text:
## Wing analysis case
using AeroMDAO

## Surfaces

# Wing
wing = Wing(foils     = fill(naca4(2,4,1,2), 2),
            chords    = [1.0, 0.6],
            twists    = [2.0, 0.0],
            spans     = [8.0],
            dihedrals = [5.],
            sweeps    = [5.]);

x_w, y_w, z_w = wing_mac = mean_aerodynamic_center(wing)
S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing);

# Horizontal tail
htail = Wing(foils     = fill(naca4(0,0,1,2), 2),
             chords    = [0.7, 0.42],
             twists    = [0.0, 0.0],
             spans     = [2.5],
             dihedrals = [0.],
             sweeps    = [6.39],
             position  = [4., 0, -0.1],
             angle     = -2.,
             axis      = [0., 1., 0.])

# Vertical tail
vtail = HalfWing(foils     = fill(naca4(0,0,0,9), 2),
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [1.0],
                 dihedrals = [0.],
                 sweeps    = [7.97],
                 position  = [4., 0, 0],
                 angle     = 90.,
                 axis      = [1., 0., 0.])

# Print info
print_info(wing, "Wing")
print_info(htail, "Horizontal Tail")
print_info(vtail, "Vertical Tail")

## WingMesh type
wing_mesh  = WingMesh(wing, [12], 6,
                      span_spacing = Cosine()
                     )
htail_mesh = WingMesh(htail, [12], 6,
                      span_spacing = Cosine()
                     )
vtail_mesh = WingMesh(vtail, [12], 6,
                      span_spacing = Cosine()
                     )

aircraft =  (
                wing  = make_horseshoes(wing_mesh),
                htail = make_horseshoes(htail_mesh),
                vtail = make_horseshoes(vtail_mesh)
            );

## Case
fs  = Freestream(alpha = 0.0,
                 beta  = 0.0,
                 omega = [0., 0., 0.]);

ref = References(speed     = 1.0,
                 density   = 1.225,
                 viscosity = 1.5e-5,
                 area      = projected_area(wing),
                 span      = span(wing),
                 chord     = mean_aerodynamic_chord(wing),
                 location  = mean_aerodynamic_center(wing))

##
@time begin
    system = solve_case(aircraft, fs, ref;
                        # print            = true, # Prints the results for only the aircraft
                        # print_components = true, # Prints the results for all components
                       );

    # Compute dynamics
    ax       = Geometry() # Geometry, Stability(), Body()
    CFs, CMs = surface_coefficients(system; axes = ax)
    Fs, Ms   = surface_dynamics(system; axes = ax)
    # Fs       = surface_forces(system; axes = ax)
    # vels     = surface_velocities(system)

    nfs = nearfield_coefficients(system)
    ffs = farfield_coefficients(system)
end;

## Viscous drag prediction

# Equivalent flat-plate skin friction estimation
CDv_wing  = profile_drag_coefficient(wing,  [0.8, 0.8], system.reference)
CDv_htail = profile_drag_coefficient(htail, [0.6, 0.6], system.reference)
CDv_vtail = profile_drag_coefficient(vtail, [0.6, 0.6], system.reference)

CDv_plate = CDv_wing + CDv_htail + CDv_vtail

## Local dissipation form factor friction estimation (WRONG???)
import LinearAlgebra: norm

edge_speeds = norm.(surface_velocities(system)); # Inviscid speeds on the surfaces

# Drag coefficients
CDvd_wing  = profile_drag_coefficient(wing_mesh,  [0.8, 0.8], edge_speeds.wing,  ref)
CDvd_htail = profile_drag_coefficient(htail_mesh, [0.6, 0.6], edge_speeds.htail, ref)
CDvd_vtail = profile_drag_coefficient(vtail_mesh, [0.6, 0.6], edge_speeds.vtail, ref)

CDv_diss   = CDvd_wing + CDvd_htail + CDvd_vtail

## Viscous drag coefficient
CDv = CDv_plate

## Total force coefficients with empirical viscous drag prediction
CDi_nf, CY_nf, CL_nf, Cl, Cm, Cn = nf = nearfield(system)
CDi_ff, CY_ff, CL_ff = ff = farfield(system)

nf_v = (CD = CDi_nf + CDv, CDv = CDv, nf...)
ff_v = (CD = CDi_ff + CDv, CDv = CDv, ff...)

## Spanwise forces/lifting line loads
wing_ll  = spanwise_loading(chord_panels(wing_mesh),  CFs.wing,  S)
htail_ll = spanwise_loading(chord_panels(htail_mesh), CFs.htail, S)
vtail_ll = spanwise_loading(chord_panels(vtail_mesh), CFs.vtail, S);

## Plotting
using CairoMakie
CairoMakie.activate!()

set_theme!(
            # theme_black()
            # theme_light()
          )

using LaTeXStrings
const LS = LaTeXString

## Streamlines
# Spanwise distribution
span_points = 30
init        = chop_leading_edge(wing, span_points)
dx, dy, dz  = 0, 0, 1e-1
seed        = init # [ init .+ Ref([dx, dy,  dz])
                   #   init .+ Ref([dx, dy, -dz]) ];

distance = 5
num_stream_points = 100
streams = streamlines(system, seed, distance, num_stream_points);

wing_cam_connec  = triangle_connectivities(LinearIndices(wing_mesh.camber_mesh))
htail_cam_connec = triangle_connectivities(LinearIndices(htail_mesh.camber_mesh))
vtail_cam_connec = triangle_connectivities(LinearIndices(vtail_mesh.camber_mesh));

## Surface velocities
vels = surface_velocities(system);
sps  = norm.(vels)

wing_sp_points  = extrapolate_point_mesh(sps.wing)
htail_sp_points = extrapolate_point_mesh(sps.htail)
vtail_sp_points = extrapolate_point_mesh(sps.vtail)

## Surface pressure coefficients
cps  = norm.(CFs) * S

wing_cp_points  = extrapolate_point_mesh(cps.wing)
htail_cp_points = extrapolate_point_mesh(cps.htail)
vtail_cp_points = extrapolate_point_mesh(cps.vtail)

## Figure plot
fig1  = Figure(resolution = (1280, 720))

scene = LScene(fig1[1:4,1])
ax    = fig1[1:4,2] = GridLayout()

ax_cd = Axis(ax[1,1:2], ylabel = L"C_{D_i}", title = LS("Spanwise Loading"))
ax_cy = Axis(ax[2,1:2], ylabel = L"C_Y",)
ax_cl = Axis(ax[3,1:2], xlabel = L"y", ylabel = L"C_L")

# Spanload plot
function plot_spanload!(ax, ll_loads, name = "Wing")
    @views lines!(ax[1,1:2], ll_loads[:,1], ll_loads[:,2], label = name,)
    @views lines!(ax[2,1:2], ll_loads[:,1], ll_loads[:,3], label = name,)
    @views lines!(ax[3,1:2], ll_loads[:,1], ll_loads[:,4], label = name,)

    nothing
end

plot_spanload!(ax, wing_ll, LS("Wing"))
plot_spanload!(ax, htail_ll, LS("Horizontal Tail"))
plot_spanload!(ax, vtail_ll, LS("Vertical Tail"))

# Legend
Legend(ax[4,1:2], ax_cl)

# Surface pressure meshes
m1 = poly!(scene, vec(wing_mesh.camber_mesh),  wing_cam_connec,  color = vec(wing_cp_points))
m2 = poly!(scene, vec(htail_mesh.camber_mesh), htail_cam_connec, color = vec(htail_cp_points))
m3 = poly!(scene, vec(vtail_mesh.camber_mesh), vtail_cam_connec, color = vec(vtail_cp_points))

Colorbar(fig1[4,1], m1.plots[1], label = L"Pressure Coefficient, $C_p$", vertical = false)

fig1[0, :] = Label(fig1, LS("Vortex Lattice Analysis"), textsize = 20)

# Airfoil meshes
# wing_surf = surface_coordinates(wing_mesh, wing_mesh.n_span, 60)
# surf_connec = triangle_connectivities(LinearIndices(wing_surf))
# wing_surf_mesh = mesh(vec(wing_surf), surf_connec)
# w1 = wireframe!(scene, wing_surf_mesh.plot[1][], color = :grey, alpha = 0.1)

# Borders
lines!(scene, plot_planform(wing))
lines!(scene, plot_planform(htail))
lines!(scene, plot_planform(vtail))

# l1 = [ lines!(scene, pts, color = :grey) for pts in plot_panels(camber_panels(wing_mesh))  ]
# l2 = [ lines!(scene, pts, color = :grey) for pts in plot_panels(camber_panels(htail_mesh)) ]
# l3 = [ lines!(scene, pts, color = :grey) for pts in plot_panels(camber_panels(vtail_mesh)) ]

# Streamlines
[ lines!(scene, stream[:], color = :green) for stream in eachcol(streams) ]

fig1.scene

## Save figure
# save("plots/VortexLattice.pdf", fig1, px_per_unit = 1.5)

## Animation settings
pts = [ Node(Point3f0[stream]) for stream in streams[1,:] ]

[ lines!(scene, pts[i], color = :green, axis = (; type = Axis3)) for i in eachindex(pts) ]

# Recording
fps     = 30
nframes = length(streams[:,1])

# record(fig, "plots/vlm_animation.mp4", 1:nframes) do i
#     for j in eachindex(streams[1,:])
#         pts[j][] = push!(pts[j][], Point3f0(streams[i,j]))
#     end
#     sleep(1/fps) # refreshes the display!
#     notify(pts[i])
# end

## Arrows
# hs_pts = vec(Tuple.(bound_leg_center.(horses)))
# arrows!(scene, getindex.(hs_pts, 1), getindex.(hs_pts, 2), getindex.(hs_pts, 3),
#                 vec(CDis),vec(CYs),vec( Ls),
#                 arrowsize = Vec3f.(0.3, 0.3, 0.4),
#                 lengthscale = 10,
#                 label = "Forces (Exaggerated)")

## VARIABLE ANALYSES
#=========================================================#

using Setfield
using Base.Iterators: product
using Plots

pyplot()

## Speed sweep
Vs = 1.0:10:300
res_Vs = permutedims(combinedimsview(
    map(Vs) do V
        ref1 = @set ref.speed = V
        sys = solve_case(aircraft, fs, ref1)
        [ V; farfield(sys)...; nearfield(sys)... ]
    end
))

##
Plots.plot(
    res_Vs[:,1], res_Vs[:,2:end],
    layout = (3,3), size = (900, 800),
    xlabel = "V",
    labels = ["CD_ff" "CY_ff" "CL_ff" "CD" "CY" "CL" "Cl" "Cm" "Cn"]
)

##
f2 = Figure()
ax1 = Axis(f2[1,1], xlabel = L"V", ylabel = L"C_{D_{ff}}")
lines!(res_Vs[:,1], res_Vs[:,2]) # CD_ff
ax2 = Axis(f2[1,2], xlabel = L"V", ylabel = L"C_{Y_{ff}}")
lines!(res_Vs[:,1], res_Vs[:,3])
ax3 = Axis(f2[1,3], xlabel = L"V", ylabel = L"CL_{ff}")
lines!(res_Vs[:,1], res_Vs[:,4])
ax4 = Axis(f2[2,1], xlabel = L"V", ylabel = L"C_D")
lines!(res_Vs[:,1], res_Vs[:,5])
ax5 = Axis(f2[2,2], xlabel = L"V", ylabel = L"C_Y")
lines!(res_Vs[:,1], res_Vs[:,6])
ax6 = Axis(f2[2,3], xlabel = L"V", ylabel = L"C_L")
lines!(res_Vs[:,1], res_Vs[:,7])
ax7 = Axis(f2[3,1], xlabel = L"V", ylabel = L"C_\ell")
lines!(res_Vs[:,1], res_Vs[:,8])
ax8 = Axis(f2[3,2], xlabel = L"V", ylabel = L"C_m")
lines!(res_Vs[:,1], res_Vs[:,9])
ax9 = Axis(f2[3,3], xlabel = L"V", ylabel = L"C_n")
lines!(res_Vs[:,1], res_Vs[:,10])

f2

## Alpha sweep
αs = -5:0.5:5
res_αs = permutedims(combinedimsview(
    map(αs) do α
        fst = @set fs.alpha = deg2rad(α)
        sys = solve_case(aircraft, fst, ref)
        [ α; farfield(sys)...; nearfield(sys)... ]
    end
))

Plots.plot(
    res_αs[:,1], res_αs[:,2:end],
    layout = (3,3), size = (900, 800),
    xlabel = "α",
    labels = ["CD_ff" "CY_ff" "CL_ff" "CD" "CY" "CL" "Cl" "Cm" "Cn"]
)

## (Speed, alpha) sweep
res = combinedimsview(
    map(product(Vs, αs)) do (V, α)
        ref1 = @set ref.speed = V
        fst = @set fs.alpha = deg2rad(α)
        sys = solve_case(aircraft, fst, ref1)
        [ V; α; farfield(sys)...; nearfield(sys)... ]
    end
)

##
res_p = permutedims(res, (3,1,2))

# CDi
plt_CDi_ff = Plots.plot(camera = (75, 30), ylabel = L"\alpha", xlabel = "V", zlabel = L"C_{D_i}")
[ Plots.plot!(res_p[:,1,n], res_p[:,2,n], res_p[:,3,n], label = "", c = :black) for n in axes(res_p,3) ]

