Merge pull request #27 from Digiratory/23-kasdin-generator

Kasdin generator

Author: Sinitca-Aleksandr

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## 1.7.1
 
+* [GH-23](https://github.com/Digiratory/StatTools/issues/23) feat: add Kasdin generator. fix: change first arg in lfilter in LBFBm generator.
 * [GH-15](https://github.com/Digiratory/StatTools/issues/15) feat&fix: LBFBm generator update: generate with input value and return an increment instead of the absolute value of the signal.
 
 ## 1.7.0

README.md

@@ -30,7 +30,7 @@ from StatTools.filters import FilteredArray
 
 h = 0.8                 # choose Hurst parameter
 total_vectors = 1000    # total number of vectors in output
-vectors_length = 1440   # each vector's length 
+vectors_length = 1440   # each vector's length
 t = 8                   # threads in use during computation
 
 correlated_vectors = Filter(h, vectors_length).generate(n_vectors=total_vectors,
@@ -49,13 +49,32 @@ base = 1.2          # the basis for the filter optimization algorithm
 target_len = 4000   # number of generation iterations
 
 generator = LBFBmGenerator(h, filter_len, base)
-trajectory = []
-    for value in islice(generator, target_len):
-        trajectory.append(value)
+signal = []
+for value in islice(generator, target_len):
+    signal.append(value)
 ```
 
 For more information and generator validation, see [lbfbm_generator.ipynb](/research/lbfbm_generator.ipynb).
 
+It is also possible to use the method of generating increments with a given H using `KasdinGenerator`.
+
+```python
+from StatTools.generators.kasdin_generator import KasdinGenerator
+h = 0.8             # choose Hurst parameter
+target_len = 4000   # number of generation iterations
+
+generator = KasdinGenerator(h, length=target_len)
+
+# the first option
+signal = generator.get_full_sequence()
+
+# the second option
+signal_list = []
+for sample in generator:
+    signal_list.append(sample)
+```
+For more information see Kasdin, N. J. (1995). Discrete simulation of colored noise and stochastic processes and 1/f/sup /spl alpha// power law noise generation. doi:10.1109/5.381848.
+
 ### Fluctuational Analysis
 
 1. Example of Detrended Fluctuational Analysis (DFA)

StatTools/generators/kasdin_generator.py

@@ -0,0 +1,71 @@
+from itertools import islice
+from typing import Iterator, Optional
+
+import numpy as np
+from scipy.signal import lfilter
+
+
+class KasdinGenerator:
+    """
+    Generates a sequence of numbers according to the Kasdin model.
+    Based on the method proposed in the article Kasdin, N. J. (1995).
+        Discrete simulation of colored noise and stochastic processes and 1/f/sup /spl alpha// power law noise generation.
+        doi:10.1109/5.381848
+
+    Args:
+        h (float): Hurst exponent (0 < H < 2)
+        length (int): Maximum length of the sequence.
+        random_generator (Iterator[float], optional): Iterator providing random values.
+            Defaults is iter(np.random.randn(), None).
+    Raises:
+        ValueError: If length is less than 1
+        StopIteration('Sequence exhausted') : If maximum sequence length has been reached.
+
+    Example usage:
+    >>> generator = KasdinGenerator(h, length)
+    >>> trj = list(generator)
+    """
+
+    def __init__(
+        self,
+        h: float,
+        length: int,
+        random_generator: Optional[Iterator[float]] = iter(np.random.randn, None),
+    ) -> None:
+        if length is not None and length < 1:
+            raise ValueError("Length must be more than 1")
+        self.h = h
+        self.length = length
+        self.random_generator = random_generator
+
+        # init filter coefficients
+        beta = 2 * self.h - 1
+        self.filter_coefficients = np.zeros(self.length, dtype=np.float64)
+        self.filter_coefficients[0] = 1.0
+        k = np.arange(1, self.length)
+        self.filter_coefficients[1:] = np.cumprod((k - 1 - beta / 2) / k)
+
+        # generate the sequence
+        random_sequence = np.fromiter(
+            islice(random_generator, self.length), dtype=np.float64
+        )
+        self.sequence = lfilter(1, self.filter_coefficients, random_sequence)
+        self.current_index = 0
+
+    def get_filter_coefficients(self):
+        """Returns the filter coefficients."""
+        return self.filter_coefficients
+
+    def __iter__(self) -> "KasdinGenerator":
+        return self
+
+    def __next__(self) -> float:
+        """Return next value in sequence"""
+        if self.current_index >= self.length:
+            raise StopIteration("Sequence exhausted")
+        self.current_index += 1
+        return self.sequence[self.current_index - 1]
+
+    def get_full_sequence(self) -> np.ndarray:
+        """Return full generated sequence."""
+        return self.sequence

StatTools/generators/lbfbm_generator.py

@@ -1,9 +1,10 @@
 import math
 import warnings
-from typing import List, Iterator, Optional
+from typing import Iterator, List, Optional
+
 import numpy as np
-from scipy.signal import lfilter
 from numpy.typing import NDArray
+from scipy.signal import lfilter
 
 
 def signed_power(base: float, degree: float) -> float:
@@ -176,8 +177,7 @@ def _find_filter_len(self, base, length):
 
     def _calculate_step(self) -> float:
         """Applies a filter."""
-        res = lfilter(np.ones(self.filter_len), self.matrix_a, self.bins[::-1])
-        return res[-1] - res[-2]
+        return lfilter(1, self.matrix_a, self.bins[::-1])[-1]
 
     def __iter__(self) -> "LBFBmGenerator":
         return self

tests/test_filter.py

@@ -1,8 +1,8 @@
 import numpy as np
 import pytest
 
-from StatTools.generators.base_filter import Filter
 from StatTools.analysis.dfa import DFA
+from StatTools.generators.base_filter import Filter
 
 testdata = {
     "target_mean": [0.5, 0.7, 0.9],
@@ -24,21 +24,25 @@ def test_filter_generator(h, length, target_std, target_mean):
     Test that the generated data has the specified mean and standard deviation.
     """
     generator = Filter(h, length, set_mean=target_mean, set_std=target_std)
-    trajectory = list(generator)
-
-    actual_mean = np.mean(trajectory)
-    actual_std = np.std(trajectory, ddof=1)
-    actual_h = DFA(trajectory).find_h()
+    mean_difference_mean = 0
+    mean_difference_std = 0
+    mean_difference_h = 0
+    times = 5
+    for _ in range(times):
+        trajectory = list(generator)
+        mean_difference_mean += np.mean(trajectory)
+        mean_difference_std += np.std(trajectory, ddof=1)
+        mean_difference_h += DFA(trajectory).find_h()
 
     assert (
-        abs(actual_mean - target_mean) < 0.001
-    ), f"Mean deviation too large: expected {target_mean}, got {actual_mean}"
+        abs(mean_difference_mean / times - target_mean) < 0.001
+    ), f"Mean deviation too large: expected {target_mean}, got {mean_difference_mean}"
     assert (
-        abs(actual_std - target_std) < 0.001
-    ), f"Std deviation too large: expected {target_std}, got {actual_std}"
-    assert abs(actual_h - h) < (
+        abs(mean_difference_std / times - target_std) < 0.001
+    ), f"Std deviation too large: expected {target_std}, got {mean_difference_std}"
+    assert abs(mean_difference_h / times - h) < (
         h * 0.15
-    ), f"Hurst deviation too large: expected {h}, got {actual_h}"
+    ), f"Hurst deviation too large: expected {h}, got {mean_difference_h}"
 
 
 @pytest.mark.parametrize("h", testdata["h"])
@@ -50,21 +54,25 @@ def test_filter(h, length, target_std, target_mean):
     Test that the generated data has the specified mean and standard deviation.
     """
     generator = Filter(h, length, set_mean=target_mean, set_std=target_std)
-    trajectory = generator.generate(n_vectors=1)
-
-    actual_mean = np.mean(trajectory)
-    actual_std = np.std(trajectory, ddof=1)
-    actual_h = DFA(trajectory).find_h()
+    mean_difference_mean = 0
+    mean_difference_std = 0
+    mean_difference_h = 0
+    times = 3
+    for _ in range(times):
+        trajectory = generator.generate(n_vectors=1)
+        mean_difference_mean += np.mean(trajectory)
+        mean_difference_std += np.std(trajectory, ddof=1)
+        mean_difference_h += DFA(trajectory).find_h()
 
     assert (
-        abs(actual_mean - target_mean) < 0.001
-    ), f"Mean deviation too large: expected {target_mean}, got {actual_mean}"
+        abs(mean_difference_mean / times - target_mean) < 0.001
+    ), f"Mean deviation too large: expected {target_mean}, got {mean_difference_mean}"
     assert (
-        abs(actual_std - target_std) < 0.001
-    ), f"Std deviation too large: expected {target_std}, got {actual_std}"
-    assert abs(actual_h - h) < (
+        abs(mean_difference_std / times - target_std) < 0.001
+    ), f"Std deviation too large: expected {target_std}, got {mean_difference_std}"
+    assert abs(mean_difference_h / times - h) < (
         h * 0.15
-    ), f"Hurst deviation too large: expected {h}, got {actual_h}"
+    ), f"Hurst deviation too large: expected {h}, got {mean_difference_h}"
 
 
 @pytest.mark.parametrize("h", testdata["h"])
@@ -75,22 +83,27 @@ def test_filter_2d(h, length, target_std, target_mean):
     """
     Test that the generated data has the specified mean and standard deviation.
     """
+    count = 3
     generator = Filter(h, length, set_mean=target_mean, set_std=target_std)
-    trajectories = generator.generate(n_vectors=3)
+    trajectories = generator.generate(n_vectors=count)
 
-    for i in range(3):
+    mean_difference_mean = 0
+    mean_difference_std = 0
+    mean_difference_h = 0
+
+    for i in range(count):
         trajectory = trajectories[i]
 
-        actual_mean = np.mean(trajectory)
-        actual_std = np.std(trajectory, ddof=1)
-        actual_h = DFA(trajectory).find_h()
+        mean_difference_mean += np.mean(trajectory)
+        mean_difference_std += np.std(trajectory, ddof=1)
+        mean_difference_h += DFA(trajectory).find_h()
 
-        assert (
-            abs(actual_mean - target_mean) < 0.001
-        ), f"Mean deviation too large: expected {target_mean}, got {actual_mean}"
-        assert (
-            abs(actual_std - target_std) < 0.001
-        ), f"Std deviation too large: expected {target_std}, got {actual_std}"
-        assert abs(actual_h - h) < (
-            h * 0.15
-        ), f"Hurst deviation too large: expected {h}, got {actual_h}"
+    assert (
+        abs(mean_difference_mean / count - target_mean) < 0.001
+    ), f"Mean deviation too large: expected {target_mean}, got {mean_difference_mean}"
+    assert (
+        abs(mean_difference_std / count - target_std) < 0.001
+    ), f"Std deviation too large: expected {target_std}, got {mean_difference_std}"
+    assert abs(mean_difference_h / count - h) < (
+        h * 0.15
+    ), f"Hurst deviation too large: expected {h}, got {mean_difference_h}"

tests/test_kasdin_generator.py

@@ -0,0 +1,56 @@
+import numpy as np
+import pytest
+
+from StatTools.analysis.dfa import DFA
+from StatTools.generators.kasdin_generator import KasdinGenerator
+
+testdata = {
+    "h_list": [i * 0.01 for i in range(50, 200, 20)],
+    "rate_list": [14],
+}
+
+
+def get_test_h(
+    h: float,
+    target_len: int,
+) -> float:
+    """
+    Calculates the Hurst exponent for the generated trajectory.
+
+    Parameters:
+        base: The base of the number system for bins
+        filter_len: Filter length
+        h: The specified Hurst exponent
+        scales: Scales for analysis
+        step: The step for analysis
+
+    Returns:
+        Calculated Hurst exponent (h_gen)
+    """
+    generator = KasdinGenerator(h, length=target_len)
+    signal = generator.get_full_sequence()
+    dfa = DFA(signal)
+    return dfa.find_h()
+
+
+@pytest.mark.parametrize("h", testdata["h_list"])
+@pytest.mark.parametrize("rate", testdata["rate_list"])
+def test_kasdin_generator(h: float, rate: int):
+    """
+    It tests the generator for compliance with the specified Hurst exponent.
+
+    Parameters:
+        h: The specified Hurst exponent
+        base: The base of the number system for bins
+    """
+    threshold = 0.10
+    times = 3
+    mean_difference = 0
+    length = 2**rate
+    for _ in range(times):
+        h_gen = get_test_h(h, length)
+        mean_difference += abs(h_gen - h) / h
+    mean_difference /= times
+    assert (
+        mean_difference <= threshold
+    ), f"Diff between h and h_gen exceeds {threshold * 100}%: h={h}, h_gen={h_gen}, mean diff={mean_difference * 100:.2f}%"

tests/test_lbfbm_generator.py

@@ -1,6 +1,6 @@
-from scipy import stats
 import numpy as np
 import pytest
+from scipy import stats
 
 from StatTools.analysis.dpcca import dpcca
 from StatTools.generators.lbfbm_generator import LBFBmGenerator, normalize
@@ -79,7 +79,7 @@ def test_lbfbm_generator(h: float, base: float, rate: int):
         base: The base of the number system for bins
     """
     threshold = 0.10
-    times = 10
+    times = 3
     mean_difference = 0
     length = 2**rate
     scales = np.array([2**i for i in range(3, rate - 3)])