// =============================================================================
// ACAS_sim.cpp - Simulation Engine for Asymptomatic Carotid Artery Stenosis
// 
// Model Features:
// - Individual patient simulation with stenosis progression
// - Stroke risk prediction using Ferket meta-analysis model
// - Screening cascade: DUS followed by confirmatory MRA
// - Revascularization complications and outcomes
// - Modified Rankin Scale (mRS) disability tracking (0-6, 6=death)
// - Cost and quality-adjusted life year (QALY) calculations
// - Common Random Numbers (CRN) for variance reduction
// 
// Health States: H=healthy, T=post-TIA, I=post-IS, B=post-both (IS and ICH)
// Stenosis Categories: 0=<50%, 1=50-69%, 2=70-89%, 3=90-99%, 4=occluded
// =============================================================================

#include <Rcpp.h>
#include <iostream>
#include <cmath>
#include <cerrno>
#include <cfenv>
#include <cstring>
#include <algorithm>
#include <fstream>

using namespace Rcpp;

// === UTILITY FUNCTIONS ===

// [[Rcpp::export]]
double RateToProb(double rate, double t){  // Convert annual rates to probabilities
  return(1-exp(-rate*t));
}

// [[Rcpp::export]]
double ProbToRate(double prob, double t){  // Convert probabilities to annual rates
  return(-log(1-prob)/t);
}

// [[Rcpp::export]]
double ProbToProb(double prob_t, double t) {  // Convert probability over time t to annual
  double P1 = 1 - pow((1 - prob_t), 1/t);
  return P1;
}

// [[Rcpp::export]]
double OR_to_RR(double OR, double prob){  // Convert OR to RR given baseline probability
  return(OR/(1-prob + prob*OR));
}

// === mRS MAPPING FUNCTIONS ===

