#include <math.h>
#include "/misc/radon2/xiaoxiak/thesis/header/nr.h"
#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define TINY 1.0e-25
#define ITMAX 200
                
void powell(double *p, double **xi, int n, double ftol, int *iter, double *fret,
        double (*func)(double []))
{
        void linmin(double p[], double xi[], int n, double *fret,
                double (*func)(double []));
        int i,ibig,j;
        double del,fp,fptt,t,*pt,*ptt,*xit;

        pt=dvector(1,n);
        ptt=dvector(1,n);
        xit=dvector(1,n);
        *fret=(*func)(p);
        for (j=1;j<=n;j++) pt[j]=p[j];
        for (*iter=1;;++(*iter)) {
                fp=(*fret);
                ibig=0;
                del=0.0;
                for (i=1;i<=n;i++) {
                        for (j=1;j<=n;j++) xit[j]=xi[j][i];
                        fptt=(*fret);
                        linmin(p,xit,n,fret,func);
                        if (fptt-(*fret) > del) {
                                del=fptt-(*fret);
                                ibig=i;
                        }
                }
                if (2.0*(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))+TINY) 
{
                        free_dvector(xit,1,n);
                        free_dvector(ptt,1,n);
                        free_dvector(pt,1,n);
                        return;
                }
                if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");
                for (j=1;j<=n;j++) {
                        ptt[j]=2.0*p[j]-pt[j];
                        xit[j]=p[j]-pt[j];
                        pt[j]=p[j];
                }
                fptt=(*func)(ptt);
                if (fptt < fp) {
                        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
                        if (t < 0.0) {
                                linmin(p,xit,n,fret,func);
                                for (j=1;j<=n;j++) {
                                        xi[j][ibig]=xi[j][n];
                                        xi[j][n]=xit[j];
                                } 
                        }  
                }
        }
}
#undef ITMAX
#undef NRANSI
                        
                        
#include <math.h>
#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY2 1.0e-20
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
        
void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
        double (*func)(double))
{
        double ulim,u,r,q,fu,dum;
        
        *fa=(*func)(*ax);
        *fb=(*func)(*bx);
        if (*fb > *fa) {
                SHFT(dum,*ax,*bx,dum)
                SHFT(dum,*fb,*fa,dum)
        }
        *cx=(*bx)+GOLD*(*bx-*ax);
        *fc=(*func)(*cx);
        while (*fb > *fc) { 
                r=(*bx-*ax)*(*fb-*fc);
                q=(*bx-*cx)*(*fb-*fa);
                u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
                        (2.0*SIGN(FMAX(fabs(q-r),TINY2),q-r));
                ulim=(*bx)+GLIMIT*(*cx-*bx);
                if ((*bx-u)*(u-*cx) > 0.0) {
                        fu=(*func)(u);
                        if (fu < *fc) {
                                *ax=(*bx);
                                *bx=u;
                                *fa=(*fb);
                                *fb=fu;
                                return;
                        } else if (fu > *fb) {
                                *cx=u;
                                *fc=fu;
                                return;
                        }
                        u=(*cx)+GOLD*(*cx-*bx);
                        fu=(*func)(u);
                } else if ((*cx-u)*(u-ulim) > 0.0) {
                        fu=(*func)(u);
                        if (fu < *fc) {
                                SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
                                SHFT(*fb,*fc,fu,(*func)(u))
                        }
                } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
                        u=ulim;
                        fu=(*func)(u);
                } else {
                        u=(*cx)+GOLD*(*cx-*bx);
                        fu=(*func)(u);
                }
                SHFT(*ax,*bx,*cx,u)   
                SHFT(*fa,*fb,*fc,fu)  
        }
}
#undef GOLD
#undef GLIMIT
#undef TINY2
#undef SHFT
#undef NRANSI
                        

#define NRANSI
#include "nrutil.h"
#define TOL 2.0e-4

int ncom;
double *pcom,*xicom,(*nrfunc)(double []);
                        
