text:
using CDDLib

@testset "Testing Repeated Game functionality" begin

    pd_payoff = [9.0 1.0
                 10.0 3.0]

    A = Player(pd_payoff)
    B = Player(pd_payoff)
    nfg = NormalFormGame((A, B))

    # Tests construction of repeated game
    rpd = RepeatedGame(nfg, 0.75)
    C, H, Z = Games.initialize_sg_hpl(4, [0.0, 0.0], 1.0)

    #
    # Test various helper functions
    #
    @testset "Testing flow utility computations" begin
        @test abs(flow_u_1(rpd, 1, 1) - 9.0) < 1e-14
        @test abs(flow_u_2(rpd, 2, 2) - 3.0) < 1e-14
        @test maximum(abs, flow_u(rpd, 2, 2) - [3.0, 3.0]) < 1e-14
    end

    @testset "Testing best deviation computations" begin
        @test best_dev_i(rpd, 1, 1) == 2
    end

    @testset "Testing unit circle function" begin
        H = Games.unitcircle(4)
        points = [1.0 0.0
                  0.0 1.0
                  -1.0 0.0
                  0.0 -1.0]

        @test maximum(abs, H - points) < 1e-12
    end

    @testset "Testing subgradient and hyperplane level initialize" begin
        C, H, Z = Games.initialize_sg_hpl(4, [0.0, 0.0], 1.0)

        @test maximum(abs, C - ones(4)) < 1e-12
        @test maximum(abs, H - Z') < 1e-12
    end

    @testset "Testing worst value computation" begin
        @test abs(worst_value_i(rpd, H, C, 1) + 1.0) < 1e-12
    end

    #
    # Test the actual computation
    #
    @testset "Testing outer approximation" begin
        dict = Dict(:nH=>64, :maxiter=>150, :tol=>1e-9)
        plib = CDDLib.Library()
        for kwargs in [dict, merge(dict, Dict(:plib=>plib))]
            vertices = @inferred(outerapproximation(rpd; kwargs...))
            p_in_v = [vertices[i, :] for i in 1:size(vertices, 1)]

            mybools = [all(isapprox.([3.0, 3.0], p)) for p in p_in_v]
            @test any(mybools)
        end
    end

end


text:
# Functions that convert between r, c, v and r, c|v, 1 style Associative Arrays

function col2type(A,splitSep)
  # Inverse of val2col (see below for more info)
  r,c,v = find(A)
  cType, cVal = SplitStr(c,splitSep)
  return Assoc(r,Array{Union{AbstractString,Number}}(cType),Array{Union{AbstractString,Number}}(cVal));

end

function val2col(A::Assoc,splitSep)
  #val2col: Append associative array values to column keys; inverse of col2type.
  #Associative array user function.
  #  Usage:
  #    AA = val2col(A,splitSep)
  #  Inputs:
  #    A = associative array with string column keys
  #    splitSep = single character separator
  # Outputs:
  #    A = associative array with string column keys and numeric values

    r, cType, cVal, = find(A);
    c = CatStr(cType,splitSep, cVal);
    return Assoc(r,c,1);
end

ext: jl


text:
"""
`Subject` has following public fields:

* `path` - path to the subject directory
* `identifier` - identifier of subject, extracted from directory name `sub-<subject_id>`
* `sessions` - a vector of `Session`

The `Subject` can be initialized by specifying only `path`. Other behavior can also be
tweaked accordingly with the optional parameters:

* `search::Bool` - defaults to `true`. If `true` this will search over `root` to find
  `Subject`, `Session` and `Files`.
* `load_metadata::Bool` - defaults to `true`. If `true`, this expects a JSON-sidecar to
  be present alongside with every BIDS files.
* `require_modality::Bool` - defaults to `true`. If `true` this expects a modality
  presents in every file name, e.g. sub-subtest_ses-1_run-001_**T1w**.nii.gz
* `longitudinal::Bool` - defaults to `true`. If `true` this expects a session directory
  `ses-` exists under subject directory.
* `strict::Bool` - defaults to `true`. If `true`, `BIDSTools` will throw errors on
  invalid BIDS filenames, this can be turned to `false` to not parse those file names
  and display a warning instead.
* `extract_from_full_path::Bool` - defaults to `true`. If `true`, will try to extract
  `sub` and `ses` from full path and append to `entities` if they can't be found in the
  parsed filename.

# Example

```julia-repl
julia> sub = Subject("/path/to/bids/root/sub-Subtest/")
Subject:
    identifier = Subtest
    total session = 1
    total files = 8
```
"""
struct Subject
    path::String
    identifier::String
    sessions::Vector{Session}
end

function Subject(
    path::AbstractString;
    search::Bool=true,
    load_metadata::Bool=true,
    require_modality::Bool=true,
    longitudinal::Bool=true,
    strict::Bool=true,
    extract_from_full_path::Bool=true
)
    sessions = Vector{Session}()
    sep = Vector{Char}(Base.Filesystem.path_separator)
    @assert startswith(basename(rstrip(path, sep)), "sub-")
    identifier = replace(basename(rstrip(path, sep)), "sub-" => "")
    if search
        for d in readdir(path)
            if isdir(joinpath(path, d)) && startswith(d, "ses-")
                push!(
                    sessions,
                    Session(
                        joinpath(path, d),
                        search=search,
                        load_metadata=load_metadata,
                        require_modality=require_modality,
                        longitudinal=longitudinal,
                        strict=strict,
                        extract_from_full_path=extract_from_full_path
                    )
                )

            end
        end
    end
    Subject(path, identifier, sessions)
end

#-------------------------------------------------------------------------------

"""
    function total_sessions(subject::Subject)

Get number of session of a subject. The argument `Subject` can be changed into `Layout`.
"""
function total_sessions(subject::Subject)
    length(subject.sessions)
end

# Docstring in Session.jl
function total_files(subject::Subject)
    sum(total_files, subject.sessions)
end

function Base.show(io::IO, subject::Subject)
    print(io, """
        Subject:
            identifier = $(subject.identifier)
            total session = $(total_sessions(subject))
            total files = $(total_files(subject))""")
end

# Docstring of get_files is in File.jl
function get_files(
    subject::Subject; path::Union{String, Regex, Nothing}=nothing, kws...
)
    result = Vector{File}()
    for ses in subject.sessions
        files = get_files(ses; path=path, kws...)
        push!(result, files...)
    end
    return result
end


text:
# Everytime I encounter something that raises some error and I fix it, I should add that
# specific problem to this list to make sure, that this specific run then works without
# bugs.
using ProbNumDiffEq
using Test
using LinearAlgebra
using UnPack
using ParameterizedFunctions
using OrdinaryDiffEq
using Plots

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary:
    prob_ode_fitzhughnagumo, prob_ode_vanstiff, prob_ode_2Dlinear, prob_ode_linear

@testset "Problem with analytic solution" begin
    linear(u, p, t) = p .* u
    linear_analytic(u0, p, t) = @. u0 * exp(p * t)
    prob =
        ODEProblem(ODEFunction(linear, analytic=linear_analytic), [1 / 2], (0.0, 1.0), 1.01)

    @test solve(prob, EK0()) isa ProbNumDiffEq.ProbODESolution
    @test solve(ProbNumDiffEq.remake_prob_with_jac(prob), EK1()) isa
          ProbNumDiffEq.ProbODESolution
    sol = solve(prob, EK0())
    @test sol.errors isa Dict
    @test all(haskey.(Ref(sol.errors), (:l∞, :l2, :final)))
end

@testset "Matrix-Valued Problem" begin
    prob = prob_ode_2Dlinear
    my_p = similar(prob.u0)
    my_p .= 1.01
    prob = remake(prob, p=my_p)

    @testset "$alg" for alg in [EK0(), EK1()]
        @test solve(prob, alg) isa ProbNumDiffEq.ProbODESolution
        sol = solve(prob, alg)

        @test length(sol.u[1]) == length(sol.pu.μ[1])
        @test sol.u[1][:] == sol.pu.μ[1]
        @test sol.u ≈ sol.u_analytic rtol = 1e-4
        @test plot(sol) isa AbstractPlot
    end
end

@testset "scalar-valued problem" begin
    @testset "$alg" for alg in [EK0(), EK1()]
        @test_broken solve(prob_ode_linear, alg) isa ProbNumDiffEq.ProbODESolution
    end
end

@testset "Stiff Vanderpol" begin
    prob = ProbNumDiffEq.remake_prob_with_jac(prob_ode_vanstiff)
    @test solve(prob, EK1(order=3)) isa ProbNumDiffEq.ProbODESolution
end

@testset "Big Float" begin
    prob = prob_ode_fitzhughnagumo
    prob = remake(prob, u0=big.(prob.u0))
    sol = solve(prob, EK0(order=3))
    @test eltype(eltype(sol.u)) == BigFloat
    @test eltype(eltype(sol.pu.μ)) == BigFloat
    @test eltype(eltype(sol.pu.Σ)) == BigFloat
    @test sol isa ProbNumDiffEq.ProbODESolution
end

@testset "IIP problem" begin
    f(du, u, p, t) = (du[1] = p[1] * u[1] .* (1 .- u[1]))
    prob = ODEProblem(f, [1e-1], (0.0, 5), [3.0])
    @testset "without jacobian" begin
        # first without defined jac
        @test solve(prob, EK0(order=4)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1(order=4)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=1)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=2)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=3)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK0(initialization=RungeKuttaInit())) isa
              ProbNumDiffEq.ProbODESolution
    end
    @testset "with jacobian" begin
        # now with defined jac
        prob = ProbNumDiffEq.remake_prob_with_jac(prob)
        @test solve(prob, EK1(order=4)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=1)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=2)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=3)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK0(initialization=RungeKuttaInit())) isa
              ProbNumDiffEq.ProbODESolution
    end
end

@testset "OOP problem" begin
    f(u, p, t) = ([p[1] * u[1] .* (1 .- u[1])])
    prob = ODEProblem(f, [1e-1], (0.0, 5), [3.0])
    @testset "without jacobian" begin
        # first without defined jac
        @test solve(prob, EK0(order=4)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1(order=4)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=1)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=2)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=3)) isa ProbNumDiffEq.ProbODESolution
    end
    @testset "with jacobian" begin
        # now with defined jac
        prob = ProbNumDiffEq.remake_prob_with_jac(prob)
        @test solve(prob, EK0(order=4)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1(order=4)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=1)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=2)) isa ProbNumDiffEq.ProbODESolution
        @test solve(prob, EK1FDB(order=4, jac_quality=3)) isa ProbNumDiffEq.ProbODESolution
    end
end

@testset "Callback: Harmonic Oscillator with condition on E=2" begin
    u0 = ones(2)
    function harmonic_oscillator(du, u, p, t)
        du[1] = u[2]
        return du[2] = -u[1]
    end
    prob = ODEProblem(harmonic_oscillator, u0, (0.0, 100.0))

    function Callback()
        function affect!(integ)
            @unpack x_filt, Proj, PI, E0 = integ.cache

            x = x_filt

            m, P = x.μ, x.Σ

            m0, P0 = E0 * m, ProbNumDiffEq.X_A_Xt(P, E0)

            e = m0'm0
            H = 2m0'E0
            S = H * P * H'

            S_inv = inv(S)
            K = P * H' * S_inv

            mnew = m + K * (2 .- e)
            Pnew = ProbNumDiffEq.X_A_Xt(P, (I - K * H)) # + X_A_Xt(R, K)

            # @info m P e S K mnew
            copy!(m, mnew)
            return copy!(P, Pnew)
        end
        condtion = (t, u, integrator) -> true
        save_positions = (true, true)
        return DiscreteCallback(condtion, affect!, save_positions=save_positions)
    end

    @test solve(prob, EK0(order=3)) isa ProbNumDiffEq.ProbODESolution
    @test solve(prob, EK0(order=3), callback=Callback()) isa ProbNumDiffEq.ProbODESolution
end

@testset "ManifoldUpdate callback test" begin
    # Again: Harmonic Oscillator with condition on E=2
    u0 = ones(2)
    function harmonic_oscillator(du, u, p, t)
        du[1] = u[2]
        return du[2] = -u[1]
    end
    prob = ODEProblem(harmonic_oscillator, u0, (0.0, 100.0))
    appxsol = solve(prob, Vern9(), abstol=1e-12, reltol=1e-12)

    E(u) = [dot(u, u) - 2]

    sol1 = solve(prob, EK0(order=3))
    @test sol1 isa ProbNumDiffEq.ProbODESolution
    sol2 = solve(prob, EK0(order=3), callback=ManifoldUpdate(E))
    @test sol2 isa ProbNumDiffEq.ProbODESolution

    @test E(sol1[end]) .^ 2 > E(sol2[end]) .^ 2

    err1 = sol1[end] .- appxsol[end]
    err2 = sol2[end] .- appxsol[end]
    @test all(err1 .^ 2 > err2 .^ 2)
end

@testset "Problem definition with ParameterizedFunctions.jl" begin
    f = @ode_def LotkaVolterra begin
        dx = a * x - b * x * y
        dy = -c * y + d * x * y
    end a b c d
    p = [1.5, 1, 3, 1]
    tspan = (0.0, 10.0)
    u0 = [1.0, 1.0]
    prob = ODEProblem(f, u0, tspan, p)
    @test solve(prob, EK1(order=3)) isa ProbNumDiffEq.ProbODESolution
end

@testset "Robertson in mass-matrix-ODE form" begin
    function rober(du, u, p, t)
        y₁, y₂, y₃ = u
        k₁, k₂, k₃ = p
        du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
        du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
        du[3] = y₁ + y₂ + y₃ - 1
        return nothing
    end
    M = [
        1.0 0 0
        0 1.0 0
        0 0 0
    ]
    f = ODEFunction(rober, mass_matrix=M)
    prob = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e1), (0.04, 3e7, 1e4))

    sol1 = solve(prob, EK1(order=3))
    sol2 = solve(prob, RadauIIA5())
    @test sol1[end] ≈ sol2[end] rtol = 1e-5
