Method 3: Multivariate model mixing with a Gaussian process

This method uses the same framework as the previous method, but now includes a Gaussian process (GP) in the mixing.

A diagnostic tool that helps with determining whether or not our mixed model result is reasonable is the Mahalanobis distance, calculated as

\[ D^{2}_{MD} = (\mathbf{y} - \mathbf{m})^{T}\textit{K}^{-1}(\mathbf{y} - \mathbf{m}), \]

and given in the functions below.

GP(g, loworder, highorder, kernel='RBF', nu=None, ci=68, error_model='informative', new=False)

The parameter settings of the kernel will be set by the user in this initial function. This class 'wraps' the scikit learn package.

Example

GP(g=np.linspace(1e-6,1.0,100), loworder=5, highorder=2, kernel="Matern", ci=68, error_model='informative')

Parameters:

Name	Type	Description	Default
g	numpy linspace	The linspace across the coupling constant space used for the GP.	required
highorder	(ndarray, float, int)	The truncation order of the large-g expansion.	required
kernel	str	The type of kernel the user wishes to use. Default is the RBF kernel; possible choices are RBF, Matern, and Rational Quadratic.	'RBF'
nu	float	The value of the Matern kernel used, if kernel="Matern". Otherwise, default is None.	None
ci	int	The uncertainty interval to use. Must be 68 or 95.	68
error_model	str	The error model to be used in the calculation. Options are 'uninformative' and 'informative'. Default is 'informative'.	'informative'
new	bool	Control variable for additional edits being made to the code for the dissertation alterations. Default is False.	False

Returns:

Type	Description
	None.

Source code in samba/gaussprocess.py

def __init__(self, g, loworder, highorder, kernel="RBF", nu=None, ci=68, error_model='informative', new=False):

    r'''
    A class that will pull from the Models class to perform GP emulation on 
    two models from the small-g expansion region to the large-g expansion region. 
    The parameter settings of the kernel will be set by the user in this 
    initial function. This class 'wraps' the scikit learn package. 

    Example:
        GP(g=np.linspace(1e-6,1.0,100), loworder=5, highorder=2, kernel="Matern",
            ci=68, error_model='informative')

    Parameters:
        g (numpy linspace): The linspace across the coupling constant space 
            used for the GP.

        highorder (numpy.ndarray, float, int): The truncation order of the 
            large-g expansion. 

        kernel (str): The type of kernel the user wishes to use. Default is 
            the RBF kernel; possible choices are RBF, Matern, and Rational 
            Quadratic. 

        nu (float): The value of the Matern kernel used, if kernel="Matern". 
            Otherwise, default is None.

        ci (int): The uncertainty interval to use. Must be 68 or 95. 

        error_model (str): The error model to be used in the calculation. 
            Options are 'uninformative' and 'informative'. Default is 'informative'. 

        new (bool): Control variable for additional edits being made to the code
            for the dissertation alterations. Default is False.

    Returns:
        None.
    ''' 

    #set up the prediction array as a class variable for use later
    self.gpredict = np.copy(g)

    #extract uncertainty interval for later use
    self.ci = ci 

    # extract the control variable 
    self.new = new

    #check type and assign class variables
    if isinstance(loworder, float) == True or isinstance(loworder, int) == True:
        loworder = np.array([loworder])

    if isinstance(highorder, float) == True or isinstance(highorder, int) == True:
        highorder = np.array([highorder])

    self.loworder = loworder 
    self.highorder = highorder 

    # Models(), Uncertainties()
    self.m = Models(self.loworder, self.highorder)
    self.u = Uncertainties(error_model)

    # instantiate the class variable error_model for ease class crossing
    self.error_model = self.u.error_model

    # integral length
    self.gint = np.empty([])

    # kernel set-up for the rest of the class (one-dimensional)
    kconstant = kernels.ConstantKernel(1.0)

    if kernel == "RBF":
        k = kernels.RBF(length_scale=0.5, length_scale_bounds=(1e-5,1e5))
    elif kernel == "Matern":
        if nu is None:
            raise ValueError('Matern kernel must be supplied a value for nu.')
        else:
            k = kernels.Matern(length_scale=0.4, length_scale_bounds=(1e-5,1e5), nu=nu)
    elif kernel == "Rational Quadratic":
        k = kernels.RationalQuadratic(length_scale=1.0, alpha=1)
    else:
        raise ValueError('Please choose an available kernel.')

    self.kern = kconstant * k

    return None

MD_set(pts=3, plot=False)

Takes the training set of points and uses them to cut the testing set to their limits. This reduces the MD calculation to the region of interest.

Example

GP.MD_set()

Parameters:

Name	Type	Description	Default
pts	int	The number of points to use to calculate the Mahalanobis distance. Can be any number up to the size of self.gpredict.	3
plot	bool	The option to plot the MD points across the input space. Default is False.	False

Returns:

Name	Type	Description
md_g	ndarray	The input values used in the MD calculation.
md_mean	ndarray	The mean values from the GP corresponding to the md_g points.
md_sig	ndarray	The error bars corresponding to the md_g points.
md_cov	ndarray	The covariance matrix corresponding to the md_g points.

Source code in samba/gaussprocess.py

