text:
using TensorOperations
using QuadGK
using Zygote
using ForwardDiff
using StaticArrays
using LinearAlgebra

function gaussiandensity_self_energy(rtol=1e-12; tmin=0, tmax=Inf)
    _f(x) = 1/(2x^2+1)^(3//2)
    return 2π^(5//2).*quadgk(_f, tmin, tmax, rtol=rtol)
end

function gaussian_coulomb_integral(a₁=1, a₂=1,  d=0; rtol=1e-12, tmin=0, tmax=Inf)
    @assert a₁>0
    @assert a₂>0
    function _f(t; a₁= a₁, a₂=a₂, d=d)
        c = (a₁+a₂)*t^2 +  a₁*a₂
        return exp(-( a₁*a₂/c)*t^2*d^2)/c^(3//2)
    end
    return 2π^(5//2).*quadgk(_f, tmin, tmax, rtol=rtol)
end

function gaussian_coulomb_integral_grad(a₁=1, a₂=1,  d=0; rtol=1e-12, tmin=0, tmax=Inf)
    function eri_ssss(x)
        a1 = x[1]
        a2 = x[2]
        d = x[3]
        boys(x; n=0) = quadgk( t-> t^(2n)*exp(-x[1]*t^2), 0, 1; rtol=1e-13)[1]
        θ = a1*a2/(a1+a2)
        F₀=boys(θ*d^2)
        return 2π^(5//2)*sqrt(θ)*F₀/(a1*a2)^(3//2)
    end
    h(x)=Zygote.gradient( z->Zygote.forwarddiff(eri_ssss, z), x)[1]
    return h( [Float64(a₁), Float64(a₂), Float64(d)])
end

function gaussian_density_nuclear_potential(a=1, d=0; rtol=1e-12, tmin=0, tmax=Inf)
    #f(t) = exp(-a*t^2*d^2/(t^2+a)) * (t^2+a)^(-3//2)
    f(t) = (t^2+a)^(-3//2)
    return 2π.*quadgk(f, tmin, tmax; rtol=rtol)
end

"""
    test_nuclear_potential(a, ne::Int, ng::Int, nt::Int; Kwargs...) -> Dict

Test nuclear potential accuracy to Gaussien electron density.

# Arguments
- `a`        : Grid box size in bohr or unitful unit if given
- `ne::Int`  : Number of elements per dimension
- `ng::Int`  : Number of Gauss points per element per dimension
- `nt::Int`  : Number of Gauss points for t-integration

# Keywords
- `α=1`           :  Width of Gaussian
- `origin=0.`     :  Nuclear coordinate -> (origin, origin, origin)
- `tmin=0`        :  Minimum t-value
- `tmax=30`       :  Maximum t-value
- `δ=0.25`        :  Local correction parameter, if used
- `mode="normal"` :  Sets mode, options are: "normal", "log", "loglocal", "preset"
"""
function test_nuclear_potential(a, ne::Int, ng::Int, nt::Int;
                                α=1, σ=0.1, origin=zeros(3), tmin=0, tmax=30, δ=0.25, mode="normal")
    ceg = CubicElementGrid(a, ne, ng)
    if mode == "normal"
        @info "normal mode"
        npt1 = NuclearPotentialTensor(origin[1], ceg, nt; tmin=tmin, tmax=tmax)
        npt2 = NuclearPotentialTensor(origin[2], ceg, nt; tmin=tmin, tmax=tmax)
        npt3 = NuclearPotentialTensor(origin[3], ceg, nt; tmin=tmin, tmax=tmax)
    elseif mode == "log"
        @info "logarithmic mode"
        npt1 = NuclearPotentialTensorLog(origin[1], ceg, nt; tmin=tmin, tmax=tmax)
        npt2 = NuclearPotentialTensorLog(origin[2], ceg, nt; tmin=tmin, tmax=tmax)
        npt3 = NuclearPotentialTensorLog(origin[3], ceg, nt; tmin=tmin, tmax=tmax)
    elseif mode == "loglocal"
        @info "logarithmic mode with local correction δ=$δ"
        npt1 = NuclearPotentialTensorLogLocal(origin[1], ceg, nt; tmin=tmin, tmax=tmax, δ=δ)
        npt2 = NuclearPotentialTensorLogLocal(origin[2], ceg, nt; tmin=tmin, tmax=tmax, δ=δ)
        npt3 = NuclearPotentialTensorLogLocal(origin[3], ceg, nt; tmin=tmin, tmax=tmax, δ=δ)
    elseif mode == "preset"
        @info "Preset mode"
        v1 = optimal_nuclear_tensor(ceg, origin[1])
        v2 = optimal_nuclear_tensor(ceg, origin[2])
        v3 = optimal_nuclear_tensor(ceg, origin[3])
        tmin = v1.tmin
        tmax = v1.tmax
        @info "tmin is set to $tmin"
        @info "tmax is set to $tmax"
        vv = nuclear_potential(ceg, 1, v1,v2,v3)
        V = vv.vals
    elseif mode == "gaussian"
        @info "Gaussian mode"
        npt1 = NuclearPotentialTensorGaussian(origin[1], ceg, nt; tmin=tmin, tmax=tmax, σ=σ)
        npt2 = NuclearPotentialTensorGaussian(origin[2], ceg, nt; tmin=tmin, tmax=tmax, σ=σ)
        npt3 = NuclearPotentialTensorGaussian(origin[3], ceg, nt; tmin=tmin, tmax=tmax, σ=σ)
    else
        error("mode not recognized")
    end
    if mode != "preset"
        pt = PotentialTensor(npt1, npt2, npt3)
        V = Array(pt)
    end
    @info "Using using cubic box of size ($a)^3 = $(a^3)"
    @info "Using $ne^3 = $(ne^3) elements"
    @info "Using $ng^3 = $(ng^3) Gauss points per element"
    @info "Total ammount of points is $((ne*ng)^3)"
    @info "t-integration is using $nt points"
    r = position_operator(ceg)
    r0 = r - austrip.(collect(origin)).*u"bohr"
    r2 = r0⋅r0
    ρ = exp(-1u"bohr^-2"*α*r2).vals
    integral = integrate(ρ, ceg, V)
    ref = gaussian_density_nuclear_potential(α; tmin=tmin, tmax=tmax)
    ref_tot = gaussian_density_nuclear_potential(α)
    tail = gaussian_density_nuclear_potential(α; tmin=tmax)
    @info "Integral = $(integral)"
    @info "Reference = $(ref[1])"
    @info "Relative error = $( round((integral-ref[1])/ref[1]; sigdigits=2) )"
    @info "Total reference = $(ref_tot[1])"
    @info "Relative error to total = $( round((integral-ref[1])/ref_tot[1]; sigdigits=2) )"
    @info "Tail energy = $(round(tail[1]; sigdigits=2))"
    @info "Tail relative to total = $(round(tail[1]/ref_tot[1]; sigdigits=2))"
    return Dict("integral"=>integral,
                "reference"=>ref[1],
                "total reference"=>ref_tot[1],
                "tail"=>tail[1])
end

"""
    test_accuracy(a::Real, ne::Int, ng::Int, nt::Int; kwords) -> Dict

Test accuracy on Gaussian charge distribution self energy.

# Arguments
- `a::Real`   : Simulation box size
- `ne::Int`   : Number of elements per dimension
- `ng::Int`   : Number of Gausspoints per dimension for r
- `nt::Int`   : Number of Gausspoints per dimension for t

# Keywords
- `tmax=300`          : Maximum t-value for integration
- `tmin=0`            : Minimum t-value for integration
- `mode=:combination` : Integration type - options `:normal`, `:log`, `:loglocal`, `:local`, `:combination` and `:optimal`
- `δ=0.25`            : Localization parameter for local integration types
- `correction=true`   : Correction to `tmax`-> ∞ integration - `true` calculate correction, `false` do not calculate
- `tboundary=20`      : Parameter for `:combination` mode. Switch to `:loglocal` for t>`tboundary` else use `:log`
- `α1=1`              : 1st Gaussian is exp(-α1*r^2)
- `α2=1`              : 2nd Gaussian is exp(-α2*r^2)
- `d=0`               : Distance between two Gaussian centers
"""
function test_accuracy(a, ne::Int, ng::Int, nt::Int;
        tmax=300, tmin=0, mode=:optimal, δ=0.25, correction=true, tboundary=20, α1=1, α2=1, d=0u"bohr", showprogress=true)
    ceg = CubicElementGrid(a, ne, ng)
    if mode == :normal
        @info "Normal mode"
        ct = CoulombTransformation(ceg, nt; tmax=tmax, tmin=tmin)
    elseif mode == :log
        @info "Logritmic mode"
        ct = CoulombTransformationLog(ceg, nt; tmax=tmax, tmin=tmin)
    elseif mode == :local
        @info "Local mode"
        ct = CoulombTransformationLocal(ceg, nt; tmax=tmax, tmin=tmin, δ=δ)
    elseif mode == :loglocal
        @info "Log local mode"
        ct = CoulombTransformationLogLocal(ceg, nt; tmax=tmax, tmin=tmin, δ=δ)
    elseif mode == :combination
        @info "Combination mode"
        ct = optimal_coulomb_tranformation(ceg, nt; tmax=tmax, δ=δ, tboundary=tboundary)
    elseif  mode == :optimal
        @info "Optimal mode"
        ct = optimal_coulomb_tranformation(ceg, nt)
    else
        error("Mode not known")
    end
    return test_accuracy(ceg, ct; correction=correction, α1=α1, α2=α2, d=d, showprogress=showprogress)
end

function test_accuracy(ceg::CubicElementGrid, ct::AbstractCoulombTransformation;
                            α1=1, α2=1, d=0u"bohr", correction=true, showprogress=true)
    tmin = ct.tmin
    tmax = ct.tmax
    r = position_operator(ceg)
    r1 = r + [0.5, 0., 0.].*d
    r2 = r - [0.5, 0., 0.].*d
    q1 = QuantumState( exp( -0.5α1*u"bohr^-2"*(r1⋅r1) ) )
    q2 = QuantumState( exp( -0.5α1*u"bohr^-2"*(r2⋅r2) ) )
    sd = SlaterDeterminant(q1)
    V = coulomb_operator(sd, ct; correction=correction, showprogress=showprogress)
    V_corr = ScalarOperator(ceg, coulomb_correction(density_operator(q1).vals, ct.tmax); unit=u"hartree")

    E_cor = bracket(q2, V_corr, q2)
    E_int = bracket(q2, V, q2)
    E_tail = gaussian_coulomb_integral(α1, α2, austrip(d); tmin=tmax)[1] * u"hartree"
    E_true = gaussian_coulomb_integral(α1, α2, austrip(d); tmin=tmin, tmax=tmax)[1] * u"hartree"
    if correction
        E = E_int + E_cor
    else
        E = E_int
    end
    E_tot = gaussian_coulomb_integral(α1, α2, austrip(d))[1]*u"hartree"
    @info "Calculated energy = $E"
    @info "True inregration energy = $E_true"
    @info "Total Energy (reference) = $E_tot"
    @info "Integration error = $(round(u"hartree", (E_int-E_true); sigdigits=2)) ; error/E = $( round((E_int-E_true)/E_true; sigdigits=2))"
    @info "Integration error relative to total energy = $( round((E_int-E_true)/E_tot; sigdigits=2))"
    @info "Tail energy = $(round(u"hartree", E_tail; sigdigits=5)) ; E_tail/E_tot = $(round(E_tail/E_tot; sigdigits=2))"
    @info "Energy correction = $(round(u"hartree", E_cor; sigdigits=5))  ; (E_cor-E_tail)/E_tot = $(round((E_cor-E_tail)/E_tot;sigdigits=2))"
    @info "Error to total energy error/E_tot = $(round((E-E_tot)/E_tot; sigdigits=2))"
    out = Dict("calculated"=>E,
        "reference integral"=>E_true,
        "calculated integral"=>E_int,
        "correction"=>E_cor,
        "reference tail"=>E_tail,
        "total reference energy"=>E_tot
        )
    return out
end

function test_accuracy_ad(a, ne, ng, nt; tmax=300, α1=1, α2=1, d=0.5, δ=0.25, showprogress=true)
    ceg = CubicElementGrid(a, ne, ng)
    ct = optimal_coulomb_tranformation(ceg, nt; tmax=tmax, δ=δ, tboundary=20)
    return test_accuracy_ad(ceg::CubicElementGrid, ct::AbstractCoulombTransformation; α1=α1, α2=α2, d=d, showprogress=showprogress)
end

function test_accuracy_ad(ceg::CubicElementGrid, ct::AbstractCoulombTransformation;
                         α1=1., α2=1., d=0.5, showprogress=true)
    function f(x)
        tmin = ct.tmin
        tmax = ct.tmax
        r1 = SVector(0.5x[3], 0., 0.)
        r2 = SVector(-0.5x[3], 0., 0.)
        ρ1 = density_tensor(ceg; a=x[1], r=r1)
        ρ2 = density_tensor(ceg; a=x[2], r=r2)
        V = poisson_equation(ρ1, ct; tmax=tmax, showprogress=showprogress)
        return integrate(ρ2, ceg, V)
    end
    x = [Float64(α1), Float64(α2), Float64(d)]
    @info "Zygote forward mode"
    g = Zygote.gradient( z->Zygote.forwarddiff(f, z), x)[1]
    #@info "ForwardDiff"
    #gf = x -> ForwardDiff.gradient(f, x)
    #fd = gf(x)
    g_ref = gaussian_coulomb_integral_grad(α1, α2, d)
    @info "Zygote forward mode gradient $g"
    #@info "ForwardDiff $fd"
    @info "Analytic $g_ref"
    @info "Relative error  $( round.((g.-g_ref)./g_ref; sigdigits=2))"
    return (g, g_ref)
end

##  Kinetic energy tests
function test_kinetic_energy(a, ne, ng; ν=0, ω=1)
    ceg = CubicElementGrid(a, ne, ng)

    H = HamiltonOperatorFreeParticle(ceg)

    ψ = HarmonicEigenstate(ν; ω=ω)
    ϕ = ReferenceStates.harmonic_state(ceg, ψ)
    T = bracket(ϕ, H, ϕ)
    Tref = 0.5*3*ReferenceStates.energy(ψ)
    @info "Calculated kinetic energy = $T"
    @info "Reference kinetic energy = $Tref"
    @info "Error = $(T-Tref)"
    @info "Relative error = $((T-Tref)/Tref)"
    return T, Tref
end

# Wave function test
function test_atom_wavefunction(a, ne, nt; r=zeros(3), aname="H", ni=20)
    Z = elements[Symbol(aname)].number
    ceg = CubicElementGrid(a, ne, nt)
    ra = austrip.(r) .* u"bohr"
    r0 = position_operator(ceg) - ra
    rr = sqrt(r0⋅r0)

    V = -1u"e_au"*nuclear_potential(ceg, aname, r)

    phi = exp(-0.5u"bohr^-1"*rr)
    ϕ = QuantumState(ceg, phi.vals)
    normalize!(ϕ)

    H = HamiltonOperator(V)
    ψ = helmholtz_equation(ϕ, H; showprogress=true )

    for i in 1:ni
        helmholtz_equation!(ψ, H)
    end
    E = bracket(ψ, H, ψ) |> u"Ry"
    V = bracket(ψ, V, ψ) |> u"Ry"

    E_ref = Z^2 * u"Ry"

    @info "Energy = $E"
    @info "Reference energy = $E_ref"

    return Dict(
        "energy" => E,
        "E_ref"  => E_ref,
        "H"      => H,
        "ψ"      => ψ
    )
end


text:
### TRY numero 2

tile_size = 4
board_size = 8 # must be even
mini_board = [zeros(tile_size, tile_size) ones(tile_size, tile_size);
              ones(tile_size, tile_size) zeros(tile_size, tile_size)]

chessboard = repeat(mini_board, outer=(convert(Integer,(board_size/2)), convert(Integer,(board_size/2))))
chessboard = repeat(mini_board, outer=(2, 2))

img = chessboard

### test clean for u

chess_norm = chessboard
chess_norm_f = Float32.(chess_norm)

chess_warp_f = copy(chess_norm_f)
for k in 2:size(chess_warp_f)[1]
    for l in 2:size(chess_warp_f)[2]
        chess_warp_f[k-1, l-1] = chess_warp_f[k, l]
        if l == size(chess_warp_f)[2]
            chess_warp_f[k-1, l] = chess_warp_f[k, l]
        end
        if k == size(chess_warp_f)[1]
            chess_warp_f[k, l-1] = chess_warp_f[k, l]
        end
    end
end

filter_half_size = 2
filter_size = filter_half_size * 2 + 1

u_est, coeffs = LAP_julia.registration.single_lap(chess_norm_f, chess_warp_f, 3, filter_half_size, [filter_size, filter_size])

window_half_size = [1, 1]
u_out = LAP_julia.helpers.smooth_with_gaussian(u_est[3:end-2, 3:end-2], window_half_size)

LAP_julia.visualise.showflow(u_est[3:end-2, 3:end-2])

LAP_julia.visualise.showflow(u_out)

### test gen rand flow
for k in 1:2
    out = gen_tiled_flow((100,100),10)
    showflow(out, 1)
end

out = gen_tiled_flow((100,100),10)
showflow(out)
# --->works

### Tim Holy INTERPOLATION example

using Interpolations, ColorVectorSpace, ImageFiltering, Images

function imWarp(img, dx, dy)
    itp = interpolate(img, BSpline(Linear()))
    inds = indices_spatial(img)
    rng = extrema.(inds)
    imw = similar(img, eltype(itp))
    println(rng[1])
    for I in CartesianIndices(inds)
        dxi, dyi = dx[I], dy[I]
        y, x = I[1]+dyi, I[2]+dxi
        if is_in(y, rng[1]) && is_in(x,rng[2])
            imw[I] = itp(y, x)
        else
            imw[I] = 0
        end
    end
    return imw
end

function is_in(x, (low, high))
    return x >= low && x <= high
end

is_in(1.0, (1, 512))

# Usage demo
using TestImages
img = testimage("mandril");
img

# Let's use a random warp. We blur it with a Gaussian kernel to make it reasonably smooth.
using ImageFiltering

#test1:
kern = KernelFactors.IIRGaussian((10,10))  # IIRGaussian is fast even for very large σ
dx, dy = imfilter(100*randn(size(img)), kern), imfilter(100*randn(size(img)), kern);
imgw = imWarp(img, dx, dy);
imgw

#test2:
flow = gen_tiled_flow(size(img), 50);
dx, dy = real(flow), imag(flow);
imgw = imWarp(img, dx, dy);
imgw

#test3:
n = 20
grid = ones(512, 512);
@views grid[1:n:end, :] .= 0;
@views grid[:, 1:n:end] .= 0;
fig = PyPlot.figure(dpi = 400, figsize = (5, 5))
PyPlot.imshow(grid, cmap = :gray); gcf()
flow = gen_tiled_flow(size(grid), 50);
dx, dy = real(flow), imag(flow);
imgw = imWarp(grid, dx, dy);

fig = PyPlot.figure(dpi = 400, figsize = (5, 5))
PyPlot.imshow(imgw, cmap = :gray); gcf()

fig = PyPlot.figure(dpi = 400, figsize = (5, 5))
PyPlot.close_figs()
showflow(flow,200)

## pyplot disp image !!!! Does not work for TestImages

PyPlot.imshow(chessboard, cmap = :gray)
gcf()

###
using Interpolations

t = 0:.1:1
x = sin.(2π*t)
y = cos.(2π*t)
A = hcat(x,y)

itp = Interpolations.scale(interpolate(A, (BSpline(Cubic(Natural(OnCell()))), NoInterp())), t, 1:2)

tfine = 0:.01:1
xs, ys = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine];

using PyPlot

figure()
scatter(x, y, label="knots")
PyPlot.plot(xs, ys, label="spline")
gcf()

###
using Interpolations
using PyPlot

t = 0:.1:1
x = sin.(2π*t)

itp = Interpolations.interpolate(x, (BSpline(Quadratic(Natural(OnGrid())))))
itp = Interpolations.interpolate(x, (BSpline(Linear())));
sitp = scale(itp, t);
#sitp = itp

tfine = 0:.01:1;
xs = [sitp(t) for t in tfine]

figure()
scatter(t, x, label="knots")
PyPlot.plot(tfine, xs, label="spline")
gcf()

### get rid of unwanted image replication on the right and bottom of images in the interpolation

img = testimage("mandril_gray");
ix = size(img)[1]
iy = size(img)[2]

ixhalf = Int64(ix/2)
iyhalf = Int64(iy/2)

# flow all to centere
flow = Array{ComplexF64}(undef, ix, iy)
flow[1:ixhalf, 1:iyhalf] .= -20 .+ -20 .* im;
flow[1:ixhalf, iyhalf:iy] .= 20 .+ -20 .* im;
flow[ixhalf:ix, 1:iyhalf] .= -20 .+ 20 .* im;
flow[ixhalf:ix, iyhalf:iy] .= 20 .+ 20 .* im;

flow = LAP_julia.helpers.smooth_with_gaussian(flow, 60);

showflow(flow)

dx, dy = real(flow), imag(flow);
imgw = imWarp(img, dx, dy)
imgw

save("img.png", imgw)

###

ext: jl


text:
using LinearAlgebra
using SparseArrays
using Test
using ProximalOperators

@testset "Quadratic" begin

# Test with full matrices

n, k = 5, 4

A = randn(n, k)
Q = A*A'
q = randn(n)
f = Quadratic(Q, q)
@test typeof(f) <: ProximalOperators.QuadraticDirect

predicates_test(f)

@test ProximalOperators.is_smooth(f) == true
@test ProximalOperators.is_quadratic(f) == true
@test ProximalOperators.is_set(f) == false

x = randn(n)

grad_fx, fx = gradient_test(f, x)
@test fx ≈ 0.5*dot(x, Q*x) + dot(x, q)
@test all(grad_fx .≈ (Q*x + q))

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 1.5)

f = Quadratic(Q, q, iterative=true)
@test typeof(f) <: ProximalOperators.QuadraticIterative

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 2.1)

# Test with sparse matrices

Q = sparse(Q)
f = Quadratic(Q, q)
@test typeof(f) <: ProximalOperators.QuadraticDirect

grad_fx, fx = gradient_test(f, x)
@test fx ≈ 0.5*dot(x, Q*x) + dot(x, q)
@test all(grad_fx .≈ (Q*x + q))

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 0.8)

f = Quadratic(Q, q, iterative=true)
@test typeof(f) <: ProximalOperators.QuadraticIterative

call_test(f, x)
prox_test(f, x)
prox_test(f, x, 1.3)

end


text:
module TokyoVariantAnalysis4
using Logging
using CSV
using DataFrames
using Dates
using HypothesisTests: BinomialTest, confint

export variants_df
export date_to_value, value_to_date

@info "module", @__MODULE__

function has_name(df::DataFrame, name::String)
    return !isnothing(findfirst(isequal(name), names(df)))
end

@info "--------"
@info "作業ディレクトリ"
@info pwd()

@info "--------"
@info "基準日"
const epoch_date = Date("2021-01-01")
@info "epoch_date", epoch_date
#
date_to_value(d::Date) = Dates.value(d - epoch_date)
value_to_date(v)::Date = epoch_date + Day(v)

@info "========"
@info "データ読み込み・変換"

const organization_names = ["k", "p"]   # k 健安研, p 民間検査機関等
const variant_names = ["N501Y", "L452R"]
@info "organization_names", organization_names
@info "variant_names", variant_names

function load_df(csv_filepath)::DataFrame
    @info csv_filepath

    @info "--------"
    @info "CSV 読込み"
    csv_df = CSV.read(csv_filepath, DataFrame; missingstring="NA")
    @info size(csv_df)

    @info "--------"
    @info "データフレーム作成"
    df = DataFrame()
    @info "--------"
    @info "日付 (start_date, end_date, period_days, mid_date, date_value, date_error)"
    @assert has_name(csv_df, "date")
    dates = csv_df.date
    df.date_start  = dates
    df.date_end    = dates
    df.days        = date_to_value.(df.date_end) .+ 1.0 .- date_to_value.(df.date_start)
    df.date_error  = df.days ./ 2.0
    df.date_value  = date_to_value.(df.date_start) .+ df.date_error
    df.date_mid    = value_to_date.(round.(df.date_value .- 0.5))
    @assert all(df.date_start .≤ df.date_mid .≤ df.date_end)
    @info "df.date_start[end]", df.date_start[end]
    @info "df.date_mid[end]", df.date_mid[end]
    @info "df.date_end[end]", df.date_end[end]
    @info size(df)
    @info "--------"
    @info "変異株検査総数のコピー (ORG_n)"
    for orgname in organization_names
        colname = orgname*"_n"  # カラム名
        @assert has_name(csv_df, colname)
        df[!,colname] = csv_df[!,colname]
    end
    @info size(df)
    @info "--------"
    @info "変異株カウントのコピー (ORG_VAR)"
    for orgname in organization_names
        for varname in variant_names
            colname = orgname*"_"*varname   # カラム名
            @assert has_name(csv_df, colname)
            df[!,colname] = csv_df[!,colname]
        end
    end
    @info size(df)

    n_rows = nrow(df)
    @info "--------"
    @info "検査機関にわたる変異株検査総数の合計 (n)"
    ns::Vector{Union{Int,Missing}} = fill(missing, n_rows)
    for i in 1:n_rows
        n = missing
        for orgname in organization_names
            colname = orgname*"_n"  # 検査機関別カラム名
            orgn = df[i,colname]
            if !ismissing(orgn)
                if ismissing(n)
                    n  = orgn
                else
                    n += orgn
                end
            end
        end
        ns[i] = n
    end
    df.n = ns
    @info size(df)
    @info "--------"
    @info "検査機関にわたる変異株カウント合計 (VAR)"
    for varname in variant_names
        cs::Vector{Union{Int,Missing}} = fill(missing, n_rows)
        for i in 1:n_rows
            c = missing
            for orgname in organization_names
                colname = orgname*"_"*varname   # 検査機関別カラム名
                orgc = df[i,colname]
                if !ismissing(orgc)
                    if ismissing(c)
                        c  = orgc
                    else
                        c += orgc
                    end
                end
            end
            cs[i] = c
        end
        df[!,varname] = cs
    end
    @info size(df)

    @info "--------"
    @info "クラスター・渡航関連者を補正した L452R 変異株数と検査総数"
    @assert has_name(csv_df, "L452R_tr")   # 渡航者接触者
    @assert has_name(csv_df, "L452R_cn")   # クラスター数
    @assert has_name(csv_df, "L452R_cp")   # クラスター内人数
    ns_adj = copy(df.n)
    cs_adj = copy(df.L452R)
    for i in 1:n_rows
        n_adj = ns_adj[i]
        c_adj = cs_adj[i]
        tr = csv_df.L452R_tr[i]
        if !ismissing(tr)
            n_adj -= tr
            c_adj -= tr
        end
        cn = csv_df.L452R_cn[i]
        cp = csv_df.L452R_cp[i]
        if !ismissing(cp)
            n_adj -= cp
            c_adj -= cp
            if !ismissing(cn)
                n_adj += cn
                c_adj += cn
            end
        end
        @assert ismissing(n_adj) || (n_adj ≥ 0)
        @assert ismissing(c_adj) || (c_adj ≥ 0)
        ns_adj[i] = n_adj
        cs_adj[i] = c_adj
    end
    df.n     = ns_adj
    df.L452R = cs_adj
    return df
end

const tokyo_variants_csv_filename = "東京変異株4.csv"
variants_df = load_df(tokyo_variants_csv_filename)

@info "========"
@info "統計量算出"

const variant_ratio_pairs = [
    ("L452R", "N501Y"),
]

function statistics!(df::DataFrame)
    n_rows = nrow(df)
    @info "--------"
    @info "全体に対する割合を計算 (VAR_prop)"
    for varname in variant_names
        df[!,varname*"_prop"] = df[!,varname] ./ df.n
    end
    @info size(df)
    @info "--------"
    @info "1 変異株の全体に対する信頼区間を計算 (VAR_prop_cil, VAR_prop_ciu)"
    for varname in variant_names
        cs = df[!,varname]
        ns = df.n
        cils::Vector{Union{Float64,Missing}} = fill(missing, n_rows)
        cius::Vector{Union{Float64,Missing}} = fill(missing, n_rows)
        for i in 1:n_rows
            (cils[i], cius[i]) = if !ismissing(cs[i]) && !ismissing(ns[i])
                confint(BinomialTest(cs[i], ns[i]))
            else
                (missing, missing)
            end
        end
        df[!,varname*"_prop_cil"] = cils
        df[!,varname*"_prop_ciu"] = cius
    end
    @info size(df)
    @info "--------"
    @info "1 変異株のオッズを計算 (VAR_odds)"
    for varname in variant_names
        df[!,varname*"_odds"] = 1.0 ./ ((df.n ./ df[!,varname]) .- 1.0)
    end
    @info size(df)
    @info "--------"
    @info "オッズの信頼区間を計算 (VAR_odds_cil, VAR_odds_ciu)"
    for varname in variant_names
        cs = df[!,varname]
        ns = df.n
        cils::Vector{Union{Float64,Missing}} = fill(missing, n_rows)
        cius::Vector{Union{Float64,Missing}} = fill(missing, n_rows)
        for i in 1:n_rows
            (cils[i], cius[i]) = if !ismissing(cs[i]) && !ismissing(ns[i])
                1.0 ./ (1.0 ./ confint(BinomialTest(cs[i], ns[i])) .- 1.0)
            else
                (missing, missing)
            end
        end
        df[!,varname*"_odds_cil"] = cils
        df[!,varname*"_odds_ciu"] = cius
    end
    @info size(df)
    @info "--------"
    @info "オッズの対数（ロジット） (log_VAR_odds, log_VAR_odds_cil, log_VAR_odds_ciu)"
    for varname in variant_names
        oddname = varname*"_odds"
        df[!,"log_"*oddname]        = log.(df[!,oddname])
        df[!,"log_"*oddname*"_cil"] = log.(df[!,oddname*"_cil"])
        df[!,"log_"*oddname*"_ciu"] = log.(df[!,oddname*"_ciu"])
    end
    @info "--------"
    @info "特定 2 変異株間の比を計算 (VAR1_VAR2)"
    for (varname1, varname2) in variant_ratio_pairs
        @info varname1, varname2
        df[!,varname1*"_"*varname2] = df[!,varname1] ./ df[!,varname2]
    end
    @info size(df)
    @info "--------"
    @info "2 変異株間の比の信頼区間を計算 (VAR1_VAR2_cil, VAR1_VAR2_ciu)"
    @info "2 変異株のみからなる 2 項分布による近似"
    for (varname1, varname2) in variant_ratio_pairs
        @info varname1, varname2
        cs1 = df[!,varname1]
        cs2 = df[!,varname2]
        cils::Vector{Union{Float64,Missing}} = fill(missing, n_rows)
        cius::Vector{Union{Float64,Missing}} = fill(missing, n_rows)
        for i in 1:n_rows
            (cils[i], cius[i]) = if !ismissing(cs1[i]) && !ismissing(cs2[i])
                1.0 ./ (1.0 ./ confint(BinomialTest(cs1[i], cs1[i] + cs2[i])) .- 1.0)
            else
                (missing, missing)
            end
        end
        df[!,varname1*"_"*varname2*"_cil"] = cils
        df[!,varname1*"_"*varname2*"_ciu"] = cius
    end
    @info size(df)
    @info "--------"
    @info "2 変異株間比の対数 (log_VAR1_VAR2, log_VAR1_VAR2_cil, log_VAR1_VAR2_ciu)"
    for (varname1, varname2) in variant_ratio_pairs
        @info varname1, varname2
        pairname = varname1*"_"*varname2
        df[!,"log_"*pairname]        = log.(df[!,pairname])
        df[!,"log_"*pairname*"_cil"] = log.(df[!,pairname*"_cil"])
        df[!,"log_"*pairname*"_ciu"] = log.(df[!,pairname*"_ciu"])
    end
    @info size(df)
end

statistics!(variants_df)

@info "--------"
@info "解析データ書き出し"
const output_csv_filename = "tokyo_variant_analysis4.csv"
@info output_csv_filename
CSV.write(output_csv_filename, variants_df; missingstring="NA")
@info size(variants_df)

end #module

max_stars_repo_name: muimuiz/TokyoVariant

max_issues_repo_name: muimuiz/TokyoVariant

max_forks_repo_name: muimuiz/TokyoVariant


text:
f = open("2019/01/input.txt", "r")
inlines = [parse(Int, x) for x in readlines(f)]
close(f)

outlines = [[(x÷3)-2] for x in inlines]

while true
    going = false
    for l in outlines
        fuel = (last(l)÷3)-2
        if fuel > 0
            #println(fuel)
            push!(l, fuel)
            going = true
        end
    end
    if !going
        break
    end
end

println(sum([sum(l) for l in outlines]))

max_stars_repo_name: chaserobertson/advent

max_issues_repo_name: chaserobertson/advent

max_forks_repo_name: chaserobertson/advent


text:
export srDistance

"""
function srDistance

Unnormalized distance (reconstruction error) between two pairs (M,w):

min_P ||M1*diagm(w1) - M2*diagm(w2)*P ||_p

where P is a permutation matrix and || ||_p denotes the entrywise matrix pn-norm
(i.e., pn=2 corresponds to the Frobenius norm).

Before computing the error, the first dimensions of the binary matrices are
augmented from m*(p-1) to m*p so that the column sums of each block become 1.
This can be disabled by setting p=1.

If n has different values among the inputs, then the largest n is chosen and
the rest are augmented with zeros to make the dimensions compatible.

Input:
    M1 - binary array of size (m*(p-1),n) or (m*(p-1),n,N)
    w1 - frequency vector of length n or an array of size (n,N)
    M2 - as M1
    w2 - as w1

Keyword arguments:
    p - number of categories (default=2)
    pn - exponent in the matrix norm (default=2.0)

Output:
    distance(s) as a scalar or an N-vector
"""

function srDistance(M1::BitArray{2},w1::Vector,M2::BitArray{2},w2::Vector;p::Int64=2,pn::Float64=2.0)
    n = maximum((size(M1,2),length(w1),size(M2,2),length(w2)))
    M1 = [M1 falses(size(M1,1),n-size(M1,2))]
    w1 = [w1; zeros(n-length(w1))]
    M2 = [M2 falses(size(M2,1),n-size(M2,2))]
    w2 = [w2; zeros(n-length(w2))]

    perms = permutations(1:n)

    if p>1
        M1 = fullCategoricalMatrix(M1,p)
        M2 = fullCategoricalMatrix(M2,p)
    end

    distance = Inf
    for perm in perms
        M2perm = M2[:,perm]
        w2perm = w2[perm]

        mat = M1*diagm(w1) - M2perm*diagm(w2perm)
        dist = (sum(abs.(mat[:]).^pn))^(1/pn)
        if dist < distance
            distance = dist
        end
    end

    return distance
end

function srDistance(M1::BitArray{3},w1::Array{Float64,2},M2::BitArray{3},w2::Array{Float64,2};p::Int64=2,pn::Float64=2.0)
    N = size(M1,3)
    distances = SharedArray{Float64}(N)

    @sync @parallel for idx = 1:N
        distances[idx] = srDistance(M1[:,:,idx],w1[:,idx],M2[:,:,idx],w2[:,idx],p=p,pn=pn)
    end

    return distances
end

function fullCategoricalMatrix(M::BitArray{2},p::Int64)
    n = size(M,2)

    M1 = reshape(M,(p-1,:))
    blockcolsums = sum(M1,1)
    Mfull = [1-blockcolsums; M1]

    return reshape(Mfull,(:,n))
end


text:
@testset "preprocess" begin
    s = """
        @def z = [1,2,3,
            4,5,6]
        """
    tokens = F.find_tokens(s, F.MD_TOKENS, F.MD_1C_TOKENS)
    F.find_indented_blocks!(tokens, s)

    @test tokens[1].name == :MD_DEF_OPEN
    @test tokens[2].name == :LR_INDENT
    @test tokens[3].name == :LINE_RETURN
    @test tokens[4].name == :EOS

    F.preprocess_candidate_mddefs!(tokens)

    @test tokens[1].name == :MD_DEF_OPEN
    @test tokens[2].name == :LINE_RETURN
    @test tokens[3].name == :EOS

    blocks, = F.find_all_ocblocks(tokens, F.MD_OCB_ALL)

    content = F.content.(blocks) .|> String
    @test isapproxstr(content[1], "z = [1,2,3,4,5,6]")

    s = """
        @def z1 = [1,2,3,
            4,5,6]
        @def z2 = 3
        @def z3 = [1,
            2,
                3]
        """
    tokens = F.find_tokens(s, F.MD_TOKENS, F.MD_1C_TOKENS)
    F.find_indented_blocks!(tokens, s)
    seqtok = (:MD_DEF_OPEN, :LR_INDENT, :LINE_RETURN,
              :MD_DEF_OPEN, :LINE_RETURN,
              :MD_DEF_OPEN, :LR_INDENT, :LR_INDENT, :LINE_RETURN,
              :EOS)
    for i in 1:length(tokens)
        @test tokens[i].name == seqtok[i]
    end

    F.preprocess_candidate_mddefs!(tokens)
    seqtok = (:MD_DEF_OPEN, :LINE_RETURN,
              :MD_DEF_OPEN, :LINE_RETURN,
              :MD_DEF_OPEN, :LINE_RETURN,
              :EOS)
    for i in 1:length(tokens)
        @test tokens[i].name == seqtok[i]
    end

    blocks, = F.find_all_ocblocks(tokens, F.MD_OCB_ALL)
    content = F.content.(blocks) .|> String
    @test isapproxstr(content[1], "z1=[1,2,3,4,5,6]")
    @test isapproxstr(content[2], "z2=3")
    @test isapproxstr(content[3], "z3=[1,2,3]")
end

@testset "mddefs1" begin
    F.def_LOCAL_VARS!()
    s = """
        @def x = 5
        """ |> fd2html_td
    @test F.locvar("x") == 5
    s = """
        @def x = "hello"; <!-- nothing -->
        """ |> fd2html_td
    @test F.locvar("x") == "hello"
    s = """
        @def x = "hello"
        @def y = pi
        """ |> fd2html_td
    @test F.locvar("x") == "hello"
    @test F.locvar("y") == pi
    s = """
        @def z = [1,2,3,
                  4,5,6]
        A
        """ |> fd2html_td
    @test F.locvar("z") == collect(1:6)
    s = """
        @def s = \"\"\"foo
            bar
            baz etc\"\"\"
        @def s2 = "nothing"
        """ |> fd2html_td
    @test isapproxstr(F.locvar("s"), "foo bar baz etc")
    @test isapproxstr(F.locvar("s2"), "nothing")
    @test F.locvar("foobar") === nothing
    @test F.globvar("foobar") === nothing
end


text:
using YaoExtensions, Yao
using Random, Test

@testset "EchoBlock" begin
    block = EchoBlock(4, :Test)
    reg = rand_state(4)
    @test reg |> block == reg
    @test mat(block) != nothing
    @test ishermitian(block)
    @test isreflexive(block)
    @test isunitary(block)
    @test getiparams(block) == ()

    @test EchoBlock()(5) isa EchoBlock
    @test EchoBlock()(5).sym == :ECHO
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
"""
Bridges `MOI.Indicator{A, CP.DifferentFrom}` to linear constraints
(including, possibly, strict inequalities). This constraint adds one variable
to store the absolute value of the difference, and uses it for the indicator.

For `AbstractFloat` arguments (like `Float64`): equivalent to abs(x) > 0.0, i.e.
a `Strictly(GreaterThan(0.0))`.
"""
struct IndicatorDifferentFrom2PseudoMILPBridge{T <: Real, A} <: MOIBC.AbstractBridge
    var_abs::Union{Nothing, MOI.VariableIndex}
    con_abs::Union{Nothing, MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.AbsoluteValue}, MOI.ConstraintIndex{MOI.VariableIndex, CP.AbsoluteValue}}
    con_indic::Union{
        MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}, T}}},
        MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}}},
        MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.EqualTo{T}}},
    }
