Documentation for the neutron-rich-bmm code

GP BMM structure and helper functions

SigmoidChangepoint(ls1, ls2, cbar1, cbar2, changepoint=1.0, changepoint_bounds=(1e-05, 100000.0), width=1.0, width_bounds=(1e-05, 100000.0))

Designs a non-stationary changepoint kernel that inherits from the sklearn RBF Kernel class.

The kernel is given by: .. math:: k(x_i, x_j) = (1 - \sigma(x_i)) * K1(x_i,x_j) * (1 - \sigma(x_j)) + \sigma(x_i) * K2(x_i,x_j) * \sigma(x_j)

where K1 is the first kernel and K2 is the second kernel, with the changepoint defined by a chosen switching function. The current options are 'linear' and 'sigmoid'.

Source code in src/neutron_rich_bmm/custom_kernels.py

def __init__(self, ls1, ls2, cbar1, cbar2, changepoint=1.0, changepoint_bounds=(1e-5, 1e5), 
             width=1.0, width_bounds=(1e-5, 1e5)):

    # will not touch these for now
    self.ls1 = ls1
    self.ls2 = ls2
    self.cbar1 = cbar1
    self.cbar2 = cbar2

    # optimizable parameters
    self.width = width
    self.width_bounds = width_bounds
    self.changepoint = changepoint 
    self.changepoint_bounds = changepoint_bounds

    # which kernel type am I using
    self.type = 'sigmoid'

    return None

TanhChangepoint(ls1, ls2, cbar1, cbar2, changepoint=1.0, changepoint_bounds=(1e-05, 100000.0), width=1.0, width_bounds=(1e-05, 100000.0))

Designs a non-stationary changepoint kernel that inherits from the sklearn RBF Kernel class.

The kernel is given by: .. math:: k(x_i, x_j) = (1 - \sigma(x_i)) * K1(x_i,x_j) * (1 - \sigma(x_j)) + \sigma(x_i) * K2(x_i,x_j) * \sigma(x_j)

where K1 is the first kernel and K2 is the second kernel, with the changepoint defined by a chosen switching function. The current option is 'tanh'.

Source code in src/neutron_rich_bmm/custom_kernels.py

def __init__(self, ls1, ls2, cbar1, cbar2, changepoint=1.0, changepoint_bounds=(1e-5, 1e5), 
             width=1.0, width_bounds=(1e-5, 1e5)):

    # will not touch these for now
    self.ls1 = ls1
    self.ls2 = ls2
    self.cbar1 = cbar1
    self.cbar2 = cbar2

    # optimizable parameters
    self.width = width
    self.width_bounds = width_bounds
    self.changepoint = changepoint 
    self.changepoint_bounds = changepoint_bounds

    return None

boundary_conditions(dens, pres_dict, index=0)

Helper function to find boundary conditions from the pQCD results.

Parameters:

Name	Type	Description	Default
dens	ndarray	The density array as input to find the BCs.	required
pres_dict	dict	The dictionary of pressure values corresponding to the input density array.	required

Returns:

Name	Type	Description
mu_FG	ndarray	The 1-d array of chemical potentials corresponding to the values of density that were input.
mU_FG	ndarray	The array of shape [:,None] that is used in the gsum truncation error analysis.
edens_dict	dict	The energy density values at the chosen density index, used as the BCs for the speed of sound calculation.

Source code in src/neutron_rich_bmm/eos_utils.py

def boundary_conditions(dens, pres_dict, index=0):

    '''
    Helper function to find boundary conditions
    from the pQCD results. 

    Parameters:
        dens (numpy.ndarray): The density array as 
            input to find the BCs.

        pres_dict (dict): The dictionary of pressure values 
            corresponding to the input density array.

    Returns:
        mu_FG (numpy.ndarray): The 1-d array of chemical potentials
            corresponding to the values of density that
            were input. 

        mU_FG (numpy.ndarray): The array of shape [:,None] that is
            used in the gsum truncation error
            analysis. 

        edens_dict (dict): The energy density values at the chosen
            density index, used as the BCs for the 
            speed of sound calculation.
    '''

    # call pQCD class
    pqcd = PQCD(X=1, Nf=2) # classic implementation here

    # constants
    hbarc = 197.327 # Mev fm

    # unpack dictionary
    pres_FG = pres_dict['FG']
    pres_NLO = pres_dict['NLO']
    pres_N2LO = pres_dict['N2LO']

    # set up new dictionary for eps(n) BCs
    edens_FG = dict()
    edens_NLO = dict()
    edens_N2LO = dict()

    # make mu_FG array from the selected density array (no playing around)
    n_q = dens*3.0  # n_q [fm^-3]

    # convert to GeV^3 for mu_q
    conversion_fm3 = ((1000.0)**(3.0))/((197.33)**(3.0)) # [fm^-3]  (do the opposite of this)
    n_q = n_q/conversion_fm3  # [GeV^3]

    # invert to get mu
    _, _, mu_FG = pqcd.inversion(n_mu=n_q)  # [GeV] # these are quark chemical potentials
    mU_FG = mu_FG[:, None]

    # FG BCs
    edens_FG['mean'] = ((3.0 / (2 * np.pi**2.0)) * (3.0 * np.pi**2.0 * dens/2.0)**(4.0/3.0) * hbarc)[index]
    edens_FG['lower'] = (dens*3*1000.*mu_FG - (pres_dict_FG['mean']-pres_dict_FG['std_dev']))[index]
    edens_FG['upper'] = (dens*3*1000.*mu_FG - (pres_dict_FG['mean']+pres_dict_FG['std_dev']))[index]

    # NLO BCs
    edens_NLO['mean'] = ((pqcd.mu_1(mU_FG)[:,0]*1000.) * 3.0 * dens - \
                         (pres_dict_NLO['mean'] - pres_dict_FG['mean']))[index]
    edens_NLO['lower'] = ((pqcd.mu_1(mU_FG)[:,0]*1000.) * 3.0 * dens \
    - ((pres_dict_NLO['mean'] - pres_dict_FG['mean']) - \
       (pres_dict_NLO['std_dev']-pres_dict_FG['std_dev'])))[index]
    edens_NLO['upper'] = ((pqcd.mu_1(mU_FG)[:,0]*1000.) * 3.0 * dens \
    - ((pres_dict_NLO['mean'] - pres_dict_FG['mean']) + \
       (pres_dict_NLO['std_dev']-pres_dict_FG['std_dev'])))[index]

    # N2LO BCs
    edens_N2LO['mean'] = ((pqcd.mu_2(mU_FG)[:,0]*1000.) * 3.0 * dens - \
                          (pres_dict_N2LO['mean'] - pres_dict_NLO['mean']))[index]

    edens_N2LO['lower'] = ((pqcd.mu_2(mU_FG)[:,0]*1000.) * 3.0 * dens - \
                           ((pres_dict_N2LO['mean'] - pres_dict_NLO['mean']) - \
                            (pres_dict_N2LO['std_dev']-pres_dict_NLO['std_dev'])))[index]

    edens_N2LO['upper'] = ((pqcd.mu_2(mU_FG)[:,0]*1000.) * 3.0 * dens - \
                           ((pres_dict_N2LO['mean'] - pres_dict_NLO['mean']) + \
                            (pres_dict_N2LO['std_dev']-pres_dict_NLO['std_dev'])))[index]

    # add corrections to single out each order
    edens_NLO['mean'] += edens_FG['mean']
    edens_NLO['lower'] += edens_FG['lower']
    edens_NLO['upper'] += edens_FG['upper']

    edens_N2LO['mean'] += edens_NLO['mean']
    edens_N2LO['lower'] += edens_NLO['lower']
    edens_N2LO['upper'] += edens_NLO['upper']

    # combine into dictionary and return
    edens_dict = {
        'FG': edens_FG,
        'NLO': edens_NLO,
        'N2LO': edens_N2LO
    }

    return mu_FG, mU_FG, edens_dict

get_closest_mask(array, values)

Returns a mask corresponding to the locations in array that are closest to values.

array and values must be sorted

Taken from gsum, originally written by J. A. Melendez

Source code in src/neutron_rich_bmm/eos_utils.py

def get_closest_mask(array, values):
    """Returns a mask corresponding to the locations in array that are closest to values.

    array and values must be sorted

    Taken from gsum, originally written by J. A. Melendez
    """
    idxs = np.searchsorted(array, values, side="left")

    # find indexes where previous index is closer
    prev_idx_is_less = ((idxs == len(array))|(np.fabs(values - array[np.maximum(idxs-1, 0)]) \
                                              < np.fabs(values - array[np.minimum(idxs, len(array)-1)])))
    idxs[prev_idx_is_less] -= 1
    return np.isin(np.arange(len(array)), idxs)

get_linear_mask_in_log_space(x, x_min, x_max, log_x_step, base=np.e)

Mask for getting linear data in the log space. Written by J. A. Melendez.

Parameters:

Name	Type	Description	Default
x	ndarray	Array of x values.	required
x_min float)		Min x value.	required
x_max	float	Max x value.	required
log_x_step	float	The step size.	required
base	float	The base of the log we are using. Default is natural log (np.e).	e

Returns:

Type	Description
	The linear mask in the logarithmic space.

Source code in src/neutron_rich_bmm/eos_utils.py

def get_linear_mask_in_log_space(x, x_min, x_max, log_x_step, base=np.e):

    '''
    Mask for getting linear data in the log space. Written by
    J. A. Melendez.

    Parameters:
        x (numpy.ndarray): Array of x values.

        x_min float): Min x value.

        x_max (float): Max x value.

        log_x_step (float): The step size.

        base (float): The base of the log we are using. Default
            is natural log (np.e). 

    Returns:
        The linear mask in the logarithmic space.

    '''
    lin_x = np.arange(
        np.emath.logn(n=base, x=x_min),
        np.emath.logn(n=base, x=x_max),
        log_x_step
    )
    closest = get_closest_mask(np.emath.logn(n=base, x=x), lin_x)
    return (x <= x_max) & (x >= x_min) & closest

gp_data(data_xeft, data_pqcd, cutoff=40, all_orders=True, matter='SNM', kernel='rbf')

Helper function for determining training data from the Chiral EFT and pQCD full training sets. Used for the BMM when a GP is the method of choice.

Parameters:

Name	Type	Description	Default
data_xeft	dict	The dictionary of densities, means, and variances of the Chiral EFT data.	required
data_pqcd	dict	The dictionary of densities, means, and variances of the pQCD data.	required
cutoff	int	The scaled density cutoff we are using for pQCD data.	40
all_orders	bool	Toggle if data is more than one-dimensional. Default is True.	True

Returns:

Name	Type	Description
training_set	dict	The dictionary of selected training data
		concatenated from both EOSs.

Source code in src/neutron_rich_bmm/eos_utils.py

def gp_data(data_xeft, data_pqcd, cutoff=40, all_orders=True, matter='SNM', kernel='rbf'):

    '''
    Helper function for determining training data 
    from the Chiral EFT and pQCD full training sets.
    Used for the BMM when a GP is the method of choice.

    Parameters:
        data_xeft (dict): The dictionary of densities, means, 
            and variances of the Chiral EFT data.

        data_pqcd (dict): The dictionary of densities, means, 
            and variances of the pQCD data.

        cutoff (int): The scaled density cutoff we are using for
            pQCD data.

        all_orders (bool): Toggle if data is more than one-dimensional.
            Default is True.

    Returns:
        training_set (dict): The dictionary of selected training data 
        concatenated from both EOSs.
    '''

    # split into training and testing data
    n_xeft = data_xeft["density"]
    n_pqcd = data_pqcd["density"]

    # toggle
    if all_orders is True:
        p_mean_xeft = data_xeft["mean"][:, -1]
        p_stdv_xeft = data_xeft["std_dev"][:, -1]
        p_cov_xeft = data_xeft["cov"][..., -1]
    else:
        p_mean_xeft = data_xeft["mean"]
        p_stdv_xeft = data_xeft["std_dev"]
        p_cov_xeft = data_xeft["cov"]

    p_mean_pqcd = data_pqcd["mean"][:, -1]
    p_stdv_pqcd = data_pqcd["std_dev"][:, -1]
    p_cov_pqcd = data_pqcd["cov"][..., -1]

    # cut into training and testing sets
    # training data
    n_train_xeft = n_xeft[1::2]
    n_train_pqcd = n_pqcd[1::2]

    chiral_train = {
        'dens': n_train_xeft,
        'mean': p_mean_xeft[1::2],
        'std': p_stdv_xeft[1::2],
        'cov': p_cov_xeft[1::2, 1::2]
    }
    pqcd_train = {
        'dens': n_train_pqcd,
        'mean': p_mean_pqcd[1::2],
        'std': p_stdv_pqcd[1::2],
        'cov': p_cov_pqcd[1::2, 1::2]
    }

    # testing data
    n_test_xeft = n_xeft[::2]
    n_test_pqcd = n_pqcd[::2]

    chiral_test = {
        'dens': n_test_xeft,
        'mean': p_mean_xeft[::2],
        'std': p_stdv_xeft[::2],
        'cov': p_cov_xeft[::2,::2]
    }
    pqcd_test = {
        'dens': n_test_pqcd,
        'mean': p_mean_pqcd[::2],
        'std': p_stdv_pqcd[::2],
        'cov': p_cov_pqcd[::2,::2]
    }

    # store cutoff in terms of density
    pqcd_dens_cutoff = cutoff * 0.164

    # cut the training sets up
    if chiral_train['dens'][-1] > 0.34:
        chiral_cutoff = np.where([chiral_train['dens']>=0.34])[1][0]
    else:
        chiral_cutoff = -1

    chiral_tr = {}
    for key,i in chiral_train.items():
        if chiral_train[key].ndim == 1:
            chiral_tr[key] = chiral_train[key][:chiral_cutoff]
        elif chiral_train[key].ndim == 2:
            chiral_tr[key] = chiral_train[key][:chiral_cutoff, :chiral_cutoff]

    # chiral point selection
    if matter == 'ANM':   
        chiral_tr_final = {}
        for key,i in chiral_tr.items():
            if chiral_tr[key].ndim == 1:
                if kernel == 'rbf':
                    chiral_tr_final[key] = chiral_tr[key][2::20]
                else:
                    chiral_tr_final[key] = chiral_tr[key][::5]  # [2::20]
            elif chiral_tr[key].ndim == 2:
                if kernel == 'rbf':
                    chiral_tr_final[key] = chiral_tr[key][2::20, 2::20]   # [2::20]
                else:
                    chiral_tr_final[key] = chiral_tr[key][::5, ::5]   # [2::20]

    elif matter == 'SNM':
        chiral_tr_final = {}
        for key,i in chiral_tr.items():
            if chiral_tr[key].ndim == 1:
                chiral_tr_final[key] = chiral_tr[key][40::30] #[40::30] 
            elif chiral_tr[key].ndim == 2:
                chiral_tr_final[key] = chiral_tr[key][40::30, 40::30] #[40::30, 40::30]

    elif matter == 'PNM':
        chiral_tr_final = {}
        for key,i in chiral_tr.items():
            if chiral_tr[key].ndim == 1:
                chiral_tr_final[key] = chiral_tr[key][30::9] 
            elif chiral_tr[key].ndim == 2:
                chiral_tr_final[key] = chiral_tr[key][30::9, 30::9]

    print(chiral_tr_final['dens'].shape, chiral_tr_final['mean'].shape, \
          chiral_tr_final['std'].shape, chiral_tr_final['cov'].shape)

    # pqcd point selection
    if pqcd_train['dens'][-1] < pqcd_dens_cutoff:
        pqcd_cutoff = min(pqcd_train['dens'])
    else:
        pqcd_cutoff = np.where([pqcd_train['dens']>=pqcd_dens_cutoff])[1][0]

    pqcd_tr = {}
    for key,i in pqcd_train.items():
        if pqcd_train[key].ndim == 1:
            pqcd_tr[key] = pqcd_train[key][pqcd_cutoff:]
        elif pqcd_train[key].ndim == 2:
            pqcd_tr[key] = pqcd_train[key][pqcd_cutoff:, pqcd_cutoff:]

    # concatenate everything 
    log_space_pqcd = get_linear_mask_in_log_space(pqcd_tr['dens'], pqcd_tr['dens'][0],\
                                                    pqcd_tr['dens'][-1], 0.20, base=np.e)
    pqcd_tr_final = {}
    for key,i in pqcd_tr.items():
        if pqcd_tr[key].ndim == 1:
            if kernel == 'rbf':
                pqcd_tr_final[key] = pqcd_tr[key][::40]
            else:
                pqcd_tr_final[key] = pqcd_tr[key][::5] #[::40] for adding more points for Matern kernel
        elif pqcd_tr[key].ndim == 2:
            if kernel == 'rbf':
                pqcd_tr_final[key] = pqcd_tr[key][::40, ::40]
            else:
                pqcd_tr_final[key] = pqcd_tr[key][::5, ::5] #[::40, ::40]

    print(chiral_tr_final['dens'].shape, chiral_tr_final['mean'].shape, \
          chiral_tr_final['std'].shape, chiral_tr_final['cov'].shape)
    print(pqcd_tr_final['dens'].shape, pqcd_tr_final['mean'].shape, \
          pqcd_tr_final['std'].shape, pqcd_tr_final['cov'].shape)

    # now concatenate into block diagonal matrix
    training_set = {
        'dens' : np.concatenate((chiral_tr_final['dens'], pqcd_tr_final['dens'])),
        'mean' : np.concatenate((chiral_tr_final['mean'], pqcd_tr_final['mean'])),
        'std' : np.concatenate((chiral_tr_final['std'], pqcd_tr_final['std'])),
        'cov' : sp.linalg.block_diag(chiral_tr_final['cov'], pqcd_tr_final['cov'])
    }

    # print the covariance to check
    print('Cov shape:', training_set['cov'].shape)

    return chiral_tr_final, pqcd_tr_final, training_set

pal_eos(kf)

Python version of PAL (Prakash, Ainsworth, Lattimer) EOS. Coupling constants found via the FORTRAN code paleoscc.f90, not included in this function. This function is designed to be used as a mean function in the GP for chiral EFT.

Parameters:

Name	Type	Description	Default
kf	ndarray	The Fermi momentum to be used to calculate PAL for the energy per particle.	required

Returns:

Name	Type	Description
enperpart_kf	ndarray	The energy per particle, in terms of the Fermi momentum.

Source code in src/neutron_rich_bmm/eos_utils.py

def pal_eos(kf):

    '''
    Python version of PAL (Prakash, Ainsworth, Lattimer) EOS. 
    Coupling constants found via the FORTRAN code paleoscc.f90,
    not included in this function. This function is designed
    to be used as a mean function in the GP for chiral EFT.

    Parameters:
        kf (numpy.ndarray): The Fermi momentum to be used to calculate PAL for
            the energy per particle.

    Returns:
        enperpart_kf (numpy.ndarray): The energy per particle, in terms of the
            Fermi momentum.
    '''

    # extract coupling constants from cc dict
    K0 = 260. #cc['K0']       # MeV
    A = -47.83618 #cc['A']    # MeV
    B = 31.01158 #cc['B']     # MeV
    Bp = 0. #cc['Bp']       
    Sig = 1.500259  #cc['Sig'] 

    # if we want Bp != 0.0
#     A = -22.97032
#     B = 22.3410432
#     Bp = 0.2
#     Sig = 1.9999723

    # other constants
    hc = 197.33     # MeV fm
    n0 = 0.164      # fm^-3
    mn = 939.       # MeV
    kf0 = (1.5*np.pi**2.*0.164)**(1./3.)    # fm^1
    ef0 = (hc*kf0)**2./2./939.              # MeV
    sufac = (2.**(2./3.)-1.)*0.6*ef0        # MeV
    s0 = 30.                                # MeV

    # other coupling constants
    C1 = -83.841                            # MeV
    C2 = 22.999                             # MeV
    Lambda1 = 1.5*kf0                       # fm^-1
    Lambda2 = 3.*kf0                        # fm^-1

    # conversion from kf to n to solve that problem
    n = 2.0 * kf**3.0 / (3.0 * np.pi**2.0)          # fm^-3

    # write it as E/A first
    one = mn * n0 * (n/n0) + (3.0/5.0)*ef0*n0*(n/n0)**(5.0/3.0)       # MeV/fm^3
    two = 0.5*A*n0*(n/n0)**2.0 + (B*n0*(n/n0)**(Sig+1.0))/(1.0 + Bp * (n/n0)**(Sig - 1.0))   # MeV/fm^3
    sum_1 = C1 * (Lambda1/kf0)**3.0 * ((Lambda1/kf) - np.arctan(Lambda1/kf))  # MeV
    sum_2 = C2 * (Lambda2/kf0)**3.0 * ((Lambda2/kf) - np.arctan(Lambda2/kf))  # MeV 
    three = 3.0 * n0 * (n/n0) * (sum_1 + sum_2)                               # MeV/fm^3

    eps_kf = one + two + three     # MeV/fm^3

    # convert to E/A from eps
    enperpart_kf = eps_kf * n    # MeV

    # now calculate pressure using differentiation
 #   derivEA = np.gradient(enperpart_kf, kf)
 #   pressure_kf = (n * kf / 3.0) * np.asarray(derivEA)

    return enperpart_kf

pressure_pal_eos(kf)

The PAL EOS pressure calculation.

Parameters:

Name	Type	Description	Default
kf	ndarray	The Fermi momentum .	required

Returns:

Name	Type	Description
pressure_kf	ndarray	The pressure in terms of the Fermi momentum.

Source code in src/neutron_rich_bmm/eos_utils.py

def pressure_pal_eos(kf):

    '''
    The PAL EOS pressure calculation.

    Parameters:
        kf (numpy.ndarray): The Fermi momentum  .

    Returns:
        pressure_kf (numpy.ndarray): The pressure in terms of the Fermi
            momentum.
    '''

    # calculate n first again (so we don't have to pass it in)
    n = 2.0 * kf**3.0 / (3.0 * np.pi**2.0)

    # now calculate pressure using differentiation
    derivEA = ndt.Derivative(pal_eos, step=1e-4, method='central')
    pressure_kf = (n * kf / 3.0) * np.asarray(derivEA(kf))

    return pressure_kf

