text:
#=
The commands below are a subset of the file read_exact developed for matlab.  It aims to
expose similar functionality, but at this time only for reading files, not for
reading buffers directly from the Gantner system.  It uses giutility.dll on
Windows machines to do the work.

This work was inspired from https://nextcloud.gantner-instruments.com/s/fGismwHwQynM9rE
=#

# location of shared library is in local directory
const giutility = joinpath(@__DIR__, "giutility.dll")

"""
    gantOpenFile(filename::String)
Open the gantner .dat file
"""
function gantOpenFile(filename::String) #, HCLIENT, HCONNECTION)
    resultDict = Dict{Int32, String}(
        0 => "OK",
        1 => "ERROR",
        2 => "CONNECTION_ERROR",
        3 => "INIT_ERROR",
        4 => "LIMIT_ERROR",
        5 => "SYNC_CONF_ERROR",
        6 => "MULTYUSED_ERROR",
        7 => "INDEX_ERROR",
        8 => "FILE_ERROR",
        9 => "NOT_READY",
        10 => "EXLIB_MISSING",
        11 => "NOT_CONNECTED",
        12 => "NO_FILE",
        13 => "CORE_ERROR",
        14 => "POINTER_INVALID")

    HCLIENT = HCONNECTION = Ref{Int32}(0)  # Initalize values

    ret = ccall((:_CD_eGateHighSpeedPort_DecodeFile_Select, giutility),
    Cint,(Ref{Cint},Ref{Cint},Cstring),HCLIENT,HCONNECTION,filename)

    # @show(filename, HCLIENT, HCONNECTION, ret)

    if ret != 0
        error("File did not open, Return Code $ret meaning $(resultDict[ret])")
    end
    return HCLIENT[], HCONNECTION[]
end

"""
    gantCloseFile(HCLIENT::Integer, HCONNECTION::Integer)
Close the file that has been opened.  And probably data streams as well.
"""
function gantCloseFile(HCLIENT::Integer, HCONNECTION::Integer)

    channelCount = Ref{Int32}(-1)

    ret = ccall((:_CD_eGateHighSpeedPort_Close, giutility),
    Cint,(Cint,Cint),HCLIENT, HCONNECTION)

    HCLIENT = HCONNECTION = Int32(-1)  # Initalize values

    # @show(HCONNECTION, channelCount, ret)

  if ret != 0
    error("Channel numbers not read, Return Code $ret")
  end
  return HCLIENT[], HCONNECTION[]
end

"""
    gantChanNumRead(HCONNECTION::Integer)
Get the number of channels
"""
function gantChanNumRead(HCONNECTION::Integer)

    DADI_INPUT = 0   # input channels
    channelCount = Ref{Int32}(0)

    ret = ccall((:_CD_eGateHighSpeedPort_GetNumberOfChannels, giutility),
    Cint,(Cint,Cint,Ref{Cint}),HCONNECTION,DADI_INPUT,channelCount)

    # @show(HCONNECTION, channelCount, ret)

  if ret != 0
    error("Channel numbers not read, Return Code $ret")
  end
  return channelCount[]
end

"""
    gantNumSamples(HCLIENT::Integer, HCONNECTION::Integer)
Get the number of samples per channel
"""
function gantNumSamples(HCLIENT::Integer, HCONNECTION::Integer)
    # read the number of samples per channel in the file
    numValues = ccall((:_CD_eGateHighSpeedPort_GetBufferFrames_All, giutility),
    Cint, (Cint, Cint), HCONNECTION, HCLIENT)

   return numValues
end

"""
    gantSampleRate(HCONNECTION::Integer)
Get the sample rate
"""
function gantSampleRate(HCONNECTION::Integer)

    fs = Ref{Float64}(0)
    channelInfo = ""

    ret = ccall((:_CD_eGateHighSpeedPort_GetDeviceInfo, giutility),
    Cint,(Cint,Cint,Cint,Ref{Cdouble},Cstring),HCONNECTION,16,1,fs,channelInfo)

    # @show(HCONNECTION, fs, channelInfo, ret)

  if ret != 0
    error("Channel sample rate not read, Return Code $ret")
  end
  return fs[]
end

"""
    gantChanName(HCONNECTION::Integer, channelindex::Integer)
Get the channel name
"""
function gantChanName(HCONNECTION::Integer, channelindex::Integer)

    CHINFO_NAME = 0  # type of info = name
    DADI_INPUT = 0   # input channels
    info = Vector{UInt8}(undef, 1024)  # initialize info

    ret = ccall((:_CD_eGateHighSpeedPort_GetChannelInfo_String, giutility),
    Cint,(Cint,Cint,Cint,Cint,Ptr{Cchar}),
    HCONNECTION,CHINFO_NAME,DADI_INPUT,channelindex-1,info)

    # @show(HCONNECTION, channelindex, info, ret)

  if ret != 0
    error("Channel sample rate not read, Return Code $ret")
  end
  return unsafe_string(pointer(info))
end

"""
    gantChanDataRead(HCLIENT::Integer, HCONNECTION::Integer, channelindex::Integer)
Get the data, using logic of read_exact.c
This logic seems convoluted but think it is a workaround to ideosyncracies in
the giutility.dll file.
"""
function gantChanDataRead(HCLIENT::Integer, HCONNECTION::Integer, channelindex::Integer)
    HSPdict = Dict{Int,String}(
    0 => "HSP_OK",
    1 => "HSP_ERROR",
    2 => "HSP_CONNECTION_ERROR",
    3 => "HSP_INIT_ERROR",
    4 => "HSP_LIMIT_ERROR",
    5 => "HSP_SYNC_CONF_ERROR",
    6 => "HSP_MULTYUSED_ERROR",
    7 => "HSP_INDEX_ERROR",
    8 => "HSP_FILE_ERROR",
    9 => "HSP_NOT_READY")

    HSP_OK = 0

    numChan = gantChanNumRead(HCONNECTION)
    #@show(numChan)

    if 1 > channelindex || channelindex > numChan
        error("Channel $channelindex is not in the file")
    end

    # read_exact gets the sample rate and timestep

    # read the number of values in the file
    numValues = ccall((:_CD_eGateHighSpeedPort_GetBufferFrames_All, giutility),
    Cint, (Cint, Cint), HCONNECTION, HCLIENT)
    #@show(numValues)

    if numValues == 0
        error("No Data in file...")
    end
    Data = zeros(numValues)

    frameSize = Ref{Int32}(0)
    ret = ccall((:_CD_eGateHighSpeedPort_LoadBufferData, giutility),
    Cint, (Cint, Ref{Cint}), HCONNECTION, frameSize)
    index = 1
    found = 0
    finished = 0

    while ret == HSP_OK
        value = Ref{Float64}(0)

        # Load next value from buffer
        for i = 1:frameSize[]
            # read channel value
            ret = ccall((:_CD_eGateHighSpeedPort_ReadBuffer_Single, giutility),
            Cint, (Cint,Cint,Cint, Ref{Cdouble}), HCONNECTION,HCLIENT,channelindex-1, value)

            if ret != HSP_OK
                #display("error at read buffer single... CHind: $channelindex returncode: $(HSPdict[ret])")
                found = 1
                break
            end

            # @show(value[])
            Data[index] = value[]
            index += 1

            # load next value
            ret = ccall((:_CD_eGateHighSpeedPort_ReadBuffer_NextFrame, giutility),
            Cint, (Cint, Cint), HCONNECTION, HCLIENT)
            if ret != HSP_OK
                #@show(i, ret)
                #display("Error at next frame.... returncode: $ret, $(HSPdict[ret])")
                break
            end

        end
        if found == 1
            break
        end

        # load next data -> if no data available: ret = 1
        ret = ccall((:_CD_eGateHighSpeedPort_LoadBufferData, giutility),
        Cint, (Cint, Ref{Cint}), HCONNECTION, frameSize)

        if ret != HSP_OK
            finished = 1
        end
    end
    if finished == 1
        ret = 0
    end

    # rewind the file
    ret = ccall((:_CD_eGateHighSpeedPort_Rewind, giutility),
    Cint, (Cint, Cint, Cint), HCONNECTION, HCLIENT, numValues)
    return Data
end

#=
"""
    gantChanDataRead(HCLIENT::Integer, HCONNECTION::Integer, channelindex::Integer)
Get the data.  This was an attempt to make the logic from read_exact more
straightforward.  However the return codes from giutility.dll are strange so
it gives funny error messages that I have not fully figured out.
"""
function gantChanDataRead(HCLIENT::Integer, HCONNECTION::Integer, channelindex::Integer)
    HSPdict = Dict{Int,String}(
    0 => "HSP_OK",
    1 => "HSP_ERROR",
    2 => "HSP_CONNECTION_ERROR",
    3 => "HSP_INIT_ERROR",
    4 => "HSP_LIMIT_ERROR",
    5 => "HSP_SYNC_CONF_ERROR",
    6 => "HSP_MULTYUSED_ERROR",
    7 => "HSP_INDEX_ERROR",
    8 => "HSP_FILE_ERROR",
    9 => "HSP_NOT_READY")

    HSP_OK = 0

    numChan = gantChanNumRead(HCONNECTION)

    @show(numChan)

    if 1 > channelindex || channelindex > numChan
        error("Channel $channelindex is not in the file")
    end

    # read_exact gets the sample rate and timestep

    # read the number of values in the file
    numValues = ccall((:_CD_eGateHighSpeedPort_GetBufferFrames_All, "./giutility.dll"),
    Cint, (Cint, Cint), HCONNECTION, HCLIENT)

    if numValues == 0
        error("No Data in file...")
    end

    @show(numValues)

    Data = zeros(numValues)

    value = Ref{Float64}(0)
    frameSize = 0  # Initialize to keep in scope
    for i = 1:numValues
        # Load next value from buffer
        if i == 1 || i % frameSize == 2 # works because i == 1 evaluated first
            frameSize = Ref{Int32}(0)
            ret = ccall((:_CD_eGateHighSpeedPort_LoadBufferData, "./giutility.dll"),
            Cint, (Cint, Ref{Cint}), HCONNECTION, frameSize)
            if ret != HSP_OK
                display("error at load next data -> if no data available ret = 1, ret: $ret")
            end
            frameSize = frameSize[]
            @show(i,frameSize)
        else
            ret = ccall((:_CD_eGateHighSpeedPort_ReadBuffer_NextFrame, "./giutility.dll"),
            Cint, (Cint, Cint), HCONNECTION, HCLIENT)
            if ret != HSP_OK
                @show(i, ret)
                display("Error at next frame.... returncode: $ret, $(HSPdict[ret])")
            end
        end
        # read channel value
        ret = ccall((:_CD_eGateHighSpeedPort_ReadBuffer_Single, "./giutility.dll"),
        Cint, (Cint,Cint,Cint, Ref{Cdouble}), HCONNECTION,HCLIENT,channelindex-1, value)

        if ret != HSP_OK
            display("error at read buffer single... CHind: $channelindex returncode: $(HSPdict[ret])")
        end

        # @show(value[])

        Data[i] = value[]

    end
    @show(i)
    # load next data -> if no data available: ret = 1

    # rewind the file
    ret = ccall((:_CD_eGateHighSpeedPort_Rewind, "./giutility.dll"),
    Cint, (Cint, Cint, Cint), HCONNECTION, HCLIENT, numValues)
    return Data
end
=#