# CL
plt_CL_ff = Plots.plot(camera = (75,15),
ylabel = L"\alpha", xlabel = "V", zlabel = L"C_{L}",
label = "", c = :black)
[ Plots.plot!(
    res_p[:,1,n], res_p[:,2,n], res_p[:,8,n], label = "", c = :black) for n in axes(res_p,3) ]

plt_Cm_ff = Plots.plot(camera = (100,30), ylabel = L"\alpha", xlabel = "V", zlabel = L"C_m", label = "")
[ Plots.plot!(
    res_p[:,1,n], res_p[:,2,n], res_p[:,10,n], label = "", c = :black) for n in axes(res_p,3) ]

##
p = Plots.plot(plt_CDi_ff, plt_CL_ff, plt_Cm_ff, layout = (1,3), size = (1300, 400))

##
savefig(p, "coeffs.png")


text:
# BGZF blocks
#
# Code in this file should not "know about" Codecs, TranscodingStreams,
# or any of that. It should only rely on Base and LibDeflate, so that means the
# code in this file can easily be cannibalized for other packages, or repurposed.

# We must write these bytes *exactly* per the BGZF spec,
# we can't just write an empty block, as some of the fields
# (e.g. MTIME or OS) may not be these exact bytes.
const EOF_BLOCK = [
    0x1f, 0x8b, 0x08, 0x04, 0x00, 0x00, 0x00, 0x00,
    0x00, 0xff, 0x06, 0x00, 0x42, 0x43, 0x02, 0x00,
    0x1b, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00
]

# BGZF blocks are no larger than 64 KiB before and after compression.
const MAX_BLOCK_SIZE = UInt(64 * 1024)

# Maximum number of bytes to be compressed at one time. Random bytes usually end up filling
# a bit more when compressed, so we have a generous 256 byte margin of safety.
const SAFE_BLOCK_SIZE = UInt(MAX_BLOCK_SIZE - 256)

# Field descriptions are for decompressors / compressors
mutable struct Block{T}
    de_compressor::T
    outdata::Vector{UInt8}
    indata::Vector{UInt8}
    task::Task
    gzip_extra_fields::Vector{LibDeflate.GzipExtraField} # cached to avoid allocating
    crc32::UInt32    # stated checksum / calculated checksum

    # BGZF blocks can store 0:typemax(UInt16)+1 bytes
    # so unfortunately UInt16 will not suffice here.
    outlen::UInt32   # Length of decompressed payload / compressed block
    inlen::UInt32    # Length of compressed payload / total input block
end

function Block(dc::T) where T <: DE_COMPRESSOR
    outdata = Vector{UInt8}(undef, MAX_BLOCK_SIZE)
    indata = similar(outdata)

    # We initialize with a trivial, but completable task for sake of simplicity
    task = schedule(Task(() -> nothing))
    return Block{T}(dc, outdata, indata, task, LibDeflate.GzipExtraField[], 0, 0, 0)
end

function Base.empty!(block::Block)
    block.outlen = 0
    block.inlen = 0
end
Base.isempty(block::Block) = iszero(block.outlen)

nfull(::Type{Block{Decompressor}}) = MAX_BLOCK_SIZE
nfull(::Type{Block{Compressor}}) = SAFE_BLOCK_SIZE
Base.wait(b::Block) = wait(b.task)

function check_eof_block(block::Block{Decompressor})
    # We don't store the metadata, so can't compare to EOF_BLOCK. BGZF specs
    # allows us to simply check for a zero-length block.
    if !iszero(block.outlen)
        bgzferror("No EOF block. Truncated file?")
    end
end

"Get the BSIZE field from block by reading from given vector"
function bsize(block::Block{Decompressor}, vector::Vector{UInt8})::Union{UInt16, Nothing}
    fieldnum = findfirst(block.gzip_extra_fields) do field
        field.tag === (UInt8('B'), UInt8('C'))
    end
    fieldnum === nothing && return nothing
    field = @inbounds block.gzip_extra_fields[fieldnum]
    field.data === nothing && return nothing
    length(field.data) != 2 && return nothing
    return (vector[first(field.data)] % UInt16) | ((vector[last(field.data)] % UInt16) << 8)
end

"""Load data from buffer starting at index 1 into the block, and queue its decompression
in another thread. Return number of bytes consumed"""
function load_block!(block::Block{Decompressor}, buffer::Vector{UInt8}, len::Unsigned)
    # Parse header
    header_len, header = LibDeflate.parse_gzip_header(
        buffer,
        len % UInt,
        block.gzip_extra_fields
    )

    bsiz = bsize(block, buffer)
    bsiz === nothing && throw(BGZFError(NO_GZIP_EXTRA_FIELD_BSIZE))
    # By spec, BSIZE is block size -1. Include header_len bytes header, 8 byte tail
    blocksize = bsiz + 1
    block.inlen = blocksize - header_len - 8
    block.crc32 = bitload(UInt32, buffer, blocksize - 7)
    block.outlen = bitload(UInt32, buffer, blocksize - 3)
    unsafe_copyto!(block.indata, 1, buffer, header_len + 1, block.inlen)
    queue!(block)
    return blocksize
end

"Load `len_minus_one+1` bytes from `buffer` into the block and queue its compression in another thread."
function load_block!(block::Block{Compressor}, buffer::Vector{UInt8}, len_minus_one::UInt16)
    block.inlen = len_minus_one + 1
    # We zero the outlen to mark that it has no data _until_ the queue happens,
    # where the correct outlen is set
    block.outlen = 0
    unsafe_copyto!(block.indata, 1, buffer, 1, len_minus_one + 1)
    queue!(block)
    return nothing
end

"Process the block in another thread"
queue!(block::Block) = block.task = @spawn _queue!(block)

function _queue!(block::Block{Decompressor})
    unsafe_decompress!(Base.HasLength(), block.de_compressor,
                       pointer(block.outdata), block.outlen,
                       pointer(block.indata), block.inlen)

    crc32 = unsafe_crc32(pointer(block.outdata), block.outlen)
    crc32 != block.crc32 && bgzferror("CRC32 checksum does not match")
    return nothing
end

function _queue!(block::Block{Compressor})
    # Meat: The compressed data
    compress_len = unsafe_compress!(block.de_compressor,
                   pointer(block.outdata, 19), MAX_BLOCK_SIZE - 26,
                   pointer(block.indata), block.inlen)
    block.crc32 = unsafe_crc32(pointer(block.indata), block.inlen)
    block.outlen = compress_len + 26

    # Header: 18 bytes of header
    unsafe_copyto!(block.outdata, 1, BLOCK_HEADER, 1, 16)
    bitstore(UInt16(block.outlen - 1), block.outdata, 17)

    # Tail: CRC + isize
    bitstore(block.crc32, block.outdata, 18 + compress_len + 1)
    bitstore(block.inlen % UInt32, block.outdata, 18 + compress_len + 5)
    return nothing
end

ext: jl


text:
module OffsetBuffers

export
    OffsetBuffer,
    RollingBuffer,
    SlidingBuffer,
    capacity, # counfirst, count, offset,
    isfull,
    prepare_append!,
    reset!

using CustomUnitRanges
include(CustomUnitRanges.filename_for_urange)  # defines URange

include("rolling_buffer.jl")
include("sliding_buffer.jl")

const OffsetBuffer{T} = Union{RollingBuffer{T}, SlidingBuffer{T}}

end # module


text:
#@nospecialize

function determine_variables(expr)::Set{Symbol}
    variables = Set{Symbol}()
    ignored_symbols = Set{Symbol}()
    push!(ignored_symbols, :target)
    prewalk(expr) do x
        if x isa Symbol
            x ∈ ignored_symbols || push!(variables, x)
        elseif x isa Expr
            if x.head === :(=)
                LHS = first(x.args)
                if LHS isa Symbol
                    push!(ignored_symbols, LHS)
                elseif LHS isa Expr && LHS.head === :ref
                    push!(ignored_symbols, first(LHS.args))
                end
            elseif x.head === :call
                push!(ignored_symbols, first(x.args))
            end
        end
        return x
    end
    variables
end

function no_fmadd_dist_call(y, dargs)
    dist = Expr(:call, first(dargs), y)
    for i ∈ 2:length(dargs)
        push!(dist.args, dargs[i])
    end
    Expr(:(=), :target, Expr(:call, Expr(:(.), :ProbabilityModels, QuoteNode(:vadd)), :target, dist))
end
function fmadd_dist_call(y, dargs)
    farg₁ = dargs[2]
    # farg₁ must be an expr of the form (X * β + α)
    (farg₁ isa Expr && farg₁.head === :call) || return no_fmadd_dist_call(y, dargs)
    args = farg₁.args
    Nargs = length(args)
    Nargs > 2 || return no_fmadd_dist_call(y, dargs)
    if first(args) !== :+
        if first(dargs) === :Normal && first(args) === :* && length(args) == 3
            dist = Expr(:call, :Normal, y, args[2], args[3])
            for i ∈ 3:length(dargs)
                push!(dist.args, dargs[i])
            end
            Expr(:(=), :target, Expr(:call, Expr(:(.), :ProbabilityModels, QuoteNode(:vadd)), :target, dist))
        else
            return no_fmadd_dist_call(y, dargs)
        end
    end
    j = 2
    while j ≤ Nargs
        argsⱼ = args[j]
        if argsⱼ isa Expr && (first(argsⱼ.args) === :*) && length(argsⱼ.args) == 3
            break
        end
        j += 1
    end
    j > Nargs && return no_fmadd_dist_call(y, dargs)
    # if we're here, j is the index of the mul
    argsⱼ = args[j].args
    dist = Expr(:call, first(dargs), y, argsⱼ[2], argsⱼ[3])
    if Nargs == 3
        push!(dist.args, args[j == 2 ? 3 : 2])
    else #Nargs > 3
        for i ∈ 2:Nargs
            i == j || push!(dist.args, args[i])
        end
    end
    for i ∈ 3:length(dargs)
        push!(dist.args, dargs[i])
    end
    Expr(:(=), :target, Expr(:call, Expr(:(.), :ProbabilityModels, QuoteNode(:vadd)), :target, dist))
end
function dist_call(x)
    distcall::Expr = x.args[3]
    distcall.head === :call || return x
    length(distcall.args) < 2 && return x
    y = x.args[2]
    dargs = distcall.args
    if first(dargs) ∈ ProbabilityDistributions.FMADD_DISTRIBUTIONS
        fmadd_dist_call(y, dargs)
    else
        no_fmadd_dist_call(y, dargs)
    end
end

function translate_sampling_statement(x)
    x isa Expr || return x
    if x.head === :(+=)
        return Expr(:(=), x.args[1], Expr(:+, x.args[1], x.args[2]))
    elseif x.head === :call && first(x.args) === :~
        dist_call(x)
    else
        return x
    end
end

translate_sampling_statements(expr::Expr)::Expr = prewalk(translate_sampling_statement, expr)

# This does not work, because eval only looks in the current module, and doesn't accept a module argument (ie, Main)
# function interpolate_globals(expr)::Expr
#     postwalk(expr) do x
#         if x isa Expr && x.head == :$
#             return eval(first(x.args))
#         else
#             return x
#         end
#     end
# end
maybeget(d::Dict{K,V}, ex) where {K,V} = ex isa K ? get(d, ex, ex) : ex
ssa_sym(i::Int) = Symbol("##SSAValue##$(i)##")
function ssa_to_sym(expr)
    postwalk(expr) do ex
        if ex isa Core.SSAValue
            return ssa_sym(ex.id)
        elseif ex isa GlobalRef
            return Expr(:., Symbol(ex.mod), QuoteNode(ex.name))
        else
            return ex
        end
    end
end
function subssasyms!(q, ex, substitutions, i)
    ex = postwalk(ex) do x
        y = ssa_to_sym(x)
        maybeget(substitutions, y)
    end
    push!(q.args, :($(ssa_sym(i)) = $ex))
    nothing
end
function flatten_expression(expr)::Expr
    lowered_array = Meta.lower(ProbabilityModels, expr).args
    substitutions = Dict{Symbol,Expr}()
    assignments = Dict{Symbol,Symbol}()
    @assert length(lowered_array) == 1
    lowered = first(lowered_array).code
    q = quote end
    for i ∈ 1:length(lowered) - 1 # skip return statement
        ex = lowered[i]
        if ex isa Expr
            if ex.head === :call
                if first(ex.args) === :(:)
                    push!(q.args, :($(ssa_sym(i)) = ProbabilityModels.PaddedMatrices.Static{($(ex.args[2]),$(ex.args[3]))}()))
                elseif ex.args[1] isa GlobalRef && ex.args[1].name == :getproperty
                    substitutions[ssa_sym(i)] = Expr(:., ex.args[2], ex.args[3])
                else
                    subssasyms!(q, ex, substitutions, i)
                end
            elseif ex.head === :(=)
                assignments[ssa_to_sym(ex.args[2])] = ex.args[1]
            else
                subssasyms!(q, ex, substitutions, i)
            end
        else
            subssasyms!(q, ex, substitutions, i)
        end
    end
    q = postwalk(q) do ex
        maybeget(assignments, ex)
    end
    # println(q)
    q
end

