pf_base.h

Go to the documentation of this file.

#ifndef PF_BASE_H
#define PF_BASE_H
 
#include <map>
#include <string>
#include <vector>
 
#ifdef DROPPINGTHISINRPACKAGE
    #include <RcppEigen.h>
    // [[Rcpp::depends(RcppEigen)]]
#else
    #include <Eigen/Dense>
#endif
 
 
 
namespace pf {
 
namespace bases {
 
 
 
/************************************************************************************************************/
 
 
template<typename float_t, size_t dimobs, size_t dimstate>
class pf_base {
public:
 
    /* expose float type to users of this ABCTP */
    using float_type = float_t;
 
    /* expose observation-sized vector type to users  */
    using obs_sized_vec = Eigen::Matrix<float_t,dimobs,1>;
 
    /* expose state-sized vector type to users */
    using state_sized_vec = Eigen::Matrix<float_t,dimstate,1>;
    
    /* expose state-sized vector to users */
    using dynamic_matrix = Eigen::Matrix<float_t,Eigen::Dynamic,Eigen::Dynamic>;
    
    /* a function  */
    using func = std::function<const dynamic_matrix(const state_sized_vec&)>;
    
    /* functions  */
    using func_vec = std::vector<func>;
 
    /* the dimension of each observation vector (allows indirect access to template parameters)*/
    static constexpr unsigned int dim_obs = dimobs;
 
    /* the dimension of the state vector (allows indirect access to template parameters)*/
    static constexpr unsigned int dim_state = dimstate;
 
 
    virtual void filter(const obs_sized_vec &data, const func_vec& fs = func_vec() ) = 0;
 
 
    virtual float_t getLogCondLike() const = 0;
 
 
    virtual ~pf_base(){};
};
 
 
/************************************************************************************************************/
 
 
template<typename float_t, size_t dimobs, size_t dimstate, size_t dimcov>
class pf_withcov_base {
public:
 
    /* expose float type to users of this ABCTP */
    using float_type = float_t;
 
    /* expose observation-sized vector type to users  */
    using obs_sized_vec = Eigen::Matrix<float_t,dimobs,1>;
 
    /* expose state-sized vector type to users */
    using state_sized_vec = Eigen::Matrix<float_t,dimstate,1>;
   
    using cov_sized_vec = Eigen::Matrix<float_t,dimcov,1>;
 
    /* expose state-sized vector to users */
    using dynamic_matrix = Eigen::Matrix<float_t,Eigen::Dynamic,Eigen::Dynamic>;
    
    /* a function  */
    using func = std::function<const dynamic_matrix(const state_sized_vec&, 
                                                    const cov_sized_vec&)>;
    
    /* functions  */
    using func_vec = std::vector<func>;
 
    /* the dimension of each observation vector (allows indirect access to template parameters)*/
    static constexpr unsigned int dim_obs = dimobs;
 
    /* the dimension of the state vector (allows indirect access to template parameters)*/
    static constexpr unsigned int dim_state = dimstate;
 
 
    virtual void filter(const obs_sized_vec &data, 
                        const cov_sized_vec &cov, 
                        const func_vec& fs = func_vec() ) = 0;
 
 
    virtual float_t getLogCondLike() const = 0;
 
 
    virtual ~pf_withcov_base(){};
};
 
 
/************************************************************************************************************/
 
template<typename float_t, size_t dim_s_state, size_t dim_ns_state, size_t dimobs>
class rbpf_base {
public:
 
    /* expose observation-sized vector type to users  */
    using obs_sized_vec = Eigen::Matrix<float_t,dimobs,1>;
    
    /* expose sampled-state-size vector type to users */
    using sampled_state_sized_vec  = Eigen::Matrix<float_t,dim_s_state,1>;
    
    /* expose not-sampled-state-sized vector type to users */
    using not_sampled_state_sized_vec = Eigen::Matrix<float_t,dim_ns_state,1>;
    
    /* expose state-sized vector to users */
    using dynamic_matrix = Eigen::Matrix<float_t,Eigen::Dynamic,Eigen::Dynamic>;
    
    /* a function */
    using func  = std::function<const dynamic_matrix(const not_sampled_state_sized_vec&, const sampled_state_sized_vec&)>;
    
    /* functions */
    using func_vec = std::vector<func>;
 
    /* the size of the sampled state portion */
    static constexpr unsigned int dim_sampled_state = dim_s_state;
 
    /* the size of the not sampled state portion */
    static constexpr unsigned int dim_not_sampled_state = dim_ns_state;
 
