function DensityComponents()
%
% this is a skeleton script that illustrates the algorithm
% and the gradient computations for Maximum Likelihood learning
% of Density Components, as described in Karklin and Lewicki (2005).
%
%
% yan karklin. 2006.


I = 399;  % dimensionality of the data (e.g. image patches)
J = 100;  % number of higher-order basis functions
N = 300;  % number of data samples in a single batch

epsv = 0.1;  % step size for v (might want to adjust so it tapers)
epsB = 0.01; % step size for B

% initialize B
B = .1 * randn(I,J);

% main loop - go up the gradient of the likelihood 
for i=1:1000,

  % get new batch of 300 image patches here (resample x from images)
  % W is learned separately using standard ICA.
  % x are whitened image patches; all u's should have variance 1 
  %     over the entire data ensemble
  x = getData();  % <--- resampling code not included!
  u = W * x; 
  
  % initialize v_MAP;  
  vmap = .1*randn(J,N);

  % use gradient descent to obtain MAP estimate \hat{v}
  for j=1:30,
    dv = .5 * B'*(sqrt(2)*abs(u)./ exp(B*vmap/2) - 1);
    dv = dv - sqrt(2)*sign(vmap);
    vmap = vmap + epsv .* dv;
  end;

  % update B
  dB = .5 * (sqrt(2)*abs(u)./ exp(B*vmap/2) - 1) * vmap';
  B = B + epsB .* dB;

end;




function L = calcL(B, u, v)
% computes the negative log-likelihood
% -logP(u) = -logP(u|B,v) - logP(B) - logP(v)

N = size(u,2);
Bv  = B * v;

% conditionally Laplacian prior P(u|B,v)
% with variance  sigma^2 = exp(Bv);  
lPu = sum(sum(Bv/2 + (sqrt(2)*abs(u)) ./ exp(Bv/2)));

% assume Laplacian prior on v
lPv = sqrt(2)*sum(sum(abs(v)));

L = lPu/N + lPv/N;