text:
mutable struct DemoTrainingDQN
    # network weights
    w::Vector{Array{Float32,N} where N}

    # tree params
    max_stepc::Int64
    epsilon::Float64
    MIN_EPSILON::Float64
    fac::Float64
    gamma::Float64

    # training params
    LR::Float64
    batchsize::Int64
    train_epochs::Int64
    max_states::Int64
    train_when_min::Int64

    # for validation
    val_safety::Float64

    # experience buffer
    sensors::Vector{Array{UInt8,2}}
    qs::Vector{Vector{Float32}}
    priorities::Vector{Float64}

    # for stats
    rewards::Vector{Float64}
    ranks::Vector{Int64}
    stepscum::Vector{Int64}
    episode_mean_speed::Vector{Float64}
    episode_std_speed::Vector{Float64}
    val_mean_speed::Vector{Float64}
    val_std_speed::Vector{Float64}
    val_steps::Vector{Float64}
    val_rewards::Vector{Float64}
    val_stepscum::Vector{Int64}
    avg_window::Int64
    _stepscum::Int64
    iter_of_training_start::Int64
    epsilon0::Float64
    function DemoTrainingDQN()
        new(
            Vector{Array{Float32,N} where N}(),
            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            Vector{Array{UInt8,2}}(), Vector{Vector{Float32}}(),
            Vector{Float64}(), Vector{Float64}(),
            Vector{Int64}(), Vector{Int64}(),
            Vector{Float64}(), Vector{Float64}(),
            Vector{Float64}(), Vector{Float64}(),
            Vector{Float64}(), Vector{Float64}(),
            Vector{Int64}(),
            0, 0, 0, 0,
        )
    end
end

function InitializeDTDQN()
    dtdqn = DemoTrainingDQN()
    dtdqn.w = initialize_weights(actionc_steer*actionc_accel)
    dtdqn.max_stepc = 5500
    dtdqn.epsilon = 0.5
    dtdqn.MIN_EPSILON = 0.01
    dtdqn.fac = 0.995
    dtdqn.gamma = 0.98
    dtdqn.LR = 0.01
    dtdqn.batchsize = 64
    dtdqn.train_epochs = 20
    dtdqn.max_states = 80000
    dtdqn.train_when_min = 50000
    dtdqn.val_safety = 1.0
    dtdqn.avg_window = 50
    dtdqn._stepscum = 0
    dtdqn.iter_of_training_start = -1
    dtdqn.epsilon0 = dtdqn.epsilon
    dtdqn
end

function traindtdqn(dtdqn::DemoTrainingDQN, iters::Integer, fname)

    w = map(KnetArray, dtdqn.w)
    for i in 1:iters
        #### Accumulating episode ###################################################

        println("--------------------------------------------------------------- i="*string(i))
        flush(stdout)

        added_states = 0
        episodsc = 0
        tot_episode_len = 0
        tot_passed = 0
        while true
            # accumulate episodes
            _speeds = Float64[]
            _sensors = Array{UInt8,2}[]
            _qs = Vector{Float32}[]
            _rewards = Float64[]
            rank = 0

            episodsc += 1
            state, env = initialize_simple_road();
            sensor = get_sensors(env, state)
            reward = 0.0f0

            while true
                # single episode loop
                dtdqn._stepscum += 1
                rank += 1

                sensor = get_sensors(env, state)
                q = network_predict(env, w, sensor)

                push!(_speeds, state.speed)
                push!(_sensors, sensor)
                push!(_qs, q)
                push!(_rewards, reward)

                action_ix = argmax(q)
                action = action_ix_to_action(env, Int32(action_ix))
                state, reward, done = sim!(env, state, action, 5, dtdqn.val_safety)
                if done
                    sensor = get_sensors(env, state)
                    q = network_predict(env, w, sensor)
                    q .= 0.0f0

                    push!(_speeds, state.speed)
                    push!(_sensors, sensor)
                    push!(_qs, q)
                    push!(_rewards, reward)
                end
                done && break
                rank > dtdqn.max_stepc && break
            end
            # n-step backprop
            for i in length(_qs):-1:2
                _qs[i-1][argmax(_qs[i-1])] = _rewards[i] + dtdqn.gamma*maximum(_qs[i])
            end
            append!(dtdqn.qs, _qs)
            append!(dtdqn.sensors, _sensors)
            f = if length(dtdqn.priorities)==0
                1.0
            else
                median(dtdqn.priorities)
            end
            append!(dtdqn.priorities, f.*ones(Float32, length(_qs)))

            push!(dtdqn.stepscum, dtdqn._stepscum)
            push!(dtdqn.ranks, rank)
            push!(dtdqn.episode_mean_speed, mean(_speeds))
            push!(dtdqn.episode_std_speed, std(_speeds))
            push!(dtdqn.rewards, sum(_rewards))
            if length(dtdqn.qs) > dtdqn.max_states
                dtdqn.qs = dtdqn.qs[(end-dtdqn.max_states+1):end]
                dtdqn.sensors = dtdqn.sensors[(end-dtdqn.max_states+1):end]
                dtdqn.priorities = dtdqn.priorities[(end-dtdqn.max_states+1):end]
            end
            added_states += length(_qs)
            added_states >= dtdqn.max_stepc-1 && break
            tot_episode_len += rank
            tot_passed += state.times_passed_yolocar
            print("."); flush(stdout)
        end
        mean_episode_len = tot_episode_len / episodsc
        mean_passed = tot_passed / episodsc

        #### Training #######################################################

        if length(dtdqn.qs) >= dtdqn.train_when_min
            if dtdqn.iter_of_training_start<0
                dtdqn.iter_of_training_start = div(dtdqn._stepscum, dtdqn.max_stepc)
            end
            dtdqn.epsilon = max(dtdqn.MIN_EPSILON, dtdqn.epsilon0 * dtdqn.fac^(div(dtdqn._stepscum, dtdqn.max_stepc)-dtdqn.iter_of_training_start))
            N = round(Int64, added_states * dtdqn.train_epochs / dtdqn.batchsize)
            println("\n---- training for "*string(N)*" iterations"); flush(stdout)
            w,_ = trainit(N, w, dtdqn.sensors,dtdqn.qs,dtdqn.LR,dtdqn.batchsize, dtdqn.priorities)
            dtdqn.w = map(Array, w)
        end

        #### Validating #####################################################

        valv = Float64[]
        valr = 0.0
        spd = mean(dtdqn.episode_mean_speed[max(1,end-dtdqn.avg_window):end])
        if isnan(spd)
            spd = 0.8
        end
        state, val_env = initialize_simple_road(spd);
        for _ in 1:dtdqn.max_stepc
            push!(valv, state.speed)
            _sensors = get_sensors(val_env, state)
            _vqs = network_predict(val_env, w, _sensors)
            action = action_ix_to_action(val_env,Int32(argmax(_vqs)))
            state, reward, done = sim!(val_env, state, action, 5, dtdqn.val_safety)
            valr += reward
            done && break
        end
        vpassed = state.times_passed_yolocar
        push!(dtdqn.val_rewards, valr)
        push!(dtdqn.val_mean_speed, mean(valv))
        push!(dtdqn.val_std_speed, std(valv))
        push!(dtdqn.val_steps, length(valv))
        push!(dtdqn.val_stepscum, dtdqn._stepscum)

        #### Reporting #####################################################

        println()
        @show length(dtdqn.qs), dtdqn.rewards[end], div(dtdqn._stepscum, dtdqn.max_stepc)
        println("         eps = "*string(round(dtdqn.epsilon, digits=4))*"    gamma ="*string(round(dtdqn.gamma, digits=4)))
        println("         rnk = "*string(mean_episode_len)*             "      mvel = "*string(round(dtdqn.episode_mean_speed[end], digits=4))*"    svel = "*string(round(dtdqn.episode_std_speed[end], digits=4)))
        println("        vrnk = "*string(dtdqn.val_steps[end])*         "     vmvel = "*string(round(dtdqn.val_mean_speed[end], digits=4))*"   vsvel = "*string(round(dtdqn.val_std_speed[end], digits=4)))
        println(" mean_passed = "*string(floor(Int64, mean_passed))*    "   vpassed = "*string(vpassed))
        flush(stdout)

        if i%10 == 0
            IJulia.clear_output();
        end

        if i%100 == 0
            @save fname dtdqn
        end

    end
end

function plotdtdqn(dtdqn)
    p1 = plot(dtdqn.stepscum/5500, moving_avg(dtdqn.rewards, dtdqn.avg_window), label="train reward")
    plot!(p1, dtdqn.val_stepscum/5500, moving_avg(dtdqn.val_rewards, dtdqn.avg_window), lw=3, label="val reward")

    p2 = plot(dtdqn.stepscum/5500, moving_avg(dtdqn.ranks/dtdqn.max_stepc, dtdqn.avg_window), label="train rank (perc.)")
    plot!(p2, dtdqn.val_stepscum/5500, moving_avg(dtdqn.val_steps/dtdqn.max_stepc, dtdqn.avg_window), lw=3, label="max tree rank (perc.)")

    p3 = plot(dtdqn.stepscum/5500, moving_avg(dtdqn.episode_mean_speed, dtdqn.avg_window), label="train mean speed")
    plot!(p3, dtdqn.val_stepscum/5500, moving_avg(dtdqn.val_mean_speed, dtdqn.avg_window), lw=3, label="val mean speed")

    plot(p1,p2,p3, layout=(3,1), size=(900,1100))
end


text:
module Retriever

using PyCall

export check_for_updates, dataset_names, download
export install_csv, install_mysql, install_postgres
export install_sqlite, install_msaccess, install_json
export install_xml, reset_retriever

rt = pyimport("retriever")

"""
```julia
    check_for_updates()
```

Check Retriever scripts for updates.
"""
function check_for_updates()
    rt.check_for_updates()
end

"""
```julia
    reload_scripts()
```

Retriever reload scripts.
"""
function reload_scripts()
    rt.reload_scripts()
end

"""
```julia
    reset_retriever(scope::String="all", ask_permission::Bool=true)
```

Retriever reset scripts or data or all.
"""
function reset_retriever(scope::String="all", ask_permission::Bool=true)
    rt.reset_retriever(scope, ask_permission)
end

"""
```julia
    dataset_names()
```

Return list of all available dataset names.
"""
function dataset_names()
    rt.dataset_names()
end

"""
```julia
    download(dataset; path::String="./", quite::Bool=false,
                subdir::String="", debug::Bool=false, use_cache::Bool=false)
```

Download scripts for retriever.
"""
function download(dataset; path::String="./", quite::Bool=false,
                subdir::String="", debug::Bool=false, use_cache::Bool=false)
    rt.download(dataset, path, quite, subdir, debug, use_cache)
end

"""
```julia

    install_csv(dataset; table_name::String="", compile::Bool=false,
        data_dir=pwd(), debug::Bool=false, use_cache::Bool=true)
```

Install Retriever scripts in csv format.
"""
function install_csv(dataset; table_name::String="{db}_{table}",
        data_dir=pwd(), debug::Bool=false, use_cache::Bool=true)
    rt.install_csv(dataset, table_name, data_dir, debug, use_cache)
end

"""
```julia
    install_mysql(dataset; user::String="root",
                password::String="", host::String="localhost",
                port::Int=3306, database_name::String="", table_name::String="",
                compile::Bool=false, debug::Bool=false, quite::Bool=false,
                use_cache::Bool=true)
```

Install Retriever scripts in mysql database.
"""
function install_mysql(dataset; user::String="root",
                password::String="", host::String="localhost",
                port::Int=3306, database_name::String="{db}", table_name::String="{db}.{table}",
                debug::Bool=false, use_cache::Bool=true)
    rt.install_mysql(dataset, user, password, host, port, database_name,
                     table_name, debug, use_cache)
end

"""
```julia
    install_postgres(dataset; user::String="postgres",
                password::String="", host::String="localhost",
                port::Int=5432, database::String="postgres",
                database_name::String="", table_name::String="", bbox::Any=[],
                debug::Bool=false, use_cache::Bool=true)
```

Install Retriever scripts in database.
"""
function install_postgres(dataset; user::String="postgres",
                password::String="", host::String="localhost",
                port::Int=5432, database::String="postgres",
                database_name::String="{db}", table_name::String="{db}.{table}", bbox::Any=[],
                debug::Bool=false, use_cache::Bool=true)

    rt.install_postgres(dataset, user, password, host, port, database,
                        database_name, table_name, bbox, debug, use_cache)
end