text:
if !@isdefined(TestINet)
    struct TestINet <: ANetwork; end
    struct TestIDiNet <: ADiNetwork; end
    struct TestIIdxEdge <: AIndexedEdge; end
end

@testset "$TEST $G" begin

g = TestINet()
h = TestIDiNet()
e = TestIIdxEdge()

@test_throws ErrorException idx(e)

if G<: ANetwork

g = G(10,20)
m = ones(10,10)
em = add_edge_property!(g, "rand", m)
@test typeof(em) <: EdgeMap
@test em.data === m
@test eprop(g, "rand")[1,1] == 1

m = fill(2, 10)
vm = add_vertex_property!(g, "rand", m)
@test typeof(vm) <: VertexMap
@test vm.data === m
@test vprop(g, "rand")[1] == 2

end #if

end # testset

ext: jl


text:
using HTTP
using Test
using Pluto
using Pluto: ServerSession, ClientSession, SessionActions
using Pluto.Configuration
using Pluto.Configuration: notebook_path_suggestion, from_flat_kwargs, _convert_to_flags

@testset "Configurations" begin

cd(Pluto.project_relative_path("test")) do
    @test notebook_path_suggestion() == joinpath(pwd(), "")
end

@testset "from_flat_kwargs" begin
    opt = from_flat_kwargs(;compile="min", launch_browser=false)
    @test opt.compiler.compile == "min"
    @test opt.server.launch_browser == false

    @test_throws ArgumentError from_flat_kwargs(;asdfasdf="test")
end

@testset "flag conversion" begin
    if VERSION > v"1.5.0-"
        @test _convert_to_flags(Configuration.CompilerOptions(threads="123")) ==
            ["--startup-file=no", "--history-file=no", "--threads=123"]

        @test _convert_to_flags(Configuration.CompilerOptions(threads=123)) ==
            ["--startup-file=no", "--history-file=no", "--threads=123"]

        @test _convert_to_flags(Configuration.CompilerOptions()) ⊇
            ["--startup-file=no", "--history-file=no"]
    else
        @test _convert_to_flags(Configuration.CompilerOptions()) ==
            ["--startup-file=no", "--history-file=no"]
    end
    @test _convert_to_flags(Configuration.CompilerOptions(compile="min")) ⊇
    ["--compile=min", "--startup-file=no", "--history-file=no"]
end

