From e29ec621f59efc778346f75e2cd90fae7367dfa3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoitsierro@bluewin.ch>
Date: Wed, 25 Aug 2021 12:30:55 +0200
Subject: [PATCH 1/7] added version command

---
 src/scgenerator/cli/cli.py | 3 ++-
 src/scgenerator/const.py   | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/scgenerator/cli/cli.py b/src/scgenerator/cli/cli.py
index 666a3f8..25b0320 100644
--- a/src/scgenerator/cli/cli.py
+++ b/src/scgenerator/cli/cli.py
@@ -6,7 +6,7 @@ import re
 
 import numpy as np
 
-from .. import env, io, scripts
+from .. import env, io, scripts, const
 from ..logger import get_logger
 from ..physics.fiber import dispersion_coefficients
 from ..physics.simulate import SequencialSimulations, resume_simulations, run_simulation_sequence
@@ -37,6 +37,7 @@ def create_parser():
         parser.add_argument(
             *names, **{k: v for k, v in args.items() if k not in {"short_name", "type"}}
         )
+    parser.add_argument("--version", action="version", version=const.__version__)
 
     run_parser = subparsers.add_parser("run", help="run a simulation from a config file")
     run_parser.add_argument("configs", help="path(s) to the toml configuration file(s)", nargs="+")
diff --git a/src/scgenerator/const.py b/src/scgenerator/const.py
index e650575..a8e391d 100644
--- a/src/scgenerator/const.py
+++ b/src/scgenerator/const.py
@@ -1,4 +1,4 @@
-__version__ = "0.1.1"
+__version__ = "0.1.1dev"
 
 
 from typing import Any

From bc95f6d8d57004019a44f2bc211c749f145933be Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoitsierro@bluewin.ch>
Date: Thu, 26 Aug 2021 09:21:29 +0200
Subject: [PATCH 2/7] Added Noise correlation parameter

---
 src/scgenerator/const.py           | 2 +-
 src/scgenerator/defaults.py        | 1 +
 src/scgenerator/initialize.py      | 7 +++++--
 src/scgenerator/physics/pulse.py   | 8 +++++---
 src/scgenerator/utils/parameter.py | 1 +
 5 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/src/scgenerator/const.py b/src/scgenerator/const.py
index a8e391d..bcbadaf 100644
--- a/src/scgenerator/const.py
+++ b/src/scgenerator/const.py
@@ -1,4 +1,4 @@
-__version__ = "0.1.1dev"
+__version__ = "0.1.2dev"
 
 
 from typing import Any
diff --git a/src/scgenerator/defaults.py b/src/scgenerator/defaults.py
index d3fed53..a1dfdec 100644
--- a/src/scgenerator/defaults.py
+++ b/src/scgenerator/defaults.py
@@ -20,6 +20,7 @@ default_parameters = dict(
     shape="gaussian",
     frep=40e6,
     behaviors=["spm", "ss"],
+    noise_correlation=0,
     raman_type="agrawal",
     parallel=True,
     repeat=1,
diff --git a/src/scgenerator/initialize.py b/src/scgenerator/initialize.py
index 423390f..dad4a43 100644
--- a/src/scgenerator/initialize.py
+++ b/src/scgenerator/initialize.py
@@ -81,7 +81,9 @@ class Params(BareParams):
 
         # Technical noise
         if self.intensity_noise is not None and self.intensity_noise > 0:
-            delta_int, delta_T0 = pulse.technical_noise(self.intensity_noise)
+            delta_int, delta_T0 = pulse.technical_noise(
+                self.intensity_noise, self.noise_correlation
+            )
             self.peak_power *= delta_int
             self.t0 *= delta_T0
             self.width *= delta_T0
@@ -203,7 +205,7 @@ class Config(BareConfig):
             if self.loss == "capillary":
                 for param in ["core_radius", "he_mode"]:
                     self.get_fiber(param)
-        for param in ["length", "input_transmission"]:
+        for param in ["length", "input_transmission", "n2"]:
             self.get(param)
 
     def gas_consistency(self):
@@ -244,6 +246,7 @@ class Config(BareConfig):
             "interpolation_degree",
             "ideal_gas",
             "recovery_last_stored",
+            "noise_correlation",
         ]:
             self.get(param)
 
diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py
index 63cc9e3..7f9dab0 100644
--- a/src/scgenerator/physics/pulse.py
+++ b/src/scgenerator/physics/pulse.py
@@ -139,6 +139,7 @@ def modify_field_ratio(
     target_power: float = None,
     target_energy: float = None,
     intensity_noise: float = None,
+    noise_correlation: float = None,
 ) -> float:
     """multiply a field by this number to get the desired effects
 
@@ -165,7 +166,7 @@ def modify_field_ratio(
         ratio *= np.sqrt(target_power / abs2(field).max())
 
     if intensity_noise is not None:
-        d_int, _ = technical_noise(intensity_noise)
+        d_int, _ = technical_noise(intensity_noise, noise_correlation)
         ratio *= np.sqrt(d_int)
     return ratio
 
@@ -318,6 +319,7 @@ def setup_custom_field(params: BareParams) -> bool:
             params.peak_power,
             params.energy,
             params.intensity_noise,
+            params.noise_correlation,
         )
         width, peak_power, energy = measure_field(params.t, field_0)
     else:
@@ -374,7 +376,7 @@ def pulse_energy_with_loss(spectrum, dw, alpha, h) -> float:
     return np.sum(spec2 * dw) - h * np.sum(alpha * spec2 * dw)
 
 
-def technical_noise(rms_noise, relative_factor=0.4):
+def technical_noise(rms_noise, noise_correlation=-0.4):
     """
     To implement technical noise as described in Grenier2019, we need to know the
     noise properties of the laser, summarized into the RMS amplitude noise
@@ -391,7 +393,7 @@ def technical_noise(rms_noise, relative_factor=0.4):
         delta_T0 : float
     """
     psy = np.random.normal(1, rms_noise)
-    return psy, 1 - relative_factor * (psy - 1)
+    return psy, 1 + noise_correlation * (psy - 1)
 
 
 def shot_noise(w_c, w0, T, dt):
diff --git a/src/scgenerator/utils/parameter.py b/src/scgenerator/utils/parameter.py
index 338840d..37ba1f9 100644
--- a/src/scgenerator/utils/parameter.py
+++ b/src/scgenerator/utils/parameter.py
@@ -392,6 +392,7 @@ class BareParams:
     shape: str = Parameter(literal("gaussian", "sech"))
     wavelength: float = Parameter(in_range_incl(100e-9, 3000e-9), display_info=(1e9, "nm"))
     intensity_noise: float = Parameter(in_range_incl(0, 1), display_info=(1e2, "%"))
+    noise_correlation: float = Parameter(in_range_incl(-10, 10))
     width: float = Parameter(in_range_excl(0, 1e-9), display_info=(1e15, "fs"))
     t0: float = Parameter(in_range_excl(0, 1e-9), display_info=(1e15, "fs"))
 

From 24bb7866cdeade260e7026e41a84a7ec690ea09f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoitsierro@bluewin.ch>
Date: Thu, 26 Aug 2021 14:10:17 +0200
Subject: [PATCH 3/7] marcuse equation

---
 src/scgenerator/physics/fiber.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/scgenerator/physics/fiber.py b/src/scgenerator/physics/fiber.py
index af63faf..c91a430 100644
--- a/src/scgenerator/physics/fiber.py
+++ b/src/scgenerator/physics/fiber.py
@@ -310,6 +310,7 @@ def A_eff_marcuse(wl: T, core_radius: float, numerical_aperture: float) -> T:
     """
     V = 2 * pi * core_radius * numerical_aperture / wl
     w_eff = core_radius * (0.65 + 1.619 / V ** 1.5 + 2.879 / V ** 6)
