LogNorm.m

Go to the documentation of this file.

namespace testing{


/* (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 LogNorm {
    public: /* ( Static ) */


        static function [aln , times , rowvec<integer>st_sizes ] = compareSimTransJac_FullJac(models.BaseFullModel m,rowvec<integer> st_sizes) {

            /*  Get required data */
            jtd = m.Data.JacobianTrainData;
            QFull = m.ErrorEstimator.QFull;
            f = m.System.f;

            /*  Input checks */
            if nargin < 2
                st_sizes = 1:size(QFull,2);
            end

            /*  Preps */
            stn = length(st_sizes);
            n = size(jtd.xi,2);
            aln = zeros(stn,n);
            times = zeros(stn,n);
            pi = ProcessIndicator(" Computing approximated log norms for %d full jacobians and %d sim. trans. orders ",...
                n*stn,false,n,stn);
            for nr = 1:n
                J = f.getStateJacobian(jtd.xi(:,nr),jtd.ti(nr),jtd.mui(:,nr));
                for snr = 1:stn
                    Qk = QFull(:,1:st_sizes(snr));
                    t = tic;
                    aln(snr,nr) = Utils.logNorm(Qk^t*J*Qk);
                    times(snr,nr) = toc(t);
                    pi.step;
                end
            end
            pi.stop;
            if nargout == 0
                testing.DEIM.compareSimTransJac_FullJac_plots(m, aln, times, st_sizes);
            end
        }
        static function pm = compareSimTransJac_FullJac_plots(m,aln,times,st_sizes,pm) {
            if nargin < 5
                pm = PlotManager(false,2,2);
            end
            jstd = m.Data.JacSimTransData;
            pm.nextPlot(" correct_ln ",sprintf(" Full logarithmic norms\nJacobian full dimension %dx%d ",...
                m.System.f.fDim,m.System.f.xDim));
            plot(jstd.LogNorms);

            if size(aln,1) > 1
                [X,Y] = meshgrid(1:size(aln,2),st_sizes);
                LN = repmat(jstd.LogNorms,size(aln,1),1);
            end

            if size(aln,1) == 1
                h = pm.nextPlot(" approx_ln ",...
                    sprintf(" Approx logarithmic norms\nSimilarity transformation size %d ",...
                    st_sizes));
                plot(h, aln);
            else
                h = pm.nextPlot(" approx_ln "," Approx logarithmic norms ",...
                    " training sample "," similarity transform size ");
                LogPlot.nicesurf(h, X, Y, aln);
            end

            h = pm.nextPlot(" rel_err "," Relative error "," training sample ",...
                " similarity transform size ");
            LogPlot.logsurf(h,X,Y,abs((LN-aln)./LN));

            sv = m.ErrorEstimator.QSingVals;
            s = m.ErrorEstimator.JacSimTransMaxSize;
            h = pm.nextPlot(" sing_vals ",...
                sprintf(" Singular values of SVD of eigenvector matrix\nBlack line: Max size of sim. trans., singular value= %g ",...
                sv(s))," POD Modes "," Singular values ");
            semilogy(h,sv);
            hold on;
            plot(h,[s s+eps],[min(sv) max(sv)]," k ");
            hold off;

            h = pm.nextPlot(" comp_times ",...
                sprintf(" Average computation times (over %d values)\nfor log norms of similarity transformed matrices ",...
                    size(times,2))," similarity transformation size "," mean computation time ");
            me = mean(times,2);
            semilogy(h,st_sizes,me," -s ");
            if nargout < 1
                pm.done;
            end
        }


