new; cls;
size=100;
fn df(x)=(1/2)^(x-1); /* distance function */
ngen=100; alpha=0.8; n=3; sig=0.1; etasb=0.5; d=0; s=0.01; ppert=0.05;
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                           */
lstype=7;  /* 0: random 1: elitist 2: probabilistic 3: random by distance 4: random&prob */
           /* 5: 2D dish random 6: 2D dish random&probabilistic                          */
           /* 7: 3D bulb random 8: 3D bulb random&probabilistic                          */
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=CellularRealGA(size,ngen,alpha,n,sig,etasb,d,s,ppert,ctype,lstype,range,dim,&df,&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 CellularRealGA(size,ngen,alpha,n,sig,etasb,d,s,ppert,ctype,lstype,range,dim,&df,&f);
   local size1,x,fitv,minfitv,fitv0,i,j,k,di,rorder,i1,i2,i3,ii,inc,inew,seq,p,cp,ks,f:proc;
   local min,max,xs,x1,x2,x3,fitv1,fitv2,fitv3,wi,u,betaq,min1,min2,max1,max2,e,actinc,mdist;
   local rii,v,w,cdata,df:proc;
   if rows(ctype)==1;
      ctype=ctype*ones(ngen,1);
   elseif rows(ctype)<ngen;
      ctype=ctype|ctype[rows(ctype)]*ones(ngen-rows(ctype),1);
   endif;
   if ctype<0;
      ctype=floor(6*rndu(ngen,1))+1; /* change 6 to more when we add more options */
   endif;
   if rows(d)==1;
      d=d*ones(dim,1);
   endif;
   if rows(s)==1;
      s=s*ones(dim,1);
   endif;
   size1=prodc(size);
   if lstype==5 or lstype==6;
      cdata=zeros(size1,2);
      i=1;
      do while i<=size1;
         u=2*rndu(1,1)-1;
         v=2*rndu(1,1)-1;
         if sqrt(u^2+v^2)<=1;
            cdata[i,.]=u~v;
            i=i+1;
         endif;
      endo;
   endif;
   if lstype==7 or lstype==8;
      cdata=zeros(size1,3);
      i=1;
      do while i<=size1;
         u=2*rndu(1,1)-1;
         v=2*rndu(1,1)-1;
         w=2*rndu(1,1)-1;
         if sqrt(u^2+v^2+w^2)<=1;
            cdata[i,.]=u~v~w;
            i=i+1;
         endif;
      endo;
   endif;
   x=range[1,.]+(range[2,.]-range[1,.]).*rndu(size1,dim);
   fitv=f(x);
   minfitv=minc(fitv);
   fitv=fitv-minc(minfitv|0);
   j=1;
   do while j<=ngen;
      rorder=rankindx(rndu(size1,1),1);
      i=1;
      do while i<=size1;
         if lstype==3 or lstype==4;
            e=zeros(size1,1); e[rorder[i]]=1;
            seq=seqa(1,1,size1);
            seq=delif(seq,e);
            if rows(size)==1;
               ii=seqa(1,1,size1);
               rii=rorder[i];
            elseif rows(size)==2;
               i1=ceil(seq/size[2]); 
               i2=seq-size[2]*(ceil(seq/size[2])-1);
               ii=i1~i2;
               i1=ceil(rorder[i]/size[2]); 
               i2=rorder[i]-size[2]*(ceil(rorder[i]/size[2])-1);
               rii=i1~i2;
            elseif rows(size)==3;
               i3=ceil(seq/(size[1]*size[2]));
               i1=ceil((seq-(i3-1)*size[1]*size[2])/size[2]);
               i2=(seq-(i3-1)*size[1]*size[2])-size[2]*(ceil((seq-(i3-1)*size[1]*size[2])/size[2])-1);
               ii=i1~i2~i3;
               i3=ceil(rorder[i]/(size[1]*size[2]));
               i1=ceil((rorder[i]-(i3-1)*size[1]*size[2])/size[2]);
               i2=(rorder[i]-(i3-1)*size[1]*size[2])-size[2]*(ceil((rorder[i]-(i3-1)*size[1]*size[2])/size[2])-1);
               rii=i1~i2~i3;
            endif;
            actinc=zeros(size1-1,rows(size));
            di=1;
            do while di<=rows(size);
               k=1;
               do while k<=size1-1;
                  actinc[k,di]=minc(abs(ii[k,di]-rii[di])|(minc(ii[k,di]|rii[di])+size[di]-maxc(ii[k,di]|rii[di])));
                  k=k+1;
               endo;
               di=di+1;
            endo;
            mdist=maxc(actinc');
         elseif lstype==5 or lstype==6 or lstype==7 or lstype==8;
            e=zeros(size1,1); e[rorder[i]]=1;
            seq=seqa(1,1,size1);
            seq=delif(seq,e);
            mdist=zeros(size1,1);
            k=1;
            do while k<=size1;
               mdist[k]=sqrt((cdata[k,.]-cdata[rorder[i],.])*(cdata[k,.]-cdata[rorder[i],.])');
               k=k+1;
            endo;
            mdist=delif(mdist,e);
         else;
            if rows(size)==1;
               i1=rorder[i];
               let inc[2,1]=
               -1
                1
               ;
               inew=i1+inc;
            elseif rows(size)==2;
               i1=ceil(rorder[i]/size[2]); 
               i2=rorder[i]-size[2]*(ceil(rorder[i]/size[2])-1);
               let inc[8,2]=
               -1 -1
                0 -1
                1 -1
               -1  0
                1  0
               -1  1
                0  1
                1  1
               ;
               inew=(i1~i2)+inc;
            elseif rows(size)==3;
               i3=ceil(rorder[i]/(size[1]*size[2]));
               i1=ceil((rorder[i]-(i3-1)*size[1]*size[2])/size[2]);
               i2=(rorder[i]-(i3-1)*size[1]*size[2])-size[2]*(ceil((rorder[i]-(i3-1)*size[1]*size[2])/size[2])-1);
               let inc[26,3]=
               -1 -1 -1
                0 -1 -1
                1 -1 -1
               -1  0 -1
                0  0 -1
                1  0 -1
               -1  1 -1
                0  1 -1
                1  1 -1
               -1 -1  0
                0 -1  0
                1 -1  0
               -1  0  0
                1  0  0
               -1  1  0
                0  1  0
                1  1  0
               -1 -1  1
                0 -1  1
                1 -1  1
               -1  0  1
                0  0  1
                1  0  1
               -1  1  1
                0  1  1
                1  1  1
               ;
               inew=(i1~i2~i3)+inc;
            endif;
            di=1;
            do while di<=rows(size);
               k=1;
               do while k<=rows(inc);
                  if inew[k,di]==0;
                     inew[k,di]=size[di];
                  elseif inew[k,di]==size[di]+1;
                     inew[k,di]=1;
                  endif;
                  k=k+1;
               endo;
               di=di+1;
            endo;
            if rows(size)==1;
               seq=inew;
            elseif rows(size)==2;
               seq=(size[2]*(inew[.,1]-1)+inew[.,2]);
            elseif rows(size)==3;
               seq=size[2]*(inew[.,1]-1)+inew[.,2]+(inew[.,3]-1)*size[1]*size[2];
            endif;
         endif;
         /* pick one of neighborhoods */
         if lstype==0;
            ks=seq[floor(rows(seq)*rndu(1,1))+1];
         elseif lstype==1;
            ks=seq[maxindc(fitv[seq])];
         elseif lstype==2;
            fitv0=fitv[seq];
            p=fitv0/sumc(fitv0);
            cp=cumsumc(p);
            ks=seq[sumc(rndu(1,1).>cp)+1];
         elseif lstype==3 or lstype==5 or lstype==7;
            fitv0=df(mdist);
            p=fitv0/sumc(fitv0);
            cp=cumsumc(p);
            ks=seq[sumc(rndu(1,1).>cp)+1];
         elseif lstype==4 or lstype==6 or lstype==8;
            fitv0=df(mdist).*fitv[seq];
            p=fitv0/sumc(fitv0);
            cp=cumsumc(p);
            ks=seq[sumc(rndu(1,1).>cp)+1];
         endif;
         if ctype[j]==0; /* linear crossover between rorder[i]-th and ks-th */
            x1=zeros(1,dim);
            x2=zeros(1,dim);
            x3=zeros(1,dim);
            di=1;
            do while di<=dim;
               if x[rorder[i],di]<=x[ks,di];
                  x1[di]=maxc((3/2*x[rorder[i],di]-1/2*x[ks,di])|range[1,di]);
                  x2[di]=1/2*x[rorder[i],di]+1/2*x[ks,di];
                  x3[di]=minc((-1/2*x[rorder[i],di]+3/2*x[ks,di])|range[2,di]);
               else;
                  x1[di]=minc((3/2*x[rorder[i],di]-1/2*x[ks,di])|range[2,di]);
                  x2[di]=1/2*x[rorder[i],di]+1/2*x[ks,di];
                  x3[di]=maxc((-1/2*x[rorder[i],di]+3/2*x[ks,di])|range[1,di]);
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            fitv2=f(x2);
            fitv3=f(x3);
            if f(x[rorder[i],.])<maxc(fitv1|fitv2|fitv3);
               if fitv1>fitv2 and fitv1>fitv3;
                  fitv1=fitv1-minc(minfitv|0);
                  x[rorder[i],.]=x1;
                  fitv[rorder[i]]=fitv1;
               elseif fitv2>fitv1 and fitv2>fitv3;
                  fitv2=fitv2-minc(minfitv|0);
                  x[rorder[i],.]=x2;
                  fitv[rorder[i]]=fitv2;
               else;
                  fitv3=fitv3-minc(minfitv|0);
                  x[rorder[i],.]=x3;
                  fitv[rorder[i]]=fitv3;
               endif;
            endif;
         elseif ctype[j]==1; /* blend crossover between rorder[i]-th and ks-th */
            x1=zeros(1,dim);
            di=1;
            do while di<=dim;
               min=minc(x[rorder[i],di]|x[ks,di])-alpha*abs(x[rorder[i],di]-x[ks,di]);
               max=maxc(x[rorder[i],di]|x[ks,di])+alpha*abs(x[rorder[i],di]-x[ks,di]);
               x1[di]=min+(max-min)*rndu(1,1);
               if x1[di]<range[1,di];
                  x1[di]=range[1,di];
               elseif range[2,di]<x1[di];
                  x1[di]=range[2,di];
               endif;
               di=di+1;
            endo;      
            fitv1=f(x1);
            if f(x[rorder[i],.])<fitv1;
               fitv1=fitv1-minc(minfitv|0);
               x[rorder[i],.]=x1;
               fitv[rorder[i]]=fitv1;
            endif;
         elseif ctype[j]==2; /* arithmetic crossover between rorder[i]-th and ks-th */
            x1=zeros(1,dim);
            x2=zeros(1,dim);
            di=1;
            do while di<=dim;
               wi=rndn(1,1)+0.5;  /* N(0.5,1) */
               x1[di]=wi*x[rorder[i],di]+(1-wi)*x[ks,di];
               x2[di]=(1-wi)*x[rorder[i],di]+wi*x[ks,di];
               if x1[di]<range[1,di];
                  x1[di]=range[1,di];
               elseif range[2,di]<x1[di];
                  x1[di]=range[2,di];
               endif;
               if x2[di]<range[1,di];
                  x2[di]=range[1,di];
               elseif range[2,di]<x2[di];
                  x2[di]=range[2,di];
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            fitv2=f(x2);
            if f(x[rorder[i],.])<maxc(fitv1|fitv2);
               if fitv1>fitv2;
                  fitv1=fitv1-minc(minfitv|0);
                  x[rorder[i],.]=x1;
                  fitv[rorder[i]]=fitv1;
               else;
                  fitv2=fitv2-minc(minfitv|0);
                  x[rorder[i],.]=x2;
                  fitv[rorder[i]]=fitv2;
               endif;
            endif;
         elseif ctype[j]==3; /* UNDX crossover between rorder[i]-th and ks-th */
            x1=zeros(1,dim);
            di=1;
            do while di<=dim;
               x1[di]=(x[rorder[i],di]+x[ks,di])/2+sig*(x[ks,di]-x[rorder[i],di])*rndn(1,1);
               if x1[di]<range[1,di];
                  x1[di]=range[1,di];
               elseif range[2,di]<x1[di];
                  x1[di]=range[2,di];
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            if f(x[rorder[i],.])<fitv1;
               fitv1=fitv1-minc(minfitv|0);
               x[rorder[i],.]=x1;
               fitv[rorder[i]]=fitv1;
            endif;
         elseif ctype[j]==4; /* SBX crossover between rorder[i]-th and ks-th */
            x1=zeros(1,dim);
            x2=zeros(1,dim);
            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[di]=0.5*((1+betaq)*x[rorder[i],di]+(1-betaq)*x[ks,di]);
               x2[di]=0.5*((1-betaq)*x[rorder[i],di]+(1+betaq)*x[ks,di]);
               if x1[di]<range[1,di];
                  x1[di]=range[1,di];
               elseif range[2,di]<x1[di];
                  x1[di]=range[2,di];
               endif;
               if x2[di]<range[1,di];
                  x2[di]=range[1,di];
               elseif range[2,di]<x2[di];
                  x2[di]=range[2,di];
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            fitv2=f(x2);
            if f(x[rorder[i],.])<maxc(fitv1|fitv2);
               if fitv1>fitv2;
                  fitv1=fitv1-minc(minfitv|0);
                  x[rorder[i],.]=x1;
                  fitv[rorder[i]]=fitv1;
               else;
                  fitv2=fitv2-minc(minfitv|0);
                  x[rorder[i],.]=x2;
                  fitv[rorder[i]]=fitv2;
               endif;
            endif;
         elseif ctype[j]==5; /* unimodal averaing n blend crossover between rorder[i]-th and ks-th */
            x1=zeros(1,dim);
            di=1;
            do while di<=dim;
               min=minc(x[rorder[i],di]|x[ks,di])-alpha*abs(x[rorder[i],di]-x[ks,di]);
               max=maxc(x[rorder[i],di]|x[ks,di])+alpha*abs(x[rorder[i],di]-x[ks,di]);
               x1[di]=sumc(min+(max-min)*rndu(n,1))/n;
               if x1[di]<range[1,di];
                  x1[di]=range[1,di];
               elseif range[2,di]<x1[di];
                  x1[di]=range[2,di];
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            if f(x[rorder[i],.])<fitv1;
               fitv1=fitv1-minc(minfitv|0);
               x[rorder[i],.]=x1;
               fitv[rorder[i]]=fitv1;
            endif;
         elseif ctype[j]==6; /* bimodal averaing n blend crossover between rorder[i]-th and ks-th */
            x1=zeros(1,dim);
            x2=zeros(1,dim);
            di=1;
            do while di<=dim;
               min1=minc(x[rorder[i],di]|x[ks,di])-alpha*abs(x[rorder[i],di]-x[ks,di]);
               max1=minc(x[rorder[i],di]|x[ks,di])+alpha*abs(x[rorder[i],di]-x[ks,di]);
               min2=maxc(x[rorder[i],di]|x[ks,di])-alpha*abs(x[rorder[i],di]-x[ks,di]);
               max2=maxc(x[rorder[i],di]|x[ks,di])+alpha*abs(x[rorder[i],di]-x[ks,di]);
               x1[di]=sumc(min1+(max1-min1)*rndu(n,1))/n;
               x2[di]=sumc(min2+(max2-min2)*rndu(n,1))/n;
               if x1[di]<range[1,di];
                  x1[di]=range[1,di];
               elseif range[2,di]<x1[di];
                  x1[di]=range[2,di];
               endif;
               if x2[di]<range[1,di];
                  x2[di]=range[1,di];
               elseif range[2,di]<x2[di];
                  x2[di]=range[2,di];
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            fitv2=f(x2);
            if f(x[rorder[i],.])<maxc(fitv1|fitv2);
               if fitv1>fitv2;
                  fitv1=fitv1-minc(minfitv|0);
                  x[rorder[i],.]=x1;
                  fitv[rorder[i]]=fitv1;
               else;
                  fitv2=fitv2-minc(minfitv|0);
                  x[rorder[i],.]=x2;
                  fitv[rorder[i]]=fitv2;
               endif;
            endif;
         endif;
         /* perturbation(uniform) */
         if d/=0; /* no uniform perturbation if d=0 */
            x1=zeros(1,dim);
            di=1;
            do while di<=dim;
               if rndu(1,1)<ppert;
                  x1[di]=x[rorder[i],di]-d[di]+2*d[di]*rndu(1,1);
                  if x1[di]<range[1,di];
                     x1[di]=range[1,di];
                  elseif range[2,di]<x1[di];
                     x1[di]=range[2,di];
                  endif;
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            if f(x[rorder[i],.])<fitv1;             
               fitv1=fitv1-minc(minfitv|0);
               x[rorder[i],.]=x1;
               fitv[rorder[i]]=fitv1;
            endif;
         endif;
         if s/=0; /* no normal perturbation if s=0 */
            x1=zeros(1,dim);
            di=1;
            do while di<=dim;
               if rndu(1,1)<ppert;
                  x1[di]=x[rorder[i],di]+s[di]*rndn(1,1);
                  if x1[di]<range[1,di];
                     x1[di]=range[1,di];
                  elseif range[2,di]<x1[di];
                     x1[di]=range[2,di];
                  endif;
               endif;
               di=di+1;
            endo;
            fitv1=f(x1);
            if f(x[rorder[i],.])<fitv1;             
               fitv1=fitv1-minc(minfitv|0);
               x[rorder[i],.]=x1;
               fitv[rorder[i]]=fitv1;
            endif;
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   xs=x[maxindc(fitv),.];
   retp(xs);
endp;