// P76MixConDis.sce
// Mixture estimation - static components and pointer model
// - 14 continuous + 2 disctere  variables 
// - random initialization (prepared for expert knowledge)
// - display by marginals - select v1,v2 (variables) and z (discrete)
[u,t,n]=file();                         // find working directory
chdir(dirname(n(2)));                   // set working directory
clear("u","t","n")                      // clear auxiliary data
exec("ScIntro.sce",-1),mode(0)          // intro to sesion
rand('seed',0)

function mixGr(v1,v2,z)
  // plot and print of clusters
  // v1,v2 - variables for marginal
  // z     - value of discrete variable
  disp('Centers of selected marginals - all components')
  set(scf(100),'position',[850 50 600 400])
  for i=1:nc
    m=psi2row([z,i],[nz,nc]);
    if ~isempty(yC(m))
      plot(yC(m)(v1,:),yC(m)(v2,:),tx(m),'markersize',3)
      disp(Est.Cy(m).th([v1,v2]),'Cent '+string(m)+', (z = '+string(z)+')')
    end
  end
  title('Variables '+string(v1)+', '+string(v2)+'  and  zt = '+string(z))
endfunction
// .............................................................................

nd=1350;            // number of data: all = 1350
nc=5;               // number of components
I_estCov=0;         // estimation of noise covariances ? 0|1 no|yes

// DATA 
load _data/d_con.dat d_con;         // HERE, THE DATA ARE LOADED
load _data/d_ord.dat d_ord;
load _data/d_nom.dat d_nom;

// data assignment (time runs in rows)
[dsc,mD,sD]=scal([d_con';d_ord']); // scaling
// by scalling we transform data to standard range -> init. by rand is possible
yt=dsc;                       // scaled continuous and ordinal data
yu=unscal(yt,mD,sD);          // unscaled data
[nv,nd]=size(yt);
dd=d_nom;
bn=max(dd,'r');               // nominal data
zz=dd';
for i=1:nd
  zt(1,i)=psi2row(zz(:,i),bn); // coding 2var -> 1var
end
nz=max(zt);                   // number of values of nominal variable

// MODEL INITIALIZATION
// random generation of initial centers
for j=1:nz     // components are indexed by z and c
for k=1:nc
  i=psi2row([j,k],[nz,nc]);
  Est.Cy(i).th=rand(1,nv,'u'); // HERE SET THE INITIAL CENTERS OF COMPONENTS
  // for good clustering it is necessary to set better initial centers
  Est.Cy(i).sd=0.01*eye(nv,nv); // covariances of components
  Est.Cy(i).cv=Est.Cy(i).sd**2;
  Ps=[Est.Cy(i).th 1]';       // initial parameters
  Est.Cy(i).V=Ps*Ps';         // statistics
end
end
Est.ka=ones(1,nc);            // counter
// initial discrete pars
Est.Cz.V=rand(nz,nc,'u')+1;
Est.Cz.th=fnorm(Est.Cz.V,2);
// initial pointer probabilities
Est.Cp.V=ones(1,nc);          // pointer statistics
Est.Cp.th=fnorm(ones(1,nc));  // pointer parameter
Est.ct(1)=1;                  // initial poiner value
w=fnorm(ones(1,nc));          // weights
tht=list(); for i=1:nc, tht(i)=[]; end
wt=zeros(nc*nz,nd);
ii=zeros(nc,nd);

