Diff of /trunk/gadget/hooke.cc

Parent Directory | Revision Log | Patch

	revision 11, Thu Jul 23 19:00:38 2015 UTC	revision 14, Thu Jul 30 12:36:46 2015 UTC
#	Line 156	Line 156
156	DoubleVector z(point);	DoubleVector z(point);
157
158	minf = prevbest;	minf = prevbest;
	//  for (int k=0;k<point.Size(); k++)
	//        cout << z[k] ;
	//  cout << endl;
159	for (i = 0; i < point.Size(); i++) {	for (i = 0; i < point.Size(); i++) {
	//        for (int k=0;k<point.Size(); k++)
	//                cout << z[k] << " " ;
	//cout << endl;
160	z[param[i]] = point[param[i]] + delta[param[i]];	z[param[i]] = point[param[i]] + delta[param[i]];
161	ftmp = EcoSystem->SimulateAndUpdate(z);	ftmp = EcoSystem->SimulateAndUpdate(z);
	//    cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " << ftmp << endl;
162	if (ftmp < minf) {	if (ftmp < minf) {
163	minf = ftmp;	minf = ftmp;
164	} else {	} else {
#	Line 178	Line 170
170	else	else
171	z[param[i]] = point[param[i]];	z[param[i]] = point[param[i]];
172	}	}
	//    cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " << ftmp  <<" - " << prevbest << endl;
173	}	}
174
175	for (i = 0; i < point.Size(); i++)	for (i = 0; i < point.Size(); i++)
#	Line 188	Line 179
179
180	/* given a point, look for a better one nearby, one coord at a time */	/* given a point, look for a better one nearby, one coord at a time */
181	#ifndef NO_OPENMP	#ifndef NO_OPENMP
182		/*
183		* function bestBeraby parallelized with OpenMP
184		* · 2 threads per coord to parallelize the calculation of +delta/-delta
185		* · parallelize the calculation of the best nearby of the coord
186		*/
187	double OptInfoHooke::bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {	double OptInfoHooke::bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
	//cout << "beastNearbyINI" << endl;
188	double minf;//, ftmp;	double minf;//, ftmp;
189	int i, j, k;	int i, j, k;
190	DoubleVector z(point);	DoubleVector z(point);
#	Line 215	Line 210
210	}	}
211
212	while ( i < nvars) {	while ( i < nvars) {
	//        for (int k=0;k<point.Size(); k++)
	//                cout << z[k] << " " ;
	//          cout << endl;
213	if ((i + paral_tokens -1) >= nvars)	if ((i + paral_tokens -1) >= nvars)
214	paral_tokens = nvars - i;	paral_tokens = nvars - i;
	//cout << "i:" << i << endl;
215	omp_set_dynamic(0);	omp_set_dynamic(0);
216	omp_set_nested(1);	omp_set_nested(1); //permit the nested parallelization
217	#pragma omp parallel for num_threads(paral_tokens) private(k)	#pragma omp parallel for num_threads(paral_tokens) private(k) //parallelize the parameters (numThr/2)
218	for (j = 0; j < paral_tokens; ++j) {	for (j = 0; j < paral_tokens; ++j) {
	//                cout << "thr_for:" << omp_get_num_threads()<<endl;
	//                cout << "J:" << j << endl;
219	storage[j].z = z;	storage[j].z = z;
220	storage[j].delta = delta;	storage[j].delta = delta;
221	DoubleVector v1(z);	DoubleVector v1(z);
#	Line 235	Line 224
224	v1[k] +=  delta[k];	v1[k] +=  delta[k];
225	v2[k] -=  delta[k];	v2[k] -=  delta[k];
226
227	#pragma omp parallel sections num_threads(2)	#pragma omp parallel sections num_threads(2) //parrallelize the +/- delta simulation for each parameter
228	{	{
	//                        cout << "thr_sec:" << omp_get_num_threads()<<endl;
229	#pragma omp section	#pragma omp section
230	{	{
231	storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);	storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
	//                cout << "->" << j << endl;
232	}	}
233	#pragma omp section	#pragma omp section
234	{	{
	//                                cout << "->" << j+paral_tokens << endl;
235	storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);	storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
236	}	}
237	}	}
	//cout << "-----" << endl;
	//                cout <<i+j<< "-z["<< param[i+j]<<"]:" << (storage[j].z)[param[i+j]] << " - " << storage[j].ftmp << endl;
238
239	if (storage[j].ftmp < minf) {	if (storage[j].ftmp < minf) {
240	storage[j].iters = 1;	storage[j].iters = 1;
241	storage[j].z[k] = v1[k];	storage[j].z[k] = v1[k];
	//                        storage[j].delta[k] += delta[k];
	//                      minf = ftmp;
242	} else {	} else {
243	storage[j].iters = 2;	storage[j].iters = 2;
244	storage[j].delta[k] = 0.0 - delta[k];	storage[j].delta[k] = 0.0 - delta[k];
245	if (storage[j+paral_tokens].ftmp < minf) {	if (storage[j+paral_tokens].ftmp < minf) {
	//                                minf = ftmp;
246	storage[j].ftmp = storage[j+paral_tokens].ftmp;	storage[j].ftmp = storage[j+paral_tokens].ftmp;
247	storage[j].z[k] = v2[k];	storage[j].z[k] = v2[k];
248	}	}
	//                        else
	//                                storage[j].z[k] = point[k];
249	}	}
250	}	}
251
	//        cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " ;
