feedfw.m

Nonlinear dynamical factor analysis Matlab package

function x = feedfw(s0, net, approximation)
% FEEDFW  Do feedforward
%
%    Usage:
%      x = feedfw(s, net, approximation)
%      where x will be a cell array with intermediate results
%      from several phases of computation with x{4} containing the
%      final result.

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

% Identity map ignores feedforward
if isfield(net, 'identity') & net.identity,
  x{4} = s0;
  return;
end

% Default approximation is hermite
if nargin < 3,
  approximation = 'hermite';
end

s = probdist(s0);
x = cell(1, 4);

sdim = size(s, 1);
tdim = size(s, 2);
hdim = size(net.w1, 1);
odim = size(net.w2, 1);

% Initialisation:         x{1} = s
x{1}.e = s.e;
x{1}.var = s.var;
x{1}.multi = repmat(eye(sdim), [1 1 tdim]);

% First linear layer:     x{2} = A*x{1} + a
x{2}.e = net.w1.e * x{1}.e + repmat(net.b1.e, [1, tdim]);
x{2}.extra = net.w1.var * (x{1}.e.^2 + x{1}.var) + ...
    repmat(net.b1.var, [1, tdim]);
x{2}.multi = reshape(net.w1.e * reshape(x{1}.multi, [sdim sdim*tdim]), ...
		     [hdim sdim tdim]);
x{2}.var = x{2}.extra + net.w1.e.^2 * x{1}.var;

% Nonlinear hidden layer: x{3} = phi(x{2})
if strcmp(approximation, 'hermite'),
  [x{3}, x{5}] = nonlin_hermite(x{2}, net.nonlin);
elseif strcmp(approximation, 'taylor'),
  x{3} = nonlin_taylor(x{2}, net.nonlin);
elseif strcmp(approximation, 'adaptive'),
  [x{3}, x{5}] = nonlin_adaptive(x{2}, net.nonlin);
elseif strcmp(approximation, 'none'),
  x{3} = x{2};
  x{3}.e = feval(net.nonlin, x{2}.e);
else
  error('Unsupported approximation')
end

% Second linear layer:    x{4} = B*x{3} + b
x{4}.e = net.w2.e * x{3}.e + repmat(net.b2.e, [1, tdim]);
x{4}.extra = net.w2.var * x{3}.e.^2 + ...
    net.w2.e.^2 * x{3}.extra + net.w2.var * x{3}.var + ...
    repmat(net.b2.var, [1, tdim]);
x{4}.multi = reshape(net.w2.e * reshape(x{3}.multi, [hdim sdim*tdim]), ...
		     [odim sdim tdim]);
x{4}.var = x{4}.extra + ...
    reshape(sum(repmat(shiftdim(x{1}.var, -1), [odim 1 1]) ...
		.* (x{4}.multi .^ 2), 2), [odim tdim]);