speed_of_sound(dens, pressure, edens=None, sat=False, bounds=68, integrate='forward', sampled=False)

Function to evaluate the speed of sound of a system given the pressure, number density, and initial parameters for the energy density integration.

Parameters:

Name	Type	Description	Default
dens	numpy 1d array	The number density of the system.	required
pressure	dict	The dictonary of pressure means and standard deviations from the system.	required
edens	dict	The dictionary of energy density means and standard deviations for a specific starting point in density.	None
sat	bool	Starting at saturation density (0.16 fm^-3) or not. Default is False.	False
integrate	str	Decision to integrate forward or backward. Default is 'forward'.	'forward'
sampled	bool	If using samples from the speed of sound, run the std and mean using nanmean and nanstd from numpy instead of computing envelopes. Default is 'False'.	False

Returns:

Name	Type	Description
cs2	dict	The dictonary of results for the speed of sound (calculated using 1\mu dP/dn) and the lower and upper bounds of it at one sigma, returned when sampled is True.
edens_full	dict	The energy density dictionary of means and variances returned when sampled is True.
dens_arr	ndarray	The densities corresponding to the speed of sound calculation (if sat is True, this will reflect from saturation up), returned when sampled is False.
cs2_log	dict	The dict of speed of sound values from using the n * dlog(mu)/dn method. Returned when sampled is False.
edens_int	dict	The dict of energy densities, returned when sampled is False.
mu_dict	dict	The dict of chemical potential values, returned when sampled is False.

Source code in src/neutron_rich_bmm/eos_utils.py

def speed_of_sound(dens, pressure, edens=None, sat=False, bounds=68, integrate='forward', sampled=False):

    '''
    Function to evaluate the speed of sound of
    a system given the pressure, number density,
    and initial parameters for the energy
    density integration. 

    Parameters:
        dens (numpy 1d array): The number density of the system.

        pressure (dict): The dictonary of pressure means
            and standard deviations from the system.

        edens (dict): The dictionary of energy density 
            means and standard deviations for a 
            specific starting point in density.

        sat (bool): Starting at saturation density 
            (0.16 fm^-3) or not. Default is False.

        integrate (str): Decision to integrate forward or backward.
            Default is 'forward'.

        sampled (bool): If using samples from the speed of sound, run
            the std and mean using nanmean and nanstd from
            numpy instead of computing envelopes.
            Default is 'False'. 

    Returns:
        cs2 (dict): The dictonary of results for the 
            speed of sound (calculated using 1\mu dP/dn)
            and the lower and upper bounds of it at 
            one sigma, returned when sampled is True.

        edens_full (dict): The energy density dictionary of means and
            variances returned when sampled is True.

        dens_arr (numpy.ndarray): The densities corresponding to the 
            speed of sound calculation (if sat is True, 
            this will reflect from saturation up), returned
            when sampled is False.

        cs2_log (dict): The dict of speed of sound values from using
            the n * dlog(mu)/dn method. Returned when 
            sampled is False.

        edens_int (dict): The dict of energy densities, returned when 
            sampled is False.

        mu_dict (dict): The dict of chemical potential values, returned
            when sampled is False.
    '''

    # check for saturation point integration
    if sat is True:
        dens_arr = np.linspace(0.164, 16.4, 1200)
    else:
        dens_arr = dens

    # using samples
    if sampled is True:
        pres = np.asarray(pressure['samples'])   # (nB, n_samples) shape
        edens_0 = edens['samples']

        # huge list for all sampled curves
        edens_full = []

        # collect the function together
        dn = dens[1] - dens[0]    # equally spaced

        # interpolation and integration for each sample 
        for i in range(len(pres.T)):

            # empty list for storing (re-initialize to dump old data)
            en_samples = np.zeros(len(pres))

            # outer term (not changing with n)
            outer = (edens_0[i]/dens[-1])  # adding change of integration constant w/each sample

            # running integration backward from pQCD
            for j in range(len(dens)):

                # Simpson's Rule integration
                en_samples[j] = dens[j] * (outer - scint.simps((pres[j:, i]/dens[j:]**2.0), dens[j:]))

            edens_full.append(en_samples)   # shape (n_samples, nB)

        # now calculate chemical potential and derivative
        mu_samples = np.asarray([((np.asarray(edens_full)[i,:] + pres[:,i]))/dens for \
                                 i in range(len(edens_full))])   # samples, nB

        # get the results using 1/mu dP/dn instead (more stable)
        dpdn_samples = np.gradient(pres, dn, axis=0, edge_order=2)

        cs2_samples = np.asarray([(mu_samples[i,:])**(-1.0) * dpdn_samples[:,i] \
                                  for i in range(len(edens_full))])

        # get mean, std_dev estimations out, store and return
        cs2_mean = np.nanmean(cs2_samples, axis=0)  # run over samples
        cs2_std = np.nanstd(cs2_samples, axis=0)

        cs2 = {
            'mean': cs2_mean,
            'std': cs2_std,
            'samples': cs2_samples
        }

        return cs2, edens_full

    # first, make bounds flexible
    zscore = stats.norm.ppf((1.0 + (bounds/100.0))/2.0)

    # extract the necessary information
    p_mean = pressure['mean']
    p_low = pressure['mean'] - zscore*pressure['std_dev']
    p_high = pressure['mean'] + zscore*pressure['std_dev']

    # extract the parameters for edens (for pqcd these will be full curves)
    e_mean = edens['mean'][0]
    e_low = edens['lower'][0]
    e_high = edens['upper'][0]

    # calculate the interpolants
    p_mean_interp = interp1d(dens, (p_mean), kind='cubic', \
                            fill_value='extrapolate')
    p_lower_interp = interp1d(dens, (p_low), kind='cubic', \
                            fill_value='extrapolate')
    p_upper_interp = interp1d(dens, (p_high), kind='cubic', \
                            fill_value='extrapolate')

    # define internal functions for integration
    def pres_mean(n):
        return p_mean_interp(n) / (n)**2.0
    def pres_lower(n):
        return p_lower_interp(n) / (n)**2.0
    def pres_upper(n):
        return p_upper_interp(n) / (n)**2.0

    # perform integration
    en_mean = []
    en_lower = []
    en_upper = []

    # integrating forwards
    if integrate == 'forward':

        for n in dens_arr:
            en_mean.append(n*(e_mean/dens_arr[0] + \
                            scint.quad(lambda x : pres_mean(x), dens_arr[0], n, epsabs=1e-10, epsrel=1e-10)[0]))

            en_lower.append(n*(e_low/dens_arr[0] + \
                            scint.quad(lambda x : pres_lower(x), dens_arr[0], n, epsabs=1e-10, epsrel=1e-10)[0]))
            en_upper.append(n*(e_high/dens_arr[0] + \
                            scint.quad(lambda x : pres_upper(x), dens_arr[0], n, epsabs=1e-10, epsrel=1e-10)[0]))

    # try integrating backwards
    elif integrate == 'backward':

        tol = 1e-08

        for n in dens_arr:
            en_mean.append(n*(e_mean/dens_arr[-1] - \
                            scint.quad(lambda x : pres_mean(x), n, dens_arr[-1], epsabs=tol, epsrel=tol)[0]))
            en_lower.append(n*(e_low/dens_arr[-1] - \
                            scint.quad(lambda x : pres_lower(x), n, dens_arr[-1], epsabs=tol, epsrel=tol)[0]))
            en_upper.append(n*(e_high/dens_arr[-1] - \
                            scint.quad(lambda x : pres_upper(x), n, dens_arr[-1], epsabs=tol, epsrel=tol)[0]))

    # dict of energy densities
    edens_int = {
        'mean': np.asarray(en_mean),
        'lower': np.asarray(en_lower),
        'upper': np.asarray(en_upper)
    }

    # calculate deriv of pressure
    dpdn_mean = ndt.Derivative(p_mean_interp, step=1e-6, method='central')
    dpdn_lower = ndt.Derivative(p_lower_interp, step=1e-6, method='central')
    dpdn_upper = ndt.Derivative(p_upper_interp, step=1e-6, method='central')

    # calculate the chemical potential
    mu_mean = (en_mean + p_mean_interp(dens_arr))/dens_arr
    mu_lower = (en_lower + p_lower_interp(dens_arr))/dens_arr
    mu_upper = (en_upper + p_upper_interp(dens_arr))/dens_arr

    # calculate the log of the chemical potential
    log_mu_mean = np.log(mu_mean)
    log_mu_lower = np.log(mu_lower)
    log_mu_upper = np.log(mu_upper)

    # calculate speed of sound using chemical potential
    # at desired density array
    cs2_mu_mean = dpdn_mean(dens_arr) / mu_mean
    cs2_mu_lower = dpdn_lower(dens_arr) / mu_upper
    cs2_mu_upper = dpdn_upper(dens_arr) / mu_lower

    # calculate speed of sound using log(mu)
    # at desired density array
    cs2_log_mean = dens_arr * np.gradient(log_mu_mean, dens_arr, edge_order=2)
    cs2_log_lower = dens_arr * np.gradient(log_mu_lower, dens_arr, edge_order=2)
    cs2_log_upper = dens_arr * np.gradient(log_mu_upper, dens_arr, edge_order=2)

    # collect into dict and return
    cs2 = {
        'mean' : cs2_mu_mean,
        'lower' : cs2_mu_lower,
        'upper' : cs2_mu_upper
    }

    # collect log method and return
    cs2_log = {
        'mean': cs2_log_mean,
        'lower': cs2_log_lower,
        'upper': cs2_log_upper
    }

    # collect mu into dict and return
    mu_dict = {
        'mean': mu_mean,
        'lower': mu_lower,
        'upper': mu_upper
    }

    return dens_arr, cs2, cs2_log, edens_int, mu_dict

causality_stability(cs2, edens, pressure)

Function to eliminate any draws from the EOS that are unstable or acausal.

Example

causality_stability(cs2=array, edens=array, pressure=array)

Parameters:

Name	Type	Description	Default
cs2	array	The array of speed of sound draws. Shape should be [density, draw].	required
edens	array	The array of corresponding energy densities, in the shape [density, draw].	required
pressure	array	The array of corresponding pressures in the shape of [density, draw].	required

Returns:

Name	Type	Description
cs2_reduced	array	The 2D array of reduced samples for the speed of sound.
edens_reduced	array	The 2D array of reduced samples for the energy density.
pres_reduced	array	The 2D array of reduced samples for the pressure.

Source code in src/neutron_rich_bmm/tov_utils.py

def causality_stability(cs2, edens, pressure):

    r"""
    Function to eliminate any draws from the EOS that are unstable
    or acausal.

    Example:
        causality_stability(cs2=array, edens=array, pressure=array)

    Parameters:
        cs2 (array): The array of speed of sound draws. Shape should
            be [density, draw].

        edens (array): The array of corresponding energy densities,
            in the shape [density, draw].

        pressure (array): The array of corresponding pressures in the
            shape of [density, draw]. 

    Returns:
        cs2_reduced (array): The 2D array of reduced samples for 
            the speed of sound.

        edens_reduced (array): The 2D array of reduced samples for 
            the energy density.

        pres_reduced (array): The 2D array of reduced samples 
            for the pressure.
    """

    # check if reshaping is needed (probably unnecessary)
    if cs2.ndim == 1:
        cs2 = cs2.reshape(-1,1)

    # run over the draws (must be given in [density, draws]!)
    ind_list = []
    for i in range(len(cs2)):
        for j in range(len(cs2.T)):
            if cs2[i,j] > 1.0 or cs2[i,j] < 0.0:   # could do max > 1.0 or min < 0.0 or something
                ind_list.append(j)   # cut out the draws (make sure this works right)

    # cut these out of the draws in the TOV results below
    cs2_reduced = np.delete(cs2, ind_list, axis=1)
    pressure_reduced = np.delete(pressure, ind_list, axis=1)
    edens_reduced = np.delete(edens, ind_list, axis=1)

    return cs2_reduced, edens_reduced, pressure_reduced

tov_data(edens_full, pres_dict, cs2_data=None, save=False, filepath=None)

Function to collect the data needed for the TOV input and to organize it such that it can be read correctly by the TOV solver we are using. It also attaches the low-density EOS information up to 0.5n0 from Bethe, Pethick, and Sutherland, and Negele and Vautherin's work.

Example

tov_data(edens_full=array, pres_dict=dict, cs2_data=array)

Parameters:

Name	Type	Description	Default
edens_full	array	The array of energy densities in MeV/fm^3.	required
pres_dict	dict	The dictionary of pressures and their standard deviations in MeV/fm^3, as well as the number density in fm^-3.	required
cs2_data	array	The array of corresponding speed of sound values in dimensionless units. This is optional, as these values are only needed if the TOV solver is selected to calculate tidal deformability.	None
save	bool	If the user wishes to save the data file, supply 'True'.	False
filepath	str	If the user wishes to save the data file, supply a proper file path, including any external folders.	None

Returns:

Name	Type	Description
tov_dict	dict	The full dictionary of dens, edens, pres, and cs2 (optionally) to be used in the TOV solver.

Source code in src/neutron_rich_bmm/tov_utils.py

def tov_data(edens_full, pres_dict, cs2_data=None, save=False, filepath=None):

    r"""
    Function to collect the data needed for the TOV input and to organize
    it such that it can be read correctly by the TOV solver we are using. 
    It also attaches the low-density EOS information up to 0.5n0 from Bethe,
    Pethick, and Sutherland, and Negele and Vautherin's work. 

    Example:
        tov_data(edens_full=array, pres_dict=dict, cs2_data=array)

    Parameters:
        edens_full (array): The array of energy densities in MeV/fm^3.

        pres_dict (dict): The dictionary of pressures and their standard deviations
            in MeV/fm^3, as well as the number density in fm^-3.

        cs2_data (array): The array of corresponding speed of sound values in 
            dimensionless units. This is optional, as these values are only 
            needed if the TOV solver is selected to calculate tidal deformability.

        save (bool): If the user wishes to save the data file, supply 'True'.

        filepath (str): If the user wishes to save the data file, supply a proper
            file path, including any external folders.

    Returns:
        tov_dict (dict): The full dictionary of dens, edens, pres, and cs2 
            (optionally) to be used in the TOV solver. 
    """

    # number of samples
    samples = len(pres_dict['samples'].T)

    # organize the samples into the correct pairing (same style as a file would be)
    low_den_file = np.loadtxt("../data/NSM_data/MFT_ns6p.dat", skiprows=1)

    # turn edens_full into an array  [density, draws] is needed
#     edens = np.asarray(edens_full)  #.T
#     cs2 = np.asarray(cs2_data['samples'])  #.T

    # rename dicts
    edens = edens_full
    cs2 = cs2_data

    # get organized data
    tov_index = (np.where([pres_dict['dens'][i] <= 0.08 for i in range(len(pres_dict['dens']))])[0][-1] + 1)
    edens_final = edens[tov_index:, :]
    gp_final = np.asarray([pres_dict['samples'][tov_index:, i] for i in range(samples)]).T
    density_final = pres_dict['dens'][tov_index:]

    if cs2_data is not None:
        cs2_final = cs2[tov_index:, :]

    # run through and append the low density data to these arrays and then save to file
    edens_tov = np.asarray([np.concatenate((low_den_file[::-1,0], edens_final[:,i]))\
                            for i in range(samples)]).T
    pres_tov = np.asarray([np.concatenate((low_den_file[::-1,1], gp_final[:,i])) \
                            for i in range(samples)]).T
    dens_tov = np.concatenate((low_den_file[::-1,2], density_final)).reshape(-1,1)

    if cs2_data is not None:
        # solve for the speed of sound using edens and pres at low densities
        dpdeps = np.gradient(low_den_file[::-1,1], low_den_file[::-1,0], edge_order=2)  # might need tweaking
        cs2_tov = np.asarray([np.concatenate((dpdeps, cs2_final[:,i])) for i in range(samples)]).T
    else:
        cs2_tov = np.zeros(len(density_final) + len(low_den_file[:,0]))

    # eliminate the information above the core of the star (assuming 10 times saturation high enough)
    index_10 = np.where([dens_tov[i] <= 1.64 for i in range(len(dens_tov))])[0][-1] + 1

    # cut the arrays
    dens_end = dens_tov  #[:index_10]
    edens_end = edens_tov  #[:index_10]
    pres_end = pres_tov  #[:index_10]
    cs2_end = cs2_tov  #[:index_10]

    # save the result if desired
    if save is True:
        np.savez(filepath, density=dens_end, \
                edens=edens_end, pres=pres_end, cs2=cs2_end)
    else: 
        print('Data not saved.')

    # make a dict out of the results and return
    tov_dict = {
        'dens': dens_end,
        'edens': edens_end,
        'pres': pres_end,
        'cs2': cs2_end
    }

    print('''I'm done!''')

    return tov_dict

Chiral EFT EOS code

Chiral_model(density=None, Lambda=500, high_density=True)

Parameters:

Name	Type	Description	Default
density	array	The density specified for calculation of the chiral EOS. Default is None, which will make the chiral EOS default to the original data density range.	None
Lambda	int	The value of the cutoff chosen. Can be either 450 MeV or 500 MeV.	500
high_density	bool	Whether we want to use high-density data or not. Default is True. Sets up the data for interpolation in the chiral EOS class.	True

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/chiral_model.py

def __init__(self, density=None, Lambda=500, high_density=True):

    '''
    Model class to give Taweret to mix. 

    Parameters:
        density (numpy.array): The density specified for calculation of the chiral
            EOS. Default is None, which will make the chiral EOS
            default to the original data density range. 

        Lambda (int): The value of the cutoff chosen. Can be either
            450 MeV or 500 MeV. 

        high_density (bool): Whether we want to use high-density data or not.
            Default is True. Sets up the data for interpolation
            in the chiral EOS class. 

    Returns:
        None.
    '''

    # set class variables
    self.Lambda = Lambda
    self.density = density

    # call Chiral class and instantiate
    self.chiral = Chiral(density_input=self.density, Lambda=self.Lambda, \
                         high_density=high_density) 

    return None 

evaluate(input_space=None, N3LO=True, scaled=True, extend=False)

Returns the mean and standard deviation of the chiral EFT EOS in terms of pressure wrt baryon chemical potential. The mean is calculated from the chiral EFT EOS formalism of C. Drischler et al. (2021). The standard deviation is calculated via the truncation error models in the gsum package, used in the chiral EOS paper.

Parameters:

Name	Type	Description	Default
input_space	array	The input space array. Not actually necessary for this function but necessary for Taweret.	None
N3LO	bool	If True, returns only the N3LO results for mean and std_dev. Otherwise will return all results up to and through N3LO.	True
scaled	bool	If the data is scaled, then this is True. Else, it is False. Default is True.	True
extend	bool	Extends the data to higher truncation values. Default is False.	False

Returns:

Name	Type	Description
mean	numpy.1darray	The mean of the pressure.
std_dev	numpy.1darray	The standard deviation of the pressure.

Source code in src/neutron_rich_bmm/chiral_model.py

    def evaluate(self, input_space=None, N3LO=True, scaled=True, extend=False):

        '''
        Returns the mean and standard deviation of the chiral EFT EOS
        in terms of pressure wrt baryon chemical potential. The mean is calculated 
        from the chiral EFT EOS formalism of C. Drischler et al. (2021). The 
        standard deviation is calculated via the truncation error models in the gsum 
        package, used in the chiral EOS paper.

        Parameters:
            input_space (numpy.array): The input space array. Not actually necessary 
                for this function but necessary for Taweret. 

            N3LO (bool): If True, returns only the N3LO results for mean and std_dev.
                Otherwise will return all results up to and through N3LO.

            scaled (bool): If the data is scaled, then this is True. Else, it is False.
                Default is True.

            extend (bool): Extends the data to higher truncation values. Default is
                False.

        Returns:
            mean (numpy.1darray): The mean of the pressure. 

            std_dev (numpy.1darray): The standard deviation of the pressure.
        '''

        # correct for input_space (for now)
        if input_space is not None:
            input_space = None

        # set up the data containers for interpolation
        self.obs_neutron, self.obs_nuclear, self.obs_sym_energy = \
            self.chiral.data_interpolation(density_int=None, kf_s_int=None, extend=extend)

        # energy per particle
        self.energies_s, self.energy_s_stds = \
            self.chiral.energy_per_particle(add_rest_mass=True, orders='all')

        # chemical potential (calls pressure and eps already)
        self.mu_s, self.mu_s_stds = \
            self.chiral.chemical_potential(method=1, add_rest_mass=True)

        # inversion for n(mu)
#         self.density_mu_N3LO, self.mu_N3LO_array, self.kf_N3LO = \
#             self.chiral.inversion()

        # call the container again with n(mu) and kf(n(mu))
#         self.obs_neutron, self.obs_nuclear, self.obs_sym_energy = \
#             self.chiral.data_interpolation(density_int=self.density_mu_N3LO, \
#                                            kf_s_int=self.kf_N3LO)

        # calculate pressure
        self.pressures_s, self.pressure_s_stds, self.pressure_s_cov = self.chiral.pressure(\
            orders='all')

        # N3LO results only
        if N3LO is True:
            mean = self.pressures_s[:,3]
            std_dev = self.pressure_s_stds[:,3]

            if scaled is True:
                pressure_free = np.zeros([len(self.density), 4])
                hbarc = 197.327 # MeV fm
                for i in range(4):
                    pressure_free[:,i] = ((1.0/(2.0*np.square(np.pi)))*(self.chiral.kf_s_all**4.0))*(hbarc)
                mean = mean/pressure_free[:,3]
                std_dev = std_dev/pressure_free[:,3]

            else:
                mean = mean
                std_dev = std_dev