// ESTIMATION 
printf(' ')
kk=ceil(nd/10); printf('\n  2 4 6 8 |\n ')
for t=1:nd                    // time loop
if t/kk==fix(t/kk), printf('.'); end
  // proximities continuous and ordinal
  for k=1:nc
    i=psi2row([k,zt(t)],[nc,nz]); 
    [xxx,G(k)]=GaussN(yt(:,t),Est.Cy(i).th,Est.Cy(i).cv); // likelihood
    G(k)=G(k)*Est.Cz.th(zt(t));
    ii(k,t)=i;
  end
  Lq=G-max(G);
  q=exp(Lq(:)');
  
  ww=q'.*Est.Cp.th'; w=ww/sum(ww);           // generation of weights
  wt(ii(:,t),t)=w;

  // Update of statistic 
  Ps=[yt(:,t)' 1];                           // extended reg.vec.
  for k=1:nc
    i=psi2row([k,zt(t)],[nc,nz]);            // selected components (zt)
    Est.Cy(i).V=Est.Cy(i).V+w(k)*Ps'*Ps;     // information matrix
    Est.ka(k)=Est.ka(k)+w(k);                // counter
    Est.Cz.V(zt(t),k)=Est.Cz.V(zt(t),k)+w(k);            // discrete stat
    Est.Cp.V(k)=Est.Cp.V(k)+w(k);            // pointer statistics

    // Point estimates
    Vyy=Est.Cy(i).V(1:nv,1:nv);              // part Vyy - y.y'
    Vy=Est.Cy(i).V($,1:nv);                  // part Vy  - psi.y
    V1=Est.Cy(i).V($,$);                     // part V1  - psi.psi'
    Est.Cy(i).th=inv(V1+1e-8*eye(V1))*Vy;    // pt.est. - reg.coef.
    if t>nd+1, I_estCov=1; end               // DELAYED ESTIMATION OF COV
    if I_estCov~=0
      // pt.est. of noise covariance - used or not
      Est.Cy(i).cv=(Vyy-Vy'*inv(V1+1e-8*eye(V1))*Vy)/Est.ka(k); 
    end
  tht(k)=[tht(k) Est.Cy(k).th'];
  end
  Est.Cz.th(zt(t),:)=fnorm(Est.Cz.V(zt(t),:)); // nominal variables
      
  Est.Cp.th=fnorm(Est.Cp.V,2);     // pt.est. of pointer parameter
  [ss,Est.ct(1,t)]=max(w);         // store of c
end
[xxx,cE]=max(wt,'r');
EstO=Est;
save estO.dat EstO;                // save for repetitive running

// RESULTS
disp ' '; disp(vals(Est.ct),'Frequencies of pointer values')

// plot of weights (last 100 weights)
set(scf(1),'position',[50 50 400 300]);
plot(wt(:,($-100):$)')
title('Weights - check of functionality')

// plot of of evolution of parameter estimates
if 0                // activate by 1
for i=1:nc
  scf(10+i);
  plot(tht(i)')
  title('comp '+string(i))
end
end

// Detection of centers of components in original scale
// - for each z, nc-points from 14-dimensional space is stored
// - type Th + Enter to see the centers
Th=list();
for j=1:nz
  zz=[];
  for k=1:nc
    i=psi2row([j,k],[nz,nc]);
    zp=unscal(Est.Cy(i).th',mD,sD)';
    zz=[zz; zp];
  end
  jz=row2psi(j,bn);
  Cz(jz(1),jz(2)).c=zz;
  Th(j)=zz;
end

tx1=['bo';'ro';'go';'mo';'ko';'co';'yo'];
tx2=['bx';'rx';'gx';'mx';'kx';'cx';'yx'];
tx3=['b+';'r+';'g+';'m+';'k+';'c+';'y+'];
tx=[tx1;tx2;tx3];

yC=list();
for k=1:nz
  for i=1:nc
    m=psi2row([k,i],[nz,nc]);
    j=find((cE==i)&(zt==k));
    yC(m)=yt(:,j);
  end
end

v1=1; v2=4; z=1;        // set marginal variables and discrete z
mixGr(v1,v2,z)

// Parameters of all nc*nz components (for scaled data)
set(scf(111),'position',[200 10 600 600]);
title('Components -> down, discrete variable -> right')
for i=1:nc*nz
  subplot(nc,nz,i)
  bar(Est.Cy(i).th)
end

  
  