    /* the size of the observations */
    static constexpr unsigned int dim_obs = dimobs;
 
 
    virtual void filter(const obs_sized_vec &data, const func_vec& fs = func_vec() ) = 0;
 
 
    virtual ~rbpf_base(){};
};
 
 
/************************************************************************************************************/
 
 
template<size_t dimx, size_t dimy, typename float_t>
class ForwardMod {
private:
 
    /* expose observation-sized vector type to users (private because they clash with pf_base) */
    using obs_sized_vec = Eigen::Matrix<float_t,dimy,1>;
 
    /* expose state-sized vector type to users (private because they clash with pf_base) */
    using state_sized_vec = Eigen::Matrix<float_t,dimx,1>;
 
public:
 
    /* a pair of paths (xts first, yts second) */
    using aPair = std::pair<std::vector<state_sized_vec>, std::vector<obs_sized_vec> >;
 
 
    aPair sim_forward(unsigned int T);
 
 
    virtual state_sized_vec muSamp() = 0;
 
 
    virtual state_sized_vec fSamp (const state_sized_vec &xtm1) = 0;
 
 
    virtual obs_sized_vec gSamp (const state_sized_vec &xt)   = 0;
};
 
 
template<size_t dimx, size_t dimy, typename float_t>
auto ForwardMod<dimx,dimy,float_t>::sim_forward(unsigned int T) -> aPair {
 
    std::vector<state_sized_vec> xs;
    std::vector<obs_sized_vec> ys;
 
    // time 1
    xs.push_back(this->muSamp());
    ys.push_back(this->gSamp(xs[0]));
 
    // time > 1
    for(size_t i = 1; i < T; ++i) {
        
        auto xt = this->fSamp(xs[i-1]);
        xs.push_back(xt);
        ys.push_back(this->gSamp(xt));
    }
 
    return aPair(xs, ys); 
}
 
 
/************************************************************************************************************/
 
 
template<size_t dimx, size_t dimy, typename float_t>
class GenForwardMod {
private:
 
    /* expose observation-sized vector type to users (private because they clash with pf_base) */
    using obs_sized_vec = Eigen::Matrix<float_t,dimy,1>;
 
    /* expose state-sized vector type to users (private because they clash with pf_base) */
    using state_sized_vec = Eigen::Matrix<float_t,dimx,1>;
 
public:
 
    /* a pair of paths (xts first, yts second) */
    using aPair = std::pair<std::vector<state_sized_vec>, std::vector<obs_sized_vec> >;
 
 
    aPair sim_forward(unsigned int T);
 
 
    virtual state_sized_vec muSamp() = 0;
 
 
    virtual state_sized_vec fSamp (const state_sized_vec &xtm1, const obs_sized_vec &ytm1) = 0;
 
 
    virtual obs_sized_vec gSamp (const state_sized_vec &xt)   = 0;
};
 
 
template<size_t dimx, size_t dimy, typename float_t>
auto GenForwardMod<dimx,dimy,float_t>::sim_forward(unsigned int T) -> aPair {
 
    std::vector<state_sized_vec> xs;
    std::vector<obs_sized_vec> ys;
 
    // time 1
    xs.push_back(this->muSamp());
    ys.push_back(this->gSamp(xs[0]));
 
    // time > 1
    for(size_t i = 1; i < T; ++i) {
        
        auto xt = this->fSamp(xs[i-1], ys[i-1]);
        xs.push_back(xt);
        ys.push_back(this->gSamp(xt));
    }
 
    return aPair(xs, ys); 
}
 
 
/************************************************************************************************************/
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
class FutureSimulator {
private:
 
    /* expose observation-sized vector type to users  (private because they clash with pf_base) */
    using obs_sized_vec = Eigen::Matrix<float_t,dimy,1>;
 
    /* expose state-sized vector type to users */
    using state_sized_vec = Eigen::Matrix<float_t,dimx,1>;
 
public:
    /* one time point's pair of of nparts states and nparts observations  */
    using timePair = std::pair<std::array<state_sized_vec, nparts>, std::array<obs_sized_vec, nparts> >;
 
    /* many path pairs. time, (state/obs), particle */
    using manyPairs = std::vector<timePair>;
 
    /* observation paths (time, particle) */
    using obsPaths = std::vector<std::array<obs_sized_vec, nparts> >;
 
