/* Fuzzy regression mixed with FCM */
new; cls;
x={1,2,3,5,6,8,9};
y={2,3,2,4,3,5,6};
alpha=0.5; /* 0<=alpha<=1 */
nc=2;
w=2;
n=rows(x);
xx1=ones(n,1)~x;

m=10;
lambda=100;
ngen=10000;
dim=n*(nc-1);
s=rndu(dim,1);
range=zeros(1,dim)|ones(1,dim);
ntimes=10;
yhats=wfregmcMOCMAESm(alpha,nc,w,y,xx1,m,lambda,ngen,s,range,dim,ntimes,&f);

library pgraph;
graphset;
_plctrl=-1|zeros(nc,1);
_pltype=6;
xy(x,y~yhats);

/*
**  fregfcm.txt - Fuzzy regression mixed with FCM.
** (C) Copyright 2011 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Estimates fuzzy regression coefficients considering weighted FCM in terms of common membership matrix.
**              We mix fuzzy regression with fuzzy c-means with the ratio of alpha and 1-alpha.
**
**  Format:     yhats=wfregmcMOCMAESm(alpha,nc,w,y,xx1,m,lambda,ngen,s,range,dim,ntimes,&f);
**
**  Input:      alpha    scalar,    weight or ratio of fuzzy regression over fuzzy c-means (0<=alpha<=1)
**
**              nc       scalar,    number of cluster (2,3,...,nc)
**
**              w        scalar,    parameter of power (1<w)
**
**              y        vector,    n dependent variables
**
**              x1       matrix,    n x k independent variables
**
**              m        scalar,    number of parents (1,2,...,m)
**
**              lambda   scalar,    number of children (1,2,...,lambda)
**
**              ngen     scalar,    number of generations (1,2,...,ngen)
**
**              s        vector,    dim standard deviations (0<s)
**
**              range    matrix,    2 x dim matrix consisting of zeros for the first row and ones for the second row(default)
**
**              dim      scalar,    n*(nc-1) (default)
**
**              times    scalar,    number of times to repeat the calculation (1,2,...,times)
**
**              &f       function,  objective function to be maximized
**
**  Output:     yhats    matrix,    n x nc matrix of resulting estimates of y
**
**  Notice:     It takes lots of memory to run. Try to use a small number of m for Light version of GAUSS.
**
*/

