text:
function allConsequentAreSet(productionRule::ProductionRule)
    check = true
    for i=1:length(productionRule._consequent)
        if productionRule._consequent[i]._depth == -Inf
            check = false
            break
        end
    end
    return check
end

function setFourthDepth!(gramm::Grammar)

    function isRecursive(symbol::UInt16)

        rule = getProductionRule(gramm, symbol)
        recursive = false
        for i=1:length(rule._consequent)
            if rule._consequent[i]._recursive
                recursive = true
            else
                recursiveInternal = false
                if rule._consequent[i]._depth != 0
                    for j=1:length(rule._consequent[i]._symbols)
                        if !isTerminalSymbol(gramm, rule._consequent[i]._symbols[j])
                            if isRecursive(rule._consequent[i]._symbols[j])
                                recursive = true
                                recursiveInternal = true
                            end
                        end
                    end
                    rule._consequent[i]._recursive = recursiveInternal
                end
            end
        end
        return recursive
    end

    for i=1:length(gramm._P)
        productionRule = gramm._P[i]
        isRecursive(productionRule._antecedent)
    end

end

"""
    Profundidad para aquellos consecuentes que sean Inf, en este punto ya deberia conocer cual es
    su profundidad.
"""

function setThirdDepth!(gramm::Grammar)

    for i=1:length(gramm._P)
        productionRule = gramm._P[i]
        for j=1:length(productionRule._consequent)
            antecedent = productionRule._antecedent
            consequent = productionRule._consequent[j]
            depth = 0;
            count = 0;

            if consequent._depth == 32000
                max = 0
                sum = 0
                for g=1:length(consequent._symbols)
                    symbol = consequent._symbols[g]
                    if !isTerminalSymbol(gramm, symbol)
                        current = getDepth(gramm, symbol)
                        sum = sum + current
                        if current > max
                            max = current
                        end
                    end
                end
                consequent._minProductions = sum
                consequent._depth = max
            end
        end
    end

end

"""
    Esta profundidad es para aquellos consecuentes que sean -Inf (que no esten seteados),
    y para setear los antecedentes que tenga como consecuente Inf y algun valor establecido
    de profundidad en otros, cogiendo el minimo de estos otros como profundidad.
"""

function setSecondDepth!(gramm::Grammar)

    count = length(gramm._P)
    nIter = 0
    nMaxIter = getNRules(gramm) + getNProductionRules(gramm)

    while (count != 0 && nIter < nMaxIter)
        count = 0
        nIter = nIter + 1
        for i=1:length(gramm._P)
            productionRule = gramm._P[i]
            nInf = 0
            if !allConsequentAreSet(productionRule) || !hasDepth(productionRule)
                nConsequentSet = 0
                min = 32000

                for j=1:length(productionRule._consequent)
                    antecedent = productionRule._antecedent
                    consequent = productionRule._consequent[j]
                    if consequent._depth == 32000
                        nInf = nInf + 1
                        nConsequentSet = nConsequentSet + 1
                    elseif consequent._depth == -32000
                        max = -32000
                        sum = 0
                        for g=1:length(consequent._symbols)
                            symbol = consequent._symbols[g]
                            if !isTerminalSymbol(gramm, symbol)

                                current = getDepth(gramm, symbol)
                                sum = sum + current
                                if current == -32000
                                    max = -32000
                                    break
                                end
                                if current > max
                                    max = current
                                end
                            end
                        end
                        if sum != -32000
                            consequent._minProductions = sum
                            consequent._depth = max
                        end
                    else
                        nConsequentSet = nConsequentSet + 1
                        current = consequent._depth
                        if current < min
                            min = current
                        end
                    end

                end
                if nInf == length(productionRule._consequent)
                    error("The production rule of <$(productionRule._antecedent)> doesn't lead to any terminal or it appears in all of its consequents")
                end
                if nConsequentSet == length(productionRule._consequent)
                    productionRule._depth = min + 1
                else
                    count = count + 1
                end
            end
        end
    end
    if nIter == nMaxIter
        error("The grammar has a loop. For example: <ex> -> <var> ; <var> -> <ex>")
    end
end

""" Setea la primera profundidad es decir, aquellos consecuentes que sean un no terminal se setean a 0 profundidad,
    y las reglas de produccion que los contienen a 1 de profundidad. Por otro lado, aquellas reglas con recursividad
    directa se setean como recursivas e Inf.
"""
function setFirstDepth!(gramm::Grammar)

    for i=1:length(gramm._P)
        productionRule = gramm._P[i]
        for j=1:length(productionRule._consequent)
            antecedent = productionRule._antecedent
            consequent = productionRule._consequent[j]
            depth = 0;
            count = 0;
            for g=1:length(consequent._symbols)
                symbol = consequent._symbols[g]

                if !isTerminalSymbol(gramm, symbol)
                    if symbol == antecedent
                        consequent._depth = 32000
                        consequent._recursive = true
                        break
                    end
                else
                    count = count + 1
                end
            end
            if count == length(consequent._symbols)
                consequent._depth = 0
                consequent._minProductions = 0
                productionRule._depth = 1
            end
        end
    end
end

function setDepth!(gramm::Grammar)
    setFirstDepth!(gramm)
    setSecondDepth!(gramm)
    setThirdDepth!(gramm)
    setFourthDepth!(gramm)
end

function reorderNonTerminals!(N::Array{String}, productionRules::Array{ProductionRule})
    len = length(N)

    newN = similar(Array{String}, len)
    reference = similar(Array{UInt16}, len)

    for i in eachindex(productionRules)
        reference[i] = productionRules[i]._antecedent
    end
    aux = [(reference[i], N[i]) for i=1:len]
    aux = sort(aux)

    for i in eachindex(productionRules)
        productionRules[i]._antecedent = i
        N[i] = aux[i][2]
    end

    for i in eachindex(productionRules)
        productionRule = productionRules[i]
        for j in eachindex(productionRule._consequent)
            antecedent = productionRule._antecedent
            consequent = productionRule._consequent[j]
            depth = 0;
            count = 0;
            for g in eachindex(consequent._symbols)
                if consequent._symbols[g] < 1000
                    consequent._symbols[g] = findfirst(x->consequent._symbols[g]==x, reference)
                end
            end
        end
    end

end

function readSymbol(str::String)
    nullChar(element) = findnext(x->x==element, [' ','\n','\r','\t'], 1) != nothing
    count = 0
    len = length(str)
    i = 2
    symbol = ""
    while (str[i]!='>')
        if nullChar(str[i])
            error("There shouldn't be neither whitespaces nor line breaks inside a <symbol>")
        else
            symbol=symbol*str[i]
            i=i+1
        end
        if i > len
            error("While reading a symbol, the string ends prematurely")
        end
    end
    if i == 1
        error("Symbol must have, at least, one character")
    end

    return symbol, i
end

function createGrammar(N::Array{String}, T::Array{String}, P::Array{String}, S::String)
    if findnext(x->x==S,N,1) == nothing
        error("The initial symbol S ($S) must be in the non-terminal set, N")
    end

    nullChar(element) = findnext(x->x==element, [' ','\n','\r','\t'], 1) != nothing
    function findNextNotNullChar(str, i)
        index = findnext(!nullChar, str, i)
        index == nothing ? nothing : collect(index)[1]
    end
    function findNextNull(str, i)
        index = findnext(nullChar, str, i)
        index == nothing ? nothing : collect(index)[1]
    end

    productionRules = Array{ProductionRule}(undef, 0)

    for productionRule in P
        consequents = Array{Consequent}(undef,0)
        symbol = ""
        antecedent = false

        len = length(productionRule)
        antecedentSymbol = ""
        symbols = Array{UInt16}(undef, 0)

        i = 1
        while (i<=len)
            if !nullChar(productionRule[i])
                # If the antecedent is not found yet, the first letter must be '<'
                if !antecedent
                    if productionRule[i] != '<'
                        error("The first character of the production rule must be '<' which represents the begining of the antecedent symbol")
                    else
                        antecedentSymbol, sum = readSymbol(productionRule[i:end])
                        symbolIndex = findnext(x->x==antecedentSymbol, N, 1)

                        if symbolIndex == nothing
                            error("Nonterminal symbol ($symbol) is in production rules, but it isn't in Non-terminal set, N")
                        end

                        antecedentSymbol = symbolIndex

                        i = i + sum
                        i = findNextNotNullChar(productionRule, i)
                        if productionRule[i] != '-'
                            error("The character expected after the antecedent symbol is '-' ")
                        end
                        if productionRule[i+1] != '>'
                            error("The character expected after the antecedent symbol and '-' is '>'")
                        end
                        i = i+2
                        antecedent = true
                    end

                # Inside the consequent of the rule
                else
                    if productionRule[i] == '<'
                        symbol, sum = readSymbol(productionRule[i:end])
                        symbolIndex = findnext(x->x==symbol, N, 1)
                        if symbolIndex == nothing
                            error("Non-terminal symbol ($symbol) is in production rules, but it isn't in Non-terminal set, N")
                        end

                        push!(symbols, symbolIndex)
                        i = i + sum - 1
                        if findnext(x->!nullChar(x), productionRule, i+1) == nothing
                            push!(consequents, Consequent(symbols))
                            break
                        end
                    elseif productionRule[i] == '|'
                        if length(symbols) == 0
                            error("A '|' found without any symbol before")
                        end

                        push!(consequents, Consequent(symbols))
                        symbols = Array{UInt16}(undef, 0)
                        if productionRule[i:end] == "| <var> "
                            activate = true
                        end
                    else

                        next = findnext(x->nullChar(x) || x=='<', productionRule, i)
                        if next == nothing
                            symbol = productionRule[i:end]
                            symbolIndex = findnext(x->x==symbol, T, 1)
                            if symbolIndex == nothing
                                error("Terminal symbol ($symbol) is in production rules, but it isn't in Terminal set, T")
                            end
                            push!(symbols, symbolIndex+1000)
                            push!(consequents, Consequent(symbols))
                            break
                        else
                            next = collect(next)[1]
                            symbol = productionRule[i:next-1]
                            if symbol == "SP"
                                if findfirst(" ", T) == nothing
                                    push!(T, " ")
                                end
                                symbol = " "
                            end
                            symbolIndex = findnext(x->x==symbol, T, 1)
                            if symbolIndex == nothing
                                error("Terminal symbol ($symbol) is in production rules, but it isn't in Terminal set, T")
                            end
                            push!(symbols, symbolIndex+1000)
                            i = next-1
                        end
                    end
                end
            end
            i = i+1
        end
        if antecedent == -1 || antecedent == 0
            error("The antecendent must be followed by '->' and there should be an antecedent")
        end

        push!(productionRules, ProductionRule(antecedentSymbol, consequents))
    end
    for i in eachindex(N)
        antecedents = [x._antecedent for x in productionRules]
        if findnext(x->x==i, antecedents, 1) == nothing
            error("Non-terminal symbol ($(N[i])) isn't an antecedent of any production rule")
        end
    end

    reorderNonTerminals!(N, productionRules)
    gramm = Grammar(N, T, productionRules, S)

    setDepth!(gramm)

    return gramm
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using Test
import LokiLogger, HTTP, Sockets, JSON3, Logging, URIs

const ch = Channel{Any}(Inf)
function fivexx(http, c=500)
    HTTP.setstatus(http, c)
    HTTP.setheader(http, "Content-Length" => "0")
    HTTP.startwrite(http)
    HTTP.close(http.stream)
end
const server = Sockets.listen(Sockets.InetAddr(Sockets.ip"0.0.0.0", 9100))
tsk = @async HTTP.listen("0.0.0.0", 9100; server=server) do http
    if http.message.method != "POST"
        return fivexx(http)
    end
    ctype = HTTP.header(http.message, "Content-Type")
    if ctype != "application/json"
        return fivexx(http)
    end
    bytes = read(http)
    clength = HTTP.header(http.message, "Content-Length")
    if length(bytes) != parse(Int, clength)
        return fivexx(http)
    end
    try
        put!(ch, JSON3.read(bytes))
    catch
        return fivexx(http)
    end
    # All good
    HTTP.setstatus(http, 204)
    HTTP.startwrite(http)
end
atexit(() -> (close(server); try fetch(tsk) catch e end))
const lokiserver = "http://localhost:9100"
const starttime = time() * 1e9

@testset "LokiLogger.jl" begin
    # Defaults
    logger = LokiLogger.Logger(lokiserver)
    Logging.with_logger(logger) do
        @info "hello, world"
        @error "hello, world" error("hello, error")
    end
    json = take!(ch)
    @test keys(json) == Set([:streams])
    @test length(json.streams) == 1
    @test keys(json.streams[1]) == Set([:stream, :values])
    stream = json.streams[1].stream
    @test stream.host == Sockets.gethostname()
    @test stream.app == "LokiLogger.jl"
    values = json.streams[1].values
    @test length(values) == 1
    @test length(values[1]) == 2
    @test starttime < parse(Int, values[1][1]) < time() * 1e9
    msg = values[1][2]
    @test occursin("level=info msg=\"hello, world\" module=", msg)
    # Error capturing
    json = take!(ch)
    msg = json.streams[1].values[1][2]
    @test occursin("level=error msg=\"Exception while generating log record in module", msg)
    # JSON format
    logger = LokiLogger.Logger(LokiLogger.json, lokiserver; labels=Dict("region" => "eu-central"))
    Logging.with_logger(logger) do
        @error "hello, world" _group=:group
    end
    line = (@__LINE__) - 2
    json = take!(ch)
    @test keys(json) == Set([:streams])
    @test length(json.streams) == 1
    @test keys(json.streams[1]) == Set([:stream, :values])
    stream = json.streams[1].stream
    @test stream.region == "eu-central"
    values = json.streams[1].values
    @test length(values) == 1
    @test length(values[1]) == 2
    @test starttime < parse(Int, values[1][1]) < time() * 1e9
    msg = values[1][2]
    @test occursin("{\"level\":\"error\",\"msg\":\"hello, world\",\"module\":\"", msg)
    msgjson = JSON3.read(msg)
    @test keys(msgjson) == Set([:level, :msg, :module, :file, :line, :group, :id])
    @test msgjson.level == "error"
    @test msgjson.msg == "hello, world"
    @test msgjson.module == "Main"
    @test msgjson.file == @__FILE__
    @test msgjson.line == line
    @test msgjson.group == "group"
    # Custom format
    logger = LokiLogger.Logger(URIs.URI(lokiserver)) do io, args
        println(io, args.level, "| ", args.message)
    end
    Logging.with_logger(logger) do
        @debug "hello, world"
    end
    json = take!(ch)
    msg = json.streams[1].values[1][2]
    @test msg == "Debug| hello, world"
end

ext: jl


text:
using KernelAbstractions
using CUDA
# One can write the weights as a matrix (μⱼ/(μⱼ + μⱼ)); needs to be checked!
# qp_matrix = qp_μ ./ (qp_μ .+ qp_μ');

# qp_matrix2 = (1 ./ qp_μ .+ 1 ./ qp_μ')
# Random tau:
# τ_nSpec = 4 .+ randn(1000, )

# and exp(1/μᵢ + 1/μⱼ)
# mat2   = repeat(qp_matrix2, outer=[1, 1, size(τ_nSpec, 1)])
# tau_3d = repeat(τ_nSpec, inner=size(qp_matrix, 1), outer=size(qp_matrix, 1))

@kernel function get_r!(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w, ndoubl)
    i, j, n = @index(Global, NTuple)
    # @show i, j, n
    r⁻⁺[i,j,n] = ϖ[n] * Z⁻⁺[i,j] * (μ[j] / (μ[i] + μ[j])) * (1 - exp(-dτ[n] * ((1 / μ[i]) + (1 / μ[j])))) * w[j]
    if μ[i] == μ[j]
        # 𝐓⁺⁺(μᵢ, μᵢ) = (exp{-τ/μᵢ} + ϖ ̇𝐙⁺⁺(μᵢ, μᵢ) ̇(τ/μᵢ) ̇exp{-τ/μᵢ}) ̇𝑤ᵢ
        if i == j
            t⁺⁺[i,j,n] = exp(-dτ[n] / μ[i]) + ϖ[n] * Z⁺⁺[i,i] * (dτ[n] / μ[i]) * exp(-dτ[n] / μ[i]) * w[i]
        else
            t⁺⁺[i,j,n] = ϖ[n] * Z⁺⁺[i,i] * (dτ[n] / μ[i]) * exp(-dτ[n] / μ[i]) .* w[i]
        end
    else
        # 𝐓⁺⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁺⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ-μⱼ)) ̇(exp{-τ/μᵢ} - exp{-τ/μⱼ}) ̇𝑤ⱼ
        # (𝑖 ≠ 𝑗)
        t⁺⁺[i,j,n] = ϖ[n] * Z⁺⁺[i,j] .* (μ[j] / (μ[i] - μ[j])) * (exp(-dτ[n] / μ[i]) - exp.(-dτ[n] / μ[j])) * w[j]
    end
end

n = 40
nSpec = 10000
r⁻⁺_CPU, t⁺⁺_CPU = rand(n, n, nSpec), rand(n, n, nSpec);
ϖ, dτ    = rand(nSpec), rand(nSpec);
μ, w     = rand(n, n), rand(n, n);
Z⁻⁺, Z⁺⁺ = rand(n, n), rand(n, n);

device = KernelAbstractions.CPU()
kernel! = get_r!(device)
kernel!(r⁻⁺_CPU, t⁺⁺_CPU, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w,0, ndrange=size(r⁻⁺_CPU));

function test2(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w)
    event = kernel!(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w, 0, ndrange=size(r⁻⁺));
    wait(device, event)
    synchronize();
end
@time test2(r⁻⁺_CPU, t⁺⁺_CPU, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w)
@time test2(r⁻⁺_CPU, t⁺⁺_CPU, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w)

r⁻⁺, t⁺⁺ = CuArray(rand(n, n, nSpec)), CuArray(rand(n, n, nSpec));
ϖ, dτ    = CuArray(rand(nSpec)), CuArray(rand(nSpec));
μ, w     = CuArray(rand(n, n)), CuArray(rand(n, n));
Z⁻⁺, Z⁺⁺ = CuArray(rand(n, n)), CuArray(rand(n, n));

device = KernelAbstractions.CUDADevice()
kernel! = get_r!(device)
kernel!(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w,0, ndrange=size(r⁻⁺));
function test2(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w)
    event = kernel!(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w, 0, ndrange=size(r⁻⁺));
    wait(device, event)
    synchronize();
end
@time test2(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w)
@time test2(r⁻⁺, t⁺⁺, ϖ, dτ, Z⁻⁺, Z⁺⁺, μ, w)


text:
module BinanceBot

using Binance, CoinMarketCap, DataFrames, CSV

export write_symbol_map, read_symbol_map, write_market_cap_list, read_market_cap_list,
    write_price_data, read_price_data, sorted_roundtrip

function write_symbol_map(file::AbstractString="symbolmap.csv", client=Client())
    df = get_symbol_map(DataFrame, client)
    CSV.write(file, df)
    df
end
function read_symbol_map(file::AbstractString="symbolmap.csv")
    get_symbol_map(Dict, nothing, CSV.read(file))  #Dict{String, Tuple{String,String}}
end
function write_market_cap_list(file::AbstractString="marketcaplist.csv",
                               market=Market();
                               thresh::Float64=0.5e9)
    df = ticker_by_cap(DataFrame, market, thresh)
    CSV.write(file, df[[:symbol]])
end
function read_market_cap_list(file::AbstractString="marketcaplist.csv")
    d = CSV.read(file)
    convert(Vector, d[:symbol])  #Vector{String}
end
function write_price_data(file::AbstractString="prices.csv", client=Client())
    df = get_all_tickers(DataFrame, client)
    CSV.write(file, df)
end
function read_price_data(file::AbstractString="prices.csv")
    get_all_tickers(Dict, nothing, CSV.read(file))
end

function sorted_roundtrip(client, sym_list::Vector{String},
                          prices::Dict{String,Float64}=get_all_tickers(Dict, client),
                          base_syms::Vector{String}=base_symbols(),
                          derived_syms::Vector{String}=setdiff(sym_list,base_symbols());
                          cost_per_trade::Float64=0.001)
    slip = 1.0 - cost_per_trade
    v = []
    for b1 in base_syms, d in derived_syms, b2 in base_syms
        if b1 != b2
            if !haskey(prices, "$d$b1") || !haskey(prices, "$d$b2")
                continue
            end
            if  !haskey(prices, "$b1$b2") && !haskey(prices, "$b2$b1")
                continue
            end
            p1 = prices["$d$b1"] * slip
            p2 = 1 / prices["$d$b2"] * slip
            if haskey(prices, "$b1$b2")
                p3 = prices["$b1$b2"] * slip
            elseif haskey(prices, "$b2$b1")
                p3 = 1 / prices["$b2$b1"] * slip
            else
                error("shouldn't have gotten here... $b2$b1")
            end
            p = p1 * p2 * p3
            push!(v, (b1, d, b2, p1, p2, p3, p))
        end
    end
    sort!(v, by=x->x[7], rev=true)
end

end # module


text:
"""
    solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, ::ExtendedMethod, δ::Real, vertex_list;...)

Solve NORBiP using the extended formulation.
Creates a disjunctive constraint over vertices of each dual polyhedron using indicator constraints.

`vertex_list` is a vector of vectors of length `bp.mu`, with each
individual sub-vector containing a list of tuples `(α, β)` with α a vector of `bp.ml` elements and β a scalar.
"""
function solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, ::ExtendedMethod, δ::Real, vertex_list; verbose=false, resvar=false)
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds) = bilevel_optimality(m, bp, upperlevel=false, resvar=resvar)
    # early termination if relaxed or bilevel optimistic infeasible
    optimize!(m)
    st = JuMP.termination_status(m)
    if st != MOI.OPTIMAL && st != MOI.DUAL_INFEASIBLE
        return (m, st)
    end
    w = Vector{Vector{VariableRef}}(undef, bp.mu)
    for k in 1:bp.mu
        w[k] = @variable(m, [l=1:length(vertex_list[k])], Bin, base_name="w_$(k)_")
        @constraint(m, sum(w[k]) >= 1)
        for l in eachindex(vertex_list[k])
            (α, β) = vertex_list[k][l]
            first_part = if resvar
                dot(α, m[:r])
            else
                dot(α, bp.b - bp.A * x)
            end
            @constraint(m, [w[k][l], first_part + β * dot(bp.d, v) + dot(bp.G[k,:], x)]
                in MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(bp.q[k] - β * δ))
            )
        end
    end
    return (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds, w)
end

function bilevel_optimality(m::JuMP.Model, bp::BilevelProblem; upperlevel=true, resvar=false)
    (_, x, v, λ, σ, s, upperfeas, lowerfeas, kkt) = high_point_relaxation(m, bp, upperlevel=upperlevel, resvar=resvar)
    @constraint(m, kkt2_bounds[i=1:bp.ml],
        [λ[i], s[i]] in MOI.SOS1([0.4,0.6])
    )
    @constraint(m, kkt2_var[j=1:bp.nl],
        [v[j], σ[j]] in MOI.SOS1([0.4,0.6])
    )
    return (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds)
end