#             for i in range(len(mean)):
#                 if mean[i] <= mean[i-1]:
#                     std_dev[i] = 1e10

        else:
            mean = self.pressures_s
            std_dev = self.pressure_s_stds

        return mean, std_dev

log_likelihood_elementwise()

Calculates the log likelihood for the model. Not needed for this model.

Source code in src/neutron_rich_bmm/chiral_model.py

def log_likelihood_elementwise(self):

    '''
    Calculates the log likelihood for the model.
    Not needed for this model.
    '''

    return None

set_prior()

Sets the prior on the model. Not needed for this model.

Source code in src/neutron_rich_bmm/chiral_model.py

def set_prior(self):

    '''
    Sets the prior on the model. Not needed for
    this model. 
    '''

    return None

Chiral(density_input=None, Lambda=500, high_density=True)

function imports the necessary data to begin the calculation, and sets up the density and data for interpolation. Choice of Lambda (450, 500) is also made here.

Parameters:

Name	Type	Description	Default
density_input	linspace	The density regime for the chiral EOS. Default is set below in the code.	None
Lambda	int	The desired cutoff to use. Default is 500; either 450 or 500 MeV may currently be chosen.	500
high_density	bool	Decision to use high density data or not. Default is True.	True

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/chiral.py

def __init__(self, density_input=None, Lambda=500, high_density=True):

    '''
    The class for the chiral EFT EOS model. Will be used
    in the Taweret model class for the chiral EOS. This init
    function imports the necessary data to begin the calculation, and
    sets up the density and data for interpolation. Choice of Lambda 
    (450, 500) is also made here. 

    Parameters:
        density_input (numpy.linspace): The density regime for the chiral 
            EOS. Default is set below in the code. 

        Lambda (int): The desired cutoff to use. Default is 500; either
            450 or 500 MeV may currently be chosen. 

        high_density (bool): Decision to use high density data or not.
            Default is True. 

    Returns:
        None.
    '''

    # set cutoff
    self.Lambda = Lambda

    if high_density is True:
        filename = '../nuclear-matter-convergence/data/all_matter_data_high_density.csv'
    else:
        filename = '../nuclear-matter-convergence/data/all_matter_data.csv'

    # import the data
    data = InputData(filename, Lambda)

    # set orders and breakdown scale
    self.orders = np.array([0, 2, 3, 4])  
    self.breakdown = 600 # MeV  

    # set up all of the necessary parameters
    self.kf_n = data.kf_n
    self.Kf_n = data.Kf_n
    self.kf_s = data.kf_s
    self.Kf_s = data.Kf_s
    self.kf_d = data.kf_avg
    self.Kf_d = data.Kf_avg

    self.density = data.density

    self.ref_2bf = data.ref_2bf
    self.ref_n_3bf = data.ref_n_3bf
    self.ref_s_3bf = data.ref_s_3bf
    self.ref_d_3bf = data.ref_avg_3bf

    self.y_s_2_plus_3bf = data.y_s_2_plus_3bf
    self.y_n_2_plus_3bf = data.y_n_2_plus_3bf
    self.y_d_2_plus_3bf = data.y_d_2_plus_3bf

    self.y_s_2bf = data.y_s_2bf
    self.y_n_2bf = data.y_n_2bf
    self.y_d_2bf = data.y_d_2bf

    self.y_s_3bf = data.y_s_3bf
    self.y_n_3bf = data.y_n_3bf
    self.y_d_3bf = data.y_d_3bf

    self.fit_n2lo = data.fit_n2lo
    self.fit_n3lo = data.fit_n3lo

    self.data = data 

    # set up density
    if density_input is not None:
        self.density_all = density_input
    else:
        self.density_all = np.arange(self.density[0], self.density[-1], 0.005)

    self.N_all = len(self.density_all)

    #neutron matter
    self.kf_n_all = fermi_momentum(self.density_all, degeneracy=2)
    self.Kf_n_all = self.kf_n_all[:, None]

    #symmetric matter
    self.kf_s_all = fermi_momentum(self.density_all, degeneracy=4)
    self.Kf_s_all = self.kf_s_all[:, None]

    #symmetry energy
    self.kf_d_all = (self.kf_n_all + self.kf_s_all) / 2.
    self.Kf_d_all = self.kf_d_all[:, None]

    # error band setup for MBPT calculations
    min_uncertainty = 0.02  # Twenty keV
    uncertainty_factor = 0.001  # 0.1%

    self.err_y_n = np.abs(self.y_n_2_plus_3bf[:, -1]) * uncertainty_factor
    self.err_y_n[np.abs(self.err_y_n) < min_uncertainty] = min_uncertainty

    self.err_y_s = np.abs(self.y_s_2_plus_3bf[:, -1]) * uncertainty_factor
    self.err_y_s[np.abs(self.err_y_s) < min_uncertainty] = min_uncertainty

    self.err_y_d = np.sqrt(self.err_y_n**2 + self.err_y_s**2)

    # momenta and std dev setup
    self.kf0_n = fermi_momentum(0.164, 2)
    self.kf0_s = fermi_momentum(0.164, 4)

    self.ref_neutron = 16 / self.kf0_n**2
    self.ref_nuclear = 16 / self.kf0_s**2

    if self.Lambda == 500:
        self.std_neutron = 1.00
        self.ls_neutron = 0.973
        self.std_nuclear = 2.95
        self.ls_nuclear = 0.484

        self.rho = None   # letting it calculate rho
        self.ls_n_sym = self.ls_neutron
        self.ls_s_sym = self.ls_nuclear

    elif self.Lambda == 450:
        self.std_neutron = 0.8684060649936118
        self.ls_neutron = 0.7631421388401067
        self.std_nuclear = 2.6146499024837073
        self.ls_nuclear = 0.46603268529311087

        self.rho = 0.95
        self.ls_n_sym = (self.ls_neutron + self.ls_nuclear) / 2
        self.ls_s_sym = None

    return None

chemical_potential(method=1, add_rest_mass=False)

Calculation of the chemical potential depending on which type is desired: either type=1 ((P+eps)/n) or type=2 (deps/dn). Both should be equivalent if thermodynamic consistency is preserved.

Parameters:

Name	Type	Description	Default
method	int	The type of calculation to obtain the chemical potential. Possibilities include: 1 = (P+eps)/n and 2 = d(eps)/dn. Default is 1.	1
add_rest_mass	bool	Adds the rest mass to the energy per particle if desired. Default is False.	False

Returns:

Type	Description
	self.mu_s, self.mu_s_stds (numpy.ndarray): The chemical potential and std dev for the chiral EOS, wrt n.

Source code in src/neutron_rich_bmm/chiral.py

def chemical_potential(self, method=1, add_rest_mass=False):

    '''
    Calculation of the chemical potential depending on which
    type is desired: either type=1 ((P+eps)/n) or type=2 (deps/dn).
    Both should be equivalent if thermodynamic consistency is 
    preserved. 

    Parameters:
        method (int): The type of calculation to obtain the chemical 
            potential. Possibilities include: 1 = (P+eps)/n and 
            2 = d(eps)/dn. Default is 1. 

        add_rest_mass (bool): Adds the rest mass to the energy per 
            particle if desired. Default is False. 

    Returns:
        self.mu_s, self.mu_s_stds (numpy.ndarray): The chemical 
            potential and std dev for the chiral EOS, wrt n. 
    '''

    # call the energy_density function and tell it rest mass choice
    self.energy_density_s, self.energy_density_s_stds = \
        self.energy_dens(add_rest_mass=add_rest_mass)

    # type 1
    if method == 1:

        # call the pressure function
        self.pressures_s, self.pressure_s_stds, self.pressure_s_cov = self.pressure()

        ti_upper = self.pressures_s + self.energy_density_s
        ti_upper_stds = self.pressure_s_stds + self.energy_density_s_stds
        mu_s = np.zeros([len(ti_upper), 4])
        mu_s_stds = np.zeros([len(ti_upper), 4])

        for i in range(4):
            mu_s[:,i] = ti_upper[:,i]/self.density_all
            mu_s_stds[:,i] = ti_upper_stds[:,i]/self.density_all

        # class variable
        self.mu_s = mu_s 
        self.mu_s_stds = mu_s_stds

        return self.mu_s, self.mu_s_stds

    elif method == 2:

        # compute the derivative of the energy density
        k_F = (3.0/2.0 * np.pi**2.0 * self.density_all)**(1.0/3.0)

        dE = np.zeros([len(self.density_all), 4])
        dE_cov = np.zeros([len(dE), len(dE)])
        dE_std = np.zeros([len(dE), 4])

        mu_s_deriv = np.zeros([len(self.density_all), 4])
        mu_s_deriv_stds = np.zeros([len(self.density_all), 4])

        for i, n in enumerate(self.orders):
            dE[:,i] = self.obs_nuclear.get_pred(order=n, deriv=1) 
            dE_cov = self.obs_nuclear.get_cov(order=n, deriv1=1, deriv2=1)
            dE_std[:,i] = np.sqrt(np.diag(dE_cov))

        for i in range(4):
            mu_s_deriv[:,i] = (k_F * dE[:,i])/3.0 + \
                (self.energy_density_s[:,i]/self.density_all)
            mu_s_deriv_stds[:,i] = (k_F * dE_std[:,i])/3.0 + \
                  (self.energy_density_s_stds[:,i]/self.density_all)

        # class variable
        self.mu_s = mu_s_deriv
        self.mu_s_stds = mu_s_deriv_stds

        return self.mu_s, self.mu_s_stds

    else:
        raise ValueError('Method must be either 1 or 2.')

data_interpolation(density_int=None, kf_s_int=None, extend=False)

The interpolation of the data using the observable containers from nuclear-matter-convergence for neutron matter, symmetric nuclear matter, and the symmetry energy calculation.

Parameters:

Name	Type	Description	Default
density_int	numpy.1darray	The interpolated number density values from the mu(n) inversion. Default is None; this will leave the container to use the original interpolated density_all.	None
kf_s_int	numpy.1darray	The interpolated Fermi momenta for SNM. This will be replaced when calculating the inverted density and chemical potential. Default is None; this will leave the container to use the original interpolated kf_s.	None
extend	bool	Whether or not to extend the data using the mean function estimation and extending the truncation error. Default is False.	False

Returns:

Type	Description
	self.obs_neutron, self.obs_nuclear, self.obs_sym_energy (objects): The objects that contain the important interpolation information for the chiral EFT EOS for PNM, SNM, and the symmetry energy.

Source code in src/neutron_rich_bmm/chiral.py

def data_interpolation(self, density_int=None, kf_s_int=None, extend=False):

    '''
    The interpolation of the data using the observable
    containers from nuclear-matter-convergence for neutron matter,
    symmetric nuclear matter, and the symmetry energy calculation.

    Parameters:
        density_int (numpy.1darray): The interpolated number density values 
            from the mu(n) inversion.
            Default is None; this will leave the container to use the 
            original interpolated density_all. 

        kf_s_int (numpy.1darray): The interpolated Fermi momenta for SNM. This will be replaced
            when calculating the inverted density and chemical potential.
            Default is None; this will leave the container to use the
            original interpolated kf_s. 

        extend (bool): Whether or not to extend the data using the mean 
            function estimation and extending the truncation error. Default is False.

    Returns:
        self.obs_neutron, self.obs_nuclear, self.obs_sym_energy (objects):
            The objects that contain the important interpolation information
            for the chiral EFT EOS for PNM, SNM, and the symmetry energy. 
    '''

    setup_time_start = time.time()
    verbose = True

    # alter the density and kf_s in the obs_nuclear container
    # if the routine is called for the inversion calculation
    if density_int is not None:
        self.dens_interp = density_int
    else: 
        self.dens_interp = self.density_all 

    if kf_s_int is not None:
        self.kf_s_all = kf_s_int

    print('Setting up neutron matter...', flush=True)
    self.obs_neutron = ObservableContainer(
        density=self.density,
        kf=self.kf_n,
        y=self.y_n_2_plus_3bf,
        orders=self.orders,
        density_interp=self.density_all,
        kf_interp=self.kf_n_all,
        std=self.std_neutron,
        ls=self.ls_neutron,
        ref=self.ref_neutron,
        breakdown=self.breakdown,
        err_y=self.err_y_n,
        include_3bf=False,
        derivs=[0, 1, 2],
        verbose=verbose,
        extend=extend,
    )

    print('Setting up nuclear matter...', flush=True)
    self.obs_nuclear = ObservableContainer(
        density=self.density,
        kf=self.kf_s,
        y=self.y_s_2_plus_3bf,
        orders=self.orders,
        density_interp=self.dens_interp, 
        kf_interp=self.kf_s_all,
        std=self.std_nuclear,
        ls=self.ls_nuclear,
        ref=self.ref_nuclear,
        breakdown=self.breakdown,
        err_y=self.err_y_s,
        include_3bf=False,
        derivs=[0, 1, 2],
        verbose=verbose,
        extend=extend,
    )

    print('Setting up symmetry energy...', flush=True)
    self.obs_sym_energy = SymmetryEnergyContainer(
        density=self.density,
        y=self.y_d_2_plus_3bf,
        orders=self.orders,
        density_interp=self.density_all,
        std_n=self.std_neutron,
        ls_n=self.ls_n_sym,
        std_s=self.std_nuclear,
        ls_s=self.ls_s_sym,
        ref_n=self.ref_neutron,
        ref_s=self.ref_nuclear,
        breakdown=self.breakdown,
        err_y=self.err_y_d,
        include_3bf=False,
        derivs=[0, 1],
        verbose=verbose,
        rho=self.rho,
    )

    print('Setup time:', time.time() - setup_time_start)

    return self.obs_neutron, self.obs_nuclear, self.obs_sym_energy

energy_dens(add_rest_mass=False, orders='all')

Computes the energy density and standard deviation of the energy density of the chiral EOS. Adds rest mass if desired.

Parameters:

Name	Type	Description	Default
add_rest_mass	bool	The option to include the rest mass in the E/A (and in the energy density).	False
orders	str	Option to either calculate all orders ('all') or only N3LO ('N3LO'). Default is 'all'.	'all'

Returns:

Type	Description
	self.energy_density, self.energy_density_s_stds (numpy.ndarray): The energy density and standard deviation of the energy density.

Source code in src/neutron_rich_bmm/chiral.py

def energy_dens(self, add_rest_mass=False, orders='all'):

    '''
    Computes the energy density and standard deviation of the energy
    density of the chiral EOS. Adds rest mass if desired. 

    Parameters:
        add_rest_mass (bool): The option to include the rest mass in 
            the E/A (and in the energy density). 

        orders (str): Option to either calculate all orders ('all') or 
            only N3LO ('N3LO'). Default is 'all'.

    Returns:
        self.energy_density, self.energy_density_s_stds (numpy.ndarray):
            The energy density and standard deviation of the energy 
            density.
    '''

    # if rest_mass is true, add it and compute
    if add_rest_mass is True:

        if orders == 'all':
            # call energy per particle function
            self.energies_s_mn, self.energy_s_stds = \
                self.energy_per_particle(add_rest_mass=True, orders='all')

            self.energy_density_s = np.zeros([len(self.energies_s_mn), 4])
            self.energy_density_s_stds = np.zeros([len(self.energy_s_stds), 4])

            for i in range(4):
                self.energy_density_s[:,i] = self.energies_s_mn[:,i] \
                    * self.density_all 
                self.energy_density_s_stds[:,i] = \
                    self.energy_s_stds[:,i] * self.density_all 

        elif orders == 'N3LO':

            # call energy per particle function
            self.energies_s_mn, self.energy_s_stds = \
                self.energy_per_particle(add_rest_mass=True, orders='N3LO')

            self.energy_density_s = self.energies_s_mn[:,0] * self.dens_interp
            self.energy_density_s_stds = self.energy_s_stds[:,0] * self.dens_interp

    # otherwise calculate eps and std dev as normal
    else:

        if orders == 'all':

            self.energies_s, self.energy_s_stds = \
                self.energy_per_particle(add_rest_mass=False, orders='all')

            self.energy_density_s = np.zeros([len(self.energies_s), 4])
            self.energy_density_s_stds = np.zeros([len(self.energy_s_stds), 4])

            for i in range(4):
                self.energy_density_s[:,i] = self.energies_s[:,i] \
                    * self.density_all 
                self.energy_density_s_stds[:,i] = self.energy_s_stds[:,i] \
                    * self.density_all 

        elif orders == 'N3LO':

            # call energy per particle function
            self.energies_s, self.energy_s_stds = \
                self.energy_per_particle(add_rest_mass=False, orders='N3LO')

            self.energy_density_s = self.energies_s[:,0] * self.dens_interp
            self.energy_density_s_stds = self.energy_s_stds[:,0] * self.dens_interp

    return self.energy_density_s, self.energy_density_s_stds

energy_per_particle(add_rest_mass=False, case='SNM', orders='all')

A function that ouputs the energy per particle with truncation errors from gsum.

Parameters:

Name	Type	Description	Default
add_rest_mass	bool	If desired, will add the rest mass to the energy per particle.	False
orders	str	Command to determine whether we use all orders for the EFT or just one of them. Default is 'all', but 'N3LO' will convert to only one EFT at N3LO.	'all'

Returns:

Type	Description
	self.energies_s, self.energy_s_stds (numpy.ndarray): The energy per particle (E/A) and standard deviation of E/A.
	self.energies_s_mn, self.energy_s_stds (numpy.ndarray): E/A (inclusive of
	the rest mass) and the standard deviation of E/A. This is returned if
	add_rest_mass is True.

Source code in src/neutron_rich_bmm/chiral.py

def energy_per_particle(self, add_rest_mass=False, case='SNM', orders='all'):

    '''
    A function that ouputs the energy per particle with truncation 
    errors from gsum. 

    Parameters:
        add_rest_mass (bool): If desired, will add the rest mass to the energy per
            particle. 

        orders (str): Command to determine whether we use all orders for the EFT or
            just one of them. Default is 'all', but 'N3LO' will convert
            to only one EFT at N3LO. 

    Returns:
        self.energies_s, self.energy_s_stds (numpy.ndarray): The energy per 
            particle (E/A) and standard deviation of E/A.

        self.energies_s_mn, self.energy_s_stds (numpy.ndarray): E/A (inclusive of 
        the rest mass) and the standard deviation of E/A. This is returned if
        add_rest_mass is True.
    '''

    if case == 'SNM':

        # adding the rest mass in (average p,n)
        self.rest_mass = 938.91875434       # MeV

        # energy per particle
        if orders == 'all':
            self.energies_s = np.array([self.obs_nuclear.get_pred(order=n, deriv=0) \
                                for n in self.orders]).T
            self.energy_s_stds = np.array([self.obs_nuclear.get_std(order=n, deriv=0, \
                                    include_trunc=True) for n in self.orders]).T
        elif orders == 'N3LO':
            self.energies_s = np.array([self.obs_nuclear.get_pred(order=4,\
                                         deriv=0)]).T
            self.energy_s_stds = np.array([self.obs_nuclear.get_std(order=4, \
                                         deriv=0, include_trunc=True)]).T

        # if rest_mass is True, add it to E/A
        if add_rest_mass is True:
            self.energies_s_mn = self.energies_s + self.rest_mass
            return self.energies_s_mn, self.energy_s_stds

        else:
            return self.energies_s, self.energy_s_stds

    elif case == 'PNM':

        # adding the rest mass in (n only here)
        self.rest_mass = 939.6       # MeV

        # energy per particle
        if orders == 'all':
            self.energies_n = np.array([self.obs_neutron.get_pred(order=n, deriv=0) \
                                for n in self.orders]).T
            self.energy_n_stds = np.array([self.obs_neutron.get_std(order=n, deriv=0, \
                                    include_trunc=True) for n in self.orders]).T
        elif orders == 'N3LO':
            self.energies_n = np.array([self.obs_neutron.get_pred(order=4,\
                                         deriv=0)]).T
            self.energy_n_stds = np.array([self.obs_neutron.get_std(order=4, \
                                         deriv=0, include_trunc=True)]).T

        # if rest_mass is True, add it to E/A
        if add_rest_mass is True:
            self.energies_n_mn = self.energies_n + self.rest_mass
            return self.energies_n_mn, self.energy_n_stds

        else:
            return self.energies_n, self.energy_n_stds

f_n(mu, mu_func, guess=0.33)

fsolve function to invert mu(n) to n(mu).

Parameters:

Name	Type	Description	Default
mu	float	The current value in the chemical potential to be inverted.	required
mu_func	function	The function from interp1d that is to be solved for n.	required
guess	float	The guess for fsolve to find a root. Default is 0.33.	0.33

Returns:

Type	Description
	The value of n with respect to the given values of the
	chemical potential.

Source code in src/neutron_rich_bmm/chiral.py

def f_n(self, mu, mu_func, guess=0.33):

    '''
    fsolve function to invert mu(n) to n(mu).

    Parameters:
        mu (float): The current value in the chemical potential 
            to be inverted.

        mu_func (function): The function from interp1d that is 
            to be solved for n.

        guess (float): The guess for fsolve to find a root. 
            Default is 0.33. 

    Returns:
        The value of n with respect to the given values of the 
        chemical potential. 
    '''

    return optimize.fsolve(lambda n : mu - mu_func(n), x0 = guess)[0]

inversion(guess=0.33)

Inverts the function mu(n) in favour of n(mu). For Lambda = 500 MeV, able to use the fsolve function to do this. For Lambda = 450 MeV, must manually invert. ***Note: because of this, the range of mu is not the same for both cases of Lambda.

Parameters:

Name	Type	Description	Default
guess	float	The guess for the fsolve function to invert for Lambda = 500 MeV. If using Lambda = 450 MeV, this argument is ignored. Default is 0.33.	0.33

Returns:

