text:
using CSV, DataFrames, Statistics
using Printf
df = CSV.read(joinpath(@__DIR__, "results.csv"))

## This bit is a little annoying:
# To check for OOMs, we need to manually look through the log (see `log.md`)
# Which was numbered according to the algorithm and problem index.
# However, later two more algorithms were added (`MISDP_dB` and `MISDP (dB^2)`)
# which can mess up the numbering. So first we will filter out those to mimick
# the state of the first run, and add flags for OOMs.
include("common.jl") # need to load the list of algorithms (`algos`)

algos_pre_dB = filter(algos) do a
    a.algo != "MISDP_dB" && a.algo != "MISDP (dB^2)"
end

# Looking through `log.md`, we find two issues:
# out of memory (OOM) errors during 4 solves, and unknown errors in 2 solves.
OOM = [(problem=10, algo=4), (problem=10, algo=5), (problem=10, algo=9),
       (problem=10, algo=10)]
unknown_error = [(problem=12, algo=4), (problem=12, algo=9)]

# Let us record the OOM errors in the dataframe.
inds = falses(size(df, 1))
for (p_idx, a_idx) in OOM
    inds .= inds .| ((df.problem .== problems[p_idx].problem) .& (df.algo .== algos_pre_dB[a_idx].algo) .& (df.settings .== algos_pre_dB[a_idx].settings))
end
df.OOM = inds
## End of problem-index specific content.

# Now we will filter out the `MISDP` algorithm (which had `dB^2 + 1` entries, which
# is 1 more than needed, and is replaced by the `MISDP (dB^2)` version.
df = filter(df) do row
    row.algo != "MISDP"
end

df.timeout = isnan.(df.optval) .& (!).(df.errors)

df.error_not_solved_not_timeout = df.errors .& isnan.(df.optval) .& (!).(df.timeout)
# Now we will caculate discrepencies in `optval`
# First we will filter to rows that claim to have a correct answer
row_should_be_correct(row) = !startswith(row.algo, "guesswork_upper_bound") && !isnan(row.optval) && row.algo != "MISDP_dB"

df_right = filter(row_should_be_correct, df)
gdf_right = groupby(df_right, "problem")

# Check for BB84(1)
bb84_true_val = (big(1) / big(4)) * (10 - sqrt(big(10)))

@assert all(gdf_right[4].problem .== Ref("BB84(1)")) # check we've got the right group
@assert maximum((gdf_right[4].optval .- bb84_true_val ) ./ bb84_true_val) < 1e-7 # not too big of an error

mean_solutions_df = combine(gdf_right, :optval => mean)
mean_solutions = Dict(mean_solutions_df.problem .=> mean_solutions_df.optval_mean)

# Now that we've calculated the mean solutions, we calculate the discrepencies

df.discrepency_from_mean = [ isnan(row.optval) ? NaN : !haskey(mean_solutions, row.problem) ? missing : abs(row.optval - mean_solutions[row.problem]) for row in eachrow(df)]

df.relative_discrepency_from_mean = [ isnan(row.optval) ? NaN : !haskey(mean_solutions, row.problem) ? missing : abs(row.optval - mean_solutions[row.problem])/mean_solutions[row.problem] for row in eachrow(df)]

# Let us check the maximum discrepency is not too bad
df_right = filter(row_should_be_correct, df)
findmax(df_right.relative_discrepency_from_mean) # (1.7202621117336732e-6, 45)

# Now let us generate the final table.
## Rows in table: algorithm with settings
## Columns in table: # solved, # timeout, # errored, avg relative_discrepency_from_mean, average time

df.algo_settings = df.algo .* Ref(" ") .* df.settings

function pretty_algo(row)
    if startswith(row.algo, "guesswork_upper_bound")
        return "Upper bound"
    end
    if row.algo == "MISDP_dB"
        return "MISDP (\$d_B\$)"
    elseif row.algo == "MISDP"
        return "MISDP (\$d_B^2 + 1\$)"
    elseif row.algo == "MISDP (dB^2)"
        return "MISDP (\$d_B^2\$)"
    elseif row.algo == "dual_SDP"
        return "SDP (dual)"
    end
    return row.algo
end

function pretty_settings(row)
    n = length("guesswork_upper_bound")

    if startswith(row.algo, "guesswork_upper_bound")
        j = findfirst(==('='), row.algo)
        max_time = row.algo[j+1:end-1]
        return row.settings * ", \$t_\\text{max}=$(max_time)\$"
    end
    if row.settings == "Pajarito(Mosek, Gurobi, MSD=true)"
        return "Pajarito (c1)"
    elseif row.settings == "Pajarito(Mosek, Gurobi, MSD=false)"
        return "Pajarito (c2)"
    elseif row.settings == "Pajarito(SCS, Cbc, MSD=false)"
        return "Pajarito (o)"
    end
    return row.settings
end

df.algo2 = [ pretty_algo(row) for row in eachrow(df) ]
df.settings2 =  [ pretty_settings(row) for row in eachrow(df) ]

gdf = groupby(df, [:algo2, :settings2])

function mean_no_nan(collection)
    c = [x for x in collection if  !ismissing(x) && !isnan(x)]
    isempty(c) ? missing : mean(c)
end
float_fmt(f) = @sprintf("%2.2f", f)
pct_fmt(f) =  float_fmt(f*100) * "\\,\\%"
time_fmt(f) =  float_fmt(f) * "\\,s"

table1 = combine(gdf, :relative_discrepency_from_mean => pct_fmt ∘ mean_no_nan => "average relative error", :elapsed_seconds => time_fmt ∘ mean_no_nan => "average time", :optval => (x -> sum((!isnan).(x))) => "solved", :timeout => sum => "timeouts", :error_not_solved_not_timeout => sum => "errors")

sort!(table1, [:algo2, :settings2])
table1 = rename(table1, :algo2 => :Algorithm, :settings2 => :Parameters)

df.time_status = ifelse.(df.timeout, Ref(:timeout), ifelse.(df.error_not_solved_not_timeout, Ref(:error), df.elapsed_seconds))

df2 = copy(df)

function tol_time(row)
    disc = ismissing(row.relative_discrepency_from_mean) ? "(?\\,\\%)" : isnan(row.relative_discrepency_from_mean) ? "" : "("* pct_fmt(row.relative_discrepency_from_mean)*")"

    time = row.time_status isa Symbol ? string(row.time_status) : time_fmt(row.elapsed_seconds)
    disc == "" ? time : "$time $disc"
end

df2.value = [tol_time(row) for row in eachrow(df2) ]
df_times = unstack(df2, [:algo2, :settings2], :problem, :value)
problem_stems =  unique([ str[1:end-3] for str in df2.problem])

nms = names(df_times, Not(["algo2", "settings2"]))

for stem in problem_stems
    df_times[:, stem] =  ["$(row[1]), $(row[2])" for row in eachrow(df_times[:, nms[findall(startswith(stem), nms)]])]
end
table2 = df_times[:, ["algo2", "settings2", problem_stems...]]
sort!(table2, [:algo2, :settings2])
table2 = rename(table2, :algo2 => :Algorithm, :settings2 => :Parameters)

using Latexify
latexify(table1, env=:table, latex=false)
latexify(table2[:,["Algorithm", "Parameters", "2qubits", "2qutrits", "Y"] ], env=:table, latex=false)
latexify(table2[:,["Algorithm", "Parameters", "3qubits", "3qutrits", "BB84"] ], env=:table, latex=false)


text:
using WebBasedWorkspace
using Test
using URIs
using Pkg

EXAMPLE_URI = URI("https://github.com/MarkNahabedian/WebBasedWorkspace")

@testset "WebBasedWorkspace.jl" begin
    pspecs = package_specs(github_raw(EXAMPLE_URI))
    @test length(pspecs) == 4
    Pkg.add(pspecs)
end

ext: jl

max_stars_repo_name: MarkNahabedian/WebBasedWorkspace

max_issues_repo_name: MarkNahabedian/WebBasedWorkspace

max_forks_repo_name: MarkNahabedian/WebBasedWorkspace


text:
@eval Discord struct Baz
    id::Int
    x::Int
    y::Union{Int, Missing}
end
@eval Discord @boilerplate Baz :merge
using Discord: Baz

@testset "State" begin
    @testset "insert_or_update!" begin
        d = Dict()
        v = Baz[]

        # Inserting a new entry works fine.
        b = Baz(1, 2, 3)
        insert_or_update!(d, b.id, b)
        @test d[b.id] == b
        # But when we insert a value with the same ID, it's merged.
        b = Baz(1, 3, 3)
        insert_or_update!(d, b.id, b)
        @test d[b.id] == b
        b = Baz(1, 10, missing)
        insert_or_update!(d, b.id, b)
        @test d[b.id] == Baz(b.id, b.x, 3)

        # We can also insert/update without specifying the key, :id is inferred.
        b = Baz(2, 2, 1)
        insert_or_update!(d, b)
        @test d[b.id] == b
        b = Baz(2, 1, missing)
        insert_or_update!(d, b)
        @test d[b.id] == Baz(b.id, b.x, 1)

        # We can also do this on lists.
        b = Baz(1, 2, 3)
        insert_or_update!(v, b.id, b)
        @test length(v) == 1 && first(v) == b
        b = Baz(1, 10, missing)
        insert_or_update!(v, b.id, b)
        @test length(v) == 1 && first(v) == Baz(b.id, b.x, 3)
        b = Baz(1, 2, 2)
        insert_or_update!(v, b)
        @test length(v) == 1 && first(v) == b

        # If we need something fancier to index into a list, we can use a key function.
        b = Baz(1, 4, 5)
        insert_or_update!(v, b.id, b; key=x -> x.y)
        # We inserted a new element because we looked for an element with y == b.id.
        @test length(v) == 2 && last(v) == b
        # If we leave out the insert key, the key function is used on the value too.
        b = Baz(1, 0, 5)
        # We updated the value with y == b.y.
        insert_or_update!(v, b; key=x -> x.y)
        @test length(v) == 2 && last(v) == b

        # The updated value is returned.
        empty!(d)
        empty!(v)
        b = Baz(1, 2, 3)
        @test insert_or_update!(d, b) == b
        @test insert_or_update!(v, b) == b
        b = Baz(1, 3, missing)
        @test insert_or_update!(d, b) == Baz(1, 3, 3)
        @test insert_or_update!(v, b) == Baz(1, 3, 3)
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.4
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# # Hamiltonian Monte Carlo with leapfrog
#
# Scalar version: https://github.com/genkuroki/public/blob/main/0018/HMC%20leapfrog.ipynb

# %% tags=[]
module My

using ConcreteStructs: @concrete
using Parameters: @unpack

using LinearAlgebra: dot
using ForwardDiff: gradient
using Random: default_rng, randn!
using StaticArrays: SVector, MVector

@concrete struct LFProblem{dim} ϕ; H; F; dt; nsteps end

"""Assume ϕ is a potential function."""
function LFProblem(dim, ϕ; dt = 1.0, nsteps = 40)
    H(x, v, param) = dot(v, v)/2 + ϕ(x, param)
    F(x, param) = -gradient(x -> ϕ(x, param), x)
    LFProblem{dim}(ϕ, H, F, dt, nsteps)
end

"""Assume ϕ is a potential function and ∇ϕ its gradient."""
function LFProblem(dim, ϕ, ∇ϕ; dt = 1.0, nsteps = 40)
    H(x, v, param) = dot(v, v)/2 + ϕ(x, param)
    F(x, param) = -∇ϕ(x, param)
    LFProblem{dim}(ϕ, H, F, dt, nsteps)
end

"""Numerically solve Hamilton's equation of motion with leapfrog method"""
function solve(lf::LFProblem, x, v, param)
    @unpack F, dt, nsteps = lf
    v = v + F(x, param)*dt/2
    x = x + v*dt
    for _ in 2:nsteps
        v = v + F(x, param)*dt
        x = x + v*dt
    end
    v = v + F(x, param)*dt/2
    x, v
end

@inline function _update!(lf::LFProblem{dim}, x, vtmp, param, rng) where dim
    @unpack H = lf
    v = SVector{dim}(randn!(rng, vtmp))
    xnew, vnew = solve(lf, x, v, param)
    dH = H(xnew, vnew, param) - H(x, v, param)
    rand(rng) ≤ exp(-dH) ? xnew : x
end

"""Hamiltonian Monte Carlo"""
function HMC(lf::LFProblem{dim}, param = nothing;
        niters = 10^5, thin = 1, nwarmups = 0, rng = default_rng(),
        init = SVector{dim}(randn(rng, dim))) where dim
    vtmp = MVector{dim}(zeros(eltype(init), dim))
    x = init
    for _ in 1:nwarmups
        x = _update!(lf, x, vtmp, param, rng)
    end
    sample = Vector{typeof(init)}(undef, niters)
    for i in 1:niters
        for _ in 1:thin
            x = _update!(lf, x, vtmp, param, rng)
        end
        @inbounds sample[i] = x
    end
    sample
end

end

# %%
using Plots
using BenchmarkTools
using StaticArrays
using LinearAlgebra
using KernelDensity
using Statistics
using QuadGK
using Distributions
using Symbolics

# %% [markdown]
# ## 2-dimensional normal distribution

# %%
A = @SMatrix [
     1  1/2
    1/2  1
]
param = (; A = A)
ϕ(x, param) = dot(x, param.A, x)/2
lf = My.LFProblem(2, ϕ)

# %%
@time sample = My.HMC(lf, param)
@time sample = My.HMC(lf, param)
@time sample = My.HMC(lf, param);

# %%
@btime My.HMC($lf, $param);

# %%
X, Y = first.(sample), last.(sample)
d = InterpKDE(kde((X, Y)))
x, y = range(extrema(X)...; length=201), range(extrema(Y)...; length=201)
heatmap(x, y, (x, y) -> pdf(d, x, y); size=(450, 400), right_margin=3Plots.mm)

