## Author: M. Turhan ÇOBAN
## Created: 2021-06-06
function retval = nelder_opt2(f,a,da,maxiteration,tolerance) 
 def  nelder(fnelder,a,da,maxiteration,tolerance,printlist): 
    #print("a",a)
    #print("da",da)
    n=len(a)
    m=n+1
    x=[[0.0 for j in range(n)] for i in range(m)]
    p=[[0.0 for j in range(m)] for i in range(m)]
    s=""
    NDIMS = len(x)-1
    NPTS = len(x)
    FUNC = NDIMS
    ncalls = 0;
    for i in range(m):
        for j in range(n):
            if i==j:
                x[i][j]=a[j]+da[j]
                p[i][j]=x[i][j]
            else:
                x[i][j]=a[j]
                p[i][j]=x[i][j]
        p[i][FUNC] = fnelder.func(p[i])
        #print("x[",i,"]=",x[i])    
        #print("p[",i,"]=",p[i])    
    print("p",p)
    #Inlet variable definitions
    #fnelder : abstract multivariable function f(x)
    #x : independent variable set of n+1 simplex elements
    #maxiteration : maximum iteration number
    #tolerance : 

    #///// construct the starting simplex //////////////////
    z=[0.0 for i in range(NDIMS)]
    best = 1.0E99;
    #/////////////// calculate the first function values for the simplex ////////////////
    iter=0;
    for iter in range(1,maxiteration): 
        #////////  define lo, nhi, hi (low high next_to_high //////////////
        ilo=0
        ihi=0
        inhi = -1  
        # -1 means missing
        flo = p[0][FUNC];
        fhi = flo
        #print("flo=",flo,"fhi=",fhi)
        #double pavg,sterr;
        for i in range(1,NPTS): 
            if p[i][FUNC] < flo: 
                flo=p[i][FUNC] 
                ilo=i  
            if  p[i][FUNC] > fhi:
                fhi=p[i][FUNC] 
                ihi=i
        print("flo=",flo,"fhi=",fhi,"ilo=",ilo,"ihi=",ihi)
        fnhi = flo
        inhi = ilo

        for i in range(NPTS): 
            if (i != ihi) and (p[i][FUNC] > fnhi): 
                fnhi=p[i][FUNC] 
                inhi=i
        print("fnhi=",fnhi,"inhi=",inhi)
        #///////// exit criteria //////////////
        xx=iter % 4*NDIMS
        print("xx=",xx,"iter=",iter,"NDIMS=",NDIMS)
        if (iter % 4*NDIMS) == 0:
            #calculate the avarage (including maximum value)
            pavg=0
            for i in range(NPTS): 
                 pavg=pavg+p[i][FUNC]
            pavg=pavg/NPTS
            print("pavg",pavg)
            tot=0.0
            if printlist!=0:
                print(iter)
                for j in range(NDIMS): 
                    print(p[ilo][j]," ")
            for i in range(NPTS): 
                tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg)
            sterr=sqrt(tot/NPTS)
            print("sterr",sterr)
            for j in range(NDIMS): 
                z[j]=p[ilo][j]
            best = p[ilo][FUNC];
            print("best",best)
        # end of if iter%4*NDMS)==0
        #//// calculate avarage without maximum value //////
        ave=[0.0 for i in range(NDIMS)]
        for j in range(NDIMS): 
            ave[j] = 0
        for i in range(NPTS): 
            if i != ihi:
                for j in range(NDIMS): 
                    ave[j] = ave[j]+p[i][j];
        for j in range(NDIMS): 
            ave[j]=ave[j]/(NPTS-1)
        #print("ave[",j,"]=",ave[j])
        #//////// reflect ////////////////
        #print("reflect")
        r=[0.0 for i in range(NDIMS)]
        for j in range(NDIMS): 
            r[j] = 2.0*ave[j] - p[ihi][j]
        fr = fnelder.func(r)
        if (flo <= fr) and (fr < fnhi): #in zone: accept
            for j in range(NDIMS):
                p[ihi][j] = r[j]
            p[ihi][FUNC] = fr
            continue
        #///(//////////// expand
        if fr < flo: 
            #print("expand")
            e=[0.0 for i in range(NDIMS)]
            for j in range(NDIMS): 
                e[j] = 3.0*ave[j] - 2.0*p[ihi][j];
            fe = fnelder.func(e);
            if fe < fr: 
                for j in range(NDIMS): 
                    p[ihi][j] = e[j]
                p[ihi][FUNC] = fe
                continue
            else: 
                for j in range(NDIMS): 
                    p[ihi][j] = r[j]
                p[ihi][FUNC] = fr
                continue
        #//////////// shrink:
        if fr < fhi:
            #print("srink")
            c=[0.0 for i in range(NDIMS)]
            for j in range(NDIMS): 
                c[j] = 1.5*ave[j] - 0.5*p[ihi][j]
            fc = fnelder.func(c);
            if fc <= fr: 
                for j in range(NDIMS): 
                    p[ihi][j] = c[j];
                p[ihi][FUNC] = fc
                continue
            else:   #////// shrink
                for i in range(NPTS): 
                    if i != ilo: 
                        for j in range(NDIMS): 
                            p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j]
                        p[i][FUNC] = fnelder.func(p[i])
                continue
        #/////////////// contract 
        if fr >= fhi:  
            #print("contract")
            cc=[0.0 for i in range(NDIMS)]
            for j in range(NDIMS): 
                cc[j] = 0.5*ave[j] + 0.5*p[ihi][j];
            fcc = fnelder.func(cc)
            if fcc < fhi:
                for j in range(NDIMS): 
                    p[ihi][j] = cc[j];
                p[ihi][FUNC] = fcc
                continue
            else:
                for i in range(NPTS): 
                    if i != ilo:
                        for j in range(NDIMS): 
                            p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j]
                        p[i][FUNC] = fnelder.func(p[i]);
    return z  
   retval =z;
endfunction