text:
#############################################################################################
#                             COPY TO BLOCK MATRICES                                       #
#############################################################################################

"Copy vectors from a vector of vectors to a concatenated vector"
function copy_inds(dest, src, inds)
    for i in eachindex(inds)
        dest[inds[i]] = src[i]
    end
end

"Copy constraints to a single concatenated vector"
function copy_constraints!(d, solver::DirectSolver)
    conSet = get_constraints(solver)
    for (i,con) in enumerate(conSet.constraints)
        copy_inds(d, con.vals, solver.con_inds[i])
    end
    return nothing
end

"Copy active set to a single concatenated vector"
function copy_active_set!(a, solver::DirectSolver)
    conSet = get_constraints(solver)
    for i = 1:length(conSet.constraints)
        copy_inds(solver.active_set, conSet.constraints[i].active, solver.con_inds[i])
    end
end

"""Copy constraint Jacobians to given indices in a sparse array
Dispatches on bandedness of the constraint"""
function copy_jacobian!(D, con::ConstraintVals{T,Stage}, cinds, xinds, uinds) where T
    for (i,k) in enumerate(con.inds)
        zind = [xinds[k]; uinds[k]]
        D[cinds[i], zind] .= con.∇c[i]
    end
end

function copy_jacobian!(D, con::ConstraintVals{T,State}, cinds, xinds, uinds) where T
    for (i,k) in enumerate(con.inds)
        D[cinds[i], xinds[k]] .= con.∇c[i]
    end
end

function copy_jacobian!(D, con::ConstraintVals{T,Control}, cinds, xinds, uinds) where T
    for (i,k) in enumerate(con.inds)
        D[cinds[i], uinds[k]] .= con.∇c[i]
    end
end

function copy_jacobian!(D, con::ConstraintVals{T,Coupled}, cinds, xinds, uinds) where T
    for (i,k) in enumerate(con.inds)
        zind = [xinds[k]; uinds[k]; xinds[k+1]; uinds[k+1]]
        D[cinds[i], zind] .= con.∇c[i]
    end
end

function copy_jacobian!(D, con::Union{ConstraintVals{T,Dynamical}, ConstraintVals{T,Coupled,<:DynamicsConstraint{Q}}},
		cinds, xinds, uinds) where {T,Q<:Implicit}
    for (i,k) in enumerate(con.inds)
        zind = [xinds[k]; uinds[k]; xinds[k+1]]
        D[cinds[i], zind] .= con.∇c[i]
    end
end

"Copy all constraint Jacobians to a sparse matrix"
function copy_jacobians!(D, solver::DirectSolver)
    conSet = get_constraints(solver)
    xinds, uinds = primal_partition(solver)
    cinds = solver.con_inds

    for i = 1:length(conSet.constraints)
        copy_jacobian!(D, conSet.constraints[i], cinds[i], xinds, uinds)
    end
    return nothing
end

"Copy constraint Jacobians to linear indices of a vector"
function copy_jacobian!(d::AbstractVector{<:Real}, con::ConstraintVals, linds)
	for (j,k) in enumerate(con.inds)
		inds = linds[j]
		d[inds] = con.∇c[j]
	end
end

"Copy all constraint Jacobians to linear indices of a vector"
function copy_jacobians!(jac::AbstractVector{<:Real}, solver::DirectSolver)
    conSet = get_constraints(solver)
    xinds, uinds = primal_partition(solver)
    cinds = solver.con_inds
    linds = jacobian_linear_inds(solver)

    for (i,con) in enumerate(conSet.constraints)
		copy_jacobian!(jac, con, linds[i])
    end
    return nothing
end

function copy_gradient!(grad_f, E::Vector{<:SizedCostExpansion}, xinds, uinds)
    for k = 1:length(uinds)
        grad_f[xinds[k]] = E[k].x
        grad_f[uinds[k]] = E[k].u
    end
    if length(xinds) != length(uinds)
        grad_f[xinds[end]] = E[end].x
    end
end


text:
# A collection of various physical constants

module constants

export el,me,mp,kB,ep0

const el = 1.60218e-19           # Electron charge in C
const me = 9.094e-31             # Electron mass in kg
const mp = 1.6726e-27            # Proton mass in kg
const kB = 1.3806488e-23         # Boltzmann's constant in SI units (Joule/degK)
const ep0 = 8.854188e-12         # Dielectric constant

end

ext: jl

max_stars_repo_name: gjwilkie/PlasmaTools

max_issues_repo_name: gjwilkie/PlasmaTools

max_forks_repo_name: gjwilkie/PlasmaTools


text:
module SWCs
include("PointObjs.jl")
using .PointObjs
const ONE_UINT32 = UInt32(1)
export SWC

const SWC = Vector{PointObj}

function SWC( swcString::AbstractString )
    swc = SWC()
    for line in split(swcString, "\n")
        if isempty(line)
            continue
        end
        try
            numbers = map(Meta.parse, split(line))
            # construct a point object
            pointObj = PointObj( numbers[2:7]... )
            push!(swc, pointObj)
        catch err
            if occursin(r"^\s*#", line)
                println("comment in swc file: $line")
            elseif isempty(line) || line == " " || line == "\n"
                continue
            else
                @warn("invalid line: $line")
            end
        end
    end
    swc
end

################## properties #######################
function get_node_num(self::SWC) length(self) end

function get_edge_num(self::SWC)
    num_edges = 0
    for pointObj in self
        if PointObjs.get_parent( pointObj ) != -1
            num_edges += 1
        end
    end
    num_edges
end

"""
    get_edges(self::SWC)

get the edges represented as a Vector{NTuple{2,UInt32}}
"""
function get_edges(self::SWC)
    edges = Vector{NTuple{2,UInt32}}()
    for (index, pointObj) in enumerate(self)
        if PointObjs.get_parent( pointObj ) != -1
            push!(edges, ( PointObjs.get_parent( pointObj ), index ))
        end
    end
    edges
end

"""
    get_total_length( self::SWC )
accumulate all the euclidean distance of edges
"""
function get_total_length( self::SWC )
    total_length = Float64(0)
    edges = get_edges( self )
    for e in edges
        total_length += PointObjs.euclidean_distance( self[e[1]], self[e[2]] )
    end
    total_length
end

################## IO ###############################

"""
get binary buffer formatted as neuroglancer nodeNet.

# Binary format
    UInt32: number of vertex
    UInt32: number of edges
    Array{Float32,2}: Nx3 array, xyz coordinates of vertex
    Array{UInt32,2}: Mx2 arrray, node index pair of edges
reference:
https://github.com/seung-lab/neuroglancer/wiki/Skeletons

since there is no node type and radius information in precomputed format,
    we won't have a function to read/load precomputed binary data to SWC.
"""
function get_neuroglancer_precomputed(self::SWC)
    @show get_node_num(self)
    @show get_edge_num(self)
    # total number of bytes
    num_bytes = 4 + 4 + 4*3*get_node_num(self) + 4*2*get_edge_num(self)
    buffer = IOBuffer( read=false, write=true, maxsize=num_bytes )
    # write the number of vertex, and edges
    write(buffer, UInt32(get_node_num(self)))
    write(buffer, UInt32(get_edge_num(self)))
    # write the node coordinates
    for pointObj in self
        write(buffer, pointObj.x)
        write(buffer, pointObj.y)
        write(buffer, pointObj.z)
    end
    # write the edges
    for edge in get_edges( self )
        # neuroglancer index is 0-based
        write(buffer, UInt32( edge[1]-ONE_UINT32 ))
        write(buffer, UInt32( edge[2]-ONE_UINT32 ))
    end
    bin = Vector{UInt8}(take!(buffer))
    close(buffer)
    return bin
end

function Base.String(self::SWC)
    io = IOBuffer()
    for (index, pointObj) in enumerate(self)
        write(io, "$index $(String(pointObj)) \n")
    end
    str = String(take!(io))
    close(io)
    str
end

function deserialize(data::Vector{UInt8})
    # a pointObj is 21 byte
    @assert mod(length(data), 21) == 0 "the binary file do not match the byte layout of pointObj."
    nodeNum = div(length(data), 21)
    swc = Vector{PointObj}()
    sizehint!(swc, nodeNum)
    for i in 1:nodeNum
        pointObj = PointObj( data[(i-1)*21+1:i*21] )
        push!(swc, pointObj)
    end
    swc
end

"""
    load_swc_bin( fileName::AbstractString )
load binary swc file
"""
function load_swc_bin( fileName::AbstractString )
    read( fileName ) |> deserialize
end

function load_swc(fileName::AbstractString)
    swcString = read( fileName , String)
    SWC( swcString )
end

function load(fileName::AbstractString)
    if endswith(fileName, ".swc")
        return load_swc(fileName)
    elseif endswith(fileName, ".swc.bin")
        return load_swc_bin(fileName)
    else
        error("only support .swc and .swc.bin, but this file is: ", fileName)
    end
end

function serialize(self::SWC)
    io = IOBuffer( read=false, write=true, maxsize=length(self)*21 )
    for pointObj in self
        write(io, PointObjs.serialize(pointObj))
    end
    data = take!(io)
end

"""
    save_swc_bin( self::SWC, fileName::AbstractString )
represent swc file as binary file. the data structure is the same with swc.
"""
function save_swc_bin( self::SWC, fileName::AbstractString )
    data = serialize(self)
    write(fileName, data)
end

function save_swc(self::SWC, file_name::AbstractString)
    f = open(file_name, "w")
    for (index, pointObj) in enumerate(self)
        write(f, "$index $(String(pointObj)) \n")
    end
    close(f)
end

function save(self::SWC, fileName::AbstractString)
    if endswith(fileName, ".swc")
        save_swc(self, fileName)
    elseif endswith(fileName, ".swc.bin")
        save_swc_bin(self, fileName)
    else
        error("only support .swc and .swc.bin, but this file is: ", fileName)
    end
end

#################### manipulate ######################

"""
note that only stretch the coordinates here, not including the radius
since radius was already adjusted in the neighborhood weights
"""
function stretch_coordinates!(self::SWC, expansion::Tuple)
    @assert length(expansion) == 3
    for i in 1:length( self )
        self[i].x       *= expansion[1]
        self[i].y       *= expansion[2]
        self[i].z       *= expansion[3]
        self[i].radius  *= (prod(expansion))^(1/3)
    end
end

"""
stretch the coordinate according to the mip level
normally, we only build mip level at XY plane, not Z
"""
function stretch_coordinates!(self::SWC, mip::Integer)
    stretch_coordinates!(self, (2^(mip-1), 2^(mip-1), 1))
end

function add_offset!(self::SWC, offset::Tuple)
    @assert length(offset) == 3
    for in in length( self )
        self[i].x += offset[1]
        self[i].y += offset[2]
        self[i].z += offset[3]
    end
end

end # module of SWCs


text:
# using Revise
using TopOpt, LinearAlgebra, StatsFuns
using StatsFuns: logsumexp

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 3.0

# ### Define the problem
problem = PointLoadCantilever(Val{:Linear}, (160, 40), (1.0, 1.0), E, v, f)

# ### Parameter settings
V = 0.5 # volume fraction
xmin = 0.0001 # minimum density
steps = 40 # maximum number of penalty steps, delta_p0 = 0.1

# ### Continuation SIMP
x0 = fill(1.0, 160 * 40) # initial design
x = copy(x0)
for p in [1.0, 2.0, 3.0]
    global penalty, stress, filter, result, stress, x
    penalty = TopOpt.PowerPenalty(p)
    solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)
    stress = TopOpt.MicroVonMisesStress(solver)
    filter = DensityFilter(solver; rmin=rmin)
    volfrac = TopOpt.Volume(problem, solver)

    obj = x -> volfrac(filter(x)) - V
    thr = 10 # stress threshold
    constr = x -> begin
        s = stress(filter(x))
        vcat((s .- thr) / 100, logsumexp(s) - log(length(s)) - thr)
    end
    alg = PercivalAlg()
    options = PercivalOptions()
    optimizer = Optimizer(obj, constr, x, alg; options=options)
    simp = SIMP(optimizer, solver, p)
    result = simp(x)
    x = result.topology
end

maximum(stress(filter(x0))) # 0.51
maximum(stress(filter(x))) # 10.01

# ### (Optional) Visualize the result using Makie.jl
# Need to run `using Pkg; Pkg.add(Makie)` first
# ```julia
# using Makie
# using TopOpt.TopOptProblems.Visualization: visualize
# fig = visualize(
#    problem; topology = penalty.(filter(result.topology)), default_exagg_scale = 0.07,
#    scale_range = 10.0, vector_linewidth = 3, vector_arrowsize = 0.5,
# )
# Makie.display(fig)
# ```

#md # ## [Plain Program](@id local-stress-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [local-stress.jl](local_stress.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
module TestLDict

using FactCheck
using MultidimensionalTables

