monitor.hpp

Go to the documentation of this file.

//============================================================================
// vSMC/include/vsmc/core/monitor.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_MONITOR_HPP
#define VSMC_CORE_MONITOR_HPP

#include <vsmc/internal/common.hpp>

#define VSMC_RUNTIME_ASSERT_CORE_MONITOR_ID(func)                             \
    VSMC_RUNTIME_ASSERT(                                                      \
        (id < dim()), "**Monitor::" #func "** INVALID ID NUMBER ARGUMENT")

#define VSMC_RUNTIME_ASSERT_CORE_MONITOR_ITER(func)                           \
    VSMC_RUNTIME_ASSERT((iter < iter_size()),                                 \
        "**Monitor::" #func "** INVALID ITERATION NUMBER ARGUMENT")

#define VSMC_RUNTIME_ASSERT_CORE_MONITOR_EVAL                                 \
    VSMC_RUNTIME_ASSERT(static_cast<bool>(eval_),                             \
        "**Monitor::eval** INVALID EVALUATION OBJECT")

namespace vsmc
{

enum MonitorStage {
    MonitorMove,     
    MonitorResample, 
    MonitorMCMC      
};                   // enum MonitorStage

template <typename T>
class Monitor
{
    public:
    using value_type = T;
    using eval_type =
        std::function<void(std::size_t, std::size_t, Particle<T> &, double *)>;

    Monitor(std::size_t dim, const eval_type &eval, bool record_only = false,
        MonitorStage stage = MonitorMCMC)
        : dim_(dim)
        , eval_(eval)
        , recording_(true)
        , record_only_(record_only)
        , stage_(stage)
        , name_(dim)
    {
    }

    std::size_t dim() const { return dim_; }

    bool record_only() const { return record_only_; }

    MonitorStage stage() const { return stage_; }

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

    void reserve(std::size_t num)
    {
        index_.reserve(num);
        record_.reserve(dim_ * num);
    }

    bool empty() const { return !static_cast<bool>(eval_); }

    std::string &name(std::size_t id)
    {
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ID(name);

        return name_[id];
    }

    const std::string &name(std::size_t id) const
    {
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ID(name);

        return name_[id];
    }

    std::size_t index(std::size_t iter) const
    {
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ITER(index);

        return index_[iter];
    }

    double record(std::size_t id) const
    {
        std::size_t iter = iter_size() ? iter_size() - 1 : iter_size();
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ID(record);
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ITER(record);

        return record_[iter * dim_ + id];
    }

    double record(std::size_t id, std::size_t iter) const
    {
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ID(record);
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ITER(record);

        return record_[iter * dim_ + id];
    }

    template <typename OutputIter>
    void read_index(OutputIter first) const
    {
        std::copy(index_.begin(), index_.end(), first);
    }

    const std::size_t *index_data() const { return index_.data(); }

    const double *record_data() const { return record_.data(); }

    const double *record_data(std::size_t iter) const
    {
        return record_.data() + iter * dim_;
    }

    template <typename OutputIter>
    void read_record(std::size_t id, OutputIter first) const
    {
        const std::size_t N = iter_size();
        const double *riter = record_.data() + id;
        for (std::size_t i = 0; i != N; ++i, ++first, riter += dim_)
            *first = *riter;
    }

    template <typename OutputIterIter>
    void read_record_matrix(OutputIterIter first) const
    {
        for (std::size_t d = 0; d != dim_; ++d, ++first)
            read_record(d, *first);
    }

    template <MatrixLayout Layout, typename OutputIter>
    void read_record_matrix(OutputIter first) const
    {
        const std::size_t N = iter_size();
        if (Layout == ColMajor) {
            for (std::size_t d = 0; d != dim_; ++d) {
                const double *riter = record_.data() + d;
                for (std::size_t i = 0; i != N; ++i, ++first, riter += dim_)
                    *first = *riter;
            }
        }

        if (Layout == RowMajor)
            std::copy(record_.begin(), record_.end(), first);
    }

    void set_eval(const eval_type &new_eval) { eval_ = new_eval; }

    void eval(std::size_t iter, Particle<T> &particle, MonitorStage stage)
    {
        if (!recording_)
            return;

        if (stage != stage_)
            return;

        VSMC_RUNTIME_ASSERT_CORE_MONITOR_EVAL;

        result_.resize(dim_);
        if (record_only_) {
            eval_(iter, dim_, particle, result_.data());
            push_back(iter);

            return;
        }

        const std::size_t N = static_cast<std::size_t>(particle.size());
        buffer_.resize(N * dim_);
        eval_(iter, dim_, particle, buffer_.data());
        ::cblas_dgemv(::CblasColMajor, ::CblasNoTrans,
            static_cast<VSMC_CBLAS_INT>(dim_), static_cast<VSMC_CBLAS_INT>(N),
            1.0, buffer_.data(), static_cast<VSMC_CBLAS_INT>(dim_),
            particle.weight().data(), 1, 0.0, result_.data(), 1);
        push_back(iter);
    }

    void clear()
    {
        index_.clear();
        record_.clear();
    }

    bool recording() const { return recording_; }

    void turn_on() { recording_ = true; }

    void turn_off() { recording_ = false; }

    private:
    std::size_t dim_;
    eval_type eval_;
    bool recording_;
    bool record_only_;
    MonitorStage stage_;
    Vector<std::string> name_;
    Vector<std::size_t> index_;
    Vector<double> record_;
    Vector<double> result_;
    Vector<double> buffer_;

    void push_back(std::size_t iter)
    {
        index_.push_back(iter);
        record_.insert(record_.end(), result_.begin(), result_.end());
    }
}; // class Monitor

} // namespace vsmc

#endif // VSMC_CORE_MONITOR_HPP