252	for (j = 0; j < paral_tokens; ++j) {	for (j = 0; j < paral_tokens; ++j) {
253	i++;	i++;
254	iters += storage[j].iters;	iters += storage[j].iters;
#	Line 278	Line 256
256	minf = storage[j].ftmp;	minf = storage[j].ftmp;
257	z = storage[j].z;	z = storage[j].z;
258	delta = storage[j].delta;	delta = storage[j].delta;
	//                        cout << storage[j].ftmp;
259	break;	break;
260	}	}
261	}	}
262	}	}
263
	//  omp_set_num_threads(numThr);
264	for (i = 0; i < nvars; ++i)	for (i = 0; i < nvars; ++i)
265	point[i] = z[i];	point[i] = z[i];
266	return minf;	return minf;
#	Line 316	Line 291
291	EcoSystem->scaleVariables();	EcoSystem->scaleVariables();
292	#ifndef NO_OPENMP	#ifndef NO_OPENMP
293	int numThr = omp_get_max_threads ( );	int numThr = omp_get_max_threads ( );
294	for (i = 0; i < numThr; i++)	for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
295	EcoSystems[i]->scaleVariables();	EcoSystems[i]->scaleVariables();
296	#endif	#endif
297	EcoSystem->getOptScaledValues(x);	EcoSystem->getOptScaledValues(x);
#	Line 337	Line 312
312	bestx[i] = x[i];	bestx[i] = x[i];
313	trialx[i] = x[i];	trialx[i] = x[i];
314	param[i] = i;	param[i] = i;
	//FIXME rand_r
