BaseModel.m

Go to the documentation of this file.

namespace models{


/* (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 BaseModel
  :public KerMorObject {
    public: /* ( setObservable ) */

        models.BaseFirstOrderSystem System = "[]";
        char Name = "";
        matrix<double> G = 1;
        double RealTimePlottingMinPause = .1;
        logical isStatic = false;
        integer DefaultInput = "[]";
    public: /* ( Dependent ) */

        Times;
        rowvec<double> scaledTimes;
        double Tscaled;
    public: /* ( Dependent, setObservable ) */

        double tau;
        double T;
        double dt;
        solvers.BaseSolver ODESolver;
        logical RealTimePlotting;
        colvec<double> DefaultMu = "[]";
    public:

        dtscaled = .1;
        char gitRefOnSave = "";
    private:

        struct ftau = 1;
        fdt = .1;

        fT = 1;

        frtp = false;

        fODEs;

        steplistener;

        ctime;

        fDefMu = "[]";

        
    public:


        BaseModel() {
            
            /*  Call superclass constructor first
             * (not necessary in this version as automatically called first, but one never knows..) */
            this = this@KerMorObject;

            /*  Init defaults */
            this.ODESolver= solvers.MLWrapper(@ode23);

            /*  Register default properties */
            this.registerProps(" System "," T "," ODESolver "," dt "," G ",...
                " tau "," RealTimePlottingMinPause "," RealTimePlotting ");
        }


        function  delete() {
           this.ODESolver= [];
           this.System= [];
        }


        function  initDefaultParameter() {
            mu = zeros(this.System.ParamCount,1);
            /* p = struct; */
            for k = 1:this.System.ParamCount
                mu(k) = this.System.Params(k).Default;
                /* name = regexprep(this.System.Params(k).Name,'[^A-Za-z0-9_]', '_');
                 *p.(name) = k; */
            end
            this.fDefMu= mu;
            /* this.ParamIdx = p; */
        }
        function [rowvec<double>t , matrix<double>y , doublesec , x ] = simulate(colvec<double> mu,integer inputidx) {
            if nargin < 3
                inputidx = this.DefaultInput;
                if nargin < 2
                    mu = this.DefaultMu;
                end
            end
            /*  Transpose if necessary */
            if size(mu,2) > 1 && size(mu,1) == 1
                mu =mu^t;
            end
            this.WorkspaceVariableName= inputname(1);

            if this.RealTimePlotting
                this.ctime= tic;
            end
            if this.isStatic
                /*  solve static equation */
                [t, x, time] = this.solveStatic(mu, inputidx);
            else

                /*  Get scaled state trajectory */
                [t, x, time] = this.computeTrajectory(mu, inputidx);
            end
            /*  Measure rest of time */
            starttime = tic;

            /*  Convert to output */
            y = this.System.computeOutput(x);

            /*  Scale times back to real units */
            t = t*this.tau;

            sec = toc(starttime) + time;

            this.WorkspaceVariableName= ;
        }
        function [handlef , handleax ] = plot(rowvec<double> t,matrix<double> y,varargin) {
            if isempty(varargin)
                f = figure;
                ax = gca(f);
            else
                ax = varargin[1];
                f = get(ax," Parent ");
            end
            y = Utils.preparePlainPlot(y);
            plot(ax,t,y);
            title(ax,sprintf(" Output plot for model '%s' ", this.Name));
            xlabel(ax," Time "); ylabel(ax," Output values ");
        }
        function [handlef , handleax ] = plotState(rowvec t,matrix<double> x,varargin) {
            if isempty(varargin)
                f = figure;
                ax = gca(f);
            else
                ax = varargin[1];
                f = get(ax," Parent ");
            end
            x = Utils.preparePlainPlot(x);
            plot(ax,t,x);
            title(ax,sprintf(" State space plot for model '%s' ", this.Name));
            xlabel(ax," Time "); ylabel(ax," State space values ");
        }
        function [f , ax ] = plotSingle(double t,colvec<double> y,varargin) {
            [f,ax] = this.plot(t, y, varargin[:]);
        }
        function  plotInputs(pm) {
            if nargin < 2
                pm = PlotManager;
                pm.LeaveOpen= true;
            end
            h = pm.nextPlot(" inputs "," Model inputs "," t "," u(t) ");
            hold(h," on ");
            leg = [];
            cg = LineSpecIterator;
            for k=1:this.System.InputCount
                plot(h,this.Times,this.System.Inputs[k](this.scaledTimes)," Color ",cg.nextColor);
                leg[k] = sprintf(" u_%d(t) ",k);/* #ok */

            end
            legend(h,leg[:]);
            pm.done;
        }


