text:
#############################################################################################

print("   running RosCtx test...")

before_time = time()

@context RosCtx

function rosenbrock(x::Vector{Float64})
    a = 1.0
    b = 100.0
    result = 0.0
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end

x = rand(2)
@inferred(overdub(RosCtx(), rosenbrock, x))

messages = String[]
Cassette.prehook(::RosCtx, f, args...) = push!(messages, string("calling ", f, args))
@test overdub(RosCtx(), rosenbrock, x) == rosenbrock(x)
@test length(messages) > 90

Cassette.prehook(ctx::RosCtx, f, args...) = push!(ctx.metadata, string("calling ", f, args))
messages2 = String[]
@test overdub(RosCtx(metadata=messages2), rosenbrock, x) == rosenbrock(x)
@test messages == messages2

Cassette.prehook(::RosCtx, f, args...) = nothing
Cassette.prehook(ctx::RosCtx, f, args::Number...) = push!(ctx.metadata, args)
argslog = Any[]
@test @overdub(RosCtx(metadata=argslog), rosenbrock(x)) == rosenbrock(x)
for args in argslog
    @test all(x -> isa(x, Number), args)
end

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running HookCtx test...")

before_time = time()

@context HookCtx

Cassette.prehook(ctx::HookCtx, f, args...) = push!(ctx.metadata[1], (f, args))
Cassette.posthook(ctx::HookCtx, out, f, args...) = push!(ctx.metadata[2], (out, f, args))

pres, posts = Any[], Any[]
ctx = HookCtx(metadata=(pres, posts))
x1, x2, x3 = rand(Int), rand(Int), rand(Int)

@overdub(ctx, x1 + x2 * x3)
@test pres == [(*, (x2, x3)),
               (Base.mul_int, (x2, x3)),
               (+, (x1, x2*x3)),
               (Base.add_int, (x1, x2*x3))]
@test posts == [(Base.mul_int(x2, x3), Base.mul_int, (x2, x3)),
                (*(x2, x3), *, (x2, x3)),
                (Base.add_int(x1, x2*x3), Base.add_int, (x1, x2*x3)),
                (+(x1, x2*x3), +, (x1, x2*x3))]
empty!(pres)
empty!(posts)

Cassette.overdub(::HookCtx, ::typeof(+), args...) = +(args...)

@overdub(ctx, x1 + x2 * x3)
@test pres == [(*, (x2, x3)),
               (Base.mul_int, (x2, x3)),
               (+, (x1, x2*x3))]
@test posts == [(Base.mul_int(x2, x3), Base.mul_int, (x2, x3)),
                (*(x2, x3), *, (x2, x3)),
                (+(x1, x2*x3), +, (x1, x2*x3))]
empty!(pres)
empty!(posts)

@overdub(ctx, Core._apply(+, (x1, x2), (x2 * x3, x3)))
@test pres == [(tuple, (x1, x2)),
               (*, (x2, x3)),
               (Base.mul_int, (x2, x3)),
               (tuple, (x2*x3, x3)),
               (+, (x1, x2, x2*x3, x3))]
@test posts == [((x1, x2), tuple, (x1, x2)),
                (Base.mul_int(x2, x3), Base.mul_int, (x2, x3)),
                (*(x2, x3), *, (x2, x3)),
                ((x2*x3, x3), tuple, (x2*x3, x3)),
                (+(x1, x2, x2*x3, x3), +, (x1, x2, x2*x3, x3))]

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running SinCtx test...")

before_time = time()

@context SinCtx
x = rand()
sin_plus_cos(x) = sin(x) + cos(x)
@test @overdub(SinCtx(), sin_plus_cos(x)) === sin_plus_cos(x)
Cassette.overdub(::SinCtx, ::typeof(sin), x) = cos(x)
@test @overdub(SinCtx(), sin_plus_cos(x)) === (2 * cos(x))

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running FoldCtx test...")

before_time = time()

@context FoldCtx
x = 2
foldmul(x, args...) = Core._apply(Base.afoldl, (*, x), args...)
@test @overdub(FoldCtx(), foldmul(x)) === foldmul(x)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running CountCtx test...")

before_time = time()

@context CountCtx
count1 = Ref(0)
Cassette.prehook(ctx::CountCtx, f, args::Number...) = (ctx.metadata[] += 1)
@overdub(CountCtx(metadata=count1), sin(1))
Cassette.prehook(ctx::CountCtx, f, args::Number...) = (ctx.metadata[] += 2)
count2 = Ref(0)
@overdub(CountCtx(metadata=count2), sin(1))
@test (2 * count1[]) === count2[]

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running CountCtx2 test...")

before_time = time()

@context CountCtx2
mutable struct Count{T}
    count::Int
end
function Cassette.prehook(ctx::CountCtx2{Count{T}}, f, arg::T, args::T...) where {T}
    ctx.metadata.count += 1
end
mapstr(x) = map(string, x)
c = Count{Union{String,Int}}(0)
@test @overdub(CountCtx2(metadata=c), mapstr(1:10)) == mapstr(1:10)
@test c.count > 1000

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running SqrCtx test...")

before_time = time()

@context SqrCtx
square_closure(x) = (y -> y * x)(x)
x = rand()
@test square_closure(x) == @overdub(SqrCtx(), square_closure(x))

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running CompCtx test...")

before_time = time()

@context CompCtx
comprehension1(x) = [i for i in x]
comprehension2(f, x, y) = [f(x, i) for i in y]
f, x, y = hypot, rand(), rand(2)
@test comprehension1(x) == @overdub(CompCtx(), comprehension1(x))
@test comprehension1(y) == @overdub(CompCtx(), comprehension1(y))
@test comprehension2(f, x, y) == @overdub(CompCtx(), comprehension2(f, x, y))

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running PassCtx test...")

before_time = time()

@context PassCtx
sig_collection = DataType[]
mypass = @pass (ctx, ref) -> begin
    # TODO: check !in(sig, sig_collection) first to mimic a caching technique that will be used in practice
    push!(sig_collection, ref.signature)
    return ref.code_info
end
@overdub(PassCtx(pass=mypass), sum(rand(3)))
@test !isempty(sig_collection) && all(T -> T <: Tuple, sig_collection)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running PassFallbackCtx test...")

before_time = time()

@context PassFallbackCtx
fallbackpass = @pass (ctx, ref) -> begin
    if ref.signature <: Tuple{typeof(sin),Any}
        return :(cos($(Cassette.OVERDUB_ARGUMENTS_NAME)[2]))
    end
    return ref.code_info
end
x = rand(30)
sin_kernel(i) = i > 0.5 ? sin(i) : i
y = @inferred(overdub(PassFallbackCtx(pass=fallbackpass), sum, sin_kernel, x))
@test sum(i -> i > 0.5 ? cos(i) : i, x) === y

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running WorldCtx test...")

before_time = time()

@context WorldCtx
worldtest = 0
oldctx = WorldCtx()
overdub(oldctx, sin, 1)

Cassette.prehook(::WorldCtx, args...) = (global worldtest += 1)
overdub(WorldCtx(), sin, 1)
@test worldtest > 100

tmp = worldtest
overdub(oldctx, sin, 1)
@test tmp < worldtest

tmp = worldtest
Cassette.prehook(::WorldCtx, args...) = nothing
overdub(WorldCtx(), sin, 1)
@test tmp === worldtest

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running TraceCtx test...")

before_time = time()

x, y, z = rand(3)
trtest(x, y, z) = x*y + y*z
trkwtest(x; _y = 1.0, _z = 2.0) = trtest(x, _y, _z)

@context TraceCtx

function Cassette.overdub(ctx::TraceCtx, args...)
    subtrace = Any[]
    push!(ctx.metadata, args => subtrace)
    if canrecurse(ctx, args...)
        newctx = similarcontext(ctx, metadata = subtrace)
        return recurse(newctx, args...)
    else
        return fallback(ctx, args...)
    end
end
trace = Any[]
@test @overdub(TraceCtx(metadata = trace), trtest(x, y, z)) == trtest(x, y, z)
@test trace == Any[
    (trtest,x,y,z) => Any[
        (*,x,y) => Any[(Base.mul_float,x,y)=>Any[]]
        (*,y,z) => Any[(Base.mul_float,y,z)=>Any[]]
        (+,x*y,y*z) => Any[(Base.add_float,x*y,y*z)=>Any[]]
    ]
]

# jrevels/Cassette.jl#48
@context HookTraceCtx

mutable struct HookTrace
    current::Vector{Any}
    stack::Vector{Any}
    HookTrace() = new(Any[], Any[])
end

tracekw = Any[]
@overdub(TraceCtx(metadata = tracekw), trkwtest(x, _y = y, _z = z)) == trtest(x, y, z)
subtracekw = first(Iterators.filter(t -> t[1] === (Core.kwfunc(trkwtest), (_y = y, _z = z), trkwtest, x), tracekw))[2]
@test subtracekw == trace

function enter!(t::HookTrace, args...)
    pair = args => Any[]
    push!(t.current, pair)
    push!(t.stack, t.current)
    t.current = pair.second
    return nothing
end

function exit!(t::HookTrace)
    t.current = pop!(t.stack)
    return nothing
end

Cassette.prehook(ctx::HookTraceCtx, args...) = enter!(ctx.metadata, args...)
Cassette.posthook(ctx::HookTraceCtx, args...) = exit!(ctx.metadata)

htrace = HookTrace()
@overdub(HookTraceCtx(metadata = htrace), trtest(x, y, z))
@test htrace.current == trace

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running NestedReflectCtx test...")

before_time = time()

@context NestedReflectCtx
r_pre = Cassette.reflect((typeof(sin), Int))
r_post = Cassette.reflect((typeof(overdub), typeof(NestedReflectCtx()), typeof(sin), Int))
@test isa(r_pre, Cassette.Reflection) && isa(r_post, Cassette.Reflection)
Cassette.overdub_pass!(r_pre, typeof(NestedReflectCtx()))
@test r_pre.code_info.code == r_post.code_info.code

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running GemvCtx test...")

before_time = time()

# issue #51
@context GemvCtx
α, β = 1.0, 2.0
A, X, Y = rand(Float64, 10, 10), rand(Float64, 10), rand(Float64, 10)
Y_copy = copy(Y)
Y_out = overdub(GemvCtx(), LinearAlgebra.BLAS.gemv!, 'T', α, A, X, β, Y)
Y_copy_out = LinearAlgebra.BLAS.gemv!('T', α, A, X, β, Y_copy)
@test Y_out == Y_copy_out
@test Y_out === Y
@test Y_copy_out === Y_copy

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running InferCtx test...")

before_time = time()

@context InferCtx

dispatchtupletest(::Type{T}) where {T} = Base.isdispatchtuple(Tuple{T}) ? T : Any
relu(x) = max(zero(x), x)
relulayer(W, x, b) = relu.(W*x .+ b)
kwargtest(foobar; foo = 1, bar = 2) = nothing

@inferred(overdub(InferCtx(), typejoin, Float32, Float32, Float32))
@inferred(overdub(InferCtx(), dispatchtupletest, Float32))
@inferred(overdub(InferCtx(), (a, b) -> Core.apply_type(a, b), AbstractVector, Int))
@inferred(overdub(InferCtx(), eltype, rand(1)))
@inferred(overdub(InferCtx(), *, rand(1, 1), rand(1, 1)))
@inferred(overdub(InferCtx(), *, rand(Float32, 1, 1), rand(Float32, 1, 1)))
if VERSION <= v"1.3"
    @inferred(overdub(InferCtx(), *, rand(Float32, 1, 1), rand(Float32, 1)))
    @inferred(overdub(InferCtx(), relulayer, rand(Float64, 1, 1), rand(Float32, 1), rand(Float32, 1)))
    @inferred(overdub(InferCtx(), rand, Float32, 1))
else
    # test depends on constant propagation
    @test_throws Exception @inferred(overdub(InferCtx(), *, rand(Float32, 1, 1), rand(Float32, 1)))
    # test depends on M*v which is the test above
    @test_throws Exception @inferred(overdub(InferCtx(), relulayer, rand(Float64, 1, 1), rand(Float32, 1), rand(Float32, 1)))
    # XXX: Figure out why this broke
    @test_throws Exception @inferred(overdub(InferCtx(), rand, Float32, 1))
end
@inferred(overdub(InferCtx(), broadcast, +, rand(1), rand(1)))
@inferred(overdub(InferCtx(), () -> kwargtest(42; foo = 1, bar = 2)))

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running InModuleCtx test...")

