ResponseSurfaceApprox.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 ResponseSurfaceApprox
  :public handle {
public:

    EC;

    XYRange = -5:.1:5;

    PM;

    TrainPerc = .2;

    NumSteps = 10;

    LoopCallback;


public:

    X;

    Y;

    Z;

    atd;

    
public:

    ResponseSurfaceApprox() {
        ec = kernels.config.ExpansionConfig;
        ec.Prototype.Kernel= kernels.Wendland;
        ec.StateConfig= kernels.config.WendlandConfig(" G ",10," S ",2," Dim ",2);
        this.EC= ec;
        pm = PlotManager(false,3,2);
        pm.LeaveOpen= true;
        this.PM= pm;            
    }


    function Z = getSurface(X,Y) {
        Z = sin(X/2).*Y.^2+5*X+Y;
        this.Z= Z;
    }


    function alg = setupVKOGA() {
        alg = approx.algorithms.VKOGA;
        alg.MaxExpansionSize= this.NumSteps;
        alg.UsefPGreedy= 0;
        alg.ExpConfig= this.EC;
    }


    function kernels.KernelExpansionkexp = compute() {
        this.init;
        alg = this.setupVKOGA;
        kexp = alg.computeApproximation(this.atd);
    }


    function  init(sel) {
        [X,Y] = meshgrid(this.XYRange);
        if nargin < 2
            tot = numel(X);
            sel = randperm(tot);
            n = round(this.TrainPerc*tot);
        end
        Z = this.getSurface(X,Y);
        this.Z= Z;
        fxi = Z(sel(1:n));
        xi = [X(sel(1:n)); Y(sel(1:n))];
        atd = data.ApproxTrainData(xi,[],[]);
        atd.fxi= fxi;
        this.X= X;
        this.Y= Y;
        this.atd= atd;
    }


    function  plotOriginal() {
        h = this.PM.nextPlot(" fun "," Original response surface "," x_1 "," x_2 ");
        surfl(h, this.X, this.Y, this.Z);
        zlabel(" f(x) ");
        shading interp;
    }


    function  plotOriginalWithTD() {
        h = this.PM.nextPlot(" fun_atd "," Original response surface with training data "," x "," f(x) ");
        surf(h, this.X, this.Y, this.Z);
        shading interp;
        hold(h," on ");
        xi = this.atd.xi.toMemoryMatrix;
        plot3(h,xi(1,:),xi(2,:),fxi," r. ");
    }


    function  plotTrainData() {
        h = this.PM.nextPlot(" traindata "," Training data "," x_1 "," x_2 ");
        xi = this.atd.xi.toMemoryMatrix;
        plot3(h,xi(1,:),xi(2,:),fxi," r. ");
    }


    function  iterationPlots(kernels.KernelExpansion kexp,pm) {
        if nargin < 3
            pm = this.PM;
        end
        allxi = [this.X(:) this.Y(:)]^t;

        /*          aZ = kexp.evaluate(allxi);
         *         aZ = reshape(aZ,size(X,1),[]);
         *         h = pm.nextPlot('approx','Approx','x','f(x)');
         *         surf(h, X, Y, aZ);
         *         h = pm.nextPlot('error','Error','x','f(x)');
         *         surf(h, X, Y, abs(Z-aZ)); */
        args = [" EdgeColor "," interp "," FaceAlpha ",.2," EdgeColor "," none "," FaceColor "," red "];

        centerpos = Utils.findVecInMatrix(allxi,kexp.Centers.xi);

        /* % Iteration plots */
        zlim = [];
        for s = 0:this.NumSteps
            if s == 0
                fxi = zeros(size(this.X));
            else
                k = kexp.getSubExpansion(s);
                fxiflat = k.evaluate(allxi);
                fxi = reshape(fxiflat,size(this.X,1),[]);
            end
            if isa(pm," PlotManager ")
                h = pm.nextPlot(sprintf(" iter%d ",s),...
                    sprintf(" Approximation at step %d ",s)," x_1 "," x_2 ");
            else
                h = pm;
                cla(h);
            end
            surf(h, this.X, this.Y, fxi," EdgeColor "," interp ");
            zlabel(" f(x) ");
            hold(h," on ");
            surf(h, this.X, this.Y, this.Z,args[:]);
            if s > 0
                plot3(allxi(1,centerpos(s)),allxi(2,centerpos(s)),...
                fxiflat(1,centerpos(s))," y* "," MarkerSize ",20);
                plot3(allxi(1,centerpos(s)),allxi(2,centerpos(s)),...
                    fxiflat(1,centerpos(s))," k. "," MarkerSize ",20);
                if s > 1
                    plot3(allxi(1,centerpos(1:s-1)),allxi(2,centerpos(1:s-1)),...
                        fxiflat(1,centerpos(1:s-1))," r. "," MarkerSize ",16);
                end
            end
            /*      h = pm.nextPlot(sprintf('err_iter%d',s),...
             *         sprintf('Error at step %d',s),'x','f(x)');
             *     surf(h, X, Y, abs(Z-fxi),'EdgeColor','interp');
             *axis(h,'tight'); */
            if ~isempty(this.LoopCallback)
                this.LoopCallback(this,s,h);
            end
            if isempty(zlim)
                zlim = get(h," ZLim ");
            else
                set(h," ZLim ",zlim);
            end
        end

        if nargin < 2
            pm.done;
        end
    }



public: /* ( Static ) */

    static function  run() {

        rsa = demos.ResponseSurfaceApprox;
        kexp = rsa.compute;
        rsa.iterationPlots(kexp);
    }
};
}