"""
```julia
    install_sqlite(dataset; file::String="", table_name::String="",
                data_dir=pwd(), debug::Bool=false,
                use_cache::Bool=true)
```

Install Retriever scripts in database.
"""
function install_sqlite(dataset; file::String="", table_name::String="{db}.{table}",
                data_dir=pwd(), debug::Bool=false, use_cache::Bool=true)
    return rt.install_sqlite(dataset, file, table_name, data_dir, debug, use_cache)
end

"""
```julia
    install_msaccess(dataset; file::String="", table_name::String="",
                data_dir=pwd(), debug::Bool=false, use_cache::Bool=true)
```

Install Retriever scripts in msacces.
"""
function install_msaccess(dataset; file::String="", table_name::String="",
                data_dir=pwd(), debug::Bool=false, use_cache::Bool=true)
    rt.install_msaccess(dataset, file, table_name, data_dir, debug, quite,
                        use_cache)
end

"""
```julia
    install_json(dataset; table_name::String="", data_dir=pwd(),
    debug::Bool=false, use_cache::Bool=true)
```

Install Retriever scripts in json format.
"""
function install_json(dataset; table_name::String="{db}_{table}.json",
                data_dir=pwd(), debug::Bool=false, use_cache::Bool=true)
    rt.install_json(dataset, table_name, data_dir, debug, use_cache)

end

"""
```julia
    install_xml(dataset; table_name::String="", compile::Bool=false, data_dir=pwd(),
            debug::Bool=false, use_cache::Bool=true)
```

Install Retriever scripts in xml format.
"""
function install_xml(dataset; table_name::String="{db}_{table}.xml",
                data_dir=pwd(), debug::Bool=false, use_cache::Bool=true)
    rt.install_xml(dataset, table_name, data_dir, debug, use_cache)

end

"""
```julia
    reset_retriever(; scope::String="all")
```

Remove stored information on scripts, data, and connection
"""
function reset_retriever(; scope::String="all")
    rt.reset_retriever(scope)
end

end # module

ext: jl


text:
#TODO: do binary...although this is exactly what BigInt does

function serialize(s::AbstractSerializer, t::fmpz)
  Serialization.serialize_type(s, fmpz)
  return serialize(s, base(t, 62))
end

function deserialize(s::AbstractSerializer, ::Type{fmpz})
  return parse(fmpz, deserialize(s), 62)
end

function serialize(s::AbstractSerializer, t::fmpq)
  Serialization.serialize_type(s, fmpq)
  serialize(s, base(numerator(t), 62))
  return serialize(s, base(denominator(t), 62))
end

function deserialize(s::AbstractSerializer, ::Type{fmpq})
  n = parse(fmpz, deserialize(s), 62)
  d = parse(fmpz, deserialize(s), 62)
  return fmpq(n, d)
end

function serialize(s::AbstractSerializer, t::PolyElem{T}) where T
  Serialization.serialize_type(s, PolyElem{T})
  serialize(s, length(t))
  for i=0:length(t)
    serialize(s, coeff(t, i))
  end
end

function deserialize(s::AbstractSerializer, ::Type{PolyElem{T}}) where T
  L = T[]
  l = deserialize(s)
  for i=0:l
    push!(L, deserialize(s))
  end
  R = parent(L[1])
  Rx, x = PolynomialRing(R, :S, cached=false)
  return Rx(L)
end

function serialize(s::AbstractSerializer, t::AnticNumberField)
  Serialization.serialize_type(s, AnticNumberField)
  return serialize(s, t.pol)
end

function deserialize(s::AbstractSerializer, ::Type{AnticNumberField})
  return number_field(deserialize(s), cached=false)[1]
end

function serialize(s::AbstractSerializer, t::NfToNfMor)
  Serialization.serialize_type(s, NfToNfMor)
  return serialize(s, (domain(t), codomain(t), t.prim_img))
end

function deserialize(s::AbstractSerializer, ::Type{NfToNfMor})
  K, L, a = deserialize(s)
  return NfToNfMor(K, L, a)
end

add_verbose_scope(:Par)

function _bizarre(a::Int, b::Int)
  return length(Hecke.class_group(Hecke.wildanger_field(a, b)[1])[1])
end

function parallel_work(f::Function, r::UnitRange)
  l = length(workers())
  @vprint :Par 1 "using $l workers\n"
  cache = []
  st = start(r)
  res = []
  c = Channel(2*l+1)
  no_job = 0
  tw = 0.0
  job_stat = Dict(w => 0 for w = workers())
  jobs = Dict{Int, Any}()
  while !done(r, st)
    while no_job < l
      if done(r, st)
        break
      end
      v, st = next(r, st)
      no_job += 1
      d = Distributed.nextproc()
      job_stat[d] += 1
      jobs[d] = time()
      @vtime :Par 2 @async put!(c, remotecall_fetch(f, d, v, d))
    end
    tw -= time()
    @vtime :Par 2 wait(c)
    tw += time()
    while isready(c)
      @vtime :Par 2 t = take!(c)
      push!(res, (t, time() - jobs[t[4]]))
      no_job -= 1
    end
  end
  @vprint :Par 2 "Finishing off..\n"
  push!(res, ((1,2,3.0,4), 5.0))
  while no_job > 0
    tw -= time()
    @vtime :Par 2 wait(c)
    tw += time()
    t = take!(c)
    push!(res, (t, time() - jobs[t[4]]))
    no_job -= 1
  end
  println("total waiting time was $tw")
  for (w,v) = job_stat
    println("$w -> $v")
  end
  return res
end

function addprocs(machines::AbstractVector; tunnel=false, sshflags=``, max_parallel=10, kwargs...)
  new = Base.Distributed.addprocs(Base.Distributed.SSHManager(machines); tunnel=tunnel, sshflags=sshflags, max_parallel=max_parallel, kwargs...)
  for d = new
    eval(Main, :(@spawnat $d eval(Main, :(include("/home/fieker/.julia/v0.6/Nemo/src/Nemo.jl");include("/home/fieker/.julia/v0.6/Hecke/src/Hecke.jl")))))
  end
  return new
end


text:
srand(0) #to make the tests repeatable

unittests = [
    "policies",
    "parameters",
    "data",
    "samples",
    "distributions",
    "proposals",
    "samplers",
    "indicators",
    "models",
    "testutil",
    "tuners",
    "chains",
    "jobsegments",
    "runners",
    "mhrunner",
    "smhrunner",
    "gmhrunner"]

println("===================")
println("Running unit tests:")
println("===================")

for t in unittests
    tfile = t*".jl"
    println("  * $(tfile) *")
    include(string("unit/",tfile))
    println()
    println()
end