// [[Rcpp::export]]
int map_OTT_mRS_tPA(int mRS_in, double t_OTT){  // Map OTT to mRS for tPA patients
  
  double OR_mRS0       = 1.735*exp(-0.0009*t_OTT);
  double odds_mRS0     = OR_mRS0 * 0.174;
  double prob_mRS0     = odds_mRS0/(1 + odds_mRS0);
  
  double OR_mRS1       = 2.12312*exp(-0.00199*t_OTT);
  double odds_mRS1     = OR_mRS1 * 0.53987;
  double prob_mRS0_1   = odds_mRS1/(1 + odds_mRS1);
  //double prob_mRS1   = prob_mRS0_1 - prob_mRS0;
  
  double OR_mRS2       = 1.8749*exp(-0.00184*t_OTT);
  double odds_mRS2     = OR_mRS2 * 0.9075;
  double prob_mRS0_2   = odds_mRS2/(1 + odds_mRS2);
  //double prob_mRS2   = prob_mRS0_2 - prob_mRS0_1;
  
  double OR_mRS3       = 1.9405*exp(-0.00194*t_OTT);
  double odds_mRS3     = OR_mRS3 * 1.6006;
  double prob_mRS0_3   = odds_mRS3/(1 + odds_mRS3);
  //double prob_mRS3   = prob_mRS0_3 - prob_mRS0_2;
  
  double OR_mRS4       = 1.53133*exp(-0.00185*t_OTT);
  double odds_mRS4     = OR_mRS4 * 3.918;
  double prob_mRS0_4   = odds_mRS4/(1 + odds_mRS4);
  //double prob_mRS4   = prob_mRS0_4 - prob_mRS0_3;
  
  double OR_mRS5       = 1.35215*exp(-0.00158*t_OTT);
  double odds_mRS5     = OR_mRS5 * 7.537;
  double prob_mRS0_5   = odds_mRS5/(1 + odds_mRS5);
  //double prob_mRS5   = prob_mRS0_5 - prob_mRS0_4;
  
  //double prob_mRS6   = 1 - (prob_mRS0 + prob_mRS1 + prob_mRS2 + prob_mRS3 + prob_mRS4 + prob_mRS5);
  
  int mRS_out          = -999;
  double randnum       = R::unif_rand();
  
  //adjust for input mRS
  NumericVector delta_prob = NumericVector::create(0, prob_mRS0, prob_mRS0_1, prob_mRS0_2, prob_mRS0_3, prob_mRS0_4, prob_mRS0_5);
  
  if(randnum < (prob_mRS0 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in)) ){
    mRS_out = 0;
  }else if(randnum < (prob_mRS0_1 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 1;
  }else if(randnum < (prob_mRS0_2 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 2;
  }else if(randnum < (prob_mRS0_3 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 3;
  }else if(randnum < (prob_mRS0_4 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 4;
  }else if(randnum < (prob_mRS0_5 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 5;
  }else{
    mRS_out = 6;
  }
  
  return(mRS_out);
}

// [[Rcpp::export]]
int map_OTT_mRS_notPA(int mRS_in, NumericVector p_vec_mRS_IS_notx){  // mRS for non-tPA patients
  double prob_mRS0   = p_vec_mRS_IS_notx[0];
  double prob_mRS1   = p_vec_mRS_IS_notx[1];
  double prob_mRS2   = p_vec_mRS_IS_notx[2];
  double prob_mRS3   = p_vec_mRS_IS_notx[3];
  double prob_mRS4   = p_vec_mRS_IS_notx[4];
  double prob_mRS5   = p_vec_mRS_IS_notx[5];
  
  double prob_mRS0_1 = prob_mRS0   + prob_mRS1;
  double prob_mRS0_2 = prob_mRS0_1 + prob_mRS2;
  double prob_mRS0_3 = prob_mRS0_2 + prob_mRS3;
  double prob_mRS0_4 = prob_mRS0_3 + prob_mRS4;
  double prob_mRS0_5 = prob_mRS0_4 + prob_mRS5;
  
  int mRS_out        = -999;
  double randnum     = R::unif_rand();
  
  //adjust for input mRS
  NumericVector delta_prob = NumericVector::create(0, prob_mRS0, prob_mRS0_1, prob_mRS0_2, prob_mRS0_3, prob_mRS0_4, prob_mRS0_5);
  
  if(randnum < (prob_mRS0 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in)) ){
    mRS_out = 0;
  }else if(randnum < (prob_mRS0_1 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 1;
  }else if(randnum < (prob_mRS0_2 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 2;
  }else if(randnum < (prob_mRS0_3 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 3;
  }else if(randnum < (prob_mRS0_4 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 4;
  }else if(randnum < (prob_mRS0_5 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 5;
  }else{
    mRS_out = 6;
  }
  
  return(mRS_out);
}

// [[Rcpp::export]]
int map_OTT_mRS_ICH(int mRS_in, NumericVector p_vec_mRS_ICH){  // mRS after ICH
  double prob_mRS0   = p_vec_mRS_ICH[0];
  double prob_mRS1   = p_vec_mRS_ICH[1];
  double prob_mRS2   = p_vec_mRS_ICH[2];
  double prob_mRS3   = p_vec_mRS_ICH[3];
  double prob_mRS4   = p_vec_mRS_ICH[4];
  double prob_mRS5   = p_vec_mRS_ICH[5];
  
  double prob_mRS0_1 = prob_mRS0   + prob_mRS1;
  double prob_mRS0_2 = prob_mRS0_1 + prob_mRS2;
  double prob_mRS0_3 = prob_mRS0_2 + prob_mRS3;
  double prob_mRS0_4 = prob_mRS0_3 + prob_mRS4;
  double prob_mRS0_5 = prob_mRS0_4 + prob_mRS5;
  
  int mRS_out        = -999;
  double randnum     = R::unif_rand();
  
  //adjust for input mRS
  NumericVector delta_prob = NumericVector::create(0, prob_mRS0, prob_mRS0_1, prob_mRS0_2, prob_mRS0_3, prob_mRS0_4, prob_mRS0_5);
  
  if(randnum < (prob_mRS0 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in)) ){
    mRS_out = 0;
  }else if(randnum < (prob_mRS0_1 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 1;
  }else if(randnum < (prob_mRS0_2 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 2;
  }else if(randnum < (prob_mRS0_3 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 3;
  }else if(randnum < (prob_mRS0_4 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 4;
  }else if(randnum < (prob_mRS0_5 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 5;
  }else{
    mRS_out = 6;
  }
  
  return(mRS_out);
}

// [[Rcpp::export]]
int map_OTT_mRS_ICH_CRN(int mRS_in, NumericVector p_vec_mRS_ICH, double randnum){  // ICH mRS with CRN
  double prob_mRS0   = p_vec_mRS_ICH[0];
  double prob_mRS1   = p_vec_mRS_ICH[1];
  double prob_mRS2   = p_vec_mRS_ICH[2];
  double prob_mRS3   = p_vec_mRS_ICH[3];
  double prob_mRS4   = p_vec_mRS_ICH[4];
  double prob_mRS5   = p_vec_mRS_ICH[5];
  
  double prob_mRS0_1 = prob_mRS0   + prob_mRS1;
  double prob_mRS0_2 = prob_mRS0_1 + prob_mRS2;
  double prob_mRS0_3 = prob_mRS0_2 + prob_mRS3;
  double prob_mRS0_4 = prob_mRS0_3 + prob_mRS4;
  double prob_mRS0_5 = prob_mRS0_4 + prob_mRS5;
  
  int mRS_out        = -999;
  
  //adjust for input mRS
  NumericVector delta_prob = NumericVector::create(0, prob_mRS0, prob_mRS0_1, prob_mRS0_2, prob_mRS0_3, prob_mRS0_4, prob_mRS0_5);
  
  if(randnum < (prob_mRS0 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in)) ){
    mRS_out = 0;
  }else if(randnum < (prob_mRS0_1 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 1;
  }else if(randnum < (prob_mRS0_2 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 2;
  }else if(randnum < (prob_mRS0_3 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 3;
  }else if(randnum < (prob_mRS0_4 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 4;
  }else if(randnum < (prob_mRS0_5 - delta_prob(mRS_in))/(1 - delta_prob(mRS_in))){
    mRS_out = 5;
  }else{
    mRS_out = 6;
  }
  
  return(mRS_out);
}

// === ACUTE ISCHEMIC STROKE MODULE ===

// [[Rcpp::export]]
List mod_acuteIS(double age, int mRS_in,  // Simulate acute IS event with tPA eligibility
                 double p_OTD_doc, double OTD_meanlog, double OTD_sdlog,
                 double DTN_tPA_shape, double DTN_tPA_scale, double p_tPA_contra, double p_tPA_other,
                 NumericVector p_vec_mRS_IS_notx,
                 double r_recur_postIS_1yr, NumericVector hr_vec_recur_postIS,
                 NumericVector c_vec_acuteIS_notx_age18_44, NumericVector c_vec_acuteIS_notx_age45_64,
                 NumericVector c_vec_acuteIS_notx_age65_80, NumericVector c_vec_acuteIS_notx_age81up,
                 double c_tPA,
                 NumericVector c_vec_mRS,
                 NumericVector u_vec_mRS
){
  //Output
  int mRS;
  double c_acuteIS = 0;
  double u_acuteIS = 0;
  bool event_rstroke = 0;
  
  //Onset to door
  double t_OTD = R::rlnorm(OTD_meanlog, OTD_sdlog);
  
  //Door to needle
  double t_DTN = R::rweibull(DTN_tPA_shape, DTN_tPA_scale);
  
  //Onset to treatment
  double t_OTT = t_OTD + t_DTN;
  
  //Getting tPA or not
  bool tPA_elig = 0;
  if(R::unif_rand() < p_OTD_doc * (1 - p_tPA_contra)){tPA_elig = 1;}
  
  bool tPA_luck = 0;
  if(R::unif_rand() < 1 - p_tPA_other){tPA_luck = 1;}
  
  bool tPA = 0;
  if((t_OTT <= 4.5*60) & tPA_elig & tPA_luck){tPA = 1;}
  
  if(tPA){
    mRS    = map_OTT_mRS_tPA(mRS_in, t_OTT);
    c_acuteIS += c_tPA;
  }else{
    mRS    = map_OTT_mRS_notPA(mRS_in, p_vec_mRS_IS_notx);
  }
  
  //Add hospitalization cost excluding tPA
  if(age < 45){
    c_acuteIS += c_vec_acuteIS_notx_age18_44(mRS);
  }else if(age >= 45 & age < 65){
    c_acuteIS += c_vec_acuteIS_notx_age45_64(mRS);
  }else if(age >= 65 & age < 80){
    c_acuteIS += c_vec_acuteIS_notx_age65_80(mRS);
  }else{
    c_acuteIS += c_vec_acuteIS_notx_age81up(mRS);
  }
  
  if(mRS < 6){
    //Recurrent IS within a year?
    if(R::unif_rand() < RateToProb(r_recur_postIS_1yr * hr_vec_recur_postIS[mRS] / hr_vec_recur_postIS[3], 1)){
      //Yes
      event_rstroke = 1;
      //Cycle rewards before recurrent stroke
      u_acuteIS += 0.5 *u_vec_mRS(mRS);
      c_acuteIS += 0.25*c_vec_mRS(mRS);
      
      //----- Recurrent stroke
      
      //Onset to door
      double t_OTD = R::rlnorm(OTD_meanlog, OTD_sdlog);
      //Door to needle
      double t_DTN = R::rweibull(DTN_tPA_shape, DTN_tPA_scale);
      //OTT
      double t_OTT = t_OTD + t_DTN;
      
      //Getting tPA or not
      bool tPA_elig = 0;
      if(R::unif_rand() < p_OTD_doc * (1 - p_tPA_contra)){tPA_elig = 1;}
      bool tPA_luck = 0;
      if(R::unif_rand() < 1 - p_tPA_other){tPA_luck = 1;}
      bool tPA = 0;
      if((t_OTT <= 4.5*60) & tPA_elig & tPA_luck){tPA = 1;}
      
      if(tPA){
        mRS    = map_OTT_mRS_tPA(mRS, t_OTT);
        c_acuteIS += c_tPA;
      }else{
        mRS    = map_OTT_mRS_notPA(mRS, p_vec_mRS_IS_notx);
      }
      
      //Add hospitalization cost ex tPA
      if(age < 45){
        c_acuteIS += c_vec_acuteIS_notx_age18_44(mRS);
      }else if(age >= 45 & age < 65){
        c_acuteIS += c_vec_acuteIS_notx_age45_64(mRS);
      }else if(age >= 65 & age < 80){
        c_acuteIS += c_vec_acuteIS_notx_age65_80(mRS);
      }else{
        c_acuteIS += c_vec_acuteIS_notx_age81up(mRS);
      }
      
      //Cycle Rewards
      u_acuteIS += 0.5 *u_vec_mRS(mRS);
      c_acuteIS += 0.25*c_vec_mRS(mRS);
    }else{
      //No recurrent stroke
      c_acuteIS += 0.75 *c_vec_mRS(mRS);
      u_acuteIS += u_vec_mRS(mRS);
    }
  }
  
  
  return(List::create(
      Named("mRS") = mRS,
      Named("c_acuteIS") = c_acuteIS,
      Named("u_acuteIS") = u_acuteIS,
      Named("event_rstroke") = event_rstroke
  )
  );
}


// [[Rcpp::export]]
List mod_acuteIS_all(int n_pop, double p_OTD_doc, double OTD_meanlog, double OTD_sdlog,  // Pre-compute IS outcomes
                     double DTN_tPA_shape, double DTN_tPA_scale, double p_tPA_contra, double p_tPA_other,
                     NumericVector p_vec_mRS_IS_notx,
                     double r_recur_postIS_1yr, NumericVector hr_vec_recur_postIS,
                     NumericVector c_vec_acuteIS_notx_age18_44, NumericVector c_vec_acuteIS_notx_age45_64,
                     NumericVector c_vec_acuteIS_notx_age65_80, NumericVector c_vec_acuteIS_notx_age81up,
                     double c_tPA,
                     NumericVector c_vec_mRS,
                     NumericVector u_vec_mRS
){
  // output
  NumericVector age_init_vec(n_pop*240);
  NumericVector mRS_init_vec(n_pop*240);
  NumericVector mRS_output_vec(n_pop*240);
  NumericVector c_acuteIS_vec(n_pop*240);
  NumericVector u_acuteIS_vec(n_pop*240);
  NumericVector event_rstroke_vec(n_pop*240);
  
  for(double age_mod = 40; age_mod < 80; age_mod++){
    for (int mRS = 0; mRS < 6; mRS++){
      for (int id = 0; id < n_pop; id++){
        int id_all = n_pop * (40 * mRS + (age_mod - 40)) + id;
        age_init_vec(id_all) = age_mod;
        mRS_init_vec(id_all) = mRS;
        
        List output_acuteIS = mod_acuteIS(age_mod, mRS,
                                          p_OTD_doc, OTD_meanlog, OTD_sdlog,
                                          DTN_tPA_shape, DTN_tPA_scale,
                                          p_tPA_contra, p_tPA_other,
                                          p_vec_mRS_IS_notx,
                                          r_recur_postIS_1yr, hr_vec_recur_postIS,
                                          c_vec_acuteIS_notx_age18_44, c_vec_acuteIS_notx_age45_64,
                                          c_vec_acuteIS_notx_age65_80, c_vec_acuteIS_notx_age81up,
                                          c_tPA, c_vec_mRS, u_vec_mRS);
        
        mRS_output_vec(id_all) = output_acuteIS["mRS"];
        c_acuteIS_vec(id_all) = output_acuteIS["c_acuteIS"];
        u_acuteIS_vec(id_all) = output_acuteIS["u_acuteIS"];
        event_rstroke_vec(id_all) = output_acuteIS["event_rstroke"];
        
      }
    }
  }
  
  return(List::create(
      Named("age_init")      = age_init_vec,
      Named("mRS_init")      = mRS_init_vec,
      Named("mRS_output")    = mRS_output_vec,
      Named("c_acuteIS")     = c_acuteIS_vec,
      Named("u_acuteIS")     = u_acuteIS_vec,
      Named("event_rstroke") = event_rstroke_vec
  )
  );
  
}

// [[Rcpp::export]]
List mod_acuteIS_select(double prob_loc, double age, int mRS_in,  // Select pre-computed IS outcome
                        NumericVector c_acuteIS_vec, 
                        NumericVector u_acuteIS_vec,
                        NumericVector mRS_prop_acuteIS_vec
){
  // output
  double c_acuteIS = 0;
  double u_acuteIS = 0;
  
  double age_group = 0;
  if (age >= 45) {age_group = 1;}
  if (age >= 65) {age_group = 2;}
  
  // every age group: 27. 
  // for x = mRS_in,
  // start at c(0,7,13,18,22,25) -- ((7+8-x)*x)/2, get 7-x
  
  //double prob_loc = R::unif_rand();
  int final_count = 0;
  
  for (int count = 0; count < 7 - mRS_in; count++) {
    if (prob_loc <= mRS_prop_acuteIS_vec[27*age_group + (15 - mRS_in)*mRS_in / 2 + count]){
      final_count = count;
      break;
    }
  }
  
  c_acuteIS = c_acuteIS_vec[27*age_group + (15 - mRS_in)*mRS_in / 2 + final_count];
  u_acuteIS = u_acuteIS_vec[27*age_group + (15 - mRS_in)*mRS_in / 2 + final_count];
  //mRS = mRS_out_acuteIS_vec[27*age_group + (15 - mRS_in)*mRS_in / 2 + final_count];
  
  return(List::create(
      Named("mRS") = final_count + mRS_in,
      Named("c_acuteIS") = c_acuteIS,
      Named("u_acuteIS") = u_acuteIS
  )
  );
}

// === STROKE RISK PREDICTION ===

// [[Rcpp::export]]
double predict_stroke_update(NumericVector coefficients,  // Ferket model with time-varying coefficients 
                             NumericVector coef_Ferket_stroke_risk_time_regression_vec,
                             double years_elapsed,
                             double male, double RIDAGEYR, 
                             double race, double BPXSY1, 
                             double SMQ040, double SBPtxt,
                             double diabetes_rollup, double vascular_history){
  
  // Ensure there are enough coefficients
  if (coefficients.size() != 9 || coef_Ferket_stroke_risk_time_regression_vec.size() != 20) {
    stop("Incorrect number of coefficients provided.");
  }
  
  // regression coef
  // double coef_Ferket_intercept = coef_Ferket_stroke_risk_time_regression_vec(0);
  double coef_Ferket_years_elapsed = coef_Ferket_stroke_risk_time_regression_vec(1);
  // double coef_Ferket_age_baseline = coef_Ferket_stroke_risk_time_regression_vec(2);
  // double coef_Ferket_age_baseline_squared = coef_Ferket_stroke_risk_time_regression_vec(3);
  // double coef_Ferket_male = coef_Ferket_stroke_risk_time_regression_vec(4);
  // double coef_Ferket_sbp_baseline = coef_Ferket_stroke_risk_time_regression_vec(5);
  // double coef_Ferket_sbp_baseline_squared = coef_Ferket_stroke_risk_time_regression_vec(6);
  // double coef_Ferket_cvd_history_baseline = coef_Ferket_stroke_risk_time_regression_vec(7);
  // double coef_Ferket_htnmed_baseline = coef_Ferket_stroke_risk_time_regression_vec(8);
  // double coef_Ferket_dm_baseline = coef_Ferket_stroke_risk_time_regression_vec(9);
  // double coef_Ferket_cursmk_baseline = coef_Ferket_stroke_risk_time_regression_vec();
  // double coef_Ferket_race_black = coef_Ferket_stroke_risk_time_regression_vec(10);
  double coef_Ferket_years_elapsed_age_baseline = coef_Ferket_stroke_risk_time_regression_vec(11);
  double coef_Ferket_years_elapsed_age_baseline_squared = coef_Ferket_stroke_risk_time_regression_vec(12);
  double coef_Ferket_years_elapsed_male = coef_Ferket_stroke_risk_time_regression_vec(13);
  double coef_Ferket_years_elapsed_sbp_baseline = coef_Ferket_stroke_risk_time_regression_vec(14);
  double coef_Ferket_years_elapsed_sbp_baseline_squared = coef_Ferket_stroke_risk_time_regression_vec(15);
  double coef_Ferket_years_elapsed_cvd_history_baseline = coef_Ferket_stroke_risk_time_regression_vec(16);
  double coef_Ferket_years_elapsed_htnmed_baseline = coef_Ferket_stroke_risk_time_regression_vec(17);
  double coef_Ferket_years_elapsed_dm_baseline = coef_Ferket_stroke_risk_time_regression_vec(18);
  //double coef_Ferket_years_elapsed_cursmk_baseline = coef_Ferket_stroke_risk_time_regression_vec();
  double coef_Ferket_years_elapsed_race_black = coef_Ferket_stroke_risk_time_regression_vec(19);
  
  // Calculate log odds
  double log_odds = coefficients[0] + male * coefficients[1] + RIDAGEYR * coefficients[2] + 
    race * coefficients[3] + BPXSY1 * coefficients[4] + 
    SMQ040 * coefficients[5] + SBPtxt * coefficients[6] + 
    diabetes_rollup * coefficients[7] + vascular_history * coefficients[8];
  
  // Convert log odds to probability
  double probability = 1 / (1 + exp(-log_odds));
  
  // Convert probability to log prob
  double log_prob = log(probability);
  
  // adding Ferket into the probability
  log_prob = log_prob + coef_Ferket_years_elapsed * years_elapsed + 
    coef_Ferket_years_elapsed_age_baseline * years_elapsed * RIDAGEYR +
    coef_Ferket_years_elapsed_age_baseline_squared * years_elapsed * RIDAGEYR * RIDAGEYR +
    coef_Ferket_years_elapsed_male * years_elapsed * male +
    coef_Ferket_years_elapsed_sbp_baseline * years_elapsed * BPXSY1 +
    coef_Ferket_years_elapsed_sbp_baseline_squared * years_elapsed *BPXSY1 * BPXSY1 +
    coef_Ferket_years_elapsed_cvd_history_baseline * years_elapsed * vascular_history +
    coef_Ferket_years_elapsed_htnmed_baseline * years_elapsed * SBPtxt +
    coef_Ferket_years_elapsed_dm_baseline * years_elapsed * diabetes_rollup +
    //coef_Ferket_years_elapsed_cursmk_baseline * years_elapsed * SMQ040 +
    coef_Ferket_years_elapsed_race_black * years_elapsed * race;
  
  // Convert log prob back to prob
  probability = exp(log_prob);
  
  return probability;
  
}

// [[Rcpp::export]]
double predict_stroke_sten(double p_stroke_avg, int stenosis_cat,  // Adjust risk by stenosis severity
                           double or_stroke_cat1,
                           double or_stroke_cat2,
                           double or_stroke_cat3,
                           double hr_stroke_cat4){
  
  double p_stroke = p_stroke_avg;
  
  if(stenosis_cat == 1){
    p_stroke = OR_to_RR(or_stroke_cat1, p_stroke_avg)*p_stroke_avg;
  }else if(stenosis_cat == 2){
    p_stroke = OR_to_RR(or_stroke_cat2, p_stroke_avg)*p_stroke_avg;
  }else if(stenosis_cat == 3){
    p_stroke = OR_to_RR(or_stroke_cat3, p_stroke_avg)*p_stroke_avg;
  }else if(stenosis_cat == 4){
    p_stroke = RateToProb(hr_stroke_cat4*ProbToRate(p_stroke_avg,1), 1);
  }
  
  return p_stroke;
}

// === SCREENING AND STENOSIS MANAGEMENT ===

// [[Rcpp::export]]
List mod_stenosis_screening(NumericVector rn_once, int stenosis_cat, bool is_TIA,  // DUS + MRA cascade
                            double p_DUS_sensitivity_70, double p_DUS_specificity_70,
                            double p_MRA_sensitivity, double p_MRA_specificity,
                            double c_DUS, double c_MRA
){
  // CRN
  int r_DUS = 0;
  int r_MRA = 1;
  int r_MRA_TIA = 5;
  
  // output
  bool revasc = 0;
  double c_screening = 0;
  
  if (is_TIA){
    c_screening += c_MRA;
    if (stenosis_cat != 4){
      if ((stenosis_cat > 1 & rn_once[r_MRA_TIA] < p_MRA_sensitivity) || 
          (stenosis_cat <= 1 & rn_once[r_MRA_TIA] < 1 - p_MRA_specificity)) {
        revasc = 1;
      }
    }
  } else {
    c_screening += c_DUS;
    if (stenosis_cat != 4){
      if ((stenosis_cat > 1 & rn_once[r_DUS] < p_DUS_sensitivity_70) || 
          (stenosis_cat <= 1 & rn_once[r_DUS] < 1 - p_DUS_specificity_70)) {
        c_screening += c_MRA;
        if ((stenosis_cat > 1 & rn_once[r_MRA] < p_MRA_sensitivity) || 
            (stenosis_cat <= 1 & rn_once[r_MRA] < 1 - p_MRA_specificity)) {
          revasc = 1;
        }
      }
    }
  }
  
  return(List::create(
      Named("revasc") = revasc,
      Named("c_screening") = c_screening));
}


// [[Rcpp::export]]
int mod_stenosis_dev(double rn_stenosis_dev,  // Annual stenosis progression/regression 
                     int stenosis_cat, 
                      bool revasc, 
                      double p_restenosis_revasc, 
                      double p_stenosis_prog_cat0, 
                      double p_stenosis_prog,
                      double p_stenosis_prog_2cats, 
                      double p_stenosis_prog_3cats,
                      double p_stenosis_reg){
  if (stenosis_cat == 0) {
    if(revasc){
      //with revasc history
      if (rn_stenosis_dev < p_restenosis_revasc) {
        stenosis_cat++;
      }
    }else{
      //no revasc history
      if (rn_stenosis_dev < p_stenosis_prog_cat0) {
        stenosis_cat++;
      }
    }
    
  } else {
    if(stenosis_cat != 4){
      //double rn_stenosis_dev = R::unif_rand();
      if (rn_stenosis_dev < p_stenosis_prog * p_stenosis_prog_2cats * p_stenosis_prog_3cats){
        stenosis_cat = std::min(stenosis_cat + 3, 4);
      } else if (rn_stenosis_dev < p_stenosis_prog * p_stenosis_prog_2cats) {
        stenosis_cat = std::min(stenosis_cat + 2, 4);
      } else if (rn_stenosis_dev < p_stenosis_prog) {
        stenosis_cat = std::min(stenosis_cat + 1, 4);
      } else if (rn_stenosis_dev > 1 - p_stenosis_reg){
        stenosis_cat = std::max(stenosis_cat - 1, 0);
      }
    }
  }
  return(stenosis_cat);
}





// === MAIN SIMULATION ENGINE ===

// [[Rcpp::export]]
List mod_main_skipIS(int n_pop, bool screening_strategy,  // Main Markov cohort simulation
              NumericVector c_acuteIS_vec, NumericVector u_acuteIS_vec,
              NumericVector mRS_prop_acuteIS_vec,
              NumericVector age_vec, NumericVector male_vec, NumericVector race_vec,
              NumericVector SBP_vec, NumericVector smoking_vec, NumericVector SBPtxt_vec,
              NumericVector diabetes_vec, NumericVector vasc_hist_vec,
              NumericVector category_vec, 
              NumericVector coefficients,
              NumericVector coef_Ferket_stroke_risk_time_regression_vec,
              NumericVector pdie_m, NumericVector pdie_f,
              double disc_rate,
              double r_recur_postIS_1yr, double r_recur_postIS_post1yr, double p_recur_ICH,
              NumericVector hr_vec_recur_postIS, NumericVector hr_vec_mort_postIS,
              NumericVector p_vec_mRS_ICH,
              double r_IS_postICH, double r_ICH_postICH, double hr_mort_postICH,
              double c_ICH,
              NumericVector u_vec_mRS,
              NumericVector c_vec_mRS,
              double p_DUS_sensitivity_70, double p_DUS_specificity_70,
              double p_MRA_sensitivity, double p_MRA_specificity,
              double p_stenosis_prog_cat0,
              double p_stenosis_prog,
              double p_stenosis_prog_2cats,
              double p_stenosis_prog_3cats,
              double p_stenosis_reg,
              double p_revasc_complication,
              double p_stroke_revasc_complication,
              double p_death_revasc_complication,
              NumericVector p_complication_mRS_vec,
              double p_restenosis_revasc,
              double hr_stroke_after_revasc,
              NumericVector or_rr_stroke_vec,
              double c_DUS, double c_MRA, double c_revasc_asymp, double c_revasc_symp,
              NumericVector c_healthy_coef_vec, 
              double u_healthy, double u_revasc, double t_revasc,
              NumericVector coef_TIA_female,
              NumericVector coef_TIA_male, 
              double hr_stroke_after_TIA,
              double c_TIA, double u_TIA, double c_med_TIA,
              NumericVector randomnums
){
  //Output
  NumericVector LE_vec(n_pop);
  NumericVector LE_disc_vec(n_pop);
  NumericVector QALE_vec(n_pop);
  NumericVector QALE_disc_vec(n_pop);
  NumericVector cost_vec(n_pop);
  NumericVector cost_disc_vec(n_pop);
  NumericVector t_rstroke_vec(n_pop);
  NumericVector stroke_vec(n_pop);
  NumericVector true_positive_vec(n_pop);
  NumericVector false_negative_vec(n_pop);
  NumericVector false_positive_vec(n_pop);
  NumericVector true_negative_vec(n_pop);
  
  //life table
  double pdie[2][101];
  for(int i = 0; i < 100; i++){
    pdie[0][i] = pdie_f[i];
    pdie[1][i] = pdie_m[i];
  }
  
  double or_stroke_cat1 = or_rr_stroke_vec(0); 
  double or_stroke_cat2 = or_rr_stroke_vec(1);
  double or_stroke_cat3 = or_rr_stroke_vec(2);
  double hr_stroke_cat4 = or_rr_stroke_vec(3);
  
  //Common Random Number
  int age_min = 50;
  int tmax = 100 - age_min;
  // type of events: indice
  // initialize
  int r_complication = 2;
  int r_type_complication = 3;
  int r_mRS_complication = 4;
  
  int r_complication_TIA = 6;
  int r_type_complication_TIA = 7;
  int r_mRS_complication_TIA = 8;
  int length_events_once = 9;
  
  // looping
  int r_bg_mortality = 0;
  int r_ICH = 1;
  int r_IS = 2;
  int r_mRS = 3;
  int r_TIA = 4;
  int r_stenosis_dev = 5;
  int length_events_loop = 6;
  
  //looping over population
  for(int id = 0; id < n_pop; id++){
    //Read in individual characteristics
    int age_ini = age_vec(id); //initial age
    int age_mod = age_ini;
    bool male   = male_vec(id);
    int stenosis_cat = category_vec(id);
    bool race = race_vec(id);
    double SBP = SBP_vec(id);
    bool smoking = smoking_vec(id);
    bool SBPtxt = SBPtxt_vec(id);
    bool diabetes = diabetes_vec(id);
    bool vasc_hist = vasc_hist_vec(id);
    
    //Define tracker variables
    char state = 'H';
    double cycle = 0; // may be not integer
    double cycle_offset = 0;
    double p_stroke_updating;
    double p_stroke;
    int mRS = 0;
    bool revasc = 0;
    
    //Define randomnums
    NumericVector rn_once(length_events_once);
    for (int i = 0; i < length_events_once; i++){
      rn_once[i] = randomnums[id + length_events_once + i];
    }
    
    //Screening or not (depends on the strategies)
    if (screening_strategy) {
      // use a function
      List output_screening  = mod_stenosis_screening(rn_once, stenosis_cat, 0,
                                                      p_DUS_sensitivity_70, 
                                                      p_DUS_specificity_70,
                                                      p_MRA_sensitivity, 
                                                      p_MRA_specificity,
                                                      c_DUS, c_MRA);
      
      double c_screening = output_screening["c_screening"];
      revasc = output_screening["revasc"];
      cost_vec(id) += c_screening;
      cost_disc_vec(id) += c_screening /pow(1.0 + disc_rate, cycle);
      
      // tp, fp, tn, fn
      if (revasc){
        //positive
        if (stenosis_cat > 1) {
          true_positive_vec(id) = 1;
        } else {
          false_positive_vec(id) = 1;
        }
      } else {
        // negative
        if (stenosis_cat > 1 & stenosis_cat < 4){
          false_negative_vec(id) = 1;
        } else {
          true_negative_vec(id) = 1;
        }
      }
      
      if (revasc) {
        
        if(rn_once[r_complication] > p_revasc_complication){
          //no complication
          stenosis_cat = 0; //return to 0-49
          
          //revasc for asymp cost
          cost_vec(id) += c_revasc_asymp;
          cost_disc_vec(id) += c_revasc_asymp /pow(1.0 + disc_rate, cycle);
          
          cycle_offset = t_revasc;
          
          //Utility
          QALE_vec(id) += u_revasc*cycle_offset;
          QALE_disc_vec(id) += u_revasc*cycle_offset /pow(1.0 + disc_rate, cycle);
          
          // LE
          LE_vec(id) += cycle_offset;
          LE_disc_vec(id) += cycle_offset /pow(1.0+disc_rate, cycle);
          
        } else { 
          if (rn_once[r_type_complication] <= p_stroke_revasc_complication){
            // Complication == stroke, keep the stenosis_cat the same
            // the complication will finish this cycle, 
            // not considering the 2nd stroke right now.
            state = 'I';
            double rn_complication_mRS = rn_once[r_mRS_complication];
            if (rn_complication_mRS < p_complication_mRS_vec[0]) {
              mRS = 0;
            } else if (rn_complication_mRS < p_complication_mRS_vec[1]) {
              mRS = 1;
            } else if (rn_complication_mRS < p_complication_mRS_vec[2]) {
              mRS = 2;
            } else if (rn_complication_mRS < p_complication_mRS_vec[3]) {
              mRS = 3;
            } else if (rn_complication_mRS < p_complication_mRS_vec[4]) {
              mRS = 4;
            } else {
              mRS = 5;
            }
            
            //revasc for asymp cost
            cost_vec(id) += c_revasc_symp;
            cost_disc_vec(id) += c_revasc_symp /pow(1.0 + disc_rate, cycle);
            
            //Utility
            QALE_vec(id) += u_vec_mRS(mRS);
            QALE_disc_vec(id) += u_vec_mRS(mRS) /pow(1.0 + disc_rate, cycle);
            
            //LE
            LE_vec(id) += 1;
            LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
            
            age_mod++;
            cycle++;
            
          } else {
            //complication == death
            //revasc for symp cost
            cost_vec(id) += c_revasc_symp;
            cost_disc_vec(id) += c_revasc_symp /pow(1.0 + disc_rate, cycle);
            
            mRS = 6;
          }
        } // revasc finished
      }// screening finished
    }
    // else if not screening, nothing happened
    
    // loops of years
    while(mRS != 6 & age_mod < 100){
      
      if (state == 'H'){
        //1. Die from background mortality
        double p_die_bg = pdie[male][age_mod];
        
        if(randomnums[id + length_events_once + tmax * r_bg_mortality + (age_mod - age_min)] < p_die_bg){
          mRS = 6;
          break;
        }
        
        // update stroke risk every 10 years
        if (std::remainder(cycle, 10.0) < 1 & std::remainder(cycle, 10.0) >= 0) {
          //initialize or update stroke risk
          p_stroke_updating = ProbToProb(predict_stroke_update(
            coefficients, coef_Ferket_stroke_risk_time_regression_vec, cycle, male, age_ini,
            race, SBP, smoking, SBPtxt, diabetes, vasc_hist), 10);
          }
        
        p_stroke = predict_stroke_sten(p_stroke_updating,
                                                stenosis_cat,
                                                or_stroke_cat1,
                                                or_stroke_cat2,
                                                or_stroke_cat3,
                                                hr_stroke_cat4);
          
        if(revasc){
          p_stroke = p_stroke*hr_stroke_after_revasc;
        }
        
        
        //2. IS happening?
        if (randomnums[id + length_events_once + tmax * r_IS + (age_mod - age_min)] < p_stroke) {
          state = 'I';
          //Enter acute IS module
          //Run acute stroke module
          List output_acuteIS = mod_acuteIS_select(randomnums[id + length_events_once + tmax * r_mRS + (age_mod - age_min)],
                                                   age_mod, mRS,
                                                   c_acuteIS_vec, 
                                                   u_acuteIS_vec,
                                                   mRS_prop_acuteIS_vec);
          
          double c_acuteIS = output_acuteIS["c_acuteIS"];
          double u_acuteIS = output_acuteIS["u_acuteIS"];
          
          //mRS
          mRS = output_acuteIS["mRS"];

          //Acute IS cost
          cost_vec(id) += c_acuteIS;
          cost_disc_vec(id) += c_acuteIS /pow(1.0 + disc_rate, cycle);
          //Utility
          QALE_vec(id) += u_acuteIS*(1 - cycle_offset);
          QALE_disc_vec(id) += u_acuteIS*(1 - cycle_offset) /pow(1.0 + disc_rate, cycle);
          
          //LE
          if(mRS == 6){
            cycle_offset = 0;
            break;
          } else {
            LE_vec(id) += 1 - cycle_offset;
            LE_disc_vec(id) += (1 - cycle_offset) /pow(1.0+disc_rate, cycle);
            
            cycle_offset = 0;
            ++age_mod;
            ++cycle;
            continue;
          }
        }
        
        //3. TIA happening?
        // p_TIA needs to be defined
        double p_TIA;
        
        if (male == 0){
          p_TIA = exp(coef_TIA_female[0] + coef_TIA_female[1] * age_mod);
        } else {
          p_TIA = exp(coef_TIA_male[0] + coef_TIA_male[1] * age_mod);
        }
        
        if (randomnums[id + length_events_once + tmax * r_TIA + (age_mod - age_min)]< p_TIA){
          // TIA event
          state = 'T';
          // <- + TIA cost 
          cost_vec(id) += c_TIA;
          // <- ACAS screen
          List output_screening_TIA = mod_stenosis_screening(rn_once, stenosis_cat, 1,
                                                             p_DUS_sensitivity_70, 
                                                             p_DUS_specificity_70,
                                                             p_MRA_sensitivity, 
                                                             p_MRA_specificity,
                                                             c_DUS, c_MRA);
            
          double c_screening = output_screening_TIA["c_screening"];
          bool revasc_TIA = output_screening_TIA["revasc"];
          if(revasc_TIA){revasc = 1;}
          cost_vec(id) += c_screening;
          cost_disc_vec(id) += c_screening /pow(1.0 + disc_rate, cycle);
          
          if (revasc_TIA) {
            if(rn_once[r_complication_TIA] > p_revasc_complication){
              //no complication
              stenosis_cat = 0; //return to 0-49
              
              //revasc for symp cost
              cost_vec(id) += c_revasc_symp;
              cost_disc_vec(id) += c_revasc_symp /pow(1.0 + disc_rate, cycle);
                
              cycle_offset = t_revasc;
                
              //Utility
              QALE_vec(id) += u_revasc*cycle_offset;
              QALE_disc_vec(id) += u_revasc*cycle_offset /pow(1.0 + disc_rate, cycle);
                
              // LE
              LE_vec(id) += cycle_offset;
              LE_disc_vec(id) += cycle_offset /pow(1.0+disc_rate, cycle);
                
            } else { 
              if (rn_once[r_type_complication_TIA] <= p_stroke_revasc_complication){
                // Complication == stroke, keep the stenosis_cat the same
                // the complication will finish this cycle, 
                // not considering the 2nd stroke right now.
                  state = 'I';
                  double rn_complication_mRS = rn_once[r_mRS_complication_TIA];
                  if (rn_complication_mRS < p_complication_mRS_vec[0]) {
                    mRS = 0;
                  } else if (rn_complication_mRS < p_complication_mRS_vec[1]) {
                    mRS = 1;
                  } else if (rn_complication_mRS < p_complication_mRS_vec[2]) {
                    mRS = 2;
                  } else if (rn_complication_mRS < p_complication_mRS_vec[3]) {
                    mRS = 3;
                  } else if (rn_complication_mRS < p_complication_mRS_vec[4]) {
                    mRS = 4;
                  } else {
                    mRS = 5;
                  }
                  
                  //revasc for symp cost
                  cost_vec(id) += c_revasc_symp;
                  cost_disc_vec(id) += c_revasc_symp /pow(1.0 + disc_rate, cycle);
                  
                  //Utility
                  QALE_vec(id) += u_vec_mRS(mRS);
                  QALE_disc_vec(id) += u_vec_mRS(mRS) /pow(1.0 + disc_rate, cycle);
                  
                  //LE
                  LE_vec(id) += 1;
                  LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
                  
                  age_mod++;
                  cycle++;
                  
                } else {
                  //complication == death
                  //revasc for symp cost
                  cost_vec(id) += c_revasc_symp;
                  cost_disc_vec(id) += c_revasc_symp /pow(1.0 + disc_rate, cycle);
                  
                  mRS = 6;
                }
              } // complication
            }// revasc finished
          

        } // TIA finished
        
        //4. stenosis development(only being considered if state == H)
        if(state == 'H'){
          stenosis_cat = mod_stenosis_dev(randomnums[id + length_events_once + tmax * r_stenosis_dev + (age_mod - age_min)],
                                          stenosis_cat, revasc, 
                                          p_restenosis_revasc, 
                                          p_stenosis_prog_cat0, 
                                          p_stenosis_prog,
                                          p_stenosis_prog_2cats, 
                                          p_stenosis_prog_3cats,
                                          p_stenosis_reg);
        }
        
        double c_healthy = std::max(((c_healthy_coef_vec[2] - c_healthy_coef_vec[4] * c_healthy_coef_vec[0])*age_mod +
                                    (c_healthy_coef_vec[1] - c_healthy_coef_vec[3] * c_healthy_coef_vec[0])) / 
                                    (1 - c_healthy_coef_vec[3] - c_healthy_coef_vec[4]*age_mod), c_vec_mRS[0]);
        
        //Cycle reward
        cost_vec(id) += c_healthy;
        cost_disc_vec(id) += c_healthy /pow(1.0 + disc_rate, cycle);
        QALE_vec(id) += u_healthy * (1 - cycle_offset);
        QALE_disc_vec(id) += u_healthy * (1 - cycle_offset) /pow(1.0 + disc_rate, cycle);
        LE_vec(id) += 1 - cycle_offset;
        LE_disc_vec(id) += (1 - cycle_offset) /pow(1.0+disc_rate, cycle);
        cycle_offset = 0;
        
      } else if (state == 'T'){
        //1. Die from background mortality
        double p_die_bg = pdie[male][age_mod];
        if(randomnums[id + length_events_once + tmax * r_bg_mortality + (age_mod - age_min)] < p_die_bg){
          mRS = 6;
          break;
        }
        
        // stroke risk
        if (std::remainder(cycle, 10.0) < 1 & std::remainder(cycle, 10.0) >= 0) {
          //initialize or update stroke risk
          p_stroke_updating = ProbToProb(predict_stroke_update(
            coefficients, coef_Ferket_stroke_risk_time_regression_vec, cycle, male, age_ini,
            race, SBP, smoking, SBPtxt, diabetes, vasc_hist), 10);
        }
        
        p_stroke = predict_stroke_sten(p_stroke_updating,
                                       stenosis_cat,
                                       or_stroke_cat1,
                                       or_stroke_cat2,
                                       or_stroke_cat3,
                                       hr_stroke_cat4)*hr_stroke_after_TIA;
        
        if(revasc){
          p_stroke = p_stroke*hr_stroke_after_revasc;
        }
        
        
        //2. IS happening?
        if (randomnums[id + length_events_once + tmax * r_IS + (age_mod - age_min)] < p_stroke) {
          state = 'I';
          //Enter acute IS module
          //Run acute stroke module
          List output_acuteIS = mod_acuteIS_select(randomnums[id + length_events_once + tmax * r_mRS + (age_mod - age_min)], 
                                                   age_mod, mRS,
                                                   c_acuteIS_vec, 
                                                   u_acuteIS_vec,
                                                   mRS_prop_acuteIS_vec);
          
          double c_acuteIS = output_acuteIS["c_acuteIS"];
          double u_acuteIS = output_acuteIS["u_acuteIS"];
          
          //mRS
          mRS = output_acuteIS["mRS"];
          
          //Acute IS cost
          cost_vec(id) += c_acuteIS;
          cost_disc_vec(id) += c_acuteIS /pow(1.0 + disc_rate, cycle);
          //Utility
          QALE_vec(id) += u_acuteIS*(1 - cycle_offset);
          QALE_disc_vec(id) += u_acuteIS*(1 - cycle_offset) /pow(1.0 + disc_rate, cycle);
          
          //LE
          if(mRS == 6){
            cycle_offset = 0;
            break;
          } else {
            LE_vec(id) += 1 - cycle_offset;
            LE_disc_vec(id) += (1 - cycle_offset) /pow(1.0+disc_rate, cycle);
            
            cycle_offset = 0;
            ++age_mod;
            ++cycle;
            continue;
          
          }
        }
        
        //3. stenosis development(only being considered if state == H or T)
        if(state == 'T'){
          stenosis_cat = mod_stenosis_dev(randomnums[id + length_events_once + tmax * r_stenosis_dev + (age_mod - age_min)],
                                          stenosis_cat, revasc, 
                                          p_restenosis_revasc, 
                                          p_stenosis_prog_cat0, 
                                          p_stenosis_prog,
                                          p_stenosis_prog_2cats, 
                                          p_stenosis_prog_3cats,
                                          p_stenosis_reg);
        }
        
        double c_healthy = std::max(((c_healthy_coef_vec[2] - c_healthy_coef_vec[4] * c_healthy_coef_vec[0])*age_mod +
                                    (c_healthy_coef_vec[1] - c_healthy_coef_vec[3] * c_healthy_coef_vec[0])) / 
                                    (1 - c_healthy_coef_vec[3] - c_healthy_coef_vec[4]*age_mod), c_vec_mRS[0]);
        
        //Cycle reward
        cost_vec(id) += c_healthy + c_med_TIA;
        cost_disc_vec(id) += (c_healthy + c_med_TIA) /pow(1.0 + disc_rate, cycle);
        QALE_vec(id) += u_TIA * (1 - cycle_offset);
        QALE_disc_vec(id) += u_TIA * (1 - cycle_offset) /pow(1.0 + disc_rate, cycle);
        LE_vec(id) += 1 - cycle_offset;
        LE_disc_vec(id) += (1 - cycle_offset) /pow(1.0+disc_rate, cycle);
        cycle_offset = 0;
        
      } else if (state == 'I'){
        //Post IS
        
        //1. Die from background mortality
        double p_die_bg = RateToProb(ProbToRate(pdie[male][age_mod], 1) * hr_vec_mort_postIS(mRS) / hr_vec_mort_postIS(3), 1);
        if(randomnums[id + length_events_once + tmax * r_bg_mortality + (age_mod - age_min)] < p_die_bg){
          mRS = 6;
          break;
        }
        
        //2. Recurrent Stroke?
        
        //ICH
        if(randomnums[id + length_events_once + tmax * r_ICH + (age_mod - age_min)] < p_recur_ICH *
           RateToProb(r_recur_postIS_post1yr * hr_vec_recur_postIS(mRS) / hr_vec_recur_postIS(3), 1)){
          //2.1 ICH?
          state = 'B';
          mRS = map_OTT_mRS_ICH_CRN(mRS, p_vec_mRS_ICH, randomnums[id + length_events_once + tmax * r_mRS + (age_mod - age_min)]);
          if(t_rstroke_vec(id) == 0){
            t_rstroke_vec(id) = 0.5 + cycle;
          }
          
          //Acute ICH cost
          cost_vec(id) += c_ICH + 0.75*c_vec_mRS[mRS];
          cost_disc_vec(id) += (c_ICH + 0.75*c_vec_mRS[mRS]) /pow(1.0 + disc_rate, cycle);
          //Utility
          QALE_vec(id) += u_vec_mRS[mRS];
          QALE_disc_vec(id) += u_vec_mRS[mRS] /pow(1.0 + disc_rate, cycle);
          
          //LE
          if(mRS == 6){
            break;
          }else{
            LE_vec(id) += 1;
            LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
            
            ++age_mod;
            ++cycle;
            continue;
          }
        }
        
        //IS
        if(randomnums[id + length_events_once + tmax * r_IS + (age_mod - age_min)] < (1 - p_recur_ICH) *
           RateToProb(r_recur_postIS_post1yr * hr_vec_recur_postIS(mRS) / hr_vec_recur_postIS(3), 1)){
          //2.2 IS?
          if(t_rstroke_vec(id) == 0){
            t_rstroke_vec(id) = 0.5 + cycle;
          }
          //Enter acute IS module
          //Run acute stroke module
          List output_acuteIS = mod_acuteIS_select(randomnums[id + length_events_once + tmax * r_mRS + (age_mod - age_min)], 
                                                   age_mod, mRS,
                                                   c_acuteIS_vec, 
                                                   u_acuteIS_vec,
                                                   mRS_prop_acuteIS_vec);
          
          double c_acuteIS = output_acuteIS["c_acuteIS"];
          double u_acuteIS = output_acuteIS["u_acuteIS"];
          //mRS
          mRS = output_acuteIS["mRS"];
          
          //Acute IS cost
          cost_vec(id) += c_acuteIS;
          cost_disc_vec(id) += c_acuteIS /pow(1.0 + disc_rate, cycle);
          //Utility
          QALE_vec(id) += u_acuteIS;
          QALE_disc_vec(id) += u_acuteIS /pow(1.0 + disc_rate, cycle);
          
          //LE
          if(mRS == 6){
            break;
          }else{
            LE_vec(id) += 1;
            LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
            
            ++age_mod;
            ++cycle;
            continue;
          }
          
        }
        
        //Cycle reward
        cost_vec(id) += c_vec_mRS(mRS);
        cost_disc_vec(id) += c_vec_mRS(mRS) /pow(1.0 + disc_rate, cycle);
        QALE_vec(id) += u_vec_mRS(mRS);
        QALE_disc_vec(id) += u_vec_mRS(mRS) /pow(1.0 + disc_rate, cycle);
        LE_vec(id) += 1;
        LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
        
      } else if (state == 'B'){
        //Post both IS and ICH
        //Post IS
        
        //1. Die from background mortality
        double p_die_bg = std::max(RateToProb(ProbToRate(pdie[male][age_mod], 1) * hr_mort_postICH, 1),
                                   RateToProb(ProbToRate(pdie[male][age_mod], 1) * hr_vec_mort_postIS(mRS) / hr_vec_mort_postIS(3), 1));
        if(randomnums[id + length_events_once + tmax * r_bg_mortality + (age_mod - age_min)] < p_die_bg){
          mRS = 6;
          break;
        }
        
        //2. Recurrent Stroke?
        
        if(randomnums[id + length_events_once + tmax * r_ICH + (age_mod - age_min)] < RateToProb(r_ICH_postICH, 1)){
          //2.1 ICH?
          if(t_rstroke_vec(id) == 0){
            t_rstroke_vec(id) = 0.5 + cycle;
          }
          mRS = map_OTT_mRS_ICH_CRN(mRS, p_vec_mRS_ICH, randomnums[id + length_events_once + tmax * r_mRS + (age_mod - age_min)]);
          
          //Acute ICH cost
          cost_vec(id) += c_ICH + 0.75*c_vec_mRS[mRS];
          cost_disc_vec(id) += (c_ICH + 0.75*c_vec_mRS[mRS]) /pow(1.0 + disc_rate, cycle);
          //Utility
          QALE_vec(id) += u_vec_mRS[mRS];
          QALE_disc_vec(id) += u_vec_mRS[mRS] /pow(1.0 + disc_rate, cycle);
          
          //LE
          if(mRS == 6){
            break;
          }else{
            LE_vec(id) += 1;
            LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
            
            ++age_mod;
            ++cycle;
            continue;
          }
        }
        
        if(randomnums[id + length_events_once + tmax * r_IS + (age_mod - age_min)] < std::max(RateToProb(r_IS_postICH, 1),
                        RateToProb(r_recur_postIS_post1yr * hr_vec_recur_postIS(mRS) / hr_vec_recur_postIS(3), 1)*(1-p_recur_ICH))){
          //2.2 IS?
          if(t_rstroke_vec(id) == 0){
            t_rstroke_vec(id) = 0.5 + cycle;
          }
          //Enter acute IS module
          //Run acute stroke module
          List output_acuteIS = mod_acuteIS_select(randomnums[id + length_events_once + tmax * r_mRS + (age_mod - age_min)],
                                                   age_mod, mRS,
                                                   c_acuteIS_vec, 
                                                   u_acuteIS_vec,
                                                   mRS_prop_acuteIS_vec);
          
          double c_acuteIS = output_acuteIS["c_acuteIS"];
          double u_acuteIS = output_acuteIS["u_acuteIS"];
          //mRS
          mRS = output_acuteIS["mRS"];
          
          //Acute IS cost
          cost_vec(id) += c_acuteIS;
          cost_disc_vec(id) += c_acuteIS /pow(1.0 + disc_rate, cycle);
          //Utility
          QALE_vec(id) += u_acuteIS;
          QALE_disc_vec(id) += u_acuteIS /pow(1.0 + disc_rate, cycle);
          
          //LE
          if(mRS == 6){
            break;
          }else{
            LE_vec(id) += 1;
            LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
            
            ++age_mod;
            ++cycle;
            continue;
          }
        }
        
        //Cycle reward
        cost_vec(id) += c_vec_mRS(mRS);
        cost_disc_vec(id) += c_vec_mRS(mRS) /pow(1.0 + disc_rate, cycle);
        QALE_vec(id) += u_vec_mRS(mRS);
        QALE_disc_vec(id) += u_vec_mRS(mRS) /pow(1.0 + disc_rate, cycle);
        LE_vec(id) += 1;
        LE_disc_vec(id) += 1 /pow(1.0+disc_rate, cycle);
        
      } else {
        Rcerr << "Error, bad input, quitting\n";
      }
      
      ++age_mod;
      ++cycle;

    }
    
    if (state != 'H') {
      stroke_vec(id) = 1;
    }
  }
  
  return(List::create(
      Named("LE")        = LE_vec,
      Named("LE_disc")   = LE_disc_vec,
      Named("QALE")      = QALE_vec,
      Named("QALE_disc") = QALE_disc_vec,
      Named("cost")      = cost_vec,
      Named("cost_disc") = cost_disc_vec,
      Named("t_rstroke") = t_rstroke_vec,
      Named("stroke") = stroke_vec,
      Named("true_positive") = true_positive_vec,
      Named("false_positive") = false_positive_vec,
      Named("false_negative") = false_negative_vec,
      Named("true_negative") = true_negative_vec
  )
  );
}

// === HELPER FUNCTIONS ===

// [[Rcpp::export]]
double compute_c_healthy(NumericVector c_healthy_coef_vec, double age_mod) {  // Age-adjusted costs
  if(c_healthy_coef_vec.length() != 5) {
    stop("c_healthy_coef_vec should be of length 5");
  }
  
  double denominator = 1 - c_healthy_coef_vec[3] - c_healthy_coef_vec[4]*age_mod;
  
  if(denominator == 0) {
    stop("Division by zero error: denominator is zero");
  }
  
  double c_healthy = ((c_healthy_coef_vec[2] - c_healthy_coef_vec[4] * c_healthy_coef_vec[0])*age_mod +
                      (c_healthy_coef_vec[1] - c_healthy_coef_vec[3] * c_healthy_coef_vec[0])) / 
                      denominator;
  
  return c_healthy;
}