function SymModel = sbml2Symbolic(SimbioModel, opts)
%SBML2SYMBOLIC converts a Simbiology model into a structure of symbolic
%   variables so that differentiation and other mathematical manipulation
%   can be done on the system.
% 
%   SymModel = sbml2Symbolic(SimbioModel, opts)
%
%   Inputs
%       SimbioModel - A Simbiology model object
%       opts        - Optional function options
%           .Verbose [ true | {false} ]
%               Print progress to command window
%           .Order [ 0 | 1 | {2} | 3 ]
%               Determines how deep the derivatives should be taken with
%               respect to x and p. Each level increases the cost
%               exponentially, but may be required for certain
%               applications.
%
%   Outputs
%       SymModel
%       	.XNames    - nX by 1 cell vector of strings containing the
%                        natural names of the species
%           .Knames    - nK by 1 cell vector of strings containing the
%                        natural names of the rate constants
%           .x0        - nX by 1 vector defining the initial conditions
%           .k         - nK by 1 vector defining the values of the rate
%                        constants
%
%           .Xsyms     - nX by 1 vector of symbolics containing new
%                        differentiable names of the species
%           .Ksyms     - nK by 1 vector of symbolics containing new
%                        differentiable names of the rate constants
%           Xsyms and Ksyms are used in the symbolics below
%
%           .Rsyms     - nR by 1 vector of symbolic expressions
%                        representing the rates of each reaction
%           .S         - nX by nR stoichiometry matrix representing how
%                        each reaction changes each species
%           .f         - nX by 1 vector of symbolics containing expressions
%                        for the odes defining the system
%           .dfdx      - nX by nX matrix of symbolics containing
%                        expressions for the partial derivatives on each
%                        ode (rows) with respect to each species (columns)
%           .dfdk      - nX by nK matrix of symbolics containing
%                        expressions for the partial derivatives on each
%                        ode (rows) with respect to each rate constant
%                        (columns)
%           .df2dx2    - nX*nX by nX matrix of symbolics containing
%                        expressions for the double partial derivative on
%                        vectorized dfdx (rows) with respect to the species
%                        (columns)
%           .df2dk2    - nX*nK by nK matrix of symbolics containing
%                        expressions for the double partial derivative on
%                        vectorized dfdk (rows) with respect to the rate
%                        constants (columns)
%           .df2dkdx   - nX*nX by nK matrix of symbolics containing
%                        expressions for the double partial derivative on
%                        vectorized dfdx (rows) with respect to the species
%                        (columns)
%           .df2dxdk   - nX*nK by nX matrix of symbolics containing
%                        expressions for the double partial derivative on
%                        vectorized dfdx (rows) with respect to the rate
%                        constants (columns)
%
%   Limitations:
%   Not all Simbiology features are compatible with this kind of analysis.
%   This function ignores any events, rules, and functions of the model.
%   Reaction.Reversible is ignored
%   Reaction.Active is ignored
%   Unit conversion is completely ignored.
%   All species are assumed to be in dimensions given by DefaultSpeciesDimension.
%   All reactions are assumed to be in dimensions of amount/time.
%   Warning: there may be other features missing as well. Do careful
%   comparisons to simbiology before trusting this code.

% (c) 2010 David R Hagen & Bruce Tidor
% This work is released under the MIT license.

%% Options
% Resolve missing inputs
if nargin < 2
    opts = [];
end

% Copy model object
SimbioModel = copyobj(SimbioModel);

%Options for displaying progress
defaultOpts.Verbose = 0;
defaultOpts.Order   = 2;

opts = mergestruct(defaultOpts, opts);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Part 1: Extracting the Model Variables %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if opts.Verbose; fprintf('Extracting model components...'); end
%% Model name
Name = SimbioModel.Name;

%% Build up the table of compartments
compartments = SimbioModel.Compartments;

nC = length(compartments);
Cnames = cell(nC,1);
Ccap = zeros(nC,1);

%Get compartment names
for i = 1:nC
    Cnames{i}   = compartments(i).Name;
    Ccap(i)     = compartments(i).Capacity;
end

%Cnames will be used to remove the prefixes off species. Cnames and Ccap
%will also be appended to Knames and k, respectively, to represent when the
%capacity is used as a constant.

%% Build up the table of species
species = SimbioModel.Species;

nX = length(species);
Xnames = cell(nX,1);
x0 = zeros(nX,1);
xCap = zeros(nX,1);

%Get species names and compartment volume
for i = 1:nX %Foreach species
    Xnames{i} = species(i).Name;
    x0(i)     = species(i).InitialAmount;
    xCap(i)   = species(i).Parent.Capacity;