        function [doublet , colvec<double>x , doublectime ] = solveStatic(colvec<double> mu,integer inputidx) {
            if nargin < 3
                inputidx = [];
                if nargin < 2
                    mu = [];
                end
            end
            sys = this.System;
            /*  Stop the time */
            st = tic;

            /*  Prepare the system by setting mu and inputindex. */
            sys.prepareSimulation(mu, inputidx);
            t = this.scaledTimes;

            /*  Prepare the right-hand side */
            rhs = zeros(size(sys.A.evaluate(1,0),1), length(t));
            if ~isempty(sys.B)
                if ~sys.B.TimeDependent
                    rhs = sys.B.evaluate(0, sys.mu)*sys.u(t);
                else
                    for tdx = 1:length(t)
                        rhs(:,tdx) = sys.B.evaluate(t(tdx), sys.mu)*sys.u(t(tdx));
                    end
                end
            end
            if ~isempty(sys.f)
                rhs = rhs + sys.f.evaluateMulti(zeros(sys.f.xDim,0), t, repmat(sys.mu,1,length(t)));
            end

            /*  solve the system A*x + rhs = 0 for x */
            if ~sys.A.TimeDependent
                /*  precompute lu decomposition */
                [l,u] = lu(sys.A.evaluate(1,0,sys.mu));
                x = -u\(l\rhs);
            else
                x = zeros(sys.A.xDim, length(t));
                for tdx = 1:length(t) 
                    x(:,tdx) = -sys.A.evaluate(1,t(tdx),sys.mu)\rhs(:,tdx);
                end
            end

            ctime = toc(st);
        }
        function [rowvec<double>t , matrix<double>x , doublectime ] = computeTrajectory(colvec<double> mu,integer inputidx) {
            if nargin < 3
                inputidx = this.DefaultInput;
                if nargin < 2
                    mu = this.DefaultMu;
                end
            end
            sys = this.System;

            /*  Stop the time */
            st = tic;

            /*  Prepare the system by setting mu and inputindex. */
            sys.prepareSimulation(mu, inputidx);

            /*  Check explicit solvers */
            if isempty(this.System.MaxTimestep) && this.fODEs.SolverType ~= solvers.SolverTypes.Implicit
                warning(" Attention: Using an non-implicit solver without System.MaxTimestep set. Please check. ");
            end

            /* % Solve ODE */
            slv = this.ODESolver;
            slv.MaxStep= []; slv.InitialStep= [];
            if ~isempty(sys.MaxTimestep)
                /*  Remember: When scaling is used, these are the  */
                slv.MaxStep= sys.MaxTimestep/this.tau;
                slv.InitialStep= .5*sys.MaxTimestep/this.tau;
            end

            /*  Assign jacobian information if available */
            if isa(slv," solvers.AJacobianSolver ")
                slv.JPattern= sys.SparsityPattern;
                slv.JacFun= @sys.getJacobian;
            end

            /*  Assign mass matrix to solver if present */
            slv.M= [];
            if ~isempty(sys.M)
                slv.M= sys.getMassMatrix;
            end

            /*  Call solver */
            [t, x] = slv.solve(@sys.ODEFun, this.scaledTimes, sys.getX0(mu));

            if length(this.scaledTimes) == 2
                t = [t(1), t(end)];
                x = [x(:,1), x(:,end)];
            end

            /*  Get used time */
            ctime = toc(st);
        }
        function matrix<double>mu = getRandomParam(integer num,integer seed) {
            if nargin < 3
                seed = round(cputime*100);
                if nargin < 2
                    num = 1;
                end
            end
            s = this.System;
            r = RandStream(" mt19937ar "," Seed ",seed);
            if s.ParamCount > 0
                pmin = [s.Params(:).MinVal]^t;
                pmax = [s.Params(:).MaxVal]^t;
                mu = r.rand(s.ParamCount,num) .* repmat(pmax-pmin,1,num) + repmat(pmin,1,num);
            else
                mu = [];
            end
        }
    protected:


        function this = saveobj() {
            this.gitRefOnSave= KerMor.getGitBranch;
        }
    protected: /* ( Static ) */

        static function this = loadobj(this,s) {
            if ~isa(this," models.BaseModel ") && nargin < 2
                error(" The model class has changed but the loadobj method does not implement backwards-compatible loading behaviour.\nPlease implement the loadobj-method in your subclass and pass the loaded object struct as second argument. ");
            end
            if nargin == 2
                this.System= s.System;
                this.Name= s.Name;
                this.G= s.G;
                this.RealTimePlottingMinPause= s.RealTimePlottingMinPause;
                this.ftau= s.ftau;
                this.fdt= s.fdt;
                this.fT= s.fT;
                this.frtp= s.frtp;
                this.fODEs= s.fODEs;
                this.steplistener= s.steplistener;
                this.ctime= s.ctime;
                this.dtscaled= s.dtscaled;
                this.gitRefOnSave= s.gitRefOnSave;                
                if isfield(s," DefaultMu ")
                    this.DefaultMu= s.DefaultMu;
                    this.DefaultInput= s.DefaultInput;
                end
            end
            this = loadobj@DPCMObject(this);
        }

    
    protected: /* ( Sealed ) */

