Annotation of /trunk/gadget/predatoraggregator.cc

Parent Directory | Revision Log

---

Revision 1 - (view) (download)

1 :	agomez	1	#include "predatoraggregator.h"
2 :			#include "prey.h"
3 :			#include "predator.h"
4 :			#include "stockpredator.h"
5 :			#include "mathfunc.h"
6 :			#include "popinfovector.h"
7 :			#include "errorhandler.h"
8 :			#include "gadget.h"
9 :			#include "global.h"
10 :
11 :			PredatorAggregator::PredatorAggregator(const PredatorPtrVector& Predators,
12 :			const PreyPtrVector& Preys, const IntMatrix& Areas,
13 :			const LengthGroupDivision* const predLgrpDiv, const LengthGroupDivision* const preyLgrpDiv)
14 :			: predators(Predators), preys(Preys), areas(Areas),
15 :			doeseat(Predators.Size(), Preys.Size(), 0), dptr(0), alk(0), usepredages(0) {
16 :
17 :			int i, j;
18 :			for (i = 0; i < predators.Size(); i++) {
19 :			if (!checkLengthGroupStructure(predators[i]->getLengthGroupDiv(), predLgrpDiv))
20 :			handle.logMessage(LOGFAIL, "Error in predatoraggregator - invalid predator length group structure");
21 :			predConv.AddRows(1, predators[i]->getLengthGroupDiv()->numLengthGroups(), -1);
22 :			for (j = 0; j < predConv.Ncol(i); j++)
23 :			predConv[i][j] = predLgrpDiv->numLengthGroup(predators[i]->getLengthGroupDiv()->meanLength(j));
24 :			}
25 :
26 :			for (i = 0; i < preys.Size(); i++) {
27 :			if (!checkLengthGroupStructure(preys[i]->getLengthGroupDiv(), preyLgrpDiv))
28 :			handle.logMessage(LOGFAIL, "Error in predatoraggregator - invalid prey length group structure");
29 :			preyConv.AddRows(1, preys[i]->getLengthGroupDiv()->numLengthGroups(), -1);
30 :			for (j = 0; j < preyConv.Ncol(i); j++)
31 :			preyConv[i][j] = preyLgrpDiv->numLengthGroup(preys[i]->getLengthGroupDiv()->meanLength(j));
32 :			}
33 :
34 :			for (i = 0; i < predators.Size(); i++)
35 :			for (j = 0; j < preys.Size(); j++)
36 :			if (predators[i]->doesEat(preys[j]->getName()))
37 :			doeseat[i][j] = 1;
38 :
39 :			for (i = 0; i < areas.Nrow(); i++)
40 :			total.resize(new DoubleMatrix(predLgrpDiv->numLengthGroups(), preyLgrpDiv->numLengthGroups(), 0.0));
41 :			}
42 :
43 :			PredatorAggregator::PredatorAggregator(const PredatorPtrVector& Predators,
44 :			const PreyPtrVector& Preys, const IntMatrix& Areas,
45 :			const IntMatrix& predAges, const LengthGroupDivision* const preyLgrpDiv)
46 :			: predators(Predators), preys(Preys), areas(Areas),
47 :			doeseat(Predators.Size(), Preys.Size(), 0), dptr(0), alk(0), usepredages(1) {
48 :
49 :			int i, j, k, l, minage, maxage;
50 :			for (i = 0; i < predators.Size(); i++) {
51 :			if (predators[i]->getType() != STOCKPREDATOR)
52 :			handle.logMessage(LOGFAIL, "Error in predatoraggregator - predator is not age structured", predators[i]->getName());
53 :
54 :			minage = ((StockPredator*)predators[i])->minAge();
55 :			maxage = ((StockPredator*)predators[i])->maxAge();
56 :			predConv.AddRows(1, maxage + 1, -1);
57 :			for (j = minage; j <= maxage; j++)
58 :			for (k = 0; k < predAges.Nrow(); k++)
59 :			for (l = 0; l < predAges.Ncol(k); l++)
60 :			if (j == predAges[k][l])
61 :			predConv[i][j] = k;
62 :			}
63 :
64 :			for (i = 0; i < preys.Size(); i++) {
65 :			if (!checkLengthGroupStructure(preys[i]->getLengthGroupDiv(), preyLgrpDiv))
66 :			handle.logMessage(LOGFAIL, "Error in predatoraggregator - invalid prey length group structure");
67 :			preyConv.AddRows(1, preys[i]->getLengthGroupDiv()->numLengthGroups(), -1);
68 :			for (j = 0; j < preyConv.Ncol(i); j++)
69 :			preyConv[i][j] = preyLgrpDiv->numLengthGroup(preys[i]->getLengthGroupDiv()->meanLength(j));
70 :			}
71 :
72 :			for (i = 0; i < predators.Size(); i++)
73 :			for (j = 0; j < preys.Size(); j++)
74 :			if (predators[i]->doesEat(preys[j]->getName()))
75 :			doeseat[i][j] = 1;
76 :
77 :			for (i = 0; i < areas.Nrow(); i++)
78 :			total.resize(new DoubleMatrix(predAges.Nrow(), preyLgrpDiv->numLengthGroups(), 0.0));
79 :			}
80 :
81 :			PredatorAggregator::~PredatorAggregator() {
82 :			int i;
83 :			for (i = 0; i < total.Size(); i++)
84 :			delete total[i];
85 :			}
86 :
87 :			void PredatorAggregator::Print(ofstream& outfile) const {
88 :			int i, j, k;
89 :			for (i = 0; i < total.Size(); i++) {
90 :			outfile << "\t\tInternal areas " << i << endl;
91 :			for (j = 0; j < total[i]->Nrow(); j++) {
92 :			outfile << TAB << TAB;
93 :			for (k = 0; k < total[i]->Ncol(j); k++)
94 :			outfile << setw(smallwidth) << (*total[i])[j][k] << sep;
95 :			outfile << endl;
96 :			}
97 :			}
98 :			outfile.flush();
99 :			}
100 :
101 :			void PredatorAggregator::Reset() {
102 :			int i;
103 :			for (i = 0; i < total.Size(); i++)
104 :			(*total[i]).setToZero();
105 :			}
106 :
107 :			void PredatorAggregator::Sum() {
108 :			int g, h, i, j, k, l, m;
109 :			double consum;
110 :			DoubleVector agesum, lensum;
111 :
112 :			this->Reset();
113 :			//sum over the appropriate preys, predators, areas and lengths
114 :			for (g = 0; g < predators.Size(); g++) {
115 :			for (h = 0; h < preys.Size(); h++) {
116 :			if (doeseat[g][h]) {
117 :			for (l = 0; l < areas.Nrow(); l++) {
118 :			for (j = 0; j < areas.Ncol(l); j++) {
119 :			if (predators[g]->isInArea(areas[l][j]) && preys[h]->isPreyArea(areas[l][j])) {
120 :			dptr = &predators[g]->getConsumption(areas[l][j], preys[h]->getName());
121 :
122 :			if (usepredages) {
123 :			//need to convert from length groups to age groups
124 :			alk = &((StockPredator*)predators[g])->getCurrentALK(areas[l][j]);
125 :
126 :			//first calculate how many predators there are in each age and length group
127 :			agesum.Reset();
128 :			lensum.Reset();
129 :			agesum.resize(alk->maxAge() + 1, 0.0);
130 :			lensum.resize(predators[g]->getLengthGroupDiv()->numLengthGroups(), 0.0);
131 :			for (k = alk->minAge(); k <= alk->maxAge(); k++) {
132 :			for (m = alk->minLength(k); m < alk->maxLength(k); m++) {
133 :			agesum[k] += (*alk)[k][m].N;
134 :			lensum[m] += (*alk)[k][m].N;
135 :			}
136 :			}
137 :
138 :			//then calculate the total consumption by the predators
139 :			for (k = alk->minAge(); k <= alk->maxAge(); k++) {
140 :			if (predConv[g][k] >= 0) {
141 :			for (i = 0; i < dptr->Ncol(k); i++) {
142 :			if (preyConv[h][i] >= 0) {
143 :			consum = 0.0;
144 :			for (m = 0; m < dptr->Nrow(); m++)
145 :			if (!isZero(lensum[m]))
146 :			consum += (*dptr)[m][i] * (*alk)[k][m].N / lensum[m];
147 :
148 :			(*total[l])[predConv[g][k]][preyConv[h][i]] += consum;
149 :			}
150 :			}
151 :			}
152 :			}
153 :
154 :			//finally convert this total consumption to a per-predator consumption
155 :			for (k = alk->minAge(); k <= alk->maxAge(); k++)
156 :			if ((predConv[g][k] >= 0) && (!isZero(agesum[k])))
157 :			for (i = 0; i < (*total[l])[predConv[g][k]].Size(); i++)
158 :			(*total[l])[predConv[g][k]][i] /= agesum[k];
159 :
160 :			} else {
161 :			for (k = 0; k < dptr->Nrow(); k++)
162 :			if (predConv[g][k] >= 0)
163 :			for (i = 0; i < dptr->Ncol(k); i++)
164 :			if (preyConv[h][i] >= 0)
165 :			(*total[l])[predConv[g][k]][preyConv[h][i]] += (*dptr)[k][i];
166 :
167 :			}
168 :			}
169 :			}
170 :			}
171 :			}
172 :			}
173 :			}
174 :			}
175 :
176 :			void PredatorAggregator::NumberSum() {
177 :			if (usepredages)
178 :			handle.logMessage(LOGFAIL, "Error in predatoraggregator - cannot sum numbers for age structured predators");
179 :
180 :			int g, h, i, j, k, l;
181 :			const PopInfoVector* preymeanw;
182 :
183 :			this->Reset();
184 :			for (g = 0; g < predators.Size(); g++) {
185 :			for (h = 0; h < preys.Size(); h++) {
186 :			if (doeseat[g][h]) {
187 :			for (l = 0; l < areas.Nrow(); l++) {
188 :			for (j = 0; j < areas.Ncol(l); j++) {
189 :			if (predators[g]->isInArea(areas[l][j]) && preys[h]->isPreyArea(areas[l][j])) {
190 :			dptr = &predators[g]->getConsumption(areas[l][j], preys[h]->getName());
191 :			preymeanw = &predators[g]->getConsumptionPopInfo(areas[l][j], preys[h]->getName());
192 :			for (k = 0; k < dptr->Nrow(); k++)
193 :			if (predConv[g][k] >= 0)
194 :			for (i = 0; i < dptr->Ncol(k); i++)
195 :			if (preyConv[h][i] >= 0 && (!(isZero((*preymeanw)[i].W))))
196 :			(*total[l])[predConv[g][k]][preyConv[h][i]] += (*dptr)[k][i] / (*preymeanw)[i].W;
197 :
198 :			}
199 :			}
200 :			}
201 :			}
202 :			}
203 :			}
204 :			}

---

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