KernelInterpol.m

Go to the documentation of this file.

namespace general{
namespace interpolation{


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

        logical UseNewtonBasis = true;
        double RelTol = 1e-13;
    private:

        data.FileMatrix K;
    private:

        kexp;
    public:


        KernelInterpol() {
            this = this@KerMorObject;
            this.MultiTargetComputation= true;
            this.registerProps(" UseNewtonBasis ");
        }


        function copy = clone() {
            copy = general.interpolation.KernelInterpol;
            copy = clone@IKernelCoeffComp(this, copy);
            copy.UseNewtonBasis= this.UseNewtonBasis;
            copy.K= this.K;
            copy.kexp= this.kexp;
        }


        function [matrix<double>ai , matrix<double>nbase , usedcenteridx ] = interpolate(fxi) {
            N = size(fxi,2);
            if this.UseNewtonBasis   
                NV = zeros(N,N);
                c = zeros(size(fxi,1),N);
                fxin = Norm.L2(fxi);
                sumNsq = zeros(1,N);
                fresidual = fxi;
                usedcenteridx = zeros(1,N);
                for m = 1:N
                    res = sum(fresidual.^2,1);
                    [~, maxidx] = max(res);
                    col = this.kexp.getKernelMatrixColumn(maxidx);
                    tN = col - NV(:,1:m-1)*NV(maxidx,1:m-1)^t;
                    tNnorm = sqrt(col(maxidx)-sumNsq(maxidx));
                    if max(res./fxin) < this.RelTol || ~isreal(tNnorm) || tNnorm <= 0
                        m = m-1;/* #ok */

                        break;
                    end
                    NV(:,m) = tN/tNnorm;
                    c(:,m) = fresidual(:,maxidx)./tNnorm;                    
                    fresidual = fresidual - c(:,m)*(NV(:,m))^t;
                    sumNsq = sumNsq + (NV(:,m).^2)^t;
                    usedcenteridx(m) = maxidx;
                end
                usedcenteridx = usedcenteridx(1:m);
                nbase = NV(:,1:m);
                ai = c(:,1:m);

                usedcenteridx = sort(usedcenteridx," ascend ");
                nbase = nbase(usedcenteridx,:);
                if length(usedcenteridx) < N && nargout < 3
                    error(" Less basis translates required than given, but their indices are not considered as output. Please receive the 'used' output and select the indicated subset as centers. ");
                end
            else
                nbase = 1;
                ai = (this.K\fxi" ) ";
                usedcenteridx = 1:N;
            end
        }
        function  init(kernels.KernelExpansion kexp) {
            this.K= [];
            this.kexp= [];
            if this.UseNewtonBasis
                this.kexp= kexp;
            else
                this.K= kexp.getKernelMatrix;
            end
        }
        function [rowvecci , integersvidx , sf ] = computeKernelCoefficients(rowvec<double> fxi,unused1) {

            [ci, base, svidx] = this.interpolate(fxi);
            /*  Here: Transform the coefficients back to direct translate
             * basis as the Base cannot be set here and might be different
             * for each component function */
            if this.UseNewtonBasis
                ci = ci / base;
            end
            sf = StopFlag.SUCCESS;
        }
    public: /* ( Static ) */


        static function  test_KernelInterpolation() {

            pm = PlotManager(false,2,2);
            pm.LeaveOpen= true;

/*              x = -5:.05:5;
 *             fx = sinc(x); */

            x = 0:.05:10;
            fx = x.^2.*sin(2*x);

            n = length(x);
            samp = 1:15:n;

            ke = kernels.KernelExpansion;

            ke.Kernel= kernels.GaussKernel(.5);
            internal_test(x,fx,false);
            internal_test(x,fx,true);

            internal_test(x,ones(size(x))*5,false);
            internal_test(x,ones(size(x))*5,true);

            ke.Kernel= kernels.InvMultiquadrics(-1.4,2);
            internal_test(x,fx,true);
            internal_test(x,fx,false);

            ke.Kernel= kernels.InvMultiquadrics(-4,5);
            internal_test(x,fx,true);
            internal_test(x,fx,false);

            /*  2D test */
            fx = [fx; x.^3.*sin(.4*x)];
            internal_test(x,fx,true);
            internal_test(x,fx,false);

            pm.done;

            function internal_test(x,fx,lu)
                xi = x(samp);
                fxi = fx(:,samp);

                ki = general.interpolation.KernelInterpol;

                ke.Centers.xi= xi;
                ki.UseNewtonBasis= lu;
                ki.init(ke);

                [ke.Ma, ke.Base] = ki.interpolate(fxi);

                fsvr = ke.evaluate(x);

                h = pm.nextPlot(" xx ",sprintf(" Interpolation test with kernel %s, newton=%d ",class(ke.Kernel),lu));
                hold(h," on ");
                for k=1:size(fx,1)
                    plot(h,x,fx(k,:)," r ");
                    /*  Plot approximated function */
                    plot(h,x,fsvr(k,:)," b-- ");
                    /* skipped = setdiff(1:length(x),svidx); */
                    plot(h,xi,fxi(k,:)," .r ");
                end
            end
        }
        static function  test_KernelInterpolationError() {
            realprec = 20;
            range = [-4 4];

            pm = PlotManager(false,2,3);
            pm.LeaveOpen= true;

            ke = kernels.KernelExpansion;
            ke.Kernel= kernels.GaussKernel(1);
            testfun = @(x)sinc(x) * 2 .*sin(3*x);
            ke.Kernel= kernels.InvMultiquadrics(-1,1);

            ki = general.interpolation.KernelInterpol;

            hsteps = [.5:.25:2 2:2:10];
            err = zeros(size(hsteps));
            for n = 1:length(hsteps)
                h = 1/hsteps(n);
                xi = range(1):h:range(2);
                xi = xi + rand(size(xi))*max(xi)/100+.1;
                fxi = testfun(xi);

                x = range(1):h/realprec:range(2);
                fx = testfun(x);

                ke.Centers.xi= xi;
                ki.init(ke);
                [ke.Ma, ke.Base, used] = ki.interpolate(fxi);
                ke.Centers.xi= xi(:,used);
                fint = ke.evaluate(x);

                ax = pm.nextPlot([" step " n],sprintf(" h=%g ",h)," x "," fx ");
                plot(ax,xi,fxi," r. ",x,fx," r ",x,fint," b ");
                err(n) = max(abs(fx-fint));
            end
            pm = PlotManager; pm.LeaveOpen= true;
            ax = pm.nextPlot(" errors "," Error development over h sizes "," 1/h "," errors ");
            h = 1./hsteps;
            semilogy(ax,hsteps,exp(log(h)./h)," b ",hsteps,exp(-1./h)," b-- ",hsteps,err," r ");
            legend(ax," exp(log(h)/h) "," exp(-1/h) "," ||f-sfx||_\infty ");
            pm.done;
        }


};
}
}