# %%
f(n) = mean(x -> x*x', @view sample[1:n])

n = 1:1000
S = f.(n)
S11 = (S -> S[1,1]).(S)
S22 = (S -> S[2,2]).(S)
S12 = (S -> S[1,2]).(S)

ymin = min(-1.5, minimum(S11), minimum(S22), minimum(S12))
ymax = max(2.5, maximum(S11), maximum(S22), maximum(S12))

plot(ylim = (ymin, ymax))
plot!(S11; label="s11", c=1)
hline!([inv(A)[1,1]]; label="", c=1, ls=:dash)
plot!(S22; label="s22", c=2)
hline!([inv(A)[2,2]]; label="", c=2, ls=:dash)
plot!(S12; label="s12", c=3)
hline!([inv(A)[1,2]]; label="", c=3, ls=:dash)

# %% [markdown]
# ## φ(x) = a(x² - 1)²

# %%
ϕ4(x, a) = a * (x[1]^2 - 1)^2
a = [3, 4, 5, 6, 7, 8]
XX = Vector{Float64}[]
ZZ = Float64[]
PP = []
for i in eachindex(a)
    Z = quadgk(x -> exp(-ϕ4((x,), a[i])), -Inf, Inf)[1]
    push!(ZZ, Z)
    lf = My.LFProblem(1, ϕ4; dt = 0.05, nsteps = 100)
    @time X = first.(My.HMC(lf, a[i]))
    flush(stdout)
    push!(XX, X)
    P = plot()
    histogram!(X; norm=true, alpha=0.3, label="HMC LF sample", bin=100, c=i)
    plot!(x -> exp(-ϕ4(x, a[i]))/Z, -2, 2; label="exp(-ϕ2(x))/Z", lw=2, c=i)
    plot!(; legend=false, xtick=-2:0.5:2)
    title!("ϕ(x) = a(x² - 1)²  for  a = $(a[i])", titlefontsize=9)
    push!(PP, P)
end
plot(PP...; size=(800, 450), layout=(3, 2))

# %%
QQ = []
for i in eachindex(a)
    Q = plot(XX[i][1:10000]; ylim=(-1.5, 1.5), label="", c=i, lw=0.5)
    title!("ϕ(x) = a(x² - 1)²  for  a = $(a[i])", titlefontsize=9)
    push!(QQ, Q)
end
plot(QQ...; size=(800, 900), layout=(length(a), 1))

# %% [markdown]
# ## Baysian inference for a sample of the standard normal distribution

# %%
n = 10
sample_normal = randn(n)
f(y, m, s) = (y - m)^2/(2s^2) + log(s)
negloglik(w, sample) = sum(y -> f(y, w[1], exp(w[2])), sample)
lf = My.LFProblem(2, negloglik; dt = 0.1, nsteps = 30)

# %%
@time sample = My.HMC(lf, sample_normal; init = SVector(0.0, 0.0))
@time sample = My.HMC(lf, sample_normal; init = SVector(0.0, 0.0))
@time sample = My.HMC(lf, sample_normal; init = SVector(0.0, 0.0))

# %%
m, logs = first.(sample), last.(sample)
d = InterpKDE(kde((m, logs)))
x, y = range(extrema(m)...; length=201), range(extrema(logs)...; length=201)
heatmap(x, y, (x, y) -> pdf(d, x, y); size=(450, 400), xlabel="μ", ylabel="log(σ)")

# %% [markdown]
# ## Symbolics.jl example

# %%
dim = 2
@variables a[1:dim, 1:dim], x[1:dim]
aa, xx = collect.((a, x))

# %%
phi_sym = dot(xx, aa, xx)/2 |> expand |> simplify

# %%
dphi_sym = Symbolics.gradient(phi_sym, xx)

# %%
phi_rgf = build_function(phi_sym, xx, aa; expression = Val(false))
dphi_rgf = build_function(dphi_sym, xx, aa; expression = Val(false))[1]
lf = My.LFProblem(dim, phi_rgf, dphi_rgf);

# %%
param = A = @SMatrix [
     1  1/2
    1/2  1
]

# %%
@time sample = My.HMC(lf, param)
@time sample = My.HMC(lf, param)
@time sample = My.HMC(lf, param);

# %%
@btime My.HMC($lf, $param);

# %%
X, Y = first.(sample), last.(sample)
d = InterpKDE(kde((X, Y)))
x, y = range(extrema(X)...; length=201), range(extrema(Y)...; length=201)
heatmap(x, y, (x, y) -> pdf(d, x, y); size=(450, 400), right_margin=3Plots.mm)

# %%
f(n) = mean(x -> x*x', @view sample[1:n])

n = 1:1000
S = f.(n)
S11 = (S -> S[1,1]).(S)
S22 = (S -> S[2,2]).(S)
S12 = (S -> S[1,2]).(S)

ymin = min(-1.5, minimum(S11), minimum(S22), minimum(S12))
ymax = max(2.5, maximum(S11), maximum(S22), maximum(S12))

plot(ylim = (ymin, ymax))
plot!(S11; label="s11", c=1)
hline!([inv(A)[1,1]]; label="", c=1, ls=:dash)
plot!(S22; label="s22", c=2)
hline!([inv(A)[2,2]]; label="", c=2, ls=:dash)
plot!(S12; label="s12", c=3)
hline!([inv(A)[1,2]]; label="", c=3, ls=:dash)

# %%

max_stars_repo_path: 0020/HMC leapfrog Symbolics.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0020/HMC leapfrog Symbolics.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0020/HMC leapfrog Symbolics.jl

max_forks_repo_name: genkuroki/public


text:
module PSFDistiller

using IndexFunArrays
using Images # find_local_maxima
using NearestNeighbors # to discard them in PSF extraction
using InverseModeling # for gauss_fit
using FourierTools # for shift
using FFTW # for rfft
using NDTools # for select_region

include("psf_extraction.jl")
include("utils.jl")

end # module


text:
using JuMP
using StructJuMP

include("Benders_pmap.jl")

#===================================================
 The following function is taken from JuMP src/solvers.jl:conicdata().
 Modified for extension StructJuMP/Benders_pmap.jl.
===================================================#
function conicconstraintdata(m::Model)
    # TODO this generates redundant constraints for variables that are already bounded (one-sided)
    stoch = getStructure(m)
    parent = stoch.parent
    numMasterCols = 0
    if parent !== nothing
        numMasterCols = parent.numCols
    end
    v = Symbol[]

    var_cones = Any[cone for cone in m.varCones]
    con_cones = Any[]
    nnz = 0

    numSDPRows, numSymRows, nnz = JuMP.getSDrowsinfo(m)

    linconstr::Vector{LinearConstraint} = m.linconstr
    numLinRows = length(linconstr)

    numBounds = JuMP.variable_range_to_cone!(var_cones, m)

    nnz += numBounds
    for c in 1:numLinRows
        nnz += length(linconstr[c].terms.coeffs)
    end

    numSOCRows = JuMP.getNumSOCRows(m)
    numNormRows = length(m.socconstr)

    numRows = numLinRows + numBounds + numSOCRows + numSDPRows + numSymRows

    # constr_to_row is not used but fill_bounds_constr! and fillconstr! for SDP needs them
    constr_to_row = Array{Vector{Int}}(numBounds + 2*length(m.sdpconstr))

    b = Array{Float64}(numRows)

    I_m = Int[]
    J_m = Int[]
    V_m = Float64[]
    I_s = Int[]
    J_s = Int[]
    V_s = Float64[]

    # Fill it up
    if numMasterCols > 0
        tmprow_m = JuMP.IndexedVector(Float64, parent.numCols)
    end
    tmprow_s = JuMP.IndexedVector(Float64, m.numCols)

    JuMP.fillconstrRHS!(b, con_cones, 0, m.linconstr)
    if numMasterCols > 0
        JuMP.fillconstrLHS!(I_m, J_m, V_m, tmprow_m, 0, m.linconstr, parent, true)
    end
    c = JuMP.fillconstrLHS!(I_s, J_s, V_s, tmprow_s, 0, m.linconstr, m, true)

    for idx in 1:m.numCols
        # identify integrality information
        push!(v, m.colCat[idx])
    end
    c, d = JuMP.fill_bounds_constr!(I_s, J_s, V_s, b, con_cones, constr_to_row, c, 0, m)

    @assert c == numLinRows + numBounds
    @assert d == numBounds

    JuMP.fillconstrRHS!(b, con_cones, c, m.socconstr)
    if numMasterCols > 0
        JuMP.fillconstrLHS!(I_m, J_m, V_m, tmprow_m, c, m.socconstr, parent, true)
    end
    c = JuMP.fillconstrLHS!(I_s, J_s, V_s, tmprow_s, c, m.socconstr, m, true)

    @assert c == numLinRows + numBounds + numSOCRows

    if numMasterCols > 0
        c, d = JuMP.fillconstr!(I_m, J_m, V_m, b, con_cones, tmprow_m, constr_to_row, c, d, m.sdpconstr, m, true)
    end
    c, d = JuMP.fillconstr!(I_s, J_s, V_s, b, con_cones, tmprow_s, constr_to_row, c, d, m.sdpconstr, m, true)

    if c < length(b)
        # This happens for example when symmetry constraints are dropped with SDP
        resize!(b, c)
    end

    f_s = JuMP.prepAffObjective(m)

    # The conic MPB interface defines conic problems as
    # always being minimization problems, so flip if needed
    m.objSense == :Max && scale!(f_s, -1.0)

    if numMasterCols > 0
        JuMP.rescaleSDcols!(spzeros(numMasterCols), J_m, V_m, parent)
    end
    JuMP.rescaleSDcols!(f_s, J_s, V_s, m)

    A = sparse(I_m, J_m, V_m, numRows, numMasterCols)
    B = sparse(I_s, J_s, V_s, numRows, m.numCols)

    return f_s, A, B, b, var_cones, con_cones, v
end

function BendersBridge(m::Model, master_solver, sub_solver)

    c_all = Vector{Float64}[]
    A_all = SparseMatrixCSC{Float64}[]
    B_all = SparseMatrixCSC{Float64}[]
    b_all = Vector{Float64}[]
    K_all = Any[]
    C_all = Any[]
    v_all = Symbol[]

    (c,A,B,b,var_cones, constr_cones, v) = conicconstraintdata(m)
    push!(c_all, c)
    push!(A_all, B)
    push!(b_all, b)
    push!(K_all, constr_cones)
    push!(C_all, var_cones)
    append!(v_all, v)

    for i = 1:num_scenarios(m)
        (c,A,B,b,var_cones, constr_cones, v) = conicconstraintdata(getchildren(m)[i])
        push!(c_all, getprobability(m)[i] * c)
        push!(A_all, A)
        push!(B_all, B)
        push!(b_all, b)
        push!(K_all, constr_cones)
        push!(C_all, var_cones)
        append!(v_all, v)
    end

    return Benders_pmap(c_all,A_all,B_all,b_all,K_all,C_all,v_all,master_solver,sub_solver)
end

function DLP(m::Model, solver)

    c_all = Vector{Float64}[]
    A_all = SparseMatrixCSC{Float64}[]
    B_all = SparseMatrixCSC{Float64}[]
    b_all = Vector{Float64}[]
    K_all = Any[]
    C_all = Any[]
    v_all = Symbol[]

    (c,A,B,b,var_cones, constr_cones, v) = conicconstraintdata(m)
    push!(c_all, c)
    push!(A_all, B)
    push!(b_all, b)
    push!(K_all, constr_cones)
    push!(C_all, var_cones)
    append!(v_all, v)

    for i = 1:num_scenarios(m)
        (c,A,B,b,var_cones, constr_cones, v) = conicconstraintdata(getchildren(m)[i])
        push!(c_all, getprobability(m)[i] * c)
        push!(A_all, A)
        push!(B_all, B)
        push!(b_all, b)
        push!(K_all, constr_cones)
        push!(C_all, var_cones)
        append!(v_all, v)
    end
    # v_all is unused for now

    dlp_rows = sum(length(b) for b in b_all)
    dlp_cols = sum(length(c) for c in c_all)

    c_dlp = zeros(dlp_cols)
    A_dlp = spzeros(dlp_rows, dlp_cols)
    b_dlp = zeros(dlp_rows)
    K_dlp = Any[] # constr_cones
    C_dlp = Any[] # var_cones

    rows = 1
    cols = 1
    for i = 1:length(c_all)
        rows_end = rows + length(b_all[i]) - 1
        cols_end = cols + length(c_all[i]) - 1

        c_dlp[cols:cols_end] = c_all[i]
        b_dlp[rows:rows_end] = b_all[i]

        A_dlp[rows:rows_end, 1:size(A_all[i], 2)] = A_all[i]
        if i != 1
            A_dlp[rows:rows_end, cols:cols_end] = B_all[i-1]
        end

        append!(K_dlp, [(cone, rows-1 + idx) for (cone, idx) in K_all[i]])
        append!(C_dlp, [(cone, cols-1 + idx) for (cone, idx) in C_all[i]])
        rows = rows_end+1
        cols = cols_end+1
    end

    model = MathProgBase.ConicModel(solver)
    MathProgBase.loadproblem!(model, c_dlp, A_dlp, b_dlp, K_dlp, C_dlp)

    # solve conic model
    MathProgBase.optimize!(model)
    status = MathProgBase.status(model)

    # return status, objective value, and solution
    return status, MathProgBase.getobjval(model), MathProgBase.getsolution(model)[1:length(c_all[1])]

end


text:
""" Julia Program to implement the Cutting a Rod problem

We are given a rod of length 'n' and the cost of the rod of different length from
1 to n. We need to find the maximum sell price that can be obtained by cutting the rod
at any lengths.

This problem can be solved by using Dynamic Programming
"""

# Build a dp Table in Botom-Up manner, such that the n'th element stores
# the maximum price obtained by a rod of length 'n'
function max_sell_price(m, cost)
    dp = zeros(Int, n + 1)
    for i in 1:n
        curr_max = -1
        for j in 0:(i - 1)
            # Maximum Price obtained by cutting a rod of length zero is zero
            if ((i - j - 1) == 0)
                curr_max = max(curr_max, cost[j + 1])
            else
                curr_max = max(curr_max, cost[j + 1] + dp[i - j - 1])
            end
        end
        #  Stores the maximum cost obtainable for a rod of length i
        dp[i] = curr_max
    end
    return dp[n]
end

print("What is the length of the rod? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("The given rod length is Invalid!!!")
    exit()
end
cost = Int[]
print("Enter the cost of the rod for length 1 to $n :  ")
cost = [parse(Int, num) for num in split(readline())]
max_price = max_sell_price(n, cost)
print("The maximized sell price of the: $max_price")

"""
Time Complexity- O(n^2), where 'n' is the length of the rod
Space Complexity- O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE 1

What is the length of the rod? 5
Enter the cost of the rod for length 1 to 5: 2 5 7 8 10
The maximized sell price of the: 12

SAMPLE 2

What is the length of the rod? -3
The given rod length is Invalid!!!

"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
function make_results_df(times::AbstractVector{T}, rvs_clean::AbstractVector{T}, σ_rvs_clean::AbstractVector{T},
              rvs_dirty::AbstractVector{T} = Vector{T}(undef,0), σ_rvs_dirty::AbstractVector{T} = Vector{T}(undef,0),
              scores::AbstractArray{T,2} = Array{T,2}(undef,0,0), σ_scores::AbstractArray{T,2} = Array{T,2}(undef,0,0),
              #unit_scores::AbstractVector{String} = fill("nondim",size(scores,2))
              ) where { T<:Real }
  @assert length(times) == length(rvs_clean) == length(σ_rvs_clean)
  @assert length(rvs_dirty) == length(σ_rvs_dirty)
  output_df = DataFrame("Time [MJD]"=>times, "RV_C [m/s]"=>rvs_clean, "e_RV_C [m/s]"=>σ_rvs_clean,
                    "RV_D [m/s]"=>rvs_dirty, "e_RV_D [m/s]"=>σ_rvs_dirty )
  if length(rvs_dirty) == length(times)
    output_df[!,"RV_D [m/s]"] = rvs_dirty
    output_df[!,"e_RV_D [m/s]"] = σ_rvs_dirty
  end
  if size(scores,1) >= 1
    @assert length(times) == size(scores,1)
    if size(σ_scores,2) == 0
      σ_scores = fill(NaN,length(times),size(scores,2)  )
    end
    @assert length(times) == size(σ_scores,1)
    for i in 1:size(scores,2)
      output_df[!,"Score_" * string(i) * " [nondim]"] = vec(scores[:,i])
      output_df[!,"e_Score_" * string(i) * " [nondim]"] = vec(σ_scores[:,i])
    end
  end
  return output_df
end


text:
@testset "@ivp for a continuous system" begin
    P1 = @ivp(x' = -x, x(0) ∈ Interval(-1.0, 1.0))
    testSystem = InitialValueProblem(LinearContinuousSystem(I(-1.0, 1)), Interval(-1, 1))
    @test P1 == testSystem
    # throw error if x(0) is not defined
    @test_throws ArgumentError @ivp(x' = -x)

    sys = @system(x' = -x)
    P2 = @ivp(sys, x(0) ∈ Interval(-1.0, 1.0))
    @test P2 == testSystem
    P3 = @ivp(@system(x' = -x), x(0) ∈ Interval(-1.0, 1.0))
    @test P3 == testSystem
end


text:
function evaluate(p::MPoly{T}, es...) where {T}
    r = zero(T)
    for (m, c) in p
        t = c
        for i = 1:length(m)
            t *= es[i]^m[i]
        end
        r += t
    end
    r
end

function evaluate_basis(p::MPoly{T}, es...) where {T}
    r = Array{T}(undef, 0)
    for (m, c) in p
        r1 = c
        for i = 1:length(m)
            r1 *= es[i]^m[i]
        end
        push!(r, r1)
    end
    r
end


text:
### Quantum.jl
# Core functionality for setting up quantum Hilbert spaces
# of operators, handling tensor products, and computing
# partial traces
###

###
# Quantum Hilbert Space Factor
###
"""
    QFactor(dim::Int, name::QName[, ops::QOps])

Single Hilbert space factor for a quantum space.

### Fields:
  - dim : Hilbert space dimension
  - name : String naming the factor for clarity
  - ops  : Dict of named operators belonging to space

### Bases:
  - Indexing as an array with [] indexes from 1
  - Indexing as a function with () indexes from 0
"""
struct QFactor <: QObj
    dim  :: Int
    name :: QName
    ops  :: QOps
    # Input sanity check
    function QFactor(d::Int, n::QName, ops::QOps)
        (d > 0) || error("Dimension must be positive.")
        new(d,n,ops)
    end
end
# Simplified constructor that auto-creates identity operator
QFactor(dim::Int, name::QName) =
    QFactor(dim, name,
            QOps('i' => sparse(QComp(1.0)SparseArrays.I, dim, dim))
            )

# Helper functions
size(s :: QFactor) = s.dim
length(s :: QFactor) = s.dim
name(s :: QFactor) = s.name

# Overload function call of string to get operator
(q::QFactor)(arg::QOpName) = q.ops[arg]

###
# Quantum Tensor Product of Hilbert Spaces
###
"""
    QSpace(dim::Int, name::QName)
    QSpace(f::QFactor)
    QSpace(fs::Vector{QFactor}, ops::QOps)

Tensor product containing quantum Hilbert space factors.

### Fields:
  - factors : Vector of Hilbert space factors
  - ops     : Dict of named operators belonging to space

### Bases:
  - Indexing as an array with [] indexes from 1
  - Indexing as a function with () indexes from 0
"""
struct QSpace <: QObj
    factors :: Vector{QFactor}
    ops     :: QsOps
end
# Convenient constructors
QSpace(f::QFactor) = QSpace(vec([f]),
                            QsOps([(string(k),v) for (k,v) in f.ops]))
QSpace(dim::Int, name::QName) = QSpace(QFactor(dim, name))

# Helper functions
factors(s :: QSpace) = s.factors
@inline size(s :: QSpace) = tuple(map(size, factors(s))...)
length(s :: QSpace) = prod(size(s))
name(s :: QSpace) = join(map(name, factors(s)), " ⊗ ")

# Overload function call
# Memoize creation of tensor product of composite operators for efficiency
(q::QSpace)(arg::AbstractString) = begin
    if haskey(q.ops, arg)
        q.ops[arg]
    elseif length(arg) == length(factors(q))
        q.ops[arg] = ⊗([q.factors[i](c) for (i,c) in enumerate(arg)]...)
    else
        error("QSpace does not contain named operator")
    end
end

###
# Pretty printing
###
function show(io::IO, q::QObj)
    println(io, "$(typeof(q)): ", name(q))
    println(io, "Dims  : ", size(q))
    println(io, "Ops   : '", join(keys(q.ops), "', '"), "'")
end

###
# Subviews
###
"""
    subview(s::QObj, op::QOp)

Return a reshaped view on quantum operator that makes
subsystem indices transparent.

`t[i,j]` is equivalent to ``T_{i;j} = a_{i} ⊗ b_{j}``.

`t[i,k,j,l]` is equivalent to ``T_{ik;jl} = a_{ij} ⊗ b_{kl}``.

**Example:**
```julia
a = ["a11" "a12" ; "a21" "a22"];
b = ["b11" "b12" ; "b21" "b22"];
c = ["c11" "c12" "c13" ; "c21" "c22" "c23"; "c31" "c32" "c33"];
abc = kron(a,b,c);
sys = qubit() ⊗ qubit() ⊗ osc(3);
abcv = subview(sys, abc);
"a21b12c31" == abcv[2,1,3,1,2,1] # True
```
"""
@inline function subview(s :: QSpace, op :: AbstractArray)
    # Compute correct permutation of dimensions
    # for the reshaped view of the array in memory
    dim = size(s)
    d = length(dim)
    rdim = reverse(dim)
    l = length(size(op))
    if l==1
        # Note: for a dimension list (2, 3, 4)
        #       the reshape should be: reshape(o, 4, 3, 2)
        #       with permutation: [3, 2, 1]
        PermutedDimsArray(reshape(op, rdim...), reverse(1:d))
    elseif l==2
        # Note: for a dimension list (2, 3, 4)
        #       the reshape should be: reshape(o, 4, 3, 2, 4, 3, 2)
        #       with permutation: [3, 2, 1, 6, 5 ,4]
        perm = vcat(reverse(1:d), reverse((d+1):(2*d)))
        PermutedDimsArray(reshape(op, rdim..., rdim...), perm)
    else
        error("subview assumes an initial 1D or 2D array")
    end
end
subview(s::QFactor, op::AbstractArray) = subview(QSpace(s), op)

"""
    superket(op::QOp)

Convert a square matrix operator into a ket vector in the superoperator space.
"""
superket(op::QOp) = reshape(op, length(op), 1)

"""
    unsuperket(op::AbstractArray)

Convert a superket back into a square matrix operator.
"""
function unsuperket(op::QOp)
    d = Int(sqrt(length(op)))
    d^2 == length(op) || error("Not a superket")
    reshape(op, d, d)
end

"""
    superopl(op::QOp)

Convert `op` into a superoperator that performs left multiplication
by `op`. `superopl(op)*superket(v) == superket(op * v)`
"""
superopl(op::QOp) = sparse(QComp(1.0)SparseArrays.I, size(op)) ⊗ op

"""
    superopr(op::QOp)

Convert `op` into a superoperator that performs right multiplication
by `op`. `superopr(op)*superket(v) == superket(v * op)`
"""
superopr(op::QOp) = transpose(op) ⊗ sparse(QComp(1.0)SparseArrays.I, size(op))

###################################################
# Workhorse functions
###

###
# Tensor Product
###
"""
    ⊗(q::AbstractArray, s::AbstractArray)
    ⊗(q::QObj, s::QObj)

Tensor product.

  - equivalent to `kron` for arrays
  - constructs tensor-product QSpaces from QFactors

"""
@inline ⊗(q::AbstractArray...) = kron(q...)
@inline function ⊗(q::QSpace...)
    factors = vcat(map(s->s.factors, q)...)
    QSpace(factors, QsOps())
end
@inline ⊗(q::QFactor...) = ⊗(map(QSpace,q)...)
⊗(q::QFactor, s::QSpace) = ⊗(QSpace(q),s)
⊗(s::QSpace, q::QFactor) = ⊗(s,QSpace(q))

###
# Lift to Joint Space
###

"""
    lift(q::QSpace, o::QOp, i::Int)

Lift an operator of a single factor into a joint space,
assuming that the factor is at position `i` of the tensor
product.

### Returns:
  - QOp : result of tensor product with appropriate identities
"""
@inline function lift(q::QSpace, o::QOp, i::Int)
    is = map(f -> f('i'), factors(q))
    splice!(is, i, [o])
    ⊗(is...)
end

###
# Partial Trace
###

"""
    ptr(s::QObj, o::QOp, n::Int...)

Partial trace of o, over subsystems in positions n, inside quantum space s.

### Returns:
  - (snew, onew)
  - snew : reduced system space
  - onew : reduced operator

"""
@inline ptr(s::QFactor, o::QOp, n::Int...) = tr(o)
@inline function ptr(s::QSpace, o::QOp, n::Int)
    v = subview(s,o)
    # Check that valid subsystem was specified
    dim = size(v)
    nsys = Int(fld(length(dim),2))
    (n > 0 && n <= nsys) || error("Partial trace over nonexistent system.")
    # Find second index matching n
    n2 = n+nsys
    # Helper function to truncate indices
    function trunc(l)
        l2 = [l...]
        deleteat!(l2, (n, n2))
        tuple(l2...)
    end

    # Remove indicated subsystem from system list
    sysind = [trunc(collect(1:length(dim)))[1:(nsys-1)]...]
    # Create new QSpace without indicated subsystem
    newsys = foldl(⊗, s.factors[sysind])
    if typeof(newsys)==QFactor
        newsys = QSpace(newsys)
    end
    # Create empty sparse matrix to hold result of trace
    m = spzeros(QComp, size(newsys(repeat("i",length(sysind))))...)
    # Create new view with easy indexing
    mv = subview(newsys, m)

    # Extract nonempty indices from original view
    vinds = CartesianIndices(v)
    oinds = LinearIndices(o)
    indices = []
    rows = rowvals(o)
    for c in 1:last(size(o))
        for r in nzrange(o, c)
            push!(indices, Tuple(vinds[oinds[rows[r], c]]))
        end
    end
    # Group indices that match when truncated
    for i in Set(map(trunc, indices)) # Set ensures no duplicates
        # Isolate diagonal elements within subgrouping, keep duplicates
        diagi = filter(e -> e[n]==e[n2], filter(j->trunc(j)==i, indices))
        if !isempty(diagi)
            # Set new sparse matrix elements as sums of diagonals
            mv[i...] = sum(map(j->v[j...], diagi))
        end
    end
    # Return both the new space, and the partially traced matrix
    newsys, m
end
@inline function ptr(s::QSpace, o::QOp, n::Int...)
    isempty(n) && return tr(o)
    n = sort([n...])
    s2, o2 = ptr(s,o,last(n))
    ptr(s2, o2, n[1:end-1]...)
end

#####################################################################
# Convenience functions for generating operators and states
###

###
# Sparse matrix basis generators by index
###

# Using (i) index for a 1D (vector) representation
Base.@propagate_inbounds Base.getindex(A::QFactor, i::Int) =
    sparsevec([i], [QComp(1)], length(A))
# Using (i,j) indices for the 2D representation
Base.@propagate_inbounds Base.getindex(A::QFactor, i::Int, j::Int) =
    sparse([i], [j], [QComp(1)], length(A), length(A))

# Using (i1,j1,k1,...,i2,j2,k2,...) indices for subsystem representation
@inline Base.@propagate_inbounds function Base.getindex(A::QSpace, inds::Int...)
    l = length(size(A))
    fs = factors(A)
    if length(inds) == l  # 1D (vector)
        ⊗([f[i] for (f,i) in zip(fs,inds)]...)
    elseif length(inds) == 2*l  # 2D (matrix)
        subinds = reshape(collect(inds), l, 2)
        ⊗([f[subinds[i,:]...] for (i,f) in enumerate(fs)]...)
    else
        error("Incorrect indexing.")
    end
end

# Overload function calls to match indexing from 0
@inline (q::QFactor)(i::Int...) = q[(1 .+ [i...])...]
@inline (q::QSpace)(i::Int...) = q[(1 .+ [i...])...]

###
# Ground states
###

"""
    groundvec(o::QObj)

Ground state vector of quantum space `o`.
"""
groundvec(o::QObj) = o[1]

"""
    groundvec(o::QObj)

Ground state density matrix of quantum space `o`.
"""
ground(o::QObj) = o[1,1]

"""
    projector(ψ::QKet)

Convert ket vector ψ into a projection operator.
"""
function projector(ψ::QKet)
    first(size(ψ)) == 1 && error("Not a column vector.")
    ρ = ψ * ψ'
    sparse(ρ / tr(ρ))
end
transition(ψl,ψr) = sparse(ψl * ψr') / (norm(ψl)*norm(ψr))


text:
module RomeoApp

using ArgParse
using MriResearchTools
using ROMEO

include("argparse.jl")
include("caller.jl")

function julia_main()::Cint
    try
        unwrapping_main(ARGS)
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

export unwrapping_main

end # module


text:
module WordleSolver
using Random, ProgressMeter, Counters, Random

export hist_check, wordle_solver, wordle_play, wordle_score

export fast_wordle_score

const five = 5

"""
    wordle_score(guess, answer)
Given two five letter words, return a five-tuple with the following meanings:
* `0` means that letter in `guess` does not appear in `answer`.
* `1` means that letter in `guess` appears in `answer`, but in another position.
* `2` means that latter in `guess` appears in `answer` in that position.

This assumes that `guess` and `answer` are valid five-letter words.
"""
function wordle_score(guess::String, code::String)::NTuple{five,Int}
    result = zeros(Int, five)
    used = zeros(Bool, five)

    # find letters that match in position
    for i = 1:five
        if guess[i] == code[i]
            result[i] = 2
            used[i] = true
        end
    end

    # find letters that are correct, but in wrong position
    for i = 1:five
        for j = 1:five
            if i == j
                continue
            end
            if guess[i] == code[j] && result[i] == 0 && !used[j]
                result[i] = 1
                used[j] = true
            end
        end
    end
    return Tuple(result)
end

include("startup.jl")

function fast_wordle_score(guess::Int, answer::Int)::Byte
    return SCORE_MATRIX[guess, answer]
end

include("evaluation.jl")

"""
    read_tuple
Read in characters from the keyboard and parse them into a 5-tuple of integers.
If the input isn't parsed as a 5-tuple of `Int`s, then throw an error.
Returns the 5-tuple.
"""
function read_tuple()::NTuple{five,Int}
    str = readline()
    parsed = Base.Meta.parse(str)
    result = eval(parsed)

    for t in result
        if t < 0 || t > 2
            error("Invalid tuple $result")
        end
    end

    return result
end

"""
    wordle_solver(scorer = min_max_score)
The computer guesses five-letter words and the human scores them.
Scores are entered as a list of five comma-separated numbers with
* 0 letter is not in the word
* 1 letter is in the word but in the wrong position
* 2 letter is in the word and in the right position

The optional parameter `scorer` is the name of a scoring function.
Choices are:
    * `min_max_score` --- default
    * `entropy_score`
    * `equi_score`
"""
function wordle_solver(scorer::Function = min_max_score)
    history = Dict{Int,Byte}()
    stop = code_2_byte((2, 2, 2, 2, 2))

    bad_message = "\nYour entry is invalid. Please try again."

    first_guess::Bool = true

    while true
        alist = filter_answers(history)
        if length(alist) == 0    # something wrong; no possible answer
            println("I give up")
            print("What was your word? ")
            w = readline()
            w = uppercase(w)
            try
                x = ANS_DICT[w]
            catch
                println("$w is not a valid word")
                return
            end
            history_check(history, w)
            return
        end

        shuffle!(alist)
        g = 1

        if first_guess
            g = FIRST_GUESS[scorer]
            first_guess = false
        else
            g = best_guess(alist, scorer)
        end
        guess = GUESS_LIST[g]

        # get user input
        result = (0, 0, 0, 0, 0)
        while true
            println("\nI guess the word is $guess")
            print("Enter the result --> ")
            try
                result = read_tuple()
            catch
                println(bad_message)
                continue
            end
            break
        end
        x = code_2_byte(result)

        history[g] = code_2_byte(result)
        if x == stop
            steps = length(history)
            println("\nSuccess in $steps guesses!!")
            return
        end
    end
end

"""
    history_check(history,answer::String)
Check the user inputs saved in the dictionary `d` with the code word `code`.
"""
function history_check(history::Dict{Int,Byte}, answer::String)
    answer = uppercase(answer)
    a = ANS_DICT[answer]  # code for answer
    for g in keys(history)
        x = fast_wordle_score(g, a)
        if x != history[g]
            guess = GUESS_LIST[g]
            println("\bWhen I guessed $guess you entered $(byte_2_code(history[g]))")
            println("     but you should have entered $(byte_2_code(x))")
        end
    end
end

"""
    read_word
Read a word from the keyboard. Return the pair `(w,valid)`
where `w` is the uppercase version of the input and `valid`
is `true` if `w` is a five-letter string of characters from
`A` to `Z`. We don't check if `w` is in a dictionary.
"""
function read_word()::Tuple{String,Bool}
    w = readline()
    valid = true
    w = uppercase(w)

    if w == "?"
        return w, true     # special case to escape
    end

    if length(w) != five
        valid = false
    end
    for c in w
        if c < 'A' || c > 'Z'
            valid = false
        end
    end
    return w, valid
end

"""
    wordle_play
The human tries to guess the secret word.
"""
function wordle_play(code::String)
    code = uppercase(code)
    count = 0
    while true
        count += 1
        w = "ABCDE"
        while true
            print("\nEnter your guess: ")
            w, valid = read_word()
            if valid
                break
            end
            println("Your input $w is not a valid $five-letter word. Please try again.")
        end

        if w == "?"
            println("The code word is $code")
            return
        end

        result = wordle_score(w, code)

        wordle_print(w, result)

        println("Score = $result")
        if w == code
            break
        end
    end
    println("\nSolved in $count steps")
    return
end

function wordle_play()
    idx = mod1(rand(Int), NA)
    code = ANS_LIST[idx]
    wordle_play(code)
end

function present_score(word::String, score::NTuple{five,Int})
    for c in word
        print("$c ")
    end
    println()

    symbols = "-+*"
    for t in score
        print(symbols[t+1], " ")
    end
    println()
end
include("wordle_print.jl")
include("auto_play.jl")

include("solve_tools.jl")

@info "Compilation completed"

end # module


text:
function jaccard_similarity(set1::Array, set2::Array)

    intersec = 0.0
    set1 = unique(set1)
    set2 = unique(set2)

    for i in set1
        for j in set2
            if i == j
                intersec += 1.0
            end
        end
    end
    return intersec / (length(set1) + length(set2) - intersec)
end


text:
"""
    inner_eltype(A::Type{<:AbstractArray{<:AbstractArray{V,M},N}})
    inner_eltype(A::AbstractArray{<:AbstractArray{V,M},N})

Returns `T`, the common `eltype` of the elements of `A`.
"""
inner_eltype(::Type{<:NestedArray{V}}) where {V} = @isdefined(V) ? V : Any
inner_eltype(A::NestedArray) = inner_eltype(typeof(A))

"""
    inner_ndims(A::Type{<:AbstractArray{<:AbstractArray{V,M},N}})
    inner_ndims(A::AbstractArray{<:AbstractArray{V,M},N})

Returns `M`, the common dimensionality of the elements of `A`.
"""
inner_ndims(::Type{<:NestedArray{<:Any,M}}) where {M} = M
inner_ndims(A::NestedArray) = inner_ndims(typeof(A))

"""
    inner_axes(A::AbstractArray{<:AbstractArray{V,M},N}, [dim])

Returns the common axes of the elements of `A`.
Throws an error if the elements of `A` do not have equal axes or if `A` is empty.
"""
@inline function inner_axes(A::NestedArray{<:Any,M}) where {M}
    # TODO would this be wrong for offset arrays?
    return isempty(A) ? ntuple(_ -> Base.OneTo(0), M) : _scan_inner(axes, A)
end

@inline function inner_axes(A::NestedArray{M}, d::Integer) where {M}
    d <= M ? inner_axes(A)[d] : Base.OneTo(0)
end

"""
    inner_size(A::AbstractArray{<:AbstractArray{V,M},N}, [dim])

Returns the common size of the elements of `A`.
Throws an error if the elements of `A` do not have equal axes or if `A` is empty.
"""
@inline inner_size(A::NestedArray) = map(Base.unsafe_length, inner_axes(A))
@inline inner_size(A::NestedArray, d::Integer) = Base.unsafe_length(inner_axes(A, d))

"""
    inner_length(A::AbstractArray{<:AbstractArray{V,M},N})

Returns the common length of the elements of `A`.
Throws an error if the elements of `A` do not have equal length or if `A` is empty.
"""
@inline inner_length(A::AbstractArray) = _scan_inner(length, A) # faster than prod(inner_size(A))

function _scan_inner(f::F, A::AbstractArray) where {F}
    if isempty(A)
        argerror("Cannot apply $f to an empty array")
    else
        x = f(first(A))
        for I in eachindex(A)[2:end]
            f(A[I]) == x || argerror("The element arrays of A do not have matching $f")
        end
        return x
    end
end


text:
@testset "Oblivious routing: uncertainty in the demand" begin
    opt = ECOS.Optimizer

    g = MetaDiGraph()
    add_vertex!(g, :name, "A")
    add_vertex!(g, :name, "B")
    add_vertex!(g, :name, "C")
    add_vertex!(g, :name, "D")
    add_vertex!(g, :name, "E")
    add_vertex!(g, :name, "F")
    add_edge!(g, 1, 3, :capacity, 2.0)
    add_edge!(g, 2, 3, :capacity, 2.0)
    add_edge!(g, 2, 6, :capacity, 2.0)
    add_edge!(g, 3, 4, :capacity, 20.0) # High capacity, shared.
    add_edge!(g, 4, 5, :capacity, 2.0)
    add_edge!(g, 1, 5, :capacity, 2.0) # Direct links.
    add_edge!(g, 4, 6, :capacity, 2.0)
    paths_d1 = [[Edge(1, 5)], [Edge(1, 3), Edge(3, 4), Edge(4, 5)]]
    paths_d2 = [[Edge(2, 6)], [Edge(2, 3), Edge(3, 4), Edge(4, 6)]]

    k = SimpleDiGraph(6)
    add_edge!(k, 1, 5)
    add_edge!(k, 2, 6)
    d1 = Edge(1, 5)
    d2 = Edge(2, 6)

    @testset "$type$(cg ? " with column generation" : "")" for (type, cg) in [(FlowFormulation, false), (PathFormulation, false), (PathFormulation, true)]
        @testset "$algo (with OPT_D computation: $compute_opt_d)" for (algo, compute_opt_d) in [(CuttingPlane, false), (CuttingPlane, true), (DualReformulation, false)]
            n_paths = (cg ? 1 : 5) # Only two paths for each demand, force column generation to generate a new one.
            mt = ModelType(Load(), MinimumMaximum(), type(), cg, algo(), ObliviousUncertainty(), UncertainDemand())
            rd = RoutingData(g, k, opt, mt, npaths=n_paths, model_exact_opt_d=compute_opt_d, logmessage=(x)->nothing)
            r = compute_routing(rd)

            @test r !== nothing
            @test r.data == rd
            @test r.n_iter >= 0
            if algo == CuttingPlane
                @test r.n_matrices >= 1
            else
                # Matrix output is disabled in this test.
                @test algo == DualReformulation
                @test r.n_matrices == 0
            end
            if cg
                @test r.n_columns >= 0
            else
                @test r.n_columns == 0
            end
            if algo == CuttingPlane
                @test r.n_cuts > 0
            else
                @test algo == DualReformulation
                @test r.n_cuts == 0
            end
            @test r.time_precompute_ms > 0.0
            @test r.total_time_solve_ms > 0.0
            @test r.total_time_export_ms >= 0.0
            @test r.total_time_ms > 0.0
            @test length(r.objectives) >= 1
            if algo == CuttingPlane
                @test length(r.matrices) >= 1
            else
                @test algo == DualReformulation
                # Matrix output is disabled in this test.
                if length(r.matrices) > 0
                    println("=== Check me: ")
                    println(r.matrices)
                end
                @test length(r.matrices) == 0
            end
            @test length(r.routings) >= 1
            @test r.master_model !== nothing

            @test r.objectives[r.n_iter] ≈ 1.0 atol=1.e-6
        end
    end
end


text:
__precompile__(true)
module IntervalMatrices

using LinearAlgebra

import Random
using Random: AbstractRNG, GLOBAL_RNG, seed!

using Reexport
@reexport using IntervalArithmetic

# ================================
# Type defining an interval matrix
# ================================
include("matrix.jl")
include("affine.jl")

# =================================
# Operations for interval matrices
# =================================
include("operations/arithmetic.jl")
include("operations/mult.jl")
include("operations/norm.jl")
include("operations/numops.jl")
include("operations/random.jl")
include("operations/setops.jl")

# ==================================================
# Methods to compute the power of an interval matrix
# ==================================================
include("operations/power.jl")

# =======================================================
# Methods to handle the exponential of an interval matrix
# =======================================================
include("exponential.jl")

# ==============================================================
# Methods to compute a correction matrix from an interval matrix
# ==============================================================
include("correction_matrices.jl")

# ========
# Exports
# ========
export AbstractIntervalMatrix,
       IntervalMatrix,
       set_multiplication_mode, get_multiplication_mode

export inf, sup,
       rand, sample,
       opnorm,
       diam_norm,
       scale, scale!,
       square,
       quadratic_expansion,
       exp_overapproximation,
       exp_underapproximation

export IntervalMatrixPower,
       increment!, increment,
       base,
       index

export AffineIntervalMatrix,
       AffineIntervalMatrix1

set_multiplication_mode(config[:multiplication])

end # module


text:
export get_sol2array, flatten, unique

"""Get the vector of vectors of results. `sol` can be `:EnsembleSolution` or other solutions from `DifferentialEquations.solve`
"""
get_sol2array(sol) = getfield.(sol.u, :u)

flatten(arr_arrs) = collect(Iterators.flatten(arr_arrs))

"""
Multiplication of series.
"""
function mul(array1D)
    m = 1
    for i in array1D
        m = m * i
    end
    return m
end

"""
Multipication of series with type converting.
"""
function mul(type, array1D)
    m = 1
    for i in array1D
        m = m * i
    end
    return convert(type, m)
end

function Base.unique(sols :: T; tol_digit=4) where T<: EnsembleSolution
    sols_ = map( u->trunc.(u, digits=tol_digit), sols.u  )
    unique!(sols_)
    return sols_
end


text:
using SparseTransforms
using Test

function test_signal()
    println("Testing signal retrieval")
    n, k = 8, 3
    locs = sample(0:2^n-1, 2^k, replace=false)
    strengths = Float64.(rand(Uniform(0.1, 10), 2^k)) .* (-1) .^ rand(Bool, 2^k)
    signal = TestSignal(n, locs, strengths)

    lazy_signal = LazySignal(n, locs, strengths)

    println("Testing lazy signal")
    for i = 1:2^n
        @test signal[i] ≈ lazy_signal[i]
    end

    println("Testing multiple index subsampling")
    for _ = 1:20
        idxs = sample(1:2^n-1, 2^k, replace=false)
        @test signal[idxs] ≈ lazy_signal[idxs]
    end
end


text:
using AWSCore
using XMLDict
using JSON

aws = aws_config()

struct Service
    request::Function
    service_name::String
    api_version::String
    service_endpoint::String
end

function rest_xml(service::Service, request_method::String, request_uri::String, args...)
    AWSCore.service_rest_xml(
        aws;
        service=service.service_name,
        version=service.api_version,
        verb=request_method,
        resource=request_uri,
        args=args
    )
end

function parse_metadata()
    return JSON.parsefile("metadata.json")
end

function service(service_name::String)
    service = metadata[service_name]

    request_function = getfield(Main, Symbol(service["protocol"]))
    api_version = service["api_version"]

    service_endpoint = ""

    if haskey(service, "endpoint")
        service_endpoint = service["endpoint"]
    end

    return Service(request_function, service_name, api_version, service_endpoint)
end

function main()
    global metadata = parse_metadata()
    s3 = service("s3")

    buckets = s3.request(s3, "GET", "/")["Buckets"].x
    buckets = xml_dict(parse_xml(string(buckets)))

    println(json(buckets, 4))
end

main()

max_stars_repo_name: mattBrzezinski/AWSCore-Prototype

max_issues_repo_name: mattBrzezinski/AWSCore-Prototype

max_forks_repo_name: mattBrzezinski/AWSCore-Prototype


text:
module Parameters

"""
    ThermodynamicsParameters

Parameters for Thermodynamics.jl.

# Example
```
import CLIMAParameters as CP
import Thermodynamics.Parameters as TP
FT = Float64;
toml_dict = CP.create_toml_dict(FT; dict_type = "alias");
aliases = string.(fieldnames(TP.ThermodynamicsParameters));
param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics");
param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...);
```
"""
Base.@kwdef struct ThermodynamicsParameters{FT}
    T_0::FT
    MSLP::FT
    cp_v::FT
    cp_l::FT
    cp_i::FT
    LH_v0::FT
    LH_s0::FT
    press_triple::FT
    T_triple::FT
    T_freeze::FT
    T_min::FT
    T_max::FT
    entropy_reference_temperature::FT
    entropy_dry_air::FT
    entropy_water_vapor::FT
    kappa_d::FT
    gas_constant::FT
    molmass_dryair::FT
    molmass_water::FT
    T_surf_ref::FT
    T_min_ref::FT
    grav::FT
    T_icenuc::FT
    pow_icenuc::FT
end

const ATP = ThermodynamicsParameters

Base.broadcastable(ps::ATP) = Ref(ps)

# wrappers
for fn in fieldnames(ATP)
    @eval $(fn)(ps::ATP) = ps.$(fn)
end

# Derived parameters
R_d(ps::ATP) = ps.gas_constant / ps.molmass_dryair
R_v(ps::ATP) = ps.gas_constant / ps.molmass_water
molmass_ratio(ps::ATP) = ps.molmass_dryair / ps.molmass_water
LH_f0(ps::ATP) = ps.LH_s0 - ps.LH_v0
e_int_v0(ps::ATP) = ps.LH_v0 - R_v(ps) * ps.T_0
e_int_i0(ps::ATP) = LH_f0(ps)
cp_d(ps::ATP) = R_d(ps) / ps.kappa_d
cv_d(ps::ATP) = cp_d(ps) - R_d(ps)
cv_v(ps::ATP) = ps.cp_v - R_v(ps)
cv_l(ps::ATP) = ps.cp_l
cv_i(ps::ATP) = ps.cp_i

end


text:
# this is an example of how to augment Wavs
using WAV
using SpeechAugment

##------------------------------------##
##---------- exmaple one -------------##
##------------------------------------##

# 1. read a wav file as a speech example
batchsize = 8;
data,fs = wavread("/home/xxx/ASpeechExample.wav");

# 2. init all the augmentation functions you want
echo  = initAddEcho(fs, (0.05,0.4), (3.0,3.2,2.5,3.5,2.0,3.0));
noise = initAddNoise("XXPathFullOfNoiseWAVs", 2, (5,15));
clip  = initClipWav((0.5,2.0));
drop  = initDropWav(fs, (0.09,0.15));
far   = initFarfieldWav(fs, (0.5,0.9));
speed = initSpeedWav((0.8,1.2));

# 3. make a function list
fnlist = [echo noise clip drop far speed];

# 4. augment one wav into #batchSize audios
wavlist = augmentWav(fnlist, data, batchsize)
for i = 1:8
    wavwrite(wavlist[i], "A$i.wav",Fs=16000,nbits=32)
end

##------------------------------------##
##---------- exmaple two -------------##
##------------------------------------##

# 1. prepare #batchSize audios
wavs = Vector(undef, batchsize);
for i = 1:batchsize
    wavs[i] = copy(data)
end
# 2. augment #batchSize audios
wavs = augmentWavs(fnlist, wavs)
for i = 1:batchsize
    wavwrite(wavs[i], "B$i.wav",Fs=16000,nbits=32)
end

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
# This file was generated, do not modify it. # hide
import MLJFlux
import MLJ
import DataFrames: DataFrame
import Statistics
import Flux
using Random
using PyPlot

MLJ.color_off() # hide
Random.seed!(11)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
const act_fn = relu

convx(ni, nf; ks=3, stride=1) = Conv(
        (ks, ks), ni => nf, stride = stride, pad = ks ÷ 2, init = Flux.kaiming_normal)

function convxlayer(ni, nf; ks=3, stride=1, zero_bn=false, act=true)
    bn = BatchNorm(nf, act ? act_fn : identity)
    fill!(bn.γ, zero_bn ? 0 : 1)
    return Chain(convx(ni, nf; ks = ks, stride = stride), bn)
end

struct ResBlock
    convs
    idconv
    pool
end

function ResBlock(expansion::Int, ni::Int, nh::Int; stride::Int = 1)
    nf, ni = nh * expansion, ni * expansion
    if expansion == 1
        layers = [
            convxlayer(ni, nh; stride = stride),
            convxlayer(nh, nf; zero_bn = true, act = false),
        ]
    else
        layers = [
            convxlayer(ni, nh; ks = 1),
            convxlayer(nh, nh; stride = stride),
            convxlayer(nh, nf; ks = 1, zero_bn = true, act = false)
        ]
    end

    return ResBlock(
        Chain(layers...),
        ni == nf ? identity : convxlayer(ni, nf; ks = 1, stride = 1),
        stride == 1 ? identity : MeanPool((2, 2))
    )
end
Flux.@functor ResBlock

(r::ResBlock)(x) = act_fn.(r.convs(x) .+ r.idconv(r.pool(x)))

function make_layer(expansion, ni, nf, n_blocks, stride)
    return Chain([
            ResBlock(expansion, i == 1 ? ni : nf, nf; stride =(i == 1 ? stride : 1))
            for i in 1:n_blocks
            ]...)
end

function XResNet(expansion, layers; c_in = 3)
    nfs = [c_in, (c_in+1)*8, 64, 64]
    stem = [convxlayer(nfs[i], nfs[i+1]; stride = i == 1 ? 2 : 1) for i in 1:3]

    nfs = [64 ÷ expansion, 64, 128, 256, 512]
    res_layers = [
        make_layer(expansion, nfs[i], nfs[i+1], l, i == 1 ? 1 : 2)
        for (i, l) in enumerate(layers)
    ]

    return Chain(
        stem...,
        MaxPool((3, 3); pad = 1, stride = 2),
        res_layers...,
    )
end

xresnet18(; kwargs...) = XResNet(1, [2, 2, 2, 2]; kwargs...)
xresnet50(; kwargs...) = XResNet(4, [3, 4, 6, 3]; kwargs...)


text:
using Xpress

mpb_path = ""
if VERSION < v"0.7"
    mpb_path = Pkg.dir("MathProgBase")
else
    import MathProgBase
    mpb_path = joinpath(dirname(pathof(MathProgBase)),"..")
end

include(joinpath(mpb_path,"test","linproginterface.jl"))
println("Testing linproginterface with solver Xpress.XpressSolver")
linprogsolvertest(Xpress.XpressSolver())
println("Done")

include(joinpath(mpb_path,"test","linprog.jl"))
linprogtest(Xpress.XpressSolver(PRESOLVE = 0))

include(joinpath(mpb_path,"test","mixintprog.jl"))
mixintprogtest(Xpress.XpressSolver())

include(joinpath(mpb_path,"test","quadprog.jl"))
quadprogtest(Xpress.XpressSolver(BARGAPSTOP = 1e-9))
socptest(Xpress.XpressSolver())
#qpdualtest(Xpress.XpressSolver()) #fails

include(joinpath(mpb_path,"test","conicinterface.jl"))
println("Testing coniclineartest with solver Xpress.XpressSolver")
coniclineartest(Xpress.XpressSolver())
println("Done")
println("Testing conicSOCtest with solver Xpress.XpressSolver")
conicSOCtest(Xpress.XpressSolver())
println("Done")


text:
centroid(V::Array{Float64,2}) = sum(V,2)/size(V,2)
centroid(V::Array{Int64,2}) = sum(V,2)/size(V,2)

	function larView(model::Tuple{Array{Float64, 2}, Array{Array{Int64, 1}, 1}})
		p["VIEW"](Plasm.lar2hpc(model...))
	end
	function larView(model::Tuple{Array{Array{Float64,N} where
	N,1},Array{Array{Int64,1},1}})
		V,CV = model
		V = hcat(V...)
		p["VIEW"](Plasm.lar2hpc(V,CV))
	end
	function larView(model::Array{Any,1})
		HPC_value_array = [Plasm.lar2hpc(item[1],item[2]) for item in model]
		p["VIEW"](p["STRUCT"](HPC_value_array))
	end
	function larView(V::Array{Float64, 2}, CV::Array{Array{Int64, 1}, 1})
		p["VIEW"](Plasm.lar2hpc(V, CV))
	end
	function larView(V::Array{Int64, 2}, CV::Array{Array{Int64, 1}, 1})
		W = convert(Array{Float64,2}, V)
		p["VIEW"](Plasm.lar2hpc(W, CV))
	end

	# LAR `model` ->  numbered `HPC` obyect
	function lar2numbered_hpc(larmodel,scaling=1.0)
		V,cells = larmodel
		VV,EV,FV,CV = cells

		Z = hcat(V[:,1],V)
		W = PyCall.PyObject([Any[Z[h,k] for h=1:size(Z,1)] for k=1:size(Z,2)])

		VV,EV,FV,CV = map(Plasm.doublefirst, [VV+1,EV+1,FV+1,CV+1])
		WW,EW,FW,CW = map(Plasm.array2list,[VV,EV,FV,CV])
		PyCall.PyObject([WW,EW,FW,CW])
		wire = p["MKPOL"](PyCall.PyObject([W,EW,[]]))

		VV,EV,FV,CV = VV-1,EV-1,FV-1,CV-1
		WW,EW,FW,CW = map(Plasm.array2list,[VV,EV,FV,CV])
		hpc = p.larModelNumbering(1,1,1)(W,PyCall.PyObject([WW,EW,FW,CW]),wire,scaling)
	end

	# Display a numbered `HPC` object from a `LAR` model with the `PyPlasm` viewer
	viewnumbered(larmodel,scaling=1.0) = p["VIEW"](lar2numbered_hpc(larmodel,scaling))

	# translate the columns of `V` matrix by sum with `t` vector
	function translate( t, V )
		broadcast(+,t,V)
	end

	# scale the columns of `V` matrix by product times `s` vector
	function scale( s, V )
		broadcast(*,s,V)
	end

	# rotate the columns of `V` matrix by properly using the `args` parameters
	function rotate(args,V)
		n = length(args)
		if n == 1 # rotation in 2D
			angle = args[1]; Cos = cos(angle); Sin = sin(angle)
			mat = eye(2)
			mat[1,1] = Cos;    mat[1,2] = -Sin;
			mat[2,1] = Sin;    mat[2,2] = Cos;
		elseif n == 3 # rotation in 3D
			mat = eye(3)
			angle = norm(collect(args)); axis = normalize(collect(args))
			Cos = cos(angle); Sin = sin(angle)
			# elementary rotations (in 3D)
			if axis[2]==axis[3]==0.0    # rotation about x
				mat[2,2] = Cos;    mat[2,3] = -Sin;
				mat[3,2] = Sin;    mat[3,3] = Cos;
			elseif axis[1]==axis[3]==0.0    # rotation about y
				mat[1,1] = Cos;    mat[1,3] = Sin;
				mat[3,1] = -Sin;    mat[3,3] = Cos;
			elseif axis[1]==axis[2]==0.0    # rotation about z
				mat[1,1] = Cos;    mat[1,2] = -Sin;
				mat[2,1] = Sin;    mat[2,2] = Cos;
			# general rotations (in 3D)
			else  # general 3D rotation (Rodrigues' rotation formula)
				I = eye(3) ; u = axis;
				Ux = Array([
					 0        -u[3]      u[2];
					 u[3]        0      -u[1];
					-u[2]      u[1]        0 ])
				UU = Array([
					 u[1]*u[1]    u[1]*u[2]    u[1]*u[3];
					 u[2]*u[1]    u[2]*u[2]    u[2]*u[3];
					 u[3]*u[1]    u[3]*u[2]    u[3]*u[3] ])
				mat = Cos*I + Sin*Ux + (1.0-Cos)*UU
			end
		end
		mat*V
	end

	# Visualize solid cells
	function viewsolidcells(sx=1.2, sy=1.2, sz=1.2)
		scaling = [sx; sy; sz]
		function viewsolidcells0(V,CV,FV,EV,cscCF,cscFE)
			local3cells = Lar.map_3cells_to_localbases(V,CV,FV,EV,cscCF,cscFE)
			hpcs = Any[]
			for local3cell in local3cells
				v,tv = local3cell
				centroid = sum(v,2)/size(v,2)
				scaledcentroid = scaling.*centroid
				translation = scaledcentroid - centroid
				w = v .+ translation
				hpc = p["SOLIDIFY"](Plasm.lar2hpc(w,tv))
				append!(hpcs, [hpc])
			end
			p["VIEW"](p["STRUCT"](hpcs))
			return viewsolidcells0
		end
	end

end


text:
"""
    ℞X(X)

Ensure X in the base unit U (assuming unit(X)=U).

If X has units of U or equivalent, or if it has no units,
it returns it in U. If X has any units other than K,
it then rasies an ArgumentError.

# Notes
It accepts vectors, but will throw an error if elements have different units.

# Examples
```julia-repl
julia> ℞T(300)
300 K
```
"""

using Unitful

"""
    ℞T(T)

Ensure T is in K.
"""
function ℞T(T)
    if unit.(upreferred.(T)) == unit.(ustrip.(T) .* u"K")
        return T .|> upreferred
    elseif unit.(upreferred.(T)) == unit.(ustrip.(T) .* Unitful.NoUnits)
        return T .* u"K"
    else
        throw(ArgumentError("Input: in K (or equivalent) or unitless"))
    end
end

"""
    ℞µ(µ)

Ensure µ is in kg/m/s.
"""
function ℞µ(µ)
    if unit.(upreferred.(µ)) == unit.(ustrip.(µ) .* u"kg" ./ u"m" ./ u"s")
        return µ .|> upreferred
    elseif unit.(upreferred.(µ)) == unit.(ustrip.(µ) .* Unitful.NoUnits)
        return µ .* u"kg" ./ u"m" ./ u"s"
    else
        throw(ArgumentError("Input: in kg/m/s (or equivalent) or unitless"))
    end
end

function ℞p(p)
    if unit.(upreferred.(p)) == unit.(ustrip.(p) .* u"kg/m/s^2")
        return p .|> upreferred
    elseif unit.(upreferred.(p)) == unit.(ustrip.(p) .* Unitful.NoUnits)
        return p .* u"kg/m/s^2"
    else
        throw(ArgumentError("Input: in kg/m/s^2 (or equivalent) or unitless"))
    end
end

function ℞M_wt(M_wt)
    if unit.(upreferred.(M_wt)) == unit.(ustrip.(M_wt) .* u"kg/mol")
        return M_wt .|> upreferred
    elseif unit.(upreferred.(M_wt)) == unit.(ustrip.(M_wt) .* Unitful.NoUnits)
        return M_wt .* u"kg/mol"
    else
        throw(ArgumentError("Input: in kg/mol (or equivalent) or unitless"))
    end
end


text:
using BinaryBuilder

# zlib version
version = v"1.2.11"

# Collection of sources required to build zlib
sources = [
    "https://zlib.net/zlib-$(version).tar.gz" =>
    "c3e5e9fdd5004dcb542feda5ee4f0ff0744628baf8ed2dd5d66f8ca1197cb1a1",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/zlib-*

# On windows platforms, our ./configure and make invocations differ a bit
if [[ ${target} == *-w64-mingw* ]]; then
    EXTRA_CONFIGURE_FLAGS="--sharedlibdir=${prefix}/bin"
    EXTRA_MAKE_FLAGS="SHAREDLIB=libz.dll SHAREDLIBM=libz-1.dll SHAREDLIBV=libz-1.2.11.dll LDSHAREDLIBC= "
fi

./configure ${EXTRA_CONFIGURE_FLAGS} --prefix=${prefix}
make install ${EXTRA_MAKE_FLAGS} -j${nproc}

# Materialize DLL files (https://github.com/bicycle1885/CodecZlib.jl/issues/24).
if [[ ${target} == *-w64-mingw* ]]; then
    mkdir ${WORKSPACE}/destdir/tmp
    cp -L ${WORKSPACE}/destdir/bin/* ${WORKSPACE}/destdir/tmp
    rm -r ${WORKSPACE}/destdir/bin
    mv ${WORKSPACE}/destdir/tmp ${WORKSPACE}/destdir/bin
fi
"""

# Build for ALL THE PLATFORMS!
platforms = supported_platforms()
#platforms = expand_gcc_versions(platforms)

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "libz", :libz),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

build_tarballs(ARGS, "Zlib", version, sources, script, platforms, products, dependencies)


text:
# use of Symbols as node is deprecated
src_dst_symbol_as_node = (:DirectedEdge, :UndirectedEdge)
for f in src_dst_symbol_as_node
    @eval @deprecate $f(src::Symbol, dst::Symbol) $f(Node(src), Node(dst))
end

dst_symbol_as_node = (:*, :^, src_dst_symbol_as_node...)
for f in dst_symbol_as_node
    @eval @deprecate $f(src, dst::Symbol) $f(src, Node(dst))
end

src_symbol_as_node = (:ModifyingNode, dst_symbol_as_node...)
for f in src_symbol_as_node
    @eval @deprecate $f(src::Symbol, dst) $f(Node(src), dst)
end

@deprecate keep(x::Symbol, y, z) keep(Node(x), y, z)
@deprecate keep(x::Symbol, y) keep(Node(x), y)


text:
# This package is based on the shrinkage estimator from
# https://cran.r-project.org/web/packages/entropy/

module ShrinkageEstimators

include("ExportedFunctions.jl")
include("Discretization.jl")
include("Shrinkage.jl")
include("Formulae.jl")

end # module

max_stars_repo_name: Tchanders/ShrinkageEstimators

max_issues_repo_name: Tchanders/ShrinkageEstimators

max_forks_repo_name: Tchanders/ShrinkageEstimators


text:
σ = [2,2,1]
phase_unwrap = :laplacian
phase_scaling_type = :linear
phase_scaling_strength = 4
tmp_folder = "tmp"
mag_combine = (:CNR=>(:gm, :wm))
o = Options(mag_combine=mag_combine, phase_hp_σ=σ, phase_unwrap=phase_unwrap, phase_scaling_type=phase_scaling_type, phase_scaling_strength=phase_scaling_strength, writesteps=tmp_folder)
@test o.phase_hp_σ == σ
@test o.phase_unwrap == phase_unwrap
@test o.phase_scaling_type == phase_scaling_type
@test o.phase_scaling_strength == phase_scaling_strength
@test o.mag_combine == (:CNR=>(:gm, :wm))
@test o.mag_sens === nothing
@test o.writesteps == "tmp"
@test o.mag_softplus == true

mkpath(tmp_folder)
saveconfiguration(o)
text = """mag_combine: :CNR => (:gm, :wm)
mag_sens: nothing
mag_softplus: true
phase_unwrap: laplacian
phase_hp_σ: [2, 2, 1]
phase_scaling_type: linear
phase_scaling_strength: 4
writesteps: tmp
"""
file = "$tmp_folder/settings_swi.txt"
@test read(file, String) == text
rm(file)
try rm(tmp_folder) catch end

@test isdir(CLEARSWI.dir("test", "testData", "small"))

printstyled("Tests Passed!"; color=:green)


text:
@testset "non_convex" begin

    @testset "non_convex_small_sigma" begin

        # Setup problem
        P = sparse([2. 5.; 5. 1.])
        q = [3.; 4]
        A = sparse([-1 0; 0 -1; -1 -3; 2 5; 3 4])
        u = [0.; 0.; -15; 100; 80]
        l = -Inf * ones(length(u))
        options = Dict(:verbose => false, :sigma => 1e-06)
        model = OSQP.Model()
        try
            # Setup should fail due to (P + sigma I) having a negative eigenvalue
            global test_setup = 1
            OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)
        catch
            global test_setup = 0
        end

        @test test_setup == 0

    end

    @testset "non_convex_big_sigma" begin

        # Setup workspace with new sigma
        P = sparse([2. 5.; 5. 1.])
        q = [3.; 4]
        A = sparse([-1 0; 0 -1; -1 -3; 2 5; 3 4])
        u = [0.; 0.; -15; 100; 80]
        l = -Inf * ones(length(u))
        options = Dict(:verbose => false, :sigma => 5.)
        model = OSQP.Model()
        OSQP.setup!(model; P=P, q=q, A=A, l=l, u=u, options...)

        # Solve problem
        results = OSQP.solve!(model)

        @test isnan(results.info.obj_val)
        @test results.info.status == :Non_convex

    end

end


text:
struct AlwaysAssertionError
    msg::String
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

function randombubbles()
	store = []
	mycircle(r,n) = Lar.circle(r)([n])
	for k=1:50
		r = rand()
		n = 8
		while true
			n = abs(rand(Int8)+1)
			if n>2 break end
		end
		scale = Lar.s(0.25,0.25)
		transl = Lar.t(rand(2)...)
		s = Lar.Struct([ transl, scale, mycircle(r,n) ])
		push!(store, Lar.struct2lar(s))
	end
	s = Lar.Struct(store)
	V,EV = Lar.struct2lar(s)
	V = GL.normalize2(V)
	GL.VIEW([ GL.GLLines(V,EV) ])

	W,EW = Lar.fragmentlines((V,EV))
	V,FVs,EVs = Lar.arrange2D(W,EW)
end

# ////////////////////////////////////////////////////////////
# generation of 2D arrangement
V,FVs,EVs = randombubbles()

# native OpenGL visualization
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,1,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,3,1));
GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,1));

GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,1,1));
GL.VIEW(GL.GLExplode(V,EVs,1,1,1,1,1));


text:
"""
        CalEqOddsWrapper

Calibrated equalized odds is a post-processing technique that optimizes over calibrated classifier score outputs to find probabilities with which to change output labels with an equalized odds objective
"""

struct CalEqOddsWrapper{M<:MLJBase.Model} <: Deterministic
    grp::Symbol
    classifier::M
    fp_rate::Int64
    fn_rate::Int64
    alpha::Float64
end

"""
        CalEqOddsWrapper(classifier=nothing, grp=:class, fp_rate=1, fn_rate=1)

Instantiates CalEqOddsWrapper which wraps the classifier
"""
function CalEqOddsWrapper(; classifier::MLJBase.Model=nothing,grp::Symbol=:class, fp_rate=1, fn_rate=1, alpha=1.0)
    model =  CalEqOddsWrapper(grp, classifier, fp_rate, fn_rate, alpha)
    message = MLJBase.clean!(model)
    isempty(message) || @warn message
    return model
end

function MLJBase.clean!(model::CalEqOddsWrapper)
    warning = ""
	model.classifier!=nothing || (warning *= "No classifier specified in model\n")
    target_scitype(model) <: AbstractVector{<:Finite} || (warning *= "Only Binary Classifiers are supported\n")
	(model.alpha>=0 && model.alpha<=1) || (warning*="alpha should be between 0 and 1 (inclusive)\n")
    return warning