+    return np.pi * w_eff ** 2
 
 
 def HCPCF_find_with_given_ZDW(

From 9375968125270b86946d415b12a57fb7083982e0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoitsierro@bluewin.ch>
Date: Thu, 26 Aug 2021 14:45:28 +0200
Subject: [PATCH 4/7] alert at the end of simulation

---
 src/scgenerator/cli/cli.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/src/scgenerator/cli/cli.py b/src/scgenerator/cli/cli.py
index 25b0320..ff66523 100644
--- a/src/scgenerator/cli/cli.py
+++ b/src/scgenerator/cli/cli.py
@@ -3,7 +3,8 @@ import os
 from collections import ChainMap
 from pathlib import Path
 import re
-
+import sys
+import subprocess
 import numpy as np
 
 from .. import env, io, scripts, const
@@ -148,6 +149,16 @@ def run_sim(args):
 
     method = prep_ray()
     run_simulation_sequence(*args.configs, method=method)
+    if sys.platform == "darwin" and sys.stdout.isatty():
+        subprocess.run(
+            [
+                "osascript",
+                "-e",
+                'tell app "System Events" to display dialog "simulation finished !"',
+            ],
+            stdout=subprocess.DEVNULL,
+            stderr=subprocess.DEVNULL,
+        )
 
 
 def merge(args):

From f8d2f53083d0395662114212f281d5bd37c910fa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoitsierro@bluewin.ch>
Date: Thu, 26 Aug 2021 15:00:49 +0200
Subject: [PATCH 5/7] fixed energy conservation

---
 src/scgenerator/physics/simulate.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/scgenerator/physics/simulate.py b/src/scgenerator/physics/simulate.py
index a6107b3..28719eb 100644
--- a/src/scgenerator/physics/simulate.py
+++ b/src/scgenerator/physics/simulate.py
@@ -114,12 +114,12 @@ class RK4IP:
             elif self.alpha is not None:
                 self.logger.debug("Conserved quantity : energy with loss")
                 self.conserved_quantity_func = lambda spectrum, h: pulse.pulse_energy_with_loss(
-                    spectrum, self.dw, self.alpha, h
+                    self.C_to_A_factor * spectrum, self.dw, self.alpha, h
                 )
             else:
                 self.logger.debug("Conserved quantity : energy without loss")
                 self.conserved_quantity_func = lambda spectrum, h: pulse.pulse_energy(
-                    spectrum, self.dw
+                    self.C_to_A_factor * spectrum, self.dw
                 )
         else:
             self.conserved_quantity_func = lambda spectrum, h: 0.0

From 5236beedd6a06ddfd521de9364fbab20135abb55 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoitsierro@bluewin.ch>
Date: Fri, 27 Aug 2021 08:37:58 +0200
Subject: [PATCH 6/7] iterable RK4IP solver, also in parallel

---
 src/scgenerator/initialize.py       |   5 ++
 src/scgenerator/physics/simulate.py | 124 +++++++++++++++++++++++-----
 2 files changed, 107 insertions(+), 22 deletions(-)

diff --git a/src/scgenerator/initialize.py b/src/scgenerator/initialize.py
index dad4a43..18a642f 100644
--- a/src/scgenerator/initialize.py
+++ b/src/scgenerator/initialize.py
@@ -361,6 +361,11 @@ class ParamSequence:
     def count_variations(self) -> int:
         return count_variations(self.config)
 
+    @property
+    def first(self) -> Params:
+        for _, params in self:
+            return params
+
 
 class ContinuationParamSequence(ParamSequence):
     def __init__(self, prev_sim_dir: os.PathLike, new_config: BareConfig):
diff --git a/src/scgenerator/physics/simulate.py b/src/scgenerator/physics/simulate.py
index 28719eb..a253cc0 100644
--- a/src/scgenerator/physics/simulate.py
+++ b/src/scgenerator/physics/simulate.py
@@ -1,9 +1,10 @@
 import multiprocessing
+import multiprocessing.connection
 import os
 import random
 from datetime import datetime
 from pathlib import Path
-from typing import Dict, List, Tuple, Type, Union
+from typing import Dict, Generator, List, Tuple, Type, Any
 
 import numpy as np
 
@@ -53,14 +54,18 @@ class RK4IP:
 
         self.job_identifier = job_identifier
         self.id = task_id
+        self.save_data = save_data
 
-        self.sim_dir = io.get_sim_dir(self.id)
-        self.sim_dir.mkdir(exist_ok=True)
-        self.data_dir = self.sim_dir / self.job_identifier
+        if self.save_data:
+            self.sim_dir = io.get_sim_dir(self.id)
+            self.sim_dir.mkdir(exist_ok=True)
+            self.data_dir = self.sim_dir / self.job_identifier
+        else:
+            self.sim_dir = None
+            self.data_dir = None
 
         self.logger = get_logger(self.job_identifier)
         self.resuming = False
-        self.save_data = save_data
 
         self.w_c = params.w_c
         self.w = params.w
@@ -144,9 +149,6 @@ class RK4IP:
         ]
         self.size_fac = 2 ** (1 / 5)
 
