Optimize

To optimize the rate constants, a model needs to be created that can create a prediction of the concentrations, compare it with the experimental data, and subsequently calculate the corresponding errors. This model can then:

Create a guess of the rate constants.
Compare the result of the prediction with the experimental data
Weigh the errors
Save the intermediate results
It will then repeat 1-4 with new rate constants until the error no longer improves and the rate constants have converged to a stable value or the maximum number of iterations has been reached. The method that is used to find the next set of rate constants is the Nelder-Mead algorithm with adaptive parameters.

The optimize.RateConstantOptimizerTemplate is an abstract base class which has implemented the numbered steps above. However, the user must define the error functions and the exact methodology of creating a prediction. The definition of the error functions consists out of two parts. First the data (either experimental or predicted) must be converted to a curve of interest, such as: \(A / (A + A_{labeled})\). Subsequently a metric wil be used to evaluate the error between each curve of both datasets.

class optimize.RateConstantOptimizerTemplate(experimental, metric, raw_weights=None)

Enables easy optimization of a model, which must be semi-implemented by the user.

Note

The user must implement the following abstract functions:

create_prediction()
calculate_curves()

Parameters:

experimental (pd.DataFrame) – The experimental data
metric (Callable[[np.ndarray, np.ndarray], float]) – An error metric which takes as input two np.ndarrays for the keywords y_pred and y_true and returns a float. Lower values should indicate a better prediction.
raw_weights (Optional[dict[str, float]]) – A dictionary containing patterns and weight. Each pattern will be searched for in the errors. Upon a match the error yielded will be multiplied by its respective weight. If an error is matched with multiple patterns, the weight will be decreased in a multiplicative manner. If None (default), no weights will be applied.

Variables:

weights (Optional[np.ndarray]) – The final weight per error type.

abstract static create_prediction(x, x_description)

Creates a prediction of the system, given a set of parameters. For DRL experiments DRL.predict_concentration does most of the required calculations.

Parameters:

x (np.ndarray) – Contains all parameters, which are to be optimized.
x_description (list[str]) – The description of each parameter.

Returns:

A DataFrame that contains the predicted concentrations as a function of time.

Return type:

pd.DataFrame

abstract static calculate_curves(data)

Calculates the curves corresponding to the data (either experimental or predicted). The experimental curves will only be calculated upon initialization and are stored for subsequent use.

Parameters:: data (pd.DataFrame) – The data from which the curves should be calculated. Either experimental or predicted by create_prediction().
Returns:: A dictionary containing a description of each curve, and the corresponding curve.
Return type:: dict[str, np.ndarray]

calculate_errors(prediction)

Calculates the (unweighted) error caused by each error function. This is done by calculating the curves corresponding to the prediction, and comparing them with the experimental curves using the metric.

Parameters:: prediction (DataFrame) – The predicted concentrations as a function of time
Returns:: The unweighted errors of each error function.
Return type:: pd.Series
Raises:: ValueError – When the metric returns a nan value, a ValueError will be raised. The error message will detail information on which curves caused the metric to return nan.

weigh_errors(errors)

Weighs the errors.

Parameters:: errors (Series) – The unweighted errors
Returns:: The weighted errors
Return type:: pd.Series

calculate_total_error(errors)

Weighs and sums the errors.

Parameters:: errors (Series) – The unweighted errors
Returns:: The total error in the model.
Return type:: float

optimize(x0, x_description, x_bounds, path, metadata=None, maxiter=50000, resume_from_simplex=None, show_pbar=True, _overwrite_log=False)

Optimizes the system, utilizing a nelder-mead algorithm.

Parameters:

x0 (np.ndarray) – Parameters which are to be optimized. Always contains the rate constants.
x_description (list[str]) – Description of each parameter.
x_bounds (Bounds) – The scipy.optimize.bounds of each parameter.
path (str | pathlib.Path) – The path to the folder where the optimization progress should be stored. The entire path to the folder will be created if it did not exist.
metadata (Optional[dict]) – The metadata that should be saved alongside the solution. This data will be stored in the settings_info.json file.
maxiter (float) – The maximum number of iterations.
resume_from_simplex (np.ndarray) – When a simplex is given of size [N+1, N] where N is the number of parameters, the solution starts here. This can be used to resume the optimization process.
show_pbar (bool) – If True, shows a progress bar.
_overwrite_log (bool) – If True, the logs will be overwritten. Should only be used in test scripts to avoid accidental loss of data.

Returns:

All relevant metadata and progress on each iteration will be stored in path. It can be loaded and analyzed either by using the visualize.VisualizeModel, or load_optimization_progress which returns an optimize.OptimizedModel instance.

Return type:

None

optimize_multiple(path, n_runs, x_description, x_bounds, x0_bounds=None, x0_min=1e-06, n_jobs=1, backend='loky', **optimize_kwargs)

