Python API

Calibration toolbox

This module contains the Bayesian calibration class.

class BayesianCalibration(system, calibration, num_iter=1, curr_iter=0, save_fig=-1, callback=None)

Bases: object

This is the Bayesian calibration class.

A Bayesian calibration requires the following input:

  1. The DynamicSystem that encapsulates the observation data and simulation data,

  2. The inference method, for example, IterativeBayesianFilter,

  3. The number of iterations

4. The current iteration number if the user simply wants to process their data with GrainLearning for one iteration # or continue from that iteration (TODO)

# 5. A tolerance to stop the iterations if the maximum uncertainty is below the tolerance

  1. and the flag for skipping (-1), showing (0), or saving (1) the figures.

There are two ways of initializing a calibration toolbox class.

Method 1 - dictionary style (recommended)

bayesian_calibration = BayesianCalibration.from_dict(
    {
        "num_iter": 8,
        "callback": run_sim,
        "system": {
            "system_type": DynamicSystem,
            "model_name": "test",
            "param_names": ["a", "b"],
            "param_min": [0, 0],
            "param_max": [1, 10],
            "num_samples": 10,
            "obs_data": [2,4,8,16],
            "ctrl_data": [1,2,3,4],
        },
        "calibration": {
            "inference": {"ess_target": 0.3},
            "sampling": {"max_num_components": 1},
        },
        "save_fig": -1,
    }
)

or

Method 2 - class style

bayesian_calibration = BayesianCalibration(
    num_iter = 8,
    system = DynamicSystem(...),
    calibration = IterativeBayesianFilter(...)
    save_fig = -1
)
Parameters:

  • system (Type[DynamicSystem]) – A dynamic system whose observables and hidden states evolve dynamically over “time”

  • calibration (Type[IterativeBayesianFilter]) – An iterative Bayesian Filter that iteratively sample the parameter space

  • num_iter (int) – Number of iteration steps

  • curr_iter (int) – Current iteration step

  • save_fig (int) – Flag for skipping (-1), showing (0), or saving (1) the figures

  • callback (Optional[Callable]) – A callback function that runs the external software and passes the parameter sample to generate outputs

classmethod from_dict(obj)

An alternative constructor to allow choosing a system type (e.g., dynamic system or IO dynamic system)

Parameters:

obj (Dict) – a dictionary containing the keys and values to construct a BayesianCalibration object

Returns:

a BayesianCalibration object

get_most_prob_params()

Return the most probable set of parameters

Returns:

Estimated parameter values

increase_curr_iter()

Increase the current iteration step by one

load_all()

Simply load all previous iterations of Bayesian calibration

load_and_process(sigma=0.1)

Load existing simulation data and compute the posterior distribution using an assumed sigma

Parameters:

sigma (float) – assumed uncertainty coefficient, defaults to 0.1

load_and_run_one_iteration()

Load existing simulation data and run Bayesian calibration for one iteration Note the maximum uncertainty sigma_max is solved to reach a certain effective sample size ess_target, unlike being assumed as an input for load_and_process(…)

load_system()

Load existing simulation data from disk into the dynamic system

plot_uq_in_time()

Plot the evolution of uncertainty moments and distribution over time

resample()

Learn and resample from a proposal distribution todo this should be refactored

Returns:

Combinations of resampled parameter values

run()

This is the main calibration loop which does the following steps 1. First iteration of Bayesian calibration starts with a Halton sequence 2. Iterations continue by resampling the parameter space until certain criteria are met.

run_callback()

Run the callback function

run_one_iteration(index=-1)

Run Bayesian calibration for one iteration.

Parameters:

index (int) – iteration step, defaults to -1

Dynamic systems

This module contains various classes of state-space systems that describe the time evolution of the system’s hidden state based on partial observations from the real world.

class DynamicSystem(obs_data, num_samples, param_min, param_max, ctrl_data=None, obs_names=None, ctrl_name=None, inv_obs_weight=None, sim_name='sim', sim_data=None, curr_iter=0, param_data=None, param_names=None, sigma_max=1000000.0, sigma_tol=0.001)

Bases: object

This is the dynamical system class.

A dynamical system (also known as a state-space system) describes the time evolution of the system’s (hidden) state and observation using the following equations:

\[ \begin{align}\begin{aligned}x_t & = f(x_{t−1}) + q_{t−1}\\y_t & = h(x_t) + r_t\end{aligned}\end{align} \]

where \(x_t\) is the hidden state, \(y_t\) is the observation, both represented as random processes, \(f\) is the state transition function, \(h\) is the observation function, \(q_{t−1}\) is the process noise, and \(r_t\) is the observation noise.

In the context of Bayesian parameter estimation, \(f\) is the model that describe the physical process and \(h\) is the model that describe the relationship between the hidden state and observation. In the simplest case, \(h\) is an identity matrix which indicate the one-to-one relationship between \(x_t\) and \(y_t\).

Therefore, the DynamicSystem class is used to encapsulate the observation data and the simulation data, which require specifying the number of samples and the lower and upper bound of the parameters.

There are two ways of initializing the class.

Method 1 - dictionary style

system_cls = DynamicSystem.from_dict(
    {
        "param_names": ["a", "b"],
        "param_min": [0, 0],
        "param_max": [1, 10],
        "num_samples": 14,
        "obs_data": y_obs,
        "ctrl_data": y_ctrl,
    }
)

or

Method 2 - class style

system_cls = DynamicSystem(
    param_names = ["a", "b"],
    param_min = [0, 0],
    param_max = [1, 10],
    num_samples = 14,
    obs_data = y_obs,
    ctrl_data = y_ctrl,
)

