#include "keeper.h"
#include "errorhandler.h"
#include "runid.h"
#include "ecosystem.h"
#include "gadget.h"
#include "global.h"
extern Ecosystem* EcoSystem;
#ifdef _OPENMP
extern Ecosystem** EcoSystems;
#endif
Keeper::Keeper() {
stack = new StrStack();
boundsgiven = 0;
fileopen = 0;
numoptvar = 0;
bestlikelihood = 0.0;
}
void Keeper::keepVariable(double& value, Parameter& attr) {
int i, index = -1;
for (i = 0; i < switches.Size(); i++)
if (switches[i] == attr)
index = i;
if (index == -1) {
//attr was not found -- add it to switches and values
index = switches.Size();
switches.resize(attr);
values.resize(1, value);
bestvalues.resize(1, value);
lowerbds.resize(1, -9999.0); // default lower bound
upperbds.resize(1, 9999.0); // default upper bound
opt.resize(1, 0);
scaledvalues.resize(1, 1.0);
initialvalues.resize(1, 1.0);
address.resize();
address[index].resize();
address[index][0] = &value;
if (stack->getSize() != 0)
address[index][0] = stack->sendAll();
} else {
if (value != values[index]) {
handle.logFileMessage(LOGFAIL, "read repeated switch name but different initial value", switches[index].getName());
} else {
i = address[index].Size();
address[index].resize();
address[index][i] = &value;
if (stack->getSize() != 0)
address[index][i] = stack->sendAll();
}
}
}
Keeper::~Keeper() {
delete stack;
if (fileopen) {
handle.Close();
outfile.close();
outfile.clear();
}
}
void Keeper::deleteParameter(const double& var) {
int i, j, check;
check = 0;
for (i = 0; i < address.Nrow(); i++) {
for (j = 0; j < address[i].Size(); j++) {
if (address[i][j] == &var) {
check++;
address[i].Delete(j);
if (address[i].Size() == 0) {
//the variable we deleted was the only one with this switch
address.Delete(i);
switches.Delete(i);
values.Delete(i);
bestvalues.Delete(i);
opt.Delete(i);
lowerbds.Delete(i);
upperbds.Delete(i);
scaledvalues.Delete(i);
initialvalues.Delete(i);
i--;
}
}
}
}
if (check != 1)
handle.logMessage(LOGFAIL, "Error in keeper - failed to delete parameter");
}
void Keeper::changeVariable(const double& pre, double& post) {
int i, j, check;
check = 0;
for (i = 0; i < address.Nrow(); i++) {
for (j = 0; j < address.Ncol(i); j++) {
if (address[i][j] == &pre) {
check++;
address[i][j] = &post;
}
}
}
if (check != 1)
handle.logMessage(LOGFAIL, "Error in keeper - failed to change variables");
}
void Keeper::clearLast() {
stack->clearString();
}
void Keeper::clearAll() {
stack->clearStack();
}
void Keeper::setString(const char* str) {
stack->clearStack();
stack->storeString(str);
}
void Keeper::addString(const char* str) {
stack->storeString(str);
}
void Keeper::addString(const string str) {
this->addString(str.c_str());
}
int Keeper::numVariables() const {
return switches.Size();
}
void Keeper::getCurrentValues(DoubleVector& val) const {
int i;
for (i = 0; i < values.Size(); i++)
val[i] = values[i];
}
void Keeper::getInitialValues(DoubleVector& val) const {
int i;
for (i = 0; i < initialvalues.Size(); i++)
val[i] = initialvalues[i];
}
void Keeper::getScaledValues(DoubleVector& val) const {
int i;
for (i = 0; i < scaledvalues.Size(); i++)
val[i] = scaledvalues[i];
}
void Keeper::getOptScaledValues(DoubleVector& val) const {
int i, j = 0;
if (val.Size() != numoptvar)
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
for (i = 0; i < scaledvalues.Size(); i++) {
if (opt[i]) {
val[j] = scaledvalues[i];
j++;
}
}
}
void Keeper::getOptInitialValues(DoubleVector& val) const {
int i, j = 0;
if (val.Size() != numoptvar)
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
for (i = 0; i < initialvalues.Size(); i++) {
if (opt[i]) {
val[j] = initialvalues[i];
j++;
}
}
}
void Keeper::resetVariables() {
int i;
for (i = 0; i < values.Size(); i++) {
initialvalues[i] = 1.0;
scaledvalues[i] = values[i];
}
}
void Keeper::scaleVariables() {
int i;
for (i = 0; i < values.Size(); i++) {
if (isZero(values[i])) {
if (opt[i])
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[i].getName());
initialvalues[i] = 1.0;
scaledvalues[i] = values[i];
//JMB - do we want to scale the bounds here as well??
} else {
initialvalues[i] = values[i];
scaledvalues[i] = 1.0;
}
}
}
void Keeper::Update(const DoubleVector& val) {
int i, j;
if (val.Size() != values.Size())
handle.logMessage(LOGFAIL, "Error in keeper - received wrong number of variables to update");
for (i = 0; i < address.Nrow(); i++) {
for (j = 0; j < address.Ncol(i); j++)
*address[i][j].addr = val[i];
values[i] = val[i];
if (isZero(initialvalues[i])) {
if (opt[i])
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[i].getName());
scaledvalues[i] = val[i];
} else
scaledvalues[i] = val[i] / initialvalues[i];
}
}
void Keeper::Update(int pos, double& value) {
int i;
if (pos <= 0 && pos >= address.Nrow())
handle.logMessage(LOGFAIL, "Error in keeper - received invalid variable to update");
for (i = 0; i < address.Ncol(pos); i++)
*address[pos][i].addr = value;
values[pos] = value;
if (isZero(initialvalues[pos])) {
if (opt[pos])
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[pos].getName());
scaledvalues[pos] = value;
} else
scaledvalues[pos] = value / initialvalues[pos];
}
void Keeper::writeBestValues() {
int i, j = 0;
DoubleVector tmpvec(numoptvar, 0.0);
for (i = 0; i < bestvalues.Size(); i++) {
if (opt[i]) {
tmpvec[j] = bestvalues[i];
j++;
}
}
handle.logMessage(LOGINFO, tmpvec);
}
void Keeper::openPrintFile(const char* const filename) {
if (fileopen)
handle.logMessage(LOGFAIL, "Error in keeper - cannot open output file");
fileopen = 1;
outfile.open(filename, ios::out);
handle.checkIfFailure(outfile, filename);
handle.Open(filename);
outfile << "; ";
RUNID.Print(outfile);
}
void Keeper::writeInitialInformation(const LikelihoodPtrVector& likevec) {
if (!fileopen)
handle.logMessage(LOGFAIL, "Error in keeper - cannot write to output file");
int i, j;
outfile << "; Listing of the switches used in the current Gadget run\n";
for (i = 0; i < address.Nrow(); i++) {
outfile << switches[i].getName() << TAB;
for (j = 0; j < address.Ncol(i); j++)
outfile << address[i][j].getName() << TAB;
outfile << endl;
}
outfile << ";\n; Listing of the likelihood components used in the current Gadget run\n;\n";
outfile << "; Component\tType\tWeight\n";
for (i = 0; i < likevec.Size(); i++)
outfile << likevec[i]->getName() << TAB << likevec[i]->getType() << TAB << likevec[i]->getWeight() << endl;
outfile << ";\n; Listing of the output from the likelihood components for the current Gadget run\n;\n";
}
void Keeper::writeValues(const LikelihoodPtrVector& likevec, int prec) {
if (!fileopen)
handle.logMessage(LOGFAIL, "Error in keeper - cannot write to output file");
//JMB - print the number of function evaluations at the start of the line
int iters = EcoSystem->getFuncEval();
#ifdef _OPENMP
if(EcoSystems != NULL) {
int numThr = omp_get_max_threads ( );
for (int i = 0; i < numThr; i++)
iters += EcoSystems[i]->getFuncEval();
}
#endif
outfile << iters << TAB;
int i, p, w;
p = prec;
if (prec == 0)
p = printprecision;
w = p + 4;
for (i = 0; i < values.Size(); i++)
outfile << setw(w) << setprecision(p) << values[i] << sep;
if (prec == 0)
p = smallprecision;
w = p + 4;
outfile << TAB << TAB;
for (i = 0; i < likevec.Size(); i++)
outfile << setw(w) << setprecision(p) << likevec[i]->getUnweightedLikelihood() << sep;
if (prec == 0)
p = fullprecision;
w = p + 4;
outfile << TAB << TAB << setw(w) << setprecision(p) << EcoSystem->getLikelihood() << endl;
}
void Keeper::Update(const StochasticData* const Stoch) {
int i, j;
if (Stoch->numSwitches() > 0) {
if (Stoch->isOptGiven())
boundsgiven = 1;
else
numoptvar = switches.Size();
IntVector match(Stoch->numVariables(), 0);
IntVector found(switches.Size(), 0);
for (i = 0; i < Stoch->numVariables(); i++) {
for (j = 0; j < switches.Size(); j++) {
if (Stoch->getSwitch(i) == switches[j]) {
values[j] = Stoch->getValue(i);
bestvalues[j] = Stoch->getValue(i);
if (!boundsgiven) {
//JMB we are going to optimise all variables
opt[j] = 1;
} else {
lowerbds[j] = Stoch->getLowerBound(i);
upperbds[j] = Stoch->getUpperBound(i);
opt[j] = Stoch->getOptFlag(i);
if (opt[j])
numoptvar++;
}
if (isZero(initialvalues[j])) {
if (opt[j])
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[j].getName());
scaledvalues[j] = values[j];
} else
scaledvalues[j] = values[j] / initialvalues[j];
match[i]++;
found[j]++;
}
}
}
if (handle.getLogLevel() >= LOGWARN) {
for (i = 0; i < Stoch->numVariables(); i++)
if (!match[i])
handle.logMessage(LOGWARN, "Warning in keeper - failed to match switch", Stoch->getSwitch(i).getName());
for (i = 0; i < switches.Size(); i++)
if (!found[i])
handle.logMessage(LOGWARN, "Warning in keeper - using default values for switch", switches[i].getName());
}
} else {
if (this->numVariables() != Stoch->numVariables())
handle.logMessage(LOGFAIL, "Error in keeper - received wrong number of variables to update");
numoptvar = Stoch->numVariables();
for (i = 0; i < numoptvar; i++) {
opt[i] = 1; //JMB we are going to optimise all variables
values[i] = Stoch->getValue(i);
bestvalues[i] = values[i];
if (isZero(initialvalues[i])) {
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[i].getName());
scaledvalues[i] = values[i];
} else
scaledvalues[i] = values[i] / initialvalues[i];
}
}
for (i = 0; i < address.Nrow(); i++)
for (j = 0; j < address.Ncol(i); j++)
*address[i][j].addr = values[i];
}
void Keeper::getOptFlags(IntVector& optimise) const {
int i;
for (i = 0; i < optimise.Size(); i++)
optimise[i] = opt[i];
}
void Keeper::getSwitches(ParameterVector& sw) const {
int i;
for (i = 0; i < sw.Size(); i++)
sw[i] = switches[i];
}
void Keeper::writeParams(const OptInfoPtrVector& optvec, const char* const filename, int prec, int interrupt) {
int i, p, w, check;
ofstream paramfile;
paramfile.open(filename, ios::out);
handle.checkIfFailure(paramfile, filename);
handle.Open(filename);
p = prec;
if (prec == 0)
p = largeprecision;
w = p + 4;
paramfile << "; ";
RUNID.Print(paramfile);
if (interrupt) {
int iters = EcoSystem->getFuncEval();
#ifdef _OPENMP
int numThr = omp_get_max_threads ( );
for (i = 0; i < numThr; i++)
iters += EcoSystems[i]->getFuncEval();
#endif
paramfile << "; Gadget was interrupted after a total of " << iters
<< " function evaluations\n; the best likelihood value found so far is "
<< setprecision(p) << bestlikelihood << endl;
} else if (EcoSystem->getFuncEval() == 0) {
paramfile << "; a simulation run was performed giving a likelihood value of "
<< setprecision(p) << EcoSystem->getLikelihood() << endl;
} else {
for (i = 0; i < optvec.Size(); i++)
optvec[i]->Print(paramfile, p);
}
paramfile << "switch\tvalue\t\tlower\tupper\toptimise\n";
for (i = 0; i < bestvalues.Size(); i++) {
//JMB - if a switch is outside the bounds, we need to reset this back to the bound
//note that the simulation should have used the value of the bound anyway ...
check = 0;
if (lowerbds[i] > bestvalues[i]) {
check++;
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << lowerbds[i];
handle.logMessage(LOGWARN, "Warning in keeper - parameter has a final value", bestvalues[i]);
handle.logMessage(LOGWARN, "which is lower than the corresponding lower bound", lowerbds[i]);
} else if (upperbds[i] < bestvalues[i]) {
check++;
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << upperbds[i];
handle.logMessage(LOGWARN, "Warning in keeper - parameter has a final value", bestvalues[i]);
handle.logMessage(LOGWARN, "which is higher than the corresponding upper bound", upperbds[i]);
} else
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << bestvalues[i];
paramfile << TAB << setw(smallwidth) << setprecision(smallprecision) << lowerbds[i]
<< sep << setw(smallwidth) << setprecision(smallprecision) << upperbds[i]
<< sep << setw(smallwidth) << opt[i];
if (check)
paramfile << " ; warning - parameter has been reset to bound";
paramfile << endl;
}
handle.Close();
paramfile.close();
paramfile.clear();
}
void Keeper::getLowerBounds(DoubleVector& lbs) const {
int i;
for (i = 0; i < lbs.Size(); i++)
lbs[i] = lowerbds[i];
}
void Keeper::getUpperBounds(DoubleVector& ubs) const {
int i;
for (i = 0; i < ubs.Size(); i++)
ubs[i] = upperbds[i];
}
void Keeper::getOptLowerBounds(DoubleVector& lbs) const {
int i, j = 0;
if (lbs.Size() != numoptvar)
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
for (i = 0; i < lowerbds.Size(); i++) {
if (opt[i]) {
lbs[j] = lowerbds[i];
j++;
}
}
}
void Keeper::getOptUpperBounds(DoubleVector& ubs) const {
int i, j = 0;
if (ubs.Size() != numoptvar)
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
for (i = 0; i < upperbds.Size(); i++) {
if (opt[i]) {
ubs[j] = upperbds[i];
j++;
}
}
}
void Keeper::checkBounds(const LikelihoodPtrVector& likevec) const {
if (!boundsgiven)
return;
int i, count;
//check that we have a boundlikelihood component
count = 0;
for (i = 0; i < likevec.Size(); i++)
if (likevec[i]->getType() == BOUNDLIKELIHOOD)
count++;
if ((count == 0) && (values.Size() != 0))
handle.logMessage(LOGWARN, "Warning in keeper - no boundlikelihood component found\nNo penalties will be applied if any of the parameter bounds are exceeded");
if (count > 1)
handle.logMessage(LOGWARN, "Warning in keeper - repeated boundlikelihood components found");
//check the values of the switches are within the bounds to start with
for (i = 0; i < values.Size(); i++) {
if ((lowerbds[i] > values[i]) || (upperbds[i] < values[i]))
handle.logMessage(LOGFAIL, "Error in keeper - initial value outside bounds for parameter", switches[i].getName());
if (upperbds[i] < lowerbds[i])
handle.logMessage(LOGFAIL, "Error in keeper - upper bound lower than lower bound for parameter", switches[i].getName());
if ((lowerbds[i] < 0.0) && (upperbds[i] > 0.0) && (opt[i]))
handle.logMessage(LOGWARN, "Warning in keeper - bounds span zero for parameter", switches[i].getName());
}
}
void Keeper::storeVariables(double likvalue, const DoubleVector& point) {
int i, j = 0;
bestlikelihood = likvalue;
for (i = 0; i < bestvalues.Size(); i++) {
if (opt[i]) {
bestvalues[i] = point[j];
j++;
}
}
}