diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..9498e50 Binary files /dev/null and b/.DS_Store differ diff --git a/.codacy.yaml b/.codacy.yaml deleted file mode 100644 index ccca973..0000000 --- a/.codacy.yaml +++ /dev/null @@ -1,8 +0,0 @@ ---- -engines: - cppcheck: - language: c++ - -# https://stackoverflow.com/questions/36825903/cppcheck-claims-that-a-field-is-not-used-while-it-is-in-another-file -exclude_paths: - - "src/**.hpp" diff --git a/.github/workflows/pip.yaml b/.github/workflows/pip.yaml index 774ba84..48bae90 100644 --- a/.github/workflows/pip.yaml +++ b/.github/workflows/pip.yaml @@ -14,8 +14,8 @@ jobs: strategy: fail-fast: false matrix: - platform: [ubuntu-latest] - python-version: ["3.7", "3.11"] + platform: [ubuntu-22.04] + python-version: ["3.10.4", "3.11.0"] steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml new file mode 100644 index 0000000..e3b6fcb --- /dev/null +++ b/.github/workflows/publish.yaml @@ -0,0 +1,34 @@ +name: "Publish to PyPI" + +on: + release: + types: [published] + +jobs: + publish: + name: Publish to PyPI + runs-on: ubuntu-latest + + steps: + - name: Check out code + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.11" + + - name: Install build tools + run: | + python -m pip install --upgrade pip + pip install build twine + + - name: Build package + run: pip install . + + - name: Publish to PyPI + env: + TWINE_USERNAME: "__token__" + TWINE_PASSWORD: "${{ secrets.PYPI_API_KEY }}" + run: | + twine upload dist/* diff --git a/.github/workflows/sphinx.yaml b/.github/workflows/sphinx.yaml index 6eaf20b..1c82281 100644 --- a/.github/workflows/sphinx.yaml +++ b/.github/workflows/sphinx.yaml @@ -9,12 +9,6 @@ on: types: [published] workflow_dispatch: -# Sets permissions of the GITHUB_TOKEN to allow deployment to GitHub Pages -permissions: - contents: read - pages: write - id-token: write - jobs: build: runs-on: ubuntu-latest @@ -59,16 +53,38 @@ jobs: run: mv README_files docs/build/html # Upload - - name: Upload artifacts + - name: Upload artifact uses: actions/upload-pages-artifact@v3 with: name: github-pages path: docs/build/html/ + post-page-artifact: + if: ${{ github.event_name == 'pull_request' && github.event.pull_request.head.repo.full_name == github.repository }} + needs: build + runs-on: ubuntu-latest + permissions: + contents: read + pull-requests: write + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Post comment preview + uses: CDCgov/cfa-actions/post-artifact@v1.0.0 + with: + artifact-name: github-pages + gh-token: ${{ secrets.GITHUB_TOKEN}} + message: "Thanks for your contribution, ${{ github.actor }}; your `{ artifact-name }` is ready for download [here]({ artifact-url })." + deploy: # Deploy to the github-pages environment # but not on PRs if: ${{ github.event_name != 'pull_request' }} + permissions: + pages: write + id-token: write + environment: name: github-pages url: ${{ steps.deployment.outputs.page_url }} diff --git a/.gitignore b/.gitignore index 99aefa9..f0fc1ff 100644 --- a/.gitignore +++ b/.gitignore @@ -1,17 +1,19 @@ .cache/ .ipynb_checkpoints/ .pytest_cache/ +.ropeproject/ +.nocheck/ _build/ dist/ build/ +dist/ __pycache__/ *.egg-info/ .ropeproject/ compile_commands.json + *.html *.ipynb -.nocheck/ - *_files !README_files diff --git a/CMakeLists.txt b/CMakeLists.txt index b64bc87..b18e872 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,7 @@ include_directories(include/epiworld) python_add_library(_core MODULE src/database.cpp + src/entity.cpp src/main.cpp src/model.cpp src/saver.cpp diff --git a/Vagrantfile b/Vagrantfile new file mode 100644 index 0000000..2061eb7 --- /dev/null +++ b/Vagrantfile @@ -0,0 +1,34 @@ +Vagrant.configure("2") do |config| + config.vm.box = "ubuntu/jammy64" + + config.vm.provider "virtualbox" do |vb| + vb.memory = "2048" + vb.cpus = 2 + end + + config.vm.network "private_network", type: "dhcp" + + config.vm.provision "shell", inline: <<-SHELL + # Update and install prerequisites + sudo apt-get update + sudo apt-get install -y python3 python3-pip python3-venv cmake build-essential + + # Create a Python virtual environment + python3 -m venv /home/vagrant/epiworldPy-env + + # Activate virtual environment and upgrade pip + source /home/vagrant/epiworldPy-env/bin/activate + pip install --upgrade pip setuptools wheel + + # Install dependencies from pyproject.toml (if it exists) + if [ -f /vagrant/pyproject.toml ]; then + pip install . + fi + + # Deactivate the virtual environment + deactivate + SHELL + + # Synchronize the project directory with the VM + config.vm.synced_folder ".", "/vagrant" +end diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml deleted file mode 100644 index c23d550..0000000 --- a/conda.recipe/meta.yaml +++ /dev/null @@ -1,40 +0,0 @@ -package: - name: epiworldpy - version: 0.0.1 - -source: - path: .. - -build: - number: 0 - script: {{ PYTHON }} -m pip install . -vv - -requirements: - build: - - python - - {{ compiler('cxx') }} - - host: - - python - - pip - - scikit-build-core - - pybind11 >=2.10.0 - - run: - - python - - numpy - - -test: - imports: - - epiworldpy - requires: - - pytest - source_files: - - tests - commands: - - pytest tests - -about: - summary: An example project built with pybind11 and scikit-build. - license_file: LICENSE diff --git a/docs/conf.py b/docs/conf.py index 5f6bad7..19da662 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -67,9 +67,9 @@ # built documents. # # The short X.Y version. -version = '0.0.1' +version = '0.6.0-0' # The full version, including alpha/beta/rc tags. -release = '0.0.1' +release = '0.6.0-0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/include/epiworld/agent-meat-virus-sampling.hpp b/include/epiworld/agent-meat-virus-sampling.hpp index d728993..bfe668e 100644 --- a/include/epiworld/agent-meat-virus-sampling.hpp +++ b/include/epiworld/agent-meat-virus-sampling.hpp @@ -20,7 +20,7 @@ namespace sampler { * @return Virus* of the selected virus. If none selected (or none * available,) returns a nullptr; */ -template +template inline std::function*,Model*)> make_update_susceptible( std::vector< epiworld_fast_uint > exclude = {} ) @@ -179,7 +179,7 @@ inline std::function*,Model*)> make_update_susceptible( * @return Virus* of the selected virus. If none selected (or none * available,) returns a nullptr; */ -template +template inline std::function*(Agent*,Model*)> make_sample_virus_neighbors( std::vector< epiworld_fast_uint > exclude = {} ) @@ -347,7 +347,7 @@ inline std::function*(Agent*,Model*)> make_sample_virus_ * @return Virus* of the selected virus. If none selected (or none * available,) returns a nullptr; */ -template +template inline Virus * sample_virus_single(Agent * p, Model * m) { diff --git a/include/epiworld/epiworld.hpp b/include/epiworld/epiworld.hpp index cfdee7a..6483465 100644 --- a/include/epiworld/epiworld.hpp +++ b/include/epiworld/epiworld.hpp @@ -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; @@ -88,4 +88,4 @@ namespace epiworld { } -#endif \ No newline at end of file +#endif diff --git a/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp b/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp index 6c494b9..46146a9 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp +++ b/include/epiworld/math/lfmcmc/lfmcmc-bones.hpp @@ -23,15 +23,15 @@ using LFMCMCKernelFun = std::function 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* m ); @@ -60,32 +60,32 @@ inline LFMCMCProposalFun 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 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* 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 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* m ); @@ -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 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* m ); @@ -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); @@ -128,77 +131,89 @@ 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 simulation_fun; - LFMCMCSummaryFun summary_fun; - LFMCMCProposalFun proposal_fun = proposal_fun_normal; - LFMCMCKernelFun kernel_fun = kernel_fun_uniform; + LFMCMCSimFun m_simulation_fun; + LFMCMCSummaryFun m_summary_fun; + LFMCMCProposalFun m_proposal_fun = proposal_fun_normal; + LFMCMCKernelFun m_kernel_fun = kernel_fun_uniform; // 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 time_start; - std::chrono::time_point time_end; + std::chrono::time_point m_start_time; + std::chrono::time_point m_end_time; + // Timing // std::chrono::milliseconds - std::chrono::duration time_elapsed = + std::chrono::duration m_elapsed_time = std::chrono::duration::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 fun); void set_simulation_fun(LFMCMCSimFun fun); void set_summary_fun(LFMCMCSummaryFun fun); void set_kernel_fun(LFMCMCKernelFun fun); + + void set_params_names(std::vector< std::string > names); + void set_stats_names(std::vector< std::string > names); /** * @name Random number generation @@ -206,8 +221,8 @@ class LFMCMC { * @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(); @@ -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 & verbose_off(); + LFMCMC & verbose_on(); + void print(size_t burnin = 0u) const; }; diff --git a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp index 38cbf56..6a658ab 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp +++ b/include/epiworld/math/lfmcmc/lfmcmc-meat-print.hpp @@ -2,61 +2,83 @@ #define LFMCMC_MEAT_PRINT_HPP template -inline void LFMCMC::print() +inline void LFMCMC::print(size_t burnin) const { - std::vector< epiworld_double > summ_params(n_parameters * 3, 0.0); - std::vector< epiworld_double > summ_stats(n_statistics * 3, 0.0); - for (size_t k = 0u; k < n_parameters; ++k) + // For each statistic or parameter in the model, we print three values: + // - mean, the 2.5% quantile, and the 97.5% quantile + std::vector< epiworld_double > summ_params(m_n_params * 3, 0.0); + std::vector< epiworld_double > summ_stats(m_n_stats * 3, 0.0); + + // Compute the number of samples to use based on burnin rate + size_t n_samples_print = m_n_samples; + if (burnin > 0) + { + if (burnin >= m_n_samples) + throw std::length_error( + "The burnin is greater than or equal to the number of samples." + ); + + n_samples_print = m_n_samples - burnin; + + } + + epiworld_double n_samples_dbl = static_cast< epiworld_double >( + n_samples_print + ); + + // Compute parameter summary values + for (size_t k = 0u; k < m_n_params; ++k) { // Retrieving the relevant parameter - std::vector< epiworld_double > par_i(n_samples); - for (size_t i = 0u; i < n_samples; ++i) + std::vector< epiworld_double > par_i(n_samples_print); + for (size_t i = burnin; i < m_n_samples; ++i) { - par_i[i] = accepted_params[i * n_parameters + k]; - summ_params[k * 3] += par_i[i]/n_samples; + par_i[i-burnin] = m_all_accepted_params[i * m_n_params + k]; + summ_params[k * 3] += par_i[i-burnin]/n_samples_dbl; } // Computing the 95% Credible interval std::sort(par_i.begin(), par_i.end()); summ_params[k * 3 + 1u] = - par_i[std::floor(.025 * static_cast(n_samples))]; + par_i[std::floor(.025 * n_samples_dbl)]; summ_params[k * 3 + 2u] = - par_i[std::floor(.975 * static_cast(n_samples))]; + par_i[std::floor(.975 * n_samples_dbl)]; } - for (size_t k = 0u; k < n_statistics; ++k) + // Compute statistics summary values + for (size_t k = 0u; k < m_n_stats; ++k) { // Retrieving the relevant parameter - std::vector< epiworld_double > stat_k(n_samples); - for (size_t i = 0u; i < n_samples; ++i) + std::vector< epiworld_double > stat_k(n_samples_print); + for (size_t i = burnin; i < m_n_samples; ++i) { - stat_k[i] = accepted_stats[i * n_statistics + k]; - summ_stats[k * 3] += stat_k[i]/n_samples; + stat_k[i-burnin] = m_all_accepted_stats[i * m_n_stats + k]; + summ_stats[k * 3] += stat_k[i-burnin]/n_samples_dbl; } // Computing the 95% Credible interval std::sort(stat_k.begin(), stat_k.end()); - summ_stats[k * 3 + 1u] = - stat_k[std::floor(.025 * static_cast(n_samples))]; + stat_k[std::floor(.025 * n_samples_dbl)]; summ_stats[k * 3 + 2u] = - stat_k[std::floor(.975 * static_cast(n_samples))]; + stat_k[std::floor(.975 * n_samples_dbl)]; } printf_epiworld("___________________________________________\n\n"); printf_epiworld("LIKELIHOOD-FREE MARKOV CHAIN MONTE CARLO\n\n"); - printf_epiworld("N Samples : %ld\n", n_samples); + printf_epiworld("N Samples (total) : %zu\n", m_n_samples); + printf_epiworld("N Samples (after burn-in period) : %zu\n", m_n_samples - burnin); std::string abbr; epiworld_double elapsed; - get_elapsed("auto", &elapsed, &abbr, false); + get_elapsed_time("auto", &elapsed, &abbr, false); printf_epiworld("Elapsed t : %.2f%s\n\n", elapsed, abbr.c_str()); //////////////////////////////////////////////////////////////////////////// @@ -78,10 +100,10 @@ inline void LFMCMC::print() nchar_par_num += 5; // 1 for neg padd, 2 for decimals, 1 the decimal point, and one b/c log(<10) < 1. std::string charlen = std::to_string(nchar_par_num); - if (names_parameters.size() != 0u) + if (m_param_names.size() != 0u) { int nchar_par = 0; - for (auto & n : names_parameters) + for (auto & n : m_param_names) { int tmp_nchar = n.length(); if (nchar_par < tmp_nchar) @@ -96,15 +118,15 @@ inline void LFMCMC::print() std::string(".2f] (initial : % ") + charlen + std::string(".2f)\n"); - for (size_t k = 0u; k < n_parameters; ++k) + for (size_t k = 0u; k < m_n_params; ++k) { printf_epiworld( fmt_params.c_str(), - names_parameters[k].c_str(), + m_param_names[k].c_str(), summ_params[k * 3], summ_params[k * 3 + 1u], summ_params[k * 3 + 2u], - params_init[k] + m_initial_params[k] ); } @@ -117,7 +139,7 @@ inline void LFMCMC::print() std::string(".2f] (initial : % ") + charlen + std::string(".2f)\n"); - for (size_t k = 0u; k < n_parameters; ++k) + for (size_t k = 0u; k < m_n_params; ++k) { printf_epiworld( @@ -126,7 +148,7 @@ inline void LFMCMC::print() summ_params[k * 3], summ_params[k * 3 + 1u], summ_params[k * 3 + 2u], - params_init[k] + m_initial_params[k] ); } @@ -150,10 +172,10 @@ inline void LFMCMC::print() // Figuring out format std::string fmt_stats; - if (names_statistics.size() != 0u) + if (m_stat_names.size() != 0u) { int nchar_stats = 0; - for (auto & n : names_statistics) + for (auto & n : m_stat_names) { int tmp_nchar = n.length(); if (nchar_stats < tmp_nchar) @@ -168,15 +190,15 @@ inline void LFMCMC::print() std::string(".2f] (Observed: % ") + nchar_char + std::string(".2f)\n"); - for (size_t k = 0u; k < n_statistics; ++k) + for (size_t k = 0u; k < m_n_stats; ++k) { printf_epiworld( fmt_stats.c_str(), - names_statistics[k].c_str(), + m_stat_names[k].c_str(), summ_stats[k * 3], summ_stats[k * 3 + 1u], summ_stats[k * 3 + 2u], - observed_stats[k] + m_observed_stats[k] ); } @@ -190,7 +212,7 @@ inline void LFMCMC::print() std::string(".2f] (Observed: % ") + nchar_char + std::string(".2f)\n"); - for (size_t k = 0u; k < n_statistics; ++k) + for (size_t k = 0u; k < m_n_stats; ++k) { printf_epiworld( fmt_stats.c_str(), @@ -198,7 +220,7 @@ inline void LFMCMC::print() summ_stats[k * 3], summ_stats[k * 3 + 1u], summ_stats[k * 3 + 2u], - observed_stats[k] + m_observed_stats[k] ); } diff --git a/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp b/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp index 688f4fd..239f5fa 100755 --- a/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp +++ b/include/epiworld/math/lfmcmc/lfmcmc-meat.hpp @@ -5,20 +5,20 @@ /** * @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 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* m ) { - for (size_t p = 0u; p < m->get_n_parameters(); ++p) - params_now[p] = params_prev[p] + m->rnorm(); + for (size_t p = 0u; p < m->get_n_params(); ++p) + new_params[p] = old_params[p] + m->rnorm(); return; } @@ -44,20 +44,20 @@ inline LFMCMCProposalFun make_proposal_norm_reflective( LFMCMCProposalFun fun = [scale,lb,ub]( - 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* m ) { // Making the proposal - for (size_t p = 0u; p < m->get_n_parameters(); ++p) - params_now[p] = params_prev[p] + m->rnorm() * scale; + for (size_t p = 0u; p < m->get_n_params(); ++p) + new_params[p] = old_params[p] + m->rnorm() * scale; // Checking boundaries epiworld_double d = ub - lb; int odd; epiworld_double d_above, d_below; - for (auto & p : params_now) + for (auto & p : new_params) { // Correcting if parameter goes above the upper bound @@ -84,7 +84,7 @@ inline LFMCMCProposalFun make_proposal_norm_reflective( } #ifdef EPI_DEBUG - for (auto & p : params_now) + for (auto & p : new_params) if (p < lb || p > ub) throw std::range_error("The parameter is out of bounds."); #endif @@ -103,20 +103,20 @@ inline LFMCMCProposalFun 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 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* m ) { - for (size_t p = 0u; p < m->get_n_parameters(); ++p) - params_now[p] = (params_prev[p] + m->runif(-1.0, 1.0)); + for (size_t p = 0u; p < m->get_n_params(); ++p) + new_params[p] = (old_params[p] + m->runif(-1.0, 1.0)); return; } @@ -124,24 +124,24 @@ inline void proposal_fun_unif( /** * @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 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* m ) { epiworld_double ans = 0.0; - for (size_t p = 0u; p < m->get_n_parameters(); ++p) - ans += std::pow(stats_obs[p] - stats_now[p], 2.0); + for (size_t p = 0u; p < m->get_n_params(); ++p) + ans += std::pow(observed_stats[p] - simulated_stats[p], 2.0); return std::sqrt(ans) < epsilon ? 1.0 : 0.0; @@ -153,21 +153,24 @@ 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 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* m ) { epiworld_double ans = 0.0; - for (size_t p = 0u; p < m->get_n_parameters(); ++p) - ans += std::pow(stats_obs[p] - stats_now[p], 2.0); + for (size_t p = 0u; p < m->get_n_params(); ++p) + ans += std::pow(observed_stats[p] - simulated_stats[p], 2.0); return std::exp( -.5 * (ans/std::pow(1 + std::pow(epsilon, 2.0)/3.0, 2.0)) @@ -179,25 +182,25 @@ inline epiworld_double kernel_fun_gaussian( template inline void LFMCMC::set_proposal_fun(LFMCMCProposalFun fun) { - proposal_fun = fun; + m_proposal_fun = fun; } template inline void LFMCMC::set_simulation_fun(LFMCMCSimFun fun) { - simulation_fun = fun; + m_simulation_fun = fun; } template inline void LFMCMC::set_summary_fun(LFMCMCSummaryFun fun) { - summary_fun = fun; + m_summary_fun = fun; } template inline void LFMCMC::set_kernel_fun(LFMCMCKernelFun fun) { - kernel_fun = fun; + m_kernel_fun = fun; } @@ -205,7 +208,8 @@ template inline void LFMCMC::run( std::vector< epiworld_double > params_init_, size_t n_samples_, - epiworld_double epsilon_ + epiworld_double epsilon_, + int seed ) { @@ -213,100 +217,128 @@ inline void LFMCMC::run( chrono_start(); // Setting the baseline parameters of the model - n_samples = n_samples_; - epsilon = epsilon_; - params_init = params_init_; - n_parameters = params_init_.size(); + m_n_samples = n_samples_; + m_epsilon = epsilon_; + m_initial_params = params_init_; + m_n_params = params_init_.size(); + + if (seed >= 0) + this->seed(seed); - params_now.resize(n_parameters); - params_prev.resize(n_parameters); + m_current_proposed_params.resize(m_n_params); + m_current_accepted_params.resize(m_n_params); - if (sampled_data != nullptr) - sampled_data->resize(n_samples); + if (m_simulated_data != nullptr) + m_simulated_data->resize(m_n_samples); - params_prev = params_init; - params_now = params_init; + m_current_accepted_params = m_initial_params; + m_current_proposed_params = m_initial_params; // Computing the baseline sufficient statistics - summary_fun(observed_stats, *observed_data, this); - n_statistics = observed_stats.size(); + m_summary_fun(m_observed_stats, m_observed_data, this); + m_n_stats = m_observed_stats.size(); // Reserving size - drawn_prob.resize(n_samples); - sampled_accepted.resize(n_samples, false); - sampled_stats.resize(n_samples * n_statistics); - sampled_stats_prob.resize(n_samples); + m_current_proposed_stats.resize(m_n_stats); + m_current_accepted_stats.resize(m_n_stats); + m_all_sample_drawn_prob.resize(m_n_samples); + m_all_sample_acceptance.resize(m_n_samples, false); + m_all_sample_params.resize(m_n_samples * m_n_params); + m_all_sample_stats.resize(m_n_samples * m_n_stats); + m_all_sample_kernel_scores.resize(m_n_samples); - accepted_params.resize(n_samples * n_parameters); - accepted_stats.resize(n_samples * n_statistics); - accepted_params_prob.resize(n_samples); + m_all_accepted_params.resize(m_n_samples * m_n_params); + m_all_accepted_stats.resize(m_n_samples * m_n_stats); + m_all_accepted_kernel_scores.resize(m_n_samples); - TData data_i = simulation_fun(params_init, this); + TData data_i = m_simulation_fun(m_initial_params, this); - std::vector< epiworld_double > proposed_stats_i; - summary_fun(proposed_stats_i, data_i, this); - accepted_params_prob[0u] = kernel_fun(proposed_stats_i, observed_stats, epsilon, this); + m_summary_fun(m_current_proposed_stats, data_i, this); + m_all_accepted_kernel_scores[0u] = m_kernel_fun( + m_current_proposed_stats, m_observed_stats, m_epsilon, this + ); // Recording statistics - for (size_t i = 0u; i < n_statistics; ++i) - sampled_stats[i] = proposed_stats_i[i]; + for (size_t i = 0u; i < m_n_stats; ++i) + m_all_sample_stats[i] = m_current_proposed_stats[i]; + + m_current_accepted_stats = m_current_proposed_stats; - for (size_t k = 0u; k < n_statistics; ++k) - accepted_params[k] = proposed_stats_i[k]; + for (size_t k = 0u; k < m_n_params; ++k) + m_all_accepted_params[k] = m_initial_params[k]; + + for (size_t k = 0u; k < m_n_params; ++k) + m_all_sample_params[k] = m_initial_params[k]; - for (size_t i = 1u; i < n_samples; ++i) + // Init progress bar + progress_bar = Progress(m_n_samples, 80); + if (verbose) { + progress_bar.next(); + } + + // Run LFMCMC + for (size_t i = 1u; i < m_n_samples; ++i) { - // Step 1: Generate a proposal and store it in params_now - proposal_fun(params_now, params_prev, this); + // Step 1: Generate a proposal and store it in m_current_proposed_params + m_proposal_fun(m_current_proposed_params, m_current_accepted_params, this); - // Step 2: Using params_now, simulate data - TData data_i = simulation_fun(params_now, this); + // Step 2: Using m_current_proposed_params, simulate data + TData data_i = m_simulation_fun(m_current_proposed_params, this); // Are we storing the data? - if (sampled_data != nullptr) - sampled_data->operator[](i) = data_i; + if (m_simulated_data != nullptr) + m_simulated_data->operator[](i) = data_i; // Step 3: Generate the summary statistics of the data - summary_fun(proposed_stats_i, data_i, this); + m_summary_fun(m_current_proposed_stats, data_i, this); // Step 4: Compute the hastings ratio using the kernel function - epiworld_double hr = kernel_fun(proposed_stats_i, observed_stats, epsilon, this); - sampled_stats_prob[i] = hr; + epiworld_double hr = m_kernel_fun( + m_current_proposed_stats, m_observed_stats, m_epsilon, this + ); + + m_all_sample_kernel_scores[i] = hr; // Storing data - for (size_t k = 0u; k < n_statistics; ++k) - sampled_stats[i * n_statistics + k] = proposed_stats_i[k]; + for (size_t k = 0u; k < m_n_params; ++k) + m_all_sample_params[i * m_n_params + k] = m_current_proposed_params[k]; + + for (size_t k = 0u; k < m_n_stats; ++k) + m_all_sample_stats[i * m_n_stats + k] = m_current_proposed_stats[k]; // Running Hastings ratio - epiworld_double r = runif(); - drawn_prob[i] = r; + epiworld_double r = runif(); + m_all_sample_drawn_prob[i] = r; // Step 5: Update if likely - if (r < std::min(static_cast(1.0), hr / accepted_params_prob[i - 1u])) + if (r < std::min(static_cast(1.0), hr / m_all_accepted_kernel_scores[i - 1u])) { - accepted_params_prob[i] = hr; - sampled_accepted[i] = true; + m_all_accepted_kernel_scores[i] = hr; + m_all_sample_acceptance[i] = true; - for (size_t k = 0u; k < n_statistics; ++k) - accepted_stats[i * n_statistics + k] = - proposed_stats_i[k]; - - params_prev = params_now; + for (size_t k = 0u; k < m_n_stats; ++k) + m_all_accepted_stats[i * m_n_stats + k] = + m_current_proposed_stats[k]; + m_current_accepted_params = m_current_proposed_params; + m_current_accepted_stats = m_current_proposed_stats; } else { - for (size_t k = 0u; k < n_statistics; ++k) - accepted_stats[i * n_statistics + k] = - accepted_stats[(i - 1) * n_statistics + k]; + for (size_t k = 0u; k < m_n_stats; ++k) + m_all_accepted_stats[i * m_n_stats + k] = + m_all_accepted_stats[(i - 1) * m_n_stats + k]; - accepted_params_prob[i] = accepted_params_prob[i - 1u]; + m_all_accepted_kernel_scores[i] = m_all_accepted_kernel_scores[i - 1u]; } - for (size_t k = 0u; k < n_parameters; ++k) - accepted_params[i * n_parameters + k] = params_prev[k]; + for (size_t k = 0u; k < m_n_params; ++k) + m_all_accepted_params[i * m_n_params + k] = m_current_accepted_params[k]; + if (verbose) { + progress_bar.next(); + } } // End timing @@ -318,7 +350,7 @@ inline void LFMCMC::run( template inline epiworld_double LFMCMC::runif() { - return runifd->operator()(*engine); + return runifd->operator()(*m_engine); } template @@ -327,13 +359,13 @@ inline epiworld_double LFMCMC::runif( epiworld_double ub ) { - return runifd->operator()(*engine) * (ub - lb) + lb; + return runifd->operator()(*m_engine) * (ub - lb) + lb; } template inline epiworld_double LFMCMC::rnorm() { - return rnormd->operator()(*engine); + return rnormd->operator()(*m_engine); } template @@ -342,13 +374,13 @@ inline epiworld_double LFMCMC::rnorm( epiworld_double sd ) { - return (rnormd->operator()(*engine)) * sd + mean; + return (rnormd->operator()(*m_engine)) * sd + mean; } template inline epiworld_double LFMCMC::rgamma() { - return rgammad->operator()(*engine); + return rgammad->operator()(*m_engine); } template @@ -362,7 +394,7 @@ inline epiworld_double LFMCMC::rgamma( rgammad->param(std::gamma_distribution<>::param_type(alpha, beta)); - epiworld_double ans = rgammad->operator()(*engine); + epiworld_double ans = rgammad->operator()(*m_engine); rgammad->param(old_param); @@ -373,14 +405,14 @@ inline epiworld_double LFMCMC::rgamma( template inline void LFMCMC::seed(epiworld_fast_uint s) { - this->engine->seed(s); + this->m_engine->seed(s); } template -inline void LFMCMC::set_rand_engine(std::mt19937 & eng) +inline void LFMCMC::set_rand_engine(std::shared_ptr< std::mt19937 > & eng) { - engine = ŋ + m_engine = eng; } template @@ -390,9 +422,9 @@ inline void LFMCMC::set_rand_gamma(epiworld_double alpha, epiworld_double } template -inline std::mt19937 & LFMCMC::get_rand_endgine() +inline std::shared_ptr< std::mt19937 > & LFMCMC::get_rand_endgine() { - return *engine; + return m_engine; } // Step 1: Simulate data @@ -405,16 +437,16 @@ inline std::mt19937 & LFMCMC::get_rand_endgine() #define DURCAST(tunit,txtunit) {\ elapsed = std::chrono::duration_cast(\ - time_end - time_start).count(); \ + m_end_time - m_start_time).count(); \ abbr_unit = txtunit;} template -inline void LFMCMC::get_elapsed( +inline void LFMCMC::get_elapsed_time( std::string unit, epiworld_double * last_elapsed, std::string * unit_abbr, bool print -) { +) const { // Preparing the result epiworld_double elapsed; @@ -425,7 +457,7 @@ inline void LFMCMC::get_elapsed( { size_t tlength = std::to_string( - static_cast(floor(time_elapsed.count())) + static_cast(floor(m_elapsed_time.count())) ).length(); if (tlength <= 1) @@ -470,46 +502,46 @@ inline void LFMCMC::get_elapsed( template inline void LFMCMC::chrono_start() { - time_start = std::chrono::steady_clock::now(); + m_start_time = std::chrono::steady_clock::now(); } template inline void LFMCMC::chrono_end() { - time_end = std::chrono::steady_clock::now(); - time_elapsed += (time_end - time_start); + m_end_time = std::chrono::steady_clock::now(); + m_elapsed_time += (m_end_time - m_start_time); } template -inline void LFMCMC::set_par_names(std::vector< std::string > names) +inline void LFMCMC::set_params_names(std::vector< std::string > names) { - if (names.size() != n_parameters) + if (names.size() != m_n_params) throw std::length_error("The number of names to add differs from the number of parameters in the model."); - names_parameters = names; + m_param_names = names; } template inline void LFMCMC::set_stats_names(std::vector< std::string > names) { - if (names.size() != n_statistics) + if (names.size() != m_n_stats) throw std::length_error("The number of names to add differs from the number of statistics in the model."); - names_statistics = names; + m_stat_names = names; } template -inline std::vector< epiworld_double > LFMCMC::get_params_mean() +inline std::vector< epiworld_double > LFMCMC::get_mean_params() { - std::vector< epiworld_double > res(this->n_parameters, 0.0); + std::vector< epiworld_double > res(this->m_n_params, 0.0); - for (size_t k = 0u; k < n_parameters; ++k) + for (size_t k = 0u; k < m_n_params; ++k) { - for (size_t i = 0u; i < n_samples; ++i) - res[k] += (this->accepted_params[k + n_parameters * i])/ - static_cast< epiworld_double >(n_samples); + for (size_t i = 0u; i < m_n_samples; ++i) + res[k] += (this->m_all_accepted_params[k + m_n_params * i])/ + static_cast< epiworld_double >(m_n_samples); } return res; @@ -517,19 +549,33 @@ inline std::vector< epiworld_double > LFMCMC::get_params_mean() } template -inline std::vector< epiworld_double > LFMCMC::get_stats_mean() +inline std::vector< epiworld_double > LFMCMC::get_mean_stats() { - std::vector< epiworld_double > res(this->n_statistics, 0.0); + std::vector< epiworld_double > res(this->m_n_stats, 0.0); - for (size_t k = 0u; k < n_statistics; ++k) + for (size_t k = 0u; k < m_n_stats; ++k) { - for (size_t i = 0u; i < n_samples; ++i) - res[k] += (this->accepted_stats[k + n_statistics * i])/ - static_cast< epiworld_double >(n_samples); + for (size_t i = 0u; i < m_n_samples; ++i) + res[k] += (this->m_all_accepted_stats[k + m_n_stats * i])/ + static_cast< epiworld_double >(m_n_samples); } return res; } +template +inline LFMCMC & LFMCMC::verbose_off() +{ + verbose = false; + return *this; +} + +template +inline LFMCMC & LFMCMC::verbose_on() +{ + verbose = true; + return *this; +} + #endif \ No newline at end of file diff --git a/include/epiworld/misc.hpp b/include/epiworld/misc.hpp index dbeb73d..6e0376b 100644 --- a/include/epiworld/misc.hpp +++ b/include/epiworld/misc.hpp @@ -124,7 +124,7 @@ inline bool IN(const Ta & a, const std::vector< Ta > & b) noexcept * @return int If -1 then it means that none got sampled, otherwise the index * of the entry that got drawn. */ -template +template inline int roulette( const std::vector< TDbl > & probs, Model * m diff --git a/include/epiworld/model-bones.hpp b/include/epiworld/model-bones.hpp index cf454d4..d405509 100644 --- a/include/epiworld/model-bones.hpp +++ b/include/epiworld/model-bones.hpp @@ -140,7 +140,7 @@ class Model { std::vector< Entity > entities = {}; std::vector< Entity > entities_backup = {}; - std::mt19937 engine; + std::shared_ptr< std::mt19937 > engine = std::make_shared< std::mt19937 >(); std::uniform_real_distribution<> runifd = std::uniform_real_distribution<> (0.0, 1.0); @@ -286,8 +286,8 @@ class Model { * @param s Seed */ ///@{ - 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(size_t s); void set_rand_norm(epiworld_double mean, epiworld_double sd); void set_rand_unif(epiworld_double a, epiworld_double b); diff --git a/include/epiworld/model-meat.hpp b/include/epiworld/model-meat.hpp index f04f79d..4916016 100644 --- a/include/epiworld/model-meat.hpp +++ b/include/epiworld/model-meat.hpp @@ -674,12 +674,6 @@ inline void Model::agents_empty_graph( } -// template -// inline void Model::set_rand_engine(std::mt19937 & eng) -// { -// engine = std::make_shared< std::mt19937 >(eng); -// } - template inline void Model::set_rand_gamma(epiworld_double alpha, epiworld_double beta) { @@ -720,7 +714,7 @@ template inline epiworld_double & Model::operator()(std::string pname) { if (parameters.find(pname) == parameters.end()) - throw std::range_error("The parameter "+ pname + "is not in the model."); + throw std::range_error("The parameter '"+ pname + "' is not in the model."); return parameters[pname]; @@ -820,7 +814,7 @@ inline void Model::set_backup() // } template -inline std::mt19937 & Model::get_rand_endgine() +inline std::shared_ptr< std::mt19937 > & Model::get_rand_endgine() { return engine; } @@ -828,86 +822,86 @@ inline std::mt19937 & Model::get_rand_endgine() template inline epiworld_double Model::runif() { // CHECK_INIT() - return runifd(engine); + return runifd(*engine); } template inline epiworld_double Model::runif(epiworld_double a, epiworld_double b) { // CHECK_INIT() - return runifd(engine) * (b - a) + a; + return runifd(*engine) * (b - a) + a; } template inline epiworld_double Model::rnorm() { // CHECK_INIT() - return rnormd(engine); + return rnormd(*engine); } template inline epiworld_double Model::rnorm(epiworld_double mean, epiworld_double sd) { // CHECK_INIT() - return rnormd(engine) * sd + mean; + return rnormd(*engine) * sd + mean; } template inline epiworld_double Model::rgamma() { - return rgammad(engine); + return rgammad(*engine); } template inline epiworld_double Model::rgamma(epiworld_double alpha, epiworld_double beta) { auto old_param = rgammad.param(); rgammad.param(std::gamma_distribution<>::param_type(alpha, beta)); - epiworld_double ans = rgammad(engine); + epiworld_double ans = rgammad(*engine); rgammad.param(old_param); return ans; } template inline epiworld_double Model::rexp() { - return rexpd(engine); + return rexpd(*engine); } template inline epiworld_double Model::rexp(epiworld_double lambda) { auto old_param = rexpd.param(); rexpd.param(std::exponential_distribution<>::param_type(lambda)); - epiworld_double ans = rexpd(engine); + epiworld_double ans = rexpd(*engine); rexpd.param(old_param); return ans; } template inline epiworld_double Model::rlognormal() { - return rlognormald(engine); + return rlognormald(*engine); } template inline epiworld_double Model::rlognormal(epiworld_double mean, epiworld_double shape) { auto old_param = rlognormald.param(); rlognormald.param(std::lognormal_distribution<>::param_type(mean, shape)); - epiworld_double ans = rlognormald(engine); + epiworld_double ans = rlognormald(*engine); rlognormald.param(old_param); return ans; } template inline int Model::rbinom() { - return rbinomd(engine); + return rbinomd(*engine); } template inline int Model::rbinom(int n, epiworld_double p) { auto old_param = rbinomd.param(); rbinomd.param(std::binomial_distribution<>::param_type(n, p)); - epiworld_double ans = rbinomd(engine); + epiworld_double ans = rbinomd(*engine); rbinomd.param(old_param); return ans; } template inline void Model::seed(size_t s) { - this->engine.seed(s); + this->engine->seed(s); } template @@ -1297,7 +1291,7 @@ inline Model & Model::run( this->ndays = ndays; if (seed >= 0) - engine.seed(seed); + engine->seed(seed); array_double_tmp.resize(std::max( size(), diff --git a/include/epiworld/models/surveillance.hpp b/include/epiworld/models/surveillance.hpp index 6dca7e4..8b1a288 100644 --- a/include/epiworld/models/surveillance.hpp +++ b/include/epiworld/models/surveillance.hpp @@ -227,7 +227,7 @@ inline ModelSURV::ModelSURV( // How many will we find std::binomial_distribution<> bdist(m->size(), m->par("Surveilance prob.")); - int nsampled = bdist(m->get_rand_endgine()); + int nsampled = bdist(*m->get_rand_endgine()); int to_go = nsampled + 1; diff --git a/include/epiworld/random_graph.hpp b/include/epiworld/random_graph.hpp index dfe6883..ee13b5f 100644 --- a/include/epiworld/random_graph.hpp +++ b/include/epiworld/random_graph.hpp @@ -5,7 +5,7 @@ class RandGraph { private: std::shared_ptr< std::mt19937 > engine; - std::shared_ptr< std::uniform_real_distribution<> > unifd; + std::uniform_real_distribution<> unifd; int N = 0; bool initialized = false; @@ -14,7 +14,7 @@ class RandGraph { RandGraph(int N_) : N(N_) {}; void init(int s); - void set_rand_engine(std::mt19937 & e); + void set_rand_engine(std::shared_ptr< std::mt19937 > & e); epiworld_double runif(); }; @@ -27,18 +27,15 @@ inline void RandGraph::init(int s) { engine->seed(s); - if (!unifd) - unifd = std::make_shared< std::uniform_real_distribution<> >(0, 1); - initialized = true; } -inline void RandGraph::set_rand_engine(std::mt19937 & e) +inline void RandGraph::set_rand_engine(std::shared_ptr< std::mt19937 > & e) { - engine = std::make_shared< std::mt19937 >(e); + engine = e; } @@ -47,7 +44,7 @@ inline epiworld_double RandGraph::runif() { if (!initialized) throw std::logic_error("The object has not been initialized"); - return (*unifd)(engine); + return unifd(*engine); } diff --git a/noxfile.py b/noxfile.py deleted file mode 100644 index 050d60a..0000000 --- a/noxfile.py +++ /dev/null @@ -1,23 +0,0 @@ -from __future__ import annotations - -import nox - -nox.options.sessions = ["lint", "tests"] - - -@nox.session -def lint(session: nox.Session) -> None: - """ - Run the linter. - """ - session.install("pre-commit") - session.run("pre-commit", "run", "--all-files", *session.posargs) - - -@nox.session -def tests(session: nox.Session) -> None: - """ - Run the unit and regular tests. - """ - session.install(".[test]") - session.run("pytest", *session.posargs) diff --git a/pyproject.toml b/pyproject.toml index 57e6993..77338b9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,10 +1,10 @@ [build-system] -requires = ["scikit-build-core>=0.3.3", "pybind11"] +requires = ["scikit-build-core>=0.10", "pybind11"] build-backend = "scikit_build_core.build" [project] name = "epiworldpy" -version = "0.0.1" +version = "0.6.0-0" description = "Python bindings for epiworld" readme = "README.md" requires-python = ">=3.7" diff --git a/src/entity.cpp b/src/entity.cpp new file mode 100644 index 0000000..b31ab89 --- /dev/null +++ b/src/entity.cpp @@ -0,0 +1,32 @@ +#include "entity.hpp" + +#include +#include +#include + +using namespace epiworld; +using namespace epiworldpy; +using namespace pybind11::literals; +namespace py = pybind11; + +static epiworld::Entity new_entity(std::string name, + epiworld::EntityToAgentFun fun) { + Entity entity(name, fun); + + return entity; +} + +static epiworld::EntityToAgentFun new_entity_to_agent_fun( + const std::function &, Model *)> &fun) { + return EntityToAgentFun(fun); +} + +void epiworldpy::export_entity( + pybind11::class_, + std::shared_ptr>> &c) { + c.def(py::init(&new_entity), "Add a list of user data.", py::arg("name"), + py::arg("fun") = nullptr) + .def_static("new_entity_to_agent_fun", &new_entity_to_agent_fun, + "Create a new entity to agent function from a lambda.", + py::arg("fun")); +} diff --git a/src/entity.hpp b/src/entity.hpp new file mode 100644 index 0000000..864894b --- /dev/null +++ b/src/entity.hpp @@ -0,0 +1,11 @@ +#ifndef EPIWORLDPY_ENTITY_HPP +#define EPIWORLDPY_ENTITY_HPP + +#include "epiworld-common.hpp" + +namespace epiworldpy { +void export_entity(pybind11::class_, + std::shared_ptr>> &c); +}; + +#endif /* EPIWORLDPY_ENTITY_HPP */ diff --git a/src/epiworldpy/__init__.py b/src/epiworldpy/__init__.py index ad24a5b..e3e7624 100644 --- a/src/epiworldpy/__init__.py +++ b/src/epiworldpy/__init__.py @@ -2,9 +2,10 @@ from ._core import __doc__, __version__, Model, ModelDiffNet, ModelSEIR, \ ModelSEIRCONN, ModelSEIRD, ModelSIR, ModelSIRCONN, ModelSIRD, \ - ModelSIRDCONN, ModelSIS, ModelSISD, ModelSURV, Saver, Tool, UpdateFun, \ - Virus + ModelSIRDCONN, ModelSIS, ModelSISD, ModelSURV, Entity, Saver, Tool, \ + UpdateFun, Virus __all__ = ["__doc__", "__version__", "ModelDiffNet", "ModelSEIR", "ModelSEIRCONN", "ModelSEIRD", "ModelSIR", "ModelSIRCONN", "ModelSIRD", "ModelSIRDCONN", - "ModelSIS", "ModelSISD", "ModelSURV", "Saver", "Tool", "UpdateFun", "Virus"] + "ModelSIS", "ModelSISD", "ModelSURV", "Entity" "Saver", "Tool", "UpdateFun", + "Virus"] diff --git a/src/main.cpp b/src/main.cpp index 3240e3b..8c53f54 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,19 +1,15 @@ -#include -#include #include -#include -#include #include -#include // for py::array_t +#include #include #include // silently fails when removed. #include "database.hpp" -#include "epiworld-common.hpp" #include "misc.hpp" #include "model-bones.hpp" #include "model.hpp" +#include "entity.hpp" #include "saver.hpp" #include "tool.hpp" #include "virus.hpp" @@ -55,6 +51,8 @@ PYBIND11_MODULE(_core, m) { m, "Model", "A generic model of some kind; a parent class."); auto database = py::class_, std::shared_ptr>>( m, "DataBase", "A container for data generated by a model run."); + auto entity = py::class_, std::shared_ptr>>( + m, "Entity", "Unknown."); auto saver = py::class_>( m, "Saver", "Saves the result of multiple runs."); @@ -66,6 +64,7 @@ PYBIND11_MODULE(_core, m) { epiworldpy::export_model(model); epiworldpy::export_all_models(m); epiworldpy::export_database(database); + epiworldpy::export_entity(entity); epiworldpy::export_saver(saver); epiworldpy::export_tool(tool); epiworldpy::export_virus(virus); diff --git a/src/model.cpp b/src/model.cpp index cce604e..79c1f0b 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -130,6 +130,8 @@ void epiworldpy::export_model(py::class_> &c) { "Adds a virus to the model", py::arg("virus")) .def("add_tool", &Model::add_tool, "Adds a tool to modify the model.", py::arg("tool")) + .def("add_entity", &Model::add_entity) + .def("get_entity", &Model::get_entity) .def("agents_smallworld", &Model::agents_smallworld, "Populate the model without an edgelist.", py::arg("n"), py::arg("k"), py::arg("d"), py::arg("p")) diff --git a/src/virus.cpp b/src/virus.cpp index 2c9004c..0b80f0f 100644 --- a/src/virus.cpp +++ b/src/virus.cpp @@ -1,5 +1,4 @@ #include "virus.hpp" -#include "virus-distribute-meat.hpp" #include #include diff --git a/tests/test_basic.py b/tests/test_basic.py index a51982d..e2a5e46 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -1,4 +1,4 @@ -import epiworldpy as epiworld +from importlib.metadata import version def test_version(): - assert epiworld.__version__ == "0.0.1" + assert version('epiworldpy') == "0.6.0.post0" diff --git a/tests/test_distribution.py b/tests/test_distribution.py new file mode 100644 index 0000000..6e35643 --- /dev/null +++ b/tests/test_distribution.py @@ -0,0 +1,13 @@ +import epiworldpy as epiworld + +def test_distribution(): + hypothetical = epiworld.ModelSIR( + name = 'hypothetical', + prevalence = 0.01, + transmission_rate = 0.1, + recovery_rate = 0.14 + ) + + hypothetical.agents_smallworld(100000, 10, False, 0.01) + + hypothetical.run(100, 223) diff --git a/tests/test_entity.py b/tests/test_entity.py new file mode 100644 index 0000000..486a000 --- /dev/null +++ b/tests/test_entity.py @@ -0,0 +1,30 @@ +import epiworldpy as epiworld + +def dist_factory(start, end): + def entity_to_agent_fun(entity, model): + agents = model.get_agents() + + for i in range(start, end): + entity.add_agent(agents[i], model) + + return entity_to_agent_fun + +def test_entity(): + hypothetical = epiworld.ModelSIRCONN( + name = 'hypothetical', + n = 10000, + prevalence = 0.01, + contact_rate = 10, + transmission_rate = 0.1, + recovery_rate = 0.14 + ) + + entity_1 = epiworld.Entity("Entity 1", dist_factory(0, 3000)) + entity_2 = epiworld.Entity("Entity 2", dist_factory(3000, 6000)) + entity_3 = epiworld.Entity("Entity 3", dist_factory(6000, 10000)) + + hypothetical.add_entity(entity_1) + hypothetical.add_entity(entity_2) + hypothetical.add_entity(entity_3) + + hypothetical.run(100, 223) diff --git a/tests/test_models.py b/tests/test_models.py index 0ad0a93..7bf487b 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -1,75 +1,38 @@ import epiworldpy as epiworld -# TODO: Decide the right way to validate data; how often is the underlying C++ -# simulation code changed in a way that will change outputs here? What's the -# tradeoff? - -def test_diffnet_simple(): - # TODO: How do we invoke this model? - pass - -def test_seir_simple(): - # TODO: Implement `agents_from_adjlist` or something similar. - pass - def test_seirconn_simple(): - covid19 = epiworld.ModelSEIRCONN( - name = 'covid-19', + hypothetical = epiworld.ModelSEIRCONN( + name = 'hypothetical', n = 10000, - prevalence = .01, + prevalence = 0.01, contact_rate = 2.0, - transmission_rate = .1, + transmission_rate = 0.1, incubation_days = 7.0, recovery_rate = 0.14 ) - covid19.run(100, 223) - -def test_seird_simple(): - # TODO: Implement `agents_from_adjlist` or something similar. - pass - -def test_sir_simple(): - # TODO: Implement `agents_from_adjlist` or something similar. - pass + hypothetical.run(100, 223) def test_sirconn_simple(): - covid19 = epiworld.ModelSIRCONN( - name = 'covid-19', + hypothetical = epiworld.ModelSIRCONN( + name = 'hypothetical', n = 10000, - prevalence = .01, + prevalence = 0.01, contact_rate = 2.0, - transmission_rate = .1, + transmission_rate = 0.1, recovery_rate = 0.14 ) - covid19.run(100, 223) - -def test_sird_simple(): - # TODO: Implement `agents_from_adjlist` or something similar. - pass + hypothetical.run(100, 223) -def test_sirdconn_simple(): - covid19 = epiworld.ModelSIRDCONN( - name = 'covid-19', - n = 10000, - prevalence = .01, - contact_rate = 2.0, - transmission_rate = .1, - recovery_rate = 0.14, - death_rate = 0.1, +def test_smallworld(): + hypothetical = epiworld.ModelSIR( + name = 'hypothetical', + prevalence = 0.01, + transmission_rate = 0.1, + recovery_rate = 0.14 ) - covid19.run(100, 223) - -def test_sis_simple(): - # TODO: Implement `agents_from_adjlist` or something similar. - pass - -def test_sisd_simple(): - # TODO: Implement `agents_from_adjlist` or something similar. - pass + hypothetical.agents_smallworld(100000, 10, False, 0.01) -def test_surv_simple(): - # TODO: Implement `agents_from_adjlist` or something similar. - pass + hypothetical.run(100, 223) diff --git a/tests/test_saver.py b/tests/test_saver.py deleted file mode 100644 index a234b47..0000000 --- a/tests/test_saver.py +++ /dev/null @@ -1,18 +0,0 @@ -import epiworldpy as epiworld - -def test_saver_basic(): - covid19 = epiworld.ModelSEIRCONN( - name = 'covid-19', - n = 10000, - prevalence = .01, - contact_rate = 2.0, - transmission_rate = .1, - incubation_days = 7.0, - recovery_rate = 0.14 - ) - - saver = epiworld.Saver("total_hist", "virus_hist") - - saver.run_multiple(covid19, 100, 4, nthreads=1) - - # TODO: Verify things worked correctly, as is the point of tesing.