proc f(x);
   local times,U,Uwj,bj,Sj1,cj,dj,Sj2,sum,out,i,j,label,yj,oldval;
   times=cols(x);
   out=zeros(times,1);
   i=1;
   do while i<=times;
      U=zeros(nc,n);
      j=1;
      do while j<=nc-1;
         U[j,.]=x[n*(j-1)+1:n*j,i]';
         j=j+1;
      endo;
      j=2;
      do while j<=nc-1;
         U[j,.]=U[j,.].*(ones(1,n)-sumc(U[1:j-1,.])');
         j=j+1;
      endo;
      U[nc,.]=ones(1,n)-sumc(U[1:nc-1,.])';
      U=U';
      sum=0;
      j=1;
      do while j<=nc;
         Uwj=diagrv(zeros(n,n),U[.,j]^w);
         oldval=trapchk(1);
         trap 1,1;
         yj=inv(xx1'Uwj*xx1);
         trap oldval,1;
         if scalerr(yj);
            out[i]=-1e+256;
            goto label;
         endif;
         bj=inv(xx1'Uwj*xx1)*xx1'Uwj*y;
         if bj==miss(zeros(rows(bj),1),0);
            out[i]=-1e+256;
            goto label;
         endif;
         Sj1=(y-xx1*bj)'Uwj*(y-xx1*bj);
         cj=sumc((U[.,j]^w).*xx1)./sumc(U[.,j]^w);
         dj=sqrt(sumc(((xx1-cj')^2)'));
         Sj2=dj'Uwj*dj;
         sum=sum+alpha*Sj1+(1-alpha)*Sj2;
         j=j+1;
      endo;
      out[i]=-sum;
    label:
      i=i+1;
   endo;
   retp(out);
endp;

proc wfregmcMOCMAESm(alpha,nc,w,y,xx1,m,lambda,ngen,s,range,dim,ntimes,&f);
   local xs,f0,i,xi,fi,n,U,j,yhats,Uwj,bj,Sj1,cj,dj,Sj2,sum,f:proc;
   xs=0;
   f0=-1e+256;
   i=1;
   do while i<=ntimes;
      xi=wfregmcMOCMAES(alpha,nc,w,y,xx1,m,lambda,ngen,s,range,dim,&f);
      fi=f(xi);
      if f0<fi;
         xs=xi;
         f0=fi;
      endif;
      i=i+1;
   endo;
   n=rows(xx1);
   U=zeros(nc,n);
   j=1;
   do while j<=nc-1;
      U[j,.]=xs[n*(j-1)+1:n*j]';
      j=j+1;
   endo;
   j=2;
   do while j<=nc-1;
      U[j,.]=U[j,.].*(ones(1,n)-sumc(U[1:j-1,.])');
      j=j+1;
   endo;
   U[nc,.]=ones(1,n)-sumc(U[1:nc-1,.])';
   U=U';
   print "U=" U;
   yhats=zeros(n,nc);
   sum=0;
   j=1;
   do while j<=nc;
      Uwj=diagrv(zeros(n,n),U[.,j]^w);
      bj=inv(xx1'Uwj*xx1)*xx1'Uwj*y;
      print " bj=" bj';
      yhats[.,j]=xx1*bj;
      Sj1=(y-xx1*bj)'Uwj*(y-xx1*bj);
      cj=sumc((U[.,j]^w).*xx1)./sumc(U[.,j]^w);
      dj=sqrt(sumc(((xx1-cj')^2)'));
      Sj2=dj'Uwj*dj;
      sum=sum+alpha*Sj1+(1-alpha)*Sj2;
      j=j+1;
   endo;
   print "sum=" sum;
   print "yhats=" yhats;
   retp(yhats);
endp;

proc wfregmcMOCMAES(alpha,nc,w,y,xx1,m,lambda,ngen,s,range,dim,&f);
   local d,ptarget,cp,cc,ccov,pthresh,pcs,ss,Cs,pc,C,x,x_1,x1,i,j,k,pbar,p,count,fxk,fx1,fx,f:proc;
   d=1+dim/(2*lambda);
   ptarget=1/(5+sqrt(lambda)/2);
   cp=ptarget*lambda/(2+ptarget*lambda);
   cc=2/(dim+2);
   ccov=2/(dim^2+6);
   pthresh=0.44;
   pcs=zeros(m*dim,1);
   if rows(s)==1;
      ss=s*ones(dim,1);
   else;
      ss=s;
   endif;
   Cs=eye(dim);
   i=2;
   do while i<=m;
      ss=ss|s;
      Cs=Cs|eye(dim);
      i=i+1;
   endo;
   x=range[1,.]+(range[2,.]-range[1,.]).*rndu(m,dim);
   pbar=0;
   count=0;
   k=1;
   do while k<=m;
      pc=pcs[(k-1)*dim+1:k*dim];
      s=ss[(k-1)*dim+1:k*dim];
      C=Cs[(k-1)*dim+1:k*dim,.];
      j=1;
      do while j<=ngen;
         x1=rndn(lambda,dim);
         if dim<lambda;
            x1=x[k,.]+(x1-meanc(x1)')*inv(chol(vcx(x1-meanc(x1)')))*(s'.*chol(corrvc(C)));
         else;
            x1=x[k,.]+x1*(s'.*chol(corrvc(C)));
         endif;
         i=1;
         do while i<=lambda;
            x1[i,.]=maxc(range[1,.]|x1[i,.])';
            x1[i,.]=minc(x1[i,.]|range[2,.])';
            i=i+1;
         endo;
         count=sumc(f(x[k,.]').<f(x1'));
         x_1=x[k,.];
         if count>=1;
            x[k,.]=x1[maxindc(f(x1')),.];
         endif;
         p=count/lambda;
         pbar=(1-cp)*pbar+cp*p;
         if pbar<pthresh;
            pc=(1-cc)*pc+sqrt(cc*(2-cc))*((x[k,.]'-x_1')./s);
            C=(1-ccov)*C+ccov*pc*pc';
         else;
            pc=(1-cc)*pc;
            C=(1-ccov)*C+ccov*(pc*pc'+cc*(2-cc)*C);
         endif;
         C=(C+C')/2;
         s=s*exp(1/d*((pbar-ptarget)/(1-ptarget)));
         trap 1;
         if s<1e-10 or scalerr(inv(chol(vcx(x1-meanc(x1)')))) or scalerr(chol(corrvc(C)));
            break;
         endif;
         j=j+1;
      endo;
      k=k+1;
   endo;
   x=x[maxindc(f(x')),.];
   retp(x');
endp;