software/kermor/_bell_function_8m_source.html

 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 BellFunction

   :public kernels.ARBFKernel {

     public: /* ( Dependent, setObservable ) */


         r0;

     public: /* ( setObservable ) */


         NewtonTolerance = 1e-7;

         MaxNewtonIterations = 5000;

     public:


         rm;

     private: /* ( Transient ) */


         p;


         sc;


         d1;


         d2;


     private:


         fr0;


     public:


         BellFunction() {

             this = this@kernels.ARBFKernel;

             this.registerProps(" r0 "," NewtonTolerance "," MaxNewtonIterations ");

         }


         function c = getGlobalLipschitz() {

             c = abs(this.evaluateD1(this.r0));

         }

         function copy = clone(copy) {

             copy = clone@kernels.ARBFKernel(this, copy);

             copy.NewtonTolerance= this.NewtonTolerance;

             copy.MaxNewtonIterations= this.MaxNewtonIterations;

             copy.fr0= this.fr0;

             copy.p= this.p;

             copy.sc= this.sc;

             copy.d1= this.d1;

             copy.d2= this.d2;

             copy.rm= this.rm;

             copy.r0= this.r0;

         }


         function rtmp = ModifiedNewton(rstart,s) {

             rtmp = rstart+2*this.NewtonTolerance;

             r = rstart;


             /*  Add eps at division as for zero nominator the expression is zero. (no error made) */

             n = zeros(3,length(s));

             n(1,:) = this.d1(1) - (this.sc(1)-this.evaluateScalar(s))./(eps-s);

             n(2,:) = this.d1(2) - (this.sc(2)-this.evaluateScalar(s))./(this.p(2)-s+eps);

             n(3,:) = this.d1(3) - (this.sc(3)-this.evaluateScalar(s))./(this.p(3)-s+eps);


             /*  Lim x->y n'(x) = \phi''(x) / 2 (De L'hospital) */

             dn = zeros(3,length(s));

             dn(1,:) = this.d2(1) - n(1,:)./-s;

             dn(1,isnan(dn(1,:))) = this.evaluateD2(0)/2;

             dn(2,:) = this.d2(2) - n(2,:)./(this.p(2)-s);

             dn(2,isinf(dn(2,:))) = 0; /*  == ddf(x0)/2; */


             dn(3,:) = this.d2(3) - n(3,:)./(this.p(3)-s);

             dn(3,isinf(dn(3,:))) = this.d2(3)/2;


             cnt = 0;


             /* conv = 1:length(x);

              *finished = zeros(size(x));

              *while ~isempty(conv) && cnt < this.MaxNewtonIterations */

             while any(abs(rtmp-r) > this.NewtonTolerance) && cnt < this.MaxNewtonIterations


                 /*  Experimental setting: thinning out the finished values does not improve the

                  * overall performance.

                  * Find indices of finished values

                  *                 fidx = abs(xtmp-x) <= this.NewtonTolerance;

                  *                 if any(fidx)

                  *                     % Store finished values

                  *                     finished(conv(fidx)) = xtmp(fidx);

                  *                     conv(fidx) = [];

                  *                     x(fidx) = []; xtmp(fidx) = []; y(fidx) = [];

                  *                     n0(fidx) = []; dn0(fidx) = [];

                  *                     nx0(fidx) = []; dnx0(fidx) = [];

                  *                     nxr(fidx) = []; dnxr(fidx) = [];

                  *                 end */


                 [g,dg] = this.optFun(r,s,n,dn);


                 rtmp = r;

                 r = r - g./dg;

                 cnt = cnt + 1;


                 /* if cnt > this.MaxNewtonIterations*.5

                  *showNewton;

                  *end */

             end

             if cnt == this.MaxNewtonIterations

                 error(" Bellfunction->ModifiedNewton: Max iterations of %d reached ",this.MaxNewtonIterations);

             end


             /*              function showNewton %#ok<*DEFNU>

              *                 % Fix an identical y for all iterations

              *                 si = 200;

              *                 for i=1:5

              *

              *                     X = linspace(.5*min(min([x; xtmp],[],1)),2*max(max([x; xtmp],[],1)),si);

              *                     gs = zeros(length(y),si);

              *                     for idx = 1:length(y)

              *                         gs(idx,:) = error.lipfun.ImprovedLocalSecantLipschitz.optFun(X,repmat(y(idx),size(X)),f,df,ddf,this.x0);

              *                     end

              *

              *                     [g,dg] = error.lipfun.ImprovedLocalSecantLipschitz.optFun(x,y,f,df,ddf,this.x0);

              *                     xtmp = x;

              *                     x = x - g./dg;

              *                     %plot(xs,gs,'r');%,'LineWidth',2

              *                     figure(1);

              *                     subplot(1,2,1);

              *                     plot(X,gs,'b',X,0,'black-');

              *                     hold on;

              *                     plot([xtmp; x],[g; zeros(size(g))],x,0,'rx',xtmp,g,'rx');

              *                     hold off;

              *                     axis tight;

              *

              *                     subplot(1,2,2);

              *                     X = linspace(0-this.xR,2*this.xR,si);

              *                     gs = zeros(length(y),si);

              *                     for idx = 1:length(y)

              *                         gs(idx,:) = kernels.BellFunction.optFun(X,repmat(y(idx),size(X)),f,df,ddf,this.x0,this.PenaltyFactor);

              *                     end

              *                     plot(X,gs);

              *                     hold on;

              *                     plot(X,0,'black-');

              *                     hold off;

              *                     pause;

              *                 end

              *             end */

         }


 #if 0 //mtoc++: 'set.r0'

 function  r0(value) {

             if ~isreal(value) || value < 0 || ~isscalar(value)

                 error(" r0 must be a scalar greater than zero. ");

             end

             this.fr0= value;

             this.setConstants;

         }


 #endif


 #if 0 //mtoc++: 'get.r0'

 function value = r0() {

             value = this.fr0;

         }


 #endif


     private:


         function  setConstants() {


             /*  Set outer bound for 'r_s' positions */

             this.rm= this.evaluateScalar(0)*this.r0 /...

                 (this.evaluateScalar(0)-this.evaluateScalar(this.r0));

             /*  Set values needed for Newton iteration */

             q = [0 this.r0 this.rm];

             this.sc= this.evaluateScalar(q);

             this.d1= this.evaluateD1(q);

             this.d2= this.evaluateD2(q);

             this.p= q;

         }

     protected:


         function [g , dg , pi , pl , pr ] = optFun(r,s,n,dn) {

             g = zeros(size(r));

             d = g; c = d; dg = g;


             a = s < this.r0;

             pi = a & r < this.r0 | ~a & r > this.r0;

             pr = a & r > this.rm;

             pl = ~a & r < 0;

             std = ~(pi | pl | pr);


             /*  Standard case */

             rs = r(std); ss = s(std);

             g(std) = this.evaluateD1(rs) - (this.evaluateScalar(rs)-this.evaluateScalar(ss))./(rs-ss);

             g(isnan(g)) = 0;

             dg(std) = this.evaluateD2(rs) - g(std)./(rs-ss);

             ninf = isinf(dg);

             dg(ninf) = this.evaluateD2(rs(ninf))/2;


             /*  Penalty case */

             penalty = ~std;

             c(pi) = dn(2,pi).^2 ./ (4*n(2,pi));

             d(pi) = 2*n(2,pi)./dn(2,pi) - this.r0;

             c(pl) = dn(1,pl).^2 ./ (4*n(1,pl));

             d(pl) = 2*n(1,pl)./dn(1,pl);

             c(pr) = dn(3,pr).^2 ./ (4*n(3,pr));

             d(pr) = 2*n(3,pr)./dn(3,pr) - this.rm;


             g(penalty) = c(penalty) .* (r(penalty) + d(penalty)).^2;

             dg(penalty) = 2 * c(penalty) .* (r(penalty) + d(penalty));

             /*  Avoid exact matches! */

             g(dg == 0) = 0;

             dg(dg == 0) = 1;

         }


     public: /* ( Abstract ) */


         virtual function dr = evaluateD1() = 0;

         virtual function ddr = evaluateD2() = 0;

     public: /* ( Static ) */


     protected: /* ( Static ) */


         static function this = loadobj(this,from) {

             if nargin > 1

                 /*  NOTE: r0 must be set via setter in subclass at some stage. */

                 if isfield(from," fr0 ") && ~isempty(from.fr0)

                     this.fr0= from.fr0;

                 end

                 this.MaxNewtonIterations= from.MaxNewtonIterations;

                 this.NewtonTolerance= from.NewtonTolerance;

                 this = loadobj@KerMorObject(this, from);

             elseif ~isa(this, " kernels.BellFunction ")

                 error(" Object passed is not a kernels.BellFunction instance and no init struct is given. ");

             end

             if isempty(this.fr0)

                 error(" Internal r0 field not persisted yet and not initialized in subclass loadobj method. ");

             end

             this.setConstants;

         }

 };

 }