"""
    high_point_relaxation(m::JuMP.Model, bp::BilevelProblem)

Build the relaxed high-point relaxation model, along with dual feasibility (no complementarity constraint).
This problem is purely linear & continuous.
"""
function high_point_relaxation(m::JuMP.Model, bp::BilevelProblem; upperlevel = true, resvar=false)
    @variable(m, x[1:bp.nu] >= 0)
    @variable(m, v[1:bp.nl] >= 0)
    @variable(m, σ[1:bp.nl] >= 0)
    @variable(m, λ[1:bp.ml] >= 0)
    @variable(m, s[1:bp.ml] >= 0)
    if upperlevel
        @constraint(m, bp.G * x + bp.H * v .<= bp.q)
    end
    if resvar
        @variable(m, r[i=1:bp.ml])
        @constraint(m, bp.b - bp.A * x .== r)
        @constraint(m, lowerfeas, bp.B * v  .+ s .== r)
    else
        @constraint(m, lowerfeas, bp.A * x + bp.B * v  .+ s .== bp.b)
    end
    @constraint(m, kkt, bp.d .+ bp.B' * λ .- σ .== 0)
    return (m, x, v, λ, σ, s, nothing, lowerfeas, kkt)
end

"""
Attach the upper-level objective to the model
"""
function add_upper_level(m::JuMP.Model, bp::BilevelProblem, x, v)
    @objective(m, Min, dot(bp.cx, x) + dot(bp.cy, v))
    return nothing
end

"""
    solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, ::LazyExtended, δ::Real; optimizer, poly_lib, verbose=false)
"""
function solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, ::LazyExtended, δ::Real; optimizer, poly_lib, verbose=false)
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds) = bilevel_optimality(m, bp, upperlevel=true)
    add_upper_level(m, bp, x, v)
    # early termination if relaxed or bilevel optimistic infeasible
    optimize!(m)
    st = JuMP.termination_status(m)
    w = [VariableRef[] for _ in 1:bp.mu]
    problem_ref = (model = m, x = x, v = v, λ = λ, σ = σ, s = s, w = w)
    if st != MOI.OPTIMAL && st != MOI.DUAL_INFEASIBLE
        return (Symbol(st), problem_ref)
    end
    unexplored_subproblems = BitSet(1:bp.mu)
    valid_subproblems = BitSet()
    res = bp.b - bp.A * x
    while length(unexplored_subproblems) > 0
        optimize!(m)
        empty!(valid_subproblems)
        k = pop!(unexplored_subproblems)
        subresult = dual_polyhedron(bp.B, bp.d, bp.H, k, x, v, bp.A, bp.G, bp.q, bp.b, δ; poly_lib=poly_lib, optimizer=optimizer)
        if subresult[1] == :INFEASIBLE
            return (:INFEASIBLE,)
        elseif subresult[1] == :OPTIMAL
            push!(valid_subproblems, k)
        else subresult[1] == :VERTICES
            vertex_list = subresult[2]
            w[k] = @variable(m, [l=1:length(vertex_list)], Bin, base_name = "w_$k")
            @constraint(m, sum(w[k]) >= 1)
            for l in eachindex(vertex_list)
                (α, β) = vertex_list[l]
                @constraint(m, [w[k][l], dot(α, res) + β * dot(bp.d, v) + dot(bp.G[k,:], x)]
                    in MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(bp.q[k] - β * δ))
                )
            end
            optimize!(m)
            if st != MOI.OPTIMAL && st != MOI.DUAL_INFEASIBLE
                return (Symbol(st), problem_ref)
            end
            # put the optimal back in the subproblem pool
            for v in valid_subproblems
                push!(unexplored_subproblems, v)
            end
        end
    end
    return (:OPTIMAL, problem_ref)
end

"""
    solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, ::LazyBatched, δ::Real; optimizer, poly_lib)
"""
function solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, lb::LazyBatched, δ::Real; optimizer, poly_lib, verbose=false)
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds) = bilevel_optimality(m, bp, upperlevel=true)
    add_upper_level(m, bp, x, v)
    # early termination if relaxed or bilevel optimistic infeasible
    optimize!(m)
    st = JuMP.termination_status(m)
    w = [VariableRef[] for _ in 1:bp.mu]
    problem_ref = (model = m, x = x, v = v, λ = λ, σ = σ, s = s, w = w)
    if st != MOI.OPTIMAL && st != MOI.DUAL_INFEASIBLE
        return (Symbol(st), problem_ref)
    end
    unexplored_subproblems = BitSet(1:bp.mu)
    optimize!(m)
    valid_subproblems = BitSet()
    current_batch = BitSet()
    Γ = lb.batch_size
    while length(unexplored_subproblems) > 0
        verbose && @info "new batch, $(length(unexplored_subproblems)) unexplored subproblems"
        empty!(current_batch)
        while length(unexplored_subproblems) > 0 && length(current_batch) < Γ
            k = pop!(unexplored_subproblems)
            verbose && @info "subproblem $k"
            optimize!(m)
            subresult = dual_polyhedron(bp.B, bp.d, bp.H, k, x, v, bp.A, bp.G, bp.q, bp.b, δ; poly_lib=poly_lib, optimizer=optimizer)
            if subresult[1] == :INFEASIBLE
                verbose && @info "infeasible subproblem"
                return (:INFEASIBLE,)
            elseif subresult[1] == :OPTIMAL
                push!(valid_subproblems, k)
                verbose && @info "$k optimal"
            else subresult[1] == :VERTICES
                push!(current_batch, k)
                vertex_list = subresult[2]
                @info "Adding $(length(vertex_list)) vertices at subproblem $k"
                w[k] = @variable(m, [l=eachindex(vertex_list)], Bin, base_name = "w_$k")
                @constraint(m, sum(w[k]) >= 1)
                res = bp.b - bp.A * x
                for l in eachindex(vertex_list)
                    (α, β) = vertex_list[l]
                    @constraint(m, [w[k][l], dot(α, res) + β * dot(bp.d, v) + dot(bp.G[k,:], x)]
                        in MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(bp.q[k] - β * δ))
                    )
                end
                for vsp in valid_subproblems
                    push!(unexplored_subproblems, vsp)
                end
                empty!(valid_subproblems)
            end
        end
        optimize!(m)
        st = JuMP.termination_status(m)
        if st != MOI.OPTIMAL && st != MOI.DUAL_INFEASIBLE
            return (Symbol(st), problem_ref)
        end
    end
    return (:OPTIMAL, problem_ref)
end

"""
    solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, ::BilinearTerms, δ::Real; verbose=false)
"""
function solve_near_optimal(m::JuMP.Model, bp::BilevelProblem, ::BilinearTerms, δ::Real; verbose=false)
    (m, x, v, λ, σ, s, upperfeas, lowerfeas, kkt, kkt2_var, kkt2_bounds) = bilevel_optimality(m, bp, upperlevel=true)
    add_upper_level(m, bp, x, v)
    @variable(m, α[1:bp.mu, 1:bp.ml] >= 0)
    @variable(m, β[1:bp.mu] >= 0)
    @constraint(m, dual_polyhedron_feasible[k=1:bp.mu],
        bp.B' * α[k,:] .+ β[k] * bp.d .≥ bp.H[k,:]
    )
    @constraint(m, bilinear_cons[k=1:bp.mu],
        α[k,:] ⋅ (bp.b - bp.A * x) .+ β[k] * (bp.d ⋅ v + δ) .≤ bp.q[k] - bp.G[k,:] ⋅ x
    )
    return (model = m, x = x, v = v, λ = λ, σ = σ, s = s, α = α, β = β)
end


text:
using ClinicalTrialSampleSize
using Base.Test

# test for power calculations
@test_throws ErrorException power(OneSampleMean(0, 9), n = 1, std = 1)
@test_throws ErrorException power(OneSampleMean(0, 9), std = 1)
power(TwoSampleProp(0.9, 0.5, 1), std = 1, n = 20, alpha = 0.05)
@test_throws ErrorException power(TwoSampleMean(0.9, 0.5, 1), n = 20, std = -1,  alpha = 0.05)
@test_throws ErrorException power(TwoSampleMean(0.9, 0.5, 1), n = 20, alpha = 0.05)
@test_throws ErrorException power(TwoSampleMean(0.9, 0.5, 1), n = 20, std = (-0.02, -0.03), alpha = 0.05)
@test_throws ErrorException power(OneSampleMean(0, 9), n = 25, std = 1, alpha = 1.5)
@test_throws ErrorException power(OneSampleMean(0, 9), n = 10, std = 1, alpha = 0.05, side = "greater")

# test for sample size calculations
@test_throws ErrorException sample_size(OneSampleMean(0, 9), power = 1.2, std = 1)
@test_throws ErrorException sample_size(OneSampleMean(0, 9), std = 1)
sample_size(TwoSampleProp(0.9, 0.5, 1), std = 1, power = 0.8, alpha = 0.05)
@test_throws ErrorException  sample_size(TwoSampleMean(0.9, 0.5, 1), power = 0.8, alpha = 0.05)
@test_throws ErrorException sample_size(TwoSampleMean(0, 9, 1), power = 0.9, std = (-2, -3), alpha = 0.05)
@test_throws ErrorException sample_size(TwoSampleMean(0, 9, 1), power = 0.9, std = -1, alpha = 0.5)
@test_throws ErrorException sample_size(TwoSampleMean(0, 9, 1), power = 0.9, std = (2, 3), alpha = 1.5)
@test_throws ErrorException sample_size(OneSampleProp(0.2, 0.5), power = 0.8, std = 1, alpha = 0.05, side = "greater")

return true


text:
# using Revise
using POMDPs
using BasicPOMCP
using POMCPOW
using VOOTreeSearch
using ParticleFilters
using CPUTime
using DataFrames
using ARDESPOT
using Distributions
using VDPTag2
using Dates
using CSV
using POMDPSimulators
using VDPTag2Experiments
using LinearAlgebra: diagm, PosDefException, det, eigvals, eigvecs
using Random: MersenneTwister, seed!
using Distances: euclidean

using MCTS
using StaticArrays
using POMDPModelTools: GenerativeBeliefMDP
include("./vdp_next_action.jl")

max_time = 0.1
max_depth = 10

@show start_mean = [110.0, 30.0, 30.0, 5.0, 100.0, 0.8]
@show start_cov = diagm([100.0^2, 10.0^2, 10.0^2, 10.0^2, 20.0^2, 0.1^2])
d = MvNormal(start_mean, start_cov)
on_cluster = false
save_data = false
parallelize = true

if on_cluster
  max_query = 100_000
  K = 150 # number of parameter samples
  n = 20 # number of evaluation simulations
  m = 30 # number of elite samples
  max_iters = 50
else
  max_query = 100
  K = 15 # number of parameter samples
  n = 1 # number of evaluation simulations
  m = 10 # number of elite samples
  max_iters = 3
end

# VDPTag2
@show pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 2.8)))

println(stderr, "Max time: ", max_time)
println(stderr, "Max depth: ", max_depth)
println(stderr, "Max query: ", max_query)
println(stderr, "Start mean: ", start_mean)
println(stderr, "Start cov: ", start_cov)
println(stderr, "POMDP: ", pomdp)

function gen_sims(x::Vector{Float64}, n, k, seed)
    c = x[1]
    @assert c >= 0.0
    k_act = x[2]
    @assert k_act >= 1.0
    inv_alpha_act = x[3]
    @assert inv_alpha_act >= 0.1
    k_obs = x[4]
    @assert k_obs >= 1.0
    inv_alpha_obs = x[5]
    @assert inv_alpha_obs >= 0.1
    p_voo = x[6]
    @assert p_voo >= 0.1
    @assert p_voo <= 1.0

    sims = []

    for i in 1:n
        rng = MersenneTwister(7)

        sig = diagm([0.1])
        epsilon = 0.1/10.0
        voo_opt = VOOActionGenerator(p_voo, Distributions.Uniform(), MvNormal(sig), rng, epsilon, true, 20)

        m = 10
        node_updater = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                           LowVarianceResampler(m),
                                           0.05, m, rng)
        ro = ToNextML(mdp(pomdp), rng)
        solver = DPWSolver(n_iterations=max_query,
                           exploration_constant=c,
                           depth=max_depth,
                           max_time=max_time,
                           k_action = k_act,
                           alpha_action = 1/inv_alpha_act,
                           k_state = k_obs,
                           alpha_state = 1/inv_alpha_obs,
                           estimate_value=RolloutEstimator(ro),
                           next_action=RootToNextMLFirstVOO(voo_opt),
                           rng=rng
                          )
        belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), node_updater)
        planner = solve(solver, belief_mdp)
        filter = SIRParticleFilter(deepcopy(pomdp),
                                   100_000,
                                   rng=MersenneTwister(i+90_000))

        seed!(planner, i+40000*k)
        sim = Sim(deepcopy(pomdp),
                  planner,
                  filter,
                  rng=MersenneTwister(i+50_000*k),
                  max_steps=50,
                  metadata=Dict(:i=>i, :k=>k)
                 )
        push!(sims, sim)
    end

    return sims
end

