
xrange = -5:0.2:5;
[Xin,Yin] =  meshgrid(xrange);      % create grid of test points
response_inputs = [Xin(:) Yin(:)];  % convert to N x 2 matrix
[Nresponse, dummy]  = size(response_inputs);

% create training data
N = 300;
Xtrain = randn(N,2)*1.5;
Ytrain = Xtrain(:,1).^2 + Xtrain(:,2)*2; %+ randn(N,1)*0.01;

% sinusoidal training set
% Xtrain = randn(N,2)*1.5;
% Ytrain = sin(Xtrain(:,1)) + 2*cos(Xtrain(:,2));

% input transformations
% choose transformation type
% 1 - linear, no transformation
% 2 - quadratic expansion
% 3 - exponential terms (X exp(abs(X))
% 4 - frequency transformation

trans_type = 3;
% calculate response to training data
X = transform(Xtrain, trans_type);
% estimate optimal weighting matrix
W = pinv(X)*Ytrain;

% calculate output response
X = transform(response_inputs, trans_type);
response_output = X*W;

mesh(Xin,Yin,reshape(response_output,sqrt(Nresponse),sqrt(Nresponse)));
hold on
plot3(Xtrain(:,1),Xtrain(:,2),Ytrain,'o');
hold off

nin = length(X(1,:));
figure(2);
for i = 1:nin
   subplot(ceil(sqrt(nin)),ceil(sqrt(nin)),i);
   mesh(Xin,Yin,reshape(X(:,i),sqrt(Nresponse),sqrt(Nresponse)));
end