Added the tested code for synth_and_ana

synth_and_ana.py

@@ -0,0 +1,99 @@
+#  Synthetise new data and analyse it directly to return fit parametersi
+
+import numpy as np
+import panda as pd
+
+from fsl_mrs.core import MRS
+from fsl_mrs.utils import plotting, mrs_io
+from fsl_mrs.utils.synthetic import synthetic_from_basis as synth
+from fsl_mrs.utils.misc import parse_metab_groups
+from fsl_mrs.utils.fitting import fit_FSLModel
+
+
+def synth_and_ana(noise_cov,
+                  fit_parameters, fit_snr, fit_sdnoise, fit_varnoise,
+                  df_parameter_synth, basis_path, output_path,
+                  verbose=False):
+    """Synthetise spectra with given noise-covariance, analyse the data and
+    return table of fitting parameters with a fit-plot.
+    We ignore Gly. Independent scaling of Macro Molecules."""
+
+    # Load input
+    synth_parameter = df_parameter_synth.mean().to_dict()
+    basis, names, Bheader = mrs_io.read_basis(basis_path)
+
+    # Set the adjustment lists from our input
+    broadening = [(synth_parameter["gamma_0"], synth_parameter["sigma_0"]),
+                  (synth_parameter["gamma_1"], synth_parameter["sigma_1"])]
+    shifting = [(synth_parameter["eps_0"]), (synth_parameter["eps_1"])]
+    baseline = [synth_parameter["B_real_0"], synth_parameter["B_imag_0"],
+                synth_parameter["B_real_1"], synth_parameter["B_imag_1"],
+                synth_parameter["B_real_2"], synth_parameter["B_imag_2"]]
+    coilphase = [synth_parameter["Phi0"]]
+
+    # Generate synthetic data
+    fidS, headerS, concentrationsS = synth.syntheticFromBasisFile(basis_path,
+                                                                  concentrations=synth_parameter,
+                                                                  ignore=['Gly'], ind_scaling=['mm'],
+                                                                  metab_groups='mm', broadening=broadening,  shifting=shifting,
+                                                                  # correct for complex noise
+                                                                  noisecovariance=np.divide(
+                                                                      noise_cov, 2),
+                                                                  # CAVE: baseline chosen manually
+                                                                  baseline=baseline, baseline_ppm=(
+                                                                      .2, 4.2),
+                                                                  coilphase=coilphase,
+                                                                  bandwidth=6000)
+    # Create mrs object for further use
+    mrsA = MRS(FID=fidS, header=headerS, basis=basis,
+               basis_hdr=Bheader[0], names=names)
+
+    mrsA.ignore(['Gly'])
+    mrsA.processForFitting(ind_scaling=['mm'])
+
+    # Scale it to Input data
+    # rescaled_FID, __ = rescale_FID(mrsA.FID, 1/mrsA.scaling["FID"]*100)
+    mrsA.set_FID(fidS)
+
+    metab_groups = parse_metab_groups(mrsA, 'mm')
+    # Use Voigt line broadening, fit between .2 and 4.2,
+    #  with a 2nd order polynomial baseline
+    FitArgs = {
+        'model': 'voigt',
+        'metab_groups': metab_groups,
+        'ppmlim': (.2, 4.2),
+        'baseline_order': 2}
+    res = fit_FSLModel(mrsA, **FitArgs)
+
+    # Plot fitting results.
+    if verbose:
+        fig = plotting.plot_fit(mrsA, pred=res.pred,
+                                baseline=res.baseline, out="name")
+
+    # Combine highly correlated metabolites
+    combinationList = [['Glu', 'Gln'],
+                       ['GPC', 'PCho'],
+                       ['Cr', 'PCr'],
+                       ['Glc', 'Tau'],
+                       ["NAA", "NAAG"]]
+    res.combine(combinationList)
+
+    # Store parameters from fit
+    fit_parameters.append(res.fitResults)
+    fit_snr.append(res.SNR.spectrum)
+    # Old version of noise
+    # noise_sd=np.max(np.real(mrsA.get_spec(ppmlim=(.2,4.2))))/res.SNR.spectrum
+    noise_sd = np.std(mrsA.FID[1000:1600])
+    fit_sdnoise.append(noise_sd)
+    noise_var = np.var(mrsA.FID[1000:1600])
+    fit_varnoise.append(noise_var)
+
+    df_params = pd.concat(fit_parameters, ignore_index=True)
+    df_params['SNR'] = fit_snr
+    df_params['noise_sd'] = fit_sdnoise
+    df_params['noise_var'] = fit_varnoise
+
+    if verbose:
+        return df_params, fig
+
+    return df_params