97 static const realnum bnsfac[2][3] = { {-1.0f,-2.0f,0.0f}, {1.0f,0.0f,2.0f} };
195 if( scale[0] >= 0.0f )
197 for( i=1; i <= n; i++ )
200 xpscl = x[i_] + scale[i_];
207 scl[0] = (
realnum)fabs(scale[0]);
208 for( i=1; i <= n; i++ )
211 xpscl = x[i_] + scl[0];
216 if( ((mode/2)%2) == 0 )
220 else if(
usubc.alpha <= 0.0f )
224 else if(
usubc.beta <= 0.0f ||
usubc.beta >= 1.0f )
228 else if(
usubc.gamma <= 1.0f )
232 else if(
usubc.delta <= 0.0f ||
usubc.delta >= 1.0f )
236 else if(
usubc.psi <= 0.0f ||
usubc.psi >= 1.0f )
240 else if(
usubc.omega <= 0.0f ||
usubc.omega >= 1.0f )
249 else if( n < ((n - 1)/
usubc.nsmax + 1)*
usubc.nsmin )
261 else if(
usubc.bonus < 0.0f )
265 else if(
usubc.nfstop < 0 )
274 ifsptr = isptr +
usubc.nsmax*(
usubc.nsmax + 3);
276 if( scale[0] > 0.0f )
279 cdcopy(n,scale,1,&work[istptr-1],1);
284 cdcopy(n,scl,0,&work[istptr-1],1);
286 for( i=1; i <= n; i++ )
293 if(
usubc.irepl == 0 )
297 else if(
usubc.minf )
305 if(
usubc.nfstop == 0 )
309 else if(
usubc.minf )
328 else if( *iflag == 2 )
341 else if( *iflag == -1 )
346 else if( *iflag == 0 )
360 for( i=0; i < n; i++ )
362 work[i] = (
realnum)fabs(work[i]);
366 partx(n,iwork,work,&nsubs,&iwork[insptr-1]);
369 insfnl = insptr + nsubs - 1;
379 simplx(n,&work[istptr-1],ns,&iwork[ipptr-1],maxnfe,cmode,x,&sfx,
380 nfe,&work[isptr-1],&work[ifsptr-1],iflag);
395 for( i=0; i < n; i++ )
397 work[i] = x[i] - work[i];
409 for( i=0; i < n; i++ )
412 xstop +=
POW2(work[i]);
413 ttemp = (
realnum)fabs(work[istep]);
416 xstop2 =
MAX2( xstop2 , ttemp );
420 if( sqrt(xstop) > tol || xstop2*
usubc.psi > (
realnum)tol )
422 setstp(nsubs,n,work,&work[istptr-1]);
643 if( incx == 1 && incy == 1 )
652 for( i=0; i < m; i++ )
662 for( i=m; i < n; i += 7 )
682 ix = (-n + 1)*incx + 1;
684 iy = (-n + 1)*incy + 1;
685 for( i=0; i < n; i++ )
762 for( i=1; i <= ns; i++ )
765 x[ips[i_]-1] = xs[i_];
770 if(
usubc.irepl == 0 )
781 else if(
isubc.IntNew )
786 newbst = fx <
usubc.ftest;
791 newbst = fx >
usubc.ftest;
793 if(
usubc.initx || newbst )
795 if(
usubc.irepl == 1 )
803 if(
usubc.irepl == 1 )
892 for( i=1; i <= n; i++ )
895 if( deltax[i_] == 0.f )
897 step[i_] = -step[i_];
956 while( ifirst <= ilast )
958 for( i=ifirst; i <= ilast; i++ )
963 if( xkey[ixi-1] < xkey[ixip1-1] )
971 for( i=ilast; i >= ifirst; i-- )
976 if( xkey[ixi-1] < xkey[ixip1-1] )
1052 for( i=1; i < n; i++ )
1061 for( i=0; i < (
usubc.nsmin - 1); i++ )
1063 as1 += absdx[ip[nused+i]-1];
1068 for( ns1=
usubc.nsmin; ns1 <= limit; ns1++ )
1071 as1 += absdx[ip[nused+ns1_]-1];
1075 if( ns2 >= ((ns2 - 1)/
usubc.nsmax + 1)*
usubc.nsmin )
1078 gap = as1/ns1 - as2/ns2;
1082 nsvals[*nsubs-1] = ns1;
1087 else if( as1/ns1 > gapmax )
1092 nused += nsvals[*nsubs-1];
1098 nsvals[*nsubs-1] = ns1;
1105#define S(I_,J_) (*(s+(I_)*(ns)+(J_)))
1220 start(n,x,step,ns,ips,s,&small);
1228 if(
usubc.irepl > 0 )
1230 isubc.IntNew =
false;
1231 evalf(ns,ips,&
S(0,0),n,x,&fs[0],nfe);
1237 isubc.IntNew =
true;
1238 for( j=2; j <= npts; j++ )
1241 evalf(ns,ips,&
S(j_,0),n,x,&fs[j_],nfe);
1244 order(npts,fs,&il,&is,&ih);
1251 calcc(ns,s,ih,inew,updatc,&
S(icent-1,0));
1257 newpt(ns,
usubc.alpha,&
S(icent-1,0),&
S(ih-1,0),
true,&
S(itemp-1,0),
1261 evalf(ns,ips,&
S(itemp-1,0),n,x,&fr,nfe);
1267 newpt(ns,-
usubc.gamma,&
S(icent-1,0),&
S(itemp-1,0),
true,
1271 evalf(ns,ips,&
S(ih-1,0),n,x,&
fe,nfe);
1278 cdcopy(ns,&
S(itemp-1,0),1,&
S(ih-1,0),1);
1282 else if( fr < fs[is-1] )
1287 cdcopy(ns,&
S(itemp-1,0),1,&
S(ih-1,0),1);
1298 &
S(itemp-1,0),&small);
1302 newpt(ns,-
usubc.beta,&
S(icent-1,0),&
S(itemp-1,0),
false,
1307 evalf(ns,ips,&
S(itemp-1,0),n,x,&fc,nfe);
1310 cdcopy(ns,&
S(itemp-1,0),1,&
S(ih-1,0),1);
1318 for( j=1; j <= npts; j++ )
1327 evalf(ns,ips,&
S(j_,0),n,x,&fs[j_],nfe);
1333 order(npts,fs,&il,&is,&ih);
1340 if(
usubc.irepl == 0 )
1353 if( *nfe >= maxnfe )
1355 if( !(
dist(ns,&
S(ih-1,0),&
S(il-1,0)) <= tol || small) )
1371 for( i=0; i < ns; i++ )
1373 x[ips[i]-1] =
S(il-1,i);
1383 static long int nsv;
1425 usubc.fxstat[3] = 0.0f;
1430 f1sv =
usubc.fxstat[0];
1431 usubc.nfxe += ifxwt;
1437 fscale)+nsv*
POW2((
usubc.fxstat[0]-f1sv)/fscale)+ifxwt*
1478 for( i=0; i < m; i++ )
1480 dtemp += (
realnum)fabs(dx[i]);
1486 for( i=m; i < n; i += 6 )
1488 dtemp += (
realnum)(fabs(dx[i]) + fabs(dx[i+1]) + fabs(dx[i+2]) +
1489 fabs(dx[i+3]) + fabs(dx[i+4]) + fabs(dx[i+5]));
1502 ix = (-n + 1)*incx + 1;
1504 for( i=0; i < n; i++ )
1506 dtemp += (
realnum)fabs(dx[ix-1]);
1537 if( !(n <= 0 || incx <= 0) )
1550 for( i=1; i <= m; i++ )
1561 for( i=mp1; i <= n; i += 5 )
1579 for( i=0; i<nincx; i=i+incx)
1592#define S(I_,J_) (*(s+(I_)*(ns)+(J_)))
1650 for( i=1; i <= ns; i++ )
1653 S(0,i_) = x[ips[i_]-1];
1656 for( j=2; j <= (ns + 1); j++ )
1660 S(j_,j_-1) =
S(0,j_-1) + step[ips[j_-1]-1];
1665 for( j=2; j <= (ns + 1); j++ )
1668 if( (
double)(
S(j_,j_-1)) == (
double)(
S(0,j_-1)) )
1731 if( fs[j-1] >= fs[*il-1] )
1743 for( i=il0 + 1; i <= (il0 + npts - 2); i++ )
1746 if( fs[j-1] >= fs[*ih-1] )
1751 else if( fs[j-1] > fs[*is-1] )
1755 else if( fs[j-1] < fs[*il-1] )
1802 absxmy = (
realnum)fabs(x[0]-y[0]);
1803 if( absxmy <= 1.0f )
1805 sum = absxmy*absxmy;
1814 for( i=2; i <= n; i++ )
1817 absxmy = (
realnum)fabs(x[i_]-y[i_]);
1818 if( absxmy <= scale )
1820 sum +=
POW2(absxmy/scale);
1824 sum = 1.0f + sum*
POW2(scale/absxmy);
1828 dist_v = (
realnum)(scale*sqrt(sum));
1835#define S(I_,J_) (*(s+(I_)*(ns)+(J_)))
1848 realnum xNothing[1] = { 0.0f };
1896 cdcopy(ns,xNothing,0,c,1);
1897 for( j=1; j <= (ns + 1); j++ )
1901 cdaxpy(ns,1.0f,&
S(j_,0),1,c,1);
1905 else if( ih != inew )
1907 for( i=0; i < ns; i++ )
1909 c[i] += (
S(inew-1,i) -
S(ih-1,i))/ns;
1990 for( i=0; i < ns; i++ )
1992 xnew[i] = (
realnum)(xbase[i] + coef*(xbase[i] - xold[i]));
1993 eqbase = eqbase && ((double)(xnew[i]) == (double)(xbase[i]));
1994 eqold = eqold && ((double)(xnew[i]) == (double)(xold[i]));
1999 for( i=0; i < ns; i++ )
2002 xold[i] = (
realnum)(xbase[i] + coef*(xbase[i] - xold[i]));
2003 eqbase = eqbase && ((double)(xold[i]) == (double)(xbase[i]));
2004 eqold = eqold && ((double)(xold[i]) == (double)(xoldi));
2007 *small = eqbase || eqold;
2035 if( incx == 1 && incy == 1 )
2046 for( i=1; i <= m; i++ )
2057 for( i=m; i < n; i += 4 )
2074 ix = (-n + 1)*incx + 1;
2076 iy = (-n + 1)*incy + 1;
2077 for( i=0; i < n; i++ )
STATIC double da(double z, double temp, double eden)
bool fp_equal(sys_float x, sys_float y, int n=3)
#define DEBUG_ENTRY(funcname)
chi2_type optimize_func(const realnum param[], int grid_index=-1)
STATIC double dist(long, realnum[], realnum[])
STATIC void sortd(long, realnum[], long[])
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
STATIC void cdcopy(long, realnum[], long, realnum[], long)
STATIC void cdaxpy(long, double, realnum[], long, realnum[], long)
STATIC void newpt(long, double, realnum[], realnum[], int, realnum[], int *)
STATIC void csscal(long, double, realnum[], long)
STATIC void evalf(long, long[], realnum[], long, realnum[], realnum *, long *)
STATIC void fstats(double, long, int)
STATIC void simplx(long, realnum[], long, long[], long, int, realnum[], realnum *, long *, realnum *, realnum[], long *)
STATIC void calcc(long, realnum *, long, long, int, realnum[])
void optimize_subplex(long int n, double tol, long int maxnfe, long int mode, realnum scale[], realnum x[], realnum *fx, long int *nfe, realnum work[], long int iwork[], long int *iflag)
STATIC void setstp(long, long, realnum[], realnum[])
STATIC double cdasum(long, realnum[], long)
STATIC void partx(long, long[], realnum[], long *, long[])
STATIC void order(long, realnum[], long *, long *, long *)