TPWLApprox.m

Go to the documentation of this file.

namespace approx{


/* (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 TPWLApprox
  :public approx.BaseApprox {
    public: /* ( setObservable ) */

        MinWeightValue = 1e-10;
    public: /* ( Dependent, setObservable ) */

        Beta;
    public:

        xi;
        Ai;
        bi;

        GaussWeight;

        fBeta = 25;

        
    public:

        
#if 0 //mtoc++: 'set.MinWeightValue'
function  MinWeightValue(value) {
            if ~isposrealscalar(value)
                error(" MinWeightValue must be a positive real scalar ");
            end
            this.MinWeightValue= value;
        }

#endif


        
#if 0 //mtoc++: 'set.Beta'
function  Beta(value) {
            if ~isposrealscalar(value)
                error(" Beta must be a positive real scalar. ");
            end
            this.fBeta= value;
            this.GaussWeight.Gamma= 1/sqrt(value);
        }

#endif


        
#if 0 //mtoc++: 'get.Beta'
function value = Beta() {
            value = this.fBeta;
        }

#endif


        TPWLApprox(sys) {
            this = this@approx.BaseApprox(sys);

            /*  Beta from Paper is e^{-25 d/m}, for gauss: e^{d/(m*g^2)}, so g=.2
             * Initialize the gaussian weight with what ever Beta is set to */
            this.GaussWeight= kernels.GaussKernel(1/sqrt(this.Beta));
            this.registerProps(" Beta "," MinWeightValue ");

            /*  Set training data selection to epsselector (default strategy for original TPWL) */
            this.TrainDataSelector= data.selection.EpsSelector;
            this.CustomProjection= true;
        }


        function copy = clone() {
            copy = approx.TPWLApprox(this.System);
            copy = clone@approx.BaseApprox(this, copy);
            copy.Ai= this.Ai;
            copy.bi= this.bi;
            copy.xi= this.xi;
            copy.Beta= this.Beta;
            copy.MinWeightValue= this.MinWeightValue;
        }


        function projected = project(V,W) {
            projected = project@approx.BaseApprox(this, V, W, this.clone);
            for i=1:length(this.Ai)
                projected.Ai[i] = W^t*this.Ai[i]*V;
            end
            projected.bi= W^t*this.bi;
            projected.xi= W^t*this.xi;
        }
        function matrix<double>y = evaluate(colvec<double> x,double t,colvec<double> mu) {
            as = length(this.Ai);
            y = zeros(size(x));
            /*  Single argument evaluation */
            if (size(x,2) == 1)
                /*  Compute all distances */
                di = sqrt(sum((this.xi - repmat(x,1,as)).^2));
                di(di == 0) = eps;
                w = this.GaussWeight.evaluateScalar(sqrt(di/min(di)));
                /*  Normalize local weights */
                w = w / sum(w);
                idx = 1:as;
                idx = idx(w > this.MinWeightValue);
                for i=idx
                    y = y + w(i)*(this.Ai[i]*x + this.bi(:,i));
                end
            /*  Multi argument evaluation */
            else
                n = size(x,2);
                xisel = reshape(meshgrid(1:n,1:as),[],1)^t;
                /*  Compute all distances */
                di = sqrt(sum((repmat(this.xi,1,n)-x(:,xisel)).^2));
                di = reshape(di" ,as,[]) ";
                di(di == 0) = eps;
                dimin = min(di,[],2);
                w = this.GaussWeight.evaluateScalar(sqrt(di./repmat(dimin,1,as)));
                w = w ./ repmat(sum(w,2),1,as);
                y = zeros(size(x));
                for i=1:as
                    /*  Find relevant weights */
                    rel = w(:,i) > this.MinWeightValue;
                    /*  Normalize local weights
                     *wl = w(rel,i) ./ s(rel);
                     * Only update entries for relevant weights */
                    hlp = this.Ai[i]*x(:,rel) + this.bi(:,i*ones(1,sum(rel)));
                    y(:,rel) = y(:,rel) + hlp*diag(w(rel,i));
                end
            end
        }
        function  approximateSystemFunction(models.BaseFullModel model) {
            atd = model.Data.ApproxTrainData;
            this.xi= atd.xi;
            ti = atd.ti;
            mui = atd.mui;
            fxi = atd.fxi;

            as = size(this.xi,2);
            this.Ai= cell(0,as);
            for i = 1:as
                if isempty(mui)
                    mu = [];
                else
                    mu = mui(:,i);
                end
                model.System.setConfig(mu,model.DefaultInput);
                this.Ai[i] = model.System.f.getStateJacobian(this.xi(:,i),ti(i));
                this.bi(:,i) = fxi(:,i) - this.Ai[i]*this.xi(:,i);
            end
        }


        function matrix<double>y = evaluateCoreFun() {
        }
    public: /* ( Static ) */

        static function res = test_TWPLApprox() {
            m = models.synth.KernelTest(100,false);
            m.ErrorEstimator= [];
            m.Approx= approx.TPWLApprox(m.System);
            m.Approx.Beta= 25;
            m.dt= 0.01;
            m.Sampler= sampling.ManualSampler(m.getRandomParam);

            m.offlineGenerations;
            mu = m.getRandomParam;
            [~,y] = m.simulate(mu);
            r = m.buildReducedModel;
            [~,yr] = r.simulate(mu);
            fprintf(" Max absolute error: %g\n ",max(Norm.L2(y-yr)));
            res = true;
        }

    
};
}