Log In | Get Help   
Home My Page Projects Code Snippets Project Openings Mareframe
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files
[mareframe] Annotation of /trunk/gadget/keeper.cc
[mareframe] / trunk / gadget / keeper.cc Repository:
ViewVC logotype

Annotation of /trunk/gadget/keeper.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (view) (download)

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

root@forge.cesga.es
ViewVC Help
Powered by ViewVC 1.0.0  

Powered By FusionForge