#ifndef MIX_METHOD_H_
#define MIX_METHOD_H_

#include <sampler/Metropolis.h>
#include <sampler/StepAdapter.h>
#include <vector>

class StepAdapter;

/**
 * @short Metropolis-Hastings sampling method with tempered transitions
 *
 * TemperedMetropolis implements a random-walk Metropolis-Hastings
 * method with tempered transitions (Neal, R. Sampling from multimodal
 * distributions using tempered transitions, Statistics and Computing,
 * vol 6, pp 353-355, 1996).
 *
 * A tempered log density at temperature T is derived from the target
 * log density L by replacing L with L/T.  For large T, the tempered
 * density is flattened, converging to the prior density as T tends to
 * infinity.
 *
 * In a tempered transition, a new Metropolis-Hasting proposal is
 * generated by a sequence of Metropolis-Hastings updates that are
 * stationary with respect to a tempered density.  In this sequence, T
 * starts close to 1 increases to a maximum (e.g. 1000) and then
 * decreases to 1 again. Tempering allows the proposal to jump to a
 * new local mode of the target density by tunnelling, at high
 * temperature, through a region of low density.  
 *
 * Three elements are essential for successful tempered updates. The
 * first is a sufficiently high maximum temperature. Without this,
 * tunnelling can never occur between well-separated modes of the
 * target density. The second element is a sufficiently small
 * temperature increment between levels.  If the increase is too
 * large, then the proposal generated by the tempered updates has a
 * low acceptance probability. The third element is a sufficient
 * number of replicates at each temperature level.  In
 * TemperedMetropolis the updates at each level are tuned to have an
 * acceptance probability converging to 0.234.  Several replicate
 * updates are therefore required to guarantee that a M-H jump takes
 * place. Without this, the proposal risks getting stuck in a region
 * of low probability. These elements are controlled by the parameters
 * max_level, max_temp, and nrep respectively.
 */
class TemperedMetropolis : public Metropolis
{
    const int _max_level;
    const unsigned int _nrep;
    const std::vector<double> _pwr;
    int _t;
    int _tmax;
    std::vector<StepAdapter*> _step_adapter;
    double _pmean;
    unsigned int _niter;
    void temperedUpdate(RNG *rng, double &log_prior0, double &log_likelihood0, 
                        std::vector<double> &value0);
public:
    /**
     * Constructor.
     *
     * @param nlevel Number of levels between minimum and maximum temperature
     *
     * @param max_temp  Maximum temperature
     *
     * @param nrep Number of Metropolis-Hastings updates to do at each
     * level
     */
    TemperedMetropolis(std::vector<double> const &value, 
                       int nlevel, double max_temp, unsigned int nrep);
    ~TemperedMetropolis();
    /**
     * Updates the current value using tempered transitions.
     */
    void update(RNG *rng);
    /**
     * Modifies the step size at each temperature level to achieve the
     * target acceptance probability using a noisy gradient algorithm
     *
     * When the mean acceptance probability is within the target
     * range, the temperature is increased up to max_temp.
     *
     * @param prob acceptance probability at current update
     */
    void rescale(double prob);
    /**
     * Checks whether the maximum temperature has been reached.
     */
    bool checkAdaptation() const;
    /**
     * The target density is assumed to be the product of a prior density
     * and a likelihood. Only the likelihood part of the density is 
     * tempered. This may be necessary in order to ensure that the 
     * tempered distribution can be normalized.
     *
     * This function returns the log of the prior density function. 
     */    
    virtual double logPrior() const = 0;
    /**
     * Returns the log of the likelihood.
     */
    virtual double logLikelihood() const = 0;
    /**
     * Modifies the given value vector in place by adding an
     * independent normal increment to each element.  It can be
     * overridden to provide non-normal increments, or a random walk
     * on some transformed scale (but see RMetropolis#logJacobian).
     *
     * Note that this function does not modify the value of the
     * RWMetropolis object.
     */
    virtual void step(std::vector<double> &value, double s, RNG *rng) const;
    /**
     * If the random walk takes place on a transformed scale
     * (e.g. log, logistic), then the log density of the target
     * distribution must be penalized by the log Jacobian of the
     * transformation.
     *
     * This function calculates the log Jacobian at the given value.
     */
    virtual double logJacobian(std::vector<double> const &value) const;
};

#endif /* MIX_METHOD_H_ */
