From aa0a36e7e4b75651c8fa40984e9390e1e113d59a Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Wed, 17 Apr 2024 22:42:38 +0900 Subject: [PATCH 1/9] Add MSE --- examples/model_comparison.py | 105 +++++++++++++---------------------- 1 file changed, 38 insertions(+), 67 deletions(-) diff --git a/examples/model_comparison.py b/examples/model_comparison.py index 21fae1a..edd66b6 100644 --- a/examples/model_comparison.py +++ b/examples/model_comparison.py @@ -4,97 +4,68 @@ import matplotlib.pyplot as plt import numpy as np +from scipy.ndimage.interpolation import zoom + + +def mse(x, y): + n = len(x) + assert n == len(y) + return 1 / n * sum((x - y) ** 2) def main(): - 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") - ] + path = Path("./examples/data/fad_kinetics") + start = 19 + raw_data = np.genfromtxt(path / "FADpH21_data.txt")[start:] + raw_data_time = np.genfromtxt(path / "FADpH21_data_time.txt")[start:] + kinetics_data = np.genfromtxt(path / "FAD_pH21_kinetics.txt") + kinetics_time = np.genfromtxt(path / "FAD_pH21_kinetics_time.txt") + semiclassical_data = np.genfromtxt(path / "semiclassical.txt") + semiclassical_time = np.genfromtxt(path / "semiclassical_time.txt") + semiclassical_kinetics_data = np.genfromtxt(path / "semiclassical_kinetics.txt") + semiclassical_kinetics_time = np.genfromtxt( + path / "semiclassical_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") - ] + print(f"{len(kinetics_data)=}") # 500 + print(f"{len(semiclassical_data)=}") # 600 + print(f"{len(semiclassical_kinetics_data)=}") # 600 + + raw_size = len(raw_data) + raw_kinetics_zoom = zoom(raw_data, len(kinetics_data) / raw_size) + raw_semiclassical_zoom = zoom(raw_data, len(semiclassical_data) / raw_size) + raw_semiclassical_kinetics_zoom = zoom( + raw_data, len(semiclassical_kinetics_data) / raw_size ) - # semiclassical_kinetics_time3 = np.array( - # [ - # np.genfromtxt(file_path) - # for file_path in Path(path).glob("semiclassical_kinetics_time3.txt") - # ] - # ) + print(f"{mse(kinetics_data, raw_kinetics_zoom)=}") + print(f"{mse(semiclassical_data, raw_semiclassical_zoom)=}") + print(f"{mse(semiclassical_kinetics_data, raw_semiclassical_kinetics_zoom)=}") plt.figure(1) plt.plot( - raw_data_time[0, :] - raw_data_time[0, 0], - raw_data[0, :] / raw_data[0, :].max(), + raw_data_time - raw_data_time[0], + raw_data / raw_data.max(), "k", linewidth=3, ) plt.plot( - kinetics_time[0, :] * 1e6, - kinetics_data[0, :] / kinetics_data[0, :].max(), + kinetics_time * 1e6, + kinetics_data / kinetics_data.max(), "b", linewidth=3, ) plt.plot( - semiclassical_time[0, :], - semiclassical_data[0, :] / semiclassical_data[0, :].max(), + semiclassical_time, + semiclassical_data / semiclassical_data.max(), "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 / semiclassical_kinetics_data.max(), "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]) From a0d56b1570142033e51330e618f239946527a27b Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Wed, 17 Apr 2024 22:58:00 +0900 Subject: [PATCH 2/9] Normalise --- examples/model_comparison.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/examples/model_comparison.py b/examples/model_comparison.py index edd66b6..4c465e7 100644 --- a/examples/model_comparison.py +++ b/examples/model_comparison.py @@ -18,12 +18,17 @@ def main(): path = Path("./examples/data/fad_kinetics") start = 19 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" ) @@ -42,27 +47,28 @@ def main(): print(f"{mse(semiclassical_kinetics_data, raw_semiclassical_kinetics_zoom)=}") plt.figure(1) + # plt.plot(raw_data_time) plt.plot( - raw_data_time - raw_data_time[0], - raw_data / raw_data.max(), + raw_data_time, + raw_data, "k", linewidth=3, ) plt.plot( kinetics_time * 1e6, - kinetics_data / kinetics_data.max(), + kinetics_data, "b", linewidth=3, ) plt.plot( semiclassical_time, - semiclassical_data / semiclassical_data.max(), + semiclassical_data, "g", linewidth=3, ) plt.plot( semiclassical_kinetics_time * 1e6, - semiclassical_kinetics_data / semiclassical_kinetics_data.max(), + semiclassical_kinetics_data, "r", linewidth=3, ) From 052d96db7ad3d94a79917fd0b7cec1e7df640571 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 18 Apr 2024 11:38:28 +0900 Subject: [PATCH 3/9] Fix MSE --- README.md | 2 ++ examples/model_comparison.py | 46 ++++++++++++++++-------------------- 2 files changed, 22 insertions(+), 26 deletions(-) diff --git a/README.md b/README.md index 5a8c435..d697b03 100644 --- a/README.md +++ b/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 diff --git a/examples/model_comparison.py b/examples/model_comparison.py index 4c465e7..c0c4a1f 100644 --- a/examples/model_comparison.py +++ b/examples/model_comparison.py @@ -7,16 +7,15 @@ from scipy.ndimage.interpolation import zoom -def mse(x, y): - n = len(x) - assert n == len(y) - return 1 / n * sum((x - y) ** 2) +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) -def main(): +def main(start=19): path = Path("./examples/data/fad_kinetics") - start = 19 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:] @@ -32,28 +31,16 @@ def main(): semiclassical_kinetics_time = np.genfromtxt( path / "semiclassical_kinetics_time.txt" ) - print(f"{len(kinetics_data)=}") # 500 - print(f"{len(semiclassical_data)=}") # 600 - print(f"{len(semiclassical_kinetics_data)=}") # 600 - - raw_size = len(raw_data) - raw_kinetics_zoom = zoom(raw_data, len(kinetics_data) / raw_size) - raw_semiclassical_zoom = zoom(raw_data, len(semiclassical_data) / raw_size) - raw_semiclassical_kinetics_zoom = zoom( - raw_data, len(semiclassical_kinetics_data) / raw_size - ) - print(f"{mse(kinetics_data, raw_kinetics_zoom)=}") - print(f"{mse(semiclassical_data, raw_semiclassical_zoom)=}") - print(f"{mse(semiclassical_kinetics_data, raw_semiclassical_kinetics_zoom)=}") + 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, raw_data, "k", linewidth=3) plt.plot( kinetics_time * 1e6, kinetics_data, @@ -86,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)=}") From c6842e7ec23a65a0042c0692ea1f56c712770d8c Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 18 Apr 2024 11:51:20 +0900 Subject: [PATCH 4/9] Fix docs --- .readthedocs.yaml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.readthedocs.yaml b/.readthedocs.yaml index b59a964..cde2cf4 100644 --- a/.readthedocs.yaml +++ b/.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 From ce3c1c36f9b855265f4bd3236bc0fabf2f597fb0 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 18 Apr 2024 11:54:36 +0900 Subject: [PATCH 5/9] Add `docs/requirements.txt` --- docs/requirements.txt | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/requirements.txt diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..6c5d5d4 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1 @@ +sphinx-rtd-theme From 5c59e0d50a035685535449cc4b3f24b9f2200b29 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 18 Apr 2024 11:56:44 +0900 Subject: [PATCH 6/9] Add typehints --- docs/requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/requirements.txt b/docs/requirements.txt index 6c5d5d4..2fae7cf 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1 +1,2 @@ sphinx-rtd-theme +sphinx-autodoc-typehints From 9ddbcc3d4b5e8efd963748c9a5ca665475b87f55 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 18 Apr 2024 12:01:00 +0900 Subject: [PATCH 7/9] Enable typehints --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 0e6335b..4d16483 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -28,7 +28,7 @@ "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", ] From 910c62317f08cde174dc75877b6bf41fa644fc72 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 18 Apr 2024 17:48:15 +0900 Subject: [PATCH 8/9] Rename parameter --- radicalpy/classical.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/radicalpy/classical.py b/radicalpy/classical.py index 86b7046..fb830a7 100644 --- a/radicalpy/classical.py +++ b/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]) From 87a05f2e12169faff566634533db6270a1b368d2 Mon Sep 17 00:00:00 2001 From: Emil VATAI Date: Thu, 18 Apr 2024 23:23:58 +0900 Subject: [PATCH 9/9] Intersphinx --- docs/conf.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/conf.py b/docs/conf.py index 4d16483..c69cfc9 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -24,6 +24,7 @@ "sphinx.ext.coverage", "sphinx.ext.githubpages", "sphinx.ext.ifconfig", + "sphinx.ext.intersphinx", "sphinx.ext.mathjax", "sphinx.ext.napoleon", "sphinx.ext.todo", @@ -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