@testset "authentication" begin
    port = 1238
    options = Pluto.Configuration.from_flat_kwargs(; port=port, launch_browser=false, workspace_use_distributed=false)
    🍭 = Pluto.ServerSession(; options=options)
    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient
    host = 🍭.options.server.host
    secret = 🍭.secret
    println("Launching test server...")
    server_task = @async Pluto.run(🍭)
    sleep(2)

    local_url(suffix) = "http://$host:$port/$suffix"
    withsecret(url) = occursin('?', url) ? "$url&secret=$secret" : "$url?secret=$secret"
    @test HTTP.get(local_url("favicon.ico")).status == 200

    function requeststatus(url, method)
        r = HTTP.request(method, url; status_exception=false, redirect=false)
        r.status
    end

    nb = SessionActions.open(🍭, Pluto.project_relative_path("sample", "Basic.jl"); as_sample=true)

    simple_routes = [
        ("", "GET"),
        ("edit?id=$(nb.notebook_id)", "GET"),
        ("notebookfile?id=$(nb.notebook_id)", "GET"),
        ("notebookexport?id=$(nb.notebook_id)", "GET"),
        ("statefile?id=$(nb.notebook_id)", "GET"),
    ]

    function tempcopy(x)
        p = tempname()
        Pluto.readwrite(x, p)
        p
    end
    @assert isfile(Pluto.project_relative_path("sample", "Basic.jl"))

    effect_routes = [
        ("new", "GET"),
        ("new", "POST"),
        ("open?url=$(HTTP.URIs.escapeuri("https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.14.5/sample/Basic.jl"))", "GET"),
        ("open?url=$(HTTP.URIs.escapeuri("https://raw.githubusercontent.com/fonsp/Pluto.jl/v0.14.5/sample/Basic.jl"))", "POST"),
        ("open?path=$(HTTP.URIs.escapeuri(Pluto.project_relative_path("sample", "Basic.jl") |> tempcopy))", "GET"),
        ("open?path=$(HTTP.URIs.escapeuri(Pluto.project_relative_path("sample", "Basic.jl") |> tempcopy))", "POST"),
        ("sample/Basic.jl", "GET"),
        ("sample/Basic.jl", "POST"),
        ("notebookupload", "POST"),
    ]

    for (suffix, method) in simple_routes ∪ effect_routes
        url = local_url(suffix)
        @test requeststatus(url, method) == 403
    end

    # no notebooks were opened
    @test length(🍭.notebooks) == 1

    for (suffix, method) in simple_routes
        url = local_url(suffix) |> withsecret
        @test requeststatus(url, method)  ∈ 200:299
    end

    for (suffix, method) in setdiff(effect_routes, [("notebookupload", "POST")])
        url = local_url(suffix) |> withsecret
        @test requeststatus(url, method) ∈ 200:399 # 3xx are redirects
    end

    @async schedule(server_task, InterruptException(); error=true)
end

end # testset


text:
export TraceTestSystem

"""
    TraceTestSystem(F, p₀, v) <: AbstractSystem

The system `G(x,l) := [F(x, p₀+x[end]v); l ⋅[x[1:end-1];1])`.
"""
struct TraceTestSystem{S<:AbstractSystem, T} <: AbstractSystem
    F::S
    p₀::Vector{T}
    v::Vector{T}
end

function TraceTestSystem(F, p::AbstractVector, v::AbstractVector)
    TraceTestSystem(F, promote(p, v)...)
end

struct TraceTestSystemCache{SC<:AbstractSystemCache, T1, T2} <: AbstractSystemCache
    u::Vector{T1}
    U::Matrix{T1}
    Up::Matrix{T1}
    p::Vector{T2}
    cache::SC
end

function cache(F::TraceTestSystem, x, l)
    c = cache(F.F, x, l)
    u = evaluate(F.F, x, l, c)
    U = similar(u, size(F.F))
    Up = similar(u, size(F.F)[1], length(F.p₀))
    p = F.p₀ + 0.1 * F.v
    TraceTestSystemCache(u, U, Up, p, c)
end

Base.size(F::TraceTestSystem) = size(F.F) .+ 1

function evaluate!(u, F::TraceTestSystem, x, l, c::TraceTestSystemCache)
    λ = x[end]
    for i in eachindex(c.p)
        c.p[i] = F.p₀[i] + λ * F.v[i]
    end
    evaluate!(c.u, F.F, x, c.p, c.cache)
    n = length(x)
    for i in 1:length(c.u)
        u[i] = c.u[i]
    end
    u[end] = l[n]
    for i in 1:n-1
        u[end] += l[i] * x[i]
    end

    u
end

function jacobian!(U, F::TraceTestSystem, x, l, c::TraceTestSystemCache)
    λ = x[end]
    for i in eachindex(c.p)
        c.p[i] = F.p₀[i] + λ * F.v[i]
    end
    jacobian!(c.U, F.F, x, c.p, c.cache)
    differentiate_parameters!(c.Up, F.F, x, c.p, c.cache)
    n, m = length(x), size(U, 1)
    for j in 1:n-1, i in 1:m-1
        U[i,j] = c.U[i,j]
    end
    for i in 1:m
        U[i,n] = zero(eltype(U))
    end
    for k in 1:size(c.Up, 2), i in 1:m-1
        U[i, n] += c.Up[i,k] * F.v[k]
    end

    for j in 1:n-1
        U[m,j] = l[j]
    end
    U
end

function evaluate_and_jacobian!(u, U, F::TraceTestSystem, x, l, c::TraceTestSystemCache)
    λ = x[end]
    for i in eachindex(c.p)
        c.p[i] = F.p₀[i] + λ * F.v[i]
    end
    n, m = length(x), size(U, 1)

    evaluate_and_jacobian!(c.u, c.U, F.F, x, c.p, c.cache)
    for i in 1:length(c.u)
        u[i] = c.u[i]
    end
    u[end] = l[n]
    for i in 1:n-1
        u[end] += l[i] * x[i]
    end

    differentiate_parameters!(c.Up, F.F, x, c.p, c.cache)

    for j in 1:n-1, i in 1:m-1
        U[i,j] = c.U[i,j]
    end
    for i in 1:m
        U[i,n] = zero(eltype(U))
    end
    for k in 1:size(c.Up, 2), i in 1:m-1
        U[i, n] += c.Up[i,k] * F.v[k]
    end

    for j in 1:n-1
        U[m,j] = l[j]
    end

    nothing
end

function differentiate_parameters!(U, F::TraceTestSystem, x, l, c::TraceTestSystemCache)
    U .= zero(eltype(U))
    for i in 1:length(l)-1
        U[end, i] = x[i]
    end
    U[end,length(l)] = one(eltype(U))
    U
end

function differentiate_parameters(F::TraceTestSystem, x, l, c::TraceTestSystemCache)
    differentiate_parameters!(similar(c.U, size(F)[1], length(l)), F, x, l, c)
end

function evaluate(F::TraceTestSystem, x, l, c::TraceTestSystemCache)
    evaluate!(similar(c.u, size(F)[1]), F, x, l, c)
end
function jacobian(F::TraceTestSystem, x, l, c::TraceTestSystemCache)
    jacobian!(similar(c.U, size(F)), F, x, l, c)
end


text:
module LowStorageRungeKuttaMethod
export LowStorageRungeKutta, updatedt!

using Requires

@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
  using .CuArrays
  using .CuArrays.CUDAnative
  using .CuArrays.CUDAnative.CUDAdrv

  include("LowStorageRungeKuttaMethod_cuda.jl")
end

using ...ODESolvers
ODEs = ODESolvers
using ...SpaceMethods

"""
    LowStorageRungeKutta(f, Q; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

    Q̇ = f(Q)

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This uses the fourth-order, low-storage, Runge--Kutta scheme of Carpenter
and Kennedy (1994) (in their notation (5,4) 2N-Storage RK scheme).

### References

    @TECHREPORT{CarpenterKennedy1994,
      author = {M.~H. Carpenter and C.~A. Kennedy},
      title = {Fourth-order {2N-storage} {Runge-Kutta} schemes},
      institution = {National Aeronautics and Space Administration},
      year = {1994},
      number = {NASA TM-109112},
      address = {Langley Research Center, Hampton, VA},
    }
"""
struct LowStorageRungeKutta{T, AT, Nstages} <: ODEs.AbstractODESolver
  "time step"
  dt::Array{T,1}
  "time"
  t::Array{T,1}
  "rhs function"
  rhs!::Function
  "Storage for RHS during the LowStorageRungeKutta update"
  dQ::AT
  "low storage RK coefficient vector A (rhs scaling)"
  RKA::NTuple{Nstages, T}
  "low storage RK coefficient vector B (rhs add in scaling)"
  RKB::NTuple{Nstages, T}
  "low storage RK coefficient vector C (time scaling)"
  RKC::NTuple{Nstages, T}
  function LowStorageRungeKutta(rhs!::Function, Q::AT; dt=nothing,
                                t0=0) where {AT<:AbstractArray}

    @assert dt != nothing

    T = eltype(Q)
    dt = [T(dt)]
    t0 = [T(t0)]
    # FIXME: Add reference
    RKA = (T(0),
           T(-567301805773)  / T(1357537059087),
           T(-2404267990393) / T(2016746695238),
           T(-3550918686646) / T(2091501179385),
           T(-1275806237668) / T(842570457699 ))

    RKB = (T(1432997174477) / T(9575080441755 ),
           T(5161836677717) / T(13612068292357),
           T(1720146321549) / T(2090206949498 ),
           T(3134564353537) / T(4481467310338 ),
           T(2277821191437) / T(14882151754819))

    RKC = (T(0),
           T(1432997174477) / T(9575080441755),
           T(2526269341429) / T(6820363962896),
           T(2006345519317) / T(3224310063776),
           T(2802321613138) / T(2924317926251))

    new{T, AT, length(RKA)}(dt, t0, rhs!, similar(Q), RKA, RKB, RKC)
  end
end

function LowStorageRungeKutta(spacedisc::AbstractSpaceMethod, Q; dt=nothing,
                              t0=0)
  LowStorageRungeKutta((x...) -> SpaceMethods.odefun!(spacedisc, x...),
                       Q; dt=dt, t0=t0)
end

"""
    updatedt!(lsrk::LowStorageRungeKutta, dt)

Change the time step size to `dt` for `lsrk.
"""
updatedt!(lsrk::LowStorageRungeKutta, dt) = lsrk.dt[1] = dt

function ODEs.dostep!(Q, lsrk::LowStorageRungeKutta, timeend,
                      adjustfinalstep)
  time, dt = lsrk.t[1], lsrk.dt[1]
  if adjustfinalstep && time + dt > timeend
    dt = timeend - time
    @assert dt > 0
  end
  RKA, RKB, RKC = lsrk.RKA, lsrk.RKB, lsrk.RKC
  rhs!, dQ = lsrk.rhs!, lsrk.dQ
  for s = 1:length(RKA)
    rhs!(dQ, Q, time)

    # update solution and scale RHS
    # FIXME: GPUify
    update!(Val(size(Q,2)), Val(size(Q,1)), dQ.Q, Q.Q, Q.realelems,
            RKA[s%length(RKA)+1], RKB[s], dt)

    time += RKC[s] * dt
  end
  if dt == lsrk.dt[1]
    lsrk.t[1] += dt
  else
    lsrk.t[1] = timeend
  end

end

# {{{ Update solution (for all dimensions)
function update!(::Val{nstates}, ::Val{Np}, rhs::Array{T, 3}, Q, elems, rka,
                 rkb, dt) where {nstates, Np, T}
  @inbounds for e = elems, s = 1:nstates, i = 1:Np
    Q[i, s, e] += rkb * dt * rhs[i, s, e]
    rhs[i, s, e] *= rka

  end
end
# }}}

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
export getBoundaryIndices, getBoundaryIndicesFEM

"""
function getBoundaryIndices(M::RegularMesh, bc=(:dir,:neu,:neu,:neu))

gets indices of boundary points, neighbors of boundary points and interior
points in extended mesh

Input:

	M::RegularMesh - mesh
	bc             - chooses Dirichlet or Neuman BC for each boundary
	                 Ex1: dim=2, domain=[0,1,0,1] --> bc is a vector of length 4
					      edge 1: x1=0, x2 in [0,1] (left)
						  edge 2: x2=0, x1 in [0,1] (bottom)
						  edge 3: x1=1, x2 in [0,1] (right)
						  edge 4: x2=1, x1 in [0,1] (top)
		             Ex1: dim=3, domain=[0,1,0,1,0,1] --> bc is a vector of length 6
		 			      face 1: x1=0, x2,x3 in [0,1]
		 			 	  face 2: x2=0, x1,x3 in [0,1]
		 			 	  face 3: x3=0, x1,x2 in [0,1]
		 			      face 4: x1=1, x2,x3 in [0,1]
		 			 	  face 5: x2=2, x1,x3 in [0,1]
		 			 	  face 6: x3=3, x1,x2 in [0,1]

Outputs:
	iddir     - indices of ghost points associated with Dirichlet conditions
	idneu     - indices of ghost points associated with Neuman conditions
	iddirn    - indices of neighbors of ghost points associated with Dirichlet conditions
	idneun    - indices of neighbors of ghost points associated with Neuman conditions
	idint     - indices of interior nodes

"""
function getBoundaryIndices(M::RegularMesh, bc=(:dir,:neu,:neu,:neu,:neu,:neu))

	# get boundary conditions
	idg,idgn,idint = getGhostIndices(M.n)
	iddir  = zeros(Int,0)
	idneu  = zeros(Int,0)
	iddirn = zeros(Int,0)
	idneun = zeros(Int,0)
	for k=1:length(bc)
		if bc[k]==:dir
			iddir  = vcat(iddir,idg[k])
			iddirn = vcat(iddirn, idgn[k])
		elseif bc[k]==:neu
			idneu  = vcat(idneu,idg[k])
			idneun = vcat(idneun,idgn[k])
		end
	end

	return iddir, idneu, iddirn, idneun, idint
end

"""
function getBoundaryIndicesFEM(M::RegularMesh, bc=(:dir,:neu,:neu,:neu))

Input:

	M::RegularMesh - mesh
	bc             - chooses Dirichlet or Neuman BC for each boundary
	                 Ex1: dim=2, domain=[0,1,0,1] --> bc is a vector of length 4
					      edge 1: x1=0, x2 in [0,1] (left)
						  edge 2: x2=0, x1 in [0,1] (bottom)
						  edge 3: x1=1, x2 in [0,1] (right)
						  edge 4: x2=1, x1 in [0,1] (top)
		             Ex1: dim=3, domain=[0,1,0,1,0,1] --> bc is a vector of length 6
		 			      face 1: x1=0, x2,x3 in [0,1]
		 			 	  face 2: x2=0, x1,x3 in [0,1]
		 			 	  face 3: x3=0, x1,x2 in [0,1]
		 			      face 4: x1=1, x2,x3 in [0,1]
		 			 	  face 5: x2=2, x1,x3 in [0,1]
		 			 	  face 6: x3=3, x1,x2 in [0,1]

Outputs:
	iddir     - indices of FEM nodes associated with Dirichlet conditions

"""
function getBoundaryIndicesFEM(M::RegularMesh, bc=(:dir,:neu,:neu,:neu,:neu,:neu))

	# get all nodes
	ids    = reshape(1:prod(M.n.+1),tuple(M.n.+1...))
	iddir  = zeros(Int,0)

	if M.dim==2
		if bc[1]==:dir
			iddir = vcat(iddir,vec(ids[1,:]));
		end
		if bc[2]==:dir
			iddir = vcat(iddir,vec(ids[:,1]));
		end
		if bc[3]==:dir
			iddir = vcat(iddir,vec(ids[end,:]));
		end
		if bc[4]==:dir
			iddir = vcat(iddir,vec(ids[:,end]));
		end
	elseif M.dim==3
		if bc[1]==:dir
			iddir = vcat(iddir,vec(ids[1,:,:]));
		end
		if bc[2]==:dir
			iddir = vcat(iddir,vec(ids[:,1,:]));
		end
		if bc[3]==:dir
			iddir = vcat(iddir,vec(ids[:,:,1]));
		end
		if bc[4]==:dir
			iddir = vcat(iddir,vec(ids[end,:,:]));
		end
		if bc[5]==:dir
			iddir = vcat(iddir,vec(ids[:,end,:]));
		end
		if bc[6]==:dir
			iddir = vcat(iddir,vec(ids[:,:,end]));
		end
	end
	return unique(iddir)
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
h = nothing
lsatx = LSAtT(0.0)

# Create the stopping object to test
stop = LS_Stopping(h, (x,y)-> armijo(x,y), lsatx, max_iter = 10)

# We tests different functions of stopping
OK = update_and_start!(stop, x = 1.0, g₀ = NaN, h₀ = NaN, ht = NaN)
@test OK == true
@test status(stop) == :DomainError
@test stop.current_state.x == 1.0

reinit!(stop)
@test stop.meta.domainerror == false
update!(lsatx, g₀ = 0.0, h₀ = 0.0, ht = 10.0)
@test (stop.current_state.g₀ == lsatx.g₀) && (lsatx.g₀ == 0.0)
@test (stop.current_state.h₀ == lsatx.h₀) && (lsatx.h₀ == 0.0)
@test (stop.current_state.ht == lsatx.ht) && (lsatx.ht == 10.0)
@test start!(stop) == false
@test status(stop) == :Unknown
#
@test update_and_stop!(stop, ht = -10.0) == true
@test status(stop) == :Optimal

reinit!(stop)
@test stop.meta.optimal == false
update!(stop.current_state, ht = 10.0)
# Check if _tired_check works
@test isnan(stop.meta.start_time) #default time value is NaN
start!(stop) #start initializes the start_time
@test !isnan(stop.meta.start_time)
@test !stop!(stop) #shall we stop?Not yet
#If we check stop! without a given start_time:
stop.meta.start_time = NaN
@test !(stop!(stop))
@test stop.meta.tired == false
@test status(stop) == :Unknown

# we can check _tired_check for max evals with NLPStopping

# Check if _unbounded_check! works
update!(stop.current_state, x = 1e100)
@test stop!(stop)

reinit!(stop, rstate = true, x = 1.0)
@test stop.current_state.x == 1.0
@test stop.current_state.ht == nothing

## _optimality_check and _null_test are tested with NLP
try
armijo(stop.pb, stop.current_state)
@test false #nothing entries in the stop
catch
@test true
end
try
wolfe(stop.pb, stop.current_state)
@test false #nothing entries in the stop
catch
@test true
end
try
armijo_wolfe(stop.pb, stop.current_state)
@test false #nothing entries in the stop
catch
@test true
end
update!(stop.current_state, h₀ = 1.0, ht = 0.0, g₀ = 1.0, gt = 0.0)
@test wolfe(stop.pb, stop.current_state) == 0.0
@test armijo_wolfe(stop.pb, stop.current_state) == 0.0
mutable struct Tpb
    d :: Number
end
stop.pb = Tpb(0.0)
@test shamanskii_stop(stop.pb, stop.current_state) == 0.0 #specific LineModel
@test goldstein(stop.pb, stop.current_state) >= 0.0

stop.optimality_check = (x,y) -> 0.0
stop.pb = ADNLPModel(x -> 0.0, [1.0])
stop.meta.max_f = -1
reinit!(stop.current_state, 0.0)
@test stop.current_state.ht == nothing
@test stop!(stop) == true
@test stop.current_state.ht == nothing
@test stop.meta.resources == true


text:
using ProgressMeter
using MATLAB
using Printf
include("PSS_function.jl")
# H=range(.12,stop=0.144281,length=211);
# H=[0.112012,.125789, .135]
H=range(.1,stop=0.15,length=10);
max_hit=100
width=.03
height=.4
h_start=replace(@sprintf("%.3f",H[1]),"."=>"_")
h_end=replace(@sprintf("%.3f",H[end]),"."=>"_")

location="~/Desktop/PSS_MOVIE/"
# video_name=location*"PSS_MOVIE_"*h_start*"to"*h_end*".avi"
# mat"writerObj = VideoWriter($video_name);"

# mat"""
# writerObj.FrameRate=15;
# writerObj.Quality = 100;
# open(writerObj);
mat"""
figure('position',[0, 0, 1200, 900]);
hold on;
pbaspect([10 5 1])
"""
mat"axis([-$width, $width, -$height, $height])"

n_iter_Q=11;#50
Q_start=-.25
Q_end=.25
n_iter_P=10;#51
P_start=-.015
P_end=.015
t_end=1e2

H_differences=zeros(0)

ArrP=range(P_start,stop=P_end,length=n_iter_P)
ArrQ=range(Q_start,stop=Q_end,length=n_iter_Q)

#defining colors for PSS
COLOR=["b", "g", "c", "m", "y", "k", "r"]

@showprogress 1 "Computing..." for Energy in H
    h=replace(@sprintf("%.8f",Energy),"."=>"_")
    file_name=location*"PSS_"*h*".fig"
    n=0
    mat"hold on;"
    mat"pbaspect([10 5 1])"
    mat"axis([-$width, $width, -$height, $height])"
    for j=1:n_iter_P
        for k=1:n_iter_Q
            Q,P,dH=PSS_function(ArrQ[k], ArrP[j], Energy, t_end, max_hit)
            if dH<5e-5 #make sure dH is ok
                current_color=COLOR[mod(k,length(COLOR))+1]
                # mat"plot($P_PSS,$Q_PSS,'.','MarkerSize',.1,'color',$current_color); hold on;"
                mat"plot($P,$Q,'k.','MarkerSize',.1);"
            end
        end
    end
    mat"savefig($file_name)"
    #
    # mat"""
    # try
    #     thisFrame = getframe;
    # end
    # myMovie($n) = thisFrame;
    # """
    # if n>1
    #     mat"""
    #     try
    #         writeVideo(writerObj, thisFrame);
    #     end
    #     """
    # end
    mat"cla reset"
    # n=n+1
end
# mat"save('test.mat',myMovie)"
# mat"close(writerObj)"


text:
"""
    MathOptLMO{OT <: MOI.Optimizer} <: LinearMinimizationOracle

Linear minimization oracle with feasible space defined through a MathOptInterface.Optimizer.
The oracle call sets the direction and reruns the optimizer.

The `direction` vector has to be set in the same order of variables as the `MOI.ListOfVariableIndices()` getter.

The Boolean `use_modify` determines if the objective in`compute_extreme_point` is updated with
`MOI.modify(o, ::MOI.ObjectiveFunction, ::MOI.ScalarCoefficientChange)` or with `MOI.set(o, ::MOI.ObjectiveFunction, f)`.
`use_modify = true` decreases the runtime and memory allocation for models created as an optimizer object and defined directly
with MathOptInterface. `use_modify = false` should be used for CachingOptimizers.
"""
struct MathOptLMO{OT<:MOI.AbstractOptimizer} <: LinearMinimizationOracle
    o::OT
    use_modify::Bool
    function MathOptLMO(o, use_modify=true)
        MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        return new{typeof(o)}(o, use_modify)
    end
end

function compute_extreme_point(lmo::MathOptLMO{OT}, direction::AbstractVector{T}; kwargs...) where {OT,T<:Real}
    variables = MOI.get(lmo.o, MOI.ListOfVariableIndices())
    if lmo.use_modify
        for i in eachindex(variables)
            MOI.modify(
                lmo.o,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
                MOI.ScalarCoefficientChange(variables[i], direction[i]),
            )
        end
    else
        terms = [MOI.ScalarAffineTerm(d, v) for (d, v) in zip(direction, variables)]
        obj = MOI.ScalarAffineFunction(terms, zero(T))
        MOI.set(lmo.o, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    end
    return _optimize_and_return(lmo, variables)
end

function compute_extreme_point(lmo::MathOptLMO{OT}, direction::AbstractMatrix{T}; kwargs...) where {OT,T<:Real}
    n = size(direction, 1)
    v = compute_extreme_point(lmo, vec(direction))
    return reshape(v, n, n)
end

function Base.copy(lmo::MathOptLMO{OT}; ensure_identity=true) where {OT}
    opt = OT() # creates the empty optimizer
    index_map = MOI.copy_to(opt, lmo.o)
    if ensure_identity
        for (src_idx, des_idx) in index_map.var_map
            if src_idx != des_idx
                error("Mapping of variables is not identity")
            end
        end
    end
    return MathOptLMO(opt)
end

function Base.copy(lmo::MathOptLMO{OT}; ensure_identity=true) where {OTI, OT <: MOIU.CachingOptimizer{OTI}}
    opt = MOIU.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        OTI(),
    )
    index_map = MOI.copy_to(opt, lmo.o)
    if ensure_identity
        for (src_idx, des_idx) in index_map.var_map
            if src_idx != des_idx
                error("Mapping of variables is not identity")
            end
        end
    end
    return MathOptLMO(opt)
end

function compute_extreme_point(
    lmo::MathOptLMO{OT},
    direction::AbstractVector{MOI.ScalarAffineTerm{T}};
    kwargs...
) where {OT,T}
    if lmo.use_modify
        for d in direction
            MOI.modify(
                lmo.o,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
                MOI.ScalarCoefficientChange(d.variable,d.coefficient),
            )
        end

        variables = MOI.get(lmo.o, MOI.ListOfVariableIndices())
        variables_to_zero = setdiff(variables, [dir.variable for dir in direction])

        terms = [MOI.ScalarAffineTerm(d, v) for (d, v) in zip(zeros(length(variables_to_zero)), variables_to_zero)]

        for t in terms
            MOI.modify(
                lmo.o,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
                MOI.ScalarCoefficientChange(t.variable,t.coefficient),
            )
        end
    else
        variables = [d.variable for d in direction]
        obj = MOI.ScalarAffineFunction(direction, zero(T))
        MOI.set(lmo.o, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    end
    return _optimize_and_return(lmo, variables)
end

function _optimize_and_return(lmo, variables)
    MOI.optimize!(lmo.o)
    term_st = MOI.get(lmo.o, MOI.TerminationStatus())
    if term_st ∉ (MOI.OPTIMAL, MOI.ALMOST_OPTIMAL, MOI.SLOW_PROGRESS)
        @error "Unexpected termination: $term_st"
        return MOI.get.(lmo.o, MOI.VariablePrimal(), variables)
    end
    return MOI.get.(lmo.o, MOI.VariablePrimal(), variables)
end

"""
    convert_mathopt(lmo::LMO, optimizer::OT; kwargs...) -> MathOptLMO{OT}

Converts the given LMO to its equivalent MathOptInterface representation using `optimizer`.
Must be implemented by LMOs.
"""
function convert_mathopt end


text:
### A Pluto.jl notebook ###
# v0.9.11

using Markdown
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.peek, el) ? Base.peek(el) : missing
        el
    end
end

# ╔═╡ 5d93901c-ba6d-11ea-354f-658440cfd6fb
using Plots

# ╔═╡ e40d8e1c-ba5d-11ea-0057-f3ac3801c047
# Simple model of infection

# ╔═╡ c2c32a88-ba5e-11ea-0c75-b17710c4421b
## Each individual infects c new people each day

# ╔═╡ 121f4b78-ba5f-11ea-1abe-1b53c89d2475
I₀ = 1 # \_0 <TAB>

# ╔═╡ b7ebc04c-ba5f-11ea-0456-05cb3654c89d
c = 0.01 # average number of people each individual infects per day

# ╔═╡ eecb614a-ba5f-11ea-268f-e3879450b557
λ = 1 + c # \lambda <TAB>

# ╔═╡ 308091de-ba60-11ea-0491-efc040c0b2fe
I₁ = λ * I₀

# ╔═╡ 625fe844-ba60-11ea-35f6-a13d8310aa9b
I₂ = λ * I₁

# ╔═╡ 7a23dcb0-ba60-11ea-2dd6-ade51d455d8e
# Arrays

# ╔═╡ 91603fa4-ba60-11ea-139f-e3983d0cd261
begin
	T = 10
	I = zeros(T)
end

# ╔═╡ f18523d6-ba60-11ea-298b-e77c44717b3c
I[1] = I₀

# ╔═╡ b17a49d2-ba61-11ea-3dd8-e1a0175fd8ee
for n in 1:T-1
	I[n+1] = λ * I[n]
end

# ╔═╡ 494cc8ba-ba6d-11ea-2c4d-0d9bd08659f5
I

# ╔═╡ a8cae8e4-ba6d-11ea-1e0a-cfed38c2f685
plot(I, m=:o, label="I[n]", leg=:topleft)

# ╔═╡ 10683e04-ba6e-11ea-18a8-e11794ea52b1
# Functions

# ╔═╡ fc894d30-ba6e-11ea-3af0-5f1289515858
function run_infection(T=20)
	I = zeros(T)
	I[1] = I₀

	for n ∈ 1:T-1
		I[n+1] = λ * I[n]
	end

	return I
end

# ╔═╡ 6da88666-ba6f-11ea-07e3-c990bcaf93f1
methods(run_infection)

# ╔═╡ 80d5f08c-ba71-11ea-198d-d392e0143e10
max = 1000

# ╔═╡ 2f2ba61a-ba70-11ea-34ad-f30ca470d61b
@bind range_num html"<input type='range' min=1 max=1000>"

# ╔═╡ 812bd2d8-ba6f-11ea-3424-319f7b22b82e
I_result = run_infection(range_num)

# ╔═╡ a397fb74-ba6f-11ea-1586-e19d4898924e
begin
	plot(I_result, m=:o)

	xlims!(0, max)
	ylims!(0, 10)
end

# ╔═╡ 0e22838a-ba70-11ea-0473-49bd2932698f
# Exponential growth

# ╔═╡ d917e7aa-ba71-11ea-0b45-a791fde442f8
plot(run_infection(1000))

# ╔═╡ e9e54262-ba71-11ea-2a5d-add1bc1230de
plot(run_infection(1000), yscale=:log10)

# ╔═╡ 2f692bfa-ba72-11ea-28a3-558290ff0b73
begin
	I_result_log = run_infection(1000)

	plot(log.(I_result_log))
end

# ╔═╡ a95a44a8-ba72-11ea-06a4-b108eb9b0549
begin
	p = 0.02
	α = 0.01
	N = 1000

	β(I, S) = p * α * (N - I)

	function run_infection_logistic(T=20)
		I = zeros(T)
		I[1] = I₀

		for n ∈ 1:T-1
			I[n+1] = I[n] + β(I[n], N - I[n]) * I[n]
		end

		return I
	end
end

# ╔═╡ 5b0010bc-ba77-11ea-1c2f-d71f73bbc813
ril = run_infection_logistic()

# ╔═╡ 8376e7dc-ba77-11ea-0d7f-d35b8a611b23
plot(ril, m=:o)

# ╔═╡ a4dc902a-ba77-11ea-32bf-4767238e9828
plot(ril, m=:o, yscale=:log10)

# ╔═╡ b7f5591c-ba77-11ea-117f-cfbab8eb8f3d
plot(run_infection_logistic(100), m=:o) #sigmoid shape

# ╔═╡ df4c9c8a-ba77-11ea-2d73-9f3a35f7e21c
plot(run_infection_logistic(100), m=:o, yscale=:log10)

# ╔═╡ 388a4f9e-ba78-11ea-16a7-3f597ef95fef
# More realism: Hetreogeneity of individuals

# ╔═╡ 99e9557c-ba7c-11ea-2afd-e1767f6ee8fc
rand()

# ╔═╡ f8f18c9c-ba7c-11ea-18a0-c98ac848c314
randn()

# ╔═╡ 0d855362-ba7a-11ea-3f82-4bc5f6c382c8
begin
	function run_infection_realistic(T=20, c_average=1.1)
		I = zeros(T)
		I[1] = I₀

		for n ∈ 1:T-1
			c = c_average + 0.1 * rand()
			I[n+1] = I[n] +c * I[n]
		end

		return I
	end
end

# ╔═╡ c718542a-ba7d-11ea-0532-f10f96230741
c_average = 1.1

# ╔═╡ e51af200-ba7d-11ea-2cf7-39b1b93b872b
cs = [c_average + 0.1 * randn() for _ ∈ 1:100]

# ╔═╡ f5d2ca8e-ba7d-11ea-3d0b-89123de86872
scatter(cs)

# ╔═╡ 6dc1c28e-ba7e-11ea-0826-39089532e11e
plot(run_infection_realistic(100, 0.02))

# ╔═╡ Cell order:
# ╠═e40d8e1c-ba5d-11ea-0057-f3ac3801c047
# ╠═c2c32a88-ba5e-11ea-0c75-b17710c4421b
# ╠═121f4b78-ba5f-11ea-1abe-1b53c89d2475
# ╠═b7ebc04c-ba5f-11ea-0456-05cb3654c89d
# ╠═eecb614a-ba5f-11ea-268f-e3879450b557
# ╠═308091de-ba60-11ea-0491-efc040c0b2fe
# ╠═625fe844-ba60-11ea-35f6-a13d8310aa9b
# ╠═7a23dcb0-ba60-11ea-2dd6-ade51d455d8e
# ╠═91603fa4-ba60-11ea-139f-e3983d0cd261
# ╠═f18523d6-ba60-11ea-298b-e77c44717b3c
# ╠═b17a49d2-ba61-11ea-3dd8-e1a0175fd8ee
# ╠═494cc8ba-ba6d-11ea-2c4d-0d9bd08659f5
# ╠═5d93901c-ba6d-11ea-354f-658440cfd6fb
# ╠═a8cae8e4-ba6d-11ea-1e0a-cfed38c2f685
# ╠═10683e04-ba6e-11ea-18a8-e11794ea52b1
# ╠═fc894d30-ba6e-11ea-3af0-5f1289515858
# ╠═6da88666-ba6f-11ea-07e3-c990bcaf93f1
# ╠═80d5f08c-ba71-11ea-198d-d392e0143e10
# ╠═2f2ba61a-ba70-11ea-34ad-f30ca470d61b
# ╠═812bd2d8-ba6f-11ea-3424-319f7b22b82e
# ╠═a397fb74-ba6f-11ea-1586-e19d4898924e
# ╠═0e22838a-ba70-11ea-0473-49bd2932698f
# ╠═d917e7aa-ba71-11ea-0b45-a791fde442f8
# ╠═e9e54262-ba71-11ea-2a5d-add1bc1230de
# ╠═2f692bfa-ba72-11ea-28a3-558290ff0b73
# ╠═a95a44a8-ba72-11ea-06a4-b108eb9b0549
# ╠═5b0010bc-ba77-11ea-1c2f-d71f73bbc813
# ╠═8376e7dc-ba77-11ea-0d7f-d35b8a611b23
# ╠═a4dc902a-ba77-11ea-32bf-4767238e9828
# ╠═b7f5591c-ba77-11ea-117f-cfbab8eb8f3d
# ╠═df4c9c8a-ba77-11ea-2d73-9f3a35f7e21c
# ╠═388a4f9e-ba78-11ea-16a7-3f597ef95fef
# ╠═99e9557c-ba7c-11ea-2afd-e1767f6ee8fc
# ╠═f8f18c9c-ba7c-11ea-18a0-c98ac848c314
# ╠═0d855362-ba7a-11ea-3f82-4bc5f6c382c8
# ╠═c718542a-ba7d-11ea-0532-f10f96230741
# ╠═e51af200-ba7d-11ea-2cf7-39b1b93b872b
# ╠═f5d2ca8e-ba7d-11ea-3d0b-89123de86872
# ╠═6dc1c28e-ba7e-11ea-0826-39089532e11e

max_stars_repo_name: elihugarret/my-notebooks


text:
import ForwardDiff: Dual, partials, value

@static if isdefined(LinearAlgebra, :gradient)
    import LinearAlgebra.gradient
end

######################
# Extraction methods #
######################

# Extractions are supposed to unpack the value and the partials
# The partials should be put into a tensor of higher order.
# The extraction methods need to know the input type to the function
# that generated the result. The reason for this is that there is no
# difference in the output type (except the number of partials) for
# norm(v) and det(T) where v is a vector and T is a second order tensor.

####################
# Value extraction #
####################

# Scalar output -> Scalar value
@inline function _extract_value(v::Dual)
    return value(v)
end
# AbstractTensor output -> AbstractTensor gradient
@generated function _extract_value(v::AbstractTensor{<:Any,<:Any,<:Dual})
    TensorType = get_base(v)
    ex = Expr(:tuple)
    for i in 1:n_components(TensorType)
        # Can use linear indexing even for SymmetricTensor
        # when indexing the underlying tuple
        push!(ex.args, :(value(get_data(v)[$i])))
    end
    quote
        $(Expr(:meta, :inline))
        @inbounds return $TensorType($ex)
    end
end

#######################
# Gradient extraction #
#######################

# Scalar output, Scalar input -> Scalar gradient
@inline function _extract_gradient(v::Dual, ::Number)
    return @inbounds partials(v)[1]
end
# Vec, Tensor{2/4}, SymmetricTensor{2/4} output, Scalar input -> Vec, Tensor{2/4}, SymmetricTensor{2/4} gradient
@generated function _extract_gradient(v::AbstractTensor{<:Any,<:Any,<:Dual}, ::Number)
    TensorType = get_base(v)
    ex = Expr(:tuple)
    for i in 1:n_components(TensorType)
        # Can use linear indexing even for SymmetricTensor
        # when indexing the underlying tuple
        push!(ex.args, :(partials(get_data(v)[$i])[1]))
    end
    quote
        $(Expr(:meta, :inline))
        @inbounds return $TensorType($ex)
    end
end
# Scalar output, Vec input -> Vec gradient
@inline function _extract_gradient(v::Dual, ::Vec{N}) where {N}
    return Vec{N}(partials(v).values)
end
# Scalar output, SecondOrderTensor input -> SecondOrderTensor gradient
@inline function _extract_gradient(v::Dual, ::SymmetricTensor{2, dim}) where {dim}
    return SymmetricTensor{2, dim}(partials(v).values)
end
@inline function _extract_gradient(v::Dual, ::Tensor{2, dim}) where {dim}
    return Tensor{2, dim}(partials(v).values)
end

# Vec output, Vec input -> Tensor{2} gradient
@inline function _extract_gradient(v::Vec{1, <: Dual}, ::Vec{1})
    @inbounds begin
        p1 = partials(v[1])
        ∇f = Tensor{2, 1}((p1[1],))
    end
    return ∇f
end
# Vec output, Vec input -> Tensor{2} gradient
@inline function _extract_gradient(v::Vec{2, <: Dual}, ::Vec{2})
    @inbounds begin
        p1, p2 = partials(v[1]), partials(v[2])
        ∇f = Tensor{2, 2}((p1[1], p2[1], p1[2], p2[2]))
    end
    return ∇f
end
# Vec output, Vec input -> Tensor{2} gradient
@inline function _extract_gradient(v::Vec{3, <: Dual}, ::Vec{3})
    @inbounds begin
        p1, p2, p3 = partials(v[1]), partials(v[2]), partials(v[3])
        ∇f = Tensor{2, 3}((p1[1], p2[1], p3[1], p1[2], p2[2], p3[2], p1[3], p2[3], p3[3]))
    end
    return ∇f
end

# Tensor{2} output, Tensor{2} input -> Tensor{4} gradient
@inline function _extract_gradient(v::Tensor{2, 1, <: Dual}, ::Tensor{2, 1})
    @inbounds begin
        p1 = partials(v[1,1])
        ∇f = Tensor{4, 1}((p1[1],))
    end
    return ∇f
end
# SymmetricTensor{2} output, SymmetricTensor{2} input -> SymmetricTensor{4} gradient
@inline function _extract_gradient(v::SymmetricTensor{2, 1, <: Dual}, ::SymmetricTensor{2, 1})
    @inbounds begin
        p1 = partials(v[1,1])
        ∇f = SymmetricTensor{4, 1}((p1[1],))
    end
    return ∇f
end
# Tensor{2} output, Tensor{2} input -> Tensor{4} gradient
@inline function _extract_gradient(v::Tensor{2, 2, <: Dual}, ::Tensor{2, 2})
    @inbounds begin
        p1, p2, p3, p4 = partials(v[1,1]), partials(v[2,1]), partials(v[1,2]), partials(v[2,2])
        ∇f = Tensor{4, 2}((p1[1], p2[1], p3[1], p4[1],
                           p1[2], p2[2], p3[2], p4[2],
                           p1[3], p2[3], p3[3], p4[3],
                           p1[4], p2[4], p3[4], p4[4]))
    end
    return ∇f
end
# SymmetricTensor{2} output, SymmetricTensor{2} input -> SymmetricTensor{4} gradient
@inline function _extract_gradient(v::SymmetricTensor{2, 2, <: Dual}, ::SymmetricTensor{2, 2})
    @inbounds begin
        p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[2,2])
        ∇f = SymmetricTensor{4, 2}((p1[1], p2[1], p3[1],
                                    p1[2], p2[2], p3[2],
                                    p1[3], p2[3], p3[3]))
    end
    return ∇f