You can pass the simulation data to the class using the DynamicSystem.set_sim_data() method.

system_cls.set_sim_data(x)

The simulation data is a numpy array of shape (num_samples, num_obs, num_steps).

Parameters:

  • obs_data (ndarray) – observation or reference data

  • num_samples (int) – Sample size of the ensemble of model evaluations

  • num_steps – Number of steps or sequence size in the dataset

  • num_obs – Number of observations in the dataset

  • num_ctrl – Number of control data (identical between simulation and observation)

  • param_min (List[float]) – List of parameter lower bounds

  • param_max (List[float]) – List of parameter Upper bounds

  • curr_iter (int) – current iteration ID, defaults to 0

  • ctrl_data (Optional[ndarray]) – control data (e.g, time), defaults to None, optional

  • obs_names (Optional[List[str]]) – Column names of the observation data, defaults to None, optional

  • ctrl_name (Optional[str]) – Column name of the control data, defaults to None, optional

  • inv_obs_weight (Optional[List[float]]) – Inverse of the observation weight, defaults to None, optional

  • param_data (Optional[ndarray]) – Parameter data, defaults to None, optional

  • param_names (Optional[List[str]]) – Parameter names, defaults to None, optional

  • sim_data (Optional[ndarray]) – Simulation data, defaults to None, optional

  • sigma_max (float) – Maximum uncertainty, defaults to 1.0e6, optional

  • sigma_tol (float) – Tolerance of the estimated uncertainty, defaults to 1.0e-3, optional

  • sim_name (str) – Name of the simulation, defaults to ‘sim’, optional

  • sigma_min – Minimum uncertainty, defaults to 1.0e-6, optional

  • _inv_normalized_sigma – Calculated normalized sigma to weigh the covariance matrix

  • estimated_params – Estimated parameter as the first moment of the distribution (\(x_\mu = \sum_i w_i * x_i\)), defaults to None, optional

  • estimated_params_cv – Estimated parameter coefficient of variation as the second moment of the distribution (\(x_\sigma = \sqrt{\sum_i w_i * (x_i - x_\mu)^2} / x_\mu\)), defaults to None, optional

classmethod backup_sim_data()

Virtual function to backup simulation data

compute_estimated_params(posteriors)

Compute the estimated means and uncertainties for the parameters.

This function is vectorized for all time steps

Parameters:

posteriors (array) – Posterior distribution of shape (num_steps, num_samples)

compute_inv_normalized_sigma()

Compute the inverse of the matrix that apply different weights on the observables

classmethod from_dict(obj)

Initialize the class using a dictionary style

Parameters:

obj (dict) – Dictionary object

Returns:

DynamicSystem: DynamicSystem object

get_inv_normalized_sigma()

Get the inverse of the matrix that apply different weights on the observables

classmethod get_sim_data_files()

Virtual function to get simulation data files from disk

classmethod load_param_data()

Virtual function to load param data from disk

classmethod load_sim_data()

Virtual function to load simulation data

classmethod move_data_to_sim_dir()

Virtual function to move data into simulation directory

reset_inv_normalized_sigma()

Reset the inverse of the weighting matrix to None

set_ctrl_data(data)

Set the control data :type data: list :param data: control data of shape (num_ctrl, num_steps)

Parameters:

data (list) –

set_obs_data(data)

Set the observation data :type data: list :param data: observation data of shape (num_obs, num_steps)

Parameters:

data (list) –

set_param_data(data)

Set the simulation data

Parameters:

data (list) – parameter data of shape (num_samples, num_params)

set_sim_data(data)

Set the simulation data

Parameters:

data (list) – simulation data of shape (num_samples, num_obs, num_steps)

classmethod set_up_sim_dir()

Virtual function to set up simulation directory

classmethod write_params_to_table()

Write the parameter data into a text file

class IODynamicSystem(sim_name, sim_data_dir, sim_data_file_ext, obs_data_file, obs_names, ctrl_name, num_samples, param_min, param_max, obs_data=None, ctrl_data=None, inv_obs_weight=None, sim_data=None, curr_iter=0, param_data_file='', param_data=None, param_names=None)

Bases: DynamicSystem

This is the I/O dynamic system class derived from the dynamic system class. Extra functionalities are added to handle I/O operations.

There are two ways of initializing the class.

Method 1 - dictionary style

system_cls = IODynamicSystem.from_dict(
    {
        "system_type": IODynamicSystem,
        "param_min": [0, 0],
        "param_max": [1, 10],
        "param_names": ['a', 'b'],
        "num_samples": 14,
        "obs_data_file": 'obs_data.txt',
        "obs_names": ['y_obs'],
        "ctrl_name": 'y_ctrl
        "sim_name": 'linear',
        "sim_data_file_ext": '.txt',
    }
)

or

Method 2 - class style

