sampler.hpp

Go to the documentation of this file.

//============================================================================
// vSMC/include/vsmc/core/sampler.hpp
//----------------------------------------------------------------------------
//                         vSMC: Scalable Monte Carlo
//----------------------------------------------------------------------------
// Copyright (c) 2013-2016, Yan Zhou
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//   Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
//
//   Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//============================================================================

#ifndef VSMC_CORE_SAMPLER_HPP
#define VSMC_CORE_SAMPLER_HPP

#include <vsmc/internal/common.hpp>
#include <vsmc/core/monitor.hpp>
#include <vsmc/core/particle.hpp>

#define VSMC_RUNTIME_ASSERT_CORE_SAMPLER_MONITOR_NAME(iter, map, func)        \
    VSMC_RUNTIME_ASSERT(                                                      \
        (iter != map.end()), "**Sampler::" #func "** INVALID MONITOR NAME")

#define VSMC_RUNTIME_ASSERT_CORE_SAMPLER_FUNCTOR(func, caller, name)          \
    VSMC_RUNTIME_ASSERT(                                                      \
        (func), "**Sampler::" #caller "** INVALID " #name " OBJECT")

namespace vsmc
{

enum SamplerStage {
    SamplerInit = 1 << 0, 
    SamplerMove = 1 << 1, 
    SamplerMCMC = 1 << 2  
};                        // enum SamplerStage

inline constexpr SamplerStage operator&(SamplerStage s1, SamplerStage s2)
{
    return static_cast<SamplerStage>(
        static_cast<int>(s1) & static_cast<int>(s2));
}

inline constexpr SamplerStage operator|(SamplerStage s1, SamplerStage s2)
{
    return static_cast<SamplerStage>(
        static_cast<int>(s1) | static_cast<int>(s2));
}

inline constexpr SamplerStage operator^(SamplerStage s1, SamplerStage s2)
{
    return static_cast<SamplerStage>(
        static_cast<int>(s1) ^ static_cast<int>(s2));
}

inline constexpr SamplerStage operator~(SamplerStage s)
{
    return static_cast<SamplerStage>(~static_cast<int>(s));
}

inline SamplerStage &operator&=(SamplerStage &s1, SamplerStage s2)
{
    return s1 = s1 & s2;
}

inline SamplerStage &operator|=(SamplerStage &s1, SamplerStage s2)
{
    return s1 = s1 | s2;
}

inline SamplerStage &operator^=(SamplerStage &s1, SamplerStage s2)
{
    return s1 = s1 ^ s2;
}

template <typename T>
class Sampler
{
    public:
    using size_type = typename Particle<T>::size_type;
    using eval_type = std::function<std::size_t(std::size_t, Particle<T> &)>;

    template <typename... Args>
    explicit Sampler(Args &&... args)
        : particle_(std::forward<Args>(args)...)
        , iter_num_(0)
        , resample_threshold_(resample_threshold_never())
    {
    }

    Sampler<T> clone() const
    {
        Sampler<T> sampler(*this);
        sampler.particle().rng_set().seed();
        Seed::instance()(sampler.particle().rng());

        return sampler;
    }

    size_type size() const { return particle_.size(); }

    void reserve(std::size_t num)
    {
        num += iter_num_;
        size_history_.reserve(num);
        ess_history_.reserve(num);
        resampled_history_.reserve(num);
        status_history_.reserve(num);
        for (auto &m : monitor_)
            m.second.reserve(num);
    }

    std::size_t iter_size() const { return size_history_.size(); }

    std::size_t iter_num() const { return iter_num_; }

    Sampler<T> &eval(
        const eval_type &new_eval, SamplerStage stage, bool append = true)
    {
        if (!append)
            eval_.clear();
        eval_.push_back(std::make_pair(new_eval, stage));

        return *this;
    }

    Sampler<T> &resample_method(
        ResampleScheme scheme, double threshold = resample_threshold_always())
    {
        switch (scheme) {
            case Multinomial:
                resample_eval_ = ResampleEval<T>(ResampleMultinomial());
                break;
            case Residual:
                resample_eval_ = ResampleEval<T>(ResampleResidual());
                break;
            case Stratified:
                resample_eval_ = ResampleEval<T>(ResampleStratified());
                break;
            case Systematic:
                resample_eval_ = ResampleEval<T>(ResampleSystematic());
                break;
            case ResidualStratified:
                resample_eval_ = ResampleEval<T>(ResampleResidualStratified());
                break;
            case ResidualSystematic:
                resample_eval_ = ResampleEval<T>(ResampleResidualSystematic());
                break;
        }
        resample_threshold(threshold);

        return *this;
    }

    Sampler<T> &resample_method(const eval_type &res_eval,
        double threshold = resample_threshold_always())
    {
        resample_eval_ = res_eval;
        resample_threshold(threshold);

        return *this;
    }

    double resample_threshold() const { return resample_threshold_; }

    Sampler<T> &resample_threshold(double threshold)
    {
        resample_threshold_ = threshold;
        return *this;
    }

    static double resample_threshold_never()
    {
        return -std::numeric_limits<double>::infinity();
    }

    static double resample_threshold_always()
    {
        return std::numeric_limits<double>::infinity();
    }

    Sampler<T> &monitor(const std::string &name, const Monitor<T> &mon)
    {
        monitor_.insert(std::make_pair(name, mon));

        return *this;
    }

    Monitor<T> &monitor(const std::string &name)
    {
        typename monitor_map_type::iterator iter = monitor_.find(name);

        VSMC_RUNTIME_ASSERT_CORE_SAMPLER_MONITOR_NAME(iter, monitor_, monitor);

        return iter->second;
    }

    const Monitor<T> &monitor(const std::string &name) const
    {
        typename monitor_map_type::const_iterator citer = monitor_.find(name);

        VSMC_RUNTIME_ASSERT_CORE_SAMPLER_MONITOR_NAME(
            citer, monitor_, monitor);

        return citer->second;
    }

    bool monitor_clear(const std::string &name)
    {
        return monitor_.erase(name) ==
            static_cast<typename monitor_map_type::size_type>(1);
    }

    Sampler<T> &monitor_clear()
    {
        monitor_.clear();

        return *this;
    }

    Sampler<T> &initialize()
    {
        do_initialize();

        return *this;
    }

    Sampler<T> &iterate(std::size_t num = 1)
    {
        if (num > 1)
            reserve(num);
        for (std::size_t i = 0; i != num; ++i)
            do_iterate();

        return *this;
    }

    Particle<T> &particle() { return particle_; }

    const Particle<T> &particle() const { return particle_; }

    size_type size_history(std::size_t iter) const
    {
        return size_history_[iter];
    }

    template <typename OutputIter>
    OutputIter read_size_history(OutputIter first) const
    {
        return std::copy(size_history_.begin(), size_history_.end(), first);
    }

    double ess_history(std::size_t iter) const { return ess_history_[iter]; }

    template <typename OutputIter>
    OutputIter read_ess_history(OutputIter first) const
    {
        return std::copy(ess_history_.begin(), ess_history_.end(), first);
    }

    bool resampled_history(std::size_t iter) const
    {
        return resampled_history_[iter];
    }

    template <typename OutputIter>
    OutputIter read_resampled_history(OutputIter first) const
    {
        return std::copy(
            resampled_history_.begin(), resampled_history_.end(), first);
    }

    std::size_t status_history(std::size_t iter, std::size_t id) const
    {
        return status_history_[iter][id];
    }

    std::map<std::string, Vector<double>> summary() const
    {
        const double missing_data = std::numeric_limits<double>::quiet_NaN();

        std::map<std::string, Vector<double>> df;
        Vector<double> data(iter_size());

        std::copy(size_history_.begin(), size_history_.end(), data.begin());
        df["Size"] = data;

        std::copy(resampled_history_.begin(), resampled_history_.end(),
            data.begin());
        df["Resampled"] = data;

        df["ESS"] = ess_history_;

        std::size_t ssize = 0;
        for (std::size_t i = 0; i != iter_size(); ++i)
            if (ssize < status_history_[i].size())
                ssize = status_history_[i].size();
        for (std::size_t j = 0; j != ssize; ++j) {
            for (std::size_t i = 0; i != iter_size(); ++i) {
                data[i] = status_history_[i].size() > j ?
                    status_history_[i][j] :
                    missing_data;
            }
            df["Status." + std::to_string(j)] = data;
        }

        for (const auto &m : monitor_) {
            if (m.second.iter_size() > 0) {
                for (std::size_t d = 0; d != m.second.dim(); ++d) {
                    std::size_t miter = 0;
                    for (std::size_t i = 0; i != iter_size(); ++i) {
                        if (miter != m.second.iter_size() &&
                            i == m.second.index(miter)) {
                            data[i] = m.second.record(d, miter++);
                        } else {
                            data[i] = missing_data;
                        }
                    }
                    if (m.second.name(d).empty())
                        df[m.first + "." + std::to_string(d)] = data;
                    else
                        df[m.second.name(d)] = data;
                }
            }
        }

        return df;
    }

    template <typename CharT, typename Traits>
    std::basic_ostream<CharT, Traits> &print(
        std::basic_ostream<CharT, Traits> &os, char sepchar = '\t') const
    {
        if (!os || iter_size() == 0)
            return os;

        const std::map<std::string, Vector<double>> df = summary();

        for (const auto &v : df)
            os << v.first << sepchar;
        os << '\n';
        for (std::size_t i = 0; i != iter_size(); ++i) {
            for (const auto &v : df)
                os << v.second[i] << sepchar;
            os << '\n';
        }

        return os;
    }

    private:
    using monitor_map_type = std::map<std::string, Monitor<T>>;

    Particle<T> particle_;
    std::size_t iter_num_;

    Vector<std::pair<eval_type, SamplerStage>> eval_;
    eval_type resample_eval_;
    double resample_threshold_;

    Vector<size_type> size_history_;
    Vector<double> ess_history_;
    Vector<bool> resampled_history_;
    Vector<Vector<std::size_t>> status_history_;
    Vector<std::size_t> status_;
    monitor_map_type monitor_;

    void do_reset()
    {
        iter_num_ = 0;
        particle_.weight().set_equal();
        for (auto &m : monitor_)
            m.second.clear();

        size_history_.clear();
        ess_history_.clear();
        resampled_history_.clear();
        status_history_.clear();
        status_.clear();
    }

    void do_initialize()
    {
        do_reset();
        status_.clear();
        do_init();
        do_common();
    }

    void do_iterate()
    {
        ++iter_num_;
        status_.clear();
        do_move();
        do_common();
    }

    void do_common()
    {
        do_monitor(MonitorMove);
        do_resample(resample_threshold_);
        do_monitor(MonitorResample);
        do_mcmc();
        do_monitor(MonitorMCMC);
        status_history_.push_back(std::move(status_));
    }

    void do_init()
    {
        for (auto &e : eval_)
            if ((e.second & SamplerInit) != 0)
                status_.push_back(e.first(iter_num_, particle_));
    }

    void do_move()
    {
        for (auto &e : eval_)
            if ((e.second & SamplerMove) != 0)
                status_.push_back(e.first(iter_num_, particle_));
    }

    void do_mcmc()
    {
        for (auto &e : eval_)
            if ((e.second & SamplerMCMC) != 0)
                status_.push_back(e.first(iter_num_, particle_));
    }

    void do_resample(double threshold)
    {
        size_history_.push_back(size());
        ess_history_.push_back(particle_.weight().ess());

        if (ess_history_.back() < size() * threshold) {
            resampled_history_.push_back(true);
            resample_eval_(iter_num_, particle_);
        } else {
            resampled_history_.push_back(false);
        }
    }

    void do_monitor(MonitorStage stage)
    {
        for (auto &m : monitor_)
            if (!m.second.empty())
                m.second(iter_num_, particle_, stage);
    }
}; // class Sampler

template <typename CharT, typename Traits, typename T>
inline std::basic_ostream<CharT, Traits> &operator<<(
    std::basic_ostream<CharT, Traits> &os, const Sampler<T> &sampler)
{
    return sampler.print(os);
}

} // namespace vsmc

#endif // VSMC_CORE_SAMPLER_HPP