Annotation of /trunk/paramin-beta/paraminsimann.cc

Parent Directory | Revision Log

---

Revision 1 - (view) (download)

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 :			}

---

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