system_cls = IODynamicSystem(
        param_min = [0, 0],
        param_max = [1, 10],
        param_names = ['a', 'b'],
        num_samples = 14,
        obs_data_file = 'obs_data.txt',
        obs_names = ['y_obs'],
        ctrl_name = 'y_ctrl',
        sim_name = 'linear',
        sim_data_file_ext = '.txt',
)
Parameters:

  • param_min (List[float]) – List of parameter lower bounds

  • param_max (List[float]) – List of parameter Upper bounds

  • param_names (Optional[List[str]]) – Parameter names, defaults to None

  • num_samples (int) – Sample size of the ensemble of model evaluations

  • obs_data_file (str) – Observation data file, defaults to None

  • obs_names (List[str]) – Column names of the observation data, defaults to None

  • ctrl_name (str) – Column name of the control data, defaults to None

  • sim_name (str) – Name of the simulation, defaults to ‘sim’

  • sim_data_dir (str) – Simulation data directory, defaults to ‘./sim_data’

  • sim_data_file_ext (str) – Simulation data file extension, defaults to ‘.npy’

  • curr_iter (int) – Current iteration ID, defaults to 0

  • param_data_file (str) – Parameter data file, defaults to None, optional

  • obs_data (Optional[ndarray]) – observation or reference data, optional

  • ctrl_data (Optional[ndarray]) – Control data (e.g, time), defaults to None, optional

  • inv_obs_weight (Optional[List[float]]) – Inverse of the observation weight, defaults to None, optional

  • param_data (Optional[ndarray]) – Parameter data, defaults to None, optional

  • sim_data (Optional[ndarray]) – Simulation data, defaults to None, optional

  • sim_data_files – List of simulation data files (num_samples), defaults to None, optional

backup_sim_data()

Backup simulation data files to a backup directory.

compute_estimated_params(posteriors)

Compute the estimated means and uncertainties for the parameters.

This function is vectorized for all time steps

Parameters:

posteriors (array) – Posterior distribution of shape (num_steps, num_samples)

compute_inv_normalized_sigma()

Compute the inverse of the matrix that apply different weights on the observables

classmethod from_dict(obj)

Initialize the class using a dictionary style

Parameters:

obj (dict) – Dictionary object

Return IODynamicSystem:

IODynamicSystem object

get_inv_normalized_sigma()

Get the inverse of the matrix that apply different weights on the observables

get_obs_data()

Get the observation data from the observation data file.

Separate the control data from the observation data if the name of control variable is given. Otherwise, the observation data is the entire data in the observation data file.

get_sim_data_files()

Get the simulation data files from the simulation data directory.

load_param_data()

Load parameter data from a table written in a text file.

load_sim_data()

Load the simulation data from the simulation data files.

The function does the following: 1. Load simulation data into an IO dynamic system object 2. Check if parameter values read from the table matches those used to creat the simulation data

move_data_to_sim_dir()

Move simulation data files and corresponding parameter table into the directory defined per iteration

reset_inv_normalized_sigma()

Reset the inverse of the weighting matrix to None

set_ctrl_data(data)

Set the control data :type data: list :param data: control data of shape (num_ctrl, num_steps)

Parameters:

data (list) –

set_obs_data(data)

Set the observation data :type data: list :param data: observation data of shape (num_obs, num_steps)

Parameters:

data (list) –

set_param_data(data)

Set the simulation data

Parameters:

data (list) – parameter data of shape (num_samples, num_params)

set_sim_data(data)

Set the simulation data

Parameters:

data (list) – simulation data of shape (num_samples, num_obs, num_steps)

set_up_sim_dir()

Create a directory to store simulation data and write the parameter data into a text file

write_params_to_table()

Write the parameter data into a text file.

Return param_data_file:

The name of the parameter data file

Iterative Bayesian Filter

This module contains various Bayesian filtering classes by mixing the inference and sampling methods from the inference and sampling modules.

class IterativeBayesianFilter(inference=None, sampling=None, ess_tol=0.01, initial_sampling='halton', proposal=None, proposal_data_file=None)

Bases: object

This is the Iterative Bayesian Filter class.

The idea is to solve an inverse problem repeatedly by refining the previous knowledge of the dynamical system. An iterative Bayesian filter is essentially a sequential Monte Carlo filter that iteratively samples the parameter space, using a proposal distribution trained with non-Gaussian mixture models. It consists of the inference class and the (re)sampling class. The inference class is a sequential Monte Carlo filter, and the resampling class is a Gaussian mixture model. The Gaussian mixture model is trained with the previous ensemble (i.e., samples and associated weights).

The steps of the iterative Bayesian filter are as follows:

  1. Initialize the ensemble of model (including parameter) states using a low-discrepancy sequence.

  2. Compute the posterior distribution of model states such that the target effective sample size is reached.

  3. Generate new samples from a proposal density that is trained with the previous ensemble (i.e., samples and associated weights).

  4. Run the model with the new samples and update the ensemble of model states using DynamicSystem.run() and BayesianCalibration.load_system().

  5. Repeat steps 2-4 until convergence. The iterations are done in the Bayesian Calibration class by calling BayesianCalibration.run_one_iteration().

There are two ways of initializing the class.

Method 1 - dictionary style

ibf_cls = IterativeBayesianFilter.from_dict(
    {
        "inference":{
            "ess_target": 0.3,
            "scale_cov_with_max": True
        },
        "sampling":{
            "max_num_components": 2
        }
    }
)

or

Method 2 - class style

system_cls = IterativeBayesianFilter(
        inference = SMC(...),
        sampling = GaussianMixtureModel(...)
)
Parameters:

  • inference (Optional[Type[SMC]]) – Sequential Monte Carlo class (SMC)

  • sampling (Optional[Type[GaussianMixtureModel]]) – Gaussian Mixture Model class (GMM)

  • initial_sampling (str) – The initial sampling method, defaults to Halton

  • ess_tol (float) – Tolerance for the target effective sample size to converge, defaults to 1.0e-2

  • proposal (Optional[ndarray]) – A proposal distribution to sample the state-parameter space, defaults to None

  • proposal_data_file (Optional[str]) – Pickle that stores the previously trained proposal distribution, defaults to None

  • param_data_list – List of the parameter samples of shape (num_samples, num_params) generated in all iterations

  • sigma_list – List of sigma values optimized to satisfy the target effective sample size in all iterations

  • posterior – The posterior distribution of model states at the last time step

