splinend(n, varargin)

function zi = splinend(n, varargin)
% Complete argument lists:
% function zi = splinend(n, x1, x2, ..., xn, z, xi1, xi2, ..., xin, ...
%                        action, options)
% Tensorial spline interpolation in n-dimensional
%
% Short calling:
% function s = splinend(n, x1, ..., xn, z, [],...,[],'InitOnly') OR
%          zi = splinend(n, x1, x2, ..., xn, z, xi1, xi2, ..., xin) OR
%          zi = splinend(..., 'ForceReset') OR
%          zi = splinend(..., s)
%
%          zi = splinend(..., options); % fifth parameters
%
% input
%    n: integer, dimension
%    x1, ...xn: vectors (required to be monotonics), OR
%    they could be the an ND-array of same dimension than z
%                                                (as generated by ndgrid)
%    z: data matrix z(length(x1),...,length(xn))
%    xi1, ..., xin: vectors/arrays.
%       location of the points where the interpolation will performed.
%        a) If they have a same size, they contains the scattered
%                   coordinates of interpolation points
%        b) If they are vectors, their cartesian product grid will be
%                   interpolated on.
%          In case of cartesian product needs to be generated on
%          same-length linear grid data, we only need to transpose one of
%          the input vector {xik} to different orientation.
%    options: optional structure with following field
%      zndgrid: integer flag [0]: ndgrid or meshgrid on output z style for n<2,
%      economy [0], 1: Economy or Recursive spline structure
%      dorder  [zeros(1,n)]: order of the derivative
%      nogrid [false]: reserved field name, not documented 
%      tensorindex: reserved field name, not documented
%
% output zi: result of the interpolation, see remark a) b) on input {xik}
%            for of the size of output zi.
%
% tensorial bi-cubic-spline interpolation (natural conditions by default)
%
% Author: Bruno Luong
% Last update: 19/Feb/2009
%    Call new find_idx with mex dichotomy search for speed
%
% HOW TO CALL:
% - If the last input parameter action is 'InitOnly', splinend returns
% the structure
% "s" that contains initialization from the coarse data (x1, ... xn,z)
% - The structure "s" can be passed then as last argument for successive
% interpolation calls (in this case inputs (x1, ... xn,z) are ignored).
% - The initialization structure is remembered by the function splinend,
% unless until the function is called with 'InitOnly' or 'ForceReset'.
%
%

persistent splinestruct;

if ~isscalar(n)
    error('splinend: The first input should be the dimension');
end

xlist=varargin(1:n);
z=varargin{n+1};
xilist=varargin(n+2:2*n+1);

varargin(end+1:2*n+3)={''};
[action options]=deal(varargin{2*n+2:2*n+3});

actionpos=2*n+3;

if (nargin>=actionpos) && isstruct(action) && isfield(action, 'Tag') && ...
        ischar(action.Tag) && ...
        ~isempty(regexp(action.Tag,'splinend|spline1d', 'once'))
    
    splinestruct=action;
    
elseif (nargin>=actionpos && strcmpi(action, 'InitOnly')) || ...
       (nargin>=actionpos && strcmpi(action, 'ForceReset'))
    
    % format the input vectors
    for d=1:n
        if ~isvector(xlist{d})
            if ndims(xlist{d})==n % extract the vectors and the dimension d
                subs=num2cell(ones(1,n));
                subs{d}=':';
                sr = struct('type','()','subs', {subs});
                xlist{d} = reshape(subsref(xlist{d},sr),[],1);
            else
                error('splinend: input x must be vector or %d-dimensional grid', n);
            end
        else % put them in column
            xlist{d} = xlist{d}(:);
        end
    end
    
    % Check for meshgrid style in dimension <=2
    if n<=2 && getoption(options, 'zndgrid', 0)==0
        z = z.'; % transpose
    end
    
    % Switch to ndgrid mode
    options.zndgrid=1;
    
    % Basic dimensional check
    dimmatch = arrayfun(@(k) length(xlist{k})==size(z,k), 1:n);
    if ~all(dimmatch)
        error('splinend: dimensions not match');
    end
    
    if getoption(options, 'economy', 0)
        % Building tensorial spline structure, recursive, economy size
        splinestruct = spline_tens_ecom(xlist, z, options);
    else
        % Building tensorial spline structure, recursive, large size
        splinestruct = spline_tensorial(xlist, z, options);
    end
    
    griddata = struct('x', {xlist}, ...
        'data', z, ...
        'options', options);
    splinestruct.griddata = griddata;
    
    % Return the structure
    if strcmpi(action, 'InitOnly')
        zi=splinestruct;
        return
    end
     
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate the splinend
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(splinestruct)
    error('Splinend: ForceReset must be switched on');
