KernelExpansion.m

Go to the documentation of this file.

namespace kernels{


/* (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 KernelExpansion
  :public KerMorObject,
   public general.AProjectable,
   public dscomponents.IGlobalLipschitz {
    public: /* ( Dependent, setObservable ) */

        kernels.BaseKernel Kernel;
    public: /* ( Dependent ) */

        rowvec Ma_norms;
        double NativeNorm;
        colvec<double> ComponentNorms;
        MBnd;
        logical HasCustomBase;
        Size;
    public: /* ( setObservable ) */

        struct Centers = struct("'xi',[]");
        matrix<double> Ma = "[]";
        matrix<double> Base = 1;
    protected:

        kernels.BaseKernel fSK;
    public:

        KernelExpansion() {

            /*  Set custom projection to true as the project method is
             * overridden */
            this = this@KerMorObject;

            this.fSK= kernels.GaussKernel;

            /*  DONT CHANGE THIS LINE unless you know what you are doing or
             * you are me :-)
 *             this.updateRotInv; */

            this.registerProps(" Kernel "," Centers "," Base ");
        }
        function fx = evaluate(matrix x,varargin) {
            fx = this.Ma * (this.Base \ this.getKernelVector(x)^t);
        }
        function phi = getKernelVector(matrix x,varargin) {
            phi = this.fSK.evaluate(x, this.Centers.xi);
        }
        function row = getKernelMatrixColumn(integer idx,matrix x,varargin) {
            if nargin < 3
                x = this.Centers.xi;
            end
            row = this.fSK.evaluate(x, x(:,idx));
        }
        function J = getStateJacobian(x,varargin) {
            J = this.Ma * (this.Base \ this.fSK.getNabla(x, this.Centers.xi)^t);
        }
        function K = getKernelMatrix() {
            K = [];
            if ~isempty(this.Centers)
                K = this.fSK.evaluate(this.Centers.xi,[]);
/*                  cent = this.Centers.xi;
 *                 if isa(cent,'data.FileMatrix')
 *                     cent = cent.toMemoryMatrix;
 *                 end
 *                 K = this.fSK.evaluate(cent,[]); */
            end
        }
        function v = scalarProductWith(f) {
            if ~isa(f," kernels.KernelExpansion ")
                error(" f argument must be another KernelExpansion ");
            end
            if size(this.Centers.xi,1) ~= size(f.Centers.xi,1)
                error(" Function input dimension mismatch ");
            end
            if size(this.Ma,1) ~= size(f.Ma,1)
                error(" Function output dimension mismatch ");
            end
            v = sum(sum(this.Ma .* f.evaluate(this.Centers.xi),2));
        }


        function  normalize() {
            this.Ma= this.Ma / this.NativeNorm;
        }


        function  setCentersFromATD(data.ApproxTrainData atd,rowvec<integer> idx) {
            this.Centers.xi= atd.xi(:,idx);
        }
        function kernels.KernelExpansionkexp = toTranslateBase() {
            kexp = this.clone;
            if kexp.HasCustomBase
                kexp.Ma= kexp.Ma/kexp.Base;
                kexp.Base= 1;
            end
        }
        function subexp = getSubExpansion(arg) {
            subexp = this.clone;
            if isscalar(arg)
                sel = 1:arg;
            else
                sel = arg;
            end
            subexp.Ma= subexp.Ma(:,sel);
            subexp.Centers.xi= subexp.Centers.xi(:,sel);
            if isfield(this.Centers," ti ") && ~isempty(this.Centers.ti)
                subexp.Centers.ti= this.Centers.ti(sel);
            end
            if isfield(this.Centers," mui ") && ~isempty(this.Centers.mui)
                subexp.Centers.mui= this.Centers.mui(:,sel);
            end
            if subexp.HasCustomBase
                subexp.Base= subexp.Base(sel,sel);
            end
        }


        function c = getGlobalLipschitz(double t,colvec<double> mu) {
            c = sum(this.Ma_norms) * this.Kernel.getGlobalLipschitz;
        }