datestring = Dates.format(now(), "e_d_u_Y_HH_MM")
results_csv = DataFrame(c=Float64[], k_act=Float64[], inv_alpha_act=Float64[], k_obs=Float64[], inv_alpha_obs=Float64[], p=Float64[])
try
  for i in 1:max_iters
      sims = []
      params = Vector{Float64}[]
      println(stderr,"creating $K simulation sets")
      for k in 1:K
          p = rand(d)
          p[1] = max(0.0, p[1])
          p[2] = max(1.0, p[2])
          p[3] = max(0.1, p[3])
          p[4] = max(1.0, p[4])
          p[5] = max(0.1, p[5])
          p[6] = max(0.1, min(1.0, p[6]))
          push!(params, p)
          k_sims = gen_sims(p, n, k, i)
          append!(sims, k_sims)
      end
      @assert length(params) == K
      if parallelize
        results = run_parallel(sims)
      else
        results = run(sims)
      end
      combined = by(results, :k) do df
          DataFrame(mean_reward=mean(df[:reward]))
      end
      order = sortperm(combined[:mean_reward])
      elite = params[combined[:k][order[K-m:end]]]
      elite_matrix = Matrix{Float64}(undef, length(start_mean), m)
      for k in 1:m
          elite_matrix[:,k] = elite[k]
      end
      try
          global d = fit(typeof(d), elite_matrix)
      catch ex
          if ex isa PosDefException
              println(stderr,"pos def exception")
              global d = fit(typeof(d), elite_matrix += 0.01*randn(size(elite_matrix)))
          else
              rethrow(ex)
          end
      end
      println(stderr, "Iteration $i")
      println(stderr, "Reward (mean): ", mean(combined[:mean_reward]))
      println(stderr, "Mean: ", mean(d))
      println(stderr, "Cov (det): ", det(cov(d)))
      ev = eigvals(cov(d))
      println(stderr, "Cov (eig): ", ev)
      for j in 1:length(ev)
          println(stderr, "Eigvecs: ", eigvecs(cov(d))[:,j])
      end
      push!(results_csv, mean(d))
      if save_data
        fname_csv = joinpath(@__DIR__, "../data", "VDP_pft_vpw_opt"*"_"*datestring*".csv")
        CSV.write(fname_csv, results_csv)
        println(stderr,"Saving parameters to ", fname_csv)
      end
  end
finally
  Distributed.rmprocs(worker_ids)
end


text:
using Base.Test

module test_heatcapacity_resources
import WaterData

func = exp
exponential = WaterData.TFuncHeatCapacity(func)
func2d(x, y) = exp(x) + 2exp(y)
exponential2d = WaterData.PTFuncHeatCapacity(func2d)
constant = WaterData.ConstantHeatCapacity(1)
end

@testset "Heat capacity handling" begin
    res = test_heatcapacity_resources

    @testset "Values are as expected" begin
        cₚ_exp = res.exponential
        cₚ_exp2 = res.exponential2d
        cₚ_const = res.constant

        @testset "called directly" begin
            # constant values are always equal
            @test cₚ_const(0) == cₚ_const(10)
            @test cₚ_const(0, 10) == 1

            # calling a single-input heat capacity
            @test cₚ_exp(1) == e^1
            @test cₚ_exp(0) == e^0
            @test cₚ_exp(1, 2) == e^2

            # calling a temp- and pressure-dependent heat capacity
            @test cₚ_exp2(0, 1) == e^0 + 2e^1
            @test cₚ_exp2(1, 1) == e^1 + 2e^1
            @test cₚ_exp2(1, 3) == e^1 + 2e^3
            @test cₚ_exp2(2, 1) == e^2 + 2e^1

            # needs both temperature and pressure
            @test_throws MethodError cₚ_exp2(1)
        end
    end
end


text:
# These are some functions to allow static type-checking of Julia programs

module TypeCheck
export checkreturntypes, checklooptypes, checkmethodcalls,
  methodswithdescendants

## Modifying functions from Base

# return the type-inferred AST for each method of a generic function
function Base.code_typed(f::Function)
  Expr[code_typed(m) for m in f.env]
end

# return the type-inferred AST for one method of a generic function
function Base.code_typed(m::Method)
 linfo = m.func.code
 (tree,ty) = Base.typeinf(linfo,m.sig,())
 if !isa(tree,Expr)
     ccall(:jl_uncompress_ast, Any, (Any,Any), linfo, tree)
 else
    tree
 end
end

# specify a method via a method, function, or function+argument-type-tuple
# prints out the local variables of that method, and their types
function Base.whos(f,args...)
  for e in code_typed(f,args...)
    display(MethodSignature(e))
    for x in sort(e.args[2][2];by=x->x[1])
      println("\t",x[1],"\t",x[2])
    end
    println("")
  end
end

## Basic Operations on Function Exprs

# given an Expr representing a method, return its inferred return type
returntype(e::Expr) =  e.args[3].typ

# given an Expr representing a method, return an Array of Exprs representing its body
body(e::Expr) = e.args[3].args

# given an Expr representing a method, return all of the return statement in its body
returns(e::Expr) = filter(x-> typeof(x) == Expr && x.head==:return,body(e))

# given an Expr representing a method,
# return all function all Exprs contained in return statements in the method body
function extractcallsfromreturns(e::Expr)
  rs = returns(e)
  rs_with_calls = filter(x->typeof(x.args[1]) == Expr && x.args[1].head == :call,rs)
  Expr[expr.args[1] for expr in rs_with_calls]
end

# A type that covers both Types and TypeVars
AType = Union(Type,TypeVar)

# given an Expr representing a method, get the type of each of the arguments in the signature
function argumenttypes(e::Expr)
  argnames = Symbol[typeof(x) == Symbol ? x : x.args[1] for x in e.args[1]]
  argtuples = filter(x->x[1] in argnames, e.args[2][2]) #only arguments, no local vars
  AType[t[2] for t in argtuples]
end

# given an Expr, determine if it is calling a TopNode
# (this affects how we should handle resolving the callee name)
istop(e) = Base.is_expr(e,:call) && typeof(e.args[1]) == TopNode

# given a call Expr (:call, :call1, :new), determine its return type
function returntype(e::Expr,context::Expr) #must be :call,:new,:call1
  if Base.is_expr(e,:new); return e.typ; end
  if Base.is_expr(e,:call1) && isa(e.args[1], TopNode); return e.typ; end
  if !Base.is_expr(e,:call); error("Expected :call Expr"); end

  if istop(e)
    return e.typ
  end

  callee = e.args[1]
  if istop(callee)
    return returntype(callee,context)
  elseif isa(callee,SymbolNode) # only seen (func::F), so non-generic function
    return Any
  elseif is(callee,Symbol)
    if e.typ != Any || any([isa(x,LambdaStaticData) for x in e.args[2:end]])
      return e.typ
    end

    if isdefined(Base,callee)
      f = eval(Base,callee)
      if !isa(f,Function) || !isgeneric(f)
        return e.typ
      end
      fargtypes = tuple([argumenttype(ea,context) for ea in e.args[2:end]])
      return Union([returntype(ef) for ef in code_typed(f,fargtypes)]...)
    else
      return @show e.typ
    end
  end

  return e.typ
end

# for an Expr `e` used as an argument to a call in function Expr `context`,
# determine the type of `e`
function argumenttype(e::Expr,context::Expr)
 if Base.is_expr(e,:call) || Base.is_expr(e,:new) || Base.is_expr(e,:call1)
   return returntype(e,context)
 end

 @show e
 return Any
end

# for a Symbol `s` used as an argument to a call in a function Expr `e`,
# determine the type of `s`.
function argumenttype(s::Symbol,e::Expr)
  vartypes = [x[1] => x[2] for x in e.args[2][2]]
  s in vartypes ? (vartypes[@show s]) : Any
end

# as above, but for different call argument types
argumenttype(s::SymbolNode,e::Expr) = s.typ
argumenttype(t::TopNode,e::Expr) = Any
argumenttype(l::LambdaStaticData,e::Expr) = Function
argumenttype(q::QuoteNode,e::Expr) = argumenttype(q.value,e)

# as above, but for various literal values
argumenttype(n::Number,e::Expr) = typeof(n)
argumenttype(c::Char,e::Expr) = typeof(c)
argumenttype(s::String,e::Expr) = typeof(s)
argumenttype(i,e::Expr) = typeof(i) #catch all, hopefully for more literals

# start, next, and done are the functions for-loops use to iterate
# having implementations of them makes a type iterable
# this defines iterating over a DataType to mean iterating over its descendants
# (breadth-first-search ordering)
Base.start(t::DataType) = [t]
function Base.next(t::DataType,arr::Vector{DataType})
  c = pop!(arr)
  append!(arr,[x for x in subtypes(c)])
  (c,arr)
end
Base.done(t::DataType,arr::Vector{DataType}) = length(arr) == 0

# this is like the `methodswith` function from Base,
# but it considers all descendants of the type
# rather than the type itself (or the type + super types)
# it produces a list of function names and statistics
# about how many descendants implement that function
# this is good for discovering implicit interfaces
function methodswithdescendants(t::DataType;onlyleaves::Bool=false,lim::Int=10)
  d = Dict{Symbol,Int}()
  count = 0
  for s in t
    if !onlyleaves || (onlyleaves && isleaftype(s))
      count += 1
      fs = Set{Symbol}()
      for m in methodswith(s)
        push!(fs,m.func.code.name)
      end
      for sym in fs
        d[sym] = get(d,sym,0) + 1
      end
    end
  end
  l = [(k,v/count) for (k,v) in d]
  sort!(l,by=(x->x[2]),rev=true)
  l[1:min(lim,end)]
end

# check all the generic functions in a module
function checkallmodule(m::Module;test=checkreturntypes,kwargs...)
  score = 0
  for n in names(m)
    f = eval(m,n)
    if isgeneric(f) && typeof(f) == Function
      fm = test(f;mod=m,kwargs...)
      score += length(fm.methods)
      display(fm)
    end
  end
  println("The total number of failed methods in $m is $score")
end

# use checkallmodule to implement the Module version of other checks
checkreturntypes(m::Module;kwargs...) = checkallmodule(m;test=checkreturntypes,kwargs...)
checklooptypes(m::Module) = checkallmodule(m;test=checklooptypes)
checkmethodcalls(m::Module) = checkallmodule(m;test=checkmethodcalls)

## Checking that return values are base only on input *types*, not values.

type MethodSignature
  typs::Vector{AType}
  returntype::Union(Type,TypeVar) # v0.2 has TypeVars as returntypes; v0.3 does not
end
MethodSignature(e::Expr) = MethodSignature(argumenttypes(e),returntype(e))
function Base.writemime(io, ::MIME"text/plain", x::MethodSignature)
  println(io,"(",join([string(t) for t in x.typs],","),")::",x.returntype)
end

type FunctionSignature
  methods::Vector{MethodSignature}
  name::Symbol
end

function Base.writemime(io, ::MIME"text/plain", x::FunctionSignature)
  for m in x.methods
    print(io,string(x.name))
    display(m)
  end
end

# given a function, run checkreturntypes on each method
function checkreturntypes(f::Function;kwargs...)
  results = MethodSignature[]
  for e in code_typed(f)
    (ms,b) = checkreturntype(e;kwargs...)
    if b push!(results,ms) end
  end
  FunctionSignature(results,f.env.name)
end

# given an Expr representing a Method,
# determine whether its return type is based
# only on the arugment types or whether it is
# also influenced by argument values
# (the Method fails the check if the return type depends on values)
function checkreturntype(e::Expr;kwargs...)
  (typ,b) = isreturnbasedonvalues(e;kwargs...)
  (MethodSignature(argumenttypes(e),typ),b)
end

# Determine whether this method's return type might change based on input values rather than input types
function isreturnbasedonvalues(e::Expr;mod=Base)
  rt = returntype(e)
  ts = argumenttypes(e)
  if isleaftype(rt) || rt == None return (rt,false) end

  for t in ts
    if !isleaftype(t)
      return (rt,false)
    end
  end

  cs = [returntype(c,e) for c in extractcallsfromreturns(e)]
  for c in cs
    if rt == c
       return (rt,false)
    end
  end

  return (rt,true) # return is not concrete type; all args are concrete types
end

## Checking that variables in loops have concrete types

type LoopResult
  msig::MethodSignature
  lines::Vector{(Symbol,Type)} #TODO should this be a specialized type? SymbolNode?
  LoopResult(ms::MethodSignature,ls::Vector{(Symbol,Type)}) = new(ms,unique(ls))
end

function Base.writemime(io, ::MIME"text/plain", x::LoopResult)
  display(x.msig)
  for (s,t) in x.lines
    println(io,"\t",string(s),"::",string(t))
  end
end

type LoopResults
  name::Symbol
  methods::Vector{LoopResult}
