#include "quotapredator.h" #include "keeper.h" #include "prey.h" #include "mathfunc.h" #include "errorhandler.h" #include "gadget.h" #include "global.h" QuotaPredator::QuotaPredator(CommentStream& infile, const char* givenname, const IntVector& Areas, const TimeClass* const TimeInfo, Keeper* const keeper, Formula multscaler) : LengthPredator(givenname, Areas, keeper, multscaler) { type = QUOTAPREDATOR; keeper->addString("predator"); keeper->addString(givenname); //first read in the suitability parameters this->readSuitability(infile, TimeInfo, keeper); int i, j; double tmp; char text[MaxStrLength]; strncpy(text, "", MaxStrLength); functionnumber = 0; functionname = new char[MaxStrLength]; strncpy(functionname, "", MaxStrLength); //the next entry in the input file should be the name of the quota function infile >> functionname >> ws; if (strcasecmp(functionname, "simple") == 0) functionnumber = 1; else if (strcasecmp(functionname, "simplesum") == 0) functionnumber = 2; else if (strcasecmp(functionname, "simpleselect") == 0) functionnumber = 3; else if (strcasecmp(functionname, "annual") == 0) functionnumber = 4; else if (strcasecmp(functionname, "annualsum") == 0) functionnumber = 5; else if (strcasecmp(functionname, "annualselect") == 0) functionnumber = 6; else handle.logFileMessage(LOGFAIL, "\nError in quotapredator - unrecognised function", functionname); //now read in the stock biomass levels infile >> text >> ws; if (strcasecmp(text, "biomasslevel") != 0) handle.logFileUnexpected(LOGFAIL, "biomasslevel", text); while (isdigit(infile.peek()) && !infile.eof()) { infile >> tmp >> ws; biomasslevel.resize(1, tmp); } if (biomasslevel.Size() == 0) handle.logMessage(LOGFAIL, "Error in quotapredator - missing quota data"); //check the data to make sure that it is continuous and positive if (biomasslevel[0] < 0.0) handle.logFileMessage(LOGFAIL, "biomass levels must be positive"); if (biomasslevel.Size() != 1) for (i = 1; i < biomasslevel.Size(); i++) if ((biomasslevel[i] - biomasslevel[i - 1]) < verysmall) handle.logFileMessage(LOGFAIL, "biomass levels must be strictly increasing"); //finally read in the quota parameters infile >> text >> ws; if (strcasecmp(text, "quotalevel") != 0) handle.logFileUnexpected(LOGFAIL, "quotalevel", text); infile >> ws; keeper->addString("quota"); quotalevel.resize(biomasslevel.Size() + 1, keeper); for (i = 0; i < quotalevel.Size(); i++) if (!(infile >> quotalevel[i])) handle.logFileMessage(LOGFAIL, "invalid format of quotalevel vector"); quotalevel.Inform(keeper); keeper->clearLast(); calcquota.resize(preference.Size(), 0.0); selectprey.resize(preference.Size(), 0); //default to not using the prey //if we need to select the stocks used to calculate the fishing quota infile >> text >> ws; if ((functionnumber == 3) || (functionnumber == 6)) { if (strcasecmp(text, "selectstocks") != 0) handle.logFileUnexpected(LOGFAIL, "selectstocks", text); i = 0; CharPtrVector preynames; infile >> text >> ws; while (!infile.eof() && (strcasecmp(text, "amount") != 0)) { preynames.resize(new char[strlen(text) + 1]); strcpy(preynames[i++], text); infile >> text >> ws; } if (preynames.Size() == 0) handle.logFileMessage(LOGFAIL, "no stocks selected to calculate quota"); //check that the prey names are unique for (i = 0; i < preynames.Size(); i++) for (j = 0; j < preynames.Size(); j++) if ((strcasecmp(preynames[i], preynames[j]) == 0) && (i != j)) handle.logFileMessage(LOGFAIL, "repeated stock", preynames[i]); //work out which preys are needed to calculate the fishing quota int check = 0; for (i = 0; i < preynames.Size(); i++) { check = 0; for (j = 0; j < preference.Size(); j++) { if (strcasecmp(preynames[i], this->getPreyName(j)) == 0) { selectprey[j] = 1; check = 1; } } if (check == 0) handle.logFileMessage(LOGFAIL, "failed to match stock", preynames[i]); } //finally free the memory since the prey names are no longer needed for (i = 0; i < preynames.Size(); i++) delete[] preynames[i]; } if (strcasecmp(text, "amount") != 0) handle.logFileUnexpected(LOGFAIL, "amount", text); keeper->clearLast(); keeper->clearLast(); } QuotaPredator::~QuotaPredator() { delete[] functionname; } void QuotaPredator::Eat(int area, const AreaClass* const Area, const TimeClass* const TimeInfo) { int inarea = this->areaNum(area); int prey, preyl; int predl = 0; //JMB there is only ever one length group ... double tmp, bio; totalcons[inarea][predl] = 0.0; tmp = prednumber[inarea][predl].N * multi * TimeInfo->getTimeStepSize() / TimeInfo->numSubSteps(); if (isZero(tmp)) { //JMB no predation takes place on this timestep for (prey = 0; prey < this->numPreys(); prey++) (*predratio[inarea])[prey][predl] = 0.0; return; } switch (functionnumber) { case 1: //Calculate the fishing level based on the available biomass of each stock for (prey = 0; prey < this->numPreys(); prey++) { if (this->getPrey(prey)->isPreyArea(area)) { (*predratio[inarea])[prey][predl] = tmp * calcQuota(this->getPrey(prey)->getTotalBiomass(area)); if ((*predratio[inarea])[prey][predl] > 10.0) //JMB arbitrary value here ... handle.logMessage(LOGWARN, "Warning in quotapredator - excessive consumption required"); } else (*predratio[inarea])[prey][predl] = 0.0; } break; case 2: //Calculate the fishing level based on the available biomass of all stocks bio = 0.0; for (prey = 0; prey < this->numPreys(); prey++) if (this->getPrey(prey)->isPreyArea(area)) bio += this->getPrey(prey)->getTotalBiomass(area); tmp *= calcQuota(bio); if (tmp > 10.0) //JMB arbitrary value here ... handle.logMessage(LOGWARN, "Warning in quotapredator - excessive consumption required"); for (prey = 0; prey < this->numPreys(); prey++) { if (this->getPrey(prey)->isPreyArea(area)) (*predratio[inarea])[prey][predl] = tmp; else (*predratio[inarea])[prey][predl] = 0.0; } break; case 3: //Calculate the fishing level based on the available biomass of the selected stocks bio = 0.0; for (prey = 0; prey < this->numPreys(); prey++) if ((this->getPrey(prey)->isPreyArea(area)) && (selectprey[prey])) bio += this->getPrey(prey)->getTotalBiomass(area); tmp *= calcQuota(bio); if (tmp > 10.0) //JMB arbitrary value here ... handle.logMessage(LOGWARN, "Warning in quotapredator - excessive consumption required"); for (prey = 0; prey < this->numPreys(); prey++) { if (this->getPrey(prey)->isPreyArea(area)) (*predratio[inarea])[prey][predl] = tmp; else (*predratio[inarea])[prey][predl] = 0.0; } break; case 4: //Calculate the annual fishing level based on the available biomass of each stock for (prey = 0; prey < this->numPreys(); prey++) { if (this->getPrey(prey)->isPreyArea(area)) { if (TimeInfo->getStep() == 1) calcquota[prey] = calcQuota(this->getPrey(prey)->getTotalBiomass(area)); //JMB this doesnt work if there is more than one stock since the value of quota will have changed (*predratio[inarea])[prey][predl] = calcquota[prey] * tmp; if ((*predratio[inarea])[prey][predl] > 10.0) //JMB arbitrary value here ... handle.logMessage(LOGWARN, "Warning in quotapredator - excessive consumption required"); } else (*predratio[inarea])[prey][predl] = 0.0; } break; case 5: //Calculate the annual fishing level based on the available biomass of all stocks if (TimeInfo->getStep() == 1) { bio = 0.0; for (prey = 0; prey < this->numPreys(); prey++) if (this->getPrey(prey)->isPreyArea(area)) bio += this->getPrey(prey)->getTotalBiomass(area); calcquota[0] = calcQuota(bio); //JMB all quotas are the same } tmp *= calcquota[0]; if (tmp > 10.0) //JMB arbitrary value here ... handle.logMessage(LOGWARN, "Warning in quotapredator - excessive consumption required"); for (prey = 0; prey < this->numPreys(); prey++) { if (this->getPrey(prey)->isPreyArea(area)) (*predratio[inarea])[prey][predl] = tmp; else (*predratio[inarea])[prey][predl] = 0.0; } break; case 6: //Calculate the annual fishing level based on the available biomass of the selected stocks if (TimeInfo->getStep() == 1) { bio = 0.0; for (prey = 0; prey < this->numPreys(); prey++) if ((this->getPrey(prey)->isPreyArea(area)) && (selectprey[prey])) bio += this->getPrey(prey)->getTotalBiomass(area); calcquota[0] = calcQuota(bio); //JMB all quotas are the same } tmp *= calcquota[0]; if (tmp > 10.0) //JMB arbitrary value here ... handle.logMessage(LOGWARN, "Warning in quotapredator - excessive consumption required"); for (prey = 0; prey < this->numPreys(); prey++) { if (this->getPrey(prey)->isPreyArea(area)) (*predratio[inarea])[prey][predl] = tmp; else (*predratio[inarea])[prey][predl] = 0.0; } break; default: handle.logMessage(LOGWARN, "Warning in quotapredator - unrecognised function", functionname); break; } for (prey = 0; prey < this->numPreys(); prey++) { if (isZero((*predratio[inarea])[prey][predl])) { for (preyl = 0; preyl < (*cons[inarea][prey])[predl].Size(); preyl++) (*cons[inarea][prey])[predl][preyl] = 0.0; } else { for (preyl = 0; preyl < (*cons[inarea][prey])[predl].Size(); preyl++) { (*cons[inarea][prey])[predl][preyl] = (*predratio[inarea])[prey][predl] * this->getSuitability(prey)[predl][preyl] * this->getPrey(prey)->getBiomass(area, preyl); totalcons[inarea][predl] += (*cons[inarea][prey])[predl][preyl]; } //inform the preys of the consumption this->getPrey(prey)->addBiomassConsumption(area, (*cons[inarea][prey])[predl]); } } } void QuotaPredator::adjustConsumption(int area, const TimeClass* const TimeInfo) { int prey, preyl; int inarea = this->areaNum(area); int predl = 0; //JMB there is only ever one length group ... overcons[inarea][predl] = 0.0; if (isZero(totalcons[inarea][predl])) //JMB no predation takes place on this timestep return; double maxRatio, tmp; maxRatio = TimeInfo->getMaxRatioConsumed(); for (prey = 0; prey < this->numPreys(); prey++) { if (this->getPrey(prey)->isOverConsumption(area)) { hasoverconsumption[inarea] = 1; DoubleVector ratio = this->getPrey(prey)->getRatio(area); for (preyl = 0; preyl < (*cons[inarea][prey])[predl].Size(); preyl++) { if (ratio[preyl] > maxRatio) { tmp = maxRatio / ratio[preyl]; overcons[inarea][predl] += (1.0 - tmp) * (*cons[inarea][prey])[predl][preyl]; (*cons[inarea][prey])[predl][preyl] *= tmp; (*usesuit[inarea][prey])[predl][preyl] *= tmp; } } } } if (hasoverconsumption[inarea]) { totalcons[inarea][predl] -= overcons[inarea][predl]; overconsumption[inarea][predl] += overcons[inarea][predl]; } totalconsumption[inarea][predl] += totalcons[inarea][predl]; for (prey = 0; prey < this->numPreys(); prey++) if (this->getPrey(prey)->isPreyArea(area)) for (preyl = 0; preyl < (*cons[inarea][prey])[predl].Size(); preyl++) (*consumption[inarea][prey])[predl][preyl] += (*cons[inarea][prey])[predl][preyl]; } void QuotaPredator::Print(ofstream& outfile) const { outfile << "QuotaPredator\n"; PopPredator::Print(outfile); } double QuotaPredator::calcQuota(double biomass) { int i; double quota = 0.0; if (biomass < biomasslevel[0]) { quota = quotalevel[0]; } else if (biomass > biomasslevel[biomasslevel.Size() - 1]) { quota = quotalevel[biomasslevel.Size()]; } else { for (i = 1; i < biomasslevel.Size(); i++) if ((biomasslevel[i - 1] < biomass) && (biomass < biomasslevel[i])) quota = quotalevel[i]; } if (quota < 0.0) handle.logMessage(LOGWARN, "Warning in quotapredator - negative quota", quota); return quota; }