sample_vector.h

Go to the documentation of this file.

/* $Id: sample_vector.h 17393 2014-08-23 20:19:14Z predoehl $ */
/* =========================================================================== *
   |
   |  Copyright (c) 1994-2010 by Kobus Barnard (author)
   |
   |  Personal and educational use of this code is granted, provided that this
   |  header is kept intact, and that the authorship is not misrepresented, that
   |  its use is acknowledged in publications, and relevant papers are cited.
   |
   |  For other use contact the author (kobus AT cs DOT arizona DOT edu).
   |
   |  Please note that the code in this file has not necessarily been adequately
   |  tested. Naturally, there is no guarantee of performance, support, or fitness
   |  for any particular task. Nonetheless, I am interested in hearing about
   |  problems that you encounter.
   |
   |  Author:  Ernesto Brau
 * =========================================================================== */

#ifndef SAMPLE_VECTOR_H_INCLUDED
#define SAMPLE_VECTOR_H_INCLUDED

#include <algorithm>
#include "sample_cpp/sample_step.h"
#include "sample_cpp/sample_default.h"

#include <wrap_lapack/wrap_lapack.h>
#include <prob_cpp/prob_distribution.h>
#include <m_cpp/m_matrix.h>
#include <n_cpp/n_cholesky.h>

template<class Model>
inline
double get_vector_model_parameter(const Model& m, int i)
{
    return m[i];
}

template<class Model>
inline
void set_vector_model_parameter(Model& m, int i, double x)
{
    m[i] = x;
}

template<class Model>
inline
void move_vector_model_parameter(Model& m, int i, double x)
{
    Move_model_parameter_as_plus<Model> mmpap(set_vector_model_parameter<Model>,
                                              get_vector_model_parameter<Model>);
    mmpap(m, i, x);
}


template<class Model>
inline
int get_vector_model_dimension(const Model& m)
{
    return m.size();
}

template<typename Model>
class Vector_numerical_gradient : public Numerical_gradient<Model>
{
public:
    typedef Numerical_gradient<Model> Parent;

    typedef typename Parent::Target_distribution Target_distribution;
    typedef typename Parent::Get_neighborhood Get_neighborhood;

    Vector_numerical_gradient
    (
        const Target_distribution& log_target,
        const Get_neighborhood&    get_neighborhood
    ) :
        Parent
        (
            log_target,
            move_vector_model_parameter<Model>,
            get_neighborhood,
            get_vector_model_dimension<Model>
        )
    {}
};

template<typename Model, bool CONSTRAIN_TARGET = false, bool INCLUDE_ACCEPT_STEP = true, bool REVERSIBLE = true>
class Vector_hmc_step : public Basic_hmc_step<Model, CONSTRAIN_TARGET, INCLUDE_ACCEPT_STEP, REVERSIBLE>
{
public:
    typedef Basic_hmc_step<Model, CONSTRAIN_TARGET, INCLUDE_ACCEPT_STEP, REVERSIBLE> Parent;

    typedef typename Parent::Target_distribution Target_distribution;
    typedef typename Parent::Gradient Gradient;
    typedef typename Parent::Get_step_size Get_step_size;
    typedef typename Parent::Get_lower_bounds Get_lower_bounds;
    typedef typename Parent::Get_upper_bounds Get_upper_bounds;

    Vector_hmc_step
    (
        const Target_distribution& log_target, 
        int                        num_dynamics_steps,
        const Gradient&            gradient,
        const Get_step_size&       get_step_size,
        double                     alpha = 0.0
    ) :
        Parent(log_target,
               num_dynamics_steps,
               gradient,
               move_vector_model_parameter<Model>, 
               get_step_size,
               get_vector_model_dimension<Model>,
               alpha)
    {}