end
# Tensor{2} output, Tensor{2} input -> Tensor{4} gradient
@inline function _extract_gradient(v::Tensor{2, 3, <: Dual}, ::Tensor{2, 3})
    @inbounds begin
        p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[3,1])
        p4, p5, p6 = partials(v[1,2]), partials(v[2,2]), partials(v[3,2])
        p7, p8, p9 = partials(v[1,3]), partials(v[2,3]), partials(v[3,3])
        ∇f = Tensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1], p7[1], p8[1], p9[1],
                           p1[2], p2[2], p3[2], p4[2], p5[2], p6[2], p7[2], p8[2], p9[2],  #    ###  #
                           p1[3], p2[3], p3[3], p4[3], p5[3], p6[3], p7[3], p8[3], p9[3],  #    # #  #
                           p1[4], p2[4], p3[4], p4[4], p5[4], p6[4], p7[4], p8[4], p9[4],  ###  ###  ###
                           p1[5], p2[5], p3[5], p4[5], p5[5], p6[5], p7[5], p8[5], p9[5],
                           p1[6], p2[6], p3[6], p4[6], p5[6], p6[6], p7[6], p8[6], p9[6],
                           p1[7], p2[7], p3[7], p4[7], p5[7], p6[7], p7[7], p8[7], p9[7],
                           p1[8], p2[8], p3[8], p4[8], p5[8], p6[8], p7[8], p8[8], p9[8],
                           p1[9], p2[9], p3[9], p4[9], p5[9], p6[9], p7[9], p8[9], p9[9]))
    end
    return ∇f
