Merge pull request #156 from Spin-Chemistry-Labs/155-add-mse

Add MSE

.readthedocs.yaml

@@ -6,7 +6,7 @@ version: 2
 
 # Set the version of Python and other tools you might need
 build:
-  os: ubuntu-20.04
+  os: ubuntu-22.04
   tools:
     python: "3.10"
     # You can also specify other tool versions:
@@ -25,4 +25,5 @@ sphinx:
 # Optionally declare the Python requirements required to build your docs
 python:
    install:
+   - requirements: docs/requirements.txt
    - requirements: requirements.txt

README.md

@@ -1,3 +1,5 @@
+[![Documentation Status](https://readthedocs.org/projects/radicalpy/badge/?version=latest)](https://radicalpy.readthedocs.io/en/latest/?badge=latest)
+
 # RadicalPy: a toolbox for radical pair spin dynamics
 
 RadicalPy in an intuitive (object-oriented) open-source Python

docs/conf.py

@@ -24,11 +24,12 @@
     "sphinx.ext.coverage",
     "sphinx.ext.githubpages",
     "sphinx.ext.ifconfig",
+    "sphinx.ext.intersphinx",
     "sphinx.ext.mathjax",
     "sphinx.ext.napoleon",
     "sphinx.ext.todo",
     "sphinx.ext.viewcode",
-    # "sphinx_autodoc_typehints", # did't work when I tried :(
+    "sphinx_autodoc_typehints",  # did't work when I tried :(
     "sphinx_rtd_theme",
 ]
 
@@ -37,6 +38,10 @@
 exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
 # autoclass_content = "both"
 todo_include_todos = True
+intersphinx_mapping = {
+    "python": ("https://docs.python.org/3", None),
+    "numpy": ("http://docs.scipy.org/doc/numpy", None),
+}
 
 # -- Options for HTML output -------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output

docs/requirements.txt

@@ -0,0 +1,2 @@
+sphinx-rtd-theme
+sphinx-autodoc-typehints

examples/model_comparison.py

@@ -4,97 +4,61 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+from scipy.ndimage.interpolation import zoom
 
 
-def main():
+def mse(raw, raw_time, data):
+    N = len(data) - 1
+    indices = N * raw_time / raw_time[-1]
+    return sum((data[indices.astype(int)] - raw) ** 2) / len(raw)
 
-    path = "./examples/data/fad_kinetics"
-    raw_data = np.array(
-        [np.genfromtxt(file_path) for file_path in Path(path).glob("FADpH21_data.txt")]
-    )
-    raw_data_time = np.array(
-        [
-            np.genfromtxt(file_path)
-            for file_path in Path(path).glob("FADpH21_data_time.txt")
-        ]
-    )
-    kinetics_data = np.array(
-        [
-            np.genfromtxt(file_path)
-            for file_path in Path(path).glob("FAD_pH21_kinetics.txt")
-        ]
-    )
-    kinetics_time = np.array(
-        [
-            np.genfromtxt(file_path)
-            for file_path in Path(path).glob("FAD_pH21_kinetics_time.txt")
-        ]
-    )
-    semiclassical_data = np.array(
-        [np.genfromtxt(file_path) for file_path in Path(path).glob("semiclassical.txt")]
-    )
-    semiclassical_time = np.array(
-        [
-            np.genfromtxt(file_path)
-            for file_path in Path(path).glob("semiclassical_time.txt")
-        ]
-    )
-    semiclassical_kinetics_data = np.array(
-        [
-            np.genfromtxt(file_path)
-            for file_path in Path(path).glob("semiclassical_kinetics.txt")
-        ]
-    )
-    # semiclassical_kinetics3 = np.array(
-    #     [
-    #         np.genfromtxt(file_path)
-    #         for file_path in Path(path).glob("semiclassical_kinetics_new3.txt")
-    #     ]
-    # )
-    semiclassical_kinetics_time = np.array(
-        [
-            np.genfromtxt(file_path)
-            for file_path in Path(path).glob("semiclassical_kinetics_time.txt")
-        ]
+
+def main(start=19):
+
+    path = Path("./examples/data/fad_kinetics")
+    raw_data = np.genfromtxt(path / "FADpH21_data.txt")[start:]
+    raw_data /= raw_data.max()
+    raw_data_time = np.genfromtxt(path / "FADpH21_data_time.txt")[start:]
+    raw_data_time -= raw_data_time[0]
+    kinetics_data = np.genfromtxt(path / "FAD_pH21_kinetics.txt")
+    kinetics_data /= kinetics_data.max()
+    kinetics_time = np.genfromtxt(path / "FAD_pH21_kinetics_time.txt")
+    semiclassical_data = np.genfromtxt(path / "semiclassical.txt")
+    semiclassical_data /= semiclassical_data.max()
+    semiclassical_time = np.genfromtxt(path / "semiclassical_time.txt")
+    semiclassical_kinetics_data = np.genfromtxt(path / "semiclassical_kinetics.txt")
+    semiclassical_kinetics_data /= semiclassical_kinetics_data.max()
+    semiclassical_kinetics_time = np.genfromtxt(
+        path / "semiclassical_kinetics_time.txt"
     )
-    # semiclassical_kinetics_time3 = np.array(
-    #     [
-    #         np.genfromtxt(file_path)
-    #         for file_path in Path(path).glob("semiclassical_kinetics_time3.txt")
-    #     ]
-    # )
+    mse_semiclassical = mse(raw_data, raw_data_time, semiclassical_data)
+    mse_kinetics = mse(raw_data, raw_data_time, kinetics_data)
+    mse_new_method = mse(raw_data, raw_data_time, semiclassical_kinetics_data)
+    print(f"{mse_semiclassical=}")
+    print(f"{mse_kinetics=}")
+    print(f"{mse_new_method=}")
 
     plt.figure(1)
+    # plt.plot(raw_data_time)
+    plt.plot(raw_data_time, raw_data, "k", linewidth=3)
     plt.plot(
-        raw_data_time[0, :] - raw_data_time[0, 0],
-        raw_data[0, :] / raw_data[0, :].max(),
-        "k",
-        linewidth=3,
-    )
-    plt.plot(
-        kinetics_time[0, :] * 1e6,
-        kinetics_data[0, :] / kinetics_data[0, :].max(),
+        kinetics_time * 1e6,
+        kinetics_data,
         "b",
         linewidth=3,
     )
     plt.plot(
-        semiclassical_time[0, :],
-        semiclassical_data[0, :] / semiclassical_data[0, :].max(),
+        semiclassical_time,
+        semiclassical_data,
         "g",
         linewidth=3,
     )
     plt.plot(
-        semiclassical_kinetics_time[0, :] * 1e6,
-        semiclassical_kinetics_data[0, :] / semiclassical_kinetics_data[0, :].max(),
+        semiclassical_kinetics_time * 1e6,
+        semiclassical_kinetics_data,
         "r",
         linewidth=3,
     )
-    # plt.plot(
-    #     semiclassical_kinetics_time3[0, :] * 1e6,
-    #     semiclassical_kinetics3[0, :, 16] / semiclassical_kinetics3[0, :, 16].max(),
-    #     "m",
-    #     linewidth=3,
-    # )
     plt.xlabel("Time / $\mu s$", size=24)
     plt.ylabel("Normalised $\Delta \Delta A$ / a.u.", size=24)
     plt.ylim([-0.1, 1.1])
@@ -109,7 +73,14 @@ def main():
 
     path = __file__[:-3] + f"_{0}.png"
     plt.savefig(path, dpi=300, bbox_inches="tight")
+    return mse_semiclassical, mse_kinetics, mse_new_method
 
 
 if __name__ == "__main__":
-    main()
+    results = []
+    for start in [19]:
+        # for start in range(0, 30):
+        print(f"raw data cut off: {start=}")
+        results.append(main(start))
+    results = np.array(results)
+    print(f"{results.min(axis=0)=}")

radicalpy/classical.py

@@ -169,19 +169,18 @@ def reaction_scheme(path: str, rate_equations: dict):
     Path(path).write_text(texcode)
 
 
-def random_theta_phi(n: int = 1) -> ArrayLike:
+def random_theta_phi(num_samples: int = 1) -> ArrayLike:
     """Random sampling of theta and phi.
 
     Args:
 
-            n_samples (int): The number of samples generated.
+            num_samples: The number of samples generated.
 
     Returns:
             Theta and phi (radians).
-
     """
-    phi = np.random.uniform(0, 2 * np.pi, size=n)
-    theta = np.arccos(np.random.uniform(-1, 1, size=n))
+    phi = np.random.uniform(0, 2 * np.pi, size=num_samples)
+    theta = np.arccos(np.random.uniform(-1, 1, size=num_samples))
     return np.array([theta, phi])