TEST: Basic tests for dynamic MRS fitting

fsl_mrs/tests/test_utils_dynamic_dyn_results.py

+'''FSL-MRS test script
+
+Test the core dynamic MRS class
+
+Copyright Will Clarke, University of Oxford, 2021'''
+
+import pytest
+import numpy as np
+import pandas as pd
+import fsl_mrs.utils.synthetic as syn
+from fsl_mrs.core import MRS, basis
+import fsl_mrs.utils.dynamic as dyn
+
+
+# Fixture returning two MRS objects (with basis) linked by a concentration scaling of 2x.
+@pytest.fixture
+def fixed_ratio_mrs():
+    FID_basis1 = syn.syntheticFID(chemicalshift=[1, ], amplitude=[1], noisecovariance=[[0]], damping=[3])
+    FID_basis2 = syn.syntheticFID(chemicalshift=[3, ], amplitude=[1], noisecovariance=[[0]], damping=[3])
+    FID_basis1[1]['fwhm'] = 3 * np.pi
+    FID_basis2[1]['fwhm'] = 3 * np.pi
+    b = basis.Basis(
+        np.stack((FID_basis1[0][0], FID_basis2[0][0]), axis=1),
+        ['Met1', 'Met2'],
+        [FID_basis1[1], FID_basis2[1]])
+
+    FID1 = syn.syntheticFID(chemicalshift=[1, 3], amplitude=[1, 1], noisecovariance=[[0.01]])
+    FID2 = syn.syntheticFID(chemicalshift=[1, 3], amplitude=[2, 2], noisecovariance=[[0.01]])
+
+    mrs1 = MRS(FID=FID1[0][0], header=FID1[1], basis=b)
+    mrs2 = MRS(FID=FID2[0][0], header=FID2[1], basis=b)
+
+    mrs1.check_FID(repair=True)
+    mrs1.check_Basis(repair=True)
+    mrs2.check_FID(repair=True)
+    mrs2.check_Basis(repair=True)
+
+    mrs_list = [mrs1, mrs2]
+    return mrs_list
+
+
+def test_dynRes(fixed_ratio_mrs):
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+    res = dyn_obj.fit()
+
+    res_obj = res['result']
+
+    import plotly
+    fig = res_obj.plot_spectra(init=True, fit_to_init=True)
+    assert isinstance(fig, plotly.graph_objs._figure.Figure)
+
+    import matplotlib.pyplot
+    fig = res_obj.plot_mapped(fit_to_init=True)
+    assert isinstance(fig, matplotlib.pyplot.Figure)
+
+    assert isinstance(res_obj.data_frame, pd.DataFrame)
+    assert isinstance(res_obj.x, np.ndarray)
+    assert res_obj.x.shape[0] == res_obj.data_frame.shape[1]
+
+    assert isinstance(res_obj.mapped_parameters, np.ndarray)
+    assert res_obj.mapped_parameters.shape == (1, len(mrs_list), len(dyn_obj.mapped_names))
+
+    assert isinstance(res_obj.mapped_parameters_init, np.ndarray)
+    assert res_obj.mapped_parameters_init.shape == (len(mrs_list), len(dyn_obj.mapped_names))
+
+    assert isinstance(res_obj.free_parameters_init, np.ndarray)
+    assert res_obj.free_parameters_init.shape == (len(dyn_obj.free_names),)
+
+    assert isinstance(res_obj.init_dataframe, pd.DataFrame)
+    assert res_obj.init_dataframe.shape == (1, len(dyn_obj.free_names),)
+
+    assert isinstance(res_obj.mapped_parameters_fitted_init, np.ndarray)
+    assert res_obj.mapped_parameters_fitted_init.shape == (len(mrs_list), len(dyn_obj.mapped_names))
+
+    assert res_obj.mapped_names == dyn_obj.mapped_names
+    assert res_obj.free_names == dyn_obj.free_names
+
+
+def test_dynRes_newton(fixed_ratio_mrs):
+    """Test newton optimiser specific components"""
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+    res = dyn_obj.fit()
+
+    res_obj = res['result']
+    assert isinstance(res_obj.cov, pd.DataFrame)
+    assert res_obj.cov.shape == (10, 10)
+
+    assert isinstance(res_obj.corr, pd.DataFrame)
+    assert res_obj.corr.shape == (10, 10)
+
+    assert isinstance(res_obj.std, pd.Series)
+    assert res_obj.std.shape == (10,)
+
+    assert np.allclose(res_obj.std, np.sqrt(np.diagonal(res_obj.cov)))
+
+
+def test_dynRes_mcmc(fixed_ratio_mrs):
+    """Test mcmc optimiser specific components"""
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+
+    res = dyn_obj.fit(method='MH', mh_jumps=20)
+
+    res_obj = res['result']
+    assert isinstance(res_obj.cov, pd.DataFrame)
+    assert res_obj.cov.shape == (10, 10)
+
+    assert isinstance(res_obj.corr, pd.DataFrame)
+    assert res_obj.corr.shape == (10, 10)
+
+    assert isinstance(res_obj.std, pd.Series)
+    assert res_obj.std.shape == (10,)
+
+    assert np.allclose(res_obj.std, np.sqrt(np.diagonal(res_obj.cov)))

