TPWLLocalLipEstimator.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 TPWLLocalLipEstimator
  :public error.BaseEstimator {
    public:


        UseTimeDiscreteC = true;
    private:

        KernelLipschitzFcn;

        A;

        Ab;

        b;

        AB;

        bB;

        B;

        AV;

        neg_e1 = false;

        Ainorms;

        xi;

        
    public:

        TPWLLocalLipEstimator() {
            this.ExtraODEDims= 2;
        }


        function copy = clone() {
            error(" not yet implemented ");
        }


        function p = prepareForReducedModel(models.ReducedModel rmodel) {

            /*  Call superclass method to perform standard estimator
             * computations */
            p = prepareForReducedModel@error.BaseEstimator(this, rmodel);

            p.KernelLipschitzFcn= rmodel.System.f.GaussWeight.getLipschitzFunction;

            /*  Get full model */
            fm = rmodel.FullModel;
            if ~isempty(fm.System.B)
                try
                    fB = fm.System.B.evaluate([],[]);
                catch ME/* #ok */

                    fB = fm.System.B.evaluate(0,rmodel.getRandomParam);
                    warning(" Some:Id "," Error estimator for current system will not work correctly! (B is not linear and mu-independent! ");
                end
            end

            p.xi= fm.Approx.xi;

            Ai = fm.Approx.Ai;
            bi = fm.Approx.bi;
            G = rmodel.GScaled;
            n = length(Ai);
            /*  Perform any offline computations/preparations
             * Only prepare matrices if projection is used */
            if ~isempty(rmodel.V) && ~isempty(rmodel.W)
                /*  Compute projected versions (tilde..) */
                for i=1:n
                    Ai[i] = (Ai[i] - rmodel.V*(rmodel.W^t*Ai[i]))*rmodel.V;
                end
                bi = bi - rmodel.V*(rmodel.W^t*bi);

                p.A= cell(n,n);
                p.Ab= cell(n,n);
                p.b= zeros(n,n);
                [x,y] = meshgrid(1:n);
                all = [reshape(x,1,[]); reshape(y,1,[])];
                for k = 1:n*n
                    i = all(1,k);
                    j = all(2,k);
                    p.A[i,j] = Ai[i]^t*(G*Ai[j]);
                    p.Ab[i,j] = Ai[i]^t*(G*bi(:,j));
                    /* p.b(i,j) = bi(:,i)'*G*bi(:,j); */
                end
                p.b= bi^t*(G*bi);

                if ~isempty(fB)
                    p.AB= cell(1,n);
                    for i=1:n
                        p.AB[i] = Ai[i]^t*(G*fB);
                    end
                    p.bB= bi^t*(G*fB);
                    p.B= fB^t*(G*fB);
                end

                /*  Compute AV_i`s for \beta(s) comp */
                p.AV= cell(1,n);
                p.Ainorms= zeros(1,n);
                for i=1:n
                    p.AV[i] = Ai[i]^t*(G*Ai[i]);
                    p.Ainorms(i) = norm(Ai[i]);
                end
            else
                error(" \todo ");
/*                  % No projection means no projection error!
 *                 n = size(Ma,2);
 *                 
 *                 if ~isempty(B)
 *                     b = size(B,2);
 *                     p.M2 = zeros(n,b);
 *                     p.M3 = zeros(b,b);
 *                 end */
            end
        }
        function e = evalODEPart(colvec x,double t,double ut) {

            eold = x(end-this.ExtraODEDims+1:end);
            e = zeros(this.ExtraODEDims,1);
            wf = this.ReducedModel.System.f.GaussWeight;

            /*  Get actual z(s) */
            z = x(1:end-this.ExtraODEDims);

            /*  Get distances */
            di = this.xi - repmat(this.ReducedModel.V*z,1,size(this.xi,2));
            di = sqrt(sum(di.^2,1));
            di(di == 0) = eps;
            w = wf.evaluateScalar(di/min(di));
            w = w / sum(w);
            idx = 1:length(w);
            /*  Use same criterion as also used in normal approx evaluation */
            idx = idx(w > this.ReducedModel.System.f.MinWeightValue);
            m = length(idx);
            [x,y] = meshgrid(idx);
            all = [reshape(x,1,[]); reshape(y,1,[])];
            e(1) = 0;
            /*  Run semi-quadratic loop as many weights hopefully are zero */
            for k = 1:m*m
                i = all(1,k);
                j = all(2,k);
                e(1) = e(1) + w(i)*w(j)*(z" *this.A{i,j}*z + 2*z "*this.Ab[i,j] + this.b(i,j));
            end

            if nargin == 5 /*  An input function u is set */

                for i=1:m
                    e(1) = e(1) + w(i)*(z^t*this.AB[i]*ut + this.bB(i,:)*ut);
                end
                e(1) = e(1) + ut^t*this.B*ut;
            end

            if ~this.neg_e1 && e(1) < 0
                this.neg_e1= true;
            end
            /* e(1) = sqrt(abs(e(1))); */
            e(1) = sqrt(max(e(1),0));

            /* % Normal computations
             * Standard (worst-) Case */
            Ct = Inf;
            /*  Time-discrete computation */
            if this.UseTimeDiscreteC
                Ct = eold(1) + eold(2);
            end

            /*  Part one of beta */
            beta = wf.evaluateScalar(max(0,di-Ct)) * this.Ainorms^t;
            /*  Part two of beta */
            av = cellfun(@(A)sqrt(z^t*A*z),this.AV);
            beta = beta + av * this.KernelLipschitzFcn(di,Ct,t,this.mu)^t;

            e(2) = beta*(eold(1) + eold(2));

/*              % Iteration stuff
 *             if this.Iterations > 0
 *                 this.times(end+1) = t;
 *                 this.e1vals(end+1) = e(1);
 *                 this.divals(end+1,:) = di;
 *             end */
        }
        function ct = postProcess(double t,colvec<double> x,integer inputidx) {
            ct = postProcess@error.BaseEstimator(this, t, x, mu, inputidx);
            ts = tic;
            eint = x(end-this.ExtraODEDims+1:end,:);
            if all(eint == 0)
                warning(" CompWiseErrorEstimator:process "," Integral part is all zero. Attention! ");
            end
            this.StateError= eint(1,:) + eint(2,:);

            if this.neg_e1
                disp(" TPWLLocalLipEstimator: Negative alpha(t) norms occurred. Used zero instead. ");
                this.neg_e1= false;
            end

/*              % Iteration stuff
 *             if this.Iterations > 0
 *                 this.times(end+1) = t(end);
 *                 if this.UseTimeDiscreteC
 *                     warning('error:TPWLLocalLipEstimator','Using Iterations together with TimeDiscreteC will yield no improvement. Not performing iterations.');
 *                 else
 *                     this.performIterations(t, mu);
 *                 end
 *             end */

            /*  Tranform to output error estimation */
            C = this.ReducedModel.FullModel.System.C.evaluate(0,[]);
            this.StateError= norm(C)*this.StateError;

            ct = ct + toc(ts);
        }


        function e0 = init(colvec<double> mu) {
            e0 = [this.ReducedModel.getExo(mu); 0];
        }
        function  clear() {
            clear@error.BaseEstimator(this);
            this.neg_e1= false;
        }

    
   
    
    public: /* ( Static ) */

        static function errmsg = validModelForEstimator(models.ReducedModel rmodel) {
            errmsg = ;
            if ~isa(rmodel.FullModel.System.C," dscomponents.LinearOutputConv ")
                errmsg = " Local Lipschitz estimators work only for constant linear output conversion. ";
            end
            if isempty(errmsg) && ~isa(rmodel.FullModel.Approx," approx.TPWLApprox ")
                errmsg = " The Approx class of the full model must be a subclass of approx.TPWLApprox for this error estimator. "; 
            end
        }
};
}