Annotation of /trunk/gadget/maturity.cc

Parent Directory | Revision Log

---

Revision 1 - (view) (download)

1 :	agomez	1	#include "maturity.h"
2 :			#include "stock.h"
3 :			#include "mathfunc.h"
4 :			#include "readfunc.h"
5 :			#include "readword.h"
6 :			#include "conversionindex.h"
7 :			#include "errorhandler.h"
8 :			#include "gadget.h"
9 :			#include "global.h"
10 :
11 :			// ********************************************************
12 :			// Functions for base Maturity
13 :			// ********************************************************
14 :			Maturity::Maturity(const IntVector& tmpareas, int minage, int numage,
15 :			const LengthGroupDivision* const lgrpdiv, const char* givenname)
16 :			: HasName(givenname), LivesOnAreas(tmpareas) {
17 :
18 :			int i;
19 :			istagged = 0;
20 :			tmpratio = 1.0;
21 :			ratioscale = 1.0; //JMB used to scale the ratios to ensure that they sum to 1
22 :			LgrpDiv = new LengthGroupDivision(*lgrpdiv);
23 :			if (LgrpDiv->Error())
24 :			handle.logMessage(LOGFAIL, "Error in maturity - failed to create length group");
25 :
26 :			IntVector lower(numage, 0);
27 :			IntVector agesize(numage, LgrpDiv->numLengthGroups());
28 :			Storage.resize(areas.Size(), minage, lower, agesize);
29 :			for (i = 0; i < Storage.Size(); i++)
30 :			Storage[i].setToZero();
31 :			}
32 :
33 :			Maturity::~Maturity() {
34 :			int i;
35 :			for (i = 0; i < matureStockNames.Size(); i++)
36 :			delete[] matureStockNames[i];
37 :			for (i = 0; i < CI.Size(); i++)
38 :			delete CI[i];
39 :			delete LgrpDiv;
40 :			}
41 :
42 :			void Maturity::setStock(StockPtrVector& stockvec) {
43 :			int i, j, index;
44 :
45 :			for (i = 0; i < matureStockNames.Size(); i++)
46 :			for (j = 0; j < matureStockNames.Size(); j++)
47 :			if ((strcasecmp(matureStockNames[i], matureStockNames[j]) == 0) && (i != j))
48 :			handle.logMessage(LOGFAIL, "Error in maturity - repeated stock", matureStockNames[i]);
49 :
50 :			for (i = 0; i < stockvec.Size(); i++)
51 :			for (j = 0; j < matureStockNames.Size(); j++)
52 :			if (strcasecmp(stockvec[i]->getName(), matureStockNames[j]) == 0)
53 :			matureStocks.resize(stockvec[i]);
54 :
55 :			if (matureStocks.Size() != matureStockNames.Size()) {
56 :			handle.logMessage(LOGWARN, "Error in maturity - failed to match mature stocks");
57 :			for (i = 0; i < stockvec.Size(); i++)
58 :			handle.logMessage(LOGWARN, "Error in maturity - found stock", stockvec[i]->getName());
59 :			for (i = 0; i < matureStockNames.Size(); i++)
60 :			handle.logMessage(LOGWARN, "Error in maturity - looking for stock", matureStockNames[i]);
61 :			handle.logMessage(LOGFAIL, ""); //JMB this will exit gadget
62 :			}
63 :
64 :			//JMB ensure that the ratio vector is indexed in the right order
65 :			ratioindex.resize(matureStocks.Size(), 0);
66 :			for (i = 0; i < matureStocks.Size(); i++)
67 :			for (j = 0; j < matureStockNames.Size(); j++)
68 :			if (strcasecmp(matureStocks[i]->getName(), matureStockNames[j]) == 0)
69 :			ratioindex[i] = j;
70 :
71 :			for (i = 0; i < matureStocks.Size(); i++) {
72 :			CI.resize(new ConversionIndex(LgrpDiv, matureStocks[i]->getLengthGroupDiv()));
73 :			if (CI[i]->Error())
74 :			handle.logMessage(LOGFAIL, "Error in maturity - error when checking length structure");
75 :
76 :			index = 0;
77 :			for (j = 0; j < areas.Size(); j++)
78 :			if (!matureStocks[i]->isInArea(areas[j]))
79 :			index++;
80 :
81 :			if (index != 0)
82 :			handle.logMessage(LOGWARN, "Warning in maturity - mature stock isnt defined on all areas");
83 :			}
84 :			}
85 :
86 :			void Maturity::Print(ofstream& outfile) const {
87 :			int i;
88 :			outfile << "\nMaturity\n\tNames of mature stocks:";
89 :			for (i = 0; i < matureStockNames.Size(); i++)
90 :			outfile << sep << matureStockNames[i];
91 :			outfile << "\n\tRatio maturing into each stock:";
92 :			for (i = 0; i < matureRatio.Size(); i++)
93 :			outfile << sep << (matureRatio[ratioindex[i]] * ratioscale);
94 :			outfile << "\n\tStored numbers:\n";
95 :			for (i = 0; i < areas.Size(); i++) {
96 :			outfile << "\tInternal area " << areas[i] << endl;
97 :			Storage[i].printNumbers(outfile);
98 :			}
99 :			}
100 :
101 :			void Maturity::Reset(const TimeClass* const TimeInfo) {
102 :			//JMB check that the sum of the ratios is 1
103 :			if (TimeInfo->getTime() == 1) {
104 :			int i;
105 :			ratioscale = 0.0;
106 :			for (i = 0; i < matureRatio.Size(); i++ )
107 :			ratioscale += matureRatio[i];
108 :
109 :			if (isZero(ratioscale)) {
110 :			handle.logMessage(LOGWARN, "Warning in maturity - specified ratios are zero");
111 :			ratioscale = 1.0;
112 :			} else if (isEqual(ratioscale, 1.0)) {
113 :			// do nothing
114 :			} else {
115 :			handle.logMessage(LOGWARN, "Warning in maturity - scaling ratios using", ratioscale);
116 :			ratioscale = 1.0 / ratioscale;
117 :			}
118 :			}
119 :			}
120 :
121 :			void Maturity::Move(int area, const TimeClass* const TimeInfo) {
122 :			if (!(this->isMaturationStep(TimeInfo)))
123 :			handle.logMessage(LOGFAIL, "Error in maturity - maturity requested on wrong timestep");
124 :			int i, inarea = this->areaNum(area);
125 :			for (i = 0; i < matureStocks.Size(); i++) {
126 :			if (!matureStocks[i]->isInArea(area))
127 :			handle.logMessage(LOGFAIL, "Error in maturity - mature stock doesnt live on area", area);
128 :
129 :			tmpratio = matureRatio[ratioindex[i]] * ratioscale;
130 :			matureStocks[i]->Add(Storage[inarea], CI[i], area, tmpratio);
131 :			if (istagged && tagStorage.numTagExperiments() > 0)
132 :			matureStocks[i]->Add(tagStorage, CI[i], area, tmpratio);
133 :			}
134 :
135 :			Storage[inarea].setToZero();
136 :			if (istagged && tagStorage.numTagExperiments() > 0)
137 :			tagStorage[inarea].setToZero();
138 :			}
139 :
140 :			void Maturity::storeMatureStock(int area, int age, int length, double number, double weight) {
141 :			if (isZero(number) || isZero(weight)) {
142 :			Storage[this->areaNum(area)][age][length].setToZero();
143 :			} else {
144 :			Storage[this->areaNum(area)][age][length].N = number;
145 :			Storage[this->areaNum(area)][age][length].W = weight;
146 :			}
147 :			}
148 :
149 :			void Maturity::storeMatureTagStock(int area, int age, int length, double number, int id) {
150 :			if (!istagged)
151 :			handle.logMessage(LOGFAIL, "Error in maturity - invalid tagging experiment");
152 :			if ((id >= tagStorage.numTagExperiments()) || (id < 0))
153 :			handle.logMessage(LOGFAIL, "Error in maturity - invalid tagging experiment");
154 :			if (isZero(number))
155 :			*(tagStorage[this->areaNum(area)][age][length][id].N) = 0.0;
156 :			else
157 :			*(tagStorage[this->areaNum(area)][age][length][id].N) = number;
158 :			}
159 :
160 :			const StockPtrVector& Maturity::getMatureStocks() {
161 :			return matureStocks;
162 :			}
163 :
164 :			void Maturity::setTagged() {
165 :			istagged = 1;
166 :			//resize tagStorage to be the same size as Storage
167 :			int i, minage, maxage;
168 :			minage = Storage[0].minAge();
169 :			maxage = Storage[0].maxAge();
170 :			IntVector lower(maxage - minage + 1, 0);
171 :			IntVector size(maxage - minage + 1, LgrpDiv->numLengthGroups());
172 :			tagStorage.resize(areas.Size(), minage, lower, size);
173 :			for (i = 0; i < tagStorage.Size(); i++)
174 :			tagStorage[i].setToZero();
175 :			}
176 :
177 :			void Maturity::addMaturityTag(const char* tagname) {
178 :			if (!istagged)
179 :			handle.logMessage(LOGFAIL, "Error in maturity - invalid tagging experiment", tagname);
180 :			tagStorage.addTag(tagname);
181 :			}
182 :
183 :			void Maturity::deleteMaturityTag(const char* tagname) {
184 :			if (!istagged)
185 :			handle.logMessage(LOGFAIL, "Error in maturity - invalid tagging experiment", tagname);
186 :
187 :			int minage, maxage, age, len, a;
188 :			int id = tagStorage.getTagID(tagname);
189 :
190 :			if (id >= 0) {
191 :			minage = tagStorage[0].minAge();
192 :			maxage = tagStorage[0].maxAge();
193 :			//free memory allocated for tagging experiment
194 :			for (a = 0; a < tagStorage.Size(); a++) {
195 :			for (age = minage; age <= maxage; age++) {
196 :			for (len = tagStorage[a].minLength(age); len < tagStorage[a].maxLength(age); len++) {
197 :			delete[] (tagStorage[a][age][len][id].N);
198 :			(tagStorage[a][age][len][id].N) = NULL;
199 :			}
200 :			}
201 :			}
202 :			tagStorage.deleteTag(tagname);
203 :
204 :			} else
205 :			handle.logMessage(LOGWARN, "Warning in maturity - failed to delete tagging experiment", tagname);
206 :			}
207 :
208 :			// ********************************************************
209 :			// Functions for MaturityA
210 :			// ********************************************************
211 :			MaturityA::MaturityA(CommentStream& infile, const TimeClass* const TimeInfo,
212 :			Keeper* const keeper, int minage, int numage, const IntVector& tmpareas,
213 :			const char* givenname, const LengthGroupDivision* const lgrpdiv)
214 :			: Maturity(tmpareas, minage, numage, lgrpdiv, givenname), minStockAge(minage) {
215 :
216 :			//JMB store the length of the timestep
217 :			timesteplength = TimeInfo->getTimeStepSize();
218 :
219 :			char text[MaxStrLength];
220 :			strncpy(text, "", MaxStrLength);
221 :			int i;
222 :
223 :			keeper->addString("maturity");
224 :			infile >> text >> ws;
225 :			if ((strcasecmp(text, "nameofmaturestocksandratio") != 0) && (strcasecmp(text, "maturestocksandratios") != 0))
226 :			handle.logFileUnexpected(LOGFAIL, "maturestocksandratios", text);
227 :
228 :			i = 0;
229 :			infile >> text >> ws;
230 :			while (strcasecmp(text, "coefficients") != 0 && !infile.eof()) {
231 :			matureStockNames.resize(new char[strlen(text) + 1]);
232 :			strcpy(matureStockNames[i], text);
233 :			matureRatio.resize(1, keeper);
234 :			if (!(infile >> matureRatio[i]))
235 :			handle.logFileMessage(LOGFAIL, "invalid format for mature ratio");
236 :			matureRatio[i].Inform(keeper);
237 :
238 :			infile >> text >> ws;
239 :			i++;
240 :			}
241 :
242 :			if (infile.eof())
243 :			handle.logFileEOFMessage(LOGFAIL);
244 :			maturityParameters.resize(4, keeper);
245 :			maturityParameters.read(infile, TimeInfo, keeper);
246 :			preCalcMaturation.AddRows(numage, LgrpDiv->numLengthGroups(), 0.0);
247 :
248 :			infile >> ws;
249 :			if (!infile.eof()) {
250 :			infile >> text >> ws;
251 :			handle.logFileUnexpected(LOGFAIL, "<end of file>", text);
252 :			}
253 :			handle.logMessage(LOGMESSAGE, "Read maturity data file");
254 :			keeper->clearLast();
255 :			}
256 :
257 :			void MaturityA::Reset(const TimeClass* const TimeInfo) {
258 :			Maturity::Reset(TimeInfo);
259 :
260 :			if (TimeInfo->didStepSizeChange())
261 :			timesteplength = TimeInfo->getTimeStepSize();
262 :
263 :			maturityParameters.Update(TimeInfo);
264 :			if (maturityParameters.didChange(TimeInfo)) {
265 :			if (maturityParameters[1] < LgrpDiv->minLength())
266 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - l50 less than minimum length for stock", this->getName());
267 :			if (maturityParameters[1] > LgrpDiv->maxLength())
268 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - l50 greater than maximum length for stock", this->getName());
269 :
270 :			int age, len;
271 :			for (age = 0; age < preCalcMaturation.Nrow(); age++) {
272 :			for (len = 0; len < LgrpDiv->numLengthGroups(); len++) {
273 :			tmpratio = exp(-maturityParameters[0] * (LgrpDiv->meanLength(len) - maturityParameters[1]) -
274 :			maturityParameters[2] * (age + minStockAge - maturityParameters[3]));
275 :			preCalcMaturation[age][len] = 1.0 / (1.0 + tmpratio);
276 :			}
277 :			}
278 :
279 :			if (handle.getLogLevel() >= LOGMESSAGE)
280 :			handle.logMessage(LOGMESSAGE, "Reset maturity data for stock", this->getName());
281 :			}
282 :			}
283 :
284 :			void MaturityA::setStock(StockPtrVector& stockvec) {
285 :			Maturity::setStock(stockvec);
286 :
287 :			int i;
288 :			minMatureAge = 9999;
289 :			double minlength = 9999.0;
290 :			for (i = 0; i < matureStocks.Size(); i++) {
291 :			minMatureAge = min(matureStocks[i]->minAge(), minMatureAge);
292 :			minlength = min(matureStocks[i]->getLengthGroupDiv()->minLength(), minlength);
293 :			}
294 :			minMatureLength = LgrpDiv->numLengthGroup(minlength);
295 :			if (minMatureAge < minStockAge)
296 :			handle.logMessage(LOGFAIL, "Error in maturity - minimum mature age is less than stock age for stock", this->getName());
297 :			}
298 :
299 :			double MaturityA::calcMaturation(int age, int length, int growth, double weight) {
300 :
301 :			if ((age >= minMatureAge) && ((length + growth) >= minMatureLength)) {
302 :			tmpratio = preCalcMaturation[age - minStockAge][length] *
303 :			(maturityParameters[0] * growth * LgrpDiv->dl() + maturityParameters[2] * timesteplength);
304 :			return (min(max(0.0, tmpratio), 1.0));
305 :			}
306 :			return 0.0;
307 :			}
308 :
309 :			void MaturityA::Print(ofstream& outfile) const {
310 :			Maturity::Print(outfile);
311 :			outfile << "\tPrecalculated maturity:\n";
312 :			preCalcMaturation.Print(outfile);
313 :			}
314 :
315 :			int MaturityA::isMaturationStep(const TimeClass* const TimeInfo) {
316 :			return 1;
317 :			}
318 :
319 :			// ********************************************************
320 :			// Functions for MaturityB
321 :			// ********************************************************
322 :			MaturityB::MaturityB(CommentStream& infile, const TimeClass* const TimeInfo,
323 :			Keeper* const keeper, int minage, int numage, const IntVector& tmpareas,
324 :			const char* givenname, const LengthGroupDivision* const lgrpdiv)
325 :			: Maturity(tmpareas, minage, numage, lgrpdiv, givenname) {
326 :
327 :			//JMB set the default value for currenttimestep
328 :			currentmaturitystep = 0;
329 :
330 :			char text[MaxStrLength];
331 :			strncpy(text, "", MaxStrLength);
332 :			int i, tmpint = 0;
333 :
334 :			keeper->addString("maturity");
335 :			infile >> text >> ws;
336 :			if ((strcasecmp(text, "nameofmaturestocksandratio") != 0) && (strcasecmp(text, "maturestocksandratios") != 0))
337 :			handle.logFileUnexpected(LOGFAIL, "maturestocksandratios", text);
338 :
339 :			i = 0;
340 :			infile >> text >> ws;
341 :			while (strcasecmp(text, "maturitysteps") != 0 && !infile.eof()) {
342 :			matureStockNames.resize(new char[strlen(text) + 1]);
343 :			strcpy(matureStockNames[i], text);
344 :			matureRatio.resize(1, keeper);
345 :			if (!(infile >> matureRatio[i]))
346 :			handle.logFileMessage(LOGFAIL, "invalid format for mature ratio");
347 :			matureRatio[i].Inform(keeper);
348 :
349 :			infile >> text >> ws;
350 :			i++;
351 :			}
352 :
353 :			if (infile.eof())
354 :			handle.logFileEOFMessage(LOGFAIL);
355 :
356 :			infile >> ws;
357 :			while (isdigit(infile.peek()) && !infile.eof()) {
358 :			infile >> tmpint >> ws;
359 :			maturitystep.resize(1, tmpint);
360 :			}
361 :			if (infile.eof())
362 :			handle.logFileEOFMessage(LOGFAIL);
363 :
364 :			i = 0;
365 :			infile >> text;
366 :			if (strcasecmp(text, "maturitylengths") != 0)
367 :			handle.logFileUnexpected(LOGFAIL, "maturitylengths", text);
368 :			while (i < maturitystep.Size() && !infile.eof()) {
369 :			maturitylength.resize(1, keeper);
370 :			maturitylength[i++].read(infile, TimeInfo, keeper);
371 :			}
372 :
373 :			for (i = 0; i < maturitystep.Size(); i++)
374 :			if (maturitystep[i] < 1 || maturitystep[i] > TimeInfo->numSteps())
375 :			handle.logFileMessage(LOGFAIL, "invalid maturity step", maturitystep[i]);
376 :
377 :			if (maturitylength.Size() != maturitystep.Size())
378 :			handle.logFileMessage(LOGFAIL, "number of maturitysteps does not equal number of maturitylengths");
379 :
380 :			infile >> ws;
381 :			if (!infile.eof()) {
382 :			infile >> text >> ws;
383 :			handle.logFileUnexpected(LOGFAIL, "<end of file>", text);
384 :			}
385 :			handle.logMessage(LOGMESSAGE, "Read maturity data file");
386 :			keeper->clearLast();
387 :			}
388 :
389 :			void MaturityB::Print(ofstream& outfile) const {
390 :			int i;
391 :			Maturity::Print(outfile);
392 :			outfile << "\tMaturity timesteps";
393 :			for (i = 0; i < maturitystep.Size(); i++)
394 :			outfile << sep << maturitystep[i];
395 :			outfile << "\n\tMaturity lengths";
396 :			for (i = 0; i < maturitylength.Size(); i++)
397 :			outfile << sep << maturitylength[i];
398 :			outfile << endl;
399 :			}
400 :
401 :			void MaturityB::setStock(StockPtrVector& stockvec) {
402 :			Maturity::setStock(stockvec);
403 :			}
404 :
405 :			void MaturityB::Reset(const TimeClass* const TimeInfo) {
406 :			Maturity::Reset(TimeInfo);
407 :
408 :			int i;
409 :			maturitylength.Update(TimeInfo);
410 :			if (maturitylength.didChange(TimeInfo)) {
411 :			for (i = 0; i < maturitylength.Size(); i++) {
412 :			if (maturitylength[i] < LgrpDiv->minLength())
413 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - length less than minimum stock length for stock", this->getName());
414 :			if (maturitylength[i] > LgrpDiv->maxLength())
415 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - length greater than maximum stock length for stock", this->getName());
416 :			}
417 :
418 :			if (handle.getLogLevel() >= LOGMESSAGE)
419 :			handle.logMessage(LOGMESSAGE, "Reset maturity data for stock", this->getName());
420 :			}
421 :			}
422 :
423 :			double MaturityB::calcMaturation(int age, int length, int growth, double weight) {
424 :			if ((LgrpDiv->meanLength(length) > maturitylength[currentmaturitystep]) ||
425 :			(isEqual(LgrpDiv->meanLength(length), maturitylength[currentmaturitystep])))
426 :			return 1.0;
427 :			return 0.0;
428 :			}
429 :
430 :			int MaturityB::isMaturationStep(const TimeClass* const TimeInfo) {
431 :			int i;
432 :			for (i = 0; i < maturitystep.Size(); i++) {
433 :			if (maturitystep[i] == TimeInfo->getStep()) {
434 :			currentmaturitystep = i;
435 :			return 1;
436 :			}
437 :			}
438 :			return 0;
439 :			}
440 :
441 :			// ********************************************************
442 :			// Functions for MaturityC
443 :			// ********************************************************
444 :			MaturityC::MaturityC(CommentStream& infile, const TimeClass* const TimeInfo,
445 :			Keeper* const keeper, int minage, int numage, const IntVector& tmpareas,
446 :			const char* givenname, const LengthGroupDivision* const lgrpdiv, int numMatConst)
447 :			: Maturity(tmpareas, minage, numage, lgrpdiv, givenname), minStockAge(minage) {
448 :
449 :			char text[MaxStrLength];
450 :			strncpy(text, "", MaxStrLength);
451 :			int i, tmpint = 0;
452 :
453 :			keeper->addString("maturity");
454 :			infile >> text >> ws;
455 :			if ((strcasecmp(text, "nameofmaturestocksandratio") != 0) && (strcasecmp(text, "maturestocksandratios") != 0))
456 :			handle.logFileUnexpected(LOGFAIL, "maturestocksandratios", text);
457 :
458 :			i = 0;
459 :			infile >> text >> ws;
460 :			while (strcasecmp(text, "coefficients") != 0 && !infile.eof()) {
461 :			matureStockNames.resize(new char[strlen(text) + 1]);
462 :			strcpy(matureStockNames[i], text);
463 :			matureRatio.resize(1, keeper);
464 :			if (!(infile >> matureRatio[i]))
465 :			handle.logFileMessage(LOGFAIL, "invalid format for mature ratio");
466 :			matureRatio[i].Inform(keeper);
467 :
468 :			infile >> text >> ws;
469 :			i++;
470 :			}
471 :
472 :			if (infile.eof())
473 :			handle.logFileEOFMessage(LOGFAIL);
474 :			maturityParameters.resize(numMatConst, keeper);
475 :			maturityParameters.read(infile, TimeInfo, keeper);
476 :			preCalcMaturation.AddRows(numage, LgrpDiv->numLengthGroups(), 0.0);
477 :
478 :			infile >> text >> ws;
479 :			if ((strcasecmp(text, "maturitystep") != 0) && (strcasecmp(text, "maturitysteps") != 0))
480 :			handle.logFileUnexpected(LOGFAIL, "maturitysteps", text);
481 :
482 :			while (isdigit(infile.peek()) && !infile.eof()) {
483 :			infile >> tmpint >> ws;
484 :			maturitystep.resize(1, tmpint);
485 :			i++;
486 :			}
487 :
488 :			for (i = 0; i < maturitystep.Size(); i++)
489 :			if (maturitystep[i] < 1 || maturitystep[i] > TimeInfo->numSteps())
490 :			handle.logFileMessage(LOGFAIL, "invalid maturity step", maturitystep[i]);
491 :
492 :			infile >> ws;
493 :			if (!infile.eof()) {
494 :			infile >> text >> ws;
495 :			handle.logFileUnexpected(LOGFAIL, "<end of file>", text);
496 :			}
497 :
498 :			handle.logMessage(LOGMESSAGE, "Read maturity data file");
499 :			keeper->clearLast();
500 :			}
501 :
502 :			void MaturityC::setStock(StockPtrVector& stockvec) {
503 :			Maturity::setStock(stockvec);
504 :
505 :			int i;
506 :			minMatureAge = 9999;
507 :			double minlength = 9999.0;
508 :			for (i = 0; i < matureStocks.Size(); i++) {
509 :			minMatureAge = min(matureStocks[i]->minAge(), minMatureAge);
510 :			minlength = min(matureStocks[i]->getLengthGroupDiv()->minLength(), minlength);
511 :			}
512 :			minMatureLength = LgrpDiv->numLengthGroup(minlength);
513 :			if (minMatureAge < minStockAge)
514 :			handle.logMessage(LOGFAIL, "Error in maturity - minimum mature age is less than stock age for stock", this->getName());
515 :			}
516 :
517 :			void MaturityC::Reset(const TimeClass* const TimeInfo) {
518 :			Maturity::Reset(TimeInfo);
519 :
520 :			maturityParameters.Update(TimeInfo);
521 :			if (maturityParameters.didChange(TimeInfo)) {
522 :			if (maturityParameters[1] < LgrpDiv->minLength())
523 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - l50 less than minimum length for stock", this->getName());
524 :			if (maturityParameters[1] > LgrpDiv->maxLength())
525 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - l50 greater than maximum length for stock", this->getName());
526 :
527 :			int age, len;
528 :			for (age = 0; age < preCalcMaturation.Nrow(); age++) {
529 :			for (len = 0; len < LgrpDiv->numLengthGroups(); len++) {
530 :			if ((age + minStockAge >= minMatureAge) && (len >= minMatureLength)) {
531 :			tmpratio = exp(-1.0 * maturityParameters[0] * (LgrpDiv->meanLength(len) - maturityParameters[1])
532 :			- maturityParameters[2] * (age + minStockAge - maturityParameters[3]));
533 :			preCalcMaturation[age][len] = 1.0 / (1.0 + tmpratio);
534 :			if (preCalcMaturation[age][len] < 0.0)
535 :			preCalcMaturation[age][len] = 0.0;
536 :			if (preCalcMaturation[age][len] > 1.0)
537 :			preCalcMaturation[age][len] = 1.0;
538 :
539 :			} else
540 :			preCalcMaturation[age][len] = 0.0;
541 :			}
542 :			}
543 :
544 :			if (handle.getLogLevel() >= LOGMESSAGE)
545 :			handle.logMessage(LOGMESSAGE, "Reset maturity data for stock", this->getName());
546 :			}
547 :			}
548 :
549 :			double MaturityC::calcMaturation(int age, int length, int growth, double weight) {
550 :			if (age >= minStockAge)
551 :			return preCalcMaturation[age - minStockAge][length];
552 :			return 0.0;
553 :			}
554 :
555 :			int MaturityC::isMaturationStep(const TimeClass* const TimeInfo) {
556 :			int i;
557 :			for (i = 0; i < maturitystep.Size(); i++)
558 :			if (maturitystep[i] == TimeInfo->getStep())
559 :			return 1;
560 :			return 0;
561 :			}
562 :
563 :			void MaturityC::Print(ofstream& outfile) const {
564 :			int i;
565 :			Maturity::Print(outfile);
566 :			outfile << "\tPrecalculated maturity:\n";
567 :			preCalcMaturation.Print(outfile);
568 :			outfile << "\tMaturity timesteps:";
569 :			for (i = 0; i < maturitystep.Size(); i++)
570 :			outfile << sep << maturitystep[i];
571 :			outfile << endl;
572 :			}
573 :
574 :			// ********************************************************
575 :			// Functions for MaturityD
576 :			// ********************************************************
577 :			MaturityD::MaturityD(CommentStream& infile, const TimeClass* const TimeInfo,
578 :			Keeper* const keeper, int minage, int numage, const IntVector& tmpareas,
579 :			const char* givenname, const LengthGroupDivision* const lgrpdiv, int numMatConst, const char* refWeightFile)
580 :			: MaturityC(infile, TimeInfo, keeper, minage, numage, tmpareas, givenname, lgrpdiv, numMatConst) {
581 :
582 :			//read information on reference weights.
583 :			ifstream subweightfile;
584 :			subweightfile.open(refWeightFile, ios::in);
585 :			handle.checkIfFailure(subweightfile, refWeightFile);
586 :			handle.Open(refWeightFile);
587 :			CommentStream subweightcomment(subweightfile);
588 :			DoubleMatrix tmpRefW;
589 :			readRefWeights(subweightcomment, tmpRefW);
590 :			handle.Close();
591 :			subweightfile.close();
592 :			subweightfile.clear();
593 :
594 :			//aggregate the reference weight data to be the same length groups
595 :			if (LgrpDiv->meanLength(0) < tmpRefW[0][0] ||
596 :			LgrpDiv->meanLength(LgrpDiv->numLengthGroups() - 1) > tmpRefW[tmpRefW.Nrow() - 1][0])
597 :			handle.logFileMessage(LOGFAIL, "lengths for reference weights must span the range of growth lengths");
598 :
599 :			refWeight.resize(LgrpDiv->numLengthGroups(), 0.0);
600 :			int i, j, pos = 0;
601 :			double tmplen;
602 :			for (j = 0; j < LgrpDiv->numLengthGroups(); j++) {
603 :			tmplen = LgrpDiv->meanLength(j);
604 :			for (i = pos; i < tmpRefW.Nrow() - 1; i++) {
605 :			if (((tmplen > tmpRefW[i][0]) || (isEqual(tmplen, tmpRefW[i][0]))) &&
606 :			((tmplen < tmpRefW[i + 1][0]) || (isEqual(tmplen, tmpRefW[i + 1][0])))) {
607 :
608 :			tmpratio = (tmplen - tmpRefW[i][0]) / (tmpRefW[i + 1][0] - tmpRefW[i][0]);
609 :			refWeight[j] = tmpRefW[i][1] + tmpratio * (tmpRefW[i + 1][1] - tmpRefW[i][1]);
610 :			pos = i;
611 :			}
612 :			}
613 :			}
614 :			}
615 :
616 :			void MaturityD::setStock(StockPtrVector& stockvec) {
617 :			Maturity::setStock(stockvec);
618 :
619 :			int i;
620 :			minMatureAge = 9999;
621 :			double minlength = 9999.0;
622 :			for (i = 0; i < matureStocks.Size(); i++) {
623 :			minMatureAge = min(matureStocks[i]->minAge(), minMatureAge);
624 :			minlength = min(matureStocks[i]->getLengthGroupDiv()->minLength(), minlength);
625 :			}
626 :			minMatureLength = LgrpDiv->numLengthGroup(minlength);
627 :			if (minMatureAge < minStockAge)
628 :			handle.logMessage(LOGFAIL, "Error in maturity - minimum mature age is less than stock age for stock", this->getName());
629 :			}
630 :
631 :			void MaturityD::Reset(const TimeClass* const TimeInfo) {
632 :			Maturity::Reset(TimeInfo);
633 :
634 :			maturityParameters.Update(TimeInfo);
635 :			if (maturityParameters.didChange(TimeInfo)) {
636 :			if (maturityParameters[1] < LgrpDiv->minLength())
637 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - l50 less than minimum length for stock", this->getName());
638 :			if (maturityParameters[1] > LgrpDiv->maxLength())
639 :			handle.logMessage(LOGWARN, "Warning in maturity calculation - l50 greater than maximum length for stock", this->getName());
640 :
641 :			if (handle.getLogLevel() >= LOGMESSAGE)
642 :			handle.logMessage(LOGMESSAGE, "Reset maturity data for stock", this->getName());
643 :			}
644 :			}
645 :
646 :			double MaturityD::calcMaturation(int age, int length, int growth, double weight) {
647 :
648 :			if ((age >= minMatureAge) && (length >= minMatureLength)) {
649 :			double tmpweight, my;
650 :
651 :			if ((length >= refWeight.Size()) || (isZero(refWeight[length])))
652 :			tmpweight = maturityParameters[5];
653 :			else
654 :			tmpweight = weight / refWeight[length];
655 :
656 :			my = exp(-1.0 * maturityParameters[0] * (LgrpDiv->meanLength(length) - maturityParameters[1])
657 :			- maturityParameters[2] * (age + minStockAge - maturityParameters[3])
658 :			- maturityParameters[4] * (tmpweight - maturityParameters[5]));
659 :			tmpratio = 1.0 / (1.0 + my);
660 :			return (min(max(0.0, tmpratio), 1.0));
661 :			}
662 :			return 0.0;
663 :			}
664 :
665 :			int MaturityD::isMaturationStep(const TimeClass* const TimeInfo) {
666 :			int i;
667 :			for (i = 0; i < maturitystep.Size(); i++)
668 :			if (maturitystep[i] == TimeInfo->getStep())
669 :			return 1;
670 :			return 0;
671 :			}
672 :
673 :			void MaturityD::Print(ofstream& outfile) const {
674 :			int i;
675 :			Maturity::Print(outfile);
676 :			outfile << "\tMaturity timesteps:";
677 :			for (i = 0; i < maturitystep.Size(); i++)
678 :			outfile << sep << maturitystep[i];
679 :			outfile << endl;
680 :			}

---

root@forge.cesga.es	ViewVC Help
Powered by ViewVC 1.0.0