src package

Submodules

src.ACBICI.calibrator module

Created on Fri Feb 17 16:09:02 2023

@author: Ignacio Romero, Christina Schenk

class src.ACBICI.calibrator.ACBICImodel(*args, **kwargs)

Bases: ABC

addParameter(*, label=None, prior)

Appends a parameter to the model that will be calibrated.

Parameters:

• label (a string with the name of the parameter.)

• prior (a prior distribution for the parameter)

• Example

• heta$' (self.addParameter(label=r'$)

• prior=Uniform(a=1.0

• b=2.0))

Return type:

None

getInputDimension()

Returns the dimension of the input variables of the model.

Parameters:

None

Return type:

The dimension of the input variables

getNParam()

Returns the number of parameters in the model.

Parameters:

None

Return type:

The number of parameters.

getOutputDimension()

Returns the dimension of the output variables of the model = number of tasks.

Parameters:

None

Return type:

The dimension of the output variables

abstractmethod symbolicModel(x, t)

Evaluates the model f(x,t), the expression that we want to calibrate.

Parameters:

• x (an np array. On each row, a vector of observable inputs)

• t (an np array. On each row, a vector of parameters)

Return type:

An np array, value f(x,t) at each of the pairs (x_i, t_i)

class src.ACBICI.calibrator.Calibrator(*args, **kwargs)

Bases: ABC

PHEbounds()

Compute lower and upper bounds for the parameters to be calibrated

Parameters:

None

Returns:

• lower (array with lower bounds. Can be -np.inf)

• upper (array with upper boudns. Can be +np.in)

PHEguess()

Compute lower and upper bounds guesses for the parameters to be calibrated

Parameters:

None

Returns:

• lower (array with lower guess. Can be -np.inf)

• upper (array with upper guess. Can be +np.in)

addHyperparameter(name, prior)

Appends a hyperparameter to the calibration object.

Parameters:

• name ((string) the name of the hyperparameter)

• prior (a prior probability)

Return type:

None

abstractmethod calibrate(*, method='mcmc', nsteps=10000, burn=0.2, thin=1, nwalkers=16, nsynthetic=30)

This is the main function of the calibrator. It taks care of the whole calibration process.

Parameters:

• (optional) (nsteps)

• burn ((optional): the ration of burn in samples)

• thin (1 if thinning is to be performed.)

• nwalkers (number of walkers for parallel MCMC)

• nsynthetic (Number of synthetic points to be generated if a) – surrogate is built. If the calibrator does not u use a surrogate, the option is ignored.

Return type:

None

gaussianLogLikelihood(x, mu, sigma, sparse=False)

Calculate the log-likelihood of a Gaussian variable at a value, given its mean and covariance.

Parameters:

• x (the value where the log-likelihood is to be evaluated)

• mu (the mean vector. dim(mu) = dim(x))

• sigma (the covariance matrix. dim(sigma) = dim(x) x dim(x))

• sparse (flag for using sparse Cholesky, requires scikit-sparse package,)

• Cholesky (default=False uses dense)

Returns:

• log-likelihood or the log-likelihood minus a constant that does

• not depend on x.

generateSyntheticData(*, npoints=30)

Generate synthetic (random) input/parameters pairs, run the model and store the synthetic (x, theta, y) triplets that will be used to feed the surrogate.

Parameters:

npoints (number of synthetic points to be generated.)

Return type:

None

genericCalibration(*, method='mcmc', nsteps=10000, burn=0.2, thin=1, nwalkers=16, fposterior='samples', flogprob='logprob', flogprior='logprior')

This function takes care of the calibration process, once all the parameters and required functions have been taken care of. It should not be called directly by the user but from a derived calibrator.

Parameters:

• method ("mcmc": Markov Chain Monte Carlo (default)) – “vbmc”: Variational Bayesian Monte Carlo using pyvbmc “vbmcsimple”: a simplified Variational Bayes Monte Carlo solver

• nsteps (number of MCMC steps)

• burn (ratio of burn in samples in MCMC)

• thin (1 if thinning is to be performed in MCMC)

• nwalkers (number of walkers for parallel MCMC)

• fposterior (filename where the posteriors will be written)

• flogprob (filename where the log-pdf will be written)

• flogprior (filename where the priors will be written)

Return type:

None

getAllLabels()

Retrieves the labels (strings) of all the variables that need to be calibrated: the model parameters, the calibrator hyperparameters (if any) and the experimental error (if unknown).

Parameters:

--

Return type:

Array of labels.

getExperimentalSTD()

Parameters:

--

Return type:

The value of the experimental error.

getNHyper()

Returns the number of hyperparameters in the GP (if any).

Parameters:

--

Return type:

The number hyperparameters in the calibration (>=0)

getNPHE()

Returns the number of variables that need to be calibrated in the problem. This includes the parameters in the model, the hyperparameters in the GP (if any) and the experimental error (if unknown).

Parameters:

--

Returns:

• The number of parameters + hyperparameters (if any) + experimental

• error (if any) in the calibration.

getNParam()

Returns the number of parameters in the model that is being calibrated.

Parameters:

--

Return type:

The number of parameters.

abstractmethod logLikelihood_vect()

Vectorized function that calculates the log-likelihood

Parameters:

None

Return type:

None

logPosterior(PHE)

Log posterior computation

Parameters:

PHE (an np array with the parameters, hyperparameters (if any), expError (if any))

Return type:

logPosterior / -inf if the Hyperpar-Par-Exp have invalid values

logPriors()

Logarithm of priors samples for the PHE.

Parameters:

None

Return type:

Array of log-priors of params + hyperparams + exp

mcmcCalibration(*, nsteps=10000, burn=0.2, thin=1, nwalkers=16, fposterior='samples', flogprob='logprob', flogprior='logprior')

Markov chain Monte Carlo calibration.

Parameters:

• nsteps (number of MCMC steps)

• burn (ratio of burn in samples in MCMC)

• thin (1 if thinning is to be performed in MCMC)

• nwalkers (number of walkers for parallel MCMC)

• fposterior (filename where the posteriors will be written)

• flogprob (filename where the log-pdf will be written)

• flogprior (filename where the priors will be written)

Return type:

None

abstractmethod plot(*, compare=True, discrepancy=True, prior=True, trace=True, corner=True, pearson=True, prerror=True, onscreen=False, dumpfiles=True)

Abstract function for the plot generation.

Parameters:

• Compare (generate data vs. predictions plot if True)

• Discrepancy (generate a plot wit the discrepancy function if True)

• prior (generate figure with prior probabilities for parameters if True)

• trace (generate MCMC trace plots if True)

• corner (generate corner plots for posterior if True)

• pearson (generate correlation plot)

• prerror (generate prediction error plot)

• dumpfiles (write files to disk if True)

Return type:

None

plotCorner(*, onscreen=True, dumpfiles=True)

Creates a plot with the pairwise and single variable marginalized probability distributions of parameters and hyperparameters. Show the MAP and mean.

Parameters:

• onscreen (if True, plot the figure on the screen)

• dumpfiles (if True, write the figure on a file)

Return type:

None

plotCorrelations(*, onscreen=True, dumpfiles=True)

Writes a short report on the correlations between variables in the log file.

Parameters:

• onscreen (if True, plot the figure on the screen)

• dumpfiles (if True, write the figure on a file)

Return type:

None

plotPredictionErrors(*, onscreen=True, dumpfiles=True)

Generate a figure illustrating the errors made by the calibrated model.

Parameters:

• onscreen (if True, plot the figure on the screen)

• dumpfiles (if True, write the figure on a file)

Return type:

None

plotPriorsPosteriors(*, onscreen=True, dumpfiles=True)

Generate a figure illustrating the prior probability distributions for the parameters and hyperparameters.

Parameters:

• onscreen (if True, plot the figure on the screen)

• dumpfiles (if True, write the figure on a file)

Return type:

None

plotTrace(*, onscreen=True, dumpfiles=True)

Creates a plot with the MCMC sample trace. It thins the number of points to NPOINTS, so as to be able to distinguish the trace.

Parameters:

• onscreen (if True, plot the figure on the screen)

• dumpfiles (if True, write the figure on a file)

Return type:

None

predictionError()

Once a model has been calibrated, compute an array of errors made comparing the experimental values with the model predictions.

Parameters:

None

Return type:

An array of errors of shape nexp x ydim , or None if it can not be calculated

printInfo()

Prints a summary of the calibration process in the log file

Parameters:

None

Return type:

None

printReport()

Writes a file with the statistics for the posterior distribution.

Parameters:

None

Return type:

None

printlog(*args, **kwargs)

This is a wrapper that allows use the command ‘print’ directly dumping on the log file of the calibration.

Parameters:

• *args (this is the standard way of)

• **kwargs (this is the standard way of)

• function. (declaring the arguments for the print)

Return type:

None

randomize()

Compute random value of vars to calibrate according to prior distributions.

Parameters:

Nonen

Return type:

Array with random, valid values for the PHE.

replaceHyperparameterPrior(*, label, prior)

Replace the prior distribution for a hyperparameter in the calibration. All calibrations, when created, assign a default prior distribution to all their hyperparameters. The user might want to use this function to replace the default with a more appropriate one. Careful: the label provided to the function has to match one of the existing labels for the calibrator hyperparameters. Example: theCalibrator.replaceHyperparameterPrior(label=r’\(eta_x\)’,prior=Cauchy(mu=1.0, sigma=0.1))

Parameters:

• label (string with the name of the hyperparameter to be replaced.)

• prior (prior probability distribution (of ACBICI type))

Return type:

None

setDefaultExperimentalSTDPrior()

Set a default prior distribution for the standard deviation of the experimental error, itself being a normal variable. A desirable mean and std are first calculated in terms of the available experimental data. Then, the prior parameters are set so that the selected distribution has this mean and std. Note that the prior support must be the positive numbers

Parameters:

None

Return type:

None

setExperimentalSTDPrior(prior)

All calibrators assume that the experimental measurements have a certain error which is Gaussian and has zero mean. With this command, the user sets the standard deviation of the Gaussian to be a random variable itself, and provides a prior probability distribution for it.

Parameters:

prior (prior probability distribution (of ACBICI type))

Return type:

None

setExperimentalSTDValue(val)

All calibrators assume that the experimental measurements have a certain error which is Gaussian and has zero mean. With this command, the user fixes the standard deviation of the Gaussian.

Parameters:

val (standard deviation of the distribution of the experimental error.)

Return type:

None

storeExperimentalData(data)

This function transfers the experimental data, a numpy array, to the calibrator for its later user. The data must have as many rows as desired but xdim+1 columns, where xdim is the dimension of the input space. Saves xExp in (nsamples, xdim) and yExp in (nsamples,) format.

Parameters:

data (numpy array with the (input, output) data from experiments.)

Return type:

None

storePHE(phe)

Given an array that contains the model parameters, the surrogate hyperparameters (if any) and the uknown experimental error (if any), store each data in its position.

storeSyntheticData(data)

This function transfers the synthetic data, a numpy array, to the calibrator for its later user. The data must have as many rows as desired but xdim+pdim+1 columns, where xdim is the dimension of the input space, pdim is the number of parameters in the model. Saves xSyn in (nsamples, xdim) and ySyn in (nsamples,) format.

Parameters:

data (numpy array with the (x_i, theta_i, output_i) from model runs.)

Return type:

None

class src.ACBICI.calibrator.KOHCalibrator(*args, **kwargs)

Bases: Calibrator

KOH model class. Features of the calibrator: - surrogate gaussian process for the model. - gaussian process for the discrepancy (optional) - experimental error is gaussian, with known variance

calibrate(*, method='mcmc', nsteps=10000, burn=0.2, thin=1, nwalkers=16, nsynthetic=30)

This is the main function of the calibrator. It taks care of the whole calibration process.

Parameters:

• (optional) (nsteps)

• burn ((optional): the ration of burn in samples)

• thin (1 if thinning is to be performed.)

• nwalkers (number of walkers for parallel MCMC)

• nsynthetic (Number of synthetic points to be generated if a) – surrogate is built. If the calibrator does not u use a surrogate, the option is ignored.

Return type:

None

covd(x1, x2)

covt(t1, t2)

covx(x1, x2)

getOutputDimension()

Returns the dimension of the output variables of the model = number of tasks.

Parameters:

None

Return type:

The dimension

logLikelihood_vect()

Log-likelihood function for KOH ansatz with GP surrogates for the model and discrepancy.

This is the log-likelihood of z being sampled from a KOH ansatz for which we have x datapoints, and hyperparameters beta, lam, …

Parameters:

None

Returns:

• log-likelihood of the expensive model with the value of the parameters

• and hyperparameters

plot(*, compare=True, discrepancy=True, prior=True, trace=True, corner=True, pearson=True, prerror=True, onscreen=False, dumpfiles=True)

Abstract function for the plot generation.

Parameters:

• Compare (generate data vs. predictions plot if True)

• Discrepancy (generate a plot wit the discrepancy function if True)

• prior (generate figure with prior probabilities for parameters if True)

• trace (generate MCMC trace plots if True)

• corner (generate corner plots for posterior if True)

• pearson (generate correlation plot)

• prerror (generate prediction error plot)

• dumpfiles (write files to disk if True)

Return type:

None

plotCompare(*, onscreen=False, dumpfiles=True)

plotCompare_multix(*, onscreen=False, dumpfiles=True)

plotDiscrepancy(*, onscreen=False, dumpfiles=True)

plotDiscrepancy_multix(*, onscreen=False, dumpfiles=True)

remove_extrarows(theta, num_tasks)

Reverse the augmentation process by removing the extra rows and restoring the original theta.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

remove_task_index(augmented_xExp, num_tasks)

Reverse the augmentation process by removing the task index column and restoring original xExp.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

setDefaultHyperparametersPriors()

When the user does not provide a prior distribution for the hyperparameters, ACBICI will provide the ones in this function. For the lengthscale beta, we propose to use the average pairwise distance of the data in the numerator; for the variance lambda, we use the variance of the output.

Parameters:

None

Return type:

None

src.ACBICI.calibrator.average_pairwise_distance(data)

Given an two dimensional array of size n x m, this function assume that it corresponds to n vectors of data, each of them of length m. Then, it finds the average Euclidean distance among all data

Parameters:

data (a numpy array of dim n x m, containing n vectors of size m.)

Return type:

The average mean Euclidean distance between all pairs of data.

class src.ACBICI.calibrator.classicalCalibrator(*args, **kwargs)

Bases: Calibrator

A class that takes care of a classical Bayesian calibration process: given an analytical expression (or function) that depends on a controlable input and a collection of parameters, finds the best values for the paramters in the sense that they explain the experimental results best and are consistent with our prior knowledge. The experimental data must have a known experimental error assumed to be Gaussian.

calibrate(*, method='mcmc', nsteps=10000, burn=0.2, thin=1, nwalkers=16, nsynthetic=30)

This is the main function of the calibrator. It taks care of the whole calibration process.

Parameters:

• (optional) (nsteps)

• burn ((optional): the ration of burn in samples)