Type	Description
	self.density_mu_N3LO (numpy.ndarray): The array in density that we obtain from inversion.
	self.mu_array_N3LO (numpy.ndarray): The inverted chemical potential array.
	self.kf_N3LO (numpy.ndarray): The new array in kf for the observable container once mu has been inverted.

Source code in src/neutron_rich_bmm/chiral.py

def inversion(self, guess=0.33):

    '''
    Inverts the function mu(n) in favour of n(mu). For Lambda = 500 MeV,
    able to use the fsolve function to do this. For Lambda = 450 MeV,
    must manually invert. ***Note: because of this, the range of mu is 
    not the same for both cases of Lambda. 

    Parameters:
        guess (float): The guess for the fsolve function to invert for 
            Lambda = 500 MeV. If using Lambda = 450 MeV, this argument 
            is ignored. Default is 0.33. 

    Returns:
        self.density_mu_N3LO (numpy.ndarray): The array in density that 
            we obtain from inversion. 

        self.mu_array_N3LO (numpy.ndarray): The inverted chemical potential 
            array. 

        self.kf_N3LO (numpy.ndarray): The new array in kf for the observable 
            container once mu has been inverted. 
    '''

    # set interpolation space outside of cutoff dependence
    nnew = np.linspace(0.05, 0.339, 500)   # should this be changed to 16?
    #nnew = np.linspace(min(self.density_all), max(self.density_all), 500)

    if self.Lambda == 500:

        # N3LO results
        mu_inversion_N3LO = np.asarray(self.invert(self.density_all, self.mu_s[:,3], guess=guess, nnew=nnew))

        # send in the densities in terms of the mu we want
        self.mu_N3LO_array = mu_inversion_N3LO[0] # mu values we now are working with in ALL cases below
        self.density_mu_N3LO = mu_inversion_N3LO[1] # densities found given a specific mu array

        # figure out the density cutoff in mu by finding nearest value
        n_cutoff = mu_inversion_N3LO[1][0]
        nearest = self.density_mu_N3LO.flat[np.abs(self.density_mu_N3LO - n_cutoff).argmin()]
        index = np.where(self.density_mu_N3LO==nearest)[0][0]
        mu_cutoff = self.mu_N3LO_array[index]

        # now use this cutoff to stop from calculating beyond this point 
        index = np.where(self.mu_N3LO_array > mu_cutoff)
        self.mu_N3LO_array = self.mu_N3LO_array[index[0][0]:]
        self.density_mu_N3LO = self.density_mu_N3LO[index[0][0]:]

        # print the range in density to check
        print(min(self.density_mu_N3LO), max(self.density_mu_N3LO))

        # also recalculate kf for the observable container
        self.kf_N3LO = (3.0 * np.square(np.pi) * self.density_mu_N3LO / 2.0)**(1.0/3.0)

        return self.density_mu_N3LO, self.mu_N3LO_array, self.kf_N3LO

    elif self.Lambda == 450:

        # try just inverting the interpolation
        mu_N3LO = self.mu_s[:,3]

        # interpolate the array in mu
        mu_N3LO_interp = interpolate.interp1d(self.density_all, mu_N3LO, kind='cubic')

        # set a new density linspace and mu for tight interpolation
        mu_new = mu_N3LO_interp(nnew)

        # determine region where two solutions exist to avoid it
        mu_min = min(mu_new)
        mu_max = mu_new[0]
        index_mu = np.where(mu_new==mu_min)[0][0]

        # write new density array without the bottom half of mu
        mu_upper = mu_new[index_mu:]
        n_upper = nnew[index_mu:]

        # use interpolation to get more points in density (function)
        n_interp = interpolate.interp1d(mu_upper, n_upper, kind='cubic')

        # use the function to generate density curve in mu
        nearest = mu_N3LO.flat[np.abs(mu_N3LO - mu_max).argmin()]
        index = np.where(mu_N3LO==nearest)[0][0]  

        # chemical potential above density cutoff
        mu_slice = mu_N3LO[index+1:]

        # cut this array to fit the interpolation length
        nearest = mu_slice.flat[np.abs(mu_slice - mu_upper[-1]).argmin()]
        index = np.where(mu_slice==nearest)[0][0]  
        mu_slice = mu_slice[:index-1]

        # new density array from interpolant and mu values
        nmu_new = n_interp(mu_slice)

        # collect arrays from before
        self.mu_N3LO_array = mu_slice
        self.density_mu_N3LO = nmu_new

        # print the density limits now
        print(min(self.density_mu_N3LO), max(self.density_mu_N3LO))

        # also recalculate kf for the observable container
        self.kf_N3LO = (3.0 * np.square(np.pi) * self.density_mu_N3LO / 2.0)**(1.0/3.0)

        return self.density_mu_N3LO, self.mu_N3LO_array, self.kf_N3LO

invert(n, mu, guess=0.33, nnew=None)

Inversion function that uses fsolve. At present, only used for Lambda = 500 MeV.

Parameters:

Name	Type	Description	Default
n	ndarray	The array of densities that we currently possess.	required
mu	ndarray	The array of chemical potential values for the range of densities in n.	required
guess	float	The guess for fsolve. Default is 0.33.	0.33
nnew	ndarray	The values for the new array in density. Default is None.	None

Returns:

Name	Type	Description
mu_new	ndarray	Array of mu results at new points in density.
f_n_result	ndarray	The results of the root finding wrt the new chemical potential array.

Source code in src/neutron_rich_bmm/chiral.py

def invert(self, n, mu, guess=0.33, nnew=None):

    '''
    Inversion function that uses fsolve. At present, only used
    for Lambda = 500 MeV. 

    Parameters:
        n (numpy.ndarray): The array of densities that we currently
            possess.

        mu (numpy.ndarray): The array of chemical potential values 
            for the range of densities in n.

        guess (float): The guess for fsolve. Default is 0.33.

        nnew (numpy.ndarray): The values for the new array in density. 
            Default is None.

    Returns:
        mu_new (numpy.ndarray): Array of mu results at new points 
            in density. 

        f_n_result (numpy.ndarray): The results of the root finding wrt 
            the new chemical potential array.
    '''

    # interpolate the arrays
    mu_interp = interpolate.interp1d(n, mu, kind='cubic')

    # set a new density linspace and mu for tight interpolation
    if nnew is not None:
        mu_new = mu_interp(nnew)
    else:
        mu_new = mu_interp(n)

    # call root finder
    f_n_result = []
    for i in mu_new:
        f_n_result.append(self.f_n(i, mu_interp, guess=guess))

    return mu_new, f_n_result

pressure(orders='all', matter='SNM')

The pressure of the chiral EOS. Note: Rest mass will not affect this calculation, so no option for adding the rest mass is included here.

Parameters:

Name	Type	Description	Default
orders	str	Command to determine how many orders in the EFT expansion to calculate. Default is 'all', but can be set to 'N3LO' to only return that one.	'all'

Returns:

Type	Description
	self.pressure_s, self.pressure_s_stds (numpy.ndarray): The pressure and standard deviation of the pressure.

Source code in src/neutron_rich_bmm/chiral.py

def pressure(self, orders='all', matter='SNM'):

    '''
    The pressure of the chiral EOS. 
    Note: Rest mass will not affect this calculation, so no option for 
    adding the rest mass is included here. 

    Parameters:
        orders (str): Command to determine how many orders in the EFT 
            expansion to calculate. Default is 'all', but can be set to 
            'N3LO' to only return that one.

    Returns:
        self.pressure_s, self.pressure_s_stds (numpy.ndarray): The 
            pressure and standard deviation of the pressure. 
    '''

    if matter == 'SNM':
        # set up the pressure and std dev
        pressures_s = []
        pressure_s_stds = []
        pressure_s_cov_n = {}
        pressure_s_cov_arr = np.zeros([len(self.obs_nuclear.kf_interp), len(self.obs_nuclear.kf_interp), len(self.orders)])

        if orders == 'all':

            for i, n in enumerate(self.orders):
                pressure_s = compute_pressure(
                    self.obs_nuclear.density_interp,
                    self.obs_nuclear.kf_interp,
                    dE=self.obs_nuclear.get_pred(order=n, deriv=1)
                )
                pressure_s_cov = compute_pressure_cov(
                    self.obs_nuclear.density_interp,
                    self.obs_nuclear.kf_interp,
                    dE_cov=self.obs_nuclear.get_cov(order=n, deriv1=1, deriv2=1)
                )

                pressure_s_cov_n[n] = pressure_s_cov

                pressure_s_std = np.sqrt(np.diag(pressure_s_cov))

                pressures_s.append(pressure_s)
                pressure_s_stds.append(pressure_s_std)

            pressures_s = np.array(pressures_s).T
            pressure_s_stds = np.array(pressure_s_stds).T

            for i,value in zip(range(len(self.orders)), pressure_s_cov_n.values()):
                pressure_s_cov_arr[:,:,i] = np.array(value).T

            # transform into class variables
            self.pressures_s = pressures_s 
            self.pressure_s_stds = pressure_s_stds 
            self.pressure_s_cov = pressure_s_cov_arr

            return self.pressures_s, self.pressure_s_stds, self.pressure_s_cov 

        elif orders == 'N3LO':

            pressures_s = compute_pressure(
                self.obs_nuclear.density_interp,
                self.obs_nuclear.kf_interp,
                dE=self.obs_nuclear.get_pred(order=4, deriv=1)
            )

            pressure_s_cov = compute_pressure_cov(
                self.obs_nuclear.density_interp,
                self.obs_nuclear.kf_interp,
                dE_cov=self.obs_nuclear.get_cov(order=4, deriv1=1, deriv2=1)
            )

            pressure_s_stds = np.sqrt(np.diag(pressure_s_cov))

            pressures_s = np.array(pressures_s).T
            pressure_s_stds = np.array(pressure_s_stds).T
            pressure_s_cov = np.array(pressure_s_cov).T

            # transform into class variables
            self.pressures_s = pressures_s 
            self.pressure_s_stds = pressure_s_stds 
            self.pressure_s_cov = pressure_s_cov

            return self.pressures_s, self.pressure_s_stds, self.pressure_s_cov 

    if matter == 'PNM':
        # set up the pressure and std dev
        pressures_n = []
        pressure_n_stds = []
        pressure_n_cov_n = {}
        pressure_n_cov_arr = np.zeros([len(self.obs_neutron.kf_interp), len(self.obs_neutron.kf_interp), len(self.orders)])

        if orders == 'all':

            for i, n in enumerate(self.orders):
                pressure_n = compute_pressure(
                    self.obs_neutron.density_interp,
                    self.obs_neutron.kf_interp,
                    dE=self.obs_neutron.get_pred(order=n, deriv=1)
                )
                pressure_n_cov = compute_pressure_cov(
                    self.obs_neutron.density_interp,
                    self.obs_neutron.kf_interp,
                    dE_cov=self.obs_neutron.get_cov(order=n, deriv1=1, deriv2=1)
                )

                pressure_n_cov_n[n] = pressure_n_cov

                pressure_n_std = np.sqrt(np.diag(pressure_n_cov))

                pressures_n.append(pressure_n)
                pressure_n_stds.append(pressure_n_std)

            pressures_n = np.array(pressures_n).T
            pressure_n_stds = np.array(pressure_n_stds).T

            for i,value in zip(range(len(self.orders)), pressure_n_cov_n.values()):
                pressure_n_cov_arr[:,:,i] = np.array(value).T

            # transform into class variables
            self.pressures_n = pressures_n 
            self.pressure_n_stds = pressure_n_stds 
            self.pressure_n_cov = pressure_n_cov_arr

            return self.pressures_n, self.pressure_n_stds, self.pressure_n_cov 

Perturbative QCD EOS code

Gorda(mu, X, Nf, mu_FG=None)

function since it is already calibrated.

Parameters:

Name	Type	Description	Default
mu	linspace	The quark chemical potential needed to generate the mean and standard deviation of the pressure in the pQCD EOS model.	required
X	int	The value of the renormalization scale parameter.	required
Nf	int	The number of flavours of quarks considered.	required
mu_FG	ndarray	The FG chemical potential array.	None

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/gorda_model.py

def __init__(self, mu, X, Nf, mu_FG=None):

    '''
    Sets up the chemical potential for the evaluate() 
    function. pQCD EOS model only possesses the evaluate()
    function since it is already calibrated. 

    Parameters:
        mu (numpy.linspace): The quark chemical potential needed 
            to generate the mean and standard deviation of the pressure
            in the pQCD EOS model. 

        X (int): The value of the renormalization scale parameter.

        Nf (int): The number of flavours of quarks considered.

        mu_FG (numpy.ndarray): The FG chemical potential array.

    Returns:
        None. 
    '''

    # set the chemical potential (change later if needed)
    self.mu = mu
    self.X = X
    self.Nf = Nf
    if mu_FG is not None:
        self.mu_FG = mu_FG
        self.mU_FG = mu_FG[:,None]

    # instantiate the Gorda class
    self.gorda = PQCD(X=self.X, Nf=self.Nf)

    # instantiate number of orders
    self.orders = 3

    # set up coefficients and Truncation class
    coeffs = np.array([self.gorda.c_0(self.mu), self.gorda.c_1(self.mu), self.gorda.c_2(self.mu)]).T
    self.trunc = Truncation(x=self.mu, x_FG=self.mu_FG, norders=3, \
                            yref=self.gorda.yref, expQ=self.gorda.expQ, coeffs=coeffs)

    return None

evaluate(input_space=None, N2LO=True, scaled=True)

The evaluation function for Taweret to obtain the calibrated mean and variance of the pQCD model.

Parameters:

Name	Type	Description	Default
input_space	ndarray	The number density array.	None
N2LO	bool	If we only want the data from the pQCD pressure at N2LO. Default is True.	True
scaled	bool	If the data is scaled, this is True. Else, it is False. Default is True.	True

Returns:

Type	Description
	mean, std_dev (numpy.ndarray): The mean and standard deviations at the selected points in the input space of the pQCD pressure.

Source code in src/neutron_rich_bmm/gorda_model.py

def evaluate(self, input_space=None, N2LO=True, scaled=True):

    '''
    The evaluation function for Taweret to obtain the 
    calibrated mean and variance of the pQCD model.

    Parameters:
        input_space (numpy.ndarray): The number density array.

        N2LO (bool): If we only want the data from the pQCD pressure 
            at N2LO. Default is True.

        scaled (bool): If the data is scaled, this is True. 
            Else, it is False. Default is True.

    Returns:
        mean, std_dev (numpy.ndarray): The mean and standard 
            deviations at the selected points in the input space 
            of the pQCD pressure.
    '''

    # KLW formalism for the pQCD EOS pressure
    conversion = (1000)**4.0/(197.327)**3.0

    # correct input space
    if input_space is not None:
        input_space = None

    # call interpolation and work through
    _, _, _ = self.trunc.gp_interpolation(center=0.0, sd=1.0)

    # coeffs and data solved at mu_FG
    coeffs_FG = np.array([self.gorda.c_0(self.mu_FG), self.gorda.c_1(self.mu_FG), self.gorda.c_2(self.mu_FG)]).T
    data_FG = gm.partials(coeffs_FG, ratio=self.gorda.expQ(self.mU_FG), \
                          ref=self.gorda.yref(self.mU_FG), orders=[range(3)])
    _, coeffs_trunc, std_dev = \
        self.trunc.uncertainties(data=data_FG, expQ=self.gorda.expQ, yref=self.gorda.yref)

    # fix lower densities when run below n = 0.16 fm^-3
    for j in range(self.orders):
        for i in range(len(std_dev)):
            if np.isnan(std_dev[i,j]) == True or np.isinf(std_dev[i,j]) == True:
                std_dev[i,j] = 1e10

    # KLW pressure call
    pressure_n = self.gorda.pressure_KLW(self.mu_FG)

    # put these into an array
    mean = np.array([pressure_n["LO"], pressure_n["NLO"], pressure_n["N2LO"]]).T

    if N2LO is True:
        mean = mean[:,2]
        std_dev = std_dev[:,2]

    if scaled is True:
        mean = mean/self.gorda.yref(self.mU_FG)
        std_dev = std_dev/self.gorda.yref(self.mU_FG)
        return mean, std_dev 
    else:
        return mean, std_dev

log_likelihood_elementwise()

The log likelihood function that would calculate this quantity in Taweret.

Source code in src/neutron_rich_bmm/gorda_model.py

def log_likelihood_elementwise(self):
    '''
    The log likelihood function that would calculate
    this quantity in Taweret.
    '''
    return None

set_prior()

The prior function to set a prior for Taweret to take in for this model.

Source code in src/neutron_rich_bmm/gorda_model.py

def set_prior(self):
    '''
    The prior function to set a prior for 
    Taweret to take in for this model.
    '''
    return None

PQCD(X=1, Nf=2)

the rest of the code.

:Example: PQCD(X=1, Nf=2)

Parameters:

Name	Type	Description	Default
X	int	The value of the coefficient multiplying mu to determine the renormalization scale. Default is 1, as in Gorda et al. (2023).	1
Nf	int	The number of different quark flavours being considered. Default is 2, for SNM. NSM is 3.	2

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def __init__(self, X=1, Nf=2):

    '''
    Initialize the class so that mu and lambda_bar
    have been set and constants have been defined for
    the rest of the code.

    :Example:
        PQCD(X=1, Nf=2)

    Parameters:        
        X (int): The value of the coefficient multiplying mu
            to determine the renormalization scale. Default is 1,
            as in Gorda et al. (2023).

        Nf (int): The number of different quark flavours being 
            considered. Default is 2, for SNM. NSM is 3.

    Returns:
        None.
    '''

    # initialize constants
    self.X = X
    self.Nc = 3
    self.Nf = Nf
    self.Ca = self.Nc
    self.Cf = (self.Nc**2.0 - 1.0)/(2.0*self.Nc)
    self.lambda_MS = 0.38 #[GeV]

    self.beta0 = (11.0*self.Ca - 2.0*self.Nf)/3.0   # for alpha_s

    # QCD constants 
    self.b0 = (11.0/6.0) * self.Nc * self.Nf**(-1.0) - (1.0/3.0)
    self.b1 = (17.0/12.0) * self.Nc**2 * self.Nf**(-2) + \
    (self.Nc**(-1.0) * self.Nf**(-1.0))/8.0 - \
    (13.0/24.0) * self.Nc * self.Nf**(-1.0)
    self.delta = -0.85638

    # pressure constants
    self.dA = 8.0
    self.a11 = -0.75 * self.dA * self.Nc**(-1.0)
    self.a21 = -(3.0/8.0) * self.dA * self.Nc**(-1.0)  
    self.a22 = self.b0 * self.a11 
    self.a23 = self.dA * (((11.0/6.0) - (np.pi**2)/8.0 - (5.0 * np.log(2))/4.0 + \
                np.log(2)**2  - (3.0/16.0)*self.delta) * self.Nc**(-1.0) - \
                (415.0/192.0)*self.Nf**(-1.0) - \
                (51.0/64.0)*self.Nc**(-2) * self.Nf**(-1.0))

    return None

alpha_s(mu, loop=2)

The function for alpha_s with an option for either first order or second order (loops) inclusion.

:Example: PQCD.alpha_s(mu, loop=1)

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential.	required
loop	int	The order at which alpha_s is calculated. Default is 2.	2

Returns:

Name	Type	Description
alpha_s	ndarray	The values of alpha_s at the order requested.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def alpha_s(self, mu, loop=2):

    '''
    The function for alpha_s with an option for either 
    first order or second order (loops) inclusion.

    :Example:
        PQCD.alpha_s(mu, loop=1)

    Parameters:
        mu (numpy.ndarray): The quark chemical potential. 

        loop (int): The order at which alpha_s is calculated.
            Default is 2.

    Returns:
        alpha_s (numpy.ndarray): The values of alpha_s at 
            the order requested.
    '''

    lambda_bar = 2. * self.X * mu   # X range [1/2, 2]

    # renormalization scale
    ell = np.log((lambda_bar**2.0)/self.lambda_MS**2.0)

    if (loop == 1):
        beta1 = 0.0
        # running coupling constant
        alpha_s = ((4.0*np.pi)/(self.beta0*ell))
    elif (loop == 2):
        beta1 = (17.0/3.0)*self.Ca**2.0 - self.Cf*self.Nf - (5.0/3.0)*(self.Ca*self.Nf) 
        # running coupling constant
        alpha_s = ((4.0*np.pi)/(self.beta0*ell)) * (1.0 - (2.0*beta1*np.log(ell))/(ell*self.beta0**2.0))
    else:
        raise ValueError('Loop selected must be 1 or 2.')

    return alpha_s

c_0(x)

The c0 coefficient of pQCD.

Parameters:

Name	Type	Description	Default
x	array	The quark chemical potential.	required

Returns:

Type	Description
	The value of the coefficient at each point
	in the array.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def c_0(self, x):

    '''
    The c0 coefficient of pQCD. 

    Parameters:
        x (numpy.array): The quark chemical potential.

    Returns:
        The value of the coefficient at each point
        in the array.
    '''
    return np.ones(len(x))

c_1(x)

The c1 coefficient of pQCD.

Parameters:

Name	Type	Description	Default
x	array	The quark chemical potential.	required

Returns:

Type	Description
	The value of the coefficient at each point
	in the array.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def c_1(self, x):

    '''
    The c1 coefficient of pQCD. 

    Parameters:
        x (numpy.array): The quark chemical potential.

    Returns:
        The value of the coefficient at each point
        in the array.
    '''
    return (self.a11*np.ones(len(x))/self.Nf)

c_2(x)

The c2 coefficient of pQCD.

Parameters:

Name	Type	Description	Default
x	array	The quark chemical potential.	required

Returns:

Type	Description
	The value of the coefficient at each point
	in the array.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def c_2(self, x):

    '''
    The c2 coefficient of pQCD. 

    Parameters:
        x (numpy.array): The quark chemical potential.

    Returns:
        The value of the coefficient at each point
        in the array.
    '''
    lambda_bar = 2.0 * self.X * x
    one = self.a21 * np.log(self.Nf * self.alpha_s(x, loop=2)/np.pi)
    two = self.a22 * np.log(lambda_bar/(2.0*x))
    three = self.a23
    return (one + two + three)/self.Nf

