//
// Genetic Algorithm with an approximated evaluation principle 
// Only a decreasing percent of the population is evaluated exactly
// The approximation is build with a local RBF interpolation
//

function y=exactfunction(x)                //    the function to optimize
n=length(x);
y=n+sum(x.^2-cos(2*%pi*x));               // Rastrigin                           
endfunction
//
//----------------------------------------------------------------
//
function y=approxfunction(x,Xexact,A,b)  // the approximated function by RBF interpolation
[N1,N2]=size(Xexact);n=N2-1;
u=linsolve(A,-b);
y=0;
for i=1:N1 
   y=y+u(i)*phi(x-Xexact(i,1:n))          
end
endfunction
//
//----------------------------------------------------------------
//
function y=phi(x)                         // RBF function
y=norm(x)
endfunction
//
//----------------------------------------------------------------
//
function [A,b]=construct(Xexact)        // construction of the matrix for interpolation
//
[N1,N2]=size(Xexact);n=N2-1;
A=zeros(N1,N1);b=zeros(N1,1);
for i=1:N1
  for j=1:N1
    A(i,j)=phi(Xexact(i,1:n)-Xexact(j,1:n));
  end
  b(i)=Xexact(i,n+1);
end
endfunction
//
//---------------------------------------------------------------------
function Xexactsel=selection(x,Xexact,Nsize)  // selection of the centers for interpolation
//
[N1,N2]=size(Xexact);n=N2-1;normx=[];
for i=1:N1
  normx(i)=norm(x-Xexact(i,1:n))
end
[x1,x2]=gsort(-normx);
Xexactsel=Xexact(x2(1:Nsize),:)
endfunction
//----------------------------------------------------------------
function evalpop=evaluation(X);              // exact evaluation of the whole population
[Npop,N2]=size(X);n=N2-1;
evalpop=zeros(Npop,1);
for i=1:Npop;
  evalpop(i)=exactfunction(X(i,1:n));
end
endfunction
//
//----------------------------------------------------------------
//
function evalpop=evaluationapprox(X,Xexact,Nsize);  // approximated evaluation of the whole population
[Npop,N2]=size(X);n=N2-1;
evalpop=zeros(Npop,1);
for i=1:Npop;
  Xexactsel=selection(X(i,1:n),Xexact,Nsize);
  [A,b]=construct(Xexactsel);
  evalpop(i)=approxfunction(X(i,1:n),Xexactsel,A,b);
end
endfunction
//
//------------------------------------------------------------------
//
function [Xpar,bestpar]=elitism(Xnew,Xpar,Xtot,Npop,gen) // Elitism
[minval,index]=min(Xtot(:,n+1));  
bestpar=Xtot(index(1),:);
[minval,index]=min(Xnew(:,n+1));
bestnewpop=Xnew(index(1),:);
if (bestnewpop($)>bestpar($))&(gen>1) then
    nrand=int(Npop*rand())+1;
    Xpar(nrand,:)=bestpar;
else 
    bestpar=bestnewpop;
end
endfunction
//
//--------------------------------------------------------------
//
function Xsel=rankselect(X)                           // ranking selection
Xsel=X;
[Npop,n]=size(X);n=n-1;
[Fgsort,ind]=gsort(X(:,n+1));
X=X(ind,:);
p=1:Npop;roulette=cumsum(p)/sum(p);
for i=1:Npop;
  u=rand();
  index=find(roulette<u);Nselect=max(index)+1;
  Xsel(i,:)=X(Nselect,:);
end
endfunction
//
//--------------------------------------------------------------
//
function xoff=blcross(xpar)                             // blend crossover
[np,n]=size(xpar);n=n-1;
xoff=xpar;
for i=1:n
   u=rand();
   xoff(1,:)=u*xpar(1,:)+(1-u)*xpar(2,:);
   xoff(2,:)=u*xpar(2,:)+(1-u)*xpar(1,:);
end
endfunction
//
//--------------------------------------------------------------
//
function xoff=mutate(xpar,ngen,Ngen,xmin,xmax)            // non uniform mutation
[np,n]=size(xpar);n=n-1;
b=5;
xoff=xpar;
for i=1:n
  u1=rand();u2=rand();
  if (u1<1/2) then
    xoff(1,i)=xpar(1,i)+(xmax(i)-xpar(i))*u2*(1-(ngen-1)/Ngen)^b;
   else
    xoff(1,i)=xpar(1,i)-(xpar(i)-xmin(i))*u2*(1-(ngen-1)/Ngen)^b;
   end;
end;
endfunction
//
//---------------------------------------------------------------------
//------------------------------------------------------------------
//
//
irep=evstr(x_dialog('new computation (0) or resume (1)','0'));
n=evstr(x_dialog('number of parameters (Rastrigin function)','2'));
Xup=evstr(x_dialog('upper bound of the domain','5'));
Xlo=evstr(x_dialog('lower bound of the domain','-5'));
Npop=evstr(x_dialog('population size','10'));
Ngen=evstr(x_dialog('number of generations','30'));
Ngeninit=evstr(x_dialog('number of generations with only exact evaluations','2'));
Nsize=evstr(x_dialog('number of centers used for RBF reconstruction','10'));
Nexacttot=evstr(x_dialog('range of allowable exact evaluation number','70'));
//
deff('val=valalpha(alpha)','Npercent0=[ones(1,Ngeninit),linspace(alpha,1/Npop,Ngen-Ngeninit)];val=sum(int(Npercent0*Npop))+1-Nexacttot');
alpha0=fsolve(0.4,valalpha);alpha0=min(1,alpha0);alpha0=max(0,alpha0);
Npercent0=[ones(1,Ngeninit),linspace(alpha0,1/Npop,Ngen-Ngeninit)];        // percent of exact evaluations at generation gen
//
//------------------------------------------------------------------
//----------------------------------------------------------------------
//
xmin=Xlo*ones(1,n);
xmax=Xup*ones(1,n);
Xmin=ones(Npop,1)*xmin;
Xmax=ones(Npop,1)*xmax;
u=rand(Npop,n);
pop=Xmin+(Xmax-Xmin).*u;  // random initialisation of the population
Xpar=[pop,zeros(Npop,1)];  
Xoff=Xpar;Xoldpar=Xpar;mineval=0;newval=zeros(Npop,1);Traj=[];
Xtot=[];Nexact=[];fmin=[];
//
//
if (irep==1) then,load('session.sav');gen0=gen+1;else gen0=1;end
//
for gen=gen0:Ngen;
Npercent=Npercent0(gen);
//
  Xpar=Xoff;
