diff --git a/fsl_mrs/tests/test_utils_dynamic_dyn_results.py b/fsl_mrs/tests/test_utils_dynamic_dyn_results.py new file mode 100644 index 0000000000000000000000000000000000000000..424105c5740e8f27e40dde5b745caebce5b9a805 --- /dev/null +++ b/fsl_mrs/tests/test_utils_dynamic_dyn_results.py @@ -0,0 +1,140 @@ +'''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))) diff --git a/fsl_mrs/tests/test_utils_dynamic_dynmrs.py b/fsl_mrs/tests/test_utils_dynamic_dynmrs.py new file mode 100644 index 0000000000000000000000000000000000000000..bae4442b94290237eefa7a79decb84c76e7ae761 --- /dev/null +++ b/fsl_mrs/tests/test_utils_dynamic_dynmrs.py @@ -0,0 +1,173 @@ +'''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) diff --git a/fsl_mrs/tests/testdata/dynamic/simple_linear_model.py b/fsl_mrs/tests/testdata/dynamic/simple_linear_model.py new file mode 100644 index 0000000000000000000000000000000000000000..2f32f21dc564c6b1602ccb87b38cb4802247647f --- /dev/null +++ b/fsl_mrs/tests/testdata/dynamic/simple_linear_model.py @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:97bfd4a0ba3bbab4f000e7edea26449a358e55b2e2c170c718855984a1b56c0a +size 863 diff --git a/fsl_mrs/utils/dynamic/dyn_results.py b/fsl_mrs/utils/dynamic/dyn_results.py index 2b005d22f4cd47605054424fcc25204f7304d969..005eb0b5a464edf685c629a51f01db1a0843789a 100644 --- a/fsl_mrs/utils/dynamic/dyn_results.py +++ b/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) diff --git a/fsl_mrs/utils/dynamic/dynmrs.py b/fsl_mrs/utils/dynamic/dynmrs.py index 901fc3b1ed2ee6d90d46b488018f9250944f847f..26f136adf0858db602453eb835b9790040fafac2 100644 --- a/fsl_mrs/utils/dynamic/dynmrs.py +++ b/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 diff --git a/release_instructions.md b/release_instructions.md index b16cf7c699474cd651b88a08bc0c6cb20f780902..97dbc65b46afacd95057648100b6d3e93955098d 100644 --- a/release_instructions.md +++ b/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