function drop_const_assignment!(first_pass, tracked_vars, ex, LHS, RHS, verbose)
    if RHS isa Expr && RHS.head === :call
        f = first(RHS.args)
        A = @view(RHS.args[2:end])
        if f ∈ ProbabilityDistributions.DISTRIBUTION_DIFF_RULES
            track_tup = Expr(:tuple,)
            for a ∈ A
                if a ∈ tracked_vars
                    push!(track_tup.args, true)
                    push!(tracked_vars, LHS)
                else
                    push!(track_tup.args, false)
                end
            end
            if verbose
                printstring = "distribution $f (ret: $LHS): "
                push!(first_pass, :(println($printstring)))
            end
            push!(first_pass, :($LHS = ProbabilityModels.ProbabilityDistributions.$f(Val{$track_tup}(), $(A...))))
            verbose && push!(first_pass, :(println($LHS)))
        else
            for a ∈ A
                if a ∈ tracked_vars
                    push!(tracked_vars, LHS)
                    break
                end
            end
            if f === :add
                push!(first_pass, :( $LHS = ProbabilityModels.vadd($(A...))))
            else
                push!(first_pass, ex)
            end
        end
    else
        A = RHS.args
        for a ∈ A
            if a ∈ tracked_vars
                push!(tracked_vars, out)
                break
            end
        end
        push!(first_pass, ex)
    end
end

"""
This pass is for when we aren't taking partial derivatives.
"""
function constant_drop_pass!(first_pass::Vector{Any}, expr, tracked_vars, verbose = false)
    for x ∈ expr.args
        if !(x isa Expr)
            push!(first_pass, x)
            continue
        end
        ex::Expr = x
        if ex.head === :(=)
            drop_const_assignment!(first_pass, tracked_vars, ex, ex.args[1], ex.args[2], verbose)
        else
            push!(first_pass, ex)
        end
    end
    nothing
end

types_to_vals(::Type{T}) where {T} = Val{T}()
types_to_vals(v) = v
types_to_vals(A::AbstractArray{Union{Missing,T}}) where {T} = DistributionParameters.maybe_missing(A)
extract_typeval(::Type{Val{T}}) where {T} = T
extract_typeval(::Type{Val{Type{T}}}) where {T} = T

function load_and_constrain_quote(ℓ, model_name, variables, variable_type_names, θ, Tθ, T)
    Nparam = gensym(:number_parameters)
    logjac = gensym(:logjac)
    θq = quote
        @generated function ProbabilityModels.DistributionParameters.constrain($ℓ::$(model_name){$Nparam, $(variable_type_names...)}, $θ::AbstractVector{$T}, ::Val{$logjac} = Val{false}()) where {$Nparam, $T, $logjac, $(variable_type_names...)}
            return_partials = false
            first_pass = quote end
            push!(first_pass.args, Expr(:(=), Symbol("##θparameter##"), Expr(:call, :(ProbabilityModels.vectorizable), $(QuoteNode(θ)))))
            second_pass = quote end
            transformed_params = Expr(:tuple,)
            return_expr = $logjac ? Expr(:tuple, transformed_params, :target) : transformed_params
        end
    end
    v = Symbol("##storage_vector##")
    θq_vec = quote
        @generated function ProbabilityModels.DistributionParameters.constrain!(
            $v::AbstractVector{$T}, $ℓ::$(model_name){$Nparam, $(variable_type_names...)}, $θ::AbstractVector{$T}
        ) where {$Nparam, $T, $(variable_type_names...)}
            return_partials = false
            first_pass = quote end
            push!(first_pass.args, Expr(:(=), Symbol("##stack_pointer##"), Expr(:call, :(ProbabilityModels.StackPointer), Expr(:call, :pointer, Symbol("##storage_vector##")))))
            push!(first_pass.args, Expr(:(=), Symbol("##θparameter##"), Expr(:call, :(ProbabilityModels.vectorizable), $(QuoteNode(θ)))))
            second_pass = quote end
        end
    end
    param_names = gensym(:param_names)
    param_names_quote = quote
        function ProbabilityModels.DistributionParameters.parameter_names($ℓ::$(model_name){$Nparam, $(variable_type_names...)}) where {$Nparam, $(variable_type_names...)}
            $param_names = String[]
            $([quote
               if $v <: Val
               # @show $v
               append!($param_names, ProbabilityModels.DistributionParameters.parameter_names(ProbabilityModels.extract_typeval($v), $(QuoteNode(v)))::Vector{String})
               elseif $v <: ProbabilityModels.DistributionParameters.MissingDataArray
               #append!($param_names, ProbabilityModels.DistributionParameters.parameter_names($v, $(QuoteNode(v)))::Vector{String})
               append!($param_names, ProbabilityModels.DistributionParameters.parameter_names($ℓ.$(variables[i]), $(QuoteNode(v)))::Vector{String})
               end
               end for (i,v) ∈ enumerate(variable_type_names)]...)
            $param_names
        end
    end
    Nconstrainedparam = Symbol("##number_of_parameters##")
    constrained_length_quote = quote
        @generated function ProbabilityModels.DistributionParameters.constrained_length($ℓ::$(model_name){$Nparam, $(variable_type_names...)}) where {$Nparam, $(variable_type_names...)}
            $Nconstrainedparam = 0
            $([quote
               if $v <: Val
                   $Nconstrainedparam += ProbabilityModels.PaddedMatrices.type_length(ProbabilityModels.extract_typeval($v))
               elseif $v <: ProbabilityModels.DistributionParameters.MissingDataArray
                   $Nconstrainedparam += ProbabilityModels.PaddedMatrices.type_length($v)
               end
               end for v ∈ variable_type_names]...)
            $Nconstrainedparam
        end
    end
    θq_body = θq.args[end].args[end].args[end].args;
    θq_vec_body = θq_vec.args[end].args[end].args[end].args;
    for i ∈ eachindex(variables)
        load_data = Expr(:quote, :($(variables[i]) = $ℓ.$(variables[i])))
        load_incomplete_data = Expr(:quote, Expr(:(=), Symbol("##incomplete##", variables[i]), :($ℓ.$(variables[i]))))
        push!(θq_body, quote
          if $(variable_type_names[i]) <: Val
              ProbabilityModels.DistributionParameters.load_parameter!(
                  first_pass.args, second_pass.args, $(QuoteNode(variables[i])),
                  ProbabilityModels.extract_typeval($(variable_type_names[i])),
                  false, ProbabilityModels, nothing, $logjac, true
              )
              push!(transformed_params.args, Expr(:(=), $(QuoteNode(variables[i])), $(QuoteNode(variables[i]))))
          elseif $(variable_type_names[i]) <: ProbabilityModels.DistributionParameters.MissingDataArray
              push!(first_pass.args, $load_incomplete_data)
              ProbabilityModels.DistributionParameters.load_missing_as_vector!(
                  first_pass.args, second_pass.args, $(QuoteNode(variables[i])), ($(variable_type_names[i])),
                  false, ProbabilityModels, nothing, $logjac, true
              )
              push!(transformed_params.args, Expr(:(=), $(QuoteNode(variables[i])), $(QuoteNode(variables[i]))))
          end
        end)
        push!(θq_vec_body, quote
            if $(variable_type_names[i]) <: Val
                ProbabilityModels.DistributionParameters.load_parameter!(first_pass.args, second_pass.args, $(QuoteNode(variables[i])), ProbabilityModels.extract_typeval($(variable_type_names[i])), false, ProbabilityModels, Symbol("##stack_pointer##"), false, true)
            #     if ProbabilityModels.extract_typeval($(variable_type_names[i])) <: ProbabilityModels.DistributionParameters.RealFloat
            #         push!(first_pass.args, Expr(:call, :(ProbabilityModels.store!),
            #                               Expr(:call, :pointer, Symbol("##stack_pointer##"), $(QuoteNode(T))), $(QuoteNode(variables[i]))))
            #         push!(first_pass.args, Expr(:(+=), Symbol("##stack_pointer##"), Expr(:call, :sizeof, $(QuoteNode(T)) )))
            #   end
              elseif $(variable_type_names[i]) <: ProbabilityModels.DistributionParameters.MissingDataArray
                  push!(first_pass.args, $load_incomplete_data)
                  ProbabilityModels.DistributionParameters.load_missing_as_vector!(
                      first_pass.args, second_pass.args, $(QuoteNode(variables[i])), ($(variable_type_names[i])),
                      false, ProbabilityModels, Symbol("##stack_pointer##"), false, true
                  )
            end
       end)
    end
    push!(θq_body, :(push!(first_pass.args, return_expr)))
    push!(θq_vec_body, Expr(:call, :push!, :(first_pass.args), QuoteNode(Symbol("##storage_vector##"))))
    push!(θq_body, :(first_pass))
    # push!(θq_body, :(QuoteNode(first_pass)))
    push!(θq_vec_body, :(first_pass))
    θq, θq_vec, param_names_quote, constrained_length_quote
end

