text:
@generated function transform!(context::MD5_CTX)
    ret = quote
        pbuf = buffer_pointer(context)
    end
    ex  = quote
        @inbounds A = context.state[1]
        @inbounds B = context.state[2]
        @inbounds C = context.state[3]
        @inbounds D = context.state[4]
    end
    push!(ret.args, ex)
    for i in 0:63
        if 0 ≤ i ≤ 15
            ex = :(F = (B & C) | ((~B) & D))
            g = i
        elseif 16 ≤ i ≤ 31
            ex = :(F = (D & B) | ((~D) & C))
            g = 5i + 1
        elseif 32 ≤ i ≤ 47
            ex = :(F = B ⊻ C ⊻ D)
            g = 3i + 5
        elseif 48 ≤ i ≤ 63
            ex = :(F = C ⊻ (B | (~D)))
            g = 7i
        end
        push!(ret.args, ex)
        g = (g % 16) + 1
        ex = quote
            temp = D
            D = C
            C = B
            inner = A + F + $(kk[i+1]) + unsafe_load(pbuf, $g)
            rot_inner = lrot($(ss[i+1]), inner, 32)
            B = B + rot_inner
            A = temp
        end
        push!(ret.args, ex)
    end

    ex = quote
        @inbounds context.state[1] += A
        @inbounds context.state[2] += B
        @inbounds context.state[3] += C
        @inbounds context.state[4] += D
    end
    push!(ret.args, ex)
    quote
        $ret
    end
end

function digest!(context::T) where {T<:MD5_CTX}
    pad_remainder!(context)

    bitcount_idx = div(short_blocklen(T), sizeof(context.bytecount)) + 1
    pbuf = Ptr{typeof(context.bytecount)}(pointer(context.buffer))
    unsafe_store!(pbuf, 8context.bytecount, bitcount_idx)

    # Final transform:
    transform!(context)

    # ctx has been mutated
    reinterpret(UInt8, context.state)
end


text:
#
### Point / vector wrappers
#
@enum ReprType ambiant_repr manifold_repr
mutable struct Point{Tambiantpoint, Tmanifoldpoint}
    amb_repr::Tambiantpoint
    man_repr::Tmanifoldpoint
    M::Manifold
    repr::ReprType
end

function copy(p::Point)
    return Point(deepcopy(p.amb_repr), deepcopy(p.man_repr), deepcopy(p.M), p.repr)
end

embed(x) = x
function embed(x::Point)
    return x.repr == ambiant_repr ? x.amb_repr : embed(x.M, x.man_repr)
end

get_repr(x) = x
function get_repr(x::Point)
    return x.repr == ambiant_repr ? x.amb_repr : x.man_repr
end

function convert_point_repr!(ostate, optimizer, pb)
    convert_point_repr!(ostate.x, ostate.selected_update)
    return
end

ext: jl


text:
# see also the macro versions in datasets.jl

using .CSV

export load_boston, load_ames, load_iris
export load_reduced_ames
export load_crabs

datadir = joinpath(@__DIR__, "..", "data")

load_boston() = CSV.read(joinpath(datadir, "Boston.csv"), copycols=true,
                         categorical=true)

function load_reduced_ames()
    df = CSV.read(joinpath(datadir, "reduced_ames.csv"), copycols=true,
                  categorical=true)
    return coerce(df, :OverallQual => OrderedFactor,
                  :GarageCars => Count,
                  :YearBuilt => Continuous,
                  :YearRemodAdd => Continuous)
end

load_ames() = CSV.read(joinpath(datadir, "ames.csv"), copycols=true,
                  categorical=true)

load_iris() = CSV.read(joinpath(datadir, "iris.csv"), pool=true, copycols=true,
                  categorical=true)

load_crabs() = CSV.read(joinpath(datadir, "crabs.csv"), pool=true,
                  copycols=true, categorical=true)


text:
function main
    println("Test")
end

max_stars_repo_name: Templarghost/TravelingWhithamSolver

max_issues_repo_name: Templarghost/TravelingWhithamSolver

max_forks_repo_name: Templarghost/TravelingWhithamSolver


text:
function DBI.columninfo(
    db::SQLiteDatabaseHandle,
    table::String,
    column::String
)
    datatype = Array(Ptr{Uint8}, 1)
    collseq = Array(Ptr{Uint8}, 1)
    notnull = Array(Cint, 1)
    primarykey = Array(Cint, 1)
    autoinc = Array(Cint, 1)

    status = sqlite3_table_column_metadata(
        db.ptr,
        table,
        column,
        datatype,
        collseq,
        notnull,
        primarykey,
        autoinc
    )
    db.status = status
    if status != SQLITE_OK
        error(errstring(db))
    end

    datatype, length = sql2jltype(bytestring(datatype[1]))

    return DBI.DatabaseColumn(
        column,
        datatype,
        length,
        bytestring(collseq[1]),
        bool(1 - notnull[1]),
        bool(primarykey[1]),
        bool(autoinc[1])
    )
end

function Base.connect(
    ::Type{SQLite3},
    path::String,
    accessmode = convert(Cint, SQLITE_OPEN_READWRITE),
    vfs::Ptr{Void} = C_NULL
)
    dbptrptr = Array(Ptr{Void}, 1)
    # TODO: Check whether path should be UTF8
    status = sqlite3_open_v2(utf8(path), dbptrptr, accessmode, vfs)
    db = SQLiteDatabaseHandle(dbptrptr[1], status)
    db.status = status
    if status != SQLITE_OK
        msg = @sprintf("Failed to connect to database at '%s'", path)
        error(msg)
    end
    return db
end

function DBI.disconnect(db::SQLiteDatabaseHandle)
    status = sqlite3_close_v2(db.ptr)
    db.status = status
    if status != SQLITE_OK
        error(errstring(db))
    end
    return
end

function DBI.errcode(db::SQLiteDatabaseHandle)
    return sqlite3_extended_errcode(db.ptr)
end

function DBI.errstring(db::SQLiteDatabaseHandle)
    errmsg1 = bytestring(sqlite3_errstr(db.status))
    errmsg2 = bytestring(sqlite3_errmsg(db.ptr))
    io = IOBuffer()
    @printf(io, "Error code %d: %s\n", db.status, errmsg1)
    @printf(io, "Error message: %s\n", errmsg2)
    return bytestring(io)
end

function DBI.execute(stmt::SQLiteStatementHandle)
    status = sqlite3_step(stmt.ptr)
    stmt.db.status = status
    if !(status in [SQLITE_OK, SQLITE_ROW, SQLITE_DONE])
        error(errstring(db))
    else
        stmt.executed += 1
    end
    return
end

# Bind parameters, then step and reset
function DBI.execute(stmt::SQLiteStatementHandle, parameters::Vector)
    for index in 1:length(parameters)
        parameter = parameters[index]
        # TODO: Error check these bind calls
        if isa(parameter, FloatingPoint)
            sqlite3_bind_double(stmt.ptr, index, parameter)
        elseif isa(parameter, Integer)
            if WORD_SIZE == 64
                sqlite3_bind_int64(stmt.ptr, index, parameter)
            else
                sqlite3_bind_int(stmt.ptr, index, parameter)
            end
        elseif isa(parameter, NAtype)
            sqlite3_bind_null(stmt.ptr, index)
        # TODO: Allow blob binding
        else
            s = utf8(parameter)
            sqlite3_bind_text(stmt.ptr, index, s, length(s), C_NULL)
        end
    end

    status = sqlite3_step(stmt.ptr)
    stmt.db.status = status
    if !(status in [SQLITE_OK, SQLITE_ROW, SQLITE_DONE])
        error(errstring(stmt.db))
    else
        stmt.executed += 1
    end

    # TODO: Check status of reset call?
    sqlite3_reset(stmt.ptr)

    return
end

# Assumes statement has been executed once already
function DBI.fetchrow(stmt::SQLiteStatementHandle)
    # Unline sqlite3_column_count, this returns number of values
    # actually available at this moment
    ncols = sqlite3_data_count(stmt.ptr)
    row = Array(Any, ncols)
    if ncols == 0
        return row
    end
    for i = 1:ncols
        t = sqlite3_column_type(stmt.ptr, i - 1)
        if t == SQLITE_INTEGER
            if WORD_SIZE == 64
                row[i] = sqlite3_column_int64(stmt.ptr, i - 1)
            else
                row[i] = sqlite3_column_int(stmt.ptr, i - 1)
            end
        elseif t == SQLITE_FLOAT
            row[i] = sqlite3_column_double(stmt.ptr, i - 1)
        elseif t == SQLITE3_TEXT
            row[i] = bytestring(sqlite3_column_text(stmt.ptr, i - 1))
        # TODO: Support blobs
        else
            row[i] = NA
        end
    end
    status = sqlite3_step(stmt.ptr)
    stmt.db.status = status
    if !(status in [SQLITE_OK, SQLITE_ROW, SQLITE_DONE])
        error(errstring(stmt.db))
    end
    return row
end

# Assumes statement has been executed once already
function DBI.fetchall(stmt::SQLiteStatementHandle)
    rows = Array(Any, 0)
    allfetched = false
    while !allfetched
        row = fetchrow(stmt)
        if !isempty(row)
            push!(rows, row)
        else
            allfetched = true
        end
    end
    return rows
end

# Assumes statement has been executed once already
function DBI.fetchdf(stmt::SQLiteStatementHandle)
    ncols = sqlite3_column_count(stmt.ptr)
    cols = Array(Any, ncols)
    names = Array(Symbol, ncols)
    ntypedcols = 0
    for i in 1:ncols
        names[i] = symbol(bytestring(sqlite3_column_name(stmt.ptr, i - 1)))
        t = sqlite3_column_type(stmt.ptr, i - 1)
        if t == SQLITE3_TEXT
            cols[i] = DataArray(UTF8String, 0)
            ntypedcols += 1
        elseif t == SQLITE_FLOAT
            cols[i] = DataArray(Float64, 0)
            ntypedcols += 1
        elseif t == SQLITE_INTEGER
            cols[i] = DataArray(Int, 0)
            ntypedcols += 1
        else
            cols[i] = DataArray(Any, 0)
        end
    end

    while stmt.db.status != SQLITE_DONE
        for i = 1:ncols
            t = sqlite3_column_type(stmt.ptr, i - 1)
            if t == SQLITE3_TEXT
                # TODO: Make this more like raw API
                val = bytestring(sqlite3_column_text(stmt.ptr, i - 1))
            elseif t == SQLITE_FLOAT
                val = sqlite3_column_double(stmt.ptr, i - 1)
            elseif t == SQLITE_INTEGER
                if WORD_SIZE == 64
                    val = sqlite3_column_int64(stmt.ptr, i - 1)
                else
                    val = sqlite3_column_int(stmt.ptr, i - 1)
                end
            else
                val = NA
            end
            push!(cols[i], val)
        end
        execute(stmt)
    end

    # If the type of any column was not determined by its first row,
    # we tighten the type here by finding the first non-NA value
    # and using its type
    if ntypedcols != ncols
        nrows = length(cols[1])
        for j in 1:ncols
            if nrows > 0 && isna(cols[j][1])
                da = cols[j]
                i = 2
                while i <= nrows && isna(da[j])
                    i += 1
                end
                if i <= nrows
                    T = typeof(da[i])
                    # Need to perform convert safely
                    # cols[j] = convert(DataVector{T}, da)
                    tmp = DataArray(T, length(da))
                    for ind in 1:length(da)
                        tmp[ind] = da[ind]
                    end
                    cols[j] = tmp
                end
            end
        end
    end

    return DataFrame(cols, DataFrames.Index(names))
end

function DBI.finish(stmt::SQLiteStatementHandle)
    status = sqlite3_finalize(stmt.ptr)
    stmt.db.status = status
    if status != SQLITE_OK
        error(errstring(db))
    end
    return
end

function DBI.lastinsertid(db::SQLiteDatabaseHandle)
    return sqlite3_last_insert_rowid(db.ptr)
end

function DBI.prepare(db::SQLiteDatabaseHandle, sql::String)
    stmtptrptr = Array(Ptr{Void}, 1)
    status = sqlite3_prepare_v2(
        db.ptr,
        utf8(sql),
        stmtptrptr,
        [C_NULL]
    )
    db.status = status
    if status != SQLITE_OK
        error(errstring(db))
    end
    return SQLiteStatementHandle(db, stmtptrptr[1])
end

# TODO: Put in DBI, export or delete
function sql(stmt::SQLiteStatementHandle)
    return bytestring(sqlite3_sql(stmt.ptr))
end

function DBI.tableinfo(db::SQLiteDatabaseHandle, table::String)
    stmt = prepare(db, "PRAGMA table_info($table)")
    execute(stmt)
    columnlist = fetchall(stmt)
    ncols = size(columnlist, 1)
    finish(stmt)
    cols = Array(DBI.DatabaseColumn, ncols)
    # NB: We are discarding information about default values
    #     provided by PRAGMA table_info().
    for i in 1:ncols
        column = columnlist[i][2]
        cols[i] = columninfo(db, table, column)
    end
    return DBI.DatabaseTable(table, cols)
end


text:
export cycles

function _find_cycle(f::DiscreteFunction, s::Int)
    hits = Set{Int}()
    push!(hits, s)
    current = s

    while true
        current = f(current)
        if in(current, hits)  # we found a cycle!
            result = Array{Int,1}()
            push!(result, current)
            while true
                current = f(current)
                if current == result[1]
                    return result
                end
                push!(result, current)
            end
        end
        push!(hits, current)
    end
    error("This can't happen")
end