facts("LDict tests") do
  context("constructor tests") do
    d1 = LDict(Dict(:a=>1, :c=>3, :b=>2), [:c, :b, :a])
    d2 = LDict([:c, :b, :a], [3, 2, 1])
    d3 = LDict(:c=>3, :b=>2, :a=>1)
    d4 = LDict((:c,3), (:b,2), (:a,1))
    @fact d1 --> d2
    @fact d2 --> d3
    @fact d3 --> d4
  end
  context("method tests") do
    @fact merge(LDict(:a=>1, :b=>2, :c=>3), LDict(:a=>10, :d=>15)) -->
      LDict(:a=>10, :b=>2, :c=>3, :d=>15)
    @fact merge(LDict(:a=>1, :b=>2, :c=>3), LDict(:a=>10, :d=>"15")) -->
      LDict(:a=>10, :b=>2, :c=>3, :d=>"15")
    @fact delete(LDict(:a=>[10,11,12], 10=>3, :c=>5.0, :d=>1), :a, 10) -->
      LDict(:c=>5.0, :d=>1)
    @fact pick(LDict(:a=>[10,11,12], 10=>3, :c=>5.0, :d=>1), [:d, :c]) -->
      LDict(:d=>1, :c=>5.0)
  end
end

end


text:
# This duplicates some functionality from Dolo.jl. Should be moved to Dolang.

sanitize(s, dynvars::Array{Symbol,1}) = s
sanitize(s::Symbol, dynvars::Array{Symbol,1})= s in dynvars ? :($s(0)) : s
function sanitize(eq::Expr, dynvars::Array{Symbol,1})
    if eq.head == :(=)
        return sanitize( :( $(eq.args[1])==$( eq.args[2] ) ), dynvars )
    else
        if eq.head == :call
            if eq.args[1] in dynvars
                return eq
            else
                return Expr(eq.head, [sanitize(e, dynvars) for e in eq.args]...)
            end
        else
            head = eq.head
            args = eq.args
            return Expr(eq.head, [sanitize(e, dynvars) for e in args]...)
        end
    end
end


text:
using CoordinateDescent

include(joinpath(@__DIR__, "..", "benchmark", "locpoly_bench.jl"))

tests = [
  "atom_iterator",
  "coordinate_descent",
  "lasso",
  "varying_coefficient_lasso"
]

for t in tests
	f = "$t.jl"
	println("* running $f ...")
	include(f)
end


text:
module MaxHelpingHandDisableIcePhysicsTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/DisableIcePhysicsTrigger" DisableIcePhysicsTrigger(x::Integer, y::Integer,
    width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, disableIcePhysics::Bool=true)

const placements = Ahorn.PlacementDict(
    "Disable Ice Physics (max480's Helping Hand)" => Ahorn.EntityPlacement(
        DisableIcePhysicsTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# Integer arithmetic functions

for op ∈ (:(~), :(+), :(-))
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline $rename(s1::IntegerTypes) = $op(s1)
        @vectordef $rename function Base.$op(v1) where {W,T<:IntegerTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1))
        end
    end
end
# @inline vnot(s1::Bool) = !s1
# @inline vnot(v1::Vec{N,Bool}) where {N} = vbitwise_not(v1)
# @inline vnot(v1::SVec{N,Bool}) where {N} = SVec(vbitwise_not(extract_data(v1)))
# @inline Base.:!(v1::SVec{N,Bool}) where {N} = SVec(vbitwise_not(extract_data(v1)))

@inline vabs(s1::IntTypes) = abs(s1)
@inline function vabs(v1::_Vec{W,T}) where {W,T<:IntTypes}
    # s = -vbroadcast(Vec{W,T}, signbit(v1))
    s = vright_bitshift(v1, Val{8*sizeof(T)}())
    # Note: -v1 == ~v1 + 1
    vsub(vxor(s, v1), s)
end
@inline vabs(v1::SVec{W,T}) where {W,T<:IntTypes} = SVec(vabs(extract_data(v1)))
@inline Base.abs(v1::SVec{W,T}) where {W,T<:IntTypes} = SVec(vabs(extract_data(v1)))

@inline vabs(v1::_Vec{W,T}) where {W,T<:UIntTypes} = v1
@inline vabs(v1::AbstractStructVec{W,T}) where {W,T<:UIntTypes} = v1
@inline Base.abs(v1::SVec{W,T}) where {W,T<:UIntTypes} = v1

# TODO: Try T(v1>0) - T(v1<0)
#       use a shift for v1<0
#       evaluate v1>0 as -v1<0 ?

@inline vsign(s1::IntTypes) = sign(s1)
@inline vsign(v1::_Vec{W,T}) where {W,T<:IntTypes} =
    vifelse(visequal(v1, vbroadcast(_Vec{W,T},0)), vbroadcast(_Vec{W,T},0),
        vifelse(vless(v1, vbroadcast(_Vec{W,T},0)), vbroadcast(_Vec{W,T},-1), vbroadcast(_Vec{W,T},1)))
@inline vsign(v1::_Vec{W,T}) where {W,T<:UIntTypes} =
    vifelse(visequal(v1, vbroadcast(_Vec{W,T},0)), vbroadcast(_Vec{W,T},0), vbroadcast(_Vec{W,T},1))
@inline vsign(v1::SVec) = SVec(vsign(extract_data(v1)))
@inline Base.sign(v1::SVec) = SVec(vsign(extract_data(v1)))

@inline vsignbit(s1::IntegerTypes) = signbit(s1)
@inline vsignbit(v1::_Vec{W,T}) where {W,T<:IntTypes} = vless(v1, vbroadcast(_Vec{W,T}, 0))
@inline vsignbit(v1::SVec{W,T}) where {W,T<:IntTypes} = SVec(vsignbit(extract_data(v1)))
@inline Base.signbit(v1::SVec{W,T}) where {W,T<:IntTypes} = SVec(vsignbit(extract_data(v1)))

@inline vsignbit(v1::_Vec{W,T}) where {W,T<:UIntTypes} = vsignbit(SVec(v1))
@inline vsignbit(v1::SVec{W,T}) where {W,T<:UIntTypes} = zero(Mask{W})
@inline Base.signbit(v1::SVec{W,T}) where {W,T<:UIntTypes} = zero(Mask{W})

for op ∈ (:(&), :(|), :(⊻), :(+), :(-), :(*), :(÷), :(%) )
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline $rename(s1::IntegerTypes, s2::IntegerTypes) = $op(s1, s2)

        @vectordef $rename function Base.$op(v1, v2) where {W,T<:IntegerTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2))
        end

        # @inline $rename(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypes} =
        #     llvmwrap(Val{$(QuoteNode(op))}(), v1, v2)
    end
end

@inline Base.:%(v::SVec{W,I}, ::Type{I}) where {W,I} = v
@inline Base.:%(v::SVec{W,I}, ::Type{Vec{W,I}}) where {W,I} = v
@inline Base.:%(v::SVec{W,I}, ::Type{SVec{W,I}}) where {W,I} = v
# @inline Base.:%(v::Vec{W,I}, ::Type{SVec{W,I}}) where {W,I} = SVec(v)
Base.:%(v::SVec{W,I1}, ::Type{I2}) where {W,I1,I2} =  SVec(vconvert(Vec{W,I2}, extract_data(v)))
@inline Base.:%(v::SVec{W,I1}, ::Type{SVec{W,I2}}) where {W,I1,I2} = v % I2
@inline Base.:%(v::SVec{W,I1}, ::Type{Vec{W,I2}}) where {W,I1,I2} = v % I2

@inline vrem(v1::V, ::Type{V}) where {V <: AbstractSIMDVector} = v1

@inline vcopysign(s1::IntegerTypes, s2::IntegerTypes) = copysign(s1, s2)
@inline vcopysign(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T<:IntTypes} =
    vifelse(vsignbit(v2), -abs(v1), abs(v1))
# @inline function vcopysign(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T}
    # SVec(vcopysign(extract_data(v1), extract_data(v2)))
# end
@inline function Base.copysign(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T}
    SVec(vcopysign(extract_data(v1), extract_data(v2)))
end
@inline vcopysign(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T<:UIntTypes} = v1
@inline vcopysign(v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}) where {W,T<:UIntTypes} = v1
@inline Base.copysign(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T<:UIntTypes} = v1

@inline vflipsign(s1::IntegerTypes, s2::IntegerTypes) = flipsign(s1, s2)
@inline vflipsign(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T<:IntTypes} = vifelse(vsignbit(v2), -v1, v1)
@inline function vflipsign(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T<:IntTypes}
    SVec(vflipsign(extract_data(v1), extract_data(V2)))
end
@inline function Base.flipsign(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T<:IntTypes}
    SVec(vflipsign(extract_data(v1), extract_data(V2)))
end

@inline vflipsign(v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}) where {W,T<:UIntTypes} = v1
@inline vflipsign(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T<:UIntTypes} = v1
@inline Base.flipsign(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T<:UIntTypes} = v1

@inline vmax(s1::IntegerTypes, s2::IntegerTypes) = max(s1, s2)
@inline vmin(s1::IntegerTypes, s2::IntegerTypes) = min(s1, s2)
@inline vmax(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T<:IntegerTypes} = vifelse(vless(v2, v1), v1, v2)
@inline vmin(v1::_Vec{W,T}, v2::_Vec{W,T}) where {W,T<:IntegerTypes} = vifelse(vless(v1, v2), v1, v2)
@inline function Base.max(v1::AbstractStructVec, v2::AbstractStructVec)
    V = promote_vtype(typeof(v1), typeof(v2))
    SVec(vmax(extract_data(vconvert(V, v1)), extract_data(vconvert(V, v2))))
end
@inline function Base.max(v1, v2::AbstractStructVec)
    V = promote_vtype(typeof(v1), typeof(v2))
    SVec(vmax(extract_data(vconvert(V, v1)), extract_data(vconvert(V, v2))))
end
@inline function Base.max(v1::AbstractStructVec, v2)
    V = promote_vtype(typeof(v1), typeof(v2))
    SVec(vmax(extract_data(vconvert(V, v1)), extract_data(vconvert(V, v2))))
end
@inline function Base.min(v1::AbstractStructVec, v2::AbstractStructVec)
    V = promote_vtype(typeof(v1), typeof(v2))
    SVec(vmin(extract_data(vconvert(V, v1)), extract_data(vconvert(V, v2))))
end
@inline function Base.min(v1, v2::AbstractStructVec)
    V = promote_vtype(typeof(v1), typeof(v2))
    SVec(vmin(extract_data(vconvert(V, v1)), extract_data(vconvert(V, v2))))
end
@inline function Base.min(v1::AbstractStructVec, v2)
    V = promote_vtype(typeof(v1), typeof(v2))
    SVec(vmin(extract_data(vconvert(V, v1)), extract_data(vconvert(V, v2))))
end

vmuladd(s1::ScalarTypes, s2::ScalarTypes, s3::ScalarTypes) = muladd(s1,s2,s3)
@inline function vmuladd(v1::_Vec{W,T}, v2::_Vec{W,T}, v3::_Vec{W,T}) where {W,T<:IntegerTypes}
    vadd(vmul(v1,v2),v3)
end
@inline function vmuladd(
    v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}, v3::AbstractStructVec{W,T}
) where {W,T<:IntegerTypes}
    SVec(vadd(vmul(extract_data(v1),extract_data(v2)),extract_data(v3)))
end
@inline function Base.muladd(
    v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}, v3::AbstractStructVec{W,T}
) where {W,T<:IntegerTypes}
    SVec(vadd(vmul(extract_data(v1),extract_data(v2)),extract_data(v3)))
end

vfma(s1::ScalarTypes, s2::ScalarTypes, s3::ScalarTypes) = fma(s1,s2,s3)
@inline function vfma(v1::_Vec{W,T}, v2::_Vec{W,T}, v3::_Vec{W,T}) where {W,T<:IntegerTypes}
    vadd(vmul(v1,v2),v3)
end
@inline function vfma(
    v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}, v3::AbstractStructVec{W,T}
) where {W,T<:IntegerTypes}
    SVec(vadd(vmul(extract_data(v1),extract_data(v2)),extract_data(v3)))
end
@inline function Base.fma(
    v1::AbstractStructVec{W,T}, v2::AbstractStructVec{W,T}, v3::AbstractStructVec{W,T}
) where {W,T<:IntegerTypes}
    SVec(vadd(vmul(extract_data(v1),extract_data(v2)),extract_data(v3)))
