Correct the simulation module

The simulation module was broken: depending on the values given in the
fake profiles, that are necessary to run the simulation, the
back-calculated R2 values were offset differently. This is now fixed.

The calculated profiles are also available in the `Data` directory.

chemex/chemex.py

@@ -21,7 +21,7 @@
 from chemex.models import model
 from chemex.models.loader import register_kinetic_settings
 from chemex.optimize.fitting import run_methods
-from chemex.optimize.helper import execute_post_fit
+from chemex.optimize.helper import execute_simulation
 from chemex.parameters import database
 
 
@@ -46,8 +46,9 @@ def run_sim(args: Namespace, experiments: Experiments):
     plot = args.plot == "normal"
 
     database.fix_all_parameters()
+    experiments.back_calculate()
 
-    execute_post_fit(experiments, path, plot=plot)
+    execute_simulation(experiments, path, plot=plot)
 
 
 def run(args: Namespace):

chemex/containers/data.py

@@ -45,6 +45,15 @@ def scale_calc(self) -> None:
         self.calc *= scale
         self.scale = scale
 
+    def prepare_for_simulation(self) -> None:
+        self.exp = self.calc.copy()
+        self.err *= 0.0
+        if any(self.refs):
+            scale = float(100.0 / np.mean(self.exp[self.refs]))
+            self.exp *= scale
+            self.calc *= scale
+            self.scale *= scale
+
     def monte_carlo(self) -> Data:
         """Generate a data set to run Monte-Carlo simulation"""
         data = deepcopy(self)

chemex/containers/experiment.py

@@ -38,6 +38,9 @@ def residuals(self, params: ParametersLF) -> list[float]:
     def plot(self, path: Path):
         self.plotter.plot(path, self.profiles)
 
+    def plot_simulation(self, path: Path):
+        self.plotter.plot_simulation(path, self.profiles)
+
     def write(self, path: Path):
         filename = (path / self.filename.name).with_suffix(".dat")
         with filename.open("w") as file_dat:
@@ -93,6 +96,10 @@ def estimate_noise(self, kind: str):
         for profile in self.profiles:
             profile.set_noise(noise_mean)
 
+    def prepare_for_simulation(self) -> None:
+        for profile in self.profiles:
+            profile.prepare_for_simulation()
+
     def monte_carlo(self) -> Experiment:
         profiles = [profile.monte_carlo() for profile in self.profiles]
         return Experiment(

chemex/containers/experiments.py

@@ -39,6 +39,15 @@ def residuals(self, params: Parameters) -> np.ndarray:
             )
         )
 
+    def back_calculate(self) -> None:
+        params_lf = database.build_lmfit_params(self.param_ids)
+        self.residuals(params_lf)
+
+    def prepare_for_simulation(self) -> None:
+        self.back_calculate()
+        for experiment in self:
+            experiment.prepare_for_simulation()
+
     def write(self, path: Path):
         path_dat = path / "Data"
         path_dat.mkdir(parents=True, exist_ok=True)
@@ -49,6 +58,10 @@ def plot(self, path: Path):
         for experiment in self:
             experiment.plot(path)
 
+    def plot_simulation(self, path: Path):
+        for experiment in self:
+            experiment.plot_simulation(path)
+
     def select(self, selection: Selection):
         if selection.include is None and selection.exclude is None:
             return

chemex/containers/profile.py

@@ -94,6 +94,9 @@ def filter(self, params: ParametersLF) -> None:
     def set_noise(self, value: float):
         self.data.err[:] = value
 
+    def prepare_for_simulation(self) -> None:
+        self.data.prepare_for_simulation()
+
     def monte_carlo(self) -> Profile:
         profile = deepcopy(self)
         profile.data = profile.data.monte_carlo()

chemex/optimize/gridding.py

@@ -21,6 +21,7 @@
 from chemex.optimize.grouping import create_groups
 from chemex.optimize.grouping import Group
 from chemex.optimize.helper import calculate_statistics
+from chemex.optimize.helper import execute_post_fit
 from chemex.optimize.helper import execute_post_fit_groups
 from chemex.optimize.helper import print_header
 from chemex.optimize.helper import print_values