fsl_mrs/tests/test_utils_dynamic_dynmrs.py

+'''FSL-MRS test script
+
+Test the core dynamic MRS class
+
+Copyright Will Clarke, University of Oxford, 2021'''
+
+import pytest
+import numpy as np
+import fsl_mrs.utils.synthetic as syn
+from fsl_mrs.core import MRS, basis
+import fsl_mrs.utils.dynamic as dyn
+
+
+@pytest.fixture
+def fixed_ratio_mrs():
+    FID_basis1 = syn.syntheticFID(chemicalshift=[1, ], amplitude=[1], noisecovariance=[[0]], damping=[3])
+    FID_basis2 = syn.syntheticFID(chemicalshift=[3, ], amplitude=[1], noisecovariance=[[0]], damping=[3])
+    FID_basis1[1]['fwhm'] = 3 * np.pi
+    FID_basis2[1]['fwhm'] = 3 * np.pi
+    b = basis.Basis(
+        np.stack((FID_basis1[0][0], FID_basis2[0][0]), axis=1),
+        ['Met1', 'Met2'],
+        [FID_basis1[1], FID_basis2[1]])
+
+    FID1 = syn.syntheticFID(chemicalshift=[1, 3], amplitude=[1, 1], noisecovariance=[[0.01]])
+    FID2 = syn.syntheticFID(chemicalshift=[1, 3], amplitude=[2, 2], noisecovariance=[[0.01]])
+
+    mrs1 = MRS(FID=FID1[0][0], header=FID1[1], basis=b)
+    mrs2 = MRS(FID=FID2[0][0], header=FID2[1], basis=b)
+
+    mrs1.check_FID(repair=True)
+    mrs1.check_Basis(repair=True)
+    mrs2.check_FID(repair=True)
+    mrs2.check_Basis(repair=True)
+
+    mrs_list = [mrs1, mrs2]
+    return mrs_list
+
+
+def test_fixtures(fixed_ratio_mrs):
+    assert fixed_ratio_mrs[0].names == ['Met1', 'Met2']
+    assert fixed_ratio_mrs[1].names == ['Met1', 'Met2']
+
+
+def test_dynMRS_setup(fixed_ratio_mrs):
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+
+    # Test properties
+    assert dyn_obj.metabolite_names == ['Met1', 'Met2']
+    assert dyn_obj.free_names == [
+        'conc_c_0_Met1',
+        'conc_c_g_Met1',
+        'conc_c_0_Met2',
+        'conc_c_g_Met2',
+        'gamma_0',
+        'eps_0',
+        'Phi_0_0',
+        'Phi_1_0',
+        'baseline_0',
+        'baseline_1']
+    assert dyn_obj.mapped_names == [
+        'conc_Met1',
+        'conc_Met2',
+        'gamma_00',
+        'eps_00',
+        'Phi_0_00',
+        'Phi_1_00',
+        'baseline_00',
+        'baseline_01']
+
+    assert isinstance(dyn_obj.vm, dyn.VariableMapping)
+
+
+def test_process_mrs_list(fixed_ratio_mrs):
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+
+    assert dyn_obj.mrs_list[0].scaling['FID'] == 1.0
+    assert dyn_obj.mrs_list[0].scaling['basis'] == 1.0
+
+    dyn_obj_scaled = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=True)
+
+    assert dyn_obj_scaled.mrs_list[0].scaling['FID'] > 1.0
+    assert dyn_obj_scaled.mrs_list[0].scaling['basis'] > 1.0
+    assert dyn_obj_scaled.mrs_list[0].scaling['FID'] == dyn_obj_scaled.mrs_list[1].scaling['FID']
+
+
+# Test Utility methods
+def test_get_constants(fixed_ratio_mrs):
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+
+    consts = dyn_obj._get_constants(
+        mrs_list[0],
+        (0.2, 4.2),
+        0,
+        [0, 0])
+
+    assert len(consts) == 8
+    assert np.allclose(consts[0], mrs_list[0].frequencyAxis)
+    assert np.allclose(consts[1], mrs_list[0].timeAxis)
+    assert np.allclose(consts[2], mrs_list[0].basis)
+    assert consts[3].shape == (2048, 2)
+    assert consts[4] == [0, 0]
+    assert consts[5] == 1
+    assert consts[6] == mrs_list[0].ppmlim_to_range((0.2, 4.2), True)[0]
+    assert consts[7] == mrs_list[0].ppmlim_to_range((0.2, 4.2), True)[1]
+
+
+def test_dynMRS_fit(fixed_ratio_mrs):
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+    res = dyn_obj.fit()
+
+    concs = res['result'].data_frame.filter(like='conc').to_numpy()
+    assert np.allclose(concs, [1, 1, 1, 1], atol=0.1)
+
+
+def test_dynMRS_fit_mcmc(fixed_ratio_mrs):
+    mrs_list = fixed_ratio_mrs
+
+    dyn_obj = dyn.dynMRS(
+        mrs_list,
+        [0, 1],
+        'fsl_mrs/tests/testdata/dynamic/simple_linear_model.py',
+        model='lorentzian',
+        baseline_order=0,
+        metab_groups=[0, 0],
+        rescale=False)
+    res = dyn_obj.fit(method='MH', mh_jumps=50)
+
+    concs = res['result'].data_frame.filter(like='conc').to_numpy()
+    assert np.allclose(concs, [1, 1, 1, 1], atol=0.1)

