diff --git a/CHANGELOG b/CHANGELOG index bcae4a13c..2fc121e64 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -8,6 +8,10 @@ What's new - Support gamma angle averaging. +**Simulator** + +- New instance method for the `Simulator` class -- `.optimize()` -- which pre-computes transition pathways before least-squares minimization. This improves the efficiency of least-squares minimizations. + v0.7.0 ------ diff --git a/docs/introduction/fitting_example.rst b/docs/introduction/fitting_example.rst index 21f63436b..b2cf8440e 100644 --- a/docs/introduction/fitting_example.rst +++ b/docs/introduction/fitting_example.rst @@ -520,7 +520,15 @@ minimization, and print the :context: close-figs from lmfit import Minimizer - minner = Minimizer(sf.LMFIT_min_function, fit_parameters, fcn_args=(sim, processor, sigma)) + # Optimize fitting by pre-computing transition pathways + opt = sim.optimize() + + minner = Minimizer( + sf.LMFIT_min_function, + fit_parameters, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt} + ) result = minner.minimize() result diff --git a/docs/introduction/isotopomers_example.rst b/docs/introduction/isotopomers_example.rst index 788345361..780ca4247 100644 --- a/docs/introduction/isotopomers_example.rst +++ b/docs/introduction/isotopomers_example.rst @@ -215,7 +215,7 @@ spin systems and the list of your two methods, and run the simulations. sim = Simulator( spin_systems = [isotopomer1, isotopomer2, isotopomer3], - methods=[method_H, method_C] + methods = [method_H, method_C] ) sim.run() diff --git a/fitting_source/1D_fitting/plot_1_31P_Na2PO4_MAS.py b/fitting_source/1D_fitting/plot_1_31P_Na2PO4_MAS.py index c3290b439..c596b32b0 100644 --- a/fitting_source/1D_fitting/plot_1_31P_Na2PO4_MAS.py +++ b/fitting_source/1D_fitting/plot_1_31P_Na2PO4_MAS.py @@ -99,14 +99,6 @@ experiment=experiment, # experimental dataset ) -# A method object queries every spin system for a list of transition pathways that are -# relevant to the given method. Since the method and the number of spin systems remains -# unchanged during the least-squares analysis, a one-time query is sufficient. To avoid -# querying for the transition pathways at every iteration in a least-squares fitting, -# evaluate the transition pathways once and store it as follows -for sys in spin_systems: - sys.transition_pathways = MAS.get_transition_pathways(sys) - # %% # **Step 3:** Create the Simulator object and add the method and spin system objects. sim = Simulator(spin_systems=spin_systems, methods=[MAS]) @@ -169,13 +161,28 @@ print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"])) # %% -# **Step 7:** Perform the least-squares minimization. For the user's convenience, we +# **Step 7:** Perform the least-squares minimization. +# A method object queries every spin system for a list of transition pathways that are +# relevant to the given method. Since the method and the number of spin systems remains +# unchanged during the least-squares analysis, a one-time query is sufficient. To avoid +# querying for the transition pathways at every iteration in a least-squares fitting, +# call the :py:mth:~`mrsimulator.Simulator.optimize()` method to pre-compute the +# pathways. +# +# For the user's convenience, we # also provide a utility function, # :func:`~mrsimulator.utils.spectral_fitting.LMFIT_min_function`, for evaluating the # difference vector between the simulation and experiment, based on -# the parameters update. You may use this function directly as the argument of the -# LMFIT Minimizer class, as follows, -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +# the parameters update. You may use this function directly to instantiate the +# LMFIT Minimizer class where `fcn_args` and `fcn_kws` are arguments passed to the +# function, as follows, +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/1D_fitting/plot_1_31p_Na2PO4_static.py b/fitting_source/1D_fitting/plot_1_31p_Na2PO4_static.py index 50a6f4fa7..533f456d9 100644 --- a/fitting_source/1D_fitting/plot_1_31p_Na2PO4_static.py +++ b/fitting_source/1D_fitting/plot_1_31p_Na2PO4_static.py @@ -73,11 +73,6 @@ experiment=experiment, # experimental dataset ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = static1D.get_transition_pathways(sys) - # %% # **Guess Model Spectrum** @@ -122,7 +117,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result 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 86b0efa49..f18c67a97 100644 --- a/fitting_source/1D_fitting/plot_2_13C_glycine_960Hz.py +++ b/fitting_source/1D_fitting/plot_2_13C_glycine_960Hz.py @@ -76,12 +76,6 @@ spectral_dimensions=spectral_dims, experiment=experiment, # experimental dataset ) - -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = MAS.get_transition_pathways(sys) - # %% # **Guess Model Spectrum** @@ -126,7 +120,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result 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 index 0fbab9a0f..f0038c580 100644 --- 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 @@ -196,7 +196,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processors, sigmas)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processors, sigmas), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/1D_fitting/plot_2_PASS_cross_sections.py b/fitting_source/1D_fitting/plot_2_PASS_cross_sections.py index b15f35651..d4c3bb9e8 100644 --- a/fitting_source/1D_fitting/plot_2_PASS_cross_sections.py +++ b/fitting_source/1D_fitting/plot_2_PASS_cross_sections.py @@ -76,10 +76,6 @@ experiment=pass_cross_section, # also add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = PASS.get_transition_pathways(sys) # %% # **Guess Spectrum** @@ -122,7 +118,13 @@ # %% # Run the minimization using LMFIT -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/1D_fitting/plot_3_Na2SiO3.py b/fitting_source/1D_fitting/plot_3_Na2SiO3.py index 19794624c..1c3a37bbd 100644 --- a/fitting_source/1D_fitting/plot_3_Na2SiO3.py +++ b/fitting_source/1D_fitting/plot_3_Na2SiO3.py @@ -104,13 +104,6 @@ experiment=experiment, # experimental dataset ) -# A method object queries every spin system for a list of transition pathways that are -# relevant for the given method. Since the method and the number of spin systems remain -# the same during the least-squares fit, a one-time query is sufficient. To avoid -# querying for the transition pathways at every iteration in a least-squares fitting, -# evaluate the transition pathways once and store it as follows -for sys in spin_systems: - sys.transition_pathways = MAS_CT.get_transition_pathways(sys) # %% # **Step 3:** Create the Simulator object and add the method and spin system objects. @@ -177,9 +170,16 @@ # provide a utility function, # :func:`~mrsimulator.utils.spectral_fitting.LMFIT_min_function`, for evaluating the # difference vector between the simulation and experiment, based on -# the parameters update. You may use this function directly as the argument of the -# LMFIT Minimizer class, as follows, -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +# the parameters update. You may use this function directly to instantiate the +# LMFIT Minimizer class where `fcn_args` and `fcn_kws` are arguments passed to the +# function, as follows, +opt = sim.optimize() +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/1D_fitting/plot_4_11B_Quad_NMR.py b/fitting_source/1D_fitting/plot_4_11B_Quad_NMR.py index dfe4af7c4..2756c8b61 100644 --- a/fitting_source/1D_fitting/plot_4_11B_Quad_NMR.py +++ b/fitting_source/1D_fitting/plot_4_11B_Quad_NMR.py @@ -70,11 +70,6 @@ experiment=experiment, # add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = MAS_CT.get_transition_pathways(sys) - # %% # **Guess Model Spectrum** @@ -119,7 +114,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/1D_fitting/plot_4_27Al_YAG.py b/fitting_source/1D_fitting/plot_4_27Al_YAG.py index 58b62cc38..35d591f29 100644 --- a/fitting_source/1D_fitting/plot_4_27Al_YAG.py +++ b/fitting_source/1D_fitting/plot_4_27Al_YAG.py @@ -83,11 +83,6 @@ experiment=experiment, # add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = MAS.get_transition_pathways(sys) - # %% # **Guess Spectrum** @@ -132,7 +127,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/1D_fitting/plot_4_2H_quad.py b/fitting_source/1D_fitting/plot_4_2H_quad.py index cdc5519c1..7353f3de5 100644 --- a/fitting_source/1D_fitting/plot_4_2H_quad.py +++ b/fitting_source/1D_fitting/plot_4_2H_quad.py @@ -71,10 +71,6 @@ experiment=experiment, # experimental dataset ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = MAS.get_transition_pathways(sys) # %% # **Guess Model Spectrum** @@ -120,7 +116,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/1D_fitting/plot_5_119Sn_sideband.py b/fitting_source/1D_fitting/plot_5_119Sn_sideband.py index 516a6e0bc..63ddc0b4b 100644 --- a/fitting_source/1D_fitting/plot_5_119Sn_sideband.py +++ b/fitting_source/1D_fitting/plot_5_119Sn_sideband.py @@ -88,10 +88,6 @@ experiment=experiment, # add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = MAS.get_transition_pathways(sys) # %% # **Guess Spectrum** @@ -155,7 +151,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py b/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py index cb3c8d46b..54c9ae171 100644 --- a/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py +++ b/fitting_source/2D_fitting/plot_1_LHistidine_PASS.py @@ -89,11 +89,6 @@ experiment=mat_dataset, # add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = PASS.get_transition_pathways(sys) - # %% # **Guess Spectrum** @@ -137,7 +132,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/2D_fitting/plot_1_Rb2SO4_QMAT.py b/fitting_source/2D_fitting/plot_1_Rb2SO4_QMAT.py index a3e4ac521..5e1e1d680 100644 --- a/fitting_source/2D_fitting/plot_1_Rb2SO4_QMAT.py +++ b/fitting_source/2D_fitting/plot_1_Rb2SO4_QMAT.py @@ -85,11 +85,6 @@ experiment=qmat_dataset, # add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = PASS.get_transition_pathways(sys) - # %% # **Guess Spectrum** @@ -135,7 +130,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/2D_fitting/plot_2_Coesite_DAS.py b/fitting_source/2D_fitting/plot_2_Coesite_DAS.py index a01dcccd9..40e2dfa1b 100644 --- a/fitting_source/2D_fitting/plot_2_Coesite_DAS.py +++ b/fitting_source/2D_fitting/plot_2_Coesite_DAS.py @@ -115,11 +115,6 @@ experiment=experiment, # also add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the das method. -for sys in spin_systems: - sys.transition_pathways = DAS.get_transition_pathways(sys) - # %% # **Guess Spectrum** @@ -166,11 +161,16 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize(method="powell") result - # %% # The best fit solution # --------------------- diff --git a/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py b/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py index fb3fd6d52..2fafef84d 100644 --- a/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py +++ b/fitting_source/2D_fitting/plot_3_RbNO3_MQMAS.py @@ -84,11 +84,6 @@ experiment=experiment, # add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the das method. -for sys in spin_systems: - sys.transition_pathways = MQMAS.get_transition_pathways(sys) - # %% # **Guess Spectrum** @@ -135,7 +130,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/fitting_source/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.py b/fitting_source/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.py index e25e6d0e1..619ecf086 100644 --- a/fitting_source/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.py +++ b/fitting_source/2D_fitting/plot_4_NiCl2.2D2O_shifting-d.py @@ -123,11 +123,6 @@ experiment=experiment, # also add the measurement to the method. ) -# Optimize the script by pre-setting the transition pathways for each spin system from -# the method. -for sys in spin_systems: - sys.transition_pathways = shifting_d.get_transition_pathways(sys) - # %% # **Guess Spectrum** @@ -173,7 +168,13 @@ # %% # **Solve the minimizer using LMFIT** -minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +opt = sim.optimize() # Pre-compute transition pathways +minner = Minimizer( + sf.LMFIT_min_function, + params, + fcn_args=(sim, processor, sigma), + fcn_kws={"opt": opt}, +) result = minner.minimize() result diff --git a/src/c_lib/base/base_model.pyx b/src/c_lib/base/base_model.pyx index 574868c6c..9c138f10b 100644 --- a/src/c_lib/base/base_model.pyx +++ b/src/c_lib/base/base_model.pyx @@ -15,6 +15,8 @@ clib.generate_tables() @cython.wraparound(False) def core_simulator(method, list spin_systems, + list transition_pathways, # Same length as spin_systems list and + list transition_weights, # elements coorespond to each system int verbose=0, # for debug purpose only. unsigned int number_of_sidebands=90, unsigned int integration_density=72, @@ -165,9 +167,7 @@ def core_simulator(method, affine_matrix_c[3] -= affine_matrix_c[1]*affine_matrix_c[2] # sites _______________________________________________________________________________ - p_isotopes = None - - cdef int number_of_sites, number_of_couplings, p_number_of_sites=0 + cdef int number_of_sites, number_of_couplings cdef ndarray[int] spin_index_ij cdef ndarray[float] spin_i cdef ndarray[double] gyromagnetic_ratio_i @@ -206,7 +206,11 @@ def core_simulator(method, # ------------------------------------------------------------------------- # sample __________________________________________________________________ - for spin_sys in spin_systems: + for idx in range(len(spin_systems)): + spin_sys = spin_systems[idx] + segments = transition_pathways[idx] + weights = transition_weights[idx] + abundance = spin_sys.abundance isotopes = [site.isotope.symbol for site in spin_sys.sites] if channel not in isotopes: @@ -393,28 +397,11 @@ def core_simulator(method, # amp_individual.append([]) # continue - 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, dtype=np.float32) - transition_pathway_c = transition_pathway.ravel() - transition_pathway_weight_c = weights.view(dtype=np.float64) - else: - # convert transition objects to list - weights = [(item.weight.real, item.weight.imag) for item in transition_pathway] - transition_pathway_weight_c = np.asarray(weights, dtype=np.float64).ravel() - # transition_pathway_weight_c = weights - - transition_pathway = np.asarray(transition_pathway) - lst = [item.tolist() for item in transition_pathway.ravel()] - transition_pathway_c = np.asarray(lst, dtype=np.float32).ravel() - - pathway_count, transition_count_per_pathway = transition_pathway.shape[:2] - pathway_increment = 2*number_of_sites*transition_count_per_pathway + transition_pathway_c = np.asarray(segments, dtype=np.float32).ravel() # NOTE: Why 32 then 64 bits? + transition_pathway_weight_c = weights.view(dtype=np.float64) - p_number_of_sites = number_of_sites - p_isotopes = isotopes + pathway_count, transition_count_per_pathway = segments.shape[:2] + pathway_increment = 2*number_of_sites*transition_count_per_pathway # if spin_sys.transitions is not None: # transition_pathway_c = np.asarray( diff --git a/src/mrsimulator/method/__init__.py b/src/mrsimulator/method/__init__.py index f6e6d4900..0de4e30a0 100644 --- a/src/mrsimulator/method/__init__.py +++ b/src/mrsimulator/method/__init__.py @@ -479,12 +479,12 @@ def _get_transition_pathways_np(self, spin_system): 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 - ] + return np.array( + [[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 == []: + def _get_transition_pathway_weights_np(self, pathways, spin_system): + if np.asarray(pathways).size == 0: return np.asarray([]) symbol = [item.isotope.symbol for item in spin_system.sites] @@ -546,7 +546,7 @@ def _calculate_transition_connect_weight(trans1, trans2, spins, alpha, beta, gam 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) + weights = self._get_transition_pathway_weights_np(segments, spin_system) return segments, weights # indexes = np.where(weights != 0)[0] # return np.asarray(segments)[indexes], weights[indexes] diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py index 6375c6105..75558b0a2 100644 --- a/src/mrsimulator/simulator/__init__.py +++ b/src/mrsimulator/simulator/__init__.py @@ -132,7 +132,6 @@ class Simulator(Parseable): spin_systems: List[SpinSystem] = [] methods: List[Method] = [] config: ConfigSimulator = ConfigSimulator() - # indexes = [] class Config: validate_assignment = True @@ -320,11 +319,41 @@ def export_methods(self, filename: str): mth, outfile, ensure_ascii=False, sort_keys=False, allow_nan=False ) + def optimize(self) -> None: + """Pre-computes transition pathways and associated weights for each of the + :py:class:`~mrsimulator.method.Method` and + :py:class:`~mrsimulator.spin_system.SpinSystem` objects held by the simulator. + This increases the efficiency during least-squared minimization since pathways + are not re-computed during every iteration. + + Example + ------- + + >>> sim = Simulator() + >>> # Add spin systems and methods + >>> optimization = sim.optimize() + >>> sim.run(opt=optimization) + """ + opt = {"precomputed_pathways": [], "precomputed_weights": []} + for mth in self.methods: + pathways = [] + weights = [] + for sys in self.spin_systems: + p, w = mth._get_transition_pathway_and_weights_np(sys) + pathways.append(p) + weights.append(w) + + opt["precomputed_pathways"].append(pathways) + opt["precomputed_weights"].append(weights) + + return opt + def run( self, method_index: list = None, n_jobs: int = 1, pack_as_csdm: bool = True, + opt: dict = None, **kwargs, ): """Run the simulation and compute spectrum. @@ -337,11 +366,15 @@ def run( bool pack_as_csdm: If true, the simulation results are stored as a `CSDM `_ object, otherwise, as a `ndarray - `_ + `_ object. The simulations are stored as the value of the :attr:`~mrsimulator.Method.simulation` attribute of the corresponding method. + dict opt: An optional optimization dictionary storing pre-computed + transition pathways and transition weights for the given methods and + spin systems in the simulator. The default is None, that is, the + pathways and weights are calculated in the run method. Example ------- @@ -349,54 +382,50 @@ def run( >>> sim.run() # doctest:+SKIP """ verbose = 0 - if method_index is None: + + if opt is None: + opt = self.optimize() + + if method_index is None: # Simulate for all methods method_index = np.arange(len(self.methods)) - elif isinstance(method_index, int): + elif isinstance(method_index, int): # Package single method index as list method_index = [method_index] + for index in method_index: method = self.methods[index] spin_sys = get_chunks(self.spin_systems, n_jobs) - kwargs_dict = self.config.get_int_dict() + + # Chunking transition pathways and weights form the optimization dictionary + pathways = get_chunks(opt["precomputed_pathways"][index], n_jobs) + weights = get_chunks(opt["precomputed_weights"][index], n_jobs) + + config_dict = self.config.get_int_dict() jobs = ( delayed(core_simulator)( - method=method, spin_systems=sys, **kwargs_dict, **kwargs + method=method, + spin_systems=sys, + transition_pathways=pth, + transition_weights=wht, + **config_dict, + **kwargs, ) - for sys in spin_sys + for sys, pth, wht in zip(spin_sys, pathways, weights) ) amp = Parallel( n_jobs=n_jobs, verbose=verbose, backend="loky", - # **{ - # "backend": { - # "threads": "threading", - # "processes": "multithreading", - # None: None, - # }["threads"] - # }, )(jobs) - # self.indexes.append(indexes) - gyromagnetic_ratio = method.channels[0].gyromagnetic_ratio B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density larmor_freq = np.abs(B0 * gyromagnetic_ratio * 1e6) for seq in method.spectral_dimensions: seq.origin_offset = larmor_freq + seq.reference_offset - if isinstance(amp[0], list): - simulated_dataset = [] - for item in amp: - simulated_dataset += item - if isinstance(amp[0], np.ndarray): - simulated_dataset = [np.asarray(amp).sum(axis=0)] - - method.simulation = ( - self._as_csdm_object(simulated_dataset, method) - if pack_as_csdm - else np.asarray(simulated_dataset) + self._package_amp_after_simulation( + method=method, amp=amp, pack_as_csdm=pack_as_csdm ) - amp = None def save(self, filename: str, with_units: bool = True): """Serialize the simulator object to a JSON file. @@ -532,6 +561,31 @@ def _get_dv_metadata(self, index): } return obj + def _package_amp_after_simulation(self, method, amp, pack_as_csdm): + """Helper function for packaging frequency array (amp) into a numpy array or + CSDM object after a simulation. Puts the packaged simulation into the passed + method object and does not return anything + + Args: + Method method: The :py:class:`~mrsimulator.Method` object for the simulated + spectrum. + np.ndarray amp: Calculated spectrum as a numpy array + bool pack_as_csdm: Packages the simulated spectrum as a CSDM object if true. + Otherwise kept as a numpy array. + """ + if isinstance(amp[0], list): + simulated_dataset = [] + for item in amp: + simulated_dataset += item + if isinstance(amp[0], np.ndarray): + simulated_dataset = [np.asarray(amp).sum(axis=0)] + + method.simulation = ( + self._as_csdm_object(simulated_dataset, method) + if pack_as_csdm + else np.asarray(simulated_dataset) + ) + class Sites(AbstractList): """A list of unique :ref:`site_api` objects within a simulator object.""" diff --git a/src/mrsimulator/simulator/tests/test_simulator.py b/src/mrsimulator/simulator/tests/test_simulator.py index 1721ad49f..290649424 100644 --- a/src/mrsimulator/simulator/tests/test_simulator.py +++ b/src/mrsimulator/simulator/tests/test_simulator.py @@ -392,3 +392,13 @@ def get_blocks(indexes): lst[i] += 1 items_list = np.arange(85).tolist() check_chunks(items_list, -1, lst) + + +def test_sim_optimize(): + sim = Simulator() + sim.spin_systems = [SpinSystem(sites=[Site(isotope="1H"), Site(isotope="23Na")])] + sim.methods = [BlochDecaySpectrum(channels=["1H"])] + + opt = sim.optimize() + + assert set(opt.keys()) == {"precomputed_pathways", "precomputed_weights"} diff --git a/src/mrsimulator/utils/spectral_fitting.py b/src/mrsimulator/utils/spectral_fitting.py index 416a7e7fd..9cac3e1a8 100644 --- a/src/mrsimulator/utils/spectral_fitting.py +++ b/src/mrsimulator/utils/spectral_fitting.py @@ -430,7 +430,11 @@ def _check_for_experiment_data(methods_list: list): def LMFIT_min_function( - params: Parameters, sim: Simulator, processors: list = None, sigma: list = None + params: Parameters, + sim: Simulator, + processors: list = None, + sigma: list = None, + **sim_kwargs, ): """The simulation routine to calculate the vector difference between simulation and experiment based on the parameters update. @@ -438,10 +442,12 @@ def LMFIT_min_function( Args: params: Parameters object containing parameters for OLS minimization. sim: Simulator object. - processors: A list of PostSimulator objects corresponding to the methods in the - Simulator object. + processors: A list of :py:class:~`mrsimulator.signal_processor.Processor` + objects to apply post-simulation processing to the simulated spectra. sigma: A list of standard deviations corresponding to the experiments in the - Simulator.methods attribute + :py:attr:~`mrsimulator.Simulator.methods` attribute. + sim_kwargs: Keyword arguments to pass to the + :py:mth:~`mrsimulator.Simulator.run()` method. Returns: Array of the differences between the simulation and the experimental datasets. """ @@ -452,7 +458,7 @@ def LMFIT_min_function( _check_for_experiment_data(sim.methods) update_mrsim_obj_from_params(params, sim, processors) - sim.run() + sim.run(**sim_kwargs) processed_dataset = [ item.apply_operations(dataset=data.simulation)