#include "migration.h"
#include "readfunc.h"
#include "readword.h"
#include "errorhandler.h"
#include "mathfunc.h"
#include "migrationarea.h"
#include "gadget.h"
#include "global.h"
// ********************************************************
// Functions for base Migration
// ********************************************************
Migration::Migration(const IntVector& Areas, const char* givenname)
: HasName(givenname), LivesOnAreas(Areas) {
if (areas.Size() == 1)
handle.logMessage(LOGWARN, "Warning in migration - only one area defined");
}
// ********************************************************
// Functions for MigrationNumbers
// ********************************************************
MigrationNumbers::MigrationNumbers(CommentStream& infile, const IntVector& Areas,
const AreaClass* const Area, const TimeClass* const TimeInfo,
const char* givenname, Keeper* const keeper)
: Migration(Areas, givenname) {
int i, j;
ifstream subfile;
CommentStream subcomment(subfile);
char text[MaxStrLength];
char filename[MaxStrLength];
strncpy(text, "", MaxStrLength);
strncpy(filename, "", MaxStrLength);
keeper->addString("migrationmatrix");
// Open and read year and step information
readWordAndValue(infile, "yearstepfile", filename);
subfile.open(filename, ios::in);
handle.checkIfFailure(subfile, filename);
handle.Open(filename);
this->readTimeStepData(subcomment, TimeInfo);
handle.Close();
subfile.close();
subfile.clear();
// Open and read information about how the matrices are defined
infile >> text >> filename >> ws;
subfile.open(filename, ios::in);
handle.checkIfFailure(subfile, filename);
handle.Open(filename);
if (strcasecmp(text, "definematrices") == 0) {
this->readGivenMatrices(subcomment, keeper);
} else if (strcasecmp(text, "defineratios") == 0) {
this->readGivenRatios(subcomment, keeper, Area);
} else
handle.logFileUnexpected(LOGFAIL, "definematrices or defineratios", text);
handle.Close();
subfile.close();
subfile.clear();
this->checkMatrixIndex();
handle.logMessage(LOGMESSAGE, "Read migration numbers file - number of migration matrices", readMigration.Size());
keeper->clearLast();
}
MigrationNumbers::~MigrationNumbers() {
int i;
for (i = 0; i < timeindex.Size(); i++)
if (timeindex[i] != -1)
delete[] allmatrixnames[i];
for (i = 0; i < usedmatrixnames.Size(); i++) {
delete[] usedmatrixnames[i];
delete readMigration[i];
delete calcMigration[i];
}
}
void MigrationNumbers::readTimeStepData(CommentStream& infile, const TimeClass* const TimeInfo) {
int year, step, timeid, count;
char text[MaxStrLength];
strncpy(text, "", MaxStrLength);
// check the number of columns in the inputfile
infile >> ws;
if (countColumns(infile) != 3)
handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 3");
year = step = count = 0;
allmatrixnames.resizeBlank(TimeInfo->numTotalSteps() + 1);
timeindex.resize(TimeInfo->numTotalSteps() + 1, -1);
while (!infile.eof()) {
infile >> year >> step >> text >> ws;
if (TimeInfo->isWithinPeriod(year, step)) {
timeid = TimeInfo->calcSteps(year, step);
timeindex[timeid] = 0;
allmatrixnames[timeid] = new char[strlen(text) + 1];
strcpy(allmatrixnames[timeid], text);
count++;
}
}
if (count == 0)
handle.logMessage(LOGWARN, "Warning in migration - found no migration data");
if (count != TimeInfo->numTotalSteps())
handle.logMessage(LOGWARN, "Warning in migration - migration data doesnt span time range");
handle.logMessage(LOGMESSAGE, "Read migration time data - number of entries", count);
}
void MigrationNumbers::readGivenRatios(CommentStream& infile, Keeper* const keeper, const AreaClass* const Area) {
int innerrow, innercol, fromArea, toArea;
int i, j, k, count;
int numAreas = areas.Size();
char text[MaxStrLength];
char name[MaxStrLength];
strncpy(text, "", MaxStrLength);
strncpy(name, "", MaxStrLength);
/* File format:
[migrationmatrix]
name matrixname
fromarea toarea ratio
[migrationmatrix]
...
*/
toArea = fromArea = 0;
// Read ratio information from file
while (!infile.eof() && !infile.fail()) {
infile >> text >> ws;
if (strcasecmp(text, "[migrationmatrix]") != 0)
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix]", text);
readWordAndValue(infile, "name", name);
if (this->useMatrix(name)) {
this->setMatrixName(name);
count = readMigration.Size();
checkvalues.AddRows(1, numAreas, 0);
readMigration.resize(new FormulaMatrix(numAreas, numAreas, 0.0));
infile >> ws;
while (!infile.eof() && !infile.fail() && (infile.peek() != '[')) {
if (!(isdigit(infile.peek()))) {
infile >> text >> ws;
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix] or area number", text);
}
infile >> fromArea >> toArea >> ws;
innerrow = Area->getInnerArea(toArea);
innercol = Area->getInnerArea(fromArea);
if (!this->isInArea(innerrow) || !this->isInArea(innercol))
handle.logMessage(LOGFAIL, "Error in migration - invalid area for matrix", name);
toArea = this->areaNum(innerrow);
fromArea = this->areaNum(innercol);
infile >> (*readMigration[count])[toArea][fromArea] >> ws;
checkvalues[count][fromArea]++;
}
} else {
handle.logMessage(LOGWARN, "Warning in migration - matrix not used in the simulation", name);
infile >> ws;
while (!infile.eof() && !infile.fail() && infile.peek() != '[')
infile >> text >> ws;
}
}
// Finished reading from file, now check the values
Formula missing;
missing.setValue(1.0);
for (i = 0; i < readMigration.Size(); i++) {
for (j = 0; j < readMigration[i]->Nrow(); j++) {
if (checkvalues[i][j] == 0) {
// no data read for migration from this area
(*readMigration[i])[j][j].setValue(1.0);
} else if (checkvalues[i][j] == numAreas) {
// nothing to do - data read for migration from this area to every other area
} else if (isZero((*readMigration[i])[j][j])) {
// read data for migration from this area to other areas
// so simply subtract the migration from 1 to find those that are left
count = 0;
vector<Formula*> checklist;
for (k = 0; k < readMigration[i]->Ncol(); k++) {
if (!(isZero((*readMigration[i])[k][j]))) {
checklist.push_back(&(*readMigration[i])[k][j]);
count++;
}
}
if (count == 0)
handle.logFileMessage(LOGFAIL, "missing data for migration matrix", usedmatrixnames[i]);
checklist.insert(checklist.begin(), &missing);
Formula* missingValue = new Formula(MINUS, checklist);
(*readMigration[i])[j][j] = *missingValue;
delete missingValue;
} else if (checkvalues[i][j] == (numAreas - 1)) {
// read migration data from this area to all other areas but 1
count = -1;
vector<Formula*> checklist;
for (k = 0; k < readMigration[i]->Ncol(); k++) {
if (isZero((*readMigration[i])[k][j]))
count = k;
else
checklist.push_back(&(*readMigration[i])[k][j]);
}
if (count == -1)
handle.logFileMessage(LOGFAIL, "missing data for migration matrix", usedmatrixnames[i]);
checklist.insert(checklist.begin(), &missing);
Formula* missingValue = new Formula(MINUS, checklist);
(*readMigration[i])[count][j] = *missingValue;
delete missingValue;
} else {
// havent read enough information to reliably construct matrix
handle.logFileMessage(LOGFAIL, "missing data for migration matrix", usedmatrixnames[i]);
}
}
}
// Inform keeper of the values and resize
for (i = 0; i < readMigration.Size(); i++) {
(*readMigration[i]).Inform(keeper);
calcMigration.resize(new DoubleMatrix(numAreas, numAreas, 0.0));
}
}
void MigrationNumbers::readGivenMatrices(CommentStream& infile, Keeper* const keeper) {
int i, j, count;
int numAreas = areas.Size();
char text[MaxStrLength];
char name[MaxStrLength];
strncpy(text, "", MaxStrLength);
strncpy(name, "", MaxStrLength);
/* File format:
[migrationmatrix]
name matrixname
x11 ... x1n
...
xn1 ... xnn
[migrationmatrix]
...
*/
// Read matrix information from file
while (!infile.eof() && !infile.fail()) {
infile >> text >> ws;
if (strcasecmp(text, "[migrationmatrix]") != 0)
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix]", text);
readWordAndValue(infile, "name", name);
if (this->useMatrix(name)) {
this->setMatrixName(name);
count = readMigration.Size();
readMigration.resize(new FormulaMatrix(numAreas, numAreas, 0.0));
infile >> ws;
if (infile.peek() == '[' || infile.eof()) {
for (i = 0; i < numAreas; i++)
(*readMigration[count])[i][i].setValue(1.0);
} else if (isdigit(infile.peek()) || infile.peek() == '#' || infile.peek() == '(') {
// look for a number or a formula as the next entry to be read
for (i = 0; i < numAreas; i++)
for (j = 0; j < numAreas; j++)
if (!(infile >> (*readMigration[count])[i][j]))
handle.logMessage(LOGFAIL, "Error in migration - failed to read matrix", name);
infile >> ws; //JMB strip whitespace
} else {
infile >> text;
handle.logFileUnexpected(LOGFAIL, "[migrationmatrix] or matrix value", text);
}
} else {
handle.logMessage(LOGWARN, "Warning in migration - matrix not used in the simulation", name);
infile >> ws;
while (!infile.eof() && !infile.fail() && infile.peek() != '[')
infile >> text >> ws;
}
}
// Inform keeper of the values and resize
for (i = 0; i < readMigration.Size(); i++) {
(*readMigration[i]).Inform(keeper);
calcMigration.resize(new DoubleMatrix(numAreas, numAreas, 0.0));
}
}
void MigrationNumbers::checkMatrixIndex() {
int i, j, check;
for (i = 0; i < timeindex.Size(); i++) {
if (timeindex[i] == 0) {
check = -1;
for (j = 0; j < usedmatrixnames.Size(); j++)
if (strcasecmp(allmatrixnames[i], usedmatrixnames[j]) == 0)
check = j;
if (check == -1)
handle.logMessage(LOGFAIL, "Error in migration - failed to read matrix", allmatrixnames[i]);
timeindex[i] = check;
} else if (timeindex[i] != -1) {
handle.logMessage(LOGFAIL, "Error in migration - repeated migration matrix", i);
}
}
}
void MigrationNumbers::Print(ofstream& outfile) {
int i, j, k;
outfile << "\nMigration\n\tNames of migration matrices:\n\t";
for (i = 0; i < timeindex.Size(); i++)
if (timeindex[i] != -1)
outfile << allmatrixnames[i] << sep;
outfile << "\n\n\tMigration matrices";
for (i = 0; i < calcMigration.Size(); i++) {
outfile << "\n\tMatrix name: " << usedmatrixnames[i] << "\n\t";
for (j = 0; j < calcMigration[i]->Nrow(); j++) {
for (k = 0; k < calcMigration[i]->Ncol(j); k++)
outfile << setw(smallwidth) << (*calcMigration[i])[j][k] << sep;
outfile << "\n\t";
}
}
outfile.flush();
}
const DoubleMatrix& MigrationNumbers::getMigrationMatrix(const TimeClass* const TimeInfo) {
return (*calcMigration[timeindex[TimeInfo->getTime()]]);
}
void MigrationNumbers::setMatrixName(char* name) {
int i;
// check to ensure that this matrix name is unique
for (i = 0; i < usedmatrixnames.Size(); i++)
if (strcasecmp(usedmatrixnames[i], name) == 0)
handle.logMessage(LOGFAIL, "Error in migration - repeated matrix", name);
// store the new matrix name
char* tempname = new char[strlen(name) + 1];
strcpy(tempname, name);
usedmatrixnames.resize(tempname);
}
void MigrationNumbers::Reset() {
//JMB need to reset the penalty vector first
penalty.Reset();
int i, j, k;
double colsum;
for (i = 0; i < readMigration.Size(); i++) {
for (k = 0; k < readMigration[i]->Ncol(); k++) {
colsum = 0.0;
for (j = 0; j < readMigration[i]->Nrow(); j++) {
if (isZero((*readMigration[i])[j][k])) {
(*calcMigration[i])[j][k] = 0.0;
} else if ((*readMigration[i])[j][k] < 0.0) {
penalty.resize(1, (*readMigration[i])[j][k]);
(*calcMigration[i])[j][k] = 0.0;
handle.logMessage(LOGWARN, "Warning in migration - value outside bounds", (*readMigration[i])[j][k]);
} else if ((*readMigration[i])[j][k] > 1.0) {
penalty.resize(1, (*readMigration[i])[j][k]);
(*calcMigration[i])[j][k] = 1.0;
handle.logMessage(LOGWARN, "Warning in migration - value outside bounds", (*readMigration[i])[j][k]);
} else
(*calcMigration[i])[j][k] = (*readMigration[i])[j][k];
colsum += (*calcMigration[i])[j][k];
}
if (isZero(colsum)) {
handle.logMessage(LOGWARN, "Warning in migration - column sum is zero");
// just putting anything here ... AJ 23.06.05, Needs to be changed!!!!!
penalty.resize(1, 1.0);
for (j = 0; j < readMigration[i]->Nrow(); j++) {
if (k == j)
(*calcMigration[i])[j][k] = 1.0;
else
(*calcMigration[i])[j][k] = 0.0;
}
} else if (isEqual(colsum, 1.0)) {
// everything is OK ...
} else {
penalty.resize(1, colsum);
colsum = 1.0 / colsum;
for (j = 0; j < readMigration[i]->Nrow(); j++)
(*calcMigration[i])[j][k] *= colsum;
}
}
}
if (handle.getLogLevel() >= LOGMESSAGE)
handle.logMessage(LOGMESSAGE, "Reset migration data for stock", this->getName());
}
int MigrationNumbers::useMatrix(char* name) {
int i;
for (i = 0; i < timeindex.Size(); i++)
if ((timeindex[i] != -1) && (strcasecmp(allmatrixnames[i], name) == 0))
return 1;
return 0;
}
int MigrationNumbers::isMigrationStep(const TimeClass* const TimeInfo) {
if (timeindex[TimeInfo->getTime()] == -1)
return 0;
return 1;
}
// ********************************************************
// Functions for MigrationFunction
// ********************************************************
MigrationFunction::MigrationFunction(CommentStream& infile, const IntVector& Areas,
const AreaClass* const Area, const TimeClass* const TimeInfo,
const char* givenname, Keeper* const keeper)
: Migration(Areas, givenname) {
int i;
ifstream subfile;
CommentStream subcomment(subfile);
char text[MaxStrLength];
char filename[MaxStrLength];
strncpy(text, "", MaxStrLength);
strncpy(filename, "", MaxStrLength);
keeper->addString("migrationfunction");
readWordAndModelVariable(infile, "diffusion", diffusion, TimeInfo, keeper);
readWordAndModelVariable(infile, "driftx", driftx, TimeInfo, keeper);
readWordAndModelVariable(infile, "drifty", drifty, TimeInfo, keeper);
readWordAndVariable(infile, "lambda", lambda);
delta = TimeInfo->getTimeStepLength() / TimeInfo->numSubSteps();
// Open and read given year and step information
readWordAndValue(infile, "areadefinition", filename);
subfile.open(filename, ios::in);
handle.checkIfFailure(subfile, filename);
handle.Open(filename);
this->readAreaData(subcomment, Area, TimeInfo, keeper);
handle.Close();
subfile.close();
subfile.clear();
calcMigration.AddRows(oceanareas.Size(), oceanareas.Size(), 0.0);
if (oceanareas.Size() != areas.Size())
handle.logMessage(LOGWARN, "Warning in migration - migration data doesnt cover all areas");
handle.logMessage(LOGMESSAGE, "Read migration function file - number of migration areas", oceanareas.Size());
keeper->clearLast();
}
void MigrationFunction::readAreaData(CommentStream& infile, const AreaClass* const Area,
const TimeClass* const TimeInfo, Keeper* const keeper) {
int area, inarea;
ifstream subfile;
CommentStream subcomment(subfile);
char text[MaxStrLength];
char filename[MaxStrLength];
char areaname[MaxStrLength];
strncpy(text, "", MaxStrLength);
strncpy(filename, "", MaxStrLength);
strncpy(areaname, "", MaxStrLength);
/* File format:
[area]
name areaname
number areanumber
rectangles filename containing the rectangle information
[area]
...
*/
// Read area information from file
infile >> ws;
while (!infile.eof() && !infile.fail()) {
infile >> text >> ws;
if (strcasecmp(text, "[area]") != 0)
handle.logFileUnexpected(LOGFAIL, "[area]", text);
readWordAndValue(infile, "name", areaname);
readWordAndVariable(infile, "number", area);
readWordAndValue(infile, "rectangles", filename);
inarea = Area->getInnerArea(area);
if (!this->isInArea(inarea))
handle.logMessage(LOGFAIL, "Error in migration - invalid area", area);
subfile.open(filename, ios::in);
handle.checkIfFailure(subfile, filename);
handle.Open(filename);
oceanareas.resize(new MigrationArea(subcomment, areaname, inarea));
handle.Close();
subfile.close();
subfile.clear();
}
}
MigrationFunction::~MigrationFunction() {
int i;
for (i = 0; i < oceanareas.Size(); i++)
delete oceanareas[i];
}
void MigrationFunction::Print(ofstream& outfile) {
int i, j;
outfile << "\nMigration\n\n\tMigration matrix\n\t";
for (i = 0; i < calcMigration.Nrow(); i++) {
for (j = 0; j < calcMigration.Ncol(i); j++)
outfile << setw(smallwidth) << calcMigration[i][j] << sep;
outfile << "\n\t";
}
outfile.flush();
}
const DoubleMatrix& MigrationFunction::getMigrationMatrix(const TimeClass* const TimeInfo) {
if (this->updateVariables(TimeInfo))
this->recalcMatrix();
return calcMigration;
}
void MigrationFunction::recalcMatrix() {
int from, to, idfrom, idto, i, j;
double mig, sa;
double colsum, sum;
for (from = 0; from < oceanareas.Size(); from++) {
idfrom = this->areaNum(oceanareas[from]->getAreaID());
colsum = 0.0;
for (to = 0; to < oceanareas.Size(); to++) {
idto = this->areaNum(oceanareas[to]->getAreaID());
sum = 0.0;
for (i = 0; i < oceanareas[from]->getNumRectangles(); i++) {
sa = (oceanareas[from]->getRectangles()[i])->getArea();
if (!(isZero(sa))) {
for (j = 0; j < oceanareas[to]->getNumRectangles(); j++) {
mig = this->getMigrationFunction(oceanareas[from]->getRectangles()[i], oceanareas[to]->getRectangles()[j]);
sum += (mig * sa);
}
}
}
calcMigration[idto][idfrom] = sum / oceanareas[from]->getArea();
colsum += calcMigration[idto][idfrom];
}
if (isZero(colsum)) {
handle.logMessage(LOGWARN, "Warning in migration - column sum is zero");
penalty.resize(1, 1.0);
calcMigration[idfrom][idfrom] = 1.0;
} else if (isEqual(colsum, 1.0)) {
// everything is OK ...
} else {
handle.logMessage(LOGWARN, "Warning in migration - column sum", colsum);
penalty.resize(1, colsum);
colsum = 1.0 / colsum;
for (i = 0; i < calcMigration.Ncol(); i++)
calcMigration[i][idfrom] *= colsum;
}
}
}
void MigrationFunction::Reset() {
//JMB need to reset the penalty vector
penalty.Reset();
if (handle.getLogLevel() >= LOGMESSAGE)
handle.logMessage(LOGMESSAGE, "Reset migration data for stock", this->getName());
}
int MigrationFunction::isMigrationStep(const TimeClass* const TimeInfo) {
return 1;
}
int MigrationFunction::updateVariables(const TimeClass* const TimeInfo) {
//update the values of the variables that can change
delta = TimeInfo->getTimeStepLength() / TimeInfo->numSubSteps();
diffusion.Update(TimeInfo);
driftx.Update(TimeInfo);
drifty.Update(TimeInfo);
if ((TimeInfo->didStepSizeChange()) || (diffusion.didChange(TimeInfo))
|| (driftx.didChange(TimeInfo)) || (drifty.didChange(TimeInfo)))
return 1;
return 0;
}
double MigrationFunction::getMigrationFunction(Rectangle* fromRec, Rectangle* toRec) {
double fx, fy;
double dx, dy;
double xiL, yiL, xiU, yiU;
double xfL, yfL, xfU, yfU;
if (isZero(diffusion) || isZero(lambda)) // prevent divide by zero errors ...
return 0.0;
xiL = fromRec->getLowerX();
yiL = fromRec->getLowerY();
xiU = fromRec->getUpperX();
yiU = fromRec->getUpperY();
xfL = toRec->getLowerX();
yfL = toRec->getLowerY();
xfU = toRec->getUpperX();
yfU = toRec->getUpperY();
dx = delta * diffusion;
dy = dx * lambda;
// functions from Violeta ...
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);
if (fx < verysmall)
return 0.0;
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);
if (fy < verysmall)
return 0.0;
return (0.5 * fx * fy);
}
double MigrationFunction::f1x(double w, double u, double D, double beta) {
return ((2/sqrt(pivalue))*sqrt(D)*exp(-(w-(beta*delta)-u)*(w-(beta*delta)-u)/(4*D)));
}
double MigrationFunction::f2x(double w, double u, double D, double beta) {
return ((w-(beta*delta)-u)*erf((w-(beta*delta)-u)/(2*sqrt(D))));
}
Parent Directory
Revision Log