#include <admodel.h>

  extern "C"  {
    void ad_boundf(int i);
  }
#include <epiglm2.htp>

model_data::model_data(int argc,char * argv[]) : ad_comm(argc,argv)
{
  cgr_phase.allocate("cgr_phase");
  ctr_phase.allocate("ctr_phase");
  n.allocate("n");
  pgr.allocate("pgr");
  ptr.allocate("ptr");
  sumh.allocate(1,n,"sumh");
  tau.allocate(1,n,"tau");
  k.allocate(1,n,"k");
  maxk.allocate(1,n,"maxk");
  Zgr.allocate(1,n,1,pgr,"Zgr");
  Ztr.allocate(1,n,1,ptr,"Ztr");
  Y.allocate(1,n,1,k,"Y");
}

model_parameters::model_parameters(int sz,int argc,char * argv[]) : 
 model_data(argc,argv) , function_minimizer(sz)
{
  initializationfunction();
  ctr.allocate(1,ptr,-10,10,ctr_phase,"ctr");
  cgr.allocate(1,pgr,-10,10,cgr_phase,"cgr");
  survivorvec.allocate(1,n,"survivorvec");
  #ifndef NO_AD_INITIALIZE
    survivorvec.initialize();
  #endif
  betavec.allocate(1,n,"betavec");
  #ifndef NO_AD_INITIALIZE
    betavec.initialize();
  #endif
  extras.allocate(1,n,"extras");
  #ifndef NO_AD_INITIALIZE
    extras.initialize();
  #endif
  f.allocate("f");
}

void model_parameters::userfunction(void)
{
  f = 0;
  extras = mfexp(Zgr * cgr);
  betavec = mfexp(Ztr * ctr);
  for (int i=1; i <=n; i++)
  {
    survivorvec(i) = extras(i) + maxk(i) - k(i);
  }
  // negative log-likelihood
  for (int i=1; i<=n; i++)
  {
    f -= gammln(survivorvec(i) + k(i) + 1) - gammln(survivorvec(i) + 1);
    if (k(i) >= 1)
    {
      f -= -betavec(i) * sumh(i);
      for (int j=1; j <= k(i); j++)
      {
        f -= log(betavec(i) * Y(i,j));
      }
    }
    f -= survivorvec(i) * -betavec(i) * tau(i);
  }
}

void model_parameters::preliminary_calculations(void){
  admaster_slave_variable_interface(*this);
}

model_data::~model_data()
{}

model_parameters::~model_parameters()
{}

void model_parameters::report(void){}

void model_parameters::final_calcs(void){}

void model_parameters::set_runtime(void){}

#ifdef _BORLANDC_
  extern unsigned _stklen=10000U;
#endif


#ifdef __ZTC__
  extern unsigned int _stack=10000U;
#endif

  long int arrmblsize=0;

int main(int argc,char * argv[])
{
    ad_set_new_handler();
  ad_exit=&ad_boundf;
    gradient_structure::set_NO_DERIVATIVES();
    gradient_structure::set_YES_SAVE_VARIABLES_VALUES();
  #if defined(__GNUDOS__) || defined(DOS386) || defined(__DPMI32__)  || \
     defined(__MSVC32__)
      if (!arrmblsize) arrmblsize=150000;
  #else
      if (!arrmblsize) arrmblsize=25000;
  #endif
    model_parameters mp(arrmblsize,argc,argv);
    mp.iprint=10;
    mp.preliminary_calculations();
    mp.computations(argc,argv);
    return 0;
}

extern "C"  {
  void ad_boundf(int i)
  {
    /* so we can stop here */
    exit(i);
  }
}
