text:
using Test, Distributions, ForwardDiff

@testset "Normal" begin
    @test isa(convert(Normal{Float64}, Float16(0), Float16(1)),
              Normal{Float64})
    @test Inf === logpdf(Normal(0, 0), 0)
    @test -Inf === logpdf(Normal(), Inf)
    @test iszero(logcdf(Normal(0, 0), 0))
    @test iszero(logcdf(Normal(), Inf))
    @test -Inf === logccdf(Normal(0, 0), 0)
    @test iszero(logccdf(Normal(eps(), 0), 0))
    @test -Inf === quantile(Normal(), 0)
    @test iszero(quantile(Normal(), 0.5))
    @test Inf === quantile(Normal(), 1)
    @test -Inf === quantile(Normal(0, 0), 0)
    @test iszero(quantile(Normal(0, 0), 0.75))
    @test Inf === quantile(Normal(0, 0), 1)
    @test -Inf === quantile(Normal(0.25, 0), 0)
    @test 0.25 == quantile(Normal(0.25, 0), 0.95)
    @test Inf === quantile(Normal(0.25, 0), 1)
    @test Inf === cquantile(Normal(), 0)
    @test iszero(cquantile(Normal(), 0.5))
    @test -Inf === cquantile(Normal(), 1)
    @test Inf === cquantile(Normal(0, 0), 0)
    @test iszero(cquantile(Normal(0, 0), 0.75))
    @test -Inf === cquantile(Normal(0, 0), 1)
    @test Inf === cquantile(Normal(0.25, 0), 0)
    @test 0.25 == cquantile(Normal(0.25, 0), 0.95)
    @test -Inf === cquantile(Normal(0.25, 0), 1)
    @test -Inf === invlogcdf(Normal(), -Inf)
    @test isnan(invlogcdf(Normal(), NaN))
    @test Inf === invlogccdf(Normal(), -Inf)
    @test isnan(invlogccdf(Normal(), NaN))
    # test for #996 being fixed
    let d = Normal(0, 1), x = 1.0, ∂x = 2.0
        @inferred cdf(d, ForwardDiff.Dual(x, ∂x)) ≈ ForwardDiff.Dual(cdf(d, x), ∂x * pdf(d, x))
    end
end

