new; cls;
tp=1;      /* top-percent rate 0<tp<=1 */
nel=2;     /* number of elitists(suppress when nel=0) <=size */
ng=0;      /* number of replacements in Genitor algorithm(no Genitor algorithm when ng=0) */
gr=0;      /* Gray code(1 or nonzero) or N-ary code(0)       */
sigmas=0;  /* 0(no sharing) or length,length+1,... Usually it can work in Gray code  */
cr=0;      /* 0<=cr<=1, no initialization when cr=0 */
w={1.5,1.3,1,0.7,0.4}; /* weight for each age from age 0 to age rows(w) and the rest is zero. */
d=0;       /* heterogeneous recombination(no recombination if d=0) */
sc=0.5;    /* parameter for power in sharing function. It is linear when sc=1.       */
oc=0;      /* 0:neither polygyny nor parthenogenesis 1:allowing both 2:allowing polygyny only */
size=50;
length=15; base=4; length=31; base=2;
ngen=1000; pc=0.8; npc=1; npc=1; pm=0.1; npm=1; pinv=0; pscr=0; ptlo=0; pdl=0;
ct=-1;      
           /* 1:altogether 2:independently 3:the same pairs differently 4:single at random   */
           /* negative: mix at random between 1 and 4                                        */
ctype=-1;  
           /* 1: one-point crossover 2: n-point crossover 3: uniform crossover 4: arithmetic */
           /* 5: one-point crossover with different lengths 6:n-point with different lengths */
           /* 7: uniform crossover with different points                                     */
           /* negative: mix at random between 1 and 7                                        */
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.5: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                        */
           /* 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.5: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   */
tsize=4; /* tournament size for stype=2 */
tr={0.9,0};  /* 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.2,
0.5 0.3,
0.6 0.35,
0.7 0.4,
0.8 0.5,
0.9 0.7
}; 
fr={
5.5 0.1,
4.4 0.2,
3.3 0.6,
2.2 0.8,
1.1   1};
mu=150;  /* mu=1,2,3,...,size,size+1,size+2,...,rows(fitv)  */
c=2.0;   /* 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)                */
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 */
            /* 5:with reversed candidates 6:augmengted 7:age scaling                     */
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=NGAmD(tp,ng,gr,sigmas,d,cr,w,sc,size,musize,lambda,etype,length,base,ngen,pc,oc,npc,pm,npm,pinv,pscr,ptlo,pdl,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);                   /* a user-defined function */
   local r,nn;
   r=rankindx(fitv,1);  /* in ascending order */
   nn=rows(fitv);
   fitv=r/(2*(c-1)).*(c-sqrt(c^2-4*(c-1)*rndu(nn,1)));
   retp( fitv );
endp;

