#include /include/license.stan

functions {

  #include /include/linear_interpolation.stan
  #include /include/loglogistic.stan

  vector TKTD_varIT( real t,
                     vector y,
                     vector theta,
                     vector expo_x,
                     vector expo_y) {
    // - parameters
    real kd = theta[1];
    // - new variables
    vector[1] dy_dt; //
    // interpolation of exposure
    real expo_t = interpolate(t, expo_x, expo_y) ;
    // - model
    dy_dt[1] =  kd * ( expo_t - y[1]);
    return(dy_dt);
  }
}

data {
   #include /include/data_guts.stan
   real alpha_meanlog10;
   real alpha_sdlog10;
   real beta_minlog10;
   real beta_maxlog10;
}
parameters {
    real kd_log10, hb_log10, alpha_log10, beta_log10 ;
}
transformed parameters{
    // Priors
    real kd, hb, alpha, beta ;
    // ODE
    vector<lower=0>[1] y0;
    // remove 1e-12 to start before min time !
    real t0 = 0.0 - 1e-12;
    vector[1] theta ;
    array[n_N] vector[1] y_hat;
    // model GUTS
    vector<lower=0, upper=1>[n_N] Psurv_hat;
    vector<lower=0, upper=1>[n_N] Conditional_Psurv_hat;
    
    kd = 10^kd_log10 ;
    if(HB_FIXED == 0){
        hb = 10^hb_log10 ;
    } else{
        hb = hb_value;
    }
    alpha = 10^alpha_log10 ;
    beta  = 10^beta_log10 ;
    
    theta[1] = kd;

    for (g in 1:n_group){
        array[GROUP_N[3,g]] real ts ;
        vector[GROUP_X[3,g]] exposure_t ;
        vector[GROUP_X[3,g]] exposure_v ;
        // add 1e-12 to be > 0!
        ts = time_N[GROUP_N[1,g]:GROUP_N[2,g]];
        exposure_t = to_vector(time_X[GROUP_X[1,g]:GROUP_X[2,g]]);
        exposure_v = to_vector(conc[GROUP_X[1,g]:GROUP_X[2,g]]);
        y0 = [0]';
        
        y_hat[GROUP_N[1,g]:GROUP_N[2,g]]  = ode_rk45(
          TKTD_varIT, y0, t0, ts, theta, 
          exposure_t, exposure_v
        );
          
        for(i in GROUP_N[1,g]:GROUP_N[2,g]){
          // Psurv_hat[i] = exp(- hb * time_N[i]) * (1-exp(loglogistic_lcdf(max(y_hat[GROUP_N[1,g]:i, 1]) | alpha, beta)));
          Psurv_hat[i] = exp(- hb * time_N[i]) * (1-exp(lognormal_lcdf(max(y_hat[GROUP_N[1,g]:i, 1]) | alpha, beta)));
          Conditional_Psurv_hat[i] = i == GROUP_N[1,g] ? Psurv_hat[i] : Psurv_hat[i] / Psurv_hat[i-1] ;
        }
    }
}
model {
    // PRIORS
    target += normal_lpdf(kd_log10 | kd_meanlog10, kd_sdlog10);
    target += normal_lpdf(hb_log10 | hb_meanlog10, hb_sdlog10);
    target += normal_lpdf(alpha_log10 | alpha_meanlog10, alpha_sdlog10);
    // USING LOGLOGISTIC, MEAN IS ALPHA, VARIANCE IS SIMILAR TO HB!
    target += normal_lpdf(beta_log10 | alpha_meanlog10 , hb_sdlog10 );
    // target += uniform_lpdf(beta_log10 | beta_minlog10 , beta_maxlog10 );
    for (g in 1:n_group){
      target += binomial_lpmf(
          Nsurv[GROUP_N[1,g]:GROUP_N[2,g]] | 
          Nprec[GROUP_N[1,g]:GROUP_N[2,g]], 
          Conditional_Psurv_hat[GROUP_N[1,g]:GROUP_N[2,g]]
        );
    }
}
generated quantities {
  #include /include/generated_quantities.stan
}
