IterationCompLemmaEstimator.m

Go to the documentation of this file.

namespace error{


/* (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 IterationCompLemmaEstimator
  :public error.BaseCompLemmaEstimator {
    public:

        error.lipfun.Base LocalLipschitzFcn;
        Iterations;
        logical UseTimeDiscreteC = true;
    public: /* ( Transient ) */

        d_iValues;
    private:

        Ma_norms;

                c;
fIterations = 0;

        fTDC = true;

        G;

        smallz_flag = false;

        
    private: /* ( Transient ) */

        errEst;
        tstep;

        lstPreSolve;

        lstPostSolve;

        
    public: /* ( Transient ) */

        function copy = clone() {
            copy = error.IterationCompLemmaEstimator;
            /*  ExtraODEDims is set in constructor! */
            copy = clone@error.BaseCompLemmaEstimator(this, copy);
            /*  Clone local lipschitz function */
            copy.LocalLipschitzFcn= this.LocalLipschitzFcn.clone;
            copy.Iterations= this.Iterations;
            copy.UseTimeDiscreteC= this.UseTimeDiscreteC;
            copy.Ma_norms= this.Ma_norms;
            copy.c= this.c;
            copy.d_iValues= this.d_iValues;
            copy.errEst= this.errEst;
            copy.tstep= this.tstep;
            if ~isempty(copy.ReducedModel)
                copy.lstPreSolve= addlistener(copy.ReducedModel.ODESolver," PreSolve ",@copy.cbPreSolve);
                copy.lstPreSolve.Enabled= this.lstPreSolve.Enabled;
                copy.lstPostSolve= addlistener(copy.ReducedModel.ODESolver," PostSolve ",@copy.cbPostSolve);
                copy.lstPostSolve.Enabled= this.lstPostSolve.Enabled;
            end
            copy.G= this.G;
            copy.smallz_flag= this.smallz_flag;
        }
        function  offlineComputations(models.BaseFullModel fm) {

            offlineComputations@error.BaseCompLemmaEstimator(this, fm);

            /*  Obtain the correct snapshots */
            f = [];
            if ~isempty(fm.Approx)
                f = fm.Approx.Expansion;
            end
            if isempty(f)
                /*  This is the also possible case that the full core
                 * function of the system is a KernelExpansion.
                 *
                 * Get centers of full core function */
                f = fm.System.f.Expansion;
            end
            this.c= f.Centers;
            this.Ma_norms= f.Ma_norms;

            /*  Assign Lipschitz function */
            lfcn = error.lipfun.ImprovedLocalSecantLipschitz(f.Kernel);
            /*  Pre-Compute bell function xfeats if applicable */
            [x,X] = Utils.getBoundingBox(this.c.xi);
            d = norm(X-x);
            lfcn.precompMaxSecants(d*2,200);
            this.LocalLipschitzFcn= lfcn;
        }
        function prepared = prepareForReducedModel(models.ReducedModel rm) {
            prepared = prepareForReducedModel@error.BaseCompLemmaEstimator(this, rm);

            prepared.lstPreSolve= addlistener(rm.ODESolver," PreSolve ",@prepared.cbPreSolve);
            prepared.lstPostSolve= addlistener(rm.ODESolver," PostSolve ",@prepared.cbPostSolve);
            prepared.G= rm.G;

            f = rm.FullModel.Approx;
            if isempty(f)
                f = rm.FullModel.System.f;
            end
            if f.Expansion.Kernel.IsRBF && ~isempty(rm.V)
                /*  Use the projected centers z_i from the reduces system with x_i = Vz_i in this
                 * case. */
                hlp = f.project(rm.V, rm.W);
                prepared.c= hlp.Expansion.Centers;
                clear hlp;
                /*  The norm is then ||Vz - x_i||_G = ||z-z_i||_V'GV */
                prepared.G= rm.V^t*(prepared.G*rm.V);
                /*  Set flag for small z state vectors to true */
                prepared.smallz_flag= true;
            end
        }
        function b = getBeta(colvec xfull,t) {

            /*  Project x variable back to full space if centers are not within the space spanned by V
             * (indicated by smallz_flag) */
            x = xfull(1:end-1);
            if ~this.smallz_flag && ~isempty(this.ReducedModel.V)
                x = this.ReducedModel.V*x;
            end
            di = this.c.xi - repmat(x,1,size(this.c.xi,2));
            di = sqrt(sum(di.*(this.G*di),1));

            /* % Normal computations
             * Standard (worst-) Case */
            Ct = Inf;
            /*  Time-discrete computation */
            if this.UseTimeDiscreteC
                Ct = xfull(end);
                /*  Keep track of distances when iterations are used */
            elseif this.Iterations > 0
                this.d_iValues(this.StepNr,:) = di;
            end

            /*  Check if time and param kernels are used */
            if (~isempty(this.mu))
                f = this.ReducedModel.System.f;
                hlp = this.Ma_norms ...
                    .* f.Expansion.TimeKernel.evaluate(t, this.c.ti) ...
                    .* f.Expansion.ParamKernel.evaluate(this.mu, this.c.mui);
            else
                hlp = this.Ma_norms;
            end

            b = hlp * this.LocalLipschitzFcn.evaluate(di, Ct)^t;
        }
        function  clear() {
            clear@error.BaseCompLemmaEstimator(this);
            this.d_iValues= [];
        }