end

% Check id the psstruct structure is correctly initialized and
% it's compatible
if ~isstruct(splinestruct) || ~isfield(splinestruct,'Tag') || ...
        isempty(regexp(splinestruct.Tag,'splinend|spline1d', 'once'))
    error('splinend: not valid structure');
end

% Make sure they have the same dimension
[xilist{:}] = duplicate(xilist{:});

% Get the derivative order
if isfield(splinestruct,'dorder')
    spline_dorder=splinestruct.dorder;
    % Keep only the n-last (recursive reason)
    spline_dorder=spline_dorder(end-n+1:end);
else
    spline_dorder=zeros(1,n);
end

% Additional derivative order
dorder = getoption(options, 'dorder', zeros(1,n));
options.dorder = dorder;

if getoption(options, 'nogrid', false)
    odim = 0;
else % check if interpolated point are grid
    xigrid = cell(size(xilist));
    [odim xigrid{:}] = isgrid(xilist{:});
end

if isfield(splinestruct,'economy')
    economy=splinestruct.economy;
else
    economy=0;
end

if all(odim>0) % Plaid data
    
    % Retreive grid data    
    if economy && isfield(options,'tensorindex')
        % Compute the dimension indices
        d = length(splinestruct.xgrid)-n+1;
        index = [options.tensorindex options.dflag];
        % Get the whole the n-last dimensions
        index(d:end/2) = NaN;
        % Get the function value only
        index(end/2+d:end) = 1;
        [xlist zi]=splineinfo_nocheck(splinestruct, 'grid', 'data', index);
        % Truncate the first part of the list
        xlist = xlist(d:end);
    else
        % Retreive full grid data
        [xlist zi]=splineinfo_nocheck(splinestruct, 'grid', 'data', 0);
    end
    
    if isfield(splinestruct,'griddata') && ...
       isfield(splinestruct.griddata,'options')
        % Get the options that has been used to build the spline function
        orgopt = splinestruct.griddata.options;
    end
    
    sizez = size(zi);
    % Initialize orgopt if not exist
    orgopt.zndgrid = 1; % ndgrid mode, working dimension is the first
    % Loop on dimension, circular rotation of the current
    % array with successive tensorial-interpolation
    for d=1:n
        orgopt.dorder = options.dorder(d)+spline_dorder(d); % setup derivative
        zi = spline1d(xlist{d}, reshape(zi,sizez(d),[]), ...
                      xigrid{d}, 'ForceReset', orgopt).';
    end
    zi = reshape(zi,cellfun(@length,xigrid));
    if any(odim~=(1:n))
        zi = permute(zi,odim);
    end
    
    return % Done
end

% Not Plaid data -------------------------