function generate_generated_funcs_expressions(model_name, expr)
    # Determine the set of variables that are either parameters or data, after interpolating any globals inserted via `$`
    # expr = interpolate_globals(expr)]
    variable_set = determine_variables(expr)
    variables = sort!([v for v ∈ variable_set]) # ensure order is constant
    variable_type_names = [Symbol("##Type##", v) for v ∈ variables]

    Nparam = Symbol("##number_of_parameters##")
    stride = gensym(:LDA)
    L1 = gensym(:L)
    L2 = gensym(:L)
    T = Symbol("##element_type##")

    var_vartype_pairs = [:( $(variables[i])::$(variable_type_names[i]) ) for i ∈ eachindex(variables)]
    struct_quote = quote
        struct $model_name{$Nparam, $(variable_type_names...)} <: ProbabilityModels.AbstractProbabilityModel{$Nparam}
            $(var_vartype_pairs...)
        end
    end
    precomp = gensym(:precompile)
    # precomp_quote = if VERSION >= v"1.3.0-alpha.0" && ProbabilityModels.NTHREADS[] > 1
    #     tlag = gensym(:tlag); tl = gensym(:tl)
    #     quote
    #         $tlag = Threads.@spawn precompile(ProbabilityModels.logdensity_and_gradient!, (ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, $model_name{$Nparam,$(variable_type_names...)}, ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, ProbabilityModels.StackPointer))
    #         $tl = Threads.@spawn precompile(ProbabilityModels.logdensity, ($model_name{$Nparam,$(variable_type_names...)}, ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, ProbabilityModels.StackPointer))
    #         Threads.@spawn precompile(ProbabilityModels.InplaceDHMC.mcmc_with_warmup!, (
    #             ProbabilityModels.VectorizedRNG.PtrPCG{4}, ProbabilityModels.StackPointer, ProbabilityModels.PaddedMatrices.DynamicPtrVector{Float64}, ProbabilityModels.PaddedMatrices.DynamicPtrVector{ProbabilityModels.InplaceDHMC.TreeStatisticsNUTS}, $model_name{$Nparam, $(variable_type_names...)}, Int))
    #         wait($tlag); wait($tl)
    #     end
    # else
    #     quote
    #         precompile(ProbabilityModels.logdensity_and_gradient!, (ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, $model_name{$Nparam,$(variable_type_names...)}, ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, ProbabilityModels.StackPointer))
    #         precompile(ProbabilityModels.logdensity, ($model_name{$Nparam,$(variable_type_names...)}, ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, ProbabilityModels.StackPointer))
    #     end
    # end
    precomp_quote = quote
        precompile(ProbabilityModels.logdensity_and_gradient!, (ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, $model_name{$Nparam,$(variable_type_names...)}, ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, ProbabilityModels.StackPointer))
        precompile(ProbabilityModels.logdensity, ($model_name{$Nparam,$(variable_type_names...)}, ProbabilityModels.PtrVector{$Nparam, Float64, $Nparam, false}, ProbabilityModels.StackPointer))
    end
    struct_kwarg_quote = quote
        function $model_name{$Nparam}( $(var_vartype_pairs...)) where {$Nparam, $(variable_type_names...)}
            $model_name{$Nparam,$(variable_type_names...)}($(variables...))
        end
        function $model_name{$Nparam}($precomp::Bool, $(var_vartype_pairs...)) where {$Nparam, $(variable_type_names...)}
            if $precomp
                $precomp_quote
            end
            $model_name{$Nparam,$(variable_type_names...)}($(variables...))
        end
        function $model_name($precomp::Bool = true; $(var_vartype_pairs...)) where {$(variable_type_names...)}
            $Nparam = 0
            # @show $(variables...)
            # @show $(variable_type_names...)
            $([quote
                if isa(ProbabilityModels.types_to_vals($v), Val) || isa($v, ProbabilityModels.DistributionParameters.MissingDataArray)
                   $Nparam += ProbabilityModels.PaddedMatrices.param_type_length($v)
               elseif isa($v, AbstractArray{Union{Missing,T}} where T)
                    $v = ProbabilityModels.DistributionParameters.maybe_missing($v)
                   $Nparam += ProbabilityModels.PaddedMatrices.param_type_length($v)
                end
            end for v ∈ variables]...)
            $model_name{$Nparam}($precomp, $([:(ProbabilityModels.types_to_vals($v)) for v ∈ variables]...))
        end
    end
    # Translate the sampling statements, and then flatten the expression to remove nesting.
    expr = translate_sampling_statements(expr) |>
                flatten_expression# |>
    # The plan in this definition is to make each keyword arg default to the appropriate field of ℓ
    # This allows one to optionally override in the case of a single evaluation.
    # Additionally, it will throw an error if one of the required fields was not defined.
    V = gensym(:V)
    ℓ = Symbol("####ℓ#data####")# gensym(:ℓ)
    θ = Symbol("##θ_parameter_vector##")#gensym(:θ)
    Tθ = gensym(:Tθ)
    constrain_quote, cvq, pn, clq  = load_and_constrain_quote(ℓ, model_name, variables, variable_type_names, θ, Tθ, T)
    base_stack_pointer = STACK_POINTER_REF[]# + 9VectorizationBase.REGISTER_SIZE
    stack_pointer_expr = NTHREADS[] == 1 ? base_stack_pointer : :($base_stack_pointer + (Threads.threadid()-1)*$(LOCAL_STACK_SIZE[]))
    # we have to split these, because of dispatch ambiguity errors
    θq_value = quote
        @generated function ProbabilityModels.logdensity(
                    $ℓ::$(model_name){$Nparam, $(variable_type_names...)},
                    $θ::ProbabilityModels.PtrVector{$Nparam, $T, $Nparam, false},
                    $(Symbol("##stack_pointer##"))::ProbabilityModels.StackPointer = $stack_pointer_expr
        ) where {$Nparam, $T, $(variable_type_names...)}

            return_partials = false
            TLθ = $Nparam
            model_parameters = Symbol[]
            first_pass = Any[]
            second_pass = Any[]
        end
    end
    θq_valuegradient = quote
        @generated function ProbabilityModels.logdensity_and_gradient!(
                        $(Symbol("##∂θ_parameter_vector##"))::ProbabilityModels.PtrVector{$Nparam, $T, $Nparam, false},
                        $ℓ::$(model_name){$Nparam, $(variable_type_names...)},
                        $θ::ProbabilityModels.PtrVector{$Nparam, $T, $Nparam, false},
                        $(Symbol("##initial_stack_pointer##"))::ProbabilityModels.StackPointer = $stack_pointer_expr
        ) where {$Nparam, $T, $(variable_type_names...)}

            return_partials = true
            TLθ = $Nparam
            model_parameters = Symbol[]
            first_pass = Any[]
            second_pass = Any[]
        end
    end
    for (θq, return_partials) ∈ ((θq_value, false), (θq_valuegradient, true))
        θq_body = θq.args[end].args[end].args[end].args;
        push!(θq_body, :(expr = $(Expr(:quote, deepcopy(expr)))))
        # Now, we work on assembling the function.
        for i ∈ eachindex(variables)
            load_data = Expr(:quote, :($(variables[i]) = $ℓ.$(variables[i])))
              load_data_dynamic_ptr_array = Expr(:quote,
                                         :($(variables[i]) = ProbabilityModels.PaddedMatrices.DynamicPtrArray(
                                             pointer( $ℓ.$(variables[i])),
                                             size($ℓ.$(variables[i])),
                                             size($ℓ.$(variables[i]),1)) )
                                         )
            load_data_ptr_array = Expr(:quote, :($(variables[i]) = ProbabilityModels.PaddedMatrices.PtrArray( $ℓ.$(variables[i]) ) ) )
            missingvar = Symbol("##missing##", variables[i])
            load_incomplete_data = Expr(:quote, Expr(:(=), Symbol("##incomplete##", variables[i]), :($ℓ.$(variables[i]))))
            push!(θq_body, quote
                if $(variable_type_names[i]) <: Val
                    push!(model_parameters, $(QuoteNode(variables[i])))
                  ProbabilityModels.DistributionParameters.load_parameter!(
                      first_pass, second_pass,
                      $(QuoteNode(variables[i])), ProbabilityModels.extract_typeval($(variable_type_names[i])),
                      $return_partials, ProbabilityModels, Symbol("##stack_pointer##"), true, false
                  )
                elseif $(variable_type_names[i]) <: ProbabilityModels.DistributionParameters.MissingDataArray
                    push!(model_parameters, $(QuoteNode(variables[i])))
                    push!(first_pass, $load_incomplete_data)
                  ProbabilityModels.DistributionParameters.load_parameter!(
                      first_pass, second_pass, $(QuoteNode(variables[i])), ($(variable_type_names[i])),
                      $return_partials, ProbabilityModels, Symbol("##stack_pointer##"), true, false
                  )
                elseif $(variable_type_names[i]) <: Union{Array,SubArray{<:Any,<:Any,<:Array}}
                  push!(first_pass, $load_data_dynamic_ptr_array)
                elseif $(variable_type_names[i]) <: ProbabilityModels.PaddedMatrices.AbstractMutableFixedSizeArray
                  push!(first_pass, $load_data_ptr_array)
                else
                    push!(first_pass, $load_data)
                end
            end)
        end
        if return_partials
            processing = quote
                $(verbose_models() > 0 ? :(println("Before differentiating, the model is: \n", ProbabilityModels.PaddedMatrices.simplify_expr(expr), "\n\n")) : nothing)
                ProbabilityModels.ReverseDiffExpressions.reverse_diff_pass!(first_pass, second_pass, expr, tracked_vars, :ProbabilityModels, $(verbose_models() > 1))
                expr_out = quote
                    target = ProbabilityModels.initialize_target($(Symbol("##element_type##")))
                    # $(Symbol("##stack_pointer##")) = $(Symbol("##initial_stack_pointer##"))
                    $(Symbol("##stack_pointer##")) = ProbabilityModels.StackPointer(pointer($(Symbol("##initial_stack_pointer##")), $(Symbol("##element_type##"))))
                    $(Symbol("##θparameter##")) = ProbabilityModels.vectorizable($(Symbol("##θ_parameter_vector##")))
                    # $(Symbol("##θparameter##")) = ProbabilityModels.vectorizable(ProbabilityModels.SIMDPirates.noalias!(pointer($(Symbol("##θ_parameter_vector##")))))
                    # $(Symbol("##∂θparameter##")) = ProbabilityModels.vectorizable($(Symbol("##∂θ_parameter_vector##")))
                    $(Symbol("##∂θparameter##")) = ProbabilityModels.vectorizable(pointer($(Symbol("##∂θ_parameter_vector##"))))
                end
                append!(expr_out.args, first_pass)
                append!(expr_out.args, second_pass)
                push!(expr_out.args, Expr(:(=), Symbol("##scalar_target##"), :(ProbabilityModels.vsum($(name_dict[:target])))))
                push!(expr_out.args, :((isfinite($(Symbol("##scalar_target##"))) && all(isfinite, $(Symbol("##∂θ_parameter_vector##")))) || ($(Symbol("##scalar_target##")) = typemin($(Symbol("##element_type##"))))))
                push!(expr_out.args, Symbol("##scalar_target##"))
            end
        else
            processing = quote
                ProbabilityModels.constant_drop_pass!(first_pass, expr, tracked_vars, $(verbose_models() > 1))
                pushfirst!(first_pass, Expr(:(=), Symbol("##initial_stack_pointer##"), Expr(:call, :pointer, Symbol("##stack_pointer##"))))
                expr_out = quote
                    target = ProbabilityModels.initialize_target($(Symbol("##element_type##")))
                    $(Symbol("##θparameter##")) = ProbabilityModels.vectorizable($(Symbol("##θ_parameter_vector##")))
                end
                append!(expr_out.args, first_pass)
                push!(expr_out.args, Expr(:(=), Symbol("##scalar_target##"), :(ProbabilityModels.vsum($(name_dict[:target])))))
                push!(expr_out.args, Expr(:call, :(ProbabilityModels.lifetime_end!), Symbol("##initial_stack_pointer##"), Val{8192}()))
                push!(expr_out.args, :(isfinite($(Symbol("##scalar_target##"))) ? $(Symbol("##scalar_target##")) : typemin($(Symbol("##element_type##")))))
            end
        end
        push!(θq_body, quote
            tracked_vars = Set(model_parameters)
              name_dict = Dict(:target => :target)
              $processing
              final_quote = ProbabilityModels.StackPointers.stack_pointer_pass(
                  expr_out, $(QuoteNode(Symbol("##stack_pointer##"))),
                  nothing, $(verbose_models() > 1), :ProbabilityModels
              ) |> ProbabilityModels.PaddedMatrices.simplify_expr
              $(verbose_models() > 0) && println(final_quote)
              quote
                  @inbounds begin
                      $final_quote
                  end
              end
        end)
    end
    # Now, we would like to apply
    PaddedMatrices.simplify_expr.((struct_quote, struct_kwarg_quote, θq_value, θq_valuegradient, constrain_quote, cvq, pn, clq, variables))
    #struct_quote, struct_kwarg_quote, θq_value, θq_valuegradient, constrain_quote, cvq, pn, clq, variables))
end

macro model(model_name, expr)
    # struct_quote, struct_kwarg_quote, θq_value, θq_valuegradient, constrain_q, dim_q, variables = generate_generated_funcs_expressions(model_name, expr)
    struct_quote, struct_kwarg_quote, θq_value, θq_valuegradient, constrain_q, cvq, pn, clq, variables = generate_generated_funcs_expressions(model_name, expr)
    printstring = """
        Defined model: $model_name.
        Unknowns: $(variables[1])$([", " * string(variables[i]) for i ∈ 2:length(variables)]...).
    """
    esc(quote
        # $struct_quote; $dim_q; $struct_kwarg_quote; $θq_value; $θq_valuegradient; $constrain_q;
        $struct_quote; $struct_kwarg_quote; $θq_value; $θq_valuegradient; $constrain_q; $cvq; $pn; $clq;
#        ProbabilityModels.Distributed.myid() == 1 && println($printstring)
        println($printstring)
    end)
end
# macro show_model_expr(model_name, expr)
#     struct_quote, struct_kwarg_quote, θq_value, θq_valuegradient, constrain_q, cvq, pn, clq, variables = generate_generated_funcs_expressions(model_name, expr)
#     @show struct_quote
#     @show struct_kwarg_quote
#     @show θq_value
#     @show θq_valuegradient
#     @show constrain_q
#     @show cvq
#     @show, pn
#     @show clq
#     @show variables
#     nothing
# end

# @specialize


text:
using Test
using SafeTestsets

@safetestset "Sample Entropy test" begin include("test_sample_entropy.jl") end


text:
#=
Derivatives of (r, l, P, T) with respect to the mass coordinate.
=#

include("constants.jl")
include("physics.jl")
include("opacity.jl")

# TODO: implement an opacity that changes with composition
const logkappa_spl = opal_test_spl()

function drdm(m, r, l, P, T, rho, mu)
    #=
    Derivative of r with respect to m, from mass conservation.
    =#
    return 1 / (4 * pi * r .^ 2 .* rho)
end # drdm

function dldm(m, r, l, P, T, rho, X, Y)
    #=
    Derivative of l with respect to m, from energy conservation.
    =#
    return epsilon_pp(rho, T, X) + epsilon_cno(rho, T, X, Y)
end # dldm

function dPdm(m, r, l, P, T)
    #=
    Derivative of P with respect to m, from hydrostatic equilibrium.
    =#
    return -G * m / (4 * pi * r .^ 4)
end # dPdm

function dTdm(m, r, l, P, T, mu, kappa)
    #=
    Derivative of T with respect to m, from... thermodynamics.
    =#
    grad_rad = del_rad(m, l, P, T, kappa)
    grad_ad = del_ad()
    stable =  grad_rad < grad_ad
    if stable
        # @debug("dTdm: stable to convection")
        grad = grad_rad
    else
        # @debug("dTdm: unstable to convection")
        grad = grad_ad
    end # if
    factor = - G .* m .* T ./ (4 * pi .* r .^ 4 .* P)
     return factor .* grad
end # dTdm

function deriv(m, r, l, P, T, X, Y)
    #=
    Subroutine for shootf.
    Calulates the derivatives of (r, l, P, T) with respect to mass.
    See Kippenhahn chapter 10, page 89.
    =#
    # @debug("deriv: m, y = ", m, ", ", [r, l, P, T])
    mu = mu_from_composition(X, Y)
    rho = rho_ideal(P, T, mu)
    kappa = try
        opacity(logkappa_spl, rho, T)
    catch
        @debug("deriv: Opacity interpolation failed!")
        @debug("deriv: y = ", [r, l, P, T])
        @debug("deriv: logT = ", log10(T))
        @debug("deriv: logR = ", log10(rho / (T / 1e6)^3))
        throw(DomainError)
    end
    dy = [drdm(m, r, l, P, T, rho, mu),
          dldm(m, r, l, P, T, rho, X, Y),
          dPdm(m, r, l, P, T),
          dTdm(m, r, l, P, T, mu, kappa)]
    # @debug("deriv: dy = ", dy)
    return dy
end # deriv


text:
#jl # EFTfitter.jl - Empty Template

#!jl # ### Parameters
#jl #~============= Parameters =============================================#
parameters = BAT.NamedTupleDist(
    p1 = -2..2,
)

#!jl # ### Observables
#jl #~============= Observables ============================================#
function observable1(params)
    return params.p1
end

#!jl # ### Measurements
#jl #~============= Measurements ===========================================#
measurements = (
    Meas1 = Measurement(observable1, 0.0, uncertainties = (unc1 = 0.1,), active=true),

    #MeasDist = MeasurementDistribution(obs_array, values_array, uncertainties = (unc1 = unc1_array,), active=false),
)