end

%% Build up the table of constants
constants = SimbioModel.Parameters;

nK = length(constants); %number of total parameters
nKm = nK; %length of model parameters
Knames = cell(nK,1);
k = zeros(nK,1);

%Get model parameters
for i = 1:nK %Foreach parameter
    Knames{i}   = constants(i).Name;
    k(i,1)      = constants(i).Value;
end

%Get kinetic law parameters
reactions = SimbioModel.Reactions;

nR = length(reactions);

for i = 1:nR
    kineticlaw = reactions(i).KineticLaw; %will be empty if no kinetic law parameters exist
    if ~isempty(kineticlaw)
        constantskl = kineticlaw.Parameters; %will only fetch parameters unique to this kinetic law
        nKkl = length(constantskl);
        for j = 1:nKkl
            nK = nK + 1; %add one more parameter
            Knames{nK,1} = constantskl(j).Name;
            k(nK,1) = constantskl(j).Value;
        end
    end
end

if opts.Verbose; fprintf('done.\n'); end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Part 2: Renaming Everything to Allow for Symbolic Handling %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if opts.Verbose; fprintf('Renaming variables...'); end
%% Convert compartment names to differentiable names
%The names to be used will be "co1x", "co2x" etc...
Csyms = sym(zeros(nC,1));
Cnicestrs = cell(nC,1);