• thin (1 if thinning is to be performed.)

• nwalkers (number of walkers for parallel MCMC)

• nsynthetic (Number of synthetic points to be generated if a) – surrogate is built. If the calibrator does not u use a surrogate, the option is ignored.

Return type:

None

logLikelihood_vect()

Log-likelihood function for known model and normal experimental error.

This is the log-likelihood of z being sampled from a KOH ansatz for which we have x datapoints, and hyperparameters beta, lam, …

Parameters:

None

Returns:

• log-likelihood of the KOH model with the value of the parameters

• and hyperparameters

plot(*, compare=True, discrepancy=True, prior=True, trace=True, corner=True, pearson=True, prerror=True, onscreen=False, dumpfiles=True)

Abstract function for the plot generation.

Parameters:

• Compare (generate data vs. predictions plot if True)

• Discrepancy (generate a plot wit the discrepancy function if True)

• prior (generate figure with prior probabilities for parameters if True)

• trace (generate MCMC trace plots if True)

• corner (generate corner plots for posterior if True)

• pearson (generate correlation plot)

• prerror (generate prediction error plot)

• dumpfiles (write files to disk if True)

Return type:

None

plotCompare(*, onscreen=True, dumpfiles=True)

Creates a plot that compares the provided data and the predictions of the calibrated model.

Parameters:

• onscreen (if True, plot the figure on the screen)

• dumpfiles (if True, write the figure on a file)

Return type:

None

plotCompare_multix(*, onscreen=True, dumpfiles=True)

Creates a plot for each x dimension that compares the provided data and the predictions of the calibrated model.

Parameters:

• onscreen (if True, plot the figure on the screen)

• dumpfiles (if True, write the figure on a file)

Return type:

None

remove_extrarows(theta, num_tasks)

Reverse the augmentation process by removing the extra rows and restoring the original theta.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

remove_task_index(augmented_xExp, num_tasks)

Reverse the augmentation process by removing the task index column and restoring original xExp.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

setDefaultHyperparametersPriors()

When the user does not provide a prior distribution for the hyperparameters, ACBICI will provide the ones in this function. For the classical calibrator, only the std of the experimental error needs to be initialized.

Parameters:

None

Return type:

None

storeSamples(xExperimental, yExperimental)

Function to store experimental data

Parameters:

• xExperimental (x values for experimental datapoints)

• yExperimental (y values for experimental datapoints)

Return type:

None

transformed_model(xExp, mtheta)

class src.ACBICI.calibrator.discrepancyCalibrator(*args, **kwargs)

Bases: Calibrator

Features of the calibrator:

• the model is known exactly and it is easy to evaluate.

• gaussian process for the discrepancy (optional)

• experimental error is gaussian. Its variance may or may not be known.

calibrate(*, method='mcmc', nsteps=10000, burn=0.2, thin=1, nwalkers=16, nsynthetic=30, lambda_x_std=None)

This is the main function of the calibrator. It taks care of the whole calibration process.

Parameters:

• (optional) (nsteps)

• burn ((optional): the ration of burn in samples)

• thin (1 if thinning is to be performed.)

• nwalkers (number of walkers for parallel MCMC)

• nsynthetic (Number of synthetic points to be generated if a) – surrogate is built. If the calibrator does not u use a surrogate, the option is ignored.

Return type:

None

covd(x1, x2)

Covariance discrepancy kernel for pairs of input vectors.

Parameters:

• x1 (first input vector. Dimensions xdim+pdim)

• x2 (second input vector. Dimensions xdim+pdim)

Return type:

Covariance matrix for discrepancy

getOutputDimension()

Returns the dimension of the output variables of the model = number of tasks.

Parameters:

None

Return type:

The dimension

logLikelihood_vect()

Vectorized version Log-likelihood function for KOH ansatz with GP for the discrepancy error.

This is the log-likelihood of z being sampled from a KOH ansatz for which we have x datapoints, and hyperparameters beta, lam, …

Parameters:

None

Return type:

log-likelihood of the KOH model with the value of the parameters and hyperparameters

plot(*, compare=True, discrepancy=True, prior=True, trace=True, corner=True, pearson=True, prerror=True, onscreen=False, dumpfiles=True)

Abstract function for the plot generation.

Parameters:

• Compare (generate data vs. predictions plot if True)

• Discrepancy (generate a plot wit the discrepancy function if True)

• prior (generate figure with prior probabilities for parameters if True)

• trace (generate MCMC trace plots if True)

• corner (generate corner plots for posterior if True)

• pearson (generate correlation plot)

• prerror (generate prediction error plot)

• dumpfiles (write files to disk if True)

Return type:

None

plotCompare(*, onscreen=False, dumpfiles=True)

plotCompare_multix(*, onscreen=False, dumpfiles=True)

plotDiscrepancy(*, onscreen=False, dumpfiles=True)

plotDiscrepancy_multix(*, onscreen=False, dumpfiles=True)

remove_extrarows(theta, num_tasks)

Reverse the augmentation process by removing the extra rows and restoring the original theta.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

remove_task_index(augmented_xExp, num_tasks)