-        if self.save_data:
-            self._save_current_spectrum(0)
-
         # Initial step size
         if self.adapt_step_size:
             self.initial_h = (self.z_targets[0] - self.z) / 2
@@ -165,6 +167,16 @@ class RK4IP:
         self._save_data(self.cons_qty, f"cons_qty")
         self.step_saved()
 
+    def get_current_spectrum(self) -> tuple[int, np.ndarray]:
+        """returns the current spectrum
+
+        Returns
+        -------
+        np.ndarray
+            spectrum
+        """
+        return self.C_to_A_factor * self.current_spectrum
+
     def _save_data(self, data: np.ndarray, name: str):
         """calls the appropriate method to save data
 
@@ -178,6 +190,24 @@ class RK4IP:
         io.save_data(data, self.data_dir, name)
 
     def run(self):
+        time_start = datetime.today()
+
+        for step, num, _ in self.irun():
+            if self.save_data:
+                self._save_current_spectrum(num)
+
+        self.logger.info(
+            "propagation finished in {} steps ({} seconds)".format(
+                step, (datetime.today() - time_start).total_seconds()
+            )
+        )
+
+        if self.save_data:
+            self._save_data(self.z_stored, "z.npy")
+
+        return self.stored_spectra
+
+    def irun(self):
 
         # Print introduction
         self.logger.debug(
@@ -189,7 +219,8 @@ class RK4IP:
         h_taken = self.initial_h
         h_next_step = self.initial_h
         store = False  # store a spectrum
-        time_start = datetime.today()
+
+        yield step, len(self.stored_spectra) - 1, self.get_current_spectrum()
 
         while self.z < self.z_final:
             h_taken, h_next_step, self.current_spectrum = self.take_step(
@@ -205,8 +236,8 @@ class RK4IP:
                 self.logger.debug("{} steps, z = {:.4f}, h = {:.5g}".format(step, self.z, h_taken))
 
                 self.stored_spectra.append(self.current_spectrum)
-                if self.save_data:
-                    self._save_current_spectrum(len(self.stored_spectra) - 1)
+
+                yield step, len(self.stored_spectra) - 1, self.get_current_spectrum()
 
                 self.z_stored.append(self.z)
                 del self.z_targets[0]
@@ -225,17 +256,6 @@ class RK4IP:
                 store = True
                 h_next_step = self.z_targets[0] - self.z
 
-        self.logger.info(
-            "propagation finished in {} steps ({} seconds)".format(
-                step, (datetime.today() - time_start).total_seconds()
-            )
-        )
-
-        if self.save_data:
-            self._save_data(self.z_stored, "z.npy")
-
-        return self.stored_spectra
-
     def take_step(
         self, step: int, h_next_step: float, current_spectrum: np.ndarray
     ) -> Tuple[float, float, np.ndarray]:
@@ -731,5 +751,65 @@ def resume_simulations(sim_dir: Path, method: Type[Simulations] = None) -> Simul
     return Simulations.new(param_seq, task_id, method)
 
 
+def __parallel_RK4IP_worker(
+    worker_id: int,
+    msq_queue: multiprocessing.connection.Connection,
+    data_queue: multiprocessing.Queue,
+    params: utils.BareParams,
+):
+    logger = get_logger(__name__)
+    logger.debug(f"workder {worker_id} started")
+    for out in RK4IP(params).irun():
+        logger.debug(f"worker {worker_id} waiting for msg")
+        msq_queue.recv()
+        logger.debug(f"worker {worker_id} got msg")
+        data_queue.put((worker_id, out))
+        logger.debug(f"worker {worker_id} sent data")
+
+
+def parallel_RK4IP(
+    config,
+) -> Generator[
+    tuple[tuple[list[tuple[str, Any]], initialize.Params, int, int, np.ndarray], ...], None, None
+]:
+    logger = get_logger(__name__)
+    params = list(initialize.ParamSequence(config))
+    n = len(params)
+    z_num = params[0][1].z_num
+
+    cpu_no = multiprocessing.cpu_count()
+    if len(params) < cpu_no:
+        cpu_no = len(params)
+
+    pipes = [multiprocessing.Pipe(duplex=False) for i in range(n)]
+    data_queue = multiprocessing.Queue()
+    workers = [
+        multiprocessing.Process(target=__parallel_RK4IP_worker, args=(i, pipe[0], data_queue, p[1]))
+        for i, (pipe, p) in enumerate(zip(pipes, params))
+    ]
+    try:
+        [w.start() for w in workers]
+        logger.debug("pool started")
+        for i in range(z_num):
+            for q in pipes:
+                q[1].send(0)
+                logger.debug("msg sent")
+            computed_dict = {}
+            for j in range(n):
+                w_id, computed = data_queue.get()
+                computed_dict[w_id] = computed
+            computed_dict = list(computed_dict.items())
+            computed_dict.sort()
+            yield tuple((*p, *c) for p, c in zip(params, [el[1] for el in computed_dict]))
+        print("finished")
+    finally:
+        for w, cs in zip(workers, pipes):
+            w.join()
+            w.close()
+            cs[0].close()
+            cs[1].close()
+        data_queue.close()
+
+
 if __name__ == "__main__":
     pass

From 46999f15f7bbf26250153ec2d4044d172ea928e8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= <benoitsierro@bluewin.ch>
Date: Sat, 28 Aug 2021 14:15:10 +0200
Subject: [PATCH 7/7] miscelanious

---
 src/scgenerator/physics/fiber.py    | 7 +++++++
 src/scgenerator/physics/pulse.py    | 6 ++----
 src/scgenerator/physics/simulate.py | 2 +-
 3 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/scgenerator/physics/fiber.py b/src/scgenerator/physics/fiber.py
index c91a430..d8334ad 100644
--- a/src/scgenerator/physics/fiber.py
+++ b/src/scgenerator/physics/fiber.py
@@ -711,6 +711,13 @@ def compute_loss(params: BareParams) -> Optional[np.ndarray]:
         out = np.zeros_like(params.l)
         out[mask] = alpha
         return out
+    # else:
+        # return np.where(
+        #     (params.l < params.lower_wavelength_interp_limit)
+        #     | (params.l > params.upper_wavelength_interp_limit),
+        #     100,
+        #     0,
+        # )
     return None
 
 
diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py
index 7f9dab0..d1b086b 100644
--- a/src/scgenerator/physics/pulse.py
+++ b/src/scgenerator/physics/pulse.py
@@ -1019,10 +1019,8 @@ def rin_curve(spectra: np.ndarray) -> np.ndarray:
     rin_curve : np.ndarray
         RIN curve
     """
-    spec2 = abs2(spectra)
-    # return np.std(spec, axis=0) / np.mean(spec, axis=0)
-    m = np.mean(spec2, axis=0)
-    return np.sqrt(np.mean((spec2 - m) ** 2)) / m
+    A2 = abs2(spectra)
+    return np.std(A2, axis=0) / np.mean(A2, axis=0)
 
 
 def measure_field(t: np.ndarray, field: np.ndarray) -> Tuple[float, float, float]:
diff --git a/src/scgenerator/physics/simulate.py b/src/scgenerator/physics/simulate.py
index a253cc0..b10298d 100644
--- a/src/scgenerator/physics/simulate.py
+++ b/src/scgenerator/physics/simulate.py
@@ -207,7 +207,7 @@ class RK4IP:
 
         return self.stored_spectra
 
-    def irun(self):
+    def irun(self) -> Generator[tuple[int, int, np.ndarray], None, None]:
 
         # Print introduction
         self.logger.debug(