expQ(x)

The expansion parameter function for pQCD. Default here is Nf * alpha_s / pi, but another may be substituted.

Parameters:

Name	Type	Description	Default
x	array	The quark chemical potential.	required

Returns:

Type	Description
	The value of the expansion parameter at each
	point in the quark chemical potential array.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def expQ(self, x): 

    '''
    The expansion parameter function for pQCD. 
    Default here is Nf * alpha_s / pi, but another may
    be substituted.

    Parameters:
        x (numpy.array): The quark chemical potential.

    Returns:
        The value of the expansion parameter at each
        point in the quark chemical potential array.
    '''
    return (self.Nf * self.alpha_s(x, loop=2)/np.pi)[:,0]

inversion(n_mu=None)

Function to invert n(mu) to obtain mu(n).

Parameters:

Name	Type	Description	Default
n_mu	array	Linspace over n_q for the inversion.	None

Returns:

Name	Type	Description
n_mu	ndarray	The number density.
f_mu_2_new	array	The array corresponding to the inverted function values for mu(n).
f_mu_FG_new	ndarray	The values of the FG chemical potential.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def inversion(self, n_mu=None):

    '''
    Function to invert n(mu) to obtain mu(n). 

    Parameters:
        n_mu (numpy.array): Linspace over n_q for the inversion.

    Returns:
        n_mu (numpy.ndarray): The number density. 

        f_mu_2_new (numpy.array): The array corresponding to the 
            inverted function values for mu(n).

        f_mu_FG_new (numpy.ndarray): The values of the FG chemical potential. 
    '''

    # write the root finding function
    def f_mu_2(n, guess):  
        return opt.fsolve(lambda mu : n - self.n_mu(mu), x0 = guess)[0]

    # FG chemical potential (mu_FG)
    def f_mu_FG(n, guess):
        return opt.fsolve(lambda mu_FG : n - self.n_FG_mu(mu_FG), x0 = guess, xtol=1e-10)[0]

    if n_mu is None:
        n_mu = np.linspace(0.01, 0.8, 1000)  # n_q right now -> 3*n_B = (0.0, 1.25) GeV^3

    # invert to get mu(n^2)
    f_mu_2_new = []
    f_mu_FG_new = []
    for i in n_mu:
        f_mu_2_new.append(f_mu_2(i, guess=2.0))
        f_mu_FG_new.append(f_mu_FG(i, guess=2.0))

    # convert to n_B here to not make mistake later
    self.nB = n_mu/3.0   # n_q/3

    # convert the mu arrays over (still muB/3 -> mu_q)
    f_mu_2_new = np.asarray(f_mu_2_new)
    self.f_mu_FG_new = np.asarray(f_mu_FG_new)

    return n_mu, f_mu_2_new, self.f_mu_FG_new

mask_array(array, neg=False, fill_value=None) staticmethod

Returns a masked array from an original array that contains unwanted nan values. Can also fill values in the place of the mask, if desired.

:Example: PQCD.mask_array(array=mu, fill_value=0)

Parameters:

Name	Type	Description	Default
array	ndarray	The array with values to mask.	required
neg	bool	If False, will not also mask negative values. If True, will check for negative values and mask using current fill value.	False
fill_value	(int, float)	The value with which to fill the mask. Default is None.	None

Returns:

Name	Type	Description
masked_array	ndarray	The masked array with or without filled values.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

@staticmethod
def mask_array(array, neg=False, fill_value=None):

    '''
    Returns a masked array from an original array that contains 
    unwanted nan values. Can also fill values in the place of 
    the mask, if desired.

    :Example:
        PQCD.mask_array(array=mu, fill_value=0)

    Parameters:
        array (numpy.ndarray): The array with values to mask. 

        neg (bool): If False, will not also mask negative values. If True,
            will check for negative values and mask using current
            fill value. 

        fill_value (int, float): The value with which to fill the mask. 
            Default is None.

    Returns:
        masked_array (numpy.ndarray): The masked array with or without 
            filled values.
    '''

    # initialize the mask and new array
    mask = np.zeros(len(array), dtype=bool)
    masked_array = np.zeros(len(array))

    # check for neg bool and mask if needed
    if neg is True:
        for i in range(len(array)):
            if array[i] < 0.0:
                array[i] = fill_value

    # check for nan values and replace in mask
    for i in range(len(array)):
        if np.isnan(array[i]):
            mask[i] = True

    # mask the original array
    masked_array = np.ma.array(array, mask=mask)

    # if there were given fill values, fill the array
    if fill_value is not None:
        masked_array = np.ma.filled(masked_array, fill_value=fill_value)

    return masked_array

mu_1(mu_FG)

The mu1 term of pQCD.

Parameters:

Name	Type	Description	Default
mu_FG	array	The FG quark chemical potential.	required

Returns:

Type	Description
	The value of mu1 at each point
	in the array.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def mu_1(self, mu_FG): 

    '''
    The mu1 term of pQCD. 

    Parameters:
        mu_FG (numpy.array): The FG quark chemical potential.

    Returns:
        The value of mu1 at each point
        in the array.
    '''

    mU_FG = mu_FG[:, None]
    numerator = self.c_1(mu_FG) * self.expQ(mU_FG) * (self.Nf * mu_FG**3.0/(np.pi**2.0))
    #numerator = self.c_1(mu_FG) * self.alpha_s(mu_FG, loop=2) * self.n_FG_mu(mu_FG) #checked 
    denominator = self.c_0(mu_FG) * (self.Nf * 3.0 * mu_FG**2.0 / (np.pi**2.0)) #checked
    return -numerator/denominator

mu_2(mu_FG)

The mu2 term of pQCD.

Parameters:

Name	Type	Description	Default
mu_FG	array	The FG quark chemical potential.	required

Returns:

Type	Description
	The value of mu2 at each point
	in the array.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def mu_2(self, mu_FG): 

    '''
    The mu2 term of pQCD. 

    Parameters:
        mu_FG (numpy.array): The FG quark chemical potential.

    Returns:
        The value of mu2 at each point
        in the array.
    '''

    mU_FG = mu_FG[:, None]
    self.derivalpha = - self.beta0 * self.alpha_s(mu_FG, loop=2)**2.0 / (2.0 * mu_FG * np.pi)
    self.derivQ = self.derivalpha * (self.expQ(mU_FG)/self.alpha_s(mu_FG, loop=2))
    secderivP0 = self.Nf * 3.0 * mu_FG**2.0 / (np.pi**2.0)

    first = 0.5 * self.c_0(mu_FG) * self.mu_1(mu_FG)**2.0 * (self.Nf * 6.0 * mu_FG / (np.pi**2.0)) #checked
    #third = self.c_1(mu_FG) * derivalpha * self.pressure_FG(mu_FG) #checked
    third = self.c_1(mu_FG) * self.derivQ * self.pressure_FG(mu_FG)

  #  second = self.mu_1(mu_FG) * self.c_1(mu_FG) * self.alpha_s(mu_FG, loop=2) * secderivP0 #checked
    second = self.mu_1(mu_FG) * self.c_1(mu_FG) * self.expQ(mU_FG) * secderivP0
  #  fourth = self.c_2(mu_FG) * self.alpha_s(mu_FG, loop=2)**2.0 * (self.Nf * mu_FG**3.0 / (np.pi**2.0)) #checked
    fourth = self.c_2(mu_FG) * self.expQ(mU_FG)**2.0 * (self.n_FG_mu(mu_FG))

    mu2 = -(first + second + third + fourth)/(self.c_0(mu_FG) * secderivP0)

    return mu2

n_FG_mu(mu)

Free quark number density calculation from P_FG(mu) derivative. Note that this will yield a different n value for the same mu value as n_mu yields, so scaling must be taken into account.

Parameters:

Name	Type	Description	Default
mu	ndarray	The input chemical potential.	required

Returns:

Name	Type	Description
n_FG	ndarray	The result of the number density at the input chemical potential.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def n_FG_mu(self, mu):

    '''
    Free quark number density calculation from P_FG(mu)
    derivative. Note that this will yield a different n 
    value for the same mu value as n_mu yields, so 
    scaling must be taken into account.

    Parameters:
        mu (numpy.ndarray): The input chemical potential.

    Returns:
        n_FG (numpy.ndarray): The result of the number density at the input 
            chemical potential.
    '''

    n_FG = self.Nf * mu**3.0 / (np.pi**2.0)

    return n_FG

n_convert_mu(density)

The function that converts from a desired number density to quark chemical potential.

Parameters:

Name	Type	Description	Default
density	linspace	The density array.	required

Returns:

Name	Type	Description
mu_FG	ndarray	The FG quark chemical potential.
mu_n	ndarray	The quark chemical potential array.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def n_convert_mu(self, density):

    '''
    The function that converts from a desired number 
    density to quark chemical potential. 

    Parameters:
        density (numpy.linspace): The density array.

    Returns:
        mu_FG (numpy.ndarray): The FG quark chemical potential. 

        mu_n (numpy.ndarray): The quark chemical potential array.
    '''

    n0 = 0.164

    # take the density and invert to obtain the corresponding mu
    n_q = density*3.0  # n_q [fm^-3]

    # convert to GeV^3 for mu_q
    conversion_fm3 = ((1000.0)**(3.0))/((197.33)**(3.0)) # [fm^-3]  (do the opposite of this)
    n_q = n_q/conversion_fm3

    # invert to obtain mu_n and mu_FG (returns quark chemical potentials)
    _, mu_n, mu_FG = self.inversion(n_mu=n_q)

    return mu_FG, mu_n

n_mu(mu)

Simple first derivative of the pressure (up to second order in P(mu) equation) with respect to the chemical potential to obtain the number density with respect to mu.

NOTE: this derivative is to second order in P, so the equation will be: n(mu) = dP0/dmu + dP1/mu + dP2/mu.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential.	required

Returns:

Name	Type	Description
n	ndarray	The array of the number density with respect to the chemical potential.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def n_mu(self, mu):

    '''
    Simple first derivative of the pressure (up to second
    order in P(mu) equation) with respect to the chemical 
    potential to obtain the number density with respect to mu.

    NOTE: this derivative is to second order in P, so the
    equation will be: n(mu) = dP0/dmu + dP1/mu + dP2/mu.

    Parameters:
        mu (numpy.ndarray): The quark chemical potential. 

    Returns:
        n (numpy.ndarray): The array of the number density with 
            respect to the chemical potential.
    '''

    n_mu = ndt.Derivative(self.pressure_mu, step=0.00001, method='central')
    n = n_mu(mu)

    return n

pressure_FG(mu)

The FG contribution to the pressure in terms of mu.

:Example: PQCD.pressure_FG(mu=np.linspace())

Parameters:

Name	Type	Description	Default
mu	linspace	The original range in mu.	required

Returns:

Name	Type	Description
p_FG	linspace	The zeroth order (LO) contribution the pressure.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def pressure_FG(self, mu):

    '''
    The FG contribution to the pressure in terms of mu.

    :Example:
        PQCD.pressure_FG(mu=np.linspace())

    Parameters:
        mu (numpy.linspace): The original range in mu.

    Returns:
        p_FG (numpy.linspace): The zeroth order (LO) contribution
            the pressure. 
    '''

    p_FG = (self.Nf * self.Nc * mu**4.0) / (12.0 * np.pi**2.0)  # extremely general form

    return p_FG

pressure_KLW(mu_FG)

The pressure equation written out for the KLW inversion.

Parameters:

Name	Type	Description	Default
mu_FG	array	The FG quark chemical potential.	required

Returns:

Name	Type	Description
pressure_n	dict	The values of the pressure from the KLW inversion at LO, NLO, and N2LO.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def pressure_KLW(self, mu_FG):

    '''
    The pressure equation written out for the KLW
    inversion. 

    Parameters:
        mu_FG (numpy.array): The FG quark chemical potential.

    Returns:
        pressure_n (dict): The values of the pressure from the KLW
            inversion at LO, NLO, and N2LO.
    '''

    # mu_FG assignment
    mU_FG = mu_FG[:,None]

    pressure_0 = self.c_0(mu_FG) * self.yref(mU_FG)

    pressure_1 = self.c_1(mu_FG)*self.expQ(mU_FG)*self.yref(mU_FG) + \
    self.mu_1(mu_FG) * self.c_0(mu_FG) * self.n_FG_mu(mu_FG)

    pressure_2 = self.c_2(mu_FG)*self.expQ(mU_FG)**2.0*self.yref(mU_FG) + \
    (self.mu_2(mu_FG) * self.c_0(mu_FG) * self.n_FG_mu(mu_FG) + \
     0.5 * self.mu_1(mu_FG)**2.0 * (self.Nf * 3.0 * mu_FG**2.0 / np.pi**2.0) + \
     self.mu_1(mu_FG) * self.c_1(mu_FG) * self.n_FG_mu(mu_FG) * self.expQ(mU_FG))

    # organise so each order contains the last
    pressure_LO = pressure_0
    pressure_NLO = pressure_0 + pressure_1
    pressure_N2LO = pressure_0 + pressure_1 + pressure_2

    # stash in dict
    pressure_n = {
        "LO": pressure_LO,
        "NLO": pressure_NLO,
        "N2LO": pressure_N2LO
    }

    return pressure_n

pressure_mu(mu, order=2)

The pressure with respect to mu for NLO or N2LO. Note: alpha_s has been used up to second order here to maintain the validity of the value of alpha_s at 2 GeV. From Gorda et al. (2023), supplemental material.

:Example: PQCD.pressure_mu(mu=np.linspace(), order=1)

Parameters:

Name	Type	Description	Default
mu	linspace	The original range of mu.	required
order	int	The order at which the pressure is calculated. Default is 2, options are 1 and 2.	2

Returns:

Name	Type	Description
pressure	ndarray	The pressure as a function of mu.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def pressure_mu(self, mu, order=2):

    '''
    The pressure with respect to mu for NLO or N2LO.
    Note: alpha_s has been used up to second order here to
    maintain the validity of the value of alpha_s at 2 GeV.
    From Gorda et al. (2023), supplemental material.

    :Example:
        PQCD.pressure_mu(mu=np.linspace(), order=1)

    Parameters:
        mu (numpy.linspace): The original range of mu.

        order (int): The order at which the pressure is calculated. Default is 2,
            options are 1 and 2.

    Returns:
        pressure (numpy.ndarray): The pressure as a function of mu.
    '''

    lambda_bar = 2. * self.X * mu     # X range [1/2, 2]

    # log terms
    logQ = np.log(self.Nf*self.alpha_s(mu, loop=2) / np.pi)
    logL = np.log(lambda_bar/(2.0*mu))

    if order == 1:
        pressure = self.pressure_FG(mu) * (1.0 + self.a11*(self.alpha_s(mu, loop=2)/np.pi))

    elif order == 2:

        # N2LO terms
        first = self.a21 * logQ
        second = self.a22 * logL
        third = self.a23

        # full expression
        pressure_NLO = self.a11*(self.alpha_s(mu, loop=2)/np.pi)
        pressure_N2LO = ((first + second + third)*self.Nf \
                         *(self.alpha_s(mu,loop=2)/np.pi)**2)

        pressure = self.pressure_FG(mu) * (1.0 + pressure_NLO + pressure_N2LO)

    return pressure

pressure_old(mu, order=2)

Old version of the pressure, using Nf=3 implicitly. Not used in the paper.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential.	required

Returns:

Name	Type	Description
pressure	ndarray	The value of the pressure at the chemical potentials.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def pressure_old(self, mu, order=2):

    '''
    Old version of the pressure, using Nf=3
    implicitly. Not used in the paper.

    Parameters:
        mu (numpy.ndarray): The quark chemical potential.

    Returns:
        pressure (numpy.ndarray): The value of the pressure at the
            chemical potentials. 
    '''

    lambda_bar = 2.0 * self.X * mu

    if order == 1:

        pressure = self.pressure_FG(mu) * (1.0 - 2.0 * \
            self.alpha_s(mu, loop=2)/np.pi)

        return pressure

    elif order == 2:

        first = 2.0 * self.alpha_s(mu, loop=2) / np.pi
        second = 0.303964 * self.alpha_s(mu, loop=2)**2.0 \
            * np.log(self.alpha_s(mu, loop=2))
        second_2 = self.alpha_s(mu, loop=2)**2.0 * (0.874355 \
            + 0.911891 * np.log(lambda_bar/mu))
        pressure = self.pressure_FG(mu) * (1.0 - first - \
            second - second_2)

        return pressure

yref(x)

The function to evaluate yref from pQCD.

Parameters:

Name	Type	Description	Default
x	array	The quark chemical potential.	required

Returns:

Name	Type	Description
yref	numpy.2darray	The [:,None] array for yref.

Source code in src/neutron_rich_bmm/pqcd_reworked.py

def yref(self, x):

    '''
    The function to evaluate yref from pQCD.

    Parameters:
        x (numpy.array): The quark chemical potential.

    Returns:
        yref (numpy.2darray): The [:,None] array for yref.
    ''' 

    yref = ((self.Nf * self.Nc * x**4.0) / (12.0 * np.pi**2.0))  # FG pressure for any flavour, colour
    yref = yref[:,0]
    return yref

Kurkela(X=2, Nf=3)

A simple class to find the roots and perform the integration needed for the Kurkela et al. (2010) CQM paper formulation of the pressure of massless quarks.

Also converts from chemical potential to number density if needed.

Source code in src/neutron_rich_bmm/pqcd.py

def __init__(self, X=2, Nf=3):

    # initialise PQCD class to use below
    self.X = X
    self.Nf = Nf
    self.pqcd_root = PQCD(X=self.X, Nf=self.Nf)

    return None

c_0(mu)

LO coefficent for the pQCD EOS using Kurkela formalism.

Parameters:

Name	Type	Description	Default
mu	ndarray	The chemical potential for the pQCD EOS.	required

Returns:

Type	Description
	np.ones(len(mu)) (numpy.ndarray): The leading order coefficient.

Source code in src/neutron_rich_bmm/pqcd.py

def c_0(self, mu):

    '''
    LO coefficent for the pQCD EOS using Kurkela formalism.

    Parameters:
        mu (numpy.ndarray): The chemical potential for the pQCD EOS.

    Returns:
        np.ones(len(mu)) (numpy.ndarray): The leading order coefficient. 
    '''

    return np.ones(len(mu))

c_1(mu)

The NLO coefficient for the pQCD EOS using Kurkela formalism.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential array.	required

Returns:

Type	Description
	The value of the coefficient.

Source code in src/neutron_rich_bmm/pqcd.py

def c_1(self, mu):

    '''
    The NLO coefficient for the pQCD EOS using Kurkela
    formalism. 

    Parameters:
        mu (numpy.ndarray): The quark chemical potential array.

    Returns:
        The value of the coefficient. 
    '''

    # reshape the mu array for yref
    mu_arr = mu[:, None]

    # calculate the pressure 
    _, p_1, _ = self.pressure(mu, scaled=False)

    # take only needed pressure (fix later)
    p_1 = p_1[1:]

    one = p_1/self.yref(mu_arr)
    two = self.c_0(mu)

    return (one - two)/self.pqcd_root.alpha_s(mu, loop=2)

c_2(mu)

The N2LO coefficent for the pQCD EOS using Kurkela formalism.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential.	required

Returns:

Type	Description
	The value of the coefficient.

Source code in src/neutron_rich_bmm/pqcd.py

def c_2(self, mu):

    '''
    The N2LO coefficent for the pQCD EOS using Kurkela
    formalism. 

    Parameters:
        mu (numpy.ndarray): The quark chemical potential. 

    Returns:
        The value of the coefficient. 
    '''

    # reshape the mu array again
    mu_arr = mu[:, None]

    # calculate the pressure
    _, _, p_2 = self.pressure(mu, scaled=False)

    # take only needed pressure (fix later)
    p_2 = p_2[1:]

    one = p_2/self.yref(mu_arr)
    two = self.c_1(mu) * self.pqcd_root.alpha_s(mu, loop=2)
    three = self.c_0(mu)

    return (one - two - three)/np.square(self.pqcd_root.alpha_s(mu,loop=2))

expQ(mu)

The expansion parameter, Q, of the pQCD EOS using the Kurkela formalism.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potentials.	required

Returns:

Type	Description
	The value of expQ at the chosen potentials.

Source code in src/neutron_rich_bmm/pqcd.py

def expQ(self, mu):

    '''
    The expansion parameter, Q, of the pQCD EOS using
    the Kurkela formalism. 

    Parameters:
        mu (numpy.ndarray): The quark chemical potentials. 

    Returns:
        The value of expQ at the chosen potentials. 
    '''

    return self.pqcd_root.alpha_s(mu, loop=2)[:,0]

find_roots()

Function to determine the value of the chemical potential when the number density reaches zero, so that we can define the limits of our integration for the pressure.

Returns:

Type	Description
	[root1, root2] (numpy.array): The roots of the first and second order number densities in the chemical potential.

Source code in src/neutron_rich_bmm/pqcd.py

def find_roots(self):

    '''
    Function to determine the value of the chemical potential
    when the number density reaches zero, so that we can define
    the limits of our integration for the pressure.

    Parameters:
        None.

    Returns:
        [root1, root2] (numpy.array): The roots of the first and 
            second order number densities in the chemical potential. 
    '''

    # set the linspace needed to solve precisely for the root
    mu_root = np.linspace(0.2, 2.5, 10000)

    # set up the root finder for both orders in n
    root_first = np.vstack((mu_root, self.n_first(mu_root)))
    index1 = np.where(root_first[1,:] > 0.0)[0][0]
    root1 = mu_root[index1]

    root_second = np.vstack((mu_root, self.n_second(mu_root)))
    index2 = np.where((root_second[1,:] > 0.0))[0][0]
    root2 = mu_root[index2]

    # print statement for roots if desired
    # print('The roots found for X = {} are:'.format(self.X))
    # print('n_1 = {}, n_2 = {}'.format(root1, root2))

    return [root1, root2]

