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(eval_, "**Monitor::eval** INVALID EVALUATION OBJECT")

namespace vsmc
{

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

template <typename T>
class Monitor
{
    public:
    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)
        , record_only_(record_only)
        , stage_(stage)
        , name_(dim)
    {
        internal::size_check<VSMC_BLAS_INT>(dim_, "Monitor::Monitor");
    }

    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 !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() const
    {
        std::size_t iter = iter_size() > 0 ? iter_size() - 1 : iter_size();
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ITER(index);

        return index_[iter];
    }

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

        return index_[iter];
    }

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

    double record(std::size_t id) const
    {
        std::size_t iter = iter_size() > 0 ? 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>
    OutputIter read_record(std::size_t id, OutputIter first) const
    {
        VSMC_RUNTIME_ASSERT_CORE_MONITOR_ID(record);

        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;

        return first;
    }

    template <typename OutputIter>
    OutputIter read_record_matrix(MatrixLayout layout, OutputIter first) const
    {
        if (layout == RowMajor)
            return std::copy(record_.begin(), record_.end(), first);

        if (layout == ColMajor) {
            for (std::size_t d = 0; d != dim_; ++d)
                for (std::size_t i = 0; i != iter_size(); ++i)
                    *first++ = record(d, i);
        }

        return first;
    }

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

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

        for (std::size_t d = 0; d != dim_; ++d) {
            for (std::size_t i = 0; i != iter_size(); ++i)
                data[i] = record(d, i);
            if (name_[d].empty())
                df["X." + std::to_string(d)] = data;
            else
                df[name_[d]] = data;
        }

        return df;
    }

    void eval(const eval_type &new_eval, bool record_only = false,
        MonitorStage stage = MonitorMCMC)
    {
        eval_ = new_eval;
        record_only_ = record_only;
        stage_ = stage;
    }

    void operator()(
        std::size_t iter, Particle<T> &particle, MonitorStage stage)
    {
        internal::size_check<VSMC_BLAS_INT>(particle.size(), "Monitor::eval");

        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());
        internal::cblas_dgemv(internal::CblasColMajor, internal::CblasNoTrans,
            static_cast<VSMC_BLAS_INT>(dim_), static_cast<VSMC_BLAS_INT>(N),
            1.0, buffer_.data(), static_cast<VSMC_BLAS_INT>(dim_),
            particle.weight().data(), 1, 0.0, result_.data(), 1);
        push_back(iter);
    }

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

    private:
    std::size_t dim_;
    eval_type eval_;
    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