%
% cwms.m
% Cluster-Weighted State Analysis
% (c) Neil Gershenfeld 1/28/97
%
clear all
%
min_var = 1e-6;				% minimum cluster variance
dim = 2;				% number of measurements per record
data_file_name='testdata.dat';	% input data file
%
% read the data file
%    lines with comments start with spaces or %
%    lines with data have the format
%       record state data(1) ... data(dim) optional_comment
% 
data_file = fopen(data_file_name,'r');
point = 0;
while (1)
   c = fscanf(data_file,'%c',1);
   if feof(data_file)
      break
      end
   position = ftell(data_file);
   fseek(data_file,(position-1),'bof');
   if ((c == '%') | (c == ' ') | (c == 13))
      fgetl(data_file);  
   else
      point = point + 1;
      record(point) = fscanf(data_file,'%d',1);
      states(point) = fscanf(data_file,'%d',1);
      data(point,1:dim) = (fscanf(data_file,'%f',dim))';
      fgetl(data_file);
      end
   end
fclose(data_file);
npts = point;
sort_states = sort(states);
unique_states = states([1 1+find(sort_states(1:(npts-1)) ~= ...
   sort_states(2:npts))]);
nstates = length(unique_states);
fprintf('   read %d points in %dD with %d states from %s\n', ...
   npts, dim, nstates, data_file_name);
%
% make data zero mean unit variance
%
%data = data - kron(ones(npts,1),mean(data));
%data = data ./ kron(ones(npts,1),std(data));
%
% get input info
%
nclusters = input('Number of clusters? ');
nsteps = input('Number of iterations? ');
momentum = input('Momentum? ');
%
% generate random starting clusters with constant weights
%
pc = ones(1,nclusters)/nclusters;
mu = kron(ones(nclusters,1),(max(data)-min(data))) ...
     .* rand(nclusters,dim) + kron(ones(nclusters,1),min(data));
dmu = zeros(nclusters,dim);
var = kron(ones(nclusters,1),(max(data)-min(data)).^2);
var_err = [];
var_mistakes = [];
meanvar_mistakes = [];
mean_mistakes = [];
covar = [];
covar_err = [];
for cluster = 1:nclusters
   covar = [covar; diag((max(data)-min(data)).^2)];
   end
ps = ones(nstates,nclusters)/nstates;
%
% E-M iteration
%
step = 1;
while (step <= nsteps)
   %
   % plot the first two axes for data and clusters
   %
   for cluster = 1:nclusters
      plot(mu(cluster,1),mu(cluster,2),'.');
      hold on
      [value state] = max(ps(:,cluster));
      h=text(mu(cluster,1),mu(cluster,2),sprintf('%d',state), ...
         'FontSize',30,'FontWeight','bold','Color','Red');
      [vec,val] = eig(covar((1:dim)+dim*(cluster-1),:));
      for d = 1:dim
         %
	     % plot using covariances
	     %
         plot([mu(cluster,1)-sqrt(val(d,d))*vec(1,d) ...
	       mu(cluster,1)+sqrt(val(d,d))*vec(1,d)], ...
              [mu(cluster,2)-sqrt(val(d,d))*vec(2,d) ...
	       mu(cluster,2)+sqrt(val(d,d))*vec(2,d)]);
 	    end
	 %
	 % plot using variances
	 %
