POD.m

Go to the documentation of this file.

namespace general{


/* (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 POD
  :public KerMorObject {
    public: /* ( setObservable ) */

        Mode = "rel";
        Value = .3;
        UseSVDS = false;
    public: /* ( setObservable ) */


        POD() {
            this = this@KerMorObject;

            this.registerProps(" Mode "," Value "," UseSVDS ");
        }


        function [matrix<double>podvec , rowvec<double>s ] = computePOD(matrix<double> data,matrix<double> Vexclude,colvec<integer> targetdims) {

            /* % Dimension checks
             * This is for the "explicit" modes where the target dimension
             * is either a fixed value or a fraction of the full dimension. */
            if nargin < 4
                targetdims = " : ";
                if nargin < 3
                    Vexclude = [];
                end
            end
            target_dim=[];
            if strcmpi(this.Mode," rel ")
                target_dim = ceil(min(size(data))*this.Value);
            elseif strcmpi(this.Mode," abs ")
                target_dim = this.Value;
            end
            if ~isempty(target_dim) && target_dim > min(size(data))
                warning(" KerMor:general:POD "," Targeted reduced space dimension (%d) has to be <= the smallest matrix dimension (%d). Using size %d. ",...
                    target_dim,min(size(data)),min(size(data)));
                target_dim = min(size(data));
            end

            /* % Block SVD case */
            if isa(data," data.ABlockedData ")
                if ~any(strcmpi(this.Mode,[" abs "," rel "]))
                    target_dim = min(size(data));
                    if max(size(data)) > 10000
                        warning(" POD on data.ABlockedData with dimension %d and mode '%s' is possibly very very slow. Try to use  "" abs ""  or  "" rel ""  modes. ",max(size(data)),this.Mode);
                    end
                end
                if KerMor.App.Verbose > 2
                    [n,m] = size(data);
                    fprintf(" Starting POD on ABlockedData(%s, %dx%d, %d blocks) with mode  "" %s ""  and value %g\n ",...
                        class(data),n,m,data.getNumBlocks,this.Mode,this.Value);
                end
                [podvec, s] = data.getSVD(target_dim, Vexclude, targetdims);
                s = diag(s);

                sig = 1:size(podvec,2);
                if strcmpi(this.Mode," sign ")
                    sig = s >= s(1)*this.Value;
                elseif strcmpi(this.Mode," eps ")
                    sig = s >= this.Value;
                end
                podvec = podvec(:,sig);

            /* % Matrix argument case */
            else
                /*  Reduce the data to the targeted dimensions first */
                data = data(targetdims,:);

                /*  Subtract the space spanned by Vexclude if given */
                if ~isempty(Vexclude)
                    data = data - Vexclude*(Vexclude^t*data);
                    target_dim = min(target_dim, min(size(data))-size(Vexclude,2));
                end

                /*  Create full matrix out of sparse fxi sets to enable use of
                 * svd instead of svds */
                nonzero = [];
                if issparse(data)
                    nonzero = find(sum(abs(data),2) ~= 0);
                    /*  Check if same sparsity pattern holds for each fxi column */
                    if ~isempty(nonzero)
                        if KerMor.App.Verbose > 2
                            fprintf(" Reducing data dimension from %d to %d as sparse matrix is given. ",...
                                size(data,1),length(find(nonzero)));
                        end
                        fulldim = size(data,1);
                        data = full(data(nonzero,:));
                    end
                end

                if KerMor.App.Verbose > 2
                    fprintf(" Starting POD on %dx%d matrix with mode  "" %s ""  and value %g (UseSVDS=%d) ...  ",...
                        size(data),this.Mode,this.Value,this.UseSVDS);
                end
                /* % Reduction for modes 'sign' and 'eps'
                 * For these "dynamic" modes the full singular values have
                 * to be computed in order to determine how many to use. */
                if any(strcmpi(this.Mode,[" sign "," eps "]))
                    [u,s] = svd(data," econ ");
                    s = diag(s);
                    if strcmpi(this.Mode," sign ")
                        sig = s >= s(1)*this.Value;
                    elseif strcmpi(this.Mode," eps ")
                        sig = s >= this.Value;
                    end
                    if KerMor.App.Verbose > 2
                        fprintf(" selecting %d singular values ...  ",length(find(sig)));
                    end
                    /*  Select wanted subspace */
                    u = u(:,sig);

                    err = sum(s(~sig));
                    rerr = 100*err/sum(s);
                else
                    /* % Reduction for modes 'abs' and 'rel'
                     * For cases 'abs' or 'rel': fixed target dimension.
                     * So just let svds extract the wanted components! */

                    /*  As tests showed that the svds method is less reliable
                     * computing correct subspaces, an option is added that
                     * explicitly allows to choose if svds is used rather than
                     * svd. */
                    if this.UseSVDS || issparse(data)
                        [u, s] = svds(data, target_dim);
                        s = diag(s);
                    else
                        [u,s] = svd(data," econ ");
                        s = diag(s);
                        err = sum(s(target_dim+1:end));
                        rerr = 100*err/sum(s);
                        u(:,target_dim+1:end) = [];
                    end
                end
                if KerMor.App.Verbose > 2
                    fprintf(" error %g (%g%% relative)\n ",err,rerr);
                end

                /*  Safety for zero singular values. */
                z = find(s(1:size(u,2)) == 0);
                if ~isempty(z)
                    warning(" KerMor:POD "," %d of %d selected singular values are zero. Assuming output size %d ",...
                        length(z),length(s),length(s)-length(z));
                end
                u(:,z) = [];

                /*  Re-create the sparse matrix structure after POD computation if input argument
                 * was a sparse matrix */
                if ~isempty(nonzero)
                    u = MatUtils.toSparse(u, nonzero, fulldim);
                end

                podvec = u;    
            end
            /*  Security checks (only for true matrices so far) */
            if any(any(round(podvec^t*podvec) ~= eye(size(podvec,2))))
                warning(" KerMor:POD ",[" Resulting POD vectors not  " ...
                    " sufficiently orthogonal. Orthonormalizing. "]);
                /*  Sometimes the u vectors generated by svds are not
                 * properly orthonormal */
                o = general.Orthonormalizer;
                podvec = o.orthonormalize(podvec);
            end
        }
#if 0 //mtoc++: 'set.UseSVDS'
function  UseSVDS(value) {
            if ~islogical(value)
                error(" Value should be logical ");
            end
            this.UseSVDS= value;
        }

#endif


        
#if 0 //mtoc++: 'set.Mode'
function  Mode(value) {
            if ~any(strcmpi(value, [" sign "," eps "," abs "," rel "]))
                error([" Unknown POD reduction mode:  " value " \nAllowed: sign, eps, abs, rel "]);
            end
            this.Mode= lower(value);
        }

#endif


        
#if 0 //mtoc++: 'set.Value'
function  Value(value) {
            if ~isreal(value) || ~isscalar(value)
                error(" Value property must be a positive real scalar. ");
            end
            this.Value= value;
        }

#endif

    
    public: /* ( Static ) */ /* ( setObservable ) */