% Economy strategy
if economy
    
    % Compute the dimension indice
    d = length(splinestruct.xgrid)-n+1;
    % Get the breaks at the specified dimension
    xgrid=splinestruct.xgrid{d};
    nx=length(xgrid);
    dx=diff(xgrid);
    % Get the interpolation points
    xi = xilist{1};
    
    if n==1 % Last dimension is reached
        
        index = [options.tensorindex options.dflag];
        % Get the whole last dimension & both value and 2nd derivative
        index([d end]) = NaN; 
        zdata = splineinfo_nocheck(splinestruct,'data',index);
        
        s1 = struct('Tag', 'spline1d structure', ...
            'xgrid', reshape(xgrid,1,[]), ...
            'nx', nx, ...
            'dx', reshape(dx,1,[]), ...
            'zz', zdata(:,1), ...
            'zzxx', zdata(:,2), ...
            'dorder', spline_dorder);
        
        zi = spline1d([], [], xi, s1, options);
    else % Recursive call
        
        %
        % Put this flag to 0 (independent searching for index)
        %               or 1 (MATLAB searching for index)
        %
        UsingMatlabInterp1=0;
        
        if UsingMatlabInterp1
            xi=interp1(xgrid,(1:nx),xi,'linear','extrap');
            %xi=min(max(xi,1),nx);
        else
            xi=find_idx(xi, xgrid, 'extrap'); %#ok
        end
        
        %
        idx=max(min(fix(xi),nx-1),1);
        %
        % Basis of the cubic polynome in x-direction evaluate
        % at the points of interpolation
        %
        % sqrdx=dx(idx).^2';
        % bx=(xi-idx);
        % ax=(1-bx);
        % cx=sqrdx.*(ax.^3-ax);
        % dx=sqrdx.*(bx.^3-bx);
        
        h=reshape(dx(idx),size(idx));
        
        switch round(spline_dorder(1)+dorder(1)) % Cumulative Derivative
            case 0, % interpolated spline function
                bx=(xi-idx);
                ax=(1-bx);
                g=-(h.^2).*ax.*bx;
                cx=(ax+1).*g;
                dx=(bx+1).*g;
            case 1, % the first derivative
                b=(xi-idx); a=(1-b);
                ax=-1./h+zeros(size(xi));
                bx=+1./h+zeros(size(xi));
                cx=h.*(+3*b.*(a+1)-2);
                dx=h.*(-3*a.*(b+1)+2);
            case 2, % the second derivative
                b=(xi-idx); a=(1-b);
                ax=zeros(size(xi));
                bx=zeros(size(xi));
                cx=6*a;
                dx=6*b;
            case 3, % the third derivative (discontinuous)
                ax=zeros(size(xi));
                bx=zeros(size(xi));
                cx=-6./h+zeros(size(xi));
                dx=+6./h+zeros(size(xi));
            otherwise
                if dorder(1)>3
                    error('spline1d: 4th derivative is ill-defined (Dirac)');
                else
                    error('spline1d: not valid derivative');
                end
        end
        
        %
        % These linear-index will be used to access to two brakets indices
        % where the points of interpolation belong.
        %
        idx0=idx; % left bracket
        % idx1=idx0+1;
        
        % Allocate memory
        zz0 = zeros(size(xi));
        zz1 = zeros(size(xi));
        zzxx0 = zeros(size(xi));
        zzxx1 = zeros(size(xi));
        
        % Sort the points within their brakets
        % This is usually fater than accumarray
        [isort p]=sort(idx0(:));
        istart=find(diff([0; isort(:); Inf])); % 0 and Inf must be indices
        % that fall outside
        
        % Compute the dimension indice
        d = length(splinestruct.xgrid)-n+1;
        
        if d==1 %~isfield(options,'tensorindex')
            options.tensorindex = nan(1,n);
            options.dflag = ones(1,n);
        end        
        
        % recursive of spline nd
        % Loop on break points
        
        % Set the derivative for next loop on dimension
        options.dorder = options.dorder(2:end);
        for k=1:length(istart)-1
            
            % number of the braket
            j=isort(istart(k));
            % indices of points falling on the breaks;
            idxj = p(istart(k):istart(k+1)-1);
            
            xionbreak = arrayfun(@(d) xilist{d}(idxj), 2:n, ...
                                 'UniformOutput', false);
            
            % left values
            options.tensorindex(d)=j;
            options.dflag(d) = 1;
            zz0(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, splinestruct, options);
            
            % right value
            options.tensorindex(d)=j+1;
            options.dflag(d) = 1;
            zz1(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, splinestruct, options);
            
            % left second derivative
            options.tensorindex(d)=j;
            options.dflag(d) = 2;
            zzxx0(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, splinestruct, options);
            
            % right second derivative
            options.tensorindex(d)=j+1;
            options.dflag(d) = 2;
            zzxx1(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, splinestruct, options);
        end
        
        zi = bsxfun2006(@times,ax,zz0) + ...
             bsxfun2006(@times,bx,zz1) + ...
             bsxfun2006(@times,cx,zzxx0) + ...
             bsxfun2006(@times,dx,zzxx1);
    end
    