add_curr_param_data_to_list(param_data)

Add the current parameter samples to the list of parameter samples.

Parameters:

param_data (ndarray) – Current parameter samples of shape (num_samples, num_params)

classmethod from_dict(obj)

Initialize the class using a dictionary style

Parameters:

obj (dict) – a dictionary containing the keys and values to construct an iBF object

Returns:

an iBF object

initialize(system)

Initialize the ensemble of model (including parameter) states using a low-discrepancy sequence.

Parameters:

system (Type[DynamicSystem]) – Dynamic system class

load_proposal_from_file(system)

Load the proposal density from a file.

Parameters:

system (Type[IODynamicSystem]) – Dynamic system class

run_inference(system)

Compute the posterior distribution of model states such that the target effective sample size is reached.

Parameters:

system (Type[DynamicSystem]) – Dynamic system class

run_sampling(system)

Generate new samples from a proposal density.

The proposal density can be trained with the previous ensemble (i.e., samples and associated weights).

Parameters:

system (Type[DynamicSystem]) – Dynamic system class

save_proposal_to_file(system)

Save the proposal density to a file.

Parameters:

system (Type[IODynamicSystem]) – Dynamic system class

solve(system)

Run both inference and sampling for a dynamic system

The iterations are done at the level of the Bayesian Calibration class by calling BayesianCalibration.run_one_iteration().

Parameters:

system (Type[DynamicSystem]) – Dynamic system class

Inference

This module contains various methods for performing statistical and Bayesian inference

class SMC(ess_target, scale_cov_with_max=False, cov_matrices=None)

Bases: object

This is the Sequential Monte Carlo class that recursively update model states based on sequential observations using the Bayes’ theorem

There are two ways of initializing the class.

Method 1 - dictionary style

system_cls = SMC.from_dict(
    {
        "ess_target": 0.3,
        "scale_cov_with_max": False
    }
)

or

Method 2 - class style

system_cls = SMC(
        ess_target = 0.3,
        scale_cov_with_max = False
)
Parameters:

  • ess_target (float) – Target effective sample size (w_0 / sum_i w_i^2) where w_0 = 1 / N_p, defaults to 0.3

  • scale_cov_with_max (bool) – True if the covariance matrix is scaled with the maxima of the observations, defaults to True

  • cov_matrices (Optional[array]) – Covariance matrices of shape (num_steps, num_obs, num_obs), defaults to None, Optional

  • likelihoods – Likelihood distributions of shape (num_steps, num_samples)

  • posteriors – Posterior distributions of shape (num_steps, num_samples)

  • ess – Time evolution of the effective sample size

compute_effective_sample_size()

Compute the effective sample size

data_assimilation_loop(sigma, system, proposal=None)

Perform data assimilation loop

Parameters:

  • sigma (float) – Uncertainty

  • proposal (Optional[ndarray]) – Proposal distribution from which the samples are sampled from, defaults to None, optional

  • system (Type[DynamicSystem]) – Dynamic system class

Returns:

Result of the objective function which converges to a user defined effective sample size

classmethod from_dict(obj)

Initialize the class using a dictionary style

Parameters:

obj (dict) – a dictionary containing the keys and values to construct an SMC object

Returns:

an SMC object

get_covariance_matrices(sigma, system)

Compute the (diagonal) covariance matrices from an uncertainty that is either assumed by the user or updated by SMC to satisfy the target effective sample size.

This function is vectorized for all time steps

Parameters:

  • sigma (float) – Uncertainty

  • system (Type[DynamicSystem]) – Dynamic system

Returns:

Covariance matrices for all time steps

static get_likelihoods(system, cov_matrices)

Compute the likelihood distributions of simulation data as a multivariate normal centered around the observation.

This function is vectorized for all time steps

Parameters:

  • system (Type[DynamicSystem]) – Dynamic system class

  • cov_matrices (array) – Covariance matrices of shape (num_steps, num_obs, num_obs)

Returns:

Likelihood matrices of shape (num_steps, num_samples) considering all time steps

get_posterior_at_time(time_step=-1)

Get the posterior distribution at a certain time step

Parameters:

time_step (int) – input time step, defaults to -1 (the last step in time), optional

Returns:

Posterior distribution at a certain time step

static get_posteriors(system, likelihoods, proposal=None)

Compute the posterior distributions from the likelihood for all the time steps

This function is vectorized for all time steps

Parameters:

  • system (Type[DynamicSystem]) – Dynamic system class

  • likelihoods (array) – Likelihood matrices of shape (num_steps, num_samples)

  • proposal (Optional[array]) – Proposal distribution at the initial time step, defaults to None, optional

Returns:

Posterior matrices of shape (num_steps, num_samples) considering all time steps

set_posteriors(posteriors=None)

Set the posterior distribution

Parameters:

posteriors (Optional[ndarray]) –

Sampling

This module contains various methods to sample the state-parameter space of a dynamic system.

class GaussianMixtureModel(max_num_components=0, weight_concentration_prior=0.0, covariance_type='tied', n_init=1, tol=1e-05, max_iter=100, random_state=None, init_params='k-means++', warm_start=False, expand_factor=10, slice_sampling=False, deviation_factor=0.0)

Bases: object

This is the wrapper class of the Bayesian Gaussian Mixture. from scikit-learn.

The BayesianGaussianMixture. class is extended to include the following functionalities:

There are two ways of initializing the class.

