PODGreedy.m

Go to the documentation of this file.

namespace spacereduction{


/* (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 PODGreedy
  :public spacereduction.BaseSpaceReducer,
   public IParallelizable {
    public: /* ( setObservable ) */

        double Eps = 1e-6;
        double MinRelImprovement = .005;
        integer MaxSubspaceSize = 100;
        matrix<double> InitialSpace = "[]";
    public:

        ErrorImprovements;

        
    public:

        PODGreedy() {
            this.registerProps(" Eps "," MinRelImprovement "," MaxSubspaceSize "," InitialSpace ");
        }


        function  plotSummary(pm,context) {
            if nargin < 2
                pm = PlotManager;
                pm.LeaveOpen= true;
                context = " POD-Greedy ";
            end
            plotSummary@spacereduction.BaseSpaceReducer(this, pm, context);
            if ~isempty(this.ErrorImprovements)
                str = sprintf(" %s: Error decay over training data\nEps:%g, MaxSize:%d, FXI: %d, FD: %d, BSpan: %d ",...
                    context,this.Eps,this.MaxSubspaceSize,this.IncludeTrajectoryFxiData,...
                    this.IncludeFiniteDifferences,this.IncludeBSpan);
                h = pm.nextPlot(" podgreedy_gains ",str," subspace size "," Gain ");
                semilogy(h,this.ErrorImprovements," LineWidth ",2);
            else
                warning(" spacereduction:PODGreedy ",...
                    " Error data empty. Not providing summary. ");
            end
        }


    
    protected:


        function [V , W ] = generateReducedSpaceImpl(models.BaseFullModel model,subset) {
            bdata = model.Data.TrajectoryData;

            /*  Wrap in finite difference adder */
            if this.IncludeFiniteDifferences
                bdata = data.JoinedBlockData(bdata, data.FinDiffBlockData(bdata));
            end
            if this.IncludeAxData
                bdata = data.JoinedBlockData(bdata, ...
                    data.AxBlockData(bdata, model.System.A, model.scaledTimes));
            end
            /*  Augment block data with fxi values */
            if this.IncludeTrajectoryFxiData
                if isempty(model.Data.TrajectoryFxiData)
                    error(" No training fxi data found in ModelData. ");
                end
                bdata = data.JoinedBlockData(bdata, model.Data.TrajectoryFxiData);
            end

            if KerMor.App.Verbose > 2
                fprintf(" POD-Greedy: Starting subspace computation using %d trajectories...\n ",model.Data.TrajectoryData.getNumTrajectories);
            end

            pod = general.POD;
            pod.Value= 1; /*  MUST stay at 1 or getInitialSpace will fail. */

            pod.Mode= " abs ";

            o = general.Orthonormalizer;
            o.Algorithm= " gs ";

            /*  Compute initial space */
            if KerMor.App.Verbose > 2
                fprintf(" POD-Greedy: Computing initial space...\n ");
            end
            if ~isempty(this.InitialSpace)
                V = this.InitialSpace(subset,:);
                if this.IncludeBSpan
                    V = o.orthonormalize([V bdata.InputSpaceSpan(subset,:)]);
                end
            elseif this.IncludeBSpan
                V = model.Data.InputSpaceSpan(subset,:);
            else
                V = this.getInitialSpace(bdata, pod, subset);
            end
            /*  Compute initial space errors */
            err = [];
            for k=1:size(V,2)
               [err(end+1), idx] = this.getMaxErr(V(:,1:k), bdata, subset);/* #ok */

            end

            cnt = 1; 
            /*  Maximum possible subspace size */
            ss = size(V,1);
            olderr = err;
            impr = 1;
            while (err(end) > this.Eps) && cnt < this.MaxSubspaceSize && size(V,2) < ss && impr(end) > this.MinRelImprovement
                x = bdata.getBlock(idx); /*  get trajectory */

                x = x(subset,:);
                e = x - V*(V^t*x);
                Vn = pod.computePOD(e);
                V = o.orthonormalize([V Vn]);
                [err(end+1), idx] = this.getMaxErr(V, bdata, subset);/* #ok */

                impr(end+1) = (olderr-err(end))/olderr;/* #ok */

                olderr = err(end);
                cnt = cnt+1;
            end
            this.ErrorImprovements= impr;
            this.ProjectionError= err;
            if cnt == this.MaxSubspaceSize
                fprintf(" POD-Greedy: Maximum predefined subspace size %d reached. Aborting.\n ",this.MaxSubspaceSize);
            end
            if size(V,2) == ss
                fprintf(" POD-Greedy: Maximum possible subspace size %d reached. Aborting.\n ",ss);
            end
            if impr(end) <= this.MinRelImprovement
                fprintf(" POD-Greedy: Minimum relative improvement of %f not satisfied. Aborting.\n ",this.MinRelImprovement);
            end
            if KerMor.App.Verbose > 1
                fprintf(" POD-Greedy: Finished with subspace size %d and error %e.\n ",size(V,2),err(end));
                if KerMor.App.Verbose > 2
                    figure; subplot(1,2,1);
                    semilogy(err);
                    xlabel(" POD-Greedy step "); ylabel(" Error norms ");
                    title(" L_2 norm sum over whole projection error trajectory vectors ");
                    subplot(1,2,2);
                    plot(impr);
                    xlabel(" POD-Greedy step "); ylabel(" Relative improvement ");
                    title(" Error improvement relative to previous error ");
                end
            end
            /*  Galerkin projection! */
            W = [];
        }



    
    private:


