Skip to content

Commit

Permalink
[core] update, fix #82
Browse files Browse the repository at this point in the history
  • Loading branch information
@milobanks
milobanks committed Jan 15, 2025
1 parent 07037fa commit 98b2add
Show file tree
Hide file tree
Showing 10 changed files with 361 additions and 277 deletions.
6 changes: 3 additions & 3 deletions include/epiworld/agent-meat-virus-sampling.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ namespace sampler {
* @return Virus<TSeq>* of the selected virus. If none selected (or none
* available,) returns a nullptr;
*/
template<typename TSeq>
template<typename TSeq = EPI_DEFAULT_TSEQ>
inline std::function<void(Agent<TSeq>*,Model<TSeq>*)> make_update_susceptible(
std::vector< epiworld_fast_uint > exclude = {}
)
Expand Down Expand Up @@ -179,7 +179,7 @@ inline std::function<void(Agent<TSeq>*,Model<TSeq>*)> make_update_susceptible(
* @return Virus<TSeq>* of the selected virus. If none selected (or none
* available,) returns a nullptr;
*/
template<typename TSeq = int>
template<typename TSeq = EPI_DEFAULT_TSEQ>
inline std::function<Virus<TSeq>*(Agent<TSeq>*,Model<TSeq>*)> make_sample_virus_neighbors(
std::vector< epiworld_fast_uint > exclude = {}
)
Expand Down Expand Up @@ -347,7 +347,7 @@ inline std::function<Virus<TSeq>*(Agent<TSeq>*,Model<TSeq>*)> make_sample_virus_
* @return Virus<TSeq>* of the selected virus. If none selected (or none
* available,) returns a nullptr;
*/
template<typename TSeq = int>
template<typename TSeq = EPI_DEFAULT_TSEQ>
inline Virus<TSeq> * sample_virus_single(Agent<TSeq> * p, Model<TSeq> * m)
{

Expand Down
4 changes: 2 additions & 2 deletions include/epiworld/epiworld.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@

/* Versioning */
#define EPIWORLD_VERSION_MAJOR 0
#define EPIWORLD_VERSION_MINOR 4
#define EPIWORLD_VERSION_MINOR 6
#define EPIWORLD_VERSION_PATCH 0

static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR;
Expand Down Expand Up @@ -88,4 +88,4 @@ namespace epiworld {

}

#endif
#endif
177 changes: 101 additions & 76 deletions include/epiworld/math/lfmcmc/lfmcmc-bones.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,15 +23,15 @@ using LFMCMCKernelFun = std::function<epiworld_double(const std::vector< epiworl

/**
* @brief Proposal function
* @param params_now Vector where to save the new parameters.
* @param params_prev Vector of reference parameters.
* @param new_params Vector where to save the new parameters.
* @param old_params Vector of reference parameters.
* @param m LFMCMC model.
* @tparam TData
*/
template<typename TData>
inline void proposal_fun_normal(
std::vector< epiworld_double >& params_now,
const std::vector< epiworld_double >& params_prev,
std::vector< epiworld_double >& new_params,
const std::vector< epiworld_double >& old_params,
LFMCMC<TData>* m
);

Expand Down Expand Up @@ -60,32 +60,32 @@ inline LFMCMCProposalFun<TData> make_proposal_norm_reflective(
* Proposals are made within a radious 1 of the current
* state of the parameters.
*
* @param params_now Where to write the new parameters
* @param params_prev Reference parameters
* @param new_params Where to write the new parameters
* @param old_params Reference parameters
* @tparam TData
* @param m LFMCMC model.
*/
template<typename TData>
inline void proposal_fun_unif(
std::vector< epiworld_double >& params_now,
const std::vector< epiworld_double >& params_prev,
std::vector< epiworld_double >& new_params,
const std::vector< epiworld_double >& old_params,
LFMCMC<TData>* m
);

/**
* @brief Uses the uniform kernel with euclidean distance
*
* @param stats_now Vector of current statistics based on
* simulated data.
* @param stats_obs Vector of observed statistics
* @param simulated_stats Vector of statistics based on
* simulated data
* @param observed_stats Vector of observed statistics
* @param epsilon Epsilon parameter
* @param m LFMCMC model.
* @param m LFMCMC model
* @return epiworld_double
*/
template<typename TData>
inline epiworld_double kernel_fun_uniform(
const std::vector< epiworld_double >& stats_now,
const std::vector< epiworld_double >& stats_obs,
const std::vector< epiworld_double >& simulated_stats,
const std::vector< epiworld_double >& observed_stats,
epiworld_double epsilon,
LFMCMC<TData>* m
);
Expand All @@ -94,14 +94,17 @@ inline epiworld_double kernel_fun_uniform(
* @brief Gaussian kernel
*
* @tparam TData
* @param epsilon
* @param m
* @param simulated_stats Vector of statistics based on
* simulated data
* @param observed_stats Vector of observed statistics
* @param epsilon Epsilon parameter
* @param m LFMCMC model
* @return epiworld_double
*/
template<typename TData>
inline epiworld_double kernel_fun_gaussian(
const std::vector< epiworld_double >& stats_now,
const std::vector< epiworld_double >& stats_obs,
const std::vector< epiworld_double >& simulated_stats,
const std::vector< epiworld_double >& observed_stats,
epiworld_double epsilon,
LFMCMC<TData>* m
);
Expand All @@ -116,7 +119,7 @@ class LFMCMC {
private:

// Random number sampling
std::mt19937 * engine = nullptr;
std::shared_ptr< std::mt19937 > m_engine = nullptr;

std::shared_ptr< std::uniform_real_distribution<> > runifd =
std::make_shared< std::uniform_real_distribution<> >(0.0, 1.0);
Expand All @@ -128,86 +131,98 @@ class LFMCMC {
std::make_shared< std::gamma_distribution<> >();

// Process data
TData * observed_data;

// Information about the size of the problem
size_t n_samples;
size_t n_statistics;
size_t n_parameters;
TData m_observed_data;
std::vector< TData > * m_simulated_data = nullptr;

epiworld_double epsilon;
// Information about the size of the process
size_t m_n_samples;
size_t m_n_stats;
size_t m_n_params;

std::vector< epiworld_double > params_now;
std::vector< epiworld_double > params_prev;
std::vector< epiworld_double > params_init;
epiworld_double m_epsilon;

std::vector< epiworld_double > observed_stats; ///< Observed statistics
std::vector< epiworld_double > m_initial_params; ///< Initial parameters
std::vector< epiworld_double > m_current_proposed_params; ///< Proposed parameters for the next sample
std::vector< epiworld_double > m_current_accepted_params; ///< Most recently accepted parameters (current state of MCMC)
std::vector< epiworld_double > m_current_proposed_stats; ///< Statistics from simulation run with proposed parameters
std::vector< epiworld_double > m_current_accepted_stats; ///< Statistics from simulation run with most recently accepted params

std::vector< epiworld_double > sampled_params; ///< Sampled Parameters
std::vector< epiworld_double > sampled_stats; ///< Sampled statistics
std::vector< epiworld_double > sampled_stats_prob; ///< Sampled statistics
std::vector< bool > sampled_accepted; ///< Indicator of accepted statistics
std::vector< epiworld_double > m_observed_stats; ///< Observed statistics

std::vector< epiworld_double > accepted_params; ///< Posterior distribution (accepted samples)
std::vector< epiworld_double > accepted_stats; ///< Posterior distribution (accepted samples)
std::vector< epiworld_double > accepted_params_prob; ///< Posterior probability
std::vector< epiworld_double > m_all_sample_params; ///< Parameter samples
std::vector< epiworld_double > m_all_sample_stats; ///< Statistic samples
std::vector< bool > m_all_sample_acceptance; ///< Indicator if sample was accepted
std::vector< epiworld_double > m_all_sample_drawn_prob; ///< Drawn probabilities (runif()) for each sample
std::vector< epiworld_double > m_all_sample_kernel_scores; ///< Kernel scores for each sample

std::vector< epiworld_double > drawn_prob; ///< Drawn probabilities (runif())
std::vector< TData > * sampled_data = nullptr;
std::vector< epiworld_double > m_all_accepted_params; ///< Posterior distribution of parameters from accepted samples
std::vector< epiworld_double > m_all_accepted_stats; ///< Posterior distribution of statistics from accepted samples
std::vector< epiworld_double > m_all_accepted_kernel_scores; ///< Kernel scores for each accepted sample

// Functions
LFMCMCSimFun<TData> simulation_fun;
LFMCMCSummaryFun<TData> summary_fun;
LFMCMCProposalFun<TData> proposal_fun = proposal_fun_normal<TData>;
LFMCMCKernelFun<TData> kernel_fun = kernel_fun_uniform<TData>;
LFMCMCSimFun<TData> m_simulation_fun;
LFMCMCSummaryFun<TData> m_summary_fun;
LFMCMCProposalFun<TData> m_proposal_fun = proposal_fun_normal<TData>;
LFMCMCKernelFun<TData> m_kernel_fun = kernel_fun_uniform<TData>;

// Misc
std::vector< std::string > names_parameters;
std::vector< std::string > names_statistics;
std::vector< std::string > m_param_names;
std::vector< std::string > m_stat_names;

std::chrono::time_point<std::chrono::steady_clock> time_start;
std::chrono::time_point<std::chrono::steady_clock> time_end;
std::chrono::time_point<std::chrono::steady_clock> m_start_time;
std::chrono::time_point<std::chrono::steady_clock> m_end_time;

// Timing
// std::chrono::milliseconds
std::chrono::duration<epiworld_double,std::micro> time_elapsed =
std::chrono::duration<epiworld_double,std::micro> m_elapsed_time =
std::chrono::duration<epiworld_double,std::micro>::zero();

inline void get_elapsed(
inline void get_elapsed_time(
std::string unit,
epiworld_double * last_elapsed,
std::string * unit_abbr,
bool print
);
) const;

void chrono_start();
void chrono_end();

// Progress
bool verbose = true;
Progress progress_bar;

public:

void run(
std::vector< epiworld_double > param_init,
std::vector< epiworld_double > params_init_,
size_t n_samples_,
epiworld_double epsilon_
epiworld_double epsilon_,
int seed = -1
);

LFMCMC() {};
LFMCMC(TData & observed_data_) : observed_data(&observed_data_) {};
LFMCMC(const TData & observed_data_) : m_observed_data(observed_data_) {};
~LFMCMC() {};

void set_observed_data(TData & observed_data_) {observed_data = &observed_data_;};
// Setting LFMCMC variables
void set_observed_data(const TData & observed_data_) {m_observed_data = observed_data_;};

void set_proposal_fun(LFMCMCProposalFun<TData> fun);
void set_simulation_fun(LFMCMCSimFun<TData> fun);
void set_summary_fun(LFMCMCSummaryFun<TData> fun);
void set_kernel_fun(LFMCMCKernelFun<TData> fun);

void set_params_names(std::vector< std::string > names);
void set_stats_names(std::vector< std::string > names);

/**
* @name Random number generation
*
* @param eng
*/
///@{
void set_rand_engine(std::mt19937 & eng);
std::mt19937 & get_rand_endgine();
void set_rand_engine(std::shared_ptr< std::mt19937 > & eng);
std::shared_ptr< std::mt19937 > & get_rand_endgine();
void seed(epiworld_fast_uint s);
void set_rand_gamma(epiworld_double alpha, epiworld_double beta);
epiworld_double runif();
Expand All @@ -219,28 +234,38 @@ class LFMCMC {
///@}

// Accessing parameters of the function
size_t get_n_samples() const {return n_samples;};
size_t get_n_statistics() const {return n_statistics;};
size_t get_n_parameters() const {return n_parameters;};
epiworld_double get_epsilon() const {return epsilon;};

const std::vector< epiworld_double > & get_params_now() {return params_now;};
const std::vector< epiworld_double > & get_params_prev() {return params_prev;};
const std::vector< epiworld_double > & get_params_init() {return params_init;};
const std::vector< epiworld_double > & get_statistics_obs() {return observed_stats;};
const std::vector< epiworld_double > & get_statistics_hist() {return sampled_stats;};
const std::vector< bool > & get_statistics_accepted() {return sampled_accepted;};
const std::vector< epiworld_double > & get_posterior_lf_prob() {return accepted_params_prob;};
const std::vector< epiworld_double > & get_drawn_prob() {return drawn_prob;};
std::vector< TData > * get_sampled_data() {return sampled_data;};

void set_par_names(std::vector< std::string > names);
void set_stats_names(std::vector< std::string > names);
size_t get_n_samples() const {return m_n_samples;};
size_t get_n_stats() const {return m_n_stats;};
size_t get_n_params() const {return m_n_params;};
epiworld_double get_epsilon() const {return m_epsilon;};

const std::vector< epiworld_double > & get_initial_params() const {return m_initial_params;};
const std::vector< epiworld_double > & get_current_proposed_params() const {return m_current_proposed_params;};
const std::vector< epiworld_double > & get_current_accepted_params() const {return m_current_accepted_params;};
const std::vector< epiworld_double > & get_current_proposed_stats() const {return m_current_proposed_stats;};
const std::vector< epiworld_double > & get_current_accepted_stats() const {return m_current_accepted_stats;};

const std::vector< epiworld_double > & get_observed_stats() const {return m_observed_stats;};

const std::vector< epiworld_double > & get_all_sample_params() const {return m_all_sample_params;};
const std::vector< epiworld_double > & get_all_sample_stats() const {return m_all_sample_stats;};
const std::vector< bool > & get_all_sample_acceptance() const {return m_all_sample_acceptance;};
const std::vector< epiworld_double > & get_all_sample_drawn_prob() const {return m_all_sample_drawn_prob;};
const std::vector< epiworld_double > & get_all_sample_kernel_scores() const {return m_all_sample_kernel_scores;};

const std::vector< epiworld_double > & get_all_accepted_params() const {return m_all_accepted_params;};
const std::vector< epiworld_double > & get_all_accepted_stats() const {return m_all_accepted_stats;};
const std::vector< epiworld_double > & get_all_accepted_kernel_scores() const {return m_all_accepted_kernel_scores;};

std::vector< TData > * get_simulated_data() const {return m_simulated_data;};

std::vector< epiworld_double > get_params_mean();
std::vector< epiworld_double > get_stats_mean();
std::vector< epiworld_double > get_mean_params();
std::vector< epiworld_double > get_mean_stats();

void print() ;
// Printing
LFMCMC<TData> & verbose_off();
LFMCMC<TData> & verbose_on();
void print(size_t burnin = 0u) const;

};

Expand Down
Loading

0 comments on commit 98b2add

Please sign in to comment.