Issue 107

src/nh99/mod_hess.cpp

@@ -60,6 +60,8 @@ void function_minimizer::hess_routine(void)
 
 void function_minimizer::hess_routine_noparallel(void)
 {
+  between_phases_calculations();
+
   int nvar=initial_params::nvarcalc(); // get the number of active parameters
   //if (adjm_ptr) set_labels_for_hess(nvar);
   independent_variables x(1,nvar);

tests/GNUmakefile

@@ -392,7 +392,7 @@ endif
 	$(MAKE) --directory=vectorize clean
 	$(MAKE) --directory=vonmises clean
 	$(MAKE) --directory=issue158 clean
-	$(MAKE) --directory=dd2 clean
+	$(MAKE) --directory=issue-107 clean
 
 dist-clean:
 ifeq ($(CMDSHELL),cmd)

tests/issue-107/GNUmakefile

@@ -0,0 +1,86 @@
+all:
+	../../admb pella_t_modified
+	-./pella_t_modified > out
+	../../admb pella_t_modified2
+	cp pella_t_modified.dat pella_t_modified2.dat
+	-./pella_t_modified2 > out2
+
+nohess:
+	../../admb pella_t_modified
+	-./pella_t_modified -nohess > out
+	../../admb pella_t_modified2
+	cp pella_t_modified.dat pella_t_modified2.dat
+	-./pella_t_modified2 -nohess > out2
+
+hess_step:
+	../../admb pella_t_modified
+	-./pella_t_modified -hess_step > out
+	../../admb pella_t_modified2
+	cp pella_t_modified.dat pella_t_modified2.dat
+	-./pella_t_modified2 -hess_step > out2
+
+
+diff:
+	sed 's/-0\.0000/0\.0000/g' pella_t_modified.cor > pella_t_modified.cor2
+	sed 's/-0\.0000/0\.0000/g' pella_t_modified2.cor > pella_t_modified2.cor2
+	-diff -w --side-by-side pella_t_modified.cor2 pella_t_modified2.cor2
+	-diff -w --side-by-side pella_t_modified.eva pella_t_modified2.eva
+	-diff -w --side-by-side pella_t_modified.log pella_t_modified2.log
+	-diff -w --side-by-side pella_t_modified.p01 pella_t_modified2.p01
+	-diff -w --side-by-side pella_t_modified.p02 pella_t_modified2.p02
+	-diff -w --side-by-side pella_t_modified.p03 pella_t_modified2.p03
+	-diff -w --side-by-side pella_t_modified.par pella_t_modified2.par
+	-diff -w --side-by-side pella_t_modified.std pella_t_modified2.std
+	-diff -w --side-by-side --width=180 -w out out2
+clean:
+	@rm -vf admodel.cov
+	@rm -vf admodel.dep
+	@rm -vf admodel.hes
+	@rm -vf cmpdiff.tmp
+	@rm -vf gradfil1.tmp
+	@rm -vf gradfil2.tmp
+	@rm -vf varssave.tmp
+	@rm -vf fmin.log out out2
+	@rm -vf pella_t_modified
+	@rm -vf pella_t_modified.b01
+	@rm -vf pella_t_modified.b02
+	@rm -vf pella_t_modified.b03
+	@rm -vf pella_t_modified.bar
+	@rm -vf pella_t_modified.cor
+	@rm -vf pella_t_modified.cor2
+	@rm -vf pella_t_modified.cpp
+	@rm -vf pella_t_modified.eva
+	@rm -vf pella_t_modified.htp
+	@rm -vf pella_t_modified.log
+	@rm -vf pella_t_modified.obj
+	@rm -vf pella_t_modified.p01
+	@rm -vf pella_t_modified.p02
+	@rm -vf pella_t_modified.p03
+	@rm -vf pella_t_modified.par
+	@rm -vf pella_t_modified.std
+	@rm -vf pella_t_modified2
+	@rm -vf pella_t_modified2.dat
+	@rm -vf pella_t_modified2.b01
+	@rm -vf pella_t_modified2.b02
+	@rm -vf pella_t_modified2.b03
+	@rm -vf pella_t_modified2.bar
+	@rm -vf pella_t_modified2.cor
+	@rm -vf pella_t_modified2.cor2
+	@rm -vf pella_t_modified2.cpp
+	@rm -vf pella_t_modified2.eva
+	@rm -vf pella_t_modified2.htp
+	@rm -vf pella_t_modified2.log
+	@rm -vf pella_t_modified2.obj
+	@rm -vf pella_t_modified2.p01
+	@rm -vf pella_t_modified2.p02
+	@rm -vf pella_t_modified2.p03
+	@rm -vf pella_t_modified2.par
+	@rm -vf pella_t_modified2.std
+	@rm -vf pella_t_modified.r01
+	@rm -vf pella_t_modified.r02
+	@rm -vf pella_t_modified.r03
+	@rm -vf pella_t_modified.rep
+	@rm -vf pella_t_modified2.r01
+	@rm -vf pella_t_modified2.r02
+	@rm -vf pella_t_modified2.r03
+	@rm -vf pella_t_modified2.rep

tests/issue-107/changes

@@ -0,0 +1,33 @@
+diff --git a/src/nh99/model5.cpp b/src/nh99/model5.cpp
+index 18a38a11..ae7177f9 100644
+--- a/src/nh99/model5.cpp
++++ b/src/nh99/model5.cpp
+@@ -18,10 +18,10 @@ void param_init_bounded_dev_vector::set_value(const dvar_vector& x,
+     {
+       ::set_value(*this,x,ii,minb,maxb,pen);
+     }
+-    dvariable s=mean(*this);
+-    pen+=10000.0*s*s;
+     if (!initial_params::mc_phase)
+     {
++      dvariable s=mean(*this);
++      //pen+=10000.0*s*s;
+       (*this)-=s;
+     }
+   }
+diff --git a/src/nh99/xmodelm3.cpp b/src/nh99/xmodelm3.cpp
+index 2db14566..7dc2f739 100644
+--- a/src/nh99/xmodelm3.cpp
++++ b/src/nh99/xmodelm3.cpp
+@@ -529,7 +529,10 @@ void tracing_message(int traceflag,const char *s);
+       report(g);
+       // in case the user changes some initial_params in the report section
+       // call reset again
+-      initial_params::reset(dvar_vector(x));
++      if (initial_params::current_phase < initial_params::max_number_phases)
++      {
++        initial_params::reset(dvar_vector(x));
++      }
+       report_function_minimizer_stats();
+       if (quit_flag=='Q') break;
+       if (!quit_flag || quit_flag == 'N')

tests/issue-107/pella_t_modified.dat

@@ -0,0 +1,63 @@
+# number of years
+#   of data
+       56
+ 32            44.45399      59.49697
+ 33            46.795        62.28313 
+ 34            47.49         65.33897 
+ 35            47.34299      75.53572 
+ 36            48.923        70.38986 
+ 37            49.53899      79.26112 
+ 38            49.55299      87.31544 
+ 39            50.903        79.23817 
+ 40            53.381        81.68614 
+ 41            52.231        85.15594 
+ 42            50.38799      89.97632 
+ 43            53.699        95.99872 
+ 44            53.435        110.9247 
+ 45            53.395        101.9471 
+ 46            60.266        100.8948 
+ 47            55.7          99.15089 
+ 48            55.564        100.0702 
+ 49            55.02499      95.47215 
+ 50            57.23399      95.66711 
+ 51            56.045        95.57533 
+ 52            62.26199      107.7313 
+ 53            59.837        131.5304 
+ 54            70.58298      133.3808 
+ 55            57.521        119.6114 
+ 56            66.681        130.005 
+ 57            60.854        109.6818 
+ 58            64.51399      120.868 
+ 59            71.24298      130.0228 
+ 60            86.62286      133.8806 
+ 61            92.56155      127.611 
+ 62            110.0678      115.1604 
+ 63            103.5499      106.4599 
+ 64            94.96712      101.7176 
+ 65            100.529       99.29117 
+ 66            94.47342      100.3763 
+ 67            88.63864      100.822 
+ 68            79.19131      103.2759 
+ 69            85.24725      95.86926 
+ 70            82.78796      93.34418 
+ 71            81.88229      89.67853 
+ 72            76.74669      80.95505 
+ 73            62.64147      64.91022 
+ 74            52.49203      62.3805 
+ 75            48.93383      63.43256 
+ 76            50.60567      53.82018 
+ 77            42.80426      60.51315 
+ 78            41.80861      67.21815 
+ 79            47.3311       70.05229 
+ 80            52.06058      92.7 
+ 81            48.50106      115.7 
+ 82            46.63046      114.4 
+ 83            55.82103      124.4 
+ 84            59.72166      135.6 
+ 85            70.6          140.7 
+ 86            86.6          132.6 
+ 87            86.4          125.4 
+-1,-1,-1
+0,0,0
+0,0,0
+Pacific Halibut supplied by Rick Deriso

tests/issue-107/pella_t_modified.tpl

@@ -0,0 +1,157 @@
+// Author: David Fournier
+// Copyright (c) 2009-2018 ADMB Foundation
+
+DATA_SECTION
+  init_int nobs;
+  init_matrix data(1,nobs,1,3)
+  vector obs_catch(1,nobs);
+  vector cpue(1,nobs);
+  vector effort(1,nobs);
+  number avg_effort
+
+  // set up the mcmc chain counter
+  int mceval_counter;
+  !! mceval_counter = 0;
+  int mcmc_counter;
+  !! mcmc_counter = -2;
+  int count;
+  !! count = 0;
+
+INITIALIZATION_SECTION
+  m 2.
+  beta 1.
+  r 1.
+PARAMETER_SECTION
+  init_bounded_number r(0.,5,2)
+  init_bounded_number beta(0.,5.)
+  init_number log_binit(2)
+  init_bounded_number q(0.,1.)
+  init_bounded_number m(1,10.,4)
+  init_bounded_dev_vector effort_devs(1,nobs,-5.,5.,3)
+  init_bounded_vector k_devs(2,nobs,-5.,5.,4)
+  number binit
+  vector pred_catch(1,nobs)
+  vector biomass(1,nobs)
+  vector f(1,nobs)
+  vector k(1,nobs)
+  vector k_trend(1,nobs)
+  sdreport_number k_1
+  sdreport_number k_last
+  sdreport_number k_change
+  sdreport_number k_ratio
+  sdreport_number B_projected
+  number tmp_mort;
+  number bio_tmp;
+  number c_tmp;
+  objective_function_value ff;
+PRELIMINARY_CALCS_SECTION
+  // get the data out of the data matrix into
+  obs_catch=column(data,2);
+  cpue=column(data,3);
+  // divide the catch by the cpue to get the effort
+  effort=elem_div(obs_catch,cpue);
+  // normalize the effort and save the average
+  double avg_effort=mean(effort);
+  effort/=avg_effort;
+  cout << " beta" << beta << endl;
+PROCEDURE_SECTION
+  // calculate the fishing mortality
+  calculate_fishing_mortality();
+  // calculate the biomass and predicted catch
+  calculate_biomass_and_predicted_catch();
+  // calculate the objective function
+  calculate_the_objective_function();
+
+  // MCMC output
+  if (initial_params::mc_phase==1){
+    cout << " MCMC: " << ++mcmc_counter;
+  }
+  if (mceval_phase()){
+    cout << " MCeval: " << ++mceval_counter;
+  }
+
+FUNCTION calculate_fishing_mortality
+  // calculate the fishing mortality
+  f=q*effort;
+  if (active(effort_devs)) f=elem_prod(f,exp(effort_devs));
+
+FUNCTION calculate_biomass_and_predicted_catch
+  // calculate the biomass and predicted catch
+  if (!active(log_binit))
+  {
+    log_binit=log(obs_catch(1)/(q*effort(1)));
+  }
+  binit=exp(log_binit);
+  biomass[1]=binit;
+  if (active(k_devs))
+  {
+    k(1)=binit/beta;
+    for (int i=2;i<=nobs;i++)
+    {
+      k(i)=k(i-1)*exp(k_devs(i));
+    }
+  }
+  else
+  {
+    // set the whole vector equal to the constant k value
+    k=binit/beta;
+  }
+  // only calculate these for the standard deviation report
+  if (sd_phase())
+  {
+    k_1=k(1);
+    k_last=k(nobs);
+    k_ratio=k(nobs)/k(1);
+    k_change=k(nobs)-k(1);
+  }
+  // two time steps per year desired
+  int nsteps=2;
+  double delta=1./nsteps;
+  // Integrate the logistic dynamics over n time steps per year
+  for (int i=1; i<=nobs; i++)
+  {
+    bio_tmp=1.e-20+biomass[i];
+    c_tmp=0.;
+    for (int j=1; j<=nsteps; j++)
+    {
+      //This is the new biomass after time step delta
+      bio_tmp=bio_tmp*(1.+r*delta)/
+	(1.+ (r*pow(bio_tmp/k(i),m-1.)+f(i))*delta );
+      // This is the catch over time step delta
+      c_tmp+=f(i)*delta*bio_tmp;
+    }
+    pred_catch[i]=c_tmp;        // This is the catch for this year
+    if (i<nobs)
+    {
+      biomass[i+1]=bio_tmp;// This is the biomass at the begining of the next
+    }			   // year
+    else
+    {
+      B_projected=bio_tmp; // get the projected biomass for std dev report
+    }
+  }
+
+FUNCTION calculate_the_objective_function
+  if (!active(effort_devs))
+  {
+    ff=nobs/2.*log(norm2(log(obs_catch)-log(1.e-10+pred_catch)));
+  }
+  else if(!active(k_devs))
+  {
+    ff= .5*(size_count(obs_catch)+size_count(effort_devs))*
+      log(norm2(log(obs_catch)-log(1.e-10+pred_catch))
+			      +0.1*norm2(effort_devs));
+  }
+  else
+  {
+    ff= .5*( size_count(obs_catch)+size_count(effort_devs)
+	     +size_count(k_devs) )*
+      log(norm2(log(obs_catch)-log(pred_catch))
+	   + 0.1*norm2(effort_devs)+10.*norm2(k_devs));
+  }
+  // Bayesian contribution for Pella Tomlinson m
+  ff+=2.*square(log(m-1.));
+  if (initial_params::current_phase<3)
+  {
+    ff+=1000.*square(log(mean(f)/.4));
+  }