Merge pull request #46 from grburgess/dev

Dev

.flake8

@@ -1,5 +1,5 @@
 [flake8]
-ignore = E127,E201,E202,E203,E231,E252,E266,E402,E999,F841,W503,W605,F403,F401
+ignore = E127,E201,E202,E203,E231,E252,E266,E402,E999,F841,W503,W605,F403,F401, E731
 select = C,E,F,W,B,B950
 extend-ignore = E203, E501
 #extend-ignore = F401 # ignore unused

.github/workflows/build_test.yml

@@ -24,7 +24,8 @@ jobs:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
-        pip install --upgrade flake8 numba numpy ipython scipy h5py matplotlib hypothesis pystan betagen coverage pytest pytest-cov cython codecov
+        python -m pip install --upgrade pip wheel
+        pip install --upgrade flake8 numba numpy ipython ipython_genutils scipy h5py matplotlib hypothesis pystan betagen coverage pytest pytest-cov cython codecov
         python setup.py install
     - name: Lint with flake8
       run: |

popsynth/aux_samplers/viewing_angle_sampler.py

@@ -36,3 +36,9 @@ def true_sampler(self, size: int) -> None:
         )
 
         self._true_values = np.arccos(1.0 - theta_inverse)
+
+    def _compute_probability(self):
+
+        return np.sin(self._true_values) / (
+            1 - np.cos(np.deg2rad(self.max_angle))
+        )

popsynth/auxiliary_sampler.py

@@ -27,6 +27,7 @@ def __init__(
         true_values: ArrayLike,
         obs_values: ArrayLike,
         selection: ArrayLike,
+        probability: ArrayLike,
     ) -> None:
         """
         A container for secondary properties that adds dict
@@ -40,14 +41,17 @@ def __init__(
         :type obs_values: ArrayLike
         :param selection:
         :type selection: ArrayLike
+        :param probability:
+        :type probability: ArrayLike
+        :returns:
         :returns:
 
         """
 
         self._true_values: ArrayLike = true_values
         self._obs_values: ArrayLike = obs_values
         self._selection: ArrayLike = selection
-
+        self._probability: ArrayLike = probability
         self._name: str = name
 
     @property
@@ -84,6 +88,16 @@ def selection(self) -> ArrayLike:
         """
         return self._selection
 
+    @property
+    def probability(self) -> ArrayLike:
+        """
+        The the probability of the draws
+
+        :returns:
+
+        """
+        return self._probability
+
     def __getitem__(self, key):
 
         if key == "selection":
@@ -95,6 +109,9 @@ def __getitem__(self, key):
         elif key == "obs_values":
             return self._obs_values
 
+        elif key == "probability":
+            return self._probability
+
         else:
 
             log.error("trying to access something that does not exist")
@@ -177,22 +194,23 @@ def __init__(
         :type uses_sky_position: bool
         """
 
-        self._parameter_storage = {}  # type: Dict[str, float]
-        self._name = name  # type: str
-        self._obs_name = "%s_obs" % name  # type: str
+        self._parameter_storage: Dict[str, float] = {}
+        self._name: str = name
+        self._obs_name: str = "%s_obs" % name
 
-        self._obs_values = None  # type: ArrayLike
-        self._true_values = None  # type: ArrayLike
-        self._is_observed = observed  # type: bool
-        self._secondary_samplers = {}  # type: SamplerDict
-        self._is_secondary = False  # type: bool
+        self._obs_values: ArrayLike = None
+        self._true_values: ArrayLike = None
+        self._is_observed: bool = observed
+        self._secondary_samplers: SamplerDict = {}
+        self._is_secondary: bool = False
         self._parent_names = []
-        self._has_secondary = False  # type: bool
-        self._is_sampled = False  # type: bool
-        self._selector = UnitySelection()  # type: SelectionProbability
-        self._uses_distance = uses_distance  # type: bool
-        self._uses_luminosity = uses_luminosity  # type: bool
-        self._uses_sky_position = uses_sky_position  # type: bool
+        self._has_secondary: bool = False
+        self._is_sampled: bool = False
+        self._selector: SelectionProbability = UnitySelection()
+        self._uses_distance: bool = uses_distance
+        self._uses_luminosity: bool = uses_luminosity
+        self._uses_sky_position: bool = uses_sky_position
+        self._probability: Optional[np.ndarray] = None
 
     def display(self):
 