fsl_mrs/utils/dynamic/dyn_results.py

@@ -16,6 +16,14 @@ from fsl_mrs.utils.misc import calculate_lap_cov
 
 # Plotting functions:
 def subplot_shape(plots_needed):
+    """Calculates the number and shape of subplots needed for
+    given number of plots.
+
+    :param plots_needed: Number of plots needed
+    :type plots_needed: int
+    :return: Number of columns and number of rows
+    :rtype: tuple
+    """
     col = int(np.floor(np.sqrt(plots_needed)))
     row = int(np.floor(np.sqrt(plots_needed)))
     counter = 0
@@ -29,8 +37,21 @@ def subplot_shape(plots_needed):
 
 
 class dynRes:
+    """Base dynamic results class. Not intended to be created directly, but via child classes
+    dynRes_newton and dynRes_mcmc."""
 
     def __init__(self, samples, dyn, init):
+        """Initilisation of dynRes class object.
+
+        Typically called from __init__ of dynRes_newton and dynRes_mcmc
+
+        :param samples: Array of free parameters returned by optimiser, can be 2D in mcmc case.
+        :type samples: numpy.ndarray
+        :param dyn: Copy of dynMRS class object.
+        :type dyn: fsl_mrs.utils.dynamic.dynMRS
+        :param init: Results of the initilisation optimisation, containing 'resList' and 'x'.
+        :type init: dict
+        """
         self._data = pd.DataFrame(data=samples, columns=dyn.free_names)
         self._data.index.name = 'samples'
 
@@ -41,10 +62,12 @@ class dynRes:
 
     @property
     def data_frame(self):
+        """Return the pandas dataframe view of the results."""
         return self._data
 
     @property
     def x(self):
+        """Return the (mcmc: mean) results as a numpy array."""
         return self.data_frame.mean().to_numpy()
 
     @staticmethod
@@ -218,7 +241,7 @@ class dynRes:
             ax.set_title(paramname)
             handles, labels = ax.get_legend_handles_labels()
         fig.legend(handles, labels, loc='right')
-        plt.show()
+        return fig
 
     def plot_spectra(self, init=False, fit_to_init=False):
         """Plot individual spectra as fitted using the dynamic model
@@ -327,26 +350,66 @@ class dynRes:
 
 
 class dynRes_mcmc(dynRes):
+    """Results class for MCMC optimised dynamic fitting.
 
