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

        dscomponents.ACoreFun f;
        dscomponents.AffLinCoreFun,dscomponents.LinearCoreFun A = "[]";
        dscomponents.AInputConv B = "[]";
        dscomponents.LinearOutputConv C;
        dscomponents.AInitialValue x0 = "[]";
        dscomponents.AMassMatrix M = "[]";
        dscomponents.ACoreFun g = "[]";
        Inputs = {""};
        Params = data.ModelParam.empty;
        integer DependentParamIndices = -1;
        double MaxTimestep = "[]";
        colvec StateScaling = 1;
        SparsityPattern;
    public: /* ( Transient ) */

        mu = "[]";
        u = "[]";
        inputidx = "[]";
    public:

        NumStateDofs = "[]";

        NumAlgebraicDofs = 0;

        NumTotalDofs = "[]";

        
    public: /* ( Dependent ) */

        InputCount;
        ParamCount;
    public:

        Model;
    public:

        BaseFirstOrderSystem(models.BaseFullModel model) {
            this.validateModel(model);
            this.Model= model;
            this.C= dscomponents.LinearOutputConv(1);

            /*  Register default properties */
            this.registerProps(" A "," f "," B "," C "," x0 "," Inputs "," Params "," MaxTimestep "," StateScaling ");
        }
        function rsys = getReducedSystemInstance(models.ReducedModel rmodel) {
            rsys = models.ReducedSystem(rmodel);
        }
        function  setConfig(mu,inputidx) {
            if isempty(mu) && ~isempty(this.Params)
                error(" A model with parameters cannot be simulated without a parameter argument. ");
            end
            if size(mu,1) ~= this.ParamCount
                error(" The mu vector size mismatches the defined parameter number. ");    
            end
            this.mu= mu;

            if isempty(inputidx) && ~isempty(this.Inputs)
                cprintf(" red "," Attention: Starting simulations without input, but the system has configured some.\n ");
                /* inputidx = 1; */
            elseif inputidx > length(this.Inputs)
                error(" Invalid input index '%d'. There are %d possible inputs configured. ",inputidx,length(this.Inputs));
            end
            this.u= @(t)[];
            this.inputidx= [];
            if ~isempty(inputidx)
                this.u= this.Inputs[inputidx];
                this.inputidx= inputidx;
            end
        }
        function  prepareSimulation(colvec<double> mu,integer inputidx) {
            this.setConfig(mu, inputidx);

            if isempty(this.NumStateDofs)
                this.NumStateDofs= size(this.x0.evaluate(mu),1);
                fprintf(2," NumStateDofs not set. Guess from initial value: %d\n ",this.NumStateDofs);
                this.updateDimensions;
            end

            /*  Forward preparation call to nonlinearity, if present */
            if ~isempty(this.A)
                this.A.prepareSimulation(mu);
            end

            /*  Forward preparation call to nonlinearity, if present */
            if ~isempty(this.f)
                this.f.prepareSimulation(mu);
            end

            /*  Forward preparation call to nonlinearity, if present */
            if ~isempty(this.B)
                this.B.prepareSimulation(mu);
            end
        }


        function dx = ODEFun(double t,colvec<double> x) {
            dx = zeros(this.NumTotalDofs,1);
            if ~isempty(this.A)
                dx = dx + this.A.evaluate(x, t);
            end
            if ~isempty(this.f)
                dx = dx + this.f.evaluate(x, t);
            end
            if ~isempty(this.B) && ~isempty(this.inputidx)
                dx = dx + this.B.evaluate(t, this.mu)*this.u(t);
            end
            if ~isempty(this.g)
                dx = [dx; this.g.evaluate(x,t)];
            end
        }
        function y = computeOutput(matrix<double> x,colvec<double> mu) {

            if nargin < 3
                mu = this.mu;
            end

            /*  Re-scale state variable */
            if ~isequal(this.StateScaling,1)
                if isscalar(this.StateScaling)
                    x = this.StateScaling*x;
                else
                    x = bsxfun(@times,x,this.StateScaling);
                end
            end
            y = x;
            if ~isempty(this.C)
                if this.C.TimeDependent
                    /*  Evaluate the output conversion at each time t
                     * Figure out resulting size of C*x evaluation */
                    t = this.Model.scaledTimes;
                    hlp = this.C.evaluate(t(1),mu)*x(:,1);
                    y = zeros(size(hlp,1),length(t));
                    y(:,1) = hlp;
                    for idx=2:length(t)
                        y(:,idx) = this.C.evaluate(t(idx),mu)*x(:,idx);
                    end
                else
                    /*  otherwise it's a constant matrix so multiplication
                     * can be preformed much faster. */
                    y = this.C.evaluate([],mu)*x;
                end
            end
        }
        function  updateSparsityPattern() {
            td = this.NumTotalDofs;
            sd = this.NumStateDofs;
            JP = logical(sparse(td,td));
            if ~isempty(this.A) && ~isempty(this.A.JSparsityPattern)
                JP(1:sd,1:sd) = JP(1:sd,1:sd) | this.A.JSparsityPattern;
            end
            if ~isempty(this.f) && ~isempty(this.f.JSparsityPattern)
                JP(1:sd,:) = JP(1:sd,:) | this.f.JSparsityPattern;
            end
            if ~isempty(this.g)
                JP(sd+1:end,:) = this.g.JSparsityPattern;
            end
            this.SparsityPattern= JP;
        }


        function J = getJacobian(double t,xc) {
            td = this.NumTotalDofs;
            sd = this.NumStateDofs;
            x = xc(1:sd);
            J = sparse(td,td);
            if isempty(this.g)
                if ~isempty(this.A)
                    J = J + this.A.getStateJacobian(x, t);
                end
                if ~isempty(this.f)
                    J = J + this.f.getStateJacobian(xc, t);
                end
            else
                if ~isempty(this.A)
                    J(1:sd,1:sd) = J(1:sd,1:sd) + this.A.getStateJacobian(x, t);
                elseif ~isempty(this.f)
                    J(1:sd,:) = J(1:sd,:) + this.f.getStateJacobian(xc, t);
                end
                J(sd+1:end,:) = this.g.getStateJacobian(xc,t);
            end
        }
        function matrix<double>x0 = getX0(rowvec<double> mu) {
            ss = this.StateScaling;
            if (size(mu,2) > 1) && ~isscalar(ss)
                ss = repmat(ss,1,size(mu,2));
            end
            x0 = this.x0.evaluate(mu) ./ ss;

            /*  Append DoFs for algebraic conditions below (if set) */
            if this.NumAlgebraicDofs > 0
                x0 = [x0; this.getAlgebraicDofsInitialConditions];
            end
        }
        function M = getMassMatrix() {
            M = this.M;
            if isempty(M)
                M = dscomponents.ConstMassMatrix(speye(this.NumStateDofs));
            else
                if this.NumAlgebraicDofs > 0
                    error(" First order with AlgebraicDofs not yet implemented. ");
                end
            end
        }
        function ModelParamp = addParam(char name,default,varargin) {
            if ~isempty(this.Params) && ~isempty(find(strcmpi(name,[this.Params(:).Name]),1))
                error(" Parameter with name %s already exists. Use setParam instead. ",name);
            end
            p = data.ModelParam(name, default, varargin[:]);
            this.Params(end+1) = p;
        }
        function pidx = getParamIndexFromName(paramname) {
            if ~ischar(paramname)
                error(" Parameter paramname must be a char array ");
            end
            pidx = find(strcmpi(paramname,[this.Params(:).Name]),1);
        }
        function pt = getParamInfo(colvec<double> mu) {
            if nargin < 2
                mu = this.Model.DefaultMu;
            end
            /* str = ''; */
            pt = PrintTable;
            pt.HasRowHeader= true;
            pt.HasHeader= true;
            pt.addRow(" Parameter name "," Value "," Index ");
            for i=1:this.ParamCount
               /*  str = sprintf('%s%d - %s: %f\n',str,i,this.Params(i).Name,mu(i)); */
                pt.addRow(this.Params(i).Name,mu(i),i);
            end
            if nargout < 1
                pt.print;
            end
        }


        function  plotInputs(pm) {
            if nargin < 2
                pm = PlotManager;
                pm.LeaveOpen= true;
            end
            for i=1:this.InputCount
                ui = this.Inputs[i];
                h = pm.nextPlot(sprintf(" u%d ",i),sprintf(" u_%d = %s ",i,func2str(ui)),...
                    " t ",sprintf(" u_%d ",i));
                plot(h,this.Model.Times,ui(this.Model.scaledTimes)," LineWidth ",2);
            end
            pm.done;
        }

    
    protected:


        function  updateDimensions() {
            this.NumTotalDofs= this.NumStateDofs + this.NumAlgebraicDofs;
        }


        function ad_ic = getAlgebraicDofsInitialConditions() {
            ad_ic = zeros(this.NumAlgebraicDofs,1);
        }
        function  validateModel(models.BaseFullModel model) {
            if ~isa(model, " models.BaseModel ")
                error(" The Model property has to be a child of models.BaseModel ");
            end
        }
        /* % Getter & Setter */
    public:

        
