'This function determines the BLUP parameters and evaluates the reduced
likelihood function for the given autocorrelation parameters theta.
Maximizing this function wrt the autocorrelation parameters theta is
equivalent to maximizing the likelihood of the assumed joint Gaussian
distribution of the observations y evaluated onto the design of
experiments X.
Parameters
theta : array_like, optional
An array containing the autocorrelation parameters at which the
Gaussian Process model parameters should be determined.
Default uses the built-in autocorrelation parameters
(ie ``theta = self.theta_``).
Returns
reduced_likelihood_function_value : double
The value of the reduced likelihood function associated to the
given autocorrelation parameters theta.
par : dict
A dictionary containing the requested Gaussian Process model
parameters:
- ``sigma2`` is the Gaussian Process variance.
- ``beta`` is the generalized least-squares regression weights for
Universal Kriging or given beta0 for Ordinary Kriging.
- ``gamma`` is the Gaussian Process weights.
- ``C`` is the Cholesky decomposition of the correlation
matrix [R].
- ``Ft`` is the solution of the linear equation system
[R] x Ft = F
- ``G`` is the QR decomposition of the matrix Ft.'

def reduced_likelihood_function(self, theta=None):

    check_is_fitted(self, 'X')
    if (theta is None):
        theta = self.theta_
    reduced_likelihood_function_value = (- np.inf)
    par = {}
    n_samples = self.X.shape[0]
    D = self.D
    ij = self.ij
    F = self.F
    if (D is None):
        (D, ij) = l1_cross_distances(self.X)
        if ((np.min(np.sum(D, axis=1)) == 0.0) and (self.corr != correlation.pure_nugget)):
            raise Exception('Multiple    X    are    not    allowed')
        F = self.regr(self.X)
    r = self.corr(theta, D)
    R = (np.eye(n_samples) * (1.0 + self.nugget))
    R[(ij[:, 0], ij[:, 1])] = r
    R[(ij[:, 1], ij[:, 0])] = r
    try:
        C = linalg.cholesky(R, lower=True)
    except linalg.LinAlgError:
        return (reduced_likelihood_function_value, par)
    Ft = linalg.solve_triangular(C, F, lower=True)
    (Q, G) = linalg.qr(Ft, mode='economic')
    sv = linalg.svd(G, compute_uv=False)
    rcondG = (sv[(-1)] / sv[0])
    if (rcondG < 1e-10):
        sv = linalg.svd(F, compute_uv=False)
        condF = (sv[0] / sv[(-1)])
        if (condF > 1000000000000000.0):
            raise Exception('F    is    too    ill    conditioned.    Poor    combination    of    regression    model    and    observations.')
        else:
            return (reduced_likelihood_function_value, par)
    Yt = linalg.solve_triangular(C, self.y, lower=True)
    if (self.beta0 is None):
        beta = linalg.solve_triangular(G, np.dot(Q.T, Yt))
    else:
        beta = np.array(self.beta0)
    rho = (Yt - np.dot(Ft, beta))
    sigma2 = ((rho ** 2.0).sum(axis=0) / n_samples)
    detR = (np.diag(C) ** (2.0 / n_samples)).prod()
    reduced_likelihood_function_value = ((- sigma2.sum()) * detR)
    par['sigma2'] = (sigma2 * (self.y_std ** 2.0))
    par['beta'] = beta
    par['gamma'] = linalg.solve_triangular(C.T, rho)
    par['C'] = C
    par['Ft'] = Ft
    par['G'] = G
    return (reduced_likelihood_function_value, par)
