new; cls;
ct=-1;      
           /* 2: independently 3: the same pairs differently 4: single at random */
           /* negative: mix at random bewteen 2 and 4                            */
ctype=-1;   
           /* 0: linear crossover 1: blend crossover(alpha) 2: arithmetic crossover     */ 
           /* 3: UNDX crossover(sig) 4: SBX crossover(etasb)                            */
           /* 5: unimodal averaging crossovers(alpha,n) 6: bimodal averaging(alpha,n)   */
           /* negative: mix at random between 1 and 6 except 0                          */
stype=5;   
           /* 0:generational 1:elite 2:tournament 3:ranking 4:roulette 5:expected-value     */
           /* 6:top-percent(random) 7:scaled selection 8:replicative selection              */
           /* 9:inverse tournament                                                          */
           /* 0.1:steady-state 0.2:(mu,lambda) 0.3:(mu+lambda)                              */
           /* 2.1:stochastic for Tsize=2 2.2:stochastic for any Tsize 2.3:EP tournament     */
           /* 2.4:EP tournament with score setting 2.5:stochastic EP tournament             */
           /* 2.6:post unbiased tournament 2.7:pre unbiased tournament 2.8:unbiased elitist */
           /* 2.9:stochastic unbiased tournament(based on 2.6)                              */
           /* 3.1:Baker's linear ranking with mu  3.2:uniform ranking with mu               */
           /* 3.3:exponential ranking with eta  3.4: Whitley's linear ranking               */
           /* 4.1:roulette with linear scaling 5.1:expected-value with linear scaling       */
           /* 4.2:roulette with sigma truncation 5.2:expected-value with sigma truncation   */
           /* 4.3:roulette with power-law scaling 5.3:expected-value with power-law scaling */
           /* 4.4:roulette with Boltzmann scaling 5.4:expected-value with Boltzmann scaling */
           /* 4.5:roulette with window scaling 5.5:expected-value with window scaling       */
           /* 4.6:roulette with custom scaling 5.6:expected-value with custom scaling       */
           /* 4.7:roulette with linear ranking 5.7:expected-value with linear ranking       */
           /* 4.8:roulette with exponential ranking 5.8:expected-value with exp ranking     */
           /* 4.9:roulette with linear normalization 5.9:expected-value with linear norm    */
           /* 7.1:cutting linear 7.2:kinked linear 7.3:power function 7.4:percentile        */
           /* 8.1:with frequencies 8.2:percentile 8.3:flexible percentile(roulette)         */
           /* 9.1:stochastic for T=2 9.2:deterministic for any T 9.3:stochastic for any T   */
tp=1;      /* top-percent rate 0<tp<=1 */
nel=0;     /* number of elitists(suppress when nel=0) <=size */
ng=0;      /* number of replacements in Genitor algorithm(no Genitor algorithm when ng=0) */
sigmas=0;  /* no sharing if sigmas=0 */
cr=0;      /* 0<=cr<=1, no initialization when cr=0 */
w=1-0.2*seqa(0,1,50+1)/50;
sc=0.5;    /* parameter for power in sharing function. It is linear when sc=1.         */
tsize=4; /* tournament size for stype=2 */
tr={0.8,0.1};  /* probability for larger one in stochastic tournament */
ts={3,1};/* s[1]:score for win  s[2]:score for draw */
eta=2;   /* 1<=eta<=2 for stype=3 and 0<eta for stype=3.3 and eta/=1 for stype=3.4  */
pp={
0.3  0.1,
0.4  0.2,
0.5  0.3,
0.6  0.4,
0.7  0.5,
0.8  0.6,
0.85 0.7,
0.9  0.8,
0.95 0.9
};   
fr={
6.3 0.1,
5.2 0.2,
3.9 0.6,
2.2 0.8,
1.5   1};
mu=150;  /* mu=1,2,3,...,size,size+1,size+2,...,rows(fitv) */
c=0.1;   /* constant for stype=4.1&5.1(c=2) and 4.2&5.2(c=1~3) and 4.3&5.3(c[1],[2]>0)      */
         /* and 4.4&5.4(c[1],[2]>0) and 4.5&5.5(c>=0) and 4.6&5.6(c's depend on function)   */ 
         /* and 4.7&5.7(0<=c[1]<c[2]) and 4.8&5.8(0<c<1) and 4.9&5.9(none)                  */
oc=2;    /* 0:neither polygyny nor parthenogenesis 1:allowing both 2:allowing polygyny only */
size=100;
ngen=100; pc=0.8; alpha=0.5; n=5; sig=0.1; etasb=2; d=0.1; s=0; ppert=0.01;
musize=80;  /* musize=size(default) or 0<<musize<=size for (mu,lambda) or (mu+lambda)  */
lambda=150; /* even integer where lambda>size>musize (activated only if musize/=size)  */
etype=0;    /* 0:default 1:steady-state 2:(mu,lambda) 3:(mu+lambda) 4:cross generational 7:age */
range1={0,3}; range2={0,3}; range=range1~range2;
dim=2;
xx={4.2,4.8,6.0,3.7,4.5,3.2,5.3,4.9,4.8,6.8};
y={5.4,5.9,5.9,6.7,3.7,2.3,4.2,3.2,7.3,6.6};
data=xx~y;
xs=RealGA(tp,nel,ng,sigmas,cr,w,sc,size,musize,lambda,etype,ngen,pc,oc,alpha,n,sig,etasb,d,s,ppert,ct,ctype,stype,tsize,tr,ts,eta,pp,fr,mu,c,range,dim,&f);
print xs;





