new; cls;
gr=1;         /* Gray code(1 or nonzero) or binary(0)  */
size=100;
fn df(x)=(1/2)^(x-1); /* distance function */
length=31;
ngen=100;
npc=3;
pm=0.05; pinv=0; pscr=0; ptlo=0;
ctype=3;   /* 1: one-point crossover 2: n-point crossover 3: uniform crossover           */
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=CellularBGA(gr,size,length,ngen,npc,pm,pinv,pscr,ptlo,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 CellularBGA(gr,size,length,ngen,npc,pm,pinv,pscr,ptlo,ctype,lstype,range,dim,&df,&f);
   local size1,x,fitv,minfitv,i,j,k,di,rorder,i1,i2,i3,inc,inew,seq,p,cp,ks,f:proc;
   local min,max,x01,ii,xs,xsc,xs1,xs2,l,xc,x1,x2,fitv0,fitvc,fitv1,fitv2,y,xp;
   local mask,li,sli,rs,e,actinc,mdist,rii,u,v,w,cdata,df:proc;
   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=round(rndu(size1,length*dim));
   if gr==0;                           /* binary code */
      xp=x;
   else;
      xp=(cumsumc(x[.,1:length]')%2)'; /*  Gray code  */
      ii=2;
      do while ii<=dim;
         xp=xp~(cumsumc(x[.,(ii-1)*length+1:ii*length]')%2)'; 
         ii=ii+1;
      endo;
   endif;
   x01=xp[.,1:length]*((1/2)^seqa(1,1,length));
   ii=2;
   do while ii<=dim;
      x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/2)^seqa(1,1,length));
      ii=ii+1;
   endo;
   xs=range[1,.]+(range[2,.]-range[1,.]).*x01;
   fitv=f(xs);
   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;
         xc=x[rorder[i],.];
         if ctype==1;    /* one-point crossover */
            l=floor((length*dim-1)*rndu(1,1))+1; 
            x1=x[rorder[i],1:l]~x[ks,l+1:length*dim];
            x2=x[ks,1:l]~x[rorder[i],l+1:length*dim];
         elseif ctype==2; /* n-point crossover */
            l=rankindx((length*dim-1)*rndu(length*dim-1,1),1);
            l=l[1:npc];
            l=0|sortc(l,1)|length*dim;
            x1=zeros(1,length*dim);
            x2=zeros(1,length*dim);
            k=1;
            do while k<=npc+1;
               if k%2==1;
                  x1[.,l[k]+1:l[k+1]]=x[rorder[i],l[k]+1:l[k+1]];
                  x2[.,l[k]+1:l[k+1]]=x[ks,l[k]+1:l[k+1]];
               else;
                  x1[.,l[k]+1:l[k+1]]=x[ks,l[k]+1:l[k+1]];
                  x2[.,l[k]+1:l[k+1]]=x[rorder[i],l[k]+1:l[k+1]];
               endif;
               k=k+1;
            endo;
         elseif ctype==3; /* uniform crossover */
            x1=zeros(1,length*dim);
            x2=zeros(1,length*dim);
            mask=round(rndu(length*dim,1));
            k=1;
            do while k<=length*dim;
               if mask[k]==1;
                  x1[.,k]=x[rorder[i],k];
                  x2[.,k]=x[ks,k];
               else;
                  x1[.,k]=x[ks,k];
                  x2[.,k]=x[rorder[i],k];
               endif;
               k=k+1;
            endo;
         endif;
         if gr==0;                             /* binary code */
            xp=xc;
         else;
            xp=(cumsumc(xc[.,1:length]')%2)';  /*  Gray code  */
            ii=2;
            do while ii<=dim;
               xp=xp~(cumsumc(xc[.,(ii-1)*length+1:ii*length]')%2)'; /*  Gray code  */
               ii=ii+1;
            endo;
         endif;
         x01=xp[.,1:length]*((1/2)^seqa(1,1,length));
         ii=2;
         do while ii<=dim;
            x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/2)^seqa(1,1,length));
            ii=ii+1;
         endo;
         xsc=range[1,.]+(range[2,.]-range[1,.]).*x01;
         if gr==0;                             /* binary code */
            xp=x1;
         else;
            xp=(cumsumc(x1[.,1:length]')%2)';  /*  Gray code  */
            ii=2;
            do while ii<=dim;
               xp=xp~(cumsumc(x1[.,(ii-1)*length+1:ii*length]')%2)'; /*  Gray code  */
               ii=ii+1;
            endo;
         endif;
         x01=xp[.,1:length]*((1/2)^seqa(1,1,length));
         ii=2;
         do while ii<=dim;
            x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/2)^seqa(1,1,length));
            ii=ii+1;
         endo;
         xs1=range[1,.]+(range[2,.]-range[1,.]).*x01;
         if gr==0;                             /* binary code */
            xp=x2;
         else;
            xp=(cumsumc(x2[.,1:length]')%2)';  /*  Gray code  */
            ii=2;
            do while ii<=dim;
               xp=xp~(cumsumc(x2[.,(ii-1)*length+1:ii*length]')%2)'; /*  Gray code  */
               ii=ii+1;
            endo;
         endif;
         x01=xp[.,1:length]*((1/2)^seqa(1,1,length));
         ii=2;
         do while ii<=dim;
            x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/2)^seqa(1,1,length));
            ii=ii+1;
         endo;
         xs2=range[1,.]+(range[2,.]-range[1,.]).*x01;
         fitvc=f(xsc);
         fitv1=f(xs1);
         fitv2=f(xs2);
         if fitvc<maxc(fitv1|fitv2);
            if fitv1>fitv2;
               fitvc=fitv1;
               fitv1=fitv1-minc(minfitv|0);
               x[rorder[i],.]=x1;
               fitv[rorder[i]]=fitv1;
            else;
               fitvc=fitv2;
               fitv2=fitv2-minc(minfitv|0);
               x[rorder[i],.]=x2;
               fitv[rorder[i]]=fitv2;
            endif;
         endif;
         /* mutation */
         y=(rndu(1,length*dim).<pm);
         x1=(x[rorder[i],.]+y)%2;
         /* inversion */
         if rndu(1,1)<pinv;
            li=floor(length*dim*rndu(2,1))+1;
            li=sortc(li,1);
            x1[.,li[1]:li[2]]=rev(x1[.,li[1]:li[2]]')';
         endif;
         /* scramble */
         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);
            x1[.,li[1]:li[2]]=x1[.,sli[rankindx(rndu(li[2]-li[1]+1,1),1)]];
         endif;
         /* translocation */
         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;
            x1[.,rs:rs+li[2]-li[1]]=x1[.,li[1]:li[2]];
         endif;
         if gr==0;                             /* binary code */
            xp=x1;
         else;
            xp=(cumsumc(x1[.,1:length]')%2)';  /*  Gray code  */
            ii=2;
            do while ii<=dim;
               xp=xp~(cumsumc(x1[.,(ii-1)*length+1:ii*length]')%2)'; /*  Gray code  */
               ii=ii+1;
            endo;
         endif;
         x01=xp[.,1:length]*((1/2)^seqa(1,1,length));
         ii=2;
         do while ii<=dim;
            x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/2)^seqa(1,1,length));
            ii=ii+1;
         endo;
         xs1=range[1,.]+(range[2,.]-range[1,.]).*x01;
         fitv1=f(xs1);
         if fitvc<fitv1;             
            fitv1=fitv1-minc(minfitv|0);
            x[i,.]=x1;
            fitv[i]=fitv1;
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   if gr==0;                           /* binary code */
      xp=x;
   else;
      xp=(cumsumc(x[.,1:length]')%2)'; /*  Gray code  */
      ii=2;
      do while ii<=dim;
         xp=xp~(cumsumc(x[.,(ii-1)*length+1:ii*length]')%2)'; 
         ii=ii+1;
      endo;
   endif;
   x01=xp[.,1:length]*((1/2)^seqa(1,1,length));
   ii=2;
   do while ii<=dim;
      x01=x01~xp[.,(ii-1)*length+1:ii*length]*((1/2)^seqa(1,1,length));
      ii=ii+1;
   endo;
   xs=range[1,.]+(range[2,.]-range[1,.]).*x01;
   xs=xs[maxindc(fitv),.];
   retp(xs);
endp;