@@ -378,6 +396,22 @@ def draw(self, size: int = 1):
 
             self._apply_selection()
 
+            self._probability = self._compute_probability()
+
+    def _compute_probability(self) -> Optional[np.ndarray]:
+        return None
+
+    @property
+    def probability(self) -> np.ndarray:
+
+        if self._probability is None:
+
+            return np.ones_like(self._true_values)
+
+        else:
+
+            return self._probability
+
     def reset(self):
         """
         Reset all the selections.
@@ -457,7 +491,11 @@ def get_secondary_properties(
 
         recursive_secondaries.add_secondary(
             SecondaryContainer(
-                self._name, self._true_values, self._obs_values, self._selector
+                self._name,
+                self._true_values,
+                self._obs_values,
+                self._selector,
+                self.probability,
             )
         )
 
@@ -670,6 +708,9 @@ def true_sampler(self, size: int = 1):
 
         pass
 
+    def _probability(self, true_values: np.ndarray) -> np.ndarray:
+        pass
+
     def observation_sampler(self, size: int = 1) -> np.ndarray:
 
         return self._true_values

popsynth/distribution.py

@@ -1,9 +1,10 @@
 import abc
 from dataclasses import dataclass
-from typing import Dict, Union
+from typing import Dict, Optional, Union
 
 import numpy as np
 import pandas as pd
+import scipy.integrate as integrate
 from class_registry import AutoRegister
 from IPython.display import Markdown, Math, display
 from numpy.typing import ArrayLike
@@ -42,11 +43,23 @@ def __init__(self, name: str, seed: int, form: str) -> None:
         """
 
         self._parameter_storage: Dict[str, float] = {}
+        self._normalization_parameter: Optional[str] = None
 
         self._seed = seed  # type: int
         self._name = name  # type: str
         self._form = form  # type: str
 
+        self._probability: Optional[np.ndarray] = None
+
+    @property
+    def normalization_parameter(self) -> Optional[float]:
+
+        return self._parameter_storage[self._normalization_parameter]
+
+    @normalization_parameter.setter
+    def normalization_parameter(self, value):
+        self._parameter_storage[self._normalization_parameter] = value
+
     @property
     def name(self) -> str:
         """
@@ -91,6 +104,11 @@ def truth(self) -> Dict[str, float]:
 
         return out
 
+    @property
+    def probability(self) -> np.ndarray:
+
+        return self._probability
+
     def display(self):
         """
         use ipython pretty display to
@@ -274,7 +292,7 @@ def spatial_values(self):
 
     def draw_sky_positions(self, size: int) -> None:
         """
-        Draw teh sky positions of the objects
+        Draw the sky positions of the objects
 
         :param size:
         :type size: int
@@ -283,7 +301,7 @@ def draw_sky_positions(self, size: int) -> None:
         """
         self._theta, self._phi = sample_theta_phi(size)
 
-    def draw_distance(self, size: int) -> None:
+    def draw_distance(self, size: int, normalize: bool = False) -> None:
         """
         Draw the distances from the specified dN/dr model.
 
@@ -298,11 +316,32 @@ def draw_distance(self, size: int) -> None:
             / self.time_adjustment(r)
         )
 
+        if normalize:
+
+            # set the prefactor to unity
+
+            if self._normalization_parameter is not None:
+
+                old_value = self.normalization_parameter
+                self.normalization_parameter = 1
+
+            integral = integrate.quad(dNdr, 0.0, self.r_max)[0]
+
+            dNdr_norm = lambda r: dNdr(r) / integral
+
+            if self._normalization_parameter is not None:
+
+                self.normalization_parameter = old_value
+
+        else:
+
+            dNdr_norm = dNdr
+
         # find the maximum point
         tmp = np.linspace(
             0.0, self.r_max, 500, dtype=np.float64
         )  # type: ArrayLike
-        ymax = np.max(dNdr(tmp))  # type: float
+        ymax = np.max(dNdr_norm(tmp))  # type: float
 
         # rejection sampling the distribution
         r_out = []
@@ -324,7 +363,7 @@ def draw_distance(self, size: int) -> None:
 
                     # compare them
 
-                    if y < dNdr(r):
+                    if y < dNdr_norm(r):
                         r_out.append(r)
                         flag = False
         else:
@@ -335,6 +374,20 @@ def draw_distance(self, size: int) -> None:
 
         self._distances = np.array(r_out)  # type: ArrayLike
 
+        # now compute the differential probability
+        # of the distance draws
+
+        if self._normalization_parameter is not None:
+
+            old_value = self.normalization_parameter
+            self.normalization_parameter = 1
+
+        self._probability = dNdr_norm(self._distances)
+
+        if self._normalization_parameter is not None:
+
+            self.normalization_parameter = old_value
+
 
 class LuminosityDistribution(Distribution):
     _distribution_name = "LuminosityDistribtuion"
@@ -368,7 +421,7 @@ def phi(self, luminosity):
         raise RuntimeError("Must be implemented in derived class")
 
     @abc.abstractmethod
-    def draw_luminosity(self, size):
+    def draw_luminosity(self, size: int) -> np.ndarray:
         """
         function to draw the luminosity via an alternative method
         must be implemented in child class
