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