GaussKernel.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 GaussKernel
  :public kernels.BellFunction {
    public:

        GaussKernel(double Gamma) {
            this = this@kernels.BellFunction;
            this.registerProps(" Gamma ");/* ,'Sigma' */


            if nargin == 1
                this.Gamma= Gamma;
            end
            this.updateGammaDependants;
            this.addlistener(" Gamma "," PostSet ",@this.updateGammaDependants);
        }
        function K = evaluate(matrix<double> x,matrix<double> y) {
            K = exp(-this.getSqDiffNorm(x, y)/this.Gamma^2);
        }
        function Nablax = getNabla(colvec<double> x,matrix<double> y) {
            if size(x,2) > 1 && size(y,2) > 1
                error(" One argument must be a vector. ");
            end
            if ~isempty(this.fP)
                error(" Not yet implemented correctly. ");
/*                  xl = x(this.P,:);
 *                 yl = y(this.P,:);
 *             else
 *                 xl = x; yl = y; */
            end
            hlp = bsxfun(@minus,x,y);
            hlp = -2*hlp/this.Gamma^2;
            Nablax = bsxfun(@times,hlp,this.evaluate(x, y));
        }
        function dx = evaluateD1(r) {
            dx = -2*r/this.Gamma^2 .* exp(-r.^2/this.Gamma^2);
        }
        function ddx = evaluateD2(r) {
            ddx = (2/this.Gamma^2) * (2*r.^2/this.Gamma^2-1) .* exp(-r.^2/this.Gamma^2);
        }
        function phi = evaluateScalar(r) {
            phi = exp(-r.^2/this.Gamma^2);
        }
        function g = setGammaForDistance(dist,ep) {
            if nargin == 2
                ep = eps;
            end
            g = dist/sqrt(-log(ep));
            this.Gamma= g;

            if KerMor.App.Verbose > 3
                fprintf(" Setting Gamma = %12.20f\n ",g);
            end
        }
        function copy = clone() {
            copy = kernels.GaussKernel;
            copy.Gamma= this.Gamma;
            copy = clone@kernels.BellFunction(this, copy);
        }

    
    private:

        function  updateGammaDependants(unused1,unused2) {
            this.r0= this.Gamma/sqrt(2);
        }
    protected: /* ( Static ) */

        static function obj = loadobj(obj) {
            if ~isa(obj, " kernels.GaussKernel ")
                newinst = kernels.GaussKernel;
                newinst.Gamma= obj.Gamma;
                newinst.updateGammaDependants;
                obj = loadobj@kernels.BellFunction(newinst, obj);
            else
                obj = loadobj@kernels.BellFunction(obj);
            end
            obj.addlistener(" Gamma "," PostSet ",@this.updateGammaDependants);
        }

        
    public: /* ( Static ) */

        static function res = test_InterpolGamma() {
            ki = general.interpolation.KernelInterpol;
            ki.UseNewtonBasis= false;
            kexp = kernels.KernelExpansion;
            k = kernels.GaussKernel;
            kexp.Kernel= k;
            dx = .2;
            x = -3:dx:3;
            fx = sin(x*pi);
            plot(x,fx);
            epsteps = 0.05:.05:.95;
            dlog = zeros(3,length(epsteps));
            for epidx=1:length(epsteps)
                ep = epsteps(epidx);
                k.setGammaForDistance(dx,ep);
                for idx = 1:length(x)
                    x2 = x;
                    x2(idx) = [];
                    fx2 = fx;
                    fx2(idx) = [];

                    kexp.Centers.xi= x2;
                    ki.init(kexp);
                    kexp.Ma= ki.interpolate(fx2);

                    fxi = kexp.evaluate(x);
                    diff(idx) = abs(fx(idx) - fxi(idx));/* #ok */

                end
                dlog(1,epidx) = ep;
                dlog(2,epidx) = min(diff);
                dlog(3,epidx) = max(diff);
            end            
            disp(dlog);
            [~, idx] = min(dlog(2,:));
            fprintf(" Min distance: %f at ep=%f\n ",dlog(2,idx),dlog(1,idx));
            res = true;
        }


        static function [res , pm ] = test_GaussKernel(pm) {
            if nargin < 1
                pm = PlotManager(false,2,3);
                pm.LeaveOpen= true;
            end
            c = 0;
            x = (-1.2:.01:1.2)+c;
            [X,Y] = meshgrid(x);
            x2 = [X(:)" ; Y(:) "];
            k = kernels.GaussKernel;
            kexp = kernels.KernelExpansion;
            kexp.Kernel= k;
            kexp.Centers.xi= c;
            kexp.Ma= 1;
            conf = [1 2 3];
            for n = 1:length(conf)
                d = conf(n);
                g = k.setGammaForDistance(d,eps);
                tag = strrep(sprintf(" g_1d_d%d_g%g ",d,g)," . "," _ ");
                h = pm.nextPlot(tag,sprintf(" Gauss kernel with dist=%d,\\gamma=%g on 1D data ",d,g));
                plot(h,x,kexp.evaluate(x));
            end
            kexp.Centers.xi= [c; c];
            for n = 1:length(conf)
                d = conf(n);
                g = k.setGammaForDistance(d,eps);
                tag = strrep(sprintf(" g_2d_d%d_g%g ",d,g)," . "," _ ");
                h = pm.nextPlot(tag,sprintf(" Gauss kernel with dist=%d,\\gamma=%g on 2D data ",d,g));
                surf(h,X,Y,reshape(kexp.evaluate(x2),length(x),[])," EdgeColor "," none ");
            end
            if nargin < 1
                pm.done;
            end
            res = true;
        }

    

};
}