function MI = mi(cb, p)

%Fitness(cb,p)     -  Returns the conditional prob of a given codebook cb with priors p.
%
%  The fitness calculated is a lower bound on the mutual information assuming 
%     multiplicative covariant Gaussian noise.
%  The codebook 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. 


M = size(cb, 1);
N = size(cb, 2);

global C
global SIGMASQUARED
global MAXF
global MINF

f = 0;
K = zeros(M,N,N);
Kinv = zeros(M,N,N);

for s = 1:M
  K(s,:,:) = SIGMASQUARED*(C*cb(s,:)'*cb(s,:) + (1-C)*diag(cb(s,:).*cb(s,:)));
  Knumer = diag((N*C + 1 - C) * ones(N,1)) - C;
  Kdenom = cb(s,:)' * cb(s,:) * SIGMASQUARED * (1-C) * (N*C + 1 - C);
  Kinv(s,:,:) = Knumer ./ Kdenom;
end

condProb = zeros(M,MAXF-MINF+1,MAXF-MINF+1);
respRange = [MINF:MAXF];
piFact = (2*pi)^(N/2);

for s=1:M
  Kinvs = squeeze(Kinv(s,:,:));
  Ks = squeeze(K(s,:,:)); 
  for rind1 = 1:length(respRange)
    for rind2 = 1:length(respRange)
      rDelt = [respRange(rind1) respRange(rind2)] - cb(s,:);
      condProb(s,rind1,rind2) = exp(- .5 * rDelt * Kinvs * rDelt');
    end
  end
  condProb(s,:,:) = condProb(s,:,:) ./ (piFact*sqrt(det(Ks)));
  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

% now 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);

% now 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

xentR = 0.;
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));
      end
    end
  end
end

MI = xentR + entR;