end

function MOIBC.bridge_constraint(
    ::Type{IndicatorDifferentFrom2PseudoMILPBridge{T, A}},
    model,
    f::MOI.VectorOfVariables,
    s::MOI.Indicator{A, CP.DifferentFrom{T}},
) where {T <: Real, A}
    return MOIBC.bridge_constraint(
        IndicatorDifferentFrom2PseudoMILPBridge{T, A},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{IndicatorDifferentFrom2PseudoMILPBridge{T, A}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::MOI.Indicator{A, CP.DifferentFrom{T}},
) where {T <: Real, A}
    f_scalars = MOIU.scalarize(f)

    # Add the absolute value.
    var_abs = MOI.add_variable(model)
    con_abs = MOI.add_constraint(
        model,
        MOIU.vectorize([
            MOI.ScalarAffineFunction(
                [MOI.ScalarAffineTerm(one(T), var_abs)],
                zero(T),
            ),
            f_scalars[2] - s.set.value,
        ]),
        CP.AbsoluteValue(),
    )

    # New indicator. For floats, the difference may take any real value:
    # > 0 must be encoded as-is.
    con_indic = MOI.add_constraint(
        model,
        MOIU.vectorize(
            [
                f_scalars[1],
                one(T) * var_abs
            ]
        ),
        MOI.Indicator{A}(CP.Strictly(MOI.GreaterThan(zero(T))))
    )

    return IndicatorDifferentFrom2PseudoMILPBridge{T, A}(var_abs, con_abs, con_indic)
end

function MOIBC.bridge_constraint(
    ::Type{IndicatorDifferentFrom2PseudoMILPBridge{T, A}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::MOI.Indicator{A, CP.DifferentFrom{T}},
) where {T <: Integer, A}
    f_scalars = MOIU.scalarize(f)

    # Add the absolute value.
    var_abs = MOI.add_variable(model)
    con_abs = MOI.add_constraint(
        model,
        MOIU.vectorize([one(T) * var_abs, f_scalars[2] - s.set.value]),
        CP.AbsoluteValue(),
    )

    # New indicator. For integers, the difference takes an integer value:
    # > 0 is equivalent to >= 1.
    con_indic = MOI.add_constraint(
        model,
        MOIU.vectorize(
            [
                f_scalars[1],
                one(T) * var_abs
            ]
        ),
        MOI.Indicator{A}(MOI.GreaterThan(one(T)))
    )

    return IndicatorDifferentFrom2PseudoMILPBridge{T, A}(var_abs, con_abs, con_indic)
end

function MOIBC.bridge_constraint(
    ::Type{IndicatorDifferentFrom2PseudoMILPBridge{Bool, A}},
    model,
    f::MOI.VectorAffineFunction{Bool},
    s::MOI.Indicator{A, CP.DifferentFrom{Bool}},
) where {A}
    # New indicator. For Booleans, no need to have an absolute value:
    # != 0 is equivalent to == 1, and vice-versa.
    con_indic = MOI.add_constraint(
        model,
        f,
        MOI.Indicator{A}(MOI.EqualTo{Bool}(1 - s.set.value))
    )

    return IndicatorDifferentFrom2PseudoMILPBridge{Bool, A}(nothing, nothing, con_indic)
end

function MOI.supports_constraint(
    ::Type{IndicatorDifferentFrom2PseudoMILPBridge{T, A}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{MOI.Indicator{A, CP.DifferentFrom{T}}},
) where {T <: Real, A}
    return true
end

function MOIB.added_constrained_variable_types(::Type{IndicatorDifferentFrom2PseudoMILPBridge{T, A}}) where {T <: Real, A}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{IndicatorDifferentFrom2PseudoMILPBridge{T, A}}) where {T <: AbstractFloat, A}
    return [
        (MOI.VectorAffineFunction{T}, CP.AbsoluteValue),
        (MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}}}),
    ]
