Annotation of /trunk/gadget/grower.cc

Parent Directory | Revision Log

---

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

---

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