function m = postProcessModel(m)
%
% Finalize the data structure
%    
% 
% This work is licensed under the Creative Commons Attribution-Noncommercial-No Derivative 
% Works 3.0 United States License. To view a copy of this license, visit 
% http://creativecommons.org/licenses/by-nc-nd/3.0/us/ or send a letter to Creative Commons,
% 171 Second Street, Suite 300, San Francisco, California, 94105, USA.
    m.S = sparse(m.nX, m.nR);
    [m.spinds.i_S, m.spinds.j_S, m.spinds.x_S, m.spinds.n_S, m.S] = ...
        vec2sparse(m.spinds.i_S, m.spinds.j_S, m.spinds.x_S, m.spinds.n_S, m.S);
    
    [m.spinds.i_a1, m.spinds.j_a1, m.spinds.x_a1, m.spinds.n_a1, m.a1] = ...
        vec2sparse(m.spinds.i_a1, m.spinds.j_a1, m.spinds.x_a1, m.spinds.n_a1, m.a1);

    [m.spinds.i_b1, m.spinds.j_b1, m.spinds.x_b1, m.spinds.n_b1, m.b1] = ...
        vec2sparse(m.spinds.i_b1, m.spinds.j_b1, m.spinds.x_b1, m.spinds.n_b1, m.b1);

    [m.spinds.i_a2, m.spinds.j_a2, m.spinds.x_a2, m.spinds.n_a2, m.a2] = ...
        vec2sparse(m.spinds.i_a2, m.spinds.j_a2, m.spinds.x_a2, m.spinds.n_a2, m.a2);

    [m.spinds.i_b2, m.spinds.j_b2, m.spinds.x_b2, m.spinds.n_b2, m.b2] = ...
        vec2sparse(m.spinds.i_b2, m.spinds.j_b2, m.spinds.x_b2, m.spinds.n_b2, m.b2);

    [m.spinds.i_a1P, m.spinds.j_a1P, m.spinds.x_a1P, m.spinds.n_a1P, m.a1P] = ...
        vec2sparse(m.spinds.i_a1P, m.spinds.j_a1P, m.spinds.x_a1P, m.spinds.n_a1P, m.a1P);
    
    [m.spinds.i_b1P, m.spinds.j_b1P, m.spinds.x_b1P, m.spinds.n_b1P, m.b1P] = ...
        vec2sparse(m.spinds.i_b1P, m.spinds.j_b1P, m.spinds.x_b1P, m.spinds.n_b1P, m.b1P);

    [m.spinds.i_a2P, m.spinds.j_a2P, m.spinds.x_a2P, m.spinds.n_a2P, m.a2P] = ...
        vec2sparse(m.spinds.i_a2P, m.spinds.j_a2P, m.spinds.x_a2P, m.spinds.n_a2P, m.a2P);
    
    [m.spinds.i_b2P, m.spinds.j_b2P, m.spinds.x_b2P, m.spinds.n_b2P, m.b2P] = ...
        vec2sparse(m.spinds.i_b2P, m.spinds.j_b2P, m.spinds.x_b2P, m.spinds.n_b2P, m.b2P);

    Ix	= speye(m.nX);
    Iu	= speye(m.nU);
    IxkronIx = speye(m.nX*m.nX);
    
    pRkronx = kron(m.p.',  Ix);
    m.a1    = pRkronx*(m.a1P);
    m.a2    = pRkronx*(m.a2P);
    m.k     = pRkronx*(m.kP);
    m.b1    = pRkronx*(m.b1P);
    m.b2    = pRkronx*(m.b2P);


    % Find the columns of A2 that contain at least one entry
    [rows, cols] = find(m.a2P);
    cols         = unique(cols);
        
    ixx = mod(cols-1, m.nX)+1;
    jxx = ceil(cols / m.nX);
    lxx = (ixx-1).*m.nX +jxx;
    
    kxx = jxx;
    jxx = ixx;
    ixx = lxx;


    ixI = cols;
    jxI = mod(ixI-1, m.nX) + 1;
    ixI_i = cols;
    ixI_j = jxI;
    
    kxI = ceil(ixI/m.nX);
    lxI = (jxI-1)*(m.nX*m.nX)+ixI;

    ixI = lxI;
    jxI = kxI;


    iIx = cols;
    jIx = ceil(iIx/m.nX);
    iIx_i = cols;
    iIx_j = jIx;
    
    kIx = mod(iIx-1, m.nX) + 1;
    lIx = ((jIx-1)*(m.nX*m.nX))+iIx;

    iIx = lIx;
    jIx = kIx;
    
    %Preallocate space for frequently used terms
    xkronx  = sparse(m.nX*m.nX, 1);
    xkronI  = sparse(m.nX*m.nX, m.nX);
    Ikronx  = sparse(m.nX*m.nX, m.nX);
    
%     indsXP = kron(0:(m.nX-1), ones(1, m.nX))'*ones(1,m.nP)*m.nX*m.nP+kron(ones(m.nX, 1), reshape(1:m.nX*m.nP, m.nX, m.nP));
%     indsXX = kron(0:(m.nX-1), ones(1, m.nX))'*ones(1,m.nX)*m.nX^2+kron(ones(m.nX, 1), reshape(1:m.nX^2, m.nX, m.nX));
%     indsXPp(indsXP) = 1:numel(indsXP);
%     indsXXp(indsXX) = 1:numel(indsXX);
%     indsShuffle = reshape((ones(m.nX, 1)*(m.nX*(0:m.nP-1))+(1:m.nX)'*ones(1, m.nP))', [], 1);
    
%     % preprocess shuffled A2 matrix with any properly squared species multiplied
%     a2Deriv           = m.a2';
%     a2Deriv           = a2Deriv(indsXX);
% ******* (THIS IS RIGHT)
%     a2Deriv           = 2*a2Deriv;
% ******* END(THIS IS RIGHT)

% ******* (THIS IS WRONG)
%     indsQuad           = find(kron(ones(m.nX, 1), eye(m.nX, m.nX)));
%     a2Deriv(indsQuad) = 2*a2Deriv(indsQuad);
% ******* END(THIS IS WRONG)
    
	m.f         = @f;
    m.dfdx      = @dfdx;
    m.dfdu      = @dfdu;
    m.dfdp      = @dfdp;
    m.d2fdx2    = @d2fdx2;
    m.d2fdp2    = @d2fdp2;
    m.d2fdpdx   = @d2fdpdx;
    m.d2fdxdp   = @d2fdxdp;
    m.update    = @update;
    
    mfields = {'name'
               'nX'
               'nP'
               'nU'
               'nY'
               'nR'
               'xNames'
               'yNames'
               'pNames'
               'uNames'
               'spinds'
               'a1'
               'a2'
               'b1'
               'b2'
               'k'
               'a1P'
               'a2P'
               'b1P'
               'b2P'
               'kP'
               'c'
               'ic'
               'p'
               'S'
               'f'
               'dfdx'
               'dfdp'
               'dfdu'
               'd2fdx2'
               'd2fdp2'
               'd2fdpdx'
               'd2fdxdp'
               'update'
              };
    m = orderfields(m, mfields);
    
    function val = f(t, x, u)
        % f is a column with length nX
        xkronx(ixx) = x(jxx).*x(kxx); % For Raw Speed!
        val = m.a1*x +  m.a2*xkronx + m.b1*u + m.b2*kroncol(x, u) + m.k;
%         val = m.a1*x +  m.a2*kron(x, x) + m.b1*u + m.b2*kron(x, u) + m.k;   % non-RAW
    end
    
    function val = dfdx(t, x, u)
        % dfdx is f by x
        xkronI = sparse(ixI_i, ixI_j, x(jxI), m.nX^2, m.nX);
        Ikronx = sparse(iIx_i, iIx_j, x(jIx), m.nX^2, m.nX);
        val = m.a1 + m.a2 * (xkronI + Ikronx) + m.b2*kron(Ix, u);
%         val = m.a1 + m.a2 * (kron(x, Ix) + kron(Ix, x)) + m.b2*kron(Ix, u); % non-RAW
    end
    
    function val = dfdu(t, x, u)
        % dfdu is f by u
        val = m.b1 + m.b2*kron(x, Iu);
    end

    function val = dfdp(t, x, u)
        % dfdp is f by p
        xkronx(ixx) = x(jxx).*x(kxx);
        val = m.a1P * x + m.a2P * xkronx + m.b1P*u + m.b2P*kroncol(x, u) + m.kP;
%         val = m.a1P * x + m.a2P * kron(x,x) + m.b1P * u + m.b2P * kron(x,u) + m.kP; % non-RAW
        val = reshape(val, m.nX, m.nP);
    end
    
    function val = d2fdx2(t, x, u)
        % d2fdx2 is f*x by x
        val = 2*m.a2*IxkronIx;
        val = reshape(val, m.nX*m.nX, m.nX);
    end

    function val = d2fdp2(t, x, u)
        % d2fdp2 is f*p by p
        val = sparse(m.nX*m.nP, m.nP);
    end

    function val = d2fdxdp(t, x, u)
        % d2fdxdp is f*p by x
        xkronI = sparse(ixI_i, ixI_j, x(jxI), m.nX^2, m.nX);
        Ikronx = sparse(iIx_i, iIx_j, x(jIx), m.nX^2, m.nX);
%         xkronI(ixI) = x(jxI); % For Double Raw Speed!
%         Ikronx(iIx) = x(jIx); % For Double Raw Speed!
        val = m.a1P + m.a2P * (xkronI + Ikronx) + m.b2P*kron(Ix, u);
%         val = m.a1P + m.a2P * (kron(x,Ix) + kron(I,x)) + m.b2P*kron(Ix, u); % non-RAW
    end

    function val = d2fdpdx(t, x, u)
        % d2fdpdx is f*x by p
        xkronI = sparse(ixI_i, ixI_j, x(jxI), m.nX^2, m.nX);
        Ikronx = sparse(iIx_i, iIx_j, x(jIx), m.nX^2, m.nX);
%         xkronI(ixI) = x(jxI); % For Double Raw Speed!
%         Ikronx(iIx) = x(jIx); % For Double Raw Speed!
        val = m.a1P + m.a2P * (xkronI + Ikronx) + m.b2P*kron(Ix, u);
%         val = m.a1P + m.a2P * (kron(x,Ix) + kron(I,x)) + m.b2P*kron(Ix, u); % non-RAW
        val = full(val);
        val = reshape(val, m.nX, m.nP, m.nX); % Create 3D matrix
        val = permute(val, [1,3,2]); % Flip the 2nd and 3rd axes
        val = reshape(val, m.nX*m.nX, m.nP); % Reshape to f*x by p
        val = sparse(val); % Return as sparse matrix
    end

    function mout = update(p, ic)
        mout = m;
        if nargin>=1 && exists(p),  mout.p  = p;  end
        if nargin>=2 && exists(ic), mout.ic = ic; end
        
        %% update system matrices
        if(m.nP > 0)
            pRkronx = kron(p.',  Ix);
            mout.a1    = pRkronx*(m.a1P);
            mout.a2    = pRkronx*(m.a2P);
            mout.k     = pRkronx*(m.kP);
            mout.b1    = pRkronx*(m.b1P);
            mout.b2    = pRkronx*(m.b2P);
        end
        mout = postProcessModel(mout);
    end
end