        static function res = test_POD() {
            pod = general.POD;
            res = true;

            /*  d < N */
            vec = rand(50,1000);
            res = res && general.POD.internalPODTest(pod, vec);
            res = res && general.POD.internalPODTest(pod, vec, 1:30);
            res = res && general.POD.internalPODTest(pod, vec, 30:50);

            /*  d > N */
            vec = rand(1000,50);
            res = res && general.POD.internalPODTest(pod, vec);
            res = res && general.POD.internalPODTest(pod, vec, 1:300);
            res = res && general.POD.internalPODTest(pod, vec, 800:900);
        }


    
    private: /* ( Static ) */ /* ( setObservable ) */

        static function res = internalPODTest(pod,vec,selection) {
            if nargin < 3
                selection = " : ";
            end
            res = true;
            pod.Mode= " eps ";
            pod.Value= 7;
            V = pod.computePOD(vec,[],selection);
            res = res && isequal(round(V^t*V),eye(size(V,2)));

            pod.Mode= " sign ";
            pod.Value= .3;
            V = pod.computePOD(vec,[],selection);
            res = res && isequal(round(V^t*V),eye(size(V,2)));

            pod.Mode= " abs ";
            pod.Value= 10;
            V = pod.computePOD(vec,[],selection);
            res = res && isequal(round(V^t*V),eye(size(V,2)));

            /*  Vexclude mode */
            o = general.Orthonormalizer;

            exclu = o.orthonormalize(rand(length(vec(selection,1)),round(pod.Value/2)));
            Vex = pod.computePOD(vec,exclu,selection);
            res = res && isequal(round(Vex^t*Vex),eye(size(Vex,2)));
            res = res && norm(exclu^t*Vex) < 1e-12;

            pod.Mode= " abs ";
            pod.Value= 10;
            Vf = pod.computePOD(data.FileMatrix(vec),[],selection);
            res = res && isequal(round(Vf" *Vf),eye(size(Vf,2))) && norm(abs(V)-abs(Vf), "fro^t) < 1e-10;

            pod.Mode= " rel ";
            pod.Value= .3;
            V = pod.computePOD(vec,[],selection);
            res = res && isequal(round(V^t*V),eye(size(V,2)));

            pod.Mode= " rel ";
            pod.Value= .3;
            Vf = pod.computePOD(data.FileMatrix(vec),[],selection);
            res = res && isequal(round(Vf" *Vf),eye(size(Vf,2))) && norm(abs(V)-abs(Vf), "fro^t) < 1e-10;
        }

    
};
}