# Synthetise new data and analyse it directly to return fit parameters import numpy as np import pandas 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): """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) # 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 return df_params