@testset "Normal logpdf & pdf type inference" begin
    @test @inferred(pdf(Normal(0.0, 0.0), 0.0))           === Inf
    @test @inferred(pdf(Normal(0.0, 0.0), -1.0))          === 0.0
    @test @inferred(pdf(Normal(0.0, 0.0), 0.0f0))         === Inf
    @test @inferred(pdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(pdf(Normal(0.0f0, 0.0f0), 0.0))       === Inf
    @test @inferred(pdf(Normal(0.0f0, 0.0f0), 0.0f0))     === Inf32
    @test @inferred(pdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(pdf(Normal(NaN, 0.0), 0.0f0))         === NaN
    @test @inferred(pdf(Normal(NaN32, 0.0f0), 0.0f0))     === NaN32
    @test @inferred(pdf(Normal(0 // 1, 0 // 1), 0 // 1))  === Inf
    @test @inferred(pdf(Normal(0 // 1, 0 // 1), NaN))     === NaN
    @test @inferred(pdf(Normal(0.0, 0.0), BigInt(1)))     == big(0.0)
    @test @inferred(pdf(Normal(0.0, 0.0), BigFloat(1)))   == big(0.0)
    @test isequal(@inferred(pdf(Normal(0.0, 0.0), BigFloat(NaN))), big(NaN))

    @test @inferred(logpdf(Normal(0.0, 0.0), 0.0))           === Inf
    @test @inferred(logpdf(Normal(0.0, 0.0), -1.0))          === -Inf
    @test @inferred(logpdf(Normal(0.0, 0.0), 0.0f0))         === Inf
    @test @inferred(logpdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(logpdf(Normal(0.0f0, 0.0f0), 0.0))       === Inf
    @test @inferred(logpdf(Normal(0.0f0, 0.0f0), 0.0f0))     === Inf32
    @test @inferred(logpdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(logpdf(Normal(NaN, 0.0), 0.0f0))         === NaN
    @test @inferred(logpdf(Normal(NaN32, 0.0f0), 0.0f0))     === NaN32
    @test @inferred(logpdf(Normal(0 // 1, 0 // 1), 0 // 1))  === Inf
    @test @inferred(logpdf(Normal(0 // 1, 0 // 1), NaN))     === NaN
    @test @inferred(logpdf(Normal(0.0, 0.0), BigInt(1)))     == big(-Inf)
    @test @inferred(logpdf(Normal(0.0, 0.0), BigFloat(1)))   == big(-Inf)
    @test isequal(@inferred(logpdf(Normal(0.0, 0.0), BigFloat(NaN))), big(NaN))

    @test @inferred(cdf(Normal(0.0, 0.0), 0.0))           === 1.0
    @test @inferred(cdf(Normal(0.0, 0.0), -1.0))          === 0.0
    @test @inferred(cdf(Normal(0.0, 0.0), 0.0f0))         === 1.0
    @test @inferred(cdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(cdf(Normal(0.0f0, 0.0f0), 0.0))       === 1.0
    @test @inferred(cdf(Normal(0.0f0, 0.0f0), 0.0f0))     === 1.0f0
    @test @inferred(cdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(cdf(Normal(NaN, 0.0), 0.0f0))         === NaN
    @test @inferred(cdf(Normal(NaN32, 0.0f0), 0.0f0))     === NaN32
    @test @inferred(cdf(Normal(0 // 1, 0 // 1), 0 // 1))  === 1.0
    @test @inferred(cdf(Normal(0 // 1, 0 // 1), NaN))     === NaN
    @test @inferred(cdf(Normal(0.0, 0.0), BigInt(1)))     == big(1.0)
    @test @inferred(cdf(Normal(0.0, 0.0), BigFloat(1)))   == big(1.0)
    @test isequal(@inferred(cdf(Normal(0.0, 0.0), BigFloat(NaN))), big(NaN))

    @test @inferred(logcdf(Normal(0.0, 0.0), 0.0))           === -0.0
    @test @inferred(logcdf(Normal(0.0, 0.0), -1.0))          === -Inf
    @test @inferred(logcdf(Normal(0.0, 0.0), 0.0f0))         === -0.0
    @test @inferred(logcdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(logcdf(Normal(0.0f0, 0.0f0), 0.0))       === -0.0
    @test @inferred(logcdf(Normal(0.0f0, 0.0f0), 0.0f0))     === -0.0f0
    @test @inferred(logcdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(logcdf(Normal(NaN, 0.0), 0.0f0))         === NaN
    @test @inferred(logcdf(Normal(NaN32, 0.0f0), 0.0f0))     === NaN32
    @test @inferred(logcdf(Normal(0 // 1, 0 // 1), 0 // 1))  === -0.0
    @test @inferred(logcdf(Normal(0 // 1, 0 // 1), NaN))     === NaN
    @test @inferred(logcdf(Normal(0.0, 0.0), BigInt(1)))     == big(0.0)
    @test @inferred(logcdf(Normal(0.0, 0.0), BigFloat(1)))   == big(0.0)
    @test isequal(@inferred(logcdf(Normal(0.0, 0.0), BigFloat(NaN))), big(NaN))

    @test @inferred(ccdf(Normal(0.0, 0.0), 0.0))           === 0.0
    @test @inferred(ccdf(Normal(0.0, 0.0), -1.0))          === 1.0
    @test @inferred(ccdf(Normal(0.0, 0.0), 0.0f0))         === 0.0
    @test @inferred(ccdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(ccdf(Normal(0.0f0, 0.0f0), 0.0))       === 0.0
    @test @inferred(ccdf(Normal(0.0f0, 0.0f0), 0.0f0))     === 0.0f0
    @test @inferred(ccdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(ccdf(Normal(NaN, 0.0), 0.0f0))         === NaN
    @test @inferred(ccdf(Normal(NaN32, 0.0f0), 0.0f0))     === NaN32
    @test @inferred(ccdf(Normal(0 // 1, 0 // 1), 0 // 1))  === 0.0
    @test @inferred(ccdf(Normal(0 // 1, 0 // 1), NaN))     === NaN
    @test @inferred(ccdf(Normal(0.0, 0.0), BigInt(1)))     == big(0.0)
    @test @inferred(ccdf(Normal(0.0, 0.0), BigFloat(1)))   == big(0.0)
    @test isequal(@inferred(ccdf(Normal(0.0, 0.0), BigFloat(NaN))), big(NaN))

    @test @inferred(logccdf(Normal(0.0, 0.0), 0.0))           === -Inf
    @test @inferred(logccdf(Normal(0.0, 0.0), -1.0))          === -0.0
    @test @inferred(logccdf(Normal(0.0, 0.0), 0.0f0))         === -Inf
    @test @inferred(logccdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(logccdf(Normal(0.0f0, 0.0f0), 0.0))       === -Inf
    @test @inferred(logccdf(Normal(0.0f0, 0.0f0), 0.0f0))     === -Inf32
    @test @inferred(logccdf(Normal(0.0, 0.0), NaN))           === NaN
    @test @inferred(logccdf(Normal(NaN, 0.0), 0.0f0))         === NaN
    @test @inferred(logccdf(Normal(NaN32, 0.0f0), 0.0f0))     === NaN32
    @test @inferred(logccdf(Normal(0 // 1, 0 // 1), 0 // 1))  === -Inf
    @test @inferred(logccdf(Normal(0 // 1, 0 // 1), NaN))     === NaN
    @test @inferred(logccdf(Normal(0.0, 0.0), BigInt(1)))     == big(-Inf)
    @test @inferred(logccdf(Normal(0.0, 0.0), BigFloat(1)))   == big(-Inf)
    @test isequal(@inferred(logccdf(Normal(0.0, 0.0), BigFloat(NaN))), big(NaN))

    @test @inferred(quantile(Normal(1.0, 0.0), 0.0f0))     === -Inf
    @test @inferred(quantile(Normal(1.0, 0.0f0), 1.0))     ===  Inf
    @test @inferred(quantile(Normal(1.0f0, 0.0), 0.5))     ===  1.0
    @test @inferred(quantile(Normal(1.0f0, 0.0), NaN))     ===  NaN
    @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.0f0)) === -Inf32
    @test @inferred(quantile(Normal(1.0f0, 0.0f0), 1.0f0)) ===  Inf32
    @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.5f0)) ===  1.0f0
    @test @inferred(quantile(Normal(1.0f0, 0.0f0), NaN32)) ===  NaN32
    @test @inferred(quantile(Normal(1//1, 0//1), 1//2))    ===  1.0

    @test @inferred(cquantile(Normal(1.0, 0.0), 0.0f0))     ===  Inf
    @test @inferred(cquantile(Normal(1.0, 0.0f0), 1.0))     === -Inf
    @test @inferred(cquantile(Normal(1.0f0, 0.0), 0.5))     ===  1.0
    @test @inferred(cquantile(Normal(1.0f0, 0.0), NaN))     ===  NaN
    @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.0f0)) ===  Inf32
    @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 1.0f0)) === -Inf32
    @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.5f0)) ===  1.0f0
    @test @inferred(cquantile(Normal(1.0f0, 0.0f0), NaN32)) ===  NaN32
    @test @inferred(cquantile(Normal(1//1, 0//1), 1//2))    ===  1.0
end


text:
"""
The Multibreak package provides the `@multibreak` macro for breaking
out of, and optionally continuing, several nested loops at once.

The simplest use is of the form
```
using Multibreak

@multibreak for i = 1:5
    for j = 1:5
        break; break
    end
end
```
Read the tutorial at $(abspath(dirname(@__DIR__), "test",
"tutorial.jl")) or
https://github.com/GunnarFarneback/Multibreak.jl/blob/master/test/tutorial.jl
for the full documentation.
"""
module Multibreak
export @multibreak

function multibreak_transform_break_and_continue(args, break_labels,
                                                 continue_labels)
    out = Any[]
    nbreak = 0
    ncontinue = 0
    for arg in args
        if Meta.isexpr(arg, :break)
            if ncontinue > 0
                push!(out, :(error("multibreak: continue cannot precede a break")))
            end
            nbreak += 1
        elseif Meta.isexpr(arg, :continue)
            if ncontinue > 0
                push!(out, :(error("multibreak: multiple continue not allowed")))
            end
            ncontinue += 1
        elseif typeof(arg) == LineNumberNode
            push!(out, arg)
        else
            push!(out, arg)
        end
    end
    n = nbreak + ncontinue
    if  n > length(break_labels)
        push!(out, :(error("multibreak: not enough nested loops for requested multiple break/continue")))
    elseif n > 0
        push!(out, Expr(:symbolicgoto,
                        ifelse(ncontinue == 0,
                               break_labels[end - n + 1],
                               continue_labels[end - n + 1])))
    end
    return out
end

function multibreak_transform_ast(ast, loop_counter = [1],
                                  break_labels = Symbol[],
                                  continue_labels = Symbol[])
    if typeof(ast) != Expr
        return ast
    end

    if ast.head == :for || ast.head == :while
        n = loop_counter[1]
        break_labels = push!(copy(break_labels),
                             gensym(Symbol("loop", n, "break")))
        continue_labels = push!(copy(continue_labels),
                                gensym(Symbol("loop", n, "continue")))
        loop_counter[1] += 1
    end

    args = [multibreak_transform_ast(arg, loop_counter, break_labels,
                                     continue_labels) for arg in ast.args]

    if ast.head == :for || ast.head == :while
        arg2 = Expr(:block,
                    args[2],
                    Expr(:symboliclabel, last(continue_labels)))
        return Expr(:block,
                    Expr(ast.head, args[1], arg2),
                    Expr(:symboliclabel, last(break_labels)))
    elseif ast.head == :block
        return Expr(:block,
                    multibreak_transform_break_and_continue(args, break_labels,
                                                            continue_labels)...)
    end

    return Expr(ast.head, args...)
end

# TODO: This string repeats parts of the module documentation
# string. The common parts should be reused but first find out how to
# do that without running into scoping issues.
"""
`@multibreak` allows breaking out of, and optionally continuing,
several nested loops at once.

The simplest use is of the form
```
using Multibreak

@multibreak for i = 1:5
    for j = 1:5
        break; break
    end
end
```
Read the tutorial at $(abspath(dirname(@__DIR__), "test",
"tutorial.jl")) or
https://github.com/GunnarFarneback/Multibreak.jl/blob/master/test/tutorial.jl
for the full documentation.
"""
macro multibreak(blk)
    esc(multibreak_transform_ast(blk))
end

end

ext: jl


text:
"""
    ProjectedPointDistribution(m::Manifold, d, proj!, x)

Generates a random point in ambient space of `m` and projects it to `m`
using function `proj!`. Generated arrays are of type `TResult`, which can be
specified by providing the `x` argument.
"""
struct ProjectedPointDistribution{TResult,TM<:Manifold,TD<:Distribution,TProj} <:
       MPointDistribution{TM}
    manifold::TM
    d::TD
    proj!::TProj
end

function ProjectedPointDistribution(
    M::Manifold,
    d::Distribution,
    proj!,
    x::TResult,
) where {TResult}
    return ProjectedPointDistribution{TResult,typeof(M),typeof(d),typeof(proj!)}(
        M,
        d,
        proj!,
    )
end

function rand(rng::AbstractRNG, d::ProjectedPointDistribution{TResult}) where {TResult}
    x = convert(TResult, rand(rng, d.d))
    return d.proj!(d.manifold, x)
end

function _rand!(rng::AbstractRNG, d::ProjectedPointDistribution, x::AbstractArray{<:Number})
    _rand!(rng, d.d, x)
    return d.proj!(d.manifold, x)
end

support(d::ProjectedPointDistribution) = MPointSupport(d.manifold)

"""
    ProjectedFVectorDistribution(type::VectorBundleFibers, x, d, project_vector!)

Generates a random vector from ambient space of manifold `type.manifold`
at point `x` and projects it to vector space of type `type` using function
`project_vector!`, see [`project_vector`](@ref) for documentation.
Generated arrays are of type `TResult`.
"""
struct ProjectedFVectorDistribution{
    TResult,
    TSpace<:VectorBundleFibers,
    TX,
    TD<:Distribution,
    TProj,
} <: FVectorDistribution{TSpace,TX}
    type::TSpace
    x::TX
    d::TD
    project_vector!::TProj
end

function ProjectedFVectorDistribution(
    type::VectorBundleFibers,
    x,
    d::Distribution,
    project_vector!,
    xt::TResult,
) where {TResult}
    return ProjectedFVectorDistribution{
        TResult,
        typeof(type),
        typeof(x),
        typeof(d),
        typeof(project_vector!),
    }(
        type,
        x,
        d,
        project_vector!,
    )
end

function rand(rng::AbstractRNG, d::ProjectedFVectorDistribution{TResult}) where {TResult}
    v = convert(TResult, reshape(rand(rng, d.d), size(d.x)))
    return d.project_vector!(d.type, v, d.x, v)
end

function _rand!(
    rng::AbstractRNG,
    d::ProjectedFVectorDistribution,
    v::AbstractArray{<:Number},
)
    # calling _rand!(rng, d.d, v) doesn't work for all arrays types
    return copyto!(v, rand(rng, d))
end

support(tvd::ProjectedFVectorDistribution) = FVectorSupport(tvd.type, tvd.x)


text:
#special casing script tags so that they get removed and re-added on a change
#this ensures any js updates get run. We need this for Interact to work.
#see: https://github.com/shashi/Patchwork.jl/issues/14
println("Mathbox.jl needs to overwrite Patchwork's diff!(Elem{ns, tag}, Elem{ns, tag}, Any, Any) - you can ignore the warning below.\n")

import Patchwork: Patch, DictDiff, Overwrite, diff!

function diff!{ns, tag}(a::Elem{ns, tag}, b::Elem{ns, tag}, index, patches)
    if a === b return patches end

    patch = get(patches, index, Patch[])

    proppatch = Patchwork.diff(properties(a), properties(b))
    if !is(proppatch, nothing)
        patch = push!(patch, DictDiff(proppatch))
    end

    if tag == :script && ns == :xhtml
        scriptpatch = Dict()
        diff!(a.children, b.children, 0, scriptpatch)
        if !isempty(scriptpatch)
            return patches[index] = Overwrite(b)
        end
    else
        diff!(a.children, b.children, index, patches, parentpatch=patch)
    end

    if !isempty(patch)
        patches[index] = patch
    end
end


text:
push!(LOAD_PATH, "../src")

excludedfiles = [ "test_excluded.jl"];

function runtests()
    exename      = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir      = pwd()
    istest(f)    = endswith(f, ".jl") && startswith(basename(f), "test_")
    ispretest(f) = endswith(f, ".jl") && startswith(basename(f), "pretest_")    # NOTE: pretest_JustSayIt1.jl, pretest_JustSayIt2.jl test the Python related installations (restart needed after first installation!).
    testfiles    = sort(filter(istest,    vcat([joinpath.(root, files) for (root, dirs, files) in walkdir(testdir)]...)))
    pretestfiles = sort(filter(ispretest, vcat([joinpath.(root, files) for (root, dirs, files) in walkdir(testdir)]...)))

    nfail = 0
    printstyled("Testing package JustSayIt.jl\n"; bold=true, color=:white)
    for f in [pretestfiles; testfiles]
        println("")
        if f ∈ excludedfiles
            println("Test Skip:")
            println("$f")
            continue
        end
        try
            run(`$exename -O3 --startup-file=no $(joinpath(testdir, f))`)
        catch ex
            nfail += 1
        end
    end
    return nfail
end
exit(runtests())


text:
# precise values were computed with Maple 2018.2.1
@testset "Type stability" begin
  @testset for F = [Float32, Float64, BigFloat]
    x, y = SimpleBirthDeathProcess.gradient_hessian([F(1), F(1)], 1, 0, F(1))
    @test isa(x, Vector{F})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(1), F(1)], 1, 0, F(0))
    @test isa(x, Vector{F})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(0), F(1)], 1, 0, F(1))
    @test isa(x, Vector{F})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(1), F(0)], 1, 0, F(1))
    @test isa(x, Vector{F})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(0), F(0)], 1, 0, F(1))
    @test isa(x, Vector{F})
  end
end

@testset "Parameters are zero" begin
  x, y = SimpleBirthDeathProcess.gradient_hessian([1.0, 1.0], 1, 0, 0.0)
  @test x[1] === 0.0
  @test x[2] === 0.0

  x, y = SimpleBirthDeathProcess.gradient_hessian([0.0, 0.0], 1, 0, 1.0)
  @test x[1] === NaN
  @test x[2] === NaN

  x, y = SimpleBirthDeathProcess.gradient_hessian([0.0, 0.0], 2, 2, 1.0)
  @test x[1] === -2.0
  @test x[2] === -2.0
end

@testset "Death rate μ = 0" begin
  u = ((1, 0), (1, 1), (1, 100), (1000, 0), (1000, 500), (1000, 998), (1000, 999), (1000, 1000), (1000, 2000))
  v = ((0.001, 0.001), (0.001, 0.1), (0.001, 1.0), (0.001, 10.0), (0.001, 100.0),
       (  0.1, 0.001), (  0.1, 0.1), (  0.1, 1.0), (  0.1, 10.0), (  0.1, 100.0),
       (  1.0, 0.001), (  1.0, 0.1), (  1.0, 1.0), (  1.0, 10.0), (  1.0, 100.0),
       ( 10.0, 0.001), ( 10.0, 0.1), ( 10.0, 1.0), ( 10.0, 10.0),
       (100.0, 0.001), (100.0, 0.1), (100.0, 1.0))

  # case i = 1, j = 0
  λ1 = [-1000.0005000000833, -10.000500008333333, -1.0005000833333319, -0.10050083333194445, -0.10508331944775050e-1, -1000.0500008333333, -10.050083333194445, -1.0508331944775050, -0.15819767068693264, -0.10000454019910097, -1000.5000833333319, -10.508331944775050, -1.5819767068693264, -1.0000454019910097, -1.0000000000000000, -1005.0083333194445, -15.819767068693264, -10.000454019910097, -10.000000000000000, -1050.8331944775050, -100.00454019910097, -100.00000000000000]
  μ1 = fill(Float64(NaN), 22)

  # case i = 1, j = 1
  λ2 = [-0.10000000000000000e-2, -0.10000000000000000e-2, -0.10000000000000000e-2, -0.10000000000000000e-2, -0.10000000000000000e-2, -0.10000000000000000, -0.10000000000000000, -0.10000000000000000, -0.10000000000000000, -0.10000000000000000, -1.0, -1.0, -1.0, -1.0, -1.0, -10.0, -10.0, -10.0, -10.0, -100.0, -100.0, -100.0]
  μ2 = [-0.99999900000033333e-3, -0.99990000333325000e-3, -0.99900033325001666e-3, -0.99003325016638929e-3, -0.90325163928080854e-3, -0.99990000333325000e-1, -0.99003325016638929e-1, -0.90325163928080854e-1, -0.26424111765711536e-1, .80000907998595250e-1, -0.99900033325001666, -0.90325163928080854, -0.26424111765711536, .80000907998595250, .98000000000000000, -9.9003325016638929, -2.6424111765711536, 8.0000907998595250, 9.8000000000000000, -90.325163928080854, 80.000907998595250, 98.000000000000000]

  # case i = 1, j = 100
  λ3 = [98999.949500008250, 989.94950082500000, 98.949508249999863, 9.8495824998625003, 0.94032486253272991, 98994.950082500000, 984.95824998625003, 94.032486253272991, 5.6615693980063316, -0.99550520289004091e-1, 98949.508249999863, 940.32486253272991, 56.615693980063316, -0.99550520289004091, -1.0000000000000000, 98495.824998625003, 566.15693980063316, -9.9550520289004091, -10.000000000000000, 94032.486253272991, -99.550520289004091, -100.00000000000000]
  μ3 = [-0.50499957750016833e-1, -0.50495775168328988e-1, -0.50457766828988342e-1, -0.50079178995902986e-1, -0.46439070316410744e-1, -5.0495775168328988, -5.0079178995902986, -4.6439070316410744, -2.1459870421747642, 0.79596374218064201e-1, -50.457766828988342, -46.439070316410744, -21.459870421747642, 0.79596374218064201, 0.98000000000000000, -500.79178995902986, -214.59870421747642, 7.9596374218064201, 9.8000000000000000, -4643.9070316410744, 79.596374218064201, 98.000000000000000]

  # case i = 1000, j = 0
  λ4 = [-1000000.5000000833, -10000.500008333333, -1000.5000833333319, -100.50083333194445, -10.508331944775050, -1000050.0008333333, -10050.083333194445, -1050.8331944775050, -158.19767068693264, -100.00454019910097, -1000500.0833333319, -10508.331944775050, -1581.9767068693264, -1000.0454019910097, -1000.0000000000000, -1005008.3333194445, -15819.767068693264, -10000.454019910097, -10000.000000000000, -1050833.1944775050, -100004.54019910097, -100000.00000000000]
  μ4 = [-0.50000041666683333, -0.50004166833337639, -0.50041683337639722, -0.50418337647235800, -0.54337723598142662, -50.004166833337639, -50.418337647235800, -54.337723598142662, -113.63051215897188, -220164.65340786806, -500.41683337639722, -543.37723598142662, -1136.3051215897188, -2201646.5340786806, -0.26881171418161354e45, -5041.8337647235800, -11363.051215897188, -22016465.340786806, -0.26881171418161354e46, -54337.723598142662, -220164653.40786806, -0.26881171418161354e47]

  # case i = 1000, j = 500
  λ5 = [-500000.75000004167, -5000.7500041666667, -500.75004166666597, -50.750416665972224, -5.7541659723875248, -500075.00041666667, -5075.0416665972224, -575.41659723875248, -129.09883534346632, -100.00227009955048, -500750.04166666597, -5754.1659723875248, -1290.9883534346632, -1000.0227009955048, -1000.0000000000000, -507504.16665972224, -12909.883534346632, -10000.227009955048, -10000.000000000000, -575416.59723875248, -100002.27009955048, -100000.00000000000]
  μ5 = [-0.75000121033759135, -0.75012103623423997, -0.75121058744638998, -0.76212847856695186, -0.87364004716317359, -75.012103623423997, -76.212847856695186, -87.364004716317359, -233.16916999601915, -440850.68591396041, -751.21058744638998, -873.64004716317359, -2331.6916999601915, -4408506.8591396041, -0.53816212919325036e45, -7621.2847856695186, -23316.916999601915, -44085068.591396041, -0.53816212919325036e46, -87364.004716317359, -440850685.91396041, -0.53816212919325036e47]

  # case i = 1000, j = 998
  λ6 = [-2000.9990000001667, -20.999000016666667, -2.9990001666666639, -1.1990016666638889, -1.0190166638895501, -2099.9000016666667, -119.90016666638889, -101.90166638895501, -100.11639534137387, -100.00000908039820, -2999.0001666666639, -1019.0166638895501, -1001.1639534137387, -1000.0000908039820, -1000.0000000000000, -11990.016666638889, -10011.639534137387, -10000.000908039820, -10000.000000000000, -101901.66638895501, -100000.00908039820, -100000.00000000000]
  μ6 = [-0.99999800083366642, -1.0988000867465001, -1.9970012495001055, -10.979124800529834, -100.88560947870789, -109.88000867465001, -1097.9124800529834, -10088.560947870789, -108533.89212839828, -220024333.74274823, -1997.0012495001055, -100885.60947870789, -1085338.9212839828, -2200243337.4274823, -0.26854290246743193e48, -109791.24800529834, -10853389.212839828, -22002433374.274823, -0.26854290246743193e49, -10088560.947870789, -220024333742.74823, -0.26854290246743193e50]

  # case i = 1000, j = 999
  λ7 = [-1000.9995000000833, -10.999500008333333, -1.9995000833333319, -1.0995008333319444, -1.0095083319447750, -1099.9500008333333, -109.95008333319444, -100.95083319447750, -100.05819767068693, -100.00000454019910, -1999.5000833333319, -1009.5083319447750, -1000.5819767068693, -1000.0000454019910, -1000.0000000000000, -10995.008333319444, -10005.819767068693, -10000.000454019910, -10000.000000000000, -100950.83319447750, -100000.00454019910, -100000.00000000000]
  μ7 = fill(Float64(NaN), 22)

  # case i = 1000, j = 1000
  λ8 = [-1.0, -1.0, -1.0, -1.0, -1.0, -100.0, -100.0, -100.0, -100.0, -100.0, -1000.0, -1000.0, -1000.0, -1000.0, -1000.0, -10000.0, -10000.0, -10000.0, -10000.0, -100000.0, -100000.0, -100000.0]
  μ8 = [-0.99900000000033325, -0.90000000325000000, -0.25000001388888924e-6, 9.0000500001111112, 99.080026115675099, -90.000000325000000, 900.00500001111112, 9908.0026115675099, 108481.08672432000, 220024493.74457239, -0.25000001388888924e-3, 99080.026115675099, 1084810.8672432000, 2200244937.4457239, 0.26854290246743193e48, 90000.500001111112, 10848108.672432000, 22002449374.457239, 0.26854290246743193e49, 9908002.6115675099, 220024493744.57239, 0.26854290246743193e50]

  # case i = 1000, j = 2000
  λ9 = [999998.50000008333, 9998.5000083333333, 998.50008333333194, 98.500833331944448, 8.5083319447750496, 999850.00083333333, 9850.0833331944448, 850.83319447750496, -41.802329313067358, -99.995459800899031, 998500.08333333194, 8508.3319447750496, -418.02329313067358, -999.95459800899031, -1000.0000000000000, 985008.33331944448, -4180.2329313067358, -9999.5459800899031, -10000.000000000000, 850833.19447750496, -99995.459800899031, -100000.00000000000]
  μ9 = [-1.4999975853318353, -1.4997585381333236, -1.4975858311258006, -1.4759031040469013, -1.2633260093356671, -149.97585381333236, -147.59031040469013, -126.33260093356671, 60.565273388861400, 219884.60595632866, -1497.5858311258006, -1263.3260093356671, 605.65273388861400, 2198846.0595632866, 0.26827462783959234e45, -14759.031040469013, 6056.5273388861400, 21988460.595632866, 0.26827462783959234e46, -126332.60093356671, 219884605.95632866, 0.26827462783959234e47]

  λs = [λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9]
  μs = [μ1 μ2 μ3 μ4 μ5 μ6 μ7 μ8 μ9]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, λ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([λ, 0.0], i, j, t)
      @test x[1] ≈ λs[a, b]

      if (b == 1) || (b == 7)
        @test isnan(x[2])
      else
        @test x[2] ≈ μs[a, b]
      end
    end
  end
end

@testset "Birth rate λ = 0" begin
  u = ((1, 0), (1, 1), (1, 2), (1, 3), (1, 100),
       (1000, 0), (1000, 500), (1000, 1000), (1000, 1001), (1000, 1002), (1000, 2000))
  v = ((0.001, 0.001), (0.001, 0.1), (0.001, 1.0), (0.001, 10.0), (0.001, 100.0),
       (  0.1, 0.001), (  0.1, 0.1), (  0.1, 1.0), (  0.1, 10.0), (  0.1, 100.0),
       (  1.0, 0.001), (  1.0, 0.1), (  1.0, 1.0), (  1.0, 10.0), (  1.0, 100.0),
       ( 10.0, 0.001), ( 10.0, 0.1), ( 10.0, 1.0), ( 10.0, 10.0),
       (100.0, 0.001), (100.0, 0.1), (100.0, 1.0))

  # case i = 1, j = 0
  λ1 = [-.49999958333350000e-3, -.49995833499995695e-3, -.49958349995695278e-3, -.49584995702764239e-3, -.45995776441545389e-3, -.49995833499995695e-1, -.49584995702764239e-1, -.45995776441545389e-1, -.21409726569788410e-1, -.40861998033439283e-5, -.49958349995695278, -.45995776441545389, -.21409726569788410, -.40861998033439283e-4, -.36828752162606276e-43, -4.9584995702764239, -2.1409726569788410, -.40861998033439283e-3, -.36828752162606276e-42, -45.995776441545389, -.40861998033439283e-2, -.36828752162606276e-41]
  μ1 = [999.99950000008333, 9.9995000083333333, .99950008333333194, .99500833331944448e-1, .95083319447750496e-2, 999.95000083333333, 9.9500833331944448, .95083319447750496, .58197670686932642e-1, .45401991009687768e-5, 999.50008333333194, 9.5083319447750496, .58197670686932642, .45401991009687768e-4, .37200759760208360e-43, 995.00833331944448, 5.8197670686932642, .45401991009687768e-3, .37200759760208360e-42, 950.83319447750496, .45401991009687768e-2, .37200759760208360e-41]

  # case i = 1, j = 1
  λ2 = [-.99999900000033333e-3, -.99990000333325000e-3, -.99900033325001666e-3, -.99003325016638929e-3, -.90325163928080854e-3, -.99990000333325000e-1, -.99003325016638929e-1, -.90325163928080854e-1, -.26424111765711536e-1, .80000907998595250e-1, -.99900033325001666, -.90325163928080854, -.26424111765711536, .80000907998595250, .98000000000000000, -9.9003325016638929, -2.6424111765711536, 8.0000907998595250, 9.8000000000000000, -90.325163928080854, 80.000907998595250, 98.000000000000000]
  μ2 = [-.10000000000000000e-2, -.10000000000000000e-2, -.10000000000000000e-2, -.10000000000000000e-2, -.10000000000000000e-2, -.10000000000000000, -.10000000000000000, -.10000000000000000, -.10000000000000000, -.10000000000000000, -1., -1., -1., -1., -1., -10., -10., -10., -10., -100., -100., -100.]

  # case i = 1, j = 2
  λ3 = fill(Float64(NaN), 22)
  μ3 = [-1000.0015000000833, -10.001500008333333, -1.0015000833333319, -.10150083333194445, -.11508331944775050e-1, -1000.1500008333333, -10.150083333194445, -1.1508331944775050, -.25819767068693264, -.20000454019910097, -1001.5000833333319, -11.508331944775050, -2.5819767068693264, -2.0000454019910097, -2.0000000000000000, -1015.0083333194445, -25.819767068693264, -20.000454019910097, -20.000000000000000, -1150.8331944775050, -200.00454019910097, -200.00000000000000]

  # case i = 1, j = 3
  λ4 = [-.20000038333340000e-2, -.20003833400003361e-2, -.20038340003361445e-2, -.20384003364455497e-2, -.23903395557079497e-2, -.20003833400003361, -.20384003364455497, -.23903395557079497, -.68814964393585032, -1321.2279395160034, -2.0038340003361445, -2.3903395557079497, -6.8814964393585032, -13212.279395160034, -.16128702850896813e43, -20.384003364455497, -68.814964393585032, -132122.79395160034, -.16128702850896813e44, -239.03395557079497, -1321227.9395160034, -.16128702850896813e45]
  μ4 = [-2000.0020000001667, -20.002000016666667, -2.0020001666666639, -.20200166666388890, -.22016663889550099e-1, -2000.2000016666667, -20.200166666388890, -2.2016663889550099, -.41639534137386528, -.30000908039820194, -2002.0001666666639, -22.016663889550099, -4.1639534137386528, -3.0000908039820194, -3.0000000000000000, -2020.0166666388890, -41.639534137386528, -30.000908039820194, -30.000000000000000, -2201.6663889550099, -300.00908039820194, -300.00000000000000]

  # case i = 1, j = 100
  λ5 = [-.50500042270424997e-1, -.50504227209154007e-1, -.50542287245844974e-1, -.50924391770930393e-1, -.54899807241247986e-1, -5.0504227209154007, -5.0924391770930393, -5.4899807241247986, -11.495293725082332, -22241.204792222054, -50.542287245844974, -54.899807241247986, -114.95293725082332, -222412.04792222054, -.27155469085693613e44, -509.24391770930393, -1149.5293725082332, -2224120.4792222054, -.27155469085693613e45, -5489.9807241247986, -22241204.792222054, -.27155469085693613e46]
  μ5 = [-99000.050500008250, -990.05050082500000, -99.050508249999863, -9.9505824998625003, -1.0413248625327299, -99005.050082500000, -995.05824998625003, -104.13248625327299, -15.761569398006332, -10.000449479710996, -99050.508249999863, -1041.3248625327299, -157.61569398006332, -100.00449479710996, -100.00000000000000, -99505.824998625003, -1576.1569398006332, -1000.0449479710996, -1000.0000000000000, -104132.48625327299, -10000.449479710996, -10000.000000000000]

  # case i = 1000, j = 0
  λ6 = [-.49999958333350000, -.49995833499995695, -.49958349995695278, -.49584995702764239, -.45995776441545389, -49.995833499995695, -49.584995702764239, -45.995776441545389, -21.409726569788410, -.40861998033439283e-2, -499.58349995695278, -459.95776441545389, -214.09726569788410, -.40861998033439283e-1, -.36828752162606276e-40, -4958.4995702764239, -2140.9726569788410, -.40861998033439283, -.36828752162606276e-39, -45995.776441545389, -4.0861998033439283, -.36828752162606276e-38]
  μ6 = [999999.50000008333, 9999.5000083333333, 999.50008333333194, 99.500833331944448, 9.5083319447750496, 999950.00083333333, 9950.0833331944448, 950.83319447750496, 58.197670686932642, .45401991009687768e-2, 999500.08333333194, 9508.3319447750496, 581.97670686932642, .45401991009687768e-1, .37200759760208360e-40, 995008.33331944448, 5819.7670686932642, .45401991009687768, .37200759760208360e-39, 950833.19447750496, 4.5401991009687768, .37200759760208360e-38]

  # case i = 1000, j = 500
  λ7 = [-.74999879366292465, -.74987936876736038, -.74879391256996842, -.73796152217638516, -.63176278847795166, -74.987936876736038, -73.796152217638516, -63.176278847795166, 30.174345657161533, 109722.71751532172, -748.79391256996842, -631.76278847795166, 301.74345657161533, 1097227.1751532172, .13386930676309896e45, -7379.6152217638516, 3017.4345657161533, 10972271.751532172, .13386930676309896e46, -63176.278847795166, 109722717.51532172, .13386930676309896e47]
  μ7 = [499999.25000004167, 4999.2500041666667, 499.25004166666597, 49.250416665972224, 4.2541659723875248, 499925.00041666667, 4925.0416665972224, 425.41659723875248, -20.901164656533679, -49.997729900449516, 499250.04166666597, 4254.1659723875248, -209.01164656533679, -499.97729900449516, -500.00000000000000, 492504.16665972224, -2090.1164656533679, -4999.7729900449516, -5000.0000000000000, 425416.59723875248, -49997.729900449516, -50000.000000000000]

  # case i = 1000, j = 1000
  λ8 = [-.99900000000033325, -.90000000325000000, -.25000001388888924e-6, 9.0000500001111112, 99.080026115675099, -90.000000325000000, 900.00500001111112, 9908.0026115675099, 108481.08672432000, 220024493.74457239, -.25000001388888924e-3, 99080.026115675099, 1084810.8672432000, 2200244937.4457239, .26854290246743193e48, 90000.500001111112, 10848108.672432000, 22002449374.457239, .26854290246743193e49, 9908002.6115675099, 220024493744.57239, .26854290246743193e50]
  μ8 = [-1., -1., -1., -1., -1., -100., -100., -100., -100., -100., -1000., -1000., -1000., -1000., -1000., -10000., -10000., -10000., -10000., -100000., -100000., -100000.]

  # case i = 1000, j = 1001
  λ9 = fill(Float64(NaN), 22)
  μ9 = [-1001.0005000000833, -11.000500008333333, -2.0005000833333319, -1.1005008333319444, -1.0105083319447750, -1100.0500008333333, -110.05008333319444, -101.05083319447750, -100.15819767068693, -100.10000454019910, -2000.5000833333319, -1010.5083319447750, -1001.5819767068693, -1001.0000454019910, -1001.0000000000000, -11005.008333319444, -10015.819767068693, -10010.000454019910, -10010.000000000000, -101050.83319447750, -100100.00454019910, -100100.00000000000]

  # case i = 1000, j = 1002
  λ10 = [-1.0020020008336671, -1.1012000867535001, -2.0030012505001056, -11.021125200531279, -101.28794959317297, -110.12000867535001, -1102.1125200531279, -10128.794959317297, -108968.62671672793, -220905752.70567148, -2003.0012505001056, -101287.94959317297, -1089686.2671672793, -2209057527.0567148, -.26961868694758675e48, -110211.25200531279, -10896862.671672793, -22090575270.567148, -.26961868694758675e49, -10128794.959317297, -220905752705.67148, -.26961868694758675e50]
  μ10 = [-2001.0010000001667, -21.001000016666667, -3.0010001666666639, -1.2010016666638889, -1.0210166638895501, -2100.1000016666667, -120.10016666638889, -102.10166638895501, -100.31639534137387, -100.20000908039820, -3001.0001666666639, -1021.0166638895501, -1003.1639534137387, -1002.0000908039820, -1002.0000000000000, -12010.016666638889, -10031.639534137387, -10020.000908039820, -10020.000000000000, -102101.66638895501, -100200.00908039820, -100200.00000000000]

  # case i = 1000, j = 2000
  λ11 = [-1.5000024186691687, -1.5002418718670768, -1.5024191688785827, -1.5242368968264347, -1.7470793256997726, -150.02418718670768, -152.42368968264347, -174.70793256997726, -466.12045451745941, -881259.55795312004, -1502.4191688785827, -1747.0793256997726, -4661.2045451745941, -8812595.5795312004, -.10757850183164074e46, -15242.368968264347, -46612.045451745941, -88125955.795312004, -.10757850183164074e47, -174707.93256997726, -881259557.95312004, -.10757850183164074e48]
  μ11 = [-1000001.5000000833, -10001.500008333333, -1001.5000833333319, -101.50083333194445, -11.508331944775050, -1000150.0008333333, -10150.083333194445, -1150.8331944775050, -258.19767068693264, -200.00454019910097, -1001500.0833333319, -11508.331944775050, -2581.9767068693264, -2000.0454019910097, -2000.0000000000000, -1015008.3333194445, -25819.767068693264, -20000.454019910097, -20000.000000000000, -1150833.1944775050, -200004.54019910097, -200000.00000000000]

  λs = [λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10 λ11]
  μs = [μ1 μ2 μ3 μ4 μ5 μ6 μ7 μ8 μ9 μ10 μ11]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, μ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([0.0, μ], i, j, t)
      @test x[2] ≈ μs[a, b]

      if (b == 3) || (b == 9)
        @test isnan(x[1])
      else
        @test x[1] ≈ λs[a, b]
      end
    end
  end
end

@testset "Equal rates" begin
  u = ((1, 0), (1, 1), (100, 0), (100, 1), (100, 50), (100, 100), (100, 200))
  v = ((0.001, 0.001), (0.001, 0.1), (0.001, 1.0), (0.001, 10.0), (0.001, 100.0),
       (  0.1, 0.001), (  0.1, 0.1), (  0.1, 1.0), (  0.1, 10.0), (  0.1, 100.0),
       (  1.0, 0.001), (  1.0, 0.1), (  1.0, 1.0), (  1.0, 10.0), (  1.0, 100.0),
       ( 10.0, 0.001), ( 10.0, 0.1), ( 10.0, 1.0), ( 10.0, 10.0),
       (100.0, 0.001), (100.0, 0.1), (100.0, 1.0))

  # case i = 1, j = 0
  λ1 = [-.49999950000050000e-3, -.49995000499950005e-3, -.49950049950049950e-3, -.49504950495049505e-3, -.45454545454545455e-3, -.49995000499950005e-1, -.49504950495049505e-1, -.45454545454545455e-1, -.25000000000000000e-1, -.45454545454545455e-2, -.49950049950049950, -.45454545454545455, -.25000000000000000, -.45454545454545455e-1, -.49504950495049505e-2, -4.9504950495049505, -2.5000000000000000, -.45454545454545455, -.49504950495049505e-1, -45.454545454545455, -4.5454545454545455, -.49504950495049505]
  μ1 = [999.99950000050000, 9.9995000499950005, .99950049950049950, .99504950495049505e-1, .95454545454545455e-2, 999.95000499950005, 9.9504950495049505, .95454545454545455, .75000000000000000e-1, .54545454545454545e-2, 999.50049950049950, 9.5454545454545455, .75000000000000000, .54545454545454545e-1, .50495049504950495e-2, 995.04950495049505, 7.5000000000000000, .54545454545454545, .50495049504950495e-1, 954.54545454545455, 5.4545454545454545, .50495049504950495]

  # case i = 1, j = 1
  λ2 = [-.99999900000100000e-3, -.99990000999900010e-3, -.99900099900099900e-3, -.99009900990099010e-3, -.90909090909090909e-3, -.99990000999900010e-1, -.99009900990099010e-1, -.90909090909090909e-1, -.50000000000000000e-1, -.90909090909090909e-2, -.99900099900099900, -.90909090909090909, -.50000000000000000, -.90909090909090909e-1, -.99009900990099010e-2, -9.9009900990099010, -5., -.90909090909090909, -.99009900990099010e-1, -90.909090909090909, -9.0909090909090909, -.99009900990099010]
  μ2 = [-.99999900000100000e-3, -.99990000999900010e-3, -.99900099900099900e-3, -.99009900990099010e-3, -.90909090909090909e-3, -.99990000999900010e-1, -.99009900990099010e-1, -.90909090909090909e-1, -.50000000000000000e-1, -.90909090909090909e-2, -.99900099900099900, -.90909090909090909, -.50000000000000000, -.90909090909090909e-1, -.99009900990099010e-2, -9.9009900990099010, -5., -.90909090909090909, -.99009900990099010e-1, -90.909090909090909, -9.0909090909090909, -.99009900990099010]

  # case i = 100, j = 0
  λ3 = [-.49999950000050000e-1, -.49995000499950005e-1, -.49950049950049950e-1, -.49504950495049505e-1, -.45454545454545455e-1, -4.9995000499950005, -4.9504950495049505, -4.5454545454545455, -2.5000000000000000, -.45454545454545455, -49.950049950049950, -45.454545454545455, -25., -4.5454545454545455, -.49504950495049505, -495.04950495049505, -250., -45.454545454545455, -4.9504950495049505, -4545.4545454545455, -454.54545454545455, -49.504950495049505]
  μ3 = [99999.950000050000, 999.95000499950005, 99.950049950049950, 9.9504950495049505, .95454545454545455, 99995.000499950005, 995.04950495049505, 95.454545454545455, 7.5000000000000000, .54545454545454545, 99950.049950049950, 954.54545454545455, 75., 5.4545454545454545, .50495049504950495, 99504.950495049505, 750., 54.545454545454545, 5.0495049504950495, 95454.545454545455, 545.45454545454545, 50.495049504950495]

  # case i = 100, j = 1
  λ4 = [-.50499949500050500e-1, -.50494950504949505e-1, -.50449550449550450e-1, -.50000000000000000e-1, -.45909090909090909e-1, -5.0494950504949505, -5., -4.5909090909090909, -2.5250000000000000, -.45909090909090909, -50.449550449550450, -45.909090909090909, -25.250000000000000, -4.5909090909090909, -.50000000000000000, -500., -252.50000000000000, -45.909090909090909, -5., -4590.9090909090909, -459.09090909090909, -50.]
  μ4 = [98999.949500050500, 989.94950504949505, 98.949550449550450, 9.8500000000000000, .94409090909090909, 98994.950504949505, 985., 94.409090909090909, 7.3750000000000000, .53090909090909091, 98949.550449550450, 944.09090909090909, 73.750000000000000, 5.3090909090909091, .49000000000000000, 98500., 737.50000000000000, 53.090909090909091, 4.9000000000000000, 94409.090909090909, 530.90909090909091, 49.]

  # case i = 100, j = 50
  λ5 = [-.74999876960859314e-1, -.74987696828491296e-1, -.74877035844303681e-1, -.73777168447797450e-1, -.63507697938872812e-1, -7.4987696828491296, -7.3777168447797450, -6.3507697938872812, -2.1057046979865772, -.24835227282408946, -74.877035844303681, -63.507697938872812, -21.057046979865772, -2.4835227282408946, -.25481471624824203, -737.77168447797450, -210.57046979865772, -24.835227282408946, -2.5481471624824203, -6350.7697938872812, -248.35227282408946, -25.481471624824203]
  μ5 = [49999.925000123039, 499.92501230317151, 49.925122964155696, 4.9262228315522025, .43649230206112719, 49992.501230317151, 492.62228315522025, 43.649230206112719, 2.8942953020134228, .25164772717591054, 49925.122964155696, 436.49230206112719, 28.942953020134228, 2.5164772717591054, .24518528375175797, 49262.228315522025, 289.42953020134228, 25.164772717591054, 2.4518528375175797, 43649.230206112719, 251.64772717591054, 24.518528375175797]

  # case i = 100, j = 100
  λ6 = [-.99990000000149985e-1, -.99000050981668072e-1, -.90049754058324478e-1, -.30223650198890085e-1, -.25338869475809811e-2, -9.9000050981668072, -3.0223650198890085, -.25338869475809811, -.25125628140703518e-1, -.25990862764163180e-2, -90.049754058324478, -2.5338869475809811, -.25125628140703518, -.25990862764163180e-1, -.43313158689553725e-2, -302.23650198890085, -2.5125628140703518, -.25990862764163180, -.43313158689553725e-1, -253.38869475809811, -2.5990862764163180, -.43313158689553725]
  μ6 = [-.99990000000149985e-1, -.99000050981668072e-1, -.90049754058324478e-1, -.30223650198890085e-1, -.25338869475809811e-2, -9.9000050981668072, -3.0223650198890085, -.25338869475809811, -.25125628140703518e-1, -.25990862764163180e-2, -90.049754058324478, -2.5338869475809811, -.25125628140703518, -.25990862764163180e-1, -.43313158689553725e-2, -302.23650198890085, -2.5125628140703518, -.25990862764163180, -.43313158689553725e-1, -253.38869475809811, -2.5990862764163180, -.43313158689553725]

  # case i = 100, j = 200
  λ7 = [99999.850000248020, 999.85002480048006, 99.850247869667363, 9.8524650621838242, .87316511644498909, 99985.002480048006, 985.24650621838242, 87.316511644498909, 5.8110367892976589, .50601018116488633, 99850.247869667363, 873.16511644498909, 58.110367892976589, 5.0601018116488633, .49724365359025618, 98524.650621838242, 581.10367892976589, 50.601018116488633, 4.9724365359025618, 87316.511644498909, 506.01018116488633, 49.724365359025618]
  μ7 = [-.14999975198034802, -.14997519951993646, -.14975213033263712, -.14753493781617583, -.12683488355501091, -14.997519951993646, -14.753493781617583, -12.683488355501091, -4.1889632107023411, -.49398981883511367, -149.75213033263712, -126.83488355501091, -41.889632107023411, -4.9398981883511367, -.50275634640974382, -1475.3493781617583, -418.89632107023411, -49.398981883511367, -5.0275634640974382, -12683.488355501091, -493.98981883511367, -50.275634640974382]

  λs = [λ1 λ2 λ3 λ4 λ5 λ6 λ7]
  μs = [μ1 μ2 μ3 μ4 μ5 μ6 μ7]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, λ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([λ, λ], i, j, t)
      @test x[1] ≈ λs[a, b]
      @test x[2] ≈ μs[a, b]
    end
  end
end

@testset "Unequal rates" begin
  u = ((1, 0), (1, 1), (100, 0), (100, 1), (100, 2), (100, 50), (100, 100), (100, 200))

  # θ = (λ - μ) * t
  # ω = log(λ / μ)
  # compute the gradient at (approximate) points (θ, ω) at different 't':
  # (-20, -10), (-10, -7), (-5, -2), (-1, -5), (-0.5, -0.1), (-0.0001, -0.5)
  # (0.0001, 0.5), (0.5, 0.1), (1, 5), (5, 2), (10, 7), (20, 10)
  v = ((0.001,      0.5, 20_001.0), (0.001,    10.0, 10_010.0), (0.001,   800.0, 5_800.0), (0.001,     7.0, 1_007.0), (0.001,  4_700.0, 5_200.0), (0.001,     0.15,     0.25),
       (0.001,     0.25,     0.15), (0.001, 5_200.0,  4_700.0), (0.001, 1_007.0,     7.0), (0.001, 5_800.0,   800.0), (0.001, 10_010.0,    10.0), (0.001, 20_001.0,      0.5),
       (  1.0,   0.0005,   20.005), (  1.0,   0.005,   10.005), (  1.0,     0.8,     5.8), (  1.0,   0.007,   1.007), (  1.0,      4.7,     5.2), (  1.0, 0.000154, 0.000254),
       (  1.0, 0.000254, 0.000154), (  1.0,     5.2,      4.7), (  1.0,   1.007,   0.007), (  1.0,     5.8,     0.8), (  1.0,   10.005,   0.005), (  1.0,   20.005,   0.0005),
       (100.0,   9.0e-6, 0.200009), (100.0,  9.0e-5,  0.10009), (100.0,  0.0078,  0.0578), (100.0,  7.0e-5, 0.01007), (100.0,   0.0475,  0.0525), (100.0, 1.542e-6, 2.542e-6),
       (100.0, 2.542e-6, 1.542e-6), (100.0,  0.0525,   0.0475), (100.0, 0.01007,  0.7e-4), (100.0,  0.0578,  0.0078), (100.0,  0.10009,  0.9e-4), (100.0, 0.200009,   9.0e-6))

  # case i = 1, j = 0
  λ1 = [-.19570707918544786e-11, -.40821178708009919e-7, -.46906227010616650e-5, -.21315409211113093e-3, -.69899957089946077e-4, -.49988335624552032e-3, -.49991668291348118e-3, -.97583272912637388e-4, -.41618175383103213e-3, -.16672076034240454e-3, -.99854747794594254e-4, -.49997498064924518e-4, -.19492786674769594e-8, -.40841578171455945e-4, -.46906227010616650e-2, -.21315409211113093, -.69899957089946077e-1, -.49988135711991863, -.49991468365460000, -.97583272912637388e-1, -.41618175383103213, -.16672076034240454, -.99904647952845889e-1, -.49987501072373520e-1, -.19580078352125812e-6, -.40825256970283447e-2, -.47067585408845584, -21.315409211113093, -6.9319464697466986, -49.988125716405827, -49.991458369207570, -9.6772881176829805, -41.618175383103213, -16.729428323168471, -9.9864723838002413, -4.9997748040189202]
  μ1 = [.20600692320685258e-11, .45352105505645797e-7, .56930327610437369e-5, .57686690555327770e-3, .94724419395054920e-4, 3.9995000833170865, 6.6661667833104211, .14286600035686243e-3, .14264398876503173, .12453093772989383e-2, .99999959178821292e-1, 1.9999999999980429, .20518454250070254e-8, .45377034660357587e-4, .56930327610437369e-2, .57686690555327770, .94724419395054920e-1, 3936.5079593320934, 6493.0066121493736, .14286600035686243, 142.64398876503173, 1.2453093772989383, 199.99995915842183, 1999.9999999980507, .20610562417692556e-6, .45357089163137749e-2, .57160973911526626, 57.686690555327770, 9.3703309299360672, 393341.03922613119, 648458.44248388152, 14.120685109200670, 14264.398876503173, 127.73445235103975, 11111.107028585414, 111111.11111091531]

  # case i = 1, j = 1
  λ2 = [.90000250014360659e-3, .80000916162460336e-3, .60418606441048513e-3, -.26235477048360901e-3, -.56811749106295586e-4, -.99975004582437675e-3, -.99985003249362625e-3, -.27815471089887134e-3, -.99631692140071710e-3, -.94700883049741754e-3, -.99980029957120788e-3, -.99999750018759721e-3, .90002249514388002, .80000912082567646, .60418606441048513, -.26235477048360901, -.56811749106295586e-1, -.99974604757317337, -.99984603397586390, -.27815471089887134, -.99631692140071710, -.94700883049741754, -.99990009988771115, -.99999750118718439, 90.000000020629158, 80.000915346807881, 60.415379273083935, -26.235477048360901, -5.5650764321337403, -99.974584766145266, -99.984583405081529, -27.653392742725618, -99.631692140071710, -94.815587627597788, -99.982025165802420, -99.999550020268566]
  μ2 = [-.99999750018759721e-3, -.99980029957120788e-3, -.94700883049741754e-3, -.99631692140071710e-3, -.27815471089887134e-3, -.99985003249362625e-3, -.99975004582437675e-3, -.56811749106295586e-4, -.26235477048360901e-3, .60418606441048513e-3, .80000916162460336e-3, .90000250014360659e-3, -.99999750118718439, -.99990009988771115, -.94700883049741754, -.99631692140071710, -.27815471089887134, -.99984603397586390, -.99974604757317337, -.56811749106295586e-1, -.26235477048360901, .60418606441048513, .80000912082567646, .90002249514388002, -99.999550020268566, -99.982025165802420, -94.815587627597788, -99.631692140071710, -27.653392742725618, -99.984583405081529, -99.974584766145266, -5.5650764321337403, -26.235477048360901, 60.415379273083935, 80.000915346807881, 90.000000020629158]

  # case i = 100, j = 0
  λ3 = [-.19570707918544786e-9, -.40821178708009919e-5, -.46906227010616650e-3, -.21315409211113093e-1, -.69899957089946077e-2, -.49988335624552032e-1, -.49991668291348118e-1, -.97583272912637388e-2, -.41618175383103213e-1, -.16672076034240454e-1, -.99854747794594254e-2, -.49997498064924518e-2, -.19492786674769594e-6, -.40841578171455945e-2, -.46906227010616650, -21.315409211113093, -6.9899957089946077, -49.988135711991863, -49.991468365460000, -9.7583272912637388, -41.618175383103213, -16.672076034240454, -9.9904647952845889, -4.9987501072373520, -.19580078352125812e-4, -.40825256970283447, -47.067585408845584, -2131.5409211113093, -693.19464697466986, -4998.8125716405827, -4999.1458369207570, -967.72881176829805, -4161.8175383103213, -1672.9428323168471, -998.64723838002413, -499.97748040189202]
  μ3 = [.20600692320685258e-9, .45352105505645797e-5, .56930327610437369e-3, .57686690555327770e-1, .94724419395054920e-2, 399.95000833170865, 666.61667833104211, .14286600035686243e-1, 14.264398876503173, .12453093772989383, 9.9999959178821292, 199.99999999980429, .20518454250070254e-6, .45377034660357587e-2, .56930327610437369, 57.686690555327770, 9.4724419395054920, 393650.79593320934, 649300.66121493736, 14.286600035686243, 14264.398876503173, 124.53093772989383, 19999.995915842183, 199999.99999980507, .20610562417692556e-4, .45357089163137749, 57.160973911526626, 5768.6690555327770, 937.03309299360672, 39334103.922613119, 64845844.248388152, 1412.0685109200670, 1426439.8876503173, 12773.445235103975, 1111110.7028585414, 11111111.111091531]

  # case i = 100, j = 1
  λ4 = [.90000230639359820e-3, .79596786493251037e-3, .13981441700538030e-3, -.21364609889485571e-1, -.69769075010109572e-2, -.50488202314130889e-1, -.50491601640928263e-1, -.99388987292499727e-2, -.42198310550672898e-1, -.17452364104395467e-1, -.10885420331236039e-1, -.59497498086151245e-2, .90002230216529194, .79596580458670233, .13981441700538030, -21.364609889485571, -6.9769075010109572, -50.488000402445118, -50.491399715781264, -9.9388987292499727, -42.198310550672898, -17.452364104395467, -10.890460247219454, -5.9487601073521629, 89.999980636351590, 79.596745302802075, 13.818469718326807, -2136.4609889485571, -691.82777693705690, -5048.7990306903222, -5049.1389619566310, -985.70491639334068, -4219.8310550672898, -1751.0289916212764, -1088.6427911620263, -594.97725561814167]
  μ4 = [-.99999729624074324e-3, -.99531044112614895e-3, -.38339858715408759e-3, .56113506728373776e-1, .90995628092115657e-2, 395.94950839835907, 659.94951179768587, .14086922286223085e-1, 14.121492532967657, .12388981441700538, 9.9007959678649325, 198.00090000230639, -.99999729805448732, -.99540777345633575, -.38339858715408759, 56.113506728373776, 9.0995628092115657, 389713.28812784327, 642806.65485674041, 14.086922286223085, 14121.492532967657, 123.88981441700538, 19800.795965804587, 198000.90002230217, -99.999529615811772, -99.532989983087356, -38.226223455186428, 5611.3506728373776, 900.00936932094503, 38940662.898803582, 64197285.831319504, 1392.3827493787326, 1412149.2532967657, 12706.126162026019, 1100079.5967453028, 11000089.999980636]

  # case i = 100, j = 2
  λ5 = [1.9927787122357217, .32163692475737808e-1, .97425715959198324e-3, -.21391743089932451e-1, -.69203236347768731e-2, -.50988063953204707e-1, -.50991531960205365e-1, -.10027203421813526e-1, -.42778266788065129e-1, -.17954082799035675e-1, -.11508040738856832e-1, -.59536091560667918e-2, 1992.8094101534447, 37.861954491218859, .97425715959198324, -21.391743089932451, -6.9203236347768731, -50.987859961585254, -50.991327954991404, -10.027203421813526, -42.778266788065129, -17.954082799035675, -11.626596058432330, -5.9526049167402596, 110949.97705808668, 3315.8745224852283, 97.958885888721212, -2139.1743089932451, -686.12894646219536, -5098.7849762047094, -5099.1317754773521, -994.39253155393727, -4277.8266788065129, -1801.7526974316092, -1152.8870791430617, -595.19214856651272]
  μ5 = [-.10538541438174042e-2, -.17178309486470416e-2, -.10575310748977438e-2, .54540501831597234e-1, .88189504243403200e-2, 391.94900846803979, 653.28234526938013, .13930740195010361e-1, 13.978608256910068, .12347425715959198, 9.8321636924757378, 197.99277871223572, -1.0538296105668030, -1.8314936096567176, -1.0575310748977438, 54.540501831597234, 8.8189504243403200, 385775.78032558832, 636312.64850367478, 13.930740195010361, 13978.608256910068, 123.47425715959198, 19637.861954491219, 197992.80941015344, -105.21419757430737, -173.76828605683934, -106.25096732780302, 5454.0501831597234, 872.27413511272940, 38547221.875305561, 63548727.414764392, 1377.0289482746467, 1397860.8256910068, 12662.061449991285, 1092204.7634113741, 10999838.865946976]

  # case i = 100, j = 50
  λ6 = [97.626946975214449, 3.7588588469969708, .54473827961662456e-1, .27612412855101704e-1, -.13920642282632803e-2, -.74969660300858135e-1, -.74981129888725301e-1, -.83002754798428564e-2, -.70208264569878180e-1, -.25430656334196357e-1, -.21243401325429413e-1, -.61274573649451976e-2, 97628.277448243412, 6747.6436818595847, 54.473827961662456, 27.612412855101704, -1.3920642282632803, -74.969168275231885, -74.980637843082738, -8.3002754798428564, -70.208264569878180, -25.430656334196357, -24.683222535856405, -6.1258372797256635, 5432659.2333512069, 411665.24643513353, 5561.5356395185165, 2761.2412855101704, -137.80601377839591, -7496.9143674013607, -7498.0613240863636, -822.22312733527382, -7020.8264569878180, -2559.9435412482824, -2172.3816223275053, -605.13875780174958]
  μ6 = [-.36275823586955101e-2, -.16248396330424418e-1, -.16809966679023943e-1, -.20555831600662688e-1, .13151091355417589e-2, 199.92501887011127, 333.25836367303248, .92462336440771452e-2, 7.1704695557122446, .11697382796166246, 8.7588588469969708, 197.62694697521445, -3.6264621235147163, -19.685721286481093, -16.809966679023943, -20.555831600662688, 1.3151091355417589, 196775.41306294432, 324600.35550704944, 9.2462336440771452, 7170.4695557122446, 116.97382796166246, 16747.643681859585, 197628.27744824341, -355.15000729552236, -1672.8312176916776, -1694.8916381340956, -2055.5831600662688, 130.15782504567856, 19662053.472900933, 32417924.616112495, 914.82556516897251, 717046.95557122446, 11971.792049774927, 967220.80199068908, 10988214.788906762]

  # case i = 100, j = 100
  λ7 = [197.28754575851234, 8.2329736217213995, .11365599575369508, 1.0856508893970181, .50988841320189165e-2, -.97500473079758485e-1, -.98500284337856984e-1, -.49478791660518938e-2, -.90616161433438351e-1, -.28953817235122780e-1, -.26012207465436708e-1, -.62887454204118196e-2, 197289.95555009261, 15337.684217738644, 113.65599575369508, 1085.6508893970181, 5.0988841320189165, -97.460501257372059, -98.460304414813390, -4.9478791660518938, -90.616161433438351, -28.953817235122780, -31.051535137038095, -6.2866141221186399, 10976912.709213279, 906054.78376089458, 11620.339023013812, 108565.08893970181, 504.67383903857071, -9745.8502693947145, -9845.8305442163355, -489.97846233890028, -9061.6161433438351, -2917.2121763139237, -2671.1918928900612, -614.80983636429526]
  μ7 = [-.62887454204118196e-2, -.26012207465436708e-1, -.28953817235122780e-1, -.90616161433438351e-1, -.49478791660518938e-2, -.98500284337856984e-1, -.97500473079758485e-1, .50988841320189165e-2, 1.0856508893970181, .11365599575369508, 8.2329736217213995, 197.28754575851234, -6.2866141221186399, -31.051535137038095, -28.953817235122780, -90.616161433438351, -4.9478791660518938, -98.460304414813390, -97.460501257372059, 5.0988841320189165, 1085.6508893970181, 113.65599575369508, 15337.684217738644, 197289.95555009261, -614.80983636429526, -2671.1918928900612, -2917.2121763139237, -9061.6161433438351, -489.97846233890028, -9845.8305442163355, -9745.8502693947145, 504.67383903857071, 108565.08893970181, 11620.339023013812, 906054.78376089458, 10976912.709213279]

  # case i = 100, j = 200
  λ8 = [396.70077304860189, 17.543356145022943, .23411301648660350, 14.342944230834941, .18526112914933737e-1, 666.51672783140746, 399.85003803142806, .27015902793601083e-2, -.41095405051067463e-1, -.33415716510125988e-1, -.32264187884264475e-1, -.65675877322286981e-2, 396704.82629451725, 33545.920039238453, 234.11301648660350, 14342.944230834941, 18.526112914933737, 649200.71150720687, 393550.82642485962, 2.7015902793601083, -41.095405051067463, -33.415716510125988, -39.142670002576738, -6.5646783716526818, 22067442.846711107, 1937285.2595546954, 23960.940006289331, 1434294.4230834941, 1832.9593262562349, 64835849.281574616, 39324106.975737793, 267.40900733191310, -4109.5405051067463, -3369.3612921885534, -3322.4617386149191, -632.42077145229062]
  μ8 = [-.11567337744728073e-1, -.42254197874274465e-1, -.50657095820470816e-1, -.14040027098949845, -.16529178951409123e-1, -.14996196857194033, -.14993883525920239, -.27504828297471143e-2, .57229945120655245e-1, .10911301648660350, 7.5433561450229433, 196.70077304860189, -11.563428684074576, -49.137672501327363, -50.657095820470816, -140.40027098949845, -16.529178951409123, -149.96097671517856, -149.93784344247571, -2.7504828297471143, 57.229945120655245, 109.11301648660350, 13545.920039238453, 196704.82629451725, -1132.3982724647451, -4321.5625478865746, -5099.4650984169271, -14040.027098949845, -1637.3528974299917, -14996.092712246653, -14993.779385176546, -272.30383163850191, 5722.9945120655245, 11140.427185776511, 826174.14844358434, 10956331.735599996]

  λs = [λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8]
  μs = [μ1 μ2 μ3 μ4 μ5 μ6 μ7 μ8]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, λ, μ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([λ, μ], i, j, t)
      @test x[1] ≈ λs[a, b]
      @test x[2] ≈ μs[a, b]
    end
  end
end


text:
# Datasets are not exported to avoid colliding with names in the Main scope

@doc raw"""
    NILE

The absolute path for the `NILE` dataset stored inside StateSpaceModels.jl.
This dataset provides measurements of the annual flow of the Nile river at Aswan from 1871 to 1970,
in ``10^8 m^3``.

See more on [Nile river annual flow](@ref)

# References
 * Durbin, James, & Siem Jan Koopman. (2012). "Time Series Analysis by State Space Methods: Second Edition." Oxford University Press. (Chapter 2)
"""
const NILE = joinpath(dirname(@__DIR__()), "datasets", "nile.csv")

@doc raw"""
    AIR_PASSENGERS

The absolute path for the `AIR_PASSENGERS` dataset stored inside StateSpaceModels.jl.
This dataset provides monthly totals of a US airline passengers from 1949 to 1960.

See more on [Airline passengers](@ref)

# References
 * https://www.stata-press.com/data/r12/ts.html
"""
const AIR_PASSENGERS = joinpath(dirname(@__DIR__()), "datasets", "airpassengers.csv")

@doc raw"""
    INTERNET

The absolute path for the `INTERNET` dataset stored inside StateSpaceModels.jl.
This dataset provides the number of users logged on to an Internet server each minute over 100 minutes.

# References
 * Durbin, James, & Siem Jan Koopman. (2012). "Time Series Analysis by State Space Methods:
    Second Edition." Oxford University Press. (Chapter 9)
"""
const INTERNET = joinpath(dirname(@__DIR__()), "datasets", "internet.csv")

@doc raw"""
    WHOLESALE_PRICE_INDEX

TODO
"""
const WHOLESALE_PRICE_INDEX = joinpath(
    dirname(@__DIR__()), "datasets", "wholesalepriceindex.csv"
)

@doc raw"""
    VEHICLE_FATALITIES

The absolute path for the `VEHICLE_FATALITIES` dataset stored inside StateSpaceModels.jl.
This dataset provides the number of annual road traffic fatalities in Norway and Finland.

See more on [Finland road traffic fatalities](@ref)

# References
 * Commandeur, Jacques J.F. & Koopman, Siem Jan, 2007. "An Introduction to State Space Time
    Series Analysis," OUP Catalogue, Oxford University Press (Chapter 3)
 * http://staff.feweb.vu.nl/koopman/projects/ckbook/OxCodeAll.zip
"""
const VEHICLE_FATALITIES = joinpath(
    dirname(@__DIR__()), "datasets", "vehicle_fatalities.csv"
)

@doc raw"""
    FRONT_REAR_SEAT_KSI

The absolute path for the `FRONT_REAR_SEAT_KSI` dataset stored inside StateSpaceModels.jl.
This dataset provides the log of british people killed or serious injuried in road accidents accross
UK.

# References
 * Commandeur, Jacques J.F. & Koopman, Siem Jan, 2007. "An Introduction to State Space Time
    Series Analysis," OUP Catalogue, Oxford University Press (Chapter 3)
 * http://staff.feweb.vu.nl/koopman/projects/ckbook/OxCodeAll.zip
"""
const FRONT_REAR_SEAT_KSI = joinpath(dirname(@__DIR__()), "datasets", "front_rear_seat.csv")

@doc raw"""
    RJ_TEMPERATURE

Weekly mean temperature in Rio de Janeiro in Kelvin (K).
"""
const RJ_TEMPERATURE = joinpath(dirname(@__DIR__()), "datasets", "rj_temperature.csv")

@doc raw"""
    US_CHANGE

Percentage changes in quarterly personal consumption expenditure, personal disposable income,
production, savings and the unemployment rate for the US, 1960 to 2016.

Federal Reserve Bank of St Louis.

# References
 * Hyndman, Rob J., Athanasopoulos, George. "Forecasting: Principles and Practice"
"""
const US_CHANGE = joinpath(dirname(@__DIR__()), "datasets", "uschange.csv")

@doc raw"""
    SUNSPOTS_YEAR

Yearly numbers of sunspots from 1700 to 1988 (rounded to one digit).
Federal Reserve Bank of St Louis.

# References
 * H. Tong (1996) Non-Linear Time Series. Clarendon Press, Oxford, p. 471.
"""
const SUNSPOTS_YEAR = joinpath(dirname(@__DIR__()), "datasets", "sunspot_year.csv")


text:
function test_cic( nstep, nbpart)

    @show nstep, nbpart

    nx     = 128
    ny     = 128

    xp  = zeros(Float64, nbpart)
    yp  = zeros(Float64, nbpart)
    vx  = zeros(Float64, (nx+1,ny+1))
    vy  = zeros(Float64, (nx+1,ny+1))

    xmin, xmax = -5.0, 5.0
    ymin, ymax = -5.0, 5.0
    dx = (xmax-xmin)/nx
    dy = (ymax-ymin)/ny
    dt = 0.1

    xp .= rand(nbpart)
    yp .= rand(nbpart)

    for j = 1:ny+1, i = 1:nx+1
        vx[i,j] =   (j-1)*dy + xmin
        vy[i,j] = - (i-1)*dx + ymin
    end

    for istep = 1:nstep

        for k=1:nbpart

           xp[k] = xmin + mod(xp[k]-xmin,xmax - xmin)
           yp[k] = ymin + mod(yp[k]-ymin,ymax - ymin)

           xc = (xp[k] - xmin)/dx
           yc = (yp[k] - ymin)/dy

           i = trunc(Int64,xc)
           j = trunc(Int64,yc)

           dpx = xc - i
           dpy = yc - j

           i += 1
           j += 1

           a1 = (1.0-dpx) * (1.0-dpy)
           a2 = (    dpx) * (1.0-dpy)
           a3 = (    dpx) * (    dpy)
           a4 = (1.0-dpx) * (    dpy)

           vpx = ( a1 * vx[i  ,j  ] + a2 * vx[i+1,j  ]
                 + a3 * vx[i+1,j+1] + a4 * vx[i  ,j+1] )
           vpy = ( a1 * vy[i  ,j  ] + a2 * vy[i+1,j  ]
                 + a3 * vy[i+1,j+1] + a4 * vy[i  ,j+1] )

           xp[k] += vpx * dt
           yp[k] += vpy * dt

        end

        #println(sum(xp)/nbpart, " - ", sum(yp)/nbpart)

    end

end

test_cic( 1, 1) #trigger comoilation
@time test_cic( 1000, 409600)


text:
using ADCME
using PyPlot
using MLDatasets
using ProgressMeter

function encoder(x, n_hidden, n_output, rate)
    local μ, σ
    variable_scope("encoder") do
        y = dense(x, n_hidden, activation = "elu")
        y = dropout(y, rate, ADCME.options.training.training)
        y = dense(y, n_hidden, activation = "tanh")
        y = dropout(y, rate, ADCME.options.training.training)
        y = dense(y, 2n_output)
        μ = y[:, 1:n_output]
        σ = 1e-6 + softplus(y[:,n_output+1:end])
    end
    return μ, σ
end

function decoder(z, n_hidden, n_output, rate)
    local y
    variable_scope("decoder") do
        y = dense(z, n_hidden, activation="tanh")
        y = dropout(y, rate, ADCME.options.training.training)
        y = dense(y, n_hidden, activation="elu")
        y = dropout(y, rate, ADCME.options.training.training)
        y = dense(y, n_output, activation="sigmoid")
    end
    return y
end

function autoencoder(xh, x, dim_img, dim_z, n_hidden, rate)
    μ, σ = encoder(xh, n_hidden, dim_z, rate)
    z = μ + σ .* tf.random_normal(size(μ), 0, 1, dtype=tf.float64)
    y = decoder(z, n_hidden, dim_img, rate)
    y = clip(y, 1e-8, 1-1e-8)

    marginal_likelihood = sum(x .* log(y) + (1-x).*log(1-y), dims=2)
    KL_divergence = 0.5 * sum(μ^2 + σ^2 - log(1e-8 + σ^2) - 1, dims=2)

    marginal_likelihood = mean(marginal_likelihood)
    KL_divergence = mean(KL_divergence)

    ELBO = marginal_likelihood - KL_divergence
    loss = -ELBO
    return y, loss, -marginal_likelihood, KL_divergence
end

function step(epoch)
    tx = train_x[1:batch_size,:]
    @showprogress for i = 1:div(60000, batch_size)
        idx = Array((i-1)*batch_size+1:i*batch_size)
        run(sess, opt, x=>train_x[idx,:])
    end
    y_, loss_, ml_, kl_ = run(sess, [y, loss, ml, KL_divergence],
            feed_dict = Dict(
                ADCME.options.training.training=>false,
                x => tx
            ))
    println("epoch $epoch: L_tot = $(loss_), L_likelihood = $(ml_), L_KL = $(kl_)")

    close("all")
    for i = 1:3
        for j = 1:3
            k = (i-1)*3 + j
            img = reshape(y_[k,:], 28, 28)'|>Array
            subplot(3,3,k)
            imshow(img)
        end
    end
    savefig("result$epoch.png")
end

n_hidden = 500
rate = 0.1
dim_z = 20
dim_img = 28^2
batch_size = 128
ADCME.options.training.training = placeholder(true)
x = placeholder(Float64, shape = [128, 28^2])
xh = x
y, loss, ml, KL_divergence = autoencoder(xh, x, dim_img, dim_z, n_hidden, rate)
opt = AdamOptimizer(1e-3).minimize(loss)

# prepare data
train_x = MNIST.traintensor(Float64);
train_x = Array(reshape(train_x, :, 60000)');

sess = Session(); init(sess)
for i = 1:100
    step(i)
end


text:
#' # Minimizing functions - how a computer learns

#' In the previous notebooks, we've seen that by changing **parameters** in a
#' function, we could find a "best fit" of that function to some data. We use
#' a **loss function** to quantify the "goodness" of a set of parameters and
#' we look for those parameters that **optimize**, in fact **minimize**,
#' the loss function. If we teach a machine how to minimize the loss function
#' on its own, we say that machine is able to **learn** how to model our data.
#'
#' In the last notebook, we threw around terms like **derivative**,
#' **gradient**, and **gradient descent** to give you a rough sense of how we
#' minimize a function on a computer. In this notebook, we will step through
#' these concepts more carefully, with the aim of being able to implement them
#' using Julia.

#' ## Minimizing a 1D function using calculus
#' Let's draw the 1D loss function $L_1$ from a previous notebook again:

using Plots
#gr()
#pyplot()
plotly()

sigma(x) = 1 / (1 + exp(-x))
f(x,w) = sigma(w * x)

x1 = 2
y1 = 0.8

L1(w) = (y1 - f(x1, w))^2

plot(L1, -2, 1.5, xlabel="w", ylabel="L1(w)", legend=false)

#' By eye, we can see that the minimum is around $w=0.6$. But how can we get
#' the computer to work this out on its own?

#' In the previous notebook, we thought of this function plot as a hill,
#' viewed from the side. We could find the minimum by making the hill sticky,
#' and letting a ball roll down it. The ball will find and settle in the
#' minimum of the function.  Now let's see how to teach a computer to do this.

#' We need to find the downhill direction along the hill, which is related to
#' its *slope* (how steep it is). Calculus provides us with tools to calculate
#' that slope!

#' Namely, the slope of a curve $L_1(w)$ at $w$ is given by its **derivative**
#' $L_1'(w)$; geometrically, this is the slope of the **tangent line** to the
#' curve at that point, i.e. the straight line which touches the curve at that
#' point.

#' Calculus provides us with some rules to calculate an analytical formula for
#' the derivative, and we will see later how to apply these rules, albeit
#' indirectly, for machine learning.
#' To gain understanding, however, we will see here how to get the computer
#' to help us out by calculating the derivatives **numerically** instead!

#' ## Approximating derivatives

#' Let's recall that the derivative $L_1'(w)$ of a function is defined as
#'
#' $$L_1'(w) \simeq \frac{L_1(w+h) - L_1(w)}{h},$$
#'
#' for a small step size $h$. (Strictly speaking, we must take the limit
#' when $h$ tends to $0$ to obtain the exact value of the derivative.)

#' #### Exercise 1
#'
#' Write a function to calculate the derivative of a function at a given
#' point. Note that in Julia, we can easily pass functions as arguments
#' to other functions!
#' The function should take the function, the point $w$ at which to calculate
#' the derivative, and the value of $h$, which should have a default value
#' of 0.001.
#'
#' *Note*: A function input argument can have a *default value* if we set
#' the input argument equal to that default value when we define the function.
#' For example,
#'
#' ```julia
#' f(x, a = 3) = a * x^2
#' ```
#'
#' The function `f` will square the input we give it and multiply by `a`.
#' However,  if we choose to call `f(x)` *without* passing it an `a` input,
#' it will assume `a` is `3` and return `3*x^2`.

#' #### Solution
#'
#' We could define the `derivative` function as
function derivative(f, w, h=0.001)
    Ans = f(w + h) - f(w)
    Ans /= h
    return Ans
end

#' **Tests**
@assert isapprox(derivative(x -> x^2, 0.2), 0.4, atol=0.01)
@assert isapprox(derivative(x -> x^2, 0.2, 0.5), 0.9, atol=0.5)
@assert isapprox(derivative(x -> x^3, 0.2), 0.12, atol=0.01)

#' #### Exercise 2
#'
#' Write an interactive visualization of the tangent line to the graph of
#' $L_1$, so that we can visualize the tangent at any point on $L_1$.
#' Include the current value of the derivative in the title.
#'
#' *Hint*: Recall that the straight line through the point $(x_0, y_0)$ with
#' slope $m$ is given by
#'
#' $$\frac{y - y_0}{x - x_0} = m,$$
#'
#' so
#'
#' $$y = y_0 + m*(x - x_0).$$

wmin = -2
wmax = 1
w = 0.5
pl = plot(L1, wmin, wmax, label="L_1")

x0 = w
y0 = L1(w)
scatter!(pl, [x0], [y0], label="")

m = derivative(L1, w)

plot!(pl, x -> y0 + m * (x - x0), ls=:dash, label="tangent line")

ylims!(pl, -0.1, 0.6)
xlims!(pl, wmin, wmax)

#' #### Exercise 3
#'
#' What is the value of the derivative (slope of the tangent line) at a
#' minimum? Can this happen anywhere else?

#' #### Solution
#' The derivative is `0` at a minimum. This happens at any minimum, maximum or
#' "point of inflection" of a smooth function.

#' #### Exercise 4
#'
#' When the derivative $L_1'(w)$ is positive, that means that $L_1$ increases
#' from left to right at point $w$.
#'
#' If the derivative $L_1'(w)$ is positive (> 0), in which direction should
#' we move $w$ to *decrease* $L_1$?

#' #### Solution
#'
#' We should move $w$ left, i.e. decrease it.

#' #### Exercise 5
#'
#' If the derivative $L_1'(w)$ is negative (< 0), in which direction should we
#' move $w$ to *decrease* $L_1$?

#' We can use this information to tell the computer which way to take a step.
# This constitutes the numerical algorithm called **gradient descent**;
#' it is called this since we are descending (moving downwards) along the
#' graph of the function by using information about its gradient (slope).

#' #### Exercise 6
#'
#' Implement gradient descent by following this prescription for an
#' **iterative (repetitive) algorithm**:
#'
#' 1. Start at an initial guess $w_0$ for $w$.
#'
#' 2. At each step, calculate the derivative, $L_1'(w_n)$ at the current value
#' of $w_n$ using the function that you created above.
#'
#' 3. Modify the value of $w$ by a small multiple (for example, $\eta=0.01$)
#' of the value of the derivative you just created, via
#' $w_{n+1} = w_n - \eta L_1'(w_n)$.
#'
#' For this problem, start with $w_0 = -2.0$. Repeat steps 2 and 3 a total
#' of `2000` times.
#'
#' Package this code into a function called `gradient_descent` that takes two
#' inputs, a function and a range for values of $w$, and returns the final
#' value of $w$ and $L_1(w)$.
#'
#' Using `L1` and `-2:0.01:1.5` as your inputs to `gradient_descent`,
#' for what value of $w$ is $L_1$ at a minimum?

#' #### Solution
function gradient_descent(f) # wrange not used?
    w = -2.0
    for i in 1:2000
        fprime = derivative(f, w)
        step = 0.01 * fprime
        w = w - step
    end

    return w, f(w)
end

(w, fw) = gradient_descent(L1)

#' #### Exercise 7
#'
#' Modify your code for gradient descent to return the result once it has
#' found an answer within some *tolerance*, rather than taking a set number
#' of steps. The new prescription for this algorithm is:
#'
#' 1. Start at an initial guess $w_0$ for $w$.
#'
#' 2. At each step, calculate the derivative, $L_1'(w_n)$ at the current value
#' of $w_n$, using the function that you created above.
#'
#' 3. Modify the value of $w$ by a small multiple (for example, $\eta=0.01$)
#' of the value of the derivative you just created, via
#' $w_{n+1} = w_n - \eta L_1'(w_n)$.
#'
#' 4. Check how different $w_{n+1}$ is from $w_n$. If you're satisfied that
#' $L_1(w_{n+1})$ is minimized, return $w_{n+1}$ and $L_1(w_{n+1})$.
#' Otherwise, go to step (2) and continue.
#'
#' Edit `gradient_descent` so that it takes three inputs: a function, a range
#' for values of $w$, and a tolerance that tells you how close $w_{n+1}$ must
#' be to $w_n$ before you can stop iterating.
#'
#' Using `L1`, `-2:0.01:1.5`, and `.000001` as your inputs to
#' `gradient_descent`, for what value of $w$ is $L_1$ at a minimum?

#' #### Solution
#'
#'
#' Your code should look something like this
#'
#' $L_1$ is at a minimum when `gradient_descent(L1, -2:0.01:1.5, .000001)[1]`,
#' or $w$, is around `0.69216`.

function gradient_descent(f, tol)
    w = -2.0
    ccount = 0

    while ccount < 2000000   # do no more than 2000 times
        fprime = derivative(f, w)
        step = 0.01 * fprime
        w = w - step
        if abs(step) < tol
            return w, f(w)
        end

        ccount += 1
    end
end

(w, fw) = gradient_descent(L1, 0.000001)
@assert 0.69215 < w < 0.69217

#' After you've completed this step, you will have found a good approximation
#' to the value $w^*$ where $L_1$ reaches its minimum, and the minimum
#' value $L_1(w^*)$, which in this case is (almost) $0$.

#' #### Exercise 8
#'
#' Alter the function `gradient_descent` so that it stores the results
#' `(w, L1(w))` at each step of the algorithm as an array and returns
#' this array. How many steps does the algorithm take for input parameters
#' `(L1, -2:0.01:1.5, .000001)` before terminating? You should count your
#' starting $w_0$ as your first step.

"""
    gradient_descent(f, wrange, tol)

Use gradient descent method to find the minimum of a function.

# Examples

# Arguments
- `f`: Find the minimum of the function `f`.
- `wrange`: User proivded initial guess range.
- `tol`: Tolerance of consequtive step size.  Returns when the difference of
        consequtive steps is less than `tol`.
"""
function gradient_descent(f, wrange; tol=0.000001, w_init=undef)
    # find initial guess from wrange
    if w_init == undef
        w_vec = Vector{Float64}(undef, 0)
        fw_vec = Vector{Float64}(undef, 0)

        for w = wrange
            fw = f(w)
            push!(w_vec, w)
            push!(fw_vec, fw)
        end

        min_ind = argmin(fw_vec)
        w_init = w_vec[min_ind]
        fw_init = fw_vec[min_ind]
        println("Initial guess: $w_init, $fw_init")
    # return (w_inti, fw_init)
    end

    # gradient descent
    w = w_init
    w_vec = Vector{Float64}(undef, 0)
    fw_vec = Vector{Float64}(undef, 0)
    ccount = 0
    fw = f(w)
    push!(w_vec, w)
    push!(fw_vec, fw)
    while ccount < 1000000
        fprime = derivative(f, w)
        step = 0.01 * fprime

        w = w - step
        fw = f(w)
        push!(w_vec, w)
        push!(fw_vec, fw)
        if abs(step) < tol
            println("iteration = $ccount")
            return w_vec, fw_vec
        end

        ccount += 1
    end

    error("Exceeded maximum iterations")
end
(w, fw) = gradient_descent(L1, -2:0.1:1.5; tol=0.000001)
println("Iteration count: $(length(w))")

#' #### Exercise 9
#'
#' Overlay the steps taken in the last exercise with a plot of
#' $L_1(w)$ vs. $w$.
#'
#' Where does our algorithm take the largest steps?
#' (Where does the ball move fastest down the hill?)
#'
#' A) Between w = -2:-1

#' B) Between w = -1:0

#' C) Between w = 0:.6
#'
#' *Hint*: It may be easier to see what's going on if you only plot,
#' for example, every 15th step.

# Initial guess = 0.7 so the plot range needs to be smaller to see
(w, fw) = gradient_descent(L1, -2:0.1:1.5; tol=0.000001)
pl = plot(L1, 0.692, 0.70, xlabel="w", ylabel="L1(w)", legend=false)
ws = reverse(w[end:-50:1]) # Make sure end point is included
fws = reverse(fw[end:-50:1])
scatter!(pl, ws,fws)

# Initial guess = -2.0
(w, fw) = gradient_descent(L1, -2:0.1:1.5; tol=0.000001, w_init=-2)
pl = plot(L1, -2, 1.5, xlabel="w", ylabel="L1(w)", legend=false)
ws = reverse(w[end:-50:1]) # Make sure end point is included
fws = reverse(fw[end:-50:1])
scatter!(pl, ws,fws)

#' ## Functions of 2 variables and their derivatives
#' So far, we have been looking at the minimizing a loss function
#' $L_1(w)$ that depends on a single parameter, $w$. Now let's turn to
#' the cost function $L_2(w, b)$ from a previous notebook, that is a
#' function of *two* parameters, $w$ and $b$, and try to minimize it.
#'
#' As we've seen, we get a **surface**, instead of a curve,
#' when we graph $L_2$ as a function of both of its parameters.

#' **Exercise 10**
#'
#' Draw a surface plot of $L_2$, given by
#'
#' $$L_{2}(w, b) = \sum_i(y_i - g(x_i, w, b))^2$$
#'
#' using the `surface` function from `Plots.jl`. For this plot, use the
#' values of `xs` and `ys` from notebook 5:
#'
#' We can get a nice interactive 3D plot by using the Plotly backend
#' of `Plots.jl` by executing the command
#'
#'     plotly()
# if !(@isdefined WEAVE_ARGS)
#     plotly()
# end

xs = [2, -3, -1, 1]
ys = [0.8, 0.3, 0.4, 0.4]

sigma(x) = 1/(1 + exp(-x))
g(x,w,b) = sigma(w*x) + b

L2(w,b) = sum(abs2, ys .- g.(xs,w,b))

pl = plot()
surface!(pl, -2:0.02:2, -2:0.02:2, L2, alpha=0.6, xlabel="w",
        ylabel="b", zlabel="L2(w,b)",
        size=(600,600))
#'

#+ echo=false
if !@isdefined WEAVE_ARGS # exe when not from Weave
    gui(pl)
end

#' ### Finding the minimum
#' We can just about see, by rotating the graph, that $L_2$ has a single
#' minimum. We want to find the values of $w$ and $b$ where this
#' minimum value is reached.

#' Following what we did for the function $L_1$ above, we expect that
#' we will need to calculate derivatives of $L_2$. Since the function is
#' more complicated, though, the derivatives are too!
#'
#' It turns out that the right concept is that of the
#' [**gradient**](https://en.wikipedia.org/wiki/Gradient) of
#' $L_2$, denoted $\nabla L_2(w, b)$. This is a **vector** consisting
#' of $2$ numbers if there are $2$ parameters
#'[or, in general, $n$ numbers if there are $n$ parameters].
#'
#' The numbers that form the gradient $\nabla L_2(w, b)$ are called the
#' **partial derivatives** of $L_2$ with respect to $w$ and $b$,  written as
#'
#' $$\frac{\partial L_2}{\partial w} \quad \text{and} \quad \frac{\partial L_2}{\partial b}.$$
#'
#' Although this notation might look complicated, all it means is that we
#' calculate derivatives just like we did before, except that we fix the value
#' of the other variable.
#' For example, to calculate $\frac{\partial L_2}{\partial w}$,
#' the partial derivative of $L_2$ with respect to $w$, we fix the value of
#' $b$ and think of the resulting function as a function of the single
#' variable $w$; then we use the formula for derivatives of functions of a
#' single variable.
#'
#' [Note that $\frac{\partial L_2}{\partial w}$ is itself a function of $w$ and $b$; we could write $\frac{\partial L_2}{\partial w}(w, b)$.]

#' #### Exercise 11
#'
#' Write functions that will allow you to calculate the partial derivatives of
#' a function $f$ of two variables.
#'
#' In particular, declare functions called `partial_w` and `partial_b`.
#' Each should take four inputs - a function $f$ of two variables, the first
#' input argument to $f$, the second input argument to $f$, and a step size
#' `h` with default value `0.001`. `partial_w` should return the partial
#' derivative of $f$ with respect to its first input argument and `partial_b`
#' should return the partial derivative of $f$ with respect to its second
#' input argument.

#' #### Solution
#'
#' The partial derivatives are
partial_w(f, w, b, h=0.001) = (f(w+h,b) - f(w,b))/h
partial_b(f, w, b, h=0.001) = (f(w,b+h) - f(w,b))/h

#' **Test**
partial_x = partial_w
partial_y = partial_b
fn(x,y) = 3*x*y
#' Test default step
@assert isapprox(partial_x(fn, 0.3, 0.3), 0.9; atol=1.e-6)
@assert isapprox(partial_y(fn, 0.3, 0.3), 0.9; atol=1.e-6)
#' Test custom step
fn(x,y) = 3 * x^2 * y^3
@assert isapprox(partial_x(fn, 2.3, 10.2, 1), 17828.2944; atol=1.e-2)
@assert isapprox(partial_y(fn, 2.3, 10.2, 1), 5454.84; atol=1.e-2)

#' #### Exercise 12
#'
#' Use `partial_b` from the last exercise to find the partial derivative of
#' $L_2$ with respect to $w$ at b = 0.3,
#' $\frac{\partial L_2}{\partial w}|_{b = 0.3}$ for `w = -2:0.01:1`

wrange = -2:0.01:1
partialw = partial_w.(L2, wrange, 0.3)

#' #### Exercise 13
#'
#' Plot the cross section of the surface of $L_2(w, b)$ at $b = 0.3$.
#' Make this plot interactive to show that the function `partial_w` gives
#' the slope of the tangent to this cross section for any point `w` in
#' the range `-2:0.01:1`.
#'
#' For what value of $w$ in this range is the slope of the cross section
#' closest to -1?

L2_b(w) = L2(w, 0.3) # b = 0.3

w = 0.0 # Change me
pl = plot(L2_b, wrange) #plot(wrange, L2_b)
scatter!([w], [L2_b(w)])

slope = partial_w(L2, w, 0.3)
wrange_value = slope .* wrange .- slope*w .+ L2.(w,0.3)
plot!(pl, wrange, wrange_value)
ylims!(0.3, 2.0)
title!("Slope = $slope @ $(w)")
#'
#' ## **Optional**: Functions with $n$ inputs
#' If a function $f$ takes $n$ input arguments, we can write them as
#' $p_1, \ldots, p_n$, where $p_i$ means the "$i$th parameter".
#' In Julia, we can wrap them up into a single **vector**. Now we can
#' calculate the partial derivative $\frac{\partial L_2}{\partial p_i}$ with
#' respect to the $i$th variable.

#' #### Exercise 14
#'
#' For the next exercise, you will need to use the splat command, `...`. You can use this command to "open up" a collection and pass all the elements of that collection as inputs to a function.
#'
#' For example, say you have an array, `numbers`,
#'
#' ```julia
#' numbers = [4, 3, 2]
#' ```
#'
#' and you want to use `numbers` to create a $4\times3\times3$ randomly populated array via `rand`. `rand(numbers)` will not do what you want. You could index into `numbers` to grab the values you want and pass them to `rand`, as in
#'
#' ```julia
#' rand(numbers[1], numbers[2], numbers[3])
#' ```
#'
#' or you could use a splat:
#'
#' ```julia
#' rand(numbers...)
#' ```
#'
#' Use `...` to pass the contents of `inputs`
#'
#' ```julia
#' inputs = [30, 12, "cats"]
#' ```
#'
#' to the function `dreams`
#'
#' ```julia
#' dreams(i, j, perfect_mammal) = "I wish I had $(i + j) $perfect_mammal."
#' ```

inputs = [30, 12, "cats"]
dreams(i, j, perfect_mammal) = "I wish I had $(i + j) $perfect_mammal."
dreams(inputs...)

# **Tests**

@assert dreams(inputs...) == "I wish I had 42 cats."

#' #### Exercise 15:
#'
#' Write a function, `partial`, to calculate the $i$th partial derivative of a
#' function. This function should have four inputs
#'
#' * a function, $f$, for which you want to compute the partial derivative
#' * an array, *p*, specifying the values of all input arguments to $f$ at
#'   the point where you want $\frac{\partial f}{\partial p_i}$ computed
#' * the index, $i$, of the variable with respect to which you want to calculate the partial derivative of $f$
#' * a step size with default value 0.001
#'
#' Hint: you will need to `copy` and modify `p` within `partial`.

function partial(fn, p, i, h=0.001)
    p = copy(p)
    a = fn(p...)
    p[i] += 1*h
    b = fn(p...)
    Ans = (b - a)/(1*h)

    return Ans
end

fn1(x,y,z) = x^2 * y * z
fn1(xyz) = fn1(xyz...)
#' $\frac{\partial fn_1}{\partial x}$
fn1px(x,y,z)  = 2*x*y*z
fn1py(x,y,z) = x^2*1*z
fn1pz(x,y,z) = x^2*y*1

x = 1.0
y = 2.0
z = 1.0
xyz = [x, y, z]

actual = partial(fn1, xyz, 1, 1.e-6)
expected = fn1px(xyz...)
@assert isapprox(actual, expected; atol=1.e-4)

using FiniteDifferences
# what the?  Why do I have to use 2nd order central difference for
# first derivate?
actual = central_fdm(2,1)(x->fn1(x,y,z), 1.0)
@assert isapprox(actual, expected; atol=1.e-6)
actual = forward_fdm(2,1)(x->fn1(x,y,z), 1.0)
@assert isapprox(actual, expected; atol=1.e-6)

using FiniteDiff
# This is right
actual = FiniteDiff.finite_difference_derivative(x->fn1(x,y,z), 1.0)
@assert isapprox(actual, expected; atol=1.e-6)

# GradientFiniteDiff.finite_difference_gradient(fn1, xyz)
(gx,gy,gz) = FiniteDiff.finite_difference_gradient(fn1, xyz)

w = 1.1
x = 2.2
y = 3.3
z = 4.4
fn2(w,x,y,z) = w * x^2 * y^3 * z^4
fn2pw(w,x,y,z) = x^2 * y^3 * z^4
fn2px = w * 2 * x * y^3 * z^4
fn2py = w * x^2 * 3 * y^2 * z^4
fn2pz = w * x^2 * y^3 * 4 * z^3

# Test partial()
expected = fn2pw(w,x,y,z)

actual = partial(fn2, [w,x,y,z], 1, 1.e-6)
@assert isapprox(actual, expected; rtol=1.e-6)

actual = FiniteDiff.finite_difference_derivative(
            w->fn2(w,x,y,z), w)
@assert isapprox(actual, expected; rtol=1.e-9)

#' ## Gradient descent in 2 dimensions
#' It turns out that the gradient vector of the function $L_2(w, b)$ gives the
#' direction in the plane $(w, b)$ in which the function $L_2$
#' **increases fastest**.
#'
#' In other words, if you start at the position $(w_0, b_0)$ and take a small
#' step of length $\eta$ to new values $(w', b')$, the value of the function
#' will change to a new value $L_2(w', b')$. How do we get to the minimum of
#' $L_2(w, b)$ as fast as possible? We want to step in the direction where
#' $L_2$ decreases fastest!
#'
#' In multivariable calculus courses, it is shown that $L_2$ will *increase*
#' fastest if you take a step **in the direction of the gradient
#' $\nabla L_2(w, b)$**! To decrease $L_2$ the fastest, we should take a step
#' in the *opposite* direction, $-\nabla L_2(w, b)$.
#'
#' Let's now generalize the gradient descent algorithm that we wrote previously
#' to work our 2-dimensional function.

#' #### Exercise 16
#'
#' Extend your 1D implementation of `gradient_descent` so that it will minimize the function $L_2$.
#'
#' Requirements:
#'
#' * Your new method for `gradient_descent` will take four input arguments: the function $f$ for which you seek the minimum, the range of values for $f$'s first input argument that you will consider, the range of values for $f$'s second input argument that you will consider, and a tolerance that will specify the maximum allowable step size, $\sum_i \eta \frac{\partial f}{\partial p_i}$
#'
#' * Use $\eta = .01$. For example, for a function $f(w, b)$, update $w$ such that $w_{n+1} = w_n - 0.01 * \frac{\partial f}{\partial w_n}$
#'
#' * Seed `gradient_descent` with the starting coordinates [-2.0, -2.0], i.e. $w_0 = -2.0$ and $b_0 = -2.0$.
#'
#' * Return all steps (their coordinates) taken during gradient descent and the values of the loss function at these coordinates.
#'
#' Once you have done this, execute
#'
#' ```julia
#' gradient_descent(L2, -2:0.02:2, -2:0.02:2, .001)
#' ```
#'
#' How many steps were taken by gradient descent?
#'
#' Hint: Do not count your starting coordinates `[-2.0, -2.0]` as a step.

function gradient_descent(fn, wrange, brange; tol=0.001, eta=0.01)
    w = wrange[1]
    b = brange[1]

    w_steps = []
    b_steps = []
    f_steps = []

    while true
        (gfw,gfb) = FiniteDiff.finite_difference_gradient(fn, [w,b])
        w = w - eta * gfw
        b = b - eta * gfb
        f = fn(w,b)

        push!(w_steps, w)
        push!(b_steps, b)
        push!(f_steps, f)

        if abs(eta*gfw + eta*gfb) < tol
            return (w_steps,b_steps,f_steps)
        end

    end

    #return (w_steps, b_steps, f_steps)
end

xs = [2, -3, -1, 1]
ys = [0.8, 0.3, 0.4, 0.4]

sigma(x) = 1/(1 + exp(-x))
g(x,w,b) = sigma(w*x) + b

L2(w,b) = sum(abs2, ys .- g.(xs,w,b))
L2(xb) = L2(xb...)
wrange = -2.0:0.02:2.0
brange = -2.0:0.02:2.0
(w,b,fv) = gradient_descent(L2, wrange, brange)

#' #### Exercise 17
#'
#' Use the `surface` and `scatter!` commands to illustrate the path taken by
#' `gradient_descent` from [-2.0, -2.0] to wherever the algorithm terminates.
#'
#' Where do the scattered points representing the steps of gradient descent
#' appear the most dense?
#'
#' A) Near the starting point at [-2.0, -2.0]

#' B) Near the point [-1.8, 0] where $C_2$ appears nearly flat

#' C) Near the minimum of $C_2$

#' #### Solution
#'
#' B) The following code will plot the function `C2` with an overlay of the
#' steps taken in gradient descent
#'
#' You can change the frequency with which steps are plotted by changing the
#' iterator in the `for` loop from `1:length(values)` to
#' `1:n:length(values)` for varying values of `n`. This will allow you to
#' see more clearly where steps are densely or sparsely packed.

fig = surface(-2:0.02:2, -2:0.02:2, L2, alpha=0.6, xlabel = "w",
             ylabel = "b", zlabel = "C2(w, b)", size=(800,800))
for i in 1:10:length(fv)
   scatter!(fig, [w[i]], [b[i]], [fv[i]], markersize=1)
end
fig