Optimizes the system, utilizing a nelder-mead algorithm, for a given number of runs. Each run has random starting positions for each parameter, which is distributed according to a loguniform distribution. This choice of distribution allows the model to sample each parameter equally a logarithmic scale, for example k1 will be sampled just as often between 10e-5 and 10e-4 as between 0.1 and 1, allowing a wide range of values to be examined. The bounds of the starting position (x0_bounds) can be separately controlled from the bounds the system is allowed to explore (x_bounds). If the given path already has an existing directory called ‘optimization_multiple_guess’, the optimization will be resumed from that point onwards.

Parameters:

path (str | pathlib.Path) – The path to the folder where the optimization progress should be stored.
n_runs (int) – The number of runs which are to be computed.
x_description (list[str]) – Description of each parameter.
x_bounds (Bounds) – The scipy.optimize.bounds of each parameter.
x0_bounds (Optional[Bounds]) – The scipy.optimize.bounds for the starting value of each parameter. By default, it is identical to the x_bounds. Lower bounds smaller than x0_min are set to x0_min. When the upper bound is 0, the corresponding x0 will also be set to 0. This disables the reaction.
x0_min (float) – The minimum value the lower bound of x0_bounds can take. Any values lower than it is set to x0_min.
n_jobs (int) – The number of processes which should be used, allowing negative indexing, i.e. n_jobs=-1 would use up all available cores, n_jobs=-2 would use all available cores except one, etc.
backend (str) – The backend that is used by Joblib. Loky (default) works on all platforms.
**optimize_kwargs – The key word arguments that will be passed to self.optimize.

Returns:

All data of run n will be stored at ‘path/guess_n/’.

Return type:

None

static load_optimized_model(path)

Loads in the data from the log files.

Parameters:: path (str | pathlib.Path) – The path to the folder where the optimization progress has been stored.
Returns:: A structured OptimizerProgress instance which contains all information that was logged.
Return type:: OptimizedModel

class optimize.OptimizedModel(path)

Formatted data structure of an optimized model.

Parameters:

path (str | pathlib.Path) – The path to the folder where the optimization progress has been stored.

Variables:

metadata (dict[str, any]) – The stored metadata.
x_description (list[str]) – The description of each parameter.
n_dimensions (int) – The number of parameters.
n_function_evaluations (int) – The number of computed iterations.
all_x (pd.DataFrame) – The applied parameters for each iteration.
all_errors (pd.Series) – The error for each iteration.
all_times (pd.Series) – The timestamp at which the computation of each iteration was finished.
optimal_x (pd.Series) – The parameters for the iteration with the lowest error.
optimal_error (float) – The error corresponding to the iteration with the lowest error.
simplex (np.ndarray) – An array of size [N + 1, N] corresponding to the N parameters for each of the N + 1 best iterations

class optimize.OptimizedMultipleModels(path)

Formatted data structure for multiple optimized models. The data will be formatted such that the best run starts at index 0, and the error ascends with each index number.

Parameters:

path (str | pathlib.Path) – The path to the directory, which contains subdirectories. Each subdirectory in turn should contain its optimization progress.

Variables:

x_description (list[str]) – The description of each parameter.
all_initial_x (pd.DataFrame) – The initial parameters, x, in each model.
all_optimal_x (pd.DataFrame) – The optimal parameters, x, in each model.
all_optimal_error (np.ndarray) – The error when using the optimal parameters in each model.
best (OptimizedModel) – The model with the lowest error.

get_model(n)

Returns the n-th model, when sorted by the error (the best model at index 0).

Parameters:: n (int) – The index of the model that is to be retrieved.
Return type:: OptimizedModel

Example

To optimize the rate constants we will analyze a simple system. We can describe the model and the experimental conditions as follows:

\[\require{mhchem} \ce{catalyst + A <=>[\ce{k1}][\ce{k_{-1}}] B ->[\ce{{k2}}] C + catalyst}\]

import numpy as np
reactions = [
    ('k1', ['A', 'cat'], ['B'],),
    ('k-1', ['B'], ['A', 'cat'],),
    ('k2', ['B'], ['C', 'cat']),

    # labeled
    ('k1', ['A-d10', 'cat'], ['B-d10'],),
    ('k-1', ['B-d10'], ['A-d10', 'cat'],),
    ('k2', ['B-d10'], ['C-d10', 'cat'])
]

# look at as simple of a system as possible.
concentration_initial = {'A': 1, 'cat': 1 / 5}
concentration_labeled = {'A-d10': 1}
dilution_factor = 1
time_pre = np.linspace(0, 10, 50)
time_post = np.linspace(10, 90, 8 * 50)

To be able to fit to this hypothetical reaction fake data was generated by first creating a prediction using the DRL class. To this data we than add noise based on its intensity and a base level of noise to mimic real experimental data.

from delayed_reactant_labeling.predict import DRL
import pandas as pd

