intgrad1.m

The inverse of the gradient function. I ve provided versions that work on 1-d vectors, or 2-d or 3-d

function fhat = intgrad1(fx,dx,f1,method)
% intgrad: generates a vector, integrating derivative information.
% usage: fhat = intgrad1(dfdx)
% usage: fhat = intgrad1(dfdx,dx)
% usage: fhat = intgrad1(dfdx,dx,f1)
% usage: fhat = intgrad1(dfdx,dx,f1,method)
%
% arguments: (input)
%  dfdx - vector of length nx, as gradient would have produced.
%
%    dx - (OPTIONAL) scalar or vector - denotes the spacing in x
%         if dx is a scalar, then spacing in x (the column index
%         of fx and fy) will be assumed to be constant = dx.
%         if dx is a vector, it denotes the actual coordinates
%         of the points in x (i.e., the column dimension of fx
%         and fy.) length(dx) == nx
%
%         DEFAULT: dx = 1
%
%    f1 - (OPTIONAL) scalar - defines the first eleemnt of fhat
%         after integration. This is just the constant of integration.
%
%         DEFAULT: f1 = 0
%
%  method - (OPTIONAL) scalar - either 0, 1, 2, or 3. Defines
%         the integration scheme used.
%
%         method = 0 --> cumtrapz
%
%         method = 1 --> solves central finite difference
%                        approximation using linear algebra
%                        A second order fda. At least 3 points
%                        are necessary.
%
%         method = 2 --> integrated spline model
%                        This will almost always be the most
%                        accurate among the alternative methods.
%
%         method = 3 --> integrated pchip model
%
%         method = 4 --> higher order finite difference approximation
%                        A 4th order fda. At least 5 points are
%                        necessary.
%
%         DEFAULT: method = 2
%
%         Note: method = 0 (cumtrapz) will generally be the fastest,
%         and method = 2 (spline integral) will be the most accurate
%         of the four methods.
%         Methods 1, 3, and 4 were put in there mainly for fun on my
%         part, lthough for equally spaced points, the 4th order fda
%         should also be quite accurate.
%
%         Data series with noise in them may be best integrated using
%         a lower order method to avoid noise amplification.
%
% arguments: (output)
%   fhat - vector of length nx, containing the integrated function
%
% Example usage: 
%  x = 0:.001:1;
%  f = exp(x) + exp(-x);
%  dfdx = exp(x) - exp(-x);
%  tic,fhat = intgrad1(dfdx,.001,2,2);toc

% size 
if (length(fx)~=numel(fx))
  error 'dfdx must be a vector.'
end
sx = size(fx);
fx = fx(:);
nx = length(fx);
if nx<2
  error 'dfdx must be a vector of length >= 2'
end

% supply defaults if needed
if (nargin<2) || isempty(dx)
  % default x spacing is 1
  dx = 1;
end
if (nargin<3) || isempty(f1)
  % default integration constant is 0
  f1 = 0;
end
if (nargin<4) || isempty(method)
  % default integration method is 2
  method = 2;
end

% if scalar spacings, expand them to be vectors
dx=dx(:);
uniflag = 1;
if length(dx) == 1
  mdx = dx; % mean of dx
  xp = (0:(nx-1))'*dx;
  dx = repmat(dx,nx-1,1);
elseif length(dx)==nx
  % dx was a vector, use diff to get the spacing
  xp = dx;
  dx = diff(dx);
  mdx = mean(dx);
  ddx = diff(dx);
  if any(dx<=0)
    error 'x points must be monotonic increasing'
  elseif any(abs(ddx)>(xp(end)*1.e-15))
    uniflag = 0;
  end
else
  error 'dx is not a scalar or of length == nx'
end

if (length(f1) > 1) || ~isnumeric(f1) || isnan(f1) || ~isfinite(f1)
  error 'f1 must be a finite scalar numeric variable.'
end

if ~ismember(method,[0 1 2 3 4])
  error 'Method must be one of: [0, 1, 2, 3 4]'
end

