cloudy trunk
Loading...
Searching...
No Matches
optimize_phymir.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*optimize_phymir Peter van Hoof's optimization routine */
4
5#include "cddefines.h"
6#include "version.h"
7#include "optimize.h"
8#if defined(__unix) || defined(__APPLE__)
9#include <cstddef>
10#include <sys/types.h>
11#include <sys/wait.h>
12#include <unistd.h>
13#else
14#define pid_t int
15#define fork() TotalInsanityAsStub<pid_t>()
16#define wait(X) TotalInsanityAsStub<pid_t>()
17#endif
18#include "mpi_utilities.h"
19
20
23const char* STATEFILE = "continue.pmr";
24const char* STATEFILE_BACKUP = "continue.bak";
25
26/* The optimization algorithm Phymir and its subsidiary routines have been
27 * written by Peter van Hoof. */
28
29inline void wr_block(const void*,size_t,const char*);
30inline void rd_block(void*,size_t,const char*);
31
32
33template<class X, class Y, int NP, int NSTR>
35{
36 DEBUG_ENTRY( "p_clear1()" );
37
38 // first zero out the entire struct so that even the padding gets initialized
39 // this prevents complaints about writing uninitialized bytes by valgrind
40 memset( this, 0, sizeof(*this) );
41
42 p_xmax = numeric_limits<X>::max();
43 p_ymax = numeric_limits<Y>::max() / Y(2.);
44 for( int i=0; i < 2*NP+1; ++i )
45 {
46 for( int j=0; j < NP; ++j )
47 p_xp[i][j] = -p_xmax;
48 p_yp[i] = -p_ymax;
49 }
50 for( int i=0; i < NP; ++i )
51 {
52 p_absmin[i] = -p_xmax;
53 p_absmax[i] = p_xmax;
54 p_varmin[i] = p_xmax;
55 p_varmax[i] = -p_xmax;
56 for( int j=0; j < NP; ++j )
57 p_a2[i][j] = -p_xmax;
58 p_c1[i] = -p_xmax;
59 p_c2[i] = -p_xmax;
60 p_xc[i] = -p_xmax;
61 p_xcold[i] = -p_xmax;
62 }
63 p_vers = VRSNEW;
64 p_toler = -p_xmax;
65 p_dmax = X(0.);
66 p_dold = X(0.);
67 p_ymin = p_ymax;
68 p_dim = int32(NP);
69 p_sdim = int32(NSTR);
70 p_nvar = int32(0);
71 p_noptim = int32(0);
72 p_maxiter = int32(0);
73 p_jmin = int32(-1);
74 p_maxcpu = int32(1);
75 p_curcpu = int32(0);
77 // p_chState, and p_chStr[123] have already been initialized with memset above
78
79 // calculate the size of the struct from the start upto, but not including p_size
80 // this section of the struct will be written in binary mode to a state file
81 // cast pointers to size_t so that we get the size in bytes for fwrite / fread
82 // make p_size an uint32 so that the width of the field is platform independent
83 p_size = static_cast<uint32>(reinterpret_cast<size_t>(&p_size) - reinterpret_cast<size_t>(this));
84 p_func = NULL;
85}
86
87template<class X, class Y, int NP, int NSTR>
88void phymir_state<X,Y,NP,NSTR>::p_wr_state(const char *fnam) const
89{
90 DEBUG_ENTRY( "p_wr_state()" );
91
92 if( cpu.i().lgMaster() && strlen(fnam) > 0 )
93 {
94 FILE *fdes = open_data( fnam, "wb", AS_LOCAL_ONLY_TRY );
95 bool lgErr = ( fdes == NULL );
96 lgErr = lgErr || ( fwrite( &p_size, sizeof(uint32), 1, fdes ) != 1 );
97 lgErr = lgErr || ( fwrite( this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
98 lgErr = lgErr || ( fclose(fdes) != 0 );
99 if( lgErr )
100 {
101 printf( "p_wr_state: error writing file: %s\n", fnam );
102 remove(fnam);
103 }
104 }
105 return;
106}
107
108template<class X, class Y, int NP, int NSTR>
110{
111 DEBUG_ENTRY( "p_rd_state()" );
112
113 FILE *fdes = open_data( fnam, "rb", AS_LOCAL_ONLY );
114 uint32 wrsize;
115 bool lgErr = ( fread( &wrsize, sizeof(uint32), 1, fdes ) != 1 );
116 lgErr = lgErr || ( p_size != wrsize );
117 lgErr = lgErr || ( fread( this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
118 lgErr = lgErr || ( fclose(fdes) != 0 );
119 if( lgErr )
120 {
121 printf( "p_rd_state: error reading file: %s\n", fnam );
123 }
124 return;
125}
126
127template<class X, class Y, int NP, int NSTR>
129 int jj,
130 int runNr)
131{
132 DEBUG_ENTRY( "p_execute_job()" );
133
134 pid_t pid;
135 switch( p_mode )
136 {
137 case PHYMIR_SEQ:
138 return p_lgLimitExceeded(x) ? p_ymax : p_func(x,runNr);
139 case PHYMIR_FORK:
140 p_curcpu++;
141 if( p_curcpu > p_maxcpu )
142 {
143 /* too many processes, wait for one to finish */
144 (void)wait(NULL);
145 p_curcpu--;
146 }
147 /* flush all open files */
148 fflush(NULL);
149 pid = fork();
150 if( pid < 0 )
151 {
152 fprintf( ioQQQ, "creating the child process failed\n" );
154 }
155 else if( pid == 0 )
156 {
157 /* this is child process so execute */
158 p_execute_job_parallel( x, jj, runNr );
159 /* at this point ioQQQ points to the main output of the parent process,
160 * the child should not close that in cdEXIT(), so wipe the handle here */
161 ioQQQ = NULL;
163 }
164 // the real y-value is not available yet...
165 return p_ymax;
166 case PHYMIR_MPI:
167 if( (jj%cpu.i().nCPU()) == cpu.i().nRANK() )
168 p_execute_job_parallel( x, jj, runNr );
169 // the real y-value is not available yet...
170 return p_ymax;
171 default:
173 }
174}
175
176// p_execute_job_parallel MUST be const, otherwise changes by child processes may be lost!
177template<class X, class Y, int NP, int NSTR>
179 int jj,
180 int runNr) const
181{
182 DEBUG_ENTRY( "p_execute_job_parallel()" );
183
184 char fnam1[20], fnam2[20];
185 sprintf(fnam1,"yval_%d",jj);
186 sprintf(fnam2,"output_%d",jj);
187 /* each child should have its own output file */
188 FILE *ioQQQ_old = ioQQQ;
189 ioQQQ = open_data( fnam2, "w", AS_LOCAL_ONLY );
190 // fail-safe in case p_func crashes without being caught...
191 Y yval = p_ymax;
192 wr_block(&yval,sizeof(yval),fnam1);
193 if( !p_lgLimitExceeded(x) )
194 {
195 try
196 {
197 yval = p_func(x,runNr);
198 }
199 catch( ... )
200 {
201 // make sure output is complete since files may not have been closed properly...
202 fflush(NULL);
203 yval = p_ymax;
204 }
205 wr_block(&yval,sizeof(yval),fnam1);
206 }
207 fclose( ioQQQ );
208 ioQQQ = ioQQQ_old;
209}
210
211template<class X, class Y, int NP, int NSTR>
213 int jhi)
214{
215 DEBUG_ENTRY( "p_barrier()" );
216
217 switch( p_mode )
218 {
219 case PHYMIR_SEQ:
220 // nothing to do...
221 break;
222 case PHYMIR_FORK:
223 /* wait for child processes to terminate */
224 while( p_curcpu > 0 )
225 {
226 (void)wait(NULL);
227 p_curcpu--;
228 }
229 p_process_output( jlo, jhi );
230 break;
231 case PHYMIR_MPI:
232 // wait for all function evaluations to finish so that output is complete
233 MPI::COMM_WORLD.Barrier();
234 p_process_output( jlo, jhi );
235 // y values are known now, so broadcast to other ranks
236 MPI::COMM_WORLD.Bcast( &p_yp[jlo], jhi-jlo+1, MPI::type(p_yp), 0 );
237 break;
238 default:
240 }
241}
242
243template<class X, class Y, int NP, int NSTR>
245 int jhi)
246{
247 DEBUG_ENTRY( "p_process_output()" );
248
249 if( cpu.i().lgMaster() )
250 {
251 char fnam[20];
252 for( int jj=jlo; jj <= jhi; jj++ )
253 {
254 sprintf(fnam,"yval_%d",jj);
255 rd_block(&p_yp[jj],sizeof(p_yp[jj]),fnam);
256 remove(fnam);
257 sprintf(fnam,"output_%d",jj);
258 append_file(ioQQQ,fnam);
259 remove(fnam);
260 }
261 }
262}
263
264template<class X, class Y, int NP, int NSTR>
266{
267 DEBUG_ENTRY( "p_evaluate_hyperblock()" );
268
269 int jlo = 1, jhi = 0;
270 for( int j=0; j < p_nvar; j++ )
271 {
272 X sgn = X(1.);
273 for( int jj=2*j+1; jj <= 2*j+2; jj++ )
274 {
275 sgn = -sgn;
276 for( int i=0; i < p_nvar; i++ )
277 {
278 p_xp[jj][i] = p_xc[i] + sgn*p_dmax*p_c2[j]*p_a2[j][i];
279 p_varmax[i] = max(p_varmax[i],p_xp[jj][i]);
280 p_varmin[i] = min(p_varmin[i],p_xp[jj][i]);
281 }
282 if( !lgMaxIterExceeded() )
283 {
284 // p_execute_job will not return the correct chi^2 in parallel mode
285 // only after p_barrier() has finished will p_yp[jj] contain the correct value
286 p_yp[jj] = p_execute_job( p_xp[jj], jj, p_noptim++ );
287 jhi = jj;
288 }
289 }
290 }
291 /* wait for jobs to terminate */
292 p_barrier( jlo, jhi );
293}
294
295template<class X, class Y, int NP, int NSTR>
297{
298 DEBUG_ENTRY( "p_setup_next_hyperblock()" );
299
300 const Y F0 = Y(1.4142136);
301 const X F1 = X(0.7071068);
302 const X F2 = X(0.1);
303
304 /* find best model */
305 for( int jj=1; jj <= 2*p_nvar; jj++ )
306 {
307 if( p_yp[jj] < p_ymin )
308 {
309 p_ymin = p_yp[jj];
310 p_jmin = jj;
311 }
312 }
313 bool lgNewCnt = p_jmin > 0;
314
315 /* determine minimum and relative uncertainties */
316 bool lgNegd2 = false;
317 X xnrm = X(0.);
318 X xhlp[NP];
319 for( int i=0; i < p_nvar; i++ )
320 {
321 Y d1 = p_yp[2*i+2] - p_yp[2*i+1];
322 Y d2 = Y(0.5)*p_yp[2*i+2] - p_yp[0] + Y(0.5)*p_yp[2*i+1];
323 if( d2 <= Y(0.) )
324 lgNegd2 = true;
325 xhlp[i] = -p_dmax*p_c2[i]*(static_cast<X>(max(min((Y(0.25)*d1)/max(d2,Y(1.e-10)),F0),-F0)) -
326 p_delta(2*i+1,p_jmin) + p_delta(2*i+2,p_jmin));
327 xnrm += pow2(xhlp[i]);
328 }
329 xnrm = static_cast<X>(sqrt(xnrm));
330 /* set up new transformation matrix */
331 int imax = 0;
332 X amax = X(0.);
333 for( int j=0; j < p_nvar; j++ )
334 {
335 for( int i=0; i < p_nvar; i++ )
336 {
337 if( xnrm > X(0.) )
338 {
339 if( j == 0 )
340 {
341 p_a2[0][i] *= xhlp[0];
342 }
343 else
344 {
345 p_a2[0][i] += xhlp[j]*p_a2[j][i];
346 p_a2[j][i] = p_delta(i,j);
347 if( j == p_nvar-1 && abs(p_a2[0][i]) > amax )
348 {
349 imax = i;
350 amax = abs(p_a2[0][i]);
351 }
352 }
353 }
354 else
355 {
356 p_a2[j][i] = p_delta(i,j);
357 }
358 }
359 }
360 /* this is to assure maximum linear independence of the base vectors */
361 if( imax > 0 )
362 {
363 p_a2[imax][0] = X(1.);
364 p_a2[imax][imax] = X(0.);
365 }
366 /* apply Gram-Schmidt procedure to get orthogonal base */
367 p_phygrm( p_a2, p_nvar );
368
369 /* set up hyperblock dimensions in new base and choose new center */
370 for( int i=0; i < p_nvar; i++ )
371 {
372 p_c2[i] = X(0.);
373 for( int j=0; j < p_nvar; j++ )
374 {
375 p_c2[i] += pow2(p_a2[i][j]/p_c1[j]);
376 }
377 p_c2[i] = static_cast<X>(1./sqrt(p_c2[i]));
378 p_xc[i] = p_xp[p_jmin][i];
379 p_xp[0][i] = p_xc[i];
380 }
381 p_yp[0] = p_yp[p_jmin];
382 p_jmin = 0;
383 /* choose size of next hyperblock */
384 X r1, r2;
385 if( lgNegd2 )
386 {
387 /* this means that the hyperblock either is bigger than the scale
388 * on which the function varies, or is so small that we see noise.
389 * in both cases make the hyperblock smaller. */
390 r1 = F1;
391 r2 = F1;
392 }
393 else
394 {
395 r1 = F2;
396 if( lgNewCnt )
397 {
398 /* when we make progress, p_dmax may become bigger */
399 r2 = sqrt(1.f/F1);
400 }
401 else
402 {
403 /* when there is no progress force p_dmax smaller, otherwise there
404 * may never be an end */
405 r2 = sqrt(F1);
406 }
407 }
408 p_dmax = min(max(xnrm/p_c2[0],p_dmax*r1),p_dmax*r2);
409 /* fail-safe against excessive behaviour */
411}
412
413template<class X, class Y, int NP, int NSTR>
415{
416 DEBUG_ENTRY( "p_reset_hyperblock()" );
417
418 if( !lgConvergedRestart() )
419 {
420 /* reset hyperblock so that we can restart the optimization */
421 for( int i=0; i < p_nvar; i++ )
422 {
423 p_xcold[i] = p_xc[i];
424 p_c2[i] = p_c1[i];
425 }
426 p_dmax = p_dold;
428 }
429}
430
431template<class X, class Y, int NP, int NSTR>
433 int n)
434{
435 DEBUG_ENTRY( "p_phygrm()" );
436
437 /* apply modified Gram-Schmidt procedure to a */
438 for( int i=0; i < n; i++ )
439 {
440 X ip = X(0.);
441 for( int k=0; k < n; k++ )
442 ip += pow2(a[i][k]);
443 ip = sqrt(ip);
444 for( int k=0; k < n; k++ )
445 a[i][k] /= ip;
446 for( int j=i+1; j < n; j++ )
447 {
448 X ip = X(0.);
449 for( int k=0; k < n; k++ )
450 ip += a[i][k]*a[j][k];
451 for( int k=0; k < n; k++ )
452 a[j][k] -= ip*a[i][k];
453 }
454 }
455 return;
456}
457
458template<class X, class Y, int NP, int NSTR>
460{
461 DEBUG_ENTRY( "p_lgLimitExceeded()" );
462
463 for( int i=0; i < p_nvar; i++ )
464 {
465 if( x[i] < p_absmin[i] )
466 return true;
467 if( x[i] > p_absmax[i] )
468 return true;
469 }
470 return false;
471}
472
473template<class X, class Y, int NP, int NSTR>
475{
476 DEBUG_ENTRY( "p_reset_transformation_matrix()" );
477
478 /* initialize transformation matrix to unity */
479 for( int i=0; i < p_nvar; i++ )
480 for( int j=0; j < p_nvar; j++ )
481 p_a2[j][i] = p_delta(i,j);
482}
483
484template<class X, class Y, int NP, int NSTR>
485void phymir_state<X,Y,NP,NSTR>::init_minmax(const X pmin[], // pmin[nvar]: minimum values for params
486 const X pmax[], // pmax[nvar]: maximum values for params
487 int nvar) // will not be initialized yet
488{
489 DEBUG_ENTRY( "init_minmax()" );
490
491 ASSERT( !lgInitialized() );
492
493 for( int i=0; i < nvar; i++ )
494 {
495 p_absmin[i] = pmin[i];
496 p_absmax[i] = pmax[i];
497 }
498}
499
500template<class X, class Y, int NP, int NSTR>
501void phymir_state<X,Y,NP,NSTR>::init_strings(const char* date, // date string for inclusion in state file
502 const char* version, // version string for inclusion in state file
503 const char* host_name) // host name for inclusion in state file
504{
505 DEBUG_ENTRY( "init_strings()" );
506
507 // use NSTR-1 so that last 0 byte is not overwritten...
508 if( date != NULL )
509 strncpy( p_chStr1, date, NSTR-1 );
510 if( version != NULL )
511 strncpy( p_chStr2, version, NSTR-1 );
512 if( host_name != NULL )
513 strncpy( p_chStr3, host_name, NSTR-1 );
514}
515
516template<class X, class Y, int NP, int NSTR>
517void phymir_state<X,Y,NP,NSTR>::init_state_file_name(const char* fnam) // name of the state file that will be written
518{
519 DEBUG_ENTRY( "init_state_file_name()" );
520
521 // use NSTR-1 so that last 0 byte is not overwritten...
522 strncpy( p_chState, fnam, NSTR-1 );
523}
524
525template<class X, class Y, int NP, int NSTR>
526void phymir_state<X,Y,NP,NSTR>::initial_run(Y (*func)(const X[],int), // function to be optimized
527 int nvar, // number of parameters to be optimized
528 const X start[], // start[n]: initial estimates for params
529 const X del[], // del[n]: initial stepsizes for params
530 X toler, // absolute tolerance on parameters
531 int maxiter, // maximum number of iterations
532 phymir_mode mode, // execution mode: sequential, fork, mpi
533 int maxcpu) // maximum no. of CPUs to use
534{
535 DEBUG_ENTRY( "initial_run()" );
536
537 ASSERT( nvar > 0 && nvar <= NP );
538
539 p_func = func;
540 p_nvar = nvar;
541 p_toler = toler;
542 p_maxiter = maxiter;
543 p_mode = mode;
544 p_maxcpu = maxcpu;
545 p_noptim = 0;
546
547 /* initialize hyperblock dimensions and center */
548 p_dmax = X(0.);
549 for( int i=0; i < p_nvar; i++ )
550 p_dmax = max(p_dmax,abs(del[i]));
551
552 p_dold = p_dmax;
553 for( int i=0; i < p_nvar; i++ )
554 {
555 p_xc[i] = start[i];
556 p_xcold[i] = p_xc[i] + X(10.)*p_toler;
557 p_c1[i] = abs(del[i])/p_dmax;
558 p_c2[i] = p_c1[i];
559 p_xp[0][i] = p_xc[i];
560 p_varmax[i] = max(p_varmax[i],p_xc[i]);
561 p_varmin[i] = min(p_varmin[i],p_xc[i]);
562 }
563
564 // p_execute_job will not return the correct chi^2 in parallel mode
565 // only after p_barrier() has finished will p_yp[0] contain the correct value
566 p_yp[0] = p_execute_job( p_xc, 0, p_noptim++ );
567 p_barrier( 0, 0 );
568
569 p_ymin = p_yp[0];
570 p_jmin = 0;
571
573
575}
576
577template<class X, class Y, int NP, int NSTR>
578void phymir_state<X,Y,NP,NSTR>::continue_from_state(Y (*func)(const X[],int), // function to be optimized
579 int nvar, // number of parameters to be optimized
580 const char* fnam, // file name of the state file
581 X toler, // absolute tolerance on parameters
582 int maxiter, // maximum number of iterations
583 phymir_mode mode, // execution mode: sequential, fork, mpi
584 int maxcpu) // maximum no. of CPUs to use
585{
586 DEBUG_ENTRY( "continue_from_state()" );
587
588 p_rd_state( fnam );
589 // sanity checks
590 if( !fp_equal( p_vers, VRSNEW ) )
591 {
592 printf( "optimize continue - file has incompatible version, sorry\n" );
594 }
595 if( p_dim != NP )
596 {
597 printf( "optimize continue - arrays have wrong dimension, sorry\n" );
599 }
600 if( p_sdim != NSTR )
601 {
602 printf( "optimize continue - strings have wrong length, sorry\n" );
604 }
605 if( p_nvar != nvar )
606 {
607 printf( "optimize continue - wrong number of free parameters, sorry\n" );
609 }
610 // this pointer was not part of the state file, it would be useless...
611 p_func = func;
612 // Option to override the tolerance set in the state file. This is useful
613 // if you want to refine an already finished run to a smaller tolerance.
614 p_toler = toler;
615 // you ran out of iterations, but wanted to continue...
616 p_maxiter = maxiter;
617 // it is OK to continue in a different mode
618 p_mode = mode;
619 p_maxcpu = maxcpu;
620}
621
622template<class X, class Y, int NP, int NSTR>
624{
625 DEBUG_ENTRY( "optimize()" );
626
628
629 while( !lgConverged() )
630 {
632 if( lgMaxIterExceeded() )
633 break;
636 }
637}
638
639template<class X, class Y, int NP, int NSTR>
641{
642 DEBUG_ENTRY( "optimize_with_restart()" );
643
645
646 while( !lgConvergedRestart() )
647 {
648 optimize();
649 if( lgMaxIterExceeded() )
650 break;
652 }
653}
654
655template<class X, class Y, int NP, int NSTR>
657{
658 DEBUG_ENTRY( "lgConvergedRestart()" );
659
660 if( lgConverged() )
661 {
662 X dist = X(0.);
663 for( int i=0; i < p_nvar; i++ )
664 dist += pow2(p_xc[i] - p_xcold[i]);
665 dist = static_cast<X>(sqrt(dist));
666 return ( dist <= p_toler );
667 }
668 return false;
669}
670
672 const realnum del[],
673 long int nvarPhymir,
674 chi2_type *ymin,
675 realnum toler)
676{
677 DEBUG_ENTRY( "optimize_phymir()" );
678
679 if( nvarPhymir > LIMPAR )
680 {
681 fprintf( ioQQQ, "optimize_phymir: too many parameters are varied, increase LIMPAR\n" );
683 }
684
686
687 // first check if a statefile already exists, back it up in that case
688 (void)remove(STATEFILE_BACKUP);
689 FILE *tmp = open_data( STATEFILE, "r", AS_LOCAL_ONLY_TRY );
690 if( tmp != NULL )
691 {
692 fclose( tmp );
693 // create backup copy of statefile
694 FILE *dest = open_data( STATEFILE_BACKUP, "wb", AS_LOCAL_ONLY_TRY );
695 if( dest != NULL )
696 {
697 append_file( dest, STATEFILE );
698 fclose( dest );
699 }
700 }
701
702 phymir_mode mode = optimize.lgParallel ? ( cpu.i().lgMPI() ? PHYMIR_MPI : PHYMIR_FORK ) : PHYMIR_SEQ;
703 long nCPU = optimize.lgParallel ? ( cpu.i().lgMPI() ? cpu.i().nCPU() : optimize.useCPU ) : 1;
704 if( optimize.lgOptCont )
705 {
706 phymir.continue_from_state( optimize_func, nvarPhymir, STATEFILE, toler,
707 optimize.nIterOptim, mode, nCPU );
708 }
709 else
710 {
712 phymir.init_strings( t_version::Inst().chDate, t_version::Inst().chVersion,
713 cpu.i().host_name() );
714 phymir.initial_run( optimize_func, nvarPhymir, xc, del, toler,
715 optimize.nIterOptim, mode, nCPU );
716 }
717
718 phymir.optimize_with_restart();
719
720 if( phymir.lgMaxIterExceeded() )
721 {
722 fprintf( ioQQQ, " Optimizer exceeding maximum iterations.\n" );
723 fprintf( ioQQQ, " This can be reset with the OPTIMIZE ITERATIONS command.\n" );
724 }
725
726 // updates to optimize.nOptimiz and optimize.varmin/max in child processes are lost, so repair here...
727 optimize.nOptimiz = phymir.noptim();
728 for( int i=0; i < nvarPhymir; i++ )
729 {
730 xc[i] = phymir.xval(i);
731 optimize.varmax[i] = min(phymir.xmax(i),optimize.varang[i][1]);
732 optimize.varmin[i] = max(phymir.xmin(i),optimize.varang[i][0]);
733 }
734 *ymin = phymir.yval();
735 return;
736}
737
739inline void wr_block(const void *ptr,
740 size_t len,
741 const char *fnam)
742{
743 DEBUG_ENTRY( "wr_block()" );
744
745 FILE *fdes = open_data( fnam, "wb", AS_LOCAL_ONLY );
746 if( fwrite(ptr,len,size_t(1),fdes) != 1 ) {
747 printf( "error writing on file: %s\n",fnam );
748 fclose(fdes);
750 }
751 fclose(fdes);
752 return;
753}
754
756inline void rd_block(void *ptr,
757 size_t len,
758 const char *fnam)
759{
760 DEBUG_ENTRY( "rd_block()" );
761
762 FILE *fdes = open_data( fnam, "rb", AS_LOCAL_ONLY );
763 if( fread(ptr,len,size_t(1),fdes) != 1 ) {
764 printf( "error reading on file: %s\n",fnam );
765 fclose(fdes);
767 }
768 fclose(fdes);
769 return;
770}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
T pow2(T a)
Definition cddefines.h:931
#define MIN2
Definition cddefines.h:761
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
long max(int a, long b)
Definition cddefines.h:775
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
static t_version & Inst()
Definition cddefines.h:175
bool lgConverged() const
Definition optimize.h:139
void continue_from_state(Y(*)(const X[], int), int, const char *, X, int, phymir_mode, int)
int32 noptim() const
Definition optimize.h:145
X p_absmin[NP]
Definition optimize.h:77
void init_strings(const char *, const char *, const char *)
char p_chStr2[NSTR]
Definition optimize.h:102
bool lgInitialized() const
Definition optimize.h:138
void p_rd_state(const char *)
X p_c2[NP]
Definition optimize.h:83
bool p_lgLimitExceeded(const X[]) const
void p_reset_transformation_matrix()
void p_barrier(int, int)
X p_absmax[NP]
Definition optimize.h:78
X xval(int i) const
Definition optimize.h:141
char p_chStr1[NSTR]
Definition optimize.h:101
void p_phygrm(X[][NP], int)
int32 p_maxcpu
Definition optimize.h:97
Y yval() const
Definition optimize.h:144
int32 p_maxiter
Definition optimize.h:95
int32 p_curcpu
Definition optimize.h:98
void initial_run(Y(*)(const X[], int), int, const X[], const X[], X, int, phymir_mode, int)
Y p_execute_job(const X[], int, int)
void p_evaluate_hyperblock()
int32 p_dim
Definition optimize.h:91
void p_setup_next_hyperblock()
int32 p_nvar
Definition optimize.h:93
void init_state_file_name(const char *)
X p_xp[2 *NP+1][NP]
Definition optimize.h:75
X p_delta(int i, int j) const
Definition optimize.h:123
X xmax(int i) const
Definition optimize.h:143
Y p_yp[2 *NP+1]
Definition optimize.h:76
void p_process_output(int, int)
void optimize_with_restart()
char p_chState[NSTR]
Definition optimize.h:100
int32 p_noptim
Definition optimize.h:94
char p_chStr3[NSTR]
Definition optimize.h:103
Y(* p_func)(const X[], int)
Definition optimize.h:107
int32 p_jmin
Definition optimize.h:96
bool lgConvergedRestart() const
X p_xcold[NP]
Definition optimize.h:85
X p_varmax[NP]
Definition optimize.h:80
int32 p_sdim
Definition optimize.h:92
void init_minmax(const X[], const X[], int)
X p_varmin[NP]
Definition optimize.h:79
bool lgMaxIterExceeded() const
Definition optimize.h:137
void p_execute_job_parallel(const X[], int, int) const
uint32 p_size
Definition optimize.h:106
phymir_mode p_mode
Definition optimize.h:99
X xmin(int i) const
Definition optimize.h:142
X p_c1[NP]
Definition optimize.h:82
void p_wr_state(const char *) const
X p_a2[NP][NP]
Definition optimize.h:81
X p_xc[NP]
Definition optimize.h:84
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
static t_cpu cpu
Definition cpu.h:355
@ AS_LOCAL_ONLY
Definition cpu.h:208
@ AS_LOCAL_ONLY_TRY
Definition cpu.h:207
void append_file(FILE *dest, const char *source)
t_MPI COMM_WORLD
t_optimize optimize
Definition optimize.cpp:5
double chi2_type
Definition optimize.h:49
const realnum VRSNEW
Definition optimize.h:47
chi2_type optimize_func(const realnum param[], int grid_index=-1)
phymir_mode
Definition optimize.h:65
@ PHYMIR_MPI
Definition optimize.h:65
@ PHYMIR_FORK
Definition optimize.h:65
@ PHYMIR_ILL
Definition optimize.h:65
@ PHYMIR_SEQ
Definition optimize.h:65
const long LIMPAR
Definition optimize.h:61
void optimize_phymir(realnum xc[], const realnum del[], long int nvarPhymir, chi2_type *ymin, realnum toler)
#define pid_t
const char * STATEFILE
void rd_block(void *, size_t, const char *)
#define fork()
void wr_block(const void *, size_t, const char *)
#define wait(X)
const char * STATEFILE_BACKUP
STATIC double dist(long, realnum[], realnum[])
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
#define F1(x, y, z)
#define F2(x, y, z)