function [pi,t]=fpfsp(N,rho,alpha,traj,n,m,t,z)
% [pi,t]=fpfsp(N,rho,alpha,traj,n,m,t,z)
% "fast" permutation flowshop problem
% solved via CE with "fast" traj. gen.
% Outputs -
% pi: the output permutation
% t:  the cost (in time) matrix used
% Inputs -
% N: Number of samples to generate each round
% rho: fraction of best samples to take
% alpha: smoothing parameter
% traj: 0 for alias technique (default) uses x% = 70%
%       1 for composition method
% n: number of jobs
% m: number of machines
% t: t(i,j) is the cost of job i on mach j
% z: number of succesive rho-th quantiles the same before stop

xx=.7;

     if nargin<=4,
	       z=10; % additional stopping criterion
	       n=5; % number of jobs
     m=3; % number of machines;
t=[1,.3,.1;.4,.2,.6;.1,.5,.8;.3,.8,.9;.1,.2,.3];
end
if nargin<=3,
	  traj=0;          
end
if nargin==0,
	       N=100;
	       rho=.05;
alpha=.5;
end
tol=0.005;

kk=(1:n)';

% initial P

st=zeros(z,1);

Pold=zeros(n);
if (traj==1),
P=ones(n)*1/(n-1);
P=P-diag(diag(P));
else
P=ones(n)*1/n;
end

count=1;
while max(max(abs(P-Pold)))>tol,
%X=[];
%S=[];
for i=1:N,
if traj==1,
	X(i,:)=fgpfsp2(P);
else
X(i,:)=fgpfsp(P,xx);
end
if(traj==1),
S(i)=spfsp(X(i,:),t);
else

y=X(i,:);

for zi=1:n,
        ki=y(zi);
x(ki)=zi;
end
%size(x)
%size(Y)
Y(i,:)=x;

%end

S(i)=spfsp(x,t);

end

end
[S,I]=sort(S);
g=floor(rho*N);

if(traj==1),

for i=1:n,
	for j=1:n,
% portion of elite samples having a transition from i to j
                ss=0;
for k=1:g,
                ii=kk(X(I(k),:)==i);
                if(ii<n),
                    if(X(I(k),ii+1)==j),
ss=ss+1;
end
else
if(X(I(k),1)==j),
ss=ss+1;
end
end
end

		w(i,j)=ss/g;
end
end

else %traj==0, ie node placement

for i=1:n,
for j=1:n,

%w(i,j)=sum(S(I(1:g)).*(Y(:,i)==j))/sum(S(I)
ss=0;
for k=1:g,
  ss=ss+(Y(I(k),i)==j);
end

w(i,j)=ss/g;

end
end

end

Pold=P;
P=alpha*w+(1-alpha)*P;
st(1:(z-1))=st(2:z);
st(z)=S(g);
[count,S(g),S(1)]
if (sum(st==st(z))==z),
break;
end
count=count+1;
end
if(traj==1),
pi=X(I(1),:);
else
pi=Y(I(1),:);
end