text:
@testset "chainrules" begin
    Random.seed!(1)

    @testset "general: single input" begin
        for x in (1.0, -1.0, 0.0, 0.5, 10.0, -17.1, 1.5 + 0.7im)
            test_scalar(erf, x)
            test_scalar(erfc, x)
            test_scalar(erfcx, x)
            test_scalar(erfi, x)

            test_scalar(airyai, x)
            test_scalar(airyaiprime, x)
            test_scalar(airybi, x)
            test_scalar(airybiprime, x)

            test_scalar(dawson, x)

            if x isa Real
                test_scalar(logerfc, x)
                test_scalar(logerfcx, x)

                test_scalar(invdigamma, x)
            end

            if x isa Real && 0 < x < 1
                test_scalar(erfinv, x)
                test_scalar(erfcinv, x)
            end

            if x isa Real && x > 0 || x isa Complex
                test_scalar(gamma, x)
                test_scalar(digamma, x)
                test_scalar(trigamma, x)

                if x isa Real
                    test_scalar(airyaix, x)
                    test_scalar(airyaiprimex, x)
                end
            end
        end
    end

    @testset "Bessel functions" begin
        for x in (1.5, 2.5, 10.5, -0.6, -2.6, -3.3, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
            test_scalar(besselj0, x)
            test_scalar(besselj1, x)

            isreal(x) && x < 0 && continue

            test_scalar(bessely0, x)
            test_scalar(bessely1, x)

            for nu in (-1.5, 2.2, 4.0)
                test_frule(besseli, nu, x)
                test_rrule(besseli, nu, x)
                test_frule(besselix, nu, x) # derivative is `NotImplemented`
                test_frule(besselix, nu ⊢ NoTangent(), x) # derivative is a number
                test_rrule(besselix, nu, x)

                test_frule(besselj, nu, x)
                test_rrule(besselj, nu, x)
                test_frule(besseljx, nu, x) # derivative is `NotImplemented`
                test_frule(besseljx, nu ⊢ NoTangent(), x) # derivative is a number
                test_rrule(besseljx, nu, x)

                test_frule(besselk, nu, x)
                test_rrule(besselk, nu, x)
                test_frule(besselkx, nu, x)
                test_rrule(besselkx, nu, x)

                test_frule(bessely, nu, x)
                test_rrule(bessely, nu, x)
                test_frule(besselyx, nu, x) # derivative is `NotImplemented`
                test_frule(besselyx, nu ⊢ NoTangent(), x) # derivative is a number
                test_rrule(besselyx, nu, x)

                test_frule(hankelh1, nu, x)
                test_rrule(hankelh1, nu, x)
                test_frule(hankelh1x, nu, x)
                test_rrule(hankelh1x, nu, x)

                test_frule(hankelh2, nu, x)
                test_rrule(hankelh2, nu, x)
                test_frule(hankelh2x, nu, x)
                test_rrule(hankelh2x, nu, x)
            end
        end
    end

    @testset "erf, beta, and logbeta" begin
        test_points = (1.5, 2.5, 10.5, 1.6 + 1.6im, 1.6 - 1.6im, 4.6 + 1.6im)
        for x in test_points, y in test_points
            test_frule(beta, x, y)
            test_rrule(beta, x, y)

            test_frule(logbeta, x, y)
            test_rrule(logbeta, x, y)

            if x isa Real && y isa Real
                test_frule(erf, x, y)
                test_rrule(erf, x, y)
            end
        end
    end

    @testset "log gamma and co" begin
        # It is important that we have negative numbers with both odd and even integer parts
        test_points = (1.5, 2.5, 10.5, -0.6, -2.6, -3.3, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
        for x in test_points
            for m in (0, 1, 2, 3)
                test_frule(polygamma, m, x; rtol=1e-8)
                test_rrule(polygamma, m, x; rtol=1e-8)
            end

            isreal(x) && x < 0 && continue
            test_scalar(loggamma, x)
            for a in test_points
                test_frule(gamma, a, x; rtol=1e-8)
                test_rrule(gamma, a, x; rtol=1e-8)

                test_frule(loggamma, a, x)
                test_rrule(loggamma, a, x)
            end

            isreal(x) || continue
            test_frule(logabsgamma, x)
            test_rrule(logabsgamma, x; output_tangent=(randn(), randn()))
            for a in test_points
                isreal(a) && a > 0 || continue
                test_frule(gamma_inc, a, x, 0)
                test_rrule(gamma_inc, a, x, 0; output_tangent=(randn(), randn()))
            end
        end
    end

    @testset "exponential integrals" begin
        for x in (1.5, 2.5, 10.5, 1.6 + 1.6im, 1.6 - 1.6im, -4.6 + 1.6im)
            test_scalar(expint, x)
            test_scalar(expintx, x)

            for nu in (-1.5, 2.2, 4.0)
                test_frule(expint, nu, x)
                test_rrule(expint, nu, x)

                test_frule(expintx, nu, x)
                test_rrule(expintx, nu, x)
            end

            isreal(x) || continue
            test_scalar(expinti, x)
            test_scalar(sinint, x)
            test_scalar(cosint, x)
        end
    end

    # https://github.com/JuliaMath/SpecialFunctions.jl/issues/307
    @testset "promotions" begin
        # one argument
        for f in (erf, erfc, logerfc, erfcinv, erfcx, logerfcx, erfi, erfinv, sinint)
            _, ẏ = frule((NoTangent(), 1f0), f, 1f0)
            @test ẏ isa Float32
            _, back = rrule(f, 1f0)
            _, x̄ = back(1f0)
            @test x̄ isa Float32
        end

        # two arguments
        _, ẏ = frule((NoTangent(), 1f0, 1f0), erf, 1f0, 1f0)
        @test ẏ isa Float32
        _, back = rrule(erf, 1f0, 1f0)
        _, x̄ = back(1f0)
        @test x̄ isa Float32
    end
end


text:
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIF = MOI.FileFormats

# Use this as cache
MOIU.@model(ClassicConeOptimizer,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
    (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction,),
    true
)

cache = ClassicConeOptimizer{Float64}()
bridge = MOI.Bridges.full_bridge_optimizer(cache, Float64)

# Read initial instance
finst = ARGS[1]
cbf = MOIF.CBF.Model()
MOI.read_from_file(cbf, finst)
MOI.copy_to(bridge, cbf)

#
cbf2 = MOIF.CBF.Model()
MOI.copy_to(cbf2, cache)
MOI.write_to_file(cbf2, ARGS[2])

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
# Implementation of the streaming LogSumExp trick
# http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html

mutable struct LogSumExpAcc{T<:Number}
    m::T # Current maximum value
    s::T # Sum
end

LogSumExpAcc() = LogSumExpAcc(-Inf, 0.0)

function add!(acc::LogSumExpAcc{T}, val::T) where T
    if val === -Inf
        return
    elseif val <= acc.m
        acc.s += exp(val - acc.m)
    else
        acc.s *= exp(acc.m - val)
        acc.s += 1.0
        acc.m = val
    end
end

function add!(acc::LogSumExpAcc{T}, vals::Vector{T}) where T
    for val in vals
        add!(acc, val)
    end
end

sum(acc::LogSumExpAcc{T}) where T = log(acc.s) + acc.m


text:
"""
    TypeQ()

Lorem ipsum dolor sit amet.
"""
struct TypeQ

end


text:
module BenchmarkAStarSearchCollatz

using BenchmarkTools
using AStarSearch

neighbours(n) = n % 2 == 0 ? [n ÷ 2, 3n + 1] : [3n + 1]

# The best case is if the number is a power of 2, then the moves are just the log2, otherwise this is a very optimistic heuristic (it doesn't overestimate)
heuristic(n, goal) = floor(Int64, log2(n))

astar_collatz(x; kwargs...) = astar(neighbours, x, 1; heuristic, kwargs...)

suite = BenchmarkGroup()

numbers = [1, 34, 523, 1345, 8765]
for number in numbers
  suite[[string(number)]] = @benchmarkable astar_collatz($number)
end

end

BenchmarkAStarSearchCollatz.suite


text:
@testset "16.3sum-closest.jl" begin
    @test three_sum_closest([-1, 2, 1, -4], 1) == 2
    @test three_sum_closest([0, 0, 0], 1) == 0
end


text:
module CthulhuVision

include("random.jl")
include("math.jl")
include("light.jl")
include("image.jl")
include("materials.jl")
include("camera.jl")
include("triangle.jl")
include("world.jl")
include("bvh.jl")
include("scene.jl")
include("rendering.jl")

end


text:
### Generalized Metropolis-Hastings Monte Carlo runner
immutable GMHRunner <: AbstractMHRunner
    numburnin::Int
    numiterations::Int
    numproposals::Int
    numindicatorsamples::Int
    policy::MHRuntimePolicy #policies determining runtime behaviour of the MCMC

    function GMHRunner(nburnin::Int,niterations::Int,nproposals::Int,nindsamples::Int,policy::MHRuntimePolicy)
        @assert 0 <= nburnin && 0 < niterations && 0 < nproposals && 0 < nindsamples
        new(nburnin,niterations,nproposals,nindsamples,policy)
    end
end

###Factory functions
_runner(::Type{Val{:generalized}},policy::MHRuntimePolicy,niterations::Int,nproposals::Int;numburnin =0,numindicatorsamples =nproposals) = GMHRunner(numburnin,niterations,nproposals,numindicatorsamples,policy)

function run!(runner_::GMHRunner,model_::AbstractModel,sampler_::AbstractSampler,tuner_::AbstractTuner)
    indicatorstate_ = samplerstate(sampler_,1,runner_.policy.sampletype,runner_.policy.calculationtype)
    segments_ = remotesegments(runner_.policy,model_,sampler_,runner_.numproposals)
    show(segments_)
    indicator_,tunerstate_,chain_ = createcommon(runner_,tuner_,numparas(model_),numtotalproposals(segments_),runner_.numindicatorsamples)
    tic()
    indicatorstate_ = burnin!(runner_,model_,sampler_,indicatorstate_,segments_,tuner_,tunerstate_,indicator_,chain_)
    indicatorstate_ = main!(runner_,model_,sampler_,indicatorstate_,segments_,tuner_,tunerstate_,indicator_,chain_)
    chain_.runtime = toq()
    chain_
end

function burnin!(runner_::GMHRunner,model_::AbstractModel,sampler_::AbstractSampler,indicatorstate_::AbstractSamplerState,segments_::AbstractRemoteSegments,
                 tuner_::AbstractTuner,tunerstate_::AbstractTunerState,indicator_::AbstractIndicatorMatrix,chain_::AbstractChain)
    println("==============")
    println("Burnin started")
    println("==============")
    showsamplerstatevars(indicatorstate_,"indicator state")
    println("==================================")
    storeduring = _storeduring(:burnin,runner_.policy)
    initialize!(runner_,model_,indicatorstate_,chain_,storeduring)
    for i=1:runner_.numburnin
        iterate!(runner_,model_,indicatorstate_,segments_,indicator_,chain_,storeduring)
        accepted!(tunerstate_,indicator_)
        needstuning(tuner_,i)?tune!(runner_,sampler_,indicatorstate_,segments_,tuner_,tunerstate_):nothing
        indicatorstate_ = preparenext!(runner_,model_,indicatorstate_,segments_,indicator_) #prepare for the next iteration
    end
    println("================")
    println("Burnin completed")
    println("================")
    showsamplerstatevars(indicatorstate_,"indicator state")
    println("==================================")
    indicatorstate_
end

function main!(runner_::GMHRunner,model_::AbstractModel,sampler_::AbstractSampler,indicatorstate_::AbstractSamplerState,segments_::AbstractRemoteSegments,
               tuner_::AbstractTuner,tunerstate_::AbstractTunerState,indicator_::AbstractIndicatorMatrix,chain_::AbstractChain)
    println("================")
    println("Main run started")
    println("================")
    showsamplerstatevars(indicatorstate_,"indicator state")
    println("==================================")
    storeduring = _storeduring(:main,runner_.policy)
    for i=1:runner_.numiterations
        iterate!(runner_,model_,indicatorstate_,segments_,indicator_,chain_,storeduring)
        needstuning(tuner_,i)?println("Iteration $i/$(runner_.numiterations)"):nothing
        indicatorstate_ = preparenext!(runner_,model_,indicatorstate_,segments_,indicator_) #prepare for the next iteration
    end
    println("==================")
    println("Main run completed")
    println("==================")
    showsamplerstatevars(indicatorstate_,"indicator state")
    println("==================================")
    indicatorstate_
end

function iterate!(runner_::GMHRunner,model_::AbstractModel,indicatorstate_::AbstractSamplerState,segments_::AbstractRemoteSegments,
                  indicator_::AbstractIndicatorMatrix,chain_::AbstractChain,storeduring::Bool)
    auxiliary!(runner_,model_,indicatorstate_)
    indicatoracceptance = acceptance!(indicatorstate_)
    segmentacceptances = iterate!(segments_,indicatorstate_)
    transitionprobability!(indicator_,indicatoracceptance,segmentacceptances)
    sampleindicator!(indicator_)
    storeduring?store!(runner_,indicatorstate_,segments_,indicator_,chain_):nothing
end

function preparenext!(runner_::GMHRunner,model_::AbstractModel,indicatorstate_::AbstractSamplerState,
                  segments_::AbstractRemoteSegments,indicator_::AbstractIndicatorMatrix)
    indicatorend = indicatorsamples(indicator_)[end]
    if indicatorend != numproposals(indicator_) + 1
        indicatorstate_ = prepare!(segments_,indicatorstate_,indicatorend)
    else
        indicatorstate_ = prepareindicator!(indicatorstate_)
    end
    traitvalue(runner_.policy.model)==:stochastic?geometry!(model_,from(indicatorstate_)):nothing
    indicatorstate_
end

function auxiliary!(runner_::GMHRunner,model_::AbstractModel,indicatorstate_::AbstractSamplerState)
    auxcounter = 0
    while ~isfinite((propose!(indicatorstate_) ; geometry!(model_,proposals(indicatorstate_))).logprior[1])
        if (auxcounter+=1) > maxinitialcounter
            error("Problems generating auxiliary proposal, giving up. Please check the parameter priors.")
        end
    end
    indicatorstate_
end

function store!(runner_::GMHRunner,indicatorstate_::AbstractSamplerState,segments_::AbstractRemoteSegments,
                indicator_::AbstractIndicatorMatrix,chain_::AbstractChain)
    nproposals = numproposals(indicator_)
    nsamples = numsamples(indicator_)
    indicatorsamples_ = indicatorsamples(indicator_)
    retrievesamples!(segments_,indicatorsamples_)
    from_ = from(indicatorstate_)
    for i=2:nsamples+1
        @inbounds if indicatorsamples_[i] == nproposals + 1
            store!(chain_,from_,1)
        else
            @inbounds store!(segments_,chain_,indicatorsamples_[i])
        end
    end
    accepted!(chain_,indicator_)
end

function tune!(runner_::AbstractRunner,sampler_::AbstractSampler,indicatorstate_::AbstractSamplerState,
               segments_::AbstractRemoteSegments,tuner_::AbstractTuner,tunerstate_::AbstractTunerState)
    tvals = tune(tuner_,tunerstate_)
    showstep(tuner_,tunerstate_)
    tune!(indicatorstate_,tvals...)
    tune!(segments_,tvals...)
    nextindex!(tunerstate_)
end

function show(io::IO,r::GMHRunner)
    println(io,"Generalized Metropolis-Hastings runner with:")
    println(io," numburnin: ",r.numburnin)
    println(io," numiterations: ",r.numiterations)
    println(io," numproposals: ",r.numproposals)
    println(io," numindicatorsamples: ",r.numindicatorsamples)
    print(io," policy: ")
    show(io,r.policy)
    println(io)
    nothing
end


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

using UpROOT
using Test

using ArraysOfArrays
using TypedTables

@testset "testfiles" begin
    @testset "leaflist" begin
        file = TFile(UpROOT.testfiles["leaflist"])
        @test "tree;1" in keys(file)
        tree = file["tree"]

        @test tree[1] isa NamedTuple
        @test tree[1].leaflist isa NamedTuple
        @test tree[1:5] isa Table
        @test tree[1:5].leaflist isa Table
        @test tree[:] isa Table
        @test tree[:].leaflist isa Table

        @test tree.leaflist[1] isa NamedTuple
        @test tree.leaflist[1:5] isa Table
        @test tree.leaflist[:] isa Table

        @test tree[1].leaflist == tree.leaflist[1]
    end

    @testset "HZZ" begin
        file = TFile(UpROOT.testfiles["HZZ"])
        @test "events;1" in keys(file)
        tree = file["events"]

        @test tree[1] isa NamedTuple
        @test tree[1:5] isa Table
        @test tree[:] isa Table

        @test tree.Jet_E[1] isa AbstractVector
        @test tree.Jet_E[1:5] isa VectorOfVectors
        @test tree.Jet_E[:] isa VectorOfVectors
    end
end


text:
# Example from "Bounds on the reachable sets of nonlinear control systems" Constant Temperature Profile
id = :ScottB2013a

pL = [420.0]
pU = [480.0]
ScottB2013a_x0(p::Vector{T}) where T = [1.5*one(T); 0.5*one(T); zero(T)]
function ScottB2013a_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T

    R = 8.134
    A1 = 2400
    A2 = 8800
    E1 = 6.9*10^3
    E2 = 1.69*10^4

    du[1] = -A1*exp(E1/(R*p[1]))*u[1]
    du[2] = A1*exp(E1/(R*p[1]))*u[1] - A2*exp(E2/(R*p[1]))*u[2]
    du[3] = A2*exp(E2/(R*p[1]))*u[2]
    return
end
tspan = (0.0, 0.08)

prob = ODERelaxProb(ScottB2013a_f!, tspan, ScottB2013a_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/pii/S0005109812004839"
source = "ScottB2013"
desc = "3 Series Rxn Network"


text:
#=
    Cheating model.

    From
    https://discourse.julialang.org/t/fitting-a-observed-value-to-a-binomial-distribution-turing/66619
    """
    Fitting a observed value to a Binomial Distribution Turing

    I am new to Turing and trying to learn it by trying to replicate the Chapters from the book
    https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers 2 .

    Here is my problem (from Chapter 2) :

    1.) Have a probability p ~ Uniform(0,1) which I have to infer the “number_of_cheaters”
    2.) p_skewed = p*(0.5) + (0.25) (deterministic value)
    3.) with model:
    yes_responses = pm.Binomial(“number_of_cheaters”, 100, p_skewed, observed=35)

    How do I write this in Turing ?

    ...

    @model function model(yes_responses,N)
        p ~ Uniform(0,1)
        p_skewed = p*(0.5)+0.25
        rv_yes_responses ~ Binomial(N,p_skewed)

        // This line on how to include yes_responses i am unable to figure out

        return p_skewed

    end
   """

   The model refered to is in
   https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter2_MorePyMC/Ch2_MorePyMC_PyMC3.ipynb

    Output of this Turing.jl model:
        Summary Statistics
    parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

               p    0.2042    0.0993     0.0031    0.0040   446.4728    0.9990      394.4106
        p_skewed    0.3521    0.0497     0.0016    0.0020   446.4728    0.9990      394.4106

    Quantiles
    parameters      2.5%     25.0%     50.0%     75.0%     97.5%
        Symbol   Float64   Float64   Float64   Float64   Float64

               p    0.0349    0.1325    0.2020    0.2749    0.3981
        p_skewed    0.2674    0.3162    0.3510    0.3875    0.4490

    Summary Stats:
    Length:         1000
    Missing Count:  0
    Mean:           0.352105
    Minimum:        0.250395
    1st Quartile:   0.316228
    Median:         0.350993
    3rd Quartile:   0.387453
    Maximum:        0.494718

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function cheating_model2(yes_responses,N)
    p ~ Uniform(0,1)
    # Using Dirac(.) makes p_skewed in the chain
    p_skewed ~ Dirac(p*(0.5)+0.25)
    # p_skewed = p*(0.5)+0.25
    yes_responses ~ Binomial(N,p_skewed)

    # An alternative is to return the variable (see below how to handle this)
    return p_skewed
end

yes_response = 35
N = 100
model = cheating_model2(yes_response,N)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5),  1_000)
# chns = sample(model, PG(5),  MCMCThreads(), 10_000, 4)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)
# display(plot(chns))

chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
genq = generated_quantities(model, chains_params)
display(summarystats(vcat(genq...)))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
##file to include in ComplexModeler.jl which contains all files associated to Matpower

include("MatpowerGenerator.jl")
include("MatpowerLoad.jl")
include("MatpowerROPFShunt.jl")
include("MatpowerLine_π_Imax.jl")
include("MatpowerVolt.jl")
include("read_matpower_simple_Imax.jl")


text:
"""
    NonlinearKoopman(K, B, C, basis, Q, P, discrete)

An approximation of the Koopman operator which is nonlinear in the states.

`K` is the array representing the operator, `B` is the (possible present) array
representing the influence of exogenous inputs on the evolution.
`C` is the array mapping from the Koopman space to the original state space. `basis` is a
[Basis]@ref(Basis), mapping the state space to the Koopman space.

`Q` and `P` are matrices used for updating the operator with new measurements.
`discrete` indicates if the operator is discrete or continuous.

The Koopman operator is callable with the typical signature of `f(u,p,t)` and `f(du,u,p,t)`, respectively.

# Example

```julia
k = EDMD(X, basis)

u = k([2.0; 0.5], nothing, nothing)
du = similar(u)
k(du, u, nothing, nothing)
```
"""
mutable struct NonlinearKoopman <: AbstractKoopmanOperator
    operator::AbstractArray
    input::AbstractArray
    output::AbstractArray

    basis::Basis

    Q::AbstractArray
    P::AbstractArray

    discrete::Bool
end

(k::NonlinearKoopman)(u, p::DiffEqBase.NullParameters, t) = k(u, [], t)
(k::NonlinearKoopman)(du, u, p::DiffEqBase.NullParameters, t) = k(du, u, [], t)

function (k::NonlinearKoopman)(u,  p::AbstractArray = [], t = nothing)
    return k.output*k.operator*k.basis(u, p, t)
end

function (k::NonlinearKoopman)(du, u, p::AbstractArray = [], t = nothing)
    mul!(du, k.output*k.operator, k.basis(u, p, t))
end

"""
    update!(k, X, Y; p = [], t = [], threshold = eps())

Update the Koopman `k` given new data `X` and `Y`. The operator is updated in place if
the L2 error of the prediction exceeds the `threshold`.

`p` and `t` are the parameters of the basis and the vector of timepoints, if necessary.
"""
function update!(k::NonlinearKoopman, X::AbstractArray, Y::AbstractArray; p::AbstractArray = [], t::AbstractVector = [], threshold::T = eps()) where {T <: Real}
    @assert updatable(k) "Linear Koopman is not updatable."

    Ψ₀ = k.basis(X, p, t)
    Ψ₁ = k.basis(Y, p, t)

    ϵ = norm(Ψ₁-operator(k)*Ψ₀, 2)

    if ϵ < threshold
        return
    end

    k.Q += Ψ₁*Ψ₀'
    k.P += Ψ₀*Ψ₀'
    k.operator .= k.Q / k.P

    if norm(Y - outputmap(k)*k.operator*Ψ₀) < threshold
        return
    end

    # TODO Make this a proper rank 1 update
    k.output .= X / Ψ₀

    return
end

"""
    reduce_basis(k; threshold)

Reduces the `basis` of the nonlinear Koopman using the 1-norm of each row
of the matrix `C*K`. Rows where the threshold is not reached are deleted.
"""
function reduce_basis(k::NonlinearKoopman; threshold = 1e-5, kwargs...)
    b = k.output*k.operator
    inds = vec(sum(abs, b, dims = 1) .> threshold)
    return Basis(map(x->x.rhs, k.basis[inds]), variables(k.basis), parameters = parameters(k.basis), iv = independent_variable(k.basis), kwargs...)
end

"""
    ODESystem(k; threshold = eps())

Convert a `NonlinearKoopman` into an `ODESystem`. `threshold` determines the cutoff
for the entries of the matrix representing the state space evolution of the system.
"""
function ModelingToolkit.ODESystem(k::NonlinearKoopman; threshold = eps(), kwargs...)
    @assert threshold > zero(threshold) "Threshold must be greater than zero"

    eqs = Any[]
    A = outputmap(k)*k.operator
    A[abs.(A) .< threshold] .= zero(eltype(A))
    eqs = hcat([x.rhs for x in k.basis.eqs])
    b = Basis(simplify.(A*eqs)[:,1], variables(k.basis), parameters = parameters(k.basis), iv = independent_variable(k.basis))
    return ODESystem(b, kwargs...)
end


text:
# ---
# title: 909. Snakes and Ladders
# id: problem909
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/snakes-and-ladders/description/>
# hidden: true
# ---
#
# On an N x N `board`, the numbers from `1` to `N*N` are written
# _boustrophedonically_   **starting from the bottom  left of the board**, and
# alternating direction each row.  For example, for a 6 x 6 board, the numbers
# are written as follows:
#
#
#
#     ![](https://assets.leetcode.com/uploads/2018/09/23/snakes.png)
#
#
# You start on square `1` of the board (which is always in the last row and
# first column).  Each move, starting from square `x`, consists of the
# following:
#
#   * You choose a destination square `S` with number `x+1`, `x+2`, `x+3`, `x+4`, `x+5`, or `x+6`, provided this number is `<= N*N`.
#     * (This choice simulates the result of a standard 6-sided die roll: ie., there are always **at most 6 destinations, regardless of the size of the board**.)
#   * If `S` has a snake or ladder, you move to the destination of that snake or ladder.  Otherwise, you move to `S`.
#
# A board square on row `r` and column `c` has a "snake or ladder" if
# `board[r][c] != -1`.  The destination of that snake or ladder is
# `board[r][c]`.
#
# Note that you only take a snake or ladder at most once per move: if the
# destination to a snake or ladder is the start of another snake or ladder, you
# do **not** continue moving.  (For example, if the board is `[[4,-1],[-1,3]]`,
# and on the first move your destination square is `2`, then you finish your
# first move at `3`, because you do **not** continue moving to `4`.)
#
# Return the least number of moves required to reach square N*N.  If it is not
# possible, return `-1`.
#
# **Example 1:**
#
#
#
#     Input: [
#     [-1,-1,-1,-1,-1,-1],
#     [-1,-1,-1,-1,-1,-1],
#     [-1,-1,-1,-1,-1,-1],
#     [-1,35,-1,-1,13,-1],
#     [-1,-1,-1,-1,-1,-1],
#     [-1,15,-1,-1,-1,-1]]
#     Output: 4
#     Explanation:
#     At the beginning, you start at square 1 [at row 5, column 0].
#     You decide to move to square 2, and must take the ladder to square 15.
#     You then decide to move to square 17 (row 3, column 5), and must take the snake to square 13.
#     You then decide to move to square 14, and must take the ladder to square 35.
#     You then decide to move to square 36, ending the game.
#     It can be shown that you need at least 4 moves to reach the N*N-th square, so the answer is 4.
#
#
# **Note:**
#
#   1. `2 <= board.length = board[0].length <= 20`
#   2. `board[i][j]` is between `1` and `N*N` or is equal to `-1`.
#   3. The board square with number `1` has no snake or ladder.
#   4. The board square with number `N*N` has no snake or ladder.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
get_union_type_2(::Type{Union{Missing, T}}) where T = T

function create_time(time_str, ampm, tz)
    if ismissing(time_str) || ismissing(ampm)
        return missing
    end
    dt = DateTime(time_str, dateformat"yyyy/mm/dd H:MM")
    dt += hour(dt) == 12 ? (ampm == "AM" ? Hour(12) : Hour(0)) : (ampm == "PM" ? Hour(12) : Hour(0))
    ZonedDateTime(dt, tz)
end


text:
module BayesTesting

using Distributions, KernelDensity, LinearAlgebra, StatsBase, Plots, StatsPlots
# using PlotRecipes
export update_mean, marginal_posterior_mu, mcodds, todds, bayespval, pdr_pval
export blinreg, bayesregNIG, bayesreg, gsreg, compare_means, compare_proportions
export beta_posterior, beta_update, equiv_test, Bayesian_ttest, correlation_ttest
# package code goes here
include("bayes_meta_analysis.jl")
include("bayes_odds_pvalue.jl")
include("bayesreg.jl")
include("gsreg.jl")
include("compare_means.jl")
include("bayesian_ttest.jl")
include("correlation_coefficient_v1.jl")
end # module


text:
"""
    ComputeAdimensionalHubbleFactor(z::Float64, w0waCDMCosmology::w0waCDMCosmology)

This function, given the value of the cosmological parameters, evaluate the
Adimensional Hubble Factor for ``w_0 w_a``CDM cosmologies.
The analitycal expression is given by:
```math
E(z)=\\sqrt{\\Omega_M(1+z)^3+\\Omega_R(1+z)^4+
\\Omega_{DE}(1+z)^{3(1+w_0+w_a)}\\exp(-3w_a \\frac{z}{1+z})+\\Omega_k(1+z)^2}
```

!!! warning
    This expression is valid only for the
    [CPL parameterization](https://arxiv.org/abs/astro-ph/0208512)
    of the Dark Energy Equation of State.
"""
function ComputeAdimensionalHubbleFactor(z::Float64, w0waCDMCosmology::w0waCDMCosmology)
    E_z = sqrt(w0waCDMCosmology.ΩM*(1+z)^3 + w0waCDMCosmology.Ωr*(1+z)^4+
    w0waCDMCosmology.Ωk*(1+z)^2
    +w0waCDMCosmology.ΩDE*(1+z)^(3*(1+w0waCDMCosmology.w0+
    w0waCDMCosmology.wa))*exp(-3*w0waCDMCosmology.wa*z/(1+z)) )
    return E_z
end

function ComputeAdimensionalHubbleFactor(z::Float64,
    flatw0waCDMCosmology::Flatw0waCDMCosmology)
    E_z = sqrt(flatw0waCDMCosmology.ΩM*(1+z)^3
    +(1-flatw0waCDMCosmology.ΩM)*(1+z)^(3*(1+flatw0waCDMCosmology.w0+
    flatw0waCDMCosmology.wa))*exp(-3*flatw0waCDMCosmology.wa*z/(1+z)) )
    return E_z
end

"""
    ComputeHubbleFactor(z::Float64, AbstractCosmology::AbstractCosmology)

This function, given the value of the cosmological parameters, evaluate the
Hubble Factor for ``w_0 w_a``CDM cosmologies, whose expression is given by
```math
H(z)=H_0\\sqrt{\\Omega_M(1+z)^3+\\Omega_R(1+z)^4+
\\Omega_{DE}(1+z)^{3(1+w_0+w_a)}\\exp(-3w_a \\frac{z}{1+z})+\\Omega_k(1+z)^2}
```

"""
function ComputeHubbleFactor(z::Float64, AbstractCosmology::AbstractCosmology)
    H_z = AbstractCosmology.H0*ComputeAdimensionalHubbleFactor(z,
    AbstractCosmology)
    return H_z
end

"""
    Computeχ(z::Float64, AbstractCosmology::AbstractCosmology)

This function, given the value of the cosmological parameters, evaluate the
Comoving Distance χ. It is evaluated as:
```math
\\chi(z)=\\frac{c}{H_0}\\int_0^z \\frac{dz'}{E(z')}
```
"""
function Computeχ(z::Float64,
    AbstractCosmology::AbstractCosmology)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    integral, err = QuadGK.quadgk(x -> 1 /
    ComputeAdimensionalHubbleFactor(x,AbstractCosmology), 0, z, rtol=1e-12)
    return integral*c_0/AbstractCosmology.H0
end

"""
    ComputeBackgroundQuantitiesOverGrid(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    AbstractCosmology::AbstractCosmology)

This function evaluate the Hubble factor and the comoving distance over the
[`CosmologicalGrid`](@ref).
"""
function ComputeBackgroundQuantitiesGrid!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    AbstractCosmology::AbstractCosmology)
    for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
        BackgroundQuantities.HZArray[idx_ZArray] = ComputeHubbleFactor(
        CosmologicalGrid.ZArray[idx_ZArray], AbstractCosmology)
        BackgroundQuantities.rZArray[idx_ZArray] = Computeχ(
        CosmologicalGrid.ZArray[idx_ZArray], AbstractCosmology)
    end
