'The Gaussian Process model fitting method.
Parameters
X : double array_like
An array with shape (n_samples, n_features) with the input at which
observations were made.
y : double array_like
An array with shape (n_samples, ) or shape (n_samples, n_targets)
with the observations of the output to be predicted.
Returns
gp : self
A fitted Gaussian Process model object awaiting data to perform
predictions.'

def fit(self, X, y):

    self._check_params()
    self.random_state = check_random_state(self.random_state)
    (X, y) = check_X_y(X, y, multi_output=True, y_numeric=True)
    self.y_ndim_ = y.ndim
    if (y.ndim == 1):
        y = y[:, np.newaxis]
    (n_samples, n_features) = X.shape
    (_, n_targets) = y.shape
    self._check_params(n_samples)
    if self.normalize:
        X_mean = np.mean(X, axis=0)
        X_std = np.std(X, axis=0)
        y_mean = np.mean(y, axis=0)
        y_std = np.std(y, axis=0)
        X_std[(X_std == 0.0)] = 1.0
        y_std[(y_std == 0.0)] = 1.0
        X = ((X - X_mean) / X_std)
        y = ((y - y_mean) / y_std)
    else:
        X_mean = np.zeros(1)
        X_std = np.ones(1)
        y_mean = np.zeros(1)
        y_std = np.ones(1)
    (D, ij) = l1_cross_distances(X)
    if ((np.min(np.sum(D, axis=1)) == 0.0) and (self.corr != correlation.pure_nugget)):
        raise Exception('Multiple    input    features    cannot    have    the    same    target    value.')
    F = self.regr(X)
    n_samples_F = F.shape[0]
    if (F.ndim > 1):
        p = F.shape[1]
    else:
        p = 1
    if (n_samples_F != n_samples):
        raise Exception('Number    of    rows    in    F    and    X    do    not    match.    Most    likely    something    is    going    wrong    with    the    regression    model.')
    if (p > n_samples_F):
        raise Exception(('Ordinary    least    squares    problem    is    undetermined    n_samples=%d    must    be    greater    than    the    regression    model    size    p=%d.' % (n_samples, p)))
    if (self.beta0 is not None):
        if (self.beta0.shape[0] != p):
            raise Exception('Shapes    of    beta0    and    F    do    not    match.')
    self.X = X
    self.y = y
    self.D = D
    self.ij = ij
    self.F = F
    (self.X_mean, self.X_std) = (X_mean, X_std)
    (self.y_mean, self.y_std) = (y_mean, y_std)
    if ((self.thetaL is not None) and (self.thetaU is not None)):
        if self.verbose:
            print('Performing    Maximum    Likelihood    Estimation    of    the    autocorrelation    parameters...')
        (self.theta_, self.reduced_likelihood_function_value_, par) = self._arg_max_reduced_likelihood_function()
        if np.isinf(self.reduced_likelihood_function_value_):
            raise Exception('Bad    parameter    region.    Try    increasing    upper    bound')
    else:
        if self.verbose:
            print('Given    autocorrelation    parameters.    Computing    Gaussian    Process    model    parameters...')
        self.theta_ = self.theta0
        (self.reduced_likelihood_function_value_, par) = self.reduced_likelihood_function()
        if np.isinf(self.reduced_likelihood_function_value_):
            raise Exception('Bad    point.    Try    increasing    theta0.')
    self.beta = par['beta']
    self.gamma = par['gamma']
    self.sigma2 = par['sigma2']
    self.C = par['C']
    self.Ft = par['Ft']
    self.G = par['G']
    if (self.storage_mode == 'light'):
        if self.verbose:
            print('Light    storage    mode    specified.    Flushing    autocorrelation    matrix...')
        self.D = None
        self.ij = None
        self.F = None
        self.C = None
        self.Ft = None
        self.G = None
    return self