gp_interpolation(mu, kernel=None, center=0.0, sd=1.0)

The function responsible for fitting the coefficients with a GP and predicting at new points. This information will be used in constructing our truncated GP in the function 'Uncertainties'.

Parameters:

Name	Type	Description	Default
mu	ndarray	The chemical potential linspace needed.	required
kernel	obj	The kernel needed for the interpolation GP. Can be fed in from the outside for specific parameter alterations.	None
center	float	Value for the center of the prior.	0.0
sd	float	The scale of the prior.	1.0

Returns:

Name	Type	Description
pred	ndarray	An array of predictions from the GP.
std	ndarray	The standard deviation at the points in 'pred'.
underlying_std	ndarray	The underlying standard deviation of the GP.

Source code in src/neutron_rich_bmm/pqcd.py

    def gp_interpolation(self, mu, kernel=None, center=0.0, sd=1.0):

        '''
        The function responsible for fitting the coefficients with a GP
        and predicting at new points. This information will be used in 
        constructing our truncated GP in the function 'Uncertainties'. 

        Parameters:
            mu (numpy.ndarray): The chemical potential linspace needed.

            kernel (obj): The kernel needed for the interpolation GP. Can be fed in 
                from the outside for specific parameter alterations.

            center (float): Value for the center of the prior. 

            sd (float): The scale of the prior. 

        Returns:
            pred (numpy.ndarray): An array of predictions from the GP.

            std (numpy.ndarray): The standard deviation at the points in 'pred'.

            underlying_std (numpy.ndarray): The underlying standard deviation of the GP.
        '''

#         # grab the correct mu_array from the pressure
#         mu_array, _, _ = self.pressure(mu, scaled=True)

        # reshape the mu_array to what we need
        #mu = mu_array["x2n2"][0]    # need the quark chemical potential
        mu_arr = mu[:, None]

        # interpolate the coefficents using GPs and gsum 
        np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) 
        self.mask = self.gp_mask(mu)

        # Set up gp objects with fixed mean and standard deviation 
        if kernel is None:
            self.kernel = self.gp_kernel()  # in case kernel has not been yet called
        else:
            self.kernel = kernel
            self.gp_interp = gm.ConjugateGaussianProcess(
                kernel=self.kernel, center=center, disp=0, df=3.0, scale=sd, nugget=0) 

        # fit and predict using the interpolated GP
        self.gp_interp.fit(mu_arr[self.mask], self.coeffs[self.mask])
        pred, std = self.gp_interp.predict(mu_arr, return_std=True)
        underlying_std = np.sqrt(self.gp_interp.cov_factor_)

        # print the kernel parameters for viewing outside the function
        print(self.gp_interp.kernel_)
        print(self.gp_interp.cov_factor_)

        return pred, std, underlying_std

gp_kernel()

The kernel that we will use both for interpolating the coefficients and for predicting the truncation error bands.

Returns:

Type	Description
	self.kernel (sklearn object): The kernel needed for the GPs in both 'uncertainties' and 'gp_interpolation'.

Source code in src/neutron_rich_bmm/pqcd.py

def gp_kernel(self):

    '''
    The kernel that we will use both for interpolating the 
    coefficients and for predicting the truncation error bands.

    Parameters:
        None.

    Returns:
        self.kernel (sklearn object): The kernel needed for the GPs 
            in both 'uncertainties' and 'gp_interpolation'. 
    '''

    self.ls = 0.496    # starting guess; can get really close if we set 0.25 and fix it
    self.sd = 0.2    # makes a difference on the band of the regression curve for c_2 
    self.center = 0
    self.nugget = 1e-10  # nugget goes to Cholesky decomp, not the kernel (kernel has own nugget)
    self.kernel = RBF(length_scale=self.ls, length_scale_bounds='fixed') + \
    WhiteKernel(noise_level=self.nugget, noise_level_bounds='fixed') # letting this vary 

    return self.kernel

gp_mask(mu)

The mask array needed to correctly separate our training and testing data.

Parameters:

Name	Type	Description	Default
mu	ndarray	The chemical potential linspace needed.	required

Returns:

Type	Description
	self.mask (numpy.ndarray): The mask for use when interpolating or using the truncated GP.

Source code in src/neutron_rich_bmm/pqcd.py

def gp_mask(self, mu):

    '''
    The mask array needed to correctly separate our training
    and testing data.

    Parameters:
        mu (numpy.ndarray): The chemical potential linspace needed.

    Returns:
        self.mask (numpy.ndarray): The mask for use when interpolating 
            or using the truncated GP.
    '''

    # mask for values above 40*n0 only 
    low_bound = next(i for i, val in enumerate(mu)
                              if val > 0.88616447)
    mu_mask = mu[low_bound:]
    mask_true = np.array([(i) % 25 == 0 for i in range(len(mu_mask))])

    # concatenate with a mask over the other elements of mu before low_bound
    mask_false = np.full((1,len(mu[:low_bound])), fill_value=False)
    self.mask = np.concatenate((mask_false[0], mask_true))

    # old mask
    #self.mask = np.array([(i-3) % 25 == 0 for i in range(len(mu))]) #np.array([(i-6) % 25 == 0 for i in range(len(mu))])

    return self.mask 

inversion(n_mu=None)

Function to invert n(mu) to obtain mu(n).

Parameters:

Name	Type	Description	Default
n_mu	array	Linspace over n_q for the inversion.	None

Returns:

Type	Description
	f_mu_1_result, f_mu_2_result (numpy.array): The two arrays corresponding to the inverted function values for mu(n^(1)) and mu(n^(2)).

Source code in src/neutron_rich_bmm/pqcd.py

def inversion(self, n_mu=None):

    '''
    Function to invert n(mu) to obtain mu(n). 

    Parameters:
        n_mu (numpy.array): Linspace over n_q for the inversion.

    Returns:
        f_mu_1_result, f_mu_2_result (numpy.array): The two arrays 
            corresponding to the inverted function values
            for mu(n^(1)) and mu(n^(2)). 
    '''

    # write the root finding function
    def f_mu_1(n, guess):  
        return opt.fsolve(lambda mu : n - self.pqcd_root.n_1(mu), x0 = guess)[0]

    def f_mu_2(n, guess):  
        return opt.fsolve(lambda mu : n - self.pqcd_root.n_2(mu), x0 = guess)[0]

    # call the function
    f_mu_1_result = []
    f_mu_2_result = []

    if n_mu is None:
        n_mu = np.linspace(0.01, 0.8, 1000)  # n_q right now -> 3*n_B = (0.0, 1.25) GeV^3
    else:
        n_mu = n_mu

    # invert to get mu(n^1) and mu(n^2)
    for i in n_mu:
        f_mu_1_result.append(f_mu_1(i, guess=2.0))   # these results will be in muB/3 (because of n(mu) function)
        f_mu_2_result.append(f_mu_2(i, guess=2.0))

    # convert to n_B here to not make mistake later
    self.nB = n_mu/3.0   # n_q/3

    # convert the mu arrays over (still muB/3)
    f_mu_1_result = np.asarray(f_mu_1_result)
    f_mu_2_result = np.asarray(f_mu_2_result)

    return n_mu, f_mu_1_result, f_mu_2_result

n_first(mu)

The first order in the number density, to be solved to find the root.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potentials.	required

Returns:

Name	Type	Description
n	ndarray	The number density array.

Source code in src/neutron_rich_bmm/pqcd.py

def n_first(self, mu):

    '''The first order in the number density, to be solved
    to find the root.

    Parameters:
        mu (numpy.ndarray): The quark chemical potentials. 

    Returns:
        n (numpy.ndarray): The number density array. 
    '''

    # set up the function
    n = self.pqcd_root.n_1(mu)/self.pqcd_root.n_FG(mu)

    return n

n_second(mu)

The second order in the number density, to be solved to find the root.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potentials.	required

Returns:

Name	Type	Description
n	ndarray	The number density.

Source code in src/neutron_rich_bmm/pqcd.py

def n_second(self, mu):

    '''The second order in the number density, to be solved
    to find the root.

    Parameters:
        mu (numpy.ndarray): The quark chemical potentials. 

    Returns:
        n (numpy.ndarray): The number density. 
    '''

    # set up the function
    n = self.pqcd_root.n_2(mu)/self.pqcd_root.n_FG(mu)

    return n

pressure(mu, n_mu=None, scaled=True)

Calculation of the pressure with respect to the chemical potential using the method of Kurkela et al. (2010). Employs the integration of the number density with respect to mu, using the integration limits determined by the find_roots function.

Parameters:

Name	Type	Description	Default
mu	ndarray	The linspace in chemical potential with which we started.	required
n_mu	ndarray	The corresponding number density.	None
scaled	bool	Toggle to return either scaled pressure (wrt free quark gas) or to return without scaling. Default is True.	True

Returns:

Name	Type	Description
mu_array	dict	Dictionary of the new mu arrays spanning the integration limits of the pressure integration.
pressure_1	ndarray	The mean value of the pressure at first order for the chosen chemical potential (or number density) array.
pressure_2	ndarray	The mean value of the pressure at second order for the chosen chemical potential (or number density) array.

Source code in src/neutron_rich_bmm/pqcd.py

    def pressure(self, mu, n_mu=None, scaled=True):

        '''
        Calculation of the pressure with respect to the chemical
        potential using the method of Kurkela et al. (2010). 
        Employs the integration of the number density with respect
        to mu, using the integration limits determined by the find_roots
        function. 

        Parameters:
            mu (numpy.ndarray): The linspace in chemical potential with 
                which we started.

            n_mu (numpy.ndarray): The corresponding number density.

            scaled (bool): Toggle to return either scaled pressure 
                (wrt free quark gas) or to return without scaling. 
                Default is True.

        Returns:
            mu_array (dict): Dictionary of the new mu arrays spanning 
                the integration limits of the pressure integration. 

            pressure_1 (numpy.ndarray): The mean value of the pressure at 
                first order for the chosen chemical potential 
                (or number density) array. 

            pressure_2 (numpy.ndarray): The mean value of the pressure at 
                second order for the chosen chemical potential 
                (or number density) array. 

        '''

        ### In this version, only using the second order limits 
        ### since they are higher in mu than the first order ones.
        ### A necessary condition for BUQEYE truncation errors at
        ### the present time. 

        # set integration limits first
        roots = self.find_roots()

        if self.X == 2:
            mu_0 = {
                "2" : roots,
            }

        # bag constant
        B = 0.0
        hbarc = 0.197327**3.0

        # set up dicts
        p0 = defaultdict(list)
        P = defaultdict(list)
        mu_array = defaultdict(list)

        # X values and orders keys for dicts
        if self.X == 2:
            mu_array_keys = ["x2n1", "x2n2"]
            mu_keys = ["2", "2"]

            # create each array for each X and order 
            # (using second order chemical potential array)
            mu_array_i = np.concatenate((np.array([mu_0["2"][1]]), mu)) 
            for i,j in zip(mu_array_keys, [0,1]):
                #mu_array[i].append(np.linspace(mu_0["2"][1], max(mu), len(mu))) # concatenate mu_0 with mu_n (OK b/c cutoff always lower)
                mu_array[i].append(mu_array_i)

            # first order, X=2
            for i in mu_array["x2n1"][0]:   # now we are in mu_B
                p0["x2n1"].append(self.pqcd_root.pressure_FG(i)/hbarc)
                P["x2n1"].append(-B + sc.quad((lambda x : (1.0/3.0) * \
                                self.pqcd_root.n_1(x/3.0)), 3.0 * mu_0["2"][0], 3.0 * i)[0])

            # second order, X=2
            for i in mu_array["x2n2"][0]: 
                p0["x2n2"].append(self.pqcd_root.pressure_FG(i)/hbarc)
                P["x2n2"].append(-B + sc.quad((lambda x : (1.0/3.0) * \
                                 self.pqcd_root.n_2(x/3.0)), 3.0 * mu_0["2"][1], 3.0 * i)[0])

            # class variables for p0
            #p0 = {k: v / hbarc for k, v in p0.items()}
            self.p0 = p0    #unscaled, be careful here

#             # re-set mu_array to 100 points we need in density
#             mu_array["x2n1"][0] = mu_array["x2n1"][0][1:]
#             mu_array["x2n2"][0] = mu_array["x2n2"][0][1:]

            if scaled == True:
                p_1 = np.asarray(P["x2n1"])/hbarc/np.asarray(p0["x2n1"])#[1:])
                p_2 = np.asarray(P["x2n2"])/hbarc/np.asarray(p0["x2n2"])#[1:])
                return mu_array, p_1, p_2
            else:
                p_1 = np.asarray(P["x2n1"])/hbarc#[1:])
                p_2 = np.asarray(P["x2n2"])/hbarc#[1:])
                return mu_array, p_1, p_2

        else:
            return None

speed_sound(mu, n)

The speed of sound calculation. Not used in the paper.

Parameters:

Name	Type	Description	Default
mu	ndarray	The chemical potential.	required
n	ndarray	The number density.	required

Returns:

Name	Type	Description
cs2	ndarray	The speed of sound.

Source code in src/neutron_rich_bmm/pqcd.py

def speed_sound(self, mu, n):

    # speed of sound for the Kurkela EOS using
    # inversion as for the pressure 

    '''
    The speed of sound calculation. Not used in 
    the paper. 

    Parameters:
        mu (numpy.ndarray): The chemical potential. 

        n (numpy.ndarray): The number density.

    Returns:
        cs2 (numpy.ndarray): The speed of sound. 
    '''
    import numdifftools as ndt

    # terms
    one = n/mu
    dndmu = ndt.Derivative(self.pqcd_root.n_2, n=1)
    two = dndmu(mu)

    cs2 = one * (two)**(-1.0)   # handles inverse of derivative

    return cs2

uncertainties(mu, n_orders=3, kernel=None, test=None)

Calculation of the truncation error bands for the pQCD EOS, using the Kurkela et al. (2010) formulation for the pressure. This function uses techniques from the gsum package.

Parameters:

Name	Type	Description	Default
mu	ndarray	The linspace of chemical potential needed.	required
n_orders	int	The highest order to which the pressure EOS is calculated.	3
kernel	obj	The kernel needed for the interpolation and truncation GP. Can be fed in from the outside to change parameters.	None
test	ndarray	Testing array.	None

Returns:

Name	Type	Description
data	ndarray	The data array.
		self.coeffs (numpy.ndarray): The values of the coefficents at the chemical potential mu.
std_trunc	ndarray	The arrays of truncation errors per each order.

Source code in src/neutron_rich_bmm/pqcd.py

    def uncertainties(self, mu, n_orders=3, kernel=None, test=None):

        '''
        Calculation of the truncation error bands for the pQCD EOS,
        using the Kurkela et al. (2010) formulation for the pressure.
        This function uses techniques from the gsum package. 

        Parameters:
            mu (numpy.ndarray): The linspace of chemical potential needed.

            n_orders (int): The highest order to which the pressure EOS 
                is calculated.

            kernel (obj): The kernel needed for the interpolation and 
                truncation GP. Can be fed in from the outside to change 
                parameters.

            test (numpy.ndarray): Testing array. 

        Returns:
            data (numpy.ndarray): The data array.

            self.coeffs (numpy.ndarray): The values of the coefficents 
                at the chemical potential mu.

            std_trunc (numpy.ndarray): The arrays of truncation errors 
                per each order.
        '''

#         # grab the correct mu_array from the pressure
#         mu_array, _, _ = self.pressure(mu, scaled=True)

#         # reshape the mu_array to what we need
#         mu = mu_array["x2n2"][0]    # need the quark chemical potential
        mu_arr = mu[:, None]
        self.mu_arr = mu_arr          # for use outside

        # orders
        orders = np.arange(0, n_orders)

        # construct the mask
        self.mask = self.gp_mask(mu)

        # construct the needed data and coefficients
        if test is not None:
            coeffs_all = test
        else:
            coeffs_all = np.array([self.c_0(mu), self.c_1(mu), self.c_2(mu)]).T
        data_all = gm.partials(coeffs_all, ratio=self.expQ(mu_arr), \
                               ref=self.yref(mu_arr), orders=[0,1,2])

        # shape for correct format
        self.coeffs = coeffs_all[:, :n_orders]
        data = data_all[:, :n_orders]

        # call the kernel function 
        if kernel is None:
            self.kernel = self.gp_kernel()
        else:
            self.kernel = kernel

        # set up the truncation GP
        trunc_gp = gm.TruncationGP(kernel=self.kernel, ref=self.yref, \
                            ratio=self.expQ, disp=0, df=3, scale=1, optimizer=None) # disp is variance of the mean
        trunc_gp.fit(mu_arr[self.mask], data[self.mask], orders=orders)

        std_trunc = np.zeros([len(mu_arr), n_orders])
        for i, n in enumerate(orders):
            # Only get the uncertainty due to truncation (kind='trunc')
            _, std_trunc[:,n] = trunc_gp.predict(mu_arr, order=n, return_std=True, kind='trunc')

        return data, self.coeffs, std_trunc

yref(mu)

The reference for the expansion.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential.	required

Returns:

Name	Type	Description
yref	ndarray	The values of yref at the given chemical potentials.

Source code in src/neutron_rich_bmm/pqcd.py

def yref(self, mu):

    '''
    The reference for the expansion. 

    Parameters:
        mu (numpy.ndarray): The quark chemical potential.

    Returns:
        yref (numpy.ndarray): The values of yref at the given
            chemical potentials. 
    '''

    # scaling
    hbarc = 0.197327**3.0

    yref = ((3.0 * mu**4.0) / (4.0 * np.pi**2.0))
    yref = yref[:,0]/hbarc

    return yref

PQCD(X=2, Nf=3)

the rest of the code.

:Example: PQCD(mu=np.linspace(), X=1)

Parameters:

Name	Type	Description	Default
mu	linspace	Range of values to use for mu. Units: GeV.	required
X	int	The value of the coefficient multiplying mu to determine the renormalization scale. Default is 2, as in Kurkela et al. (2010).	2
Nf	int	The number of different quark flavours being considered. Default is 3.	3

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/pqcd.py

def __init__(self, X=2, Nf=3):

    '''
    Initialize the class so that mu and lambda_bar
    have been set and constants have been defined for
    the rest of the code.

    :Example:
        PQCD(mu=np.linspace(), X=1)

    Parameters:
        mu (numpy.linspace): Range of values to use for mu. Units: GeV. 

        X (int): The value of the coefficient multiplying mu
            to determine the renormalization scale. Default is 2,
            as in Kurkela et al. (2010).

        Nf (int): The number of different quark flavours being 
            considered. Default is 3. 

    Returns:
        None.
    '''

    # initialize constants and lambda_bar
    self.X = X
    self.Nc = 3
    self.Nf = Nf
    self.Ca = self.Nc
    self.Cf = (self.Nc**2.0 - 1.0)/(2.0*self.Nc)
    self.lambda_MS = 0.38 #[GeV]

    self.beta0 = (11.0*self.Nc - 2.0*self.Nf)/3.0

    return None

alpha_s(mu, loop=2)

The function for alpha_s with an option for either first order or second order (loops) inclusion.

:Example: PQCD.alpha_s(mu, loop=1)

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential.	required
loop	int	The order at which alpha_s is calculated. Default is 2.	2

Returns:

Name	Type	Description
alpha_s	ndarray	The values of alpha_s for the input chemical potential.

Source code in src/neutron_rich_bmm/pqcd.py

def alpha_s(self, mu, loop=2):

    '''
    The function for alpha_s with an option for either 
    first order or second order (loops) inclusion.

    :Example:
        PQCD.alpha_s(mu, loop=1)

    Parameters:
        mu (numpy.ndarray): The quark chemical potential. 

        loop (int): The order at which alpha_s is calculated. 
            Default is 2.

    Returns:
        alpha_s (numpy.ndarray): The values of alpha_s for the input
            chemical potential.
    '''

    lambda_bar = self.X * mu

    if (loop == 1):
        beta1 = 0.0
    elif (loop == 2):
        beta1 = (17.0/3.0)*self.Ca**2.0 - self.Cf*self.Nf - (5.0/3.0)*(self.Ca*self.Nf)
    else:
        raise ValueError('Loop selected must be 1 or 2.')

    # renormalization scale
    ell = np.log((lambda_bar**2.0)/self.lambda_MS**2.0)

    # running coupling constant
    alpha_s = ((4.0*np.pi)/(self.beta0*ell)) * (1.0 - (2.0*beta1*np.log(ell))/(ell*self.beta0**2.0))      

    return alpha_s

cs2(n, mu)

Calculaton of the speed of sound. Not used in this paper.

Parameters:

Name	Type	Description	Default
n	ndarray	Number density array.	required
mu	ndarray	Corresponding chemical potential array.	required

Returns:

Name	Type	Description
cs2	ndarray	The speed of sound.

Source code in src/neutron_rich_bmm/pqcd.py

def cs2(self, n, mu):

    '''
    Calculaton of the speed of sound. Not used in this
    paper.

    Parameters:
        n (numpy.ndarray): Number density array.

        mu (numpy.ndarray): Corresponding chemical 
            potential array. 

    Returns:
        cs2 (numpy.ndarray): The speed of sound.
    '''

    # mu_FG again
    mu_FG = ((n*np.pi**2.0)/3.0)**(1.0/3.0)

    # determine denominator in mu
    denom = (1.0 - 2.0 * self.alpha_s(mu=mu_FG, loop=2)/np.pi)**(1.0/3.0)

    # now calculate derivative of mu wrt n
    dmudn = (1.0/3.0) * (np.pi**2.0 / (3.0 * denom))**(1.0/3.0) * n**(-2.0/3.0)

    # calculate cs2
    cs2 = (n / mu) * dmudn 

    return cs2

