function [MI, miDeriv] = miPoisSamp(cb, p, rSamples)

%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

global MAXF
global MINF

M = size(cb, 1); % number of symbols
N = size(cb, 2); % number of neurons
S = size(rSamples,1);

% conditional probability computed for every response vector
%   assumes noise for each neuron is conditionally independent
%   also normalizes cond prob over the response vectors
condProb = zeros(M,S);
for s=1:M
    for samp=1:S
        condProb(s,samp) = prod(poisspdf(rSamples(samp,:), cb(s,:)));
    end
    condProb(s,:) = condProb(s,:) ./ sum(condProb(s,:));
end

% clf;
% hold on;
% s=1;
% plot3(rSamples(:,1), rSamples(:,2), condProb(s,:)','bd');

% compute p(r) = sum_s p(r|s) p(s)
probR = zeros(S,1);
for samp=1:S
    for s=1:M
        probR(samp) = probR(samp) + condProb(s,samp)*p(s);
    end
end
probR = probR ./ sum(probR);

% compute entropy of p(r)
entR = 0.;
for samp=1:S
    if probR(samp)
      entR = entR - probR(samp) * log(probR(samp));
    end
end

% compute cross-entropy, and mean rate for each symbol
xentR = 0.;
meanR = zeros(M,N);
for samp=1:S
    for s = 1:M
      if condProb(s,samp)
        xentR = xentR + condProb(s,samp) * p(s) * ...
	        log(condProb(s,samp));
        meanR(s,:) = meanR(s,:) + ...
	        rSamples(samp,:) * condProb(s,samp);
      end
    end
end

MI = xentR + entR;

% compute partials of MI wrt cb
miDeriv = zeros(M,N);
for samp=1:S
    for s = 1:M
      if condProb(s,samp)
        miDeriv(s,:) = miDeriv(s,:) + ...
            (rSamples(samp,:) - meanR(s,:)) * ...
            condProb(s,samp) * ...
            log(condProb(s,samp)/probR(samp));
      end
    end
end
for s = 1:M
  miDeriv(s,:) = miDeriv(s,:) ./ cb(s,:);
end



