funcApprox class

Defines a class to represent an approximated function

Objects created by this class are a subclass of basis, adding fields to identify a function and methods to interpolate, compute Jacobian and Hessian.

Apart from the properties inherited from basis, objects of class funcApprox have the following properties:

df: scalar, number of interpolated functions
fnodes: value of interpolated function at basis nodes
coef: interpolation coefficients
Phi: interpolation matrix, evaluated at basis nodes
Phiinv: inverse of Phi
fnames: cell of scalars, optional names for interpolated functions

Object of class funcApprox have the following methods:

funcApprox: class constructor
updateCoef: computes interpolation coefficients if fnodes is modified
Interpolate: interpolates the function
Jacobian: computes the Jacobian of the function
Hessian: computes the Hessian of the function

To date, only basisChebyshev has been implemented, so calling the function with type 'spli' or 'lin' returns an error.

Last updated: October 4, 2014.

Licensed under the MIT license, see LICENSE.txt

funcApprox
set.fnodes
updateCoef
Interpolate
Jacobian
Hessian

classdef funcApprox < basis

    properties
        fnodes  % value of functions at nodes
        fnames  % names for functions
    end


    properties (SetAccess = protected)
        Phi     % interpolation matrix
        Phiinv  % inverse of interpolation matrix
        df      % number of functions
        coef    % interpolation coefficients
    end


    methods

funcApprox

        function F = funcApprox(B,fnodes,fnames)

Constructor for funcApprox

F = funcApprox(B,fnodes,fnames)

The inputs are:

