baranov.cpp
Go to the documentation of this file.
```00001 #include "statsLib.h"
00002
00051 double get_ft(const double& ct, const double& m, const dvector& va, const dvector& ba)
00052 {
00053   double ft;
00054   //initial guess for ft
00055   ft=ct/(va*(ba*exp(-m/2.)));
00056
00057   for(int i=1;i<=50;i++)
00058   {
00059     dvector f = ft*va;
00060     dvector z = m+f;
00061     dvector s = exp(-z);
00062     dvector o = (1.-s);
00063
00064     dvector t1 = elem_div(f,z);
00065     dvector t2 = elem_prod(t1,o);
00066     dvector t3 = elem_prod(o,ba);
00067     //predicted catch
00068     double pct = t2*ba;
00069
00070     //derivative of catch wrt ft
00071     double dct = sum(
00072       elem_div(t3,z)
00073       - elem_div(elem_prod(f,t3),square(z))
00074       + elem_prod(elem_prod(t1,s),ba));
00075
00076     ft -= (pct-ct)/dct;  //newton step
00077     //if(fabs(pct-ct)<1.e-9) break; //do not use for dvariables
00078   }
00079
00080   return(ft);
00081 }
00082
00094 dvector get_ft( dvector& ct,const double& m, const dmatrix& V,const dvector& ba )
00095 {
00096   /*  ARGUMENTS:
00097      ct is the observed catch for each fleet
00098      m is the natural mortality rate
00099      va is the vulnerability row fleet, col age va(1,ngear,1,nage)
00100      ba is the start of year biomass at age
00101   */
00102
00103   int i,a,A;
00104   double minsurv = 0.05;
00105   int ng=size_count(ct);  //number of gears
00106   a=ba.indexmin(); A=ba.indexmax();
00107   dvector ft(1,ng); ft=0;
00108   dvector ctmp(1,ng);
00109   dvector ct_hat(1,ng); //predicted catch
00110   dvector dct_hat(1,ng);  //derivative
00111   dvector zage(a,A);
00112   dvector sage(a,A);
00113   dvector ominus(a,A);
00114   dmatrix F(1,ng,a,A);
00115
00116
00117   //initial guess for ft=ct/(0.98 Bt);
00118   ctmp=ct;
00119
00120   for(i=1;i<=ng;i++)
00121   {
00122     ft(i) = ctmp(i)/(0.98*ba*V(i)*exp(-0.5*m));
00123     if(1.-ft(i)<minsurv)
00124     {
00125       ft(i)=1.-minsurv;
00126       ctmp=ft(i)*ba*V(i)*exp(-0.5*m);
00127     }
00128   }
00129   ct=ctmp;  //don't do this for the differentiable version.
00130
00131   //now solve baranov catch equation iteratively.
00132   for(int iter=1; iter<=17; iter++)
00133   {
00134     for(i=1;i<=ng;i++)F(i)=ft(i)*V(i);
00135     zage=m+colsum(F);
00136     sage=mfexp(-zage);
00137     ominus=(1.-sage);
00138
00139     for(i=1;i<=ng;i++)
00140     {
00141       dvector t1 = elem_div(F(i),zage);
00142       dvector t2 = elem_prod(t1,ominus);
00143       dvector t3 = elem_prod(ominus,ba);
00144
00145       ct_hat(i) = t2*ba;
00146
00147       dct_hat(i) = sum(
00148             elem_div(t3,zage)
00149             - elem_div(elem_prod(F(i),t3),square(zage))
00150             + elem_prod(elem_prod(t1,sage),ba));
00151
00152         ft(i) -= (ct_hat(i)-ctmp(i))/dct_hat(i);
00153     }
00154     //cout<<iter<<"\t"<<ft<<"\t"<<ct_hat-ct<<endl;
00155   }
00156   //cout<<ft<<"\t\t"<<ct<<"\t\t"<<ctmp<<endl;
00157
00158   return(ft);
00159 }
00160
00174 dvector get_ft( dvector& ct,const double& m, const dmatrix& V,const dvector& na, const dvector& wa )
00175 {
00176   /*  ARGUMENTS:
00177      ct is the observed catch for each fleet
00178      m is the natural mortality rate
00179      va is the vulnerability row fleet, col age va(1,ngear,1,nage)
00180      na is the start of year numbers at age
00181      wa is the mean weight-at-age
00182   */
00183
00184   int i,a,A;
00185   double minsurv = 0.05;
00186   int ng=size_count(ct);  //number of gears
00187   a=na.indexmin(); A=na.indexmax();
00188   dvector ft(1,ng); ft=0;
00189   dvector ctmp(1,ng);
00190   dvector ct_hat(1,ng); //predicted catch
00191   dvector dct_hat(1,ng);  //derivative
00192   dvector ba(a,A);    //biomass at age
00193   dvector ca(a,A);    //catch-at-age in numbers
00194   dvector zage(a,A);
00195   dvector sage(a,A);
00196   dvector ominus(a,A);
00197   dmatrix F(1,ng,a,A);
00198
00199
00200   //initial guess for ft=ct/(0.98 Bt);
00201   ctmp=ct;
00202   ba = elem_prod(na,wa);
00203   for(i=1;i<=ng;i++)
00204   {
00205     ft(i) = ctmp(i)/(0.98*ba*V(i)*exp(-0.5*m));
00206     if(1.-ft(i)<minsurv)
00207     {
00208       ft(i)=1.-minsurv;
00209       ctmp(i)=ft(i)*ba*V(i)*exp(-0.5*m);
00210     }
00211   }
00212   ct=ctmp;  //don't do this for the differentiable version.
00213
00214   //now solve baranov catch equation iteratively.
00215   for(int iter=1; iter<=17; iter++)
00216   {
00217     for(i=1;i<=ng;i++)F(i)=ft(i)*V(i);
00218     zage=m+colsum(F);
00219     sage=mfexp(-zage);
00220     ominus=(1.-sage);
00221
00222     for(i=1;i<=ng;i++)
00223     {
00224       dvector t1 = elem_div(F(i),zage);
00225       dvector t2 = elem_prod(t1,ominus);
00226       dvector t3 = elem_prod(ominus,na);
00227       ca = elem_prod(t2,na);
00228
00229       ct_hat(i) = ca*wa;
00230
00231       dvector t4 = ft(i)*square(V(i));
00232
00233       dvector t5 = elem_div(elem_prod(elem_prod(V(i),ominus),na),zage)
00234         - elem_div(elem_prod(elem_prod(t4,ominus),na),square(zage))
00235         + elem_div(elem_prod(elem_prod(t4,sage),na),zage);
00236       dct_hat(i) = t5*wa;
00237
00238         ft(i) -= (ct_hat(i)-ctmp(i))/dct_hat(i);
00239     }
00240     //cout<<iter<<"\t"<<ft<<"\t"<<ct_hat-ctmp<<endl;
00241     //SJDM, this algorithm does converge niceley for multiple fleets
00242   }
00243   //cout<<ft<<"\t\t"<<ct<<"\t\t"<<ctmp<<endl;
00244   //cout<<ct<<endl;
00245   return(ft);
00246 }
00247
```