mask_array(array, neg=False, fill_value=None) staticmethod

Returns a masked array from an original array that contains unwanted nan values. Can also fill values in the place of the mask, if desired.

:Example: PQCD.mask_array(array=mu, fill_value=0)

Parameters:

Name	Type	Description	Default
array	ndarray	The array with values to mask.	required
neg	bool	If False, will not also mask negative values. If True, will check for negative values and mask using current fill value.	False
fill_value	(int, float)	The value with which to fill the mask. Default is None.	None

Returns:

Name	Type	Description
masked_array	ndarray	The masked array with or without filled values.

Source code in src/neutron_rich_bmm/pqcd.py

@staticmethod
def mask_array(array, neg=False, fill_value=None):

    '''
    Returns a masked array from an original array that contains 
    unwanted nan values. Can also fill values in the place of 
    the mask, if desired.

    :Example:
        PQCD.mask_array(array=mu, fill_value=0)

    Parameters:
        array (numpy.ndarray): The array with values to mask. 

        neg (bool): If False, will not also mask negative values. If True,
            will check for negative values and mask using current
            fill value. 

        fill_value (int, float): The value with which to fill the mask. 
            Default is None.

    Returns:
        masked_array (numpy.ndarray): The masked array with or without 
            filled values.
    '''

    # initialize the mask and new array
    mask = np.zeros(len(array), dtype=bool)
    masked_array = np.zeros(len(array))

    # check for neg bool and mask if needed
    if neg is True:
        for i in range(len(array)):
            if array[i] < 0.0:
                array[i] = fill_value

    # check for nan values and replace in mask
    for i in range(len(array)):
        if np.isnan(array[i]):
            mask[i] = True

    # mask the original array
    masked_array = np.ma.array(array, mask=mask)

    # if there were given fill values, fill the array
    if fill_value is not None:
        masked_array = np.ma.filled(masked_array, fill_value=fill_value)

    return masked_array

n_1(mu)

The first order (NLO) term in the number density of Eq. (59) in Kurkela et al. (2010).

:Example: PQCD.n_1(self, mu=np.linspace())

Parameters:

Name	Type	Description	Default
mu	linspace	The original mu range.	required

Returns:

Name	Type	Description
n_1	ndarray	The first order contribution to the number density.

Source code in src/neutron_rich_bmm/pqcd.py

def n_1(self, mu):

    '''
    The first order (NLO) term in the number density of Eq. (59) in
    Kurkela et al. (2010).

    :Example:
        PQCD.n_1(self, mu=np.linspace())

    Parameters:
        mu (numpy.linspace): The original mu range.

    Returns:
        n_1 (numpy.ndarray): The first order contribution to the 
            number density.
    '''

    n_1 = self.n_FG(mu) * (1.0 - (2.0/np.pi) * self.alpha_s(mu, loop=2))  

    return n_1

n_2(mu)

The second order (NNLO) term in the number density in Eq. (61) of Kurkela et al. (2010).

:Example: PQCD.n_2(mu=np.linspace())

Parameters:

Name	Type	Description	Default
mu	linspace	The original mu range.	required

Returns:

Name	Type	Description
n_2	ndarray	The second order contribution to the number density.

Source code in src/neutron_rich_bmm/pqcd.py

def n_2(self, mu):

    '''
    The second order (NNLO) term in the number density in Eq. (61)
    of Kurkela et al. (2010). 

    :Example:
        PQCD.n_2(mu=np.linspace())

    Parameters:
        mu (numpy.linspace): The original mu range.

    Returns:
        n_2 (numpy.ndarray): The second order contribution to 
            the number density.
    '''

    lambda_bar = self.X * mu

    outer = self.n_FG(mu)
    one = 1.0
    two = 2.0*self.alpha_s(mu, loop=2)/np.pi
    three = (self.alpha_s(mu, loop=2)/np.pi)**2.0
    four = (61.0/4.0) - 11.0*np.log(2) - 0.369165*self.Nf 
    five = self.Nf*np.log(self.Nf*self.alpha_s(mu, loop=2)/np.pi) 
    six = self.beta0*np.log(lambda_bar/mu)

    n_2 = outer * (one - two - three*(four + five + six))

    return n_2

n_FG(mu)

The zeroth order (LO) term in the number density, from Eq. (59) in Kurkela et al. (2010). Aka the FG contribution.

:Example: PQCD.n_FG(mu=np.linspace())

Parameters:

Name	Type	Description	Default
mu	linspace	The original mu range.	required

Returns:

Name	Type	Description
n_FG	ndarray	The FG contribution to the number density.

Source code in src/neutron_rich_bmm/pqcd.py

def n_FG(self, mu):

    '''
    The zeroth order (LO) term in the number density, from Eq. (59)
    in Kurkela et al. (2010). Aka the FG contribution.

    :Example:
        PQCD.n_FG(mu=np.linspace())

    Parameters:
        mu (numpy.linspace): The original mu range. 

    Returns:
        n_FG (numpy.ndarray): The FG contribution to the number 
            density.
    '''

    n_FG = (3.0 * mu**3.0) / np.pi**2.0

    return n_FG

pressure_FG(mu)

The FG contribution to the pressure in terms of mu.

:Example: PQCD.pressure_FG(mu=np.linspace())

Parameters:

Name	Type	Description	Default
mu	linspace	The original range in mu.	required

Returns:

Name	Type	Description
p_FG	linspace	The zeroth order (LO) contribution the pressure.

Source code in src/neutron_rich_bmm/pqcd.py

def pressure_FG(self, mu):

    '''
    The FG contribution to the pressure in terms of mu.

    :Example:
        PQCD.pressure_FG(mu=np.linspace())

    Parameters:
        mu (numpy.linspace): The original range in mu.

    Returns:
        p_FG (numpy.linspace): The zeroth order (LO) contribution 
            the pressure. 
    '''

    p_FG = (3.0 * mu**4.0) / (4.0 * np.pi**2.0)

    return p_FG

pressure_mu(mu, order=2)

The pressure with respect to mu for either NLO or NNLO. Note: alpha_s has been used up to second order here to maintain the validity of the value of alpha_s at 2 GeV.

:Example: PQCD.pressure_mu(mu=np.linspace(), order=1)

Parameters:

Name	Type	Description	Default
mu	linspace	The original range of mu.	required
order	int	The order at which the pressure is calculated. Default is 2, options are 1 and 2.	2

Returns:

Name	Type	Description
pressure	ndarray	The pressure as a function of mu.

Source code in src/neutron_rich_bmm/pqcd.py

def pressure_mu(self, mu, order=2):

    '''
    The pressure with respect to mu for either NLO or NNLO.
    Note: alpha_s has been used up to second order here to
    maintain the validity of the value of alpha_s at 2 GeV.

    :Example:
        PQCD.pressure_mu(mu=np.linspace(), order=1)

    Parameters:
        mu (numpy.linspace): The original range of mu.

        order (int): The order at which the pressure is calculated. 
            Default is 2, options are 1 and 2.

    Returns:
        pressure (numpy.ndarray): The pressure as a function of mu.
    '''

    lambda_bar = self.X * mu 

    if order == 1:

        pressure = self.pressure_FG(mu) * (1.0 - 2.0 * \
            self.alpha_s(mu, loop=2)/np.pi)

        return pressure

    elif order == 2:

        first = 2.0 * self.alpha_s(mu, loop=2) / np.pi
        second = 0.303964 * self.alpha_s(mu, loop=2)**2.0 \
            * np.log(self.alpha_s(mu, loop=2))
        second_2 = self.alpha_s(mu, loop=2)**2.0 * (0.874355 \
            + 0.911891 * np.log(lambda_bar/mu))
        pressure = self.pressure_FG(mu) * (1.0 - first - \
            second - second_2)

        return pressure

pressure_n(n, alpha_s_mu=None)

The pressure as a function of number density, to first order.

:Example:
PQCD.pressure_n(n=np.linspace(), alpha_s_mu=None)

Parameters:

Name	Type	Description	Default
n	linspace	The range in number density being considered.	required
alpha_s_mu	(linspace, ndarray)	A possible different chemical potential range to send to alpha_s. Default is None.	None

Returns:

Name	Type	Description
p_n	ndarray	The pressure as a function of density, at NLO.

Source code in src/neutron_rich_bmm/pqcd.py

def pressure_n(self, n, alpha_s_mu=None):

    '''
    The pressure as a function of number density, to first
    order.

    :Example:  
        PQCD.pressure_n(n=np.linspace(), alpha_s_mu=None)

    Parameters:
        n (numpy.linspace): The range in number density being 
            considered.

        alpha_s_mu (numpy.linspace, numpy.ndarray): A possible 
            different chemical potential range to send to alpha_s. 
            Default is None.

    Returns:
        p_n (numpy.ndarray): The pressure as a function of density, 
            at NLO.
    '''

    if alpha_s_mu is not None:
        p_n = self.pressure_n_FG(n) * \
            (1.0 + 2.0 * self.alpha_s(mu=alpha_s_mu, loop=2)/(3.0 * np.pi))**(4.0) \
            * (1.0 - 2.0 * self.alpha_s(mu=alpha_s_mu, loop=2))
    else:
        raise ValueError('Supply a value of mu for alpha_s.')

    return p_n

pressure_n3lo(mu, scaled=False)

A first attempt at coding the results from Gorda et al. (2023) including N3LO.

Parameters:

Name	Type	Description	Default
mu	ndarray	The quark chemical potential.	required
scaled	bool	Whether the data is scaled or not. Default is False.	False

Returns:

Name	Type	Description
p_n3lo	ndarray	The pressure at N3LO.

Source code in src/neutron_rich_bmm/pqcd.py

def pressure_n3lo(self, mu, scaled=False):

    '''
    A first attempt at coding the results from 
    Gorda et al. (2023) including N3LO. 

    Parameters:
        mu (numpy.ndarray): The quark chemical potential. 

        scaled (bool): Whether the data is scaled or not. 
            Default is False. 

    Returns:
        p_n3lo (numpy.ndarray): The pressure at N3LO. 
    '''

    # define constants and pieces
    lambda_bar = self.X * mu 
    c0 = -23.0 # 68% uncertainty of +/- 10.0

    # define constants
    c_32 = 11.0/12.0
    c_31 = -6.5968 - 3.0 * np.log(lambda_bar/(2.0*mu))
    c_30 = 5.1342 + (2.0/3.0)*c0 - 18.284 * np.log(lambda_bar/(2.0*mu)) \
        - (9.0/2.0) * np.log(lambda_bar/(2.0*mu))**2.0

    p_free = 3.0 * mu**4.0 / (4.0 * np.pi**2.0)

    # define the equation (with and without scaling)
    one = np.log(3.0 * self.alpha_s(mu,loop=2)/np.pi) + \
            3.0 * np.log(lambda_bar/(2.0*mu)) + 5.0021 
    two = c_32 * np.log(3.0*self.alpha_s(mu,loop=2)/np.pi)**2.0 + \
            c_31 * np.log(3.0*self.alpha_s(mu,loop=2)/np.pi) + c_30
    pressure = 1.0 - 2.0 * (self.alpha_s(mu, loop=2)/np.pi) - \
            3.0 * (self.alpha_s(mu, loop=2)/np.pi)**2.0 * one + 9.0 * (self.alpha_s(mu,loop=2)/np.pi)**3.0 * two

    if scaled is True:
        p_n3lo = pressure 
    else:
        p_n3lo = p_free * pressure

    return p_n3lo

pressure_n_FG(n)

FG pressure with respect to number density n.

Parameters:

Name	Type	Description	Default
n	ndarray	The baryon number density.	required

Returns:

Name	Type	Description
p_n_FG	ndarray	The FG pressure as a function of density.

Source code in src/neutron_rich_bmm/pqcd.py

def pressure_n_FG(self, n):

    '''
    FG pressure with respect to number density n.

    Parameters:
        n (numpy.ndarray): The baryon number density. 

    Returns:
        p_n_FG (numpy.ndarray): The FG pressure as a function of density.
    '''

    p_n_FG = (3.0 / 4.0 * np.pi**2.0) * (n * np.pi**2.0 / 3.0)**(4.0/3.0)

    return p_n_FG

Truncation(x, x_FG, norders, orders, yref, expQ, coeffs, mask=None)

The truncation error class that wraps the gsum package. For use on the pQCD EOS results.

Parameters:

Name	Type	Description	Default
x	array	Quark chemical potential.	required
x_FG	2D array	FG quark chemical potential.	required
norders	int	Number of orders used.	required
orders	list	The list of orders ([0 1 2], etc.)	required
yref	function	Functional form for yref.	required
expQ	function	Functional form for the expansion parameter.	required
coeffs	array	The coefficient array of arrays. Must be transposed when sent in to this class.	required
mask	bool array	If using a mask, send the mask in here.	None

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/truncation_error.py

def __init__(self, x, x_FG, norders, orders, yref, expQ, coeffs, mask=None):

    # immediately create this 
    self.orders_mask = None
    self.orders = orders
    self.n_orders = norders

    # define class variables for pressure to all orders
    self.x = x
    self.X = x[:,None]
    if x_FG is not None:
        self.x_FG = x_FG
        self.X_FG = x_FG[:,None]
    else:
        self.x_FG = None
        self.X_FG = None
    self.norders = norders

    # declare yref, Q
    self.yref = yref(self.X) 
    self.expQ = expQ(self.X)

    ### formalism from gsum docs ###

    # separate the coeffs
    self.coeffs_list = []

    # if masking, create separate instance of coeffs
    if mask is not None:
        self.coeffs_list_trunc = []
        self.orders_mask = mask
        coeffs_trunc = coeffs[:,self.orders_mask]

        for i in range(len(coeffs_trunc.T)):
            self.coeffs_list_trunc.append(coeffs_trunc[:,i])

    for i in range(len(coeffs.T)):
        self.coeffs_list.append(coeffs[:,i])

    # construct total arrays of each quantity
    self.coeffs_all = np.array(self.coeffs_list).T
    self.data_all = gm.partials(self.coeffs_all, ratio=self.expQ, ref=self.yref, orders=[range(norders)])
    self.diffs_all = np.array([self.data_all[:, 0], *np.diff(self.data_all, axis=1).T]).T

    # save different coeff array for the GP to fit
    if mask is not None:
        self.coeffs_all_trunc = np.array(self.coeffs_list_trunc).T

    # get the "all-orders" curve
    self.data_true = self.data_all[:, -1]

    # specify range
    if mask is not None:
        self.coeffs_trunc = self.coeffs_all_trunc

    self.coeffs = self.coeffs_all[:, :norders]
    self.data = self.data_all[:, :norders]
    self.diffs = self.diffs_all[:, :norders]

    return None

diagnostics(dx_train=30, dx_test=15)

The diagnostic function to check the validity of the truncation error obtained via gsum. Uses gsum to perform Mahalanobis distance and pivoted Cholesky calculations. Plots the results.

Parameters:

Name	Type	Description	Default
dx_train	int	The number to use as a step size for the training data.	30
dx_test	int	The number to use as a step size for the testing data.	15

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/truncation_error.py

def diagnostics(self, dx_train=30, dx_test=15):

    r'''
    The diagnostic function to check the validity of 
    the truncation error obtained via gsum. Uses
    gsum to perform Mahalanobis distance and pivoted
    Cholesky calculations. Plots the results.

    Parameters:
        dx_train (int): The number to use as a step size for the
            training data.

        dx_test (int): The number to use as a step size for the
            testing data.

    Returns:
        None.
    '''

    # set the plot labels
    MD_label = r'$\mathrm{D}_{\mathrm{MD}}^2$'
    PC_label = r'$\mathrm{D}_{\mathrm{PC}}$'

    # set up plotting tools
    mpl.rcParams['text.usetex'] = True
    mpl.rcParams['figure.dpi'] = 150   # change for paper plots
    mpl.rcParams['font.size'] = 8
    mpl.rcParams['ytick.direction'] = 'in'
    mpl.rcParams['xtick.direction'] = 'in'
    mpl.rcParams['xtick.labelsize'] = 8
    mpl.rcParams['ytick.labelsize'] = 8
    WIDE_IMG_WIDTH = 800
    NARROW_IMG_WIDTH = 400

    cmaps = [plt.get_cmap(name) for name in ['Oranges', 'Greens', 'Blues', 'Reds']]
    colors = [cmap(0.55 - 0.1 * (i==0)) for i, cmap in enumerate(cmaps)]
    light_colors = [cmap(0.25) for cmap in cmaps]

    edgewidth = 0.6
    text_bbox = dict(boxstyle='round', fc=(1, 1, 1, 0.6), ec='k', lw=0.8)

    # call the masking function for diagnostics
    x_train_mask, x_valid_mask = self.regular_train_test_split(self.x, dx_train, dx_test, \
                                                               offset_train=0, \
                                                               offset_test=0, xmin=0.88616447)

    # print the values for checking (use only ones after 40 n0 again)
    print('Number of points in training set:', np.shape(self.X[self.mask])[0])
    print('Number of points in validation set:', np.shape(self.X[x_valid_mask])[0])

    print('\nTraining set: \n', self.X[self.mask])
    print('\nValidation set: \n', self.X[x_valid_mask])

    # overwrite training mask with the original mask for keeping range correct
    x_train_mask = self.mask 

    # check if the two arrays have equal elements
    for i in self.X[x_train_mask]:
        for j in self.X[x_valid_mask]:
            if i == j:
                print('Found an equal value!')

    # already fit the GP kernel, do diagnostics directly now 
    underlying_std = np.sqrt(self.trunc_gp.coeffs_process.cov_factor_)
    print(underlying_std)
    print(self.trunc_gp.coeffs_process.nugget)

    gp_diagnostic = self.trunc_gp.coeffs_process  # set equal here (same object)
    print('Calculated value :', gp_diagnostic.df_ * gp_diagnostic.scale_**2 / (gp_diagnostic.df_ + 2))

    # Print out the kernel of the fitted GP
    print('Trained kernel: ', self.trunc_gp.coeffs_process.kernel_)
    print('Scale: ', gp_diagnostic.scale_)

    # MD diagnostic plotting
    mean_underlying = gp_diagnostic.mean(self.X[x_valid_mask])
    cov_underlying = gp_diagnostic.cov(self.X[x_valid_mask])
    print('Condition number:', np.linalg.cond(cov_underlying))

    # coeffs coming from initial set up
    gdgn = gm.GraphicalDiagnostic(self.coeffs[x_valid_mask], mean_underlying, cov_underlying, \
                                  colors=colors,
                                  gray='gray', black='k')

    def offset_xlabel(ax):
        ax.set_xticks([0])
        ax.set_xticklabels(labels=[0], fontdict=dict(color='w'))
        ax.tick_params(axis='x', length=0)
        return ax

    fig, ax = plt.subplots(figsize=(1, 3.2))
    ax = gdgn.md_squared(type='box', trim=False, title=None, xlabel=MD_label)
    offset_xlabel(ax)
    ax.set_ylim(0, 20)

    # Pivoted Cholesky as well
    with plt.rc_context({"text.usetex": True}):
        fig, ax = plt.subplots(figsize=(3.2, 3.2))
        gdgn.pivoted_cholesky_errors(ax=ax, title=None)
        ax.text(0.04, 0.967, PC_label, bbox=text_bbox, transform=ax.transAxes, va='top', ha='left')
        plt.show()

    return None

gp_interpolation(center=0.0, sd=1.0)

The function responsible for fitting the coefficients with a GP and predicting at new points. This information will be used in constructing our truncated GP in the function 'Uncertainties'.

Parameters:

Name	Type	Description	Default
center	float	The center value for the prior.	0.0
sd	float	The scale for the prior.	1.0

Returns:

Name	Type	Description
pred	array	An array of predictions from the GP.
std	array	The standard deviation at the points in 'pred'.
underlying_std	array	The underlying standard deviation of the GP.

Source code in src/neutron_rich_bmm/truncation_error.py

def gp_interpolation(self, center=0.0, sd=1.0):

    r'''
    The function responsible for fitting the coefficients with a GP
    and predicting at new points. This information will be used in 
    constructing our truncated GP in the function 'Uncertainties'. 

    Parameters:
        center (float): The center value for the prior.

        sd (float): The scale for the prior.

    Returns:
        pred (array): An array of predictions from the GP.

        std (array): The standard deviation at the points in 'pred'.

        underlying_std (array): The underlying standard deviation of
            the GP.
    '''

    # interpolate the coefficents using GPs and gsum 
    np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) 

    # call the mask function to set this for interpolation training
    self.mask = self.gp_mask(self.x)

    # Set up gp objects with fixed mean and standard deviation 
    kernel = self.gp_kernel(ls=3.0, sd=sd, center=0)  
    self.gp_interp = gm.ConjugateGaussianProcess(
        kernel=kernel, center=center, disp=0, df=3.0, scale=sd, nugget=0) 

    # fit and predict using the interpolated GP
    if self.orders_mask is not None:
        self.gp_interp.fit(self.X[self.mask], self.coeffs_trunc[self.mask])
    else:
        self.gp_interp.fit(self.X[self.mask], self.coeffs[self.mask])
    pred, std = self.gp_interp.predict(self.X, return_std=True)
    underlying_std = np.sqrt(self.gp_interp.cov_factor_)

    # print the kernel parameters for viewing outside the function
    print(self.gp_interp.kernel_)
    print(self.gp_interp.cov_factor_)

    return pred, std, underlying_std

gp_kernel(ls=3.0, sd=0.5, center=0, nugget=1e-10)

The kernel that we will use both for interpolating the coefficients and for predicting the truncation error bands. This one is unfixed, so the value of the ls obtained here will be used to fix the second run when calling params attribute.

Parameters:

Name	Type	Description	Default
ls	float	The lengthscale guess for the kernel.	3.0
sd	float	The scale for the prior.	0.5
center	float	The center value for the prior.	0
nugget	(int, float)	The value of the nugget to send to the Cholesky decomposition.	1e-10

