function [K,Mx,xe] = covGrid(cov, xg, hyp, x, z, i)

% covGrid - Covariance function based on a grid.
%
% A grid covariance function k(x,z) is composed as a product
% k(x,z) = k1(x(i1),z(i1)) * .. * kp(x(ip),z(ip)) of p covariance functions
% operating on mutually disjoint components of the data x,z.
% The resulting covariance matrix is given by the Kronecker product
%   K = kron(kron(kron(...,K3),K2),K1) = K_p x .. x K_2 x K_1.
%
% The covariance function grid contains p grid factors xg = {x1,..,xp}
% so that for each of the factors xi a covariance function ki can be defined.
% A factor xi is of size (ni,di) so that the Kronecker grid g has size
% (n1,n2,..,np,D), where D=d1+d2+..+dp. Hence, the grid g contains N=n1*n2*..*np
% data points overall. Note that the factors xi do neither need to
% be sorted nor do they need to be 1d a priori. If a factor xi contains unevenly
% spaced values, we require di=1.
%
% A factor xi can contain:
%  a) a univariate axis xi of size (ni,1),
%  b) a multivariate axis xi of size (ni,di),
%     where the values need to be equispaced, or
%  c) a stationary equispaced subgrid composed of di univariate equispaced axes
%     {xi_1,xi_2,xi_di}, where each axis xij is of size (nij,1) so that
%     ni=ni_1*..*ni_di intended to be used with stationary covariance functions
%     When a stationary equispaced subgrid is specified, we silently assume
%     the covariance function to be stationary.
%
% For fast computations, we exploit two kinds of structure:
%  1) The Kronecker structure of the covariance matrix induced BY the p factors.
%     Hence, for p=1, there is nothing to gain here.
%  2) The Toeplitz or BTTB (Block-Toeplitz with Toeplitz Blocks) WITHIN a factor
%     if a grid factor xi is equispaced (or multivariate and equispaced).
%     Note that empirically, only for matrices of sizes above 500x500, the
%     FFT-based MVMs are faster than dense matrix operations.
%
% Some examples with sizes and domain:
%  - a single factor with a univariate axis, case a)
%     xg = { 5*rand(150,1) }                => N = 150, D = 1, dom = [0,5]
%     xg = { linspace(0,3,400)' }           => N = 400, D = 1, dom = [0,3]
%  - a single factor with a multivariate axis (equispaced is mandatory), case b)
%     xg = { [linspace(0,3,100)',...
%             linspace(1,4,100)'] }         => N = 100, D = 2, dom = [0,3]x[1,4]
%  - a single factor with a univariate equispaced&stationary subgrid, case c)
%     where we assume a stationary covariance and exploit Toeplitz structure
%     xg = { {linspace(0,10,175)'} }        => N = 175, D = 1, dom = [0,10]
%  - a single factor with a bivariate equispaced&stationary subgrid, case c)
%     where we assume a stationary covariance and exploit BTTB structure
%     xg = { {linspace(0,3,20)',...
%             linspace(1,7,30)'} }          => N = 600, D = 2, dom = [0,3]x[1,7]
% - a two-factor grid of 2 univariate axes, case a)
%     xg = { 2*rand(40,1), 5*rand(20,1) }   => N = 800, D = 2, dom = [0,2]x[0,5]
% - a two-factor grid of 2 univariate equispaced and stationary axes, case c)
%     where we assume two stationary covariances and exploit Toeplitz structure
%     in both factors
%     xg = { {linspace(0,2,25)'}, ...
%            {linspace(1,3,25)'} }          => N = 625, D = 2, dom = [0,2]x[1,3]
% - a four-factor grid with a Toeplitz factor, two ordinary Kronecker
%   factors and a 2d BTTB factor
%     xg = { {linspace(0,1,50)'}, ...       => N = 4e6, D = 5, dom = [0,1]^5
%            rand(20,1), ...
%            linspace(0,1,40)', ... 
%            {linspace(0,1,10)',linspace(0,1,10)'} }
%
% The covGrid function can be used to expand the (nested) cell array xg into a
% multivariate grid xe of size (N,D) via:
%     [xe,nx,Dx]  = covGrid('expand',xg);                             => mode 1)
% The operation can be reverted (if no subgrids are used) by:
%     xg = covGrid('factor',{xe,ng,Dg});                              => mode 2)
%
% Given scattered data x of size (n,D), we can create a grid xg covering
% the support of x using:
%     xg = covGrid('create',x,eq,k);                                  => mode 3)
% The flag eq (default value 1) can be used to enforce an equispaced
% grid. The integer k, indicates the number of grid points per dimension. If
% k is a real number from (0,1], then the number of grid points equals
% k*numel(unique(x(:,1))).
% We require at least two different components per dimension.
%
% The variables v={x,z} can either be a) grid indices or b) data points.
% a) The variable v has size (nv,1) and contains integers from [1,N]. Then
%    the datapoints are obtained as g2 = reshape(g,N,D); v = g2(v,:).
% b) The variable v has size (nv,D) and directly represents the data points.
% The mechanism works for x and z separately.
%
% An arbitrary data point x -- be it an index vector of size (n,1) or a data
% point of size (n,D) -- is converted into a regular data point xx of 
% size (n,D) by:
%     [xx,ng,Dg] = covGrid('idx2dat',xg,x);                           => mode 4)
% If x is already of size (n,D), xx will simply equal x.
%
% Given a grid xg and given arbitrary data points x, the interpolation
% matrix Mx can directly be computed without computing cross-covariances via:
%     [Mx,dMx] = covGrid('interp',xg,x);                              => mode 5)
% The cell array dMx contains the derivatives d Mx / d xi with respect to
% the i=1..p grid components.
%
% Given a nested grid xg, we can compute a flattened grid xf by:
%     xf = covGrid('flatten',xg);                                     => mode 6)
% without nesting and containing p axis elements.
%
% The hyperparameters are:
% hyp = [ hyp_1
%         hyp_2
%          ..
%         hyp_p ],
%
% Copyright (c) by Hannes Nickisch and Andrew Wilson 2015-09-26.
%
% See also COVFUNCTIONS.M, INFGRID.M.