before_time = time()

module DefineStuffInModule
    using Cassette
    Cassette.@context InModuleCtx
    x = 0
    incrpass = Cassette.@pass (ctx, ref) -> (global x += 1; ref.code_info)
    f(x) = Cassette.overdub(InModuleCtx(pass = incrpass), sin, x)
end

@test DefineStuffInModule.f(1) === sin(1)
@test DefineStuffInModule.x > 0

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running InvokeCtx test...")

before_time = time()

invoked(x::Int) = x + x
invoked(x::Number) = x^2
invoker(x) = invoke(invoked, Tuple{Number}, x)

@context InvokeCtx
Cassette.prehook(ctx::InvokeCtx, f, args...) = push!(ctx.metadata, f)
ctx = InvokeCtx(metadata=Any[])

@test overdub(ctx, invoker, 3) === 9
# This is kind of fragile and may break for unrelated reasons - the main thing
# we're testing here is that we properly trace through the `invoke` call.
@test ctx.metadata == Any[Core.apply_type, Core.invoke, Core.apply_type,
                          Val{2}, Core.apply_type, Base.literal_pow, *,
                          Base.mul_int]

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running Invoke2Ctx test...")

before_time = time()

# issue #89
invoked2(x::Int) = x + x
invoked2(x::Number, y=4) = x^2 + y
invoker2(x) = invoke(invoked2, Tuple{Number}, x)

@context Invoke2Ctx
Cassette.prehook(ctx::Invoke2Ctx, f, args...) = push!(ctx.metadata, f)
ctx = Invoke2Ctx(metadata=Any[])

@test overdub(ctx, invoker2, 3) === invoker2(3) === 13

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

# taken from https://stackoverflow.com/questions/52050262/how-to-do-memoization-or-memoisation-in-julia-1-0/52062639#52062639

print("   running MemoizeCtx test...")

before_time = time()

fib(x) = x < 3 ? 1 : fib(x - 2) + fib(x - 1)
fibtest(n) = fib(2 * n) + n

@context MemoizeCtx

function Cassette.overdub(ctx::MemoizeCtx, ::typeof(fib), x)
    result = get(ctx.metadata, x, 0)
    result === 0 && return recurse(ctx, fib, x)
    return result
end

Cassette.posthook(ctx::MemoizeCtx, fibx, ::typeof(fib), x) = (ctx.metadata[x] = fibx)

ctx = MemoizeCtx(metadata = Dict{Int,Int}())
n = 10
result = Cassette.overdub(ctx, fibtest, n)

@test result == fibtest(n)
@test length(ctx.metadata) == 2 * n
@test all(fib(k) == v for (k, v) in ctx.metadata)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

# ref https://github.com/jrevels/Cassette.jl/issues/73

if VERSION >= v"1.1-"
    print("   running NoOpCtx test...")

    before_time = time()

    @context NoOpCtx

    function loop73(x, n)
        r = x / x
        while n > 0
            r *= sin(x)
            n -= 1
        end
        return r
    end

    f73(x, n) = overdub(NoOpCtx(), loop73, x, n)
    ff73(x, n) = overdub(NoOpCtx(), f73, x, n)
    fff73(x, n) = overdub(NoOpCtx(), ff73, x, n)

    f73(2, 50) # warm up
    ff73(2, 50) # warm up
    fff73(2, 50) # warm up

    @test @allocated(f73(2, 50)) == 0
    @test @allocated(ff73(2, 50)) == 0
    @test_broken @allocated(fff73(2, 50)) == 0

    println("done (took ", time() - before_time, " seconds)")
end

#############################################################################################

print("   running DisableHooksCtx test...")

before_time = time()

@context DisableHooksCtx
Cassette.prehook(ctx::DisableHooksCtx, args...) = ctx.metadata[1] += 1
Cassette.posthook(ctx::DisableHooksCtx, args...) = ctx.metadata[2] += 1
ctx = disablehooks(DisableHooksCtx(metadata = [0, 0]))
@test overdub(ctx, sin, 1.0) == sin(1.0)
@test all(ctx.metadata .== 0)
@test overdub(ctx, overdub, ctx, sin, 1.0) == sin(1.0)
@test all(ctx.metadata .== 0)

println("done (took ", time() - before_time, " seconds)")

#############################################################################################

print("   running SliceCtx test...")

before_time = time()

using Core: CodeInfo, SlotNumber, SSAValue

@context SliceCtx

function Cassette.overdub(ctx::SliceCtx, f, callback, args...)
    if canrecurse(ctx, f, args...)
        _ctx = similarcontext(ctx, metadata = callback)
        return recurse(_ctx, f, args...) # return result, callback
    else
        return fallback(ctx, f, args...), callback
    end
end

const global_test_cache = Any[]

push_to_global_test_cache!(x) = push!(global_test_cache, x)

function Cassette.overdub(ctx::SliceCtx, ::typeof(push_to_global_test_cache!), callback, x)
    return nothing, () -> (callback(); push_to_global_test_cache!(x))
end

function sliceprintln(::Type{<:SliceCtx}, reflection::Cassette.Reflection)
    ir = reflection.code_info
    callbackslotname = gensym("callback")
    push!(ir.slotnames, callbackslotname)
    push!(ir.slotflags, 0x00)
    callbackslot = SlotNumber(length(ir.slotnames))
    getmetadata = Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), Expr(:contextslot), QuoteNode(:metadata))

    # insert the initial `callbackslot` assignment into the IR.
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                (stmt, i) -> i == 1 ? 2 : nothing,
                                (stmt, i) -> [Expr(:(=), callbackslot, getmetadata), stmt])

    # replace all calls of the form `f(args...)` with `f(callback, args...)`, taking care to
    # properly handle Core._apply calls and destructure the returned `(result, callback)`
    # into the appropriate statements
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                (stmt, i) -> begin
                                    i > 1 || return nothing # don't slice the callback assignment
                                    stmt = Base.Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
                                    if Base.Meta.isexpr(stmt, :call)
                                        isapply = Cassette.is_ir_element(stmt.args[1], GlobalRef(Core, :_apply), ir.code)
                                        return 3 + isapply
                                    end
                                    return nothing
                                end,
                                (stmt, i) -> begin
                                    items = Any[]
                                    callstmt = Base.Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
                                    callssa = SSAValue(i)
                                    if Cassette.is_ir_element(callstmt.args[1], GlobalRef(Core, :_apply), ir.code)
                                        push!(items, Expr(:call, Expr(:nooverdub, GlobalRef(Core, :tuple)), callbackslot))
                                        push!(items, Expr(:call, callstmt.args[1], callstmt.args[2], SSAValue(i), callstmt.args[3:end]...))
                                        callssa = SSAValue(i + 1)
                                    else
                                        push!(items, Expr(:call, callstmt.args[1], callbackslot, callstmt.args[2:end]...))
                                    end
                                    push!(items, Expr(:(=), callbackslot, Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), callssa, 2)))
                                    result = Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), callssa, 1)
                                    if Base.Meta.isexpr(stmt, :(=))
                                        result = Expr(:(=), stmt.args[1], result)
                                    end
                                    push!(items, result)
                                    return items
                                end)

    # replace return statements of the form `return x` with `return (x, callback)`
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                  (stmt, i) -> Base.Meta.isexpr(stmt, :return) ? 2 : nothing,
                                  (stmt, i) -> begin
                                      return [
                                          Expr(:call, Expr(:nooverdub, GlobalRef(Core, :tuple)), stmt.args[1], callbackslot)
                                          Expr(:return, SSAValue(i))
                                      ]
                                  end)
    return ir
end

const sliceprintlnpass = Cassette.@pass sliceprintln

a = rand(3)
b = rand(3)

function add(a, b)
    push_to_global_test_cache!(a)
    push_to_global_test_cache!(b)
    c = a + b
    push_to_global_test_cache!(c)
    return c
end

ctx = SliceCtx(pass=sliceprintlnpass, metadata = () -> nothing)

result, callback = Cassette.recurse(ctx, add, a, b)

@test result == a + b
@test isempty(global_test_cache)
callback()
@test global_test_cache == [a, b, result]

println("done (took ", time() - before_time, " seconds)")

# Test overdubbing of a call overload invoke

using LinearAlgebra

struct Dense{F,S,T}
  W::S
  b::T
  σ::F
end

Dense(W, b) = Dense(W, b, identity)

function (a::Dense)(x::AbstractArray)
  W, b, σ = a.W, a.b, a.σ
  σ.(W*x .+ b)
end

function (a::Dense{<:Any,W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}}
  invoke(a, Tuple{AbstractArray}, x)
end

function (a::Dense{<:Any,W})(x::AbstractArray{<:Real}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}}
  a(T.(x))
end

Cassette.@context CtxCallOverload;
let d = Dense(3,3)
    data = rand(3)
    Cassette.overdub(CtxCallOverload(), d, data)
end

#############################################################################################

print("   running OverdubOverdubCtx test...")

# Fixed in PR #148
Cassette.@context OverdubOverdubCtx
overdub_overdub_me() = 2
Cassette.overdub(OverdubOverdubCtx(), Cassette.overdub, OverdubOverdubCtx(), overdub_overdub_me)

#############################################################################################

print("   running NukeCtx test...")

@Cassette.context NukeContext
struct Silo; end

Base.iterate(x::Silo) = (println("Launching Nukes"); error("What's the point?"))

function Cassette.overdub(ctx::NukeContext, ::typeof(iterate), x::Silo)
    nothing
end

@test Cassette.overdub(NukeContext(), iterate, Silo()) === nothing

launch(s::Silo) = (s...,)

if VERSION >= v"1.4.0-DEV.304"
    @test Cassette.overdub(NukeContext(), launch, Silo()) === ()
else
    @test_broken Cassette.overdub(NukeContext(), launch, Silo()) === ()
end

if VERSION >= v"1.4.0-DEV.304"
    Cassette.@context ApplyIterateCtx;

    const instructions = []
    function Cassette.prehook(ctx::ApplyIterateCtx,
                              op::Any,
                              a::T1, b::T2) where {T1, T2}
        push!(instructions, (op, T1, T2))
    end

    Cassette.overdub(ApplyIterateCtx(), ()->pi*2.0)
    @test instructions[end] === (Core.Intrinsics.mul_float, Float64, Float64)
end


text:
using HttpServer, DataFrames

function loadScoredFile(dfd::DataFrame)
    println("loadScoredFile")
    return true
end

function rest(ip::String)
    http = HttpHandler() do req::Request, res::Response
        if ismatch(r"^/optimize/", req.resource)
            try
                dfd = readtable(IOBuffer(convert(String,req.data) ) )
                println(dfd)
                if loadScoredFile(dfd)
                    return Response("done")
                else
                    return Response("error - could not upload file")
                end
            catch
                return Response(422, "Error: coudn't solve puzzle.")
            end
        elseif ismatch(r"^/tst/", req.resource)
            return Response("ooo hello!")
        else
            return Response("unrecognized URL")
        end
    end
    server = Server( http )
    run(server, host=IPv4(ip), port=8000)
end

rest("10.63.36.18")

#wget 10.63.36.18:8000/hello/name/  --post-data "email=abc@abc.com&file1=@FILE_HERE&file2=@FILE_HERE"
#wget 10.63.36.18:8000/hello/name/  --post-file=file
# wget 10.63.36.18:8000/optimize/  --post-file=/mnt/resource/analytics/models/rev/scored.csv

ext: jl

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
## Incidence matrix
ei,ej = copy.(eachcol(Int.(readdlm("undir-graph-14.edges"))))
xy = readdlm("undir-graph-14.xy")'
n = max(maximum(ei),maximum(ej))
A = sparse(ei,ej,1,n,n)
A = A+A'
##
m = length(ei) # number of edges
B = sparse(1:m,ei,1,m,n) - sparse(1:m,ej,1,m,n)
## Output this for Latex
using Printf
@printf("\\left[ \\begin{array}{*{14}{c}} \n ")
for i=1:m
  for j=1:n
    if B[i,j] == 0
      @printf("\\textcolor{gray}{0} & ")
    else
      @printf("\\mathbf{%i} & ", B[i,j])
    end
  end
  println("\\\\")
end
@printf("\\end{array} \\right] \n ")