B: basis object
fnodes: matrix, values of the interpolant at basis B nodes
fnames: cell of strings, names of functions (optional)

            if nargin==0
                return
            end

            for s = fieldnames(B)'
                s1 = char(s);
                F.(s1) = B.(s1);
            end

            F.Phi = B.Interpolation;
            F.Phiinv = (F.Phi'*F.Phi)\F.Phi';
            F.fnodes = fnodes;
            F.df = size(fnodes,2);

            if nargin>2
                F.fnames = fnames;
            end

        end

set.fnodes

        function set.fnodes(F,value)

Setting new values for fnodes also update the interpolation coefficients

            if size(value,1) ~= size(F.nodes,1)
                error('New value for fnodes must have %d rows (1 per node)',size(F.nodes,1))
            end

            F.fnodes = value;
            F.updateCoef;
        end

updateCoef

        function updateCoef(F)

To update coefficients, the inverse interpolation matrix is premultiplied by the value of the function at the nodes

            F.coef = F.Phiinv * F.fnodes;
        end

Interpolate

        function y = Interpolate(F,...
                varargin... inputs for interpolation method
                )

y = F.Interpolate(x,order,integrate)

Interpolates function f: R^d --> R^m at values x; optionally computes derivative/integral. It obtains the interpolation matrices by calling Interpolation on its basis.

Its inputs are

x, k.d matrix of evaluation points.
order, h.d matrix, for order of derivatives. Defaults to zeros(1,d).
integrate, logical, integrate if true, derivative if false. Defaults to false

Output y returns the interpolated functions as a k.m.h array, where k is the number of evaluation points, m the number of functions (number of columns in f.coef), and h is number of order derivatives.

            if nargin <2
                y = F.fnodes;
                return
            end

            Phix = F.Interpolation(varargin{:});

            nx = size(Phix,1);  % number of evaluation points
            no = size(Phix,3);  % number of order evaluations

            y = zeros(nx,F.df,no);

            for h = 1:no
                y(:,:,h) = Phix(:,:,h) * F.coef;
            end

            if F.df==1  % only one function
                y = squeeze(y);
            end

        end %Evaluate

Jacobian

        function [DY,Y] = Jacobian(F, x, index)

[DY,Y] = F.Jacobian(x, index)

Computes the Jacobian of the approximated function f: R^d --> R^m.

Inputs:

x, k.d matrix of evaluation points.
index, 1.d boolean, take partial derivative only wrt variables with index=true. Defaults to true(1,d).

Outputs

DY, k.m.d1 Jacobian matrix evaluated at x, where d1 = sum(index)
Y, k.m interpolated function (same as Y = B.Evaluate(coef,x), provided for speed if both Jacobian and funcion value are required).

Solve for the one-dimensional basis

            if F.d==1
                if nargout == 1
                    Phix  = F.Interpolation(x,1,false);
                    DY = Phix * F.coef;
                else
                    Phix = F.Interpolation(x,[1;0],false);
                    DY = Phix(:,:,1) * F.coef;
                    Y = Phix(:,:,2) * F.coef;
                end

                return
            end

Solve for the multi-dimensional basis

            % Check validity of input x
            if size(x,2)~=F.d, error('In Jacobian, class basis: x must have d columns'); end

Keep track of required derivatives: Required is logical with true indicating derivative is required

            if nargin<3
                Required = true(1,F.d);
                index = 1:F.d;
            elseif numel(index) < F.d % assume that index have scalars of the desired derivatives
                Required = false(1,F.d);
                Required(index) = true;
            else % assume that index is nonzero for desired derivatives
                Required = logical(index);
                index = find(index);
            end

            nRequired = sum(Required);

HANDLE NODES DIMENSION

            Nrows = size(x,1);

HANDLE POLYNOMIALS DIMENSION

            c = F.opts.validPhi;
            Ncols = size(c,1);

Compute interpolation matrices for each dimension

            Phi0 = zeros(Nrows,Ncols,F.d);
            Phi1 = zeros(Nrows,Ncols,F.d);

            for k = 1:F.d
                if Required(k)
                    PhiOneDim = F.B1(k).Interpolation(x(:,k),...
                        [0 1],...
                        false);
                    Phi01 = PhiOneDim(:,c(:,k),:);
                    Phi0(:,:,k) = Phi01(:,:,1); % function value
                    Phi1(:,:,k) = Phi01(:,:,2); % its derivative
                else
                    PhiOneDim = F.B1(k).Interpolation(x(:,k),...
                        0,...
                        false);
                    Phi0(:,:,k) = PhiOneDim(:,c(:,k));
                end
            end

Compute the Jacobian Preallocate memory

            DY  = zeros(Nrows,F.df,nRequired);

            % Multiply the 1-dimensional bases

            for k=1:nRequired
                Phik = Phi0;
                Phik(:,:,index(k)) = Phi1(:,:,index(k));
                Phix = prod(Phik,3);
                DY(:,:,k) = Phix * F.coef;
            end

Compute the function if requested

            if nargout > 1
                Y = prod(Phi0,3) * F.coef;
            end

        end % Jacobian

Hessian

        function Hy = Hessian(F,x)

Hy = F.Hessian(x,coef)

Computes the Hessian of a function approximated by basis B and coefficients coef.

Its input is

x, k.d matrix of evaluation points.

Its output Hy returns the k.m.d.d Hessian evaluated at x.

            order = repmat({[0 1 2]'},1,F.d);
            order = gridmake(order{:});
            order = order(sum(order,2)==2,:);

            Phix = F.Interpolation(x,order,false);


            nx = size(x,1);     % number of evaluated points

            %Dy = squeeze(Dy);

            Hy = zeros(nx,F.df,F.d,F.d);

            for k = 1:size(order,1)
                i = find(order(k,:));
                if numel(i)==1
                    Hy(:,:,i,i) = Phix(:,:,k) * F.coef;%  Dy(:,k);
                else
                    Hy(:,:,i(1),i(2)) = Phix(:,:,k) * F.coef; %Dy(:,k);
                    Hy(:,:,i(2),i(1)) = Hy(:,:,i(1),i(2)); %Dy(:,k);
                end
            end
        end % Hessian

    end

end