end

function ComputeLogSpacedχGrid!(CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities)
    Zχ = Dierckx.Spline1D(BackgroundQuantities.rZArray, CosmologicalGrid.ZArray)
    χArray = CosmoCentral.LogSpaced(BackgroundQuantities.rZArray[1],
    last(BackgroundQuantities.rZArray), length(BackgroundQuantities.rZArray))
    CosmologicalGrid.ZArray = Zχ.(χArray)
end

function ComputeBackgroundQuantitiesOverLogSpacedχGrid(
    CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    AbstractCosmology::AbstractCosmology)
    ComputeBackgroundQuantitiesGrid!(CosmologicalGrid, BackgroundQuantities,
    AbstractCosmology)
    ComputeLogSpacedχGrid!(CosmologicalGrid, BackgroundQuantities)
    ComputeBackgroundQuantitiesGrid!(CosmologicalGrid, BackgroundQuantities,
    AbstractCosmology)
end

function ExtractGrowthFactor!(BackgroundQuantities::BackgroundQuantities,
    PowerSpectrum::PowerSpectrum)
    BackgroundQuantities.DZArray = zeros(length(PowerSpectrum.PowerSpectrumLinArray[1,:]))
    for zidx in 1:length(BackgroundQuantities.DZArray)
        BackgroundQuantities.DZArray[zidx] =
        sqrt(PowerSpectrum.PowerSpectrumLinArray[1,zidx]./
        PowerSpectrum.PowerSpectrumLinArray[1,1])
    end