text:
mutable struct PiCamera
    cam::PyObject
    framerate::Int
    resolution::Tuple{Int, Int}
    resize::Tuple{Int, Int}
    function PiCamera(framerate::Int, W::Int, H::Int, r::Int)
        cam = py"PC"()
        resolution = (16W, 16H)
        cam.resolution = resolution
        cam.framerate = framerate
        resize = resolution .÷ r
        new(cam, framerate, resolution, resize)
    end
end

Base.isopen(a::PiCamera) = !a.cam.closed
Base.close(a::PiCamera) = isopen(a) && a.cam.close()
function Base.open(a::PiCamera)
    close(a)
    cam = py"PC"()
    cam.resolution = a.resolution
    cam.framerate = a.framerate
    a.cam = cam
end
restart(a::PiCamera) = (close(a); open(a))

function snap(camera::PiCamera)
    stream = py"PIO"()
    camera.cam.capture(stream, splitter_port = 0, resize = camera.resize, format = "jpeg", use_video_port = true)
    stream.seek(0)
    buffer = IOBuffer(stream.read())
    load(buffer)
end


text:
@testset "5.4 Inverse cotangent" begin
    include("5.4.1 Inverse cotangent functions.jl")
    include("5.4.2 Exponentials of inverse cotangent.jl")
end

max_stars_repo_path: test/5 Inverse trig functions/5.4 Inverse cotangent/.jl


text:
@testset "Fourier Rotate" begin

    @testset "Compare with ImageTransformations" begin

        function f(θ)
            x = 1.0 .* range(0.0, 1.0, length=256)' .* range(0.0, 1.0, length=256)
            f(x) = sin(x * 20) + tan(1.2 * x) + sin(x) + cos(1.1323 * x) * x^3 + x^3 + 0.23 * x^4 + sin(1/(x+0.1))
            img = 5 .+ abs.(f.(x))
            img ./= maximum(img)
            img[20:40, 100:200] .= 1
            img[20:200, 20:90] .= 0.3
            img[20:200, 100:102] .= 0.7

            m = sum(img) / length(img)

            img_1 = parent(ImageTransformations.imrotate(img, deg2rad(θ), m))
            z = ones(Float32, size(img_1))
            z .*= m
            FourierTools.center_set!(z, img)
            img_2 = FourierTools.rotate(z, θ)
            img_3 = real(FourierTools.rotate(z .+ 0im, θ))
            img_4 = FourierTools.rotate!(z, θ)

            @test all(.≈(img_1, img_2, rtol=0.6))
            @test ≈(img_1, img_2, rtol=0.03)
            @test ≈(img_3, img_2, rtol=0.01)
            @test ==(img_4, z)

            img_1c = FourierTools.center_extract(img_1, (100, 100))
            img_2c = FourierTools.center_extract(img_2, (100, 100))
            @test all(.≈(img_1c, img_2c, rtol=0.3))
            @test ≈(img_1c, img_2c, rtol=0.05)
        end

        f(-54.31)
        f(-32.31)
        f(32.31)
        f(0)
    end
end


text:
"""
    AbstractStepRange

Supertype for `StepSRange` and `StepMRange`. It's subtypes should behave
identically to `StepRange`.
"""
abstract type AbstractStepRange{T,S} <: OrdinalRange{T,S} end

function Base.isempty(r::AbstractStepRange)
    (first(r) != last(r)) & ((step(r) > zero(step(r))) != (last(r) > first(r)))
end

"""
    StepSRange

A static range with elements of type `T` with spacing of type `S`. The step
between each element is constant, and the range is defined in terms of a
`start` and `stop` of type `T` and a `step` of type `S`. Neither `T` nor `S`
should be floating point types.
"""
struct StepSRange{T,Ts,F,S,L} <: AbstractStepRange{T,Ts}

    function StepSRange{T,Ts}(start::T, step::Ts, stop::T) where {T,Ts}
        return new{T,Ts,start,step,Base.steprange_last(start, step, stop)}()
    end
end

function Base.getproperty(r::StepSRange, s::Symbol)
    if s === :start
        return first(r)
    elseif s === :step
        return step(r)
    elseif s === :stop
        return last(r)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

Base.first(r::StepSRange{T,Ts,F,S,L}) where {T,Ts,F,S,L} = F

Base.last(r::StepSRange{T,Ts,F,S,L}) where {T,Ts,F,S,L} = L

Base.length(r::StepSRange) = StaticArrays.get(Length(r))

function (::Type{<:StepSRange{T1,T2} where T1})(r::AbstractRange) where {T2}
    return StepSRange{eltype(r),T2}(r)
end

"""
    StepMRange

A mutable range with elements of type `T` with spacing of type `S`. The step
between each element is constant, and the range is defined in terms of a
`start` and `stop` of type `T` and a `step` of type `S`. Neither `T` nor `S`
should be floating point types.
"""
mutable struct StepMRange{T,S} <: AbstractStepRange{T,S}
    start::T
    step::S
    stop::T

    function StepMRange{T,S}(start::T, step::S, stop::T) where {T,S}
        return new(start, step, Base.steprange_last(start,step,stop))
    end
end

function (::Type{StepMRange{T1,T2} where T1})(r::AbstractRange) where {T2}
    return StepMRange{eltype(r),T2}(r)
end

Base.first(r::StepMRange) = getfield(r, :start)

Base.last(r::StepMRange) = getfield(r, :stop)

function Base.length(r::StepMRange{T}) where {T}
    return start_step_stop_to_length(T, first(r), step(r), last(r))
end

function Base.setproperty!(r::StepMRange, s::Symbol, val)
    if s === :start
        return set_first!(r, val)
    elseif s === :step
        return set_step!(r, val)
    elseif s === :stop
        return set_last!(r, val)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

for (F,f) in ((:M,:m), (:S,:s))
    SR = Symbol(:Step, F, :Range)
    frange = Symbol(f, :range)
    @eval begin
        function Base.getindex(r::$(SR), s::AbstractRange{<:Integer})
            Base.@_inline_meta
            @boundscheck checkbounds(r, s)
            st = oftype(first(r), first(r) + (first(s)-1)*step(r))
            return $(frange)(st, step=step(r)*step(s), length=length(s))
        end

        $(SR)(r::AbstractUnitRange{T}) where {T} = $(SR){T,T}(first(r), step(r), last(r))

        $(SR)(start::T, step::S, stop::T) where {T,S} = $(SR){T,S}(start, step, stop)

        $(SR){T1,T2}(r::$(SR){T1,T2}) where {T1,T2} = r
        function $(SR){T1,T2}(r::AbstractRange) where {T1,T2}
            return $(SR){T1,T2}(
                convert(T1, first(r)),
                convert(T2, step(r)),
                convert(T1, last(r))
               )
        end

        function Base.:(-)(r::$(SR))
            return $(frange)(-first(r), step=-step(r), length=length(r))
        end
    end
end


text:
module SomaticEvolution

using Distributions
using Statistics
using Random
using StatsBase
using DataFrames
using RecipesBase
using GLM
using LaTeXStrings
using Printf
using CSV
using Plots
using JSON
using DelimitedFiles
using AbstractTrees

export Cell,
MultiSimulation,
SimulationTracker,
SimulationInput,
MultilevelInput,
BranchingMoranInput,
BranchingInput,
MoranInput,
MultilevelInput,
Simulation,
SampledData,
PlotVAF,
PlotInverseVAF,
VAFResult,

multilevel_simulation,
multilevel_simulation_fast,
multiplesimulations,
run1simulation,
getallelefreq,
sampledallelefreq,
fitinverse,
gethist,
plotvaf,
plotvaf!,
plotinversevaf,
plotinversevaf!,
getstats,
saveinput,
branchingprocess,
moranprocess,
get_simulation,
pairwise_fixed_differences,
pairwise_fixed_differences_matrix,
pairwise_fixed_differences_statistics,
shared_fixed_mutations,
newmoduletimes,
cellpopulationsize,
average_mutations,
average_mutations_per_module,
mutations_per_cell,
mutation_ids_by_cell,
clonal_mutation_ids,
clonal_mutations,
subclonefreq,
getVAFresult,
age,
saveinput,
loadinput,
run_multilevel_from_file

include("definetypes.jl")
include("simulations.jl")
include("multilevel.jl")
include("process.jl")
include("sampling.jl")
include("analyse.jl")
include("multirun.jl")
include("plot.jl")
include("util.jl")
include("statistics.jl")
include("io.jl")
include("trees.jl")

end


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, Langevin, LinearAlgebra, Plots
using KitBase.ProgressMeter: @showprogress

cd(@__DIR__)

begin
    x0 = -1
    x1 = 1
    ncell = 40
    nface = ncell + 1
    dx = (x1 - x0) / ncell
    deg = 2 # polynomial degree
    nsp = deg + 1
    γ = 5 / 3
    cfl = 0.03
    dt = cfl * dx
    t = 0.0
end

ps = FRPSpace1D(x0, x1, ncell, deg)

begin
    uqMethod = "galerkin"
    nr = 5
    nRec = 10
    opType = "uniform"
    parameter1 = 0.9
    parameter2 = 1.1
end

uq = UQ1D(nr, nRec, parameter1, parameter2, opType, uqMethod)

u = zeros(ncell, nsp, uq.nm+1)
for i = 1:ncell, j = 1:nsp
    #u[i, j, 1] = 1.0
    #u[i, j, :] .+= uq.pce * 0.1 * sin(2π * ps.xpg[i, j])
    u[i, j, :] .= uq.pce * sin(π * ps.xpg[i, j])
end

a = zeros(uq.nm+1)
a[1] = 1.0

function dudt!(du, u, p, t)
    du .= 0.0
    J, ll, lr, lpdm, dgl, dgr, Δx, γ, uq, a = p

    nm = uq.nm
    nq = uq.nq

    ncell = size(u, 1)
    nsp = size(u, 2)

    u_ran = zeros(ncell, nsp, nq)
    for i = 1:ncell, j = 1:nsp
        u_ran[i, j, :] .= chaos_ran(u[i, j, :], uq)
    end
    a_ran = chaos_ran(a, uq)

    f = zeros(ncell, nsp, nm+1)
    for i = 1:ncell, j = 1:nsp
        _f = zeros(nq)
        for k = 1:nq
            _f[k] = u_ran[i, j, k] * a_ran[k] / J[i]
        end

        f[i, j, :] .= ran_chaos(_f, uq)
    end

    u_face = zeros(ncell, nm+1, 2)
    f_face = zeros(ncell, nm+1, 2)
    for i = 1:ncell, k = 1:nm+1
        u_face[i, k, 1] = dot(u[i, :, k], lr)
        f_face[i, k, 1] = dot(f[i, :, k], lr)
        u_face[i, k, 2] = dot(u[i, :, k], ll)
        f_face[i, k, 2] = dot(f[i, :, k], ll)
    end

    f_interaction = zeros(ncell + 1, nm+1)
    for i = 2:ncell
        f_interaction[i, :] .= (f_face[i-1, :, 1] .+ f_face[i, :, 2]) ./ 2 -
            (Δx[i-1] + Δx[i]) / 2 * (u_face[i, :, 2] .- u_face[i-1, :, 1])
    end
    f_interaction[1, :] .= (f_face[ncell, :, 1] .+ f_face[1, :, 2]) ./ 2 -
        (Δx[ncell] + Δx[1]) / 2 * (u_face[1, :, 2] .- u_face[ncell, :, 1])
    f_interaction[ncell+1, :] .= (f_face[ncell, :, 1] .+ f_face[1, :, 2]) ./ 2 -
        (Δx[ncell] + Δx[1]) / 2 * (u_face[1, :, 2] .- u_face[ncell, :, 1])

    rhs1 = zero(u)
    for i = 1:ncell, ppp1 = 1:nsp, l = 1:nm+1
        rhs1[i, ppp1, l] = dot(f[i, :, l], lpdm[ppp1, :])
    end

    idx = 1:ncell
    for i in idx, ppp1 = 1:nsp, l = 1:nm+1
        du[i, ppp1, l] =
            -(
                rhs1[i, ppp1, l] +
                (f_interaction[i, l] - f_face[i, l, 2]) * dgl[ppp1] +
                (f_interaction[i+1, l] - f_face[i, l, 1]) * dgr[ppp1]
            )
    end
end

tspan = (0.0, 50.0)
p = (ps.J, ps.ll, ps.lr, ps.dl, ps.dhl, ps.dhr, ps.dx, γ, uq, a)
prob = ODEProblem(dudt!, u, tspan, p)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

