.. _program_listing_file_lib_prediction.cpp: Program Listing for File prediction.cpp ======================================= |exhale_lsh| :ref:`Return to documentation for file ` (``lib/prediction.cpp``) .. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS .. code-block:: cpp #include "AnalysisGraph.hpp" #include "ExperimentStatus.hpp" #include #include #include #include using namespace std; using Eigen::VectorXd, Eigen::MatrixXd; using fmt::print, fmt::format; namespace rs = ranges; using boost::adaptors::transformed; using fmt::print; /* ============================================================================ Private: Prediction ============================================================================ */ void AnalysisGraph::check_bounds() { int num_verts = this->num_vertices(); for (int v = 0; v < num_verts; v++) { Node &n = (*this)[v]; if (n.has_min && this->current_latent_state[2 * v] < n.min_val) { this->current_latent_state[2 * v] = n.min_val; } if (n.has_max && this->current_latent_state[2 * v] > n.max_val) { this->current_latent_state[2 * v] = n.max_val; } } } void AnalysisGraph::generate_latent_state_sequences( double initial_prediction_step) { // Allocate memory for prediction_latent_state_sequences this->predicted_latent_state_sequences.clear(); this->predicted_latent_state_sequences = vector>( this->res, vector(this->pred_timesteps, VectorXd(this->num_vertices() * 2))); // configure monitoring of experiment progress ExperimentStatus es(this->experiment_id, this->id); double progress_step = 1.0/101.0; // finish with progress at 0.99 es.enter_working_state(); cout << "\nPredicting for " << this->pred_timesteps << " time steps..." << endl; for (int samp : tq::trange(this->res)) { // The sampled transition matrices would be either of matrix exponential // (continuous) version or discretized version depending on whether the // matrix exponential (continuous) version or the discretized transition // matrix version had been used to train the model. This allows us to // make the prediction routine common for both the versions, except for // the exponentiation of the matrix exponential (continuous) transition // matrices. MatrixXd A; MatrixXd &A_samp = this->transition_matrix_collection[samp]; if (this->head_node_model == HeadNodeModel::HNM_FOURIER) { this->merge_concept_LDS_into_seasonal_head_node_modeling_LDS( A_samp, this->initial_latent_state_collection[samp]); } else { this->current_latent_state = this->initial_latent_state_collection[samp]; this->generate_head_node_latent_sequences( samp, initial_prediction_step + this->pred_timesteps); } if (this->continuous) { if (this->head_node_model == HeadNodeModel::HNM_FOURIER) { // Here A = Ac = this->transition_matrix_collection[samp] (continuous) // Evolving the system till the initial_prediction_step A = (this->A_fourier_base * initial_prediction_step).exp(); //this->predicted_latent_state_sequences[samp][0] = this->current_latent_state = A * this->current_latent_state; // After jumping to time step ips - 1, we take one step of length // Δt at a time. So compute the transition matrix for a single // step. Computing the matrix exponential for a Δt time step. // By default we are using Δt = 1 A = e^{Ac * Δt) A = (this->A_fourier_base * this->delta_t).exp(); } else { A = A_samp.exp(); // Evolving the system till the initial_prediction_step for (int ts = 0; ts < initial_prediction_step; ts++) { this->update_latent_state_with_generated_derivatives(ts, ts + 1); // Set derivatives for frozen nodes for (const auto& [v, deriv_func] : this->external_concepts) { const Indicator& ind = this->graph[v].indicators[0]; this->current_latent_state[2 * v + 1] = deriv_func(ts, ind.mean); } this->current_latent_state = A * this->current_latent_state; } this->update_latent_state_with_generated_derivatives( initial_prediction_step, initial_prediction_step + 1); // Set derivatives for frozen nodes for (const auto& [v, deriv_func] : this->external_concepts) { const Indicator& ind = this->graph[v].indicators[0]; this->current_latent_state[2 * v + 1] = deriv_func(initial_prediction_step, ind.mean); } } } else { // Here A = Ad = this->transition_matrix_collection[samp] (discrete) // This is the discrete transition matrix to take a single step of // length Δt // Evolving the system till the initial_prediction_step if (this->head_node_model == HeadNodeModel::HNM_FOURIER) { A = this->A_fourier_base; this->current_latent_state = A.pow(initial_prediction_step) * this->current_latent_state; } else { A = A_samp; for (int ts = 0; ts < initial_prediction_step; ts++) { this->update_latent_state_with_generated_derivatives(ts, ts + 1); // Set derivatives for frozen nodes for (const auto& [v, deriv_func] : this->external_concepts) { const Indicator& ind = this->graph[v].indicators[0]; this->current_latent_state[2 * v + 1] = deriv_func(ts, ind.mean); } this->current_latent_state = A * this->current_latent_state; } this->update_latent_state_with_generated_derivatives( initial_prediction_step, initial_prediction_step + 1); // Set derivatives for frozen nodes for (const auto& [v, deriv_func] : this->external_concepts) { const Indicator& ind = this->graph[v].indicators[0]; this->current_latent_state[2 * v + 1] = deriv_func(initial_prediction_step, ind.mean); } } } this->check_bounds(); this->predicted_latent_state_sequences[samp][0] = this->current_latent_state; // Clear out perpetual constraints residual from previous sample this->perpetual_constraints.clear(); if (this->clamp_at_derivative) { // To clamp a latent state value to x_c at prediction step 1 via // clamping the derivative, we have to perturb the derivative at // prediction step 0, before evolving it to prediction time step 1. // So we have to look one time step ahead whether we have to clamp // at 1. // if (delphi::utils::in(this->one_off_constraints, 1)) { this->perturb_predicted_latent_state_at(1, samp); } } // Since we used ts = 0 for the initial_prediction_step - 1, this loop is // one ahead of the prediction time steps. That is, index 1 in the // prediction data structures is the 0th index for the requested // prediction sequence. Hence, when we are at time step ts, we should // check whether there are any constconstraints // time step index is 1 for (int ts = 1; ts < this->pred_timesteps; ts++) { // When continuous: The standard matrix exponential equation is, // s_{t+Δt} = e^{Ac * Δt } * s_t // Since vector indices are integral values, and in // the implementation s is the vector, to index into // the vector we uses consecutive integers. Thus in // the implementation, the matrix exponential // equation becomes, // s_{t+1} = e^{Ac * Δt } * s_t // What this equation says is that although vector // indices advance by 1, the duration between two // predictions stored in two adjacent vector cells // need not be 1. The time duration is actually Δt. // The actual line of code represents, // s_t = e^{Ac * Δt } * s_{t-1} // When discrete : s_t = Ad * s_{t-1} this->current_latent_state = A * this->predicted_latent_state_sequences[samp][ts - 1]; if (this->head_node_model == HeadNodeModel::HNM_NAIVE) { this->update_latent_state_with_generated_derivatives( initial_prediction_step + ts, initial_prediction_step + ts + 1); } this->check_bounds(); this->predicted_latent_state_sequences[samp][ts] = this->current_latent_state; // Set derivatives for frozen nodes for (const auto & [ v, deriv_func ] : this->external_concepts) { const Indicator& ind = this->graph[v].indicators[0]; this->predicted_latent_state_sequences[samp][ts][2 * v + 1] = deriv_func(initial_prediction_step + ts, ind.mean); } if (this->clamp_at_derivative ) { if (ts == this->rest_derivative_clamp_ts) { // We have perturbed the derivative at ts - 1. If we do not // take any action, that clamped derivative will be in effect // until another clamping or end of prediction. This is where // we take the necessary actions. if (is_one_off_constraints) { // We should revert the derivative to its original value // at this time step so that clamping does not affect // time step ts + 1. for (auto constraint : this->one_off_constraints.at(ts)) { int node_id = constraint.first; this->predicted_latent_state_sequences[samp][ts] (2 * node_id + 1) = this->initial_latent_state_collection[samp] (2 * node_id + 1); } } else { for (auto [node_id, value]: this->one_off_constraints.at(ts)) { this->predicted_latent_state_sequences[samp][ts] (2 * node_id + 1) = 0; } } } // To clamp a latent state value to x_c at prediction step ts + 1 // via clamping the derivative, we have to perturb the derivative // at prediction step ts, before evolving it to prediction time // step ts + 1. So we have to look one time step ahead whether we // have to clamp at ts + 1 and clamp the derivative now. // if (delphi::utils::in(this->one_off_constraints, ts + 1) || !this->perpetual_constraints.empty()) { this->perturb_predicted_latent_state_at(ts + 1, samp); } } else { // Apply constraints to latent state if any // Logic of this condition: // Initially perpetual_constraints = ∅ // // one_off_constraints.at(ts) = ∅ => Unconstrained ∀ ts // // one_off_constraints.at(ts) ‡ ∅ // => ∃ some constraints (But we do not know what kind) // => Perturb latent state // Call perturb_predicted_latent_state_at() // one_off_constraints == true // => We are applying One-off constraints // one_off_constraints == false // => We are applying perpetual constraints // => We add constraints to perpetual_constraints // => perpetual_constraints ‡ ∅ // => The if condition is true ∀ subsequent time steps // after ts (the first time step s.t. // one_off_constraints.at(ts) ‡ ∅ // => Constraints are perpetual // if (delphi::utils::in(this->one_off_constraints, ts) || !this->perpetual_constraints.empty()) { this->perturb_predicted_latent_state_at(ts, samp); } } } // update experiment progress es.increment_progress(progress_step); } // finish experiment progress monitoring es.enter_writing_state(); } /* * Applying constraints (interventions) to latent states * Check the data structure definition to get more descriptions * * Clamping at time step ts for value v * Two methods work in tandem to achieve this. Let us label them as follows: * glss : generate_latent_state_sequences() * ppls@ : perturb_predicted_latent_state_at() * ┍━━━━━━━━━━━━━━━━━━━━━━━━┑ * │ How ┃ * ┝━━━━━━━━━━┳━━━━━━━━━━━━━┫ * │ One-off ┃ Perpetual ┃ * ━━━━━━┳━━━━━━━━━━━━┯━━━━━━━━━━━━┿━━━━━━━━━━┻━━━━━━━━━━━━━┫ * ┃ │ Clamp at │ v - x_{ts-1} ┃ * ┃ │ (by ppls@) │ ts-1 to ──────────── ┃ * ┃ │ │ Δt ┃ * ┃ Derivative ├┈┈┈┈┈┈┈┈┈┈┈┈┼┈┈┈┈┈┈┈┈┈┈┰┈┈┈┈┈┈┈┈┈┈┈┈┈┨ * Where ┃ │ Reset at │ ts to ẋ₀ ┃ ts to 0 ┃ * ┃ │ (by glss) │ from S₀ ┃ ┃ * ┣━━━━━━━━━━━━┿━━━━━━━━━━━━┿━━━━━━━━━━╋━━━━━━━━━━━━━┫ * ┃ Value │ Clamp at │ ts to v ┃ ∀ t≥ts to v ┃ * ┃ │ (by ppls@) │ ┃ ┃ * ━━━━━━┻━━━━━━━━━━━━┷━━━━━━━━━━━━┷━━━━━━━━━━┻━━━━━━━━━━━━━┛ */ void AnalysisGraph::perturb_predicted_latent_state_at(int timestep, int sample_number) { // Let vertices of the CAG be v = 0, 1, 2, 3, ... // Then, // indices 2*v keeps track of the state of each variable v // indices 2*v+1 keeps track of the state of ∂v/∂t if (this->clamp_at_derivative) { for (auto [node_id, value]: this->one_off_constraints.at(timestep)) { // To clamp the latent state value to x_c at time step t via // clamping the derivative, we have to clamp the derivative // appropriately at time step t-1. // Example: // Say we want to clamp the latent state at t=6 to value // x_c (i.e. we want x₆ = x_c. So we have to set the // derivative at t=6-1=5, ẋ₅, as follows: // x_c - x₅ // ẋ₅ = --------- ........... (1) // Δt // x₆ = x₅ + (ẋ₅ × Δt) // = x_c // Thus clamping ẋ₅ (at t = 5) as described in (1) gives // us the desired clamping at t = 5 + 1 = 6 double clamped_derivative = (value - this->predicted_latent_state_sequences[sample_number] [timestep - 1](2 * node_id)) / this->delta_t; this->predicted_latent_state_sequences[sample_number] [timestep - 1](2 * node_id + 1) = clamped_derivative; } // Clamping the derivative at t-1 changes the value at t. // According to our model, derivatives never chance. So if we // do not revert it, clamped derivative stays till another // clamping or the end of prediction. Since this is a one-off // clamping, we have to return the derivative back to its // original value at time step t, before we use it to evolve // time step t + 1. Thus we remember the time step at which we // have to perform this. this->rest_derivative_clamp_ts = timestep; return; } if (this->is_one_off_constraints) { for (auto [node_id, value]: this->one_off_constraints.at(timestep)) { this->predicted_latent_state_sequences[sample_number][timestep](2 * node_id) = value; } } else { // Perpetual constraints if (delphi::utils::in(this->one_off_constraints, timestep)) { // Update any previous perpetual constraints for (auto [node_id, value]: this->one_off_constraints.at(timestep)) { this->perpetual_constraints[node_id] = value; } } // Apply perpetual constraints for (auto [node_id, value]: this->perpetual_constraints) { this->predicted_latent_state_sequences[sample_number][timestep](2 * node_id) = value; } } } void AnalysisGraph::generate_observed_state_sequences() { using rs::to, rs::views::transform; // Allocate memory for observed_state_sequences this->predicted_observed_state_sequences.clear(); this->predicted_observed_state_sequences = vector( this->res, PredictedObservedStateSequence(this->pred_timesteps, vector>())); for (int samp = 0; samp < this->res; samp++) { vector& sample = this->predicted_latent_state_sequences[samp]; this->predicted_observed_state_sequences[samp] = sample | transform([this](VectorXd latent_state) { return this->generate_observed_state(latent_state); }) | to(); } } vector> AnalysisGraph::generate_observed_state(VectorXd latent_state) { using rs::to, rs::views::transform; int num_verts = this->num_vertices(); vector> observed_state(num_verts); for (int v = 0; v < num_verts; v++) { vector& indicators = (*this)[v].indicators; observed_state[v] = vector(indicators.size()); observed_state[v] = indicators | transform([&](Indicator ind) { return ind.mean * latent_state[2 * v]; }) | to(); } return observed_state; } FormattedPredictionResult AnalysisGraph::format_prediction_result() { // NOTE: To facilitate clamping derivatives, we start prediction one time // step before the requested prediction start time. We are omitting // that additional time step from the results returned to the user. this->pred_timesteps--; // Access // [ sample ][ time_step ][ vertex_name ][ indicator_name ] auto result = FormattedPredictionResult( this->res, vector>>( this->pred_timesteps)); for (int samp = 0; samp < this->res; samp++) { // NOTE: To facilitate clamping derivatives, we start prediction one time // step before the requested prediction start time. We are omitting // that additional time step from the results returned to the user. for (int ts = 1; ts <= this->pred_timesteps; ts++) { for (auto [vert_name, vert_id] : this->name_to_vertex) { for (auto [ind_name, ind_id] : (*this)[vert_id].nameToIndexMap) { result[samp][ts - 1][vert_name][ind_name] = this->predicted_observed_state_sequences[samp][ts][vert_id] [ind_id]; } } } } return result; } void AnalysisGraph::run_model(int start_year, int start_month, int end_year, int end_month) { if (!this->trained) { print("Passed untrained Causal Analysis Graph (CAG) Model. \n", "Try calling .train_model(...) first!"); throw "Model not yet trained"; } // Check for sensible prediction time step ranges. // NOTE: To facilitate clamping derivatives, we start prediction one time // step before the requested prediction start time. Therefor, the // possible valid earliest prediction time step is one time step after // the training start time. if (start_year < this->training_range.first.first || (start_year == this->training_range.first.first && start_month <= this->training_range.first.second)) { print("The initial prediction date can't be before the " "initial training date. Defaulting initial prediction date " "to initial training date + 1 month."); start_month = this->training_range.first.second + 1; if (start_month == 13) { start_year = this->training_range.first.first + 1; start_month = 1; } else { start_year = this->training_range.first.first; } } /* * total_timesteps * ____________________________________________ * | | * v v * start training end prediction * |--------------------------------------------| * : |--------------------------------| * : start prediction : * ^ ^ ^ * |___________|________________________________| * diff pred_timesteps */ int total_timesteps = this->calculate_num_timesteps(this->training_range.first.first, this->training_range.first.second, end_year, end_month); this->pred_timesteps = this->calculate_num_timesteps( start_year, start_month, end_year, end_month); int pred_init_timestep = total_timesteps - pred_timesteps; int year = start_year; int month = start_month; this->pred_range.clear(); this->pred_range = vector(this->pred_timesteps); for (int t = 0; t < this->pred_timesteps; t++) { this->pred_range[t] = to_string(year) + "-" + to_string(month); if (month == 12) { year++; month = 1; } else { month++; } } // NOTE: To facilitate clamping derivatives, we start prediction one time // step before the requested prediction start time. When we are // returning results, we have to remove the predictions at the 0th // index of each predicted observed state sequence. // Adding that additional time step. // t = Requested prediction start time step // = pred_init_timestep // t - 1 = Prediction start time step // = pred_init_timestep - 1 pred_init_timestep--; this->pred_timesteps++; this->generate_latent_state_sequences(pred_init_timestep); this->generate_observed_state_sequences(); } void AnalysisGraph::add_constraint(int step, string concept_name, string indicator_name, double indicator_clamp_value) { // When constraining latent state derivatives, to reach the // prescribed value for an observation at projection step t, we have // to clamp the derivative at projection step t-1 appropriately. // Therefor, to constrain at projection step 0, we have to clamp the // derivative at projection step -1. // To facilitate this, we start projection one time step before the // requested projection start time. // To correct for that and make internal vector indexing easier, // nicer and less error prone, we shift all the constraints by one // step up. step++; // Check whether this concept is in the CAG if (!delphi::utils::in(this->name_to_vertex, concept_name)) { print("Concept \"{0}\" not in CAG!\n", concept_name); return; } int concept_id = this->name_to_vertex.at(concept_name); Node& n = (*this)[concept_id]; if (n.indicators.size() <= 0) { // This concept does not have any indicators attached. print("Concept \"{0}\" does not have any indicators attached!\n", concept_name); print("\tCannot set constraint\n"); return; } // If the indicator name is empty we clamp the first indicator attached to // this concept int ind_id = 0; if (!indicator_name.empty()) { if (!delphi::utils::in(this->indicators_in_CAG, indicator_name)) { print("Indicator \"{0}\" is not in CAG!\n", indicator_name); return; } if (!delphi::utils::in(n.nameToIndexMap, indicator_name)) { print("Indicator \"{0}\" is not attached to {1} in CAG!\n", indicator_name, concept_name); return; } ind_id = n.nameToIndexMap.at(indicator_name); } Indicator& ind = n.indicators[ind_id]; // We have to clamp the latent state value corresponding to this // indicator such that the probability where the emission Gaussian // emitting the requested indicator value is the highest. For a // Gaussian emission function, the highest probable value is its // mean. So we have to set: // μ = ind_value // latent_clamp_value * scaling_factor = ind_value // latent_clamp_value = ind_value / scaling_factor // NOTE: In our code we incorrectly call scaling_factor as // indicator mean. To avoid confusion here, I am using the // correct terminology. // (Gosh, when we have incorrect terminology and we know it // and we have not fixed it, I have to type a lot of // comments) double latent_clamp_value = indicator_clamp_value / ind.get_mean(); if (this->head_nodes.find(concept_id) == this->head_nodes.end()) { if (!delphi::utils::in(this->one_off_constraints, step)) { this->one_off_constraints[step] = vector>(); } this->one_off_constraints[step].push_back( make_pair(concept_id, latent_clamp_value)); } else { step += this->pred_start_timestep; if (!delphi::utils::in(this->head_node_one_off_constraints, step)) { this->head_node_one_off_constraints[step] = vector>(); } this->head_node_one_off_constraints[step].push_back( make_pair(concept_id, latent_clamp_value)); } } /* ============================================================================ Public: Prediction ============================================================================ */ Prediction AnalysisGraph::generate_prediction(int start_year, int start_month, int end_year, int end_month, ConstraintSchedule constraints, bool one_off, bool clamp_deri) { this->is_one_off_constraints = one_off; this->clamp_at_derivative = clamp_deri; this->one_off_constraints.clear(); this->head_node_one_off_constraints.clear(); for (auto [step, const_vec] : constraints) { for (auto constraint : const_vec) { string concept_name = get<0>(constraint); string indicator_name = get<1>(constraint); double value = get<2>(constraint); this->add_constraint(step, concept_name, indicator_name, value); } } this->run_model(start_year, start_month, end_year, end_month); return make_tuple( this->training_range, this->pred_range, this->format_prediction_result()); } void AnalysisGraph::generate_prediction(int pred_start_timestep, int pred_timesteps, ConstraintSchedule constraints, bool one_off, bool clamp_deri) { this->is_one_off_constraints = one_off; this->clamp_at_derivative = clamp_deri; this->one_off_constraints.clear(); this->head_node_one_off_constraints.clear(); this->pred_start_timestep = pred_start_timestep - 1; this->pred_timesteps = pred_timesteps + 1; for (auto [step, const_vec] : constraints) { for (auto constraint : const_vec) { string concept_name = get<0>(constraint); string indicator_name = get<1>(constraint); double value = get<2>(constraint); this->add_constraint(step, concept_name, indicator_name, value); } } this->generate_latent_state_sequences(this->pred_start_timestep); this->generate_observed_state_sequences(); } vector> AnalysisGraph::prediction_to_array(string indicator) { int vert_id = -1; int ind_id = -1; auto result = vector>(this->res, vector(this->pred_timesteps)); // Find the vertex id the indicator is attached to and // the indicator id of it. // TODO: We can make this more efficient by making indicators_in_CAG // a map from indicator names to vertices they are attached to. // This is just a quick and dirty implementation for (auto [v_name, v_id] : this->name_to_vertex) { for (auto [i_name, i_id] : (*this)[v_id].nameToIndexMap) { if (indicator.compare(i_name) == 0) { vert_id = v_id; ind_id = i_id; goto indicator_found; } } } // Program will reach here only if the indicator is not found throw IndicatorNotFoundException(format( "AnalysisGraph::prediction_to_array - indicator \"{}\" not found!\n", indicator)); indicator_found: for (int samp = 0; samp < this->res; samp++) { for (int ts = 0; ts < this->pred_timesteps; ts++) { result[samp][ts] = this->predicted_observed_state_sequences[samp][ts][vert_id][ind_id]; } } return result; }