Method 1 - dictionary style

system_cls = GaussianMixtureModel.from_dict(
    {
        "max_num_components": 2,
        "covariance_type": "full",
        "random_state": 0,
        slice_sampling: False,
    }
)

or

Method 2 - class style

system_cls = GaussianMixtureModel(
    max_num_components = 2,
    covariance_type = "full",
    random_state = 0,
    slice_sampling = False,
)
Parameters:

  • max_num_components – Maximum number of components

  • weight_concentration_prior (float) – The dirichlet concentration of each component on the weight distribution (Dirichlet), default to None, optional This is commonly called gamma in the literature. The higher concentration puts more mass in the center and will lead to more components being active, while a lower concentration parameter will lead to more mass at the edge of the mixture weights simplex. The value of the parameter must be greater than 0. If it is None, it’s set to 1. / n_components.

  • covariance_type (str) – {‘full’, ‘tied’, ‘diag’, ‘spherical’}, defaults to “full”, optional String describing the type of covariance parameters to use. Must be one of: 'full': each component has its own general covariance matrix. 'tied': all components share the same general covariance matrix. 'diag': each component has its own diagonal covariance matrix. 'spherical': each component has its own single variance.

  • n_init (int) – number of initialization to perform, defaults to 1, optional The result with the highest lower bound value on the likelihood is kept.

  • tol (float) – tolerance threshold, defaults to 1.0e-3, optional EM iterations will stop when the lower bound average gain on the likelihood (of the training data with respect to the model) is below this threshold.

  • max_iter (int) – maximum number of EM iterations to perform, defaults to 100, optional

  • random_state (Optional[int]) – random seed given to the method chosen to initialize the parameters, defaults to None, optional

  • init_params (str) – {‘kmeans’, ‘k-means++’, ‘random’, ‘random_from_data’}, default=’kmeans’, optional The method used to initialize the weights, the means and the covariances. String must be one of: 'kmeans': responsibilities are initialized using kmeans. k-means++: responsibilities are initialized using a k-means++ algorithm. 'random': responsibilities are initialized randomly. 'random_from_data': responsibilities are initialized randomly from the data.

  • warm_start (bool) – flag to use warm start, defaults to False, optional If ‘warm_start’ is True, the solution of the last fitting is used as initialization for the next call of fit(). This can speed up convergence when fit is called several times on similar problems. See the Glossary.

  • expand_factor (int) – factor used when converting the ensemble from weighted to unweighted, defaults to 10, optional

  • slice_sampling (bool) – flag to use slice sampling, defaults to False, optional

  • deviation_factor (float) – factor used to determine the threshold, scores.mean() - deviation_factor * scores.std(), for slice sampling, defaults to 0.0, optional

  • gmm – The class of the Gaussian Mixture Model

  • max_params – Current maximum values of the parameters

draw_samples_within_bounds(system, num=1)

Draw new parameter samples within the user-defined upper and lower bounds

Parameters:

  • system (Type[DynamicSystem]) – Dynamic system class

  • num (int) – Number of samples to draw

Returns:

New parameter samples

expand_and_normalize_weighted_samples(weights, system)

Converting an ensemble of samples with importance weights into an ensemble of samples with equal weights.

For instance, a sample is duplicated a few times depending on the importance weight associated to it.

Parameters:

  • weights (ndarray) – Importance weights associated with the ensemble

  • system (Type[DynamicSystem]) – Dynamic system class

Returns:

Expanded unweighted samples

classmethod from_dict(obj)

Initialize the class using a dictionary style

Parameters:

obj (dict) – a dictionary containing the keys and values to construct a GMM object

Returns:

a GMM object

load_gmm_from_file(proposal_data_file='proposal_density.pkl')

Load the Gaussian mixture model from a file.

Parameters:

proposal_data_file (str) – Name of the file that stores the trained Gaussian mixture model

regenerate_params(weight, system)

Regenerating new samples by training and sampling from a Gaussian mixture model.

Parameters:

  • weight (ndarray) – Posterior found by the data assimilation

  • system (Type[DynamicSystem]) – Dynamic system class

Returns:

Resampled parameter data

save_gmm_to_file(proposal_data_file='proposal_density.pkl')

Save the Gaussian mixture model to a file.

Parameters:

proposal_data_file (str) –

train(weight, system)

Train the Gaussian mixture model.

Parameters:

  • weight (ndarray) – Posterior found by the data assimilation

  • system (Type[DynamicSystem]) – Dynamic system class

generate_params_qmc(system, num_samples, method='halton', seed=None)

This is the class to uniformly draw samples in n-dimensional space from a low-discrepancy sequence or a Latin hypercube.

See Quasi-Monte Carlo.

Parameters:

  • system (Type[DynamicSystem]) – Dynamic system class

  • num_samples (int) – Number of samples to draw

  • method (str) – Method to use for Quasi-Monte Carlo sampling. Options are “halton”, “sobol”, and “LH”

  • seed – Seed for the random number generator

RNN Module

data pre-processing

class Preprocessor

Bases: ABC

classmethod get_dimensions(data)

Extract dimensions of sample from a tensorflow dataset.

Parameters:

data (DatasetV2) – The dataset to extract from.

Returns:

Dictionary containing: sequence_length, num_load_features, num_contact_params, num_labels

classmethod get_split(dataset, inds)

Given a dataset and the indexes, return the sub-dataset containing only the elements of indexes in inds. The returned dataset respects the format of inputs and outputs (see more info in return).