end
# SymmetricTensor{2} output, SymmetricTensor{2} input -> SymmetricTensor{4} gradient
@inline function _extract_gradient(v::SymmetricTensor{2, 3, <: Dual}, ::SymmetricTensor{2, 3})
    @inbounds begin
        p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[3,1])
        p4, p5, p6 = partials(v[2,2]), partials(v[3,2]), partials(v[3,3])
        ∇f = SymmetricTensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1],
                                    p1[2], p2[2], p3[2], p4[2], p5[2], p6[2],
                                    p1[3], p2[3], p3[3], p4[3], p5[3], p6[3],
                                    p1[4], p2[4], p3[4], p4[4], p5[4], p6[4],
                                    p1[5], p2[5], p3[5], p4[5], p5[5], p6[5],
                                    p1[6], p2[6], p3[6], p4[6], p5[6], p6[6]))
    end
    return ∇f
end

# for non dual variable
@inline function _extract_value(v::Any)
    return v
end
for TensorType in (Tensor, SymmetricTensor)
    @eval begin
        @inline function _extract_gradient(v::T, x::$TensorType{order, dim}) where {T<:Real, order, dim}
            zero($TensorType{order, dim, T})
        end
        @generated function _extract_gradient(v::$TensorType{order, dim, T}, ::$TensorType{order, dim}) where {T<:Real, order, dim}
            RetType = $TensorType{order+order, dim, T}
            return quote
                $(Expr(:meta, :inline))
                zero($RetType)
            end
        end
    end