@@ -379,3 +432,7 @@ def draw_luminosity(self, size):
 
         """
         pass
+
+    def compute_probability(self, L: np.ndarray) -> None:
+
+        self._probability = self.phi(L)

popsynth/distributions/__init__.py

@@ -16,6 +16,10 @@
     ZPowerSphericalDistribution,
 )
 from .spiral_galaxy_distribution import SpiralGalaxyDistribution
+from .uniform_distribution import (
+    LogUniLuminiosityDistribution,
+    UniformCosmoDistribution,
+)
 
 __all__ = [
     "SphericalDistribution",
@@ -33,4 +37,6 @@
     "DeltaDistribution",
     "FlatlandDistribution",
     "SpiralGalaxyDistribution",
+    "LogUniLuminiosityDistribution",
+    "UniformCosmoDistribution",
 ]

popsynth/distributions/bpl_distribution.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy.integrate as integrate
 import scipy.stats as stats
 
 from popsynth.distribution import DistributionParameter, LuminosityDistribution
@@ -45,7 +46,15 @@ def __init__(self, seed: int = 1234, name: str = "bpl"):
 
     def phi(self, L):
 
-        return bpl(L, self.Lmin, self.Lbreak, self.Lmax, self.alpha, self.beta)
+        f = lambda ll: bpl(
+            ll, self.Lmin, self.Lbreak, self.Lmax, self.alpha, self.beta
+        )
+
+        integrand = integrate.quad(f, self.Lmin, self.Lmax)[0]
+
+        # normalize
+
+        return f(L) / integrand
 
     def draw_luminosity(self, size=1):
 
@@ -98,20 +107,25 @@ def bpl(x, x0, x1, x2, a1, a2):
     """
 
     # creatre a holder for the values
-    out = np.empty_like(x)
+    x = np.atleast_1d(x)
+
+    out = np.zeros_like(x)
 
     # get the total integral to compute the normalization
     _, _, C = integrate_pl(x0, x1, x2, a1, a2)
     norm = 1.0 / C
 
     # create an index to select each piece of the function
-    idx = x < x1
+    idx = (x > x0) & (x < x1)
 
     # compute the lower power law
     out[idx] = np.power(x[idx], a1)
 
     # compute the upper power law
-    out[~idx] = np.power(x[~idx], a2) * np.power(x1, a1 - a2)
+
+    idx = (x >= x1) & (x < x2)
+
+    out[idx] = np.power(x[idx], a2) * np.power(x1, a1 - a2)
 
     return out * norm
 

popsynth/distributions/cosmological_distribution.py

@@ -91,7 +91,7 @@ def time_adjustment(self, z):
 class SFRDistribution(CosmologicalDistribution):
     _distribution_name = "SFRDistribution"
 
-    r0 = DistributionParameter(vmin=0)
+    r0 = DistributionParameter(vmin=0, is_normalization=True)
     a = DistributionParameter(vmin=0)
     rise = DistributionParameter()
     decay = DistributionParameter()
@@ -157,7 +157,7 @@ def _sfr_dndv(z, r0, a, rise, decay, peak):
 class ZPowerCosmoDistribution(CosmologicalDistribution):
     _distribution_name = "ZPowerCosmoDistribution"
 
-    Lambda = DistributionParameter(default=1, vmin=0)
+    Lambda = DistributionParameter(default=1, vmin=0, is_normalization=True)
     delta = DistributionParameter()
 
     def __init__(