begin
    x = zeros(ncell * nsp)
    w = zeros(ncell * nsp, uq.nm+1)
    for i = 1:ncell
        idx0 = (i - 1) * nsp

        for j = 1:nsp
            idx = idx0 + j
            x[idx] = ps.xpg[i, j]

            w[idx, :] .= itg.u[i, j, :]
        end
    end

    sol = zeros(ncell*nsp, 2)
    for i in axes(sol, 1)
        sol[i, 1] = mean(w[i, :], uq.op)
        sol[i, 2] = std(w[i, :], uq.op)
    end

    pic1 = plot(x, sol[:, 1], label="u", xlabel="x", ylabel="mean")
    pic2 = plot(x, sol[:, 2], label="u", xlabel="x", ylabel="std")
    plot(pic1, pic2)
end
sol0 = deepcopy(sol)

@showprogress for iter = 1:nt
    step!(itg)
end

begin
    x = zeros(ncell * nsp)
    w = zeros(ncell * nsp, uq.nm+1)
    for i = 1:ncell
        idx0 = (i - 1) * nsp

        for j = 1:nsp
            idx = idx0 + j
            x[idx] = ps.xpg[i, j]

            w[idx, :] .= itg.u[i, j, :]
        end
    end

    sol = zeros(ncell*nsp, 2)
    for i in axes(sol, 1)
        sol[i, 1] = mean(w[i, :], uq.op)
        sol[i, 2] = std(w[i, :], uq.op)
    end

    pic1 = plot(x, sol[:, 1], label="u", xlabel="x", ylabel="mean")
    pic2 = plot(x, sol[:, 2], label="u", xlabel="x", ylabel="std")
    plot(pic1, pic2)
end

plot(x, sol[:, 1], label="Numerical", lw=1.5, xlabel="x", ylabel="expected value")
plot!(x, sol0[:, 1], label="Exact", lw=1.5, line=:dash)
savefig("wave_mean.pdf")

plot(x, sol[:, 2], label="Numerical", lw=1.5, xlabel="x", ylabel="standard deviation")
plot!(x, sol0[:, 2], label="Exact", lw=1.5, line=:dash)
savefig("wave_std.pdf")

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
# Extends PowerModels/src/form/bf.jl

## Variables

""
function variable_ne_branch_current(pm::_PM.AbstractBFModel; kwargs...)
    variable_ne_buspair_current_magnitude_sqr(pm; kwargs...)
end

""
function variable_ne_buspair_current_magnitude_sqr(pm::_PM.AbstractBFAModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true)
end

## Constraint templates

"""
This constraint captures problem agnostic constraints that are used to link
the model's current variables together, in addition to the standard problem
formulation constraints.  The network expansion name (ne) indicates that the
currents in this constraint can be set to zero via an indicator variable.
"""
function constraint_ne_model_current(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default)
    constraint_ne_model_current(pm, nw)
end

""
function constraint_ne_power_losses(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    branch = _PM.ref(pm, nw, :ne_branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    r = branch["br_r"]
    x = branch["br_x"]
    tm = branch["tap"]
    g_sh_fr = branch["g_fr"]
    g_sh_to = branch["g_to"]
    b_sh_fr = branch["b_fr"]
    b_sh_to = branch["b_to"]

    vad_min = _PM.ref(pm, nw, :off_angmin)
    vad_max = _PM.ref(pm, nw, :off_angmax)

    constraint_ne_power_losses(pm, nw, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, g_sh_to, b_sh_fr, b_sh_to, tm, vad_min, vad_max)
end

"""
Defines voltage drop over a branch, linking from and to side voltage magnitude
"""
function constraint_ne_voltage_magnitude_difference(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    branch = _PM.ref(pm, nw, :ne_branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    r = branch["br_r"]
    x = branch["br_x"]
    g_sh_fr = branch["g_fr"]
    b_sh_fr = branch["b_fr"]
    tm = branch["tap"]

    constraint_ne_voltage_magnitude_difference(pm, nw, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, b_sh_fr, tm)
end

## Actual constraints

"""
do nothing, most models do not require any model-specific network expansion current constraints
"""
function constraint_ne_model_current(pm::_PM.AbstractPowerModel, n::Int)
end

""
function constraint_ne_voltage_magnitude_difference(pm::_PM.AbstractBFAModel, n::Int, i, f_bus, t_bus, f_idx, t_idx, r, x, g_sh_fr, b_sh_fr, tm)
    branch = _PM.ref(pm, n, :ne_branch, i)
    fr_bus = _PM.ref(pm, n, :bus, f_bus)
    to_bus = _PM.ref(pm, n, :bus, t_bus)
    M_hi   =  fr_bus["vmax"]^2/tm^2 - to_bus["vmin"]^2
    M_lo   = -fr_bus["vmin"]^2/tm^2 + to_bus["vmax"]^2
    p_fr = _PM.var(pm, n, :p_ne, f_idx)
    q_fr = _PM.var(pm, n, :q_ne, f_idx)
    w_fr = _PM.var(pm, n, :w, f_bus)
    w_to = _PM.var(pm, n, :w, t_bus)
    z    = _PM.var(pm, n, :branch_ne, i)

    JuMP.@constraint(pm.model, (w_fr/tm^2) - w_to <= 2*(r*p_fr + x*q_fr) + M_hi*(1-z) )
    JuMP.@constraint(pm.model, (w_fr/tm^2) - w_to >= 2*(r*p_fr + x*q_fr) - M_lo*(1-z) )
end


text:
import Distributions: quantile, Normal

"""
```julia
p₁ = binom_confidence(n::Int,c::Real)
```
Returns the probability `p₁` of successes in the binomial distribution (p=1/2) of
`n` trials with confidence `c`.

# Example
At c=0.95, i.e. 95% confidence, n=1000 tosses of
a coin will yield not more than
```julia
julia> p₁ = BitInformation.binom_confidence(1000,0.95)
0.5309897516152281
```
about 53.1% heads (or tails)."""
function binom_confidence(n::Int,c::Real)
    return 0.5 + quantile(Normal(),1-(1-c)/2)/(2*sqrt(n))
end

"""
```julia
Hf = binom_free_entropy(n::Int,c::Real,base::Real=2)
```
Returns the free entropy `Hf` associated with `binom_confidence`."""
function binom_free_entropy(n::Int,c::Real,base::Real=2)
    p = binom_confidence(n,c)
    return 1 - entropy([p,1-p],base)
end


text:
using Polynomials
using Jacobi
using Test

@testset "Jacobi.jl" begin
    include("testjacobi.jl")
    include("testgauss.jl")
end


text:
# These examples are documented by the LCP notebooks from
# https://github.com/QuantEcon/SimpleDifferentialOperators.jl
# In general, these are optimal stopping problems solved
# as linear complementarity problems.

using Parameters, DiffEqOperators, Random, Suppressor, Test, PATHSolver

# Seed
Random.seed!(1793)

# Setup
StoppingProblem = @with_kw (μ_bar = -0.01, # 1D Brownian motion with drift
                            σ_bar = 0.01,
                            S_bar = 10.0,
                            γ = 0.5, # u(x) = x^γ
                            ρ = 0.05, # discount rate
                            x_min = 0.01,
                            x_max = 5.0,
                            M = 15) # num of grid points

function LCP_objects(sp)
    # setup
    @unpack μ_bar, σ_bar, S_bar, γ, ρ, x_min, x_max, M = sp
    grid = range(x_min, x_max, length = M)
    grid_extended = [grid[1] - diff(grid)[1]; grid; grid[end] + diff(grid)[end]]
    μ = μ_bar
    S(x) = S_bar
    u(x) = x^γ
    σ = σ_bar

    # construct operator and non-operator LCP objects
    BC1 = RobinBC((0., 1., 0.), (0., 1., 0.), grid[2] - grid[1], 1) # for first-derivative
    BC2 = RobinBC((0., 1., 0.), (0., 1., 0.), grid[2] - grid[1], 1) # for second-derivative
    L1 = UpwindDifference(1, 1, grid[2] - grid[1], M, μ) * BC1
    L2 = σ^2/2 * CenteredDifference(2, 2, grid[2] - grid[1], M) * BC2
    S_vec = S.(grid)
    u_vec = u.(grid)

    return (L1 = Array(L1)[1], L2 = Array(L2)[1], S = S_vec, u = u_vec, ρ = ρ)
end

function LCPsolve(sp)
    @unpack L1, L2, S, u, ρ = LCP_objects(sp)
    f = z -> ρ .* z - L1 * z - L2 * z -
        u + ρ .* S - L1 * S - L2 * S
    n = sp.M
    lb = zeros(n)
    ub = 300 .* ones(n) # a reasonable guess
    options(convergence_tolerance = 1e-15, output = :no,
            time_limit = 600) # 10 minute budget
    exit_code, sol_z, sol_f = @suppress solveLCP(f, lb, ub)
end

@testset "Backward Case, Nonzero S" begin
  sp = StoppingProblem()
  code, sol, f = LCPsolve(sp)
  @test code == :Solved
  @test sol[1] ≈ 0.0 atol = 1e-5
  @test sol[5] ≈ 12.066957816758809 atol = 1e-5
  @test f[1] ≈ 0.39946441597137766 atol = 1e-5
end

@testset "Backward Case, Zero S" begin
  code, sol, f = LCPsolve(StoppingProblem(S_bar = 0.))
  @test code == :Solved
  @test sol[1] ≈ 2.050665004133949 atol = 1e-5
  @test sol[8] ≈ 30.258918534086924 atol = 1e-5
  @test f[1] ≈ 0.0 atol = 1e-2
end

nothing


text:
using AoC2021
using Test

@testset "AoC2021.jl" begin
    # Write your tests here.
end


text:
# Should be run from EchelleCCFs/deps or EchelleCCFs/examples or EchelleCCFs/test
# Since this May be called before EchelleCCFs is installed, so this doesn't use pkgdir.

include("download.jl")

download_url = "https://zenodo.org/record/3753254/files/res-1000-1years_full_id1.h5?download=1"
download_filename = joinpath("..","data","spectra","res-1000-1years_full_id1.h5")
download_md5 = "5659082144cd093d617bb54dca937ad9"

@warn "This is a large download.  Be prepared to be patient."
download_and_check_md5sum(download_url, download_filename, md5_goal=download_md5)

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# ============================================================
# Boundary Methods
# ============================================================

include("bc_extra.jl")
include("bc_balance.jl")
include("bc_period.jl")
include("bc_mirror.jl")
include("bc_isothermal.jl")


text:
@testset "Karger Minimum Cut" begin

    gx = path_graph(5)

    #Assumes cut[1] = 1
    for g in testgraphs(gx)
        cut = @inferred(karger_min_cut(g))
        @test findfirst(isequal(2), cut) ==  findlast(isequal(1), cut)+1
        @test karger_cut_cost(g, cut) == 1
        @test karger_cut_edges(g, cut) == [Edge(findlast(isequal(1), cut), findfirst(isequal(2), cut))]
    end

    add_vertex!(gx)

    for g in testgraphs(gx)
        cut = @inferred(karger_min_cut(g))
        @test cut == [1, 1, 1, 1, 1, 2]
        @test karger_cut_cost(g, cut) == 0
        @test karger_cut_edges(g, cut) == Vector{Edge}()

    end

    gx = star_graph(5)
    for g in testgraphs(gx)
        cut = @inferred(karger_min_cut(g))
        @test count(isequal(2), cut) == 1
        @test karger_cut_cost(g, cut) == 1
        @test karger_cut_edges(g, cut) == [Edge(1, findfirst(isequal(2), cut))]
    end

    gx = SimpleGraph(1)
    for g in testgraphs(gx)
        cut = @inferred(karger_min_cut(g))
        @test cut == zeros(Int, 1)
    end

end


text:
function ADforward!(net;returnf=false,debug=false,AllocateMemory=false)
#=   calculate the value and auxiliary values of a net function
=#
    # (c) David Barber, University College London 2015

    if AllocateMemory & net.gpu
        TransformToGPU=true
        net=convert(net,"CPU")
    else
        TransformToGPU=false
    end

    if debug; println("Get value:"); end

    if isempty(net.auxvalue)
        net.auxvalue=similar(net.value)
    end

    # forward pass:
    for i in net.ForwardPassList
        thisnode=net.node[i]
        if debug
            println("node $i: $(thisnode.f)($(thisnode.parents))");
        end
        if AllocateMemory
            if debug
                @time net.value[i],net.auxvalue[i]=thisnode.f(net.value[thisnode.parents]...) ## not in place
            else
                net.value[i],net.auxvalue[i]=thisnode.f(net.value[thisnode.parents]...) ## not in place
            end
        else
            if debug
                println("in place")
                @time thisnode.f_inplace(net.value[i],net.auxvalue[i],net.value[thisnode.parents]...) ## in place
            else
                thisnode.f_inplace(net.value[i],net.auxvalue[i],net.value[thisnode.parents]...) ## in place
            end
        end
    end

    if TransformToGPU
        net=convert(net,"GPU")
    end

    if returnf
        if net.gpu
            return to_host(net.value[net.FunctionNode])[1]
        else
            return net.value[net.FunctionNode][1]
        end
#        @gpu if net.gru
#        @cpu return net.value[net.FunctionNode][1]
    end

    if debug; println("done"); end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module SequentialTests

using Test

@testset "SequentialBackend" begin include("SequentialBackendTests.jl") end

@testset "Interfaces" begin include("InterfacesTests.jl") end

@testset "FDM" begin include("FDMTests.jl") end

@testset "FEMSA" begin include("FEMSATests.jl") end

@testset "PTimers" begin include("PTimersTests.jl") end

end # module


text:
"""
Example taken from:
https://github.com/JuliaDiffEq/ModelingToolkit.jl/issues/36#issuecomment-536221300
"""

include("NewNamedViewVector.jl")

using DifferentialEquations
using Plots

# theme(:juno)

function lorenz!(du, u, (p, f), t)
    x, y, z = u.x, u.y, u.z
    σ, ρ, β = p.σ, p.ρ, p.β

    du.x = σ*(y - x)
    du.y = x*(ρ - z) - y + f
    du.z = x*y - β*z
    return nothing
end

function lotka!(du, u, (p, f), t)
    x, y = u.x, u.y
    α, β, γ, δ = p.α, p.β, p.γ, p.δ

    du.x = α*x - β*x*y + f
    du.y = -γ*y + δ*x*y
    return nothing
end

function composed!(du, u, p, t)
    lorenz!(du.lorenz, u.lorenz, (p, u.lotka.x), t)
    lotka!(du.lotka, u.lotka, (p, u.lorenz.x), t)
    return nothing
end

p = (α=1.0, β=1.0, γ=3.1, δ=0.5, ρ=1.0, σ=1.0)

prange = 0 : 0.01 : 10

tspan = (0.0, 20.0)
lorenz_ic = (x=0.0, y=0.0, z=0.0)
lotka_ic = (x=1.0, y=1.0)

u₀ = NamedViewVector{Float64}((lorenz=lorenz_ic, lotka=lotka_ic))

prob = ODEProblem(composed!, u₀, tspan, p)
sol = solve(prob, Tsit5())

plot(sol)

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
using Optim, LinearAlgebra

# Augmented Lagrange Method code for inequality contrataint
function Augmented_Lagrange_Mult_Inequality_Box(f, g, lower, upper, x, k_max; α=1000., γ=1.01, tol=1e-1000)
    λ, nite = zeros(length(g(x))), 0
    for k in 1 : k_max; nite = k
        ψ(x) = max.(g(x),-λ./(2α))
        F(x) = f(x) + λ ⋅ g(x) + α * sum(ψ(x).^2)
        #x = Optim.minimizer(optimize(F, lower, upper, x, Fminbox(NelderMead()))) # Ordem zero
        x = Optim.minimizer(optimize(F, lower, upper, x, Fminbox(BFGS());autodiff = :forward)) # Ordem 1
        λ .+= 2*α*ψ(x)
        if norm(ψ(x),Inf)<=tol; break; end # Ver Arora pg 481
        α *= γ;
    end
    return x, nite, α, λ
end;

Sy=350.0;
r1=0.4;
E=2000.0;
P=200.0;
B=1.0;

f(x) = pi*(x[1]^2)*(x[2]^2+B^2)^(1/2);

g1(x) = P*B/(pi*x[2]*(r1^2)*Sy) - 1;
g2(x) = P*((B^2+x[2]^2)^(1/2))/(x[2]*pi*x[1]^2*Sy)-1;
g3(x) = 4*P*(B^2+x[2]^2)^(3/2)/(x[2]*pi^3*x[1]^4*E)-1;

g(x) = [g1(x), g2(x), g3(x)];

lb = [0.4, 1.0];
ub = [0.8, 6.0];

x0 = [0.5, 1.1]

@time xₑ, k, α, λ = Augmented_Lagrange_Mult_Inequality_Box(f, g, lb,ub, x0, 100)
println("x* : ", xₑ, ", iterations: ", k, ", α : ", α, ", λ : ", λ)


text:
using GeneralizedEigenvalueMinimization
using Test

@testset "GeneralizedEigenvalueMinimization.jl" begin
    # Write your tests here.
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()
config_with_basis = MOIT.TestConfig(basis = true)

@testset "Split" begin
    T = Float64
    bridged_mock = MOIB.Constraint.SplitInterval{T}(mock)
    MOIT.basic_constraint_tests(
        bridged_mock,
        config,
        include = [
            (MOI.SingleVariable, MOI.Interval{T}),
            (MOI.ScalarAffineFunction{T}, MOI.Interval{T}),
            (MOI.ScalarQuadraticFunction{T}, MOI.Interval{T}),
            (MOI.SingleVariable, MOI.EqualTo{T}),
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
            (MOI.ScalarQuadraticFunction{T}, MOI.EqualTo{T}),
            (MOI.VectorOfVariables, MOI.Zeros),
            (MOI.VectorAffineFunction{T}, MOI.Zeros),
            (MOI.VectorQuadraticFunction{T}, MOI.Zeros),
        ],
    )
end

@testset "my test" begin
    bridged_mock = MOIB.Constraint.SplitInterval{Float64}(mock)

    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, -1.0, 1.0],
        ),
    )

    # variables: x, y, z
    # maxobjective: x + -1y + z
    # c1: x in Interval(-Inf, 1.0)
    # c2: y in Interval(-1.0, Inf)
    # c3: z in Interval(-Inf, Inf)
    # c4: 1z <= 1.0
    MOIT.solve_one_sided_intervals(bridged_mock, config_with_basis)

    println("====== ****** ======")
    for (i, b) in bridged_mock.map.single_variable_constraints
        @show MOI.get(bridged_mock.model, MOI.ConstraintSet(), b.lower)
        @show MOI.get(bridged_mock.model, MOI.ConstraintSet(), b.upper)
    end

    println("====== ****** ======")

    @show MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.SingleVariable,
                MOI.Interval{Float64},
            }(),
        )

