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.

def point_eval(M, C, x):
    """
    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.
    """

    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)