✨Add basic example to evaluate doas feature (#44)

* Basic DOAS example (beta)

Add a basic example to test doas feature branch, see PR pyglotaran#764: glotaran/pyglotaran#764

* Test doas example on CI

* 🩹 Added doas data  and made an exception in the .gitignore

* 🧹 Removed old style 'type' in the model

* 👌 Added coherent artifact for testing

'vary: False' and 'maximum_number_function_evaluations=1' are only set for performance reasons

* ✨Added data file with more values

Co-authored-by: s-weigand <[email protected]>

.gitignore

@@ -135,3 +135,6 @@ dmypy.json
 
 # plot results when run headless
 plot_results
+
+# DOAS data file
+!pyglotaran_examples/ex_doas_beta/data/2008Polli_betacar_chex_sim.nc

.pre-commit-config.yaml

@@ -39,4 +39,4 @@ repos:
     hooks:
       - id: codespell
         files: ".py|.rst"
-        args: [--skip="*.ipynb"]
+        args: ["--skip='*.ipynb'", "-L doas,DOAS"]

pyglotaran_examples/ex_doas_beta/ex_doas_beta.py

@@ -0,0 +1,80 @@
+# To use this script as a 'notebook' in VS Code add '# %%'
+
+# %% Imports
+from datetime import datetime
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from glotaran.analysis.optimize import optimize
+from glotaran.io import load_dataset
+from glotaran.io import load_model
+from glotaran.io import load_parameters
+from glotaran.io import save_result
+from glotaran.project.scheme import Scheme
+from pyglotaran_extras.io.boilerplate import setup_case_study
+from pyglotaran_extras.plotting.plot_overview import plot_overview
+from pyglotaran_extras.plotting.style import PlotStyle
+
+DATA_PATH = "data/2008Polli_betacar_chex_sim.nc"
+MODEL_PATH = "models/model.yml"
+PARAMETERS_FILE_PATH = "models/parameters.yml"
+
+
+# %% Define function
+def run_doas_model(show_plot=False, block_plot=False):
+
+    results_folder, script_folder = setup_case_study(
+        output_folder_name="pyglotaran_examples_results"
+    )
+    results_folder = results_folder / "target_analysis"
+    print(f"- Using folder {results_folder} to read/write files for this run")
+
+    dataset = load_dataset(script_folder.joinpath(DATA_PATH))
+    model = load_model(script_folder.joinpath(MODEL_PATH))
+    parameters = load_parameters(script_folder.joinpath(PARAMETERS_FILE_PATH))
+
+    print(f"\n{'#'*10} DOAS Model {'#'*10}\n")
+    print(model.validate(parameters=parameters))
+
+    scheme = Scheme(
+        model=model,
+        parameters=parameters,
+        data={"dataset1": dataset},
+        non_negative_least_squares=False,
+        optimization_method="TrustRegionReflection",
+        # maximum_number_function_evaluations=3,
+        maximum_number_function_evaluations=1,
+    )
+    result = optimize(scheme)
+
+    print(f"\n{'#'*3} DOAS Model - Optimization Result {'#'*3}\n")
+    print(result)
+    print(f"\n{'#'*3} DOAS Model - Optimized Parameters {'#'*3}\n")
+    print(result.optimized_parameters)
+
+    save_result(result, results_folder, format_name="legacy", allow_overwrite=True)
+
+    plot_style = PlotStyle()
+    plt.rc("axes", prop_cycle=plot_style.cycler)
+
+    if show_plot:
+        plot_style = PlotStyle()
+        plt.rc("axes", prop_cycle=plot_style.cycler)
+
+        fig = plot_overview(result.data["dataset1"], linlog=True)
+        plt.rcParams["figure.figsize"] = (21, 14)
+
+        timestamp = datetime.today().strftime("%y%m%d_%H%M")
+        fig.savefig(results_folder.joinpath(f"plot_overview_{timestamp}.pdf"), bbox_inches="tight")
+
+        plt.show(block=block_plot)
+    script_dir = Path(__file__).resolve().parent
+    print(f"Script folder: {script_dir}")
+    script_dir.cwd()
+
+
+# %% Main
+if __name__ == "__main__":
+    run_doas_model(show_plot=True, block_plot=True)
+else:
+    run_doas_model(show_plot=True, block_plot=False)

pyglotaran_examples/ex_doas_beta/models/model.yml

@@ -0,0 +1,99 @@
+# type: doas
+
+initial_concentration:
+  j1:
+    compartments: [S2FC, S2, hotS1, S1]
+    parameters: [j.1, j.0, j.0, j.0]
+
+k_matrix:
+    k1:
+      matrix:
+        '(S2, S2FC)': kinetic.1
+        '(hotS1, S2)': kinetic.2
+        '(S1, hotS1)': kinetic.3
+        '(S1, S1)': kinetic.4
+
+megacomplex:
+  doas:
+    type: damped-oscillation
+    labels: [
+      osc1,
+      osc2,
+      osc3,
+      osc4,
+      osc5,
+      osc6,
+      osc7,
+      osc8,
+      osc9,
+      osc10,
+      osc11,
+      osc12,
+      osc13,
+      osc14,
+      osc15,
+      osc16,
+      osc17,
+    ]
+    frequencies: [
+      osc.freq.1,
+      osc.freq.2,
+      osc.freq.3,
+      osc.freq.4,
+      osc.freq.5,
+      osc.freq.6,
+      osc.freq.7,
+      osc.freq.8,
+      osc.freq.9,
+      osc.freq.10,
+      osc.freq.11,
+      osc.freq.12,
+      osc.freq.13,
+      osc.freq.14,
+      osc.freq.15,
+      osc.freq.16,
+      osc.freq.17,
+    ]
+    rates: [
+      osc.rate.1,
+      osc.rate.2,
+      osc.rate.3,
+      osc.rate.4,
+      osc.rate.5,
+      osc.rate.6,
+      osc.rate.7,
+      osc.rate.8,
+      osc.rate.9,
+      osc.rate.10,
+      osc.rate.11,
+      osc.rate.12,
+      osc.rate.13,
+      osc.rate.14,
+      osc.rate.15,
+      osc.rate.16,
+      osc.rate.17,
+    ]
+  decay:
+    type: decay
+    k_matrix: [k1]
+  baseline:
+    type: baseline
+    dimension: time
+  artifact:
+    type: coherent-artifact
+    order: 3
+
+
+irf:
+  irf1:
+    type: spectral-gaussian
+    center: irf.center
+    width: irf.width
+    dispersion_center: irf.dispcenter
+    center_dispersion: [irf.disp1, irf.disp2, irf.disp3]
+
+dataset:
+  dataset1:
+    initial_concentration: j1
+    megacomplex: [doas, decay, artifact, baseline]
+    irf: irf1

pyglotaran_examples/ex_doas_beta/models/parameters.yml

@@ -0,0 +1,30 @@
+osc:
+  freq: [
+1523.39905,       1178.15955 ,     1158.34766  , 1002.55316 , 967.485352  , 804.297241  ,721.601562  ,502.533569,
+395.366425    ,  288.248444    ,   242.915039       ,94.8479767     , 35.0097275  , 698.158386 , 2142.97974  , 2221.10693 , 3981.01416
+  ,
+  {vary: False}]
+  rate: [
+  1.54323947,  2.81027126,   1.07498121 , 1.65556943  ,5.93284607,
+1.19999278    ,  0.204981461  ,    0.886841118  ,    0.801028907  ,   1.83283675  , 1.09177792 , 0.317926913  ,2.02648878,  49.9548454  ,98.0651779,
+35.9500198 ,     65.5325775 ,
+  {vary: False}]
+
+j:
+  - ['1', 1.0, {non-negative: False, vary: False}]
+  - ['0', 0.0, {non-negative: False, vary: False}]
+
+kinetic: [
+  50.0,
+  5.9171597633136095,
+  2.808988764044944,
+  0.10905125408942203,
+]
+
+irf:
+  - ['center', -2.35207728E-03, {non-negative: False, vary: True}]
+  - ['width', 1.12892091E-02, {non-negative: True, vary: True}]
+  - ['dispcenter', 600.0, {non-negative: False, vary: False}]
+  - ['disp1', 2.75756628E-03, {non-negative: False, vary: True}]
+  - ['disp2', -2.39017978E-02, {non-negative: False, vary: True}]
+  - ['disp3', -3.09822452E-03, {non-negative: False, vary: True}]

pyglotaran_examples/ex_doas_beta/scheme.yml

@@ -0,0 +1,4 @@
+model: model.yml
+parameters: parameters.yml
+data:
+  dataset1: 2008Polli_betacar_chex_sim.nc

scripts/run_examples.py

@@ -159,6 +159,14 @@ def sim_6d_disp(*, headless=False, raise_on_deprecation=False):
     from pyglotaran_examples.test.simultaneous_analysis_6d_disp import sim_analysis_script_6d_disp
 
 
+@script_run_wrapper
+def doas_beta(*, headless=False, raise_on_deprecation=False):
+    """Runs ex_doas_beta.py
+    from pyglotaran_examples/ex_doas_beta"""
+    # The whole script is run at import.
+    from pyglotaran_examples.ex_doas_beta import ex_doas_beta
+
+
 all_funcs = [
     quick_start,
     fluorescence,
@@ -171,6 +179,7 @@ def sim_6d_disp(*, headless=False, raise_on_deprecation=False):
     sim_3d_nodisp,
     sim_3d_weight,
     sim_6d_disp,
+    doas_beta,
 ]