//#######################################################################
// A C++ program to provide examples of use of methods
// in the cpp_deSolve class to numerically
// solve an SIR model and compare the results to influenza
// data from the Midwest 2007-2008 B influenza season
// The data are contained in the FluData class, and were
// filled by the R script midwest_flu.R using a file from the
// CDC website.
//
// Author: Sherry Towers
//         smtowers@asu.edu
// Created: Mar 22nd, 2013
//
// Copyright Sherry Towers, 2013
//
// This program is not guaranteed to be free of bugs and/or errors.
//
// This program can be freely used and shared, as long as the author information
// and copyright in the header remains intact.
//#######################################################################
#include <algorithm>
#include <cstdlib>
#include <vector>
#include <iostream>
#include "UsefulUtils.h"
#include "FluData.h"
#include "cpp_deSolve.h"

//################################################################################
//################################################################################
using namespace std;
int main(int argc, char **argv){
  //##############################################################################
  // set the random seed
  // (either use process ID, or an integer argument to the program)
  //##############################################################################
  int myint = unsigned(getpid());
  //int myint = atoi(argv[1]);
  srand(unsigned(time(0))+myint);

  //##############################################################################
  // the class FluData contains the midwest influenza B data from flu
  // season 2007-08, as obtained from the CDC
  //##############################################################################
  UsefulUtils myutils;
  FluData myflu;
  vector<double> mydata;
  vector<double> time_lo;
  vector<double> time_hi;
  myflu.GetData(mydata
               ,time_lo
               ,time_hi);

  //##############################################################################
  // We will fit an SIR model
  // The model parameters for which we need to fit are the transmission, beta
  // and the time of introduction of the virus to the population, t0
  // Assume that we already know the recovery rate, gamma, from observational
  // studies of sick people
  //##############################################################################
  double gamma = 1.0/3.0;  // recovery rate for influenza

  for (int iter=0;iter<1000;iter++){
    //############################################################################
    // Randomly sample the parameters
    // Note that the time of introduction, t0, cannot be after the minimum
    // time in the data.
    // For an SIR or SEIR model, it makes more intuitive sense to sample
    // R0, and then calculate beta from that.
    //############################################################################
    double R0   = myutils.RandomUniform(1.0,2.0);
    double beta = R0*gamma;
    double t0   = double(int(myutils.RandomUniform(0.0,myutils.Min(time_lo))));

    //############################################################################
    // Let's solve the model for times t=t0 to the maximum
    // time in the data (plus a bit)  in steps of 1/100th of a day
    // Note that you need to make the time step quite small in order for this 
    // method to work!
    //############################################################################
  
    double mytime = t0;
    double delta_t = 0.01;
    vector<double> vtime;
    while (mytime<=((myutils.Max(time_hi)+1.0)+delta_t)){
       vtime.push_back(mytime);
       mytime = mytime + delta_t;
    }
  
    //############################################################################
    // fill vectors with the parameters and the initial values
    //############################################################################
    vector<double> Parameters;
    Parameters.push_back(beta);
    Parameters.push_back(gamma);
  
    double npop = 52000000.0; //this is approximately the population of IL IN MI MN OH WI (CDC region 5);
    double I_0 = 1.0;
    double S_0 = npop-I_0;
    double R_0 = 0.0;
    vector<double> InitialValues;
    InitialValues.push_back(S_0);
    InitialValues.push_back(I_0);
    InitialValues.push_back(R_0);
    
    //############################################################################
    // ModelType==1 is an SIR model
    //############################################################################
    int ModelType = 1;
    cpp_deSolve mysir(InitialValues
                     ,Parameters
                     ,ModelType);
  
    //############################################################################
    // now solve the model using the Runge Kutta method
    //############################################################################
    vector<vector<double> > ModelEstimates;
    vector<vector<double> > ModelEstimatesTranspose;
    mysir.SimulateModel(vtime
                       ,ModelEstimates 
                       ,ModelEstimatesTranspose
                       );
  
    //############################################################################
    // I've written the code in the cpp_deSolve SimulateModel method
    // such that when ModelType==1 then ModelEstimatesTranspose[3] is equal
    // to number of new infections per unit time in an SIR model 
    // beta*S*I/N
    //############################################################################
    vector<double> vnewI = ModelEstimatesTranspose[3];

    //############################################################################
    // now accumulate the vnewI vector in the same time steps as the data
    // vector
    //############################################################################
    vector<double> accumulated_model = myutils.AccumulateModel(mydata
                                                              ,time_hi
                                                              ,vnewI  
                                                              ,vtime);
    double mypears =  myutils.PearsonChiSquared(mydata
                                               ,accumulated_model);
    
    cout << t0 << " "
         << R0 << " "
         << beta << " "
         << gamma << " " 
         << mydata.size() << " "
         << vnewI.size() << " "
         << mypears << " " 
         << endl;
   
  } // end loop over iterations
 
  return 0;

} // end program