Computer Science COMP 199 (Winter term)
Excursions in Computing Science

---

%MATLABpak09cI:
%
%	histogArith.m				mutInfo.m
%
% histogArith.m			THM				090824
% function hist = histogArith(hist1,op,hist2)
% arithmetic operations on histograms (discrete absolute distributions)
% e.g.	  1   2   3   2   1
%     +	 -2  -1   0   1   2
%   1 -1 -3  -2  -1   0   1
%   1  0 -2  -1   0   1   2  =>	-3  -2  -1   0   1   2   3	hist(1,k)
%   1  1 -1   0   1   2   3	 1   3   6   7   6   3   1	hist(2,k)
% these are not equi-step histograms
% e.g.	  1   1   1
%     +	 -1   0   1
%   1 -4 -5  -4  -3
%   1  0 -1   0   1  =>	-5  -4  -3  -1   0   1   3   4   5    NB missing -2   2
%   1  4  3   4   5	 1   1   1   1   1   1   1   1   1		  0   0
% usage: print and plot (histogArith plots bighist)
% hist1 = [-1,0,1;1,1,1]
% hist = histogArith(hist1,'+',hist1)
% hold on
% hist = histogArith(hist1,'+',hist)
%First attemt
% histograms conceptually indexed by range[from,by,to]: for k=1:1:(to-from)/by+1
%function [hist,range] = histogArith(hist1,range1,op,hist2,range2)
%if range1(1,2)==range2(1,2)		% initial implementation: by1 = by2
%  switch op
%  case {'+','-'}
%    
%  otherwise disp('histogArith ',op,' not implemented')
%  end %switch op
%else disp('histogArith implementation requires range1(2)==range2(2): ',
%	  range1(1,2), range2(1,2))
%end %if range1(1,2)==range2(1,2)
function hist = histogArith(hist1,op,hist2)
% find range for bighist:
n1 = min(hist1(1,:));
n2 = min(hist2(1,:));
x1 = max(hist1(1,:));
x2 = max(hist2(1,:));
switch op
case '+'
  from = n1 + n2;
  to   = x1 + x2;
case '-'
  from = n1 - x2;
  to   = x1 - n2;
case '*'
  xtremes = [n1*n2,n1*x2,x1*n2,x1*x2];
  from = min(xtremes);
  to   = max(xtremes);
case '/'    % careful: bighist below needs integers; not guaran even after mult
  xtremes = [n1*n2,n1*x2,x1*n2,x1*x2];		% same size as mult:will expand
  from = min(xtremes);
  to   = max(xtremes);
case {'x','n'}	% mxx, min
  from = min(n1,n2);
  to = max(x1,x2);
otherwise
  ['histogArith operator ' op ' not implemented']
% quit	NEVER USE quit OR exit EXCEPT TO TERMINATE MATLAB SESSION
% I don't know how to terminate _my_ program, so this will fail on   bighist =
end %switch op
sz1 = size(hist1);
sz2 = size(hist2);
bighist = zeros(1,to-from+1);
for k1 = 1:sz1(2)
  a1 = hist1(1,k1);
  a2 = hist1(2,k1);
  for k2 = 1:sz2(2)
    b1 = hist2(1,k2);
    b2 = hist2(2,k2);
    switch op
    case '+'
      ab = a1 + b1;
    case '-'
      ab = a1 - b1;
    case '*'
      ab = a1 * b1;
    case '/'
      ab = a1 / b1;
    case 'x'
      ab = max(a1,b1);
    case 'n'
      ab = min(a1,b1);
    end %switch op
    bighist(1,ab-from+1) = bighist(1,ab-from+1) + a2*b2;
  end %for k2 = 1:sz2(2)
end %for k1 = 1:sz1(2)
plot([from:to],bighist)
hist = [];
kk = 0;
for k = 1:to-from+1
  if bighist(1,k)~=0,
    histvec = [k+from-1; bighist(1,k)];
    hist = [hist histvec];
  end %if bighist(1,k)~=0 
end %for k = 1:to-from+1

% function [mI,tot] = mutInfo(M)		THM			120107
%	see ../matlab/flowMutInfo.m
% calculate mutual information for unnormalized flow matrix M
% uses log2z()
% (don't confuse with mutInfo(), which assumes normalized distribution)
function [mI,tot] = mutInfo(M)
marg1 = sum(M);
marg2 = sum(M')';
tot = sum(marg1);
expanded = [[M;marg1],[marg2;tot]]	% remove ; to display matrix with sums
pjk = M/tot;
pj = marg2/tot;
pk = marg1/tot;
Igpjk = -sum(sum(pjk.*log2z(pjk)));
Ignpj = -sum(pj.*log2z(pj));
Ignpk = -sum(pk.*log2z(pk));
mI = Ignpj + Ignpk - Igpjk;