end


text:
@testset "CST" begin
  naca = NACA"0015"
  NP = 65
  x0, y0 = gen_airfoil(naca, NP, orient=:CCW)
  # original airfoil
  cst = fit(:CST, x0, y0, N=6)
  xtest = x0[NP:end]
  yuaf = y_upper(cst, xtest)
  ylaf = y_lower(cst, xtest)
  @test cst.a0 > 0
  @test all(isapprox.(y0[NP:-1:1], yuaf, atol=2e-4))
  @test all(isapprox.(y0[NP:end], ylaf, atol=2e-4))

  dynaca = dy_upper(cst, xtest)
  dycst = dy_upper(cst, xtest)
  @test all(isapprox.(dynaca, dycst, rtol=.01, atol=.01))
  dynaca = dy_lower(cst, xtest)
  dycst = dy_lower(cst, xtest)
  @test all(isapprox.(dynaca, dycst, rtol=.01, atol=.01))
  cstr = fit(:CST, reverse(x0), reverse(y0), N=6)
  @test get_param(cst) ≈ get_param(cstr)
  # fit cambered airfoil
  cst = CST(undef, Float64, 5)
  a = get_param(cst)
  set_param!(cst, a)
  naca = NACA"4412"
  x0, y0 = gen_airfoil(naca, NP)
  fit!(cst, x0, y0)
  xx, yy = gen_airfoil(cst, NP)
  @test xx[end] == xx[1] == 1
  @test issorted(xx[NP:-1:1]) && issorted(xx[NP:end])
  @test y_upper(cst, 1.0) > y_lower(cst, 1.0)
  modify_te!(cst, T_TE=0) # modify to sharp TE
  @test isapprox(y_upper(cst, 1.0), y_lower(cst, 1.0), atol=1e-6)
  end


text:
@testset "115.distinct-subsequences.jl" begin
    @test num_distinct("rabbbit", "rabbit") == 3
    @test num_distinct("babgbag", "bag") == 5
    @test num_distinct("", "") == 1
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/nonlinear.jl
# Test general nonlinear
#############################################################################
using Compat, Compat.LinearAlgebra, Compat.SparseArrays, Compat.Test
using MathProgBase, JuMP
# If solvers not loaded, load them (i.e running just these tests)
!isdef(:nlp_solvers) && error("You must declare the solvers to test")

mutable struct DummyNLPSolver <: MathProgBase.AbstractMathProgSolver
end
mutable struct DummyNLPModel <: MathProgBase.AbstractNonlinearModel
end