end

##################
# Load functions #
##################

# Loaders are supposed to take a tensor of real values and convert it
# into a tensor of dual values where the seeds are correctly defined.

@inline function _load(v::Number)
    return Dual(v, one(v))
end

@inline function _load(v::Vec{1, T}) where {T}
    @inbounds v_dual = Vec{1}((Dual(v[1], one(T)),))
    return v_dual
end

@inline function _load(v::Vec{2, T}) where {T}
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Vec{2}((Dual(v[1], o, z),
                               Dual(v[2], z, o)))
    return v_dual
end

@inline function _load(v::Vec{3, T}) where {T}
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Vec{3}((Dual(v[1], o, z, z),
                               Dual(v[2], z, o, z),
                               Dual(v[3], z, z, o)))
    return v_dual
end

# Second order tensors
@inline function _load(v::Tensor{2, 1, T}) where {T}
    @inbounds v_dual = Tensor{2, 1}((Dual(get_data(v)[1], one(T)),))
    return v_dual
end

@inline function _load(v::SymmetricTensor{2, 1, T}) where {T}
    @inbounds v_dual = SymmetricTensor{2, 1}((Dual(get_data(v)[1], one(T)),))
    return v_dual
end

@inline function _load(v::Tensor{2, 2, T}) where {T}
    data = get_data(v)
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Tensor{2, 2}((Dual(data[1], o, z, z, z),
                                     Dual(data[2], z, o, z, z),
                                     Dual(data[3], z, z, o, z),
                                     Dual(data[4], z, z, z, o)))
    return v_dual
end

@inline function _load(v::SymmetricTensor{2, 2, T}) where {T}
    data = get_data(v)
    o = one(T)
    o2 = convert(T, 1/2)
    z = zero(T)
    @inbounds v_dual = SymmetricTensor{2, 2}((Dual(data[1], o, z, z),
                                              Dual(data[2], z, o2, z),
                                              Dual(data[3], z, z, o)))
    return v_dual
end

@inline function _load(v::Tensor{2, 3, T}) where {T}
    data = get_data(v)
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Tensor{2, 3}((Dual(data[1], o, z, z, z, z, z, z, z, z),
                                     Dual(data[2], z, o, z, z, z, z, z, z, z),
                                     Dual(data[3], z, z, o, z, z, z, z, z, z),
                                     Dual(data[4], z, z, z, o, z, z, z, z, z),
                                     Dual(data[5], z, z, z, z, o, z, z, z, z),
                                     Dual(data[6], z, z, z, z, z, o, z, z, z),
                                     Dual(data[7], z, z, z, z, z, z, o, z, z),
                                     Dual(data[8], z, z, z, z, z, z, z, o, z),
                                     Dual(data[9], z, z, z, z, z, z, z, z, o)))
    return v_dual
end

@inline function _load(v::SymmetricTensor{2, 3, T}) where {T}
    data = get_data(v)
    o = one(T)
    o2 = convert(T, 1/2)
    z = zero(T)
    @inbounds v_dual = SymmetricTensor{2, 3}((Dual(data[1], o, z, z, z, z, z),
                                              Dual(data[2], z, o2, z, z, z, z),
                                              Dual(data[3], z, z, o2, z, z, z),
                                              Dual(data[4], z, z, z, o, z, z),
                                              Dual(data[5], z, z, z, z, o2, z),
                                              Dual(data[6], z, z, z, z, z, o)))
    return v_dual
end

"""
    gradient(f::Function, v::Union{SecondOrderTensor, Vec, Number})
    gradient(f::Function, v::Union{SecondOrderTensor, Vec, Number}, :all)

Computes the gradient of the input function. If the (pseudo)-keyword `all`
is given, the value of the function is also returned as a second output argument.

# Examples
```jldoctest
julia> A = rand(SymmetricTensor{2, 2});

julia> ∇f = gradient(norm, A)
2×2 SymmetricTensor{2,2,Float64,3}:
 0.434906  0.56442
 0.56442   0.416793

julia> ∇f, f = gradient(norm, A, :all);
```
"""
function gradient(f::F, v::Union{SecondOrderTensor, Vec, Number}) where {F}
    v_dual = _load(v)
    res = f(v_dual)
    return _extract_gradient(res, v)
end
function gradient(f::F, v::Union{SecondOrderTensor, Vec, Number}, ::Symbol) where {F}
    v_dual = _load(v)
    res = f(v_dual)
    return _extract_gradient(res, v), _extract_value(res)
end
const ∇ = gradient

"""
    hessian(f::Function, v::Union{SecondOrderTensor, Vec, Number})
    hessian(f::Function, v::Union{SecondOrderTensor, Vec, Number}, :all)

Computes the hessian of the input function. If the (pseudo)-keyword `all`
is given, the lower order results (gradient and value) of the function is
also returned as a second and third output argument.

# Examples
```jldoctest
julia> A = rand(SymmetricTensor{2, 2});

julia> ∇∇f = hessian(norm, A)
2×2×2×2 SymmetricTensor{4,2,Float64,9}:
[:, :, 1, 1] =
  0.596851  -0.180684
 -0.180684  -0.133425

[:, :, 2, 1] =
 -0.180684   0.133546
  0.133546  -0.173159

[:, :, 1, 2] =
 -0.180684   0.133546
  0.133546  -0.173159

[:, :, 2, 2] =
 -0.133425  -0.173159
 -0.173159   0.608207

julia> ∇∇f, ∇f, f = hessian(norm, A, :all);
```
"""
function hessian(f::F, v::Union{SecondOrderTensor, Vec, Number}) where {F}
    gradf = y -> gradient(f, y)
    return gradient(gradf, v)
end

function hessian(f::F, v::Union{SecondOrderTensor, Vec, Number}, ::Symbol) where {F}
    gradf = y -> gradient(f, y)
    return gradient(gradf, v), gradient(f, v, :all)...
end
const ∇∇ = hessian

"""
    divergence(f, x)

Calculate the divergence of the vector field `f`, in the point `x`.

# Examples
```jldoctest
julia> f(x) = 2x;

julia> x = rand(Vec{3});

julia> divergence(f, x)
6.0
```
"""
divergence(f::F, v::Vec) where {F<:Function} = tr(gradient(f, v))

"""
    curl(f, x)

Calculate the curl of the vector field `f`, in the point `x`.

# Examples
```jldoctest
julia> f(x) = Vec{3}((x[2], x[3], -x[1]));

julia> x = rand(Vec{3});

julia> curl(f, x)
3-element Tensor{1,3,Float64,3}:
 -1.0
  1.0
 -1.0
```
"""
function curl(f::F, v::Vec{3}) where F
    @inbounds begin
        ∇f = gradient(f, v)
        c = Vec{3}((∇f[3,2] - ∇f[2,3], ∇f[1,3] - ∇f[3,1], ∇f[2,1] - ∇f[1,2]))
    end
    return c
end
curl(f::F, v::Vec{1, T}) where {F, T} = curl(f, Vec{3}((v[1], T(0), T(0))))
curl(f::F, v::Vec{2, T}) where {F, T} = curl(f, Vec{3}((v[1], v[2], T(0))))

"""
    laplace(f, x)

Calculate the laplacian of the field `f`, in the point `x`.
If `f` is a vector field, use broadcasting.

# Examples
```jldoctest
julia> x = rand(Vec{3});

julia> f(x) = norm(x);

julia> laplace(f, x)
1.7833701103136868

julia> g(x) = x*norm(x);

julia> laplace.(g, x)
3-element Tensor{1,3,Float64,3}:
 2.107389336871036
 2.7349658311504834
 2.019621767876747
```
"""
function laplace(f::F, v) where F
    return divergence(x -> gradient(f, x), v)
end
const Δ = laplace

function Broadcast.broadcasted(::typeof(laplace), f::F, v::Vec{3}) where {F}
    @inbounds begin
        vdd = _load(_load(v))
        res = f(vdd)
        v1 = res[1].partials[1].partials[1] + res[1].partials[2].partials[2] + res[1].partials[3].partials[3]
        v2 = res[2].partials[1].partials[1] + res[2].partials[2].partials[2] + res[2].partials[3].partials[3]
        v3 = res[3].partials[1].partials[1] + res[3].partials[2].partials[2] + res[3].partials[3].partials[3]
    end
    return Vec{3}((v1, v2, v3))
end


text:
using Base.Dates

const FredKeyTerms = [(:Nominal, :Real, :Spot),
                      (:ConstantMaturity, (:Nominal, :Real), (:Daily, :Weekly, :Monthly)),
                      (:Forward, (:InstantRate, :TermPremium), (:Daily,)),
                      (:BreakevenInflation, (:InstantRate, :TermPremium), (:Daily, :Monthly)),
                      (:ZeroCoupon, (:InstantRate, :InstantYield), (:Daily)) ]

"""
    TreasuryRates
       (:ConstantMaturity   ==> (:Nominal, :Real) => (:Daily, :Weekly, :Monthly) -> (Day(n), Month(n)),
       (:Forward            ==> ()),
       (:BreakevenInflation ==> (:InstantRate, :TermPremium)  => (:Daily, :Monthly) : (Day(n), Month(n), Year(n)),
       (:ZeroCoupon         ==> (:InstantRate, :InstantYield) => (:Daily))
    IndexSeries
        (:TradeWeighted  => (:Daily, :Weekly, :Monthly))
    SpotSeries
        (:Gold => (:Morning[:USD,:EURO], :Afternoon[:USD.:EURO]) ),
        (:CrudeOil => (:WestTexas[:USD], :Brent[:USD]) ),
        (:NaturalGas => :USD) )
"""

