Merge pull request #51 from uit-cosmo/one_dim

Add 1D option to ignore 2D functionality

blobmodel/blobs.py

@@ -45,9 +45,10 @@ def discretize_blob(
         t: NDArray,
         Ly: float,
         periodic_y: bool = False,
+        one_dimensional: bool = False,
     ) -> NDArray:
         """
-        Discretize blob on grid
+        Discretize blob on grid. If one_dimensional the perpendicular pulse shape is ignored.
         The following blob shapes are implemented:
                 gauss: 2D gaussian function
                 exp: one sided exponential in x and gaussian in y
@@ -56,14 +57,19 @@ def discretize_blob(
                 -------
                 discretized blob on 3d array with dimensions x,y and t : np.array
         """
+        # If one_dimensional, then Ly should be 0.
+        assert (one_dimensional and Ly == 0) or not one_dimensional
+
         if (self.width_perp > 0.1 * Ly or self.width_prop > 0.1 * Ly) and periodic_y:
             warnings.warn("blob width big compared to Ly")
 
         x_perp, y_perp = self._rotate(
             origin=(self.pos_x, self.pos_y), x=x, y=y, angle=-self.theta
         )
-        if not periodic_y:
-            return self._single_blob(x_perp, y_perp, t, Ly, periodic_y)
+        if not periodic_y or one_dimensional:
+            return self._single_blob(
+                x_perp, y_perp, t, Ly, periodic_y, one_dimensional=one_dimensional
+            )
         if np.sin(self.theta) == 0:
             _x_border = Ly - self.pos_y
             _adjusted_Ly = Ly
@@ -92,6 +98,7 @@ def discretize_blob(
                 _number_of_y_propagations,
                 x_offset=Ly * np.sin(self.theta),
                 y_offset=Ly * np.cos(self.theta),
+                one_dimensional=one_dimensional,
             )
             + self._single_blob(
                 x_perp,
@@ -102,6 +109,7 @@ def discretize_blob(
                 _number_of_y_propagations,
                 x_offset=-Ly * np.sin(self.theta),
                 y_offset=-Ly * np.cos(self.theta),
+                one_dimensional=one_dimensional,
             )
         )
 
