%
%  macm 202 -- 20 jan 03 -- djm
%
%  w4nm3.m:  given payoff matrix, find optimal "system"
%               probabilities by nelder-mead algorithm
%    - matrix "Pmatrix" must be previously defined
%    - to run, type "w4nm3" at the MATLAB prompt

%  set search length
Niter = 42;

%  simplex data (3 guesses as column vectors)
Asimp = eye(3,3);

Avals = zeros(1,3);
for j=1:3
    Avals(j) = min(transpose(Asimp(:,j))*Pmatrix);
end

%  start simplex loop (with plotting)
for j=1:Niter
 
%  simplex plotting (for mRPS only)
    Mat = [1/sqrt(3)  2/sqrt(3) ; 1 0];
    if (mod(j,3)==1)
        Aplot = (Mat*[Asimp(1:2,1:3) Asimp(1:2,1)])' ; 
        plot(Aplot(:,1),Aplot(:,2),'k')
        hold on;  axis equal;  drawnow;  pause(1)
    end

%  find best/worst
    [Abest Ibest] = max(Avals);
    [Awrst Iwrst] = min(Avals);

%  generate pivot point by AVERAGING ALL BUT WORST
    Apivot = zeros(3,1);
    for k=1:3
        Apivot = Apivot + (k~=Iwrst) * Asimp(:,k)/2;
    end
    
%  calculate four nelder-mead points 
    Anew1 =  2 *Apivot(:,1) -     Asimp(:,Iwrst); 
    Anew2 =  3 *Apivot(:,1) -  2 *Asimp(:,Iwrst); 
    Anew3 = 1.5*Apivot(:,1) - 0.5*Asimp(:,Iwrst); 
    Anew4 = 0.5*Apivot(:,1) + 0.5*Asimp(:,Iwrst); 

%  nullify nelder-mead points if not in [0,1]
    if (sum(Anew2>1)+sum(Anew2<0) ~= 0)
        Anew2 = Asimp(:,Iwrst);
    end
    if (sum(Anew1>1)+sum(Anew1<0) ~= 0)
        Anew1 = Asimp(:,Iwrst);
    end
    if (sum(Anew3>1)+sum(Anew3<0) ~= 0)
        Anew3 = Asimp(:,Iwrst);
    end

%  nelder-mead points & data
    Anew  = [Anew1 Anew2 Anew3 Anew4]; 

    Atest = zeros(1,4);
    for k=1:4
        Atest(k) = min(transpose(Anew(:,k))*Pmatrix);
    end

%  find nelder-mead best
    [Anmb Inew] = max(Atest);

%  replace if better, shrink if not
    if (Anmb>Awrst)
        Asimp(:,Iwrst) = Anew(:,Inew);
        Avals(  Iwrst) = Anmb;
        disp([j Inew Anmb-Awrst])
    else
        for k=1:3
            Asimp(:,k) = (Asimp(:,k) + Asimp(:,Ibest))/2;
            Avals(  k) = min(transpose(Asimp(:,k))*Pmatrix);
        end
        disp([j 0])
    end
end

Aprob = Asimp(:,Ibest)'
Abest = Abest