135 lines
4.5 KiB
Python
135 lines
4.5 KiB
Python
# mc_gen.py - Generating data for Monte-Carlo style simulations
|
|
#
|
|
# Author: Konstantin E Bosbach <konstantin.bosbach@mars.uni-freiburg.de>
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
import time
|
|
|
|
from fsl_mrs.core import MRS
|
|
from fsl_mrs.utils import 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):
|
|
"""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
|
|
|
|
|
|
def mc(
|
|
n, noise_sd, df_parameter_synth, basis_path, output_path,
|
|
fit_parameters=[], fit_snr=[], fit_sdnoise=[], fit_varnoise=[]
|
|
):
|
|
"""Function for calling synth_and_ana repeatedly,
|
|
as in Monte-Carlo approach"""
|
|
runtime = time.time()
|
|
|
|
# Define workspace output path
|
|
file_out_path = str(
|
|
output_path + "noise_sd_" +
|
|
str(round(noise_sd, 3)) + "_runs_"+str(n) + ".csv"
|
|
)
|
|
|
|
print(
|
|
"Starting noise_sd", round(noise_sd, 2), " with ",
|
|
round(n, 2), "repetitions"
|
|
)
|
|
|
|
# Call function generation the desired amount of times
|
|
for k in range(0, n):
|
|
noise_fit = synth_and_ana(
|
|
[[np.square(noise_sd)]],
|
|
fit_parameters=fit_parameters,
|
|
fit_snr=fit_snr, fit_sdnoise=fit_sdnoise,
|
|
fit_varnoise=fit_varnoise,
|
|
df_parameter_synth=df_parameter_synth, basis_path=basis_path
|
|
)
|
|
|
|
# Print results to file
|
|
noise_fit.to_csv(file_out_path)
|
|
|
|
print(str("Finishing noise_sd", round(noise_sd, 2), " with ",
|
|
str(round(n, 2)), "repetitions, Runtime took ",
|
|
round(time.time()-runtime, 2), "[s]"))
|
|
|
|
return noise_fit, file_out_path
|