text:
using DataInterpolations, Test
using StableRNGs

@testset "Linear Interpolation" begin
    u = 2.0collect(1:10)
    t = 1.0collect(1:10)
    A = LinearInterpolation(u,t)

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    @test A(0)        == 0.0
    @test A(5.5)      == 11.0
    @test A(11)       == 22

    u = vcat(2.0collect(1:10)', 3.0collect(1:10)')
    A = LinearInterpolation(u,t)

    for (_t, _u) in zip(t, eachcol(u))
        @test A(_t) == _u
    end
    @test A(0)        == [0.0, 0.0]
    @test A(5.5)      == [11.0, 16.5]
    @test A(11)       == [22, 33]
end

@testset "Quadratic Interpolation" begin
    u = [1.0, 4.0, 9.0, 16.0]
    t = [1.0, 2.0, 3.0, 4.0]
    A = QuadraticInterpolation(u,t)

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    @test A(0.0) == 0.0
    @test A(1.5) == 2.25
    @test A(2.5) == 6.25
    @test A(3.5) == 12.25
    @test A(5.0) == 25

    u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
    A = QuadraticInterpolation(u,t)

  for (_t, _u) in zip(t, eachcol(u))
        @test A(_t) == _u
    end
    @test A(0.0) == [0.0  ,  0.0 ]
    @test A(1.5) == [2.25 ,  2.25]
    @test A(2.5) == [6.25 ,  6.25]
    @test A(3.5) == [12.25, 12.25]
    @test A(5.0) == [25.0 , 25.0 ]
end

@testset "Lagrange Interpolation" begin
    u = [1.0, 4.0, 9.0]
    t = [1.0, 2.0, 3.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == 4.0
    @test A(1.5) == 2.25

    u = [1.0, 8.0, 27.0, 64.0]
    t = [1.0, 2.0, 3.0, 4.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == 8.0
    @test A(1.5) ≈ 3.375
    @test A(3.5) ≈ 42.875

    u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == [4.0,4.0]
    @test A(1.5) ≈ [2.25,2.25]
    @test A(3.5) ≈ [12.25,12.25]
end

@testset "Akima Interpolation" begin
    u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
    t = collect(0.0:10.0)
    A = AkimaInterpolation(u, t)

    @test A(0.0) ≈ 0.0
    @test A(0.5) ≈ 1.375
    @test A(1.0) ≈ 2.0
    @test A(1.5) ≈ 1.5
    @test A(2.5) ≈ 1.953125
    @test A(3.5) ≈ 2.484375
    @test A(4.5) ≈ 4.1363636363636366866103344
    @test A(5.1) ≈ 5.9803623910336236590978842
    @test A(6.5) ≈ 5.5067291516462386624652936
    @test A(7.2) ≈ 5.2031367459745245795943447
    @test A(8.6) ≈ 4.1796554159017080820603951
    @test A(9.9) ≈ 3.4110386597938129327189927
    @test A(10.0) ≈ 3.0
end

@testset "ConstantInterpolation" begin

    t = [1.0, 2.0, 3.0, 4.0]

    @testset "Vector case" for u in
        [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[2]
        @test A(2.0) == u[2]
        @test A(2.5) == u[3]
        @test A(3.0) == u[3]
        @test A(3.5) == u[1]
        @test A(4.0) == u[1]
        @test A(4.5) == u[1]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[1]
        @test A(2.0) == u[2]
        @test A(2.5) == u[2]
        @test A(3.0) == u[3]
        @test A(3.5) == u[3]
        @test A(4.0) == u[1]
        @test A(4.5) == u[1]
    end

    @testset "Matrix case" for u in
        [[1.0 2.0 0.0 1.0; 1.0 2.0 0.0 1.0], ["B" "C" "A" "B"; "B" "C" "A" "B"]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[:,1]
        @test A(1.0) == u[:,1]
        @test A(1.5) == u[:,2]
        @test A(2.0) == u[:,2]
        @test A(2.5) == u[:,3]
        @test A(3.0) == u[:,3]
        @test A(3.5) == u[:,1]
        @test A(4.0) == u[:,1]
        @test A(4.5) == u[:,1]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[:,1]
        @test A(1.0) == u[:,1]
        @test A(1.5) == u[:,1]
        @test A(2.0) == u[:,2]
        @test A(2.5) == u[:,2]
        @test A(3.0) == u[:,3]
        @test A(3.5) == u[:,3]
        @test A(4.0) == u[:,1]
        @test A(4.5) == u[:,1]
    end
end

@testset "QuadraticSpline Interpolation" begin
    u = [0.0, 1.0, 3.0]
    t = [-1.0, 0.0, 1.0]

    A = QuadraticSpline(u,t)

    # Solution
    P₁ = x -> (x + 1)^2 # for x ∈ [-1, 0]
    P₂ = x -> 2*x + 1   # for x ∈ [ 0, 1]

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    @test A(-2.0) == P₁(-2.0)
    @test A(-0.5) == P₁(-0.5)
    @test A(0.7)  == P₂( 0.7)
    @test A(2.0)  == P₂( 2.0)
end

@testset "CubicSpline Interpolation" begin
    u = [0.0, 1.0, 3.0]
    t = [-1.0, 0.0, 1.0]

    A = CubicSpline(u,t)

    # Solution
    P₁ = x -> 1 + 1.5x + x^2 + 0.5x^3 # for x ∈ [-1.0, 0.0]
    P₂ = x -> 1 + 1.5x + x^2 - 0.5x^3 # for x ∈ [0.0, 1.0]

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    for x in (-1.5, -0.5, -0.7)
        @test A(x) ≈ P₁(x)
    end
    for x in (0.3, 0.5, 1.5)
        @test A(x) ≈ P₂(x)
    end
end

# BSpline Interpolation and Approximation
t = [0,62.25,109.66,162.66,205.8,252.3]
u = [14.7,11.51,10.41,14.95,12.24,11.22]

A = BSplineInterpolation(u,t,2,:Uniform,:Uniform)

@test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845]
@test [A(190.0), A(225.0)] == [14.07428439395079, 11.057784141519251]
@test [A(t[1]), A(t[end])] == [u[1], u[end]]

A = BSplineInterpolation(u,t,2,:ArcLen,:Average)

@test [A(25.0), A(80.0)] == [13.363814458968486, 10.685201117692609]
@test [A(190.0), A(225.0)] == [13.437481084762863, 11.367034741256463]
@test [A(t[1]), A(t[end])] == [u[1], u[end]]

A = BSplineApprox(u,t,2,4,:Uniform,:Uniform)

@test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178]
@test [A(190.0), A(225.0)] ≈ [13.851245975109263, 12.963685868886575]
@test [A(t[1]), A(t[end])] ≈ [u[1], u[end]]

# Curvefit Interpolation
rng = StableRNG(12345)
model(x, p) = @. p[1]/(1 + exp(x - p[2]))
t = range(-10, stop=10, length=40)
u = model(t, [1.0, 2.0]) + 0.01*randn(rng, length(t))
p0 = [0.5, 0.5]

A = Curvefit(u, t, model, p0, LBFGS())

ts = [-7.0, -2.0, 0.0, 2.5, 5.0]
vs = [1.0013468217936277, 0.9836755196317837, 0.8833959853995836, 0.3810348276782708, 0.048062978598861855]
us = A.(ts)

@test vs ≈ us

# missing values handling tests
u = [1.0, 4.0, 9.0, 16.0, 25.0, missing, missing]
t = [1.0, 2.0, 3.0, 4.0, missing, 6.0, missing]
A = QuadraticInterpolation(u,t)

@test A(2.0) == 4.0
@test A(1.5) == 2.25
@test A(3.5) == 12.25
@test A(2.5) == 6.25

u = copy(hcat(u, u)')
A = QuadraticInterpolation(u,t)

@test A(2.0) == [4.0, 4.0]
@test A(1.5) == [2.25, 2.25]
@test A(3.5) == [12.25,12.25]
@test A(2.5) == [6.25, 6.25]


text:
"""
    get_target(x::Parser)
Return the [`TargetInfo`](@ref) in the [`Parser`](@ref).
"""
get_target(x::Parser) = getTargetInfo(x)

"""
    getSema(x::Parser)
Return the [`Sema`](@ref) in the [`Parser`](@ref).
"""
getSema(x::Parser) = getActions(x)

ext: jl


text:
using ArgCheck

# arithmetic comparisons between scalar tensors and scalars
Base.isless(t::Tensor{T,0}, val::S) where {T,S<:Number} = t.data[1] < val
Base.isless(val::S, t::Tensor{T,0}) where {T,S<:Number} = val < t.data[1]
Base.:(==)(t::Tensor{T,0}, val::S) where {T,S<:Number} = t.data[1] == val
Base.:(==)(val::S, t::Tensor{T,0}) where {T,S<:Number} = val == t.data[1]

# arithmetic comparisons between scalar tensors
Base.isless(t1::Tensor{T,0}, t2::Tensor{S,0}) where {T,S} = t1.data[1] < t2.data[1]

# arithmetic operations between scalar tensors and scalars
Base.:+(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1] + val)
Base.:+(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val + t.data[1])
Base.:-(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1] - val)
Base.:-(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val - t.data[1])
Base.:*(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1] * val)
Base.:*(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val * t.data[1])
Base.:/(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1] / val)
Base.:/(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val / t.data[1])
Base.:÷(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1] ÷ val)
Base.:÷(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val ÷ t.data[1])
Base.:\(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1] \ val)
Base.:\(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val \ t.data[1])
Base.:^(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1]^val)
Base.:^(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val^t.data[1])
Base.:%(t::Tensor{T,0}, val::S) where {T,S<:Number} = Tensor(t.data[1] % val)
Base.:%(val::S, t::Tensor{T,0}) where {T,S<:Number} = Tensor(val % t.data[1])

# arithmetic operations between scalar tensors
Base.:*(t1::Tensor{T,0}, t2::Tensor{S,0}) where {T,S} = Tensor(t1.data[1] * t2.data[1])
Base.:/(t1::Tensor{T,0}, t2::Tensor{S,0}) where {T,S} = Tensor(t1.data[1] / t2.data[1])
Base.:÷(t1::Tensor{T,0}, t2::Tensor{S,0}) where {T,S} = Tensor(t1.data[1] ÷ t2.data[1])
Base.:\(t1::Tensor{T,0}, t2::Tensor{S,0}) where {T,S} = Tensor(t1.data[1] \ t2.data[1])
Base.:^(t1::Tensor{T,0}, t2::Tensor{S,0}) where {T,S} = Tensor(t1.data[1]^t2.data[1])
Base.:%(t1::Tensor{T,0}, t2::Tensor{S,0}) where {T,S} = Tensor(t1.data[1] % t2.data[1])

# arithmetic operations for arbitrary tensors
Base.:+(t::Tensor) = t
function Base.:+(t1::Tensor, t2::Tensor)
    @argcheck size(t1) == size(t2) DimensionMismatch
    return Tensor(t1.data + t2.data)
end

Base.:-(t::Tensor) = Tensor(-t.data)
function Base.:-(t1::Tensor, t2::Tensor)
    @argcheck size(t1) == size(t2) DimensionMismatch
    return Tensor(t1.data - t2.data)
end

# NOTE: Scalar tensors behaving like scalars for these operations might
# be removed in future versions
Base.:*(t::Tensor, val::S) where {S<:Number} = Tensor(t.data * val)
Base.:*(val::S, t::Tensor) where {S<:Number} = Tensor(val * t.data)
Base.:*(t1::Tensor, t2::Tensor{T,0}) where {T} = Tensor(t1.data * t2.data[1])
Base.:*(t1::Tensor{T,0}, t2::Tensor) where {T} = Tensor(t1.data[1] * t2.data)

Base.:/(t::Tensor, val::S) where {S<:Number} = Tensor(t.data / val)
Base.:/(t1::Tensor, t2::Tensor{T,0}) where {T} = Tensor(t1.data / t2.data[1])

# allow typical array operations like sum or maximum for tensors
function Base.sum(f::Union{Function,Type}, t::Tensor; kwargs...)
    return Tensor(sum(f, t.data; kwargs...))
end
Base.sum(t::Tensor; kwargs...) = sum(identity, t; kwargs...)

function Base.prod(f::Union{Function,Type}, t::Tensor; kwargs...)
    return Tensor(prod(f, t.data; kwargs...))
end
Base.prod(t::Tensor; kwargs...) = prod(identity, t; kwargs...)