@testset "Nonlinear" begin

    @testset "getvalue on arrays" begin
        m = Model()
        @variable(m, x, start = π/2)
        @NLexpression(m, f1, sin(x))
        @NLexpression(m, f2, sin(2x))
        @NLexpression(m, f3[i=1:2], sin(i*x))

        @test isapprox(getvalue(f1), 1, atol=1e-5)
        @test isapprox(getvalue(f2), 0, atol=1e-5)

        @test getvalue([f1, f2]) == getvalue(f3)

        v = [1.0, 2.0]
        @NLparameter(m, vparam[i=1:2] == v[i])
        @test getvalue(vparam) == v
        v[1] = 3.0
        setvalue(vparam, v)
        @test getvalue(vparam) == v
    end

    @testset "HS071 solves correctly with $nlp_solver" for nlp_solver in nlp_solvers
        # hs071
        # Polynomial objective and constraints
        # min x1 * x4 * (x1 + x2 + x3) + x3
        # st  x1 * x2 * x3 * x4 >= 25
        #     x1^2 + x2^2 + x3^2 + x4^2 = 40
        #     1 <= x1, x2, x3, x4 <= 5
        # Start at (1,5,5,1)
        # End at (1.000..., 4.743..., 3.821..., 1.379...)
        m = Model(solver=nlp_solver)
        initval = [1,5,5,1]
        @variable(m, 1 <= x[i=1:4] <= 5, start=initval[i])
        @NLobjective(m, Min, x[1]*x[4]*(x[1]+x[2]+x[3]) + x[3])
        @NLconstraint(m, x[1]*x[2]*x[3]*x[4] >= 25)
        @NLconstraint(m, sum(x[i]^2 for i=1:4) == 40)
        @test JuMP.numnlconstr(m) == 2
        @test MathProgBase.numconstr(m) == 2
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getvalue(x),
            [1.000000, 4.742999, 3.821150, 1.379408], atol=1e-3)
    end

    @testset "HS071 solves correctly (no macros) with $nlp_solver" for nlp_solver in nlp_solvers
        # hs071
        # Polynomial objective and constraints
        # min x1 * x4 * (x1 + x2 + x3) + x3
        # st  x1 * x2 * x3 * x4 >= 25
        #     x1^2 + x2^2 + x3^2 + x4^2 = 40
        #     1 <= x1, x2, x3, x4 <= 5
        # Start at (1,5,5,1)
        # End at (1.000..., 4.743..., 3.821..., 1.379...)
        m = Model(solver=nlp_solver)
        initval = [1,5,5,1]
        @variable(m, 1 <= x[i=1:4] <= 5, start=initval[i])
        JuMP.setNLobjective(m, :Min, :($(x[1])*$(x[4])*($(x[1])+$(x[2])+$(x[3])) + $(x[3])))
        JuMP.addNLconstraint(m, :($(x[1])*$(x[2])*$(x[3])*$(x[4]) >= 25))
        JuMP.addNLconstraint(m, :($(x[1])^2+$(x[2])^2+$(x[3])^2+$(x[4])^2 == 40))
        @test MathProgBase.numconstr(m) == 2
        @test_throws ErrorException JuMP.addNLconstraint(m, :(x[1]^2+x[2]^2+x[3]^2+x[4]^2 == 40))
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getvalue(x),
            [1.000000, 4.742999, 3.821150, 1.379408], atol=1e-3)
    end

    @testset "HS071 solves correctly (epigraph) with $nlp_solver" for nlp_solver in nlp_solvers
        # hs071, with epigraph formulation
        # Linear objective, nonlinear constraints
        # min t
        # st  t >= x1 * x4 * (x1 + x2 + x3) + x3
        #     ...
        m = Model(solver=nlp_solver)
        start = [1.0, 5.0, 5.0, 1.0]
        @variable(m, 1 <= x[i=1:4] <= 5, start = start[i])
        @variable(m, t, start = 100)
        @objective(m, Min, t)
        @NLconstraint(m, t >= x[1]*x[4]*(x[1]+x[2]+x[3]) + x[3])
        @NLconstraint(m, x[1]*x[2]*x[3]*x[4] >= 25)
        @NLconstraint(m, sum{x[i]^2,i=1:4} == 40)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getvalue(x),
            [1.000000, 4.742999, 3.821150, 1.379408], atol=1e-3)
    end

    @testset "ifelse with $nlp_solver" for nlp_solver in nlp_solvers
        (occursin("OsilSolver", "$(typeof(nlp_solver))") || occursin("NLoptSolver", "$(typeof(nlp_solver))") || occursin("BaronSolver", "$(typeof(nlp_solver))")) && continue
        # XXX baron doesn't support ifelse
        occursin("baron", nlp_solver.solver_name) && continue
        m = Model(solver=nlp_solver)
        @variable(m, x, start = 2)
        # minimizer at smooth point, solvers should be okay
        @NLobjective(m, Min, ifelse( x <= 1, x^2, x) )
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getvalue(x), 0.0, atol=1e-5)
    end

    @testset "Accepting fixed variables with $nlp_solver (simplify = $simplify)" for nlp_solver in convex_nlp_solvers, simplify in [true,false]
        m = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
        @variable(m, x == 0)
        @variable(m, y ≥ 0)
        @objective(m, Min, y)
        @NLconstraint(m, y ≥ x^2)
        for α in 1:4
            JuMP.fix(x, α)
            solve(m)
            @test isapprox(getvalue(y), α^2, atol=1e-6)
        end
    end

    @testset "QP solve through NL pathway with $nlp_solver" for nlp_solver in nlp_solvers
        # Solve a problem with quadratic objective with linear
        # constraints, but force it to use the nonlinear code.
        m = Model(solver=nlp_solver)
        @variable(m, 0.5 <= x <=  2)
        @variable(m, 0.0 <= y <= 30)
        @NLparameter(m, param == 1.0)
        @objective(m, Min, (x+y)^2)
        @NLconstraint(m, x + y >= param)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), 1.0, atol=1e-6)
        @NLexpression(m, lhs, x+y)
        @test isapprox(getvalue(x)+getvalue(y), 1.0, atol=1e-6)
        @test isapprox(getvalue(lhs), 1.0, atol=1e-6)

        setvalue(param,10)
        @test m.internalModelLoaded == true
        status = solve(m)
        @test isapprox(getobjectivevalue(m), 10.0^2, atol=1e-6)
        @test isapprox(getvalue(x)+getvalue(y), 10.0, atol=1e-6)
        @test isapprox(getvalue(lhs), 10.0, atol=1e-6)
    end

    @testset "Quad con solve through NL pathway" for nlp_solver in convex_nlp_solvers
        # Solve a problem with linear objective with quadratic
        # constraints, but force it to use the nonlinear code.
        m = Model(solver=nlp_solver)
        @variable(m, -2 <= x <= 2)
        @variable(m, -2 <= y <= 2)
        @NLobjective(m, Min, x - y)
        @constraint(m, x + x^2 + x*y + y^2 <= 1)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), -1-4/sqrt(3), atol=1e-6)
        @test isapprox(getvalue(x) + getvalue(y), -1/3, atol=1e-3)
    end

    @testset "Resolve with parameter with $nlp_solver (simplify = $simplify)" for nlp_solver in convex_nlp_solvers, simplify in [true,false]
        m = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
        @variable(m, z)
        @NLparameter(m, x == 1.0)
        @NLobjective(m, Min, (z-x)^2)
        status = solve(m)
        @test status == :Optimal
        @test isapprox(getvalue(z), 1.0, atol=1e-3)

        setvalue(x, 5.0)
        status = solve(m)
        @test status == :Optimal
        @test isapprox(getvalue(z), 5.0, atol=1e-3)
    end

    @testset "Two-sided nonlinear constraints with $nlp_solver" for nlp_solver in convex_nlp_solvers
        # TODO support that
        println("Doubly bounded constaints are not yet supported")
        continue
        m = Model(solver=nlp_solver)
        @variable(m, x)
        @NLobjective(m, Max, x)
        l = -1
        u = 1
        @NLconstraint(m, l <= x <= u)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), u, atol=1e-6)

        @NLobjective(m, Min, x)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), l, atol=1e-6)
    end

    @testset "Two-sided nonlinear constraints (no macros) with $nlp_solver" for nlp_solver in convex_nlp_solvers
        println("Doubly bounded constaints are not yet supported")
        continue
        m = Model(solver=nlp_solver)
        @variable(m, x)
        JuMP.setNLobjective(m, :Max, x)
        l = -1
        u = 1
        JuMP.addNLconstraint(m, :($l <= $x <= $u))
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), u, atol=1e-6)

        JuMP.setNLobjective(m, :Min, x)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), l, atol=1e-6)
    end

    @testset "Quadratic equality constraints with $nlp_solver" for nlp_solver in nlp_solvers
        m = Model(solver=nlp_solver)
        @variable(m, 0 <= x[1:2] <= 1)
        @constraint(m, x[1]^2 + x[2]^2 == 1/2)
        @NLobjective(m, Max, x[1] - x[2])
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getvalue(x), [sqrt(1/2), 0], atol=1e-6)
    end

    @testset "Mixed integer nonlinear problems with $minlp_solver" for minlp_solver in minlp_solvers
        # Solve test problem 1 (Synthesis of processing system) in
        # M. Duran & I.E. Grossmann, "An outer approximation algorithm for
        # a class of mixed integer nonlinear programs", Mathematical
        # Programming 36, pp. 307-339, 1986.  The problem also appears as
        # problem synthes1 in the MacMINLP test set.
        m = Model(solver=minlp_solver)
        x_U = [2,2,1]
        @variable(m, x_U[i] >= x[i=1:3] >= 0)
        @variable(m, y[4:6], Bin)
        @NLobjective(m, Min, 10 + 10*x[1] - 7*x[3] + 5*y[4] + 6*y[5] + 8*y[6] - 18*log(x[2]+1) - 19.2*log(x[1]-x[2]+1))
        @NLconstraints(m, begin
            0.8*log(x[2] + 1) + 0.96*log(x[1] - x[2] + 1) - 0.8*x[3] >= 0
            log(x[2] + 1) + 1.2*log(x[1] - x[2] + 1) - x[3] - 2*y[6] >= -2
            x[2] - x[1] <= 0
            x[2] - 2*y[4] <= 0
            x[1] - x[2] - 2*y[5] <= 0
            y[4] + y[5] <= 1
        end)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), 6.00976, atol=1e-5)
        @test isapprox(getvalue(x), [1.30098, 0.0, 1.0], atol=1e-5)
        @test isapprox(getvalue(y)[:], [0.0, 1.0, 0.0], atol=1e-5)
    end

    @testset "Continuous relaxation of minlp test problem with $nlp_solver" for nlp_solver in nlp_solvers
        # Solve continuous relaxation of test problem 1 (Synthesis of processing system) in
        # M. Duran & I.E. Grossmann, "An outer approximation algorithm for
        # a class of mixed integer nonlinear programs", Mathematical
        # Programming 36, pp. 307-339, 1986.  The problem also appears as
        # problem synthes1 in the MacMINLP test set.
        # Introduce auxiliary nonnegative variable for the x[1]-x[2]+1 term
        m = Model(solver=nlp_solver)
        x_U = [2,2,1]
        @variable(m, x_U[i] >= x[i=1:3] >= 0)
        @variable(m, 1 >= y[4:6] >= 0)
        @variable(m, z >= 0, start=1)
        @NLobjective(m, Min, 10 + 10*x[1] - 7*x[3] + 5*y[4] + 6*y[5] + 8*y[6] - 18*log(x[2]+1) - 19.2*log(z))
        @NLconstraints(m, begin
            0.8*log(x[2] + 1) + 0.96*log(z) - 0.8*x[3] >= 0
            log(x[2] + 1) + 1.2*log(z) - x[3] - 2*y[6] >= -2
            x[2] - x[1] <= 0
            x[2] - 2*y[4] <= 0
            x[1] - x[2] - 2*y[5] <= 0
            y[4] + y[5] <= 1
            x[1] - x[2] + 1 == z
        end)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getobjectivevalue(m), 0.7593, atol=5e-5)
        @test isapprox(getvalue(x), [1.1465, 0.54645, 1.0], atol=2e-4)
        @test isapprox(getvalue(y)[:], [0.2732, 0.3, 0.0], atol=2e-4)
        @test isapprox(getvalue(z), 1.6, atol=2e-4)
    end

    @testset "Maximization objective with $nlp_solver" for nlp_solver in convex_nlp_solvers
        # Solve a simple problem with a maximization objective
        # XXX baron does seem to be able to solve it
        m = Model(solver=nlp_solver)
        @variable(m, -2 <= x <= 2); setvalue(x, -1.8)
        @variable(m, -2 <= y <= 2); setvalue(y,  1.5)
        @NLobjective(m, Max, y - x)
        @constraint(m, x + x^2 + x*y + y^2 <= 1)

        @test solve(m) == :Optimal
        @test isapprox(getobjectivevalue(m), 1+4/sqrt(3), atol=1e-6)
        @test isapprox(getvalue(x) + getvalue(y), -1/3, atol=1e-3)
    end

    @testset "Maximization objective (embedded expressions) with $nlp_solver (simplify = $simplify)" for nlp_solver in convex_nlp_solvers, simplify in [true,false]
        # XXX baron does seem to be able to solve it
        m = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
        @variable(m, -2 <= x <= 2); setvalue(x, -1.8)
        @variable(m, -2 <= y <= 2); setvalue(y,  1.5)
        @NLobjective(m, Max, y - x)
        @NLexpression(m, quadexpr, x + x^2 + x*y + y^2)
        @NLconstraint(m, quadexpr <= 1)

        @test solve(m) == :Optimal
        @test isapprox(getobjectivevalue(m), 1+4/sqrt(3), atol=1e-6)
        @test isapprox(getvalue(x) + getvalue(y), -1/3, atol=1e-3)
        @test isapprox(getvalue(quadexpr), 1, atol=1e-5)
        @NLexpression(m, quadexpr2, x + x^2 + x*y + y^2)
        @test isapprox(getvalue(quadexpr2), 1, atol=1e-5)
    end

    @testset "Infeasibility detection with $nlp_solver" for nlp_solver in convex_nlp_solvers
        occursin("NLoptSolver", string(typeof(nlp_solver))) && continue
        occursin("MosekSolver", string(typeof(nlp_solver))) && continue
        # (Attempt to) solve an infeasible problem
        m = Model(solver=nlp_solver)
        n = 10
        @variable(m, 0 <= x[i=1:n] <= 1)
        @NLobjective(m, Max, x[n])
        for i in 1:n-1
            @NLconstraint(m, x[i+1]-x[i] == 0.15)
        end
        @test solve(m, suppress_warnings=true) == :Infeasible
    end

    @testset "Unboundedness detection with $nlp_solver" for nlp_solver in convex_nlp_solvers
        occursin("NLoptSolver", string(typeof(nlp_solver))) && continue
        # (Attempt to) solve an unbounded problem
        m = Model(solver=nlp_solver)
        @variable(m, x >= 0)
        @NLobjective(m, Max, x)
        @NLconstraint(m, x >= 5)
        @test solve(m, suppress_warnings=true) == :Unbounded
    end

    @testset "Entropy maximization with $nlp_solver" for nlp_solver in convex_nlp_solvers
        (occursin("conopt", nlp_solver.solver_name) || occursin("snopt", nlp_solver.solver_name) || nlp_solver.solver_name == "" || nlp_solver.solver_name == "minos") && continue
        m = Model(solver=nlp_solver)
        N = 3
        @variable(m, x[1:N] >= 0, start = 1)
        @NLexpression(m, entropy[i=1:N], -x[i]*log(x[i]))
        @NLobjective(m, Max, sum(entropy[i] for i = 1:N))
        @constraint(m, sum(x) == 1)

        @test solve(m) == :Optimal
        @test isapprox(getvalue(x), [1/3,1/3,1/3], atol=1e-3)
    end

    @testset "Entropy maximization (reformulation) with $nlp_solver" for nlp_solver in convex_nlp_solvers
        (occursin("conopt", nlp_solver.solver_name) || occursin("snopt", nlp_solver.solver_name) || nlp_solver.solver_name == "" || nlp_solver.solver_name == "minos") && continue

        m = Model(solver=nlp_solver)
        idx = [1,2,3,4]
        @variable(m, x[idx] >= 0, start = 1)
        @variable(m, z[1:4], start = 0)
        @NLexpression(m, entropy[i=idx], -x[i]*log(x[i]))
        @NLobjective(m, Max, sum(z[i] for i = 1:2) + sum(z[i]/2 for i=3:4))
        @NLconstraint(m, z_constr1[i=1], z[i] <= entropy[i])
        @NLconstraint(m, z_constr1_dup[i=2], z[i] <= entropy[i]) # duplicate expressions
        @NLconstraint(m, z_constr2[i=3:4], z[i] <= 2*entropy[i])
        @constraint(m, sum(x) == 1)

        @test solve(m) == :Optimal
        @test isapprox([getvalue(x[i]) for i in idx], [1/4,1/4,1/4,1/4], atol=1e-3)
        @test isapprox(getvalue(entropy[1]), -(1/4)*log(1/4), atol=1e-4)
        @NLexpression(m, zexpr[i=1:4], z[i])
        @test isapprox(getvalue(zexpr[1]), -(1/4)*log(1/4), atol=1e-4)
    end

    @testset "Derivatives of x^4, x < 0 with $nlp_solver" for nlp_solver in convex_nlp_solvers
        m = Model(solver=nlp_solver)
        @variable(m, x >= -1, start = -0.5)
        @NLobjective(m, Min, x^4)
        status = solve(m)

        @test status == :Optimal
        @test isapprox(getvalue(x), 0.0, atol=1e-2)
    end

    @testset "Nonlinear duals with $nlp_solver (simplify = $simplify)" for nlp_solver in nlp_solvers, simplify in [true,false]
        applicable(MathProgBase.getconstrduals, MathProgBase.NonlinearModel(nlp_solver)) || continue
        modA = Model(solver=nlp_solver, simplify_nonlinear_expressions=simplify)
        @variable(modA, x >= 0)
        @variable(modA, y <= 5)
        @variable(modA, 2 <= z <= 4)
        @variable(modA, 0 <= r[i=3:6] <= i)
        @NLobjective(modA, Min, -((x + y)/2.0 + 3.0)/3.0 - z - r[3])
        @constraint(modA, cons1, x+y >= 2)
        @constraint(modA, cons2, sum(r[i] for i=3:5) <= (2 - x)/2.0)
        @NLconstraint(modA, cons3, 7.0*y <= z + r[6]/1.9)

        # Getter/setters
        @test MathProgBase.numconstr(modA) == 3
        @test JuMP.numnlconstr(modA) == 1

        # Solution
        @test solve(modA) == :Optimal
        @test isapprox(getobjectivevalue(modA), -5.8446115, atol=1e-6)
        @test isapprox(getvalue(x), 0.9774436, atol=1e-6)
        @test isapprox(getvalue(y), 1.0225563, atol=1e-6)
        @test isapprox(getvalue(z), 4.0, atol=1e-6)
        @test isapprox(getvalue(r)[3], 0.5112781, atol=1e-6)
        @test isapprox(getvalue(r)[4], 0.0, atol=1e-6)
        @test isapprox(getvalue(r)[5], 0.0, atol=1e-6)
        # knitro has issue with atol=1e-6 :(
        @test isapprox(getvalue(r)[6], 6.0, atol=1e-5)

        # Reduced costs
        @test isapprox(getdual(x), 0.0, atol=1e-6)
        @test isapprox(getdual(y), 0.0, atol=1e-6)
        @test isapprox(getdual(z), -1.0714286, atol=1e-6)
        @test isapprox(getdual(r)[3], 0.0, atol=1e-6)
        @test isapprox(getdual(r)[4], 1.0, atol=1e-6)
        @test isapprox(getdual(r)[5], 1.0, atol=1e-6)
        @test isapprox(getdual(r)[6], -0.03759398, atol=1e-6)

        # Row duals
        @test isapprox(getdual(cons1), 0.333333, atol=1e-6)
        @test isapprox(getdual(cons2), -1.0, atol=1e-6)
        @test isapprox(getdual(cons3), -0.0714286, atol=1e-6)
    end

    @testset "Nonlinear duals (Max) with $nlp_solver" for nlp_solver in nlp_solvers
        applicable(MathProgBase.getconstrduals, MathProgBase.NonlinearModel(nlp_solver)) || continue
        modA = Model(solver=nlp_solver)
        @variable(modA, x >= 0)
        @variable(modA, y <= 5)
        @variable(modA, 2 <= z <= 4)
        @variable(modA, 0 <= r[i=3:6] <= i)
        @NLobjective(modA, Max, ((x + y)/2.0 + 3.0)/3.0 + z + r[3])
        @constraint(modA, cons1, x+y >= 2)
        @constraint(modA, cons2, sum{r[i],i=3:5} <= (2 - x)/2.0)
        cons3 = @NLconstraint(modA, 7.0*y <= z + r[6]/1.9)

        # Solution
        @test solve(modA) == :Optimal
        @test isapprox(getobjectivevalue(modA), 5.8446115, atol=1e-6)
        @test isapprox(getvalue(x), 0.9774436, atol=1e-6)
        @test isapprox(getvalue(y), 1.0225563, atol=1e-6)
        @test isapprox(getvalue(z), 4.0, atol=1e-6)
        @test isapprox(getvalue(r)[3], 0.5112781, atol=1e-6)
        @test isapprox(getvalue(r)[4], 0.0, atol=1e-6)
        @test isapprox(getvalue(r)[5], 0.0, atol=1e-6)
        # knitro has issue with atol=1e-6 :(
        @test isapprox(getvalue(r)[6], 6.0, atol=1e-5)

        # Reduced costs
        @test isapprox(getdual(x), 0.0, atol=1e-6)
        @test isapprox(getdual(y), 0.0, atol=1e-6)
        @test isapprox(getdual(z), 1.0714286, atol=1e-6)
        @test isapprox(getdual(r)[3], 0.0, atol=1e-6)
        @test isapprox(getdual(r)[4], -1.0, atol=1e-6)
        @test isapprox(getdual(r)[5], -1.0, atol=1e-6)
        @test isapprox(getdual(r)[6], 0.03759398, atol=1e-6)

        # Row duals
        @test isapprox(getdual(cons1), -0.333333, atol=1e-6)
        @test isapprox(getdual(cons2), 1.0, atol=1e-6)
        @test isapprox(getdual(cons3), 0.0714286, atol=1e-6)
    end

    @testset "Changing objectives with $nlp_solver" for nlp_solver in nlp_solvers
        m = Model(solver=nlp_solver)
        @variable(m, x >= 0)
        @variable(m, y >= 0)
        @objective(m, Max, x+y)
        @NLconstraint(m, x+2y <= 1)
        @test solve(m) == :Optimal
        @test isapprox(getvalue(x), 1.0, atol=1e-4)
        @test isapprox(getvalue(y), 0.0, atol=1e-4)
        @test isapprox(getobjectivevalue(m), 1.0, atol=1e-4)

        @objective(m, Max, 2x+y)
        @test solve(m) == :Optimal
        @test isapprox(getvalue(x), 1.0, atol=1e-4)
        @test isapprox(getvalue(y), 0.0, atol=1e-4)
        @test isapprox(getobjectivevalue(m), 2.0, atol=1e-4)
    end

    @testset "Hessian chunking code with $nlp_solver" for nlp_solver in nlp_solvers
        m = Model(solver=nlp_solver)
        @variable(m, x[1:18] >= 1, start = 1.2)
        @NLobjective(m, Min, prod(x[i] for i=1:18))
        @test solve(m) == :Optimal
        @test isapprox(getvalue(x), ones(18), atol=1e-4)
    end

    #############################################################################
    # Test that output is produced in correct MPB form
    MathProgBase.NonlinearModel(s::DummyNLPSolver) = DummyNLPModel()
    function MathProgBase.loadproblem!(m::DummyNLPModel, numVar, numConstr, x_l, x_u, g_lb, g_ub, sense, d::MathProgBase.AbstractNLPEvaluator)
        MathProgBase.initialize(d, [:ExprGraph])
        objexpr = MathProgBase.obj_expr(d)
        @testset "NL MPB interface ($objexpr)" begin
            @test objexpr == :(x[1]^x[2]) || objexpr == :(-1.0*x[1]+1.0*x[2])
            @test MathProgBase.isconstrlinear(d,1)
            @test MathProgBase.isconstrlinear(d,3)
            @test MathProgBase.constr_expr(d,1) == :(2.0*x[1] + 1.0*x[2] <= 1.0)
            @test MathProgBase.constr_expr(d,2) == :(2.0*x[1] + 1.0*x[2] <= -0.0)
            @test MathProgBase.constr_expr(d,3) == :(-5.0 <= 2.0*x[1] + 1.0*x[2] <= 5.0)
            if numConstr > 3
                @test MathProgBase.constr_expr(d,4) == :(2.0*x[1]*x[1] + 1.0*x[2] + -2.0 >= 0)
                @test MathProgBase.constr_expr(d,5) == :(sin(x[1]) * cos(x[2]) - 5 == 0.0)
                @test MathProgBase.constr_expr(d,6) == :(1.0*x[1]^2 - 1.0 == 0.0)
                @test MathProgBase.constr_expr(d,7) == :(2.0*x[1]^2 - 2.0 == 0.0)
                @test MathProgBase.constr_expr(d,8) == :(-0.5 <= sin(x[1]) <= 0.5)
                @test MathProgBase.constr_expr(d,9) == :(ψ(x[1]) + Ψ(x[1],x[2]) - 3.0 <= 0.0)
            end
        end
    end
    MathProgBase.setwarmstart!(m::DummyNLPModel,x) = nothing
    MathProgBase.optimize!(m::DummyNLPModel) = nothing
    MathProgBase.status(m::DummyNLPModel) = :Optimal
    MathProgBase.getobjval(m::DummyNLPModel) = NaN
    MathProgBase.getsolution(m::DummyNLPModel) = [1.0,1.0]
    MathProgBase.setvartype!(m::DummyNLPModel,vartype) = @test any(vartype .== :Fixed) == false

    function test_nl_mpb()
        m = Model(solver=DummyNLPSolver())
        @variable(m, x == 1)
        @variable(m, y, Bin)
        @objective(m, Min, -x+y)
        @constraint(m, 2x+y <= 1)
        @constraint(m, 2x+y <= 0)
        @constraint(m, -5 <= 2x+y <= 5)
        #solve(m) # FIXME maybe?
        lb,ub = JuMP.constraintbounds(m)
        @test lb == [-Inf,-Inf,-5.0]
        @test ub == [1.0,-0.0,5.0]

        @constraint(m, 2x^2+y >= 2)
        ψ(x) = 1
        Ψ(x,y) = 2
        JuMP.register(m, :ψ, 1, ψ, autodiff=true)
        JuMP.register(m, :Ψ, 2, Ψ, autodiff=true)

        @NLconstraint(m, sin(x)*cos(y) == 5)
        @NLconstraint(m, nlconstr[i=1:2], i*x^2 == i)
        @NLconstraint(m, -0.5 <= sin(x) <= 0.5)
        @NLconstraint(m, ψ(x) + Ψ(x,y) <= 3)
        solve(m)

        @NLobjective(m, Min, x^y)
        solve(m)

        lb,ub = JuMP.constraintbounds(m)
        @test lb == [-Inf,-Inf,-5.0,0.0,0.0,0.0,0.0,-0.5,-Inf]
        @test ub == [1.0,-0.0,5.0,Inf,0.0,0.0,0.0,0.5,0.0]
    end
    test_nl_mpb()

    @testset "Simplified expression graphs" begin
        m = Model(simplify_nonlinear_expressions=true)
        # this is expert behavior, expression simplification is experimental
        @variable(m, x == 2)
        @variable(m, y)
        @NLobjective(m, Min, x^2 + y^2)
        @NLexpression(m, ex, exp(x))
        @NLconstraint(m, ex - y == 0)
        @NLconstraint(m, ex + 1 == 0)

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:ExprGraph])
        @test MathProgBase.obj_expr(d) == :(4.0 + x[2] ^ 2.0)
        @test MathProgBase.constr_expr(d,1) == :(($(exp(2)) - x[2]) - 0.0 == 0.0)
        @test MathProgBase.constr_expr(d,2) == :($(exp(2) + 1) == 0.0)
    end

    @testset "Expression graph for linear problem" begin
        m = Model()
        @variable(m, x)
        @constraint(m, 0 <= x <= 1)
        @objective(m, Max, x)
        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:ExprGraph])
        @test MathProgBase.obj_expr(d) == :(+(1.0 * x[1]))
    end

    @testset "Expression graph for ifelse" begin
        m = Model()
        @variable(m, x, start = 2)
        @NLobjective(m, Min, ifelse( x <= 1, x^2, x) )
        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:ExprGraph])
        @test MathProgBase.obj_expr(d) == :(ifelse( x[1] <= 1, x[1]^2, x[1]))
    end

    @testset "Expression graphs for corner cases" begin
        m = Model()
        @variable(m, x, start = 2)
        @constraint(m, 0 <= 1)
        @NLconstraint(m, x <= sum(0 for i in []) + prod(1 for i in []))
        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:ExprGraph])
        @test MathProgBase.constr_expr(d,1) == :(0 <= 1.0)
        @test MathProgBase.constr_expr(d,2) == :(x[1] - (0 + 1) <= 0.0)
    end

    @testset "Hessians through MPB" begin
        # Issue 435
        m = Model()
        @variable(m, a, start = 1)
        @variable(m, b, start = 2)
        @variable(m, c, start = 3)

        @NLexpression(m, foo, a * b + c^2)

        @NLobjective(m, Min, foo)
        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Hess])
        I,J = MathProgBase.hesslag_structure(d)
        V = zeros(length(I))
        MathProgBase.eval_hesslag(d, V, m.colVal, 1.0, Float64[])
        hess_raw = sparse(I,J,V)
        # Convert from lower triangular
        hess_sparse = hess_raw + hess_raw' - sparse(Diagonal(diag(hess_raw)))
        @test isapprox(hess_sparse, [0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 2.0])

        # make sure we don't get NaNs in this case
        @NLobjective(m, Min, a * b + 3*c^2)
        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:Hess])
        setvalue(c, -1.0)
        V = zeros(length(I))
        MathProgBase.eval_hesslag(d, V, m.colVal, 1.0, Float64[])
        hess_raw = sparse(I,J,V)
        hess_sparse = hess_raw + hess_raw' - sparse(Diagonal(diag(hess_raw)))
        @test isapprox(hess_sparse, [0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 6.0])

        # Initialize again
        MathProgBase.initialize(d, [:Hess])
        V = zeros(length(I))
        MathProgBase.eval_hesslag(d, V, m.colVal, 1.0, Float64[])
        hess_raw = sparse(I,J,V)
        hess_sparse = hess_raw + hess_raw' - sparse(Diagonal(diag(hess_raw)))
        @test isapprox(hess_sparse, [0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 6.0])
    end

    @testset "Hess-vec through MPB" begin
        m = Model()
        @variable(m, a, start = 1)
        @variable(m, b, start = 2)
        @variable(m, c, start = 3)

        @NLobjective(m, Min, a*b + c^2)
        @constraint(m, c*b <= 1)
        @NLconstraint(m, a^2/2 <= 1)
        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:HessVec])
        h = ones(3) # test that input values are overwritten
        v = [2.4,3.5,1.2]
        MathProgBase.eval_hesslag_prod(d, h, m.colVal, v, 1.0, [2.0,3.0])
        correct = [3.0 1.0 0.0; 1.0 0.0 2.0; 0.0 2.0 2.0]*v
        @test isapprox(h, correct)
    end

    @testset "Hess-vec through MPB with subexpressions" begin
        m = Model()
        @variable(m, a, start = 1)
        @variable(m, b, start = 2)
        @variable(m, c, start = 3)

        @NLexpression(m, ab, a*b)
        @NLobjective(m, Min, ab + c^2)
        @constraint(m, c*b <= 1)
        @NLconstraint(m, a^2/2 <= 1)
        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:HessVec])
        h = ones(3) # test that input values are overwritten
        v = [2.4,3.5,1.2]
        MathProgBase.eval_hesslag_prod(d, h, m.colVal, v, 1.0, [2.0,3.0])
        correct = [3.0 1.0 0.0; 1.0 0.0 2.0; 0.0 2.0 2.0]*v
        @test isapprox(h, correct)
    end

    @testset "NaN corner case (#695)" begin

        m = Model()
        x0 = 0.0
        y0 = 0.0

        @variable(m, x >= -1, start = 1.0)
        @variable(m, y, start = 2.0)

        @NLobjective(m, Min, (x - x0) /(sqrt(y0) + sqrt(y)))

        d = JuMP.NLPEvaluator(m)
        MathProgBase.initialize(d, [:HessVec])
        h = ones(2)
        v = [2.4,3.5]
        MathProgBase.eval_hesslag_prod(d, h, m.colVal, v, 1.0, Float64[])
        correct = [0.0 -1/(2*2^(3/2)); -1/(2*2^(3/2)) 3/(4*2^(5/2))]*v
        @test isapprox(h, correct)
    end

    mysquare(x) = x^2
    function myf(x,y)
        return (x-1)^2+(y-2)^2
    end

    println("User-defined functions not yet supported")
    if length(convex_nlp_solvers) > 0 && false
        @testset "User-defined functions" begin
            m = Model(solver=convex_nlp_solvers[1])

            JuMP.register(m, :myf, 2, myf, autodiff=true)
            JuMP.register(m, :myf_2, 2, myf, (g,x,y) -> (g[1] = 2(x-1); g[2] = 2(y-2)))
            JuMP.register(m, :mysquare, 1, mysquare, autodiff=true)
            JuMP.register(m, :mysquare_2, 1, mysquare, x-> 2x, autodiff=true)
            JuMP.register(m, :mysquare_3, 1, mysquare, x-> 2x, x -> 2.0)

            @variable(m, x[1:2] >= 0.5)
            @NLobjective(m, Min, myf(x[1],mysquare(x[2])))

            d = JuMP.NLPEvaluator(m)
            MathProgBase.initialize(d, [:Grad])
            gradout = zeros(2)
            xval = [1,sqrt(2.0)]
            @test isapprox(MathProgBase.eval_f(d, xval), 0.0, atol=1e-10)
            MathProgBase.eval_grad_f(d, gradout, xval)
            @test isapprox(gradout, [0.0,0.0], atol=1e-10)

            @test solve(m) == :Optimal

            @test isapprox(getvalue(x), xval)

            @NLobjective(m, Min, myf_2(x[1],mysquare_2(x[2])))

            d = JuMP.NLPEvaluator(m)
            MathProgBase.initialize(d, [:Grad])
            gradout = zeros(2)
            xval = [1,sqrt(2.0)]
            @test isapprox(MathProgBase.eval_f(d, xval), 0.0, atol=1e-10)
            MathProgBase.eval_grad_f(d, gradout, xval)
            @test isapprox(gradout, [0.0,0.0], atol=1e-10)
            setvalue(x[1],0.5)
            setvalue(x[2],0.5)
            @test solve(m) == :Optimal

            @test isapprox(getvalue(x), xval)

            # Test just univariate functions because hessians are disabled
            # if any multivariate functions are present.
            @NLobjective(m, Min, mysquare(x[1]-1) + mysquare_2(x[2]-2) + mysquare_3(x[1]))
            @test solve(m) == :Optimal
            @test isapprox(getvalue(x), [0.5,2.0], atol=1e-4)

            # Test #927
            m = Model(solver=convex_nlp_solvers[1])
            JuMP.register(m, :myf, 2, myf, autodiff=true)
            @variable(m, x)
            @NLobjective(m, Min, myf(x,x))
            @test solve(m) == :Optimal
            @test isapprox(getvalue(x), 1.5, atol=1e-4)

        end

        @testset "Anonymous nonlinear expression" begin
            m = Model(solver=convex_nlp_solvers[1])
            @variable(m, -1 <= x <= 1)
            obj = @NLexpression(m, x^4 + 1)
            @NLobjective(m, Min, obj)
            @test solve(m) == :Optimal
            @test isapprox(getvalue(x), 0, atol=1e-4)
        end
    end