end
# TODO: Handle negative shift counts
#       use vifelse
#       ensure vifelse is efficient
for op ∈ (:(<<), :(>>), :(>>>))
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline $rename(s1::IntegerTypes, s2::IntegerTypes) = $op(s1, s2)

        @vectordef $rename function Base.$op(v1, ::Val{I}) where {W,T<:IntegerTypes,I}
            llvmwrapshift(Val{$(QuoteNode(op))}(), extract_data(v1), Val{I}())
        end
        @vectordef $rename function Base.$op(v1, x2::Unsigned) where {W,T<:IntegerTypes}
            llvmwrapshift(Val{$(QuoteNode(op))}(), extract_data(v1), x2)
        end
        @vectordef $rename function Base.$op(v1, x2::Int) where {W,T<:IntegerTypes}
            llvmwrapshift(Val{$(QuoteNode(op))}(), extract_data(v1), x2)
        end
        @vectordef $rename function Base.$op(v1, x2::Integer) where {W,T<:IntegerTypes}
            llvmwrapshift(Val{$(QuoteNode(op))}(), extract_data(v1), x2)
        end
        @vectordef $rename function Base.$op(v1, v2) where {W,T<:IntegerTypes}
            llvmwrapshift(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2))
        end
        @vectordef $rename function Base.$op(x1::T, v2) where {W,T<:IntegerTypes}
            $rename(vbroadcast(Vec{W,I}, x1), extract_data(v2))
        end
        @vectordef $rename function Base.$op(x1::I, v2) where {W,I<:Integer,T<:IntegerTypes}
            $rename(vbroadcast(Vec{W,I}, x1), vconvert(Vec{W,I},extract_data(v2)))
        end
        @eval @inline Base.$op(v1::SVec{W,I1}, v2::SVec{W,I2}) where {W,I1,I2<:IntegerTypes} = SVec(llvmwrapshift(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2)))
    end
end
@inline function Base.:(<<)(v1::AbstractStructVec{W,I1}, v2::AbstractStructVec{W,I2}) where {W, I1 <: Unsigned, I2}
    v2pos = v2 > 0
    v2abs = vconvert(SVec{W, I1}, abs(v2))
    vifelse(v2pos, v1 << v2abs, v1 >> v2abs)
end

# Promote scalars of all IntegerTypes to vectors of IntegerTypes, leaving the
# vector type unchanged

for op ∈ (
        :(==), :(!=), :(<), :(<=), :(>), :(>=),
        :(&), :(|), :(⊻), :(+), :(-), :(*),
        :copysign, :div, :flipsign, :max, :min, :rem
    )
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline $rename(s1::ScalarTypes, s2::ScalarTypes) = $op(s1, s2)

        @vectordef $rename function Base.$op(s1::IntegerTypes, v2) where {W, T <: IntegerTypes}
            $rename(vbroadcast(typeof(v2), s1), extract_data(v2))
        end
        @vectordef $rename function Base.$op(v1, s2::IntegerTypes) where {W, T <: IntegerTypes}
            $rename(extract_data(v1), vbroadcast(typeof(v1), s2))
        end
    end
end
@inline evadd(v1,v2) = vadd(v1,v2)
@inline evsub(v1,v2) = vsub(v1,v2)
@inline evmul(v1,v2) = vmul(v1,v2)
@generated function vpmaddwd(a::NTuple{W,Core.VecElement{Int16}}, b::NTuple{W,Core.VecElement{Int16}}) where {W}
    Wh = W >>> 1
    @assert 2Wh == W
    @assert (VectorizationBase.REGISTER_SIZE >> 1) ≥ W
    S = W * 16
    # decl = "@llvm.x86.avx512.pmaddw.d.512"
    instr = "@llvm.x86.avx512.pmaddw.d.$S"
    decl = "declare <$Wh x i32> $instr(<32 x i16>, <32 x i16>)"
    instrs = String[
        "%res = call <$Wh x i32> $instr(<$W x i16> %0, <$W x i16> %1)",
        "ret <$Wh x i32> %res"
    ]
    quote
        $(Expr(:meta,:inline))
        llvmcall(
            $((decl,join(instrs,"\n"))),
            NTuple{$Wh,Core.VecElement{Int32}},
            Tuple{NTuple{$W,Core.VecElement{Int16}},NTuple{$W,Core.VecElement{Int16}}},
            a, b
        )
    end
end

@generated function vleading_zeros(v::_Vec{_W,I}) where {_W, I <: Integer}
    W = _W + 1
    typ = "i$(8sizeof(I))"
    vtyp = "<$W x $typ>"
    instr = "@llvm.ctlz.v$(W)$(typ)"
    decl = "declare $vtyp $instr($vtyp, i1)"
    instrs = String[
        "%res = call $vtyp $instr($vtyp %0, i1 1)"
        "ret $vtyp %res"
    ]
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $((decl, join(instrs, "\n"))),
            Vec{$W,$I}, Tuple{Vec{$W,$I}}, v
        )
    end
end
@inline vleading_zeros(v::SVec{<:Any,<:Integer}) = SVec(vleading_zeros(extract_data(v)))
@inline Base.leading_zeros(v::SVec{<:Any,<:Integer}) = SVec(vleading_zeros(extract_data(v)))

@generated function vtrailing_zeros(v::_Vec{_W,I}) where {_W, I <: Integer}
    W = _W + 1
    typ = "i$(8sizeof(I))"
    vtyp = "<$W x $typ>"
    instr = "@llvm.cttz.v$(W)$(typ)"
    decl = "declare $vtyp $instr($vtyp, i1)"
    instrs = String[
        "%res = call $vtyp $instr($vtyp %0, i1 1)"
        "ret $vtyp %res"
    ]
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $((decl, join(instrs, "\n"))),
            Vec{$W,$I}, Tuple{Vec{$W,$I}}, v
        )
    end
end
@inline vtrailing_zeros(v::SVec{<:Any,<:Integer}) = SVec(vtrailing_zeros(extract_data(v)))
@inline Base.trailing_zeros(v::SVec{<:Any,<:Integer}) = SVec(vtrailing_zeros(extract_data(v)))

# @generated function vadd(m::Mask{W,U}, v::Vec{W,I}) where {W,U,I<:Integer}
#     vityp = "<$W x i$(8sizeof(I))>"
#     instrs = String[]
#     if 8sizeof(U) == W
#         push!(instrs, "%bitvector = bitcast i$(W) %0 to <$W x i1>")
#     else
#         push!(instrs, "%truncint = trunc i$(8sizeof(U)) %0 to i$(W)")
#         push!(instrs, "%bitvector = bitcast i$(W) %truncint to <$W x i1>")
#     end
#     push!(instrs, "%zextv = zext <$W x i1> %bitvector to $vityp")
#     push!(instrs, "%res = add $vityp %zextv, %1")
#     push!(instrs, "ret $vityp %res")
#     quote
#         $(Expr(:meta, :inline))
#         SVec(llvmcall(
#             $(join(instrs, "\n")),
#             Vec{$W,$I}, Tuple{$U, Vec{$W,$I}}, m.u, v
#         ))
#     end
# end
# @inline vadd(m::Mask{W}, v::AbstractStructVec{W,I}) where {W,I<:IntegerTypes} = vadd(m, extract_data(v))

@inline Base.:(<<)(i::_MM{W}, j::SVec{W,T}) where {W,T<:Integer} = svrange(_MM{W}(i.i %T)) << j
@inline Base.:(>>)(i::_MM{W}, j::SVec{W,T}) where {W,T<:Integer} = svrange(_MM{W}(i.i %T)) >> j
@inline Base.:(>>>)(i::_MM{W}, j::SVec{W,T}) where {W,T<:Integer} = svrange(_MM{W}(i.i %T)) >>> j
@inline Base.:(<<)(i::SVec{W,T}, j::_MM{W}) where {W,T<:Integer} = i << svrange(_MM{W}(j.i %T))
@inline Base.:(>>)(i::SVec{W,T}, j::_MM{W}) where {W,T<:Integer} = i >> svrange(_MM{W}(j.i %T))
@inline Base.:(>>>)(i::SVec{W,T}, j::_MM{W}) where {W,T<:Integer} = i >>> svrange(_MM{W}(j.i %T))
@inline Base.:(<<)(i::_MM{W}, j::_MM{W}) where {W} = svrange(i) << svrange(j)
@inline Base.:(>>)(i::_MM{W}, j::_MM{W}) where {W} = svrange(i) >> svrange(j)
@inline Base.:(>>>)(i::_MM{W}, j::_MM{W}) where {W} = svrange(i) >>> svrange(j)

makesvec(v::AbstractSIMDVector) = SVec(v)
makesvec(s) = s
@inline function funnel_shift_left(v1, v2, v3)
    vp1, vp2, vp3 = vpromote(v1, v2, v3)
    makesvec(llvmwrap(Val{:funnel_shift_left}(), extract_data(vp1), extract_data(vp2), extract_data(vp3)))
end
@inline function funnel_shift_right(v1, v2, v3)
    vp1, vp2, vp3 = vpromote(v1, v2, v3)
    makesvec(llvmwrap(Val{:funnel_shift_right}(), extract_data(vp1), extract_data(vp2), extract_data(vp3)))
end
@inline rotate_left(a, b) = funnel_shift_left(a, a, b)
@inline rotate_right(a, b) = funnel_shift_right(a, a, b)


text:
module NamedDims
using Base: @propagate_inbounds
using Base.Broadcast:
    Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle, Unknown
using LinearAlgebra
using Pkg
using Requires
using Statistics

export NamedDimsArray, dim, rename, unname

function __init__()
    # NOTE: NamedDims is only compatible with Tracker v0.2.2; but no nice way to enforce that.
    @require Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("tracker_compat.jl")
end

# We use CoVector to workout if we are taking the tranpose of a tranpose etc
const CoVector = Union{Adjoint{<:Any, <:AbstractVector}, Transpose{<:Any, <:AbstractVector}}

include("name_core.jl")
include("wrapper_array.jl")
include("broadcasting.jl")
include("functions.jl")
include("functions_dims.jl")
include("functions_math.jl")

end # module


text:
# ---
# title: 1556. Thousand Separator
# id: problem1556
# author: Indigo
# date: 2021-01-22
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/thousand-separator/description/>
# hidden: true
# ---
#
# Given an integer `n`, add a dot (".") as the thousands separator and return it
# in string format.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 987
#     Output: "987"
#
#
# **Example 2:**
#
#
#
#     Input: n = 1234
#     Output: "1.234"
#
#
# **Example 3:**
#
#
#
#     Input: n = 123456789
#     Output: "123.456.789"
#
#
# **Example 4:**
#
#
#
#     Input: n = 0
#     Output: "0"
#
#
#
#
# **Constraints:**
#
#   * `0 <= n < 2^31`
#
#
## @lc code=start
using LeetCode

function thousand_separator(n::Int)
    s = string(n)
    r = length(s) % 3
    ed = 0
    if r == 0
        r += 3
        ed += 1
    end
    res = s[1:r]
    for i in 1 : length(s) ÷ 3 - ed
        res *= "." * s[r + i * 3 - 2 : r + i * 3]
    end
    res
end
## @lc code=end


text:
module HyrodiumHome

using Documenter

"""
Generate language list
The following script is a modified version of the Documenter.jl script.
"""
function Documenter.Writers.HTMLWriter.generate_version_file(versionfile::AbstractString, entries, symlinks = [])
    script = """
    var DOC_VERSIONS = [
      "en",
      "ja",
    ];
    var DOCUMENTER_NEWEST = "en";
    var DOCUMENTER_STABLE = "en";
    """
    write(versionfile, script)
end

end


text:
using Test, ChooseOptimizer, Cbc

@test get_solver_name() == Cbc


text:
using TriMatrices
using TriMatrices.Testing

@testset "make_test_matrix" begin
	@test make_test_matrix(TriLower{true}(), 3) == (
		[10, 20, 30, 40, 50, 60],
		[10  0  0;
		 20 30  0;
		 40 50 60]
	)
	@test make_test_matrix(TriLower{false}(), 3) == (
		[10, 20, 30],
		[-1  0  0;
		 10 -1  0;
		 20 30 -1]
	)

	@test make_test_matrix(TriUpper{true}(), 3) == (
		[10, 20, 30, 40, 50, 60],
		[10 20 40;
		 0 30 50;
		 0  0 60]
	)
	@test make_test_matrix(TriUpper{false}(), 3) == (
		[10, 20, 30],
		[-1 10 20;
		 0 -1 30;
		 0  0 -1]
	)

	@test make_test_matrix(TriSymmetric{true}(), 3) == (
		[10, 20, 30, 40, 50, 60],
		[10 20 40;
		 20 30 50;
		 40 50 60]
	)
	@test make_test_matrix(TriSymmetric{false}(), 3) == (
		[10, 20, 30],
		[-1 10 20;
		 10 -1 30;
		 20 30 -1]
	)
end


