MarkertLaw.m

Go to the documentation of this file.

namespace models{
namespace muscle{
namespace functions{


/* (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 MarkertLaw
  :public general.functions.AFunGen {
    public: /* ( Constant ) */

        static const LatestLinearizationLambda = 5;

        
    public:

        b;

        d;

        max_modulus;

        t0 = "[]";

        ft0 = "[]";

        
    public:

        MarkertLaw(b,d,max_modulus) {
            if nargin < 3
                max_modulus = [];
                /*  Take data from fit to original function over [1 1.3]
                 * -> general.functions.MarkertLaw.test_FitToOriginal(general.functions.MarkertLawOriginal,[1 1.3]); */
                if nargin < 2
                    d = 0.1660;
                    if nargin < 1
                        b = 1.3720e8;
                    end
                end
            end
            this.b= b;
            this.d= d;
            if ~isempty(max_modulus) /* && d > 1 */

                [f,df] = this.getFunction;
                rightval = 2;
                found = false;
                while ~found
                    pts = linspace(1,rightval,3);
                    /*  Compute maxit to the numer of iterations where the
                     * interval length divided by 2^maxit < eps */
                    maxit = ceil(log((pts(3)-pts(1))/eps)/log(2));
                    it = 0; diff = Inf;
                    while diff > sqrt(eps) && it < maxit;
                        fpts = df(pts)-max_modulus;
                        if fpts(2) < 0
                            pts = [pts(2) (pts(2)+pts(3))/2 pts(3)];
                        else
                            pts = [pts(1) (pts(1)+pts(2))/2 pts(2)];
                        end
                        diff = abs(fpts(2));
                        it = it+1;
                        if pts(2) > this.LatestLinearizationLambda
                            found = true;
                            break;
                        end
                    end
                    if rightval-pts(2) < 1e-10
                        rightval = rightval*2;
                        /* fprintf('Increasing right interval value to %d\n',rightval); */
                    else
                        found = true;
                    end
                end
                this.t0= pts(2);
                fprintf(" Found t0=%5.3g with diff %g to max modulus/derivative of %g after %d iterations (b=%g, d=%g).\n ",...
                    this.t0,diff,max_modulus,it,b,d);
                this.ft0= f(this.t0);
                this.max_modulus= max_modulus;
            end
        }


        function [fhandle , dfhandle , fbdhandle , dfbdhandle ] = getFunction() {
            b1 = this.b;
            d1 = this.d;
            tlin = this.t0;
            if ~isempty(tlin)
                mm = this.max_modulus;
                ftlin = this.ft0;
                /* fhandle = @(t)(t<tlin).*max(0,b1.*(t.^d1-1)) + (t>=tlin).*(ftlin + mm*(t-tlin));
                 *dfhandle = @(t)(t<tlin).*max(0,b1.*((d1-2).*t.^d1+2)) + (t>=tlin)*mm; */
                fhandle = @(t)(t >= 1 & t<tlin).*(b1./t.^2.*(t.^d1-1).*(t-1).^2)...
                    + (t>=tlin).*(ftlin + mm*(t-tlin));
                dfhandle = @(t)(t >= 1 & t<tlin).*(b1.*(t-1)./t.^3.*((d1*(t-1)+2).*t.^d1-2)) ...
                    + (t>=tlin)*mm;
                fbdhandle = @(t,b1,d1)(t >= 1 & t<tlin).*(b1./t.^2.*(t.^d1-1).*(t-1).^2)...
                    + (t>=tlin).*(ftlin + mm*(t-tlin));
                dfbdhandle = @(t,b1,d1)(t >= 1 & t<tlin).*(b1.*(t-1)./t.^3.*((d1*(t-1)+2).*t.^d1-2)) ...
                    + (t>=tlin)*mm;
            else
/*                  fhandle = @(t)max(0,b1.*(t.^d1-1));
 *                 dfhandle = @(t)max(0,b1.*((d1-2).*t.^d1+2)); */
                fhandle = @(t)(t>=1).*(b1./t.^2.*(t.^d1-1).*(t-1).^2);
                dfhandle = @(t)(t>=1).*b1.*(t-1)./t.^3.*((d1*(t-1)+2).*t.^d1-2);
                fbdhandle = @(t,b1,d1)(t>=1).*(b1./t.^2.*(t.^d1-1).*(t-1).^2);
                dfbdhandle = @(t,b1,d1)(t>=1).*b1.*(t-1)./t.^3.*((d1*(t-1)+2).*t.^d1-2);
            end
        }