end


text:
Tables.istable(::Type{<:File}) = true
Tables.columnaccess(::Type{<:File}) = true
Tables.schema(f::File)  = Tables.Schema(getnames(f), _eltype.(gettypes(f)))
Tables.columns(f::File) = f
Base.propertynames(f::File) = getnames(f)

struct Column{T, P} <: AbstractVector{T}
    file::File
    col::Int
end

_eltype(::Type{T}) where {T} = T
_eltype(::Type{PooledString}) = String
_eltype(::Type{Union{PooledString, Missing}}) = Union{String, Missing}

function Column(f::File, i::Int)
    @inbounds T = gettypes(f)[i]
    return Column{_eltype(T), T}(f, i)
end

Base.size(c::Column) = (Int(getrows(c.file)),)
Base.IndexStyle(::Type{<:Column}) = Base.IndexLinear()
metaind(x) = 2 * x - 1
valind(x) = 2 * x

Base.setindex!(c::Column, x, i::Int) = throw(ArgumentError("CSV.Column is read-only; to get a mutable vector, do `copy(col)` or to make all columns mutable do `CSV.read(file; copycols=true)` or `CSV.File(file) |> DataFrame`"))

function Base.copy(c::Column{T}) where {T}
    len = length(c)
    A = Vector{T}(undef, len)
    @simd for i = 1:len
        @inbounds A[i] = c[i]
    end
    return A
end

function Base.copy(c::Column{T, T}) where {T <: Union{String, Union{String, Missing}}}
    len = length(c)
    A = StringVector{T}(undef, len)
    @simd for i = 1:len
        @inbounds A[i] = c[i]
    end
    return A
end