void linmin(double p[], double xi[], int n, double *fret, double (*func)(double []))
{
        double brent(double ax, double bx, double cx,
                double (*f)(double), double tol, double *xmin);
        double f1dim(double x);
        void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
                double *fc, double (*func)(double));
        int j;
        double xx,xmin,fx,fb,fa,bx,ax;

        ncom=n;
        pcom=dvector(1,n);
        xicom=dvector(1,n);
        nrfunc=func;            
        for (j=1;j<=n;j++) {
                pcom[j]=p[j];
                xicom[j]=xi[j];
        }
        ax=0.0;
        xx=1.0;
        mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
        *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
        for (j=1;j<=n;j++) {
                xi[j] *= xmin;
                p[j] += xi[j];
        }
        free_dvector(xicom,1,n);
        free_dvector(pcom,1,n);
}
#undef TOL
#undef NRANSI

                        
#define NRANSI
#include "nrutil.h"
#define TOL 2.0e-4
                        
void (*nrdfun)(double [], double []);
                
void dlinmin(double p[], double xi[], int n, double *fret, double (*func)(double []),
        void (*dfunc)(double [], double []))
{
        double dbrent(double ax, double bx, double cx,
                double (*f)(double), double (*df)(double), double tol, double *xmin);
        double f1dim(double x);   
        double df1dim(double x); 
        void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
                double *fc, double (*func)(double));
        int j;
        double xx,xmin,fx,fb,fa,bx,ax;
                        
        ncom=n;
        pcom=dvector(1,n);
        xicom=dvector(1,n);
        nrfunc=func;
        nrdfun=dfunc;
        for (j=1;j<=n;j++) {
                pcom[j]=p[j];
                xicom[j]=xi[j];
        }
        ax=0.0; 
        xx=1.0; 
        mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
        *fret=dbrent(ax,xx,bx,f1dim,df1dim,TOL,&xmin);
        for (j=1;j<=n;j++) {
                xi[j] *= xmin;
                p[j] += xi[j];
        }               
        free_dvector(xicom,1,n);
        free_dvector(pcom,1,n);
}
#undef TOL      
#undef NRANSI
        
 
#define NRANSI   
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
        
extern int ncom;
extern double *pcom,*xicom,(*nrfunc)(double []);
        
double f1dim(double x)
{
        int j;
        double f,*xt;
                
        xt=dvector(1,ncom);
        for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
        f=(*nrfunc)(xt);
        free_dvector(xt,1,ncom);
        return f;
}
#undef NRANSI

#define NRANSI
extern void (*nrdfun)(double *,double *);

double df1dim(double x)
{
   int j;
   double df1=0.0;
   double *xt,*df;
   xt=dvector(1,ncom);
   df=dvector(1,ncom);
   for(j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
   (*nrdfun)(xt,df);
   for(j=1;j<=ncom;j++) df1 += df[j]*xicom[j];
   free_dvector(df,1,ncom);
   free_dvector(xt,1,ncom);
   return df1;
}
#undef NRANSI


#include <math.h>
#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define ITMAX 100
#define CGOLD 0.3819660
#define ZEPS 1.0e-10
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
double brent(double ax, double bx, double cx, double (*f)(double), double tol,
        double *xmin)
{                               
        int iter;
        double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
        double e=0.0;

        a=(ax < cx ? ax : cx);
        b=(ax > cx ? ax : cx);
        x=w=v=bx;   
        fw=fv=fx=(*f)(x);
        for (iter=1;iter<=ITMAX;iter++) {
                xm=0.5*(a+b);
                tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
                if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
                        *xmin=x;
                        return fx;
                }
                if (fabs(e) > tol1) {
                        r=(x-w)*(fx-fv);
                        q=(x-v)*(fx-fw);
                        p=(x-v)*q-(x-w)*r;
                        q=2.0*(q-r);
                        if (q > 0.0) p = -p;
                        q=fabs(q);
                        etemp=e;
                        e=d;
                        if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
                                d=CGOLD*(e=(x >= xm ? a-x : b-x));
                        else {
                                d=p/q;
                                u=x+d;
                                if (u-a < tol2 || b-u < tol2)
                                        d=SIGN(tol1,xm-x);
                        }
                } else {        
                        d=CGOLD*(e=(x >= xm ? a-x : b-x));
                }
                u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
                fu=(*f)(u);
                if (fu <= fx) {
                        if (u >= x) a=x; else b=x;
                        SHFT(v,w,x,u)
                        SHFT(fv,fw,fx,fu)
                } else {
                        if (u < x) a=u; else b=u;
                        if (fu <= fw || w == x) {
                                v=w;
                                w=u;
                                fv=fw;
                                fw=fu;
                        } else if (fu <= fv || v == x || v == w) {
                                v=u;
                                fv=fu;  
                        }
                }
        }
        nrerror("Too many iterations in brent");
        *xmin=x;
        return fx;
}
#undef ITMAX
#undef CGOLD
#undef ZEPS
#undef SHFT
#undef NRANSI
                                        
                         
#include <math.h>
#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define ITMAX 100
#define ZEPS 1.0e-10
#define MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f);
                        
