JacCompEvalWrapper.m

Go to the documentation of this file.

namespace general{


/* (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 JacCompEvalWrapper
  :public dscomponents.ACompEvalCoreFun {
    public:

        f;

        trafo = {""};

        
    public:

        fullxDim;

        
    public:

        JacCompEvalWrapper(f) {

          /*  Checks */
          if ~isa(f," dscomponents.ACompEvalCoreFun ")
              error(" The system "" s nonlinearity must be a component-wise evaluable function. ");
          end

          /*  Call superclass constructor */
          this = this@dscomponents.ACompEvalCoreFun(f.System);

          /*  Original full system function */
          this.f= f;

          /*  ACoreFun settings */
          this.xDim= f.xDim;
          /*  output dim is size of jacobian */
          if ~isempty(f.JSparsityPattern)
              fdim = length(find(f.JSparsityPattern));
          else
              fdim = f.fDim * f.xDim;
          end
          this.fDim= fdim;
          /*  So sparsity pattern for the wrapper */
          this.JSparsityPattern= [];

          /*  Extra: Store original full dimension, as f gets possibly
           * projected */
          this.fullxDim= f.xDim;
        }
        function colvec<double>fx = evaluate(colvec<double> x,double t) {

            if ~isempty(this.f.JSparsityPattern)
                nonzero = find(this.f.JSparsityPattern);
            else
                nonzero = 1:this.f.fDim*this.f.xDim;
            end
            J = this.f.getStateJacobian(x,t,this.mu);
            fx = J(nonzero);
        }
        function colvec<double>fx = evaluateMulti(colvec<double> x,double t,colvec<double> mu) {

            if ~isempty(this.f.JSparsityPattern)
                nonzero = find(this.f.JSparsityPattern);
            else
                nonzero = 1:this.f.fDim*this.f.xDim;
            end
            fx = zeros(length(nonzero),size(x,2));
            singlemu = size(mu,2) == 1;
            if singlemu
                oldmu = this.f.mu;
                this.f.prepareSimulation(mu);
            else
                for i=1:size(x,2)
                    if ~singlemu
                        this.f.prepareSimulation(mu(:,i));
                    end
                    J = this.f.getStateJacobian(x(:,i),t(i));
                    fx(:,i) = J(nonzero);
                end
            end
            if singlemu
                this.f.prepareSimulation(oldmu);
            end
        }
        function  prepareSimulation(colvec<double> mu) {
            prepareSimulation@dscomponents.ACompEvalCoreFun(this,mu);
            this.f.prepareSimulation(mu);
        }


        function  evaluateCoreFun() {
            error(" Boo. Do not call me. (directly implemented evaluate) ");
        }
        function  setPointSet(nr,pts) {

            /*  Transform indices to jacobian matrix indices
             * (the jac wrapper only "works" on nonzero jacobian entries) */
            if ~isempty(this.f.JSparsityPattern)
                jidx = find(this.f.JSparsityPattern);
                pts = jidx(pts)^t;
            end

            /*  Ensure row vectors */
            if size(pts,1) > 1
                pts = reshape(pts,1,[]);
            end
            if length(unique(pts)) ~= length(pts)
                error(" Points have to be unique. ");
            end

            /*  Get matrix indexing of the desired points
             * ind2sub direct replacement! */
            i = rem(pts-1, this.f.fDim)+1;
            j = (pts-i)/this.f.fDim+1;
            /*  Only use the i-th effective components of f, the j's
             * correspond to the derivatives of the i-th components with
             * respect to the j-th variable. */
            ui = unique(i);

            deriv = [];
            pos = [];
            for k = 1:length(ui)
                didx = find(i == ui(k));
                deriv[k] = j(didx); /* #ok<*AGROW> */

                pos = [pos didx];
            end
            /*  Can directly set point sets as f instance is clone of
             * original function */
            this.f.setPointSet(nr+2, ui, deriv);

            l = length(pos);
            this.trafo[nr] = sparse(pos, 1:l, ones(l,1),l,l);
        }
        function fx = evaluateComponentSet(nr,colvec<double> x,double t) {
            dfx = this.f.evaluateJacobianSet(nr+2, x, t);
            fx = this.trafo[nr]*dfx;
        }


        function fx = evaluateComponentSetMulti(nr,colvec<double> x,double t,colvec<double> mu) {
            dfx = this.f.evaluateJacobianSetMulti(nr+2, x, t, mu);
            fx = this.trafo[nr]*dfx;
        }


        function copy = clone() {
            copy = general.JacCompEvalWrapper(this.f);
            copy = clone@dscomponents.ACompEvalCoreFun(this, copy);
            copy.trafo= this.trafo;
            copy.fullxDim= this.fullxDim;
        }


        function proj = project(V,W) {
            proj = project@dscomponents.ACompEvalCoreFun(this, V, W, this.clone);
            /*  Project the wrapped function, too (ignores the extra copy.f
             * created in clone as project creates a clone anyways) */
            proj.f= this.f.project(V, W);
            proj.xDim= proj.f.xDim;
            /*  vec-op leads to product of jacobian size dimension */
            proj.fDim= proj.f.fDim*proj.f.xDim;
        }

    
    protected:

        function  evaluateComponents() {
            error(" This should not be called as this is a wrapper for ACompEvalCoreFuns ");
        }
};
}