Merge pull request #27 from JuliaDiffEq/myb/sub

Use submodules

Author: ChrisRackauckas

src/DiffEqProblemLibrary.jl

@@ -1,59 +1,25 @@
-__precompile__(false)
+__precompile__(true)
 
 module DiffEqProblemLibrary
 
-using DiffEqBase, ParameterizedFunctions, DiffEqOperators, DiffEqBiological,
-      LinearAlgebra
-
-include("ode/ode_premade_problems.jl")
-include("dae_premade_problems.jl")
-include("dde_premade_problems.jl")
-include("sde_premade_problems.jl")
-include("jump_premade_problems.jl")
-
-#ODE Example Problems
-export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear,
-       prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody,
-       prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff,
-       prob_ode_lorenz, prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades,
-       prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego,
-       prob_ode_hires, prob_ode_pollution, prob_ode_filament,
-       SolverDiffEq, filament_prob
-
- #SDE Example Problems
- export prob_sde_wave, prob_sde_linear, prob_sde_cubic, prob_sde_2Dlinear,
-        prob_sde_lorenz, prob_sde_2Dlinear, prob_sde_additive,
-        prob_sde_additivesystem, oval2ModelExample, prob_sde_bistable,
-        prob_sde_bruss, prob_sde_oscilreact
+module ODEProblemLibrary
+importodeproblems() = include(joinpath(@__DIR__, "ode/ode_premade_problems.jl"))
+end # module
 
- #SDE Stratonovich Example Problems
- export prob_sde_linear_stratonovich, prob_sde_2Dlinear_stratonovich
+module DAEProblemLibrary
+importdaeproblems() = include(joinpath(@__DIR__, "dae_premade_problems.jl"))
+end # module
 
- #DAE Example Problems
- export prob_dae_resrob
+module DDEProblemLibrary
+importddeproblems() = include(joinpath(@__DIR__, "dde_premade_problems.jl"))
+end # module
 
-# DDE Example Problems
-# examples with constant delays
-export prob_dde_1delay, prob_dde_1delay_notinplace, prob_dde_1delay_scalar_notinplace,
-       prob_dde_2delays, prob_dde_2delays_notinplace, prob_dde_2delays_scalar_notinplace,
-       prob_dde_1delay_long, prob_dde_1delay_long_notinplace,
-       prob_dde_1delay_long_scalar_notinplace, prob_dde_2delays_long,
-       prob_dde_2delays_long_notinplace, prob_dde_2delays_long_scalar_notinplace,
-       prob_dde_mackey, prob_dde_wheldon, prob_dde_qs,
-# examples with vanishing time dependent delays
-       prob_ddde_neves1, prob_dde_neves_thompson,
-# examples with state dependent delays
-       prob_dde_paul1, prob_dde_paul2, prob_dde_mahaffy1, prob_dde_mahaffy2,
-# examples with vanishing state dependent delays
-       prob_neves2, prob_dde_gatica
+module SDEProblemLibrary
+importsdeproblems() = include(joinpath(@__DIR__, "sde_premade_problems.jl"))
+end # module
 
-# Jump Example Problems
-export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs,
-# examples mixing mass action and constant rate jumps
-       prob_jump_osc_mixed_jumptypes,
-# examples used in published benchmarks / comparisions
-       prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor,
-# examples approximating diffusion by continuous time random walks
-       prob_jump_diffnetwork
+module JumpProblemLibrary
+importjumpproblems() = include(joinpath(@__DIR__, "jump_premade_problems.jl"))
+end # module
 
 end # module

src/dae_premade_problems.jl

@@ -1,3 +1,8 @@
+using DiffEqBase
+
+#DAE Example Problems
+export prob_dae_resrob
+
 ### DAE Problems
 
 f = function (r, yp, y, p, tres)

src/dde_premade_problems.jl

