ACompEvalCoreFun.m

Go to the documentation of this file.

namespace dscomponents{


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

        PointSets;

        
    private:

        pts = {""};


        jrow = {""};
        jend = {""};

        jself = {""};

        deriv = {""};

        dfxsel = {""};

        T = {""};

        
    protected:

        matrix<double> S = {""};
    public:


        ACompEvalCoreFun(sys) {
            this = this@dscomponents.ACoreFun(sys);
        }


        function fx = evaluateComponentSet(integer nr,colvec<double> x,double t) {
            fx = this.evaluateComponents(this.PointSets[nr],...
                this.jend[nr], this.jrow[nr], this.jself[nr], this.S[nr]*x, t);
        }
        function fx = evaluateComponentSetMulti(integer nr,matrix<double> x,rowvec<double> t,matrix<double> mu) {
            fx = this.evaluateComponentsMulti(this.PointSets[nr],...
                this.jend[nr], this.jrow[nr], this.jself[nr], this.S[nr]*x, t, mu);
        }
        function dfx = evaluateComponentSetGradientsAt(integer nr,colvec<double> x,double t) {
            dfx = this.evaluateComponentGradientsAt(this.PointSets[nr],...
                    this.jend[nr], this.jrow[nr], this.jself[nr],...
                    this.S[nr]*x, t) * this.S[nr];

            /*  the multiplication by S{nr} from the right is due to the
             * inner derivative (chain rule).
             * this is identity if no projection occured.
             * If this function has been projected, the argument passed to
             * it is Vz, so the component set gradient is right-multiplied
             * by the components of V (which are in S in this case, see
             * setPointSet). */
        }
        function J = evaluateJacobianSet(integer nr,colvec<double> x,double t) {

            J = this.evaluateComponentPartialDerivatives(this.PointSets[nr],...
                this.jend[nr], this.jrow[nr], this.deriv[nr], ...
                this.jself[nr], this.S[nr]*x, t, this.dfxsel[nr]);
            /*  The deriv{nr} contains only derivative indices for non-zero
             * jacobian elements (determined in setPointSet using the
             * JSparsityPattern). Thus, the return values of
             * evaluateComponentPartialDerivatives only are of the size of
             * the actually non-zero jacobian values.
             * The T matrix is a sparse transformation, that places those
             * values into the correct spots with all other values that have
             * been demanded to be evaluated but are actually zero.
             *
             * See setPointSet */
            J = this.T[nr] * J;
        }
        function J = evaluateJacobianSetMulti(integer nr,matrix<double> x,rowvec<double> t,colvec<double> mu) {
            singlemu = size(mu,2) == 1;
            if isempty(mu)
                mu = [];
                singlemu = true;
            end
            if singlemu
                this.prepareSimulation(mu);
            end
            if isempty(t)
                t = double.empty(0,size(x,2));
            elseif length(t) == 1
                t = ones(1,size(x,2))*t;
            end
            m = size(x,2);
            if m == 0
                m = length(t);
            end
            J = zeros(length(this.deriv[nr]), m);
            for idx = 1:size(x,2)
                if ~singlemu
                    this.prepareSimulation(mu(:, idx));
                end
                J(:,idx) = this.evaluateComponentPartialDerivatives(this.PointSets[nr],...
                this.jend[nr], this.jrow[nr], this.deriv[nr], ...
                this.jself[nr], this.S[nr]*x(:,idx), t(idx), this.dfxsel[nr]);
            end

            /*  The deriv{nr} contains only derivative indices for non-zero
             * jacobian elements (determined in setPointSet using the
             * JSparsityPattern). Thus, the return values of
             * evaluateComponentPartialDerivatives only are of the size of
             * the actually non-zero jacobian values.
             * The T matrix is a sparse transformation, that places those
             * values into the correct spots with all other values that have
             * been demanded to be evaluated but are actually zero.
             *
             * See setPointSet */
            J = this.T[nr] * J;
        }
        function  setPointSet(nr,pts,jpd) {

            if nr > 4
                warning(" KerMor:Devel ",[" Point set numbering seems to  " ...
                    " be used in a different context. 1-4 are safe to work  "...
                    " with at this stage, make sure nothing gets overridden. "]);
            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

            this.pts[nr] = pts;

            jr = []; js = logical.empty(0,1);
            je = zeros(1,length(pts));
            SP = this.JSparsityPattern;
            if ~isempty(SP)
                /*  Jacobian extras */
                if nargin == 4
                    deri = [];
                    full_mapping = [];
                    requested_len = 0;
                    dfx_sel = sparse(length(pts),0);
                end
                for i=1:length(pts)
                    sprow = SP(pts(i),:);
                    inew = find(sprow);
                    if isempty(inew)
                        inew = 1;
                        warning(" setPointSet: Found all-zero Jacobian row for point %d. Using first state-space dof as 'dependency' ",pts(i));
                    end
                    jr = [jr inew];/* #ok */

                    js = [js inew == pts(i)];/* #ok */

                    je(i) = length(jr);
                    /*  Jacobian extras */
                    if nargin == 4
                        des_der = jpd[i];
                        /*  Select elements of sparsity pattern that have been
                         * requested to be evaluated on evaluateJacobianSet */
                        full_mapping = [full_mapping sprow(des_der)];/* #ok
                         * deri contains the accumulated indices of 
                         * jacobian entries that are nonzero and thus make sense
                         * to evaluate. */

                        [~, dpos, matchidx] = intersect(des_der, inew);
                        /*  Take the match indices, but sort them in the order
                         * the elements occur in jpd{i} to maintain correct
                         * ordering. */
                        [~, sidx] = sort(dpos);
                        pos = reshape(matchidx(sidx),1,[]);

    /*                      pos1 = jpd{i}(deriv_elem); */
                        if ~isempty(pos)
                            if i==1
                                offs = 0;
                            else
                                offs = je(i-1);
                            end    
                            deri = [deri pos+offs];/* #ok */

                        end
                        /*  Sum up the total length for later transformation */
                        requested_len = requested_len + length(des_der);

                        /*  Augment dfx selection matrix (only needed for default
                         * finite differences implementation) */
                        dfx_sel(i,(end+1):(end+length(inew))) = 1;/* #ok */

                    end
                end
            else
                warning(" Empty JSparsityPattern. DEIM approximation will not use any state space arguments. ");
            end
            this.jrow[nr] = jr;
            this.jself[nr] = js;
            this.jend[nr] = je;

            /*  Compose x-entry selection matrix */
            len = je(end);
            sel = jr(1:len);
            if ~isempty(this.V)
                this.S[nr] = this.V(sel,:);
            else
                this.S[nr] = sparse(1:len,sel,ones(len,1),len,this.xDim);
            end

            /*  Extras for jacobian evaluations */
            this.deriv[nr] = [];
            this.T[nr] = [];
            if nargin == 4
                this.deriv[nr] = deri^t;
                /*  Find all the indices of the full mapping record that are
                 * nonzero */
                at = find(full_mapping);
                /*  Create a sparse transformation matrix, that assigns the
                 * possibly nonzero jacobian entries returned by
                 * evaluateComponentPartialDerivatives to their correct spot
                 * according to the initially demanded derivatives (this
                 * inserts zeros everywhere). If no zero jacobian entries
                 * have been specified in setPointSet this is the identity
                 * matrix. */
                this.T[nr] = sparse(at,1:length(at),ones(length(at),1),...
                    requested_len,length(at));
                /*  Store convenience dfx selection matrix */
                this.dfxsel[nr] = logical(dfx_sel);
            end
        }
        function target = project(V,W,target) {
            if nargin < 4
                target = this.clone;
            end
            target = project@dscomponents.ACoreFun(this, V, W, target);
            /*  In case the point sets are not created newly: update the
             * current selection matrices */
            for i=1:length(this.S)
                if ~isempty(this.S[i])
                    target.S[i] = this.S[i]*V;
                end
            end
        }