text:
function KomaUI(;dark=true,frame=true, phantom_mode="2D", sim=Dict{String,Any}(), rec=Dict{Symbol,Any}())
## ASSETS
path = @__DIR__
assets = AssetRegistry.register(dirname(path*"/ui/assets/"))
scripts = AssetRegistry.register(dirname(path*"/ui/scripts/"))
css = AssetRegistry.register(dirname(path*"/ui/css/"))
# Assets
background = assets*"/spiral-bg.svg" #In Windows joinpath causes problems "/assetserver/...-assets\Logo.png"
logo = joinpath(assets, "Logo_dark.svg")
loading = joinpath(assets, "Loading.gif")
# JS
bsjs = joinpath(scripts, "bootstrap.bundle.min.js") #this already has Popper
bscss = joinpath(css,"bootstrap.min.css")
bsiconcss = joinpath(css,"bootstrap-icons.css")
jquery = joinpath(scripts,"jquery-3.4.1.slim.min.js")
# mathjaxsetup = joinpath(scripts, "mathjaxsetup.js")
# KaTeX
katexrender = joinpath(scripts, "auto-render.min.js")
katexjs = joinpath(scripts,"katex.min.js")
katexcss = joinpath(css,"katex.min.css")
# User defined JS and CSS
customcss = joinpath(css,"custom.css")
customjs = joinpath(scripts,"custom.js")
customjs2 = joinpath(scripts,"custom2.js")
sidebarcss = joinpath(css,"sidebars.css")
# Custom icons
icons = joinpath(css,"icons.css")
## WINDOW
global w = Blink.Window(Dict(
    "title"=>"KomaUI",
    "autoHideMenuBar"=>true,
    "frame"=>frame, #removes title bar
    "node-integration" => true,
    :icon=>path*"/ui/assets/Logo_icon.png",
    "width"=>1200,
    "height"=>800
    ),async=false);
## LOADING BAR
buffericon = """

"""
progressbar = """

0%

"""
## NAV BAR
sidebar = open(f->read(f, String), path*"/ui/html/sidebar.html")
sidebar = replace(sidebar, "LOGO"=>logo)
## CONTENT
index = open(f->read(f, String), path*"/ui/html/index.html")
index = replace(index, "BACKGROUND_IMAGE"=>background)
## CSS
loadcss!(w, bscss)
loadcss!(w, bsiconcss)
loadcss!(w, customcss)
loadcss!(w, sidebarcss)
# KATEX
loadcss!(w, katexcss)
loadjs!(w, katexjs)
loadjs!(w, katexrender)
# JQUERY, BOOSTRAP JS
loadjs!(w, customjs)    #must be before jquery
loadjs!(w, jquery)
loadjs!(w, bsjs)        #after jquery
loadjs!(w, customjs2)   #must be after jquery
# LOAD ICONS
loadcss!(w, icons)
## INITAL VARIABLES
#PHANTOM init
@info "Loading Phantom (default)"
if phantom_mode == "3D"
    global phantom = brain_phantom3D()
else
    global phantom = brain_phantom2D()
end
phantom.Δw .*= 0
println("Phantom object \"$(phantom.name)\" successfully loaded!")
#SCANNER init
@info "Loading Scanner (default)"
global sys = Scanner()
println("B0 = $(sys.B0) T")
println("Gmax = $(round(sys.Gmax*1e3,digits=2)) mT/m")
println("Smax = $(sys.Smax) mT/m/ms")
#SEQ init
@info "Loading Sequence (default) "
B1 = sys.B1; durRF = π/2/(2π*γ*B1) #90-degree hard excitation pulse
EX = PulseDesigner.RF_hard(B1, durRF, sys; G=[0,0,0])
N = 101
FOV = 23e-2
EPI = PulseDesigner.EPI(FOV, N, sys)
TE = 30e-3
d1 = TE-dur(EPI)/2-dur(EX)
if d1 > 0 DELAY = Delay(d1) end
global seq = d1 > 0 ? EX + DELAY + EPI : EX + EPI
seq.DEF["TE"] = round(d1 > 0 ? TE : TE - d1, digits=4)*1e3
println("EPI successfully loaded! (TE = $(seq.DEF["TE"]) ms)")
#Init
global darkmode = dark
global raw_ismrmrd = RawAcquisitionData(Dict("systemVendor"=>"","encodedSize"=>[2,1,1],
"reconSize"=>[2,1,1],"number_of_samples"=>2,"encodedFOV"=>[100,1,1]), [Profile(AcquisitionHeader(),[0im; 0;;],[0im; 0;;])])
global rawfile = ""
global image =  [0.0im 0.; 0. 0.]
global kspace = [0.0im 0.; 0. 0.]
#Reco
default = Dict{Symbol,Any}(:reco=>"direct") #, :iterations=>10, :λ=>1e-5,:solver=>"admm",:regularization=>"TV")
global recParams = merge(default, rec)
#Simulation
Ncores = Hwloc.num_physical_cores()
default = Dict{String,Any}("Δt"=>1e-3,"Δt_rf"=>1e-4,"gpu"=>has_cuda(),"Nthreads"=>Ncores)
global simParams = merge(default, sim)
#GPUs
if simParams["gpu"]
    @info "Loading GPUs"
    print_gpus()
    device!(0) #By default it uses first GPU, multiGPU not supported yet
end
#OBERSVABLES
global seq_obs = Observable{Sequence}(seq)
global pha_obs = Observable{Phantom}(phantom)
global sig_obs = Observable{RawAcquisitionData}(raw_ismrmrd)
global img_obs = Observable{Any}(image)
## MENU FUNCTIONS
handle(w, "index") do args...
    content!(w, "div#content", index)
 end
handle(w, "pulses_seq") do args...
    include(path*"/ui/PulsesGUI_seq.jl")
end
handle(w, "pulses_kspace") do args...
    include(path*"/ui/PulsesGUI_kspace.jl")
end
handle(w, "pulses_M0") do args...
    include(path*"/ui/PulsesGUI_M0.jl")
end
handle(w, "phantom") do args...
    include(path*"/ui/PhantomGUI.jl")
end
handle(w, "sig") do args...
    include(path*"/ui/SignalGUI.jl")
end
handle(w, "reconstruction_absI") do args...
    include(path*"/ui/ReconGUI_absI.jl")
end
handle(w, "reconstruction_angI") do args...
    include(path*"/ui/ReconGUI_angI.jl")
end
handle(w, "reconstruction_absK") do args...
    include(path*"/ui/ReconGUI_absK.jl")
end
handle(w, "sim_params") do args...
    include(path*"/ui/SimParams_view.jl")
end
handle(w, "rec_params") do args...
    include(path*"/ui/RecParams_view.jl")
end
handle(w, "simulate") do args...
    # @js_ w document.getElementById("simulate!").prop("disabled", true); #Disable button during SIMULATION
    @js_ w (@var progressbar = $progressbar; document.getElementById("simulate!").innerHTML=progressbar)
    #To SequenceGUI
    signal = simulate(phantom, seq, sys; simParams, w)
    #After simulation go to RECON
    @js_ w document.getElementById("simulate!").innerHTML="Simulate!"
    #EXPORT to ISMRMRD -> To SignalGUI
    println("")
    @info "Exporting to ISMRMRD file"
    global raw_ismrmrd = rawSignalToISMRMRD([signal;;],seq;phantom,sys,simParams)
    global rawfile = tempdir()*"/Koma_signal_ismrmrd.h5"
    global sig_obs[] = raw_ismrmrd
    fout = ISMRMRDFile(rawfile)
    save(fout, raw_ismrmrd)
    println("$rawfile successfully saved!")
    #Message
    @js_ w Toasty("1", "Simulation successfull" ,"""

            <button class="btn btn-dark btn-circle btn-circle-sm m-1" onclick="Blink.msg('sig', 1)"></button>
            Updating Raw signal plots ...

            <button class="btn btn-success btn-circle btn-circle-sm m-1" onclick="Blink.msg('recon', 1)"></button>
            Ready to reconstruct?

    """)
    # @js_ w document.getElementById("simulate!").prop("disabled", false); #Re-enable button
    # @js_ w (@var button = document.getElementById("recon!"); @var bsButton = @new bootstrap.Button(button); vsButton.toggle())
end
handle(w, "recon") do args...
    @js_ w (@var buffericon = $buffericon; document.getElementById("recon!").innerHTML=buffericon)
    #IMPORT ISMRMRD raw data
    raw_ismrmrd.profiles = raw_ismrmrd.profiles[getproperty.(getproperty.(raw_ismrmrd.profiles, :head), :flags) .!= 268435456] #Extra profile in JEMRIS simulations
    acqData = AcquisitionData(raw_ismrmrd)
    acqData.traj[1].circular = false #Removing circular window
    acqData.traj[1].nodes = acqData.traj[1].nodes[1:2,:] ./ maximum(2*abs.(acqData.traj[1].nodes[:])) #Normalize k-space to -.5 to .5 for NUFFT
    Nx, Ny = raw_ismrmrd.params["reconSize"][1:2]
    recParams[:reconSize] = (Nx, Ny)
    recParams[:densityWeighting] = true
    #Reconstruction
    println("")
    @info "Running reconstruction of $rawfile"
    aux = @time MRIReco.reconstruction(acqData, recParams)
    global image  = reshape(aux.data,Nx,Ny,:)
    global kspace = fftc(reshape(aux.data,Nx,Ny,:))
    # global img_obs[] = image
    println("Reconstruction successfull!")
    #After Recon go to Image
    @js_ w document.getElementById("recon!").innerHTML="Reconstruct!"
    @js_ w Toasty("2", "Reconstruction successfull" ,"""

            <button class="btn btn-dark btn-circle btn-circle-sm m-1" onclick="Blink.msg('reconstruction_absI', 1)"></button>
            Updating Reconstruction plots ...

    """
    )
end
handle(w, "close") do args...
    global darkmode = nothing

    global phantom = nothing
    global seq = nothing
    global sys = nothing

    global raw_ismrmrd = nothing
    global rawfile = nothing

    global image = nothing
    global kspace = nothing

    global simParams = nothing
    global recParams = nothing

    global seq_obs = nothing
    global pha_obs = nothing
    global sig_obs = nothing
    global img_obs = nothing
    close(w)
end
#Update GUI's home
w = body!(w, *(sidebar,index), async=false)
if darkmode
    @js_ w document.getElementById("main").style="background-color:rgb(13,16,17);"
end
#Sequence observable
load_seq = filepicker(".seq (Pulseq)/.seqk (Koma)"; accept=".seq,.seqk")
map!(f->if f!="" #Assigning function of data when load button (filepicker) is changed
            if splitext(f)[end]==".seqk" #Koma
                global seq = JLD2.load(FileIO.File{FileIO.DataFormat{:JLD2}}(f),"seq")
            elseif splitext(f)[end]==".seq" #Pulseq
                global seq = read_seq(f) #Pulseq read
            end
            @js_ w (@var name = $(basename(f));
            Toasty("0", "Loaded "+name+" successfully", """

                    <button class="btn btn-dark btn-circle btn-circle-sm m-1" onclick="Blink.msg('pulses_seq', 1)"></button>
                    Updating Sequence plots ...

                    <button class="btn btn-primary btn-circle btn-circle-sm m-1" onclick="Blink.msg('simulate', 1)"></button>
                    Ready to simulate?

            """))
            seq
        else
            seq #default sequence
        end
    , seq_obs, load_seq)
w = content!(w, "#seqfilepicker", load_seq, async=false)
#Phantom observable
load_pha = filepicker(".phantom (Koma)/.h5 (JEMRIS)"; accept=".phantom,.h5")
map!(f->if f!="" #Assigning function of data when load button (filepicker) is changed
            if splitext(f)[end]==".phantom" #Koma
                global phantom = JLD2.load(FileIO.File{FileIO.DataFormat{:JLD2}}(f),"phantom")
            elseif splitext(f)[end]==".h5" #JEMRIS
                global phantom = read_phantom_jemris(f)
            end
            @js_ w (@var name = $(basename(f));
            Toasty("0", "Loaded "+name+" successfully", """

                    <button class="btn btn-dark btn-circle btn-circle-sm m-1" onclick="Blink.msg('phantom', 1)"></button>
                    Updating Phantom plots ...

                    <button class="btn btn-primary btn-circle btn-circle-sm m-1" onclick="Blink.msg('simulate', 1)"></button>
                    Ready to simulate?

            """))
            phantom
        else
            phantom #default
        end
    , pha_obs, load_pha)
w = content!(w, "#phafilepicker", load_pha, async=false)
#Signal observable
load_sig = filepicker(".h5 (ISMRMRD)"; accept=".h5")
map!(f->if f!="" #Assigning function of data when load button (filepicker) is changed
            fraw = ISMRMRDFile(f)
            global rawfile = f
            global raw_ismrmrd = RawAcquisitionData(fraw)

            not_Koma = raw_ismrmrd.params["systemVendor"] != "Koma.jl"
            if not_Koma
                @warn "ISMRMRD files generated externally could cause problems during the reconstruction. We are currently improving compatibility."
            end

            @js_ w (@var name = $(basename(f));
            Toasty("0", "Loaded "+name+" successfully", """

                    <button class="btn btn-dark btn-circle btn-circle-sm m-1" onclick="Blink.msg('sig', 1)"></button>
                    Updating Raw data plots ...

                    <button class="btn btn-success btn-circle btn-circle-sm m-1" onclick="Blink.msg('recon', 1)"></button>
                    Ready to reconstruct?

            """))
            raw_ismrmrd
        else
            raw_ismrmrd #default
        end
    , sig_obs, load_sig)