switch method
  case 0
    % cumtrapz
    fhat = f1 + cumtrapz(xp,fx);
  case 1
    % builds a system using finite differences, then solves using \
    if nx<3
      error 'Method == 2 requires at least 3 points'
    end
    
    % build gradient design matrix, sparsely. Use a central difference
    % in the body of the array, and forward/backward differences along
    % the edges.
    
    % A will be the final design matrix. it will be sparse.
    Af = zeros(nx,9);
    
    % non-central difference at first point
    d = igfda([dx(1),dx(1)+dx(2)],1);
    Af(1,:) = [1 1 1 1 2 3,d];
    
    % interior df/dx, central difference
    i = (2:(nx-1))';
    if uniflag
      d = repmat([0 -1 1]./(2*mdx),nx-2,1);
    else
      dxi = [-dx(i-1),dx(i)];
      d = igfda(dxi,1);
    end
    Af(i,:) = [i, i, i, i, i-1, i+1,d];
    
    % non-central difference at last point
    d = igfda([-dx(end-1)-dx(end),-dx(end)],1);
    Af(nx,:) = [nx, nx, nx, nx, nx-2, nx-1, d];
    
    % finally, we can build the rest of A itself, in its sparse form.
    A = sparse(Af(:,1:3),Af(:,4:6),Af(:,7:9),nx,nx);
    
    % Finish up with f11, the constant of integration.
    % eliminate the first unknown, as f11 is given.
    fhat = A(:,2:end)\(fx - A(:,1)*f1);
    fhat = [f1;fhat];
    
  case 2
    % integrate a spline model
    pp = spline(xp,fx);
    c = pp.coefs;
    fhat = dx.*(c(:,4) + dx.*(c(:,3)./2 + dx.*(c(:,2)./3 + dx.*c(:,1)./4)));
    fhat = f1+[0;cumsum(fhat)];
    
  case 3
    % integrate a pchip model
    pp = pchip(xp,fx);
    c = pp.coefs;
    fhat = dx.*(c(:,4) + dx.*(c(:,3)./2 + dx.*(c(:,2)./3 + dx.*c(:,1)./4)));
    fhat = f1+[0;cumsum(fhat)];
    
  case 4
    % builds a system using finite differences, then solves using \
    if nx<5
      error 'Method == 4 requires at least 5 points.'
    end
    
    % build gradient design matrix, sparsely. Use a central difference
    % in the body of the array, and forward/backward differences along
    % the edges.
    
    % A will be the final design matrix. it will be sparse.
    Af = zeros(nx,15);
    
    % non-central difference at start
    d = igfda(cumsum(dx(1:4)'),1);
    Af(1,:) = [1 1 1 1 1, 1 2 3 4 5, d];
    d = igfda([-dx(1),cumsum(dx(2:4)')],1);
    Af(2,:) = [2 2 2 2 2, 2 1 3 4 5, d];
    
    % interior df/dx, central difference
    i = (3:(nx-2))';
    if uniflag
      mdx = mean(dx);
      %       d = fda(mdx*[-2 -1 1 2],1);
      d = [0 1/12 -2/3 2/3 -1/12]./mdx;
      d = repmat(d,nx-4,1);
    else
      dxi = [-dx(i-2)-dx(i-1), -dx(i-1), dx(i), dx(i)+dx(i+1)];
      d = igfda(dxi,1);
    end
    Af(i,:) = [i i i i i, i i-2 i-1 i+1 i+2, d];
    
    % non-central difference at end point
    d = igfda([-dx(end-3)-dx(end-2)-dx(end-1), ...
             -dx(end-2)-dx(end-1),-dx(end-1),dx(end)],1);
    Af(nx-1,:) = [repmat(nx-1,1,5), nx+[-1 -4 -3 -2 0], d];
    
    d = igfda([-dx(end-3)-dx(end-2)-dx(end-1)-dx(end), ...
             -dx(end-2)-dx(end-1)-dx(end), ...
             -dx(end-1)-dx(end),-dx(end)],1);
    Af(nx,:) = [repmat(nx,1,5), nx+[0 -4 -3 -2 -1], d];
    
    % finally, we can build the rest of A itself, in its sparse form.
    A = sparse(Af(:,1:5),Af(:,6:10),Af(:,11:15),nx,nx);
    
    % Finish up with f11, the constant of integration.
    % eliminate the first unknown, as f11 is given.
    fhat = A(:,2:end)\(fx - A(:,1)*f1);
    fhat = [f1;fhat];
    
end

% restore row or column shape of fhat
fhat = reshape(fhat,sx);

%=============================================
%   begin subfunctions
%=============================================
function coef = igfda(dxlist,estimator)
% fda: computes finite difference coefficients
% usage: coef = igfda(dxlist,estimator)
%
% Will compute one fda for each row of dxlist
%
% arguments: (input)
%  dxlist - (nxp) numeric arrray - list of deltas from the
%         central node of the fda.
%
%  estimator - (OPTIONAL) scalar - denotes the derivative order
%         to be estimated
%
%         estimator = k --> estimate the k'th derivative
%
%         DEFAULT: estimator = 1

if (nargin<2) || isempty(estimator)
  estimator = 1;
end

[n,p] = size(dxlist);
nfda = p + 1;

if (estimator < 0) || (estimator > (nfda-1))
  error 'Invalid derivative order requested'
end

% tolerance for consolidator to cluster the deltas
tol = max(abs(dxlist(:))*1e-14);
if size(n,1)>1
  [dxu,ind] = consolidator(dxlist,[],[],tol);
else
  dxu = dxlist;
  ind = 1:n;
end
nu = size(dxu,1);

% loop over the distinct rows of dxlist
coef = zeros(nu,nfda);
E = zeros(nfda,1);
E(nfda - estimator) = 1;
pwr = repmat(p:-1:0,nfda,1);
f = repmat(1./factorial((nfda-1):-1:0),nfda,1);
for i = 1:nu
  M = f.*(repmat([0;dxu(i,:)'],1,nfda).^pwr);
  coef(i,:) = (M'\E)';
end

% expand any consolidated sets
coef = coef(ind,:);