end

function Base.writemime(io, ::MIME"text/plain", x::LoopResults)
  for lr in x.methods
    print(io,string(x.name))
    display(lr)
  end
end

# for a given Function, run checklooptypes on each Method
function checklooptypes(f::Function;kwargs...)
  lrs = LoopResult[]
  for e in code_typed(f)
    lr = checklooptypes(e)
    if length(lr.lines) > 0 push!(lrs,lr) end
  end
  LoopResults(f.env.name,lrs)
end

# for an Expr representing a Method,
# check that the type of each variable used in a loop
# has a concrete type
checklooptypes(e::Expr;kwargs...) = LoopResult(MethodSignature(e),loosetypes(loopcontents(e)))

# This is a function for trying to detect loops in the body of a Method
# Returns lines that are inside one or more loops
function loopcontents(e::Expr)
  b = body(e)
  loops = Int[]
  nesting = 0
  lines = {}
  for i in 1:length(b)
    if typeof(b[i]) == LabelNode
      l = b[i].label
      jumpback = findnext(
        x-> (typeof(x) == GotoNode && x.label == l) || (Base.is_expr(x,:gotoifnot) && x.args[end] == l),
        b, i)
      if jumpback != 0
        push!(loops,jumpback)
        nesting += 1
      end
    end
    if nesting > 0
      push!(lines,(i,b[i]))
    end

    if typeof(b[i]) == GotoNode && in(i,loops)
      splice!(loops,findfirst(loops,i))
      nesting -= 1
    end
  end
  lines
end

# given `lr`, a Vector of expressions (Expr + literals, etc)
# try to find all occurances of a variables in `lr`
# and determine their types
# `method` is only used as part of constructing the return type
function loosetypes(lr::Vector)
  lines = (Symbol,Type)[]
  for (i,e) in lr
    if typeof(e) == Expr
      es = copy(e.args)
      while !isempty(es)
        e1 = pop!(es)
        if typeof(e1) == Expr
          append!(es,e1.args)
        elseif typeof(e1) == SymbolNode && !isleaftype(e1.typ) && typeof(e1.typ) == UnionType
          push!(lines,(e1.name,e1.typ))
        end
      end
    end
  end
  return lines
end

## Check method calls

type CallSignature
  name::Symbol
  argumenttypes::Vector{AType}
end
function Base.writemime(io, ::MIME"text/plain", x::CallSignature)
  println(io,string(x.name),"(",join([string(t) for t in x.argumenttypes],","),")")
end

type MethodCalls
  m::MethodSignature
  calls::Vector{CallSignature}
end

function Base.writemime(io, ::MIME"text/plain", x::MethodCalls)
  display(x.m)
  for c in x.calls
    print(io,"\t")
    display(c)
  end
end

type FunctionCalls
  name::Symbol
  methods::Vector{MethodCalls}
end

function Base.writemime(io, ::MIME"text/plain", x::FunctionCalls)
  for mc in x.methods
    print(io,string(x.name))
    display(mc)
  end
end

# given a Function, run `checkmethodcalls` on each Method
# and collect the results into a FunctionCalls
function checkmethodcalls(f::Function;kwargs...)
  calls = MethodCalls[]
  for m in f.env
    e = code_typed(m)
    mc = checkmethodcalls(e,m;kwargs...)
    if !isempty(mc.calls)
      push!(calls, mc)
    end
  end
  FunctionCalls(f.env.name,calls)
end

# given an Expr representing a Method,
# and the Method it represents,
# check the Method body for calls to non-existant Methods
function checkmethodcalls(e::Expr,m::Method;kwargs...)
  if Base.arg_decl_parts(m)[3] == symbol("deprecated.jl")
    CallSignature[]
  end
  nomethoderrors(e,methodcalls(e);kwargs...)
end

# Find any methods that match the given CallSignature
function hasmatches(mod::Module,cs::CallSignature)
  if isdefined(mod,cs.name)
    f = eval(mod,cs.name)
    if isgeneric(f)
      opts = methods(f,tuple(cs.argumenttypes...))
      if isempty(opts)
        return false
      end
    end
  else
    #println("$mod.$(cs.name) is undefined")
  end
  return true
end

# Find any CallSignatures that indicate potential NoMethodErrors
function nomethoderrors(e::Expr,cs::Vector{CallSignature};mod=Base)
  output = CallSignature[]
  for callsig in cs
    if !hasmatches(mod,callsig)
      push!(output,callsig)
    end
  end
  MethodCalls(MethodSignature(e),output)
end

# Look through the body of the function for `:call`s
function methodcalls(e::Expr)
  b = body(e)
  lines = CallSignature[]
  for s in b
    if typeof(s) == Expr
      if s.head == :return
        append!(b, s.args)
      elseif s.head == :call
        if typeof(s.args[1]) == Symbol
          push!(lines,CallSignature(s.args[1], [argumenttype(e1,e) for e1 in s.args[2:end]]))
        end
      end
    end
  end
  lines
end

end  #end module


text:
export atomdofs, variablecell, fixedcell, dofmgr_resetref!, variablecell!,
       fixedcell!, set_clamp!, set_free!, set_mask!, reset_clamp!,
       inplane!

using SparseArrays: SparseMatrixCSC, nnz, sparse, findnz

using LinearAlgebra: rmul!, det

# ========================================================================
#          AUXILIARY FUNCTIONS
# ========================================================================

function zeros_free(n::Integer, x::AbstractVector, free)
   z = zeros(eltype(x), n)
   z[free] = x
   return z
end

function insert_free!(p::AbstractArray, x::AbstractVector, free)
   p[free] = x
   return p
end

# a helper function to get a valid positions array from a dof-vector
positions(at::AbstractAtoms,
          ifree::AbstractVector{<:Integer},
          dofs::Dofs) =
   insert_free!(positions(at) |> mat, dofs, ifree) |> vecs

"""
`_pos_to_dof` : a helper function that will convert a positions-based
block-hessian into a classical dof-based hessian with the standard JuLIP
ordering of dofs.
"""
function _pos_to_dof(Hpos::SparseMatrixCSC, at::AbstractAtoms{T}) where {T}
   I, J, Z = Int[], Int[], T[]
   for C in (I, J, Z); sizehint!(C, 9 * nnz(Hpos)); end
   Nat = length(at)
   @assert Nat == size(Hpos, 2)
   # TODO: this findnz creates an extra copy of all data, which we should avoid
   for (iat, jat, zat) in zip(findnz(Hpos)...)
      for a = 1:3, b = 1:3
         push!(I, 3 * (iat-1) + a)
         push!(J, 3 * (jat-1) + b)
         push!(Z, zat[a,b])
      end
   end
   return sparse(I, J, Z, 3*Nat, 3*Nat)
end

# TODO: this is a temporary hack, and I think we need to
#       figure out how to do this for more general constraints
#       maybe not too terrible
projectxfree(at::Atoms, A::SparseMatrixCSC) =
      A[at.dofmgr.xfree, at.dofmgr.xfree]

"""
`analyze_mask` : helper function to generate list of dof indices from
lists of atom indices indicating free and clamped atoms
"""
function analyze_mask(at, free, clamp, mask)
   if length(findall((free != nothing, clamp != nothing, mask != nothing))) > 1
      error("FixedCell: only one of `free`, `clamp`, `mask` may be provided")
   elseif all( (free == nothing, clamp == nothing, mask == nothing) )
      # in this case (default) all atoms are free
      return collect(1:3*length(at))
   end
   # determine free dof indices
   Nat = length(at)
   if clamp != nothing
      # revert to setting free
      free = setdiff(1:Nat, clamp)
   end
   if free != nothing
      # revert to setting mask
      mask = fill(false, 3, Nat)
      if !isempty(free)
         mask[:, free] .= true
      end
   end
   return findall(mask[:])
end

analyze_mask(at, free::Colon, clamp::Nothing, mask::Nothing) =
      collect(1:3*length(at))

analyze_mask(at, free::Nothing, clamp::Colon, mask::Nothing) =
      Int[]

# ========================================================================
#          MAIN DOF MANAGER INTERFACE
# ========================================================================

"""
`VariableCell`: both atom positions and cell shape are free;

**WARNING:** (1) before manipulating the dof-vectors returned by a `VariableCell`
constraint, read *meaning of dofs* instructions at bottom of help text!

(2) The `volume` field is a signed volume.

Constructor:
```julia
VariableCell(at::AbstractAtoms; free=..., clamp=..., mask=..., fixvolume=false)
```
Set at most one of the kwargs:
* no kwarg: all atoms are free
* `free` : list of free atom indices (not dof indices)
* `clamp` : list of clamped atom indices (not dof indices)
* `mask` : 3 x N Bool array to specify individual coordinates to be clamped

### Meaning of dofs

On call to the constructor, `VariableCell` stored positions and deformation
`X0, F0`, dofs are understood *relative* to this initial configuration.

`dofs(at, cons::VariableCell)` returns a vector that represents a pair
`(Y, F1)` of a displacement and a deformation matrix. These are to be understood
*relative* to the reference `X0, F0` stored in `cons` as follows:
* `F = F1`   (the cell is then `F'`)
* `X = [F1 * (F0 \\ y)  for y in Y)]`

One aspect of this definition is that clamped atom positions still change via
`F`
"""

DofManager(at::AbstractAtoms) = DofManager(length(at), eltype(at))

DofManager(Nat::Integer, T = Float64) =
      DofManager( false,                      # variablecell
                  collect(1:3*Nat),    # xfree
                  LinearConstraint{T}[],      # lincons
                  JVec{T}[],                  # X0
                  zero(JMat{T}) )             # F0

import Base.==
==(d1::DofManager, d2::DofManager) =
      (d1.variablecell == d2.variablecell) && (d1.xfree == d2.xfree)
      # TODO: potential need to add equality about the constraints?
      #       is xfree even relevant for equality?

variablecell(at::Atoms) = at.dofmgr.variablecell
fixedcell(at::Atoms) = !variablecell(at::Atoms)

function dofmgr_resetref!(at)
   at.dofmgr.X0 = copy(positions(at))
   at.dofmgr.F0 = cell(at)'
   return at
end

# convenience function to return DoFs associated with a particular atom
atomdofs(at::Atoms, i::Integer) =
      findall(in(3*i-2:3*i), at.dofmgr.xfree)

"""
`variablecell!(at)` : make the cell-shape variable
"""
variablecell!(at) = set_variablecell!(at, true)

"""
`fixedcell!(at)` : freeze the cell shape
"""
fixedcell!(at) = set_variablecell!(at, false)

"""
`set_variablecell!` : specify whether cell shape is variable or fixed
"""
function set_variablecell!(at::AbstractAtoms, tf::Bool)
   at.dofmgr.variablecell = tf
   if tf
      dofmgr_resetref!(at)
   end
   return at
end

set_clamp!(at::Atoms, Iclamp::AbstractVector{<: Integer}) =
      set_clamp!(at; clamp = Iclamp)

set_free!(at::Atoms, Ifree::AbstractVector{<: Integer}) =
      set_clamp!(at; free = Ifree)

set_mask!(at::Atoms, mask) =
      set_clamp!(at; mask=mask)

reset_clamp!(at) = set_clamp!(at; free = :)

function set_clamp!(at::Atoms; free=nothing, clamp=nothing, mask=nothing)
   at.dofmgr.xfree = analyze_mask(at, free, clamp, mask)
   return at
end

# reverse map:
#   F -> F
#   X[n] = F * F^{-1} X0[n]

function position_dofs(at::Atoms)
   if fixedcell(at)
      return mat(at.X)[at.dofmgr.xfree]
   end
   # variable cell case:
   X = positions(at)
   F = cell(at)'
   A = at.dofmgr.F0 * inv(F)
   U = [A * x for x in X]   # switch to broadcast!
   return [mat(U)[at.dofmgr.xfree]; Matrix(F)[:]]
end

posdofs(x) = x[1:end-9]

celldofs(x) = x[end-8:end]