w = content!(w, "#sigfilepicker", load_sig, async=false)
#Update Koma version
version =  string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]))
content!(w, "#version", version, async=false)
@info "Currently using Koma v$version"

nothing
end


text:
using Images
using ImageView
using Color
using FixedPointNumbers

function nonzero_indexes{T <: Number}(mat::Matrix{T})
    idxs = (Int, Int)[]
    for i=1:size(mat, 1)
        for j=1:size(mat, 2)
            if mat[i, j] != 0
                push!(idxs, (i, j))
            end
        end
    end
    return idxs
end

function collect_nonzeros{T <: Number}(mat::Matrix{T}, idxs::Array{(Int, Int)})
    arr = zeros(T, length(idxs))
    k = 1
    for (i, j) in idxs
        arr[k] = mat[i, j]
        k += 1
    end
    return arr
end

function map_nonzeros{T <: Number}(matsize::(Int, Int), arr::Array{T, 1},
                      idxs::Array{(Int, Int)})
    mat = zeros(T, matsize)
    k = 1
    for (i, j) in idxs
        mat[i, j] = arr[k]
        k += 1
    end
    return mat
end

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
function mpccmodel_setup_symdiff_sparse(config::MPCCModelConfig)

    @unpack dimspec, defn, fns = config
    @unpack n, q, l, me, mi, r, s = config.dimspec

    s_x = config.x
    s_pr = config.pr
    s_ps = config.ps

    # We need to setup the num arrays for the sparse jacobian and hessians
    # Note that although SparseMatrixCSC is used consistently below, most of these
    # are vectors; Symbolics spits out the mtrx type and can't be bothered converting

    # grad entries
    gradf_sparse_num = transpose( Symbolics.sparsejacobian([ defn.f ], s_x) )
    gradce_sparse_num = transpose( Symbolics.sparsejacobian( defn.ce, s_x ) )
    gradci_sparse_num = transpose( Symbolics.sparsejacobian( defn.ci, s_x ) )
    gradF_sparse_num = Matrix{SparseMatrixCSC}(undef, l, q)
    for lp_q=1:q
        for lp_l=1:l
            gradF_sparse_num[lp_l, lp_q] = transpose( Symbolics.sparsejacobian( [ defn.F[lp_l, lp_q] ], s_x ) )
        end
    end

    println("gradf: ", gradf_sparse_num)
    println("gradce: ", gradce_sparse_num)
    println("gradci: ", gradci_sparse_num)
    println("gradF: ", gradF_sparse_num)

    # hess entries
    hessf_sparse_num = Symbolics.sparsehessian(defn.f, s_x)
    hessce_sparse_num = Vector{SparseMatrixCSC}(undef, me)
    for lp_ce=1:me
        hessce_sparse_num[lp_ce] = Symbolics.sparsehessian(defn.ce[lp_ce], s_x)
    end
    hessci_sparse_num = Vector{SparseMatrixCSC}(undef, mi)
    for lp_ci=1:mi
        hessci_sparse_num[lp_ci] = Symbolics.sparsehessian(defn.ci[lp_ci], s_x)
    end
    hessF_sparse_num = Matrix{SparseMatrixCSC}(undef, l, q)
    for lp_q=1:q
        for lp_l=1:l
            hessF_sparse_num[lp_l, lp_q] = Symbolics.sparsehessian(defn.F[lp_l, lp_q], s_x)
        end
    end

    println("hessf: ", hessf_sparse_num)
    println("hessce: ", hessce_sparse_num)
    println("hessci: ", hessci_sparse_num)
    println("hessF: ", hessF_sparse_num)

    # dp entries
    fdp_sparse_num = Vector{Num}(undef, r)
    for lp_r=1:r
        fdp_sparse_num[lp_r] = Symbolics.derivative(defn.f, s_pr[lp_r]; simplify=true)
    end
    cedp_sparse_num = Vector{Vector{Num}}(undef, r)
    for lp_r=1:r
        cedp_sparse_num[lp_r] = Symbolics.derivative(defn.ce, s_pr[lp_r]; simplify=true)
    end
    cidp_sparse_num = Vector{Vector{Num}}(undef, r)
    for lp_r=1:r
        cidp_sparse_num[lp_r] = Symbolics.derivative(defn.ci, s_pr[lp_r]; simplify=true)
    end
    Fdp_sparse_num = Vector{Matrix{Num}}(undef, r)
    for lp_r=1:r
        Fdp_sparse_num[lp_r] = Symbolics.derivative(defn.F, s_pr[lp_r]; simplify=true)
    end

    println("fdp: ", fdp_sparse_num)
    println("cedp: ", cedp_sparse_num)
    println("cidp: ", cidp_sparse_num)
    println("Fdp: ", Fdp_sparse_num)

    # graddp entries
    gradfdp_sparse_num = Vector{SparseMatrixCSC}(undef, r)
    for lp_r=1:r
        gradfdp_sparse_num[lp_r] = transpose( Symbolics.sparsejacobian([fdp_sparse_num[lp_r]], s_x) )
    end
    gradcedp_sparse_num = Vector{SparseMatrixCSC}(undef, r)
    for lp_r=1:r
        gradcedp_sparse_num[lp_r] = transpose( Symbolics.sparsejacobian( cedp_sparse_num[lp_r], s_x ) )
    end
    gradcidp_sparse_num = Vector{SparseMatrixCSC}(undef, r)
    for lp_r=1:r
        gradcidp_sparse_num[lp_r] = transpose( Symbolics.sparsejacobian( cidp_sparse_num[lp_r], s_x ) )
    end
    gradFdp_sparse_num = Vector{Matrix{SparseMatrixCSC}}(undef, r)
    for lp_r=1:r
        gradFdp_sparse_num[lp_r] = Matrix{SparseMatrixCSC}(undef, l, q)
        for lp_q=1:q
            for lp_l=1:l
                gradFdp_sparse_num[lp_r][lp_l, lp_q] = transpose( Symbolics.sparsejacobian( [ Fdp_sparse_num[lp_r][lp_l, lp_q] ], s_x ) )
            end
        end
    end

    println("gradfdp: ", gradfdp_sparse_num)
    println("gradcedp: ", gradcedp_sparse_num)
    println("gradcidp: ", gradcidp_sparse_num)
    println("gradFdp: ", gradFdp_sparse_num)

    # f, ce, ci, F all return usual dense vectors/matrices because they are, well, dense.

    # Full f functions (really just converting a 1 element vector to a scalar)
    function local_f(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}) where {S <: Real, T <: Real}
        return fns.f(x, pr, ps)[1]
    end

    # Full ce functions
    function local_ce(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}) where {S <: Real, T <: Real}   # ::Vector{S}
        return fns.ce(x, pr, ps)
    end

    # Indexed ce functions
    function local_ce(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}, idxs_ce::Vector{Int64}) where {S <: Real, T <: Real} # ::Vector{S}
        return mm_fd_dn_ce_i(dimspec, fns.ce_i, x, pr, ps, idxs_ce)
    end

    # Full ci functions
    function local_ci(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}) where {S <: Real, T <: Real}  # ::Vector{S}
        return fns.ci(x, pr, ps)
    end

    # Indexed ci functions
    function local_ci(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}, idxs_ci::Vector{Int64}) where {S <: Real, T <: Real}  # ::Vector{S}
        return mm_fd_dn_ci_i(dimspec, fns.ci_i, x, pr, ps, idxs_ci)
    end

    # Full F functions
    function local_F(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}) where {S <: Real, T <: Real}  # ::Matrix{S}
        return fns.F(x, pr, ps)
    end

    # Indexed F functions; NOTE, these yeild a vector in the order that the (l,q) index tuples were in!
    function local_F(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}, idxs_F::Vector{Tuple{Int64,Int64}}) where {S <: Real, T <: Real}  # ::Vector{S}
        return mm_fd_dn_F_i(dimspec, fns.F_i, x, pr, ps, idxs_F)
    end

    # gradf functions
    gradf_sparse_fn = build_function(gradf_sparse_num, x, pr, ps; expression=Val{false})

    function local_gradf(x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}) where {S <: Real, T <: Real}
        return gradf_sparse_fn[1](x0, pr0, ps0)
    end

    # Development placeholder functions
    # local_f(x0::Vector) = sum(x)
    # function local_f!(out_f::AbstractVector, x0::Vector)
    #     out_f[1] = sum(x)
    #     return nothing
    # end

    # local_ce(x0::Vector) = sum(x)
    # function local_ce!(out_ce::AbstractVector, x0::Vector)
    #     out_ce[1] = sum(x)
    #     return nothing
    # end

    # local_ci(x0::Vector) = sum(x)
    # function local_ci!(out_ci::AbstractVector, x0::Vector)
    #     out_ci[1] = sum(x)
    #     return nothing
    # end

    # local_F(x0::Vector) = sum(x)
    # function local_F!(out_F::AbstractVector, x0::Vector)
    #     out_F[1] = sum(x)
    #     return nothing
    # end

    # local_gradf(x0::Vector) = sum(x)
    # function local_gradf!(out_gradf::AbstractVector, x0::Vector)
    #     out_gradf[1] = sum(x)
    #     return nothing
    # end

    local_gradce(x0::Vector) = sum(x)
    function local_gradce!(out_gradce::AbstractVector, x0::Vector)
        out_gradce[1] = sum(x)
        return nothing
    end

    local_gradci(x0::Vector) = sum(x)
    function local_gradci!(out_gradci::AbstractVector, x0::Vector)
        out_gradci[1] = sum(x)
        return nothing
    end

    local_gradF(x0::Vector) = sum(x)
    function local_gradF!(out_gradF::AbstractVector, x0::Vector)
        out_gradF[1] = sum(x)
        return nothing
    end

    local_hessf(x0::Vector) = sum(x)
    function local_hessf!(out_hessf::AbstractVector, x0::Vector)
        out_hessf[1] = sum(x)
        return nothing
    end

    local_hessce(x0::Vector) = sum(x)
    function local_hessce!(out_hessce::AbstractVector, x0::Vector)
        out_hessce[1] = sum(x)
        return nothing
    end

    local_hessci(x0::Vector) = sum(x)
    function local_hessci!(out_hessci::AbstractVector, x0::Vector)
        out_hessci[1] = sum(x)
        return nothing
    end

    local_hessF(x0::Vector) = sum(x)
    function local_hessF!(out_hessF::AbstractVector, x0::Vector)
        out_hessF[1] = sum(x)
        return nothing
    end

    local_fdp(x0::Vector) = sum(x)
    function local_fdp!(out_fdp::AbstractVector, x0::Vector)
        out_fdp[1] = sum(x)
        return nothing
    end

    local_cedp(x0::Vector) = sum(x)
    function local_cedp!(out_cedp::AbstractVector, x0::Vector)
        out_cedp[1] = sum(x)
        return nothing
    end

    local_cidp(x0::Vector) = sum(x)
    function local_cidp!(out_cidp::AbstractVector, x0::Vector)
        out_cidp[1] = sum(x)
        return nothing
    end

    local_Fdp(x0::Vector) = sum(x)
    function local_Fdp!(out_Fdp::AbstractVector, x0::Vector)
        out_Fdp[1] = sum(x)
        return nothing
    end

    local_gradfdp(x0::Vector) = sum(x)
    function local_gradfdp!(out_gradfdp::AbstractVector, x0::Vector)
        out_gradfdp[1] = sum(x)
        return nothing
    end

    local_gradcedp(x0::Vector) = sum(x)
    function local_gradcedp!(out_gradcedp::AbstractVector, x0::Vector)
        out_gradcedp[1] = sum(x)
        return nothing
    end

    local_gradcidp(x0::Vector) = sum(x)
    function local_gradcidp!(out_gradcidp::AbstractVector, x0::Vector)
        out_gradcidp[1] = sum(x)
        return nothing
    end

    local_gradFdp(x0::Vector) = sum(x)
    function local_gradFdp!(out_gradFdp::AbstractVector, x0::Vector)
        out_gradFdp[1] = sum(x)
        return nothing
    end

    # TODO need to use the [1] format
    return MPCCModel(
            config,
            local_f, missing,
            local_ce, missing,
            local_ci, missing,
            local_F, missing,
            local_gradf, missing,
            local_gradce, missing,
            local_gradci, missing,
            local_gradF, missing,
            local_hessf, missing,
            local_hessce, missing,
            local_hessci, missing,
            local_hessF, missing,
            local_fdp, missing,
            local_cedp, missing,
            local_cidp, missing,
            local_Fdp, missing,
            local_gradfdp, missing,
            local_gradcedp, missing,
            local_gradcidp, missing,
            local_gradFdp, missing
        )