end

function MOIB.added_constraint_types(::Type{IndicatorDifferentFrom2PseudoMILPBridge{T, A}}) where {T <: Integer, A}
    return [
        (MOI.VectorAffineFunction{T}, CP.AbsoluteValue),
        (MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}}),
    ]
end

function MOIB.added_constraint_types(::Type{IndicatorDifferentFrom2PseudoMILPBridge{Bool, A}}) where {A}
    return [
        (MOI.VectorAffineFunction{Bool}, MOI.Indicator{A, MOI.EqualTo{Bool}}),
    ]
end

function MOI.get(::IndicatorDifferentFrom2PseudoMILPBridge{T, A}, ::MOI.NumberOfVariables) where {T <: Real, A}
    return 1
end

function MOI.get(::IndicatorDifferentFrom2PseudoMILPBridge{Bool, A}, ::MOI.NumberOfVariables) where {A}
    return 0
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.AbsoluteValue,
    },
) where {T <: Real, A}
    return 1
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{Bool, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{Bool}, CP.AbsoluteValue,
    },
) where {A}
    return 0
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}}},
    },
) where {T <: AbstractFloat, A}
    return 1
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}}},
    },
) where {T <: Real, A}
    return 0
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}},
    },
) where {T <: AbstractFloat, A}
    return 0
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}},
    },
) where {T <: Real, A}
    return 1
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{Bool, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{Bool}, MOI.Indicator{A, MOI.GreaterThan{Bool}},
    },
) where {A}
    return 0
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.EqualTo{T}},
    },
) where {T <: Real, A}
    return 0
