% program for INDCLUS, according to Chaturvedi & Carroll (1994)
%
% min sum(k) SS( S_k -P W_k P')
% 
% P binary matrix (overlap allowed); W_k diagonal
%
% input: S=(S(1)...S(K)) matrices;
%	 r: dimensionality of solution
%
%
%	solution: W (rows containing W(k)'s) and P (n x r)
%
% 	start: 0 = default; 1 = random; 2 = current values 
 
[n,mn]=size(S);
m=mn/n;
conv=[];
func=[];func2=[];dysm=[];
conv=input(' convergence crit. (def. 1e-6): ');
if isempty(conv),conv=1e-6;end;
method=input('1 = original algorithm; 2= elementwise algorithm ');
start=input(' 0 = default; 1 = random; 2 = current values ');
if isempty(start),start=0;end;
if start~=1,nanal=1;end;
if start==1
  nanal=input(' How many starts do you want to use ');
  ss=input(' Which number of random start to start with');
end;
offdia=input(' Only offdiagonal ? 1=Yes (default), 0 = No  ');
addit=input('Fit additive term too ? 1=Yes (default), 0=No  ');
if isempty(offdia),offdia=1;end;
if isempty(addit),addit=1;end;


if offdia==1, disp('Offdiagonal fitting !');end;
if addit==1,disp('Additive terms fitted as well !');end;
if method==1,disp('Original algorithm');end;
if method==2,disp('Elementwise algorithm');end;

Csol=[];Psol=[];Wsol=[];Qsol=[];IT=[];
for ana=1:nanal

% initialize P
 if start==0
   U=zeros(n,n);
   for k=1:m
     U=U+S(:,(k-1)*n+1:k*n); 
   end;
   [p,d,q]=svd(U,0);
   X=p(:,1:r);
   P=X>mean(mean(X));
   W=ones(m,r);
 end;
 if start==1
   rand('seed',(ss+ana-1)^2)
   P=round(rand(n,r));
   W=rand(m,r);
 end;
 Q=P;	% splitting 
 
 P=logical(P);
 Q=logical(Q);

 I=eye(n);
 if addit==0,c=zeros(1,m);end;
 if addit==1
   c=[];
   for k=1:m
     if offdia==1
       c(k)=sum(sum( (S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*Q').*(1-I) ))/sum(sum(1-I));
     else
       c(k)=mean(mean(S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*Q'));
     end;
   end;
 end;			% additive terms
 D1=1-[ones(m,1)*I(:)'];		% indicates offdiag terms in T1
 D2=1-[kron(ones(1,m),I)];	% indicates offdiag terms in T2=SS and T3

 ff=0;
 for k=1:m
    ff=ff+ssq(S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*P');
 end;
 f=ff;
 fprintf('\n Function value at start is %12.8f \n',f);
 fold=f+2*conv*f;
 iter=0;
 dysm(ana)=0;
 while fold-f>conv*f | iter<3
  fold=f;
  iter=iter+1;
  for i=1:r
    v=ones(1,r);
    v=logical(v);
    v(i)=0;SS=[];
    if r==1 
      for k=1:m
        SS=[SS S(:,(k-1)*n+1:k*n)-c(k)];
      end;
    end;
    if r>1 
      for k=1:m
        SS=[SS S(:,(k-1)*n+1:k*n)-P(:,v)*diag(W(k,v))*Q(:,v)'-c(k)];
      end;
    end;

  % updating W(:,i)

    T1=[];
    for k=1:m
      Sk=SS(:,(k-1)*n+1:k*n);
      T1=[T1;Sk(:)'];
    end;

    % correct Pi or Qi in case column is zero: replace it by 1
    while sum(P(:,i))==0
      P(:,i)=round(rand(n,1)); 
      disp('correction of P');
      if method==2,Q(:,i)=P(:,i);end;
    end;
    while sum(Q(:,i))==0
      Q(:,i)=round(rand(n,1));
      disp('correction of Q');
    end;

    if offdia==1
      W(:,i)=T1/((kron(Q(:,i),P(:,i)))'.*D1(1,:));
    else
      W(:,i)=T1/(kron(Q(:,i),P(:,i)))';
    end;

    % set negative elements to 0
    v=W(:,i)<0;
    if length(v)==1,
      if v==1,W(1,i)=0;end;
    end;
    if length(v)>1,W(v,i)=zeros(sum(v),1);end; 


   if method==1	 % original SINDCLUS algorithm 
 
  % updating P(:,i)
    
    if offdia==1
      V1=sum(((SS-ones(n,1)*kron(W(:,i),Q(:,i))').^2)'.*D2')';  
      V2=sum(((SS).^2)'.*D2')';
    else
      V1=sum(((SS-ones(n,1)*kron(W(:,i),Q(:,i))').^2)')';
      V2=sum(((SS).^2)')';
    end;
    P(:,i)=V1<V2;	% choose 1 if V1<V2 ; 0 otherwise

  % updating Q(:,i)
    
    T3=[];
    for k=1:m
      Sk=SS(:,(k-1)*n+1:k*n);
      T3=[T3 Sk'];
    end;
    if offdia==1
      V1=sum(((T3-ones(n,1)*kron(W(:,i),P(:,i))').^2)'.*D2')';
      V2=sum(((T3).^2)'.*D2')';
    else
      V1=sum(((T3-ones(n,1)*kron(W(:,i),P(:,i))').^2)')';
      V2=sum(((T3).^2)')';
    end;
    Q(:,i)=V1<V2;	% choose 1 if V1<V2 ; 0 otherwise

   end;

   if method==2	 % elementwise  update  
 
  % updating P(:,i) and Q(:,i) simultaneously
    
    if sum(W(:,i))==0
      P(:,i)=round(rand(n,1));
      disp('arbitrary update Pi because of zero Wi');
    end;	
	% in this case any P(:,i) will do
    if sum(W(:,i))>0
     if offdia==1  
       G=sum(W(:,i).^2)*(ones(n,n)-eye(n));
       for k=1:m
         G=G-2*W(k,i)*(SS(:,(k-1)*n+1:k*n).*(ones(n,n)-eye(n)));
       end;
     else 
       G=sum(W(:,i).^2)*ones(n,n);
       for k=1:m
         G=G-2*W(k,i)*SS(:,(k-1)*n+1:k*n);
       end;
     end;
     for j=1:n
       P(j,i)=0.5;	% Pi with element j replaced by 1
       P(j,i)=G(j,:)*P(:,i)<0;
     end;
    end; 
    Q(:,i)=P(:,i);
   end;	% end of "if method==2"

  end; % end if "for i=.." loop 

 % updating additive terms

  if addit==1
    for k=1:m
      if offdia==1
        c(k)=sum(sum( (S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*Q').*(1-I) ))/sum(sum(1-I));
      else
        c(k)=mean(mean(S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*Q'));
      end;
    end;
  end;

  f=0;f2=0;
  for k=1:m
    if offdia==1
      f=f+ssq((S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*Q'-c(k)).*(1-I));
      if method==1,f2=f2+ssq((S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*P'-c(k)).*(1-I));end;
    else
      f=f+ssq(S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*Q'-c(k));
      if method==1,f2=f2+ssq(S(:,(k-1)*n+1:k*n)-P*diag(W(k,:))*P'-c(k));end;
    end;
  end;
  if fold<f,dysm(ana)=1;end;
  fprintf(' Function after %g iters is: %12.8f (f_symm = %12.8f ) \n',iter,f,f2);
end;    
func(ana)=f;
func2(ana)=f2;
IT(ana)=iter;
Csol(:,ana)=c';
Psol(:,ana)=P(:);
Qsol(:,ana)=Q(:);
Wsol(:,ana)=W(:);
end;
if method==1,disp('Original algorithm');end;
if method==2,disp('Elementwise algorithm');end;
if offdia==1, disp('Offdiagonal fitting !');end;
if addit==1,disp('Additive terms fitted as well !');end;
disp(' Minimal function value is ');
[mm,mi]=min(func)
c=Csol(:,mi)';
P(:)=Psol(:,mi)
Q(:)=Qsol(:,mi)
W(:)=Wsol(:,mi)


