all dynamic models for prelim experiment 1

interfere/dynamics/coupled_map_lattice.py

@@ -184,7 +184,7 @@ def __init__(
         sigma: float = 0.0,
         x_min: float = 0.0,
         x_max: float = 1.0,
-        boundry_condition: Optional[str] = "none",
+        boundary_condition: Optional[str] = "none",
         tsteps_btw_obs: int = 1,
         measurement_noise_std: Optional[np.ndarray] = None
     ):
@@ -223,7 +223,7 @@ def __init__(
             x_max (float): Optional maximum bound (applied to state
                 elementwise) to ensure that the noise does not peturb the system
                 out of it's domain.
-            boundry_condition (string): One of ["none", "truncate", "mod"]. If
+            boundary_condition (string): One of ["none", "truncate", "mod"]. If
                 "truncate", is selected, entries in `x` are reset to be the
                 closer of `x_min` or `x_max` when they leave the domain. If
                 "mod" is selected, and and entry in `x`, `x_i` is not in
@@ -242,7 +242,7 @@ def __init__(
         self.sigma = sigma
         self.x_max = x_max
         self.x_min = x_min
-        self.boundry_condition = boundry_condition
+        self.boundary_condition = boundary_condition
         super().__init__(
             adjacency_matrix, eps, f, f_params, tsteps_btw_obs, measurement_noise_std)
 
@@ -266,10 +266,10 @@ def step(self, x: np.ndarray, t: float, rng: np.random.mtrand.RandomState):
             x_next += rng.normal(0, self.sigma, size=self.dim)
 
             # See if the state is within boundary and apply appropriate transform. 
-            if self.boundry_condition == "mod":
+            if self.boundary_condition == "mod":
                 x_next = mod_interval(x_next, self.x_min, self.x_max)
 
-            elif self.boundry_condition == "truncate":
+            elif self.boundary_condition == "truncate":
                 x_next[x_next < self.x_min] = self.x_min
                 x_next[x_next > self.x_max] = self.x_max
 
@@ -280,6 +280,7 @@ def coupled_logistic_map(
     adjacency_matrix: np.array,
     logistic_param=3.72,
     eps=0.5,
+    sigma=0.0,
     measurement_noise_std: Optional[np.ndarray] = None
 ):
     """Initializes an N-dimensional coupled logistic map model.
@@ -299,6 +300,8 @@ def coupled_logistic_map(
             r, where the map is f(x) = rx(1-x).
         eps (float): A parameter that controls the relative strenth of    
             coupling. High epsilon means greater connection to neighbors.
+        sigma (float): The standard deviation of the additive gaussian noise
+            in the model.
         measurement_noise_std (ndarray): None, or a vector with shape (n,)
             where each entry corresponds to the standard deviation of the
             measurement noise for that particular dimension of the dynamic
@@ -307,8 +310,15 @@ def coupled_logistic_map(
             independent gaussian noise with standard deviation 1 and 10
             will be added to x1 and x2 respectively at each point in time.
     """
-    return CoupledMapLattice(
-        adjacency_matrix, eps, logistic_map, (logistic_param,), measurement_noise_std=measurement_noise_std)
+    return StochasticCoupledMapLattice(
+        adjacency_matrix,
+        eps,
+        logistic_map,
+        (logistic_param,),
+        sigma=sigma,
+        measurement_noise_std=measurement_noise_std, 
+        boundary_condition="truncate"
+    )
 
 
 
@@ -357,7 +367,7 @@ def stochastic_coupled_map_1dlattice(
         sigma=sigma,
         x_min=-1,
         x_max=1,
-        boundry_condition="truncate",
+        boundary_condition="truncate",
         tsteps_btw_obs=tsteps_btw_obs, measurement_noise_std=measurement_noise_std
     )
 

interfere/dynamics/kuramoto.py

@@ -11,6 +11,7 @@
 neighbors.
 """
 from typing import Optional, Callable
+from warnings import warn
 
 import numpy as np
 from pyclustering.nnet.fsync import fsync_network
@@ -104,8 +105,11 @@ def simulate(
     ) -> np.ndarray:        
         # Check initial condition.
         if np.any(np.abs(initial_condition) > 1):
-            raise ValueError("Kuramoto Models require initial conditions in "
-                             "the interval (-1, 1).")
+            warn("Kuramoto Models require initial conditions in "
+                 "the interval (-1, 1). Initial conditions will be thresholded at -1 and 1.")
+            initial_condition[initial_condition > 1] = 1
+            initial_condition[initial_condition < -1] = -1
+
         # Extract phase of the initial condition.
         theta0 = np.arcsin(initial_condition)
 

interfere/dynamics/lotka_voltera.py

@@ -24,6 +24,12 @@ def __init__(
         where r_i and k_i are the growth rates and carrying capacities of
         species i, A is the matrix of interspecies interactions.
 
+        See:
+        * Vadillo 2019, "Comparing stochastic Lotka–Volterra predator-prey
+        models"
+        * https://en.wikipedia.org/wiki/Competitive_Lotka%E2%80%93Volterra_equations        
+        * https://github.com/netsiphd/netrd/blob/master/netrd/dynamics/lotka_volterra.py
+
         Args:
             growth_rates (ndarray): A length n vector of growth rates (r_i's).
             capacities (ndarray): A length n vector of carrying capacities 
@@ -52,16 +58,17 @@ def __init__(
                             )
 
         # Assign parameters.
+        dim = len(growth_rates)
         self.growth_rates = growth_rates
         self.capacities = capacities
         self.interaction_mat = interaction_mat
         # Set dimension of the system.
-        super().__init__(len(growth_rates), measurement_noise_std)
+        super().__init__(dim, measurement_noise_std)
 
     def dXdt(self, x: np.ndarray, t: Optional[float] = None):
         """Coputes derivative of a generalized Lotka Voltera model.
 
-        dx_i/dt = r_i * x_i * (1 - x_i / k_i +  [A x]_i / k_i)
+        dx_i/dt = r_i * x_i * (1 - (x_i +  [A x]_i) / k_i)
 
         Args:
             x (ndarray): The current state of the system.
@@ -72,9 +79,7 @@ def dXdt(self, x: np.ndarray, t: Optional[float] = None):
             The derivative of the system at x and t with respect to time.
         """
         return self.growth_rates * x * (
-            1 - x / self.capacities + self.interaction_mat @ (
-                x / self.capacities
-            )
+            1 -  (x + self.interaction_mat @ x) / self.capacities
         )
 
 
@@ -90,12 +95,18 @@ def __init__(
     ):
         """Initializes class for simulating Lotka Voltera dynamics.
 
-            dx_i/dt = r_i * x_i * (1 - x_i / k_i +  [A x]_i / k_i) + sigma * dW
+            dx_i/dt = r_i * x_i * (1 - x_i / k_i +  [A x]_i / k_i) + sigma * x_i
+           * dW
 
         where r_i and k_i are the growth rates and carrying capacities of
         species i, A is the matrix of interspecies interactions and sigma
         is the magnitude of the effect of the Weiner process.
 
+        See:
+        * Vadillo 2019, "Comparing stochastic Lotka–Volterra predator-prey
+        models"
+        * https://en.wikipedia.org/wiki/Competitive_Lotka%E2%80%93Volterra_equations
+
 
         Args:
             growth_rates (ndarray): A length n vector of growth rates (r_i's).
@@ -123,6 +134,6 @@ def drift(self, x, t):
         return self.dXdt(x, t)
 
     def noise(self, x, t):
-        return self.sigma * np.eye(self.dim)
+        return self.sigma * np.diag(x)
 
 

interfere/dynamics/neuroscience.py

@@ -237,5 +237,94 @@ def noise(self, x: np.ndarray, t: float):
 
 
 class LEGIONPyclustering(DiscreteTimeDynamics):
-    pass
 
+
+    def __init__(
+        self,
+        num_neurons: int,
+        sigma: float=0.0,
+        parameters: legion_parameters = DEFAULT_LEGION_PARAMETERS,
+        type_conn: str = "all_to_all",
+        measurement_noise_std: Optional[np.ndarray] = None
+    ):
+        """LEGION (local excitatory global inhibitory oscillatory network).
+        
+        Args:
+            num_neurons (int): Number of neurons in the model. Must be an even  
+                number.
+            sigma (float): Scale of the independent stochastic noise added to
+                the system.
+            parameters (hhn_parameters): A pyclustering.nnet.hhn.hhn_paramerers 
+                object.
+            type_conn (str): Type of connection between oscillators. One
+                of ["all_to_all", "grid_four", "grid_eight", "list_bdir",
+                "dynamic"]. See pyclustering.nnet.__init__::conn_type for
+                details.
+            measurement_noise_std (ndarray): None, or a vector with shape (n,)
+                where each entry corresponds to the standard deviation of the
+                measurement noise for that particular dimension of the dynamic
+                model. For example, if the dynamic model had two variables x1
+                and x2 and measurement_noise_std = [1, 10], then independent
+                gaussian noise with standard deviation 1 and 10 will be added to
+                x1 and x2 respectively at each point in time. 
+        """
+        if num_neurons % 2 == 1:
+            raise ValueError("LEGION model requires an even number of neurons.")
+
+        self.num_excite = num_neurons // 2
+        self.parameters = parameters
+        self.sigma = sigma
+        self.type_conn = type_conn
+        self.Sigma = sigma * np.diag(np.ones(num_neurons))  # Noise covariance.
+        super().__init__(num_neurons, measurement_noise_std)
+
+
+    @copy_doc(DiscreteTimeDynamics.simulate)
+    def simulate(
+        self,
+        initial_condition: np.ndarray,
+        time_points: np.ndarray,
+        intervention: Optional[Callable[[np.ndarray, float], np.ndarray]]= None,
+        rng: np.random.mtrand.RandomState = DEFAULT_RANGE,
+    ):
+        self.legion_model = legion_network(
+            self.dim // 2,
+            self.parameters,
+            CONN_TYPE_MAP[self.type_conn],
+            ccore=False
+        )
+        # Assumes equally spaced time points.
+        self.dt = (time_points[-1] - time_points[0]) / len(time_points)
+
+
+        X_do = super().simulate(
+            initial_condition, time_points, intervention, rng)
+        return X_do
+
+
+    def step(
+        self,
+        x: np.ndarray,
+        t: float = None,
+        rng: np.random.RandomState = None,
+    ):
+
+        # Unpack the state of the excitatory and inhibitory neurons
+        x_excite = x[:self.num_excite]
+        x_inhib = x[self.num_excite:]
+
+        # Overwrite the states in the legion model 
+        self.legion_model._excitatory = list(x_excite)
+        self.legion_model._global_inhibitor = list(x_inhib)
+
+        # Calulate next states and extract them.
+        self.legion_model._calculate_states(0, t, self.dt, self.dt/10)
+        x_next = np.hstack([
+            self.legion_model._excitatory, self.legion_model._global_inhibitor
+        ])
+
+        # Add stochastic system noise:
+        if self.sigma != 0:
+            x_next +=  self.Sigma @ rng.normal(0.0, np.sqrt(self.dt))
+
+        return x_next

interfere/dynamics/ornstein_uhlenbeck.py

@@ -38,8 +38,8 @@ def __init__(
             Sigma.shape[1] != mu.shape[0]
         ]):
             raise ValueError(
-                "Parameters for Lotka Voltera must have matching dimensions. "
-                "Argument shapes: "
+                "Parameters for OrnsteinUhlenback must have matching "
+                "dimensions. Argument shapes: "
                 f"\n\ttheta.shape = {theta.shape}"
                 f"\n\tmu.shape = {mu.shape}"
                 f"\n\tSigma.shape = {Sigma.shape}"

interfere/dynamics/simple_linear_sdes.py

@@ -203,7 +203,7 @@ def imag_roots_4d_linear_sde(
     """
     Sigma = sigma * np.eye(MATRIX_ALL_IMAG_EIGS.shape[0])
     A = MATRIX_ALL_IMAG_EIGS
-    return LinearSDE(A, Sigma)
+    return LinearSDE(A, Sigma, measurement_noise_std)
 
 
 def attracting_fixed_point_4d_linear_sde(

interfere/dynamics/statespace_models.py

@@ -1,6 +1,7 @@
 """State space dynamic models such as vector autoregression or VARIMAX.
 """
 from typing import Callable, List, Optional
+from warnings import warn
 
 import numpy as np
 
@@ -94,8 +95,18 @@ def simulate(
                 `time_points`. The first p rows contain the initial condition/
                 history of the system and count towards n.
         """
-        p, m = initial_condition.shape
+        p = len(self.phi_matrices)
         q = len(self.theta_matrices)
+        m = self.phi_matrices[0].shape
+
+        if len(initial_condition.shape) == 1:
+            warn("Historic timesteps not found in initial condition. Replacing with zeros")
+            initial_condition = np.vstack([
+                np.zeros(m) for i in range(max(p, q))
+            ] + [initial_condition]
+            )
+
+        p, m = initial_condition.shape
         n = len(time_points)
         X = np.zeros((n, m))
 
@@ -105,21 +116,6 @@ def simulate(
         # Initialize noise
         noise_vecs = rng.multivariate_normal(
             np.zeros(m), self.sigma, n - p + q)
-
-        # Simulate 
-        for i in range(n - p - 1):
-            x_AR = np.sum([
-                phi @ x_i
-                for x_i, phi in zip(X[i:(p+i)], self.phi_matrices[::-1])
-            ])
-            x_MA = np.sum([
-                theta @ a_i
-                for a_i, theta in zip(
-                    noise_vecs[i:(q+i)], self.theta_matrices[::-1])
-            ])
-            x_next = noise_vecs[q+i+1] + x_AR + x_MA
-
-            X[p+i+1, :] = x_next
 
         # Simulate
         for i in range(n - p):