% create training data
N = 300;
[Xin,Yin] =  meshgrid(0.1:0.1:.9);      % create grid of training points
Xtrain = [Xin(:) Yin(:)];
Ytrain = Xtrain(:,1).^2 + Xtrain(:,2)*2; %+ randn(N,1)*0.01;
%[Xtrain, Ytrain] = getdata;  % read in classification data


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


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

% set up the RBF centres
centre_range = 0:0.15:1;
[centresX, centresY] = meshgrid(centre_range);
centres = [centresX(:) centresY(:)];
widths = 0.01*ones(size(centres));

% calculate response to training data
X = rbftransform(Xtrain, centres, widths);

% estimate optimal weighting matrix
W = pinv(X)*Ytrain;

% calculate output response
X = rbftransform(response_inputs, centres, widths);
response_output = X*W;

figure(1)
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);
clf
for i = 1:nin
   mesh(Xin,Yin,reshape(X(:,i),sqrt(Nresponse),sqrt(Nresponse)));
   hold on
end
hold off