proc NGAmD(tp,ng,gr,sigmas,d,cr,w,sc,size,musize,lambda,etype,length,base,ngen,pc,oc,npc,pm,npm,pinv,pscr,ptlo,pdl,ct,ctype,stype,tsize,tr,ts,eta,pp,fr,mu,c,range,dim,&f);
   local na,x,y,rp,i,j,l,index,mask,k,fitv,x01,xs,sindex,li,sli,rs,sums,dj,sj,xindex1,xindex2,f:proc;
   local nn,xindex,vk,r,p,cp,minfitv,freq,maxfreq,count,a,b,cs,fitv1,fitv1_1,fitv_1,xp,minx,bestx;
   local xnew,xold,id,size0,nusize,pindex,u,u1,beta,nn1,ri,e,nn_1p,x_1p,fitv_1p,xr,xq,tempbit;
   local tol,ngens,fitvstars,xc,age,bestage,w1,ii,xmul,xtemp1,xtemp2,l1,l2;
   tol=1e-16; ngens=1000; fitvstars=1|zeros(ngens-1,1); /* convergence setting */
   if base^length>2^31;
      errorlog "ERROR: Reduce either length or base.";
      retp(".");
   endif;
   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(4*rndu(rows(stype),1))+1; 
   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(7*rndu(rows(stype),1))+1; /* change 7 to more when we add more options */
   endif;
   if rows(d)==1;
      d=d*ones(rows(stype),1);
   endif;
   if rows(d)<rows(stype);
      d=d|d[rows(d)]*ones(rows(stype)-rows(d),1);
   endif;
   if etype==1;
      size0=size[1];
      nusize=musize;
      size=size-nusize;
   endif;
   if etype/=2 and etype/=3;
      musize=size;
   endif;
   x=floor(base*rndu(size[1],length*dim));
   if gr==0;                           /* n-ary code */
      xp=x;
      x01=xp[.,1:length]*((1/base)^seqa(1,1,length));
      ii=2;
      do while ii<=dim;
         x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/base)^seqa(1,1,length));
         ii=ii+1;
      endo;
   else;
      nn=rows(x);                      /* n-ary Gray code */
      x01=zeros(rows(x),dim);
      k=1;
      do while k<=nn;
         ii=1;
         do while ii<=dim;
            x01[k,ii]=baryGraytodec(x[k,(ii-1)*length+1:ii*length]',base)/(base^length);
            ii=ii+1;
         endo;
         k=k+1;
      endo;
   endif;
   xs=range[1,.]+(range[2,.]-range[1,.]).*x01;
   minfitv=minc(f(xs));
   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;
         if gr==0;                            /* n-ary code */
            xp=x;
            x01=xp[.,1:length]*((1/base)^seqa(1,1,length));
            ii=2;
            do while ii<=dim;
               x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/base)^seqa(1,1,length));
               ii=ii+1;
            endo;
         else;
            nn=rows(x);                       /* n-ary Gray code */
            x01=zeros(rows(x),dim);
            k=1;
            do while k<=nn;
               ii=1;
               do while ii<=dim;
                  x01[k,ii]=baryGraytodec(x[k,(ii-1)*length+1:ii*length]',base)/(base^length);
                  ii=ii+1;
               endo;
               k=k+1;
            endo;
         endif;
         xs=range[1,.]+(range[2,.]-range[1,.]).*x01;
         x=selif(x,rankindx(-f(xs),1).<=musize[j]);
      endif;
      if ctype[j]==1; /* one-point 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,length*dim);
         if ct[j]==1;
            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;
               l=floor((length*dim-1)*rndu(1,1))+1; 
               x[musize[j]+i,.]=x[id[1],1:l]~x[id[2],l+1:length*dim];
               x[musize[j]+i+1,.]=x[id[2],1:l]~x[id[1],l+1:length*dim];
               i=i+2;
            endo;
         elseif ct[j]==2;
            ii=1;
            do while ii<=dim;
               index=rankindx(rndu(musize[j],1),1);
               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=index[i|(i+1)];
                  endif;
                  l=floor((length-1)*rndu(1,1))+1; 
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=x[id[1],(ii-1)*length+1:(ii-1)*length+l]~x[id[2],(ii-1)*length+l+1:ii*length];
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=x[id[2],(ii-1)*length+1:(ii-1)*length+l]~x[id[1],(ii-1)*length+l+1:ii*length];
                  i=i+2;
               endo;
               ii=ii+1;
            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;
               ii=1;
               do while ii<=dim;
                  l=floor((length-1)*rndu(1,1))+1; 
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=x[id[1],(ii-1)*length+1:(ii-1)*length+l]~x[id[2],(ii-1)*length+l+1:ii*length];
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=x[id[2],(ii-1)*length+1:(ii-1)*length+l]~x[id[1],(ii-1)*length+l+1:ii*length];
                  ii=ii+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;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               l=floor((length-1)*rndu(1,1))+1; 
               x[musize[j]+i,.]=x[id[1],.];
               x[musize[j]+i+1,.]=x[id[2],.];
               x[musize[j]+i,(ii-1)*length+1:ii*length]=x[id[1],(ii-1)*length+1:(ii-1)*length+l]~x[id[2],(ii-1)*length+l+1:ii*length];
               x[musize[j]+i+1,(ii-1)*length+1:ii*length]=x[id[2],(ii-1)*length+1:(ii-1)*length+l]~x[id[1],(ii-1)*length+l+1:ii*length];
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==2; /* n-point 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,length*dim);
         if ct[j]==1;
            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;
               l=rankindx((length*dim-1)*rndu(length*dim-1,1),1);
               l=l[1:npc];
               l=0|sortc(l,1)|length*dim;
               k=1;
               do while k<=npc+1;
                  if k%2==1;
                     x[musize[j]+i,l[k]+1:l[k+1]]=x[id[1],l[k]+1:l[k+1]];
                     x[musize[j]+i+1,l[k]+1:l[k+1]]=x[id[2],l[k]+1:l[k+1]];
                  else;
                     x[musize[j]+i,l[k]+1:l[k+1]]=x[id[2],l[k]+1:l[k+1]];
                     x[musize[j]+i+1,l[k]+1:l[k+1]]=x[id[1],l[k]+1:l[k+1]];
                  endif;
                  k=k+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==2;
            ii=1;
            do while ii<=dim;
               index=rankindx(rndu(musize[j],1),1);
               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=index[i|(i+1)];
                  endif;
                  l=rankindx((length-1)*rndu(length-1,1),1);
                  l=l[1:npc];
                  l=0|sortc(l,1)|length;
                  k=1;
                  do while k<=npc+1;
                     if k%2==1;
                        x[musize[j]+i,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[1],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                        x[musize[j]+i+1,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[2],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                     else;
                        x[musize[j]+i,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[2],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                        x[musize[j]+i+1,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[1],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                     endif;
                     k=k+1;
                  endo;
                  i=i+2;
               endo;
               ii=ii+1;
            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;
               ii=1;
               do while ii<=dim;
                  l=rankindx((length-1)*rndu(length-1,1),1);
                  l=l[1:npc];
                  l=0|sortc(l,1)|length;
                  k=1;
                  do while k<=npc+1;
                     if k%2==1;
                        x[musize[j]+i,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[1],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                        x[musize[j]+i+1,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[2],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                     else;
                        x[musize[j]+i,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[2],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                        x[musize[j]+i+1,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[1],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                     endif;
                     k=k+1;
                  endo;
                  ii=ii+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;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               l=rankindx((length-1)*rndu(length-1,1),1);
               l=l[1:npc];
               l=0|sortc(l,1)|length;
               x[musize[j]+i,.]=x[id[1],.];
               x[musize[j]+i+1,.]=x[id[2],.];
               k=1;
               do while k<=npc+1;
                  if k%2==1;
                     x[musize[j]+i,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[1],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                     x[musize[j]+i+1,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[2],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                  else;
                     x[musize[j]+i,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[2],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                     x[musize[j]+i+1,(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]]=x[id[1],(ii-1)*length+l[k]+1:(ii-1)*length+l[k+1]];
                  endif;
                  k=k+1;
               endo;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==3; /* uniform 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,length*dim);
         if ct[j]==1;
            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;
               mask=round(rndu(length*dim,1));
               k=1;
               do while k<=length*dim;
                  if mask[k]==1;
                     x[musize[j]+i,k]=x[id[1],k];
                     x[musize[j]+i+1,k]=x[id[2],k];
                  else;
                     x[musize[j]+i,k]=x[id[2],k];
                     x[musize[j]+i+1,k]=x[id[1],k];
                  endif;
                  k=k+1;
               endo;
               i=i+2;
            endo;
         elseif ct[j]==2;
            ii=1;
            do while ii<=dim;
               index=rankindx(rndu(musize[j],1),1);
               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=index[i|(i+1)];
                  endif;
                  mask=round(rndu(length,1));
                  k=1;
                  do while k<=length;
                     if mask[k]==1;
                        x[musize[j]+i,(ii-1)*length+k]=x[id[1],(ii-1)*length+k];
                        x[musize[j]+i+1,(ii-1)*length+k]=x[id[2],(ii-1)*length+k];
                     else;
                        x[musize[j]+i,(ii-1)*length+k]=x[id[2],(ii-1)*length+k];
                        x[musize[j]+i+1,(ii-1)*length+k]=x[id[1],(ii-1)*length+k];
                     endif;
                     k=k+1;
                  endo;
                  i=i+2;
               endo;
               ii=ii+1;
            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;
               ii=1;
               do while ii<=dim;
                  mask=round(rndu(length,1));
                  k=1;
                  do while k<=length;
                     if mask[k]==1;
                        x[musize[j]+i,(ii-1)*length+k]=x[id[1],(ii-1)*length+k];
                        x[musize[j]+i+1,(ii-1)*length+k]=x[id[2],(ii-1)*length+k];
                     else;
                        x[musize[j]+i,(ii-1)*length+k]=x[id[2],(ii-1)*length+k];
                        x[musize[j]+i+1,(ii-1)*length+k]=x[id[1],(ii-1)*length+k];
                     endif;
                     k=k+1;
                  endo;
                  ii=ii+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;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               mask=round(rndu(length,1));
               x[musize[j]+i,.]=x[id[1],.];
               x[musize[j]+i+1,.]=x[id[2],.];
               k=1;
               do while k<=length;
                  if mask[k]==1;
                     x[musize[j]+i,(ii-1)*length+k]=x[id[1],(ii-1)*length+k];
                     x[musize[j]+i+1,(ii-1)*length+k]=x[id[2],(ii-1)*length+k];
                  else;
                     x[musize[j]+i,(ii-1)*length+k]=x[id[2],(ii-1)*length+k];
                     x[musize[j]+i+1,(ii-1)*length+k]=x[id[1],(ii-1)*length+k];
                  endif;
                  k=k+1;
               endo;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==4; /* 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/2,length*dim);
         if ct[j]==1;
            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;
               xmul=zeros(1,length*dim);
               k=1;
               do while k<=length*dim;
                  xmul[k]=mul(x[id[1],k],x[id[2],k],base);
                  k=k+1;
               endo;
               x[musize[j]+(i-1)/2+1,.]=xmul;
               i=i+2;
            endo;
         elseif ct[j]==2;
            ii=1;
            do while ii<=dim;
               index=rankindx(rndu(musize[j],1),1);
               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=index[i|(i+1)];
                  endif;
                  xmul=zeros(1,length);
                  k=1;
                  do while k<=length;
                     xmul[k]=mul(x[id[1],(ii-1)*length+k],x[id[2],(ii-1)*length+k],base);
                     k=k+1;
                  endo;
                  x[musize[j]+(i-1)/2+1,(ii-1)*length+1:ii*length]=xmul;
                  i=i+2;
               endo;
               ii=ii+1;
            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;
               ii=1;
               do while ii<=dim;
                  xmul=zeros(1,length);
                  k=1;
                  do while k<=length;
                     xmul[k]=mul(x[id[1],(ii-1)*length+k],x[id[2],(ii-1)*length+k],base);
                     k=k+1;
                  endo;
                  x[musize[j]+(i-1)/2+1,(ii-1)*length+1:ii*length]=xmul;
                  ii=ii+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;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               xmul=zeros(1,length);
               x[musize[j]+(i-1)/2+1,.]=x[id[1],.];
               k=1;
               do while k<=length;
                  xmul[k]=mul(x[id[1],(ii-1)*length+k],x[id[2],(ii-1)*length+k],base);
                  k=k+1;
               endo;
               x[musize[j]+(i-1)/2+1,(ii-1)*length+1:ii*length]=xmul;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==5; /* one-point crossover with different lengths */
         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,length*dim);
         if ct[j]==1;
            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;
               l=floor((length*dim-1)*rndu(2,1))+1; 
               xtemp1=x[id[1],1:l[1]]~x[id[2],l[2]+1:length*dim];
               xtemp2=x[id[2],1:l[2]]~x[id[1],l[1]+1:length*dim];
               if cols(xtemp1)>=length*dim;
                  xtemp1=xtemp1[1,1:length*dim];
               else;
                  xtemp1=xtemp1~floor(base*rndu(1,length*dim-cols(xtemp1)));
               endif;
               if cols(xtemp2)>=length*dim;
                  xtemp2=xtemp2[1,1:length*dim];
               else;
                  xtemp2=xtemp2~floor(base*rndu(1,length*dim-cols(xtemp2)));
               endif;
               x[musize[j]+i,.]=xtemp1;
               x[musize[j]+i+1,.]=xtemp2;
               i=i+2;
            endo;
         elseif ct[j]==2;
            ii=1;
            do while ii<=dim;
               index=rankindx(rndu(musize[j],1),1);
               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=index[i|(i+1)];
                  endif;
                  l=floor((length-1)*rndu(2,1))+1; 
                  xtemp1=x[id[1],(ii-1)*length+1:(ii-1)*length+l[1]]~x[id[2],(ii-1)*length+l[2]+1:ii*length];
                  xtemp2=x[id[2],(ii-1)*length+1:(ii-1)*length+l[2]]~x[id[1],(ii-1)*length+l[1]+1:ii*length];
                  if cols(xtemp1)>=length;
                     xtemp1=xtemp1[1,1:length];
                  else;
                     xtemp1=xtemp1~floor(base*rndu(1,length-cols(xtemp1)));
                  endif;
                  if cols(xtemp2)>=length;
                     xtemp2=xtemp2[1,1:length];
                  else;
                     xtemp2=xtemp2~floor(base*rndu(1,length-cols(xtemp2)));
                  endif;
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=xtemp1;
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=xtemp2;
                  i=i+2;
               endo;
               ii=ii+1;
            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;
               ii=1;
               do while ii<=dim;
                  l=floor((length-1)*rndu(2,1))+1; 
                  xtemp1=x[id[1],(ii-1)*length+1:(ii-1)*length+l[1]]~x[id[2],(ii-1)*length+l[2]+1:ii*length];
                  xtemp2=x[id[2],(ii-1)*length+1:(ii-1)*length+l[2]]~x[id[1],(ii-1)*length+l[1]+1:ii*length];
                  if cols(xtemp1)>=length;
                     xtemp1=xtemp1[1,1:length];
                  else;
                     xtemp1=xtemp1~floor(base*rndu(1,length-cols(xtemp1)));
                  endif;
                  if cols(xtemp2)>=length;
                     xtemp2=xtemp2[1,1:length];
                  else;
                     xtemp2=xtemp2~floor(base*rndu(1,length-cols(xtemp2)));
                  endif;
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=xtemp1;
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=xtemp2;
                  ii=ii+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;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               l=floor((length-1)*rndu(2,1))+1; 
               x[musize[j]+i,.]=x[id[1],.];
               x[musize[j]+i+1,.]=x[id[2],.];
               xtemp1=x[id[1],(ii-1)*length+1:(ii-1)*length+l[1]]~x[id[2],(ii-1)*length+l[2]+1:ii*length];
               xtemp2=x[id[2],(ii-1)*length+1:(ii-1)*length+l[2]]~x[id[1],(ii-1)*length+l[1]+1:ii*length];
               if cols(xtemp1)>=length;
                  xtemp1=xtemp1[1,1:length];
               else;
                  xtemp1=xtemp1~floor(base*rndu(1,length-cols(xtemp1)));
               endif;
               if cols(xtemp2)>=length;
                  xtemp2=xtemp2[1,1:length];
               else;
                  xtemp2=xtemp2~floor(base*rndu(1,length-cols(xtemp2)));
               endif;
               x[musize[j]+i,(ii-1)*length+1:ii*length]=xtemp1;
               x[musize[j]+i+1,(ii-1)*length+1:ii*length]=xtemp2;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==6; /* n-point crossover with different lengths */
         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,length*dim);
         if ct[j]==1;
            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;
               l1=rankindx((length*dim-1)*rndu(length*dim-1,1),1);
               l1=l1[1:npc];
               l1=0|sortc(l1,1)|length*dim;
               l2=rankindx((length*dim-1)*rndu(length*dim-1,1),1);
               l2=l2[1:npc];
               l2=0|sortc(l2,1)|length*dim;
               k=1;
               xtemp1=x[id[1],l1[k]+1:l1[k+1]];
               xtemp2=x[id[2],l2[k]+1:l2[k+1]];
               k=2;
               do while k<=npc+1;
                  if k%2==1;
                     xtemp1=xtemp1~x[id[1],l1[k]+1:l1[k+1]];
                     xtemp2=xtemp2~x[id[2],l2[k]+1:l2[k+1]];
                  else;
                     xtemp1=xtemp1~x[id[2],l2[k]+1:l2[k+1]];
                     xtemp2=xtemp2~x[id[1],l1[k]+1:l1[k+1]];
                  endif;
                  k=k+1;
               endo;
               if cols(xtemp1)>=length*dim;
                  xtemp1=xtemp1[1,1:length*dim];
               else;
                  xtemp1=xtemp1~floor(base*rndu(1,length*dim-cols(xtemp1)));
               endif;
               if cols(xtemp2)>=length*dim;
                  xtemp2=xtemp2[1,1:length*dim];
               else;
                  xtemp2=xtemp2~floor(base*rndu(1,length*dim-cols(xtemp2)));
               endif;
               x[musize[j]+i,.]=xtemp1;
               x[musize[j]+i+1,.]=xtemp2;
               i=i+2;
            endo;
         elseif ct[j]==2;
            ii=1;
            do while ii<=dim;
               index=rankindx(rndu(musize[j],1),1);
               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=index[i|(i+1)];
                  endif;
                  l1=rankindx((length-1)*rndu(length-1,1),1);
                  l1=l1[1:npc];
                  l1=0|sortc(l1,1)|length;
                  l2=rankindx((length-1)*rndu(length-1,1),1);
                  l2=l2[1:npc];
                  l2=0|sortc(l2,1)|length;
                  k=1;
                  xtemp1=x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                  xtemp2=x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                  k=2;
                  do while k<=npc+1;
                     if k%2==1;
                        xtemp1=xtemp1~x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                        xtemp2=xtemp2~x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                     else;
                        xtemp1=xtemp1~x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                        xtemp2=xtemp2~x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                     endif;
                     k=k+1;
                  endo;
                  if cols(xtemp1)>=length;
                     xtemp1=xtemp1[1,1:length];
                  else;
                     xtemp1=xtemp1~floor(base*rndu(1,length-cols(xtemp1)));
                  endif;
                  if cols(xtemp2)>=length;
                     xtemp2=xtemp2[1,1:length];
                  else;
                     xtemp2=xtemp2~floor(base*rndu(1,length-cols(xtemp2)));
                  endif;
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=xtemp1;
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=xtemp2;
                  i=i+2;
               endo;
               ii=ii+1;
            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;
               ii=1;
               do while ii<=dim;
                  l1=rankindx((length-1)*rndu(length-1,1),1);
                  l1=l1[1:npc];
                  l1=0|sortc(l1,1)|length;
                  l2=rankindx((length-1)*rndu(length-1,1),1);
                  l2=l2[1:npc];
                  l2=0|sortc(l2,1)|length;
                  k=1;
                  xtemp1=x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                  xtemp2=x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                  k=2;
                  do while k<=npc+1;
                     if k%2==1;
                        xtemp1=xtemp1~x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                        xtemp2=xtemp2~x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                     else;
                        xtemp1=xtemp1~x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                        xtemp2=xtemp2~x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                     endif;
                     k=k+1;
                  endo;
                  if cols(xtemp1)>=length;
                     xtemp1=xtemp1[1,1:length];
                  else;
                     xtemp1=xtemp1~floor(base*rndu(1,length-cols(xtemp1)));
                  endif;
                  if cols(xtemp2)>=length;
                     xtemp2=xtemp2[1,1:length];
                  else;
                     xtemp2=xtemp2~floor(base*rndu(1,length-cols(xtemp2)));
                  endif;
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=xtemp1;
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=xtemp2;
                  ii=ii+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;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               l1=rankindx((length-1)*rndu(length-1,1),1);
               l1=l1[1:npc];
               l1=0|sortc(l1,1)|length;
               l2=rankindx((length-1)*rndu(length-1,1),1);
               l2=l2[1:npc];
               l2=0|sortc(l2,1)|length;
               x[musize[j]+i,.]=x[id[1],.];
               x[musize[j]+i+1,.]=x[id[2],.];
               k=1;
               xtemp1=x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
               xtemp2=x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
               k=2;
               do while k<=npc+1;
                  if k%2==1;
                     xtemp1=xtemp1~x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                     xtemp2=xtemp2~x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                  else;
                     xtemp1=xtemp1~x[id[2],(ii-1)*length+l2[k]+1:(ii-1)*length+l2[k+1]];
                     xtemp2=xtemp2~x[id[1],(ii-1)*length+l1[k]+1:(ii-1)*length+l1[k+1]];
                  endif;
                  k=k+1;
               endo;
               if cols(xtemp1)>=length;
                  xtemp1=xtemp1[1,1:length];
               else;
                  xtemp1=xtemp1~floor(base*rndu(1,length-cols(xtemp1)));
               endif;
               if cols(xtemp2)>=length;
                  xtemp2=xtemp2[1,1:length];
               else;
                  xtemp2=xtemp2~floor(base*rndu(1,length-cols(xtemp2)));
               endif;
               x[musize[j]+i,(ii-1)*length+1:ii*length]=xtemp1;
               x[musize[j]+i+1,(ii-1)*length+1:ii*length]=xtemp2;
               i=i+2;
            endo;
         endif;
      elseif ctype[j]==7; /* uniform crossover with different points */
         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,length*dim);
         if ct[j]==1;
            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;
               mask=round(rndu(length*dim,1));
               l=rankindx(rndu(length*dim,1),1);
               x[musize[j]+i,.]=x[id[1],.];
               x[musize[j]+i+1,.]=x[id[2],.];
               if sumc(mask)/=length*dim;
                  l=l[1:length*dim-sumc(mask)];
                  l=sortc(l,1);
                  count=1;
                  k=1;
                  do while k<=length*dim;
                     if mask[k]==0;
                        x[musize[j]+i,k]=x[id[2],l[count]];
                        x[musize[j]+i+1,l[count]]=x[id[1],k];
                        count=count+1;
                     endif;
                     k=k+1;
                  endo;
               endif;
               i=i+2;
            endo;
         elseif ct[j]==2;
            ii=1;
            do while ii<=dim;
               index=rankindx(rndu(musize[j],1),1);
               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=index[i|(i+1)];
                  endif;
                  mask=round(rndu(length,1));
                  l=rankindx(rndu(length,1),1);
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=x[id[1],(ii-1)*length+1:ii*length];
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=x[id[2],(ii-1)*length+1:ii*length];
                  if sumc(mask)/=length;
                     l=l[1:length-sumc(mask)];
                     l=sortc(l,1);
                     count=1;
                     k=1;
                     do while k<=length;
                        if mask[k]==0;
                           x[musize[j]+i,(ii-1)*length+k]=x[id[2],(ii-1)*length+l[count]];
                           x[musize[j]+i+1,(ii-1)*length+l[count]]=x[id[1],(ii-1)*length+k];
                           count=count+1;
                        endif;
                        k=k+1;
                     endo;
                  endif;
                  i=i+2;
               endo;
               ii=ii+1;
            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;
               ii=1;
               do while ii<=dim;
                  mask=round(rndu(length,1));
                  l=rankindx(rndu(length,1),1);
                  x[musize[j]+i,(ii-1)*length+1:ii*length]=x[id[1],(ii-1)*length+1:ii*length];
                  x[musize[j]+i+1,(ii-1)*length+1:ii*length]=x[id[2],(ii-1)*length+1:ii*length];
                  if sumc(mask)/=length;
                     l=l[1:length-sumc(mask)];
                     l=sortc(l,1);
                     count=1;
                     k=1;
                     do while k<=length;
                        if mask[k]==0;
                           x[musize[j]+i,(ii-1)*length+k]=x[id[2],(ii-1)*length+l[count]];
                           x[musize[j]+i+1,(ii-1)*length+l[count]]=x[id[1],(ii-1)*length+k];
                           count=count+1;
                        endif;
                        k=k+1;
                     endo;
                  endif;
                  ii=ii+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;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               mask=round(rndu(length,1));
               l=rankindx(rndu(length,1),1);
               x[musize[j]+i,.]=x[id[1],.];
               x[musize[j]+i+1,.]=x[id[2],.];
               if sumc(mask)/=length;
                  l=l[1:length-sumc(mask)];
                  l=sortc(l,1);
                  count=1;
                  k=1;
                  do while k<=length;
                     if mask[k]==0;
                        x[musize[j]+i,(ii-1)*length+k]=x[id[2],(ii-1)*length+l[count]];
                        x[musize[j]+i+1,(ii-1)*length+l[count]]=x[id[1],(ii-1)*length+k];
                        count=count+1;
                     endif;
                     k=k+1;
                  endo;
               endif;
               i=i+2;
            endo;
         endif;
      endif;
      if ct[j]==1;
         /* mutation */
         if npm<=0;     /* mutation for parents */
            y=(rndu(musize[j],length*dim).<pm).*(floor((base-1)*rndu(musize[j],length*dim))+1);
            x[1:musize[j],.]=(x[1:musize[j],.]+y)%base;
         else;          /* n-point mutation for parents */
            y=zeros(musize[j],length*dim);
            k=1;
            do while k<=musize[j];
               if rndu(1,1)<pm;
                  rp=rankindx(rndu(length*dim,1),1);
                  rp=rp[1:npm];
                  y[k,rp]=ones(1,npm);
               endif;
               k=k+1;
            endo;
            y=y.*(floor((base-1)*rndu(musize[j],length*dim))+1);
            x[1:musize[j],.]=(x[1:musize[j],.]+y)%base;
         endif;
         /* inversion */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<pinv;
               li=floor(length*dim*rndu(2,1))+1;
               li=sortc(li,1);
               x[k,li[1]:li[2]]=rev(x[k,li[1]:li[2]]')';
            endif;
            k=k+1;
         endo;
         /* scramble */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<pscr;
               li=floor(length*dim*rndu(2,1))+1;
               li=sortc(li,1);
               sli=seqa(li[1],1,li[2]-li[1]+1);
               x[k,li[1]:li[2]]=x[k,sli[rankindx(rndu(li[2]-li[1]+1,1),1)]];
            endif;
            k=k+1;
         endo;
         /* translocation */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<ptlo;
               li=floor(length*dim*rndu(2,1))+1;
               li=sortc(li,1);
               rs=floor((length*dim-(li[2]-li[1]+1)+1)*rndu(1,1))+1;
               x[k,rs:rs+li[2]-li[1]]=x[k,li[1]:li[2]];
            endif;
            k=k+1;
         endo;
         /* deletion */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<pdl;
               li=floor(length*dim*rndu(2,1))+1;
               li=sortc(li,1);
               e=zeros(length*dim,1); 
               e[li[1]:li[2]]=ones(li[2]-li[1]+1,1);
               if li[1]==1 and li[2]==length*dim;
                  x[k,.]=floor(base*rndu(1,length*dim));
               else;
                  x[k,.]=delif(x[k,.]',e)'~floor(base*rndu(1,li[2]-li[1]+1));
               endif;
            endif;
            k=k+1;
         endo;
      elseif ct[j]==2 or ct[j]==3;
         ii=1;
         do while ii<=dim;
           /* mutation */
            if npm<=0;     /* mutation for parents */
               y=(rndu(musize[j],length).<pm).*(floor((base-1)*rndu(musize[j],length))+1);
               x[1:musize[j],(ii-1)*length+1:ii*length]=(x[1:musize[j],(ii-1)*length+1:ii*length]+y)%base;
            else;          /* n-point mutation for parents */
               y=zeros(musize[j],length);
               k=1;
               do while k<=musize[j];
                  if rndu(1,1)<pm;
                     rp=rankindx(rndu(length,1),1);
                     rp=rp[1:npm];
                     y[k,rp]=ones(1,npm);
                  endif;
                  k=k+1;
               endo;
               y=y.*(floor((base-1)*rndu(musize[j],length))+1);
               x[1:musize[j],(ii-1)*length+1:ii*length]=(x[1:musize[j],(ii-1)*length+1:ii*length]+y)%base;
            endif;
            /* inversion */
            k=1;
            do while k<=musize[j];
               if rndu(1,1)<pinv;
                  li=floor(length*rndu(2,1))+1;
                  li=sortc(li,1);
                  x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]]=rev(x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]]')';
               endif;
               k=k+1;
            endo;
            /* scramble */
            k=1;
            do while k<=musize[j];
               if rndu(1,1)<pscr;
                  li=floor(length*rndu(2,1))+1;
                  li=sortc(li,1);
                  sli=seqa(li[1],1,li[2]-li[1]+1);
                  x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]]=x[k,(ii-1)*length+sli[rankindx(rndu(li[2]-li[1]+1,1),1)]];
               endif;
               k=k+1;
            endo;
            /* translocation */
            k=1;
            do while k<=musize[j];
               if rndu(1,1)<ptlo;
                  li=floor(length*rndu(2,1))+1;
                  li=sortc(li,1);
                  rs=floor((length-(li[2]-li[1]+1)+1)*rndu(1,1))+1;
                  x[k,(ii-1)*length+rs:(ii-1)*length+rs+li[2]-li[1]]=x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]];
               endif;
               k=k+1;
            endo;
            /* deletion */
            k=1;
            do while k<=musize[j];
               if rndu(1,1)<pdl;
                  li=floor(length*rndu(2,1))+1;
                  li=sortc(li,1);
                  e=zeros(length,1); 
                  e[li[1]:li[2]]=ones(li[2]-li[1]+1,1);
                  if li[1]==1 and li[2]==length;
                     x[k,(ii-1)*length+1:ii*length]=floor(base*rndu(1,length));
                  else;
                     x[k,(ii-1)*length+1:ii*length]=delif(x[k,(ii-1)*length+1:ii*length]',e)'~floor(base*rndu(1,li[2]-li[1]+1));
                  endif;
               endif;
               k=k+1;
            endo;
            ii=ii+1;
         endo;
      elseif ct[j]==4;
         /* mutation */ 
         if npm<=0;     /* mutation for parents */
            k=1;
            do while k<=musize[j];
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               y=(rndu(1,length).<pm).*(floor((base-1)*rndu(1,length))+1);
               x[k,(ii-1)*length+1:ii*length]=(x[k,(ii-1)*length+1:ii*length]+y)%base;
               k=k+1;
            endo;
         else;          /* n-point mutation for parents */
            k=1;
            do while k<=musize[j];
               if rndu(1,1)<pm;
                  ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
                  y=zeros(1,length);
                  rp=rankindx(rndu(length,1),1);
                  rp=rp[1:npm];
                  y[1,rp]=ones(1,npm);
                  y=y.*(floor((base-1)*rndu(1,length))+1);
                  x[k,(ii-1)*length+1:ii*length]=(x[k,(ii-1)*length+1:ii*length]+y)%base;
               endif;
               k=k+1;
            endo;
         endif;
         /* inversion */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<pinv;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               li=floor(length*rndu(2,1))+1;
               li=sortc(li,1);
               x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]]=rev(x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]]')';
            endif;
            k=k+1;
         endo;
         /* scramble */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<pscr;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               li=floor(length*rndu(2,1))+1;
               li=sortc(li,1);
               sli=seqa(li[1],1,li[2]-li[1]+1);
               x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]]=x[k,(ii-1)*length+sli[rankindx(rndu(li[2]-li[1]+1,1),1)]];
            endif;
            k=k+1;
         endo;
         /* translocation */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<ptlo;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               li=floor(length*rndu(2,1))+1;
               li=sortc(li,1);
               rs=floor((length-(li[2]-li[1]+1)+1)*rndu(1,1))+1;
               x[k,(ii-1)*length+rs:(ii-1)*length+rs+li[2]-li[1]]=x[k,(ii-1)*length+li[1]:(ii-1)*length+li[2]];
            endif;
            k=k+1;
         endo;
         /* deletion */
         k=1;
         do while k<=musize[j];
            if rndu(1,1)<pdl;
               ii=floor(dim*rndu(1,1))+1;   /* randomly chosen from 1 to dim */
               li=floor(length*rndu(2,1))+1;
               li=sortc(li,1);
               e=zeros(length,1); 
               e[li[1]:li[2]]=ones(li[2]-li[1]+1,1);
               if li[1]==1 and li[2]==length;
                  x[k,(ii-1)*length+1:ii*length]=floor(base*rndu(1,length));
               else;
                  x[k,(ii-1)*length+1:ii*length]=delif(x[k,(ii-1)*length+1:ii*length]',e)'~floor(base*rndu(1,li[2]-li[1]+1));
               endif;
            endif;
            k=k+1;
         endo;
      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;
      /* heterogeneous recombination */
      if d[j]/=0;
         nn=rows(x);
         x=x[rankindx(rndu(nn,1),1),.];
         i=1;
         do while i<=nn;
            mask=abs(x[i,.]-x[i+1,.]);
            dj=sumc(mask');
            if dj>d[j];
               k=1;
               do while k<=length*dim;
                  if mask[k]>0 and rndu(1,1)<0.5;
                     tempbit=x[i,k];
                     x[i,k]=x[i+1,k];
                     x[i+1,k]=tempbit;
                  endif;
                  k=k+1;
               endo;    
            endif;  
            i=i+2;
         endo;
      endif;
      /* with reversed setting */
      if etype==5;
         xr=(x+(base-1))%base;
         x=x|xr;
      endif;
      /* augmented setting with random n-point crossover */
      if etype==6;
         nn=rows(x);
         xr=(x+(base-1))%base;
         xq=zeros(rows(x)+1,cols(x));
         i=1;
         do while i<=nn;
            id=floor(nn*rndu(2,1))+1;
            l=rankindx((length*dim-1)*rndu(length*dim-1,1),1);
            npc=floor((length*dim-1)*rndu(1,1))+1; /* random npc between 1 and length-1 */
            l=l[1:npc];
            l=0|sortc(l,1)|length*dim;
            k=1;
            do while k<=npc+1;
               if k%2==1;
                  xq[i,l[k]+1:l[k+1]]=x[id[1],l[k]+1:l[k+1]];
                  xq[i+1,l[k]+1:l[k+1]]=xr[id[2],l[k]+1:l[k+1]];
               else;
                  xq[i,l[k]+1:l[k+1]]=x[id[2],l[k]+1:l[k+1]];
                  xq[i+1,l[k]+1:l[k+1]]=xr[id[1],l[k]+1:l[k+1]];
               endif;
               k=k+1;
            endo;
            i=i+2;
         endo;
         xq=xq[1:nn,.];
         x=x|xr|xq;
      endif;
      /* calculate fitness */
      if gr==0;                           /* n-ary code */
         xp=x;
         x01=xp[.,1:length]*((1/base)^seqa(1,1,length));
         ii=2;
         do while ii<=dim;
            x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/base)^seqa(1,1,length));
            ii=ii+1;
         endo;
      else;
         nn=rows(x);                      /* n-ary Gray code */
         x01=zeros(rows(x),dim);
         k=1;
         do while k<=nn;
            ii=1;
            do while ii<=dim;
               x01[k,ii]=baryGraytodec(x[k,(ii-1)*length+1:ii*length]',base)/(base^length);
               ii=ii+1;
            endo;
            k=k+1;
         endo;
      endif;
      xs=range[1,.]+(range[2,.]-range[1,.]).*x01;
      fitv=zeros(rows(xs),1);
      fitv=f(xs);
      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=sumc(abs(x-x[k,.])');                  /* in terms of original x */
            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);
      e=(sindex.<=nel);
      bestx=selif(x,e);
      if etype==7;
         bestage=selif(age,e);
      endif;
      /* 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 floor(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 or deterministic */
         if stype[j]==0.1;     /* steady-state */
            xnew=x[size[j]+1:rows(x)];
            sindex=rankindx(-fitv[1:size[j]],1);
            xold=selif(x[1:size[j]],sindex.<=(size[j]-(rows(x)-size[j])));
            x=xold|xnew;
         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;              /* generational */
            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[j],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 stype[j]==3.5; /* original Whitley's ranking */
            nn=rows(fitv);
            xindex=floor(nn/(2*(eta-1)).*(eta-sqrt(eta^2-4*(eta-1)*rndu(size[j],1))))+1;
            x=rev(sortc(fitv~x,1)); /* sort them out in descending order */
            x=x[.,2:cols(x)];
            x=x[xindex,.];
         else;
            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,.];
         endif;
         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; /* rulette 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;
         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,.];
      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.";
      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 */
      if gr==0;                           /* n-ary code */
         xp=x;
         x01=xp[.,1:length]*((1/base)^seqa(1,1,length));
         ii=2;
         do while ii<=dim;
            x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/base)^seqa(1,1,length));
            ii=ii+1;
         endo;
      else;
         nn=rows(x);                      /* n-ary Gray code */
         x01=zeros(rows(x),dim);
         k=1;
         do while k<=nn;
            ii=1;
            do while ii<=dim;
               x01[k,ii]=baryGraytodec(x[k,(ii-1)*length+1:ii*length]',base)/(base^length);
               ii=ii+1;
            endo;
            k=k+1;
         endo;
      endif;
      xs=range[1,.]+(range[2,.]-range[1,.]).*x01;
      fitv=f(xs);
      xs=xs[maxindc(fitv),.];
      /* save maxc(fitv)'s to maxfitv */
      fitvstars=maxc(fitv)|fitvstars[1:ngens-1];
      /* cataclysmic mutation */
      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),.];
         x=delif(x,e);
         y=(rndu(rows(x),cols(x)).<cr).*(floor((base-1)*rndu(rows(x),cols(x)))+1);
         x=(x+y)%2;
         x=xc|x;
         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;