def MD_set(self, pts=3, plot=False):

    r'''
    Takes the training set of points and uses them to cut the
    testing set to their limits. This reduces the MD calculation
    to the region of interest.  

    Example:
        GP.MD_set()

    Parameters:
        pts (int): The number of points to use to calculate the Mahalanobis
            distance. Can be any number up to the size of self.gpredict. 

        plot (bool): The option to plot the MD points across the input space. 
            Default is False. 

    Returns:
        md_g (numpy.ndarray): The input values used in the MD calculation.

        md_mean (numpy.ndarray): The mean values from the GP corresponding 
            to the md_g points.

        md_sig (numpy.ndarray): The error bars corresponding to the md_g 
            points.

        md_cov (numpy.ndarray): The covariance matrix corresponding to the md_g 
            points.
    '''

    #import the GP mean, cov, and errors for the prediction set
    GP_mean = self.meanp
    GP_err = self.sigp
    GP_cov = self.cov

    #calculate the variance at each expansion point from the next term
    lowvar = self.u.variance_low(self.gpredict, self.loworder[0])
    lowerr = np.sqrt(lowvar)
    highvar = self.u.variance_high(self.gpredict, self.highorder[0])
    hierr = np.sqrt(highvar)

    #compare the values and choose where the gap is
    for i in range(len(lowerr)):
        if GP_err[i] < lowerr[i]:
            index_lowerr = i
            break

    for i in range(len(hierr)-1, -1, -1):
        if GP_err[i] < hierr[i]: 
            index_hierr = i 
            break

    #cut the GP array into the gap
    md_g = self.gpredict[index_lowerr:index_hierr]
    self.gint = md_g.copy()
    md_mean = GP_mean[index_lowerr:index_hierr]
    md_sig = GP_err[index_lowerr:index_hierr]
    md_cov = GP_cov[index_lowerr:index_hierr, index_lowerr:index_hierr]

    #select points in g
    self.lenpts = pts
    points = self.create_points(int(self.lenpts), md_g[0], md_g[-1])
    #print('Location of MD points in g: ', points)

    #find the indices
    indices = np.zeros([self.lenpts])
    for i in range(self.lenpts):
        indices[i] = self.nearest_value(md_g, points[i])

    #convert to integer array
    indices = indices.astype(int)

    #pick the points out of the arrays
    md_g = md_g[indices]
    md_mean = md_mean[indices]
    md_sig = md_sig[indices]
    md_cov = md_cov[np.ix_(indices, indices)]

    #plot the check the location of the points
    if plot is True:
        plt.xlim(0.,1.)
        plt.plot(md_g, np.ones(len(md_g)), 'k.')

    return md_g, md_mean, md_sig, md_cov

create_points(N, a, b) staticmethod

A code to create a given number of points from a linspace evenly from points a to b.

Example

GP.create_points(N=3, a=0.0, b=1.0)

Parameters:

Name	Type	Description	Default
N	int	The number of points desired.	required
a	(float, int)	The left endpoint of the region of interest.	required
b	(float, int)	The right endpoint of the region of interest.	required

Returns:

Name	Type	Description
pts	ndarray	The resulting array of points.

Source code in samba/gaussprocess.py

@staticmethod
def create_points(N, a, b):

    r'''
    A code to create a given number of points from a 
    linspace evenly from points a to b. 

    Example:
        GP.create_points(N=3, a=0.0, b=1.0)

    Parameters:
        N (int): The number of points desired.

        a (float, int): The left endpoint of the region of 
            interest. 

        b (float, int): The right endpoint of the region of 
            interest. 

    Returns:
        pts (numpy.ndarray): The resulting array of points. 
    '''

    #create the linspace with endpoints
    pts_array = np.linspace(a, b, N+2)

    #remove the first and last point
    pts = pts_array[1:-1]

    return pts

mahalanobis(y, mean, inv=None, chol=False, svd=False) staticmethod

A diagnostic testing function that can calculate the Mahalanobis distance for a given set of mean, covariance data and a vector.

1). Calculate the MD of the predictions of the GP using

the inverse covariance matrix (usual method);

  2). Calculate the MD of the predictions to construct a 
      reference distribution using the inverse covariance
      matrix (usual method);
  3). Calculate the Cholesky decomposition of the MD
      information;
  4). Perform an SVD analysis and send back the MD 
      calculated via SVD.

Example

GP.MD(y=np.array([]), mean=np.array([]), inv=numpy.ndarray([]), chol=False, svd=False)

Parameters:

Name	Type	Description	Default
y	ndarray	An array of predicted values from the emulator.	required
mean	ndarray	An array of true values from the true model (simulator).	required
inv	ndarray	The covariance matrix to be inverted in the MD calculation.	None
chol	bool	The option to calculate the Cholesky decomposition of the data.	False
svd	bool	An option to perform the SVD analysis of the MD data. To use, must also have a covariance matrix sent to inv.	False

Returns:

Name	Type	Description
md	float	(if calculating MD) The Mahalanobis distance.
chol_decomp	ndarray	(if calculating Cholesky decomposition) The Cholesky decomposition results.
svderrs	ndarray	(if calculating SVD) The SVD errors at each point in the MD testing set.
		svd_md (float) (if calculating SVD) The Mahalanobis distance.

Source code in samba/gaussprocess.py

