Annotation of /trunk/gadget/stockdistribution.cc

Parent Directory | Revision Log

---

Revision 1 - (view) (download)

1 :	agomez	1	#include "stockdistribution.h"
2 :			#include "readfunc.h"
3 :			#include "readword.h"
4 :			#include "readaggregation.h"
5 :			#include "errorhandler.h"
6 :			#include "areatime.h"
7 :			#include "fleet.h"
8 :			#include "stock.h"
9 :			#include "mathfunc.h"
10 :			#include "stockprey.h"
11 :			#include "gadget.h"
12 :			#include "global.h"
13 :
14 :			StockDistribution::StockDistribution(CommentStream& infile,
15 :			const AreaClass* const Area, const TimeClass* const TimeInfo,
16 :			double weight, const char* name)
17 :			: Likelihood(STOCKDISTRIBUTIONLIKELIHOOD, weight, name), alptr(0) {
18 :
19 :			int i, j;
20 :			char text[MaxStrLength];
21 :			strncpy(text, "", MaxStrLength);
22 :			int numarea = 0, numage = 0, numlen = 0;
23 :
24 :			char datafilename[MaxStrLength];
25 :			char aggfilename[MaxStrLength];
26 :			strncpy(datafilename, "", MaxStrLength);
27 :			strncpy(aggfilename, "", MaxStrLength);
28 :			ifstream datafile;
29 :			CommentStream subdata(datafile);
30 :
31 :			timeindex = 0;
32 :			yearly = 0;
33 :			functionname = new char[MaxStrLength];
34 :			strncpy(functionname, "", MaxStrLength);
35 :			readWordAndValue(infile, "datafile", datafilename);
36 :			readWordAndValue(infile, "function", functionname);
37 :			if (strcasecmp(functionname, "multinomial") == 0) {
38 :			MN = Multinomial();
39 :			functionnumber = 1;
40 :			} else if ((strcasecmp(functionname, "sumofsquares") == 0) || (strcasecmp(functionname, "stratified") == 0)) {
41 :			functionnumber = 2;
42 :			} else
43 :			handle.logFileMessage(LOGFAIL, "\nError in stockdistribution - unrecognised function", functionname);
44 :
45 :			infile >> ws;
46 :			char c = infile.peek();
47 :			if ((c == 'a') || (c == 'A')) {
48 :			//we have found either aggregationlevel or areaaggfile ...
49 :			streampos pos = infile.tellg();
50 :
51 :			infile >> text >> ws;
52 :			if ((strcasecmp(text, "aggregation_level") == 0) || (strcasecmp(text, "aggregationlevel") == 0))
53 :			infile >> yearly >> ws;
54 :			else if (strcasecmp(text, "areaaggfile") == 0)
55 :			infile.seekg(pos);
56 :			else
57 :			handle.logFileUnexpected(LOGFAIL, "areaaggfile", text);
58 :
59 :			//JMB - peek at the next char
60 :			c = infile.peek();
61 :
62 :			if (yearly != 0 && yearly != 1)
63 :			handle.logFileMessage(LOGFAIL, "\nError in stockdistribution - aggregationlevel must be 0 or 1");
64 :			}
65 :
66 :			//JMB - changed to make the reading of overconsumption optional
67 :			if ((c == 'o') || (c == 'O')) {
68 :			readWordAndVariable(infile, "overconsumption", overconsumption);
69 :			infile >> ws;
70 :			c = infile.peek();
71 :			} else
72 :			overconsumption = 0;
73 :
74 :			if (overconsumption != 0 && overconsumption != 1)
75 :			handle.logFileMessage(LOGFAIL, "\nError in stockdistribution - overconsumption must be 0 or 1");
76 :
77 :			//JMB - changed to make the reading of minimum probability optional
78 :			if ((c == 'm') || (c == 'M'))
79 :			readWordAndVariable(infile, "minimumprobability", epsilon);
80 :			else if ((c == 'e') || (c == 'E'))
81 :			readWordAndVariable(infile, "epsilon", epsilon);
82 :			else
83 :			epsilon = 10.0;
84 :
85 :			if (epsilon < verysmall) {
86 :			handle.logFileMessage(LOGWARN, "epsilon should be a positive integer - set to default value 10");
87 :			epsilon = 10.0;
88 :			}
89 :
90 :			//read in area aggregation from file
91 :			readWordAndValue(infile, "areaaggfile", aggfilename);
92 :			datafile.open(aggfilename, ios::in);
93 :			handle.checkIfFailure(datafile, aggfilename);
94 :			handle.Open(aggfilename);
95 :			numarea = readAggregation(subdata, areas, areaindex);
96 :			handle.Close();
97 :			datafile.close();
98 :			datafile.clear();
99 :
100 :			//read in age aggregation from file
101 :			readWordAndValue(infile, "ageaggfile", aggfilename);
102 :			datafile.open(aggfilename, ios::in);
103 :			handle.checkIfFailure(datafile, aggfilename);
104 :			handle.Open(aggfilename);
105 :			numage = readAggregation(subdata, ages, ageindex);
106 :			handle.Close();
107 :			datafile.close();
108 :			datafile.clear();
109 :
110 :			//read in length aggregation from file
111 :			readWordAndValue(infile, "lenaggfile", aggfilename);
112 :			datafile.open(aggfilename, ios::in);
113 :			handle.checkIfFailure(datafile, aggfilename);
114 :			handle.Open(aggfilename);
115 :			numlen = readLengthAggregation(subdata, lengths, lenindex);
116 :			handle.Close();
117 :			datafile.close();
118 :			datafile.clear();
119 :
120 :			//Must change from outer areas to inner areas.
121 :			for (i = 0; i < areas.Nrow(); i++)
122 :			for (j = 0; j < areas.Ncol(i); j++)
123 :			areas[i][j] = Area->getInnerArea(areas[i][j]);
124 :
125 :			//Must create the length group division
126 :			LgrpDiv = new LengthGroupDivision(lengths);
127 :			if (LgrpDiv->Error())
128 :			handle.logMessage(LOGFAIL, "Error in stockdistribution - failed to create length group");
129 :
130 :			//read in the fleetnames
131 :			i = 0;
132 :			infile >> text >> ws;
133 :			if (strcasecmp(text, "fleetnames") != 0)
134 :			handle.logFileUnexpected(LOGFAIL, "fleetnames", text);
135 :			infile >> text >> ws;
136 :			while (!infile.eof() && (strcasecmp(text, "stocknames") != 0)) {
137 :			fleetnames.resize(new char[strlen(text) + 1]);
138 :			strcpy(fleetnames[i++], text);
139 :			infile >> text >> ws;
140 :			}
141 :			if (fleetnames.Size() == 0)
142 :			handle.logFileMessage(LOGFAIL, "\nError in stockdistribution - failed to read fleets");
143 :			handle.logMessage(LOGMESSAGE, "Read fleet data - number of fleets", fleetnames.Size());
144 :
145 :			//read in the stocknames
146 :			i = 0;
147 :			if (strcasecmp(text, "stocknames") != 0)
148 :			handle.logFileUnexpected(LOGFAIL, "stocknames", text);
149 :			infile >> text;
150 :			while (!infile.eof() && (strcasecmp(text, "[component]") != 0)) {
151 :			infile >> ws;
152 :			stocknames.resize(new char[strlen(text) + 1]);
153 :			strcpy(stocknames[i++], text);
154 :			infile >> text;
155 :			}
156 :			if (stocknames.Size() == 0)
157 :			handle.logFileMessage(LOGFAIL, "\nError in stockdistribution - failed to read stocks");
158 :			handle.logMessage(LOGMESSAGE, "Read stock data - number of stocks", stocknames.Size());
159 :
160 :			//We have now read in all the data from the main likelihood file
161 :			//But we have to read in the statistics data from datafilename
162 :			datafile.open(datafilename, ios::in);
163 :			handle.checkIfFailure(datafile, datafilename);
164 :			handle.Open(datafilename);
165 :			readStockData(subdata, TimeInfo, numarea, numage, numlen);
166 :			handle.Close();
167 :			datafile.close();
168 :			datafile.clear();
169 :
170 :			switch (functionnumber) {
171 :			case 2:
172 :			for (i = 0; i < numarea; i++) {
173 :			modelYearData.resize(new DoubleMatrix(stocknames.Size(), (numage * numlen), 0.0));
174 :			obsYearData.resize(new DoubleMatrix(stocknames.Size(), (numage * numlen), 0.0));
175 :			}
176 :			break;
177 :			case 1:
178 :			if (yearly)
179 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - yearly aggregation is ignored for function", functionname);
180 :			yearly = 0;
181 :			break;
182 :			default:
183 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - unrecognised function", functionname);
184 :			break;
185 :			}
186 :			}
187 :
188 :			void StockDistribution::readStockData(CommentStream& infile,
189 :			const TimeClass* TimeInfo, int numarea, int numage, int numlen) {
190 :
191 :			double tmpnumber;
192 :			char tmparea[MaxStrLength], tmpstock[MaxStrLength];
193 :			char tmpage[MaxStrLength], tmplen[MaxStrLength];
194 :			strncpy(tmparea, "", MaxStrLength);
195 :			strncpy(tmpstock, "", MaxStrLength);
196 :			strncpy(tmpage, "", MaxStrLength);
197 :			strncpy(tmplen, "", MaxStrLength);
198 :			int i, j, year, step, count, reject;
199 :			int keepdata, timeid, stockid, ageid, areaid, lenid;
200 :			int numstock = stocknames.Size();
201 :
202 :			//Check the number of columns in the inputfile
203 :			infile >> ws;
204 :			if (countColumns(infile) != 7)
205 :			handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 7");
206 :
207 :			year = step = count = reject = 0;
208 :			while (!infile.eof()) {
209 :			keepdata = 1;
210 :			infile >> year >> step >> tmparea >> tmpstock >> tmpage >> tmplen >> tmpnumber >> ws;
211 :
212 :			//crude check to see if something has gone wrong and avoid infinite loops
213 :			if (strlen(tmparea) == 0)
214 :			handle.logFileMessage(LOGFAIL, "failed to read data from file");
215 :
216 :			//if tmpstock is in stocknames find stockid, else dont keep the data
217 :			stockid = -1;
218 :			for (i = 0; i < numstock; i++)
219 :			if (strcasecmp(stocknames[i], tmpstock) == 0)
220 :			stockid = i;
221 :
222 :			if (stockid == -1)
223 :			keepdata = 0;
224 :
225 :			//if tmparea is in areaindex find areaid, else dont keep the data
226 :			areaid = -1;
227 :			for (i = 0; i < numarea; i++)
228 :			if (strcasecmp(areaindex[i], tmparea) == 0)
229 :			areaid = i;
230 :
231 :			if (areaid == -1)
232 :			keepdata = 0;
233 :
234 :			//if tmpage is in ageindex find ageid, else dont keep the data
235 :			ageid = -1;
236 :			for (i = 0; i < numage; i++)
237 :			if (strcasecmp(ageindex[i], tmpage) == 0)
238 :			ageid = i;
239 :
240 :			if (ageid == -1)
241 :			keepdata = 0;
242 :
243 :			//if tmplen is in lenindex find lenid, else dont keep the data
244 :			lenid = -1;
245 :			for (i = 0; i < numlen; i++)
246 :			if (strcasecmp(lenindex[i], tmplen) == 0)
247 :			lenid = i;
248 :
249 :			if (lenid == -1)
250 :			keepdata = 0;
251 :
252 :			//check if the year and step are in the simulation
253 :			timeid = -1;
254 :			if ((TimeInfo->isWithinPeriod(year, step)) && (keepdata == 1)) {
255 :			//if this is a new timestep, resize to store the data
256 :			for (i = 0; i < Years.Size(); i++)
257 :			if ((Years[i] == year) && (Steps[i] == step))
258 :			timeid = i;
259 :
260 :			if (timeid == -1) {
261 :			Years.resize(1, year);
262 :			Steps.resize(1, step);
263 :			timeid = (Years.Size() - 1);
264 :
265 :			obsDistribution.resize();
266 :			modelDistribution.resize();
267 :			likelihoodValues.AddRows(1, numarea, 0.0);
268 :			for (i = 0; i < numarea; i++) {
269 :			obsDistribution[timeid].resize(new DoubleMatrix(numstock, (numage * numlen), 0.0));
270 :			modelDistribution[timeid].resize(new DoubleMatrix(numstock, (numage * numlen), 0.0));
271 :			}
272 :			}
273 :
274 :			} else
275 :			keepdata = 0;
276 :
277 :			if (keepdata == 1) {
278 :			//stock distribution data is required, so store it
279 :			count++;
280 :			i = ageid + (numage * lenid);
281 :			//JMB - this should be stored as [time][area][stock][age][length]
282 :			(*obsDistribution[timeid][areaid])[stockid][i] = tmpnumber;
283 :			} else
284 :			reject++; //count number of rejected data points read from file
285 :			}
286 :
287 :			AAT.addActions(Years, Steps, TimeInfo);
288 :			if (count == 0)
289 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - found no data in the data file for", this->getName());
290 :			if (reject != 0)
291 :			handle.logMessage(LOGMESSAGE, "Discarded invalid stockdistribution data - number of invalid entries", reject);
292 :			handle.logMessage(LOGMESSAGE, "Read stockdistribution data file - number of entries", count);
293 :			}
294 :
295 :			StockDistribution::~StockDistribution() {
296 :			int i, j;
297 :			for (i = 0; i < stocknames.Size(); i++) {
298 :			delete[] stocknames[i];
299 :			delete aggregator[i];
300 :			}
301 :			delete[] aggregator;
302 :			delete[] functionname;
303 :			delete LgrpDiv;
304 :
305 :			for (i = 0; i < fleetnames.Size(); i++)
306 :			delete[] fleetnames[i];
307 :			for (i = 0; i < areaindex.Size(); i++)
308 :			delete[] areaindex[i];
309 :			for (i = 0; i < ageindex.Size(); i++)
310 :			delete[] ageindex[i];
311 :			for (i = 0; i < lenindex.Size(); i++)
312 :			delete[] lenindex[i];
313 :			for (i = 0; i < modelYearData.Size(); i++) {
314 :			delete modelYearData[i];
315 :			delete obsYearData[i];
316 :			}
317 :			for (i = 0; i < obsDistribution.Nrow(); i++) {
318 :			for (j = 0; j < obsDistribution.Ncol(i); j++) {
319 :			delete obsDistribution[i][j];
320 :			delete modelDistribution[i][j];
321 :			}
322 :			}
323 :			}
324 :
325 :			void StockDistribution::Reset(const Keeper* const keeper) {
326 :			Likelihood::Reset(keeper);
327 :			if (isZero(weight))
328 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - zero weight for", this->getName());
329 :			int i, j;
330 :			for (i = 0; i < modelDistribution.Nrow(); i++)
331 :			for (j = 0; j < modelDistribution.Ncol(i); j++)
332 :			(*modelDistribution[i][j]).setToZero();
333 :			if (yearly)
334 :			for (i = 0; i < modelYearData.Size(); i++) {
335 :			(*modelYearData[i]).setToZero();
336 :			(*obsYearData[i]).setToZero();
337 :			}
338 :
339 :			switch (functionnumber) {
340 :			case 1:
341 :			MN.setValue(epsilon);
342 :			break;
343 :			case 2:
344 :			break;
345 :			default:
346 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - unrecognised function", functionname);
347 :			break;
348 :			}
349 :			if (handle.getLogLevel() >= LOGMESSAGE)
350 :			handle.logMessage(LOGMESSAGE, "Reset stockdistribution component", this->getName());
351 :			}
352 :
353 :			void StockDistribution::Print(ofstream& outfile) const {
354 :			int i;
355 :			outfile << "\nStock Distribution " << this->getName() << " - likelihood value " << likelihood
356 :			<< "\n\tFunction " << functionname << "\n\tStock names:";
357 :			for (i = 0; i < stocknames.Size(); i++)
358 :			outfile << sep << stocknames[i];
359 :			outfile << "\n\tFleet names:";
360 :			for (i = 0; i < fleetnames.Size(); i++)
361 :			outfile << sep << fleetnames[i];
362 :			outfile << endl;
363 :			for (i = 0; i < stocknames.Size(); i++)
364 :			aggregator[i]->Print(outfile);
365 :			outfile.flush();
366 :			}
367 :
368 :			void StockDistribution::setFleetsAndStocks(FleetPtrVector& Fleets, StockPtrVector& Stocks) {
369 :			int s, i, j, k, found, minage, maxage;
370 :			FleetPtrVector fleets;
371 :			StockPtrVector stocks;
372 :			StockPtrVector checkstocks;
373 :			aggregator = new FleetPreyAggregator*[stocknames.Size()];
374 :
375 :			for (i = 0; i < fleetnames.Size(); i++) {
376 :			found = 0;
377 :			for (j = 0; j < Fleets.Size(); j++) {
378 :			if (strcasecmp(fleetnames[i], Fleets[j]->getName()) == 0) {
379 :			found++;
380 :			fleets.resize(Fleets[j]);
381 :			}
382 :			}
383 :			if (found == 0)
384 :			handle.logMessage(LOGFAIL, "Error in stockdistribution - unrecognised fleet", fleetnames[i]);
385 :			}
386 :
387 :			for (i = 0; i < fleets.Size(); i++)
388 :			for (j = 0; j < fleets.Size(); j++)
389 :			if ((strcasecmp(fleets[i]->getName(), fleets[j]->getName()) == 0) && (i != j))
390 :			handle.logMessage(LOGFAIL, "Error in stockdistribution - repeated fleet", fleets[i]->getName());
391 :
392 :			//check fleet areas
393 :			if (handle.getLogLevel() >= LOGWARN) {
394 :			for (j = 0; j < areas.Nrow(); j++) {
395 :			found = 0;
396 :			for (i = 0; i < fleets.Size(); i++)
397 :			for (k = 0; k < areas.Ncol(j); k++)
398 :			if (fleets[i]->isInArea(areas[j][k]))
399 :			found++;
400 :			if (found == 0)
401 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - fleet not defined on all areas");
402 :			}
403 :			}
404 :
405 :			for (s = 0; s < stocknames.Size(); s++) {
406 :			found = 0;
407 :			stocks.Reset();
408 :			for (j = 0; j < Stocks.Size(); j++) {
409 :			if (Stocks[j]->isEaten()) {
410 :			if (strcasecmp(stocknames[s], Stocks[j]->getName()) == 0) {
411 :			found++;
412 :			stocks.resize(Stocks[j]);
413 :			checkstocks.resize(Stocks[j]);
414 :			}
415 :			}
416 :			}
417 :			if (found == 0)
418 :			handle.logMessage(LOGFAIL, "Error in stockdistribution - unrecognised stock", stocknames[i]);
419 :
420 :			aggregator[s] = new FleetPreyAggregator(fleets, stocks, LgrpDiv, areas, ages, overconsumption);
421 :			}
422 :
423 :			for (i = 0; i < checkstocks.Size(); i++)
424 :			for (j = 0; j < checkstocks.Size(); j++)
425 :			if ((strcasecmp(checkstocks[i]->getName(), checkstocks[j]->getName()) == 0) && (i != j))
426 :			handle.logMessage(LOGFAIL, "Error in stockdistribution - repeated stock", checkstocks[i]->getName());
427 :
428 :			//check areas, ages and lengths
429 :			if (handle.getLogLevel() >= LOGWARN) {
430 :			for (j = 0; j < areas.Nrow(); j++) {
431 :			found = 0;
432 :			for (i = 0; i < checkstocks.Size(); i++)
433 :			for (k = 0; k < areas.Ncol(j); k++)
434 :			if (checkstocks[i]->isInArea(areas[j][k]))
435 :			found++;
436 :			if (found == 0)
437 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - stock not defined on all areas");
438 :			}
439 :
440 :			minage = 9999;
441 :			maxage = 0;
442 :			for (i = 0; i < ages.Nrow(); i++) {
443 :			for (j = 0; j < ages.Ncol(i); j++) {
444 :			minage = min(ages[i][j], minage);
445 :			maxage = max(ages[i][j], maxage);
446 :			}
447 :			}
448 :
449 :			found = 0;
450 :			for (i = 0; i < checkstocks.Size(); i++)
451 :			if (minage >= checkstocks[i]->minAge())
452 :			found++;
453 :			if (found == 0)
454 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - minimum age less than stock age");
455 :
456 :			found = 0;
457 :			for (i = 0; i < checkstocks.Size(); i++)
458 :			if (maxage <= checkstocks[i]->maxAge())
459 :			found++;
460 :			if (found == 0)
461 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - maximum age greater than stock age");
462 :
463 :			found = 0;
464 :			for (i = 0; i < checkstocks.Size(); i++)
465 :			if (LgrpDiv->maxLength(0) > checkstocks[i]->getLengthGroupDiv()->minLength())
466 :			found++;
467 :			if (found == 0)
468 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - minimum length group less than stock length");
469 :
470 :			found = 0;
471 :			for (i = 0; i < checkstocks.Size(); i++)
472 :			if (LgrpDiv->minLength(LgrpDiv->numLengthGroups()) < checkstocks[i]->getLengthGroupDiv()->maxLength())
473 :			found++;
474 :			if (found == 0)
475 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - maximum length group greater than stock length");
476 :			}
477 :			}
478 :
479 :			void StockDistribution::addLikelihood(const TimeClass* const TimeInfo) {
480 :
481 :			if ((!(AAT.atCurrentTime(TimeInfo))) || (isZero(weight)))
482 :			return;
483 :
484 :			if ((handle.getLogLevel() >= LOGMESSAGE) && ((!yearly) || (TimeInfo->getStep() == TimeInfo->numSteps())))
485 :			handle.logMessage(LOGMESSAGE, "Calculating likelihood score for stockdistribution component", this->getName());
486 :
487 :			int i;
488 :			timeindex = -1;
489 :			for (i = 0; i < Years.Size(); i++)
490 :			if ((Years[i] == TimeInfo->getYear()) && (Steps[i] == TimeInfo->getStep()))
491 :			timeindex = i;
492 :			if (timeindex == -1)
493 :			handle.logMessage(LOGFAIL, "Error in stockdistribution - invalid timestep");
494 :
495 :			for (i = 0; i < stocknames.Size(); i++) {
496 :			aggregator[i]->Sum();
497 :			if ((handle.getLogLevel() >= LOGWARN) && (aggregator[i]->checkCatchData()))
498 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - zero catch found");
499 :			}
500 :
501 :			double l = 0.0;
502 :			switch (functionnumber) {
503 :			case 1:
504 :			l = calcLikMultinomial();
505 :			break;
506 :			case 2:
507 :			l = calcLikSumSquares(TimeInfo);
508 :			break;
509 :			default:
510 :			handle.logMessage(LOGWARN, "Warning in stockdistribution - unrecognised function", functionname);
511 :			break;
512 :			}
513 :
514 :			if ((!yearly) || (TimeInfo->getStep() == TimeInfo->numSteps())) {
515 :			likelihood += l;
516 :			if (handle.getLogLevel() >= LOGMESSAGE)
517 :			handle.logMessage(LOGMESSAGE, "The likelihood score for this component on this timestep is", l);
518 :			}
519 :			}
520 :
521 :			//The code here is probably unnessecarily complicated because
522 :			//there is always only one length group with this class.
523 :			double StockDistribution::calcLikMultinomial() {
524 :			int age, len, area, s, i;
525 :			//JMB - number of age and length groups is constant for all stocks
526 :			int numage = ages.Nrow();
527 :			int numlen = LgrpDiv->numLengthGroups();
528 :			int numstock = stocknames.Size();
529 :			DoubleVector obsdata(numstock, 0.0);
530 :			DoubleVector moddata(numstock, 0.0);
531 :
532 :			MN.Reset();
533 :			//the object MN does most of the work, accumulating likelihood
534 :			for (area = 0; area < areas.Nrow(); area++) {
535 :			likelihoodValues[timeindex][area] = 0.0;
536 :			for (s = 0; s < numstock; s++) {
537 :			alptr = &aggregator[s]->getSum();
538 :			for (age = (*alptr)[area].minAge(); age <= (*alptr)[area].maxAge(); age++)
539 :			for (len = (*alptr)[area].minLength(age); len < (*alptr)[area].maxLength(age); len++)
540 :			(*modelDistribution[timeindex][area])[s][age + (numage * len)] = ((*alptr)[area][age][len]).N;
541 :			}
542 :
543 :			for (i = 0; i < (numage * numlen); i++) {
544 :			for (s = 0; s < numstock; s++) {
545 :			obsdata[s] = (*obsDistribution[timeindex][area])[s][i];
546 :			moddata[s] = (*modelDistribution[timeindex][area])[s][i];
547 :			}
548 :			likelihoodValues[timeindex][area] += MN.calcLogLikelihood(obsdata, moddata);
549 :			}
550 :			}
551 :			return MN.getLogLikelihood();
552 :			}
553 :
554 :			double StockDistribution::calcLikSumSquares(const TimeClass* const TimeInfo) {
555 :			double temp, totalmodel, totaldata, totallikelihood;
556 :			int age, len, area, s, i;
557 :			//JMB - number of age and length groups is constant for all stocks
558 :			int numage = ages.Nrow();
559 :			int numlen = LgrpDiv->numLengthGroups();
560 :			int numstock = stocknames.Size();
561 :
562 :			totallikelihood = 0.0;
563 :			for (area = 0; area < areas.Nrow(); area++) {
564 :			likelihoodValues[timeindex][area] = 0.0;
565 :
566 :			for (s = 0; s < numstock; s++) {
567 :			alptr = &aggregator[s]->getSum();
568 :			for (age = (*alptr)[area].minAge(); age <= (*alptr)[area].maxAge(); age++)
569 :			for (len = (*alptr)[area].minLength(age); len < (*alptr)[area].maxLength(age); len++)
570 :			(*modelDistribution[timeindex][area])[s][age + (numage * len)] = ((*alptr)[area][age][len]).N;
571 :			}
572 :
573 :			if (!yearly) { //calculate likelihood on all steps
574 :			for (i = 0; i < (numage * numlen); i++) {
575 :			totalmodel = 0.0;
576 :			totaldata = 0.0;
577 :			for (s = 0; s < numstock; s++) {
578 :			totalmodel += (*modelDistribution[timeindex][area])[s][i];
579 :			totaldata += (*obsDistribution[timeindex][area])[s][i];
580 :			}
581 :
582 :			if (!(isZero(totalmodel)))
583 :			totalmodel = 1.0 / totalmodel;
584 :			if (!(isZero(totaldata)))
585 :			totaldata = 1.0 / totaldata;
586 :
587 :			for (s = 0; s < numstock; s++) {
588 :			temp = (((*obsDistribution[timeindex][area])[s][i] * totaldata)
589 :			- ((*modelDistribution[timeindex][area])[s][i] * totalmodel));
590 :			likelihoodValues[timeindex][area] += (temp * temp);
591 :			}
592 :			}
593 :			totallikelihood += likelihoodValues[timeindex][area];
594 :
595 :			} else { //calculate likelihood on year basis
596 :
597 :			if (TimeInfo->getStep() == 1) { //start of a new year
598 :			(*modelYearData[area]).setToZero();
599 :			(*obsYearData[area]).setToZero();
600 :			}
601 :
602 :			for (s = 0; s < numstock; s++) {
603 :			alptr = &aggregator[s]->getSum();
604 :			for (age = (*alptr)[area].minAge(); age <= (*alptr)[area].maxAge(); age++) {
605 :			for (len = (*alptr)[area].minLength(age); len < (*alptr)[area].maxLength(age); len++) {
606 :			i = age + (numage * len);
607 :			(*modelYearData[area])[s][i] += (*modelDistribution[timeindex][area])[s][i];
608 :			(*obsYearData[area])[s][i] += (*obsDistribution[timeindex][area])[s][i];
609 :			}
610 :			}
611 :			}
612 :
613 :			if (TimeInfo->getStep() < TimeInfo->numSteps())
614 :			likelihoodValues[timeindex][area] = 0.0;
615 :			else { //last step in year, so need to calc likelihood contribution
616 :			for (i = 0; i < (numage * numlen); i++) {
617 :			totalmodel = 0.0;
618 :			totaldata = 0.0;
619 :			for (s = 0; s < numstock; s++) {
620 :			totalmodel += (*modelYearData[area])[s][i];
621 :			totaldata += (*obsYearData[area])[s][i];
622 :			}
623 :
624 :			if (!(isZero(totalmodel)))
625 :			totalmodel = 1.0 / totalmodel;
626 :			if (!(isZero(totaldata)))
627 :			totaldata = 1.0 / totaldata;
628 :
629 :			for (s = 0; s < numstock; s++) {
630 :			temp = (((*obsYearData[area])[s][i] * totaldata)
631 :			- ((*modelYearData[area])[s][i] * totalmodel));
632 :			likelihoodValues[timeindex][area] += (temp * temp);
633 :			}
634 :			}
635 :			totallikelihood += likelihoodValues[timeindex][area];
636 :			}
637 :			}
638 :			}
639 :			return totallikelihood;
640 :			}
641 :
642 :			void StockDistribution::printLikelihood(ofstream& outfile, const TimeClass* const TimeInfo) {
643 :
644 :			if (!AAT.atCurrentTime(TimeInfo))
645 :			return;
646 :
647 :			int area, s, i, age, len;
648 :			int numage = ages.Nrow();
649 :			int numlen = LgrpDiv->numLengthGroups();
650 :
651 :			timeindex = -1;
652 :			for (i = 0; i < Years.Size(); i++)
653 :			if ((Years[i] == TimeInfo->getYear()) && (Steps[i] == TimeInfo->getStep()))
654 :			timeindex = i;
655 :			if (timeindex == -1)
656 :			handle.logMessage(LOGFAIL, "Error in stockdistribution - invalid timestep");
657 :
658 :			for (area = 0; area < areas.Nrow(); area++) {
659 :			for (s = 0; s < stocknames.Size(); s++) {
660 :			for (i = 0; i < (numage * numlen); i++) {
661 :			// need to calculate the age and length index from i
662 :			// i = ageid + (numage * lenid);
663 :			age = i % numage;
664 :			len = (i - age) / numage;
665 :
666 :			outfile << setw(lowwidth) << Years[timeindex] << sep << setw(lowwidth)
667 :			<< Steps[timeindex] << sep << setw(printwidth) << areaindex[area] << sep
668 :			<< setw(printwidth) << stocknames[s] << sep << setw(printwidth)
669 :			<< ageindex[age] << sep << setw(printwidth) << lenindex[len]
670 :			<< sep << setprecision(largeprecision) << setw(largewidth);
671 :
672 :			//JMB crude filter to remove the 'silly' values from the output
673 :			if ((*modelDistribution[timeindex][area])[s][i] < rathersmall)
674 :			outfile << 0 << endl;
675 :			else
676 :			outfile << (*modelDistribution[timeindex][area])[s][i] << endl;
677 :			}
678 :			}
679 :			}
680 :			}
681 :
682 :			void StockDistribution::printSummary(ofstream& outfile) {
683 :			int year, area;
684 :
685 :			for (year = 0; year < likelihoodValues.Nrow(); year++) {
686 :			for (area = 0; area < likelihoodValues.Ncol(year); area++) {
687 :			if (!yearly) {
688 :			outfile << setw(lowwidth) << Years[year] << sep << setw(lowwidth)
689 :			<< Steps[year] << sep << setw(printwidth) << areaindex[area] << sep
690 :			<< setw(largewidth) << this->getName() << sep << setprecision(smallprecision)
691 :			<< setw(smallwidth) << weight << sep << setprecision(largeprecision)
692 :			<< setw(largewidth) << likelihoodValues[year][area] << endl;
693 :			} else {
694 :			if (isZero(likelihoodValues[year][area])) {
695 :			// assume that this isnt the last step for that year and ignore
696 :			} else {
697 :			outfile << setw(lowwidth) << Years[year] << " all "
698 :			<< setw(printwidth) << areaindex[area] << sep
699 :			<< setw(largewidth) << this->getName() << sep << setprecision(smallprecision)
700 :			<< setw(smallwidth) << weight << sep << setprecision(largeprecision)
701 :			<< setw(largewidth) << likelihoodValues[year][area] << endl;
702 :			}
703 :			}
704 :			}
705 :			}
706 :			outfile.flush();
707 :			}

---

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