Parameters:

  • dataset (tuple) – Full dataset to split on, in form of a tuple (inputs, outputs), where inputs is a dictionary and its keys and outputs are numpy arrays.

  • inds (array) – Indexes of the elements in dataset that are going to be gathered in this specific split.

Returns:

tuple of the split, 2 dimensions: * inputs: dictionary containing 'load_sequence' and 'contact_params'. * outputs: Labels

classmethod make_splits(dataset, train_frac, val_frac, seed)

Split data into training, validation, and test sets. The split is done on a sample by sample basis, so sequences are not broken up. It is done randomly across all pressures and experiment types present.

Parameters:

  • dataset (tuple) – Full dataset to split on, in form of a tuple (inputs, outputs), where inputs is a dictionary and its keys and the outputs are numpy arrays.

  • train_frac (float) – Fraction of data used for training set.

  • val_frac (float) – Fraction of data used for validation set. Test fraction is the remaining.

  • seed (int) – Random seed used to make the split.

Returns:

Dictionary containing 'train', 'val', and 'test' datasets.

classmethod pad_initial(array, pad_length, axis=1)

Add pad_length copies of the initial state in the sequence to the start. This is used to be able to predict also the first timestep from a window of the same size.

Parameters:

  • array (array) – Array that is going to be modified

  • pad_lenght – number of copies of the initial state to be added at the beggining of array.

  • pad_length (int) –

Returns:

Modified array

abstract prepare_datasets(**kwargs)

Abstract class to be implemented in each of the child classes. Must return a tuple containing:

  • data_split: Dictionary with keys 'train', 'val', 'test', and values of the corresponding tensorflow Datasets.

  • train_stats: Dictionary containing the shape of the data:

    sequence_length, num_input_features, num_contact_params, num_labels, and 'mean' and 'std' of the training set, in case standardize_outputs is True.

abstract prepare_single_dataset()

Abstract class to be implemented in each of the child classes. Must return a Tensorflow dataset containing a Tuple (inputs, outputs)

  • inputs: Dictionary with keys 'load_sequence', 'contact_parameters', that are the corresponding tensorflow Datasets to input to an rnn model.

  • outputs: tensorflow Dataset containing the outputs or labels.

classmethod standardize_outputs(split_data)

Standardize outputs of split_data using the mean and std of the training data taken over both the samples and the timesteps. The 3 datasets 'train', 'val', 'test' will be standardized.

Parameters:

split_data (dict) – dictionary containing 'train', 'val', 'test' keys and the respective datasets.

Returns:

Tuple containing:

  • standardized_splits: Dictionary containing the standardized datasets.

  • train_stats: Dictionary with the metrics (mean and standard deviation) used to standardize the data.

classmethod warning_config_field(key, config, default, add_default_to_config=False)

Raises a warning if key is not included in config dictionary. Also informs the default value that will be used. If add_default_to_config=True, then it adds the key and its default value to config.

Parameters:

  • key (str) – string of the key in dictionary config to look for.

  • config (dict) – Dictionary containing keys and values

  • default – default value associated to the key in config.

  • add_default_to_config (bool) – Wether to add the default value and the key to config in case it does not exists.

Returns:

Updated config dictionary.

class PreprocessorTriaxialCompression(**kwargs)

Bases: Preprocessor

Class to Preprocess data of triaxial compression experiments, inheriting from abstract class Preprocessor

Attributes: see PreprocessorTriaxialCompression.get_default_config()

classmethod check_config(config)

Checks that values requiring an input from the user would be specified in config.

Parameters:

config (dict) – Dictionary containing the values of different arguments.

Returns:

Updated config dictionary.

classmethod get_default_config()

Returns a dictionary with default values for the configuration of data preparation. Possible fields are:

  • 'raw_data': Path to hdf5 file generated using parse_data_YADE.py

  • 'pressure' and 'experiment_type': Name of the subfield of dataset to consider. It can also be ‘All’.

  • 'standardize_outputs': If True transform the data labels to have zero mean and unit variance. Also, in train_stats the mean and variance of each label will be stored, so that can be used to transform predicitons. (This is very usful if the labels are not between [0,1])

  • 'add_e0': Whether to add the initial void ratio (output) as a contact parameter.

  • 'add_pressure': Wether to add the pressure to contact_parameters.

  • 'add_experiment_type': Wether to add the experiment type to contact_parameters.

  • 'train_frac': Fraction of the data used for training, between [0,1].

  • 'val_frac': Fraction of the data used for validation, between [0,1]. The fraction of the data used for test is then 1 - train_frac - val_frac.

  • 'window_size': int, number of steps composing a window.

  • 'window_step': int, number of steps between consecutive windows (default = 1).

  • 'pad_length': int, equals to window_size. Length of the sequence that with be pad at the start.

Returns:

Dictionary containing default values of the arguments that the user can set.

classmethod get_dimensions(data)

Extract dimensions of sample from a tensorflow dataset.

Parameters:

data (DatasetV2) – The dataset to extract from.

Returns:

Dictionary containing: sequence_length, num_load_features, num_contact_params, num_labels

classmethod get_split(dataset, inds)

Given a dataset and the indexes, return the sub-dataset containing only the elements of indexes in inds. The returned dataset respects the format of inputs and outputs (see more info in return).

Parameters:

  • dataset (tuple) – Full dataset to split on, in form of a tuple (inputs, outputs), where inputs is a dictionary and its keys and outputs are numpy arrays.

  • inds (array) – Indexes of the elements in dataset that are going to be gathered in this specific split.

Returns:

tuple of the split, 2 dimensions: * inputs: dictionary containing 'load_sequence' and 'contact_params'. * outputs: Labels

