Merge pull request #37 from ReactionMechanismGenerator/jacobians

Enable manipulation of parameters and reverse mode differentiation

Author: mjohnson541

Manifest.toml

@@ -256,9 +256,9 @@ version = "1.4.0"
 
 [[FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
-git-tree-sha1 = "80becc23db3942d9e5eb8d0d42e725b9c46616f8"
+git-tree-sha1 = "43b397bf66d07d113cc2cc012dd2029d336c3894"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.5.0"
+version = "2.5.1"
 
 [[FixedPointNumbers]]
 git-tree-sha1 = "4aaea64dd0c30ad79037084f8ca2b94348e65eaa"
@@ -537,6 +537,18 @@ git-tree-sha1 = "ea172c86745810136d744fc67650d2e2de669c4f"
 uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 version = "4.4.0"
 
+[[NNPACK_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "c3d1a616362645754b18e12dbba96ec311b0867f"
+uuid = "a6bfbf70-4841-5cb9-aa18-3a8ad3c413ee"
+version = "2018.6.22+0"
+
+[[NNlib]]
+deps = ["Libdl", "LinearAlgebra", "LoopVectorization", "NNPACK_jll", "Pkg", "Requires", "Statistics"]
+git-tree-sha1 = "1d8128735fdf3ab1643dd8bc9499e4b34ccb718d"
+uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+version = "0.7.3"
+
 [[NaNMath]]
 git-tree-sha1 = "c84c576296d0e2fbb3fc134d3e09086b3ea617cd"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
@@ -679,6 +691,12 @@ git-tree-sha1 = "d37400976e98018ee840e0ca4f9d20baa231dc6b"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.0.1"
 
+[[ReverseDiff]]
+deps = ["DiffResults", "DiffRules", "ForwardDiff", "FunctionWrappers", "LinearAlgebra", "MacroTools", "NaNMath", "Random", "SpecialFunctions", "StaticArrays", "Statistics"]
+git-tree-sha1 = "d10f33d434e920442cc6e88fdd5bbc8a5b54494f"
+uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+version = "1.4.1"
+
 [[Roots]]
 deps = ["Printf"]
 git-tree-sha1 = "7bd9d7eee602f0413f24c394386a18cb0d515f36"
@@ -702,9 +720,9 @@ version = "0.8.16"
 
 [[SLEEFPirates]]
 deps = ["Libdl", "SIMDPirates", "VectorizationBase"]
-git-tree-sha1 = "c750d618b7c8268a97e55c70e8c88e56080d30fa"
+git-tree-sha1 = "67ae90a18aa8c22bf159318300e765fbd89ddf6e"
 uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
-version = "0.5.4"
+version = "0.5.5"
 
 [[SentinelArrays]]
 deps = ["Dates", "Random"]
@@ -822,6 +840,12 @@ git-tree-sha1 = "98693daea9bb49aba71eaad6b168b152d2310358"
 uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 version = "0.2.4"
 
+[[Tracker]]
+deps = ["Adapt", "DiffRules", "ForwardDiff", "LinearAlgebra", "MacroTools", "NNlib", "NaNMath", "Printf", "Random", "Requires", "SpecialFunctions", "Statistics", "Test"]
+git-tree-sha1 = "d8173d93a7b4e74d7ab898874f2a4ca04f0cae1a"
+uuid = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
+version = "0.2.8"
+
 [[TreeViews]]
 deps = ["Test"]
 git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"

Project.toml

@@ -17,12 +17,14 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SmoothingSplines = "102930c3-cf33-599f-b3b1-9a29a5acab30"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 
@@ -39,6 +41,7 @@ Parameters = "<1"
 PyCall = "1.91.4"
 PyPlot = "2.9.0"
 QuartzImageIO = "<1"
+ReverseDiff = "<2"
 SmoothingSplines = "<1"
 SpecialFunctions = "<1"
 StaticArrays = "<1"

src/Calculators/Rate.jl

@@ -13,8 +13,8 @@ export AbstractFalloffRate
         Ea::Q
         unc::P = EmptyRateUncertainty()
 end
-@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))::Q
-@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))::Q
+@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
+@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
 export Arrhenius
 
 @with_kw struct PdepArrhenius{T<:Real,Q<:AbstractRateUncertainty,Z<:AbstractRate} <: AbstractRate
@@ -28,10 +28,10 @@ PdepArrhenius(Ps::Array{Q,1},arrs::Array{Z,1}) where {Q<:Real,Z<:AbstractRate} =
     inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
 
     if inds[2] == -1
-        return @inbounds parr.arrs[inds[1]](T=T)::Q
+        return @inbounds parr.arrs[inds[1]](T=T)
     else
-        @inbounds highk = parr.arrs[inds[2]](T=T)::Q
-        @inbounds lowk = parr.arrs[inds[1]](T=T)::Q
+        @inbounds highk = parr.arrs[inds[2]](T=T)
+        @inbounds lowk = parr.arrs[inds[1]](T=T)
         @inbounds Plow = parr.Ps[inds[1]]
         @inbounds Phigh = parr.Ps[inds[2]]
         return @inbounds @fastmath lowk*10^(log10(P/Plow)/log10(Phigh/Plow)*log10(highk/lowk))
@@ -44,10 +44,10 @@ export PdepArrhenius
     unc::Q = EmptyRateUncertainty()
 end
 
