1 : |
ulcessvp |
19 |
#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 : |
|
|
#ifdef _OPENMP
|
11 : |
|
|
extern Ecosystem** EcoSystems;
|
12 : |
|
|
#endif
|
13 : |
|
|
|
14 : |
|
|
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 : |
|
|
int iters = EcoSystem->getFuncEval();
|
293 : |
|
|
#ifdef _OPENMP
|
294 : |
|
|
if(EcoSystems != NULL) {
|
295 : |
|
|
int numThr = omp_get_max_threads ( );
|
296 : |
|
|
for (int i = 0; i < numThr; i++)
|
297 : |
|
|
iters += EcoSystems[i]->getFuncEval();
|
298 : |
|
|
}
|
299 : |
|
|
#endif
|
300 : |
|
|
outfile << iters << TAB;
|
301 : |
|
|
|
302 : |
|
|
int i, p, w;
|
303 : |
|
|
p = prec;
|
304 : |
|
|
if (prec == 0)
|
305 : |
|
|
p = printprecision;
|
306 : |
|
|
w = p + 4;
|
307 : |
|
|
for (i = 0; i < values.Size(); i++)
|
308 : |
|
|
outfile << setw(w) << setprecision(p) << values[i] << sep;
|
309 : |
|
|
|
310 : |
|
|
if (prec == 0)
|
311 : |
|
|
p = smallprecision;
|
312 : |
|
|
w = p + 4;
|
313 : |
|
|
outfile << TAB << TAB;
|
314 : |
|
|
for (i = 0; i < likevec.Size(); i++)
|
315 : |
|
|
outfile << setw(w) << setprecision(p) << likevec[i]->getUnweightedLikelihood() << sep;
|
316 : |
|
|
|
317 : |
|
|
if (prec == 0)
|
318 : |
|
|
p = fullprecision;
|
319 : |
|
|
w = p + 4;
|
320 : |
|
|
outfile << TAB << TAB << setw(w) << setprecision(p) << EcoSystem->getLikelihood() << endl;
|
321 : |
|
|
}
|
322 : |
|
|
|
323 : |
|
|
void Keeper::Update(const StochasticData* const Stoch) {
|
324 : |
|
|
|
325 : |
|
|
int i, j;
|
326 : |
|
|
if (Stoch->numSwitches() > 0) {
|
327 : |
|
|
if (Stoch->isOptGiven())
|
328 : |
|
|
boundsgiven = 1;
|
329 : |
|
|
else
|
330 : |
|
|
numoptvar = switches.Size();
|
331 : |
|
|
|
332 : |
|
|
IntVector match(Stoch->numVariables(), 0);
|
333 : |
|
|
IntVector found(switches.Size(), 0);
|
334 : |
|
|
for (i = 0; i < Stoch->numVariables(); i++) {
|
335 : |
|
|
for (j = 0; j < switches.Size(); j++) {
|
336 : |
|
|
if (Stoch->getSwitch(i) == switches[j]) {
|
337 : |
|
|
values[j] = Stoch->getValue(i);
|
338 : |
|
|
bestvalues[j] = Stoch->getValue(i);
|
339 : |
|
|
|
340 : |
|
|
if (!boundsgiven) {
|
341 : |
|
|
//JMB we are going to optimise all variables
|
342 : |
|
|
opt[j] = 1;
|
343 : |
|
|
} else {
|
344 : |
|
|
lowerbds[j] = Stoch->getLowerBound(i);
|
345 : |
|
|
upperbds[j] = Stoch->getUpperBound(i);
|
346 : |
|
|
opt[j] = Stoch->getOptFlag(i);
|
347 : |
|
|
if (opt[j])
|
348 : |
|
|
numoptvar++;
|
349 : |
|
|
}
|
350 : |
|
|
|
351 : |
|
|
if (isZero(initialvalues[j])) {
|
352 : |
|
|
if (opt[j])
|
353 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[j].getName());
|
354 : |
|
|
|
355 : |
|
|
scaledvalues[j] = values[j];
|
356 : |
|
|
} else
|
357 : |
|
|
scaledvalues[j] = values[j] / initialvalues[j];
|
358 : |
|
|
|
359 : |
|
|
match[i]++;
|
360 : |
|
|
found[j]++;
|
361 : |
|
|
}
|
362 : |
|
|
}
|
363 : |
|
|
}
|
364 : |
|
|
|
365 : |
|
|
if (handle.getLogLevel() >= LOGWARN) {
|
366 : |
|
|
for (i = 0; i < Stoch->numVariables(); i++)
|
367 : |
|
|
if (!match[i])
|
368 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - failed to match switch", Stoch->getSwitch(i).getName());
|
369 : |
|
|
|
370 : |
|
|
for (i = 0; i < switches.Size(); i++)
|
371 : |
|
|
if (!found[i])
|
372 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - using default values for switch", switches[i].getName());
|
373 : |
|
|
}
|
374 : |
|
|
|
375 : |
|
|
} else {
|
376 : |
|
|
if (this->numVariables() != Stoch->numVariables())
|
377 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received wrong number of variables to update");
|
378 : |
|
|
|
379 : |
|
|
numoptvar = Stoch->numVariables();
|
380 : |
|
|
for (i = 0; i < numoptvar; i++) {
|
381 : |
|
|
opt[i] = 1; //JMB we are going to optimise all variables
|
382 : |
|
|
values[i] = Stoch->getValue(i);
|
383 : |
|
|
bestvalues[i] = values[i];
|
384 : |
|
|
if (isZero(initialvalues[i])) {
|
385 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - cannot scale switch with initial value zero", switches[i].getName());
|
386 : |
|
|
|
387 : |
|
|
scaledvalues[i] = values[i];
|
388 : |
|
|
} else
|
389 : |
|
|
scaledvalues[i] = values[i] / initialvalues[i];
|
390 : |
|
|
}
|
391 : |
|
|
}
|
392 : |
|
|
|
393 : |
|
|
for (i = 0; i < address.Nrow(); i++)
|
394 : |
|
|
for (j = 0; j < address.Ncol(i); j++)
|
395 : |
|
|
*address[i][j].addr = values[i];
|
396 : |
|
|
}
|
397 : |
|
|
|
398 : |
|
|
void Keeper::getOptFlags(IntVector& optimise) const {
|
399 : |
|
|
int i;
|
400 : |
|
|
for (i = 0; i < optimise.Size(); i++)
|
401 : |
|
|
optimise[i] = opt[i];
|
402 : |
|
|
}
|
403 : |
|
|
|
404 : |
|
|
void Keeper::getSwitches(ParameterVector& sw) const {
|
405 : |
|
|
int i;
|
406 : |
|
|
for (i = 0; i < sw.Size(); i++)
|
407 : |
|
|
sw[i] = switches[i];
|
408 : |
|
|
}
|
409 : |
|
|
|
410 : |
|
|
void Keeper::writeParams(const OptInfoPtrVector& optvec, const char* const filename, int prec, int interrupt) {
|
411 : |
|
|
|
412 : |
|
|
int i, p, w, check;
|
413 : |
|
|
ofstream paramfile;
|
414 : |
|
|
paramfile.open(filename, ios::out);
|
415 : |
|
|
handle.checkIfFailure(paramfile, filename);
|
416 : |
|
|
handle.Open(filename);
|
417 : |
|
|
|
418 : |
|
|
p = prec;
|
419 : |
|
|
if (prec == 0)
|
420 : |
|
|
p = largeprecision;
|
421 : |
|
|
w = p + 4;
|
422 : |
|
|
|
423 : |
|
|
paramfile << "; ";
|
424 : |
|
|
RUNID.Print(paramfile);
|
425 : |
|
|
|
426 : |
|
|
if (interrupt) {
|
427 : |
|
|
int iters = EcoSystem->getFuncEval();
|
428 : |
|
|
#ifdef _OPENMP
|
429 : |
|
|
int numThr = omp_get_max_threads ( );
|
430 : |
|
|
for (i = 0; i < numThr; i++)
|
431 : |
|
|
iters += EcoSystems[i]->getFuncEval();
|
432 : |
|
|
#endif
|
433 : |
|
|
paramfile << "; Gadget was interrupted after a total of " << iters
|
434 : |
|
|
<< " function evaluations\n; the best likelihood value found so far is "
|
435 : |
|
|
<< setprecision(p) << bestlikelihood << endl;
|
436 : |
|
|
} else if (EcoSystem->getFuncEval() == 0) {
|
437 : |
|
|
paramfile << "; a simulation run was performed giving a likelihood value of "
|
438 : |
|
|
<< setprecision(p) << EcoSystem->getLikelihood() << endl;
|
439 : |
|
|
|
440 : |
|
|
} else {
|
441 : |
|
|
for (i = 0; i < optvec.Size(); i++)
|
442 : |
|
|
optvec[i]->Print(paramfile, p);
|
443 : |
|
|
}
|
444 : |
|
|
|
445 : |
|
|
paramfile << "switch\tvalue\t\tlower\tupper\toptimise\n";
|
446 : |
|
|
for (i = 0; i < bestvalues.Size(); i++) {
|
447 : |
|
|
//JMB - if a switch is outside the bounds, we need to reset this back to the bound
|
448 : |
|
|
//note that the simulation should have used the value of the bound anyway ...
|
449 : |
|
|
check = 0;
|
450 : |
|
|
if (lowerbds[i] > bestvalues[i]) {
|
451 : |
|
|
check++;
|
452 : |
|
|
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << lowerbds[i];
|
453 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - parameter has a final value", bestvalues[i]);
|
454 : |
|
|
handle.logMessage(LOGWARN, "which is lower than the corresponding lower bound", lowerbds[i]);
|
455 : |
|
|
} else if (upperbds[i] < bestvalues[i]) {
|
456 : |
|
|
check++;
|
457 : |
|
|
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << upperbds[i];
|
458 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - parameter has a final value", bestvalues[i]);
|
459 : |
|
|
handle.logMessage(LOGWARN, "which is higher than the corresponding upper bound", upperbds[i]);
|
460 : |
|
|
} else
|
461 : |
|
|
paramfile << switches[i].getName() << TAB << setw(w) << setprecision(p) << bestvalues[i];
|
462 : |
|
|
|
463 : |
|
|
paramfile << TAB << setw(smallwidth) << setprecision(smallprecision) << lowerbds[i]
|
464 : |
|
|
<< sep << setw(smallwidth) << setprecision(smallprecision) << upperbds[i]
|
465 : |
|
|
<< sep << setw(smallwidth) << opt[i];
|
466 : |
|
|
|
467 : |
|
|
if (check)
|
468 : |
|
|
paramfile << " ; warning - parameter has been reset to bound";
|
469 : |
|
|
paramfile << endl;
|
470 : |
|
|
}
|
471 : |
|
|
handle.Close();
|
472 : |
|
|
paramfile.close();
|
473 : |
|
|
paramfile.clear();
|
474 : |
|
|
}
|
475 : |
|
|
|
476 : |
|
|
void Keeper::getLowerBounds(DoubleVector& lbs) const {
|
477 : |
|
|
int i;
|
478 : |
|
|
for (i = 0; i < lbs.Size(); i++)
|
479 : |
|
|
lbs[i] = lowerbds[i];
|
480 : |
|
|
}
|
481 : |
|
|
|
482 : |
|
|
void Keeper::getUpperBounds(DoubleVector& ubs) const {
|
483 : |
|
|
int i;
|
484 : |
|
|
for (i = 0; i < ubs.Size(); i++)
|
485 : |
|
|
ubs[i] = upperbds[i];
|
486 : |
|
|
}
|
487 : |
|
|
|
488 : |
|
|
void Keeper::getOptLowerBounds(DoubleVector& lbs) const {
|
489 : |
|
|
int i, j = 0;
|
490 : |
|
|
if (lbs.Size() != numoptvar)
|
491 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
|
492 : |
|
|
|
493 : |
|
|
for (i = 0; i < lowerbds.Size(); i++) {
|
494 : |
|
|
if (opt[i]) {
|
495 : |
|
|
lbs[j] = lowerbds[i];
|
496 : |
|
|
j++;
|
497 : |
|
|
}
|
498 : |
|
|
}
|
499 : |
|
|
}
|
500 : |
|
|
|
501 : |
|
|
void Keeper::getOptUpperBounds(DoubleVector& ubs) const {
|
502 : |
|
|
int i, j = 0;
|
503 : |
|
|
if (ubs.Size() != numoptvar)
|
504 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - received invalid number of optimising variables");
|
505 : |
|
|
|
506 : |
|
|
for (i = 0; i < upperbds.Size(); i++) {
|
507 : |
|
|
if (opt[i]) {
|
508 : |
|
|
ubs[j] = upperbds[i];
|
509 : |
|
|
j++;
|
510 : |
|
|
}
|
511 : |
|
|
}
|
512 : |
|
|
}
|
513 : |
|
|
|
514 : |
|
|
void Keeper::checkBounds(const LikelihoodPtrVector& likevec) const {
|
515 : |
|
|
if (!boundsgiven)
|
516 : |
|
|
return;
|
517 : |
|
|
|
518 : |
|
|
int i, count;
|
519 : |
|
|
//check that we have a boundlikelihood component
|
520 : |
|
|
count = 0;
|
521 : |
|
|
for (i = 0; i < likevec.Size(); i++)
|
522 : |
|
|
if (likevec[i]->getType() == BOUNDLIKELIHOOD)
|
523 : |
|
|
count++;
|
524 : |
|
|
|
525 : |
|
|
if ((count == 0) && (values.Size() != 0))
|
526 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - no boundlikelihood component found\nNo penalties will be applied if any of the parameter bounds are exceeded");
|
527 : |
|
|
if (count > 1)
|
528 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - repeated boundlikelihood components found");
|
529 : |
|
|
|
530 : |
|
|
//check the values of the switches are within the bounds to start with
|
531 : |
|
|
for (i = 0; i < values.Size(); i++) {
|
532 : |
|
|
if ((lowerbds[i] > values[i]) || (upperbds[i] < values[i]))
|
533 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - initial value outside bounds for parameter", switches[i].getName());
|
534 : |
|
|
if (upperbds[i] < lowerbds[i])
|
535 : |
|
|
handle.logMessage(LOGFAIL, "Error in keeper - upper bound lower than lower bound for parameter", switches[i].getName());
|
536 : |
|
|
if ((lowerbds[i] < 0.0) && (upperbds[i] > 0.0) && (opt[i]))
|
537 : |
|
|
handle.logMessage(LOGWARN, "Warning in keeper - bounds span zero for parameter", switches[i].getName());
|
538 : |
|
|
}
|
539 : |
|
|
}
|
540 : |
|
|
|
541 : |
|
|
void Keeper::storeVariables(double likvalue, const DoubleVector& point) {
|
542 : |
|
|
int i, j = 0;
|
543 : |
|
|
bestlikelihood = likvalue;
|
544 : |
|
|
for (i = 0; i < bestvalues.Size(); i++) {
|
545 : |
|
|
if (opt[i]) {
|
546 : |
|
|
bestvalues[i] = point[j];
|
547 : |
|
|
j++;
|
548 : |
|
|
}
|
549 : |
|
|
}
|
550 : |
|
|
}
|