class Sampler:
    """
    The basic sampler object that performs the actual nested sampling.

    Parameters
    ----------
    loglikelihood : function
        Function returning ln(likelihood) given parameters as a 1-d `~numpy`
        array of length `ndim`.

    prior_transform : function
        Function transforming a sample from the a unit cube to the parameter
        space of interest according to the prior.

    npdim : int, optional
        Number of parameters accepted by `prior_transform`.

    live_points : list of 3 `~numpy.ndarray` each with shape (nlive, ndim)
        Initial set of "live" points. Contains `live_u`, the coordinates
        on the unit cube, `live_v`, the transformed variables, and
        `live_logl`, the associated loglikelihoods.

    update_interval : int
        Only update the bounding distribution every `update_interval`-th
        likelihood call.

    first_update : dict
        A dictionary containing parameters governing when the sampler should
        first update the bounding distribution from the unit cube to the one
        specified by the user.

    rstate : `~numpy.random.Generator`
        `~numpy.random.Generator` instance.

    queue_size: int
        Carry out likelihood evaluations in parallel by queueing up new live
        point proposals using (at most) this many threads/members.

    pool: pool
        Use this pool of workers to execute operations in parallel.

    use_pool : dict, optional
        A dictionary containing flags indicating where the provided `pool`
        should be used to execute operations in parallel.

    """
    def __init__(self, loglikelihood, prior_transform, npdim, live_points,
                 update_interval, first_update, rstate, queue_size, pool,
                 use_pool, ncdim):

        # distributions
        self.loglikelihood = loglikelihood
        self.prior_transform = prior_transform
        self.npdim = npdim
        self.ncdim = ncdim

        # live points
        self.live_u, self.live_v, self.live_logl = live_points
        self.nlive = len(self.live_u)
        self.live_bound = np.zeros(self.nlive, dtype='int')
        self.live_it = np.zeros(self.nlive, dtype='int')

        # bounding updates
        self.update_interval = update_interval
        self.ubound_ncall = first_update.get('min_ncall', 2 * self.nlive)
        self.ubound_eff = first_update.get('min_eff', 10.)
        self.logl_first_update = None

        # random state
        self.rstate = rstate

        # set to none just for qa
        self.scale = None
        self.method = None
        self.kwargs = {}

        # parallelism
        self.pool = pool  # provided pool
        if self.pool is None:
            self.M = map
        else:
            self.M = pool.map
        self.use_pool = use_pool  # provided flags for when to use the pool
        self.use_pool_ptform = use_pool.get('prior_transform', True)
        self.use_pool_logl = use_pool.get('loglikelihood', True)
        self.use_pool_evolve = use_pool.get('propose_point', True)
        self.use_pool_update = use_pool.get('update_bound', True)
        if self.use_pool_evolve:
            self.queue_size = queue_size  # size of the queue
        else:
            self.queue_size = 1
        self.queue = []  # proposed live point queue
        self.nqueue = 0  # current size of the queue
        self.unused = 0  # total number of proposals unused
        self.used = 0  # total number of proposals used

        # sampling
        self.it = 1  # current iteration
        self.since_update = 0  # number of calls since the last update
        self.ncall = self.nlive  # number of function calls
        self.dlv = math.log((self.nlive + 1.) / self.nlive)  # shrinkage/iter
        self.bound = [UnitCube(self.ncdim)]  # bounding distributions
        self.nbound = 1  # total number of unique bounding distributions
        self.added_live = False  # whether leftover live points were used
        self.eff = 0.  # overall sampling efficiency
        self.cite = ''  # Default empty
        # results
        self.saved_run = RunRecord()

    def propose_point(self, *args):
        raise RuntimeError('Should be overriden')

    def evolve_point(self, *args):
        raise RuntimeError('Should be overriden')

    def update_proposal(self, *args):
        raise RuntimeError('Should be overriden')

    def update(self):
        raise RuntimeError('Should be overriden')

    def __getstate__(self):
        """Get state information for pickling."""

        state = self.__dict__.copy()

        # attempt to remove from internal sampler, if not dealt within
        # DynamicSampler class
        try:
            # remove random module
            del state['rstate']

            # deal with pool
            if state['pool'] is not None:
                del state['pool']  # remove pool
                del state['M']  # remove `pool.map` function hook
        except AttributeError:
            pass

        return state

    def reset(self):
        """Re-initialize the sampler."""

        # live points
        self.live_u = self.rstate.uniform(size=(self.nlive, self.npdim))
        if self.use_pool_ptform:
            # Use the pool to compute the prior transform.
            self.live_v = np.array(
                list(self.M(self.prior_transform, np.asarray(self.live_u))))
        else:
            # Compute the prior transform using the default `map` function.
            self.live_v = np.array(
                list(map(self.prior_transform, np.asarray(self.live_u))))
        self.live_logl = self.loglikelihood.map(np.asarray(self.live_v))

        self.live_bound = np.zeros(self.nlive, dtype='int')
        self.live_it = np.zeros(self.nlive, dtype='int')

        # parallelism
        self.queue = []
        self.nqueue = 0
        self.unused = 0
        self.used = 0

        # sampling
        self.it = 1
        self.since_update = 0
        self.ncall = self.nlive
        self.bound = [UnitCube(self.ncdim)]
        self.nbound = 1
        self.added_live = False

        # results
        self.saved_run = RunRecord()

    @property
    def results(self):
        """Saved results from the nested sampling run. If bounding
        distributions were saved, those are also returned."""

        d = {}
        for k in [
                'nc', 'v', 'id', 'it', 'u', 'logwt', 'logl', 'logvol', 'logz',
                'logzvar', 'h'
        ]:
            d[k] = np.array(self.saved_run.D[k])

        # Add all saved samples to the results.
        if self.save_samples:
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                results = [('nlive', self.nlive), ('niter', self.it - 1),
                           ('ncall', d['nc']), ('eff', self.eff),
                           ('samples', d['v'])]
                for k in ['id', 'it', 'u']:
                    results.append(('samples_' + k, d[k]))
                for k in ['logwt', 'logl', 'logvol', 'logz']:
                    results.append((k, d[k]))
                results.append(('logzerr', np.sqrt(d['logzvar'])))
                results.append(('information', d['h']))
        else:
            raise ValueError("You didn't save any samples!")

        # Add any saved bounds (and ancillary quantities) to the results.
        if self.save_bounds:
            results.append(('bound', copy.deepcopy(self.bound)))
            results.append(('bound_iter',
                            np.array(self.saved_run.D['bounditer'],
                                     dtype='int')))
            results.append(('samples_bound',
                            np.array(self.saved_run.D['boundidx'],
                                     dtype='int')))
            results.append(('scale', np.array(self.saved_run.D['scale'])))

        return Results(results)

    @property
    def n_effective(self):
        """
        Estimate the effective number of posterior samples using the Kish
        Effective Sample Size (ESS) where `ESS = sum(wts)^2 / sum(wts^2)`.
        Note that this is `len(wts)` when `wts` are uniform and
        `1` if there is only one non-zero element in `wts`.

        """

        if (len(self.saved_run.D['logwt']) == 0) or (np.max(
                self.saved_run.D['logwt']) > 0.01 * np.nan_to_num(-np.inf)):
            # If there are no saved weights, or its -inf return 0.
            return 0
        else:
            # Otherwise, compute Kish ESS.
            logwts = np.array(self.saved_run.D['logwt'])
            logneff = logsumexp(logwts) * 2 - logsumexp(logwts * 2)
            return np.exp(logneff)

    @property
    def citations(self):
        """
        Return list of papers that should be cited given the specified
        configuration of the sampler.

        """

        return self.cite

    def _beyond_unit_bound(self, loglstar):
        """Check whether we should update our bound beyond the initial
        unit cube."""

        if self.logl_first_update is None:
            # If we haven't already updated our bounds, check if we satisfy
            # the provided criteria for establishing the first bounding update.
            check = (self.ncall > self.ubound_ncall
                     and self.eff < self.ubound_eff)
            if check:
                # Save the log-likelihood where our first update took place.
                self.logl_first_update = loglstar
            return check
        else:
            # If we've already update our bounds, check if we've exceeded the
            # saved log-likelihood threshold. (This is useful when sampling
            # within `dynamicsampler`).
            return loglstar >= self.logl_first_update

    def _empty_queue(self):
        """Dump all live point proposals currently on the queue."""

        while True:
            try:
                # Remove unused points from the queue.
                self.queue.pop()
                self.unused += 1  # add to the total number of unused points
                self.nqueue -= 1
            except IndexError:
                # If the queue is empty, we're done!
                self.nqueue = 0
                break

    def _fill_queue(self, loglstar):
        """Sequentially add new live point proposals to the queue."""

        # Add/zip arguments to submit to the queue.
        point_queue = []
        axes_queue = []

        # All the samplers should have have a starting point
        # satisfying a strict logl>loglstar criterion
        # The slice sampler will just fail if it's not the case
        # therefore we provide those subsets of points to choose from.

        if self.method != 'unif':
            args = (np.nonzero(self.live_logl > loglstar)[0], )
            if len(args[0]) == 0:
                raise RuntimeError(
                    'No live points are above loglstar. '
                    'Do you have a likelihood plateau ? '
                    'It is also possible that you are trying to sample '
                    'excessively around the very peak of the posterior')
        else:
            args = ()
        while self.nqueue < self.queue_size:
            if self._beyond_unit_bound(loglstar):
                # Propose points using the provided sampling/bounding options.
                point, axes = self.propose_point(*args)
                evolve_point = self.evolve_point
            else:
                # Propose/evaluate points directly from the unit cube.
                point = self.rstate.uniform(size=self.npdim)
                axes = np.identity(self.ncdim)
                evolve_point = sample_unif
            point_queue.append(point)
            axes_queue.append(axes)
            self.nqueue += 1
        loglstars = [loglstar for i in range(self.queue_size)]
        scales = [self.scale for i in range(self.queue_size)]
        ptforms = [self.prior_transform for i in range(self.queue_size)]
        logls = [self.loglikelihood for i in range(self.queue_size)]
        kwargs = [self.kwargs for i in range(self.queue_size)]
        seeds = get_seed_sequence(self.rstate, self.queue_size)
        args = zip(point_queue, loglstars, axes_queue, scales, ptforms, logls,
                   seeds, kwargs)

        if self.use_pool_evolve:
            # Use the pool to propose ("evolve") a new live point.
            self.queue = list(self.M(evolve_point, args))
        else:
            # Propose ("evolve") a new live point using the default `map`
            # function.
            self.queue = list(map(evolve_point, args))

    def _get_point_value(self, loglstar):
        """Grab the first live point proposal in the queue."""

        # If the queue is empty, refill it.
        if self.nqueue <= 0:
            self._fill_queue(loglstar)

        # Grab the earliest entry.
        u, v, logl, nc, blob = self.queue.pop(0)

        self.used += 1  # add to the total number of used points
        self.nqueue -= 1

        return u, v, logl, nc, blob

    def _new_point(self, loglstar):
        """Propose points until a new point that satisfies the log-likelihood
        constraint `loglstar` is found."""

        ncall, nupdate = 0, 0
        while True:
            # Get the next point from the queue
            u, v, logl, nc, blob = self._get_point_value(loglstar)
            ncall += nc

            # Bounding checks.
            ucheck = ncall >= self.update_interval * (1 + nupdate)
            bcheck = self._beyond_unit_bound(loglstar)

            # If our queue is empty, update any tuning parameters associated
            # with our proposal (sampling) method.
            if blob is not None and self.nqueue <= 0 and bcheck:
                self.update_proposal(blob)

            # If we satisfy the log-likelihood constraint, we're done!
            if logl > loglstar:
                break

            # If there has been more than `update_interval` function calls
            # made *and* we satisfy the criteria for moving beyond sampling
            # from the unit cube, update the bound.
            if ucheck and bcheck:
                bound = self.update()
                if self.save_bounds:
                    self.bound.append(bound)
                self.nbound += 1
                nupdate += 1
                self.since_update = -ncall  # ncall will be added back later

        return u, v, logl, ncall

    def add_live_points(self):
        """Add the remaining set of live points to the current set of dead
        points. Instantiates a generator that will be called by
        the user. Returns the same outputs as :meth:`sample`."""

        # Check if the remaining live points have already been added
        # to the output set of samples.
        if self.added_live:
            raise ValueError("The remaining live points have already "
                             "been added to the list of samples!")
        else:
            self.added_live = True

        # After N samples have been taken out, the remaining volume is
        # `e^(-N / nlive)`. The remaining points are distributed uniformly
        # within the remaining volume so that the expected volume enclosed
        # by the `i`-th worst likelihood is
        # `e^(-N / nlive) * (nlive + 1 - i) / (nlive + 1)`.
        logvols = self.saved_run.D['logvol'][-1]
        logvols += np.log(1. - (np.arange(self.nlive) + 1.) /
                          (self.nlive + 1.))
        logvols_pad = np.concatenate(
            ([self.saved_run.D['logvol'][-1]], logvols))

        # Defining change in `logvol` used in `logzvar` approximation.
        dlvs = logvols_pad[:-1] - logvols_pad[1:]

        # Sorting remaining live points.
        lsort_idx = np.argsort(self.live_logl)
        loglmax = max(self.live_logl)

        # Grabbing relevant values from the last dead point.
        logz = self.saved_run.D['logz'][-1]
        logzvar = self.saved_run.D['logzvar'][-1]
        h = self.saved_run.D['h'][-1]
        loglstar = self.saved_run.D['logl'][-1]
        if self._beyond_unit_bound(loglstar):
            bounditer = self.nbound - 1
        else:
            bounditer = 0

        # Add contributions from the remaining live points in order
        # from the lowest to the highest log-likelihoods.
        for i in range(self.nlive):

            # Grab live point with `i`-th lowest log-likelihood along with
            # ancillary quantities.
            idx = lsort_idx[i]
            logvol, dlv = logvols[i], dlvs[i]
            # we are doing copies here, because live_u/live_v are
            # updated in place
            ustar = self.live_u[idx].copy()
            vstar = self.live_v[idx].copy()
            loglstar_new = self.live_logl[idx]
            boundidx = self.live_bound[idx]
            point_it = self.live_it[idx]

            (logwt, logz, logzvar,
             h) = progress_integration(loglstar, loglstar_new, logz, logzvar,
                                       logvol, dlv, h)
            loglstar = loglstar_new
            logz_remain = loglmax + logvol  # remaining ln(evidence)
            delta_logz = np.logaddexp(logz, logz_remain) - logz  # dlogz

            # Save results.
            if self.save_samples:
                self.saved_run.append(
                    dict(id=idx,
                         u=ustar,
                         v=vstar,
                         logl=loglstar,
                         logvol=logvol,
                         logwt=logwt,
                         logz=logz,
                         logzvar=logzvar,
                         h=h,
                         nc=1,
                         boundidx=boundidx,
                         it=point_it,
                         bounditer=bounditer,
                         scale=self.scale))
            self.eff = 100. * (self.it + i) / self.ncall  # efficiency

            # Return our new "dead" point and ancillary quantities.
            yield IteratorResult(worst=idx,
                                 ustar=ustar,
                                 vstar=vstar,
                                 loglstar=loglstar,
                                 logvol=logvol,
                                 logwt=logwt,
                                 logz=logz,
                                 logzvar=logzvar,
                                 h=h,
                                 nc=1,
                                 worst_it=point_it,
                                 boundidx=boundidx,
                                 bounditer=bounditer,
                                 eff=self.eff,
                                 delta_logz=delta_logz)

    def _remove_live_points(self):
        """Remove the final set of live points if they were
        previously added to the current set of dead points."""

        if self.added_live:
            self.added_live = False
            if self.save_samples:
                for k in [
                        'id', 'u', 'v', 'logl', 'logvol', 'logwt', 'logz',
                        'logzvar', 'h', 'nc', 'boundidx', 'it', 'bounditer',
                        'scale'
                ]:
                    del self.saved_run.D[k][-self.nlive:]
        else:
            raise ValueError("No live points were added to the "
                             "list of samples!")

    def sample(self,
               maxiter=None,
               maxcall=None,
               dlogz=0.01,
               logl_max=np.inf,
               n_effective=np.inf,
               add_live=True,
               save_bounds=True,
               save_samples=True):
        """
        **The main nested sampling loop.** Iteratively replace the worst live
        point with a sample drawn uniformly from the prior until the
        provided stopping criteria are reached. Instantiates a generator
        that will be called by the user.

        Parameters
        ----------
        maxiter : int, optional
            Maximum number of iterations. Iteration may stop earlier if the
            termination condition is reached. Default is `sys.maxsize`
            (no limit).

        maxcall : int, optional
            Maximum number of likelihood evaluations. Iteration may stop
            earlier if termination condition is reached. Default is
            `sys.maxsize` (no limit).

        dlogz : float, optional
            Iteration will stop when the estimated contribution of the
            remaining prior volume to the total evidence falls below
            this threshold. Explicitly, the stopping criterion is
            `ln(z + z_est) - ln(z) < dlogz`, where `z` is the current
            evidence from all saved samples and `z_est` is the estimated
            contribution from the remaining volume. Default is `0.01`.

        logl_max : float, optional
            Iteration will stop when the sampled ln(likelihood) exceeds the
            threshold set by `logl_max`. Default is no bound (`np.inf`).

        n_effective: int, optional
            Minimum number of effective posterior samples. If the estimated
            "effective sample size" (ESS) exceeds this number,
            sampling will terminate. Default is no ESS (`np.inf`).

        add_live : bool, optional
            Whether or not to add the remaining set of live points to
            the list of samples when calculating `n_effective`.
            Default is `True`.

        save_bounds : bool, optional
            Whether or not to save past distributions used to bound
            the live points internally. Default is `True`.

        save_samples : bool, optional
            Whether or not to save past samples from the nested sampling run
            (along with other ancillary quantities) internally.
            Default is `True`.

        Returns
        -------
        worst : int
            Index of the live point with the worst likelihood. This is our
            new dead point sample.

        ustar : `~numpy.ndarray` with shape (npdim,)
            Position of the sample.

        vstar : `~numpy.ndarray` with shape (ndim,)
            Transformed position of the sample.

        loglstar : float
            Ln(likelihood) of the sample.

        logvol : float
            Ln(prior volume) within the sample.

        logwt : float
            Ln(weight) of the sample.

        logz : float
            Cumulative ln(evidence) up to the sample (inclusive).

        logzvar : float
            Estimated cumulative variance on `logz` (inclusive).

        h : float
            Cumulative information up to the sample (inclusive).

        nc : int
            Number of likelihood calls performed before the new
            live point was accepted.

        worst_it : int
            Iteration when the live (now dead) point was originally proposed.

        boundidx : int
            Index of the bound the dead point was originally drawn from.

        bounditer : int
            Index of the bound being used at the current iteration.

        eff : float
            The cumulative sampling efficiency (in percent).

        delta_logz : float
            The estimated remaining evidence expressed as the ln(ratio) of the
            current evidence.

        """

        # Initialize quantities.
        if maxcall is None:
            maxcall = sys.maxsize
        if maxiter is None:
            maxiter = sys.maxsize
        self.save_samples = save_samples
        self.save_bounds = save_bounds
        ncall = 0

        # Check whether we're starting fresh or continuing a previous run.
        if self.it == 1 or len(self.saved_run.D['logl']) == 0:
            # Initialize values for nested sampling loop.
            h = 0.  # information, initially *0.*
            logz = -1.e300  # ln(evidence), initially *0.*
            logzvar = 0.  # var[ln(evidence)], initially *0.*
            logvol = 0.  # initially contains the whole prior (volume=1.)
            loglstar = -1.e300  # initial ln(likelihood)
            delta_logz = 1.e300  # ln(ratio) of total/current evidence

            # Check if we should initialize a different bounding distribution
            # instead of using the unit cube.
            if self._beyond_unit_bound(loglstar):
                bound = self.update()
                if self.save_bounds:
                    self.bound.append(bound)
                    self.nbound += 1
                self.since_update = 0
        else:
            # Remove live points (if added) from previous run.
            if self.added_live:
                self._remove_live_points()

            # Get final state from previous run.
            h = self.saved_run.D['h'][-1]  # information
            logz = self.saved_run.D['logz'][-1]  # ln(evidence)
            logzvar = self.saved_run.D['logzvar'][-1]  # var[ln(evidence)]
            logvol = self.saved_run.D['logvol'][-1]  # ln(volume)
            loglstar = min(self.live_logl)  # ln(likelihood)
            delta_logz = np.logaddexp(
                logz,
                np.max(self.live_logl) + logvol) - logz  # log-evidence ratio

        stop_iterations = False
        # The main nested sampling loop.
        for it in range(sys.maxsize):

            # Stopping criterion 1: current number of iterations
            # exceeds `maxiter`.
            if it > maxiter:
                stop_iterations = True

            # Stopping criterion 2: current number of `loglikelihood`
            # calls exceeds `maxcall`.
            if ncall > maxcall:
                stop_iterations = True

            # Stopping criterion 3: estimated (fractional) remaining evidence
            # lies below some threshold set by `dlogz`.
            logz_remain = np.max(self.live_logl) + logvol
            delta_logz = np.logaddexp(logz, logz_remain) - logz
            if dlogz is not None and delta_logz < dlogz:
                stop_iterations = True

            # Stopping criterion 4: last dead point exceeded the upper
            # `logl_max` bound.
            if loglstar > logl_max:
                stop_iterations = True

            # Stopping criterion 5: the number of effective posterior
            # samples has been achieved.
            if n_effective is not None:
                if self.n_effective > n_effective:
                    if add_live:
                        self.add_final_live(print_progress=False)
                        neff = self.n_effective
                        self._remove_live_points()
                        self.added_live = False
                    else:
                        neff = self.n_effective
                    if neff > n_effective:
                        stop_iterations = True

            if stop_iterations:
                if not self.save_samples:
                    # If dumping past states, save only the required quantities
                    # TODO I don't quite understand why we do this
                    add_info = dict(logz=logz,
                                    logzvar=logzvar,
                                    h=h,
                                    logvol=logvol,
                                    logl=loglstar)
                    self.saved_run.append(add_info)
                break

            # Expected ln(volume) shrinkage.
            logvol -= self.dlv

            # After `update_interval` interations have passed *and* we meet
            # the criteria for moving beyond sampling from the unit cube,
            # update the bound using the current set of live points.
            ucheck = self.since_update >= self.update_interval
            bcheck = self._beyond_unit_bound(loglstar)
            if ucheck and bcheck:
                bound = self.update()
                if self.save_bounds:
                    self.bound.append(bound)
                self.nbound += 1
                self.since_update = 0

            # Locate the "live" point with the lowest `logl`.
            worst = np.argmin(self.live_logl)  # index
            worst_it = self.live_it[worst]  # when point was proposed
            boundidx = self.live_bound[worst]  # associated bound index

            # Set our new worst likelihood constraint.
            # Notice we are doing copies here because live_u and live_v
            # are updated in-place
            ustar = self.live_u[worst].copy()  # unit cube position
            vstar = self.live_v[worst].copy()  # transformed position
            loglstar_new = self.live_logl[worst]  # new likelihood

            # Sample a new live point from within the likelihood constraint
            # `logl > loglstar` using the bounding distribution and sampling
            # method from our sampler.
            u, v, logl, nc = self._new_point(loglstar_new)
            ncall += nc
            self.ncall += nc
            self.since_update += nc

            (logwt, logz, logzvar,
             h) = progress_integration(loglstar, loglstar_new, logz, logzvar,
                                       logvol, self.dlv, h)
            loglstar = loglstar_new

            # Compute bound index at the current iteration.
            if self._beyond_unit_bound(loglstar):
                bounditer = self.nbound - 1
            else:
                bounditer = 0

            # Save the worst live point. It is now a "dead" point.
            if self.save_samples:
                self.saved_run.append(
                    dict(id=worst,
                         u=ustar,
                         v=vstar,
                         logl=loglstar,
                         logvol=logvol,
                         logwt=logwt,
                         logz=logz,
                         logzvar=logzvar,
                         h=h,
                         nc=nc,
                         it=worst_it,
                         bounditer=bounditer,
                         scale=self.scale))

            # Update the live point (previously our "worst" point).
            self.live_u[worst] = u
            self.live_v[worst] = v
            self.live_logl[worst] = logl
            self.live_bound[worst] = bounditer
            self.live_it[worst] = self.it

            # Compute our sampling efficiency.
            self.eff = 100. * self.it / self.ncall

            # Increment total number of iterations.
            self.it += 1

            # Return dead point and ancillary quantities.
            yield IteratorResult(worst=worst,
                                 ustar=ustar,
                                 vstar=vstar,
                                 loglstar=loglstar,
                                 logvol=logvol,
                                 logwt=logwt,
                                 logz=logz,
                                 logzvar=logzvar,
                                 h=h,
                                 nc=nc,
                                 worst_it=worst_it,
                                 boundidx=boundidx,
                                 bounditer=bounditer,
                                 eff=self.eff,
                                 delta_logz=delta_logz)

    def run_nested(self,
                   maxiter=None,
                   maxcall=None,
                   dlogz=None,
                   logl_max=np.inf,
                   n_effective=None,
                   add_live=True,
                   print_progress=True,
                   print_func=None,
                   save_bounds=True):
        """
        **A wrapper that executes the main nested sampling loop.**
        Iteratively replace the worst live point with a sample drawn
        uniformly from the prior until the provided stopping criteria
        are reached.

        Parameters
        ----------
        maxiter : int, optional
            Maximum number of iterations. Iteration may stop earlier if the
            termination condition is reached. Default is `sys.maxsize`
            (no limit).

        maxcall : int, optional
            Maximum number of likelihood evaluations. Iteration may stop
            earlier if termination condition is reached. Default is
            `sys.maxsize` (no limit).

        dlogz : float, optional
            Iteration will stop when the estimated contribution of the
            remaining prior volume to the total evidence falls below
            this threshold. Explicitly, the stopping criterion is
            `ln(z + z_est) - ln(z) < dlogz`, where `z` is the current
            evidence from all saved samples and `z_est` is the estimated
            contribution from the remaining volume. If `add_live` is `True`,
            the default is `1e-3 * (nlive - 1) + 0.01`. Otherwise, the
            default is `0.01`.

        logl_max : float, optional
            Iteration will stop when the sampled ln(likelihood) exceeds the
            threshold set by `logl_max`. Default is no bound (`np.inf`).

        n_effective: int, optional
            Minimum number of effective posterior samples. If the estimated
            "effective sample size" (ESS) exceeds this number,
            sampling will terminate. Default is no ESS (`np.inf`).

        add_live : bool, optional
            Whether or not to add the remaining set of live points to
            the list of samples at the end of each run. Default is `True`.

        print_progress : bool, optional
            Whether or not to output a simple summary of the current run that
            updates with each iteration. Default is `True`.

        print_func : function, optional
            A function that prints out the current state of the sampler.
            If not provided, the default :meth:`results.print_fn` is used.

        save_bounds : bool, optional
            Whether or not to save past bounding distributions used to bound
            the live points internally. Default is *True*.

        """

        # Define our stopping criteria.
        if dlogz is None:
            if add_live:
                dlogz = 1e-3 * (self.nlive - 1.) + 0.01
            else:
                dlogz = 0.01

        # Run the main nested sampling loop.
        pbar, print_func = get_print_func(print_func, print_progress)
        try:
            ncall = self.ncall
            for it, results in enumerate(
                    self.sample(maxiter=maxiter,
                                maxcall=maxcall,
                                dlogz=dlogz,
                                logl_max=logl_max,
                                save_bounds=save_bounds,
                                save_samples=True,
                                n_effective=n_effective,
                                add_live=add_live)):
                ncall += results.nc

                # Print progress.
                if print_progress:
                    i = self.it - 1
                    print_func(results,
                               i,
                               ncall,
                               dlogz=dlogz,
                               logl_max=logl_max)

            # Add remaining live points to samples.
            if add_live:
                it = self.it - 1
                for i, results in enumerate(self.add_live_points()):
                    ncall += results.nc

                    # Print progress.
                    if print_progress:
                        print_func(results,
                                   it,
                                   ncall,
                                   add_live_it=i + 1,
                                   dlogz=dlogz,
                                   logl_max=logl_max)
        finally:
            if pbar is not None:
                pbar.close()
            self.loglikelihood.history_save()

    def add_final_live(self, print_progress=True, print_func=None):
        """
        **A wrapper that executes the loop adding the final live points.**
        Adds the final set of live points to the pre-existing sequence of
        dead points from the current nested sampling run.

        Parameters
        ----------
        print_progress : bool, optional
            Whether or not to output a simple summary of the current run that
            updates with each iteration. Default is `True`.

        print_func : function, optional
            A function that prints out the current state of the sampler.
            If not provided, the default :meth:`results.print_fn` is used.

        """

        if print_func is None:
            print_func = print_fn

        # Add remaining live points to samples.
        pbar, print_func = get_print_func(print_func, print_progress)
        try:
            ncall = self.ncall
            it = self.it - 1
            for i, results in enumerate(self.add_live_points()):

                # Print progress.
                if print_progress:
                    print_func(results,
                               it,
                               ncall,
                               add_live_it=i + 1,
                               dlogz=0.01)
        finally:
            if pbar is not None:
                pbar.close() class UnitCube:
    """
    An N-dimensional unit cube.

    Parameters
    ----------
    ndim : int
        The number of dimensions of the unit cube.

    """
    def __init__(self, ndim):
        self.n = ndim  # dimension
        self.vol = 1.  # volume
        self.funit = 1.  # overlap with the unit cube

    def contains(self, x):
        """Checks if unit cube contains the point `x`."""

        return unitcheck(x)

    def randoffset(self, rstate=None):
        """Draw a random offset from the center of the unit cube."""

        return self.sample(rstate=rstate) - 0.5

    def sample(self, rstate=None):
        """
        Draw a sample uniformly distributed within the unit cube.

        Returns
        -------
        x : `~numpy.ndarray` with shape (ndim,)
            A coordinate within the unit cube.

        """

        return rstate.uniform(size=self.n)

    def samples(self, nsamples, rstate=None):
        """
        Draw `nsamples` samples randomly distributed within the unit cube.

        Returns
        -------
        x : `~numpy.ndarray` with shape (nsamples, ndim)
            A collection of coordinates within the unit cube.

        """

        xs = np.array([self.sample(rstate=rstate) for i in range(nsamples)])

        return xs

    def update(self, points, rstate=None, bootstrap=0, pool=None):
        """Filler function."""

        pass class Ellipsoid:
    """
    An N-dimensional ellipsoid defined by::

        (x - v)^T A (x - v) = 1

    where the vector `v` is the center of the ellipsoid and `A` is a
    symmetric, positive-definite `N x N` matrix.

    Parameters
    ----------
    ctr : `~numpy.ndarray` with shape (N,)
        Coordinates of ellipsoid center.

    cov : `~numpy.ndarray` with shape (N, N)
        Covariance matrix describing the axes.

    """
    def __init__(self, ctr, cov, am=None, axes=None):
        self.n = len(ctr)  # dimension
        self.ctr = np.asarray(ctr)  # center coordinates
        self.cov = np.asarray(cov)  # covariance matrix
        if axes is None:
            self.axes = lalg.cholesky(cov, lower=True)  # transformation axes
        else:
            self.axes = axes

        # The eigenvalues (l) of `a` are (a^-2, b^-2, ...) where
        # (a, b, ...) are the lengths of principle axes.
        # The eigenvectors (v) are the normalized principle axes.
        l, v = lalg.eigh(self.cov)
        if np.all((l > 0.) & (np.isfinite(l))):
            self.axlens = np.sqrt(l)
            # Volume of ellipsoid is the volume of an n-sphere
            # is a product of squares of eigen values
            self.logvol = logvol_prefactor(self.n) + 0.5 * np.log(l).sum()
        else:
            raise ValueError("The input precision matrix defining the "
                             "ellipsoid {0} is apparently singular with "
                             "l={1} and v={2}.".format(self.cov, l, v))
        if am is None:
            self.am = v @ np.diag(1. / l) @ v.T
            # precision matrix (inverse of covariance)
        else:
            self.am = am

        # Scaled eigenvectors are the principle axes, where `paxes[:,i]` is the
        # i-th axis. Multiplying this matrix by a vector will transform a
        # point in the unit n-sphere to a point in the ellipsoid.
        self.paxes = np.dot(v, np.diag(self.axlens))

        # Amount by which volume was increased after initialization (i.e.
        # cumulative factor from `scale_to_vol`).
        self.expand = 1.
        self.funit = 1

    def scale_to_logvol(self, logvol):
        """Scale ellipsoid to a target volume."""

        logf = (logvol - self.logvol)
        # log of the maxium axis length of the ellipsoid
        max_log_axlen = np.log(np.sqrt(self.n) / 2)
        log_axlen = np.log(self.axlens)
        if log_axlen.max() < max_log_axlen - logf / self.n:
            # we are safe to inflate the ellipsoid isothropically
            # without hitting boundaries
            f = np.exp(logf / self.n)
            self.cov *= f**2
            self.am *= 1. / f**2
            self.axlens *= f
            self.axes *= f
        else:
            logfax = np.zeros(self.n)
            curlogf = logf  # how much we have left to inflate
            curn = self.n  # how many dimensions left
            l, v = lalg.eigh(self.cov)

            # here we start from largest and go to smallest
            for curi in np.argsort(l)[::-1]:
                delta = max(
                    min(max_log_axlen - log_axlen[curi], curlogf / curn), 0)
                logfax[curi] = delta
                curlogf -= delta
                curn -= 1
            fax = np.exp(logfax)  # linear inflation of each dimension
            l1 = l * fax**2  # eigen values are squares of axes
            self.cov = v @ np.diag(l1) @ v.T
            self.am = v @ np.diag(1 / l1) @ v.T
            self.axlens *= fax
            self.axes = lalg.cholesky(self.cov, lower=True)
            # I don't quite know how to scale axes, so I rerun cholesky
        self.logvol = logvol

    def major_axis_endpoints(self):
        """Return the endpoints of the major axis."""

        i = np.argmax(self.axlens)  # find the major axis
        v = self.paxes[:, i]  # vector from center to major axis endpoint

        return self.ctr - v, self.ctr + v

    def distance(self, x):
        """Compute the normalized distance to `x` from the center of the
        ellipsoid."""

        d = x - self.ctr

        return np.sqrt(np.dot(np.dot(d, self.am), d))

    def distance_many(self, x):
        """Compute the normalized distance to `x` from the center of the
        ellipsoid."""

        d = x - self.ctr[None, :]

        return np.sqrt(np.einsum('ij,jk,ik->i', d, self.am, d))

    def contains(self, x):
        """Checks if ellipsoid contains `x`."""

        return self.distance(x) <= 1.0

    def randoffset(self, rstate=None):
        """Return a random offset from the center of the ellipsoid."""

        return np.dot(self.axes, randsphere(self.n, rstate=rstate))

    def sample(self, rstate=None):
        """
        Draw a sample uniformly distributed within the ellipsoid.

        Returns
        -------
        x : `~numpy.ndarray` with shape (ndim,)
            A coordinate within the ellipsoid.

        """

        return self.ctr + self.randoffset(rstate=rstate)

    def samples(self, nsamples, rstate=None):
        """
        Draw `nsamples` samples uniformly distributed within the ellipsoid.

        Returns
        -------
        x : `~numpy.ndarray` with shape (nsamples, ndim)
            A collection of coordinates within the ellipsoid.

        """

        xs = np.array([self.sample(rstate=rstate) for i in range(nsamples)])

        return xs

    def unitcube_overlap(self, ndraws=10000, rstate=None):
        """Using `ndraws` Monte Carlo draws, estimate the fraction of
        overlap between the ellipsoid and the unit cube."""

        samples = [self.sample(rstate=rstate) for i in range(ndraws)]
        nin = sum([unitcheck(x) for x in samples])

        return 1. * nin / ndraws

    def update(self,
               points,
               rstate=None,
               bootstrap=0,
               pool=None,
               mc_integrate=False):
        """
        Update the ellipsoid to bound the collection of points.

        Parameters
        ----------
        points : `~numpy.ndarray` with shape (npoints, ndim)
            The set of points to bound.

        rstate : `~numpy.random.Generator`, optional
            `~numpy.random.Generator` instance.

        bootstrap : int, optional
            The number of bootstrapped realizations of the ellipsoid. The
            maximum distance to the set of points "left out" during each
            iteration is used to enlarge the resulting volumes.
            Default is `0`.

        pool : user-provided pool, optional
            Use this pool of workers to execute operations in parallel.

        mc_integrate : bool, optional
            Whether to use Monte Carlo methods to compute the effective
            overlap of the final ellipsoid with the unit cube.
            Default is `False`.

        """

        # Compute new bounding ellipsoid.
        ell = bounding_ellipsoid(points)
        self.n = ell.n
        self.ctr = ell.ctr
        self.cov = ell.cov
        self.am = ell.am
        self.logvol = ell.logvol
        self.axlens = ell.axlens
        self.axes = ell.axes
        self.paxes = ell.paxes
        self.expand = ell.expand

        # Use bootstrapping to determine the volume expansion factor.
        if bootstrap > 0:

            # If provided, compute bootstraps in parallel using a pool.
            if pool is None:
                M = map
            else:
                M = pool.map
            multis = [False for it in range(bootstrap)]
            ps = [points for it in range(bootstrap)]
            seeds = get_seed_sequence(rstate, bootstrap)
            args = zip(multis, ps, seeds)
            expands = list(M(_ellipsoid_bootstrap_expand, args))

            # Conservatively set the expansion factor to be the maximum
            # factor derived from our set of bootstraps.
            expand = max(expands)

            # If our ellipsoid is over-constrained, expand it.
            if expand > 1.:
                lv = self.logvol + self.n * np.log(expand)
                self.scale_to_logvol(lv)

        # Estimate the fractional overlap with the unit cube using
        # Monte Carlo integration.
        if mc_integrate:
            self.funit = self.unitcube_overlap(rstate=rstate) class MultiEllipsoid:
    """
    A collection of M N-dimensional ellipsoids.

    Parameters
    ----------
    ells : list of `Ellipsoid` objects with length M, optional
        A set of `Ellipsoid` objects that make up the collection of
        N-ellipsoids. Used to initialize :class:`MultiEllipsoid` if provided.

    ctrs : `~numpy.ndarray` with shape (M, N), optional
        Collection of coordinates of ellipsoid centers. Used to initialize
        :class:`MultiEllipsoid` if :data:`ams` is also provided.

    covs : `~numpy.ndarray` with shape (M, N, N), optional
        Collection of matrices describing the axes of the ellipsoids. Used to
        initialize :class:`MultiEllipsoid` if :data:`ctrs` also provided.

    """
    def __init__(self, ells=None, ctrs=None, covs=None):
        if ells is not None:
            # Try to initialize quantities using provided `Ellipsoid` objects.
            if (ctrs is None) and (covs is None):
                self.nells = len(ells)
                self.ells = ells
            else:
                raise ValueError("You cannot specific both `ells` and "
                                 "(`ctrs`, `covs`)!")
        else:
            # Try to initialize quantities using provided `ctrs` and `covs`.
            if (ctrs is None) and (covs is None):
                raise ValueError("You must specify either `ells` or "
                                 "(`ctrs`, `covs`).")
            else:
                self.nells = len(ctrs)
                self.ells = [
                    Ellipsoid(ctrs[i], covs[i]) for i in range(self.nells)
                ]
        self.__update_arrays()

        # Compute quantities.
        self.expands = np.ones(self.nells)
        self.logvol_tot = logsumexp(self.logvols)
        self.expand_tot = 1.
        self.funit = 1

    def __update_arrays(self):
        """
        Update internal arrays to ensure that in sync with ells
        """
        self.ctrs = np.array([ell.ctr for ell in self.ells])
        self.covs = np.array([ell.cov for ell in self.ells])
        self.ams = np.array([ell.am for ell in self.ells])
        self.logvols = np.array([ell.logvol for ell in self.ells])

    def scale_to_logvol(self, logvols):
        """Scale ellipoids to a corresponding set of
        target volumes.
        """
        for i in range(self.nells):
            self.ells[i].scale_to_logvol(logvols[i])

        # IMPORTANT We must also update arrays ams, covs
        self.__update_arrays()

        self.expands = np.array(
            [self.ells[i].expand for i in range(self.nells)])
        logvol_tot = logsumexp(logvols)
        self.expand_tot *= np.exp(logvol_tot - self.logvol_tot)
        self.logvol_tot = logvol_tot

    def major_axis_endpoints(self):
        """Return the endpoints of the major axis of each ellipsoid."""

        return np.array([ell.major_axis_endpoints() for ell in self.ells])

    def within(self, x, j=None):
        """Checks which ellipsoid(s) `x` falls within, skipping the `j`-th
        ellipsoid if need be."""

        delt = x[None, :] - self.ctrs
        mask = np.einsum('ai,aij,aj->a', delt, self.ams, delt) < 1
        if j is not None:
            mask[j] = False
        return np.nonzero(mask)[0]

    def overlap(self, x, j=None):
        """Checks how many ellipsoid(s) `x` falls within, skipping the `j`-th
        ellipsoid."""

        q = len(self.within(x, j=j))

        return q

    def contains(self, x):
        """Checks if the set of ellipsoids contains `x`."""
        delt = x[None, :] - self.ctrs
        return np.any(np.einsum('ai,aij,aj->a', delt, self.ams, delt) < 1)

    def sample(self, rstate=None, return_q=False):
        """
        Sample a point uniformly distributed within the *union* of ellipsoids.

        Returns
        -------
        x : `~numpy.ndarray` with shape (ndim,)
            A coordinate within the set of ellipsoids.

        idx : int
            The index of the ellipsoid `x` was sampled from.

        q : int, optional
            The number of ellipsoids `x` falls within.

        """

        # If there is only one ellipsoid, sample from it.
        if self.nells == 1:
            x = self.ells[0].sample(rstate=rstate)
            idx = 0
            q = 1
            if return_q:
                return x, idx, q
            else:
                return x, idx

        probs = np.exp(self.logvols - self.logvol_tot)
        while True:
            # Select an ellipsoid at random proportional to its volume.
            idx = rand_choice(probs, rstate)

            # Select a point from the chosen ellipsoid.
            x = self.ells[idx].sample(rstate=rstate)

            # Check how many ellipsoids the point lies within
            delts = (x[None, :] - self.ctrs)
            q = (np.einsum('ai,aij,aj->a', delts, self.ams, delts) < 1).sum()

            assert q > 0  # Should never fail

            if return_q:
                # If `q` is being returned, assume the user wants to
                # explicitly apply the `1. / q` acceptance criterion to
                # properly sample from the union of ellipsoids.
                return x, idx, q
            else:
                # If `q` is not being returned, assume the user wants this
                # done internally so we repeat the loop if needed
                if q == 1 or rstate.uniform() < (1. / q):
                    return x, idx

    def samples(self, nsamples, rstate=None):
        """
        Draw `nsamples` samples uniformly distributed within the *union* of
        ellipsoids.

        Returns
        -------
        xs : `~numpy.ndarray` with shape (nsamples, ndim)
            A collection of coordinates within the set of ellipsoids.

        """

        xs = np.array([self.sample(rstate=rstate)[0] for i in range(nsamples)])

        return xs

    def monte_carlo_logvol(self,
                           ndraws=10000,
                           rstate=None,
                           return_overlap=True):
        """Using `ndraws` Monte Carlo draws, estimate the log volume of the
        *union* of ellipsoids. If `return_overlap=True`, also returns the
        estimated fractional overlap with the unit cube."""

        # Estimate volume using Monte Carlo integration.
        samples = [
            self.sample(rstate=rstate, return_q=True) for i in range(ndraws)
        ]
        qsum = sum([q for (x, idx, q) in samples])
        logvol = np.log(ndraws * 1. / qsum) + self.logvol_tot

        if return_overlap:
            # Estimate the fractional amount of overlap with the
            # unit cube using the same set of samples.
            qin = sum([q * unitcheck(x) for (x, idx, q) in samples])
            overlap = 1. * qin / qsum
            return logvol, overlap
        else:
            return logvol

    def update(self,
               points,
               rstate=None,
               bootstrap=0,
               pool=None,
               mc_integrate=False):
        """
        Update the set of ellipsoids to bound the collection of points.

        Parameters
        ----------
        points : `~numpy.ndarray` with shape (npoints, ndim)
            The set of points to bound.

        rstate : `~numpy.random.Generator`, optional
            `~numpy.random.Generator` instance.

        bootstrap : int, optional
            The number of bootstrapped realizations of the ellipsoids. The
            maximum distance to the set of points "left out" during each
            iteration is used to enlarge the resulting volumes.
            Default is `0`.

        pool : user-provided pool, optional
            Use this pool of workers to execute operations in parallel.

        mc_integrate : bool, optional
            Whether to use Monte Carlo methods to compute the effective
            volume and fractional overlap of the final union of ellipsoids
            with the unit cube. Default is `False`.

        """

        npoints, ndim = points.shape
        if npoints == 1:
            raise RuntimeError('Cannot compute the bounding ellipsoid of '
                               'a single point.')

        # Calculate the bounding ellipsoid for the points, possibly
        # enlarged to a minimum volume.
        firstell = bounding_ellipsoid(points)

        # Recursively split the bounding ellipsoid
        ells = _bounding_ellipsoids(points, firstell)

        # Update the set of ellipsoids.
        self.nells = len(ells)
        self.ells = ells
        self.__update_arrays()
        # Sanity check: all points must be contained in some ellipsoid
        if not all(self.contains(p) for p in points):
            # refuse to update
            raise RuntimeError('Rejecting invalid MultiEllipsoid region')
        self.logvol_tot = logsumexp(self.logvols)

        # Compute expansion factor.
        expands = np.array([ell.expand for ell in self.ells])
        logvols_orig = self.logvols - np.log(expands)
        logvol_tot_orig = logsumexp(logvols_orig)
        self.expand_tot = np.exp(self.logvol_tot - logvol_tot_orig)

        # Use bootstrapping to determine the volume expansion factor.
        if bootstrap > 0:

            # If provided, compute bootstraps in parallel using a pool.
            if pool is None:
                M = map
            else:
                M = pool.map
            multis = [True for it in range(bootstrap)]
            ps = [points for it in range(bootstrap)]
            seeds = get_seed_sequence(rstate, bootstrap)
            args = zip(multis, ps, seeds)
            expands = list(M(_ellipsoid_bootstrap_expand, args))

            # Conservatively set the expansion factor to be the maximum
            # factor derived from our set of bootstraps.
            expand = max(expands)

            # If our ellipsoids are overly constrained, expand them.
            if expand > 1.:
                lvs = self.logvols + ndim * np.log(expand)
                self.scale_to_logvol(lvs)

        # Estimate the volume and fractional overlap with the unit cube
        # using Monte Carlo integration.
        if mc_integrate:
            self.logvol_tot, self.funit = self.monte_carlo_logvol(
                rstate=rstate, return_overlap=True) class RadFriends:
    """
    A collection of N-balls of identical size centered on each live point.

    Parameters
    ----------
    ndim : int
        The number of dimensions of each ball.

    cov : `~numpy.ndarray` with shape `(ndim, ndim)`, optional
        Covariance structure (correlation and size) of each ball.

    """
    def __init__(self, ndim, cov=None):
        self.n = ndim

        if cov is None:
            cov = np.identity(self.n)
        self.cov = cov
        self.am = lalg.pinvh(self.cov)
        self.axes = lalg.sqrtm(self.cov)
        self.axes_inv = lalg.pinvh(self.axes)

        detsign, detln = linalg.slogdet(self.am)
        assert detsign > 0
        self.logvol_ball = (logvol_prefactor(self.n) - 0.5 * detln)
        self.expand = 1.
        self.funit = 1

    def scale_to_logvol(self, logvol):
        """Scale ball to encompass a target volume."""

        f = np.exp((logvol - self.logvol_ball) * (1.0 / self.n))
        # linear factor
        self.cov *= f**2
        self.am /= f**2
        self.axes *= f
        self.axes_inv /= f
        self.logvol_ball = logvol

    def within(self, x, ctrs):
        """Check which balls `x` falls within."""

        # Execute a brute-force search over all balls.
        idxs = np.where(
            lalg.norm(np.dot(ctrs - x, self.axes_inv), axis=1) <= 1.)[0]

        return idxs

    def overlap(self, x, ctrs):
        """Check how many balls `x` falls within."""

        q = len(self.within(x, ctrs))

        return q

    def contains(self, x, ctrs):
        """Check if the set of balls contains `x`."""

        return self.overlap(x, ctrs) > 0

    def sample(self, ctrs, rstate=None, return_q=False):
        """
        Sample a point uniformly distributed within the *union* of balls.

        Returns
        -------
        x : `~numpy.ndarray` with shape (ndim,)
            A coordinate within the set of balls.

        q : int, optional
            The number of balls `x` falls within.

        """

        nctrs = len(ctrs)  # number of balls

        # If there is only one ball, sample from it.
        if nctrs == 1:
            ds = randsphere(self.n, rstate=rstate)
            dx = np.dot(ds, self.axes)
            x = ctrs[0] + dx
            if return_q:
                return x, 1
            else:
                return x

        # Select a ball at random.
        idx = rstate.integers(nctrs)

        # Select a point from the chosen ball.
        ds = randsphere(self.n, rstate=rstate)
        dx = np.dot(ds, self.axes)
        x = ctrs[idx] + dx

        # Check how many balls the point lies within, passing over
        # the `idx`-th ball `x` was sampled from.
        q = self.overlap(x, ctrs)

        if return_q:
            # If `q` is being returned, assume the user wants to
            # explicitly apply the `1. / q` acceptance criterion to
            # properly sample from the union of balls.
            return x, q
        else:
            # If `q` is not being returned, assume the user wants this
            # done internally.
            while rstate.uniform() > (1. / q):
                idx = rstate.integers(nctrs)
                ds = randsphere(self.n, rstate=rstate)
                dx = np.dot(ds, self.axes)
                x = ctrs[idx] + dx
                q = self.overlap(x, ctrs)
            return x

    def samples(self, nsamples, ctrs, rstate=None):
        """
        Draw `nsamples` samples uniformly distributed within the *union* of
        balls.

        Returns
        -------
        xs : `~numpy.ndarray` with shape (nsamples, ndim)
            A collection of coordinates within the set of balls.

        """

        xs = np.array(
            [self.sample(ctrs, rstate=rstate) for i in range(nsamples)])

        return xs

    def monte_carlo_logvol(self,
                           ctrs,
                           ndraws=10000,
                           rstate=None,
                           return_overlap=True):
        """Using `ndraws` Monte Carlo draws, estimate the log volume of the
        *union* of balls. If `return_overlap=True`, also returns the
        estimated fractional overlap with the unit cube."""

        # Estimate volume using Monte Carlo integration.
        samples = [
            self.sample(ctrs, rstate=rstate, return_q=True)
            for i in range(ndraws)
        ]
        qsum = sum([q for (x, q) in samples])
        logvol = np.log(1. * ndraws / qsum * len(ctrs)) + self.logvol_ball

        if return_overlap:
            # Estimate the fractional amount of overlap with the
            # unit cube using the same set of samples.
            qin = sum([q * unitcheck(x) for (x, q) in samples])
            overlap = 1. * qin / qsum
            return logvol, overlap
        else:
            return logvol

    def update(self,
               points,
               rstate=None,
               bootstrap=0,
               pool=None,
               mc_integrate=False,
               use_clustering=True):
        """
        Update the radii of our balls.

        Parameters
        ----------
        points : `~numpy.ndarray` with shape (npoints, ndim)
            The set of points to bound.

        rstate : `~numpy.random.Generator`, optional
            `~numpy.random.Generator` instance.

        bootstrap : int, optional
            The number of bootstrapped realizations of the ellipsoids. The
            maximum distance to the set of points "left out" during each
            iteration is used to enlarge the resulting volumes.
            Default is `0`.

        pool : user-provided pool, optional
            Use this pool of workers to execute operations in parallel.

        mc_integrate : bool, optional
            Whether to use Monte Carlo methods to compute the effective
            volume and fractional overlap of the final union of balls
            with the unit cube. Default is `False`.

        use_clustering : bool, optional
            Whether to use clustering to avoid issues with widely-seperated
            modes. Default is `True`.

        """

        # If possible, compute bootstraps in parallel using a pool.
        if pool is None:
            M = map
        else:
            M = pool.map

        # Get new covariance.
        if use_clustering:
            self.cov = self._get_covariance_from_clusters(points)
        else:
            self.cov = self._get_covariance_from_all_points(points)
        self.am = lalg.pinvh(self.cov)
        self.axes = lalg.sqrtm(self.cov)
        self.axes_inv = lalg.pinvh(self.axes)

        # Decorrelate and re-scale points.
        points_t = np.dot(points, self.axes_inv)

        if bootstrap == 0.:
            # Construct radius using leave-one-out if no bootstraps used.
            radii = _friends_leaveoneout_radius(points_t, 'balls')
        else:
            # Bootstrap radius using the set of live points.
            ps = [points_t for it in range(bootstrap)]
            ftypes = ['balls' for it in range(bootstrap)]
            seeds = get_seed_sequence(rstate, bootstrap)
            args = zip(ps, ftypes, seeds)
            radii = list(M(_friends_bootstrap_radius, args))

        # Conservatively set radius to be maximum of the set.
        rmax = max(radii)

        # Re-scale axes.
        self.cov *= rmax**2
        self.am /= rmax**2
        self.axes *= rmax
        self.axes_inv /= rmax

        # Compute volume.
        detsign, detln = linalg.slogdet(self.am)
        assert detsign > 0
        self.logvol_ball = (logvol_prefactor(self.n) - 0.5 * detln)
        self.expand = 1.

        # Estimate the volume and fractional overlap with the unit cube
        # using Monte Carlo integration.
        if mc_integrate:
            self.funit = self.monte_carlo_logvol(points,
                                                 return_overlap=True,
                                                 rstate=rstate)[1]

    def _get_covariance_from_all_points(self, points):
        """Compute covariance using all points."""

        return np.cov(points, rowvar=False)

    def _get_covariance_from_clusters(self, points):
        """Compute covariance from re-centered clusters."""

        # Compute pairwise distances.
        distances = spatial.distance.pdist(points,
                                           metric='mahalanobis',
                                           VI=self.am)

        # Identify conglomerates of points by constructing a linkage matrix.
        linkages = cluster.hierarchy.single(distances)

        # Cut when linkage between clusters exceed the radius.
        clusteridxs = cluster.hierarchy.fcluster(linkages,
                                                 1.0,
                                                 criterion='distance')
        nclusters = np.max(clusteridxs)
        if nclusters == 1:
            return self._get_covariance_from_all_points(points)
        else:
            i = 0
            overlapped_points = np.empty_like(points)
            for idx in np.unique(clusteridxs):
                group_points = points[clusteridxs == idx, :]
                group_mean = group_points.mean(axis=0).reshape((1, -1))
                j = i + len(group_points)
                overlapped_points[i:j, :] = group_points - group_mean
                i = j
            return self._get_covariance_from_all_points(overlapped_points) class SupFriends:
    """
    A collection of N-cubes of identical size centered on each live point.

    Parameters
    ----------
    ndim : int
        The number of dimensions of the cube.

    cov : `~numpy.ndarray` with shape `(ndim, ndim)`, optional
        Covariance structure (correlation and size) of each cube.

    """
    def __init__(self, ndim, cov=None):
        self.n = ndim

        if cov is None:
            cov = np.identity(self.n)
        self.cov = cov
        self.am = lalg.pinvh(self.cov)
        self.axes = lalg.sqrtm(self.cov)
        self.axes_inv = lalg.pinvh(self.axes)

        detsign, detln = linalg.slogdet(self.am)
        assert detsign > 0
        self.logvol_cube = self.n * np.log(2.) - 0.5 * detln
        self.expand = 1.
        self.funit = 1

    def scale_to_logvol(self, logvol):
        """Scale cube to encompass a target volume."""

        f = np.exp((logvol - self.logvol_cube) * (1.0 / self.n))
        # linear factor
        self.cov *= f**2
        self.am /= f**2
        self.axes *= f
        self.axes_inv /= f
        self.logvol_cube = logvol

    def within(self, x, ctrs):
        """Checks which cubes `x` falls within."""

        # Execute a brute-force search over all cubes.
        idxs = np.where(
            np.max(np.abs(np.dot(ctrs - x, self.axes_inv)), axis=1) <= 1.)[0]

        return idxs

    def overlap(self, x, ctrs):
        """Checks how many cubes `x` falls within, skipping the `j`-th
        cube."""

        q = len(self.within(x, ctrs))

        return q

    def contains(self, x, ctrs):
        """Checks if the set of cubes contains `x`."""

        return self.overlap(x, ctrs) > 0

    def sample(self, ctrs, rstate=None, return_q=False):
        """
        Sample a point uniformly distributed within the *union* of cubes.

        Returns
        -------
        x : `~numpy.ndarray` with shape (ndim,)
            A coordinate within the set of cubes.

        q : int, optional
            The number of cubes `x` falls within.

        """

        nctrs = len(ctrs)  # number of cubes

        # If there is only one cube, sample from it.
        if nctrs == 1:
            ds = (2. * rstate.uniform(size=self.n) - 1.)
            dx = np.dot(ds, self.axes)
            x = ctrs[0] + dx
            if return_q:
                return x, 1
            else:
                return x

        # Select a cube at random.
        idx = rstate.integers(nctrs)

        # Select a point from the chosen cube.
        ds = (2. * rstate.uniform(size=self.n) - 1.)
        dx = np.dot(ds, self.axes)
        x = ctrs[idx] + dx

        # Check how many cubes the point lies within, passing over
        # the `idx`-th cube `x` was sampled from.
        q = self.overlap(x, ctrs)

        if return_q:
            # If `q` is being returned, assume the user wants to
            # explicitly apply the `1. / q` acceptance criterion to
            # properly sample from the union of balls.
            return x, q
        else:
            # If `q` is not being returned, assume the user wants this
            # done internally.
            while rstate.uniform() > (1. / q):
                idx = rstate.integers(nctrs)
                ds = (2. * rstate.uniform(size=self.n) - 1.)
                dx = np.dot(ds, self.axes)
                x = ctrs[idx] + dx
                q = self.overlap(x, ctrs)
            return x

    def samples(self, nsamples, ctrs, rstate=None):
        """
        Draw `nsamples` samples uniformly distributed within the *union* of
        cubes.

        Returns
        -------
        xs : `~numpy.ndarray` with shape (nsamples, ndim)
            A collection of coordinates within the set of cubes.

        """

        xs = np.array(
            [self.sample(ctrs, rstate=rstate) for i in range(nsamples)])

        return xs

    def monte_carlo_logvol(self,
                           ctrs,
                           ndraws=10000,
                           rstate=None,
                           return_overlap=False):
        """Using `ndraws` Monte Carlo draws, estimate the log volume of the
        *union* of cubes. If `return_overlap=True`, also returns the
        estimated fractional overlap with the unit cube."""

        # Estimate the volume using Monte Carlo integration.
        samples = [
            self.sample(ctrs, rstate=rstate, return_q=True)
            for i in range(ndraws)
        ]
        qsum = sum([q for (x, q) in samples])
        logvol = np.log(1. * ndraws / qsum * len(ctrs)) + self.logvol_cube

        if return_overlap:
            # Estimate the fractional overlap with the unit cube using
            # the same set of samples.
            qin = sum([q * unitcheck(x) for (x, q) in samples])
            overlap = 1. * qin / qsum
            return logvol, overlap
        else:
            return logvol

    def update(self,
               points,
               rstate=None,
               bootstrap=0,
               pool=None,
               mc_integrate=False,
               use_clustering=True):
        """
        Update the half-side-lengths of our cubes.

        Parameters
        ----------
        points : `~numpy.ndarray` with shape (npoints, ndim)
            The set of points to bound.

        rstate : `~numpy.random.Generator`, optional
            `~numpy.random.Generator` instance.

        bootstrap : int, optional
            The number of bootstrapped realizations of the ellipsoids. The
            maximum distance to the set of points "left out" during each
            iteration is used to enlarge the resulting volumes.
            Default is `0`.

        pool : user-provided pool, optional
            Use this pool of workers to execute operations in parallel.

        mc_integrate : bool, optional
            Whether to use Monte Carlo methods to compute the effective
            volume and fractional overlap of the final union of cubes
            with the unit cube. Default is `False`.

        use_clustering : bool, optional
            Whether to use clustering to avoid issues with widely-seperated
            modes. Default is `True`.

        """

        # If possible, compute bootstraps in parallel using a pool.
        if pool is None:
            M = map
        else:
            M = pool.map

        # Get new covariance.
        if use_clustering:
            self.cov = self._get_covariance_from_clusters(points)
        else:
            self.cov = self._get_covariance_from_all_points(points)
        self.am = lalg.pinvh(self.cov)
        self.axes = lalg.sqrtm(self.cov)
        self.axes_inv = lalg.pinvh(self.axes)

        # Decorrelate and re-scale points.
        points_t = np.dot(points, self.axes_inv)

        if bootstrap == 0.:
            # Construct radius using leave-one-out if no bootstraps used.
            hsides = _friends_leaveoneout_radius(points_t, 'cubes')
        else:
            # Bootstrap radius using the set of live points.
            ps = [points_t for it in range(bootstrap)]
            ftypes = ['cubes' for it in range(bootstrap)]
            seeds = get_seed_sequence(rstate, bootstrap)
            args = zip(ps, ftypes, seeds)
            hsides = list(M(_friends_bootstrap_radius, args))

        # Conservatively set half-side-length to be maximum of the set.
        hsmax = max(hsides)

        # Re-scale axes.
        self.cov *= hsmax**2
        self.am /= hsmax**2
        self.axes *= hsmax
        self.axes_inv /= hsmax

        detsign, detln = linalg.slogdet(self.am)
        self.logvol_cube = (self.n * np.log(2.) - 0.5 * detln)
        self.expand = 1.

        # Estimate the volume and fractional overlap with the unit cube
        # using Monte Carlo integration.
        if mc_integrate:
            self.funit = self.monte_carlo_logvol(points,
                                                 return_overlap=True,
                                                 rstate=rstate)[1]

    def _get_covariance_from_all_points(self, points):
        """Compute covariance using all points."""

        return np.cov(points, rowvar=False)

    def _get_covariance_from_clusters(self, points):
        """Compute covariance from re-centered clusters."""

        # Compute pairwise distances.
        distances = spatial.distance.pdist(points,
                                           metric='mahalanobis',
                                           VI=self.am)

        # Identify conglomerates of points by constructing a linkage matrix.
        linkages = cluster.hierarchy.single(distances)

        # Cut when linkage between clusters exceed the radius.
        clusteridxs = cluster.hierarchy.fcluster(linkages,
                                                 1.0,
                                                 criterion='distance')
        nclusters = np.max(clusteridxs)
        if nclusters == 1:
            return self._get_covariance_from_all_points(points)
        else:
            i = 0
            overlapped_points = np.empty_like(points)
            for idx in np.unique(clusteridxs):
                group_points = points[clusteridxs == idx, :]
                group_mean = group_points.mean(axis=0).reshape((1, -1))
                j = i + len(group_points)
                overlapped_points[i:j, :] = group_points - group_mean
                i = j
            return self._get_covariance_from_all_points(overlapped_points) def sample_unif(args):
    """
    Evaluate a new point sampled uniformly from a bounding proposal
    distribution. Parameters are zipped within `args` to utilize
    `pool.map`-style functions.

    Parameters
    ----------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the initial sample.

    loglstar : float
        Ln(likelihood) bound. **Not applicable here.**

    axes : `~numpy.ndarray` with shape (ndim, ndim)
        Axes used to propose new points. **Not applicable here.**

    scale : float
        Value used to scale the provided axes. **Not applicable here.**

    prior_transform : function
        Function transforming a sample from the a unit cube to the parameter
        space of interest according to the prior.

    loglikelihood : function
        Function returning ln(likelihood) given parameters as a 1-d `~numpy`
        array of length `ndim`.

    kwargs : dict
        A dictionary of additional method-specific parameters.
        **Not applicable here.**

    Returns
    -------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the final proposed point within the unit cube. **For
        uniform sampling this is the same as the initial input position.**

    v : `~numpy.ndarray` with shape (ndim,)
        Position of the final proposed point in the target parameter space.

    logl : float
        Ln(likelihood) of the final proposed point.

    nc : int
        Number of function calls used to generate the sample. For uniform
        sampling this is `1` by construction.

    blob : dict
        Collection of ancillary quantities used to tune :data:`scale`. **Not
        applicable for uniform sampling.**

    """

    # Unzipping.
    (u, loglstar, axes, scale, prior_transform, loglikelihood, rseed,
     kwargs) = args

    # Evaluate.
    v = prior_transform(np.asarray(u))
    logl = loglikelihood(np.asarray(v))
    nc = 1
    blob = None

    return u, v, logl, nc, blob def sample_rwalk(args):
    """
    Return a new live point proposed by random walking away from an
    existing live point.

    Parameters
    ----------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the initial sample. **This is a copy of an existing live
        point.**

    loglstar : float
        Ln(likelihood) bound.

    axes : `~numpy.ndarray` with shape (ndim, ndim)
        Axes used to propose new points. For random walks new positions are
        proposed using the :class:`~dynesty.bounding.Ellipsoid` whose
        shape is defined by axes.

    scale : float
        Value used to scale the provided axes.

    prior_transform : function
        Function transforming a sample from the a unit cube to the parameter
        space of interest according to the prior.

    loglikelihood : function
        Function returning ln(likelihood) given parameters as a 1-d `~numpy`
        array of length `ndim`.

    kwargs : dict
        A dictionary of additional method-specific parameters.

    Returns
    -------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the final proposed point within the unit cube.

    v : `~numpy.ndarray` with shape (ndim,)
        Position of the final proposed point in the target parameter space.

    logl : float
        Ln(likelihood) of the final proposed point.

    nc : int
        Number of function calls used to generate the sample.

    blob : dict
        Collection of ancillary quantities used to tune :data:`scale`.

    """

    # Unzipping.
    (u, loglstar, axes, scale, prior_transform, loglikelihood, rseed,
     kwargs) = args
    rstate = get_random_generator(rseed)
    return generic_random_walk(u, loglstar, axes, scale, prior_transform,
                               loglikelihood, rstate, kwargs) def sample_slice(args):
    """
    Return a new live point proposed by a series of random slices
    away from an existing live point. Standard "Gibs-like" implementation where
    a single multivariate "slice" is a combination of `ndim` univariate slices
    through each axis.

    Parameters
    ----------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the initial sample. **This is a copy of an existing live
        point.**

    loglstar : float
        Ln(likelihood) bound.

    axes : `~numpy.ndarray` with shape (ndim, ndim)
        Axes used to propose new points. For slices new positions are
        proposed along the arthogonal basis defined by :data:`axes`.

    scale : float
        Value used to scale the provided axes.

    prior_transform : function
        Function transforming a sample from the a unit cube to the parameter
        space of interest according to the prior.

    loglikelihood : function
        Function returning ln(likelihood) given parameters as a 1-d `~numpy`
        array of length `ndim`.

    kwargs : dict
        A dictionary of additional method-specific parameters.

    Returns
    -------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the final proposed point within the unit cube.

    v : `~numpy.ndarray` with shape (ndim,)
        Position of the final proposed point in the target parameter space.

    logl : float
        Ln(likelihood) of the final proposed point.

    nc : int
        Number of function calls used to generate the sample.

    blob : dict
        Collection of ancillary quantities used to tune :data:`scale`.

    """

    # Unzipping.
    (u, loglstar, axes, scale, prior_transform, loglikelihood, rseed,
     kwargs) = args
    rstate = get_random_generator(rseed)
    # Periodicity.
    nonperiodic = kwargs.get('nonperiodic', None)

    # Setup.
    n = len(u)
    assert (axes.shape[0] == n)
    slices = kwargs.get('slices', 5)  # number of slices
    nc = 0
    nexpand = 0
    ncontract = 0
    fscale = []

    # Modifying axes and computing lengths.
    axes = scale * axes.T  # scale based on past tuning

    # Slice sampling loop.
    for it in range(slices):

        # Shuffle axis update order.
        idxs = np.arange(n)
        rstate.shuffle(idxs)

        # Slice sample along a random direction.
        for idx in idxs:

            # Select axis.
            axis = axes[idx]
            (u_prop, v_prop, logl_prop, nc1, nexpand1, ncontract1,
             fscale1) = generic_slice_step(u, axis, nonperiodic, loglstar,
                                           loglikelihood, prior_transform,
                                           rstate)
            u = u_prop
            nc += nc1
            nexpand += nexpand1
            ncontract += ncontract1
            fscale.append(fscale1)

    blob = {
        'fscale': np.mean(fscale),
        'nexpand': nexpand,
        'ncontract': ncontract
    }

    return u_prop, v_prop, logl_prop, nc, blob def sample_rslice(args):
    """
    Return a new live point proposed by a series of random slices
    away from an existing live point. Standard "random" implementation where
    each slice is along a random direction based on the provided axes.

    Parameters
    ----------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the initial sample. **This is a copy of an existing live
        point.**

    loglstar : float
        Ln(likelihood) bound.

    axes : `~numpy.ndarray` with shape (ndim, ndim)
        Axes used to propose new slice directions.

    scale : float
        Value used to scale the provided axes.

    prior_transform : function
        Function transforming a sample from the a unit cube to the parameter
        space of interest according to the prior.

    loglikelihood : function
        Function returning ln(likelihood) given parameters as a 1-d `~numpy`
        array of length `ndim`.

    kwargs : dict
        A dictionary of additional method-specific parameters.

    Returns
    -------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the final proposed point within the unit cube.

    v : `~numpy.ndarray` with shape (ndim,)
        Position of the final proposed point in the target parameter space.

    logl : float
        Ln(likelihood) of the final proposed point.

    nc : int
        Number of function calls used to generate the sample.

    blob : dict
        Collection of ancillary quantities used to tune :data:`scale`.

    """

    # Unzipping.
    (u, loglstar, axes, scale, prior_transform, loglikelihood, rseed,
     kwargs) = args
    rstate = get_random_generator(rseed)
    # Periodicity.
    nonperiodic = kwargs.get('nonperiodic', None)

    # Setup.
    n = len(u)
    assert (axes.shape[0] == n)
    slices = kwargs.get('slices', 5)  # number of slices
    nc = 0
    nexpand = 0
    ncontract = 0
    fscale = []

    # Slice sampling loop.
    for it in range(slices):

        # Propose a direction on the unit n-sphere.
        drhat = rstate.standard_normal(size=n)
        drhat /= linalg.norm(drhat)

        # Transform and scale based on past tuning.
        direction = np.dot(axes, drhat) * scale

        (u_prop, v_prop, logl_prop, nc1, nexpand1, ncontract1,
         fscale1) = generic_slice_step(u, direction, nonperiodic, loglstar,
                                       loglikelihood, prior_transform, rstate)
        u = u_prop
        nc += nc1
        nexpand += nexpand1
        ncontract += ncontract1
        fscale.append(fscale1)

    blob = {
        'fscale': np.mean(fscale),
        'nexpand': nexpand,
        'ncontract': ncontract
    }

    return u_prop, v_prop, logl_prop, nc, blob def sample_hslice(args):
    """
    Return a new live point proposed by "Hamiltonian" Slice Sampling
    using a series of random trajectories away from an existing live point.
    Each trajectory is based on the provided axes and samples are determined
    by moving forwards/backwards in time until the trajectory hits an edge
    and approximately reflecting off the boundaries.
    Once a series of reflections has been established, we propose a new live
    point by slice sampling across the entire path.

    Parameters
    ----------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the initial sample. **This is a copy of an existing live
        point.**

    loglstar : float
        Ln(likelihood) bound.

    axes : `~numpy.ndarray` with shape (ndim, ndim)
        Axes used to propose new slice directions.

    scale : float
        Value used to scale the provided axes.

    prior_transform : function
        Function transforming a sample from the a unit cube to the parameter
        space of interest according to the prior.

    loglikelihood : function
        Function returning ln(likelihood) given parameters as a 1-d `~numpy`
        array of length `ndim`.

    kwargs : dict
        A dictionary of additional method-specific parameters.

    Returns
    -------
    u : `~numpy.ndarray` with shape (npdim,)
        Position of the final proposed point within the unit cube.

    v : `~numpy.ndarray` with shape (ndim,)
        Position of the final proposed point in the target parameter space.

    logl : float
        Ln(likelihood) of the final proposed point.

    nc : int
        Number of function calls used to generate the sample.

    blob : dict
        Collection of ancillary quantities used to tune :data:`scale`.

    """

    # Unzipping.
    (u, loglstar, axes, scale, prior_transform, loglikelihood, rseed,
     kwargs) = args
    rstate = get_random_generator(rseed)
    # Periodicity.
    nonperiodic = kwargs.get('nonperiodic', None)

    # Setup.
    n = len(u)
    assert (axes.shape[0] == len(u))
    slices = kwargs.get('slices', 5)  # number of slices
    grad = kwargs.get('grad', None)  # gradient of log-likelihood
    max_move = kwargs.get('max_move', 100)  # limit for `ncall`
    compute_jac = kwargs.get('compute_jac', False)  # whether Jacobian needed
    jitter = 0.25  # 25% jitter
    nc = 0
    nmove = 0
    nreflect = 0
    ncontract = 0

    # Slice sampling loop.
    for it in range(slices):
        # Define the left, "inner", and right "nodes" for a given chord.
        # We will plan to slice sampling using these chords.
        nodes_l, nodes_m, nodes_r = [], [], []

        # Propose a direction on the unit n-sphere.
        drhat = rstate.standard_normal(size=n)
        drhat /= linalg.norm(drhat)

        # Transform and scale based on past tuning.
        axis = np.dot(axes, drhat) * scale * 0.01

        # Create starting window.
        vel = np.array(axis)  # current velocity
        u_l = u.copy()
        u_r = u.copy()
        u_l -= rstate.uniform(1. - jitter, 1. + jitter) * vel
        u_r += rstate.uniform(1. - jitter, 1. + jitter) * vel
        nodes_l.append(np.array(u_l))
        nodes_m.append(np.array(u))
        nodes_r.append(np.array(u_r))

        # Progress "right" (i.e. "forwards" in time).
        reverse, reflect = False, False
        u_r = np.array(u)
        ncall = 0
        while ncall <= max_move:

            # Iterate until we can bracket the edge of the distribution.
            nodes_l.append(np.array(u_r))
            u_out, u_in = None, []
            while True:
                # Step forward.
                u_r += rstate.uniform(1. - jitter, 1. + jitter) * vel
                # Evaluate point.
                if unitcheck(u_r, nonperiodic):
                    v_r = prior_transform(np.asarray(u_r))
                    logl_r = loglikelihood(np.asarray(v_r))
                    nc += 1
                    ncall += 1
                    nmove += 1
                else:
                    logl_r = -np.inf
                # Check if we satisfy the log-likelihood constraint
                # (i.e. are "in" or "out" of bounds).
                if logl_r < loglstar:
                    if reflect:
                        # If we are out of bounds and just reflected, we
                        # reverse direction and terminate immediately.
                        reverse = True
                        nodes_l.pop()  # remove since chord does not exist
                        break
                    else:
                        # If we're already in bounds, then we're safe.
                        u_out = np.array(u_r)
                        logl_out = logl_r
                    # Check if we could compute gradients assuming we
                    # terminated with the current `u_out`.
                    if np.isfinite(logl_out):
                        reverse = False
                    else:
                        reverse = True
                else:
                    reflect = False
                    u_in.append(np.array(u_r))
                # Check if we've bracketed the edge.
                if u_out is not None:
                    break
            # Define the rest of our chord.
            if len(nodes_l) == len(nodes_r) + 1:
                if len(u_in) > 0:
                    u_in = u_in[rstate.choice(
                        len(u_in))]  # pick point randomly
                else:
                    u_in = np.array(u)
                    pass
                nodes_m.append(np.array(u_in))
                nodes_r.append(np.array(u_out))
            # Check if we have turned around.
            if reverse:
                break

            # Reflect off the boundary.
            u_r, logl_r = u_out, logl_out
            if grad is None:
                # If the gradient is not provided, we will attempt to
                # approximate it numerically using 2nd-order methods.
                h = np.zeros(n)
                for i in range(n):
                    u_r_l, u_r_r = np.array(u_r), np.array(u_r)
                    # right side
                    u_r_r[i] += 1e-10
                    if unitcheck(u_r_r, nonperiodic):
                        v_r_r = prior_transform(np.asarray(u_r_r))
                        logl_r_r = loglikelihood(np.asarray(v_r_r))
                    else:
                        logl_r_r = -np.inf
                        reverse = True  # can't compute gradient
                    nc += 1
                    # left side
                    u_r_l[i] -= 1e-10
                    if unitcheck(u_r_l, nonperiodic):
                        v_r_l = prior_transform(np.asarray(u_r_l))
                        logl_r_l = loglikelihood(np.asarray(v_r_l))
                    else:
                        logl_r_l = -np.inf
                        reverse = True  # can't compute gradient
                    if reverse:
                        break  # give up because we have to turn around
                    nc += 1
                    # compute dlnl/du
                    h[i] = (logl_r_r - logl_r_l) / 2e-10
            else:
                # If the gradient is provided, evaluate it.
                h = grad(v_r)
                if compute_jac:
                    jac = []
                    # Evaluate and apply Jacobian dv/du if gradient
                    # is defined as d(lnL)/dv instead of d(lnL)/du.
                    for i in range(n):
                        u_r_l, u_r_r = np.array(u_r), np.array(u_r)
                        # right side
                        u_r_r[i] += 1e-10
                        if unitcheck(u_r_r, nonperiodic):
                            v_r_r = prior_transform(np.asarray(u_r_r))
                        else:
                            reverse = True  # can't compute Jacobian
                            v_r_r = np.array(v_r)  # assume no movement
                        # left side
                        u_r_l[i] -= 1e-10
                        if unitcheck(u_r_l, nonperiodic):
                            v_r_l = prior_transform(np.asarray(u_r_l))
                        else:
                            reverse = True  # can't compute Jacobian
                            v_r_r = np.array(v_r)  # assume no movement
                        if reverse:
                            break  # give up because we have to turn around
                        jac.append((v_r_r - v_r_l) / 2e-10)
                    jac = np.array(jac)
                    h = np.dot(jac, h)  # apply Jacobian
                nc += 1
            # Compute specular reflection off boundary.
            vel_ref = vel - 2 * h * np.dot(vel, h) / linalg.norm(h)**2
            dotprod = np.dot(vel_ref, vel)
            dotprod /= linalg.norm(vel_ref) * linalg.norm(vel)
            # Check angle of reflection.
            if dotprod < -0.99:
                # The reflection angle is sufficiently small that it might
                # as well be a reflection.
                reverse = True
                break
            else:
                # If the reflection angle is sufficiently large, we
                # proceed as normal to the new position.
                vel = vel_ref
                u_out = None
                reflect = True
                nreflect += 1

        # Progress "left" (i.e. "backwards" in time).
        reverse, reflect = False, False
        vel = -np.array(axis)  # current velocity
        u_l = np.array(u)
        ncall = 0
        while ncall <= max_move:

            # Iterate until we can bracket the edge of the distribution.
            # Use a doubling approach to try and locate the bounds faster.
            nodes_r.append(np.array(u_l))
            u_out, u_in = None, []
            while True:
                # Step forward.
                u_l += rstate.uniform(1. - jitter, 1. + jitter) * vel
                # Evaluate point.
                if unitcheck(u_l, nonperiodic):
                    v_l = prior_transform(np.asarray(u_l))
                    logl_l = loglikelihood(np.asarray(v_l))
                    nc += 1
                    ncall += 1
                    nmove += 1
                else:
                    logl_l = -np.inf
                # Check if we satisfy the log-likelihood constraint
                # (i.e. are "in" or "out" of bounds).
                if logl_l < loglstar:
                    if reflect:
                        # If we are out of bounds and just reflected, we
                        # reverse direction and terminate immediately.
                        reverse = True
                        nodes_r.pop()  # remove since chord does not exist
                        break
                    else:
                        # If we're already in bounds, then we're safe.
                        u_out = np.array(u_l)
                        logl_out = logl_l
                    # Check if we could compute gradients assuming we
                    # terminated with the current `u_out`.
                    if np.isfinite(logl_out):
                        reverse = False
                    else:
                        reverse = True
                else:
                    reflect = False
                    u_in.append(np.array(u_l))
                # Check if we've bracketed the edge.
                if u_out is not None:
                    break
            # Define the rest of our chord.
            if len(nodes_r) == len(nodes_l) + 1:
                if len(u_in) > 0:
                    u_in = u_in[rstate.choice(
                        len(u_in))]  # pick point randomly
                else:
                    u_in = np.array(u)
                    pass
                nodes_m.append(np.array(u_in))
                nodes_l.append(np.array(u_out))
            # Check if we have turned around.
            if reverse:
                break

            # Reflect off the boundary.
            u_l, logl_l = u_out, logl_out
            if grad is None:
                # If the gradient is not provided, we will attempt to
                # approximate it numerically using 2nd-order methods.
                h = np.zeros(n)
                for i in range(n):
                    u_l_l, u_l_r = np.array(u_l), np.array(u_l)
                    # right side
                    u_l_r[i] += 1e-10
                    if unitcheck(u_l_r, nonperiodic):
                        v_l_r = prior_transform(np.asarray(u_l_r))
                        logl_l_r = loglikelihood(np.asarray(v_l_r))
                    else:
                        logl_l_r = -np.inf
                        reverse = True  # can't compute gradient
                    nc += 1
                    # left side
                    u_l_l[i] -= 1e-10
                    if unitcheck(u_l_l, nonperiodic):
                        v_l_l = prior_transform(np.asarray(u_l_l))
                        logl_l_l = loglikelihood(np.asarray(v_l_l))
                    else:
                        logl_l_l = -np.inf
                        reverse = True  # can't compute gradient
                    if reverse:
                        break  # give up because we have to turn around
                    nc += 1
                    # compute dlnl/du
                    h[i] = (logl_l_r - logl_l_l) / 2e-10
            else:
                # If the gradient is provided, evaluate it.
                h = grad(v_l)
                if compute_jac:
                    jac = []
                    # Evaluate and apply Jacobian dv/du if gradient
                    # is defined as d(lnL)/dv instead of d(lnL)/du.
                    for i in range(n):
                        u_l_l, u_l_r = np.array(u_l), np.array(u_l)
                        # right side
                        u_l_r[i] += 1e-10
                        if unitcheck(u_l_r, nonperiodic):
                            v_l_r = prior_transform(np.asarray(u_l_r))
                        else:
                            reverse = True  # can't compute Jacobian
                            v_l_r = np.array(v_l)  # assume no movement
                        # left side
                        u_l_l[i] -= 1e-10
                        if unitcheck(u_l_l, nonperiodic):
                            v_l_l = prior_transform(np.asarray(u_l_l))
                        else:
                            reverse = True  # can't compute Jacobian
                            v_l_r = np.array(v_l)  # assume no movement
                        if reverse:
                            break  # give up because we have to turn around
                        jac.append((v_l_r - v_l_l) / 2e-10)
                    jac = np.array(jac)
                    h = np.dot(jac, h)  # apply Jacobian
                nc += 1
            # Compute specular reflection off boundary.
            vel_ref = vel - 2 * h * np.dot(vel, h) / linalg.norm(h)**2
            dotprod = np.dot(vel_ref, vel)
            dotprod /= linalg.norm(vel_ref) * linalg.norm(vel)
            # Check angle of reflection.
            if dotprod < -0.99:
                # The reflection angle is sufficiently small that it might
                # as well be a reflection.
                reverse = True
                break
            else:
                # If the reflection angle is sufficiently large, we
                # proceed as normal to the new position.
                vel = vel_ref
                u_out = None
                reflect = True
                nreflect += 1

        # Initialize lengths of chords.
        if len(nodes_l) > 1:
            # remove initial fallback chord
            nodes_l.pop(0)
            nodes_m.pop(0)
            nodes_r.pop(0)
        nodes_l, nodes_m, nodes_r = (np.array(nodes_l), np.array(nodes_m),
                                     np.array(nodes_r))
        Nchords = len(nodes_l)
        axlen = np.zeros(Nchords, dtype='float')
        for i, (nl, nm, nr) in enumerate(zip(nodes_l, nodes_m, nodes_r)):
            axlen[i] = linalg.norm(nr - nl)

        # Slice sample from all chords simultaneously. This is equivalent to
        # slice sampling in *time* along our trajectory.
        axlen_init = np.array(axlen)
        while True:
            # Safety check.
            if np.any(axlen < 1e-5 * axlen_init):
                raise RuntimeError("Hamiltonian slice sampling appears to be "
                                   "stuck! Some useful output quantities:\n"
                                   "u: {0}\n"
                                   "u_left: {1}\n"
                                   "u_right: {2}\n"
                                   "loglstar: {3}.".format(
                                       u, u_l, u_r, loglstar))

            # Select chord.
            axprob = axlen / np.sum(axlen)
            idx = rstate.choice(Nchords, p=axprob)
            # Define chord.
            u_l, u_m, u_r = nodes_l[idx], nodes_m[idx], nodes_r[idx]
            u_hat = u_r - u_l
            rprop = rstate.uniform()
            u_prop = u_l + rprop * u_hat  # scale from left
            if unitcheck(u_prop, nonperiodic):
                v_prop = prior_transform(np.asarray(u_prop))
                logl_prop = loglikelihood(np.asarray(v_prop))
            else:
                logl_prop = -np.inf
            nc += 1
            ncontract += 1
            # If we succeed, move to the new position.
            if logl_prop > loglstar:
                u = u_prop
                break
            # If we fail, check if the new point is to the left/right of
            # the point interior to the bounds (`u_m`) and update
            # the bounds accordingly.
            else:
                s = np.dot(u_prop - u_m, u_hat)  # check sign (+/-)
                if s < 0:  # left
                    nodes_l[idx] = u_prop
                    axlen[idx] *= 1 - rprop
                elif s > 0:  # right
                    nodes_r[idx] = u_prop
                    axlen[idx] *= rprop
                else:
                    raise RuntimeError("Slice sampler has failed to find "
                                       "a valid point. Some useful "
                                       "output quantities:\n"
                                       "u: {0}\n"
                                       "u_left: {1}\n"
                                       "u_right: {2}\n"
                                       "u_hat: {3}\n"
                                       "u_prop: {4}\n"
                                       "loglstar: {5}\n"
                                       "logl_prop: {6}.".format(
                                           u, u_l, u_r, u_hat, u_prop,
                                           loglstar, logl_prop))

    blob = {'nmove': nmove, 'nreflect': nreflect, 'ncontract': ncontract}

    return u_prop, v_prop, logl_prop, nc, blob def unitcheck(u, nonbounded=None):
    """Check whether `u` is inside the unit cube. Given a masked array
    `nonbounded`, also allows periodic boundaries conditions to exceed
    the unit cube."""

    if nonbounded is None:
        # No periodic boundary conditions provided.
        return np.min(u) > 0 and np.max(u) < 1
    else:
        # Alternating periodic and non-periodic boundary conditions.
        unb = u[nonbounded]
        # pylint: disable=invalid-unary-operand-type
        ub = u[~nonbounded]
        return (unb.min() > 0 and unb.max() < 1 and ub.min() > -0.5
                and ub.max() < 1.5) def get_enlarge_bootstrap(sample, enlarge, bootstrap):
    """
    Determine the enlarge, bootstrap for a given run
    """
    # we should make it dimension dependent I think...
    DEFAULT_ENLARGE = 1.25
    DEFAULT_UNIF_BOOTSTRAP = 5
    if enlarge is not None and bootstrap is None:
        """If enlarge is specified and bootstrap is not we just use enlarge
        with no nootstrapping"""
        assert enlarge >= 1
        return enlarge, 0
    elif enlarge is None and bootstrap is not None:
        """
        If bootstrap is specified but enlarge is not we just use bootstrap
        And if we allow zero bootstrap if we want to force no bootstrap
        """
        assert ((bootstrap > 1) or (bootstrap == 0))
        return 1, bootstrap
    elif enlarge is None and bootstrap is None:
        """
        If neither enlarge or bootstrap are specified we are doing
        things in auto-mode. I.e. use enlarge unless the uniform
        sampler is selected
        """
        if sample == 'unif':
            return 1, DEFAULT_UNIF_BOOTSTRAP
        else:
            return DEFAULT_ENLARGE, 0
    else:
        """Both enlarge and bootstrap were specified"""
        if bootstrap == 0 or enlarge == 1:
            return enlarge, bootstrap
        else:
            raise ValueError('Enlarge and bootstrap together do not make'
                             'sense unless bootstrap=1 or enlarge = 1')