        function target = project(V,W) {
            target = this.clone;
            k = target.Kernel;
            /*  For rotation invariant kernel expansions the snapshots can be
             * transferred into the subspace without loss. */
            if k.IsRBF
                target.Centers.xi= W^t * this.Centers.xi;
                PG = V^t*(k.G*V);
                if all(all(abs(PG - eye(size(PG,1))) < 100*eps))
                    PG = 1;
                end
                target.Kernel.G= PG;
            elseif k.IsScProd
                target.Centers.xi= V^t * (k.G * this.Centers.xi);
                target.Kernel.G= 1;
            end
            target.Ma= W^t*this.Ma;
        }


        function copy = clone(copy) {
            if nargin == 1
                copy = kernels.KernelExpansion;
            end
            /*  Copy local variables */
            copy.Centers= this.Centers;
            copy.Ma= this.Ma;
            copy.Base= this.Base;
            copy.fSK= this.fSK.clone;
        }


        function  clear() {
            this.Ma= [];
            this.Centers.xi= [];
        }
        function sum = plus(B) {
            sum = compose(A, B, 1);
        }
        function diff = minus(B) {
            diff = compose(A, B, -1);
        }
        function neg = uminus() {
            neg = this.clone;
            neg.Ma= -this.Ma;
        }


        function  exportToDirectory(dir) {
            if nargin < 2
                dir = pwd;
            end
            Utils.ensureDir(dir);
            k = this.Kernel;
            if isa(k," kernels.GaussKernel ")
                kdata = [1 k.Gamma];
            elseif isa(k," kernels.Wendland ")
                kdata = [2 k.Gamma k.d k.k];
            else
                error(" Cannot export a kernel of type %s yet. ",class(this.Kernel));
            end
            Util.saveRealVector(kdata," kernel.bin ",dir);
            Util.saveRealMatrix(this.Centers.xi," centers.bin ",dir);
            Util.saveRealMatrix(this.Ma," coeffs.bin ",dir);
        }


        function json = toJSON(filename) {
            json = " { ";
            /* sprintf('{%s,"centers":%s}',coef,cent); */
            json = [json " 'coeffs': " Util.matrixToJSON(this.Ma)];
            json = [json " , 'centers': " Util.matrixToJSON(this.Centers.xi)];
            if this.HasCustomBase
                json = [json " , 'base': " Util.matrixToJSON(this.Base)];
            end
            json = [json " } "];
            if nargin > 1
                fh = fopen(filename," w+ ");
                fprintf(fh," %s ",json);
                fclose(fh);
            end
        }


        function res = test_getStateJacobianInstance() {
            n = 5;
            xdim = size(this.Centers.xi,1);
            meanx = mean(this.Centers.xi,2);
            rs = RandStream(" mt19937ar "," Seed ",round(rand*1000));
            allX = bsxfun(@plus,meanx,rs.rand(xdim,n));
            allX = [this.Centers.xi(:,1) allX];
            res = true;
            for k = 1:n
                x = allX(:,k);
                dx = ones(xdim,1)*sqrt(eps(class(x))).*max(abs(x),1);
                X = repmat(x,1,xdim);
                Jc = (this.evaluate(X+diag(dx))-this.evaluate(X))*diag(1./dx);
                J = this.getStateJacobian(x);
                maxrel = max(max(abs((J-Jc)./Jc)));
                if maxrel > 1e-2
                    fprintf(" %s: Max relative error: %g\n ",class(this.Kernel),maxrel);
                    res = false;
                end
            end
        }

    
        /* % Getter & Setter */
    public:

        
#if 0 //mtoc++: 'set.Centers'
function  Centers(value) {
            C = [" xi "," ti "," mui "];            
            if isempty(value) || any(isfield(value, C))
                this.Centers= value;
/*                  if ~isfield(value,'xi')
 *                     warning('REQUIRED_FIELD:Empty','xi is a required field, which is left empty');
 *                 end */
            else
                error(" Value passed is not a valid struct. Only the field names xi, ti, mui are allowed ");
            end            
        }

#endif


        
#if 0 //mtoc++: 'get.Kernel'
function k = Kernel() {
            k = this.fSK;
        }

#endif


        
#if 0 //mtoc++: 'get.Size'
function value = Size() {
            value = 0;
            if ~isempty(this.Centers)
                value = size(this.Centers.xi,2);
            end
        }

#endif


        
#if 0 //mtoc++: 'set.Kernel'
function  Kernel(value) {
            if isa(value," kernels.BaseKernel ")
                this.fSK= value;
/*                  this.updateRotInv; */
            else
                error(" Kernel must be a subclass of kernels.BaseKernel. ");
            end
        }

#endif


        
#if 0 //mtoc++: 'get.Ma_norms'
function m = Ma_norms() {
            m = sqrt(sum(this.Ma.^2,1));
        }

#endif


        
#if 0 //mtoc++: 'get.ComponentNorms'
function n = ComponentNorms() {
            if this.HasCustomBase
                c = this.Ma/this.Base;
            else
                c = this.Ma;
            end
            n = sqrt(sum(c"  .* (this.getKernelMatrix * c "),1))^t;
        }