end
exit()

@testset "Interval" begin
    bridged_mock = MOIB.Constraint.SplitInterval{Float64}(mock)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [5.0, 5.0],
            con_basis = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.NONBASIC],
                (MOI.SingleVariable, MOI.GreaterThan{Float64}) =>
                    [MOI.BASIC, MOI.BASIC],
            ],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [2.5, 2.5],
            con_basis = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.NONBASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.BASIC],
                (MOI.SingleVariable, MOI.GreaterThan{Float64}) =>
                    [MOI.BASIC, MOI.BASIC],
            ],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 1.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
            con_basis = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.NONBASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.BASIC],
                (MOI.SingleVariable, MOI.GreaterThan{Float64}) =>
                    [MOI.BASIC, MOI.BASIC],
            ],
        ),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [6.0, 6.0],
            con_basis = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.NONBASIC],
                (MOI.SingleVariable, MOI.GreaterThan{Float64}) =>
                    [MOI.BASIC, MOI.BASIC],
            ],
        ),
    )
    MOIT.linear10test(bridged_mock, config_with_basis)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [0.0, 0.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1],
            con_basis = [
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [MOI.BASIC],
                (MOI.SingleVariable, MOI.GreaterThan{Float64}) =>
                    [MOI.NONBASIC, MOI.NONBASIC],
            ],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [0],
        ),
    )
    MOIT.linear10btest(bridged_mock, config_with_basis)

    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.Interval{Float64},
            }(),
        ),
    )
    vis = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    newf =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], vis), 0.0)
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, newf)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ newf

    MOI.modify(bridged_mock, ci, MOI.ScalarCoefficientChange(vis[2], 1.0))
    modified_f =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(2), vis), 0.0)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ modified_f

    @testset "$attr" for attr in [
        MOI.ConstraintPrimalStart(),
        MOI.ConstraintDualStart(),
    ]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, 2.0)
        @test MOI.get(bridged_mock, attr, ci) == 2.0
        bridge = MOIB.bridge(bridged_mock, ci)
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(bridged_mock, attr, bridge.lower) == 2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == 2.0
        else
            @test MOI.get(bridged_mock, attr, bridge.lower) == 2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == 0.0
        end
        MOI.set(bridged_mock, attr, ci, -2.0)
        @test MOI.get(bridged_mock, attr, ci) == -2.0
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(bridged_mock, attr, bridge.lower) == -2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == -2.0
        else
            @test MOI.get(bridged_mock, attr, bridge.lower) == 0.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == -2.0
        end
    end

    test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
        ),
    )
end
@testset "EqualTo{$T}" for T in [Float64, Int]
    bridged_mock = MOIB.Constraint.SplitInterval{T}(mock)
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [one(T), one(T)],
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}) =>
                zeros(T, 2),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => zeros(T, 1),
        ),
    )
    MOIT.linear13test(bridged_mock, MOIT.TestConfig{T}())

    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{T},
                MOI.EqualTo{T},
            }(),
        ),
    )

    test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}, 1),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}, 0),
        ),
    )
end
@testset "Zeros" begin
    bridged_mock = MOIB.Constraint.SplitInterval{Float64}(mock)
    mock.optimize! =
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            [1.0, 0.0, 2.0],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[0.0, 0.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) =>
                [[-3, -1]],
        )
    MOIT.lin1vtest(bridged_mock, config)

    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ),
    )

    test_delete_bridge(
        bridged_mock,
        ci,
        3,
        (
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives, 0),
        ),
    )
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3e849bdf-f13a-4ea0-aa54-b6ae89894565
md"""
# Álgebra lineal básica en Julia
Autor: Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/)
(con edición de Jane Herriman)
"""

# ╔═╡ a880f5b4-1f0f-4569-a59c-0e87d21cd203
md"""
Julia's syntax is very similar to other languages but there are some important differences. Define a matrix of random normal variates
"""

# ╔═╡ 05eff53f-64c5-4833-9eff-11c7d09b9684
A = rand(1:4,3,3)

# ╔═╡ 16ca2d04-2f46-4de7-9481-daebf5112c2f
md"""
Definir un vector de unos
"""

# ╔═╡ 018e028c-43ff-4dbe-bc83-e5f0982d6bb6
x = fill(1.0, (3))

# ╔═╡ 5cc332db-a229-4182-b860-4d276e3042a9
md"""
Notamos que $A$ tiene el tipo Array{Int64,2} pero $x$ tiene tipo Array{Int64,1}. Julia define los alias a Vector{Type}=Array{Type,1} y Matrix{Type}=Array{Type,2}.

Muchas de las operaciones básicas son idénticas a otros lenguajes

#### Multiplicación
"""

# ╔═╡ 383658c8-9d47-4a1e-ba63-66fa6e7f3311
b = A*x

# ╔═╡ e2cc7ff7-781c-4102-892c-da24291f0388
md"""
#### Traspuestas
Como en otros lenguajes, `A'` es la transpuesta conjugada mientras que `A.'` es sólo la traspuesta
"""

# ╔═╡ 5c3d3b48-6b7c-4a3b-a9b6-d83f7bc03168
Asym = A + A'

# ╔═╡ 88f5c70e-b5b1-4cf0-8f50-d1241f0a769b
md"""
#### Multiplicación traspuesta
Julia nos permite escribir esto sin *
"""

# ╔═╡ 0691212e-d372-43be-b66e-61f842df09eb
Apd = A'A

# ╔═╡ 1445f674-cfa2-44ab-920f-fe4d85d00c2e
md"""
#### Resolviendo sistemas lineales
El problema $Ax=b$ para $A$ cuadrada se resulve con la función \.
"""

# ╔═╡ c42fc4ae-bf1d-4eb7-9863-a498fb9116b8
A\b

