diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index bd198bac8..f1e30e5bf 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -72,7 +72,7 @@ jobs: run: pytest --cov=./ --cov-report=xml - name: Upload coverage - uses: codecov/codecov-action@v2.0.2 + uses: codecov/codecov-action@v1 # build_windows_test: # needs: [code_lint] @@ -112,7 +112,7 @@ jobs: # run: pytest --cov=./ --cov-report=xml # - name: Upload coverage - # uses: codecov/codecov-action@v2.0.2 + # uses: codecov/codecov-action@v1 testing_windows: needs: [code_lint] @@ -148,4 +148,4 @@ jobs: shell: pwsh run: pytest --cov=./ --cov-report=xml - name: Upload coverage - uses: codecov/codecov-action@v2.0.2 + uses: codecov/codecov-action@v1 diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index e4766faac..bc7105d13 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -33,7 +33,7 @@ jobs: run: pytest --cov=./ --cov-report=xml - name: Upload coverage - uses: codecov/codecov-action@v2.0.2 + uses: codecov/codecov-action@v1 source-distribution: needs: [build-test] diff --git a/.gitignore b/.gitignore index cc3f58d94..4fcf98c94 100644 --- a/.gitignore +++ b/.gitignore @@ -11,7 +11,6 @@ __pycache__/* .Python build/ .vscode/* -.idea/* develop-eggs/ dist/ downloads/ @@ -106,7 +105,6 @@ tests/simpson/quad/** version/ docs/examples/ docs/fitting/ -_temp/ # files class-diagram* diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6f4706029..caedcb08a 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -48,17 +48,6 @@ repos: hooks: - id: flake8 language: python - args: - [ - "--ignore=E402", - "--exclude=.eggs, *.egg,build, src/mrsimulator/__init__.py", - "--filename=*.pyx, *py", - "--max-line-length=88", - "--max-complexity=8", - "--select=B,C,E,F,W,T,N8", - "--count", - "--statistics", - ] - repo: https://github.com/asottile/reorder_python_imports rev: v2.4.0 diff --git a/.pylintrc b/.pylintrc new file mode 100644 index 000000000..83b6f246b --- /dev/null +++ b/.pylintrc @@ -0,0 +1 @@ +issue code: redefined-builtin diff --git a/CHANGELOG b/CHANGELOG index d7aede704..9c27e499b 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,23 +1,6 @@ Upcoming Changes ---------------- -What's new -'''''''''' - -**Method** - -- New Event class---``SpectralEvent``, ``ConstantTimeEvent``, and ``MixingEvent``. The MixingEvent - controls the mixing of transitions in a multi-event method. -- A ``TransitionPathway`` object now includes a weight attribute, which holds the probability of - the transition pathway based on the mixing events defined within the method. - -**SpinSystem** - -- Add a new function ``simplify()`` to the spin system object, which reduces a spin system to - a list of irreducible spin systems. -- Add a new function ``site_generator()`` to the utility collection sub-module, which simplify - the process of creating Site objects in bulk. - Changes ''''''' diff --git a/docs/api_py/method.rst b/docs/api_py/method.rst index 2a7085a9c..cf6e4bdf1 100644 --- a/docs/api_py/method.rst +++ b/docs/api_py/method.rst @@ -1,13 +1,10 @@ Method ====== -Method is the root user-level object that may be used in creating custom NMR simulation -methods. - .. toctree:: method/method method/spectral_dimension method/event method/freq_contrib - method/query + method/transition_query diff --git a/docs/api_py/method/event.rst b/docs/api_py/method/event.rst index 2fb8ae1ac..280a9be04 100644 --- a/docs/api_py/method/event.rst +++ b/docs/api_py/method/event.rst @@ -15,7 +15,3 @@ Events :show-inheritance: :members: :inherited-members: BaseModel - -.. autoclass:: MixingEvent - :show-inheritance: - :members: diff --git a/docs/api_py/method/query.rst b/docs/api_py/method/transition_query.rst similarity index 93% rename from docs/api_py/method/query.rst rename to docs/api_py/method/transition_query.rst index 0b87aa1a5..47f8ac6c1 100644 --- a/docs/api_py/method/query.rst +++ b/docs/api_py/method/transition_query.rst @@ -1,3 +1,4 @@ +.. _transition_query_api: .. currentmodule:: mrsimulator.method.query diff --git a/docs/api_py/methods-summary.rst b/docs/api_py/methods-summary.rst index ba09e95ab..ca27ff1d5 100644 --- a/docs/api_py/methods-summary.rst +++ b/docs/api_py/methods-summary.rst @@ -48,11 +48,6 @@ Summary .. ~Cosy -UML Diagram ------------ - -.. figure:: ../_static/classes_methods.* - Table of contents ----------------- diff --git a/docs/api_py/utils.rst b/docs/api_py/utils.rst index 00725d8de..f615beafc 100644 --- a/docs/api_py/utils.rst +++ b/docs/api_py/utils.rst @@ -7,7 +7,6 @@ Utility functions .. currentmodule:: mrsimulator.utils.collection .. autofunction:: single_site_system_generator -.. autofunction:: site_generator .. currentmodule:: mrsimulator.utils diff --git a/docs/getting_started-objects.rst b/docs/getting_started-objects.rst index c29c2df27..13907f6a3 100644 --- a/docs/getting_started-objects.rst +++ b/docs/getting_started-objects.rst @@ -468,11 +468,7 @@ see the list of transition pathways, for example, >>> from pprint import pprint >>> pprint(sim.methods[0].get_transition_pathways(system_4)) # 17O - [|-2.5⟩⟨-1.5|, weight=(1+0j), - |-1.5⟩⟨-0.5|, weight=(1+0j), - |-0.5⟩⟨0.5|, weight=(1+0j), - |0.5⟩⟨1.5|, weight=(1+0j), - |1.5⟩⟨2.5|, weight=(1+0j)] + [|-2.5⟩⟨-1.5|, |-1.5⟩⟨-0.5|, |-0.5⟩⟨0.5|, |0.5⟩⟨1.5|, |1.5⟩⟨2.5|] Notice, there are five transition pathways for the :math:`^{17}\text{O}` site, one associated with the central-transition, two with the inner-satellites, and two with @@ -497,7 +493,7 @@ the outer-satellites. For central transition selective simulation, use the ... ) >>> # the transition pathways >>> print(sim.methods[0].get_transition_pathways(system_4)) # 17O - [|-0.5⟩⟨0.5|, weight=(1+0j)] + [|-0.5⟩⟨0.5|] Now, you may simulate the central transition selective spectrum. :numref:`fig8_using_obj` depicts a central transition selective spectrum. diff --git a/docs/installation/requirements.rst b/docs/installation/requirements.rst index 3334c649f..5f3e93d1d 100644 --- a/docs/installation/requirements.rst +++ b/docs/installation/requirements.rst @@ -13,7 +13,7 @@ Package dependencies - typing-extensions>=3.7 - `matplotlib>=3.3.3 `_ for figures and visualization, - monty>=2.0.4 -- `csdmpy>=0.4.1 `_ +- `csdmpy>=0.4 `_ - `pydantic>=1.0 `_ - monty>=2.0.4 diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_0_Wollastonite.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_0_Wollastonite.ipynb index 39addf550..de0c11073 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_0_Wollastonite.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_0_Wollastonite.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "method = BlochDecaySpectrum(\n channels=[\"29Si\"],\n magnetic_flux_density=14.1, # in T\n rotor_frequency=1500, # in Hz\n spectral_dimensions=[\n {\n \"count\": 2048,\n \"spectral_width\": 25000, # in Hz\n \"reference_offset\": -10000, # in Hz\n \"label\": r\"$^{29}$Si resonances\",\n }\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(4, 2))\nmethod.plot()\nplt.show()" + "method = BlochDecaySpectrum(\n channels=[\"29Si\"],\n magnetic_flux_density=14.1, # in T\n rotor_frequency=1500, # in Hz\n spectral_dimensions=[\n {\n \"count\": 2048,\n \"spectral_width\": 25000, # in Hz\n \"reference_offset\": -10000, # in Hz\n \"label\": r\"$^{29}$Si resonances\",\n }\n ],\n)" ] }, { @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_1_PotassiumSulfate.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_1_PotassiumSulfate.ipynb index d19024c5e..2999ae5c0 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_1_PotassiumSulfate.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_1_PotassiumSulfate.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "method = BlochDecayCTSpectrum(\n channels=[\"33S\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=14000, # in Hz\n spectral_dimensions=[\n {\n \"count\": 2048,\n \"spectral_width\": 5000, # in Hz\n \"reference_offset\": 22500, # in Hz\n \"label\": r\"$^{33}$S resonances\",\n }\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(4, 3))\nmethod.plot()\nplt.show()" + "method = BlochDecayCTSpectrum(\n channels=[\"33S\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=14000, # in Hz\n spectral_dimensions=[\n {\n \"count\": 2048,\n \"spectral_width\": 5000, # in Hz\n \"reference_offset\": 22500, # in Hz\n \"label\": r\"$^{33}$S resonances\",\n }\n ],\n)" ] }, { @@ -150,7 +150,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_2_Coesite.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_2_Coesite.ipynb index 08daab795..21322b1ab 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_2_Coesite.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_2_Coesite.ipynb @@ -175,7 +175,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_quad_csa.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_quad_csa.ipynb index 55ffd8bb4..ab96fff45 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_quad_csa.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_quad_csa.ipynb @@ -132,7 +132,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_satellite_transition_sim.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_satellite_transition_sim.ipynb index 2262cf908..074601747 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_satellite_transition_sim.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_3_satellite_transition_sim.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "method = Method1D(\n name=\"Inner Satellite Spectrum\",\n channels=[\"27Al\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=1e9, # in Hz\n spectral_dimensions=[\n {\n \"count\": 1024,\n \"spectral_width\": 1e4, # in Hz\n \"reference_offset\": 1e4, # in Hz\n \"events\": [\n {\n \"transition_query\": [\n {\"P\": [-1], \"D\": [2]} # <-- select inner satellite transitions\n ]\n }\n ],\n }\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3))\nmethod.plot()\nplt.show()" + "method = Method1D(\n channels=[\"27Al\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=1e9, # in Hz\n spectral_dimensions=[\n {\n \"count\": 1024,\n \"spectral_width\": 1e4, # in Hz\n \"reference_offset\": 1e4, # in Hz\n \"events\": [\n {\n \"transition_query\": [\n {\"P\": [-1], \"D\": [2]} # <-- select inner satellite transitions\n ]\n }\n ],\n }\n ],\n)" ] }, { @@ -123,7 +123,7 @@ }, "outputs": [], "source": [ - "method2 = Method1D(\n name=\"Satellite Spectrum\",\n channels=[\"27Al\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=1e9, # in Hz\n spectral_dimensions=[\n {\n \"count\": 1024,\n \"spectral_width\": 1e4, # in Hz\n \"reference_offset\": 1e4, # in Hz\n \"events\": [\n {\n \"transition_query\": [\n {\"P\": [-1], \"D\": [2]}, # <-- select inter satellite transitions\n {\"P\": [-1], \"D\": [4]}, # <-- select outer satellite transitions\n ]\n }\n ],\n }\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3))\nmethod2.plot()\nplt.show()" + "method2 = Method1D(\n channels=[\"27Al\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=1e9, # in Hz\n spectral_dimensions=[\n {\n \"count\": 1024,\n \"spectral_width\": 1e4, # in Hz\n \"reference_offset\": 1e4, # in Hz\n \"events\": [\n {\n \"transition_query\": [\n {\"P\": [-1], \"D\": [2]}, # <-- select inter satellite transitions\n {\"P\": [-1], \"D\": [4]}, # <-- select outer satellite transitions\n ]\n }\n ],\n }\n ],\n)" ] }, { @@ -179,7 +179,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.ipynb index 4783bbc69..f84cfd5ea 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.ipynb @@ -58,7 +58,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Selecting the triple-quantum transition\n\nFor single-site spin-5/2 spin system, there are three triple-quantum transition\n\n- $|1/2\\rangle\\rightarrow|-5/2\\rangle$ ($P=-3, D=6$)\n- $|3/2\\rangle\\rightarrow|-3/2\\rangle$ ($P=-3, D=0$)\n- $|5/2\\rangle\\rightarrow|-1/2\\rangle$ ($P=-3, D=-6$)\n\nTo select one or more triple-quantum transitions, assign the respective value of P and\nD to the `transition_query`. Here, we select the symmetric triple-quantum transition.\n\n" + "## Selecting the triple-quantum transition\n\nFor spin-site spin-5/2 spin system, there are three triple-quantum transition\n\n- $|1/2\\rangle\\rightarrow|-5/2\\rangle$ ($P=-3, D=6$)\n- $|3/2\\rangle\\rightarrow|-3/2\\rangle$ ($P=-3, D=0$)\n- $|5/2\\rangle\\rightarrow|-1/2\\rangle$ ($P=-3, D=-6$)\n\nTo select one or more triple-quantum transitions, assign the respective value of P and\nD to the `transition_query`.\n\n" ] }, { @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "method = Method1D(\n name=\"Arbitrary Transition Method\",\n channels=[\"27Al\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=1e9, # in Hz\n spectral_dimensions=[\n {\n \"count\": 1024,\n \"spectral_width\": 5e3, # in Hz\n \"reference_offset\": 2.5e4, # in Hz\n \"events\": [\n { # symmetric triple quantum transitions\n \"transition_query\": [{\"P\": [-3], \"D\": [0]}]\n }\n ],\n }\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3))\nmethod.plot()\nplt.show()" + "method = Method1D(\n channels=[\"27Al\"],\n magnetic_flux_density=21.14, # in T\n rotor_frequency=1e9, # in Hz\n spectral_dimensions=[\n {\n \"count\": 1024,\n \"spectral_width\": 5e3, # in Hz\n \"reference_offset\": 2.5e4, # in Hz\n \"events\": [\n { # symmetric triple quantum transitions\n \"transition_query\": [{\"P\": [-3], \"D\": [0]}]\n }\n ],\n }\n ],\n)" ] }, { @@ -107,7 +107,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_5_dipolar_coupled.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_5_dipolar_coupled.ipynb index 3e75b1f91..e9ebed514 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_5_dipolar_coupled.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_5_dipolar_coupled.ipynb @@ -136,7 +136,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_6_coupled_spin_system.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_6_coupled_spin_system.ipynb index 9be0a6d1b..b4f017969 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_6_coupled_spin_system.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_6_coupled_spin_system.ipynb @@ -136,7 +136,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_7_CAS_dipolar_J.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_7_CAS_dipolar_J.ipynb index 192724a33..cfecd3c67 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_7_CAS_dipolar_J.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_7_CAS_dipolar_J.ipynb @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(crystalline)/plot_8_Ethanol.ipynb b/docs/notebooks/examples/1D_simulation(crystalline)/plot_8_Ethanol.ipynb index e4b4737c3..cca6c031c 100644 --- a/docs/notebooks/examples/1D_simulation(crystalline)/plot_8_Ethanol.ipynb +++ b/docs/notebooks/examples/1D_simulation(crystalline)/plot_8_Ethanol.ipynb @@ -222,7 +222,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_0_protein_GB1.ipynb b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_0_protein_GB1.ipynb index 8902a1f7b..972629811 100644 --- a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_0_protein_GB1.ipynb +++ b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_0_protein_GB1.ipynb @@ -179,7 +179,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_3_amorphous_like.ipynb b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_3_amorphous_like.ipynb index d994b5305..a71b90a11 100644 --- a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_3_amorphous_like.ipynb +++ b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_3_amorphous_like.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "spin_systems = single_site_system_generator(\n isotope=\"29Si\",\n isotropic_chemical_shift=iso,\n shielding_symmetric={\"zeta\": zeta, \"eta\": eta},\n abundance=pdf,\n)" + "spin_systems = single_site_system_generator(\n isotopes=\"29Si\",\n isotropic_chemical_shifts=iso,\n shielding_symmetric={\"zeta\": zeta, \"eta\": eta},\n abundance=pdf,\n)" ] }, { @@ -258,7 +258,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.ipynb b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.ipynb index f021f209e..8ed89c8c5 100644 --- a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.ipynb +++ b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.ipynb @@ -80,7 +80,7 @@ }, "outputs": [], "source": [ - "spin_systems = single_site_system_generator(\n isotope=\"27Al\",\n isotropic_chemical_shift=iso,\n quadrupolar={\"Cq\": Cq * 1e6, \"eta\": eta}, # Cq in Hz\n abundance=pdf,\n)" + "spin_systems = single_site_system_generator(\n isotopes=\"27Al\",\n isotropic_chemical_shifts=iso,\n quadrupolar={\"Cq\": Cq * 1e6, \"eta\": eta}, # Cq in Hz\n abundance=pdf,\n)" ] }, { @@ -208,7 +208,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.ipynb b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.ipynb index 825b937d2..53159aa9a 100644 --- a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.ipynb +++ b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "systems = single_site_system_generator(\n isotope=\"13C\", shielding_symmetric={\"zeta\": z_dist, \"eta\": e_dist}, abundance=amp\n)" + "systems = single_site_system_generator(\n isotopes=\"13C\", shielding_symmetric={\"zeta\": z_dist, \"eta\": e_dist}, abundance=amp\n)" ] }, { @@ -159,7 +159,7 @@ }, "outputs": [], "source": [ - "systems = single_site_system_generator(\n isotope=\"71Ga\", quadrupolar={\"Cq\": cq_dist * 1e6, \"eta\": e_dist}, abundance=amp\n)" + "systems = single_site_system_generator(\n isotopes=\"71Ga\", quadrupolar={\"Cq\": cq_dist * 1e6, \"eta\": e_dist}, abundance=amp\n)" ] }, { @@ -215,7 +215,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.ipynb b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.ipynb index 7b06fc483..e41c9a9cd 100644 --- a/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.ipynb +++ b/docs/notebooks/examples/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "systems = single_site_system_generator(\n isotope=\"13C\", shielding_symmetric={\"zeta\": z_dist, \"eta\": e_dist}, abundance=amp\n)\nprint(len(systems))" + "systems = single_site_system_generator(\n isotopes=\"13C\", shielding_symmetric={\"zeta\": z_dist, \"eta\": e_dist}, abundance=amp\n)\nprint(len(systems))" ] }, { @@ -177,7 +177,7 @@ }, "outputs": [], "source": [ - "systems = single_site_system_generator(\n isotope=\"71Ga\", quadrupolar={\"Cq\": cq_dist * 1e6, \"eta\": e_dist}, abundance=amp\n)" + "systems = single_site_system_generator(\n isotopes=\"71Ga\", quadrupolar={\"Cq\": cq_dist * 1e6, \"eta\": e_dist}, abundance=amp\n)" ] }, { @@ -269,7 +269,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.ipynb index f83d6ffed..75d601341 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "method = ThreeQ_VAS(\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4, # in T\n spectral_dimensions=[\n {\n \"count\": 128,\n \"spectral_width\": 7e3, # in Hz\n \"reference_offset\": -7e3, # in Hz\n \"label\": \"Isotropic dimension\",\n },\n {\n \"count\": 256,\n \"spectral_width\": 1e4, # in Hz\n \"reference_offset\": -4e3, # in Hz\n \"label\": \"MAS dimension\",\n },\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\nmethod.plot()\nplt.show()" + "method = ThreeQ_VAS(\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4, # in T\n spectral_dimensions=[\n {\n \"count\": 128,\n \"spectral_width\": 7e3, # in Hz\n \"reference_offset\": -7e3, # in Hz\n \"label\": \"Isotropic dimension\",\n },\n {\n \"count\": 256,\n \"spectral_width\": 1e4, # in Hz\n \"reference_offset\": -4e3, # in Hz\n \"label\": \"MAS dimension\",\n },\n ],\n)" ] }, { @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_MQMAS_albite.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_MQMAS_albite.ipynb index e3dcc7342..75138710b 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_MQMAS_albite.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_MQMAS_albite.ipynb @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.ipynb index 58405e3d3..3aa4b8311 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "angles = [54.7359, 54.73]\nmethod = []\nfor angle in angles:\n method.append(\n ST1_VAS(\n channels=[\"87Rb\"],\n magnetic_flux_density=7, # in T\n rotor_angle=angle * 3.14159 / 180, # in rad (magic angle)\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 3e3, # in Hz\n \"reference_offset\": -2.4e3, # in Hz\n \"label\": \"Isotropic dimension\",\n },\n {\n \"count\": 512,\n \"spectral_width\": 5e3, # in Hz\n \"reference_offset\": -4e3, # in Hz\n \"label\": \"MAS dimension\",\n },\n ],\n )\n )\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\nmethod[0].plot()\nplt.show()" + "angles = [54.7359, 54.73]\nmethod = []\nfor angle in angles:\n method.append(\n ST1_VAS(\n channels=[\"87Rb\"],\n magnetic_flux_density=7, # in T\n rotor_angle=angle * 3.14159 / 180, # in rad (magic angle)\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 3e3, # in Hz\n \"reference_offset\": -2.4e3, # in Hz\n \"label\": \"Isotropic dimension\",\n },\n {\n \"count\": 512,\n \"spectral_width\": 5e3, # in Hz\n \"reference_offset\": -4e3, # in Hz\n \"label\": \"MAS dimension\",\n },\n ],\n )\n )" ] }, { @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.ipynb index 359073587..dfed8c7cd 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "sas = Method2D(\n name=\"Switched Angle Spinning\",\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4, # in T\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 3.5e4, # in Hz\n \"reference_offset\": 1e3, # in Hz\n \"label\": \"90 dimension\",\n \"events\": [\n {\n \"rotor_angle\": 90 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ], # in radians\n },\n {\n \"count\": 256,\n \"spectral_width\": 22e3, # in Hz\n \"reference_offset\": -4e3, # in Hz\n \"label\": \"MAS dimension\",\n \"events\": [\n {\n \"rotor_angle\": 54.74 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ], # in radians\n },\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\nsas.plot()\nplt.show()" + "sas = Method2D(\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4, # in T\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 3.5e4, # in Hz\n \"reference_offset\": 1e3, # in Hz\n \"label\": \"90 dimension\",\n \"events\": [\n {\n \"rotor_angle\": 90 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ], # in radians\n },\n {\n \"count\": 256,\n \"spectral_width\": 22e3, # in Hz\n \"reference_offset\": -4e3, # in Hz\n \"label\": \"MAS dimension\",\n \"events\": [\n {\n \"rotor_angle\": 54.74 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ], # in radians\n },\n ],\n)" ] }, { @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_3_SAS_Rb2CrO4.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_3_SAS_Rb2CrO4.ipynb index c1fa4748c..ae5ab414f 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_3_SAS_Rb2CrO4.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_3_SAS_Rb2CrO4.ipynb @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_3QMAS_Coesite.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_3QMAS_Coesite.ipynb index 2b4bd10a5..b83707e1c 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_3QMAS_Coesite.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_3QMAS_Coesite.ipynb @@ -132,7 +132,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_DAS_Coesite.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_DAS_Coesite.ipynb index 42db7af59..32be6f646 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_DAS_Coesite.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_4_DAS_Coesite.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "das = Method2D(\n name=\"Dynamic Angle Spinning\",\n channels=[\"17O\"],\n magnetic_flux_density=11.74, # in T\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 5e3, # in Hz\n \"reference_offset\": 0, # in Hz\n \"label\": \"DAS isotropic dimension\",\n \"events\": [\n {\n \"fraction\": 0.5,\n \"rotor_angle\": 37.38 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n },\n {\n \"fraction\": 0.5,\n \"rotor_angle\": 79.19 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n },\n ],\n },\n # The last spectral dimension block is the direct-dimension\n {\n \"count\": 256,\n \"spectral_width\": 2e4, # in Hz\n \"reference_offset\": 0, # in Hz\n \"label\": \"MAS dimension\",\n \"events\": [\n {\n \"rotor_angle\": 54.735 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ],\n },\n ],\n)\nsim.methods = [das] # add the method\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\ndas.plot()\nplt.show()" + "das = Method2D(\n channels=[\"17O\"],\n magnetic_flux_density=11.74, # in T\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 5e3, # in Hz\n \"reference_offset\": 0, # in Hz\n \"label\": \"DAS isotropic dimension\",\n \"events\": [\n {\n \"fraction\": 0.5,\n \"rotor_angle\": 37.38 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n },\n {\n \"fraction\": 0.5,\n \"rotor_angle\": 79.19 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n },\n ],\n },\n # The last spectral dimension block is the direct-dimension\n {\n \"count\": 256,\n \"spectral_width\": 2e4, # in Hz\n \"reference_offset\": 0, # in Hz\n \"label\": \"MAS dimension\",\n \"events\": [\n {\n \"rotor_angle\": 54.735 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ],\n },\n ],\n)\nsim.methods = [das] # add the method." ] }, { @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.ipynb index db6d44899..3289381e4 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "coaster = Method2D(\n name=\"COASTER\",\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4, # in T\n rotor_angle=70.12 * 3.14159 / 180, # in rads\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 4e4, # in Hz\n \"reference_offset\": -8e3, # in Hz\n \"label\": \"3Q dimension\",\n \"events\": [{\"transition_query\": [{\"P\": [3], \"D\": [0]}]}],\n },\n # The last spectral dimension block is the direct-dimension\n {\n \"count\": 256,\n \"spectral_width\": 2e4, # in Hz\n \"reference_offset\": -3e3, # in Hz\n \"label\": \"70.12 dimension\",\n \"events\": [{\"transition_query\": [{\"P\": [-1], \"D\": [0]}]}],\n },\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\ncoaster.plot()\nplt.show()" + "coaster = Method2D(\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4, # in T\n rotor_angle=70.12 * 3.14159 / 180, # in rads\n spectral_dimensions=[\n {\n \"count\": 256,\n \"spectral_width\": 4e4, # in Hz\n \"reference_offset\": -8e3, # in Hz\n \"label\": \"3Q dimension\",\n \"events\": [{\"transition_query\": [{\"P\": [3], \"D\": [0]}]}],\n },\n # The last spectral dimension block is the direct-dimension\n {\n \"count\": 256,\n \"spectral_width\": 2e4, # in Hz\n \"reference_offset\": -3e3, # in Hz\n \"label\": \"70.12 dimension\",\n \"events\": [{\"transition_query\": [{\"P\": [-1], \"D\": [0]}]}],\n },\n ],\n)" ] }, { @@ -168,7 +168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.ipynb index c00f65cd7..7a13ced8b 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.ipynb @@ -69,25 +69,7 @@ }, "outputs": [], "source": [ - "PASS = SSB2D(\n channels=[\"13C\"],\n magnetic_flux_density=11.74,\n rotor_frequency=2000,\n spectral_dimensions=[\n {\n \"count\": 20 * 4,\n \"spectral_width\": 2000 * 20, # value in Hz\n \"label\": \"Anisotropic dimension\",\n },\n {\n \"count\": 1024,\n \"spectral_width\": 3e4, # value in Hz\n \"reference_offset\": 1.1e4, # value in Hz\n \"label\": \"Isotropic dimension\",\n },\n ],\n)\nsim.methods = [PASS] # add the method.\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\nPASS.plot()\nplt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For 2D spinning sideband simulation, set the number of spinning sidebands in the\nSimulator.config object to `spectral_width/rotor_frequency` along the sideband\ndimension.\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "sim.config.number_of_sidebands = 20\n\n# run the simulation.\nsim.run()" + "PASS = SSB2D(\n channels=[\"13C\"],\n magnetic_flux_density=11.74,\n rotor_frequency=2000,\n spectral_dimensions=[\n {\n \"count\": 20 * 4,\n \"spectral_width\": 2000 * 20, # value in Hz\n \"label\": \"Anisotropic dimension\",\n },\n {\n \"count\": 1024,\n \"spectral_width\": 3e4, # value in Hz\n \"reference_offset\": 1.1e4, # value in Hz\n \"label\": \"Isotropic dimension\",\n },\n ],\n)\nsim.methods = [PASS] # add the method.\n\n# For 2D spinning sideband simulation, set the number of spinning sidebands in the\n# Simulator.config object to `spectral_width/rotor_frequency` along the sideband\n# dimension.\nsim.config.number_of_sidebands = 20\n\n# run the simulation.\nsim.run()" ] }, { @@ -150,7 +132,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.ipynb index 8f393698c..890ed9498 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "qmat = SSB2D(\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4,\n rotor_frequency=2604,\n spectral_dimensions=[\n {\n \"count\": 32 * 4,\n \"spectral_width\": 2604 * 32, # in Hz\n \"label\": \"Anisotropic dimension\",\n },\n {\n \"count\": 512,\n \"spectral_width\": 50000, # in Hz\n \"label\": \"Fast MAS dimension\",\n },\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\nqmat.plot()\nplt.show()" + "qmat = SSB2D(\n channels=[\"87Rb\"],\n magnetic_flux_density=9.4,\n rotor_frequency=2604,\n spectral_dimensions=[\n {\n \"count\": 32 * 4,\n \"spectral_width\": 2604 * 32, # in Hz\n \"label\": \"Anisotropic dimension\",\n },\n {\n \"count\": 512,\n \"spectral_width\": 50000, # in Hz\n \"label\": \"Fast MAS dimension\",\n },\n ],\n)" ] }, { @@ -132,7 +132,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_8_MAF.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_8_MAF.ipynb index 52d0cd3fe..e888ebd3a 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_8_MAF.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_8_MAF.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "maf = Method2D(\n name=\"Magic Angle Spinning\",\n channels=[\"29Si\"],\n magnetic_flux_density=14.1, # in T\n spectral_dimensions=[\n {\n \"count\": 128,\n \"spectral_width\": 2e4, # in Hz\n \"label\": \"Anisotropic dimension\",\n \"events\": [\n {\n \"rotor_angle\": 90 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ],\n },\n {\n \"count\": 128,\n \"spectral_width\": 3e3, # in Hz\n \"reference_offset\": -1.05e4, # in Hz\n \"label\": \"Isotropic dimension\",\n \"events\": [\n {\n \"rotor_angle\": 54.735 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ],\n },\n ],\n affine_matrix=[[1, -1], [0, 1]],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 3.5))\nmaf.plot()\nplt.show()" + "maf = Method2D(\n channels=[\"29Si\"],\n magnetic_flux_density=14.1, # in T\n spectral_dimensions=[\n {\n \"count\": 128,\n \"spectral_width\": 2e4, # in Hz\n \"label\": \"Anisotropic dimension\",\n \"events\": [\n {\n \"rotor_angle\": 90 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ],\n },\n {\n \"count\": 128,\n \"spectral_width\": 3e3, # in Hz\n \"reference_offset\": -1.05e4, # in Hz\n \"label\": \"Isotropic dimension\",\n \"events\": [\n {\n \"rotor_angle\": 54.735 * 3.14159 / 180,\n \"transition_query\": [{\"P\": [-1], \"D\": [0]}],\n }\n ],\n },\n ],\n affine_matrix=[[1, -1], [0, 1]],\n)" ] }, { @@ -150,7 +150,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(crystalline)/plot_9_shifting-d.ipynb b/docs/notebooks/examples/2D_simulation(crystalline)/plot_9_shifting-d.ipynb index 18e22d30e..9f04834d6 100644 --- a/docs/notebooks/examples/2D_simulation(crystalline)/plot_9_shifting-d.ipynb +++ b/docs/notebooks/examples/2D_simulation(crystalline)/plot_9_shifting-d.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "shifting_d = Method2D(\n name=\"Shifting-d\",\n channels=[\"2H\"],\n magnetic_flux_density=9.395, # in T\n spectral_dimensions=[\n {\n \"count\": 512,\n \"spectral_width\": 2.5e5, # in Hz\n \"label\": \"Quadrupolar frequency\",\n \"events\": [\n {\n \"rotor_frequency\": 0,\n \"transition_query\": {\"P\": [-1]},\n \"freq_contrib\": [\"Quad1_2\"],\n }\n ],\n },\n {\n \"count\": 256,\n \"spectral_width\": 2e5, # in Hz\n \"reference_offset\": 2e4, # in Hz\n \"label\": \"Paramagnetic shift\",\n \"events\": [\n {\n \"rotor_frequency\": 0,\n \"transition_query\": {\"P\": [-1]},\n \"freq_contrib\": [\"Shielding1_0\", \"Shielding1_2\"],\n }\n ],\n },\n ],\n)\n\n# A graphical representation of the method object.\nplt.figure(figsize=(5, 2.5))\nshifting_d.plot()\nplt.show()" + "shifting_d = Method2D(\n channels=[\"2H\"],\n magnetic_flux_density=9.395, # in T\n spectral_dimensions=[\n {\n \"count\": 512,\n \"spectral_width\": 2.5e5, # in Hz\n \"label\": \"Quadrupolar frequency\",\n \"events\": [\n {\n \"rotor_frequency\": 0,\n \"transition_query\": {\"P\": [-1]},\n \"freq_contrib\": [\"Quad1_2\"],\n }\n ],\n },\n {\n \"count\": 256,\n \"spectral_width\": 2e5, # in Hz\n \"reference_offset\": 2e4, # in Hz\n \"label\": \"Paramagnetic shift\",\n \"events\": [\n {\n \"rotor_frequency\": 0,\n \"transition_query\": {\"P\": [-1]},\n \"freq_contrib\": [\"Shielding1_0\", \"Shielding1_2\"],\n }\n ],\n },\n ],\n)" ] }, { @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.ipynb b/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.ipynb index f2df2ebc0..3e92603c9 100644 --- a/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.ipynb +++ b/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "spin_systems = single_site_system_generator(\n isotope=\"87Rb\",\n isotropic_chemical_shift=iso,\n quadrupolar={\"Cq\": Cq * 1e6, \"eta\": eta}, # Cq in Hz\n abundance=pdf,\n)\nlen(spin_systems)" + "spin_systems = single_site_system_generator(\n isotopes=\"87Rb\",\n isotropic_chemical_shifts=iso,\n quadrupolar={\"Cq\": Cq * 1e6, \"eta\": eta}, # Cq in Hz\n abundance=pdf,\n)\nlen(spin_systems)" ] }, { @@ -161,7 +161,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_1_I=2.5.ipynb b/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_1_I=2.5.ipynb index 1ff7046f6..fff96ba6c 100644 --- a/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_1_I=2.5.ipynb +++ b/docs/notebooks/examples/2D_simulation(macro_amorphous)/plot_1_I=2.5.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "spin_systems = single_site_system_generator(\n isotope=\"27Al\",\n isotropic_chemical_shift=iso,\n quadrupolar={\"Cq\": Cq * 1e6, \"eta\": eta}, # Cq in Hz\n abundance=pdf,\n)\nlen(spin_systems)" + "spin_systems = single_site_system_generator(\n isotopes=\"27Al\",\n isotropic_chemical_shifts=iso,\n quadrupolar={\"Cq\": Cq * 1e6, \"eta\": eta}, # Cq in Hz\n abundance=pdf,\n)\nlen(spin_systems)" ] }, { @@ -161,7 +161,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_1_29Si_cuspidine.ipynb b/docs/notebooks/fitting/1D_fitting/plot_1_29Si_cuspidine.ipynb index d04ce6e09..45478dfbb 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_1_29Si_cuspidine.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_1_29Si_cuspidine.ipynb @@ -240,7 +240,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_1_31P_Na2PO4_MAS.ipynb b/docs/notebooks/fitting/1D_fitting/plot_1_31P_Na2PO4_MAS.ipynb index fa54ff234..5f5fa9af2 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_1_31P_Na2PO4_MAS.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_1_31P_Na2PO4_MAS.ipynb @@ -240,7 +240,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_1_31p_Na2PO4_static.ipynb b/docs/notebooks/fitting/1D_fitting/plot_1_31p_Na2PO4_static.ipynb index b6a41078f..d066e3409 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_1_31p_Na2PO4_static.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_1_31p_Na2PO4_static.ipynb @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_1940Hz.ipynb b/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_1940Hz.ipynb index 6085865ee..63a5c8db5 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_1940Hz.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_1940Hz.ipynb @@ -51,7 +51,7 @@ }, "outputs": [], "source": [ - "host = \"https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/\"\nfilename = \"13C MAS 1940Hz - Glycine.csdf\"\nexperiment = cp.load(host + filename)\n\n# standard deviation of noise from the dataset\nsigma = 3.890551\n\n# For spectral fitting, we only focus on the real part of the complex dataset\nexperiment = experiment.real\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n# plot of the dataset.\nplt.figure(figsize=(8, 4))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, color=\"black\", linewidth=0.5, label=\"Experiment\")\nax.set_xlim(280, -10)\nplt.grid()\nplt.tight_layout()\nplt.show()" + "host = \"https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/\"\nfilename = \"13C MAS 1940Hz - Glycine.csdf\"\nexperiment = cp.load(host + filename)\n\n# standard deviation of noise from the dataset\nsigma = 3.822249\n\n# For spectral fitting, we only focus on the real part of the complex dataset\nexperiment = experiment.real\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n# plot of the dataset.\nplt.figure(figsize=(8, 4))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, color=\"black\", linewidth=0.5, label=\"Experiment\")\nax.set_xlim(280, -10)\nplt.grid()\nplt.tight_layout()\nplt.show()" ] }, { @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_5000Hz.ipynb b/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_5000Hz.ipynb index bebff05c0..ed9c23616 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_5000Hz.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_5000Hz.ipynb @@ -51,7 +51,7 @@ }, "outputs": [], "source": [ - "host = \"https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/\"\nfilename = \"13C MAS 5000Hz - Glycine.csdf\"\nexperiment = cp.load(host + filename)\n\n# standard deviation of noise from the dataset\nsigma = 1.97637\n\n# For spectral fitting, we only focus on the real part of the complex dataset\nexperiment = experiment.real\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n# plot of the dataset.\nplt.figure(figsize=(8, 4))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, color=\"black\", linewidth=0.5, label=\"Experiment\")\nax.set_xlim(280, -10)\nplt.grid()\nplt.tight_layout()\nplt.show()" + "host = \"https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/\"\nfilename = \"13C MAS 5000Hz - Glycine.csdf\"\nexperiment = cp.load(host + filename)\n\n# standard deviation of noise from the dataset\nsigma = 3.822249\n\n# For spectral fitting, we only focus on the real part of the complex dataset\nexperiment = experiment.real\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n# plot of the dataset.\nplt.figure(figsize=(8, 4))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, color=\"black\", linewidth=0.5, label=\"Experiment\")\nax.set_xlim(280, -10)\nplt.grid()\nplt.tight_layout()\nplt.show()" ] }, { @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_960Hz.ipynb b/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_960Hz.ipynb index b408db144..fd3b35bcd 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_960Hz.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_960Hz.ipynb @@ -51,7 +51,7 @@ }, "outputs": [], "source": [ - "host = \"https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/\"\nfilename = \"13C MAS 960Hz - Glycine.csdf\"\nexperiment = cp.load(host + filename)\n\n# standard deviation of noise from the dataset\nsigma = 3.982936\n\n# For spectral fitting, we only focus on the real part of the complex dataset\nexperiment = experiment.real\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n# plot of the dataset.\nplt.figure(figsize=(8, 4))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, color=\"black\", linewidth=0.5, label=\"Experiment\")\nax.set_xlim(280, -10)\nplt.grid()\nplt.tight_layout()\nplt.show()" + "host = \"https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/\"\nfilename = \"13C MAS 960Hz - Glycine.csdf\"\nexperiment = cp.load(host + filename)\n\n# standard deviation of noise from the dataset\nsigma = 3.822249\n\n# For spectral fitting, we only focus on the real part of the complex dataset\nexperiment = experiment.real\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n# plot of the dataset.\nplt.figure(figsize=(8, 4))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, color=\"black\", linewidth=0.5, label=\"Experiment\")\nax.set_xlim(280, -10)\nplt.grid()\nplt.tight_layout()\nplt.show()" ] }, { @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_multi_spectra_fit.ipynb b/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_multi_spectra_fit.ipynb deleted file mode 100644 index 2ad7ddee5..000000000 --- a/docs/notebooks/fitting/1D_fitting/plot_2_13C_glycine_multi_spectra_fit.ipynb +++ /dev/null @@ -1,194 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "# This cell is added by sphinx-gallery\n!pip install mrsimulator --quiet\n\n\n%matplotlib inline\n\nimport mrsimulator\nprint(f'You are using mrsimulator v{mrsimulator.__version__}')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n# \u00b9\u00b3C MAS NMR of Glycine (CSA) multi-spectra fit\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The following is a multi-data least-squares fitting example of\n$^{13}\\text{C}$ MAS NMR spectrum of Glycine spinning at 5 kHz, 1.94 kHz, and\n960 Hz. Before trying multi-data fitting, we recommend that you first try individual\nfits. The experimental datasets are part of DMFIT [#f1]_ examples.\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "import csdmpy as cp\nimport matplotlib.pyplot as plt\nfrom lmfit import Minimizer\n\nfrom mrsimulator import Simulator, SpinSystem, Site\nfrom mrsimulator.methods import BlochDecaySpectrum\nfrom mrsimulator import signal_processing as sp\nfrom mrsimulator.utils import spectral_fitting as sf\nfrom mrsimulator.utils import get_spectral_dimensions" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Import the datasets\nImport the datasets and assign the standard deviation of noise for each dataset. Here,\n``sigma1``, ``sigma2``, and ``sigma3`` are the noise standard deviation for the\ndataset acquired at 5 kHz, 1.94 kHz, and 960 Hz spinning speeds, respectively.\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "host = \"https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/\"\nfilename1 = \"13C MAS 5000Hz - Glycine.csdf\"\nfilename2 = \"13C MAS 1940Hz - Glycine.csdf\"\nfilename3 = \"13C MAS 960Hz - Glycine.csdf\"\n\nexperiment1 = cp.load(host + filename1).real\nexperiment2 = cp.load(host + filename2).real\nexperiment3 = cp.load(host + filename3).real\nexperiments = [experiment1, experiment2, experiment3]\n\n# standard deviation of noise from the dataset\nsigma1 = 1.97637\nsigma2 = 3.822249\nsigma3 = 3.982936\nsigmas = [sigma1, sigma2, sigma3]\n\n# Convert the coordinates along each dimension from Hz to ppm for each dataset\nfig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={\"projection\": \"csdm\"})\nfor i, experiment in enumerate(experiments):\n _ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n # plot of the dataset.\n ax[i].plot(experiment, color=\"black\", linewidth=0.5, label=\"Experiment\")\n ax[i].set_xlim(280, -10)\n ax[i].grid()\nplt.tight_layout()\nplt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Create a fitting model\n**Spin System**: The objective of a multi-data fitting is to optimize the spin\nsystem parameters using multiple datasets. In this example, we create two single-site\nspin systems, which are then shared by three method objects.\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "C1 = Site(\n isotope=\"13C\",\n isotropic_chemical_shift=176.0, # in ppm\n shielding_symmetric={\"zeta\": 60, \"eta\": 0.6}, # zeta in Hz\n)\nC2 = Site(\n isotope=\"13C\",\n isotropic_chemical_shift=43.0, # in ppm\n shielding_symmetric={\"zeta\": 30, \"eta\": 0.5}, # zeta in Hz\n)\n\nspin_systems = [SpinSystem(sites=[C1], name=\"C1\"), SpinSystem(sites=[C2], name=\"C2\")]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Method**: Create the three MAS method objects with respective MAS spinning speeds.\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "# Get the spectral dimension parameters from the respective experiment and setup the\n# corresponding method.\n\n# Method for dataset 1\nspectral_dims1 = get_spectral_dimensions(experiment1)\nMAS1 = BlochDecaySpectrum(\n channels=[\"13C\"],\n magnetic_flux_density=7.05, # in T\n rotor_frequency=5000, # in Hz\n spectral_dimensions=spectral_dims1,\n experiment=experiment1, # add experimental dataset 1\n)\n\n# Method for dataset 2\nspectral_dims2 = get_spectral_dimensions(experiment2)\nMAS2 = BlochDecaySpectrum(\n channels=[\"13C\"],\n magnetic_flux_density=7.05, # in T\n rotor_frequency=1940, # in Hz\n spectral_dimensions=spectral_dims2,\n experiment=experiment2, # add experimental dataset 2\n)\n\n# Method for dataset 3\nspectral_dims3 = get_spectral_dimensions(experiment3)\nMAS3 = BlochDecaySpectrum(\n channels=[\"13C\"],\n magnetic_flux_density=7.05, # in T\n rotor_frequency=960, # in Hz\n spectral_dimensions=spectral_dims3,\n experiment=experiment3, # add experimental dataset 3\n)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Guess Model Spectrum**\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "# Simulation\n# ----------\n# Add the spin systems and the three methods to the simulator object.\nsim = Simulator(spin_systems=spin_systems, methods=[MAS1, MAS2, MAS3])\nsim.config.decompose_spectrum = \"spin_system\"\nsim.run()\n\n# Post Simulation Processing\n# --------------------------\n# Add signal processing to simulation dataset from the three methods.\n\n# Processor for dataset 1\nprocessor1 = sp.SignalProcessor(\n operations=[\n sp.IFFT(),\n sp.apodization.Exponential(FWHM=\"20 Hz\", dv_index=0), # spin system 0\n sp.apodization.Exponential(FWHM=\"200 Hz\", dv_index=1), # spin system 1\n sp.FFT(),\n sp.Scale(factor=10), # dataset is scaled independently using scale factor.\n ]\n)\n\n# Processor for dataset 2\nprocessor2 = sp.SignalProcessor(\n operations=[\n sp.IFFT(),\n sp.apodization.Exponential(FWHM=\"30 Hz\", dv_index=0), # spin system 0\n sp.apodization.Exponential(FWHM=\"300 Hz\", dv_index=1), # spin system 1\n sp.FFT(),\n sp.Scale(factor=100), # dataset is scaled independently using scale factor.\n ]\n)\n\n# Processor for dataset 3\nprocessor3 = sp.SignalProcessor(\n operations=[\n sp.IFFT(),\n sp.apodization.Exponential(FWHM=\"10 Hz\", dv_index=0), # spin system 0\n sp.apodization.Exponential(FWHM=\"150 Hz\", dv_index=1), # spin system 1\n sp.FFT(),\n sp.Scale(factor=50), # dataset is scaled independently using scale factor.\n ]\n)\nprocessors = [processor1, processor2, processor3]\n\nprocessed_data = []\nfor i, proc in enumerate(processors):\n processed_data.append(proc.apply_operations(data=sim.methods[i].simulation).real)\n\n\n# Plot of the guess Spectrum\n# --------------------------\n\nfig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={\"projection\": \"csdm\"})\nfor i, exp_data in enumerate(experiments):\n ax[i].plot(exp_data, color=\"black\", linewidth=0.5, label=\"Experiment\")\n ax[i].plot(processed_data[i], linewidth=2, alpha=0.6)\n ax[i].set_xlim(280, -10)\n ax[i].grid()\n\nplt.legend()\nplt.tight_layout()\nplt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Least-squares minimization with LMFIT\nUse the :func:`~mrsimulator.utils.spectral_fitting.make_LMFIT_params` for a quick\nsetup of the fitting parameters. Note, the first two arguments of this function is\nthe simulator object and a list of signal processing objects, ``processors``. The\nfitting parameters corresponding to the signal processor objects are generated using\n``SP_i_operation_j_FunctionName_FunctionArg``, where *i* is the *ith* signal\nprocessor within the list, *j* is the operation index of the *ith* processor, and\n*FunctionName* and *FunctionArg* are the operation function name and function\nargument, respectively.\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "params = sf.make_LMFIT_params(sim, processors, include={\"rotor_frequency\"})\nprint(params.pretty_print(columns=[\"value\", \"min\", \"max\", \"vary\", \"expr\"]))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Solve the minimizer using LMFIT**\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processors, sigmas))\nresult = minner.minimize()\nresult" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## The best fit solution\n\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "all_best_fit = sf.bestfit(sim, processors) # a list of best fit simulations\nall_residuals = sf.residuals(sim, processors) # a list of residuals\n\n# Plot the spectrum\nfig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={\"projection\": \"csdm\"})\nfor i, proc in enumerate(processors):\n ax[i].plot(experiments[i], color=\"black\", linewidth=0.5, label=\"Experiment\")\n ax[i].plot(all_residuals[i], color=\"gray\", linewidth=0.5, label=\"Residual\")\n ax[i].plot(all_best_fit[i], linewidth=2, alpha=0.6)\n ax[i].set_xlim(280, -10)\n ax[i].grid()\n\nplt.legend()\nplt.tight_layout()\nplt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - ".. [#f1] D.Massiot, F.Fayon, M.Capron, I.King, S.Le Calv\u00e9, B.Alonso, J.O.Durand,\n B.Bujoli, Z.Gan, G.Hoatson, 'Modelling one and two-dimensional solid-state NMR\n spectra.', Magn. Reson. Chem. **40** 70-76 (2002)\n `DOI: 10.1002/mrc.984 `_\n\n" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.11" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} \ No newline at end of file diff --git a/docs/notebooks/fitting/1D_fitting/plot_2_PASS_cross_sections.ipynb b/docs/notebooks/fitting/1D_fitting/plot_2_PASS_cross_sections.ipynb index e8f67dfdb..e1615ef34 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_2_PASS_cross_sections.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_2_PASS_cross_sections.ipynb @@ -179,7 +179,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_3_Na2SiO3.ipynb b/docs/notebooks/fitting/1D_fitting/plot_3_Na2SiO3.ipynb index eaaf11d75..32efa1a09 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_3_Na2SiO3.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_3_Na2SiO3.ipynb @@ -240,7 +240,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_4_11B_Quad_NMR.ipynb b/docs/notebooks/fitting/1D_fitting/plot_4_11B_Quad_NMR.ipynb index 696f24edf..ef1274181 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_4_11B_Quad_NMR.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_4_11B_Quad_NMR.ipynb @@ -179,7 +179,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_4_23Na_Nasicon.ipynb b/docs/notebooks/fitting/1D_fitting/plot_4_23Na_Nasicon.ipynb index f21ab4695..bcbc83d89 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_4_23Na_Nasicon.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_4_23Na_Nasicon.ipynb @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_4_27Al_YAG.ipynb b/docs/notebooks/fitting/1D_fitting/plot_4_27Al_YAG.ipynb index a05769a00..ab80aca0d 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_4_27Al_YAG.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_4_27Al_YAG.ipynb @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_4_2H_quad.ipynb b/docs/notebooks/fitting/1D_fitting/plot_4_2H_quad.ipynb index b67e74fba..db0d1a263 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_4_2H_quad.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_4_2H_quad.ipynb @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/1D_fitting/plot_5_119Sn_sideband.ipynb b/docs/notebooks/fitting/1D_fitting/plot_5_119Sn_sideband.ipynb index 20df5a42d..ee53b2d14 100644 --- a/docs/notebooks/fitting/1D_fitting/plot_5_119Sn_sideband.ipynb +++ b/docs/notebooks/fitting/1D_fitting/plot_5_119Sn_sideband.ipynb @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/2D_fitting/plot_1_LHistidine_PASS.ipynb b/docs/notebooks/fitting/2D_fitting/plot_1_LHistidine_PASS.ipynb index 0df8f9fc2..a0b731b5b 100644 --- a/docs/notebooks/fitting/2D_fitting/plot_1_LHistidine_PASS.ipynb +++ b/docs/notebooks/fitting/2D_fitting/plot_1_LHistidine_PASS.ipynb @@ -87,7 +87,7 @@ }, "outputs": [], "source": [ - "shifts = [120, 128, 135, 175, 55, 25] # in ppm\nzeta = [-70, -65, -60, -60, -10, -10] # in ppm\neta = [0.8, 0.4, 0.9, 0.3, 0.0, 0.0]\n\nspin_systems = single_site_system_generator(\n isotope=\"13C\",\n isotropic_chemical_shift=shifts,\n shielding_symmetric={\"zeta\": zeta, \"eta\": eta},\n abundance=100 / 6,\n)" + "shifts = [120, 128, 135, 175, 55, 25] # in ppm\nzeta = [-70, -65, -60, -60, -10, -10] # in ppm\neta = [0.8, 0.4, 0.9, 0.3, 0.0, 0.0]\n\nspin_systems = single_site_system_generator(\n isotopes=\"13C\",\n isotropic_chemical_shifts=shifts,\n shielding_symmetric={\"zeta\": zeta, \"eta\": eta},\n abundance=100 / 6,\n)" ] }, { @@ -204,7 +204,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/2D_fitting/plot_1_Rb2SO4_QMAT.ipynb b/docs/notebooks/fitting/2D_fitting/plot_1_Rb2SO4_QMAT.ipynb index a84145ec3..6279386f9 100644 --- a/docs/notebooks/fitting/2D_fitting/plot_1_Rb2SO4_QMAT.ipynb +++ b/docs/notebooks/fitting/2D_fitting/plot_1_Rb2SO4_QMAT.ipynb @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/2D_fitting/plot_2_Coesite_DAS.ipynb b/docs/notebooks/fitting/2D_fitting/plot_2_Coesite_DAS.ipynb index a4c9c51fe..8650d53b0 100644 --- a/docs/notebooks/fitting/2D_fitting/plot_2_Coesite_DAS.ipynb +++ b/docs/notebooks/fitting/2D_fitting/plot_2_Coesite_DAS.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "shifts = [29, 39, 54.8, 51, 56] # in ppm\nCq = [6.1e6, 5.4e6, 5.5e6, 5.5e6, 5.1e6] # in Hz\neta = [0.1, 0.2, 0.15, 0.15, 0.3]\nabundance_ratio = [1, 1, 2, 2, 2]\nabundance = np.asarray(abundance_ratio) / 8 * 100 # in %\n\nspin_systems = single_site_system_generator(\n isotope=\"17O\",\n isotropic_chemical_shift=shifts,\n quadrupolar={\"Cq\": Cq, \"eta\": eta},\n abundance=abundance,\n)" + "shifts = [29, 39, 54.8, 51, 56] # in ppm\nCq = [6.1e6, 5.4e6, 5.5e6, 5.5e6, 5.1e6] # in Hz\neta = [0.1, 0.2, 0.15, 0.15, 0.3]\nabundance_ratio = [1, 1, 2, 2, 2]\nabundance = np.asarray(abundance_ratio) / 8 * 100 # in %\n\nspin_systems = single_site_system_generator(\n isotopes=\"17O\",\n isotropic_chemical_shifts=shifts,\n quadrupolar={\"Cq\": Cq, \"eta\": eta},\n abundance=abundance,\n)" ] }, { @@ -204,7 +204,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/2D_fitting/plot_3_RbNO3_MQMAS.ipynb b/docs/notebooks/fitting/2D_fitting/plot_3_RbNO3_MQMAS.ipynb index 1db008eba..4ffcec662 100644 --- a/docs/notebooks/fitting/2D_fitting/plot_3_RbNO3_MQMAS.ipynb +++ b/docs/notebooks/fitting/2D_fitting/plot_3_RbNO3_MQMAS.ipynb @@ -69,7 +69,7 @@ }, "outputs": [], "source": [ - "shifts = [-26.4, -28.5, -31.3] # in ppm\nCq = [1.7e6, 2.0e6, 1.7e6] # in Hz\neta = [0.2, 1.0, 0.6]\nabundance = [33.33, 33.33, 33.33] # in %\n\nspin_systems = single_site_system_generator(\n isotope=\"87Rb\",\n isotropic_chemical_shift=shifts,\n quadrupolar={\"Cq\": Cq, \"eta\": eta},\n abundance=abundance,\n)" + "shifts = [-26.4, -28.5, -31.3] # in ppm\nCq = [1.7e6, 2.0e6, 1.7e6] # in Hz\neta = [0.2, 1.0, 0.6]\nabundance = [33.33, 33.33, 33.33] # in %\n\nspin_systems = single_site_system_generator(\n isotopes=\"87Rb\",\n isotropic_chemical_shifts=shifts,\n quadrupolar={\"Cq\": Cq, \"eta\": eta},\n abundance=abundance,\n)" ] }, { @@ -197,7 +197,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.ipynb b/docs/notebooks/fitting/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.ipynb index 3af60eb4a..0b2306511 100644 --- a/docs/notebooks/fitting/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.ipynb +++ b/docs/notebooks/fitting/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.ipynb @@ -204,7 +204,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/docs/notebooks/fitting/2D_fitting/plot_5_CoCl2.2D2O_shifting-d.ipynb b/docs/notebooks/fitting/2D_fitting/plot_5_CoCl2.2D2O_shifting-d.ipynb index 2a4e973e2..a507c5375 100644 --- a/docs/notebooks/fitting/2D_fitting/plot_5_CoCl2.2D2O_shifting-d.ipynb +++ b/docs/notebooks/fitting/2D_fitting/plot_5_CoCl2.2D2O_shifting-d.ipynb @@ -204,7 +204,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.11" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/environment-dev.yml b/environment-dev.yml index 087f8ea71..abfb73254 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -14,7 +14,7 @@ dependencies: - pip - pip: - monty>=2.0.4 - - csdmpy>=0.4.1 + - csdmpy>=0.4 - pydantic>=1.0 - typing-extensions>=3.7 - flake8 diff --git a/environment.yml b/environment.yml index c86c6ce25..175c06249 100644 --- a/environment.yml +++ b/environment.yml @@ -16,7 +16,7 @@ dependencies: - joblib>=1.0.0 - pandas>=1.1.3 - monty>=2.0.4 - - csdmpy>=0.4.1 + - csdmpy>=0.4 - pydantic>=1.0 - typing-extensions>=3.7 - lmfit>=1.0.2 diff --git a/examples_source/1D_simulation(crystalline)/plot_0_Wollastonite.py b/examples_source/1D_simulation(crystalline)/plot_0_Wollastonite.py index ac5f0315d..ef80cc2e7 100644 --- a/examples_source/1D_simulation(crystalline)/plot_0_Wollastonite.py +++ b/examples_source/1D_simulation(crystalline)/plot_0_Wollastonite.py @@ -60,11 +60,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(4, 2)) -method.plot() -plt.show() - # %% # **Step 4:** Create the Simulator object and add the method and spin system objects. sim = Simulator() diff --git a/examples_source/1D_simulation(crystalline)/plot_1_PotassiumSulfate.py b/examples_source/1D_simulation(crystalline)/plot_1_PotassiumSulfate.py index 1d36e8aa2..da967b2e4 100644 --- a/examples_source/1D_simulation(crystalline)/plot_1_PotassiumSulfate.py +++ b/examples_source/1D_simulation(crystalline)/plot_1_PotassiumSulfate.py @@ -16,7 +16,7 @@ from mrsimulator.methods import BlochDecayCTSpectrum from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 2 # %% # **Step 1:** Create the spin system @@ -44,11 +44,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(4, 3)) -method.plot() -plt.show() - # %% # **Step 3:** Create the Simulator object and add method and spin system objects. sim = Simulator() diff --git a/examples_source/1D_simulation(crystalline)/plot_3_satellite_transition_sim.py b/examples_source/1D_simulation(crystalline)/plot_3_satellite_transition_sim.py index 69f37d669..e5e8b1941 100644 --- a/examples_source/1D_simulation(crystalline)/plot_3_satellite_transition_sim.py +++ b/examples_source/1D_simulation(crystalline)/plot_3_satellite_transition_sim.py @@ -59,7 +59,6 @@ # For illustrative purposes, let's look at the infinite speed spectrum from this # satellite transition. method = Method1D( - name="Inner Satellite Spectrum", channels=["27Al"], magnetic_flux_density=21.14, # in T rotor_frequency=1e9, # in Hz @@ -79,11 +78,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3)) -method.plot() -plt.show() - # %% # Create the Simulator object and add the method and the spin system object. sim = Simulator() @@ -111,7 +105,6 @@ # - :math:`|-1/2\rangle\rightarrow|-3/2\rangle` (:math:`P=-1, D=2`) # - :math:`|-3/2\rangle\rightarrow|-5/2\rangle` (:math:`P=-1, D=4`) method2 = Method1D( - name="Satellite Spectrum", channels=["27Al"], magnetic_flux_density=21.14, # in T rotor_frequency=1e9, # in Hz @@ -132,11 +125,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3)) -method2.plot() -plt.show() - # %% # Update the method object in the Simulator object. sim.methods[0] = method2 # add the method diff --git a/examples_source/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.py b/examples_source/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.py index 3008d85d6..bdc3e7383 100644 --- a/examples_source/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.py +++ b/examples_source/1D_simulation(crystalline)/plot_4_multi-quantum_spectrum.py @@ -13,7 +13,7 @@ from mrsimulator import Simulator, SpinSystem, Site from mrsimulator.methods import Method1D -# sphinx_gallery_thumbnail_number = 2 +# sphinx_gallery_thumbnail_number = 1 # %% # Create a single-site arbitrary spin system. @@ -29,16 +29,15 @@ # Selecting the triple-quantum transition # --------------------------------------- # -# For single-site spin-5/2 spin system, there are three triple-quantum transition +# For spin-site spin-5/2 spin system, there are three triple-quantum transition # # - :math:`|1/2\rangle\rightarrow|-5/2\rangle` (:math:`P=-3, D=6`) # - :math:`|3/2\rangle\rightarrow|-3/2\rangle` (:math:`P=-3, D=0`) # - :math:`|5/2\rangle\rightarrow|-1/2\rangle` (:math:`P=-3, D=-6`) # # To select one or more triple-quantum transitions, assign the respective value of P and -# D to the `transition_query`. Here, we select the symmetric triple-quantum transition. +# D to the `transition_query`. method = Method1D( - name="Arbitrary Transition Method", channels=["27Al"], magnetic_flux_density=21.14, # in T rotor_frequency=1e9, # in Hz @@ -56,11 +55,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3)) -method.plot() -plt.show() - # %% # Create the Simulator object and add the method and the spin system object. sim = Simulator() diff --git a/examples_source/1D_simulation(crystalline)/plot_9_custom_method_Hahnecho.py b/examples_source/1D_simulation(crystalline)/plot_9_custom_method_Hahnecho.py deleted file mode 100644 index 9ba96bee2..000000000 --- a/examples_source/1D_simulation(crystalline)/plot_9_custom_method_Hahnecho.py +++ /dev/null @@ -1,140 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -""" -Writing Custom methods (HahnEcho) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Writing custom methods using the Event objects. -""" -# %% -# Import the relevant modules -import matplotlib.pyplot as plt -import numpy as np - -from mrsimulator import Simulator, SpinSystem, Site, Coupling -from mrsimulator.methods import Method1D -from mrsimulator.method.event import MixingEvent, SpectralEvent - -# sphinx_gallery_thumbnail_number = 1 - -# %% -# For demonstration, we will create two spin systems, one with a single site and other -# with two spin 1/2 sites. - -S1 = Site( - isotope="1H", - isotropic_chemical_shift=10, # in ppm - shielding_symmetric={"zeta": -80, "eta": 0.25}, # zeta in ppm -) -S2 = Site(isotope="1H", isotropic_chemical_shift=-10) -S12 = Coupling( - site_index=[0, 1], isotropic_j=100, dipolar={"D": 2000, "eta": 0, "alpha": 0} -) - -spin_system_1 = SpinSystem(sites=[S1], label="Uncoupled system") -spin_system_2 = SpinSystem(sites=[S1, S2], couplings=[S12], label="Coupled system") - - -# %% -# **Create a custom method** -# -# Writing a custom method is simply specifying an appropriate list of event objects per -# spectral dimension. In this example, we are interested in a one-dimensional Hahnecho -# method, and we use the generic `Method1D` class as a template. For a Hahnecho, we will -# use two types of Event objects---SpectralEvent and MixingEvent. -# -# A SpectralEvent object is where we sample the frequency contributions. The net -# frequency along a given spectral dimension is a weighted average of the frequencies -# from all SpectralEvent objects within a given SpectralDimension, i.e., -# -# .. math:: -# f = \sum_j w_j \times \nu_j, -# -# where :math:`f` is the net averaged frequency along a spectral dimension, :math:`w_j` -# is the weight (the attribute ``fraction``), and :math:`\nu_j` is the frequency from -# the :math:`j^\text{th}` SpectralEvent. The index :math:`j` runs over all spectral -# events within a spectral dimension. -# -# In the case of a one-dimensional Hahnecho method, the frequency is equally averaged -# over two spectral events, corresponding to the symmetry pathway, -# -# .. math:: -# p = 1 \rightarrow -1. -# -# In the following code, we define the two SpectralEvent objects with fraction 0.5 and -# the transition_query on channel-1 of P=[1] and P=[-1], respectively. Notice, the value -# for the ``P`` attribute is a list. Here, it is a list with a single integer. The list -# notation, ``[1]``, implies that the query selects all transitions where exactly one -# spin is undergoing a :math:`p=+1` transition with the remaining spin at :math:`p=0`. -# A similar argument holds for ``[-1]`` query. By implementing query objects, we -# decouple the method from the spin system, i.e., once a method is defined, it can be -# used to simulate spectra from any given spin system. We will demonstrate this -# momentarily by simulating a Hahnecho spectrum from single and two-site spin systems. -# -# Besides the SpectralEvent, you may also notice a MixingEvent sandwitched in-between -# the two SpectralEvent. A MixingEvent does not directly contribute to the frequencies. -# As the name suggests, a mixing event is used for the mixing of transitions in a -# multi-event method such as HahnEcho. In the following code, we define a mixing query -# on channel-1 by setting the attributes ``tip_angle`` and ``phase`` to :math:`\pi` and -# 0, respectively. There two parameters are analogous to the pulse angle and phase. -hahn_echo = Method1D( - channels=["1H"], - magnetic_flux_density=9.4, # in T - spectral_dimensions=[ - { - "count": 512, - "spectral_width": 2e4, # in Hz - "events": [ - SpectralEvent(fraction=0.5, transition_query=[{"ch1": {"P": [1]}}]), - MixingEvent(mixing_query={"ch1": {"tip_angle": np.pi, "phase": 0}}), - SpectralEvent(fraction=0.5, transition_query=[{"ch1": {"P": [-1]}}]), - ], - } - ], -) - -# %% -# You may also visualize the method using the `plot` function. -# hahn_echo.plot() - -# %% -# As mentioned before, a method object is decoupled from the spin system object. Notice, -# when we get the transition pathways from this method for a single-site spin system, we -# get a single transition pathway. -print(hahn_echo.get_transition_pathways(spin_system_1)) - -# %% -# In the case of a homonuclear two-site spin 1/2 spin system, the same method returns -# four transition pathways. -print(hahn_echo.get_transition_pathways(spin_system_2)) - -# %% -# Create the Simulator object, add the method and spin system objects, and run the -# simulation. -sim = Simulator() -sim.spin_systems = [spin_system_1, spin_system_2] # add the spin systems -sim.methods = [hahn_echo] # add the method -sim.config.decompose_spectrum = "spin_system" - -sim.run() - -# %% -# The simulation from each spin system is stored as a dependent variable within the -# CSDM object. Use the `split` function to split the list of the dependent variables -# into a list of CSDM objects. -simulation_results = sim.methods[0].simulation.split() - -# The plot of the two simulations. -fig, ax = plt.subplots(1, 2, figsize=(8.0, 3.0), subplot_kw={"projection": "csdm"}) -for i in range(2): - ax[i].plot(simulation_results[i].real, color="black", linewidth=1) - ax[i].invert_xaxis() -plt.tight_layout() -plt.show() - -# %% -# Notice, in the single-site spin system, the hahn echo refocuses the isotropic chemical -# shits and chemical shift anisotropies. The end result is a resonance at zero -# frequency. In the case of the two homonuclear spin 1/2 coupled spin system, the Hahn -# echo refocuses the isotropic chemical shits and chemical shift anisotropies, but not -# the dipolar and `J` couplings. diff --git a/examples_source/1D_simulation(macro_amorphous)/plot_3_amorphous_like.py b/examples_source/1D_simulation(macro_amorphous)/plot_3_amorphous_like.py index e7896d82f..b1f6c473c 100644 --- a/examples_source/1D_simulation(macro_amorphous)/plot_3_amorphous_like.py +++ b/examples_source/1D_simulation(macro_amorphous)/plot_3_amorphous_like.py @@ -89,8 +89,8 @@ # are the array of tensor parameter coordinates, and ``pdf`` is the array of the # corresponding amplitudes. spin_systems = single_site_system_generator( - isotope="29Si", - isotropic_chemical_shift=iso, + isotopes="29Si", + isotropic_chemical_shifts=iso, shielding_symmetric={"zeta": zeta, "eta": eta}, abundance=pdf, ) diff --git a/examples_source/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.py b/examples_source/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.py index 4c0e7d985..aea4c294f 100644 --- a/examples_source/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.py +++ b/examples_source/1D_simulation(macro_amorphous)/plot_4_amorphous_like_quad.py @@ -71,8 +71,8 @@ # Use the :func:`~mrsimulator.utils.collection.single_site_system_generator` utility # function to generate single-site spin systems. spin_systems = single_site_system_generator( - isotope="27Al", - isotropic_chemical_shift=iso, + isotopes="27Al", + isotropic_chemical_shifts=iso, quadrupolar={"Cq": Cq * 1e6, "eta": eta}, # Cq in Hz abundance=pdf, ) diff --git a/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py b/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py index 2ef03f1c7..d07a39c66 100644 --- a/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py +++ b/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py @@ -58,7 +58,7 @@ # :math:`\eta` parameters, use the # :func:`~mrsimulator.utils.collection.single_site_system_generator` utility function. systems = single_site_system_generator( - isotope="13C", shielding_symmetric={"zeta": z_dist, "eta": e_dist}, abundance=amp + isotopes="13C", shielding_symmetric={"zeta": z_dist, "eta": e_dist}, abundance=amp ) # %% @@ -107,7 +107,7 @@ # # Create the spin systems. systems = single_site_system_generator( - isotope="71Ga", quadrupolar={"Cq": cq_dist * 1e6, "eta": e_dist}, abundance=amp + isotopes="71Ga", quadrupolar={"Cq": cq_dist * 1e6, "eta": e_dist}, abundance=amp ) # %% diff --git a/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py b/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py index 4feb075c8..9b9c0056b 100644 --- a/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py +++ b/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py @@ -53,7 +53,7 @@ # # Create the spin systems from the above :math:`\zeta` and :math:`\eta` parameters. systems = single_site_system_generator( - isotope="13C", shielding_symmetric={"zeta": z_dist, "eta": e_dist}, abundance=amp + isotopes="13C", shielding_symmetric={"zeta": z_dist, "eta": e_dist}, abundance=amp ) print(len(systems)) @@ -107,7 +107,7 @@ # **Static spectrum** # Create the spin systems. systems = single_site_system_generator( - isotope="71Ga", quadrupolar={"Cq": cq_dist * 1e6, "eta": e_dist}, abundance=amp + isotopes="71Ga", quadrupolar={"Cq": cq_dist * 1e6, "eta": e_dist}, abundance=amp ) # %% diff --git a/examples_source/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.py b/examples_source/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.py index a8b026252..3d9a6226c 100644 --- a/examples_source/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.py +++ b/examples_source/2D_simulation(crystalline)/plot_0_MQMAS_RbNO3.py @@ -17,7 +17,7 @@ from mrsimulator.methods import ThreeQ_VAS from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 2 # %% # Generate the site and spin system objects. @@ -62,11 +62,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -method.plot() -plt.show() - # %% # Create the Simulator object, add the method and spin system objects, and # run the simulation. diff --git a/examples_source/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.py b/examples_source/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.py index 623dd52da..ab6ac7460 100644 --- a/examples_source/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.py +++ b/examples_source/2D_simulation(crystalline)/plot_1_STMAS_RbNO3.py @@ -15,7 +15,7 @@ from mrsimulator.methods import ST1_VAS from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 2 # %% # Generate the site and spin system objects. @@ -70,11 +70,6 @@ ) ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -method[0].plot() -plt.show() - # %% # Create the Simulator object, add the method and spin system objects, and # run the simulation. diff --git a/examples_source/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.py b/examples_source/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.py index 3b0138a0c..20576d80e 100644 --- a/examples_source/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.py +++ b/examples_source/2D_simulation(crystalline)/plot_2_SAS_Rb2SO4.py @@ -16,7 +16,7 @@ from mrsimulator.methods import Method2D from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 2 # %% # Generate the site and spin system objects. @@ -39,7 +39,6 @@ # method parameters, as shown below. Note, the Method2D method simulates an infinite # spinning speed spectrum. sas = Method2D( - name="Switched Angle Spinning", channels=["87Rb"], magnetic_flux_density=9.4, # in T spectral_dimensions=[ @@ -70,11 +69,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -sas.plot() -plt.show() - # %% # Create the Simulator object, add the method and spin system objects, and # run the simulation. diff --git a/examples_source/2D_simulation(crystalline)/plot_4_DAS_Coesite.py b/examples_source/2D_simulation(crystalline)/plot_4_DAS_Coesite.py index 169cd3c73..cab70503c 100644 --- a/examples_source/2D_simulation(crystalline)/plot_4_DAS_Coesite.py +++ b/examples_source/2D_simulation(crystalline)/plot_4_DAS_Coesite.py @@ -16,7 +16,7 @@ from mrsimulator.methods import Method2D from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 2 # %% # Create the Simulator object and load the spin systems database or url address. @@ -31,7 +31,6 @@ # method parameters, as shown below. Note, the Method2D method simulates an infinite # spinning speed spectrum. das = Method2D( - name="Dynamic Angle Spinning", channels=["17O"], magnetic_flux_density=11.74, # in T spectral_dimensions=[ @@ -68,12 +67,7 @@ }, ], ) -sim.methods = [das] # add the method - -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -das.plot() -plt.show() +sim.methods = [das] # add the method. # %% # Run the simulation diff --git a/examples_source/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.py b/examples_source/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.py index 4cae98ce5..621da35b2 100644 --- a/examples_source/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.py +++ b/examples_source/2D_simulation(crystalline)/plot_5_COASTER_Rb2CrO4.py @@ -18,7 +18,7 @@ from mrsimulator.methods import Method2D from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 2 # %% # Generate the site and spin system objects. @@ -41,7 +41,6 @@ # the method parameters, as shown below. Note, the Method2D method simulates an infinite # spinning speed spectrum. coaster = Method2D( - name="COASTER", channels=["87Rb"], magnetic_flux_density=9.4, # in T rotor_angle=70.12 * 3.14159 / 180, # in rads @@ -64,11 +63,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -coaster.plot() -plt.show() - # %% # Create the Simulator object, add the method and spin system objects, and # run the simulation. diff --git a/examples_source/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.py b/examples_source/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.py index d6de00e2b..694c8ecbf 100644 --- a/examples_source/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.py +++ b/examples_source/2D_simulation(crystalline)/plot_6_PASS_itraconazole_drug.py @@ -17,7 +17,7 @@ from mrsimulator.methods import SSB2D from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 2 +# sphinx_gallery_thumbnail_number = 1 # %% # There are 41 :math:`^{13}\text{C}` single-site spin systems partially describing the @@ -53,12 +53,6 @@ ) sim.methods = [PASS] # add the method. -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -PASS.plot() -plt.show() - -# %% # For 2D spinning sideband simulation, set the number of spinning sidebands in the # Simulator.config object to `spectral_width/rotor_frequency` along the sideband # dimension. diff --git a/examples_source/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.py b/examples_source/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.py index ed3795159..d5eeead0b 100644 --- a/examples_source/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.py +++ b/examples_source/2D_simulation(crystalline)/plot_7_QSSB_Rb2SO4.py @@ -16,7 +16,7 @@ from mrsimulator import Simulator, SpinSystem, Site from mrsimulator.methods import SSB2D -# sphinx_gallery_thumbnail_number = 2 +# sphinx_gallery_thumbnail_number = 1 # %% # Generate the site and spin system objects. @@ -58,11 +58,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -qmat.plot() -plt.show() - # %% # Create the Simulator object, add the method and spin system objects, and # run the simulation. diff --git a/examples_source/2D_simulation(crystalline)/plot_8_MAF.py b/examples_source/2D_simulation(crystalline)/plot_8_MAF.py index 840ae5bbb..d0787a7d3 100644 --- a/examples_source/2D_simulation(crystalline)/plot_8_MAF.py +++ b/examples_source/2D_simulation(crystalline)/plot_8_MAF.py @@ -17,7 +17,7 @@ from mrsimulator.methods import Method2D from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 2 +# sphinx_gallery_thumbnail_number = 1 # %% # Create the sites and spin systems @@ -46,7 +46,6 @@ # spectrum by customizing the method parameters, as shown below. Note, the Method2D # method simulates an infinite spinning speed spectrum. maf = Method2D( - name="Magic Angle Spinning", channels=["29Si"], magnetic_flux_density=14.1, # in T spectral_dimensions=[ @@ -77,11 +76,6 @@ affine_matrix=[[1, -1], [0, 1]], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 3.5)) -maf.plot() -plt.show() - # %% # Create the Simulator object, add the method and spin system objects, and run the # simulation. diff --git a/examples_source/2D_simulation(crystalline)/plot_9_shifting-d.py b/examples_source/2D_simulation(crystalline)/plot_9_shifting-d.py index b96a73939..cf84e1671 100644 --- a/examples_source/2D_simulation(crystalline)/plot_9_shifting-d.py +++ b/examples_source/2D_simulation(crystalline)/plot_9_shifting-d.py @@ -18,7 +18,7 @@ from mrsimulator.methods import Method2D from mrsimulator import signal_processing as sp -# sphinx_gallery_thumbnail_number = 3 +# sphinx_gallery_thumbnail_number = 2 # %% # Generate the site and spin system objects. @@ -108,7 +108,6 @@ # the first-order shielding with zeroth and second-rank tensor contributions, # respectively. See :ref:`freq_contrib_api` for details. shifting_d = Method2D( - name="Shifting-d", channels=["2H"], magnetic_flux_density=9.395, # in T spectral_dimensions=[ @@ -140,11 +139,6 @@ ], ) -# A graphical representation of the method object. -plt.figure(figsize=(5, 2.5)) -shifting_d.plot() -plt.show() - # %% # Create the Simulator object, add the method and spin system objects, and # run the simulation. diff --git a/examples_source/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.py b/examples_source/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.py index c9fee013f..b6f734ec0 100644 --- a/examples_source/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.py +++ b/examples_source/2D_simulation(macro_amorphous)/plot_0_crystalline_disorder.py @@ -90,8 +90,8 @@ def get_prob_dist(iso, Cq, eta, eps, cov): # ---------------- # Generate spin systems from the above probability distribution. spin_systems = single_site_system_generator( - isotope="87Rb", - isotropic_chemical_shift=iso, + isotopes="87Rb", + isotropic_chemical_shifts=iso, quadrupolar={"Cq": Cq * 1e6, "eta": eta}, # Cq in Hz abundance=pdf, ) diff --git a/examples_source/2D_simulation(macro_amorphous)/plot_1_I=2.5.py b/examples_source/2D_simulation(macro_amorphous)/plot_1_I=2.5.py index 497deceb0..873f696be 100644 --- a/examples_source/2D_simulation(macro_amorphous)/plot_1_I=2.5.py +++ b/examples_source/2D_simulation(macro_amorphous)/plot_1_I=2.5.py @@ -77,8 +77,8 @@ # :func:`~mrsimulator.utils.collection.single_site_system_generator` utility function to # generate single-site spin systems. spin_systems = single_site_system_generator( - isotope="27Al", - isotropic_chemical_shift=iso, + isotopes="27Al", + isotropic_chemical_shifts=iso, quadrupolar={"Cq": Cq * 1e6, "eta": eta}, # Cq in Hz abundance=pdf, ) diff --git a/fitting_source/1D_fitting/plot_2_13C_glycine_1940Hz.py b/fitting_source/1D_fitting/plot_2_13C_glycine_1940Hz.py index 017e1d116..cbf5d58b8 100644 --- a/fitting_source/1D_fitting/plot_2_13C_glycine_1940Hz.py +++ b/fitting_source/1D_fitting/plot_2_13C_glycine_1940Hz.py @@ -29,7 +29,7 @@ experiment = cp.load(host + filename) # standard deviation of noise from the dataset -sigma = 3.890551 +sigma = 3.822249 # For spectral fitting, we only focus on the real part of the complex dataset experiment = experiment.real diff --git a/fitting_source/1D_fitting/plot_2_13C_glycine_5000Hz.py b/fitting_source/1D_fitting/plot_2_13C_glycine_5000Hz.py index 500ebb21c..5aa739ae5 100644 --- a/fitting_source/1D_fitting/plot_2_13C_glycine_5000Hz.py +++ b/fitting_source/1D_fitting/plot_2_13C_glycine_5000Hz.py @@ -29,7 +29,7 @@ experiment = cp.load(host + filename) # standard deviation of noise from the dataset -sigma = 1.97637 +sigma = 3.822249 # For spectral fitting, we only focus on the real part of the complex dataset experiment = experiment.real diff --git a/fitting_source/1D_fitting/plot_2_13C_glycine_960Hz.py b/fitting_source/1D_fitting/plot_2_13C_glycine_960Hz.py index de035b2a8..f33f3a260 100644 --- a/fitting_source/1D_fitting/plot_2_13C_glycine_960Hz.py +++ b/fitting_source/1D_fitting/plot_2_13C_glycine_960Hz.py @@ -29,7 +29,7 @@ experiment = cp.load(host + filename) # standard deviation of noise from the dataset -sigma = 3.982936 +sigma = 3.822249 # For spectral fitting, we only focus on the real part of the complex dataset experiment = experiment.real diff --git a/fitting_source/1D_fitting/plot_2_13C_glycine_multi_spectra_fit.py b/fitting_source/1D_fitting/plot_2_13C_glycine_multi_spectra_fit.py deleted file mode 100644 index f39925a74..000000000 --- a/fitting_source/1D_fitting/plot_2_13C_glycine_multi_spectra_fit.py +++ /dev/null @@ -1,225 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -""" -¹³C MAS NMR of Glycine (CSA) multi-spectra fit -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -""" -# %% -# The following is a multi-data least-squares fitting example of -# :math:`^{13}\text{C}` MAS NMR spectrum of Glycine spinning at 5 kHz, 1.94 kHz, and -# 960 Hz. Before trying multi-data fitting, we recommend that you first try individual -# fits. The experimental datasets are part of DMFIT [#f1]_ examples. -import csdmpy as cp -import matplotlib.pyplot as plt -from lmfit import Minimizer - -from mrsimulator import Simulator, SpinSystem, Site -from mrsimulator.methods import BlochDecaySpectrum -from mrsimulator import signal_processing as sp -from mrsimulator.utils import spectral_fitting as sf -from mrsimulator.utils import get_spectral_dimensions - -# sphinx_gallery_thumbnail_number = 3 - -# %% -# Import the datasets -# ------------------- -# Import the datasets and assign the standard deviation of noise for each dataset. Here, -# ``sigma1``, ``sigma2``, and ``sigma3`` are the noise standard deviation for the -# dataset acquired at 5 kHz, 1.94 kHz, and 960 Hz spinning speeds, respectively. -host = "https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/" -filename1 = "13C MAS 5000Hz - Glycine.csdf" -filename2 = "13C MAS 1940Hz - Glycine.csdf" -filename3 = "13C MAS 960Hz - Glycine.csdf" - -experiment1 = cp.load(host + filename1).real -experiment2 = cp.load(host + filename2).real -experiment3 = cp.load(host + filename3).real -experiments = [experiment1, experiment2, experiment3] - -# standard deviation of noise from the dataset -sigma1 = 1.97637 -sigma2 = 3.822249 -sigma3 = 3.982936 -sigmas = [sigma1, sigma2, sigma3] - -# Convert the coordinates along each dimension from Hz to ppm for each dataset -fig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={"projection": "csdm"}) -for i, experiment in enumerate(experiments): - _ = [item.to("ppm", "nmr_frequency_ratio") for item in experiment.dimensions] - - # plot of the dataset. - ax[i].plot(experiment, color="black", linewidth=0.5, label="Experiment") - ax[i].set_xlim(280, -10) - ax[i].grid() -plt.tight_layout() -plt.show() - -# %% -# Create a fitting model -# ---------------------- -# **Spin System**: The objective of a multi-data fitting is to optimize the spin -# system parameters using multiple datasets. In this example, we create two single-site -# spin systems, which are then shared by three method objects. -C1 = Site( - isotope="13C", - isotropic_chemical_shift=176.0, # in ppm - shielding_symmetric={"zeta": 60, "eta": 0.6}, # zeta in Hz -) -C2 = Site( - isotope="13C", - isotropic_chemical_shift=43.0, # in ppm - shielding_symmetric={"zeta": 30, "eta": 0.5}, # zeta in Hz -) - -spin_systems = [SpinSystem(sites=[C1], name="C1"), SpinSystem(sites=[C2], name="C2")] - -# %% -# **Method**: Create the three MAS method objects with respective MAS spinning speeds. - -# Get the spectral dimension parameters from the respective experiment and setup the -# corresponding method. - -# Method for dataset 1 -spectral_dims1 = get_spectral_dimensions(experiment1) -MAS1 = BlochDecaySpectrum( - channels=["13C"], - magnetic_flux_density=7.05, # in T - rotor_frequency=5000, # in Hz - spectral_dimensions=spectral_dims1, - experiment=experiment1, # add experimental dataset 1 -) - -# Method for dataset 2 -spectral_dims2 = get_spectral_dimensions(experiment2) -MAS2 = BlochDecaySpectrum( - channels=["13C"], - magnetic_flux_density=7.05, # in T - rotor_frequency=1940, # in Hz - spectral_dimensions=spectral_dims2, - experiment=experiment2, # add experimental dataset 2 -) - -# Method for dataset 3 -spectral_dims3 = get_spectral_dimensions(experiment3) -MAS3 = BlochDecaySpectrum( - channels=["13C"], - magnetic_flux_density=7.05, # in T - rotor_frequency=960, # in Hz - spectral_dimensions=spectral_dims3, - experiment=experiment3, # add experimental dataset 3 -) - -# %% -# **Guess Model Spectrum** - -# Simulation -# ---------- -# Add the spin systems and the three methods to the simulator object. -sim = Simulator(spin_systems=spin_systems, methods=[MAS1, MAS2, MAS3]) -sim.config.decompose_spectrum = "spin_system" -sim.run() - -# Post Simulation Processing -# -------------------------- -# Add signal processing to simulation dataset from the three methods. - -# Processor for dataset 1 -processor1 = sp.SignalProcessor( - operations=[ - sp.IFFT(), - sp.apodization.Exponential(FWHM="20 Hz", dv_index=0), # spin system 0 - sp.apodization.Exponential(FWHM="200 Hz", dv_index=1), # spin system 1 - sp.FFT(), - sp.Scale(factor=10), # dataset is scaled independently using scale factor. - ] -) - -# Processor for dataset 2 -processor2 = sp.SignalProcessor( - operations=[ - sp.IFFT(), - sp.apodization.Exponential(FWHM="30 Hz", dv_index=0), # spin system 0 - sp.apodization.Exponential(FWHM="300 Hz", dv_index=1), # spin system 1 - sp.FFT(), - sp.Scale(factor=100), # dataset is scaled independently using scale factor. - ] -) - -# Processor for dataset 3 -processor3 = sp.SignalProcessor( - operations=[ - sp.IFFT(), - sp.apodization.Exponential(FWHM="10 Hz", dv_index=0), # spin system 0 - sp.apodization.Exponential(FWHM="150 Hz", dv_index=1), # spin system 1 - sp.FFT(), - sp.Scale(factor=50), # dataset is scaled independently using scale factor. - ] -) -processors = [processor1, processor2, processor3] - -processed_data = [] -for i, proc in enumerate(processors): - processed_data.append(proc.apply_operations(data=sim.methods[i].simulation).real) - - -# Plot of the guess Spectrum -# -------------------------- - -fig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={"projection": "csdm"}) -for i, exp_data in enumerate(experiments): - ax[i].plot(exp_data, color="black", linewidth=0.5, label="Experiment") - ax[i].plot(processed_data[i], linewidth=2, alpha=0.6) - ax[i].set_xlim(280, -10) - ax[i].grid() - -plt.legend() -plt.tight_layout() -plt.show() - - -# %% -# Least-squares minimization with LMFIT -# ------------------------------------- -# Use the :func:`~mrsimulator.utils.spectral_fitting.make_LMFIT_params` for a quick -# setup of the fitting parameters. Note, the first two arguments of this function is -# the simulator object and a list of signal processing objects, ``processors``. The -# fitting parameters corresponding to the signal processor objects are generated using -# ``SP_i_operation_j_FunctionName_FunctionArg``, where *i* is the *ith* signal -# processor within the list, *j* is the operation index of the *ith* processor, and -# *FunctionName* and *FunctionArg* are the operation function name and function -# argument, respectively. -params = sf.make_LMFIT_params(sim, processors, include={"rotor_frequency"}) -print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"])) - -# %% -# **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processors, sigmas)) -result = minner.minimize() -result - -# %% -# The best fit solution -# --------------------- -all_best_fit = sf.bestfit(sim, processors) # a list of best fit simulations -all_residuals = sf.residuals(sim, processors) # a list of residuals - -# Plot the spectrum -fig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={"projection": "csdm"}) -for i, proc in enumerate(processors): - ax[i].plot(experiments[i], color="black", linewidth=0.5, label="Experiment") - ax[i].plot(all_residuals[i], color="gray", linewidth=0.5, label="Residual") - ax[i].plot(all_best_fit[i], linewidth=2, alpha=0.6) - ax[i].set_xlim(280, -10) - ax[i].grid() - -plt.legend() -plt.tight_layout() -plt.show() - -# %% -# -# .. [#f1] D.Massiot, F.Fayon, M.Capron, I.King, S.Le Calvé, B.Alonso, J.O.Durand, -# B.Bujoli, Z.Gan, G.Hoatson, 'Modelling one and two-dimensional solid-state NMR -# spectra.', Magn. Reson. Chem. **40** 70-76 (2002) -# `DOI: 10.1002/mrc.984 `_ diff --git a/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py b/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py index 5c9358724..b7a9d4ccc 100644 --- a/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py +++ b/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py @@ -67,8 +67,8 @@ eta = [0.8, 0.4, 0.9, 0.3, 0.0, 0.0] spin_systems = single_site_system_generator( - isotope="13C", - isotropic_chemical_shift=shifts, + isotopes="13C", + isotropic_chemical_shifts=shifts, shielding_symmetric={"zeta": zeta, "eta": eta}, abundance=100 / 6, ) diff --git a/fitting_source/2D_fitting/plot_2_Coesite_DAS.py b/fitting_source/2D_fitting/plot_2_Coesite_DAS.py index 2d03c514a..f8a769412 100644 --- a/fitting_source/2D_fitting/plot_2_Coesite_DAS.py +++ b/fitting_source/2D_fitting/plot_2_Coesite_DAS.py @@ -66,8 +66,8 @@ abundance = np.asarray(abundance_ratio) / 8 * 100 # in % spin_systems = single_site_system_generator( - isotope="17O", - isotropic_chemical_shift=shifts, + isotopes="17O", + isotropic_chemical_shifts=shifts, quadrupolar={"Cq": Cq, "eta": eta}, abundance=abundance, ) diff --git a/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py b/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py index b221b3ca9..9d8390ab1 100644 --- a/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py +++ b/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py @@ -63,8 +63,8 @@ abundance = [33.33, 33.33, 33.33] # in % spin_systems = single_site_system_generator( - isotope="87Rb", - isotropic_chemical_shift=shifts, + isotopes="87Rb", + isotropic_chemical_shifts=shifts, quadrupolar={"Cq": Cq, "eta": eta}, abundance=abundance, ) diff --git a/requirements-dev.txt b/requirements-dev.txt index 06fe6cb63..377f57117 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,6 +1,6 @@ numpy>=1.17 matplotlib>=3.3.3 -csdmpy>=0.4.1 +csdmpy>=0.4 pydantic>=1.0 monty>=2.0.4 typing-extensions>=3.7 diff --git a/requirements.txt b/requirements.txt index ffcf220ef..84b1c515e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ numpy>=1.17 matplotlib>=3.3.3 -csdmpy>=0.4.1 +csdmpy>=0.4 pydantic>=1.0 monty>=2.0.4 typing-extensions>=3.7 diff --git a/setup.py b/setup.py index fe6dcb388..2468c8536 100644 --- a/setup.py +++ b/setup.py @@ -200,11 +200,11 @@ def get_location(self, dict_info): getattr(self, item).extend(dict_info[item]) def message(self, lib_centos, lib_ubuntu): - print(f"Warning: {lib_ubuntu} might not be installed. See below.\n") stat = f"yum install {lib_centos}" print(f"For CentOS users\n\t{stat}") stat = f"sudo apt-get update\n\tsudo apt-get install {lib_ubuntu}" print(f"For Ubuntu users\n\t{stat}") + sys.exit(1) class MacOSSetup(Setup): @@ -421,7 +421,7 @@ def fftw_info(self): setup_requires=["numpy>=1.17"], install_requires=[ "numpy>=1.17", - "csdmpy>=0.4.1", + "csdmpy>=0.4", "pydantic>=1.0", "monty>=2.0.4", "typing-extensions>=3.7", diff --git a/src/c_lib/base/base_model.pxd b/src/c_lib/base/base_model.pxd index b04daa560..9358f20a4 100644 --- a/src/c_lib/base/base_model.pxd +++ b/src/c_lib/base/base_model.pxd @@ -145,7 +145,6 @@ cdef extern from "simulation.h": site_struct *sites, coupling_struct *couplings, float *transition_pathway, # Pointer to a list of transitions. - double transition_pathway_weight, # The weight of transition pathway. int n_dimension, # the number of dimensions. MRS_dimension *dimensions, # the dimensions within method. MRS_fftw_scheme *fftw_scheme, # the fftw scheme diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index d4a97ee8a..2873760ac 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -68,7 +68,6 @@ def one_d_spectrum(method, # transitions of the observed spin cdef int transition_increment cdef ndarray[float, ndim=1] transition_pathway_c - cdef ndarray[double, ndim=1] transition_pathway_weight_c cdef int number_of_transitions # transition_pathway_c = np.asarray([-0.5, 0.5]).ravel() # number_of_transitions = int(transition_pathway_c.size/2) @@ -98,7 +97,6 @@ def one_d_spectrum(method, max_n_sidebands = number_of_sidebands total_n_points = 1 - cdef int n_ev cdef ndarray[int] n_event cdef ndarray[double] magnetic_flux_density_in_T, frac cdef ndarray[double] srfiH @@ -119,34 +117,31 @@ def one_d_spectrum(method, prev_n_sidebands = 0 for i, dim in enumerate(method.spectral_dimensions): - n_ev = 0 for event in dim.events: - if event.__class__.__name__ != "MixingEvent": - freq_contrib = np.append(freq_contrib, event._freq_contrib_flags()) - if event.rotor_frequency < 1.0e-3: - rotor_frequency_in_Hz = 1.0e9 - rotor_angle_in_rad = 0.0 - number_of_sidebands = 1 - if prev_n_sidebands == 0: prev_n_sidebands = 1 - else: - rotor_frequency_in_Hz = event.rotor_frequency - rotor_angle_in_rad = event.rotor_angle - if prev_n_sidebands == 0: prev_n_sidebands = number_of_sidebands - - if prev_n_sidebands != number_of_sidebands: - raise ValueError( - ( - 'The library does not support spectral dimensions containing ' - 'both zero and non-zero rotor frequencies. Consider using a ' - 'smaller value instead of zero.' - ) + freq_contrib = np.append(freq_contrib, event._freq_contrib_flags()) + if event.rotor_frequency < 1.0e-3: + rotor_frequency_in_Hz = 1.0e9 + rotor_angle_in_rad = 0.0 + number_of_sidebands = 1 + if prev_n_sidebands == 0: prev_n_sidebands = 1 + else: + rotor_frequency_in_Hz = event.rotor_frequency + rotor_angle_in_rad = event.rotor_angle + if prev_n_sidebands == 0: prev_n_sidebands = number_of_sidebands + + if prev_n_sidebands != number_of_sidebands: + raise ValueError( + ( + 'The library does not support spectral dimensions containing ' + 'both zero and non-zero rotor frequencies. Consider using a ' + 'smaller value instead of zero.' ) + ) - fr.append(event.fraction) # fraction - Bo.append(event.magnetic_flux_density) # in T - vr.append(rotor_frequency_in_Hz) # in Hz - th.append(rotor_angle_in_rad) # in rad - n_ev +=1 + fr.append(event.fraction) # fraction + Bo.append(event.magnetic_flux_density) # in T + vr.append(rotor_frequency_in_Hz) # in Hz + th.append(rotor_angle_in_rad) # in rad total_n_points *= dim.count @@ -154,7 +149,7 @@ def one_d_spectrum(method, offset = dim.spectral_width / 2.0 coordinates_offset.append(-dim.reference_offset * factor - offset) increment.append(dim.spectral_width / dim.count) - event_i.append(n_ev) + event_i.append(len(dim.events)) if dim.origin_offset is None: dim.origin_offset = np.abs(Bo[0] * gyromagnetic_ratio * 1e6) @@ -168,15 +163,6 @@ def one_d_spectrum(method, coord_off = np.asarray(coordinates_offset, dtype=np.float64) n_event = np.asarray(event_i, dtype=np.int32) - # print('n_events', n_event) - # print('frac', frac) - # print('magnetic_flux_density_in_T', magnetic_flux_density_in_T) - # print('srfiH', srfiH) - # print('rair', rair) - # print('cnt', cnt) - # print('incre', incre) - # print('coord_off', coord_off) - # create spectral_dimensions dimensions = clib.MRS_create_dimensions(the_averaging_scheme, &cnt[0], &coord_off[0], &incre[0], &frac[0], &magnetic_flux_density_in_T[0], @@ -454,16 +440,11 @@ def one_d_spectrum(method, if number_of_sites != p_number_of_sites and isotopes != p_isotopes: transition_pathway = spin_sys.transition_pathways if transition_pathway is None: - segments, weights = method._get_transition_pathway_and_weights_np(spin_sys) - transition_pathway = np.asarray(segments) + transition_pathway = np.asarray(method._get_transition_pathways_np(spin_sys)) transition_pathway_c = np.asarray(transition_pathway, dtype=np.float32).ravel() - transition_pathway_weight_c = np.asarray(weights.real, dtype=np.float64).ravel() else: - # convert transition objects to list - weight = [item.weight.real for item in transition_pathway] - transition_pathway_weight_c = np.asarray(weight, dtype=np.float64).ravel() - transition_pathway = np.asarray(transition_pathway) + # convert transition objects to list lst = [item.tolist() for item in transition_pathway.ravel()] transition_pathway_c = np.asarray(lst, dtype=np.float32).ravel() @@ -485,9 +466,6 @@ def one_d_spectrum(method, # transition_increment = 2*number_of_sites # number_of_transitions = int((transition_pathway_c.size)/transition_increment) - # print('pathway', transition_pathway_c) - # print('weight', transition_pathway_weight_c) - # print('pathway_count, inc', pathway_count, pathway_increment) for trans__ in range(pathway_count): clib.__mrsimulator_core( # spectrum information and related amplitude @@ -495,7 +473,6 @@ def one_d_spectrum(method, &sites_c, &couplings_c, &transition_pathway_c[pathway_increment*trans__], - transition_pathway_weight_c[trans__], n_dimension, # The total number of spectroscopic dimensions. dimensions, # Pointer to MRS_dimension structure the_fftw_scheme, # Pointer to the fftw scheme. @@ -560,67 +537,12 @@ def get_zeeman_states(sys): @cython.wraparound(False) def transition_connect_factor(float l, float m1_f, float m1_i, float m2_f, float m2_i, double theta, double phi): - """Evaluate the probability of connecting two transitions driven by an external rf - pulse of phase phi and tip_angle theta. The connected transitions are - | m1_f >< m1_i | --> | m2_f > < m2_i |. - - Args: - float l: The angular momentum quantum number of the spin involved in the transition. - float m1_f Final quantum number of the starting transition. - float m1_i Initial quantum number of the starting transition. - float m2_f Final quantum number of the connecting transition. - float m2_i Initial quantum number of the connecting transition. - float theta The tip-angle of the rf pulse. - float phi The phase of the rf pulse. - - Return: A complex amplitude. - """ - cdef ndarray[double] factor = np.asarray([1, 0], dtype=np.float64) + cdef ndarray[double] factor = np.zeros(2, dtype=np.float64) clib.transition_connect_factor(l, m1_f, m1_i, m2_f, m2_i, theta, phi, &factor[0]) factor = np.around(factor, decimals=12) return complex(factor[0], factor[1]) -@cython.profile(False) -@cython.boundscheck(False) -@cython.wraparound(False) -def calculate_transition_connect_weight( - ndarray[float, ndim=2] trans1, - ndarray[float, ndim=2] trans2, - ndarray[float, ndim=1] spin, - ndarray[double, ndim=1] theta, - ndarray[double, ndim=1] phi - ): - """Evaluate the probability of connecting two transitions driven by an external rf - pulse of phase phi and tip_angle theta. The connected transitions are - | m1_f >< m1_i | --> | m2_f > < m2_i |. - - Args: - float l: The angular momentum quantum number of the spin involved in the transition. - float m1_f Final quantum number of the starting transition. - float m1_i Initial quantum number of the starting transition. - float m2_f Final quantum number of the connecting transition. - float m2_i Initial quantum number of the connecting transition. - float theta The tip-angle of the rf pulse. - float phi The phase of the rf pulse. - - Return: A complex amplitude. - """ - cdef ndarray[double] factor = np.asarray([1, 0], dtype=np.float64) - cdef int i, n_sites = spin.size - cdef float m1_f, m1_i, m2_f, m2_i - for i in range(n_sites): - m1_f = trans1[1][i] # starting transition final state - m1_i = trans1[0][i] # starting transition initial state - m2_f = trans2[1][i] # landing transition final state - m2_i = trans2[0][i] # landing transition initial state - - clib.transition_connect_factor( - spin[i], m1_f, m1_i, m2_f, m2_i, theta[i], phi[i], &factor[0] - ) - return complex(factor[0], factor[1]) - - # @cython.profile(False) # @cython.boundscheck(False) # @cython.wraparound(False) diff --git a/src/c_lib/include/angular_momentum/wigner_element.h b/src/c_lib/include/angular_momentum/wigner_element.h index 1d2b1d1e1..edd479244 100644 --- a/src/c_lib/include/angular_momentum/wigner_element.h +++ b/src/c_lib/include/angular_momentum/wigner_element.h @@ -27,11 +27,11 @@ extern double wigner_d_element(const float l, const float m1, const float m2, * rf pulse. The connected transitions are |m1_f >< m1_i | --> | m2_f > < m2_i |. * * @param l The angular momentum quantum number of the spin involved in the transition. - * @param m1_f Final quantum number of the starting transition. - * @param m1_i Initial quantum number of the starting transition. - * @param m2_f Final quantum number of the connecting transition. - * @param m2_i Initial quantum number of the connecting transition. - * @param theta The tip-angle of the rf pulse. + * @param m1_f Final quantum number of the starting tranition. + * @param m1_i Initial quantum number of the starting tranition. + * @param m2_f Final quantum number of the connecting tranition. + * @param m2_i Initial quantum number of the connecting tranition. + * @param theta The tip-anlge of the rf pulse. * @param phi The phase of the rf pulse. * @param factor The probability of connection. */ diff --git a/src/c_lib/include/frequency_averaging.h b/src/c_lib/include/frequency_averaging.h index eb69238e2..5946e2cee 100644 --- a/src/c_lib/include/frequency_averaging.h +++ b/src/c_lib/include/frequency_averaging.h @@ -10,10 +10,8 @@ #include "simulation.h" void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - MRS_fftw_scheme *fftw_scheme, double *spec, - double transition_pathway_weight); + MRS_fftw_scheme *fftw_scheme, double *spec); void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, MRS_fftw_scheme *fftw_scheme, double *spec, - double transition_pathway_weight, unsigned int number_of_sidebands, double *affine_matrix); diff --git a/src/c_lib/include/simulation.h b/src/c_lib/include/simulation.h index 08ca84cae..54107d057 100644 --- a/src/c_lib/include/simulation.h +++ b/src/c_lib/include/simulation.h @@ -28,8 +28,8 @@ extern void mrsimulator_core( double rotor_frequency_in_Hz, // The rotor spin frequency. double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis. - float *transition_pathway, // Pointer to a list of transitions. - double transition_pathway_weight, // The weight of transition pathway. + float *transition_pathway, // Pointer to a list of transitions. + // powder orientation average int integration_density, // The number of triangle along the edge of octahedron. unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. @@ -45,13 +45,12 @@ extern void __mrsimulator_core( // transition is a list of quantum numbers packed as quantum numbers from the // initial energy state followed by the quantum numbers from the final energy state. // The energy states are given in Zeeman basis. - float *transition_pathway, // Pointer to a list of transitions. - double transition_pathway_weight, // The weight of transition pathway. - int n_dimension, // The total number of spectroscopic dimensions. - MRS_dimension *dimensions, // Pointer to MRS_dimension structure. - MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. - MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. - bool interpolation, // If true, perform a 1D interpolation. + float *transition_pathway, + int n_dimension, // The total number of spectroscopic dimensions. + MRS_dimension *dimensions, // Pointer to MRS_dimension structure. + MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. + MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. + bool interpolation, // If true, perform a 1D interpolation. /** * Each event consists of the following freq contrib ordered as diff --git a/src/c_lib/lib/angular_momentum/wigner_element.c b/src/c_lib/lib/angular_momentum/wigner_element.c index 15d96206f..5003c094d 100644 --- a/src/c_lib/lib/angular_momentum/wigner_element.c +++ b/src/c_lib/lib/angular_momentum/wigner_element.c @@ -282,7 +282,7 @@ void transition_connect_factor(const float l, const float m1_f, const float m1_i const double phi, double *restrict factor) { unsigned int delta_p = (unsigned int)((m2_f - m2_i) - (m1_f - m1_i)); double scale, phase = (double)delta_p * phi; - double cx = cos(phase), sx = sin(phase), re, im, temp; + double cx = cos(phase), sx = sin(phase); delta_p %= 4; scale = __wigner_d_element(l, m2_f, m1_f, theta); @@ -290,25 +290,22 @@ void transition_connect_factor(const float l, const float m1_f, const float m1_i switch (delta_p) { case 0: // * 1 - re = scale * cx; - im = scale * sx; + factor[0] = scale * cx; + factor[1] = scale * sx; break; case 1: // * -I - re = scale * sx; - im = -scale * cx; + factor[0] = scale * sx; + factor[1] = -scale * cx; break; case 2: // * -1 - re = -scale * cx; - im = -scale * sx; + factor[0] = -scale * cx; + factor[1] = -scale * sx; break; case 3: // * I - re = -scale * sx; - im = scale * cx; + factor[0] = -scale * sx; + factor[1] = scale * cx; break; } - temp = factor[0] * re - factor[1] * im; - factor[1] = factor[0] * im + factor[1] * re; - factor[0] = temp; } // void transition_connect_factor(const float *l, const float *transition_inital, diff --git a/src/c_lib/lib/frequency_averaging.c b/src/c_lib/lib/frequency_averaging.c index 8a59d9e9c..9088adc25 100644 --- a/src/c_lib/lib/frequency_averaging.c +++ b/src/c_lib/lib/frequency_averaging.c @@ -21,8 +21,7 @@ static inline void get_multi_event_amplitudes(int n_events, MRS_event *restrict } static inline void averaging_1D(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - MRS_fftw_scheme *fftw_scheme, double *spec, - double transition_pathway_weight) { + MRS_fftw_scheme *fftw_scheme, double *spec) { unsigned int i, j, k1, address; unsigned int nt = scheme->integration_density, npts = scheme->octant_orientations; @@ -56,9 +55,6 @@ static inline void averaging_1D(MRS_dimension *dimensions, MRS_averaging_scheme } } - cblas_dscal(dimensions->events->plan->size, transition_pathway_weight, - dimensions->events->freq_amplitude, 1); - offset_0 = dimensions->normalize_offset + dimensions->R0_offset; if (fabs(*freq - freq[nt]) < TOL && fabs(*freq - freq[npts - 1]) < TOL) delta_interpolation = true; @@ -99,19 +95,17 @@ static inline void averaging_1D(MRS_dimension *dimensions, MRS_averaging_scheme } void one_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, - MRS_fftw_scheme *fftw_scheme, double *spec, - double transition_pathway_weight) { + MRS_fftw_scheme *fftw_scheme, double *spec) { // multiply amplitudes from all events to the amplitude array from the first event. if (dimensions->n_events != 1) { get_multi_event_amplitudes(dimensions->n_events, dimensions->events, dimensions->events->plan->size); } - averaging_1D(dimensions, scheme, fftw_scheme, spec, transition_pathway_weight); + averaging_1D(dimensions, scheme, fftw_scheme, spec); } void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme *scheme, MRS_fftw_scheme *fftw_scheme, double *spec, - double transition_pathway_weight, unsigned int number_of_sidebands, double *affine_matrix) { unsigned int i, k, j, evt; @@ -177,7 +171,6 @@ void two_dimensional_averaging(MRS_dimension *dimensions, MRS_averaging_scheme * cblas_dscal(planA->n_octants * number_of_sidebands, planA->norm_amplitudes[j], &freq_ampB[j], scheme->octant_orientations); } - cblas_dscal(size, transition_pathway_weight, freq_ampA, 1); for (i = 0; i < number_of_sidebands; i++) { offsetA = offset0 + planA->vr_freq[i] * dimensions[0].inverse_increment; diff --git a/src/c_lib/lib/mrsimulator.c b/src/c_lib/lib/mrsimulator.c index 4b74e8914..e7d702097 100644 --- a/src/c_lib/lib/mrsimulator.c +++ b/src/c_lib/lib/mrsimulator.c @@ -330,7 +330,6 @@ void MRS_get_amplitudes_from_plan(MRS_averaging_scheme *scheme, MRS_plan *plan, * operation again updates the values of the array, `vector`. */ fftw_execute(fftw_scheme->the_fftw_plan); - /** ONLY VALID FOR SINGLE EVENT **/ /** * Evaluate the absolute value square of the `vector` array. The absolute value square * is stores as the real part of the `vector` array. The imaginary part is now @@ -338,6 +337,15 @@ void MRS_get_amplitudes_from_plan(MRS_averaging_scheme *scheme, MRS_plan *plan, vm_double_square_inplace(2 * plan->size, (double *)fftw_scheme->vector); cblas_daxpy(plan->size, 1.0, (double *)fftw_scheme->vector + 1, 2, (double *)fftw_scheme->vector, 2); + + /* Scaling the absolute value square with the powder scheme weights. Only the real + * part is scaled and the imaginary part is left as is. + */ + // for (i = 0; i < scheme->octant_orientations; i++) { + // cblas_dscal(plan->n_octants * plan->number_of_sidebands, + // plan->norm_amplitudes[i], (double *)&fftw_scheme->vector[i], + // 2 * scheme->octant_orientations); + // } } /** @@ -604,7 +612,7 @@ void MRS_rotate_components_from_PAS_to_common_frame( * = (2π / m ωr) (exp(I m ωr t) - 1) * |--scale--| * = scale (exp(I m ωr t) - 1) - * = scale [[cos(m ωr t) -1] + I sin(m ωr t)], + * = scale [[cos(m ωr t) -1] +Isin(m ωr t)], * * where ωr is the sample spinning frequency in Hz, m goes from -4 to 4, t is a vector * of length `number_of_sidebands` given as diff --git a/src/c_lib/lib/simulation.c b/src/c_lib/lib/simulation.c index b334b0ca3..6ca109cf0 100644 --- a/src/c_lib/lib/simulation.c +++ b/src/c_lib/lib/simulation.c @@ -44,13 +44,12 @@ void __mrsimulator_core( // transition is a list of quantum numbers packed as quantum numbers from the // initial energy state followed by the quantum numbers from the final energy state. // The energy states are given in Zeeman basis. - float *transition_pathway, // Pointer to the transition pathway, - double transition_pathway_weight, // The weight of transition pathway. - int n_dimension, // The total number of spectroscopic dimensions. - MRS_dimension *dimensions, // Pointer to MRS_dimension structure. - MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. - MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. - bool interpolation, // If true, perform a 1D interpolation. + float *transition_pathway, + int n_dimension, // The total number of spectroscopic dimensions. + MRS_dimension *dimensions, // Pointer to MRS_dimension structure. + MRS_fftw_scheme *fftw_scheme, // Pointer to the fftw scheme. + MRS_averaging_scheme *scheme, // Pointer to the powder averaging scheme. + bool interpolation, // If true, perform a 1D interpolation. /** * Each event consists of the following freq contrib ordered as @@ -163,13 +162,11 @@ void __mrsimulator_core( switch (n_dimension) { case 1: - one_dimensional_averaging(dimensions, scheme, fftw_scheme, spec, - transition_pathway_weight); + one_dimensional_averaging(dimensions, scheme, fftw_scheme, spec); break; case 2: two_dimensional_averaging(dimensions, scheme, fftw_scheme, spec, - transition_pathway_weight, plan->number_of_sidebands, - affine_matrix); + plan->number_of_sidebands, affine_matrix); break; } } @@ -192,7 +189,8 @@ void mrsimulator_core( double rotor_angle_in_rad, // The rotor angle relative to lab-frame z-axis // Pointer to the a list of transitions. - float *transition_pathway, double transition_pathway_weight, + float *transition_pathway, + // powder orientation average int integration_density, // The number of triangle along the edge of octahedron unsigned int integration_volume, // 0-octant, 1-hemisphere, 2-sphere. @@ -231,8 +229,9 @@ void mrsimulator_core( sites, // Pointer to a list of sites within the spin system. couplings, // Pointer to a list of couplings within a spin system. transition_pathway, // Pointer to a list of transition. - transition_pathway_weight, n_dimension, dimensions, fftw_scheme, scheme, - interpolation, freq_contrib, affine_matrix); + + n_dimension, dimensions, fftw_scheme, scheme, interpolation, freq_contrib, + affine_matrix); // gettimeofday(&end, NULL); // clock_time = (double)(end.tv_usec - begin.tv_usec) / 1000000. + diff --git a/src/mrsimulator/__init__.py b/src/mrsimulator/__init__.py index 538bb3978..b95b805a7 100644 --- a/src/mrsimulator/__init__.py +++ b/src/mrsimulator/__init__.py @@ -22,7 +22,7 @@ __license__ = "BSD License" __maintainer__ = "Deepansh J. Srivastava" __status__ = "Beta" -__version__ = "0.7.0dev1" +__version__ = "0.7.0dev0" import os diff --git a/src/mrsimulator/benchmark.py b/src/mrsimulator/benchmark.py index 8d226466c..011ecdd79 100644 --- a/src/mrsimulator/benchmark.py +++ b/src/mrsimulator/benchmark.py @@ -5,41 +5,41 @@ import numpy as np from mrsimulator import __version__ from mrsimulator import Simulator +from mrsimulator import Site +from mrsimulator import SpinSystem from mrsimulator.methods import BlochDecayCentralTransitionSpectrum from mrsimulator.methods import BlochDecaySpectrum -from mrsimulator.utils.collection import single_site_system_generator - -# from mrsimulator import SpinSystem -# from mrsimulator import Site # import platform # os_system = platform.system() -__author__ = ["Deepansh Srivastava", "Matthew D. Giammar"] -__email__ = ["srivastava.89@osu.edu", "giammar.7@buckeyemail.osu.edu"] - - def generate_spin_half_spin_system(n=1000): iso = np.random.normal(loc=0.0, scale=10.0, size=n) zeta = np.random.normal(loc=50.0, scale=15.12, size=n) eta = np.random.normal(loc=0.25, scale=0.012, size=n) - return single_site_system_generator( - isotope="29Si", - isotropic_chemical_shift=iso, - shielding_symmetric={"zeta": zeta, "eta": eta}, - ) + spin_systems = [] + for i, z, e in zip(iso, zeta, eta): + site = Site( + isotope="29Si", + isotropic_chemical_shift=i, + shielding_symmetric={"zeta": z, "eta": e}, + ) + spin_systems.append(SpinSystem(sites=[site])) + return spin_systems def generate_spin_half_int_quad_spin_system(n=1000): iso = np.random.normal(loc=0.0, scale=10.0, size=n) Cq = np.random.normal(loc=5.0e6, scale=1e5, size=n) eta = np.random.normal(loc=0.1, scale=0.012, size=n) - return single_site_system_generator( - isotope="17O", - isotropic_chemical_shift=iso, - quadrupolar={"Cq": Cq, "eta": eta}, - ) + spin_systems = [] + for i, c, e in zip(iso, Cq, eta): + site = Site( + isotope="17O", isotropic_chemical_shift=i, quadrupolar={"Cq": c, "eta": e} + ) + spin_systems.append(SpinSystem(sites=[site])) + return spin_systems def generate_simulator(spin_systems, method): diff --git a/src/mrsimulator/method/__init__.py b/src/mrsimulator/method/__init__.py index 673f4df3e..2684b91fa 100644 --- a/src/mrsimulator/method/__init__.py +++ b/src/mrsimulator/method/__init__.py @@ -6,10 +6,7 @@ from typing import Union import csdmpy as cp -import matplotlib as mpl import numpy as np -from mrsimulator.base_model import transition_connect_factor -import pandas as pd from mrsimulator.spin_system.isotope import Isotope from mrsimulator.transition import SymmetryPathway from mrsimulator.transition import Transition @@ -19,16 +16,12 @@ from pydantic import PrivateAttr from pydantic import validator -from .plot import plot as _plot from .spectral_dimension import CHANNELS from .spectral_dimension import SpectralDimension from .utils import cartesian_product -from .utils import mixing_query_connect_map -from .utils import tip_angle_and_phase_list - -__author__ = ["Deepansh J. Srivastava", "Matthew D. Giammar"] -__email__ = ["srivastava.89@osu.edu", "giammar.7@buckeyemail.osu.edu"] +__author__ = "Deepansh J. Srivastava" +__email__ = "srivastava.89@osu.edu" class Method(Parseable): @@ -278,17 +271,19 @@ def json(self, units=True) -> dict: Returns: dict """ - mth = {k: getattr(self, k) for k in ["name", "label", "description"]} + # mth = super().json(units=unit) + mth = {_: self.__getattribute__(_) for _ in ["name", "label", "description"]} mth["channels"] = [item.json() for item in self.channels] mth["spectral_dimensions"] = [ item.json(units=units) for item in self.spectral_dimensions ] # add global parameters + evt_d = self.property_units.items() global_ = ( - {k: f"{getattr(self, k)} {u}" for k, u in self.property_units.items()} + {k: f"{self.__getattribute__(k)} {u}" for k, u in evt_d} if units - else {k: getattr(self, k) for k in self.property_units} + else {k: self.__getattribute__(k) for k, u in evt_d} ) mth.update(global_) @@ -404,7 +399,7 @@ def get_symmetry_pathways(self, symmetry_element: str) -> List[SymmetryPathway]: dim._get_symmetry_pathways(symmetry_element) for dim in self.spectral_dimensions ] - sp_idx = np.arange(len(sym_path)) + sp_indexes = np.arange(len(sym_path)) indexes = [np.arange(len(item)) for item in sym_path] products = cartesian_product(*indexes) @@ -412,7 +407,11 @@ def get_symmetry_pathways(self, symmetry_element: str) -> List[SymmetryPathway]: SymmetryPathway( channels=self.channels, **{ - ch: [_ for sp, i in zip(sp_idx, item) for _ in sym_path[sp][i][ch]] + ch: [ + _ + for sp, i in zip(sp_indexes, item) + for _ in sym_path[sp][i][ch] + ] for ch in CHANNELS }, ) @@ -431,74 +430,17 @@ def _get_transition_pathways_np(self, spin_system): evt.filter_transitions(all_transitions, isotopes, channels) for dim in self.spectral_dimensions for evt in dim.events - if evt.__class__.__name__ != "MixingEvent" ] + # if segments == []: + # return [] + segments_index = [np.arange(item.shape[0]) for item in segments] cartesian_index = cartesian_product(*segments_index) return [ [segments[i][j] for i, j in enumerate(item)] for item in cartesian_index ] - def _get_transition_pathway_weights(self, pathways, spin_system): - if pathways == []: - return [] - - symbol = [item.isotope.symbol for item in spin_system.sites] - spins = np.asarray([item.isotope.spin for item in spin_system.sites]) - channels = [item.symbol for item in self.channels] - weights = np.ones(len(pathways), dtype=complex) - - mapping = mixing_query_connect_map(self.spectral_dimensions) - for obj in mapping: - theta_, phi_ = tip_angle_and_phase_list( - symbol, channels, obj["mixing_query"] - ) - map_ = obj["near_index"] - for j, path in enumerate(pathways): - # if weights[j] == 0: - # continue - weights[j] *= self._calculate_transition_connect_weight( - path[map_[0]], path[map_[1]], spins, theta_, phi_ - ) - return np.round(weights, decimals=6) - - @staticmethod - def _calculate_transition_connect_weight(trans1, trans2, spins, theta, phi): - """Return the transition connection weight of transitions `trans1` and `trans2`. - - Args: - trans1: ndarray of shape (2, n_sites) representing the starting transition, - where `n_sites` is the number of sites within the spin system. The two - entries at axis 0 corresponds to initial (index 0) and final (index 1) - Zeeman energy states. - trans2: ndarray of shape (2, n_sites) representing the target transition, - where `n_sites` is the number of sites within the spin system. The two - entries at axis 0 corresponds to initial (index 0) and final (index 1) - Zeeman energy states. - spins: ndarray of spin quantum numbers of the sites. - theta: ndarray of rotation tip_angle per site from the mixing event. - phi: ndarray of rotation phase per site from the mixing event. - """ - amp = 1 - for i, spin in enumerate(spins): - m1_f = trans1[1][i] # starting transition final state - m1_i = trans1[0][i] # starting transition initial state - m2_f = trans2[1][i] # landing transition final state - m2_i = trans2[0][i] # landing transition initial state - - amp *= transition_connect_factor( - spin, m1_f, m1_i, m2_f, m2_i, theta[i], phi[i] - ) - return amp - - def _get_transition_pathway_and_weights_np(self, spin_system): - segments = self._get_transition_pathways_np(spin_system) - weights = self._get_transition_pathway_weights(segments, spin_system) - return segments, weights - # indexes = np.where(weights != 0)[0] - # return np.asarray(segments)[indexes], weights[indexes] - def get_transition_pathways(self, spin_system) -> List[TransitionPathway]: """Return a list of transition pathways from the given spin system that satisfy the query selection criterion of the method. @@ -516,255 +458,21 @@ def get_transition_pathways(self, spin_system) -> List[TransitionPathway]: >>> sys = SpinSystem(sites=[{'isotope': '27Al'}, {'isotope': '29Si'}]) >>> method = ThreeQ_VAS(channels=['27Al']) >>> pprint(method.get_transition_pathways(sys)) - [|1.5, -0.5⟩⟨-1.5, -0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, weight=(1+0j), - |1.5, -0.5⟩⟨-1.5, -0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(1+0j), - |1.5, 0.5⟩⟨-1.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, weight=(1+0j), - |1.5, 0.5⟩⟨-1.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(1+0j)] + [|1.5, -0.5⟩⟨-1.5, -0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, + |1.5, -0.5⟩⟨-1.5, -0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, + |1.5, 0.5⟩⟨-1.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, + |1.5, 0.5⟩⟨-1.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|] """ - segments, weights = self._get_transition_pathway_and_weights_np(spin_system) + segments = self._get_transition_pathways_np(spin_system) return [ TransitionPathway( [ Transition(initial=tr[0].tolist(), final=tr[1].tolist()) for tr in item - ], - weight=w, - ) - for item, w in zip(segments, weights) - if w != 0 - ] - - def _add_simple_props_to_df(self, df, prop_dict, required, drop_constant_columns): - """Helper method for summary to reduce complexity""" - # Iterate through property and valid Event subclass for property - for prop, valid in prop_dict.items(): - lst = [ - getattr(ev, prop) if ev.__class__.__name__ in valid else np.nan - for dim in self.spectral_dimensions - for ev in dim.events - ] - if prop not in required and drop_constant_columns: - # NOTE: np.isnan() cannot be passed an object (freq_contrib) - lst_copy = np.asarray(lst)[~np.isnan(np.asarray(lst))] - if np.unique(lst_copy).size < 2: - continue - df[prop] = lst - - def summary(self, drop_constant_columns=True) -> pd.DataFrame: - r"""Returns a DataFrame giving a summary of the Method. A user can specify - optional attributes to include which appear as columns in the DataFrame. A user - can also ask to leave out attributes which remain constant throughout the - method. Invalid attributes for an Event will be replaced with NAN. - - Args: - (bool) drop_constant_columns: - Removes constant properties if True. Default is True. - - Returns: - pd.DataFrame df: - Event number as row and property as column. Invalid properties for an - event type are filled with np.nan - - **Columns** - - - (str) type: Event type - - (int) spec_dim_index: Index of spectral dimension which event belongs to - - (str) label: Event label - - (float) duration: Duration of the ConstantDurationEvent - - (float) fraction: Fraction of the SpectralEvent - - (MixingQuery) mixing_query: MixingQuery object of the MixingEvent - - (float) magnetic_flux_density: Magnetic flux density during event in Tesla - - (float) rotor_frequency: Rotor frequency during event in Hz - - (float) rotor_angle: Rotor angle during event converted to Degrees - - (FrequencyEnum) freq_contrib: - - Example: - **User Defined Method2D Example** - - >>> from mrsimulator.methods import Method2D - >>> method = Method2D( - ... channels=['1H'], - ... spectral_dimensions=[ - ... { - ... "events": [ - ... { - ... "fraction": 0.7, - ... "transition_query": [{"ch1": {"P": [1]}}] - ... }, - ... { - ... "duration": 1.8, - ... "transition_query": [{"ch1": {"P": [0], "D": [0]}}] - ... } - ... ], - ... }, - ... { - ... "events": [ - ... { - ... "fraction": 0.3, - ... "transition_query": [{"ch1": {"P": [-1]}}] - ... }, - ... ], - ... } - ... ] - ... ) - >>> df = method.summary() - >>> # Columns are reduced to fit within 80 lines - >>> drop = ["freq_contrib", "spec_dim_label", "label", "mixing_query"] - >>> df.drop(drop, axis=1, inplace=True) - >>> pprint(df) - type spec_dim_index duration fraction p d - 0 SpectralEvent 0 NaN 0.7 [1.0] [nan] - 1 ConstantDurationEvent 0 1.8 NaN [0.0] [0.0] - 2 SpectralEvent 1 NaN 0.3 [-1.0] [nan] - - **All Possible Columns** - - >>> from mrsimulator.methods import ThreeQ_VAS - >>> method = ThreeQ_VAS(channels=["17O"]) - >>> df = method.summary(drop_constant_columns=False) - >>> pprint(list(df.columns)) - ['type', - 'spec_dim_index', - 'spec_dim_label', - 'label', - 'duration', - 'fraction', - 'mixing_query', - 'magnetic_flux_density', - 'rotor_frequency', - 'rotor_angle', - 'freq_contrib', - 'p', - 'd'] - """ - # Make sure spectral_dimensions has at least one SpectralDimension - if len(self.spectral_dimensions) == 0: - raise AttributeError( - "Method has empty spectral_dimensions. At least one SpectralDimension " - "is needed with at least one Event." + ] ) - - # Make sure there is at least one event within spectral_dimensions - if sum([len(spec_dim.events) for spec_dim in self.spectral_dimensions]) == 0: - raise AttributeError( - "Method has no Events. At least one SpectralDimension " - "is needed with at least one Event." - ) - - CD = "ConstantDurationEvent" - SP = "SpectralEvent" - MX = "MixingEvent" - - # Columns required to be present - required = [ - "type", - "label", - "duration", - "fraction", - "mixing_query", - "spec_dim_index", - "spec_dim_label", - "freq_contrib", - "p", - "d", - ] - - # Properties which can accessed by getattr() - prop_dict = { - "label": (CD, SP, MX), - "duration": CD, - "fraction": SP, - "mixing_query": MX, - "magnetic_flux_density": (CD, SP), - "rotor_frequency": (CD, SP), - "rotor_angle": (CD, SP), - "freq_contrib": (CD, SP), - } - - # Create the DataFrame - df = pd.DataFrame() - - # Populate columns which cannot be calculated from iteration - df["type"] = [ - ev.__class__.__name__ - for dim in self.spectral_dimensions - for ev in dim.events + for item in segments ] - df["spec_dim_index"] = [ - i for i, dim in enumerate(self.spectral_dimensions) for ev in dim.events - ] - df["spec_dim_label"] = [ - dim.label for dim in self.spectral_dimensions for ev in dim.events - ] - self._add_simple_props_to_df(df, prop_dict, required, drop_constant_columns) - - # Add p and d symmetry pathways to dataframe - df["p"] = np.transpose( - [sym.total for sym in self.get_symmetry_pathways("P")] - ).tolist() - df["d"] = np.transpose( - [sym.total for sym in self.get_symmetry_pathways("D")] - ).tolist() - - # Convert rotor_angle to degrees - if "rotor_angle" in df.columns: - df["rotor_angle"] = df["rotor_angle"] * 180 / np.pi - - # Convert rotor_frequency to kHz - if "rotor_frequency" in df.columns: - df["rotor_frequency"] = df["rotor_frequency"] / 1000 - - return df - - def plot(self, df=None, include_legend=False) -> mpl.pyplot.figure: - """Creates a diagram representing the method. By default, only parameters which - vary throughout the method are plotted. Figure can be finley adjusted using - matplotlib rcParams. - - Args: - DataFrame df: - DataFrame to plot data from. By default DataFrame is calculated from - summary() and will show only parameters which vary throughout the - method plus 'p' symmetry pathway and 'd' symmetry pathway if it is not - none or defined - - bool include_legend: - Optional argument to include a key for event colors. Default is False - and no key will be included in figure - - Returns: - matplotlib.pyplot.figure - - Example: - >>> from mrsimulator.methods import BlochDecaySpectrum - >>> method = BlochDecaySpectrum(channels=["13C"]) - >>> fig = method.plot() - - **Adjusting Figure Size rcParams** - - >>> import matplotlib as mpl - >>> from mrsimulator.methods import FiveQ_VAS - >>> mpl.rcParams["figure.figsize"] = [14, 10] - >>> mpl.rcParams["font.size"] = 14 - >>> method = FiveQ_VAS(channels=["27Al"]) - >>> fig = method.plot(include_legend=True) - - **Plotting all Parameters, including Constant** - - >>> from mrsimulator.methods import FiveQ_VAS - >>> method = FiveQ_VAS(channels=["27Al"]) - >>> df = method.summary(drop_constant_columns=False) - >>> fig = method.plot(df=df) - """ - fig = mpl.pyplot.gcf() - - if df is None: - df = self.summary() - - _plot(fig=fig, df=df, include_legend=include_legend) - - fig.suptitle(t=self.name if self.name is not None else "") - return fig def shape(self) -> tuple: """The shape of the method's spectral dimension array. diff --git a/src/mrsimulator/method/event.py b/src/mrsimulator/method/event.py index a41b578c6..9f1f57e57 100644 --- a/src/mrsimulator/method/event.py +++ b/src/mrsimulator/method/event.py @@ -22,24 +22,23 @@ class BaseEvent(Parseable): - """Base BaseEvent class. If the value of the attribute is None, the value of the - corresponding global attribute will be used instead. + """Base BaseEvent class. Attributes ---------- magnetic_flux_density: The macroscopic magnetic flux density, :math:`H_0`, of the applied external - magnetic field during the event in units of T. The default value is ``None``. + magnetic field during the event in units of T. The default value is ``9.4``. rotor_frequency: The sample spinning frequency :math:`\nu_r`, during the event in units of Hz. - The default value is ``None``. + The default value is ``0``. rotor_angle: The angle between the sample rotation axis and the applied external magnetic field vector, :math:`\theta`, during the event in units of rad. - The default value is ``None``, i.e. the magic angle. + The default value is ``0.9553166``, i.e. the magic angle. freq_contrib: A list of FrequencyEnum enumeration. The default is all frequency enumerations. @@ -112,11 +111,10 @@ def filter_transitions(self, all_transitions, isotopes, channels): """ symmetry_permutations = self.permutation(isotopes, channels) - # print("symmetry_permutations", symmetry_permutations) + segment = [] for item in symmetry_permutations: st = all_transitions[:] - # print(item["P"].size, item["D"].size) st = st[P_symmetry_indexes(st, item["P"])] if item["P"].size > 0 else st st = st[D_symmetry_indexes(st, item["D"])] if item["D"].size > 0 else st segment += [st] diff --git a/src/mrsimulator/method/frequency_contrib.py b/src/mrsimulator/method/frequency_contrib.py index 75d8f4959..cbfc8d12c 100644 --- a/src/mrsimulator/method/frequency_contrib.py +++ b/src/mrsimulator/method/frequency_contrib.py @@ -73,7 +73,7 @@ class FrequencyEnum(str, Enum): # values[item] = item # return None - def json(self, **kwargs) -> str: + def json(self, **kwargs) -> dict: """Parse the class object to a JSON compliant python dictionary object.""" return self.value diff --git a/src/mrsimulator/method/plot.py b/src/mrsimulator/method/plot.py deleted file mode 100644 index 5291c88a6..000000000 --- a/src/mrsimulator/method/plot.py +++ /dev/null @@ -1,728 +0,0 @@ -# -*- coding: utf-8 -*- -from itertools import groupby -from math import ceil - -import matplotlib.projections as proj -import matplotlib.pyplot as plt -import numpy as np -from matplotlib.patches import Patch -from matplotlib.patches import Rectangle -from matplotlib.ticker import FixedLocator -from matplotlib.ticker import MaxNLocator - -__author__ = "Matthew D. Giammar" -__email__ = "giammar.7@buckeyemail.osu.edu" - - -# NOTE: Matplotlib should automatically generate new colors when none specified -DURATION_WIDTH = 0.5 # Width of one ConstantDurationEvent -SPECTRAL_MULTIPLIER = 0.8 # Width multiplier for all SpectralEvents -MIXING_WIDTH = 0.25 # tip_angle of 360 degrees -DEFAULTFONTSIZE = 9 -MAX_SYMMETRY_TICKS = 5 # Maximum number of ticks allowed on symmetry plot, always odd -COLORS = { - "ConstantDurationEvent": "orange", - "SpectralEvent": "g", - "MixingEvent": "b", - "inf_speed": "y", - "undef": "k", -} -LABLES = { - "rotor_angle": r"$\theta_r$ / °", - "rotor_frequency": r"$\nu_r$ / kHz", - "magnetic_flux_density": r"$B_0$ / T", -} -DEFAULT_ANNO_KWARGS = { - "annotation_clip": False, - "color": "#202020", - "ha": "center", - "va": "center", - "fontsize": 10, -} - - -class CustomAxes(plt.Axes): - """Subclass of matplotlib Axes holding specific functions and data - - Attributes - ---------- - - name: - String of the matplotlib projection name of this Axes subclass. This value - is not intended to be changed - - x_data: - Unformatted x-axis data to be plotted - - y_data: - Unformatted y-axis data to be plotted - - col_name: - Column which the y_data came from - - xmax: - Maximum x value in x_data. Used to set xlim to (0, xmax) - - mix_ev: - True where event is MixingEvent, false otherwise - - """ - - name = "custom_axes" - x_data = None - y_data = None - col_name = None - xmax = None - mix_ev = None - - def make_plot(self, x_data, y_data, col_name, mix_ev, format_kwargs, plot_kwargs): - """Main workflow function to format and plot data on Axes""" - self.x_data = np.array(x_data, dtype=float) - self.y_data = np.array(y_data, dtype=float) - self.col_name = col_name - self.xmax = max(x_data) - self.mix_ev = np.array(mix_ev) - - self._format(**format_kwargs) - self.plot(x=self.x_data, y=self.y_data, **plot_kwargs) - - def plot(self, x, y, _format=True, labels=True, **kwargs): - """Plot formatted data""" - if _format: - # matplotlib automatically blanks out nan values. y may contain nans - # mask only differs from y when plotting "rotor_frequency" - x, y, y_mask = self._format_x_and_y_data(x=x, y=y) - else: - y_mask = y - - # 'x' and 'y' cannot be passed as keyword args to plt.Axes.plot() - super().plot(x, y_mask, **kwargs) - - # labels must come after plot to setup y-axis limits - if labels: - self._add_blank_space_labels(x=x, y=y) - - def _add_blank_space_labels(self, x, y): - """Adds labels for blank spaces in plot""" - if self.col_name == "rotor_frequency": - # Locate places with infinite spinning speed (1e9 Hz) - regions = self._locate_regions_with_val(x=x, y=y, val=1e9) - for region in regions: - self._add_rect_with_label( - x0=region[0], - x1=region[1], - label="∞ speed", - rect_kwargs=dict(color=COLORS["inf_speed"]), - ) - # Locate places with undefined parameters - regions = self._locate_regions_with_val(x=x, y=y) - for region in regions: - self._add_rect_with_label( - x0=region[0], - x1=region[1], - label="undefined", - rect_kwargs=dict(color=COLORS["undef"]), - ) - - def _locate_regions_with_val(self, x, y, val=np.nan): - """Locates regions and returns a list of tuples denoting range with value""" - # Locate value in y - if val is np.nan: - loc = np.isnan(y) - else: - loc = y >= val - - # Find indexes where loc_k and loc_k+1 are both true - left_x = np.argwhere(np.logical_and(loc[1:], loc[:-1])).flatten() - - # Create tuples of left and right x points - return [(x[i], x[i + 1]) for i in left_x] - - def _format_x_and_y_data(self, x, y): - """Removes invalid event data (nans from MixingEvents) while keeping valid but - undefined event data in y. Extends y double to achieve stair-step pattern - """ - y_data = y[~self.mix_ev] # Keep data from events that are not MixingEvents - x_data = x - y_data = np.asarray([n for num in y_data for n in (num, num)]) # Extend double - - # Mask pseudo-infinite spinning speed with nans for plotting "rotor_frequency" - if self.col_name == "rotor_frequency": - mask = np.where(y_data >= 1e9, np.nan, y_data) - return x_data, y_data, mask - - # Insert zero or remove first x point for p and d - if self.col_name in ["p", "d"]: - if self.mix_ev[0]: - y_data = np.insert(y_data, 0, 0) - else: - x_data = x[1:] - - return x_data, y_data, y_data - - def _format( - self, - locator=None, - grid={ - "axis": "y", - "color": "black", - "linestyle": "--", - "linewidth": 0.5, - "alpha": 0.5, - }, - ymargin=0.1, - show_xaxis=False, - show_yaxis=True, - **kwargs, - ): - """Format Axes helper function - - Args: - (matplotlib.ticker) locator: - Object which determines the ticks on the y-axis - - (dict) grid: - Grid kwargs which determine the appearance of the grid - - (float) ymargin: - Value to be passed to matplotlib.axes.Axes.margin for padding - - (bool) show_xaxis: - Shows x-axis and ticks if True, otherwise hidden - - (bool) show_yaxis: - Shows y-axis and ticks if True, otherwise hidden - """ - label = LABLES[self.col_name] if self.col_name in LABLES else self.col_name - - if locator is None: - locator = MaxNLocator(nbins=4, integer=True, min_n_ticks=3) - - if self.col_name == "rotor_frequency": - locator.set_params(integer=False, steps=[1.5, 5, 7]) - - self.get_xaxis().set_visible(show_xaxis) - self.set_xlim(0, self.xmax) - - self.get_yaxis().set_visible(show_yaxis) - self.get_yaxis().set_major_locator(locator) - self.set_ylabel(label, fontsize=DEFAULTFONTSIZE) - - self.spines["top"].set_visible(show_xaxis) - self.spines["right"].set_visible(show_xaxis) - self.spines["bottom"].set_visible(show_xaxis) - - self.grid(**grid) - self.margins(y=ymargin) - self.tick_params(axis="both", labelsize=DEFAULTFONTSIZE) - - def _add_rect_with_label(self, x0, x1, label, rect_kwargs, anno_kwargs={}): - """Add a rectangle between x0 and x1 on ax representing event""" - # [height, color, alpha] required in rect_kwargs - bottom, top = self.get_ylim() - if "height" not in rect_kwargs: - rect_kwargs["height"] = top - bottom - - if "alpha" not in rect_kwargs: - rect_kwargs["alpha"] = 0.2 - - if "color" not in rect_kwargs: - raise ValueError("No color in `rect_kwargs`. A color must be specified") - - if anno_kwargs == {}: - anno_kwargs = DEFAULT_ANNO_KWARGS - - rect = Rectangle(xy=(x0, bottom), width=x1 - x0, **rect_kwargs) - self.add_patch(rect) - if label is not None: - y_mid = (top + bottom) / 2 - self.annotate(text=label, xy=((x1 + x0) / 2, y_mid), **anno_kwargs) - - -class MultiLineAxes(CustomAxes): - """Axes subclass for multiline plots such as 'p' or 'd'""" - - name = "multi_line_axes" - - def make_plot(self, x_data, y_data, col_name, mix_ev, format_kwargs, plot_kwargs): - """Main workflow function to format and plot data on Axes""" - self.x_data = x_data - self.y_data = np.stack(y_data.values).transpose() # each row is symm pathway - - if np.asarray(self.y_data).ndim != 2: - raise ValueError("Symmetry pathway data is misshapen. Data must be 2d") - - self.col_name = col_name - self.xmax = max(x_data) - self.mix_ev = mix_ev - - self._format(**format_kwargs) - self.plot(x=self.x_data, y=self.y_data, **plot_kwargs) - - def plot(self, x, y, do_offset=True, **kwargs): - if do_offset: - y = self._offset_overlaps(y) - - for row in y: - super().plot(x=x, y=row, labels=False, **kwargs) - - def _format(self, ymargin=0.15, **kwargs): - """Determines locator for axes and calls super()._format()""" - ticks = self._make_ticks() - margin = ymargin + 1 - - # After set_ylim() is called, Axes.margins() had no effect - self.set_ylim(ticks[0] * margin, ticks[-1] * margin) - - super()._format(locator=FixedLocator(ticks), **kwargs) - - def _make_ticks(self): - """Logic for deciding where to place the ticks for a symmetry query""" - # Find maximum distance from origin in dataset - dist = max( - 1, np.nanmax(self.y_data.flatten()), -np.nanmin(self.y_data.flatten()) - ) - - # Calculate tick step so at most MAX_SYMMETRY_TICKS ticks - step = ceil(2.3 * dist / MAX_SYMMETRY_TICKS) - - # Round up d pathway to next even int if odd - if self.col_name == "d" and step % 2 == 1: - step += 1 - - # Round dist up to next greatest multiple of step to force symmetry - dist = (dist + step - 1) - ((dist + step - 1) % step) - - return np.arange(-dist, dist + 1, step, dtype=int) - - def _offset_overlaps(self, y, offset_pct=0.03): - """Offsets y at overlapping values""" - # calculate raw offset based on percent of data range [max(y) - min(y)] - offset = np.nanmax(y.flatten()) - np.nanmin(y.flatten()) - if offset == 0: - offset = 2 - offset *= offset_pct - - # Check for only one symmetry pathway present - if y.shape[0] == 1: - return y - data = y.transpose() # Each row represents all symmetry pathways for an event - - # Loop over rows and find where pathways overlap - for row in data: - unique = np.unique(row) - - # No symmetry pathway overlap - if unique.size == row.size: - continue - - # Find indices of overlap for each number and add offsets - for num in unique: - where = np.where(row == num)[0] - np.add.at( - row, - where, - np.linspace( - -where.size * offset / 2, - where.size * offset / 2, - num=where.size, - ), - ) - - return data.transpose() - - -class SequenceDiagram(CustomAxes): - """Axes subclass holding the sequence diagram of events""" - - name = "sequence_axes" - ylim = None - - def make_plot(self, df, x_data, ylim=[0, 1]): - """Main workflow function to plot sequence diagram""" - self.x_data = x_data - self.ylim = ylim - self.xmax = max(x_data) - - self._format(grid={}, ylim=self.ylim, margin=None, axis="off") - self._plot_sequence_diagram(df=df) - - def _format(self, **kwargs): - # Turn of x and y axis and call super(),_format - self.axis("off") - self.set_ylim(*self.ylim) - - super()._format(**kwargs) - - def _format_mix_label(self, ta, p): - """Helper method to format label for mixing events. Returns (str, float)""" - # Format to one decimal place if float, otherwise no decimal places - tip_angle = "{1:0.{0}f}".format(int(not float(ta).is_integer()), ta) - phase = "{1:0.{0}f}".format(int(not float(p).is_integer()), p) - return (f"({tip_angle}, {phase})", ta / 360 * MIXING_WIDTH) - - def _plot_spec_dim_lines(self, df, ev_groups): - """Adds lines and labels denoting spectral dimensions""" - x_data = self.x_data - ylim = self.ylim - plot_kwargs = dict(_format=False, labels=False) - - self.plot(x=[0, 0], y=ylim, color="black", **plot_kwargs) - - spec_dim_idx = 0 - df_idx = 0 - x_idx = 0 - for _type, num in ev_groups: - for j in range(num): - if x_idx < len(x_data) - 1 and x_data[x_idx] == x_data[x_idx + 1]: - x_idx += 1 - if df["spec_dim_index"][df_idx + j] != spec_dim_idx: # Next spec dim - x0 = x_data[x_idx] - if df["type"][df_idx + j] == "MixingEvent": - x0 += sum(df["tip_angle"][df_idx:j]) / 360.0 * MIXING_WIDTH - self.plot(x=[x0, x0], y=ylim, color="black", **plot_kwargs) - spec_dim_idx += 1 - x_idx += 1 - df_idx += num - - # Plot last spectral dimension - self.plot(x=[max(x_data)] * 2, y=ylim, color="black", **plot_kwargs) - - def _plot_sequence_diagram(self, df): - """Plots sequence diagram (order of events)""" - ev_groups = [(_type, sum(1 for _ in gp)) for _type, gp in groupby(df["type"])] - - # Last event is MixingEvent - if ev_groups[-1][0] == "MixingEvent": - self.x_data = np.append(self.x_data, self.x_data[-1]) - n_end_mix = ev_groups[-1][1] - # Total angle / 360 * MIXING_WIDTH - # Get last 'n_end_mix' events - offset = sum(df["tip_angle"][-n_end_mix:]) / 360 * MIXING_WIDTH - self.x_data[-2] -= offset - - self._plot_spec_dim_lines(df=df, ev_groups=ev_groups) - - df_idx = 0 - x_idx = 0 - for _type, num in ev_groups: - if _type == "MixingEvent": - # Leftmost x point of next MixingEvent rectangle - x0 = self.x_data[x_idx] - # Iterate over each MixingEvent in group and plot rectangle - for j in range(num): - label, width = self._format_mix_label( - ta=df["tip_angle"][df_idx + j], - p=df["phase"][df_idx + j], - ) - super()._add_rect_with_label( - x0=x0, - x1=x0 + width, - label=label, - rect_kwargs=dict( - color=COLORS["MixingEvent"], - alpha=0.2, - ), - anno_kwargs=dict( - color="black", - ha="center", - va="center", - rotation=90, - ), - ) - x0 += width - x_idx += 1 - else: - for j in range(num): - # Increment x_idx if no mixing event between events - if self.x_data[x_idx] == self.x_data[x_idx + 1]: - x_idx += 1 - super()._add_rect_with_label( - x0=self.x_data[x_idx], - x1=self.x_data[x_idx + 1], - label=df["label"][df_idx + j], - rect_kwargs=dict( - color=COLORS[df["type"][df_idx]], - alpha=0.2, - ), - ) - x_idx += 1 - - df_idx += num - - -class SpectralDimensionLabels(CustomAxes): - """Axes subclass with""" - - name = "label_axes" - - def plot_labels(self, df, x_data, ylim=[0, 1], style="italic", format_kwargs={}): - """Plot labels of spectral dimensions""" - self.xmax = max(x_data) - self._format(ylim, **format_kwargs) - - ev_groups = [(_type, sum(1 for _ in gp)) for _type, gp in groupby(df["type"])] - - last_spec_dim_x = 0 - spec_dim_idx = 0 - df_idx = 0 - x_idx = 0 - for _type, num in ev_groups: - for j in range(num): - if x_idx < len(x_data) - 1 and x_data[x_idx] == x_data[x_idx + 1]: - x_idx += 1 - if df["spec_dim_index"][df_idx + j] != spec_dim_idx: # Next spec dim - x0 = x_data[x_idx] - if df["type"][df_idx + j] == "MixingEvent": - x0 += sum(df["tip_angle"][df_idx:j]) / 360.0 * MIXING_WIDTH - self.annotate( - text=df["spec_dim_label"][df_idx + j - 1], - xy=((last_spec_dim_x + x0) / 2, ylim[0] + 0.17), - **DEFAULT_ANNO_KWARGS, - style=style, - ) - last_spec_dim_x = x0 - spec_dim_idx += 1 - x_idx += 1 - df_idx += num - - self.annotate( - text=df["spec_dim_label"][df_idx - 1], - xy=((last_spec_dim_x + max(x_data)) / 2, ylim[0] + 0.17), - **DEFAULT_ANNO_KWARGS, - style=style, - ) - - def _format(self, ylim, **kwargs): - self.axis("off") - self.set_ylim(*ylim) - super()._format(**kwargs) - - -def _check_columns(df): - """Helper method to ensure required columns are present. Returns list of included - optional parameters in dataframe columns. - """ - required = [ - "type", - "label", - "spec_dim_index", - "spec_dim_label", - "duration", - "fraction", - "mixing_query", - "p", - ] - - # Required columns are not present in df - if not (set(required)).issubset(set(df.columns)): - raise ValueError("Some required columns were not present in the DataFrame.") - - # Remove 'freq_contrib' since it cannot be plotted - if "freq_contrib" in df.columns: - required += ["freq_contrib"] - - return list(df.columns.drop(required)) - - -def _add_tip_angle_and_phase(df): - """Add tip_angle and phase columns to dataframe from mixing_query""" - # NOTE Only columns for ch1 are created (1 ch only for now) - # NOTE What should empty MixingQueries add? (MixingQuery default?) - # BUG: Empty mixing queries (i.e. "mixing_query": {}) gives query of `None` causing - # error to be thrown - df["tip_angle"] = [ - query.ch1.tip_angle * 180 / np.pi - if query.__class__.__name__ == "MixingQuery" - else np.nan - for query in df["mixing_query"] - ] - df["phase"] = [ - query.ch1.phase * 180 / np.pi - if query.__class__.__name__ == "MixingQuery" - else np.nan - for query in df["mixing_query"] - ] - - -def _format_df(df): - """Formats dataframe and returns params to plot""" - params = _check_columns(df) - - # Check if d symmetry pathway is all nan, removing if true - if np.all(np.isnan(df["d"].tolist())): - params.remove("d") - - _add_tip_angle_and_phase(df) - return params - - -def _make_x_data(df): - """Returns list of x points to use in plotting""" - points = [0] - - for i, row in df.iterrows(): - if row["type"] == "SpectralEvent": - next_x = points[-1] + (row["fraction"] * SPECTRAL_MULTIPLIER) - points.extend((next_x, next_x)) - elif row["type"] == "ConstantDurationEvent": - next_x = points[-1] + DURATION_WIDTH - points.extend((next_x, next_x)) - - points.pop() - return points - - -def _offset_x_data(df, x_data): - """Offsets x_data based on MixingEvents""" - offset_x = np.array([0] + x_data) - ev_groups = [(_type, sum(1 for _ in group)) for _type, group in groupby(df["type"])] - - # Remove MixingEvents from end of sequence - if ev_groups[-1][0] == "MixingEvent": - ev_groups.pop() - - df_idx = 0 - # Extend first jump if first event(s) are MixingEvent - if ev_groups[0][0] == "MixingEvent": - gp__ = ev_groups[0][1] - offset = sum(df["tip_angle"][0:gp__]) / 360.0 * MIXING_WIDTH - offset_x[1] += offset - # Increment event indexer by number of MixingEvents in first group - df_idx += gp__ - ev_groups.pop(0) - - x_idx = 1 - for _type, num in ev_groups: - if _type == "MixingEvent": - up_lim__ = df_idx + num - offset = sum(df["tip_angle"][df_idx:up_lim__]) / 360.0 * MIXING_WIDTH - offset_x[x_idx] -= offset / 2 - offset_x[x_idx + 1] += offset / 2 - x_idx += 1 - else: - if offset_x[x_idx] == offset_x[x_idx + 1]: - x_idx += 1 - x_idx += (num * 2) - 1 - # Increment event indexer by number of Events in group - df_idx += num - - return offset_x - - -def _make_normal_and_offset_x_data(df): - """Calculates proper x_data and returns normal and offset x_data""" - x_data = _make_x_data(df) - - if len(x_data) == 0: - raise ValueError( - "The DataFrame does not contain any SpectralEvents or " - "ConstandDurationEvents. At least one must be present to construct a plot" - ) - - return x_data, _offset_x_data(df, x_data) - - -def _add_legend(fig): - """Adds legend for event color and type""" - fig.legend( - handles=[ - Patch(facecolor=COLORS["MixingEvent"], alpha=0.2, label="MixingEvent"), - Patch(facecolor=COLORS["SpectralEvent"], alpha=0.2, label="SpectralEvent"), - Patch( - facecolor=COLORS["ConstantDurationEvent"], - alpha=0.2, - label="ConstantDurationEvent", - ), - ], - loc="upper left", - fontsize="small", - ) - - -def _calculate_n_channels(df): - """Calculates the number of channels present in the method DataFrame""" - # TODO: (future) implement functionality for calculation - # Currently hardcoded to 1 - # Maybe move this into _format_df? - return 1 - - -def plot(fig, df, include_legend) -> plt.figure: - """Plot symmetry pathways and other requested parameters on figure""" - # TODO: (future) add functionality for multiple channels - mix_ev = np.array(df["type"] == "MixingEvent") - params = _format_df(df) - x_data, x_offset = _make_normal_and_offset_x_data(df) - - # Calculations and setup gridspec object for number of subplots - n_channels = _calculate_n_channels(df) - # spec_dim label + channels + p pathway + params - nrows = n_channels + 2 + len(params) - height_ratios = np.ones(nrows) - height_ratios[0] = 0.11 - height_ratios[1] = 0.75 - gs = fig.add_gridspec( - nrows=nrows, ncols=1, height_ratios=height_ratios, right=0.95, bottom=0.05 - ) # nrows for multiple channels - gs.update(hspace=0.17) - gs_row_idx = 0 - - # Register custom matplotlib projections - proj.register_projection(CustomAxes) - proj.register_projection(MultiLineAxes) - proj.register_projection(SequenceDiagram) - proj.register_projection(SpectralDimensionLabels) - - # Spectral Dimension labels Axes - spec_dim_label_ax = fig.add_subplot(gs[gs_row_idx, 0], projection="label_axes") - spec_dim_label_ax.plot_labels(df=df, x_data=x_offset) - gs_row_idx += 1 - - # Sequence diagram Axes - seq_ax = fig.add_subplot(gs[gs_row_idx, 0], projection="sequence_axes") - seq_ax.make_plot(df=df, x_data=x_offset) - gs_row_idx += 1 - - # p and d Axes - p_ax = fig.add_subplot(gs[gs_row_idx, 0], projection="multi_line_axes") - p_ax.make_plot( - x_data=x_offset, - y_data=df["p"], - col_name="p", - mix_ev=mix_ev, - format_kwargs={}, - plot_kwargs={"alpha": 0.6, "linewidth": 2.5}, - ) - gs_row_idx += 1 - - if "d" in params: - params.remove("d") - d_ax = fig.add_subplot(gs[gs_row_idx, 0], projection="multi_line_axes") - d_ax.make_plot( - x_data=x_offset, - y_data=df["d"], - col_name="d", - mix_ev=mix_ev, - format_kwargs={}, - plot_kwargs={"alpha": 0.6, "linewidth": 2.5}, - ) - gs_row_idx += 1 - - # params Axes - for i, param in enumerate(params, gs_row_idx): - ax = fig.add_subplot(gs[i, 0], projection="custom_axes") - ax.make_plot( - x_data=x_data, - y_data=df[param], - col_name=param, - mix_ev=mix_ev, - format_kwargs={}, - plot_kwargs={"alpha": 0.6, "linewidth": 2.5}, - ) - - # Add legend for event colors - if include_legend: - _add_legend(fig) - - return fig diff --git a/src/mrsimulator/method/query.py b/src/mrsimulator/method/query.py index d5da01f0d..2e942122d 100644 --- a/src/mrsimulator/method/query.py +++ b/src/mrsimulator/method/query.py @@ -260,7 +260,6 @@ class Config: class MixingQuery(Parseable): - # TODO: Fix example code """MixingQuery class for quering transition mixing between events. Attributes @@ -278,7 +277,7 @@ class MixingQuery(Parseable): Example ------- - >>> query = MixingQuery(ch1={"tip_angle": 1.570796, "phase": 3.141593}) + >>> query = TransitionQuery(ch1={'P': [1], 'D': [0]}, ch2={'P': [-1]}) """ ch1: RFRotation = None diff --git a/src/mrsimulator/method/spectral_dimension.py b/src/mrsimulator/method/spectral_dimension.py index 0a36b91dc..183b6650a 100644 --- a/src/mrsimulator/method/spectral_dimension.py +++ b/src/mrsimulator/method/spectral_dimension.py @@ -22,16 +22,6 @@ CHANNELS = ["ch1", "ch2", "ch3"] -class Reciprocal(Parseable): - """Reciprocal dimension from CSDM object.""" - - coordinates_offset: float = 0.0 - - property_unit_types: ClassVar[Dict] = {"coordinates_offset": "time"} - property_default_units: ClassVar[Dict] = {"coordinates_offset": "s"} - property_units: Dict = {"coordinates_offset": "s"} - - class SpectralDimension(Parseable): r"""Base SpectralDimension class defines a spectroscopic dimension of the method. @@ -66,11 +56,11 @@ class SpectralDimension(Parseable): events: A list of :ref:`event_api` or equivalent dict objects (optional). The value describes a series of events along the spectroscopic dimension. """ + count: int = Field(1024, gt=0) spectral_width: float = Field(default=25000.0, gt=0) reference_offset: float = Field(default=0.0) origin_offset: float = None - reciprocal: Reciprocal = None events: List[Union[MixingEvent, ConstantDurationEvent, SpectralEvent]] = [] property_unit_types: ClassVar[Dict] = { @@ -96,8 +86,8 @@ class Config: @classmethod def parse_dict_with_units(cls, py_dict: dict): - """Parse the physical quantities of a SpectralDimension object from a python - dictionary object. + """Parse the physical quantities of a SpectralDimension object from a + python dictionary object. Args: dict py_dict: Dict object @@ -120,7 +110,8 @@ def coordinates_Hz(self) -> np.ndarray: x_\text{Hz} = \left([0, 1, ... N-1] - T\right) \frac{\Delta x}{N} + x_0 where :math:`T=N/2` and :math:`T=(N-1)/2` for even and odd values of - :math:`N`, respectively.""" + :math:`N`, respectively. + """ n = self.count Tk = int(n / 2) increment = self.spectral_width / self.count @@ -133,7 +124,8 @@ def coordinates_ppm(self) -> np.ndarray: .. math:: x_\text{ppm} = \frac{x_\text{Hz}} {x_0 + \omega_0} - where :math:`\omega_0` is the Larmor frequency.""" + where :math:`\omega_0` is the Larmor frequency. + """ if self.origin_offset is None: warnings.warn( UserWarning( @@ -151,12 +143,6 @@ def to_csdm_dimension(self) -> cp.Dimension: increment = self.spectral_width / self.count label = "" if self.label is None else self.label description = "" if self.description is None else self.description - - default_reciprocal = {"coordinates_offset": f"{-1/(2*increment)} s"} - reciprocal = ( - default_reciprocal if self.reciprocal is None else self.reciprocal.json() - ) - dim = cp.Dimension( type="linear", count=self.count, @@ -165,26 +151,12 @@ def to_csdm_dimension(self) -> cp.Dimension: label=label, description=description, complex_fft=True, - reciprocal=reciprocal, + reciprocal={"coordinates_offset": f"{-1/(2*increment)} s"}, ) if self.origin_offset is not None: dim.origin_offset = f"{self.origin_offset} Hz" return dim - # def events_to_dataframe(self) -> pd.DataFrame: - # """Returns events list as DataFrame with event number as columns""" - # attributes = list(Event().property_units.keys()) - # attributes.append("fraction") - # rows = attributes.copy() - # rows.extend(["p", "d"]) - # df = pd.DataFrame(index=rows) - # for i in range(len(self.events)): - # _lst = [getattr(self.events[i], att) for att in attributes] - # _lst.append(self.events[i].transition_query.get_p()) - # _lst.append(self.events[i].transition_query.get_d()) - # df[i] = _lst - # return df - def _get_symmetry_pathways(self, symmetry_element: str) -> list: """Generate a list of symmetry pathways for the event. diff --git a/src/mrsimulator/method/tests/test_method_utils.py b/src/mrsimulator/method/tests/test_method_utils.py index 02ba729cc..57f948a84 100644 --- a/src/mrsimulator/method/tests/test_method_utils.py +++ b/src/mrsimulator/method/tests/test_method_utils.py @@ -1,80 +1,10 @@ # -*- coding: utf-8 -*- from mrsimulator import Site from mrsimulator import SpinSystem -from mrsimulator.method import SpectralDimension -from mrsimulator.method.event import MixingEvent -from mrsimulator.method.query import MixingQuery -from mrsimulator.method.utils import mixing_query_connect_map -from mrsimulator.method.utils import nearest_nonmixing_event -from mrsimulator.method.utils import tip_angle_and_phase_list from mrsimulator.methods import Method1D -__author__ = "Deepansh Srivastava" -__email__ = "srivastava.89@osu.edu" - def test_warnings(): s = SpinSystem(sites=[Site(isotope="23Na")]) m = Method1D(channels=["1H"]) assert m.get_transition_pathways(s) == [] - - -ME = "MixingEvent" -SE = "SpectralEvent" -CE = "ConstantDurationEvent" - - -def test_nearest_mixing_query(): - events = [SE, ME, SE, SE, CE, ME, SE, ME, ME, CE] - assert nearest_nonmixing_event(events, 1) == [0, 2] - assert nearest_nonmixing_event(events, 5) == [4, 6] - assert nearest_nonmixing_event(events, 7) == [6, 9] - - -def test_tip_angle_and_phase_list(): - symbol = ["29Si", "17O", "29Si", "29Si", "27Al"] - channels = ["29Si", "27Al"] - mixing_query = MixingQuery(ch1={"tip_angle": 1.55, "phase": 0.5}) - lst = tip_angle_and_phase_list(symbol, channels, mixing_query) - assert lst == ([1.55, 0, 1.55, 1.55, 0], [0.5, 0, 0.5, 0.5, 0]) - - mixing_query = MixingQuery( - ch1={"tip_angle": 1.55, "phase": 0.5}, - ch2={"tip_angle": 3.1415, "phase": 1.570}, - ) - lst = tip_angle_and_phase_list(symbol, channels, mixing_query) - assert lst == ([1.55, 0, 1.55, 1.55, 3.1415], [0.5, 0, 0.5, 0.5, 1.57]) - - -def test_mixing_query_connect_map(): - MX1 = MixingEvent(mixing_query={"ch1": {"tip_angle": 0.12}}) - MX2 = MixingEvent(mixing_query={"ch2": {"tip_angle": 1.12}}) - spectral_dimmensions = [ - SpectralDimension( - events=[ - {"fraction": 0.5}, # 0 - MX1, - {"duration": 0.5}, # 1 - {"fraction": 0.5}, # 2 - MX2, - {"duration": 0.5}, # 3 - ] - ), - SpectralDimension( - events=[ - MX1, - MX2, - {"duration": 0.5}, # 4 - MX2, - {"duration": 0.5}, # 5 - ] - ), - ] - res = mixing_query_connect_map(spectral_dimmensions) - assert res == [ - {"mixing_query": MX1.mixing_query, "near_index": [0, 1]}, - {"mixing_query": MX2.mixing_query, "near_index": [2, 3]}, - {"mixing_query": MX1.mixing_query, "near_index": [3, 4]}, - {"mixing_query": MX2.mixing_query, "near_index": [3, 4]}, - {"mixing_query": MX2.mixing_query, "near_index": [4, 5]}, - ] diff --git a/src/mrsimulator/method/tests/test_plot.py b/src/mrsimulator/method/tests/test_plot.py deleted file mode 100644 index 61bc85b02..000000000 --- a/src/mrsimulator/method/tests/test_plot.py +++ /dev/null @@ -1,412 +0,0 @@ -# -*- coding: utf-8 -*- -import matplotlib.projections as proj -import numpy as np -import pandas as pd -import pytest -from matplotlib import pyplot as plt -from mrsimulator.method import Method -from mrsimulator.method.plot import _add_tip_angle_and_phase -from mrsimulator.method.plot import _check_columns -from mrsimulator.method.plot import _format_df -from mrsimulator.method.plot import _make_normal_and_offset_x_data -from mrsimulator.method.plot import _make_x_data -from mrsimulator.method.plot import _offset_x_data -from mrsimulator.method.plot import CustomAxes -from mrsimulator.method.plot import MultiLineAxes -from mrsimulator.method.plot import SequenceDiagram - - -__author__ = "Matthew D. Giammar" -__email__ = "giammar.7@buckeyemail.osu.edu" - - -def test_check_columns(): - df1 = method1_df() - df2 = method2_df() - - # Test required columns present in DataFrame - bad_df1 = df1.drop("fraction", axis=1) - bad_df2 = df2.drop(["p"], axis=1) - error = r".*Some required columns were not present in the DataFrame.*" - with pytest.raises(ValueError, match=error): - _check_columns(bad_df1) - with pytest.raises(ValueError, match=error): - _check_columns(bad_df2) - - # Test correct columns dropped from DataFrame - assert _check_columns(df1) == ["d"] - assert _check_columns(df2) == [ - "magnetic_flux_density", - "rotor_frequency", - "rotor_angle", - "d", - ] - - -def test_add_tip_angle_and_phase(): - df1 = method1_df() - df2 = method2_df() - - # Modify the dataframes - _add_tip_angle_and_phase(df1) - _add_tip_angle_and_phase(df2) - - # Check 'tip_angle' and 'phase' are columns now in DataFrame - assert "tip_angle" in df1.columns - assert "phase" in df1.columns - assert "tip_angle" in df2.columns - assert "phase" in df2.columns - - # Check correct calculations for tip angle and phase - ta_should_be = np.array([np.nan, np.nan, 90.0, np.nan]) - p_should_be = np.array([np.nan, np.nan, 0.0, np.nan]) - - assert np.allclose(np.asarray(df1["tip_angle"]), ta_should_be, equal_nan=True) - assert np.allclose(np.asarray(df1["phase"]), p_should_be, equal_nan=True) - - ta_should_be = np.array( - [ - 90.0, - np.nan, - 180.0, - 270.0, - np.nan, - np.nan, - np.nan, - np.nan, - np.nan, - 68.75493541569878, - 42.97183463481174, - 18.90760723931717, - np.nan, - np.nan, - np.nan, - ] - ) - p_should_be = np.array( - [ - 17.188733853924695, - np.nan, - 8.594366926962348, - 0.0, - np.nan, - np.nan, - np.nan, - np.nan, - np.nan, - 180.0, - 45.0, - 60.0, - np.nan, - np.nan, - np.nan, - ] - ) - - assert np.allclose(np.asarray(df2["tip_angle"]), ta_should_be, equal_nan=True) - assert np.allclose(np.asarray(df2["phase"]), p_should_be, equal_nan=True) - - -def test_format_df(): - df1 = method1_df() - df2 = method2_df() - params1 = _format_df(df1) - params2 = _format_df(df2) - - # 'd' should not be in params1 since 'd' is never defined - assert "d" not in params1 - - # Check 'tip_angle' and 'phase' are columns now in DataFrame - assert "tip_angle" in df1.columns - assert "phase" in df1.columns - assert "tip_angle" in df2.columns - assert "phase" in df2.columns - - # Check expected params returned - assert isinstance(params1, list) - assert isinstance(params2, list) - - -def test_make_x_data(): - df1 = method1_df() - df2 = method2_df() - x1 = _make_x_data(df1) - x2 = _make_x_data(df2) - - # Check return type - assert isinstance(x1, list) - assert isinstance(x2, list) - - # Check expected x_data returned - # NOTE: Should arrays be hardcoded? Or should be calculated in similar way - x1_should_be = [0, 0.8, 0.8, 1.3, 1.3, 1.7] - x2_should_be = [ - 0, - 0.8, - 0.8, - 1.3, - 1.3, - 1.46, - 1.46, - 2.1, - 2.1, - 2.6, - 2.6, - 3.1, - 3.1, - 3.5, - 3.5, - 4.0, - 4.0, - 4.4, - ] - - assert np.allclose(x1, x1_should_be) - assert np.allclose(x2, x2_should_be) - - -def test_offset_x_data(): - # Setup DataFrames - df1 = method1_df() - df2 = method2_df() - df3 = method3_df() - _add_tip_angle_and_phase(df1) - _add_tip_angle_and_phase(df2) - - # Calculate needed x_data - x1 = _make_x_data(df1) - x2 = _make_x_data(df2) - x3 = _make_x_data(df3) - off_x1 = _offset_x_data(df1, x1) - off_x2 = _offset_x_data(df2, x2) - off_x3 = _offset_x_data(df3, x3) - - # Check return type - assert isinstance(off_x1, np.ndarray) - assert isinstance(off_x2, np.ndarray) - assert isinstance(off_x3, np.ndarray) - - # Check expected offset_x returned - # NOTE: Should arrays be hardcoded? Or should be calculated in similar way - off_x1_should_be = [0.0, 0.0, 0.8, 0.8, 1.26875, 1.33125, 1.7] - off_x2_should_be = [ - 0.0, - 0.0625, - 0.64375, - 0.95625, - 1.3, - 1.3, - 1.46, - 1.46, - 2.1, - 2.1, - 2.6, - 2.6, - 3.0546408412188097, - 3.1453591587811904, - 3.5, - 3.5, - 4.0, - 4.0, - 4.4, - ] - off_x3_should_be = [0, 0, 0.08, 0.08, 0.8] - - assert np.allclose(off_x1, off_x1_should_be) - assert np.allclose(off_x2, off_x2_should_be) - assert np.allclose(off_x3, off_x3_should_be) - - -def test_make_normal_and_offset_x_data(): - # Setup DataFrames - df1 = method1_df() - df2 = method2_df() - _add_tip_angle_and_phase(df1) - _add_tip_angle_and_phase(df2) - - # Test when x_data is length zero - no_events_df1 = df1.drop(df1.index, axis=0) - no_events_df2 = df2.drop(df2.index[1:], axis=0) - - error = ( - r".*The DataFrame does not contain any SpectralEvents or " - r"ConstandDurationEvents. At least one must be present to construct a plot.*" - ) - with pytest.raises(ValueError, match=error): - _make_normal_and_offset_x_data(no_events_df1) - with pytest.raises(ValueError, match=error): - _make_normal_and_offset_x_data(no_events_df2) - - # Check return types - x1, off_x1 = _make_normal_and_offset_x_data(df1) - x2, off_x2 = _make_normal_and_offset_x_data(df2) - - assert isinstance(x1, list) - assert isinstance(x2, list) - assert isinstance(off_x1, np.ndarray) - assert isinstance(off_x2, np.ndarray) - - -def test_SequenceDiagram(): - df3 = method3_df() - _add_tip_angle_and_phase(df3) - x3 = [0.0, 0.0, 0.08, 0.08, 0.8] - x_data_should_be = [0, 0, 0.08, 0.08, 0.675, 0.8] - - # Setup matplotlib objects - fig = plt.figure() - proj.register_projection(SequenceDiagram) - axes_obj = fig.add_subplot(projection="sequence_axes") - axes_obj.make_plot(df3, x3) - - assert np.allclose(axes_obj.x_data, x_data_should_be) - - -def test_MultilineAxes(): - x_data = [0, 0, 1, 1, 2, 2, 3, 3] - y_data = pd.Series([[1, 1, 1], [1, 0, 0], [1, 0, -1], [1, 1, 0]]) - - # Setup matplotlib objects - fig = plt.figure() - proj.register_projection(MultiLineAxes) - axes_obj = fig.add_subplot(projection="multi_line_axes") - - # make_plot - # Test y_data dimensions check - bad_y_data = pd.Series([0, 1, 2, 3]) # not 2 dimensional - error = r".*Symmetry pathway data is misshapen. Data must be 2d.*" - with pytest.raises(ValueError, match=error): - axes_obj.make_plot(x_data, bad_y_data, "", [], {}, {}) - - # _offset_overlaps - # Test only one symmetry pathway - y_data = np.zeros((1, 10)) - assert np.array_equal(axes_obj._offset_overlaps(y_data), y_data) - - # Test overlapping symmetry pathways - # Each column represents a symmetry pathway - y_data = np.array( - [ - np.array([1, 1, 1, 1], dtype=float), - np.array([1, 0, 0, 1], dtype=float), - np.array([1, 0, -1, 0], dtype=float), - ] - ) - offset_should_be = [ - [0.91, 0.97, 1.0, 0.94], - [1.0, -0.06, 0.0, 1.06], - [1.09, 0.06, -1.0, -0.03], - ] - assert np.allclose(axes_obj._offset_overlaps(y_data), offset_should_be) - - -def test_CustomAxes(): - # Setup matplotlib objects - fig = plt.figure() - proj.register_projection(CustomAxes) - axes_obj = fig.add_subplot(projection="custom_axes") - - # _add_rect_with_labels - # Test error thrown when no color supplied - error = r".*No color in `rect_kwargs`. A color must be specified.*" - with pytest.raises(ValueError, match=error): - axes_obj._add_rect_with_label(0, 1, "", {"height": 1}) - - # _add_blank_space_labels - axes_obj.make_plot( - [1, 2, 2, 3, 3, 4], [0, 1e12, np.nan], "rotor_frequency", [False] * 3, {}, {} - ) - - -def method1_df(): - method1 = Method( - channels=["1H"], - spectral_dimensions=[ - { - "label": "Spectral and Constant Duration Event", - "events": [ - {"fraction": 1}, - {"duration": 1.5}, - ], - }, - { - "label": "Mixing and Spectral Event", - "events": [ - {"mixing_query": {"ch1": {"tip_angle": np.pi / 2, "phase": 0}}}, - {"fraction": 0.5}, - ], - }, - ], - ) - - return method1.summary() - - -def method2_df(): - method2 = Method( - channels=["1H"], - spectral_dimensions=[ - { - "label": "Mixing, Spectral, Mixing, Mixing, ConstantDuration ", - "events": [ - {"mixing_query": {"ch1": {"tip_angle": np.pi / 2, "phase": 0.3}}}, - {"fraction": 1, "transition_query": [{"P": [1], "D": [0]}]}, - {"mixing_query": {"ch1": {"tip_angle": np.pi, "phase": 0.15}}}, - {"mixing_query": {"ch1": {"tip_angle": 3 * np.pi / 2, "phase": 0}}}, - {"duration": 1.5, "transition_query": [{"P": [-1], "D": [2]}]}, - ], - }, - { - "label": "Spectral, Spectral, ConstantDuration, ConstantDuration", - "events": [ - {"fraction": 0.2}, - {"fraction": 0.8}, - {"duration": 1.3}, - {"duration": 0.6}, - ], - }, - { - "label": "Mixing, Mixing, Mixing, Spectral, ConstantDuration, Spectral", - "events": [ - {"mixing_query": {"ch1": {"tip_angle": 1.2, "phase": np.pi}}}, - {"mixing_query": {"ch1": {"tip_angle": 0.75, "phase": np.pi / 4}}}, - {"mixing_query": {"ch1": {"tip_angle": 0.33, "phase": np.pi / 3}}}, - {"fraction": 0.5}, - {"duration": 1.5}, - {"fraction": 0.5}, - ], - }, - ], - ) - - return method2.summary(drop_constant_columns=False) - - -def method3_df(): - method3 = Method( - channels=["1H"], - spectral_dimensions=[ - { - "label": "Spectral Spectral Mixing", - "events": [ - { - "fraction": 0.1, - "transition_query": [ - {"ch1": {"P": [1, 1]}}, - {"ch1": {"P": [2]}}, - ], - }, - { - "fraction": 0.9, - "transition_query": [ - {"ch1": {"P": [-1, -1]}}, - {"ch1": {"P": [-2]}}, - ], - }, - {"mixing_query": {"ch1": {"tip_angle": np.pi, "phase": 0}}}, - ], - } - ], - ) - - return method3.summary() diff --git a/src/mrsimulator/method/tests/test_summary.py b/src/mrsimulator/method/tests/test_summary.py deleted file mode 100644 index abf4fe2ad..000000000 --- a/src/mrsimulator/method/tests/test_summary.py +++ /dev/null @@ -1,303 +0,0 @@ -# -*- coding: utf-8 -*- -import numpy as np -import pandas as pd -import pytest -from matplotlib.pyplot import Figure -from mrsimulator.method import Method -from mrsimulator.method import SpectralDimension - -__author__ = "Matthew D. Giammar" -__email__ = "giammar.7@buckeyemail.osu.edu" - - -ME = "MixingEvent" -CDE = "ConstantDurationEvent" -SE = "SpectralEvent" -REQUIRED = [ - "type", - "label", - "duration", - "fraction", - "mixing_query", - "spec_dim_index", - "spec_dim_label", - "p", - "d", -] -ALL_PARAMS = [ - "type", - "spec_dim_index", - "spec_dim_label", - "label", - "duration", - "fraction", - "mixing_query", - "magnetic_flux_density", - "rotor_frequency", - "rotor_angle", - "freq_contrib", - "p", - "d", -] - - -def test_empty_spec_dims(): - # Empty list - empty_method = Method(channels=["1H"], spectral_dimensions=[]) - error = ( - r".*Method has empty spectral_dimensions. At least one SpectralDimension " - r"is needed with at least one Event..*" - ) - with pytest.raises(AttributeError, match=error): - empty_method.summary() - - # No events - empty_method.spectral_dimensions = [SpectralDimension(events=[])] - error = ( - r".*Method has no Events. At least one SpectralDimension " - r"is needed with at least one Event..*" - ) - with pytest.raises(AttributeError, match=error): - empty_method.summary() - - -def basic_summary_tests(the_method): - def check_col_equal(col, arr): - check = [] - for col_item, arr_item in zip(col, arr): - if np.isnan(col_item): - check.append(np.isnan(arr_item)) - else: - check.append(np.isclose(col_item, arr_item)) - # Returns True if all items are equal, False otherwise - return np.all(check) - - def check_col_equal_2d(col, arr): - check = [] - for col_list, arr_list in zip(col, arr): - for col_item, arr_item in zip(col_list, arr_list): - if np.isnan(col_item): - check.append(np.isnan(arr_item)) - else: - check.append(np.isclose(col_item, arr_item)) - # Returns True if all items are equal, False otherwise - return np.all(check) - - df = the_method.summary(drop_constant_columns=False) - - # Check correct return type - assert isinstance(df, pd.DataFrame) - - # Check for correct number of event and properties (df.shape) - assert df.shape[0] == 6 - assert df.shape[1] == len(ALL_PARAMS) - - # Check for correct columns - assert set(df.columns) == set(ALL_PARAMS) - - # Check type - assert np.all(df["type"] == [ME, CDE, SE, ME, CDE, SE]) - - # Check label - assert np.all(df["label"] == ["Mix0", "Dur0", "Spec0", "Mix1", "Dur1", "Spec1"]) - - # NOTE: Need to run loops so NANs can be checked since `nan == nan` evals False - # Check duration - temp = [np.nan, 1.0, np.nan, np.nan, 2.0, np.nan] - assert check_col_equal(df["duration"], temp) - - # Check fraction - temp = [np.nan, np.nan, 0.3, np.nan, np.nan, 0.7] - assert check_col_equal(df["fraction"], temp) - - # Check p - temp = [[np.nan], [0.0], [1.0], [np.nan], [3.0], [4.0]] - assert check_col_equal_2d(df["p"], temp) - - # BUG: column d is 2d array of correct size but all nan - # Check d - temp = [[np.nan], [0.0], [-2.0], [np.nan], [-4.0], [-6.0]] - assert check_col_equal_2d(df["d"], temp) - - # Check magnetic_flux_density - temp = [np.nan, 1.0, 2.0, np.nan, 3.0, 4.0] - assert check_col_equal(df["magnetic_flux_density"], temp) - - # Check rotor_frequency - temp = [np.nan, 10.0, 20.0, np.nan, 30.0, 40.0] - assert check_col_equal(df["rotor_frequency"], temp) - - # Check rotor_angle (degrees) - temp = [ - np.nan, - 5.729577951308232, - 11.459155902616464, - np.nan, - 17.188733853924695, - 22.91831180523293, - ] - assert check_col_equal(df["rotor_angle"], temp) - - # TODO: Check mixing_query - # TODO: Check freq_contrib - - -def args_summary_tests(the_method): - # Pass drop_constant_columns=False - df = the_method.summary(drop_constant_columns=False) - - assert df.shape[0] == 6 - assert df.shape[1] == len(ALL_PARAMS) - assert set(df.columns) == set(ALL_PARAMS) - - # Pass drop_constant_columns=True (default) - df = the_method.summary(drop_constant_columns=True) - props_should_be = ["rotor_frequency", "freq_contrib"] - - assert df.shape[0] == 6 - assert df.shape[1] == len(REQUIRED + props_should_be) - assert set(df.columns) == set(REQUIRED + props_should_be) - - # Make sure columns always present are not dropped even if constant - event_ = [{"label": "all labels the same", "fraction": 0.2}] - all_SpectralEvents_spec_dims = [{"events": event_ * 5}] - - the_method = Method( - name="all-spectral-method-1", - channels=["1H", "13C"], - spectral_dimensions=all_SpectralEvents_spec_dims, - ) - - df = the_method.summary(drop_constant_columns=True) - - assert df.shape[0] == 5 - assert df.shape[1] == len(REQUIRED + ["freq_contrib"]) - assert set(df.columns) == set(REQUIRED + ["freq_contrib"]) - - -def test_summary(): - # TODO: Make one mixing_query be empty - all_defined_no_constant_spec_dims = [ - { - "events": [ - { - "label": "Mix0", - "mixing_query": { - "ch1": {"tip_angle": np.pi / 4, "phase": np.pi / 2} - }, - }, - { - "label": "Dur0", - "duration": 1, - "magnetic_flux_density": 1, - "rotor_frequency": 10000, # in kHz - "rotor_angle": 0.1, - "transition_query": [{"ch1": {"P": [0], "D": [0]}}], - }, - { - "label": "Spec0", - "fraction": 0.3, - "magnetic_flux_density": 2, - "rotor_frequency": 20000, - "rotor_angle": 0.2, - "transition_query": [{"ch1": {"P": [1], "D": [-2]}}], - }, - ] - }, - { - "events": [ - { - "label": "Mix1", - "mixing_query": {"ch1": {"tip_angle": np.pi / 2, "phase": np.pi}}, - }, - { - "label": "Dur1", - "duration": 2, - "magnetic_flux_density": 3, - "rotor_frequency": 30000, - "rotor_angle": 0.3, - "transition_query": [{"ch1": {"P": [3], "D": [-4]}}], - }, - { - "label": "Spec1", - "fraction": 0.7, - "magnetic_flux_density": 4, - "rotor_frequency": 40000, - "rotor_angle": 0.4, - "transition_query": [{"ch1": {"P": [4], "D": [-6]}}], - }, - ] - }, - ] - - const_mfd_rotor_ang_spec_dims = [ - { - "events": [ - { - "label": "Mix0", - "mixing_query": {"ch1": {"tip_angle": np.pi / 2, "phase": np.pi}}, - }, - { - "label": "Dur0", - "duration": 1, - "rotor_frequency": 10000, - "rotor_angle": 0.5, - "transition_query": [{"ch1": {"P": [0], "D": [0]}}], - }, - { - "label": "Spec0", - "fraction": 0.3, - "rotor_frequency": 20000, - "rotor_angle": 0.5, - "transition_query": [{"ch1": {"P": [1], "D": [0]}}], - }, - ] - }, - { - "events": [ - { - "label": "Mix1", - "mixing_query": {"ch1": {"tip_angle": np.pi / 2, "phase": np.pi}}, - }, - { - "label": "Dur1", - "duration": 2, - "rotor_frequency": 30000, - "rotor_angle": 0.5, - "transition_query": [{"ch1": {"P": [3], "D": [0]}}], - }, - { - "label": "Spec1", - "fraction": 0.7, - "rotor_frequency": 40000, - "rotor_angle": 0.5, - "transition_query": [{"ch1": {"P": [4], "D": [0]}}], - }, - ] - }, - ] - - method1 = Method( - name="test-method-1", - channels=["1H", "13C"], - spectral_dimensions=all_defined_no_constant_spec_dims, - ) - - # Constant 'magnetic_flux_density' and 'rotor_frequency' - method2 = Method( - name="test-method-2", - channels=["1H", "13C"], - spectral_dimensions=const_mfd_rotor_ang_spec_dims, - ) - - assert method1.name == "test-method-1" - assert len(method1.spectral_dimensions) == 2 - - basic_summary_tests(method1) - - assert method2.name == "test-method-2" - assert len(method2.spectral_dimensions) == 2 - - args_summary_tests(method2) - - assert isinstance(method2.plot(), Figure) diff --git a/src/mrsimulator/method/tests/test_transitions.py b/src/mrsimulator/method/tests/test_transitions.py index 5c23d2b4c..a1a44f57b 100644 --- a/src/mrsimulator/method/tests/test_transitions.py +++ b/src/mrsimulator/method/tests/test_transitions.py @@ -1,9 +1,6 @@ # -*- coding: utf-8 -*- -import numpy as np from mrsimulator import SpinSystem from mrsimulator.method import Method -from mrsimulator.methods import Method1D -from mrsimulator.methods import Method2D from mrsimulator.transition import TransitionPathway __author__ = "Deepansh J. Srivastava" @@ -21,7 +18,7 @@ def check_transition_set(got, expected): - assert len(got) == len(expected), "Inconsistent transition pathway count" + assert len(got) == len(expected), "Inconsisten transition pathway count" for item in got: assert item in expected, f"Transition pathways not found: {item}" @@ -86,100 +83,3 @@ def test_04(): TransitionPathway([{"final": [0.5, -0.5, 1], "initial": [0.5, 0.5, 1]}]), ] check_transition_set(tr, expected) - - -def test_hahn(): - system = SpinSystem(sites=[{"isotope": "13C"}, {"isotope": "13C"}]) - hahn = Method1D( - channels=["13C"], - spectral_dimensions=[ - { - "events": [ - {"fraction": 0.5, "transition_query": {"P": [1]}}, - {"mixing_query": {"ch1": {"tip_angle": np.pi, "phase": 0}}}, - {"fraction": 0.5, "transition_query": {"P": [-1]}}, - ] - }, - ], - ) - tr = hahn.get_transition_pathways(system) - - weights = np.asarray([1, 1, 1, 1]) - transition_pathways = 0.5 * np.asarray( - [ - [[[-1, -1], [-1, 1]], [[1, 1], [1, -1]]], - [[[-1, -1], [1, -1]], [[1, 1], [-1, 1]]], - [[[-1, 1], [1, 1]], [[1, -1], [-1, -1]]], - [[[1, -1], [1, 1]], [[-1, 1], [-1, -1]]], - ] - ) - - expected = [ - TransitionPathway( - pathway=[ - {"initial": list(states[0]), "final": list(states[1])} - for states in transitions - ], - weight=w, - ) - for transitions, w in zip(transition_pathways, weights) - ] - - assert tr == expected - - -def test_cosy(): - system = SpinSystem(sites=[{"isotope": "1H"}, {"isotope": "1H"}]) - cosy = Method2D( - channels=["1H"], - spectral_dimensions=[ - { - "events": [ - {"fraction": 0.5, "transition_query": {"P": [-1]}}, - {"mixing_query": {"ch1": {"tip_angle": np.pi / 2, "phase": 0}}}, - ], - }, - { - "events": [{"fraction": 0.5, "transition_query": {"P": [-1]}}], - }, - ], - ) - tr = cosy.get_transition_pathways(system) - - weights = np.asarray([1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1, -1, 1, 1, -1, 1]) * 0.25 - transition_pathways = 0.5 * np.asarray( - [ - [[[-1, 1], [-1, -1]], [[-1, 1], [-1, -1]]], - [[[-1, 1], [-1, -1]], [[1, -1], [-1, -1]]], - [[[-1, 1], [-1, -1]], [[1, 1], [-1, 1]]], - [[[-1, 1], [-1, -1]], [[1, 1], [1, -1]]], - # - [[[1, -1], [-1, -1]], [[-1, 1], [-1, -1]]], - [[[1, -1], [-1, -1]], [[1, -1], [-1, -1]]], - [[[1, -1], [-1, -1]], [[1, 1], [-1, 1]]], - [[[1, -1], [-1, -1]], [[1, 1], [1, -1]]], - # - [[[1, 1], [-1, 1]], [[-1, 1], [-1, -1]]], - [[[1, 1], [-1, 1]], [[1, -1], [-1, -1]]], - [[[1, 1], [-1, 1]], [[1, 1], [-1, 1]]], - [[[1, 1], [-1, 1]], [[1, 1], [1, -1]]], - # - [[[1, 1], [1, -1]], [[-1, 1], [-1, -1]]], - [[[1, 1], [1, -1]], [[1, -1], [-1, -1]]], - [[[1, 1], [1, -1]], [[1, 1], [-1, 1]]], - [[[1, 1], [1, -1]], [[1, 1], [1, -1]]], - ] - ) - - expected = [ - TransitionPathway( - pathway=[ - {"initial": list(states[0]), "final": list(states[1])} - for states in transitions - ], - weight=w, - ) - for transitions, w in zip(transition_pathways, weights) - ] - - assert tr == expected diff --git a/src/mrsimulator/method/utils.py b/src/mrsimulator/method/utils.py index b561029fb..139cedb24 100644 --- a/src/mrsimulator/method/utils.py +++ b/src/mrsimulator/method/utils.py @@ -52,113 +52,6 @@ def get_iso_dict(channels, isotopes): return {item: (np.where(isotopes == item))[0] for item in intersection} -def nearest_nonmixing_event(event_name, i): - """Return the indexes of the nearest non mixing events (SpectralEvent and - ConstantDurationEvent) about a mixing event at index `i`. - - Args: - event_name: List of event class names. - i: Int index of the mixing event. - """ - options = ["SpectralEvent", "ConstantDurationEvent"] - low_range = event_name[:i] - high_range = event_name[i:] - upper = [high_range.index(item) for item in options if item in high_range] - lower = [low_range[::-1].index(item) for item in options if item in low_range] - return [ - i - 1 - (min(lower) if lower != [] else 0), - i + (min(upper) if upper != [] else 0), - ] - - -def tip_angle_and_phase_list(symbol, channels, mixing_query): - """Return a list of tip_angles and phase of size equal to the number of sites within - the spin system, corresponding to a mixing_query from a MixingEvent. - - If the site matches the channel, append the tip_angle and phase of the corresponding - channel to the list, otherwise append 0. - - Args: - symbols: List of site symbols. - channels: List of method channel symbols. - mixing_query: Mixing query object of the MixingEvent. - """ - angle_mappable = map_mix_query_attr_to_ch(mixing_query) - tip_angle_ = [ - angle_mappable["tip_angle"][channels.index(sym)] if sym in channels else 0 - for sym in symbol - ] - phase_ = [ - angle_mappable["phase"][channels.index(sym)] if sym in channels else 0 - for sym in symbol - ] - return tip_angle_, phase_ - - -def get_mixing_query(spectral_dimensions, index): - """Return the mixing query object corresponding to the event at index `index`. The - indexing is over flattened list of events from all spectral dimensions. - - Args: - spectral_dimension: A list SpectralDimension objects. - index: The index of the event from a flatten event list. - """ - n_events = len(spectral_dimensions[0].events) - sp = 0 - while index >= n_events: - index -= n_events - sp += 1 - n_events = len(spectral_dimensions[sp].events) - return spectral_dimensions[sp].events[index].mixing_query - - -def map_mix_query_attr_to_ch(mixing_query): - """Map the mixing query attributes (tip_angle and phase) to the channel index. - If the attribute is defined for the channel use the defined value else set it to 0. - - Args: - spectral_dimension: A list SpectralDimension objects. - index: The index of the event from a flatten event list. - """ - attributes = ["tip_angle", "phase"] - return { - item: { - i: getattr(getattr(mixing_query, f"ch{i+1}"), item) or 0 - if getattr(mixing_query, f"ch{i+1}") is not None - else 0 - for i in range(3) - } - for item in attributes - } - - -def mixing_query_connect_map(spectral_dimensions): - """Return a list of mappables corresponding to each mixing event. The mappable - corresponds to mixing event and the index of next nearest transition indexes. - - Args: - spectral_dimensions: A list of SpectralDimension objects.""" - mapping = {} - event_names = [ - evt.__class__.__name__ for dim in spectral_dimensions for evt in dim.events - ] - non_mix_index = [i for i, ev in enumerate(event_names) if ev != "MixingEvent"] - non_mix_index_map = {index: i for i, index in enumerate(non_mix_index)} - - mapping = [ - { - "mixing_query": get_mixing_query(spectral_dimensions, i), - "near_index": [ - non_mix_index_map[k] for k in nearest_nonmixing_event(event_names, i) - ], - } - for i, name in enumerate(event_names) - if name == "MixingEvent" - ] - - return mapping - - # Deprecated # def query_permutations(query, isotope, channel, transition_symmetry="P"): # """ diff --git a/src/mrsimulator/methods/correlation.py b/src/mrsimulator/methods/correlation.py index c36a6508d..5d0e2107c 100644 --- a/src/mrsimulator/methods/correlation.py +++ b/src/mrsimulator/methods/correlation.py @@ -4,8 +4,6 @@ - Cosy - Inadequate """ -import numpy as np - from .base import BaseNamedMethod2D __author__ = "Deepansh J. Srivastava" @@ -50,22 +48,22 @@ class Cosy(BaseNamedMethod2D): ... couplings=[Coupling(site_index=[0, 1], isotropic_j=200)] ... ) >>> pprint(cosy.get_transition_pathways(sys)) - [|-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, weight=(0.25+0j), - |-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, weight=(-0.25-0j), - |-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(0.25+0j), - |-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, weight=(0.25+0j), - |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, weight=(-0.25-0j), - |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, weight=(0.25+0j), - |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(0.25+0j), - |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, weight=(0.25+0j), - |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, weight=(0.25+0j), - |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, weight=(0.25+0j), - |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(0.25+0j), - |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, weight=(-0.25-0j), - |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, weight=(0.25+0j), - |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, weight=(0.25+0j), - |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(-0.25-0j), - |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, weight=(0.25+0j)] + [|-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, + |-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, + |-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, + |-0.5, -0.5⟩⟨-0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, + |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, + |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, + |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, + |-0.5, -0.5⟩⟨0.5, -0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, + |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, + |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, + |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, + |-0.5, 0.5⟩⟨0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, + |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, + |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, + |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, + |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|] """ name: str = "Cosy" @@ -75,10 +73,7 @@ class Cosy(BaseNamedMethod2D): @classmethod def update(self, **kwargs): - event_0 = [ - {"transition_query": [{"ch1": {"P": [-1]}}]}, - {"mixing_query": {"ch1": {"tip_angle": np.pi / 2}}}, - ] + event_0 = [{"transition_query": [{"ch1": {"P": [-1]}}]}] event_1 = [{"transition_query": [{"ch1": {"P": [-1]}}]}] return { @@ -125,10 +120,10 @@ class Inadequate(BaseNamedMethod2D): ... couplings=[Coupling(site_index=[0, 1], isotropic_j=200)] ... ) >>> pprint(inadequate.get_transition_pathways(sys)) - [|-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, weight=(1+0j), - |-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, weight=(1+0j), - |-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(1+0j), - |-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|, weight=(1+0j)] + [|-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨-0.5, 0.5|, + |-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, -0.5⟩⟨0.5, -0.5|, + |-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, + |-0.5, -0.5⟩⟨0.5, 0.5| ⟶ |0.5, -0.5⟩⟨0.5, 0.5|] """ name: str = "Inadequate" diff --git a/src/mrsimulator/methods/mqvas.py b/src/mrsimulator/methods/mqvas.py index 0a0ec0b09..343059c10 100644 --- a/src/mrsimulator/methods/mqvas.py +++ b/src/mrsimulator/methods/mqvas.py @@ -104,7 +104,7 @@ class ThreeQ_VAS(MQ_VAS): ... ) >>> sys = SpinSystem(sites=[Site(isotope='87Rb')]) >>> method.get_transition_pathways(sys) - [|-1.5⟩⟨1.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j)] + [|-1.5⟩⟨1.5| ⟶ |-0.5⟩⟨0.5|] """ @@ -139,7 +139,7 @@ class FiveQ_VAS(MQ_VAS): ... ) >>> sys = SpinSystem(sites=[Site(isotope='17O')]) >>> method.get_transition_pathways(sys) - [|-2.5⟩⟨2.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j)] + [|-2.5⟩⟨2.5| ⟶ |-0.5⟩⟨0.5|] """ @@ -174,5 +174,5 @@ class SevenQ_VAS(MQ_VAS): ... ) >>> sys = SpinSystem(sites=[Site(isotope='51V')]) >>> method.get_transition_pathways(sys) - [|-3.5⟩⟨3.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j)] + [|-3.5⟩⟨3.5| ⟶ |-0.5⟩⟨0.5|] """ diff --git a/src/mrsimulator/methods/ssb2d.py b/src/mrsimulator/methods/ssb2d.py index 60c0af24e..91b7c8b09 100644 --- a/src/mrsimulator/methods/ssb2d.py +++ b/src/mrsimulator/methods/ssb2d.py @@ -35,7 +35,7 @@ class SSB2D(BaseNamedMethod2D): ... ) >>> sys = SpinSystem(sites=[Site(isotope='13C')]) >>> method.get_transition_pathways(sys) - [|-0.5⟩⟨0.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j)] + [|-0.5⟩⟨0.5| ⟶ |-0.5⟩⟨0.5|] """ name: str = "SSB2D" diff --git a/src/mrsimulator/methods/stvas.py b/src/mrsimulator/methods/stvas.py index 37e048514..e2d882816 100644 --- a/src/mrsimulator/methods/stvas.py +++ b/src/mrsimulator/methods/stvas.py @@ -102,8 +102,7 @@ class ST1_VAS(ST_VAS): ... ) >>> sys = SpinSystem(sites=[Site(isotope='87Rb')]) >>> pprint(method.get_transition_pathways(sys)) - [|-1.5⟩⟨-0.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j), - |0.5⟩⟨1.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j)] + [|-1.5⟩⟨-0.5| ⟶ |-0.5⟩⟨0.5|, |0.5⟩⟨1.5| ⟶ |-0.5⟩⟨0.5|] """ @@ -139,6 +138,5 @@ class ST2_VAS(ST_VAS): ... ) >>> sys = SpinSystem(sites=[Site(isotope='17O')]) >>> pprint(method.get_transition_pathways(sys)) - [|-2.5⟩⟨-1.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j), - |1.5⟩⟨2.5| ⟶ |-0.5⟩⟨0.5|, weight=(1+0j)] + [|-2.5⟩⟨-1.5| ⟶ |-0.5⟩⟨0.5|, |1.5⟩⟨2.5| ⟶ |-0.5⟩⟨0.5|] """ diff --git a/src/mrsimulator/simulator/tests/test_simulator.py b/src/mrsimulator/simulator/tests/test_simulator.py index 3a00c54af..9ac13a116 100644 --- a/src/mrsimulator/simulator/tests/test_simulator.py +++ b/src/mrsimulator/simulator/tests/test_simulator.py @@ -327,8 +327,8 @@ def test_sites_to_pandas_df(): eta_q = [None, None, None, 0.34] spin_systems = single_site_system_generator( - isotope=isotopes, - isotropic_chemical_shift=shifts, + isotopes=isotopes, + isotropic_chemical_shifts=shifts, shielding_symmetric={"zeta": zeta, "eta": eta_n}, quadrupolar={"Cq": Cq, "eta": eta_q}, abundance=1, diff --git a/src/mrsimulator/spin_system/__init__.py b/src/mrsimulator/spin_system/__init__.py index 776d3f030..75def4ec8 100644 --- a/src/mrsimulator/spin_system/__init__.py +++ b/src/mrsimulator/spin_system/__init__.py @@ -121,7 +121,7 @@ class SpinSystem(Parseable): ... ] ... ] >>> print(sys1.transition_pathways) - [|2.5, 0.5⟩⟨-2.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(1+0j)] + [|2.5, 0.5⟩⟨-2.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|] .. note:: From any given spin system, the list of relevant transition pathways is diff --git a/src/mrsimulator/spin_system/site.py b/src/mrsimulator/spin_system/site.py index eb6a30f23..304dd0c59 100644 --- a/src/mrsimulator/spin_system/site.py +++ b/src/mrsimulator/spin_system/site.py @@ -183,8 +183,6 @@ def spin_must_be_at_least_one(cls, v, values): @validator("shielding_symmetric", "shielding_antisymmetric") def shielding_symmetric_must_not_contain_Cq_and_D(cls, v, values): - if v is None: - return v _ = [ v.property_units.pop(item) if item in v.property_units else None for item in ["Cq", "D"] diff --git a/src/mrsimulator/spin_system/zeemanstate.py b/src/mrsimulator/spin_system/zeemanstate.py index 8fc019e14..974bdaa47 100644 --- a/src/mrsimulator/spin_system/zeemanstate.py +++ b/src/mrsimulator/spin_system/zeemanstate.py @@ -9,19 +9,22 @@ class ZeemanState: def __init__(self, n_sites, *args): self.n_sites = n_sites - _ = [setattr(self, f"m{i}", args[i]) for i in range(n_sites)] + for i in range(n_sites): + self.__setattr__(f"m{i}", args[i]) def __repr__(self): return self.__str__() # def _repr_html_(self): - # lst = [int(2 * getattr(self, f"m{i}")) for i in range(self.n_sites)] + # lst = [int(2 * self.__getattribute__(f"m{i}")) for i in range(self.n_sites)] # string = " ".join([rf"m_{i}={{{item} \over 2}}" for i, item in enumerate(lst)]) # return rf"$|{string}⟩$" def __str__(self): - lst = ", ".join([f"{getattr(self, f'm{i}')}" for i in range(self.n_sites)]) + lst = ", ".join( + [f"{self.__getattribute__(f'm{i}')}" for i in range(self.n_sites)] + ) return f"|{lst}⟩" def tolist(self): - return [getattr(self, f"m{i}") for i in range(self.n_sites)] + return [self.__getattribute__(f"m{i}") for i in range(self.n_sites)] diff --git a/src/mrsimulator/transition/pathway.py b/src/mrsimulator/transition/pathway.py index 555fd4d64..01e130cd3 100644 --- a/src/mrsimulator/transition/pathway.py +++ b/src/mrsimulator/transition/pathway.py @@ -97,36 +97,28 @@ class TransitionPathway(TransitionList): >>> t2 = Transition(initial=[0.5, 0.5], final=[-0.5, 0.5]) >>> path = TransitionPathway([t1, t2]) >>> path - |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5|, weight=(1+0j) + |0.5, -0.5⟩⟨0.5, 0.5| ⟶ |-0.5, 0.5⟩⟨0.5, 0.5| """ - def __init__(self, pathway: list = [], weight=(1.0 + 0j)): - super().__init__(pathway) - self.weight = weight + # @property + # def weight(self): + # pass def __str__(self): return self.__repr__() - def __eq__(self, other): - return self._list == other._list and self.weight == other.weight - def __repr__(self): - path = " ⟶ ".join([repr(item) for item in self._list]) - return path + f", weight={self.weight}" + return " ⟶ ".join([repr(item) for item in self._list]) def json(self, **kwargs) -> dict: """Parse the class object to a JSON compliant python dictionary object. Example: >>> pprint(path.json()) - {'pathway': [{'final': [0.5, -0.5], 'initial': [0.5, 0.5]}, - {'final': [-0.5, 0.5], 'initial': [0.5, 0.5]}], - 'weight': (1+0j)} + [{'final': [0.5, -0.5], 'initial': [0.5, 0.5]}, + {'final': [-0.5, 0.5], 'initial': [0.5, 0.5]}] """ - return { - "pathway": [item.json() for item in self._list], - "weight": self.weight, - } + return [item.json() for item in self._list] def tolist(self): """Expand TransitionPathway to a Python list. diff --git a/src/mrsimulator/transition/tests/test_transition_pathway.py b/src/mrsimulator/transition/tests/test_transition_pathway.py index 0617497ab..56e7a38b8 100644 --- a/src/mrsimulator/transition/tests/test_transition_pathway.py +++ b/src/mrsimulator/transition/tests/test_transition_pathway.py @@ -21,7 +21,8 @@ def test_transition_pathway(): assert trans_path[2].initial == [0.5, 0.5] assert trans_path[2].final == [-0.5, -0.5] - assert trans_path.json() == {"pathway": [a, b, c], "weight": 1} + # to dict with unit + assert trans_path.json() == [a, b, c] assert trans_path.tolist() == [ 0.5, @@ -41,4 +42,4 @@ def test_transition_pathway(): string = " ⟶ ".join( ["|-0.5, 1.5⟩⟨0.5, 1.5|", "|0.5, 1.5⟩⟨0.5, -1.5|", "|-0.5, -0.5⟩⟨0.5, 0.5|"] ) - assert str(trans_path) == string + ", weight=(1+0j)" + assert str(trans_path) == string diff --git a/src/mrsimulator/utils/collection.py b/src/mrsimulator/utils/collection.py index 4bb99fef6..42da0bc81 100644 --- a/src/mrsimulator/utils/collection.py +++ b/src/mrsimulator/utils/collection.py @@ -1,319 +1,248 @@ # -*- coding: utf-8 -*- -from typing import Dict -from typing import List -from typing import Union - import numpy as np from mrsimulator import Site from mrsimulator import SpinSystem -__author__ = ["Deepansh Srivastava", "Matthew D. Giammar"] -__email__ = ["srivastava.89@osu.edu", "giammar.7@buckeyemail.osu.edu"] +__author__ = "Deepansh Srivastava" +__email__ = "srivastava.89@osu.edu" + SHIELDING_SYM_PARAMS = ["zeta", "eta", "alpha", "beta", "gamma"] QUADRUPOLAR_PARAMS = ["Cq", "eta", "alpha", "beta", "gamma"] -LIST_LEN_ERROR_MSG = ( - "All arguments must be the same size. If one attribute is a type list of length n, " - "then all attributes with list types must also be of length n, and all remaining " - "attributes must be scalar (singular float, int, or str)." -) def single_site_system_generator( - isotope: Union[str, List[str]], - isotropic_chemical_shift: Union[float, List[float], np.ndarray] = 0, - shielding_symmetric: Dict = None, - shielding_antisymmetric: Dict = None, - quadrupolar: Dict = None, - abundance: Union[float, List[float], np.ndarray] = None, - site_name: Union[str, List[str]] = None, - site_label: Union[str, List[str]] = None, - site_description: Union[str, List[str]] = None, - rtol: float = 1e-3, -) -> List[SpinSystem]: + isotopes, + isotropic_chemical_shifts=0, + shielding_symmetric=None, + quadrupolar=None, + abundance=None, + site_labels=None, + site_names=None, + site_descriptions=None, + rtol=1e-3, +): r"""Generate and return a list of single-site spin systems from the input parameters. Args: - isotope: - A required string or a list of site isotopes. - isotropic_chemical_shift: - A float or a list/ndarray of isotropic chemical shifts per site per spin - system. The default is 0. + + isotopes: + A string or a list of site isotopes. + isotropic_chemical_shifts: + A float or a list/ndarray of values. The default value is 0. shielding_symmetric: - A shielding symmetric dict object, where the keyword value can either + A shielding symmetric like dict object, where the keyword value can either be a float or a list/ndarray of floats. The default value is None. The allowed keywords are ``zeta``, ``eta``, ``alpha``, ``beta``, and ``gamma``. - shielding_antisymmetric: - A shielding antisymmetric dict object, where the keyword value can either - be a float or a list/ndarray of floats. The default value is None. The - allowed keywords are ``zeta``, ``alpha``, and ``beta``. quadrupolar: - A quadrupolar dict object, where the keyword value can either be a float or - a list/ndarray of floats. The default value is None. The allowed keywords - are ``Cq``, ``eta``, ``alpha``, ``beta``, and ``gamma``. + A quadrupolar like dict object, where the keyword value can either be a + float or a list/ndarray of floats. The default value is None. The allowed + keywords are ``Cq``, ``eta``, ``alpha``, ``beta``, and ``gamma``. abundance: A float or a list/ndarray of floats describing the abundance of each spin system. - site_name: - A string or a list of strings with site names per site per spin system. The - default is None. - site_label: - A string or a list of strings with site labels per site per spin system. The - default is None. - site_description: - A string or a list of strings with site descriptions per site per spin - system. The default is None. + site_labels: + A string or a list of strings each with a site label. The default is None. + site_names: + A string or a list of strings each with a site name. The default is None. + site_descriptions: + A string or a list of strings each with a site description. Default is None. rtol: - The relative tolerance used in determining the cutoff abundance, given as, + The relative tolerance. This value is used in determining the cutoff + abundance given as :math:`\tt{abundance}_{\tt{cutoff}} = \tt{rtol} * \tt{max(abundance)}.` The spin systems with abundance below this threshold are ignored. - Returns: - List of :ref:`spin_sys_api` objects with a single :ref:`site_api` - - Example: - **Single spin system:** - - >>> sys1 = single_site_system_generator( - ... isotope=["1H"], - ... isotropic_chemical_shift=10, - ... site_name="Single Proton", - ... ) - >>> print(len(sys1)) - 1 - - **Multiple spin system:** - - >>> sys2 = single_site_system_generator( - ... isotope="1H", - ... isotropic_chemical_shift=[10] * 5, - ... site_name="5 Protons", - ... ) - >>> print(len(sys2)) - 5 - - **Multiple spin system with dictionary arguments:** - - >>> Cq = [4.2e6] * 12 - >>> sys3 = single_site_system_generator( - ... isotope="17O", - ... isotropic_chemical_shift=60.0, # in ppm, - ... quadrupolar={"Cq": Cq, "eta": 0.5}, # Cq in Hz - ... ) - >>> print(len(sys3)) - 12 - .. note:: The parameter value can either be a float or a list/ndarray. If the parameter value is a float, the given value is assigned to the respective parameter in all - the spin systems. If the parameter value is a list or ndarray, its `ith` value - is assigned to the respective parameter of the `ith` spin system. When multiple + the spin systems. If the parameter value is a list or ndarray, its ith value is + assigned to the respective parameter of the ith spin system. When multiple parameter values are given as lists/ndarrays, the length of all the lists must be the same. """ - sites = site_generator( - isotope=isotope, - isotropic_chemical_shift=isotropic_chemical_shift, - shielding_symmetric=shielding_symmetric, - shielding_antisymmetric=shielding_antisymmetric, - quadrupolar=quadrupolar, - name=site_name, - label=site_label, - description=site_description, + isotopes = _fix_item(isotopes) + n_isotopes = _get_length(isotopes) + + isotropic_chemical_shifts = _fix_item(isotropic_chemical_shifts) + n_iso = _get_length(isotropic_chemical_shifts) + + abundance = _fix_item(abundance) + n_abd = _get_length(abundance) + + site_labels = _fix_item(site_labels) + n_site_labels = _get_length(site_labels) + + site_names = _fix_item(site_names) + n_site_names = _get_length(site_names) + + site_descriptions = _fix_item(site_descriptions) + n_site_descriptions = _get_length(site_descriptions) + + n_ss, shield_keys, shielding_symmetric = _get_shielding_info(shielding_symmetric) + n_q, quad_keys, quadrupolar = _get_quad_info(quadrupolar) + + n_len = _check_size( + np.asarray( + [ + n_isotopes, + n_iso, + *n_ss, + *n_q, + n_abd, + n_site_labels, + n_site_names, + n_site_descriptions, + ] + ) ) - n_sites = len(sites) - abundance = 1 / n_sites if abundance is None else abundance - abundance = _extend_to_nparray(_fix_item(abundance), n_sites) - n_abd = abundance.size - - if n_sites == 1: - sites = np.asarray([sites[0] for _ in range(n_abd)]) - n_sites = sites.size - - if n_sites != n_abd: - raise ValueError( - "Number of sites does not match the number of abundances. " - f"{LIST_LEN_ERROR_MSG}" + if abundance is None: + abundance = 1 / n_len + + # system with only isotope and isotropic shifts parameters + isotopes_ = _get_default_lists(isotopes, n_len) + iso_chemical_shifts_ = _get_default_lists(isotropic_chemical_shifts, n_len) + abundance_ = _get_default_lists(abundance, n_len) + site_labels_ = _get_default_lists(site_labels, n_len) + site_names_ = _get_default_lists(site_names, n_len) + site_descriptions_ = _get_default_lists(site_descriptions, n_len) + + index = np.where(abundance_ > rtol * abundance_.max())[0] + + sys = [ + SpinSystem( + sites=[ + Site( + isotope=ist__, + isotropic_chemical_shift=iso__, + label=lab__, + name=name__, + description=dis__, + ) + ], + abundance=abd__, + ) + for ist__, iso__, lab__, name__, dis__, abd__ in zip( + isotopes_[index], + iso_chemical_shifts_[index], + site_labels_[index], + site_names_[index], + site_descriptions_[index], + abundance_[index], ) - - keep_idxs = np.where(abundance > rtol * abundance.max())[0] - - return [ - SpinSystem(sites=[site], abundance=abd) - for site, abd in zip(sites[keep_idxs], abundance[keep_idxs]) ] - -def site_generator( - isotope: Union[str, List[str]], - isotropic_chemical_shift: Union[float, List[float], np.ndarray] = 0, - shielding_symmetric: Dict = None, - shielding_antisymmetric: Dict = None, - quadrupolar: Dict = None, - name: Union[str, List[str]] = None, - label: Union[str, List[str]] = None, - description: Union[str, List[str]] = None, -) -> List[Site]: - r"""Generate a list of Site objects from lists of site attributes. - - Args: - isotope: - A required string or a list of site isotopes. - isotropic_chemical_shift: - A float or a list/ndarray of isotropic chemical shifts per site. The default - is 0. - shielding_symmetric: - A shielding symmetric dict object, where the keyword value can either - be a float or a list/ndarray of floats. The default value is None. The - allowed keywords are ``zeta``, ``eta``, ``alpha``, ``beta``, and ``gamma``. - shielding_antisymmetric: - A shielding antisymmetric dict object, where the keyword value can either - be a float or a list/ndarray of floats. The default value is None. The - allowed keywords are ``zeta``, ``alpha``, and ``beta``. - quadrupolar: - A quadrupolar dict object, where the keyword value can either be a float or - a list/ndarray of floats. The default value is None. The allowed keywords - are ``Cq``, ``eta``, ``alpha``, ``beta``, and ``gamma``. - name: - A string or a list of strings with site names per site. The default is None. - label: - A string or a list of strings with site labels per site. The default is - None. - description: - A string or a list of strings with site descriptions per site. The default - is None. - - Returns: - sites: List of :ref:`site_api` objects - - Example: - **Generating 10 hydrogen sites:** - - >>> sites1 = site_generator( - ... isotope=["1H"] * 10, - ... isotropic_chemical_shift=-15, - ... name="10 Protons", - ... ) - >>> print(len(sites1)) - 10 - - **Generating 10 hydrogen sites with different shifts:** - - >>> shifts = np.arange(-10, 10, 2) - >>> sites2 = site_generator( - ... isotope=["1H"] * 10, - ... isotropic_chemical_shift=shifts, - ... name="10 Proton", - ... ) - >>> print(len(sites2)) - 10 - - **Generating multiple sites with dictionary arguments:** - - >>> Cq = [4.2e6] * 12 - >>> sys3 = site_generator( - ... isotope="17O", - ... isotropic_chemical_shift=60.0, # in ppm, - ... quadrupolar={"Cq": Cq, "eta": 0.5}, # Cq in Hz - ... ) - >>> print(len(sys3)) - 12 - """ - lst = [isotope, isotropic_chemical_shift, name, label, description] - attributes = [_fix_item(item) for item in lst] - - n_sites = _check_lengths(attributes) - - lst_extend = [shielding_symmetric, shielding_antisymmetric, quadrupolar] - for obj in lst_extend: - if obj is not None: - obj_ext, n_dict = _extend_dict_values(obj, n_sites) - n_sites = max(n_sites, n_dict) - attributes.append(obj_ext) - else: - attributes.append(None) - - # Attributes order is same as below in list comprehension - attributes = [_extend_to_nparray(attr, n_sites) for attr in attributes] - - return np.asarray( - [ - Site( - isotope=iso, - isotropic_chemical_shift=shift, - name=name, - label=label, - description=desc, - shielding_symmetric=symm, - shielding_antisymmetric=antisymm, - quadrupolar=quad, - ) - for iso, shift, name, label, desc, symm, antisymm, quad in zip(*attributes) + # system with additional shielding symmetric parameters + if shielding_symmetric is not None: + lst = [ + _get_default_lists(shielding_symmetric[item], n_len)[index] + if item in shield_keys + else [None for _ in range(index.size)] + for item in SHIELDING_SYM_PARAMS ] - ) - - -def _fix_item(item): - """Flattens multidimensional arrays into 1d array.""" - if isinstance(item, (list, np.ndarray)): - return np.asarray(item).ravel() - return item + _populate_shielding(sys, lst) + + # system with additional quadrupolar parameters + if quadrupolar is not None: + lst = [ + _get_default_lists(quadrupolar[item], n_len)[index] + if item in quad_keys + else [None for _ in range(index.size)] + for item in QUADRUPOLAR_PARAMS + ] + _populate_quadrupolar(sys, lst) + return sys -def _extend_to_nparray(item, n): - """If item is already list/array return np.array, otherwise extend to length n.""" - data = item if isinstance(item, (list, np.ndarray)) else [item for _ in range(n)] - return np.asarray(data) +def _get_shielding_info(shielding_symmetric): + n_ss, shield_keys = [], [] + if shielding_symmetric is not None: + shield_keys = shielding_symmetric.keys() + shielding_symmetric = { + item: _fix_item(shielding_symmetric[item]) + for item in SHIELDING_SYM_PARAMS + if item in shield_keys + } + n_ss = [ + _get_length(shielding_symmetric[item]) + for item in SHIELDING_SYM_PARAMS + if item in shield_keys + ] + return n_ss, shield_keys, shielding_symmetric -def _extend_dict_values(_dict, n_sites): - """Checks and extends dict values. Returns dict or list of dicts and max length.""" - _dict = {key: _fix_item(val) for key, val in _dict.items()} - n_sites_dict = _check_lengths(list(_dict.values())) - if n_sites != 1 and n_sites_dict != 1 and n_sites != n_sites_dict: - raise ValueError(f"A list in a dictionary was misshapen. {LIST_LEN_ERROR_MSG}") - if n_sites_dict == 1: - _dict = { - key: val[0] if isinstance(val, (list, np.ndarray)) else val - for key, val in _dict.items() +def _get_quad_info(quadrupolar): + n_q, quad_keys = [], [] + if quadrupolar is not None: + quad_keys = quadrupolar.keys() + quadrupolar = { + item: _fix_item(quadrupolar[item]) + for item in QUADRUPOLAR_PARAMS + if item in quad_keys + } + n_q = [ + _get_length(quadrupolar[item]) + for item in QUADRUPOLAR_PARAMS + if item in quad_keys + ] + return n_q, quad_keys, quadrupolar + + +def _populate_quadrupolar(sys, items): + n = len(items[0]) + for i in range(n): + if sys[i].sites[0].isotope.spin > 0.5: + sys[i].sites[0].quadrupolar = { + "Cq": items[0][i], + "eta": items[1][i], + "alpha": items[2][i], + "beta": items[3][i], + "gamma": items[4][i], + } + + +def _populate_shielding(sys, items): + n = len(items[0]) + for i in range(n): + sys[i].sites[0].shielding_symmetric = { + "zeta": items[0][i], + "eta": items[1][i], + "alpha": items[2][i], + "beta": items[3][i], + "gamma": items[4][i], } - return _dict, 1 - _dict = {key: _extend_to_nparray(val, n_sites_dict) for key, val in _dict.items()} - return _zip_dict(_dict), n_sites_dict +def _get_default_lists(item, n): + if isinstance(item, (list, np.ndarray)): + return np.asarray(item) -def _check_lengths(attributes): - """Ensures all attribute lengths are 1 or maximum attribute length.""" - lengths = np.array([np.asarray(attr).size for attr in attributes]) + return np.asarray([item for _ in range(n)]) - if np.all(lengths == 1): - return 1 - lengths = lengths[np.where(lengths != 1)] - if np.unique(lengths).size == 1: - return lengths[0] +def _fix_item(item): + if isinstance(item, (list, np.ndarray)): + return np.asarray(item).ravel() + return item - raise ValueError( - f"An array or list was either too short or too long. {LIST_LEN_ERROR_MSG}" - ) +def _get_length(item): + """Return the length of item it item is a list, else 0.""" + if isinstance(item, (list, np.ndarray)): + return np.asarray(item).size + return 0 -# BUG: doctest fails on example code -def _zip_dict(_dict): - """Makes list of dicts with the same keys and scalar values from dict of lists. - Single dictionaries of only None will return None. - Example: - >>> foo = {'k1': [1, None, 3, 4], 'k2': [5, None, 7, 8], 'k3': [9, None, 11, 12]} - >>> pprint(_zip_dict(foo)) - [{'k1': 1, 'k2': 5, 'k3': 9}, - None, - {'k1': 3, 'k2': 7, 'k3': 11}, - {'k1': 4, 'k2': 8, 'k3': 12}] - """ - lst = [dict(zip(_dict.keys(), v)) for v in zip(*(_dict[k] for k in _dict.keys()))] - lst = [None if np.all([item is None for item in d.values()]) else d for d in lst] - return lst +def _check_size(n_list): + index = np.where(n_list > 0) + n_list_reduced = n_list[index] + first_item = n_list_reduced[0] + if np.all(n_list_reduced == first_item): + return first_item + raise ValueError( + "Each entry can either be a single item or a list of items. If an entry is a " + "list, it's length must be equal to the length of other lists present in the " + "system." + ) diff --git a/src/mrsimulator/utils/emcee.py b/src/mrsimulator/utils/emcee.py new file mode 100644 index 000000000..846d36b7c --- /dev/null +++ b/src/mrsimulator/utils/emcee.py @@ -0,0 +1,200 @@ +import os +import csdmpy as cp +from mrsimulator.utils import get_spectral_dimensions +import matplotlib.pyplot as plt +import pandas as pd +import numpy as np +import emcee +from mrsimulator import Simulator, SpinSystem, Site +from mrsimulator.methods import BlochDecaySpectrum +from mrsimulator import signal_processing as sp +from lmfit import Minimizer, Parameters, fit_report, minimize +import json +from mrsimulator.utils.spectral_fitting import get_correct_data_order +from multiprocessing import Pool +import copy + +class mrsim_emcee: + ''' + A utility class to sample the posterior distribution function (PDF) of NMR paramters with Markov Chain Monte Carlo (MCMC) method. + This is a class depends on the emcee (https://emcee.readthedocs.io/en/stable/#) package. + + Parameters + -------------------------- + params: LMFIT parameters + Initial guess of the parameters to start MCMC walkers. User can use mrsimulator.utils.spectral_fitting.make_LMFIT_params to generate LMFIT parameters for mrsimulator. + simulator: MRsimulator simulator + The simulator object constructed with NMR paramters. + processor: MRsimulator processor + The processor object storing NMR signal processing parameters. + sigma: float + The noise level of the spectrum baseline. + + ''' + def __init__(self,params, simulator, processor, sigma = None): + self.params = self._vary_only(params) + self.nvarys = len(self.params) + self.simulator = simulator + self.processor = processor + self.sigma = sigma + + @staticmethod + def _vary_only(params): + ''' + Filter the parameter list, leave only variables. + + Parameters + ------------------------------ + params: Parameters + LMFIT parameters obj that store the variables. + ''' + for k, v in list(params.items()): + if not v.vary: + params.pop(k) + return(params) + + def mcmc(self,steps = 1000, nwalkers = 100, burn = 100, thin = 10, progress = True): + ''' + Bayesian sampling of the posterior distribution. + + This method uses the emcee package 'https://emcee.readthedocs.io/en/stable/' to sample the posterior distribution of NMR parameters with Marcov Chain Monte Carlo (MCMC) method. + + Parameters + ----------------------- + steps: int + The total number of samples to draw from the posterior distribution for each walkers. + nwalkers: int + The numbers of walkers in the emnsemble. + burn: + Number of steps to discard at the begining of the sampling process. The first few steps when the chain is making its way towards the center of the PDF, + are called ‘burn in’ and are removed since they do not represent the final PDF. + thin: int + Only accept 1 in every `thin` samples. + progress: bool + If True, show a progress bar while running. + + ''' + # randstate + seed = 100 + rand = np.random.RandomState(seed) + # initial guess + var = np.zeros(self.nvarys) + for i, v in enumerate(self.params): + param = self.params[v] + var[i] = param.value + p0 = 1 + rand.randn(nwalkers, self.nvarys) * 1.e-4 + p0*=var + + #run mcmc sampler + with Pool() as pool: + self.sampler = emcee.EnsembleSampler(nwalkers,self.nvarys,self._log_probability, + args = (self.params, self.simulator, self.processor, self.sigma), + pool = pool) + + self.sampler.run_mcmc(p0, steps, progress=progress) + + result = [] + result['flat_chain'] = self.sampler.get_chain(thin=thin, discard=burn, flat = True) + result['log_prob'] = self.sampler.get_log_prob(thin=thin, discard = burn) + quantiles = np.percentile(result['flat_chain'],[15.87,50,84.13],axis=0) + result['params'] = copy.deepcopy(self.params) + + for i,k in enumerate(result['params']): + std_l, median, std_u = quantiles[:,i] + result['params'][k].value = median + result['params'][k].stderr = 0.5*(std_u-std_l) + + return(result) + + @staticmethod + def _log_prior(theta,bounds): + ''' + Calculate the log prior for the parameters given. + If the values in params are out side its bound, return -inf, else return 0. + ''' + if np.any(theta>bounds[0,:]) or np.any(theta