<commit_before>//
//  pfilter.cpp
//  EpiGenMCMC
//
//  Created by <NAME> on 05/05/2016.
//  Copyright (c) 2016 <NAME>, Imperial College London. All rights reserved.
//

#include "pfilter.h"
#include <iostream>



namespace EpiGenPfilter {
    double pfilter(Model & sim_model, Parameter & model_params, MCMCoptions & options, Particle &particles, Trajectory & output_traj, TimeSeriesData &epi_data, TreeData &tree_data, MultiTreeData &multitree_data) {
        double loglik = 0.0;
        int num_groups = options.num_groups;
        int num_particles = options.particles;
        int init_seed = options.seed;
        int total_dt = options.total_dt;
        double sim_dt = options.sim_dt;
        int total_steps = ceil((double)total_dt/(double)options.pfilter_every);
        int add_dt = 0;
        double ESS_threshold = options.pfilter_threshold*(double)num_particles;
        Likelihood likelihood_calc;
//        std::vector <Parameter> values;// (options.num_threads, model_params);
//        for (int i=0; i!=options.num_threads; ++i) values.push_back(model_params);
//        for (int i=0; i!=model_params.get_total_params(); ++i) values.push_back(model_params.get(i));
        std::vector <std::vector<double> > values(options.num_threads, std::vector<double>(model_params.get_total_params(), 0.0));
        for (int i=0; i!=options.num_threads; ++i) {
            for (int j=0; j!=model_params.get_total_params(); ++j) {
                values[i][j] = model_params.get(j);
            }
        }
//        printf("Size of values = %d\n",values.size());
        double reporting_rate = 1.0;
        if (model_params.param_exists("reporting")) {
            reporting_rate = model_params.get("reporting");
        }
        std::vector <std::string> param_names = model_params.get_names_vector();
        if (model_params.param_exists("time_before_data")) {
            add_dt = model_params.get("time_before_data");
        }
        if (options.save_traj) {
            if (add_dt > 0) {
                particles.start_particle_tracing(add_dt+total_dt, num_groups);
            }
            else if (add_dt < 0) {
                particles.start_particle_tracing(add_dt+total_dt, num_groups);
                total_steps = ceil((double)(total_dt+add_dt)/(double)options.pfilter_every);
            }
            else {
                particles.start_particle_tracing(total_dt, num_groups);
            }
        }
        std::vector <Model> models;
        for (int i=0; i<options.num_threads; ++i) {
            models.push_back(sim_model);
        }
        // Simulate model and calculate likelihood assuming no observed data
        if (model_params.param_exists("time_before_data")) {
            if (add_dt > 0) {
                omp_set_num_threads(options.num_threads);
//                std::vector <Trajectory *> curr_trajs;
//                for (int i=0; i!=num_particles; ++i) {
//                    curr_trajs.push_back(particles.get_traj(i));
//                }
#pragma omp parallel for shared(particles, values) schedule(static, 1)
                for (int tn=0; tn<options.num_threads; tn++) {
                    for (int i=tn; i<num_particles; i+=options.num_threads) {
                        // Adjust length of trajectory
                        particles.get_traj(i)->resize(add_dt, num_groups);
                        models[tn].simulate(values[tn], param_names, particles.get_traj(i), 0, add_dt, sim_dt, total_dt, options.rng[tn]);
                        if (options.which_likelihood<2) {
                            double w = likelihood_calc.binomial_lik(reporting_rate, particles.get_traj(i)->get_total_traj(), add_dt+total_dt, 0, add_dt, num_groups, false);
                            particles.set_weight(w, i, false);
                        }
                        if (options.save_traj) {
                            particles.save_traj_to_matrix(i, 0, add_dt);
                            particles.save_ancestry(i, 0, add_dt);
                        }
                    }
                }
            }
        }
        init_seed += num_particles;
        int t=0;
        int start_dt;
        int end_dt;
        for (t=0; t!=total_steps; ++t) {
//            std::vector<double> we(options.particles, 0.0), wg(options.particles, 0.0);
            start_dt = t*options.pfilter_every;
            end_dt = std::min(total_dt, (t+1)*options.pfilter_every);
            omp_set_num_threads(options.num_threads);
#pragma omp parallel for shared (particles, values) schedule(static,1)
            for (int tn=0; tn<options.num_threads; tn++) {
                for (int i=tn; i<num_particles; i+=options.num_threads) {
                    // Adjust length of trajectory
                    particles.get_traj(i)->resize(end_dt-start_dt, options.num_groups);
                    models[tn].simulate(values[tn], param_names, particles.get_traj(i), start_dt, end_dt, sim_dt, total_dt, options.rng[tn]);<commit_msg>latest working version<commit_after>//
//  pfilter.cpp
//  EpiGenMCMC
//
//  Created by <NAME> on 05/05/2016.
//  Copyright (c) 2016 <NAME>, Imperial College London. All rights reserved.
//

#include "pfilter.h"
#include <iostream>



namespace EpiGenPfilter {
    double pfilter(Model & sim_model, Parameter & model_params, MCMCoptions & options, Particle &particles, Trajectory & output_traj, TimeSeriesData &epi_data, TreeData &tree_data, MultiTreeData &multitree_data) {
        double loglik = 0.0;
        int num_groups = options.num_groups;
        int num_particles = options.particles;
        int init_seed = options.seed;
        int total_dt = options.total_dt;
        double sim_dt = options.sim_dt;
        int total_steps = ceil((double)total_dt/(double)options.pfilter_every);
        int add_dt = 0;
        double ESS_threshold = options.pfilter_threshold*(double)num_particles;
        Likelihood likelihood_calc;
//        std::vector <Parameter> values;// (options.num_threads, model_params);
//        for (int i=0; i!=options.num_threads; ++i) values.push_back(model_params);
//        for (int i=0; i!=model_params.get_total_params(); ++i) values.push_back(model_params.get(i));
        std::vector <std::vector<double> > values(options.num_threads, std::vector<double>(model_params.get_total_params(), 0.0));
        for (int i=0; i!=options.num_threads; ++i) {
            for (int j=0; j!=model_params.get_total_params(); ++j) {
                values[i][j] = model_params.get(j);
            }
        }
//        printf("Size of values = %d\n",values.size());
        double reporting_rate = 1.0;
        if (model_params.param_exists("reporting")) {
            reporting_rate = model_params.get("reporting");
        }
        std::vector <std::string> param_names = model_params.get_names_vector();
        if (model_params.param_exists("time_before_data")) {
            add_dt = model_params.get("time_before_data");
        }
        if (options.save_traj) {
            if (add_dt > 0) {
                particles.start_particle_tracing(add_dt+total_dt, num_groups);
            }
            else if (add_dt < 0) {
                particles.start_particle_tracing(add_dt+total_dt, num_groups);
                total_steps = ceil((double)(total_dt+add_dt)/(double)options.pfilter_every);
            }
            else {
                particles.start_particle_tracing(total_dt, num_groups);
            }
        }
        std::vector <Model> models;
        for (int i=0; i<options.num_threads; ++i) {
            models.push_back(sim_model);
        }
        // Simulate model and calculate likelihood assuming no observed data
        if (model_params.param_exists("time_before_data")) {
            if (add_dt > 0) {
                omp_set_num_threads(options.num_threads);
//                std::vector <Trajectory *> curr_trajs;
//                for (int i=0; i!=num_particles; ++i) {
//                    curr_trajs.push_back(particles.get_traj(i));
//                }
#pragma omp parallel for shared(particles, values) schedule(static, 1)
                for (int tn=0; tn<options.num_threads; tn++) {
                    for (int i=tn; i<num_particles; i+=options.num_threads) {
                        // Adjust length of trajectory
                        particles.get_traj(i)->resize(add_dt, num_groups);
                        models[tn].simulate(values[tn], param_names, particles.get_traj(i), 0, add_dt, sim_dt, total_dt, options.rng[tn]);
                        if (options.which_likelihood<2) {
                            double w = likelihood_calc.binomial_lik(reporting_rate, particles.get_traj(i)->get_total_traj(), add_dt+total_dt, 0, add_dt, num_groups, false);
                            particles.set_weight(w, i, false);
                        }
                        if (options.save_traj) {
                            particles.save_traj_to_matrix(i, 0, add_dt);
                            particles.save_ancestry(i, 0, add_dt);
                        }
                    }
                }
            }
        }
        init_seed += num_particles;
        int t=0;
        int start_dt;
        int end_dt;
        for (t=0; t!=total_steps; ++t) {
//            std::vector<double> we(options.particles, 0.0), wg(options.particles, 0.0);
            start_dt = t*options.pfilter_every;
            end_dt = std::min(total_dt, (t+1)*options.pfilter_every);
            omp_set_num_threads(options.num_threads);
#pragma omp parallel for shared (particles, values) schedule(static,1)
            for (int tn=0; tn<options.num_threads; tn++) {
                for (int i=tn; i<num_particles; i+=options.num_threads) {
                    // Adjust length of trajectory
//                    if (tn==0) std::cout << i << ' ' << std::endl;
                    particles.get_traj(i)->resize(end_dt-start_dt, options.num_groups);
                    models[tn].simulate(values[tn], param_names, particles.get_traj(i), start_dt, end_dt, sim_dt, total_dt, options.rng[tn]);
{'fraction_non_alphanumeric': 0.1229802513464991, 'fraction_numerical': 0.008876920007979254, 'mean_word_length': 2.130502653762098, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1017303', 'n_tokens_mistral': 3202, 'n_tokens_neox': 3076, 'n_words': 572}