Wendland.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 Wendland
  :public kernels.ARBFKernel {
    public: /* ( Dependent ) */

        d;
        k;
    private:

        fd = 1;

        fk = 0;

        expo;

        co;

        polystr;

        polyfun;

        polystrD1 = {"   \
                '\@(r)(l+1)'    \
                '\@(r)2*(l^2+4*l+3)*r/3+(l+2)'    \
                '\@(r)(3*(l^3+9*l^2+23*l+15)*r.^2 + 2*(6*l^2+36*l+45)*r)/15+(l+3)'"};

        polyfunD1;

        
    public:

        Wendland() {
            this.polystr= [...
                " @(r)(l+1)*r+1 " ...
                " @(r)(l^2+4*l+3)*r.^2/3+(l+2)*r+1 " ...
                " @(r)((l^3+9*l^2+23*l+15)*r.^3 + (6*l^2+36*l+45)*r.^2)/15+(l+3)*r+1 "];
            this.updateCoeffs;
        }


        function Kxy = evaluate(colvec<double> x,matrix<double> y) {
            r = sqrt(this.getSqDiffNorm(x, y))/this.Gamma;
            rp = max(1-r, 0);
            p = 1;
            if (this.fk > 0)
                p = this.polyfun(r);
            end
            Kxy=(rp.^this.expo).*p;
        }


        function Kxy = evaluateScalar(r) {
            error(" not implemented ");
/*              rp = max(1-r,0);
 *             %rp = (1-r).*(r <= 1);
 *             
 *             p = 1;
 *             if (this.fk > 0)
 *                 p = this.polyfun(r);
 *             end
 *             Kxy=(rp.^this.expo).*p; */
        }



        function Nabla = getNabla(colvec<double> x,matrix<double> y) {
            if ~isempty(this.fP)
                error(" Not yet implemented correctly. ");
            end
            r = sqrt(this.getSqDiffNorm(x, y))/this.Gamma;
            rp = max(1-r, 0);
            p = 1;
            dp = 0;
            if (this.fk > 0)
                p = this.polyfun(r);
                dp = this.polyfunD1(r);
            end
            hlp = (-this.expo*(rp.^(this.expo-1)).*p + (rp.^this.expo).*dp)./(r*this.Gamma^2);
            /*  For exact center matches r=0, so we would have inf/-infs */
            hlp(~isfinite(hlp)) = 0;
            Nabla = bsxfun(@times,hlp,bsxfun(@minus,x,y));
        }
         function dx = evaluateD1(r) {
            error(" not implemented ");
        end

        function ddx = evaluateD2(this, r)
            /*  Method for second derivative evaluation */
            error(" not implemented ");
        end

        /*  Returns the global lipschitz constant of this kernel.
         *
         * Exprimental state as not implemented & checked for all kernels. */
        function c = getGlobalLipschitz(this)
            error(" not implemented ");
        end

        function copy = clone(this)
            copy = clone@kernels.ARBFKernel(this, kernels.Wendland);
            /*  Triggers update of coeffs */
            copy.d= this.fd;
            copy.k= this.fk;
        end
    end

    methods(Access=private)
        function updateCoeffs(this)
            /*  calculates coefficients and exponent for polynomial
             * part p_{d,k}(r) of the Wendland kernel function. */
            l = floor(this.fd/2) + this.fk + 1;
            this.expo= l + this.fk;
            this.polyfun= [];
            if (this.fk > 0)
                this.polyfun= eval(this.polystr[this.fk]);
                this.polyfunD1= eval(this.polystrD1[this.fk]);
            end
        end
    end

    /* % Getter & Setter */
    methods
        function set.k(this, value)
            if value < 0 || value > 3
                error(" Only k values in [0,3] are allowed ");
            elseif round(value) ~= value
                error(" Only the k values 0,1,2,3 are allowed ");
            end
            this.fk= value;
            this.updateCoeffs;
        end

        function set.d(this, value)
            if round(value) ~= value
                error(" The dimension must be an integer ");
            end
            this.fd= value;
            this.updateCoeffs;
        end

        function v = get.d(this)
            v = this.fd;
        end

        function v = get.k(this)
            v = this.fk;
        end
    end

    methods(Static)
        function res = test_WendlandKernel(pm)
            if nargin < 1
                pm = PlotManager(false,3,3);
                pm.LeaveOpen= true;
            end
            c = 0;
            x = (-1.2:.01:1.2)+c;
            [X,Y] = meshgrid(x);
            x2 = [X(:)" ; Y(:) "];
            k = kernels.Wendland;
            kexp = kernels.KernelExpansion;
            kexp.Kernel= k;
            kexp.Centers.xi= c;
            kexp.Ma= 1;
            conf = Utils.createCombinations([1 2 3 4 5],[0 1 2 3]);
            for n = 1:length(conf)
                k.d= conf(1,n);
                k.k= conf(2,n);
                tag = sprintf(" w_1d_d%d_k%d ",k.d,k.k);
                h = pm.nextPlot(tag,sprintf(" Wendland kernel with d=%d,k=%d on 1D data ",k.d,k.k));
                plot(h,x,kexp.evaluate(x));
            end
            kexp.Centers.xi= [c; c];
            for n = 1:length(conf)
                k.d= conf(1,n);
                k.k= conf(2,n);
                tag = sprintf(" w_2d_d%d_k%d ",k.d,k.k);
                h = pm.nextPlot(tag,sprintf(" Wendland kernel with d=%d,k=%d on 2D data ",k.d,k.k));
                surf(h,X,Y,reshape(kexp.evaluate(x2),length(x),[])," EdgeColor "," none ");
            end
            if nargin < 1
                pm.done;
            end
            res = true;
        end
    end

    methods(Static,Access=protected)
        function this = loadobj(this)
            this = loadobj@kernels.ARBFKernel(this);
            this.updateCoeffs;
        end
    end
end
         }
};