end

function MOI.get(
    ::IndicatorDifferentFrom2PseudoMILPBridge{Bool, A},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{Bool}, MOI.Indicator{A, MOI.EqualTo{Bool}},
    },
) where {A}
    return 1
end

function MOI.get(
    b::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.ListOfVariableIndices,
) where {T <: Real, A}
    return [b.var_abs]
end

function MOI.get(
    b::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.AbsoluteValue,
    },
) where {T <: Real, A}
    return [b.con_abs]
end

function MOI.get(
    b::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, CP.Strictly{MOI.GreaterThan{T}}},
    },
) where {T <: Real, A}
    return [b.con_indic]
end

function MOI.get(
    b::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.GreaterThan{T}},
    },
) where {T <: Real, A}
    return [b.con_indic]
end

function MOI.get(
    b::IndicatorDifferentFrom2PseudoMILPBridge{T, A},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{A, MOI.EqualTo{T}},
    },
) where {T <: Real, A}
    return [b.con_indic]
end


text:
include("utils.jl")
include("../Applications/Smuggler.jl")

#ENV["JULIA_DEBUG"] = Main
# for n in 20:20:100
#     for seed in 1:10
#         generate_data(n, seed, 0.5)
#     end
# end

n = 20
k = 2
inst_id = 1
density = 0.5
heuristicmode = "HEURCG"
const constraint_uncertainty = false

instance = "$(inst_id)_$(n)_$(density).txt"
data = read_data(instance)
K = 1:k
N = 1:data.n
B = 5
Ω = 1
ρ = 0.5

# The following excutions avoid warming-up issues
@info "Executing functions to warm up ..."
oldstd = stdout
redirect_stdout(open("/dev/null", "w"))
reset_timer!(to::TimerOutput)
STARTTIME = time()
scenario_generation(false,true)
redirect_stdout(oldstd) # recover original stdout
@info "warm-up finished!"

#Now comes the instance that is tested, passed by the arguments
n = parse(Int, ARGS[1])
k = parse(Int, ARGS[2])
inst_id = parse(Int, ARGS[3])

instance = "$(inst_id)_$(n)_$(density).txt"
outputfilename = "/../results/Smuggler_scenario_generation_$(n)_$(k)_$(B)_$(density).csv"
instance_being_solved = "Solving instance $inst_id with k=$k and density=$density and ϵ=$ϵ ..."
@info instance_being_solved
data = read_data(instance)
K = 1:k
N = 1:data.n

reset_timer!(to::TimerOutput)
STARTTIME = time()
runtime = @elapsed lb, ub, solution, gap, NUMBER_GENSC, NUMITER, NUMBER_GENSOL = scenario_generation(false,false)
outputfile = open(string(dirname(@__FILE__),outputfilename), "a")
output = "$n $inst_id $k $NUMBER_GENSC $(data.nX_before_dom) $NUMBER_GENSOL $ϵ $runtime $NUMITER $lb $ub $gap $(TimerOutputs.time(to["heuristic"])/10^9) $(TimerOutputs.time(to["gen_sol"])/10^9) $(data.heuristic_ub)"
println(output)
println()
println(outputfile, output)
close(outputfile)
show(to)

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
mutable struct AppUser
    id::Union{Missing,String}
    login::Union{Missing,String}
    password::Union{Missing,String}
    lastname::Union{Missing,String}
    firstname::Union{Missing,String}
    email::Union{Missing,String}

    # Convenience constructor that allows us to create a vector of instances
    #   from a JuliaDB.table using the dot syntax: `Myclass.(a_JuliaDB_table)`
    AppUser(args::NamedTuple) = AppUser(;args...)
    AppUser(;id = missing,
             login = missing,
             password = missing,
             lastname = missing,
             firstname = missing,
             email = missing) = (
                  # First call the default constructor with missing values only so that
                  #   there is no risk that we don't assign an argument to the wrong attribute
                  x = new(missing,missing,missing,missing,missing,
                          missing);
                  x.id = id;
                  x.login = login;
                  x.password = password;
                  x.lastname = lastname;
                  x.firstname = firstname;
                  x.email = email;

                  return x )

end


text:
@resumable function retireProcess( sim::Simulation, mpSim::MPsim )

    retirement = mpSim.retirement

    if ( retirement.retirementAge == 0.0 ) &&
        ( retirement.maxCareerLength == 0.0 )
        return
    end  # if ( retirement.retirementAge == 0.0 ) && ...

    # Timing of the process.
    processTime = Dates.Millisecond( 0 )
    tStart = now()
    retName = string( "Default retirement process " )

    # Preparatory steps.
    priority = typemax( Int )
    timeToWait = retirement.offset

    # If an initial population snapshot is uploaded, and it contains zero-time
    #   events, don't execute zero-time retirement.
    if ( timeToWait == 0 ) && !mpSim.isVirgin
        timeToWait += retirement.freq
    end  # if ( timeToWait == 0 ) && ...

    # Process loop.
    while now( sim ) + timeToWait <= mpSim.simLength
        processTime += now() - tStart
        @yield timeout( sim, timeToWait, priority = priority )
        tStart = now()
        retirementCycle( retirement, mpSim )
        timeToWait = retirement.freq
    end  # while now( sim ) + timeToWait <= mpSim.simLength

    # Timing of the process.
    processTime += now() - tStart

    if mpSim.showInfo
        println( retName, "took ", processTime.value / 1000, " seconds." )
    end  # if mpSim.showInfo

end  # @resumable function retireProcess( sim, mpSim )

function retirementCycle( retirement::Retirement, mpSim::MPsim )

    currentTime = now( mpSim )
    queryCmd = string( "SELECT id, currentNode FROM `", mpSim.persDBname,
        "` WHERE",
        "\n    status IS 'active' AND" )

    if retirement.retirementAge == 0.0
        queryCmd = string( queryCmd,
            "\n    ", currentTime, " - timeEntered >= ",
            retirement.maxCareerLength )
    elseif retirement.maxCareerLength == 0.0
        queryCmd = string( queryCmd,
            "\n    ", currentTime, " - timeEntered + ageAtRecruitment >= ",
            retirement.retirementAge )
    else
        queryCmd = string( queryCmd,
            "\n    (", currentTime, " - timeEntered >= ",
            retirement.maxCareerLength, retirement.isEither ? " OR " : " AND ",
            "\n        ", currentTime, " - timeEntered + ageAtRecruitment >= ",
            retirement.retirementAge, ")" )
    end  # if retirement.retirementAge == 0.0

    queryCmd = string( queryCmd,
        "\n    ORDER BY currentNode" )
    ids = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )

    if isempty( ids )
        return
    end  # if isempty( ids )

    currentNodes = Vector{String}( ids[:, :currentNode] )
    ids = Vector{String}( ids[:, Symbol( mpSim.idKey )] )
    removePersons( ids, currentNodes, "retirement", mpSim )

end  # retirementCycle( retirement, mpSim )

function removePersons( ids::Vector{String}, currentNodes::Vector{String},
    simTimes::Union{Float64, Vector{Float64}}, reason::String, mpSim::MPsim )

    # Add the retirement transition to the transition records.
    sqliteCmd = string.( "\n    ('", ids, "', ", simTimes, ", '", reason,
        "', '", currentNodes, "')" )
    sqliteCmd = string( "INSERT INTO `", mpSim.transDBname, "` (`", mpSim.idKey,
        "`, timeIndex, transition, sourceNode) VALUES", join( sqliteCmd, "," ) )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    # Clear the ids from the inNodeSince field of their current node.
    for name in unique( currentNodes )
        currentNode = mpSim.baseNodeList[name]
        delete!.( Ref( currentNode.inNodeSince ), ids[currentNodes .== name] )
    end  # for name in unique( currentNodes )

    # Update the personnel records.
    sqliteCmd = string.( "UPDATE `", Ref(mpSim.persDBname), "`",
        "\n    SET status = '", Ref(reason), "',",
        "\n        timeExited = ", simTimes, ",",
        "\n        currentNode = NULL,",
        "\n        inNodeSince = NULL",
        "\n    WHERE `", mpSim.idKey, "` IS '", ids, "'" )
    DBInterface.execute.( Ref(mpSim.simDB), sqliteCmd )

end  # removePersons( ids, currentNodes, reason, mpSim )

removePersons( ids::Vector{String}, currentNode::String,
    simTimes::Union{Float64, Vector{Float64}}, reason::String, mpSim::MPsim ) =
    removePersons( ids, fill( currentNode, length(ids) ), simTimes, reason,
    mpSim )

removePersons( ids::Vector{String},
    currentNodes::Union{String, Vector{String}}, reason::String,
    mpSim::MPsim ) = removePersons( ids, currentNodes, now(mpSim), reason, mpSim )

removePerson( id::String, currentNode::String, reason::String, mpSim::MPsim ) =
    removePersons( [id], currentNode, reason, mpSim )


text:
module DeepReshapes

export describe, deep_length, deep_reshape, flatten, pack

# describe(x) recursively describes the dimensions of x
#   1. so that some deep reshapes can be reversed and
#   2. to determine the number of scalar values contained within (from
#      deep_length()).
#
# At the moment, the only recognized scalars are Numbers.
#
# In the untyped version (typed = false) this is a generalization of
# Base.size(): the specification of
#   - a scalar is an empty tuple,
#   - an Array of scalars is the Tuple of its dimensions,
#   - any other Array is the Array{Any} of specifications of its elements,
#   - a non-empty Tuple is the Tuple of specification of its elements,
#   - anything else is not defined.
#
# When typed = true, the specification of
#   - a scalar is its type,
#   - an Array of scalars is the Tuple of its element type concatenated with its
#     dimensions,
#   - and of all other objects is the same as in the untyped case.

describe{T <: Number}(x::T; typed::Bool = true) = typed ? T : ()

describe(x::Tuple{Any, Vararg{Any}}; typed::Bool = true) =
  tuple([describe(e, typed = typed) for e in x]...)

