path.hpp

Go to the documentation of this file.

//============================================================================
// vSMC/include/vsmc/core/path.hpp
//----------------------------------------------------------------------------
//                         vSMC: Scalable Monte Carlo
//----------------------------------------------------------------------------
// Copyright (c) 2013,2014, 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>
#include <vsmc/utility/aligned_memory.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_FUNCTOR(func, caller, name) \
    VSMC_RUNTIME_ASSERT(static_cast<bool>(func),                             \
            ("**Path::"#caller"** INVALID "#name" OBJECT"))                  \

namespace vsmc {

template <typename T>
class Path
{
    public :

    typedef T value_type;
    typedef cxx11::function<double (
            std::size_t, const Particle<T> &, double *)> eval_type;

    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_[0];}

    const std::size_t *integrand_data () const {return &integrand_[0];}

    const std::size_t *grid_data () const {return &grid_[0];}

    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, const Particle<T> &particle)
    {
        if (!recording_)
            return;

        VSMC_RUNTIME_ASSERT_CORE_PATH_FUNCTOR(eval_, eval, EVALUATION);

        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 *const bptr = &buffer_[0];
        const double *const wptr = particle.weight_set().weight_data();
        double grid = eval_(iter, particle, bptr);
        double integrand = math::dot(N, wptr, bptr);
        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_;
    std::vector<std::size_t> index_;
    std::vector<double, AlignedAllocator<double> > integrand_;
    std::vector<double, AlignedAllocator<double> > grid_;
    std::vector<double, AlignedAllocator<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