end

# returns mixrates
function MMI.fit(model::CalEqOddsWrapper, verbosity::Int,
	X, y)
	grps = X[:, model.grp]
	n = length(levels(grps)) # Number of different values for sensitive attribute
	# As equalized odds is a postprocessing algorithm, the model needs to be fitted a_
	mch = machine(model.classifier, X, y)
	fit!(mch)
	ŷ = MMI.predict(mch, X)
	if typeof(ŷ[1]) <: MLJBase.UnivariateFinite
		ŷ = MLJBase.mode.(ŷ)
	end

	labels = levels(y)
	favLabel = labels[2]
	unfavLabel = labels[1]
	y = y.==favLabel
	ŷ = ŷ.==favLabel

    # Compares the first class against all others
	a_Class = levels(grps)[1]
	a_Grp = grps .== a_Class

	a_const = ŷ[a_Grp]

	b_Grp = grps .!= a_Class

	b_const = ŷ[b_Grp]

	br1 = mean(y[a_Grp])
	br2 = mean(y[b_Grp])
	fp_rate = model.fp_rate
	fn_rate = model.fn_rate
	mix_rate1,mix_rate2 = CalEqOdds(ŷ[a_Grp],ŷ[b_Grp],y[a_Grp],y[b_Grp],fp_rate,fn_rate)

	# ft = fair_tensor(categorical(ŷ), categorical(y), categorical(grps))
	fitresult = [mix_rate1, mix_rate2, br1, br2, mch.fitresult, labels]
	return fitresult, nothing, nothing