if nargin<2, error('Not enough parameters provided.'), end
dense_max = 0; % 500 if larger grid we use FFT-algebra rather than dense algebra

% mode 1) expand axes xg representation into full grid x
if     strcmp(cov,'expand')          % call: [xe,nx,Dx]  = covGrid('expand',xg);
  [K,Mx,xe] = expandgrid(xg); return

% mode 2) factor full x grid into axes representation xg
elseif strcmp(cov,'factor')           % call: xg = covGrid('factor',{xe,ng,Dg});
  K = factorgrid(xg{:}); return

% mode 3) create axes representation xg from scattered data
elseif strcmp(cov,'create')               % call: xg = covGrid('create',x,eq,k);
  if nargin<3, eq = 1; else eq = hyp; end                   % set default values
  if nargin<4, k = 1; else k = x; end, x = xg;                % set input params
  K = creategrid(x,eq,k); return

% mode 4) convert possible index vector into data space
elseif strcmp(cov,'idx2dat')               % call: xx = covGrid('idx2dat',xg,x);
  [ng,Dg] = gridsize(xg); N = prod(ng); 
  if isidx(hyp,N), xe = expandgrid(xg); K = xe(hyp,:); else K = hyp; end
  Mx = ng; xe = Dg; return

% mode 5) compute interpolation matrix
elseif strcmp(cov,'interp')            % call [Mx,dMx] = covGrid('interp',xg,x);
  if nargout>1, [K,Mx] = interpgrid(xg,hyp); else K = interpgrid(xg,hyp); end
  return

% mode 6) provide flattened interpolation grid without nesting
elseif strcmp(cov,'flatten')                  % call xf = covGrid('flatten',xg);
  K = flattengrid(xg); return
end

% mode 0) regular covariance function computations
p = numel(xg); [ng,Dg] = gridsize(xg);             % number of Kronecker factors
if numel(cov)~=p, error('We require p factors.'), end
for ii = 1:p                                 % iterate over covariance functions
  f = cov(ii); if iscell(f{:}), f = f{:}; end   % expand cell array if necessary
  D = Dg(ii); j(ii) = cellstr(num2str(eval(feval(f{:}))));  % collect nbr hypers
end

if nargin<4                                        % report number of parameters
  K = char(j(1)); for ii=2:length(cov), K = [K, '+', char(j(ii))]; end, return