Returns:

Name	Type	Description
kernel	sklearn object	The kernel needed for the GPs.

Source code in src/neutron_rich_bmm/truncation_error.py

def gp_kernel(self, ls=3.0, sd=0.5, center=0, nugget=1e-10):

    r'''
    The kernel that we will use both for interpolating the 
    coefficients and for predicting the truncation error bands.
    This one is unfixed, so the value of the ls obtained here will 
    be used to fix the second run when calling params attribute.

    Parameters:
        ls (float): The lengthscale guess for the kernel.

        sd (float): The scale for the prior.

        center (float): The center value for the prior.

        nugget (int, float): The value of the nugget to send 
            to the Cholesky decomposition.

    Returns:
        kernel (sklearn object): The kernel needed for the GPs. 
    '''

    self.ls = ls    # starting guess; can get really close if we set 0.25 and fix it
    self.sd = sd    # makes a difference on the band of the regression curve for c_2 
    self.center = center
    self.nugget = nugget # nugget for the Cholesky, not the kernel 
    kernel = RBF(length_scale=self.ls)# + \
   # WhiteKernel(noise_level=self.nugget, noise_level_bounds='fixed')

    return kernel

gp_mask(mu)

The mask array needed to correctly separate our training and testing data.

Parameters:

Name	Type	Description	Default
mu	array	The chemical potential linspace needed.	required

Returns:

Type	Description
	self.mask (array): The mask for use when interpolating or using the truncated GP.

Source code in src/neutron_rich_bmm/truncation_error.py

    def gp_mask(self, mu):

        r'''
        The mask array needed to correctly separate our training
        and testing data.

        Parameters:
            mu (array): The chemical potential linspace needed.

        Returns:
            self.mask (array): The mask for use when interpolating 
                or using the truncated GP.
        '''

        # mask for values above 40*n0 only (good)
        low_bound = next(i for i, val in enumerate(mu)
                                  if val > 0.946136639) # value of mean of mu_FG and mu_n at 40*n0
        mu_mask = mu[low_bound:]
#        mu_mask = mu    # no mask applied

        # set the mask s.t. it picks the same training point number each time (this needs to be an input from outside)
        mask_num = len(mu_mask) // 1.5   # 3.5 for the n(mu) calculation
        mask_true = np.array([(i) % mask_num == 0 for i in range(len(mu_mask))])  # i-3 for <40n0

        # concatenate with a mask over the other elements of mu before low_bound
        mask_false = np.full((1,len(mu[:low_bound])), fill_value=False)
        self.mask = np.concatenate((mask_false[0], mask_true))
#        self.mask = mask_true     #no mask applied

        return self.mask 

uncertainties(data=None, expQ=None, yref=None, sd=0.5, nugget=1e-10, excluded=None)

Calculation of the truncation error bands for the pQCD EOS, using the Gorda et al. (2021) formulation for the pressure. This function uses techniques from the gsum package.

Parameters:

Name	Type	Description	Default
data	array	The data given in an array of arrays for each order-by-order result.	None
expQ	function	The functional form of the expansion parameter for gsum to use.	None
yref	function	The functional form of yref for gsum to use.	None
sd	float	The scale for the prior.	0.5
nugget	(int, float)	The nugget for the Cholesky decomposition.	1e-10
excluded	list	The orders we wish to exclude from training on in the coefficient arrays. Default is None.	None

Returns:

Name	Type	Description
data	array	The data array, containing partials at each order.
		self.coeffs (array): The values of the coefficents at x.
std_trunc	array	The arrays of truncation errors per each order.

Source code in src/neutron_rich_bmm/truncation_error.py

def uncertainties(self, data=None, expQ=None, yref=None, sd=0.5, nugget=1e-10, excluded=None):

    r'''
    Calculation of the truncation error bands for the pQCD EOS, using 
    the Gorda et al. (2021) formulation for the pressure.
    This function uses techniques from the gsum package. 

    Parameters:
        data (array): The data given in an array of arrays for 
            each order-by-order result.

        expQ (function): The functional form of the expansion 
            parameter for gsum to use.

        yref (function): The functional form of yref for gsum 
            to use.

        sd (float): The scale for the prior.

        nugget (int, float): The nugget for the Cholesky 
            decomposition.

        excluded (list): The orders we wish to exclude from training
            on in the coefficient arrays. Default is None.

    Returns:
        data (array): The data array, containing partials at each 
            order.

        self.coeffs (array): The values of the coefficents at x.

        std_trunc (array): The arrays of truncation errors per each 
            order.
    '''

    # construct mask
    self.mask = self.gp_mask(self.x)

    # get correct data shape
    if data is None:         
        data = self.data_all[:, :self.n_orders]
    else:
        data = data[:, :self.n_orders]

    # save nugget
    self.nugget = nugget

    # set up the truncation GP
    self.kernel = self.gp_kernel(ls=3.0, sd=sd, center=0, nugget=self.nugget) 
    self.trunc_gp = gm.TruncationGP(kernel=self.kernel, ref=yref, \
                        ratio=expQ, disp=0, df=3.0, scale=sd, excluded=excluded, nugget=self.nugget)

    self.trunc_gp.fit(self.X[self.mask], data[self.mask], orders=self.orders)

    std_trunc = np.zeros([len(self.X), self.n_orders])
    cov_trunc = np.zeros([len(self.X), len(self.X), self.n_orders])
    for i, n in enumerate(self.orders):
        # Only get the uncertainty due to truncation (kind='trunc')
        _, std_trunc[:,n] = self.trunc_gp.predict(self.X, order=n, return_std=True, \
                                                  kind='trunc', pred_noise=True)
        _, cov_trunc[:,:,n] = self.trunc_gp.predict(self.X, order=n, return_std=False, \
                                                    return_cov=True, kind='trunc', pred_noise=True)

    # external access without altering return
    self.cov_trunc = cov_trunc

    # check the kernel hyperparameters
    print(self.trunc_gp.coeffs_process.kernel_)
    print(self.trunc_gp.coeffs_process.nugget)

    return data, self.coeffs, std_trunc, cov_trunc

TruncationDens(nb, norders, orders, yref, expQ, coeffs, mask=None)

The truncation error class that wraps the gsum package. For use on the pQCD EOS results, in density (KLW).

Parameters:

Name	Type	Description	Default
nb	array	Baryon number density.	required
norders	int	Number of orders used.	required
orders	list	The list of orders ([0 1 2], etc.)	required
yref	function	Functional form for yref, dependent on the number density.	required
expQ	function	Functional form for the expansion parameter, dependent on the number density.	required
coeffs	array	The coefficient array of arrays. Must be transposed when sent in to this class. Dependent on number density.	required
mask	boolean array	If using a mask over some orders, send the mask in here.	None

Returns:

Type	Description
	None.

Source code in src/neutron_rich_bmm/truncation_error_dens.py

def __init__(self, nb, norders, orders, yref, expQ, coeffs, mask=None):

    # immediately create this 
    self.orders_mask = None
    self.orders = orders
    self.n_orders = norders

    # define class variables for pressure to all orders
    self.nb = nb
    self.nB = nb[:,None]
    self.norders = norders

    # declare yref, Q
    self.yref = yref(self.nB) 
    self.expQ = expQ(self.nB)

    # separate the coeffs
    self.coeffs_list = []

    # if masking, create separate instance of coeffs
    if mask is not None:
        self.coeffs_list_trunc = []
        self.orders_mask = mask
        coeffs_trunc = coeffs[:,self.orders_mask]

        for i in range(len(coeffs_trunc.T)):
            self.coeffs_list_trunc.append(coeffs_trunc[:,i])

    for i in range(len(coeffs.T)):
        self.coeffs_list.append(coeffs[:,i])

    # construct total arrays of each quantity
    self.coeffs_all = np.array(self.coeffs_list).T
    self.data_all = gm.partials(self.coeffs_all, ratio=self.expQ, ref=self.yref, orders=[range(norders)])
    self.diffs_all = np.array([self.data_all[:, 0], *np.diff(self.data_all, axis=1).T]).T

    # save different coeff array for the GP to fit
    if mask is not None:
        self.coeffs_all_trunc = np.array(self.coeffs_list_trunc).T

    # get the "all-orders" curve
    self.data_true = self.data_all[:, -1]

    # specify range
    if mask is not None:
        self.coeffs_trunc = self.coeffs_all_trunc

    # "known" orders
    self.coeffs = self.coeffs_all[:, :norders]
    self.data = self.data_all[:, :norders]
    self.diffs = self.diffs_all[:, :norders]

    return None

diagnostics(dx_train=30, dx_test=15)

The diagnostic function to check the validity of the truncation error obtained via gsum. Uses gsum to perform Mahalanobis distance and pivoted Cholesky calculations. Plots the results.

Parameters:

Name	Type	Description	Default
dx_train	int	The number to use as a step size for the training data.	30
dx_test	int	The number to use as a step size for the testing data.	15

Returns:

Type	Description
	None

Source code in src/neutron_rich_bmm/truncation_error_dens.py

def diagnostics(self, dx_train=30, dx_test=15):

    '''
    The diagnostic function to check the validity of 
    the truncation error obtained via gsum. Uses
    gsum to perform Mahalanobis distance and pivoted
    Cholesky calculations. Plots the results.

    Parameters:
        dx_train (int): The number to use as a step size for the
            training data.

        dx_test (int): The number to use as a step size for the
            testing data.

    Returns:
        None
    '''

    # set the plot labels
    MD_label = r'$\mathrm{D}_{\mathrm{MD}}^2$'
    PC_label = r'$\mathrm{D}_{\mathrm{PC}}$'

    # set up plotting tools
    mpl.rcParams['text.usetex'] = True
    mpl.rcParams['figure.dpi'] = 150   # change for paper plots
    mpl.rcParams['font.size'] = 8
    mpl.rcParams['ytick.direction'] = 'in'
    mpl.rcParams['xtick.direction'] = 'in'
    mpl.rcParams['xtick.labelsize'] = 8
    mpl.rcParams['ytick.labelsize'] = 8
    WIDE_IMG_WIDTH = 800
    NARROW_IMG_WIDTH = 400

    cmaps = [plt.get_cmap(name) for name in ['Oranges', 'Greens', 'Blues', 'Reds']]
    colors = [cmap(0.55 - 0.1 * (i==0)) for i, cmap in enumerate(cmaps)]
    light_colors = [cmap(0.25) for cmap in cmaps]

    edgewidth = 0.6
    text_bbox = dict(boxstyle='round', fc=(1, 1, 1, 0.6), ec='k', lw=0.8)

    # call the masking function for diagnostics
    nb_train_mask, nb_valid_mask = self.regular_train_test_split(self.nb, dx_train, dx_test, \
                                                               offset_train=0, \
                                                               offset_test=0, xmin=0.88616447)

    # print the values for checking (use only ones after 40 n0 again)
    print('Number of points in training set:', np.shape(self.nB[self.mask])[0])
    print('Number of points in validation set:', np.shape(self.nB[nb_valid_mask])[0])

    print('\nTraining set: \n', self.nB[self.mask])
    print('\nValidation set: \n', self.nB[nb_valid_mask])

    # overwrite training mask with the original mask for keeping range correct
    nb_train_mask = self.mask 

    # check if the two arrays have equal elements
    for i in self.nB[nb_train_mask]:
        for j in self.nB[nb_valid_mask]:
            if i == j:
                print('Found an equal value!')

    # already fit the GP kernel, do diagnostics directly now 
    underlying_std = np.sqrt(self.trunc_gp.coeffs_process.cov_factor_)
    print(underlying_std)
    print(self.trunc_gp.coeffs_process.nugget)

    gp_diagnostic = self.trunc_gp.coeffs_process  # set equal here (same object)
    print('Calculated value :', gp_diagnostic.df_ * gp_diagnostic.scale_**2 / (gp_diagnostic.df_ + 2))

    # Print out the kernel of the fitted GP
    print('Trained kernel: ', self.trunc_gp.coeffs_process.kernel_)
    print('Scale: ', gp_diagnostic.scale_)

    # MD diagnostic plotting
    mean_underlying = gp_diagnostic.mean(self.nB[nb_valid_mask])
    cov_underlying = gp_diagnostic.cov(self.nB[nb_valid_mask])
    print('Condition number:', np.linalg.cond(cov_underlying))

    # coeffs coming from initial set up
    gdgn = gm.GraphicalDiagnostic(self.coeffs[nb_valid_mask], mean_underlying, cov_underlying, \
                                  colors=colors,
                                  gray='gray', black='k')

    def offset_xlabel(ax):
        ax.set_xticks([0])
        ax.set_xticklabels(labels=[0], fontdict=dict(color='w'))
        ax.tick_params(axis='x', length=0)
        return ax

    fig, ax = plt.subplots(figsize=(1, 3.2))
    ax = gdgn.md_squared(type='box', trim=False, title=None, xlabel=MD_label)
    offset_xlabel(ax)
    ax.set_ylim(0, 20)

    # Pivoted Cholesky as well
    with plt.rc_context({"text.usetex": True}):
        fig, ax = plt.subplots(figsize=(3.2, 3.2))
        gdgn.pivoted_cholesky_errors(ax=ax, title=None)
        ax.text(0.04, 0.967, PC_label, bbox=text_bbox, transform=ax.transAxes, va='top', ha='left')
        plt.show()

    return None

gp_interpolation(center=0.0, dof=3.0, sd=1.0)

The function responsible for fitting the coefficients with a GP and predicting at new points. This information will be used in constructing our truncated GP in the function 'Uncertainties'.

Parameters:

Name	Type	Description	Default
center	float	The center value for the prior.	0.0
sd	float	The scale for the prior.	1.0

Returns:

Name	Type	Description
pred	ndarray	An array of predictions from the GP.
std	ndarray	The standard deviation at the points in 'pred'.
underlying_std	ndarray	The underlying standard deviation of the GP.

Source code in src/neutron_rich_bmm/truncation_error_dens.py

def gp_interpolation(self, center=0.0, dof=3.0, sd=1.0):

    '''
    The function responsible for fitting the coefficients with a GP
    and predicting at new points. This information will be used in 
    constructing our truncated GP in the function 'Uncertainties'. 

    Parameters:
        center (float): The center value for the prior.

        sd (float): The scale for the prior.

    Returns:
        pred (numpy.ndarray): An array of predictions from the GP.

        std (numpy.ndarray): The standard deviation at the points in 
            'pred'.

        underlying_std (numpy.ndarray): The underlying standard deviation 
            of the GP.
    '''

    # interpolate the coefficents using GPs and gsum 
    np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) 

    # call the mask function to set this for interpolation training
    self.mask = self.gp_mask(self.nb)

    # Set up gp objects with fixed mean and standard deviation 
    kernel = self.gp_kernel(ls=3.0, sd=sd, center=0)  
    self.gp_interp = gm.ConjugateGaussianProcess(
        kernel=kernel, center=center, disp=0, df=dof, scale=sd, nugget=0) 

    # fit and predict using the interpolated GP (on number density now)
    if self.orders_mask is not None:
        self.gp_interp.fit(self.nB[self.mask], self.coeffs_trunc[self.mask])
    else:
        self.gp_interp.fit(self.nB[self.mask], self.coeffs[self.mask])
    pred, std = self.gp_interp.predict(self.nB, return_std=True)
    underlying_std = np.sqrt(self.gp_interp.cov_factor_)

    # print the kernel parameters for viewing outside the function
    print('Interpolation parameters: ', self.gp_interp.kernel_)
    print(self.gp_interp.cov_factor_)

    return pred, std, underlying_std

gp_kernel(ls=3.0, sd=0.5, center=0, nugget=1e-10)

The kernel that we will use both for interpolating the coefficients and for predicting the truncation error bands. This one is unfixed, so the value of the ls obtained here will be used to fix the second run when calling params attribute.

Parameters:

Name	Type	Description	Default
ls	float	The lengthscale guess for the kernel.	3.0
sd	float	The scale for the prior.	0.5
center	float	The center value for the prior.	0
nugget	(int, float)	The value of the nugget to send to the Cholesky decomposition.	1e-10

Returns:

Name	Type	Description
kernel	sklearn object	The kernel needed for the GPs.

Source code in src/neutron_rich_bmm/truncation_error_dens.py

def gp_kernel(self, ls=3.0, sd=0.5, center=0, nugget=1e-10):

    '''
    The kernel that we will use both for interpolating the 
    coefficients and for predicting the truncation error bands.
    This one is unfixed, so the value of the ls obtained here will 
    be used to fix the second run when calling params attribute.

    Parameters:
        ls (float): The lengthscale guess for the kernel.

        sd (float): The scale for the prior.

        center (float): The center value for the prior.

        nugget (int, float): The value of the nugget to send to the 
            Cholesky decomposition.

    Returns:
        kernel (sklearn object): The kernel needed for the GPs. 
    '''

    self.ls = ls    # starting guess; can get really close if we set 0.25 and fix it
    self.sd = sd    # makes a difference on the band of the regression curve for c_2 
    self.center = center
    self.nugget = nugget # nugget for the Cholesky, not the kernel 
    kernel = RBF(length_scale=self.ls)

    return kernel

gp_mask(nb)

The mask array needed to correctly separate our training and testing data.

Parameters:

Name	Type	Description	Default
nb	ndarray	The number density linspace needed.	required

Returns:

Type	Description
	self.mask (numpy.ndarray): The mask for use when interpolating or using the truncated GP.

Source code in src/neutron_rich_bmm/truncation_error_dens.py

def gp_mask(self, nb):

    '''
    The mask array needed to correctly separate our training
    and testing data.

    Parameters:
        nb (numpy.ndarray): The number density linspace needed.

    Returns:
        self.mask (numpy.ndarray): The mask for use when 
            interpolating or using the truncated GP.
    '''

    # mask for values above 40*n0 only (good)
    low_bound = next(i for i, val in enumerate(nb)
                              if val > 6.56)  # fm^-3 
    dens_mask = nb[low_bound:]

    # set the mask s.t. it picks the same training point number each time
    mask_num = len(dens_mask) // 2.5  
    mask_true = np.array([(i) % mask_num == 0 for i in range(len(dens_mask))]) 

    # concatenate with a mask over the other elements of nb before low_bound
    mask_false = np.full((1,len(nb[:low_bound])), fill_value=False)
    self.mask = np.concatenate((mask_false[0], mask_true))

    return self.mask 

uncertainties(data=None, expQ=None, yref=None, dof=3.0, sd=0.5, nugget=1e-10, excluded=None)

Calculation of the truncation error bands for the pQCD EOS, using the Gorda et al. (2023) formulation for the pressure. This function uses techniques from the gsum package.

Parameters:

Name	Type	Description	Default
data	ndarray	The data given in an array of arrays for each order-by-order result.	None
expQ	function	The functional form of the expansion parameter for gsum to use.	None
yref	function	The functional form of yref for gsum to use.	None
sd	float	The scale for the prior.	0.5
nugget	(int, float)	The nugget for the Cholesky decomposition.	1e-10
excluded	list	The orders we wish to exclude from training on in the coefficient arrays. Default is None.	None

Returns:

Name	Type	Description
data	ndarray	The data array, containing partials at each order.
		self.coeffs (numpy.ndarray): The values of the coefficents at nb.
std_trunc	ndarray	The arrays of truncation errors per each order.

Source code in src/neutron_rich_bmm/truncation_error_dens.py

def uncertainties(self, data=None, expQ=None, yref=None, dof=3.0, sd=0.5, nugget=1e-10, excluded=None):

    '''
    Calculation of the truncation error bands for the pQCD EOS, using 
    the Gorda et al. (2023) formulation for the pressure.
    This function uses techniques from the gsum package. 

    Parameters:
        data (numpy.ndarray): The data given in an array of arrays for 
            each order-by-order result.

        expQ (function): The functional form of the expansion parameter
            for gsum to use.

        yref (function): The functional form of yref for gsum to use.

        sd (float): The scale for the prior.

        nugget (int, float): The nugget for the Cholesky decomposition.

        excluded (list): The orders we wish to exclude from training
            on in the coefficient arrays. Default is None.

    Returns:
        data (numpy.ndarray): The data array, containing partials at each order.

        self.coeffs (numpy.ndarray): The values of the coefficents at nb.

        std_trunc (numpy.ndarray): The arrays of truncation errors per each order.
    '''

    # construct mask
    self.mask = self.gp_mask(self.nb)

    # get correct data shape
    if data is None:         
        data = self.data_all[:, :self.n_orders]
    else:
        data = data[:, :self.n_orders]

    # save nugget
    self.nugget = nugget

    # set up the truncation GP
    self.kernel = self.gp_kernel(ls=3.0, sd=sd, center=0, nugget=self.nugget) 
    self.trunc_gp = gm.TruncationGP(kernel=self.kernel, ref=yref, \
                        ratio=expQ, disp=0, df=dof, scale=sd, excluded=excluded, nugget=self.nugget)

    self.trunc_gp.fit(self.nB[self.mask], data[self.mask], orders=self.orders)

    std_trunc = np.zeros([len(self.nB), self.n_orders])
    cov_trunc = np.zeros([len(self.nB), len(self.nB), self.n_orders])
    for i, n in enumerate(self.orders):
        # Only get the uncertainty due to truncation (kind='trunc')
        _, std_trunc[:,n] = self.trunc_gp.predict(self.nB, order=n, return_std=True, \
                                                  kind='trunc', pred_noise=True)
        _, cov_trunc[:,:,n] = self.trunc_gp.predict(self.nB, order=n, return_std=False, \
                                                    return_cov=True, kind='trunc', pred_noise=True)

    # external access without altering return
    self.cov_trunc = cov_trunc

    # check the kernel hyperparameters
    print(self.trunc_gp.coeffs_process.kernel_)
    print(self.trunc_gp.coeffs_process.nugget)

    # save the GP object from this GP if possible
    self.trunc_gp_obj = self.trunc_gp

    return data, self.coeffs, std_trunc, cov_trunc