118 lines
2.6 KiB
C
118 lines
2.6 KiB
C
|
|
#include <math.h>
|
|
#include "nrutil.h"
|
|
#define PFAC 0.1
|
|
#define MNPT 15
|
|
#define MNBS 60
|
|
#define TINY 1.0e-30
|
|
#define BIG 1.0e30
|
|
|
|
static long iran=0;
|
|
|
|
void miser(func,regn,ndim,npts,dith,ave,var)
|
|
float (*func)(),*ave,*var,dith,regn[];
|
|
int ndim;
|
|
unsigned long npts;
|
|
{
|
|
void ranpt();
|
|
float *regn_temp;
|
|
unsigned long n,npre,nptl,nptr;
|
|
int j,jb;
|
|
float avel,varl;
|
|
float fracl,fval;
|
|
float rgl,rgm,rgr,s,sigl,siglb,sigr,sigrb;
|
|
float sum,sumb,summ,summ2;
|
|
float *fmaxl,*fmaxr,*fminl,*fminr;
|
|
float *pt,*rmid;
|
|
|
|
pt=vector(1,ndim);
|
|
if (npts < MNBS) {
|
|
summ=summ2=0.0;
|
|
for (n=1;n<=npts;n++) {
|
|
ranpt(pt,regn,ndim);
|
|
fval=(*func)(pt);
|
|
summ += fval;
|
|
summ2 += fval * fval;
|
|
}
|
|
*ave=summ/npts;
|
|
*var=FMAX(TINY,(summ2-summ*summ/npts)/(npts*npts));
|
|
}
|
|
else {
|
|
rmid=vector(1,ndim);
|
|
npre=LMAX((unsigned long)(npts*PFAC),MNPT);
|
|
fmaxl=vector(1,ndim);
|
|
fmaxr=vector(1,ndim);
|
|
fminl=vector(1,ndim);
|
|
fminr=vector(1,ndim);
|
|
for (j=1;j<=ndim;j++) {
|
|
iran=(iran*2661+36979) % 175000;
|
|
s=SIGN(dith,(float)(iran-87500));
|
|
rmid[j]=(0.5+s)*regn[j]+(0.5-s)*regn[ndim+j];
|
|
fminl[j]=fminr[j]=BIG;
|
|
fmaxl[j]=fmaxr[j] = -BIG;
|
|
}
|
|
for (n=1;n<=npre;n++) {
|
|
ranpt(pt,regn,ndim);
|
|
fval=(*func)(pt);
|
|
for (j=1;j<=ndim;j++) {
|
|
if (pt[j]<=rmid[j]) {
|
|
fminl[j]=FMIN(fminl[j],fval);
|
|
fmaxl[j]=FMAX(fmaxl[j],fval);
|
|
}
|
|
else {
|
|
fminr[j]=FMIN(fminr[j],fval);
|
|
fmaxr[j]=FMAX(fmaxr[j],fval);
|
|
}
|
|
}
|
|
}
|
|
sumb=BIG;
|
|
jb=0;
|
|
siglb=sigrb=1.0;
|
|
for (j=1;j<=ndim;j++) {
|
|
if (fmaxl[j] > fminl[j] && fmaxr[j] > fminr[j]) {
|
|
sigl=FMAX(TINY,pow(fmaxl[j]-fminl[j],2.0/3.0));
|
|
sigr=FMAX(TINY,pow(fmaxr[j]-fminr[j],2.0/3.0));
|
|
sum=sigl+sigr;
|
|
if (sum<=sumb) {
|
|
sumb=sum;
|
|
jb=j;
|
|
siglb=sigl;
|
|
sigrb=sigr;
|
|
}
|
|
}
|
|
}
|
|
free_vector(fminr,1,ndim);
|
|
free_vector(fminl,1,ndim);
|
|
free_vector(fmaxr,1,ndim);
|
|
free_vector(fmaxl,1,ndim);
|
|
if (!jb) jb=1+(ndim*iran)/175000;
|
|
rgl=regn[jb];
|
|
rgm=rmid[jb];
|
|
rgr=regn[ndim+jb];
|
|
fracl=fabs((rgm-rgl)/(rgr-rgl));
|
|
nptl=(unsigned long)(MNPT+(npts-npre-2*MNPT)*fracl*siglb
|
|
/(fracl*siglb+(1.0-fracl)*sigrb));
|
|
nptr=npts-npre-nptl;
|
|
regn_temp=vector(1,2*ndim);
|
|
for (j=1;j<=ndim;j++) {
|
|
regn_temp[j]=regn[j];
|
|
regn_temp[ndim+j]=regn[ndim+j];
|
|
}
|
|
regn_temp[ndim+jb]=rmid[jb];
|
|
miser(func,regn_temp,ndim,nptl,dith,&avel,&varl);
|
|
regn_temp[jb]=rmid[jb];
|
|
regn_temp[ndim+jb]=regn[ndim+jb];
|
|
miser(func,regn_temp,ndim,nptr,dith,ave,var);
|
|
free_vector(regn_temp,1,2*ndim);
|
|
*ave=fracl*avel+(1-fracl)*(*ave);
|
|
*var=fracl*fracl*varl+(1-fracl)*(1-fracl)*(*var);
|
|
free_vector(rmid,1,ndim);
|
|
}
|
|
free_vector(pt,1,ndim);
|
|
}
|
|
#undef MNPT
|
|
#undef MNBS
|
|
#undef TINY
|
|
#undef BIG
|
|
#undef PFAC
|