## Author: M. Turhan ÇOBAN
## Created: 2021-06-19
function retval = nelder_opt(f,a,da,maxiteration,tolerance) 
    ## reflection
    alpha=1.0;
    ## contraction
    beta=0.5;
    ## expansion
    gamma=2.0;
    ## shrinkage
    delta=0.5;
    ## change basic multiplication parameters if needed
    n=length(a);
    m=n+1;
    NDIMS = n;
    NPTS = n+1;
    FUNC = n+1;
    ncalls = 0;
    p=zeros(m,m);
    x=zeros(m,n);
    for i=1:m
        for j=1:n
            if i==j
                x(i,j)=a(j)+da(j);               
            else
                x(i,j)=a(j);
            endif
            p(i,j)=x(i,j);
        endfor   
    endfor
    for i=1:m
      p(i,m)=f(x(i,:));
    endfor
    #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 //////////////////
    best = 1.0E99;
    z=zeros(1,NDIMS);
    #/////////////// calculate the first function values for the simplex ////////////////
    iter=0;
    for iter=1:maxiteration
        #////////  define lo, nhi, hi (low high next_to_high //////////////
        ilo=1;
        ihi=1;
        inhi = 0 ; 
        # 0 means missing
        flo = p(1,FUNC);
        fhi = flo;
        for i=1:NPTS
            if p(i,FUNC) < flo 
                flo=p(i,FUNC); 
                ilo=i;
            endif  
            if  p(i,FUNC) > fhi
                fhi=p(i,FUNC) ;
                ihi=i;
            endif            
        endfor ## i=1:NPTS
        ##fprintf("fhi = %i  ihi = %i flo = %i ili= %i ",fhi,ihi,flo,ilo);
        fnhi = flo;
        inhi = ilo;
        for i=1:NPTS 
            if (i != ihi) && (p(i,FUNC) > fnhi) 
                fnhi=p(i,FUNC); 
                inhi=i;
            endif ##  (i != ihi) && (p(i,FUNC) > fnhi) 
        endfor ## i=1:NPTS 
        ##fprintf("fnhi= %i inhi = %i\n",fnhi,inhi)
                #///////// exit criteria //////////////
        xx1=mod(iter,4);
        xx2=xx1*NDIMS;
        bx= (xx2==0);
        if bx
          pavg=0.0;
          for i=1:NPTS
            pavg=pavg+p(i,FUNC);
          endfor
          pavg=pavg/NPTS;
          tot=0.0;
          for i=1:NPTS
            tot=(p(i,FUNC)-pavg)*(p(i,FUNC)-pavg);
          endfor
          sterr=sqrt(tot/NPTS);
          for j=1:NDIMS 
                z(j)=p(ilo,j);
          endfor
          best = p(ilo,FUNC);
        endif
        # end of if iter%4*NDMS)==0
        #//// calculate avarage without maximum value //////
        for j=1:NDIMS 
            ave(j) = 0;
        endfor ## j=1:NDIMS 
        for i=1:NPTS 
            if i != ihi
                for j=1:NDIMS 
                    ave(j) = ave(j)+p(i,j);
                endfor ## j=1:NDIMS 
            endif # i != ihi
        endfor ## i=1:NPTS 
        for j=1:NDIMS 
            ave(j)=ave(j)/(NPTS-1);
        end ##j=1:NDIMS 
        #//////// reflect ////////////////
        xihi=ihi;
        for j=1:NDIMS 
            r(j) = ave(j)+alpha*(ave(j) - p(ihi,j));
        endfor ## j=1:NDIMS 
        fr = f(r);
        if (flo <= fr) && (fr < fnhi) # in zone accept reflection
            ##printf("reflect");
            for j=1:NDIMS
                p(ihi,j) = r(j);
            endfor ## j=1:NDIMS
            p(ihi,FUNC) = fr;
            continue;
        endif ## flo <= fr) && (fr < fnhi)
        #///(//////////// expand //////////////
        if fr < flo 
            for j=1:NDIMS 
                e(j) = ave(j)+gamma*(r(j)-ave(j));
            endfor ## j=1:NDIMS  
            fe = f(e);
            if fe < fr 
                ##fprintf("expand fe < fr");
                for j=1:NDIMS 
                    p(ihi,j) = e(j);
                endfor
                p(ihi,FUNC) = fe;
                continue;
            else 
                ##fprintf("expand fe >= fr");
                for j=1:NDIMS 
                    p(ihi,j) = r(j);
                endfor ## j=1:NDIMS 
                p(ihi,FUNC) = fr;
                continue;
            endif ## fe < fr
        endif ## fr < flo  
        #/////////////// contract 

        if fr >= fnhi 
          if fr<fhi 
            for j=1:NDIMS 
              cc(j) = ave(j) + beta*(r(j)-ave(j));
            endfor ## j=1:NDIMS 
            fcc = f(cc);
            if fcc <= fr ##outside
                ##fprintf("contract outside");
                for j=1:NDIMS 
                    p(ihi,j) = cc(j);
                endfor ## j=1:NDIMS 
                p(ihi,FUNC) = fcc;
                continue;
            endif ## fcc <= fr
          else ## fr>=fhi
            for j=1:NDIMS 
              cc(j) = ave(j) + beta*(p(ihi,j)-ave(j));
            endfor ## j=1:NDIMS 
            fcc = f(cc);
            if fcc<fhi
                ##fprintf("contract outside");
                for j=1:NDIMS 
                    p(ihi,j) = cc(j);
                endfor ## j=1:NDIMS 
                p(ihi,FUNC) = fcc;
                continue;
            endif ## fcc<fhi
          endif ## if fr<fhi else fr>=fhi 
         endif ##  fr >= fnhi 
        #//////////// shrink
         ##fprintf("shrink");
                for i=1:NPTS 
                    if i != ilo 
                        for j=1:NDIMS 
                            p(i,j) = p(i,j)+delta*(p(i,j)-p(ilo,j));
                        endfor
                        p(i,FUNC) = f(p(i,:));
                    endif
                continue;
                endfor    
   endfor  ## iter=1:maxiteration
   retval =z;
endfunction