# ╔═╡ 462e5a23-f323-4158-bf27-cf36bbc5cec5
md"""
#### Sistemas sobredeterminados
Cuando nuestra matriz es alta (número de renglones mayores al número de columnas), tenemos un sistema lineal sobredeterminado.

En este caso \ calcula la de mínimos cuadrados
"""

# ╔═╡ 8c87b3e1-b737-4179-9955-defcf0905b1b
Atall = rand(3, 2)
display(Atall)
Atall\b

# ╔═╡ 7dea4ca8-27e0-429e-bbf3-dccf2ca67df3
md"""
La función \ también sirve ocn problemas deficientes de rango de mínimos cuadrados. En este caso, la solución no es única y Julia regresa el valor con la menor norma.

Para crear un problema de rango deficiente de mínimos cuadrados, vamos a crear una matriz deficiente en rango con columnas linealmente dependientes
"""

# ╔═╡ 962df670-f9d1-441d-bf86-f4619f29ab4f
v = randn(3)

# ╔═╡ 306db6d0-193e-11eb-3bee-cbb7d8c4c5c6
rankdef = [v v]

# ╔═╡ 92b4cb6a-47c9-459b-8694-aee9158439dd
rankdef\b

# ╔═╡ 4c836857-08fc-411a-8b65-6340c6e94c3b
md"""
#### Sistemas indeterminados
cuando A es corta (número de columnas mayor al número de renglones), tenemos un sistema indeterminado

En este caso \ regresa la solución con la norma mínima
"""

# ╔═╡ 0784963e-8842-43f6-af98-781aa7e214b6
Ashort = rand(2, 3)

# ╔═╡ 3346487c-193e-11eb-15dd-4dd8241ac876
display(Ashort)

# ╔═╡ 3346b974-193e-11eb-2c77-9133e1f9b17b
Ashort\b[1:2]

# ╔═╡ 86dc3af8-2f80-48a5-b0d3-caeddc5ae8a8
md"""
### Ejercicios

```
A = [
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 1  2  3  4  5  6  7  8  9  10
 ]
```

Quieres obtener

```
A = [
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 7  8  9  10  1  2  3  4  5  6
 ]
```
"""

# ╔═╡ f86a1d32-1e1d-461e-ad40-f5b906b04cf9

# ╔═╡ 22e5c662-6bf6-48d4-b7f9-93d8ac33ffdd
md"""
10.2 Toma el producto de un vector `v` con sí mismo.
"""

# ╔═╡ 403752d3-4a6b-41bf-b4ed-fd6fe7b9ee95

# ╔═╡ 6b98b27f-f671-4c5c-af5e-6d290f0d158b
md"""
10.3 Toma el producto de un vector `v` con sí mismo.
"""

# ╔═╡ fd670b40-f44e-440e-90b0-b27522c7adbd

# ╔═╡ Cell order:
# ╟─3e849bdf-f13a-4ea0-aa54-b6ae89894565
# ╟─a880f5b4-1f0f-4569-a59c-0e87d21cd203
# ╠═05eff53f-64c5-4833-9eff-11c7d09b9684
# ╟─16ca2d04-2f46-4de7-9481-daebf5112c2f
# ╠═018e028c-43ff-4dbe-bc83-e5f0982d6bb6
# ╟─5cc332db-a229-4182-b860-4d276e3042a9
# ╠═383658c8-9d47-4a1e-ba63-66fa6e7f3311
# ╟─e2cc7ff7-781c-4102-892c-da24291f0388
# ╠═5c3d3b48-6b7c-4a3b-a9b6-d83f7bc03168
# ╟─88f5c70e-b5b1-4cf0-8f50-d1241f0a769b
# ╠═0691212e-d372-43be-b66e-61f842df09eb
# ╟─1445f674-cfa2-44ab-920f-fe4d85d00c2e
# ╠═c42fc4ae-bf1d-4eb7-9863-a498fb9116b8
# ╟─462e5a23-f323-4158-bf27-cf36bbc5cec5
# ╠═8c87b3e1-b737-4179-9955-defcf0905b1b
# ╟─7dea4ca8-27e0-429e-bbf3-dccf2ca67df3
# ╠═962df670-f9d1-441d-bf86-f4619f29ab4f
# ╠═306db6d0-193e-11eb-3bee-cbb7d8c4c5c6
# ╠═92b4cb6a-47c9-459b-8694-aee9158439dd
# ╟─4c836857-08fc-411a-8b65-6340c6e94c3b
# ╠═0784963e-8842-43f6-af98-781aa7e214b6
# ╠═3346487c-193e-11eb-15dd-4dd8241ac876
# ╠═3346b974-193e-11eb-2c77-9133e1f9b17b
# ╟─86dc3af8-2f80-48a5-b0d3-caeddc5ae8a8
# ╠═f86a1d32-1e1d-461e-ad40-f5b906b04cf9
# ╟─22e5c662-6bf6-48d4-b7f9-93d8ac33ffdd
# ╠═403752d3-4a6b-41bf-b4ed-fd6fe7b9ee95
# ╟─6b98b27f-f671-4c5c-af5e-6d290f0d158b
# ╠═fd670b40-f44e-440e-90b0-b27522c7adbd

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
# double run_P(int N1,int N2,double dt,double *spike_times_1,double *spike_times_2){
function run_P(spike_times_1::Vector{<:Real},
                spike_times_2::Vector{<:Real},
                dt::Real,
            )::Float64
    Nab = 0
    N1 = length(spike_times_1)
    N2 = length(spike_times_2)

    spike_index2 = 1
    for spike_index1 = 1:N1
        while (spike_index2 < N2)
            # check every spike in train 1 to see if there's a spike in train 2 within dt  (don't count spike pairs)
            # don't need to search all spike_index2 each iteration
            if abs(spike_times_1[spike_index1] - spike_times_2[spike_index2]) <=
               dt
                Nab = Nab + 1
                break
            elseif spike_times_2[spike_index2] > spike_times_1[spike_index1]
                break
            else
                spike_index2 = spike_index2 + 1
            end # if
        end # while
    end # for

    return Nab
end

# double run_T(int N1v,double dtv,double startv, double endv, double *spike_times_1)
"""
dt is the window size used for searching for spike cooccurance

TODO: what is missing now is the how the vectors should be defined, because this would remove ambiguity about the time resolutions
"""
function run_T(spike_times_1::Vector{<:Real},
                dt::Real,
                start_time::Real,
                stop_time::Real,
            )::Float64
    N1 = length(spike_times_1)
    vec_index = 1

    # maximum
    time_A = 2 * N1 * dt

    #  if just one spike in train
    if N1 == 1
        first_spike_time = spike_times_1[0]
        if first_spike_time - start_time < dt
            time_A = time_A - start_time + first_spike_time - dt
        elseif first_spike_time + dt > stop_time
            time_A = time_A - first_spike_time - dt + stop_time
        end
    else # if more than one spike in train
        while vec_index < N1
            dif_index = vec_index + 1
            diff = spike_times_1[dif_index] - spike_times_1[vec_index]

            if diff < 2 * dt # subtract overlap
                time_A = time_A - 2 * dt + diff
            end
            vec_index += 1
        end

        # check if spikes are within dt of the start_time and/or end, if so just need to subract
        # overlap of first and/or last spike as all within-train overlaps have been accounted for
        if spike_times_1[1] - start_time < dt
            time_A = time_A - start_time + spike_times_1[1] - dt
        end

        if stop_time - spike_times_1[N1-1] < dt
            time_A = time_A - spike_times_1[N1] - dt + stop_time
        end
    end

    return time_A
end

# void run_sttc(int *N1v,int *N2v,double *dtv,double *Time,double *index,double *spike_times_1,double *spike_times_2){
function tailing_coefficient(spike_times_1::Vector{<:Real},
                                spike_times_2::Vector{<:Real},
                                dt::Real,
                                start_time::Real,
                                stop_time::Real,
                            )

    N1 = length(spike_times_1)
    N2 = length(spike_times_2)

    if (N1 == 0 || N2 == 0)
        index = NaN
    else
        T = stop_time - start_time
        # TA=run_T(N1,dt,Time[0],Time[1], spike_times_1) / T
        # TB=run_T(N2,dt,Time[0],Time[1], spike_times_2) / T
        # PA=run_P(N1,N2,dt, spike_times_1, spike_times_2) / N1
        # PB=run_P(N2,N1,dt, spike_times_2, spike_times_1) / N2

        TA = run_T(spike_times_1, dt, start_time, stop_time) / T
        TB = run_T(spike_times_2, dt, start_time, stop_time) / T
        PA = run_P(spike_times_1, spike_times_2, dt) / N1
        PB = run_P(spike_times_2, spike_times_1, dt) / N2

        index =
            0.5 * (PA - TB) / (1 - TB * PA) + 0.5 * (PB - TA) / (1 - TA * PB)
    end

    return index
end


text:
#This is a modified version of FixedSizeStrings.jl
#adapted so that appended zero Chars are omitted from
#the string. So most credits got to authors of
#https://github.com/JuliaComputing/FixedSizeStrings.jl
# Maybe this should be moved to FizedSizeStrings.jl,
#but until then let's keep it here...
module MaxLengthStrings
import Base: iterate, lastindex, getindex, sizeof, length, ncodeunits, codeunit, isvalid, read, write
export MaxLengthString

struct MaxLengthString{N,T} <: AbstractString
    data::NTuple{N,T}
    function MaxLengthString{N,T}(itr) where {N,T}
      new(totuple_appendzero(NTuple{N,T},itr))
    end
end
import Base: tuple_type_head, tuple_type_tail
#totuple_appendzero(::Type{Tuple{}}, itr, s...) = ()
tuple_append(::Type{NTuple{N,T}}) where {N,T} = (zero(T), tuple_append(NTuple{N-1,T})...)
tuple_append(::Type{Tuple{}}) = ()
function totuple_appendzero(::Type{Tuple{}},itr,s...)
  iterate(itr, s...)===nothing || error("String is too long to fit into MaxLengthString")
  ()
end
function totuple_appendzero(::Type{NTuple{N,T}}, itr, s...)::Tuple{Vararg{T,N}} where {N,T}
    y = iterate(itr, s...)
    if y === nothing
      tuple_append(NTuple{N,T})
    else
      (convert(T, y[1]), totuple_appendzero(NTuple{N-1,T}, itr, y[2])...)
    end
end

function MaxLengthString(s::AbstractString,N=length(s),T=UInt8)
  MaxLengthString{N,T}(rpad(s,N,'\0'))
end

function iterate(s::MaxLengthString{N}, i::Int = 1) where N
    i > N && return nothing
    c = s.data[i]
    iszero(c) && return nothing
    return (Char(c), i+1)
end

lastindex(s::MaxLengthString{N}) where {N} = findlast(!iszero,s.data)

function getindex(s::MaxLengthString, i::Int)
  checkbounds(s,i)
  Char(s.data[i])
end

sizeof(s::MaxLengthString) = sizeof(s.data)

length(s::MaxLengthString) = findlast(!iszero,s.data)

ncodeunits(s::MaxLengthString) = length(s)

codeunit(::MaxLengthString{<:Any,T}) where T = T
function codeunit(s::MaxLengthString, i::Integer)
  checkbounds(s,i)
  s.data[i]
end

isvalid(s::MaxLengthString{<:Any,UInt8}, i::Int) = checkbounds(Bool, s, i)
isvalid(s::MaxLengthString, i::Int) = checkbounds(Bool, s, i) && isvalid(Char,s.data[i])

function read(io::IO, T::Type{<:MaxLengthString{N}}) where N
    return read!(io, Ref{T}())[]::T
end

function write(io::IO, s::MaxLengthString{N}) where N
    return write(io, Ref(s))
end
end


text:
using Test

struct MyType
end

@testset "IterableTables" begin

include("test_integration_dataframes.jl")
# include("test_integration_datastreams.jl")
# TODO Reenable once GLM is 0.7 compatible
# include("test_integration_statsmodels.jl")
# include("test_integration_timeseries.jl")
# include("test_integration_temporal.jl")
# include("test_integration_juliadb.jl")
# include("test_integration_gadfly.jl")

end


text:
"""
    validate(bt::XLSXTable)
    validate(jws::JSONWorksheet)

validate
"""
function validate(tb::XLSXTable)
    for sheet in sheetnames(tb)
        errors = validate(tb, sheet)
        if !isempty(errors)
            msgs = pretty_schemaerror(tb, sheet, errors)
            for (msg_error, msg_sol) in msgs
                print_section(msg_error[1], msg_error[2]; color=:yellow)
                printstyled(msg_sol, "\n"; color=:red)
            end
        end
    end
    nothing