        function copy = clone(copy) {
            copy = clone@dscomponents.ACoreFun(this, copy);
            copy.pts= this.pts;
            copy.S= this.S;
            copy.jrow= this.jrow;
            copy.jend= this.jend;
            copy.jself= this.jself;
            copy.deriv= this.deriv;
            copy.dfxsel= this.dfxsel;
            copy.T= this.T;
        }


        
#if 0 //mtoc++: 'get.PointSets'
function psets = PointSets() {
            psets = this.pts;
        }

#endif


        function res = test_ComponentEvalMatch(xsize) {

            if ~isempty(this.V)
                error(" Cannot run test for projected ACompEvalCoreFuns. ");
            end

            if nargin < 2
                xsize = 5;
            end
            x = rand(this.xDim,xsize);
            mu = rand(20,xsize); /*  simply assume param dim<=20 */

            t = rand(1,xsize);
            fx = this.evaluateMulti(x, t, mu);
            oldpts = [];
            if ~isempty(this.PointSets)
                oldpts = this.PointSets[1];
            end

            nsets = 3;
            s = RandStream(" mt19937ar "," Seed ",2);
            /*  Limit set sizes to 5000
             *setsizes = s.randi(min(size(fx,1),5000), nsets, 1); @temp */
            setsizes = s.randi(200, nsets, 1);
            sets = cell(nsets,1);
            for i = 1:length(setsizes)
                sets[i] = unique(s.randi(size(fx,1),1,setsizes(i)));
            end
            /*  Add an extra set with full size (only for functions with dim less than 10000) */
            if size(fx,1) <= 10000
                sets[end+1] = 1:size(fx,1);
            end
            res = true;

            for idx=1:length(sets)
                set = sets[idx];
                this.setPointSet(1, set);
                fxc = this.evaluateComponentSetMulti(1, x, t, mu);
                tmp = fx(set,:);
                err = abs((tmp-fxc)./tmp);
                err = err(tmp ~= 0);
                d = max(err);
                lres = isempty(d) || d < 1e-3;
                if ~lres || d > 10*sqrt(eps)
                    fprintf(2," Warning! ACompEvalCoreFun evaluation test: Max rel. diff. for component wise vs. full evaluation: %e\n ",full(d));
                end
                res = res && lres;
            end
            if ~isempty(oldpts)
                this.setPointSet(1,oldpts);
            end
        }
    protected:

        function matrix<double>dfx = evaluateComponentGradientsAt(rowvec<integer> pts,rowvec<integer> ends,rowvec<integer> idx,rowvec<integer> self,colvec<double> x,double t) {
            dt = sqrt(eps);
            d = size(x,1);
            X = repmat(x,1,d); T = repmat(t,1,d); MU = repmat(this.mu,1,d);
            I = speye(d,d)*dt;
            hlp = this.evaluateComponents(pts, ends, idx, self, x, t);
            dfx = (this.evaluateComponentsMulti(pts, ends, idx, self, X+I, T, MU) ...
                - repmat(hlp,1,d))/dt;
        }
        function dfx = evaluateComponentPartialDerivatives(rowvec<integer> pts,rowvec<integer> ends,rowvec<integer> idx,rowvec<integer> deriv,rowvec<integer> self,colvec<double> x,double t,dfxsel) {
            dt = sqrt(eps);
            d = length(deriv);
            xd = size(x,2);

            X = repmat(x,1,d); T = repmat(t,1,d); /* #ok<*PROP> */

            I = sparse(deriv,1:d,ones(1,d),size(x,1),d)*dt;
            dfx = (this.evaluateComponentsMulti(pts, ends, idx, self, X+I, T, this.mu) ...
                - this.evaluateComponentsMulti(pts, ends, idx, self, X, T, this.mu))/dt;
            dfx = dfx(repmat(dfxsel(:,deriv),1,xd));
        }
        function dfx = evaluateComponentPartialDerivativesMulti(pts,ends,idx,deriv,self,colvec<double> x,double t,colvec<double> mu,dfxsel) {
            dt = sqrt(eps);
            d = length(deriv);
            xd = size(x,2);
            el = reshape(repmat(1:xd,d,1),1,[]);
            X = x(:,el); T = t(el); MU = mu(:,el);
            I = repmat(sparse(deriv,1:d,ones(1,d),size(x,1),d)*dt,1,xd);
            dfx = (this.evaluateComponentsMulti(pts, ends, idx, self, X+I, T, MU) ...
                - this.evaluateComponentsMulti(pts, ends, idx, self, X, T, MU))/dt;
            dfx = dfx(repmat(dfxsel(:,deriv),1,xd));
            dfx = reshape(dfx,[],xd);
        }
        function fx = evaluateComponentsMulti(pts,ends,idx,self,colvec<double> x,double t,colvec<double> mu) {
            fx = zeros(length(pts),size(x,2));
            for k=1:size(x,2)
                this.prepareSimulation(mu(:,k));
                fx(:,k) = this.evaluateComponents(pts, ends, idx, self, x(:,k), t(k));
            end
        }
    protected: /* ( Abstract ) */

        virtual function fx = evaluateComponents(rowvec<integer> pts,rowvec<integer> ends,rowvec<integer> idx,rowvec<integer> self,matrix<double> x,rowvec<double> t,colvec<double> mu) = 0;
    protected: /* ( Static ) */

        static function obj = loadobj(obj,from) {
            if nargin == 2
                obj.pts= from.pts;
                obj.S= from.S;
                obj.jrow= from.jrow;
                obj.jend= from.jend;
                obj.jself= from.jself;
                obj.deriv= from.deriv;
                obj.dfxsel= from.dfxsel;
                obj.T= from.T;
                obj = loadobj@dscomponents.ACoreFun(obj, from);
            else
                obj = loadobj@dscomponents.ACoreFun(obj);
            end
        }

};
}