kernels.BellFunction.getGlobalLipschitz
function   c  = getGlobalLipschitz()
Computes the absolute value of the first derivative at x0 Implements the template method from BaseKer...
Definition: BellFunction.m:145

kernels.BellFunction.ModifiedNewton
function   rtmp  = ModifiedNewton(rstart, s)
Definition: BellFunction.m:170

kernels.BellFunction.clone
function   copy  = clone(copy)
The interface method with returns a copy of the current class instance.
Definition: BellFunction.m:156

kernels.BellFunction
BELLFUNCTION Summary of this class goes here Detailed explanation goes here.
Definition: BellFunction.m:18

kernels.BellFunction.loadobj
static function   this  = loadobj(this, from)
As the constant properties are transient, they have to be re-computed upon loading.
Definition: BellFunction.m:371

DPCMObject.registerProps
function  registerProps(varargin)
Call this method at any class that defines DPCM observed properties.
Definition: DPCMObject.m:125

kernels.BellFunction.evaluateD2
virtual function   ddr  = evaluateD2()
Method for second derivative evaluation.

kernels.BellFunction.rm
rm
The maximum ("right") value for any .
Definition: BellFunction.m:109

kernels.BellFunction.MaxNewtonIterations
MaxNewtonIterations
Hard break for iteration count of modified newton algorithm.
Definition: BellFunction.m:92

KerMorObject.KerMorObject
KerMorObject()
Constructs a new KerMor object.
Definition: KerMorObject.m:55

kernels.BellFunction.r0
r0
Point of maximum first derivative on scalar evaluation.
Definition: BellFunction.m:61

kernels.BellFunction.NewtonTolerance
NewtonTolerance
Error tolerance for modified newton iteration.
Definition: BellFunction.m:77

kernels.ARBFKernel.evaluateScalar
virtual function   phir  = evaluateScalar(matrix< double > r)
Allows the evaluation of the function  for scalar  directly.

kernels.BellFunction.optFun
function [   g ,   dg ,   pi ,   pl ,   pr  ] = optFun(r, s, n, dn)
Definition: BellFunction.m:310

kernels.BellFunction.evaluateD1
virtual function   dr  = evaluateD1()
Method for first derivative evaluation.

kernels.ARBFKernel
Abstract class for radial basis function / rotation- and translation invariant kernels.
Definition: ARBFKernel.m:18

kernels.BellFunction.BellFunction
BellFunction()
Definition: BellFunction.m:139