classmethod make_splits(dataset, train_frac, val_frac, seed)

Split data into training, validation, and test sets. The split is done on a sample by sample basis, so sequences are not broken up. It is done randomly across all pressures and experiment types present.

Parameters:

  • dataset (tuple) – Full dataset to split on, in form of a tuple (inputs, outputs), where inputs is a dictionary and its keys and the outputs are numpy arrays.

  • train_frac (float) – Fraction of data used for training set.

  • val_frac (float) – Fraction of data used for validation set. Test fraction is the remaining.

  • seed (int) – Random seed used to make the split.

Returns:

Dictionary containing 'train', 'val', and 'test' datasets.

classmethod pad_initial(array, pad_length, axis=1)

Add pad_length copies of the initial state in the sequence to the start. This is used to be able to predict also the first timestep from a window of the same size.

Parameters:

  • array (array) – Array that is going to be modified

  • pad_lenght – number of copies of the initial state to be added at the beggining of array.

  • pad_length (int) –

Returns:

Modified array

prepare_datasets(seed=42)

Convert raw data into preprocessed split datasets. First split the data into train, val and test datasets and then apply the Sliding windows transformation. This is to avoid having some parts of a dataset in train and some in val and/or in test (i.e. data leak).

Parameters:

seed (int) – Random seed used to split the datasets.

Returns:

Tuple (split_data, train_stats)

  • split_data: Dictionary with keys 'train', 'val', 'test', and values the corresponding tensorflow Datasets.

  • train_stats: Dictionary containing the shape of the data: sequence_length, num_load_features, num_contact_params, num_labels, and 'mean' and 'std' of the training set, in case standardize_outputs is True.

prepare_single_dataset()

Convert raw data into a tensorflow dataset with compatible format to predict and evaluate a rnn model.

Returns:

Tensorflow dataset containing a Tuple (inputs, outputs)

  • inputs: Dictionary with keys 'load_sequence', 'contact_parameters', that are the corresponding tensorflow Datasets to input to an rnn model.

  • outputs: tensorflow Dataset containing the outputs or labels.

classmethod standardize_outputs(split_data)

Standardize outputs of split_data using the mean and std of the training data taken over both the samples and the timesteps. The 3 datasets 'train', 'val', 'test' will be standardized.

Parameters:

split_data (dict) – dictionary containing 'train', 'val', 'test' keys and the respective datasets.

Returns:

Tuple containing:

  • standardized_splits: Dictionary containing the standardized datasets.

  • train_stats: Dictionary with the metrics (mean and standard deviation) used to standardize the data.

classmethod warning_config_field(key, config, default, add_default_to_config=False)

Raises a warning if key is not included in config dictionary. Also informs the default value that will be used. If add_default_to_config=True, then it adds the key and its default value to config.

Parameters:

  • key (str) – string of the key in dictionary config to look for.

  • config (dict) – Dictionary containing keys and values

  • default – default value associated to the key in config.

  • add_default_to_config (bool) – Wether to add the default value and the key to config in case it does not exists.

Returns:

Updated config dictionary.

models

Module containing a function that creates a RNN model.

rnn_model(input_shapes, window_size=20, lstm_units=50, dense_units=20, seed=42, **_)

Neural network with an LSTM layer.

Takes in a load sequence and contact parameters, and outputs the macroscopic responses. The contact parameters are used to initialize the hidden state of the LSTM.

Parameters:

  • input_shapes (dict) – Dictionary containing ‘num_load_features’, ‘num_contact_params’, ‘num_labels’. It can contain other keys but hese are the ones used here.

  • window_size (int) – Length of time window.

  • lstm_units (int) – Number of units of the hidden state of the LSTM.

  • dense_units (int) – Number of units used in the dense layer after the LSTM.

  • seed (int) – The random seed used to initialize the weights.

Returns:

A Keras model.

windows

extract_tensors(data)

Given a tensorflow Dataset extract all tensors.

Parameters:

data (DatasetV2) – Tensorflow dataset.

Returns:

3 numpy arrays: inputs, contacts, outputs.

predict_over_windows(inputs, model, window_size, sequence_length)

Take a batch of full sequences, iterate over windows making predictions. It splits up the sequence into windows of given length, each offset by one timestep, uses the model to make predictions on all of those windows, and concatenates the result into a whole sequence again. Note the length of the output sequence will be shorter by the window_size than the input sequence.

Parameters:

  • inputs (dict) – dict containing inputs: ‘load_sequence’ and ‘contact_parameters’, both being tensorflow.Tensor.

  • model (Model) – The model to predict with.

  • window_size (int) – Number of timesteps in a single window.

  • sequence_length (int) – Number of timesteps in a full sequence.

Returns:

tensorflow.Tensor of predicted sequences.

windowize_single_dataset(data, window_size=1, window_step=1, seed=42, **_)

Take a dataset ((inputs, contact_params),outputs) of sequences of shape N, S, L and output another dataset of shorter sequences of size window_size, taken at intervals window_step. The resulting dataset will have M=((sequence_length - window_size)/window_step) + 1) samples, corresponding to independent windows in the given dataset. For a given window taken from inputs, the output is the element of outputs sequence at the last position of the window.

Note

For a clear picture of what goes on here check Sliding windows section in the documentation.

Data is shuffled.

Parameters:

  • data (DatasetV2) – dataset of sequences of shape N, S, L.

  • window_size (int) – Size of the window.

  • window_step (int) – Offset between subsequent windows.

  • seed (int) – Random seed.

