Ready to fit models! ✨

Author: echedey-ls

Irradiances_ratios/E_ratio_script.py

@@ -11,8 +11,7 @@
 
 import numpy as np
 import pandas as pd
-
-from datetime import datetime
+# import matplotlib.pyplot as plt
 
 # Matrix of values to test
 # Atmosphere characterization required params
@@ -27,11 +26,21 @@
 }
 
 # what do we want to plot E_λ<λ₀/E against? (None = default behaviour)
-plot_keys = "datetime"
+plot_keys = None
+
+# model function
+model_inputs = ["relative_airmass", "aerosol_turbidity_500nm"]
+p0 = None  # [0.2, 0.1]
+
+
+def model(xdata, c0, c1):  # use this func as model template
+    r_am, aod500 = xdata
+    return c0 * r_am + c1 * aod500
+
 
 bench = MR_E_ratio(
     datetimes=pd.date_range(
-        "2023-11-27T00", "2023-11-28T00", freq=pd.Timedelta(minutes=1)
+        "2023-11-27T00", "2023-11-28T00", freq=pd.Timedelta(minutes=15)
     )
 )
 
@@ -40,14 +49,20 @@
 bench.cutoff_lambda = LAMBDA0["monosi"]  # == polysi
 bench.simulate_from_product(**spectrl2_generator_input)
 bench.plot_results(plot_keys=plot_keys)
+optim_result = bench.optimization_from_model(
+    model=model, model_inputs=model_inputs, p0=p0
+)
 bench.times_summary()
 
+# TODO: PLOT RESULTS & MODEL PREDICTION
+
 # %%
 # Test with asi cutoff wavelength
 bench.reset_simulation_state()
 bench.cutoff_lambda = LAMBDA0["asi"]
 bench.simulate_from_product(**spectrl2_generator_input)
 bench.plot_results(plot_keys=plot_keys)
+bench.optimization_from_model(model=model, model_inputs=model_inputs, p0=p0)
 bench.times_summary()
 
 # %%

Irradiances_ratios/spectrl2_E_ratio_bench.py

@@ -15,11 +15,13 @@
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
 
 from itertools import product
 from functools import partial
 from datetime import datetime
 from time import time
+from typing import Callable
 
 
 class MR_E_ratio:
@@ -118,11 +120,11 @@ def simulation_prerun(self):
             self.solpos["azimuth"],
         )
         self.time_params = {
-            "apparent_zenith": self.solpos["apparent_zenith"],
-            "aoi": self.aoi,
+            "apparent_zenith": self.solpos["apparent_zenith"].to_numpy(),
+            "aoi": self.aoi.to_numpy(),
             "relative_airmass": self.locus.get_airmass(solar_position=self.solpos)[
                 "airmass_relative"
-            ],
+            ].to_numpy(),
             "dayofyear": np.fromiter(
                 map(day_of_year, self.datetimes), dtype=np.float64
             ),
@@ -306,3 +308,72 @@ def plot_results(
         plt.close()
 
         self.processing_time["plot_results"] = time() - start_time
+
+    def optimization_from_model(
+        self, model: Callable = None, model_inputs: tuple = None, **kwargs
+    ):
+        """
+        Optimize a model to fit generated data.
+
+        Parameters
+        ----------
+        model : Callable
+            Function with the model to be optimised.
+        model_inputs : str or iterable of str
+            Order and parameters of ``model``. Must be any of:
+                * ``datetime``
+                * ``apparent_zenith``, ``aoi``, ``relative_airmass`` or ``dayofyear``
+                * any parameter name provided to ``simulate_from_product``
+        **kwargs :
+            Redirected to ``scipy.optimize.curve_fit``.
+
+        Returns
+        -------
+        ``scipy.optimize.curve_fit``'s return values
+        """
+        start_time = time()  # Initialize start time of block
+
+        if isinstance(model_inputs, str):
+            model_inputs = (model_inputs,)
+
+        # number of inputs from user: n-left-most columns
+        n_inputs = len(self.input_keys)
+        # Fitting data
+        ydata = self.results.iloc[:, n_inputs:].to_numpy().flatten()
+        # Prepare input vector
+        xdata = []  # length of each value must be
+        dates_len = len(self.datetimes)
+        try:
+            for var_name in model_inputs:
+                # broadcast all inputs to match ydata
+                if var_name in self.input_keys:
+                    xdata.append(self.results[var_name].to_numpy().repeat(dates_len))
+                elif var_name in self.time_params.keys():
+                    xdata.append(
+                        np.tile(self.time_params[var_name], self.results.shape[0])
+                    )
+                elif var_name in {"datetime"}:
+                    xdata.append(
+                        np.tile(self.datetimes.to_numpy(), self.results.shape[0])
+                    )
+                else:
+                    raise ValueError(f"'{var_name}' is not a valid parameter name!")
+
+        except TypeError:
+            raise TypeError(
+                "Provide a valid model input names vector. Must be iterable"
+                + " of strings, and that input will be provided to 'model'"
+                + f" in the same order.\nYou provided {model_inputs}"
+            )
+
+        ## This is kept here for debug purposes: check valid representation as 1D arrays
+        # fig, axs = plt.subplots(len(model_inputs))
+        # for i, name in enumerate(model_inputs):
+        #     axs[i].set_title(name)
+        #     axs[i].scatter(xdata[i], ydata)
+        # plt.show()
+
+        curve_fit_results = curve_fit(model, xdata, ydata, nan_policy="omit", **kwargs)
+
+        self.processing_time["optimization_from_model"] = time() - start_time
+        return curve_fit_results