#!jl # ### Correlations
#jl #~============= Correlations ===========================================#
correlations = (
    unc1 = NoCorrelation(active=true),
)

#corr_matrix = to_correlation_matrix(measurements,
##  (:Meas1, :Meas2, 0.1),
#)


text:
"""
    baryweights(x, y = nothing, csy = nothing) -> w

Compute the barycentric weights for the interpolation points `x`
and poles `[y; csy; conj.(csy)]`.
"""
function baryweights(x, y = nothing)
    n = length(x)

    wx = prodpot(x)
    wy = prodpot(x,y)
    s = LogNumber(true, mean(logabs.(wx)) - mean(logabs.(wy)))
    return float.(s.*wy./wx)
end

"""
    prodpot(x[,y::AbstractVector=x])

Evaluate the product potential `prod(x .- y)` at the scalar `x`,
or evaluate the potential pointwise if `x` is a vector. The result
is returned as a (vector of) `LogNumber`s to avoid over- or
underflow.

Pointwise evaluation is provided as a separate function rather than
through the dot syntax `prodpot.(x,Ref(y))` to allow this function
to be overloaded for point sets where `prodpot(x,y)` can be evaluated
more efficiently (e.g. Chebyshev points).
"""
prodpot(x) = default_prodpot(x)
prodpot(x,y) = default_prodpot(x,y)
prodpot(x,y::Nothing) = true

function default_prodpot(x)
    n = length(x)
    return (i -> prodpot(x[i], @view x[[1:i-1; i+1:n]])).(1:n)
end
default_prodpot(x,y) = default_prodpot.(x,Ref(y))
default_prodpot(x::Number,y) = mapreduce(yi->LogNumber(x - yi), *, y, init=LogNumber(one(x-one(eltype(y)))))

"""
    baryfactor(x̂,x,w) -> î,fac

If `all(x̂.!=x)`, returns `fac = sum(w./(x̂ .- x))` and an undefined `î`.
Otherwise, returns an `î` such that `x̂ == x[î]` and an undefined `fac`.
"""
function baryfactor(x̂,x,w)
    @assert length(x) == length(w)
    n = length(x)
    T = typeof(one(eltype(w))/(x̂-one(eltype(x))))
    fac = zero(T)
    for i = 1:n
        x[i] == x̂ && return i,fac
        fac += w[i]/(x̂-x[i])
    end
    return 0,fac
end

"""
    Barycentric(x, y = nothing, w = baryweights(x,y))

Barycentric basis for rational interpolation with interpolation points `x` and poles `y`.

If `isnothing(y)`, rational interpolation reduced to polynomial interpolation.
"""
struct Barycentric{X,Y,W} <: AbstractBasis
    points::X
    poles::Y
    weights::W
end
Barycentric(x, y = nothing) = Barycentric(x,y, baryweights(x,y))

Base.length(B::Barycentric) = length(B.points)

approxpoints(B::Barycentric) = B.points
approxtransform(B::Barycentric,f) = f

function iterate_basis(B::Barycentric, x̂, (i,î,fac) = (1,baryfactor(x̂,B.points,B.weights)...))
    i > length(B) && return nothing
    x = B.points
    w = B.weights
    î == 0 && return w[i]/((x̂ - x[i]) * fac), (i+1,î,fac)
    T = typeof(iterate_basis(B,x̂, (i,0,fac))[1])
    return (i == î ? one(T) : zero(T)), (i+1,î,fac)
end

evaluate_basis(B::Barycentric,i,x̂) = iterate_basis(B,x̂,(i,baryfactor(x̂,B.points,B.weights)...))[1]

# function evaltransform(B::Barycentric, x̂::Number, c::AbstractVector)
#     # TODO: implement this operation with single pass over data
# end


text:
module Compat

using ..LyceumDevTools: gitauthurl, create_git_cmd
using ..LyceumDevTools: with_tempdir, issubpath

using GitHub: Repo, Authorization, PullRequest, Branch
using GitHub: repo, authenticate, create_pull_request, edit_issue, pull_requests

using Dates
using Markdown: Markdown, Paragraph, Bold, Table, List, Header, HorizontalRule

using Pkg
using Pkg.Types: PackageEntry, PackageSpec, VersionSpec, VersionBound, semver_spec, is_stdlib
using Pkg.Types: Context, EnvCache, write_env, manifest_info
using Pkg.Types: Project, projectfile_path, read_package, read_project, write_project
using Pkg.Types: Manifest, manifestfile_path, read_manifest, write_manifest
using Pkg.Operations: abspath!, with_temp_env

using UUIDs
using UnPack

export RepoSpec, update_tomls!, update_tomls_from_test!, ghactions

# "magic bytes" embedded in pull requests to identify if a preexisiting PR
const COMPAT_UUID = "0062fa4e-0639-437d-8ed2-9da17d9c0af2"
# For TagBot
const PR_LABELS = ["compat", "nochangelog"]

const MASTER_BRANCH = "master"
const COMPAT_BRANCH = "compat"
const KEEP_OLD_COMPAT = true
const DROP_PATCH = true
const UPDATE_JULIA_COMPAT = false
const UPDATE_MANIFEST = true

include("types.jl")
include("ci.jl")
include("util.jl")
include("update.jl")

function update_tomls!(
    env_dir::AbstractString;
    keep_old_compat::Bool = KEEP_OLD_COMPAT,
    drop_patch::Bool = DROP_PATCH,
    update_julia_compat::Bool = UPDATE_JULIA_COMPAT,
    update_manifest::Bool = UPDATE_MANIFEST,
)
    # only update the registry once
    @info "Updating registries"
    Pkg.Registry.update()

    @info "Resolving project"
    ctx = Context(env = EnvCache(realpath(projectfile_path(env_dir))))
    Pkg.instantiate(ctx)
    Pkg.resolve(ctx)

    vers = get_updated_versions(ctx)

    result = CompatResult(project_file = ctx.env.project_file)
    _update_compat!(ctx, result, vers, keep_old_compat, drop_patch, update_julia_compat)
    #&& write_project(ctx.env.project, ctx.env.project_file)
    update_manifest && _update_manifest!(ctx, result)

    # Sanity check: the versions one would get after a Pkg.update() match those used to generate
    # the compat entries.
    with_sandbox(ctx) do tmp
        Pkg.update(update_registry = false)
        _check_versions(Context().env.manifest, vers)
    end

    results = [result]
    msg = format_message(env_dir, results)
    println()
    display(msg)
end

function update_tomls_from_test!(
    pkg_dir::AbstractString;
    keep_old_compat::Bool = KEEP_OLD_COMPAT,
    drop_patch::Bool = DROP_PATCH,
    update_julia_compat::Bool = UPDATE_JULIA_COMPAT,
    update_manifest::Bool = UPDATE_MANIFEST,
    remove_redundant::Bool = true,
)
    # only update the registry once
    @info "Updating registries"
    Pkg.Registry.update()

    pkg_ctx = Context(env = EnvCache(realpath(projectfile_path(pkg_dir))))
    pkg_info = Pkg.project(pkg_ctx)
    pkg_info.ispackage || error("Not a package: $pkg_dir")

    @info "Resolving test project"
    test_ctx = Context(env = EnvCache(realpath(projectfile_path(joinpath(pkg_dir, "test")))))
    Pkg.instantiate(test_ctx)
    Pkg.resolve(test_ctx)
    #was_deved = test_ctx.env.manifest TODO

    vers = get_updated_versions_from_test(pkg_ctx, test_ctx)

    pkg_result = CompatResult(project_file = pkg_ctx.env.project_file)
    _update_compat!(pkg_ctx, pkg_result, vers, keep_old_compat, drop_patch, update_julia_compat)
    write_project(pkg_ctx.env.project, pkg_ctx.env.project_file)

    test_result = CompatResult(project_file = test_ctx.env.project_file)
    remove = remove_redundant ? [pkg_info.uuid, values(pkg_ctx.env.project.deps)...] : UUID[]
    _update_compat!(test_ctx, test_result, vers, keep_old_compat, drop_patch, update_julia_compat, remove = remove)
    write_project(test_ctx.env.project, test_ctx.env.project_file)
    # TODO resolve
    if update_manifest
        Pkg.update(test_ctx, update_registry = false)
        test_result.manifest_updated = test_ctx.env.manifest != test_ctx.env.original_manifest
    end

    # Sanity check: the versions one would get after a Pkg.update() match those used to generate
    # the compat entries.
    with_sandbox(test_ctx) do tmp
        Pkg.develop(PackageSpec(uuid=pkg_info.uuid, path=pkg_ctx.env.manifest_file))
        Pkg.update(update_registry = false)
        _check_versions(Context().env.manifest, vers)
    end

    results = [pkg_result, test_result]
    msg = format_message(pkg_dir, results)
    println()
    display(msg)
    return results, msg
end

function _check_versions(man::Manifest, vers::Dict{UUID,VersionNumber})
    for (uuid, entry) in man
        if haskey(vers, uuid) && entry.version != vers[uuid]
            vpre, vpost = vers[uuid], entry.version
            error("Version mismatch after compat update for $(entry.name): $vpre != $vpost")
        end
    end
end

function _update_compat!(
    ctx::Context,
    result::CompatResult,
    vers::Dict{UUID,VersionNumber},
    keep_old_compat::Bool,
    drop_patch::Bool,
    update_julia_compat::Bool;
    remove::Vector{UUID} = UUID[],
)
    deps = ctx.env.project.deps
    old_compat = ctx.env.original_project.compat
    new_compat = ctx.env.project.compat
    for (name, uuid) in deps
        (is_stdlib(uuid) || isjll(name) || !haskey(vers, uuid)) && continue

        ver = vers[uuid]

        if haskey(old_compat, name)
            old_entry = old_compat[name]
            if uuid in remove
                delete!(new_compat, name)
                push!(result.updated, (name, old_entry, nothing))
                continue
            end

            if keep_old_compat
                new_entry = format_compat(old_entry, ver, drop_patch)
            else
                new_entry = format_compat(ver, drop_patch)
            end

            old_spec = semver_spec(old_entry)
            new_spec = semver_spec(new_entry)
            if lowerbound(new_spec) < lowerbound(old_spec) || new_spec == old_spec
                # Don't update the compat entry if it would result in a downgrade
                # and if the specs are equal, keep the old one (however it was formatted)
                # TODO format entry to canonical/compressed form.
                push!(result.unchanged, (name, old_entry))
            else
                new_compat[name] = new_entry
                push!(result.updated, (name, old_entry, new_entry))
            end
        else
            if uuid in remove
                push!(result.unchanged, (name, nothing))
            else
                new_entry = format_compat(ver, drop_patch)
                new_compat[name] = new_entry
                push!(result.new, (name, new_entry))
            end
        end
    end

    if update_julia_compat
        new_compat["julia"] = format_compat(VERSION, drop_patch)
        if new_compat["julia"] != get(old_compat, "julia", nothing)
            result.julia_compat = new_compat["julia"]
        end
    end

    project_updated(result)

    return result
end

end # module


text:
using CUDA
using MOCNeutronTransport
using BenchmarkTools
using Test

# Number of lines to use in vectors
N = 2^20
println("Using a LineSegment_2D vector of length $N")

# Check num threads and give warning
nthreads = Threads.nthreads()
if nthreads === 1
    @warn "Only using single-thread for cpu. Try restarting julia with 'julia --threads n'"
else
    println("Using $nthreads threads for CPU multi-threading")
end

# Single threads CPU intersection
function sequential_intersection!(out, l, lines)
    for i in eachindex(lines)
        @inbounds out[i] = l ∩ lines[i]
    end
    return nothing
end

# Multithreaded CPU intersection
function parallel_intersection!(out, l, lines)
    Threads.@threads for i in eachindex(lines)
        @inbounds out[i] = l ∩ lines[i]
    end
    return nothing
end

# Single threaded GPU intersection
function gpu_1_thread_intersection!(out, l, lines)
    for i = 1:length(lines)
        @inbounds out[i] = l ∩ lines[i]
    end
    return nothing
end

function bench_gpu_1_thread_intersection!(out, l, lines)
    CUDA.@sync begin
        @cuda gpu_1_thread_intersection!(out, l, lines)
    end
end

# Single block GPU intersection
function gpu_1_block_intersection!(out, l, lines)
    index = threadIdx().x    # this example only requires linear indexing, so just use `x`
    stride = blockDim().x
    for i = index:stride:length(lines)
        @inbounds out[i] = l ∩ lines[i]
    end
    return nothing
end

function bench_gpu_1_block_intersection!(out, l, lines)
    CUDA.@sync begin
        @cuda threads=512 gpu_1_block_intersection!(out, l, lines)
    end
end

# Multiple blocks GPU intersection
function gpu_multiblock_intersection!(out, l, lines)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x
    for i = index:stride:length(lines)
        @inbounds out[i] = l ∩ lines[i]
    end
    return nothing
end

function bench_gpu_multiblock_intersection!(out, l, lines)
    numblocks = ceil(Int64, length(lines)/512)
    CUDA.@sync begin
        @cuda threads=512 blocks=numblocks gpu_multiblock_intersection!(out, l, lines)
    end
end

