Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add 1D option to ignore 2D functionality #51

Merged
merged 7 commits into from
Oct 5, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 25 additions & 12 deletions blobmodel/blobs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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,
)
)

Expand All @@ -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
Expand All @@ -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)
)
Expand Down Expand Up @@ -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"
Expand All @@ -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
Expand Down
16 changes: 16 additions & 0 deletions blobmodel/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from .stochasticality import BlobFactory, DefaultBlobFactory
from .geometry import Geometry
from nptyping import NDArray
import warnings


class Model:
Expand All @@ -26,6 +27,7 @@ def __init__(
blob_factory: BlobFactory = DefaultBlobFactory(),
labels: str = "off",
label_border: float = 0.75,
one_dimensional: bool = False,
) -> None:
"""
Attributes
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down
14 changes: 12 additions & 2 deletions blobmodel/stochasticality.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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"
3 changes: 3 additions & 0 deletions examples/custom_blobfactory.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
3 changes: 3 additions & 0 deletions examples/single_blob.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@ def sample_blobs(
for i in range(num_blobs)
]

def is_one_dimensional(self) -> bool:
return False


bf = CustomBlobFactory()

Expand Down
10 changes: 5 additions & 5 deletions tests/test_blob.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
pos_x=0,
pos_y=0,
t_init=0,
t_drain=10 ** 10,
t_drain=10**10,
)


Expand All @@ -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"
Expand All @@ -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))

Expand All @@ -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)
Expand All @@ -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"
Expand Down
3 changes: 3 additions & 0 deletions tests/test_labels.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
44 changes: 44 additions & 0 deletions tests/test_one_dim.py
Original file line number Diff line number Diff line change
@@ -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()
6 changes: 6 additions & 0 deletions tests/test_vx_or_vy=0.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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()

Expand Down