1 : |
agomez |
1 |
#include "migration.h" |
2 : |
|
|
#include "readfunc.h" |
3 : |
|
|
#include "readword.h" |
4 : |
|
|
#include "errorhandler.h" |
5 : |
|
|
#include "mathfunc.h" |
6 : |
|
|
#include "migrationarea.h" |
7 : |
|
|
#include "gadget.h" |
8 : |
|
|
#include "global.h" |
9 : |
|
|
|
10 : |
|
|
// ********************************************************
|
11 : |
|
|
// Functions for base Migration
|
12 : |
|
|
// ********************************************************
|
13 : |
|
|
Migration::Migration(const IntVector& Areas, const char* givenname)
|
14 : |
|
|
: HasName(givenname), LivesOnAreas(Areas) {
|
15 : |
|
|
|
16 : |
|
|
if (areas.Size() == 1)
|
17 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - only one area defined");
|
18 : |
|
|
}
|
19 : |
|
|
|
20 : |
|
|
// ********************************************************
|
21 : |
|
|
// Functions for MigrationNumbers
|
22 : |
|
|
// ********************************************************
|
23 : |
|
|
MigrationNumbers::MigrationNumbers(CommentStream& infile, const IntVector& Areas,
|
24 : |
|
|
const AreaClass* const Area, const TimeClass* const TimeInfo,
|
25 : |
|
|
const char* givenname, Keeper* const keeper)
|
26 : |
|
|
: Migration(Areas, givenname) {
|
27 : |
|
|
|
28 : |
|
|
int i, j;
|
29 : |
|
|
ifstream subfile;
|
30 : |
|
|
CommentStream subcomment(subfile);
|
31 : |
|
|
char text[MaxStrLength];
|
32 : |
|
|
char filename[MaxStrLength];
|
33 : |
|
|
strncpy(text, "", MaxStrLength);
|
34 : |
|
|
strncpy(filename, "", MaxStrLength);
|
35 : |
|
|
keeper->addString("migrationmatrix");
|
36 : |
|
|
|
37 : |
|
|
// Open and read year and step information
|
38 : |
|
|
readWordAndValue(infile, "yearstepfile", filename);
|
39 : |
|
|
subfile.open(filename, ios::in);
|
40 : |
|
|
handle.checkIfFailure(subfile, filename);
|
41 : |
|
|
handle.Open(filename);
|
42 : |
|
|
this->readTimeStepData(subcomment, TimeInfo);
|
43 : |
|
|
handle.Close();
|
44 : |
|
|
subfile.close();
|
45 : |
|
|
subfile.clear();
|
46 : |
|
|
|
47 : |
|
|
// Open and read information about how the matrices are defined
|
48 : |
|
|
infile >> text >> filename >> ws;
|
49 : |
|
|
subfile.open(filename, ios::in);
|
50 : |
|
|
handle.checkIfFailure(subfile, filename);
|
51 : |
|
|
handle.Open(filename);
|
52 : |
|
|
|
53 : |
|
|
if (strcasecmp(text, "definematrices") == 0) {
|
54 : |
|
|
this->readGivenMatrices(subcomment, keeper);
|
55 : |
|
|
} else if (strcasecmp(text, "defineratios") == 0) {
|
56 : |
|
|
this->readGivenRatios(subcomment, keeper, Area);
|
57 : |
|
|
} else
|
58 : |
|
|
handle.logFileUnexpected(LOGFAIL, "definematrices or defineratios", text);
|
59 : |
|
|
|
60 : |
|
|
handle.Close();
|
61 : |
|
|
subfile.close();
|
62 : |
|
|
subfile.clear();
|
63 : |
|
|
|
64 : |
|
|
this->checkMatrixIndex();
|
65 : |
|
|
handle.logMessage(LOGMESSAGE, "Read migration numbers file - number of migration matrices", readMigration.Size());
|
66 : |
|
|
keeper->clearLast();
|
67 : |
|
|
}
|
68 : |
|
|
|
69 : |
|
|
MigrationNumbers::~MigrationNumbers() {
|
70 : |
|
|
int i;
|
71 : |
|
|
for (i = 0; i < timeindex.Size(); i++)
|
72 : |
|
|
if (timeindex[i] != -1)
|
73 : |
|
|
delete[] allmatrixnames[i];
|
74 : |
|
|
for (i = 0; i < usedmatrixnames.Size(); i++) {
|
75 : |
|
|
delete[] usedmatrixnames[i];
|
76 : |
|
|
delete readMigration[i];
|
77 : |
|
|
delete calcMigration[i];
|
78 : |
|
|
}
|
79 : |
|
|
}
|
80 : |
|
|
|
81 : |
|
|
void MigrationNumbers::readTimeStepData(CommentStream& infile, const TimeClass* const TimeInfo) {
|
82 : |
|
|
|
83 : |
|
|
int year, step, timeid, count;
|
84 : |
|
|
char text[MaxStrLength];
|
85 : |
|
|
strncpy(text, "", MaxStrLength);
|
86 : |
|
|
|
87 : |
|
|
// check the number of columns in the inputfile
|
88 : |
|
|
infile >> ws;
|
89 : |
|
|
if (countColumns(infile) != 3)
|
90 : |
|
|
handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 3");
|
91 : |
|
|
|
92 : |
|
|
year = step = count = 0;
|
93 : |
|
|
allmatrixnames.resizeBlank(TimeInfo->numTotalSteps() + 1);
|
94 : |
|
|
timeindex.resize(TimeInfo->numTotalSteps() + 1, -1);
|
95 : |
|
|
while (!infile.eof()) {
|
96 : |
|
|
infile >> year >> step >> text >> ws;
|
97 : |
|
|
if (TimeInfo->isWithinPeriod(year, step)) {
|
98 : |
|
|
timeid = TimeInfo->calcSteps(year, step);
|
99 : |
|
|
timeindex[timeid] = 0;
|
100 : |
|
|
allmatrixnames[timeid] = new char[strlen(text) + 1];
|
101 : |
|
|
strcpy(allmatrixnames[timeid], text);
|
102 : |
|
|
count++;
|
103 : |
|
|
}
|
104 : |
|
|
}
|
105 : |
|
|
|
106 : |
|
|
if (count == 0)
|
107 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - found no migration data");
|
108 : |
|
|
if (count != TimeInfo->numTotalSteps())
|
109 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - migration data doesnt span time range");
|
110 : |
|
|
handle.logMessage(LOGMESSAGE, "Read migration time data - number of entries", count);
|
111 : |
|
|
}
|
112 : |
|
|
|
113 : |
|
|
void MigrationNumbers::readGivenRatios(CommentStream& infile, Keeper* const keeper, const AreaClass* const Area) {
|
114 : |
|
|
|
115 : |
|
|
int innerrow, innercol, fromArea, toArea;
|
116 : |
|
|
int i, j, k, count;
|
117 : |
|
|
int numAreas = areas.Size();
|
118 : |
|
|
char text[MaxStrLength];
|
119 : |
|
|
char name[MaxStrLength];
|
120 : |
|
|
strncpy(text, "", MaxStrLength);
|
121 : |
|
|
strncpy(name, "", MaxStrLength);
|
122 : |
|
|
|
123 : |
|
|
/* File format:
|
124 : |
|
|
[migrationmatrix]
|
125 : |
|
|
name matrixname
|
126 : |
|
|
fromarea toarea ratio
|
127 : |
|
|
[migrationmatrix]
|
128 : |
|
|
...
|
129 : |
|
|
*/
|
130 : |
|
|
|
131 : |
|
|
toArea = fromArea = 0;
|
132 : |
|
|
// Read ratio information from file
|
133 : |
|
|
while (!infile.eof() && !infile.fail()) {
|
134 : |
|
|
infile >> text >> ws;
|
135 : |
|
|
if (strcasecmp(text, "[migrationmatrix]") != 0)
|
136 : |
|
|
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix]", text);
|
137 : |
|
|
|
138 : |
|
|
readWordAndValue(infile, "name", name);
|
139 : |
|
|
if (this->useMatrix(name)) {
|
140 : |
|
|
this->setMatrixName(name);
|
141 : |
|
|
count = readMigration.Size();
|
142 : |
|
|
checkvalues.AddRows(1, numAreas, 0);
|
143 : |
|
|
readMigration.resize(new FormulaMatrix(numAreas, numAreas, 0.0));
|
144 : |
|
|
|
145 : |
|
|
infile >> ws;
|
146 : |
|
|
while (!infile.eof() && !infile.fail() && (infile.peek() != '[')) {
|
147 : |
|
|
if (!(isdigit(infile.peek()))) {
|
148 : |
|
|
infile >> text >> ws;
|
149 : |
|
|
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix] or area number", text);
|
150 : |
|
|
}
|
151 : |
|
|
|
152 : |
|
|
infile >> fromArea >> toArea >> ws;
|
153 : |
|
|
innerrow = Area->getInnerArea(toArea);
|
154 : |
|
|
innercol = Area->getInnerArea(fromArea);
|
155 : |
|
|
if (!this->isInArea(innerrow) || !this->isInArea(innercol))
|
156 : |
|
|
handle.logMessage(LOGFAIL, "Error in migration - invalid area for matrix", name);
|
157 : |
|
|
|
158 : |
|
|
toArea = this->areaNum(innerrow);
|
159 : |
|
|
fromArea = this->areaNum(innercol);
|
160 : |
|
|
infile >> (*readMigration[count])[toArea][fromArea] >> ws;
|
161 : |
|
|
checkvalues[count][fromArea]++;
|
162 : |
|
|
}
|
163 : |
|
|
|
164 : |
|
|
} else {
|
165 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - matrix not used in the simulation", name);
|
166 : |
|
|
infile >> ws;
|
167 : |
|
|
while (!infile.eof() && !infile.fail() && infile.peek() != '[')
|
168 : |
|
|
infile >> text >> ws;
|
169 : |
|
|
}
|
170 : |
|
|
}
|
171 : |
|
|
|
172 : |
|
|
// Finished reading from file, now check the values
|
173 : |
|
|
Formula missing;
|
174 : |
|
|
missing.setValue(1.0);
|
175 : |
|
|
for (i = 0; i < readMigration.Size(); i++) {
|
176 : |
|
|
for (j = 0; j < readMigration[i]->Nrow(); j++) {
|
177 : |
|
|
if (checkvalues[i][j] == 0) {
|
178 : |
|
|
// no data read for migration from this area
|
179 : |
|
|
(*readMigration[i])[j][j].setValue(1.0);
|
180 : |
|
|
|
181 : |
|
|
} else if (checkvalues[i][j] == numAreas) {
|
182 : |
|
|
// nothing to do - data read for migration from this area to every other area
|
183 : |
|
|
|
184 : |
|
|
} else if (isZero((*readMigration[i])[j][j])) {
|
185 : |
|
|
// read data for migration from this area to other areas
|
186 : |
|
|
// so simply subtract the migration from 1 to find those that are left
|
187 : |
|
|
count = 0;
|
188 : |
|
|
vector<Formula*> checklist;
|
189 : |
|
|
for (k = 0; k < readMigration[i]->Ncol(); k++) {
|
190 : |
|
|
if (!(isZero((*readMigration[i])[k][j]))) {
|
191 : |
|
|
checklist.push_back(&(*readMigration[i])[k][j]);
|
192 : |
|
|
count++;
|
193 : |
|
|
}
|
194 : |
|
|
}
|
195 : |
|
|
|
196 : |
|
|
if (count == 0)
|
197 : |
|
|
handle.logFileMessage(LOGFAIL, "missing data for migration matrix", usedmatrixnames[i]);
|
198 : |
|
|
checklist.insert(checklist.begin(), &missing);
|
199 : |
|
|
Formula* missingValue = new Formula(MINUS, checklist);
|
200 : |
|
|
(*readMigration[i])[j][j] = *missingValue;
|
201 : |
|
|
delete missingValue;
|
202 : |
|
|
|
203 : |
|
|
} else if (checkvalues[i][j] == (numAreas - 1)) {
|
204 : |
|
|
// read migration data from this area to all other areas but 1
|
205 : |
|
|
count = -1;
|
206 : |
|
|
vector<Formula*> checklist;
|
207 : |
|
|
for (k = 0; k < readMigration[i]->Ncol(); k++) {
|
208 : |
|
|
if (isZero((*readMigration[i])[k][j]))
|
209 : |
|
|
count = k;
|
210 : |
|
|
else
|
211 : |
|
|
checklist.push_back(&(*readMigration[i])[k][j]);
|
212 : |
|
|
}
|
213 : |
|
|
|
214 : |
|
|
if (count == -1)
|
215 : |
|
|
handle.logFileMessage(LOGFAIL, "missing data for migration matrix", usedmatrixnames[i]);
|
216 : |
|
|
checklist.insert(checklist.begin(), &missing);
|
217 : |
|
|
Formula* missingValue = new Formula(MINUS, checklist);
|
218 : |
|
|
(*readMigration[i])[count][j] = *missingValue;
|
219 : |
|
|
delete missingValue;
|
220 : |
|
|
|
221 : |
|
|
} else {
|
222 : |
|
|
// havent read enough information to reliably construct matrix
|
223 : |
|
|
handle.logFileMessage(LOGFAIL, "missing data for migration matrix", usedmatrixnames[i]);
|
224 : |
|
|
|
225 : |
|
|
}
|
226 : |
|
|
}
|
227 : |
|
|
}
|
228 : |
|
|
|
229 : |
|
|
// Inform keeper of the values and resize
|
230 : |
|
|
for (i = 0; i < readMigration.Size(); i++) {
|
231 : |
|
|
(*readMigration[i]).Inform(keeper);
|
232 : |
|
|
calcMigration.resize(new DoubleMatrix(numAreas, numAreas, 0.0));
|
233 : |
|
|
}
|
234 : |
|
|
}
|
235 : |
|
|
|
236 : |
|
|
void MigrationNumbers::readGivenMatrices(CommentStream& infile, Keeper* const keeper) {
|
237 : |
|
|
|
238 : |
|
|
int i, j, count;
|
239 : |
|
|
int numAreas = areas.Size();
|
240 : |
|
|
char text[MaxStrLength];
|
241 : |
|
|
char name[MaxStrLength];
|
242 : |
|
|
strncpy(text, "", MaxStrLength);
|
243 : |
|
|
strncpy(name, "", MaxStrLength);
|
244 : |
|
|
|
245 : |
|
|
/* File format:
|
246 : |
|
|
[migrationmatrix]
|
247 : |
|
|
name matrixname
|
248 : |
|
|
x11 ... x1n
|
249 : |
|
|
...
|
250 : |
|
|
xn1 ... xnn
|
251 : |
|
|
[migrationmatrix]
|
252 : |
|
|
...
|
253 : |
|
|
*/
|
254 : |
|
|
|
255 : |
|
|
// Read matrix information from file
|
256 : |
|
|
while (!infile.eof() && !infile.fail()) {
|
257 : |
|
|
infile >> text >> ws;
|
258 : |
|
|
if (strcasecmp(text, "[migrationmatrix]") != 0)
|
259 : |
|
|
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix]", text);
|
260 : |
|
|
|
261 : |
|
|
readWordAndValue(infile, "name", name);
|
262 : |
|
|
if (this->useMatrix(name)) {
|
263 : |
|
|
this->setMatrixName(name);
|
264 : |
|
|
count = readMigration.Size();
|
265 : |
|
|
readMigration.resize(new FormulaMatrix(numAreas, numAreas, 0.0));
|
266 : |
|
|
|
267 : |
|
|
infile >> ws;
|
268 : |
|
|
if (infile.peek() == '[' || infile.eof()) {
|
269 : |
|
|
for (i = 0; i < numAreas; i++)
|
270 : |
|
|
(*readMigration[count])[i][i].setValue(1.0);
|
271 : |
|
|
|
272 : |
|
|
} else if (isdigit(infile.peek()) || infile.peek() == '#' || infile.peek() == '(') {
|
273 : |
|
|
// look for a number or a formula as the next entry to be read
|
274 : |
|
|
for (i = 0; i < numAreas; i++)
|
275 : |
|
|
for (j = 0; j < numAreas; j++)
|
276 : |
|
|
if (!(infile >> (*readMigration[count])[i][j]))
|
277 : |
|
|
handle.logMessage(LOGFAIL, "Error in migration - failed to read matrix", name);
|
278 : |
|
|
infile >> ws; //JMB strip whitespace
|
279 : |
|
|
|
280 : |
|
|
} else {
|
281 : |
|
|
infile >> text;
|
282 : |
|
|
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix] or matrix value", text);
|
283 : |
|
|
}
|
284 : |
|
|
|
285 : |
|
|
} else {
|
286 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - matrix not used in the simulation", name);
|
287 : |
|
|
infile >> ws;
|
288 : |
|
|
while (!infile.eof() && !infile.fail() && infile.peek() != '[')
|
289 : |
|
|
infile >> text >> ws;
|
290 : |
|
|
}
|
291 : |
|
|
}
|
292 : |
|
|
|
293 : |
|
|
// Inform keeper of the values and resize
|
294 : |
|
|
for (i = 0; i < readMigration.Size(); i++) {
|
295 : |
|
|
(*readMigration[i]).Inform(keeper);
|
296 : |
|
|
calcMigration.resize(new DoubleMatrix(numAreas, numAreas, 0.0));
|
297 : |
|
|
}
|
298 : |
|
|
}
|
299 : |
|
|
|
300 : |
|
|
void MigrationNumbers::checkMatrixIndex() {
|
301 : |
|
|
int i, j, check;
|
302 : |
|
|
|
303 : |
|
|
for (i = 0; i < timeindex.Size(); i++) {
|
304 : |
|
|
if (timeindex[i] == 0) {
|
305 : |
|
|
check = -1;
|
306 : |
|
|
for (j = 0; j < usedmatrixnames.Size(); j++)
|
307 : |
|
|
if (strcasecmp(allmatrixnames[i], usedmatrixnames[j]) == 0)
|
308 : |
|
|
check = j;
|
309 : |
|
|
|
310 : |
|
|
if (check == -1)
|
311 : |
|
|
handle.logMessage(LOGFAIL, "Error in migration - failed to read matrix", allmatrixnames[i]);
|
312 : |
|
|
timeindex[i] = check;
|
313 : |
|
|
|
314 : |
|
|
} else if (timeindex[i] != -1) {
|
315 : |
|
|
handle.logMessage(LOGFAIL, "Error in migration - repeated migration matrix", i);
|
316 : |
|
|
}
|
317 : |
|
|
}
|
318 : |
|
|
}
|
319 : |
|
|
|
320 : |
|
|
void MigrationNumbers::Print(ofstream& outfile) {
|
321 : |
|
|
int i, j, k;
|
322 : |
|
|
outfile << "\nMigration\n\tNames of migration matrices:\n\t";
|
323 : |
|
|
for (i = 0; i < timeindex.Size(); i++)
|
324 : |
|
|
if (timeindex[i] != -1)
|
325 : |
|
|
outfile << allmatrixnames[i] << sep;
|
326 : |
|
|
outfile << "\n\n\tMigration matrices";
|
327 : |
|
|
for (i = 0; i < calcMigration.Size(); i++) {
|
328 : |
|
|
outfile << "\n\tMatrix name: " << usedmatrixnames[i] << "\n\t";
|
329 : |
|
|
for (j = 0; j < calcMigration[i]->Nrow(); j++) {
|
330 : |
|
|
for (k = 0; k < calcMigration[i]->Ncol(j); k++)
|
331 : |
|
|
outfile << setw(smallwidth) << (*calcMigration[i])[j][k] << sep;
|
332 : |
|
|
outfile << "\n\t";
|
333 : |
|
|
}
|
334 : |
|
|
}
|
335 : |
|
|
outfile.flush();
|
336 : |
|
|
}
|
337 : |
|
|
|
338 : |
|
|
const DoubleMatrix& MigrationNumbers::getMigrationMatrix(const TimeClass* const TimeInfo) {
|
339 : |
|
|
return (*calcMigration[timeindex[TimeInfo->getTime()]]);
|
340 : |
|
|
}
|
341 : |
|
|
|
342 : |
|
|
void MigrationNumbers::setMatrixName(char* name) {
|
343 : |
|
|
int i;
|
344 : |
|
|
// check to ensure that this matrix name is unique
|
345 : |
|
|
for (i = 0; i < usedmatrixnames.Size(); i++)
|
346 : |
|
|
if (strcasecmp(usedmatrixnames[i], name) == 0)
|
347 : |
|
|
handle.logMessage(LOGFAIL, "Error in migration - repeated matrix", name);
|
348 : |
|
|
|
349 : |
|
|
// store the new matrix name
|
350 : |
|
|
char* tempname = new char[strlen(name) + 1];
|
351 : |
|
|
strcpy(tempname, name);
|
352 : |
|
|
usedmatrixnames.resize(tempname);
|
353 : |
|
|
}
|
354 : |
|
|
|
355 : |
|
|
void MigrationNumbers::Reset() {
|
356 : |
|
|
//JMB need to reset the penalty vector first
|
357 : |
|
|
penalty.Reset();
|
358 : |
|
|
|
359 : |
|
|
int i, j, k;
|
360 : |
|
|
double colsum;
|
361 : |
|
|
for (i = 0; i < readMigration.Size(); i++) {
|
362 : |
|
|
for (k = 0; k < readMigration[i]->Ncol(); k++) {
|
363 : |
|
|
colsum = 0.0;
|
364 : |
|
|
for (j = 0; j < readMigration[i]->Nrow(); j++) {
|
365 : |
|
|
if (isZero((*readMigration[i])[j][k])) {
|
366 : |
|
|
(*calcMigration[i])[j][k] = 0.0;
|
367 : |
|
|
} else if ((*readMigration[i])[j][k] < 0.0) {
|
368 : |
|
|
penalty.resize(1, (*readMigration[i])[j][k]);
|
369 : |
|
|
(*calcMigration[i])[j][k] = 0.0;
|
370 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - value outside bounds", (*readMigration[i])[j][k]);
|
371 : |
|
|
} else if ((*readMigration[i])[j][k] > 1.0) {
|
372 : |
|
|
penalty.resize(1, (*readMigration[i])[j][k]);
|
373 : |
|
|
(*calcMigration[i])[j][k] = 1.0;
|
374 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - value outside bounds", (*readMigration[i])[j][k]);
|
375 : |
|
|
} else
|
376 : |
|
|
(*calcMigration[i])[j][k] = (*readMigration[i])[j][k];
|
377 : |
|
|
|
378 : |
|
|
colsum += (*calcMigration[i])[j][k];
|
379 : |
|
|
}
|
380 : |
|
|
|
381 : |
|
|
if (isZero(colsum)) {
|
382 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - column sum is zero");
|
383 : |
|
|
// just putting anything here ... AJ 23.06.05, Needs to be changed!!!!!
|
384 : |
|
|
penalty.resize(1, 1.0);
|
385 : |
|
|
for (j = 0; j < readMigration[i]->Nrow(); j++) {
|
386 : |
|
|
if (k == j)
|
387 : |
|
|
(*calcMigration[i])[j][k] = 1.0;
|
388 : |
|
|
else
|
389 : |
|
|
(*calcMigration[i])[j][k] = 0.0;
|
390 : |
|
|
}
|
391 : |
|
|
|
392 : |
|
|
} else if (isEqual(colsum, 1.0)) {
|
393 : |
|
|
// everything is OK ...
|
394 : |
|
|
|
395 : |
|
|
} else {
|
396 : |
|
|
penalty.resize(1, colsum);
|
397 : |
|
|
colsum = 1.0 / colsum;
|
398 : |
|
|
for (j = 0; j < readMigration[i]->Nrow(); j++)
|
399 : |
|
|
(*calcMigration[i])[j][k] *= colsum;
|
400 : |
|
|
}
|
401 : |
|
|
}
|
402 : |
|
|
}
|
403 : |
|
|
|
404 : |
|
|
if (handle.getLogLevel() >= LOGMESSAGE)
|
405 : |
|
|
handle.logMessage(LOGMESSAGE, "Reset migration data for stock", this->getName());
|
406 : |
|
|
}
|
407 : |
|
|
|
408 : |
|
|
int MigrationNumbers::useMatrix(char* name) {
|
409 : |
|
|
int i;
|
410 : |
|
|
for (i = 0; i < timeindex.Size(); i++)
|
411 : |
|
|
if ((timeindex[i] != -1) && (strcasecmp(allmatrixnames[i], name) == 0))
|
412 : |
|
|
return 1;
|
413 : |
|
|
return 0;
|
414 : |
|
|
}
|
415 : |
|
|
|
416 : |
|
|
int MigrationNumbers::isMigrationStep(const TimeClass* const TimeInfo) {
|
417 : |
|
|
if (timeindex[TimeInfo->getTime()] == -1)
|
418 : |
|
|
return 0;
|
419 : |
|
|
return 1;
|
420 : |
|
|
}
|
421 : |
|
|
|
422 : |
|
|
// ********************************************************
|
423 : |
|
|
// Functions for MigrationFunction
|
424 : |
|
|
// ********************************************************
|
425 : |
|
|
MigrationFunction::MigrationFunction(CommentStream& infile, const IntVector& Areas,
|
426 : |
|
|
const AreaClass* const Area, const TimeClass* const TimeInfo,
|
427 : |
|
|
const char* givenname, Keeper* const keeper)
|
428 : |
|
|
: Migration(Areas, givenname) {
|
429 : |
|
|
|
430 : |
|
|
int i;
|
431 : |
|
|
ifstream subfile;
|
432 : |
|
|
CommentStream subcomment(subfile);
|
433 : |
|
|
char text[MaxStrLength];
|
434 : |
|
|
char filename[MaxStrLength];
|
435 : |
|
|
strncpy(text, "", MaxStrLength);
|
436 : |
|
|
strncpy(filename, "", MaxStrLength);
|
437 : |
|
|
keeper->addString("migrationfunction");
|
438 : |
|
|
|
439 : |
|
|
readWordAndModelVariable(infile, "diffusion", diffusion, TimeInfo, keeper);
|
440 : |
|
|
readWordAndModelVariable(infile, "driftx", driftx, TimeInfo, keeper);
|
441 : |
|
|
readWordAndModelVariable(infile, "drifty", drifty, TimeInfo, keeper);
|
442 : |
|
|
readWordAndVariable(infile, "lambda", lambda);
|
443 : |
|
|
delta = TimeInfo->getTimeStepLength() / TimeInfo->numSubSteps();
|
444 : |
|
|
|
445 : |
|
|
// Open and read given year and step information
|
446 : |
|
|
readWordAndValue(infile, "areadefinition", filename);
|
447 : |
|
|
subfile.open(filename, ios::in);
|
448 : |
|
|
handle.checkIfFailure(subfile, filename);
|
449 : |
|
|
handle.Open(filename);
|
450 : |
|
|
this->readAreaData(subcomment, Area, TimeInfo, keeper);
|
451 : |
|
|
handle.Close();
|
452 : |
|
|
subfile.close();
|
453 : |
|
|
subfile.clear();
|
454 : |
|
|
|
455 : |
|
|
calcMigration.AddRows(oceanareas.Size(), oceanareas.Size(), 0.0);
|
456 : |
|
|
if (oceanareas.Size() != areas.Size())
|
457 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - migration data doesnt cover all areas");
|
458 : |
|
|
|
459 : |
|
|
handle.logMessage(LOGMESSAGE, "Read migration function file - number of migration areas", oceanareas.Size());
|
460 : |
|
|
keeper->clearLast();
|
461 : |
|
|
}
|
462 : |
|
|
|
463 : |
|
|
void MigrationFunction::readAreaData(CommentStream& infile, const AreaClass* const Area,
|
464 : |
|
|
const TimeClass* const TimeInfo, Keeper* const keeper) {
|
465 : |
|
|
|
466 : |
|
|
int area, inarea;
|
467 : |
|
|
ifstream subfile;
|
468 : |
|
|
CommentStream subcomment(subfile);
|
469 : |
|
|
char text[MaxStrLength];
|
470 : |
|
|
char filename[MaxStrLength];
|
471 : |
|
|
char areaname[MaxStrLength];
|
472 : |
|
|
strncpy(text, "", MaxStrLength);
|
473 : |
|
|
strncpy(filename, "", MaxStrLength);
|
474 : |
|
|
strncpy(areaname, "", MaxStrLength);
|
475 : |
|
|
|
476 : |
|
|
/* File format:
|
477 : |
|
|
[area]
|
478 : |
|
|
name areaname
|
479 : |
|
|
number areanumber
|
480 : |
|
|
rectangles filename containing the rectangle information
|
481 : |
|
|
[area]
|
482 : |
|
|
...
|
483 : |
|
|
*/
|
484 : |
|
|
|
485 : |
|
|
// Read area information from file
|
486 : |
|
|
infile >> ws;
|
487 : |
|
|
while (!infile.eof() && !infile.fail()) {
|
488 : |
|
|
infile >> text >> ws;
|
489 : |
|
|
if (strcasecmp(text, "[area]") != 0)
|
490 : |
|
|
handle.logFileUnexpected(LOGFAIL, "[area]", text);
|
491 : |
|
|
|
492 : |
|
|
readWordAndValue(infile, "name", areaname);
|
493 : |
|
|
readWordAndVariable(infile, "number", area);
|
494 : |
|
|
readWordAndValue(infile, "rectangles", filename);
|
495 : |
|
|
|
496 : |
|
|
inarea = Area->getInnerArea(area);
|
497 : |
|
|
if (!this->isInArea(inarea))
|
498 : |
|
|
handle.logMessage(LOGFAIL, "Error in migration - invalid area", area);
|
499 : |
|
|
|
500 : |
|
|
subfile.open(filename, ios::in);
|
501 : |
|
|
handle.checkIfFailure(subfile, filename);
|
502 : |
|
|
handle.Open(filename);
|
503 : |
|
|
oceanareas.resize(new MigrationArea(subcomment, areaname, inarea));
|
504 : |
|
|
handle.Close();
|
505 : |
|
|
subfile.close();
|
506 : |
|
|
subfile.clear();
|
507 : |
|
|
}
|
508 : |
|
|
}
|
509 : |
|
|
|
510 : |
|
|
MigrationFunction::~MigrationFunction() {
|
511 : |
|
|
int i;
|
512 : |
|
|
for (i = 0; i < oceanareas.Size(); i++)
|
513 : |
|
|
delete oceanareas[i];
|
514 : |
|
|
}
|
515 : |
|
|
|
516 : |
|
|
void MigrationFunction::Print(ofstream& outfile) {
|
517 : |
|
|
int i, j;
|
518 : |
|
|
outfile << "\nMigration\n\n\tMigration matrix\n\t";
|
519 : |
|
|
for (i = 0; i < calcMigration.Nrow(); i++) {
|
520 : |
|
|
for (j = 0; j < calcMigration.Ncol(i); j++)
|
521 : |
|
|
outfile << setw(smallwidth) << calcMigration[i][j] << sep;
|
522 : |
|
|
outfile << "\n\t";
|
523 : |
|
|
}
|
524 : |
|
|
outfile.flush();
|
525 : |
|
|
}
|
526 : |
|
|
|
527 : |
|
|
const DoubleMatrix& MigrationFunction::getMigrationMatrix(const TimeClass* const TimeInfo) {
|
528 : |
|
|
if (this->updateVariables(TimeInfo))
|
529 : |
|
|
this->recalcMatrix();
|
530 : |
|
|
return calcMigration;
|
531 : |
|
|
}
|
532 : |
|
|
|
533 : |
|
|
void MigrationFunction::recalcMatrix() {
|
534 : |
|
|
int from, to, idfrom, idto, i, j;
|
535 : |
|
|
double mig, sa;
|
536 : |
|
|
double colsum, sum;
|
537 : |
|
|
|
538 : |
|
|
for (from = 0; from < oceanareas.Size(); from++) {
|
539 : |
|
|
idfrom = this->areaNum(oceanareas[from]->getAreaID());
|
540 : |
|
|
colsum = 0.0;
|
541 : |
|
|
for (to = 0; to < oceanareas.Size(); to++) {
|
542 : |
|
|
idto = this->areaNum(oceanareas[to]->getAreaID());
|
543 : |
|
|
sum = 0.0;
|
544 : |
|
|
for (i = 0; i < oceanareas[from]->getNumRectangles(); i++) {
|
545 : |
|
|
sa = (oceanareas[from]->getRectangles()[i])->getArea();
|
546 : |
|
|
if (!(isZero(sa))) {
|
547 : |
|
|
for (j = 0; j < oceanareas[to]->getNumRectangles(); j++) {
|
548 : |
|
|
mig = this->getMigrationFunction(oceanareas[from]->getRectangles()[i], oceanareas[to]->getRectangles()[j]);
|
549 : |
|
|
sum += (mig * sa);
|
550 : |
|
|
}
|
551 : |
|
|
}
|
552 : |
|
|
}
|
553 : |
|
|
|
554 : |
|
|
calcMigration[idto][idfrom] = sum / oceanareas[from]->getArea();
|
555 : |
|
|
colsum += calcMigration[idto][idfrom];
|
556 : |
|
|
}
|
557 : |
|
|
|
558 : |
|
|
if (isZero(colsum)) {
|
559 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - column sum is zero");
|
560 : |
|
|
penalty.resize(1, 1.0);
|
561 : |
|
|
calcMigration[idfrom][idfrom] = 1.0;
|
562 : |
|
|
|
563 : |
|
|
} else if (isEqual(colsum, 1.0)) {
|
564 : |
|
|
// everything is OK ...
|
565 : |
|
|
|
566 : |
|
|
} else {
|
567 : |
|
|
handle.logMessage(LOGWARN, "Warning in migration - column sum", colsum);
|
568 : |
|
|
penalty.resize(1, colsum);
|
569 : |
|
|
colsum = 1.0 / colsum;
|
570 : |
|
|
for (i = 0; i < calcMigration.Ncol(); i++)
|
571 : |
|
|
calcMigration[i][idfrom] *= colsum;
|
572 : |
|
|
}
|
573 : |
|
|
}
|
574 : |
|
|
}
|
575 : |
|
|
|
576 : |
|
|
void MigrationFunction::Reset() {
|
577 : |
|
|
//JMB need to reset the penalty vector
|
578 : |
|
|
penalty.Reset();
|
579 : |
|
|
if (handle.getLogLevel() >= LOGMESSAGE)
|
580 : |
|
|
handle.logMessage(LOGMESSAGE, "Reset migration data for stock", this->getName());
|
581 : |
|
|
}
|
582 : |
|
|
|
583 : |
|
|
int MigrationFunction::isMigrationStep(const TimeClass* const TimeInfo) {
|
584 : |
|
|
return 1;
|
585 : |
|
|
}
|
586 : |
|
|
|
587 : |
|
|
int MigrationFunction::updateVariables(const TimeClass* const TimeInfo) {
|
588 : |
|
|
//update the values of the variables that can change
|
589 : |
|
|
delta = TimeInfo->getTimeStepLength() / TimeInfo->numSubSteps();
|
590 : |
|
|
diffusion.Update(TimeInfo);
|
591 : |
|
|
driftx.Update(TimeInfo);
|
592 : |
|
|
drifty.Update(TimeInfo);
|
593 : |
|
|
|
594 : |
|
|
if ((TimeInfo->didStepSizeChange()) || (diffusion.didChange(TimeInfo))
|
595 : |
|
|
|| (driftx.didChange(TimeInfo)) || (drifty.didChange(TimeInfo)))
|
596 : |
|
|
return 1;
|
597 : |
|
|
return 0;
|
598 : |
|
|
}
|
599 : |
|
|
|
600 : |
|
|
double MigrationFunction::getMigrationFunction(Rectangle* fromRec, Rectangle* toRec) {
|
601 : |
|
|
double fx, fy;
|
602 : |
|
|
double dx, dy;
|
603 : |
|
|
double xiL, yiL, xiU, yiU;
|
604 : |
|
|
double xfL, yfL, xfU, yfU;
|
605 : |
|
|
|
606 : |
|
|
if (isZero(diffusion) || isZero(lambda)) // prevent divide by zero errors ...
|
607 : |
|
|
return 0.0;
|
608 : |
|
|
|
609 : |
|
|
xiL = fromRec->getLowerX();
|
610 : |
|
|
yiL = fromRec->getLowerY();
|
611 : |
|
|
xiU = fromRec->getUpperX();
|
612 : |
|
|
yiU = fromRec->getUpperY();
|
613 : |
|
|
xfL = toRec->getLowerX();
|
614 : |
|
|
yfL = toRec->getLowerY();
|
615 : |
|
|
xfU = toRec->getUpperX();
|
616 : |
|
|
yfU = toRec->getUpperY();
|
617 : |
|
|
|
618 : |
|
|
dx = delta * diffusion;
|
619 : |
|
|
dy = dx * lambda;
|
620 : |
|
|
|
621 : |
|
|
// functions from Violeta ...
|
622 : |
|
|
fx = -f1x(xfU,xiU,dx,driftx) + f1x(xfU,xiL,dx,driftx) + f1x(xfL,xiU,dx,driftx) - f1x(xfL,xiL,dx,driftx) - f2x(xfU,xiU,dx,driftx) + f2x(xfU,xiL,dx,driftx) + f2x(xfL,xiU,dx,driftx) - f2x(xfL,xiL,dx,driftx);
|
623 : |
|
|
|
624 : |
|
|
if (fx < verysmall)
|
625 : |
|
|
return 0.0;
|
626 : |
|
|
|
627 : |
|
|
fy = -f1x(yfU,yiU,dy,drifty) + f1x(yfU,yiL,dy,drifty) + f1x(yfL,yiU,dy,drifty) - f1x(yfL,yiL,dy,drifty) - f2x(yfU,yiU,dy,drifty) + f2x(yfU,yiL,dy,drifty) + f2x(yfL,yiU,dy,drifty) - f2x(yfL,yiL,dy,drifty);
|
628 : |
|
|
|
629 : |
|
|
if (fy < verysmall)
|
630 : |
|
|
return 0.0;
|
631 : |
|
|
|
632 : |
|
|
return (0.5 * fx * fy);
|
633 : |
|
|
}
|
634 : |
|
|
|
635 : |
|
|
double MigrationFunction::f1x(double w, double u, double D, double beta) {
|
636 : |
|
|
return ((2/sqrt(pivalue))*sqrt(D)*exp(-(w-(beta*delta)-u)*(w-(beta*delta)-u)/(4*D)));
|
637 : |
|
|
}
|
638 : |
|
|
|
639 : |
|
|
double MigrationFunction::f2x(double w, double u, double D, double beta) {
|
640 : |
|
|
return ((w-(beta*delta)-u)*erf((w-(beta*delta)-u)/(2*sqrt(D))));
|
641 : |
|
|
}
|