@@ -115,6 +123,7 @@ def _single_blob(
         number_of_y_propagations: NDArray = 0,
         x_offset: NDArray = 0,
         y_offset: NDArray = 0,
+        one_dimensional: bool = False,
     ) -> NDArray:
         return (
             self.amplitude
@@ -126,11 +135,15 @@ def _single_blob(
                 periodic_y,
                 number_of_y_propagations=number_of_y_propagations,
             )
-            * self._perpendicular_direction_shape(
-                y_perp + y_offset,
-                Ly,
-                periodic_y,
-                number_of_y_propagations=number_of_y_propagations,
+            * (
+                1
+                if one_dimensional
+                else self._perpendicular_direction_shape(
+                    y_perp + y_offset,
+                    Ly,
+                    periodic_y,
+                    number_of_y_propagations=number_of_y_propagations,
+                )
             )
             * self._blob_arrival(t)
         )
@@ -161,9 +174,9 @@ def _propagation_direction_shape(
             x_diffs = x - self._prop_dir_blob_position(t)
 
         if self.blob_shape == "gauss":
-            return 1 / np.sqrt(np.pi) * np.exp(-(x_diffs ** 2 / self.width_prop ** 2))
+            return 1 / np.sqrt(np.pi) * np.exp(-(x_diffs**2 / self.width_prop**2))
         elif self.blob_shape == "exp":
-            return np.exp(x_diffs) * np.heaviside(-1.0 * (x_diffs), 1)
+            return np.exp(x_diffs / self.width_prop) * np.heaviside(-1.0 * (x_diffs), 1)
         else:
             raise NotImplementedError(
                 f"{self.__class__.__name__}.blob_shape not implemented"
@@ -184,10 +197,10 @@ def _perpendicular_direction_shape(
             )
         else:
             y_diffs = y - self._perp_dir_blob_position()
-        return 1 / np.sqrt(np.pi) * np.exp(-(y_diffs ** 2) / self.width_perp ** 2)
+        return 1 / np.sqrt(np.pi) * np.exp(-(y_diffs**2) / self.width_perp**2)
 
     def _prop_dir_blob_position(self, t: NDArray) -> NDArray:
-        return self.pos_x + (self.v_x ** 2 + self.v_y ** 2) ** 0.5 * (t - self.t_init)
+        return self.pos_x + (self.v_x**2 + self.v_y**2) ** 0.5 * (t - self.t_init)
 
     def _perp_dir_blob_position(self) -> NDArray:
         return self.pos_y

blobmodel/model.py

@@ -6,6 +6,7 @@
 from .stochasticality import BlobFactory, DefaultBlobFactory
 from .geometry import Geometry
 from nptyping import NDArray
+import warnings
 
 
 class Model:
@@ -26,6 +27,7 @@ def __init__(
         blob_factory: BlobFactory = DefaultBlobFactory(),
         labels: str = "off",
         label_border: float = 0.75,
+        one_dimensional: bool = False,
     ) -> None:
         """
         Attributes
@@ -56,7 +58,20 @@ def __init__(
         label_border: float, optional
             defines region of blob as region where density >= label_border * amplitude of Blob
             only used if labels = "same" or "individual"
+        one_dimensional: bool, optional
+            If True, the perpendicular shape of the blobs will be discarded. Parameters
+            for the y-component (Ny and Ly) will be overwritten to Ny=1, Ly=0.
         """
+        self._one_dimensional = one_dimensional
+        if self._one_dimensional and (Ny != 1 or Ly != 0):
+            warnings.warn("Overwritting Ny=1 and Ly=0 to allow one dimensional model")
+            Ny = 1
+            Ly = 0
+            if not blob_factory.is_one_dimensional():
+                warnings.warn(
+                    "Using a one dimensional model with a with a blob factory that is not one-dimensional. Are you sure you know what you are doing?"
+                )
+
         self._geometry: Geometry = Geometry(
             Nx=Nx,
             Ny=Ny,
@@ -184,6 +199,7 @@ def _sum_up_blobs(
             t=self._geometry.t_matrix[:, :, _start:_stop],
             periodic_y=self._geometry.periodic_y,
             Ly=self._geometry.Ly,
+            one_dimensional=self._one_dimensional,
         )
 
         self._density[:, :, _start:_stop] += _single_blob

blobmodel/stochasticality.py

@@ -16,6 +16,12 @@ def sample_blobs(
         """creates list of Blobs used in Model."""
         raise NotImplementedError
 
+    @abstractmethod
+    def is_one_dimensional(self) -> bool:
+        """returns True if the BlobFactory is compatible with a one_dimensional
+        model."""
+        raise NotImplementedError
+
 
 class DefaultBlobFactory(BlobFactory):
     """Default implementation of BlobFactory.
@@ -42,7 +48,7 @@ def __init__(
         normal: normal distribution with zero mean and `free_parameter` as scale parameter
         uniform: uniorm distribution with mean 1 and `free_parameter` as width
         ray: rayleight distribution with mean 1
-        deg: array on ones
+        deg: degenerate distribution at `free_parameter`
         zeros: array of zeros
         """
         self.A_dist = A_dist
@@ -76,7 +82,7 @@ def _draw_random_variables(
         elif dist_type == "ray":
             return np.random.rayleigh(scale=np.sqrt(2.0 / np.pi), size=num_blobs)
         elif dist_type == "deg":
-            return np.ones(num_blobs)
+            return free_parameter * np.ones(num_blobs)
         elif dist_type == "zeros":
             return np.zeros(num_blobs)
         else:
@@ -119,3 +125,7 @@ def sample_blobs(
 
         # sort blobs by amplitude
         return np.array(_Blobs)[np.argsort(_amp)]
+
+    def is_one_dimensional(self) -> bool:
+        # Perpendicular width parameters are irrelevant since perp shape should be ignored by the Bolb class.
+        return self.vy_dist == "zeros"

examples/custom_blobfactory.py

@@ -41,6 +41,9 @@ def sample_blobs(
             for i in range(num_blobs)
         ]
 
+    def is_one_dimensional(self) -> bool:
+        return False
+
 
 bf = CustomBlobFactory()
 tmp = Model(

examples/single_blob.py

@@ -37,6 +37,9 @@ def sample_blobs(
             for i in range(num_blobs)
         ]
 
+    def is_one_dimensional(self) -> bool:
+        return False
+
 
 bf = CustomBlobFactory()
 

tests/test_blob.py

@@ -12,7 +12,7 @@
     pos_x=0,
     pos_y=0,
     t_init=0,
-    t_drain=10 ** 10,
+    t_drain=10**10,
 )
 
 
@@ -23,7 +23,7 @@ def test_initial_blob():
     mesh_x, mesh_y = np.meshgrid(x, y)
     blob_values = blob.discretize_blob(x=mesh_x, y=mesh_y, t=0, periodic_y=False, Ly=10)
 
-    expected_values = 1 / np.pi * np.exp(-(mesh_x ** 2)) * np.exp(-(mesh_y ** 2))
+    expected_values = 1 / np.pi * np.exp(-(mesh_x**2)) * np.exp(-(mesh_y**2))
     error = np.max(abs(expected_values - blob_values))
 
     assert error < 1e-10, "Numerical error too big"
@@ -39,7 +39,7 @@ def test_periodicity():
     )
 
     expected_values = 1 / np.pi * np.exp(-((mesh_x - 2) ** 2)) * np.exp(
-        -(mesh_y ** 2)
+        -(mesh_y**2)
     ) + 1 / np.pi * np.exp(-((mesh_x - 2) ** 2)) * np.exp(-((mesh_y - 10) ** 2))
     error = np.max(abs(expected_values - blob_values))
 
@@ -58,7 +58,7 @@ def test_single_point():
         pos_x=0,
         pos_y=6,
         t_init=0,
-        t_drain=10 ** 100,
+        t_drain=10**100,
     )
 
     x = np.arange(0, 10, 0.1)
@@ -69,7 +69,7 @@ def test_single_point():
         x=mesh_x, y=mesh_y, t=0, periodic_y=True, Ly=10
     )
 
-    expected_values = 1 / np.pi * np.exp(-(mesh_x ** 2)) * np.exp(-(4 ** 2))
+    expected_values = 1 / np.pi * np.exp(-(mesh_x**2)) * np.exp(-(4**2))
     error = np.max(abs(expected_values - blob_values))
 
     assert error < 1e-10, "Numerical error too big"

tests/test_labels.py

@@ -39,6 +39,9 @@ def sample_blobs(
             for i in range(num_blobs)
         ]
 
+    def is_one_dimensional(self) -> bool:
+        return False
+
 
 def test_bloblabels_speedup():
     warnings.filterwarnings("ignore")

tests/test_one_dim.py

@@ -0,0 +1,44 @@
+from blobmodel import Model, DefaultBlobFactory
+import xarray as xr
+import numpy as np
+
+
+# use DefaultBlobFactory to define distribution functions fo random variables
+bf = DefaultBlobFactory(A_dist="deg", W_dist="deg", vx_dist="deg", vy_dist="zeros")
+
+one_dim_model = Model(
+    Nx=100,
+    Ny=1,
+    Lx=10,
+    Ly=0,
+    dt=1,
+    T=1000,
+    blob_shape="exp",
+    t_drain=2,
+    periodic_y=False,
+    num_blobs=10000,
+    blob_factory=bf,
+    one_dimensional=True,
+)
+
+
+def test_one_dim():
+    ds = one_dim_model.make_realization(speed_up=True, error=1e-2)
+    model_profile = ds.n.isel(y=0).mean(dim=("t"))
+
+    x = np.linspace(0, 10, 100)
+    t_p = 1
+    t_w = 1 / 10
+    amp = 1
+    v_p = 1.0
+    t_loss = 2.0
+    t_d = t_loss * t_p / (t_loss + t_p)
+
+    analytical_profile = t_d / t_w * amp * np.exp(-x / (v_p * t_loss))
+
+    error = np.mean(abs(model_profile.values - analytical_profile))
+
+    assert error < 0.1, "Numerical error too big"
+
+
+test_one_dim()

tests/test_vx_or_vy=0.py

@@ -38,6 +38,9 @@ def sample_blobs(
             for i in range(num_blobs)
         ]
 
+    def is_one_dimensional(self) -> bool:
+        return False
+
 
 class CustomBlobFactoryVx0(BlobFactory):
     def __init__(self) -> None:
@@ -74,6 +77,9 @@ def sample_blobs(
             for i in range(num_blobs)
         ]
 
+    def is_one_dimensional(self) -> bool:
+        return False
+
 
 bf_vy_0 = CustomBlobFactoryVy0()