@staticmethod
def mahalanobis(y, mean, inv=None, chol=False, svd=False):

    r'''
    A diagnostic testing function that can calculate the Mahalanobis 
    distance for a given set of mean, covariance data and a vector. 

    Uses: 1). Calculate the MD of the predictions of the GP using
              the inverse covariance matrix (usual method);
          2). Calculate the MD of the predictions to construct a 
              reference distribution using the inverse covariance
              matrix (usual method);
          3). Calculate the Cholesky decomposition of the MD
              information;
          4). Perform an SVD analysis and send back the MD 
              calculated via SVD. 

    Example:
        GP.MD(y=np.array([]), mean=np.array([]), inv=numpy.ndarray([]),
            chol=False, svd=False)

    Parameters:
        y (numpy.ndarray): An array of predicted values from the emulator.

        mean (numpy.ndarray): An array of true values from the true model 
            (simulator).

        inv (numpy.ndarray): The covariance matrix to be inverted in the 
            MD calculation.

        chol (bool): The option to calculate the Cholesky decomposition
            of the data. 

        svd (bool): An option to perform the SVD analysis of the MD data.
            To use, must also have a covariance matrix sent to inv. 

    Returns:
        md (float): (if calculating MD) The Mahalanobis distance. 

        chol_decomp (numpy.ndarray): (if calculating Cholesky decomposition) 
            The Cholesky decomposition results. 

        svderrs (numpy.ndarray): (if calculating SVD) The SVD errors at each   
            point in the MD testing set. 

        svd_md (float) (if calculating SVD) The Mahalanobis distance. 
    '''

    y = np.atleast_2d(y)

    #cholesky option (solves for Cholesky decomposition)
    if (inv is not None) and (chol is True):

        chol = cholesky(inv)
        errs = scl.solve_triangular(chol, (y-mean).T, lower=True).T
        chol_decomp = np.linalg.norm(errs, axis=-1)

        return chol_decomp 

    #SVD option
    if (svd is True) and (inv is not None):

        #perform SVD
        _, s, vh = np.linalg.svd(inv)
        print('Eigenvalues: ',s)
        sinv = np.linalg.inv(np.diag(s))   #inverse of eigenvalue matrix
        one = vh @ (y-mean).T
        svd_md = np.squeeze(one.T @ sinv @ one)
        print('MD^2 (SVD): ', svd_md)

        #SVD errors
        svderrs = np.zeros([len(s)])
        for i in range(len(s)):
            svderrs[i] = np.square(1.0/np.sqrt(s[i]) * np.dot(vh[i,:],(y-mean).T))

        return svderrs, svd_md

    #inverse option (normal MD calculation)
    if (chol is False) and (svd is False) and (inv is not None):

        md = np.squeeze(np.sqrt(np.diag((y - mean) @ np.linalg.inv(inv) @ (y - mean).T)))

        return md

    #if nothing is selected
    if (inv is None):
        raise ValueError('Please input a covariance matrix.')

md_plotter(md_gp, md_ref, md_mean=None, md_cov=None, hist=True, box=False)

A plotting function that allows the Mahalanobis distance to be plotted using either a histogram or a box and whisker plot, or both.

Box and whisker plot code heavily drawn from J. Melendez' gsum code (https://github.com/buqeye/gsum).

Example

GP.md_plotter(md_gp=np.array([]), md_ref=np.array([]), hist=False, box=True)

Parameters:

Name	Type	Description	Default
md_gp	float	The MD^2 value for the GP curve.	required
md_ref	ndarray	The array of MD^2 values for the reference distribution.	required
md_mean	ndarray	The values of the GP mean at the md_g points. Only used for box and whisker option; default is None.	None
md_cov	ndarray	The values of the GP covariance matrix at the md_g points. Only used for box and whisker option; default is None.	None
hist	bool	Toggle for plotting a histogram. Default is True.	True
box	bool	Toggle for plotting a box plot. Default is False.	False

Returns:

Type	Description
	None.

Source code in samba/gaussprocess.py

def md_plotter(self, md_gp, md_ref, md_mean=None, md_cov=None, hist=True, box=False):

    r'''
    A plotting function that allows the Mahalanobis distance
    to be plotted using either a histogram or a box and whisker
    plot, or both. 

    Box and whisker plot code heavily drawn from J. Melendez' gsum
    code (https://github.com/buqeye/gsum).

    Example:
        GP.md_plotter(md_gp=np.array([]), md_ref=np.array([]),
        hist=False, box=True)

    Parameters:
        md_gp (float): The MD^2 value for the GP curve. 

        md_ref (numpy.ndarray): The array of MD^2 values for the reference
            distribution.

        md_mean (numpy.ndarray): The values of the GP mean at the md_g points. 
            Only used for box and whisker option; default is None. 

        md_cov (numpy.ndarray): The values of the GP covariance matrix at the 
            md_g points. Only used for box and whisker option; default is None.

        hist (bool): Toggle for plotting a histogram. Default is True. 

        box (bool): Toggle for plotting a box plot. Default is False. 

    Returns:
        None.
    '''

    title = 'Mahalanobis Distance'
    xlabel = r'$\mathrm{D}_{\mathrm{MD}}^{2}$'

    #histogram option
    if hist is True:
        fig = plt.figure(figsize=(8,6), dpi=600)
        fig.patch.set_facecolor('white')
        ax = plt.axes()
        ax.set_xlabel(xlabel, fontsize=18)
        ax.set_title(title, fontsize=22)
        ax.set_xlim(0.0, max(md_ref))
        ax.hist(md_ref, bins=50, density=True, histtype='bar', facecolor='black', \
                ec='white', label='Reference distribution')
        ax.plot(md_gp, 0.0, marker='o', color='r', markersize=10)

        #add chi-squared to histogram
        n = 200
        x = np.linspace(0.0, max(md_ref), n)
        ax.plot(x, stats.chi2.pdf(x, df=self.lenpts), 'r', linewidth=2, label=r'$\chi^2$ (df={})'.format(self.lenpts))

        #include legend
        legend = True

    #box-and-whisker option
    if box is True:

        dist = self.ref_dist(md_mean, md_cov)

        legend = False

        #set up the figure
        fig = plt.figure(figsize=(8,6), dpi=100)
        ax = plt.axes()
        ax.set_xlabel(xlabel, fontsize=18)

        #reference distribution (using chi2, NOT md_ref)
        boxartist = self.ref_boxplot(dist, ax=ax, patch_artist=True, widths=0.8)
        gray = 'gray'
        for box in boxartist['boxes']:
            box.update(dict(facecolor='lightgrey', edgecolor=gray))
        for whisk in boxartist["whiskers"]:
            whisk.update(dict(color=gray))
        for cap in boxartist["caps"]:
            cap.update(dict(color=gray))
        for med in boxartist["medians"]:
            med.update(dict(color=gray))

        #ax.boxplot(md_ref, showfliers=False)
        ax.get_xaxis().set_ticks([])
        ax.tick_params(direction='in')
        ax.set_ylim(0,20)
        ax.set_aspect(0.25)
        sns.despine(offset=0, bottom=True, ax=ax)

        #plot the individual GP MD value
        ax.plot(1.0, md_gp, color='red', marker='o', markersize=10)

    #finish up plot
    if legend is True:
        ax.legend(loc='upper right', fontsize=18)

    plt.show()

    return None