end

# Modifies the predictions of the model based on the mixrates
function MMI.predict(model::CalEqOddsWrapper, fitresult, Xnew)
	mix_rate1, mix_rate2, y1br, y2br,classifier_fitresult, labels = fitresult
	ŷ = MMI.predict(model.classifier, classifier_fitresult, Xnew)
	if typeof(ŷ[1]) <: MLJBase.UnivariateFinite
		ŷ = MLJBase.mode.(ŷ)
	end
	favLabel = labels[2]
	unfavLabel = labels[1]

	grps = Xnew[:, model.grp]
	a_Class = levels(grps)[1]
	a_Grp = grps .== a_Class
	b_Grp = grps .!= a_Class
	a_const = shuffle(findall((grps.== a_Class)))
	b_const = shuffle(findall((grps.!= a_Class)))
	indices1 = a_const[1:convert(Int64, floor((mix_rate1)*(length(a_const))))]
	indices2 =  b_const[1:convert(Int64, floor((mix_rate2)*(length(b_const))))]
	new_ŷ = deepcopy(ŷ)
	new_ŷ[indices1] .= convert(Int64,round(y1br))
	new_ŷ[indices2] .= convert(Int64,round(y2br))
	return new_ŷ
end

MMI.input_scitype(::Type{<:CalEqOddsWrapper{M}}) where M = input_scitype(M)
MMI.target_scitype(::Type{<:CalEqOddsWrapper{M}}) where M = AbstractVector{<:Finite{2}}