    /* many state paths (time, particle) */
    using statePaths = std::vector<std::array<state_sized_vec, nparts> >;
 
 
    manyPairs sim_future(unsigned int num_time_steps);
 
 
    obsPaths sim_future_obs(unsigned int num_time_steps);
 
 
    statePaths sim_future_states(unsigned int num_time_steps);
 
 
    virtual std::array<state_sized_vec,nparts> get_uwtd_samps() const = 0;
 
 
    virtual state_sized_vec fSamp (const state_sized_vec &xtm1) = 0;
 
 
    virtual obs_sized_vec gSamp (const state_sized_vec &xt) = 0;
};
 
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
auto FutureSimulator<dimx,dimy,float_t,nparts>::sim_future(unsigned int num_time_steps) -> manyPairs {
 
    // this gets returned
    manyPairs allFutures;
 
    // stuff that gets changed every time loop
    std::array<state_sized_vec, nparts> states;
    std::array<obs_sized_vec, nparts> observations;
    std::array<state_sized_vec,nparts> past_states;
    std::array<state_sized_vec,nparts> first_states = this->get_uwtd_samps();
 
    // iterate over time
    for(unsigned int time = 0; time < num_time_steps; ++time){
 
        // use particle filter's most recent information if you're going one step ahead into the future
        // otherwise use output generated in a previous iteration of this time loop
        if(time == 0){
            past_states = first_states;
        }
 
        // go particle by particle
        for(size_t j = 0; j < nparts; ++j){
            states[j] = this->fSamp( past_states[j] );
            observations[j] = this->gSamp( states[j] );
        }
           
        // add time period content 
        allFutures.push_back(timePair(states, observations)); //TODO make sure deep copy
        past_states = states;
    }
 
    return allFutures; 
}
 
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
auto FutureSimulator<dimx,dimy,float_t,nparts>::sim_future_obs(unsigned int num_time_steps) -> obsPaths {
 
    // this gets returned
    manyPairs obs_and_state_paths = this->sim_future(num_time_steps);
    obsPaths future_obs;
    for(size_t time = 0; time < obs_and_state_paths.size(); ++time){
        future_obs.push_back( std::get<1>(obs_and_state_paths[time]) );
    }
    
    return future_obs;
}
 
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
auto FutureSimulator<dimx,dimy,float_t,nparts>::sim_future_states(unsigned int num_time_steps) -> statePaths {
 
    // this gets returned
    statePaths future_states;
    manyPairs obs_and_state_paths = sim_future(num_time_steps);
    for(size_t time = 0; time < obs_and_state_paths.size(); ++time){
        future_states.push_back( std::get<0>(obs_and_state_paths[time]) );
    }
    
    return future_states;
}
 
 
/************************************************************************************************************/
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
class GenFutureSimulator {
private:
 
    /* expose observation-sized vector type to users  (private because they clash with pf_base) */
    using obs_sized_vec = Eigen::Matrix<float_t,dimy,1>;
 
    /* expose state-sized vector type to users */
    using state_sized_vec = Eigen::Matrix<float_t,dimx,1>;
 
public:
    /* one time point's pair of of nparts states and nparts observations  */
    using timePair = std::pair<std::array<state_sized_vec, nparts>, std::array<obs_sized_vec, nparts> >;
 
    /* many path pairs. time, (state/obs), particle */
    using manyPairs = std::vector<timePair>;
 
    /* observation paths (time, particle) */
    using obsPaths = std::vector<std::array<obs_sized_vec, nparts> >;
 
