DEIM.m

Go to the documentation of this file.

namespace general{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class DEIM
  :public KerMorObject,
   public general.AProjectable,
   public IReductionSummaryPlotProvider {
    public: /* ( setObservable ) */

        integer MaxOrder = 40;
    public: /* ( Dependent ) */

        rowvec<integer> Order;
    public:

        rowvec<double> SingularValues = "[]";
    protected:

        matrix<double> u;
        rowvec<double> pts;
    private:

        rowvec<integer> fOrder = "[]";
    public:

        matrix<double> U;
        matrix<double> U_nonproj;
        matrix<double> Uerr1;
        Uerr2;


        M1;

        M2;

        
    public:

        dscomponents.ACompEvalCoreFun f;
        rowvec<double> Residuals;
    public:

        function  computeDEIM(dscomponents.ACompEvalCoreFun f,matrix<double> fxi) {

            this.f= f;
            if ~isa(this.f," dscomponents.ACompEvalCoreFun ");
                error(" Cannot use DEIM with non ACompEvalCoreFun-implementing functions. ");
            elseif ~f.test_ComponentEvalMatch
                error(" Component evaluation does not match direct evaluation. ");
            elseif ~f.test_ComponentEvalMatch(10)
                error(" Component evaluation does not match direct evaluation for multi-arguments. ");
            end

            if size(fxi,1) < this.MaxOrder
                fprintf(2," Warning: MaxOrder %d of DEIM approximation is larger than actual f dimension %d. Setting MaxOrder=%d\n ",...
                    this.MaxOrder,size(fxi,1),size(fxi,1));
                this.MaxOrder= size(fxi,1);
            end

            /* % Generate u_1 ... u_m base */
            p = general.POD;
            p.Mode= " abs ";
            p.Value= this.MaxOrder;

            /* p.UseSVDS = size(fxi,1) > 10000; */
            [Utmp, this.SingularValues] = p.computePOD(fxi);
            if isa(Utmp," data.FileMatrix ")
                Utmp = Utmp.toMemoryMatrix;
            end
            this.u= Utmp;

            if size(this.u,2) < this.MaxOrder
                fprintf(" POD returned less (=%d) than MaxOrder (=%d) basis vectors. Setting MaxOrder=%d.\n ",...
                    size(this.u,2),this.MaxOrder,size(this.u,2));
                this.MaxOrder= size(this.u,2);
            end

            /*  Get interpolation points */
            this.pts= this.getInterpolationPoints(this.u);

            /*  Trigger computation of U matrix.
             * If no order has been set yet ("first computation"), use
             * predefined order. Otherwise, just ensure the order data
             * corresponds to the current u/pts. */
            if isempty(this.fOrder)
                this.Order= ceil(this.MaxOrder/10);
            else
                this.updateOrderData;
            end
        }
        function fx = evaluate(colvec<double> x,double t) {
            fx = this.U * this.f.evaluateComponentSet(1, x, t);
        }


        function fx = evaluateMulti(colvec<double> x,double t,colvec<double> mu) {
            fx = this.U * this.f.evaluateComponentSetMulti(1, x, t, mu);
        }


        function J = getStateJacobian(colvec<double> x,double t) {
            hlp = this.U * this.f.evaluateComponentSetGradientsAt(1, x, t);
            if isempty(this.V)
                p = logical(this.f.JSparsityPattern);
                J = sparse(this.f.fDim,this.f.xDim);
                J(p) = hlp(p);
            else
                J = hlp;
            end
        }


        function pts = getInterpolationPoints(matrix<double> u) {

            n = size(u,1);
            m = size(u,2);
            pts = zeros(1, m);
            v = pts;
            [v(1), pts(1)] = max(abs(u(:,1)));
            P = zeros(n,1);
            P(pts(1)) = 1;
            for i=2:m
                c = (P" *u(:,1:(i-1))) \ (P "*u(:,i));
                [v(i), pts(i)] = max(abs(u(:,i) - u(:,1:(i-1))*c));
                P = sparse(pts(1:i),1:i,ones(i,1),n,i);
            end
            this.Residuals= v;
            if KerMor.App.Verbose > 2
                fprintf(" DEIM interpolation points [%s] with values [%s]\n ",...
                    Utils.implode(pts,"   "," %d "),Utils.implode(v,"   "," %2.2e "));
            end
        }
        function err = getEstimatedError(colvec<double> x,double t,colvec<double> mu) {
            if this.fOrder(2) == 0
                error(" No error estimation possible with zero ErrorOrder property ");
            end
            err = this.Uerr1 * this.f.evaluateComponentSetMulti(1, x, t, mu) ...
                - this.Uerr2 * this.f.evaluateComponentSetMulti(2, x, t, mu);
        }


        function DEIMtarget = project(matrix<double> V,matrix<double> W,DEIM target) {

            if nargin < 4
                target = this.clone;
            end
            target = project@general.AProjectable(this, V, W, target);
            /*  Important: Project the component evaluable function, too! */
            target.f= this.f.project(V, W);
            target.updateOrderData;
        }
        function copy = clone(copy) {
            if nargin < 2
                copy = general.DEIM;
            end
            copy = clone@general.AProjectable(this, copy);
            /*  Clone associated f as different orders for the cloned object
             * lead to different PointSets for the component evaluable
             * function. */
            copy.f= this.f.clone;
            copy.fOrder= this.fOrder;
            copy.u= this.u;
            copy.pts= this.pts;
            copy.U= this.U;
            copy.U_nonproj= this.U_nonproj;
            copy.MaxOrder= this.MaxOrder;
            copy.Uerr1= this.Uerr1;
            copy.Uerr2= this.Uerr2;
            copy.M1= this.M1;
            copy.M2= this.M2;
            copy.SingularValues= this.SingularValues;
        }


        function  plotSummary(pm,context) {
            if ~isempty(this.SingularValues)
                str = sprintf(" %s: DEIM POD singular value decay, MaxOrder: %d ",context,this.MaxOrder);
                h = pm.nextPlot(" deim_singvals ",...
                    str," POD size "," singular values ");
                semilogy(h,this.SingularValues," LineWidth ",2);
            end
            if ~isempty(this.Residuals)
                str = sprintf(" %s: maximum basis residual, MaxOrder: %d ",context,this.MaxOrder);
                h = pm.nextPlot(" max_residuals ",...
                    str," Point number "," residual ");
                semilogy(h,this.Residuals," LineWidth ",2);
            end
        }

    
    protected:

        function  updateOrderData() {

            if KerMor.App.Verbose > 3
                fprintf(" general.DEIM.updateOrderData: Updating order data of DEIM (%s, #%s) to [%d %d]\n ",class(this),this.ID,this.fOrder);
            end

            o = this.fOrder(1);
            n = size(this.u,1);

            P = sparse(this.pts(1:o),1:o,ones(o,1),n,o);
            this.U_nonproj= this.u(:,1:o) / (P^t*this.u(:,1:o));

            /*  Set primary point set in ACompEvalCoreFun */
            this.f.setPointSet(1, this.pts(1:o));

            om = this.fOrder(2);
            if om > 0
                /*  Use second point set in ACompEvalCoreFun */
                this.f.setPointSet(2, this.pts(o+1:o+om));

                Perr = sparse(this.pts(o+1:o+om),1:om,ones(om,1),n,om);

                Um = this.u(:,1:o);
                Umd = this.u(:,o+1:o+om);
                A = P^t*Um;
                B = P^t*Umd;
                C = Perr^t*Um;
                D = Perr^t*Umd;
                E = C / A;
                F = D - E*B;
                this.Uerr2= ((Um/A)*B - Umd) / F;
                this.Uerr1= this.Uerr2 * E;
            else
                this.Uerr1= [];
                this.Uerr2= [];
            end

            if ~isempty(this.W)
                this.U= this.W^t*this.U_nonproj;
                if om > 0
                    /*  Compute values for error estimator */
                    this.M1= this.Uerr1 + (Um - this.V*(this.W^t*Um))/A;
                    this.M2= this.Uerr2;

                    /*  Project DEIM error estimation values */
                    this.Uerr1= this.W^t*this.Uerr1;
                    this.Uerr2= this.W^t*this.Uerr2;
                end
            else
                this.U= this.U_nonproj;
            end

            /*  Fire order updated event */
            notify(this, " OrderUpdated ", event.EventData);
        }
        /* % Getter & Setter */
    public:

        
#if 0 //mtoc++: 'get.Order'
function o = Order() {
            o = this.fOrder;
        }

#endif


        
#if 0 //mtoc++: 'set.Order'
function  Order(value) {
            if isempty(this.u) || isempty(this.pts)
                error(" Cannot set DEIM order as approximateSystemFunction has not been called yet. ");
            end
            if numel(value) ~= 2
                if numel(value) ~= 1
                    error(" Order must be either a scalar or two element integer vector ");
                else
                    value(2) = 0;
                end
            end
            if any(value) > size(this.u,2) || any(value) < 1
                error(" Invalid Order/ErrOrder value. Allowed are integers in [1, %d] ",size(this.u,2));
            elseif sum(value) > this.MaxOrder
                warning(" KerMor:DEIM ",sprintf(...
                    " Order (%d) and ErrOrder (%d) values may not exceed MaxOrder (%d). Using [%d 0] ",...
                    value,this.MaxOrder,this.MaxOrder));/* #ok */

                value = [this.MaxOrder 0];
            end
            /*  Only re-compute error order matrices if an actual change
             * happenend */
            if ~isequal(this.fOrder, value)
                this.fOrder= value;
                this.updateOrderData;
            end
        }

#endif


        
#if 0 //mtoc++: 'set.MaxOrder'
function  MaxOrder(value) {
            if ~isposintscalar(value)
                error(" MaxOrder has to be a positive integer scalar. ");
            end
            this.MaxOrder= value;
        }

#endif

    
    public:

        EVENT OrderUpdated;
    protected: /* ( Static ) */

        static function obj = loadobj(obj,varargin) {
            obj = loadobj@general.AProjectable(obj, varargin[:]);
            obj = loadobj@KerMorObject(obj, varargin[:]);
        }

};
}