end

function validate(tb::XLSXTable, sheet)
    schema = tb.schemas[sheet]
    if !ismissing(schema)
        update!(schema)
        return validate(tb[sheet], schema)
    end
    return []
end

function validate(jws::JSONWorksheet, schema::SchemaData)
    err = OrderedDict()
    @inbounds for (i, row) in enumerate(jws)
        val = JSONSchema.validate(row, schema.data)
        if !isnothing(val)
            marker = "$i"
            err[marker] = val
        end
    end
    return err
end

function pretty_schemaerror(tb::XLSXTable, sheet, err::AbstractDict)
    paths = map(el -> el.path, values(err))

    msgs = []
    for p in unique(paths)
        # error가 난 데이터 내용
        cause = []
        for el in err
            if el[2].path == p
        x = (el[1], el[2].x)
                push!(cause, x)
        end
        end
        # error 원인, 값
        errors = filter(el -> el.path == p, collect(values(err)))
        reason = unique(map(el -> el.reason, values(errors)))
        schemaval = unique(map(el -> el.val, values(errors)))

        error_info = [
        """
        $(tb) Validation failed from {key: $reason, summary: $(summary(schemaval[1]))}
        sheet:        $sheet
        column:       $p
        instance:     $cause
        """,
        "$(tb) Validation failed from {key: $reason, summary: $(summary(schemaval[1]))}"]

        solution = get_schema_description(tb, sheet, p)
        if !ismissing(solution)
            sol_info = "해결방법\n  ↳ $(solution)"
        else
            if length(reason) == 1
                sol_info = "\"$(reason[1])\": $(schemaval[1])"
            else
                sol_info = "$(reason): $(schemaval)"
            end
            if length(sol_info) > 42
                sol_info = sol_info[1:40] * "......"
            end
        end
        push!(msgs, (error_info, sol_info))
    end
    return msgs
end

function get_schema_description(tb::XLSXTable, sheet, path)
    schema = tb.schemas[sheet]
    desc = get_schema_description(schema, path)
end
function get_schema_description(schema::SchemaData, path)
    d = get_schemaelement(schema.data, path)
    if isa(d, AbstractDict)
        get(d, "description", missing)
    else
        missing
    end
end
function get_schemaelement(schema, path)
    # patternProperties는 element를 찾지 않는다
    if !haskey(schema.data, "properties")
        return missing
    end
    wind = schema.data["properties"]
    paths = filter(!isempty, replace.(split(path, "]"), "[" => ""))
    for (i, p) in enumerate(paths)
        if ismissing(wind)
            break
        end
        if i == 1
            wind = get(wind, p, missing)
        elseif occursin(r"^\d+$", p)
            wind = get(wind, "items", missing)
            if isa(wind, AbstractArray)
                idx = parse(Int, p)
                if idx > lastindex(wind)
                    throw(BoundsError(schema, paths))
                end
                wind = wind[parse(Int, p)]
            end
        else
            if haskey(wind, "\$ref")
                if haskey(wind["\$ref"], "properties")
                    wind = get(wind["\$ref"]["properties"], p, missing)
                elseif haskey(wind["\$ref"], "patternProperties")
                    wind = get(wind["\$ref"]["patternProperties"], p, missing)
                else
                    wind = missing
                end
            else
                wind = get(wind["properties"], p, missing)
            end
        end
    end
    return wind
end


text:
module DualArray

import TensorOperations: IndexTuple
import TensorOperations: contract!
import ForwardDiff: Dual
import ForwardDiff

struct DualTensor{N, T}
    A::Array{T}
    ∂::NTuple{N, Array{T}}

    "Construct an instance of `DualTensor` from `Array` components."
    DualTensor(A::Array{T},
               ∂::NTuple{N, Array{T}}) where {N, T} = begin
        all(map(sz -> sz == size(A), size.(∂))) ||
            throw(DimensionMismatch("Tensor data and its derivatives must have the same size."))

        new{N, T}(A, ∂)
    end
end

"Convert `Array` of `Dual`s into `DualTensor`."
DualTensor(A::Array{Dual{Tag, T, N}}) where {N, T, Tag} = begin
    AA = ForwardDiff.value.(A)
    ∂A = (( (x -> x[i∂]).(ForwardDiff.partials.(A)) for i∂=1:N )...,)

    DualTensor(AA, ∂A)
end

include("tensor_opr.jl")

end


text:
"""
    knapsack_solver(array::Array{Int64,1}, weight::Array{Int64,1}, capacity::Int64)

The knapsack problem is linked to the combinatorial optimization. In more detail all
array-elements are characterized by a given weight. The excersive is to maximize the the
value of subsequence of array-elements by not ooutbiding the `capacity`. For this purpose a
recurisve comparsion between old and new added items is designed wich continously update the
maximum reacahbele value sum. For more information see: [https://en.wikipedia.org/wiki/Knapsack_problem](https://en.wikipedia.org/wiki/Knapsack_problem)

# Arguments
- `array::Array{Int64,1}`: Array with the current values of each items.
- `weight::Array{Int64,1}`: Array with the current weoght of each items.
- `capacity::Int64`: Maximum weight to carry.

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> val = [60, 100, 120, 230, 220]
julia> wt = [10, 20, 25, 30, 25]
julia> W = 50
julia> ClassicAlgorithmsCollections.knapsack_solver(val, wt, W)
340
```
"""
function knapsack_solver(array::Array{Int64,1}, weight::Array{Int64,1}, capacity::Int64)
    n = length(array)

    if n == 1 || capacity == 0
        return 0
    end

    if (weight[n] > capacity)
        return knapsack_solver(array[1:n-1], weight[1:n], capacity)

    else
        item_new =
            array[n] + knapsack_solver(array[1:n-1], weight[1:n], capacity - weight[n])
        item_old = knapsack_solver(array[1:n-1], weight[1:n], capacity)
        if item_new > item_old
            return item_new
        end
        return item_old
    end
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using OrdinaryDiffEq, Flux, Optim, Random, Plots
using DiffEqSensitivity
using Zygote
using ForwardDiff
using LinearAlgebra, Statistics
using ProgressBars, Printf
using Flux.Optimise: update!, ExpDecay
using Flux.Losses: mae, mse
using BSON: @save, @load

Random.seed!(1234);

###################################
# Arguments
is_restart = true;
n_epoch = 10000;
n_plot = 50;
datasize = 50;
tstep = 1;
n_exp_train = 20;
n_exp_test = 10;
n_exp = n_exp_train + n_exp_test;
noise = 0.05;
ns = 6;
nr = 3;
alg = AutoTsit5(Rosenbrock23(autodiff=false));
atol = 1e-6;
rtol = 1e-3;

# opt = ADAMW(5.f-3, (0.9, 0.999), 1.f-6);
opt = Flux.Optimiser(ExpDecay(5e-3, 0.5, 500 * n_exp_train, 1e-4),
                     ADAMW(0.005, (0.9, 0.999), 1.f-6));

lb = 1.f-6;
ub = 1.f1;
####################################

function trueODEfunc(dydt, y, k, t)
    # TG(1),ROH(2),DG(3),MG(4),GL(5),R'CO2R(6)
    r1 = k[1] * y[1] * y[2];
    r2 = k[2] * y[3] * y[2];
    r3 = k[3] * y[4] * y[2];
    dydt[1] = - r1;  # TG
    dydt[2] = - r1 - r2 - r3;  # TG
    dydt[3] = r1 - r2;  # DG
    dydt[4] = r2 - r3;  # MG
    dydt[5] = r3;  # GL
    dydt[6] = r1 + r2 + r3;  # R'CO2R
    dydt[7] = 0.f0;
end

logA = Float32[18.60f0, 19.13f0, 7.93f0];
Ea = Float32[14.54f0, 14.42f0, 6.47f0];  # kcal/mol

function Arrhenius(logA, Ea, T)
    R = 1.98720425864083f-3
    k = exp.(logA) .* exp.(-Ea ./ R ./ T)
    return k
end

# Generate datasets
u0_list = rand(Float32, (n_exp, ns + 1));
u0_list[:, 1:2] = u0_list[:, 1:2] .* 2.0 .+ 0.2;
u0_list[:, 3:ns] .= 0.0;
u0_list[:, ns + 1] = u0_list[:, ns + 1] .* 20.0 .+ 323.0;  # T[K]
tspan = Float32[0.0, datasize * tstep];
tsteps = range(tspan[1], tspan[2], length=datasize);

ode_data_list = zeros(Float32, (n_exp, ns, datasize));
yscale_list = [];
function max_min(ode_data)
    return maximum(ode_data, dims=2) .- minimum(ode_data, dims=2) .+ lb
end
for i in 1:n_exp
    u0 = u0_list[i, :]
    k = Arrhenius(logA, Ea, u0[end])
    prob_trueode = ODEProblem(trueODEfunc, u0, tspan, k)
    ode_data = Array(solve(prob_trueode, alg, saveat=tsteps))[1:end - 1, :]
    ode_data += randn(size(ode_data)) .* ode_data .* noise
    ode_data_list[i, :, :] = ode_data
    push!(yscale_list, max_min(ode_data))
end
yscale = maximum(hcat(yscale_list...), dims=2);

np = nr * (ns + 2) + 1;
p = randn(Float32, np) .* 0.1;
p[1:nr] .+= 0.8;
p[nr * (ns + 1) + 1:nr * (ns + 2)] .+= 0.8;
p[end] = 0.1;