double dbrent(double ax, double bx, double cx, double (*f)(double),
        double (*df)(double), double tol, double *xmin)
{
        int iter,ok1,ok2;
        double a,b,d,d1,d2,du,dv,dw,dx,e=0.0;
        double fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;
                
        a=(ax < cx ? ax : cx); 
        b=(ax > cx ? ax : cx);
        x=w=v=bx;
        fw=fv=fx=(*f)(x);
        dw=dv=dx=(*df)(x);
        for (iter=1;iter<=ITMAX;iter++) {
                xm=0.5*(a+b);
                tol1=tol*fabs(x)+ZEPS;
                tol2=2.0*tol1;
                if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
                        *xmin=x;
                        return fx;
                }
                if (fabs(e) > tol1) {
                        d1=2.0*(b-a);
                        d2=d1;
                        if (dw != dx) d1=(w-x)*dx/(dx-dw);
                        if (dv != dx) d2=(v-x)*dx/(dx-dv);
                        u1=x+d1;
                        u2=x+d2;
                        ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
                        ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
                        olde=e;
                        e=d;
                        if (ok1 || ok2) {
                                if (ok1 && ok2)
                                        d=(fabs(d1) < fabs(d2) ? d1 : d2);
                                else if (ok1)
                                        d=d1;
                                else
                                        d=d2;
                                if (fabs(d) <= fabs(0.5*olde)) {
                                        u=x+d;
                                        if (u-a < tol2 || b-u < tol2)
                                                d=SIGN(tol1,xm-x);
                                } else {
                                        d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
                                }
                        } else {
                                d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
                        }
                } else {
                        d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
                }
                if (fabs(d) >= tol1) {
                        u=x+d;
                        fu=(*f)(u);  
                } else {
                        u=x+SIGN(tol1,d);
                        fu=(*f)(u);
                        if (fu > fx) {
                                *xmin=x;
                                return fx;
                        }
                }
                du=(*df)(u);   
                if (fu <= fx) {
                        if (u >= x) a=x; else b=x;
                        MOV3(v,fv,dv, w,fw,dw) 
                        MOV3(w,fw,dw, x,fx,dx)
                        MOV3(x,fx,dx, u,fu,du)
                } else {
                        if (u < x) a=u; else b=u;
                        if (fu <= fw || w == x) {
                                MOV3(v,fv,dv, w,fw,dw)
                                MOV3(w,fw,dw, u,fu,du)
                        } else if (fu < fv || v == x || v == w) {
                                MOV3(v,fv,dv, u,fu,du)
                        }
                }
        }
        nrerror("Too many iterations in routine dbrent");
        return 0.0;
}
#undef ITMAX
#undef ZEPS
#undef MOV3
#undef NRANSI
                        
                        
#include <math.h>
#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define ITMAX 200
#define EPS 1.0e-10
#define FREEALL free_dvector(xi,1,n);free_dvector(h,1,n);free_dvector(g,1,n);
                
void frprmn(double p[], int n, double ftol, int *iter, double *fret,
        double (*func)(double []), void (*dfunc)(double [], double []))
{
        void linmin(double p[], double xi[], int n, double *fret,
                double (*func)(double []));
        int j,its;
        double gg,gam,fp,dgg;
        double *g,*h,*xi;
         
        g=dvector(1,n);
        h=dvector(1,n);
        xi=dvector(1,n);
        fp=(*func)(p);
        (*dfunc)(p,xi);
        for (j=1;j<=n;j++) {
                g[j] = -xi[j];
                xi[j]=h[j]=g[j];
        }
        for (its=1;its<=ITMAX;its++) {
                *iter=its;
                linmin(p,xi,n,fret,func);
                if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {
                        FREEALL
                        return;
                }
                fp= *fret;
                (*dfunc)(p,xi);
                dgg=gg=0.0;
                for (j=1;j<=n;j++) {
                        gg += g[j]*g[j];
                        dgg += (xi[j]+g[j])*xi[j];
                }
                if (gg == 0.0) {
                        FREEALL
                        return;
                }
                gam=dgg/gg;
                for (j=1;j<=n;j++) {
                        g[j] = -xi[j];
                        xi[j]=h[j]=g[j]+gam*h[j];
                }
        }
        nrerror("Too many iterations in frprmn");
}
#undef ITMAX
#undef EPS
#undef FREEALL
#undef NRANSI
        
        
#include <math.h>
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
#define GET_PSUM  for (n=1;n<=ndim;n++) {\
                for (sum=0.0,m=1;m<=mpts;m++) sum += p[m][n];\
                psum[n]=sum;}