end


text:
module External
using  ArgParse
using Logging
import Dates
import ..Gadgetron

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table! s begin
	"port"
		help="Port to use for connection"
		required=true
		arg_type = Int
	"module"
		help="Either the name of a module or a Julia file"
		required = true
	"target"
		help="Function name to call with connection (module only)"
		required= false
	end
	return parse_args(s; as_symbols=true)
end

function load_function(module_name, target)

	is_julia_file = module_name[end-2:end] == ".jl"
	if is_julia_file
		throw(ErrorException("Using julia files directly is currently unsupported"))
	end
	loaded_module = Base.require(Main,Symbol(module_name))
	func_expr = Expr(:.,loaded_module,QuoteNode(Symbol(target)))
	@debug "Created expression $func_expr"

	return Meta.eval(func_expr)
end

function log_formatter(level, _module , group, id, file, line)
	timestamp = Dates.format(Dates.now(),"m-d H:M:S.s")
	prefix = "$timestamp $level [$file:$line] "
	suffix = ""

	return (:default, prefix, suffix )
end

function main()
	args = parse_commandline()

	Base.global_logger(ConsoleLogger(stdout,Logging.Info ))

	@info  "Starting external Julia module $(args[:module]) in state: [ACTIVE]"
	@info "Connection to parent on port $(args[:port])"

	func = load_function(args[:module],args[:target])
	connection = Gadgetron.connect("localhost",args[:port])

	try
		Base.invokelatest(func,connection)
	finally
		close(connection)
	end

end

Base.precompile(Tuple{typeof(main)})
Base.precompile(Tuple{typeof(load_function),String,String})   # time: 0.10844494

end


text:
using CairoMakie
using DelimitedFiles

QX = readdlm("QX.dat")
QY = readdlm("QY.dat")
FT_real = readdlm("FT_real.dat")
FT_imag = readdlm("FT_imag.dat")

FT_res = FT_real + 1im .* FT_imag

fig = CairoMakie.Figure(resolution = (1800, 1800))
ax =
        fig[1, 1] = Axis(
                fig,
                xlabel = "qx",
                ylabel = "qy",
                xlabelpadding = 0,
                ylabelpadding = 0,
                xlabelsize = 12,
                ylabelsize = 12,
                xticklabelsize = 12,
                yticklabelsize = 12,
                aspect = AxisAspect(1),
                xticklabelfont = "serif-roman",
                yticklabelfont = "serif-roman",
                xlabelfont = "serif-italic",
                ylabelfont = "serif-italic",
        )

hm = heatmap!(qx, qy, angle.(FT_res), colormap = :oslo)
Colorbar(fig[1, 2], hm, width = 40)
fig
save("FT_phase.pdf", fig)

fig2 = CairoMakie.Figure(resolution = (1800, 1800))
ax2 =
        fig2[1, 1] = Axis(
                fig2,
                xlabel = "qx",
                ylabel = "qy",
                xlabelpadding = 0,
                ylabelpadding = 0,
                xlabelsize = 12,
                ylabelsize = 12,
                xticklabelsize = 12,
                yticklabelsize = 12,
                aspect = AxisAspect(1),
                xticklabelfont = "serif-roman",
                yticklabelfont = "serif-roman",
                xlabelfont = "serif-italic",
                ylabelfont = "serif-italic",
        )

hm2 = heatmap!(qx, qy, abs.(FT_res), colormap = :oslo)
Colorbar(fig2[1, 2], hm2, width = 40)
fig2
save("FT_amplitude.pdf", fig)


text:
using MLStyle
using MLStyle.MatchCore: mangle
using MLStyle.Infras: @format

def_pattern(Main,
    predicate = (@λ begin
        :[$_ for $_ in $_] -> true
        _ -> false
    end),
    rewrite = (tag, case, mod) -> begin
        @match case begin
            :[$expr for $iter in $seq] => begin
                function (body)
                    iter_var = mangle(mod)
                    produced_elt_var = mangle(mod)
                    inner_test_var = mangle(mod)
                    seq_var = mangle(mod)
                    decons_elt = mk_pattern(iter_var, expr, mod)(iter)
                    decons_seq = mk_pattern(seq_var, seq, mod)(body)
                    @format [
                        seq_var, iter_var, inner_test_var, produced_elt_var,
                        body, tag, decons_elt, decons_seq, push!] quote
                        if tag isa Vector
                            let seq_var = [], inner_test_var = true
                                for iter_var in tag
                                    produced_elt_var = decons_elt
                                    if produced_elt_var === failed
                                        inner_test_var = false
                                        break
                                    end
                                    push!(seq_var, produced_elt_var)
                                end
                                if inner_test_var
                                    decons_seq
                                else
                                    failed
                                end
                            end
                        else
                            failed
                        end
                    end
                end
            end
        end
    end
)

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: thautwarm/MLStyle-Playground

max_forks_repo_name: thautwarm/MLStyle-Playground


text:
abstract type Hexacopter <: Multicopter end

include("LeeHexacopter.jl")


text:
abstract type StartValueGenerator end

const two_pi = convert(Float64, 2π)

struct TorusStartValueGenerator{N} <: StartValueGenerator
end

startvalue(::TorusStartValueGenerator{2}, w=nothing) = SVector(rand() * two_pi, rand() * two_pi)
startvalue(::TorusStartValueGenerator{3}, w=nothing) = SVector(rand() * two_pi, rand() * two_pi, rand() * two_pi)

struct DomainStartValueGenerator2D <: StartValueGenerator
    xmin::Float64
    ymin::Float64

    dx::Float64
    dy::Float64
end

function DomainStartValueGenerator2D(lims::Tuple{<:Real, <:Real, <:Real, <:Real})
    xmin, xmax, ymin, ymax = convert.(Float64, lims)
    DomainStartValueGenerator2D(xmin, ymin, xmax - xmin, ymax - ymin)
end
DomainStartValueGenerator2D(xlims, ylims) = DomainStartValueGenerator2D((xlims..., ylims...))

# function DomainStartValueGenerator2D(f::MP.AbstractPolynomial; apply_exp=false, kwargs...)
#     xmin, xmax, ymin, ymax = amoeba_carcase_domain_heuristic(f; kwargs...)
#     DomainStartValueGenerator2D(xmin, ymin, xmax - xmin, xmax - ymin)
# end

function startvalue(c::DomainStartValueGenerator2D, y=nothing)
    w1 = rand() * c.dx + c.xmin
    w2 = rand() * c.dy + c.ymin

    SVector(w1, w2)
end

struct DomainStartValueGenerator3D <: StartValueGenerator
    xmin::Float64
    ymin::Float64
    zmin::Float64

    dx::Float64
    dy::Float64
    dz::Float64
end

function DomainStartValueGenerator3D(lims::Tuple{<:Real, <:Real, <:Real, <:Real, <:Real, <:Real})
    xmin, xmax, ymin, ymax, zmin, zmax = convert.(Float64, lims)
    DomainStartValueGenerator3D(
        xmin, ymin, zmin,
        xmax - xmin, ymax - ymin, zmax - zmin)
end
DomainStartValueGenerator3D(xlims, ylims, zlims) = DomainStartValueGenerator3D(xlims..., ylims..., zlims...)

function startvalue(c::DomainStartValueGenerator3D, w=nothing)
    w1 = rand() * c.dx + c.xmin
    w2 = rand() * c.dy + c.ymin
    w3 = rand() * c.dz + c.zmin

    SVector(w1, w2, w3)
end

startvalue_generator(::AmoebaFiber2D) = TorusStartValueGenerator{2}()
startvalue_generator(::AmoebaFiber3D) = TorusStartValueGenerator{3}()


text:
using Documenter, SeisML

makedocs(;
    modules=[SeisML],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tclements/SeisML.jl/blob/{commit}{path}#L{line}",
    sitename="SeisML.jl",
    authors="Tim Clements",
    assets=String[],
)

deploydocs(;
    repo="github.com/tclements/SeisML.jl",
)


text:
#=
test_src:
- Author: Berk
- Date: 2020-06-16
This tests lse approximations.
=#

using PyCall

function compile_lse_constraints(n_samples = 300)
    sagemarks = pyimport("sagebenchmarks.literature.solved");
    ineq_constraints = BlackBoxLearner[]
    for idx=1:25
        md = OCTHaGOn.sagemark_to_GlobalModel(idx, false);
        update_bounds!(md, lbs = Dict(md.vks[end] => -200), ubs= Dict(md.vks[end] => 200))
        if !any(isinf.(values(md.lbs))) && !any(isinf.(values(md.ubs)))
            println("Adding constraints from test ", string(idx));
            for i=1:length(md.bbls)
                set_param(md.bbls[i], :n_samples, n_samples)
                md.bbls[i].name = Float64(idx) + i/10.;
                push!(ineq_constraints, md.bbls[i]);
            end
        end
    end
    println("Compiled all multivariate LSE inequalities!")
    return ineq_constraints
end

ineqs = compile_lse_constraints();
uniform_sample_and_eval!(ineqs);

learn_constraint!(ineqs)

# Plotting all accuracies (non-convex approximation)
ineqs = ineqs[findall(x -> length(x.accuracies) > 0, ineqs)]
plot_accuracies(ineqs)

convex_idxs = findall(x -> "logconvex" in x.tags, ineqs);
nonconvex_idxs = findall(x -> x ∉ convex_idxs, [1:length(ineqs)]);

plot_accuracies(ineqs[convex_idxs])
plot_accuracies(ineqs[nonconvex_idxs])

# nonconvex_scores = bar!(nonconvex_idxs, [ineq.accuracies[end] for ineq in ineqs[nonconvex_idxs]],
#                   xlabel="Constraint number",
#                         title="Accuracy over all constraints", xticks=1:length(ineqs),
#                         label="non-convex", legend=:topright, ylims=[0.9, 1.])