proc dectobaryvec(y,NN,nb);
   local x,i;
   x=zeros(nb,1);
   i=1;
   do while i<=nb;
      x[i]=floor(y/(NN^(nb-i)));
      y=y-(NN^(nb-i))*x[i];
      i=i+1;
   endo;
   retp(x);
endp;

proc baryGraytodec(x,n);
   local k,i,y;
   k=rows(x);
   y=x[1]*n^(k-1);
   i=2;
   do while i<=k;
      if sumc(x[1:i-1])%2==1;
         y=y+(n-1-x[i])*n^(k-i);
      else;
         y=y+x[i]*n^(k-i);
      endif;
      i=i+1;
   endo;
   retp(y);
endp;

proc mul(a,b,base);
   local table;
   if base==2;
      let table[2,2]=
      0  0
      0  1
      ;
   elseif base==3;
      let table[3,3]=
      0  0  0
      0  1  2
      0  2  1
      ;
   elseif base==4;
      let table[4,4]=
      0  0  0  0
      0  1  2  3
      0  2  3  1
      0  3  1  2
      ;
   elseif base==5;
      let table[5,5]=
      0  0  0  0  0
      0  1  2  3  4 
      0  2  4  1  3 
      0  3  1  4  2
      0  4  3  2  1
      ;
   elseif base==7;
      let table[7,7]=
      0  0  0  0  0  0  0
      0  1  2  3  4  5  6
      0  2  4  6  1  3  5
      0  3  6  2  5  1  4
      0  4  1  5  2  6  3
      0  5  3  1  6  4  2
      0  6  5  4  3  2  1
      ;
   elseif base==8;
      let table[8,8]=
      0  0  0  0  0  0  0  0
      0  1  2  3  4  5  6  7
      0  2  4  6  3  1  7  5 
      0  3  6  5  7  4  1  2 
      0  4  3  7  6  2  5  1 
      0  5  1  4  2  7  3  6 
      0  6  7  1  5  3  2  4 
      0  7  5  2  1  6  4  3 
      ;
   elseif base==9;
      let table[9,9]=
      0  0  0  0  0  0  0  0  0
      0  1  2  3  4  5  6  7  8
      0  2  1  6  8  7  3  5  4
      0  3  6  2  5  8  1  4  7
      0  4  8  5  6  1  7  2  3
      0  5  7  8  1  3  4  6  2
      0  6  3  1  7  4  2  8  5
      0  7  5  4  2  6  8  3  1
      0  8  4  7  3  2  5  1  6
      ;
   elseif base==11;
      let table[11,11]=
      0  0  0  0  0  0  0  0  0  0  0 
      0  1  2  3  4  5  6  7  8  9 10 
      0  2  4  6  8 10  1  3  5  7  9 
      0  3  6  9  1  4  7 10  2  5  8 
      0  4  8  1  5  9  2  6 10  3  7 
      0  5 10  4  9  3  8  2  7  1  6 
      0  6  1  7  2  8  3  9  4 10  5 
      0  7  3 10  6  2  9  5  1  8  4 
      0  8  5  2 10  7  4  1  9  6  3 
      0  9  7  5  3  1 10  8  6  4  2 
      0 10  9  8  7  6  5  4  3  2  1
      ;
   elseif base==13;
      let table[13,13]=
      0  0  0  0  0  0  0  0  0  0  0  0  0 
      0  1  2  3  4  5  6  7  8  9 10 11 12 
      0  2  4  6  8 10 12  1  3  5  7  9 11 
      0  3  6  9 12  2  5  8 11  1  4  7 10 
      0  4  8 12  3  7 11  2  6 10  1  5  9 
      0  5 10  2  7 12  4  9  1  6 11  3  8 
      0  6 12  5 11  4 10  3  9  2  8  1  7 
      0  7  1  8  2  9  3 10  4 11  5 12  6 
      0  8  3 11  6  1  9  4 12  7  2 10  5 
      0  9  5  1 10  6  2 11  7  3 12  8  4 
      0 10  7  4  1 11  8  5  2 12  9  6  3 
      0 11  9  7  5  3  1 12 10  8  6  4  2 
      0 12 11 10  9  8  7  6  5  4  3  2  1 
      ;
   elseif base==16;
      let table[16,16]=
      0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 
      0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 
      0  2  4  6  8 10 12 14  3  1  7  5 11  9 15 13 
      0  3  6  5 12 15 10  9 11  8 13 14  7  4  1  2 
      0  4  8 12  3  7 11 15  6  2 14 10  5  1 13  9 
      0  5 10 15  7  2 13  8 14 11  4  1  9 12  3  6 
      0  6 12 10 11 13  7  1  5  3  9 15 14  8  2  4 
      0  7 14  9 15  8  1  6 13 10  3  4  2  5 12 11 
      0  8  3 11  6 14  5 13 12  4 15  7 10  2  9  1 
      0  9  1  8  2 11  3 10  4 13  5 12  6 15  7 14 
      0 10  7 13 14  4  9  3 15  5  8  2  1 11  6 12 
      0 11  5 14 10  1 15  4  7 12  2  9 13  6  8  3 
      0 12 11  7  5  9 14  2 10  6  1 13 15  3  4  8 
      0 13  9  4  1 12  8  5  2 15 11  6  3 14 10  7 
      0 14 15  1 13  3  2 12  9  7  6  8  4 10 11  5 
      0 15 13  2  9  6  4 11  1 14 12  3  8  7  5 10
      ;
   else;
      errorlog "ERROR: There is no such a base.";
   endif;
   retp(table[a+1,b+1]);
endp;