long idum=-3856; 
double tt;
                
void amebsa(double **p, double y[], int ndim, double pb[], double *yb, double ftol,
        double (*funk)(double []), int *iter, double temptr)
{
        double amotsa(double **p, double y[], double psum[], int ndim, double pb[],
                double *yb, double (*funk)(double []), int ihi, double *yhi, double fac);
        float ran1(long *idum);  
        int i,ihi,ilo,j,m,n,mpts=ndim+1;
        double rtol,sum,swap,yhi,ylo,ynhi,ysave,yt,ytry,*psum;
                 
        psum=dvector(1,ndim);
        tt = -temptr;
        GET_PSUM
        for (;;) {
                ilo=1;
                ihi=2;   
                ynhi=ylo=y[1]+tt*log(ran1(&idum));
                yhi=y[2]+tt*log(ran1(&idum));
                if (ylo > yhi) {
                        ihi=1;
                        ilo=2;
                        ynhi=yhi;
                        yhi=ylo;
                        ylo=ynhi;
                }
                for (i=3;i<=mpts;i++) {
                        yt=y[i]+tt*log(ran1(&idum));
                        if (yt <= ylo) {
                                ilo=i;  
                                ylo=yt;
                        }
                        if (yt > yhi) {
                                ynhi=yhi;
                                ihi=i;
                                yhi=yt;
                        } else if (yt > ynhi) {
                                ynhi=yt;
                        }
                }
                rtol=2.0*fabs(yhi-ylo)/(fabs(yhi)+fabs(ylo));
                if (rtol < ftol || *iter < 0) {
                        swap=y[1];
                        y[1]=y[ilo];
                        y[ilo]=swap;
                        for (n=1;n<=ndim;n++) {
                                swap=p[1][n];
                                p[1][n]=p[ilo][n];
                                p[ilo][n]=swap;
                        }
                        break;
                }
                *iter -= 2;
                ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,-1.0);
                if (ytry <= ylo) {
                        ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,2.0);
                } else if (ytry >= ynhi) {
                        ysave=yhi;
                        ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,0.5);
                        if (ytry >= ysave) {
                                for (i=1;i<=mpts;i++) {
                                        if (i != ilo) {
                                                for (j=1;j<=ndim;j++) {
                                                        psum[j]=0.5*(p[i][j]+p[ilo][j]);
                                                        p[i][j]=psum[j];
                                                }
                                                y[i]=(*funk)(psum);
                                        }
                                }
                                *iter -= ndim; 
                                GET_PSUM
                        }
                } else ++(*iter);
        }
        free_dvector(psum,1,ndim);
}
#undef GET_PSUM
                        
                                
#include <math.h>
#define NRANSI
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
                 
extern long idum;
extern double tt;
                
double amotsa(double **p, double y[], double psum[], int ndim, double pb[],
        double *yb, double (*funk)(double []), int ihi, double *yhi, double fac)
{
        float ran1(long *idum);
        int j;
        double fac1,fac2,yflu,ytry,*ptry;
                                        
        ptry=dvector(1,ndim);
        fac1=(1.0-fac)/ndim;
        fac2=fac1-fac;
        for (j=1;j<=ndim;j++)
                ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
        ytry=(*funk)(ptry);
        if (ytry <= *yb) {
                for (j=1;j<=ndim;j++) pb[j]=ptry[j];
                *yb=ytry;
        }
        yflu=ytry-tt*log(ran1(&idum));
        if (yflu < *yhi) {
                y[ihi]=ytry;
                *yhi=yflu;
                for (j=1;j<=ndim;j++) {
                        psum[j] += ptry[j]-p[ihi][j];
                        p[ihi][j]=ptry[j];
                }               
        }                       
        free_dvector(ptry,1,ndim);
        return yflu;
}
#undef NRANSI