else % Not economy
    
    xgrid=splinestruct.xgrid;
    nx=splinestruct.nx;
    dx=splinestruct.dx;
    zz=splinestruct.zz;
    zzxx=splinestruct.zzxx;
    
    if isstruct(zz) && ~isempty(strfind(zz(1).Tag,'spline'))
        
        xi = xilist{1};
        
        %
        % Put this flag to 0 (independent searching for index)
        %               or 1 (MATLAB searching for index)
        %
        UsingMatlabInterp1=0;
        
        if UsingMatlabInterp1
            xi=interp1(xgrid,(1:nx),xi,'linear','extrap');
            %xi=min(max(xi,1),nx);
        else
            xi=find_idx(xi, xgrid, 'extrap'); %#ok
        end
        
        %
        idx=max(min(fix(xi),nx-1),1);
        %
        % Basis of the cubic polynome in x-direction evaluate
        % at the points of interpolation
        %
        % sqrdx=dx(idx).^2';
        % bx=(xi-idx);
        % ax=(1-bx);
        % cx=sqrdx.*(ax.^3-ax);
        % dx=sqrdx.*(bx.^3-bx);
        
        h=reshape(dx(idx),size(idx));
        
        switch round(spline_dorder(1)+dorder(1)) % Cumulative Derivative
            case 0, % interpolated spline function
                bx=(xi-idx);
                ax=(1-bx);
                g=-(h.^2).*ax.*bx;
                cx=(ax+1).*g;
                dx=(bx+1).*g;
            case 1, % the first derivative
                b=(xi-idx); a=(1-b);
                ax=-1./h+zeros(size(xi));
                bx=+1./h+zeros(size(xi));
                cx=h.*(+3*b.*(a+1)-2);
                dx=h.*(-3*a.*(b+1)+2);
            case 2, % the second derivative
                b=(xi-idx); a=(1-b);
                ax=zeros(size(xi));
                bx=zeros(size(xi));
                cx=6*a;
                dx=6*b;
            case 3, % the third derivative (discontinuous)
                ax=zeros(size(xi));
                bx=zeros(size(xi));
                cx=-6./h+zeros(size(xi));
                dx=+6./h+zeros(size(xi));
            otherwise
                if dorder(1)>3
                    error('spline1d: 4th derivative is ill-defined (Dirac)');
                else
                    error('spline1d: not valid derivative');
                end
        end
        
        %
        % These linear-index will be used to access to two brakets indices
        % where the points of interpolation belong.
        %
        idx0=idx; % left bracket
        % idx1=idx0+1;
        
        % Allocate memory
        zz0 = zeros(size(xi));
        zz1 = zeros(size(xi));
        zzxx0 = zeros(size(xi));
        zzxx1 = zeros(size(xi));
        
        % Sort the points within their brakets
        [isort p]=sort(idx0(:));
        istart=find(diff([0; isort(:); Inf])); % 0 and Inf must be indices
        % that fall outside
        
        % recursive of spline nd
        % Loop on break points
        options.dorder=options.dorder(2:end)+spline_dorder(2:end);
        for k=1:length(istart)-1
            
            % number of the braket
            j=isort(istart(k));
            % indices of points falling on the breaks;
            idxj = p(istart(k):istart(k+1)-1);
            
            xionbreak = arrayfun(@(d) xilist{d}(idxj), 2:n, ...
                'UniformOutput', false);
            % left values
            zz0(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, zz(j), options);
            % right value
            zz1(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, zz(j+1), options);
            % left second derivative
            zzxx0(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, zzxx(j), options);
            % right second derivative
            zzxx1(idxj) = splinend(n-1, xlist{2:end}, [], ...
                xionbreak{:}, zzxx(j+1), options);
        end
        
        zi = bsxfun2006(@times,ax,zz0) + ...
             bsxfun2006(@times,bx,zz1) + ...
             bsxfun2006(@times,cx,zzxx0) + ...
             bsxfun2006(@times,dx,zzxx1);
        
    elseif  ~isempty(strfind(splinestruct.Tag,'spline1d')) % 1D spline
        zi = spline1d([], [], xilist{1}, splinestruct, options);
    else % other type of function, not used right now
        switch round(dorder(1))
            case 0,
                zi = evalfct(splinestruct, xilist{1});
            case 1,
                zi = evalder(splinestruct, 1, xilist{1});
            otherwise
                error('splinend: generic derivative not implemented');
        end
    end
    