Reverse the augmentation process by removing the task index column and restoring original xExp.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

setDefaultHyperparametersPriors()

When the user does not provide a prior distribution for the hyperparameters, ACBICI will provide the ones in this function. For the lengthscale beta, we propose to use the average pairwise distance of the data in the numerator; for the variance lambda, we use the variance of the output.

Parameters:

None

Return type:

None

transformed_model(xExp, mtheta)

class src.ACBICI.calibrator.expensiveCalibrator(*args, **kwargs)

Bases: Calibrator

Features of the calibrator:

• surrogate gaussian process for the model.

• experimental error is gaussian, with known variance

calibrate(*, method='mcmc', nsteps=10000, burn=0.2, thin=1, nwalkers=16, nsynthetic=30)

This is the main function of the calibrator. It taks care of the whole calibration process.

Parameters:

• (optional) (nsteps)

• burn ((optional): the ration of burn in samples)

• thin (1 if thinning is to be performed.)

• nwalkers (number of walkers for parallel MCMC)

• nsynthetic (Number of synthetic points to be generated if a) – surrogate is built. If the calibrator does not u use a surrogate, the option is ignored.

Return type:

None

covt(t1, t2)

Computes the vallue of the covariance matrix for two arrays of synthetic data. The two arrays could be the same.

Parameters:

• t1 (an numpy array of input points. Each row must be xdim+pdim size.)

• t2 (an numpy array of input points. Each row must be xdim+pdim size.)

Return type:

Covariance matrix.

covx(x1, x2)

Computes the value of the covariance matrix for two arrays of input_data + parameters. The two arrays x1, x2 could be the same one.

Parameters:

• x1 (an numpy array of input points. Each row must be xdim+pdim size.)

• x2 (an numpy array of input points. Each row must be xdim+pdim size.)

Return type:

Covariance matrix.

getOutputDimension()

Returns the dimension of the output variables of the model = number of tasks.

Parameters:

None

Return type:

The dimension

logLikelihood_vect()

Log-likelihood function for KOH ansatz with GP surrogate for the model.

This is the log-likelihood of z being sampled from a KOH ansatz for which we have x datapoints, and hyperparameters beta, lam, …

Parameters:

None

Returns:

• log-likelihood of the expensive model with the value of the parameters

• and hyperparameters

plot(*, compare=True, discrepancy=True, prior=True, trace=True, corner=True, pearson=True, prerror=True, onscreen=False, dumpfiles=True)

Abstract function for the plot generation.

Parameters:

• Compare (generate data vs. predictions plot if True)

• Discrepancy (generate a plot wit the discrepancy function if True)

• prior (generate figure with prior probabilities for parameters if True)

• trace (generate MCMC trace plots if True)

• corner (generate corner plots for posterior if True)

• pearson (generate correlation plot)

• prerror (generate prediction error plot)

• dumpfiles (write files to disk if True)

Return type:

None

plotCompare(*, onscreen=False, dumpfiles=True)

plotCompare_multix(*, onscreen=False, dumpfiles=True)

remove_extrarows(theta, num_tasks)

Reverse the augmentation process by removing the extra rows and restoring the original theta.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

remove_task_index(augmented_xExp, num_tasks)

Reverse the augmentation process by removing the task index column and restoring original xExp.

Parameters:

• augmented_xExp (numpy array of shape (n * num_tasks, xdim + 1)) – The augmented data with repeated xExp rows and an additional task index column.

• num_tasks (int) – The number of task repetitions.

Returns:

original_xExp – The restored xExp without task indices.

Return type:

numpy array of shape (n, xdim)

setDefaultHyperparametersPriors()

When the user does not provide a prior distribution for the hyperparameters, ACBICI will provide the ones in this function. For the lengthscale beta, we propose to use the average pairwise distance of the data in the numerator; for the variance lambda, we use the variance of the output.

Parameters:

None

Return type:

None

src.ACBICI.calibrator.findApproximateMAP(data)

Obtains the approximate value of the MAP estimate in each column of x.

Parameters:

• data (an array of samples. Each row is a vector or samples of the)

• the (various variables. Each column corresponds to many samples of)

• variable. (same)

Returns:

MAP

Return type:

the map of each column

src.ACBICI.calibrator.is_function_effectively_empty(func)

src.ACBICI.calibrator.log_multivariate_normal_pdf(x, mu, cov)

Calculates the log likelihood with dense Cholesky decomposition

Parameters:

• x (the value where the log-likelihood is to be evaluated)

• mu (the mean vector. dim(mu) = dim(x))

• cov (the covariance matrix. dim(cov) = dim(x) x dim(x))

Returns:

• log-likelihood or the log-likelihood minus a constant that does

• not depend on x.

src.ACBICI.calibrator.log_multivariate_normal_pdf_sparse(x, mu, cov)

Calculates the log likelihood with sparse Cholesky decomposition

Parameters:

• x (the value where the log-likelihood is to be evaluated)

• mu (the mean vector. dim(mu) = dim(x))

• cov (the covariance matrix. dim(cov) = dim(x) x dim(x))

Returns:

• log-likelihood or the log-likelihood minus a constant that does