%      plot([mu(cluster,1)-sqrt(var(cluster,1)) ...
%         mu(cluster,1)+sqrt(var(cluster,1))], ...
%         [mu(cluster,2) mu(cluster,2)]);
%      plot([mu(cluster,1) mu(cluster,1)], ...
%         [mu(cluster,2)-sqrt(var(cluster,2)) ...
%         mu(cluster,2)+sqrt(var(cluster,2))]);
      end
   for state = 1:nstates
      index = find(states == state);
      plot(data(index,1),data(index,2),'.');
      text(data(index,1),data(index,2),sprintf('%d',state),'FontSize',12);
      end
   grid off
   hold off
   drawnow
   %
   % calculate forward probability
   %
   for cluster = 1:nclusters
      %
      % px = p(x_n|c_m)
      % p = p(s_n,x_n,c_m)
      %
      % calculate using covariances
      %
      cinv = inv(covar((1:dim)+dim*(cluster-1),:));
      px(:,cluster) = (sqrt(det(cinv))/((2*pi)^(dim/2))) ...
         * exp(-dot((data-kron(ones(npts,1),mu(cluster,:)))', ...
	     (cinv*(data-kron(ones(npts,1),mu(cluster,:)))'))'/2);
      p(:,cluster) = ps(states,cluster) .* px(:,cluster) ...
         .* kron(ones(npts,1),pc(cluster));
      %
      % calculate using variances
      %
      pxvar(:,cluster) = prod((exp(-(data-kron(ones(npts,1),mu(cluster,:))).^2 ...
         ./ (2*kron(ones(npts,1),var(cluster,:)))) ...
         ./ kron(ones(npts,1),sqrt(2*pi*var(cluster,:))))')';
      pvar(:,cluster) = ps(states,cluster) .* pxvar(:,cluster) ...
         .* kron(ones(npts,1),pc(cluster));
      end
   %
   % calculate posterior probability
   % pp = p(c_m|s_n,x_n)
   %
   pp = p ./ kron(sum(p')',ones(1,nclusters));
   %
   % update cluster weights
   % pc = p(c_m)
   %
   pc = sum(pp) / npts;
   %
   % update input parameters
   %
   for cluster = 1:nclusters
      munew = momentum * dmu(cluster,:) + ...
         sum(data .* kron(pp(:,cluster),ones(1,dim))) ...
	     / (npts * pc(cluster));
      dmu(cluster,:) = munew - mu(cluster,:);
      mu(cluster,:) = munew;
      var(cluster,:) = min_var + sum((data - ...
         kron(ones(npts,1),mu(cluster,:))).^2 ...
         .* kron(pp(:,cluster),ones(1,dim))) / (npts * pc(cluster));
      for d = 1:dim
         covar(d+(cluster-1)*dim,:) = sum( ...
	    (data - kron(ones(npts,1),mu(cluster,:))) .* ...
	    (kron(data(:,d),ones(1,dim)) - ones(npts,dim)*mu(cluster,d)) ...
            .* kron(pp(:,cluster),ones(1,dim))) / (npts * pc(cluster));
         end
      covar(((1:dim)+(cluster-1)*dim),:) = covar(((1:dim)+(cluster-1)*dim),:) ...
         + min_var*diag(ones(1,dim));
      end
   %
   % update output parameters
   %
   for state = 1:nstates
      index = find(states == state);
      ps(state,:) = sum(p(index,:)) ./ sum(p);
      end 
   %
   % calculate covariance sensitivity
   % dpx = partial_x p(x_n,c_m))
   % dp = partial_x p(s_n,x_n,c_m)
   % dps = partial_x p(s_n|x_n)
   %
   dpx = zeros(npts,dim);
   dp = zeros(npts,dim);
   for cluster = 1:nclusters
      cinv = inv(covar((1:dim)+dim*(cluster-1),:));
	  dpx = dpx - px(:,cluster*ones(1,dim)) .* ...
		    (pc(cluster)*ones(npts,dim)) ...
	        .* (cinv*(data-kron(ones(npts,1),mu(cluster,:)))')';
      dp = dp - ps(states,cluster*ones(1,dim)) .* ...
		    px(:,cluster*ones(1,dim)) .* (pc(cluster)*ones(npts,dim)) ...
	        .* (cinv*(data-kron(ones(npts,1),mu(cluster,:)))')';
      end
   dps = (kron(sum(px')',ones(1,dim)) .* dp ...
      - kron(sum(p')',ones(1,dim)) .* dpx) ...
      ./ kron((sum(px')').^2,ones(1,dim));
   %
   % calculate variance confidence per axis
   % ppredvar = p(s_n|x_{nd})
   %
   for d = 1:dim
      for cluster = 1:nclusters
         pxvard(:,cluster) = exp(-(data(:,d)-mu(cluster,d)).^2 ...
            / (2*var(cluster,d))) ...
            / sqrt(2*pi*var(cluster,d));
         end
      for state = 1:nstates
         ppredvar(:,state) = zeros(npts,1);
         for cluster = 1:nclusters
            ppredvar(:,state) = ppredvar(:,state) + ...
               kron(ones(npts,1),ps(state,cluster)) .* pxvard(:,cluster) ...
	           .* kron(ones(npts,1),pc(cluster));
            end
         end
      ppredvar = ppredvar ./ kron(sum(ppredvar')',ones(1,nstates));
	  varconfs(:,d) = diag(ppredvar(:,states));
      varconf(d) = (nstates*sum(diag(ppredvar(:,states))) - npts) ...
	     /(npts*(nstates-1));
	  end
   %
   % predict states for data set
   %
   % using covariance PDF
   % ppred = p(s|x_n)
   % covar_err = N - sum_n p(s_n|x_n)
   %
   for state = 1:nstates
      ppred(:,state) = zeros(npts,1);
      for cluster = 1:nclusters
         ppred(:,state) = ppred(:,state) + ...
            kron(ones(npts,1),ps(state,cluster)) .* px(:,cluster) ...
	    .* kron(ones(npts,1),pc(cluster));
         end
      end
   ppred = ppred ./ kron(sum(ppred')',ones(1,nstates));
   covar_err = [covar_err npts-sum(diag(ppred(:,states)))];
   %
   % using variance PDF
   %
   for state = 1:nstates
      ppredvar(:,state) = zeros(npts,1);
      for cluster = 1:nclusters
         ppredvar(:,state) = ppredvar(:,state) + ...
            kron(ones(npts,1),ps(state,cluster)) .* pxvar(:,cluster) ...
	    .* kron(ones(npts,1),pc(cluster));
         end
      end
   ppredvar = ppredvar ./ kron(sum(ppredvar')',ones(1,nstates));
   var_err = [var_err npts-sum(diag(ppredvar(:,states)))];
   [value varindex] = max(ppredvar');
   var_mistakes = [var_mistakes length(find(~(states == varindex)))];
   %
   % using means and variances
   %
   distances = [];
   for cluster = 1:nclusters
      distances(:,cluster) = sum(((data - kron(ones(npts,1),mu(cluster,:))).^2 ...
         ./ kron(ones(npts,1),var(cluster,:)))')';
      end
   [value, cluster_state] = max(ps);
   [value, meanvar_index] = min(distances');
   meanvar_mistakes = [meanvar_mistakes ...
      length(find(~(cluster_state(meanvar_index) == states)))];
   %
   % using means
   %
   distances = [];
   for cluster = 1:nclusters
      distances(:,cluster) = sum(((data - kron(ones(npts,1),mu(cluster,:))).^2)')';
      end
   [value, mean_index] = min(distances');
   mean_mistakes = [mean_mistakes ...
      length(find(~(cluster_state(mean_index) == states)))];
   %
   % print and plot errors
   %
   fprintf('step %d  covar conf: %.1f var conf: %.1f err: %i  mean/var err %i:  mean err: %i\n', ...
      step, covar_err(step), var_err(step), var_mistakes(step), meanvar_mistakes(step), ...
      mean_mistakes(step));
   if (step == nsteps)
      input('CR to plot convergence ');
      plot(1:step,covar_err,1:step,var_err,1:step,var_mistakes,1:step,meanvar_mistakes, ...
         1:step,mean_mistakes)
	  input('CR to plot sensitivity  ');
	  hist(abs(dps));
	  input('CR to plot axis dependence  ');
	  hist(varconfs);
	  fprintf('axis: ');
	  for d = 1:dim
	     fprintf('%f ',varconf(d));
		 end
      fprintf('\n');
      answer = input('Total number of iterations? (CR to end) ');
      if ~isempty(answer)
         nsteps = answer;
         end
      end
   step = step + 1;
   end % iterate step loop
clear state;
