Annotation of /trunk/gadget/stomachcontent.cc

Parent Directory | Revision Log

---

Revision 1 - (view) (download)

1 :	agomez	1	#include "stomachcontent.h"
2 :			#include "areatime.h"
3 :			#include "predator.h"
4 :			#include "stockpredator.h"
5 :			#include "prey.h"
6 :			#include "predatoraggregator.h"
7 :			#include "readfunc.h"
8 :			#include "readword.h"
9 :			#include "readaggregation.h"
10 :			#include "multinomial.h"
11 :			#include "errorhandler.h"
12 :			#include "gadget.h"
13 :			#include "global.h"
14 :
15 :			// ********************************************************
16 :			// Functions for main likelihood component StomachContent
17 :			// ********************************************************
18 :			StomachContent::StomachContent(CommentStream& infile,
19 :			const AreaClass* const Area, const TimeClass* const TimeInfo,
20 :			Keeper* const keeper, double weight, const char* name)
21 :			: Likelihood(STOMACHCONTENTLIKELIHOOD, weight, name) {
22 :
23 :			functionname = new char[MaxStrLength];
24 :			strncpy(functionname, "", MaxStrLength);
25 :			readWordAndValue(infile, "function", functionname);
26 :
27 :			char datafilename[MaxStrLength];
28 :			strncpy(datafilename, "", MaxStrLength);
29 :			readWordAndValue(infile, "datafile", datafilename);
30 :
31 :			if (strcasecmp(functionname, "scnumbers") == 0) {
32 :			StomCont = new SCNumbers(infile, Area, TimeInfo, keeper, datafilename, this->getName());
33 :
34 :			} else if (strcasecmp(functionname, "scratios") == 0) {
35 :			char numfilename[MaxStrLength];
36 :			strncpy(numfilename, "", MaxStrLength);
37 :			readWordAndValue(infile, "numberfile", numfilename);
38 :			StomCont = new SCRatios(infile, Area, TimeInfo, keeper, datafilename, numfilename, this->getName());
39 :
40 :			} else if (strcasecmp(functionname, "scamounts") == 0) {
41 :			char numfilename[MaxStrLength];
42 :			strncpy(numfilename, "", MaxStrLength);
43 :			readWordAndValue(infile, "numberfile", numfilename);
44 :			StomCont = new SCAmounts(infile, Area, TimeInfo, keeper, datafilename, numfilename, this->getName());
45 :
46 :			} else if (strcasecmp(functionname, "scsimple") == 0) {
47 :			StomCont = new SCSimple(infile, Area, TimeInfo, keeper, datafilename, this->getName());
48 :
49 :			} else {
50 :			handle.logFileMessage(LOGFAIL, "\nError in stomachcontent - unrecognised function", functionname);
51 :			}
52 :			}
53 :
54 :			StomachContent::~StomachContent() {
55 :			delete StomCont;
56 :			delete[] functionname;
57 :			}
58 :
59 :			void StomachContent::addLikelihood(const TimeClass* const TimeInfo) {
60 :			if (isZero(weight))
61 :			return;
62 :			likelihood += StomCont->calcLikelihood(TimeInfo);
63 :			}
64 :
65 :			void StomachContent::Reset(const Keeper* const keeper) {
66 :			Likelihood::Reset(keeper);
67 :			if (isZero(weight))
68 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - zero weight for", this->getName());
69 :			StomCont->Reset();
70 :			}
71 :
72 :			void StomachContent::Print(ofstream& outfile) const {
73 :			outfile << "\nStomach Content " << this->getName() << " - likelihood value " << likelihood
74 :			<< "\n\tFunction " << functionname << endl;
75 :			StomCont->Print(outfile);
76 :			}
77 :
78 :			// ********************************************************
79 :			// Functions for base SC
80 :			// ********************************************************
81 :			SC::SC(CommentStream& infile, const AreaClass* const Area, const TimeClass* const TimeInfo,
82 :			Keeper* const keeper, const char* datafilename, const char* givenname)
83 :			: HasName(givenname), aggregator(0), preyLgrpDiv(0), predLgrpDiv(0), dptr(0) {
84 :
85 :			int i, j;
86 :			char text[MaxStrLength];
87 :			strncpy(text, "", MaxStrLength);
88 :			int numpred = 0;
89 :			int numarea = 0;
90 :
91 :			timeindex = 0;
92 :			usepredages = 0;
93 :
94 :			char aggfilename[MaxStrLength];
95 :			strncpy(aggfilename, "", MaxStrLength);
96 :
97 :			ifstream datafile;
98 :			CommentStream subdata(datafile);
99 :
100 :			//JMB - changed to make the reading of minimum probability optional
101 :			infile >> ws;
102 :			char c = infile.peek();
103 :			if ((c == 'm') || (c == 'M'))
104 :			readWordAndVariable(infile, "minimumprobability", epsilon);
105 :			else if ((c == 'e') || (c == 'E'))
106 :			readWordAndVariable(infile, "epsilon", epsilon);
107 :			else
108 :			epsilon = 10.0;
109 :
110 :			if (epsilon < verysmall) {
111 :			handle.logFileMessage(LOGWARN, "epsilon should be a positive integer - set to default value 10");
112 :			epsilon = 10.0;
113 :			}
114 :
115 :			//read in area aggregation from file
116 :			readWordAndValue(infile, "areaaggfile", aggfilename);
117 :			datafile.open(aggfilename, ios::in);
118 :			handle.checkIfFailure(datafile, aggfilename);
119 :			handle.Open(aggfilename);
120 :			numarea = readAggregation(subdata, areas, areaindex);
121 :			handle.Close();
122 :			datafile.close();
123 :			datafile.clear();
124 :
125 :			//Must change from outer areas to inner areas.
126 :			for (i = 0; i < areas.Nrow(); i++)
127 :			for (j = 0; j < areas.Ncol(i); j++)
128 :			areas[i][j] = Area->getInnerArea(areas[i][j]);
129 :
130 :			//read in the predators
131 :			i = 0;
132 :			infile >> text >> ws;
133 :			if ((strcasecmp(text, "predators") != 0) && (strcasecmp(text, "predatornames") != 0))
134 :			handle.logFileUnexpected(LOGFAIL, "predatornames", text);
135 :			infile >> text >> ws;
136 :			while (!infile.eof() && ((strcasecmp(text, "predatorlengths") != 0)
137 :			&& (strcasecmp(text, "predatorages") != 0))) {
138 :			predatornames.resize(new char[strlen(text) + 1]);
139 :			strcpy(predatornames[i++], text);
140 :			infile >> text >> ws;
141 :			}
142 :			if (predatornames.Size() == 0)
143 :			handle.logFileMessage(LOGFAIL, "\nError in stomachcontent - failed to read predators");
144 :			handle.logMessage(LOGMESSAGE, "Read predator data - number of predators", predatornames.Size());
145 :
146 :			if (strcasecmp(text, "predatorlengths") == 0) { //read predator lengths
147 :			usepredages = 0; //predator is length structured
148 :			readWordAndValue(infile, "lenaggfile", aggfilename);
149 :			datafile.open(aggfilename, ios::in);
150 :			handle.checkIfFailure(datafile, aggfilename);
151 :			handle.Open(aggfilename);
152 :			numpred = readLengthAggregation(subdata, predatorlengths, predindex);
153 :			handle.Close();
154 :			datafile.close();
155 :			datafile.clear();
156 :			} else if (strcasecmp(text, "predatorages") == 0) { //read predator ages
157 :			usepredages = 1; //predator is age structured
158 :			readWordAndValue(infile, "ageaggfile", aggfilename);
159 :			datafile.open(aggfilename, ios::in);
160 :			handle.checkIfFailure(datafile, aggfilename);
161 :			handle.Open(aggfilename);
162 :			numpred = readAggregation(subdata, predatorages, predindex);
163 :			handle.Close();
164 :			datafile.close();
165 :			datafile.clear();
166 :			} else
167 :			handle.logFileUnexpected(LOGFAIL, "predatorlengths", text);
168 :
169 :			//read in the preys
170 :			readWordAndValue(infile, "preyaggfile", aggfilename);
171 :			datafile.open(aggfilename, ios::in);
172 :			handle.checkIfFailure(datafile, aggfilename);
173 :			handle.Open(aggfilename);
174 :			i = readPreyAggregation(subdata, preynames, preylengths, digestioncoeff, preyindex, keeper);
175 :			handle.Close();
176 :			datafile.close();
177 :			datafile.clear();
178 :
179 :			//prepare for next likelihood component
180 :			infile >> ws;
181 :			if (!infile.eof()) {
182 :			infile >> text >> ws;
183 :			if (strcasecmp(text, "[component]") != 0)
184 :			handle.logFileUnexpected(LOGFAIL, "[component]", text);
185 :			}
186 :			}
187 :
188 :			void SC::aggregate(int i) {
189 :			aggregator[i]->Sum();
190 :			}
191 :
192 :			double SC::calcLikelihood(const TimeClass* const TimeInfo) {
193 :
194 :			if (!(AAT.atCurrentTime(TimeInfo)))
195 :			return 0.0;
196 :
197 :			int i, a, k, p;
198 :			if (handle.getLogLevel() >= LOGMESSAGE)
199 :			handle.logMessage(LOGMESSAGE, "Calculating likelihood score for stomachcontent component", this->getName());
200 :
201 :			timeindex = -1;
202 :			for (i = 0; i < Years.Size(); i++)
203 :			if ((Years[i] == TimeInfo->getYear()) && (Steps[i] == TimeInfo->getStep()))
204 :			timeindex = i;
205 :			if (timeindex == -1)
206 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - invalid timestep");
207 :
208 :			//Get the consumption from aggregator, indexed the same way as in obsConsumption
209 :			int numprey = 0;
210 :			for (i = 0; i < preyindex.Size(); i++) {
211 :			this->aggregate(i);
212 :			for (a = 0; a < areas.Nrow(); a++) {
213 :			dptr = aggregator[i]->getSum()[a];
214 :			for (k = 0; k < dptr->Nrow(); k++)
215 :			for (p = 0; p < dptr->Ncol(k); p++)
216 :			(*modelConsumption[timeindex][a])[k][numprey + p] = (*dptr)[k][p] * digestion[i][p];
217 :
218 :			}
219 :			numprey += preylengths[i].Size() - 1;
220 :			}
221 :
222 :			//Now calculate likelihood score
223 :			double l = calcLikelihood();
224 :			if (handle.getLogLevel() >= LOGMESSAGE)
225 :			handle.logMessage(LOGMESSAGE, "The likelihood score for this component on this timestep is", l);
226 :			return l;
227 :			}
228 :
229 :			void SC::Reset() {
230 :			//JMB - calculate the digestion coefficiant matrix
231 :			int i, j;
232 :			if (digestion.Nrow() != digestioncoeff.Nrow())
233 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - missing digestion coefficient data");
234 :
235 :			for (i = 0; i < digestion.Nrow(); i++)
236 :			for (j = 0; j < digestion.Ncol(i); j++)
237 :			digestion[i][j] = digestioncoeff[i][0] + digestioncoeff[i][1] *
238 :			pow(preyLgrpDiv[i]->meanLength(j), digestioncoeff[i][2]);
239 :
240 :			for (i = 0; i < modelConsumption.Nrow(); i++)
241 :			for (j = 0; j < modelConsumption.Ncol(i); j++)
242 :			(*modelConsumption[i][j]).setToZero();
243 :
244 :			if (handle.getLogLevel() >= LOGMESSAGE)
245 :			handle.logMessage(LOGMESSAGE, "Reset stomachcontent component", this->getName());
246 :			}
247 :
248 :			void SC::printSummary(ofstream& outfile, double weight) {
249 :			int year, area;
250 :
251 :			for (year = 0; year < likelihoodValues.Nrow(); year++) {
252 :			for (area = 0; area < likelihoodValues.Ncol(year); area++) {
253 :			outfile << setw(lowwidth) << Years[year] << sep << setw(lowwidth)
254 :			<< Steps[year] << sep << setw(printwidth) << areaindex[area] << sep
255 :			<< setw(largewidth) << this->getName() << sep << setprecision(smallprecision)
256 :			<< setw(smallwidth) << weight << sep << setprecision(largeprecision)
257 :			<< setw(largewidth) << likelihoodValues[year][area] << endl;
258 :			}
259 :			}
260 :			outfile.flush();
261 :			}
262 :
263 :			void SC::printLikelihood(ofstream& outfile, const TimeClass* const TimeInfo) {
264 :
265 :			if (!AAT.atCurrentTime(TimeInfo))
266 :			return;
267 :
268 :			int i, area, pred, prey;
269 :			timeindex = -1;
270 :			for (i = 0; i < Years.Size(); i++)
271 :			if ((Years[i] == TimeInfo->getYear()) && (Steps[i] == TimeInfo->getStep()))
272 :			timeindex = i;
273 :			if (timeindex == -1)
274 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - invalid timestep");
275 :
276 :			for (area = 0; area < modelConsumption.Ncol(timeindex); area++) {
277 :			for (pred = 0; pred < modelConsumption[timeindex][area]->Nrow(); pred++) {
278 :			for (prey = 0; prey < modelConsumption[timeindex][area]->Ncol(pred); prey++) {
279 :			outfile << setw(lowwidth) << Years[timeindex] << sep << setw(lowwidth)
280 :			<< Steps[timeindex] << sep << setw(printwidth) << areaindex[area] << sep
281 :			<< setw(printwidth) << predindex[pred] << sep << setw(printwidth)
282 :			<< preyindex[prey] << sep << setprecision(largeprecision) << setw(largewidth);
283 :
284 :			//JMB crude filter to remove the 'silly' values from the output
285 :			if ((*modelConsumption[timeindex][area])[pred][prey] < rathersmall)
286 :			outfile << 0 << endl;
287 :			else
288 :			outfile << (*modelConsumption[timeindex][area])[pred][prey] << endl;
289 :			}
290 :			}
291 :			}
292 :			}
293 :
294 :			SC::~SC() {
295 :			int i, j;
296 :			for (i = 0; i < obsConsumption.Nrow(); i++) {
297 :			for (j = 0; j < obsConsumption[i].Size(); j++) {
298 :			delete obsConsumption[i][j];
299 :			delete modelConsumption[i][j];
300 :			}
301 :			}
302 :
303 :			for (i = 0; i < preyindex.Size(); i++) {
304 :			delete aggregator[i];
305 :			delete preyLgrpDiv[i];
306 :			for (j = 0; j < preynames[i].Size(); j++)
307 :			delete[] preynames[i][j];
308 :			}
309 :
310 :			if (aggregator != 0) {
311 :			delete[] aggregator;
312 :			aggregator = 0;
313 :			}
314 :			if (preyLgrpDiv != 0) {
315 :			delete[] preyLgrpDiv;
316 :			preyLgrpDiv = 0;
317 :			}
318 :
319 :			if (!usepredages)
320 :			delete predLgrpDiv;
321 :
322 :			for (i = 0; i < predatornames.Size(); i++)
323 :			delete[] predatornames[i];
324 :			for (i = 0; i < areaindex.Size(); i++)
325 :			delete[] areaindex[i];
326 :			for (i = 0; i < predindex.Size(); i++)
327 :			delete[] predindex[i];
328 :			for (i = 0; i < preyindex.Size(); i++)
329 :			delete[] preyindex[i];
330 :			}
331 :
332 :			void SC::setPredatorsAndPreys(PredatorPtrVector& Predators, PreyPtrVector& Preys) {
333 :			int i, j, k, l, found;
334 :			int minage, maxage;
335 :			PredatorPtrVector predators;
336 :			aggregator = new PredatorAggregator*[preyindex.Size()];
337 :
338 :			for (i = 0; i < predatornames.Size(); i++) {
339 :			found = 0;
340 :			for (j = 0; j < Predators.Size(); j++)
341 :			if (strcasecmp(predatornames[i], Predators[j]->getName()) == 0) {
342 :			found++;
343 :			predators.resize(Predators[j]);
344 :			}
345 :
346 :			if (found == 0)
347 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - failed to match predator", predatornames[i]);
348 :			}
349 :
350 :			//check predator areas
351 :			if (handle.getLogLevel() >= LOGWARN) {
352 :			for (j = 0; j < areas.Nrow(); j++) {
353 :			found = 0;
354 :			for (i = 0; i < predators.Size(); i++)
355 :			for (k = 0; k < areas.Ncol(j); k++)
356 :			if (predators[i]->isInArea(areas[j][k]))
357 :			found++;
358 :			if (found == 0)
359 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - predator not defined on all areas");
360 :			}
361 :			}
362 :
363 :			//check that the predators are stocks and not fleets
364 :			for (i = 0; i < predators.Size(); i++)
365 :			if (predators[i]->getType() != STOCKPREDATOR)
366 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - cannot aggregate predator", predators[i]->getName());
367 :
368 :			preyLgrpDiv = new LengthGroupDivision*[preyindex.Size()];
369 :			if (!usepredages) {
370 :			predLgrpDiv = new LengthGroupDivision(predatorlengths);
371 :			if (predLgrpDiv->Error())
372 :			handle.logFileMessage(LOGFAIL, "\nError in stomachcontent - failed to create length group");
373 :			}
374 :
375 :			if (handle.getLogLevel() >= LOGWARN) {
376 :			if (usepredages) {
377 :			//check predator ages
378 :			minage = 9999;
379 :			maxage = -1;
380 :			for (j = 0; j < predatorages.Nrow(); j++) {
381 :			for (k = 0; k < predatorages.Ncol(j); k++) {
382 :			minage = min(predatorages[j][k], minage);
383 :			maxage = max(predatorages[j][k], maxage);
384 :			}
385 :			}
386 :
387 :			found = 0;
388 :			for (j = 0; j < predators.Size(); j++)
389 :			if (minage >= ((StockPredator*)predators[j])->minAge())
390 :			found++;
391 :			if (found == 0)
392 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - minimum age less than predator age");
393 :
394 :			found = 0;
395 :			for (j = 0; j < predators.Size(); j++)
396 :			if (maxage <= ((StockPredator*)predators[j])->maxAge())
397 :			found++;
398 :			if (found == 0)
399 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - maximum age greater than predator age");
400 :
401 :			} else {
402 :			//check predator lengths
403 :			found = 0;
404 :			for (j = 0; j < predators.Size(); j++)
405 :			if (predLgrpDiv->maxLength(0) > predators[j]->getLengthGroupDiv()->minLength())
406 :			found++;
407 :			if (found == 0)
408 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - minimum length group less than predator length");
409 :
410 :			found = 0;
411 :			for (j = 0; j < predators.Size(); j++)
412 :			if (predLgrpDiv->minLength(predLgrpDiv->numLengthGroups()) < predators[j]->getLengthGroupDiv()->maxLength())
413 :			found++;
414 :			if (found == 0)
415 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - maximum length group greater than predator length");
416 :			}
417 :			}
418 :
419 :			for (i = 0; i < preynames.Nrow(); i++) {
420 :			PreyPtrVector preys;
421 :			for (j = 0; j < preynames.Ncol(i); j++) {
422 :			found = 0;
423 :			for (k = 0; k < Preys.Size(); k++)
424 :			if (strcasecmp(preynames[i][j], Preys[k]->getName()) == 0) {
425 :			found++;
426 :			preys.resize(Preys[k]);
427 :			}
428 :
429 :			if (found == 0)
430 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - failed to match prey", preynames[i][j]);
431 :			}
432 :
433 :			//check prey areas
434 :			if (handle.getLogLevel() >= LOGWARN) {
435 :			for (j = 0; j < areas.Nrow(); j++) {
436 :			found = 0;
437 :			for (k = 0; k < preys.Size(); k++)
438 :			for (l = 0; l < areas.Ncol(j); l++)
439 :			if (preys[k]->isInArea(areas[j][l]))
440 :			found++;
441 :			if (found == 0)
442 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - prey not defined on all areas");
443 :			}
444 :			}
445 :
446 :			//resize the digestion matrix
447 :			digestion.AddRows(1, preylengths[i].Size(), 0.0);
448 :			preyLgrpDiv[i] = new LengthGroupDivision(preylengths[i]);
449 :			if (preyLgrpDiv[i]->Error())
450 :			handle.logFileMessage(LOGFAIL, "\nError in stomachcontent - failed to create length group");
451 :
452 :			if (usepredages)
453 :			aggregator[i] = new PredatorAggregator(predators, preys, areas, predatorages, preyLgrpDiv[i]);
454 :			else
455 :			aggregator[i] = new PredatorAggregator(predators, preys, areas, predLgrpDiv, preyLgrpDiv[i]);
456 :			}
457 :			}
458 :
459 :			void SC::Print(ofstream& outfile) const {
460 :			int i, j;
461 :
462 :			outfile << "\tPredators:\n\t\t";
463 :			for (i = 0; i < predatornames.Size(); i++)
464 :			outfile << predatornames[i] << sep;
465 :
466 :			if (usepredages) {
467 :			outfile << "\n\t\tages:";
468 :			for (i = 0; i < predatorages.Nrow(); i++) {
469 :			outfile << "\n\t\t\t";
470 :			for (j = 0; j < predatorages.Ncol(i); j++)
471 :			outfile << predatorages[i][j] << sep;
472 :			}
473 :			outfile << endl;
474 :			} else {
475 :			outfile << "\n\t\tlengths: ";
476 :			for (i = 0; i < predatorlengths.Size(); i++)
477 :			outfile << predatorlengths[i] << sep;
478 :			outfile << endl;
479 :			}
480 :
481 :			outfile << "\n\tPreys:";
482 :			for (i = 0; i < preyindex.Size(); i++) {
483 :			outfile << "\n\t\t" << preyindex[i] << "\n\t\t";
484 :			for (j = 0; j < preynames[i].Size(); j++)
485 :			outfile << preynames[i][j] << sep;
486 :			outfile << "\n\t\tlengths: ";
487 :			for (j = 0; j < preylengths[i].Size(); j++)
488 :			outfile << preylengths[i][j] << sep;
489 :			outfile << endl;
490 :			aggregator[i]->Print(outfile);
491 :			}
492 :			outfile.flush();
493 :			}
494 :
495 :			// ********************************************************
496 :			// Functions for SCNumbers
497 :			// ********************************************************
498 :			SCNumbers::SCNumbers(CommentStream& infile, const AreaClass* const Area,
499 :			const TimeClass* const TimeInfo, Keeper* const keeper,
500 :			const char* datafilename, const char* givenname)
501 :			: SC(infile, Area, TimeInfo, keeper, datafilename, givenname) {
502 :
503 :			ifstream datafile;
504 :			CommentStream subdata(datafile);
505 :			//read in stomach content from file
506 :			datafile.open(datafilename, ios::in);
507 :			handle.checkIfFailure(datafile, datafilename);
508 :			handle.Open(datafilename);
509 :			readStomachNumberContent(subdata, TimeInfo);
510 :			handle.Close();
511 :			datafile.close();
512 :			datafile.clear();
513 :
514 :			MN = Multinomial();
515 :			MN.setValue(epsilon);
516 :			mndist.resize(likelihoodValues.Nrow(), 0.0);
517 :			mndata.resize(likelihoodValues.Nrow(), 0.0);
518 :			}
519 :
520 :			void SCNumbers::readStomachNumberContent(CommentStream& infile, const TimeClass* const TimeInfo) {
521 :
522 :			int i, year, step, count, reject;
523 :			double tmpnumber;
524 :			char tmparea[MaxStrLength], tmppred[MaxStrLength], tmpprey[MaxStrLength];
525 :			strncpy(tmparea, "", MaxStrLength);
526 :			strncpy(tmppred, "", MaxStrLength);
527 :			strncpy(tmpprey, "", MaxStrLength);
528 :			int keepdata, timeid, areaid, predid, preyid;
529 :
530 :			if (usepredages) //age structured predator
531 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - age based predators cannot be used with scnumbers");
532 :
533 :			int numpred = predatorlengths.Size() - 1;
534 :			int numarea = areas.Nrow();
535 :			int numprey = 0;
536 :			for (i = 0; i < preylengths.Nrow(); i++)
537 :			numprey += preylengths[i].Size() - 1;
538 :			if (numprey == 0)
539 :			handle.logMessage(LOGWARN, "Warning in stomachcontents - no prey found for", this->getName());
540 :
541 :			//Check the number of columns in the inputfile
542 :			infile >> ws;
543 :			if (countColumns(infile) != 6)
544 :			handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 6");
545 :
546 :			year = step = count = reject = 0;
547 :			while (!infile.eof()) {
548 :			keepdata = 1;
549 :			infile >> year >> step >> tmparea >> tmppred >> tmpprey >> tmpnumber >> ws;
550 :
551 :			//crude check to see if something has gone wrong and avoid infinite loops
552 :			if (strlen(tmparea) == 0)
553 :			handle.logFileMessage(LOGFAIL, "failed to read data from file");
554 :
555 :			//if tmparea is in areaindex find areaid, else dont keep the data
556 :			areaid = -1;
557 :			for (i = 0; i < areaindex.Size(); i++)
558 :			if (strcasecmp(areaindex[i], tmparea) == 0)
559 :			areaid = i;
560 :
561 :			if (areaid == -1)
562 :			keepdata = 0;
563 :
564 :			//if tmppred is in predindex find predid, else dont keep the data
565 :			predid = -1;
566 :			for (i = 0; i < predindex.Size(); i++)
567 :			if (strcasecmp(predindex[i], tmppred) == 0)
568 :			predid = i;
569 :
570 :			if (predid == -1)
571 :			keepdata = 0;
572 :
573 :			//if tmpprey is in preyindex find preyid, else dont keep the data
574 :			preyid = -1;
575 :			for (i = 0; i < preyindex.Size(); i++)
576 :			if (strcasecmp(preyindex[i], tmpprey) == 0)
577 :			preyid = i;
578 :
579 :			if (preyid == -1)
580 :			keepdata = 0;
581 :
582 :			//check if the year and step are in the simulation
583 :			timeid = -1;
584 :			if ((TimeInfo->isWithinPeriod(year, step)) && (keepdata == 1)) {
585 :			//if this is a new timestep, resize to store the data
586 :			for (i = 0; i < Years.Size(); i++)
587 :			if ((Years[i] == year) && (Steps[i] == step))
588 :			timeid = i;
589 :
590 :			if (timeid == -1) {
591 :			Years.resize(1, year);
592 :			Steps.resize(1, step);
593 :			timeid = Years.Size() - 1;
594 :
595 :			obsConsumption.resize();
596 :			modelConsumption.resize();
597 :			likelihoodValues.AddRows(1, numarea, 0.0);
598 :			for (i = 0; i < numarea; i++) {
599 :			obsConsumption[timeid].resize(new DoubleMatrix(numpred, numprey, 0.0));
600 :			modelConsumption[timeid].resize(new DoubleMatrix(numpred, numprey, 0.0));
601 :			}
602 :			}
603 :
604 :			} else
605 :			keepdata = 0;
606 :
607 :			if (keepdata == 1) {
608 :			//stomach content data is required, so store it
609 :			count++;
610 :			(*obsConsumption[timeid][areaid])[predid][preyid] = tmpnumber;
611 :			} else
612 :			reject++; //count number of rejected data points read from file
613 :			}
614 :
615 :			AAT.addActions(Years, Steps, TimeInfo);
616 :			if (count == 0)
617 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - found no data in the data file for", this->getName());
618 :			if (reject != 0)
619 :			handle.logMessage(LOGMESSAGE, "Discarded invalid stomachcontent data - number of invalid entries", reject);
620 :			handle.logMessage(LOGMESSAGE, "Read stomachcontent data file - number of entries", count);
621 :			}
622 :
623 :			double SCNumbers::calcLikelihood() {
624 :			int a, pred, prey;
625 :			MN.Reset();
626 :			for (a = 0; a < areas.Nrow(); a++) {
627 :			likelihoodValues[timeindex][a] = 0.0;
628 :			for (prey = 0; prey < obsConsumption[timeindex][a]->Ncol(0); prey++) {
629 :			for (pred = 0; pred < mndata.Size(); pred++) {
630 :			mndata[pred] = (*obsConsumption[timeindex][a])[pred][prey];
631 :			mndist[pred] = (*modelConsumption[timeindex][a])[pred][prey];
632 :			}
633 :			likelihoodValues[timeindex][a] += MN.calcLogLikelihood(mndata, mndist);
634 :			}
635 :			}
636 :			return MN.getLogLikelihood();
637 :			}
638 :
639 :			void SCNumbers::aggregate(int i) {
640 :			aggregator[i]->NumberSum();
641 :			}
642 :
643 :			// ********************************************************
644 :			// Functions for SCAmounts
645 :			// ********************************************************
646 :			SCAmounts::SCAmounts(CommentStream& infile, const AreaClass* const Area,
647 :			const TimeClass* const TimeInfo, Keeper* const keeper,
648 :			const char* datafilename, const char* numfilename, const char* givenname)
649 :			: SC(infile, Area, TimeInfo, keeper, datafilename, givenname) {
650 :
651 :			ifstream datafile;
652 :			CommentStream subdata(datafile);
653 :			//read in stomach content amounts from file
654 :			datafile.open(datafilename, ios::in);
655 :			handle.checkIfFailure(datafile, datafilename);
656 :			handle.Open(datafilename);
657 :			readStomachAmountContent(subdata, TimeInfo);
658 :			handle.Close();
659 :			datafile.close();
660 :			datafile.clear();
661 :
662 :			//read in stomach content sample size from file
663 :			datafile.open(numfilename, ios::in);
664 :			handle.checkIfFailure(datafile, numfilename);
665 :			handle.Open(numfilename);
666 :			readStomachSampleContent(subdata, TimeInfo);
667 :			handle.Close();
668 :			datafile.close();
669 :			datafile.clear();
670 :			}
671 :
672 :			void SCAmounts::readStomachAmountContent(CommentStream& infile, const TimeClass* const TimeInfo) {
673 :
674 :			int i, year, step, count, reject;
675 :			double tmpnumber, tmpstddev;
676 :			char tmparea[MaxStrLength], tmppred[MaxStrLength], tmpprey[MaxStrLength];
677 :			strncpy(tmparea, "", MaxStrLength);
678 :			strncpy(tmppred, "", MaxStrLength);
679 :			strncpy(tmpprey, "", MaxStrLength);
680 :			int keepdata, timeid, areaid, predid, preyid;
681 :
682 :			int numpred = 0;
683 :			if (usepredages) //age structured predator
684 :			numpred = predatorages.Nrow();
685 :			else
686 :			numpred = predatorlengths.Size() - 1;
687 :
688 :			int numarea = areas.Nrow();
689 :			int numprey = 0;
690 :			for (i = 0; i < preylengths.Nrow(); i++)
691 :			numprey += preylengths[i].Size() - 1;
692 :			if (numprey == 0)
693 :			handle.logMessage(LOGWARN, "Warning in stomachcontents - no prey found for", this->getName());
694 :
695 :			//Check the number of columns in the inputfile
696 :			infile >> ws;
697 :			if (countColumns(infile) != 7)
698 :			handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 7");
699 :
700 :			year = step = count = reject = 0;
701 :			while (!infile.eof()) {
702 :			keepdata = 1;
703 :			infile >> year >> step >> tmparea >> tmppred >> tmpprey >> tmpnumber >> tmpstddev >> ws;
704 :
705 :			//crude check to see if something has gone wrong and avoid infinite loops
706 :			if (strlen(tmparea) == 0)
707 :			handle.logFileMessage(LOGFAIL, "failed to read data from file");
708 :
709 :			//if tmparea is in areaindex find areaid, else dont keep the data
710 :			areaid = -1;
711 :			for (i = 0; i < areaindex.Size(); i++)
712 :			if (strcasecmp(areaindex[i], tmparea) == 0)
713 :			areaid = i;
714 :
715 :			if (areaid == -1)
716 :			keepdata = 0;
717 :
718 :			//if tmppred is in predindex find predid, else dont keep the data
719 :			predid = -1;
720 :			for (i = 0; i < predindex.Size(); i++)
721 :			if (strcasecmp(predindex[i], tmppred) == 0)
722 :			predid = i;
723 :
724 :			if (predid == -1)
725 :			keepdata = 0;
726 :
727 :			//if tmpprey is in preyindex find preyid, else dont keep the data
728 :			preyid = -1;
729 :			for (i = 0; i < preyindex.Size(); i++)
730 :			if (strcasecmp(preyindex[i], tmpprey) == 0)
731 :			preyid = i;
732 :
733 :			if (preyid == -1)
734 :			keepdata = 0;
735 :
736 :			//check if the year and step are in the simulation
737 :			timeid = -1;
738 :			if ((TimeInfo->isWithinPeriod(year, step)) && (keepdata == 1)) {
739 :			//if this is a new timestep, resize to store the data
740 :			for (i = 0; i < Years.Size(); i++)
741 :			if ((Years[i] == year) && (Steps[i] == step))
742 :			timeid = i;
743 :
744 :			if (timeid == -1) {
745 :			Years.resize(1, year);
746 :			Steps.resize(1, step);
747 :			timeid = Years.Size() - 1;
748 :
749 :			obsConsumption.resize();
750 :			modelConsumption.resize();
751 :			stddev.resize();
752 :			likelihoodValues.AddRows(1, numarea, 0.0);
753 :			for (i = 0; i < numarea; i++) {
754 :			obsConsumption[timeid].resize(new DoubleMatrix(numpred, numprey, 0.0));
755 :			modelConsumption[timeid].resize(new DoubleMatrix(numpred, numprey, 0.0));
756 :			stddev[timeid].resize(new DoubleMatrix(numpred, numprey, 0.0));
757 :			}
758 :			}
759 :
760 :			} else
761 :			keepdata = 0;
762 :
763 :			if (keepdata == 1) {
764 :			//stomach content data is required, so store it
765 :			count++;
766 :			(*obsConsumption[timeid][areaid])[predid][preyid] = tmpnumber;
767 :			(*stddev[timeid][areaid])[predid][preyid] = tmpstddev;
768 :			} else
769 :			reject++; //count number of rejected data points read from file
770 :			}
771 :
772 :			AAT.addActions(Years, Steps, TimeInfo);
773 :			if (count == 0)
774 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - found no data in the data file for", this->getName());
775 :			if (reject != 0)
776 :			handle.logMessage(LOGMESSAGE, "Discarded invalid stomachcontent data - number of invalid entries", reject);
777 :			handle.logMessage(LOGMESSAGE, "Read stomachcontent data file - number of entries", count);
778 :			}
779 :
780 :			void SCAmounts::readStomachSampleContent(CommentStream& infile, const TimeClass* const TimeInfo) {
781 :
782 :			int i, year, step, count, reject;
783 :			double tmpnumber;
784 :			int keepdata, timeid, areaid, predid;
785 :			char tmparea[MaxStrLength], tmppred[MaxStrLength];
786 :			strncpy(tmparea, "", MaxStrLength);
787 :			strncpy(tmppred, "", MaxStrLength);
788 :
789 :			int numpred = 0;
790 :			if (usepredages) //age structured predator
791 :			numpred = predatorages.Nrow();
792 :			else
793 :			numpred = predatorlengths.Size() - 1;
794 :
795 :			//We know the size that numbers[] will be from obsConsumption
796 :			int numarea = areas.Nrow();
797 :			for (i = 0; i < obsConsumption.Nrow(); i++)
798 :			number.resize(new DoubleMatrix(numarea, numpred, 0.0));
799 :
800 :			//Check the number of columns in the inputfile
801 :			infile >> ws;
802 :			if (countColumns(infile) != 5)
803 :			handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 5");
804 :
805 :			year = step = count = reject = 0;
806 :			while (!infile.eof()) {
807 :			keepdata = 1;
808 :			infile >> year >> step >> tmparea >> tmppred >> tmpnumber >> ws;
809 :
810 :			//crude check to see if something has gone wrong and avoid infinite loops
811 :			if (strlen(tmparea) == 0)
812 :			handle.logFileMessage(LOGFAIL, "failed to read data from file");
813 :
814 :			//check if the year and step are in the simulation
815 :			timeid = -1;
816 :			if (TimeInfo->isWithinPeriod(year, step))
817 :			//find the timeid from Years and Steps
818 :			for (i = 0; i < Years.Size(); i++)
819 :			if ((Years[i] == year) && (Steps[i] == step))
820 :			timeid = i;
821 :
822 :			if (timeid == -1)
823 :			keepdata = 0;
824 :
825 :			//if tmparea is in areaindex find areaid, else dont keep the data
826 :			areaid = -1;
827 :			for (i = 0; i < areaindex.Size(); i++)
828 :			if (strcasecmp(areaindex[i], tmparea) == 0)
829 :			areaid = i;
830 :
831 :			if (areaid == -1)
832 :			keepdata = 0;
833 :
834 :			//if tmppred is in predindex find predid, else dont keep the data
835 :			predid = -1;
836 :			for (i = 0; i < predindex.Size(); i++)
837 :			if (strcasecmp(predindex[i], tmppred) == 0)
838 :			predid = i;
839 :
840 :			if (predid == -1)
841 :			keepdata = 0;
842 :
843 :			if (keepdata == 1) {
844 :			//stomach content data is required, so store it
845 :			count++;
846 :			(*number[timeid])[areaid][predid] = tmpnumber;
847 :			} else
848 :			reject++; //count number of rejected data points read from file
849 :			}
850 :
851 :			if (count == 0)
852 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - found no data in the data file for", this->getName());
853 :			if (reject != 0)
854 :			handle.logMessage(LOGMESSAGE, "Discarded invalid stomachcontent data - number of invalid entries", reject);
855 :			handle.logMessage(LOGMESSAGE, "Read stomachcontent data file - number of entries", count);
856 :			}
857 :
858 :			SCAmounts::~SCAmounts() {
859 :			int i, j;
860 :			for (i = 0; i < stddev.Nrow(); i++) {
861 :			delete number[i];
862 :			for (j = 0; j < stddev[i].Size(); j++)
863 :			delete stddev[i][j];
864 :			}
865 :			}
866 :
867 :			//JMB - note this ignores the number of samples ...
868 :			void SCAmounts::printLikelihood(ofstream& outfile, const TimeClass* const TimeInfo) {
869 :
870 :			if (!AAT.atCurrentTime(TimeInfo))
871 :			return;
872 :
873 :			int i, area, pred, prey;
874 :			timeindex = -1;
875 :			for (i = 0; i < Years.Size(); i++)
876 :			if ((Years[i] == TimeInfo->getYear()) && (Steps[i] == TimeInfo->getStep()))
877 :			timeindex = i;
878 :			if (timeindex == -1)
879 :			handle.logMessage(LOGFAIL, "Error in stomachcontent - invalid timestep");
880 :
881 :			for (area = 0; area < modelConsumption.Ncol(timeindex); area++) {
882 :			for (pred = 0; pred < modelConsumption[timeindex][area]->Nrow(); pred++) {
883 :			for (prey = 0; prey < modelConsumption[timeindex][area]->Ncol(pred); prey++) {
884 :			outfile << setw(lowwidth) << Years[timeindex] << sep << setw(lowwidth)
885 :			<< Steps[timeindex] << sep << setw(printwidth) << areaindex[area] << sep
886 :			<< setw(printwidth) << predindex[pred] << sep << setw(printwidth)
887 :			<< preyindex[prey] << sep << setprecision(largeprecision) << setw(largewidth);
888 :
889 :			//JMB crude filter to remove the 'silly' values from the output
890 :			if ((*modelConsumption[timeindex][area])[pred][prey] < rathersmall)
891 :			outfile << 0;
892 :			else
893 :			outfile << (*modelConsumption[timeindex][area])[pred][prey];
894 :
895 :			outfile << sep << setprecision(largeprecision) << setw(largewidth)
896 :			<< (*stddev[timeindex][area])[pred][prey] << endl;
897 :			}
898 :			}
899 :			}
900 :			}
901 :
902 :			double SCAmounts::calcLikelihood() {
903 :			int a, pred, prey;
904 :			double tmplik, lik = 0.0;
905 :
906 :			for (a = 0; a < areas.Nrow(); a++) {
907 :			likelihoodValues[timeindex][a] = 0.0;
908 :			for (pred = 0; pred < obsConsumption[timeindex][a]->Nrow(); pred++) {
909 :			if (!(isZero((*number[timeindex])[a][pred]))) {
910 :			tmplik = 0.0;
911 :			for (prey = 0; prey < obsConsumption[timeindex][a]->Ncol(pred); prey++) {
912 :			if (!(isZero((*stddev[timeindex][a])[pred][prey])))
913 :			tmplik += ((*modelConsumption[timeindex][a])[pred][prey] -
914 :			(*obsConsumption[timeindex][a])[pred][prey]) *
915 :			((*modelConsumption[timeindex][a])[pred][prey] -
916 :			(*obsConsumption[timeindex][a])[pred][prey]) /
917 :			((*stddev[timeindex][a])[pred][prey] * (*stddev[timeindex][a])[pred][prey]);
918 :			}
919 :			tmplik *= (*number[timeindex])[a][pred];
920 :			likelihoodValues[timeindex][a] += tmplik;
921 :			}
922 :			}
923 :			lik += likelihoodValues[timeindex][a];
924 :			}
925 :			return lik;
926 :			}
927 :
928 :			// ********************************************************
929 :			// Functions for SCRatios
930 :			// ********************************************************
931 :			void SCRatios::setPredatorsAndPreys(PredatorPtrVector& Predators, PreyPtrVector& Preys) {
932 :			int i, j, k, l;
933 :			double tmpdivide, scale;
934 :			SC::setPredatorsAndPreys(Predators, Preys);
935 :			//Scale each row such that it sums up to 1
936 :			for (i = 0; i < obsConsumption.Nrow(); i++) {
937 :			for (j = 0; j < obsConsumption.Ncol(i); j++) {
938 :			for (k = 0; k < obsConsumption[i][j]->Nrow(); k++) {
939 :			scale = 0.0;
940 :			for (l = 0; l < obsConsumption[i][j]->Ncol(k); l++)
941 :			scale += (*obsConsumption[i][j])[k][l];
942 :
943 :			if (!(isZero(scale))) {
944 :			tmpdivide = 1.0 / scale;
945 :			for (l = 0; l < obsConsumption[i][j]->Ncol(k); l++)
946 :			(*obsConsumption[i][j])[k][l] *= tmpdivide;
947 :			}
948 :			}
949 :			}
950 :			}
951 :			}
952 :
953 :			double SCRatios::calcLikelihood() {
954 :			int a, pred, prey;
955 :			double scale, tmplik, tmpdivide;
956 :			double lik = 0.0;
957 :
958 :			for (a = 0; a < areas.Nrow(); a++) {
959 :			likelihoodValues[timeindex][a] = 0.0;
960 :			for (pred = 0; pred < obsConsumption[timeindex][a]->Nrow(); pred++) {
961 :			scale = 0.0;
962 :			for (prey = 0; prey < modelConsumption[timeindex][a]->Ncol(pred); prey++)
963 :			scale += (*modelConsumption[timeindex][a])[pred][prey];
964 :
965 :			if (!(isZero(scale))) {
966 :			tmpdivide = 1.0 / scale;
967 :			for (prey = 0; prey < obsConsumption[timeindex][a]->Ncol(pred); prey++)
968 :			(*modelConsumption[timeindex][a])[pred][prey] *= tmpdivide;
969 :
970 :			if (!(isZero((*number[timeindex])[a][pred]))) {
971 :			tmplik = 0.0;
972 :			for (prey = 0; prey < obsConsumption[timeindex][a]->Ncol(pred); prey++) {
973 :			if (!(isZero((*stddev[timeindex][a])[pred][prey])))
974 :			tmplik += ((*modelConsumption[timeindex][a])[pred][prey] -
975 :			(*obsConsumption[timeindex][a])[pred][prey]) *
976 :			((*modelConsumption[timeindex][a])[pred][prey] -
977 :			(*obsConsumption[timeindex][a])[pred][prey]) /
978 :			((*stddev[timeindex][a])[pred][prey] * (*stddev[timeindex][a])[pred][prey]);
979 :			}
980 :			tmplik *= (*number[timeindex])[a][pred];
981 :			likelihoodValues[timeindex][a] += tmplik;
982 :			}
983 :			}
984 :			}
985 :			lik += likelihoodValues[timeindex][a];
986 :			}
987 :			return lik;
988 :			}
989 :
990 :			// ********************************************************
991 :			// Functions for SCSimple
992 :			// ********************************************************
993 :			SCSimple::SCSimple(CommentStream& infile, const AreaClass* const Area,
994 :			const TimeClass* const TimeInfo, Keeper* const keeper,
995 :			const char* datafilename, const char* givenname)
996 :			: SC(infile, Area, TimeInfo, keeper, datafilename, givenname) {
997 :
998 :			ifstream datafile;
999 :			CommentStream subdata(datafile);
1000 :			//read in stomach content from file
1001 :			datafile.open(datafilename, ios::in);
1002 :			handle.checkIfFailure(datafile, datafilename);
1003 :			handle.Open(datafilename);
1004 :			readStomachSimpleContent(subdata, TimeInfo);
1005 :			handle.Close();
1006 :			datafile.close();
1007 :			datafile.clear();
1008 :			}
1009 :
1010 :			void SCSimple::readStomachSimpleContent(CommentStream& infile, const TimeClass* const TimeInfo) {
1011 :
1012 :			int i, year, step, count, reject;
1013 :			double tmpnumber;
1014 :			char tmparea[MaxStrLength], tmppred[MaxStrLength], tmpprey[MaxStrLength];
1015 :			strncpy(tmparea, "", MaxStrLength);
1016 :			strncpy(tmppred, "", MaxStrLength);
1017 :			strncpy(tmpprey, "", MaxStrLength);
1018 :			int keepdata, timeid, areaid, predid, preyid;
1019 :
1020 :			int numpred = 0;
1021 :			if (usepredages) //age structured predator
1022 :			numpred = predatorages.Nrow();
1023 :			else
1024 :			numpred = predatorlengths.Size() - 1;
1025 :
1026 :			int numarea = areas.Nrow();
1027 :			int numprey = 0;
1028 :			for (i = 0; i < preylengths.Nrow(); i++)
1029 :			numprey += preylengths[i].Size() - 1;
1030 :			if (numprey == 0)
1031 :			handle.logMessage(LOGWARN, "Warning in stomachcontents - no prey found for", this->getName());
1032 :
1033 :			//Check the number of columns in the inputfile
1034 :			infile >> ws;
1035 :			if (countColumns(infile) != 6)
1036 :			handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 6");
1037 :
1038 :			year = step = count = reject = 0;
1039 :			while (!infile.eof()) {
1040 :			keepdata = 1;
1041 :			infile >> year >> step >> tmparea >> tmppred >> tmpprey >> tmpnumber >> ws;
1042 :
1043 :			//crude check to see if something has gone wrong and avoid infinite loops
1044 :			if (strlen(tmparea) == 0)
1045 :			handle.logFileMessage(LOGFAIL, "failed to read data from file");
1046 :
1047 :			//if tmparea is in areaindex find areaid, else dont keep the data
1048 :			areaid = -1;
1049 :			for (i = 0; i < areaindex.Size(); i++)
1050 :			if (strcasecmp(areaindex[i], tmparea) == 0)
1051 :			areaid = i;
1052 :
1053 :			if (areaid == -1)
1054 :			keepdata = 0;
1055 :
1056 :			//if tmppred is in predindex find predid, else dont keep the data
1057 :			predid = -1;
1058 :			for (i = 0; i < predindex.Size(); i++)
1059 :			if (strcasecmp(predindex[i], tmppred) == 0)
1060 :			predid = i;
1061 :
1062 :			if (predid == -1)
1063 :			keepdata = 0;
1064 :
1065 :			//if tmpprey is in preyindex find preyid, else dont keep the data
1066 :			preyid = -1;
1067 :			for (i = 0; i < preyindex.Size(); i++)
1068 :			if (strcasecmp(preyindex[i], tmpprey) == 0)
1069 :			preyid = i;
1070 :
1071 :			if (preyid == -1)
1072 :			keepdata = 0;
1073 :
1074 :			//check if the year and step are in the simulation
1075 :			timeid = -1;
1076 :			if ((TimeInfo->isWithinPeriod(year, step)) && (keepdata == 1)) {
1077 :			//if this is a new timestep, resize to store the data
1078 :			for (i = 0; i < Years.Size(); i++)
1079 :			if ((Years[i] == year) && (Steps[i] == step))
1080 :			timeid = i;
1081 :
1082 :			if (timeid == -1) {
1083 :			Years.resize(1, year);
1084 :			Steps.resize(1, step);
1085 :			timeid = Years.Size() - 1;
1086 :
1087 :			obsConsumption.resize();
1088 :			modelConsumption.resize();
1089 :			likelihoodValues.AddRows(1, numarea, 0.0);
1090 :			for (i = 0; i < numarea; i++) {
1091 :			obsConsumption[timeid].resize(new DoubleMatrix(numpred, numprey, 0.0));
1092 :			modelConsumption[timeid].resize(new DoubleMatrix(numpred, numprey, 0.0));
1093 :			}
1094 :			}
1095 :
1096 :			} else
1097 :			keepdata = 0;
1098 :
1099 :			if (keepdata == 1) {
1100 :			//stomach content data is required, so store it
1101 :			count++;
1102 :			(*obsConsumption[timeid][areaid])[predid][preyid] = tmpnumber;
1103 :			} else
1104 :			reject++; //count number of rejected data points read from file
1105 :			}
1106 :
1107 :			AAT.addActions(Years, Steps, TimeInfo);
1108 :			if (count == 0)
1109 :			handle.logMessage(LOGWARN, "Warning in stomachcontent - found no data in the data file for", this->getName());
1110 :			if (reject != 0)
1111 :			handle.logMessage(LOGMESSAGE, "Discarded invalid stomachcontent data - number of invalid entries", reject);
1112 :			handle.logMessage(LOGMESSAGE, "Read stomachcontent data file - number of entries", count);
1113 :			}
1114 :
1115 :			void SCSimple::setPredatorsAndPreys(PredatorPtrVector& Predators, PreyPtrVector& Preys) {
1116 :			int i, j, k, l;
1117 :			double tmpdivide, scale;
1118 :			SC::setPredatorsAndPreys(Predators, Preys);
1119 :			//Scale each row such that it sums up to 1
1120 :			for (i = 0; i < obsConsumption.Nrow(); i++) {
1121 :			for (j = 0; j < obsConsumption.Ncol(i); j++) {
1122 :			for (k = 0; k < obsConsumption[i][j]->Nrow(); k++) {
1123 :			scale = 0.0;
1124 :			for (l = 0; l < obsConsumption[i][j]->Ncol(k); l++)
1125 :			scale += (*obsConsumption[i][j])[k][l];
1126 :
1127 :			if (!(isZero(scale))) {
1128 :			tmpdivide = 1.0 / scale;
1129 :			for (l = 0; l < obsConsumption[i][j]->Ncol(k); l++)
1130 :			(*obsConsumption[i][j])[k][l] *= tmpdivide;
1131 :			}
1132 :			}
1133 :			}
1134 :			}
1135 :			}
1136 :
1137 :			double SCSimple::calcLikelihood() {
1138 :			int a, pred, prey;
1139 :			double scale, tmplik, tmpdivide;
1140 :			double lik = 0.0;
1141 :
1142 :			for (a = 0; a < areas.Nrow(); a++) {
1143 :			likelihoodValues[timeindex][a] = 0.0;
1144 :			for (pred = 0; pred < obsConsumption[timeindex][a]->Nrow(); pred++) {
1145 :			scale = 0.0;
1146 :			for (prey = 0; prey < modelConsumption[timeindex][a]->Ncol(pred); prey++)
1147 :			scale += (*modelConsumption[timeindex][a])[pred][prey];
1148 :
1149 :			if (!(isZero(scale))) {
1150 :			tmpdivide = 1.0 / scale;
1151 :			tmplik = 0.0;
1152 :			for (prey = 0; prey < obsConsumption[timeindex][a]->Ncol(pred); prey++) {
1153 :			(*modelConsumption[timeindex][a])[pred][prey] *= tmpdivide;
1154 :			tmplik += ((*modelConsumption[timeindex][a])[pred][prey] -
1155 :			(*obsConsumption[timeindex][a])[pred][prey]) *
1156 :			((*modelConsumption[timeindex][a])[pred][prey] -
1157 :			(*obsConsumption[timeindex][a])[pred][prey]);
1158 :			}
1159 :			likelihoodValues[timeindex][a] += tmplik;
1160 :			}
1161 :			}
1162 :			lik += likelihoodValues[timeindex][a];
1163 :			}
1164 :			return lik;
1165 :			}

---

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