Typo

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,51 @@ 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=[],
+    shortage=False
+):
+    """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"
+    )
+
+    if shortage:
+        try:
+            noise_fit = pd.read_csv(file_out_path)
+            print("Found existing file: ", file_out_path)
+            return noise_fit, file_out_path
+        except:
+            None
+
+    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("Finishing noise_sd", round(noise_sd, 2), " with ",
+          round(n, 2), "repetitions, Runtime took ",
+          round(time.time()-runtime, 2), "[s]")
+
+    return noise_fit, file_out_path

fsl_mrs_mce/mc_in.py

@@ -0,0 +1,124 @@
+# mc_in.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
+from configparser import ConfigParser
+
+
+def load_config(input_folder="input/", section="default"):
+    """Loads configuration file.
+    If given option section="": returns list of config sections"""
+# Information from config txt
+    try:
+        # Load txt
+        configFilePath = str(input_folder + "config.txt")
+        configParser = ConfigParser()
+        configParser.read(configFilePath)
+        try:
+            if section == "":
+                return configParser.sections()
+            else:
+                return configParser[section]
+        except input:
+            print(configParser.sections())
+            print(f"No {section} section.")
+    except input:
+        print(f"No {configFilePath} provided.")
+
+
+def conc_df_from_file(input_path="input/", foundation="vivo_average",
+                      include_molecules=[], save_df=True):
+    sections = load_config(input_path, "")
+    """Generates concentration dataframe for mc method.
+    Supports different variation types (absolute/delta/relative)
+    Returns data_frame of different concentrations.
+    include_molecules needed for option some_molecules"""
+    # Load all relevant sections
+    list_parser_parameter = []
+
+    for i in sections:
+        if i[0:20] == "generation_parameter":
+            list_parser_parameter.append(load_config(input_path, i))
+
+    # Load foundation file, a line with e.g. in-vivo averages
+    # The code expects a one-line-foundation.
+    molecule_names = ['Ala', 'Asc', 'Asp', 'Cr', 'GABA', 'GPC', 'GSH',
+                      'Glc', 'Gln', 'Glu', 'Ins', 'Lac', 'NAA', 'NAAG',
+                      'PCho', 'PCr', 'PE', 'Scyllo', 'Tau',
+                      'mm', 'Glu+Gln', 'GPC+PCho', 'Cr+PCr',
+                      'Glc+Tau', 'NAA+NAAG']
+
+    if foundation == "vivo_average":
+        try:
+            df_foundation = pd.read_csv(
+                "basis/fit_conc_result.csv").mean().to_frame().transpose()
+        except:
+            print("No basis/fit_conc_result.csv file supplied")
+    # Sets all initial molecules, including some groups, to 0
+    if foundation == "flat":
+        try:
+            df_foundation = pd.read_csv(
+                "basis/fit_conc_result.csv").mean().to_frame().transpose()
+            for molecule in molecule_names:
+                df_foundation[molecule] = 0
+        except:
+            print("No basis/fit_conc_result.csv file supplied")
+    if foundation == "some_molecules":
+        try:
+            df_foundation = pd.read_csv(
+                "basis/fit_conc_result.csv").mean().to_frame().transpose()
+            for molecule in molecule_names:
+                if molecule not in include_molecules:
+                    df_foundation[molecule] = 0
+        except:
+            print("No basis/fit_conc_result.csv file supplied")
+    else:
+        print("Only foundation mode vivo_average not chosen")
+
+    # Load values, dependent on type absolute/delta/relative
+    list_parameter_values = []
+    list_parameter_names = []
+    for i in list_parser_parameter:
+        string_values = i["amount"][1:-1]
+        list_values = list(map(float, string_values.split(",")))
+        if i["type"] == "absolute":
+            # set absolute values
+            None
+
+        elif i["type"] == "delta":
+            # set absolute value offset from foundation
+            list_values = np.add(list_values, df_foundation[i["param"]].mean())
+
+        elif i["type"] == "relative":
+            # set relative value from foundation
+            list_values = np.multiply(
+                list_values, df_foundation[i["param"]].mean())
+
+        list_parameter_values.append(list_values)
+        list_parameter_names.append(i["param"])
+
+    # Create parameter tensor
+    mesh_parameter_values = np.meshgrid(*list_parameter_values)
+
+    # Create 'empty' standard data
+    df_conc = pd.concat(
+        [df_foundation]*int(np.array(mesh_parameter_values[0]).size),
+        ignore_index=True
+    )
+
+    # Write parameter values into dataframe
+    i = 0
+    while i <= len(list_parameter_values)-1:
+        parameter_values = mesh_parameter_values[i].ravel()
+        parameter_name = list_parameter_names[i]
+        print(parameter_name, parameter_values)
+        df_conc[parameter_name] = parameter_values
+        i = i+1
+
+    if save_df:
+        df_conc.to_csv(str(input_path+"df_conc.csv"))
+
+    return df_conc

fsl_mrs_mce/mc_sim.py

@@ -0,0 +1,49 @@
+# 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 data results from last two
+    files contain are within convergence criteria.
+    Returns Boolean"""
+
+    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)
+
+    # check if results within convergence criteria
+    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",
+        round(measure_mean, 4), round(measure_std, 4),
+        "convergence: ", str(convergence)
+    )
+
+    return convergence

setup.py

@@ -1,10 +1,10 @@
 from setuptools import setup
 
 setup(name='fsl_mrs_mce',
-      version='0.0.1',
+      version='0.0.22',
       description='A fsl_mrs Moncte Carlo estimation approach',
       url='https://git.thoffbauer.de/konstantin.bosbach/fsl_mrs_mce.git',
-      author='Konstantin Bosbach',
+      author='Konstantin E Bosbach',
       author_email='konstantin.bosbach@mars.uni-freiburg.de',
       packages=['fsl_mrs_mce'],
       zip_safe=False)