8 #if defined(__unix) || defined(__APPLE__)
10 #include <sys/types.h>
15 #define fork() TotalInsanityAsStub<pid_t>()
16 #define wait(X) TotalInsanityAsStub<pid_t>()
28 inline void wr_block(
const void*,
size_t,
const char*);
29 inline void rd_block(
void*,
size_t,
const char*);
32 template<
class X,
class Y,
int NP,
int NSTR>
39 memset(
this, 0,
sizeof(*
this) );
43 for(
int i=0; i < 2*NP+1; ++i )
45 for(
int j=0; j < NP; ++j )
49 for(
int i=0; i < NP; ++i )
51 p_absmin[i] = -p_xmax;
54 p_varmax[i] = -p_xmax;
55 for(
int j=0; j < NP; ++j )
82 p_size =
static_cast<uint32
>(
reinterpret_cast<size_t>(&p_size) - reinterpret_cast<size_t>(
this));
86 template<
class X,
class Y,
int NP,
int NSTR>
94 bool lgErr = ( fdes == NULL );
95 lgErr = lgErr || ( fwrite( &p_size,
sizeof(uint32), 1, fdes ) != 1 );
96 lgErr = lgErr || ( fwrite(
this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
97 lgErr = lgErr || ( fclose(fdes) != 0 );
100 printf(
"p_wr_state: error writing file: %s\n", fnam );
107 template<
class X,
class Y,
int NP,
int NSTR>
114 bool lgErr = ( fread( &wrsize,
sizeof(uint32), 1, fdes ) != 1 );
115 lgErr = lgErr || ( p_size != wrsize );
116 lgErr = lgErr || ( fread(
this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
117 lgErr = lgErr || ( fclose(fdes) != 0 );
120 printf(
"p_rd_state: error reading file: %s\n", fnam );
126 template<
class X,
class Y,
int NP,
int NSTR>
137 return p_lgLimitExceeded(x) ? p_ymax : p_func(x,runNr);
140 if( p_curcpu > p_maxcpu )
151 fprintf(
ioQQQ,
"creating the child process failed\n" );
157 p_execute_job_parallel( x, jj, runNr );
167 p_execute_job_parallel( x, jj, runNr );
176 template<
class X,
class Y,
int NP,
int NSTR>
183 char fnam1[20], fnam2[20];
184 sprintf(fnam1,
"yval_%d",jj);
185 sprintf(fnam2,
"output_%d",jj);
187 FILE *ioQQQ_old =
ioQQQ;
192 if( !p_lgLimitExceeded(x) )
196 yval = p_func(x,runNr);
210 template<
class X,
class Y,
int NP,
int NSTR>
223 while( p_curcpu > 0 )
228 p_process_output( jlo, jhi );
233 p_process_output( jlo, jhi );
235 MPI_Bcast( &p_yp[jlo], jhi-jlo+1, MPI_type(p_yp), 0, MPI_COMM_WORLD );
242 template<
class X,
class Y,
int NP,
int NSTR>
251 for(
int jj=jlo; jj <= jhi; jj++ )
253 sprintf(fnam,
"yval_%d",jj);
254 rd_block(&p_yp[jj],
sizeof(p_yp[jj]),fnam);
256 sprintf(fnam,
"output_%d",jj);
263 template<
class X,
class Y,
int NP,
int NSTR>
268 int jlo = 1, jhi = 0;
269 for(
int j=0; j < p_nvar; j++ )
272 for(
int jj=2*j+1; jj <= 2*j+2; jj++ )
275 for(
int i=0; i < p_nvar; i++ )
277 p_xp[jj][i] = p_xc[i] + sgn*p_dmax*p_c2[j]*p_a2[j][i];
278 p_varmax[i] =
max(p_varmax[i],p_xp[jj][i]);
279 p_varmin[i] =
min(p_varmin[i],p_xp[jj][i]);
281 if( !lgMaxIterExceeded() )
285 p_yp[jj] = p_execute_job( p_xp[jj], jj, p_noptim++ );
291 p_barrier( jlo, jhi );
294 template<
class X,
class Y,
int NP,
int NSTR>
299 const Y F0 = Y(1.4142136);
300 const X
F1 = X(0.7071068);
304 for(
int jj=1; jj <= 2*p_nvar; jj++ )
306 if( p_yp[jj] < p_ymin )
312 bool lgNewCnt = p_jmin > 0;
315 bool lgNegd2 =
false;
318 for(
int i=0; i < p_nvar; i++ )
320 Y d1 = p_yp[2*i+2] - p_yp[2*i+1];
321 Y d2 = Y(0.5)*p_yp[2*i+2] - p_yp[0] + Y(0.5)*p_yp[2*i+1];
324 xhlp[i] = -p_dmax*p_c2[i]*(
static_cast<X
>(
max(
min((Y(0.25)*d1)/
max(d2,Y(1.e-10)),F0),-F0)) -
325 p_delta(2*i+1,p_jmin) + p_delta(2*i+2,p_jmin));
326 xnrm +=
pow2(xhlp[i]);
328 xnrm =
static_cast<X
>(sqrt(xnrm));
332 for(
int j=0; j < p_nvar; j++ )
334 for(
int i=0; i < p_nvar; i++ )
340 p_a2[0][i] *= xhlp[0];
344 p_a2[0][i] += xhlp[j]*p_a2[j][i];
345 p_a2[j][i] = p_delta(i,j);
346 if( j == p_nvar-1 && abs(p_a2[0][i]) > amax )
349 amax = abs(p_a2[0][i]);
355 p_a2[j][i] = p_delta(i,j);
362 p_a2[imax][0] = X(1.);
363 p_a2[imax][imax] = X(0.);
366 p_phygrm( p_a2, p_nvar );
369 for(
int i=0; i < p_nvar; i++ )
372 for(
int j=0; j < p_nvar; j++ )
374 p_c2[i] +=
pow2(p_a2[i][j]/p_c1[j]);
376 p_c2[i] =
static_cast<X
>(1./sqrt(p_c2[i]));
377 p_xc[i] = p_xp[p_jmin][i];
378 p_xp[0][i] = p_xc[i];
380 p_yp[0] = p_yp[p_jmin];
407 p_dmax =
min(
max(xnrm/p_c2[0],p_dmax*r1),p_dmax*r2);
409 p_dmax =
MIN2(p_dmax,p_dold);
412 template<
class X,
class Y,
int NP,
int NSTR>
417 if( !lgConvergedRestart() )
420 for(
int i=0; i < p_nvar; i++ )
422 p_xcold[i] = p_xc[i];
426 p_reset_transformation_matrix();
430 template<
class X,
class Y,
int NP,
int NSTR>
437 for(
int i=0; i < n; i++ )
440 for(
int k=0; k < n; k++ )
443 for(
int k=0; k < n; k++ )
445 for(
int j=i+1; j < n; j++ )
448 for(
int k=0; k < n; k++ )
449 ip += a[i][k]*a[j][k];
450 for(
int k=0; k < n; k++ )
451 a[j][k] -= ip*a[i][k];
457 template<
class X,
class Y,
int NP,
int NSTR>
462 for(
int i=0; i < p_nvar; i++ )
464 if( x[i] < p_absmin[i] )
466 if( x[i] > p_absmax[i] )
472 template<
class X,
class Y,
int NP,
int NSTR>
478 for(
int i=0; i < p_nvar; i++ )
479 for(
int j=0; j < p_nvar; j++ )
480 p_a2[j][i] = p_delta(i,j);
483 template<
class X,
class Y,
int NP,
int NSTR>
490 ASSERT( !lgInitialized() );
492 for(
int i=0; i < nvar; i++ )
494 p_absmin[i] = pmin[i];
495 p_absmax[i] = pmax[i];
499 template<
class X,
class Y,
int NP,
int NSTR>
502 const char* host_name)
508 strncpy( p_chStr1, date, NSTR-1 );
509 if( version != NULL )
510 strncpy( p_chStr2, version, NSTR-1 );
511 if( host_name != NULL )
512 strncpy( p_chStr3, host_name, NSTR-1 );
515 template<
class X,
class Y,
int NP,
int NSTR>
521 strncpy( p_chState, fnam, NSTR-1 );
524 template<
class X,
class Y,
int NP,
int NSTR>
536 ASSERT( nvar > 0 && nvar <= NP );
548 for(
int i=0; i < p_nvar; i++ )
549 p_dmax =
max(p_dmax,abs(del[i]));
552 for(
int i=0; i < p_nvar; i++ )
555 p_xcold[i] = p_xc[i] + X(10.)*p_toler;
556 p_c1[i] = abs(del[i])/p_dmax;
558 p_xp[0][i] = p_xc[i];
559 p_varmax[i] =
max(p_varmax[i],p_xc[i]);
560 p_varmin[i] =
min(p_varmin[i],p_xc[i]);
565 p_yp[0] = p_execute_job( p_xc, 0, p_noptim++ );
571 p_reset_transformation_matrix();
573 p_wr_state( p_chState );
576 template<
class X,
class Y,
int NP,
int NSTR>
591 printf(
"optimize continue - file has incompatible version, sorry\n" );
596 printf(
"optimize continue - arrays have wrong dimension, sorry\n" );
601 printf(
"optimize continue - strings have wrong length, sorry\n" );
606 printf(
"optimize continue - wrong number of free parameters, sorry\n" );
621 template<
class X,
class Y,
int NP,
int NSTR>
626 ASSERT( lgInitialized() );
628 while( !lgConverged() )
630 p_evaluate_hyperblock();
631 if( lgMaxIterExceeded() )
633 p_setup_next_hyperblock();
634 p_wr_state( p_chState );
638 template<
class X,
class Y,
int NP,
int NSTR>
643 ASSERT( lgInitialized() );
645 while( !lgConvergedRestart() )
648 if( lgMaxIterExceeded() )
650 p_reset_hyperblock();
654 template<
class X,
class Y,
int NP,
int NSTR>
662 for(
int i=0; i < p_nvar; i++ )
663 dist +=
pow2(p_xc[i] - p_xcold[i]);
664 dist =
static_cast<X
>(sqrt(dist));
665 return ( dist <= p_toler );
680 fprintf(
ioQQQ,
"optimize_phymir: too many parameters are varied, increase LIMPAR\n" );
722 fprintf(
ioQQQ,
" Optimizer exceeding maximum iterations.\n" );
723 fprintf(
ioQQQ,
" This can be reset with the OPTIMIZE ITERATIONS command.\n" );
728 for(
int i=0; i < nvarPhymir; i++ )
730 xc[i] = phymir.
xval(i);
734 *ymin = phymir.
yval();
746 if( fwrite(ptr,len,
size_t(1),fdes) != 1 ) {
747 printf(
"error writing on file: %s\n",fnam );
763 if( fread(ptr,len,
size_t(1),fdes) != 1 ) {
764 printf(
"error reading on file: %s\n",fnam );
void p_rd_state(const char *)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
void set_used_nCPU(long n)
NORETURN void TotalInsanity(void)
void rd_block(void *, size_t, const char *)
Y p_execute_job(const X[], int, int)
chi2_type optimize_func(const realnum param[], int grid_index=-1)
realnum varang[LIMPAR][2]
void p_wr_state(const char *) const
void optimize_phymir(realnum[], const realnum[], long, chi2_type *, realnum)
void p_execute_job_parallel(const X[], int, int) const
static t_version & Inst()
void init_strings(const char *, const char *, const char *)
void continue_from_state(Y(*)(const X[], int), int, const char *, X, int, phymir_mode, int)
void wr_block(const void *, size_t, const char *)
bool fp_equal(sys_float x, sys_float y, int n=3)
void init_minmax(const X[], const X[], int)
bool lgConvergedRestart() const
STATIC double dist(long, realnum[], realnum[])
void p_reset_transformation_matrix()
void append_file(FILE *dest, const char *source)
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
void init_state_file_name(const char *)
#define DEBUG_ENTRY(funcname)
bool p_lgLimitExceeded(const X[]) const
bool lgMaxIterExceeded() const
void p_setup_next_hyperblock()
int fprintf(const Output &stream, const char *format,...)
void p_reset_hyperblock()
void optimize_with_restart()
void p_evaluate_hyperblock()
void p_process_output(int, int)
const char * STATEFILE_BACKUP
const char * host_name() const
void initial_run(Y(*)(const X[], int), int, const X[], const X[], X, int, phymir_mode, int)
#define MPI_Bcast(V, W, X, Y, Z)
void p_phygrm(X[][NP], int)