1 : |
agomez |
1 |
#include "paraminsimann.h" |
2 : |
|
|
|
3 : |
|
|
// ********************************************************
|
4 : |
|
|
// functions for class ParaminSimann
|
5 : |
|
|
// ********************************************************
|
6 : |
|
|
ParaminSimann::ParaminSimann(NetInterface* netInt) : ParaminSearch(netInt) {
|
7 : |
|
|
type = OPTSIMANN;
|
8 : |
|
|
maxim = 0;
|
9 : |
|
|
T = 100.0;
|
10 : |
|
|
cs = 2.0;
|
11 : |
|
|
|
12 : |
|
|
// Vector tempVec(numvar);
|
13 : |
|
|
// vm = tempVec;
|
14 : |
|
|
vm.resize(numvar, 1.0);
|
15 : |
|
|
// vm.setValue(1.0);
|
16 : |
|
|
// xp = tempVec;
|
17 : |
|
|
xp.resize(numvar, 0.0);
|
18 : |
|
|
// changed to iters...
|
19 : |
|
|
// nfcnev = 0;
|
20 : |
|
|
nt = 2;
|
21 : |
|
|
ns = 5;
|
22 : |
|
|
check = 4;
|
23 : |
|
|
uratio = 0.7;
|
24 : |
|
|
lratio = 0.3;
|
25 : |
|
|
rt = 0.85;
|
26 : |
|
|
eps = 1e-4;
|
27 : |
|
|
maxiterations = 2000;
|
28 : |
|
|
/// NEED TO check if 0 is OK or maybe -1????
|
29 : |
|
|
// ID = new int[numvar];
|
30 : |
|
|
// nacp = new int[numvar];
|
31 : |
|
|
// acpPointID = new int[numvar];
|
32 : |
|
|
ID.resize(numvar, 0);
|
33 : |
|
|
nacp.resize(numvar, 0);
|
34 : |
|
|
acpPointID.resize(numvar, 0);
|
35 : |
|
|
// total number of processes initiated at beginning
|
36 : |
|
|
NumberOfHosts = net->getTotalNumProc();
|
37 : |
|
|
// converged = 0;
|
38 : |
|
|
}
|
39 : |
|
|
|
40 : |
|
|
ParaminSimann::~ParaminSimann() {
|
41 : |
|
|
// delete[] ID;
|
42 : |
|
|
// delete[] nacp;
|
43 : |
|
|
// delete[] acpPointID;
|
44 : |
|
|
}
|
45 : |
|
|
|
46 : |
|
|
void ParaminSimann::read(CommentStream& infile, char* text) {
|
47 : |
|
|
int i = 0;
|
48 : |
|
|
int j = 0;
|
49 : |
|
|
double temp;
|
50 : |
|
|
|
51 : |
|
|
while (!infile.eof() && strcasecmp(text, "seed") && strcasecmp(text, "[hooke]") && strcasecmp(text, "[bfgs]")) {
|
52 : |
|
|
infile >> ws;
|
53 : |
|
|
if (strcasecmp(text, "ns") == 0) {
|
54 : |
|
|
infile >> ns;
|
55 : |
|
|
|
56 : |
|
|
} else if (strcasecmp(text, "nt") == 0) {
|
57 : |
|
|
infile >> nt;
|
58 : |
|
|
|
59 : |
|
|
} else if (strcasecmp(text, "check") == 0) {
|
60 : |
|
|
infile >> check;
|
61 : |
|
|
|
62 : |
|
|
} else if (strcasecmp(text, "rt") == 0) {
|
63 : |
|
|
infile >> rt;
|
64 : |
|
|
|
65 : |
|
|
} else if ((strcasecmp(text, "eps") == 0) || (strcasecmp(text, "simanneps") == 0)) {
|
66 : |
|
|
infile >> eps;
|
67 : |
|
|
|
68 : |
|
|
} else if (strcasecmp(text, "t") == 0) {
|
69 : |
|
|
infile >> T;
|
70 : |
|
|
|
71 : |
|
|
} else if (strcasecmp(text, "cstep") == 0) {
|
72 : |
|
|
infile >> cs;
|
73 : |
|
|
|
74 : |
|
|
} else if (strcasecmp(text, "uratio") == 0) {
|
75 : |
|
|
infile >> uratio;
|
76 : |
|
|
|
77 : |
|
|
} else if (strcasecmp(text, "lratio") == 0) {
|
78 : |
|
|
infile >> lratio;
|
79 : |
|
|
|
80 : |
|
|
} else if (strcasecmp(text, "vm") == 0) {
|
81 : |
|
|
infile >> temp;
|
82 : |
|
|
for (j = 0; j < vm.Size(); j++)
|
83 : |
|
|
vm[j] = temp;
|
84 : |
|
|
|
85 : |
|
|
} else if ((strcasecmp(text, "maxiterations") == 0) || (strcasecmp(text, "simanniter") == 0)) {
|
86 : |
|
|
infile >> maxiterations;
|
87 : |
|
|
|
88 : |
|
|
} else {
|
89 : |
|
|
cerr << "Error while reading optinfo for Simulated Annealing - unknown option " << text << endl;
|
90 : |
|
|
exit(EXIT_FAILURE);
|
91 : |
|
|
}
|
92 : |
|
|
infile >> text;
|
93 : |
|
|
i++;
|
94 : |
|
|
}
|
95 : |
|
|
|
96 : |
|
|
if (i == 0)
|
97 : |
|
|
cerr << "Warning - no optinfo give for Simulated Annealing in file - using default parameters" << endl;
|
98 : |
|
|
|
99 : |
|
|
//check the values specified in the optinfo file ...
|
100 : |
|
|
if ((uratio < 0.5) || (uratio > 1)) {
|
101 : |
|
|
cerr << "\nError in value of uratio - setting to default value of 0.7\n";
|
102 : |
|
|
uratio = 0.7;
|
103 : |
|
|
}
|
104 : |
|
|
if ((lratio < 0) || (lratio > 0.5)) {
|
105 : |
|
|
cerr << "\nError in value of lratio - setting to default value of 0.3\n";
|
106 : |
|
|
lratio = 0.3;
|
107 : |
|
|
}
|
108 : |
|
|
if ((rt < 0) || (rt > 1)) {
|
109 : |
|
|
cerr << "\nError in value of rt - setting to default value of 0.85\n";
|
110 : |
|
|
rt = 0.85;
|
111 : |
|
|
}
|
112 : |
|
|
cs = cs / lratio;
|
113 : |
|
|
}
|
114 : |
|
|
|
115 : |
|
|
void ParaminSimann::OptimiseLikelihood() {
|
116 : |
|
|
int i, numtoset;
|
117 : |
|
|
int numloops_ns, numloops_nt;
|
118 : |
|
|
int numset_nsloop; // 0 < numset_nsloop <= numvar
|
119 : |
|
|
int rock = 0;
|
120 : |
|
|
|
121 : |
|
|
bestx = net->getInitialX();
|
122 : |
|
|
xstart = bestx;
|
123 : |
|
|
for (i = 0; i < numvar; i++) {
|
124 : |
|
|
if (xstart[i] > upperbound[i] || xstart[i] < lowerbound[i]) {
|
125 : |
|
|
cerr << "Error in Simmulated Annealing - x is not within bounds\n";
|
126 : |
|
|
exit(EXIT_FAILURE);
|
127 : |
|
|
}
|
128 : |
|
|
}
|
129 : |
|
|
initialVM = vm;
|
130 : |
|
|
bestf = net->getScore();
|
131 : |
|
|
|
132 : |
|
|
if (!maxim)
|
133 : |
|
|
bestf = -bestf;
|
134 : |
|
|
|
135 : |
|
|
// Vector tempVec(check);
|
136 : |
|
|
// tempVec.setValue(bestf);
|
137 : |
|
|
fstar.Reset();
|
138 : |
|
|
fstar.resize(check, bestf);
|
139 : |
|
|
// fstar = tempVec;
|
140 : |
|
|
fstart = bestf;
|
141 : |
|
|
|
142 : |
|
|
for (i = 0; i < numvar; i++) {
|
143 : |
|
|
ID[i] = i;
|
144 : |
|
|
nacp[i] = 0;
|
145 : |
|
|
acpPointID[i] = -1;
|
146 : |
|
|
}
|
147 : |
|
|
// Find out how many values to set at beginning of each ns loop.
|
148 : |
|
|
if (numvar > NumberOfHosts)
|
149 : |
|
|
numtoset = NumberOfHosts;
|
150 : |
|
|
else
|
151 : |
|
|
numtoset = numvar;
|
152 : |
|
|
|
153 : |
|
|
// start the main loop. Note that it terminates if (i) the algorithm
|
154 : |
|
|
// succesfully otimizes the function or (ii) there are too many func. eval.
|
155 : |
|
|
while ((rock == 0) && (iters < maxiterations)) {
|
156 : |
|
|
numloops_nt = 0;
|
157 : |
|
|
net->startNewDataGroup((ns * nt * (numvar + 1)));
|
158 : |
|
|
while ((numloops_nt < nt) && (iters < maxiterations)) {
|
159 : |
|
|
numloops_ns = 0;
|
160 : |
|
|
naccepted_nsloop = 0;
|
161 : |
|
|
randomOrder(ID);
|
162 : |
|
|
|
163 : |
|
|
while ((numloops_ns < ns) && (iters < maxiterations)) {
|
164 : |
|
|
if (numloops_ns == 0) {
|
165 : |
|
|
for (i = 0; i < numtoset; i++)
|
166 : |
|
|
SetXP(ID[i]);
|
167 : |
|
|
numset_nsloop = numtoset;
|
168 : |
|
|
} else
|
169 : |
|
|
numset_nsloop = 0;
|
170 : |
|
|
|
171 : |
|
|
// sends all available data to all free hosts.
|
172 : |
|
|
net->sendToAllIdleHosts();
|
173 : |
|
|
while ((numset_nsloop < numvar) && (iters < maxiterations)) {
|
174 : |
|
|
// get a function value and do any update that's necessary.
|
175 : |
|
|
ReceiveValue();
|
176 : |
|
|
if (iters < maxiterations) {
|
177 : |
|
|
SetXP(ID[numset_nsloop]);
|
178 : |
|
|
numset_nsloop++;
|
179 : |
|
|
}
|
180 : |
|
|
}
|
181 : |
|
|
|
182 : |
|
|
numloops_ns++;
|
183 : |
|
|
}
|
184 : |
|
|
|
185 : |
|
|
if (iters < maxiterations) {
|
186 : |
|
|
while ((net->getNumNotAns() > 0) && (iters < maxiterations)) {
|
187 : |
|
|
// get a function value and do any update that's necessary
|
188 : |
|
|
ReceiveValue();
|
189 : |
|
|
}
|
190 : |
|
|
// adjust vm so approximately half of all evaluations are accepted.
|
191 : |
|
|
UpdateVM();
|
192 : |
|
|
|
193 : |
|
|
} else
|
194 : |
|
|
i = net->receiveAll(); //Why? this takes time
|
195 : |
|
|
|
196 : |
|
|
numloops_nt++;
|
197 : |
|
|
}
|
198 : |
|
|
|
199 : |
|
|
net->stopUsingDataGroup();
|
200 : |
|
|
for (i = 0; i < numvar; i++)
|
201 : |
|
|
acpPointID[i] = -1;
|
202 : |
|
|
|
203 : |
|
|
// check termination criteria.
|
204 : |
|
|
for (i = check - 1; i > 0; i--)
|
205 : |
|
|
fstar[i] = fstar[i - 1];
|
206 : |
|
|
fstar[0] = fstart;
|
207 : |
|
|
|
208 : |
|
|
rock = 0;
|
209 : |
|
|
if (fabs(bestf - fstart) < eps) {
|
210 : |
|
|
rock = 1;
|
211 : |
|
|
for (i = 0; i < check - 1; i++)
|
212 : |
|
|
if (fabs(fstar[i + 1] - fstar[i]) > eps)
|
213 : |
|
|
rock = 0;
|
214 : |
|
|
}
|
215 : |
|
|
|
216 : |
|
|
cout << "\nChecking convergence criteria after " << iters << " function evaluations ...\n";
|
217 : |
|
|
|
218 : |
|
|
// if termination criteria is not met, prepare for another loop.
|
219 : |
|
|
if (rock == 0) {
|
220 : |
|
|
T *= rt;
|
221 : |
|
|
cout << "Reducing the temperature to " << T << endl;
|
222 : |
|
|
bestf = -bestf;
|
223 : |
|
|
cout << "simann best result so far " << bestf << endl;
|
224 : |
|
|
bestf = -bestf;
|
225 : |
|
|
fstart = bestf;
|
226 : |
|
|
xstart = bestx;
|
227 : |
|
|
}
|
228 : |
|
|
}
|
229 : |
|
|
|
230 : |
|
|
// Either reached termination criteria or maxiterations or both
|
231 : |
|
|
if (!maxim)
|
232 : |
|
|
bestf = -bestf;
|
233 : |
|
|
|
234 : |
|
|
net->setInitialScore(bestx, bestf);
|
235 : |
|
|
score = bestf;
|
236 : |
|
|
// Breytti rock == 0 í stað rock == 1
|
237 : |
|
|
if ((iters >= maxiterations) && (rock == 0)) {
|
238 : |
|
|
cout << "\nSimulated Annealing optimisation completed after " << iters
|
239 : |
|
|
<< " iterations (max " << maxiterations << ")\nThe model terminated "
|
240 : |
|
|
<< "because the maximum number of iterations was reached\n";
|
241 : |
|
|
} else {
|
242 : |
|
|
cout << "\nStopping Simulated Annealing\n\nThe optimisation stopped after " << iters
|
243 : |
|
|
<< " function evaluations (max " << maxiterations << ")\nThe optimisation stopped "
|
244 : |
|
|
<< "because an optimum was found for this run\n";
|
245 : |
|
|
converge = 1;
|
246 : |
|
|
}
|
247 : |
|
|
}
|
248 : |
|
|
|
249 : |
|
|
// generate xp, the trial value of x - note use of vm to choose xp.
|
250 : |
|
|
void ParaminSimann::SetXP(int k) {
|
251 : |
|
|
int i;
|
252 : |
|
|
DoubleVector temp;
|
253 : |
|
|
int id;
|
254 : |
|
|
for (i = 0; i < numvar; i++) {
|
255 : |
|
|
if (i == k) {
|
256 : |
|
|
xp[k] = xstart[k] + ((randomNumber() * 2.0) - 1.0) * vm[k];
|
257 : |
|
|
if (xp[k] < lowerbound[k] || xp[k] > upperbound[k]) {
|
258 : |
|
|
while((xp[k] < lowerbound[k] || xp[k] > upperbound[k])) {
|
259 : |
|
|
xp[k] = xstart[k] + ((randomNumber() * 2.0) - 1.0) * vm[k];
|
260 : |
|
|
}
|
261 : |
|
|
|
262 : |
|
|
}
|
263 : |
|
|
} else {
|
264 : |
|
|
if (acpPointID[i] >= 0) {
|
265 : |
|
|
// Use parameter from x with id = nacp_ns[i] which was accepted earlier
|
266 : |
|
|
temp = net->getX(acpPointID[i]);
|
267 : |
|
|
xp[i] = temp[i];
|
268 : |
|
|
} else
|
269 : |
|
|
xp[i] = xstart[i];
|
270 : |
|
|
}
|
271 : |
|
|
}
|
272 : |
|
|
|
273 : |
|
|
if (!net->dataGroupFull()) {
|
274 : |
|
|
// net->setX(xp); // this uses a FIFO stack
|
275 : |
|
|
net->setXFirstToSend(xp); // this uses a LIFO stack
|
276 : |
|
|
} else {
|
277 : |
|
|
cerr << "During Simulated Annealing have set too many values" << endl;
|
278 : |
|
|
exit(EXIT_FAILURE);
|
279 : |
|
|
}
|
280 : |
|
|
net->sendToAllIdleHosts();
|
281 : |
|
|
}
|
282 : |
|
|
|
283 : |
|
|
void ParaminSimann::AcceptPoint() {
|
284 : |
|
|
int i;
|
285 : |
|
|
DoubleVector temp;
|
286 : |
|
|
xstart = xp;
|
287 : |
|
|
fstart = fp;
|
288 : |
|
|
nacp[ID[returnID % numvar]]++;
|
289 : |
|
|
acpPointID[ID[returnID % numvar]] = returnID;
|
290 : |
|
|
naccepted_nsloop++;
|
291 : |
|
|
|
292 : |
|
|
// if better than any other point record as new optimum.
|
293 : |
|
|
if (fp > bestf) {
|
294 : |
|
|
bestx = xp;
|
295 : |
|
|
cout << "\nNew optimum after " << iters << " function evaluations, f(x) = " << -fp << " at\n";
|
296 : |
|
|
temp = net->unscaleX(bestx);
|
297 : |
|
|
for (i = 0; i < temp.Size(); i++)
|
298 : |
|
|
cout << temp[i] << " ";
|
299 : |
|
|
bestf = fp;
|
300 : |
|
|
}
|
301 : |
|
|
}
|
302 : |
|
|
|
303 : |
|
|
void ParaminSimann::UpdateVM() {
|
304 : |
|
|
int i;
|
305 : |
|
|
double ratio;
|
306 : |
|
|
|
307 : |
|
|
// adjust vm so that approximately half of all evaluations are accepted.
|
308 : |
|
|
for (i = 0; i < numvar; i++) {
|
309 : |
|
|
ratio = (double) nacp[i] / ns;
|
310 : |
|
|
nacp[i] = 0;
|
311 : |
|
|
if (ratio > uratio) {
|
312 : |
|
|
vm[i] = vm[i] * (1.0 + cs * (ratio - uratio));
|
313 : |
|
|
} else if (ratio < lratio) {
|
314 : |
|
|
vm[i] = vm[i] / (1.0 + cs * (lratio - ratio));
|
315 : |
|
|
}
|
316 : |
|
|
|
317 : |
|
|
if (vm[i] < rathersmall)
|
318 : |
|
|
vm[i] = rathersmall;
|
319 : |
|
|
if (vm[i] > (upperbound[i] - lowerbound[i]))
|
320 : |
|
|
vm[i] = upperbound[i] - lowerbound[i];
|
321 : |
|
|
|
322 : |
|
|
}
|
323 : |
|
|
}
|
324 : |
|
|
|
325 : |
|
|
// get a function value and do any updates that are necessary.
|
326 : |
|
|
void ParaminSimann::ReceiveValue() {
|
327 : |
|
|
double p, pp;
|
328 : |
|
|
int receive;
|
329 : |
|
|
|
330 : |
|
|
receive = net->receiveOne();
|
331 : |
|
|
if (receive == net->netSuccess()) {
|
332 : |
|
|
returnID = net->getReceiveID();
|
333 : |
|
|
if (returnID >= 0) {
|
334 : |
|
|
// received data belonging to correct datagroup
|
335 : |
|
|
iters++;
|
336 : |
|
|
fp = net->getY(returnID);
|
337 : |
|
|
xp = net->getX(returnID);
|
338 : |
|
|
if (!maxim)
|
339 : |
|
|
fp = -fp;
|
340 : |
|
|
// accept the new point if the function value increases.
|
341 : |
|
|
if (fp >= fstart) {
|
342 : |
|
|
AcceptPoint();
|
343 : |
|
|
//cout << "accepted point" << endl;
|
344 : |
|
|
//cout << "fp is " << fp << endl;
|
345 : |
|
|
}
|
346 : |
|
|
else {
|
347 : |
|
|
p = expRep((fp - fstart) / T);
|
348 : |
|
|
pp = randomNumber();
|
349 : |
|
|
|
350 : |
|
|
// accept the new point if Metropolis condition is satisfied.
|
351 : |
|
|
if (pp < p)
|
352 : |
|
|
AcceptPoint();
|
353 : |
|
|
}
|
354 : |
|
|
}
|
355 : |
|
|
|
356 : |
|
|
} else {
|
357 : |
|
|
cerr << "Trying to receive value during Simulated Annealing but failed" << endl;
|
358 : |
|
|
exit(EXIT_FAILURE);
|
359 : |
|
|
}
|
360 : |
|
|
}
|
361 : |
|
|
void ParaminSimann::Print(ofstream& outfile, int prec) {
|
362 : |
|
|
outfile << "; Simmulated annealing algorithm ran for " << iters
|
363 : |
|
|
<< " function evaluations\n; and stopped when the likelihood value was "
|
364 : |
|
|
<< setprecision(prec) << score;
|
365 : |
|
|
if (converge == -1)
|
366 : |
|
|
outfile << "\n; because an error occured during the optimisation\n";
|
367 : |
|
|
else if (converge == 1)
|
368 : |
|
|
outfile << "\n; because the convergence criteria were met\n";
|
369 : |
|
|
else
|
370 : |
|
|
outfile << "\n; because the maximum number of function evaluations was reached\n";
|
371 : |
|
|
}
|