added example file to do iterative optimization

llepe/__init__.py

@@ -1,2 +1,4 @@
 from .llepe import LLEPE
-from .utils import get_xml_value, set_size
+from .utils import get_xml_value, set_size
+from .optimizers import *
+from .objectives import *

llepe/objectives.py

@@ -0,0 +1,59 @@
+import numpy as np
+import pandas as pd
+
+
+def lmse_perturbed_obj(predicted_dict,
+                       measured_df,
+                       species_list,
+                       epsilon=1e-100):
+    meas_aq = np.concatenate([measured_df['{0}_aq_eq'.format(species)].values
+                              for species in species_list])
+    pred_aq = np.concatenate([
+        predicted_dict['{0}_aq_eq'.format(species)]
+        for species in species_list])
+
+    meas_d = np.concatenate([measured_df['{0}_d_eq'.format(species)].values
+                             for species in species_list])
+    pred_d = np.concatenate([
+        predicted_dict['{0}_d_eq'.format(species)]
+        for species in species_list])
+
+    meas_org = meas_aq * meas_d
+    pred_org = np.concatenate([
+        predicted_dict['{0}_org_eq'.format(species)]
+        for species in species_list])
+
+    perturbed_pred_d = (pred_org + epsilon) / (pred_aq + epsilon)
+    perturbed_meas_d = (meas_org + epsilon) / (meas_aq + epsilon)
+    log_pred_d = np.log10(perturbed_pred_d)
+    log_meas_d = np.log10(perturbed_meas_d)
+
+    fun4 = (log_meas_d - log_pred_d) ** 2
+    obj = np.mean(fun4)
+    return obj
+
+
+def ind_lmse_perturbed_obj(predicted_dict,
+                           measured_df,
+                           species_list,
+                           epsilon=1e-100):
+    pred_df = pd.DataFrame(predicted_dict)
+    objectives = []
+    for i in range(len(pred_df)):
+        pred_aq = pred_df[['{0}_aq_eq'.format(species)
+                           for species in species_list]].values[i]
+        pred_org = pred_df[['{0}_org_eq'.format(species)
+                            for species in species_list]].values[i]
+        meas_aq = measured_df[['{0}_aq_eq'.format(species)
+                               for species in species_list]].values[i]
+        meas_d = measured_df[['{0}_d_eq'.format(species)
+                              for species in species_list]].values[i]
+        meas_org = meas_aq * meas_d
+        perturbed_pred_d = (pred_org + epsilon) / (pred_aq + epsilon)
+        perturbed_meas_d = (meas_org + epsilon) / (meas_aq + epsilon)
+        log_pred_d = np.log10(perturbed_pred_d)
+        log_meas_d = np.log10(perturbed_meas_d)
+        fun1 = (log_meas_d - log_pred_d) ** 2
+        objectives.append(np.mean(fun1))
+    objectives = np.array(objectives)
+    return objectives

llepe/optimizers.py

@@ -0,0 +1,39 @@
+import scipy.optimize as scipy_opt
+from scipy.optimize import minimize
+import skopt
+
+def dual_anneal_optimizer(objective, x_guess):
+    bounds = [(1e-1, 1e1)] * len(x_guess)
+    bounds[1] = (1e-1, 2)
+    res = scipy_opt.dual_annealing(objective,
+                                   [(1e-1, 1e1)]*len(x_guess),
+                                   x0=x_guess)
+    est_parameters = res.x
+    return est_parameters, res.fun
+
+
+def diff_evo_optimizer(objective, x_guess):
+    bounds = [(1e-1, 1e1)] * len(x_guess)
+    bounds[1] = (1e-1, 2)
+    res = scipy_opt.differential_evolution(objective,
+                                           bounds)
+    est_parameters = res.x
+    return est_parameters, res.fun
+
+
+def forest_lbfgsb_optimizer(objective, x_guess):
+    x_guess = list(x_guess)
+    bounds = [(1e-1, 1e1)]*len(x_guess)
+    bounds[1] = (1e-1, 2)
+    res = skopt.forest_minimize(objective,
+                                bounds,
+                                random_state=1,
+                                acq_func='LCB',
+                                n_random_starts=30,
+                                x0=x_guess,
+                                xi=1e-4)
+    x_guess = res.x
+    optimizer_kwargs = {"method": 'l-bfgs-b',
+                        "bounds": bounds}
+    res = minimize(objective, x_guess, **optimizer_kwargs)
+    return res.x, res.fun