function Base.copy(c::Column{T, S}) where {T <: Union{String, Union{String, Missing}}, S <: Union{PooledString, Union{PooledString, Missing}}}
    len = length(c)
    catg = getcategorical(c.file)
    tape = gettape(c.file, c.col)
    if S === PooledString
        refs = Dict{String, UInt32}()
        foreach(x->setindex!(refs, UInt32(x[1]), x[2]), enumerate(getrefs(c.file, c.col)))
        values = Vector{UInt32}(undef, len)
        @simd for i = 1:len
            @inbounds values[i] = ref(tape[valind(i)])
        end
    else
        if catg
            refs = Dict{String, UInt32}()
            foreach(x->setindex!(refs, UInt32(x[1]), x[2]), enumerate(getrefs(c.file, c.col)))
            missingref = UInt32(0)
            values = Vector{UInt32}(undef, len)
        else
            refs = Dict{Union{String, Missing}, UInt32}()
            foreach(x->setindex!(refs, UInt32(x[1]), x[2]), enumerate(getrefs(c.file, c.col)))
            missingref = UInt32(length(refs) + 1)
            refs[missing] = missingref
            values = Vector{UInt32}(undef, len)
        end
        @simd for i = 1:len
            @inbounds offlen = tape[metaind(i)]
            @inbounds values[i] = ifelse(missingvalue(offlen), missingref, ref(tape[valind(i)]))
        end
    end
    if catg
        pool = CategoricalPool(refs)
        levels!(pool, sort(levels(pool)))
        A = CategoricalArray{T, 1}(values, pool)
    else
        A = PooledArray(PooledArrays.RefArray(values), refs)
    end
    return A
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{Missing}, row::Int)
    @boundscheck checkbounds(c, row)
    return missing
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{T}, row::Int) where {T}
    @boundscheck checkbounds(c, row)
    @inbounds x = reinterp_func(T)(gettape(c.file, c.col)[valind(row)])
    return x
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{Union{T, Missing}}, row::Int) where {T}
    @boundscheck checkbounds(c, row)
    @inbounds offlen = gettape(c.file, c.col)[metaind(row)]
    @inbounds x = ifelse(missingvalue(offlen), missing, reinterp_func(T)(gettape(c.file, c.col)[valind(row)]))
    return x
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{Float64}, row::Int)
    @boundscheck checkbounds(c, row)
    @inbounds offlen = gettape(c.file, c.col)[metaind(row)]
    @inbounds v = gettape(c.file, c.col)[valind(row)]
    @inbounds x = ifelse(intvalue(offlen), Float64(int64(v)), float64(v))
    return x
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{Union{Float64, Missing}}, row::Int)
    @boundscheck checkbounds(c, row)
    @inbounds offlen = gettape(c.file, c.col)[metaind(row)]
    @inbounds v = gettape(c.file, c.col)[valind(row)]
    @inbounds x = ifelse(missingvalue(offlen), missing, ifelse(intvalue(offlen), Float64(int64(v)), float64(v)))
    return x
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{String, PooledString}, row::Int)
    @boundscheck checkbounds(c, row)
    @inbounds x = getrefs(c.file, c.col)[gettape(c.file, c.col)[valind(row)]]
    return x
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{Union{String, Missing}, Union{PooledString, Missing}}, row::Int)
    @boundscheck checkbounds(c, row)
    @inbounds offlen = gettape(c.file, c.col)[metaind(row)]
    if missingvalue(offlen)
        return missing
    else
        @inbounds x = getrefs(c.file, c.col)[gettape(c.file, c.col)[valind(row)]]
        return x
    end
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{String}, row::Int)
    @boundscheck checkbounds(c, row)
    @inbounds offlen = gettape(c.file, c.col)[metaind(row)]
    s = PointerString(pointer(getbuf(c.file), getpos(offlen)), getlen(offlen))
    return escapedvalue(offlen) ? unescape(s, gete(c.file)) : String(s)
end

@inline Base.@propagate_inbounds function Base.getindex(c::Column{Union{String, Missing}}, row::Int)
    @boundscheck checkbounds(c, row)
    @inbounds offlen = gettape(c.file, c.col)[metaind(row)]
    if missingvalue(offlen)
        return missing
    else
        s = PointerString(pointer(getbuf(c.file), getpos(offlen)), getlen(offlen))
        return escapedvalue(offlen) ? unescape(s, gete(c.file)) : String(s)
    end
end

function Base.getproperty(f::File, col::Symbol)
    i = findfirst(==(col), getnames(f))
    i === nothing && return getfield(f, col)
    return Column(f, i)
end


text:
# Map onto SE(3) embedded in R12
if !(@isdefined LinkFrameTransform)
    @computed struct LinkFrameTransform{Rn,R12,S,JC} <: TaskMap{Rn,R12,S}
        state_cache::StateCache{S,JC}
        link_cache::fulltype(LinkFrameCache{dim(Rn),S})
    end
end
struct LinkFrameTransformT{Rn,R12,S,JC} <: TaskMapT{Rn,R12,S}
    base_map::LinkFrameTransform{Rn,R12,S,JC}
end
LinkFrameTransform{ℝ{n},R12}(state_cache::StateCache{S,JC},
        link_cache::LinkFrameCache{n,S}) where {n,R12,S,JC} =
    LinkFrameTransform{ℝ{n},R12,S,JC}(state_cache, link_cache)

task_map_emb(::EmbRep, ::EmbRep, xme, task_map::LinkFrameTransform{ℝ{n},R12}) where n =
    link_frame_transform(xme, task_map.state_cache, task_map.link_cache.link)

function task_jacobian_emb(xme, task_map::LinkFrameTransform{ℝ{n},R12}, arg...) where n
    state_cache, link_cache = task_map.state_cache, task_map.link_cache
    Jfc = link_cache.transform_jacobian
    update_link_transform_jacobian!(Jfc, xme, state_cache, link_cache.link)
    Jf = Jfc.data[1]
end

task_jacobian_chart(xm, task_map::LinkFrameTransform{ℝ{n},R12},
    CM::Chart{1,ℝ{n}}, CN::Chart{1,R12}) where n = task_jacobian_emb(xm, task_map)

function task_jacobian_emb_dot(xme, vme, task_map::LinkFrameTransform{ℝ{n},R12}, arg...) where n
    state_cache, link_cache = task_map.state_cache, task_map.link_cache
    Jfc_dot = link_cache.transform_jacobian_dot
    update_link_transform_jacobian_dot!(Jfc_dot, xme, vme, state_cache, link_cache.link)
    Jf_dot = Jfc_dot.data[1]
end

task_jacobian_chart_dot(xm, vm, task_map::LinkFrameTransform{ℝ{n},R12},
    CM::Chart{1,ℝ{n}}, CN::Chart{1,R12}) where n = task_jacobian_emb_dot(xm, vm, task_map)


text:
function blocksize end
function blockgetindex end

macro hierarchical(HierarchicalType, Types...)
    blocks = :(begin end)
    push!(blocks.args, :($(Symbol("$(HierarchicalType)blocks"))::Matrix{$HierarchicalType{T}}))
    for i in 1:length(Types)
        if typeof(Types[i]) == Symbol
            push!(blocks.args, :($(Symbol("$(Types[i])blocks"))::Matrix{$(Types[i]){T}}))
        elseif typeof(Types[i]) == Expr
            push!(blocks.args, :($(Symbol("$(Types[i].args[1])blocks"))::Matrix{$(Types[i].args[1]){T,$(Types[i].args[2]){T}}}))
        else
            error("Invalid hierarchical symbol.")
        end
    end
    AbstractHierarchicalType = Meta.parse("Abstract"*string(HierarchicalType))
    Factorization = Meta.parse(string(HierarchicalType)*"Factorization")
    return esc(quote
        import Base: +, -, *, /, \, ==
        import Base: size, getindex, setindex!
        import LinearAlgebra: Factorization

        # import HierarchicalMatrices:

        AbstractSuperType = promote_type([typeof($Types[i]) == Symbol ? eval($Types[i]) : eval($Types[i].args[1]) for i in 1:length($Types)]...)

        export $AbstractHierarchicalType, $HierarchicalType, $Factorization

        abstract type $AbstractHierarchicalType{T} <: AbstractSuperType{T} end

        HierarchicalMatrices.blocksize(H::$AbstractHierarchicalType) = size(H.assigned)

        function size(H::$AbstractHierarchicalType)
            M, N = HierarchicalMatrices.blocksize(H)

            p = 0
            for m = 1:M
                p += HierarchicalMatrices.blocksize(H, m, N, 1)
            end

            q = 0
            for n = 1:N
                q += HierarchicalMatrices.blocksize(H, 1, n, 2)
            end

            p, q
        end

        struct $HierarchicalType{T} <: $AbstractHierarchicalType{T}
            $blocks
            assigned::Matrix{Int}
        end

        @generated function $HierarchicalType(::Type{T}, M::Int, N::Int) where T
            L = length(fieldnames($HierarchicalType))
            HM = $HierarchicalType
            str = "$HM(Matrix{$HM{T}}(undef, M, N), "
            for l in 2:L-1
                S = $Types[l-1]
                if typeof(S) == Symbol
                    str *= "Matrix{$S{T}}(undef, M, N), "
                else
                    str *= "Matrix{$(S.args[1]){T, $(S.args[2]){T}}}(undef, M, N), "
                end
            end
            str *= "zeros(Int, M, N))"
            return Meta.parse(str)
        end
        $HierarchicalType(M::Int, N::Int) = $HierarchicalType(Float64, M, N)

        struct $Factorization{T} <: Factorization{T}
            $HierarchicalType::$HierarchicalType{T}
            factors::Matrix{Matrix{T}}
        end

        HierarchicalMatrices.blocksize(H::$HierarchicalType, m::Int, n::Int) = HierarchicalMatrices.blocksize(H, m, n, 1), HierarchicalMatrices.blocksize(H, m, n, 2)

        @generated function HierarchicalMatrices.blocksize(H::$HierarchicalType, m::Int, n::Int, k::Int)
            L = length(fieldnames(H))-1
            T = fieldname(H, 1)
            str = "
            begin
                Hmn = H.assigned[m,n]
                if Hmn == 1
                    return size(getindex(H.$T, m, n), k)"
            for l in 2:L
                T = fieldname(H, l)
                str *= "
                elseif Hmn == $l
                    return size(getindex(H.$T, m, n), k)"
            end
            str *= "
                end
                return 0
            end"
            return Meta.parse(str)
        end

        @generated function HierarchicalMatrices.blockgetindex(H::$HierarchicalType{S}, m::Int, n::Int, i::Int, j::Int) where S
            L = length(fieldnames(H))-1
            T = fieldname(H, 1)
            str = "
            begin
                Hmn = H.assigned[m,n]
                if Hmn == 1
                    return getindex(getindex(H.$T, m, n), i, j)"
            for l in 2:L
                T = fieldname(H, l)
                str *= "
                elseif Hmn == $l
                    return getindex(getindex(H.$T, m, n), i, j)"
            end
            str *= "
                end
                return zero(S)
            end"
            return Meta.parse(str)
        end

        function getindex(H::$HierarchicalType, i::Int, j::Int)
            p, q = size(H)
            M, N = HierarchicalMatrices.blocksize(H)

            m = 1
            while m ≤ M
                r = HierarchicalMatrices.blocksize(H, m, N, 1)
                if i > r
                    i -= r
                    m += 1
                else
                    break
                end
            end

            n = 1
            while n ≤ N
                s = HierarchicalMatrices.blocksize(H, 1, n, 2)
                if j > s
                    j -= s
                    n += 1
                else
                    break
                end
            end

            HierarchicalMatrices.blockgetindex(H, m, n, i, j)
        end

        @generated function setindex!(H::$HierarchicalType{S}, A::AbstractMatrix{S}, B1::Block, B2::Block) where S
            L = length(fieldnames(H))-1
            T = fieldname(H, 1)
            str = "
            begin
                m, n = B1.K, B2.K
                if typeof(H.$T) == Matrix{typeof(A)}
                    setindex!(H.$T, A, m, n)
                    H.assigned[m,n] = 1
                    return H"
            for l in 2:L
                T = fieldname(H, l)
                str *= "
                elseif typeof(H.$T) == Matrix{typeof(A)}
                    setindex!(H.$T, A, m, n)
                    H.assigned[m,n] = $l
                    return H"
            end
            str *= "
                end
                return H
            end"
            return Meta.parse(str)
        end

        @generated function (+)(G::$HierarchicalType, H::$HierarchicalType)
            L = length(fieldnames(H))-1
            T = fieldname(H, 1)
            str = "
            begin
                F = deepcopy(G)
                M, N = HierarchicalMatrices.blocksize(H)
                for m = 1:M
                    for n = 1:N
                        Gmn = G.assigned[m,n]
                        Hmn = H.assigned[m,n]
                        if Hmn == Gmn == 1
                            F[Block(m), Block(n)] = getindex(G.$T, m, n) + getindex(H.$T, m, n)"
            for l in 2:L
                T = fieldname(H, l)
                str *= "
                        elseif Hmn == Gmn == $l
                            F[Block(m), Block(n)] = getindex(G.$T, m, n) + getindex(H.$T, m, n)"
            end
            str *= "
                        end
                    end
                end
                return F
            end"
            return Meta.parse(str)
        end

        @generated function (-)(G::$HierarchicalType, H::$HierarchicalType)
            L = length(fieldnames(H))-1
            T = fieldname(H, 1)
            str = "
            begin
                F = deepcopy(G)
                M, N = HierarchicalMatrices.blocksize(H)
                for m = 1:M
                    for n = 1:N
                        Gmn = G.assigned[m,n]
                        Hmn = H.assigned[m,n]
                        if Hmn == Gmn == 1
                            F[Block(m), Block(n)] = getindex(G.$T, m, n) - getindex(H.$T, m, n)"
            for l in 2:L
                T = fieldname(H, l)
                str *= "
                        elseif Hmn == Gmn == $l
                            F[Block(m), Block(n)] = getindex(G.$T, m, n) - getindex(H.$T, m, n)"
            end
            str *= "
                        end
                    end
                end
                return F
            end"
            return Meta.parse(str)
        end
    end)