    Vector_hmc_step
    (
        const Target_distribution& log_target, 
        int                        num_dynamics_steps,
        const Gradient&            gradient,
        const Get_step_size&       get_step_size,
        const Get_lower_bounds&     get_lower_bounds,
        const Get_upper_bounds&     get_upper_bounds,
        double                     alpha = 0.0
    ) :
        Parent(log_target,
               num_dynamics_steps,
               gradient,
               move_vector_model_parameter<Model>, 
               get_step_size,
               get_vector_model_dimension<Model>,
               get_vector_model_parameter<Model>,
               get_lower_bounds,
               get_upper_bounds,
               alpha)
    {}

    virtual ~Vector_hmc_step()
    {}
};

template<class Model, bool REVERSIBLE = true>
struct Vector_sd_step
{
    typedef Vector_hmc_step<Model, false, REVERSIBLE> Type;
};

template<typename SimpleVector>
class Independent_gaussian_proposer
{
private:
    kjb::Vector scale_;

    Independent_gaussian_proposer(const kjb::Vector& covariance) : 
        scale_(covariance.size())
    {
        std::transform(covariance.begin(), covariance.end(), scale_.begin(), static_cast<double(*)(double)>(sqrt));
    }

    Mh_proposal_result operator()(const SimpleVector& m, SimpleVector& mp) const
    {
        if(m.size() != scale_.size())
        {
            KJB_THROW(kjb::Dimension_mismatch);
        }

        mp.resize(m.size());
        for(size_t i = 0; i < mp.size(); ++i)
        {
            mp[i] = m[i] + kjb::sample(kjb::STD_NORMAL) * scale_[i];
        }

        return Mh_proposal_result();
    }
};

template<typename VectorModel>
class Mv_gaussian_proposer
{
private:
    typedef Mv_gaussian_proposer Self;
    kjb::Matrix chol_cov_;
    double log_lambda_;
    double sqrt_t_;

public:
    Mv_gaussian_proposer(const kjb::Matrix& covariance) :
        chol_cov_(kjb::cholesky_decomposition(covariance)),
        log_lambda_(0.0),
        sqrt_t_(1.0)
    {}

    void set_temperature(double t)
    {
        sqrt_t_ = sqrt(t);
    }

    Mh_proposal_result operator()(const VectorModel& m, VectorModel& mp) const
    {
        size_t N = m.size();
        if(N != chol_cov_.get_num_rows() ||
           N != mp.size())
        {
            KJB_THROW(kjb::Dimension_mismatch);
        }

        kjb::Vector delta(N);
        std::generate(delta.begin(), delta.end(), boost::bind(Self::gauss_sample_));
        // scale by temperature, global scale, and covariance matrix
        delta = sqrt_t_ * exp(log_lambda_/2.0) * chol_cov_ * delta;

        std::transform(
                m.begin(),
                m.end(),
                delta.begin(),
                mp.begin(),
                std::plus<double>());

        return Mh_proposal_result();
    }

    const kjb::Matrix& get_unscaled_cholesky_covariance() const
    {
        return chol_cov_;
    }

    kjb::Matrix get_covariance() const
    {
        kjb::Matrix result = sqrt_t_ * exp(log_lambda_/2.0) * chol_cov_;
        result = result.transpose() * result;
        return result;
    }

    double get_global_scale() const
    {
        return sqrt(exp(log_lambda_));
    }

    template <class VectorType>
    void adapt_covariance(const VectorType& delta, double gamma)
    {
#ifdef TEST
//        kjb::Matrix ref = chol_cov_;
//        kjb::Matrix test = chol_cov_;
#endif 

        rank_1_update_2(chol_cov_, delta, gamma);
        
#ifdef TEST
        // try rank-one update a few different ways, to see if they're all equivalent
//
//        double diff;
//        rank_1_update_ref(ref, delta, gamma);
//        diff = max_abs_difference(chol_cov_,ref);
//        assert(diff < FLT_EPSILON);
//
//        assert(diff < FLT_EPSILON);
#endif /* TEST */
        
    }

    void adapt_global_scale(double x)
    {
        log_lambda_ += x;
    }

