Merge pull request #40 from FlorianPfaff/restructure/manifold_vs_type

Differentiate between manifold (e.g., torus) and type (e.g., mixture)

environment.yml

@@ -13,3 +13,4 @@ dependencies:
   - pip: # Installing them via pip is much faster than using the conda-forge channel
       - filterpy
       - healpy
+      - pyshtools

pyproject.toml

@@ -3,7 +3,7 @@ name = "pyrecest"
 description = "Framework for recursive Bayesian estimation in Python."
 readme = "README.md"
 authors = ["Florian Pfaff <pfaff@kit.edu>"]
-version = "0.1.2"
+version = "0.2.0"
 
 [tool.poetry.dependencies]
 python = "^3.9"
@@ -12,6 +12,7 @@ scipy = "*"
 matplotlib = "*"
 mpmath = "*"
 filterpy = "*"
+pyshtools = "*"
 
 [tool.poetry.extras]
 healpy_support = ["healpy"]

pyrecest/distributions/__init__.py

@@ -1,8 +1,7 @@
+from .abstract_custom_distribution import AbstractCustomDistribution
 from .abstract_dirac_distribution import AbstractDiracDistribution
 from .abstract_disk_distribution import AbstractDiskDistribution
-from .abstract_distribution import AbstractDistribution
 from .abstract_ellipsoidal_ball_distribution import AbstractEllipsoidalBallDistribution
-from .abstract_non_conditional_distribution import AbstractNonConditionalDistribution
 from .abstract_periodic_distribution import AbstractPeriodicDistribution
 from .abstract_uniform_distribution import AbstractUniformDistribution
 from .circle.abstract_circular_distribution import AbstractCircularDistribution
@@ -11,7 +10,6 @@
 from .circle.custom_circular_distribution import CustomCircularDistribution
 from .circle.von_mises_distribution import VonMisesDistribution
 from .circle.wrapped_normal_distribution import WrappedNormalDistribution