# # Comparing accuracy vs. convexity
# data = Dict(
#     "x" => ["convex", "non-convex"],
#     "y" => [sum(scores.*convex)/sum(convex),
#             sum(scores.*(1 .- convex))/(length(convex) - sum(convex))],
# );
# convexity_effect = bar(data["x"], data["y"], title="Accuracy vs. convexity", legend=false,ylims=[0.9, 1.])
#
# # Retrain non-convex with different numbers of samples
# n_samples = [100, 250, 500, 750, 1000];
# for idx in problem_idxs
#     md = OCTHaGOn.sagemark_to_GlobalModel(problem_idxs[idx], false);
#     n_dims = length(md.lbs);
#     for j in n_samples
#         plan, _ = LHCoptim(j, n_dims, 1);
#         X = scaleLHC(plan,[(md.lbs[i], md.ubs[i]) for i=1:n_dims]);
#         feas_tree = learn_constraints!(lnr, [ineqs[idx]], X);
#     IAI.write_json(string(PROJECT_ROOT, "/test/data/constraint",idx,"_",j,"samples",".tree"), feas_tree[1]);
#     end
# end
#
# # Loading trees, and scoring w.r.t. maximum of the samples
# n_samples = [100, 250, 500, 750, 1000];
# grids = [IAI.read_json(string("data/constraint",idx,"_",j,"samples",".tree")) for j in n_samples for idx=1:length(ineqs)]; # inside for loop gets executed first.
# scores = [];
# for idx=1:length(ineqs)
#     md = OCTHaGOn.sagemark_to_GlobalModel(problem_idxs[idx], false);
#     n_dims = length(md.lbs);
#     plan, _ = LHCoptim(maximum(n_samples), n_dims, 1);
#     X = scaleLHC(plan,[(md.lbs[i], md.ubs[i]) for i=1:n_dims]);
#     Y = [ineqs[idx](X[j,:]) >= 0 for j = 1:maximum(n_samples)];
#     score_list = [];
#     for j in n_samples
#         nlnr = IAI.read_json(string("data/constraint",idx,"_",j,"samples",".tree"));
#         append!(score_list, IAI.score(nlnr, X, Y))
#     end
#     append!(scores, score_list)
# end
# scores = transpose(reshape(scores, (length(n_samples), length(ineqs))));
#
# # Finally plot
# for i=1:length(convex_idxs)
#     display(plot!(n_samples, scores[i,:], xlabel="Number of samples", ylabel="Accuracy",
#                   legend=false, color=:blue))
# end
# # Plotting trend
# display(plot!(n_samples, transpose(sum(scores, dims=1))./size(scores,1), xlabel="Number of samples", ylabel="Accuracy",
#                   legend=false,linewidth=4, color=:blue))
#
# # Dimension vs. accuracy
# plot(xlabel="Number of samples", ylabel="Accuracy")
# for i in unique(dimension)
#     dim_idxs = findall(x -> x==i, dimension);
#     y = sum(scores[dim_idxs, :], dims=1)./size(dim_idxs);
#     display(plot!(n_samples, y', label=string("Dimension ",i)))
# end


text:
using Nettle
using Random: AbstractRNG

using CryptoUtils: is_quadratic_residue, sqrt_mod_prime

struct Hash
    spec::String
end

(h::Hash)(x::Vector{UInt8}) = hex2bytes(hexdigest(h.spec, x))
(h::Hash)(x::String) = h(Vector{UInt8}(x))

# Dispatching on value types seems as plausable solution
function bitlength(h::Hash)
    s = h.spec

    if s == "sha256"
        return 256
    elseif s == "sha384"
        return 384
    elseif s == "sha512"
        return 512
    else
        error("No corepsonding mapping for $x implemented")
    end
end

struct PRG <: AbstractRNG
    h::Hash
    s::Vector{UInt8}
end

(prg::PRG)(i::UInt32) = prg.h([prg.s..., reverse(reinterpret(UInt8, UInt32[i]))...])

function Base.getindex(prg::PRG, range)
    (; start, stop) = range

    a = bitlength(prg.h) ÷ 8 # outlength

    K = div(stop, a, RoundUp) - 1

    r = UInt8[]

    for i in UInt32(0):UInt32(K)
        ri = prg(i)
        append!(r, ri)
    end

    return r[range]
end

struct RO
    h::Hash
    n_out::Int
end

zerofirst(x, n) = (x << n) >> n # Puts first n bits of a number x to zero.

function (ro::RO)(d::Vector{UInt8})
    (; h, n_out) = ro

    nb = reinterpret(UInt8, UInt32[n_out])
    s = h([reverse(nb)...,d...]) # Numbers on Java are represented in reverse
    prg = PRG(h, s)

    a = prg[1:div(n_out, 8, RoundUp)]

    if mod(n_out, 8) != 0
        a[1] = zerofirst(a[1], 8 - mod(n_out, 8))
    end

    return a
end

_tobig(x) = parse(BigInt, bytes2hex(reverse(x)), base=16)
interpret(::Type{BigInt}, x::Vector{UInt8}) = _tobig(reverse(x))

function interpret(::Type{Vector{T}}, 𝐫::Vector{UInt8}, N::Int) where T <: Integer
    M = length(𝐫) ÷ N
    𝐮 = reshape(𝐫, (M, N))
    𝐭 = [interpret(T, 𝐮[:, i]) for i in 1:N]
    return 𝐭
end

function Base.rand(prg::PRG, ::Type{T}, N::Int; n = bitlength(T)) where T <: Integer

    M = div(n, 8, RoundUp) # bytes for each number

    total = M * N

    𝐫 = prg[1:total]
    𝐭 = interpret(Vector{BigInt}, 𝐫, N)

    return 𝐭
end

Base.rand(prg::PRG, n::Int, N::Int) = Base.rand(prg, BigInt, N; n)

struct ROPRG
    ρ::Vector{UInt8}
    rohash::Hash
    prghash::Hash
end

function (roprg::ROPRG)(x::Vector{UInt8})

    (; ρ, rohash, prghash) = roprg

    ns = bitlength(prghash) # outlen
    ro = RO(rohash, ns)

    d = UInt8[ρ..., x...]

    s = ro(d)
    prg = PRG(prghash, s)
    return prg
end

(roprg::ROPRG)(x::String) = roprg(Vector{UInt8}(x))
(roprg::ROPRG)(x::Symbol) = roprg(string(x))

function Base.rand(prg::PRG, spec::MODP, N::Integer; nr::Integer = 0)

    p = modulus(spec)
    q = order(spec)

    @assert !isnothing(q) "Order of the group must be known"

    np = bitlength(p)

    𝐭 = rand(prg, BigInt, N; n = np + nr)

    𝐭′ = mod.(𝐭, big(2)^(np + nr))

    𝐡 = powermod.(𝐭′, (p - 1) ÷ q, p)

    return 𝐡
end

function Base.rand(prg::PRG, spec::ECP, N::Integer; nr::Integer = 0)

    (; a, b) = spec

    p = modulus(spec)
    q = order(spec)

    np = bitlength(p) # 1

    𝐭 = rand(prg, BigInt, N*10; n = np + nr)  # OPTIMIZE

    𝐭′ = mod.(𝐭, big(2)^(np + nr))

    𝐳 = mod.(𝐭′, p)

    𝐡 = Vector{Tuple{BigInt, BigInt}}(undef, N)

    l = 1

    f(x) = x^3 + a*x + b # This assumes that I do know how to do arithmetics with fields.

    for zi in 𝐳
        y2 = mod(f(zi), p)

        if is_quadratic_residue(y2, p)

            x = zi
            y = sqrt_mod_prime(y2, p)

            # The smallest root is taken
            if p - y < y
                y = p - y
            end

            𝐡[l] = (x, y)

            if l == N
                break
            else
                l += 1
            end
        end
    end

    if l != N
        error("Not enough numbers for 𝐭 have been allocated")
    end

    return 𝐡
end


text:
__precompile__(false)

module CodeManager

using Distributed

""" TODO """
debug(s...) = println("\x1b[36mDEBUG: ",s...,"\x1b[39m\x1b[49m")

""" TODO """
warn(s...) = println("\x1b[35mWARNING: ",s...,"\x1b[39m\x1b[49m")

const EMPTY_FUNCTION = (args...) -> nothing

""" force garbage collection, free memory """
const force_gc = GC.gc
const free_memory = GC.gc
export force_gc, free_memory

""" TODO """
function print_backtrace(ex::Exception)
  print("\x1b[91m")
  Base.showerror(stderr, ex, catch_backtrace())
  println("\x1b[0m")
end

""" TODO """
function print_backtrace(f::Function)
  try
    return f()
  catch ex
    print_backtrace(ex)
  end
  nothing
end
export print_backtrace

""" TODO """
modulechild(mod::Module, sym::Symbol) = Core.eval(Base.moduleroot(mod), Meta.parse(string(mod)*"."*string(sym)))

""" TODO """
function safe_include(mod::Module, path::String)
  result = (false, nothing)
  #safe_clean!(mod)
  print_backtrace(() -> begin
    result = (true, Core.include(mod, abspath(path)))
  end)
  result
end
export safe_include

#""" TODO """
#function include_module2(mod::Module, sym::Symbol, path::String)
#  result = (false, nothing)
#  safe_clean!(mod, sym)
#  print_backtrace(() -> begin
#    result = (true, Core.eval(mod, :(module $sym; include(abspath($path)); end)))
#  end)
#  result
#end

""" TODO """
function cleanCode(code::String)

	# problem using match(): only "function (name)" will be detected!
	# for x in eachmatch(r"function \s(\w)", code) println(x.captures[1]) end

	# TODO merge code lines
	code=replace(code, "\r" => "")
	#code=replace(code, r"\#\=(?(?=\=\#)then|else)*\=\#", "")
	code=replace(code, r"(?s)(?<=\#\=)(.*?)(?=\=\#)" => "")
	code=replace(code, r"(\#[^\n]+)" => "")
	code=replace(code, r"\n\s+\n" => "\n")
	code=replace(code, r"\n+" => "\n")
  code=replace(code, r"\s*([\+\-\*\/]?\=)\s*" => s"\1")
	code=replace(code, r"(\()\s*" => s"\1")
	code=replace(code, r"\s*(\))" => s"\1")
  code=replace(code, r"\n" => ";")
	code=replace(code, r"\s*([,;])\s*" => s"\1")
end

#code = cleanCode(code)
#code = replace(code, r"^\s*module\s+\w+\s*" => "")
#code = replace(code, r"\s+end\s*$" => "")
#println(code)

#mgr=@__MODULE__
#LF=";"
#add = ""
#add *= "const self = "*name*LF
#add *= "safe_eval(x) = $mgr.safe_eval("*name*", x)"*LF
#add *= "safe_call(x) = $mgr.safe_call("*name*", x)"*LF
#add *= "include_module(x;name=:Nothing) = $mgr.include_module("*name*", x; name=name)"*LF
#add *= "safe_clean!(x) = $mgr.safe_clean!("*name*", x)"*LF
#code = "module "*name*LF*code*LF*"end; println(names("*name*")); "

""" TODO """
function include_module(mod::Module, path::String; name::Symbol=:Nothing)
  result = (false, nothing)

	code = print_backtrace(() -> begin  code = ""; open(abspath(path)) do fp	code = read(fp, String) end; code; end)
	if code != nothing
    find = r"^.*\s*module\s+(\w+)\s*.*$"is
    if occursin(find, code) name = Symbol(replace(code,find=>s"\1")) end
  end

  if name == :Nothing
    print_backtrace(()->error("No module found!"))

  elseif isdefined(mod,name) safe_clean!(mod, name) end
  print_backtrace(() -> begin
    result = (true, Core.include(mod, abspath(path)))
    mgr=@__MODULE__
    child = modulechild(mod, name)

    expr = Expr(:toplevel,
      :(free!() = $mgr.safe_clean!(@__MODULE__)),
      :(safe_invoke(link::Expr, args...) = $mgr.safe_invoke(@__MODULE__, link, args...)),
      :(safe_eval(args...) = $mgr.safe_eval(@__MODULE__, args...)),
      :(safe_call(args...) = $mgr.safe_call(@__MODULE__, args...)),
      :(safe_clean!(args...) = $mgr.safe_clean!(@__MODULE__, args...)),
      :(include_module(path::String;name=:Nothing) = $mgr.include_module(@__MODULE__, path; name=name)),
      )

    Core.eval(child, expr)
  end)

  result
end
export include_module

""" TODO """
safe_invoke(mod::Module, link::Expr, args...) = safe_invoke(safe_eval(mod, link)[2], args...)

""" TODO """
function safe_invoke(f::Function, args...)
  result = (false, nothing)
  print_backtrace(() -> begin
    result = (true, Base.invokelatest(f, args...))
  end)
  result
end
export safe_invoke

