| 1 : |
agomez |
1 |
#include "keeper.h" |
| 2 : |
|
|
#include "errorhandler.h" |
| 3 : |
|
|
#include "runid.h" |
| 4 : |
|
|
#include "ecosystem.h" |
| 5 : |
|
|
#include "gadget.h" |
| 6 : |
|
|
#include "global.h" |
| 7 : |
|
|
|
| 8 : |
|
|
extern Ecosystem* EcoSystem;
|
| 9 : |
|
|
|
| 10 : |
ulcessvp |
11 |
#ifndef NO_OPENMP
|
| 11 : |
|
|
extern Ecosystem** EcoSystems;
|
| 12 : |
|
|
#endif
|
| 13 : |
|
|
|
| 14 : |
agomez |
1 |
Keeper::Keeper() {
|
| 15 : |
|
|
stack = new StrStack();
|
| 16 : |
|
|
boundsgiven = 0;
|
| 17 : |
|
|
fileopen = 0;
|
| 18 : |
|
|
numoptvar = 0;
|
| 19 : |
|
|
bestlikelihood = 0.0;
|
| 20 : |
|
|
}
|
| 21 : |
|
|
|
| 22 : |
|
|
void Keeper::keepVariable(double& value, Parameter& attr) {
|
| 23 : |
|
|
|
| 24 : |
|
|
int i, index = -1;
|
| 25 : |
|
|
for (i = 0; i < switches.Size(); i++)
|
| 26 : |
|
|
if (switches[i] == attr)
|
| 27 : |
|
|
index = i;
|
| 28 : |
|
|
|
| 29 : |
|
|
if (index == -1) {
|
| 30 : |
|
|
//attr was not found -- add it to switches and values
|
| 31 : |
|
|
index = switches.Size();
|
| 32 : |
|
|
switches.resize(attr);
|
| 33 : |
|
|
values.resize(1, value);
|
| 34 : |
|
|
bestvalues.resize(1, value);
|
| 35 : |
|
|
lowerbds.resize(1, -9999.0); // default lower bound
|
| 36 : |
|
|
upperbds.resize(1, 9999.0); // default upper bound
|
| 37 : |
|
|
opt.resize(1, 0);
|
| 38 : |
|
|
scaledvalues.resize(1, 1.0);
|
| 39 : |
|
|
initialvalues.resize(1, 1.0);
|
| 40 : |
|
|
address.resize();
|
| 41 : |
|
|
address[index].resize();
|
| 42 : |
|
|
address[index][0] = &value;
|
| 43 : |
|
|
if (stack->getSize() != 0)
|
| 44 : |
|
|
address[index][0] = stack->sendAll();
|
| 45 : |
|
|
|
| 46 : |
|
|
} else {
|
| 47 : |
|
|
if (value != values[index]) {
|
| 48 : |
|
|
handle.logFileMessage(LOGFAIL, "read repeated switch name but different initial value", switches[index].getName());
|
| 49 : |
|
|
|
| 50 : |
|
|
} else {
|
| 51 : |
|
|
i = address[index].Size();
|
| 52 : |
|
|
address[index].resize();
|
| 53 : |
|
|
address[index][i] = &value;
|
| 54 : |
|
|
if (stack->getSize() != 0)
|
| 55 : |
|
|
address[index][i] = stack->sendAll();
|
| 56 : |
|
|
}
|
| 57 : |
|
|
}
|
| 58 : |
|
|
}
|
| 59 : |
|
|
|
| 60 : |
|
|
Keeper::~Keeper() {
|
| 61 : |
|
|
delete stack;
|
| 62 : |
|
|
if (fileopen) {
|
| 63 : |
|
|
handle.Close();
|
| 64 : |
|
|
outfile.close();
|
| 65 : |
|
|
outfile.clear();
|
| 66 : |
|
|
}
|
| 67 : |
|
|
}
|
| 68 : |
|
|
|
| 69 : |
|
|
void Keeper::deleteParameter(const double& var) {
|
| 70 : |
|
|
int i, j, check;
|
| 71 : |
|
|
check = 0;
|
| 72 : |
|
|
for (i = 0; i < address.Nrow(); i++) {
|
| 73 : |
|
|
for (j = 0; j < address[i].Size(); j++) {
|
| 74 : |
|
|
if (address[i][j] == &var) {
|
| 75 : |
|
|
check++;
|
| 76 : |
|
|
address[i].Delete(j);
|
| 77 : |
|
|
if (address[i].Size() == 0) {
|
| 78 : |
|
|
//the variable we deleted was the only one with this switch
|
| 79 : |
|
|
address.Delete(i);
|
| 80 : |
|
|
switches.Delete(i);
|
| 81 : |
|
|
values.Delete(i);
|
| 82 : |
|
|
bestvalues.Delete(i);
|
| 83 : |
|
|
opt.Delete(i);
|
| 84 : |
|
|
lowerbds.Delete(i);
|
| 85 : |
|
|
upperbds.Delete(i);
|
| 86 : |
|
|
scaledvalues.Delete(i);
|
| 87 : |
|
|
initialvalues.Delete(i);
|
| 88 : |
|
|
i--;
|
| 89 : |
|
|
}
|
| 90 : |
|
|
}
|
| 91 : |
|
|
}
|
| 92 : |
|
|
}
|
| 93 : |
|
|
if (check != 1)
|
| 94 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - failed to delete parameter");
|
| 95 : |
|
|
}
|
| 96 : |
|
|
|
| 97 : |
|
|
void Keeper::changeVariable(const double& pre, double& post) {
|
| 98 : |
|
|
int i, j, check;
|
| 99 : |
|
|
check = 0;
|
| 100 : |
|
|
for (i = 0; i < address.Nrow(); i++) {
|
| 101 : |
|
|
for (j = 0; j < address.Ncol(i); j++) {
|
| 102 : |
|
|
if (address[i][j] == &pre) {
|
| 103 : |
|
|
check++;
|
| 104 : |
|
|
address[i][j] = &post;
|
| 105 : |
|
|
}
|
| 106 : |
|
|
}
|
| 107 : |
|
|
}
|
| 108 : |
|
|
if (check != 1)
|
| 109 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - failed to change variables");
|
| 110 : |
|
|
}
|
| 111 : |
|
|
|
| 112 : |
|
|
void Keeper::clearLast() {
|
| 113 : |
|
|
stack->clearString();
|
| 114 : |
|
|
}
|
| 115 : |
|
|
|
| 116 : |
|
|
void Keeper::clearAll() {
|
| 117 : |
|
|
stack->clearStack();
|
| 118 : |
|
|
}
|
| 119 : |
|
|
|
| 120 : |
|
|
void Keeper::setString(const char* str) {
|
| 121 : |
|
|
stack->clearStack();
|
| 122 : |
|
|
stack->storeString(str);
|
| 123 : |
|
|
}
|
| 124 : |
|
|
|
| 125 : |
|
|
void Keeper::addString(const char* str) {
|
| 126 : |
|
|
stack->storeString(str);
|
| 127 : |
|
|
}
|
| 128 : |
|
|
|
| 129 : |
|
|
void Keeper::addString(const string str) {
|
| 130 : |
|
|
this->addString(str.c_str());
|
| 131 : |
|
|
}
|
| 132 : |
|
|
|
| 133 : |
|
|
int Keeper::numVariables() const {
|
| 134 : |
|
|
return switches.Size();
|
| 135 : |
|
|
}
|
| 136 : |
|
|
|
| 137 : |
|
|
void Keeper::getCurrentValues(DoubleVector& val) const {
|
| 138 : |
|
|
int i;
|
| 139 : |
|
|
for (i = 0; i < values.Size(); i++)
|
| 140 : |
|
|
val[i] = values[i];
|
| 141 : |
|
|
}
|
| 142 : |
|
|
|
| 143 : |
|
|
void Keeper::getInitialValues(DoubleVector& val) const {
|
| 144 : |
|
|
int i;
|
| 145 : |
|
|
for (i = 0; i < initialvalues.Size(); i++)
|
| 146 : |
|
|
val[i] = initialvalues[i];
|
| 147 : |
|
|
}
|
| 148 : |
|
|
|
| 149 : |
|
|
void Keeper::getScaledValues(DoubleVector& val) const {
|
| 150 : |
|
|
int i;
|
| 151 : |
|
|
for (i = 0; i < scaledvalues.Size(); i++)
|
| 152 : |
|
|
val[i] = scaledvalues[i];
|
| 153 : |
|
|
}
|
| 154 : |
|
|
|
| 155 : |
|
|
void Keeper::getOptScaledValues(DoubleVector& val) const {
|
| 156 : |
|
|
int i, j = 0;
|
| 157 : |
|
|
if (val.Size() != numoptvar)
|
| 158 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
|
| 159 : |
|
|
|
| 160 : |
|
|
for (i = 0; i < scaledvalues.Size(); i++) {
|
| 161 : |
|
|
if (opt[i]) {
|
| 162 : |
|
|
val[j] = scaledvalues[i];
|
| 163 : |
|
|
j++;
|
| 164 : |
|
|
}
|
| 165 : |
|
|
}
|
| 166 : |
|
|
}
|
| 167 : |
|
|
|
| 168 : |
|
|
void Keeper::getOptInitialValues(DoubleVector& val) const {
|
| 169 : |
|
|
int i, j = 0;
|
| 170 : |
|
|
if (val.Size() != numoptvar)
|
| 171 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
|
| 172 : |
|
|
|
| 173 : |
|
|
for (i = 0; i < initialvalues.Size(); i++) {
|
| 174 : |
|
|
if (opt[i]) {
|
| 175 : |
|
|
val[j] = initialvalues[i];
|
| 176 : |
|
|
j++;
|
| 177 : |
|
|
}
|
| 178 : |
|
|
}
|
| 179 : |
|
|
}
|
| 180 : |
|
|
|
| 181 : |
|
|
void Keeper::resetVariables() {
|
| 182 : |
|
|
int i;
|
| 183 : |
|
|
for (i = 0; i < values.Size(); i++) {
|
| 184 : |
|
|
initialvalues[i] = 1.0;
|
| 185 : |
|
|
scaledvalues[i] = values[i];
|
| 186 : |
|
|
}
|
| 187 : |
|
|
}
|
| 188 : |
|
|
|
| 189 : |
|
|
void Keeper::scaleVariables() {
|
| 190 : |
|
|
int i;
|
| 191 : |
|
|
for (i = 0; i < values.Size(); i++) {
|
| 192 : |
|
|
if (isZero(values[i])) {
|
| 193 : |
|
|
if (opt[i])
|
| 194 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[i].getName());
|
| 195 : |
|
|
|
| 196 : |
|
|
initialvalues[i] = 1.0;
|
| 197 : |
|
|
scaledvalues[i] = values[i];
|
| 198 : |
|
|
//JMB - do we want to scale the bounds here as well??
|
| 199 : |
|
|
} else {
|
| 200 : |
|
|
initialvalues[i] = values[i];
|
| 201 : |
|
|
scaledvalues[i] = 1.0;
|
| 202 : |
|
|
}
|
| 203 : |
|
|
}
|
| 204 : |
|
|
}
|
| 205 : |
|
|
|
| 206 : |
|
|
void Keeper::Update(const DoubleVector& val) {
|
| 207 : |
|
|
int i, j;
|
| 208 : |
|
|
if (val.Size() != values.Size())
|
| 209 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received wrong number of variables to update");
|
| 210 : |
|
|
|
| 211 : |
|
|
for (i = 0; i < address.Nrow(); i++) {
|
| 212 : |
|
|
for (j = 0; j < address.Ncol(i); j++)
|
| 213 : |
|
|
*address[i][j].addr = val[i];
|
| 214 : |
|
|
|
| 215 : |
|
|
values[i] = val[i];
|
| 216 : |
|
|
if (isZero(initialvalues[i])) {
|
| 217 : |
|
|
if (opt[i])
|
| 218 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[i].getName());
|
| 219 : |
|
|
|
| 220 : |
|
|
scaledvalues[i] = val[i];
|
| 221 : |
|
|
} else
|
| 222 : |
|
|
scaledvalues[i] = val[i] / initialvalues[i];
|
| 223 : |
|
|
}
|
| 224 : |
|
|
}
|
| 225 : |
|
|
|
| 226 : |
|
|
void Keeper::Update(int pos, double& value) {
|
| 227 : |
|
|
int i;
|
| 228 : |
|
|
if (pos <= 0 && pos >= address.Nrow())
|
| 229 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received invalid variable to update");
|
| 230 : |
|
|
|
| 231 : |
|
|
for (i = 0; i < address.Ncol(pos); i++)
|
| 232 : |
|
|
*address[pos][i].addr = value;
|
| 233 : |
|
|
|
| 234 : |
|
|
values[pos] = value;
|
| 235 : |
|
|
if (isZero(initialvalues[pos])) {
|
| 236 : |
|
|
if (opt[pos])
|
| 237 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[pos].getName());
|
| 238 : |
|
|
|
| 239 : |
|
|
scaledvalues[pos] = value;
|
| 240 : |
|
|
} else
|
| 241 : |
|
|
scaledvalues[pos] = value / initialvalues[pos];
|
| 242 : |
|
|
}
|
| 243 : |
|
|
|
| 244 : |
|
|
void Keeper::writeBestValues() {
|
| 245 : |
|
|
int i, j = 0;
|
| 246 : |
|
|
DoubleVector tmpvec(numoptvar, 0.0);
|
| 247 : |
|
|
for (i = 0; i < bestvalues.Size(); i++) {
|
| 248 : |
|
|
if (opt[i]) {
|
| 249 : |
|
|
tmpvec[j] = bestvalues[i];
|
| 250 : |
|
|
j++;
|
| 251 : |
|
|
}
|
| 252 : |
|
|
}
|
| 253 : |
|
|
handle.logMessage(LOGINFO, tmpvec);
|
| 254 : |
|
|
}
|
| 255 : |
|
|
|
| 256 : |
|
|
void Keeper::openPrintFile(const char* const filename) {
|
| 257 : |
|
|
if (fileopen)
|
| 258 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - cannot open output file");
|
| 259 : |
|
|
fileopen = 1;
|
| 260 : |
|
|
outfile.open(filename, ios::out);
|
| 261 : |
|
|
handle.checkIfFailure(outfile, filename);
|
| 262 : |
|
|
handle.Open(filename);
|
| 263 : |
|
|
outfile << "; ";
|
| 264 : |
|
|
RUNID.Print(outfile);
|
| 265 : |
|
|
}
|
| 266 : |
|
|
|
| 267 : |
|
|
void Keeper::writeInitialInformation(const LikelihoodPtrVector& likevec) {
|
| 268 : |
|
|
if (!fileopen)
|
| 269 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - cannot write to output file");
|
| 270 : |
|
|
|
| 271 : |
|
|
int i, j;
|
| 272 : |
|
|
outfile << "; Listing of the switches used in the current Gadget run\n";
|
| 273 : |
|
|
for (i = 0; i < address.Nrow(); i++) {
|
| 274 : |
|
|
outfile << switches[i].getName() << TAB;
|
| 275 : |
|
|
for (j = 0; j < address.Ncol(i); j++)
|
| 276 : |
|
|
outfile << address[i][j].getName() << TAB;
|
| 277 : |
|
|
outfile << endl;
|
| 278 : |
|
|
}
|
| 279 : |
|
|
|
| 280 : |
|
|
outfile << ";\n; Listing of the likelihood components used in the current Gadget run\n;\n";
|
| 281 : |
|
|
outfile << "; Component\tType\tWeight\n";
|
| 282 : |
|
|
for (i = 0; i < likevec.Size(); i++)
|
| 283 : |
|
|
outfile << likevec[i]->getName() << TAB << likevec[i]->getType() << TAB << likevec[i]->getWeight() << endl;
|
| 284 : |
|
|
outfile << ";\n; Listing of the output from the likelihood components for the current Gadget run\n;\n";
|
| 285 : |
|
|
}
|
| 286 : |
|
|
|
| 287 : |
|
|
void Keeper::writeValues(const LikelihoodPtrVector& likevec, int prec) {
|
| 288 : |
|
|
if (!fileopen)
|
| 289 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - cannot write to output file");
|
| 290 : |
|
|
|
| 291 : |
|
|
//JMB - print the number of function evaluations at the start of the line
|
| 292 : |
ulcessvp |
11 |
int iters = EcoSystem->getFuncEval();
|
| 293 : |
|
|
#ifndef NO_OPENMP
|
| 294 : |
|
|
int numThr = omp_get_max_threads ( );
|
| 295 : |
|
|
for (int i = 0; i < numThr; i++)
|
| 296 : |
|
|
iters += EcoSystems[i]->getFuncEval();
|
| 297 : |
|
|
#endif
|
| 298 : |
|
|
outfile << iters << TAB;
|
| 299 : |
agomez |
1 |
|
| 300 : |
|
|
int i, p, w;
|
| 301 : |
|
|
p = prec;
|
| 302 : |
|
|
if (prec == 0)
|
| 303 : |
|
|
p = printprecision;
|
| 304 : |
|
|
w = p + 4;
|
| 305 : |
|
|
for (i = 0; i < values.Size(); i++)
|
| 306 : |
|
|
outfile << setw(w) << setprecision(p) << values[i] << sep;
|
| 307 : |
|
|
|
| 308 : |
|
|
if (prec == 0)
|
| 309 : |
|
|
p = smallprecision;
|
| 310 : |
|
|
w = p + 4;
|
| 311 : |
|
|
outfile << TAB << TAB;
|
| 312 : |
|
|
for (i = 0; i < likevec.Size(); i++)
|
| 313 : |
|
|
outfile << setw(w) << setprecision(p) << likevec[i]->getUnweightedLikelihood() << sep;
|
| 314 : |
|
|
|
| 315 : |
|
|
if (prec == 0)
|
| 316 : |
|
|
p = fullprecision;
|
| 317 : |
|
|
w = p + 4;
|
| 318 : |
|
|
outfile << TAB << TAB << setw(w) << setprecision(p) << EcoSystem->getLikelihood() << endl;
|
| 319 : |
|
|
}
|
| 320 : |
|
|
|
| 321 : |
|
|
void Keeper::Update(const StochasticData* const Stoch) {
|
| 322 : |
|
|
|
| 323 : |
|
|
int i, j;
|
| 324 : |
|
|
if (Stoch->numSwitches() > 0) {
|
| 325 : |
|
|
if (Stoch->isOptGiven())
|
| 326 : |
|
|
boundsgiven = 1;
|
| 327 : |
|
|
else
|
| 328 : |
|
|
numoptvar = switches.Size();
|
| 329 : |
|
|
|
| 330 : |
|
|
IntVector match(Stoch->numVariables(), 0);
|
| 331 : |
|
|
IntVector found(switches.Size(), 0);
|
| 332 : |
|
|
for (i = 0; i < Stoch->numVariables(); i++) {
|
| 333 : |
|
|
for (j = 0; j < switches.Size(); j++) {
|
| 334 : |
|
|
if (Stoch->getSwitch(i) == switches[j]) {
|
| 335 : |
|
|
values[j] = Stoch->getValue(i);
|
| 336 : |
|
|
bestvalues[j] = Stoch->getValue(i);
|
| 337 : |
|
|
|
| 338 : |
|
|
if (!boundsgiven) {
|
| 339 : |
|
|
//JMB we are going to optimise all variables
|
| 340 : |
|
|
opt[j] = 1;
|
| 341 : |
|
|
} else {
|
| 342 : |
|
|
lowerbds[j] = Stoch->getLowerBound(i);
|
| 343 : |
|
|
upperbds[j] = Stoch->getUpperBound(i);
|
| 344 : |
|
|
opt[j] = Stoch->getOptFlag(i);
|
| 345 : |
|
|
if (opt[j])
|
| 346 : |
|
|
numoptvar++;
|
| 347 : |
|
|
}
|
| 348 : |
|
|
|
| 349 : |
|
|
if (isZero(initialvalues[j])) {
|
| 350 : |
|
|
if (opt[j])
|
| 351 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[j].getName());
|
| 352 : |
|
|
|
| 353 : |
|
|
scaledvalues[j] = values[j];
|
| 354 : |
|
|
} else
|
| 355 : |
|
|
scaledvalues[j] = values[j] / initialvalues[j];
|
| 356 : |
|
|
|
| 357 : |
|
|
match[i]++;
|
| 358 : |
|
|
found[j]++;
|
| 359 : |
|
|
}
|
| 360 : |
|
|
}
|
| 361 : |
|
|
}
|
| 362 : |
|
|
|
| 363 : |
|
|
if (handle.getLogLevel() >= LOGWARN) {
|
| 364 : |
|
|
for (i = 0; i < Stoch->numVariables(); i++)
|
| 365 : |
|
|
if (!match[i])
|
| 366 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - failed to match switch", Stoch->getSwitch(i).getName());
|
| 367 : |
|
|
|
| 368 : |
|
|
for (i = 0; i < switches.Size(); i++)
|
| 369 : |
|
|
if (!found[i])
|
| 370 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - using default values for switch", switches[i].getName());
|
| 371 : |
|
|
}
|
| 372 : |
|
|
|
| 373 : |
|
|
} else {
|
| 374 : |
|
|
if (this->numVariables() != Stoch->numVariables())
|
| 375 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received wrong number of variables to update");
|
| 376 : |
|
|
|
| 377 : |
|
|
numoptvar = Stoch->numVariables();
|
| 378 : |
|
|
for (i = 0; i < numoptvar; i++) {
|
| 379 : |
|
|
opt[i] = 1; //JMB we are going to optimise all variables
|
| 380 : |
|
|
values[i] = Stoch->getValue(i);
|
| 381 : |
|
|
bestvalues[i] = values[i];
|
| 382 : |
|
|
if (isZero(initialvalues[i])) {
|
| 383 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[i].getName());
|
| 384 : |
|
|
|
| 385 : |
|
|
scaledvalues[i] = values[i];
|
| 386 : |
|
|
} else
|
| 387 : |
|
|
scaledvalues[i] = values[i] / initialvalues[i];
|
| 388 : |
|
|
}
|
| 389 : |
|
|
}
|
| 390 : |
|
|
|
| 391 : |
|
|
for (i = 0; i < address.Nrow(); i++)
|
| 392 : |
|
|
for (j = 0; j < address.Ncol(i); j++)
|
| 393 : |
|
|
*address[i][j].addr = values[i];
|
| 394 : |
|
|
}
|
| 395 : |
|
|
|
| 396 : |
|
|
void Keeper::getOptFlags(IntVector& optimise) const {
|
| 397 : |
|
|
int i;
|
| 398 : |
|
|
for (i = 0; i < optimise.Size(); i++)
|
| 399 : |
|
|
optimise[i] = opt[i];
|
| 400 : |
|
|
}
|
| 401 : |
|
|
|
| 402 : |
|
|
void Keeper::getSwitches(ParameterVector& sw) const {
|
| 403 : |
|
|
int i;
|
| 404 : |
|
|
for (i = 0; i < sw.Size(); i++)
|
| 405 : |
|
|
sw[i] = switches[i];
|
| 406 : |
|
|
}
|
| 407 : |
|
|
|
| 408 : |
|
|
void Keeper::writeParams(const OptInfoPtrVector& optvec, const char* const filename, int prec, int interrupt) {
|
| 409 : |
|
|
|
| 410 : |
|
|
int i, p, w, check;
|
| 411 : |
|
|
ofstream paramfile;
|
| 412 : |
|
|
paramfile.open(filename, ios::out);
|
| 413 : |
|
|
handle.checkIfFailure(paramfile, filename);
|
| 414 : |
|
|
handle.Open(filename);
|
| 415 : |
|
|
|
| 416 : |
|
|
p = prec;
|
| 417 : |
|
|
if (prec == 0)
|
| 418 : |
|
|
p = largeprecision;
|
| 419 : |
|
|
w = p + 4;
|
| 420 : |
|
|
|
| 421 : |
|
|
paramfile << "; ";
|
| 422 : |
|
|
RUNID.Print(paramfile);
|
| 423 : |
|
|
|
| 424 : |
|
|
if (interrupt) {
|
| 425 : |
ulcessvp |
11 |
int iters = EcoSystem->getFuncEval();
|
| 426 : |
|
|
#ifndef NO_OPENMP
|
| 427 : |
|
|
int numThr = omp_get_max_threads ( );
|
| 428 : |
|
|
for (i = 0; i < numThr; i++)
|
| 429 : |
|
|
iters += EcoSystems[i]->getFuncEval();
|
| 430 : |
|
|
#endif
|
| 431 : |
|
|
paramfile << "; Gadget was interrupted after a total of " << iters
|
| 432 : |
agomez |
1 |
<< " function evaluations\n; the best likelihood value found so far is "
|
| 433 : |
|
|
<< setprecision(p) << bestlikelihood << endl;
|
| 434 : |
|
|
|
| 435 : |
|
|
} else if (EcoSystem->getFuncEval() == 0) {
|
| 436 : |
|
|
paramfile << "; a simulation run was performed giving a likelihood value of "
|
| 437 : |
|
|
<< setprecision(p) << EcoSystem->getLikelihood() << endl;
|
| 438 : |
|
|
|
| 439 : |
|
|
} else {
|
| 440 : |
|
|
for (i = 0; i < optvec.Size(); i++)
|
| 441 : |
|
|
optvec[i]->Print(paramfile, p);
|
| 442 : |
|
|
}
|
| 443 : |
|
|
|
| 444 : |
|
|
paramfile << "switch\tvalue\t\tlower\tupper\toptimise\n";
|
| 445 : |
|
|
for (i = 0; i < bestvalues.Size(); i++) {
|
| 446 : |
|
|
//JMB - if a switch is outside the bounds, we need to reset this back to the bound
|
| 447 : |
|
|
//note that the simulation should have used the value of the bound anyway ...
|
| 448 : |
|
|
check = 0;
|
| 449 : |
|
|
if (lowerbds[i] > bestvalues[i]) {
|
| 450 : |
|
|
check++;
|
| 451 : |
|
|
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << lowerbds[i];
|
| 452 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - parameter has a final value", bestvalues[i]);
|
| 453 : |
|
|
handle.logMessage(LOGWARN, "which is lower than the corresponding lower bound", lowerbds[i]);
|
| 454 : |
|
|
} else if (upperbds[i] < bestvalues[i]) {
|
| 455 : |
|
|
check++;
|
| 456 : |
|
|
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << upperbds[i];
|
| 457 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - parameter has a final value", bestvalues[i]);
|
| 458 : |
|
|
handle.logMessage(LOGWARN, "which is higher than the corresponding upper bound", upperbds[i]);
|
| 459 : |
|
|
} else
|
| 460 : |
|
|
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << bestvalues[i];
|
| 461 : |
|
|
|
| 462 : |
|
|
paramfile << TAB << setw(smallwidth) << setprecision(smallprecision) << lowerbds[i]
|
| 463 : |
|
|
<< sep << setw(smallwidth) << setprecision(smallprecision) << upperbds[i]
|
| 464 : |
|
|
<< sep << setw(smallwidth) << opt[i];
|
| 465 : |
|
|
|
| 466 : |
|
|
if (check)
|
| 467 : |
|
|
paramfile << " ; warning - parameter has been reset to bound";
|
| 468 : |
|
|
paramfile << endl;
|
| 469 : |
|
|
}
|
| 470 : |
|
|
handle.Close();
|
| 471 : |
|
|
paramfile.close();
|
| 472 : |
|
|
paramfile.clear();
|
| 473 : |
|
|
}
|
| 474 : |
|
|
|
| 475 : |
|
|
void Keeper::getLowerBounds(DoubleVector& lbs) const {
|
| 476 : |
|
|
int i;
|
| 477 : |
|
|
for (i = 0; i < lbs.Size(); i++)
|
| 478 : |
|
|
lbs[i] = lowerbds[i];
|
| 479 : |
|
|
}
|
| 480 : |
|
|
|
| 481 : |
|
|
void Keeper::getUpperBounds(DoubleVector& ubs) const {
|
| 482 : |
|
|
int i;
|
| 483 : |
|
|
for (i = 0; i < ubs.Size(); i++)
|
| 484 : |
|
|
ubs[i] = upperbds[i];
|
| 485 : |
|
|
}
|
| 486 : |
|
|
|
| 487 : |
|
|
void Keeper::getOptLowerBounds(DoubleVector& lbs) const {
|
| 488 : |
|
|
int i, j = 0;
|
| 489 : |
|
|
if (lbs.Size() != numoptvar)
|
| 490 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
|
| 491 : |
|
|
|
| 492 : |
|
|
for (i = 0; i < lowerbds.Size(); i++) {
|
| 493 : |
|
|
if (opt[i]) {
|
| 494 : |
|
|
lbs[j] = lowerbds[i];
|
| 495 : |
|
|
j++;
|
| 496 : |
|
|
}
|
| 497 : |
|
|
}
|
| 498 : |
|
|
}
|
| 499 : |
|
|
|
| 500 : |
|
|
void Keeper::getOptUpperBounds(DoubleVector& ubs) const {
|
| 501 : |
|
|
int i, j = 0;
|
| 502 : |
|
|
if (ubs.Size() != numoptvar)
|
| 503 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
|
| 504 : |
|
|
|
| 505 : |
|
|
for (i = 0; i < upperbds.Size(); i++) {
|
| 506 : |
|
|
if (opt[i]) {
|
| 507 : |
|
|
ubs[j] = upperbds[i];
|
| 508 : |
|
|
j++;
|
| 509 : |
|
|
}
|
| 510 : |
|
|
}
|
| 511 : |
|
|
}
|
| 512 : |
|
|
|
| 513 : |
|
|
void Keeper::checkBounds(const LikelihoodPtrVector& likevec) const {
|
| 514 : |
|
|
if (!boundsgiven)
|
| 515 : |
|
|
return;
|
| 516 : |
|
|
|
| 517 : |
|
|
int i, count;
|
| 518 : |
|
|
//check that we have a boundlikelihood component
|
| 519 : |
|
|
count = 0;
|
| 520 : |
|
|
for (i = 0; i < likevec.Size(); i++)
|
| 521 : |
|
|
if (likevec[i]->getType() == BOUNDLIKELIHOOD)
|
| 522 : |
|
|
count++;
|
| 523 : |
|
|
|
| 524 : |
|
|
if ((count == 0) && (values.Size() != 0))
|
| 525 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - no boundlikelihood component found\nNo penalties will be applied if any of the parameter bounds are exceeded");
|
| 526 : |
|
|
if (count > 1)
|
| 527 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - repeated boundlikelihood components found");
|
| 528 : |
|
|
|
| 529 : |
|
|
//check the values of the switches are within the bounds to start with
|
| 530 : |
|
|
for (i = 0; i < values.Size(); i++) {
|
| 531 : |
|
|
if ((lowerbds[i] > values[i]) || (upperbds[i] < values[i]))
|
| 532 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - initial value outside bounds for parameter", switches[i].getName());
|
| 533 : |
|
|
if (upperbds[i] < lowerbds[i])
|
| 534 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - upper bound lower than lower bound for parameter", switches[i].getName());
|
| 535 : |
|
|
if ((lowerbds[i] < 0.0) && (upperbds[i] > 0.0) && (opt[i]))
|
| 536 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - bounds span zero for parameter", switches[i].getName());
|
| 537 : |
|
|
}
|
| 538 : |
|
|
}
|
| 539 : |
|
|
|
| 540 : |
|
|
void Keeper::storeVariables(double likvalue, const DoubleVector& point) {
|
| 541 : |
|
|
int i, j = 0;
|
| 542 : |
|
|
bestlikelihood = likvalue;
|
| 543 : |
|
|
for (i = 0; i < bestvalues.Size(); i++) {
|
| 544 : |
|
|
if (opt[i]) {
|
| 545 : |
|
|
bestvalues[i] = point[j];
|
| 546 : |
|
|
j++;
|
| 547 : |
|
|
}
|
| 548 : |
|
|
}
|
| 549 : |
|
|
}
|
Parent Directory
Revision Log