EstimatorAnalyzer.m

Go to the documentation of this file.

/* (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 EstimatorAnalyzer
  :public handle {
    public:

        models.BaseFullModel Model;
        models.ReducedModel ReducedModel;
        rowvec EstimatorIterations = "[1 2]";
        rowvec EstimatorVersions = "[1 1 0 0 1 0 0 1]";
        logical LogarithmicPlot = true;
        char SaveTexTables = "tables.tex";
        logical SortResultTable = false;
        logical UseOutputError = true;
        MarkerSize = 8;
        NumMarkers = 5;
        LineWidth = 2;
        ErrorOrders = "[1 2 3 5 10]";
        struct Est;
        PlotStartIndex = 1;
    public:

        ModelData;
    public:


        EstimatorAnalyzer(models.BaseFullModel model) {
            this.ModelData= struct(" Name ",[]," ErrT ",[]," RelErrT ",[]);
            if nargin == 1
                this.setModel(model);
            end
        }


        function  setModel(models.BaseFullModel model) {
            if isa(model," models.ReducedModel ")
                this.Model= model.FullModel;
                this.ReducedModel= model;
            else
                this.Model= model;

                model.offlineGenerations;
                this.ReducedModel= model.buildReducedModel;
            end
            /*  Only call struct building for old, non-deim estimators */
            if ~isa(this.ReducedModel.ErrorEstimator," error.DEIMEstimator ")
                this.buildEstimatorStruct(this.ReducedModel);
            end
        }
        function est = getDefaultEstStruct() {
            est = struct.empty;
            est(end+1).Name = " True error ";
            e = error.DefaultEstimator;
            e.offlineComputations(this.ReducedModel.FullModel);
            est(end).Estimator = e.prepareForReducedModel(this.ReducedModel);
            est(end).Estimator.Enabled= true;
            est(end).MarkerStyle = " o ";
            est(end).LineStyle = " - ";
            est(end).Color = [0 0 1];
        }
        function [errs , relerrs , ctimes , doublet , PlotManagerpm ] = compute(colvec mu,integer inidx,PlotManager pm) {

            if nargin < 3
                inidx = [];
                if nargin < 2
                    mu = [];
                end
            end

            [errs, ctimes, relerrs] = this.getErrorEstimates(mu, inidx, true);


            /*  MaxErr data */
            this.ModelData(end+1).Name = this.Model.Name;
            this.ModelData(end).ErrT = errs(:,end)^t;
            this.ModelData(end).MinErr = min(errs(:));
            this.ModelData(end).OverestT = (errs(:,end)/errs(1,end))^t;
            this.ModelData(end).RelErrT = relerrs(:,end)^t;
            this.ModelData(end).MinRelErr = min(relerrs(:));
            this.ModelData(end).CTimes = ctimes;
            this.ModelData(end).mu = mu;
            this.ModelData(end).inputidx = inidx;

            /*  Table overview */
            t = this.getResultTable(errs, ctimes);
            t.display;

            if ~isempty(this.SaveTexTables)
                t.Format= " tex ";
                t.saveToFile(this.SaveTexTables);
            end

            if nargin == 4
                this.createPlots(errs, relerrs, ctimes, pm);
            end
        }
        function [errs , ctimes , varargout ] = getErrorEstimates(colvec<double> mu,inidx,withrel) {
            if nargin == 3
                withrel = false;
            end
            num = length(this.Est);
            ctimes = zeros(1,num);
            errs = zeros(num,length(this.Model.Times));

            /* % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Simulations
             *fprintf('Starting estimator demo for model %s!\n',this.Model.Name);
             * Save old estimator */
            oldest = this.ReducedModel.ErrorEstimator;
            str = ; compplot = [];
            try
                for idx = 1:num
                    el = this.Est(idx);
                    /*  Set estimator and run */
                    this.ReducedModel.ErrorEstimator= el.Estimator;

                    /*  Optionally run some commands before starting stuff */
                    if isfield(el," Callback ") && ~isempty(el.Callback)
                        el.Callback(el.Estimator);
                    end

                    fprintf(" Performing estimation %d of %d: %s...\n ",idx,num,el.Name);
                    [~, ~, ctimes(idx)] = this.ReducedModel.simulate(mu,inidx);
                    if this.UseOutputError
                        errs(idx,:) = this.ReducedModel.ErrorEstimator.OutputError;
                    else
                        errs(idx,:) = this.ReducedModel.ErrorEstimator.StateError;
                    end
                    /*  Plotting preparations */
                    if ~isa(el.Estimator," error.ExpensiveBetaEstimator ")
                        str = [str sprintf(" errs(%d,end),ctimes(%d), "" %s "" , ",idx,idx,el.MarkerStyle)]; /* #ok<*AGROW> */

                        compplot(end+1) = idx;
                    end
                end
            catch ME
                /*  Restore old estimator */
                this.ReducedModel.ErrorEstimator= oldest;
                rethrow(ME);
            end
            /*  Restore old estimator */
            this.ReducedModel.ErrorEstimator= oldest;

            /*  Sanity checks */
            if any(abs(errs(1,1) - errs(2:end,1))/errs(1,1) > 1e-4)
                warning(" KerMor:EstimatorAnalyzer ",...
                    " Relative difference in initial errors over 1e-4 detected. ");
                figure; semilogy(errs(:,1));
            end

            if withrel
                varargout[1] = this.getRelativeErrorEstimates(errs, mu, inidx);
            end
        }


        function relerrs = getRelativeErrorEstimates(errs,colvec<double> mu,inidx) {
            if this.UseOutputError
                [~,yf] = this.Model.simulate(mu, inidx);
            else
                [~,yf] = this.Model.computeTrajectory(mu, inidx);
                yf = bsxfun(@times, yf, this.ReducedModel.System.StateScaling);
            end
            yfullnorm = sqrt(sum(yf.*(this.Model.G*yf),1));
            relerrs = errs ./ repmat(yfullnorm,size(errs,1),1);
        }


        function pm = createPlots(errs,relerrs,ctimes,pm) {
            if nargin < 5
                pm = PlotManager(false, 2, 2);
            end

            /*  Absolute error plots */
            this.plotErrors(errs, pm)

            /*  Relative error plots */
            this.plotRelativeErrors(relerrs, pm);

            /*  Computation ctimes plot */
            this.plotCTimes(errs, ctimes, pm);

            if nargout < 1
                pm.done;
            end
        }


        function  plotErrors(errs,pm) {
            ax = pm.nextPlot(" abserr ",[" Error estimations for model:  " this.Model.Name],...
                " Time ", " Error estimates ");

            this.doPlots(errs, ax);

            axis(ax,[this.Model.Times(this.PlotStartIndex) this.Model.T this.getYMin(errs) 3*max(errs(:))]);
        }


        function  plotRelativeErrors(relerrs,pm) {
            ax = pm.nextPlot(" relerr ",[" Relative error estimations e(t)/||y||, \Delta x(t)/||y|| for model:  " ...
                this.Model.Name], " Time ", " Relative error estimates ");

            this.doPlots(relerrs, ax);

            re = relerrs(:);
            re(isinf(re)) = -Inf;
            axis(ax,[this.Model.Times(this.PlotStartIndex) this.Model.T this.getYMin(relerrs) 3*max(re)]);
        }


        function  plotCTimes(errs,ctimes,pm) {
            ax = pm.nextPlot(" ctimes ",[" Error estimator computation times:  " this.Model.Name], " \Delta(T) ", " Comp. time [s] ");
            hold(ax, " on ");
            ph = zeros(length(this.Est),1);
            plotfun = @plot;
            if this.LogarithmicPlot
                plotfun = @semilogx;
            end
            ci = LineSpecIterator;
            for idx = 1:length(this.Est)
                e = this.Est(idx);
                if isfield(e," Color ")
                    if isempty(e.Color)
                        error(" If a field 'Color' exists, all estimators must have assigned a color. ");
                    end
                    c = e.Color;
                else
                    c = ci.nextColor;
                end
                ph(idx) = plotfun(ax,errs(idx,end),ctimes(idx),...
                            " Marker ",e.MarkerStyle," Color ",c,...
                            " MarkerSize ",this.MarkerSize*1.5);                
            end
            /*  Fill symbols (later, so char-color specs can be used in estimator structs) */
            for i=1:length(ph)
                c = get(ph(i)," Color ");
                set(ph(i)," MarkerFaceColor ",c);
                set(ph(i)," MarkerEdgeColor ",c*.7);
            end
            /*  Add line for minimum error! */
            plot(ax,[errs(1,end) errs(1,end)],[min(ctimes) max(ctimes)]," black ");
            hold(ax, " off ");
            /*  Add legend */
            a = cell(1,length(this.Est));
            [a[:]] = this.Est(:).Name;
            legend(a," Location "," NorthEast "," Interpreter "," latex ");
            axis([.9*min(errs(:,end)) 1.1*max(errs(:,end)) .9*min(ctimes(:)) 1.1*max(ctimes(:))]);
        }


        function pt = getResultTable(errs,ctimes) {
            if this.SortResultTable
                str = " Estimator hierarchy (error*time product) ";
                [~,idx] = sort(errs(:,end).*ctimes^t);
            else
                str = " Estimator list ";
                idx = 1:size(errs,2);
            end
            pt = PrintTable(" %s for model '%s' ",str,this.Model.Name);
            pt.HasRowHeader= true;
            pt.addRow(" Name ",sprintf(" $\\Delta(%g)$ ",this.Model.T)," Time "," Efficiency ");
            sfun = @(v)Utils.getLatexStr(v,3);
            for id = 1:length(this.Est)
                  pt.addRow(this.Est(idx(id)).Name,errs(idx(id),end),ctimes(idx(id)),...
                    errs(idx(id),end)/errs(1,end),[sfun," %2.2fs ",sfun]);
/*                  pt.addRow(this.Est(idx(id)).Name,errs(idx(id),end),ctimes(idx(id)),...
 *                     errs(idx(id),end)/errs(1,end),{'$%1.3$','%2.2fs','$%1.3e$'}); */
            end
            pt.HasHeader= true;
        }


        function ts = createStatsTables(sort) {
            m = length(this.ModelData);
            if nargin < 2
                sort = 1:m;
            end
            fields = [" ErrT "," RelErrT "," OverestT "," CTimes "];
            fieldnames = [" Errors "," Relative errors "," Overestimations "," Computation times "];

            for fi = 1:length(fields)
                t = PrintTable;
                t.Format= " tex ";
                t.HasHeader= true;
                t.Caption= sprintf(" %s of estimation runs ",fieldnames[fi]);
                t.addRow(" Model / Est ", this.Est(:).Name);

                fmt = " $%1.1e$ ";
                if strcmp(fields[fi]," CTimes ")
                    fmt = " $%2.2fs$ ";
                end
                fmt = [[" $%s$ "] repmat([fmt],1,length(this.Est))];
                for it=1:m
                    md = this.ModelData(sort(it));

                    /* % Compile name and add */
                    nstr = ;
                    if ~isempty(md.mu)
                        nstr = [" , \\mu=[ " sprintf(" %2.2e  ",md.mu) " ] "];
                    end
                    if ~isempty(md.inputidx)
                        nstr = [nstr sprintf(" , u_%d ",md.inputidx)];
                    end
                    hlp = num2cell(md.(fields[fi]));
                    t.addRow(nstr,hlp[:],fmt);
                end
                /* str = strrep(strrep(strrep(strrep(t.print,'e+','e^{'),'e-','e^{-'),'{0','{'),'{-0','{-');
                 *fprintf(fid,'%s',str);
                 *t.print(fid); */
                ts(fi) = t;
            end
            /* fclose(fid); */
        }
        /* % Getter & Setter */
    public:

        
