#include "boundlikelihood.h"
#include "readfunc.h"
#include "readword.h"
#include "areatime.h"
#include "keeper.h"
#include "errorhandler.h"
#include "gadget.h"
#include "global.h"
BoundLikelihood::BoundLikelihood(CommentStream& infile, const AreaClass* const Area,
const TimeClass* const TimeInfo, const Keeper* const keeper, double weight, const char* name)
: Likelihood(BOUNDLIKELIHOOD, weight, name) {
int i, j;
Parameter tempParam;
double temp;
int count = 0;
//set flag to initialise the bounds - called in Reset
checkInitialised = 0;
infile >> ws;
if (countColumns(infile) != 4)
handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 4");
while (!infile.eof()) {
infile >> tempParam >> ws;
if (strcasecmp(tempParam.getName(), "default") == 0) {
infile >> defPower >> defLW >> defUW >> ws;
count++;
} else {
switches.resize(tempParam);
infile >> temp >> ws;
powers.resize(1, temp);
infile >> temp >> ws;
lowerweights.resize(1, temp);
infile >> temp >> ws;
upperweights.resize(1, temp);
switchnr.resize(1, -1);
count++;
}
infile >> ws;
}
handle.logMessage(LOGMESSAGE, "Read penalty file - number of entries", count);
}
void BoundLikelihood::Reset(const Keeper* const keeper) {
Likelihood::Reset(keeper);
handle.setNaNFlag(0); // reset the NaN count
if (isZero(weight))
handle.logMessage(LOGWARN, "Warning in boundlikelihood - zero weight for", this->getName());
if (!checkInitialised) {
if (!keeper->boundsGiven())
handle.logMessage(LOGWARN, "Warning in boundlikelihood - no bounds have been set in input file");
int i, j, k, numvar, numset, numfail;
numvar = keeper->numVariables();
numset = switches.Size();
if (numset != 0) {
numfail = 0;
ParameterVector sw(numvar);
keeper->getSwitches(sw);
for (i = 0; i < numset; i++)
for (j = 0; j < numvar; j++)
if (switches[i] == sw[j])
switchnr[i] = j;
for (i = 0; i < numset; i++) {
if (switchnr[i] == -1) {
handle.logMessage(LOGWARN, "Warning in boundlikelihood - failed to match switch", switches[i].getName());
numfail++;
// delete the entries for the non-existant switch
switches.Delete(i);
powers.Delete(i);
lowerweights.Delete(i);
upperweights.Delete(i);
switchnr.Delete(i);
if (numfail != numset)
i--;
}
}
numset -= numfail;
}
IntVector done(numset, 0);
// resize vectors to store data
likelihoods.resize(numvar, 0.0);
lowerbound.resize(numvar, 0.0);
upperbound.resize(numvar, 0.0);
values.resize(numvar, 0.0);
powers.resize(numvar - numset, 0.0);
lowerweights.resize(numvar - numset, 0.0);
upperweights.resize(numvar - numset, 0.0);
switchnr.resize(numvar - numset, -1);
DoubleVector lbs(numvar);
DoubleVector ubs(numvar);
keeper->getLowerBounds(lbs);
keeper->getUpperBounds(ubs);
k = 0;
for (i = 0; i < numvar; i++) {
if (switchnr[i] != -1) {
lowerbound[i] = lbs[switchnr[i]];
upperbound[i] = ubs[switchnr[i]];
if (i < numset)
done[i] = switchnr[i];
else
handle.logMessage(LOGFAIL, "Error in boundlikelihood - received invalid variable to check bounds");
} else {
for (j = 0; j < numset; j++)
if (k == done[j])
k++;
switchnr[i] = k;
lowerbound[i] = lbs[k];
upperbound[i] = ubs[k];
powers[i] = defPower;
lowerweights[i] = defLW;
upperweights[i] = defUW;
k++;
}
}
for (i = 0; i < powers.Size(); i++)
if (powers[i] < verysmall)
handle.logMessage(LOGFAIL, "Error in boundlikelihood - invalid value for power", powers[i]);
checkInitialised = 1;
}
if (handle.getLogLevel() >= LOGMESSAGE)
handle.logMessage(LOGMESSAGE, "Reset boundlikelihood component", this->getName());
}
void BoundLikelihood::addLikelihoodKeeper(const TimeClass* const TimeInfo, Keeper* const keeper) {
int i;
double temp;
keeper->getCurrentValues(values);
for (i = 0; i < switchnr.Size(); i++) {
if (values[switchnr[i]] < lowerbound[i]) {
temp = fabs(values[switchnr[i]] - lowerbound[i]);
if (temp < 1.0)
likelihoods[i] = lowerweights[i] * pow(temp, (1.0 / powers[i]));
else
likelihoods[i] = lowerweights[i] * pow(temp, powers[i]);
likelihood += likelihoods[i];
keeper->Update(switchnr[i], lowerbound[i]);
} else if (values[switchnr[i]] > upperbound[i]) {
temp = fabs(values[switchnr[i]] - upperbound[i]);
if (temp < 1.0)
likelihoods[i] = upperweights[i] * pow(temp, (1.0 / powers[i]));
else
likelihoods[i] = upperweights[i] * pow(temp, powers[i]);
likelihood += likelihoods[i];
keeper->Update(switchnr[i], upperbound[i]);
} else
likelihoods[i] = 0.0;
}
if ((handle.getLogLevel() >= LOGMESSAGE) && (isZero(likelihood)))
handle.logMessage(LOGMESSAGE, "For this model simulation, no parameters are outside the bounds");
if (handle.getNaNFlag()) {
likelihood += verybig;
if (handle.getLogLevel() >= LOGMESSAGE)
handle.logMessage(LOGMESSAGE, "For this model simulation, a NaN was found within the model");
}
if (handle.getLogLevel() >= LOGMESSAGE)
handle.logMessage(LOGMESSAGE, "Calculated likelihood score for boundlikelihood component to be", likelihood);
}
void BoundLikelihood::printSummary(ofstream& outfile) {
//JMB there is only one likelihood score here ...
if (!(isZero(likelihood))) {
outfile << "all all all" << sep << setw(largewidth) << this->getName() << sep
<< setprecision(smallprecision) << setw(smallwidth) << weight << sep
<< setprecision(largeprecision) << setw(largewidth) << likelihood << endl;
outfile.flush();
}
}
void BoundLikelihood::Print(ofstream& outfile) const {
outfile << "\nBoundlikelihood " << this->getName() << " - likelihood value "
<< likelihood << endl;
outfile.flush();
}