proc f(x);
   local n,sse,x1,y,i;
   n=rows(x);
   sse=zeros(n,1);
   x1=ones(rows(data),1)~data[.,1:cols(data)-1];
   y=data[.,cols(data)];
   i=1;
   do while i<=n;
      sse[i]=sumc((y-x1*x[i,.]')^2);
      i=i+1;
   endo;
   retp(-sse);
endp;

proc userf(fitv,c,j,ngen);                
   retp(fitv);            /* not used at this time */
endp;
   
proc RealGA(tp,nel,ng,sigmas,cr,w,sc,size,musize,lambda,etype,ngen,pc,oc,alpha,n,sig,etasb,d,s,ppert,ct,ctype,stype,tsize,tr,ts,eta,pp,fr,mu,c,range,dim,&f);
   local di,na,x,i,j,index,min,max,x1,x2,fitv,xs,sindex,xtemp,sums,dj,sj,xindex1,xindex2,f:proc;
   local nn,xindex,li,vk,k,r,p,cp,minfitv,freq,maxfreq,count,a,b,cs,fitv1,fitv1_1,fitv_1,minx,bestx;
   local xnew,xold,id,size0,nusize,pindex,u,u1,beta,nn1,ri,e,nn_1p,x_1p,fitv_1p;
   local tol,ngens,fitvstars,xc,dc,age,bestage,w1,wi,betaq,min1,min2,max1,max2;
   tol=1e-16; ngens=50; fitvstars=1|zeros(ngens-1,1); /* convergence setting */
   if rows(size)==1;
      size=size*ones(ngen,1);
   else;
      ngen=rows(size);
      if etype/=0;
         errorlog "Warning: etype is set to be 0.";
         etype=0;
      endif;
   endif;
   musize=musize*ones(ngen,1);
   if rows(stype)==1;
      stype=stype*ones(ngen,1);
   endif;
   if rows(stype)<ngen;
      stype=stype|stype[rows(stype)]*ones(ngen-rows(stype),1);
   endif;
   if rows(ct)==1;
      ct=ct*ones(rows(stype),1);
   elseif rows(ct)<rows(stype);
      ct=ct|ct[rows(ct)]*ones(rows(stype)-rows(ct),1);
   endif;
   if ct<0;
      ct=floor(3*rndu(rows(stype),1))+2; 
   endif;
   if rows(ctype)==1;
      ctype=ctype*ones(rows(stype),1);
   elseif rows(ctype)<rows(stype);
      ctype=ctype|ctype[rows(ctype)]*ones(rows(stype)-rows(ctype),1);
   endif;
   if ctype<0;
      ctype=floor(6*rndu(rows(stype),1))+1; /* change 6 to more when we add more options */
   endif;
   if etype==1;
      size0=size[1];
      nusize=musize;
      size=size-nusize;
   endif;
   if etype/=2 and etype/=3;
      musize=size;
   endif;
   x=range[1,.]+(range[2,.]-range[1,.]).*rndu(size[1],dim);
   minfitv=minc(f(x));
   age=0;
   j=1;
   do while j<=ngen and not abs(fitvstars-fitvstars[1])<=tol;
      if stype[j]==0.1;
         etype=0;
      endif;
      if stype[j]==0.2;
         etype=2;
      endif;
      if stype[j]==0.3;
         etype=3;
         lambda=maxc(lambda|size[j]-musize[j]);
      endif;
      if stype[j]==0;
         pc=1;
      endif;
      na=ceil(round(size[j]*pc)/2)*2;
      if etype==2 or etype==3;
         x=selif(x,rankindx(-f(x),1).<=musize[j]);
      endif;
      if ctype[j]==0;  /* linear crossover */
         index=rankindx(rndu(musize[j],1),1);
         if etype==2;
            na=maxc(lambda|size[j]);
         endif;
         if etype==1 or etype==3;
            na=lambda;
         endif;
         x=x[index,.]|zeros(1.5*na,dim);
         if ct[j]==2;
            i=1;
            do while i<=na;
               di=1;
               do while di<=dim;
                  if oc==1;
                     id=floor(musize[j]*rndu(2,1))+1;
                  elseif oc==2;
                     id=rankindx(rndu(musize[j],1),1);
                  else;
                     id=i|(i+1);
                  endif;
                  if x[id[1],di]<=x[id[2],di];
                     x[musize[j]+3*(i-1)/2+1,di]=maxc((3/2*x[id[1],di]-1/2*x[id[2],di])|range[1,di]);
                     x[musize[j]+3*(i-1)/2+2,di]=1/2*x[id[1],di]+1/2*x[id[2],di];
                     x[musize[j]+3*(i-1)/2+3,di]=minc((-1/2*x[id[1],di]+3/2*x[id[2],di])|range[2,di]);
                  else;
                     x[musize[j]+3*(i-1)/2+1,di]=minc((3/2*x[id[1],di]-1/2*x[id[2],di])|range[2,di]);
                     x[musize[j]+3*(i-1)/2+2,di]=1/2*x[id[1],di]+1/2*x[id[2],di];
                     x[musize[j]+3*(i-1)/2+3,di]=maxc((-1/2*x[id[1],di]+3/2*x[id[2],di])|range[1,di]);
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==3;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=1;
               do while di<=dim;
                  if x[id[1],di]<=x[id[2],di];
                     x[musize[j]+3*(i-1)/2+1,di]=maxc((3/2*x[id[1],di]-1/2*x[id[2],di])|range[1,di]);
                     x[musize[j]+3*(i-1)/2+2,di]=1/2*x[id[1],di]+1/2*x[id[2],di];
                     x[musize[j]+3*(i-1)/2+3,di]=minc((-1/2*x[id[1],di]+3/2*x[id[2],di])|range[2,di]);
                  else;
                     x[musize[j]+3*(i-1)/2+1,di]=minc((3/2*x[id[1],di]-1/2*x[id[2],di])|range[2,di]);
                     x[musize[j]+3*(i-1)/2+2,di]=1/2*x[id[1],di]+1/2*x[id[2],di];
                     x[musize[j]+3*(i-1)/2+3,di]=maxc((-1/2*x[id[1],di]+3/2*x[id[2],di])|range[1,di]);
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==4;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=floor(dim*rndu(1,1))+1;
                  x[musize[j]+3*(i-1)/2+1,.]=x[id[1],.];
                  x[musize[j]+3*(i-1)/2+2,.]=x[id[2],.];
                  x[musize[j]+3*(i-1)/2+3,.]=x[id[2],.];
                  if x[id[1],di]<=x[id[2],di];
                     x[musize[j]+3*(i-1)/2+1,di]=maxc((3/2*x[id[1],di]-1/2*x[id[2],di])|range[1,di]);
                     x[musize[j]+3*(i-1)/2+2,di]=1/2*x[id[1],di]+1/2*x[id[2],di];
                     x[musize[j]+3*(i-1)/2+3,di]=minc((-1/2*x[id[1],di]+3/2*x[id[2],di])|range[2,di]);
                  else;
                     x[musize[j]+3*(i-1)/2+1,di]=minc((3/2*x[id[1],di]-1/2*x[id[2],di])|range[2,di]);
                     x[musize[j]+3*(i-1)/2+2,di]=1/2*x[id[1],di]+1/2*x[id[2],di];
                     x[musize[j]+3*(i-1)/2+3,di]=maxc((-1/2*x[id[1],di]+3/2*x[id[2],di])|range[1,di]);
                  endif;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==1;  /* blend crossover with alpha */
         index=rankindx(rndu(musize[j],1),1);
         if etype==2;
            na=maxc(lambda|size[j]);
         endif;
         if etype==1 or etype==3;
            na=lambda;
         endif;
         x=x[index,.]|zeros(na,dim);
         if ct[j]==2;
            i=1;
            do while i<=na;
               di=1;
               do while di<=dim;
                  if oc==1;
                     id=floor(musize[j]*rndu(2,1))+1;
                  elseif oc==2;
                     id=rankindx(rndu(musize[j],1),1);
                  else;
                     id=i|(i+1);
                  endif;
                  min=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=min+(max-min)*rndu(1,1);
                  x2=min+(max-min)*rndu(1,1);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==3;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=1;
               do while di<=dim;
                  min=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=min+(max-min)*rndu(1,1);
                  x2=min+(max-min)*rndu(1,1);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==4;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=floor(dim*rndu(1,1))+1;
                  x[musize[j]+i,.]=x[id[1],.];
                  x[musize[j]+i+1,.]=x[id[2],.];
                  min=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=min+(max-min)*rndu(1,1);
                  x2=min+(max-min)*rndu(1,1);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==2; /* arithmetic crossover */
         index=rankindx(rndu(musize[j],1),1);
         if etype==2;
            na=maxc(lambda|size[j]);
         endif;
         if etype==1 or etype==3;
            na=lambda;
         endif;
         x=x[index,.]|zeros(na,dim);
         if ct[j]==2;
            i=1;
            do while i<=na;
               di=1;
               do while di<=dim;
                  if oc==1;
                     id=floor(musize[j]*rndu(2,1))+1;
                  elseif oc==2;
                     id=rankindx(rndu(musize[j],1),1);
                  else;
                     id=i|(i+1);
                 endif;
                  wi=rndn(1,1)+0.5;  /* N(0.5,1) */
                  x1=wi*x[id[1],di]+(1-wi)*x[id[2],di];
                  x2=(1-wi)*x[id[1],di]+wi*x[id[2],di];
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==3;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=1;
               do while di<=dim;
                  wi=rndn(1,1)+0.5;  /* N(0.5,1) */
                  x1=wi*x[id[1],di]+(1-wi)*x[id[2],di];
                  x2=(1-wi)*x[id[1],di]+wi*x[id[2],di];
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==4;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=floor(dim*rndu(1,1))+1;
                  x[musize[j]+i,.]=x[id[1],.];
                  x[musize[j]+i+1,.]=x[id[2],.];
                  wi=rndn(1,1)+0.5;  /* N(0.5,1) */
                  x1=wi*x[id[1],di]+(1-wi)*x[id[2],di];
                  x2=(1-wi)*x[id[1],di]+wi*x[id[2],di];
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==3; /* UNDX crossover */
         index=rankindx(rndu(musize[j],1),1);
         if etype==2;
            na=maxc(lambda|size[j]);
         endif;
         if etype==1 or etype==3;
            na=lambda;
         endif;
         x=x[index,.]|zeros(na,dim);
         if ct[j]==2;
            i=1;
            do while i<=na;
               di=1;
               do while di<=dim;
                  if oc==1;
                     id=floor(musize[j]*rndu(2,1))+1;
                  elseif oc==2;
                     id=rankindx(rndu(musize[j],1),1);
                  else;
                     id=i|(i+1);
                  endif;
                  x1=(x[id[1],di]+x[id[2],di])/2+sig*(x[id[2],di]-x[id[1],di])*rndn(1,1);
                  x2=(x[id[1],di]+x[id[2],di])/2+sig*(x[id[2],di]-x[id[1],di])*rndn(1,1);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                    x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==3;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=1;
               do while di<=dim;
                  x1=(x[id[1],di]+x[id[2],di])/2+sig*(x[id[2],di]-x[id[1],di])*rndn(1,1);
                  x2=(x[id[1],di]+x[id[2],di])/2+sig*(x[id[2],di]-x[id[1],di])*rndn(1,1);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                    x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==4;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=floor(dim*rndu(1,1))+1;
                  x[musize[j]+i,.]=x[id[1],.];
                  x[musize[j]+i+1,.]=x[id[2],.];
                  x1=(x[id[1],di]+x[id[2],di])/2+sig*(x[id[2],di]-x[id[1],di])*rndn(1,1);
                  x2=(x[id[1],di]+x[id[2],di])/2+sig*(x[id[2],di]-x[id[1],di])*rndn(1,1);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                    x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==4; /* SBX crossover */
         index=rankindx(rndu(musize[j],1),1);
         if etype==2;
            na=maxc(lambda|size[j]);
         endif;
         if etype==1 or etype==3;
            na=lambda;
         endif;
         x=x[index,.]|zeros(na,dim);
         if ct[j]==2;
            i=1;
            do while i<=na;
               di=1;
               do while di<=dim;
                  if oc==1;
                     id=floor(musize[j]*rndu(2,1))+1;
                  elseif oc==2;
                     id=rankindx(rndu(musize[j],1),1);
                  else;
                     id=i|(i+1);
                  endif;
                  u=rndu(1,1);
                  if u<=0.5;
                     betaq=(2*u)^(1/(etasb+1));
                  else;
                     betaq=(1/(2*(1-u)))^(1/(etasb+1));
                  endif;
                  x1=0.5*((1+betaq)*x[id[1],di]+(1-betaq)*x[id[2],di]);
                  x2=0.5*((1-betaq)*x[id[1],di]+(1+betaq)*x[id[2],di]);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==3;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=1;
               do while di<=dim;
                  u=rndu(1,1);
                  if u<=0.5;
                     betaq=(2*u)^(1/(etasb+1));
                  else;
                     betaq=(1/(2*(1-u)))^(1/(etasb+1));
                  endif;
                  x1=0.5*((1+betaq)*x[id[1],di]+(1-betaq)*x[id[2],di]);
                  x2=0.5*((1-betaq)*x[id[1],di]+(1+betaq)*x[id[2],di]);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==4;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=floor(dim*rndu(1,1))+1;
                  x[musize[j]+i,.]=x[id[1],.];
                  x[musize[j]+i+1,.]=x[id[2],.];
                  u=rndu(1,1);
                  if u<=0.5;
                     betaq=(2*u)^(1/(etasb+1));
                  else;
                     betaq=(1/(2*(1-u)))^(1/(etasb+1));
                  endif;
                  x1=0.5*((1+betaq)*x[id[1],di]+(1-betaq)*x[id[2],di]);
                  x2=0.5*((1-betaq)*x[id[1],di]+(1+betaq)*x[id[2],di]);
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==5; /* unimodal averaing n blend crossovers */
         index=rankindx(rndu(musize[j],1),1);
         if etype==2;
            na=maxc(lambda|size[j]);
         endif;
         if etype==1 or etype==3;
            na=lambda;
         endif;
         x=x[index,.]|zeros(na,dim);
         if ct[j]==2;
            i=1;
            do while i<=na;
               di=1;
               do while di<=dim;
                  if oc==1;
                     id=floor(musize[j]*rndu(2,1))+1;
                  elseif oc==2;
                     id=rankindx(rndu(musize[j],1),1);
                  else;
                     id=i|(i+1);
                  endif;
                  min=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=sumc(min+(max-min)*rndu(n,1))/n;
                  x2=sumc(min+(max-min)*rndu(n,1))/n;
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==3;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=1;
               do while di<=dim;
                  min=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=sumc(min+(max-min)*rndu(n,1))/n;
                  x2=sumc(min+(max-min)*rndu(n,1))/n;
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==4;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=floor(dim*rndu(1,1))+1;
                  x[musize[j]+i,.]=x[id[1],.];
                  x[musize[j]+i+1,.]=x[id[2],.];
                  min=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=sumc(min+(max-min)*rndu(n,1))/n;
                  x2=sumc(min+(max-min)*rndu(n,1))/n;
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==6; /* bimodal averaing n blend crossovers */
         index=rankindx(rndu(musize[j],1),1);
         if etype==2;
            na=maxc(lambda|size[j]);
         endif;
         if etype==1 or etype==3;
            na=lambda;
         endif;
         x=x[index,.]|zeros(na,dim);
         if ct[j]==2;
            i=1;
            do while i<=na;
               di=1;
               do while di<=dim;
                  if oc==1;
                     id=floor(musize[j]*rndu(2,1))+1;
                  elseif oc==2;
                     id=rankindx(rndu(musize[j],1),1);
                  else;
                     id=i|(i+1);
                  endif;
                  min1=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max1=minc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  min2=maxc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max2=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=sumc(min1+(max1-min1)*rndu(n,1))/n;
                  x2=sumc(min2+(max2-min2)*rndu(n,1))/n;
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==3;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=1;
               do while di<=dim;
                  min1=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max1=minc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  min2=maxc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max2=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=sumc(min1+(max1-min1)*rndu(n,1))/n;
                  x2=sumc(min2+(max2-min2)*rndu(n,1))/n;
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
                  di=di+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==4;
            i=1;
            do while i<=na;
               if oc==1;
                  id=floor(musize[j]*rndu(2,1))+1;
               elseif oc==2;
                  id=rankindx(rndu(musize[j],1),1);
               else;
                  id=i|(i+1);
               endif;
               di=floor(dim*rndu(1,1))+1;
                  x[musize[j]+i,.]=x[id[1],.];
                  x[musize[j]+i+1,.]=x[id[2],.];
                  min1=minc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max1=minc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  min2=maxc(x[id[1],di]|x[id[2],di])-alpha*abs(x[id[1],di]-x[id[2],di]);
                  max2=maxc(x[id[1],di]|x[id[2],di])+alpha*abs(x[id[1],di]-x[id[2],di]);
                  x1=sumc(min1+(max1-min1)*rndu(n,1))/n;
                  x2=sumc(min2+(max2-min2)*rndu(n,1))/n;
                  if x1<range[1,di];
                     x[musize[j]+i,di]=range[1,di];
                  elseif range[2,di]<x1;
                     x[musize[j]+i,di]=range[2,di];
                  else;
                     x[musize[j]+i,di]=x1;
                  endif;
                  if x2<range[1,di];
                     x[musize[j]+i+1,di]=range[1,di];
                  elseif range[2,di]<x2;
                     x[musize[j]+i+1,di]=range[2,di];
                  else;
                     x[musize[j]+i+1,di]=x2;
                  endif;
               i=i+2;
            endo;
         endif;
      endif;
      if ct[j]==2 or ct[j]==3;
         /* perturbation(uniform) */
         if d/=0; /* no uniform perturbation if d=0 */
            i=1;
            do while i<=musize[j];
               if rndu(1,1)<ppert;
                  di=1;
                  do while di<=dim;
                     xtemp=x[i,di]-d+2*d*rndu(1,1);
                     if xtemp<range[1,di];
                        x[i,di]=range[1,di];
                     elseif range[2,di]<xtemp;
                        x[i,di]=range[2,di];
                     else;
                        x[i,di]=xtemp;
                     endif;
                     di=di+1;
                  endo;
               endif;
               i=i+1;
            endo;
         endif;
         /* perturbation(normal) */
         if s/=0; /* no uniform perturbation if s=0 */
            i=1;
            do while i<=musize[j];
               if rndu(1,1)<ppert;
                  di=1;
                  do while di<=dim;
                     xtemp=x[i,di]+s*rndn(1,1);
                     if xtemp<range[1,di];
                        x[i,di]=range[1,di];
                     elseif range[2,di]<xtemp;
                        x[i,di]=range[2,di];
                     else;
                        x[i,di]=xtemp;
                     endif;
                     di=di+1;
                  endo;
               endif;
               i=i+1;
            endo;
         endif;
      elseif ct[j]==4;
         /* perturbation(uniform) */
         if d/=0; /* no uniform perturbation if d=0 */
            i=1;
            do while i<=musize[j];
               if rndu(1,1)<ppert;
                  di=floor(dim*rndu(1,1))+1;
                     xtemp=x[i,di]-d+2*d*rndu(1,1);
                     if xtemp<range[1,di];
                        x[i,di]=range[1,di];
                     elseif range[2,di]<xtemp;
                        x[i,di]=range[2,di];
                     else;
                        x[i,di]=xtemp;
                     endif;
               endif;
               i=i+1;
            endo;
         endif;
         /* perturbation(normal) */
         if s/=0; /* no uniform perturbation if s=0 */
            i=1;
            do while i<=musize[j];
               if rndu(1,1)<ppert;
                  di=floor(dim*rndu(1,1))+1;
                     xtemp=x[i,di]+s*rndn(1,1);
                     if xtemp<range[1,di];
                        x[i,di]=range[1,di];
                     elseif range[2,di]<xtemp;
                        x[i,di]=range[2,di];
                     else;
                        x[i,di]=xtemp;
                     endif;
               endif;
               i=i+1;
            endo;
         endif;
      endif;
      /* age */
      if etype==7;
         if age==0;
            age=ones(musize[j],1)|zeros(rows(x)-musize[j],1);
         else;
            age=(age+1)|zeros(rows(x)-musize[j],1);
         endif;
      endif;
      /* calculate fitness(the larger, the better) */
      fitv=f(x);
      if etype==7;
         if rows(w)<rows(stype);
            w1=w|zeros(rows(stype)-rows(w)+1,1);
         else;
            w1=w;
         endif;
         fitv=fitv.*w1[(age+1)];
      endif;
      /* sharing */
      if sigmas/=0;      
         nn=rows(fitv);
         sums=zeros(nn,1);
         k=1;
         do while k<=nn;
            dj=sqrt(sumc(((x-x[k,.])^2)'));
            sj=((1-dj/sigmas)^sc).*(dj.<=sigmas);
            sums[k]=sumc(sj);
            k=k+1;
         endo;
         fitv=fitv./sums;
      endif;
      /* keep the best nel candidates as bestx */
      if nel>size[j];
         errorlog "ERROR: Parameter nel<=size.";
         nel=size[j];
      endif;
      sindex=rankindx(-fitv,1);
      bestx=selif(x,sindex.<=nel);
      /* genitor algorithm */
      r=rankindx(-fitv,1);
      nn=rows(fitv);
      p=1/nn*(eta-(2*eta-2)*(r-1)/(nn-1));
      cp=cumsumc(p);
      k=1;
      do while k<=ng; 
         xindex1=sumc(rndu(1,1).>cp)+1;
         xindex2=sumc(rndu(1,1).>cp)+1;
         minx=minindc(fitv);
         if rndu(1,1)<0.5;
            x[minx,.]=x[xindex1,.];
            fitv[minx]=fitv[xindex1];
            if etype==7;
               age[minx]=age[xindex1];
            endif;
         else;
            x[minx,.]=x[xindex2,.];
            fitv[minx]=fitv[xindex2];
            if etype==7;
               age[minx]=age[xindex2];
            endif;
         endif;
         k=k+1;
      endo;
      /* top-percent pre-selection */
      if stype[j]/=0 and stype[j]/=1;
         sindex=rankindx(-fitv,1);
         e=(sindex.<=(rows(fitv)*tp));
         x=selif(x,e);
         fitv=selif(fitv,e);
         if etype==7;
            age=selif(age,e);
         endif;
      endif;
      /* for (mu,lambda) */
      if etype==2;
         x=x[musize[j]+1:rows(x),.];
         fitv=fitv[musize[j]+1:rows(fitv)];
      endif;
      /* steady-state based on the first musize of parents */
      if etype==1;
         sindex=rankindx(-fitv[1:size0],1);
         xold=selif(x[1:size0,.],sindex.<=nusize[j]);
         x=x[size0+1:rows(x),.];
         fitv=fitv[size0+1:rows(fitv)];
      endif;
      /* cross generational setting */
      if etype==4;
         nn_1p=rows(x);
         if j/=1;
            x=x|x_1p;
            fitv=fitv|fitv_1p;
         endif;
         x_1p=x[1:nn_1p,.];
         fitv_1p=fitv[1:nn_1p];
      endif;
      /* selection */
      if floor(stype[j])==0;         /* generational */
         if stype[j]==0.1; /* steady-state */
            xnew=x[size[j]+1:rows(x),.];
            sindex=rankindx(-fitv[1:size[j]],1);
            if rows(x)-size[j]==size[j];
               x=xnew;
            else;
               xold=selif(x[1:size[j],.],sindex.<=(size[j]-(rows(x)-size[j])));
               x=xold|xnew;
            endif;
         elseif stype[j]==0.2; /* (mu,lambda)  */
            x=x[1:size[j],.];
         elseif stype[j]==0.3; /* (mu+lambda)  */
            x=x[1:size[j],.];
         else;
            x=x[size[j]+1:2*size[j],.];
         endif;
      elseif stype[j]==1;             /* elite */
         sindex=rankindx(-fitv,1);
         e=(sindex.<=size[j]);
         x=selif(x,e);
         if etype==7;
            age=selif(age,e);
         endif;
      elseif floor(stype[j])==2;             /* tournament */
         if stype[j]==2.3;      /* EP tournament */
            ts={1,0.5};
            xindex=eptournament(fitv,size[j],tsize,ts);
         elseif stype[j]==2.4;  /* EP tournament with score setting */
            xindex=eptournament(fitv,size[j],tsize,ts);
         elseif stype[j]==2.5;  /* stochastic EP tournament */
            xindex=septournament(fitv,size[j],tsize,tr,ts);
         elseif stype[j]==2.6;  /* post unbiased tournament */
            nn=rows(fitv);
            pindex=calcpindex(nn);
            xindex=zeros(size[j],1);
            k=1;
            do while k<=size[j];
               i=floor(nn*rndu(1,1))+1;
               if fitv[i]>fitv[pindex[i]];
                  xindex[k]=i;
               else;
                  xindex[k]=pindex[i];
               endif;
               k=k+1;
            endo;
         elseif stype[j]==2.7;    /* pre unbiased tournament */
            /* pre-random selection */
            nn=rows(fitv);
            sindex=rankindx(rndu(nn,1),1);
            sindex=sindex[1:size[j]];
            x=x[sindex,.];
            fitv=fitv[sindex];
            /* unbiased tournament */
            pindex=calcpindex(size[j]);
            xindex=zeros(size,1);
            k=1;
            do while k<=size[j];
               if fitv[k]>fitv[pindex[k]];
                  xindex[k]=k;
               else;
                  xindex[k]=pindex[k];
               endif;
               k=k+1;
            endo;
         elseif stype[j]==2.8;    /* unbiased tournament and then elitist */
            /* unbiased tournament */
            nn=rows(fitv);
            pindex=calcpindex(nn);
            xindex=zeros(nn,1);
            k=1;
            do while k<=nn;
               if fitv[k]>fitv[pindex[k]];
                  xindex[k]=k;
               else;
                  xindex[k]=pindex[k];
               endif;
               k=k+1;
            endo;
            x=x[xindex,.];
            fitv=fitv[xindex];
            /* elitist */
            sindex=rankindx(-fitv,1);
            xindex=selif(seqa(1,1,nn),sindex.<=size[j]);
         elseif stype[j]==2.9;    /* stochastic unbiased tournament */
            nn=rows(fitv);
            pindex=calcpindex(nn);
            xindex=zeros(size[j],1);
            k=1;
            do while k<=size[j];
               i=floor(nn*rndu(1,1))+1;
               if fitv[i]>fitv[pindex[i]];
                  if rndu(1,1)<tr[1];
                     xindex[k]=i;
                  else;
                     xindex[k]=pindex[i];
                  endif;
               else;
                  if rndu(1,1)<tr[1];
                     xindex[k]=pindex[i];
                  else;
                     xindex[k]=i;
                  endif;
               endif;
               k=k+1;
            endo;
         else;
            nn=rows(fitv);
            xindex=zeros(size[j],1);
            k=1;
            do while k<=size[j]; 
               li=rankindx(rndu(nn,1),1);
               if stype[j]==2.1; /* stocastic tournament for T=2 */
                  li=li[1:2];
                  vk=fitv[li];
                  if rndu(1,1)<tr[1];
                     xindex[k]=li[maxindc(vk)];
                  else;
                     xindex[k]=li[minindc(vk)];
                  endif;
               elseif stype[j]==2.2; /* stochastic tournament for any T */
                  li=li[1:tsize];
                  vk=fitv[li];
                  xindex[k]=li[stournament(vk,tr[1])];
               else;
                  li=li[1:tsize];
                  vk=fitv[li];
                  xindex[k]=li[tournament(vk)];
               endif;
               k=k+1;
            endo;
         endif;
         x=x[xindex,.];
         if etype==7;
            age=age[xindex];
         endif;
      elseif floor(stype[j])==3;      /* ranking */
         if mu>rows(fitv);
            mu=rows(fitv);
         endif;
         if stype[j]==3.1; /* Baker's linear with mu */
            r=rankindx(-fitv,1);
            nn=rows(fitv);
            p=zeros(nn,1);
            k=1;
            do while k<=nn;
               if r[k]<=mu;
                  p[k]=1/mu*(eta-(2*eta-2)*(r[k]-1)/(mu-1));
               else;
                  p[k]=0; /* p=0 for ranking mu+1,mu+2,...,nn */
               endif;
               k=k+1;
            endo;
         elseif stype[j]==3.2; /* uniform with mu */
            r=rankindx(-fitv,1);
            nn=rows(fitv);
            p=zeros(nn,1);
            k=1;
            do while k<=nn;
               if r[k]<=mu;
                  p[k]=1/mu;
               else;
                  p[k]=0; /* p=0 for ranking mu+1,mu+2,...,nn */
               endif;
               k=k+1;
            endo;
         elseif stype[j]==3.3; /* exponential */
            r=rankindx(fitv,1);  /* in ascending order */
            nn=rows(fitv);
            p=(1-exp(-r))/sumc(1-exp(-seqa(nn,-1,nn)));
         elseif stype[j]==3.4; /* Whitley's linear */
            r=rankindx(fitv,1);  /* in ascending order */
            nn=rows(fitv);
            p=r/(2*(eta-1)).*(eta-sqrt(eta^2-4*(eta-1)*rndu(nn,1)));
            p=p/sumc(p);
         else;
            r=rankindx(-fitv,1);
            nn=rows(fitv);
            p=1/nn*(eta-(2*eta-2)*(r-1)/(nn-1));
         endif;
         cp=cumsumc(p);
         xindex=zeros(size[j],1);
         k=1;
         do while k<=size[j]; 
            xindex[k]=sumc(rndu(1,1).>cp)+1;
            k=k+1;
         endo;
         x=x[xindex,.];
         if etype==7;
            age=age[xindex];
         endif;
      elseif floor(stype[j])==4;      /* roulette */
         if minfitv<0;
            fitv=fitv-minfitv;
         endif;
         k=1;
         do while k<=rows(fitv);
            if fitv[k]<0;
               fitv[k]=0;
            endif;
            k=k+1;
         endo;
         if stype[j]==4.1;    /* roulette with linear scaling */
            if minc(fitv)>(c[1]*meanc(fitv)-maxc(fitv))/(c[1]-1);
               if maxc(fitv)-meanc(fitv)/=0;
                  a=meanc(fitv)*(c[1]-1)/(maxc(fitv)-meanc(fitv));
                  b=meanc(fitv)*(maxc(fitv)-c[1]*meanc(fitv))/(maxc(fitv)-meanc(fitv));
               endif;
            else;
               if meanc(fitv)-minc(fitv)/=0;
                  a=meanc(fitv)/(meanc(fitv)-minc(fitv));
                  b=-minc(fitv)*meanc(fitv)/(meanc(fitv)-minc(fitv));
               endif;
            endif;
            fitv=a*fitv+b;
         elseif stype[j]==4.2; /* roulette with sigma truncation */
            fitv=fitv-(meanc(fitv)-c[1]*stdc(fitv));
         elseif stype[j]==4.3; /* roulette with power-law scaling */
            cs=(c[1]+(j-1)/(ngen-1)*(c[2]-c[1]));
            fitv=fitv^cs;
         elseif stype[j]==4.4; /* roulette with Boltzmann scaling */
            cs=(c[1]+(j-1)/(ngen-1)*(c[2]-c[1]));
            fitv=exp(fitv/cs)./exp(meanc(fitv)/cs);   /* or just fitv=exp(fitv/cs) */
         elseif stype[j]==4.5; /* roulette with window scaling */
            fitv=fitv-minc(fitv)+c[1];
         elseif stype[j]==4.6; /* roulette with custom scaling */       
            fitv=userf(fitv,c,j,ngen);
         elseif stype[j]==4.7; /* roulette with linear ranking */
            r=rankindx(-fitv,1);
            nn=rows(fitv);
            fitv=c[2]-(r-1)*(c[2]-c[1])/(nn-1);
         elseif stype[j]==4.8; /* roulette with exponential ranking */
            r=rankindx(-fitv,1);
            fitv=c*(1-c)^(r-1);   /* or fitv=c^(r-1); */
         elseif stype[j]==4.9; /* roulette with linear normalization */
            r=rankindx(fitv,1); 
            nn=rows(fitv);
            fitv1=zeros(nn,1);
            fitv1[1]=1;
            k=1;
            do while k<=nn;
               if r[k]==1;
                  fitv1[k]=1;
                  fitv_1=fitv[k];
                  fitv1_1=1;
               endif;
               k=k+1;
            endo;
            i=2;
            do while i<=nn;
               k=1;
               do while k<=nn;
                  if r[k]==i;
                     if fitv_1/=fitv[k];
                        fitv1[k]=fitv1_1+1;
                        fitv1_1=fitv1_1+1;
                     else;
                        fitv1[k]=fitv1_1;
                     endif;
                     fitv_1=fitv[k];
                  endif;
                  k=k+1;
               endo;
               i=i+1;
            endo;
            fitv=fitv1;
         endif;
         k=1;
         do while k<=rows(fitv);
            if fitv[k]<0;
               fitv[k]=0;
            endif;
            k=k+1;
         endo;
         p=fitv/sumc(fitv);
         cp=cumsumc(p);
         xindex=zeros(size[j],1);
         k=1;
         do while k<=size[j]; 
            xindex[k]=sumc(rndu(1,1).>cp)+1;
            k=k+1;
         endo;
         x=x[xindex,.];
         if etype==7;
            age=age[xindex];
         endif;
      elseif floor(stype[j])==5;        /* expected-value */
         if minfitv<0;
            fitv=fitv-minfitv;
         endif;
         k=1;
         do while k<=rows(fitv);
            if fitv[k]<0;
               fitv[k]=0;
            endif;
            k=k+1;
         endo;
         if stype[j]==5.1;    /* expected-value with linear scaling */
            if minc(fitv)>(c[1]*meanc(fitv)-maxc(fitv))/(c[1]-1);
               if maxc(fitv)-meanc(fitv)/=0;
                  a=meanc(fitv)*(c[1]-1)/(maxc(fitv)-meanc(fitv));
                  b=meanc(fitv)*(maxc(fitv)-c[1]*meanc(fitv))/(maxc(fitv)-meanc(fitv));
               endif;
            else;
               if meanc(fitv)-minc(fitv)/=0;
                  a=meanc(fitv)/(meanc(fitv)-minc(fitv));
                  b=-minc(fitv)*meanc(fitv)/(meanc(fitv)-minc(fitv));
               endif;
            endif;
            fitv=a*fitv+b;
         elseif stype[j]==5.2; /* expected-value with sigma truncation */
            fitv=fitv-(meanc(fitv)-c[1]*stdc(fitv));
         elseif stype[j]==5.3; /* expected-value with power-law scaling */
            cs=(c[1]+(j-1)/(ngen-1)*(c[2]-c[1]));
            fitv=fitv^cs;
         elseif stype[j]==5.4; /* expected-value with Boltzmann scaling */
            cs=(c[1]+(j-1)/(ngen-1)*(c[2]-c[1]));
            fitv=exp(fitv/cs)./exp(meanc(fitv)/cs);   /* or just fitv=exp(fitv/cs) */
         elseif stype[j]==5.5; /* expected-value with window scaling */
            fitv=fitv-minc(fitv)+c[1];
         elseif stype[j]==5.6; /* expected-value with custom scaling */       
            fitv=userf(fitv,c,j,ngen);
         elseif stype[j]==5.7; /* expected-value with linear ranking */
            r=rankindx(-fitv,1);
            nn=rows(fitv);
            fitv=c[2]-(r-1)*(c[2]-c[1])/(nn-1);
         elseif stype[j]==5.8; /* expected-value with exponential ranking */
            r=rankindx(-fitv,1);
            fitv=c*(1-c)^(r-1);   /* or fitv=c^(r-1); */
         elseif stype[j]==5.9; /* expected-value with linear normalization */
            r=rankindx(fitv,1); 
            nn=rows(fitv);
            fitv1=zeros(nn,1);
            fitv1[1]=1;
            k=1;
            do while k<=nn;
               if r[k]==1;
                  fitv1[k]=1;
                  fitv_1=fitv[k];
                  fitv1_1=1;
               endif;
               k=k+1;
            endo;
            i=2;
            do while i<=nn;
               k=1;
               do while k<=nn;
                  if r[k]==i;
                     if fitv_1/=fitv[k];
                        fitv1[k]=fitv1_1+1;
                        fitv1_1=fitv1_1+1;
                     else;
                        fitv1[k]=fitv1_1;
                     endif;
                     fitv_1=fitv[k];
                  endif;
                  k=k+1;
               endo;
               i=i+1;
            endo;
            fitv=fitv1;
         endif;
         k=1;
         do while k<=rows(fitv);
            if fitv[k]<0;
               fitv[k]=0;
            endif;
            k=k+1;
         endo;
         p=fitv/sumc(fitv);
         nn=rows(fitv);
         freq=(p*size[j]-floor(p*size[j]))~p*size[j]~seqa(1,1,nn);
         freq=sortc(freq,1);
         maxfreq=maxc(floor(freq[.,2]));
         xindex=zeros(size[j],1);
         count=1;
         i=maxfreq;
         do while i>=0;
            k=nn;
            do while k>0; 
               if floor(freq[k,2])>=i and count<=size[j];
                  xindex[count]=freq[k,3];
                  count=count+1;
               endif;
               k=k-1;
            endo;
            i=i-1;
         endo;
         if xindex==0;
            xindex=rev(freq[rows(freq)-rows(xindex)+1:rows(freq),3]);
         endif;
         x=x[xindex,.];
         if etype==7;
            age=age[xindex];
         endif;
      elseif stype[j]==6;        /* top-percent(complete random when tp=1) */
         xindex=floor(rows(x)*rndu(size[j],1))+1;
         x=x[xindex,.];
         if etype==7;
            age=age[xindex];
         endif;
      elseif floor(stype[j])==7; /* scaled selection(original Whitley's ranking) */
         nn=rows(fitv);
         u=rndu(size[j],1);
         if stype[j]==7.1; /* cutting linear */
            if eta<0 or eta>1;
               errorlog "ERROR: Parameter eta is out of range. 0<eta<1.";
               retp(".");
            endif;
            u1=eta*u;
         elseif stype[j]==7.2; /* kinked linear */
            u1=(pp[2]/pp[1])*u.*(u.<pp[1])+(((1-pp[2])/(1-pp[1]))*(u-pp[1])+pp[2]).*(u.>=pp[1]);
         elseif stype[j]==7.3; /* power function */
            u1=u^eta;
         elseif stype[j]==7.4; /* percentile */
            if pp[1,1]/=0;
               pp=zeros(1,2)|pp;
            endif;
            if pp[rows(pp),1]/=1;
               pp=pp|ones(1,2);
            endif;
            beta=(pp[2:rows(pp),2]-pp[1:rows(pp)-1,2])./(pp[2:rows(pp),1]-pp[1:rows(pp)-1,1]);
            r=zeros(size[j],1);
            nn1=rows(pp);
            k=1;
            do while k<=size[j];
               i=1;
               do while i<=nn1;
                  r[k]=r[k]+(u[k].>=pp[i,1]);
                  i=i+1;
               endo;
               k=k+1;
            endo;
            u1=beta[r].*(u-pp[r,1])+pp[r,2];
         else;
            u1=1/(2*(eta-1))*(eta-sqrt(eta^2-4*(eta-1)*u));
         endif;
         xindex=floor(nn*u1)+1;
         x=rev(sortc(fitv~x,1)); /* sort them out in descending order */
         x=x[.,2:cols(x)];
         x=x[xindex,.];
         if etype==7;
            age=age[xindex];
         endif;
      elseif floor(stype[j])==8; /* replicative selection */
         if stype[j]==8.1; /* frequency setting */
            fr=fr[1:4,1]~((1/4)|(1/2)|(3/4)|1);
         elseif stype[j]==8.2 or stype[j]==8.3; /* percentile */
            fr=fr;
         else;
            fr=(4|3|2|1)~((1/4)|(1/2)|(3/4)|1);
         endif;
         nn=rows(fitv);
         ri=rankindx(-fitv,1);
         r=ones(nn,1); /* from ones */
         nn1=rows(fr);
         k=1;
         do while k<=nn;
            r[k]=sumc(ri[k].>(fr[.,2]*nn))+1;
            k=k+1;
         endo;
         r=fr[r,1];
         if stype[j]==8.3;
            p=r/sumc(r);
            cp=cumsumc(p);
            xindex=zeros(size[j],1);
            k=1;
            do while k<=size[j]; 
               xindex[k]=sumc(rndu(1,1).>cp)+1;
               k=k+1;
            endo;
         else;
            nn1=sumc(r);
            xindex=zeros(nn1,1);
            count=1;
            k=1;
            do while k<=nn;
               i=1;
               do while i<=r[k];
                  xindex[count]=k;
                  count=count+1;
                  i=i+1;
               endo;
               k=k+1;
            endo;
            xindex=xindex[floor(nn1*rndu(size[j],1))+1];
         endif;
         x=x[xindex,.];
         if etype==7;
            age=age[xindex];
         endif;
      elseif floor(stype[j])==9; /* inverse tournament */
         if rows(fitv)<=size[j];
            errorlog "ERROR: Increase the crossover rate.";
            retp(".");
         endif;
         do while rows(fitv)>size[j]; 
            nn=rows(fitv);
            li=rankindx(rndu(nn,1),1);
            if stype[j]==9.2 or stype[j]==9.3;
               if nn-size[j]>=tsize-1;
                  li=li[1:tsize];
               else;
                  li=li[1:nn-size[j]+1];
               endif;
               vk=fitv[li];
               if stype[j]==9.2;
                  li=delif(li,seqa(1,1,rows(li)).==tournament(vk));
               else;
                  li=delif(li,seqa(1,1,rows(li)).==stournament(vk,tr[1]));
               endif;
            else;
               li=li[1:2];
               vk=fitv[li];
               if stype[j]==9.1;
                  if rndu(1,1)<tr[1];
                     li=li[minindc(vk)];
                  else;
                     li=li[maxindc(vk)];
                  endif;
               else;
                  li=li[minindc(vk)];
               endif;
            endif;
            e=zeros(nn,1);
            e[li]=ones(rows(li),1);
            x=delif(x,e);
            fitv=delif(fitv,e);
            if etype==7;
               age=delif(age,e);
            endif;
         endo;
      else;                        /* others to be updated */
         errorlog "ERROR: The number of stype must be 0 or 1 or 2 or 3 or 4 or 5 or 6 or 7 or 8 or 9.";
         retp(".");
      endif;
      /* combine with elitists */
      if nel/=0;
         if nel==size[j];
            x=bestx;
            if etype==7;
               age=bestage;
            endif;
         else;
            x=bestx|x[1:size[j]-nel,.]; 
            if etype==7;
               age=bestage|age[1:size[j]-nel];
            endif;
         endif;
      endif;
      /* only if steady-state */
      if etype==1;
         x=xold|x;
      endif;
      /* find x to maximize fitv */
      fitv=f(x);
      xs=x[maxindc(fitv),.];
      /* save maxc(fitv)'s to maxfitv */
      fitvstars=maxc(fitv)|fitvstars[1:ngens-1];
      /* initialization */
      if cr/=0 and abs(fitvstars-fitvstars[1])<=tol;
         fitvstars=maxc(fitv)|zeros(ngens-1,1);
         e=zeros(rows(x),1);
         e[maxindc(fitv)]=1;
         xc=x[maxindc(fitv),.];
         dc=minc(abs(range-xc'))';
         x=xc|((xc-cr*dc)+2*cr*dc.*rndu(rows(e)-1,dim));
         if etype==7;
            age=age[maxindc(fitv)]|zeros(rows(e)-1,1);
         endif;
      endif;
      j=j+1;
   endo;
   retp(xs);
endp;

proc tournament(v);
   local n,index,v1,index1,i,j;
   index=seqa(1,1,rows(v));
   j=1;
   do while rows(index)>=2;
      i=1;
         if v[i]>v[i+1];
            v1=v[i];
            index1=index[i];
         else;
            v1=v[i+1];
            index1=index[i+1];
         endif;
      n=rows(index);
      i=3;
      do while i<=floor(n/2)*2;
         if v[i]>v[i+1];
            v1=v1|v[i];
            index1=index1|index[i];
         else;
            v1=v1|v[i+1];
            index1=index1|index[i+1];
         endif;
         i=i+2;
      endo;
      if i==n;
         v1=v1|v[i];
         index1=index1|index[i];
      endif;
      v=v1;
      index=index1;
      j=j+1;
   endo;
   retp(index);
endp;

proc stournament(v,tr);
   local n,index,v1,index1,i,j;
   index=seqa(1,1,rows(v));
   j=1;
   do while rows(index)>=2;
      i=1;
         if v[i]>v[i+1];
            if rndu(1,1)<tr;
               v1=v[i];
               index1=index[i];
            else;
               v1=v[i+1];
               index1=index[i+1];
            endif;
         else;
            if rndu(1,1)<tr;
               v1=v[i+1];
               index1=index[i+1];
            else;
               v1=v[i];
               index1=index[i];
            endif;
         endif;
      n=rows(index);
      i=3;
      do while i<=floor(n/2)*2;
         if v[i]>v[i+1];
            if rndu(1,1)<tr;
               v1=v1|v[i];
               index1=index1|index[i];
            else;
               v1=v1|v[i+1];
               index1=index1|index[i+1];
            endif;
         else;
            if rndu(1,1)<tr;
               v1=v1|v[i+1];
               index1=index1|index[i+1];
            else;
               v1=v1|v[i];
               index1=index1|index[i];
            endif;
         endif;
         i=i+2;
      endo;
      if i==n;
         v1=v1|v[i];
         index1=index1|index[i];
      endif;
      v=v1;
      index=index1;
      j=j+1;
   endo;
   retp(index);
endp;

proc eptournament(v,size,tsize,ts);
   local nn,index,vscore,i,j,iop,sindex;
   nn=rows(v);
   index=seqa(1,1,nn);
   vscore=zeros(nn,1);
   j=1;
   do while j<=nn;
      iop=rankindx(rndu(nn-1,1),1);
      iop=iop[1:tsize];
      iop=iop+ones(tsize,1).*(iop.>=j);
      i=1;
      do while i<=tsize;
         if v[j]>v[iop[i]];
           vscore[j]=vscore[j]+ts[1];
         elseif v[j]<v[iop[i]]; 
           vscore[iop[i]]=vscore[iop[i]]+ts[1];     
         else;
           vscore[j]=vscore[j]+ts[2];
           vscore[iop[i]]=vscore[iop[i]]+ts[2];   
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   sindex=rankindx(-vscore,1);
   index=selif(index,sindex.<=size);
   retp(index);
endp;

proc septournament(v,size,tsize,tr,ts);
   local nn,index,vscore,i,j,iop,sindex,r;
   nn=rows(v);
   index=seqa(1,1,nn);
   vscore=zeros(nn,1);
   j=1;
   do while j<=nn;
      iop=rankindx(rndu(nn-1,1),1);
      iop=iop[1:tsize];
      iop=iop+ones(tsize,1).*(iop.>=j);
      i=1;
      do while i<=tsize;
         if v[j]>v[iop[i]];
           r=rndu(1,1);
           if r<tr[1];
              vscore[j]=vscore[j]+ts[1];
           elseif r<tr[1]+tr[2];
              vscore[j]=vscore[j]+ts[2];
              vscore[iop[i]]=vscore[iop[i]]+ts[2];   
           else;
              vscore[j]=vscore[iop[i]]+ts[1]; 
           endif;
         elseif v[j]<v[iop[i]]; 
           r=rndu(1,1);
           if r<tr[1];
              vscore[j]=vscore[iop[i]]+ts[1];
           elseif r<tr[1]+tr[2];
              vscore[j]=vscore[j]+ts[2];
              vscore[iop[i]]=vscore[iop[i]]+ts[2];   
           else;
              vscore[j]=vscore[j]+ts[1]; 
           endif;   
         else;
           vscore[j]=vscore[j]+ts[2];
           vscore[iop[i]]=vscore[iop[i]]+ts[2];   
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   sindex=rankindx(-vscore,1);
   index=selif(index,sindex.<=size);
   retp(index);
endp;

proc calcpindex(n);
   local index,rindex,tempindex,r,i,j;
   i=1;
   do while i<1e+256;
      index=seqa(1,1,n);
      rindex=zeros(n,1);
      j=1;
      do while j<=n-1;
         r=floor((n-j)*rndu(1,1))+1;
         tempindex=delif(index,index.==j);
         rindex[j]=tempindex[r];
         index=delif(index,index.==rindex[j]);
         j=j+1;
      endo;
      j=n;
         rindex[j]=index;
      if rindex[n]/=n;
         break;
      endif;
      i=i+1;
   endo;
   retp(rindex);
endp;