Annotation of /trunk/gadget/grower.cc

Parent Directory | Revision Log

---

Revision 2 - (view) (download)

1 :	agomez	1	#include "grower.h"
2 :			#include "errorhandler.h"
3 :			#include "readfunc.h"
4 :			#include "readword.h"
5 :			#include "keeper.h"
6 :			#include "areatime.h"
7 :			#include "growthcalc.h"
8 :			#include "gadget.h"
9 :			#include "global.h"
10 :
11 :			Grower::Grower(CommentStream& infile, const LengthGroupDivision* const OtherLgrpDiv,
12 :			const LengthGroupDivision* const GivenLgrpDiv, const IntVector& Areas,
13 :			const TimeClass* const TimeInfo, Keeper* const keeper, const char* refWeight,
14 :			const char* givenname, const AreaClass* const Area, const CharPtrVector& lenindex)
15 :			: HasName(givenname), LivesOnAreas(Areas), LgrpDiv(0), CI(0), growthcalc(0) {
16 :
17 :			char text[MaxStrLength];
18 :			strncpy(text, "", MaxStrLength);
19 :			int i;
20 :
21 :			keeper->addString("grower");
22 :			fixedweights = 0;
23 :			functionnumber = 0;
24 :			LgrpDiv = new LengthGroupDivision(*GivenLgrpDiv);
25 :	ulcessvp	2	vector_OK = 0;
26 :	agomez	1	if (LgrpDiv->Error())
27 :			handle.logMessage(LOGFAIL, "Error in grower - failed to create length group");
28 :			CI = new ConversionIndex(OtherLgrpDiv, LgrpDiv, 1);
29 :			if (CI->Error())
30 :			handle.logMessage(LOGFAIL, "Error in grower - error when checking length structure");
31 :
32 :			readWordAndValue(infile, "growthfunction", text);
33 :			if (strcasecmp(text, "multspec") == 0) {
34 :			functionnumber = 1;
35 :			growthcalc = new GrowthCalcA(infile, areas, TimeInfo, keeper);
36 :			} else if (strcasecmp(text, "fromfile") == 0) {
37 :			functionnumber = 2;
38 :			growthcalc = new GrowthCalcB(infile, areas, TimeInfo, keeper, Area, lenindex);
39 :			} else if (strcasecmp(text, "weightvb") == 0) {
40 :			functionnumber = 3;
41 :			growthcalc = new GrowthCalcC(infile, areas, TimeInfo, LgrpDiv, keeper, refWeight);
42 :			} else if (strcasecmp(text, "weightjones") == 0) {
43 :			functionnumber = 4;
44 :			growthcalc = new GrowthCalcD(infile, areas, TimeInfo, LgrpDiv, keeper, refWeight);
45 :			} else if (strcasecmp(text, "weightvbexpanded") == 0) {
46 :			functionnumber = 5;
47 :			growthcalc = new GrowthCalcE(infile, areas, TimeInfo, LgrpDiv, keeper, refWeight);
48 :			} else if (strcasecmp(text, "lengthvb") == 0) {
49 :			functionnumber = 6;
50 :			fixedweights = 1;
51 :			growthcalc = new GrowthCalcF(infile, areas, TimeInfo, keeper, Area, lenindex);
52 :			} else if (strcasecmp(text, "lengthpower") == 0) {
53 :			functionnumber = 7;
54 :			fixedweights = 1;
55 :			growthcalc = new GrowthCalcG(infile, areas, TimeInfo, keeper, Area, lenindex);
56 :			} else if (strcasecmp(text, "lengthvbsimple") == 0) {
57 :			functionnumber = 8;
58 :			growthcalc = new GrowthCalcH(infile, areas, TimeInfo, keeper);
59 :			} else if (strcasecmp(text, "weightjonessimple") == 0) {
60 :			functionnumber = 9;
61 :			growthcalc = new GrowthCalcI(infile, areas, TimeInfo, keeper);
62 :			} else if (strcasecmp(text, "lengthvpsimplet0") == 0) {
63 :			functionnumber = 10;
64 :			growthcalc = new GrowthCalcJ(infile, areas, TimeInfo, keeper);
65 :			} else if (strcasecmp(text, "lengthgompertz") == 0) {
66 :			functionnumber = 11;
67 :			growthcalc = new GrowthCalcK(infile, areas, TimeInfo, keeper);
68 :			} else {
69 :			handle.logFileMessage(LOGFAIL, "unrecognised growth function", text);
70 :			}
71 :
72 :			infile >> ws >> text;
73 :			if (strcasecmp(text, "beta") == 0) {
74 :			//Beta binomial growth distribution code is used
75 :			infile >> beta;
76 :			beta.Inform(keeper);
77 :			readWordAndVariable(infile, "maxlengthgroupgrowth", maxlengthgroupgrowth);
78 :
79 :			if (maxlengthgroupgrowth > LgrpDiv->numLengthGroups()) {
80 :			handle.logMessage(LOGWARN, "Warning in grower - maxlengthgroupgrowth is greater than number of length groups");
81 :			maxlengthgroupgrowth = LgrpDiv->numLengthGroups();
82 :			}
83 :
84 :			part1.resize(maxlengthgroupgrowth + 1, 0.0);
85 :			part2.resize(maxlengthgroupgrowth + 1, 0.0);
86 :			part4.resize(maxlengthgroupgrowth + 1, 0.0);
87 :
88 :			} else if (strcasecmp(text, "meanvarianceparameters") == 0) {
89 :			handle.logFileMessage(LOGFAIL, "\nThe mean variance parameters implementation of the growth is no longer supported\nUse the beta-binomial distribution implementation of the growth instead");
90 :
91 :			} else
92 :			handle.logFileUnexpected(LOGFAIL, "beta", text);
93 :
94 :			//Finished reading from input files.
95 :			keeper->clearLast();
96 :			int noareas = areas.Size();
97 :			int len = LgrpDiv->numLengthGroups();
98 :			int otherlen = OtherLgrpDiv->numLengthGroups();
99 :			PopInfo nullpop;
100 :			numGrow.AddRows(noareas, len, nullpop);
101 :
102 :			//setting storage spaces for growth
103 :	ulcessvp	2	// calcLengthGrowth.AddRows(noareas, len, 0.0);
104 :			// calcWeightGrowth.AddRows(noareas, len, 0.0);
105 :			_calcLengthGrowth.AddRows(noareas, len, 0.0);
106 :			_calcWeightGrowth.AddRows(noareas, len, 0.0);
107 :	agomez	1	interpLengthGrowth.AddRows(noareas, otherlen, 0.0);
108 :			interpWeightGrowth.AddRows(noareas, otherlen, 0.0);
109 :			dummyfphi.resize(len, 0.0);
110 :	ulcessvp	2	_lgrowth = new Matrix[noareas];
111 :			_wgrowth = new Matrix[noareas];
112 :
113 :	agomez	1	for (i = 0; i < noareas; i++) {
114 :	ulcessvp	2	// lgrowth.resize(new DoubleMatrix(maxlengthgroupgrowth + 1, otherlen, 0.0));
115 :			// wgrowth.resize(new DoubleMatrix(maxlengthgroupgrowth + 1, otherlen, 0.0));
116 :			_lgrowth[i] = *new Matrix(maxlengthgroupgrowth + 1, otherlen, 0.0);
117 :			_wgrowth[i] = *new Matrix(maxlengthgroupgrowth + 1, otherlen, 0.0);
118 :	agomez	1	}
119 :			}
120 :
121 :			Grower::~Grower() {
122 :			int i;
123 :	ulcessvp	2	//TODO ojo
124 :			int size = sizeof(_lgrowth)/sizeof(_lgrowth[0]);
125 :			for (i = 0; i < size; i++) {
126 :			delete &_lgrowth[i];
127 :			delete &_wgrowth[i];
128 :	agomez	1	}
129 :			delete CI;
130 :			delete LgrpDiv;
131 :			delete growthcalc;
132 :			}
133 :
134 :			void Grower::Print(ofstream& outfile) const {
135 :			int i, j, area;
136 :
137 :			outfile << "\nGrower\n\t";
138 :			LgrpDiv->Print(outfile);
139 :			for (area = 0; area < areas.Size(); area++) {
140 :			outfile << "\tLength increase on internal area " << areas[area] << ":\n\t";
141 :	ulcessvp	2	for (i = 0; i < _calcLengthGrowth.Ncol(); i++)
142 :			outfile << sep << _calcLengthGrowth[area][i];
143 :	agomez	1	outfile << "\n\tWeight increase on internal area " << areas[area] << ":\n\t";
144 :	ulcessvp	2	for (i = 0; i < _calcWeightGrowth.Ncol(); i++)
145 :			outfile << sep << _calcWeightGrowth[area][i];
146 :	agomez	1	outfile << "\n\tDistributed length increase on internal area " << areas[area] << ":\n";
147 :	ulcessvp	2	for (i = 0; i < _lgrowth[area].Nrow(); i++) {
148 :	agomez	1	outfile << TAB;
149 :	ulcessvp	2	for (j = 0; j < _lgrowth[area].Ncol(); j++)
150 :			outfile << sep << _lgrowth[area][i][j];
151 :	agomez	1	outfile << endl;
152 :			}
153 :			outfile << "\tDistributed weight increase on internal area " << areas[area] << ":\n";
154 :	ulcessvp	2	for (i = 0; i < _wgrowth[area].Nrow(); i++) {
155 :	agomez	1	outfile << TAB;
156 :	ulcessvp	2	for (j = 0; j < _wgrowth[area].Ncol(); j++)
157 :			outfile << sep << _wgrowth[area][i][j];
158 :	agomez	1	outfile << endl;
159 :			}
160 :			}
161 :			}
162 :
163 :			void Grower::Sum(const PopInfoVector& NumberInArea, int area) {
164 :			numGrow[this->areaNum(area)].Sum(&NumberInArea, *CI);
165 :			}
166 :
167 :			void Grower::calcGrowth(int area,
168 :			const AreaClass* const Area, const TimeClass* const TimeInfo) {
169 :
170 :			this->calcGrowth(area, Area, TimeInfo, dummyfphi, dummyfphi);
171 :			}
172 :
173 :			void Grower::calcGrowth(int area,
174 :			const AreaClass* const Area, const TimeClass* const TimeInfo,
175 :			const DoubleVector& FPhi, const DoubleVector& MaxCon) {
176 :
177 :			int inarea = this->areaNum(area);
178 :	ulcessvp	2	// growthcalc->calcGrowth(area, calcLengthGrowth[inarea], calcWeightGrowth[inarea],
179 :			// numGrow[inarea], Area, TimeInfo, FPhi, MaxCon, LgrpDiv);
180 :	agomez	1
181 :	ulcessvp	2	// cout << "size::" << _calcLengthGrowth.Ncol() << endl;
182 :			growthcalc->calcGrowth(area, _calcLengthGrowth[inarea], _calcWeightGrowth[inarea],
183 :			numGrow[inarea], Area, TimeInfo, FPhi, MaxCon, LgrpDiv, _calcLengthGrowth.Ncol());
184 :			// ofstream outfile;
185 :			// outfile.open("aaa", ios::out);
186 :			// _calcLengthGrowth.Print(outfile);
187 :			// int a;
188 :			// cout << "2-----------------" << endl;
189 :			// cin >> a;
190 :			//CI->interpolateLengths(interpLengthGrowth[inarea], calcLengthGrowth[inarea]);
191 :			CI->interpolateLengths(interpLengthGrowth[inarea], _calcLengthGrowth[inarea], _calcLengthGrowth.Ncol());
192 :	agomez	1	switch (functionnumber) {
193 :			case 1:
194 :			case 2:
195 :			case 3:
196 :			case 4:
197 :			case 5:
198 :			case 6:
199 :			case 7:
200 :			case 9:
201 :	ulcessvp	2	// CI->interpolateLengths(interpWeightGrowth[inarea], calcWeightGrowth[inarea]);
202 :			CI->interpolateLengths(interpWeightGrowth[inarea], _calcWeightGrowth[inarea], _calcWeightGrowth.Ncol());
203 :	agomez	1	break;
204 :			case 8:
205 :			case 10:
206 :			case 11:
207 :			break;
208 :			default:
209 :			handle.logMessage(LOGFAIL, "Error in grower - unrecognised growth function", functionnumber);
210 :			break;
211 :			}
212 :	ulcessvp	2	// ofstream outfile1;
213 :			// outfile1.open("aaa", ios::out);
214 :			// _calcLengthGrowth.Print(outfile1);
215 :			// cout << "-----------------" << endl;
216 :			// cin >> a;
217 :	agomez	1	}
218 :
219 :			void Grower::Reset() {
220 :			int i, j, area;
221 :			double factorialx, tmppart, tmpmax;
222 :
223 :	ulcessvp	2	// calcLengthGrowth.setToZero();
224 :			// calcWeightGrowth.setToZero();
225 :			_calcLengthGrowth.setToZero();
226 :			_calcWeightGrowth.setToZero();
227 :	agomez	1	interpLengthGrowth.setToZero();
228 :			for (area = 0; area < areas.Size(); area++) {
229 :	ulcessvp	2	(_lgrowth[area]).setToZero();
230 :	agomez	1	for (i = 0; i < LgrpDiv->numLengthGroups(); i++)
231 :			numGrow[area][i].setToZero();
232 :			}
233 :
234 :			switch (functionnumber) {
235 :			case 1:
236 :			case 2:
237 :			case 3:
238 :			case 4:
239 :			case 5:
240 :			case 6:
241 :			case 7:
242 :			case 9:
243 :			interpWeightGrowth.setToZero();
244 :			for (area = 0; area < areas.Size(); area++)
245 :	ulcessvp	2	(_wgrowth[area]).setToZero();
246 :	agomez	1	break;
247 :			case 8:
248 :			case 10:
249 :			case 11:
250 :			break;
251 :			default:
252 :			handle.logMessage(LOGFAIL, "Error in grower - unrecognised growth function", functionnumber);
253 :			break;
254 :			}
255 :
256 :			tmpmax = double(maxlengthgroupgrowth);
257 :			part1[0] = 1.0;
258 :			part1[1] = tmpmax;
259 :			factorialx = 1.0;
260 :			tmppart = tmpmax;
261 :			if (maxlengthgroupgrowth > 1) {
262 :			for (i = 2; i < maxlengthgroupgrowth + 1; i++) {
263 :			tmppart *= (tmpmax - i + 1);
264 :			factorialx *= double(i);
265 :			part1[i] = (tmppart / factorialx);
266 :			}
267 :			}
268 :
269 :			part2[maxlengthgroupgrowth] = 1.0;
270 :			part2[maxlengthgroupgrowth - 1] = beta;
271 :			if (maxlengthgroupgrowth > 1)
272 :			for (i = maxlengthgroupgrowth - 2; i >= 0; i--)
273 :			part2[i] = part2[i + 1] * (beta + tmpmax - i - 1);
274 :
275 :			//JMB this will never change so we can set it once
276 :			part4[0] = 1.0;
277 :
278 :			if (handle.getLogLevel() >= LOGMESSAGE)
279 :			handle.logMessage(LOGMESSAGE, "Reset grower data for stock", this->getName());
280 :			}

---

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