315	delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign	delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign
316	}	}
317
#	Line 371	Line 345
345	/* randomize the order of the parameters once in a while */	/* randomize the order of the parameters once in a while */
346	rchange = 0;	rchange = 0;
347	while (rchange < nvars) {	while (rchange < nvars) {
	//FIXME rand_r
348	rnumber = rand() % nvars;	rnumber = rand() % nvars;
349	rcheck = 1;	rcheck = 1;
350	for (i = 0; i < rchange; i++)	for (i = 0; i < rchange; i++)
#	Line 401	Line 374
374	#ifdef NO_OPENMP	#ifdef NO_OPENMP
375	iters = EcoSystem->getFuncEval() - offset;	iters = EcoSystem->getFuncEval() - offset;
376	#endif	#endif
	//    cout << "newf:" << newf << "-" << iters<<endl;
377	if (iters > hookeiter) {	if (iters > hookeiter) {
378	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
379	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
#	Line 499	Line 471
471	#ifdef NO_OPENMP	#ifdef NO_OPENMP
472	iters = EcoSystem->getFuncEval() - offset;	iters = EcoSystem->getFuncEval() - offset;
473	#endif	#endif
	//      cout << "newf:" << newf << "-" << iters<<endl;
474	if (iters > hookeiter) {	if (iters > hookeiter) {
475	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
476	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
#	Line 552	Line 523
523	}	}
524	}	}
525
526		/* Functions to perform the parallelization of the algorithm of HJ with OpenMP*/
527	#ifdef GADGET_OPENMP	#ifdef GADGET_OPENMP
	//TODO doc
528	double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {	double OptInfoHooke::bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param) {
529	//cout << "beastNearbyINI" << endl;	double minf;
	double minf;//, ftmp;
530	int i, j, k, ii;	int i, j, k, ii;
531	DoubleVector z(point);	DoubleVector z(point);
532	int bestId = 0;	int bestId = 0;
#	Line 569	Line 539
539
540	minf = prevbest;	minf = prevbest;
541
542	int paral_tokens, numThr, nvars = point.Size();	int paral_tokens, numThr, nvars = point.Size();
543	numThr = omp_get_max_threads ( );	numThr = omp_get_max_threads ( );
544
#	Line 583	Line 552
552	for (ii=0; ii< paral_tokens; ii++) {	for (ii=0; ii< paral_tokens; ii++) {
553	i = 0;	i = 0;
554	while ( i < nvars) {	while ( i < nvars) {
	//        for (int k=0;k<point.Size(); k++)
	//                cout << z[k] << " " ;
	//          cout << endl;
555	if ((i + paral_tokens -1) >= nvars)	if ((i + paral_tokens -1) >= nvars)
556	paral_tokens = nvars - i;	paral_tokens = nvars - i;
	//cout << "i:" << i << endl;
557	omp_set_dynamic(0);	omp_set_dynamic(0);
558	omp_set_nested(1);	omp_set_nested(1);
559	#pragma omp parallel for num_threads(paral_tokens) private(k)	#pragma omp parallel for num_threads(paral_tokens) private(k)
560	for (j = 0; j < paral_tokens; ++j) {	for (j = 0; j < paral_tokens; ++j) {
	//                cout << "thr_for:" << omp_get_num_threads()<<endl;
	//                cout << "J:" << j << endl;
561	storage[j].z = z;	storage[j].z = z;
562	storage[j].delta = delta;	storage[j].delta = delta;
563	DoubleVector v1(z);	DoubleVector v1(z);
#	Line 605	Line 568
568
569	#pragma omp parallel sections num_threads(2)	#pragma omp parallel sections num_threads(2)
570	{	{
	//                        cout << "thr_sec:" << omp_get_num_threads()<<endl;
571	#pragma omp section	#pragma omp section
572	{	{
573	storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);	storage[j].ftmp = EcoSystems[j]->SimulateAndUpdate(v1);
	//                cout << "->" << j << endl;
574	}	}
575	#pragma omp section	#pragma omp section
576	{	{
	//                                cout << "->" << j+paral_tokens << endl;
577	storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);	storage[j+paral_tokens].ftmp = EcoSystems[j+paral_tokens]->SimulateAndUpdate(v2);
578	}	}
579	}	}
	//cout << "-----" << endl;
	//                cout <<i+j<< "-z["<< param[i+j]<<"]:" << (storage[j].z)[param[i+j]] << " - " << storage[j].ftmp << endl;
580	if (storage[j].ftmp < minf) {	if (storage[j].ftmp < minf) {
581	storage[j].iters = 1;	storage[j].iters = 1;
582	storage[j].z[k] = v1[k];	storage[j].z[k] = v1[k];
	//                        storage[j].delta[k] += delta[k];
	//                      minf = ftmp;
583	} else {	} else {
584	storage[j].iters = 2;	storage[j].iters = 2;
585	storage[j].delta[k] = 0.0 - delta[k];	storage[j].delta[k] = 0.0 - delta[k];
586	if (storage[j+paral_tokens].ftmp < minf) {	if (storage[j+paral_tokens].ftmp < minf) {
	//                                minf = ftmp;
587	storage[j].ftmp = storage[j+paral_tokens].ftmp;	storage[j].ftmp = storage[j+paral_tokens].ftmp;
588	storage[j].z[k] = v2[k];	storage[j].z[k] = v2[k];
589	}	}
590	else iters += 2;	else iters += 2;
	//                                storage[j].z[k] = point[k];
591	}	}
592	}	}
593
	//        cout << i <<"-z["<< param[i]<<"]:" <<z[param[i]] << " - " ;
594	bestId = 0;	bestId = 0;
595	for (j = 1; j < paral_tokens; ++j) {	for (j = 1; j < paral_tokens; ++j) {
596	if (storage[j].ftmp < storage[bestId].ftmp)	if (storage[j].ftmp < storage[bestId].ftmp)
#	Line 649	Line 601
601	minf = storage[bestId].ftmp;	minf = storage[bestId].ftmp;
602	z = storage[bestId].z;	z = storage[bestId].z;
603	delta = storage[bestId].delta;	delta = storage[bestId].delta;
	//                        cout << storage[j].ftmp;
604	}	}
605
606	i += paral_tokens;	i += paral_tokens;
607	}	}
608	}	}
609
	//  omp_set_num_threads(numThr);
610	delete[] storage;	delete[] storage;
611	for (i = 0; i < nvars; ++i)	for (i = 0; i < nvars; ++i)
612	point[i] = z[i];	point[i] = z[i];
#	Line 686	Line 635
635	IntVector trapped(nvars, 0);	IntVector trapped(nvars, 0);
636
637	EcoSystem->scaleVariables();	EcoSystem->scaleVariables();
	#ifndef NO_OPENMP
638	int numThr = omp_get_max_threads ( );	int numThr = omp_get_max_threads ( );
639	for (i = 0; i < numThr; i++)	for (i = 0; i < numThr; i++) // scale the variables for the ecosystem of every thread
640	EcoSystems[i]->scaleVariables();	EcoSystems[i]->scaleVariables();
	#endif
641	EcoSystem->getOptScaledValues(x);	EcoSystem->getOptScaledValues(x);
642	EcoSystem->getOptLowerBounds(lowerb);	EcoSystem->getOptLowerBounds(lowerb);
643	EcoSystem->getOptUpperBounds(upperb);	EcoSystem->getOptUpperBounds(upperb);
#	Line 709	Line 656
656	bestx[i] = x[i];	bestx[i] = x[i];
657	trialx[i] = x[i];	trialx[i] = x[i];
658	param[i] = i;	param[i] = i;
	//FIXME rand_r
659	delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign	delta[i] = ((2 * (rand() % 2)) - 1) * rho;  //JMB - randomise the sign
660	}	}
661
#	Line 743	Line 689
689	/* randomize the order of the parameters once in a while */	/* randomize the order of the parameters once in a while */
690	rchange = 0;	rchange = 0;
691	while (rchange < nvars) {	while (rchange < nvars) {
	//FIXME rand_r
692	rnumber = rand() % nvars;	rnumber = rand() % nvars;
693	rcheck = 1;	rcheck = 1;
694	for (i = 0; i < rchange; i++)	for (i = 0; i < rchange; i++)
#	Line 773	Line 718
718	#ifdef NO_OPENMP	#ifdef NO_OPENMP
719	iters = EcoSystem->getFuncEval() - offset;	iters = EcoSystem->getFuncEval() - offset;
720	#endif	#endif
	//    cout << "newf:" << newf << "-" << iters<<endl;
721	if (iters > hookeiter) {	if (iters > hookeiter) {
722	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
723	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
#	Line 786	Line 730
730	for (i = 0; i < nvars; i++)	for (i = 0; i < nvars; i++)
731	bestx[i] = trialx[i] * init[i];	bestx[i] = trialx[i] * init[i];
732	EcoSystem->storeVariables(score, bestx);	EcoSystem->storeVariables(score, bestx);
	//             EcoSystem->writeOptValuesOMP();
733	return;	return;
734	}	}
735
#	Line 872	Line 815
815	#ifdef NO_OPENMP	#ifdef NO_OPENMP
816	iters = EcoSystem->getFuncEval() - offset;	iters = EcoSystem->getFuncEval() - offset;
817	#endif	#endif
	//      cout << "newf:" << newf << "-" << iters<<endl;
818	if (iters > hookeiter) {	if (iters > hookeiter) {
819	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");	handle.logMessage(LOGINFO, "\nStopping Hooke & Jeeves optimisation algorithm\n");
820	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");	handle.logMessage(LOGINFO, "The optimisation stopped after", iters, "function evaluations");
#	Line 885	Line 827
827	for (i = 0; i < nvars; i++)	for (i = 0; i < nvars; i++)
828	bestx[i] = trialx[i] * init[i];	bestx[i] = trialx[i] * init[i];
829	EcoSystem->storeVariables(score, bestx);	EcoSystem->storeVariables(score, bestx);
	//               EcoSystem->writeOptValuesOMP();
830	return;	return;
831	}	}
832	}	}
#	Line 916	Line 857
857	score = bestf;	score = bestf;
858	handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);	handle.logMessage(LOGINFO, "\nHooke & Jeeves finished with a likelihood score of", score);
859	EcoSystem->storeVariables(bestf, bestx);	EcoSystem->storeVariables(bestf, bestx);
	//             EcoSystem->writeOptValuesOMP();
860	return;	return;
861	}	}
862

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.11   changed lines   Added in v.14

---

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