%% Linear Regression
global beta
global regr

Fmatrix = feval(regr,points(1:npoints,:));
beta = Fmatrix \ ptsvals(1:npoints,:);

ptsvals_regr = ptsvals(1:npoints,:) - Fmatrix*beta;

%% RBFN approximation construction
istep = istep+1;
rbfrad = (dmax*(Dim*npoints)^(-1/Dim))^2;

rbfs = zeros(npoints,npoints);
% for i=1:npoints
%     for j=1:npoints
%         delta = sum((points(i,:)-points(j,:)).^2);
%         rbfs(i,j) = exp( -delta/rbfrad);
%     end
% end

for i=1:npoints
    delta = (sum((repmat(points(i,:),npoints,1) - points(1:npoints,:)).^2,2))';
    rbfs(i,:) = exp( -delta./rbfrad);
end

for i=1:npoints
    rbfs(i,i) = rbfs(i,i) + lambda;
end

global weights
% weights = linsolve(rbfs,ptsvals(1:npoints));
weights = linsolve(rbfs,ptsvals_regr);
clear rbfs

rbfs = zeros(npoints,npoints);

% for i=1:npoints
%     for j=1:npoints
%         delta = sum((points(i,:)-points(j,:)).^2);
%         rbfs(i,j) = exp( -delta/rbfrad);
%     end
% end

for i=1:npoints
    delta = (sum((repmat(points(i,:),npoints,1) - points(1:npoints,:)).^2,2))';
    rbfs(i,:) = exp( -delta./rbfrad);
end

%% author: Anicka Kucerova, Czech Technical University in Prague
%% version 0.2.0
%% 17 Sep 2013