00001
00002
00003
00004
00005 #if defined(unix) || defined(__unix) || defined(__unix__)
00006
00007
00008
00009
00010
00011 #define _INCLUDE_POSIX_SOURCE
00012 #endif
00013 #include "cddefines.h"
00014 #include "version.h"
00015 #include "optimize.h"
00016 #ifdef __unix
00017 #include <stddef.h>
00018 #include <sys/types.h>
00019 #include <sys/wait.h>
00020 #include <unistd.h>
00021 #endif
00022
00023 #define F0 1.4142136f
00024 #define F1 0.7071068f
00025 #define F2 0.1f
00026
00027
00028
00029
00030
00031
00032
00033 static void phygrm(float a2[][LIMPAR],long);
00034 static void wr_continue(float,long,float a2[][LIMPAR],const float[],const float[],
00035 const float[],const float[],float,float,float,long,long,const float[],
00036 const float[],const char[],const char[],const char[],const char*);
00037 static void rd_continue(float*,long*,float a2[][LIMPAR],float[],float[],float[],float[],
00038 float*,float*,float*,long*,long*,float[],float[],char[],char[],char[],const char*);
00039 #ifdef __unix
00040 static void cpfile(const char*);
00041 static void wr_block(void*,size_t,const char*);
00042 static void rd_block(void*,size_t,const char*);
00043 #endif
00044
00045
00046 void optimize_phymir(float xc[],
00047 float del[],
00048 long int nvarPhymir,
00049 float *ymin,
00050 float toler)
00051 {
00052 #define DELTA(i,j) (((i) == (j)) ? 1.f : 0.f)
00053
00054 char chDum1[STDLEN],
00055 chDum2[STDLEN],
00056 chDum3[STDLEN];
00057 bool lgFinish,
00058 lgNegd2,
00059 lgNewCnt,
00060 lgRead;
00061 long int i,
00062 imax,
00063 j,
00064 jj,
00065 jmin=0,
00066 limcp,
00067 nvarcp;
00068 float a2[LIMPAR][LIMPAR],
00069 amax,
00070 c1[LIMPAR],
00071 c2[LIMPAR],
00072 d1,
00073 d2,
00074 diff,
00075 dmax=0.f,
00076 dold=0.f,
00077 r1,
00078 r2,
00079 sgn,
00080 temp,
00081 vers,
00082 vpusedPhymir,
00083 xcold[LIMPAR],
00084 xhlp[LIMPAR],
00085 xnrm,
00086 xp[LIMPAR][2*LIMPAR + 1],
00087 yp[2*LIMPAR + 1];
00088 # ifdef __unix
00089 long int currentCPU;
00090 int stat;
00091 char fnam1[20],
00092 fnam2[20];
00093 pid_t pid;
00094 # endif
00095
00096 DEBUG_ENTRY( "optimize_phymir()" );
00097
00098 if( nvarPhymir > LIMPAR )
00099 {
00100 fprintf( ioQQQ, "optimize_phymir: too many parameters are varied, increase LIMPAR\n" );
00101 puts( "[Stop]" );
00102 cdEXIT(EXIT_FAILURE);
00103 }
00104
00105
00106 memset( chDum1, '\0', STDLEN );
00107 memset( chDum2, '\0', STDLEN );
00108 memset( chDum3, '\0', STDLEN );
00109
00110 for( i=0; i < LIMPAR; ++i )
00111 {
00112 c1[i] = -FLT_MAX;
00113 c2[i] = -FLT_MAX;
00114 xcold[i] = -FLT_MAX;
00115 for( j=0; j < LIMPAR; ++j )
00116 {
00117 a2[i][j] = -FLT_MAX;
00118 }
00119 }
00120
00121 optimize.nOptimiz = 0;
00122 # ifdef __unix
00123 currentCPU = 0;
00124 # endif
00125 lgFinish = false;
00126 lgRead = optimize.lgOptCont;
00127 if( optimize.lgOptCont )
00128 goto L_20;
00129
00130
00131 dmax = 0.f;
00132 for( i=0; i < nvarPhymir; i++ )
00133 {
00134 temp = (float)fabs(del[i]);
00135 dmax = MAX2(dmax,temp);
00136 }
00137
00138 dold = dmax;
00139 for( i=0; i < nvarPhymir; i++ )
00140 {
00141 xcold[i] = xc[i] + 10.f*toler;
00142 c1[i] = (float)fabs(del[i])/dmax;
00143 c2[i] = c1[i];
00144 xp[i][0] = xc[i];
00145 optimize.vparm[0][i] = xc[i];
00146
00147 if( optimize.lgParallel )
00148 {
00149 vpusedPhymir = MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
00150 vpusedPhymir = MAX2(vpusedPhymir,optimize.varang[i][0]);
00151 optimize.varmax[i] = MAX2(optimize.varmax[i],vpusedPhymir);
00152 optimize.varmin[i] = MIN2(optimize.varmin[i],vpusedPhymir);
00153 }
00154 }
00155 # ifdef __unix
00156 if( optimize.lgParallel )
00157 {
00158
00159 fflush(NULL);
00160 pid = fork();
00161 if( pid == 0 )
00162 {
00163
00164 sprintf(fnam1,"chi2_%ld",0L);
00165 sprintf(fnam2,"output_%ld",0L);
00166
00167 ioQQQ = fopen(fnam2,"w");
00168
00169 yp[0] = FLT_MAX;
00170 wr_block(&yp[0],(size_t)sizeof(float),fnam1);
00171 yp[0] = (float)optimize_func(xc);
00172 wr_block(&yp[0],(size_t)sizeof(float),fnam1);
00173 cdEXIT(EXIT_SUCCESS);
00174 }
00175 else
00176 {
00177
00178 pid = wait(&stat);
00179 sprintf(fnam1,"chi2_%ld",0L);
00180 sprintf(fnam2,"output_%ld",0L);
00181 rd_block(&yp[0],(size_t)sizeof(float),fnam1);
00182 remove(fnam1);
00183 cpfile(fnam2);
00184 remove(fnam2);
00185 }
00186
00187 optimize.nOptimiz++;
00188 }
00189 else
00190 {
00191 yp[0] = (float)optimize_func(xc);
00192 }
00193 # else
00194
00195 yp[0] = (float)optimize_func(xc);
00196 # endif
00197 *ymin = yp[0];
00198 jmin = 0;
00199
00200 L_10:
00201 for( i=0; i < nvarPhymir; i++ )
00202 {
00203 for( j=0; j < nvarPhymir; j++ )
00204 {
00205 a2[j][i] = DELTA(i,j);
00206 }
00207 }
00208
00209 L_20:
00210 if( lgRead )
00211 {
00212 rd_continue(&vers,&limcp,a2,c1,c2,xc,xcold,&dmax,&dold,ymin,&nvarcp,&optimize.nOptimiz,
00213 optimize.varmax,optimize.varmin,chDum1,chDum2,chDum3,CNTFILE);
00214 if( nvarcp != nvarPhymir )
00215 {
00216 printf( "optimize continue - wrong number of free parameters, sorry\n" );
00217 puts( "[Stop]" );
00218 cdEXIT(EXIT_FAILURE);
00219 }
00220 for( i=0; i < nvarPhymir; i++ )
00221 {
00222 xp[i][0] = xc[i];
00223 }
00224 yp[0] = *ymin;
00225 jmin = 0;
00226 lgRead = false;
00227 }
00228 else
00229 {
00230 strcpy(chDum1,version.chDate);
00231 strcpy(chDum2,version.chVersion);
00232 strcpy(chDum3,cpu.host_name());
00233 wr_continue(VRSNEW,LIMPAR,a2,c1,c2,xc,xcold,dmax,dold,*ymin,nvarPhymir,optimize.nOptimiz,
00234 optimize.varmax,optimize.varmin,chDum1,chDum2,chDum3,CNTFILE);
00235 }
00236
00237 if( lgFinish )
00238 {
00239 DEBUG_EXIT( "optimize_phymir()" );
00240 return;
00241 }
00242
00243 if( optimize.nOptimiz >= optimize.nIterOptim )
00244 {
00245 fprintf( ioQQQ, " Optimizer exceeding maximum iterations.\n This can be reset with the OPTIMIZE ITERATIONS command.\n" );
00246 DEBUG_EXIT( "optimize_phymir()" );
00247 return;
00248 }
00249
00250 for( j=0; j < nvarPhymir; j++ )
00251 {
00252 sgn = 1.f;
00253 for( jj=2*j+1; jj <= 2*j+2; jj++ )
00254 {
00255 sgn = -sgn;
00256 for( i=0; i < nvarPhymir; i++ )
00257 {
00258 xp[i][jj] = xc[i] + sgn*dmax*c2[j]*a2[j][i];
00259 xhlp[i] = xp[i][jj];
00260 optimize.vparm[0][i] = xhlp[i];
00261
00262 if( optimize.lgParallel )
00263 {
00264 vpusedPhymir = MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
00265 vpusedPhymir = MAX2(vpusedPhymir,optimize.varang[i][0]);
00266 optimize.varmax[i] = MAX2(optimize.varmax[i],vpusedPhymir);
00267 optimize.varmin[i] = MIN2(optimize.varmin[i],vpusedPhymir);
00268 }
00269 }
00270 # ifdef __unix
00271 if( optimize.lgParallel )
00272 {
00273 currentCPU++;
00274 if( currentCPU > optimize.useCPU )
00275 {
00276
00277 pid = wait(&stat);
00278 currentCPU--;
00279 }
00280
00281 fflush(NULL);
00282 pid = fork();
00283 if( pid == 0 )
00284 {
00285
00286 sprintf(fnam1,"chi2_%ld",jj);
00287 sprintf(fnam2,"output_%ld",jj);
00288
00289 ioQQQ = fopen(fnam2,"w");
00290
00291 yp[jj] = FLT_MAX;
00292 wr_block(&yp[jj],(size_t)sizeof(float),fnam1);
00293 yp[jj] = (float)optimize_func(xhlp);
00294 wr_block(&yp[jj],(size_t)sizeof(float),fnam1);
00295 cdEXIT(EXIT_SUCCESS);
00296 }
00297
00298 optimize.nOptimiz++;
00299 }
00300 else
00301 {
00302 yp[jj] = (float)optimize_func(xhlp);
00303 }
00304 # else
00305
00306 yp[jj] = (float)optimize_func(xhlp);
00307 # endif
00308 }
00309 }
00310 # ifdef __unix
00311
00312 if( optimize.lgParallel )
00313 {
00314 while( currentCPU > 0 )
00315 {
00316 pid = wait(&stat);
00317 currentCPU--;
00318 }
00319 }
00320 # endif
00321
00322 for( jj=1; jj <= 2*nvarPhymir; jj++ )
00323 {
00324 # ifdef __unix
00325 if( optimize.lgParallel )
00326 {
00327 sprintf(fnam1,"chi2_%ld",jj);
00328 sprintf(fnam2,"output_%ld",jj);
00329 rd_block(&yp[jj],(size_t)sizeof(float),fnam1);
00330 remove(fnam1);
00331 cpfile(fnam2);
00332 remove(fnam2);
00333 }
00334 # endif
00335 if( yp[jj] < *ymin )
00336 {
00337 *ymin = yp[jj];
00338 jmin = jj;
00339 }
00340 }
00341
00342 lgNewCnt = jmin > 0;
00343
00344 lgNegd2 = false;
00345 xnrm = 0.f;
00346 for( i=0; i < nvarPhymir; i++ )
00347 {
00348 d1 = yp[2*i+2] - yp[2*i+1];
00349 d2 = 0.5f*yp[2*i+2] - yp[0] + 0.5f*yp[2*i+1];
00350 if( d2 <= 0.f )
00351 lgNegd2 = true;
00352 xhlp[i] = -dmax*c2[i]*(MAX2(MIN2((0.25f*d1)/MAX2(d2,1.e-10f),F0),-F0) - DELTA(2*i+1,jmin) + DELTA(2*i+2,jmin));
00353 xnrm += POW2(xhlp[i]);
00354 }
00355 xnrm = (float)sqrt(xnrm);
00356
00357 imax = 0;
00358 amax = 0.f;
00359 for( j=0; j < nvarPhymir; j++ )
00360 {
00361 for( i=0; i < nvarPhymir; i++ )
00362 {
00363 if( xnrm > 0.f )
00364 {
00365 if( j == 0 )
00366 {
00367 a2[0][i] *= xhlp[0];
00368 }
00369 else
00370 {
00371 a2[0][i] += xhlp[j]*a2[j][i];
00372 a2[j][i] = DELTA(i,j);
00373 if( j == nvarPhymir-1 && (float)fabs(a2[0][i]) > amax )
00374 {
00375 imax = i;
00376 amax = (float)fabs(a2[0][i]);
00377 }
00378 }
00379 }
00380 else
00381 {
00382 a2[j][i] = DELTA(i,j);
00383 }
00384 }
00385 }
00386
00387 if( imax > 0 )
00388 {
00389 a2[imax][0] = 1.f;
00390 a2[imax][imax] = 0.f;
00391 }
00392
00393 phygrm(a2,nvarPhymir);
00394
00395
00396 for( i=0; i < nvarPhymir; i++ )
00397 {
00398 c2[i] = 0.f;
00399 for( j=0; j < nvarPhymir; j++ )
00400 {
00401 c2[i] += POW2(a2[i][j]/c1[j]);
00402 }
00403 c2[i] = 1.f/(float)sqrt(c2[i]);
00404 xc[i] = xp[i][jmin];
00405 xp[i][0] = xc[i];
00406 }
00407 yp[0] = yp[jmin];
00408 jmin = 0;
00409
00410 if( lgNegd2 )
00411 {
00412
00413
00414
00415 r1 = F1;
00416 r2 = F1;
00417 }
00418 else
00419 {
00420 r1 = F2;
00421 if( lgNewCnt )
00422 {
00423
00424 r2 = (float)sqrt(1.f/F1);
00425 }
00426 else
00427 {
00428
00429
00430 r2 = (float)sqrt(F1);
00431 }
00432 }
00433 dmax = MIN2(MAX2(xnrm/c2[0],dmax*r1),dmax*r2);
00434
00435 dmax = MIN2(dmax,dold);
00436
00437 if( dmax > toler )
00438 goto L_20;
00439
00440 diff = 0.f;
00441 for( i=0; i < nvarPhymir; i++ )
00442 {
00443 diff += POW2(xc[i] - xcold[i]);
00444 xcold[i] = xc[i];
00445 c2[i] = c1[i];
00446 }
00447 diff = (float)sqrt(diff);
00448 dmax = dold;
00449 lgFinish = diff <= toler;
00450 goto L_10;
00451 # undef DELTA
00452 }
00453
00454 static void phygrm(float a[][LIMPAR],
00455 long int n)
00456 {
00457 long int i,
00458 j,
00459 k;
00460 float ip;
00461
00462 DEBUG_ENTRY( "phygrm()" );
00463
00464
00465 for( i=0; i < n; i++ )
00466 {
00467 ip = 0.f;
00468 for( k=0; k < n; k++ )
00469 ip += POW2(a[i][k]);
00470 ip = (float)sqrt(ip);
00471 for( k=0; k < n; k++ )
00472 a[i][k] /= ip;
00473 for( j=i+1; j < n; j++ )
00474 {
00475 ip = 0.f;
00476 for( k=0; k < n; k++ )
00477 ip += a[i][k]*a[j][k];
00478 for( k=0; k < n; k++ )
00479 a[j][k] -= ip*a[i][k];
00480 }
00481 }
00482
00483 DEBUG_EXIT( "phygrm()" );
00484 return;
00485 }
00486
00487 static void wr_continue(float vers,
00488 long dim,
00489 float a2[][LIMPAR],
00490 const float c1[],
00491 const float c2[],
00492 const float xc[],
00493 const float xcold[],
00494 float dmax,
00495 float dold,
00496 float ymin,
00497 long nvar,
00498 long noptim,
00499 const float varmax[],
00500 const float varmin[],
00501 const char chDum1[],
00502 const char chDum2[],
00503 const char chDum3[],
00504 const char *fnam)
00505 {
00506 bool lgErr;
00507 int res;
00508 size_t num;
00509 FILE *fdes;
00510
00511 DEBUG_ENTRY( "wr_continue()" );
00512
00513 if( (fdes = fopen(fnam,"wb")) == 0 ) {
00514 printf( "error opening file: %s\n",fnam );
00515 return;
00516 }
00517 lgErr = false;
00518 num = 1;
00519 lgErr = lgErr || ( fwrite(&vers,sizeof(float),num,fdes) != num );
00520 lgErr = lgErr || ( fwrite(&dim,sizeof(long),num,fdes) != num );
00521 num = (size_t)(LIMPAR*LIMPAR);
00522 lgErr = lgErr || ( fwrite(a2,sizeof(float),num,fdes) != num );
00523 num = (size_t)LIMPAR;
00524 lgErr = lgErr || ( fwrite(c1,sizeof(float),num,fdes) != num );
00525 lgErr = lgErr || ( fwrite(c2,sizeof(float),num,fdes) != num );
00526 lgErr = lgErr || ( fwrite(xc,sizeof(float),num,fdes) != num );
00527 lgErr = lgErr || ( fwrite(xcold,sizeof(float),num,fdes) != num );
00528 num = 1;
00529 lgErr = lgErr || ( fwrite(&dmax,sizeof(float),num,fdes) != num );
00530 lgErr = lgErr || ( fwrite(&dold,sizeof(float),num,fdes) != num );
00531 lgErr = lgErr || ( fwrite(&ymin,sizeof(float),num,fdes) != num );
00532 lgErr = lgErr || ( fwrite(&nvar,sizeof(long),num,fdes) != num );
00533 lgErr = lgErr || ( fwrite(&noptim,sizeof(long),num,fdes) != num );
00534 num = (size_t)LIMPAR;
00535 lgErr = lgErr || ( fwrite(varmax,sizeof(float),num,fdes) != num );
00536 lgErr = lgErr || ( fwrite(varmin,sizeof(float),num,fdes) != num );
00537 num = (size_t)STDLEN;
00538 lgErr = lgErr || ( fwrite(chDum1,sizeof(char),num,fdes) != num );
00539 lgErr = lgErr || ( fwrite(chDum2,sizeof(char),num,fdes) != num );
00540 lgErr = lgErr || ( fwrite(chDum3,sizeof(char),num,fdes) != num );
00541 fclose(fdes);
00542 if( lgErr )
00543 {
00544 printf( "error writing file: %s\n",fnam );
00545 res = remove(fnam);
00546
00547 if( false) fprintf(ioQQQ,"%i\n", res );
00548 }
00549
00550 DEBUG_EXIT( "wr_continue()" );
00551 return;
00552 }
00553
00554 static void rd_continue(float *vers,
00555 long *dim,
00556 float a2[][LIMPAR],
00557 float c1[],
00558 float c2[],
00559 float xc[],
00560 float xcold[],
00561 float *dmax,
00562 float *dold,
00563 float *ymin,
00564 long *nvar,
00565 long *noptim,
00566 float varmax[],
00567 float varmin[],
00568 char chDum1[],
00569 char chDum2[],
00570 char chDum3[],
00571 const char *fnam)
00572 {
00573 bool lgErr;
00574 size_t num;
00575 FILE *fdes;
00576
00577 DEBUG_ENTRY( "rd_continue()" );
00578
00579 if( (fdes = fopen(fnam,"rb")) == 0 ) {
00580 printf( "error opening file: %s\n",fnam );
00581 puts( "[Stop]" );
00582 cdEXIT(EXIT_FAILURE);
00583 }
00584 lgErr = false;
00585 num = 1;
00586 lgErr = lgErr || ( fread(vers,sizeof(float),num,fdes) != num );
00587 if( lgErr || *vers != VRSNEW )
00588 {
00589 printf( "optimize continue - file has incompatible version, sorry\n" );
00590 puts( "[Stop]" );
00591 cdEXIT(EXIT_FAILURE);
00592 }
00593 lgErr = lgErr || ( fread(dim,sizeof(long),num,fdes) != num );
00594 if( lgErr || *dim != LIMPAR )
00595 {
00596 printf( "optimize continue - arrays have wrong dimension, sorry\n" );
00597 puts( "[Stop]" );
00598 cdEXIT(EXIT_FAILURE);
00599 }
00600 num = (size_t)(LIMPAR*LIMPAR);
00601 lgErr = lgErr || ( fread(a2,sizeof(float),num,fdes) != num );
00602 num = (size_t)LIMPAR;
00603 lgErr = lgErr || ( fread(c1,sizeof(float),num,fdes) != num );
00604 lgErr = lgErr || ( fread(c2,sizeof(float),num,fdes) != num );
00605 lgErr = lgErr || ( fread(xc,sizeof(float),num,fdes) != num );
00606 lgErr = lgErr || ( fread(xcold,sizeof(float),num,fdes) != num );
00607 num = 1;
00608 lgErr = lgErr || ( fread(dmax,sizeof(float),num,fdes) != num );
00609 lgErr = lgErr || ( fread(dold,sizeof(float),num,fdes) != num );
00610 lgErr = lgErr || ( fread(ymin,sizeof(float),num,fdes) != num );
00611 lgErr = lgErr || ( fread(nvar,sizeof(long),num,fdes) != num );
00612 lgErr = lgErr || ( fread(noptim,sizeof(long),num,fdes) != num );
00613 num = (size_t)LIMPAR;
00614 lgErr = lgErr || ( fread(varmax,sizeof(float),num,fdes) != num );
00615 lgErr = lgErr || ( fread(varmin,sizeof(float),num,fdes) != num );
00616 num = (size_t)STDLEN;
00617 lgErr = lgErr || ( fread(chDum1,sizeof(char),num,fdes) != num );
00618 lgErr = lgErr || ( fread(chDum2,sizeof(char),num,fdes) != num );
00619 lgErr = lgErr || ( fread(chDum3,sizeof(char),num,fdes) != num );
00620 fclose(fdes);
00621 if( lgErr )
00622 {
00623 printf( "error reading file: %s\n",fnam );
00624 puts( "[Stop]" );
00625 cdEXIT(EXIT_FAILURE);
00626 }
00627
00628 DEBUG_EXIT( "rd_continue()" );
00629 return;
00630 }
00631
00632 #ifdef __unix
00633
00634 static void cpfile(const char *fnam)
00635 {
00636 char chr;
00637
00638 FILE *fdes;
00639
00640 DEBUG_ENTRY( "cpfile()" );
00641
00642
00643
00644 if( (fdes = fopen(fnam,"r")) == 0 )
00645 return;
00646 chr = (char)fgetc(fdes);
00647 while( ! feof(fdes) )
00648 {
00649 fputc( chr , ioQQQ );
00650 chr = (char)fgetc(fdes);
00651 }
00652 fclose(fdes);
00653
00654 DEBUG_EXIT( "cpfile()" );
00655 return;
00656 }
00657
00658 static void wr_block(void *ptr,
00659 size_t len,
00660 const char *fnam)
00661 {
00662 FILE *fdes;
00663
00664 DEBUG_ENTRY( "wr_block()" );
00665
00666
00667
00668 if( (fdes = fopen(fnam,"wb")) == 0 ) {
00669 printf( "error opening file: %s\n",fnam );
00670 puts( "[Stop]" );
00671 cdEXIT(EXIT_FAILURE);
00672 }
00673 if( fwrite(ptr,len,(size_t)1,fdes) != 1 ) {
00674 printf( "error writing on file: %s\n",fnam );
00675 fclose(fdes);
00676 puts( "[Stop]" );
00677 cdEXIT(EXIT_FAILURE);
00678 }
00679 fclose(fdes);
00680
00681 DEBUG_EXIT( "wr_block()" );
00682 return;
00683 }
00684
00685 static void rd_block(void *ptr,
00686 size_t len,
00687 const char *fnam)
00688 {
00689 FILE *fdes;
00690
00691 DEBUG_ENTRY( "rd_block()" );
00692
00693
00694
00695 if( (fdes = fopen(fnam,"rb")) == 0 ) {
00696 printf( "error opening file: %s\n",fnam );
00697 puts( "[Stop]" );
00698 cdEXIT(EXIT_FAILURE);
00699 }
00700 if( fread(ptr,len,(size_t)1,fdes) != 1 ) {
00701 printf( "error reading on file: %s\n",fnam );
00702 fclose(fdes);
00703 puts( "[Stop]" );
00704 cdEXIT(EXIT_FAILURE);
00705 }
00706 fclose(fdes);
00707
00708 DEBUG_EXIT( "rd_block()" );
00709 return;
00710 }
00711
00712 #endif