VKOGA.m

Go to the documentation of this file.

namespace demos{


/* (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 VKOGA {
    public: /* ( Static ) */


        static function structres = VKOGA_1D_nD(integer n,logical fPGreedy,integer nG) {

            /* % Data setup */
            x = -5:.1:5;
            fx = [sin(x)+.2*x; cos(x)-.1*x; exp(-x.^2).*sin(x); cos(2*x).*sin(x)];
            if nargin < 3
                nG = 3;
                if nargin < 2
                    fPGreedy = false;
                    if nargin < 1
                         n = 1;
                    elseif n > 2
                        n = size(fx,1);
                    end
                end
            end
            fx = fx(1:n,:);
            atd = data.ApproxTrainData(x,[],[]);
            atd.fxi= fx;

            /* % Algorithm setup */
            alg = approx.algorithms.VKOGA;
            alg.MaxExpansionSize= 50;
            alg.MaxRelErr= 1e-5;
            alg.UsefPGreedy= fPGreedy;
            ec = kernels.config.ExpansionConfig;            
            ec.Prototype.Kernel= kernels.GaussKernel(.8);
            ec.StateConfig= kernels.config.GaussConfig(" G ",linspace(.5,2,nG));
            alg.ExpConfig= ec;

            kexp = alg.computeApproximation(atd);

            /* % Plot approximated function */
            m = length(alg.Used);
            pm = PlotManager(false,2,2);
            pm.LeaveOpen= true;

            h1 = pm.nextPlot(" fun "," Function "," x "," f(x) ");
            plot(h1,x,[fx; kexp.evaluate(x)]^t);

            h2 = pm.nextPlot(" nfun "," Newton Basis Functions on training data "," x "," N(x) ");
            plot(h2,x,alg.bestNewtonBasisValuesOnATD); 

            alg.plotSummary(pm, " VKOGA Demo ");
            pm.done;

            alg.getApproximationSummary;

            res.kexp= kexp;
            res.atd= atd;
            res.alg= alg;
        }
        static function  IterationPlots(struct res,integer steps,PlotManager pm) {
            if nargin < 3
                pm = PlotManager(false,3,3);
                pm.LeaveOpen= true;
                if nargin < 2
                    steps = 9;
                end
            end

            if res.alg.UsefPGreedy
                fprintf(2," Warning: f/P-Greedy selection criteria was used. Error plots are for f-Greedy case.\n ");
            end
            ms = 7;
            fx = res.atd.fxi.toMemoryMatrix;
            x = res.atd.xi.toMemoryMatrix;

            /* % Zero function plot */
            k = [];
            h0 = doPlot(0,zeros(size(fx)));

            /* % Iteration plots */
            for s = 1:steps
                k = res.kexp.getSubExpansion(s);
                fxi = k.evaluate(x);
                hn = doPlot(s,fxi);
                axis(hn,axis(h0));
            end

            if nargin < 3
                pm.done;
            end

            function h1 = doPlot(s,fxi)
                h1 = pm.nextPlot(sprintf(" step%d ",s),...
                    sprintf(" Iteration %d ",s)," x "," f(x) ");
                plot(h1,x,fx" , "LineWidth^t,2); 
                hold(h1," on ");

                err = fx-fxi;
                allerr = Norm.L2(err);
                plot(h1,x,allerr," r "); 
                plot(h1,x,fxi" , "--^t); 
                if ~isempty(k)
                    plot(h1,k.Centers.xi,k.evaluate(k.Centers.xi)" , "k." , "MarkerSize^t,17);
                end

                /*  Plot (next) max errors */
                [v, idx] = max(allerr);
                plot(h1,x(idx),v," ro "," MarkerSize ",ms); 
                [~, idx] = max(abs(err),[],2);
                if s == 2
                    for l=1:length(idx)
                        plot(h1,[x(idx(l)) x(idx(l))+eps],[fx(l,idx(l)) fxi(l,idx(l))]," k-- ");
                        plot(h1,[x(idx(l)) x(idx(l))+eps],[fx(l,idx(l)) fxi(l,idx(l))]," kx "," MarkerSize ",ms);
                    end
                end
                axis(h1," tight ");
                axis(h1,axis(h1)*1.03);
            end
        }
        static function  NewtonBasis_Schaback() {

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

            /*  get cut circle */
            [X, ind, Xemesh, Yemesh] = cutcircle(.05);

            /*  assemble training data */
            atd = data.ApproxTrainData(X,[],[]);
            fxi = exp(abs(X(1,:)-X(2,:)))-1;
            atd.fxi= fxi;

            /* % Algorithm setup */
            alg = approx.algorithms.VKOGA;
            /*  Maximal number of data sites to be finally used */
            alg.MaxExpansionSize= 200;
            /*  tolerance for residual of function recovery */
            alg.MaxAbsResidualErr= 1e-6;
            /*  "disable" maxrelerr check by ridiculous small amount */
            alg.MaxRelErr= 1e-10;
            alg.UsefPGreedy= false;
            ec = kernels.config.ExpansionConfig;

            /*  Gaussian
 *             ec.StateConfig = kernels.config.RBFConfig('G',[1 2 3]);  */

            /*  Wendland functions */
            wc = kernels.config.WendlandConfig(" G ",[1 2 3]," S ",[1 1 1]);
            wc.Dimension= 2;
            ec.StateConfig= wc;

            alg.ExpConfig= ec;

            kexp = alg.computeApproximation(atd);

            h = pm.nextPlot(" circ "," Selected data points ");
            plot(h,X(1,:),X(2,:)," . ",X(1,alg.Used),X(2,alg.Used)," ro ");

            Z = zeros(size(Xemesh));
            Z(ind) = fxi;
            h = pm.nextPlot(" fun "," Original function ");
            surf(h,Xemesh,Yemesh,Z);

            h = pm.nextPlot(" afun "," Approximation ");
            aZ = Z;
            aZ(ind) = kexp.evaluate(X);
            surf(h,Xemesh,Yemesh,aZ);

            h = pm.nextPlot(" err "," Error ");
            surf(h,Xemesh,Yemesh,abs(Z-aZ));

            for n=1:8
                plotBasis(n);
            end
            pm.done;

            function plotBasis(idx)
                h = pm.nextPlot(sprintf(" bfun%d ",idx),sprintf(" Newton basis function %d ",idx)," x "," y ");
                nv = zeros(size(Xemesh));
                nv(ind) = alg.bestNewtonBasisValuesOnATD(:,idx);
                surf(h,Xemesh,Yemesh,nv);
            end

            function [Xe, ind, xm, ym]=cutcircle(he)
                /*  generates regular points Xe in circle with lower left quadrangle cut out
                 * The index list ind indicates the positions in meshgrid (xm,ym) 
                 * when the meshgrid is written in a single point list (xee, yee) */
                dx=2; /*  space dimension */

                dy=dx;
                /*  he=0.01; % 1 D stepsize for evaluation grid. 
                 * The grid can consist of very many points.
                 * he=0.01 gives N=40401 points to work on */
                [xm, ym]=meshgrid(-1:he:1,-1:he:1);
                xee=xm(:);
                yee=ym(:);
                ind=find((xee.^2+yee.^2<=1)&((xee>=0)|(yee>=0)) );
/*                  msk=zeros(size(xee));
 *                 msk(ind)=1.0; */
                Xe=[xee(ind), yee(ind)]^t;
            end
        }
};
}