5 #if defined(unix) || defined(__unix) || defined(__unix__)
11 #define _INCLUDE_POSIX_SOURCE
18 #include <sys/types.h>
36 const realnum[],
const char[],
const char[],
const char[],
const char*);
38 realnum*,
realnum*,
realnum*,
long*,
long*,
realnum[],
realnum[],
char[],
char[],
char[],
const char*);
40 STATIC void cpfile(
const char*);
41 STATIC void wr_block(
void*,
size_t,
const char*);
42 STATIC void rd_block(
void*,
size_t,
const char*);
52 #define DELTA(i,j) (((i) == (j)) ? 1.f : 0.f)
100 fprintf(
ioQQQ,
"optimize_phymir: too many parameters are varied, increase LIMPAR\n" );
105 memset( chDum1,
'\0',
STDLEN );
106 memset( chDum2,
'\0',
STDLEN );
107 memset( chDum3,
'\0',
STDLEN );
109 for( i=0; i <
LIMPAR; ++i )
114 for( j=0; j <
LIMPAR; ++j )
131 for( i=0; i < nvarPhymir; i++ )
134 dmax =
MAX2(dmax,temp);
138 for( i=0; i < nvarPhymir; i++ )
140 xcold[i] = xc[i] + 10.f*toler;
141 c1[i] = (
realnum)fabs(del[i])/dmax;
163 sprintf(fnam1,
"chi2_%ld",0L);
164 sprintf(fnam2,
"output_%ld",0L);
169 wr_block(&yp[0],(
size_t)
sizeof(
realnum),fnam1);
171 wr_block(&yp[0],(
size_t)
sizeof(
realnum),fnam1);
180 sprintf(fnam1,
"chi2_%ld",0L);
181 sprintf(fnam2,
"output_%ld",0L);
182 rd_block(&yp[0],(
size_t)
sizeof(
realnum),fnam1);
202 for( i=0; i < nvarPhymir; i++ )
204 for( j=0; j < nvarPhymir; j++ )
206 a2[j][i] =
DELTA(i,j);
213 rd_continue(&vers,&limcp,a2,c1,c2,xc,xcold,&dmax,&dold,ymin,&nvarcp,&
optimize.
nOptimiz,
215 if( nvarcp != nvarPhymir )
217 printf(
"optimize continue - wrong number of free parameters, sorry\n" );
220 for( i=0; i < nvarPhymir; i++ )
234 wr_continue(
VRSNEW,LIMPAR,a2,c1,c2,xc,xcold,dmax,dold,*ymin,nvarPhymir,
optimize.
nOptimiz,
245 fprintf(
ioQQQ,
" Optimizer exceeding maximum iterations.\n This can be reset with the OPTIMIZE ITERATIONS command.\n" );
249 for( j=0; j < nvarPhymir; j++ )
252 for( jj=2*j+1; jj <= 2*j+2; jj++ )
255 for( i=0; i < nvarPhymir; i++ )
257 xp[i][jj] = xc[i] + sgn*dmax*c2[j]*a2[j][i];
285 sprintf(fnam1,
"chi2_%ld",jj);
286 sprintf(fnam2,
"output_%ld",jj);
291 wr_block(&yp[jj],(
size_t)
sizeof(
realnum),fnam1);
293 wr_block(&yp[jj],(
size_t)
sizeof(
realnum),fnam1);
315 while( currentCPU > 0 )
323 for( jj=1; jj <= 2*nvarPhymir; jj++ )
328 sprintf(fnam1,
"chi2_%ld",jj);
329 sprintf(fnam2,
"output_%ld",jj);
330 rd_block(&yp[jj],(
size_t)
sizeof(
realnum),fnam1);
347 for( i=0; i < nvarPhymir; i++ )
349 d1 = yp[2*i+2] - yp[2*i+1];
350 d2 = 0.5f*yp[2*i+2] - yp[0] + 0.5f*yp[2*i+1];
354 xnrm +=
POW2(xhlp[i]);
360 for( j=0; j < nvarPhymir; j++ )
362 for( i=0; i < nvarPhymir; i++ )
372 a2[0][i] += xhlp[j]*a2[j][i];
373 a2[j][i] =
DELTA(i,j);
374 if( j == nvarPhymir-1 && (
realnum)fabs(a2[0][i]) > amax )
377 amax = (
realnum)fabs(a2[0][i]);
383 a2[j][i] =
DELTA(i,j);
391 a2[imax][imax] = 0.f;
397 for( i=0; i < nvarPhymir; i++ )
400 for( j=0; j < nvarPhymir; j++ )
402 c2[i] +=
POW2(a2[i][j]/c1[j]);
404 c2[i] = 1.f/sqrt(c2[i]);
434 dmax =
MIN2(
MAX2(xnrm/c2[0],dmax*r1),dmax*r2);
436 dmax =
MIN2(dmax,dold);
442 for( i=0; i < nvarPhymir; i++ )
444 diff +=
POW2(xc[i] - xcold[i]);
450 lgFinish = diff <= toler;
466 for( i=0; i < n; i++ )
469 for( k=0; k < n; k++ )
472 for( k=0; k < n; k++ )
474 for( j=i+1; j < n; j++ )
477 for( k=0; k < n; k++ )
478 ip += a[i][k]*a[j][k];
479 for( k=0; k < n; k++ )
480 a[j][k] -= ip*a[i][k];
522 lgErr = lgErr || ( fwrite(&vers,
sizeof(
realnum),num,fdes) != num );
523 lgErr = lgErr || ( fwrite(&dim,
sizeof(
long),num,fdes) != num );
524 num = (size_t)(LIMPAR*LIMPAR);
525 lgErr = lgErr || ( fwrite(
a2,
sizeof(
realnum),num,fdes) != num );
526 num = (size_t)LIMPAR;
527 lgErr = lgErr || ( fwrite(c1,
sizeof(
realnum),num,fdes) != num );
528 lgErr = lgErr || ( fwrite(c2,
sizeof(
realnum),num,fdes) != num );
529 lgErr = lgErr || ( fwrite(xc,
sizeof(
realnum),num,fdes) != num );
530 lgErr = lgErr || ( fwrite(xcold,
sizeof(
realnum),num,fdes) != num );
532 lgErr = lgErr || ( fwrite(&dmax,
sizeof(
realnum),num,fdes) != num );
533 lgErr = lgErr || ( fwrite(&dold,
sizeof(
realnum),num,fdes) != num );
534 lgErr = lgErr || ( fwrite(&ymin,
sizeof(
realnum),num,fdes) != num );
535 lgErr = lgErr || ( fwrite(&nvar,
sizeof(
long),num,fdes) != num );
536 lgErr = lgErr || ( fwrite(&noptim,
sizeof(
long),num,fdes) != num );
537 num = (size_t)LIMPAR;
538 lgErr = lgErr || ( fwrite(varmax,
sizeof(
realnum),num,fdes) != num );
539 lgErr = lgErr || ( fwrite(varmin,
sizeof(
realnum),num,fdes) != num );
541 lgErr = lgErr || ( fwrite(chDum1,
sizeof(
char),num,fdes) != num );
542 lgErr = lgErr || ( fwrite(chDum2,
sizeof(
char),num,fdes) != num );
543 lgErr = lgErr || ( fwrite(chDum3,
sizeof(
char),num,fdes) != num );
547 printf(
"error writing file: %s\n",fnam );
550 if(
false) fprintf(
ioQQQ,
"%i\n", res );
583 lgErr = lgErr || ( fread(vers,
sizeof(
realnum),num,fdes) != num );
584 if( lgErr || *vers !=
VRSNEW )
586 printf(
"optimize continue - file has incompatible version, sorry\n" );
589 lgErr = lgErr || ( fread(dim,
sizeof(
long),num,fdes) != num );
590 if( lgErr || *dim != LIMPAR )
592 printf(
"optimize continue - arrays have wrong dimension, sorry\n" );
595 num = (size_t)(LIMPAR*LIMPAR);
596 lgErr = lgErr || ( fread(
a2,
sizeof(
realnum),num,fdes) != num );
597 num = (size_t)LIMPAR;
598 lgErr = lgErr || ( fread(c1,
sizeof(
realnum),num,fdes) != num );
599 lgErr = lgErr || ( fread(c2,
sizeof(
realnum),num,fdes) != num );
600 lgErr = lgErr || ( fread(xc,
sizeof(
realnum),num,fdes) != num );
601 lgErr = lgErr || ( fread(xcold,
sizeof(
realnum),num,fdes) != num );
603 lgErr = lgErr || ( fread(dmax,
sizeof(
realnum),num,fdes) != num );
604 lgErr = lgErr || ( fread(dold,
sizeof(
realnum),num,fdes) != num );
605 lgErr = lgErr || ( fread(ymin,
sizeof(
realnum),num,fdes) != num );
606 lgErr = lgErr || ( fread(nvar,
sizeof(
long),num,fdes) != num );
607 lgErr = lgErr || ( fread(noptim,
sizeof(
long),num,fdes) != num );
608 num = (size_t)LIMPAR;
609 lgErr = lgErr || ( fread(varmax,
sizeof(
realnum),num,fdes) != num );
610 lgErr = lgErr || ( fread(varmin,
sizeof(
realnum),num,fdes) != num );
612 lgErr = lgErr || ( fread(chDum1,
sizeof(
char),num,fdes) != num );
613 lgErr = lgErr || ( fread(chDum2,
sizeof(
char),num,fdes) != num );
614 lgErr = lgErr || ( fread(chDum3,
sizeof(
char),num,fdes) != num );
618 printf(
"error reading file: %s\n",fnam );
626 STATIC void cpfile(
const char *fnam)
639 chr = (char)fgetc(fdes);
640 while( ! feof(fdes) )
642 fputc( chr ,
ioQQQ );
643 chr = (char)fgetc(fdes);
649 STATIC void wr_block(
void *ptr,
660 if( fwrite(ptr,len,(
size_t)1,fdes) != 1 ) {
661 printf(
"error writing on file: %s\n",fnam );
669 STATIC void rd_block(
void *ptr,
680 if( fread(ptr,len,(
size_t)1,fdes) != 1 ) {
681 printf(
"error reading on file: %s\n",fnam );