end


text:
"""
Aitchison GLN layer with arbitrary number of units.
"""
struct AGLNLayer{
        MF <: AbstractMatrix{<:Real},
        VF <: AbstractVector{<:Real},
        VI <: AbstractVector{Int32}}
    input_dim::Int
    output_dim::Int
    context_dim::Int
    predictor_dim::Int
    prediction_dim::Int

    learning_rate::Float64

    lsh::LocalitySensitiveHash{MF, VF, VI}

    weights::MF
end

"""
Aitchison geometry closure function.
"""
function closure(p::AbstractArray, simplex_dim::Int)
    return p ./ sum(p, dims=simplex_dim)
end

function closure!(p::AbstractArray, simplex_dim::Int)
    p ./= sum(p, dims=simplex_dim)
end

"""
Aitchison geometry powering function.
"""
function power(p::AbstractArray, w::AbstractArray, simplex_dim::Int)
    @show size(p)
    @show size(w)
    return closure(p .^ w, simplex_dim)
end

"""
Aitchison geometry perturbation function.
"""
function perturb(p::AbstractArray, perturb_dim::Int, simplex_dim::Int)
    return dropdims(closure(prod(p, dims=perturb_dim), simplex_dim), dims=perturb_dim)
end

"""
Construct a single AGLN layer.

  * `input_dim`: number of inputs (i.e. number of units in the prev layer)
  * `output_dim`: number of outputs (i.e. number of units in this layer)
  * `context_dim`: number of context functions (inducing 2^context_dim weight vectors)
  * `predictor_dim`: dimensionality of predictors
  * `prediction_dim`: dimensionality of predictions
  * `learning_rate`: controls step size
"""
function AGLNLayer(
        input_dim::Int, output_dim::Int,
        context_dim::Int, predictor_dim::Int, prediction_dim::Int,
        learning_rate::Float64=1e-2)

    lsh = LocalitySensitiveHash(output_dim, context_dim, predictor_dim)
    weights = zeros(Float32,(input_dim, output_dim*(2^context_dim)))

    return AGLNLayer(
        input_dim,
        output_dim,
        context_dim,
        predictor_dim,
        prediction_dim,
        learning_rate,
        lsh,
        weights)
end

"""
Forward pass for Aitchison GLN layer. Combine input prediction probability
vectors into a hopefully better prediction probability vector.
"""
function forward(
        layer::AGLNLayer, p::AbstractArray,
        z::AbstractMatrix, y::Union{Nothing, AbstractMatrix})

    batch_dim = size(p, 3)

    # [output_dim, batch_dim]
    cs = layer.lsh(z)

    function predict(weights)
        # [input_dim, output_dim, 1, batch_dim]
        weights_cs = weights[:,cs]
        penalty = 1f-5 * sum(weights_cs.^2)

        weights_ = exp.(reshape(
            weights_cs, (layer.input_dim, layer.output_dim, 1, batch_dim)))

        # [input_dim, 1, prediction_dim, batch_dim]
        p_ = reshape(p, (layer.input_dim, 1, layer.prediction_dim, batch_dim))

        # weight and combine inputs p_ by powering then perturbation

        # [output_dim, prediction_dim, batch_dim]
        q = perturb(power(p_, weights, 3), 1, 3)

        return q, penalty
    end

    if y !== nothing
        function loss(weights)
            # [output_dim, prediction_dim, batch_dim]
            q, penalty = predict(weights)

            y_ = reshape(y, (1, layer.prediction_dim, batch_dim))

            neg_ll = -sum(log.(q .* y_))

            return neg_ll + penalty
        end

        dloss_dweights = gradient(loss, layer.weights)[1]
        layer.weights .-= layer.learning_rate .* dloss_dweights
        clamp!(layer.weights, -10f0, 10f0)
    end

    return predict(layer.weights)
end

abstract type WeakClassifier end

"""
Simple logistic regression used as the weak learner input into the AGLN.
"""
struct BasicLogisticRegression{AF <: AbstractArray{<:Real, 3}} <: WeakClassifier
    # weights [predictor_dim, prediction_dim, 1]
    w::AF

    # bias [predictor_dim, prediction_dim, 1]
    b::AF

    # hyperparameters for the N(0, σ) prior on weights and biases
    σ_w::Float32
    σ_b::Float32

    # gradience descent learning rate
    learning_rate::Float32
end

"""
Construct simple logistic regression that makes `predictor_dim` predictions
seperately using each of the the predictor dimensions.
"""
function BasicLogisticRegression(
        predictor_dim::Int, prediction_dim::Int;
        σ_w=5.0, σ_b=5.0, learning_rate=1f-2)
    return BasicLogisticRegression(
        1f-3 * randn(Float32, (predictor_dim, prediction_dim, 1)),
        zeros(Float32, (predictor_dim, prediction_dim, 1)),
        Float32(σ_w), Float32(σ_b), Float32(learning_rate))
end

# x: [predictor_dim, batch_dim]
# y: [prediction_dim, batch_dim]
"""
Logistic regression forward pass. If `y` is not nothing, the w
"""
function forward(
        layer::BasicLogisticRegression,
        x::AbstractMatrix, y::Union{Nothing, AbstractMatrix})

    predictor_dim = size(layer.w, 1)
    prediction_dim = size(layer.w, 2)
    batch_dim = size(x, 2)

    @assert size(x, 1) == predictor_dim
    @assert y === nothing || size(y, 1) == prediction_dim
    @assert y === nothing || size(y, 2) == batch_dim

    function predict(w, b)
        # [predictor_dim, prediction_dim, batch_dim]
        logits = b .+ w .* reshape(x, (predictor_dim, 1, batch_dim))

        p = softmax(logits, dims=2)
        @show extrema(p)

        return p
    end

    if y !== nothing
        function loss(w, b)
            p = predict(w, b)

            y_ = reshape(y, 1, prediction_dim, batch_dim)
            true_label_likelihood = dropdims(sum(y_ .* p, dims=2), dims=2)

            ll = sum(log.(true_label_likelihood))
            lprior = sum((w ./ layer.σ_w).^2) + sum((b ./ layer.σ_b).^2)

            lp = ll + lprior

            @show lp

            return -lp
        end

        dloss_dw, dloss_db = gradient(loss, layer.w, layer.b)

        layer.w .-= layer.learning_rate .* dloss_dw
        layer.b .-= layer.learning_rate .* dloss_db
    end

    return predict(layer.w, layer.b)
end

struct AitchisonGLN{AF <: AbstractArray{<:Real,3}}
    x_norm::Union{Nothing, NormalizationLayer}
    weak_learner::BasicLogisticRegression
    layers::Vector{AGLNLayer}

    # [bias_dim, prediction_dim, 1]
    p_bias::AF
end

"""
Construct an untrained Aitchison gated linear network.

  * `predictor_dim`: dimensionality of predictors (x)
  * `prediction_dim`: dimensionality of predictions (y)
  * `num_layers`: number of layers in the model
  * `layer_width`: width of each layer in the model
  * `context_dim`: number of context functions (inducing 2^context_dim weight vectors)
  * `bias`: values for `bias` units.
"""
function AitchisonGLN(
        predictor_dim::Int, prediction_dim::Int,
        num_layers::Int, layer_width::Int, context_dim::Int;
        standardize::Bool=true, bias::Float32=5.0f0)

    bias_dim = prediction_dim - 1

    layers = AGLNLayer[]

    if num_layers > 0
        push!(
            layers, AGLNLayer(
                predictor_dim + bias_dim, layer_width, context_dim,
                predictor_dim, prediction_dim))
    end

    for i in 1:num_layers-1
        push!(
            layers, AGLNLayer(
                layer_width + bias_dim, layer_width, context_dim,
                predictor_dim, prediction_dim))
    end

    bias = zeros(Float32, (bias_dim, prediction_dim, 1))
    sbp_basis!(bias)

    return AitchisonGLN(
        standardize ? NormalizationLayer(predictor_dim) : nothing,
        BasicLogisticRegression(predictor_dim, prediction_dim),
        layers, bias)
end

"""
Construct a simplical basis using sequential binary partitioning.
"""
function sbp_basis!(basis::AbstractArray)
    i = Ref(1)
    sbp_basis!(basis, i, 1, size(basis, 2))
    map!(exp, basis, basis)
    closure!(basis, 2)
end

function sbp_basis!(basis::AbstractArray, i::Ref{Int}, from::Int, to::Int)
    num_leaves = to - from + 1
    if num_leaves <= 1
        return
    end

    mid = div(from + to, 2)

    left_subtree = from:mid
    right_subtree = mid+1:to

    len = length(from:to)
    left_len = length(left_subtree)
    right_len = length(right_subtree)

    basis[i[], from:mid, 1] .= sqrt(right_len) / sqrt(left_len * len)
    basis[i[], mid+1:to, 1] .= -sqrt(left_len) / sqrt(right_len * len)

    i[] += 1

    sbp_basis!(basis, i, from, mid)
    sbp_basis!(basis, i, mid+1, to)
end

"""
Forward pass for the AGLN, training if y is given, predicting otherwise.
"""
function forward(agln::AitchisonGLN, x::AbstractMatrix, y::Union{Nothing, AbstractMatrix})

    batch_dim = size(x, 2)

    @assert y === nothing || size(x, 2) == size(y, 2)
    @assert y === nothing || size(y, 1) == size(agln.p_bias, 2)

    xz = agln.x_norm === nothing ? x : forward(agln.x_norm, x, y !== nothing)

    p = forward(agln.weak_learner, xz, y)

    for layer in agln.layers
        # concatenate bias inputs
        p = cat(repeat(agln.p_bias, 1, 1, batch_dim), p, dims=1)

        # [output_dim, prediction_dim, batch_dim]
        p = forward(layer, p, xz, y)
    end

    return dropdims(mean(p, dims=1), dims=1)
end

"""
Train the model one step using a the batch of predictors `x` and predictions
`y`.

`x` should have shape [predictors dimension, batch size]
`y` should have shape [predictions dimension, batch size]
"""
function train!(agln::AitchisonGLN, x::AbstractMatrix, y::AbstractMatrix)
    return forward(agln, x, y)
end

"""
Make predictions using a trained model for predictors `x`.

`x` should have shape [predictors dimension, batch size]
"""
function predict(agln::AitchisonGLN, x::AbstractMatrix)
    return forward(agln, x, nothing)
end


text:
@doc raw"""
    LocalLevel(y::Vector{Fl}) where Fl

The local level model is defined by:
```math
\begin{gather*}
    \begin{aligned}
        y_{t} &=  \mu_{t} + \varepsilon_{t} \quad \varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_{\varepsilon})\\
        \mu_{t+1} &= \mu_{t} + \eta_{t} \quad \eta_{t} \sim \mathcal{N}(0, \sigma^2_{\eta})\\
    \end{aligned}
\end{gather*}
```

# Example
```jldoctest
julia> model = LocalLevel(rand(100))
LocalLevel model
```

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. pp. 9
"""
mutable struct LocalLevel <: StateSpaceModel
    hyperparameters::HyperParameters
    system::LinearUnivariateTimeInvariant
    results::Results

    function LocalLevel(y::Vector{Fl}) where Fl

        # Define system matrices
        Z = ones(Fl, 1)
        T = ones(Fl, 1, 1)
        R = ones(Fl, 1, 1)
        d = zero(Fl)
        c = zeros(Fl, 1)
        H = one(Fl)
        Q = ones(Fl, 1, 1)

        system = LinearUnivariateTimeInvariant{Fl}(y, Z, T, R, d, c, H, Q)

        # Define hyperparameters names
        names = ["sigma2_ε", "sigma2_η"]
        hyperparameters = HyperParameters{Fl}(names)

        return new(hyperparameters, system, Results{Fl}())
    end
end

# Obligatory methods
function default_filter(model::LocalLevel)
    Fl = typeof_model_elements(model)
    a1 = zero(Fl)
    P1 = Fl(1e6)
    steadystate_tol = Fl(1e-5)
    return ScalarKalmanFilter(a1, P1, 1, steadystate_tol)
end

function initial_hyperparameters!(model::LocalLevel)
    Fl = typeof_model_elements(model)
    observed_variance = var(model.system.y[findall(!isnan, model.system.y)])
    initial_hyperparameters = Dict{String,Fl}(
        "sigma2_ε" => observed_variance, "sigma2_η" => observed_variance
    )
    set_initial_hyperparameters!(model, initial_hyperparameters)
    return model
end

function constrain_hyperparameters!(model::LocalLevel)
    constrain_variance!(model, "sigma2_ε")
    constrain_variance!(model, "sigma2_η")
    return model
end

function unconstrain_hyperparameters!(model::LocalLevel)
    unconstrain_variance!(model, "sigma2_ε")
    unconstrain_variance!(model, "sigma2_η")
    return model
end

function fill_model_system!(model::LocalLevel)
    model.system.H = get_constrained_value(model, "sigma2_ε")
    model.system.Q[1] = get_constrained_value(model, "sigma2_η")
    return model
end

function reinstantiate(::LocalLevel, y::Vector{Fl}) where Fl
    return LocalLevel(y)
end

has_exogenous(::LocalLevel) = false


text:
#=
for variables X in R^{n, m} and Y in S^n:
maximize    tr(C*X) :
subject to  Y - X*X' in S^n_+
            Y_ij = P_ij for (i, j) in Omega
where Omega is a set of fixed indices and P is a random PSD matrix

