Fix inconsistency with cpp impl, improve design (#3)

 - fix: RNG is skipped when volatility is 0
 - RNG as an argument to higher level generators
 - Separate computation of available units and available power
 - add unit tests

Signed-off-by: Sylvain Leclerc <[email protected]>

requirements.txt

@@ -1,3 +1,4 @@
 numpy==1.24.4
 pandas
+pydantic
 PyYAML

src/cluster_import.py

@@ -14,11 +14,14 @@
 
 import numpy as np
 
-from ts_generator import ProbilityLaw, ThermalCluster
+from ts_generator import ProbabilityLaw, ThermalCluster
 
 
 def import_thermal_cluster(path: Path, days_per_year: int = 365) -> ThermalCluster:
-    law_dict = {"UNIFORM": ProbilityLaw.UNIFORM, "GEOMETRIC": ProbilityLaw.GEOMETRIC}
+    law_dict = {
+        "UNIFORM": ProbabilityLaw.UNIFORM,
+        "GEOMETRIC": ProbabilityLaw.GEOMETRIC,
+    }
     array = np.genfromtxt(path, delimiter=",", dtype=str)
     return ThermalCluster(
         unit_count=int(array[1][1]),
@@ -28,10 +31,10 @@ def import_thermal_cluster(path: Path, days_per_year: int = 365) -> ThermalClust
         fo_volatility=float(array[5][1]),
         po_law=law_dict[array[6][1]],
         po_volatility=float(array[7][1]),
-        fo_duration=array[8][1 : days_per_year + 1].astype(float),
+        fo_duration=array[8][1 : days_per_year + 1].astype(int),
         fo_rate=array[9][1 : days_per_year + 1].astype(float),
-        po_duration=array[10][1 : days_per_year + 1].astype(float),
+        po_duration=array[10][1 : days_per_year + 1].astype(int),
         po_rate=array[11][1 : days_per_year + 1].astype(float),
-        npo_min=array[12][1 : days_per_year + 1].astype(float),
-        npo_max=array[13][1 : days_per_year + 1].astype(float),
+        npo_min=array[12][1 : days_per_year + 1].astype(int),
+        npo_max=array[13][1 : days_per_year + 1].astype(int),
     )

src/cluster_resolve.py

@@ -13,15 +13,18 @@
 import pandas as pd
 
 from cluster_parsing import InputCluster
-from ts_generator import ProbilityLaw, ThermalCluster
+from ts_generator import ProbabilityLaw, ThermalCluster
 
 
 def resolve_thermal_cluster(
     parsed_yaml: InputCluster,
     parameters_ts: pd.core.frame.DataFrame,
     modulation: pd.core.frame.DataFrame,
 ) -> ThermalCluster:
-    law_dict = {"UNIFORM": ProbilityLaw.UNIFORM, "GEOMETRIC": ProbilityLaw.GEOMETRIC}
+    law_dict = {
+        "UNIFORM": ProbabilityLaw.UNIFORM,
+        "GEOMETRIC": ProbabilityLaw.GEOMETRIC,
+    }
     return ThermalCluster(
         unit_count=parsed_yaml.unit_count,
         nominal_power=parsed_yaml.nominal_power,

src/duration_generator.py

@@ -0,0 +1,115 @@
+# Copyright (c) 2024, RTE (https://www.rte-france.com)
+#
+# See AUTHORS.txt
+#
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#
+# SPDX-License-Identifier: MPL-2.0
+#
+# This file is part of the Antares project.
+from abc import ABC, abstractmethod
+from enum import Enum
+from math import log, sqrt
+from typing import List
+
+import numpy as np
+
+from random_generator import RNG
+
+
+class ProbabilityLaw(Enum):
+    UNIFORM = "UNIFORM"
+    GEOMETRIC = "GEOMETRIC"
+
+
+class DurationGenerator(ABC):
+    """
+    Logic to generate unavailability duration.
+    """
+
+    @abstractmethod
+    def generate_duration(self, day: int) -> int:
+        ...
+
+
+class GeneratorWrapper(DurationGenerator):
+    """
+    Wraps another generator to skip unnecessary random number generations.
+    Used to keep backward compat with cpp implementation.
+    """
+
+    def __init__(
+        self, delegate: DurationGenerator, volatility: float, expecs: List[int]
+    ) -> None:
+        self.volatility = volatility
+        self.expectations = expecs
+        self.delegate = delegate
+
+    def generate_duration(self, day: int) -> int:
+        """
+        generation of random outage duration
+        """
+        expectation = self.expectations[day]
+        # Logic copied from cpp implementation for results preservation
+        if self.volatility == 0 or expectation == 1:
+            return expectation
+        return self.delegate.generate_duration(day)
+
+
+class UniformDurationGenerator(DurationGenerator):
+    def __init__(self, rng: RNG, volatility: float, expecs: List[int]) -> None:
+        self.rng = rng
+        self.a = np.empty(len(expecs), dtype=float)
+        self.b = np.empty(len(expecs), dtype=float)
+        for day, expec in enumerate(expecs):
+            xtemp = volatility * (expec - 1)
+            self.a[day] = expec - xtemp
+            self.b[day] = 2 * xtemp + 1
+
+    def generate_duration(self, day: int) -> int:
+        """
+        generation of random outage duration
+        """
+        rnd_nb = self.rng.next()
+        return int(self.a[day] + rnd_nb * self.b[day])
+
+
+class GeometricDurationGenerator(DurationGenerator):
+    def __init__(self, rng: RNG, volatility: float, expecs: List[int]) -> None:
+        self.rng = rng
+        self.a = np.empty(len(expecs), dtype=float)
+        self.b = np.empty(len(expecs), dtype=float)
+        for day, expec in enumerate(expecs):
+            xtemp = volatility * volatility * expec * (expec - 1)
+            if xtemp != 0:
+                ytemp = (sqrt(4 * xtemp + 1) - 1) / (2 * xtemp)
+                self.a[day] = expec - 1 / ytemp
+                self.b[day] = 1 / log(1 - ytemp)
+            else:
+                self.a[day] = expec - 1
+                self.b[day] = 0
+
+    def generate_duration(self, day: int) -> int:
+        """
+        generation of random outage duration
+        """
+        rnd_nb = self.rng.next()
+        return min(int(1 + self.a[day] + self.b[day] * log(rnd_nb)), 1999)
+
+
+def make_duration_generator(
+    rng: RNG, law: ProbabilityLaw, volatility: float, expectations: List[int]
+) -> DurationGenerator:
+    """
+    return a DurationGenerator for the given law
+    """
+    base_rng: DurationGenerator
+    if law == ProbabilityLaw.UNIFORM:
+        base_rng = UniformDurationGenerator(rng, volatility, expectations)
+    elif law == ProbabilityLaw.GEOMETRIC:
+        base_rng = GeometricDurationGenerator(rng, volatility, expectations)
+    else:
+        raise ValueError(f"Unknown law type: {law}")
+    return GeneratorWrapper(base_rng, volatility, expectations)

src/mersenne_twister.py

@@ -27,10 +27,12 @@
 # LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
+import dataclasses
+from dataclasses import field
 from typing import List, Tuple
 
 
+@dataclasses.dataclass
 class MersenneTwister:
     periodN: int = 624
     periodM: int = 397
@@ -40,11 +42,10 @@ class MersenneTwister:
 
     MAG: Tuple[int, int] = (0, MATRIX_A)
 
-    mt: List[int] = [0] * periodN
+    mt: List[int] = field(default_factory=lambda: [0] * 624)
 
     mti: int = 0
 
-    @classmethod
     def seed(self, seed: int) -> None:
         self.mt[0] = seed & 0xFFFFFFFF
         for i in range(1, self.periodN):
@@ -53,7 +54,6 @@ def seed(self, seed: int) -> None:
 
         self.mti = self.periodN
 
-    @classmethod
     def next(self) -> float:
         if self.mti == self.periodN:
             for j in range(self.periodN - self.periodM):
@@ -87,6 +87,5 @@ def next(self) -> float:
 
         return y / 4294967295
 
-    @classmethod
     def reset(self) -> None:
         self.seed(5489)

src/random_generator.py

@@ -0,0 +1,48 @@
+# Copyright (c) 2024, RTE (https://www.rte-france.com)
+#
+# See AUTHORS.txt
+#
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+#
+# SPDX-License-Identifier: MPL-2.0
+#
+# This file is part of the Antares project.
+import random
+from abc import ABC, abstractmethod
+from typing import Optional
+
+from mersenne_twister import MersenneTwister
+
+
+class RNG(ABC):
+    """
+    Random number generator interface
+    """
+
+    @abstractmethod
+    def next(self) -> float:
+        ...
+
+
+class PythonRNG(ABC):
+    """
+    Native python RNG.
+    """
+
+    def next(self) -> float:
+        return random.random()
+
+
+class MersenneTwisterRNG(RNG):
+    """
+    Our own RNG based on Mersenne-Twister algorithm.
+    """
+
+    def __init__(self, seed: int = 5489):
+        self._rng = MersenneTwister()
+        self._rng.seed(seed)
+
+    def next(self) -> float:
+        return self._rng.next()