                        static function [aln , times , jtimes , rowvec<integer>deim_orders , rowvec<integer>st_sizes ] = compareSimTransDEIMJac_FullJac(models.BaseFullModel m,rowvec<integer> deim_orders,rowvec<integer> st_sizes) {

            /*  Get required data */
            jtd = m.Data.JacobianTrainData;
            /* zi = m.Data.W'*jtd.xi; */
            e = m.ErrorEstimator;
            jd = e.JacMDEIM;

            /*  Input checks */
            if nargin < 3
                st_sizes = [1:e.JacSimTransMaxSize 0];
                if nargin < 2
                    deim_orders = 1:jd.MaxOrder;
                end
            end

            /*  Preps */
            stn = length(st_sizes);
            n = size(jtd.xi,2);
            no = length(deim_orders);

            aln = zeros(no,stn,n);
            times = aln;
            jtimes = aln;
            pi = ProcessIndicator(" Computing approximated log norms over %d DEIM orders and %d sim. trans. sizes on %d training values ",...
                numel(aln),false,no,stn,n);
            for onr = 1:no
                jd.Order= deim_orders(onr);
                for snr = 1:stn
                    jd.setSimilarityTransform(e.QFull(:,1:st_sizes(snr)));
                    for nr = 1:n
                        t = tic;
                        J = jd.evaluate(jtd.xi(:,nr),jtd.ti(nr),jtd.mui(:,nr));
                        jtimes(onr,snr,nr) = toc(t);
                        t = tic;
                        aln(onr,snr,nr) = Utils.logNorm(J);
                        times(onr,snr,nr) = toc(t);
                        pi.step;
                    end
                end
            end
            pi.stop;
        }
                        static function pm = compareSimTransDEIMJac_FullJac_plots(m,aln,times,jtimes,deim_orders,st_sizes,pm) {
            if nargin < 7
                pm = PlotManager(false,2,2);
            end

            jstd = m.Data.JacSimTransData;
            n = size(aln,3);
            ln = reshape(jstd.LogNorms,1,1,[]);
            ln = repmat(ln,[length(deim_orders), length(st_sizes), 1]);

            abserr = abs(aln-ln);
            relerr = abs(abserr ./ ln);
            relerr(isnan(relerr)) = 0;

            /*  Aggregate errors over trajectory sample data
             *aln = sqrt(sum(aln.^2,3)); */
            abserr = mean(abserr,3);
            relerr = mean(relerr,3);

            /*  Mean the times too */
            ttimes = times + jtimes;
            times = mean(times,3);
            jtimes = mean(jtimes,3);
            ttimes = mean(ttimes,3);

            [X,Y] = meshgrid(st_sizes,deim_orders);

            h = pm.nextPlot(" abs_err ",sprintf(" Mean absolute approximation error over %d samples ",n),...
                " Similarity transformation size "," DEIM order ");
            LogPlot.logsurf(h,X,Y,abserr);

            h = pm.nextPlot(" rel_err ",sprintf(" Mean relative approximation error over %d samples ",n),...
                " Similarity transformation size "," DEIM order ");
            LogPlot.logsurf(h,X,Y,relerr);

            h = pm.nextPlot(" comp_times_j ",...
                sprintf(" Average computation times over %d values\nfor matrix DEIM jacobian evaluation ",...
                    n)," Similarity transformation size "," DEIM order ");
            LogPlot.nicesurf(h,X,Y,jtimes);

            h = pm.nextPlot(" comp_times ",...
                sprintf(" Average computation times over %d values\nfor log norm computation of sim.trans. matrix DEIM jacobian ",...
                    n)," Similarity transformation size "," DEIM order ");
            LogPlot.nicesurf(h,X,Y,times);

            h = pm.nextPlot(" comp_times_total ",...
                sprintf(" Average total computation times over %d values\n(jac comp + log norm comp) ",...
                    n)," Similarity transformation size "," DEIM order ");
            LogPlot.nicesurf(h,X,Y,ttimes);

            if nargout < 1
                pm.done;
            end
        }
        static function structres = getApproxLogNormsAtPos(models.BaseFullModel mo,colvec<double> x,double t,matrix<double> mui) {

            /*  Get required data */
            e = mo.ErrorEstimator;
            if isa(mo," models.ReducedModel ")
                fm = mo.FullModel;
            else
                fm = mo;
            end
            jd = e.JacMDEIM;

            /*  Input checks */
            if nargin < 4
                mui = fm.Data.ParamSamples;
            end

            /*  Preps */
            n = size(x,2);
            if ~isempty(t)
                if length(t) == 1
                    t = ones(1,n)*t;
                elseif length(t) ~= n
                    error(" Size of t vector must match column count of x ");
                end
            end
            m = size(mui,2);
            aln = zeros(n,m);
            times = aln;
            jtimes = aln;
            pi = ProcessIndicator(" Computing approximated log norms for %d parameter values at %d locations ",n*m,false,m,n);
            for i = 1:n
                for j = 1:m
                    tc = tic;
                    J = jd.evaluate(x(:,i),t(i),mui(:,j));
                    jtimes(i,j) = toc(tc);
                    tc = tic;
                    aln(i,j) = Utils.logNorm(J);
                    times(i,j) = toc(tc);
                    pi.step;
                end
            end
            pi.stop;
            res.aln= aln;
            res.times= times;
            res.jtimes= jtimes;
            res.mJ= e.JacMatDEIMOrder;
            res.k= e.JacSimTransSize;
            res.mui= mui;
            res.munames= [mo.System.Params(:).Name];
        }
        static function  getApproxLogNormsAtPos_plots(res,pm) {
            str = sprintf(" \nm_J = %d, k=%d, %d locations ",res.mJ, res.k, size(res.aln,1));
            if size(res.mui,1) == 1
                if size(res.aln,1) == 1
                    h = pm.nextPlot(" aln_min ",[" Log norms at given location " str],res.munames[1]);
                    LogPlot.cleverPlot(h,res.mui,res.aln," b ");
                else
                    [X,Y] = meshgrid(res.mui,1:size(res.aln,1));
                    h = pm.nextPlot(" aln_min ",[" Log norms over all locations " str],res.munames[1]," location nr ");
                    surf(h,X,Y,res.aln," EdgeColor "," none "," FaceColor "," interp ");
                end
            elseif size(res.mui,1) == 2
                tri = delaunay(res.mui(1,:),res.mui(2,:));
                if size(res.aln,1) == 1
                    h = pm.nextPlot(" aln_min ",[" Log norms at given location " str],res.munames[:]);
                    trisurf(tri,res.mui(1,:),res.mui(2,:),res.aln,...
                        " EdgeColor "," none "," FaceColor "," interp "," Parent ",h);
                else
                    /* % Animation
 *                     h = pm.nextPlot('aln','Approximated Jacobian log norms over sample parameters');
 *                     for k=1:size(res.aln,1)
 *                         trisurf(tri,res.mui(1,:),res.mui(2,:),res.aln(k,:),...
 *                             'EdgeColor','none','FaceColor','interp','Parent',h);
 *                         title(h,sprintf('Approximated Jacobian log norms over sample parameters, pos: %d',k));
 *                         view(90,0);
 *                         pause;
 *                     end
                     *% Min/Max plot */
                    alnmin = min(res.aln,[],1);
                    alnmax = max(res.aln,[],1);
                    h = pm.nextPlot(" aln_min ",[" Minimum log norms over all locations " str],res.munames[:]);
                    trisurf(tri,res.mui(1,:),res.mui(2,:),alnmin,...
                        " EdgeColor "," none "," FaceColor "," interp "," Parent ",h);
                    h = pm.nextPlot(" aln_min ",[" Max log norms over all locations " str],res.munames[:]);
                    trisurf(tri,res.mui(1,:),res.mui(2,:),alnmax,...
                        " EdgeColor "," none "," FaceColor "," interp "," Parent ",h);
                    h = pm.nextPlot(" aln_min ",[" Max-Min log norms over all locations " str],res.munames[:]);
                    trisurf(tri,res.mui(1,:),res.mui(2,:),alnmax-alnmin,...
                        " EdgeColor "," none "," FaceColor "," interp "," Parent ",h);
                end
            else
                error(" Cannot plot results for more than two parameters ");
            end
        }


