`UserDefinedPrior` fails in mcmc due to clash multiprocess and pickle

[hposborn]

UserDefinedPrior requires defining a function which can be included in the output. In my case, I wanted to use the Kipping Beta prior on eccentricity:

def clipped_beta_func(params):
    #a=0.867; b=3.03; (gamma(a) * gamma(b) / gamma(a + b))=0.4271856369315835
    return (params[0]>0.001)&(params[0]<0.999),(params[0]**(0.867-1) * (1-params[0])**(3.03-1)) / 0.4271856369315835

priors += [radvel.prior.UserDefinedPrior(['e1'], clipped_beta_func, "Kipping (2013) Beta prior on eccentricity")]

However, when running from the command line, this seems to always fail for me at the mcmc stage, giving the error: _pickle.PicklingError: Can't pickle <function clipped_beta_func at 0x11ee86670>: it's not the same object as fitHD5457_beta.clipped_beta_func

I believe this is because multiprocessing requires pickled objects and functions cannot be explicitly pickled (unless they are natively/initially defined before multiprocessing). I guess this requires a bit more testing to prove.

I would like to suggest functionality to either run MCMC using user-defined priors (though I'm not sure how that would work) or to natively include known eccentricity priors (Kipping, Van Eylen, etc), although I know that cleanly sampling e-w space is less "clean" than with ecosw-esinw.

[hposborn]

The solution was to put any UserDefinedPrior directly into the radvel/prior.py file.

class KippingEccPrior(Prior):
    """Kipping beta distribtion prior on eccentricity. 
       Also includes HardBounds at 0.001/0.999 to avoid getting close to unphysical eccentricities.

    Args:
        num_planets (int): Number of planets. Used to set eccentricity prior for each planet
    """

    def __repr__(self):
        return "ecc constrained via Beta distribution from Kipping 2013"

    def __str__(self):
        return "$e$ constrained using $\Beta$ dist"

    def __init__(self, num_planets):
        self.num_planets = num_planets

    def __call__(self, params, vector):
        def _getpar(key, num_planet):
            return vector.vector[vector.indices['{}{}'.format(key, num_planet)]][0]

        p=0
        for num_planet in range(1, self.num_planets+1):
            try:
                e = _getpar('e', num_planet)
            except KeyError:
                e = _getpar('secosw', num_planet)**2 + _getpar('sesinw', num_planet)**2
            if (e>0.001)&(e<0.999):
                p += (e**(0.867-1) * (1-e)**(3.03-1)) / 0.4271856369315835 #(gamma(a) * gamma(b) / gamma(a + b))
            else:
                p += -np.inf
        return p

Though I'm not convinced this is the correct implementation.