Program Listing for File getSyms.m

Return to documentation for file (matlab/@amifun/getSyms.m)

function [this,model] = getSyms(this,model)
    % getSyms computes the symbolic expression for the requested function
    %
    % Parameters:
    %  model: model definition object @type amimodel
    %
    % Return values:
    %  this: updated function definition object @type amifun
    %  model: updated model definition object @type amimodel

    % store often used variables for ease of notation, dependencies should
    % ensure that these variables are defined

    persistent x p sx w ndw jacw

    nx = model.nx;
    nevent = model.nevent;
    np = model.np;
    nk = model.nk;
    nz = model.nz;
    ny = model.ny;

    fprintf([this.funstr ' | '])
    switch(this.funstr)
        case 'x'
            % create cell array of same size
            xs = cell(nx,1);
            % fill cell array
            for j=1:nx
                xs{j} = sprintf('var_x_%i',j-1);
            end
            % transform into symbolic expression
            this.sym = sym(xs);
            x = this.sym;

        case 'dx'
            % create cell array of same size
            dxs = cell(nx,1);
            % fill cell array
            for j=1:nx
                dxs{j} = sprintf('var_dx_%i',j-1);
            end
            % transform into symbolic expression
            this.sym = sym(dxs);

        case 'p'
            % create cell array of same size
            ps = cell(np,1);
            % fill cell array
            for j=1:np
                ps{j} = sprintf('var_p_%i',j-1);
            end
            % transform into symbolic expression
            this.sym = sym(ps);
            p = this.sym;

        case 'k'
            % create cell array of same size
            ks = cell(nk,1);
            % fill cell array
            for j=1:nk
                ks{j} = sprintf('var_k_%i',j-1);
            end
            % transform into symbolic expression
            this.sym = sym(ks);

        case 'sx'
            % create cell array of same size
            sxs = cell(nx,1);
            % fill cell array
            for j = 1:nx
                sxs{j} = sprintf('var_sx_%i', j-1);
            end
            % transform into symbolic expression
            this.sym = repmat(sym(sxs),[1,np]);
            sx = this.sym;

        case 'sdx'
            % create cell array of same size
            sdx = cell(nx,np);
            % fill cell array
            for j = 1:nx
                for i = 1:np
                    sdx{j,i} = sprintf('var_sdx_%i', j-1);
                end
            end
            % transform into symbolic expression
            this.sym = sym(sdx);

        case 'xB'
            % create cell array of same size
            xBs = cell(nx,1);
            % fill cell array
            for j = 1:nx
                xBs{j} = sprintf('var_xB_%i', j-1);
            end
            % transform into symbolic expression
            this.sym = sym(xBs);

        case 'dxB'
            % create cell array of same size
            dxBs = cell(nx,1);
            % fill cell array
            for j = 1:nx
                dxBs{j} = sprintf('var_dxB_%i', j-1);
            end
            % transform into symbolic expression
            this.sym = sym(dxBs);

        case 'y'
            this.sym = model.sym.y;
            % replace unify symbolic expression
            this = unifySyms(this,model);
            this = makeStrSymsFull(this);


            % activate splines
            for iy = 1:ny
                if(not(all([model.splineflag,model.minflag,model.maxflag])))
                    str = char(this.sym(iy));
                    if(strfind(str,'spline'))
                        model.splineflag = true;
                    end
                    if(strfind(str,'max'))
                        model.maxflag = true;
                    end
                    if(strfind(str,'min'))
                        model.minflag = true;
                    end
                end
            end

        case 'x0'
            this.sym = model.sym.x0;
            % replace unify symbolic expression
            this = unifySyms(this,model);

        case 'dx0'
            this.sym = model.sym.dx0;
            % replace unify symbolic expression
            this = unifySyms(this,model);

        case 'sigma_y'
            this.sym = model.sym.sigma_y;
            this = makeStrSymsFull(this);
            % replace unify symbolic expression
            this = unifySyms(this,model);

        case 'sigma_z'
            this.sym = model.sym.sigma_z;
            this = makeStrSymsFull(this);
            % replace unify symbolic expression
            this = unifySyms(this,model);

        case 'M'
            this.sym = sym(model.sym.M);
            % replace unify symbolic expression
            this = unifySyms(this,model);
            this = makeStrSyms(this);

        case 'xdot'
            this.sym = model.sym.xdot;
            % replace unify symbolic expression
            this = unifySyms(this,model);

            if(strcmp(model.wtype,'iw'))
                if(size(this.sym,2)>size(this.sym,1))
                    this.sym = -transpose(model.fun.M.sym*model.fun.dx.sym)+this.sym;
                else
                    this.sym = -model.fun.M.sym*model.fun.dx.sym+this.sym;
                end
            end

            % create cell array of same size
            xdots = cell(nx,1);
            xdot_olds = cell(nx,1);
            % fill cell array
            for j=1:nx
                xdots{j} = sprintf('xdot_%i',j-1);
                xdot_olds{j} = sprintf('xdot_old_%i',j-1);
            end
            this.strsym = sym(xdots);
            this.strsym_old = sym(xdot_olds);

            % activate splines
            for ix = 1:nx
                if(not(all([model.splineflag,model.minflag,model.maxflag])))
                    str = char(this.sym(ix));
                    if(strfind(str,'spline'))
                        model.splineflag = true;
                    end
                    if(strfind(str,'max'))
                        model.maxflag = true;
                    end
                    if(strfind(str,'min'))
                        model.minflag = true;
                    end
                end
            end

        case 'w'
            optimize = getoptimized(optsym(model.fun.xdot.sym));
            tmpxdot = sym(char(optimize(end)));
            nw = (length(optimize)-1);
            model.nw = nw;
            if(nw>0)
                exprs = arrayfun(@(x) children(x),optimize(1:(end-1)),'UniformOutput',false); % extract symbolic variable
                S.subs = {2};
                S.type='()';
                C = cellfun(@(x) subsref(x,S),exprs,'UniformOutput',false); % get second element
                this.sym = [C{:}]; % transform cell to matrix
                S.subs = {1};
                C = cellfun(@(x) subsref(x,S),exprs,'UniformOutput',false);
                temps = [C{:}];
            end