//
//  ----------------------------------Evaluation phase -------------------------------------
//----------------------------------------------------------------------------------------------
  if (gen<Ngeninit+1) then 
     val=evaluation(Xpar); 
     Nexact=[Nexact,Nexact($)+Npop];
     Xpar(:,n+1)=val;
     Xnew=Xpar;
     Xtot=[Xtot;Xnew];
  else
     val2=evaluationapprox(Xpar,Xtot,Nsize);
     Xpar(:,n+1)=val2;    
     [approxvalues,u]=gsort(-val2);
     u=u(1:int(Npercent*Npop));
     val3=evaluation(Xpar(u,:));
     Xpar(u,$)=val3;
     Xnew=Xpar(u,:);
     Xtot=[Xtot;Xnew];
     Nexact=[Nexact,Nexact($)+length(u)];
  end 
   [Xpar,bestpar]=elitism(Xnew,Xpar,Xtot,Npop,gen);
   fmin=[fmin,bestpar($)];
//
//
//  ----------------------------------------------------------------------------------------
//
  //
  Xpar=rankselect(Xpar);                                // rank selection
  for i=1:2:(Npop-1)
    u1=int(Npop*rand())+1;u2=int(Npop*rand())+1;u3=rand();
    xpar=[Xpar(u1,:);Xpar(u2,:)];Xoff(i:(i+1),:)=xpar;
     xoff=blcross(xpar);                                 //crossover
     Xoff(i:(i+1),:)=xoff; 
       for j=1:2
          xoffl=mutate(xoff(j,:),gen,Ngen,xmin,xmax);   // mutation
          Xoff(i+j-1,:)=xoffl;
       end
 end
save('session.sav');
end
//
//  ----------------------------------last evaluation phase -----------------------------
//----------------------------------------------------------------------------------------
//
if (gen<Ngeninit+1) then 
     val=evaluation(Xpar); 
     Nexact=[Nexact,Nexact($)+Npop];
     Xpar(:,n+1)=val;
     Xnew=Xpar;
     Xtot=[Xtot;Xpar];
else
    val2=evaluationapprox(Xpar,Xtot,Nsize);
    Xpar(:,n+1)=val2;
    [approxvalues,u]=gsort(-val2);
    u=u(1:int(Npercent*Npop));
    val3=evaluation(Xpar(u,:));
    Xpar(u,$)=val3;
    Xnew=Xpar(u,:);
    Xtot=[Xtot;Xnew];
    Nexact=[Nexact,Nexact($)+length(u)];
end
[Xpar,bestpar]=elitism(Xnew,Xpar,Xtot,Npop,gen);
//
fmin=[fmin,bestpar($)];
save('session.sav');
//----------- results displays --------------

xset('window',0);
plot2d(Nexact,fmin,3);
xtitle('convergence history','Neval','fmin');
xset('window',1);
plot2d1('onl',Nexact,fmin,3);
xtitle('convergence history (log scale)','Neval','fmin');
//
disp('minimum obtained:');disp(bestpar(1:n));
disp('corresponding value by f:');disp(bestpar($));
disp('evaluation number:');disp(Nexact($));



/////////////////////////////// VERSION AVEC INTERFACE UNIX
//
////-----------------------------------------------------
//function writedata(str0)
//  unix('rm -f data_heart');
//  fddata=file('open','data_heart','unknown')    
//  write(fddata,str0);
//  write(fddata,'[ensight]');
//  write(fddata,'start     = 0');
//  write(fddata,'verbose   = 4');
//  write(fddata,'#multimesh = 1');
//  file('close',fddata);
//endfunction
//
//
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
//function y=exactfunctionbis(X)       
//z1=X(1);
//theta1=X(2);
//r1=aiL*sqrt(1-(z1/ciL)**2);
//x1=r1*cos(theta1);
//y1=r1*sin(theta1);
//  tx=[t1 t2 t3 f1 f1 f2 x1 x2 x3];
//  ty=[t1 t2 t3 f1 f1 f2 y1 y2 y3];
//  tz=[t1 t2 t3 f1 f1 f2 z1 z2 z3];
//  ts=[ti1 ti2 ti3 tf1 tf2 tf3 s1 s2 s3];
//str0=ts + a + string(tx) + e + string(ty) + e + string(tz) + a;
//writedata(str0);
////
//str2='/afs/inria.fr/rocq/home/reo/ldumas/lifevREO-opt/testsuite/test_bidomain_Schaf/test_bidomain_Schaf -f data_heart >output ';
////
//unix('rm -f time_dep');
//unix('rm -f time_temp');
//unix('rm -f output');
//unix(str2);
////
//y=read('time_temp',1,1);
////
////
//endfunction
////