• not depend on x.

src.ACBICI.calibrator.printStatistics(names, samples, file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Prints a nice description of the main statistics of the sample.

Parameters:

• names (an array of strings with the names of the variables in samples)

• samples (an array of samples. Each row has the same length as 'names')

Return type:

None

src.ACBICI.gprocess module

Created on Fri Feb 17 16:09:02 2023

@author: Ignacio Romero, Christina Schenk

class src.ACBICI.gprocess.GaussianProcess(X, z, mu, cov, noisesigma=0.0)

Bases: object

A Gaussian process.

covarianceMatrix()

function for covarianceMatrix

Parameters:

self (constructor)

Returns:

Sigma

Return type:

covariance matrix

eval(Xstar)

eval function to return mean and std of the GP at points Xstar, when the parameters theta and the hyperparameters have been set before.

Parameters:

• self (constructor)

• Xstar (array of input points of dimension number of points where we want to evaluate it x 1)

Returns:

• zstar (mean)

• np.sqrt(diag) (standard deviation)

info()

info function

Parameters:

self (constructor)

Return type:

print number of data points

static test()

Use this function to test the vanilla Gaussian process. A random vector of datapoints is created, and a GP is built from them, using random hyperparameters. Finally, plots are generated to see that everything works ok.

Return type:

Plot for validating plain GP with noise

src.ACBICI.kernels module

Created on Tue Jan 28

@author: Christina Schenk, Ignacio Romero

class src.ACBICI.kernels.Expo(*args, **kwargs)

Bases: kernel

Exponential kernel. Exponential growth. Rate increases as input increases.

evalkernel(dist, lamb, beta)

Evaluates the kernel for a distance ‘dist’, given the two hyperparameters lambda and beta. All the kernels must be of the form k = lambda * f(distance/beta)

Parameters:

• dist (distance)

• lamb (hyperparameter lambda)

• beta (hyperparameter beta)

Return type:

kernel evaluation

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print information about the kernel on a file

Parameters:

file (filename)

Return type:

None

class src.ACBICI.kernels.Matern32(*args, **kwargs)

Bases: kernel

Once differentiable Matern kernel.

evalkernel(dist, lamb, beta)

Evaluates the kernel for a distance ‘dist’, given the two hyperparameters lambda and beta. All the kernels must be of the form k = lambda * f(distance/beta)

Parameters:

• dist (distance)

• lamb (hyperparameter lambda)

• beta (hyperparameter beta)

Return type:

kernel evaluation

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print information about the kernel on a file

Parameters:

file (filename)

Return type:

None

class src.ACBICI.kernels.Matern52(*args, **kwargs)

Bases: kernel

Twice differentiable Matern kernel.

evalkernel(dist, lamb, beta)

Evaluates the kernel for a distance ‘dist’, given the two hyperparameters lambda and beta. All the kernels must be of the form k = lambda * f(distance/beta)

Parameters:

• dist (distance)

• lamb (hyperparameter lambda)

• beta (hyperparameter beta)

Return type:

kernel evaluation

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print information about the kernel on a file

Parameters:

file (filename)

Return type:

None

class src.ACBICI.kernels.MultiTask(*args, **kwargs)

Bases: kernel

Multi Task kernel assuming that tasks independent, just comparing for same task dimension 0 for different tasks or matern3/2 kernel if same task

evalkernel(dist, lamb, beta, x1=None, x2=None)

Evaluates the kernel for a distance ‘dist’, given the two hyperparameters lambda and beta. All the kernels must be of the form k = lambda * f(distance/beta)

Parameters:

• dist (distance)

• lamb (hyperparameter lambda)

• beta (hyperparameter beta)

Return type:

kernel evaluation

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print information about the kernel on a file

Parameters:

file (filename)

Return type:

None

class src.ACBICI.kernels.RatQuad(*args, **kwargs)

Bases: kernel

Rational quadratic kernel with alpha=1.

evalkernel(dist, lamb, beta)

Evaluates the kernel for a distance ‘dist’, given the two hyperparameters lambda and beta. All the kernels must be of the form k = lambda * f(distance/beta)

Parameters:

• dist (distance)

• lamb (hyperparameter lambda)

• beta (hyperparameter beta)

Return type:

kernel evaluation

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print information about the kernel on a file

Parameters:

file (filename)

Return type:

None

class src.ACBICI.kernels.SqExpo(*args, **kwargs)

Bases: kernel

Squared exponential kernel (RBF or Gaussian Kernel).

evalkernel(dist, lamb, beta, x1, x2)

Evaluates the kernel for a distance ‘dist’, given the two hyperparameters lambda and beta. All the kernels must be of the form k = lambda * f(distance/beta)

Parameters:

• dist (distance)

• lamb (hyperparameter lambda)

• beta (hyperparameter beta)

Return type:

kernel evaluation

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print information about the kernel on a file

Parameters:

file (filename)

Return type:

None

class src.ACBICI.kernels.kernel(*args, **kwargs)

Bases: ABC

abstractmethod evalkernel(dist, lamb, beta, x1=None, x2=None)

Evaluates the kernel for a distance ‘dist’, given the two hyperparameters lambda and beta. All the kernels must be of the form k = lambda * f(distance/beta)

Parameters:

• dist (distance)

• lamb (hyperparameter lambda)

• beta (hyperparameter beta)

Return type:

kernel evaluation

abstractmethod print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Print information about the kernel on a file

Parameters:

file (filename)

Return type:

None

src.ACBICI.randomprocess module

Created on Fri Feb 17 16:09:02 2023

@author: Ignacio Romero, Christina Schenk

class src.ACBICI.randomprocess.KOHGaussianProcess(xExperimental, yExperimental, xSynthetic, tSynthetic, ySynthetic, calibrator)

Bases: randomProcess

Gaussian processes that have the structure of x ~ (symbolic equation or GP of metamodel) + (GP of discrepancy) + (experimental error)

Reference: M.C. Kennedy and A. O’Hagan (2001): Bayesian calibration of computer models, J.R. Statist. Soc. B (2001), 63, Part 3 , pp. 425-464.

covarianceMatrix_vect()

Vectorized code Evaluate the covariance matrix when we fix the value of the experimental data, synthetic data, parameters theta and hyperparameters

Parameters:

---

Returns:

Sigma

Return type:

covariance matrix

covarianceMatrix_vect_withdisc()

Evaluate the covariance matrix with discrepancy when we fix the value of the experimental data, synthetic data, parameters theta and hyperparameters

Parameters:

---

Returns:

Sigma

Return type:

covariance matrix

eval_vect(Xstar)

Vectorized code Evaluate mean and std of the GP at points Xstar. The parameters theta and the hyperparameters have been set before.

Parameters:

Xstar (array of input points where we want to evaluate the GP) – The 2nd dimension should be equal to the dimension of the input variables x

Returns:

• ystar (mean at every point of Xstar)

• np.sqrt(diag) (Sqrt of the diagonal of the covariance matrix)

eval_vect_withdisc(Xstar)

Vectorized code Evaluate mean and std of the GP including discrepancy at points Xstar. The parameters theta and the hyperparameters have been set before.

Parameters:

Xstar (array of input points where we want to evaluate the GP) – The 2nd dimension should be equal to the dimension of the input variables x

Returns:

• ystar (mean at every point of Xstar)

• np.sqrt(diag) (Sqrt of the diagonal of the covariance matrix)

print()

Print information about the KOH GP.

Parameters:

none

Return type:

none

static test(self)

Use this function to test the KOH model.

Return type:

Plot for testing KOH class, plotting prediction vs. experimental values

class src.ACBICI.randomprocess.expensiveGaussianProcess(xExperimental, yExperimental, xSynthetic, tSynthetic, ySynthetic, calibrator)

Bases: randomProcess

Gaussian processes that have the structure of x ~ (GP of metamodel) + (experimental error)

Reference: M.C. Kennedy and A. O’Hagan (2001): Bayesian calibration of computer models, J.R. Statist. Soc. B (2001), 63, Part 3 , pp. 425-464.

covarianceMatrix_vect()

Vectorized code Evaluate the covariance matrix with discrepancy when we fix the value of the experimental data, synthetic data, parameters theta and hyperparameters

Parameters:

---

Returns:

Sigma

Return type:

covariance matrix

eval_vect(Xstar)

Vectorized code Evaluate mean and std of the GP at points Xstar. The parameters theta and the hyperparameters have been set before.

Parameters:

Xstar (array of input points where we want to evaluate the GP) – The 2nd dimension should be equal to the dimension of the input variables x

Returns:

• ystar (mean at every point of Xstar)

• np.sqrt(diag) (Sqrt of the diagonal of the covariance matrix)

print()

Print information about the GP for expensive model.

Parameters:

none

Return type:

none

class src.ACBICI.randomprocess.inexpKOHGaussianProcess(xExperimental, yExperimental, calibrator)

Bases: object

Class for Kennedy O’Hagan calibration that involves an inexpensive model.

covarianceMatrix_vect()

Vectorized version Evaluate the covariance matrix for the inexpensive KOH model.

Experimental data, parameters theta and hyperparameters are fixed. The matrix consists of the sum of the experimental (Gaussian) error and the discrepancy, if the model has the latter.

Parameters:

--

Returns:

Sigma

Return type:

covariance matrix

covarianceMatrix_vect_withdisc()

Vectorized version Evaluate the covariance matrix with discrepancy for the inexpensive KOH model.

Experimental data, parameters theta and hyperparameters are fixed. The matrix consists of the sum of the experimental (Gaussian) error and the discrepancy, if the model has the latter.

Parameters:

--

Returns:

Sigma

Return type:

covariance matrix

eval_vect(Xstar)

Vectorized code Evaluate mean and std of the GP at points Xstar.

Parameters:

Xstar (array of input points (npoints x dim))

Returns:

• ystar (mean vector)

• std (standard deviation estimate for every evaluation point)

eval_vect_withdisc(Xstar)

Vectorized code Evaluate mean and std of the GP including discrepancy at points Xstar.

Parameters:

Xstar (array of input points)

Returns:

• ystar (mean)

• std (standard deviation estimate for every evaluation point)

print()

Print information about the random process

Parameters:

none

Return type:

none

class src.ACBICI.randomprocess.randomProcess

Bases: ABC

Abstract class for all the random processes employed for calibration in ACBICI.

abstractmethod covarianceMatrix_vect()

Calculate the covariance matrix.

abstractmethod eval_vect(Xstar)

Evaluate the random process at inputs Xstar

Paramters

Xstar: a numpy array of input values. Each row is an input

rtype:

A numpy vector of evaluations.

abstractmethod print()

src.ACBICI.priors module

Created on Nov 25 2024 Class for probability priors.

@author: Christina Schenk, Ignacio Romero

class src.ACBICI.priors.Cauchy(*args, **kwargs)

Bases: prior

Prior distribution of the class Cauchy with parameters mu and sigma

The PDF is p(x) = 2/(pi sigma) * 1/(1+ (x-mu)^2/sigma^2) for any x The distribution is identical to the half-Cauchy, the difference being that the random variable might be any real number.

See: https://distribution-explorer.github.io/continuous/cauchy.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.cauchy.html#scipy.stats.cauchy

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

class src.ACBICI.priors.Gamma(*args, **kwargs)

Bases: prior

Prior distribution of the class Gamma with shape parameter alpha and rate parameter beta.

The PDF is p(x) = 1/Gamma(alpha) * beta^alpha * x^(alpha-1) * exp[-beta x] with mean alpha/beta and variance alpha/beta**2

See: https://distribution-explorer.github.io/continuous/gamma.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html#scipy.stats.gamma

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

class src.ACBICI.priors.HalfCauchy(*args, **kwargs)

Bases: prior

Prior distribution of the class half-Cauchy with parameters mu and sigma

The PDF is p(x) = 2/(pi sigma) * 1/(1+ (x-mu)^2/sigma^2) for x >= mu

See https://distribution-explorer.github.io/continuous/halfcauchy.html

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

class src.ACBICI.priors.HalfNormal(*args, **kwargs)

Bases: prior

Prior distribution of the class half-normal with parameters mu and sigma

See https://distribution-explorer.github.io/continuous/halfnormal.html

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

class src.ACBICI.priors.Normal(*args, **kwargs)

Bases: prior

Prior distribution of the class norm with parameters mu and sigma

The PDF is p(x) = 1/sqrt(2pi sigma^2) * exp[-0.5 (x-mu)^2/sigma^2]

See: https://distribution-explorer.github.io/continuous/norm.html https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html#scipy.stats.norm

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

class src.ACBICI.priors.Uniform(*args, **kwargs)

Bases: prior

Uniform prior distribution in the interval [a,b]. Probability density function is p(x) = 1/(b-a) for a <= x <= b.

See https://distribution-explorer.github.io/continuous/uniform.html

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

class src.ACBICI.priors.Weibull(*args, **kwargs)

Bases: prior

Weibull prior distribution. Probability density function is p(x) = 0 x<= 0 p(x) = (k/lambda) * (x/lambda)^(k-1) * exp(-(x/lambda)^k)

See https://distribution-explorer.github.io/continuous/weibull.html

print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

class src.ACBICI.priors.prior(*args, **kwargs)

Bases: ABC

bounds()

Return lower and upper bound values of the distribution. Either can be -np.inf or +np.inf

Parameters:

None

Returns:

• lower (lower bound)

• upper (upper bound)

guessInterval()

Return a fair interval for the expected value of the variable

Parameters:

None

Returns:

• lower (lower bound)

• upper (upper bound)

logProbability(x)

Log of the probability density function.

Parameters:

x (scalar at which the log-pdf has to be evaluated.)

Return type:

logpdf(x)

pdf(x)

Calculate the PDF at the values x.

Parameters:

x (scalar at which the pdf has to be evaluated.)

Return type:

pdf(x)

ppf(percentiles)

Calculate the percentiles point function for input(s)

Parameters:

percentiles (array of scalars of the percentiles that need to be calculated.) – values should be in [0,1]

Return type:

array with the values of the random variable where the percentiles are attained.

abstractmethod print(file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>)

Write on file information about the probability distribution.

randomSample()

Generate a single sample distributed according to the probability distribution of the prior.

src.ACBICI.vbmc module

class src.ACBICI.vbmc.DiagonalVBmc(theCalibrator, n_mc=30)

Bases: object

elbo(params)

Computes the negative of ELBO

Parameters:

params (a numpy array with the initial value of the parameters) – to be calibrated. These include the actual model parameters, the hyperparameters of the GP (if any), and the variance of the experimental error (if needed)

Return type:

Negative ELBO (evidence lower bound)

fit()

sample_posterior(n_samples=2000)

class src.ACBICI.vbmc.VBmc(theCalibrator)

Bases: object

A variational Bayes calibrator based on the package pyvbmc (https://github.com/acerbilab/pyvbmc)

fit()

This runs variational Bayes with pybvmc. The package does not allow to dump results on a file, but not on the screen, so we have to deactivate printing to the terminal.

sample_posterior(n_samples=2000)

Module contents