#endif


        
#if 0 //mtoc++: 'get.NativeNorm'
function n = NativeNorm() {

            /*  doesnt use ComponentNorms as this way we save "sqrt(x).^2" */
            if this.HasCustomBase
                c = this.Ma/this.Base;
            else
                c = this.Ma;
            end
            n = sqrt(sum(sum(c"  .* (this.getKernelMatrix * c "))));
        }

#endif

#if 0 //mtoc++: 'get.MBnd'
function M = MBnd() {
            if this.HasCustomBase
                c = this.Ma/this.Base;
            else
                c = this.Ma;
            end
            M = max(sum(abs(c),2));
        }

#endif


        
#if 0 //mtoc++: 'get.HasCustomBase'
function v = HasCustomBase() {
            v = ~isempty(this.Base) && (numel(this.Base) > 1 || this.Base ~= 1);
        }

#endif

    
    private:

        function c = compose(B,sgn) {
            if ~isa(A," kernels.KernelExpansion ") || ~isa(B," kernels.KernelExpansion ")
                error(" Both arguments must be at least kernels.KernelExpansion subclasses ");
            end
            if A.Kernel ~= B.Kernel
                error(" Cannot add kernel expansions: different state space kernels ");
            elseif ~isempty(A.Ma) && ~isempty(B.Ma) && size(A.Ma,1) ~= size(B.Ma,1)
                error(" Cannot add kernel expansions: different output dimensions ");
            elseif ~isempty(A.Centers.xi) && ~isempty(B.Centers.xi) && size(A.Centers.xi,1) ~= size(B.Centers.xi,1)
                error(" Cannot add kernel expansions: different statespace center dimensions ");
            end
            c = A.clone;
            c.Centers.xi= [A.Centers.xi B.Centers.xi];
            c.Ma= [A.Ma sgn*B.Ma];
        }

    
    public: /* ( Static ) */

        static function res = test_getStateJacobian() {
            rs = RandStream(" mt19937ar "," Seed ",0);
            ke = kernels.KernelExpansion;
            k1 = kernels.GaussKernel(1);
            k2 = kernels.Wendland;

            s = [1 1 2 2 3 3];
            d = [1 20 1 20 1 20];
            gam = 4 + (1:6);
            n = 6;
            nc = round(rs.rand(1,6)*100);

            res = true;
            for k = 1:n
                ke.Centers.xi= rs.rand(d(k),nc(k));
                ke.Ma= rs.rand(d(k)+1,nc(k));

                k1.Gamma= gam(k);
                ke.Kernel= k1;
                res = res & ke.test_getStateJacobianInstance;

                k2.Gamma= gam(k);
                k2.d= d(k);
                k2.k= s(k);
                ke.Kernel= k2;
                res = res & ke.test_getStateJacobianInstance;
            end
        }

    
    protected: /* ( Static ) */

        static function this = loadobj(this,from) {
            if nargin > 1
                this.Kernel= from.Kernel;
                this.Ma= from.Ma;
                this.Centers= from.Centers;
                if isfield(from," Base ")
                    this.Base= from.Base;
                end
                this = loadobj@KerMorObject(this, from);
            elseif ~isa(this, " kernels.KernelExpansion ")
                newinst = kernels.KernelExpansion;
                newinst.Kernel= this.fSK;
                newinst.Ma= this.Ma;
                newinst.Centers= this.Centers;
                if isfield(this," Base ")
                    newinst.Base= this.Base;
                end
                this = loadobj@KerMorObject(newinst, this);
            end
        }
};
}