1 : |
agomez |
1 |
#include "renewal.h" |
2 : |
|
|
#include "errorhandler.h" |
3 : |
|
|
#include "readfunc.h" |
4 : |
|
|
#include "keeper.h" |
5 : |
|
|
#include "readword.h" |
6 : |
|
|
#include "gadget.h" |
7 : |
|
|
#include "global.h" |
8 : |
|
|
|
9 : |
|
|
RenewalData::RenewalData(CommentStream& infile, const IntVector& Areas,
|
10 : |
|
|
const AreaClass* const Area, const TimeClass* const TimeInfo, Keeper* const keeper,
|
11 : |
|
|
const char* refWeightFile, const char* givenname, int minage, int maxage, double DL)
|
12 : |
|
|
: HasName(givenname), LivesOnAreas(Areas), CI(0), LgrpDiv(0) {
|
13 : |
|
|
|
14 : |
|
|
keeper->addString("renewaldata");
|
15 : |
|
|
ifstream subfile;
|
16 : |
|
|
CommentStream subcomment(subfile);
|
17 : |
|
|
char text[MaxStrLength];
|
18 : |
|
|
strncpy(text, "", MaxStrLength);
|
19 : |
|
|
|
20 : |
|
|
index = 0;
|
21 : |
|
|
int i, j;
|
22 : |
|
|
double minlength, maxlength, dl;
|
23 : |
|
|
readWordAndVariable(infile, "minlength", minlength);
|
24 : |
|
|
readWordAndVariable(infile, "maxlength", maxlength);
|
25 : |
|
|
|
26 : |
|
|
//JMB - changed to make the reading of dl optional
|
27 : |
|
|
//if it isnt specifed here, it will default to the dl value of the stock
|
28 : |
|
|
infile >> ws;
|
29 : |
|
|
char c = infile.peek();
|
30 : |
|
|
if ((c == 'd') || (c == 'D'))
|
31 : |
|
|
readWordAndVariable(infile, "dl", dl);
|
32 : |
|
|
else
|
33 : |
|
|
dl = DL;
|
34 : |
|
|
|
35 : |
|
|
LgrpDiv = new LengthGroupDivision(minlength, maxlength, dl);
|
36 : |
|
|
if (LgrpDiv->Error())
|
37 : |
|
|
handle.logMessage(LOGFAIL, "Error in renewal - failed to create length group");
|
38 : |
|
|
|
39 : |
|
|
infile >> text >> ws;
|
40 : |
|
|
if ((strcasecmp(text, "normalcondfile") == 0)) {
|
41 : |
|
|
//read renewal data in mean length format, using the reference weight file
|
42 : |
|
|
readoption = 0;
|
43 : |
|
|
infile >> text >> ws;
|
44 : |
|
|
subfile.open(text, ios::in);
|
45 : |
|
|
handle.checkIfFailure(subfile, text);
|
46 : |
|
|
handle.Open(text);
|
47 : |
|
|
this->readNormalConditionData(subcomment, keeper, TimeInfo, Area, minage, maxage);
|
48 : |
|
|
handle.Close();
|
49 : |
|
|
subfile.close();
|
50 : |
|
|
subfile.clear();
|
51 : |
|
|
|
52 : |
|
|
//read information on reference weights.
|
53 : |
|
|
DoubleMatrix tmpRefW;
|
54 : |
|
|
keeper->addString("referenceweights");
|
55 : |
|
|
subfile.open(refWeightFile, ios::in);
|
56 : |
|
|
handle.checkIfFailure(subfile, refWeightFile);
|
57 : |
|
|
handle.Open(refWeightFile);
|
58 : |
|
|
readRefWeights(subcomment, tmpRefW);
|
59 : |
|
|
handle.Close();
|
60 : |
|
|
subfile.close();
|
61 : |
|
|
subfile.clear();
|
62 : |
|
|
|
63 : |
|
|
//Interpolate the reference weights. First there are some error checks.
|
64 : |
|
|
if (LgrpDiv->meanLength(0) < tmpRefW[0][0] ||
|
65 : |
|
|
LgrpDiv->meanLength(LgrpDiv->numLengthGroups() - 1) > tmpRefW[tmpRefW.Nrow() - 1][0])
|
66 : |
|
|
handle.logFileMessage(LOGFAIL, "lengths for reference weights must span the range of initial condition lengths");
|
67 : |
|
|
|
68 : |
|
|
//Aggregate the reference weight data to be the same format
|
69 : |
|
|
double ratio, tmplen;
|
70 : |
|
|
refWeight.resize(LgrpDiv->numLengthGroups(), 0.0);
|
71 : |
|
|
int pos = 0;
|
72 : |
|
|
for (j = 0; j < LgrpDiv->numLengthGroups(); j++) {
|
73 : |
|
|
tmplen = LgrpDiv->meanLength(j);
|
74 : |
|
|
for (i = pos; i < tmpRefW.Nrow() - 1; i++) {
|
75 : |
|
|
if (((tmplen > tmpRefW[i][0]) || (isEqual(tmplen, tmpRefW[i][0]))) &&
|
76 : |
|
|
((tmplen < tmpRefW[i + 1][0]) || (isEqual(tmplen, tmpRefW[i + 1][0])))) {
|
77 : |
|
|
|
78 : |
|
|
ratio = (tmplen - tmpRefW[i][0]) / (tmpRefW[i + 1][0] - tmpRefW[i][0]);
|
79 : |
|
|
refWeight[j] = tmpRefW[i][1] + ratio * (tmpRefW[i + 1][1] - tmpRefW[i][1]);
|
80 : |
|
|
pos = i;
|
81 : |
|
|
}
|
82 : |
|
|
}
|
83 : |
|
|
}
|
84 : |
|
|
keeper->clearLast();
|
85 : |
|
|
|
86 : |
|
|
} else if ((strcasecmp(text, "renewaldatafile") == 0) || (strcasecmp(text, "normalparamfile") == 0)) {
|
87 : |
|
|
//read renewal data in mean length format, using a length weight relationship
|
88 : |
|
|
readoption = 1;
|
89 : |
|
|
infile >> text >> ws;
|
90 : |
|
|
subfile.open(text, ios::in);
|
91 : |
|
|
handle.checkIfFailure(subfile, text);
|
92 : |
|
|
handle.Open(text);
|
93 : |
|
|
this->readNormalParameterData(subcomment, keeper, TimeInfo, Area, minage, maxage);
|
94 : |
|
|
handle.Close();
|
95 : |
|
|
subfile.close();
|
96 : |
|
|
subfile.clear();
|
97 : |
|
|
|
98 : |
|
|
} else if ((strcasecmp(text, "numberfile") == 0)) {
|
99 : |
|
|
//read renewal data in number format
|
100 : |
|
|
readoption = 2;
|
101 : |
|
|
infile >> text >> ws;
|
102 : |
|
|
subfile.open(text, ios::in);
|
103 : |
|
|
handle.checkIfFailure(subfile, text);
|
104 : |
|
|
handle.Open(text);
|
105 : |
|
|
this->readNumberData(subcomment, keeper, TimeInfo, Area, minage, maxage);
|
106 : |
|
|
handle.Close();
|
107 : |
|
|
subfile.close();
|
108 : |
|
|
subfile.clear();
|
109 : |
|
|
|
110 : |
|
|
} else
|
111 : |
|
|
handle.logFileMessage(LOGFAIL, "unrecognised renewal data format", text);
|
112 : |
|
|
|
113 : |
|
|
keeper->clearLast();
|
114 : |
|
|
}
|
115 : |
|
|
|
116 : |
|
|
void RenewalData::readNormalParameterData(CommentStream& infile, Keeper* const keeper,
|
117 : |
|
|
const TimeClass* const TimeInfo, const AreaClass* const Area, int minage, int maxage) {
|
118 : |
|
|
|
119 : |
|
|
char c;
|
120 : |
|
|
int year, step, area, age, keepdata, count, reject;
|
121 : |
|
|
PopInfoIndexVector poptmp(LgrpDiv->numLengthGroups(), 0);
|
122 : |
|
|
|
123 : |
|
|
//We now expect to find the renewal data in the following format
|
124 : |
|
|
//year, step, area, age, number, mean, sdev, alpha, beta
|
125 : |
|
|
infile >> ws;
|
126 : |
|
|
if (countColumns(infile) != 9)
|
127 : |
|
|
handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 9");
|
128 : |
|
|
|
129 : |
|
|
year = step = area = age = count = reject = 0;
|
130 : |
|
|
while (!infile.eof()) {
|
131 : |
|
|
//crude check to see if something has gone wrong and avoid infinite loops
|
132 : |
|
|
if (!(isdigit(infile.peek())))
|
133 : |
|
|
handle.logFileMessage(LOGFAIL, "failed to read data from file");
|
134 : |
|
|
|
135 : |
|
|
keepdata = 1;
|
136 : |
|
|
infile >> year >> step >> area >> age >> ws;
|
137 : |
|
|
|
138 : |
|
|
if (!(TimeInfo->isWithinPeriod(year, step)))
|
139 : |
|
|
keepdata = 0;
|
140 : |
|
|
|
141 : |
|
|
if (!(this->isInArea(Area->getInnerArea(area))))
|
142 : |
|
|
keepdata = 0;
|
143 : |
|
|
|
144 : |
|
|
if ((age < minage) || (age > maxage))
|
145 : |
|
|
keepdata = 0;
|
146 : |
|
|
|
147 : |
|
|
if (keepdata == 1) {
|
148 : |
|
|
//renewal data is required, so store it
|
149 : |
|
|
renewalTime.resize(1, TimeInfo->calcSteps(year, step));
|
150 : |
|
|
renewalArea.resize(1, Area->getInnerArea(area));
|
151 : |
|
|
renewalAge.resize(1, age);
|
152 : |
|
|
|
153 : |
|
|
renewalMult.resize(1, keeper);
|
154 : |
|
|
meanLength.resize(1, keeper);
|
155 : |
|
|
sdevLength.resize(1, keeper);
|
156 : |
|
|
alpha.resize(1, keeper);
|
157 : |
|
|
beta.resize(1, keeper);
|
158 : |
|
|
renewalDistribution.resize(1, new AgeBandMatrix(age, poptmp));
|
159 : |
|
|
|
160 : |
|
|
infile >> renewalMult[count] >> ws;
|
161 : |
|
|
infile >> meanLength[count] >> ws;
|
162 : |
|
|
infile >> sdevLength[count] >> ws;
|
163 : |
|
|
infile >> alpha[count] >> ws;
|
164 : |
|
|
infile >> beta[count] >> ws;
|
165 : |
|
|
|
166 : |
|
|
renewalMult[count].Inform(keeper);
|
167 : |
|
|
meanLength[count].Inform(keeper);
|
168 : |
|
|
sdevLength[count].Inform(keeper);
|
169 : |
|
|
alpha[count].Inform(keeper);
|
170 : |
|
|
beta[count].Inform(keeper);
|
171 : |
|
|
count++;
|
172 : |
|
|
|
173 : |
|
|
} else { //renewal data not required - skip rest of line
|
174 : |
|
|
reject++;
|
175 : |
|
|
infile.get(c);
|
176 : |
|
|
while (c != '\n' && !infile.eof())
|
177 : |
|
|
infile.get(c);
|
178 : |
|
|
infile >> ws;
|
179 : |
|
|
}
|
180 : |
|
|
}
|
181 : |
|
|
|
182 : |
|
|
if (count == 0)
|
183 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - found no data in the data file");
|
184 : |
|
|
if (reject != 0)
|
185 : |
|
|
handle.logMessage(LOGMESSAGE, "Discarded invalid renewal data - number of invalid entries", reject);
|
186 : |
|
|
handle.logMessage(LOGMESSAGE, "Read renewal data file - number of entries", count);
|
187 : |
|
|
}
|
188 : |
|
|
|
189 : |
|
|
void RenewalData::readNormalConditionData(CommentStream& infile, Keeper* const keeper,
|
190 : |
|
|
const TimeClass* const TimeInfo, const AreaClass* const Area, int minage, int maxage) {
|
191 : |
|
|
|
192 : |
|
|
char c;
|
193 : |
|
|
int year, step, area, age, keepdata, count, reject;
|
194 : |
|
|
PopInfoIndexVector poptmp(LgrpDiv->numLengthGroups(), 0);
|
195 : |
|
|
|
196 : |
|
|
//We now expect to find the renewal data in the following format
|
197 : |
|
|
//year, step, area, age, number, mean, sdev, cond
|
198 : |
|
|
infile >> ws;
|
199 : |
|
|
if (countColumns(infile) != 8)
|
200 : |
|
|
handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 8");
|
201 : |
|
|
|
202 : |
|
|
year = step = area = age = count = reject = 0;
|
203 : |
|
|
while (!infile.eof()) {
|
204 : |
|
|
//crude check to see if something has gone wrong and avoid infinite loops
|
205 : |
|
|
if (!(isdigit(infile.peek())))
|
206 : |
|
|
handle.logFileMessage(LOGFAIL, "failed to read data from file");
|
207 : |
|
|
|
208 : |
|
|
keepdata = 1;
|
209 : |
|
|
infile >> year >> step >> area >> age >> ws;
|
210 : |
|
|
|
211 : |
|
|
if (!(TimeInfo->isWithinPeriod(year, step)))
|
212 : |
|
|
keepdata = 0;
|
213 : |
|
|
|
214 : |
|
|
if (!(this->isInArea(Area->getInnerArea(area))))
|
215 : |
|
|
keepdata = 0;
|
216 : |
|
|
|
217 : |
|
|
if ((age < minage) || (age > maxage))
|
218 : |
|
|
keepdata = 0;
|
219 : |
|
|
|
220 : |
|
|
if (keepdata == 1) {
|
221 : |
|
|
//renewal data is required, so store it
|
222 : |
|
|
renewalTime.resize(1, TimeInfo->calcSteps(year, step));
|
223 : |
|
|
renewalArea.resize(1, Area->getInnerArea(area));
|
224 : |
|
|
renewalAge.resize(1, age);
|
225 : |
|
|
|
226 : |
|
|
renewalMult.resize(1, keeper);
|
227 : |
|
|
meanLength.resize(1, keeper);
|
228 : |
|
|
sdevLength.resize(1, keeper);
|
229 : |
|
|
relCond.resize(1, keeper);
|
230 : |
|
|
renewalDistribution.resize(1, new AgeBandMatrix(age, poptmp));
|
231 : |
|
|
|
232 : |
|
|
infile >> renewalMult[count] >> ws;
|
233 : |
|
|
infile >> meanLength[count] >> ws;
|
234 : |
|
|
infile >> sdevLength[count] >> ws;
|
235 : |
|
|
infile >> relCond[count] >> ws;
|
236 : |
|
|
|
237 : |
|
|
renewalMult[count].Inform(keeper);
|
238 : |
|
|
meanLength[count].Inform(keeper);
|
239 : |
|
|
sdevLength[count].Inform(keeper);
|
240 : |
|
|
relCond[count].Inform(keeper);
|
241 : |
|
|
count++;
|
242 : |
|
|
|
243 : |
|
|
} else { //renewal data not required - skip rest of line
|
244 : |
|
|
reject++;
|
245 : |
|
|
infile.get(c);
|
246 : |
|
|
while (c != '\n' && !infile.eof())
|
247 : |
|
|
infile.get(c);
|
248 : |
|
|
infile >> ws;
|
249 : |
|
|
}
|
250 : |
|
|
}
|
251 : |
|
|
|
252 : |
|
|
if (count == 0)
|
253 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - found no data in the data file");
|
254 : |
|
|
if (reject != 0)
|
255 : |
|
|
handle.logMessage(LOGMESSAGE, "Discarded invalid renewal data - number of invalid entries", reject);
|
256 : |
|
|
handle.logMessage(LOGMESSAGE, "Read renewal data file - number of entries", count);
|
257 : |
|
|
}
|
258 : |
|
|
|
259 : |
|
|
void RenewalData::readNumberData(CommentStream& infile, Keeper* const keeper,
|
260 : |
|
|
const TimeClass* const TimeInfo, const AreaClass* const Area, int minage, int maxage) {
|
261 : |
|
|
|
262 : |
|
|
char c;
|
263 : |
|
|
double length;
|
264 : |
|
|
int i, year, step, area, age, keepdata, count, id, lengthid, reject;
|
265 : |
|
|
int numlen = LgrpDiv->numLengthGroups();
|
266 : |
|
|
PopInfoIndexVector poptmp(numlen, 0);
|
267 : |
|
|
|
268 : |
|
|
//We now expect to find the renewal data in the following format
|
269 : |
|
|
//year, step, area, age, length, number, mean weight
|
270 : |
|
|
infile >> ws;
|
271 : |
|
|
if (countColumns(infile) != 7)
|
272 : |
|
|
handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 7");
|
273 : |
|
|
|
274 : |
|
|
year = step = area = age = count = reject = 0;
|
275 : |
|
|
while (!infile.eof()) {
|
276 : |
|
|
//crude check to see if something has gone wrong and avoid infinite loops
|
277 : |
|
|
if (!(isdigit(infile.peek())))
|
278 : |
|
|
handle.logFileMessage(LOGFAIL, "failed to read data from file");
|
279 : |
|
|
|
280 : |
|
|
keepdata = 1;
|
281 : |
|
|
infile >> year >> step >> area >> age >> length >> ws;
|
282 : |
|
|
|
283 : |
|
|
if (!(TimeInfo->isWithinPeriod(year, step)))
|
284 : |
|
|
keepdata = 0;
|
285 : |
|
|
|
286 : |
|
|
if (!(this->isInArea(Area->getInnerArea(area))))
|
287 : |
|
|
keepdata = 0;
|
288 : |
|
|
|
289 : |
|
|
if ((age < minage) || (age > maxage))
|
290 : |
|
|
keepdata = 0;
|
291 : |
|
|
|
292 : |
|
|
lengthid = -1;
|
293 : |
|
|
for (i = 0; i < numlen; i++)
|
294 : |
|
|
if (isEqual(length, LgrpDiv->minLength(i)))
|
295 : |
|
|
lengthid = i;
|
296 : |
|
|
|
297 : |
|
|
//OK the length doesnt match a minimum length so find the length group it is in
|
298 : |
|
|
if ((lengthid == -1) && (length > LgrpDiv->minLength()) && (length < LgrpDiv->maxLength())) {
|
299 : |
|
|
for (i = 1; i < numlen; i++) {
|
300 : |
|
|
if (length < LgrpDiv->minLength(i)) {
|
301 : |
|
|
lengthid = i - 1;
|
302 : |
|
|
break;
|
303 : |
|
|
}
|
304 : |
|
|
}
|
305 : |
|
|
if (lengthid == -1)
|
306 : |
|
|
lengthid = numlen - 1; //then this must be the last length group
|
307 : |
|
|
}
|
308 : |
|
|
|
309 : |
|
|
if (lengthid == -1)
|
310 : |
|
|
keepdata = 0;
|
311 : |
|
|
|
312 : |
|
|
if (keepdata == 1) {
|
313 : |
|
|
//renewal data is required, so store it
|
314 : |
|
|
id = -1;
|
315 : |
|
|
for (i = 0; i < renewalTime.Size(); i++)
|
316 : |
|
|
if ((renewalTime[i] == TimeInfo->calcSteps(year, step))
|
317 : |
|
|
&& (renewalArea[i] == Area->getInnerArea(area))
|
318 : |
|
|
&& (renewalAge[i] == age))
|
319 : |
|
|
id = i;
|
320 : |
|
|
|
321 : |
|
|
if (id == -1) {
|
322 : |
|
|
//this is a new timestep/area combination
|
323 : |
|
|
renewalTime.resize(1, TimeInfo->calcSteps(year, step));
|
324 : |
|
|
renewalArea.resize(1, Area->getInnerArea(area));
|
325 : |
|
|
renewalAge.resize(1, age);
|
326 : |
|
|
id = renewalTime.Size() - 1;
|
327 : |
|
|
|
328 : |
|
|
renewalDistribution.resize(1, new AgeBandMatrix(age, poptmp));
|
329 : |
|
|
renewalNumber.resize(new FormulaMatrix(maxage - minage + 1, numlen, 0.0));
|
330 : |
|
|
}
|
331 : |
|
|
|
332 : |
|
|
renewalDistribution[id][age][lengthid].N = 0.0;
|
333 : |
|
|
infile >> (*renewalNumber[id])[age - minage][lengthid] >> ws;
|
334 : |
|
|
infile >> renewalDistribution[id][age][lengthid].W >> ws;
|
335 : |
|
|
count++;
|
336 : |
|
|
|
337 : |
|
|
} else { //renewal data not required - skip rest of line
|
338 : |
|
|
reject++;
|
339 : |
|
|
infile.get(c);
|
340 : |
|
|
while (c != '\n' && !infile.eof())
|
341 : |
|
|
infile.get(c);
|
342 : |
|
|
infile >> ws;
|
343 : |
|
|
}
|
344 : |
|
|
}
|
345 : |
|
|
|
346 : |
|
|
for (i = 0; i < renewalNumber.Size(); i++)
|
347 : |
|
|
(*renewalNumber[i]).Inform(keeper);
|
348 : |
|
|
|
349 : |
|
|
if (count == 0)
|
350 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - found no data in the data file");
|
351 : |
|
|
if (reject != 0)
|
352 : |
|
|
handle.logMessage(LOGMESSAGE, "Discarded invalid renewal data - number of invalid entries", reject);
|
353 : |
|
|
handle.logMessage(LOGMESSAGE, "Read renewal data file - number of entries", count);
|
354 : |
|
|
}
|
355 : |
|
|
|
356 : |
|
|
RenewalData::~RenewalData() {
|
357 : |
|
|
int i;
|
358 : |
|
|
delete LgrpDiv;
|
359 : |
|
|
delete CI;
|
360 : |
|
|
for (i = 0; i < renewalNumber.Size(); i++)
|
361 : |
|
|
delete renewalNumber[i];
|
362 : |
|
|
}
|
363 : |
|
|
|
364 : |
|
|
void RenewalData::setCI(const LengthGroupDivision* const GivenLDiv) {
|
365 : |
|
|
if (!checkLengthGroupStructure(GivenLDiv, LgrpDiv))
|
366 : |
|
|
handle.logMessage(LOGFAIL, "Error in renewal - invalid length group structure for stock", this->getName());
|
367 : |
|
|
//check minimum and maximum lengths
|
368 : |
|
|
if (LgrpDiv->minLength() < GivenLDiv->minLength())
|
369 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - minimum length less than stock length for stock", this->getName());
|
370 : |
|
|
if (LgrpDiv->maxLength() > GivenLDiv->maxLength())
|
371 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - maximum length greater than stock length for stock", this->getName());
|
372 : |
|
|
CI = new ConversionIndex(LgrpDiv, GivenLDiv);
|
373 : |
|
|
if (CI->Error())
|
374 : |
|
|
handle.logMessage(LOGFAIL, "Error in renewal - error when checking length structure");
|
375 : |
|
|
}
|
376 : |
|
|
|
377 : |
|
|
void RenewalData::Print(ofstream& outfile) const {
|
378 : |
|
|
int i;
|
379 : |
|
|
outfile << "\nRenewal data\n\t";
|
380 : |
|
|
LgrpDiv->Print(outfile);
|
381 : |
|
|
//need to find all areas and ages for current index value ...
|
382 : |
|
|
for (i = 0; i < renewalTime.Size(); i++) {
|
383 : |
|
|
if (renewalTime[i] == renewalTime[index]) {
|
384 : |
|
|
outfile << "\tInternal area " << renewalArea[i] << " age " << renewalAge[i];
|
385 : |
|
|
|
386 : |
|
|
if (readoption == 2)
|
387 : |
|
|
outfile << "\n\tNumbers\n";
|
388 : |
|
|
else
|
389 : |
|
|
outfile << " multiplier " << renewalMult[i] << "\n\tNumbers\n";
|
390 : |
|
|
|
391 : |
|
|
renewalDistribution[i].printNumbers(outfile);
|
392 : |
|
|
outfile << "\tMean weights\n";
|
393 : |
|
|
renewalDistribution[i].printWeights(outfile);
|
394 : |
|
|
}
|
395 : |
|
|
}
|
396 : |
|
|
outfile.flush();
|
397 : |
|
|
}
|
398 : |
|
|
|
399 : |
|
|
void RenewalData::Reset() {
|
400 : |
|
|
int i, age, l, minage;
|
401 : |
|
|
double sum, mult, dnorm;
|
402 : |
|
|
|
403 : |
|
|
index = 0;
|
404 : |
|
|
if (readoption == 0) {
|
405 : |
|
|
for (i = 0; i < renewalTime.Size(); i++) {
|
406 : |
|
|
age = renewalAge[i];
|
407 : |
|
|
|
408 : |
|
|
//JMB check that the length data is valid
|
409 : |
|
|
if (isZero(meanLength[i]) || sdevLength[i] < 0.04) {
|
410 : |
|
|
//JMB the limit has been set at 0.04 to keep the exponential calculation sane
|
411 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - invalid length data");
|
412 : |
|
|
|
413 : |
|
|
//JMB set the population to zero
|
414 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++)
|
415 : |
|
|
renewalDistribution[i][age][l].setToZero();
|
416 : |
|
|
|
417 : |
|
|
} else {
|
418 : |
|
|
//JMB check that the mean length is within the length group range
|
419 : |
|
|
if (meanLength[i] < LgrpDiv->minLength())
|
420 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - mean length is less than minimum length for stock", this->getName());
|
421 : |
|
|
if (meanLength[i] > LgrpDiv->maxLength())
|
422 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - mean length is greater than maximum length for stock", this->getName());
|
423 : |
|
|
|
424 : |
|
|
sum = 0.0;
|
425 : |
|
|
mult = 1.0 / sdevLength[i];
|
426 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++) {
|
427 : |
|
|
dnorm = (LgrpDiv->meanLength(l) - meanLength[i]) * mult;
|
428 : |
|
|
renewalDistribution[i][age][l].N = exp(-(dnorm * dnorm) * 0.5);
|
429 : |
|
|
sum += renewalDistribution[i][age][l].N;
|
430 : |
|
|
}
|
431 : |
|
|
|
432 : |
|
|
if (isZero(sum)) {
|
433 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - calculated zero recruits for stock", this->getName());
|
434 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++)
|
435 : |
|
|
renewalDistribution[i][age][l].setToZero();
|
436 : |
|
|
|
437 : |
|
|
} else {
|
438 : |
|
|
sum = 10000.0 / sum;
|
439 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++) {
|
440 : |
|
|
renewalDistribution[i][age][l].N *= sum;
|
441 : |
|
|
renewalDistribution[i][age][l].W = refWeight[l] * relCond[i];
|
442 : |
|
|
if ((handle.getLogLevel() >= LOGWARN) && (isZero(renewalDistribution[i][age][l].W)) && (renewalDistribution[i][age][l].N > 0.0))
|
443 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - zero mean weight for stock", this->getName());
|
444 : |
|
|
}
|
445 : |
|
|
}
|
446 : |
|
|
}
|
447 : |
|
|
}
|
448 : |
|
|
|
449 : |
|
|
} else if (readoption == 1) {
|
450 : |
|
|
for (i = 0; i < renewalTime.Size(); i++) {
|
451 : |
|
|
age = renewalAge[i];
|
452 : |
|
|
|
453 : |
|
|
//JMB check that the length data is valid
|
454 : |
|
|
if (isZero(meanLength[i]) || sdevLength[i] < 0.04) {
|
455 : |
|
|
//JMB the limit has been set at 0.04 to keep the exponential calculation sane
|
456 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - invalid length data");
|
457 : |
|
|
|
458 : |
|
|
//JMB set the population to zero
|
459 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++)
|
460 : |
|
|
renewalDistribution[i][age][l].setToZero();
|
461 : |
|
|
|
462 : |
|
|
} else {
|
463 : |
|
|
//JMB check that the mean length is within the length group range
|
464 : |
|
|
if (meanLength[i] < LgrpDiv->minLength())
|
465 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - mean length is less than minimum length for stock", this->getName());
|
466 : |
|
|
if (meanLength[i] > LgrpDiv->maxLength())
|
467 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - mean length is greater than maximum length for stock", this->getName());
|
468 : |
|
|
|
469 : |
|
|
sum = 0.0;
|
470 : |
|
|
mult = 1.0 / sdevLength[i];
|
471 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++) {
|
472 : |
|
|
dnorm = (LgrpDiv->meanLength(l) - meanLength[i]) * mult;
|
473 : |
|
|
renewalDistribution[i][age][l].N = exp(-(dnorm * dnorm) * 0.5);
|
474 : |
|
|
sum += renewalDistribution[i][age][l].N;
|
475 : |
|
|
}
|
476 : |
|
|
|
477 : |
|
|
if (isZero(sum)) {
|
478 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - calculated zero recruits for stock", this->getName());
|
479 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++)
|
480 : |
|
|
renewalDistribution[i][age][l].setToZero();
|
481 : |
|
|
|
482 : |
|
|
} else {
|
483 : |
|
|
sum = 10000.0 / sum;
|
484 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++) {
|
485 : |
|
|
renewalDistribution[i][age][l].N *= sum;
|
486 : |
|
|
renewalDistribution[i][age][l].W = alpha[i] * pow(LgrpDiv->meanLength(l), beta[i]);
|
487 : |
|
|
if ((handle.getLogLevel() >= LOGWARN) && (isZero(renewalDistribution[i][age][l].W)) && (renewalDistribution[i][age][l].N > 0.0))
|
488 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - zero mean weight for stock", this->getName());
|
489 : |
|
|
}
|
490 : |
|
|
}
|
491 : |
|
|
}
|
492 : |
|
|
}
|
493 : |
|
|
|
494 : |
|
|
} else if (readoption == 2) {
|
495 : |
|
|
for (i = 0; i < renewalTime.Size(); i++) {
|
496 : |
|
|
age = renewalAge[i];
|
497 : |
|
|
minage = renewalDistribution[i].minAge();
|
498 : |
|
|
for (l = renewalDistribution[i].minLength(age); l < renewalDistribution[i].maxLength(age); l++) {
|
499 : |
|
|
renewalDistribution[i][age][l].N = (*renewalNumber[i])[age - minage][l];
|
500 : |
|
|
if (handle.getLogLevel() >= LOGWARN) {
|
501 : |
|
|
if (renewalDistribution[i][age][l].N < 0.0)
|
502 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - negative number of recruits", renewalDistribution[i][age][l].N);
|
503 : |
|
|
if ((isZero(renewalDistribution[i][age][l].W)) && (renewalDistribution[i][age][l].N > 0.0))
|
504 : |
|
|
handle.logMessage(LOGWARN, "Warning in renewal - zero mean weight for stock", this->getName());
|
505 : |
|
|
}
|
506 : |
|
|
}
|
507 : |
|
|
}
|
508 : |
|
|
|
509 : |
|
|
} else
|
510 : |
|
|
handle.logMessage(LOGFAIL, "Error in renewal - unrecognised data format");
|
511 : |
|
|
|
512 : |
|
|
if (handle.getLogLevel() >= LOGMESSAGE)
|
513 : |
|
|
handle.logMessage(LOGMESSAGE, "Reset renewal data for stock", this->getName());
|
514 : |
|
|
}
|
515 : |
|
|
|
516 : |
|
|
int RenewalData::isRenewalStepArea(int area, const TimeClass* const TimeInfo) {
|
517 : |
|
|
int i;
|
518 : |
|
|
for (i = 0; i < renewalTime.Size(); i++)
|
519 : |
|
|
if ((renewalTime[i] == TimeInfo->getTime()) && (renewalArea[i] == area))
|
520 : |
|
|
return 1;
|
521 : |
|
|
return 0;
|
522 : |
|
|
}
|
523 : |
|
|
|
524 : |
|
|
void RenewalData::addRenewal(AgeBandMatrix& Alkeys, int area, const TimeClass* const TimeInfo) {
|
525 : |
|
|
int i;
|
526 : |
|
|
for (i = 0; i < renewalTime.Size(); i++) {
|
527 : |
|
|
if ((renewalTime[i] == TimeInfo->getTime()) && (renewalArea[i] == area)) {
|
528 : |
|
|
index = i;
|
529 : |
|
|
if (readoption == 2)
|
530 : |
|
|
Alkeys.Add(renewalDistribution[i], *CI);
|
531 : |
|
|
else if (renewalMult[i] > verysmall)
|
532 : |
|
|
Alkeys.Add(renewalDistribution[i], *CI, renewalMult[i]);
|
533 : |
|
|
}
|
534 : |
|
|
}
|
535 : |
|
|
}
|