function _standardize_cycle(c::Array{Int,1})
    (x,idx) = findmin(c)
    front = c[idx:end]
    back  = c[1:idx-1]
    return vcat(front,back)
end

"""
`cycles(f::DiscreteFunction)` returns a list of the cycles in `f`.
"""
function cycles(f::DiscreteFunction)::Array{Array{Int,1},1}
    result = Set{Array{Int,1}}()

    n = length(f)
    for s=1:n
        c = _find_cycle(f,s)
        c = _standardize_cycle(c)
        push!(result,c)
    end

    return sort(collect(result))
end


text:
import RNG: gen_seed

"""
```julia
ARS1x{R} <: R123Generator1x{UInt128}
ARS1x([seed, R=7])
```

ARS1x is one kind of ARS Counter-Based RNGs. It generates one `UInt128` number at a time.

`seed` is an `Integer` which will be automatically converted to `UInt128`.

`R` denotes to the Rounds which should be at least 1 and no more than 10. With 7 rounds (by default), it has
a considerable safety margin over the minimum number of rounds with no known statistical flaws, but still has
excellent performance.

Only available when [`R123_USE_AESNI`](@ref).
"""
type ARS1x{R} <: R123Generator1x{UInt128}
    x::UInt128
    key::UInt128
    ctr::UInt128
end

function ARS1x(seed::Integer=gen_seed(UInt128), R::Integer=7)
    @assert 1 <= R <= 10
    r = ARS1x{Int(R)}(0, 0, 0)
    srand(r, seed)
end

function srand(r::ARS1x, seed::Integer=gen_seed(UInt128))
    r.key = seed % UInt128
    r.ctr = 0
    random123_r(r)
    r
end

for R = 1:10
    @eval @inline function ars1xm128i(r, ::Type{Val{$R}}, ctr, key)
        p1 = Ptr{UInt128}(pointer_from_objref(ctr))
        p2 = Ptr{UInt128}(pointer_from_objref(key))
        p = Ptr{UInt128}(pointer_from_objref(r))
        ccall(($("ars1xm128i$R"), librandom123), Void, (
        Ptr{UInt128}, Ptr{UInt128}, Ptr{UInt128}
        ), p1, p2, p)
        unsafe_load(p, 1)
    end
end

@inline function random123_r{R}(r::ARS1x{R})
    ars1xm128i(r, Val{R}, r.ctr, r.key)
    (r.x,)
end

"""
```julia
ARS4x{R} <: R123Generator4x{UInt32}
ARS4x([seed, R=7])
```

ARS4x is one kind of ARS Counter-Based RNGs. It generates four `UInt32` numbers at a time.

`seed` is a `Tuple` of four `Integer`s which will all be automatically converted to `UInt32`.

`R` denotes to the Rounds which must be at least 1 and no more than 10. With 7 rounds (by default), it has a
considerable safety margin over the minimum number of rounds with no known statistical flaws, but still has
excellent performance.

Only available when [`R123_USE_AESNI`](@ref).
"""
type ARS4x{R} <: R123Generator4x{UInt32}
    x1::UInt32
    x2::UInt32
    x3::UInt32
    x4::UInt32
    key::UInt128
    ctr1::UInt128
    p::Int
end

function ARS4x(seed::NTuple{4, Integer}=gen_seed(UInt32, 4), R::Integer=7)
    @assert 1 <= R <= 10
    r = ARS4x{Int(R)}(0, 0, 0, 0, 0, 0, 0)
    srand(r, seed)
end

function srand(r::ARS4x, seed::NTuple{4, Integer}=gen_seed(UInt32, 4))
    r.key = unsafe_load(Ptr{UInt128}(pointer_from_objref(seed)), 1)
    r.ctr1 = 0
    p = 0
    random123_r(r)
    r
end

@inline function random123_r{R}(r::ARS4x{R})
    ars1xm128i(r, Val{R}, r.ctr1, r.key)
    (r.x1, r.x2, r.x3, r.x4)
end


text:
## From JuliaSeis
primitive type IBMFloat32 32 end

ieeeOfPieces(fr::UInt32, exp::Int32, sgn::UInt32) =
    reinterpret(Float32, convert(UInt32,fr >>> 9) | convert(UInt32,exp << 23) | sgn) :: Float32

import Base.convert, Base.Float32

function convert(::Type{Float32}, ibm::IBMFloat32)
    local fr::UInt32 = ntoh(reinterpret(UInt32, ibm))
    local sgn::UInt32 = fr & 0x80000000 # save sign
    fr <<= 1 # shift sign out
    local exp::Int32 = convert(Int32,fr >>> 25) # save exponent
    fr <<= 7 # shift exponent out

    if (fr == convert(UInt32,0))
        zero(Float32)
    else
        # normalize the signficand
        local norm::UInt32 = leading_zeros(fr)
        fr <<= norm
        exp = (exp << 2) - 130 - norm
        # exp <= 0 --> ieee(0,0,sgn)
        # exp >= 255 --> ieee(0,255,sgn)
        # else -> ieee(fr<<1, exp, sgn)
        local clexp::Int32 = exp & convert(Int32,0xFF)
        ieeeOfPieces(clexp == exp ? fr << 1 : convert(UInt32,0), clexp, sgn)
    end
end
Float32(ibm::IBMFloat32) = convert(Float32,ibm)
## From JuliaSeis


text:
using BenchmarkTools
using RowMajorArrays

function sum_idx1(a)
    s = zero(eltype(a))
    @inbounds for i = 1:size(a, 1), j = 1:size(a, 2)
        s += a[i, j]
    end
    s
end

function sum_idx2(a)
    s = zero(eltype(a))
    @inbounds for j = 1:size(a, 2), i = 1:size(a, 1)
        s += a[i, j]
    end
    s
end

a1 = rand(1000, 1000)
a1_r = RowMajorArray(a1)

println("Array row major sum")
@btime sum_idx1($a1)
println("Array column major sum")
@btime sum_idx2($a1)

println("RowMajorArray row major sum")
@btime sum_idx1($a1_r)
println("RowMajorArray column major sum")
@btime sum_idx2($a1_r)


text:
# magma_cgesvd wrapper
function magmaZgesvd(jobu, jobvt, m, n, A, lda, s, U,
	 ldu, VT, ldvt, work, lwork, rwork, info)
    ccall((:magma_zgesvd, libmagma),
               Cint,
               (Cint, Cint, Cint, Cint, PtrOrCuPtr{ComplexF64},
			    Cint, PtrOrCuPtr{Float64}, PtrOrCuPtr{ComplexF64},
				Cint, PtrOrCuPtr{ComplexF64}, Cint, PtrOrCuPtr{ComplexF64},
				Cint, PtrOrCuPtr{Float64}, PtrOrCuPtr{Cint}),
                jobu, jobvt, m, n, A, lda, s, U, ldu, VT, ldvt,
			    work, lwork, rwork, info)
end


text:
push!(LOAD_PATH, "../src")
using RobustOptimization
using Ipopt, JuMP, PyPlot
using ScikitLearn
using LinearAlgebra
@sk_import model_selection: train_test_split
@sk_import metrics: roc_curve
@sk_import metrics: roc_auc_score
@sk_import metrics: accuracy_score
@sk_import metrics: classification_report
path_train = "data/ionosphere_scale"
nbfeatures = 34
data_train = read_data_libsvm(path_train, nbfeatures);
# xtr, xte, ytr, yte = train_test_split(data_train[:,1:end-1], data_train[:,end], train_size = 0.15)
# ab_train = hcat(xtr, ytr)
# aux_train = ab_train[ab_train[:,end] .< 0,:]
# ll_train = vcat(aux_train,data_train[data_train[:,end].>0,:])

rs = 75
xtr, xte, ytr, yte = train_test_split(data_train[:,1:end-1], data_train[:,end], train_size = 0.6, random_state = rs)
df_train = hcat(xtr, ytr)
# low_train = hcat(xtr, ytr)
df_aux= hcat(xte, yte);

# df_train, df_aux = create_data("diabetes_scale", nbfeatures, 0.2, train_test_split)
solver = with_optimizer(Ipopt.Optimizer, print_level=2)
N = size(df_train)[1]
ϵ = 0.05
verbosity = 500
itmax = 100
sample = N

ambiguity = "KLdivergence"
robustModel = RobustModel(N, nbfeatures, ϵ, ambiguity, LogisticRegression())
α = 1/norm(robustModel.descent_direction)
# α = 0.1
projParams = ProjParams(Int(40), 1e-5, sample, para_proj=Sequential(), para_inter=Sequential())
optParams = OptParams(itmax, 1e-5, α, verbosity = verbosity);

xnormal = normal_opt(df_train, solver, robustModel.regressionModel)
fpr, tpr, thresholds = roc_curve(df_aux[:,end], positive_rate(xnormal,df_aux))

e_tab = [0.0001,0.0005,0.0006,0.0007,0.0008,0.0009,
            0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,
            0.01,0.05,0.1,0.5]
# e_tab = [0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.05,0.1,0.5,1.0]

dm_tab = []
time_tab = []
xr1_tab = []
xr2_tab = []
fpr1_tab = []
tpr1_tab = []
fpr2_tab = []
tpr2_tab = []
for ϵ in e_tab
    println(" ")
    println("espilon = ", ϵ)
    println(" ")
    robustModel = RobustModel(N, nbfeatures, ϵ, ambiguity, LogisticRegression())
    x0 = initialize(df_train, robustModel, KLConstraint())
    xalg, yalg, dm, mem, mini, t_iter = run_algo(x0, df_train, robustModel, optParams, projParams)
    xrobust1 = getsolution(xalg, ambiguity, nbfeatures)
    xrobust2 = getsolution(mem, ambiguity, nbfeatures)
    fpr1, tpr1, thresholds = roc_curve(df_aux[:,end], positive_rate(xrobust1,df_aux))
    fpr2, tpr2, thresholds = roc_curve(df_aux[:,end], positive_rate(xrobust2,df_aux))
    push!(time_tab, t_iter)
    push!(dm_tab, dm)
    push!(xr1_tab, xrobust1)
    push!(xr2_tab, xrobust2)
    push!(fpr1_tab, fpr1)
    push!(tpr1_tab, tpr1)
    push!(fpr2_tab, fpr2)
    push!(tpr2_tab, tpr2)
end

figure()
for i in 1:(size(e_tab)[1]-1)
    plot(dm_tab[i][1:end].-dm_tab[i][end], label="ϵ = $(e_tab[i])")
end
legend()

figure()
for i in 1:(size(e_tab)[1]-1)
    plot((dm_tab[i]), label="ϵ = $(e_tab[i])")
    yscale("log")
end
legend()

figure()
for i in 1:(size(e_tab)[1]-1)
    plot(cumsum(time_tab[1]), dm_tab[1], label="ϵ = $(e_tab[i])")
    yscale("log")
end
legend(loc=1)

figure()
plot(fpr,tpr, color="black", label="LR")
for i in 1:(size(e_tab)[1]-1)
    plot(fpr1_tab[i], tpr1_tab[i], label="ϵ = $(e_tab[i])")
end
legend()

figure()
plot(fpr,tpr, color="black")
for i in 1:14
    plot(fpr2_tab[i], tpr2_tab[i])
end

plt_AUC = []
print(roc_auc_score(df_aux[:,end], positive_rate(xnormal,df_aux)))
push!(plt_AUC,roc_auc_score(df_aux[:,end], positive_rate(xnormal,df_aux)))
print("  ")
println(accuracy_score(df_aux[:,end], pred(df_aux, xnormal, LogisticRegression())))
for i in 1:(size(e_tab)[1])
    push!(plt_AUC, roc_auc_score(df_aux[:,end], positive_rate(xr1_tab[i],df_aux)))
    print(roc_auc_score(df_aux[:,end], positive_rate(xr1_tab[i],df_aux)))
    print("  ")
    println(accuracy_score(df_aux[:,end], pred(df_aux, xr1_tab[i], LogisticRegression())))
end
figure()
plot(vcat([0],e_tab),plt_AUC, marker="*")
xscale("log")

plt_AUC = []
print(roc_auc_score(df_aux[:,end], positive_rate(xnormal,df_aux)))
push!(plt_AUC,roc_auc_score(df_aux[:,end], positive_rate(xnormal,df_aux)))
print("  ")
println(accuracy_score(df_aux[:,end], pred(df_aux, xnormal, LogisticRegression())))
for i in 1:(size(e_tab)[1])
    push!(plt_AUC, roc_auc_score(df_aux[:,end], positive_rate(xr1_tab[i],df_aux)))
    print(roc_auc_score(df_aux[:,end], positive_rate(xr2_tab[i],df_aux)))
    print("  ")
    println(accuracy_score(df_aux[:,end], pred(df_aux, xr2_tab[i], LogisticRegression())))
end
figure()
plot(vcat([0],e_tab),plt_AUC, marker="*")
xscale("log")

# @everywhere df_train, df_aux = create_data("diabetes_scale", nbfeatures, 0.7, train_test_split)
solver = with_optimizer(Ipopt.Optimizer, print_level=2)
N = size(df_train)[1]
ϵ = 0.05
verbosity = 1
itmax = 50
sample = 1500

ambiguity = "wasserstein"
robustModel = RobustModel(N, nbfeatures, ϵ, ambiguity, LogisticRegression())
α = 1/norm(robustModel.descent_direction)
# α = 1.0
projParams = ProjParams(Int(200), 1e-5, sample, para_proj=Sequential(), para_inter=Sequential())
optParams = OptParams(itmax, 1e-5, α, verbosity = verbosity);

