code: def point_eval(M, C, x):
      
  
      x_ = regularize_array(x)
  
      M_out = empty(x_.shape[0])
      V_out = empty(x_.shape[0])
  
      if isinstance(C, pymc.gp.BasisCovariance):
          y_size = len(C.basis)
      elif C.obs_mesh is not None:
          y_size = C.obs_mesh.shape[0]
      else:
          y_size = 1
  
      n_chunks = ceil(y_size * x_.shape[0] / float(chunksize))
      bounds = array(linspace(0, x_.shape[0], n_chunks + 1), dtype='int')
      cmin = bounds[:-1]
      cmax = bounds[1:]
      for (cmin, cmax) in zip(bounds[:-1], bounds[1:]):
          x__ = x_[cmin:cmax]
          V_out[cmin:cmax], Uo_Cxo = C(x__, regularize=False, return_Uo_Cxo=True)
          M_out[cmin:cmax] = M(x__, regularize=False, Uo_Cxo=Uo_Cxo)
  
      if len(x.shape) > 1:
          targ_shape = x.shape[:-1]
      else:
          targ_shape = x.shape
      return M_out.reshape(targ_shape), V_out.reshape(targ_shape)
docs: Evaluates M(x) and C(x).
  
      Minimizes computation; evaluating M(x) and C(x) separately would
      evaluate the off-diagonal covariance term twice, but callling
      point_eval(M,C,x) would only evaluate it once.
  
      Also chunks the evaluations if the off-diagonal term.