Annotation of /trunk/gadget/grower.cc

Parent Directory | Revision Log

---

Revision 4 - (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	4	calcLengthGrowth.AddRows(noareas, len, 0.0);
104 :			calcWeightGrowth.AddRows(noareas, len, 0.0);
105 :	agomez	1	interpLengthGrowth.AddRows(noareas, otherlen, 0.0);
106 :			interpWeightGrowth.AddRows(noareas, otherlen, 0.0);
107 :			dummyfphi.resize(len, 0.0);
108 :	ulcessvp	4	lgrowth = new Matrix[noareas];
109 :			wgrowth = new Matrix[noareas];
110 :	ulcessvp	2
111 :	agomez	1	for (i = 0; i < noareas; i++) {
112 :	ulcessvp	4	lgrowth[i] = *new Matrix(maxlengthgroupgrowth + 1, otherlen, 0.0);
113 :			wgrowth[i] = *new Matrix(maxlengthgroupgrowth + 1, otherlen, 0.0);
114 :	agomez	1	}
115 :			}
116 :
117 :			Grower::~Grower() {
118 :			int i;
119 :	ulcessvp	4	int size = sizeof(lgrowth)/sizeof(lgrowth[0]);
120 :	ulcessvp	2	for (i = 0; i < size; i++) {
121 :	ulcessvp	4	delete &lgrowth[i];
122 :			delete &wgrowth[i];
123 :	agomez	1	}
124 :			delete CI;
125 :			delete LgrpDiv;
126 :			delete growthcalc;
127 :			}
128 :
129 :			void Grower::Print(ofstream& outfile) const {
130 :			int i, j, area;
131 :
132 :			outfile << "\nGrower\n\t";
133 :			LgrpDiv->Print(outfile);
134 :			for (area = 0; area < areas.Size(); area++) {
135 :			outfile << "\tLength increase on internal area " << areas[area] << ":\n\t";
136 :	ulcessvp	4	for (i = 0; i < calcLengthGrowth.Ncol(); i++)
137 :			outfile << sep << calcLengthGrowth[area][i];
138 :	agomez	1	outfile << "\n\tWeight increase on internal area " << areas[area] << ":\n\t";
139 :	ulcessvp	4	for (i = 0; i < calcWeightGrowth.Ncol(); i++)
140 :			outfile << sep << calcWeightGrowth[area][i];
141 :	agomez	1	outfile << "\n\tDistributed length increase on internal area " << areas[area] << ":\n";
142 :	ulcessvp	4	for (i = 0; i < lgrowth[area].Nrow(); i++) {
143 :	agomez	1	outfile << TAB;
144 :	ulcessvp	4	for (j = 0; j < lgrowth[area].Ncol(); j++)
145 :			outfile << sep << lgrowth[area][i][j];
146 :	agomez	1	outfile << endl;
147 :			}
148 :			outfile << "\tDistributed weight increase on internal area " << areas[area] << ":\n";
149 :	ulcessvp	4	for (i = 0; i < wgrowth[area].Nrow(); i++) {
150 :	agomez	1	outfile << TAB;
151 :	ulcessvp	4	for (j = 0; j < wgrowth[area].Ncol(); j++)
152 :			outfile << sep << wgrowth[area][i][j];
153 :	agomez	1	outfile << endl;
154 :			}
155 :			}
156 :			}
157 :
158 :			void Grower::Sum(const PopInfoVector& NumberInArea, int area) {
159 :			numGrow[this->areaNum(area)].Sum(&NumberInArea, *CI);
160 :			}
161 :
162 :			void Grower::calcGrowth(int area,
163 :			const AreaClass* const Area, const TimeClass* const TimeInfo) {
164 :			this->calcGrowth(area, Area, TimeInfo, dummyfphi, dummyfphi);
165 :			}
166 :
167 :			void Grower::calcGrowth(int area,
168 :			const AreaClass* const Area, const TimeClass* const TimeInfo,
169 :			const DoubleVector& FPhi, const DoubleVector& MaxCon) {
170 :
171 :			int inarea = this->areaNum(area);
172 :
173 :	ulcessvp	4	growthcalc->calcGrowth(area, calcLengthGrowth[inarea], calcWeightGrowth[inarea],
174 :			numGrow[inarea], Area, TimeInfo, FPhi, MaxCon, LgrpDiv, calcLengthGrowth.Ncol());
175 :
176 :			CI->interpolateLengths(interpLengthGrowth[inarea], calcLengthGrowth[inarea], calcLengthGrowth.Ncol());
177 :	agomez	1	switch (functionnumber) {
178 :			case 1:
179 :			case 2:
180 :			case 3:
181 :			case 4:
182 :			case 5:
183 :			case 6:
184 :			case 7:
185 :			case 9:
186 :	ulcessvp	4	CI->interpolateLengths(interpWeightGrowth[inarea], calcWeightGrowth[inarea], calcWeightGrowth.Ncol());
187 :	agomez	1	break;
188 :			case 8:
189 :			case 10:
190 :			case 11:
191 :			break;
192 :			default:
193 :			handle.logMessage(LOGFAIL, "Error in grower - unrecognised growth function", functionnumber);
194 :			break;
195 :			}
196 :			}
197 :
198 :			void Grower::Reset() {
199 :			int i, j, area;
200 :			double factorialx, tmppart, tmpmax;
201 :
202 :	ulcessvp	4	calcLengthGrowth.setToZero();
203 :			calcWeightGrowth.setToZero();
204 :	agomez	1	interpLengthGrowth.setToZero();
205 :			for (area = 0; area < areas.Size(); area++) {
206 :	ulcessvp	4	(lgrowth[area]).setToZero();
207 :	agomez	1	for (i = 0; i < LgrpDiv->numLengthGroups(); i++)
208 :			numGrow[area][i].setToZero();
209 :			}
210 :
211 :			switch (functionnumber) {
212 :			case 1:
213 :			case 2:
214 :			case 3:
215 :			case 4:
216 :			case 5:
217 :			case 6:
218 :			case 7:
219 :			case 9:
220 :			interpWeightGrowth.setToZero();
221 :			for (area = 0; area < areas.Size(); area++)
222 :	ulcessvp	4	(wgrowth[area]).setToZero();
223 :	agomez	1	break;
224 :			case 8:
225 :			case 10:
226 :			case 11:
227 :			break;
228 :			default:
229 :			handle.logMessage(LOGFAIL, "Error in grower - unrecognised growth function", functionnumber);
230 :			break;
231 :			}
232 :
233 :			tmpmax = double(maxlengthgroupgrowth);
234 :			part1[0] = 1.0;
235 :			part1[1] = tmpmax;
236 :			factorialx = 1.0;
237 :			tmppart = tmpmax;
238 :			if (maxlengthgroupgrowth > 1) {
239 :			for (i = 2; i < maxlengthgroupgrowth + 1; i++) {
240 :			tmppart *= (tmpmax - i + 1);
241 :			factorialx *= double(i);
242 :			part1[i] = (tmppart / factorialx);
243 :			}
244 :			}
245 :
246 :			part2[maxlengthgroupgrowth] = 1.0;
247 :			part2[maxlengthgroupgrowth - 1] = beta;
248 :			if (maxlengthgroupgrowth > 1)
249 :			for (i = maxlengthgroupgrowth - 2; i >= 0; i--)
250 :			part2[i] = part2[i + 1] * (beta + tmpmax - i - 1);
251 :
252 :			//JMB this will never change so we can set it once
253 :			part4[0] = 1.0;
254 :
255 :			if (handle.getLogLevel() >= LOGMESSAGE)
256 :			handle.logMessage(LOGMESSAGE, "Reset grower data for stock", this->getName());
257 :			}

---

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