#if 0 //mtoc++: 'set.f'
function  f(value) {
            if ~isempty(value) && ~isa(value, " dscomponents.ACoreFun ")
                error(" The property 'f' has to be a class implementing dscomponents.ACoreFun ");
            end
/*              if (isequal(this,value))
 *                 warning('KerMor:selfReference','Careful with self-referencing classes. See BaseFullModel class documentation for details.');
 *             end */
            this.f= value;
        }

#endif


        
#if 0 //mtoc++: 'set.B'
function  B(value) {
            if ~isempty(value) && ~isa(value, " dscomponents.AInputConv ")
                error(" The property 'B' has to be a class implementing dscomponents.AInputConv ");
            end
            this.B= value;
        }

#endif


        
#if 0 //mtoc++: 'set.A'
function  A(value) {
            if ~isempty(value) && ~(isa(value, " dscomponents.LinearCoreFun ") || isa(value, " dscomponents.AffLinCoreFun "))
                error(" The property 'A' has to be a LinearCoreFun or AffLinCoreFun. ");
            end
            this.A= value;
        }

#endif


        
#if 0 //mtoc++: 'set.C'
function  C(value) {
            if isempty(value)
                value = dscomponents.LinearOutputConv(1);
                warning(" An output conversion must always exist. Choosing dscomponents.LinearOutputConv(1) for simple forwarding. ");
            elseif ~isa(value, " dscomponents.AOutputConv ")
                error(" The property 'C' has to be a class implementing dscomponents.AOutputConv ");
            end
            this.C= value;
        }

