//#######################################################################
// A C++ program to provide examples of use of a class to numerically
// solve the Lotka-volterra model for predator/prey systems
// and fit it to the Canadian Hare/Lynx data set
//
// Author: Sherry Towers
//         smtowers@asu.edu
// Created: Feb 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 <cmath>
#include <iostream>
#include "UsefulUtils.h"
#include "PredatorPrey.h"
#include "HareLynxData.h"

//################################################################################
// The first argument to the program is used to help set the random seed
// The reason we do it this way is because when we submit multiple jobs
// simultaneously to batch, they will often start running at the same time stamp,
// thus it is prudent to also pass the job ID number to get a different
// random seed for each run of the program
//################################################################################
using namespace std;
int main(int argc, char **argv){
  srand(unsigned(time(0))+atoi(argv[1]));

  int nparam = 6;    // number of parameters we are fitting for
  int niter = 10000; // the number of iterations

  UsefulUtils myutil;

  HareLynxData mydata;
  vector<double> vtime;
  vector<double> vhare;
  vector<double> vlynx;
  mydata.GetData(vtime,vhare,vlynx);

  double mintime = myutil.Min(vtime);
  double maxtime = myutil.Max(vtime);

  double time = mintime-2.0;  // the minimum time, minus two years ( the data are binned in two year increments)
  double delta_t = 0.010;
  vector<double> vtimeb;
  while (time<=maxtime){
     vtimeb.push_back(time);
     time = time + delta_t;
  }
  //##############################################################################
  // do niter iterations: randomly sample a hypothesis for the initial conditions
  // and parameter values, calculate the predicted Hare and Lynx populations,
  // then compare the model estimates to the data
  // The values from the statistical fit are
  // double prey_0 = 20000;
  // double predator_0 = 32000;
  // double growth_rate_prey = 2.252;
  // double predation_encounter_rate = 5.006e-5;
  // double death_rate_predators = 0.12216;
  // double predation_efficiency = 0.242;
  //
  // The initial values are approximate, and found by initial exploratory 
  // studies in R in hare_lynx_fit_step1.R
  //##############################################################################
  vector<double> vparam_min;
  vector<double> vparam_max;
  vparam_min.push_back( 1.0e3 ); 
  vparam_min.push_back( 500.0 ); 
  vparam_min.push_back( 0.2 ); 
  vparam_min.push_back( 1e-07 ); 
  vparam_min.push_back( 0.2 ); 
  vparam_min.push_back( 0.001 ); 

  vparam_max.push_back( 5.0e6 ); 
  vparam_max.push_back( 1.0e6 ); 
  vparam_max.push_back( 1.5 ); 
  vparam_max.push_back( 1.0e-4 ); 
  vparam_max.push_back( 2.0 ); 
  vparam_max.push_back( 2.0 ); 

  for (int iter=0;iter<niter;iter++){
     //#######################################################################
     // uniformly randomly sample values for the model parameters
     //#######################################################################
     vector<double> vparam;
     for (int i=0;i<vparam_min.size();i++){
        vparam.push_back(myutil.RandomUniform(vparam_min[i],vparam_max[i]));
     }

     //#######################################################################
     // solve the differential equations
     //#######################################################################
     double prey_0 = vparam[0];
     double predator_0 = vparam[1];
     double growth_rate_prey = vparam[2];
     double predation_encounter_rate = vparam[3];
     double death_rate_predators = vparam[4];
     double predation_efficiency = vparam[5];

     PredatorPrey mylotka(prey_0
                         ,predator_0
                         ,growth_rate_prey
                         ,predation_encounter_rate
                         ,death_rate_predators
                         ,predation_efficiency);

     vector<double> vprey;
     vector<double> vpredator;
     mylotka.SimulatePredatorPrey(vtimeb,vprey,vpredator);
     //#######################################################################
     // now accumulate the model estimates into the same time bins
     // as the data, and then calculate the Pearson chi-squared goodness
     // of fit statistic
     //#######################################################################
     vector<double> vpreyb = myutil.AccumulateModel(vhare,vtime,vprey,vtimeb);
     vector<double> vpredatorb = myutil.AccumulateModel(vlynx,vtime,vpredator,vtimeb);

     double chisq_hare = myutil.PearsonChiSquared(vhare,vpreyb);
     double chisq_lynx = myutil.PearsonChiSquared(vlynx,vpredatorb);
     double total_chisq = chisq_hare + chisq_lynx;

     //#######################################################################
     // print out the model parameters
     // and goodness-of-fit to a file 
     //#######################################################################

     if (total_chisq<5.0e6){
        cout << iter << " ";
        for (int i=0;i<nparam;i++){
           cout << vparam[i] << " ";
        }
        cout << total_chisq     << " "   
             << endl;
     }
  }

  return 0;

} // end program