import math
import copy
import numpy as np
from .sampler import Sampler
from .bounding import (UnitCube, Ellipsoid, MultiEllipsoid, RadFriends,
                       SupFriends)
from .sampling import (sample_unif, sample_rwalk, sample_slice, sample_rslice,
                       sample_hslice)
from .utils import unitcheck, get_enlarge_bootstrap
        """Update the slice proposal scale based on the relative
        size of the slices compared to our initial guess.
        For slice sampling the scale is only 'advisory' in the sense that
        the right scale will just speed up sampling as we'll have to expand
        or contract less. It won't affect the quality of the samples much.
        """
        # see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4063214/
        # also 2002.06212
        # https://www.tandfonline.com/doi/full/10.1080/10618600.2013.791193
        # and https://github.com/joshspeagle/dynesty/issues/260
        nexpand, ncontract = max(blob['nexpand'], 1), blob['ncontract']
        mult = (nexpand * 2. / (nexpand + ncontract))
        # avoid drastic updates to the scale factor limiting to factor
        # of two
        mult = np.clip(mult, 0.5, 2)
        # Remember I can't apply the rule that scale < cube diagonal
        # because scale is multiplied by axes
        self.scale = self.scale * mult

    def update_hslice(self, blob):
        """Update the Hamiltonian slice proposal scale based
        on the relative amount of time spent moving vs reflecting."""

        nmove, nreflect = blob['nmove'], blob['nreflect']
        ncontract = blob.get('ncontract', 0)
        fmove = (1. * nmove) / (nmove + nreflect + ncontract + 2)
        norm = max(self.fmove, 1. - self.fmove)
        self.scale *= math.exp((fmove - self.fmove) / norm)

    def update_user(self, blob):
        """Update the scale based on the user-defined update function."""

        if callable(self.custom_update):
            self.scale = self.custom_update(blob, self.scale)


