move from Dual and HyperDualNumbers to ForwardDiff (#8)

move from Dual and HyperDualNumbers to ForwardDiff

* add windows in travis

* remove FowardDiff from extras

* travis to allow failure

* cleaned up

Author: briochemc

.appveyor.yml

@@ -1,2 +1,3 @@
 Manifest.toml
 docs/build/
+*html

.gitignore

@@ -3,32 +3,24 @@ language: julia
 os:
   - osx
   - linux
+  - windows
 
 julia:
-  - 1.0
+  - 1
+  - 1.3
   - nightly
 
-# Uncomment the following lines to allow failures on nightly julia
-# (tests will run but not make your overall status red)
-matrix:
-  allow_failures:
-  - julia: nightly
-
 notifications:
   email: false
 
-#script: # the default script is equivalent to the following
-#  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
-#  - julia -e 'Pkg.clone(pwd()); Pkg.build("F1Method"); Pkg.test("F1Method"; coverage=true)';
-
-after_success:
-  - julia -e 'using Pkg; cd(Pkg.dir("F1Method")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())';
-  - julia -e 'using Pkg; cd(Pkg.dir("F1Method")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())';
+codecov: true
 
 jobs:
+  allow_failures:
+    - julia: nightly
   include:
     - stage: "Documentation"
-      julia: 1.0
+      julia: 1.3
       os: linux
       script:
         - julia --project=docs/ -e 'using Pkg; Pkg.instantiate();

.travis.yml

@@ -1,20 +1,20 @@
 name = "F1Method"
 uuid = "d5605bae-6d76-11e9-2fd7-71bcf42edbaa"
 authors = ["Benoit Pasquier <pasquieb@uci.edu>"]
-version = "0.2.1"
+version = "0.3"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
-HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
+DiffEqBase = "6"
+ForwardDiff = "0.10"
 julia = "1.0"
 
 [extras]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ForwardDiff"]
+test = ["Test"]

Project.toml

@@ -3,36 +3,27 @@
 
 # F-1 algorithm
 
-<p>
-  <a href="https://doi.org/10.5281/zenodo.2667835">
-    <img src="https://zenodo.org/badge/DOI/10.5281/zenodo.2667835.svg" alt="DOI">
-  </a>
-  <a href="https://github.com/briochemc/F1Method.jl/blob/master/LICENSE">
-    <img alt="License: MIT" src="https://img.shields.io/badge/License-MIT-yellow.svg">
-  </a>
-</p>
 <p>
   <a href="https://briochemc.github.io/F1Method.jl/stable/">
-    <img src=https://img.shields.io/badge/docs-stable-blue.svg>
-  </a>
-  <a href="https://briochemc.github.io/F1Method.jl/latest/">
-    <img src=https://img.shields.io/badge/docs-dev-blue.svg>
+    <img src=https://img.shields.io/badge/docs-stable-important.svg?style=flat-square&label=Documentation&logo=Read%20the%20Docs>
   </a>
 </p>
+
 <p>
-  <a href="https://travis-ci.com/briochemc/F1Method.jl">
-    <img alt="Build Status" src="https://travis-ci.com/briochemc/F1Method.jl.svg?branch=master">
+  <a href="https://doi.org/10.5281/zenodo.2667835">
+    <img src="http://img.shields.io/badge/DOI-10.5281%20%2F%20zenodo.2667835-blue.svg?&style=flat-square">
   </a>
-  <a href='https://coveralls.io/github/briochemc/F1Method.jl'>
-    <img src='https://coveralls.io/repos/github/briochemc/F1Method.jl/badge.svg' alt='Coverage Status' />
+  <a href="https://github.com/briochemc/F1Method.jl/blob/master/LICENSE">
+    <img alt="License: MIT" src="https://img.shields.io/badge/License-MIT-blue.svg?&style=flat-square">
   </a>
 </p>
+
 <p>
-  <a href="https://ci.appveyor.com/project/briochemc/f1method-jl">
-    <img alt="Build Status" src="https://ci.appveyor.com/api/projects/status/prm2xfd6q5pba1om?svg=true">
+  <a href="https://travis-ci.com/briochemc/F1Method.jl">
+    <img alt="Build Status" src="https://img.shields.io/travis/com/briochemc/F1Method.jl/master?label=OSX/Linux/Windows&logo=travis&logocolor=white&style=flat-square">
   </a>
   <a href="https://codecov.io/gh/briochemc/F1Method.jl">
-    <img src="https://codecov.io/gh/briochemc/F1Method.jl/branch/master/graph/badge.svg" />
+    <img src="https://img.shields.io/codecov/c/github/briochemc/F1Method.jl/master?label=Codecov&logo=codecov&logoColor=white&style=flat-square">
   </a>
 </p>
 
@@ -62,8 +53,7 @@ The F-1 algorithm is **easy** to use, gives **accurate** results, and is computa
 
 - **Easy** — The F-1 algorithm basically just needs the user to provide a solver (for finding the steady-state), the mismatch function, `f`, the state function, `F`, and their derivatives, `∇ₓf` and `∇ₓF` w.r.t. the state `x`. 
     (Note these derivatives can be computed numerically, via the [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package for example.) 
-- **Accurate** — Thanks to dual and hyperdual numbers, the accuracy of the gradient and Hessian, as computed by the F-1 algorithm, are close to machine precision.
-    (The F-1 algorithm uses the [DualNumbers](https://github.com/JuliaDiff/DualNumbers.jl) and [HyperDualNumbers](https://github.com/JuliaDiff/HyperDualNumbers.jl) packages.)
+- **Accurate** — Thanks to ForwardDiff's nested dual numbers implementation, the accuracy of the gradient and Hessian, as computed by the F-1 algorithm, are close to machine precision.
 - **Fast** — The F-1 algorithm is as fast as if you derived analytical formulas for every first and second derivatives *and* used those in the most efficient way.
     This is because the bottleneck of such computations is the number of matrix factorizations, and the F-1 algorithm only requires a single one. In comparison, standard autodifferentiation methods that take the steady-state solver as a black box would require order `m` or `m^2` factorizations, where `m` is the number of parameters.
 

README.md

@@ -4,4 +4,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 
 [compat]
-Documenter = "~0.22"
+Documenter = "~0.24"

docs/Project.toml

@@ -7,7 +7,7 @@ refer to the Equation numbers in the above manuscript. A bibtex
 citation file is available in the GitHub repository.
 ======================================================================#
 
-using LinearAlgebra, DualNumbers, HyperDualNumbers, DiffEqBase
+using LinearAlgebra, ForwardDiff, DiffEqBase
 
 """
     Mem
@@ -21,21 +21,22 @@ Contains
 - `p`   the parameters 𝒑
 The `Mem`-type object should be initialized with `initialize_mem`.
 """
-mutable struct Mem 
+mutable struct Mem
     s     # 𝒔(𝒑)
     A     # factors of 𝐀 = ∇ₓ𝑭(𝒔,𝒑)
     ∇s    # ∇𝒔(𝒑)
     ∇ₓf   # ∇ₓ𝑓(𝒔,𝒑)
     p     # 𝒑
 end
 
+
 function update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     if p ≠ mem.p                      # only update mem if 𝒑 has changed
         update_solution!(F, ∇ₓF, mem, p, alg; options...)
         s, m = mem.s.u, length(p)
-        ∇ₚF = reduce(hcat, [𝔇(F(s, p + ε * e(j,m))) for j in 1:m]) # (2.7)
+        ∇ₚF = ForwardDiff.jacobian(p -> F(s,p), p)
         mem.A = factorize(∇ₓF(s,p))   # update factors of ∇ₓ𝑭(𝒔,𝒑)
-        mem.∇s .= mem.A \ -∇ₚF        # update ∇𝒔 (2.2)
+        mem.∇s .= mem.A \ -∇ₚF        # update ∇𝒔
         mem.∇ₓf .= ∇ₓf(s,p)           # update ∇ₓ𝑓(𝒔,𝒑)
         mem.p = p                     # update 𝒑
     end
@@ -73,8 +74,8 @@ Returns the gradient of the `objective` function using the F-1 method.
 function gradient(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     s, ∇s, m = mem.s, mem.∇s, length(p)
-    ∇ₚf = [𝔇(f(s,p + ε * e(j,m))) for j in 1:m]'    # (2.6)
-    return mem.∇ₓf * ∇s + ∇ₚf                       # (2.1)
+    ∇ₚf = ForwardDiff.jacobian(p -> [f(s,p)], p)
+    return mem.∇ₓf * ∇s + ∇ₚf
 end
 
 """
@@ -86,14 +87,8 @@ function hessian(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     s, A, ∇s, m = mem.s, mem.A, mem.∇s, length(p)
     A⁻ᵀ∇ₓfᵀ = vec(A' \ mem.∇ₓf') # independent of (𝑗,𝑘)
-    H, xⱼₖ = zeros(m,m), Vector{Hyper{Float64}}(undef, length(s))
-    for j in 1:m, k in j:m       # loop upper triangle (symmetry)
-        pⱼₖ = p + ε₁ * e(j,m) + ε₂ * e(k,m)              # hyperdual 𝒑
-        @views xⱼₖ .= s + ε₁ * ∇s[:,j] + ε₂ * ∇s[:,k]    # hyperdual 𝒙
-        H[j,k] = ℌ(f(xⱼₖ,pⱼₖ)) - ℌ(F(xⱼₖ,pⱼₖ))' * A⁻ᵀ∇ₓfᵀ    # (2.8)
-        j ≠ k ? H[k,j] = H[j,k] : nothing # Hessian symmetry
-    end
-    return H
+    H(λ) = f(s+∇s*λ, p+λ) - F(s+∇s*λ, p+λ)' * A⁻ᵀ∇ₓfᵀ
+    return ForwardDiff.hessian(H, zeros(m))
 end
 
 """
@@ -106,8 +101,4 @@ function initialize_mem(x, p)
     return Mem(copy(x), nothing, zeros(n,m), zeros(1,n), nothing)
 end
 
-e(j, m) = [i == j for i in 1:m]      # 𝑗th basis vector of ℝᵐ
-𝔇(x) = DualNumbers.dualpart.(x)      # dual part
-ℌ(x) = HyperDualNumbers.ε₁ε₂part.(x) # hyperdual part
-
 end