function p2vec(p)
    slope = p[nr * (ns + 2) + 1] .* 100
    w_b = p[1:nr] .* slope
    w_out = reshape(p[nr + 1:nr * (ns + 1)], ns, nr)
    w_in_Ea = abs.(p[nr * (ns + 1) + 1:nr * (ns + 2)] .* slope)
    w_in = clamp.(-w_out, 0, 4)
    w_in = vcat(w_in, w_in_Ea')
    return w_in, w_b, w_out
end

function display_p(p)
    w_in, w_b, w_out = p2vec(p);
    println("species (column) reaction (row)")
    println("w_in | w_b")
    w_in_ = vcat(w_in, w_b')'
    show(stdout, "text/plain", round.(w_in_, digits=3))
    println("\nw_out")
    show(stdout, "text/plain", round.(w_out', digits=3))
    println("\n")
end
display_p(p)

inv_R = - 1 / 1.98720425864083f-3;
function crnn(du, u, p, t)
    logX = @. log(clamp(u[1:end - 1], lb, ub))
    w_in_x = w_in' * vcat(logX, inv_R / u[end])
    du .= vcat(w_out * (@. exp(w_in_x + w_b)), 0.f0)
end

u0 = u0_list[1, :];
prob = ODEProblem(crnn, u0, tspan, saveat=tsteps, atol=atol, rtol=rtol)

sense = BacksolveAdjoint(checkpointing=true; autojacvec=ZygoteVJP());
function predict_neuralode(u0, p)
    global w_in, w_b, w_out = p2vec(p)
    pred = clamp.(Array(solve(prob, alg, u0=u0, p=p, sensalg=sense)), -ub, ub)
    return pred
end
predict_neuralode(u0, p)

i_obs = [1, 2, 3, 4, 5, 6];
function loss_neuralode(p, i_exp)
    ode_data = @view ode_data_list[i_exp, i_obs, :]
    pred = predict_neuralode(u0_list[i_exp, :], p)[i_obs, :]
    loss = mae(ode_data ./ yscale[i_obs], pred ./ yscale[i_obs])
    return loss
end

cbi = function (p, i_exp)
    ode_data = ode_data_list[i_exp, :, :]
    pred = predict_neuralode(u0_list[i_exp, :], p)
    l_plt = []
    for i in 1:ns
        plt = scatter(tsteps, ode_data[i,:], markercolor=:transparent,
                      title=string(i), label=string("data_", i))
        plot!(plt, tsteps, pred[i,:], label=string("pred_", i))
        push!(l_plt, plt)
    end
    plt_all = plot(l_plt..., legend=false)
    png(plt_all, string("figs/i_exp_", i_exp))
    return false
end

l_loss_train = []
l_loss_val = []
iter = 1
cb = function (p, loss_train, loss_val)
    global l_loss_train, l_loss_val, iter
    push!(l_loss_train, loss_train)
    push!(l_loss_val, loss_val)

    if iter % n_plot == 0
        display_p(p)
        @printf("min loss train %.4e val %.4e\n", minimum(l_loss_train), minimum(l_loss_val))

        l_exp = randperm(n_exp)[1:1];
        println("update plot for ", l_exp)
        for i_exp in l_exp
            cbi(p, i_exp)
        end

        plt_loss = plot(l_loss_train, xscale=:log10, yscale=:log10,
                        framestyle=:box, label="Training")
        plot!(plt_loss, l_loss_val, label="Validation")
        plot!(xlabel="Epoch", ylabel="Loss")
        png(plt_loss, "figs/loss")

        @save "./checkpoint/mymodel.bson" p opt l_loss_train l_loss_val iter;
    end
    iter += 1;
end

if is_restart
    @load "./checkpoint/mymodel.bson" p opt l_loss_train l_loss_val iter;
    iter += 1;
end

i_exp = 1
epochs = ProgressBar(iter:n_epoch);
loss_epoch = zeros(Float32, n_exp);
grad_norm = zeros(Float32, n_exp_train);
for epoch in epochs
    global p
    for i_exp in randperm(n_exp_train)
        grad = ForwardDiff.gradient(x -> loss_neuralode(x, i_exp), p)
        grad_norm[i_exp] = norm(grad, 2)
        update!(opt, p, grad)
    end
    for i_exp in 1:n_exp
        loss_epoch[i_exp] = loss_neuralode(p, i_exp)
    end
    loss_train = mean(loss_epoch[1:n_exp_train]);
    loss_val = mean(loss_epoch[n_exp_train + 1:end]);
    set_description(epochs, string(@sprintf("Loss train %.2e val %.2e gnorm %.1e lr %.1e",
                                             loss_train, loss_val, mean(grad_norm), opt[1].eta)))
    cb(p, loss_train, loss_val);
end

for i_exp in 1:n_exp
    cbi(p, i_exp)
end


text:
# based on NLPModels

export is_unconstrained

# Base type for an optimization model.
abstract type AbstractOptiModel end

# Base type for metadata related to an optimization model.
abstract type AbstractOptiModelMeta end

struct OptiModelMeta <: AbstractOptiModelMeta

  # A composite type that represents the main features of
  # the optimization problem
  #
  #  optimize     f(x) + g(x)
  #  subject to   c(x) in S
  #
  # where x         is an nvar-dimensional vector,
  #       f         is a differentiable real-valued objective function,
  #       g         is a proximable extended-real-valued objective function,
  #       c         is a vector-valued differentiable constraint function,
  #       S         is an ncon-dimensional closed, projectable set,
  #       optimize  is either "minimize" or "maximize".

  nvar     ::Int          # number of variables
  ncon     ::Int          # number of constraints
  x0       ::Vector       # initial primal guess
  y0       ::Vector       # initial dual guess
  minimize ::Bool         # true if optimize == minimize
  name     ::String       # problem name

  function OptiModelMeta( nvar, ncon; x0=zeros(nvar,), y0=zeros(ncon,), minimize=true, name="Generic" )
    if (nvar < 1) || (ncon < 0)
      error("Nonsensical dimensions")
    end

    @lencheck nvar x0
    @lencheck ncon y0

    new( nvar, ncon, x0, y0, minimize, name )
  end
end

"""
    is_unconstrained(prob)
    is_unconstrained(meta)
"""
is_unconstrained(meta::OptiModelMeta) = meta.ncon == 0
is_unconstrained(prob::AbstractOptiModel) = is_unconstrained( prob.meta )


text:
### THIS IS LOADED INTO the Genie module

import Logging
import REPL, REPL.Terminals

### PUBLIC ###

"""
$TYPEDSIGNATURES

Creates a new `controller` file. If `pluralize` is `false`, the name of the controller is not automatically pluralized.
"""
function newcontroller(controller_name::Union{String,Symbol}; path::Union{String,Nothing} = nothing, pluralize::Bool = true) :: Nothing
  Generator.newcontroller(string(controller_name), path = path, pluralize = pluralize)
  load_resources()

  nothing
end

"""
$TYPEDSIGNATURES

Creates all the files associated with a new resource.
If `pluralize` is `false`, the name of the resource is not automatically pluralized.
"""
function newresource(resource_name::Union{String,Symbol}; path::String = ".", pluralize::Bool = true, context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  Generator.newresource(string(resource_name), path = path, pluralize = pluralize)

  try
    Core.eval(context, :(SearchLight.Generator.newresource(uppercasefirst($resource_name))))
  catch
  end

  load_resources()

  nothing
end

"""
$TYPEDSIGNATURES

Creates a new Genie `Task` file.
"""
function newtask(task_name::Union{String,Symbol}) :: Nothing
  task_name = string(task_name)
  endswith(task_name, "Task") || (task_name = task_name * "Task")
  Toolbox.new(task_name)

  nothing
end

### PRIVATE ###

"""
$TYPEDSIGNATURES

Recursively adds subfolders of `lib/` to LOAD_PATH.
The `lib/` folder, if present, is designed to host user code in the form of .jl files.
This function loads user code into the Genie app.
"""
function load_libs(root_dir::String = Genie.config.path_lib) :: Nothing
  isdir(root_dir) || return nothing

  push!(LOAD_PATH, root_dir)

  @async for (root, dirs, files) in walkdir(root_dir)
    Threads.@threads for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, p)
    end
  end

  nothing
end

"""
$TYPEDSIGNATURES

Recursively adds subfolders of `resources/` to LOAD_PATH.
"""
function load_resources(root_dir::String = Genie.config.path_resources) :: Nothing
  isdir(root_dir) || return nothing

  push!(LOAD_PATH, root_dir)

  @async for (root, dirs, files) in walkdir(root_dir)
    Threads.@threads for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, joinpath(root, dir))
    end
  end

  nothing
end

"""
$TYPEDSIGNATURES

Recursively adds subfolders of `helpers/` to LOAD_PATH.
"""
function load_helpers(root_dir::String = Genie.config.path_helpers) :: Nothing
  isdir(root_dir) || return nothing

  push!(LOAD_PATH, root_dir)

  @async for (root, dirs, files) in walkdir(root_dir)
    Threads.@threads for dir in dirs
      p = joinpath(root, dir)
      in(p, LOAD_PATH) || push!(LOAD_PATH, joinpath(root, dir))
    end
  end

  nothing
end

"""
$TYPEDSIGNATURES

Loads (includes) the framework's configuration files into the app's module `context`.
The files are set up with `Revise` to be automatically reloaded.
"""
function load_configurations(root_dir::String = Genie.config.path_config; context::Union{Module,Nothing} = nothing) :: Nothing
  secrets_path = joinpath(root_dir, Genie.SECRETS_FILE_NAME)
  isfile(secrets_path) && Revise.includet(default_context(context), secrets_path)

  # check that the secrets_path has called Genie.secret_token!
  if isempty(Genie.secret_token(false)) # do not generate a temporary token in this check
    match_deprecated = isfile(secrets_path) ? match(r"SECRET_TOKEN\s*=\s*\"(.*)\"", readline(secrets_path)) : nothing
    if match_deprecated !== nothing # does the file use the deprecated syntax?
      Genie.secret_token!(match_deprecated.captures[1]) # resolve the issue for now
      @warn "
        $(secrets_path) is using a deprecated syntax to set the secret token.
        Call Genie.Generator.migrate_secrets_file() to resolve this warning.
      "
    end

    Genie.secret_token() # emits a warning and re-generates the token if secrets_path is not valid
  end

  nothing
end

"""
$TYPEDSIGNATURES

Loads (includes) the framework's initializers.
The files are set up with `Revise` to be automatically reloaded.
"""
function load_initializers(root_dir::String = Genie.config.path_config; context::Union{Module,Nothing} = nothing) :: Nothing
  dir = joinpath(root_dir, Genie.config.initializers_folder)

  isdir(dir) || return nothing

  Threads.@threads for i in readdir(dir)
    fi = joinpath(dir, i)
    endswith(fi, ".jl") && Revise.includet(default_context(context), fi)
  end

  nothing
end

"""
$TYPEDSIGNATURES

Loads (includes) the framework's plugins initializers.
"""
function load_plugins(root_dir::String = Genie.config.path_plugins; context::Union{Module,Nothing} = nothing) :: Nothing
  isdir(root_dir) || return nothing

  Threads.@threads for i in readdir(root_dir)
    fi = joinpath(root_dir, i)
    endswith(fi, ".jl") && Revise.includet(default_context(context), fi)
  end

  nothing
end

"""
$TYPEDSIGNATURES

Loads the routes file.
"""
function load_routes_definitions(routes_file::String = Genie.ROUTES_FILE_NAME; context::Union{Module,Nothing} = nothing, additional_routes_file::String = Genie.APP_FILE_NAME) :: Nothing
  isfile(routes_file) && Revise.includet(default_context(context), routes_file)

  nothing
end

const SECRET_TOKEN = Ref{String}("") # global state

"""
$TYPEDSIGNATURES

Return the secret token used in the app for encryption and salting.

Usually, this token is defined through `Genie.secret_token!` in the `config/secrets.jl` file.
Here, a temporary one is generated for the current session if no other token is defined and
`generate_if_missing` is true.
"""
function secret_token(generate_if_missing::Bool = true; context::Union{Module,Nothing} = nothing)
  if context != nothing
    @warn "secret_token not context-dependent any more; the context argument is deprecated"
  end

  if isempty(SECRET_TOKEN[]) && generate_if_missing
    @warn "
          No secret token is defined through `Genie.secret_token!(\"token\")`. Such a token
          is needed to hash and to encrypt/decrypt sensitive data in Genie, including cookie
          and session data.

          If your app relies on cookies or sessions make sure you generate a valid token,
          otherwise the encrypted data will become unreadable between app restarts.

          You can resolve this issue by generating a valid `config/secrets.jl` file with a
          random token, calling `Genie.Generator.write_secrets_file()`.
          "
    secret_token!()
  end

  SECRET_TOKEN[]
end

"""
$TYPEDSIGNATURES

Define the secret token used in the app for encryption and salting.
"""
function secret_token!(value::AbstractString = Generator.secret_token())
  SECRET_TOKEN[] = value

  value
end

"""
$TYPEDSIGNATURES

Sets the module in which the code is loaded (the app's module)
"""
function default_context(context::Union{Module,Nothing} = nothing)
  try
    context === nothing ? Main.UserApp : context
  catch ex
    @error ex
    @__MODULE__
  end
end

"""
$TYPEDSIGNATURES

Main entry point to loading a Genie app.
"""
function load(; context::Union{Module,Nothing} = nothing) :: Nothing
  context = default_context(context)

  Genie.Configuration.isdev() && Core.eval(context, :(__revise_mode__ = :eval))

  App.bootstrap(context)

  t = Terminals.TTYTerminal("", stdin, stdout, stderr)

  load_configurations(context = context)

  replprint("initializers", t, clearline = 0, prefix = "Loading ")
  load_initializers(context = context)

  replprint("helpers", t, prefix = "Loading ")
  load_helpers()

  replprint("lib", t, prefix = "Loading ")
  load_libs()

  replprint("resources", t, prefix = "Loading ")
  load_resources()

  replprint("plugins", t, prefix = "Loading ")
  load_plugins(context = context)

  replprint("routes", t, prefix = "Loading ")
  load_routes_definitions(context = context)

  replprint("\nReady! \n", t, clearline = 2, color = :green, bold = :true)
  println()

  nothing
end

"""
$TYPEDSIGNATURES

Prints app loading progress to the console.
"""
function replprint(output::String, terminal;
                    newline::Int = 0, clearline::Int = 1, color::Symbol = :white, bold::Bool = false, sleep_time::Float64 = 0.2,
                    prefix::String = "", prefix_color::Symbol = :green, prefix_bold::Bool = true)
  if clearline > 0
    for i in 1:(clearline+1)
      REPL.Terminals.clear_line(terminal)
    end
  end

  isempty(prefix) || printstyled(prefix, color = prefix_color, bold = prefix_bold)
  printstyled(output, color = color, bold = bold)

  if newline > 0
    for i in 1:newline
      println()
    end
  end

  sleep(sleep_time)
end


text:
module Common
using Base.Test

include("test_comparison.jl")
include("test_obtain.jl")
include("test_constructor.jl")
include("test_clone.jl")

export test_constructor
export test_obtain
export test_comnparison
export test_clone

end