Management Strategy Evaluation Demo  1.0
MSED
 All Classes Files Functions Variables Enumerations Enumerator Groups Pages
Public Member Functions | List of all members
OperatingModel Class Reference
Collaboration diagram for OperatingModel:
Collaboration graph

Public Member Functions

 OperatingModel (Scenario &cScenario, EstimatorClass &cEstimator, const HarvestControlRule &cHCR)
 
double get_bmsy ()
 
double get_fmsy ()
 
double get_msy ()
 
double get_bo ()
 
double get_est_bo ()
 
dvector get_bt ()
 
dvector get_hat_ct ()
 
dvector get_aav ()
 
void runMSEscenario (const Scenario &cScenario)
 Run the Operating model scenario. More...
 
void write_data_file (const int &nyr, const dvector &ct, const dvector &it)
 

Detailed Description

Definition at line 22 of file OperatingModel.h.

Member Function Documentation

void OperatingModel::runMSEscenario ( const Scenario &  cScenario)

Run the Operating model scenario. 

    This function runs the operating model conditional on the scenario class that is 
    passed as an arguement.
Author
Steve Martell
Date
May 29, 2013
Parameters
aScenario class object.
Returns
void
See Also

Definition at line 25 of file OperatingModel.cpp.

26 {
27  // This routine reconstructs the population dynamics based on the Scenario class
28  int i;
29  double ro, reck, a, b;
30  // [ ] TO DO, clean this up and declare this as private member variables.
31  int agek = m_agek;
32  double bo = m_bo;
33  double h = m_h;
34  double s = m_s;
35  double q = m_q;
36  double sig = m_sig;
37  double tau = m_tau;
38  dvector ft = m_ft;
39  dvector wt = m_wt;
40  dvector it = m_it;
41 
42  ro = bo*(1.-s);
43  reck = 4*h/(1.-h);
44  a = reck*ro/bo;
45  b = (reck-1.0)/bo;
46 
47  // |---------------------------------------------------------------------------------|
48  // | TRUE REFERENCE POINTS
49  // |---------------------------------------------------------------------------------|
50  // |
51  msy_reference_points cRefPoints(reck,s,bo);
52  m_fmsy = cRefPoints.get_fmsy();
53  m_bmsy = cRefPoints.get_bmsy();
54  m_msy = cRefPoints.get_msy();
55 
56  dvector bt(m_syr,m_nyr+m_pyr+1);
57  dvector rt(m_syr,m_nyr+m_pyr);
58  dvector hat_ct(m_syr,m_nyr+m_pyr);
59  dvector hat_dt(m_syr,m_nyr+m_pyr); // IUU Fishing
60  dvector hat_it(m_syr,m_nyr+m_pyr);
61 
62  // |---------------------------------------------------------------------------------|
63  // | CONDITION REFERENCE MODEL BASED ON SCENARIO INFORMATION
64  // |---------------------------------------------------------------------------------|
65  // |
66  hat_dt.initialize();
67  bt.initialize();
68  bt(m_syr) = bo;
69  rt(m_syr,m_syr+agek) = ro * exp(wt(m_syr,m_syr+agek));
70  hat_it(m_syr,m_nyr) = it;
71  for(i = m_syr; i <= m_nyr; i++)
72  {
73 
74  if(i-m_syr > agek)
75  {
76  rt(i) = a*bt(i-agek)/(1.+b*bt(i-agek)) * exp(wt(i));
77  }
78  // hat_ct(i) = bt(i) * (1.-mfexp(-ft(i)));
79  hat_ct(i) = bt(i) * ft(i);
80  bt(i+1) = s*bt(i) + rt(i) - hat_ct(i);
81  }
82  hat_ct(m_syr,m_nyr) = m_ct;
83  hat_dt(m_syr,m_nyr) = 0;
84  // | END OF CONDITIONING PERIOD
85 
86  // |---------------------------------------------------------------------------------|
87  // | PROJECTION PERIOD
88  // |---------------------------------------------------------------------------------|
89  // | -1. Set Random numbers based on random number seed.
90  // | -2. Calculate TAC based on harvest control rule and biomass estimate.
91  // | -3. Implement harvest on reference population.
92  // | -4. Update reference population.
93  // | -5. Update observation models & write data files.
94  // | -6. Conduct stock assessment & update Bo, reck and s parameters.
95  // | -7. Repeat steps 2-7 for pyr's
96  // | -8. Compute performance measures.
97  // | -9. Added sin curve to mimic non-stationarity
98  // |
99  int pyr1 = m_nyr+1;
100  int pyr2 = m_nyr+m_pyr;
101  dvector pyr(pyr1,pyr2);
102  pyr.fill_seqadd(pyr1,1);
103  int seed = m_rng;
104 
105  random_number_generator rng(seed);
106  dvector rt_dev(pyr1,pyr2);
107  dvector it_dev(pyr1,pyr2);
108  dvector pdo_dev(pyr1,pyr2);
109  dvector tac_dev(pyr1,pyr2);
110 
111  // | Recruitment Scenario
112  double lambda;
113  if( m_nScenario==1 )
114  {
115  lambda = 0;
116  }
117  else if( m_nScenario==2 )
118  {
119  lambda = 0.4;
120  }
121  pdo_dev = sin((pyr-double(pyr1))/double(pyr2-pyr1)*4.0*PI);
122 
123 
124  rt_dev.fill_randn(rng);
125  it_dev.fill_randn(rng);
126  tac_dev.fill_randn(rng);
127 
128  rt_dev = rt_dev*tau - 0.5*tau*tau;
129  it_dev = it_dev*sig - 0.5*sig*sig;
130  double phi = 0.04;
131  tac_dev = tac_dev*phi -0.5*phi*phi;
132 
133  double fmsy,bmsy,msy, tac, frate;
134  double est_bo = bo;
135  double est_reck = reck;
136  double est_s = s;
137  double est_bt = bt(pyr1);
138 
139 
140  for( i = pyr1; i <= pyr2; i++ )
141  {
142  // -2. Calculate TAC based on HCR & Biomass estimate.
143  msy_reference_points cRefPoints(est_reck,est_s,est_bo);
144  fmsy = cRefPoints.get_fmsy();
145  msy = cRefPoints.get_msy();
146  bmsy = cRefPoints.get_bmsy();
147  m_bmsy = bmsy;
148  m_msy = msy;
149  m_fmsy = fmsy;
150  tac = m_cHCR.getTac(est_bt,fmsy,msy,bmsy,est_bo,0.15,hat_ct(i-1),m_mintac);
151  tac = tac*exp(tac_dev(i));
152  // -3. Implement harvest on reference population, add implentation errors
153  // Watch out here if tac > bt(i), then need to set frate to some arbitrary max
154 
155  if( tac < bt(i) )
156  {
157  frate = -log((-tac+bt(i))/bt(i));
158  }
159  else
160  {
161  frate = 0.8;
162  }
163  hat_ct(i) = bt(i) * (1.-mfexp(-frate));
164  hat_dt(i) = bt(i) * (1.-mfexp(-m_iuu_rate));
165 
166  // -4. Update reference population
167  if(i-m_syr > agek)
168  {
169  rt(i) = a*bt(i-agek)/(1.+b*bt(i-agek)) * exp(rt_dev(i) + lambda*pdo_dev(i));
170  }
171  bt(i+1) = s*bt(i) + rt(i) - hat_ct(i) - hat_dt(i);
172 
173  // -5. Update observation models and write data files.
174  // cout<<"Q = "<<q<<" "<<m_q<<endl;
175  hat_it(i) = q*bt(i)*exp(it_dev(i));
176  write_data_file(i,hat_ct(m_syr,i),hat_it(m_syr,i));
177 
178  // -6. Conduct assessment and update parameters
179  // system("./OM -ind MSE.dat -nox -est > NUL");
180  // If using perfect information, then don't run assessment.
181  if(m_flg_perfect_information)
182  {
183  m_cEstimator.runEstimator();
184 
185  ifstream ifs("mse.par");
186  ifs>>est_bo;
187  ifs>>est_reck;
188  ifs>>est_s;
189  ifs>>est_bt;
190  m_est_bo = est_bo;
191  }
192  else if( !m_flg_perfect_information )
193  {
194  cout<<"|-- PERFECT INFORMATION --|"<<endl;
195  est_bo = m_bo;
196  est_reck = m_reck;
197  est_s = m_s;
198  est_bt = bt(i+1);
199  m_est_bo = est_bo;
200  }
201 
202 
203  // - Screen dump so you can watch the progress.
204  cout<<setprecision(3);
205  cout<<"|---------------------------|"<<endl;
206  cout<<"| - Year "<<i<<endl;
207  cout<<"|---------------------------|"<<endl;
208  cout<<"| - est Bo "<<est_bo<<endl;
209  cout<<"| - est Bt "<<est_bt<<endl;
210  cout<<"| - bmsy "<<bmsy<<endl;
211  cout<<"| - bt(i) "<<bt(i)<<endl;
212  cout<<"| - fmsy "<<fmsy<<endl;
213  cout<<"| - frate "<<frate<<endl;
214  cout<<"| - msy "<<msy<<endl;
215  cout<<"| - tac "<<tac<<endl;
216  cout<<"| - hat_ct(i) "<<hat_ct(i)<<endl;
217  cout<<"|---------------------------|"<<endl;
218 
219  } // 7. repeat steps 2-6
220  m_bt = bt(m_syr,m_nyr+m_pyr);
221  m_hat_ct = hat_ct;
222 
223  // 8. Calculate performance measures
224  // AAV = |Ct - Ct-1| / Ct
225  dvector AAV(m_syr,pyr2);
226  AAV.initialize();
227  for( i = m_syr+6; i <= pyr2; i++ )
228  {
229  AAV(i) = sum(fabs(first_difference(hat_ct(i-5,i))));
230  AAV(i) /= sum(hat_ct(i-5,i));
231  // cout<<AAV(i)<<endl;
232  }
233  m_AAV = AAV;
234 
235 }

Here is the call graph for this function:

The documentation for this class was generated from the following files: