Orthonormalizer.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 Orthonormalizer
  :public KerMorObject {
    public:

        G = 1;
        Epsilon = 1e-7;
        Algorithm = "gs";
    public:

        function onvec = orthonormalize(vec) {

            /*  Check for nonempty vec */
            if isempty(vec)
                onvec = zeros(size(vec));
                return;
            end

            /*  Check on identity of vectors */
            for i=1:(size(vec,2)-1)
                for j=(i+1):size(vec,2)
                    if isequal(vec(:,i),vec(:,j))
                        vec(:,j) = 0;
                    end;
                end;
            end

            if strcmp(this.Algorithm," gs ")
                onvec = this.ortho_gs(vec);
            elseif strcmp(this.Algorithm," qr ")
                onvec = this.ortho_qr(vec);
            elseif strcmp(this.Algorithm," ch ")
                onvec = this.ortho_ch(vec);
            end
        }
#if 0 //mtoc++: 'set.G'
function  G(value) {
            if ~isa(value, " double ")
                error(" G should be a double matirx ");
            end
            this.G= value;
        }

#endif


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

#endif


        
#if 0 //mtoc++: 'set.Algorithm'
function  Algorithm(value) {
            if ~any(strcmp(value,[" gs " " qr " " ch "]))
                error(" Invalid method: %s ",value);
            end
            this.Algorithm= value;
        }

#endif

    
    private:

        function onvec = ortho_gs(vec) {

            onvec = vec;
            for i = 1:size(vec,2)
                /*  orthogonalize next vector i and assume, that it is already
                 * orthogonal to previous ones */
                n = sqrt(onvec(:,i)^t * this.G * onvec(:,i));
                if (n < this.Epsilon)
                    onvec(:,i) = 0;
                else
                    onvec(:,i) = onvec(:,i)/n;
                end

                /*  orthogonalize remaining vectors wrt this one: */
                A_mult_onvec_i = this.G*onvec(:,i);

                /*  create row-vector with projections on the created on-vector */
                coeffs= A_mult_onvec_i^t * onvec(:,i+1:end);

                /*  create matrix of scaled versions of actual orthonormalized vector */
                coeffmat = onvec(:,i) * coeffs;

                /*  perform orthogonalization */
                onvec(:,i+1:end) = onvec(:,i+1:end) - coeffmat;
            end

            /*  eliminate zero-columns */
            nsqr = sum(onvec.^2);
            onvec = onvec(:,nsqr > 0.1);
        }
        function onvec = ortho_qr(vec) {
            if 1==1 
                error(" Algorithm 'qr' not working. Need to check. "); 
            end
            onvec = vec;

            /*  incomplete cholesky of inner-product matrix: */
            R_M = ichol(sparse(this.G)," inf ");

            /*  qr decomposition of R_M * X with permutation indices E */
            [Q, R, E]= qr(R_M * onvec, 0);

            /*  sort such that linear independent first and Q * R = R_M * Xon */
            onvec = onvec(:,E);

            /*  search nonvanishing diagonal entries of R */
            ind = find(abs(diag(R)) > this.Epsilon);
            onvec = onvec(:,ind);
            Rind = R(ind,ind);
            onvec = onvec(:,ind) / Rind;
        }


        function onvec = ortho_ch(vec) {
            if 1==1 
                error(" Algorithm 'ch' not working. Need to check. "); 
            end
            onvec = vec;

            /*  get gram matrix */
            Gr = onvec^t * this.G * onvec;
            Gr = 0.5* (Gr + Gr^t);

            options.droptol= this.Epsilon;

            /* R = full(cholinc(sparse(G),'inf')); % fill small diagonal entries with inf
             *R = full(cholinc(sparse(G),options)); % fill small diagonal entries with inf */

            /*  cholesky with dropvalue epsilon
             * If an error gets thrown here, you might have an older matlab
             * version only knowing "cholinc" as command, which has been
             * replaced by "ichol" as of R2013A.
 *             R = full(cholinc(sparse(Gr),this.Epsilon));  */
            R = full(ichol(sparse(Gr),options));

            ind = find(diag(R) > this.Epsilon);
            onvec = onvec(:,ind);
            Ri = R(ind,ind);
            onvec = onvec / Ri;

            /*  eliminate too small or too large vectors: */
            norms = onvec^t * this.G * onvec;
            onvec = onvec(:,abs(diag(norms)-1)<0.5);
        }

    
    public: /* ( Static ) */

        static function res = test_Orthogonalization() {
            res = true;
            for k = 1:5
                d2 = randi(200);
                d1 = d2 + k*randi(10);

                A = rand(d1,d2);

                o = general.Orthonormalizer;
                o.Algorithm= " gs ";
                ov1 = o.orthonormalize(A);
    /*              o.Algorithm = 'ch';
     *             ov2 = o.orthonormalize(A);
     *             o.Algorithm = 'qr';
     *             ov3 = o.orthonormalize(A) */

    /*              res = norm(ov1-ov2) < sqrt(eps); */
                res = res & isequal(round(ov1^t*ov1),eye(d2));
                res = res & abs(norm(ov1^t*ov1)-1) < 1000*eps;
            end
        }

    
};
}