# Corresponds to calib_eq_odds function which calculates mixrates
function CalEqOdds(ŷ1, ŷ2, y1, y2, fp_rate, fn_rate)
	y1_triv = ones(length(y1))*mean(y1)
	y2_triv = ones(length(y2))*mean(y2)
	if fn_rate == 0
		cost1 = mean(ŷ1[y1.==0])
		cost2 = mean(ŷ2[y2.==0])
		trivial_cost1 = mean(y1_triv[y1.==0])
		trivial_cost2 = mean(y2_triv[y2.==0])
	elseif fp_rate == 0
		cost1 = 1-mean(ŷ1[y1.==1])
		cost2 = 1-mean(ŷ2[y2.==1])
		trivial_cost1 = 1-mean(y1_triv[y1.==1])
		trivial_cost2 = 1-mean(y2_triv[y2.==1])
	else
		norm_const = fp_rate+fn_rate
		cost1 = fp_rate/norm_cost * mean(ŷ1[y1.==0]) + fn_rate/norm_cost * (1-mean(ŷ1[y1.==1]))
		cost2 = fp_rate/norm_cost * mean(ŷ2[y2.==0]) + fn_rate/norm_cost * (1-mean(ŷ2[y2.==1]))
		trivial_cost1 = fp_rate/norm_cost * mean(y1_triv[y1.==0]) + fn_rate/norm_cost * (1-mean(y1_triv[y1.==1]))
		trivial_cost2 = fp_rate/norm_cost * mean(y2_triv[y2.==0]) + fn_rate/norm_cost * (1-mean(y2_triv[y2.==1]))
	end

        if cost2>cost1 && cost2<trivial_cost1
	    	mix_rate1 = (cost2-cost1)/(trivial_cost1-cost1)
			mix_rate2 = 0
		elseif cost1>cost2 && cost1<trivial_cost2
			mix_rate1 = 0
			mix_rate2 = (cost1-cost2)/(trivial_cost2-cost2)
		else
			mix_rate1 = 0
			mix_rate2 = 0
		end
	return mix_rate1, mix_rate2
