DATA_SECTION

  init_int nobs               // # of observations
  init_int nblock             // # of blocks
  init_vector killed(1,nobs)  // # killed per trial
  init_vector size(1,nobs)    // size of individuals
  init_vector initial(1,nobs) // starting density (# individuals)
  init_matrix Z(1,nobs,1,nblock) // random-effects model matrix


PARAMETER_SECTION

  init_number c
  init_number d
  init_bounded_number h(0.0,1.0,1) 
  init_number g
  init_bounded_number sigma_c(0.00001,1.0,1)
  random_effects_vector u(1,nblock,2)
  vector prob(1,nobs)
  vector cvec(1,nobs)
  objective_function_value f

PROCEDURE_SECTION
  // dvariable fpen=0.0;
  // compute vector of c values
  cvec = c + sigma_c*(Z*u);
  // power-Ricker
  prob = 1/(1/(elem_prod(cvec,elem_prod(pow(size/d,g),exp(-size/d))))+h*initial);
  // penalties: constrain 0.001 <= prob <= 0.999
  //  prob = posfun(prob,0.001,fpen);
  // f += 1000*fpen;
  // prob = 1-posfun(1-prob,0.001,fpen);
  // f += 1000*fpen;
  // binomial negative log-likelihood: include norm constant
  f -= sum( log_comb(initial,killed)+
            elem_prod(killed,log(prob))+
            elem_prod(initial-killed,log(1-prob)));


GLOBALS_SECTION

  #include <admodel.h>
  #include <df1b2fun.h>
  #include <adrndeff.h>

  df1b2vector pow(const df1b2vector& v,const df1b2variable & x)
  {
    int mmin=v.indexmin();
    int mmax=v.indexmax();
    df1b2vector tmp(mmin,mmax);
    for (int i=mmin;i<=mmax;i++)
    {
      tmp(i)=pow(v(i),x);
    }
    return tmp;
  }