        function str = getConfigStr() {
            str = sprintf(" b: %g, d: %g, MaxModulus:%g ",this.b,this.d,this.max_modulus);
        }


        function  plot(range,varargin) {
            if nargin < 2
                if ~isempty(this.t0)
                    range = [1 1.2*this.t0];
                else
                    range = [1 1.2];
                end
            end
            plot@general.functions.AFunGen(this, " R ", range, varargin[:]);
            if (range(2) > this.t0)
                ax = get(gcf," Children ");
                ax = ax(2); /*  second one is the left one - hope this is reproducible */

                hold(ax," on ");
                plot(ax,this.t0,this.ft0," rx "," MarkerSize ",16);
            end
        }

    
    public: /* ( Static ) */

        static function  test_Linarization() {
            C = Utils.createCombinations([1 5],linspace(3,40,5),linspace(100,1e6,5));
            x = linspace(.3,6,300);
            pm = PlotManager(false,1,2);
            pm.LeaveOpen= true;
            ax = pm.nextPlot(" markertfun_lintest "," Plots for various parameters "," lambda "," value ");
            axl = pm.nextPlot(" markertfun_lintest_log "," Plots for various parameters "," lambda "," value ");
            for i = 1:size(C,2)
                ml = models.muscle.functions.MarkertLaw(C(1,i),C(2,i),C(3,i));
/*                  ml.plot(x);                 */
                [f,df] = ml.getFunction;
                fx = f(x);
                plot(ax,x,fx," b ");
                semilogy(axl,x,fx," b ");
                fprintf(" t0=%g, f(t0)=%g for %s\n ",ml.t0, ml.ft0, ml.getConfigStr);

                hold(ax," on ");
                hold(axl," on ");

                ml = models.muscle.functions.MarkertLaw(C(1,i),C(2,i),[]);
/*                  ml.plot(x); */
                [f,df] = ml.getFunction;
                fx = f(x);
                plot(ax,x,fx," r ");
                semilogy(axl,x,fx," r ");
            end
            pm.done;
            ylim(ax,[0 1e6]);
        }


        static function [res , fit , fitlin ] = test_FitToOriginal(orig,range) {
            if nargin < 2
                range = [1 1.2];
                if nargin < 1
                    orig = models.muscle.functions.MarkertLawOriginal(7990, 16.6);
                end
            end
            ntest = 100;
            fac = 100;
            p = Utils.createCombinations(...
                linspace(orig.b/fac,orig.b*fac*1000,ntest),...
                linspace(orig.d/fac,orig.d*fac,ntest));
            l = linspace(range(1),range(2),100);
            of = orig.getFunction;
            n = size(p,2);
            forig = of(l);
            err = zeros(1,n);
            pi = ProcessIndicator(" Performing fit for %d cases ",n,false,n);
            for k=1:n
                m = models.muscle.functions.MarkertLaw(p(1,k),p(2,k));
                mf = m.getFunction;
                ftest = mf(l);
                err(k) = norm(forig-ftest);
                pi.step;
            end
            pi.stop;
            [~, idx] = min(err(1,:));
            fit = models.muscle.functions.MarkertLaw(p(1,idx),p(2,idx));
            fitlin = models.muscle.functions.MarkertLaw(p(1,idx),p(2,idx),1e5);
            orig.plot(l);
            fit.plottofigure(gcf,l," P ",[" r "]);
            fitlin.plottofigure(gcf,l," P ",[" g "]);
            res = true;
        }

    

};
}
}
}