@@ -56,7 +57,8 @@ def run_group_grid(
     shape = tuple(len(values) for values in group_grid.values())
     grid_size = np.prod(shape)
 
-    filename = path / "Grid" / f"{group.path}.out"
+    basename = group.path if group.path != Path(".") else Path("grid")
+    filename = path / "Grid" / f"{basename}.out"
     filename.parent.mkdir(parents=True, exist_ok=True)
 
     with filename.open("w") as fileout:
@@ -256,3 +258,5 @@ def run_grid(
 
     if len(groups) > 1:
         execute_post_fit_groups(experiments, path, plot)
+    else:
+        execute_post_fit(experiments, path, plot != "nothing")

chemex/optimize/helper.py

@@ -80,6 +80,14 @@ def _write_files(experiments: Experiments, path: Path):
     _write_statistics(experiments, path=path)
 
 
+def _write_simulation_files(experiments: Experiments, path: Path):
+    """Write the results of the simulation to output files."""
+    print_writing_results(path)
+    path.mkdir(parents=True, exist_ok=True)
+    write_parameters(experiments, path)
+    experiments.write(path)
+
+
 def _write_plots(experiments: Experiments, path: Path):
     """Plot the experimental and fitted data."""
 
@@ -94,6 +102,20 @@ def _write_plots(experiments: Experiments, path: Path):
     print("")
 
 
+def _write_simulation_plots(experiments: Experiments, path: Path):
+    """Plot the experimental and fitted data."""
+
+    print_making_plots()
+
+    path_ = path / "Plots"
+    path_.mkdir(parents=True, exist_ok=True)
+    try:
+        experiments.plot_simulation(path=path_)
+    except KeyboardInterrupt:
+        print("  - Plotting cancelled\n")
+    print("")
+
+
 def execute_post_fit(
     experiments: Experiments,
     path: Path,
@@ -104,6 +126,17 @@ def execute_post_fit(
         _write_plots(experiments, path)
 
 
+def execute_simulation(
+    experiments: Experiments,
+    path: Path,
+    plot: bool = False,
+):
+    experiments.prepare_for_simulation()
+    _write_simulation_files(experiments, path)
+    if plot:
+        _write_simulation_plots(experiments, path)
+
+
 def execute_post_fit_groups(experiments: Experiments, path: Path, plot: str) -> None:
     print_group_name("All groups")
     params_lf = database.build_lmfit_params(experiments.param_ids)

chemex/plotters/cest.py

@@ -220,3 +220,20 @@ def plot(self, path: Path, profiles: list[Profile]) -> None:
                 self._plot_profile(file_pdf, profile, data_exp, data_calc)
                 file_exp.write(self.printer.print_exp(str(profile.name), data_exp))
                 file_calc.write(self.printer.print_calc(str(profile.name), data_calc))
+
+    def plot_simulation(self, path: Path, profiles: list[Profile]) -> None:
+
+        basename = path / self.filename.name
+        name_pdf = basename.with_suffix(".pdf")
+        name_sim = basename.with_suffix(".sim")
+
+        print_plot_filename(name_pdf, extra=False)
+
+        with ExitStack() as stack:
+            file_pdf = stack.enter_context(PdfPages(str(name_pdf)))
+            file_sim = stack.enter_context(name_sim.open("w"))
+            for profile in sorted(profiles):
+                data_exp = create_plot_data_exp(profile)
+                data_calc = create_plot_data_calc(profile)
+                self._plot_profile(file_pdf, profile, data_exp, data_calc)
+                file_sim.write(self.printer.print_calc(str(profile.name), data_calc))

chemex/plotters/cpmg.py

@@ -65,10 +65,8 @@ def calculate_exp_rates(data: Data, time_t2: float) -> np.ndarray:
 
 
 def calculate_calc_rates(data: Data, time_t2: float) -> np.ndarray:
-
     intensities = data.calc[~data.refs]
     intensities0 = data.exp[data.refs]
-
     return intensities_to_rates(intensities, intensities0, time_t2)
 
 
@@ -174,3 +172,20 @@ def plot(self, path: Path, profiles: list[Profile]) -> None:
                 plot_cpmg(file_pdf, str(profile.name), data_exp, data_calc)
                 file_exp.write(self.printer.print_exp(str(profile.name), data_exp))
                 file_calc.write(self.printer.print_calc(str(profile.name), data_calc))
+
+    def plot_simulation(self, path: Path, profiles: list[Profile]) -> None:
+
+        basename = path / self.filename.name
+        name_pdf = basename.with_suffix(".pdf")
+        name_sim = basename.with_suffix(".sim")
+
+        print_plot_filename(name_pdf, extra=False)
+
+        with ExitStack() as stack:
+            file_pdf = stack.enter_context(PdfPages(str(name_pdf)))
+            file_sim = stack.enter_context(name_sim.open("w"))
+            for profile in sorted(profiles):
+                data_exp = create_plot_data_exp(profile, self.config)
+                data_calc = create_plot_data_calc(profile, self.config)
+                plot_cpmg(file_pdf, str(profile.name), data_exp, data_calc)
+                file_sim.write(self.printer.print_calc(str(profile.name), data_calc))

chemex/plotters/plotter.py

@@ -9,3 +9,6 @@
 class Plotter(Protocol):
     def plot(self, path: Path, profiles: list[Profile]) -> None:
         ...
+
+    def plot_simulation(self, path: Path, profiles: list[Profile]) -> None:
+        ...

chemex/plotters/relaxation.py

@@ -85,3 +85,20 @@ def plot(self, path: Path, profiles: list[Profile]) -> None:
                 plot_relaxation(file_pdf, str(profile.name), data_exp, data_calc)
                 file_exp.write(self.printer.print_exp(str(profile.name), data_exp))
                 file_calc.write(self.printer.print_calc(str(profile.name), data_calc))
+
+    def plot_simulation(self, path: Path, profiles: list[Profile]) -> None:
+
+        basename = path / self.filename.name
+        name_pdf = basename.with_suffix(".pdf")
+        name_sim = basename.with_suffix(".sim")
+
+        print_plot_filename(name_pdf, extra=False)
+
+        with ExitStack() as stack:
+            file_pdf = stack.enter_context(PdfPages(str(name_pdf)))
+            file_sim = stack.enter_context(name_sim.open("w"))
+            for profile in sorted(profiles):
+                data_exp = create_plot_data_exp(profile)
+                data_calc = create_plot_data_calc(profile)
+                plot_relaxation(file_pdf, str(profile.name), data_exp, data_calc)
+                file_sim.write(self.printer.print_calc(str(profile.name), data_calc))

chemex/plotters/shift.py

@@ -42,6 +42,7 @@ def create_plot_data(profiles: list[Profile]) -> Data:
 class ShiftPlotter:
     def __init__(self, filename: Path, **_extra: Any):
         self.filename = filename
+        self.plot_simulation = self.plot
 
     def plot(self, path: Path, profiles: list[Profile]) -> None:
 

chemex/printers/plot.py

@@ -31,7 +31,7 @@ def print_exp(self, name: str, data: Data) -> str:
             start = " " if mask else "#"
             end = "" if mask else " # NOT USED IN THE FIT"
             output.append(
-                f"{start} {metadata: {self.first_column_fmt}}  {exp: 17.8e} {err[0]: 17.8e}{end}"
+                f"{start} {metadata: {self.first_column_fmt}} {exp: 17.8e} {err[0]: 17.8e}{end}"
             )
         return "\n".join(output) + "\n\n"
 
@@ -42,7 +42,7 @@ def print_calc(self, name: str, data: Data) -> str:
         ]
 
         for metadata, calc in zip(data.metadata, data.calc):
-            output.append(f"  {metadata: {self.first_column_fmt}}  {calc: 17.8e}")
+            output.append(f"  {metadata: {self.first_column_fmt}} {calc: 17.8e}")
 
         return "\n".join(output) + "\n\n"
 
@@ -64,15 +64,15 @@ def print_exp(self, name: str, data: Data) -> str:
             start = " " if mask else "#"
             end = "" if mask else " # NOT USED IN THE FIT"
             output.append(
-                f"{start} {metadata: {self.first_column_fmt}}  {exp: 17.8e} {err[0]: 17.8e} {err[1]: 17.8e}{end}"
+                f"{start} {metadata: {self.first_column_fmt}} {exp: 17.8e} {err[0]: 17.8e} {err[1]: 17.8e}{end}"
             )
         return "\n".join(output) + "\n\n"
 
     def print_calc(self, name: str, data: Data) -> str:
         output = [f"[{name}]", f"# {self.first_column_name:>12s} {'R2 (CALC)':>17s}"]
 
         for metadata, calc in zip(data.metadata, data.calc):
-            output.append(f"  {metadata: {self.first_column_fmt}}  {calc: 17.8e}")
+            output.append(f"  {metadata: {self.first_column_fmt}} {calc: 17.8e}")
 
         return "\n".join(output) + "\n\n"
 

chemex/tools/pick_cest.py

@@ -31,7 +31,7 @@ def __init__(self, profile: Profile, sw: float | None = None):
         delta_ppm = 0.0
         if settings is not None:
             spectrometer = profile.spectrometer
-            cos_n = getattr(settings, "cos_n")
+            cos_n = getattr(settings, "cos_n", None)
             sw = getattr(settings, "sw")
             if cos_n is not None and cos_n % 2 == 0 and sw is not None:
                 shifts_ppm = spectrometer.offsets_to_ppms(np.array([sw / 2.0, 0.0]))