        function ct = prepareConstants(colvec<double> mu,integer inputidx) {

            if this.Iterations > 0 && this.UseTimeDiscreteC
                warning(" Ambiguous configuration. Having Iterations and UseTimeDiscreteC set; preferring UseTimeDiscreteC ");
            end

            /*  Call superclass method */
            ct = prepareConstants@error.BaseCompLemmaEstimator(this, mu, inputidx);
            st = tic;
            /*  Returns the initial error at `t=0` of the integral part. */
            if isempty(this.lstPreSolve)
                this.lstPreSolve= addlistener(this.ReducedModel.ODESolver," PreSolve ",@this.cbPreSolve);
            end
            this.lstPreSolve.Enabled= true;
            if isempty(this.lstPostSolve)
                this.lstPostSolve= addlistener(this.ReducedModel.ODESolver," PostSolve ",@this.cbPostSolve);
            end
            this.lstPostSolve.Enabled= true;
            /*  Call ISimConstants update function to compute values that are constant during a
             * simulation. */
            this.LocalLipschitzFcn.prepareConstants;
            ct = ct + toc(st);
        }
        function ct = postProcess(colvec<double> x,double t,integer inputidx) {
            st = tic;
            this.StateError= x(end,:);

            /*  PreSolve not needed during iterations */
            this.lstPreSolve.Enabled= false;

            /*  Iteration stuff */
            if ~this.UseTimeDiscreteC && this.Iterations > 0
                /*  Switch on/off listeners */

                solver = this.ReducedModel.ODESolver;
                e0 = this.getE0(this.mu);
                for it = 1:this.Iterations
                    /*  Set time-step counter to one */
                    this.tstep= 1;
                    /*  Solve */
                    [~, e] = solver.solve(@this.iterationODEPart, t, e0);
                end
                this.StateError= e;
            end

            /*  PostSolve still needed during iterations */
            this.lstPostSolve.Enabled= false;

            ct = postProcess@error.BaseCompLemmaEstimator(this, x, t, inputidx);

            ct = ct + toc(st);
        }
        /* % Getter & Setter */
    public: /* ( Transient ) */

        
#if 0 //mtoc++: 'set.LocalLipschitzFcn'
function  LocalLipschitzFcn(value) {
            if ~isa(value," error.lipfun.Base ")
                error(" LocalLipschitzFcn must be a error.lipfun.Base subclass. ");
            end
            this.LocalLipschitzFcn= value;
        }

#endif


        
#if 0 //mtoc++: 'set.Iterations'
function  Iterations(value) {
            if isempty(value)
                value = 0;
            elseif ~isscalar(value) || value < 0
                error(" Iterations value must be a non-negative integer. ");
            end
            this.Iterations= value;
        }

#endif


        
#if 0 //mtoc++: 'set.UseTimeDiscreteC'
function  UseTimeDiscreteC(value) {
            if isempty(value)
                value = false;
            elseif ~islogical(value)
                error(" The value must be a logical ");
            end
            this.UseTimeDiscreteC= value;
        }

#endif

    
    private: /* ( Transient ) */


        function e = iterationODEPart(double t,eold) {
            idx = this.tstep;

            if abs(this.EstimationData(1,idx)-t) > 100*eps
                error(" The ODE solver does not work as required by the iterative scheme. ");
            end

            if (~isempty(this.mu))
                f = this.ReducedModel.System.f;
                hlp = this.Ma_norms ...
                    .* f.Expansion.TimeKernel.evaluate(t, this.c.ti) ...
                    .* f.Expansion.ParamKernel.evaluate(this.mu, this.c.mui);
            else
                hlp = this.Ma_norms;
            end
            b = hlp * this.LocalLipschitzFcn.evaluate(...
                this.d_iValues(idx,:), this.errEst(idx))^t;

            e = b*eold + this.EstimationData(2,idx);

            this.EstimationData(3,this.tstep) = b;

            this.tstep= this.tstep+1;
        }
        function  cbPreSolve(sender,data) {
            this.d_iValues= zeros(length(data.Times),size(this.c.xi,2));
        }
        function  cbPostSolve(sender,data) {
            this.errEst= data.States(end,:);
        }
    public: /* ( Static ) */ /* ( Transient ) */

        static function errmsg = validModelForEstimator(models.BaseFullModel model) {
            errmsg = validModelForEstimator@error.BaseCompLemmaEstimator(model);
            if isempty(errmsg) && ~isempty(model.Approx) && ~isa(model.Approx," dscomponents.ParamTimeKernelCoreFun ")
                errmsg = " The model "" s approximation function must be a subclass of dscomponents.ParamTimeKernelCoreFun for this error estimator. ";
            end
            if isempty(errmsg) && isa(model.System.f," dscomponents.ParamTimeKernelCoreFun ") ...
                    && ~isa(model.System.f.Expansion.Kernel," kernels.BellFunction ")
                errmsg = " The system "" s kernel must be a kernels.BellFunction for this error estimator. ";
            end
        }
        static function res = test_IterationCompLemmaEstimator() {
            m = models.synth.KernelTest(10);
            e = error.IterationCompLemmaEstimator;
            e.UseTimeDiscreteC= false;
            e.Iterations= 4;

            m.ErrorEstimator= e;
            m.offlineGenerations;
            r = m.buildReducedModel;

            r.simulate(r.getRandomParam,[]);
            e.UseTimeDiscreteC= true;
            r.simulate(r.getRandomParam,[]);

            res = true;
        }

    
        /*      methods(Static,Access=protected)
     *         function s = loadobj(s)
     *             s = loadobj@error.BaseCompLemmaEstimator(s);
     *             this.lstPreSolve = addlistener(s.ReducedModel.ODESolver,'PreSolve',@s.cbPreSolve);
     *             this.lstPreSolve.Enabled = false;
     *             this.lstPostSolve = addlistener(s.ReducedModel.ODESolver,'PostSolve',@s.cbPostSolve);
     *             this.lstPostSolve.Enabled = false;
     *         end
     *     end */
};
}