Merge branch 'enh/dynamic_init_options' into 'master'

ENH: Dynamic init options

See merge request fsl/fsl_mrs!70

CHANGELOG.rst

 This document contains the FSL-MRS release history in reverse chronological order.
 
-2.0.8 (WIP)
------------
+2.0.8 (Monday 7th November)
+---------------------------
 - Added check (and printed suggestion) in `fsl_mrs` and `fsl_mrsi` that default MM are added to appropriate metabolite groups.
 - Added ppm range option to `fsl_mrs_preproc` alignment stage. Use option `--align_limits`.
-
+- Added new initialisation options to dynamic fitting based on free (rather than mapped) parameters.
+- Truncation step in `fsl_mrs_preproc` now happens earlier in series.
 
 2.0.7 (Thursday 20th October 2022)
 -----------------------------------

fsl_mrs/dynamic/dyn_results.py

@@ -396,11 +396,18 @@ class dynRes:
         return df_dict
 
     # Plotting
-    def plot_mapped(self, tvals=None, fit_to_init=False):
+    def plot_mapped(self, tvals=None, fit_to_init=False, ground_truth=None):
         """Plot each mapped parameter across time points
 
+        :param tvals: 'time' values on x axis, defaults to None / those stored in results object
+        :type tvals: list, optional
         :param fit_to_init: Plot the mapped parameters as per initilisation, defaults to False
         :type fit_to_init: bool, optional
+        :param ground_truth: If a ground truth exists (from simulation) plot the mapped parameters
+            as calculated from this vector, defaults to None
+        :type ground_truth: numpy.array, optional
+        :return: Figure object
+        :rtype: matplotlib.pyplot.figure.Figure
         """
 
         init_params = self.init_mapped_parameters_array
@@ -410,17 +417,30 @@ class dynRes:
         names = self.mapped_names
         if tvals is None:
             tvals = self._time_index_1d()
+        if ground_truth is not None:
+            gtval = self._dyn.vm.free_to_mapped(ground_truth)
+        else:
+            gtval = np.empty(len(names), dtype=object)
 
         # Plot the lot
         row, col = subplot_shape(len(names))
 
         fig, axes = plt.subplots(row, col, figsize=(20, 20))
-        for ax, p_init, p_fit_init, p_dyn, p_dyn_sd, paramname \
-                in zip(axes.flatten(), init_params.T, fitted_init_params.T, dyn_params.T, dyn_params_sd.T, names):
+        for ax, p_init, p_fit_init, p_dyn, p_dyn_sd, gt, paramname \
+                in zip(
+                    axes.flatten(),
+                    init_params.T,
+                    fitted_init_params.T,
+                    dyn_params.T,
+                    dyn_params_sd.T,
+                    gtval.T,
+                    names):
             ax.plot(tvals, p_init, 'o', label='init')
+            ax.errorbar(tvals, p_dyn, yerr=p_dyn_sd, fmt='-', label='dyn')
             if fit_to_init:
                 ax.plot(tvals, p_fit_init, ':', label='fit to init')
-            ax.errorbar(tvals, p_dyn, yerr=p_dyn_sd, fmt='-', label='dyn')
+            if gt is not None:
+                ax.plot(tvals, gt, 'k--', label='Ground Truth')
             ax.set_title(paramname)
             handles, labels = ax.get_legend_handles_labels()
         fig.legend(handles, labels, loc='right')

fsl_mrs/dynamic/dynmrs.py

@@ -230,6 +230,7 @@ class dynMRS(object):
             method='Newton',
             mh_jumps=600,
             init=None,
+            x0=None,
             verbose=False,
             output_opt_sol=False):
         """Fit the dynamic model
@@ -238,8 +239,10 @@ class dynMRS(object):
         :type method: str, optional
         :param mh_jumps: Number of MH jumps, defaults to 600
         :type mh_jumps: int, optional
-        :param init: Initilisation (x0), defaults to None
+        :param init: Initilisation based on independent fitting approach, defaults to None
         :type init: dict, optional
+        :param x0: Initilisation based on free parameters, defaults to None
+        :type x0: np.array, optional
         :param verbose: Verbosity flag, defaults to False
         :param output_opt_sol: Output the scipy solution object (for debugging), defaults to False
         :type output_opt_sol: bool, optional
@@ -251,9 +254,14 @@ class dynMRS(object):
             print('Start fitting')
             start_time = time.time()
 
-        if init is None:
+        if init is None \
+                and x0 is None:
             init = self.initialise(verbose=verbose)
-        x0 = self.vm.mapped_to_free(init['x'])
+
+        if x0 is None:
+            x0 = self.vm.mapped_to_free(init['x'])
+        else:
+            init = {'x': self.vm.free_to_mapped(x0)}
 
         # MCMC or Newton
         if method.lower() == 'newton':

fsl_mrs/scripts/fsl_mrs_preproc.py

@@ -277,6 +277,13 @@ def main():
     if args.quant is not None:
         quant_data = nifti_mrs_proc.ecc(quant_data, quant_data)
 
+    if args.leftshift:
+        verbose_print('... Truncation ...')
+        supp_data = nifti_mrs_proc.truncate_or_pad(supp_data, -args.leftshift, 'first', report=report_dir)
+        ref_data = nifti_mrs_proc.truncate_or_pad(ref_data, -args.leftshift, 'first')
+        if args.quant is not None:
+            quant_data = nifti_mrs_proc.truncate_or_pad(quant_data, -args.leftshift, 'first')
+
     # HLSVD
     if args.hlsvd:
         if not args.average:
@@ -286,13 +293,6 @@ def main():
         hlsvdlimits = [-0.25, 0.25]
         supp_data = nifti_mrs_proc.remove_peaks(supp_data, hlsvdlimits, limit_units='ppm', report=report_dir)
 
-    if args.leftshift:
-        verbose_print('... Truncation ...')
-        supp_data = nifti_mrs_proc.truncate_or_pad(supp_data, -args.leftshift, 'first', report=report_dir)
-        ref_data = nifti_mrs_proc.truncate_or_pad(ref_data, -args.leftshift, 'first')
-        if args.quant is not None:
-            quant_data = nifti_mrs_proc.truncate_or_pad(quant_data, -args.leftshift, 'first')
-
     # Apply shift to reference
     verbose_print('... Shifting tCr to 3.027 ...')
     supp_data = nifti_mrs_proc.shift_to_reference(