    void rank_1_update(kjb::Matrix& C, const VectorModel& delta, double gamma)
    {
#ifdef KJB_HAVE_LAPACK
        int N = delta.size();
        kjb::Matrix tmp(N, 1);
        std::copy(delta.begin(), delta.end(), &tmp(0,0));
        lapack_solve_triangular(
                C.get_c_matrix(),
                tmp.get_c_matrix(),
                &tmp.get_underlying_representation_with_guilt());

        tmp = tmp * tmp.transpose();

        // M -= I
        for(size_t i = 0; i < N; ++i)
            tmp(i,i) -= 1.0;

        lower_triangular_product_(C, tmp, tmp);
        tmp *= gamma;
        C += tmp;
#else
        KJB_THROW_2(kjb::Missing_dependency, "lapack");
#endif
    }

    void rank_1_update_2(kjb::Matrix& C, const VectorModel& delta, double gamma)
    {
        // TODO: DEBUG ME
        // I'm failing when gamma == 1.
        // Questions to answer:
        //    why am I multiplying C by sqrt(1-gamma)?
        //    check rasmussen's adaptive mh paper
        VectorModel x = sqrt(gamma) * delta;
        size_t N = delta.size();

        C *= sqrt(1-gamma);

        for(size_t k = 0; k < N; ++k)
        {
            double& C_k = C(k,k);
            const double x_k = x[k];

            double r = sqrt(C_k*C_k+x_k * x_k);
            double c = r/C_k;
            double s = x_k/C_k;
            C_k = r;
            for(size_t k2 = k+1; k2 < N; ++k2)
            {
                C(k2,k) = (C(k2,k) + s*x[k2])/c;
                x[k2] = c*x[k2] - s*C(k2, k);
            }
        }
    }


    void rank_1_update_ref(kjb::Matrix& C, const VectorModel& delta, double gamma)
    {
        kjb::Vector tmp(delta.begin(), delta.end());
        C = C * C.transpose();
        C += gamma * (kjb::outer_product(tmp, tmp)-C);
#ifdef TEST
        kjb::Matrix m_tmp = C;
#endif

        C = kjb::cholesky_decomposition(C);

#ifdef TEST
        assert(kjb::max_abs_difference(C * C.transpose(), m_tmp) < FLT_EPSILON);
#endif
    }

private:
    // utility
    static double gauss_sample_() { return kjb::sample(kjb::STD_NORMAL); }

    void lower_triangular_product_(
            const kjb::Matrix& op1,
            const kjb::Matrix& op2,
            kjb::Matrix& out)
    {
        const size_t N = op1.get_num_cols();
        kjb::Matrix tmp(N,N, 0.0);

        for(size_t r = 0; r < N; ++r)
        for(size_t c = 0; c <= r; ++c)
        for(size_t i = c; i <= r; ++i)
            tmp(r,c) += op1(r,i)*op2(i,c);
        out.swap(tmp);
    }
};

template<typename VectorModel>
class Vector_srw_step : public Basic_mh_step<VectorModel, Mv_gaussian_proposer<VectorModel> >
{
public:
    typedef Basic_mh_step<VectorModel, Mv_gaussian_proposer<VectorModel> > Base;
    typedef Vector_srw_step Self;

    typedef typename Base::Target_distribution Target_distribution;
    typedef typename Base::Proposer Proposer;

    Vector_srw_step
    (
        const Target_distribution& target,
        const kjb::Matrix& covariance
    ) :
        Base(target,
               Proposer(covariance))
    {}

    Vector_srw_step
    (
        const Target_distribution& target,
        const kjb::Vector& covariance
    ) :
        Base(target,
               Proposer(covariance))
    {}

    inline void set_covariance(const kjb::Matrix& covariance)
    {
        Proposer& proposer = Base::get_proposer_();
        proposer = Mv_gaussian_proposer<VectorModel>(covariance);
    }

    inline void set_covariance(const kjb::Vector& covariance)
    {
        Proposer& proposer = Base::get_proposer_();
        proposer = Mv_gaussian_proposer<VectorModel>(covariance);
    }

    virtual ~Vector_srw_step()
    {}
};

#endif /*SAMPLE_VECTOR_H_INCLUDED */