Management of noise for 2 fibers coupler + correction docs

optcom/components/fiber_coupler.py

@@ -383,6 +383,8 @@ def __init__(self, name: str = default_name, length: float = 1.0,
         # Policy -------------------------------------------------------
         self.add_port_policy(([0, 1], [2, 3], True))
         self.add_wait_policy([0, 1], [2, 3])
+        # temp
+        self._get_kappa_for_noise = cnlse.get_kappa_for_noise
     # ==================================================================
     def output_ports(self, input_ports: List[int]) -> List[int]:
         # 1 input ------------------------------------------------------
@@ -420,14 +422,38 @@ def __call__(self, domain: Domain, ports: List[int], fields: List[Field]
                 else:
                     fields_port_1.append(fields[i])
         output_fields = self._stepper(domain, fields_port_0, fields_port_1)
-        # Noise management - splitting assumption
-        #if (self._NOISE):
-        #    noise = np.zeros(domain.noise_samples)
-        #    for i in range(len(ports)):
-        #        noise += fields[i].noise
-        #    noise /= len(ports)
-        #    for i in range(len(ports)):
-        #        fields[i].noise = noise
+        # Noise management
+        if (self._NOISE):
+            # 1. splitting assumption
+            #noise = np.zeros(domain.noise_samples)
+            #for i in range(len(ports)):
+            #    noise += fields[i].noise
+            #noise /= len(ports)
+            #for i in range(len(ports)):
+            #    fields[i].noise = noise
+            # 2. Anal sol, only bi-fiber coupler config supported for now
+            # and make assumptions of a symmetric coupler
+            # k_{12} = k_{21}.  See Agrawal applications top page 61 for
+            # analytical solution
+            length = self._stepper._length
+            kappa = self._get_kappa_for_noise()
+            factor = np.power(np.cos(kappa*length), 2)
+            factor_ = np.power(np.sin(kappa*length), 2)
+            accu_noise_port_0 = np.zeros(domain.noise_samples)
+            accu_noise_port_1 = np.zeros(domain.noise_samples)
+            # Solution of for ecah considered alone
+            for i, field in enumerate(fields_port_0):   # len == 2
+                accu_noise_port_0 += factor_ * field.noise
+                field.noise *= factor
+            for i, field in enumerate(fields_port_1):   # len == 2
+                accu_noise_port_1 += factor_ * field.noise
+                field.noise *= factor
+            # Redistributing the noise that switched in the other fiber
+            for field in fields_port_0:
+                field.noise += accu_noise_port_1 / len(fields_port_0)
+            for field in fields_port_1:
+                field.noise += accu_noise_port_0 / len(fields_port_1)
+
         # Get storage
         if (self._stepper.save_all):
             self.storages.append(self._stepper.storage)
@@ -448,13 +474,17 @@ def __call__(self, domain: Domain, ports: List[int], fields: List[Field]
 
     import optcom as oc
 
-    lt: oc.Layout = oc.Layout(oc.Domain(bit_width=20.0, samples_per_bit=1024))
+    noise_samples: int = 100
+    lt: oc.Layout = oc.Layout(oc.Domain(bit_width=20.0, samples_per_bit=1024,
+                                        noise_samples=noise_samples))
 
     Lambda: float = 1030.0
     pulse_1: oc.Gaussian = oc.Gaussian(channels=1, peak_power=[38.5, 0.5],
-                                       fwhm=[1.], center_lambda=[Lambda])
+                                       fwhm=[1.], center_lambda=[Lambda],
+                                       noise=np.ones(noise_samples)*12)
     pulse_2: oc.Gaussian = oc.Gaussian(channels=1, peak_power=[23.5, 0.3],
-                                       fwhm=[1.], center_lambda=[1050.0])
+                                       fwhm=[1.], center_lambda=[1050.0],
+                                       noise=np.ones(noise_samples)*5)
 
     steps: int = int(100)
     alpha: List[Union[List[float], Callable, None]] = [[0.046], [0.046]]
@@ -513,17 +543,27 @@ def __call__(self, domain: Domain, ports: List[int], fields: List[Field]
     y_datas: List[np.ndarray] = [oc.temporal_power(pulse_1[0][0].channels),
                                  oc.temporal_power(pulse_2[0][0].channels),
                                  oc.temporal_power(coupler[2][0].channels),
-                                 oc.temporal_power(coupler[3][0].channels)]
+                                 oc.temporal_power(coupler[3][0].channels),
+                                 pulse_1[0][0].noise, coupler[2][0].noise,
+                                 pulse_2[0][0].noise, coupler[3][0].noise]
     x_datas: List[np.ndarray] = [pulse_1[0][0].time, pulse_2[0][0].time,
-                                 coupler[2][0].time, coupler[3][0].time]
-    plot_groups: List[int] = [0, 1, 2, 3]
-    plot_titles: List[str] = ["Original pulse", "Original pulse bis",
+                                 coupler[2][0].time, coupler[3][0].time,
+                                 pulse_1[0][0].domain.noise_omega,
+                                 coupler[2][0].domain.noise_omega,
+                                 pulse_2[0][0].domain.noise_omega,
+                                 coupler[3][0].domain.noise_omega]
+    plot_groups: List[int] = [0, 1, 2, 3, 4, 4, 5, 5]
+    plot_titles: List[str] = ["Original pulse", "Original pulse",
                               "Pulse coming out of the "
                               "coupler with Lk = {}"
                               .format(str(round(length*kappa_,2)))]
     plot_titles.append(plot_titles[-1])
 
-    line_labels: List[Optional[str]] = ["port 0", "port 1", "port 2", "port 3"]
+    line_labels: List[Optional[str]] = ["port 0", "port 1", "port 2", "port 3",
+                                        "init noise port 0",
+                                        "final noise port 2",
+                                        "init noise port 1",
+                                        "final noise port 3"]
 
     oc.plot2d(x_datas, y_datas, plot_groups=plot_groups,
               plot_titles=plot_titles, x_labels=['t'], y_labels=['P_t'],

optcom/effects/asymmetry.py

@@ -54,12 +54,17 @@ def __init__(self, delta: Union[float, Callable] = None) -> None:
         super().__init__()
         # The asymmetry coefficient ------------------------------------
         self._delta_op: np.ndarray = np.array([])
-        self._delta: Union[np.ndarray, Callable]
+        self._delta: Callable
         if (callable(delta)):
             self._delta = delta
         else:
             self._delta = lambda omega: np.ones_like(omega) * delta
     # ==================================================================
+    @property
+    def delta(self) -> Callable:
+
+        return self._delta
+    # ==================================================================
     def set(self, center_omega: np.ndarray = np.array([]),
             abs_omega: np.ndarray = np.array([])) -> None:
 

optcom/effects/coupling.py

@@ -70,7 +70,7 @@ def __init__(self, kappa: Union[List[float], Callable],
         # The coupling coefficient -------------------------------------
         self._op: np.ndarray = np.array([])
         self._kappa_op: np.ndarray = np.array([])
-        self._kappa: Union[np.ndarray, Callable]
+        self._kappa: Callable
         if (callable(kappa)):
             self._kappa = kappa
         else:
@@ -79,6 +79,11 @@ def __init__(self, kappa: Union[List[float], Callable],
             self._order_taylor = len(kappa_) - 1
     # ==================================================================
     @property
+    def kappa(self) -> Callable:
+
+        return self._kappa            
+    # ==================================================================
+    @property
     def order_taylor(self) -> int:
 
         return self._order_taylor

optcom/equations/abstract_ampnlse.py

@@ -127,15 +127,15 @@ def __init__(self, re_fiber: AbstractREFiber, gain_order: int,
         self._gain_coeff: List[DopedFiberGain] = []
         self._absorp_coeff: List[DopedFiberGain] = []
         k: int = 0  # counter
-        for i in range(4):
+        for i in range(self._nbr_eqs):
             k = i//2
             self._gain_coeff.append(DopedFiberGain(self._sigma_e[k],
                                                    self._overlap[k], 0.))
             self._absorp_coeff.append(DopedFiberGain(self._sigma_a[k],
                                                      self._overlap[k], 0.))
         self._gains: List[Gain] = []
         self._absorps: List[Attenuation] = []
-        for i in range(4):
+        for i in range(self._nbr_eqs):
             k = i//2
             self._gains.append(Gain(self._gain_coeff[i], gain_order,
                                     start_taylor=start_taylor_gain,
@@ -152,7 +152,7 @@ def __init__(self, re_fiber: AbstractREFiber, gain_order: int,
             self._sat_gains = []
             start_taylor_gain = 1
             en_sat_: Union[float, Callable]
-            for i in range(4):
+            for i in range(self._nbr_eqs):
                 k = i//2
                 crt_en_sat = en_sat[k]
                 if (crt_en_sat is None):

optcom/equations/abstract_cnlse.py

@@ -15,7 +15,7 @@
 """.. moduleauthor:: Sacha Medaer"""
 
 import copy
-from typing import Callable, List, Optional, Union
+from typing import Callable, List, Optional, overload, Union
 
 import numpy as np
 
@@ -140,11 +140,11 @@ def __init__(self, nbr_fibers: int,
                          INTER_PORT_DELAY=INTER_PORT_DELAY)
         kerr_effects_ = util.make_list(kerr_effects, nbr_fibers)
         raman_effects_ = util.make_list(raman_effects, nbr_fibers)
-        self._asym: List[List[AbstractEffect]] =\
+        self._asym: List[List[Asymmetry]] =\
             [[] for i in range(nbr_fibers)]
-        self._coup: List[List[AbstractEffect]] =\
+        self._coup: List[List[Coupling]] =\
             [[] for i in range(nbr_fibers)]
-        self._raman: List[List[AbstractEffect]] =\
+        self._raman: List[List[Union[Raman, RamanApprox]]] =\
             [[] for i in range(nbr_fibers)]
         kappa_: List[List[Union[List[float], Callable]]] =\
             [[[] for i in range(nbr_fibers)] for i in range(nbr_fibers)]
@@ -197,12 +197,24 @@ def __init__(self, nbr_fibers: int,
                         self._coup[i].append(Coupling(kappa_[i][j]))
                         self._add_ind_effect(self._coup[i][-1], i, j)
     # ==================================================================
-    def __call__(self, waves: np.ndarray, z: float, h: float) -> np.ndarray:
-        res = np.zeros_like(waves, dtype=cst.NPFT)
-        for i in range(len(waves)):
-            res[i] = self.term_ind(waves, i, z)
+    @overload
+    def __call__(self, waves: np.ndarray, z: float, h: float
+                 )-> np.ndarray: ...
+    # ------------------------------------------------------------------
+    @overload
+    def __call__(self, waves: np.ndarray, z: float, h: float, ind: int
+                 ) -> np.ndarray: ...
+    # ------------------------------------------------------------------
+    def __call__(self, *args):
+        if (len(args) == 4):
+            waves, z, h, ind = args
+            res = np.zeros_like(waves[ind], dtype=cst.NPFT)
+            res = self.term_ind(waves, ind, z)
 
-        return res
+            return res
+        else:
+
+            raise NotImplementedError()
     # ==================================================================
     def open(self, domain: Domain, *fields: List[Field]) -> None:
         super().open(domain, *fields)
@@ -216,6 +228,20 @@ def open(self, domain: Domain, *fields: List[Field]) -> None:
             for raman in self._raman[i]:
                 raman.set()
     # ==================================================================
+    def get_kappa_for_noise(self):
+        # temp harcoding for noise manangement
+        return self._coup_coeff[0]  # sym assumption, both are equal anyway
+    # ==================================================================
+    def close(self, domain: Domain, *fields: List[Field]) -> None:
+        # tem hardcoding -> must change
+        self._coup_coeff: List[Callable] = [] # len == 2
+        if (self._nbr_eqs == 2):    # hardcoding for 2 -> need to generalize
+            for i in range(len(self._coup)): # len == 2
+                for coup in self._coup[i]:  # len == 1
+                    self._coup_coeff.append(coup.kappa(self._noise_omega)[0])
+        # Call to part
+        super().close(domain, *fields)
+    # ==================================================================
     @sync_waves_decorator
     def op_non_lin(self, waves: np.ndarray, id: int,
                    corr_wave: Optional[np.ndarray] = None

optcom/equations/abstract_equation.py

@@ -15,7 +15,7 @@
 """.. moduleauthor:: Sacha Medaer"""
 
 from abc import ABCMeta
-from typing import List, Optional, Tuple
+from typing import List, Optional, overload, Tuple
 
 import numpy as np
 
@@ -29,5 +29,27 @@ def __init__(self) -> None:
 
         return None
     # ==================================================================
+    @overload
     def __call__(self, vectors: np.ndarray, z: float, h: float
+                 )-> np.ndarray: ...
+    # ------------------------------------------------------------------
+    @overload
+    def __call__(self, vectors: np.ndarray, z: float, h: float, ind: int
                  ) -> np.ndarray: ...
+    # ------------------------------------------------------------------
+    def __call__(self, *args):
+        """
+        Parameters
+        ----------
+        vectors :
+            The value of the unknowns at the current step.
+        z :
+            The current value of the variable.
+        h :
+            The step size.
+        ind :
+            The index of the considered vector.
+
+        """
+
+        raise NotImplementedError()

optcom/equations/abstract_field_equation.py

@@ -16,7 +16,8 @@
 
 from __future__ import annotations
 
-from typing import Callable, List, Optional, Tuple, TYPE_CHECKING, Union
+from typing import Callable, List, Optional, overload, Tuple, TYPE_CHECKING,\
+                   Union
 
 import numpy as np
 
@@ -520,19 +521,40 @@ def post_pass(self) -> None:
 
         return None
     # ==================================================================
+    @overload
     def calc_noise(self, noises: np.ndarray, z: float, h: float
-                   ) -> np.ndarray:
+                 )-> np.ndarray: ...
+    # ------------------------------------------------------------------
+    @overload
+    def calc_noise(self, noises: np.ndarray, z: float, h: float, ind: int
+                 ) -> np.ndarray: ...
+    # ------------------------------------------------------------------
+    def calc_noise(self, *args):
         """Compute noise."""
-        if (self._NOISE):
+        if (len(args) == 4):
+            noises, z, h, ind = args
+            if (self._NOISE):
+                noises_ = np.zeros_like(noises[ind])
+                eq_id = self.id_tracker.eq_id_of_field_id(ind)
+                noises_ = self._noise_func_per_field[eq_id](noises, z, h, ind)
+
+                return noises_
+            else:
+
+                return noises[ind]
+        elif (len(args) == 3):
+            # For now, dummy identity function for equations that do not
+            # implement any noise, might want to find another solution
+            noises, z, h = args
             noises_ = np.zeros_like(noises)
-            nbr_noises = self.id_tracker.nbr_fields
-            for i in range(nbr_noises):
+            for i in range(len(noises)):
                 eq_id = self.id_tracker.eq_id_of_field_id(i)
-                noises_[i] = self._noise_func_per_field[eq_id](noises[i], z, h)
+                noises_[i] = self._noise_func_per_field[eq_id](noises, z, h, i)
 
             return noises_
+        else:
 
-        return noises
+            raise NotImplementedError()
     # ==================================================================
     def _initiate_noise_managament(self) -> None:
         """Initiate noise manangement. Link each noise array to the
@@ -557,7 +579,8 @@ def _initiate_noise_managament(self) -> None:
                 crt_ops = ['+' for _ in range(len(crt_noise_effects)-1)]
                 noise_fct = CallableLittExpr(effects_to_add, crt_ops)
             else:
-                noise_fct = CallableLittExpr([lambda identity, z, h: identity])
+                noise_fct = CallableLittExpr([lambda noises, z, h, ind:
+                                              noises[ind]])
             self._noise_func_per_field.append(noise_fct)
     # ==================================================================
     def clear(self) -> None:

optcom/equations/abstract_nlse.py

@@ -268,6 +268,6 @@ def open(self, domain: Domain, *fields: List[Field]) -> None:
         # Noise --------------------------------------------------------
         if (self._att is not None):
             att_coeff: np.ndarray = self._att.coeff(self._noise_omega, 0)
-            noise_fct = lambda noise, z, h: (-1*noise*att_coeff)
+            noise_fct = lambda noises, z, h, ind: (-1*noises[ind]*att_coeff)
             for i in range(self._nbr_eqs):
                 self._add_noise_effect(noise_fct, i)