-from .custom_distribution import CustomDistribution
 from .disk_uniform_distribution import DiskUniformDistribution
 from .ellipsoidal_ball_uniform_distribution import EllipsoidalBallUniformDistribution
 from .hypersphere_subset.abstract_hemispherical_distribution import (
@@ -83,11 +81,9 @@
     "CustomLinearDistribution",
     "HyperhemisphericalWatsonDistribution",
     "AbstractLinearDistribution",
-    "AbstractNonConditionalDistribution",
     "AbstractCircularDistribution",
     "AbstractDiracDistribution",
     "AbstractDiskDistribution",
-    "AbstractDistribution",
     "AbstractEllipsoidalBallDistribution",
     "AbstractHyperhemisphericalDistribution",
     "AbstractHypersphereSubsetDistribution",
@@ -98,7 +94,7 @@
     "AbstractToroidalDistribution",
     "AbstractUniformDistribution",
     "BinghamDistribution",
-    "CustomDistribution",
+    "AbstractCustomDistribution",
     "CustomHemisphericalDistribution",
     "CustomHyperhemisphericalDistribution",
     "DiskUniformDistribution",

pyrecest/distributions/abstract_bounded_domain_distribution.py

@@ -1,5 +1,7 @@
-from .abstract_non_conditional_distribution import AbstractNonConditionalDistribution
+from .abstract_manifold_specific_distribution import (
+    AbstractManifoldSpecificDistribution,
+)
 
 
-class AbstractBoundedDomainDistribution(AbstractNonConditionalDistribution):
+class AbstractBoundedDomainDistribution(AbstractManifoldSpecificDistribution):
     pass

pyrecest/distributions/abstract_custom_distribution.py

@@ -0,0 +1,30 @@
+import copy
+import warnings
+from abc import abstractmethod
+
+from .abstract_distribution_type import AbstractDistributionType
+
+
+class AbstractCustomDistribution(AbstractDistributionType):
+    def __init__(self, f, scale_by=1):
+        self.f = f
+        self.scale_by = scale_by
+
+    def pdf(self, xs):
+        # Shifting is something for subclasses
+        return self.scale_by * self.f(xs)
+
+    @abstractmethod
+    def integrate(self, integration_boundaries=None):
+        pass
+
+    def normalize(self, verify=None):
+        cd = copy.deepcopy(self)
+
+        integral = self.integrate()
+        cd.scale_by = cd.scale_by / integral
+
+        if verify and abs(cd.integrate()[0] - 1) > 0.001:
+            warnings.warn("Density is not yet properly normalized.", UserWarning)
+
+        return cd

pyrecest/distributions/abstract_custom_nonperiodic_distribution.py

@@ -0,0 +1,8 @@
+from .abstract_custom_distribution import AbstractCustomDistribution
+from .abstract_nonperiodic_distribution import AbstractNonperiodicDistribution
+
+
+class AbstractCustomNonPeriodicDistribution(
+    AbstractCustomDistribution, AbstractNonperiodicDistribution
+):
+    pass

pyrecest/distributions/abstract_dirac_distribution.py

@@ -3,16 +3,11 @@
 
 import numpy as np
 
-from .abstract_distribution import AbstractDistribution
+from .abstract_distribution_type import AbstractDistributionType
 
 
-class AbstractDiracDistribution(AbstractDistribution):
+class AbstractDiracDistribution(AbstractDistributionType):
     def __init__(self, d, w=None):
-        AbstractDistribution.__init__(self, dim=d.shape[-1])
-        if self.dim > d.shape[0]:
-            warnings.warn(
-                "Not even one Dirac per dimension. If this warning is unexpected, verify d_ is shaped correctly."
-            )
         if w is None:
             w = np.ones(d.shape[0]) / d.shape[0]
 
@@ -31,11 +26,12 @@ def normalize(self):
         dist.normalize_in_place()
         return dist
 
-    def apply_function(self, f):
-        d_new = np.apply_along_axis(f, 1, self.d)
+    def apply_function(self, f, f_supports_multiple=True):
         dist = copy.deepcopy(self)
-        dist.d = d_new
-        dist.w = self.w
+        if f_supports_multiple:
+            dist.d = f(dist.d)
+        else:
+            dist.d = np.apply_along_axis(f, 1, dist.d)
         return dist
 
     def reweigh(self, f):
@@ -56,7 +52,7 @@ def reweigh(self, f):
 
     def sample(self, n):
         ids = np.random.choice(self.w.size, size=n, p=self.w)
-        return self.d[ids, :]
+        return self.d[ids] if self.d.ndim == 1 else self.d[ids, :]
 
     def entropy(self):
         warnings.warn("Entropy is not defined in a continuous sense")
@@ -71,7 +67,7 @@ def integrate(self, left=None, right=None):
     def log_likelihood(self, *args):
         raise NotImplementedError("PDF:UNDEFINED, not supported")
 
-    def pdf(self, *args):
+    def pdf(self, xs):
         raise NotImplementedError("PDF:UNDEFINED, pdf is not defined")
 
     def integrate_numerically(self, *args):
@@ -102,3 +98,17 @@ def mode_numerical(self, _=None):
 
     def entropy_numerical(self):
         raise NotImplementedError("PDF:UNDEFINED, not supported")
+
+    @classmethod
+    def from_distribution(cls, distribution, n_samples):
+        is_valid_class = False
+        # Look if distribution is of the correct type
+        for base in cls.__bases__:
+            if isinstance(distribution, base):
+                is_valid_class = True
+
+        assert is_valid_class
+        assert isinstance(n_samples, int) and n_samples > 0
+        samples = distribution.sample(n_samples)
+        weights = np.ones(n_samples) / n_samples
+        return cls(samples, weights)

pyrecest/distributions/abstract_distribution_type.py

@@ -0,0 +1,2 @@
+class AbstractDistributionType:
+    pass

pyrecest/distributions/abstract_distribution.py → pyrecest/distributions/abstract_manifold_specific_distribution.py

@@ -3,12 +3,16 @@
 import numpy as np
 
 
-class AbstractDistribution(ABC):
+class AbstractManifoldSpecificDistribution(ABC):
     """Abstract base class for all distributions."""
 
     def __init__(self, dim=None):
         self._dim = dim
 
+    @abstractmethod
+    def get_manifold_size(self):
+        pass
+
     @property
     def dim(self):
         return self._dim

pyrecest/distributions/abstract_mixture.py

@@ -1,13 +1,14 @@
 import warnings
+from abc import abstractmethod
 
 import numpy as np
 
-from .abstract_distribution import AbstractDistribution
+from .abstract_distribution_type import AbstractDistributionType
 
 
-class AbstractMixture(AbstractDistribution):
+class AbstractMixture(AbstractDistributionType):
     def __init__(self, dists, w=None):
-        AbstractDistribution.__init__(self, dim=dists[0].dim)
+        AbstractDistributionType.__init__(self)
         if w is None:
             w = np.ones(len(dists)) / len(dists)
         else:
@@ -51,6 +52,11 @@ def sample(self, n):
 
         return s
 
+    @property
+    @abstractmethod
+    def input_dim(self):
+        pass
+
     def pdf(self, xs):
         assert xs.shape[-1] == self.input_dim, "Dimension mismatch"
 

pyrecest/distributions/abstract_nonperiodic_distribution.py

@@ -0,0 +1,7 @@
+from .abstract_manifold_specific_distribution import (
+    AbstractManifoldSpecificDistribution,
+)
+
+
+class AbstractNonperiodicDistribution(AbstractManifoldSpecificDistribution):
+    pass

pyrecest/distributions/abstract_orthogonal_basis_distribution.py

@@ -0,0 +1,40 @@
+import copy
+from abc import abstractmethod
+
+import numpy as np
+
+from .abstract_distribution_type import AbstractDistributionType
+
+
+class AbstractOrthogonalBasisDistribution(AbstractDistributionType):
+    def __init__(self, coeff_mat, transformation):
+        self.transformation = transformation
+        self.coeff_mat = coeff_mat
+        self.normalize_in_place()
+
+    @abstractmethod
+    def normalize_in_place(self):
+        pass
+
+    @abstractmethod
+    def value(self, xa):
+        pass
+
+    def normalize(self):
+        result = copy.deepcopy(self)
+        return result.normalize_in_place()
+
+    def pdf(self, xa):
+        val = self.value(xa)
+        if self.transformation == "sqrt":
+            assert all(np.imag(val) < 0.0001)
+            return np.real(val) ** 2
+
+        if self.transformation == "identity":
+            return val
+
+        if self.transformation == "log":
+            print("Warning: Density may not be normalized")
+            return np.exp(val)
+
+        raise ValueError("Transformation not recognized or unsupported")

pyrecest/distributions/abstract_uniform_distribution.py

@@ -2,13 +2,13 @@
 
 import numpy as np
 
-from .abstract_distribution import AbstractDistribution
+from .abstract_distribution_type import AbstractDistributionType
 
 
-class AbstractUniformDistribution(AbstractDistribution):
+class AbstractUniformDistribution(AbstractDistributionType):
     # Uniform distribution.
     def pdf(self, xs):
-        return 1 / self.get_manifold_size() * np.ones(np.size(xs) // self.dim)
+        return 1 / self.get_manifold_size() * np.ones(xs.shape[0])
 
     @abstractmethod
     def get_manifold_size(self):

pyrecest/distributions/circle/circular_dirac_distribution.py

@@ -6,8 +6,8 @@ class CircularDiracDistribution(
     HypertoroidalDiracDistribution, AbstractCircularDistribution
 ):
     def __init__(self, d, w=None):
-        HypertoroidalDiracDistribution.__init__(self, d, w)
         AbstractCircularDistribution.__init__(self)
+        HypertoroidalDiracDistribution.__init__(self, d, w, dim=1)
         assert self.d.shape[0] == self.w.shape[0]
 
     def plot(self):

pyrecest/distributions/circle/custom_circular_distribution.py

@@ -1,16 +1,24 @@
-from .abstract_circular_distribution import AbstractCircularDistribution
-from ..custom_distribution import CustomDistribution
 import numpy as np
 
-class CustomCircularDistribution(CustomDistribution, AbstractCircularDistribution):
-    def __init__(self, f_):
+from ..abstract_custom_distribution import AbstractCustomDistribution
+from .abstract_circular_distribution import AbstractCircularDistribution
+
+
+class CustomCircularDistribution(
+    AbstractCustomDistribution, AbstractCircularDistribution
+):
+    def __init__(self, f_, scale_by=1, shift_by=0):
         """
         It is the user's responsibility to ensure that f is a valid
         circular density, i.e., 2pi-periodic, nonnegative and
         normalized.
         """
         AbstractCircularDistribution.__init__(self)
-        CustomDistribution.__init__(self, f_, 1)
-
-    def integrate(self, integration_boundaries=np.array([0, 2*np.pi])):
+        AbstractCustomDistribution.__init__(self, f_, scale_by)
+        self.shift_by = shift_by
+
+    def pdf(self, xs):
+        super().pdf(np.mod(xs + self.shift_by, 2 * np.pi))
+
+    def integrate(self, integration_boundaries=np.array([0, 2 * np.pi])):
         return AbstractCircularDistribution.integrate(self, integration_boundaries)