# Use the occupancy API to determine threads and blocks to saturate the GPU
function bench_gpu_multiblock_autooccupancy!(out, l, lines)
    kernel = @cuda launch=false gpu_multiblock_intersection!(out, l, lines)
    config = launch_configuration(kernel.fun)
    threads = min(length(lines), config.threads)
    blocks = cld(length(lines), threads)

    CUDA.@sync begin
        kernel(out, l, lines; threads, blocks)
    end
end

cpu_time = 0.0
for T = [Float64, Float32]
    println("Using LineSegment_2D of type $T")
    l = LineSegment_2D(Point_2D{T}(0, 0.2), Point_2D{T}(1, 0.8))
    lines = rand(LineSegment_2D{T}, N)
    out = fill((false, Point_2D{T}(0, 0)), N)
    ref = l .∩ lines
    time = @belapsed sequential_intersection!($out, $l, $lines)
    μs = 1e6*time
    if T == Float64
        global cpu_time = μs
    end
    speedup = cpu_time/μs
    println("    CPU: single-thread = $μs μs")
    println("         Speed up compared to single-thread CPU & Float64 = $speedup")
    @test out == ref

    fill!(out, (false, Point_2D{T}(0, 0)))
    time = @belapsed parallel_intersection!($out, $l, $lines)
    μs = 1e6*time
    speedup = cpu_time/μs
    println("    CPU: $nthreads threads = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 = $speedup")
    @test out == ref

    lines_d = CuArray(lines)
    out_d = CUDA.fill((false, Point_2D{T}(0, 0)), N)

    time = @belapsed bench_gpu_1_thread_intersection!($out_d, $l, $lines_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    println("    GPU: single-thread/block, 1 blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 = $speedup")
    out_arr = Array(out_d)
    @test getindex.(out_arr, 1) == getindex.(ref, 1)
    @test getindex.(out_arr, 2) ≈ getindex.(ref, 2)

    fill!(out_d, (false, Point_2D{T}(0, 0)))
    time = @belapsed bench_gpu_1_block_intersection!($out_d, $l, $lines_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    println("    GPU: 512 threads/block, 1 blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 = $speedup")
    out_arr = Array(out_d)
    @test getindex.(out_arr, 1) == getindex.(ref, 1)
    @test getindex.(out_arr, 2) ≈ getindex.(ref, 2)

    fill!(out_d, (false, Point_2D{T}(0, 0)))
    time = @belapsed bench_gpu_multiblock_intersection!($out_d, $l, $lines_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    numblocks = ceil(Int64, N/512)
    println("    GPU: 512 threads/block, $numblocks blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 = $speedup")
    out_arr = Array(out_d)
    @test getindex.(out_arr, 1) == getindex.(ref, 1)
    @test getindex.(out_arr, 2) ≈ getindex.(ref, 2)

    fill!(out_d, (false, Point_2D{T}(0, 0)))
    time = @belapsed bench_gpu_multiblock_autooccupancy!($out_d, $l, $lines_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    kernel = @cuda launch=false gpu_multiblock_intersection!(out_d, l, lines_d)
    config = launch_configuration(kernel.fun)
    threads = min(N, config.threads)
    blocks = cld(N, threads)
    println("    GPU: $threads threads/block, $blocks blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 Points = $speedup")
    out_arr = Array(out_d)
    @test getindex.(out_arr, 1) == getindex.(ref, 1)
    @test getindex.(out_arr, 2) ≈ getindex.(ref, 2)
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
module ExSup

    include("csv_processing.jl")
    include("experiment_support.jl") ;

end # module


text: reverse(dict::Dict) = Dict(value => key for (key, value) in dict)


text:
module STMO

using Plots, LinearAlgebra

include("helpers.jl")
export myred, mygreen, myyellow, myblue, myblack, myorange, mycolors
export trace, notrack, PathTrack, nsteps
export path, path!, plotobj, plotobj!, colorscatter
export dist
export EdgeList, WeightedEdgeList, Vertices, AdjList, edges2adjlist, adjlist2edges
export nvertices, isconnected, vertices

include("quadratic.jl")
export fquad, quadratic
export Quadratic, SignalRecovery
#export solve_quadratic, quadratic_ls, gradient_descent, plot_quadratic
#export generate_noisy_measurements, make_bookkeeping, make_connection_matrix, signalfun

include("unconstrained.jl")
export Unconstrained
#export gradient_descent, coordinate_descent, newtons_method, backtracking_line_search

include("optimaltransport.jl")
export OptimalTransport

include("mst.jl")
export MST

include("testfuns.jl")
export TestFuns

include("tickettoride.jl")
export TicketToRide

end # module

max_stars_repo_name: mferrodo/STMO

max_issues_repo_name: mferrodo/STMO

max_forks_repo_name: mferrodo/STMO


text:
push!(LOAD_PATH,"../src/")
using CorrDCA
using Documenter
#DocMeta.setdocmeta!(CorrDCA, :DocTestSetup, :(using CorrDCA); recursive=true)
makedocs(;
    modules=[CorrDCA],
    authors="Andrea Pagnani",
    repo="https://github.com/pagnani/CorrDCA/blob/{commit}{path}#L{line}",
    sitename="CorrDCA",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pagnani.github.io/CorrDCA",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    branch = "gh-pages",
    devbranch = "main",
	repo = "github.com/pagnani/CorrDCA.git",
    versions = ["stable" => "v^"]
)

max_stars_repo_name: pagnani/CorrDCA

max_issues_repo_name: pagnani/CorrDCA

max_forks_repo_name: pagnani/CorrDCA


text:
using Printf

function rectilinear(fname::AbstractString,
                     xAxis::Array, yAxis::Array, zAxis::Array,
                     data::Array)
    """
    Save data on rectilinear grid in VTK file
    @param fname file name
    @param xAxis x axis (nodes)
    @param yAxis y axis (nodes)
    @param zAxis z axis (nodes)
    @param data field data
    """
    f = open(fname, "w")
    write(f, "# vtk DataFile Version 2.0\n")
    write(f, "upwind.jl\n")
    write(f, "ASCII\n")
    write(f, "DATASET RECTILINEAR_GRID\n")
    write(f, "DIMENSIONS")

    numCells = [length(xAxis) - 1, length(yAxis) - 1, length(zAxis) - 1]
    ntot = length(data)

    # Number of nodes in each direction
    write(f, @sprintf(" %d", numCells[1] + 1))
    write(f, @sprintf(" %d", numCells[2] + 1))
    write(f, @sprintf(" %d\n", numCells[3] + 1))

    write(f,  "X_COORDINATES ")
    write(f,  @sprintf("%d double\n", numCells[1] + 1))
    for i in 1:numCells[1] + 1
        write(f,  @sprintf("%f ", xAxis[i]))
    end

    write(f,  "\nY_COORDINATES ")
    write(f, @sprintf("%d double\n", numCells[2] + 1))
    for i in 1:numCells[2] + 1
        write(f, @sprintf("%f ", yAxis[i]))
    end

    write(f,  "\nZ_COORDINATES ")
    write(f, @sprintf("%d double\n", numCells[3] + 1))
    for i in 1:numCells[3] + 1
        write(f, @sprintf("%f ", zAxis[i]))
    end

    write(f, @sprintf("\nCELL_DATA %d\n", ntot))
    write(f, "SCALARS f double 1\n")
    write(f, "LOOKUP_TABLE default\n")
    for i in eachindex(data)
        write(f, @sprintf("%f ", data[i]))
        if i % 10 == 0
            write(f, "\n")
        end
    end
    write(f, "\n")

    close(f)

end

max_stars_repo_name: pletzer/fidibench

max_issues_repo_name: pletzer/fidibench

max_forks_repo_name: pletzer/fidibench


text:
const minimal_default = Attributes(
    font = "Dejavu Sans",
    backgroundcolor = RGBAf0(1,1,1,1),
    color = :black,
    colormap = :viridis,
    resolution = reasonable_resolution(),
    visible = true,
    clear = true,
    show_axis = true,
    show_legend = false,
    scale_plot = true,
    center = true,
    axis = Attributes(),
    legend = Attributes(),
    axis_type = automatic,
    camera = automatic,
    limits = automatic,
    padding = Vec3f0(0.1),
    raw = false
)

const _current_default_theme = Attributes(; minimal_default...) # make a copy. TODO overload copy?

function current_default_theme(; kw_args...)
    new_theme, rest = merge_attributes!(Attributes(kw_args), _current_default_theme)
    merge!(new_theme, rest)
end

function set_theme!(new_theme::Attributes)
    empty!(_current_default_theme)
    new_theme, rest = merge_attributes!(new_theme, minimal_default)
    merge!(_current_default_theme, new_theme, rest)
    return
end
function set_theme!(;kw_args...)
    set_theme!(Attributes(; kw_args...))
end


text:
using Jeff, Distributions

SLD = Jeff.Parameter(Uniform(6.3, 6.4))
ISLD = Jeff.Parameter(0.000)
THICK = Jeff.Parameter(50.0)
ROUGH = Jeff.Parameter(LogNormal(5, 1))
LAYER = Jeff.Layer(THICK, SLD, ISLD, ROUGH)

@testset "parameter_struct" begin
    @test isapprox(minimum(SLD.value), 6.3)
    @test isapprox(maximum(SLD.value), 6.4)
end

@testset "parameter_struct_no_vary" begin
    @test isapprox(THICK.value, 50.0)
end

function test_layer(layer::Jeff.Layer)
    @test isapprox(layer.thick.value, 50.0)
    @test isapprox(minimum(layer.sld.value), 6.3)
    @test isapprox(maximum(layer.sld.value), 6.4)
    @test isapprox(layer.isld.value, 0.000)
    @test isapprox(meanlogx(layer.rough.value), 5.)
    @test isapprox(stdlogx(layer.rough.value), 1.)
end

@testset "layer_struct" begin
    test_layer(LAYER)
end

@testset "model_struct" begin
    scale = Jeff.Parameter(1,)
    bkg = Jeff.Parameter(Uniform(1e-8, 1e-6))
    model = Jeff.Model(scale, bkg, [LAYER, LAYER])
    @test isapprox(model.scale.value, 1)
    @test isapprox(minimum(model.bkg.value), 1e-8)
    @test isapprox(maximum(model.bkg.value), 1e-6)
    test_layer(model.layers[1])
    test_layer(model.layers[2])
end


text:
import Pkg

original_directory = pwd()
project_root = joinpath(splitpath(@__DIR__)...)
cd(project_root)

include(joinpath(project_root, "ci", "RegistryTestingTools", "types.jl",))

if length(ARGS) == 0
    push!(ARGS, "[,)",)
end

this_registry_toml_path = joinpath(project_root, "Registry.toml")
registry_configuration = Pkg.TOML.parsefile(this_registry_toml_path)
name_to_path = Dict{String, String}()
all_packages = String[]
for pair in registry_configuration["packages"]
    name = pair[2]["name"]
    path = pair[2]["path"]
    push!(all_packages, name,)
    name_to_path[name] = path
end

unique!(all_packages)
sort!(all_packages)
n = length(all_packages)
@debug("all_packages ($(n)):")
for i = 1:n
    @debug("$(i). $(all_packages[i])")
end
if _is_interval(ARGS[1])
    this_job_interval = _construct_interval(
        convert(String, strip(ARGS[1]))
        )
    _this_job_interval_contains_x(x) = _interval_contains_x(
        this_job_interval,
        x,
        )
    packages_to_compress_in_this_job_interval = all_packages[
        _this_job_interval_contains_x.(
            all_packages
            )
        ]
else
    packages_to_compress_in_this_job_interval = String[
        strip(s) for s in ARGS
        ]
end
unique!(packages_to_compress_in_this_job_interval)
sort!(packages_to_compress_in_this_job_interval)
n = length(packages_to_compress_in_this_job_interval)
@debug("packages_to_compress_in_this_job_interval ($(n)):")
for i = 1:n
    @debug("$(i). $(packages_to_compress_in_this_job_interval[i])")
end

n = length(packages_to_compress_in_this_job_interval)
for i = 1:n
    name = packages_to_compress_in_this_job_interval[i]
    @debug("Compressing \"$(name)\" (package $(i) of $(n))")
    path = name_to_path[name]
    for j = 1:3
        for filename in ["Versions.toml", "Compat.toml", "Deps.toml",]
            full_file_path = joinpath(path, filename,)
            uncompressed = Pkg.Compress.load(full_file_path,)
            Pkg.Compress.save(full_file_path,uncompressed,)
        end
    end
    @info("Compressed package \"$(name)\"")
end

temp_registry_toml = read(this_registry_toml_path, String)
temp_registry_toml_parsed = Pkg.TOML.parse(temp_registry_toml)
if !haskey(temp_registry_toml_parsed, "packages")
    temp_output = string(temp_registry_toml,
                         "\n\n\n\n",
                         "[packages]",
                         "\n\n\n\n")
    rm(this_registry_toml_path; force = true, recursive = true)
    write(this_registry_toml_path, temp_output)
end

cd(original_directory)


text:
using PowerModels
using Graphs

type Graph
    graph               ::SimpleGraph
    vertex_map          ::Dict{Int,Any}
    edge_map            ::Dict{Int,Any}
    bus_to_vertex_map   ::Dict{Int,Int}
    data                ::Dict{AbstractString,Any}
    ref                 ::Dict{Symbol,Any}
    function Graph()
        gm = new()
        gm.graph = simple_graph(0, is_directed=false)
        gm.vertex_map = Dict{Int,Any}()
        gm.edge_map = Dict{Int,Any}()
        gm.bus_to_vertex_map = Dict{Int,Int}()
        gm.data = Dict{AbstractString,Any}()
        return gm
    end
end

type VertexNeighbors
    id          ::Int
    edge_ids    ::Vector{Int}
    function VertexNeighbors()
        vn = new()
        vn.id = 0
        vn.edge_ids = Int[]
        return vn
    end
end

type Attack
    bus_ids             ::Vector{Int}
    branch_ids          ::Vector{Int}
    total_load          ::Float64
    function Attack()
        a = new()
        a.bus_ids = Int[]
        a.branch_ids = Int[]
        a.total_load = 0.0
        return a
    end
end

function populate_gm(gm::Graph, buses::Dict{Int,Any}, branches::Dict{Int,Any})
    i = 1
    for bus in buses
        add_vertex!(gm.graph)
        key = bus[1]
        value = bus[2]
        gm.vertex_map[i] = value
        gm.bus_to_vertex_map[key] = i
        i += 1
    end

    i = 1
    for (key,value) in branches
        add_edge!(gm.graph, gm.bus_to_vertex_map[value["f_bus"]], gm.bus_to_vertex_map[value["t_bus"]])
        gm.edge_map[i] = value
        i += 1
    end

    for edge in edges(gm.graph)
        @assert gm.vertex_map[edge.source]["index"] == gm.edge_map[edge.index]["f_bus"]
        @assert gm.vertex_map[edge.target]["index"] == gm.edge_map[edge.index]["t_bus"]
    end
    return
end

function get_vertices(gm::Graph; degree = 1)
    vertices_with_degree = Dict(v => out_degree(v, gm.graph) for v in vertices(gm.graph))
    vertices_with_degree = filter((v,deg) -> deg == degree, vertices_with_degree)
    neighbors = Dict(v => collect(out_edges(v, gm.graph)) for v in collect(keys(vertices_with_degree)))

    vertex_neighbors = VertexNeighbors[]
    for (v,n) in neighbors
        if gm.vertex_map[v]["bus_type"] == 2 || gm.vertex_map[v]["bus_type"] == 3
            continue
        end
        vn = VertexNeighbors()
        vn.id = gm.vertex_map[v]["index"]
        for edge in n
            push!(vn.edge_ids, gm.edge_map[edge.index]["index"])
        end
        push!(vertex_neighbors, vn)
    end

    return vertex_neighbors
end

function create_n_1_graph(gm::Graph, buses, branches, line)
    i = 1
    for bus in buses
        add_vertex!(gm.graph)
        key = bus[1]
        value = bus[2]
        gm.vertex_map[i] = value
        gm.bus_to_vertex_map[key] = i
        i += 1
    end

    i = 1
    for branch in branches
        if branch == line
            continue
        end
        key = branch[1]
        value = branch[2]
        add_edge!(gm.graph, gm.bus_to_vertex_map[value["f_bus"]], gm.bus_to_vertex_map[value["t_bus"]])
        gm.edge_map[i] = value
        i += 1
    end

    for edge in edges(gm.graph)
        @assert gm.vertex_map[edge.source]["index"] == gm.edge_map[edge.index]["f_bus"]
        @assert gm.vertex_map[edge.target]["index"] == gm.edge_map[edge.index]["t_bus"]
    end

    return
end

function check_graph_for_isolation(gm::Graph, buses, branches)
    components = connected_components(gm.graph)
    isolated_load_vertices = Attack[]
    if length(components) == 1
        return Attack[]
    end
    ref = gm.ref
    for component in components
        is_gen = false
        for v in component
            i = gm.vertex_map[v]["index"]
            if buses[i]["bus_type"] == 2 || buses[i]["bus_type"] == 3
                total_pmin = sum([abs(ref[:gen][g]["pmin"]) for g in ref[:bus_gens][i]])
                total_pmax = sum([abs(ref[:gen][g]["pmax"]) for g in ref[:bus_gens][i]])
                is_condensor = (total_pmin == 0 && total_pmax == 0) ? true : false
                if is_condensor
                    continue
                else
                    is_gen = true
                    break
                end
            end
        end
        if !(is_gen)
            load = Attack()
            for v in component
                push!(load.bus_ids, gm.vertex_map[v]["index"])
                load.total_load += abs(buses[gm.vertex_map[v]["index"]]["pd"])
            end
            push!(isolated_load_vertices, load)

        end
    end

    return isolated_load_vertices
end

function remove_redundant_attacks(attacks::Vector{Attack})
    redundant_attack_indexes = Set()
    num_attacks = length(attacks)
    for i in 1:num_attacks
        for j in 1:num_attacks
            if j == i
                continue
            end
            num_elements_i = length(attacks[i].bus_ids)
            num_elements_j = length(attacks[j].bus_ids)
            attack_set_i = Set(attacks[i].bus_ids)
            attack_set_j = Set(attacks[j].bus_ids)
            intersection = intersect(attack_set_i, attack_set_j)
            if attack_set_i == intersection
                push!(redundant_attack_indexes, i)
            end
            if attack_set_j == intersection
                push!(redundant_attack_indexes, j)
            end
        end
    end
    deleteat!(attacks, sort(collect(redundant_attack_indexes)))
    return
end

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
# run_random_set
fname = "../poets/basic_set_generator.jl"
include(fname)

var = ARGS[1]

start_random = false # true
global HOLD=0
fname = "../poets/emt_grn_model_objective_function.jl"
include(fname)

# Set parameters and pointers

#initial_state = Array{Float64,1}(parameter_array)
full_pc_array = readdlm("../poets/results_poets/cross_validation/fold_pop_obj_$var/PC.dat")

# Execute the beast
println(now())
@show size(full_pc_array)[2]
#stop()
EC = generate_set(objective_function,full_pc_array)
println(now())
# Save output

outdir= "../poets/results_poets/cross_validation"
writedlm("$outdir/EC$var.dat", EC)
#writedlm("$outdir/PC.dat", PC)
#writedlm("$outdir/RA.dat", RA)


text:
module SimpleValueGraphs

using Graphs: AbstractGraph, AbstractEdgeIter, AbstractEdge, num_self_loops

using Graphs.SimpleGraphs: AbstractSimpleGraph, AbstractSimpleEdge,
    SimpleGraph, SimpleDiGraph, SimpleEdge, IsDirected

using Base: OneTo

import Graphs

import Graphs:
    nv, ne, is_directed,
    eltype, vertices, add_vertex!, rem_vertex!, has_vertex,
    edgetype, edges, src, dst, reverse,
    weights,
    add_edge!, rem_edge!, has_edge,
    inneighbors, outneighbors,
    zero,
    adjacency_matrix,
    squash, reverse

import Base: show, iterate, length,
    getindex, size, ==, hash

import LinearAlgebra
import LinearAlgebra:
    ishermitian, issymmetric,
    adjoint, transpose

import SparseArrays:
    SparseMatrixCSC

export
    ValGraph,
    ValOutDiGraph,
    ValDiGraph,

    AbstractValEdge,
    ValEdge,
    ValDiEdge,

    vertexvals_type,
    edgevals_type,
    graphvals_type,

    hasvertexkey,
    hasedgekey,
    hasgraphkey,

    get_vertexval,
    set_vertexval!,

    get_edgeval,
    get_edgeval_or,
    set_edgeval!,

    get_graphval,
    set_graphval!,

    inedgevals,
    outedgevals,

    AdjacencyMatrix,
    ValMatrix,

    # overridden methods from Graphs.jl
    nv,
    ne,
    is_directed,

    eltype,
    vertices,

    edgetype,
    edges,
    src,
    dst,
    reverse,

    weights,

    add_edge!,
    rem_edge!,
    has_edge,

    add_vertex!,
    rem_vertex!,
    has_vertex,

    inneighbors,
    outneighbors,

    adjacency_matrix,

    squash, reverse,

    zero,

    # overridden methods from Base
    iterate,
    length,
    size,
    getindex,
    show,
    ==,
    hash,

    # overridden methods from LinearAlgebra
    ishermitian,
    adjoint,
    transpose,
    issymmetric,

    # overriden methods from SparseArrays
    SparseMatrixCSC

# ==== Includes ===========================

include("AbstractTuples.jl")
include("utils.jl")

include("abstractvaluegraph.jl")
include("graphwrappers.jl")
include("valueedge.jl")
include("valuegraph.jl")
include("matrices.jl")
include("operators.jl")

include("integrations/SimpleGraphs.jl")

include("Experimental.jl")

# ==== Various ===========================

"""
    swissmetro_graph()

A small example graph for using in documentation.

Swissmetro was a planned (but never realised) Hyperloop style project in Switzerland.
All data was taken from Wikipedia: https://en.wikipedia.org/wiki/Swissmetro
"""
function swissmetro_graph()

    cities = [
        (name = "Basel", population = 117_595),
        (name = "Bern", population = 133_791),
        (name = "Genève", population = 201_818),
        (name = "Lausanne", population = 139_111),
        (name = "St. Gallen", population = 75_833),
        (name = "Zürich", population = 415_215),
    ]

    g = ValGraph{Int8}(
        6;
        vertexval_types=(name = String, population = Int32),
        vertexval_init=v -> cities[v],
        edgeval_types=(distance=Float64,)
    )

    add_edge!(g, 1, 6, (distance=89.0,))
    add_edge!(g, 2, 4, (distance=81.0,))
    add_edge!(g, 2, 6, (distance=104.0,))
    add_edge!(g, 3, 4, (distance=68.0,))
    add_edge!(g, 5, 6, (distance=69.0,))

    return g
end

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.3.0
#   kernelspec:
#     display_name: Julia 1.3.0
#     language: julia
#     name: julia-1.3
# ---

# + [markdown] slideshow={"slide_type": "slide"}
# # Julia and DifferentialEquations.jl
#

# +
using Plots

x = LinRange(0, 2π, 10)
y = sin.(x) + cos.(x).^2

plot(x, y)

# +
using PyCall
@pyimport scipy.interpolate as interpolate

f = interpolate.interp1d(x, y, kind="cubic")
new_x = LinRange(first(x), last(x), 1000)
plot( new_x, f(new_x))

# +
using RCall

R"""
library(missMDA)
data(orange)
nb = estim_ncpPCA(orange,ncp.max=5)
resMI = MIPCA(orange,ncp=2)
plot(resMI)
"""

max_stars_repo_name: pnavaro/DiffEq

max_issues_repo_name: pnavaro/DiffEq

max_forks_repo_name: pnavaro/DiffEq


text:
module NVML

using ..APIUtils

using ..CUDA

using CEnum: @cenum

import Libdl

function libnvml()
    @memoize begin
        if Sys.iswindows()
            # the NVSMI dir isn't added to PATH by the installer
            nvsmi = joinpath(ENV["ProgramFiles"], "NVIDIA Corporation", "NVSMI")
            if isdir(nvsmi)
                joinpath(nvsmi, "nvml.dll")
            else
                # let's just hope for the best
                "nvml"
            end
        else
            "libnvidia-ml.so.1"
        end
    end::String
end

function has_nvml()
    @memoize begin
        if Libdl.dlopen(libnvml(); throw_error=false) === nothing
            return false
        end

        # JuliaGPU/CUDA.jl#860: initialization can fail on Windows
        try
            initialize_api()
        catch err
            @error "Cannot use NVML, as it failed to initialize" exception=(err, catch_backtrace())
            return false
        end

        return true
    end::Bool
end

# core library
include("libnvml_common.jl")
include("error.jl")
include("libnvml.jl")
include("libnvml_deprecated.jl")

# wrappers
include("system.jl")
include("device.jl")

end


text:
using ReplMaker
using REPL
using Serialization
using Sockets

using OpenSSH_jll

function verify_header(io, ser_version=Serialization.ser_version)
    magic = String(read(io, length(protocol_magic)))
    if magic != protocol_magic
        error("RemoteREPL protocol magic number mismatch: $(repr(magic)) != $(repr(protocol_magic))")
    end
    version = read(io, typeof(protocol_version))
    if version != protocol_version
        error("RemoteREPL protocol version number mismatch: $version != $protocol_version")
    end
    # Version 1: We rely on the standard Serialization library for simplicity;
    # that's backward but not forward compatible depending on
    # Serialization.ser_version, so we check for an exact match
    remote_ser_version = read(io, UInt32)
    if remote_ser_version != ser_version
        error("RemoteREPL Julia Serialization format version mismatch: $remote_ser_version != $ser_version")
    end
    return true
end

# Find a free port on `network_interface`
function find_free_port(network_interface)
    # listen on port 0 => kernel chooses a free port. See, for example,
    # https://stackoverflow.com/questions/44875422/how-to-pick-a-free-port-for-a-subprocess
    server = listen(network_interface, 0)
    _, free_port = getsockname(server)
    close(server)
    # The kernel can reuse free_port here after $some_time_delay, but apprently
    # this is large enough for Selenium to have used this technique for ten
    # years...
    return free_port
end

function comm_pipeline(cmd::Cmd)
    errbuf = IOBuffer()
    proc = run(pipeline(cmd, stdout=errbuf, stderr=errbuf),
               wait=false)
    atexit() do
        kill(proc)
    end
    @async begin
        # Attempt to log any connection errors to the user
        wait(proc)
        errors = String(take!(errbuf))
        if !isempty(errors) || !success(proc)
            @warn "Tunnel output" errors=Text(errors)
        end
    end
    proc
end

function ssh_tunnel(host, port, tunnel_interface, tunnel_port; ssh_opts=``)
    OpenSSH_jll.ssh() do ssh_exe
        # Tunnel binds locally to $tunnel_interface:$tunnel_port
        # The other end jumps through $host using the provided identity,
        # and forwards the data to $port on *itself* (this is the localhost:$port
        # part - "localhost" being resolved relative to $host)
        ssh_cmd = `$ssh_exe $ssh_opts -o ExitOnForwardFailure=yes -o ServerAliveInterval=60
                            -N -L $tunnel_interface:$tunnel_port:localhost:$port $host`
        @debug "Connecting SSH tunnel to remote address $host via ssh tunnel to $port" ssh_cmd
        comm_pipeline(ssh_cmd)
    end
end

function aws_tunnel(instance_id, port, tunnel_port; region=nothing)
    region = region === nothing ? `` : `--region $region`
    aws_cmd = `aws ssm start-session $region --target $instance_id --document-name AWS-StartPortForwardingSession --parameters "{\"portNumber\":[\"$port\"],\"localPortNumber\":[\"$tunnel_port\"]}"`
    @debug "Connecting AWS session manager tunnel to EC2 instance $instance_id via port $port" aws_cmd
    comm_pipeline(aws_cmd)
end

function k8s_tunnel(host, port, tunnel_port; namespace=nothing)
    namespace = namespace === nothing ? `` : `-n $namespace`
    k8s_cmd = `kubectl port-forward $namespace $host $tunnel_port:$port`
    @debug "Connecting through kubectl to resource $host via port $port" k8s_cmd
    comm_pipeline(k8s_cmd)
end

function connect_via_tunnel(host, port; retry_timeout,
        tunnel=:ssh, ssh_opts=``, region=nothing, namespace=nothing)
    # We assume the remote server is only listening for local connections.
    tunnel_interface = Sockets.localhost
    tunnel_port = find_free_port(tunnel_interface)
    comm_proc = if tunnel == :ssh
            ssh_tunnel(host, port, tunnel_interface, tunnel_port; ssh_opts=ssh_opts)
        elseif tunnel == :aws
            aws_tunnel(host, port, tunnel_port; region=region)
        elseif tunnel == :k8s
            k8s_tunnel(host, port, tunnel_port; namespace=namespace)
        end

    # Retry loop to wait for the connection.
    for i=1:retry_timeout
        try
            return connect(tunnel_interface, tunnel_port)
        catch exc
            if (exc isa Base.IOError) && process_running(comm_proc) && i < retry_timeout
                sleep(1)
            else
                kill(comm_proc)
                wait(comm_proc)
                @error "Exceeded maximum socket connection attempts"
                rethrow()
            end
        end
    end
end

function valid_input_checker(prompt_state)
    ast = Base.parse_input_line(String(take!(copy(REPL.LineEdit.buffer(prompt_state)))),
                                depwarn=false)
    return !Meta.isexpr(ast, :incomplete)
end

struct RemoteCompletionProvider <: REPL.LineEdit.CompletionProvider
    socket
end

function REPL.complete_line(c::RemoteCompletionProvider, state::REPL.LineEdit.PromptState)
    # See REPL.jl complete_line(c::REPLCompletionProvider, s::PromptState)
    partial = REPL.beforecursor(state.input_buffer)
    full = REPL.LineEdit.input_string(state)
    serialize(c.socket, (:repl_completion, (partial, full)))
    messageid, value = deserialize(c.socket)
    if messageid != :completion_result
        @warn "Completion failure" messageid
        return ([], "", false)
    end
    return value
end

function run_remote_repl_command(socket, out_stream, cmdstr)
    ast = Base.parse_input_line(cmdstr, depwarn=false)
    messageid=nothing
    value=nothing
    try
        # See REPL.jl: display(d::REPLDisplay, mime::MIME"text/plain", x)
        display_props = Dict(
            :displaysize=>displaysize(out_stream),
            :color=>get(out_stream, :color, false),
            :limit=>true,
            :module=>Main,
        )
        serialize(socket, (:display_properties, display_props))
        serialize(socket, (:eval, ast))
        flush(socket)
        response = deserialize(socket)
        messageid, value = response isa Tuple && length(response) == 2 ?
                           response : (nothing,nothing)
    catch exc
        if exc isa Base.IOError
            messageid = :error
            value = "IOError - Remote Julia exited or is inaccessible"
        else
            rethrow()
        end
    end
    if messageid == :eval_result || messageid == :error
        if !isnothing(value) && !REPL.ends_with_semicolon(cmdstr)
            println(out_stream, value)
        end
    else
        @error "Unexpected response from server" messageid
    end
end

function setup_connection(host, port;
                          tunnel = (host!=Sockets.localhost) ? :ssh : :none,
                          ssh_opts=``, region=nothing, namespace=nothing)
    socket = if tunnel == :none
            connect(host, port)
        else
            connect_via_tunnel(host, port; retry_timeout=5,
                tunnel=tunnel, ssh_opts=ssh_opts, region=region,
                namespace=namespace)
        end

    try
        verify_header(socket)
    catch exc
        close(socket)
        rethrow()
    end

    atexit() do
        close_connection(socket)
    end

    return socket
end

function close_connection(socket)
    if isopen(socket)
        serialize(socket, (:exit,nothing))
        close(socket)
    end
end

"""
    connect_repl([host=localhost,] port::Integer=27754;
                 use_ssh_tunnel = (host != localhost),
                 ssh_opts = ``)

Connect client REPL to a remote `host` on `port`. This is then accessible as a
remote sub-repl of the current Julia session.

For security, `connect_repl()` uses an ssh tunnel for remote hosts. This means
that `host` needs to be running an ssh server and you need ssh credentials set
up for use on that host. For secure networks this can be disabled by setting
`tunnel=:none`.

To provide extra options to SSH, you may use the `ssh_opts` keyword, for example an identity file may be set with  `` ssh_opts = `-i /path/to/identity.pem` ``.
Alternatively, you may want to set this up permanently using a `Host` section in your ssh config file.

You can also use the following technologies for tunneling in place of SSH:
1) AWS Session Manager: set `tunnel=:aws`. The optional `region` keyword argument can be used to specify the AWS Region of your server.
2) kubectl: set `tunnel=:k8s`. The optional `namespace` keyword argument can be used to specify the namespace of your Kubernetes resource.

See README.md for more information.
"""
function connect_repl(host=Sockets.localhost, port::Integer=27754;
                      tunnel::Symbol = host!=Sockets.localhost ? :ssh : :none,
                      ssh_opts=``, region=nothing, namespace=nothing)
    socket = setup_connection(host, port, tunnel=tunnel, ssh_opts=ssh_opts,
                 region=region, namespace=namespace)
    out_stream = stdout
    ReplMaker.initrepl(c->run_remote_repl_command(socket, out_stream, c),
                       repl         = Base.active_repl,
                       valid_input_checker = valid_input_checker,
                       prompt_text  = "remote> ",
                       prompt_color = :magenta,
                       start_key    = '>',
                       sticky_mode  = true,
                       mode_name    = "remote_repl",
                       completion_provider = RemoteCompletionProvider(socket)
                       )
                       # startup_text = false)
    nothing
end

connect_repl(port::Integer) = connect_repl(Sockets.localhost, port)

#-------------------------------------------------------------------------------
"""
    remote_eval(cmdstr)
    remote_eval(host, port, cmdstr)

Parse a string `cmdstr`, evaluate it in the remote REPL server's `Main` module,
then close the connection.

For example, to cause the remote Julia instance to exit, you could use

```
using RemoteREPL
RemoteREPL.remote_eval("exit()")
```
"""
function remote_eval(host, port::Integer, cmdstr::AbstractString;
                     tunnel::Symbol = host!=Sockets.localhost ? :ssh : :none)
    socket = setup_connection(host, port, tunnel=tunnel)
    io = IOBuffer()
    run_remote_repl_command(socket, io, cmdstr)
    close_connection(socket)
    String(take!(io))
end

function remote_eval(cmdstr::AbstractString)
    remote_eval(Sockets.localhost, 27754, cmdstr)
end


text:
using Test
using Pluto
using Pluto.Configuration
using Pluto.Configuration: overlayed, notebook_path_suggestion, from_flat_kwargs
using Pluto.WorkspaceManager: _convert_to_flags

@testset "Configurations" begin

cd(Pluto.project_relative_path("test")) do
    @test notebook_path_suggestion() == joinpath(pwd(), "")
end

@testset "overlayed" begin
    opt = Configuration.CompilerOptions()
    @test opt.compile === nothing
    @test overlayed(opt; compile="min").compile == "min"
end

@testset "from_flat_kwargs" begin
    opt = from_flat_kwargs(;compile="min", project="test")
    @test opt.compiler.compile == "min"
    @test opt.compiler.project == "test"

    @test_throws ArgumentError from_flat_kwargs(;fake_project="test")
end

@testset "flag conversion" begin
    @test _convert_to_flags(Configuration.CompilerOptions()) ==
        ["--project=@.", "--startup-file=no", "--history-file=no"]

    @test _convert_to_flags(Configuration.CompilerOptions(compile="min")) ==
    ["--compile=min", "--project=@.", "--startup-file=no", "--history-file=no"]

    @test _convert_to_flags(Configuration.CompilerOptions(compile="min", project="test")) ==
    ["--compile=min", "--project=test", "--startup-file=no", "--history-file=no"]
end

end # testset


text:
#=
Popis
  Funkcia vypocita priemernu hodnotu z najkratsich ciest v grafe.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozene ohodnotenie hrany :: String
  components : parameter pre prepinanie graf/komponent :: Bool

Navratova hodnota
  avg_distance : priemerna hodnota najkratsej cesty pre cely graf :: Float64
  avg_distance : slovnik s mapovanim (komponent_grafu => avg_distance) :: Dict{Array{ExVertex,1},Float64}

Poznamky
  Priemernu hodnotu najkratsej cesty je mozne nastavenim parametra components vypocitat bud pre cely graf
  alebo pre jednotlive komponenty.
=#

function average_shortest_distance(g::GenericGraph; dist_key::AbstractString="weight", components::Bool=false)

  epi = AttributeEdgePropertyInspector{Float64}(dist_key)

  if components

    comp = connected_components(g)
    avg_distance = [v => 0.0 for v = comp]
    paths = Float64[]

    for vector in comp

      for v in vector

        paths = [filter(x -> x != Inf, dijkstra_shortest_paths(g, epi, [v]).dists) for v = vector ]

      end

      avg_distance[vector] = sum([sum(x) for x = paths]) / (length(vector)*(length(vector) - 1))
      empty!(paths)

    end

    return avg_distance

  elseif !(components)

    avg_distance = 0.0
    paths = [filter(x -> x != Inf, dijkstra_shortest_paths(g, epi, [v]).dists) for v = vertices(g) ]
    avg_distance = sum([sum(x) for x = paths]) / (length(vertices(g))*(length(vertices(g)) - 1))

    return avg_distance

  end

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
# ---
# title: 492. Construct the Rectangle
# id: problem492
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/construct-the-rectangle/description/>
# hidden: true
# ---
#
# A web developer needs to know how to design a web page's size. So, given a
# specific rectangular web page's area, your job by now is to design a
# rectangular web page, whose length L and width W satisfy the following
# requirements:
#
#   1. The area of the rectangular web page you designed must equal to the given target area.
#   2. The width `W` should not be larger than the length `L`, which means `L >= W`.
#   3. The difference between length `L` and width `W` should be as small as possible.
#
# Return _an array`[L, W]` where `L` and `W` are the length and width of the web
# page you designed in sequence._
#
#
#
# **Example 1:**
#
#
#
#     Input: area = 4
#     Output: [2,2]
#     Explanation: The target area is 4, and all the possible ways to construct it are [1,4], [2,2], [4,1].
#     But according to requirement 2, [1,4] is illegal; according to requirement 3,  [4,1] is not optimal compared to [2,2]. So the length L is 2, and the width W is 2.
#
#
# **Example 2:**
#
#
#
#     Input: area = 37
#     Output: [37,1]
#
#
# **Example 3:**
#
#
#
#     Input: area = 122122
#     Output: [427,286]
#
#
#
#
# **Constraints:**
#
#   * `1 <= area <= 107`
#
#
## @lc code=start
using LeetCode

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


text:
module ChebyshevExpansions
    # Use README as the docstring of the module:
    @doc read(joinpath(dirname(@__DIR__), "README.md"), String) ChebyshevExpansions
    using LinearAlgebra, SparseArrays, ProgressMeter, FFTW, ArnoldiMethod
    export randomkets, momentaKPM, densityKPM, dosKPM, averageKPM
    #export momentaKPM, densityKPM, dosKPM, averageKPM
    include("KPM.jl")
end

max_stars_repo_name: fernandopenaranda/ChebyshevExpansions

max_issues_repo_name: fernandopenaranda/ChebyshevExpansions

max_forks_repo_name: fernandopenaranda/ChebyshevExpansions