function app = rbfapprox(limx, N)
% Create normalized RBF approximator
% Arguments:
% limx: a matrix with two columns and P rows, where P is the number of input variables, containing
% the limits of the domain. limx(d, 1) is the inferior limit and limx(d, 2) is the upper limit for
% the domain on dimension d.
% N: a vector with P elements, with N(d) giving the number of basis function centers on dimension d.
% 
% Outputs: app, an RBF approximator structure.
% To use the approximator:
% app.phi(app, x) returns a column vector of basis function values at point x.
% app.eval(app, x, theta) returns the value taken by the approximator at point x given parameter
% theta. This is just app.phi(app, x)' * theta.

if size(limx, 2) ~= 2 || size(limx, 1) == 0, error('Incorrect limx structure'); end;
app.p = size(limx, 1);
if length(N) ~= app.p, error('Must specify number of basis functions for each dimension'); end;

% computer centers and radii
for d = 1:app.p
    C{d} = limx(d, 1):(limx(d, 2)-limx(d, 1))/(N(d)-1):limx(d, 2);
    rad(d) = (limx(d, 2)-limx(d, 1))/(N(d)-1);
end;

% compute centers, # of RBFs, and radii
app.C = C;
c = [];
for i = 1:length(C{2}),
    c = [c [C{1}; ones(1, length(C{1}))*C{2}(i)]];
end;
app.c = c;
app.N = size(app.c, 2);
app.rad = repmat(rad(:), 1, app.N);

% core functions
app.phi = @rbfdisc_phi; app.eval = @rbfdisc_eval;

end     % TRIANG constructor

% ===== Core functions
% BF values in x
function phi = rbfdisc_phi(app, x)
a = exp(-sum( ((repmat(x, 1, app.N)  - app.c) ./ app.rad).^2, 1 ));
phi = a(:)/sum(a);
end
% Evaluate at a given parameter vector
function v = rbfdisc_eval(app, x, theta)
v = rbfdisc_phi(app, x)' * theta;
end

% ==== Other local functions
function X = flat(x)
% Gives an explicit enumeration of the set of points in the cross product of n grids
%   X = FLAT(x)
% Parameters:
%   x       - cell array of grids
% Returns:
%   X       - flat representation of the grid elements. Each grid element is explicitly stored
%       on a column of X. X has dimension n x N where n is the number of dimensions and N the
%       total number of points on the grids.

% quick return when there is nothing to flatten
if length(x) == 1, X = x{1}; 
else
    n = length(x);
    xnd = ndgridx(x);
    X = zeros(n, numel(xnd{1}));
    for i = 1:n,
        X(i, :) = xnd{i}(:);
    end;
end;

end

function X = ndgridx(x)
% Equivalent of NDGRID for which input and output are cell arrays
% instead of lists of arguments. Useful for working with variable-dimension grids.
%   X = NDGRIDX(x)
% Parameters:
%   x       - cell array of single-dimensional generating grids.
%   X       - cell array of resulting matrices
% See ndgrid().
% 
% ACKNOWLEDGMENT: This file uses portions of code from the file ndgrid.m, developed by MathWorks.

n = length(x);
for i=length(x):-1:1,
  x{i} = full(x{i}); % Make sure everything is full
  siz(i) = numel(x{i});
end
if length(siz)<n, siz = [siz ones(1,n-length(siz))]; end

X = cell(1,n);
for i=1:n,
  xi = x{i}(:); % Extract and reshape as a vector.
  s = siz; s(i) = []; % Remove i-th dimension
  xi = reshape(xi(:,ones(1,prod(s))),[length(xi) s]); % Expand
  X{i} = permute(xi,[2:i 1 i+1:n]); % Permute to i'th dimension
end
end