md_squared(md_g, md_mean, md_cov, n_curves=1000)

A wrapper for the Mahalanobis distance calculation for the reference distribution and the GP curve. To calculate the Cholesky decomposition or to perform an SVD analysis, consult GP.mahalanobis() below.

Example

GP.md_squared(md_g=np.linspace, md_mean=np.array([]), md_cov=np.array([,]), n_curves=1000)

Parameters:

Name	Type	Description	Default
md_g	linspace	The points in input space g from the GP.MD_set() function.	required
md_mean	ndarray	The values of the GP mean at the md_g points.	required
md_cov	ndarray	The values of the GP covariance matrix at the md_g points.	required
n_curves	int	The number of curves from the reference distribution that are drawn for the MD^2 calculation (md_ref).	1000

Returns:

Name	Type	Description
md_gp	float	The individual MD^2 value for the GP curve.

md_ref (numpy.ndarray): The array of MD^2 values from the reference distribution.

Source code in samba/gaussprocess.py

def md_squared(self, md_g, md_mean, md_cov, n_curves=1000):

    r'''
    A wrapper for the Mahalanobis distance calculation for the
    reference distribution and the GP curve. To calculate the 
    Cholesky decomposition or to perform an SVD analysis, consult
    GP.mahalanobis() below. 

    Example:
        GP.md_squared(md_g=np.linspace, md_mean=np.array([]), 
                      md_cov=np.array([,]), n_curves=1000)

    Parameters:
        md_g (numpy.linspace): The points in input space g from the GP.MD_set() 
            function. 

        md_mean (numpy.ndarray): The values of the GP mean at the md_g points. 

        md_cov (numpy.ndarray): The values of the GP covariance matrix at the 
            md_g points. 

        n_curves (int): The number of curves from the reference distribution that
            are drawn for the MD^2 calculation (md_ref). 

    Returns:
        md_gp (float): The individual MD^2 value for the GP curve. 

    md_ref (numpy.ndarray): The array of MD^2 values from the reference 
        distribution.
    '''

    #calculate the ref distribution MDs
    dist = self.ref_dist(md_mean, md_cov)
    y = self.sample_ref(dist, n_curves)
    md = np.ones([n_curves])
    for i in range(n_curves):
        md[i] = self.mahalanobis(y[:,i].T, md_mean, inv=md_cov, chol=False, svd=False)

    #MD^2 (ref)
    md_ref = md**2.0 

    #calculate the GP MD 
    fval = self.m.true_model(md_g)
    mdgp = self.mahalanobis(fval.T, md_mean, inv=md_cov, chol=False, svd=False)

    #MD^2 (GP)
    md_gp = mdgp**2.0

    return md_gp, md_ref

nearest_value(array, value) staticmethod

A static method to find the index of the nearest value of an array to a desired value.

Example

GP.nearest_value(array=numpy.ndarray, value=5)

Parameters:

Name	Type	Description	Default
array	ndarray	The array of values to search.	required
value	int	The desired value to search the array for.	required

index (int): The index of the nearest value of the array to the desired value.

Source code in samba/gaussprocess.py

@staticmethod
def nearest_value(array, value):

    r'''
    A static method to find the index of the nearest value
    of an array to a desired value. 

    Example:
        GP.nearest_value(array=numpy.ndarray, value=5)

    Parameters:
        array (numpy.ndarray): The array of values to search. 

        value (int): The desired value to search the array for. 

    Returns:
    index (int): The index of the nearest value of the array
        to the desired value. 
    '''

    #calculate the difference between each point
    abs_val = np.abs(array - value)

    #find the smallest difference in the array
    index = abs_val.argmin()

    return index

plot_training(gs, datas, sigmas)

A simple plotter to plot the trained GP results and models, as well as the points at which the GP was trained.

Example

GP.plot_training(gs=np.array([]), datas=np.array([]), sigmas=np.array([]))

Parameters:

Name	Type	Description	Default
gs	ndarray	Points chosen by GP.training_set() in input space g.	required
datas	ndarray	Corresponding values of the series expansions at gs.	required
sigmas	ndarray	Corresponding error model results at each training point.	required

Returns:

Type	Description
	None.

Source code in samba/gaussprocess.py