#if 0 //mtoc++: 'set.EstimatorIterations'
function  EstimatorIterations(value) {
            this.EstimatorIterations= value;
            if ~isempty(this.ReducedModel)/* #ok */

                this.buildEstimatorStruct(this.ReducedModel);/* #ok */

            end
        }

#endif


        
#if 0 //mtoc++: 'set.EstimatorVersions'
function  EstimatorVersions(value) {
            this.EstimatorVersions= value;
            if ~isempty(this.ReducedModel)/* #ok */

                this.buildEstimatorStruct(this.ReducedModel);/* #ok */

            end
        }

#endif

    
    private:


        function matrix<double>y = getYMin(err) {
            err = err(:,this.PlotStartIndex:end);
            e1 = min(err(:,1));
            er = min(reshape(err(:,2:end),1,[]));
            if log10(abs(e1 - er)) < 5
                y = .5*er;
            else
                y = .5*min(err(:));
            end
        }


        function  doPlots(data,ax) {
            times = this.Model.Times(this.PlotStartIndex:end);
            data = data(:,this.PlotStartIndex:end);
            if this.LogarithmicPlot
                ph = semilogy(ax,times,data);
            else
                ph = plot(ax,times,data);
            end
            set(ph," LineWidth ",this.LineWidth);
            set(ph(1)," LineWidth ",this.LineWidth+.5);
            hold on;
            /*  Select extra marker places */
            nt = length(times);
            sel = round(1:nt/this.NumMarkers:nt);
            ci = LineSpecIterator;
            for idx=1:length(this.Est)
                e = this.Est(idx);
                if isfield(e," Color ")
                    if isempty(e.Color)
                        error(" If a field 'Color' exists, all estimators must have assigned a color. ");
                    end
                    c = e.Color;
                else
                    c = ci.nextColor;
                end
                set(ph(idx)," LineStyle ",e.LineStyle," Color ",c);
                /*  Shift marker positions for better visual */
                pos = mod(sel+round(nt/(this.NumMarkers*length(this.Est)))*(idx-1),nt)+1;
                h = plot(ax,times(pos), data(idx,pos), e.MarkerStyle,...
                    " MarkerSize ",this.MarkerSize);
                /*  Get resulting color, so char-color specs can be used as e.Color above) */
                c = get(ph(idx)," Color ");
                set(h," MarkerFaceColor ",c," MarkerEdgeColor ",c*.7);
            end
            a = cell(1,length(this.Est));
            [a[:]] = this.Est(:).Name;
            [~,oh] = legend(a," Location "," NorthEast "," Interpreter "," latex "," FontSize ",14);
            /*  Assign markers to legend */
            oh = findobj(oh," Type "," line ");
            for idx=1:length(this.Est)
                c = get(ph(idx)," Color ");
                set(oh(2*idx)," Marker ",this.Est(idx).MarkerStyle," MarkerFaceColor ",c,...
                    " MarkerEdgeColor ",c*.7," MarkerSize ",this.MarkerSize);
            end
        }
        function  buildEstimatorStruct(r) {
            this.Est= this.buildKernelEstimatorStruct(r);
        }


        function est = buildKernelEstimatorStruct(r) {
            est = struct.empty;

            /*  Default Estimator (true error) */
            if this.EstimatorVersions(1)
                est(end+1).Name = " True error ";
                est(end).Estimator = error.DefaultEstimator(r);
                est(end).Estimator.Enabled= true;
                est(end).MarkerStyle = " o ";
                est(end).LineStyle = " - ";
            end

            /*  GLE estimator */
            if this.EstimatorVersions(2)
                msg = error.GLEstimator.validModelForEstimator(r.FullModel);
                e = error.GLEstimator;
                e.offlineComputations(r.FullModel);
                e = e.prepareForReducedModel(r);
                if isempty(msg)
                    fprintf(" Initializing Global Lipschitz estimator...\n ");
                    est(end+1).Name = " GLE ";
                    est(end).Estimator = e;
                    est(end).MarkerStyle = " s ";
                    est(end).LineStyle = " - ";
                else
                    fprintf(" Cannot use the GLEstimator for model %s:\n%s\n ",r.Name,msg);
                end
            end

            msg = error.IterationCompLemmaEstimator.validModelForEstimator(r.FullModel);
            if isempty(msg)

                if ~isempty(this.Model.Approx)
                    k = this.Model.Approx.Expansion.Kernel;
                else
                    k = this.Model.System.f.Expansion.Kernel;
                end

                reps = this.EstimatorIterations;
                fprintf(" Using iteration counts: %s\n ",num2str(this.EstimatorIterations));
                e = error.IterationCompLemmaEstimator;
                e.offlineComputations(r.FullModel);
                e = e.prepareForReducedModel(r);
                if this.EstimatorVersions(3)
                    fprintf(" Initializing LGL estimator...\n ");
                    est(end+1).Name = " LGL ";
                    est(end).Estimator = e;
                    est(end).Estimator.LocalLipschitzFcn= error.lipfun.LocalGradientLipschitz(k);
                    est(end).Estimator.UseTimeDiscreteC= false;
                    est(end).MarkerStyle = " s ";
                    est(end).LineStyle = " - ";
                    for it = reps
                        orig = est(end).Estimator;/* #ok */

                        eval(sprintf([" est(end+1).Name =  "" LGL, %d It "" ; "...
                            " est(end).Estimator = orig.clone; "...
                            " est(end).Estimator.Iterations = %d; "...
                            " est(end).MarkerStyle =  "" s "" ; "...
                            " est(end).LineStyle =  "" -. "" ; "],it,it));
                    end
                end

                if this.EstimatorVersions(4)
                    fprintf(" Initializing LGLMod (mod secant) estimator...\n ");
                    est(end+1).Name = " LGLMod ";
                    est(end).Estimator = e.clone;
                    est(end).Estimator.LocalLipschitzFcn= error.lipfun.LocalSecantLipschitz(k);
                    est(end).Estimator.UseTimeDiscreteC= false;
                    est(end).MarkerStyle = " h ";
                    est(end).LineStyle = " - ";
                    for it = reps
                        orig = est(end).Estimator;/* #ok */

                        eval(sprintf([" est(end+1).Name =  "" LGLMod, %d It "" ; "...
                            " est(end).Estimator = orig.clone; "...
                            " est(end).Estimator.Iterations = %d; "...
                            " est(end).MarkerStyle =  "" h "" ; "...
                            " est(end).LineStyle =  "" : "" ; "],it,it));
                    end
                end

                if any(this.EstimatorVersions([5 8]))
                    ilfcn = error.lipfun.ImprovedLocalSecantLipschitz(k);
                end

                if this.EstimatorVersions(5)
                    fprintf(" Initializing LSL estimator...\n ");
                    est(end+1).Name = " LSLE ";
                    est(end).Estimator = e.clone;
                    est(end).Estimator.LocalLipschitzFcn= ilfcn.clone;
                    est(end).Estimator.UseTimeDiscreteC= false;
                    est(end).MarkerStyle = " p ";
                    est(end).LineStyle = " - ";
                    for it = reps
                        orig = est(end).Estimator;/* #ok */

                        eval(sprintf([" est(end+1).Name =  "" LSLE, %d It "" ; "...
                            " est(end).Estimator = orig.clone; "...
                            " est(end).Estimator.Iterations = %d; "...
                            " est(end).MarkerStyle =  "" p "" ; "...
                            " est(end).LineStyle =  "" -- "" ; "],it,it));
                    end
                end

                /* % TD-Versions */
                if any(this.EstimatorVersions(6:8))
                    td = error.IterationCompLemmaEstimator;
                    td.offlineComputations(r.FullModel);
                    td = td.prepareForReducedModel(r);
                    td.Iterations= 0;
                    td.UseTimeDiscreteC= true;
                end

                if this.EstimatorVersions(6)
                    fprintf(" Initializing LSL TD estimators...\n ");
                    est(end+1).Name = " LGL TD ";
                    est(end).Estimator = td;
                    est(end).Estimator.LocalLipschitzFcn= error.lipfun.LocalGradientLipschitz(k);
                    est(end).MarkerStyle = " < ";
                    est(end).LineStyle = " - ";
                end

                if this.EstimatorVersions(7)
                    fprintf(" Initializing LSL TD estimators...\n ");
                    est(end+1).Name = " LGL TD ";
                    est(end).Estimator = td.clone;
                    est(end).Estimator.LocalLipschitzFcn= error.lipfun.LocalSecantLipschitz(k);
                    est(end).MarkerStyle = " < ";
                    est(end).LineStyle = " - ";
                end

                if this.EstimatorVersions(8)
                    est(end+1).Estimator = est(end).Estimator.clone;
                    est(end).Name = " LSLE TD ";
                    est(end).Estimator = td.clone;
                    est(end).Estimator.LocalLipschitzFcn= ilfcn.clone;
                    est(end).MarkerStyle = " < ";
                    est(end).LineStyle = " - ";
                end
            else
                fprintf(" Cannot use the IterationCompLemmaEstimator for model %s:\n%s\n ",r.Name,msg);
            end
        }
};