""" TODO """
safe_call(mod::Module, value::Union{Symbol, Expr}) = safe_eval(mod, :($value()))

""" TODO """
safe_call(mod::Module, value::Function) = safe_eval(mod, value)
export safe_call

""" TODO """
safe_eval(mod::Module, value::String) = safe_eval(mod, Meta.parse(value))

""" TODO """
safe_eval(mod::Module, value::Function) = safe_eval(mod, :($value()))

""" TODO """
function safe_eval(mod::Module, value::Any)
  result = (false, nothing)
  print_backtrace(() -> begin
    result = (true, Core.eval(mod, value))
  end)
  result
end
export safe_eval

""" TODO """
safe_clean!(mod::Module, sym::Symbol) = safe_clean!(modulechild(mod, sym))

""" TODO """
#Base.compilecache()
function safe_clean!(mod::Module; clean_child_modules=false, bare=false)
  mods = Module[]
  mods_done = Module[]

  print_backtrace(() -> begin
    push!(mods,mod)
    while length(mods)>0
      mod=pop!(mods)
      parent = parentmodule(mod)
      symbol = Symbol(mod)
      name = nameof(mod) #split(string(mod),".")[end] #nameof(mod)
      for sym in names(mod, all=true)
        try
        try value = getfield(mod,sym); catch z; value = sym end
        typ = typeof(value)
        if isa(sym, Symbol) && isdefined(mod, sym)
          if occursin(r"^\#",string(sym)); value = Function; end
          try
            if value == Symbol || isa(value, Symbol) #ignore
            elseif value == Module || isa(value, Module)
              if clean_child_modules && !in(mod,mods_done) && Symbol("Main.",sym) != symbol && sym != symbol && sym != name
                push!(mods,getModule(mod, sym)) #TODO: DEEPWALK
              end
            elseif value == Function || isa(value, Function) Core.eval(mod, :($sym() = nothing))
            elseif isconst(mod, sym)
              if isa(value, Module); #ignore
              elseif isa(value, Number) Core.eval(mod,:(const $sym = $typ(0))) #set to default
              #elseif isprimitivetype(typ)
              #elseif isstructtype(typ) Core.eval(mod,:(const $sym = $typ()))
              #elseif typ == DataType  debug("Custom: $mod.$sym::$typ")
              else Core.eval(mod,:(const $sym = $typ())) #set to default constructor -> when it fails it fails
              end
            else Distributed.clear!(sym;mod=mod)
            #warn("$mod.$sym::$typ is not supported!")
            end
          catch x
            warn("$mod.$sym::$typ -> error!")
          end
        else warn("$mod.$sym::$typ is not valid!")
        end
        catch xxx
          debug("Symbol $mod.$sym error!")
        end
      end
      push!(mods_done, mod)
      #println(names(mod,all=true))
      Core.eval(parent, bare ? :(baremodule $name end) : :(module $name end)) #clear module
      #println(names(eval(Meta.parse(string(Symbol(mod)))),all=true))
    end
  end)

  force_gc()
end
export safe_clean!

end #CodeManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using Documenter
using MEstimation

makedocs(
    sitename = "MEstimation",
    authors = "Ioannis Kosmidis, Nicola Lunardon",
    format = Documenter.HTML(),
    modules = [MEstimation],
    pages = [
        "Home" => "index.md",
        "Examples" => "man/examples.md",
        "Documentation" => Any[
            "Public" => "lib/public.md",
            "Internal" => "lib/internal.md",
        ]
    ],
    doctest = false
)

deploydocs(
    repo = "github.com/ikosmidis/MEstimation.jl.git",
    branch = "gh-pages",
    target = "build",
    devbranch = "develop",
    devurl = "dev",
    push_preview = true,
)


text:
# tags used in  part

# --------------------- <meta />s
"""
    tag_meta( charset::Union{String,Nothing} = "utf-8", lang::Union{String,Nothing} = "en-us", name::Union{String,Nothing} = "", scheme::Union{String,Nothing} = "", content::Union{String,Nothing} = "", http_equiv::Union{String,Nothing} = "")
"""
function tag_meta( ; charset::Union{String,Nothing} = "utf-8", lang::Union{String,Nothing} = "en-us",
    name::Union{String,Nothing} = nothing, scheme::Union{String,Nothing} = nothing, content::Union{String,Nothing} = nothing, http_equiv::Union{String,Nothing} = nothing)
    local tmppairedattrs = Dict{String,String}()
    isa(charset, Nothing) ? nothing : tmppairedattrs["charset"] = charset
    isa(lang, Nothing) ? nothing : tmppairedattrs["lang"] = lang
    isa(name, Nothing) ? nothing : tmppairedattrs["name"] = name
    isa(scheme, Nothing) ? nothing : tmppairedattrs["scheme"] = scheme
    isa(content, Nothing) ? nothing : tmppairedattrs["content"] = content
    isa(http_equiv, Nothing) ? nothing : tmppairedattrs["http-equiv"] = http_equiv
    return SingleHtmlTag( "meta",nothing,nothing,nothing,singleAttributes = String[], pairedAttributes = tmppairedattrs )
end # tag_meta

# ----------------------- linking external sources (dropped those attributes which are not supported by HTML5)
"""
    tag_link( rel::String, type::String, href::String ; hreflang::Union{String, Nothing} = nothing, media::Union{String, Nothing} = nothing )

creates a <link /> tag which is used in  part of the html document.
all attributed not supported by HTML5 are removed.
returns a SingleHtmlTag.

1. :rel is limited to be one of [ \"alternate\", \"author\", \"help\", \"icon\", \"licence\", \"next\", \"pingback\", \"prefetch\", \"prev\", \"search\", \"sidebar\", \"stylesheet\", \"tag\" ]
"""
function tag_link( rel::String, type::String, href::String ;
    hreflang::Union{String, Nothing} = nothing, media::Union{String, Nothing} = nothing )
    # prepare attributes
    local tmppairedattrs = Dict{String,String}(
        "rel" => add_doublequotations(rel),
        "type" => add_doublequotations(type),
        "href" => add_doublequotations(href),
    )
    isa(hreflang, Nothing) ? nothing : tmppairedattrs["hreflang"] = add_doublequotations(hreflang)
    isa(media   , Nothing) ? nothing : tmppairedattrs["media"] = add_doublequotations(media)
    return SingleHtmlTag( "link", nothing, nothing, nothing, singleAttributes = String[], pairedAttributes = tmppairedattrs )
end # tag_link

#

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
# Modelo de gerador eletrico
type Gerador_Eletrico
	Gerador_Eletrico()=begin
		new(
			Eficiencia(Dict{Symbol,Any}(
				:Brief=>"Eficiencia do gerador eletrico"
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Brief=>"Potencia do gerador eletrico",
				:PosX=>1,
				:PosY=>0.5
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia total das turbinas"
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia total da turbina",
				:PosX=>0,
				:PosY=>0.5
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia total da turbina",
				:PosX=>0.02,
				:PosY=>0.25
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia total da turbina",
				:PosX=>0.02,
				:PosY=>0.75
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia total da turbina",
				:PosX=>0.5,
				:PosY=>1
			)),
			[
				:(POT_TURB = POT_TURB1 + POT_TURB2 + POT_TURB3 + POT_TURB4),
				:(POT_GE = EF_GE * POT_TURB),
			],
			[
				"Potencia total das turbinas","Potencia do Gerador Eletrico",
			],
			[:EF_GE,],
			[:POT_GE,:POT_TURB,:POT_TURB1,:POT_TURB2,:POT_TURB3,:POT_TURB4,]
		)
	end
	EF_GE::Eficiencia
	POT_GE::Pot_sinal
	POT_TURB::Potencia
	POT_TURB1::Potencia
	POT_TURB2::Potencia
	POT_TURB3::Potencia
	POT_TURB4::Potencia
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Gerador_Eletrico
function setEquationFlow(in::Gerador_Eletrico)
	addEquation(1)
	addEquation(2)
end
function atributes(in::Gerador_Eletrico,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/gerador"
	drive!(fields,_)
	return fields
end
Gerador_Eletrico(_::Dict{Symbol,Any})=begin
	newModel=Gerador_Eletrico()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Test, MathOptInterface
using LinQuadOptInterface

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const LQOI = LinQuadOptInterface

# Test that using LinQuadOptInterface passes without error. We have to do this
# here because types can't be created in testsets. If it can be defined without
# error, we're okay. This is most likely to error if fields are added to
# @LinQuadOptimizerBase without proper module prefixes.
mutable struct OptimizerTest <: LinQuadOptInterface.LinQuadOptimizer
    LinQuadOptInterface.@LinQuadOptimizerBase
    OptimizerTest() = new()
end

@testset "LinQuadOptInterface" begin
    solver = LQOI.MockLinQuadOptimizer()

    @testset "Printing" begin
        @test sprint(show, solver) == "A LinQuadOptInterface model with " *
            "backend:\nMockLinQuadModel\n    Sense       minimize\n    " *
            "Variables   0\n    Constraints 0\n"
    end

    config = MOIT.TestConfig(solve=false)

    @testset "Unit Tests" begin
        config = MOIT.TestConfig(solve=false)
        MOIT.basic_constraint_tests(solver, config)
        MOIT.unittest(solver, config, [
            "solve_affine_interval",
            "solve_qp_edge_cases",
            "solve_qcp_edge_cases",
            "solve_affine_deletion_edge_cases"
        ])
    end

    @testset "Linear tests" begin
        config = MOIT.TestConfig(solve=false)
        # Exclude partial_start since VariablePrimalStart not implemented.
        MOIT.contlineartest(solver, config, ["partial_start"])
        config = MOIT.TestConfig()
        include("contlinear.jl")
        set_linear1test_solutions!(solver)
        MOIT.linear1test(solver, config)
        set_linear2test_solutions!(solver)
        MOIT.linear2test(solver, config)
        set_linear3test_solutions!(solver)
        MOIT.linear3test(solver, config)
        set_linear4test_solutions!(solver)
        MOIT.linear4test(solver, config)
        set_linear5test_solutions!(solver)
        MOIT.linear5test(solver, config)
        set_linear6test_solutions!(solver)
        MOIT.linear6test(solver, config)
        set_linear7test_solutions!(solver)
        MOIT.linear7test(solver, config)
        set_linear8atest_solutions!(solver)
        MOIT.linear8atest(solver, config)
        set_linear8btest_solutions!(solver)
        MOIT.linear8btest(solver, config)
        set_linear8ctest_solutions!(solver)
        MOIT.linear8ctest(solver, config)
        set_linear9test_solutions!(solver)
        MOIT.linear9test(solver, config)
        set_linear10test_solutions!(solver)
        MOIT.linear10test(solver, config)
        set_linear11test_solutions!(solver)
        MOIT.linear11test(solver, config)
        set_linear12test_solutions!(solver)
        MOIT.linear12test(solver, MOIT.TestConfig(atol=1e-3,rtol=1e-3))
        set_linear13test_solutions!(solver)
        MOIT.linear13test(solver, config)
        set_linear14test_solutions!(solver)
        MOIT.linear14test(solver, config)
        # set_linear15test_solutions!(solver)
        # MOIT.linear15test(solver, config)
    end

    @testset "Quadratic tests" begin
        config = MOIT.TestConfig(solve=false)
        MOIT.contquadratictest(solver, config)
        include("contquadratic.jl")
        config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals=false)
        set_qp1test_solutions!(solver)
        MOIT.qp1test(solver, config)
        set_qp2test_solutions!(solver)
        MOIT.qp2test(solver, config)
        set_qp3test_solutions!(solver)
        MOIT.qp3test(solver, config)
        set_qcp1test_solutions!(solver)
        MOIT.qcp1test(solver, config)
        set_qcp2test_solutions!(solver)
        MOIT.qcp2test(solver, config)
        set_qcp3test_solutions!(solver)
        MOIT.qcp3test(solver, config)
        set_socp1test_solutions!(solver)
        MOIT.socp1test(solver, config)
    end

    @testset "Linear Conic tests" begin
        config = MOIT.TestConfig(solve=false)
        MOIT.lintest(solver, config)
        include("contconic.jl")
        config = MOIT.TestConfig()
        set_lin1test_solutions!(solver)
        MOIT.lin1vtest(solver, config)
        set_lin1test_solutions!(solver)
        MOIT.lin1ftest(solver, config)
        set_lin2test_solutions!(solver)
        MOIT.lin2vtest(solver, config)
        set_lin2test_solutions!(solver)
        MOIT.lin2ftest(solver, config)
        set_lin3test_solutions!(solver)
        MOIT.lin3test(solver, config)
        set_lin3test_solutions!(solver)
        MOIT.lin4test(solver, config)
    end

    @testset "Integer Linear tests" begin
        config = MOIT.TestConfig(solve=false)
        MOIT.intlineartest(solver, config, ["int2"])
        include("intlinear.jl")
        config = MOIT.TestConfig()
        set_knapsacktest_solutions!(solver)
        MOIT.knapsacktest(solver, config)
        set_int1test_solutions!(solver)
        MOIT.int1test(solver, config)
        set_int3test_solutions!(solver)
        MOIT.int3test(solver, config)
    end

    @testset "ModelLike tests" begin
        config = MOIT.TestConfig(solve=false)
        @testset "default_objective_test" begin
            MOIT.default_objective_test(LQOI.MockLinQuadOptimizer())
        end
        @testset "default_status_test" begin
            MOIT.default_status_test(LQOI.MockLinQuadOptimizer())
        end
        @testset "nametest" begin
            MOIT.nametest(LQOI.MockLinQuadOptimizer())
        end
        @testset "validtest" begin
            MOIT.validtest(LQOI.MockLinQuadOptimizer())
        end
        @testset "emptytest" begin
            MOIT.emptytest(LQOI.MockLinQuadOptimizer())
        end
        @testset "orderedindicestest" begin
            MOIT.orderedindicestest(LQOI.MockLinQuadOptimizer())
        end
        @testset "copytest" begin
            MOIT.copytest(LQOI.MockLinQuadOptimizer(),
                          LQOI.MockLinQuadOptimizer())
            MOIT.failcopytestc(LQOI.MockLinQuadOptimizer())
            MOIT.failcopytestia(LQOI.MockLinQuadOptimizer())
            MOIT.failcopytestva(LQOI.MockLinQuadOptimizer())
            MOIT.failcopytestca(LQOI.MockLinQuadOptimizer())
        end
    end