def plot_training(self, gs, datas, sigmas):

    r'''
    A simple plotter to plot the trained GP results and models, 
    as well as the points at which the GP was trained. 

    Example:
        GP.plot_training(gs=np.array([]), datas=np.array([]),
            sigmas=np.array([]))

    Parameters:
        gs (numpy.ndarray): Points chosen by GP.training_set() in input 
            space g.

        datas (numpy.ndarray): Corresponding values of the series expansions 
            at gs.

        sigmas (numpy.ndarray): Corresponding error model results at each 
            training point.

    Returns:
        None.
    '''

    #set up the plot
    fig = plt.figure(figsize=(8,6), dpi=600)
    ax = plt.axes()
    fig.patch.set_facecolor('white')
    ax.tick_params(axis='x', labelsize=18)
    ax.tick_params(axis='y', labelsize=18)
    ax.locator_params(nbins=8)
    ax.xaxis.set_minor_locator(AutoMinorLocator())
    ax.yaxis.set_minor_locator(AutoMinorLocator())
    ax.set_xlim(0.0, max(self.gpredict))
    ax.set_ylim(1.0,3.0)
    ax.set_xlabel('g', fontsize=22)
    ax.set_ylabel('F(g)', fontsize=22)
    ax.set_title('F(g): training set', fontsize=22)
    ax.plot(self.gpredict, self.m.true_model(self.gpredict), 'k', label='True model')

    #plot the data
    ax.errorbar(self.gtrlow, self.datatrlow, yerr=self.lowsigma, color='red', fmt='o', markersize=4, \
                capsize=4, label=r'$f_s$ ($N_s$ = {}) data'.format(self.loworder[0]))
    ax.errorbar(self.gtrhigh, self.datatrhigh, yerr=self.highsigma, color='blue', fmt='o', markersize=4, \
                capsize=4, label=r'$f_l$ ($N_l$ = {}) data'.format(self.highorder[0]))

    #plot the chosen training points over the whole training set
    ax.errorbar(gs, datas, yerr=sigmas, color='black', fmt='o', markersize=4, capsize=4, label='Training data')

    ax.legend(fontsize=18, loc='upper right')
    plt.show()

    #save figure option
    # response = input('Would you like to save this figure? (yes/no)')

    # if response == 'yes':
    #     name = input('Enter a file name (include .jpg, .png, etc.)')
    #     fig.savefig(name, bbox_inches='tight')

    return None

plot_validate(intervals)

A simple plotter to show the results of the GP predictions at new points in g.

Example

GP.plot_validate(intervals=np.array([,]))

Parameters:

Name	Type	Description	Default
intervals	ndarray	The uncertainty band around the prediction set.	required

Returns:

Type	Description
	None.

Source code in samba/gaussprocess.py

def plot_validate(self, intervals):

    r'''
    A simple plotter to show the results of the GP 
    predictions at new points in g. 

    Example:
        GP.plot_validate(intervals=np.array([,]))

    Parameters:
        intervals (numpy.ndarray): The uncertainty band around the 
            prediction set.

    Returns:
        None.
    '''

    #plot the results
    fig = plt.figure(figsize=(8,6), dpi=600)
    ax = plt.axes()
    fig.patch.set_facecolor('white')
    ax.tick_params(axis='x', labelsize=18)
    ax.tick_params(axis='y', labelsize=18)
    ax.locator_params(nbins=8)
    ax.xaxis.set_minor_locator(AutoMinorLocator())
    ax.yaxis.set_minor_locator(AutoMinorLocator())
    ax.set_xlim(0.0, max(self.gpredict))
    ax.set_ylim(1.0,3.0)
    ax.set_xlabel('g', fontsize=22)
    ax.set_ylabel('F(g)', fontsize=22)
    ax.set_title('F(g): GP predictions', fontsize=22)
    ax.plot(self.gpredict, self.m.true_model(self.gpredict), 'k', label='True model')

    #plot the data
    if self.bars is True:
        ax.errorbar(self.gtrlow, self.datatrlow, self.lowsigma, color="red", fmt='o', markersize=4, \
            capsize=4, alpha = 0.4, label=r"$f_s$ ($N_s$ = {})".format(self.loworder[0]), zorder=1)
        ax.errorbar(self.gtrhigh, self.datatrhigh, self.highsigma, color="blue", fmt='o', markersize=4, \
            capsize=4, alpha=0.4, label=r"$f_l$ ($N_l$ = {})".format(self.highorder[0]), zorder=1)
    else:
        ax.errorbar(self.gs, self.datas, yerr=self.sigmas, color="red", fmt='o', markersize=4, \
            capsize=4, label=r"Training data", zorder=10)

        # calculate intervals (could all be better but this isn't my focus anymore)
        _, _, interval_low, interval_high = self.fdagger(self.gpredict)

        #plot the small-g expansions and error bands
        for i,j in zip(range(len(self.loworder)), self.loworder):
            ax.plot(self.gpredict, self.m.low_g(self.gpredict)[i,:], 'r--', label=r'$f_s$ ($N_s$ = {})'.format(j))

        for i in range(len(self.loworder)):
            ax.plot(self.gpredict, interval_low[i, :, 0], 'r', linestyle='dotted', \
                label=r'$f_s$ ($N_s$ = {}) {}\% CI'.format(self.loworder[i], int(self.ci)))
            ax.plot(self.gpredict, interval_low[i, :, 1], 'r', linestyle='dotted')

        #for each large-g order, calculate and plot
        for i,j in zip(range(len(self.highorder)), self.highorder):
            ax.plot(self.gpredict, self.high_g(self.gpredict)[i,:], 'b--', label=r'$f_l$ ($N_l$ = {})'.format(j))

        for i in range(len(self.highorder)):
            ax.plot(self.gpredict, interval_high[i, :, 0], 'b', linestyle='dotted', \
                label=r'$f_l$ ($N_l$ = {}) {}\% CI'.format(self.highorder[i], int(self.ci)))
            ax.plot(self.gpredict, interval_high[i, :, 1], 'b', linestyle='dotted')
        ax.set_xlim(0.0, max(self.gpredict)+0.01)
    ax.plot(self.gpred, self.meanp, 'g', label='Predictions', zorder=2)
    ax.plot(self.gpred, intervals[:,0], color='green', linestyle='dotted', label=r'{}$\%$ CI'.format(self.ci), zorder=2)
    ax.plot(self.gpred, intervals[:,1], color='green', linestyle='dotted', zorder=2)
    ax.fill_between(self.gpred[:,0], intervals[:,0], intervals[:,1], color='green', alpha=0.3, zorder=10)

    ax.legend(fontsize=14, loc='upper right')
    plt.show()

    #save figure option
    # response = input('Would you like to save this figure? (yes/no)')

    # if response == 'yes':
    #     name = input('Enter a file name (include .jpg, .png, etc.)')
    #     fig.savefig(name, bbox_inches='tight')

    return None   # again this is the true answer, use this for the thesis

