[ADMB Users] semi-implicit differential equation codes

[Saang-Yoon]

Dear all.
Following Dr. David Fournier's advice (e.g., ADMB manual version
10.0), I use a semi-implicit differential equation for the simple
logistic surplus production model; see my TPL codes below.  I
successfully did complie-link, and created the resultant executable
file (.exe). However, estimates of fishing mortality rates in years
are all zero, and I suspect something wrong in my TPL codes. Even
without DAT and PIN files, some of you may find what I have to do to
avoid the problem.  I could send you DAT and PIN files; here I don't
see the attachment function. I don't have data on fishing effort, and
my codes below are a little different from those in the ADMB manual.
My data are
time series only of commercial catch (or yield) and survey index.  I
would greatly appreciate your help.  Saang-Yoon



DATA_SECTION
  int i;  //for loop;
  int j;  //for loop;

  init_vector years(1,3);
  int nyrs;
  int ncol;  //data column
  !!nyrs=years(2)-years(1)+1;
  !!ncol=size_count(years);

  init_matrix yieldsurvey(1,nyrs,1,ncol);
  vector yrs(1,nyrs);
  vector Yt(1,nyrs);
  vector It(1,nyrs);

  !!yrs=trans(yieldsurvey)(1);
  !!Yt=trans(yieldsurvey)(2);  //yield;
  !!It=trans(yieldsurvey)(3);

  number minK;
  number maxK;
  !!minK=max(Yt)-max(Yt)*15/100;
  !!maxK=minK+minK*3000/100;

  number minY;
  number maxY;
  number meanY;
  !!minY=min(Yt);
  !!maxY=max(Yt);
  !!meanY=mean(Yt);

  !!cout<<"yrs: "<<yrs<<endl;
  !!cout<<"Yt: "<<Yt<<endl;
  !!cout<<"It: "<<It<<endl;
  !!cout<<"K range: "<<minK<<" "<<maxK<<endl;


PARAMETER_SECTION
  init_bounded_number q(0.0,1.0,1);
  init_bounded_number r(0.0,5,2);
  init_bounded_number K(minK,maxK,3);   //aka B_infinitive
  init_bounded_number b1(0.0,1.0,4);     //B1=b1*K;
  init_bounded_dev_vector F_devs(1,nyrs,-5.0,5.0,5);

  number B1;
  vector Bt(1,nyrs+1);
  vector Ft(1,nyrs);

  vector predYt(1,nyrs); //expected yield in year t
  vector predIt(1,nyrs);

  vector sig2It(1,nyrs);

  vector meanlogIt(1,nyrs);
  number varlogI;

  number varlogY;

  number Bpen;
  number bio_tmp; //biomass temporary
  number y_tmp; //yield temporary

  number Bmsy;
  number msy; //msy;
  number Fmsy;

  objective_function_value f;

PROCEDURE_SECTION
  B1=b1*K;

  f=0.0;
  Bpen=0;
  //cout<<"f1: "<<f<<endl;

  //Project B_t
  Bt(1)=B1;
  Bt(1)=posfun(Bt(1),Yt(1),Bpen); //keeps it above

  //Ft
  for(i=1;i<=nyrs;i++)
    Ft(i)=Ft(i)*exp(F_devs(i));

  //four time steps per year desired
  int nsteps=4;
  double delta=1.0/nsteps;
  //integrate the logistic dynamics over n time steps per year
  for(i=1;i<=nyrs;i++)   {
     bio_tmp=1.e-20+Bt(i);
	 y_tmp=0.0;  //yield
	 for(j=1;j<=nsteps;j++) {
	    bio_tmp=bio_tmp*(1.0+r*delta)/(1.0+(r*bio_tmp/K+Ft(i))*delta);
	    y_tmp=y_tmp+Ft(i)*bio_tmp*delta;
     };
     predYt(i)=y_tmp; //expected yield;
     Bt(i+1)=bio_tmp;
     Bt(i+1)=posfun(Bt(i+1),meanY,Bpen);  //keeps it above
  };
  f+=Bpen;

  //
  //survey as ROUTINE;
  for(i=1;i<=nyrs;i++)
	 predIt(i)=q*Bt(i);
  varlogI=norm2(log(It)-log(1.e-10+predIt))/nyrs;
  f+=0.5*nyrs*log(varlogI)+norm2(log(It)-log(1.e-10+predIt))/
(2.0*varlogI);
  cout<<"f3: "<<f<<endl;

  //
  //commerical yield as ROUTINE;
  varlogY=norm2(log(Yt)-log(1.e-10+predYt))/nyrs;
  f+=0.5*nyrs*log(varlogY)+norm2(log(Yt)-log(1.e-10+predYt))/
(2.0*varlogY);
  cout<<"f4: "<<f<<endl;


  //management references
  Bmsy=K/2;
  msy=r*K/4;  //MSY;
  Fmsy=msy/Bmsy;


REPORT_SECTION
  report<<"yr Yt It Bt predIt varlogI Ft Ft/Fmsy (Bt-Yt) f Bpen
max.grad"<<endl;
  for(i=1;i<=nyrs;i++)   {
    report<<yrs(i)<<" "<<Yt(i)<<" "<<It(i)<<" "<<Bt(i)<<"
"<<predIt(i)<<" ";
    report<<varlogI<<" "<<Ft(i)<<" "<<Ft(i)/Fmsy<<" ";
    report<<Bt(i)-Yt(i)<<" "<<f<<" "<<Bpen<<"
"<<objective_function_value::gmax<<endl;
  };


GLOBALS_SECTION
  #include <admodel.h>
  #include <math.h>
  #include <stdio.h>
  #include <stddef.h>
  #include <stdlib.h>