#endif


        
#if 0 //mtoc++: 'set.Inputs'
function  Inputs(value) {
            if ~iscell(value)
                error(" Property 'Inputs' must be a cell array. ");
            end
            for n=1:numel(value)
                if ~isempty(value[n])
                    if ~isa(value[n]," function_handle ")
                        error(" Each 'Inputs' cell must contain a function handle. ");
/*                      elseif nargin(value{n}) ~= 1
 *                        error('Each "Inputs" function must take exactly one (=time) parameter.'); */
                    end
                end
            end
            this.Inputs= value;
        }

#endif


        
#if 0 //mtoc++: 'set.Params'
function  Params(value) {
            if ~isa(value," data.ModelParam ")
                error(" Params property must be a ModelParam array. ");
            end
            this.Params= value;
        }

#endif


        
#if 0 //mtoc++: 'set.x0'
function  x0(value) {
            this.checkType(value," dscomponents.AInitialValue ");
/*              if ~isa(value,'function_handle')
 *                 error('x0 must be a function handle.');
 *             elseif nargin(value) ~= 1
 *                 error('x0 must take exactly one argument.');
 *             end */

            this.x0= value;
        }

#endif


        
#if 0 //mtoc++: 'get.ParamCount'
function value = ParamCount() {
            value = length(this.Params);
        }

#endif


        
#if 0 //mtoc++: 'get.InputCount'
function value = InputCount() {
            value = length(this.Inputs);
        }

#endif


        
#if 0 //mtoc++: 'set.MaxTimestep'
function  MaxTimestep(value) {
            if (~isscalar(value) || value <= 0) && ~isempty(value)
                error(" Value must be a positive real scalar if not empty. ");
            end
            this.MaxTimestep= value;
        }

#endif


        
#if 0 //mtoc++: 'set.StateScaling'
function  StateScaling(value) {
            if isempty(value)
                error(" StateScaling must not be empty. (Set to 1 if not needed). ");
            elseif ~isvector(value)
                error(" StateScaling must be a vector. ");
            end
            this.StateScaling= value(:);
        }

#endif

        
};
}