function Base.maximum(f::Union{Function,Type}, t::Tensor; kwargs...)
    return Tensor(maximum(f, t.data; kwargs...))
end
Base.maximum(t::Tensor; kwargs...) = maximum(identity, t; kwargs...)

function Base.minimum(f::Union{Function,Type}, t::Tensor; kwargs...)
    return Tensor(minimum(f, t.data; kwargs...))
end
Base.minimum(t::Tensor; kwargs...) = minimum(identity, t; kwargs...)

function Base.extrema(f::Union{Function,Type}, t::Tensor; kwargs...)
    return (minimum(f, t; kwargs...), maximum(f, t; kwargs...))
end
Base.extrema(t::Tensor; kwargs...) = extrema(identity, t; kwargs...)

Base.argmax(t::Tensor; kwargs...) = argmax(t.data; kwargs...)
Base.argmin(t::Tensor; kwargs...) = argmin(t.data; kwargs...)

Base.findmax(t::Tensor; kwargs...) = (maximum(t; kwargs...), argmax(t; kwargs...))
Base.findmin(t::Tensor; kwargs...) = (minimum(t; kwargs...), argmin(t; kwargs...))

Base.conj(t::Tensor) = Tensor(conj(t.data))


text:
"""
Check what average L1 distance is for two random vertices
in a grid graph of variable size.
"""

using LightGraphs
using QuNet
using Statistics
using Plots

function make_userpairs(path_length::Int)
    p1 = rand(1:path_length)
    p2 = rand(1:path_length)
    while p2 == p1
        p2 = rand(1:path_length)
    end
    return p1, p2
end

# Params
numSamples = 10000
minDim = 2
maxDim = 40
makeNewData = true

# MAIN
if makeNewData == true
    lengths = [[] for i in minDim:maxDim]
    for (dimidx, dim) in enumerate(minDim:maxDim)
        for sample in 1:numSamples
            g = LightGraphs.path_graph(dim)
            p1, p2 = make_userpairs(dim)
            shortestPath = a_star(g, p1, p2)
            len = length(shortestPath)
            push!(lengths[dimidx], len)
        end
    end

    L1 = []
    for l in lengths
        path = mean(l)
        push!(L1, path)
    end
end

HL1 = []

function HudsonL1(n::Int)::Float64
    # return (2/3 * n*(n^2 -1) + n^2*(n+1))/(2*n^2 - 1)
    return 1/3 * (n+1)
end

for n in minDim:maxDim
    push!(HL1, HudsonL1(n))
end

println("Experimental data")
println(L1)
println("Analytic data")
println(HL1)
x = collect(minDim:maxDim)
plot(x, L1, label = "L1")
plot!(x, HL1, label = "Hudson's L1")
xlabel!("Path size")
ylabel!("Average path length")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
# Euclidean Manifold definitions and arithmetic

struct Euclid2 <: Manifold
  dof::Int
  addop::Function
  diffop::Function
  getMu
  getLambda
  domain::Tuple{Tuple{Float64,Float64},Tuple{Float64, Float64}}
  Euclid2() = new(2, +, -, KDE.getEuclidMu, KDE.getEuclidLambda, ((-Inf,Inf),(-Inf,Inf)))
end

# ??
struct EuclideanManifold <: Manifold
end

get2DLambda(Lambdas::Vector{Float64})::Float64 = sum(Lambdas)

function *(PP::Vector{MB{EuclideanManifold,BallTreeDensity}})
  bds = Vector{BallTreeDensity}(undef, length(PP))
  for p in PP
    bds[i] = p.belief
  end
  *(bds)
end

#


text:
using HelloCExecutable
using Test

@testset "hello_world" begin
    @test hello_world() == "Hello, Goma!"
end


text:
@testset "HTTP.serve" begin

# test kill switch
server = HTTP.Server()
tsk = @async HTTP.serve(server)
sleep(1.0)
put!(server.in, HTTP.Nitrogen.KILL)
sleep(0.1)
@test istaskdone(tsk)

# test http vs. https

# echo response
serverlog = HTTP.FIFOBuffer()
server = HTTP.Server((req, rep) -> HTTP.Response(String(req)), serverlog)
tsk = @async HTTP.serve(server)
sleep(1.0)
r = HTTP.get("http://127.0.0.1:8081/"; readtimeout=30)

@test HTTP.status(r) == 200
@test String(take!(r)) == ""

print(String(read(serverlog)))

# invalid HTTP
sleep(2.0)
tcp = connect(ip"127.0.0.1", 8081)
write(tcp, "GET / HTP/1.1\r\n\r\n")
sleep(2.0)
log = String(read(serverlog))

print(log)
@test contains(log, "invalid HTTP version")

# bad method
sleep(2.0)
tcp = connect(ip"127.0.0.1", 8081)
write(tcp, "BADMETHOD / HTTP/1.1\r\n\r\n")
sleep(2.0)

log = String(read(serverlog))

print(log)
@test contains(log, "invalid HTTP method")

# Expect: 100-continue
sleep(2.0)
tcp = connect(ip"127.0.0.1", 8081)
write(tcp, "POST / HTTP/1.1\r\nContent-Length: 15\r\nExpect: 100-continue\r\n\r\n")
sleep(2.0)

log = String(read(serverlog))

@test contains(log, "sending 100 Continue response to get request body")
client = String(readavailable(tcp))
@test client == "HTTP/1.1 100 Continue\r\n\r\n"

write(tcp, "Body of Request")
sleep(2.0)
log = String(read(serverlog))
client = String(readavailable(tcp))

print(client)
@test contains(client, "HTTP/1.1 200 OK\r\n")
@test contains(client, "Connection: keep-alive\r\n")
@test contains(client, "Content-Length: 15\r\n")
@test contains(client, "\r\n\r\nBody of Request")

put!(server.in, HTTP.Nitrogen.KILL)

serverlog = HTTP.FIFOBuffer()
server = HTTP.Server((req, rep) -> begin
    io = HTTP.FIFOBuffer()
    @async begin
        i = 0
        while true
            println(io, "data: $(now())\n")
            sleep(2)
            i += 1
            i > 5 && break
        end
        close(io)
    end
    r = HTTP.Response(200, Dict(
        "Content-Type" => "text/event-stream",
        "Cache-Control" => "no-cache",
        "Connection" => "keep-alive"), io)
end, serverlog)

tsk = @async HTTP.serve(server, IPv4(0,0,0,0), 8082)
sleep(5.0)
r = HTTP.get("http://localhost:8082/"; readtimeout=30, verbose=true)
log = String(read(serverlog))
println(log)
@test length(String(r)) > 175

# test readtimeout, before sending anything and then mid-request

# header overflow

# upgrade request

# handler throw error

# keep-alive vs. close: issue #81
tsk = @async HTTP.serve(HTTP.Server((req, res) -> Response("Hello\n"), STDOUT), ip"127.0.0.1", 8083)
sleep(2.0)
r = HTTP.request(HTTP.Request(major=1, minor=0, uri=HTTP.URI("http://127.0.0.1:8083/"), headers=Dict("Host"=>"127.0.0.1:8083")))
@test HTTP.status(r) == 200
@test HTTP.headers(r)["Connection"] == "close"

# body too big

# other bad requests

end # @testset


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Glibmm"
version = v"2.68.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://download.gnome.org/sources/glibmm/$(version.major).$(version.minor)/glibmm-$(version).tar.xz",
                  "6664e27c9a9cca81c29e35687f49f2e0d173a2fc9e98c3428311f707db532f8c")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/glibmm*/