end
if nargin<5, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

v = [];               % v vector indicates to which covariance parameters belong
for ii = 1:p, v = [v repmat(ii, 1, eval(char(j(ii))))]; end
if nargin==6 && i>length(v), error('Unknown hyperparameter'), end

N = prod(ng); n = size(x,1); D = sum(Dg);     % expanded grid and data dimension
ix = isidx(x,N);               % determine whether x is an index or a data array
if ~ix && size(x,2)~=D, error('Grid and data dimension are different.'), end
if nargout>1 || ~(dg||xeqz||ix), Mx = interpgrid(xg,x); end  % off-grid interpol
if dg               % evaluate as full dense vector for diagonal covariance case
  K = 1;                       % xg is not assumed to form a grid for z = 'diag'
  for ii = 1:length(cov)                       % iteration over factor functions
    f = cov(ii); if iscell(f{:}), f = f{:}; end % expand cell array if necessary
    d = sum(Dg(1:ii-1))+(1:Dg(ii));                     % dimensions of interest
    if nargin<6, i = 0; vi = 0; else vi = v(i); end; % which covariance function
    if i<=length(v)
      if ix, xii = xg{ii}; else xii = x(:,d); end  % switch Kronecker/plain prod
      if ii==vi
        j = sum(v(1:i)==vi);                % which parameter in that covariance
        Kj = feval(f{:}, hyp(v==ii), xii, z, j);        % deriv Kronecker factor
      else
        Kj = feval(f{:}, hyp(v==ii), xii, z);           % plain Kronecker factor
      end
      if ix, K = kron(K,Kj); else K = K.*Kj; end   % switch Kronecker/plain prod
    else error('Unknown hyperparameter')
    end
  end
  if ix, K = K(x); end, return
end

