path.hpp

Go to the documentation of this file.

//============================================================================
// vSMC/include/vsmc/core/path.hpp
//----------------------------------------------------------------------------
//                         vSMC: Scalable Monte Carlo
//----------------------------------------------------------------------------
// Copyright (c) 2013-2015, 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_PATH_HPP
#define VSMC_CORE_PATH_HPP

#include <vsmc/internal/common.hpp>

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

#define VSMC_RUNTIME_ASSERT_CORE_PATH_EVAL                                    \
    VSMC_RUNTIME_ASSERT(                                                      \
        static_cast<bool>(eval_), "**Path::eval** INVALID EVALUAITON OBJECT")

namespace vsmc
{

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

    explicit Path(const eval_type &eval, bool record_only = false)
        : eval_(eval)
        , recording_(true)
        , record_only_(record_only)
        , log_zconst_(0)
    {
    }

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

    void reserve(std::size_t num)
    {
        index_.reserve(num);
        integrand_.reserve(num);
        grid_.reserve(num);
    }

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

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

        return index_[iter];
    }

    double integrand(std::size_t iter) const
    {
        VSMC_RUNTIME_ASSERT_CORE_PATH_ITER(integrand);

        return integrand_[iter];
    }

    double grid(std::size_t iter) const
    {
        VSMC_RUNTIME_ASSERT_CORE_PATH_ITER(grid);

        return grid_[iter];
    }

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

    const std::size_t *integrand_data() const { return integrand_.data(); }

    const std::size_t *grid_data() const { return grid_.data(); }

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

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

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

    void set_eval(const eval_type &new_eval, bool record_only = false)
    {
        eval_ = new_eval;
        record_only_ = record_only;
    }

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

        VSMC_RUNTIME_ASSERT_CORE_PATH_EVAL;

        if (record_only_) {
            double integrand = 0;
            double grid = eval_(iter, particle, &integrand);
            push_back(iter, grid, integrand);

            return;
        }

        const std::size_t N = static_cast<std::size_t>(particle.size());
        buffer_.resize(N);
        double grid = eval_(iter, particle, buffer_.data());
        double integrand =
            dot(N, particle.weight().data(), 1, buffer_.data(), 1);
        push_back(iter, grid, integrand);
    }

    double zconst() const { return std::exp(log_zconst_); }

    double log_zconst() const { return log_zconst_; }

    void clear()
    {
        log_zconst_ = 0;
        index_.clear();
        integrand_.clear();
        grid_.clear();
    }

    bool recording() const { return recording_; }

    void turn_on() { recording_ = true; }

    void turn_off() { recording_ = false; }

    private:
    eval_type eval_;
    bool recording_;
    bool record_only_;
    double log_zconst_;
    Vector<std::size_t> index_;
    Vector<double> integrand_;
    Vector<double> grid_;
    Vector<double> buffer_;

    void push_back(std::size_t iter, double grid, double integrand)
    {
        index_.push_back(iter);
        grid_.push_back(grid);
        integrand_.push_back(integrand);
        if (iter_size() > 1) {
            std::size_t i = iter_size() - 1;
            log_zconst_ += 0.5 * (grid_[i] - grid_[i - 1]) *
                (integrand_[i] + integrand_[i - 1]);
        }
    }
}; // class PathSampling

} // namespace vsmc

#endif // VSMC_CORE_PATH_HPP