added mc_sim.py, with get_convergence function for checking if experiment converges

fsl_mrs_mce/mc_gen.py

@@ -1,10 +1,13 @@
-#  Synthetise new data and analyse it directly to return fit parameters
+# 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 plotting, mrs_io
+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
@@ -31,18 +34,20 @@ def synth_and_ana(noise_cov,
     coilphase = [synth_parameter["Phi0"]]
 
     # Generate synthetic data
-    fidS, headerS, concentrationsS = synth.syntheticFromBasisFile(basis_path,
+    fidS, headerS, concentrationsS = synth.syntheticFromBasisFile(
-                                                                  concentrations=synth_parameter,
+        basis_path,
-                                                                  ignore=['Gly'], ind_scaling=['mm'],
+        concentrations=synth_parameter,
-                                                                  metab_groups='mm', broadening=broadening,  shifting=shifting,
+        ignore=['Gly'], ind_scaling=['mm'],
-                                                                  # correct for complex noise
+        metab_groups='mm', broadening=broadening,
-                                                                  noisecovariance=np.divide(
+        shifting=shifting,
-                                                                      noise_cov, 2),
+        # correct for complex noise
-                                                                  # CAVE: baseline chosen manually
+        noisecovariance=np.divide(noise_cov, 2),
-                                                                  baseline=baseline, baseline_ppm=(
+        # CAVE: baseline chosen manually
-                                                                      .2, 4.2),
+        baseline=baseline, baseline_ppm=(.2, 4.2),
-                                                                  coilphase=coilphase,
+        coilphase=coilphase,
-                                                                  bandwidth=6000)
+        bandwidth=6000
+    )
+
     # Create mrs object for further use
     mrsA = MRS(FID=fidS, header=headerS, basis=basis,
                basis_hdr=Bheader[0], names=names)
@@ -88,3 +93,42 @@ def synth_and_ana(noise_cov,
     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

fsl_mrs_mce/mc_sim.py

@@ -0,0 +1,45 @@
+# mc_sim.py - Helpers for Monte-Carlo style simulations
+#
+# Author: Konstantin E Bosbach <konstantin.bosbach@mars.uni-freiburg.de>
+
+
+import numpy as np
+import pandas as pd
+
+
+def get_convergence(
+    output_file_paths, molecule, crit_mean=0.01, crit_std=0.10
+):
+    """Function checks if last two files contain converging datasets."""
+
+    if len(output_file_paths) < 2:
+        print("One iteration, therefore no convergence.")
+        return False
+    else:
+        newer_fit_results = pd.read_csv(output_file_paths[-1])
+        older_fit_results = pd.read_csv(output_file_paths[-2])
+
+    newer_mean = newer_fit_results[molecule].mean()
+    older_mean = older_fit_results[molecule].mean()
+
+    # Normalize with mean of latest dataset, to get deviation
+    norm = newer_mean
+    measure_mean = abs(abs(np.abs(newer_mean) - abs(older_mean))/norm)
+
+    newer_std = newer_fit_results[molecule].std()
+    older_std = older_fit_results[molecule].std()
+    measure_std = abs(abs(np.abs(newer_std) - abs(older_std))/norm)
+
+    if measure_mean <= crit_mean:
+        if measure_std <= crit_std:
+            convergence = True
+    else:
+        convergence = False
+
+    print(
+        "Convergence result for ", output_file_paths[-1], " and ",
+        output_file_paths[-2], "\n\t\t\t",
+        measure_mean, measure_std, "convergence: ", str(convergence)
+    )
+
+    return convergence