Returns:

  • inputs of shape: (M, window_size, num_load_features), with M >> N.

  • outputs of shape: (M, L_outputs)

windowize_train_val_test(split_data, window_size, window_step, **_)

Convert sequences into windows of given length. Leave test set untouched.

Parameters:

  • split_data (dict) – Dictionary with keys ‘train’, ‘val’, ‘test’ pointing to tensorflow datasets.

  • window_size (int) – Number of timesteps to include in a window.

  • window_step (int) – Offset between subsequent windows.

Returns:

windows: Dictionary of dataset splits, where the training and validation splits have been modified into windows.

train

Train a model to predict macroscopic features of a DEM simulation.

get_default_config()

Returns a dictionary with default values for the configuration of RNN model and training procedure. Possible fields are:

  • RNN model

    • 'window_size': int, number of steps composing a window.

    • 'window_step': int, number of steps between consecutive windows (default = 1).

    • 'pad_length': int, equals to window_size. Length of the sequence that with be pad at the start.

    • 'lstm_units': int, number of neurons or units in LSTM layer.

    • 'dense_units': int, number of neurons or units of dense layer.

  • Training procedure

    • 'patience': patience of tf.keras.callbacks.EarlyStopping.

    • 'epochs': Maximum number of epochs.

    • 'learning_rate': double, learning_rate of tf.keras.optimizers.Adam.

    • 'batch_size': Size of the data batches per training step.

    • 'save_weights_only': Boolean

      • True: Only the weights will be saved (Recommended fro compatibility across platforms).

      • False: The whole model will be saved.

Returns:

Dictionary containing default values of the arguments that the user can set.

train(preprocessor, config=None, model=None)

Train a model and report to weights and biases.

If called in the framework of a sweep: A sweep can be created from the command line using a configuration file, for example example_sweep.yaml, as: wandb sweep example_sweep.yaml And run with the line shown subsequently in the terminal. The config is loaded from the yaml file.

Parameters:

  • preprocessor (Preprocessor) – Preprocessor object to load and prepare the data.

  • config – dictionary containing model and training configurations.

  • model (Optional[Model]) – Keras model if None is passed then an rnn_model will be created (default).

Returns:

Same as tf.keras.Model.fit(): A History object. Its History.history attribute is a record of training loss values and metrics values at successive epochs, as well as validation loss values and validation metrics values.

train_without_wandb(preprocessor, config=None, model=None)

Train a model locally: no report to wandb. Saves either the model or its weight to folder outputs.

Parameters:

  • preprocessor (Preprocessor) – Preprocessor object to load and prepare the data.

  • config – dictionary containing taining hyperparameters and some model parameters.

  • model (Optional[Model]) – Keras model if None is passed then an rnn_model will be created (default).

Returns:

Same as tf.keras.Model.fit(): A History object. Its History.history attribute is a record of training loss values and metrics values at successive epochs, as well as validation loss values and validation metrics values.

predict

Module containing functions to load a trained RNN model and make a prediction.

get_best_run_from_sweep(entity_project_sweep_id, preprocessor=None)

Load the best performing model found with a weights and biases sweep. Also load the data splits it was trained on.

Parameters:

  • entity_project_sweep_id (str) – string of form <user>/<project>/<sweep_id>

  • preprocessor (Optional[Preprocessor]) – Preprocessor object to load and prepare the data. If None is given, then a PreprocessorTriaxialCompression will be considered

Returns:

  • model: The trained model with the lowest validation loss.

  • data: The train, val, test splits as used during training.

  • stats: Some statistics on the training set and configuration.

  • config: The configuration dictionary used to train the model.

get_pretrained_model(path_to_model)

Loads configuration, training statistics and model of a pretrained model.

Reads train_stats, and creates dataset.

Parameters:

path_to_model (str) – str or pathlib.Path to the folder where is stored.

Returns:

  • model: keras model ready to use.

  • train_stats: Array containing the values used to standardize the data (if config.standardize_outputs = True), and lenghts of sequences, load_features, contact_params, labels, window_size and window_step.

  • config: dictionary with the model configuration

load_config(path_to_model)

Searches for the configuration (of the model training) file in ‘path_to_model’. Read config.yaml into a python dictionary equivalent to config. config.yaml contains information about hyperparameters and model parameters, is generated in every run of wandb. Raises FileNotFoundError if there are not files matching possible formats.

Parameters:

path_to_model (Path) – str or pathlib.Path to the folder where is stored.

Returns:

Dictionary config.

load_model(path_to_model, train_stats, config)

Searches for the file containing the saved model in ‘path_to_model’. Raises FileNotFoundError if there are not files matching possible formats.

Parameters:

  • path_to_model (Path) – str or pathlib.Path to the folder where is stored.

  • train_stats (dict) – Dictionary containing different dimensions of the data used to trained and standardization values. Its contents are saved in ‘path_to_model/train_stats.npy’.

  • config (dict) – Dictionary containing the configuration of the training for such model.

Returns:

Keras model.

predict_macroscopics(model, data, train_stats, config, batch_size=256)

Use the given model to predict the features of the given data. If ‘standardize_outputs’ in config, rescale the predictions to their original units.

Parameters:

  • model (Model) – Keras RNN model

  • data (DatasetV2) – Tensorflow dataset containing inputs: ‘load_sequence’ and ‘contact_parameters’, and outputs.

  • train_stats (dict) – Dictionary containing statistics of the training set.

  • config (dict) – Dictionary containing the configuration with which the model was trained.

  • batch_size (int) – Size of batches to use.

Returns:

predictions: tf.Tensor containing the predictions in original units.