class UnitCubeSampler(SuperSampler):
    """
    Samples conditioned on the unit N-cube (i.e. with no bounds).

    Parameters
    ----------
    loglikelihood : function
        Function returning ln(likelihood) given parameters as a 1-d `~numpy`
        array of length `ndim`.

    prior_transform : function
        Function transforming a sample from the a unit cube to the parameter
        space of interest according to the prior.

    npdim : int
        Number of parameters accepted by `prior_transform`.

    live_points : list of 3 `~numpy.ndarray` each with shape (nlive, ndim)
        Initial set of "live" points. Contains `live_u`, the coordinates
        on the unit cube, `live_v`, the transformed variables, and
        `live_logl`, the associated loglikelihoods.

    method : {`'unif'`, `'rwalk'`,
              `'slice'`, `'rslice'`, `'hslice'`}, optional
        Method used to sample uniformly within the likelihood constraint,
        conditioned on the provided bounds.

    update_interval : int
        Only update the bounding distribution every `update_interval`-th
        likelihood call.

    first_update : dict
        A dictionary containing parameters governing when the sampler should
        first update the bounding distribution from the unit cube to the one
        specified by the user.

    rstate : `~numpy.random.Generator`
        `~numpy.random.Generator` instance.

    queue_size: int
        Carry out likelihood evaluations in parallel by queueing up new live
        point proposals using (at most) this many threads/members.

    pool: pool
        Use this pool of workers to execute operations in parallel.

    use_pool : dict, optional
        A dictionary containing flags indicating where the provided `pool`
        should be used to execute operations in parallel.

    kwargs : dict, optional
        A dictionary of additional parameters.

    """
    def __init__(self,
                 loglikelihood,
                 prior_transform,
                 npdim,
                 live_points,
                 method,
                 update_interval,
                 first_update,
                 rstate,
                 queue_size,
                 pool,
                 use_pool,
                 kwargs=None,
                 ncdim=0):

        # Initialize sampler.
        super().__init__(loglikelihood,
                         prior_transform,
                         npdim,
                         live_points,
                         method,
                         update_interval,
                         first_update,
                         rstate,
                         queue_size,
                         pool,
                         use_pool,
                         ncdim=ncdim,
                         kwargs=kwargs or {})


        self.unitcube = UnitCube(self.ncdim)