function pc = gh_packed_pc(x,fmmparam)
% GH_PACKED_PC - Pack P and C for the Gauss-Hermite transformation
%
% Syntax:
% pc = GH_PACKED_PC(x,fmmparam)
%
% In:
% x - Evaluation point
% fmmparam - Array of handles and parameters to form the functions.
%
% Out:
% pc - Output values
%
% Description:
% Packs the integrals that need to be evaluated in nice function form to
% ease the evaluation. Evaluates P = (f-fm)(f-fm)' and C = (x-m)(f-fm)'.
% Copyright (c) 2010 Hartikainen, Särkkä, Solin
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.
%%
f = fmmparam{1};
m = fmmparam{2};
fm = fmmparam{3};
if length(fmmparam) >= 4
param = fmmparam{4};
end
if ischar(f) || strcmp(class(f),'function_handle')
if ~exist('param','var')
F = feval(f,x);
else
F = feval(f,x,param);
end
elseif isnumeric(f)
F = f*x;
else
if ~exist('param','var')
F = f(x);
else
F = f(x,param);
end
end
d = size(x,1);
s = size(F,1);
% Compute P = (f-fm)(f-fm)' and C = (x-m)(f-fm)'
% and form array of [vec(P):vec(C)]
pc = zeros(s^2+d*s,size(F,2));
P = zeros(s,s);
C = zeros(d,s);
for k=1:size(F,2)
for j=1:s
for i=1:s
P(i,j) = (F(i,k)-fm(i)) * (F(j,k) - fm(j));
end
for i=1:d
C(i,j) = (x(i,k)-m(i)) * (F(j,k) - fm(j));
end
end
pc(:,k) = [reshape(P,s*s,1);reshape(C,s*d,1)];
end