end


text:
# Newton with line search
function newton!(mechanism::Mechanism{T,Nn,Nb,Ne,0};
        ε = 1e-10, newtonIter = 100, lineIter = 10, warning::Bool = false
    ) where {T,Nn,Nb,Ne}

    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    graph = mechanism.graph
    ldu = mechanism.ldu

    normf0 = normf(mechanism)
    for n = Base.OneTo(newtonIter)
        setentries!(mechanism)
        factor!(graph, ldu)
        solve!(mechanism) # x̂1 for each body and constraint

        normf1 = lineSearch!(mechanism, normf0;iter = lineIter, warning = warning)
        normΔs1 = normΔs(mechanism)

        foreachactive(updatesolution!, bodies)
        foreachactive(updatesolution!, eqcs)

        if normf1 < ε && normΔs1 < ε
            # warning && (@info string("Newton iterations: ",n))
            return
        else
            normf0 = normf1
        end
    end

    warning && (@info string("newton! did not converge. n = ", newtonIter, ", tol = ", normf(mechanism), "."))
    return
end

# Newton on interior point with line search
function newton!(mechanism::Mechanism{T,Nn,Nb,Ne,Ni};
        ε = 1e-10, σ = 0.1, μ = 1.0, newtonIter = 100, lineIter = 10, warning::Bool = false
    ) where {T,Nn,Nb,Ne,Ni}

    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    ineqcs = mechanism.ineqconstraints
    graph = mechanism.graph
    ldu = mechanism.ldu

    foreachactive(resetVars!, ineqcs)
    mechanism.μ = μ

    meritf0 = meritf(mechanism)
    for n = Base.OneTo(newtonIter)
        setentries!(mechanism) #wo mein Matrix steht da, da ist auch extendDanddeltas (rechte Seite)
        factor!(graph, ldu)
        solve!(mechanism) # x̂1 for each body and constraint

        meritf1 = lineSearch!(mechanism, meritf0;iter = lineIter, warning = warning)

        foreachactive(updatesolution!, bodies)
        foreachactive(updatesolution!, eqcs)
        foreachactive(updatesolution!, ineqcs)

        if normf(mechanism) < ε
            warning && (@info string("Newton iterations: ", n))
            return
        else
            while meritf1 < mechanism.μ
                mechanism.μ = σ * mechanism.μ
                meritf1 = meritf(mechanism)
            end
            meritf0 = meritf1
        end
    end

    warning && (@info string("newton! did not converge. n = ", newtonIter, ", tol = ", normf(mechanism), "."))
    return
end

function newton!(mechanism::LinearMechanism{T,Nn,Nb,Ne,0};
        ε = 1e-10, newtonIter = 100, lineIter = 10, warning::Bool = false
    ) where {T,Nn,Nb,Ne}

    A = mechanism.A
    Bu = mechanism.Bu
    Bλ = mechanism.Bλ
    G = mechanism.G
    nc = size(G)[1]
    Z = zeros(T,nc,nc)

    normf0 = normf(mechanism)
    for n = Base.OneTo(newtonIter)
        M = [-I Bλ;G Z]
        Δsol = M\[fdynamics(mechanism); fconstraints(mechanism)]
        mechanism.Δz = Δsol[1:12*Nb]
        mechanism.Δλ = Δsol[12*Nb+1:end]

        normf1 = lineSearch!(mechanism, normf0;iter = lineIter, warning = warning)
        # normΔs1 = normΔs(mechanism)

        mechanism.zsol[1] = mechanism.zsol[2]
        mechanism.λsol[1] = mechanism.λsol[2]

        if normf1 < ε #&& normΔs1 < ε
            # warning && (@info string("Newton iterations: ",n))
            return
        else
            normf0 = normf1
        end
    end

    warning && (@info string("newton! did not converge. n = ", newtonIter, ", tol = ", normf(mechanism), "."))
    return
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# src/quadexpr.jl
# Defines all types relating to expressions with a quadratic and affine part
# - GenericQuadExpr             ∑qᵢⱼ xᵢⱼ  +  ∑ aᵢ xᵢ  +  c
#   - QuadExpr                  Alias for (Float64, Variable)
# - GenericQuadConstraint       ∑qᵢⱼ xᵢⱼ  +  ∑ aᵢ xᵢ  +  c  [≤,≥]  0
#   - QuadConstraint            Alias for (Float64, Variable)
# Operator overloads in src/operators.jl
#############################################################################

#############################################################################
# GenericQuadExpr
# ∑qᵢⱼ xᵢⱼ  +  ∑ aᵢ xᵢ  +  c
type GenericQuadExpr{CoefType,VarType} <: AbstractJuMPScalar
    qvars1::Vector{VarType}
    qvars2::Vector{VarType}
    qcoeffs::Vector{CoefType}
    aff::GenericAffExpr{CoefType,VarType}
end
coeftype{C,V}(::GenericQuadExpr{C,V}) = C

Base.isempty(q::GenericQuadExpr) = (length(q.qvars1) == 0 && isempty(q.aff))
Base.zero{C,V}(::Type{GenericQuadExpr{C,V}}) = GenericQuadExpr(V[], V[], C[], zero(GenericAffExpr{C,V}))
Base.one{C,V}(::Type{GenericQuadExpr{C,V}})  = GenericQuadExpr(V[], V[], C[],  one(GenericAffExpr{C,V}))
Base.zero(q::GenericQuadExpr) = zero(typeof(q))
Base.one(q::GenericQuadExpr)  =  one(typeof(q))
Base.copy(q::GenericQuadExpr) = GenericQuadExpr(copy(q.qvars1),copy(q.qvars2),copy(q.qcoeffs),copy(q.aff))

function Base.append!{T,S}(q::GenericQuadExpr{T,S}, other::GenericQuadExpr{T,S})
    append!(q.qvars1, other.qvars1)
    append!(q.qvars2, other.qvars2)
    append!(q.qcoeffs, other.qcoeffs)
    append!(q.aff, other.aff)
    q
end

function assert_isfinite(q::GenericQuadExpr)
    assert_isfinite(q.aff)
    for i in 1:length(q.qcoeffs)
        isfinite(q.qcoeffs[i]) || error("Invalid coefficient $(q.qcoeffs[i]) on quadratic term $(q.qvars1[i])*$(q.qvars2[i])")
    end
end

function Base.isequal{T,S}(q::GenericQuadExpr{T,S},other::GenericQuadExpr{T,S})
    isequal(q.aff,other.aff)   || return false
    length(q.qvars1) == length(other.qvars1) || return false
    for i in 1:length(q.qvars1)
        isequal(q.qvars1[i],  other.qvars1[i]) || return false
        isequal(q.qvars2[i],  other.qvars2[i]) || return false
        isequal(q.qcoeffs[i],other.qcoeffs[i]) || return false
    end
    return true
end

# Alias for (Float64, Variable)
typealias QuadExpr GenericQuadExpr{Float64,Variable}
Base.convert(::Type{QuadExpr}, v::Union{Real,Variable,AffExpr}) = QuadExpr(Variable[], Variable[], Float64[], AffExpr(v))
QuadExpr() = zero(QuadExpr)

function setobjective(m::Model, sense::Symbol, q::QuadExpr)
    m.obj = q
    if m.internalModelLoaded
        if method_exists(MathProgBase.setquadobjterms!, (typeof(m.internalModel), Vector{Cint}, Vector{Cint}, Vector{Float64}))
            verify_ownership(m, m.obj.qvars1)
            verify_ownership(m, m.obj.qvars2)
            MathProgBase.setquadobjterms!(m.internalModel, Cint[v.col for v in m.obj.qvars1], Cint[v.col for v in m.obj.qvars2], m.obj.qcoeffs)
        else
            isa(m.internalModel, MathProgBase.AbstractLinearQuadraticModel) && Base.warn_once("Solver does not support adding a quadratic objective to an existing model. Hot-start is disabled.")
            m.internalModelLoaded = false
        end
    end
    setobjectivesense(m, sense)
end

# Copy a quadratic expression to a new model by converting all the
# variables to the new model's variables
function Base.copy(q::QuadExpr, new_model::Model)
    QuadExpr(copy(q.qvars1, new_model), copy(q.qvars2, new_model),
                copy(q.qcoeffs), copy(q.aff, new_model))
end

function getvalue(a::QuadExpr)
    ret = getvalue(a.aff)
    for it in 1:length(a.qvars1)
        ret += a.qcoeffs[it] * getvalue(a.qvars1[it]) * getvalue(a.qvars2[it])
    end
    return ret
end
getvalue(arr::Array{QuadExpr}) = map(getvalue, arr)

##########################################################################
# GenericQuadConstraint
# ∑qᵢⱼ xᵢⱼ  +  ∑ aᵢ xᵢ  +  c  [≤,≥]  0
# As RHS is implicitly taken to be zero, we store only LHS and sense
type GenericQuadConstraint{QuadType} <: AbstractConstraint
    terms::QuadType
    sense::Symbol
end
Base.copy{CON<:GenericQuadConstraint}(c::CON, new_model::Model) = CON(copy(c.terms, new_model), c.sense)

# Alias for (Float64, Variable)
typealias QuadConstraint GenericQuadConstraint{QuadExpr}

function Base.copy(c::QuadConstraint, new_model::Model)
    return QuadConstraint(copy(c.terms, new_model), c.sense)
end

function addconstraint(m::Model, c::QuadConstraint)
    push!(m.quadconstr,c)
    if m.internalModelLoaded
        if method_exists(MathProgBase.addquadconstr!, (typeof(m.internalModel),
                                                       Vector{Cint},
                                                       Vector{Float64},
                                                       Vector{Cint},
                                                       Vector{Cint},
                                                       Vector{Float64},
                                                       Char,
                                                       Float64))
            if !((s = string(c.sense)[1]) in ['<', '>', '='])
                error("Invalid sense for quadratic constraint")
            end
            terms = c.terms
            verify_ownership(m, terms.qvars1)
            verify_ownership(m, terms.qvars2)
            MathProgBase.addquadconstr!(m.internalModel,
                                        Cint[v.col for v in c.terms.aff.vars],
                                        c.terms.aff.coeffs,
                                        Cint[v.col for v in c.terms.qvars1],
                                        Cint[v.col for v in c.terms.qvars2],
                                        c.terms.qcoeffs,
                                        s,
                                        -c.terms.aff.constant)
        else
            Base.warn_once("Solver does not appear to support adding quadratic constraints to an existing model. Hot-start is disabled.")
            m.internalModelLoaded = false
        end
    end
    return ConstraintRef{Model,QuadConstraint}(m,length(m.quadconstr))