describe{T <: Number}(x::Array{T}; typed::Bool = true) =
  typed ? tuple(T, Base.size(x)...) : Base.size(x)

describe(x::Array; typed::Bool = true) =
  map!(e -> describe(e, typed = typed), similar(x, Any), x)
  # ^ We use map! here because we want to define the element type of the result.

# A deep_reshape(x, specification) iterates over all the scalars (recursively)
# contained within x, and from that stream writes into a new object whose
# structure is defined by specification.
#
# This happens in two independent processes which are executed in parallel:
#   - deep_produce() recursively walks the input and produces the scalars, and
#   - consume_structured() recursively builds the result from them according to
#     specification.
#
# By default, deep_produce() descends into Arrays and Tuples, but that behavior
# can be overridden by supplying the optional Deep argument. All encountered
# values that are not of Deep type are considered scalar by deep_produce() and
# therefore produced directly.
#
# The specification is of the same format that describe() returns.

DeepDefault = Union{Array, Tuple}

deep_reshape(x, specification; Deep = DeepDefault) =
  consume_structured(@task(deep_produce(x, Deep)), specification)

deep_reshape(x, specification...) = deep_reshape(x, specification)

# Descend into all x of Deep type.
deep_produce{T}(x::T, Deep::Type{T}) =
  for e in x; deep_produce(e, Deep) end

# Produce scalars directly.
deep_produce(x, ::Type) = produce(x)

# Consume a scalar or Array of scalars.
function consume_structured(
  producer::Task,
  Member::DataType,
  dimensions::Integer...
)
  if dimensions == ()
    convert(Member, consume(producer))
  else
    result = Array(Member, dimensions...)
    for indices in CartesianRange(dimensions)
      result[indices] = consume(producer)
    end
    result
  end
end

# Consume a scalar or Array of scalars of unspecified type.
consume_structured(producer::Task, specification::Tuple{Vararg{Integer}}) =
  consume_structured(producer, Any, specification...)

# Consume a scalar or Array of scalars of specific type.
consume_structured(producer::Task, specification::Tuple{DataType, Vararg{Integer}}) =
  consume_structured(producer, specification...)

# Recursively consume any other Array.
function consume_structured(producer::Task, specification::Array)
  dimensions = size(specification)
  result = Array(Any, dimensions...)
  for indices in CartesianRange(dimensions)
    result[indices] = consume_structured(producer, specification[indices])
  end
  result
end

# Recursively consume a Tuple.
consume_structured(producer::Task, specification::Tuple) =
  tuple([consume_structured(producer, element) for element in specification]...)

# deep_length(x) returns the count of scalar values (recursively) contained
# within x. Only works where describe(x) also works.

deep_length(x) = specification_length(describe(x, typed = false))

specification_length(::Tuple{}) = 1
specification_length(x::Tuple{Vararg{Integer}}) = prod(x)
specification_length(x) = sum(specification_length, x)

# This is just a convenience wrapper to reshape any deep-reshapable into a flat
# Array, optionally with a fixed result element type.

flatten(MemberType::DataType, x...) =
  deep_reshape(x, MemberType, deep_length(x))

flatten(x...) = flatten(Any, x...)

# A convenience wrapper to ease converting between deep and flat representations
# of the same data.

pack(MemberType::DataType, data...) =
  (flatten(MemberType, data), describe(data))

pack(x...) = pack(Any, x...)

end # module


text:
# This file was generated, do not modify it. # hide
G = SearchGraph(; dist, db)
index!(G)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
# machine code generation

## libdevice

function find_libdevice(cap)
    CUDAnative.configured || return
    global libdevice

    if isa(libdevice, Dict)
        # select the most recent & compatible library
        vers = keys(CUDAnative.libdevice)
        compat_vers = Set(ver for ver in vers if ver <= cap)
        isempty(compat_vers) && error("No compatible CUDA device library available")
        ver = maximum(compat_vers)
        path = libdevice[ver]
    else
        libdevice
    end
end

const libdevices = Dict{String, LLVM.Module}()
function load_libdevice(ctx::CompilerContext)
    path = find_libdevice(ctx.cap)

    get!(libdevices, path) do
        open(path) do io
            libdevice = parse(LLVM.Module, read(io), JuliaContext())
            name!(libdevice, "libdevice")
            libdevice
        end
    end
end

function link_libdevice!(ctx::CompilerContext, mod::LLVM.Module, libdevice::LLVM.Module)
    # override libdevice's triple and datalayout to avoid warnings
    triple!(libdevice, triple(mod))
    datalayout!(libdevice, datalayout(mod))

    # 1. save list of external functions
    exports = String[]
    for f in functions(mod)
        fn = LLVM.name(f)
        if !haskey(functions(libdevice), fn)
            push!(exports, fn)
        end
    end

    # 2. link with libdevice
    link!(mod, libdevice)

    ModulePassManager() do pm
        # 3. internalize all functions not in list from (1)
        internalize!(pm, exports)

        # 4. eliminate all unused internal functions
        #
        # this isn't necessary, as we do the same in optimize! to inline kernel wrappers,
        # but it results _much_ smaller modules which are easier to handle on optimize=false
        global_optimizer!(pm)
        global_dce!(pm)
        strip_dead_prototypes!(pm)

        # 5. run NVVMReflect pass
        push!(metadata(mod), "nvvm-reflect-ftz",
              MDNode([ConstantInt(Int32(1), JuliaContext())]))

        # 6. run standard optimization pipeline
        #
        #    see `optimize!`

        run!(pm, mod)
    end
end

## PTX code generation

function machine(cap::VersionNumber, triple::String)
    InitializeNVPTXTarget()
    InitializeNVPTXTargetInfo()
    t = Target(triple)

    InitializeNVPTXTargetMC()
    cpu = "sm_$(cap.major)$(cap.minor)"
    if cuda_driver_version >= v"9.0" && v"6.0" in ptx_support
        # in the case of CUDA 9, we use sync intrinsics from PTX ISA 6.0+
        feat = "+ptx60"
    else
        feat = ""
    end
    tm = TargetMachine(t, triple, cpu, feat)
    asm_verbosity!(tm, true)

    return tm
end

function mcgen(ctx::CompilerContext, mod::LLVM.Module, f::LLVM.Function)
    tm = machine(ctx.cap, triple(mod))

    InitializeNVPTXAsmPrinter()
    return String(emit(tm, mod, LLVM.API.LLVMAssemblyFile))
end


text:
using Mustache
using Test

@testset " interpolation " begin

	## Mustache-free templates should render as-is.
tpl = """Hello from {Mustache}!
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """Hello from {Mustache}!
"""

	## Unadorned tags should interpolate content into the template.
tpl = """Hello, {{subject}}!
"""

	@test Mustache.render(tpl, Dict{Any,Any}("subject"=>"world")) == """Hello, world!
"""

	## Basic interpolation should be HTML escaped.
tpl = """These characters should be HTML escaped: {{forbidden}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("forbidden"=>"& \" < >")) == """These characters should be HTML escaped: & " < >
"""

	## Triple mustaches should interpolate without HTML escaping.
tpl = """These characters should not be HTML escaped: {{{forbidden}}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("forbidden"=>"& \" < >")) == """These characters should not be HTML escaped: & \" < >
"""

	## Ampersand should interpolate without HTML escaping.
tpl = """These characters should not be HTML escaped: {{&forbidden}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("forbidden"=>"& \" < >")) == """These characters should not be HTML escaped: & \" < >
"""

	## Integers should interpolate seamlessly.
tpl = """\"{{mph}} miles an hour!\""""

	@test Mustache.render(tpl, Dict{Any,Any}("mph"=>85)) == """\"85 miles an hour!\""""

	## Integers should interpolate seamlessly.
tpl = """\"{{{mph}}} miles an hour!\""""

	@test Mustache.render(tpl, Dict{Any,Any}("mph"=>85)) == """\"85 miles an hour!\""""

	## Integers should interpolate seamlessly.
tpl = """\"{{&mph}} miles an hour!\""""

	@test Mustache.render(tpl, Dict{Any,Any}("mph"=>85)) == """\"85 miles an hour!\""""

	## Decimals should interpolate seamlessly with proper significance.
tpl = """\"{{power}} jiggawatts!\""""

	@test Mustache.render(tpl, Dict{Any,Any}("power"=>1.21)) == """\"1.21 jiggawatts!\""""

	## Decimals should interpolate seamlessly with proper significance.
tpl = """\"{{{power}}} jiggawatts!\""""

	@test Mustache.render(tpl, Dict{Any,Any}("power"=>1.21)) == """\"1.21 jiggawatts!\""""

	## Decimals should interpolate seamlessly with proper significance.
tpl = """\"{{&power}} jiggawatts!\""""

	@test Mustache.render(tpl, Dict{Any,Any}("power"=>1.21)) == """\"1.21 jiggawatts!\""""

	## Failed context lookups should default to empty strings.
tpl = """I ({{cannot}}) be seen!"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """I () be seen!"""

	## Failed context lookups should default to empty strings.
tpl = """I ({{{cannot}}}) be seen!"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """I () be seen!"""

	## Failed context lookups should default to empty strings.
tpl = """I ({{&cannot}}) be seen!"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """I () be seen!"""

	## Dotted names should be considered a form of shorthand for sections.
tpl = """\"{{person.name}}\" == \"{{#person}}{{name}}{{/person}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("person"=>Dict{Any,Any}("name"=>"Joe"))) == """\"Joe\" == \"Joe\""""

	## Dotted names should be considered a form of shorthand for sections.
tpl = """\"{{{person.name}}}\" == \"{{#person}}{{{name}}}{{/person}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("person"=>Dict{Any,Any}("name"=>"Joe"))) == """\"Joe\" == \"Joe\""""

	## Dotted names should be considered a form of shorthand for sections.
tpl = """\"{{&person.name}}\" == \"{{#person}}{{&name}}{{/person}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("person"=>Dict{Any,Any}("name"=>"Joe"))) == """\"Joe\" == \"Joe\""""

	## Dotted names should be functional to any level of nesting.
tpl = """\"{{a.b.c.d.e.name}}\" == \"Phil\""""

	@test Mustache.render(tpl, Dict{Any,Any}("a"=>Dict{Any,Any}("b"=>Dict{Any,Any}("c"=>Dict{Any,Any}("d"=>Dict{Any,Any}("e"=>Dict{Any,Any}("name"=>"Phil"))))))) == """\"Phil\" == \"Phil\""""

	## Any falsey value prior to the last part of the name should yield ''.
tpl = """\"{{a.b.c}}\" == \"\""""

	@test Mustache.render(tpl, Dict{Any,Any}("a"=>Dict{Any,Any}())) == """\"\" == \"\""""

	## Each part of a dotted name should resolve only against its parent.
tpl = """\"{{a.b.c.name}}\" == \"\""""

	@test Mustache.render(tpl, Dict{Any,Any}("c"=>Dict{Any,Any}("name"=>"Jim"),"a"=>Dict{Any,Any}("b"=>Dict{Any,Any}()))) == """\"\" == \"\""""

	## The first part of a dotted name should resolve as any other name.
tpl = """\"{{#a}}{{b.c.d.e.name}}{{/a}}\" == \"Phil\""""

	@test Mustache.render(tpl, Dict{Any,Any}("b"=>Dict{Any,Any}("c"=>Dict{Any,Any}("d"=>Dict{Any,Any}("e"=>Dict{Any,Any}("name"=>"Wrong")))),"a"=>Dict{Any,Any}("b"=>Dict{Any,Any}("c"=>Dict{Any,Any}("d"=>Dict{Any,Any}("e"=>Dict{Any,Any}("name"=>"Phil"))))))) == """\"Phil\" == \"Phil\""""

	## Dotted names should be resolved against former resolutions.
tpl = """{{#a}}{{b.c}}{{/a}}"""

	@test Mustache.render(tpl, Dict{Any,Any}("b"=>Dict{Any,Any}("c"=>"ERROR"),"a"=>Dict{Any,Any}("b"=>Dict{Any,Any}()))) == """"""

	## Interpolation should not alter surrounding whitespace.