apk add doxygen graphviz libxslt
meson --cross-file=${MESON_TARGET_TOOLCHAIN} --libdir lib builddir .
cd builddir/
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct(["libglibmm-$(version.major).$(version.minor)", "libglibmm-$(version.major)"], :libglibmm),
    LibraryProduct(["libglibmm_generate_extra_defs-$(version.major).$(version.minor)", "libglibmm_generate_extra_defs-$(version.major)"], :libglibmm_generate_extra_defs),
    LibraryProduct(["libgiomm-$(version.major).$(version.minor)", "libgiomm-$(version.major)"], :libgiomm)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Glib_jll", uuid="7746bdde-850d-59dc-9ae8-88ece973131d"); compat=string(version))
    Dependency(PackageSpec(name="libsigcpp_jll", uuid="8855df87-3ca1-5a3d-a68b-4f0ddc198ba8"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#####
##### Read Stata DTA files using a row-wise interface, via Tables
#####

module StataDTAFiles

using ArgCheck: @argcheck
import Dates
using DocStringExtensions: SIGNATURES, TYPEDEF
using Parameters: @unpack
using StrFs: StrF
import Tables

import Base: read, seek, iterate, length, open, close, eltype, show

export DTAFile, elapsed_days

####
#### types for byteorder handling and IO wrapper
####

abstract type ByteOrder end

struct MSF <: ByteOrder end     # big-endian

struct LSF <: ByteOrder end     # little-endian

struct ByteOrderIO{B <: ByteOrder, T <: IO} <: IO
    byteorder::B
    io::T
end

read(io::ByteOrderIO, ::Type{UInt8}) = read(io.io, UInt8)

seek(io::ByteOrderIO, pos) = seek(io.io, pos)

####
#### tag verification
####

function verifytag(io::IO, tag::AbstractArray{UInt8}, closing::Bool = false)
    read(io, UInt8) == UInt8('<') || error("First byte is not a '<'.")
    if closing
        read(io, UInt8) == UInt8('/') || error("/ expected (closing tag).")
    end
    for c in tag
        read(io, UInt8) == c || error("not the expected tag.")
    end
    read(io, UInt8) == UInt8('>') || error("Tag not closed with '>'.")
    nothing
end

verifytag(io::IO, tag::String, closing::Bool = false) = verifytag(io, codeunits(tag), closing)

function verifytag(f::Function, io::IO, tag)
    verifytag(io, tag)
    result = f(io)
    verifytag(io, tag, true)
    result
end

####
#### reading primitives
####

readfixedstring(io::IO, nb) = String(read!(io, Vector{UInt8}(undef, nb)))

function readbyteorder(io::IO)
    verifytag(io, "byteorder") do io
        order = readfixedstring(io, 3)
        if order == "MSF"
            MSF()
        elseif order == "LSF"
            LSF()
        else
            error("unknown byte order $(order)")
        end
    end
end

"""
Number types is Stata that correspond to native Julia types, and are denoted by the latter.
"""
const DATANUMTYPES = Union{Int8, Int16, Int32, Float32, Float64}

"""
Number types which are read by `readnum`, in addition to `DATANUMTYPES`.
"""
const EXTRANUMTYPES = Union{Int64, UInt8, UInt16, UInt32, UInt64}

"""
All number types read by `readnum`.
"""
const READNUMTYPES = Union{DATANUMTYPES, EXTRANUMTYPES}

readnum(boio::ByteOrderIO{MSF}, T::Type{<:READNUMTYPES}) = ntoh(read(boio.io, T))

readnum(boio::ByteOrderIO{LSF}, T::Type{<:READNUMTYPES}) = ltoh(read(boio.io, T))

"""
$(SIGNATURES)

Read chars into a buffer of length `len`, find the terminating "\0" (if any) and truncate,
returning a string.
"""
function readchompedstring(boio::ByteOrderIO, len::Integer)
    buffer = Vector{UInt8}(undef, len)
    read!(boio.io, buffer)
    numchars = findfirst(isequal(0x00), buffer)
    if numchars ≢ nothing
        resize!(buffer, numchars - 1)
    end
    String(buffer)
end

"""
$(SIGNATURES)

Read length (of type T), then read and chomp the string.
"""
function readchompedstring(boio::ByteOrderIO, T::Type{<:Integer})
    len = Int(readnum(boio, T))
    readchompedstring(boio, len)
end

####
#### header
####

"""
$(TYPEDEF)

DTA file header (without the byte order, which is encoded in the corresponding
`ByteOrderIO`.
"""
struct DTAHeader
    release::Int
    variables::Int
    observations::Int
    label::String
    timestamp::Dates.DateTime
end

"Date format of Stata file timestamps."
const TIMESTAMPFMT = Dates.DateFormat("d u y H:M") # eg 04 Jul 2032 04:23

function read_header(io::IO)
    verifytag(io, "header") do io
        @assert verifytag(io -> readfixedstring(io, 3), io, "release") == "118"
        byteorder = readbyteorder(io)
        boio = ByteOrderIO(byteorder, io)
        K = verifytag(boio -> readnum(boio, Int16), boio, "K")
        N = verifytag(boio -> readnum(boio, Int64), boio, "N")
        label = verifytag(boio -> readchompedstring(boio, Int16), boio, "label")
        timestamp_str = strip(verifytag(boio -> readchompedstring(boio, Int8), boio, "timestamp"))
        timestamp = Dates.DateTime(timestamp_str, TIMESTAMPFMT)
        DTAHeader(118, Int(K), Int(N), label, timestamp), boio
    end
end

####
#### map
####

struct DTAMap
    stata_data_open::Int64
    map::Int64
    variable_types::Int64
    varnames::Int64
    sortlist::Int64
    formats::Int64
    value_label_names::Int64
    variable_labels::Int64
    characteristics::Int64
    data::Int64
    strls::Int64
    value_labels::Int64
    stata_data_close::Int64
    eof::Int64
end

function read_map(boio::ByteOrderIO)
    verifytag(boio, "map") do boio
        map = DTAMap([readnum(boio, Int64) for _ in 1:14]...)
        @assert map.stata_data_open == 0
        map
    end
end

####
#### types
####

"""
Maximum length of `str#` (aka `strfs`) strings in Stata DTA files.
"""
const STRFSMAXLEN = 2045

"""
$(SIGNATURES)

Map the numerical type code in a DTA file to a Julia type.

Note that numeric types use a sentinel for missing values, and thus are decoded as
`Union{Missing, T}` for the corresponding type `T`.
"""
function decode_variable_type(code::UInt16)
    if 1 ≤ code ≤ STRFSMAXLEN
        StrF{Int(code)}
    elseif code == 32768
        String
    elseif 65526 ≤ code ≤ 65530
        Union{Missing, (Float64, Float32, Int32, Int16, Int8)[code - 65525]}
    else
        error("unrecognized variable type code $(code)")
    end
end

"""
$(SIGNATURES)

Read variable types from a DTA file. Return a tuple of Julia types.
"""
function read_variable_types(boio::ByteOrderIO, header::DTAHeader, map::DTAMap)
    seek(boio, map.variable_types)
    verifytag(boio, "variable_types") do boio
        ((decode_variable_type(readnum(boio, UInt16)) for _ in 1:header.variables)..., )
    end
end

####
#### metadata
####

"""
$(SIGNATURES)

Read variable names, returning a tuple of `Symbol`s.
"""
function read_variable_names(boio::ByteOrderIO, header::DTAHeader, map::DTAMap)
    seek(boio, map.varnames)
    verifytag(boio, "varnames") do boio
        ntuple(_ -> Symbol(readchompedstring(boio, 129)), header.variables)
    end
end

"""
$(SIGNATURES)

Read the sortlist, a vector if integers which contains column indexes for nested sorting.
"""
function read_sortlist(boio::ByteOrderIO, header::DTAHeader, map::DTAMap)
    seek(boio, map.sortlist)
    verifytag(boio, "sortlist") do boio
        sortlist = [readnum(boio, Int16) for _ in 1:(header.variables + 1)]
        terminator = findfirst(iszero, sortlist)
        @assert terminator ≢ nothing
        sortlist[1:(terminator-1)]
    end
end

"""
$(SIGNATURES)

Read the format strings, returning a `Vector{String}`.
"""
function read_formats(boio::ByteOrderIO, header::DTAHeader, map::DTAMap)
    seek(boio, map.formats)
    verifytag(boio, "formats") do boio
        [readchompedstring(boio, 57) for _ in 1:header.variables]
    end
end

####
#### read data
####

const MAXINT8 = Int8(0x64)
const MAXINT16 = Int16(0x7fe4)
const MAXINT32 = Int32(0x7fffffe5)
const MAXFLOAT32 = Float32(0x1.fffffep126)
const MAXFLOAT64 = Float64(0x1.fffffffffffffp1022)

decode_missing(x::Int8) = x ≥ MAXINT8 ? missing : x
decode_missing(x::Int16) = x ≥ MAXINT16 ? missing : x
decode_missing(x::Int32) = x ≥ MAXINT32 ? missing : x
decode_missing(x::Float32) = x ≥ MAXFLOAT32 ? missing : x
decode_missing(x::Float64) = x ≥ MAXFLOAT64 ? missing : x

readfield(boio::ByteOrderIO, ::Type{Union{Missing,T}}) where T =
    decode_missing(read(boio, T))

readfield(boio::ByteOrderIO, ::Type{StrF{N}}) where N = read(boio, StrF{N})

_readrow(boio::ByteOrderIO, ::Type{Tuple{}}) = ()

_readrow(boio::ByteOrderIO, ::Type{T}) where T <: Tuple =
    (readfield(boio, Base.tuple_type_head(T)), _readrow(boio, Base.tuple_type_tail(T))...)

# NOTE: type assertion below is needed because of
# https://github.com/JuliaLang/julia/issues/29970
# TODO: revisit once that improves
readrow(boio::ByteOrderIO, ::Type{NamedTuple{N,T}}) where {N,T} =
    NamedTuple{N,T}(_readrow(boio, T))::NamedTuple{N,T}

####
#### API
####

"""
$(TYPEDEF)

Representation of a Stata DTA file for which everything except the data itself has been read.

The data (rows) can be read as `NamedTuple`s, using the iteration interface.
"""
struct DTAFile{T <: NamedTuple, B <: ByteOrderIO}
    boio::B
    header::DTAHeader
    map::DTAMap
    sortlist::Vector{Int16}
    formats::Vector{String}
end

function show(io::IO, dta::DTAFile{T}) where T
    @unpack header, sortlist, formats = dta
    @unpack release, variables, observations, label, timestamp = header
    COLORHEADER = :red
    COLORVAR = :blue
    COLORTYPE = :green
    print(io, "Stata DTA file $(release), ")
    printstyled(io, "$(variables) vars"; color = COLORHEADER)
    print(io, " in ")
    printstyled(io, "$(observations) rows"; color = COLORHEADER)
    println(io, ", ", timestamp)
    isempty(label) || println(io, "    label: ", label)
    varnames = fieldnames(T)
    if isempty(sortlist)
        println(io, "    not sorted")
    else
        print(io, "    sorted by ")
        printstyled(io, varnames[sortlist], "\n"; color = COLORVAR)
    end
    for (i, (varname, format)) in enumerate(zip(varnames, formats))
        i == 1 || println(io)
        print(io, "  ")
        printstyled(io, varname; color = COLORVAR)
        print(io, "::")
        printstyled(io, fieldtype(T, varname); color = COLORTYPE)
        print(io, " [", format, "]")
    end
end

open(::Type{DTAFile}, path::AbstractString) = open(DTAFile, open(path, "r"))

function open(::Type{DTAFile}, io::IO)
    seekstart(io)
    verifytag(io, "stata_dta")
    # read header and map
    header, boio = read_header(io)
    map = read_map(boio)
    # read the rest using map
    variable_names = read_variable_names(boio, header, map)
    sortlist = read_sortlist(boio, header, map)
    formats = read_formats(boio, header, map)
    variable_types = read_variable_types(boio, header, map)
    DTAFile{NamedTuple{variable_names, Tuple{variable_types...}},
            typeof(boio)}(boio, header, map, sortlist, formats)
end

function open(f::Function, ::Type{DTAFile}, args...)
    dta = open(DTAFile, args...)
    try
        f(dta)
    finally
        close(dta)
    end
end

close(dta::DTAFile) = close(dta.boio.io)

####
#### row iteration interface
####

eltype(::Type{<: DTAFile{T}}) where T = T

length(dta::DTAFile) = dta.header.observations

function iterate(dta::DTAFile{T}, index = 1) where T
    @unpack boio = dta
    if index > dta.header.observations
        nothing
    else
        if index == 1
            seek(boio, dta.map.data)
            verifytag(boio, "data")
        end
        readrow(boio, T), index + 1
    end
end

####
#### Tables interface
####

Tables.istable(::Type{<:DTAFile}) = true

Tables.rowaccess(::Type{<:DTAFile}) = true

Tables.rows(dta::DTAFile) = dta

function Tables.schema(dta::DTAFile{NamedTuple{names, types}}) where {names, types}
    Tables.Schema(names, types)
end

####
#### date
####

"""
Start of the epoch, ie "day 0" for most date handling functions in Stata.

!!! note
    Don't use directly, see [`elapsed_days`](@ref).
"""
const EPOCH = Dates.Date(1960, 1, 1)

"""
$(SIGNATURES)

Convert a Stata "elapsed date" representation into `Date`.

Corresponds to the `%td` format in Stata.
"""
elapsed_days(Δ::Integer) = EPOCH + Dates.Day(Δ)

elapsed_days(Δ::Real) = elapsed_days(convert(Int, Δ))

end # module


text:
# stuff copy pasted from Sha.jl

# Common update and digest functions which work across SHA1 and SHA2

# update! takes in variable-length data, buffering it into blocklen()-sized pieces,
# calling transform!() when necessary to update the internal hash state.
function update!(context::MD5_CTX,data)
    T = typeof(context)
    # We need to do all our arithmetic in the proper bitwidth
    UIntXXX = typeof(context.bytecount)

    # Process as many complete blocks as possible
    len = convert(UIntXXX, length(data))
    data_idx = convert(UIntXXX, 0)
    usedspace = context.bytecount % blocklen(T)
    while len - data_idx + usedspace >= blocklen(T)
        # Fill up as much of the buffer as we can with the data given us
        for i in 1:(blocklen(T) - usedspace)
            context.buffer[usedspace + i] = data[data_idx + i]
        end

        transform!(context)
        context.bytecount += blocklen(T) - usedspace
        data_idx += blocklen(T) - usedspace
        usedspace = convert(UIntXXX, 0)
    end

    # There is less than a complete block left, but we need to save the leftovers into context.buffer:
    if len > data_idx
        for i = 1:(len - data_idx)
            context.buffer[usedspace + i] = data[data_idx + i]
        end
        context.bytecount += len - data_idx
    end
end

lrot(b,x,width) = ((x << b) | (x >> (width - b)))


text:
module HTTPDownloads

import Downloads: Downloads, Curl, Downloader

using Sockets: TCPSocket

using URIs
using HTTP

abstract type LibCurlLayer{Next <: HTTP.Layer} <: HTTP.Layer{Next} end

DOWNLOAD_LOCK = ReentrantLock()
DOWNLOADER = Ref{Union{Nothing,Downloader}}(nothing)

function http_easy_hook(easy, info)
    # Disable redirects - HTTP.jl will handle this.
    Curl.setopt(easy, Curl.CURLOPT_FOLLOWLOCATION, false)
end

function get_downloader()
    lock(DOWNLOAD_LOCK) do
        yield() # let other downloads finish
        downloader = DOWNLOADER[]
        if downloader === nothing
            # TODO: Install our own easy_handle here.
            downloader = DOWNLOADER[] = Downloader()
            downloader.easy_hook = http_easy_hook
        end
        downloader
    end
end

# Monkey-patch arg_read_size
# See https://github.com/JuliaLang/Downloads.jl/issues/142
Downloads.arg_read_size(io::Base.GenericIOBuffer) = bytesavailable(io)

function HTTP.request(::Type{LibCurlLayer{Next}}, url::URI, req, body;
            response_stream=nothing,
            iofunction=nothing,

            # # TimeoutLayer
            readtimeout=0,

            # # ConnectionLayer
            socket_type::Type=TCPSocket,

            verbose::Int=0,

            # List of all HTTP.jl keywords for the connection pool and stream
            # layers which we could/should consider emulating:

            # ConnectionPool / getconnection / newconnection
            #
            # * connection_limit::Int=default_connection_limit,
            #     https://curl.se/libcurl/c/CURLOPT_MAXCONNECTS.html
            # * pipeline_limit::Int=1
            #     # likely will be deprecated in HTTP.jl - ignore
            # * idle_timeout::Int=0,
            #     # Similar to the Downloader's grace parameter

            # ConnectionPool / sslconnection
            #
            # * require_ssl_verification=NetworkOptions.verify_host(host, "SSL"),
            #     https://curl.se/libcurl/c/CURLOPT_SSL_VERIFYPEER.html
            # * sslconfig::SSLConfig=nosslconfig,
            #     ???

            # * keepalive::Bool=false,
            #     https://curl.se/libcurl/c/CURLOPT_TCP_KEEPALIVE.html
            # * connect_timeout::Int=0,
            #     https://curl.se/libcurl/c/CURLOPT_CONNECTTIMEOUT.html
            # * readtimeout::Int=0,
            #     OK - use requests `timeout` keyword

            # * sslconfig::SSLConfig=nosslconfig,
            #     FIXME ca_roots ???

            # * proxy =
            #     https://curl.se/libcurl/c/CURLOPT_PROXY.html
            # * require_ssl_verification = NetworkOptions.verify_host(host),
            #     https://curl.se/libcurl/c/CURLOPT_SSL_VERIFYHOST.html
            # * reuse_limit =
            #     No equiv? Ignored for now

            # # StreamLayer
            # * reached_redirect_limit=false,
            #     FIXME ???
            # * response_stream=nothing,
            #     Done.
            # * iofunction=nothing,
            #     TODO: Figure out callbacks?
            # * verbose::Int=0,

            kw...) where Next

    if iofunction !== nothing || socket_type !== TCPSocket
        # Fallback to pure-Julia implementation in ConnectionPool.
        # This is required until we can figure out how to expose the iofunction
        # callback interface via libcurl
        HTTP.request(Next, url, req, body;
                     response_stream=response_stream,
                     iofunction=iofunction,
                     readtimeout=readtimeout,
                     socket_type=socket_type,
                     kw...)
    end

    input = req.body === HTTP.body_is_a_stream ? body :
            length(req.body) > 0 ? IOBuffer(req.body) : nothing

    output_buf = nothing
    output = response_stream
    if response_stream === nothing
        # HTTP.jl assumes you always want the body by default.
        output_buf = IOBuffer()
        if req.method != "HEAD"
            # Workaround https://github.com/JuliaLang/Downloads.jl/pull/131
            output = output_buf
        end
    end

    # Use Downloads.jl
    response = Downloads.request(string(url);
        downloader = get_downloader(),
        method = req.method,
        headers = req.headers,
        input = input,
        output = output,
        timeout = readtimeout <= 0 ? Inf : readtimeout,
        throw = false,
        verbose = verbose > 0)

    if response isa Downloads.RequestError
        # We couldn't even get a response. There's a huge list of possible
        # return codes in libcurl's curl.h.
        #
        # On the other hand, users of HTTP.jl largely seem to match these based
        # on whether they're "retryable", and expect to see HTTP.jl's error
        # types.
        #
        # For example, AWS.jl has the following logic:
        #
        #  # Base.IOError is needed because HTTP.jl can often have errors that aren't
        #  # caught and wrapped in an HTTP.IOError
        #  # https://github.com/JuliaWeb/HTTP.jl/issues/382
        #  @delay_retry if isa(e, Sockets.DNSError) ||
        #                  isa(e, HTTP.ParseError) ||
        #                  isa(e, HTTP.IOError) ||
        #                  isa(e, Base.IOError) ||
        #                  (isa(e, HTTP.StatusError) && _http_status(e) >= 500)
        #
        if response.code == Curl.CURLE_OPERATION_TIMEDOUT
            throw(HTTP.TimeoutRequest.ReadTimeoutError(readtimeout))
        else
            # Wrap other errors in HTTP.IOError. Downstream codes will likely
            # assume all such errors are retry-able
            #
            throw(HTTP.IOError(response, response.message))
        end
    end

    # Convert output back into HTTP.Response
    body = output_buf === nothing ?
           HTTP.MessageRequest.body_was_streamed :
           take!(output_buf)

    res = HTTP.Response(response.status, HTTP.mkheaders(response.headers);
                        body=body, request=req)
end

function has_http_layer(stack, layer)
    while true
        if stack === Union{}
            return false
        elseif stack <: layer
            return true
        end
        stack = stack.parameters[1]
    end
end

function set_downloads_backend(use_downloads::Bool=true)
    stack = HTTP.stack()
    if use_downloads && !has_http_layer(stack, LibCurlLayer)
        HTTP.insert_default!(ConnectionPoolLayer, LibCurlLayer)
    elseif !use_downloads && has_http_layer(stack, LibCurlLayer)
        HTTP.remove_default!(ConnectionPoolLayer, LibCurlLayer)
    end
    nothing
end

end


text:
print("Last updated: $(read(`date`, String))")

# Load environment
import Pkg; Pkg.activate("../../../")

# Import Libraries
using Turing
using Turing: Variational
using Distributions
using AbstractGPs, KernelFunctions
using PyPlot
using StatsFuns
using JSON3
using Flux
using ProgressBars
import Random
import LinearAlgebra

# Define a kernel.
function sekernel(alpha, rho, jitter=0)
  kernel = alpha^2 * transform(SEKernel(), invsqrt2/rho)
  jitter <= 0 || (kernel += jitter * EyeKernel())
  return kernel
end

function compute_f(kernel, X, eta, beta=0, jitter=0)
  K = kernelmatrix(kernel, X, obsdim=1)
  jitter <= 0 || (K += LinearAlgebra.I * jitter)
  return LinearAlgebra.cholesky(K).L * eta .+ beta
end

@model function GPClassify(y, X, jitter=1e-6)
    # Priors.
    alpha ~ LogNormal(0, 1)
    rho ~ LogNormal(0, 1)
    beta ~ Normal(0, 1)  # intercept.
    eta ~ filldist(Normal(0, 1), length(y))

    # Latent GP
    kernel = sekernel(alpha, rho)
    f = compute_f(sekernel(alpha, rho), X, eta, beta, jitter)

    # Sampling Distribution.
    y ~ Product(Bernoulli.(logistic.(f)))
end;

# For getting quantiles along array dimensions
quantiles(X, q; dims) = mapslices(x -> quantile(x, q), X, dims=dims)

# To extract parameters from ADVI model.
function make_extractor(m, q; nsamples=300)
    qsamples = rand(q, nsamples)
    _, sym2range = Variational.bijector(m; sym_to_ranges=Val(true))
    return sym -> qsamples[collect(sym2range[sym][1]), :]
end;

function plot_data(X, y; edgecolors=nothing)
    colors = ["blue", "red"]
    plt.scatter(X[:, 1], X[:, 2], color=colors[y .+ 1], edgecolors=edgecolors)
end;

# Function for plotting parameter posterior.
function plot_post(samples, name; bins=nothing, suffix="")
    param = samples[name]
    plt.hist(param, density=true, bins=bins)
    plt.xlabel(name)
    plt.ylabel("density")
    if suffix !== ""
        suffix = "($(suffix))"
    end
    plt.title("Histogram of $(name) $(suffix)");
end;

function makegrid(xmin, xmax, ymin, ymax, ngrid; return_all=false)
  ngrid = floor(Int, sqrt(ngrid))
  grid = [[i, j]
          for i in range(xmin, xmax, length=ngrid),
          j in range(ymin, ymax, length=ngrid)]
  out =  Matrix(hcat(vec(grid)...)')

  if return_all
    return out, grid
  else
    return out
  end
end;

function gp_predict(postsamps, X, y, nnew; eps=0, prob=true, jitter=1e-6)
    # Make grid of prediction locations.
    xmin = minimum(X, dims=1)
    xmax = maximum(X, dims=1)
    Xnew, grid = makegrid(xmin[1]-eps, xmax[1]+eps,
                          xmin[2]-eps, xmax[2]+eps,
                          nnew, return_all=true)

    # Posterior samples of model parameters.
    alpha = postsamps[:alpha]
    rho = postsamps[:rho]
    beta = postsamps[:beta]
    nsamps = length(alpha)
    eta = postsamps[:eta]

    # Predict classification probabilities.
    ps = [let
            kernel = sekernel(alpha[i], rho[i], jitter)
            gp = GP(beta[i], kernel)
            f = compute_f(kernel, X, eta[:, i], beta[i])
            pfx = posterior(gp(X'), f)
            logit_p = rand(pfx(Xnew'))
            prob ? logistic.(logit_p) : logit_p
         end for i in ProgressBar(1:nsamps)]

    return hcat(ps...), Xnew, grid
end;

function plot_kernel_params(post_samples, kernel_params; bins=nothing, figsize=(8,3))
    plt.figure(figsize=figsize)
    nparams = length(kernel_params)
    for i in 1:nparams
        param = kernel_params[i]
        plt.subplot(1, nparams, i)
        plot_post(post_samples, param, bins=bins);
    end
    plt.tight_layout()
end;

function plot_uq(postsamps, X, y, algo;
                 nnew=100, eps=0.2, figsize=(10, 8),
                 kernel_params=nothing, jitter=1e-6,
                 color_res=100, color_res_sd=100, return_stuff=false,
                 data_edgecolors=nothing, vmin_sd=0, vmax_sd=nothing, prob=true)
    nsamps = postsamps[:alpha]
    ps, Xnew, grid = gp_predict(postsamps, X, y, nnew, eps=eps, prob=prob, jitter=jitter)

    plt.figure(figsize=figsize)
    gs = [[g[i] for g in grid] for i in 1:2]

    # Plot mean function.
    plt.subplot(2, 2, 1)
    if prob
        plt.contourf(gs[1], gs[2], reshape(mean(ps, dims=2), size(gs[1])...),
                     color_res, cmap=plt.get_cmap("bwr"), vmin=0, vmax=1)
    else
        plt.contourf(gs[1], gs[2], reshape(mean(ps, dims=2), size(gs[1])...),
                     color_res, cmap=plt.get_cmap("bwr"))
    end
    plt.colorbar();
    plot_data(X, y, edgecolors=data_edgecolors);
    plt.xticks([]); plt.yticks([]);
    plt.title("Posterior Mean Function ($(algo))")

    # Plot function sd.
    plt.subplot(2, 2, 2)
    plt.contourf(gs[1], gs[2], reshape(std(ps, dims=2), size(gs[1])...),
                 color_res_sd, cmap=plt.get_cmap("Oranges"),
                 vmin=vmin_sd, vmax=vmax_sd)
    plt.colorbar();
    plot_data(X, y)
    plt.xticks([]);
    plt.yticks([]);
    plt.title("Posterior Standard Deviation Function ($algo)")

    if kernel_params != nothing
        plot_kernel_params(postsamps, kernel_params, bins=30)
    end

    if return_stuff
        return ps, Xnew, grid
    end
end;

# Read data.

# Path to data.
data_path = joinpath(@__DIR__, "../data/gp-classify-data-N50.json")

# Load data in JSON format.
data = let
    x = open(f -> read(f, String), data_path)
    JSON3.read(x)
end

# Store data.
X = [data["x1"] data["x2"]]
y = Int64.(data["y"])

# Plot data.
plot_data(X, y);

# Create model.
m =  GPClassify(Float64.(y), X);
kernel_params = [:alpha, :rho, :beta];

# Fit via ADVI.
Random.seed!(7)

# NOTE: ADVI(num_elbo_samples, max_iters)

# Compile
@time q = vi(m, ADVI(1, 1))

# RUN
@time q = vi(m, ADVI(1, 1000))

# Get posterior samples
extract_gp = make_extractor(m, q, nsamples=500)
advi_samples = Dict{Symbol, Any}(sym => vec(extract_gp(sym)) for sym in kernel_params)
advi_samples[:eta] = extract_gp(:eta);

plot_uq(advi_samples, X, y, "ADVI", eps=0.5,
        nnew=500, data_edgecolors="gold",
        kernel_params=kernel_params,
        color_res=range(0, 1, length=101),
        color_res_sd=range(0, .5, length=101));

# Fit via HMC.
Random.seed!(0)
burn = 500
nsamples = 500

# Compile
@time _ = sample(m, HMC(0.01, 1), 1)

# Run
@time hmc_chain = sample(m, HMC(0.05, 20), burn + nsamples);

# Get posterior samples
hmc_samples = Dict{Symbol, Any}([
    sym => vec(group(hmc_chain, sym).value.data)[end-nsamples+1:end]
for sym in kernel_params])
hmc_samples[:eta] = Matrix(group(hmc_chain, :eta).value.data[end-nsamples+1:end, :, 1]');

# Plot HMC
ps, Xnew, grid = plot_uq(hmc_samples, X, y, "HMC", eps=0.5,
                         nnew=500, data_edgecolors="gold",
                         kernel_params=kernel_params,
                         color_res=range(0, 1, length=101),
                         return_stuff=true);

plt.plot(get(hmc_chain, :lp)[1].data[burn+1:end]);
plt.xlabel("HMC iterations");
plt.ylabel("log unnormalized joint density");

# Fit via NUTS.
Random.seed!(0)
burn = 500
nsamples = 500

# Compile
@time _ = sample(m, NUTS(1, 0.8, max_depth=2), 2)

# Run
@time nuts_chain = sample(m, NUTS(burn, 0.8, max_depth=10), burn + nsamples);

# Get posterior samples
nuts_samples = Dict{Symbol, Any}([
    sym => vec(group(nuts_chain, sym).value.data)
for sym in kernel_params])
nuts_samples[:eta] = Matrix(group(nuts_chain, :eta).value.data[:, :, 1]');

# Plot NUTS
plot_uq(nuts_samples, X, y, "NUTS", eps=0.5,
        nnew=500, data_edgecolors="gold",
        kernel_params=kernel_params,
        color_res=range(0, 1, length=101),
        color_res_sd=range(0, .5, length=101),
        return_stuff=true);

plt.plot(get(nuts_chain, :lp)[1].data);
plt.xlabel("NUTS iterations");
plt.ylabel("log unnormalized joint density");

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
export
IndexedBinaryTree,
itree,
index,
ileaves,
itraversewhen,
ileft,
iright

const EmptyIndexedNode = -1

struct IndexedNode{T}
  left::Int
  right::Int
  next::Int # next can mean parent or next free node
  data::T
  IndexedNode{T}() where T = new{T}(EmptyIndexedNode, EmptyIndexedNode, EmptyIndexedNode)
  IndexedNode{T}(l,r) where T = new{T}(l, r, EmptyIndexedNode)
  IndexedNode{T}(l,r,n) where T = new{T}(l, r, n)
  IndexedNode{T}(l,r,n,d::T) where T = new{T}(l, r, n, d)
end

IndexedNode(a::T) where T = IndexedNode{T}(EmptyIndexedNode, EmptyIndexedNode, EmptyIndexedNode, a)
IndexedNode(l,r,n,d::T) where T = IndexedNode{T}(l, r, n, d)

mutable struct NodeVector{T}
  nodes::Vector{IndexedNode{T}}
  root::Int
  nextfree::Int
  last::Int
  #nextfreeleaf::Int  # can optionally reserve first block of indexes for leaves.  Enables direct indexing
end

isleaf(n::IndexedNode{T}) where T = EmptyIndexedNode == n.left &&
                                    EmptyIndexedNode == n.right

parent(n::IndexedNode{T}) where T = n.next
setparent(n::IndexedNode{T}, iparent) where T = IndexedNode{T}(n.left, n.right, iparent, n.data)
setnext(n::IndexedNode{T}, inext) where T = setparent(n, inext)

Base.length(nv::NodeVector{T}) where T = length(nv.nodes)

struct IndexedBinaryTree{T} <: BinaryTree{T}
  nodes::NodeVector{T}
  node::Int
end

IndexedBinaryTree{T}(d::T, ::Nothing, ::Nothing) where T = itree(d)
IndexedBinaryTree(d::T) where T = itree(d)
IndexedBinaryTree{T}(d::T) where T = itree(d)

isroot(t::IndexedBinaryTree{T}) where T = t.nodes.root == t.node

function NodeVector{T}(sz::Int) where T
  nodes = [IndexedNode{T}(EmptyIndexedNode, EmptyIndexedNode, i) for i in vcat(2:sz, EmptyIndexedNode)]
  NodeVector{T}(nodes, EmptyIndexedNode, 1, sz)
end

index(tree::IndexedBinaryTree{T}) where T = tree.node
node(tree::IndexedBinaryTree{T}) where T = tree.nodes[tree.node]

root(tree::IndexedBinaryTree{T}) where T = IndexedBinaryTree{T}(tree.nodes, tree.nodes.root)

Base.getindex(A::IndexedBinaryTree{T}, i::Number) where {T} = IndexedBinaryTree{T}(A.nodes, i)

function ileaves(t::IndexedBinaryTree{T}) where T
  L = Vector{IndexedNode{T}}(undef, 16)
  resize!(L, 0);
  ileaves(t.nodes, t.node, L);
  return L
end

function ileaves( n::NodeVector{T}, i::Int, L::Vector{IndexedNode{T}} ) where T
  nc = 0
  if n.nodes[i].left != EmptyIndexedNode
    nc = 1
    ileaves(n, n.nodes[i].left, L)
  end
  if n.nodes[i].right != EmptyIndexedNode
    nc += 1
    ileaves(n, n.nodes[i].right, L)
  end

  if nc == 0
    push!(L, n.nodes[i])
  end
end

function itraversewhen(when::Function, node::IndexedBinaryTree{T}, L::Vector{R},
                       pre::Union{Function, Nothing},
                       post::Union{Function, Nothing},
                       inorder::Union{Function, Nothing}) where {T, R}
  if !when(node)
    return
  end

  if pre != nothing
    push!(L, pre(node))
  end

  node.nodes[node.node].left != EmptyIndexedNode && itraversewhen(when, left(node), L, pre, post, inorder)

  if inorder != nothing
    push!(L, inorder(node))
  end

  node.nodes[node.node].right != EmptyIndexedNode && itraversewhen(when, right(node), L, pre, post, inorder)
  if post != nothing
    push!(L, post(node))
  end
end

function itraversewhen(when::Function, t::IndexedBinaryTree{T}, returnType::Type,
                       pre, post, inorder = nothing) where T
  L = Vector{returnType}(undef, length(t.nodes))
  resize!(L, 0);
  itraversewhen(when, t, L, pre, post, inorder);
  return L
end

function itraversewhen(when::Function, node::IndexedBinaryTree{T}, L::Vector{R},
                       pre::Function) where {T, R}
  if !when(node)
    return
  end

  push!(L, pre(node))

  node.nodes[node.node].left != EmptyIndexedNode && itraversewhen(when, left(node), L, pre)
  node.nodes[node.node].right != EmptyIndexedNode && itraversewhen(when, right(node), L, pre)
end

function itraversewhen(when::Function, t::IndexedBinaryTree{T}, returnType::Type,
                       pre::Function ) where T
  L = Vector{returnType}(undef, length(t.nodes))
  resize!(L, 0);
  itraversewhen(when, t, L, pre);
  return L
end

function resize( nv::NodeVector{T}, sz::Int ) where T
  oldsz = length(nv.nodes)
  if sz < oldsz
    return nv
  end

  if oldsz*2 > sz
    sz = oldsz*2
  end
  newnodes = [IndexedNode{T}(EmptyIndexedNode, EmptyIndexedNode, i)
              for i in vcat((oldsz+2):sz, EmptyIndexedNode)]
  nodes = vcat(nv.nodes, newnodes)
  lastn = nodes[nv.last]
  nodes[nv.last] = IndexedNode{T}(lastn.left, lastn.right, oldsz+1, lastn.data)
  NodeVector{T}(nodes, nv.root, oldsz+1, sz)
end

import Base.resize!
function resize!( nv::NodeVector{T}, sz::Int ) where T
  oldsz = length(nv.nodes)
  if sz < oldsz
    return nv
  end

  if oldsz*2 > sz
    sz = oldsz*2
  end
  resize!(nv.nodes, sz)
  newnodes = [IndexedNode{T}(EmptyIndexedNode, EmptyIndexedNode, i)
              for i in vcat((oldsz+2):sz, EmptyIndexedNode)]
  nv.nodes[oldsz+1:sz] = newnodes
  lastn = nv.nodes[nv.last]
  nv.nodes[nv.last] = IndexedNode{T}(lastn.left, lastn.right, oldsz+1, lastn.data)

  if nv.nextfree == EmptyIndexedNode
    nv.nextfree = oldsz+1
  end

  nv.last = sz
  return nv
end

"""
    reserve_leaves!(nv, n)

For internal use.
reserve the next n nodes for allocating to leave in order of insertion

this allows direct indexing of leaves with implicit in-order indices
uses any exising free nodes, so if those aren't ordered then this will not be true
"""
#==
function reserve_leaves!( nv::NodeVector{T}, n:Int ) where T
  nv.nextfreeleaf = nv.last
  oldsz = length(nv.nodes)
  resize!(nv, n+oldsz+1)
  nv.nextfree = n
  nv.nodes[nv.nextfreeleaf+n-1].next = EmptyIndexedNode

end
==#

Base.size(A::NodeVector{T}) where T = size(A.nodes)
Base.getindex(A::NodeVector{T}, i::Number) where {T} = A.nodes[i]
Base.setindex!(A::NodeVector{T}, n::IndexedNode{T}, i::Number) where {T} = A.nodes[i] = n
Base.firstindex(A::NodeVector{T}) where T = A.root

itree(ibt::IndexedBinaryTree{T}, node::Int) where T = IndexedBinaryTree{T}(ibt.nodes, node)
itree(a::T) where T = IndexedBinaryTree{T}(NodeVector{T}([IndexedNode(a)], 1, EmptyIndexedNode, 1), 1)
itree(dc::V) where {T, V<:AbstractVector{T}} = reduce(insert, Iterators.drop(dc,1); init = itree(first(dc)))

key(tree::IndexedBinaryTree{T}) where T = key(tree.nodes[tree.node].data)

left(tree::IndexedBinaryTree{T}) where T = isempty_node(tree.nodes[tree.node].left) ? nothing : itree(tree, tree.nodes[tree.node].left)

right(tree::IndexedBinaryTree{T}) where T = isempty_node(tree.nodes[tree.node].right) ? nothing : itree(tree, tree.nodes[tree.node].right)

ileft(tree::IndexedBinaryTree{T}) where T = tree.nodes[tree.node].left
iright(tree::IndexedBinaryTree{T}) where T = tree.nodes[tree.node].right

parent(tree::IndexedBinaryTree{T}) where T = isempty_node(tree.nodes[tree.node].parent) ? nothing : itree(tree, parent(tree.nodes[tree.node]))

isempty_node(tree::IndexedBinaryTree{T}) where T = tree.node == EmptyIndexedNode
isempty_node(i::Int) = i == EmptyIndexedNode

data(tree::IndexedBinaryTree{T}) where T = tree.nodes[tree.node].data
keyval(tree::IndexedBinaryTree{T}) where T = data(tree)
data(n::IndexedNode{T}) where T = n.data
keyval(n::IndexedNode{T}) where T = n.data

function search(stree::IndexedBinaryTree{T}, k::T) where T
  hit = nothing
  while true
    if k == key(stree)
      hit = (stree, :found)
      break
    elseif k < key(stree)
      if isnothing(left(stree))
        hit = (stree, :left)
        break
      else
        stree = left(stree)
      end
    else
      if isnothing(right(stree))
        hit = (stree, :right)
        break
      else
        stree = right(stree)
      end
    end
  end
  return hit
end

#==
function insert(::Type{NotFullTree}, t::IndexedBinaryTree{T}, d::T) where T
  nv = t.nodes
  if t.nodes.nextfree == EmptyIndexedNode
    nv = resize(nv, length(nv)+1)
  end
  inode = nv.nextfree
  nv_nextfree = nv[inode].next

  hit = search(t, d)
  hit = isnothing(hit) ? (t, :found) : hit

  iparent = first(hit).node
  p = nv[iparent]
  if hit[2] == :left
    p = IndexedNode{T}(inode, p.right, parent(p), p.data)
  else
    p = IndexedNode{T}(p.left, inode, parent(p), p.data)
  end

  cpnodes = copy(nv.nodes)
  cpnodes[inode] = IndexedNode{T}(EmptyIndexedNode, EmptyIndexedNode, iparent, d)
  cpnodes[iparent] =  p
  IndexedBinaryTree{T}(NodeVector{T}(cpnodes, nv.root, nv_nextfree, nv.last), nv.root)
end

import Base.insert!

insert!( t::N, d::T ) where {T, N<:BinaryTree{T}} = insert!( fullness(N), t, d )

function insert!(::Type{NotFullTree}, tree::IndexedBinaryTree{T}, d::T) where T
  nv = tree.nodes
  if tree.nodes.nextfree == EmptyIndexedNode
    nv = resize!(nv, length(nv)+1)
  end
  inode = nv.nextfree
  nv_nextfree = nv[inode].next

  hit = search(tree, d);
  hit = isnothing(hit) ? (tree, :found) : hit

  iparent = first(hit).node
  p = nv[iparent]
  if hit[2] == :left
    p = IndexedNode{T}(inode, p.right, parent(p), p.data)
  else
    p = IndexedNode{T}(p.left, inode, parent(p), p.data)
  end

  cpnodes = nv.nodes
  cpnodes[inode] = IndexedNode{T}(EmptyIndexedNode, EmptyIndexedNode, iparent, d)
  cpnodes[iparent] =  p
  IndexedBinaryTree{T}(NodeVector{T}(cpnodes, nv.root, nv_nextfree, nv.last), nv.root)
end
==#

function rebuild_node!(dir::Union{Val{:both}, Val{:left}, Val{:right}},
             lc::Branch{T}, rc::Branch{T},
             n::IndexedBinaryTree{T}, d::T ) where T

  nv = n.nodes
  inode = n.node
  if isnothing(lc)
    ilc = EmptyIndexedNode
  elseif dir != Val{:right}() && isleaf(lc)
    # copy leaf into parents vector
    if nv.nextfree == EmptyIndexedNode
      nv = resize!(nv, length(nv)+2)  # adding two for case of 2 children
    end

    #==
    #!me possible orphaning? then we need something like this, but then that should looks like
    #!me else case after this
    #!me can happen for full trees (:both) or rotations we move a node from this tree to this tree?
    if !isroot(lc)
      nv[index(lc)] = setnext(node(lc), nv.nextfree)
      nv.nextfree = index(lc)
    end
    ==#

    ilc = nv.nextfree   # if not new index it will orphan current lc. else nv[index(lc)] = setparent(node(lc), nv.nextfree);  nv.nextfree = index(lc);  # maybe put in free() method
    nv.nextfree = nv[ilc].next
    nv[ilc] = setparent(node(lc), inode)
  else
    # child and parents vector should be the same
    ilc = index(lc)
    nv[ilc] = setparent(node(lc), inode)
  end

  if isnothing(rc)
    irc = EmptyIndexedNode
  elseif dir != Val{:left}() && isleaf(rc)
    # copy leaf into parents vector
    if nv.nextfree == EmptyIndexedNode
      nv = resize!(nv, length(nv)+1)
    end
    irc = nv.nextfree
    nv.nextfree = nv[irc].next
    nv[irc] = setparent(node(rc), inode)
  else
    # child and parents vector should be the same
    irc = index(rc)
    nv[irc] = setparent(node(rc), inode)
  end

  imabastard = EmptyIndexedNode
  tn = IndexedNode{T}(ilc, irc, imabastard, d)
  nv[inode] = tn
  IndexedBinaryTree{T}(nv, inode)
end

rebuild_node(dir::Union{Val{:both}, Val{:left}, Val{:right}},
             lc::Branch{T}, rc::Branch{T},
             n::IndexedBinaryTree{T}, d::T ) where T =
  rebuild_node!(dir, lc, rc, n, d)


text:
# GurobiMathProgInterface
# Standardized MILP interface

export GurobiSolver

type GurobiMathProgModel <: AbstractLinearQuadraticModel
    inner::Model
    last_op_type::Symbol  # To support arbitrary order of addVar/addCon
                          # Two possibilities :Var :Con
    changed_constr_bounds::Bool # have we updated the bounds below?
    lb::Vector{Float64} # persistent bounds on constraints to maintain
    ub::Vector{Float64} # abstraction for lb ≤ Ax ≤ ub
    lazycb
    cutcb
    heuristiccb
    infocb
    options
end
function GurobiMathProgModel(env=nothing;options...)
   finalize_env = (env == nothing)
   if env == nothing
       env = Env()
   end
   m = GurobiMathProgModel(Model(env,""; finalize_env=finalize_env), :Con, false, Float64[], Float64[], nothing, nothing, nothing, nothing, options)
   setparams!(m)
   return m
end
function setparams!(m::GurobiMathProgModel)
    # Helper to set the parameters on the model's copy of env rather than
    # modifying the global env (ref: http://www.gurobi.com/support/faqs#P)

    # Set `InfUnbdInfo` to 1 by default so infeasibility rays available
    setparam!(m.inner, "InfUnbdInfo", 1)
    for (name,value) in m.options
        setparam!(m.inner, string(name), value)
    end
end

immutable GurobiSolver <: AbstractMathProgSolver
    env
    options
end
GurobiSolver(env=nothing; kwargs...) = GurobiSolver(env, kwargs)
LinearQuadraticModel(s::GurobiSolver) = GurobiMathProgModel(s.env; s.options...)
ConicModel(s::GurobiSolver) = LPQPtoConicBridge(LinearQuadraticModel(s))

supportedcones(::GurobiSolver) = [:Free,:Zero,:NonNeg,:NonPos,:SOC,:SOCRotated]

loadproblem!(m::GurobiMathProgModel, filename::AbstractString) = read_model(m.inner, filename)

function loadproblem!(m::GurobiMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
  # throw away old model but keep env and finalize_env
  env = m.inner.env
  finalize_env = m.inner.finalize_env
  m.inner.finalize_env = false
  free_model(m.inner)
  m.inner = Model(env, "", finalize_env=finalize_env)
  # Re-set options on new model's copy of env
  setparams!(m)

  add_cvars!(m.inner, float(obj), float(collb), float(colub))
  update_model!(m.inner)

  neginf = typemin(eltype(rowlb))
  posinf = typemax(eltype(rowub))

  # check if we have any range constraints
  # to properly support these, we will need to keep track of the
  # slack variables automatically added by gurobi.
  rangeconstrs = any((rowlb .!= rowub) & (rowlb .> neginf) & (rowub .< posinf))
  if rangeconstrs
      warn("Julia Gurobi interface doesn't properly support range (two-sided) constraints. See Gurobi.jl issue #14")
      add_rangeconstrs!(m.inner, float(A), float(rowlb), float(rowub))
  else
      b = Array(Float64,length(rowlb))
      senses = Array(Cchar,length(rowlb))
      for i in 1:length(rowlb)
          if rowlb[i] == rowub[i]
              senses[i] = '='
              b[i] = rowlb[i]
          elseif rowlb[i] > neginf
              senses[i] = '>'
              b[i] = rowlb[i]
          else
              @assert rowub[i] < posinf
              senses[i] = '<'
              b[i] = rowub[i]
          end
      end
      add_constrs!(m.inner, float(A), senses, b)
  end

  m.lb, m.ub = rowlb, rowub

  update_model!(m.inner)
  setsense!(m,sense)
end

writeproblem(m::GurobiMathProgModel, filename::AbstractString) = write_model(m.inner, filename)

getvarLB(m::GurobiMathProgModel)     = get_dblattrarray( m.inner, "LB", 1, num_vars(m.inner))
setvarLB!(m::GurobiMathProgModel, l) = set_dblattrarray!(m.inner, "LB", 1, num_vars(m.inner), l)

getvarUB(m::GurobiMathProgModel)     = get_dblattrarray( m.inner, "UB", 1, num_vars(m.inner))
setvarUB!(m::GurobiMathProgModel, u) = set_dblattrarray!(m.inner, "UB", 1, num_vars(m.inner), u)

function getconstrLB(m::GurobiMathProgModel)
    sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
    ret   = get_dblattrarray(m.inner, "RHS", 1, num_constrs(m.inner))
    for i = 1:num_constrs(m.inner)
        if sense[i] == '>' || sense[i] == '='
            # Do nothing
        else
            # LEQ constraint, so LB is -Inf
            ret[i] = -Inf
        end
     end
     return ret
end

function getconstrUB(m::GurobiMathProgModel)
    sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
    ret   = get_dblattrarray(m.inner, "RHS", 1, num_constrs(m.inner))
    for i = 1:num_constrs(m.inner)
        if sense[i] == '<' || sense[i] == '='
            # Do nothing
        else
            # GEQ constraint, so UB is +Inf
            ret[i] = +Inf
        end
    end
    return ret
end

# setconstrLB!(m::GurobiMathProgModel, lb) = (m.changed_constr_bounds = true; m.last_op_type = :Con; m.lb = copy(lb))
# setconstrUB!(m::GurobiMathProgModel, ub) = (m.changed_constr_bounds = true; m.last_op_type = :Con; m.ub = copy(ub))
setconstrLB!(m::GurobiMathProgModel, lb) = (m.changed_constr_bounds = true; m.lb = copy(lb))
setconstrUB!(m::GurobiMathProgModel, ub) = (m.changed_constr_bounds = true; m.ub = copy(ub))

getobj(m::GurobiMathProgModel)     = get_dblattrarray( m.inner, "Obj", 1, num_vars(m.inner)   )
setobj!(m::GurobiMathProgModel, c) = set_dblattrarray!(m.inner, "Obj", 1, num_vars(m.inner), c)

function addvar!(m::GurobiMathProgModel, constridx, constrcoef, l, u, objcoef)
    if m.last_op_type == :Con
        updatemodel!(m)
        m.last_op_type = :Var
    end
    add_var!(m.inner, length(constridx), constridx, float(constrcoef), float(objcoef), float(l), float(u), GRB_CONTINUOUS)
end
function addvar!(m::GurobiMathProgModel, l, u, objcoef)
    if m.last_op_type == :Con
        updatemodel!(m)
        m.last_op_type = :Var
    end
    add_var!(m.inner, 0, Integer[], Float64[], float(objcoef), float(l), float(u), GRB_CONTINUOUS)
end
function addconstr!(m::GurobiMathProgModel, varidx, coef, lb, ub)
    if m.last_op_type == :Var
        updatemodel!(m)
        m.last_op_type = :Con
    end
    if lb == -Inf
        # <= constraint
        add_constr!(m.inner, varidx, coef, '<', ub)
    elseif ub == +Inf
        # >= constraint
        add_constr!(m.inner, varidx, coef, '>', lb)
    elseif lb == ub
        # == constraint
        add_constr!(m.inner, varidx, coef, '=', lb)
    else
        # Range constraint
        error("Adding range constraints not supported yet.")
    end
end

function updatemodel!(m::GurobiMathProgModel)
    update_model!(m.inner)
    if m.changed_constr_bounds
    # update lower/upper bounds on linear constraints (if they're consistent...)
        sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
        rhs   = get_dblattrarray( m.inner, "RHS",   1, num_constrs(m.inner))
        lb, ub = m.lb, m.ub
        @assert (n_rows = num_constrs(m.inner)) == length(lb) == length(ub) ==
                                                   length(sense) == length(rhs)
        for i = 1:n_rows
            if -Inf < lb[i] == ub[i] < Inf
                sense[i] = '='
                rhs[i]   = lb[i]
            elseif (-Inf < lb[i] && ub[i] < Inf)
                error("Gurobi.jl does not currently support range constraints")
            elseif (lb[i] == -Inf && ub[i] < Inf)
                sense[i] = '<'
                rhs[i]   = ub[i]
            elseif (lb[i] > -Inf && ub[i] == Inf)
                sense[i] = '>'
                rhs[i]   = lb[i]
            else
                error("Internal error.")
            end
        end
        set_charattrarray!(m.inner, "Sense", 1, num_constrs(m.inner), sense)
        set_dblattrarray!(m.inner, "RHS", 1, num_constrs(m.inner), rhs)
        m.changed_constr_bounds = false
        update_model!(m.inner)
    end
end

getconstrmatrix(m::GurobiMathProgModel) = get_constrmatrix(m.inner)

function setsense!(m::GurobiMathProgModel, sense)
  if sense == :Min
    set_sense!(m.inner, :minimize)
  elseif sense == :Max
    set_sense!(m.inner, :maximize)
  else
    error("Unrecognized objective sense $sense")
  end
end
function getsense(m::GurobiMathProgModel)
  v = get_intattr(m.inner, "ModelSense")
  if v == -1
    return :Max
  else
    return :Min
  end
end

numvar(m::GurobiMathProgModel)    = (updatemodel!(m); num_vars(m.inner))
numconstr(m::GurobiMathProgModel) = num_constrs(m.inner) + num_qconstrs(m.inner)
numlinconstr(m::GurobiMathProgModel) = num_constrs(m.inner)
numquadconstr(m::GurobiMathProgModel) = num_qconstrs(m.inner)

function optimize!(m::GurobiMathProgModel)
    # set callbacks if present
    if m.lazycb != nothing || m.cutcb != nothing || m.heuristiccb != nothing || m.infocb != nothing
        updatemodel!(m)
        if !is_mip(m.inner)
            Base.warn_once("Gurobi ignores branch-and-bound callbacks when no discrete elements are present in the model.")
        end
        setmathprogcallback!(m)
    end
    if m.lazycb != nothing
      setparam!(m.inner, "LazyConstraints", 1)
    end
    updatemodel!(m)
    optimize(m.inner)
end

function status(m::GurobiMathProgModel)
  s = get_status(m.inner)
  if s == :optimal
    return :Optimal
  elseif s == :infeasible
    return :Infeasible
  elseif s == :unbounded
    return :Unbounded
  elseif s == :inf_or_unbd
    Base.warn_once("Gurobi reported infeasible or unbounded. Set InfUnbdInfo=1 for more specific status.")
    return :InfeasibleOrUnbounded
  elseif s == :iteration_limit || s == :node_limit || s == :time_limit || s == :solution_limit
    return :UserLimit
  elseif s == :numeric
    return :Numeric
  elseif s == :suboptimal
    return :Suboptimal # not very useful status
  elseif s == :interrupted # ended by user?
    return :UserLimit
  else
    error("Unrecognized solution status: $s")
  end
end

getobjval(m::GurobiMathProgModel)   = get_objval(m.inner)
getobjbound(m::GurobiMathProgModel) = get_objbound(m.inner)
getsolution(m::GurobiMathProgModel) = get_solution(m.inner)

function getconstrsolution(m::GurobiMathProgModel)
    sense = get_charattrarray(m.inner, "Sense", 1, num_constrs(m.inner))
    rhs   = get_dblattrarray( m.inner, "RHS",   1, num_constrs(m.inner))
    slack = get_dblattrarray( m.inner, "Slack", 1, num_constrs(m.inner))
    ret   = zeros(num_constrs(m.inner))
    for i = 1:num_constrs(m.inner)
        if sense[i] == '='
            # Must be equal to RHS if feasible
            ret[i] = rhs[i]
        else
            # Slack variable is non-negative for <= constraints, and non-positive for >= constraints
            ret[i] = rhs[i] - slack[i]
        end
    end
    return ret
end

function getreducedcosts(m::GurobiMathProgModel)
    if is_qcp(m.inner) && get_int_param(m.inner, "QCPDual") == 0
        return fill(NaN, num_vars(m.inner))
    else
        return get_dblattrarray(m.inner, "RC", 1, num_vars(m.inner))
    end
end

function getconstrduals(m::GurobiMathProgModel)
    if is_qcp(m.inner) && get_int_param(m.inner, "QCPDual") == 0
        return fill(NaN, num_constrs(m.inner))
    else
        return get_dblattrarray(m.inner, "Pi", 1, num_constrs(m.inner))
    end
end

function getquadconstrduals(m::GurobiMathProgModel)
    if is_qcp(m.inner) && get_int_param(m.inner, "QCPDual") == 0
        return fill(NaN, num_qconstrs(m.inner))
    else
        return get_dblattrarray(m.inner, "QCPi", 1, num_qconstrs(m.inner))
    end
end

getinfeasibilityray(m::GurobiMathProgModel) = -get_dblattrarray(m.inner, "FarkasDual", 1, num_constrs(m.inner)) # note sign is flipped
getunboundedray(m::GurobiMathProgModel) = get_dblattrarray(m.inner, "UnbdRay", 1, num_vars(m.inner))

getbasis(m::GurobiMathProgModel) = get_basis(m.inner)

getrawsolver(m::GurobiMathProgModel) = m.inner

const var_type_map = Dict(
  'C' => :Cont,
  'B' => :Bin,
  'I' => :Int,
  'S' => :SemiCont,
  'N' => :SemiInt
)

const rev_var_type_map = Dict(
  :Cont => 'C',
  :Bin => 'B',
  :Int => 'I',
  :SemiCont => 'S',
  :SemiInt => 'N'
)

function setvartype!(m::GurobiMathProgModel, vartype::Vector{Symbol})
    # do this to make sure we deal with new columns
    updatemodel!(m)
    nvartype = map(x->rev_var_type_map[x], vartype)
    set_charattrarray!(m.inner, "VType", 1, length(nvartype), nvartype)
    updatemodel!(m) # otherwise getvartype! will return old values
end

function getvartype(m::GurobiMathProgModel)
    ret = get_charattrarray(m.inner, "VType", 1, num_vars(m.inner))
    map(x->var_type_map[x], ret)
end

function setwarmstart!(m::GurobiMathProgModel, v)
    for j = 1:length(v)
        if isnan(v[j])
            v[j] = 1e101  # GRB_UNDEFINED
        end
    end
    set_dblattrarray!(m.inner, "Start", 1, num_vars(m.inner), v)
end

addsos1!(m::GurobiMathProgModel, idx, weight) = add_sos!(m.inner, :SOS1, idx, weight)
addsos2!(m::GurobiMathProgModel, idx, weight) = add_sos!(m.inner, :SOS2, idx, weight)

# Callbacks

setlazycallback!(m::GurobiMathProgModel,f) = (m.lazycb = f)
setcutcallback!(m::GurobiMathProgModel,f) = (m.cutcb = f)
setheuristiccallback!(m::GurobiMathProgModel,f) = (m.heuristiccb = f)
setinfocallback!(m::GurobiMathProgModel,f) = (m.infocb = f)

type GurobiCallbackData <: MathProgCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    sol::Vector{Float64}  # Used for heuristic callbacks
#    model::GurobiMathProgModel # not needed?
end

function cbgetmipsolution(d::GurobiCallbackData)
    @assert d.state == :MIPSol
    return cbget_mipsol_sol(d.cbdata, d.where)
end

function cbgetmipsolution(d::GurobiCallbackData,output)
    @assert d.state == :MIPSol
    return cbget_mipsol_sol(d.cbdata, d.where, output)
end

function cbgetlpsolution(d::GurobiCallbackData)
    @assert d.state == :MIPNode
    return cbget_mipnode_rel(d.cbdata, d.where)
end

function cbgetlpsolution(d::GurobiCallbackData, output)
    @assert d.state == :MIPNode
    return cbget_mipnode_rel(d.cbdata, d.where, output)
end

# TODO: macro for these getters?
function cbgetobj(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_objbst(d.cbdata, d.where)
    elseif d.state == :Intermediate
        return cbget_mip_objbst(d.cbdata, d.where)
    elseif d.state == :MIPSol
        error("Gurobi does not implement cbgetobj when state == MIPSol")
        # https://groups.google.com/forum/#!topic/gurobi/Az_X6Ag-y6k
        # https://github.com/JuliaOpt/MathProgBase.jl/issues/23
        # A complex workaround is possible but hasn't been implemented.
        # return cbget_mipsol_objbst(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

function cbgetbestbound(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_objbnd(d.cbdata, d.where)
    elseif d.state == :MIPSol
        return cbget_mipsol_objbnd(d.cbdata, d.where)
    elseif d.state == :Intermediate
        return cbget_mip_objbnd(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

function cbgetexplorednodes(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_nodcnt(d.cbdata, d.where)
    elseif d.state == :MIPSol
        return cbget_mipsol_nodcnt(d.cbdata, d.where)
    elseif d.state == :Intermediate
        return cbget_mip_nodcnt(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

# returns :MIPNode :MIPSol :Intermediate
cbgetstate(d::GurobiCallbackData) = d.state

function cbaddcut!(d::GurobiCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode
    cbcut(d.cbdata, convert(Vector{Cint}, varidx), float(varcoef), sense, float(rhs))
end

function cbaddlazy!(d::GurobiCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode || d.state == :MIPSol
    cblazy(d.cbdata, convert(Vector{Cint}, varidx), float(varcoef), sense, float(rhs))
end

function cbaddsolution!(d::GurobiCallbackData)
    # Gurobi doesn't support adding solutions on MIPSol.
    # TODO: support this anyway
    @assert d.state == :MIPNode
    cbsolution(d.cbdata, d.sol)
    # "Wipe" solution back to GRB_UNDEFINIED
    for i in 1:length(d.sol)
        d.sol[i] = 1e101  # GRB_UNDEFINED
    end
end

function cbsetsolutionvalue!(d::GurobiCallbackData,varidx,value)
    d.sol[varidx] = value
end

# breaking abstraction, define our low-level callback to eliminatate
# a level of indirection

function mastercallback(ptr_model::Ptr{Void}, cbdata::Ptr{Void}, where::Cint, userdata::Ptr{Void})

    model = unsafe_pointer_to_objref(userdata)::GurobiMathProgModel
    grbrawcb = CallbackData(cbdata,model.inner)
    if where == CB_MIPSOL
        state = :MIPSol
    elseif where == CB_MIPNODE
        state = :MIPNode
        # skip callback if node is reported to be cut off or infeasible --
        # nothing to do.
        # TODO: users may want this information
        status = cbget_mipnode_status(grbrawcb, where)
        if status != 2
            return convert(Cint,0)
        end
    elseif where == CB_MIP
        state = :Intermediate
    else
        # State with no relevant callbacks
        return convert(Cint,0)
    end

    grbcb = GurobiCallbackData(grbrawcb, state, where, [0.0])
    if model.infocb != nothing
        ret = model.infocb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end
    if model.cutcb != nothing && state == :MIPNode
        ret = model.cutcb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end
    if model.heuristiccb != nothing && state == :MIPNode
        grbcb.sol = fill(1e101, numvar(model))  # GRB_UNDEFINED
        ret = model.heuristiccb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end
    if model.lazycb != nothing && (state == :MIPSol || state == :MIPNode)
        ret = model.lazycb(grbcb)
        if ret == :Exit
            terminate(model.inner)
        end
    end

    return convert(Cint,0)
end

# User callback function should be of the form:
# callback(cbdata::MathProgCallbackData)
# return :Exit to indicate an error

function setmathprogcallback!(model::GurobiMathProgModel)

    is_windows() && (Sys.WORD_SIZE == 64 || error("Callbacks not currently supported on Win32. Use 64-bit Julia with 64-bit Gurobi."))
    grbcallback = cfunction(mastercallback, Cint, (Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}))
    ret = @grb_ccall(setcallbackfunc, Cint, (Ptr{Void}, Ptr{Void}, Any), model.inner.ptr_model, grbcallback, model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

# QCQP

function setquadobj!(m::GurobiMathProgModel, rowidx, colidx, quadval)
    delq!(m.inner)
    update_model!(m.inner)
    scaledvals = similar(quadval)
    for i in 1:length(rowidx)
      if rowidx[i] == colidx[i]
        # rescale from matrix format to "terms" format
        scaledvals[i] = quadval[i] / 2
      else
        scaledvals[i] = quadval[i]
      end
    end
    add_qpterms!(m.inner, rowidx, colidx, scaledvals)
    update_model!(m.inner)
end

function addquadconstr!(m::GurobiMathProgModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    if m.last_op_type == :Var
        updatemodel!(m)
        m.last_op_type = :Con
    end
    add_qconstr!(m.inner, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
end

getsolvetime(m::GurobiMathProgModel) = get_runtime(m.inner)
getnodecount(m::GurobiMathProgModel) = get_node_count(m.inner)


text:
@doc raw"""
    TranslationGroup{T<:Tuple,𝔽} <: GroupManifold{Euclidean{T,𝔽},AdditionOperation}

Translation group $\mathrm{T}(n)$ represented by translation arrays.

# Constructor
    TranslationGroup(n₁,...,nᵢ; field = 𝔽)

Generate the translation group on
$𝔽^{n₁,…,nᵢ}$ = `Euclidean(n₁,...,nᵢ; field = 𝔽)`, which is isomorphic to the group itself.
"""
const TranslationGroup{T<:Tuple,𝔽} = GroupManifold{𝔽,Euclidean{T,𝔽},AdditionOperation}

function TranslationGroup(n::Int...; field::AbstractNumbers=ℝ)
    return TranslationGroup{Tuple{n...},field}(
        Euclidean(n...; field=field),
        AdditionOperation(),
    )
end

@inline function active_traits(f, M::TranslationGroup, args...)
    if is_metric_function(f)
        #pass to Euclidean by default - but keep Group Decorator for the retraction
        return merge_traits(IsGroupManifold(M.op), IsExplicitDecorator())
    else
        return merge_traits(
            IsGroupManifold(M.op),
            HasBiinvariantMetric(),
            IsDefaultMetric(EuclideanMetric()),
            active_traits(f, M.manifold, args...),
            IsExplicitDecorator(), #pass to Euclidean by default/last fallback
        )
    end
end

identity_element!(::TranslationGroup, p) = fill!(p, 0)

function Base.show(io::IO, ::TranslationGroup{N,𝔽}) where {N,𝔽}
    return print(io, "TranslationGroup($(join(N.parameters, ", ")); field = $(𝔽))")
end


text:
module cylinder_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule

function cylinder_bend()
    # Cylinder  bent by edge by enforced displacement, axially symmetric model
    println("Cylinder  bent by edge by enforced displacement, axially symmetric model")

    # Parameters:
    E1=1.0;
    E2=1.0;
    E3=3.0;
    nu12=0.29;
    nu13=0.29;
    nu23=0.19;
    G12=0.3;
    G13=0.3;
    G23=0.3;
    p= 0.15;
    rin=1.;
    rex =1.2;
    Length = 1*rex
    ua = -0.05*Length
    tolerance=rin/1000.

    ##
    # Note that the FinEtools objects needs to be created with the proper
    # model-dimension reduction at hand.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DAxisymm

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the thickness in the radial direction, and the second
    # coordinate is the thickness in the axial direction.

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] .+= rin
    bdryfes = meshboundary(fes);

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # the symmetry plane
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    # The other end
    l1 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)
    for ixxxx = 1:length(l1)
        r = fens.xyz[l1[ixxxx], 1]
        setebc!(u,[l1[ixxxx]],true, 2, (r-(rex+rin)/2)/((rex+rin)/2)*ua)
    end

    applyebc!(u)
    numberdofs!(u)
    println("Number of degrees of freedom = $(u.nfreedofs)")

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)

    K =stiffness(femm, geom, u)
    F = nzebcloadsstiffness(femm, geom, u)

    U=  K\(F)
    scattersysvec!(u,U[:])

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
    println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")

    File =  "cylinder_bend_sigmaz.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)

end # cylinder_bend

function cylinder_pull()
    # Cylinder  pulled by enforced displacement, axially symmetric model
    println("Cylinder  pulled by enforced displacement, axially symmetric model")

    # Parameters:
    E1=1.0;
    E2=1.0;
    E3=3.0;
    nu12=0.29;
    nu13=0.29;
    nu23=0.19;
    G12=0.3;
    G13=0.3;
    G23=0.3;
    p= 0.15;
    rin=1.;
    rex =1.2;
    Length = 1*rex
    ua = -0.05*Length
    tolerance=rin/1000.

    ##
    # Note that the FinEtools objects needs to be created with the proper
    # model-dimension reduction at hand.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DAxisymm

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the thickness in the radial direction, and the second
    # coordinate is the thickness in the axial direction.

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] .+= rin
    bdryfes = meshboundary(fes);

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # the symmetry plane
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    # The other end
    l1 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)
    setebc!(u,l1,true, 2, ua)

    applyebc!(u)
    numberdofs!(u)
    println("Number of degrees of freedom = $(u.nfreedofs)")

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)

    K =stiffness(femm, geom, u)
    F = nzebcloadsstiffness(femm, geom, u)
    U=  K\(F)
    scattersysvec!(u,U[:])

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
    println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")

    File =  "orthoballoon_sigmaz.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)

end # cylinder_pull

function cylinder_pull_algo()
    # Cylinder  pulled by enforced displacement, axially symmetric model

    # Parameters:
    E1=1.0;
    E2=1.0;
    E3=3.0;
    nu12=0.29;
    nu13=0.29;
    nu23=0.19;
    G12=0.3;
    G13=0.3;
    G23=0.3;
    p= 0.15;
    rin=1.;
    rex =1.2;
    Length = 1*rex
    ua = -0.05*Length
    tolerance=rin/1000.

    ##
    # Note that the FinEtools objects needs to be created with the proper
    # model-dimension reduction at hand.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DAxisymm

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the thickness in the radial direction, and the second
    # coordinate is the thickness in the axial direction.

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] .+= rin
    bdryfes = meshboundary(fes);

    # the symmetry plane
    la1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    # The other end
    la2 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)

    e1 = FDataDict("node_list"=>la1, "component"=>2, "displacement"=>x -> 0.0)
    e2 = FDataDict("node_list"=>la2, "component"=>2, "displacement"=>x -> ua)

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)

    # Make region
    region = FDataDict("femm"=>femm)

    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region],
    "essential_bcs"=>[e1, e2])

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)
    geom = modeldata["geom"]
    u = modeldata["u"]

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
    println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")

    File =  "orthoballoon_sigmaz.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
end # cylinder_pull_algo

function allrun()
    println("#####################################################")
    println("# cylinder_bend ")
    cylinder_bend()
    println("#####################################################")
    println("# cylinder_pull ")
    cylinder_pull()
    println("#####################################################")
    println("# cylinder_pull_algo ")
    cylinder_pull_algo()
    return true
end # function allrun

end # module cylinder_examples


text:
using Luxor

function drawbars(w, h, sections, dk, lt)
    Drawing(w,h)
    background("white")
    width = 1
    height = h/sections
    for y in 0:height:h-1
        setline(width)
        for x in 0:w/width
            sethue(x % 2 == 0 ? dk: lt)
            line(Point(x*width,y), Point(x*width,y+height), :stroke)
        end
        width += 1
    end
end

drawbars(1920, 1080, 4, "black", "white")
finish()
preview()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# canonical data version: 1.7.0

using Test

include("acronym.jl")

@testset "basic" begin
    @test acronym("Portable Network Graphics") == "PNG"
end

@testset "lowercase words" begin
    @test acronym("Ruby on Rails") == "ROR"
end

@testset "punctuation" begin
    @test acronym("First In, First Out") == "FIFO"
end

@testset "all caps word" begin
    @test acronym("GNU Image Manipulation Program") == "GIMP"
end

@testset "punctuation without whitespace" begin
    @test acronym("Complementary metal-oxide semiconductor") == "CMOS"
end

@testset "very long abbreviation" begin
    @test acronym("Rolling On The Floor Laughing So Hard That My Dogs Came Over And Licked Me") == "ROTFLSHTMDCOALM"
end

@testset "consecutive delimiters" begin
    @test acronym("Something - I made up from thin air") == "SIMUFTA"
end

@testset "apostrophes" begin
    @test acronym("Halley's Comet") == "HC"
end

@testset "underscore emphasis" begin
    @test acronym("The Road _Not_ Taken") == "TRNT"
end


text:
const SpinOneSite = Union{TagType"S=1", TagType"SpinOne"}

function siteinds(::SpinOneSite,
                  N::Int; kwargs...)
  conserve_qns = get(kwargs,:conserve_qns,false)
  conserve_sz = get(kwargs,:conserve_sz,conserve_qns)
  if conserve_sz
    up = QN("Sz",+2) => 1
    z0 = QN("Sz", 0) => 1
    dn = QN("Sz",-2) => 1
    return [Index(up,z0,dn;tags="Site,S=1,n=$n") for n=1:N]
  end
  return [Index(3,"Site,S=1,n=$n") for n=1:N]
end

function state(::SpinOneSite,
               st::AbstractString)
  if st == "Up" || st == "↑"
    return 1
  elseif st == "Z0" || st == "0"
    return 2
  elseif st == "Dn" || st == "↓"
    return 3
  end
  throw(ArgumentError("State string \"$st\" not recognized for SpinOne site"))
  return 0
end

function op(::SpinOneSite,
            s::Index,
            opname::AbstractString)::ITensor
  Up = s(1)
  UpP = s'(1)
  Z0 = s(2)
  Z0P = s'(2)
  Dn = s(3)
  DnP = s'(3)

  Op = emptyITensor(s',dag(s))

  if opname == "S⁺" || opname == "Splus" || opname == "S+"
    Op[Z0P, Dn] = √2
    Op[UpP, Z0] = √2
  elseif opname == "S⁻" || opname == "Sminus" || opname == "S-"
    Op[Z0P, Up] = √2
    Op[DnP, Z0] = √2
  elseif opname == "Sˣ" || opname == "Sx"
    Op[Z0P, Up] = 1.0/√2
    Op[UpP, Z0] = 1.0/√2
    Op[DnP, Z0] = 1.0/√2
    Op[Z0P, Dn] = 1.0/√2
  elseif opname == "iSʸ" || opname == "iSy"
    Op[Z0P, Up] = -1.0/√2
    Op[UpP, Z0] = +1.0/√2
    Op[DnP, Z0] = -1.0/√2
    Op[Z0P, Dn] = +1.0/√2
  elseif opname == "Sʸ" || opname == "Sy"
    Op = complex(Op)
    Op[Z0P, Up] = +1.0/√2im
    Op[UpP, Z0] = -1.0/√2im
    Op[DnP, Z0] = +1.0/√2im
    Op[Z0P, Dn] = -1.0/√2im
  elseif opname == "Sᶻ" || opname == "Sz"
    Op[UpP, Up] = 1.0
    Op[DnP, Dn] = -1.0
  elseif opname == "Sᶻ²" || opname == "Sz2"
    Op[UpP, Up] = 1.0
    Op[DnP, Dn] = 1.0
  elseif opname == "Sˣ²" || opname == "Sx2"
    Op[UpP, Up] = 0.5
    Op[DnP, Up] = 0.5
    Op[Z0P, Z0] = 1.0
    Op[UpP, Dn] = 0.5
    Op[DnP, Dn] = 0.5
  elseif opname == "Sʸ²" || opname == "Sy2"
    Op[UpP, Up] = 0.5
    Op[DnP, Up] = -0.5
    Op[Z0P, Z0] = 1.0
    Op[UpP, Dn] = -0.5
    Op[DnP, Dn] = 0.5
  elseif opname == "projUp"
    Op[UpP, Up] = 1.
  elseif opname == "projZ0"
    Op[Z0P, Z0] = 1.
  elseif opname == "projDn"
    Op[DnP, Dn] = 1.
  elseif opname == "XUp"
    xup = emptyITensor(ComplexF64,s)
    xup[Up] = 0.5
    xup[Z0] = im*√2
    xup[Dn] = 0.5
    return xup
  elseif opname == "XZ0"
    xZ0 = emptyITensor(ComplexF64,s)
    xZ0[Up] = im*√2
    xZ0[Dn] = -im*√2
    return xZ0
  elseif opname == "XDn"
    xdn = emptyITensor(ComplexF64,s)
    xdn[Up] = 0.5
    xdn[Z0] = -im*√2
    xdn[Dn] = 0.5
    return xdn
  else
    throw(ArgumentError("Operator name '$opname' not recognized for SpinOneSite"))
  end
  return Op
end


text:
################################################################################
# Create Rust crate
mkpath("AddIntegersRust/add-integers")
write("AddIntegersRust/add-integers/Cargo.toml", """
                                                   [package]
                                                   name = "add-integers"
                                                   version = "1.0.0"
                                                   edition = "2021"

                                                   [lib]
                                                   # The library name cannot be camel case
                                                   name = "add_integers_rust"
                                                   crate-type = ["cdylib"]

                                                   [dependencies]
                                                   half = "1.8.2"
                                                   """)
mkpath("AddIntegersRust/add-integers/src")

################################################################################
# Define wrapped functions
module AddIntegersRust
using CxxInterface

const libAddIntegersRust = joinpath(pwd(), "libadd_integers_rust")

eval(rustsetup())
eval(rustnewfile("AddIntegersRust/add-integers/src/lib.rs", ""))
eval(rustcode("""
              use half::prelude::*;
              """))

eval(rustfunction(FnName(:add_int, "add_int", libAddIntegersRust), FnResult(Int32, "i32", Int, expr -> :(convert(Int, $expr))),
                  [FnArg(:x, Int32, "x", "i32", Integer, identity), FnArg(:y, Int32, "y", "i32", Integer, identity)],
                  "x + y + f16::from_bits(0).to_f32() as i32"))
end

################################################################################
# Extract Rust code
using CxxInterface
using .AddIntegersRust
AddIntegersRust.rust_write_code!()

################################################################################
# Compile Rust code
# (This fails if there is no Rust compiler available)
using Libdl: dlext
dlpre = Sys.iswindows() ? "" : "lib"
run(Cmd(`cargo build --release`; dir="AddIntegersRust/add-integers"))
cp("AddIntegersRust/add-integers/target/release/$(dlpre)add_integers_rust.$dlext", "libadd_integers_rust.$dlext"; force=true)

# Please, DO NOT call a Rust compiler manually in your own Julia
# packages. This works only in very controlled environments such as on
# CI infrastructure. If you do, your package will be fragile, and will
# create lots of headaches for your users in the wild. Instead, use
# [BinaryBuilder](https://binarybuilder.org) and store your build
# recipes on [Yggdrasil](https://github.com/JuliaPackaging/Yggdrasil).

################################################################################
# Call the wrapped function
# (This fails if the Rust compiler is not compatible with the
# currently running Julia executable)

# Only test cases where the Github CI environment supports this (:i686 is not supported)
if Sys.ARCH ≡ :x86_64
    @test AddIntegersRust.add_int(2, 3) == 5
end