xnormal = normal_opt(df_train, solver, robustModel.regressionModel)
fpr, tpr, thresholds = roc_curve(df_aux[:,end], positive_rate(xnormal,df_aux))

# e_tab_was = [0.1]
# e_tab_was = [0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
e_tab_was = [0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.5,1.0]

mem_was = []
dm_was = []
time_was = []
mini_was = []
xr1_was = []
xr2_was = []
fpr1_was = []
tpr1_was = []
fpr2_was = []
tpr2_was = []
for ϵ in e_tab_was
# for ϵ in [0.001,0.005,0.01]
    println(" ")
    println("espilon = ", ϵ)
    println(" ")
    robustModel = RobustModel(N, nbfeatures, ϵ, ambiguity, LogisticRegression())
    x0 = init_proj(df_train, robustModel, projParams);
    xalg, yalg, dm, mem, mini, t_was = run_algo(x0, df_train, robustModel, optParams, projParams)
    xrobust1 = getsolution(xalg, ambiguity, nbfeatures)
    xrobust2 = getsolution(mem, ambiguity, nbfeatures)
    fpr1, tpr1, thresholds = roc_curve(df_aux[:,end], positive_rate(xrobust1,df_aux))
    fpr2, tpr2, thresholds = roc_curve(df_aux[:,end], positive_rate(xrobust2,df_aux))
    push!(dm_was, dm)
    push!(time_was, dm)
    push!(mini_was, mini)
    push!(mem_was, mem)
    push!(xr1_was, xrobust1)
    push!(xr2_was, xalg)
    push!(fpr1_was, fpr1)
    push!(tpr1_was, tpr1)
    push!(fpr2_was, fpr2)
    push!(tpr2_was, tpr2)
end

figure()
for i in 1:size(e_tab_was)[1]
    plot(dm_was[i], label="ϵ = $(e_tab[i])")
end
legend()

figure()
for i in 1:(size(e_tab_was)[1]-2)
    plot(dm_was[i][1:end].-minimum(dm_was[i]), label="ϵ = $(e_tab_was[i])")
    yscale("log")
end
# plot(dm_was[11][1:end].-minimum(dm_was[11]), label="ϵ = $(e_tab_was[end-1])")
# plot(dm_was[end][1:end].-minimum(dm_was[end]), label="ϵ = $(e_tab_was[end])")
yscale("log")
xlabel("Iterations")
legend(loc = 1, ncol=2)

figure()
for i in 1:size(e_tab_was)[1]
    plot(cumsum(time_was[i][1:end]), dm_was[i][1:end].-mini_was[i], label="ϵ = $(e_tab[i])")
    yscale("log")
end
legend()

figure()
plot(fpr,tpr, color="black", label="LR")
for i in 1:size(e_tab_was)[1]
    plot(fpr1_was[i], tpr1_was[i], label="ϵ = $(e_tab[i])")
end
legend()

figure()
plot(fpr,tpr, color="black")
for i in 1:size(e_tab_was)[1]
    plot(fpr2_was[i], tpr2_was[i])
end

plt_was = []
print(roc_auc_score(df_aux[:,end], positive_rate(xnormal,df_aux)))
push!(plt_was, roc_auc_score(df_aux[:,end], positive_rate(xnormal,df_aux)))
print("  ")
println(accuracy_score(df_aux[:,end], pred(df_aux, xnormal, LogisticRegression())))
for i in 1:size(e_tab_was)[1]
    print(roc_auc_score(df_aux[:,end], positive_rate(xr1_was[i],df_aux)))
    push!(plt_was, roc_auc_score(df_aux[:,end], positive_rate(xr1_was[i],df_aux)))
    print("  ")
    println(accuracy_score(df_aux[:,end], pred(df_aux, xr1_was[i], LogisticRegression())))
end
figure()
plot(vcat([0],e_tab_was),plt_was, marker="*")
xscale("log")

print(roc_auc_score(df_aux[:,end], positive_rate(xnormal,df_aux)))
print("  ")
println(accuracy_score(df_aux[:,end], pred(df_aux, xnormal, LogisticRegression())))
for i in 1:21
    print(roc_auc_score(df_aux[:,end], positive_rate(xr2_was[i],df_aux)))
    print("  ")
    println(accuracy_score(df_aux[:,end], pred(df_aux, xr2_was[i], LogisticRegression())))
end

println(classification_report(df_aux[:,end], pred(df_aux, xnormal, LogisticRegression())))
println(classification_report(df_aux[:,end], pred(df_aux, xr1_was[3], LogisticRegression())))

# open("io-log-was-good.txt", "w") do f
#     write(f, "xnormal = ")
#     write(f, "$xnormal \n")
# #     write(f, "fpr = ")
# #     write(f, "$fpr \n")
# #     write(f, "tpr = ")
# #     write(f, "$tpr \n")
# #     write(f, "dm_tab = ")
# #     write(f, "$dm_tab \n")
# #     write(f, "time_tab = ")
# #     write(f, "$time_tab \n")
# #     write(f, "xr1_tab = ")
# #     write(f, "$xr1_tab \n")
# #     write(f, "fpr1_tab = ")
# #     write(f, "$fpr1_tab \n")
# #     write(f, "tpr1_tab = ")
# #     write(f, "$tpr1_tab \n")
#     write(f, "dm_was = ")
#     write(f, "$dm_was \n")
#     write(f, "time_was = ")
#     write(f, "$time_was \n")
#     write(f, "xr1_was = ")
#     write(f, "$xr1_was \n")
#     write(f, "fpr1_was = ")
#     write(f, "$fpr1_was \n")
#     write(f, "tpr1_was = ")
#     write(f, "$tpr1_was \n")
# end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
#=
:(+) => "i32.add", # these consume 2 args
:(-) => "i32.sub",
:(*) => "i32.mul",
:(div) => "i32.div_s", # julia 4/2 will return float. but div(4,2) wont
:(%) => "i32.rem_s",
:(%) => "i32.rem_s",
:(rem) => "i32.rem_s",
:(&&) => "i32.and",
:(||) => "i32.or",
:(⊻) => "i32.xor",
:(xor) => "i32.xor",

:(<<) => "i32.shl",
:(>>) => "i32.shr_s",
:(==) => "i32.eq",
# :(===) => "i32.eq", #handle this in a special way.
:(!=) => "i32.ne",
:(<) => "i32.lt_s",
:(>) => "i32.gt_s",
:(<=) => "i32.le_s",
:(>=) => "i32.ge_s",
:(!) => "i32.eqz", # these consume 1 arg
:(not_int) => "i32.eqz",
:(leading_zeros) => "i32.clz",
:(trailing_zeros) => "i32.ctz",
:(count_ones) => "i32.popcnt",
:(float) => "f32.convert_i32_s",
:(Int) => "i32.trunc_f32_s",
# :(^) => ["call \$powi",2],
# ""=>["i32.load",1]
# ""=>["i32.store",2]
# ""=>["i32.const",1]
=#

#int
add(a, b) = a + b
sub(a, b) = a - b
mul(a, b) = a * b
_div(a, b) = div(a, b) #should return int
_div2(a, b) = a / b #should return float
_mod(a, b) = mod(a, b)
_mod2(a, b) = a % b
_shl(a, b) = a << b
_shr(a, b) = a >> b #sign preserving
_eq(a, b) = a == b
_egal(a, b) = a === b #identical. (for immutable: same content on bit level. for mutable: type and adress equal)
_ne(a, b) = a != b
_lt(a, b) = a < b
_gt(a, b) = a > b
_le(a, b) = a <= b
_ge(a, b) = a >= b

#bool
_and(a,b) = a && b
_or(a, b) = a || b
_xor(a, b) = xor(a, b)
_xor2(a, b) = a ⊻ b
_not(a) = !a #intended for bools (but wasm dont have bools so actually equivalent to a == 0)
_not2(a) = Base.not_int(a) #intended for bools in julia (but wasm dont have bools so actually equivalent to a == 0)
#bits
_leading_zeros(a) = leading_zeros(a)
_trailing_zeros(a) = trailing_zeros(a)
_count_ones(a) = count_ones(a)
_float(a) = float(a) #int -> float
_Int(a) = Int(a) #whole float -> int

#=
# :(^) => ["call \$powi",2],
# ""=>["i32.load",1]
# ""=>["i32.store",2]
# ""=>["i32.const",1]
=#

function ops(a, b)

  add(a, b)
  sub(a, b)
  mul(a, b)
  _div(a, b)
  _div2(a, b)
  _shl(a, b)
  _shr(a, b)
  _eq(a, b)
  _egal(a, b)
  _ne(a, b)
  _lt(a, b)
  _gt(a, b)
  _le(a, b)
  _ge(a, b)

  #bits
  _leading_zeros(a)
  _trailing_zeros(a)
  _count_ones(a)
  _float(a)

  _mod(a, b)
  _mod2(a, b)
  #_not2(true)
  #_not2(1)

  t=true
  f=false
  #bool
  _and(t, f)
  _or(t, f)
  _xor(t, f)
  _xor2(t, f)
  _not(t)

  return nothing
end

function exports()
  ops(3, 5) #int ops
  #_mod2(3, 5)
end

exports()


text:
using LinearAlgebra, LinearOperators, SparseArrays, Printf, Test

using MPCC

using JuMP
using NLPModelsJuMP
using NLPModels

include("problems.jl")
include("rosenbrock.jl")

