function da = demo_amoeba()
    
    x = zeros(3,2);
    y = zeros(3,1);

    alpha = 1.0; % expansion during reflection
    beta  = 2.0; % expansion
    gamma = 0.5; % contraction and multiple contraction

    niter = 30;
    
    x = [  0.1  -15.2;   ...
           4.7   19  ;   ...
          14.8    5.9 ];
    
    for i=1:3
        y(i) = f( x(i,1),x(i,2) );
    end;

    fileNAME = 'amoeba_points.dat';
    fileID = fopen(fileNAME,'w');

    %% find the minimum
    for k=1:niter
        
        if (y(3)<y(2)) 
            [x,y] = swap(x,y,2,3);
        end
        
        if (y(2)<y(1)) 
            [x,y] = swap(x,y,1,2);
            if (y(3)<y(2)) 
                [x,y] = swap(x,y,2,3);
            end;
        end;

        fprintf(fileID,' %+8.4e  %+8.4e  %+8.4e  %+8.4e  %+8.4e  %+8.4e  %+8.4e\n', ...
                x(1,:), x(2,:), x(3,:), y(1));

        if (k==niter)
            fprintf('\n After %4d iterations:\n',k);
            fprintf(' y: ( %8.4e , %8.4e , %8.4e )\n',y(:));
            fprintf(' x: [ %+8.4e , %+8.4e  \n',x(1,:));
            fprintf('      %+8.4e , %+8.4e  \n',x(2,:));
            fprintf('      %+8.4e , %+8.4e ]\n',x(3,:));
            %pause;
        end;
            
	for i=1:2
            c(i)  = sum(x(1:2,i)) / 2;
            xt(i) = c(i) + alpha*(c(i)-x(3,i));
        end; 

        yt = f( xt(1),xt(2) );
        
        if (yt<y(3)) 
            x(3,1:2) = xt(1:2);            
            y(3) = yt;
            if (yt<y(1))
                xe(1:2) = c(1:2) + beta*(xt(1:2)-c(1:2));
                ye = f( xe(1),xe(2) );
                if (ye<yt) 
                    x(3,1:2) = xe(1:2);
                    y(3) = ye;
                end;
            end;
        end;

        if (yt>y(2))
            xc(1:2) = c(1:2) + gamma*(x(3,1:2)-c(1:2));
            yc = f( xc(1),xc(2) );
            if (yc<y(3))
                x(3,1:2) = xc(1:2);
                y(3)=yc;
            else
                for i=2:3
                    x(i,1:2) = x(1,1:2) + gamma*(x(i,1:2)-x(1,1:2));
                end;
                y(2) = f( x(2,1),x(2,2) );
                y(3) = f( x(3,1),x(3,2) );
            end;
        end;
    end;

    fclose(fileID);

    %% plot the result
    pts = load(fileNAME);
    [niter bf] = size(pts);
    xmin = min([ min(pts(:,1)) min(pts(:,3)) min(pts(:,5)) ]);
    xmax = max([ max(pts(:,1)) max(pts(:,3)) max(pts(:,5)) ]);
    ymin = min([ min(pts(:,2)) min(pts(:,4)) min(pts(:,6)) ]);
    ymax = max([ max(pts(:,2)) max(pts(:,4)) max(pts(:,6)) ]);

    fig=figure;
    for k=1:niter
        
        disp(['Plotting iter ' num2str(k)]);

        subplot(1,2,1);
        fill(pts(k,1:2:5),pts(k,2:2:6),'r','edgecolor','r');
       
        subplot(1,2,2);
        fill(pts(k,1:2:5),pts(k,2:2:6),'r','edgecolor','r','linewidth',3);
        axis([xmin xmax ymin ymax]);

        set(gcf,'NextPlot','add'); axes; h=title('MyTitle');
        set(gca,'Visible','off');  set(h,'Visible','on'); 
        set(gca,'Fontsize',14)
        title(['Iter.' int2str(k) ':  F (' num2str(pts(k,1)) ' , ' num2str(pts(k,2)) ' ) = ' num2str(pts(k,7)) ])
        
        pause
    end
    
    
    return;
    
%%%%%%%%%%%%

function fun = f(x,y)
    
    fun = 100 * (y-x^2)^2 + (1-x)^2;
    % fun = x^2+y^2;
        
    return;

%%%%%%%%%%%%

function [x,y] = swap(x,y,i1,i2)
    
    c       = x(i1,:);
    x(i1,:) = x(i2,:);
    x(i2,:) = c(:);
    
    d     = y(i1);
    y(i1) = y(i2);
    y(i2) = d;
    
    return;

%%%%%%%%%%%%