tpl = """| {{string}} |"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """| --- |"""

	## Interpolation should not alter surrounding whitespace.
tpl = """| {{{string}}} |"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """| --- |"""

	## Interpolation should not alter surrounding whitespace.
tpl = """| {{&string}} |"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """| --- |"""

	## Standalone interpolation should not alter surrounding whitespace.
tpl = """  {{string}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """  ---
"""

	## Standalone interpolation should not alter surrounding whitespace.
tpl = """  {{{string}}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """  ---
"""

	## Standalone interpolation should not alter surrounding whitespace.
tpl = """  {{&string}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """  ---
"""

	## Superfluous in-tag whitespace should be ignored.
tpl = """|{{ string }}|"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """|---|"""

	## Superfluous in-tag whitespace should be ignored.
tpl = """|{{{ string }}}|"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """|---|"""

	## Superfluous in-tag whitespace should be ignored.
tpl = """|{{& string }}|"""

	@test Mustache.render(tpl, Dict{Any,Any}("string"=>"---")) == """|---|"""
end


text:
# ---
# title: 136. Single Number
# id: problem136
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: Hash Table, Bit Manipulation
# link: <https://leetcode.com/problems/single-number/description/>
# hidden: true
# ---
#
# Given a **non-empty**  array of integers `nums`, every element appears _twice_
# except for one. Find that single one.
#
# **Follow up:**  Could you implement a solution with a linear runtime
# complexity and without using extra memory?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,2,1]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: nums = [4,1,2,1,2]
#     Output: 4
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 3 * 104`
#   * `-3 * 104 <= nums[i] <= 3 * 104`
#   * Each element in the array appears twice except for one element which appears only once.
#
#
## @lc code=start
using LeetCode

function single_number_136(nums::Vector{Int})::Int
    a = 0
    for i in nums
        a ⊻= i
    end
    a
end
## @lc code=end

## @lc test=start
@testset "136.single-number.jl" begin
    @test single_number_136([4,1,2,1,2]) == 4
    @test single_number_136([2, 2, 1]) == 1
end
## @lc test=end

## @lc info=start
# link: [solution to 136](https://leetcode-cn.com/problems/single-number/solution/136-zhi-chu-xian-yi-ci-de-shu-zi-yi-huo-ukdwx/)
## @lc info=end


text:
using WignerFamilies
using Documenter

makedocs(;
    modules=[WignerFamilies],
    authors="xzackli <xzackli@gmail.com> and contributors",
    repo="https://github.com/xzackli/WignerFamilies.jl/blob/{commit}{path}#L{line}",
    sitename="WignerFamilies.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xzackli.github.io/WignerFamilies.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/xzackli/WignerFamilies.jl",
)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Spectrum
export Uniform, DeltaFunction, Measured

using ....OpticSim
using ...Emitters
using DataFrames
using Distributions
import Unitful: Length, ustrip
using Unitful.DefaultSymbols

const UNIFORMSHORT = 0.450 #um
const UNIFORMLONG = 0.680 #um

abstract type AbstractSpectrum{T<:Real} end

"""
    Uniform{T} <: AbstractSpectrum{T}

Encapsulates a flat spectrum range which is sampled uniformly. Unless stated diferrently, the range used will be 450nm to 680nm.

```julia
Uniform(low_end::T, high_end::T) where {T<:Real}
Uniform(::Type{T} = Float64) where {T<:Real}
```
"""
struct Uniform{T} <: AbstractSpectrum{T}
    low_end::T
    high_end::T

    # user defined range of spectrum
    function Uniform(low_end::T, high_end::T) where {T<:Real}
        return new{T}(low_end, high_end)
    end

    # with no specific range we will use the constants' values
    function Uniform(::Type{T} = Float64) where {T<:Real}
        return new{T}(UNIFORMSHORT, UNIFORMLONG)
    end
end

Emitters.generate(s::Uniform{T}) where {T<:Real} = (one(T), rand(Distributions.Uniform(s.low_end, s.high_end)))

"""
    DeltaFunction{T} <: AbstractSpectrum{T}

Encapsulates a constant spectrum.

```julia
DeltaFunction{T<:Real}
```
"""
struct DeltaFunction{T} <: AbstractSpectrum{T}
    λ::T
end

DeltaFunction(λ::Length) = DeltaFunction{Float64}(ustrip(μm, λ))

Emitters.generate(s::DeltaFunction{T}) where {T<:Real} = (one(T), s.λ)

"""
    Measured{T} <: AbstractSpectrum{T}

Encapsulates a measured spectrum to compute emitter power. Create spectrum by reading CSV files.
Evaluate spectrum at arbitrary wavelength with [`spectrumpower`](@ref) (**more technical details coming soon**)

```julia
Measured(samples::DataFrame)
```
"""
struct Measured{T} <: AbstractSpectrum{T}
    low_wave_length::Int64
    high_wave_length::Int64
    wave_length_step::Int64
    power_samples::Vector{T}

    function Measured(samples::DataFrame)
        colnames = names(samples)
        @assert "Wavelength" in colnames
        @assert "Power" in colnames
        wavelengths = samples[!, :Wavelength]
        @assert eltype(wavelengths) <: Int64

        power::Vector{T} where {T<:Real} = samples[!, :Power] # no missing values allowed and must be real numbers
        maxpower = maximum(power)
        T = eltype(power)
        p = sortperm(wavelengths)
        wavelengths = wavelengths[p]
        power = power[p] ./ maxpower #make sure power values have the same sorted permutation as wavelengths normalized with maximum power equal to 1
        λmin = wavelengths[p[1]]
        λmax = wavelengths[p[end]]
        step = wavelengths[2] - wavelengths[1]

        for i in 3:length(wavelengths)
            @assert wavelengths[i] - wavelengths[i - 1] == step  # no missing values allowed and step size must be consistent
        end

        return new{T}(λmin, λmax, step, power)
    end
end

function Emitters.generate(spectrum::Measured{T}) where {T<:Real}
    λ = rand(Distributions.Uniform(convert(T, spectrum.low_wave_length), convert(T, spectrum.high_wave_length)))
    power = spectrumpower(spectrum, λ)
    if power === nothing        #added this condition because the compiler was allocating if just returned values directly.
        return (nothing, nothing)
    else
        return (power, λ / convert(T, 1000)) #convert from nm to um
    end
end

"""expects wavelength in nm not um"""
function spectrumpower(spectrum::Measured{T}, λ::T) where {T<:Real}
    if λ < spectrum.low_wave_length || λ > spectrum.high_wave_length
        return nothing
    end

    lowindex = floor(Int64, (λ - spectrum.low_wave_length)) ÷ spectrum.wave_length_step + 1

    if lowindex == length(spectrum.power_samples)
        return convert(T, spectrum.power_samples[end])
    else
        highindex = lowindex + 1
        α = mod(λ, spectrum.wave_length_step) / spectrum.wave_length_step

        return convert(T, (1 - α) * spectrum.power_samples[lowindex] + α * spectrum.power_samples[highindex])
    end
end

end # module Spectrum


text:
using Distributions

# ρ(x,y) in integrand
ρ(x,y) = (0 < x < 1) && (0 < y < 2) ? 1/2 : 0

#=
Implementation of random walk sampler Metropolis-Hastings.
Inputs:
- n: The number of states to return
- b: The standard deviation of the normal distribution used in the proposal
     distribution
Outputs:
- Returns an array of length n denoting states S_i in Metropolis-Hastings alg.
=#
function metropolis_hastings(n, b)

    # Samples
    samples = zeros(n,2)
    samples[1,:] = rand(2)

    # Conditional Q distribution
    d = Normal(0,b)

    acceptances = 0

    for i=2:n
        # Propose Y
        ξx = rand(d)
        ξy = rand(d)
        X = samples[i-1,1] + ξx
        Y = samples[i-1,2] + ξy

        # Acceptance probability
        α = min(1, ρ(X,Y)/ρ(samples[i-1,1],samples[i-1,2]))

        # Next state
        if rand() < α
            samples[i,1] = X
            samples[i,2] = Y
        else
            samples[i,:] = samples[i-1,:]
        end
    end

    return samples
end

# Obtain 10000 samples with b = 1
S = metropolis_hastings(10000,1)

# f(x) in integrand
f(x,y) = exp(sin(x*y))

# Estimate integral
M = 1000
int_estimate = 2*mean([f(x,y) for (x,y) in zip(S[M:end,1], S[M:end,2])])

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
module SuperFrame

end

max_stars_repo_name: Mechachleopteryx/JuliaLessons

max_issues_repo_name: Mechachleopteryx/JuliaLessons

max_forks_repo_name: Mechachleopteryx/JuliaLessons


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module Consts

    const HEAD_FILE  = "HEAD"
    const FETCH_HEAD  = "FETCH_HEAD"
    const REMOTE_ORIGIN = "origin"

    # objs
    @enum(OBJECT,
          OBJ_ANY    = -2,
          OBJ_BAD    = -1,
          OBJ_COMMIT = 1,
          OBJ_TREE   = 2,
          OBJ_BLOB   = 3,
          OBJ_TAG    = 4)

    #revwalk
    const SORT_NONE        = Cint(0)
    const SORT_TOPOLOGICAL = Cint(1 << 0)
    const SORT_TIME        = Cint(1 << 1)
    const SORT_REVERSE     = Cint(1 << 2)

    # refs
    const REF_INVALID  = Cint(0)
    const REF_OID      = Cint(1)
    const REF_SYMBOLIC = Cint(2)
    const REF_LISTALL  = REF_OID | REF_SYMBOLIC

    # checkout
    const CHECKOUT_NONE                    = Cuint(0)
    const CHECKOUT_SAFE                    = Cuint(1 << 0)
    const CHECKOUT_FORCE                   = Cuint(1 << 1)
    const CHECKOUT_RECREATE_MISSING        = Cuint(1 << 2)
    const CHECKOUT_ALLOW_CONFLICTS         = Cuint(1 << 4)
    const CHECKOUT_REMOVE_UNTRACKED        = Cuint(1 << 5)
    const CHECKOUT_REMOVE_IGNORED          = Cuint(1 << 6)
    const CHECKOUT_UPDATE_ONLY             = Cuint(1 << 7)
    const CHECKOUT_DONT_UPDATE_INDEX       = Cuint(1 << 8)
    const CHECKOUT_NO_REFRESH              = Cuint(1 << 9)
    const CHECKOUT_SKIP_UNMERGED           = Cuint(1 << 10)
    const CHECKOUT_USE_OURS                = Cuint(1 << 11)
    const CHECKOUT_USE_THEIRS              = Cuint(1 << 12)
    const CHECKOUT_DISABLE_PATHSPEC_MATCH  = Cuint(1 << 13)
    const CHECKOUT_SKIP_LOCKED_DIRECTORIES = Cuint(1 << 18)
    const CHECKOUT_DONT_OVERWRITE_IGNORED  = Cuint(1 << 19)
    const CHECKOUT_CONFLICT_STYLE_MERGE    = Cuint(1 << 20)
    const CHECKOUT_CONFLICT_STYLE_DIFF3    = Cuint(1 << 21)
    const CHECKOUT_DONT_REMOVE_EXISTING    = Cuint(1 << 22)

    const CHECKOUT_UPDATE_SUBMODULES       = Cuint(1 << 16)
    const CHECKOUT_UPDATE_SUBMODULES_IF_CHANGED = Cuint(1 << 17)

    const CHECKOUT_NOTIFY_NONE      = Cuint(0)
    const CHECKOUT_NOTIFY_CONFLICT  = Cuint(1 << 0)
    const CHECKOUT_NOTIFY_DIRTY     = Cuint(1 << 1)
    const CHECKOUT_NOTIFY_UPDATED   = Cuint(1 << 2)
    const CHECKOUT_NOTIFY_UNTRACKED = Cuint(1 << 3)
    const CHECKOUT_NOTIFY_IGNORED   = Cuint(1 << 4)
    const CHECKOUT_NOTIFY_ALL       = 0x0FFFF

    # diff
    const DIFF_OPTIONS_VERSION = Cuint(1)

    const DIFF_NORMAL                     = Cuint(0)
    const DIFF_REVERSE                    = Cuint(1 << 0)
    const DIFF_INCLUDE_IGNORED            = Cuint(1 << 1)
    const DIFF_RECURSE_IGNORED_DIRS       = Cuint(1 << 2)
    const DIFF_INCLUDE_UNTRACKED          = Cuint(1 << 3)
    const DIFF_RECURSE_UNTRACKED_DIRS     = Cuint(1 << 4)
    const DIFF_INCLUDE_UNMODIFIED         = Cuint(1 << 5)
    const DIFF_INCLUDE_TYPECHANGE         = Cuint(1 << 6)
    const DIFF_INCLUDE_TYPECHANGE_TREES   = Cuint(1 << 7)
    const DIFF_IGNORE_FILEMODE            = Cuint(1 << 8)
    const DIFF_IGNORE_SUBMODULES          = Cuint(1 << 9)
    const DIFF_IGNORE_CASE                = Cuint(1 << 10)
    const DIFF_DISABLE_PATHSPEC_MATCH     = Cuint(1 << 12)
    const DIFF_SKIP_BINARY_CHECK          = Cuint(1 << 13)
    const DIFF_ENABLE_FAST_UNTRACKED_DIRS = Cuint(1 << 14)

    const DIFF_FORCE_TEXT               = Cuint(1 << 20)
    const DIFF_FORCE_BINARY             = Cuint(1 << 21)
    const DIFF_IGNORE_WHITESPACE        = Cuint(1 << 22)
    const DIFF_IGNORE_WHITESPACE_CHANGE = Cuint(1 << 23)
    const DIFF_IGNORE_WHITESPACE_EOL    = Cuint(1 << 24)
    const DIFF_SHOW_UNTRACKED_CONTENT   = Cuint(1 << 25)
    const DIFF_SHOW_UNMODIFIED          = Cuint(1 << 26)
    const DIFF_PATIENCE                 = Cuint(1 << 28)
    const DIFF_MINIMAL                  = Cuint(1 << 29)

    const DIFF_FLAG_BINARY     = Cuint(1 << 0)
    const DIFF_FLAG_NOT_BINARY = Cuint(1 << 1)
    const DIFF_FLAG_VALID_OID  = Cuint(1 << 2)

    const DIFF_FORMAT_PATCH        = Cuint(1)
    const DIFF_FORMAT_PATCH_HEADER = Cuint(2)
    const DIFF_FORMAT_RAW          = Cuint(3)
    const DIFF_FORMAT_NAME_ONLY    = Cuint(4)
    const DIFF_FORMAT_NAME_STATUS  = Cuint(5)

    @enum(DELTA_STATUS, DELTA_UNMODIFIED = Cint(0),
                        DELTA_ADDED      = Cint(1),
                        DELTA_DELETED    = Cint(2),
                        DELTA_MODIFIED   = Cint(3),
                        DELTA_RENAMED    = Cint(4),
                        DELTA_COPIED     = Cint(5),
                        DELTA_IGNORED    = Cint(6),
                        DELTA_UNTRACKED  = Cint(7),
                        DELTA_TYPECHANGE = Cint(8))

    # index
    const IDXENTRY_NAMEMASK   = (0x0fff)
    const IDXENTRY_STAGEMASK  = (0x3000)
    const IDXENTRY_EXTENDED   = (0x4000)
    const IDXENTRY_VALID      = (0x8000)
    const IDXENTRY_STAGESHIFT = Cint(12)

    const IDXENTRY_UPDATE            = Cint(1 << 0)
    const IDXENTRY_REMOVE            = Cint(1 << 1)
    const IDXENTRY_UPTODATE          = Cint(1 << 2)
    const IDXENTRY_ADDED             = Cint(1 << 3)

    const IDXENTRY_HASHED            = Cint(1 << 4)
    const IDXENTRY_UNHASHED          = Cint(1 << 5)
    const IDXENTRY_WT_REMOVE         = Cint(1 << 6)
    const IDXENTRY_CONFLICTED        = Cint(1 << 7)

    const IDXENTRY_UNPACKED          = Cint(1 << 8)
    const IDXENTRY_NEW_SKIP_WORKTREE = Cint(1 << 9)

    const INDEXCAP_IGNORE_CASE = Cuint(1)
    const INDEXCAP_NO_FILEMODE = Cuint(2)
    const INDEXCAP_NO_SYMLINKS = Cuint(4)
    const INDEXCAP_FROM_OWNER  = ~Cuint(0)

    const INDEX_ADD_DEFAULT                = Cuint(0)
    const INDEX_ADD_FORCE                  = Cuint(1 << 0)
    const INDEX_ADD_DISABLE_PATHSPEC_MATCH = Cuint(1 << 1)
    const INDEX_ADD_CHECK_PATHSPEC         = Cuint(1 << 2)

    const INDEX_STAGE_ANY = Cint(-1)

    # merge
    @enum(GIT_MERGE, MERGE_FIND_RENAMES     = 1 << 0,
                     MERGE_FAIL_ON_CONFLICT = 1 << 1,
                     MERGE_SKIP_REUC        = 1 << 2,
                     MERGE_NO_RECURSIVE     = 1 << 3)

    @enum(GIT_MERGE_FILE, MERGE_FILE_DEFAULT                  = 0,       # Defaults
                          MERGE_FILE_STYLE_MERGE              = 1 << 0,  # Create standard conflicted merge files
                          MERGE_FILE_STYLE_DIFF3              = 1 << 1,  # Create diff3-style files
                          MERGE_FILE_SIMPLIFY_ALNUM           = 1 << 2,  # Condense non-alphanumeric regions for simplified diff file
                          MERGE_FILE_IGNORE_WHITESPACE        = 1 << 3,  # Ignore all whitespace
                          MERGE_FILE_IGNORE_WHITESPACE_CHANGE = 1 << 4,  # Ignore changes in amount of whitespace
                          MERGE_FILE_IGNORE_WHITESPACE_EOL    = 1 << 5,  # Ignore whitespace at end of line
                          MERGE_FILE_DIFF_PATIENCE            = 1 << 6,  # Use the "patience diff" algorithm
                          MERGE_FILE_DIFF_MINIMAL             = 1 << 7)  # Take extra time to find minimal diff

    @enum(GIT_MERGE_FILE_FAVOR, MERGE_FILE_FAVOR_NORMAL = 0,
                                MERGE_FILE_FAVOR_OURS   = 1,
                                MERGE_FILE_FAVOR_THEIRS = 2,
                                MERGE_FILE_FAVOR_UNION  = 3)

    @enum(GIT_MERGE_PREFERENCE, MERGE_PREFERENCE_NONE             = 0,
                                MERGE_PREFERENCE_NO_FASTFORWARD   = 1,
                                MERGE_PREFERENCE_FASTFORWARD_ONLY = 2)

    @enum(GIT_MERGE_ANALYSIS, MERGE_ANALYSIS_NONE        = 0,
                              MERGE_ANALYSIS_NORMAL      = 1 << 0,
                              MERGE_ANALYSIS_UP_TO_DATE  = 1 << 1,
                              MERGE_ANALYSIS_FASTFORWARD = 1 << 2,
                              MERGE_ANALYSIS_UNBORN      = 1 << 3)

    # reset
    const RESET_SOFT  = Cint(1) # Move the head to the given commit
    const RESET_MIXED = Cint(2) # SOFT plus reset index to the commit
    const RESET_HARD  = Cint(3) # MIXED plus changes in working tree discarded

    #rebase
    @enum(GIT_REBASE_OPERATION, REBASE_OPERATION_PICK   = Cint(0),
                                REBASE_OPERATION_REWORD = Cint(1),
                                REBASE_OPERATION_EDIT   = Cint(2),
                                REBASE_OPERATION_SQUASH = Cint(3),
                                REBASE_OPERATION_FIXUP  = Cint(4),
                                REBASE_OPERATION_EXEC   = Cint(5))

    # fetch_prune
    const FETCH_PRUNE_UNSPECIFIED = Cint(0)
    const FETCH_PRUNE             = Cint(1)
    const FETCH_NO_PRUNE          = Cint(2)

    # remote_autotag
    const REMOTE_DOWNLOAD_TAGS_UNSPECIFIED = Cint(0)
    const REMOTE_DOWNLOAD_TAGS_AUTO        = Cint(1)
    const REMOTE_DOWNLOAD_TAGS_NONE        = Cint(2)
    const REMOTE_DOWNLOAD_TAGS_ALL         = Cint(3)

    # clone
    const CLONE_LOCAL_AUTO     = Cint(0)
    const CLONE_LOCAL          = Cint(1)
    const CLONE_NO_LOCAL       = Cint(2)
    const CLONE_LOCAL_NO_LINKS = Cint(3)

    # status
    const STATUS_CURRENT          = Cuint(0)
    const STATUS_INDEX_NEW        = Cuint(1 << 0)
    const STATUS_INDEX_MODIFIED   = Cuint(1 << 1)
    const STATUS_INDEX_DELETED    = Cuint(1 << 2)
    const STATUS_INDEX_RENAMED    = Cuint(1 << 3)
    const STATUS_INDEX_TYPECHANGE = Cuint(1 << 4)
    const STATUS_WT_NEW           = Cuint(1 << 7)
    const STATUS_WT_MODIFIED      = Cuint(1 << 8)
    const STATUS_WT_DELETED       = Cuint(1 << 9)
    const STATUS_WT_TYPECHANGE    = Cuint(1 << 10)
    const STATUS_WT_RENAMED       = Cuint(1 << 11)
    const STATUS_WT_UNREADABLE    = Cuint(1 << 12)
    const STATUS_IGNORED          = Cuint(1 << 14)
    const STATUS_CONFLICTED       = Cuint(1 << 15)

    # status show
    const STATUS_SHOW_INDEX_AND_WORKDIR = Cint(0)
    const STATUS_SHOW_INDEX_ONLY        = Cint(1)
    const STATUS_SHOW_WORKDIR_ONLY      = Cint(2)

    # status options
    const STATUS_OPT_INCLUDE_UNTRACKED                = Cuint(1 << 0)
    const STATUS_OPT_INCLUDE_IGNORED                  = Cuint(1 << 1)
    const STATUS_OPT_INCLUDE_UNMODIFIED               = Cuint(1 << 2)
    const STATUS_OPT_EXCLUDE_SUBMODULES               = Cuint(1 << 3)
    const STATUS_OPT_RECURSE_UNTRACKED_DIRS           = Cuint(1 << 4)
    const STATUS_OPT_DISABLE_PATHSPEC_MATCH           = Cuint(1 << 5)
    const STATUS_OPT_RECURSE_IGNORED_DIRS             = Cuint(1 << 6)
    const STATUS_OPT_RENAMES_HEAD_TO_INDEX            = Cuint(1 << 7)
    const STATUS_OPT_RENAMES_INDEX_TO_WORKDIR         = Cuint(1 << 8)
    const STATUS_OPT_SORT_CASE_SENSITIVELY            = Cuint(1 << 9)
    const STATUS_OPT_SORT_CASE_INSENSITIVELY          = Cuint(1 << 10)
    const STATUS_OPT_RENAMES_FROM_REWRITES            = Cuint(1 << 11)
    const STATUS_OPT_NO_REFRESH                       = Cuint(1 << 12)
    const STATUS_OPT_UPDATE_INDEX                     = Cuint(1 << 13)
    const STATUS_OPT_INCLUDE_UNREADABLE               = Cuint(1 << 14)
    const STATUS_OPT_INCLUDE_UNREADABLE_AS_UNTRACKED  = Cuint(1 << 15)

    @enum(GIT_SUBMODULE_IGNORE, SUBMODULE_IGNORE_UNSPECIFIED  = -1, # use the submodule's configuration
                                SUBMODULE_IGNORE_NONE         = 1,  # any change or untracked == dirty
                                SUBMODULE_IGNORE_UNTRACKED    = 2,  # dirty if tracked files change
                                SUBMODULE_IGNORE_DIRTY        = 3,  # only dirty if HEAD moved
                                SUBMODULE_IGNORE_ALL          = 4)  # never dirty

    """
Option flags for `GitRepo`.

* `REPOSITORY_OPEN_NO_SEARCH` - Only open the repository if it can be immediately found in the `path`.  Do not walk up from the `path` looking at parent directories.
* `REPOSITORY_OPEN_CROSS_FS` - Unless this flag is set, open will not continue searching across filesystem boundaries. (E.g. Searching in a user's home directory `/home/user/source/` will not return `/.git/` as the found repo if `/` is a different filesystem than `/home`.)
* `REPOSITORY_OPEN_BARE` - Open repository as a bare repo regardless of core.bare config, and defer loading config file for faster setup.
    """
    @enum(GIT_REPOSITORY_OPEN, REPOSITORY_OPEN_DEFAULT   = 0,
                               REPOSITORY_OPEN_NO_SEARCH = 1<<0,
                               REPOSITORY_OPEN_CROSS_FS  = 1<<1,
                               REPOSITORY_OPEN_BARE      = 1<<2)

    @enum(GIT_BRANCH, BRANCH_LOCAL = 1, BRANCH_REMOTE = 2)

    @enum(GIT_FILEMODE, FILEMODE_UNREADABLE          = 0o000000,
                        FILEMODE_TREE                = 0o040000,
                        FILEMODE_BLOB                = 0o100644,
                        FILEMODE_BLOB_EXECUTABLE     = 0o100755,
                        FILEMODE_LINK                = 0o120000,
                        FILEMODE_COMMIT              = 0o160000)

    @enum(GIT_CREDTYPE, CREDTYPE_USERPASS_PLAINTEXT = Cuint(1 << 0),
                        CREDTYPE_SSH_KEY            = Cuint(1 << 1),
                        CREDTYPE_SSH_CUSTOM         = Cuint(1 << 2),
                        CREDTYPE_DEFAULT            = Cuint(1 << 3),
                        CREDTYPE_SSH_INTERACTIVE    = Cuint(1 << 4),
                        CREDTYPE_USERNAME           = Cuint(1 << 5),
                        CREDTYPE_SSH_MEMORY         = Cuint(1 << 6))

    @enum(GIT_FEATURE, FEATURE_THREADS = Cuint(1 << 0),
                       FEATURE_HTTPS   = Cuint(1 << 1),
                       FEATURE_SSH     = Cuint(1 << 2),
                       FEATURE_NSEC    = Cuint(1 << 3))

if LibGit2.version() >= v"0.24.0"
    """
Priority level of a config file.

These priority levels correspond to the natural escalation logic (from higher to lower) when searching for config entries in git.

* `CONFIG_LEVEL_DEFAULT` - Open the global, XDG and system configuration files if any available.
* `CONFIG_LEVEL_PROGRAMDATA` - System-wide on Windows, for compatibility with portable git
* `CONFIG_LEVEL_SYSTEM` - System-wide configuration file; `/etc/gitconfig` on Linux systems
* `CONFIG_LEVEL_XDG` - XDG compatible configuration file; typically `~/.config/git/config`
* `CONFIG_LEVEL_GLOBAL` - User-specific configuration file (also called Global configuration file); typically `~/.gitconfig`
* `CONFIG_LEVEL_LOCAL` - Repository specific configuration file; `\$WORK_DIR/.git/config` on non-bare repos
* `CONFIG_LEVEL_APP` - Application specific configuration file; freely defined by applications
* `CONFIG_HIGHEST_LEVEL` - Represents the highest level available config file (i.e. the most specific config file available that actually is loaded)
    """
    @enum(GIT_CONFIG, CONFIG_LEVEL_DEFAULT     = 0,
                      CONFIG_LEVEL_PROGRAMDATA = 1,
                      CONFIG_LEVEL_SYSTEM      = 2,
                      CONFIG_LEVEL_XDG         = 3,
                      CONFIG_LEVEL_GLOBAL      = 4,
                      CONFIG_LEVEL_LOCAL       = 5,
                      CONFIG_LEVEL_APP         = 6,
                      CONFIG_HIGHEST_LEVEL     =-1)
else
    """
Priority level of a config file.

These priority levels correspond to the natural escalation logic (from higher to lower) when searching for config entries in git.

* `CONFIG_LEVEL_DEFAULT` - Open the global, XDG and system configuration files if any available.
* `CONFIG_LEVEL_SYSTEM` - System-wide configuration file; `/etc/gitconfig` on Linux systems
* `CONFIG_LEVEL_XDG` - XDG compatible configuration file; typically `~/.config/git/config`
* `CONFIG_LEVEL_GLOBAL` - User-specific configuration file (also called Global configuration file); typically `~/.gitconfig`
* `CONFIG_LEVEL_LOCAL` - Repository specific configuration file; `\$WORK_DIR/.git/config` on non-bare repos
* `CONFIG_LEVEL_APP` - Application specific configuration file; freely defined by applications
* `CONFIG_HIGHEST_LEVEL` - Represents the highest level available config file (i.e. the most specific config file available that actually is loaded)
    """
    @enum(GIT_CONFIG, CONFIG_LEVEL_DEFAULT     = 0,
                      CONFIG_LEVEL_SYSTEM      = 1,
                      CONFIG_LEVEL_XDG         = 2,
                      CONFIG_LEVEL_GLOBAL      = 3,
                      CONFIG_LEVEL_LOCAL       = 4,
                      CONFIG_LEVEL_APP         = 5,
                      CONFIG_HIGHEST_LEVEL     =-1)
end

    """
Global library options.

These are used to select which global option to set or get and are used in `git_libgit2_opts()`.
    """
    @enum(GIT_OPT,  GET_MWINDOW_SIZE         = 0,
                    SET_MWINDOW_SIZE         = 1,
                    GET_MWINDOW_MAPPED_LIMIT = 2,
                    SET_MWINDOW_MAPPED_LIMIT = 3,
                    GET_SEARCH_PATH          = 4,
                    SET_SEARCH_PATH          = 5,
                    SET_CACHE_OBJECT_LIMIT   = 6,
                    SET_CACHE_MAX_SIZE       = 7,
                    ENABLE_CACHING           = 8,
                    GET_CACHED_MEMORY        = 9,
                    GET_TEMPLATE_PATH        = 10,
                    SET_TEMPLATE_PATH        = 11,
                    SET_SSL_CERT_LOCATIONS   = 12)

    @enum(GIT_PROXY, PROXY_NONE,
                     PROXY_AUTO,
                     PROXY_SPECIFIED)

end

max_stars_repo_name: liveinjs/julia

max_issues_repo_name: liveinjs/julia

max_forks_repo_name: liveinjs/julia


text:
@testset "SimResult" begin
    @info "Testing SimResult"

    import ControlSystems: SimResult, ssrand

    y = randn(1,10)
    u = randn(1,10)
    x = randn(1,10)
    t = 1:10
    sys = ssrand(1, 1, 1, Ts=1)

    r = SimResult(y,t,x,u,sys)

    # test getindex
    @test r[1] == y
    @test r[2] == t
    @test r[3] == x
    @test r[4] == u
    @test r[5] == sys

    # test destructuring
    y2,t2,x2,u2 = r
    @test y2 === y
    @test t2 === t
    @test x2 === x
    @test u2 === u

    # test inferability of destructuring
    foo(r::SimResult) = (a,b,c,d) = r
    @inferred foo(r)

    # test properties
    @test r.nx == 1
    @test r.nu == 1
    @test r.ny == 1
end


text:
immutable GadflyPackage       <: PlottingPackage end
immutable ImmersePackage      <: PlottingPackage end
immutable PyPlotPackage       <: PlottingPackage end
immutable QwtPackage          <: PlottingPackage end
immutable UnicodePlotsPackage <: PlottingPackage end
immutable WinstonPackage      <: PlottingPackage end
immutable BokehPackage        <: PlottingPackage end
immutable PlotlyPackage       <: PlottingPackage end
immutable GRPackage           <: PlottingPackage end
immutable GLVisualizePackage  <: PlottingPackage end
immutable NoPackage           <: PlottingPackage end

typealias GadflyOrImmerse @compat(Union{GadflyPackage, ImmersePackage})

export
  gadfly,
  immerse,
  pyplot,
  qwt,
  unicodeplots,
  bokeh,
  plotly,
  gr,
  glvisualize
  # winston

gadfly(; kw...)        = (default(; kw...); backend(:gadfly))
immerse(; kw...)       = (default(; kw...); backend(:immerse))
pyplot(; kw...)        = (default(; kw...); backend(:pyplot))
qwt(; kw...)           = (default(; kw...); backend(:qwt))
unicodeplots(; kw...)  = (default(; kw...); backend(:unicodeplots))
bokeh(; kw...)         = (default(; kw...); backend(:bokeh))
plotly(; kw...)        = (default(; kw...); backend(:plotly))
gr(; kw...)            = (default(; kw...); backend(:gr))
glvisualize(; kw...)   = (default(; kw...); backend(:glvisualize))
# winston(; kw...)       = (default(; kw...); backend(:winston))

backend_name(::GadflyPackage)       = :gadfly
backend_name(::ImmersePackage)      = :immerse
backend_name(::PyPlotPackage)       = :pyplot
backend_name(::UnicodePlotsPackage) = :unicodeplots
backend_name(::QwtPackage)          = :qwt
backend_name(::BokehPackage)        = :bokeh
backend_name(::PlotlyPackage)       = :plotly
backend_name(::GRPackage)           = :gr
backend_name(::GLVisualizePackage)  = :glvisualize
backend_name(::NoPackage)           = :none

include("backends/supported.jl")

include("backends/qwt.jl")
include("backends/gadfly.jl")
include("backends/unicodeplots.jl")
include("backends/pyplot.jl")
include("backends/immerse.jl")
include("backends/winston.jl")

include("backends/web.jl")
include("backends/bokeh.jl")
include("backends/plotly.jl")
include("backends/gr.jl")

include("backends/glvisualize.jl")

# ---------------------------------------------------------

plot(pkg::PlottingPackage; kw...) = error("plot($pkg; kw...) is not implemented")
plot!(pkg::PlottingPackage, plt::Plot; kw...) = error("plot!($pkg, plt; kw...) is not implemented")
_update_plot(pkg::PlottingPackage, plt::Plot, d::Dict) = error("_update_plot($pkg, plt, d) is not implemented")
# Base.display(pkg::PlottingPackage, plt::Plot) = error("display($pkg, plt) is not implemented")

_update_plot_pos_size{P<:PlottingPackage}(plt::PlottingObject{P}, d::Dict) = nothing #error("_update_plot_pos_size(plt,d) is not implemented for $P")

subplot(pkg::PlottingPackage; kw...) = error("subplot($pkg; kw...) is not implemented")
subplot!(pkg::PlottingPackage, subplt::Subplot; kw...) = error("subplot!($pkg, subplt; kw...) is not implemented")
# Base.display(pkg::PlottingPackage, subplt::Subplot) = error("display($pkg, subplt) is not implemented")

# ---------------------------------------------------------

const BACKENDS = [:qwt, :gadfly, :unicodeplots, :pyplot, :immerse, :bokeh, :plotly, :gr]
const INITIALIZED_BACKENDS = Set{Symbol}()
backends() = BACKENDS

function backendInstance(sym::Symbol)
  sym == :qwt && return QwtPackage()
  sym == :gadfly && return GadflyPackage()
  sym == :unicodeplots && return UnicodePlotsPackage()
  sym == :pyplot && return PyPlotPackage()
  sym == :immerse && return ImmersePackage()
  sym == :winston && return WinstonPackage()
  sym == :bokeh && return BokehPackage()
  sym == :plotly && return PlotlyPackage()
  sym == :gr && return GRPackage()
  sym == :glvisualize && return GLVisualizePackage()
  sym == :none && return NoPackage()
  error("Unsupported backend $sym")
end

type CurrentBackend
  sym::Symbol
  pkg::PlottingPackage
end
CurrentBackend(sym::Symbol) = CurrentBackend(sym, backendInstance(sym))

# ---------------------------------------------------------

function pickDefaultBackend()
  for pkgstr in ("PyPlot", "Immerse", "Qwt", "Gadfly", "GR", "UnicodePlots", "Bokeh", "GLVisualize")
    if Pkg.installed(pkgstr) != nothing
      return backend(symbol(lowercase(pkgstr)))
    end
  end
  backend(:plotly)
end

# ---------------------------------------------------------

"""
Returns the current plotting package name.  Initializes package on first call.
"""
function backend()

  global CURRENT_BACKEND
  if CURRENT_BACKEND.sym == :none
    pickDefaultBackend()
  end

  currentBackendSymbol = CURRENT_BACKEND.sym
  if !(currentBackendSymbol in INITIALIZED_BACKENDS)

    # initialize
    println("[Plots.jl] Initializing backend: ", CURRENT_BACKEND.sym)
    if currentBackendSymbol == :qwt
      try
        @eval import Qwt
        @eval export Qwt
      catch err
        warn("Couldn't import Qwt.  Install it with: Pkg.clone(\"https://github.com/tbreloff/Qwt.jl.git\")\n  (Note: also requires pyqt and pyqwt).")
        rethrow(err)
      end

    elseif currentBackendSymbol == :gadfly
      try
        @eval import Gadfly, Compose
        @eval export Gadfly, Compose
        @eval include(joinpath(Pkg.dir("Plots"), "src", "backends", "gadfly_shapes.jl"))
      catch err
        warn("Couldn't import Gadfly.  Install it with: Pkg.add(\"Gadfly\").")
        rethrow(err)
      end

    elseif currentBackendSymbol == :unicodeplots
      try
        @eval import UnicodePlots
        @eval export UnicodePlots
      catch err
        warn("Couldn't import UnicodePlots.  Install it with: Pkg.add(\"UnicodePlots\").")
        rethrow(err)
      end

    elseif currentBackendSymbol == :pyplot
      try
        @eval import PyPlot
        @eval export PyPlot
        @eval const pycolors = PyPlot.pywrap(PyPlot.pyimport("matplotlib.colors"))
        @eval const pypath = PyPlot.pywrap(PyPlot.pyimport("matplotlib.path"))
        @eval const mplot3d = PyPlot.pywrap(PyPlot.pyimport("mpl_toolkits.mplot3d"))
        # @eval const pycolorbar = PyPlot.pywrap(PyPlot.pyimport("matplotlib.colorbar"))
        if !isa(Base.Multimedia.displays[end], Base.REPL.REPLDisplay)
          PyPlot.ioff()  # stops wierd behavior of displaying incomplete graphs in IJulia

          # # TODO: how the hell can I use PyQt4??
          # "pyqt4"=>:qt_pyqt4
          # PyPlot.backend[1] = "pyqt4"
          # PyPlot.gui[1] = :qt_pyqt4
          # PyPlot.switch_backend("Qt4Agg")

          # only turn on the gui if we want it
          if PyPlot.gui != :none
            PyPlot.pygui(true)
          end

        end
      catch err
        warn("Couldn't import PyPlot.  Install it with: Pkg.add(\"PyPlot\").")
        rethrow(err)
      end

    elseif currentBackendSymbol == :immerse
      try
        @eval import Immerse, Gadfly, Compose, Gtk
        @eval export Immerse, Gadfly, Compose, Gtk
        @eval include(joinpath(Pkg.dir("Plots"), "src", "backends", "gadfly_shapes.jl"))
      catch err
        # error("Couldn't import Immerse.  Install it with: Pkg.add(\"Immerse\").\n   Error: ", err)
        warn("Couldn't import Immerse.  Install it with: Pkg.add(\"Immerse\").")
        rethrow(err)
      end

    elseif currentBackendSymbol == :bokeh
      try
        @eval import Bokeh
        @eval export Bokeh
      catch err
        warn("Couldn't import Bokeh.  Install it with: Pkg.add(\"Bokeh\").")
        rethrow(err)
      end

    elseif currentBackendSymbol == :plotly
      try
        @eval begin
          import JSON
          JSON._print(io::IO, state::JSON.State, dt::Union{Date,DateTime}) = print(io, '"', dt, '"')

          ############################
          # borrowed from https://github.com/spencerlyon2/Plotlyjs.jl/blob/master/src/display.jl
          _js_path = joinpath(Pkg.dir("Plots"), "deps", "plotly-latest.min.js")

          # if we're in IJulia call setupnotebook to load js and css
          if isijulia()
              # the first script is some hack I needed to do in order for the notebook
              # to not complain about Plotly being undefined
              display("text/html", """

               """)
              # display("text/html", "
Plotly javascript loaded.
")
          end
          # end borrowing (thanks :)
          ###########################

          # try
          #   include(joinpath(Pkg.dir("Plots"), "src", "backends", "plotly_blink.jl"))
          # catch err
          #   warn("Error including Plotlyjs: $err\n  Note: Will fall back to built-in display.")
          # end

        end
      catch err
        warn("Couldn't setup Plotly")
        rethrow(err)
      end

    elseif currentBackendSymbol == :gr
      try
        @eval import GR
      catch err
        warn("Couldn't import GR.  Install it with: Pkg.add(\"GR\").")
      end

    elseif currentBackendSymbol == :glvisualize
      try
        @eval import GLVisualize
        @eval export GLVisualize
      catch err
        warn("Couldn't setup GLVisualize")
        rethrow(err)
      end

    elseif currentBackendSymbol == :winston
      warn("Winston support is deprecated and broken.  Try another backend: $BACKENDS")
      try
        @eval ENV["WINSTON_OUTPUT"] = "gtk"
        @eval import Winston, Gtk
        @eval export Winston, Gtk
      catch err
        warn("Couldn't import Winston.  Install it with: Pkg.add(\"Winston\").")
        rethrow(err)
      end

    else
      error("Unknown backend $currentBackendSymbol.  Choose from: $BACKENDS")
    end
    push!(INITIALIZED_BACKENDS, currentBackendSymbol)

  end
  CURRENT_BACKEND.pkg
end

"""
Set the plot backend.  Choose from:  :qwt, :gadfly, :unicodeplots, :immerse, :pyplot
"""
function backend(pkg::PlottingPackage)
  CURRENT_BACKEND.sym = backend_name(pkg)
  CURRENT_BACKEND.pkg = pkg
end

function backend(modname)

  # set the PlottingPackage
  if modname == :qwt
    CURRENT_BACKEND.pkg = QwtPackage()
  elseif modname == :gadfly
    CURRENT_BACKEND.pkg = GadflyPackage()
  elseif modname == :unicodeplots
    CURRENT_BACKEND.pkg = UnicodePlotsPackage()
  elseif modname == :pyplot
    CURRENT_BACKEND.pkg = PyPlotPackage()
  elseif modname == :immerse
    CURRENT_BACKEND.pkg = ImmersePackage()
  elseif modname == :winston
    CURRENT_BACKEND.pkg = WinstonPackage()
  elseif modname == :bokeh
    CURRENT_BACKEND.pkg = BokehPackage()
  elseif modname == :plotly
    CURRENT_BACKEND.pkg = PlotlyPackage()
  elseif modname == :gr
    CURRENT_BACKEND.pkg = GRPackage()
  elseif modname == :glvisualize
    CURRENT_BACKEND.pkg = GLVisualizePackage()
  else
    error("Unknown backend $modname.  Choose from: $BACKENDS")
  end

  # update the symbol
  CURRENT_BACKEND.sym = modname
  # println("[Plots.jl] Switched to backend: ", modname)

  # return the package
  CURRENT_BACKEND.pkg
end