end


text:
using BayesianLinearRegression

import Random
import Printf: @printf

Random.seed!(2)
n = 20;
xx = range(0, 2π, length=n);
y = sin.(xx) .+ 0.1 .* randn(n);

function compareModels(p)
    X = chebyshev(xx,p)
    m = BayesianLinReg(X,y)
    fit!(m, callback=stopAtIteration(10))
    return logEvidence(m)
end

for (p,logEv) in enumerate(compareModels.(1:20))
    @printf "p = %2d: evidence = %8.2f\n" p exp(logEv)
end

X = chebyshev(xx,6);
@time m = BayesianLinReg(X,y) |> fit!;

priorScale(m)

noiseScale(m)

posteriorWeights(m)

effectiveNumParameters(m)

yhat = predict(m,m.X)


text:
"""
    iht(filename, k, d, phenotypes=6, covariates="", summaryfile="iht.summary.txt",
        betafile="iht.beta.txt", kwargs...)

Runs IHT with sparsity level `k`.

# Arguments
- `filename`: A `String` for VCF, binary PLINK, or BGEN file. VCF files should end
    in `.vcf` or `.vcf.gz`. Binary PLINK files should exclude `.bim/.bed/.fam`
    trailings but the trio must all be present in the same directory. BGEN files
    should end in `.bgen`.
- `k`: An `Int` for sparsity parameter = number of none-zero coefficients
- `d`: Distribution of phenotypes. Specify `Normal` for quantitative traits,
    `Bernoulli` for binary traits, `Poisson` or `NegativeBinomial` for
    count traits, and `MvNormal` for multiple quantitative traits.

# Optional Arguments
- `phenotypes`: Phenotype file name (`String`), an integer, or vector of integer. Integer(s)
    coresponds to the column(s) of PLINK's `.fam` file that stores phenotypes (default `phenotypes=6`).
    Enter multiple integers for multivariate analysis (e.g. `phenotypes=[6, 7]`).
    We recognize missing phenotypes as `NA` or `-9`. For quantitative traits
    (univariate or multivariate), missing phenotypes are imputed with the mean. Binary
    and count phenotypes cannot be imputed. Phenotype files are read using `readdlm` function
    in Julia base. We require each subject's phenotype to occupy a different row, and for
    multivariate analysis, each phenotype is comma separated. The phenotype file
    should not include a header line. Each row should be listed in the same order as in
    the PLINK and (for multivariate analysis) be comma separated.
- `covariates`: Covariate file name. Default `covariates=""` (in which case an intercept
    term will be automatically included). If `covariates` file specified, it will be
    read using `readdlm` function in Julia base. We require the covariate file to be
    comma separated, and not include a header line. Each row should be listed in the
    same order as in the PLINK. The first column should be all 1s to indicate an
    intercept. All other columns not specified in `exclude_std_idx` will be standardized
    to mean 0 variance 1.
- `summaryfile`: Output file name for saving IHT's summary statistics. Default
    `summaryfile="iht.summary.txt"`.
- `betafile`: Output file name for saving IHT's estimated genotype effect sizes.
    Default `betafile="iht.beta.txt"`.
- `covariancefile`: Output file name for saving IHT's estimated trait covariance
    matrix for multivariate analysis. Default `covariancefile="iht.cov.txt"`.
- `exclude_std_idx`: Indices of non-genetic covariates that should be excluded from
    standardization.
- `dosage`: Currently only guaranteed to work for VCF files. If `true`, will read
    genotypes dosages (i.e. `X[i, j] ∈ [0, 2]` before standardizing)
- All optional arguments available in [`fit_iht`](@ref)
"""
function iht(
    filename::AbstractString,
    k::Int,
    d::UnionAll;
    phenotypes::Union{AbstractString, Int, AbstractVector{Int}} = 6,
    covariates::AbstractString = "",
    summaryfile::AbstractString = "iht.summary.txt",
    betafile::AbstractString = "iht.beta.txt",
    covariancefile::AbstractString = "iht.cov.txt",
    exclude_std_idx::AbstractVector{<:Integer} = Int[],
    dosage::Bool = false,
    kwargs...
    )
    # read genotypes
    X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt = parse_genotypes(filename, dosage)
    if typeof(X) <: SnpData
        xla = SnpLinAlg{Float64}(X.snparray, model=ADDITIVE_MODEL,
            center=true, scale=true, impute=true)
    else
        xla = X # numeric matrix from VCF or BGEN files
    end

    # read phenotypes
    y = parse_phenotypes(X, phenotypes, d())

    # read and standardize covariates
    z = covariates == "" ? ones(size(xla, 1)) :
        parse_covariates(covariates, exclude_std_idx, standardize=true)
    is_multivariate(y) && (z = convert(Matrix{Float64}, Transpose(z)))

    # run IHT
    io = open(summaryfile, "w")
    if is_multivariate(y)
        result = fit_iht(y, Transpose(xla), z, k=k, io=io; kwargs...)
    else
        l = d == NegativeBinomial ? LogLink() : canonicallink(d()) # link function
        result = fit_iht(y, xla, z, k=k, d=d(), l=l, io=io; kwargs...)
    end

    show(io, result)

    if is_multivariate(y)
        open(betafile, "w") do io
            print(io, "SNPid")
            for i in 1:size(y, 1)
                print(io, '\t', "beta_$i")
            end
            print(io, '\n')
            writedlm(io, [X_ids result.beta'])
        end
        # writedlm(betafile, DataFrame([X_ids result.beta'], ["SNPid", "Estimated_beta"]))
        writedlm(covariancefile, result.Σ)
    else
        open(betafile, "w") do io
            println(io, "SNPid", '\t', "Estimated_beta")
            writedlm(io, [X_ids result.beta])
        end
    end

    close(io)
    flush(io)

    return result
end

# adhoc constructor for empty MvNormal distribution
MvNormal() = MvNormal(Float64[])

"""
    parse_phenotypes(x, col::Union{Int, AbstractVector{Int}}, ::Distribution)

Reads phenotypes to numeric array. If `x` is a `SnpData`, columns `col` of the `.fam`
file will be parsed as phenotypes. Otherwise, will read `x` as comma-separated text
file where each sample occupies a row. We recognize missing phenotypes as `NA` or
`-9`. For quantitative traits (univariate or multivariate), missing phenotypes are
imputed with the mean. Binary and count phenotypes cannot be imputed.
"""
function parse_phenotypes end

function parse_phenotypes(x::SnpData, col::AbstractVector{Int}, ::MvNormal)
    n = x.people
    r = length(col) # number of traits
    y = Matrix{Float64}(undef, r, n)
    offset = 5

    # impute missing phenotypes "-9" by mean of observed phenotypes
    missing_idx = Int[]
    for c in col
        empty!(missing_idx)
        s = 0.0
        for i in 1:n
            if phenotype_is_missing(x.person_info[i, c])
                y[c - offset, i] = 0.0
                push!(missing_idx, i)
            else
                y[c - offset, i] = parse(Float64, x.person_info[i, c])
                s += y[c - offset, i]
            end
        end
        avg = s / (n - length(missing_idx))
        for i in missing_idx
            y[c - offset, i] = avg
        end
    end
    return y
end

parse_phenotypes(::SnpData, ::Int, ::MvNormal) =
    throw(ArgumentError("Multivariate analysis requires multiple phenotypes! Please specify " *
        "e.g. phenotypes=[6, 7] or save each sample's phenotypes in a comma-" *
        "separated file where each sample occupies a different row and each" *
        " phenotype is separated by a single comma."))

function parse_phenotypes(x::SnpData, col::Int, ::Normal)
    n = x.people
    y = Vector{Float64}(undef, n)

    # impute missing phenotypes by mean of observed phenotypes
    missing_idx = Int[]
    s = 0.0
    for i in 1:n
        if phenotype_is_missing(x.person_info[i, col])
            y[i] = 0.0
            push!(missing_idx, i)
        else
            y[i] = parse(Float64, x.person_info[i, col])
            s += y[i]
        end
    end
    avg = s / (n - length(missing_idx))
    for i in missing_idx
        y[i] = avg
    end
    return y
end

function parse_phenotypes(x::SnpData, col::Int, ::UnivariateDistribution)
    n = x.people
    y = Vector{Float64}(undef, n)

    # missing phenotypes NOT allowed for binary/count phenotypes
    for i in 1:n
        if phenotype_is_missing(x.person_info[i, col])
            throw(MissingException("Missing phenotype detected for sample $i. Automatic " *
                "phenotype imputation are only possible for quantitative traits. " *
                "Please exclude missing phenotypes or impute them first."))
        else
            y[i] = parse(Float64, x.person_info[i, col])
        end
    end
    return y
end

function parse_phenotypes(::Any, pheno_filename::AbstractString, d)
    y = readdlm(pheno_filename, ',', Float64)
    if is_multivariate(y)
        y = convert(Matrix{Float64}, Transpose(y))
    else
        y = dropdims(y, dims=2)
    end
    return y
end

"""
    parse_covariates(filename, exclude_std_idx; standardize::Bool=true)

Reads a comma separated text file `filename`. Each row should be a sample ordered the
same as in the plink file. The first column should be array of 1 (representing
intercept). Each covariate should be comma separated. If `standardize=true`,
all columns except those in `exclude_std_idx` will be standardized.
"""
function parse_covariates(filename::AbstractString, exclude_std_idx::AbstractVector{<:Integer};
    standardize::Bool=true)
    z = readdlm(filename, ',', Float64)

    if eltype(exclude_std_idx) == Bool
        std_idx = .!exclude_std_idx
    else
        std_idx = trues(size(z, 2))
        std_idx[exclude_std_idx] .= false
    end

    if all(x == 1 for x in @view(z[:, 1]))
        std_idx[1] = false # don't standardize intercept
    else
        @warn("Covariate file provided but did not detect an intercept. An intercept will NOT be included in IHT!")
    end

    standardize && standardize!(@view(z[:, std_idx]))
    return z
end

function phenotype_is_missing(s::AbstractString)
    return s == "-9" || s == "NA"
end

"""
    cross_validate(filename, d, path=1:20, phenotypes=6, covariates="",
        cv_summaryfile="cviht.summary.txt", q=5, kwargs...)

Runs cross-validation to determinal optimal sparsity level `k`. Different
sparsity levels are specified in `path`.

# Arguments
- `filename`: A `String` for input PLINK file name (without `.bim/.bed/.fam` suffixes)
- `d`: Distribution of phenotypes. Specify `Normal` for quantitative traits,
    `Bernoulli` for binary traits, `Poisson` or `NegativeBinomial` for
    count traits, and `MvNormal` for multiple quantitative traits.

# Optional Arguments
- `path`: Different values of `k` that should be tested. One can input a vector of
    `Int` (e.g. `path=[5, 10, 15, 20]`) or a range (default `path=1:20`).
- `phenotypes`: Phenotype file name (`String`), an integer, or vector of integer. Integer(s)
    coresponds to the column(s) of `.fam` file that stores phenotypes (default 6).
    We recognize missing phenotypes as `NA` or `-9`. For quantitative traits
    (univariate or multivariate), missing phenotypes are imputed with the mean. Binary
    and count phenotypes cannot be imputed. Phenotype files are read using `readdlm` function
    in Julia base. We require each subject's phenotype to occupy a different row, and for
    multivariate analysis, each phenotype is comma separated. The phenotype file
    should not include a header line. Each row should be listed in the same order as in
    the PLINK.
- `covariates`: Covariate file name. Default `covariates=""` (in which case an intercept
    term will be automatically included). If `covariates` file specified, it will be
    read using `readdlm` function in Julia base. We require the covariate file to be
    comma separated, and not include a header line. Each row should be listed in the
    same order as in the PLINK. The first column should be all 1s to indicate an
    intercept. All other columns not specified in `exclude_std_idx` will be standardized
    to mean 0 variance 1
- `cv_summaryfile`: Output file name for saving IHT's cross validation summary statistics.
    Default `cv_summaryfile="cviht.summary.txt"`.
- `q`: Number of cross validation folds. Larger means more accurate and more computationally
    intensive. Should be larger 2 and smaller than 10. Default `q=5`.
- `dosage`: Currently only guaranteed to work for VCF files. If `true`, will read
    genotypes dosages (i.e. `X[i, j] ∈ [0, 2]` before standardizing)
- All optional arguments available in [`cv_iht`](@ref)
"""
function cross_validate(
    filename::AbstractString,
    d::UnionAll;
    path::AbstractVector{<:Integer} = 1:20,
    phenotypes::Union{AbstractString, Int, AbstractVector{Int}} = 6,
    covariates::AbstractString = "",
    cv_summaryfile::AbstractString = "cviht.summary.txt",
    q::Int = 5,
    exclude_std_idx::AbstractVector{<:Integer} = Int[],
    dosage::Bool = false,
    kwargs...
    )
    start_time = time()

    # read genotypes
    X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt = parse_genotypes(filename, dosage)
    if typeof(X) <: SnpData
        x = SnpLinAlg{Float64}(X.snparray, model=ADDITIVE_MODEL,
            center=true, scale=true, impute=true)
    else
        x = X # numeric matrix from VCF or BGEN files
    end

    # read phenotypes
    y = parse_phenotypes(X, phenotypes, d())

    # read and standardize covariates
    z = covariates == "" ? ones(size(x, 1)) :
        parse_covariates(covariates, exclude_std_idx, standardize=true)
    is_multivariate(y) && (z = convert(Matrix{Float64}, Transpose(z)))

    # run cross validation
    if is_multivariate(y)
        mse = cv_iht(y, Transpose(x), z, path=path, q=q; kwargs...)
    else
        l = d == NegativeBinomial ? LogLink() : canonicallink(d()) # link function
        mse = cv_iht(y, x, z, path=path, q=q, d=d(), l=l; kwargs...)
    end

    # save results
    open(cv_summaryfile, "w") do io
        k = path[argmin(mse)]
        print_cv_results(io, mse, path, k)
        end_time = time() - start_time
        println(io, "Total cross validation time = $end_time seconds")
    end

    return mse
end

"""
    convert_bgen_gt(t::Type{T}, b::Bgen)

Imports BGEN genotypes and chr/sampleID/pos/snpID/ref/alt into numeric arrays.
Genotypes are centered and scaled to mean 0 variance 1. Missing genotypes will
be replaced with the mean. Assumes every variant is biallelic (ie only 1 alt allele).

# Input
- `b`: a `Bgen` object
- `T`: Type for genotype array

# Output
- `G`: matrix of genotypes with type `T`.
"""
function convert_bgen_gt(t::Type{T}, b::Bgen) where T <: Real
    n = n_samples(b)
    p = n_variants(b)

    # return arrays
    G = Matrix{t}(undef, n, p)
    Gchr = Vector{String}(undef, p)
    Gpos = Vector{Int}(undef, p)
    GsnpID = [String[] for _ in 1:p] # each variant can have >1 rsid, although we don't presently allow this
    Gref = Vector{String}(undef, p)
    Galt = [String[] for _ in 1:p] # each variant can have >1 alt allele, although we don't presently allow this

    # import each variant
    i = 1
    for v in iterator(b; from_bgen_start=true)
        dose = ref_allele_dosage!(b, v; T=t) # this reads REF allele as 1
        BGEN.alt_dosage!(dose, v.genotypes.preamble) # switch 2 and 0 (ie treat ALT as 1)
        copyto!(@view(G[:, i]), dose)
        # store chr/pos/snpID/ref/alt info
        Gchr[i], Gpos[i] = chrom(v), pos(v)
        push!(GsnpID[i], rsid(v))
        ref_alt_alleles = alleles(v)
        length(ref_alt_alleles) > 2 && error("Marker $i of BGEN is not biallelic!")
        Gref[i] = ref_alt_alleles[1]
        push!(Galt[i], ref_alt_alleles[2])
        i += 1
        clear!(v)
    end

    # center/scale/impute
    standardize_genotypes!(G)

    return G, b.samples, Gchr, Gpos, GsnpID, Gref, Galt
end

"""
    standardize_genotypes!(G::AbstractMatrix)

Centers and scales each column (SNP) of `G` to mean 0 variance 1. Also each
missing entry will be imputed as mean.
"""
function standardize_genotypes!(G::AbstractMatrix)
    T = eltype(G)
    @inbounds for snp in eachcol(G)
        μi, mi = zero(T), 0
        @simd for i in eachindex(snp)
            μi += isnan(snp[i]) ? zero(t) : snp[i]
            mi += isnan(snp[i]) ? 0 : 1
        end
        μi /= mi
        σi = sqrt(μi * (1 - μi / 2))
        @simd for i in eachindex(snp)
            isnan(snp[i]) && (snp[i] = μi) # impute
            snp[i] -= μi # center
            σi > 0 && (snp[i] /= σi) # scale
        end
    end
    return nothing
end

"""
    parse_genotypes(tgtfile::AbstractString, dosage=false)

Imports genotype data from `tgtfile`. If binary PLINK files are supplied, genotypes
will not be decompressed to numeric matrices (~2 bit per entry). VCF or BGEN genotypes
will be stored in double precision matrices (64 bit per entry).

# Inputs
- `tgtfile`: VCF, binary PLINK, or BGEN file. VCF files should end in `.vcf` or
    `.vcf.gz`. Binary PLINK files should exclude `.bim/.bed/.fam` trailings but
    the trio must all be present in the same directory. BGEN files should end in
    `.bgen`

# Optional Inputs
- `dosage`: Currently only guaranteed to work for VCF files. If `true`, will read
    genotypes dosages (i.e. `X[i, j] ∈ [0, 2]` before standardizing)

# Output
- `X`: a `n × p` genotype matrix of type `Float64` (VCF or BGEN inputs) or
    `SnpArrays.SnpData` (binary PLINK inputs)
- `Gchr`: Vector of `String`s holding chromosome number for each variant
- `Gpos`: Vector of `Int` holding each variant's position
- `GsnpID`: Vector of `String`s holding variant ID for each variant
- `Gref`: Vector of `String`s holding reference allele for each variant
- `Galt`: Vector of `String`s holding alterante allele for each variant
"""
function parse_genotypes(tgtfile::AbstractString, dosage=false)
    if (endswith(tgtfile, ".vcf") || endswith(tgtfile, ".vcf.gz"))
        f = dosage ? convert_ds : convert_gt
        X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt =
            f(Float64, tgtfile, trans=false,
            save_snp_info=true, msg = "Importing from VCF file...")
        # convert missing to NaN
        replace!(X, missing => NaN)
        X = convert(Matrix{Float64}, X) # drop Missing from Matrix type
        # center/scale/impute
        standardize_genotypes!(X)
    elseif endswith(tgtfile, ".bgen")
        # dosage && error("Currently don't support reading dosage data from BGEN files")
        samplefile = isfile(tgtfile[1:end-5] * ".sample") ?
            tgtfile[1:end-5] * ".sample" : nothing
        indexfile = isfile(tgtfile * ".bgi") ? tgtfile * ".bgi" : nothing
        bgen = Bgen(tgtfile; sample_path=samplefile, idx_path=indexfile)
        X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt = MendelIHT.convert_bgen_gt(Float64, bgen)
    elseif isplink(tgtfile)
        dosage && error("PLINK files detected but dosage = true!")
        X = SnpArrays.SnpData(tgtfile)
        # get other relevant information
        X_sampleID = X.person_info[!, :iid]
        X_chr = X.snp_info[!, :chromosome]
        X_pos = X.snp_info[!, :position]
        X_ids = X.snp_info[!, :snpid]
        X_ref = X.snp_info[!, :allele1]
        X_alt = X.snp_info[!, :allele2]
    else
        error("Unrecognized target file format: target file can only be VCF" *
            " files (ends in .vcf or .vcf.gz), BGEN (ends in .bgen) or PLINK" *
            " (do not include.bim/bed/fam) and all trio must exist in 1 directory)")
    end
    return X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt
end

isplink(tgtfile::AbstractString) = isfile(tgtfile * ".bed") &&
                                   isfile(tgtfile * ".fam") &&
                                   isfile(tgtfile * ".bim")


text:
function orthogonalize(x, bases)
    F = 1.0
    x_new = copy(x)
    @inbounds for j = 1:length(bases)
        cji = bases[j] ⋅ x
        F -= cji ^ 2
        x_new -= cji * bases[j]
    end
    x_new /= √F
    x_new
end

@doc raw"""All Constraints. Fixed lengths, potential energy, etc.

```math
P = 1 - \sum_{i}^{N_C}|C_i'\rangle\langle C_i'|
```
"""
function projection_matrix(∇Cs)
    N_C = length(∇Cs)
    ∇C_unit_new = normalize.(∇Cs)
    P = I
    @inbounds for i = 1:N_C
        ∇C_unit_new[i] = orthogonalize(∇C_unit_new[i], ∇C_unit_new[1:i-1])
        tmp = flatten(∇C_unit_new[i])
        P -= tmp * tmp'
    end
    P
end


text:
using Requests
using MatrixDepot
using DataFrames
using DataFramesMeta
using PETScBinaryIO
using MatrixMarket
using MAT

function max_by(fn, itr)
  v = fn(itr[1])
  e = itr[1]
  for i in itr
    vv = fn(i)
    if vv > v
      v = vv
      e = i
    end
  end

  e
end

function largest_component(L)
  cc = connected_components(L)
  lc = sort(max_by(length, cc))
  L[lc,lc]
end

function connected_components(L)
  components = Vector{Vector{Int64}}()
  notfound = Set{Int64}(1:size(L, 1))
  while !isempty(notfound)
    s = first(notfound)
    old = Set{Int64}()
    frontier = Set{Int64}(s)
    new = Set{Int64}(s)
    while !isempty(new)
      union!(old, new)
      frontier = copy(new)
      empty!(new)
      for v in frontier
        for n in findn(L[:, v])
          if !(n in old || n in new || n in frontier)
            push!(new, n)
          end
        end
      end
    end
    push!(components, collect(old))
    setdiff!(notfound, old)
  end

  components
end

"""
Get metadata information about all graphs in the UF sparse matrix collection
"""
function get_ufcollection_metadata()
    csv = String(get("https://www.cise.ufl.edu/research/sparse/matrices/UFstats.csv").data)
    i = findnext(x -> x == '\n', csv, 1)
    i = findnext(x -> x == '\n', csv, i+1)
    df = readtable(IOBuffer(csv[i+1:end]), header=false, names = [:group, :name, :rows, :cols, :nnz, :real, :binary, :geom, :posdef, :psym, :nsym, :kind])
    df[:geom] = df[:geom] .== 1
    df[:posdef] = df[:posdef] .== 1
    df[:binary] = df[:binary] .== 1
    df[:real] = df[:real] .== 1
    df[:id] = 1:nrow(df)
    df
end

graph_kinds = [ "directed weighted graph"
              , "directed graph"
              , "undirected weighted graph"
              , "undirected random graph"
              , "undirected weighted random graph"
              , "duplicate undirected random graph"
              , "directed weighted random graph"
              , "weighted bipartite graph"
              , "undirected graph"
              , "bipartite graph"
              , "directed multigraph"
              , "undirected multigraph"
              , "undirected weighted graph sequence"
              , "undirected bipartite graph"
              , "undirected graph sequence"
              , "directed weighted graph sequence"
              , "weighted undirected graph"
              , "weighted directed graph"
              , "term/document graph"
              , "random undirected graph"
              , "random unweighted graph"
              ]

"""
Filter `df` by `:kind in kinds`
"""
function select_by_kinds(df, kinds)
    mapreduce(k -> @where(df, :kind .== k), vcat, kinds)
end

function download_graphs(filepath, df, formats)
    for (g, n) in zip(df[:group], df[:name])
        if isfile("$filepath/$n.mtx")
            continue
        end
        localpath = "$filepath/$n.tmp.mat"
        println("Downloading $g/$n to $localpath")
        println("https://www.cise.ufl.edu/research/sparse/mat/$g/$n.mat")
        download("https://www.cise.ufl.edu/research/sparse/mat/$g/$n.mat", localpath)
        A = matread(localpath)["Problem"]["A"]
        rm(localpath)
        if size(A, 1) != size(A, 2)
            A = A * A.'
            A = A - spdiagm(diag(A))
        end
        if !issymmetric(A)
            A = (A + A.') ./ 2
        end
        d = sum(A, 2)
        L = spdiagm(vec(d)) - A
        L = largest_component(L)
        for f in formats
            if f == "mat"
                matwrite("$filepath/$n.mat", Dict("L"=>L))
            elseif f == "mtx"
                mmwrite("$filepath/$n.mtx", L)
            elseif f == "petsc"
                writePETSc("$filepath/$n.petsc", L)
            end
        end
    end
end

function download_all_graphs(filepath)
    df = get_ufcollection_metadata()
    gs = select_by_kinds(df, graph_kinds)
    download_graphs(gs, ["mtx"])
end

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
using Test
unstable(x) = x <= 0 ? 0 : x

@inferred unstable(1.0)
# @code_warntype unstable(1.0)


text:
# ---
# title: 209. Minimum Size Subarray Sum
# id: problem209
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Two Pointers, Binary Search
# link: <https://leetcode.com/problems/minimum-size-subarray-sum/description/>
# hidden: true
# ---
#
# Given an array of **n** positive integers and a positive integer **s** , find
# the minimal length of a **contiguous** subarray of which the sum ≥ **s**. If
# there isn't one, return 0 instead.
#
# **Example:  **
#
#
#
#     Input: s = 7, nums = [2,3,1,2,4,3]
#     Output: 2
#     Explanation: the subarray [4,3] has the minimal length under the problem constraint.
#
# **Follow up:**
#
# If you have figured out the _O_ ( _n_ ) solution, try coding another solution
# of which the time complexity is _O_ ( _n_ log _n_ ).
#
#
## @lc code=start
using LeetCode

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


text:
using Test
using ParameterTests
using Dates: Date

@testset "Parameter Tests" verbose=true begin
    @testset "Strategy Generators" verbose=true begin
        for type in (Int8, Int16, Int32, Int64, Int128)
            @paramtest "Integers - $(type)" begin
                @given a ∈ integers(type)
                @given b ∈ integers(type; min=0)
                @given c ∈ integers(type; min=-10, max=10)
                @given d ∈ integers(type; max=3)
                @test isa(a, type)
                @test isa(b, type) && b >= 0
                @test isa(c, type) && -10 <= c <= 10
                @test isa(d, type) && d <= 3
            end
        end

        @paramtest "Integers" begin
            @given a ∈ integers()
            @given b ∈ integers(; min=0)
            @given c ∈ integers(; min=-10, max=10)
            @given d ∈ integers(; max=3)
            @test isa(a, Int)
            @test isa(b, Int) && b >= 0
            @test isa(c, Int) && -10 <= c <= 10
            @test isa(d, Int) && d <= 3
        end

        @paramtest "Floats" begin
            @given a ∈ floats(Float16)
            @given b ∈ floats(Float16; min=0, max=1)
            @given c ∈ floats(; allow_nan=true, allow_inf=true)
            @given d ∈ floats(Float32; allow_nan=false, allow_inf=true)
            @test isa(a, Float16) && !isnan(a) && !isinf(a)
            @test isa(b, Float16) && 0 <= b <= 1
            @test isa(c, Float64)
            @test isa(d, Float32)
        end

        @paramtest "Dates" begin
            @given a ∈ dates()
            @given b ∈ dates(;min=Date(2020,1,1))
            @given c ∈ dates(;min=Date(2020,1,1), max=Date(2021,1,1))
            @given d ∈ dates(;max=Date(2020,1,1))
            @test isa(a, Date)
            @test isa(b, Date) && b >= Date(2020,1,1)
            @test isa(c, Date) && Date(2020,1,1) <= c <= Date(2021,1,1)
            @test isa(d, Date) && d <= Date(2020,1,1)
        end

        @paramtest "Strings" begin
            @given a ∈ strings(; minsize=0, maxsize=10)
            @given b ∈ strings(:alphanum; minsize=3000)
            @given c ∈ strings(['A','B']; maxsize=3)
            @test length(a) <= 10
            @test 3000 <= length(b) <= 5000
            @test length(c) <= 3 && Set(c) ⊆ ('A','B')
        end

        @paramtest "Symbols" begin
            @given a ∈ symbols(; minsize=0, maxsize=10)
            @given b ∈ symbols(:alphanum; minsize=3000)
            @given c ∈ symbols(['A','B']; maxsize=3)
            @test length(String(a)) <= 10 && isa(a, Symbol)
            @test 3000 <= length(String(b)) <= 5000
            @test length(String(c)) <= 3 && Set(String(c)) ⊆ ('A','B')
        end
    end

    @testset "Strategy Parameters" verbose=true begin
        # Simple demonstration that our sample testing works as expected.
        count = 0
        @paramtest "Samples" samples=500 begin
            @given a ∈ integers(), b ∈ integers()
            count += 1
            @test a + b == b + a
        end

        @test count == 500
    end

    @testset "Array Strategies" verbose=true begin
        @paramtest "Vectors" begin
            @given a ∈ vectors(integers(); maxsize=5)
            @given b ∈ vectors(dates(;min=Date(2020,1,1)); minsize=5, maxsize=10)
            @given c ∈ vectors(floats(); minsize=10)
            @given d ∈ vectors(floats(Float32;min=0, max=10))
            @given e ∈ vectors(strings(); maxsize=3)
            @given f ∈ vectors(symbols(); maxsize=3)
            @test isa(a, Vector{Int}) && length(a) <= 5
            @test isa(b, Vector{Date}) && 5 <= length(b) <= 10
            @test all(map(x -> x >= Date(2020,1,1), b))
            @test isa(c, Vector{Float64}) && length(c) >= 10
            @test isa(d, Vector{Float32})
            @test isa(e, Vector{String})
            @test isa(f, Vector{Symbol})
        end
    end

    @testset "Parameters" verbose=true begin
        # This provides a simple example of how we can select from a parameter list
        # rather than using sample strategies
        dataset = 1:5
        @paramtest "Data" begin
            @given x ∈ dataset
            @test x <= 5
        end

        @paramtest "Tuple Data" begin
            @given (x,y) ∈ ((1,2),(3,4))
            @test x + y == y + x
        end
    end
end


text:
export CircularAperture,
       CircularAnnulus

"""
    CircularAperture(x, y, r)
    CircularAperture([x, y], r)

A circular aperture.

# Examples
```jldoctest
julia> ap = CircularAperture(0, 0, 10)
CircularAperture(0, 0, r=10)

```
"""
struct CircularAperture{T <: Number} <: AbstractAperture
    x::T
    y::T
    r::T
end

CircularAperture(center::AbstractVector, r) = CircularAperture(center..., r)
CircularAperture(x, y, r) = CircularAperture(promote(x, y, r)...)

function Base.show(io::IO, c::CircularAperture)
    print(io, "CircularAperture($(c.x), $(c.y), r=$(c.r))")
end

bbox(c::CircularAperture{<:Integer}) = (c.x - c.r, c.x + c.r, c.y - c.r, c.y + c.r)

function bbox(c::CircularAperture{<:AbstractFloat})
    xmin = floor(Int, c.x - c.r + 0.5)
    xmax = ceil(Int, c.x + c.r + 0.5)
    ymin = floor(Int, c.y - c.r + 0.5)
    ymax = ceil(Int, c.y + c.r + 0.5)
    return (xmin, xmax, ymin, ymax)
end

function mask(c::CircularAperture; method = :exact)
    bounds = edges(c)
    box = bbox(c)
    ny, nx = size(c)
    return circular_overlap(bounds..., nx, ny, c.r, method = method)
end

#######################################################

"""
    CircularAnnulus(x, y, r_in, r_out)
    CircularAnnulus([x, y], r_in, r_out)

A circular aperture.

# Examples
```jldoctest
julia> ap = CircularAnnulus(0, 0, 5, 10)
CircularAnnulus(0, 0, r_in=5, r_out=10)

```
"""
struct CircularAnnulus{T <: Number} <: AbstractAperture
    x::T
    y::T
    r_in::T
    r_out::T
end

CircularAnnulus(center::AbstractVector, r_in, r_out) = CircularAnnulus(center..., r_in, r_out)
CircularAnnulus(x, y, r_in, r_out) = CircularAnnulus(promote(x, y, r_in, r_out)...)

function Base.show(io::IO, c::CircularAnnulus)
    print(io, "CircularAnnulus($(c.x), $(c.y), r_in=$(c.r_in), r_out=$(c.r_out))")
end

bbox(c::CircularAnnulus{<:Integer}) = (c.x - c.r_out, c.x + c.r_out, c.y - c.r_out, c.y + c.r_out)

function bbox(c::CircularAnnulus{<:AbstractFloat})
    xmin = floor(Int, c.x - c.r_out + 0.5)
    xmax = ceil(Int, c.x + c.r_out + 0.5)
    ymin = floor(Int, c.y - c.r_out + 0.5)
    ymax = ceil(Int, c.y + c.r_out + 0.5)
    return (xmin, xmax, ymin, ymax)
end

function mask(c::CircularAnnulus; method = :exact)
    bounds = edges(c)
    box = bbox(c)
    ny, nx = size(c)
    out = circular_overlap(bounds..., nx, ny, c.r_out, method = method)
    out .-= circular_overlap(bounds..., nx, ny, c.r_in, method = method)
end

function edges(c::Union{CircularAperture,CircularAnnulus})
    ibox = bbox(c)
    xmin = ibox[1] - c.x - 0.5
    xmax = ibox[2] - c.x + 0.5
    ymin = ibox[3] - c.y - 0.5
    ymax = ibox[4] - c.y + 0.5
    return (xmin, xmax, ymin, ymax)
end


text:
using Sub2Index
using Test

@testset "Simple Indexing" begin

A = reshape(1:100, (10, 10))

@test @_ A₁ == A[1] == 1
@test @_ A₂₃ == A[23] == 23
@test @_ A₂₊₃ == A[2 + 3] == 5
@test @_ A₍₄₎₍₅₎ == A[4, 5] == 44

end

@testset "Fibonacci Sequence" begin

F = ones(Int, 10)

@_ for n = 3:10 Fₙ = Fₙ₋₁ + Fₙ₋₂ end

@test F == [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]

end

@testset "Pascal Triangle" begin

Π = zeros(Int, (10, 10))

@_ for i = 1:10 Π₍ᵢ₎₍₁₎ = Π₍ᵢ₎₍ᵢ₎ = 1 end
@_ for i = 3:10, j = 2:i-1
    Π₍ᵢ₎₍ⱼ₎ = Π₍ᵢ₋₁₎₍ⱼ₋₁₎ + Π₍ᵢ₋₁₎₍ⱼ₎
end

@test Π == [1  0   0   0    0    0   0   0  0  0
            1  1   0   0    0    0   0   0  0  0
            1  2   1   0    0    0   0   0  0  0
            1  3   3   1    0    0   0   0  0  0
            1  4   6   4    1    0   0   0  0  0
            1  5  10  10    5    1   0   0  0  0
            1  6  15  20   15    6   1   0  0  0
            1  7  21  35   35   21   7   1  0  0
            1  8  28  56   70   56  28   8  1  0
            1  9  36  84  126  126  84  36  9  1]

end

@testset "fallbacks" begin

@test @_(123) == 123
@test @_("Fₙ = Fₙ₋₁ + Fₙ₋₂") == "Fₙ = Fₙ₋₁ + Fₙ₋₂"
@test @_(:Fₙ) == :Fₙ

end


text:
module LanguageServer
using JSON
using REPL
import URIParser
import DocumentFormat
import CSTParser
import SymbolServer
using CSTParser
using CSTParser: EXPR
import CSTParser.Tokenize.Tokens

import StaticLint

export LanguageServerInstance

include("uri2.jl")
include("protocol/protocol.jl")
include("document.jl")
include("languageserverinstance.jl")
include("jsonrpc.jl")
include("staticlint.jl")

include("requests/init.jl")
include("requests/misc.jl")
include("requests/textdocument.jl")
include("requests/features.jl")
include("requests/hover.jl")
include("requests/completions.jl")
include("requests/workspace.jl")
include("parsing.jl")
include("utilities.jl")
include("jmd.jl")
include("display.jl")

end


text:
# Tests for what happens in the `utils.jl`

fs()
write(joinpath("_layout", "head.html"), "")
write(joinpath("_layout", "foot.html"), "")
write(joinpath("_layout", "page_foot.html"), "")
write("config.md", "")

@testset "utils1" begin
    write("utils.jl", """
        x = 5
        """)
    F.process_utils()
    s = """
        {{fill x}}
        """ |> fdi
    @test isapproxstr(s, "5")
    s = """
        {{x}}
        """ |> fdi
    @test isapproxstr(s, "5")
end

@testset "utils:hfun" begin
    write("utils.jl", """
        hfun_foo() = return "blah"
        """)
    F.process_utils()
    s = """
        {{foo}}
        """ |> fdi
    @test isapproxstr(s, "blah")

    write("utils.jl", """
        function hfun_bar(vname)
            val = locvar(vname[1])
            return round(sqrt(val), digits=2)
        end
        """)
    F.process_utils()
    s = """
        @def xx = 25
        {{bar xx}}
        """ |> fdi
    @test isapproxstr(s, "5.0")
end

@testset "utils:lxfun" begin
    write("utils.jl", """
        function lx_foo(lxc::Franklin.LxCom, _)
            return uppercase(Franklin.content(lxc.braces[1]))
        end
        """)
    F.process_utils()
    s = raw"""
        \foo{bar}
        """ |> fdi
    @test isapproxstr(s, "BAR")

    write("utils.jl", """
        function lx_baz(com, _)
            brace_content = Franklin.content(com.braces[1])
            return uppercase(brace_content)
        end
        """)
    F.process_utils()
    s = raw"""
        \baz{bar}
        """ |> fdi
    @test isapproxstr(s, "BAR")
end

# 452
@testset "utils:lxfun2" begin
    write("utils.jl", raw"""
        function lx_bold(com, _)
            text = Franklin.content(com.braces[1])
            return "**$text**"
        end
        """)
    F.process_utils()
    s = raw"""
        \bold{bar}
        """ |> fdi
    @test isapproxstr(s, "bar")
end