MLWrapper.m

Go to the documentation of this file.

namespace solvers{


/* (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 MLWrapper
  :public solvers.BaseSolver,
   public solvers.AJacobianSolver {
    public: /* ( setObservable ) */

        function_handle MLSolver = @ode23;
        struct odeopts;
    private:

        fED;
    public:


        MLWrapper(solver) {
            this = this@solvers.BaseSolver;
            this.registerProps(" MLSolver ");

            this.SolverType= solvers.SolverTypes.Explicit;
            if nargin > 0
                this.MLSolver= solver;
                slvname = char(solver);
                if ~isempty(strfind(" ode15s ",slvname))
                    this.SolverType= solvers.SolverTypes.Implicit;
                end
            end
            this.odeopts= odeset;
        }


        function [doublet , matrix<double>y ] = solve(odefun,double t,x0) {
            opts = this.odeopts;
            if ~isempty(this.MaxStep)
                opts = odeset(opts, " MaxStep ", this.MaxStep);
            end
            if ~isempty(this.InitialStep)
                opts = odeset(opts, " InitialStep ", this.InitialStep);
            end

            /*  Use AJacobianSolver properties */
            opts = odeset(opts, " Jacobian ", this.JacFun);
            JP = [];
            /*  Matlab solvers only use patterns if the jacobian function is
             * not directly given */
            if isempty(this.JacFun)
                JP = this.JPattern;
            end
            opts = odeset(opts, " JPattern ", JP);

            /*  Pass Mass Matrix and initial slope to solver (if set) */
            M = []; yp0 = odefun(0,x0); MS = [];
            if ~isempty(this.M)
                if ~this.M.TimeDependent
                    M = this.M.evaluate(0);
                    Mt0 = M;
                else
                    M = @(t,notused)this.M.evaluate(t);
                    Mt0 = M(0);
                end

                /*  See if we have a "simple" singular mass matrix */
                algdims = this.M.AlgebraicEquationDims;
                if ~isempty(algdims)
                    invertible = 1:size(Mt0,1);
                    invertible(algdims) = [];
                    yp0(invertible) = Mt0(invertible,invertible)\yp0(invertible);
                else
                    /*  Compute initial slope (considering possibly
                     * rank-deficient M) */
                    if issparse(Mt0)
                        r = sprank(Mt0);
                    else
                        r = rank(Mt0);
                    end
                    if r < size(Mt0,1)
                        odeset(opts," MassSingular "," yes ");
                        warning(" off "," MATLAB:singularMatrix ");
                    end
                    yp0 = Mt0\yp0;
                    if r < size(Mt0,1)
                        warning(" on "," MATLAB:singularMatrix ");
                        yp0(isinf(yp0)) = 0;
                        if any(isnan(yp0(:)))
                            error(" Failed to automatically produce a consistent initial slope. ");
                        end
                    end
                end
                MS = " none ";
            end
            opts = odeset(opts," Mass ",M," MStateDependence ",MS,...
                    " InitialSlope ",yp0);

            if this.RealTimeMode
                opts = odeset(opts," OutputFcn ",@this.ODEOutputFcn);
                /*  Bug in Matlab 2009a: direct assignment crashes Matlab!
                 * Seems also not to work if created within the constructor. */
                ed = solvers.SolverEventData;
                this.fED= ed;
                this.solverCall(odefun, t, x0, opts);
                t = []; y = [];
            else
                [t,y] = this.solverCall(odefun, t, x0, opts);
                y = y^t;
                t = t^t;
            end
        }


        
#if 0 //mtoc++: 'set.MLSolver'
function  MLSolver(value) {
            if isa(value," function_handle ")
                this.MLSolver= value;
                this.Name= sprintf(" Matlab Solver wrapper using %s ",func2str(value));
            else
                error(" Invalid function handle! ");
            end
        }

#endif

    
    protected:

        function varargout = solverCall(function_handle odefun,rowvec<double> t,x0,struct opts) {
            [varargout[1:nargout]] = this.MLSolver(odefun, t, x0, opts);
        }
    protected: /* ( Sealed ) */

        function status = ODEOutputFcn(t,y,flag) {
            if ~strcmp(flag," init ")
                /*  For some reason the t and y args have more than one
                 * entry, so loop over all of them. */
                for idx=1:length(t)
                    this.fED.Times= t(idx); /*  when flag==init the t var is larger than one */

                    this.fED.States= y(:,idx);
                    this.notify(" StepPerformed ",this.fED);
                end
            end
            status = 0;
        }
};
}