using FiniteStateProjection
using OrdinaryDiffEq

@parameters σ d
rn = @reaction_network begin
    σ, 0 --> A
    d, A --> 0
end σ d

sys = FSPSystem(rn)

# Parameters for our system
ps = [ 10.0, 1.0 ]

# Initial values
u0 = zeros(50)
u0[1] = 1.0

prob = ODEProblem(sys, u0, (0, 10.0), ps)
sol = solve(prob, Vern7(), atol=1e-6)

using FiniteStateProjection
using OrdinaryDiffEq

@parameters ρ σ_on σ_off d
rn = @reaction_network begin
    ρ, G_on --> G_on + M
    (σ_on, σ_off), G_off <--> G_on
    d, M --> 0
end ρ σ_on σ_off d

# This automatically reduces the dimensionality of the
# network by exploiting conservation laws
ih = ReducingIndexHandler(rn)
sys = FSPSystem(rn, ih)

# There is one conserved quantity: G_on + G_off
cons = conservedquantities([1, 0, 0], sys)

# Parameters for our system
ps = [ 15.0, 0.25, 0.15, 1.0 ]

# In the reduced model, G_off = 1 - G_on does not have to be tracked
u0 = zeros(2, 50)
u0[1,1] = 1.0

prob = ODEProblem(sys, u0, (0, 10.0), (ps, cons))
sol = solve(prob, Vern7(), atol=1e-6)