end

MOIU.@model(ModelComplete,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval,
     MOI.Semicontinuous, MOI.Semiinteger),
    (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
     MOI.SecondOrderCone, MOI.RotatedSecondOrderCone, MOI.GeometricMeanCone,
     MOI.ExponentialCone, MOI.DualExponentialCone,
     MOI.PositiveSemidefiniteConeTriangle, MOI.PositiveSemidefiniteConeSquare,
     MOI.RootDetConeTriangle, MOI.RootDetConeSquare, MOI.LogDetConeTriangle,
     MOI.LogDetConeSquare),
    (MOI.PowerCone, MOI.DualPowerCone, MOI.SOS1, MOI.SOS2),
    (MOI.SingleVariable,),
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))
@testset "Copy from/to @Model" begin
    MOIT.copytest(LQOI.MockLinQuadOptimizer(), ModelComplete{Float64}())
    MOIT.copytest(ModelComplete{Float64}(), LQOI.MockLinQuadOptimizer())
end

@testset "Issue #52" begin
    @testset "scalaraffine" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        f1 =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 1.5)
        s1 = MOI.LessThan(2.0)
        c = MOI.add_constraint(model, f1, s1)
        f2 = MOI.get(model, MOI.ConstraintFunction(), c)
        @test f1 ≈ f2
        s2 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s1) == typeof(s2)
        @test s1.upper == s2.upper
    end
    @testset "vectoraffine" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        f1 =  MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1], MOI.ScalarAffineTerm.([1.0], [x])), [1.5])
        s1 = MOI.Zeros(1)
        c = MOI.add_constraint(model, f1, s1)
        f2 = MOI.get(model, MOI.ConstraintFunction(), c)
        @test f1 ≈ f2
        s2 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s1) == typeof(s2)
        @test s1.dimension == s2.dimension
    end
end

@testset "Issue #54" begin
    @testset "scalar, one-sided" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        f =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 1.5)
        s = MOI.LessThan(2.0)
        c = MOI.add_constraint(model, f, s)
        # Change the constraint set, and verify that we get the same set
        # when we retrieve it:
        s2 = MOI.LessThan(1.0)
        MOI.set(model, MOI.ConstraintSet(), c, s2)
        s3 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s2) == typeof(s3)
        @test s2.upper == s3.upper
    end
    @testset "scalar, interval" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        f =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 1.5)
        s = MOI.Interval(2.0, 3.0)
        c = MOI.add_constraint(model, f, s)
        # Change the constraint set, and verify that we get the same set
        # when we retrieve it:
        s2 = MOI.Interval(1.0, 2.0)
        MOI.set(model, MOI.ConstraintSet(), c, s2)
        s3 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s2) == typeof(s3)
        @test s2.lower == s3.lower
        @test s2.upper == s3.upper
    end
end

@testset "Issue #32" begin
    check_set(actual::MOI.GreaterThan, expected::MOI.GreaterThan) =
        @test actual.lower == expected.lower
    check_set(actual::MOI.LessThan, expected::MOI.LessThan) =
        @test actual.upper == expected.upper
    function check_set(actual::MOI.Interval, expected::MOI.Interval)
        @test actual.lower == expected.lower
        @test actual.lower == expected.lower
    end
    function check_set(actual::S, expected::S) where {S <: Union{MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros}}
        @test actual.dimension == expected.dimension
    end

    function check_row(model, constraint_index, expected_function, expected_set)
        actual_function = MOI.get(model, MOI.ConstraintFunction(), constraint_index)
        @test actual_function ≈ expected_function
        actual_set = MOI.get(model, MOI.ConstraintSet(), constraint_index)
        check_set(actual_set, expected_set)
    end

    @testset "Scalar affine function" begin
        @testset "Single variable" begin
            model = LQOI.MockLinQuadOptimizer()
            x = MOI.add_variable(model)

            # Add the constraint and verify that we can retrieve it
            f =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0], [x]), 0.1)
            s = MOI.GreaterThan(0.5)
            c = MOI.add_constraint(model, f, s)
            check_row(model, c, f, s)

            # Replace the constraint function with itself and verify that
            # the problem is unchanged
            MOI.set(model, MOI.ConstraintFunction(), c, f)
            check_row(model, c, f, s)

            # Replace the constraint function with a new function, verify
            # that the replacement occurred and that the set is unchanged
            f2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.5)
            MOI.set(model, MOI.ConstraintFunction(), c, f2)
            check_row(model, c, f2, s)

            # Replace the constraint function with a new function which is not
            # in canonical form:
            f3 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.5], [x]), 0.5)
            MOI.set(model, MOI.ConstraintFunction(), c, f3)
            check_row(model, c, MOIU.canonical(f3), s)

            # Replace the constraint function with a new function whose sparsity
            # pattern is different:
            f4 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 2.0)
            MOI.set(model, MOI.ConstraintFunction(), c, f4)
            check_row(model, c, f4, s)
        end

        @testset "Multiple variables" begin
            model = LQOI.MockLinQuadOptimizer()
            x = MOI.add_variable(model)
            y = MOI.add_variable(model)

            # Add the constraint and verify that we can retrieve it
            f =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0], [x, y]), 0.1)
            s = MOI.GreaterThan(0.5)
            c = MOI.add_constraint(model, f, s)
            check_row(model, c, f, s)

            # Replace the constraint function with itself and verify that
            # the problem is unchanged
            MOI.set(model, MOI.ConstraintFunction(), c, f)
            check_row(model, c, f, s)

            # Replace the constraint function with a new function, verify
            # that the replacement occurred and that the set is unchanged
            f2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.5, 1.0], [x, y]), 0.5)
            MOI.set(model, MOI.ConstraintFunction(), c, f2)
            check_row(model, c, f2, s)

            # Replace the constraint function with a new function which is not
            # in canonical form:
            f3 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-0.1, 0.2], [y, x]), 0.5)
            MOI.set(model, MOI.ConstraintFunction(), c, f3)
            check_row(model, c, MOIU.canonical(f3), s)

            # Replace the constraint function with a new function whose sparsity
            # pattern is different:
            f4 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0], [y]), 2.0)
            MOI.set(model, MOI.ConstraintFunction(), c, f4)
            check_row(model, c, f4, s)
        end
    end

    @testset "Vector affine function" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        y = MOI.add_variable(model)

        # Add the constraint and verify that we can retrieve it
        f = MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
             MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, y)),
             MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(-1.0, x))],
            [0.1, 0.2])
        s = MOI.Nonpositives(2)
        c = MOI.add_constraint(model, f, s)
        check_row(model, c, f, s)

        # Replace the constraint function with itself and verify that
        # the problem is unchanged
        MOI.set(model, MOI.ConstraintFunction(), c, f)
        check_row(model, c, f, s)

        # Replace the constraint function with a new function, verify
        # that the replacement occurred and that the set is unchanged
        f2 = MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(3.0, x)),
             MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-2.0, y)),
             MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(-1.5, x))],
            [0.5, 2.0])
        MOI.set(model, MOI.ConstraintFunction(), c, f2)
        check_row(model, c, f2, s)

        # Replace the constraint function with a new function which is not
        # in canonical form:
        f3 = MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x)),
             MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(-1.5, x)),
             MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-2.5, y))],
            [0.1, 0.2])
        MOI.set(model, MOI.ConstraintFunction(), c, f3)
        check_row(model, c, MOIU.canonical(f3), s)

        # Replace the constraint function with a new function whose sparsity
        # pattern is different:
        f4 = MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(-1.5, y)),
             MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x))],
            [0.1, 0.2])
        MOI.set(model, MOI.ConstraintFunction(), c, f4)
        check_row(model, c, f4, s)
    end
end

@testset "Issue #67" begin
    model = LQOI.MockLinQuadOptimizer()
    function objective_type_test(f::MOI.AbstractScalarFunction)
        MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
        @test MOI.get(model, MOI.ObjectiveFunctionType()) == typeof(f)
    end
    x = MOI.add_variable(model)
    f = MOI.SingleVariable(x)
    objective_type_test(f)
    objective_type_test(convert(MOI.ScalarAffineFunction{Float64}, f))
    objective_type_test(convert(MOI.ScalarQuadraticFunction{Float64}, f))
end

@testset "Conflicting SingleVariable constraints" begin
    @testset "ZeroOne" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semiinteger(2.0, 3.0))
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semicontinuous(2.0, 3.0))
        c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.5))
        c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(0.5))
        @test MOI.is_valid(model, c1)
        @test MOI.is_valid(model, c2)
    end
    @testset "Integer" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Interval(2.0, 3.0))
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semiinteger(2.0, 3.0))
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semicontinuous(2.0, 3.0))
    end
    @testset "Semiinteger" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semiinteger(2.0, 3.0))
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Interval(2.0, 3.0))
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semicontinuous(2.0, 3.0))
    end
    @testset "Semicontinuous" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semicontinuous(2.0, 3.0))
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Interval(2.0, 3.0))
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
        @test_throws Exception MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Semiinteger(2.0, 3.0))
    end
end

@testset "SemiXXX variables" begin
    @testset "Semiinteger" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        c = MOI.add_constraint(
            model, MOI.SingleVariable(x), MOI.Semiinteger(1.0, 4.5))
        @test model.inner.vartype[1] == Cchar('N')
        @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.Semiinteger(1.0, 4.5)
        @test MOI.get(model, MOI.ConstraintFunction(), c) == MOI.SingleVariable(x)
        MOI.delete(model, c)
        @test model.inner.vartype[1] == Cchar('C')
    end
    @testset "Semicontinuous" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.add_variable(model)
        c = MOI.add_constraint(
            model, MOI.SingleVariable(x), MOI.Semicontinuous(1.0, 4.0))
        @test model.inner.vartype[1] == Cchar('S')
        @test MOI.get(model, MOI.ConstraintSet(), c) == MOI.Semicontinuous(1.0, 4.0)
        @test MOI.get(model, MOI.ConstraintFunction(), c) == MOI.SingleVariable(x)
        MOI.delete(model, c)
        @test model.inner.vartype[1] == Cchar('C')
    end
end

@testset "VariablePrimalStart" begin
    model = LQOI.MockLinQuadOptimizer()
    x = MOI.add_variable(model)
    @test_throws MOI.UnsupportedAttribute MOI.set(model, MOI.VariablePrimalStart(), x, 1.0)
end

@testset "Copy names" begin
    model = LQOI.MockLinQuadOptimizer()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    @test MOI.get(model, MOI.VariableName(), x) == "x"
    c = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.set(model, MOI.ConstraintName(), c, "c")
    @test MOI.get(model, MOI.ConstraintName(), c) == "c"

    model2 = LQOI.MockLinQuadOptimizer()
    copy_map = MOI.copy_to(model2, model)
    @test MOI.get(model2, MOI.VariableName(), copy_map.varmap[x]) == "x"
    @test MOI.get(model2, MOI.ConstraintName(), copy_map.conmap[c]) == "c"
end