the nonlinear constraint Y - X*X' in S^n_+ is equivalent to
the conic constraint (Y, 0.5, X) in MatrixEpiPerSquareCone(),
and also to the larger conic constraint [I X'; X Y] in S^{n + m}_+

simple case of nonlinfear matrix inequality example from Lectures on Modern
Convex Optimization by Aharon Ben-Tal and Arkadi Nemirovski, pg 154
=#

using SparseArrays

struct MatrixQuadraticJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    num_rows::Int
    num_cols::Int
    use_matrixepipersquare::Bool # use matrixepipersquare cone, else PSD cone
end

function build(inst::MatrixQuadraticJuMP{T}) where {T <: Float64}
    (num_rows, num_cols) = (inst.num_rows, inst.num_cols)
    C = randn(num_cols, num_rows)
    P = randn(num_rows, num_rows)
    P = Symmetric(P * P')
    (row_idxs, col_idxs, _) = findnz(tril!(sprand(Bool, num_rows,
        num_rows, inv(sqrt(num_rows)))) + I)

    model = JuMP.Model()
    JuMP.@variable(model, X[1:num_rows, 1:num_cols])
    JuMP.@variable(model, Y[1:num_rows, 1:num_rows], Symmetric)
    JuMP.@objective(model, Max, tr(C * X))
    JuMP.@constraint(model, [(row, col) in zip(row_idxs, col_idxs)],
        Y[row, col] == P[row, col])

    if inst.use_matrixepipersquare
        U_svec = zeros(JuMP.GenericAffExpr{T, JuMP.VariableRef},
            Cones.svec_length(num_rows))
        U_svec = Cones.smat_to_svec!(U_svec, 1.0 * Y, sqrt(2))
        JuMP.@constraint(model, vcat(U_svec, 0.5, vec(X)) in
            Hypatia.MatrixEpiPerSquareCone{T, T}(num_rows, num_cols))
    else
        JuMP.@constraint(model, Symmetric(
            [Matrix(I, num_cols, num_cols) X'; X Y]) in JuMP.PSDCone())
    end

    return model
end


text:
using Base.Test

using Pebblegames, Graphtypes

#manual pebble game moves
p = basic_pebble_game(4, 2, 3)

@test add_edge!(p, 1, 2)
@test add_edge!(p, 2, 3)

@test connected(p.graph, 1, 2)
@test !connected(p.graph, 2, 1)
@test !connected(p.graph, 3, 2)
@test connected(p.graph, 2, 3)

@test p.pebbles == [1,1,2,2]

@test slide_move!(p, 2, 3)

@test p.pebbles == [1,2,1,2]

@test add_edge!(p, 2, 4)

@test p.pebbles == [1,1,1,2]

@test slide_move!(p, 1, 2)

@test p.pebbles == [2,0,1,2]
@test connected(p.graph, 2, 1)
@test connected(p.graph, 2, 4)
@test connected(p.graph, 3, 2)

@test add_edge!(p, 1, 4)
@test slide_move!(p, 2, 1)
@test slide_move!(p, 3, 2)
@test add_edge!(p, 3, 4)

@test p.pebbles == [0,0,1,2]
@test connected(p.graph, 1, 2)
@test connected(p.graph, 2, 3)
@test connected(p.graph, 3, 4)
@test connected(p.graph, 2, 4)

@test total_edges(p.graph) == 5

#fetching a pebble automatically
game = basic_pebble_game(5, 2, 3)

game.pebbles = [0, 2, 0, 1, 2]

connect!(game.graph, 1, 3)
connect!(game.graph, 3, 2)
connect!(game.graph, 3, 4)
connect!(game.graph, 1, 5)

go_get_pebbles!(game, 1, [1,2])
@test connected(game.graph, 3,1)
@test connected(game.graph, 4,3)
@test connected(game.graph, 5,1)
@test !connected(game.graph, 1,3)
@test !connected(game.graph, 3,4)
@test !connected(game.graph, 1,5)
@test game.pebbles == [2, 2, 0, 0, 1]

#fetching a pebble automatically
game = basic_pebble_game(3, 2, 3)

game.pebbles = [1, 2, 2]
connect!(game.graph, 1, 2)

go_get_pebbles!(game, 1, [1,3])

@test game.pebbles == [2, 1, 2]

#pebble game solver
m = simple_graph(3)
connect!(m, 1, 2)
connect!(m, 3, 2)
connect!(m, 3, 1)

pebbles_left, rejected_edges = solve_basic_pebble_game(m, 2, 3)

@test pebbles_left == 3 && rejected_edges == 0
@test connected(m, 1, 2) || connected(m, 2, 1)
@test connected(m, 3, 2) || connected(m, 2, 3)
@test connected(m, 1, 3) || connected(m, 3, 1)

#pebble game solver (problem as seen on http://linkage.cs.umass.edu/pg/pg.html)
g = structure_graph(6)
connect!(g, 1, 2)
connect!(g, 1, 3)
connect!(g, 2, 4)
connect!(g, 3, 4)
connect!(g, 3, 5)
connect!(g, 3, 6)
connect!(g, 4, 5)
connect!(g, 4, 6)
connect!(g, 5, 6)

pebbles_left, rejected_edges = solve_basic_pebble_game(g, 2, 3)

@test pebbles_left == 4
@test rejected_edges == 1

#pebble game solver (problem Fig. 3.a. in Lee and Streinu 2008)
g = simple_graph(6)
connect!(g, 1, 2)
connect!(g, 2, 4)
connect!(g, 2, 4)
connect!(g, 1, 3)
connect!(g, 1, 3)
connect!(g, 3, 4)
connect!(g, 1, 4)
connect!(g, 3, 4)
connect!(g, 3, 4)
connect!(g, 1, 5)
connect!(g, 1, 5)
connect!(g, 3, 5)
connect!(g, 3, 6)
connect!(g, 3, 6)
connect!(g, 3, 6)

pebbles_left, rejected_edges = solve_basic_pebble_game(g, 3, 3)

@test pebbles_left == 3
@test rejected_edges == 0

#pebble game solver (problem Fig. 3.b. in Lee and Streinu 2008)
g = simple_graph(6)
connect!(g, 1, 2)
connect!(g, 2, 4)
connect!(g, 2, 4)
connect!(g, 1, 3)
connect!(g, 1, 3)
connect!(g, 3, 4)
connect!(g, 1, 4)
connect!(g, 3, 4)
connect!(g, 3, 4)
connect!(g, 1, 5)
connect!(g, 1, 5)
connect!(g, 3, 6)
connect!(g, 3, 6)
connect!(g, 3, 6)
connect!(g, 2, 3)

pebbles_left, rejected_edges = solve_basic_pebble_game(g, 3, 3)

@test pebbles_left == 4
@test rejected_edges == 1

max_stars_repo_name: AaltoScienceInst/pebblegame

max_issues_repo_name: AaltoScienceInst/pebblegame

max_forks_repo_name: AaltoScienceInst/pebblegame


text:
using MKL

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

using OrdinaryDiffEq, LinearAlgebra, Plots, DiffEqOperators, SparseArrays, Distributions, CUDA
using Symbolics, SparseDiffTools
using PreallocationTools
using PatternFormation
using LinearSolve
using IncompleteLU
using BenchmarkTools
using LazyArrays
using FiniteDiff

BLAS.set_num_threads(16)
CUDA.allowscalar(false)

#=============================================================
Preconditioners
==============================================================#
function incompletelu(W,du,u,p,t,newW,Plprev,Prprev,solverdata)
  if newW === nothing || newW
    Pl = ilu(convert(AbstractMatrix,W), τ = 150.0f0)
  else
    Pl = Plprev
  end
  Pl,nothing
end
Base.eltype(::IncompleteLU.ILUFactorization{Tv,Ti}) where {Tv,Ti} = Tv
#==============================================================
Grid parameters
===============================================================#
dx = 1.0f0
dy = 1.0f0
N = 300
x = range(0, N*dx, length=N)
y = range(0, N*dy, length=N)
tspan = (0.0, Float32(1e3))
#==============================================================
Problem parameters
===============================================================#
println("Setting up problem parameters")
ū = 0.5f0
v̄ = 0.4f0
D = 0.5f0
#D = cu(Float32.(0.1 .+ (0.5-0.1)*exp.(-1/5000*((x.-N*dx/2).^2 .+ (y.-N*dx/2)'.^2))))
#D,_,_,_,_ = init_cond(:ErMnO3_ErOnly, N, M = 9, α=1/80, BC=:Periodic, c01=0.3, dx = dx, f=2/3)
#D = cu(D)
β = -1.0f0 #FIXED
r₁ = 50.0f0
r₂ = 10.0f0
α = -(v̄ - ū*v̄*r₂)./(ū .- ū*v̄^2*r₁)
γ = -α
δ = 1.0f0 #1/(5*pi)^2*(α + β*D)/(2D)

##########################################################
A = [10 1; 1 10]
#γ_x = kron(A, ones(N÷2,N÷2))
checker = kron(ones(N÷4,N÷4), A)
gauss = 1.0 .+ (10.0-1.0)*exp.(-1/5000*((x[1:N÷2].-N*dx/4).^2 .+ (y[1:N÷2].-N*dx/4)'.^2))
γ_x = [checker 10*ones(N÷2,N÷2); ones(N÷2,N÷2) gauss]
γ_y = [reverse(checker, dims=2) ones(N÷2,N÷2); 10*ones(N÷2,N÷2) (10.0+1.0).-gauss]
#-----------------------------
f = 1e-3
w = 1
γ_x[N÷2-w+1:N÷2+w,:] .= f
γ_x[:, N÷2-w+1:N÷2+w] .= f
γ_y[N÷2-w+1:N÷2+w,:] .= f
γ_y[:, N÷2-w+1:N÷2+w] .= f
γ_x[1:w, :] .= f
γ_x[end-w+1:end, :] .= f
γ_x[:,1:w] .= f
γ_x[:, end-w+1:end] .= f
γ_y[1:w, :] .= f
γ_y[end-w+1:end, :] .= f
γ_y[:,1:w] .= f
γ_y[:, end-w+1:end] .= f
#----------------
#B = 5 .+ (1-5)*exp.(-1/2500*((x[1:end÷2].-N*dx/4).^2 .+ (y[1:end÷2].-N*dx/4)'.^2))
#A = 5 .+ (9-5)*exp.(-1/2500*((x[1:end÷2].-N*dx/4).^2 .+ (y[1:end÷2].-N*dx/4)'.^2))
#γ_x = cu([A B; B A])
#γ_y = cu(reverse(γ_x, dims=2))
#γ_x = 10+1e-3 .+ 5*sin.(2*pi*5*x) .+ 5*cos.(2*pi*5*y')
#γ_y = 20 .- γ_x
#----------------
γ_x = 2.0f0
γ_y = 1.0f0
#-----------------
A = [10 1; 1 10]
checker = kron(A, ones(N÷20,N÷20))
γ_x = kron(ones(N÷30, N÷30), checker)
γ_y = reverse(γ_x, dims=2)
#------------------------------------
γ_x = -9(x/2 .+ y'/2)/N.+10
γ_y = reverse(γ_x')
#------------------------------------
γ_x = 5*sin.(10*pi*x.+10*pi*y') .+ (5+1e-3)
γ_y = 1/2*sin.(10*pi*x.+10*pi*y' .+ pi) .+ (1/2+1e-3)
##########################################################################################
γᵤ = 1.0f0
γᵥ = 1.0f0
#=============================================================
Initial Conditions
=============================================================#
println("Initial Conditions")
r01 = Float32.(ū .+ rand(Uniform(-ū/100, +ū/100), N, N))
r02 = Float32.(v̄ .+ rand(Uniform(-v̄/100, +v̄/100), N, N))
#r01, x, y, dx, dy = init_cond(:NoisePatches, N, c01 = ū, dx = dx)
#r02, x, y, dx, dy = init_cond(:NoisePatches, N, c01 = v̄, dx = dx)
#r01, x, y, dx, dy = init_cond(:ErMnO3_MnOnly, N, M = 9, α=1/80, BC=:Periodic, c01 =ū, dx = dx)
#r02, x, y, dx, dy = init_cond(:ErMnO3_ErOnly, N, M = 15, α=1, BC=:Neumann, c01 =v̄, dx = dx)
#r01 .= 1 .- r01
r0 = cat(cu(r01), cu(r02), dims=3)
r0_cpu = cat(r01, r02, dims=3)
#=============================================================
Finite differences matrices
==============================================================#
println("Dealing with finite differencing.")
#Q = Neumann0BC(1.0, 1)
Q = PeriodicBC(Float64)
D2 = CenteredDifference(2, 2, 1.0, N)
L2 = D2*Q

A2 = cu(sparse(L2)[1])
A2_cpu = sparse(L2)[1]
#============================================================
Caches
=============================================================#
println("Dealing with caches.")
A2u = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
uA2 = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
D2u = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
A2v = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
vA2 = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
D2v = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
A2u_i = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
u_iA2 = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
u_i = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
A2v_i = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
v_iA2 = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
v_i = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
u_ni = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))
v_ni = DiffEqBase.dualcache(cu(zeros(Float32, N,N)))

cache =(A2u, uA2, D2u, A2v, vA2, D2v, A2u_i, u_iA2, u_i, A2v_i, v_iA2, v_i)
p = (D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, cu(γ_x), cu(γ_y), dx, dy, A2, cache)
p_cpu = (D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, dx, dy, A2_cpu, cache)
#==========================================================
Advanced versoin of the PDE, for solving
===========================================================#
function advanced_version!(dr,r,p,t)
  D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, γ_x, γ_y, dx, dy, A2, cache = p
  A2u, uA2, D2u, A2v, vA2, D2v, A2u_i, u_iA2, u_i, A2v_i, v_iA2, v_i =  cache

  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]

  A2u = get_tmp(A2u, 0.f0)
  uA2 = get_tmp(uA2, 0.f0)
  D2u = get_tmp(D2u, 0.f0)
  A2v = get_tmp(A2v, 0.f0)
  vA2 = get_tmp(vA2, 0.f0)
  D2v = get_tmp(D2v, 0.f0)
  A2u_i = get_tmp(A2u_i, 0.f0)
  u_iA2 = get_tmp(u_iA2, 0.f0)
  u_i = get_tmp(u_i, 0.f0)
  A2v_i = get_tmp(A2v_i, 0.f0)
  v_iA2 = get_tmp(v_iA2, 0.f0)
  v_i = get_tmp(v_i, 0.f0)

  mul!(A2u,A2,u)
  mul!(uA2,u,cu(A2'))
  mul!(A2v,A2,v)
  mul!(vA2,v,cu(A2'))
  @. u_i = 2.0f0*u*(u-1)*(2.0f0*u-1) - γᵤ*(γ_x*1.0f0/dy^2*A2u + γ_y*1.0f0/dx^2*uA2) # 4 allocations
  @. v_i = 2.0f0*v*(v-1)*(2.0f0*v-1) - γᵥ*(γ_x*1.0f0/dy^2*A2v + γ_y*1.0f0/dx^2*vA2) # 4 allocatoins
  mul!(A2u_i,A2,u_i)
  mul!(u_iA2,u_i,cu(A2'))
  mul!(A2v_i,A2,v_i)
  mul!(v_iA2,v_i,cu(A2'))

  @. D2u = δ*D*(1.0f0/dy^2*A2u_i + 1.0f0/dx^2*u_iA2)
  @. D2v = δ*(1.0f0/dy^2*A2v_i + 1.0f0/dx^2*v_iA2)
  @. du = D2u + α*u*(1.0f0-r₁*v^2) + v*(1.0f0 - r₂*u)
  @. dv = D2v + β*v*(1.0f0+α*r₁/β*u*v) + u*(γ + r₂*v)
  nothing
end
#=============================================================
Basic Version for computing Jacobian sparsity
==============================================================#
function basic_version!(dr,r,p,t)
  D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, dx, dy, A2, cache = p

  u = @view r[:,:,1]
  v = @view r[:,:,2]

  u_n = 2u.*(u.-1).*(2u.-1) .- γᵤ*(1/dy^2*A2*u + 1/dx^2*u*A2')
  v_n = 2v.*(v.-1).*(2v.-1) .- γᵥ*(1/dy^2*A2*v + 1/dx^2*v*A2')
  D2u = D*δ*(1/dy^2*A2*u_n + 1/dx^2*u_n*A2')
  D2v = δ*(1/dy^2*A2*v_n + 1/dx^2*v_n*A2')
  dr[:,:,1] = D2u .+ α*u.*(1 .- r₁*v.^2) .+ v.*(1 .- r₂*u)
  dr[:,:,2] = D2v .+ β*v.*(1 .+ α*r₁/β*u.*v) .+ u.*(γ .+ r₂*v)
  nothing
end
#=================================================================
Compute the Jacobian sparsity pattern
==================================================================#
#println("Jacobian sparsity")
#dr0 = copy(r0_cpu)
#@time jac_sparsity = Symbolics.jacobian_sparsity((dr,r)->basic_version!(dr,r,p_cpu,0.f0),dr0,r0_cpu)

#jac_sparsity = load("j500_n0.jld", "jac_sparsity")

#jac = Float32.(sparse(jac_sparsity))
#colors = matrix_colors(jac)
#using FiniteDiff
#sparsecache = FiniteDiff.JacobianCache(r0_cpu, colorvec=colors,sparsity=jac)
#my_basic_version!(dr0, r0) = basic_version!(dr0, r0, p_cpu, 0.f0)
#FiniteDiff.finite_difference_jacobian!(jac, my_basic_version!, r0_cpu, sparsecache, colorvec=colors)
#using Arpack
#@time jac_eig = eigs(jac, nev=2, maxiter=500, ritzvec=false)
#@show jac_eig[1]

#f = ODEFunction(advanced_version!;jac=jac,jac_prototype=Float32.(jac_sparsity))

#==================================================================
Solve the system
==================================================================#
println("Solving")
prob = ODEProblem(advanced_version!,r0,tspan,p)
CUDA.allowscalar(false)
#@time sol = solve(prob, TRBDF2(linsolve=KrylovJL_GMRES(), precs=incompletelu, concrete_jac=true), saveat=range(0.0, stop=tspan[2], length=101), progress=true, progress_steps=1);
@time sol = solve(prob, ROCK2(), saveat=range(0.0, stop=tspan[2], length=101), progress=true, progress_steps=1, reltol=1e-6, abstol=1e-6);

#==================================================================
Plot
==================================================================#
#println("Plotting")
anim = @animate for i in 1:length(sol.t)
  title = "t = $(sol.t[i])"
  #u = kron(ones(3,3), sol.u[i][:,:,1])
  u = Array(sol.u[i])[:,:,1]
  heatmap(x, y, u, c=:grays, aspect_ratio=1, axis=([], false), title=title, colorbar=false)
  #heatmap(u, c=:roma, aspect_ratio=dy/dx, axis=([], false), title=title, colorbar=false)
end

gif(anim, fps=10)


text:
using BioFetch, FASTX, GenomicAnnotations, Test

nr002726 = fetchseq("NR_002726.2")
@test length(nr002726) == 1
@test typeof(nr002726[1]) == FASTA.Record

nr002726gb = fetchseq("NR_002726.2", format = gb)
@test length(nr002726gb) == 1
@test typeof(nr002726gb[1]) == GenomicAnnotations.Record{Gene}

seqs = fetchseq(["CAA41295.1", "NP_000176"], format = fasta)
@test length(seqs) == 2
@test typeof(seqs[1]) == FASTA.Record

nr002726 = fetchseq("NR_002726.2", 10:20)
nr002726rev = fetchseq("NR_002726.2", 10:20, false)
@test length(nr002726) == 1
@test typeof(nr002726[1]) == FASTA.Record
@test length(FASTX.sequence(nr002726[1])) == 11
@test length(FASTX.sequence(nr002726rev[1])) == 11
@test map(x -> x == 'G' ? 'C' : x == 'C' ? 'G' : x == 'T' ? 'A' : x == 'A' ? 'T' : x, FASTX.sequence(nr002726[1])) == FASTX.sequence(nr002726rev[1])

@test length(fetchseq("ENSG00000141510")) == 42


text:
Base.show(io::IO, ::MIME"text/html", org::Org) = (html(io, org), nothing)

const html_entity_map =
    Dict('&' => "&",
         '<' => "<",
         '>' => ">",
         '"' => """,
         '\'' => "'",
         '/' => "/",
         '`' => "`",
         '=' => "=")

html_escape(s) = replace(s, r"[&<>\"'`=\/]" =>
    e -> html_entity_map[e[1]])

function html_tagpair(tag, attrs::Pair...)
    attrstr = if length(attrs) > 0
        ' ' * join(map(a -> string(a.first, "=\"", html_escape(a.second), '"'), attrs), ' ')
    else "" end
    (string('<', tag, attrstr, '>'), string("</", tag, '>'))
end

function html_tagwrap(s, tag, escape=false, attrs::Pair...)
    open, close = html_tagpair(tag, attrs...)
    string(open, if escape html_escape(s) else s end, close)
end
html_tagwrap(s, tag, attrs::Pair...) = html_tagwrap(s, tag, false, attrs...)

# ---------------------
# Org
# ---------------------

function html(io::IO, o::Org)
    for component in o.contents
        html(io, component)
        print(io, '\n')
    end
end

html(o::Union{Org, OrgComponent}) = html(stdout, o)

function html(_::IO, component::C) where { C <: OrgComponent}
    @warn "No method for converting $C to a html representation currently exists"
end

# ---------------------
# Sections
# ---------------------

function html(io::IO, heading::Heading)
    sopen, sclose = if isnothing(heading.section); ("", "")
    else ("
", "
") end
    print(io, sopen, "<h$(heading.level)>")
    if !isnothing(heading.keyword)
        print(io, html_tagwrap(heading.keyword, "span", true, "class" => "hkeyword"), ' ')
    end
    if !isnothing(heading.priority)
        print(io, html_tagwrap(heading.priority, "span", true, "class" => "hpriority"), ' ')
    end
    html.(Ref(io), heading.title)
    if length(heading.tags) > 0
        print.(io, html_tagwrap.(heading.tags, "span", true, "class" => "htag"))
    end
    print(io, "</h$(heading.level)>")
    if !isnothing(heading.section)
        print(io, "\n")
        html(io, heading.section)
    end
    print(io, sclose)
end

function html(io::IO, section::Section)
    for component in section.contents
        html(io, component)
        print(io, '\n')
    end
end

# ---------------------
# Greater Elements
# ---------------------

# Greater Block
# Drawer
# Dynamic Block
# FootnoteDefinition
# InlineTask

function html(io::IO, item::Item)
    open, close = html_tagpair("li")
    print(io, open)
    # TODO support counterset
    # TODO support checkbox better
    if !isnothing(item.checkbox)
        print(io, " [", item.checkbox, "] ")
    end
    if !isnothing(item.tag)
        print(io, html_escape(item.tag), " :: ")
    end
    if !isnothing(item.contents)
        html.(Ref(io), item.contents)
    end
    if !isnothing(item.sublist)
        html(io, item.sublist)
    end
    print(io, close)
end

function html(io::IO, list::List)
    open, close = html_tagpair(if list isa OrderedList; "ol" else "ul" end)
    print(io, open, '\n')
    for item in list.items
        html(io, item)
        print(io, '\n')
    end
    print(io, close)
end

# PropertyDrawer
# Table

# ---------------------
# Elements
# ---------------------

# Babel Call
# Block
# Clock
# Planning
# DiarySexp
# Comment
# Fixed Width
# Horizontal Rule
# Keyword
# LaTeX Environment
# Node Property

function html(io::IO, par::Paragraph)
    print(io, "
")
    for obj in par.contents
        html(io, obj)
    end
    print(io, "
")
end

# Table Row

# ---------------------
# Objects
# ---------------------

html(io::IO, entity::Entity) = print(io, Entities[entity.name].html)

# LaTeX Fragment

html(io::IO, snippet::ExportSnippet) =
    if snippet.backend == "html" print(io, snippet.snippet) end

# Footnote Ref
# Inline Babel Call
# Inline Source Block
# Line Break
# Link
# Macro
# Radio Target

html(io::IO, target::Target) =
    print(io, html_tagwrap("", "a", "href" => string('#', hash(target.target))))

html(io::IO, statscookie::StatisticsCookiePercent) =
    print(io, '[', statscookie.percentage, "%]")
html(io::IO, statscookie::StatisticsCookieFraction) =
    print(io, '[', if !isnothing(statscookie.complete) string(statscookie.complete) else "" end,
          '/', if !isnothing(statscookie.total) string(statscookie.total) else "" end, ']')

html(io::IO, script::Superscript) =
    print(io, html_escape(script.char), html_tagwrap(script.script, "sup", true))
html(io::IO, script::Subscript) =
    print(io, html_escape(script.char), html_tagwrap(script.script, "sub", true))

# Table Cell

function html(io::IO, tsrod::TimestampRepeaterOrDelay)
    timeunits = Dict('h' => "hour",
                     'd' => "day",
                     'w' => "week",
                     'm' => "month",
                     'y' => "year")
    if tsrod.type in (:cumulative, :catchup, :restart)
        print(io, "every ", tsrod.value, ' ', timeunits[tsrod.unit],
              if tsrod.value == 1 "" else 's' end, " thereafter")
    else
        print(io, "warning ", tsrod.value, ' ', timeunits[tsrod.unit],
              if tsrod.value == 1 "" else 's' end, " before")
    end
end

function html(io::IO, ts::TimestampInstant)
    if isnothing(ts.time)
        print(io,
              html_tagwrap(Dates.format(ts.date, dateformat"e d u Y"),
                           "time", "datetime" => Dates.format(ts.date, ISODateFormat)))
    else
        print(io,
              html_tagwrap(Dates.format(ts.date, dateformat"H:M e d u Y"),
                           "time", "datetime" => Dates.format(DateTime(ts.date, ts.time),
                                                              ISODateTimeFormat)))
    end
    if !isnothing(ts.repeater)
        print(io, " and ")
        html(io, ts.repeater)
    end
    if !isnothing(ts.warning)
        printstyled(io, ", ")
        html(io, ts.warning)
    end
end

const html_markup_codes =
    Dict(:bold => html_tagpair("b"),
         :italic => html_tagpair("em"),
         :strikethrough => html_tagpair("s"),
         :underline => html_tagpair("span", "style" => "text-decoration: underline"),
         :verbatim => html_tagpair("kbd"),
         :code => html_tagpair("code"))

function html(io::IO, markup::TextMarkup)
    tagopen, tagclose = html_markup_codes[markup.formatting]
    print(io, tagopen)
    if markup.contents isa AbstractString
        print(io, html_escape(markup.contents))
    else
        html.(Ref(io), markup.contents)
    end
    print(io, tagclose)
end

function html(io::IO, text::TextPlain)
    tsub = replace(replace(replace(replace(html_escape(text.text),
                                           "..." => "…"),
                                   r"---([^-])" => s"—\1"),
                           r"--([^-])" => s"—\1"),
                   r"\\-" => "")
    print(io, tsub)
end


text:
df = CSV.File("config/min_max_adj.csv") |> DataFrame!
const MINS = Array{Float64}(df[!, :min])
const MAXS = Array{Float64}(df[!, :adj_max])
const DIFF = MAXS .- MINS
const ON_INP = .~(Bool.(df[!, :off]))

function get_inputs(features::Array{Float64})
    inputs = (features .- MINS) ./ DIFF
    inputs = min.(max.(inputs, -1.0), 1.0)
    inputs[ON_INP]
end