end
addconstraint(m::Model, c::Array{QuadConstraint}) =
    error("Vectorized constraint added without elementwise comparisons. Try using one of (.<=,.>=,.==).")

function addVectorizedConstraint(m::Model, v::Array{QuadConstraint})
    ret = Array(ConstraintRef{Model,QuadConstraint}, size(v))
    for I in eachindex(v)
        ret[I] = addconstraint(m, v[I])
    end
    ret
end


text:
module ABBAj
include("Compress.jl")
include("Digitization.jl")
include("Reconstruction.jl")
using .Compress, .Digitization, .Reconstruction
using Statistics

export fit_transform, inverse_transform, standard_data, load_sample

"""
ABBA: adaptive Brownian bridge-based symbolic aggregation of time series
    A novel time series symbolic compressed technique with mean-based clustering.

Examples:
julia> time_series = load_samples(); # load time series samples

julia> symbols, model = fit_transform(time_series, 20, 0.1); # use 20 symbols with compressed tolerance of 0.1
julia> r_time_series = inverse_transform(symbols, model, time_series[1]); # inverse transform time series

julia> symbols, model = fit_transform(time_series, 20, 0.1, 10); # use 20 symbols with compressed tolerance of 0.1
                                                                 # and run in parallel kmeans with 10 threads
julia> r_time_series = inverse_transform(symbols, model, time_series[1]); # inverse transform time series

Reference
Elsworth, S., Güttel, S. ABBA: adaptive Brownian bridge-based symbolic aggregation of time series.
Data Min Knowl Disc 34, 1175–1200 (2020). https://doi.org/10.1007/s10618-020-00689-6

"""
global ABBA_PARAM_MU, ABBA_PARAM_SCL = Ref{Float64}(0.0), Ref{Float64}(1.0)

function fit_transform(series::AbstractVector, k::Int, tol=0.5::AbstractFloat,
                       nthreads=Threads.nthreads()::Int, max_len=typemax(Int)::Int
                       )
    @assert(k > 0 && typeof(k) == Int)
    @assert(tol > 0)
    @assert(max_len >=1)
    global ABBA_PARAM_MU, ABBA_PARAM_SCL;
    scl_series, ABBA_PARAM_MU, ABBA_PARAM_SCL= standard_data(series);
    if (ABBA_PARAM_SCL==0)
       ABBA_PARAM_SCL = 1;
    end
    # scl_series = copy(series);
    compressed = compress(scl_series, tol, max_len);
    symbols, model = digitization_k(compressed, k, nthreads);

    return symbols, model;
end

function inverse_transform(symbols::Array{Char,1},
                           model::Digitization.Model,
                           start=0.0::AbstractFloat
                           )
    pieces = inverse_digitize(symbols, model);
    pieces = quantize(pieces);
    reconstructed_ts = inverse_compress(pieces, start);
    reconstructed_ts = broadcast(*, reconstructed_ts, ABBA_PARAM_SCL);
    reconstructed_ts = broadcast(+, reconstructed_ts, ABBA_PARAM_MU);
    return reconstructed_ts;
end

function standard_data(series::AbstractVector)
    PARARM_SCL = Statistics.std(series);
    PARARM_MU =  mean(series);
    scl_series = broadcast(-, copy(series), PARARM_MU);
    scl_series = broadcast(/, scl_series, PARARM_SCL);
    return scl_series, PARARM_MU, PARARM_SCL;
end

function load_sample(len=1000::Int, freq=20.0::AbstractFloat)
    sample = zeros(len);
    j = 1;
    for i in range(0, len-1, step=1)*(1/freq)
        sample[j] = sin(i);
        j = j + 1;
    end
    return sample;
end

end

max_stars_repo_name: nla-group/Julia-ABBA

max_issues_repo_name: nla-group/Julia-ABBA

max_forks_repo_name: nla-group/Julia-ABBA


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module ImageQuilting

using ImageFiltering
using ImageMorphology
using LightGraphs
using LightGraphsFlows
using Base: @nexprs, @nloops, @nref
using Primes: factor
using StatsBase: sample, weights
using ProgressMeter: Progress, next!
using FFTW: set_num_threads
using CpuId: cpucores
using RecipesBase
using SparseArrays: spzeros
using Random: shuffle!, randperm
using Statistics: mean, std

include("utils.jl")
include("plot_recipes.jl")
include("imfilter_cpu.jl")
include("relaxation.jl")
include("tau_model.jl")
include("graphcuts.jl")
include("voxel_reuse.jl")
include("iqsim.jl")

# optionally load GeoStats.jl API
using Requires
function __init__()
    @require GeoStatsBase="323cb8eb-fbf6-51c0-afd0-f8fba70507b2" include("geostats.jl")
end

export
  # functions
  iqsim,
  voxelreuse,

  # geostats solver
  ImgQuilt

end


text:
using SolidStateDetectors
using IntervalSets
using Plots
T = Float32

struct CustomChargeDensity{T, N} <: SolidStateDetectors.AbstractChargeDensity{T}
    G1::SolidStateDetectors.AbstractGeometry{T, N}
    G2::SolidStateDetectors.AbstractGeometry{T, N}
    L1::SolidStateDetectors.AbstractChargeDensity{T}
    L2::SolidStateDetectors.AbstractChargeDensity{T}
end

#compact version of the function I showed before
function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::SolidStateDetectors.AbstractCoordinatePoint{T})::T where {T}
       (pt in cdm.G1) * SolidStateDetectors.get_charge_density(cdm.L1, pt) +
       (pt in cdm.G2) * SolidStateDetectors.get_charge_density(cdm.L2, pt)
end

G1 = SolidStateDetectors.Tube{T}(0..0.1,0..2π,0..0.04,missing)     #Tube with 0 ≤ z ≤ 0.04
G2 = SolidStateDetectors.Box{T}((-1,1),(-1,1),(0.02,0.08),missing) #Box with 0.02 ≤ z ≤ 0.08
CD1 = SolidStateDetectors.CylindricalChargeDensity{T}((0.0f0, 0.0f0, 1.0f16), (1.0f17, 0.0f0, 1.0f17))
CD2 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0, 0.0f0, 1.0f15), (1.0f17, 0.0f0, 0.0f0))
ccdm = CustomChargeDensity(G1,G2,CD1,CD2)

simulation = Simulation{T}(SSD_examples[:InvertedCoax]) #adjusted phi in the grid to go from 0 to 360 to avoid 2D calculation
simulation.detector.semiconductors[1].charge_density_model = ccdm #apply CustomChargeDensity to the semiconductor
calculate_electric_potential!(simulation, max_refinements = 2)

plot(
    plot(simulation.ρ, φ = 0, resample=false),
    plot(simulation.ρ, z = 0.01),
    plot(simulation.ρ, z = 0.03, resample=false),
    plot(simulation.ρ, z = 0.05),
    size = (1200,1200)
)

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
import Test: @test, @testset, @test_throws
import Luna: Utils, @hlock
import HDF5

@testset "Utils" begin

@testset "dict->HDF5" begin
d = Dict{String, Any}()
d["float"] = 1.0
d["float[]"] = [1.0, 2.0, 3.0]
d["string"] = "foo"
d["nothing"] = nothing
d["dict"] = Dict("foo"=>5, "bar"=>[1, 2, 3], "baz"=>Dict("complex"=>5.0+2.0im))
fpath = joinpath(Utils.cachedir(), "output_test", "test.h5")
isfile(fpath) && rm(fpath)
isdir(dirname(fpath)) || mkpath(dirname(fpath))
Utils.save_dict_h5(fpath, d)
@test_throws ErrorException Utils.save_dict_h5(fpath, d, force=false)
@hlock HDF5.h5open(fpath) do file
    for k in ["float", "float[]", "string"]
        @test d[k] == read(file[k])
    end
    @test d["dict"]["foo"] == read(file["dict"]["foo"])
    @test d["dict"]["bar"] == read(file["dict"]["bar"])
    @test d["dict"]["baz"]["complex"] == read(file["dict"]["baz"]["complex"])
end
rm(fpath)

delete!(d, "nothing") # nothing values are converted to empty arrays and not re-converted
Utils.save_dict_h5(fpath, d)
dd = Utils.load_dict_h5(fpath)
@test d == dd

rm(fpath)
rm(dirname(fpath), force=true)
end

@testset "super/subscripts" begin
@test Utils.subscript(0) == Utils.subscript("0") == Utils.subscript('0') == "₀"
@test Utils.subscript(1) == Utils.subscript("1") == Utils.subscript('1') == "₁"
@test Utils.subscript(2) == Utils.subscript("2") == Utils.subscript('2') == "₂"

@test Utils.subscript(123456789) == Utils.subscript("123456789") == "₁₂₃₄₅₆₇₈₉"
@test Utils.subscript("0123456789") == "₀₁₂₃₄₅₆₇₈₉"
end

end


text:
function _read_bibtex(path)
    str = read(path, String)
    _, dict = BibTeX.parse_bibtex(str)
    return dict
end

function _bookbib_to_Ref(dict::AbstractDict)

    ref = Ref()

    # add fields
    set_bibkey!(ref, get(dict, "bibkey", ""))
    set_authors!(ref, get(dict, "author", ""))
    set_year!(ref, get(dict, "year", ""))
    set_title!(ref, get(dict, "title", ""))
    set_doi!(ref, get(dict, "doi", ""))
    set_keywords!(ref, get(dict, "keywords", ""))
    set_abstract!(ref, get(dict, "abstract", ""))
    set_refdict!(ref, dict)

    return ref
end

function _format_BOOKBIB!(_BOOKBIB)

    # clear content
    for (bibkey, dict) in _BOOKBIB

        # doi
        if haskey(dict, "doi")
            doistr = _doi_to_url(dict["doi"])
            !isempty(doistr) && (dict["doi"] = doistr)
        end

        # add bibkeys into dicts
        dict["bibkey"] = bibkey

        for (key, val) in dict
            key == "doi" && continue
            dict[key] = _clear_bibtex_entry(val)
        end
    end
end

function _bibrefs(bibs::Vector{<:AbstractString})

    _BOOKBIB = Dict{String, Any}()
    for path in bibs
        dict = _read_bibtex(path)

        merge!(_BOOKBIB, dict)
    end
    _format_BOOKBIB!(_BOOKBIB)

    # return
    refs = OrderedDict{String, Ref}()
    for (bibkey, dict) in _BOOKBIB
        refs[bibkey] = _bookbib_to_Ref(dict)
    end
    refs
end
_bibrefs(bib::AbstractString) = _bibrefs([bib])

##
# API
bibrefs() = _bibrefs(bibpaths())


text:
# Fixed Talbot method
#
# Abate, J. and Valkó, P.P.
# Multi-precision Laplace transform inversion
# International Journal for Numerical Methods in Engineering, Vol. 60 (Iss. 5-7)  2004  pp 979–993

"""
    talbot(func::Function, t::AbstractFloat, M::Integer=talbot_default_num_terms)

Evaluate the inverse Laplace transform of `func` at the point `t`. Use `M` terms in the algorithm.
For `typeof(t)` is `Float64`, the default for `M` is `32`. For `BigFloat` the default is `64`.

If `BigFloat` precision is larger than default, try increasing `M`.

# Example

```jldoctest
julia> InverseLaplace.talbot(s -> 1 / s^3, 3)
4.50000000000153
```

!!! note
    This function uses the fixed Talbot method. It evaluates `func` for complex arguments.
"""
function talbot(func, t, M)
    bM = convert(typeof(t), M)
    r = (2 * bM) / (5 * t)
    term = (1//2) * exp(r * t) * func(r)
    for i in 1:M-1
        theta = i * (pi/bM)
        s = r * theta * (complex(cot(theta), one(theta)))
        sigma = theta + (theta * cot(theta) - 1) * cot(theta)
        term += real(exp(t * s) * complex(one(t), sigma) * func(s))
    end
    return term * 2 / (5 * t)
end

talbot(func, t) = talbot(func, t, talbot_default_num_terms)
const talbot_BigFloat_default_num_terms = 64
talbot(func, t::BigFloat) = talbot(func, t, talbot_BigFloat_default_num_terms)
talbot(func, t::Integer, args...) = talbot(func, BigFloat(t), args...)
# Hmm, at some point, one of these routines actually gave a Rational result. Don't recall how.
# But, it can't be talbot.
talbot(func, t::Rational, args...) = talbot(func, BigFloat(t), args...)

# Operate on an array of values of t. A single function evaluation
# f(s) is used for all t
# This gives more inaccurate results the further values of
# t are from tmax
"""
    talbotarr(func, ta::AbstractArray, M)

Compute the inverse Laplace transform for each element in `ta`. Each evaluation
of `func(s)` is used for all elements of `ta`. This may be faster than a
broadcast application of `talbot` (i.e. `talbot.(...`) , but is in general, less
accurate.  `talbotarr` uses the "fixed" Talbot method.
"""
function talbotarr(func, t::AbstractArray, M)
    tt = typeof(t[1])
    bM = convert(tt, M)
    terms = similar(t)
    tmax = maximum(t)
    r = (2 * bM) / (5 * tmax)
    fr = (1//2) * func(r)
    terms .= exp.(r .* t) .* fr
    for i in 1:M-1
        theta = i * (pi/bM)
        s = r * theta * (complex(cot(theta), one(theta)))
        sigma = theta + (theta * cot(theta) - 1) * cot(theta)
        fs = complex(one(tt), sigma) * func(s)
        for j in 1:length(terms)
            terms[j] += real(exp(t[j] * s) * fs)
        end
    end
    terms .= terms .* 2 / (5 * tmax)
    return terms
end

talbotarr(func, t::AbstractArray) = talbotarr(func, t, talbot_default_num_terms)
talbotarr(func, t::Vector{BigFloat}) = talbotarr(func, t, talbot_BigFloat_default_num_terms)

#  LocalWords:  talbotarr func talbot