-@inline function (marr::MultiArrhenius)(;T::Q=nothing,P::R=0.0,C::S=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (marr::MultiArrhenius)(;T::Q,P::R=0.0,C::S=0.0) where {Q<:Real,R<:Real,S<:Real}
     out = 0.0
     for arr in marr.arrs
-        @fastmath out += arr(T)::Q
+        @fastmath out += arr(T)
     end
     return out
 end
@@ -58,10 +58,10 @@ export MultiArrhenius
     unc::Q = EmptyRateUncertainty()
 end
 
-@inline function (parr::MultiPdepArrhenius)(;T::Q=nothing,P::R=0.0,C::S=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (parr::MultiPdepArrhenius)(;T::Q,P::R=0.0,C::S=0.0) where {Q<:Real,R<:Real,S<:Real}
     out = 0.0
     for pdar in parr.parrs
-        @fastmath out += pdar(T=T,P=P)::Q
+        @fastmath out += pdar(T=T,P=P)
     end
     return out
 end
@@ -73,7 +73,7 @@ export MultiPdepArrhenius
     unc::Q = EmptyRateUncertainty()
 end
 
-(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T)::Q)
+(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
 export ThirdBody
 
 @with_kw struct Lindemann{N<:Integer,K<:AbstractFloat,Q<:AbstractRateUncertainty} <: AbstractFalloffRate
@@ -84,8 +84,8 @@ export ThirdBody
 end
 
 @inline function (lnd::Lindemann)(;T::Q=nothing,P::R=0.0,C::S=nothing) where {Q<:Real,R<:Real,S<:Real}
-    k0 = lnd.arrlow(T=T)::Q
-    kinf = lnd.arrhigh(T=T)::Q
+    k0 = lnd.arrlow(T=T)
+    kinf = lnd.arrhigh(T=T)
     @fastmath Pr = k0*C/kinf
     return @fastmath kinf*Pr/(1.0+Pr)
 end
@@ -102,9 +102,9 @@ export Lindemann
     unc::R = EmptyRateUncertainty()
 end
 
-@inline function (tr::Troe)(;T::Q=nothing,P::R=0.0,C::S=nothing) where {Q<:Real,R<:Real,S<:Real}
-    k0 = tr.arrlow(T=T)::Q
-    kinf = tr.arrhigh(T=T)::Q
+@inline function (tr::Troe)(;T::Q,P::R=0.0,C::S=nothing) where {Q<:Real,R<:Real,S<:Real}
+    k0 = tr.arrlow(T=T)
+    kinf = tr.arrhigh(T=T)
     @fastmath Pr = k0*C/kinf
 
     if tr.T1 == 0.0 && tr.T3 == 0.0
@@ -119,7 +119,7 @@ end
         @fastmath c = -0.4-0.67*log10(Fcent)
         @fastmath F = 10^(log10(Fcent)/(1+((log10(Pr)+c)/(n-d*(log10(Pr))))^2))
     end
-    return @fastmath kinf*(Pr/(1+Pr))*F::S
+    return @fastmath kinf*(Pr/(1+Pr))*F
 end
 export Troe
 
@@ -200,4 +200,4 @@ end
 @inline function getkineticstype(kin::Chebyshev)
     return (string(typeof(kin).name),kin.Tmin,kin.Tmax,kin.Pmin,kin.Pmax,size(kin.coefs)) #different opt functions, but always can do
 end
-export getkineticstype
+export getkineticstype

src/Calculators/Ratevec.jl

@@ -1,4 +1,6 @@
 using Parameters
+using ReverseDiff
+using ForwardDiff
 
 abstract type AbstractRatevec end
 export AbstractRatevec
@@ -19,8 +21,8 @@ function Arrheniusvec(arrs::T) where {T<:AbstractArray}
     end
     return Arrheniusvec(A=A,n=n,Ea=Ea)
 end
-@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))::Array{Q,1}
-@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))::Array{Q,1}
+@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
+@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
 export Arrheniusvec
 
 @with_kw struct Chebyshevvec{T<:AbstractArray,Q<:Real,S<:Real,V<:Real,B<:Real} <: AbstractRate
@@ -85,7 +87,7 @@ end
 export getredpress
 
 @inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0) where {N<:Real,B<:Real,Q<:Real}
-    k = zeros(size(ch.coefs)[1])
+    k = zeros(N,size(ch.coefs)[1])
     Tred = getredtemp(ch,T)
     Pred = getredpress(ch,P)
     klen,Tlen,Plen = size(ch.coefs)

src/Calculators/Thermo.jl

@@ -14,7 +14,7 @@ end
 end
 export NASApolynomial
 
-@inline function calcHSCpdless(poly::NASApolynomial,T::Float64)
+@inline function calcHSCpdless(poly::NASApolynomial,T::X) where {X<:Real}
     if length(poly.coefs) != 7
         Tpoly0 = T
         Tpoly1 = T*T

src/Calculators/Thermovec.jl

@@ -73,7 +73,7 @@ end
 end
 export selectPoly
 
-@inline function calcHSCpdless(poly::NASApolynomialvec,T::Float64)
+@inline function calcHSCpdless(poly::NASApolynomialvec,T::X) where {X<:Real}
     if size(poly.coefs)[1] != 7
         Tpoly0 = T
         Tpoly1 = T*T
@@ -118,7 +118,7 @@ export selectPoly
     return (cpdivR,hdivRT,sdivR)
 end
 
-@inline function calcHSCpdless(nasavec::NASAvec,T::Float64)
+@inline function calcHSCpdless(nasavec::NASAvec,T::X) where {X<:Real}
     poly = selectPoly(nasavec,T)
     return calcHSCpdless(poly,T)
 end