# create a "real" prediction.
rate_constants_real = {'k1': 0.3, 'k-1': 0.05, 'k2': 0.5}
drl_real = DRL(rate_constants=rate_constants_real, reactions=reactions)
real_data = drl_real.predict_concentration(
    t_eval_pre=time_pre,
    t_eval_post=time_post,
    dilution_factor=dilution_factor,
    initial_concentrations=concentration_initial,
    labeled_concentration=concentration_labeled)

# add noise
rng = np.random.default_rng(42)
fake_data = []
for col in real_data.columns[:-1]:  # the last column contains a time array, so skip that one.
    noise_dynamic = rng.normal(loc=1, scale=0.05, size=real_data[col].shape)  # fraction error
    noise_static = rng.normal(loc=0.015, scale=0.0075, size=real_data[col].shape)
    fake_col = real_data[col] * noise_dynamic + noise_static
    fake_col[fake_col < 1e-10] = 1e-10  # no negative intensity
    fake_data.append(fake_col)
fake_data.append(real_data['time'])
fake_data = pd.DataFrame(fake_data, index=real_data.columns).T

# visualize
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
for col in fake_data.columns[:-1]:
    ax.plot(time_post, real_data[col], color="tab:gray")
    ax.scatter(time_post, fake_data[col], label=col)
ax.plot(np.nan, np.nan, color="tab:gray", label="real")
ax.set_xlabel("time")
ax.set_ylabel("intensity")
ax.legend()
fig.show()

The RateConstantOptimizerTemplate contains several functions which allow it to optimize a model. However, it does require the user to define two functions. The first create_prediction should tell the class how exactly it can create a prediction from a given set of parameters. This function allows the user also to modify the parameters its given, or to modify the output of the prediction. The second function calculate_curves describes how the data should be analyzed.

from __future__ import annotations  # required for compatibility with python 3.8
from delayed_reactant_labeling.predict import DRL
from delayed_reactant_labeling.optimize import RateConstantOptimizerTemplate

class RateConstantOptimizer(RateConstantOptimizerTemplate):
    @staticmethod
    def create_prediction(x: np.ndarray, x_description: list[str]) -> pd.DataFrame:
        rate_constants = pd.Series(x, x_description)
        # The rate constants can easily be manipulated here. For example,
        # rate_constants["k1"] = 0.42 would fixate that value.
        # Because Series are mutable, the changed version will be stored in the logs!

        drl = DRL(reactions=reactions, rate_constants=rate_constants)
        pred_labeled = drl.predict_concentration(
            t_eval_pre=time_pre,
            t_eval_post=time_post,
            initial_concentrations=concentration_initial,
            labeled_concentration=concentration_labeled,
            dilution_factor=dilution_factor,
            rtol=1e-8,
            atol=1e-8, )

        # The prediction can be altered here before its analyzed.
        return pred_labeled

    @staticmethod
    def calculate_curves(data: pd.DataFrame) -> dict[str, np.ndarray]:
        curves = {}
        for chemical in ['A', 'B', 'C']:
            chemical_sum = data[[chemical, f'{chemical}-d10']].sum(axis=1)
            curves[f'ratio_{chemical}'] = (data[chemical] / chemical_sum).to_numpy()
        return curves

Internally the class compares the curves of the predicted data with the experimental data using a metric function. The function that calculates the Mean Absolute Error can for example be defined as follows:

def metric(y_true, y_pred):
    return np.average( np.abs(y_true - y_pred))

The np.average function also takes a weight keyword, so this can easily be implemented into the metric. Other functions such as np.nanmean can be used to skip NaN values. Scikit-learn.metrics implements a lot of different metrics and is a great resource. However, their functions also always check the arguments very precisely, which can lead to significant slow-downs.

To start optimizing this system we do the following:

from scipy.optimize import Bounds
description = ['k1', 'k-1', 'k2']
bounds = Bounds(np.array([1e-9, 1e-9, 1e-9]), np.array([100, 100, 100]))  #lower bound, upper bound

RCO = RateConstantOptimizer(experimental=fake_data, metric=metric)
RCO.optimize(
    x0=np.array([1, 1, 1]),
    x_description=description,
    x_bounds=bounds,
    path='./optimization/', _overwrite_log=True)

model = RCO.load_optimized_model('./optimization/')

model.best_X
# k1,   2.112804e-01
# k-1,  1.000000e-09
# k2,   6.424953e-01

However, the results of a optimization run are heavily dependent on its starting position, especially for systems with many different reactions. To analyze multiple runs from different starting positions we can use:

RCO.optimize_multiple(
    path='./optimization_multiple/',
    n_runs=50,
    x_description=description,
    x_bounds=bounds,
    n_jobs=-2,  # use all cpu's except one.
)

from delayed_reactant_labeling.optimize import OptimizedMultipleModels
models = OptimizedMultipleModels(path='./optimization_multiple/')
print(models.best.optimal_x)
# k1     2.112738e-01
# k-1    1.000000e-09
# k2     6.425392e-01