function [MI, miDeriv] = miPois(cb, p)

%miPois(cb,p) -  Returns the MI of a given codebook cb with priors p,
%      and the derivatives of MI wrt cb.
%
%  Mutual information is calculated assuming independent Poisson noise.
%
%  cb is an N by M matrix, where N is the number of neurons and 
%     M is the number of symbols.
%  p is a vector of length M. 
%
% NOTE: This code assumes N=2

M = size(cb, 1);
N = size(cb, 2);
global MAXF
global MINF

condProb = zeros(M,MAXF-MINF+1,MAXF-MINF+1);
respRange = [MINF:MAXF];

% conditional probability computed for every response vector
%   assumes noise for each neuron is conditionally independent
%   also normalizes cond prob over the response vectors
for s=1:M
  for rind1 = 1:length(respRange)
    for rind2 = 1:length(respRange)
      condProb(s,rind1,rind2) = ...
          prod(poisspdf([respRange(rind1) respRange(rind2)], cb(s,:)));
    end
  end
  condProb(s,:,:) = condProb(s,:,:) ./ sum(sum(condProb(s,:,:)));
end

% clf;
% hold on;
% for s=1:M
%   condProbs = squeeze(condProb(s,:,:));
%   surface(MINF:MAXF,MINF:MAXF,condProbs);
% end

% compute p(r) = sum_s p(r|s) p(s)
probR = zeros(MAXF-MINF+1,MAXF-MINF+1);
for rind1 = 1:length(respRange)
  for rind2 = 1:length(respRange)
    for s=1:M
      probR(rind1,rind2) = probR(rind1,rind2) + condProb(s,rind1,rind2)*p(s);
    end
  end
end
probR(:,:) = probR(:,:) ./ sum(sum(probR(:,:)));


%figure;
%clf;
%surface(MINF:MAXF,MINF:MAXF,probR);

% compute entropy of p(r)
entR = 0.;
for rind1 = 1:length(respRange)
  for rind2 = 1:length(respRange)
    if probR(rind1,rind2)
      entR = entR - probR(rind1,rind2) * log(probR(rind1,rind2));
    end
  end
end

% compute cross-entropy, and mean rate for each symbol
xentR = 0.;
meanR = zeros(M,2);
for rind1 = 1:length(respRange)
  for rind2 = 1:length(respRange)
    for s = 1:M
      if condProb(s,rind1,rind2)
        xentR = xentR + condProb(s,rind1,rind2) * p(s) * ...
	        log(condProb(s,rind1,rind2));
        meanR(s,:) = meanR(s,:) + ...
            [respRange(rind1) respRange(rind2)] * condProb(s,rind1,rind2);
      end
    end
  end
end

MI = xentR + entR;

% compute partials of MI wrt cb
miDeriv = zeros(s,2);
for rind1 = 1:length(respRange)
  for rind2 = 1:length(respRange)
    for s = 1:M
      if condProb(s,rind1,rind2)
        miDeriv(s,:) = miDeriv(s,:) + ...
            ([respRange(rind1) respRange(rind2)] - meanR(s,:)) * ...
            condProb(s,rind1,rind2) * ...
            log(condProb(s,rind1,rind2)/probR(rind1,rind2));
      end
    end
  end
end
for s = 1:M
  miDeriv(s,:) = miDeriv(s,:) ./ cb(s,:);
end