if isidx(z,N), iz = z; z = covGrid('expand',xg); z = z(iz,:); end     % expand z
K = cell(p,1); sz = [1,1];             % covariance Kronecker factors total size
for ii = 1:p                                   % iteration over factor functions
  f = cov(ii); if iscell(f{:}), f = f{:}; end   % expand cell array if necessary
  d = sum(Dg(1:ii-1))+(1:Dg(ii));                       % dimensions of interest
  if isnumeric(z) && ~isempty(z)                                   % cross terms
    zd = z(:,d);
  elseif xeqz && iscell(xg{ii})                                  % Toeplitz/BTTB
    zd = [];
  else                                                        % symmetric matrix
    zd = z;
  end
  if nargin<6, i = 0; vi = 0; else vi = v(i); end;   % which covariance function
  if i<=length(v)
    if ii==vi
      pars = {zd,sum(v(1:i)==vi)};          % which parameter in that covariance
    else
      pars = {zd};
    end
    xii = xg{ii};                                                  % grid factor
    if xeqz && iscell(xg{ii})                                    % Toeplitz/BTTB
      di = numel(xii); ni = gridsize(xii);
      wi = zeros(di,1); for j=1:di, wi(j) = xii{j}(end)-xii{j}(1); end   % width
      eqstat = true; for j=1:di, eqstat = eqstat & equi(xii,j); end
      if ~eqstat, error('Subgrid not equispaced.'), end        % stop if problem
      if prod(ni)>dense_max  % empirical thresh: FFT-based MVM with 1 rhs faster
        kii = @(n) feval(f{:}, hyp(v==ii), (n-1)*diag(wi./(ni-1)), ...  % k(int)
                                                      zeros(1,di), pars{2:end});
        xc = cell(di,1);              % generic (Strang) circular embedding grid
        for j=1:di, n2 = floor(ni(j)-1/2)+1; xc{j} = [1:n2,n2-2*ni(j)+2:0]'; end
        ci = reshape(kii(covGrid('expand',xc)),[2*ni(:)'-1,1]);
        fi = real(fftn(ci));                % precompute FFT for circular filter
        mvmKi = @(x) bttbmvmsymfft(fi,x);        % MVM with Toeplitz/BTTB matrix
        if di==1, s = 'toep'; else s = ['bttb',num2str(di)]; end
        ki = struct('descr',s, 'mvm',mvmKi, 'kii',kii, 'size',prod(ni)*[1,1]);
      else             % simply evaluate covariance matrix if prod(ni) too small
        ki = feval(f{:},hyp(v==ii),expandgrid(xii),[],pars{2:end});
      end
      K{ii} = ki; sz = sz.*K{ii}.size;
    else
      if iscell(xii), xii = expandgrid(xii); end
      K{ii} = feval(f{:}, hyp(v==ii), xii, pars{:});    % plain Kronecker factor
      sz = sz.*size(K{ii});
    end
  else error('Unknown hyperparameter')
  end
end

if xeqz                                                    % create mvm and rest
  K = struct('mvm',@(a)kronmvm(K,a),'size',sz,'kronmvm',@kronmvm,...
                                                     'kron',struct('factor',K));
else                                                        % expand cross terms
  Ks = K; K = Ks{1}; for ii = 2:p, K = kron1(Ks{ii},K); end
  if ix, if numel(x)~=N || max(abs(x-(1:N)'))>0, K = K(x,:); end
  else   K = Mx*K;
  end
end
if nargout>2, xe = covGrid('expand',xg); end

% Perform a matrix vector multiplication b = A*x with a matrix A being a
% Kronecker product given by A = kron( kron(...,As{2}), As{1} ).
function b = kronmvm(As,x,transp)
if nargin>2 && ~isempty(transp) && transp   % transposition by transposing parts
  for i=1:numel(As)
    if isnumeric(As{i})
      As{i} = As{i}';
    else
      As{i}.mvm = As{i}.mvmt;
      As{i}.size = [As{i}.size(2),As{i}.size(1)];
    end
  end
end
m = zeros(numel(As),1); n = zeros(numel(As),1);                  % extract sizes
for i=1:numel(n)
  if isnumeric(As{i})
    [m(i),n(i)] = size(As{i});
  else
    m(i) = As{i}.size(1); n(i) = As{i}.size(2);
  end
end
d = size(x,2);
b = x;
for i=1:numel(n)
  a = reshape(b,[prod(m(1:i-1)), n(i), prod(n(i+1:end))*d]);    % prepare  input
  if isnumeric(As{i})                                % apply along one dimension
    tmp = reshape(permute(a,[1,3,2]),[],n(i))*As{i}';
    b = permute(reshape(tmp,[size(a,1),size(a,3),m(i)]),[1,3,2]);
  else
    b = zeros([prod(m(1:i-1)), m(i), prod(n(i+1:end))*d]);      % prepare output
    for j=1:size(b,3), b(:,:,j) = As{i}.mvm(a(:,:,j)')'; end
  end
end
b = reshape(b,prod(m),d);                        % bring result in correct shape

% Perform MVM b = T*a with a of size (n,m) with a BTTB (Block-Toeplitz with
% Toeplitz-blocks) matrix T of size (n,n) by pointwise multiplication with the
% Fourier-transformed filter f. All variables are assumed real valued.
% Needs O(3*m*n*log(n)) time and O(n*m) space.
function b = bttbmvmsymfft(f,a)
  m = size(a,2); ng = (size(f)+1)/2; p = numel(ng);              % extract sizes
  if p==2 && ng(2)==1, p = 1; ng = ng(1); end           % detect 1d and reduce p
  b = reshape(a,[ng,m]);
  for i=1:p, b = fft(b,2*ng(i)-1,i); end                                  % fftn
  if exist('bsxfun','builtin')            % use fast multiplication if available
    b = bsxfun(@times,f,b);
  else
    b = repmat(f,[ones(1,p),m]).*b;
  end
  for i=1:p, b = ifft(b,[],i); end                                        % ifft
  for i=1:p                          % only keep the relevant part of the result
    b = reshape(b,prod(ng(1:i-1)),2*ng(i)-1,prod(2*ng(i+1:p)-1)*m);
    b = b(:,1:ng(i),:);
  end
  b = real(reshape(b,[],m));

% perform kron along first dimension only
% the code is equivalent to the following loop
%   z = zeros(size(x,1)*size(y,1),size(x,2));
%   for i=1:size(z,2), z(:,i) = kron(x(:,i),y(:,i)); end
function z = kron1(x,y)
  nx = size(x,1); ny = size(y,1);
  z = repmat(reshape(x,1,nx,[]),[ny,1,1]).*repmat(reshape(y,ny,1,[]),[1,nx,1]);
  z = reshape(z,nx*ny,[]);

function r = isidx(i,N)     % check whether i represents an integer index vector
  r = false;
  if numel(i)>0 && ~strcmp(i,'diag') && size(i,2)==1 && ndims(i)==2
    if max(abs(i-floor(i)))<1e-13
      if 0<min(i) && max(i)<=N, r = true; end
    end
  end
  
function [ng,Dg] = gridsize(xg)          % report the size of the p grid factors
  p = numel(xg); ng = zeros(p,1); Dg = zeros(p,1);      % number of grid factors
  for i=1:p                                          % iterate over grid factors
    x = xg{i};
    if iscell(x)                                 % stationary and equispace grid
      for j=1:numel(x)
        if j==1
          [ng(i),Dg(i)] = size(x{j});
        else
          ng(i) = ng(i)*size(x{j},1);
          Dg(i) = Dg(i)+size(x{j},2);
        end
      end
    else                                                        % arbitrary grid
      [ng(i),Dg(i)] = size(x);
    end
  end

% mode 1 (expand)
function [x,ng,Dg] = expandgrid(xg)                    % expand a Kronecker grid
  [ng,Dg] = gridsize(xg);                                        % original size
  xg = flattengrid(xg);                                      % remove nestedness
  p = numel(xg); x = xg{1};                                 % expanded grid data
  ngf = zeros(p,1); Dgf = zeros(p,1); [ngf(1),Dgf(1)] = size(xg{1});
  for i=2:p
    szx = size(x); [ngf(i),Dgf(i)] = size(xg{i});
    xold = repmat(reshape(x,[],1,szx(end)),[1,ngf(i),1]);
    xnew = repmat(reshape(xg{i},[1,ngf(i),Dgf(i)]),[size(xold,1),1,1]);
    x = reshape(cat(3,xold,xnew),[szx(1:end-1),ngf(i),szx(end)+Dgf(i)]);
  end
  x = reshape(x,[],size(x,ndims(x)));

% mode 2 (factor)
function xg = factorgrid(x,ng,Dg)                      % factor a Kronecker grid
  p = numel(ng); xg = cell(p,1);         % extract individual grid components xg
  for i=1:p
    x = reshape(x,[prod(ng(1:i-1)), ng(i), prod(ng(i+1:end)), sum(Dg)]);
    xg{i} = reshape(x(1,:,1,sum(Dg(1:i-1))+(1:Dg(i))), ng(i), Dg(i));
  end

% mode 3 (create)
function xg = creategrid(x,eq,k)
  if nargin<2, eq = 1; end                                  % set default values
  if nargin<3, k = 1; end                                     % set input params
  p = size(x,2); xg = cell(p,1);                               % allocate result
  if numel(k)>0, k = ones(p,1).*k(:); end              % enforce vector-valued k
  for j=1:p                                            % iterate over dimensions
    u = sort(unique(x(:,j))); if numel(u)<2, error('Two few unique points.'),end
    if isempty(k)                              % determine number of grid points
      if eq
        ngj = ceil( (u(end)-u(1))/min(abs(diff(u))) );     % use minimum spacing
      else
        ngj = numel(u);
      end
    elseif 0<=k(j) && k(j)<=1
      ngj = ceil(k(j)*numel(u));
    else
      ngj = k(j);
    end
    du = (u(end)-u(1))/ngj; bu = [u(1)-5*du, u(end)+5*du];
    if eq                                                      % equispaced grid
      xg{j} = linspace(bu(1),bu(2),max(ngj,5))';        % at least 5 grid points
    else                                                   % non-equispaced grid
      [idx,xgj] = kmeans(u,min(numel(u),ngj-2)); xgj = sort(xgj(:))';  % cluster
      nb = ngj-numel(xgj); nb1 = floor(nb/2); nb2 = nb - nb1; % size of boundary
      xg1 = linspace(bu(1),xgj(1),nb1+1); xg2 = linspace(xgj(end),bu(2),nb2+1);
      xg{j} = [xg1(1:nb1), xgj, xg2(1+(1:nb2))]';
    end
  end

% mode 5 (interp)
function [Mx,dMx] = interpgrid(xg,x)
  xg = flattengrid(xg);                                      % remove nestedness
  p = numel(xg); Dg = zeros(p,1); ng = zeros(p,1); n = size(x,1);      % dims ..
  for i=1:p, [ng(i),Dg(i)] = size(xg{i}); end, N = prod(ng);       %.. and sizes
  ix = isidx(x,N);             % determine whether x is an index or a data array
  if ix
    Mx = sparse(1:n,x,1,n,N); if nargout>1, dMx = repmat({sparse(n,N)},p,1); end
  else
    deg = 3;                     % degree of equispaced interpolation polynomial
    s = 1;                                                      % initial stride
    for i=1:p                                 % iterate over Toeplitz components
      d = sum(Dg(1:i-1))+(1:Dg(i));                     % dimensions of interest
      xt = xg{i}; it = find(abs(xt(2,:)-xt(1,:)));        % grid nonzero inc idx
      if equi(xg,i)                         % compute interpolation coefficients
        if nargout>1                                       % equispaced grid pts
          [Ji,Ci,dCi] = eqinterp(xt(:,it),x(:,d(it)),deg);
        else
          [Ji,Ci] = eqinterp(xt(:,it),x(:,d(it)),deg);
        end
      else
        if nargout>1                                   % non-equispaced grid pts
          [Ji,Ci,dCi] = neqinterp(xt(:,it),x(:,d(it)));
        else
          [Ji,Ci] = neqinterp(xt(:,it),x(:,d(it)));
        end
      end
      nc = size(Ci,2);    % number of interpolation coefficients along dimension
      if i==1
        C = Ci; J = ones(n,1);
        if nargout>1, dC = repmat({Ci},p,1); dC{1} = dCi; end
      else
        C = repmat(C,[1,1,nc]) .* repmat(reshape(Ci,n,1,nc),[1,size(C,2),1]);
        C = reshape(C,n,[]);
        if nargout>1
          for j=1:p
            if i==j, dCij = dCi; else dCij = Ci; end
            dC{j} = repmat(dC{j},[1,1,nc]) .* ...
                    repmat(reshape(dCij,n,1,nc),[1,size(dC{j},2),1]);
            dC{j} = reshape(dC{j},n,[]);
          end
        end
      end
      J = repmat(J(:),[1,nc]) + s*repmat(Ji-1,[size(C,2)/nc,1]);  % blow 2nd idx
      s = s*ng(i);                                               % update stride
    end
    I = repmat((1:n)',[1,size(C,2)]); id = 0<J&J<=N;% first index and valid flag
    Mx = sparse(I(id),J(id),C(id),n,N);
    if nargout>1
      dMx = cell(p,1); for i=1:p,dMx{i} = sparse(I(id),J(id),dC{i}(id),n,N); end
    end
  end

% mode 6 (flatten)
function xf = flattengrid(xg)               % convert nested grid into flat grid
  xf = cell(1,0);
  for i=1:numel(xg)
    x = xg{i}; if iscell(x), xf = [xf,x]; else xf = [xf,{x}]; end
  end  

function eq = equi(xg,i)                        % grid along dim i is equispaced
  ni = size(xg{i},1);
  if ni>1                              % diagnose if data is linearly increasing
    dev = abs(diff(xg{i})-ones(ni-1,1)*(xg{i}(2,:)-xg{i}(1,:)));
    eq = max(dev(:))<1e-9;
  else
    eq = true;
  end

% Compute interpolation coefficients C (nt,nc) and interpolation coefficient
% indices J (nt,nc) from a source grid s (ns,1) to a target array t (nt,1).
% The coefficient matrix C has rows summing up to 1.
function [J,C,dC] = eqinterp(s,t,d)
  gp = false;
  switch d
    case 0, k = @(x) -0.5<x & x<=0.5; it=-1:0; dk = @(x) 0*x;
    case 1, k = @(x) max(1-abs(x),0); it=-1:0; dk = @(x) -(abs(x)<=1).*sign(x);
    case 3, k = @kcub;                it=-2:1; dk = @dkcub;
    otherwise, ell = d/5; gp = true;                          % GP interpolation
      k = @(x) exp(-x.*x/(2*ell^2));  it=(0:d-1)-floor(d/2);
      dk = @(x) -k(x).*x/(ell^2);
  end
  ds = s(2)-s(1); ns = numel(s); nt = numel(t); nc = numel(it);
  if size(s,2)*size(t,2)~=1, error('Interpolation only possible for d==1.'), end
  if ns<nc, error('Interpolation only possible for ns>=nc.'), end
  j = floor((t-s(1))/ds)+1;                % index of closest smaller grid point
  w = (t-s(1))/ds-j+1;   % relative distance to closest smaller grid point [0,1]
  j = j-it(nc); C = zeros(nt,nc); dC = zeros(nt,nc);
  for i=1:nc
    C(:,i) = k(w+it(nc+1-i)); if nargout>2, dC(:,i)=dk(w+it(nc+1-i))*(ns-1); end
  end
  if gp, kn = k(sqrt(sq_dist(1:d))); C = C/kn; dC = dC/kn; end
  v = 1; id = find(j<nc+it(1)); C(id,:) = 0; dC(id,:) = 0;  % fix lower boundary
  D = abs(repmat(s(1:nc)',numel(id),1)-repmat(t(id),[1,nc]));
  [junk,jid] = min(D,[],2);          % index of closest index in boundary region
  for i=1:numel(id), C(id(i),jid(i)) = 1; dC(id(i),jid(i)) = 0; end, j(id) = v;
  v = ns-nc+1; id = find(j>v); C(id,:) = 0; dC(id,:) = 0;   % fix upper boundary
  D = abs(repmat(s(ns-nc+1:ns)',numel(id),1)-repmat(t(id),[1,nc]));
  [junk,jid] = min(D,[],2);          % index of closest index in boundary region
  for i=1:numel(id), C(id(i),jid(i)) = 1; dC(id(i),jid(i)) = 0; end, j(id) = v;
  J = zeros(nt,nc); for i=1:nc, J(:,i) = j+i-1; end    % construct index array J

% Robert G. Keys, Cubic Convolution Interpolation for Digital Image Processing,
% IEEE ASSP, 29:6, December 1981, p. 1153-1160.
function y = kcub(x)
  y = zeros(size(x)); x = abs(x);
  q = x<=1;          % Coefficients:  1.5, -2.5,  0, 1
  y(q) =            (( 1.5 * x(q) - 2.5) .* x(q)    ) .* x(q) + 1;
  q = 1<x & x<=2;    % Coefficients: -0.5,  2.5, -4, 2
  y(q) =            ((-0.5 * x(q) + 2.5) .* x(q) - 4) .* x(q) + 2;
function y = dkcub(x)
  y = sign(x); x = abs(x);
  q = x<=1;          % Coefficients:  1.5, -2.5,  0, 1
  y(q) = y(q) .*  ( 4.5 * x(q) - 5.0) .* x(q);
  q = 1<x & x<=2;    % Coefficients: -0.5,  2.5, -4, 2
  y(q) = y(q) .* ((-1.5 * x(q) + 5.0) .* x(q) - 4.0);
  y(x>2) = 0;

% Perform piecewise linear interpolation using inverse distance weighting.
% s (ns,1) source nodes, need neither be sorted nor equispaced
% t (nt,1) target nodes
% M (nt,ns) interpolation matrix, M = sparse((1:N)'*[1,1],J,C,nt,ns);
%
% z = M*y where y (ns,1) are source values and z (nt,1) are target values
function [J,C,dC] = neqinterp(s,t)
  ns = size(s,1); nc = 2;                                       % get dimensions
  if size(s,2)*size(t,2)~=1, error('Interpolation only possible for d==1.'), end
  if ns<nc, error('Interpolation only possible for ns>=nc.'), end
  [s,ord] = sort(s); ds = diff(s); 
  if min(ds)<1e-10, error('Some source points are equal.'), end
  [junk,ii] = histc(t(:),[-inf;s(2:end-1);inf]);
  d0 = t(:)-s(ii); d1 = s(ii+1)-t(:); d0n = d0<0; d1n = d1<0;
  d0(d0n) = 0; d1(d1n) = 0;                                % boundary conditions
  J = [ord(ii),ord(ii+1)]; C = [d1./(d1+d0),d0./(d1+d0)];
  nz = 1-(d1n|d0n); if nargout>2, dC = [-nz./(d1+d0),nz./(d1+d0)]; end
