HealPixPrior (Resloves Issue #419)

Hi all, This MR introduces a HealPixMapPriorDist object to be used to create a 2D or 3D joint prior in ra/dec (ra/dec/distance) according to a passed in HealPix map or .fits file to address and resolve issue #419 (closed).

This PriorObject inherits from BaseJointPriorDist and needs to overwrite three methods (_rescale(), _sample(), and _lnprob()) To illustrate this is working I sample 5000 samples from the prior and show below that these samples are distribvuted according to whats in the HealPix Map .fits file.

Here I show that the _rescale() method is also working by generating 5000 samples uniformly on [0,1] then rescale to the HealPixPrior which also matches with what is in the map file.

bilby/core/prior.py

@@ -641,22 +641,16 @@ class ConditionalPriorDict(PriorDict):
         self._check_resolved()
         self._update_rescale_keys(keys)
         result = dict()
-        for key, index in zip(self._rescale_keys, self._rescale_indexes):
+        for key, index in zip(self.sorted_keys_without_fixed_parameters, self._rescale_indexes):
             required_variables = {k: result[k] for k in getattr(self[key], 'required_variables', [])}
             result[key] = self[key].rescale(theta[index], **required_variables)
         return [result[key] for key in keys]
 
     def _update_rescale_keys(self, keys):
         if not keys == self._least_recently_rescaled_keys:
-            self._set_rescale_keys_and_indexes(keys)
+            self._rescale_indexes = [keys.index(element) for element in self.sorted_keys_without_fixed_parameters]
             self._least_recently_rescaled_keys = keys
 
-    def _set_rescale_keys_and_indexes(self, keys):
-        unconditional_keys, unconditional_idxs, _ = np.intersect1d(keys, self.unconditional_keys, return_indices=True)
-        conditional_keys, conditional_idxs, _ = np.intersect1d(keys, self.conditional_keys, return_indices=True)
-        self._rescale_keys = np.append(unconditional_keys, conditional_keys)
-        self._rescale_indexes = np.append(unconditional_idxs, conditional_idxs)
-
     def _check_resolved(self):
         if not self._resolved:
             raise IllegalConditionsException("The current set of priors contains unresolveable conditions.")
@@ -673,6 +667,10 @@ class ConditionalPriorDict(PriorDict):
     def sorted_keys(self):
         return self.unconditional_keys + self.conditional_keys
 
+    @property
+    def sorted_keys_without_fixed_parameters(self):
+        return [key for key in self.sorted_keys if not isinstance(self[key], (DeltaFunction, Constraint))]
+
     def __setitem__(self, key, value):
         super(ConditionalPriorDict, self).__setitem__(key, value)
         self._resolve_conditions()
@@ -3175,106 +3173,124 @@ class MultivariateGaussianDist(BaseJointPriorDist):
         Add a new mode.
         """
 
-        # add means
-        if mus is not None:
-            try:
-                self.mus.append(list(mus))  # means
-            except TypeError:
-                raise TypeError("'mus' must be a list")
-        else:
-            self.mus.append(np.zeros(self.num_vars))
+    def __len__(self):
+        return len(self.names)
 
-        # add the covariances if supplied
-        if cov is not None:
-            self.covs.append(np.asarray(cov))
+    def __repr__(self):
+        """Overrides the special method __repr__.
 
-            if len(self.covs[-1].shape) != 2:
-                raise ValueError("Covariance matrix must be a 2d array")
+        Returns a representation of this instance that resembles how it is instantiated.
+        Works correctly for all child classes
 
-            if (self.covs[-1].shape[0] != self.covs[-1].shape[1] or
-                    self.covs[-1].shape[0] != self.num_vars):
-                raise ValueError("Covariance shape is inconsistent")
+        Returns
+        -------
+        str: A string representation of this instance
 
-            # check matrix is symmetric
-            if not np.allclose(self.covs[-1], self.covs[-1].T):
-                raise ValueError("Covariance matrix is not symmetric")
+        """
+        dist_name = self.__class__.__name__
+        instantiation_dict = self.get_instantiation_dict()
+        args = ', '.join(['{}={}'.format(key, repr(instantiation_dict[key]))
+                          for key in instantiation_dict])
+        return "{}({})".format(dist_name, args)
 
-            self.sigmas.append(np.sqrt(np.diag(self.covs[-1])))  # standard deviations
+    def prob(self, samp):
+        """
+        Get the probability of a sample. For bounded priors the
+        probability will not be properly normalised.
+        """
 
-            # convert covariance into a correlation coefficient matrix
-            D = self.sigmas[-1] * np.identity(self.covs[-1].shape[0])
-            Dinv = np.linalg.inv(D)
-            self.corrcoefs.append(np.dot(np.dot(Dinv, self.covs[-1]), Dinv))
-        elif corrcoef is not None and sigmas is not None:
-            self.corrcoefs.append(np.asarray(corrcoef))
+        return np.exp(self.ln_prob(samp))
 
-            if len(self.corrcoefs[-1].shape) != 2:
-                raise ValueError("Correlation coefficient matrix must be a 2d "
-                                 "array.")
+    def _check_samp(self, value):
+        """
+        Get the log-probability of a sample. For bounded priors the
+        probability will not be properly normalised.
 
-            if (self.corrcoefs[-1].shape[0] != self.corrcoefs[-1].shape[1] or
-                    self.corrcoefs[-1].shape[0] != self.num_vars):
-                raise ValueError("Correlation coefficient matrix shape is "
-                                 "inconsistent")
+        Parameters
+        ----------
+        value: array_like
+            A 1d vector of the sample, or 2d array of sample values with shape
+            NxM, where N is the number of samples and M is the number of
+            parameters.
 
-            # check matrix is symmetric
-            if not np.allclose(self.corrcoefs[-1], self.corrcoefs[-1].T):
-                raise ValueError("Correlation coefficient matrix is not "
-                                 "symmetric")
+        Returns
+        -------
+        samp: array_like
+            returns the input value as a sample array
+        outbounds: array_like
+            Boolean Array that selects samples in samp that are out of given bounds
+        """
+        samp = np.array(value)
+        if len(samp.shape) == 1:
+            samp = samp.reshape(1, self.num_vars)
 
-            # check diagonal is all ones
-            if not np.all(np.diag(self.corrcoefs[-1]) == 1.):
-                raise ValueError("Correlation coefficient matrix is not"
-                                 "correct")
+        if len(samp.shape) != 2:
+            raise ValueError("Array is the wrong shape")
+        elif samp.shape[1] != self.num_vars:
+            raise ValueError("Array is the wrong shape")
 
-            try:
-                self.sigmas.append(list(sigmas))  # standard deviations
-            except TypeError:
-                raise TypeError("'sigmas' must be a list")
+        # check sample(s) is within bounds
+        outbounds = np.ones(samp.shape[0], dtype=np.bool)
+        for s, bound in zip(samp.T, self.bounds.values()):
+            outbounds = (s < bound[0]) | (s > bound[1])
+            if np.any(outbounds):
+                break
+        return samp, outbounds
 
-            if len(self.sigmas[-1]) != self.num_vars:
-                raise ValueError("Number of standard deviations must be the "
-                                 "same as the number of parameters.")
+    def ln_prob(self, value):
+        """
+        Get the log-probability of a sample. For bounded priors the
+        probability will not be properly normalised.
 
-            # convert correlation coefficients to covariance matrix
-            D = self.sigmas[-1] * np.identity(self.corrcoefs[-1].shape[0])
-            self.covs.append(np.dot(D, np.dot(self.corrcoefs[-1], D)))
-        else:
-            # set unit variance uncorrelated covariance
-            self.corrcoefs.append(np.eye(self.num_vars))
-            self.covs.append(np.eye(self.num_vars))
-            self.sigmas.append(np.ones(self.num_vars))
+        Parameters
+        ----------
+        value: array_like
+            A 1d vector of the sample, or 2d array of sample values with shape
+            NxM, where N is the number of samples and M is the number of
+            parameters.
+        """
 
-        # get eigen values and vectors
-        try:
-            evals, evecs = np.linalg.eig(self.corrcoefs[-1])
-            self.eigvalues.append(evals)
-            self.eigvectors.append(evecs)
-        except Exception as e:
-            raise RuntimeError("Problem getting eigenvalues and vectors: "
-                               "{}".format(e))
-
-        # check eigenvalues are positive
-        if np.any(self.eigvalues[-1] <= 0.):
-            raise ValueError("Correlation coefficient matrix is not positive "
-                             "definite")
-        self.sqeigvalues.append(np.sqrt(self.eigvalues[-1]))
-
-        # set the weights
-        if weight is None:
-            self.weights.append(1.)
+        samp, outbounds = self._check_samp(value)
+        lnprob = -np.inf * np.ones(samp.shape[0])
+        lnprob = self._ln_prob(samp, lnprob, outbounds)
+        if samp.shape[0] == 1:
+            return lnprob[0]
         else:
-            self.weights.append(weight)
+            return lnprob
+
+    def _ln_prob(self, samp, lnprob, outbounds):
+        """
+        Get the log-probability of a sample. For bounded priors the
+        probability will not be properly normalised. **this method needs overwritten by child class**
 
-        # set the cumulative relative weights
-        self.cumweights = np.cumsum(self.weights) / np.sum(self.weights)
+        Parameters
+        ----------
+        samp: vector
+            sample to evaluate the ln_prob at
+        lnprob: vector
+            of -inf pased in with the same shape as the number of samples
+        outbounds: array_like
+            boolean array showing which samples in lnprob vector are out of the given bounds
 
-        # add the mode
-        self.nmodes += 1
+        Returns
+        -------
+        lnprob: vector
+            array of lnprob values for each sample given
+        """
+        """
+        Here is where the subclass where overwrite ln_prob method
+        """
+        return lnprob
 
-        # add multivariate Gaussian
-        self.mvn.append(scipy.stats.multivariate_normal(mean=self.mus[-1],
-                                                        cov=self.covs[-1]))
+    def sample(self, size=1, **kwargs):
+        """
+        Draw, and set, a sample from the Dist, accompanying method _sample needs to overwritten
+
+        Parameters
+        ----------
+        size: int
+            number of samples to generate, defualts to 1
+        """
 
     def _rescale(self, samp, **kwargs):
         try:
@@ -3613,6 +3629,7 @@ def conditional_prior_factory(prior_class):
             self.condition_func = condition_func
             self._reference_params = reference_params
             self.__class__.__name__ = 'Conditional{}'.format(prior_class.__name__)
+            self.__class__.__qualname__ = 'Conditional{}'.format(prior_class.__qualname__)
 
         def sample(self, size=None, **required_variables):
             """Draw a sample from the prior