Merge branch 'enh/dynamic_improved_init' into 'master'

ENH: Add manual or mean initialisation to the dynamic fitting initilisation.

See merge request !31

fsl_mrs/tests/test_utils_dynamic_dynmrs.py

@@ -171,3 +171,32 @@ def test_dynMRS_fit_mcmc(fixed_ratio_mrs):
 
     concs = res['result'].data_frame.filter(like='conc').mean(axis=0).to_numpy()
     assert np.allclose(concs, [1, 1, 1, 1], atol=0.1)
+
+
+def test_dynMRS_mean_fit_init(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)
+
+    # check mean fit gives same results as fit of the mean mrs lsit
+    from fsl_mrs.utils.preproc.combine import combine_FIDs
+    from fsl_mrs.utils import fitting
+    from copy import deepcopy
+    mean_fid = combine_FIDs([mrs.FID for mrs in dyn_obj.mrs_list], 'mean')
+    mean_mrs = deepcopy(dyn_obj.mrs_list[0])
+    mean_mrs.FID = mean_fid
+
+    meanres = fitting.fit_FSLModel(mean_mrs, method='Newton', **dyn_obj._fit_args)
+    assert np.allclose(meanres.params, dyn_obj.fit_mean_spectrum())
+
+    # Check the init produces the same results via both interfaces
+    init1 = dyn_obj.initialise(indiv_init=meanres.params)
+    init2 = dyn_obj.initialise(indiv_init='mean')
+    assert np.allclose(np.hstack(np.hstack(init1['x'])), np.hstack(np.hstack(init2['x'])))

fsl_mrs/utils/dynamic/dynmrs.py

@@ -190,17 +190,21 @@ class dynMRS(object):
 
         return {'result': results, 'resList': res_list, 'optimisation_sol': sol}
 
-    def initialise(self, verbose=False):
+    def fit_mean_spectrum(self):
+        from fsl_mrs.utils.preproc.combine import combine_FIDs
+        from copy import deepcopy
+
+        mean_fid = combine_FIDs([mrs.FID for mrs in self.mrs_list], 'mean')
+        mean_mrs = deepcopy(self.mrs_list[0])
+        mean_mrs.FID = mean_fid
+        return fitting.fit_FSLModel(mean_mrs, method='Newton', **self._fit_args).params
+
+    def initialise(self, indiv_init=None, verbose=False):
         """Initialise the dynamic fitting using seperate fits of each spectrum.
 
-        :param model: Spectral model 'lorentzian' or 'voigt', defaults to 'voigt'
-        :type model: str, optional
-        :param metab_groups: List of metabolite groupings, defaults to None
-        :type metab_groups: List of ints, optional
-        :param ppmlim: Ppm range over which to fit, defaults to (.2, 4.2)
-        :type ppmlim: tuple, optional
-        :param baseline_order: Polynomial baseline order, defaults to 2, -1 disables.
-        :type baseline_order: int, optional
+        :param indiv_init: Optional initilisation of individual fits.
+            Can be a numpy array of mapped parameters, 'mean' (fits the mean spectrum), or None (independent).
+            Defaults to None
         :param verbose: Print information during fitting, defaults to False
         :type verbose: bool, optional
         :return: Dict containing free parameters and individual FitRes objects
@@ -209,6 +213,9 @@ class dynMRS(object):
         if verbose:
             start_time = time.time()
 
+        if isinstance(indiv_init, str) and indiv_init == 'mean':
+            indiv_init = self.fit_mean_spectrum()
+
         varNames = models.FSLModel_vars(self._fit_args['model'])
         numMetabs = self.mrs_list[0].numBasis
         numGroups = max(self._fit_args['metab_groups']) + 1
@@ -222,7 +229,7 @@ class dynMRS(object):
         for t, mrs in enumerate(self.mrs_list):
             if verbose:
                 print(f'Initialising {t + 1}/{len(self.mrs_list)}', end='\r')
-            res = fitting.fit_FSLModel(mrs, method='Newton', **self._fit_args)
+            res = fitting.fit_FSLModel(mrs, method='Newton', x0=indiv_init, **self._fit_args)
             resList.append(res)
             params = x2p(res.params, numMetabs, numGroups)
             for i, p in enumerate(params):

fsl_mrs/utils/qc.py

@@ -122,7 +122,8 @@ def noiseSD(spectrum, noisemask):
 
 def idNoiseRegion(mrs, spectrum, startingNoiseThreshold=0.001):
     """ Identify noise region in given spectrum"""
-    normspec = np.real(spectrum) / np.max(np.real(spectrum))
+    with np.errstate(divide='ignore', invalid='ignore'):
+        normspec = np.real(spectrum) / np.max(np.real(spectrum))
     noiseThreshold = startingNoiseThreshold
     noiseRegion = np.abs(normspec) < noiseThreshold
     # print(np.sum(noiseRegion))