function f = FitnessFunc(cb, p)

%Fitness(cb,p)     -  Returns the fitness 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

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

for s = 2:M
  for t = 1:(s-1)
    %  Each unique symbol pair contributes this much.
    weight = p(s) * p(t);
    cbDiff = cb(s,:) - cb(t,:);
    Ks = squeeze(K(s,:,:));          Kt = squeeze(K(t,:,:));
    Kinvs = squeeze(Kinv(s,:,:));    Kinvt = squeeze(Kinv(t,:,:));
    KLsymm1 = sum(sum(Ks .* Kinvt)) + sum(sum(Kinvs .* Kt));
    KLsymm2 = cbDiff * (Kinvs + Kinvt) * cbDiff';
    f = f + weight * (.5*(KLsymm1 + KLsymm2) - 1); 
  end
end   



