Management Strategy Evaluation Demo  1.0
MSED
 All Classes Files Functions Variables Enumerations Enumerator Groups Pages
om.cpp
1  #undef REPORT
2  #define REPORT(object) report << #object "\n" << object << endl;
3  #undef COUT
4  #define COUT(object) cout<<fixed<<#object "\n"<<object<<endl;
5  #include <iostream>
6  #include <iomanip>
7  using namespace std;
8  #include <admodel.h>
9  #include <time.h>
10  //#include <statsLib.h>
11  #include "LRGS.h"
12  #include "MSYReferencePoints.h"
13  // #include "Scenario.h"
14  #include "OperatingModel.h"
15  // #include "harvestControlRule.h"
16 
17  time_t start,finish;
18  long hour,minute,second;
19  double elapsed_time;
20  // sLRGSparameters sPars;
21  sLRGSdata data;
22 
23 #include <admodel.h>
24 #include <contrib.h>
25 
26  extern "C" {
27  void ad_boundf(int i);
28  }
29 #include <om.htp>
30 
31 model_data::model_data(int argc,char * argv[]) : ad_comm(argc,argv)
32 {
33  int opt, on;
34  do_mse = 0;
35  rseed = 999;
36  if( ( on=option_match(ad_comm::argc,ad_comm::argv,"-mse",opt) ) >-1 )
37  {
38  do_mse = 1;
39  rseed=atoi(ad_comm::argv[on+1]);
40  COUT(rseed);
41  }
42  agek.allocate("agek");
43  syr.allocate("syr");
44  nyr.allocate("nyr");
45  iyr.allocate(syr,nyr,"iyr");
46  ct.allocate(syr,nyr,"ct");
47  it.allocate(syr,nyr,"it");
48  nScenario.allocate("nScenario");
49  n_hcr.allocate("n_hcr");
50  n_pyr.allocate("n_pyr");
51  n_flg_perfect_information.allocate("n_flg_perfect_information");
52  iuu_rate.allocate("iuu_rate");
53  min_tac.allocate("min_tac");
54  sEstimator.allocate("sEstimator");
55  data.syr = syr;
56  data.nyr = nyr;
57  data.agek = agek;
58  data.ct = ct;
59  data.it = it;
60  cout<<data.it<<endl;
61 }
62 
63 void model_parameters::initializationfunction(void)
64 {
65  log_bo.set_initial_value(8.0);
66  h.set_initial_value(0.9);
67  s.set_initial_value(0.9);
68  log_sigma.set_initial_value(4.65);
69  log_tau.set_initial_value(4.65);
70 }
71 
72 model_parameters::model_parameters(int sz,int argc,char * argv[]) :
73  model_data(argc,argv) , function_minimizer(sz)
74 {
75  initializationfunction();
76  log_bo.allocate(0,10,1,"log_bo");
77  h.allocate(0.2,1.0,1,"h");
78  s.allocate(0.0,1.0,1,"s");
79  log_sigma.allocate(2,"log_sigma");
80  log_tau.allocate(3,"log_tau");
81  wt.allocate(syr,nyr,-15,15,3,"wt");
82  f.allocate("f");
83  prior_function_value.allocate("prior_function_value");
84  likelihood_function_value.allocate("likelihood_function_value");
85  bo.allocate("bo");
86  #ifndef NO_AD_INITIALIZE
87  bo.initialize();
88  #endif
89  ro.allocate("ro");
90  #ifndef NO_AD_INITIALIZE
91  ro.initialize();
92  #endif
93  sig.allocate("sig");
94  #ifndef NO_AD_INITIALIZE
95  sig.initialize();
96  #endif
97  tau.allocate("tau");
98  #ifndef NO_AD_INITIALIZE
99  tau.initialize();
100  #endif
101  a.allocate("a");
102  #ifndef NO_AD_INITIALIZE
103  a.initialize();
104  #endif
105  b.allocate("b");
106  #ifndef NO_AD_INITIALIZE
107  b.initialize();
108  #endif
109  reck.allocate("reck");
110  #ifndef NO_AD_INITIALIZE
111  reck.initialize();
112  #endif
113  q.allocate("q");
114  #ifndef NO_AD_INITIALIZE
115  q.initialize();
116  #endif
117  fpen.allocate("fpen");
118  #ifndef NO_AD_INITIALIZE
119  fpen.initialize();
120  #endif
121  bt.allocate(syr,nyr+1,"bt");
122  #ifndef NO_AD_INITIALIZE
123  bt.initialize();
124  #endif
125  rt.allocate(syr,nyr,"rt");
126  #ifndef NO_AD_INITIALIZE
127  rt.initialize();
128  #endif
129  ft.allocate(syr,nyr,"ft");
130  #ifndef NO_AD_INITIALIZE
131  ft.initialize();
132  #endif
133  epsilon.allocate(syr,nyr,"epsilon");
134  #ifndef NO_AD_INITIALIZE
135  epsilon.initialize();
136  #endif
137  nll.allocate(1,8,"nll");
138  #ifndef NO_AD_INITIALIZE
139  nll.initialize();
140  #endif
141  sd_dep.allocate("sd_dep");
142 }
143 
144 void model_parameters::userfunction(void)
145 {
146  f =0.0;
147  sLRGSparameters sPars;
148  sPars.log_bo = log_bo;
149  sPars.h = h;
150  sPars.s = s;
151  sPars.log_sigma = log_sigma;
152  sPars.log_tau = log_tau;
153  sPars.wt = wt;
154  bo = mfexp(log_bo);
155  sig = sqrt(1.0/mfexp(log_sigma));
156  tau = sqrt(1.0/mfexp(log_tau));
157  reck = 4.*h/(1.-h);
158  // Testing out a new class called LRGS for doing all of the
159  // model calculations.
160  // LRGS cLRGSmodel(syr,nyr,agek,bo,h,s,sig,tau,ct,it,wt);
161  // Test cTest;
162  LRGS cLRGSmodel(data,sPars);
163  cLRGSmodel.initialize_model();
164  cLRGSmodel.population_dynamics();
165  cLRGSmodel.observation_model();
166  epsilon = cLRGSmodel.get_epsilon();
167  sd_dep = cLRGSmodel.get_depletion();
168  bt = cLRGSmodel.get_bt();
169  ft = cLRGSmodel.get_ft();
170  q = cLRGSmodel.get_q();
171  // cout<<q<<endl;
172  // initialize_model();
173  // population_dynamics();
174  // observation_model();
175  calc_objective_function();
176 }
177 
178 void model_parameters::calc_Reference_Points()
179 {
180  int i;
181  int j;
182  dvector fe(1,100);
183  double be,ce;
184  for( i = 1; i <= 100; i++ )
185  {
186  be = value(bo);
187  fe(i) = double((i-1.)/(100.-1.))*1.2;
188  for( j = 1; j <= 100; j++ )
189  {
190  ce = be * (1.-exp(-fe(i)));
191  be = value(s*be + a*be/(1.+b*be)) - ce;
192  }
193  cout<<setprecision(5)<<fe(i)<<"\t"<<ce<<"\t"<<be<<endl;
194  }
195 }
196 
197 void model_parameters::initialize_model()
198 {
199  rt.initialize();
200  bt.initialize();
201  bo = mfexp(log_bo);
202  ro = bo*(1.-s);
203  reck = 4.*h/(1.-h);
204  a = reck*ro/bo;
205  b = (reck-1.0)/bo;
206  rt(syr,syr+agek) = ro * exp(wt(syr,syr+agek));
207  bt(syr) = bo;
208  sig = sqrt(1.0/mfexp(log_sigma));
209  tau = sqrt(1.0/mfexp(log_tau));
210 }
211 
212 void model_parameters::population_dynamics()
213 {
214  int i;
215  dvariable btmp;
216  fpen.initialize();
217  for(i=syr;i<=nyr;i++)
218  {
219  ft(i) = -log((-ct(i)+bt(i))/bt(i));
220  if(i-syr > agek)
221  {
222  rt(i) = a*bt(i-agek)/(1.+b*bt(i-agek)) * exp(wt(i));
223  }
224  btmp = s*bt(i) + rt(i) - ct(i);
225  bt(i+1) = posfun(btmp,0.1,fpen);
226  }
227  sd_dep = bt(nyr)/bo;
228 }
229 
230 void model_parameters::observation_model()
231 {
232  int i;
233  dvar_vector zt = log(it) - log(bt(syr,nyr));
234  q = exp(mean(zt));
235  epsilon = zt - mean(zt);
236 }
237 
238 void model_parameters::run_mse()
239 {
240  int j;
241  Scenario cScenario1(agek,nScenario,n_pyr,n_flg_perfect_information,rseed,value(bo),
242  value(h),value(s),iuu_rate,
243  value(q),value(sig),value(tau),value(ft),
244  value(wt),it,ct,min_tac);
245  int e_hcr = n_hcr;
246  HarvestControlRule c_hcr(e_hcr);
247  EstimatorClass cEstimator(sEstimator);
248  OperatingModel cOM(cScenario1,cEstimator,c_hcr);
249  cOM.runMSEscenario(cScenario1);
250  ofstream ofs("OM.rep",ios::app);
251  ofs<<"t_bo\n" << cOM.get_bo() <<endl;
252  ofs<<"t_est_bo\n" << cOM.get_est_bo() <<endl;
253  ofs<<"t_bmsy\n" << cOM.get_bmsy() <<endl;
254  ofs<<"t_fmsy\n" << cOM.get_fmsy() <<endl;
255  ofs<<"t_msy\n" << cOM.get_msy() <<endl;
256  ofs<<"t_bt\n" << cOM.get_bt() <<endl;
257  ofs<<"t_aav\n" << cOM.get_aav() <<endl;
258  ofs.close();
259  // |--------------------------------------------|
260  // | NOW RUN THE MODEL WITH PERFECT INFORMATION |
261  // |--------------------------------------------|
262  Scenario cScenarioP(agek,nScenario,n_pyr,0,rseed,value(bo),
263  value(h),value(s),iuu_rate,
264  value(q),value(sig),value(tau),value(ft),
265  value(wt),it,ct);
266  OperatingModel cOMP(cScenarioP,cEstimator,c_hcr);
267  cOMP.runMSEscenario(cScenarioP);
268  ofstream ofs1("OM.rep",ios::app);
269  ofs1<<"p_bo\n" << cOMP.get_est_bo() <<endl;
270  ofs1<<"p_bmsy\n"<< cOMP.get_bmsy() <<endl;
271  ofs1<<"p_fmsy\n"<< cOMP.get_fmsy() <<endl;
272  ofs1<<"p_msy\n" << cOMP.get_msy() <<endl;
273  ofs1<<"p_bt\n" << cOMP.get_bt() <<endl;
274  ofs1<<"p_ct\n" << cOMP.get_hat_ct() <<endl;
275  ofs1<<"p_aav\n" << cOMP.get_aav() <<endl;
276  ofs1.close();
277 }
278 
279 void model_parameters::calc_objective_function()
280 {
281  nll.initialize();
282  dvariable isig2 = mfexp(log_sigma);
283  dvariable itau2 = mfexp(log_tau);
284  // No comments
285  nll(1) = dnorm(epsilon,sig);
286  nll(2) = dbeta(s,30.01,10.01);
287  nll(3) = dnorm(log_bo,log(3000),1.0);
288  nll(4) = dbeta((h-0.2)/0.8,1.01,1.01);
289  nll(5) = dgamma(isig2,1.01,1.01);
290  if(active(log_tau))
291  {
292  nll(6) = dgamma(itau2,1.01,1.01);
293  nll(7) = dnorm(wt,tau);
294  }
295  else
296  {
297  nll(7) = dnorm(wt,1.0);
298  }
299  if(fpen>0 && !mc_phase()) cout<<"Fpen = "<<fpen<<endl;
300  f = sum(nll) + 100000.*fpen;
301 }
302 
303 void model_parameters::report()
304 {
305  adstring ad_tmp=initial_params::get_reportfile_name();
306  ofstream report((char*)(adprogram_name + ad_tmp));
307  if (!report)
308  {
309  cerr << "error trying to open report file" << adprogram_name << ".rep";
310  return;
311  }
312  msy_reference_points cMSY(value(reck),value(s),value(bo));
313  REPORT(bo);
314  REPORT(h);
315  REPORT(s);
316  REPORT(q);
317  REPORT(sig);
318  REPORT(tau);
319  double fmsy = cMSY.get_fmsy();
320  double bmsy = cMSY.get_bmsy();
321  double msy = cMSY.get_msy();
322  REPORT(fmsy);
323  REPORT(bmsy);
324  REPORT(msy);
325  REPORT(ft);
326  REPORT(wt);
327  REPORT(bt);
328  REPORT(ct);
329  REPORT(epsilon);
330  // print mle estimates of key parameters for MSE
331  ofstream ofs("mse.par");
332  ofs<<bo<<"\n"<<reck<<"\n"<<s<<"\n"<<bt(nyr+1)<<endl;
333 }
334 
335 void model_parameters::final_calcs()
336 {
337  if(do_mse)
338  {
339  cout<<"Running MSE"<<endl;
340  run_mse();
341  }
342  time(&finish);
343  elapsed_time=difftime(finish,start);
344  hour=long(elapsed_time)/3600;
345  minute=long(elapsed_time)%3600/60;
346  second=(long(elapsed_time)%3600)%60;
347  cout<<endl<<endl<<"*******************************************"<<endl;
348  cout<<"--Start time: "<<ctime(&start)<<endl;
349  cout<<"--Finish time: "<<ctime(&finish)<<endl;
350  cout<<"--Runtime: ";
351  cout<<hour<<" hours, "<<minute<<" minutes, "<<second<<" seconds"<<endl;
352  cout<<"*******************************************"<<endl;
353  // Check calculus: A-O-K
354  // calcReferencePoints();
355 }
356 
357 void model_parameters::preliminary_calculations(void){
358 #if defined(USE_ADPVM)
359 
360  admaster_slave_variable_interface(*this);
361 
362 #endif
363 }
364 
365 model_data::~model_data()
366 {}
367 
368 model_parameters::~model_parameters()
369 {}
370 
371 void model_parameters::set_runtime(void){}
372 
373 #ifdef _BORLANDC_
374  extern unsigned _stklen=10000U;
375 #endif
376 
377 
378 #ifdef __ZTC__
379  extern unsigned int _stack=10000U;
380 #endif
381 
382  long int arrmblsize=0;
383 
384 int main(int argc,char * argv[])
385 {
386  ad_set_new_handler();
387  ad_exit=&ad_boundf;
388  time(&start);
389  arrmblsize = 50000000;
390  gradient_structure::set_GRADSTACK_BUFFER_SIZE(1.e7);
391  gradient_structure::set_CMPDIF_BUFFER_SIZE(1.e7);
392  gradient_structure::set_MAX_NVAR_OFFSET(5000);
393  gradient_structure::set_NUM_DEPENDENT_VARIABLES(5000);
394 
395  gradient_structure::set_NO_DERIVATIVES();
396  gradient_structure::set_YES_SAVE_VARIABLES_VALUES();
397  if (!arrmblsize) arrmblsize=15000000;
398  model_parameters mp(arrmblsize,argc,argv);
399  mp.iprint=10;
400  mp.preliminary_calculations();
401  mp.computations(argc,argv);
402  return 0;
403 }
404 
405 extern "C" {
406  void ad_boundf(int i)
407  {
408  /* so we can stop here */
409  exit(i);
410  }
411 }