@@ -1,3 +1,20 @@
+using DiffEqBase
+# DDE Example Problems
+# examples with constant delays
+export prob_dde_1delay, prob_dde_1delay_notinplace, prob_dde_1delay_scalar_notinplace,
+       prob_dde_2delays, prob_dde_2delays_notinplace, prob_dde_2delays_scalar_notinplace,
+       prob_dde_1delay_long, prob_dde_1delay_long_notinplace,
+       prob_dde_1delay_long_scalar_notinplace, prob_dde_2delays_long,
+       prob_dde_2delays_long_notinplace, prob_dde_2delays_long_scalar_notinplace,
+       prob_dde_mackey, prob_dde_wheldon, prob_dde_qs,
+# examples with vanishing time dependent delays
+       prob_ddde_neves1, prob_dde_neves_thompson,
+# examples with state dependent delays
+       prob_dde_paul1, prob_dde_paul2, prob_dde_mahaffy1, prob_dde_mahaffy2,
+# examples with vanishing state dependent delays
+       prob_neves2, prob_dde_gatica
+
+
 # DDE examples with analytical solution
 
 ## Single constant delay

src/jump_premade_problems.jl

@@ -1,9 +1,19 @@
-""" 
+using DiffEqBase, DiffEqBiological
+# Jump Example Problems
+export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs,
+# examples mixing mass action and constant rate jumps
+       prob_jump_osc_mixed_jumptypes,
+# examples used in published benchmarks / comparisions
+       prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor,
+# examples approximating diffusion by continuous time random walks
+       prob_jump_diffnetwork
+
+"""
     General structure to hold JumpProblem info. Needed since
-    the JumpProblem constructor requires the algorithm, so we 
+    the JumpProblem constructor requires the algorithm, so we
     don't create the JumpProblem here.
 """
-struct JumpProblemNetwork 
+struct JumpProblemNetwork
     network         # DiffEqBiological network
     rates           # vector of rate constants or nothing
     tstop           # time to end simulation
@@ -28,7 +38,7 @@ Nsims = 8000
 expected_avg = 5.926553750000000e+02
 prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg)
 """
-    DNA negative feedback autoregulatory model. Protein acts as repressor.    
+    DNA negative feedback autoregulatory model. Protein acts as repressor.
 """
 prob_jump_dnarepressor = JumpProblemNetwork(dna_rs, rates, tf, u0, prob, prob_data)
 
@@ -50,7 +60,7 @@ prob_jump_constproduct = JumpProblemNetwork(bd_rs, rates, tf, u0, prob, prob_dat
 
 nonlin_rs = @reaction_network dtype begin
     k1, 2A --> B
-    k2, B --> 2A 
+    k2, B --> 2A
     k3, A + B --> C
     k4, C --> A + B
     k5, 3C --> 3A
@@ -132,12 +142,12 @@ u0[ findfirst(rs.syms, :S3) ] = params[3]
 tf    = 100.
 prob = DiscreteProblem(u0, (0., tf), rates)
 """
-    Multistate model from Gupta and Mendes, 
+    Multistate model from Gupta and Mendes,
     "An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems",
     Computation 2018, 6, 9; doi:10.3390/computation6010009
     Translated from supplementary data file: Models/Multi-state/fixed_multistate.xml
 """
-prob_jump_multistate = JumpProblemNetwork(rs, rates, tf, u0, prob, 
+prob_jump_multistate = JumpProblemNetwork(rs, rates, tf, u0, prob,
                 Dict("specs_to_sym_name" => specs_sym_to_name, "rates_sym_to_idx" => rates_sym_to_idx, "params" => params))
 
 
@@ -162,14 +172,14 @@ genenetwork *= "end"
 rs = eval( parse(genenetwork) )
 u0 = zeros(Int, length(rs.syms))
 for i = 1:(2*N)
-    u0[findfirst(rs.syms, Symbol("G$(i)"))] = 1    
+    u0[findfirst(rs.syms, Symbol("G$(i)"))] = 1
 end
 tf = 2000.0
 prob = DiscreteProblem(u0, (0.0, tf))
 """
-    Twenty-gene model from McCollum et al, 
+    Twenty-gene model from McCollum et al,
     "The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior"
