Initial commit. Full solver in solve_uSEIR function.

LICENSE

+"THE BEER-WARE LICENSE":
+P. Hernandez, C. Pena, J.J. Gomez Cadenas wrote these files. As long
+as you retain this notice you can do whatever you want with this
+stuff. If we meet some day, and you think this stuff is worth it, you
+can buy us a beer in return.
+

Manifest.toml

+# This file is machine-generated - editing it directly is not advised
+
+[[Arpack]]
+deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
+git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
+uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+version = "0.4.0"
+
+[[Arpack_jll]]
+deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
+git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
+uuid = "68821587-b530-5797-8361-c406ea357684"
+version = "3.5.0+3"
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[CodeTracking]]
+deps = ["InteractiveUtils", "UUIDs"]
+git-tree-sha1 = "cab4da992adc0a64f63fa30d2db2fd8bec40cab4"
+uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
+version = "0.5.11"
+
+[[CompilerSupportLibraries_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612"
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "0.3.3+0"
+
+[[DataAPI]]
+git-tree-sha1 = "176e23402d80e7743fc26c19c681bfb11246af32"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.3.0"
+
+[[DataStructures]]
+deps = ["InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "6166ecfaf2b8bbf2b68d791bc1d54501f345d314"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.17.15"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Distributions]]
+deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
+git-tree-sha1 = "c4ed10355637fcb0725dc6a27060f74df24f13cd"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.23.2"
+
+[[FileWatching]]
+uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
+
+[[FillArrays]]
+deps = ["LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "6c89d5b673e59b8173c546c84127e5f623d865f6"
+uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
+version = "0.8.9"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[JuliaInterpreter]]
+deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
+git-tree-sha1 = "12289ab2f9ff307a603a42b20ed83e20717cc475"
+uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
+version = "0.7.15"
+
+[[LibGit2]]
+deps = ["Printf"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[LoweredCodeUtils]]
+deps = ["JuliaInterpreter"]
+git-tree-sha1 = "695206e7ec13cbec6e6ee4d19af7464e80f4d1ad"
+uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
+version = "0.4.4"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "de0a5ce9e5289f27df672ffabef4d1e5861247d5"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "0.4.3"
+
+[[OpenBLAS_jll]]
+deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
+git-tree-sha1 = "1887096f6897306a4662f7c5af936da7d5d1a062"
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.9+4"
+
+[[OpenSpecFun_jll]]
+deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
+git-tree-sha1 = "d51c416559217d974a1113522d5919235ae67a87"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.3+3"
+
+[[OrderedCollections]]
+git-tree-sha1 = "12ce190210d278e12644bcadf5b21cbdcf225cd3"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.2.0"
+
+[[PDMats]]
+deps = ["Arpack", "LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"]
+git-tree-sha1 = "2fc6f50ddd959e462f0a2dbc802ddf2a539c6e35"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.9.12"
+
+[[Pkg]]
+deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[QuadGK]]
+deps = ["DataStructures", "LinearAlgebra"]
+git-tree-sha1 = "dc84e810393cfc6294248c9032a9cdacc14a3db4"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.3.1"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Revise]]
+deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "UUIDs", "Unicode"]
+git-tree-sha1 = "3185d2ee31756af9e20ce045ddfaedcd0df9e4aa"
+uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
+version = "2.6.6"
+
+[[Rmath]]
+deps = ["Random", "Rmath_jll"]
+git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.6.1"
+
+[[Rmath_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "1660f8fefbf5ab9c67560513131d4e933012fc4b"
+uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
+version = "0.2.2+0"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SortingAlgorithms]]
+deps = ["DataStructures", "Random", "Test"]
+git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "0.3.1"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[SpecialFunctions]]
+deps = ["OpenSpecFun_jll"]
+git-tree-sha1 = "a69e1eaf3397fbb5a4d718538463e5cc20ef42a4"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "0.10.2"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
+git-tree-sha1 = "a6102b1f364befdb05746f386b67c6b7e3262c45"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.33.0"
+
+[[StatsFuns]]
+deps = ["Rmath", "SpecialFunctions"]
+git-tree-sha1 = "f290ddd5fdedeadd10e961eb3f4d3340f09d030a"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "0.9.4"
+
+[[SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
+[[Test]]
+deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

Project.toml

+name = "uSEIR"
+uuid = "56d60da9-7444-4793-9579-959e59fcd452"
+authors = ["Alberto Ramos <alberto.ramos@desy.de>"]
+version = "0.1.0"
+
+[deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

README.md

+
+# uSEIR.jl
+
+uSEIR.jl is a package to solve the Non-Markovian SEIR model with 
+arbitrary distributions for exposed/infected.
+
+## Install
+
+The easiest way is to use the package manager
+```julia
+julia> import Pkg
+(v1.1) pkg> add https://gitlab.ift.uam-csic.es/alberto/uSEIR.jl
+```
+
+## Minimal example
+
+Compare the Gamma and delta distribution results for means of
+exposed/infected 4.0 and 6.0 respectively with a canonical R0 of 2.0. 
+
+```julia
+julia> using Random, Distributions, uSEIR, PyPlot
+julia> sol_g = solve_uSEIR(Gamma(4.0,1.0), Gamma(6.0, 1.0), 0.001, 2.0/6.0)
+julia> sol_d = solve_uSEIR(4.0, 6.0, 0.001, 2.0/6.0)
+julia> ioff()
+julia> fig = figure("Comparison delta vs. Gamma",figsize=(8,8))
+julia> ax = gca()
+julia> plot(1:4:length(sol_g), sol_g[3:4:end], 1:4:length(sol_d), sol_d[3:4:end])
+```
+
+## License
+
+"THE BEER-WARE LICENSE":
+P. Hernandez, C. Pena, J.J. Gomez Cadenas wrote these files. As long
+as you retain this notice you can do whatever you want with this
+stuff. If we meet some day, and you think this stuff is worth it, you
+can buy us a beer in return.
+
+

src/uSEIR.jl

+module uSEIR
+
+import Random, Distributions
+
+# Include types
+include("uSEIRsolvers.jl") 
+
+# Methods
+export solve_uSEIR
+
+end

src/uSEIRsolvers.jl

+"""
+Create a probability table 
+
+    pdf_table(d::UnivariateDistribution, eps::Float, tol::Float64)
+
+Given a distribution `d`, the interval `[0,quantile(d, tol)]` is 
+divided in `n` sub-intervals of length `eps`. The probabilities 
+P_n for `x` falling in each subinterval are returned as a vector.
+
+## Arguments
+- d: An `UnivariateDistribution`. See `Distributions.jl`.
+- eps: Discretization step size.
+- tol: The quantile that marks the end of the interval.
+
+## Examples
+```julia-repl
+julia> pdf_t = pdf_table(Normal(10.0, 0.1), 0.01, 0.9999)
+```
+"""
+function pdf_table(d::Distributions.UnivariateDistribution, eps::Float64, tol::Float64)
+
+    n::Int64 = round(Int64, Distributions.quantile(d, tol)/eps)
+
+    pdf = zeros(Float64, n)
+    cd1 = 0.0
+    for i in 1:n-1
+        cd2 = Distributions.cdf(d, i*eps)
+        pdf[i] = cd2-cd1
+        cd1 = cd2
+    end
+    pdf[end] = 1.0-sum(pdf)
+
+    return pdf
+end
+
+
+"""
+Solve uSEIR: arbitrary distributions for infected anf exposed times
+
+    solve_seir(de, di, eps, rate, dtr = 1, i0 = 1.0E-6, stp = 1.0E-6)
+
+Given the distributions for the time that people are exposed and infected, this 
+routines solver the SEIR model.
+
+## Arguments
+- de, di: `UnivariateDistribution`. The distributions for the time that people 
+  are exposed and infected respectively. If a `Float64` is introduced the program 
+  interpret this as the distribution being a delta function.
+- eps: `Float64`. The discretization parameter.
+- rate: `Float64`. The infection rate.
+- dtr: `Int64`. Optional (default 1). How often one wants the values of 
+  `s(t), e(t), i(t), r(t)` be returned
+- i0: `Float64`. Optional (default 1.0E-6). Ratio of people infected at `t=0`.
+- stp: `Float64`. Optional (default 1.0E-10). Stoping criteria: Simulation ends when 
+  e(t) and i(t) are both smaller than `stp`.
+
+## Returns
+- solution: `Vector{Float64}`. The values of `s(t), e(t), i(t), r(t)` for each time. 
+
+## Examples
+```julia-repl
+julia> slv = solve_uSEIR(Gamma(5.0, 1.0), Gamma(7.0, 1.4), 0.01, 2.0/7.0)
+julia> # Get the values of the infected
+julia> infected = slv(3:4:end)
+julia> # Get the values of the exposed
+julia> exposed  = slv(2:4:end)
+julia> # Get the values of the recovered
+julia> inf = slv(4:4:end)
+```
+"""
+function solve_uSEIR(de::Distributions.UnivariateDistribution, di::Distributions.UnivariateDistribution, eps::Float64, prob::Float64, dtr::Int64=1, i0::Float64 = 1.0E-6, stp::Float64 = 1.0E-10)
+    
+    # Build probability table
+    pn  = prob*eps
+    pdi = pdf_table(di, eps, 0.9999)
+    ni  = length(pdi)
+    pde = pdf_table(de, eps, 0.9999)
+    ne  = length(pde)
+
+    # set initial condition
+    S    = 1.0-i0
+    E    = zeros(ne)
+    E[1] = i0
+    I    = zeros(ni)
+    R    = 0.0
+
+    sol = Vector{Float64}()
+
+    n  = 0
+    sI = 0.0
+    while true
+        n = n + 1
+        R = R + I[1]
+        for k in 1:ni-1
+	    I[k] = I[k+1] + pdi[k]*E[1]
+        end
+        I[end] = pdi[end]*E[1]
+
+        for k in 1:ne-1
+	    E[k] = E[k+1] + pn*pde[k]*sI*S
+        end
+        E[end] = pn*pde[end]*sI*S
+
+        S = S - pn*sI*S
+
+        sI = sum(I)
+        sE = sum(E)
+        
+        if (mod(n-1,dtr) == 0)
+            append!(sol, [S, sE, sI, R])
+        end
+
+        if ( (sE < stp) && (sI < stp) ) 
+	    break
+        end
+
+    end
+    
+    return sol
+end
+
+function solve_uSEIR(te::Float64, ti::Float64, eps::Float64, prob::Float64, dtr::Int64=1, i0::Float64 = 1.0E-6, stp::Float64 = 1.0E-10)
+    
+    pn = prob*eps
+    ni = round(Int64, ti/eps)
+    ne = round(Int64, te/eps)
+    
+    # set initial condition
+    S = 1.0-i0
+    E = zeros(ne)
+    E[1] = i0
+    I = zeros(ni)
+    R = 0.0
+
+    sol = Vector{Float64}()
+
+    n  = 0
+    sI = 0.0
+    while true
+        n = n + 1
+        R = R + I[1]
+        for k in 1:ni-1
+	    I[k] = I[k+1]
+        end
+        I[end] = E[1]
+
+        for k in 1:ne-1
+	    E[k] = E[k+1]
+        end
+        E[end] = pn*sI*S
+
+        S = S - pn*sI*S
+
+        sI = sum(I)
+        sE = sum(E)
+        
+        if (mod(n-1,dtr) == 0)
+            append!(sol, [S, sE, sI, R])
+        end
+
+        if ( (sE < stp) && (sI < stp) ) 
+	    break
+        end
+
+    end
+    
+    return sol
+end
+
+function solve_uSEIR(de::Distributions.UnivariateDistribution, ti::Float64, eps::Float64, prob::Float64, dtr::Int64=1, i0::Float64 = 1.0E-6, stp::Float64 = 1.0E-10)
+    
+    pn  = prob*eps
+    ni  = round(Int64, ti/eps)
+    pde = pdf_table(de, eps, 0.9999)
+    ne  = length(pde)
+    
+    # set initial condition
+    S = 1.0-i0
+    E = zeros(ne)
+    E[1] = i0
+    I = zeros(ni)
+    R = 0.0
+
+    sol = Vector{Float64}()
+
+    n  = 0
+    sI = 0.0
+    while true
+        n = n + 1
+        R = R + I[1]
+        for k in 1:ni-1
+	    I[k] = I[k+1]
+        end
+        I[end] = E[1]
+
+        for k in 1:ne-1
+	    E[k] = E[k+1] + pn*pde[k]*sI*S
+        end
+        E[end] = pn*pde[end]*sI*S
+
+        S = S - pn*sI*S
+
+        sI = sum(I)
+        sE = sum(E)
+        
+        if (mod(n-1,dtr) == 0)
+            append!(sol, [S, sE, sI, R])
+        end
+
+        if ( (sE < stp) && (sI < stp) ) 
+	    break
+        end
+
+    end
+    
+    return sol
+end
+
+function solve_uSEIR(te::Float64, di::Distributions.UnivariateDistribution, eps::Float64, prob::Float64, dtr::Int64=1, i0::Float64 = 1.0E-6, stp::Float64 = 1.0E-10)
+    
+    pn = prob*eps
+    pdi = pdf_table(di, eps, 0.9999)
+    ni  = length(pdi)
+    ne = round(Int64, te/eps)
+    
+    # set initial condition
+    S = 1.0-i0
+    E = zeros(ne)
+    E[1] = i0
+    I = zeros(ni)
+    R = 0.0
+
+    sol = Vector{Float64}()
+
+    n  = 0
+    sI = 0.0
+    while true
+        n = n + 1
+        R = R + I[1]
+        for k in 1:ni-1
+	    I[k] = I[k+1] + pdi[k]*E[1]
+        end
+        I[end] = pdi[end]*E[1]
+
+        for k in 1:ne-1
+	    E[k] = E[k+1]
+        end
+        E[end] = pn*sI*S
+
+        S = S - pn*sI*S
+
+        sI = sum(I)
+        sE = sum(E)
+        
+        if (mod(n-1,dtr) == 0)
+            append!(sol, [S, sE, sI, R])
+        end
+
+        if ( (sE < stp) && (sI < stp) ) 
+	    break
+        end
+
+    end
+    
+    return sol
+end