Commit 2c7a2b83 authored by William Clarke's avatar William Clarke
Browse files

Merge branch 'dynmrs_tests' into 'master'

TEST: Basic tests for dynamic MRS fitting

1. Tests for the dynMRS and dynRes classes.
2. Fixed bug in MCMC optimisation introduced in editing (multi basis set) enhancement.
3. Removed unnecessary methods in dynMRS.

See merge request fsl/fsl_mrs!29
parents 6c9dd697 04c74d34
Pipeline #11286 failed with stages
in 13 minutes and 3 seconds
'''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 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)
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
...@@ -16,6 +16,14 @@ from fsl_mrs.utils.misc import calculate_lap_cov ...@@ -16,6 +16,14 @@ from fsl_mrs.utils.misc import calculate_lap_cov
# Plotting functions: # Plotting functions:
def subplot_shape(plots_needed): 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))) col = int(np.floor(np.sqrt(plots_needed)))
row = int(np.floor(np.sqrt(plots_needed))) row = int(np.floor(np.sqrt(plots_needed)))
counter = 0 counter = 0
...@@ -29,8 +37,21 @@ def subplot_shape(plots_needed): ...@@ -29,8 +37,21 @@ def subplot_shape(plots_needed):
class dynRes: 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): 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 = pd.DataFrame(data=samples, columns=dyn.free_names)
self._data.index.name = 'samples' self._data.index.name = 'samples'
...@@ -41,10 +62,12 @@ class dynRes: ...@@ -41,10 +62,12 @@ class dynRes:
@property @property
def data_frame(self): def data_frame(self):
"""Return the pandas dataframe view of the results."""
return self._data return self._data
@property @property
def x(self): def x(self):
"""Return the (mcmc: mean) results as a numpy array."""
return self.data_frame.mean().to_numpy() return self.data_frame.mean().to_numpy()
@staticmethod @staticmethod
...@@ -218,7 +241,7 @@ class dynRes: ...@@ -218,7 +241,7 @@ class dynRes:
ax.set_title(paramname) ax.set_title(paramname)
handles, labels = ax.get_legend_handles_labels() handles, labels = ax.get_legend_handles_labels()
fig.legend(handles, labels, loc='right') fig.legend(handles, labels, loc='right')
plt.show() return fig
def plot_spectra(self, init=False, fit_to_init=False): def plot_spectra(self, init=False, fit_to_init=False):
"""Plot individual spectra as fitted using the dynamic model """Plot individual spectra as fitted using the dynamic model
...@@ -327,26 +350,66 @@ class dynRes: ...@@ -327,26 +350,66 @@ class dynRes:
class dynRes_mcmc(dynRes): class dynRes_mcmc(dynRes):
"""Results class for MCMC optimised dynamic fitting.
Derived from parent dynRes class.
"""
def __init__(self, samples, dyn, init): 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) super().__init__(samples, dyn, init)
@property @property
def cov(self): def cov(self):
"""Returns the covariance matrix of free parameters
:return: Covariance matrix as a DataFrame
:rtype: pandas.DataFrame
"""
return self.data_frame.cov() return self.data_frame.cov()
@property @property
def corr(self): def corr(self):
"""Returns the correlation matrix of free parameters
:return: Covariance matrix as a DataFrame
:rtype: pandas.DataFrame
"""
return self.data_frame.corr() return self.data_frame.corr()
@property @property
def std(self): def std(self):
"""Returns the standard deviations of the free parameters
:return: Std as data Series
:rtype: pandas.Series
"""
return self.data_frame.std() return self.data_frame.std()
class dynRes_newton(dynRes): class dynRes_newton(dynRes):
def __init__(self, samples, dyn, init): 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) super().__init__(samples[np.newaxis, :], dyn, init)
# Calculate covariance, correlation and uncertainty # Calculate covariance, correlation and uncertainty
...@@ -361,12 +424,27 @@ class dynRes_newton(dynRes): ...@@ -361,12 +424,27 @@ class dynRes_newton(dynRes):
@property @property
def cov(self): 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) return pd.DataFrame(self._cov, self.free_names, self.free_names)
@property @property
def corr(self): 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) return pd.DataFrame(self._corr, self.free_names, self.free_names)
@property @property
def std(self): 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) return pd.Series(self._std, self.free_names)
...@@ -74,10 +74,14 @@ class dynMRS(object): ...@@ -74,10 +74,14 @@ class dynMRS(object):
numBasis, numGroups = self.mrs_list[0].numBasis, max(metab_groups) + 1 numBasis, numGroups = self.mrs_list[0].numBasis, max(metab_groups) + 1
varNames, varSizes = models.FSLModel_vars(model, numBasis, numGroups, baseline_order) 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): def _process_mrs_list(self):
"""Apply single scaling """Apply single scaling to the mrs_list
""" """
scales = [] scales = []
for mrs in self.mrs_list: for mrs in self.mrs_list:
...@@ -117,6 +121,10 @@ class dynMRS(object): ...@@ -117,6 +121,10 @@ class dynMRS(object):
full_mapped_names.extend([f'{nn}_{idx:02.0f}' for idx in range(ns)]) full_mapped_names.extend([f'{nn}_{idx:02.0f}' for idx in range(ns)])
return full_mapped_names return full_mapped_names
@property
def vm(self):
return self._vm
def fit(self, def fit(self,
method='Newton', method='Newton',
mh_jumps=600, mh_jumps=600,
...@@ -140,7 +148,7 @@ class dynMRS(object): ...@@ -140,7 +148,7 @@ class dynMRS(object):
start_time = time.time() start_time = time.time()
if init is None: 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']) x0 = self.vm.mapped_to_free(init['x'])
# MCMC or Newton # MCMC or Newton
...@@ -229,14 +237,6 @@ class dynMRS(object): ...@@ -229,14 +237,6 @@ class dynMRS(object):
g = max(metab_groups) + 1 g = max(metab_groups) + 1
return (freq, time, basis, base_poly, metab_groups, g, first, last) 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): def _prepare_data(self, ppmlim):
"""FID to Spec and slice for fitting""" """FID to Spec and slice for fitting"""
first, last = self.mrs_list[0].ppmlim_to_range(ppmlim) first, last = self.mrs_list[0].ppmlim_to_range(ppmlim)
...@@ -339,7 +339,7 @@ class dynMRS(object): ...@@ -339,7 +339,7 @@ class dynMRS(object):
n_over_2 = len(self.data[0]) / 2 n_over_2 = len(self.data[0]) / 2
for time_index in range(len(self.vm.time_variable)): for time_index in range(len(self.vm.time_variable)):
p = np.hstack(mapped[time_index, :])