function set_position_dofs!(at::AbstractAtoms{T}, x::Dofs) where {T}
   if fixedcell(at)
      return set_positions!(at, positions(at, at.dofmgr.xfree, x))
   end
   # variable cell case:
   F = JMat{T}(celldofs(x))
   A = F * inv(at.dofmgr.F0)
   Y = copy(at.dofmgr.X0)
   mat(Y)[at.dofmgr.xfree] = posdofs(x)
   for n = 1:length(Y)
      Y[n] = A * Y[n]
   end
   set_positions!(at, Y)
   set_cell!(at, F')
   return at
end

function momentum_dofs(at::AbstractAtoms)
   if fixedcell(at)
      return mat(momenta(at))[at.dofmgr.xfree]
   end
   @error("`momentum_dofs` is not yet implmement for variable cells")
   return nothing
end

function set_momentum_dofs!(at::AbstractAtoms, p::Dofs)
   if fixedcell(at)
      return set_momenta!(at, zeros_free(3 * length(at), p, at.dofmgr.xfree) |> vecs)
   end
   @error("`set_momentum_dofs!` is not yet implmement for variable cells")
   return nothing
end

# for a variation x^t_i = (F+tU) F_0^{-1} (u_i + t v_i)
#       ~ U F^{-1} F F0^{-1} u_i + F F0^{-1} v_i
# we get
#      dE/dt |_{t=0} = U : (S F^{-T}) - < (F * inv(F0))' * frc, v>
#
# this is nice because there is no contribution from the stress to
# the positions component of the gradient

"""
`sigvol` : signed volume
"""
sigvol(C::AbstractMatrix) = det(C)
sigvol(at::AbstractAtoms) = sigvol(cell(at))

"""
`sigvol_d` : derivative of signed volume
"""
sigvol_d(C::AbstractMatrix) = sigvol(C) * inv(C)'
sigvol_d(at::AbstractAtoms) = sigvol_d(cell(at))

function gradient(calc::AbstractCalculator, at::Atoms)
   if fixedcell(at)
      return rmul!(mat(forces(at))[at.dofmgr.xfree], -1.0)
   end
   F = cell(at)'
   A = F * inv(at.dofmgr.F0)
   G = forces(at)
   for n = 1:length(G)
      G[n] = - A' * G[n]
   end
   S = - virial(at) * inv(F)'                  # ∂E / ∂F
   # S += at.dofmgr.pressure * sigvol_d(at)'     # applied stress  TODO: revive this!
   return [ mat(G)[at.dofmgr.xfree]; Array(S)[:] ]
end

function hessian(calc::AbstractCalculator, at::AbstractAtoms)
   if fixedcell(at)
      H =  _pos_to_dof(hessian_pos(calculator(at), at), at)
      return projectxfree(at, H)
   end
   @error("`hessian` is not yet implmement for variable cells")
   return nothing
end

function inplane!(at::Atoms; free = 1:length(at), clamp = nothing, i1 = 1, i2 = 2)
   if clamp != nothing
      free = setdiff(1:length(at), clamp)
   end
   mask = fill(false, 3, length(at))
   mask[i1, free] .= true
   mask[i2, free] .= true
   set_mask!(at, mask)
   return at
end

# TODO:
#   - anti-plane
#   - pressure
#   - fixvolume
#   - once we add an external potential we need to think about terminology
#     should energy -> potential_energy?; total_energy a new one?
#     total_energy(calc::AbstractCalculator, at::AbstractAtoms, cons::VariableCell) =
#              energy(calc, at) + cons.pressure * sigvol(at)


text:
#### Deprecate on 0.6 (to be removed on 0.7)

@Base.deprecate expected_logdet meanlogdet

function probs(d::DiscreteUnivariateDistribution)
    Base.depwarn("probs(d::$(typeof(d))) is deprecated. Please use pdf(d) instead.", :probs)
    return pdf(d)
end

function Binomial(n::Real, p::Real)
    Base.depwarn("Binomial(n::Real, p) is deprecated. Please use Binomial(n::Integer, p) instead.", :Binomial)
    Binomial(Int(n), p)
end

function Binomial(n::Real)
    Base.depwarn("Binomial(n::Real) is deprecated. Please use Binomial(n::Integer) instead.", :Binomial)
    Binomial(Int(n))
end

function BetaBinomial(n::Real, α::Real, β::Real)
    Base.depwarn("BetaBinomial(n::Real, α, β) is deprecated. Please use BetaBinomial(n::Integer, α, β) instead.", :BetaBinomial)
    BetaBinomial(Int(n), α, β)
end

# vectorized versions
for fun in [:pdf, :logpdf,
            :cdf, :logcdf,
            :ccdf, :logccdf,
            :invlogcdf, :invlogccdf,
            :quantile, :cquantile]

    _fun! = Symbol('_', fun, '!')
    fun! = Symbol(fun, '!')

    @eval begin
        @deprecate ($_fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray) r .= ($fun).(d, X) false
        @deprecate ($fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray) r .= ($fun).(d, X) false
        @deprecate ($fun)(d::UnivariateDistribution, X::AbstractArray) ($fun).(d, X)
    end
end

@deprecate pdf(d::DiscreteUnivariateDistribution) pdf.(Ref(d), support(d))

# Wishart constructors
@deprecate Wishart(df::Real, S::AbstractPDMat, warn::Bool) Wishart(df, S)
@deprecate Wishart(df::Real, S::Matrix, warn::Bool) Wishart(df, S)
@deprecate Wishart(df::Real, S::Cholesky, warn::Bool) Wishart(df, S)

# Deprecate 3 arguments expectation and once with function in second place
@deprecate expectation(distr::DiscreteUnivariateDistribution, g::Function, epsilon::Real) expectation(g, distr; epsilon=epsilon) false
@deprecate expectation(distr::ContinuousUnivariateDistribution, g::Function, epsilon::Real) expectation(g, distr) false
@deprecate expectation(distr::Union{UnivariateDistribution,MultivariateDistribution}, g::Function; kwargs...) expectation(g, distr; kwargs...) false

# Deprecate `MatrixReshaped`
const MatrixReshaped{S<:ValueSupport,D<:MultivariateDistribution{S}} = ReshapedDistribution{2,S,D}
@deprecate MatrixReshaped(
    d::MultivariateDistribution, n::Integer, p::Integer=n
) reshape(d, (n, p))


text:
module FluxTubes
using ModelingToolkit
using ...Electrical: Pin

@parameters t
D = Differential(t)

export PositiveMagneticPort, NegativeMagneticPort, TwoPort
include("utils.jl")

export Ground, Idle, Short, Crossing, ConstantPermeance, ConstantReluctance, EddyCurrent, ElectroMagneticConverter
include("basic.jl")

export  ConstantMagneticPotentialDifference, ConstantMagneticFlux
include("sources.jl")

end  #module


text:
module TestTMapreduce

using Test
using ThreadingTools
using BenchmarkTools: @benchmark

const MAX_ALLOCS = 120

struct FreeMonoid{T}
    word::Vector{T}
end

function Base.:*(x::FreeMonoid, y::FreeMonoid)
    FreeMonoid(vcat(x.word, y.word))
end
Base.:(==)(x::FreeMonoid, y::FreeMonoid) = x.word == y.word
pure(x) = FreeMonoid([x])
Base.one(::Type{FreeMonoid{T}}) where {T} = FreeMonoid{T}(T[])

@testset "sum, prod, minimum, maximum" begin

    for (f, arr) in [
            (identity, [1]),
            (sin, randn(10)),
            (identity, randn(10)),
            (x -> 2x, rand(Int, 10000)),
        ]
        @test Base.prod(arr) ≈ @inferred tprod(arr)
        @test Base.sum(arr) ≈ @inferred tsum(arr)
        @test Base.minimum(arr) === @inferred tminimum(arr)
        @test Base.maximum(arr) === @inferred tmaximum(arr)

        @test Base.prod(f, arr) ≈ @inferred tprod(f, arr)
        @test Base.sum(f, arr) ≈ @inferred tsum(f, arr)
        @test Base.minimum(f, arr) === @inferred tminimum(f, arr)
        @test Base.maximum(f, arr) === @inferred tmaximum(f, arr)
    end

    # empty
    @test_throws ArgumentError tmaximum(Int[])
    @test_throws ArgumentError tminimum(Int[])
    @test tprod(Int[]) === 1
    @test tsum(Int[]) === 0

    @test tsum([1,2,3], batch_size=100) === 6
    @test tsum(x->2x, [1,2,3], batch_size=100) === 12
    @test_throws ArgumentError tsum([1,2,3], batch_size=-100)

    # performance
    data = randn(10^5)
    for red in [tsum, tprod, tminimum, tmaximum]
        b = @benchmark ($red)($data) samples=1 evals=1
        @test b.allocs < MAX_ALLOCS
        b = @benchmark ($red)(sin, $data) samples=1 evals=1
        @test b.allocs < MAX_ALLOCS
    end
end

@testset "treduce, tmapreduce" begin
    @test tmapreduce(pure, *, 1:1, init=pure(-1)) == pure(-1) * pure(1)
    @test tmapreduce(pure, *, 1:1, init=pure(-1)) == mapreduce(pure, *, 1:1, init=pure(-1))

    @test tmapreduce(identity, +, Int[], init=4) == 4
    @test tmapreduce(identity, +, Int[]) == 0
    @test treduce(+, Int[]) == 0

    @test treduce(+, [1,2,3], batch_size=3) === 6
    @test tmapreduce(x->2x, +, [1,2,3], batch_size=3) === 12

    setups = [
                  (f=identity,       op=+, srcs=(1:10,),                    init=0),
                  (f=identity,       op=*, srcs=(String[],),                init="Hello"),
                  (f=identity,       op=+, srcs=(1:10,),                    init=21),
                  (f=x->2x,          op=*, srcs=(collect(1:10),),           init=21),
                  (f=x->2x,          op=*, srcs=(rand(Int, 10^5),),         init=rand(Int)),
                  (f=pure,           op=*, srcs=(randn(4),),                init=pure(42.0)),
                  (f=pure,           op=*, srcs=(randn(10^4),),             init=pure(42.0)),
       ]

    # multi arg mapreduce
    if VERSION >= v"1.2-"
        for n in 1:4
            s= (f=pure∘tuple,  op=*, srcs=[Float64[1] for _ in 1:n], init=pure(tuple(1.0:n...)))
            push!(setups, s)
        end
    else
        @warn "Skipping multi arg mapreduce tests, on julia $VERSION"
    end

    for setup in setups

        res_base = @inferred Base.reduce(setup.op, map(setup.f, setup.srcs...))
        res_tt   = @inferred   treduce(setup.op, map(setup.f, setup.srcs...))
        @test res_base == res_tt

        res_base = @inferred Base.reduce(setup.op, map(setup.f, setup.srcs...), init=setup.init)
        res_tt   = @inferred   treduce(setup.op, map(setup.f, setup.srcs...), init=setup.init)
        @test res_base == res_tt

        res_base = @inferred Base.mapreduce(setup.f, setup.op, setup.srcs...)
        res_tt   = @inferred   tmapreduce(setup.f, setup.op, setup.srcs...)
        @test res_base == res_tt

        res_base = @inferred Base.mapreduce(setup.f, setup.op, setup.srcs..., init=setup.init)
        res_tt   = @inferred   tmapreduce(setup.f, setup.op, setup.srcs..., init=setup.init)
        @test res_base == res_tt
    end

    # performance
    data = randn(10^5)
    b = @benchmark tmapreduce(sin, +, $data) samples=1 evals=1
    @test b.allocs < MAX_ALLOCS

    b = @benchmark treduce(+, $data) samples=1 evals=1
    @test b.allocs < MAX_ALLOCS
end

end#module


text:
import RowEchelon
import LinearAlgebra

#=
Homework 7
Section 3.2
Problems 9
=#

A = [1 -1 -3 0 ;
     0 1 5 4 ;
     -1 2 8 5;
     3 -1 -2 3]

detA = LinearAlgebra.det(A)

println("Original Matrix:")
println(A)
println("Determinate:")
println(detA)


text:
tbl_dir_test = joinpath(pkgdir(MortalityTables), "test", "data", "CSV")
@testset "CSV equality" begin
    @testset "CSV and XTbML equality: $id" for id in [17,428,1152,3302]
        xtbml = MortalityTables.readXTbML(joinpath(soa_tbl_dir, "t$id.xml"))
        csv = MortalityTable(CSV.File(joinpath(tbl_dir_test, "t$id.csv"), header=false, silencewarnings=true))

        @test xtbml.ultimate == csv.ultimate
        if typeof(xtbml) <: MortalityTables.SelectUltimateTable
            for issue_age in eachindex(xtbml.select)
                @test xtbml.select[issue_age] == csv.select[issue_age]
            end
        end

    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PersistentVolumeClaim is a user's request for and claim to a persistent volume

    IoK8sApiCoreV1PersistentVolumeClaim(;
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        spec=nothing,
        status=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#metadata
    - spec::IoK8sApiCoreV1PersistentVolumeClaimSpec : Spec defines the desired characteristics of a volume requested by a pod author. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#persistentvolumeclaims
    - status::IoK8sApiCoreV1PersistentVolumeClaimStatus : Status represents the current information/status of a persistent volume claim. Read-only. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#persistentvolumeclaims
"""
mutable struct IoK8sApiCoreV1PersistentVolumeClaim <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PersistentVolumeClaimSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PersistentVolumeClaimStatus } # spec name: status

    function IoK8sApiCoreV1PersistentVolumeClaim(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeClaim, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaim, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaim, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaim, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaim, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeClaim

const _property_map_IoK8sApiCoreV1PersistentVolumeClaim = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiCoreV1PersistentVolumeClaim = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiCoreV1PersistentVolumeClaimSpec", Symbol("status")=>"IoK8sApiCoreV1PersistentVolumeClaimStatus")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeClaim }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeClaim))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeClaim }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeClaim[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeClaim }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeClaim[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeClaim)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PersistentVolumeClaim }, name::Symbol, val)
end


text:
@testset "1D array sampling" begin

d = FindSteadyStates.Domain(low=3.0, high=4.0)
length_of_vector = 1000

vec = rand_vecU(length_of_vector, d)
vec2 = rand_vec(length_of_vector, Uniform(d.low, d.high))

@test length(vec) == length_of_vector
@test minimum(vec) >= d.low
@test maximum(vec) <= d.high

@test length(vec2) == length_of_vector
@test minimum(vec2) >= d.low
@test maximum(vec2) <= d.high

end

@testset "sampling in uniform and log_uniform" begin

mean(vec) = sum(vec) / length(vec)

## Parameters

down, up = 1.,100.
len = 1000
# de struct
de = DEsteady(func=(x) -> x,p=0,u0=ones(len),method=Tsit5())
domain = Domain(down, up)

domains = fill!(Array{Domain}(undef, len), Domain(down,up))
## Sampling

# Create sampler functor
uni_samp = Uniform(a=down, b=up)
log_samp = Log_uniform(a=down, b=up, base=10.)

# Sampling single value
n1 = uni_samp()
n2 = log_samp()

# Sampling a vector
v1 = uni_samp(len)
v2 = log_samp(len)

# sampling with domains

v_2 = rand_vecU(domains)
# Sampling a vector
println("Sampling a vector")
@time v3_1 = rand_vec(len, uni_samp)
@time v3_2 = rand_vec(len, log_samp)

# test uniform sampler with Demeta struct
@time v4 = rand_vecU(de, domain)

## Test
# v1 v2
@test mean(v1) >= mean(v2)

# v3
@test mean(v3_1) >= mean(v3_2)

@test mean(v1) >= mean(v2)
# v4
@test length(v4) == len
@test (maximum(v4) <= domain.high) & (minimum(v4) >= domain.low)

# v_

@testset "Sampling with domains. Input:ODEmeta" begin
    err = 0
    for i in eachindex(v_2)
        num = v_2[i]
        d = domains[i]
        if !(num >= d.low)
            err += 1
        elseif !(num <= d.high)
            err += 1
        else
            continue
        end
    end
    @test err == 0
end

end

@testset "Sampling with samplers" begin

samplers = [Uniform(1,2) for i in 1:100]
vec = rand_vec(samplers)

@test length(vec) == length(samplers)

end

@testset "General method" begin
    len = 10
    general_samplers = [rand for i in 1:len]

    ## Sampling
    v1 = rand_vec(general_samplers)

    v2 = rand_vecU([[1,2], [1,2], [1,2]])

    @test len == length(v1)
end

@testset "Random generator" begin

ParameterGrid([[1,2,3],[1,3,4]])
param_ranges = [[1.,10.],[1.,20.]]
sam = Uniform(0,1, Int)
sam2  = Uniform(0,10)
l = Log_uniform(0.12,13, 10.0,Int)
l2 = Log_uniform(0.12,1,100., Float64)

len = 10

gen1 = ParameterRandom([sam, sam2, l, l2], 100000)

gen2 = ParameterRandom(param_ranges, sam,  len)
gen3 = ParameterRandom(param_ranges, l, len )

@test length(gen2) == length(gen3)
@test gen2.methods[1] == sam(param_ranges[1])

# Plotting
p1= Plots.histogram([i[1] for i in gen1], title="Uniform Integer"); # Uniform(0,1, Int)
p2 =Plots.histogram([i[2] for i in gen1], title="Uniform"); # Uniform(0,10)
p3 = Plots.histogram([i[3] for i in gen1], title="Log Integer"); # Log_uniform(0.12,13, 10.0,Int)
p4 = Plots.histogram([i[4] for i in gen1], title="Log"); # Log_uniform(0.12,1,100., Float64)

Plots.plot(p1,p2,p3,p4)

end


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

using Markdown

# ╔═╡ f906f80a-9eec-40cb-80ba-bdb06a798b1b
md"""
# Estructuras de datos

Una vez que empecemos a trabajar con muchos datos a la vez, será conveniente guardar nuestros datos en estructuras como arreglos o diccionarios (más allá que sólo variables).

Tipos de estructuras de datos que cubrimos
1. Tuplas
2. Diccionarios
3. Arreglos

Como repaso, las tuplas y los arreglos ambos son secuencias ordenadas de elementos (entonces podemos accesarlos por medio de un índice).
Los diccionarios y los arreglos son mutables.

¡Explicaremos más brevemente!
"""

# ╔═╡ 67852a29-1c6c-4ab6-be1e-35a16be10a40
md"""
## Tuplas / Tuples

Podemos crear una tupla encerrando una secuencia ordenada de elementos con `( )`.

Sintaxis:

```julia
(item1, item2, ...)```
"""

# ╔═╡ 5f86014d-54f2-497e-b6d3-529ff3e6f43e
mis_animales_favoritos = ("pingüino", "gato", "petauro_del_azúcar")

# ╔═╡ b21a8b34-0b8d-4c06-b38c-f468e6994dd7
md"""
Podemos accesar con el índice a esta tupla,
"""

# ╔═╡ 19f89894-779a-4498-bbc5-a586f26904fe
myfavoriteanimals[1]

# ╔═╡ b5b36260-b451-4537-801b-a0c90963fce0
md"""
Pero como las tuplas con inmutables, no la podemos modificar
"""

# ╔═╡ 41d46a14-3604-41d9-b8f7-3f4961dd66c8
mis_animales_favoritos[1] = "nutria"

# ╔═╡ 4ac9e892-3ac0-4f76-bc89-fee80ace061a
md"""
## Diccionarios

Si tenemos conjuntos de datos relacionados entre sí, podemos guardar los datos en un diccionario. Un buen ejemplo es una lista de contactos, donde asociamos nombres a números de teléfono.

Sintaxis:
```julia
Dict(llave1 => valor1, llave2 => valor2, ....)```
"""

# ╔═╡ 4ac82f64-1d92-4a63-ad3b-9d65df410245
miagenda = Dict("Jenny" => "867-5309", "Cazafantasmas" => "555-2368")

# ╔═╡ a4e6740c-947d-4413-9eac-d5252e6c23e3
md"""
En este ejemplo, cada nombre y número es un par de "llave" y "valor". Podemos tomar el númer de Jenny (un valor) usando la llave asociada.
"""

# ╔═╡ 0d2e6550-8ba0-4c97-a12a-911125e71237
miagenda["Jenny"]

# ╔═╡ 680c0f29-beac-49d6-9072-862c2c40c882
md"""
Podemos agregar otra entrada al diccionario de la manera siguiente
"""

# ╔═╡ 1bb288b2-522d-4846-879e-801ba9237967
miagenda["Kramer"] = "555-FILK"

# ╔═╡ 79ad3556-1c4b-40a1-bec9-26f407c39824
md"""
Veamos como se ve nuestro diccionario ahora...
"""

# ╔═╡ 4b30cdfc-8870-432e-8a00-2bcd9e6f719f
miagenda

# ╔═╡ 3276de5f-804a-46c3-b56c-dc59b4e0e5c6
md"""
Para borrar a Kramer de nuestro diccionario - y simultáneamente tomar su número - usamos  pop!
"""

# ╔═╡ f957b7d5-aaf4-456d-b820-66e1cc1b1fb2
pop!(miagenda, "Kramer")

# ╔═╡ c016c9d6-bea4-431d-b4fe-428798797e9b
miagenda

# ╔═╡ f391e572-57d7-4bbf-8cb4-8ffcf3477eb1
md"""
A diferencia de las tuplas y los arreglos, los diccionarios no están ordenados y no podemos accesarlos con un índice
"""

# ╔═╡ 6e2a49b3-076d-49cf-918a-1a3f47adc968
miagenda[1]

# ╔═╡ 4f28e674-6bf5-487e-927d-1323740e5744
md"""
En el ejemplo anterior, `julia` piensa que estás tratando de accesar a un valor asociado a la llave `1`.
"""

# ╔═╡ 285cf3f2-6d51-4a62-89c1-fa4c1fcdce6e
md"""
## Arreglos

A diferencia de las tuplas, los arreglos son mutables. A diferencia de los diccionarios, los arreglos contienen secuencias ordenadas de elementos.

Podemos crear un arreglo encapsulando esta secuencia de elementos con `[ ]`.

Sintaxis:

```julia
[item1, item2, ...]```

Por ejemplo, podemos usar un arreglo para recordar a mis amigos
"""

# ╔═╡ 2a19ba38-cbfa-4c3d-9976-4a1cd21f4313
myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

# ╔═╡ 51e4ab4e-75e7-43b0-ba5c-33742dcd013c
md"""
O guardar una secuencia de números
"""

# ╔═╡ b80b646a-71fd-4ffd-942a-7ef55f729402
fibonacci = [1, 1, 2, 3, 5, 8, 13]

# ╔═╡ 764cbdbc-44f5-4296-91ae-2e32ce3e98e0
mezcla = [1, 1, 2, 3, "Ted", "Robyn"]

# ╔═╡ a2366d3b-d945-422d-b0f0-d31d110899ef
md"""
Una vez que tenenmos un arreglo, podemos tomar los datos individualmente dentro del arreglo accesándolos por su índice. Por ejemplo, si queremos al tercer amigo en myfriends, escribimos
"""

# ╔═╡ 1b4b94c0-18e6-4d6a-bf5b-fb8a144e0a13
myfriends[3]

# ╔═╡ 82ee1763-e192-488d-a11a-5c2138e9ca54
md"""
Podemos usar el índice para mutar la entrada del arreglo
"""

# ╔═╡ 66c708d0-c7b7-47d1-8256-5c7e00cc4f51
myfriends[3] = "Baby Bop"

# ╔═╡ 8ec49414-da74-41c8-a2b1-55722966ce8f
md"""
También se puede editar con las fuciones de `push!` y `pop!`. `push!` agrega un elemento al final del arreglo y `pop!` quita al último elemento del arreglo.

Podemos agregar otro número a la secuencia fibonacci
"""

# ╔═╡ efd0045d-1f74-4ba0-988c-88897a8fea05
push!(fibonacci, 21)

# ╔═╡ 097606fc-7e9c-4c18-a5e6-027b352d0303
md"""
y quitarlo
"""

# ╔═╡ 78d6f8b1-5f07-45a4-ae60-658a3dd3fc02
pop!(fibonacci)

# ╔═╡ aaf34386-b771-4b0b-9d59-9ae24f4e8db8
fibonacci

# ╔═╡ 1f56308b-608a-4bac-b2ed-427d652d9f98
md"""
Hasta ahora hemos dado ejemplos de arreglos de escalars unidimensionales, pero los arreglos pueden tener un número arbitrario de dimensiones y también pueden guardar otros arreglos.

Por ejemplo, estos son arreglos de arreglos
"""

# ╔═╡ e47cb827-15e0-4ec0-b35f-5dc823844ac3
favoritos = [["koobideh", "chocolate", "eggs"],["penguins", "cats", "sugargliders"]]

# ╔═╡ a21f2ecb-7dbb-44df-8e81-55b7ed249161
numbers = [[1, 2, 3], [4, 5], [6, 7, 8, 9]]

# ╔═╡ eb301c3e-035f-42b7-aec5-66f53f085553
md"""
Abajo hay arreglos de 2D y 3D poblados con valores aleatorios
"""

# ╔═╡ 3a6a5316-9f10-40f1-9424-7d932c04be29
rand(4, 3)

# ╔═╡ 6fc41d26-ebc4-46ff-ac4e-075a903178e0
rand(4, 3, 2)

# ╔═╡ ea1953f7-9b6c-4698-8b50-e1ac43e04068
md"""
¡Cuidado copiando los arreglos!
"""

# ╔═╡ c6015f8e-87dd-408a-9491-8164f8ff9a40
fibonacci

# ╔═╡ eaeb74ea-6e64-4d3e-afbe-e9d5f518a618
somenumbers = fibonacci

# ╔═╡ 4cd1a54f-050c-4db6-b3de-957ee37f6c7f
somenumbers[1] = 404

# ╔═╡ 6f9e6b20-6dd7-47dc-b1c9-f599dcceaea0
fibonacci

# ╔═╡ fb09bb5b-9f80-47fe-9c75-2d9017a14271
md"""
Editar `somenumbers` causa que `fibonacci` se edite también!

En el ejemplo superior, en realidad no hicimos una copia de `fibonacci`. Sólo creamos una nueva manera de accesar las entradas del arreglo relacionado a `fibonacci`.

Si queremos hacer una copia de un arreglo amarrado a `fibonacci`, usamos la función de `copy`.
"""

# ╔═╡ a448ed00-e41b-4686-8c09-127ad496914b
# Primero restauramos a fibonnaci
fibonacci[1] = 1
fibonacci

# ╔═╡ 366a0922-c471-484d-9e3f-39b4117e705a
somemorenumbers = copy(fibonacci)

# ╔═╡ 57c6edf7-2d12-4d48-816c-3b32de841b3f
somemorenumbers[1] = 404

# ╔═╡ 636ee6ad-7010-4c50-9582-b8431dee69e0
fibonacci

# ╔═╡ 74e89a6d-f8bd-4744-aad3-7a5c109ae9f4
md"""
En el último ejemplo, no se editó a fibonacci. Entonces vemos que los arreglos amarrados a `somemorenumbers` y `fibonacci` son distintos.
"""

# ╔═╡ 9b65e734-61c3-4059-9539-79711d7cd48f
md"""
### Exercises

3.1 Crea un arreglo, `arreglo`, que sea un arreglo de 1D de 2-elementos de 1-elemento de 1D, cada uno guardando el número 0.
Accesa a `arreglo` para agregar un `1` a cada uno de los arreglos que contiene.
"""

# ╔═╡ cf0b40c8-cf5b-473d-a0c9-4248198d395e

# ╔═╡ d6cc5532-87c1-409a-9661-a4edc42a3177
md"""
3.2 Trata de agregar "Emergencia" a `miagenda` con el valor `911`. Trata de agregar `911` como un entero y no como cadena. ¿Porqué no funciona?
"""

# ╔═╡ 56fcc723-4baf-4225-88f0-de331af394e1

# ╔═╡ da4abb22-bd0a-4eaf-89df-d9c95f183f6a
md"""
3.3 Crea un nuevo diccionario que se llame `agenda_flexible` que tenga el número de Jenny guardado como cadena y el de los Cazafantasmas como entero.
"""

# ╔═╡ c2763aaa-e0af-4d36-9997-c62cb85c10d4

# ╔═╡ 66aa9ebe-09ea-4bab-8bed-90d79cdc2826
md"""
3.4 Agrega la llave de "Emergencia al valor (entero) `911` a `agenda_flexible`.
"""

# ╔═╡ 2fa4912f-569d-48ad-81ed-9609174b967b

# ╔═╡ 15636609-63bd-4ff1-88b3-497e54deb1c7
md"""
3.5 ¿Porqué podemos agregar un entero como valor a `agenda_flexible` pero no a  `miagenda`? ¿Cómo pudimos haber inicializado `miagenda` para que aceptara enteros como valores?
"""

# ╔═╡ b42e131e-6d50-40c0-9a22-762612efef0a

# ╔═╡ Cell order:
# ╟─f906f80a-9eec-40cb-80ba-bdb06a798b1b
# ╟─67852a29-1c6c-4ab6-be1e-35a16be10a40
# ╠═5f86014d-54f2-497e-b6d3-529ff3e6f43e
# ╟─b21a8b34-0b8d-4c06-b38c-f468e6994dd7
# ╠═19f89894-779a-4498-bbc5-a586f26904fe
# ╟─b5b36260-b451-4537-801b-a0c90963fce0
# ╠═41d46a14-3604-41d9-b8f7-3f4961dd66c8
# ╟─4ac9e892-3ac0-4f76-bc89-fee80ace061a
# ╠═4ac82f64-1d92-4a63-ad3b-9d65df410245
# ╟─a4e6740c-947d-4413-9eac-d5252e6c23e3
# ╠═0d2e6550-8ba0-4c97-a12a-911125e71237
# ╟─680c0f29-beac-49d6-9072-862c2c40c882
# ╠═1bb288b2-522d-4846-879e-801ba9237967
# ╟─79ad3556-1c4b-40a1-bec9-26f407c39824
# ╠═4b30cdfc-8870-432e-8a00-2bcd9e6f719f
# ╟─3276de5f-804a-46c3-b56c-dc59b4e0e5c6
# ╠═f957b7d5-aaf4-456d-b820-66e1cc1b1fb2
# ╠═c016c9d6-bea4-431d-b4fe-428798797e9b
# ╟─f391e572-57d7-4bbf-8cb4-8ffcf3477eb1
# ╠═6e2a49b3-076d-49cf-918a-1a3f47adc968
# ╟─4f28e674-6bf5-487e-927d-1323740e5744
# ╟─285cf3f2-6d51-4a62-89c1-fa4c1fcdce6e
# ╠═2a19ba38-cbfa-4c3d-9976-4a1cd21f4313
# ╟─51e4ab4e-75e7-43b0-ba5c-33742dcd013c
# ╠═b80b646a-71fd-4ffd-942a-7ef55f729402
# ╠═764cbdbc-44f5-4296-91ae-2e32ce3e98e0
# ╟─a2366d3b-d945-422d-b0f0-d31d110899ef
# ╠═1b4b94c0-18e6-4d6a-bf5b-fb8a144e0a13
# ╟─82ee1763-e192-488d-a11a-5c2138e9ca54
# ╠═66c708d0-c7b7-47d1-8256-5c7e00cc4f51
# ╟─8ec49414-da74-41c8-a2b1-55722966ce8f
# ╠═efd0045d-1f74-4ba0-988c-88897a8fea05
# ╟─097606fc-7e9c-4c18-a5e6-027b352d0303
# ╠═78d6f8b1-5f07-45a4-ae60-658a3dd3fc02
# ╠═aaf34386-b771-4b0b-9d59-9ae24f4e8db8
# ╟─1f56308b-608a-4bac-b2ed-427d652d9f98
# ╠═e47cb827-15e0-4ec0-b35f-5dc823844ac3
# ╠═a21f2ecb-7dbb-44df-8e81-55b7ed249161
# ╟─eb301c3e-035f-42b7-aec5-66f53f085553
# ╠═3a6a5316-9f10-40f1-9424-7d932c04be29
# ╠═6fc41d26-ebc4-46ff-ac4e-075a903178e0
# ╟─ea1953f7-9b6c-4698-8b50-e1ac43e04068
# ╠═c6015f8e-87dd-408a-9491-8164f8ff9a40
# ╠═eaeb74ea-6e64-4d3e-afbe-e9d5f518a618
# ╠═4cd1a54f-050c-4db6-b3de-957ee37f6c7f
# ╠═6f9e6b20-6dd7-47dc-b1c9-f599dcceaea0
# ╟─fb09bb5b-9f80-47fe-9c75-2d9017a14271
# ╠═a448ed00-e41b-4686-8c09-127ad496914b
# ╠═366a0922-c471-484d-9e3f-39b4117e705a
# ╠═57c6edf7-2d12-4d48-816c-3b32de841b3f
# ╠═636ee6ad-7010-4c50-9582-b8431dee69e0
# ╟─74e89a6d-f8bd-4744-aad3-7a5c109ae9f4
# ╟─9b65e734-61c3-4059-9539-79711d7cd48f
# ╠═cf0b40c8-cf5b-473d-a0c9-4248198d395e
# ╟─d6cc5532-87c1-409a-9661-a4edc42a3177
# ╠═56fcc723-4baf-4225-88f0-de331af394e1
# ╟─da4abb22-bd0a-4eaf-89df-d9c95f183f6a
# ╠═c2763aaa-e0af-4d36-9997-c62cb85c10d4
# ╟─66aa9ebe-09ea-4bab-8bed-90d79cdc2826
# ╠═2fa4912f-569d-48ad-81ed-9609174b967b
# ╟─15636609-63bd-4ff1-88b3-497e54deb1c7
# ╠═b42e131e-6d50-40c0-9a22-762612efef0a

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
using StochasticDiffEq, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: oval2ModelExample
Random.seed!(200)
prob = oval2ModelExample(largeFluctuations=true,useBigs=false)
quick_prob = deepcopy(prob)
quick_prob.tspan = (0.0,1.0)

using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 0.0001
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100
@benchmark begin
  Random.seed!(100)
  sol = solve(quick_prob,SRIW1(),dt=(1/2)^(18),progress_steps=Int(1e5),
        adaptivealg=:RSwM3,progress=false,qmax=4,save_everystep=false,
        saveat=0.1,abstol=1e-5,reltol=1e-3)
end

using ProfileView
Random.seed!(100)
@profile sol = solve(quick_prob,SRIW1(),dt=(1/2)^(18),progress_steps=Int(1e5),
      adaptivealg=:RSwM3,progress=false,qmax=4,save_everystep=false,
      saveat=0.1,abstol=1e-5,reltol=1e-3)

ProfileView.view()

Random.seed!(100)
@time sol = solve(quick_prob,SRIW1(),dt=(1/2)^(18),progress_steps=Int(1e5),
      progress=false,qmax=4,save_everystep=false,
      saveat=0.1,abstol=1e-5,reltol=1e-3)

println(sol.u[end])

#[0.011503,0.939809,0.00312214,0.00155873,0.0172325,0.0577331,0.237757,0.00134921,0.000238022,4.19916e-5,7.40824e-6,1.307e-6,0.0621091,1.24463,0.0483949,199.901,137.457,0.0177237,0.132583]

Random.seed!(100)
@time integrator = init(quick_prob,SRIW1(),dt=(1/2)^(18),progress_steps=Int(1e5),
      adaptivealg=:RSwM3,progress=false,qmax=4,save_everystep=false,
      saveat=0.1,abstol=1e-5,reltol=1e-3)

@time StochasticDiffEq.solve!(integrator)

@code_warntype StochasticDiffEq.generate_tildes(integrator,1.0,1.0,1.0)
@code_warntype StochasticDiffEq.perform_step!(integrator)
@code_llvm DiffEqBase.isinplace(integrator.noise)
function test_pop(integrator)
  b2 = integrator.opts.beta2
end
@code_warntype test_pop(integrator)


text:
# ============================================================================
# Example 1 from *Reachability analysis of linear systems with uncertain
# parameters and inputs*. TODO: Agregar referencia completa.
# ============================================================================

# initial set
X0 = BallInf([1.0, 1.0], 0.1)

# linear ODE: x' = Ax
A = IntervalMatrix([-1.0 ± 0.05 -4.0 ± 0.05;
                    4.0 ± 0.05 -1.0 ± 0.05])

# IVP(LCS(A), X0) TODO: remove
interval2D_linear = @ivp x' = Ax, x(0) ∈ X0

# affine ODE: x' = Ax + Bu
B = IntervalMatrix(hcat([1.0 ± 0.0; 1.0 ± 0.0])) # why hcat?
U = Interval(-0.05, 0.05)

# IVP(CLCCS(A, B, nothing, U), X0) TODO: remove
interval2D_affine = @ivp x' = Ax + Bu, x(0) ∈ X0, u ∈ U