%             model.nw = 0;
%             nw = 0;
%             this.sym = sym(zeros(0,1));



            ws = cell(nw,1);
            ts = cell(nw,1);
            % fill cell array
            for iw = 1:nw
                ws{iw} = sprintf('var_w_%i', iw-1);
            end
            % transform into symbolic expression
            this.strsym = sym(ws);
            ndw = 0;
            if(nw>0)
                tmpxdot = mysubs(tmpxdot,temps,this.strsym); % replace common expressions
                this.sym = mysubs(this.sym,temps,this.strsym);
                model.updateRHS(tmpxdot); % update rhs
                ndw = 1;
                jacw = jacobian(this.sym,this.strsym);
                vv = sym('v',[length(this.strsym),1]);
                while(sum(jacw^ndw*vv)~=0)
                    ndw = ndw+1;
                end
                ndw = ndw - 1;
            else
                model.updateRHS(tmpxdot); % update RHS anyways to generate sym_noopt for fun.xdot
            end
            w = this.strsym;

        case 'dwdx'
            if(length(model.fun.w.sym)>0)
                jacx = jacobian(model.fun.w.sym,x);
                this.sym = jacx;
                for idw = 1:ndw
                    this.sym = this.sym + (jacw^idw)*jacx; % this part is only to get the right nonzero entries
                end
                % fill cell array
                idx_w = find(this.sym);
                this.strsym = sym(zeros(size(jacx)));
                if(numel(idx_w)>0)
                    for iw = 1:length(idx_w)
                        this.strsym(idx_w(iw)) = sym(sprintf('var_dwdx_%i', iw-1));
                    end
                    model.ndwdx = length(idx_w);
                    % update dwdx with simplified expressions, here we can exploit
                    % the proper ordering of w to ensure correctness of expressions
                    tmp = jacx + jacw*this.strsym;
                    this.sym = tmp(idx_w);
                else
                    this.strsym = sym(zeros(size(jacx)));
                end
            else
                this.sym = sym(zeros(0,nx));
                this.strsym = sym(zeros(0,nx));
            end

        case 'dwdp'
            if(length(model.fun.w.sym)>0)
                jacp = jacobian(model.fun.w.sym,p);
                this.sym = jacp;
                for idw = 1:ndw
                    this.sym = this.sym + (jacw^idw)*jacp; % this part is only to get the right nonzero entries
                end
                % fill cell array
                idx_w = find(this.sym);
                this.strsym = sym(zeros(size(jacp)));
                if(numel(idx_w)>0)
                    for iw = 1:length(idx_w)
                        this.strsym(idx_w(iw)) = sym(sprintf('var_dwdp_%i', iw-1));
                    end
                    model.ndwdp = length(idx_w);
                    % update dwdx with simplified expressions, here we can exploit
                    % the proper ordering of w to ensure correctness of expressions
                    tmp = jacp + jacw*this.strsym;
                    this.sym = tmp(idx_w);
                end
            else
                this.sym = sym(zeros(0,nx));
                this.strsym = sym(zeros(0,nx));
            end

        case 'dfdx'
            this.sym=jacobian(model.fun.xdot.sym,x) + jacobian(model.fun.xdot.sym,w)*model.fun.dwdx.strsym;
            this = makeStrSyms(this);

        case 'J'
            if(strcmp(model.wtype,'iw'))
                syms cj
                this.sym = model.fun.dfdx.sym - cj*model.fun.M.sym;
            else
                if(~isempty(w))
                    this.sym = jacobian(model.fun.xdot.sym,x)  + jacobian(model.fun.xdot.sym,w)*model.fun.dwdx.strsym;
                    this.sym_noopt = jacobian(model.fun.xdot.sym_noopt,x);
                else
                    this.sym = jacobian(model.fun.xdot.sym,x);
                    this.sym_noopt = this.sym;
                end
            end

            this = makeStrSymsSparse(this);



        case 'JDiag'
            this.sym = diag(model.fun.J.sym);
            this = makeStrSyms(this);

        case 'dxdotdp'
            if(~isempty(w))
                this.sym=jacobian(model.fun.xdot.sym,p)  + jacobian(model.fun.xdot.sym,w)*model.fun.dwdp.strsym;
                this.sym_noopt = jacobian(model.fun.xdot.sym_noopt,p);
            else
                this.sym=jacobian(model.fun.xdot.sym,p);
                this.sym_noopt = this.sym;
            end

            %%
            % build short strings for reuse of dxdotdp
            % create cell array of same size
            dxdotdps = cell(nx,1);
            % fill cells with strings
            for ix = 1:nx
                dxdotdps{ix} = sprintf('tmp_dxdotdp_%i',ix-1);
            end
            % create full symbolic array
            this.strsym = sym(dxdotdps);

        case 'sx0'
            this.sym=jacobian(model.fun.x0.sym,p);

        case 'sdx0'
            this.sym=jacobian(model.fun.dx0.sym,p);

        case 'sxdot'
            if(np>0)
                if(strcmp(model.wtype,'iw'))
                    this.sym=model.fun.J.strsym*sx(:,1)-model.fun.M.strsym*model.fun.sdx.sym(:,1)+model.fun.dxdotdp.strsym;
                else
                    this.sym=model.fun.J.strsym*sx(:,1)+model.fun.dxdotdp.strsym;
                end
            else
                this.sym = sym(zeros(size(sx,1),0));
            end

        case 'dydx'
            this.sym=jacobian(model.fun.y.sym,x);
            % create cell array of same sizex
            this.strsym = sym(zeros(ny,nx));
            % fill cell array
            this = makeStrSyms(this);

        case 'dydp'
            this.sym=jacobian(model.fun.y.sym,p);
            % create cell array of same size
            this = makeStrSyms(this);

        case 'xBdot'
            if(strcmp(model.wtype,'iw'))
                syms t
                this.sym = diff(transpose(model.fun.M.sym),t)*model.fun.xB.sym + transpose(model.fun.M.sym)*model.fun.dxB.sym - transpose(model.fun.dfdx.sym)*model.fun.xB.sym;
            else
                this.sym = model.fun.JB.sym * model.fun.xB.sym;
            end

        case 'qBdot'
            % If we do second order adjoints, we have to augment
            if (model.nxtrue < nx)
                this.sym = sym(zeros(model.ng, model.np));
                this.sym(1,:) = -transpose(model.fun.xB.sym(1:model.nxtrue)) * model.fun.dxdotdp.sym(1:model.nxtrue, :);
                for ig = 2 : model.ng
                    this.sym(ig,:) = ...
                        -transpose(model.fun.xB.sym(1:model.nxtrue)) * model.fun.dxdotdp.sym((ig-1)*model.nxtrue+1 : ig*model.nxtrue, :) ...
                        -transpose(model.fun.xB.sym((ig-1)*model.nxtrue+1 : ig*model.nxtrue)) * model.fun.dxdotdp.sym(1:model.nxtrue, :);
                end
            else
                this.sym = -transpose(model.fun.xB.sym)*model.fun.dxdotdp.sym;
            end

        case 'dsigma_ydp'
            this.sym = jacobian(model.fun.sigma_y.sym,p);
            this = makeStrSyms(this);

        case 'dsigma_zdp'
            if(nz>0)
                this.sym = jacobian(model.fun.sigma_z.sym,p);
            else
                this.sym = sym(zeros(model.nz,np));
            end
            this = makeStrSyms(this);

        case 'root'
            if(nevent>0)
                this.sym = transpose([model.event.trigger]);
                this = unifySyms(this,model);
            else
                this.sym = sym(zeros(0,1));
            end
            for ir = 1:nevent
                if(not(all([model.splineflag,model.minflag,model.maxflag])))
                    str = char(this.sym(ir));
                    if(strfind(str,'spline'))
                        model.splineflag = true;
                    end
                    if(strfind(str,'max'))
                        model.maxflag = true;
                    end
                    if(strfind(str,'min'))
                        model.minflag = true;
                    end
                end
            end

        case 'drootdp'
            this.sym = jacobian(model.fun.root.sym,p);

        case 'drootdx'
            this.sym = jacobian(model.fun.root.sym,x);

        case 'drootdt'
            % noopt is important here to get derivatives right
            this.sym = diff(model.fun.root.sym,'t') + model.fun.drootdx.sym*model.fun.xdot.sym_noopt;

        case 'sroot'
            this.sym = model.fun.drootdp.sym + model.fun.drootdx.sym*sx;

        case 'srz'
            if(isfield(model.sym,'rz')) % user defined input or from augmentation
                this.sym = jacobian(model.fun.rz.sym,p) + jacobian(model.fun.rz.sym,x)*sx;
            else
                for iz = 1:length(model.z2event)
                    this.sym(iz,:) = model.fun.sroot.sym(model.z2event(iz),:);
                end
            end

        case 's2root'
            switch(model.o2flag)
                case 1
                    s2x = reshape(sx((model.nxtrue+1):end,:),[model.nxtrue,np,np]);
                    vec = sym(eye(np));
                case 2
                    s2x = reshape(sx((model.nxtrue+1):end,:),[model.nxtrue,np,1]);
                    vec = model.sym.k((end-np+1):end);
            end
            for ievent = 1:nevent

                this.sym(ievent,:,:) = (jacobian(model.fun.sroot.sym(ievent,:),p) + jacobian(model.fun.sroot.sym(ievent,:),x(1:model.nxtrue))*sx(1:model.nxtrue,:) + jacobian(model.fun.sroot.sym(ievent,:),x(1:model.nxtrue))*sx(1:model.nxtrue,:))*vec;
                for ix = 1:model.nxtrue
                    this.sym(ievent,:,:) = this.sym(ievent,:,:) + model.fun.drootdx.sym(ievent,ix)*s2x(ix,:,:);
                end
            end

        case 'dtaudp'
            this.sym = sym(zeros(nevent,np));
            for ievent = 1:nevent
                this.sym(ievent,:) = - model.fun.drootdp.sym(ievent,:)/model.fun.drootdt.sym(ievent);
            end

        case 'dtaudx'
            this.sym = sym(zeros(nevent,nx));
            for ievent = 1:nevent
                this.sym(ievent,:) = - model.fun.drootdx.sym(ievent,:)/model.fun.drootdt.sym(ievent);
            end

        case 'stau'
            this.sym = sym(zeros(nevent,np));
            for ievent = 1:nevent
                this.sym(ievent,:) = - model.fun.sroot.sym(ievent,:)/model.fun.drootdt.sym(ievent);
            end
            % create cell array of same size
            staus = cell(1,np);
            % fill cells
            for j=1:np
                staus{j} = sprintf('var_stau_%i',j-1);
            end
            % transform to symbolic variable
            staus = sym(staus);
            % multiply
            this.strsym = staus;

        case 'deltax'
            if(nevent>0)
                this.sym = [model.event.bolus];
                this = unifySyms(this,model);
            else
                this.sym = sym(zeros(0,1));
            end

        case 'deltaxdot'
            this.sym = model.fun.xdot.strsym-model.fun.xdot.strsym_old;

        case 'ddeltaxdp'
            this.sym = sym(zeros(nx,nevent,np));
            for ievent = 1:nevent
                this.sym(:,ievent,:) = jacobian(model.fun.deltax.sym(:,ievent),p);
            end

        case 'ddeltaxdx'
            this.sym = sym(zeros(nx,nevent,nx));
            if(nx>0)
                for ievent = 1:nevent
                    this.sym(:,ievent,:) = jacobian(model.fun.deltax.sym(:,ievent),x);
                end
            end

        case 'ddeltaxdt'
            this.sym = diff(model.fun.deltax.sym,'t');

        case 'deltasx'

            if(nevent>0)
                for ievent = 1:nevent

                    % dtdp  = (1/drdt)*drdp
                    dtdp = model.fun.stau.strsym; % this 1 here is correct, we explicitely do not want ievent here as the actual stau_tmp will only have dimension np

                    % if we are just non-differentiable and but not
                    % discontinuous we can ignore some of the terms!
                    if(any(logical(model.fun.deltax.sym(:,ievent)~=0)))
                        % dxdp  = dx/dt*dt/dp + dx/dp
                        dxdp = sym(zeros(nx,np));
                        for ix = 1:nx
                            dxdp(ix,:) = model.fun.xdot.sym(ix)*dtdp + sx(ix,:);
                        end

                        this.sym(:,:,ievent) = ...
                            + permute(model.fun.ddeltaxdx.sym(:,ievent,:),[1 3 2])*dxdp ...
                            + model.fun.ddeltaxdt.sym(:,ievent)*dtdp ...
                            + permute(model.fun.ddeltaxdp.sym(:,ievent,:),[1 3 2]);
                    else
                        this.sym(:,:,ievent) = sym(zeros([nx,np]));
                    end
                    if(any(model.event(ievent).hflag))
                        this.sym(model.event(ievent).hflag,:,ievent) = ...
                            this.sym(model.event(ievent).hflag,:,ievent) ...
                            - model.fun.deltaxdot.sym(model.event(ievent).hflag)*dtdp;
                    end
                end
            end

        case 'deltaqB'
            if (model.nxtrue < nx)
                ng_tmp = round(nx / model.nxtrue);
                this.sym = sym(zeros(np*ng_tmp,nevent));
            else
                this.sym = sym(zeros(np,nevent));
            end

            for ievent = 1:nevent
                this.sym(1:np,ievent) = -transpose(model.fun.xB.sym)*squeeze(model.fun.ddeltaxdp.sym(:,ievent,:));
                % This is just a very quick fix. Events in adjoint systems
                % have to be implemented in a way more rigorous way later
                % on... Some day...
            end

        case 'deltaxB'
            this.sym = sym(zeros(nx,nevent));
            for ievent = 1:nevent
                this.sym(:,ievent) = -transpose(squeeze(model.fun.ddeltaxdx.sym(:,ievent,:)))*model.fun.xB.sym;
            end


        case 'z'
            if(nevent>0)
                this.sym = transpose([model.event.z]);
                this = unifySyms(this,model);
            else
                this.sym = sym(zeros(0,1));
            end
            % construct the event identifyer, this is a vector which maps
            % events to outputs z
            model.z2event = zeros(length(this.sym),1);
            iz = 0;
            for ievent = 1:nevent
                for jz = 1:length(model.event(ievent).z)
                    iz = iz+1;
                    model.z2event(iz) = ievent;
                end
            end
            this = makeStrSymsFull(this);

        case 'rz'
            this.sym = sym(zeros(size(model.fun.z.sym)));
            if(isfield(model.sym,'rz'))
                this.sym = model.sym.rz;
            else
                for iz = 1:length(model.z2event)
                    this.sym(iz) = model.fun.root.sym(model.z2event(iz));
                end
            end
            this = unifySyms(this,model);
            this = makeStrSymsFull(this);

        case 'dzdp'
            this.sym = jacobian(model.fun.z.sym,p);

            for iz = 1:nz
                this.sym(iz,:) = this.sym(iz,:) + diff(model.fun.z.sym(iz),sym('t'))*model.fun.dtaudp.sym(model.z2event(iz),:);
            end
            % create cell array of same size
            this = makeStrSyms(this);

        case 'drzdp'
            this.sym = jacobian(model.fun.rz.sym,p);
            this = makeStrSyms(this);

        case 'dzdx'
            this.sym = jacobian(model.fun.z.sym,x);
            for iz = 1:nz
                this.sym(iz,:) = this.sym(iz,:)+ diff(model.fun.z.sym(iz),sym('t'))*model.fun.dtaudx.sym(model.z2event(iz),:);
            end
            this = makeStrSyms(this);

        case 'drzdx'
            this.sym = jacobian(model.fun.rz.sym,x);
            this = makeStrSyms(this);

        case 'dzdt'
            if(nz>0)
                this.sym = diff(model.fun.z.sym,'t')+jacobian(model.fun.z.sym,x(1:model.nxtrue))*model.fun.xdot.sym_noopt(1:model.nxtrue);
            else
                this.sym = sym.empty();
            end

        case 'sz'
            this.sym = sym(zeros(nz,np));
            tmpsym = sym(zeros(nz-model.nztrue,np));
            for iz = 1:nz
                if iz <= model.nztrue
                    this.sym(iz,:) = ...
                        + jacobian(model.fun.z.sym(iz),p) ...
                        + jacobian(model.fun.z.sym(iz),x(1:model.nxtrue))*sx(1:model.nxtrue,:) ...
                        + model.fun.dzdt.sym(iz)*model.fun.stau.sym(model.z2event(iz),:);
                else
                    % I have discovered a truly marvelous proof of these equations, which this margin is too small to
                    % contain.
                    %
                    % we need (dtau/dt)\(ddzdpdt*dtdp + dtdpddzdtdp + ddzdpdp       + dtdp*ddzdtdt*dtdp here
                    % term above:        jacx*sx+jacp   dzdt*st       jacx*sx+jacp    dzdt*st
                    % term in aug z:     drootdt        dzdp          dzdp            drootdt
                    % augmented is always drootdt\dzdp or linear combinations
                    % we cannot correctly compute ddzdpdt, but only ddzdtdp so we omit the drootdt part of jacx*sx-jacp
                    % symmetrise it and add it later on.
                    % also the dzdp part contains second order sensitivities and the drootdt part does not (this leads
                    % to 1:model.nxtrue indexing)


                    if(model.o2flag==1)
                        tmpsym(iz-model.nztrue,:) = jacobian(1/model.fun.drootdt.sym(model.z2event(iz)),p)*model.fun.z.sym(iz)*model.fun.drootdt.sym(model.z2event(iz)) ...
                            + jacobian(1/model.fun.drootdt.sym(model.z2event(iz)),x(1:model.nxtrue))*sx(1:model.nxtrue,:)*model.fun.z.sym(iz)*model.fun.drootdt.sym(model.z2event(iz));
                    else
                       error('sz for directional second order sensis was never implemented and I do not know how to, you are on your own here.');
                    end

                    this.sym(iz,:) = ...
                        + jacobian(model.fun.z.sym(iz)*model.fun.drootdt.sym(model.z2event(iz)),p)/model.fun.drootdt.sym(model.z2event(iz)) ...
                        + jacobian(model.fun.z.sym(iz)*model.fun.drootdt.sym(model.z2event(iz)),x)*sx/model.fun.drootdt.sym(model.z2event(iz)) ...
                        + model.fun.dzdt.sym(iz)*model.fun.stau.sym(model.z2event(iz),:);
                end
            end
            if(model.nz>0)
                if(model.o2flag==1)
                    % THIS IS WHAT YOU NEED TO FIX FOR model.o2flag == 2, transpose does not do the deal :(
                    % we need to pay attention here, some sensitivities are encoded as sx and some as augmented x.
                    % the sx ones are not transposable (as the parameter is encoded in the matrix column) so we need to
                    % convert all of them to augmented x.
                    for ip = 1:np
                        tmpsym(:,ip) = subs(tmpsym(:,ip),sx(1:model.nxtrue,1),x((model.nxtrue*ip+1):(model.nxtrue*(ip+1))));
                    end
                    this.sym(model.nztrue+1:end,:) = this.sym(model.nztrue+1:end,:) + tmpsym + transpose(tmpsym);
                    % you might not believe this, but this matrix should (and hopefully will) actually be symmetric ;)
                end
            end

            % create cell array of same size
            szs = cell(nz,np);
            % fill cell array
            for j = 1:nz
                for i = 1:np
                    szs{j,i} = sprintf('sz_%i', j-1);
                end
            end
            % transform into symbolic expression
            this.strsym = sym(szs);

        case 'JBand'
            %do nothing
        case 'JBandB'
            %do nothing
        case 'JSparse'
            %do nothing

        case 'Jy'
            this.sym = model.sym.Jy;
            % replace unify symbolic expression
            this = unifySyms(this,model);
            tmp = arrayfun(@(iy)sym(sprintf('y_%i',iy-1)),1:ny,'UniformOutput',false);
            this.sym = mysubs(this.sym,[tmp{:}],model.fun.y.strsym);
            tmp = arrayfun(@(iy)sym(sprintf('sigma_y_%i',iy-1)),1:ny,'UniformOutput',false);
            this.sym = mysubs(this.sym,[tmp{:}],model.fun.sigma_y.strsym);
        case 'dJydy'
            this.sym = sym(zeros(model.nytrue, model.ng, model.ny));
            for iyt = 1 : model.nytrue
                this.sym(iyt,:,:) = jacobian(model.fun.Jy.sym(iyt,:),model.fun.y.strsym);
            end
        case 'dJydsigma'
            this.sym = sym(zeros(model.nytrue, model.ng, model.ny));
            for iyt = 1 : model.nytrue
                this.sym(iyt,:,:) = jacobian(model.fun.Jy.sym(iyt,:),model.fun.sigma_y.strsym);
            end
        case 'Jz'
            this.sym = model.sym.Jz;
            this = unifySyms(this,model);
            z = arrayfun(@(iz)sym(sprintf('z_%i',iz-1)),1:nz,'UniformOutput',false);
            this.sym = mysubs(this.sym,[z{:}],model.fun.z.strsym);
            tmp = arrayfun(@(iz)sym(sprintf('sigma_z_%i',iz-1)),1:nz,'UniformOutput',false);
            this.sym = mysubs(this.sym,[tmp{:}],model.fun.sigma_z.strsym);
        case 'Jrz'
            this.sym = model.sym.Jrz;
            this = unifySyms(this,model);
            tmp = arrayfun(@(iz)sym(sprintf('rz_%i',iz-1)),1:nz,'UniformOutput',false);
            this.sym = mysubs(this.sym,[tmp{:}],model.fun.z.strsym);
            tmp = arrayfun(@(iz)sym(sprintf('sigma_z_%i',iz-1)),1:nz,'UniformOutput',false);
            this.sym = mysubs(this.sym,[tmp{:}],model.fun.sigma_z.strsym);
        case 'dJzdz'
            this.sym = sym(zeros(model.nztrue, model.ng, model.nz));
            for iz = 1 : model.nztrue
                this.sym(iz,:,:) = jacobian(model.fun.Jz.sym(iz,:),model.fun.z.strsym);
            end
            this = makeStrSyms(this);
        case 'dJrzdz'
            this.sym = sym(zeros(model.nztrue, model.ng, model.nz));
            for iz = 1 : model.nztrue
                this.sym(iz,:,:) = jacobian(model.fun.Jrz.sym(iz,:),model.fun.rz.strsym);
            end
            this = makeStrSyms(this);
        case 'dJzdsigma'
            this.sym = sym(zeros(model.nztrue, model.ng, model.nz));
            for iz = 1 : model.nztrue
                this.sym(iz,:,:) = jacobian(model.fun.Jz.sym(iz,:),model.fun.sigma_z.strsym);
            end
        case 'dJrzdsigma'
            this.sym = sym(zeros(model.nztrue, model.ng, model.nz));
            for iz = 1 : model.nztrue
                this.sym(iz,:,:) = jacobian(model.fun.Jrz.sym(iz,:),model.fun.sigma_z.strsym);
            end

        otherwise
            error('unknown function name')
    end
end

function this = unifySyms(this,model)
    % unifySyms replaces model specific variable names with general
    % ones which conforms with C naming schemes
    %
    % Parameters:
    %  model: model definition object @type amimodel
    this.sym = mysubs(this.sym,model.sym.x,model.fun.x.sym);
    this.sym = mysubs(this.sym,model.sym.p,model.fun.p.sym);
    this.sym = mysubs(this.sym,model.sym.k,model.fun.k.sym);
end

function this = makeStrSymsSparse(this)
    this.strsym = sym(zeros(size(this.sym)));
    idx = find(this.sym);
    for isym = 1:length(idx)
        this.strsym(idx(isym)) = sym(sprintf([this.cvar '_%i'], isym-1));
    end
end

function this = makeStrSyms(this)
    this.strsym = sym(zeros(size(this.sym)));
    idx = find(this.sym);
    idx = transpose(idx(:));
    for isym = idx
        this.strsym(isym) = sym(sprintf([this.cvar '_%i'], isym-1));
    end
end

function this = makeStrSymsFull(this)
    this.strsym = sym(zeros(size(this.sym)));
    for isym = 1:numel(this.strsym)
        this.strsym(isym) = sym(sprintf([this.cvar '_%i'], isym-1));
    end
end

function out = mysubs(in, old, new)
    % mysubs is a wrapper for ther subs matlab function
    %
    % Parameters:
    %  in: symbolic expression in which to replace @type symbolic
    %  old: expression to be replaced @type symbolic
    %  new: replacement expression @type symbolic
    %
    % Return values:
    %  out: symbolic expression with replacement @type symbolic
    if(~isnumeric(in) && ~isempty(old) && ~isempty(symvar(in)))
        matVer = ver('MATLAB');
        if(str2double(matVer.Version)>=8.1)
            out = subs(in, old(:), new(:));
        else
            out = subs(in, old(:), new(:), 0);
        end
    else
        out = in;
    end
end