ref_boxplot(dist, q1=0.25, q3=0.75, whislo=0.025, whishi=0.975, ax=None, **kwargs) staticmethod

Taken from the gsum code written by J. Melendez (https://github.com/buqeye/gsum).

Source code in samba/gaussprocess.py

@staticmethod 
def ref_boxplot(dist, q1=0.25, q3=0.75, whislo=0.025, whishi=0.975, ax=None, **kwargs):

    r'''
    Taken from the gsum code written by J. Melendez (https://github.com/buqeye/gsum).
    '''

    stat_dict = [{'med': dist.median(), 'q1': dist.ppf(q1), 'q3': dist.ppf(q3),
                  'whislo': dist.ppf(whislo), 'whishi': dist.ppf(whishi)}]

    return ax.bxp(stat_dict, showfliers=False, **kwargs)

ref_dist(mean, cov) staticmethod

Constructs a multivariate normal distribution to act as a reference distribution for the Mahalanobis distance calculation.

Example

Diagnostics.ref_dist(mean=np.array([]), cov=np.array([]))

Parameters:

Name	Type	Description	Default
mean	ndarray	The mean of the GP (given by the prediction set).	required
cov	ndarray	The covariance matrix of the GP (given by the prediction set).	required

Returns:

Name	Type	Description
dist	object	A multivariate normal distribution that can be used to generate samples for the reference distribution.

Source code in samba/gaussprocess.py

@staticmethod
def ref_dist(mean, cov):

    r'''
    Constructs a multivariate normal distribution to act
    as a reference distribution for the Mahalanobis distance
    calculation. 

    Example:
        Diagnostics.ref_dist(mean=np.array([]), cov=np.array([]))

    Parameters:
        mean (numpy.ndarray): The mean of the GP (given by the 
            prediction set). 

        cov (numpy.ndarray): The covariance matrix of the GP 
            (given by the prediction set). 

    Returns:
        dist (object): A multivariate normal distribution that can 
            be used to generate samples for the reference distribution. 
    '''

    dist = stats.multivariate_normal(mean=mean, cov=cov)

    return dist

sample_ref(dist, n_curves) staticmethod

Generate some sample curves from the reference distribution.

Example

Diagnostics.sample_ref(dist, n_curves=10)

Parameters:

Name	Type	Description	Default
dist	object	The reference distribution object.	required
n_curves	int	The number of draws from the reference distribution.	required

Returns:

Name	Type	Description
samples	ndarray	The array of curves from the distribution.

Source code in samba/gaussprocess.py

@staticmethod
def sample_ref(dist, n_curves):

    r'''
    Generate some sample curves from the reference distribution.

    Example:
        Diagnostics.sample_ref(dist, n_curves=10)

    Parameters:
        dist (object): The reference distribution object. 

        n_curves (int): The number of draws from the reference 
            distribution.

    Returns:
        samples (numpy.ndarray): The array of curves from the 
            distribution. 
    '''

    samples = dist.rvs(n_curves).T

    return samples

training(error=True, method=2, plot=True)

A function that links the model data and the training function in scikit learn, and plots the training data using GP.plot_training().

Example

GP.training(error=False, method=3)

Parameters:

Name	Type	Description	Default
error	bool	A boolean variable to toggle use of a truncation error model in the kernel during training. Default is True.	True
method	int	The method used for determining the training points. Options: 1,2,3. For an extensive explanation of the methods, see the paper.	2
plot	bool	Option to plot the training set with series expansions and true model. Default is True.	True

Returns:

Name	Type	Description
sk	object	The object storing all training information from the sklearn regression performed on the data.

Source code in samba/gaussprocess.py

    def training(self, error=True, method=2, plot=True):

        r'''
        A function that links the model data and the training function in 
        scikit learn, and plots the training data using GP.plot_training().

        Example:
            GP.training(error=False, method=3)

        Parameters:
            error (bool): A boolean variable to toggle use of a truncation error model 
                in the kernel during training. Default is True.

            method (int): The method used for determining the training points. Options: 
                1,2,3. For an extensive explanation of the methods, see the paper.

            plot (bool): Option to plot the training set with series expansions and 
                true model. Default is True. 

        Returns:
            sk (object): The object storing all training information from the sklearn 
                regression performed on the data.
        '''

        #first set the method
        self.method = method 

        #call the training set generator function
        gs, datas, sigmas = self.training_set()

#         ### specific test ###
#         # split up the set
#         gs = np.array([gs[1], gs[2], gs[3]])
#         datas = np.array([datas[1], datas[2], datas[3]])
#         sigmas = np.array([sigmas[1], sigmas[2], sigmas[3]])

        #make a gs class variable for weights use
        self.gs = gs 
        self.datas = datas
        self.sigmas = sigmas

        #make column vectors for the regressor
        gc = gs.reshape(-1,1)
        datac = datas.reshape(-1,1)

        #take the data point uncertainty into the kernel 
        if error == True:
            self.alpha = np.square(sigmas)
        else:
            self.alpha = 1e-12

        #use GPR and kernel to train
        m = GaussianProcessRegressor(kernel=self.kern, alpha=self.alpha, n_restarts_optimizer=20, normalize_y=True)

        #fit the GP to the training data
        self.sk = m.fit(gc, datac)

        #print the optimized parameters for the user
        print('Gaussian process parameters: {}'.format(m.kernel_))

        #plot the results
        if plot is True:
            self.plot_training(gs, datas, sigmas)

        return self.sk

training_set()

An internal function to calculate the necessary training data set from the input prediction set.

Example

GP.training_set()

Returns: gs (numpy.ndarray): The modified array of input values for the training.

datas (numpy.ndarray): The modified array of data values for the training.

sigmas (numpy.ndarray): The modified array of the truncation errors for the training.

Source code in samba/gaussprocess.py

    def training_set(self):

        r'''
        An internal function to calculate the necessary training data set from
        the input prediction set. 

        Example:
            GP.training_set() 

        Parameters:
            None. 

        Returns:
        gs (numpy.ndarray): The modified array of input values for the training. 

        datas (numpy.ndarray): The modified array of data values for the training. 

        sigmas (numpy.ndarray): The modified array of the truncation errors for 
            the training. 
        '''

        #set up the training set from the prediction set (offset by midpoint)  # this should be fine
        self.midpoint = (self.gpredict[1] - self.gpredict[0]) / 2.0
        gtrainingset = np.linspace(min(self.gpredict)+self.midpoint, max(self.gpredict)+self.midpoint, len(self.gpredict))

        #stop the training set, negative curvature
        if self.loworder[0] % 4 == 2 or self.loworder[0] % 4 == 3:
            for i in range(len(gtrainingset)):
                if self.m.low_g(gtrainingset[i]) < -1.0:
                    lowindex = i-1
                    break

        #stop the training set, positive curvature
        elif self.loworder[0] % 4 == 0 or self.loworder[0] % 4 == 1:
            for i in range(len(gtrainingset)):
                if self.m.low_g(gtrainingset[i]) > 3.0:
                    lowindex = i-1
                    break

        #stop the training set, even orders (positive curvature)
        if self.highorder[0] % 2 == 0:
            for i in range(len(gtrainingset)):
                if self.m.high_g(gtrainingset[i]) > 3.0:
                    highindex = i+1
                else:
                    break

        #stop the training set, odd orders (negative curvature)
        else:
            for i in range(len(gtrainingset)):
                if self.m.high_g(gtrainingset[i]) < -1.0:
                    highindex = i+1
                else:
                    break

        # slice the training set for the two models (this should be fine)
        self.gtrlow = gtrainingset[:lowindex]
        self.gtrhigh = gtrainingset[highindex:]

        # calculate the data at each point
        self.datatrlow = self.m.low_g(self.gtrlow)[0,:]
        self.datatrhigh = self.m.high_g(self.gtrhigh)[0,:]

        # calculate the variance at each point from the next term
        lowvariance = self.u.variance_low(self.gtrlow, self.loworder[0])
        self.lowsigma = np.sqrt(lowvariance)
        highvariance = self.u.variance_high(self.gtrhigh, self.highorder[0])
        self.highsigma = np.sqrt(highvariance)

        # find the values of g in the other set to determine location of points
        index_ghigh = (np.where(self.gtrhigh == self.gtrlow[-1])[0])[0]

        # value of g at the optimal red points
        if self.new is False:
            pt1 = 0.0656575
            pt2 = 0.1161625
        else:
            pt1 = 0.02
            pt2 = 0.08
            pt3 = 0.12

        # method 1: using g=0.6 as a training point
        pttest = 0.6  
        indexptest = self.nearest_value(self.gtrhigh, pttest) 

        # method 2: finding based on error (2%)
        for i in range(len(self.gtrhigh)-1, -1, -1):
            if self.highsigma[i] >= 0.02*self.datatrhigh[i]:
                indexerror = i
                break 

        # find the values in the training array closest to the points (this is good regardless of what happens)
        indexpt1 = self.nearest_value(self.gtrlow, pt1)
        indexpt2 = self.nearest_value(self.gtrlow, pt2)
        if self.new is True:
            indexpt3 = self.nearest_value(self.gtrlow, pt3)

        # create two points on either side (highpoint = 20)
        if self.new is False:
            glowtr = np.array([self.gtrlow[indexpt1], self.gtrlow[indexpt2]])
            datalowtr = np.array([self.datatrlow[indexpt1], self.datatrlow[indexpt2]])
            sigmalowtr = np.array([self.lowsigma[indexpt1], self.lowsigma[indexpt2]])
        else:
            glowtr = np.array([self.gtrlow[indexpt1], self.gtrlow[indexpt2], self.gtrlow[indexpt3]])
            datalowtr = np.array([self.datatrlow[indexpt1], self.datatrlow[indexpt2], self.datatrlow[indexpt3]])
            sigmalowtr = np.array([self.lowsigma[indexpt1], self.lowsigma[indexpt2], self.lowsigma[indexpt3]])

        # choose training points depending on method entered (method 1 in thesis)
        if self.method == 1:
            ghightr = np.array([self.gtrhigh[indexptest], self.gtrhigh[-1]])
            datahightr = np.array([self.datatrhigh[indexptest], self.datatrhigh[-1]])
            sigmahightr = np.array([self.highsigma[indexptest], self.highsigma[-1]])

#         elif self.method == 2:
#             if self.new is True:
#                 newpt = self.gtrhigh[-1] - self.gtrhigh[index_ghigh]/4.0
#                 indexnewpt = self.nearest_value(self.gtrhigh, newpt)
#                 ghightr = np.array([self.gtrhigh[index_ghigh], self.gtrhigh[indexnewpt], self.gtrhigh[-1]])
#                 datahightr = np.array([self.datatrhigh[index_ghigh], self.datatrhigh[indexnewpt], self.datatrhigh[-1]])
#                 sigmahightr = np.array([self.highsigma[index_ghigh], self.highsigma[indexnewpt], self.highsigma[-1]])
# #                 ghightr = np.array([self.gtrhigh[index_ghigh], self.gtrhigh[-1]])
# #                 datahightr = np.array([self.datatrhigh[index_ghigh], self.datatrhigh[-1]])
# #                 sigmahightr = np.array([self.highsigma[index_ghigh], self.highsigma[-1]])

#             else:
#                 ghightr = np.array([self.gtrhigh[index_ghigh], self.gtrhigh[-1]])
#                 datahightr = np.array([self.datatrhigh[index_ghigh], self.datatrhigh[-1]])
#                 sigmahightr = np.array([self.highsigma[index_ghigh], self.highsigma[-1]])

        # method 2 in thesis
        elif self.method == 3:
            if self.new is True:
                newpt = (self.gtrhigh[-1] - self.gtrhigh[indexerror])/2.0 + self.gtrhigh[indexerror]
                indexnewpt = self.nearest_value(self.gtrhigh, newpt)
                ghightr = np.array([self.gtrhigh[indexerror], self.gtrhigh[indexnewpt], self.gtrhigh[-1]])
                datahightr = np.array([self.datatrhigh[indexerror], self.datatrhigh[indexnewpt], self.datatrhigh[-1]])
                sigmahightr = np.array([self.highsigma[indexerror], self.highsigma[indexnewpt], self.highsigma[-1]])   
            else:
                ghightr = np.array([self.gtrhigh[indexerror], self.gtrhigh[-1]])
                datahightr = np.array([self.datatrhigh[indexerror], self.datatrhigh[-1]])
                sigmahightr = np.array([self.highsigma[indexerror], self.highsigma[-1]])   

        # somewhere in the middle of the two points (method 3 in thesis)  --> this one we use?
        elif self.method == 4:
            indexfhigh = indexerror
            fmiddle = (self.gtrhigh[-1] - self.gtrhigh[indexerror])/3.0 + self.gtrhigh[indexerror]
            fmiddle2 = (self.gtrhigh[-1] - self.gtrhigh[indexerror])/3.0 + self.gtrhigh[indexerror] \
               + (self.gtrhigh[-1] - self.gtrhigh[indexerror])/3.0
            indexfmiddle = self.nearest_value(self.gtrhigh, fmiddle) 
            indexfmiddle2 = self.nearest_value(self.gtrhigh, fmiddle2)
            ghightr = np.array([self.gtrhigh[indexfhigh], \
                                self.gtrhigh[indexfmiddle], self.gtrhigh[indexfmiddle2], self.gtrhigh[-1]])
            datahightr = np.array([self.datatrhigh[indexfhigh], \
                                   self.datatrhigh[indexfmiddle], self.datatrhigh[indexfmiddle2], self.datatrhigh[-1]])
            sigmahightr = np.array([self.highsigma[indexfhigh], \
                                    self.highsigma[indexfmiddle], self.highsigma[indexfmiddle2], self.highsigma[-1]]) 

#         elif self.method == 5:
#             firsthigh = 0.6
#             indexfhigh = self.nearest_value(self.gtrhigh, firsthigh)
#             fmiddle = self.gtrhigh[-1] - self.gtrhigh[indexfhigh]/3.0  
#             indexfmiddle = self.nearest_value(self.gtrhigh, fmiddle)
#             ghightr = np.array([self.gtrhigh[indexfhigh], \
#                                 self.gtrhigh[indexfmiddle], self.gtrhigh[-1]])
#             datahightr = np.array([self.datatrhigh[indexfhigh], \
#                                    self.datatrhigh[indexfmiddle], self.datatrhigh[-1]])
#             sigmahightr = np.array([self.highsigma[indexfhigh], \
#                                     self.highsigma[indexfmiddle], self.highsigma[-1]]) 

        #concatenate these arrays and send back
        gtr = np.concatenate((glowtr, ghightr))
        datatr = np.concatenate((datalowtr, datahightr))
        sigmatr = np.concatenate((sigmalowtr, sigmahightr))

        return gtr, datatr, sigmatr 

validate(plot=True, run_taweret=False, bars=True)

A wrapper function for scikit learn's GP prediction function. This will predict the GP results with an interval and plot against the expansions using GP.plot_validate().

Example

GP.validate()

Parameters:

Name	Type	Description	Default
plot	bool	The option to plot the GP mean and variance over the testing set and true model. Default is True.	True
bars	bool	Whether to plot and show the error bands of each training point or not. Default is True.	True

Returns:

Name	Type	Description
meanp	ndarray	The mean array of the GP prediction results.
sigp	ndarray	The standard deviation array of the GP prediction results.
cov	ndarray	The covariance matrix of the GP prediction results.

Source code in samba/gaussprocess.py

def validate(self, plot=True, run_taweret=False, bars=True):

    r'''
    A wrapper function for scikit learn's GP prediction function. This will 
    predict the GP results with an interval and plot against the expansions
    using GP.plot_validate().

    Example:
        GP.validate()

    Parameters:
        plot (bool): The option to plot the GP mean and variance over the testing
            set and true model. Default is True. 

        bars (bool): Whether to plot and show the error bands of
            each training point or not. Default is True.

    Returns:
        meanp (numpy.ndarray): The mean array of the GP prediction results.

        sigp (numpy.ndarray): The standard deviation array of the GP prediction 
            results. 

        cov (numpy.ndarray): The covariance matrix of the GP prediction results. 
    '''

    # save this control variable for later
    self.bars = bars

    #make the prediction values into a column vector
    self.gpred = self.gpredict.reshape(-1,1)

    #predict the results for the validation data
    self.meanp, self.sigp = self.sk.predict(self.gpred, return_std=True)
    _, self.cov = self.sk.predict(self.gpred, return_cov=True)

    # issues right here when running wrapped in Taweret
    if run_taweret is False:
        self.meanp = self.meanp #[:,0]

    #calculate the interval for the predictions
    if self.ci == 68:
        factor = 1.0
    elif self.ci == 95:
        factor = 1.96
    intervals = np.zeros([len(self.meanp), 2])
    intervals[:,0] = self.meanp - factor*self.sigp
    intervals[:,1] = self.meanp + factor*self.sigp

    #plot the results
    if plot is True:
        self.plot_validate(intervals)

    return self.meanp, self.sigp, self.cov   # this is it, this is the answer, not the PPD one