feedback.m

Nonlinear dynamical factor analysis Matlab package

function [dc_dsm, dc_dsv, dc_dnetm, dc_dnetv, dx] =...
  feedback(x, net, sources, data, noiseparam, status, missing)

% Copyright (C) 2002-2005 Harri Valpola, Antti Honkela and Matti Tornio.
%
% This package comes with ABSOLUTELY NO WARRANTY; for details
% see License.txt in the program package.  This is free software,
% and you are welcome to redistribute it under certain conditions;
% see License.txt for details.

noisevar = normalvar(noiseparam);

[datadim, nsampl] = size(data);
nsources = size(sources, 1);

datavars = noisevar * ones(1, nsampl);

% If constraints are not met, modify the data
if isfield(status, 'constraints'),
  max = sparse(x{4}.e >= repmat(status.constraints(1:datadim,1), 1, nsampl));
  min = sparse(x{4}.e <= repmat(status.constraints(1:datadim,2), 1, nsampl));
  missing(max) = 0;
  missing(min) = 0;
  newdata = max .* repmat(status.constraints(1:datadim,1), 1, nsampl) + ...
	    min .* repmat(status.constraints(1:datadim,2), 1, nsampl); 
  i = (newdata ~= 0);
  data = full(~i .* data + i .* newdata);
%  full(data .* ~missing + missing*inf)

end

% The variances of unknown data samples are set to infinity to effectively
% ignore these samples
datavars(find(missing)) = inf;

dx{4}.var = .5 ./ datavars;
dx{4}.e = (x{4}.e - data) ./ datavars;
dx{4}.extra = dx{4}.var;
dx{4}.multi = repmat(shiftdim(sources.var, -1), [datadim 1 1])...
    .* repmat(reshape(dx{4}.var, [datadim 1 nsampl]), [1 nsources 1])...
    .* (2 * x{4}.multi);

multivar = zeros(size(sources));

% The first layer (linear)
temp = x{4}.multi.^2;

% Somewhat more efficient way to calculate
%   multivar(:,i) = temp(:,:,i)' * multivar(:,i);
for i=1:nsources
  multivar(i,:) = sum(reshape(temp(:,i,:), size(dx{4}.var)) .* dx{4}.var, 1);
end

dx{3}.var = net.w2.var' * dx{4}.extra;
dx{3}.e = net.w2.e' * dx{4}.e + (2*net.w2.var' * dx{4}.extra) .* x{3}.e;
dx{3}.extra = net.w2.e' .^2 * dx{4}.extra;

%dx{3}.multi = zeros(size(x{3}.multi));
%  dx{3}.multi(:,:,i) = net.w2.e' * dx{4}.multi(:,:,i);

d0 = size(net.w2, 2);
[d1 d2 d3] = size(dx{4}.multi);
dx{3}.multi = ...
    reshape(net.w2.e' * reshape(dx{4}.multi, [d1 d2*d3]), [d0 d2 d3]);


[dc_dnetm.w2, dc_dnetv.w2, dc_dnetm.b2, dc_dnetv.b2] = ...
    netgrads(x{3}, dx{4}, net.w2, net.b2);

% The second layer (nonlinear)
if strcmp(status.approximation, 'hermite'),
  [dx{2}.e, dx{2}.var, dx{2}.multi, dx{2}.extra] = ...
      feedback_hermite(dx{3}.e, dx{3}.var, dx{3}.multi, dx{3}.extra, ...
                       x{2}.e, x{2}.var, x{2}.multi, x{2}.extra, ...
                       x{3}.e, x{3}.var, net.nonlin, x{5}, status);
elseif strcmp(status.approximation, 'taylor'),
  [dx{2}.e, dx{2}.var, dx{2}.multi, dx{2}.extra] = ...
      feedback_taylor(dx{3}.e, dx{3}.var, dx{3}.multi, dx{3}.extra, ...
                      x{2}.e, x{2}.var, x{2}.multi, x{2}.extra, ...
                      net.nonlin);
elseif strcmp(status.approximation, 'adaptive'),
  [dx{2}.e, dx{2}.var, dx{2}.multi, dx{2}.extra] = ...
      feedback_adaptive(dx{3}.e, dx{3}.var, dx{3}.multi, dx{3}.extra, ...
                        x{2}.e, x{2}.var, x{2}.multi, x{2}.extra, ...
                        x{3}.e, x{3}.var, net.nonlin, x{5}, ...
			x{3}.lowvar, x{3}.highvar, status);

else
  error('Unsupported approximation')
end

dx{1}.e = net.w1.e' * dx{2}.e + ...
    (2 * net.w1.var' * (dx{2}.var + dx{2}.extra)) .* x{1}.e;
dx{1}.var = (net.w1.e'.^2 + net.w1.var') * dx{2}.var + ...
    net.w1.var' * dx{2}.extra;

[dc_dnetm.w1, dc_dnetv.w1, dc_dnetm.b1, dc_dnetv.b1] = ...
    netgradstop(x{1}, dx{2}, net.w1, net.b1);

dc_dsm = dx{1}.e;
dc_dsv = dx{1}.var + multivar;