const FredTreasuryRates = Dict([
  :ConstantMaturity => Dict([
    :Nominal => Dict([
              :Daily => Dict([
                       Year(30) => "DGS30", Year(20) => "DGS20",  Year(10) => "DGS10", Year(7) => "DGS7",
                       Year(5) => "DGS5",  Year(3) => "DGS3", Year(2) => "DGS2",  Year(1) => "DGS1",
                       Month(6) => "DGS6MO", Month(3) => "DGS3MO", Month(1) => "DGS1MO",
                       Day(1) => "DGS1DA",
              ]), # :Daily
              :Weekly => Dict([
                       Year(30) => "WGS30", Year(20) => "WGS20",  Year(10) => "WGS10", Year(7) => "WGS7",
                       Year(5) => "WGS5",  Year(3) => "WGS3", Year(2) => "WGS2",  Year(1) => "WGS1",
                       Month(6) => "WGS6MO", Month(3) => "WGS3MO", Month(1) => "WGS1MO",
                       Day(1) => "WGS1DA",
              ]), # :Weekly
              :Monthly => Dict([
                       Year(30) => "MGS30", Year(20) => "MGS20",  Year(10) => "MGS10", Year(7) => "MGS7",
                       Year(5) => "MGS5",  Year(3) => "MGS3", Year(2) => "MGS2",  Year(1) => "WGS1",
                       Month(6) => "MGS6MO", Month(3) => "MGS3MO", Month(1) => "MGS1MO",
                       Day(1) => "MGS1DA",
              ]), # :Monthly
          ]),   # :Nominal
   :Real =>  Dict([
              :Daily => Dict([
                       Year(30) => "DFII30", Year(20) => "DFII20",  Year(10) => "DFII10",
                       Year(7) => "DFII7",  Year(5) => "DFII5",
              ]), # :Daily
              :Weekly => Dict([
                      Year(30) => "WFII30", Year(20) => "WFII20",  Year(10) => "WFII10",
                      Year(7) => "WFII7",  Year(5) => "WFII5",
              ]), # :Weekly
              :Monthly => Dict([
                      Year(30) => "NFII30", Year(20) => "MFII20",  Year(10) => "MFII10",
                      Year(7) => "MFII7",  Year(5) => "MFII5",
              ]), # :Monthly
          ]),  # :Real
  ]),  # ConstantMaturity

   :BreakevenInflation => Dict([
             :Daily => Dict([
               Year(30) => "T30YIEM", Year(20) => "T20YIEM",  Year(10) => "T10YIEM",  Year(5) => "T5YIEM"
             ]), # :Daily
             :Monthly => Dict([
               Year(30) => "T30YIEM", Year(20) => "T20YIEM",  Year(10) => "T10YIEM", Year(5) => "T5YIEM"
             ]), # :Monthly
           ]),   # :BreakevenInflationRate

   :Forward => Dict([
        :InstantRate => Dict([
               :Daily => Dict([
                      Year(10) => "THREEFF10", Year(9) => "THREEFF9",  Year(8) => "THREEFF8",  Year(7) => "THREEFF7",
                      Year(6) => "THREEFF6",  Year(5) => "THREEFF5",  Year(4) => "THREEFF4",  Year(3) => "THREEFF3",
                      Year(2) => "THREEFF2", Year(1) => "THREEFF1"
                ]), # :Daily
             ]), # :InstantRate
       :TermPremium => Dict([
               :Daily => Dict([
                      Year(10) => "THREEFFTP10", Year(9) => "THREEFFTP9",  Year(8) => "THREEFFTP8",  Year(7) => "THREEFFTP7",
                      Year(6) => "THREEFFTP6",  Year(5) => "THREEFFTP5",  Year(4) => "THREEFFTP4",  Year(3) => "THREEFFTP3",
                      Year(2) => "THREEFFTP2", Year(1) => "THREEFFTP1"
                ]), # :Daily
             ]), # :TermPremium
     ]),         # :Forward

  :ZeroCoupon => Dict([
        :InstantYield => Dict([
               :Daily => Dict([
                      Year(10) => "THREEFY10", Year(9) => "THREEFY9",  Year(8) => "THREEFY8",  Year(7) => "THREEFY7",
                      Year(6) => "THREEFY6",  Year(5) => "THREEFY5",  Year(4) => "THREEFY4",  Year(3) => "THREEFY3",
                      Year(2) => "THREEFY2", Year(1) => "THREEFY1" ]),
               ]), # :Daily
       :TermPremium => Dict([
               :Daily => Dict([
                      Year(10) => "THREEFYTP10", Year(9) => "THREEFYTP9",  Year(8) => "THREEFYTP8",  Year(7) => "THREEFYTP7",
                      Year(6) => "THREEFYTP6",  Year(5) => "THREEFYTP5",  Year(4) => "THREEFYTP4",  Year(3) => "THREEFYTP3",
                      Year(2) => "THREEFYTP2", Year(1) => "THREEFYTP1" ]),
               ]), # :Daily
     ]),           # :ZeroCoupon

  ]); # FredTreasuryRates

const FredIndexedSeries = Dict([

       :TradeWeightedUSD => Dict([
          :Daily => Dict([
            :Major => "DTWEXM", :Broad => "DTWEXB"
          ]), # :Daily
          :Weekly => Dict([
            :Major => "TWEXM", :Broad => "TWEXMB"
           ]), # :Weekly
          :Monthly => Dict([
             :Nominal => Dict([ :Major => "TWEXMMTH", :Broad => "TWEXBMTH" ]),
             :Real    => Dict([ :Major => "TWEXMPA",  :Broad => "TWEXBPA"  ]),
           ]), # :Monthly
       ]),    # :TradeWeightedUSD

  ]); # FredIndexedSeries

const FredSpotSeries = Dict([

  :Gold => Dict([
    :Daily => Dict([
      :Morning => Dict([
                  :USD  => Dict([ Day(1) => "GOLDAMGBD228NLBM" ]),
                  :Euro => Dict([ Day(1) => "GOLDAMGBD230NLBM" ]),
      ]), # :Morning
      :Afternoon => Dict([
                  :USD  => Dict([ Day(1) => "GOLDPMGBD228NLBM" ]),
                  :Euro => Dict([ Day(1) => "GOLDPMGBD230NLBM" ]),
      ]), # :Afternoon
    ]), # :Daily
  ]), # :Gold

 :CrudeOil => Dict([
    :Daily => Dict([
      :WestTexas => Dict([
        :USD  => Dict([ Day(1) => "DCOILWTICOM" ]),
      ]),
      :Brent => Dict([
        :USD  => Dict([ Day(1) => "DCOILBRENTEU" ]),
      ]),
    ]), # :Daily
 ]), # :CrudeOil

 :NaturalGas => Dict([
    :Daily([
        :USD  => Dict([ Day(1) => "DHHNGSP" ]),
    ]), # :Daily
 ]), # :NaturalGas

]); # FredSpotSeries


text:
using MetadataTools
pkgmeta = get_all_pkg()
pg = make_dep_graph(pkgmeta)
pkg_graph = get_pkg_dep_graph("Gadfly", pg)


text:
import Test

using DelayedErrors

@noinline function grandparent1()::Nothing
    parent1()
    return nothing
end
@noinline function grandparent2()::Nothing
    parent2()
    return nothing
end
@noinline function grandparent3()::Nothing
    parent3()
    return nothing
end
@noinline function grandparent4()::Nothing
    parent4()
    return nothing
end

@noinline function parent1()::Nothing
    child1()
    return nothing
end
@noinline function parent2()::Nothing
    child2()
    return nothing
end
@noinline function parent3()::Nothing
    child3()
    return nothing
end
@noinline function parent4()::Nothing
    child4()
    return nothing
end

@noinline function child1()::Nothing
    grandchild1()
    return nothing
end
@noinline function child2()::Nothing
    grandchild2()
    return nothing
end
@noinline function child3()::Nothing
    grandchild3()
    return nothing
end
@noinline function child4()::Nothing
    grandchild4()
    return nothing
end

@noinline function grandchild1()::Nothing
    push_delayed_error(
        ;
        a = 1,
        b = "2",
        )
    return nothing
end
@noinline function grandchild2()::Nothing
    push_delayed_error(
        "This is a test";
        a = 1,
        b = "2",
        )
    return nothing
end
@noinline function grandchild3()::Nothing
    push_delayed_error(
        "This is ",
        "a test";
        a = 1,
        b = "2",
        )
    return nothing
end
@noinline function grandchild4()::Nothing
    push_delayed_error(
        "This ",
        "is ",
        "a ",
        "test";
        a = 1,
        b = "2",
        )
    return nothing
end

pop_delayed_errors()
pop_delayed_errors()
pop_delayed_errors()
grandparent1()
Test.@test_throws(
    ErrorException,
    pop_delayed_errors(),
    )
pop_delayed_errors()
pop_delayed_errors()
pop_delayed_errors()
grandparent2()
grandparent3()
grandparent4()
Test.@test_throws(
    ErrorException,
    pop_delayed_errors(),
    )
pop_delayed_errors()
pop_delayed_errors()
pop_delayed_errors()


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Test
using MIPLearn

MIPLearn.setup_logger()
const is_gurobi_available = ("GUROBI_HOME" in keys(ENV))

@testset "MIPLearn" begin
    include("fixtures/knapsack.jl")
    include("instance/file_instance_test.jl")
    include("instance/jump_instance_test.jl")
    include("solvers/jump_solver_test.jl")
    include("solvers/learning_solver_test.jl")
    # include("utils/benchmark_test.jl")
    include("utils/parse_test.jl")
    include("bb/lp_test.jl")
end


text:
# Unit Test for `SolveQuadraticProgram()`

using LinearAlgebra;
using SparseArrays;
using Gurobi; #<! Solves problems which SCS struggles with
using OSQP;
using Random;
using StableRNGs;
using Test;

seedNumber = 1234;
Random.seed!(seedNumber);
# Random.seed!(StableRNG(seedNumber), seedNumber);

include("GenerateQuadraticProgram.jl");
include("SolveQuadraticProgram.jl");
include("LinearSystemSolvers.jl");
include("SolveQuadraticProgramJump.jl");

# @enum ProblemClass randomQp = 1 inequalityConstrainedQp equalityConstrainedQp optimalControl portfolioOptimization lassoOptimization huberFitting supportVectorMachine isotonicRegression
@enum RefSolver gurobiSolver = 1 osqpSolver

## Parameters

# Simulaion

refSolver = osqpSolver;

numSimulations  = 10;
mNumElements = [    0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;];
mNumConstraints = [ 0000 0000;
0000 0000;
0005 0050;
0000 0000;
0000 0000;
0000 0000;
0000 0000;
0000 0000;
0000 0000;]; #<! Set to 0 for OSQP Paper dimensions

# Solver
numIterations   = 50000;
ϵAbs            = 1e-7;
ϵRel            = 1e-7;
ρ               = 0.1;
adptΡ           = true;
hLinSolInit     = FacLdlInit;
hLinSol         = FacLdl!;

absDevThr = 1e-5;

hOptFun = optimizer_with_attributes(Gurobi.Optimizer, "OptimalityTol" => 1e-8, "FeasibilityTol" => 1e-8);

