function [p,y,nfunk] = amoeba(funk, p, y, ndim, ftol, NMAX)
%function [p,y,nfunk] = amoeba(funk, p, y, ndim, ftol, NMAX)
% The Downhill Simplex method of optimization minimizes a scalar objective
% or cost function of ndim independent variables or parameters using 
% function evaluation only. It does not require derivatives of 
% the function of minimization. 
%
% INPUT
% funk - the name of the function that you want to minimize
% p - the N+1-by-N matrix, row stands for the vertices, column stands for
% the variables
% y - the function values of N+1 vertices
% ndim - the number of independentvariables
% ftol - functional convergence tolerance
% NMAX - the maximum number of function evaluations
% 
% OUTPUT
%   p,y has the same meaning
%   nfunk - the number of function evaluations
%
nfunk = 0;
TINY = 1e-20;
mpts = ndim + 1;
psum = sum(p,1);
while (1==1)
    ilo = 1;
    if y(1)>y(2)
        inhi = 2;
        ihi  = 1;
    else
        inhi = 1;
        ihi  = 2;
    end
    for i=1:mpts
        if (y(i) <= y(ilo)) 
            ilo = i;
        end
        if (y(i) > y(ihi)) 
            inhi = ihi;
            ihi  = i;
        elseif (y(i) > y(inhi)) && (i ~= ihi)
            inhi = i;
        end
    end
    
    rtol = 2.0 * abs(y(ihi)-y(ilo))/(abs(y(ihi)) + abs(y(ilo))+TINY);
    if (rtol < eps)
        [y(1),y(ilo)] = swap(y(1),y(ilo));
        [p(1,:),p(ilo,:)] = swap(p(1,:),p(ilo,:));
        break;
    end
    if (nfunk >= NMAX) 
        break;
    end
    nfunk = nfunk + 2;
    [ytry, p, y, psum] = amotry(funk,p,y,psum,ndim,ihi,-1.0);
    if (ytry <= y(ilo))
        [ytry, p, y, psum] = amotry(funk,p,y,psum,ndim,ihi,2.0);
    elseif (ytry >= y(inhi))
        ysave=y(ihi);
        [ytry, p, y, psum] = amotry(funk,p,y,psum,ndim,ihi,0.5);
        if (ytry >= ysave)
            for i=1:mpts
                if (i ~= ilo)
                    psum=0.5*(p(i,:) + p(ilo,:));
                    p(i,:)=psum;
                    y(i)=feval(@funk,psum);
                end
            end
            nfunk = nfunk + ndim;
            psum = sum(p,1);
        end
    else
        nfunk = nfunk - 2;
    end
end

function [ytry, p, y, psum] = amotry(funk, p, y, psum, ndim,ihi, fac)
fac1 = (1.0-fac) / ndim;
fac2 = fac1 - fac;
ptry = psum .* fac1 - p(ihi,:) .* fac2;
ytry = feval(funk,ptry);
if (ytry < y(ihi))
    y(ihi) = ytry;
    psum = psum + ptry - p(ihi,:);
    p(ihi,1:ndim)= ptry;
end

function [a,b]=swap(a,b)
swapv=a;
a=b;
b=swapv;