        static function [res , mScale , MScale , pos , l , sel , seli ] = CompLogNorms(m,numt) {


            atd = m.Data.ApproxTrainData;
            N = size(atd.xi,2);

            if nargin < 2
                numt = N;
            else
                numt = min(numt, N);
            end
            seed = 1;
            /*  Show plots for each numt value */
            doplot = 4==5;
            /*  Sort training points after norm values (ascending) */
            dosort = 4==5;
            /*  Plot negative local logarithmic norms */
            plotneg = 4==4;

            r = RandStream(" mt19937ar "," Seed ",seed);

            /* sel = unique(r.randi(N,1,num)); */
            sel = 1:N;
            seli = unique(r.randi(N,1,numt));
            n = length(seli);

            pi = ProcessIndicator(" Computing %d local log norms, maxed over %d training points ",n,false,n,length(sel));
            ln = zeros(1,n);
            l = ln;
            pos = l;

            if doplot
                pm = PlotManager(false, 2, 2);
            end

            xi = atd.xi(:,sel);
            fxi = atd.fxi(:,sel);
            V = m.Data.V;
            W = m.Data.W;

            /* % Compute constant terms */
            xinsq = sum(xi.^2);

            /*  Optional sorting of values */
            if doplot && dosort
                [xinsq sortidx] = sort(xinsq);
                xi = xi(:,sortidx);
                fxi = fxi(:,sortidx);
            end

            xifxi = sum(xi.*fxi);
            Vfxi = V^t*fxi;
            Vxi = V^t*xi;
            VV = 1;/* V'*V;
             * if norm(VV) - 1 < 1e-12
             *     VV = 1;
             * end */


            /*  Locally selected "Vz" terms */
            zi = W^t*atd.xi(:,seli);
            MU = atd.mui(:,seli);
            fVzi = m.System.f.evaluate(V*zi,atd.ti(seli),MU);

            /*  Only needed for efficient local lipschitz constant computation
             * fxinsq = sum(fxi.*fxi); */

            /* % Main loop */
            for i = 1:n
                z = zi(:,i);
                fz = fVzi(:,i);
                nom = xifxi - z" *Vfxi + z "*(V^t*fz);
                nom = nom - fz^t*xi;
                denom = xinsq - 2*z" *Vxi + z "*(VV*z);
                Ln = nom ./ denom;
                /* denom(zer) = sum((xi(:,zer)-V*z(:,ones(1, length(zer)))).^2); */

            /*      L = sqrt(abs((fxinsq - 2*fz'*fxi + fz'*fz) ./ denom));
             *     L2 = Norm.L2(fxi - fVzi(:,ones(1,length(sel))*i)) ./ sqrt(denom); */

                [~, idx] = sort(denom);
                Ln = Ln(idx);

                zer = denom == 0;
                Ln(zer) = [];

                [ln(i), pos(i)] = max(Ln);

            /*      [~, sortidx] = sort(denom);
             *     Ln = Ln(sortidx); */
                if doplot
                    ti = sprintf(" Idx %d: loc Lip const %g, loc log norm %g ",seli(i),l(i),ln(i));
                    ax = pm.nextPlot(" blah ",ti," training point "," local [lip const (b), log norm (g)] ");
                    neg = Ln <= 0;
                    if plotneg
                        Lnn = -Ln;
                    end
                    Ln(neg) = 0;
                    hlp = [L; Ln];
                    if plotneg
                        Lnn(~neg) = 0;
                        hlp = [hlp; Lnn];/* #ok */

                    end
                    semilogy(ax, hlp^t);
                end
                pi.step;
            end
            if doplot
                pm.done;
            end
            pi.stop;

            res = data.ApproxTrainData(zi,[],[]);
            res.fxi= ln;        
            res = data.ApproxTrainData.makeUniqueXi(res);
            [res, mScale, MScale] = data.ApproxTrainData.scaleXiZeroOne(res);
        }
};
}