        function [maxerr , midx ] = getMaxErr(V,md,reducable) {
            midx = -1;
            maxerr = 0;
            if this.ComputeParallel
                n = md.getNumBlocks;
                if KerMor.App.Verbose > 3
                   fprintf(" POD-Greedy: Computing maximum error over %d trajectories on %d workers for subspace size %d...\n ",n,matlabpool(" size "),size(V,2));
                end
                err = zeros(1,n);
                parfor i=1:n
                    x = md.getBlock(i); /* #ok<PFBNS> */

                    x = x(reducable,:);
                    hlp = sum(Norm.L2(x - V*(V^t*x)));
                    if KerMor.App.Verbose > 4
                       fprintf(" POD-Greedy: Error for block %d: %e\n ",i,hlp);
                    end
                    err(i) = hlp;
                end
                [maxerr, midx] = max(err);
            else
                for k=1:md.getNumBlocks;
                    x = md.getBlock(k);
                    x = x(reducable,:);
                    e = sum(Norm.L2(x - V*(V^t*x)));
                    if maxerr < e
                        maxerr = e;
                        midx = k;
                    end
                end
            end
            if KerMor.App.Verbose > 0
                fprintf(" POD-Greedy: Max error with subspace size %d: %e\n ",size(V,2),maxerr);
            end
        }

    
    public: /* ( Static ) */

        static function  test_PODGreedy() {
            
            sv = 15;
            dim = 100;
            dt = .05;
            T = 1;

            m = models.BaseFullModel;
            m.dt= dt;
            m.T= T;

            m.Sampler= sampling.GridSampler;

            m.Approx= [];

            tg = spacereduction.PODGreedy;
            tg.Eps= 1e-5;
            m.SpaceReducer= tg;

            m.ODESolver= solvers.ExplEuler(dt);

            s = models.BaseFirstOrderSystem(m);
            s.B= [];
            s.addParam(" mu1 ",.5," Range ",[0 1]," Desired ",4);
            s.addParam(" mu2 ",.6," Range ",[0 1]," Desired ",5);
            s.x0= dscomponents.ConstInitialValue(rand(dim,1)*5);
            s.MaxTimestep= dt;
            m.System= s;

            f = dscomponents.ParamTimeKernelCoreFun(m.System);
            kexp = kernels.ParamTimeKernelExpansion;
            kexp.Kernel= kernels.GaussKernel(40);
            kexp.TimeKernel= kernels.NoKernel;
            kexp.ParamKernel= kernels.GaussKernel(2);
            kexp.Ma= rand(dim,sv);
            kexp.Centers.xi= repmat(linspace(0,10,sv),dim,1);
            kexp.Centers.ti= 1:sv;
            kexp.Centers.mui= rand(2,sv);
            f.Expansion= kexp;
            s.f= f;

            m.offlineGenerations;
        }

    

};
}