end % Not economy

return

end

% recursive call of tensorial spline
function s = spline_tensorial(x, z, options)
x1 = x{1};
nx1 = length(x1);
nd = length(x);
s = spline1d(x1, reshape(z, nx1, []), [], 'InitOnly', options);
if nd>1
    s.zz = arrayfun(@(k) spline_tensorial(x(2:end), s.zz(k,:), options), ...
        (1:nx1)');
    s.zzxx = arrayfun(@(k) spline_tensorial(x(2:end), s.zzxx(k,:), options), ...
        (1:nx1)');
    s.Tag = 'splinend structure';
    s.nvar = nd;
end
end

% direct call of tensorial to populate for second derivatives
% in all dimensions
function s = spline_tens_ecom(x, z, options)

n=length(x);
ngrid = cellfun(@length,x);

% flag = dec2bin(0:2^n-1,n);
% flag = double(flag(:,end:-1:1)-'0');

% We have 2^n arrays of the original size
% In the later part, first indice is function value
%                    second indice is second derivative
sizedata = [ngrid 2*ones(1,n)]; % real size
data = zeros([prod(ngrid) 2^n]); % two parts

for k=1:2^n
    d=nextpow2(k);
    if d==0
        % initialized the first array with interpolation values
        data(:,k) = z(:);
    else
        %
        kp = k-2^(d-1);
        cs = circshift((1:n),[0 -(d-1)]);
        % compute the array k from array kp by computing the second
        % derivative of the dtat(:,kp) along the dimension d
        % (using spline1d)
        
        % dimension kp becomes first
        zk = permute(reshape(data(:,kp),ngrid),cs);
        % call spline1d
        sk = spline1d(x{d}, reshape(zk,ngrid(d),[]), [], ...
            'InitOnly', options);
        if k==2 % initialize prototype structure
            s = sk;
        end
        % rotate back, and assign
        data(:,k) = reshape(ipermute(reshape(sk.zzxx, ngrid(cs)), cs),...
            [],1);
    end
end

s = rmfield(s, {'dx', 'nx' 'zz', 'zzxx'});
% change or add new fields
s.Tag = 'splinend structure';
s.economy = 1;
s.nvar = n;
s.xgrid = x;
s.ngrid = ngrid;
s.tensordata=reshape(data,sizedata);

end

function val = getoption(options, field, default)
if isfield(options, field)
    val = options.(field);
else
    val = default;
end
end