printstyled("MPCCMeta tests... ")
include("test-mpcc-meta.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("MPCCCounters tests... ")
include("test-mpcc-counters.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("ADMPCC tests... ")
include("run_test_admpcc.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("MPCCNLPs tests... ")
include("run_test_mpccnlps.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("NLMPCC tests... ")
include("test-mpcc-nlmpcc.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("Bilevel Programming MPCC tests... ")
include("run_test_bpmpcc.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("Ampl MPCC tests... ")
include("run_test_mpccampl.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("MPCCAtX tests... ")
include("test-mpcc-state.jl")
printstyled("passed ✓ \n", color = :green)
printstyled("MPCCStopping tests... ")
include("test-mpcc-stopping.jl")
printstyled("passed ✓ \n", color = :green)


text:
module TestTuning

# using Revise
using Test
using MLJ
using MLJBase
import Random.seed!
seed!(1234)

@load KNNRegressor

include("foobarmodel.jl")

x1 = rand(100);
x2 = rand(100);
x3 = rand(100);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.2*rand(100);

@testset "2-parameter tune, no nesting" begin

    sel = FeatureSelector()
    stand = UnivariateStandardizer()
    ridge = FooBarRegressor()
    composite = MLJ.SimpleDeterministicCompositeModel(transformer=sel, model=ridge)

    features_ = range(composite, :(transformer.features), values=[[:x1], [:x1, :x2], [:x2, :x3], [:x1, :x2, :x3]])
    lambda_ = range(composite, :(model.lambda), lower=1e-6, upper=1e-1, scale=:log10)

    ranges = [features_, lambda_]

    holdout = Holdout(fraction_train=0.8)
    grid = Grid(resolution=10)

    tuned_model = TunedModel(model=composite, tuning=grid,
                             resampling=holdout, measure=rms,
                             ranges=ranges, full_report=false)

    MLJBase.info_dict(tuned_model)

    tuned = machine(tuned_model, X, y)

    fit!(tuned)
    report(tuned)
    tuned_model.full_report=true
    fit!(tuned)
    report(tuned)

    b = fitted_params(tuned).best_model

    measurements = tuned.report.measurements
    # should be all different:
    @test length(unique(measurements)) == length(measurements)

    @test length(b.transformer.features) == 3
    @test abs(b.model.lambda - 0.027825) < 1e-6

    # get the training error of the tuned_model:
    e = rms(y, predict(tuned, X))

    # check this error has same order of magnitude as best measurement
    # during tuning:
    r = e/tuned.report.best_measurement
    @test r < 10 && r > 0.1

    # test weights:
    tuned_model.weights = rand(length(y))
    fit!(tuned)
    @test tuned.report.measurements[1] != measurements[1]

end

#@testset "multiple resolutions" begin

    forest = EnsembleModel(atom=@load(KNNRegressor))
    pipe = @pipeline MyPipe(sel=FeatureSelector(), forest=forest)

    X = MLJ.table(rand(150, 2), names=[:x1, :x2])
    y = rand(150)

    r1 = range(pipe, :(sel.features), values=[[:x1,], [:x1, :x2]])
    r2 = range(pipe, :(forest.atom.K), lower = 1, upper=50)
    r3 = range(pipe, :(forest.bagging_fraction), lower=0.5, upper=1.0)
    ranges = [r1, r2, r3]

    holdout = Holdout(fraction_train=0.8)
    grid = Grid(resolution=5)

    tuned = TunedModel(model=pipe, tuning=grid,
                       resampling=holdout, measure=rms,
                       ranges=ranges)

    mach = machine(tuned, X, y)

    function num_grid_points(resolution)
        grid.resolution = resolution
        fit!(mach)
        return report(mach).measurements |> length
    end

    @test num_grid_points(6) == 72
    @test num_grid_points([:(forest.atom.K)=>7,
                           :(forest.bagging_fraction)=>4]) == 56
    @test num_grid_points([:(forest.atom.K)=>6, ]) == 60

    @testset "one parameter tune" begin
        ridge = FooBarRegressor()
        r = range(ridge, :lambda, lower=1e-7, upper=1e-3, scale=:log)
        tuned_model = TunedModel(model=ridge, ranges=r)
        tuned = machine(tuned_model, X, y)
        fit!(tuned)
        report(tuned)
    end

#end

@testset "nested parameter tune" begin
    tree_model = KNNRegressor()
    forest_model = EnsembleModel(atom=tree_model)
    r1 = range(forest_model, :(atom.K), lower=3, upper=4)
    r2 = range(forest_model, :bagging_fraction, lower=0.4, upper=1.0);
    self_tuning_forest_model = TunedModel(model=forest_model,
                                          tuning=Grid(resolution=3),
                                          resampling=Holdout(),
                                          ranges=[r1, r2],
                                          measure=rms)
    self_tuning_forest = machine(self_tuning_forest_model, X, y)
    fit!(self_tuning_forest, verbosity=2)
    r = report(self_tuning_forest)
    @test length(unique(r.measurements)) == 6
end

## LEARNING CURVE

@testset "learning curves" begin
    atom = FooBarRegressor()
    ensemble = EnsembleModel(atom=atom, n=500)
    mach = machine(ensemble, X, y)
    r_lambda = range(ensemble, :(atom.lambda), lower=0.0001, upper=0.1, scale=:log10)
    curve = MLJ.learning_curve!(mach; range=r_lambda)
    atom.lambda=0.3
    r_n = range(ensemble, :n, lower=10, upper=100)
    curve2 = MLJ.learning_curve!(mach; range=r_n)
end
end # module
true


text:
using Test
using MultinomialRegression
const MultinomialRegression = MultinomialRegression  # Clashes with MLJLinearModels
const fit = MultinomialRegression.fit  # Clashes with MLJLinearModels

using DataFrames
using Dates
using Logging
using MLJLinearModels
using RDatasets
using Random: randperm
using StableRNGs
using Statistics

function mlj_fit(y, X; opts...)
    reg = MLJLinearModels.MultinomialRegression(; opts...)
    b   = MLJLinearModels.fit(reg, X, y)
    nx  = size(X, 2)
    c   = floor(Int, length(b) / nx)
    B   = reshape(b, nx, c)  # (p + 1) x c
end

function mlj_predict(X, B)
    n = size(X, 1)
    c = size(B, 2)
    result = exp.(X * B)  # n x c
    for i = 1:n
        s = sum(result[i, :])
        for j = 1:c
            result[i, j] /= s
        end
    end
    result
end

include("iris.jl")
include("alligators.jl")
#include("compare_to_MLJLinearModels.jl")

@info "$(now()) Finished"


text:
__precompile__()

module Sundials

import Reexport
Reexport.@reexport using DiffEqBase
import DataStructures
import Logging
import DiffEqBase
import SparseArrays
import LinearAlgebra

import Libdl
import CEnum

## TODO: pending https://github.com/JuliaLang/julia/issues/29420
# this one is suggested in the issue, but it looks like time_t and tm are two different things?
# const Ctime_t = Base.Libc.TmStruct

const Ctm = Base.Libc.TmStruct
const Ctime_t = UInt
const Cclock_t = UInt
export Ctm, Ctime_t, Cclock_t

const warnkeywords = (
    :save_idxs,
    :d_discontinuities,
    :isoutofdomain,
    :unstable_check,
    :calck,
    :internalnorm,
    :gamma,
    :beta1,
    :beta2,
    :qmax,
    :qmin,
    :qoldinit,
)

function __init__()
    global warnlist = Set(warnkeywords)
    global warnida = union(warnlist, Set((:dtmin,)))
end

using Sundials_jll

export solve,
    SundialsODEAlgorithm, SundialsDAEAlgorithm, ARKODE, CVODE_BDF, CVODE_Adams, IDA, KINSOL

# some definitions from the system C headers wrapped into the types_and_consts.jl
const DBL_MAX = prevfloat(Inf)
const DBL_MIN = nextfloat(-Inf)
const DBL_EPSILON = eps(Cdouble)

const FILE = Nothing
const __builtin_va_list = Ptr{Cvoid}

include("API/types_and_consts.jl")
include("types_and_consts_additions.jl")

include("handle.jl")
include("nvector_wrapper.jl")

include("API/nvector.jl")
include("API/libsundials.jl")
include("API/sunmatrix.jl")
include("API/sunlinsol.jl")
include("API/sunnonlinsol.jl")
include("API/cvodes.jl")
#include("API/cvode.jl")
include("API/arkode.jl")
include("API/idas.jl")
#include("API/ida.jl")
include("API/kinsol.jl")
# this must be the last one from the API folder!
include("API/api_overloads.jl")

include("simple.jl")
include("common_interface/function_types.jl")
include("common_interface/verbosity.jl")
include("common_interface/algorithms.jl")
include("common_interface/integrator_types.jl")
include("common_interface/integrator_utils.jl")
include("common_interface/solve.jl")

##################################################################
# Deprecations
##################################################################

end # module


text:
module ThJIT

using ..ThArrays

@static if Sys.islinux()
    using LibTorchCAPI_jll
elseif Sys.isapple()
    const libtorch_capi = :libtorch_capi
end

mutable struct CompilationUnit
    mod::Ptr{Nothing}
    owner::Ptr{Nothing}

    function CompilationUnit(m::Ptr{Nothing}, o::Ptr{Nothing})
        ret = new(m, o)
        finalizer(ret) do cu
            ccall((:cunit_destroy, libtorch_capi),
                  Cvoid, (Ptr{Cvoid},),
                  cu.owner)
        end
    end
end

function compile(code::AbstractString)
    fields = [Ptr{Nothing}(0), Ptr{Nothing}(0)]
    cu = ccall((:cunit_compile, libtorch_capi),
               Ptr{Cvoid}, (Ptr{Cvoid}, Cstring),
               fields, pointer(code))
    CompilationUnit(fields[1], fields[2])
end

function run_method(cu::CompilationUnit,
                    method::AbstractString,
                    args::Vararg{Tensor})
    ptrs = map(x -> x.pointer, collect(args))
    tr = ccall((:cunit_run_method, libtorch_capi),
               Ptr{Cvoid}, (Ptr{Cvoid}, Cstring, Ptr{Cvoid}, Cint),
               cu.mod, pointer(method), ptrs, length(args))
    return ThArrays.tensor_from_ptr(tr)
end

struct Function
    cu::CompilationUnit
    method::AbstractString
end

function (f::Function)(args::Vararg{Tensor})
    run_method(f.cu, f.method, args...)
end

get_method(cu::CompilationUnit, method::AbstractString) = Function(cu, method)
Base.getindex(cu::CompilationUnit, method::AbstractString) =
    Function(cu, method)
Base.getindex(cu::CompilationUnit, method::Symbol) =
    Function(cu, string(method))

function Base.getproperty(cu::CompilationUnit, p::Symbol)
    p in fieldnames(CompilationUnit) && return getfield(cu, p)
    return cu[p]
end

end


text:
include("./mt1d.jl");
using MT1D
using PyPlot

"""
Calculate the apparent resistivity and phase from a model and then
plot a comparison of the results with actual data.

## Arguments:
- `modelFile` -- The path to a file containing the resistivity model.
- `dataFile`  -- The path to a file containing a data matrix formatted
  as [T,ρ,Φ,σρ,σΦ]

## Returns:
- Nothing
"""
function plotMT1DCompare(modelFile::AbstractString, dataFile::AbstractString)
    modelFileHandle = open(modelFile);
    model = readdlm(modelFileHandle, Float64);
    close(modelFileHandle);

    dataFileHandle = open(dataFile);
    data = readdlm(dataFileHandle, Float64);
    close(dataFileHandle);

    Ts     = data[:,1];
    fs     = Ts.^-1;
    ρData  = data[:,2];
    ΦData  = data[:,3];

    ρa, Φa = mt1d(model, fs);
    σρData = data[:,4];
    σΦData = data[:,5];

    rmsΦ = rms(ΦData, Φa, σΦData);
    rmsρ = rms(ρData, ρa, σρData);
    rmsTotal = sqrt(rmsΦ^2 + rmsρ^2);
    titleStr = @sprintf("Total RMS: %.2f", rmsTotal);

    subplot(2,1,1);
    loglog(Ts, ρa);
    errorbar(Ts, ρData, yerr=σρData, fmt="x");
    xlabel("Period (S)");
    ylabel("Apparent Resistivity (Ohm-Metres)");
    title(titleStr);

    subplot(2,1,2);
    semilogx(Ts, Φa);
    errorbar(Ts, ΦData, yerr=σΦData, fmt="x");
    xlabel("Period (S)");
    ylabel("Phase (°)");

    legend(["Model", "Data"]);
end


text:
"""
    solve(prob::SingleTermFODEProblem, n, ChebSpectral())

Using Chebyshev spectral method to solve single term fractional ordinary differential equations.

!!! warning
    Chebyshev spectral method only support for linear fractional differential equations, and the time span of the problem should be ``[-1, 1]``.

### References

```tex
@thesis{
    title = {Spectral Methods for Fractional Differential Equations}
    author = {Chinenye Assumpta Nnakenyi, J. A. C. Weideman, Nicholas Hale}
    year = {2015}
}
```
"""
struct ChebSpectral <: FractionalDiffEqAlgorithm end
#FIXME: Fix the initial conditions

function solve(FODE::SingleTermFODEProblem, n, ::ChebSpectral)
    @unpack f, α, u0, tspan = FODE
    (D, x) = Cheb(n, α)
    result = D\f.(x)
    return FODESolution(x, result)
end

function VandermondeLegendre(n)
    x = [-cos(π*j/n) for j in range(0, n, step=1)]
    A = x*x'
    A[1, :] = ones(n+1)
    A[2, :] = x
    for k = 1:n-1
        for i = 0:n
            A[k+2, i+1] = ((2*k+1)*A[2, i+1]*A[k+1, i+1] - k*A[k, i+1])/(k+1)
        end
    end
    return A'
end

function VandermondeJacobi(n, mu)
    x = [-cos(π*j/n) for j in range(0, n, step=1)]
    A = x*x'
    a = -mu
    b = mu
    A0 = 0.5*(a+b)+1
    B0 = 0.5*(a-b)
    A[1, :] = ones(n+1)
    A[2, :] = A0.*x .+ B0
    for k = 1:n-1
        for i = 0:n
            An =((2*k+a+b+1)*(2*k+a+b+2))/(2*(k+1)*(k+a+b+1))
            Bn = ((a^2-b^2)*(2*k+a+b+1))/(2*(k+1)*(k+a+b+1)*(2*k+a+b))
            Cn = ((k+a)*(k+b)*(2*k+a+b+2))/((k+1)*(k+a+b+1)*(2*k+a+b))
            A[k+2, i+1] = (An*x[i+1]+Bn)*A[k+1, i+1] - Cn*A[k, i+1]
        end
    end
    return A'
end

"""
Chebyshev fractional spectral differentiation matrix
"""
function Cheb(n, mu)
    if mu == 1
        if n == 0
            @assert("Wrong")
        else
            k = LinRange(0, 1, n+1)
            D = k*k'
            D[1, 1]= -(2*n^2 +1)/6
            D[n, n]=(2*n^2 +1)/6
            D[n, 1] = (-1)^n/2
            D[1, n] = -(-1)^n/2
            for j = 0:n
                for l = 0:n
                    if j !== l
                        if k[j+1] == 0 || k[j+1]==n
                            c=2
                        else
                            c = 1
                            D[j+1, l+1] = -c*(-1)^(j+l)/(c*(cos(k[j+1]*π) - cos(k[l+1]*π)))
                        end
                    end
                end
            end

            for i = 1:n-1
                D[i+1, i+1] = -(-1)*np.cos(k[i]*pi)/(2*(1-(cos(k[i+1]*pi))^2))
                D[1, i+1] = -2*(-1)^(i)/(1 - (cos(k[i+1]*pi)))
                D[n, i+1] = 2*(-1)^(i+n)/(1 + (cos(k[i+1]*pi)))
                D[i+1, 1] = 1*(-1)^(i)/(2*(1 - (cos(k[i+1]*pi))))
                D[i+1, n] = -(-1)^(i+n)/(2*(1 + (cos(k[i+1]*pi))))
            end
            x = [-cos(π*j/n) for j in range(0, n, step=1)]
            return D, x
        end
        else
            k = zeros(n)
            A = k*k'
            A1 = k*k'
            B = []
            x = [-cos(π*j/(n-1)) for j in range(0, n-1, step=1)]
            B1 = @. 1/((x+1)^mu)
            for k=1:n
                push!(B, gamma(k+mu)/gamma(k))
            end
            for i=0:n-1
                for j=0:n-1
                    if i == j
                        A[i+1, j+1] = B[i+1]
                        A1[i+1, j+1] = B1[i+1]
                    end
                end
            end
            P = VandermondeLegendre(n-1)
            j1 = inv(VandermondeJacobi(n-1, mu))
            D = P*A*j1*A1
            D = D[2:n, 2:n]
            x = x[2:end]
            return D, x
        end
end
#=
n=30
#tspan = collect(-1:0.01:1)

right(t)=t*exp(-t)
#RHS = rightfun.(tspan)
(D, x) = Cheb(n, 0.5)
sol = D\rightfun.(x)

using Plots
plot(x, sol)
=#
#=
U = @. (1+x)^(6+9/17)
plot!(x, U, ls=:dash)
=#


text:
using DifferentialEquations
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

# u[1]=q1, u[2]=q2, u[3]=p1, u[4]=p2, u[5]=q1 hit, u[6]=p2 hit, u[7]=H


text:
using JaynesCummings
using Base.Test

N = 12
ω_q = 2*pi*6.57E9
ω_r = ω_q
g = 0.5*2*pi*19E6
t_f = 500E-9
t_samples = 300
w_samples = 10

include("comparison.jl")

# Test state generation
# Simple case, |g,1>+1im|g,3>
state1 = gen_initialstate(N,["|g,1>","1im|g,3>"])
@test ishermitian(state1)
@test sum(diag(state1)) ≈ 1
@test state1[2,2] ≈ 0.5
@test state1[4,4] ≈ 0.5
@test state1[2,4] ≈ -0.5im
@test state1[4,2] ≈ 0.5im
# More complex case
state2 = gen_initialstate(N,["|g,0>","exp(π/3im)|g,3>","0.5im|g,4>","-3.25|g,7>"])
@test ishermitian(state2)
@test sum(diag(state2)) ≈ 1
ketmag = (1+1+0.5^2+3.25^2)
@test state2[1,1] ≈ 1/ketmag
@test state2[4,4] ≈ 1/ketmag
@test state2[5,5] ≈ 0.5^2/ketmag
@test state2[8,8] ≈ 3.25^2/ketmag
@test state2[1,4] ≈ exp(π/3im)'/ketmag
@test state2[1,5] ≈ -0.5im/ketmag
@test state2[1,8] ≈ -3.25/ketmag
@test state2[4,5] ≈ (exp(π/3im)*-0.5im)/ketmag
@test state2[5,4] ≈ (exp(π/3im)'*0.5im)/ketmag

# State for next tests
ρ = gen_initialstate(N,["|g,0>","0.707106781+0.907106781im|g,3>","1.6025403-0.5im|g,4>","-1.6|g,6>"])
@test ρ ≈ good_initialstate

H = gen_hamiltonian(N,ω_q,ω_r,g)
@test H ≈ good_hamiltonian

times, U_array = gen_timeevoarray(H,t_f,t_samples)
@test times ≈ (0:t_samples-1)*t_f/t_samples
@test U_array[10] ≈ good_time_evo_array_10

e_prob = calc_qubittimeevo(ρ,U_array)
@test e_prob[15] ≈ good_excited_prob_15

photons = calc_photonnumbers(N,g,times,e_prob)
@test photons ≈ good_photons

ρ_resonator = calc_densitymatrix_resonator(ρ,g,times,U_array)
@test ρ_resonator ≈ ptrace(ρ,[2,N],1)

W = calc_wignerfunction_resonator(ρ,g,times,U_array,w_samples)
@test W ≈ calc_wignerfunction_resonator(ρ_resonator,w_samples)

println("Tests passed.")
#println("Now benchmarking...")
#include("perf/bench.jl")


text:
module ReadWriteLocks

using Base: lock, unlock
if VERSION < v"1.2.0-DEV.28"
    using Base.Threads: AbstractLock
else
    using Base: AbstractLock
end

export ReadWriteLock, read_lock, write_lock, lock!, unlock!

struct ReadLock{T<:AbstractLock}
    rwlock::T
end

struct WriteLock{T<:AbstractLock}
    rwlock::T
end

# Need Julia VERSION > v"1.2.0-DEV.28` to have `ReentrantLock <: AbstractLock`
LockTypes = Union{AbstractLock, ReentrantLock}
mutable struct ReadWriteLock{L<:LockTypes} <: AbstractLock
    readers::Int
    writer::Bool
    lock::L  # reentrant mutex
    condition::Condition
    read_lock::ReadLock
    write_lock::WriteLock

    function ReadWriteLock(
        readers::Int=0,
        writer::Bool=false,
        lock::L=ReentrantLock(),
        condition::Condition=Condition(),
    ) where L <: LockTypes
        rwlock = new{L}(readers, writer, lock, condition)
        rwlock.read_lock = ReadLock(rwlock)
        rwlock.write_lock = WriteLock(rwlock)

        return rwlock
    end
end

read_lock(rwlock::ReadWriteLock) = rwlock.read_lock
write_lock(rwlock::ReadWriteLock) = rwlock.write_lock

function lock!(read_lock::ReadLock)
    rwlock = read_lock.rwlock
    lock(rwlock.lock)

    try
        while rwlock.writer
            wait(rwlock.condition)
        end

        rwlock.readers += 1
    finally
        unlock(rwlock.lock)
    end

    return nothing
end

function unlock!(read_lock::ReadLock)
    rwlock = read_lock.rwlock
    lock(rwlock.lock)

    try
        rwlock.readers -= 1
        if rwlock.readers == 0
            notify(rwlock.condition; all=true)
        end
    finally
        unlock(rwlock.lock)
    end

    return nothing
end

function lock!(write_lock::WriteLock)
    rwlock = write_lock.rwlock
    lock(rwlock.lock)

    try
        while rwlock.readers > 0 || rwlock.writer
            wait(rwlock.condition)
        end

        rwlock.writer = true
    finally
        unlock(rwlock.lock)
    end

    return nothing
end

function unlock!(write_lock::WriteLock)
    rwlock = write_lock.rwlock
    lock(rwlock.lock)

    try
        rwlock.writer = false
        notify(rwlock.condition; all=true)
    finally
        unlock(rwlock.lock)
    end

    return nothing
end

end # module


text:
import Random: rand, seed!
import RandomNumbers: AbstractRNG

# Chacha Random Number Generator
# DOUBLEROUNDS=10 means Chacha20
# Generates 32-bit random numbers

mutable struct Chacha{DOUBLEROUNDS} <: AbstractRNG{UInt32}
    v::Vector{UInt32}
    x::Vector{UInt32}
    subCounter::Int
end

function Chacha8()
    Chacha(
        [
            0x243F6A88,
            0x85A308D3,
            0x13198A2E,
            0x03707344,
            0xA4093822,
            0x299F31D0,
            0x082EFA98,
            0xEC4E6C89,
        ],
        UInt64(0),
        4)
    end
function Chacha(
    key::Vector{UInt32} = [
        0x243F6A88,
        0x85A308D3,
        0x13198A2E,
        0x03707344,
        0xA4093822,
        0x299F31D0,
        0x082EFA98,
        0xEC4E6C89,
    ],
    skipLen::UInt64 = UInt64(0),
    DOUBLEROUNDS::Integer = 10, #Chacha20
)
    r = Chacha{DOUBLEROUNDS}(Vector{UInt32}(undef, 16), Vector{UInt32}(undef, 16), Int8(1))
    v = r.v
    v[1] = 0x61707865
    v[2] = 0x3320646e
    v[3] = 0x79622d32
    v[4] = 0x6b206574
    v[5:12] .= key[1:8]
    v[13] = skipLen
    v[14] = skipLen >> 32
    v[15] = 0
    v[16] = 0
    raw(r)
    r
end

@inline function raw(r::Chacha{DOUBLEROUNDS}) where {DOUBLEROUNDS}
    x = r.x
    @inbounds x .= r.v
    # @inbounds for i = 1:16
    #     x[i] = r.v[i]
    # end
    @inbounds for i = 1:DOUBLEROUNDS
        quarterRound(x, 1, 5, 9, 13)
        quarterRound(x, 2, 6, 10, 14)
        quarterRound(x, 3, 7, 11, 15)
        quarterRound(x, 4, 8, 12, 16)
        quarterRound(x, 1, 6, 11, 16)
        quarterRound(x, 2, 7, 12, 13)
        quarterRound(x, 3, 8, 9, 14)
        quarterRound(x, 4, 5, 10, 15)
    end
    @inbounds @simd for i = 1:16
        x[i] += r.v[i]
    end
    r.v[13] += 1
    if (r.v[13] == 0)
        r.v[14] += 1
    end
end

#@inline ROT32(word::UInt32, count::Int) = bitrotate(word,count) #(((word) << count) | ((word) >> (32 - count)))

@inline function quarterRound(v::Vector{UInt32}, a::Int, b::Int, c::Int, d::Int)
    @inbounds begin
        v[a] += v[b]
        v[d] = bitrotate(v[d] ⊻ v[a], 16)
        v[c] += v[d]
        v[b] = bitrotate(v[b] ⊻ v[c], 12)
        v[a] += v[b]
        v[d] = bitrotate(v[d] ⊻ v[a], 8)
        v[c] += v[d]
        v[b] = bitrotate(v[b] ⊻ v[c], 7)
    end
end

@inline function skip(r::Chacha{DOUBLEROUNDS}, skipLength::UInt64) where {DOUBLEROUNDS}
    actualSkip = (skipLength + r.subCounter - 1) / 16
    this.v[13] += -1 + (actualSkip & 0xffffffff) % UInt32
    this.v[14] += (actualSkip >> 32) % UInt32
    raw(r)
    r.subCounter = (((r.subCounter + (skipLength)) % 16) + 1) % Int8
end

@inline function rand(r::Chacha{R}, ::Type{UInt32}) where {R}
    if (r.subCounter >= 17)
        r.subCounter = 1
        raw(r)
    end
    value = r.x[r.subCounter]
    r.subCounter += 1
    value
end

@inline function rand(r::Chacha{R}, ::Type{UInt64}) where {R}
    if (r.subCounter >= 17)
        r.subCounter = 1
        raw(r)
    end
    value = (UInt64(r.x[r.subCounter]) << 32) | r.x[r.subCounter+1]
    r.subCounter += 2
    value
end

#Right now, use UInt32 version
#@inline rand(r::Chacha, ::Type{Float64}) =
#(Float64(rand(r, UInt64) >> 11) + 0.5) / 9007199254740992.0


text:
using AdFem
using PyPlot
using SparseArrays

function f1func(x,y)
    18.8495559215388*pi^2*sin(pi*x)^2*sin(pi*y)*cos(pi*y) - 6.28318530717959*pi^2*sin(pi*y)*cos(pi*x)^2*cos(pi*y) + pi*sin(pi*y)*cos(pi*x)
end
function f2func(x,y)
    -18.8495559215388*pi^2*sin(pi*x)*sin(pi*y)^2*cos(pi*x) + 6.28318530717959*pi^2*sin(pi*x)*cos(pi*x)*cos(pi*y)^2 + pi*sin(pi*x)*cos(pi*y)
end

m = 50
n = 50
h = 1/n
mmesh = Mesh(m, n, h, degree=2)
ν = 0.5
K = ν*constant(compute_fem_laplace_matrix(mmesh))
B = constant(compute_interaction_matrix(mmesh))
Z = [K -B'
    -B spzero(size(B,1))]

bd = bcnode(mmesh)
bd = [bd; bd .+ mmesh.ndof; 2mmesh.ndof + 1]

F1 = eval_f_on_gauss_pts(f1func, mmesh)
F2 = eval_f_on_gauss_pts(f2func, mmesh)
F = compute_fem_source_term(F1, F2, mmesh)
rhs = [F;zeros(mmesh.nelem)]
Z, rhs = impose_Dirichlet_boundary_conditions(Z, rhs, bd, zeros(length(bd)))
sol = Z\rhs

sess = Session(); init(sess)
S = run(sess, sol)

xy = fem_nodes(mmesh)
x, y = xy[:,1], xy[:,2]
U = @. 2*pi*sin(pi*x)*sin(pi*x)*cos(pi*y)*sin(pi*y)
figure(figsize=(12,5))
subplot(121)
visualize_scalar_on_fem_points(U, mmesh)
title("Reference")
subplot(122)
visualize_scalar_on_fem_points(S[1:mmesh.nnode], mmesh)
title("Computed")
savefig("stokes1.png")

U = @. -2*pi*sin(pi*x)*sin(pi*y)*cos(pi*x)*sin(pi*y)
figure(figsize=(12,5))
subplot(121)
visualize_scalar_on_fem_points(U, mmesh)
title("Reference")
subplot(122)
visualize_scalar_on_fem_points(S[mmesh.ndof+1:mmesh.ndof+mmesh.nnode], mmesh)
title("Computed")
savefig("stokes2.png")

xy = fvm_nodes(mmesh)
x, y = xy[:,1], xy[:,2]
p = @. sin(pi*x)*sin(pi*y)
figure(figsize=(12,5))
subplot(121)
visualize_scalar_on_fvm_points(p, mmesh)
title("Reference")
subplot(122)
visualize_scalar_on_fvm_points(S[2mmesh.ndof+1:end], mmesh)
title("Computed")
savefig("stokes3.png")


text:
module Distances

export cohensH

# TODO Convert to string documentation

## Cohen's H
# Measure of distance between two proportions/probabilities
# Variables
#   - p1, p2 (int) = proportions/probabilities
function cohensH(p1, p2)
  return arcsineTransfo(p1)-arcsineTransfo(p2)
end

function arcsineTransfo(p)
  return 2*asin(sqrt(p))
end

end

max_stars_repo_name: DCLukas/HypothesisModules

max_issues_repo_name: DCLukas/HypothesisModules

max_forks_repo_name: DCLukas/HypothesisModules


text:
# InternalOptions
mutable struct InternalOptions
    wrap_structs::Bool
    ismutable::Bool
end
InternalOptions() = InternalOptions(true, false)

"""
    WrapContext
Store shared information about the wrapping session.
"""
mutable struct WrapContext
    index::Index
    headers::Vector{String}
    output_file::String
    common_file::String
    clang_includes::Vector{String}               # clang include paths
    clang_args::Vector{String}                   # additional {"-Arg", "value"} pairs for clang
    header_wrapped::Function                     # called to determine header inclusion status
                                                 #   (top_header, cursor_header) -> Bool
    header_library::Function                     # called to determine shared library for given header
                                                 #   (header_name) -> library_name::AbstractString
    header_outputfile::Function                  # called to determine output file group for given header
                                                 #   (header_name) -> output_file::AbstractString
    cursor_wrapped::Function                     # called to determine cursor inclusion status
                                                 #   (cursor_name, cursor) -> Bool
    common_buf::OrderedDict{Symbol, ExprUnit}    # output buffer for common items: typedefs, enums, etc.
    output_bufs::DefaultOrderedDict{String, Array{Any}}
    options::InternalOptions
    rewriter::Function
end

### Convenience function to initialize wrapping context with defaults
function init(; headers::Vector{String}                    = String[],
                index::Union{Index,Nothing}                = nothing,
                output_dir::String                         = "",
                output_file::String                        = "",
                common_file::String                        = "",
                clang_args::Vector{String}                 = String[],
                clang_includes::Vector{String}             = String[],
                clang_diagnostics::Bool                    = true,
                header_wrapped                             = (header, cursorname) -> true,
                header_library                             = nothing,
                header_outputfile::Union{Function,Nothing} = nothing,
                cursor_wrapped                             = (cursorname, cursor) -> true,
                options                                    = InternalOptions(),
                rewriter                                   = x -> x)

    # Set up some optional args if they are not explicitly passed.
    index == nothing && (index = Index(clang_diagnostics);)

    if output_file == "" && header_outputfile == nothing
        header_outputfile = x->joinpath(output_dir, strip(splitext(basename(x))[1]) * ".jl")
    end

    common_file == "" && (common_file = output_file;)
    common_file = joinpath(output_dir, common_file)

    if header_library == nothing
        header_library = x->strip(splitext(basename(x))[1])
    elseif isa(header_library, String)
        libname = copy(header_library)
        header_library = x->libname
    end
    if header_outputfile == nothing
        header_outputfile = x->joinpath(output_dir, output_file)
    end

    # Instantiate and return the WrapContext
    global context = WrapContext(index,
                                 headers,
                                 output_file,
                                 common_file,
                                 clang_includes,
                                 clang_args,
                                 header_wrapped,
                                 header_library,
                                 header_outputfile,
                                 cursor_wrapped,
                                 OrderedDict{Symbol,ExprUnit}(),
                                 DefaultOrderedDict{String, Array{Any}}(()->Any[]),
                                 options,
                                 rewriter)
    return context
end

function Base.run(wc::WrapContext, generate_template=true)
    # parse headers
    ctx = DefaultContext(wc.index)
    parse_headers!(ctx, wc.headers, args=wc.clang_args, includes=wc.clang_includes)
    ctx.options["is_function_strictly_typed"] = false
    ctx.options["is_struct_mutable"] = wc.options.ismutable
    # Helper to store file handles
    filehandles = Dict{String,IOStream}()
    getfile(f) = (f in keys(filehandles)) ? filehandles[f] : (filehandles[f] = open(f, "w"))

    for trans_unit in ctx.trans_units
        root_cursor = getcursor(trans_unit)
        push!(ctx.cursor_stack, root_cursor)
        header = spelling(root_cursor)
        @info "wrapping header: $header ..."

        # loop over all of the child cursors and wrap them, if appropriate.
        ctx.children = children(root_cursor)
        for (i, child) in enumerate(ctx.children)
            child_name = name(child)
            child_header = filename(child)
            ctx.children_index = i
            # what should be wrapped:
            #   1. wrap if context.header_wrapped(header, child_header) == True
            #   2. wrap if context.cursor_wrapped(cursor_name, cursor) == True
            #   3. skip compiler defs and cursors already wrapped
            if !wc.header_wrapped(header, child_header) ||
               !wc.cursor_wrapped(child_name, child) ||       # client callbacks
               startswith(child_name, "__")          ||       # skip compiler definitions
               child_name in keys(ctx.common_buffer)          # already wrapped
                continue
            end
            ctx.libname = wc.header_library(child_header)

            try
                wrap!(ctx, child)
            catch err
                push!(ctx.cursor_stack, child)
                @error "error thrown. Last cursor available in context.cursor_stack."
                rethrow(err)
            end
        end

        # apply user-supplied transformation
        ctx.api_buffer = wc.rewriter(ctx.api_buffer)

        # write output
        out_file = wc.header_outputfile(header)
        @info "writing $(out_file)"

        out_stream = getfile(out_file)
        println(out_stream, "# Julia wrapper for header: $(basename(header))")
        println(out_stream, "# Automatically generated using Clang.jl\n")
        print_buffer(out_stream, ctx.api_buffer)
        ctx.api_buffer = Expr[]  # clean api_buffer for the next header
    end

    common_buf = dump_to_buffer(ctx.common_buffer)

    # apply user-supplied transformation
    common_buf = wc.rewriter(common_buf)

    # write "common" definitions: types, typealiases, etc.
    open(wc.common_file, "w") do f
        println(f, "# Automatically generated using Clang.jl\n")
        print_buffer(f, common_buf)
    end

    map(close, values(filehandles))

    if generate_template
        copydeps(dirname(wc.common_file))
        print_template(joinpath(dirname(wc.common_file), "LibTemplate.jl"))
    end
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --color=yes --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

using JuMP
import Gurobi

function alternating()
    m = direct_model(Gurobi.Optimizer())
    @variable(m, x[1:100000] >= 0, Int)
    for e in x
        unset_integer(e)
        #is_integer(e)
        #has_upper_bound(e)
        lower_bound(e)
    end
    return
end

function batch()
    m = direct_model(Gurobi.Optimizer())
    @variable(m, x[1:100000] >= 0, Int)
    unset_integer.(x)
    is_integer.(x)
    #is_integer.(x)
    #has_upper_bound.(x)
    lower_bound.(x)
    return
end

function test()
    println("-------------------- batch --------------------------")
    @time batch()
    println("-------------------- alternating --------------------")
    @time alternating()
end

test()
println("\n-------------------- second run ---------------------\n")
test()

# The output of this code in my notebook (Intel(R) Core(TM) i7-7700HQ
# CPU @ 2.80GHz) before PR https://github.com/JuliaOpt/Gurobi.jl/pull/303
# was:
#
# -------------------- batch --------------------------
# Academic license - for non-commercial use only
#   0.758784 seconds (2.34 M allocations: 99.862 MiB, 16.04% gc time)
# -------------------- alternating --------------------
# Academic license - for non-commercial use only
#  24.357266 seconds (1.20 M allocations: 47.019 MiB, 0.08% gc time)
#
# -------------------- second run ---------------------
#
# -------------------- batch --------------------------
# Academic license - for non-commercial use only
#   0.210978 seconds (1.31 M allocations: 50.068 MiB, 8.00% gc time)
# -------------------- alternating --------------------
# Academic license - for non-commercial use only
#  23.389932 seconds (1.20 M allocations: 47.019 MiB, 0.04% gc time)


text:
export ∥

@doc raw"""
# Function call

`∥(z...)`

Type `\parallel` and hit the `tab` key to autocomplete the parallel symbol ∥

# Description

This calculates the parallel connection of impedances. Optionally, the
calculation can be performed using units from the module Unitful.

**Important note:** Use the ∥ operator in parentheses, since its precedence
is equal to the (low) precedence of relation operators or use the function call
with arguments to avoid miscalculations

# Variables

`z` Vector of impedances

# Examples

```julia
julia> using Unitful,Unitful.DefaultSymbols,ElectricalEngineering
julia> ∥(4,6)
2.4
julia> ∥(4Ω,6Ω)
2.4 Ω
julia> 2Ω+(4Ω∥6Ω) # Parentheses are ALWAYS required since
4.4 Ω
julia> 2Ω+4Ω∥6Ω # Is equal to (2Ω+4Ω)∥6Ω which is usually NOT intended
3.0 Ω
julia> (2.4Ω∥(4Ω∥6Ω))
1.2 Ω
julia> ∥(4Ω,6Ω,2.4Ω)
1.2 Ω

```

Type `\Omega` and hit the `tab` key to autocomplete the parallel symbol Ω
"""
function ∥(z...)
    ypar = 0.0/unit(z[1])
    for k=1:length(z)
        ypar = ypar + 1/z[k]
    end
    return 1/ypar
end


text:
module day20

struct Tile
    id
    data
    edges
end

bitstoint(arr) = parse(Int, join(UInt8.(arr)), base=2)

function parseimage(s)
    label, data = split(s, ":\n")
    data = Array{Bool}(hcat(map(line -> '#' .== collect(line), split(data, "\n"))...)')
    edges = bitstoint.(vcat(map(f -> f.([data[:, end], data[end, :], data[:, 1], data[1, :]]), [identity, reverse])...))
    Tile(parse(Int, label[5:end]), data, edges)
end

images = parseimage.(split(strip(read("day20_input.txt", String)), "\n\n"))
adj = Dict(im.id => Set(i.id for i in images if i.id != im.id && !isempty(im.edges ∩ i.edges)) for im in images)
println("Part 1: corner product = $(prod([i.id for i in images if length(adj[i.id]) == 2]))")

# Part 2
fourconnected(i) = map(x -> i + CartesianIndex(x...), [(0,1), (1,0), (0,-1), (-1,0)])

function fillgrid(grid, ids, adj)
    region = CartesianIndices(grid)
    if isempty(ids) return true, grid end

    # Try filling the next empty grid square
    i = first(i for i in region if grid[i] == 0)
    neighbours = [grid[j] for j in fourconnected(i) if j in region && grid[j] > 0]
    candidates = intersect(ids, [adj[n] for n in neighbours]...)
    if isempty(candidates) return false, grid end
    for c in candidates
        grid[i] = c
        newids = setdiff(ids, Set(c))
        result, grid = fillgrid(grid, newids, adj)
        if result return result, grid end
    end
    grid[i] = 0
    false, grid
end

function findlayout(images, adj)
    len = Int(sqrt(length(images)))

    # Try each tile in the top left and then see if we can make progress laying out the rest
    for im in images
        if length(adj[im.id]) > 2 continue end
        grid = zeros(Int, (len, len))
        grid[1,1] = im.id
        remaining = Set(i.id for i in images if i.id != im.id)
        filled, grid = fillgrid(grid, remaining, adj)
        if filled return grid end
    end
end

flip(im) = reverse(im, dims=1)
transforms = [identity, rotr90, rot180, rotl90]
alltransforms = [transforms; flip .∘ transforms]

function findtransform(im, nedges)
    for transform in alltransforms
        imt = transform(im)
        edges = bitstoint.([imt[:, end], imt[end, :], imt[:, 1], imt[1, :]])
        if all(map((a, b) -> isnothing(a) || b ∈ a, nedges, edges)) return imt[2:end-1, 2:end-1] end
    end
    edges = bitstoint.([im[:, end], im[end, :], im[:, 1], im[1, :]])
    error("Failed to find transform: nedges = $nedges, ours = $edges")
end

function countmonsterwaves(im)
    template = ["                  # ",
                "#    ##    ##    ###",
                " #  #  #  #  #  #   "]
    monster = Array{Bool}(hcat(map(line -> '#' .== collect(line), template)...)')
    mc = count(monster)
    region = CartesianIndices(im)
    n = 0
    for i in region
        finish = i + CartesianIndex(size(monster)...) - CartesianIndex(1, 1)
        if finish ∉ region continue end
        n += count(im[i:finish] .* monster) == mc ? mc : 0
    end
    n
end

function countwaves(images, adj)
    grid = findlayout(images, adj)
    full = zeros(Bool, size(grid) .* (size(images[1].data) .- (2, 2)))
    region = CartesianIndices(grid)
    imdict = Dict(i.id => i for i in images)
    unity = CartesianIndex(1, 1)
    for i in region
        nedges = [j in region ? Set(imdict[grid[j]].edges ∩ imdict[grid[i]].edges) : nothing for j in fourconnected(i)]
        im = findtransform(imdict[grid[i]].data, nedges)
        start = (i - unity) * size(im)[1] + unity
        finish = start + CartesianIndex(size(im)...) - unity
        full[start:finish] = im
    end

    for transform in alltransforms
        m = countmonsterwaves(transform(full))
        if m > 0 return sum(count(full)) - m end
    end
end

println("Part 2: roughness = $(countwaves(images, adj))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
# This is used to store experiment details
abstract type AbstractExperiment end

mutable struct Experiment <: AbstractExperiment
    data::DataFrame
    description::String
    selected_features::Array{Int64,1}
    selected_entries::Array{Int64,1}
end

# Constructor
# WARNING: data is not copied so any transformation of the Experiment's data
# will modify the original `data` DataFrame fed to the constructor.
function Experiment(data; description = "No description provided")
    return Experiment(data, description, 1:ncol(data), 1:nrow(data))
end

function Base.show(io::IO, e::Experiment)
    compact = get(io, :compact, false)

    if compact
        show(io, "Exp.: "*string(length(e.selected_entries))*"/"*string(nrow(e.data))*
             "-"*string(length(e.selected_features))*"/"*string(ncol(e.data)))
    else
        show(io, "Experiment with "*string(length(e.selected_entries))*"/"*string(nrow(e.data))*
             " entries and "*string(length(e.selected_features))*"/"*string(ncol(e.data))*
             " features selected.")
    end
end

abstract type AbstractReduce end
abstract type AbstractFilter <: AbstractReduce end
abstract type AbstractSimpleFilter <: AbstractFilter end
abstract type AbstractCombinationFilter <: AbstractFilter end
abstract type AbstractMissingFilter <: AbstractFilter end

mutable struct Filter <: AbstractSimpleFilter
    value::Any
    feature::Symbol
    compare::Function
    description::String
end

# Constructor
function Filter(value, feature;
                compare = isequal, description = "No description provided")
    return Filter(value, feature, compare, description)
end

mutable struct CombinationFilter <: AbstractCombinationFilter
    filter1::AbstractFilter
    filter2::AbstractFilter
    operator::Function
end

mutable struct MissingFilter <: AbstractMissingFilter end

# Constructor
function MissingFilter(feature;
                        description = "No description provided")
    return Filter(1, feature, (x,y) -> !ismissing(x), description)
end

# Methods

"""Return filtered entries in an Experiment `e` based on filter `f`
"""
function filter_entries(e::AbstractExperiment, f::AbstractSimpleFilter)
    expEntries = e.data[e.selected_entries, f.feature]
    return(e.selected_entries[f.compare.(expEntries, f.value)])
end

function filter_entries(e::AbstractExperiment, f::AbstractCombinationFilter)
    e1 = filter_entries(e, f.filter1)
    e2 = filter_entries(e, f.filter2)
    return(sort(f.operator(e1, e2)))
end

function filter_entries(e::AbstractExperiment, f::AbstractMissingFilter)
    return(sort(∩([filter_entries(e, MissingFilter(f))
                   for f in Symbol.(names(e.data)[e.selected_features])]...)))
end

"""Filter entries in an Experiment `e` based on filter(s) `f`,
updating `e.selected_entries` in place accordingly.
"""
function filter_entries!(e::AbstractExperiment, f::AbstractFilter)
    # Currently returns the indices kept
    e.selected_entries = filter_entries(e,f)
end

function filter_entries!(e::AbstractExperiment, filters::Array{T,1}) where {T<:AbstractFilter}
    for f in filters
        filter_entries!(e, f)
    end
end

abstract type AbstractReduce end
abstract type AbstractSelector <: AbstractReduce end
abstract type AbstractSimpleSelector <: AbstractSelector end
abstract type AbstractNameSelector <: AbstractSelector end
abstract type AbstractCombinationSelector <: AbstractSelector end

mutable struct Selector <: AbstractSimpleSelector
    summarize::Function
    subset::Union{Function, Nothing}
    description::String
end

# Constructor
function Selector(summarize; subset = nothing,
                  description = "No description provided")
    return Selector(summarize, subset, description)
end

mutable struct NameSelector <: AbstractNameSelector
    summarize::Function
    description::String
end

# Constructor
function NameSelector(summarize; description = "No description provided")
    return NameSelector(summarize, description)
end

mutable struct CombinationSelector <: AbstractCombinationSelector
    selector1::AbstractSelector
    selector2::AbstractSelector
    operator::Function
end

# Methods

"""Return selected features in an Experiment `e` based on selectors `s`
"""
function select_features(e::AbstractExperiment, s::AbstractSimpleSelector)
    # NB: isnothing was not implemented in 1.0
    if s.subset === nothing
        data = e.data[e.selected_entries, e.selected_features]
    else
        subIndices = s.subset(e.data[e.selected_entries,:])
        data = e.data[e.selected_entries[subIndices], e.selected_features]
    end
    selectedFtDF = mapcols(s.summarize, data)
    return(e.selected_features[[x for x in selectedFtDF[1,:]]])
end

function select_features(e::AbstractExperiment, s::AbstractNameSelector)
    selectedFtDF = map(s.summarize, names(e.data[:,e.selected_features]))
    return(e.selected_features[selectedFtDF])
end

function select_features(e::AbstractExperiment, s::AbstractCombinationSelector)
    f1 = select_features(e, s.selector1)
    f2 = select_features(e, s.selector2)
    return(sort(s.operator(f1, f2)))
end

"""Return selected features in an Experiment `e` based on selectors `s`,
updating `e.selected_features` in place accordingly.
"""
function select_features!(e::AbstractExperiment, s::AbstractSelector)
    # Currently returns the indices kept
    e.selected_features = select_features(e,s)
end

function select_features!(e::AbstractExperiment, selectors::Array{T,1}) where {T<:AbstractSelector}
    for s in selectors
        select_features!(e, s)
    end
end

# Not exported in favor of the shorter "filter!" and "select!"
"""For an experiment `e`, update in place `e.selected_features` and
`e.selected_entries` based on an array `arr` of feature selectors and
entry filters. Filters and selectors are applied sequentially.
"""
function filter_experiment!(e::AbstractExperiment, arr::Array{T,1}) where {T<:AbstractReduce}
    for a in arr
        filter!(e, a)
    end
end

"""For an experiment `e`, update in place `e.selected_features` and
`e.selected_entries` based on a feature selector `s`.
"""
function filter_experiment!(e::AbstractExperiment, s::AbstractSelector)
    select_features!(e,s)
end

"""For an experiment `e`, update in place `e.selected_features` and
`e.selected_entries` based on an entry filter `f`.
"""
filter_experiment!(e::AbstractExperiment, f::AbstractFilter) = filter_entries!(e,f)

# Aliases
@doc (@doc filter_experiment!)
Base.filter!(e::AbstractExperiment,
                   arr::Array{T,1}) where {T<:AbstractReduce} = filter_experiment!(e,arr)
Base.filter!(e::AbstractExperiment, s::AbstractSelector) = filter_experiment!(e,s)
Base.filter!(e::AbstractExperiment, f::AbstractFilter) = filter_experiment!(e,f)

@doc (@doc filter_experiment!)
DataFrames.select!(e::AbstractExperiment,
                   arr::Array{T,1}) where {T<:AbstractReduce} = filter_experiment!(e,arr)
DataFrames.select!(e::AbstractExperiment, s::AbstractSelector) = filter_experiment!(e,s)
DataFrames.select!(e::AbstractExperiment, f::AbstractFilter) = filter_experiment!(e,f)

"""Return a negative Filter or Selector by inverting
the entries or features that are kept and excluded.
"""
function negation(r::Union{AbstractNameSelector,AbstractSimpleSelector})
    neg_r = deepcopy(r)
    neg_r.summarize = !neg_r.summarize
    neg_r.description = "Do not "*neg_r.description
    return(neg_r)
end

function negation(r::AbstractSimpleFilter)
    neg_r = deepcopy(r)
    neg_r.compare = !neg_r.compare
    neg_r.description = "Do not "*neg_r.description
    return(neg_r)
end

"""Return a copy of the data in Experiment `e` for its
selected entries and features.
"""
function getdata(e::Experiment)
    return(e.data[e.selected_entries, e.selected_features])
end


text:
# Precompilation statements

#=
    Reading and parsing
=#
precompile(parse_node, (EzXML.Node, Bool))

# Hacky way to get at all the types we create
for name in names(StationXML, all=true)
    name_str = String(name)
    if occursin(r"^[A-Z]", name_str) && !occursin(r"^BASE", name_str)
        name in (:M, :COMPARABLE_TYPES, :StationXML, :Numerator, :Denominator) && continue
        getfield(StationXML, name) isa DataType || continue
        @eval begin
            # Parsing
            precompile(parse_node, (Type{$name}, EzXML.Node, Bool))
            precompile(parse_node, (Type{Union{$name,Missing}}, EzXML.Node, Bool))
            # Construction
            precompile($name, ())
            # Equality
            precompile(==, ($name, $name))
            precompile(==, ($name, Missing))
            precompile(==, (Missing, $name))
            precompile(local_equals, ($name, $name))
            precompile(local_equals, ($name, Missing))
            precompile(local_equals, (Missing, $name))
        end
        # Others
        for func in (:attribute_fields, :element_fields, :has_text_field, :text_field,
                :removed_fields)
            @eval precompile($func, ($name,))
        end
    end
end

# Hacky way to get members of a Union
_union_types(x::Union) = (x.a, _union_types(x.b)...)
_union_types(x::Type) = (x,)

for T in _union_types(ParsableTypes)
    @eval begin
        precompile(parse_node, (Type{$T}, EzXML.Node, Bool))
        precompile(parse_node, (Type{Union{Missing,$T}}, EzXML.Node, Bool))
        precompile(local_parse, (Type{$T}, String))
        precompile(local_parse, (Type{Union{$T,Missing}}, String))
        precompile(local_tryparse, (Type{$T}, String))
        precompile(local_tryparse, (Type{Union{$T,Missing}}, String))
    end
end

precompile(StationXML.read, (String,))
precompile(readstring, (String,))

#=
    Merging and appending
=#
for T in (Network, Station, Channel)
    precompile(_time_ranges_overlap, (T, T))
end
precompile(_merge!, (Nothing, FDSNStationXML, FDSNStationXML, Bool))
precompile(_merge!, (FDSNStationXML, Network, Network, Bool))
precompile(_merge!, (Network, Station, Station, Bool))
precompile(Base.append!, (Network, Network))
precompile(Base.append!, (FDSNStationXML, FDSNStationXML))
precompile(Base.merge, (FDSNStationXML, FDSNStationXML))
precompile(Base.merge, (FDSNStationXML, FDSNStationXML, FDSNStationXML))
precompile(Base.merge!, (FDSNStationXML, FDSNStationXML))
precompile(Base.merge!, (FDSNStationXML, FDSNStationXML, FDSNStationXML))


text:
module AssemblyCalculus

    using Distributions: Binomial, Bernoulli
    using Graphs
    using StatsBase
    using ThreadsX

    include("types.jl")
    include("simulation_functions.jl")
    include("partial_neural_area.jl")

    export Brain, Assembly, Stimulus
    export simulate!, random_firing!, rand_stim, zero_stim

end  # module AssemblyCalculus


text:
module CrossValidators

using Statistics: mean, std

using MLBase: StratifiedKfold, LOOCV,
              RandomSub, StratifiedRandomSub

# standard included modules
using Random
using AutoMLPipeline
using AutoMLPipeline.Utils
using DataFrames: DataFrame

export crossvalidate, holdoutcrossvalidate

function holdoutcrossvalidate(pl::Machine,X::DataFrame,Y::Vector,metric::Function,ntimes=10)
  @assert size(X)[1] == length(Y)
end

"""
    crossvalidate(pl::Machine,X::DataFrame,Y::Vector,pfunc::Function,kfolds=10)

Run K-fold crossvalidation where:
- `pfunc` is a performance metric
- `X` and `Y` are input and target
"""
function crossvalidate(pl::Machine,X::DataFrame,Y::Vector,
		       pfunc::Function,nfolds::Int=10,verbose::Bool=true)
  ## flatten arrays
  @assert size(X)[1] == length(Y)
  rowsize = size(X)[1]
  folds = StratifiedKfold(Y,nfolds) |> collect
  pacc = Float64[]
  fold = 0
  error = 0
  for trndx in folds
    ppl = deepcopy(pl)
    rlist = collect(1:rowsize) # nasty bug if not collected and placed inside setdiff
    tsndx = setdiff(rlist,trndx)
    trX = X[trndx,:]
    trY = Y[trndx] |> collect
    tstX = X[tsndx,:]
    tstY = Y[tsndx] |> collect
    res = 0.0
    try
      res = pipe_accuracy(ppl,pfunc,trX,trY,tstX,tstY)
      push!(pacc,res)
      fold += 1
      if verbose == true
	println("fold: ",fold,", ",res)
      end
    catch e
      #println(e)
      error += 1
    end
  end
  if verbose == true
    println("errors: ",error)
  end
  (mean=mean(pacc),std=std(pacc),folds=nfolds,errors=error)
end

function pipe_accuracy(plearner,perf::Function,trX::DataFrame,
		       trY::Vector,tstX::DataFrame,tstY::Vector)
  learner = deepcopy(plearner)
  fit!(learner,trX,trY)
  pred = transform!(learner,tstX)
  res = perf(pred,tstY)
  return res
end

end


text:
function get_true_roadway()
    roadway = Roadway()
    nlanes = 1
    seg1 = VecSE2(-10.0,0.0,0.0)
    length1 = 10.
    add_line!(seg1,nlanes,length1,roadway)

    seg2 = VecSE2(5.0,-5.0,-pi/2)
    length2 = 10.
    add_line!(seg2,nlanes,length2,roadway)

    seg3 = VecSE2(0.0, 3.0, pi)
    length3 = 10.
    add_line!(seg3, nlanes, length3, roadway)

    seg4 = VecSE2(23.0, 3.0, pi)
    length4 = 10.
    add_line!(seg4, nlanes, length4, roadway)

    seg5 = VecSE2(13.0, 0.0, 0.0)
    length5 = 10.
    add_line!(seg5, nlanes, length5, roadway)

    seg6 = VecSE2(8.0, -15, pi/2)
    length6 = 10.
    add_line!(seg6, nlanes, length6, roadway)

    junction = Junction([Connection(1,2), Connection(1,5), Connection(6,3), Connection(4,3), Connection(6,5), Connection(4,2)])  #, Connection(3,4), Connection(1,5), Connection(6,5), Connection(4,2)])

    add_junction!(junction, roadway)

    roadway
end

function true_scene(scene::Scene, roadway::Roadway)
    new_scene = Scene
    for veh in scene
        new_veh_state = VehicleState(veh.state.posG, roadway, veh.state.v)
        new_veh = Vehicle(new_veh_state, veh.def, veh.id)
        push!(new_scene, new_veh)
    end
    return new_scene
end

# doc, r = initialize_XML()
# convert_roadway!(r, roadway)
# write("test_intersection.xodr", doc)

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Revise

using Genie.Router, Genie.Renderer.Html
using Stipple

using StippleUI, StippleUI.Table, StippleUI.Range, StippleUI.BigNumber, StippleUI.Heading, StippleUI.Dashboard
using StippleCharts, StippleCharts.Charts

import StippleUI.Range: range

using CSV, DataFrames

# configuration
const data_opts = DataTableOptions(columns = [Column("Good_Rating"), Column("Amount", align = :right),
                                              Column("Age", align = :right), Column("Duration", align = :right)])

const plot_colors = ["#72C8A9", "#BD5631"]

const bubble_plot_opts = PlotOptions(data_labels_enabled=false, fill_opacity=0.8, xaxis_tick_amount=10, chart_animations_enabled=false,
                                      xaxis_max=80, xaxis_min=17, yaxis_max=20_000, chart_type=:bubble,
                                      colors=plot_colors, plot_options_bubble_min_bubble_radius=4, chart_font_family="Lato, Helvetica, Arial, sans-serif")

const bar_plot_opts = PlotOptions(xaxis_tick_amount=10, xaxis_max=350, chart_type=:bar, plot_options_bar_data_labels_position=:top,
                                  plot_options_bar_horizontal=true, chart_height=200, colors=plot_colors, chart_animations_enabled=false,
                                  xaxis_categories = ["20-30", "30-40", "40-50", "50-60", "60-70", "70-80"], chart_toolbar_show=false,
                                  chart_font_family="Lato, Helvetica, Arial, sans-serif")

# model
data = CSV.File("data/german_credit.csv") |> DataFrame!

Base.@kwdef mutable struct Dashboard <: ReactiveModel
  credit_data::R{DataTable} = DataTable()
  credit_data_pagination::DataTablePagination = DataTablePagination(rows_per_page=100)
  credit_data_loading::R{Bool} = false

  range_data::R{RangeData{Int}} = RangeData(15:80)

  big_numbers_count_good_credits::R{Int} = 0
  big_numbers_count_bad_credits::R{Int} = 0
  big_numbers_amount_good_credits::R{Int} = 0
  big_numbers_amount_bad_credits::R{Int} = 0

  bar_plot_options::PlotOptions = bar_plot_opts
  bar_plot_data::R{Vector{PlotSeries}} = []

  bubble_plot_options::PlotOptions = bubble_plot_opts
  bubble_plot_data::R{Vector{PlotSeries}} = []
end

# functions
function creditdata(data::DataFrame, model::M) where {M<:Stipple.ReactiveModel}
  model.credit_data[] = DataTable(data, data_opts)
end

function bignumbers(data::DataFrame, model::M) where {M<:ReactiveModel}
  model.big_numbers_count_good_credits[] = data[(data[:Good_Rating] .== true), [:Good_Rating]] |> nrow
  model.big_numbers_count_bad_credits[] = data[(data[:Good_Rating] .== false), [:Good_Rating]] |> nrow
  model.big_numbers_amount_good_credits[] = data[(data[:Good_Rating] .== true), [:Amount]] |> Array |> sum
  model.big_numbers_amount_bad_credits[] = data[(data[:Good_Rating] .== false), [:Amount]] |> Array |> sum
end

function barstats(data::DataFrame, model::M) where {M<:Stipple.ReactiveModel}
  age_stats = Dict{Symbol,Vector{Int}}(:good_credit => Int[], :bad_credit => Int[])

  for x in 20:10:70
    push!(age_stats[:good_credit],
          data[(data[:Age] .∈ [x:x+10]) .& (data[:Good_Rating] .== true), [:Good_Rating]] |> nrow)
    push!(age_stats[:bad_credit],
          data[(data[:Age] .∈ [x:x+10]) .& (data[:Good_Rating] .== false), [:Good_Rating]] |> nrow)
  end

  model.bar_plot_data[] = [PlotSeries("Good credit", PlotData(age_stats[:good_credit])),
                            PlotSeries("Bad credit", PlotData(age_stats[:bad_credit]))]
end

function bubblestats(data::DataFrame, model::M) where {M<:ReactiveModel}
  selected_columns = [:Age, :Amount, :Duration]
  credit_stats = Dict{Symbol,DataFrame}()

  credit_stats[:good_credit] = data[data[:Good_Rating] .== true, selected_columns]
  credit_stats[:bad_credit] = data[data[:Good_Rating] .== false, selected_columns]

  model.bubble_plot_data[] = [PlotSeries("Good credit", PlotData(credit_stats[:good_credit])),
                              PlotSeries("Bad credit", PlotData(credit_stats[:bad_credit]))]
end

function setmodel(data::DataFrame, model::M)::M where {M<:ReactiveModel}
  creditdata(data, model)
  bignumbers(data, model)

  barstats(data, model)
  bubblestats(data, model)

  model
end

### UI
Stipple.register_components(Dashboard, StippleCharts.COMPONENTS)

model = setmodel(data, Dashboard()) |> Stipple.init

function filterdata(model::Dashboard)
  model.credit_data_loading[] = true
  model = setmodel(data[(model.range_data[].range.start .<= data[:Age] .<= model.range_data[].range.stop), :], model)
  model.credit_data_loading[] = false

  nothing
end

function ui()
  dashboard(root(model), [
    heading("German Credits by Age")

    row([
      cell(class="st-module", [
        row([
          cell(class="st-br", [
            bignumber("Bad credits",
                      :big_numbers_count_bad_credits,
                      icon="format_list_numbered",
                      color="negative")
          ])

          cell(class="st-br", [
            bignumber("Good credits",
                      :big_numbers_count_good_credits,
                      icon="format_list_numbered",
                      color="positive")
          ])

          cell(class="st-br", [
            bignumber("Bad credits total amount",
                      R"big_numbers_amount_bad_credits | numberformat",
                      icon="euro_symbol",
                      color="negative")
          ])

          cell(class="st-br", [
            bignumber("Good credits total amount",
                      R"big_numbers_amount_good_credits | numberformat",
                      icon="euro_symbol",
                      color="positive")
          ])
        ])
      ])
    ])

    row([
      cell([
        h4("Age interval filter")

        range(18:1:90,
              :range_data;
              label=true,
              labelvalueleft="'Min age: ' + range_data.min",
              labelvalueright="'Max age: ' + range_data.max")
      ])
    ])

    row([
      cell(class="st-module", [
        h4("Credits data")

        table(:credit_data;
              style="height: 400px;",
              pagination=:credit_data_pagination,
              loading=:credit_data_loading
        )
      ])
      cell(class="st-module", [
        h4("Credits by age")
        plot(:bar_plot_data; options=:bar_plot_options)
      ])
    ])

    row([
      cell(class="st-module", [
        h4("Credits by age, amount and duration")
        plot(:bubble_plot_data; options=:bubble_plot_options)
      ])
    ])

    footer(class="st-footer q-pa-md", [
      cell([
        img(class="st-logo", src="/img/st-logo.svg")
        span(" © 2020")
      ])
    ])
  ], title="German Credits by Age") |> html
end

# handlers
on(model.range_data) do _
  filterdata(model)
end

# routes
route("/", ui)

# JS deps
function __init__()
  push!(Stipple.DEPS, () -> script(src="/js/plugins/genie_autoreload/autoreload.js"))
end

# start server
up()


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iupcontrols.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static is_windows()? (const libiupcontrols_ = "iupcontrols") : (const libiupcontrols_ = "libiupcontrols")  # Name of the shared lib.

const IUP_PRIMARY = -1
const IUP_SECONDARY = -2

function IupControlsOpen()
	ccall((:IupControlsOpen, libiupcontrols_), Cint, (),)
end
function IupControlsClose()
	ccall((:IupControlsClose, libiupcontrols_), Void, (),)
end
function IupColorbar()
	ccall((:IupColorbar, libiupcontrols_), Ptr{Cint}, (),)
end
function IupCells()
	ccall((:IupCells, libiupcontrols_), Ptr{Cint}, (),)
end
function IupColorBrowser()
	ccall((:IupColorBrowser, libiupcontrols_), Ptr{Cint}, (),)
end
function IupGauge()
	ccall((:IupGauge, libiupcontrols_), Ptr{Cint}, (),)
end
function IupDial(orient::String="")
	ccall((:IupDial, libiupcontrols_), Ptr{Ihandle}, (Ptr{UInt8},), orient)
end
function IupMatrix(action::String="")
	ccall((:IupMatrix, libiupcontrols_), Ptr{Ihandle}, (Ptr{UInt8},), action)
end
function IupMatrix(action::Ptr{Void})
	ccall((:IupMatrix, libiupcontrols_), Ptr{Ihandle}, (Ptr{UInt8},), action)
end
function IupMatrixList()
	ccall((:IupMatrixList, libiupcontrols_), Ptr{Ihandle}, (),)
end
# DEPRECATED
function IupMatSetAttribute(ih::Ptr{Int}, name::String, lin::Int, col::Int, value::String)
	ccall((:IupMatSetAttribute, libiupcontrols_), Void, (Ptr{Cint}, Ptr{UInt8}, Cint, Cint, Ptr{UInt8}), ih, name, lin, col, value)
end
function IupMatStoreAttribute(ih::Ptr{Int}, name::String, lin::Int, col::Int, value::String)
	ccall((:IupMatStoreAttribute, libiupcontrols_), Void, (Ptr{Cint}, Ptr{UInt8}, Cint, Cint, Ptr{UInt8}), ih, name, lin, col, value)
end
function IupMatGetAttribute(ih::Ptr{Int}, name::String, lin::Int, col::Int)
	ccall((:IupMatGetAttribute, libiupcontrols_), Ptr{UInt8}, (Ptr{Cint}, Ptr{UInt8}, Cint, Cint), ih, name, lin, col)
end
function IupMatGetInt(ih::Ptr{Int}, name::String, lin::Int, col::Int)
	ccall((:IupMatGetInt, libiupcontrols_), Cint, (Ptr{Cint}, Ptr{UInt8}, Cint, Cint), ih, name, lin, col)
end
function IupMatGetFloat(ih::Ptr{Int}, name::String, lin::Int, col::Int)
	ccall((:IupMatGetFloat, libiupcontrols_), Cfloat, (Ptr{Cint}, Ptr{UInt8}, Cint, Cint), ih, name, lin, col)
end


text:
struct InternalArray{T,N,P} <: AbstractArray{T,N}
	parent::P
	extracted::Int
end

#TODO check the dimensions and everything is as expected
# That is: parent[i][extracted] is meaningful
function InternalArray(parent::AbstractArray, extracted::Int)
	InternalArray{eltype(parent), ndims(eltype(parent)), typeof(parent)}(parent, extracted)
end

Base.size(v::InternalArray) = size(v.parent)
Base.getindex(v::InternalArray, i) = v.parent[i][v.extracted]

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
# NOTE: This file contains many parts that are copied from the file runtests.jl from the Package MPI.jl.
push!(LOAD_PATH, "../src")

import ParallelStencil # Precompile it.
import ParallelStencil: SUPPORTED_PACKAGES, PKG_CUDA
@static if (PKG_CUDA in SUPPORTED_PACKAGES) import CUDA end

excludedfiles = [ "test_excluded.jl"];

function runtests()
    exename   = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir   = pwd()
    istest(f) = endswith(f, ".jl") && startswith(basename(f), "test_")
    testfiles = sort(filter(istest, vcat([joinpath.(root, files) for (root, dirs, files) in walkdir(testdir)]...)))

    nfail = 0
    printstyled("Testing package ParallelStencil.jl\n"; bold=true, color=:white)

    if (PKG_CUDA in SUPPORTED_PACKAGES && !CUDA.functional())
        @warn "Test Skip: All CUDA tests will be skipped because CUDA is not functional (if this is unexpected type `import CUDA; CUDA.functional(true)` to debug your CUDA installation)."
    end

    for f in testfiles
        println("")
        if f ∈ excludedfiles
            println("Test Skip:")
            println("$f")
            continue
        end
        try
            run(`$exename -O3 --startup-file=no --check-bounds=no $(joinpath(testdir, f))`)
        catch ex
            nfail += 1
        end
    end
    return nfail
end

exit(runtests())


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerMulticomponentEquations1D(gammas           = (1.4, 1.4, 1.4),
                                                       gas_constants    = (0.4, 0.4, 0.4))

"""
    initial_condition_two_interacting_blast_waves(x, t, equations::CompressibleEulerMulticomponentEquations1D)

A multicomponent two interacting blast wave test taken from
- T. Plewa & E. Müller (1999)
  The consistent multi-fluid advection method
  [arXiv: 9807241](https://arxiv.org/pdf/astro-ph/9807241.pdf)
"""
function initial_condition_two_interacting_blast_waves(x, t, equations::CompressibleEulerMulticomponentEquations1D)

  rho1        = 0.5 * x[1]^2
  rho2        = 0.5 * (sin(20 * x[1]))^2
  rho3        = 1 - rho1 - rho2

  prim_rho    = SVector{3, real(equations)}(rho1, rho2, rho3)

  v1          = 0.0

  if x[1] <= 0.1
    p = 1000
  elseif x[1] < 0.9
    p = 0.01
  else
    p = 100
  end

  prim_other  = SVector{2, real(equations)}(v1, p)

  return prim2cons(vcat(prim_other, prim_rho), equations)
end
initial_condition = initial_condition_two_interacting_blast_waves

function boundary_condition_two_interacting_blast_waves(u_inner, orientation, direction, x, t,
                                                        surface_flux_function,
                                                        equations::CompressibleEulerMulticomponentEquations1D)

  u_inner_reflect = SVector(-u_inner[1], u_inner[2], u_inner[3], u_inner[4], u_inner[5])
  # Calculate boundary flux
  if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_inner_reflect, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_inner_reflect, u_inner, orientation, equations)
  end

  return flux
end
boundary_conditions = boundary_condition_two_interacting_blast_waves

surface_flux = flux_lax_friedrichs
volume_flux  = flux_ranocha
basis = LobattoLegendreBasis(3)
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max = 0.8,
                                         alpha_min = 0.0,
                                         alpha_smooth = true,
                                         variable=pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (0.0,)
coordinates_max = (1.0,)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=9,
                n_cells_max=10_000,
                periodicity=false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.038)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100

analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=0.1)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks, maxiters=1e5);
summary_callback() # print the timer summary