-    Comp. Bio. and Chem., 30, pg. 39-49 (2006). 
+    Comp. Bio. and Chem., 30, pg. 39-49 (2006).
 """
 prob_jump_twentygenes = JumpProblemNetwork(rs, nothing, tf, u0, prob, nothing)
 
@@ -191,11 +201,11 @@ tf = 4000.
 prob = DiscreteProblem(u0, (0.0, tf), rnpar)
 """
     Negative feedback autoregulatory gene expression model. Dimer is the repressor.
-    Taken from Marchetti, Priami and Thanh, 
-    "Simulation Algorithms for Comptuational Systems Biology", 
+    Taken from Marchetti, Priami and Thanh,
+    "Simulation Algorithms for Comptuational Systems Biology",
     Springer (2017).
 """
-prob_jump_dnadimer_repressor = JumpProblemNetwork(rn, rnpar, tf, u0, prob, 
+prob_jump_dnadimer_repressor = JumpProblemNetwork(rn, rnpar, tf, u0, prob,
                                                 Dict("specs_names" => varlabels))
 
 
@@ -217,7 +227,7 @@ end
 tf = 10.
 """
     Continuous time random walk (i.e. diffusion approximation) example.
-    Here the network in the JumpProblemNetwork is a function that returns a 
+    Here the network in the JumpProblemNetwork is a function that returns a
     network given the number of lattice sites.
     u0 is a similar function that returns the initial condition vector.
 """

src/ode/brusselator_prob.jl

@@ -13,7 +13,7 @@ function brusselator_2d_loop(du, u, p, t)
   @inbounds begin
     A, B, α, xyd, dx, N = p
     α = α/dx^2
-    for I in CartesianRange((N, N))
+    for I in CartesianIndices((N, N))
       x = xyd[I[1]]
       y = xyd[I[2]]
       i = I[1]
@@ -25,7 +25,7 @@ function brusselator_2d_loop(du, u, p, t)
       du[i,j,1] = α*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
       B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
     end
-    for I in CartesianRange((N, N))
+    for I in CartesianIndices((N, N))
       i = I[1]
       j = I[2]
       ip1 = limit(i+1, N)
@@ -40,17 +40,17 @@ end
 function init_brusselator_2d(xyd)
   N = length(xyd)
   u = zeros(N, N, 2)
-  for I in CartesianRange((N, N))
+  for I in CartesianIndices((N, N))
     x = xyd[I[1]]
     y = xyd[I[2]]
     u[I,1] = 22*(y*(1-y))^(3/2)
     u[I,2] = 27*(x*(1-x))^(3/2)
   end
   u
 end
-xyd_brusselator = linspace(0,1,32)
+xyd_brusselator = range(0,stop=1,length=32)
 
-doc"""
+@doc doc"""
 2D Brusselator
 
 ```math
@@ -120,15 +120,15 @@ function brusselator_1d(du, u_, p, t)
 end
 function init_brusselator_1d(N)
   u = zeros(N, 2)
-  x = linspace(0, 1, N)
+  x = range(0, stop=1, length=N)
   for i in 1:N
     u[i, 1] = 1 + sin(2pi*x[i])
     u[i, 2] = 3.
   end
   u
 end
 
-doc"""
+@doc doc"""
 1D Brusselator
 
 ```math

src/ode/filament_prob.jl

@@ -35,14 +35,14 @@ struct NoHydroProjectionCache <: AbstractInextensibilityCache
     J         :: Matrix{T}
     P         :: Matrix{T}
     J_JT      :: Matrix{T}
-    J_JT_LDLT :: Base.LinAlg.LDLt{T, SymTridiagonal{T}}
+    J_JT_LDLT :: LinearAlgebra.LDLt{T, SymTridiagonal{T}}
     P0        :: Matrix{T}
 
     NoHydroProjectionCache(N::Int) = new(
         zeros(N, 3*(N+1)),          # J
         zeros(3*(N+1), 3*(N+1)),    # P
         zeros(N,N),                 # J_JT
-        Base.LinAlg.LDLt{T,SymTridiagonal{T}}(SymTridiagonal(zeros(N), zeros(N-1))),
+        LinearAlgebra.LDLt{T,SymTridiagonal{T}}(SymTridiagonal(zeros(N), zeros(N-1))),
         zeros(N, 3*(N+1))
     )
 end
@@ -66,7 +66,7 @@ function FilamentCache(N=20; Cm=32, ω=200, Solver=SolverDiffEq)
 end
 function stiffness_matrix!(f::AbstractFilamentCache)
     N, μ, A = f.N, f.μ, f.A
-    A[:] = eye(3*(N+1))
+    A[:] = Matrix{Float64}(I, 3*(N+1), 3*(N+1))
     for i in 1 : 3
         A[i,i] =    1
         A[i,3+i] = -2
@@ -94,7 +94,7 @@ function stiffness_matrix!(f::AbstractFilamentCache)
             A[3*j+i,3*(j+2)+i] =  1
         end
     end
-    scale!(A, -μ^4)
+    rmul!(A, -μ^4)
     nothing
 end
 function update_separate_coordinates!(f::AbstractFilamentCache, r)
@@ -126,7 +126,7 @@ end
 function initialize!(initial_conf_type::Symbol, f::AbstractFilamentCache)
     N, x, y, z = f.N, f.x, f.y, f.z
     if initial_conf_type == :StraightX
-        x[:] = linspace(0, 1, N+1)
+        x[:] = range(0, stop=1, length=N+1)
         y[:] = 0 .* x
         z[:] = 0 .* x
     else
@@ -192,14 +192,14 @@ function projection!(f::FilamentCache)
 end
 
 function subtract_from_identity!(A)
-    scale!(-1, A)
+    rmul!(-1, A)
     @inbounds for i in 1 : size(A,1)
         A[i,i] += 1
     end
     nothing
 end
 
-function LDLt_inplace!{T<:Real}(L::Base.LinAlg.LDLt{T,SymTridiagonal{T}}, A::Matrix{T})
+function LDLt_inplace!(L::LinearAlgebra.LDLt{T,SymTridiagonal{T}}, A::Matrix{T}) where {T<:Real}
     n = size(A,1)
     dv, ev = L.data.dv, L.data.ev
     @inbounds for (i,d) in enumerate(diagind(A))

src/ode/ode_premade_problems.jl

@@ -1,5 +1,17 @@
+using DiffEqBase, ParameterizedFunctions, DiffEqOperators, Random, LinearAlgebra,
+      Markdown
+
 srand(100)
 
+#ODE Example Problems
+export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear,
+       prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody,
+       prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff,
+       prob_ode_lorenz, prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades,
+       prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego,
+       prob_ode_hires, prob_ode_pollution, prob_ode_filament,
+       SolverDiffEq, filament_prob
+
 include("ode_linear_prob.jl")
 include("ode_simple_nonlinear_prob.jl")
 include("brusselator_prob.jl")

src/ode/ode_simple_nonlinear_prob.jl

@@ -175,7 +175,7 @@ pleiades = (du,u,p,t) -> begin
   w = view(u,22:28) # y′
   du[1:7] .= v
   du[8:14].= w
-  for i in 14:21
+  for i in 14:28
     du[i] = zero(eltype(u))
   end
   for i=1:7,j=1:7
@@ -241,6 +241,7 @@ Usually solved from 0 to 3.
 prob_ode_pleiades = ODEProblem(pleiades,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],(0.0,3.0))
 
 
+
 srand(100)
 const mm_A = rand(4,4)
 mm_linear = function (du,u,p,t)

src/ode/pollution_prob.jl

@@ -188,7 +188,7 @@ u0[8]  = 0.3
 u0[9]  = 0.01
 u0[17] = 0.007
 
-doc"""
+@doc doc"""
 [Pollution Problem](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/Pollution.ipynb) (Stiff)
 
 This IVP is a stiff system of 20 non-linear Ordinary Differential Equations. It is in the form of ``\\frac{dy}{dt}=f(y), \\quad y(0)=y0,`` with