%We must be careful that the names we assign the variables are not already
%names being used by the model. This code checks if a potential name (i.e.
%"compartment1x" already exists. If it does, the code simply skips that
%number and tries the next one until it finds an unused name. Note that the
%"x" on the end of the variable names prevents string replace from seeing
%"species1" in "species10".

CurrentAttempt = 1; %The numerical appendix that will be tried
CurrentCompartment = 1; %Index through the compartments
AttemptIsGood = true; %becomes true when a name already exists

while(CurrentCompartment <= nC)%Until every compartment
    CurrentName = sprintf('co%dx', CurrentAttempt);
    %Check to make sure a compartment does not have our systematic name
    for check = 1:nC
        %compartments except the one being renamed must not have 'compartment#' as name
        if strcmp(CurrentName, Cnames{check}) && i ~= check
            AttemptIsGood = false;
        end
    end
    %Also check species names (you never know)
    for check = 1:nX
        if strcmp(CurrentName, Xnames{check})
            AttemptIsGood = false;
        end
    end
    %Also check parameter names
    for check = 1:nK
        if strcmp(CurrentName, Knames{check})
            AttemptIsGood = false;
        end
    end
    
    if AttemptIsGood %change it and move on
        SimbioModel.Compartments(CurrentCompartment).rename(CurrentName);
        Csyms(CurrentCompartment) = sym(CurrentName);
        Cnicestrs{CurrentCompartment} = CurrentName;
        CurrentCompartment = CurrentCompartment + 1;
        CurrentAttempt = CurrentAttempt + 1;
    else %it failed; try again with different number
        CurrentAttempt = CurrentAttempt + 1;
        AttemptIsGood = true;
    end
end

%% Convert species names to differentiable names
%The names will be "sp1x", "sp2x", etc...
Xsyms = sym(zeros(nX,1));
Xnicestrs = cell(nX,1);

%The important part here is to give every species a unique name, regardless
%of what compartment it is in. This way the compartments can be deleted
%without consequence.

CurrentAttempt = 1; %The numerical appendix that will be tried
CurrentSpecies = 1; %Index through the species
AttemptIsGood = true; %becomes true when a name already exists

while(CurrentSpecies <= nX) %Until every species
    CurrentName = sprintf('sp%dx', CurrentAttempt);
    %Check to make sure a compartment does not have our systematic name
    for check = 1:nC
        if strcmp(CurrentName, Cnames{check})
            AttemptIsGood = false;
        end
    end
    %Check species names
    for check = 1:nX
        if strcmp(CurrentName, Xnames{check}) && i ~= check
            AttemptIsGood = false;
        end
    end
    %Also check parameter names
    for check = 1:nK
        if strcmp(CurrentName, Knames{check})
            AttemptIsGood = false;
        end
    end
    
    if AttemptIsGood %change it and move on
        SimbioModel.Species(CurrentSpecies).rename(CurrentName);
        Xsyms(CurrentSpecies) = sym(CurrentName);
        Xnicestrs{CurrentSpecies} = CurrentName;
        CurrentSpecies = CurrentSpecies + 1;
        CurrentAttempt = CurrentAttempt + 1;
    else %it failed; try again with different number
        CurrentAttempt = CurrentAttempt + 1;
        AttemptIsGood = true;
    end
end

%% Convert parameter names to things that won't confuse the symbolic toolbox
%The names will be "pa1x", "pa2x", etc...
Ksyms = sym(zeros(nK,1));
Knicestrs = cell(nK,1);

CurrentAttempt = 1; %The numerical appendix that will be tried
CurrentParameter = 1; %Index through the name
AttemptIsGood = true; %becomes true when a name already exists

while(CurrentParameter <= nK) %Until every parameter
    CurrentName = sprintf('pa%dx', CurrentAttempt);
    %Check to make sure a compartment does not have our systematic name
    for check = 1:nC
        if strcmp(CurrentName, Cnames{check})
            AttemptIsGood = false;
        end
    end
    %Check species names
    for check = 1:nX
        if strcmp(CurrentName, Xnames{check})
            AttemptIsGood = false;
        end
    end
    %Check parameter names
    for check = 1:nK
        if strcmp(CurrentName, Knames{check}) && i ~= check
            AttemptIsGood = false;
        end
    end
    
    if AttemptIsGood %change it and move on
        if CurrentParameter <= nKm %if it is a model paramter
            SimbioModel.Parameters(CurrentParameter).rename(CurrentName);
            
        else %otherwise cycle through the kinetic law parameters
            %This loop cycles through the reactions, finding how many
            %kinetic law parameters are in each. ind starts off as the
            %number of non-model parameters and is decremented by the
            %number of kinetic law parameters in each reaction until ind is
            %less than the kinetic law parameters defined in a specific
            %reaction. This means that ind now refers to the parameter that
            %needs to be changed.
            ind = CurrentParameter - nKm;
            ParameterNameNotChanged = true;
            CurrentReaction = 1;
            while(ParameterNameNotChanged)
                kineticlaw = reactions(CurrentReaction).KineticLaw; % Will be empty if this reaction does not use a kinetic law
                if ~isempty(kineticlaw)
                    nKkl = length(kineticlaw.Parameters);
                else
                    nKkl = 0;
                end
                if ind <= nKkl
                    kineticlaw.Parameters(ind).rename(CurrentName);
                    ParameterNameNotChanged = false;
                else
                    ind = ind - nKkl;
                    CurrentReaction = CurrentReaction + 1;
                end
            end
        end
        
        Ksyms(CurrentParameter) = sym(CurrentName);
        Knicestrs{CurrentParameter} = CurrentName;
        CurrentParameter = CurrentParameter + 1;
        CurrentAttempt = CurrentAttempt + 1;
    else %it failed; try again with different number
        CurrentAttempt = CurrentAttempt + 1;
        AttemptIsGood = true;
    end
end

if opts.Verbose; fprintf('done.\n'); end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Part 3: Building the Diff Eqs %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if opts.Verbose; fprintf('Building the differential equations...'); end
%% Build up rate string and stoichiometry
reactions = SimbioModel.Reactions;
speciesDimension = SimbioModel.getconfigset.CompileOptions.DefaultSpeciesDimension;

nR = length(reactions);
S = sparse(nX, nR);
Rstrs = cell(nR,1);

%Get each reaction and build stochiometry matrix
for i = 1:nR %Foreach reaction
    
    %Get the reaction rate
    Rstrs{i,1} = reactions(i).ReactionRate;
    
    %Build the stochiometry matrix
    reactants  = reactions(i).Reactants;
    products   = reactions(i).Products;
    stoichio   = reactions(i).Stoichiometry; % = [Reactants, Products] in order

    nReact     = length(reactants);
    nProd      = length(products);

    for j = 1:nReact %Foreach reactant
        name   = reactants(j).Name;
        ind    = find(strcmp(Xnicestrs, name));
        if strcmp(speciesDimension, 'substance')
            % both stoichiometry and species are in amount
            S(ind, i) = S(ind, i) + stoichio(j);
        else % speciesDimension == 'concentration'
            % stoichiometry is in concentration, reactions are in amount
            S(ind, i) = S(ind, i) + stoichio(j) / xCap(ind);
        end
    end
    
    for j = 1:nProd %Foreach product
        name   = products(j).Name;
        ind    = find(strcmp(Xnicestrs, name));
        if strcmp(speciesDimension, 'substance')
            % both stoichiometry and species are in amount
            S(ind, i) = S(ind, i) + stoichio(j+nReact);
        else % speciesDimension == 'concentration'
            % stoichiometry is in concentration, reactions are in amount
            S(ind, i) = S(ind, i) + stoichio(j+nReact) / xCap(ind);
        end
    end
end

%% Process reaction rates so that they are differentiable
%Delete compartment prefixes "compartment1."
%Add each compartment capacity as a constant
for i = 1:nC
    %'.' is a special character in regexprep, use '\.' to really mean '.'
    Rstrs = regexprep(Rstrs, [Cnicestrs{i} '\.'], '');
    
    nK = nK + 1;
    Knames{nK,1} = Cnames{i};
    k(nK,1) = Ccap(i);
    Knicestrs{nK,1} = Cnicestrs{i};
    Ksyms(nK,1) = Csyms(i);
end

%Create symbolic rate variable
Rsyms = sym(Rstrs);

%% Generate the odes and gradient of the model
%Odes of the model for all variables
f = S*Rsyms;

if opts.Order >= 1
%Gradient of r with respect to x
drdx = jacobian(Rsyms, Xsyms);

%Gradient of r with respect to k
drdk = jacobian(Rsyms, Ksyms);

%Gradient of f with respect to x
dfdx = S*drdx;

%Gradient of f with respect to k
dfdk = S*drdk;
end

if opts.Order >= 2
%Gradient of drdx with respect to x
d2rdx2 = jacobian(vec(drdx), Xsyms);

%Gradient of drdk with respect to k
d2rdk2 = jacobian(vec(drdk), Ksyms);

%Gradient of drdx with respect to k
d2rdkdx = jacobian(vec(drdx), Ksyms);

%Gradient of drdk with respect to x
d2rdxdk = jacobian(vec(drdk), Xsyms);

%Gradient of dfdx with respect to x
d2fdx2 = S*reshape(d2rdx2, nR,nX*nX);
d2fdx2 = reshape(d2fdx2, nX*nX,nX);

%Gradient of dfdk with respect to k
d2fdk2 = S*reshape(d2rdk2, nR,nK*nK);
d2fdk2 = reshape(d2fdk2, nX*nK,nK);

%Gradient of dfdx with respect to k
d2fdkdx = S*reshape(d2rdkdx, nR,nX*nK);
d2fdkdx = reshape(d2fdkdx, nX*nX,nK);

%Gradient of dfdk with respect to x
d2fdxdk = S*reshape(d2rdxdk, nR,nK*nX);
d2fdxdk = reshape(d2fdxdk, nX*nK,nX);
end

if opts.Order >= 3
%Gradient of d2rdx2 with respect to x
d3rdx3 = jacobian(vec(d2rdx2), Xsyms);

%Gradient of d2rdx2 with respect to k
d3rdkdx2 = jacobian(vec(d2rdx2), Ksyms);

%Gradient of d2fdx2 with respect to x
d3fdx3 = S*reshape(d3rdx3, nR,nX*nX*nX);
d3fdx3 = reshape(d3fdx3, nX*nX*nX,nX);

%Gradient of d2fdx2 with respect to k
d3fdkdx2 = S*reshape(d3rdkdx2, nR,nX*nX*nK);
d3fdkdx2 = reshape(d3fdkdx2, nX*nX*nX,nK);
end

if opts.Verbose; fprintf('done.\n'); end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Part 4: Build Symbolic Model %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SymModel.Name       = Name;

SymModel.Xnames     = Xnames;   %cell nXx1
SymModel.Knames     = Knames;   %cell nKx1
SymModel.x0         = x0;       %num nXx1
%SymModel.xCap       = xCap;     %num nXx1
SymModel.k          = k;        %num nXx1

SymModel.Xsyms      = Xsyms;    %sym nXx1
SymModel.Ksyms      = Ksyms;    %sym nXx1
SymModel.Rsyms      = Rsyms;    %sym nRx1
SymModel.S          = S;        %num nXxnR

SymModel.f          = f;        %sym nXx1

if opts.Order >= 1
SymModel.dfdx       = dfdx;     %sym nXxnX
SymModel.dfdk       = dfdk;     %sym nXxnK
end

if opts.Order >= 2
SymModel.d2fdx2     = d2fdx2;   %sym nX*nXxnX
SymModel.d2fdk2     = d2fdk2;   %sym nX*nKxnK
SymModel.d2fdkdx    = d2fdkdx;  %sym nX*nXxnK
SymModel.d2fdxdk    = d2fdxdk;  %sym nX*nKxnX
end

if opts.Order >= 3
SymModel.d3fdx3     = d3fdx3;   %sym nX*nX*nXxnX
SymModel.d3fdkdx2   = d3fdkdx2; %sym nX*nX*nXxnK
end