From abf3bb0c729b2e10bf779c0e99839d387bdfaffc Mon Sep 17 00:00:00 2001 From: grburgess Date: Thu, 21 Apr 2022 09:36:09 -0400 Subject: [PATCH 1/9] starting to move to adding probs per draw --- popsynth/auxiliary_sampler.py | 5 +++++ popsynth/distribution.py | 16 +++++++++++++++- popsynth/population_synth.py | 35 +++++++++++++++++++++-------------- 3 files changed, 41 insertions(+), 15 deletions(-) diff --git a/popsynth/auxiliary_sampler.py b/popsynth/auxiliary_sampler.py index 19cced5..c19176c 100644 --- a/popsynth/auxiliary_sampler.py +++ b/popsynth/auxiliary_sampler.py @@ -670,6 +670,11 @@ 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 diff --git a/popsynth/distribution.py b/popsynth/distribution.py index 6b98274..fd6b4e7 100644 --- a/popsynth/distribution.py +++ b/popsynth/distribution.py @@ -1,9 +1,10 @@ import abc from dataclasses import dataclass -from typing import Dict, Union +from typing import Dict, Union, Optional 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 @@ -47,6 +48,8 @@ def __init__(self, name: str, seed: int, form: str) -> None: self._name = name # type: str self._form = form # type: str + self._probability: Optional[np.ndarray] = None + @property def name(self) -> str: """ @@ -91,6 +94,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 @@ -335,6 +343,12 @@ def draw_distance(self, size: int) -> None: self._distances = np.array(r_out) # type: ArrayLike + # now compute the differential probability + # of the distance draws + + integral = integrate.quad(dNdr, 0.0, self.r_max)[0] + + self._probability = dNdr(r_out) / integral class LuminosityDistribution(Distribution): _distribution_name = "LuminosityDistribtuion" diff --git a/popsynth/population_synth.py b/popsynth/population_synth.py index 9a6209f..4401628 100644 --- a/popsynth/population_synth.py +++ b/popsynth/population_synth.py @@ -907,6 +907,7 @@ def draw_survey( self, flux_sigma: Optional[float] = None, log10_flux_draw: bool = True, + n_samples: Optional[int] = None, ) -> Population: """ Draw the total survey and return a :class:`Population` object. @@ -936,23 +937,29 @@ def draw_survey( np.random.seed(self._seed) - # create a callback of the integrand - dNdr = ( - lambda r: self._spatial_distribution.dNdV(r) - * self._spatial_distribution.differential_volume(r) - / self._spatial_distribution.time_adjustment(r) - ) + if n_samples is None: + + # create a callback of the integrand + dNdr = ( + lambda r: self._spatial_distribution.dNdV(r) + * self._spatial_distribution.differential_volume(r) + / self._spatial_distribution.time_adjustment(r) + ) - # integrate the population to determine the true number of - # objects - N = integrate.quad(dNdr, 0.0, self._spatial_distribution.r_max)[ - 0 - ] # type: float + # integrate the population to determine the true number of + # objects + N = integrate.quad(dNdr, 0.0, self._spatial_distribution.r_max)[ + 0 + ] # type: float - log.info("The volume integral is %f" % N) + log.info("The volume integral is %f" % N) + + # this should be poisson distributed + n: int = np.random.poisson(N) + + else: - # this should be poisson distributed - n = np.random.poisson(N) # type: np.int64 + n: int = n_samples self._spatial_distribution.draw_distance(size=n) From 02cdf788977661160043dd23a4e57bcfe46d2629 Mon Sep 17 00:00:00 2001 From: grburgess Date: Tue, 26 Apr 2022 15:39:57 +0200 Subject: [PATCH 2/9] added uniform distributions --- popsynth/distributions/__init__.py | 6 ++ .../distributions/uniform_distribution.py | 95 +++++++++++++++++++ 2 files changed, 101 insertions(+) create mode 100644 popsynth/distributions/uniform_distribution.py diff --git a/popsynth/distributions/__init__.py b/popsynth/distributions/__init__.py index 206fdbd..a1fb8e9 100644 --- a/popsynth/distributions/__init__.py +++ b/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", ] diff --git a/popsynth/distributions/uniform_distribution.py b/popsynth/distributions/uniform_distribution.py new file mode 100644 index 0000000..c986264 --- /dev/null +++ b/popsynth/distributions/uniform_distribution.py @@ -0,0 +1,95 @@ +import numpy as np +import scipy.integrate as integrate +import scipy.stats as stats + +from popsynth.distribution import DistributionParameter, LuminosityDistribution +from popsynth.distributions.cosmological_distribution import ( + CosmologicalDistribution, +) + + +class LogUniLuminiosityDistribution(LuminosityDistribution): + _distribution_name = "LogUniLuminiosityDistribution" + + Lmin = DistributionParameter(vmin=0) + + Lmax = DistributionParameter(vmin=0) + + def __init__(self, seed: int = 1234, name: str = "logunilum"): + """A broken power law luminosity distribution. + + L ~ L^``alpha`` for L <= ``Lbreak`` + L ~ L^``beta`` for L > ``Lbreak`` + + :param seed: Random seed + :type seed: int + :param name: Name of the distribution + :type name: str + :param Lmin: Minimum value of the luminosity + :type Lmin: :class:`DistributionParameter` + :param Lmax: Maximum value of the luminosity + :type Lmax: :class:`DistributionParameter` + """ + lf_form = r"\begin{cases} C L^{\alpha} & \mbox{if } L" + lf_form += r"\leq L_\mathrm{break},\\ C L^{\beta} " + lf_form += r"L_\mathrm{break}^{\alpha - \beta}" + lf_form += r" & \mbox{if } L > L_\mathrm{break}. \end{cases}" + + super(LogUniLuminiosityDistribution, self).__init__( + seed=seed, name=name, form=lf_form + ) + + def phi(self, L): + + return 1.0 / (L * np.log(self.Lmax / self.Lmin)) + + def draw_luminosity(self, size=1): + + return stats.loguniform.rvs(self.Lmin, self.Lmax, size=size) + + +class UniformCosmoDistribution(CosmologicalDistribution): + _distribution_name = "UniformCosmoDistribution" + + r0 = DistributionParameter(vmin=0, is_normalization=True) + zmin = DistributionParameter(vmin=0) + zmax = DistributionParameter(vmin=0) + + def __init__( + self, + seed: int = 1234, + name: str = "uniform_cosmo", + is_rate: bool = True, + ): + """ + A cosmological distribution where the density + evolves uniformly. + + ``Lambda`` (1+z)^``delta`` + + :param seed: Random seed + :type seed: int + :param name: Name of the distribution + :type name: str + :param is_rate: `True` if modelling a population of transient events, + `False` if modelling a population of steady-state objects. + Affects the ``time_adjustment`` method used in cosmo calculations. + Default is `True`. + :type is_rate: bool + :param Lambda: The local density in units of Gpc^-3 + :type Lambda: :class:`DistributionParameter` + :param delta: The index of the power law + :type delta: :class:`DistributionParameter` + """ + spatial_form = r"\Lambda (z+1)^{\delta}" + + super(UniformCosmoDistribution, self).__init__( + seed=seed, + name=name, + form=spatial_form, + is_rate=is_rate, + ) + + def dNdV(self, distance): + + return self.r0 / (self.zmax - self.zmin) From 0886ed573e231e6d89896be443fc78ffab4dd581 Mon Sep 17 00:00:00 2001 From: grburgess Date: Tue, 26 Apr 2022 15:40:11 +0200 Subject: [PATCH 3/9] added the ability to normalize spatial distributions and compute probability of the draws --- popsynth/distribution.py | 52 +++++++++++--- popsynth/distributions/bpl_distribution.py | 13 +++- .../cosmological_distribution.py | 4 +- .../distributions/spherical_distribution.py | 2 +- popsynth/population.py | 72 +++++++++++++++---- popsynth/population_synth.py | 67 ++++++++++++++--- popsynth/utils/meta.py | 34 +++++++-- 7 files changed, 204 insertions(+), 40 deletions(-) diff --git a/popsynth/distribution.py b/popsynth/distribution.py index fd6b4e7..99dd135 100644 --- a/popsynth/distribution.py +++ b/popsynth/distribution.py @@ -1,6 +1,6 @@ import abc from dataclasses import dataclass -from typing import Dict, Union, Optional +from typing import Dict, Optional, Union import numpy as np import pandas as pd @@ -43,6 +43,7 @@ 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 @@ -50,6 +51,15 @@ def __init__(self, name: str, seed: int, form: str) -> None: 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: """ @@ -282,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 @@ -291,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. @@ -306,11 +316,34 @@ 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 = [] @@ -332,7 +365,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: @@ -346,9 +379,8 @@ def draw_distance(self, size: int) -> None: # now compute the differential probability # of the distance draws - integral = integrate.quad(dNdr, 0.0, self.r_max)[0] + self._probability = dNdr_norm(self._distances) - self._probability = dNdr(r_out) / integral class LuminosityDistribution(Distribution): _distribution_name = "LuminosityDistribtuion" @@ -382,7 +414,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 @@ -393,3 +425,7 @@ def draw_luminosity(self, size): """ pass + + def compute_probability(self, L: np.ndarray) -> None: + + self._probability = self.phi(L) diff --git a/popsynth/distributions/bpl_distribution.py b/popsynth/distributions/bpl_distribution.py index 4460dcc..8f32f50 100644 --- a/popsynth/distributions/bpl_distribution.py +++ b/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,6 +107,8 @@ def bpl(x, x0, x1, x2, a1, a2): """ # creatre a holder for the values + x = np.atleast_1d(x) + out = np.empty_like(x) # get the total integral to compute the normalization diff --git a/popsynth/distributions/cosmological_distribution.py b/popsynth/distributions/cosmological_distribution.py index 880dcbe..492ed3a 100644 --- a/popsynth/distributions/cosmological_distribution.py +++ b/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__( diff --git a/popsynth/distributions/spherical_distribution.py b/popsynth/distributions/spherical_distribution.py index f9fe51b..9f6eb16 100644 --- a/popsynth/distributions/spherical_distribution.py +++ b/popsynth/distributions/spherical_distribution.py @@ -43,7 +43,7 @@ def transform(self, L, r): class ConstantSphericalDistribution(SphericalDistribution): _distribution_name = "ConstantSphericalDistribution" - Lambda = DistributionParameter(default=1, vmin=0) + Lambda = DistributionParameter(default=1, vmin=0 , is_normalization=True) def __init__( self, diff --git a/popsynth/population.py b/popsynth/population.py index 6deb7be..dc4b419 100644 --- a/popsynth/population.py +++ b/popsynth/population.py @@ -57,6 +57,7 @@ def __init__( theta=None, phi=None, pop_synth: Optional[Dict[str, Any]] = None, + probability: Optional[np.ndarray] = None, ) -> None: """ A population container for all the properties of the population. @@ -167,6 +168,16 @@ def __init__( self._pop_synth: Optional[Dict[str, Any]] = pop_synth + if probability is not None: + + self._probability: Optional[SimulatedVariable] = SimulatedVariable( + probability, probability, selection + ) + + else: + + self._probability = None + if self._n_detections == 0: self._no_detection = True @@ -221,6 +232,10 @@ def graph(self): return self._graph + @property + def probability(self) -> Optional[SimulatedVariable]: + return self._probability + @property def pop_synth(self) -> Dict[str, Any]: """ @@ -570,7 +585,9 @@ def _writeto(self, f): for k, v in self._spatial_params.items(): - spatial_grp.create_dataset(k, data=np.array([v]), compression="lzf") + spatial_grp.create_dataset( + k, data=np.array([v]), compression="gzip" + ) if self._lf_params is not None: @@ -578,7 +595,9 @@ def _writeto(self, f): for k, v in self._lf_params.items(): - lum_grp.create_dataset(k, data=np.array([v]), compression="lzf") + lum_grp.create_dataset( + k, data=np.array([v]), compression="gzip" + ) f.attrs["lf_form"] = np.string_(self._lf_form) @@ -594,28 +613,33 @@ def _writeto(self, f): f.attrs["r_max"] = self._r_max f.attrs["seed"] = int(self._seed) + if self._probability is not None: + f.create_dataset( + "probabilty", data=self._probability, compression="gzip" + ) + f.create_dataset( - "luminosities", data=self._luminosities, compression="lzf" + "luminosities", data=self._luminosities, compression="gzip" ) - f.create_dataset("distances", data=self._distances, compression="lzf") + f.create_dataset("distances", data=self._distances, compression="gzip") f.create_dataset( - "known_distances", data=self._known_distances, compression="lzf" + "known_distances", data=self._known_distances, compression="gzip" ) f.create_dataset( "known_distance_idx", data=self._known_distance_idx, - compression="lzf", + compression="gzip", ) f.create_dataset( "unknown_distance_idx", data=self._unknown_distance_idx, - compression="lzf", + compression="gzip", ) - f.create_dataset("fluxes", data=self._fluxes.latent, compression="lzf") - f.create_dataset("flux_obs", data=self._fluxes, compression="lzf") - f.create_dataset("selection", data=self._selection, compression="lzf") - f.create_dataset("theta", data=self._theta, compression="lzf") - f.create_dataset("phi", data=self._phi, compression="lzf") + f.create_dataset("fluxes", data=self._fluxes.latent, compression="gzip") + f.create_dataset("flux_obs", data=self._fluxes, compression="gzip") + f.create_dataset("selection", data=self._selection, compression="gzip") + f.create_dataset("theta", data=self._theta, compression="gzip") + f.create_dataset("phi", data=self._phi, compression="gzip") aux_grp = f.create_group("auxiliary_quantities") @@ -623,17 +647,17 @@ def _writeto(self, f): q_grp = aux_grp.create_group(k) q_grp.create_dataset( - "true_values", data=v["true_values"], compression="lzf" + "true_values", data=v["true_values"], compression="gzip" ) q_grp.create_dataset( - "obs_values", data=v["obs_values"], compression="lzf" + "obs_values", data=v["obs_values"], compression="gzip" ) model_grp = f.create_group("model_spaces") for k, v in self._model_spaces.items(): - model_grp.create_dataset(k, data=v, compression="lzf") + model_grp.create_dataset(k, data=v, compression="gzip") # now store the truths recursively_save_dict_contents_to_group(f, "truth", self._truth) @@ -736,6 +760,14 @@ def _loadfrom(cls, f): known_distance_idx = None unknown_distance_idx = None + try: + + probability = f["probability"][()] + + except AttributeError: + + probability = None + fluxes = f["fluxes"][()] flux_obs = f["flux_obs"][()] selection = f["selection"][()] @@ -790,6 +822,7 @@ def _loadfrom(cls, f): theta=theta, phi=phi, pop_synth=pop_synth, + probability=probability, ) def to_sub_population(self, observed: bool = True) -> "Population": @@ -854,6 +887,14 @@ def _selector(x): itr += 1 + if self._probability is not None: + + probability = _selector(self._probability) + + else: + + probability = None + return Population( luminosities=_selector(self._luminosities), distances=_selector(self._distances), @@ -878,6 +919,7 @@ def _selector(x): graph=self._graph, theta=_selector(self._theta), phi=_selector(self._phi), + probability=probability, ) def display(self): diff --git a/popsynth/population_synth.py b/popsynth/population_synth.py index 4401628..f0a1eb0 100644 --- a/popsynth/population_synth.py +++ b/popsynth/population_synth.py @@ -906,8 +906,9 @@ def draw_log_fobs(self, f, f_sigma, size=1) -> np.ndarray: def draw_survey( self, flux_sigma: Optional[float] = None, - log10_flux_draw: bool = True, + log10_flux_rdraw: bool = True, n_samples: Optional[int] = None, + normalize: bool = False, ) -> Population: """ Draw the total survey and return a :class:`Population` object. @@ -946,13 +947,52 @@ def draw_survey( / self._spatial_distribution.time_adjustment(r) ) - # integrate the population to determine the true number of - # objects - N = integrate.quad(dNdr, 0.0, self._spatial_distribution.r_max)[ - 0 - ] # type: float + if normalize: - log.info("The volume integral is %f" % N) + log.info("The population is being normalized such that") + log.info("the integral over N * dV/dz = N") + + old_value = None + + if ( + self._spatial_distribution._normalization_parameter + is not None + ): + + old_value = ( + self._spatial_distribution.normalization_parameter + ) + log.info('setting normalization to unity for integration') + self._spatial_distribution.normalization_parameter = 1 + + integral = integrate.quad( + dNdr, 0.0, self._spatial_distribution.r_max + )[0] + + if old_value is not None: + + self._spatial_distribution.normalization_parameter = ( + old_value + ) + + log.info( + f"normalization parameter set back to {self._spatial_distribution.normalization_parameter}" + ) + + dNdr_norm = lambda r: dNdr(r) / integral + + N = integrate.quad( + dNdr_norm, 0.0, self._spatial_distribution.r_max + )[0] + + else: + # integrate the population to determine the true number of + # objects + N = integrate.quad(dNdr, 0.0, self._spatial_distribution.r_max)[ + 0 + ] + + log.info(f"The volume integral is {N}") # this should be poisson distributed n: int = np.random.poisson(N) @@ -961,7 +1001,13 @@ def draw_survey( n: int = n_samples - self._spatial_distribution.draw_distance(size=n) + probability = np.ones(n) + + self._spatial_distribution.draw_distance(size=n, normalize=normalize) + + # set the draw probability + + probability *= self._spatial_distribution.probability # now draw the sky positions @@ -1085,6 +1131,10 @@ def draw_survey( size=n ) # type: np.ndarray + self._luminosity_distribution.compute_probability(luminosities) + + probability *= self._luminosity_distribution.probability + # store the truths from the luminosity distribution truth[ self.luminosity_distribution.name @@ -1332,6 +1382,7 @@ def draw_survey( theta=self._spatial_distribution.theta, phi=self._spatial_distribution.phi, pop_synth=self.to_dict(), + probability=probability, ) def display(self) -> None: diff --git a/popsynth/utils/meta.py b/popsynth/utils/meta.py index 68ea435..c89f159 100644 --- a/popsynth/utils/meta.py +++ b/popsynth/utils/meta.py @@ -12,6 +12,7 @@ def __init__( vmin: Optional[float] = None, vmax: Optional[float] = None, free: bool = True, + is_normalization: bool = False, ): """ Parameter base class. @@ -25,18 +26,27 @@ def __init__( :param free: If `True`, parameter is free """ - self.name = None # type: Union[None, str] - self._vmin = vmin # type: Optional[float] - self._vmax = vmax # type: Optional[float] - self._default = default # type: Optional[float] - self._free = free + self.name: Optional[float] = None + self._vmin: Optional[float] = vmin # type: Optional[float] + self._vmax: Optional[float] = vmax # type: Optional[float] + self._default: Optional[float] = default # type: Optional[float] + self._free: bool = free + self._is_normalization: bool = is_normalization @property def default(self) -> Optional[float]: return self._default + @property + def is_normalization(self) -> bool: + return self._is_normalization + def __get__(self, obj, type=None) -> object: + if self._is_normalization: + + obj._normalization_parameter = self.name + try: if obj._parameter_storage[self.name] is None: @@ -62,6 +72,10 @@ def __get__(self, obj, type=None) -> object: def __set__(self, obj, value) -> None: self._is_set = True + if self._is_normalization: + + obj._normalization_parameter = self.name + if not self._free: log.error(f"{self.name} is fixed and cannot be set") @@ -142,11 +156,21 @@ def __new__(mcls, name, bases, attrs, **kwargs): cls._parameter_storage = {} + norm_counter = 0 + + norm = None + for k, v in attrs.items(): if isinstance(v, Parameter): v.name = k + if v.is_normalization: + + if norm_counter > 0: + + raise AssertionError("Can only have 1 normalization") + return cls def __subclasscheck__(cls, subclass): From a651732e0e9b7a7eb0a76a23b071ff023a2f7c3c Mon Sep 17 00:00:00 2001 From: grburgess Date: Tue, 3 May 2022 11:37:36 +0200 Subject: [PATCH 4/9] fixed normalization --- popsynth/distribution.py | 14 ++++++++++++++ popsynth/distributions/bpl_distribution.py | 10 +++++++--- popsynth/population_synth.py | 4 ++-- 3 files changed, 23 insertions(+), 5 deletions(-) diff --git a/popsynth/distribution.py b/popsynth/distribution.py index 99dd135..c462f54 100644 --- a/popsynth/distribution.py +++ b/popsynth/distribution.py @@ -379,9 +379,23 @@ def draw_distance(self, size: int, normalize: bool = False) -> None: # 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" diff --git a/popsynth/distributions/bpl_distribution.py b/popsynth/distributions/bpl_distribution.py index 8f32f50..67da3b5 100644 --- a/popsynth/distributions/bpl_distribution.py +++ b/popsynth/distributions/bpl_distribution.py @@ -109,20 +109,24 @@ def bpl(x, x0, x1, x2, a1, a2): # creatre a holder for the values x = np.atleast_1d(x) - out = np.empty_like(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 np.ndarray: def draw_survey( self, flux_sigma: Optional[float] = None, - log10_flux_rdraw: bool = True, + log10_flux_draw: bool = True, n_samples: Optional[int] = None, normalize: bool = False, ) -> Population: @@ -950,7 +950,7 @@ def draw_survey( if normalize: log.info("The population is being normalized such that") - log.info("the integral over N * dV/dz = N") + log.info("the integral over N * f(z) * dV/dz = N") old_value = None From 9f153d40e4a341bfacf31cb814df5dee6453885d Mon Sep 17 00:00:00 2001 From: grburgess Date: Mon, 13 Jun 2022 08:51:14 +0200 Subject: [PATCH 5/9] fix flake8 --- .flake8 | 2 +- .../aux_samplers/viewing_angle_sampler.py | 6 ++ popsynth/auxiliary_sampler.py | 72 ++++++++++++++----- popsynth/population_synth.py | 23 ++---- 4 files changed, 66 insertions(+), 37 deletions(-) diff --git a/.flake8 b/.flake8 index 42cf9ff..c2ff51b 100644 --- a/.flake8 +++ b/.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 diff --git a/popsynth/aux_samplers/viewing_angle_sampler.py b/popsynth/aux_samplers/viewing_angle_sampler.py index 61bc189..7468821 100644 --- a/popsynth/aux_samplers/viewing_angle_sampler.py +++ b/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)) + ) diff --git a/popsynth/auxiliary_sampler.py b/popsynth/auxiliary_sampler.py index c19176c..6d39331 100644 --- a/popsynth/auxiliary_sampler.py +++ b/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,6 +41,9 @@ def __init__( :type obs_values: ArrayLike :param selection: :type selection: ArrayLike + :param probability: + :type probability: ArrayLike + :returns: :returns: """ @@ -47,7 +51,7 @@ def __init__( 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, ) ) @@ -673,8 +711,6 @@ def true_sampler(self, size: int = 1): def _probability(self, true_values: np.ndarray) -> np.ndarray: pass - - def observation_sampler(self, size: int = 1) -> np.ndarray: return self._true_values diff --git a/popsynth/population_synth.py b/popsynth/population_synth.py index d3a1256..1bdd844 100644 --- a/popsynth/population_synth.py +++ b/popsynth/population_synth.py @@ -1083,6 +1083,7 @@ def draw_survey( self._derived_luminosity_sampler.true_values, self._derived_luminosity_sampler.obs_values, self._derived_luminosity_sampler.selector, + self._derived_luminosity_sampler.probability, ) ) @@ -1191,27 +1192,10 @@ def draw_survey( auxiliary_quantities += v.get_secondary_properties() - # # append these values to a dict - # auxiliary_quantities[k] = { - # "true_values": v.true_values, - # "obs_values": v.obs_values, - # "selection": v.selector, - # } # type: dict - # collect the secondary values for k2, v2 in v.secondary_samplers.items(): - # first we tell the sampler to go and retrieve all of - # its own secondaries - - # properties = v2.get_secondary_properties() # type: dict - - # for k3, v3 in properties.items(): - - # # now attach them - # auxiliary_quantities[k3] = v3 - # store the secondary truths truth[v2.name] = v2.truth @@ -1273,7 +1257,8 @@ def draw_survey( # selection = self._flux_selector.selection - # now apply the selection from the auxilary samplers + # now apply the selection and draw probability + # from the auxilary samplers for k, v in auxiliary_quantities.items(): @@ -1285,6 +1270,8 @@ def draw_survey( continue + probability *= v["probability"] + auxiliary_selection += v["selection"] log.info( From 713e1051f62a56f8693dfc02f16c37a25578f49e Mon Sep 17 00:00:00 2001 From: grburgess Date: Mon, 13 Jun 2022 08:57:31 +0200 Subject: [PATCH 6/9] blacken --- popsynth/distribution.py | 7 ------- popsynth/distributions/bpl_distribution.py | 3 +-- popsynth/distributions/log10_normal_distribution.py | 2 +- popsynth/distributions/log_normal_distribution.py | 2 +- popsynth/distributions/pareto_distribution.py | 2 +- popsynth/distributions/schechter_distribution.py | 2 +- popsynth/distributions/spherical_distribution.py | 2 +- popsynth/distributions/spiral_galaxy_distribution.py | 2 +- popsynth/test/test_populations.py | 2 +- 9 files changed, 8 insertions(+), 16 deletions(-) diff --git a/popsynth/distribution.py b/popsynth/distribution.py index c462f54..e0afb97 100644 --- a/popsynth/distribution.py +++ b/popsynth/distribution.py @@ -325,8 +325,6 @@ def draw_distance(self, size: int, normalize: bool = False) -> 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 @@ -384,18 +382,13 @@ def draw_distance(self, size: int, normalize: bool = False) -> 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" diff --git a/popsynth/distributions/bpl_distribution.py b/popsynth/distributions/bpl_distribution.py index 67da3b5..48ed39b 100644 --- a/popsynth/distributions/bpl_distribution.py +++ b/popsynth/distributions/bpl_distribution.py @@ -123,11 +123,10 @@ def bpl(x, x0, x1, x2, a1, a2): # compute the upper power law - idx = (x>=x1) & (x= x1) & (x < x2) out[idx] = np.power(x[idx], a2) * np.power(x1, a1 - a2) - return out * norm diff --git a/popsynth/distributions/log10_normal_distribution.py b/popsynth/distributions/log10_normal_distribution.py index ed5947a..845be51 100644 --- a/popsynth/distributions/log10_normal_distribution.py +++ b/popsynth/distributions/log10_normal_distribution.py @@ -33,7 +33,7 @@ def __init__(self, seed: int = 1234, name: str = "log10norm"): def phi(self, L): return (1.0 / (self.tau * L * np.sqrt(2 * np.pi))) * np.exp( - -((np.log10(L) - self.mu) ** 2) / (2 * self.tau ** 2) + -((np.log10(L) - self.mu) ** 2) / (2 * self.tau**2) ) def draw_luminosity(self, size=1): diff --git a/popsynth/distributions/log_normal_distribution.py b/popsynth/distributions/log_normal_distribution.py index 55dfe80..1cb783a 100644 --- a/popsynth/distributions/log_normal_distribution.py +++ b/popsynth/distributions/log_normal_distribution.py @@ -35,7 +35,7 @@ def __init__(self, seed: int = 1234, name: str = "lognorm"): def phi(self, L): return (1.0 / (self.tau * L * np.sqrt(2 * np.pi))) * np.exp( - -((np.log(L) - self.mu) ** 2) / (2 * self.tau ** 2) + -((np.log(L) - self.mu) ** 2) / (2 * self.tau**2) ) def draw_luminosity(self, size=1): diff --git a/popsynth/distributions/pareto_distribution.py b/popsynth/distributions/pareto_distribution.py index 7965729..762ee4d 100644 --- a/popsynth/distributions/pareto_distribution.py +++ b/popsynth/distributions/pareto_distribution.py @@ -37,7 +37,7 @@ def phi(self, L): idx = L >= self.Lmin out[idx] = ( - self.alpha * self.Lmin ** self.alpha / L[idx] ** (self.alpha + 1) + self.alpha * self.Lmin**self.alpha / L[idx] ** (self.alpha + 1) ) return out diff --git a/popsynth/distributions/schechter_distribution.py b/popsynth/distributions/schechter_distribution.py index 90b02d3..7e84e76 100644 --- a/popsynth/distributions/schechter_distribution.py +++ b/popsynth/distributions/schechter_distribution.py @@ -37,7 +37,7 @@ def __init__(self, seed: int = 1234, name: str = "schechter"): def phi(self, L): return ( - L ** self.alpha + L**self.alpha * np.exp(-L / self.Lmin) / (self.Lmin ** (1 + self.alpha) * sf.gamma(1 + self.alpha)) ) diff --git a/popsynth/distributions/spherical_distribution.py b/popsynth/distributions/spherical_distribution.py index 9f6eb16..0396d61 100644 --- a/popsynth/distributions/spherical_distribution.py +++ b/popsynth/distributions/spherical_distribution.py @@ -43,7 +43,7 @@ def transform(self, L, r): class ConstantSphericalDistribution(SphericalDistribution): _distribution_name = "ConstantSphericalDistribution" - Lambda = DistributionParameter(default=1, vmin=0 , is_normalization=True) + Lambda = DistributionParameter(default=1, vmin=0, is_normalization=True) def __init__( self, diff --git a/popsynth/distributions/spiral_galaxy_distribution.py b/popsynth/distributions/spiral_galaxy_distribution.py index 1e56b86..7f0db1a 100644 --- a/popsynth/distributions/spiral_galaxy_distribution.py +++ b/popsynth/distributions/spiral_galaxy_distribution.py @@ -99,7 +99,7 @@ def draw_sky_positions(self, size): 0, scale=0.07 * np.abs(self._distances), size=size ) - phi = np.arccos(zpos / np.sqrt(self._distances ** 2 + zpos ** 2)) + phi = np.arccos(zpos / np.sqrt(self._distances**2 + zpos**2)) self._theta = spiraltheta diff --git a/popsynth/test/test_populations.py b/popsynth/test/test_populations.py index 438f9a9..54d8797 100644 --- a/popsynth/test/test_populations.py +++ b/popsynth/test/test_populations.py @@ -76,7 +76,7 @@ def true_sampler(self, size): dl = cosmology.luminosity_distance(self._distance) - fluxes = self._luminosity / (4 * np.pi * dl ** 2) + fluxes = self._luminosity / (4 * np.pi * dl**2) self._true_values = fluxes From fa0524e1dd95fad9f919b84c12db6ee690ec6142 Mon Sep 17 00:00:00 2001 From: grburgess Date: Mon, 13 Jun 2022 09:00:50 +0200 Subject: [PATCH 7/9] update wheel --- .github/workflows/build_test.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/build_test.yml b/.github/workflows/build_test.yml index 2132939..4183947 100644 --- a/.github/workflows/build_test.yml +++ b/.github/workflows/build_test.yml @@ -24,6 +24,7 @@ jobs: python-version: ${{ matrix.python-version }} - name: Install dependencies run: | + python -m pip install --upgrade pip wheel pip install --upgrade flake8 numba numpy ipython scipy h5py matplotlib hypothesis pystan betagen coverage pytest pytest-cov cython codecov python setup.py install - name: Lint with flake8 From 5965d145fa7f9e8b00f6aaa29fa799282cfcd989 Mon Sep 17 00:00:00 2001 From: grburgess Date: Mon, 13 Jun 2022 09:08:37 +0200 Subject: [PATCH 8/9] add weird package --- .github/workflows/build_test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build_test.yml b/.github/workflows/build_test.yml index 4183947..2b8d414 100644 --- a/.github/workflows/build_test.yml +++ b/.github/workflows/build_test.yml @@ -25,7 +25,7 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip wheel - pip install --upgrade flake8 numba numpy ipython scipy h5py matplotlib hypothesis pystan betagen coverage pytest pytest-cov cython codecov + 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: | From 11247407c2cc143466bab2fed9ffd1798e166889 Mon Sep 17 00:00:00 2001 From: grburgess Date: Mon, 13 Jun 2022 09:26:48 +0200 Subject: [PATCH 9/9] fix error checking --- popsynth/population.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/popsynth/population.py b/popsynth/population.py index dc4b419..9e7217a 100644 --- a/popsynth/population.py +++ b/popsynth/population.py @@ -754,7 +754,7 @@ def _loadfrom(cls, f): known_distance_idx = (f["known_distance_idx"][()]).astype(int) unknown_distance_idx = (f["unknown_distance_idx"][()]).astype(int) - except AttributeError: + except KeyError: known_distances = None known_distance_idx = None @@ -764,7 +764,7 @@ def _loadfrom(cls, f): probability = f["probability"][()] - except AttributeError: + except KeyError: probability = None