    /* many state paths (time, particle) */
    using statePaths = std::vector<std::array<state_sized_vec, nparts> >;
 
 
    manyPairs sim_future(unsigned int num_time_steps, const obs_sized_vec &yt);
 
 
    obsPaths sim_future_obs(unsigned int num_time_steps, const obs_sized_vec &yt);
 
 
    statePaths sim_future_states(unsigned int num_time_steps, const obs_sized_vec &yt);
 
 
    virtual std::array<state_sized_vec,nparts> get_uwtd_samps() const = 0;
 
 
    virtual state_sized_vec fSamp (const state_sized_vec &xtm1, const obs_sized_vec &ytm1) = 0;
 
 
    virtual obs_sized_vec gSamp (const state_sized_vec &xt) = 0;
};
 
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
auto GenFutureSimulator<dimx,dimy,float_t,nparts>::sim_future(unsigned int num_time_steps, const obs_sized_vec &yt) -> manyPairs {
 
    // this gets returned
    manyPairs allFutures;
 
    // stuff that gets changed every time loop
    std::array<state_sized_vec, nparts> states;
    std::array<obs_sized_vec, nparts> observations;
    std::array<state_sized_vec,nparts> past_states;
    std::array<obs_sized_vec,nparts> past_obs;
    std::array<state_sized_vec, nparts> first_states = this->get_uwtd_samps();
 
    // iterate over time
    for(unsigned int time = 0; time < num_time_steps; ++time){
 
        // use particle filter's most recent information if you're going one step ahead into the future
        // otherwise use output generated in a previous iteration of this time loop
        if(time == 0){
            past_states = first_states;
            past_obs.fill(yt);
        }
 
        // go particle by particle
        for(size_t j = 0; j < nparts; ++j){
            states[j] = this->fSamp( past_states[j], past_obs[j] );
            observations[j] = this->gSamp( states[j] );
        }
           
        // add time period content 
        allFutures.push_back(timePair(states, observations)); //TODO make sure deep copy
        past_states = states;
        past_obs = observations; 
    }
 
    return allFutures; 
}
 
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
auto GenFutureSimulator<dimx,dimy,float_t,nparts>::sim_future_obs(unsigned int num_time_steps, const obs_sized_vec &yt) -> obsPaths {
 
    // this gets returned
    manyPairs obs_and_state_paths = this->sim_future(num_time_steps, yt);
    obsPaths future_obs;
    for(size_t time = 0; time < obs_and_state_paths.size(); ++time){
        future_obs.push_back( std::get<1>(obs_and_state_paths[time]) );
    }
    
    return future_obs;
}
 
 
template<size_t dimx, size_t dimy, typename float_t, size_t nparts>
auto GenFutureSimulator<dimx,dimy,float_t,nparts>::sim_future_states(unsigned int num_time_steps, const obs_sized_vec &yt) -> statePaths {
 
    // this gets returned
    statePaths future_states;
    manyPairs obs_and_state_paths = sim_future(num_time_steps,yt);
    for(size_t time = 0; time < obs_and_state_paths.size(); ++time){
        future_states.push_back( std::get<0>(obs_and_state_paths[time]) );
    }
    
    return future_states;
}
 
/************************************************************************************************************/
 
 
 
template<size_t dimstate, size_t dimobs, typename float_t>
class cf_filter{
 
public:
    
    /* expose observation-sized vector type to users  */
    using obs_sized_vec = Eigen::Matrix<float_t,dimobs,1>;
 
    /* expose state-sized vector type to users */
    using state_sized_vec = Eigen::Matrix<float_t,dimstate,1>;
    
    virtual ~cf_filter();
    
    
    virtual float_t getLogCondLike() const = 0;
};
 
 
template<size_t dimstate, size_t dimobs, typename float_t>
cf_filter<dimstate,dimobs,float_t>::~cf_filter() {}
 
 
/************************************************************************************************************/
 
 
template<typename float_t, size_t dimobs, size_t dimstate, size_t dimu, size_t dimur, size_t numparts>
class pf_base_crn {
public:
 
    /* expose float type to users of this ABCTP */
    using float_type = float_t;
 
    /* expose observation-sized vector type to users  */
    using obs_sized_vec = Eigen::Matrix<float_t,dimobs,1>;
 
    /* expose state-sized vector type to users */
    using state_sized_vec = Eigen::Matrix<float_t,dimstate,1>;
    
    /* expose state-sized vector to users */
    using dynamic_matrix = Eigen::Matrix<float_t,Eigen::Dynamic,Eigen::Dynamic>;
    
    /* a function  */
    using func = std::function<const dynamic_matrix(const state_sized_vec&)>;
   
    /* functions  */
    using func_vec = std::vector<func>;
 
    /* type for common random numbers used for sampling state proposals*/ 
    using usv = Eigen::Matrix<float_t, dimu, 1>;
 
    /* type for common random numbers used for resampling at a given time point  */
    using usvr = Eigen::Matrix<float_t, dimur, 1>;
 
    /* the dimension of each observation vector (allows indirect access to template parameters)*/
    static constexpr unsigned int dim_obs = dimobs;
 
    /* the dimension of the state vector (allows indirect access to template parameters)*/
    static constexpr unsigned int dim_state = dimstate;
 
 
    virtual void filter(const obs_sized_vec &data, const std::array<usv,numparts>& Us, const usvr& Uresamp, const func_vec& fs = func_vec() ) = 0;
 
 
    virtual float_t getLogCondLike() const = 0;
 
 
    virtual ~pf_base_crn(){};
};
 
} // namespace bases
} // namespace pf
#endif // PF_BASE_H