+    Derived from parent dynRes class.
+    """
     def __init__(self, samples, dyn, init):
+        """Initilise MCMC dynamic fitting results object.
+
+        Simply calls parent class init.
+
+        :param samples: Array of free parameters returned by optimiser, can be 2D in mcmc case.
+        :type samples: numpy.ndarray
+        :param dyn: Copy of dynMRS class object.
+        :type dyn: fsl_mrs.utils.dynamic.dynMRS
+        :param init: Results of the initilisation optimisation, containing 'resList' and 'x'.
+        :type init: dict
+        """
         super().__init__(samples, dyn, init)
 
     @property
     def cov(self):
+        """Returns the covariance matrix of free parameters
+
+        :return: Covariance matrix as a DataFrame
+        :rtype: pandas.DataFrame
+        """
         return self.data_frame.cov()
 
     @property
     def corr(self):
+        """Returns the correlation matrix of free parameters
+
+        :return: Covariance matrix as a DataFrame
+        :rtype: pandas.DataFrame
+        """
         return self.data_frame.corr()
 
     @property
     def std(self):
+        """Returns the standard deviations of the free parameters
+
+        :return: Std as data Series
+        :rtype: pandas.Series
+        """
         return self.data_frame.std()
 
 
 class dynRes_newton(dynRes):
 
     def __init__(self, samples, dyn, init):
+        """Initilise TNC optimised dynamic fitting results object.
+
+        Calculates the covariance, correlation and standard deviations using the Fisher information matrix.
+
+        :param samples: Array of free parameters returned by optimiser, can be 2D in mcmc case.
+        :type samples: numpy.ndarray
+        :param dyn: Copy of dynMRS class object.
+        :type dyn: fsl_mrs.utils.dynamic.dynMRS
+        :param init: Results of the initilisation optimisation, containing 'resList' and 'x'.
+        :type init: dict
+        """
         super().__init__(samples[np.newaxis, :], dyn, init)
 
         # Calculate covariance, correlation and uncertainty
@@ -361,12 +424,27 @@ class dynRes_newton(dynRes):
 
     @property
     def cov(self):
+        """Returns the covariance matrix of free parameters
+
+        :return: Covariance matrix as a DataFrame
+        :rtype: pandas.DataFrame
+        """
         return pd.DataFrame(self._cov, self.free_names, self.free_names)
 
     @property
     def corr(self):
+        """Returns the correlation matrix of free parameters
+
+        :return: Covariance matrix as a DataFrame
+        :rtype: pandas.DataFrame
+        """
         return pd.DataFrame(self._corr, self.free_names, self.free_names)
 
     @property
     def std(self):
+        """Returns the standard deviations of the free parameters
+
+        :return: Std as data Series
+        :rtype: pandas.Series
+        """
         return pd.Series(self._std, self.free_names)

fsl_mrs/utils/dynamic/dynmrs.py

@@ -74,10 +74,14 @@ class dynMRS(object):
 
         numBasis, numGroups = self.mrs_list[0].numBasis, max(metab_groups) + 1
         varNames, varSizes = models.FSLModel_vars(model, numBasis, numGroups, baseline_order)
-        self.vm = self._create_vm(model, config_file, varNames, varSizes)
+        self._vm = varmap.VariableMapping(
+            param_names=varNames,
+            param_sizes=varSizes,
+            time_variable=self.time_var,
+            config_file=config_file)
 
     def _process_mrs_list(self):
-        """Apply single scaling
+        """Apply single scaling to the mrs_list
         """
         scales = []
         for mrs in self.mrs_list:
@@ -117,6 +121,10 @@ class dynMRS(object):
                 full_mapped_names.extend([f'{nn}_{idx:02.0f}' for idx in range(ns)])
         return full_mapped_names
 
+    @property
+    def vm(self):
+        return self._vm
+
     def fit(self,
             method='Newton',
             mh_jumps=600,
@@ -140,7 +148,7 @@ class dynMRS(object):
             start_time = time.time()
 
         if init is None:
-            init = self.initialise(**self._fit_args, verbose=verbose)
+            init = self.initialise(verbose=verbose)
         x0 = self.vm.mapped_to_free(init['x'])
 
         # MCMC or Newton
@@ -229,14 +237,6 @@ class dynMRS(object):
         g = max(metab_groups) + 1
         return (freq, time, basis, base_poly, metab_groups, g, first, last)
 
-    def _create_vm(self, model, config_file, varNames, varSizes):
-        """Create Variable Mapping object"""
-        vm = varmap.VariableMapping(param_names=varNames,
-                                    param_sizes=varSizes,
-                                    time_variable=self.time_var,
-                                    config_file=config_file)
-        return vm
-
     def _prepare_data(self, ppmlim):
         """FID to Spec and slice for fitting"""
         first, last = self.mrs_list[0].ppmlim_to_range(ppmlim)
@@ -339,7 +339,7 @@ class dynMRS(object):
         n_over_2 = len(self.data[0]) / 2
         for time_index in range(len(self.vm.time_variable)):
             p = np.hstack(mapped[time_index, :])
-            pred = self.forward(p)
+            pred = self.forward[time_index](p)
             ll += np.log(np.linalg.norm(pred - self.data[time_index])) * n_over_2
         return ll
 

release_instructions.md

@@ -17,4 +17,6 @@ Instructions for releasing a new version of FSL-MRS
     ```git merge upstream/master```  
     ```git push origin --tags```
 11. Trigger manual package build using the [fsl/conda/fsl-mrs:pre-conda](https://git.fmrib.ox.ac.uk/fsl/conda/fsl-mrs:pre-conda) branch. Do this by updating the version number [here](https://git.fmrib.ox.ac.uk/fsl/conda/fsl-mrs/-/blob/mnt/pre-conda/meta.yaml#L6). This will trigger a pipeline from which the conda package can be downloaded (*Download build-noarch-conda-package:archive artifact*).
-12. Send to MW for manual upload to the FSL conda channel.
\ No newline at end of file
+12. Send to MW for manual upload to the FSL conda channel.
+
+For local development installation run: ```pip install --no-deps -e .```
\ No newline at end of file