@testset verbose = true "Unit Test vs. Gurobi" begin
for problemClass ∈ instances(ProblemClass)
    @testset "Unit Test of $(problemClass) Problem Class" begin
    for iSim ∈ 1:numSimulations
        for iDim ∈ 1:size(mNumElements, 2)
            mP, vQ, mA, vL, vU = GenerateRandomQP(problemClass, mNumElements[Int(problemClass), iDim]; numConstraints = mNumConstraints[Int(problemClass), iDim]);
            # Some problems might change the actual data
            numElements     = size(mP, 1);
            numConstraints  = size(mA, 1);

            vX = zeros(numElements);

            if (refSolver == gurobiSolver)
                runTimeRefSolver = @elapsed vT = SolveQpJump(mP, vQ, mA, vL, vU; hOptFun = hOptFun);
            elseif (refSolver == osqpSolver)
                runTimeRefSolver = @elapsed begin
                    osqPModel = OSQP.Model();
                    OSQP.setup!(osqPModel; P = mP, q = vQ, A = mA, l = vL, u = vU, rho = ρ, max_iter = numIterations, eps_abs = 1e-7, eps_rel = 1e-7, scaling = 0);
                    osqpRes = OSQP.solve!(osqPModel);
                    vT = osqpRes.x;
                end
            end

            runTime     = @elapsed convFlag = SolveQuadraticProgram!(vX, mP, vQ, mA, vL, vU, hLinSolInit, hLinSol; numIterations = numIterations, ϵAbs = ϵAbs, ϵRel = ϵRel, ρ = ρ, adptΡ = adptΡ);

            maxAbsDev = norm(vT - vX, Inf);

            println("The max absolute error is: $(maxAbsDev)");
            println("The run time is: $(runTime) [Sec]");
            println("The output flag is: $(convFlag)");
            println("The reference solver ($(refSolver)) run time is: $(runTimeRefSolver) [Sec]");
            @test maxAbsDev <= absDevThr;
        end
    end
end #<! End Test Set Problem Class
end

end #<! Test Set Unit Test

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
module Jin
using MacroTools

jind = Dict()
jind_original = Dict()
references = Dict()

function reset()
    empty!(jind)
    merge!(jind, jind_original)
end

function process_expr(ex)
    ex=deepcopy(ex)
    assign_blocks=[]
    kwargs = splitdef(ex)[:kwargs]
    f_name = splitdef(ex)[:name]
    f_nameq = Meta.quot(f_name)
    for kwarg in kwargs
        name, value = kwarg.args
        nameq = Meta.quot(name)
        qualified_name = :(($f_nameq, $nameq))
        assign_block = quote
            Jin.jind[$qualified_name] = $value
            Jin.jind_original[$qualified_name] = $value
        end
        push!(assign_blocks, assign_block)
        kwarg.args[2] = :(Jin.jind[$qualified_name])
    end
    assign_blocks, ex
end

macro register(f_ex)
    name = splitdef(f_ex)[:name]
    quote
        $(esc(f_ex))
        Jin.references[$(Meta.quot(name))] = $(esc(name))
    end
end

macro configurable(f_ex)
    assign_blocks, f_ex_new = process_expr(f_ex)
    name = splitdef(f_ex)[:name]
    quote
        $(assign_blocks...)
        $(esc(f_ex_new))
        Jin.references[$(Meta.quot(name))] = $(esc(name))
    end
end

function parse_value(value)
    value = strip(value)
    if isempty(value)
        throw("value can't be empty")
    end
    # This is a reference
    if value[1] == '@'
        ref_name = value[2:end] |> Symbol
        if !haskey(references, ref_name)
            throw("No reference to $ref_name found")
        end
        return references[ref_name]
    else
        return Meta.parse(value)
    end
end

function parse_file(text)
    lines = split(text, "\n")
    config_dict = Dict()
    for line in lines
        line = strip(line)
        if isempty(line) || startswith(line, "#")
            continue
        end
        line_parts = split(line, "=")
        if length(line_parts) != 2
            throw("Invalid line: $line")
        end
        variable, value = line_parts
        parsed_value = parse_value(value)
        variable_parts = split(variable, ".")
        if length(variable_parts) != 2
            throw("Invalid variable: $variable")
        end
        function_name, var_name = variable_parts
        config_dict[(Symbol(function_name |> strip), Symbol(var_name |> strip))]= parsed_value
    end
    return config_dict
end

function load_config(filename)
    contents = open(filename) |> read |> String
    dict = parse_file(contents)
    merge!(jind, dict)
end

end # module


text:
using Test
using SymbolicUtils

# == / != syntax is nice, let's use it in tests
macro eqtest(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end
SymbolicUtils.show_simplified[] = false

include("basics.jl")
include("order.jl")
include("rewrite.jl")
include("rulesets.jl")
include("nf.jl")
include("interface.jl")
include("fuzz.jl")

if haskey(ENV, "CI")
    include("benchmark.jl")
    using Pkg
    Pkg.add(url="https://github.com/SciML/ModelingToolkit.jl.git", rev="master")
    Pkg.test("ModelingToolkit")
end


text:
####################################################################
# ML estimation
####################################################################

model_g1 = Sem(
    specification = specification_g1,
    data = dat_g1,
    imply = RAMSymbolic
)

model_g2 = Sem(
    specification = specification_g2,
    data = dat_g2,
    imply = RAM
)

model_ml_multigroup = SemEnsemble(model_g1, model_g2; diff = semdiff)

############################################################################
### test gradients
############################################################################

@testset "ml_gradients_multigroup" begin
    @test test_gradient(model_ml_multigroup, start_test; atol = 1e-9)
end

# fit
@testset "ml_solution_multigroup" begin
    solution = sem_fit(model_ml_multigroup)
    update_estimate!(partable, solution)
    @test compare_estimates(
        partable,
        solution_lav[:parameter_estimates_ml]; atol = 1e-4,
        lav_groups = Dict(:Pasteur => 1, :Grant_White => 2))
end

@testset "fitmeasures/se_ml" begin
    solution_ml = sem_fit(model_ml_multigroup)
    @test all(test_fitmeasures(
        fit_measures(solution_ml),
        solution_lav[:fitmeasures_ml]; rtol = 1e-2, atol = 1e-7))

    update_partable!(
        partable, identifier(model_ml_multigroup), se_hessian(solution_ml), :se)
    @test compare_estimates(
        partable,
        solution_lav[:parameter_estimates_ml]; atol = 1e-3,
        col = :se, lav_col = :se,
        lav_groups = Dict(:Pasteur => 1, :Grant_White => 2))
end

####################################################################
# ML estimation - user defined loss function
####################################################################

struct UserSemML <: SemLossFunction
    objective
    gradient
    hessian
end

############################################################################
### constructor
############################################################################

UserSemML(;n_par, kwargs...) = UserSemML([1.0], zeros(n_par), zeros(n_par, n_par))

############################################################################
### functors
############################################################################

import LinearAlgebra: Symmetric, cholesky, isposdef, logdet, tr
import LinearAlgebra

function (semml::UserSemML)(par, F, G, H, model)
    if G error("analytic gradient of ML is not implemented (yet)") end
    if H error("analytic hessian of ML is not implemented (yet)") end

    a = cholesky(Symmetric(model.imply.Σ); check = false)
    if !isposdef(a)
        semml.objective[1] = Inf
    else
        ld = logdet(a)
        Σ_inv = LinearAlgebra.inv(a)
        if !isnothing(F)
            prod = Σ_inv*model.observed.obs_cov
            semml.objective[1] = ld + tr(prod)
        end
    end
end

# models
model_g1 = Sem(
    specification = specification_g1,
    data = dat_g1,
    imply = RAMSymbolic
)

model_g2 = SemFiniteDiff(
    specification = specification_g2,
    data = dat_g2,
    imply = RAMSymbolic,
    loss = UserSemML
)

model_ml_multigroup = SemEnsemble(model_g1, model_g2; diff = semdiff)

@testset "gradients_user_defined_loss" begin
    @test test_gradient(model_ml_multigroup, start_test; atol = 1e-9)
end

# fit
@testset "solution_user_defined_loss" begin
    solution = sem_fit(model_ml_multigroup)
    update_estimate!(partable, solution)
    @test compare_estimates(
        partable,
        solution_lav[:parameter_estimates_ml]; atol = 1e-4,
        lav_groups = Dict(:Pasteur => 1, :Grant_White => 2))
end

####################################################################
# GLS estimation
####################################################################

model_ls_g1 = Sem(
    specification = specification_g1,
    data = dat_g1,
    imply = RAMSymbolic,
    loss = SemWLS
)

model_ls_g2 = Sem(
    specification = specification_g2,
    data = dat_g2,
    imply = RAMSymbolic,
    loss = SemWLS
)

model_ls_multigroup = SemEnsemble(model_ls_g1, model_ls_g2; diff = semdiff)

@testset "ls_gradients_multigroup" begin
    @test test_gradient(model_ls_multigroup, start_test; atol = 1e-9)
end

@testset "ls_solution_multigroup" begin
    solution = sem_fit(model_ls_multigroup)
    update_estimate!(partable, solution)
    @test compare_estimates(
        partable,
        solution_lav[:parameter_estimates_ls]; atol = 1e-4,
        lav_groups = Dict(:Pasteur => 1, :Grant_White => 2))
end

@testset "fitmeasures/se_ls" begin
    solution_ls = sem_fit(model_ls_multigroup)
    @test all(test_fitmeasures(
        fit_measures(solution_ls),
        solution_lav[:fitmeasures_ls];
        fitmeasure_names = fitmeasure_names_ls, rtol = 1e-2, atol = 1e-5))

    update_partable!(
        partable, identifier(model_ls_multigroup), se_hessian(solution_ls), :se)
    @test compare_estimates(
        partable,
        solution_lav[:parameter_estimates_ls]; atol = 1e-2,
        col = :se, lav_col = :se,
        lav_groups = Dict(:Pasteur => 1, :Grant_White => 2))
end

if !isnothing(specification_miss_g1)

####################################################################
# FIML estimation
####################################################################

model_g1 = Sem(
    specification = specification_miss_g1,
    observed = SemObsMissing,
    loss = SemFIML,
    data = dat_miss_g1,
    imply = RAM,
    diff = SemDiffEmpty()
)

model_g2 = Sem(
    specification = specification_miss_g2,
    observed = SemObsMissing,
    loss = SemFIML,
    data = dat_miss_g2,
    imply = RAM,
    diff = SemDiffEmpty()
)

model_ml_multigroup = SemEnsemble(model_g1, model_g2; diff = semdiff)

############################################################################
### test gradients
############################################################################

start_test = [
    fill(0.5, 6);
    fill(1.0, 9);
    0.05; 0.01; 0.01; 0.05; 0.01; 0.05;
    fill(0.01, 9);
    fill(1.0, 9);
    0.05; 0.01; 0.01; 0.05; 0.01; 0.05;
    fill(0.01, 9)]

@testset "fiml_gradients_multigroup" begin
    @test test_gradient(model_ml_multigroup, start_test; atol = 1e-7)
end

@testset "fiml_solution_multigroup" begin
    solution = sem_fit(model_ml_multigroup)
    update_estimate!(partable_miss, solution)
    @test compare_estimates(
        partable_miss,
        solution_lav[:parameter_estimates_fiml]; atol = 1e-4,
        lav_groups = Dict(:Pasteur => 1, :Grant_White => 2))
end

@testset "fitmeasures/se_fiml" begin
    solution = sem_fit(model_ml_multigroup)
    @test all(test_fitmeasures(
        fit_measures(solution),
        solution_lav[:fitmeasures_fiml]; rtol = 1e-3, atol = 0))

    update_partable!(
        partable_miss, identifier(model_ml_multigroup), se_hessian(solution), :se)
    @test compare_estimates(
        partable_miss,
        solution_lav[:parameter_estimates_fiml]; atol = 1e-3,
        col = :se, lav_col = :se,
        lav_groups = Dict(:Pasteur => 1, :Grant_White => 2))
end

end