        function  plotstep(src,ed) {
            y = this.System.computeOutput(ed.States);
            this.plotSingle(ed.Times * this.tau,y);
            drawnow;
            /*  Gets first set in "simulate" */
            per = toc(this.ctime);
            pause(max(0,this.RealTimePlottingMinPause-per));
            this.ctime= tic;
        }
        /* % Getter & Setter */
    public:

        
#if 0 //mtoc++: 'get.Times'
function value = Times() {
            value = 0:this.dt:this.T;
        }

#endif


        
#if 0 //mtoc++: 'get.scaledTimes'
function value = scaledTimes() {
            value = 0:this.dtscaled:this.Tscaled;
        }

#endif


        
#if 0 //mtoc++: 'get.Tscaled'
function value = Tscaled() {
            if isempty(this.tau)
                value = this.T;
            else
                value = this.T/this.tau;
            end
        }

#endif


        
#if 0 //mtoc++: 'get.dt'
function dt = dt() {
            dt = this.fdt;
        }

#endif


        
#if 0 //mtoc++: 'get.tau'
function tau = tau() {
            tau = this.ftau;
        }

#endif


        
#if 0 //mtoc++: 'set.T'
function  T(value) {
            if ~isscalar(value) || value < 0
                error(" T must be a positive real scalar. ");
            elseif value < this.fdt
                warning(" Timestep dt must be smaller or equal to T. Using dt=%g ",value/2);
                this.fdt= value/2;
            end
            if value ~= this.fT
                this.fT= value;
            end
        }

#endif


        
#if 0 //mtoc++: 'set.dt'
function  dt(value) {
            if ~isscalar(value) || value <= 0
                error(" dt must be a positive real scalar. ");
            elseif value > this.fT
                error(" Timestep dt must be smaller or equal to T. ");
            end
            if this.fdt ~= value
                this.dtscaled= value/this.ftau;
                this.fdt= value;
            end
        }

#endif


        
#if 0 //mtoc++: 'set.tau'
function  tau(value) {
            if ~isreal(value) ||~isscalar(value)
                error(" tau must be a positive real scalar. ");
            end
            if this.ftau ~= value
                this.dtscaled= this.fdt/value;
                this.ftau= value;
            end
        }

#endif


        
#if 0 //mtoc++: 'set.System'
function  System(value) {
            if ~isempty(value) && ~isa(value," models.BaseFirstOrderSystem ")
                error(" The System property must be a subclass of models.BaseFirstOrderSystem. ");
            end
            if (isequal(this,value))
                warning(" KerMor:selfReference "," Careful with self-referencing classes. See BaseFullModel class documentation for details. ");
            end
            this.System= value;
        }

#endif

#if 0 //mtoc++: 'set.G'
function  G(value) {
            if isempty(value)
                error(" G must not be empty. Use G=1 for default euclidean scalar product and norms ");
            elseif any(abs(value-value^t) > eps)
                error(" G must be symmetric. ");
            end
            this.G= value;
        }

#endif

#if 0 //mtoc++: 'set.ODESolver'
function  ODESolver(value) {
            this.checkType(value," solvers.BaseSolver ");
            /*  Add listener if new ODE solver is passed and real time
             * plotting is turned on. */
            if this.frtp && this.fODEs ~= value
                if ~isempty(this.steplistener)
                    delete(this.steplistener);
                end
                this.steplistener= value.addlistener(" StepPerformed ",@this.plotstep);
            end
            this.fODEs= value;
        }

#endif


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

#endif


        
#if 0 //mtoc++: 'set.Name'
function  Name(value) {
            if ~ischar(value)
                error(" name is acharacter field ");
            end
            this.Name= value;
        }

#endif


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

#endif


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

#endif


        
#if 0 //mtoc++: 'set.RealTimePlotting'
function  RealTimePlotting(value) {
            if ~islogical(value)
                error(" Value must be boolean. ");
            end
            if value
                if isempty(this.steplistener)
                    this.steplistener= this.ODESolver.addlistener(" StepPerformed ",@this.plotstep);
                end
            else
                if ~isempty(this.steplistener)
                    delete(this.steplistener);
                    this.steplistener= [];
                end
            end
            this.frtp= value;
        }

#endif


        
#if 0 //mtoc++: 'set.DefaultMu'
function  DefaultMu(colvec<double> mu) {
            if size(mu,2) > 1
                error(" Default param must be a column vector ");
            end
            this.fDefMu= mu;
        }

#endif


        
#if 0 //mtoc++: 'get.DefaultMu'
function colvec<double>mu = DefaultMu() {
            if isempty(this.fDefMu)
                this.initDefaultParameter;
            end
            mu = this.fDefMu;
        }

#endif

};
}