text:
using Combinatorics

@testset "Veronese Lift" begin
    @test_throws DimensionMismatch permanent([1 2], collect(Combinatorics.permutations(1:2)))
    @test_throws DimensionMismatch SwitchOnSafety.veroneselift_explicit([1 2], 3)
    @test_throws ArgumentError SwitchOnSafety.veroneselift_explicit([1 2; 3 4], 0)
    @test_throws ArgumentError veroneselift([1 2; 3 4], 0)

    @test [1.0 2 * √3 4 * √3 8.0] ≈ @inferred veroneselift([1 2], 3)

    s3 = sqrt(3)
    expected = [11^3 s3*11^2*12 s3*11*12^2 12^3;
    s3*11^2*21 11*(11*22+2*21*12) 12*(2*11*22+21*12) s3*12^2*22;
    s3*11*21^2 21*(2*11*22+21*12) 22*(11*22+2*21*12) sqrt(3)*12*22^2;
    21^3 s3*21^2*22 s3*21*22^2 22^3]
    @test expected ≈ @inferred veroneselift([11 12; 21 22], 3)
    @test expected ≈ @inferred SwitchOnSafety.veroneselift_explicit([11 12; 21 22], 3)
end


text:
function get_affine_map(ds::ArchGDAL.IDataset)
    # ArchGDAL fails hard on datasets without
    # an affinemap. GDAL documents that on fail
    # a default affinemap should be returned.
    local gt
    try
        gt = ArchGDAL.getgeotransform(ds)
    catch y
        @warn y.msg
        gt = [0.,1.,0.,0.,0.,1.]
    end
    geotransform_to_affine(gt)
end

function geotransform_to_affine(gt::SVector{6,<:AbstractFloat})
    # See https://lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html
    # for an explanation of the geotransform format
    AffineMap(SMatrix{2,2}([gt[2] gt[3]; gt[5] gt[6]]), SVector{2}([gt[1], gt[4]]))
end
geotransform_to_affine(A::Vector{<:AbstractFloat}) = geotransform_to_affine(SVector{6}(A))

function affine_to_geotransform(am::AffineMap)
    l = am.linear
    t = am.translation
    (length(l) == 4 && length(t) == 2) || error("AffineMap has wrong dimensions.")
    [t[1], l[1], l[3], t[2], l[2], l[4]]
end

"""Check wether the AffineMap of a GeoArray contains rotations."""
function is_rotated(ga::GeoArray)
    ga.f.linear[2] != 0. || ga.f.linear[3] != 0.
end

function bbox(ga::GeoArray)
    ax, ay = ga.f(SVector(0, 0))
    bx, by = ga.f(SVector(size(ga)[1:2]))
    (min_x = min(ax, bx), min_y = min(ay, by), max_x = max(ax, bx), max_y = max(ay, by))
end

function unitrange_to_affine(x::StepRangeLen, y::StepRangeLen)
    δx, δy = float(step(x)), float(step(y))
    AffineMap(
        SMatrix{2,2}(δx, 0., 0., δy),
        SVector(x[1] - δx / 2, y[1] - δy / 2)
    )
end

function bbox_to_affine(size::Tuple{Integer,Integer}, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
    AffineMap(
        SMatrix{2,2}(float(bbox.max_x - bbox.min_x) / size[1], 0, 0, float(bbox.max_y - bbox.min_y) / size[2]),
        SVector(float(bbox.min_x), float(bbox.min_y))
    )
end

"""Set geotransform of `GeoArray` by specifying a bounding box.
Note that this only can result in a non-rotated or skewed `GeoArray`."""
function bbox!(ga::GeoArray, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
    ga.f = bbox_to_affine(size(ga)[1:2], bbox)
    ga
end

"""Generate bounding boxes for GeoArray cells."""
function bboxes(ga::GeoArray)
    c = coords(ga, Vertex())
    m, n = size(c)
    cellbounds = Matrix{NamedTuple}(undef, (m - 1, n - 1))
    for j in 1:n - 1, i in 1:m - 1
        v = c[i:i + 1, j:j + 1]
        minx, maxx = extrema(first.(v))
        miny, maxy = extrema(last.(v))
        cellbounds[i, j] = (min_x = minx, max_x = maxx, min_y = miny, max_y = maxy)
    end
    cellbounds
end

# Extend CoordinateTransformations
CoordinateTransformations.compose(ga::GeoArray, t2::AffineMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::LinearMap) = CoordinateTransformations.compose(ga.f, t2)
CoordinateTransformations.compose(ga::GeoArray, t2::Translation) = CoordinateTransformations.compose(ga.f, t2)

"""Transform an GeoArray by applying a Transformation."""
function compose!(ga::GeoArray, t2::AffineMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::LinearMap)
    ga.f = compose(ga, t2)
end
function compose!(ga::GeoArray, t2::Translation)
    ga.f = compose(ga, t2)
end

"""Function to flip GeoArray upside down to adjust to GDAL ecosystem."""
function flipud!(ga::GeoArray)
    # Flip data upside down
    ga.A = reverse(ga.A, dims=2)

    # Find new corner coordinates
    ux, uy = ga.f(SVector{2}([0, size(ga)[2]]))

    # Define y mirror and compose
    lm = LinearMap(SMatrix{2,2}([1.0 0.0; 0.0 -1.0]))
    am = compose(ga.f, lm)  # AffineMap
    translate = SVector{2}([ux, uy])
    f = AffineMap(am.linear, translate)

    ga.f = f
    ga
end

"""Check bbox overlapping

Return true if two bboxes overlap.
"""
function bbox_overlap(
    bbox_a::NamedTuple{(:min_x, :min_y, :max_x, :max_y)},
    bbox_b::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})

    # Check if bboxes are valid
    if (bbox_a.min_x >= bbox_a.max_x) ||
        (bbox_a.min_y >= bbox_a.max_y)
        error("Invalid bbox (min >= max)", bbox_a)
    end
    if (bbox_b.min_x >= bbox_b.max_x) ||
        (bbox_b.min_y >= bbox_b.max_y)
        error("Invalid bbox (min >= max)", bbox_b)
    end

    # Check if bboxes overlap
    if (bbox_a.max_x < bbox_b.min_x) ||
        (bbox_a.max_y < bbox_b.min_y) ||
        (bbox_a.min_x > bbox_b.max_x) ||
        (bbox_a.min_y > bbox_b.max_y)
        return false
    else
        return true
    end
end

"""
    function crop(ga::GeoArray, cbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y),Tuple{Float64,Float64,Float64,Float64}})

Crop input GeoArray by coordinates (box or another GeoArray). If the coordinates range is larger than the GeoArray, only the overlapping part is given in the result.
"""
function crop(ga::GeoArray, cbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
    # Check if ga and cbox overlap
    if !bbox_overlap(bbox(ga), cbox)
        error("GeoArray and crop box do not overlap")
    end
    # Check extent and get bbox indices
    ga_x, ga_y, = size(ga)
    i_min_x, i_min_y = indices(ga, [cbox.min_x, cbox.max_y])
    i_max_x, i_max_y = indices(ga, [cbox.max_x, cbox.min_y])

    # Determine indices for crop area
    i_min_x = max(i_min_x, 1)
    i_max_x = min(i_max_x, ga_x)
    i_min_y = max(i_min_y, 1)
    i_max_y = min(i_max_y, ga_y)

    # Subset and return GeoArray
    return ga[i_min_x:i_max_x,i_min_y:i_max_y,:]
end

function crop(ga::GeoArray, cga::GeoArray)
    # Check if ga and crop ga CRS are the same
    if ga.crs != cga.crs
        error("GeoArrays have different CRS")
    end

    # Crop
    return crop(ga, bbox(cga))
end

ext: jl


text:
#
# The functions in this file manage the Newton-Krylov step and
# the Jacobian/preconditioner - vector products.
#
"""
Krylov_Step!(step, x, FS, FPS, ItRules, etag, delta = 0)

Take a Newton-Krylov step. This function does lots of its work
mapping nonlinear problems to linear solvers. Only then do I get
to deploy the Krylov linear solvers.
"""
function Krylov_Step!(step, x, FS, FPS, ItRules, etag, delta = 0)
    #
    # Test for too much, too soon.
    #
    lsolver = ItRules.lsolver
    lmaxit = ItRules.lmaxit
    T = eltype(FPS)
    kstep_test(FPS, step, lsolver)
#    (nk, mk) = size(FPS)
#    n = length(step)
#    n == nk || error("Krylov vectors wrong length")
#    (lsolver == "gmres") || error(lsolver, " ", "not supported")
    Jvec = ItRules.Jvec
    Pvec = ItRules.Pvec
    kl_store = ItRules.kl_store
    pdata = ItRules.pdata
    dx = ItRules.dx
    f = ItRules.f
    fixedeta = ItRules.fixedeta
    #    s0 = zeros(size(step))
    #
    #   Initial iterate for step is zero.
    #
    s0 = step
    s0 .*= 0.0
    side = ItRules.pside
    #
    # map the Jacobian-vector and preconditioner-vector products
    # from nsoli format to what the Krylov solvers want to see
    #
    kdata = (
        pdata = pdata,
        dx = dx,
        xc = x,
        f = f,
        FS = FS,
        Jvec = Jvec,
        Pvec = Pvec,
        delta = delta,
    )
    Pvecg = Pvec2
    Jvecg = Jvec2
    Pvec == nothing && (Pvecg = Pvec)
    Jvec == dirder && (Jvecg = Jvec)
    #
    #RHS=FS
    #T == Float64 || (RHS=T.(FS))
    if lsolver == "gmres"
    kout = kl_gmres( s0, FS, Jvecg, FPS, etag, Pvecg;
        kl_store = kl_store, pdata = kdata, side = side, lmaxit = lmaxit,)
    else
    kout = kl_bicgstab( s0, FS, Jvecg, FPS, etag, Pvecg;
        kl_store = kl_store, pdata = kdata, side = side, lmaxit = lmaxit,)
    end
    step .= -kout.sol
    reshist = kout.reshist
    lits = kout.lits
    idid = kout.idid
    Lstats = (reshist = reshist, lits = lits, idid = idid)
    return (step = step, Lstats = Lstats)
end

function Pvec2(v, kdata)
    F = kdata.f
    FS = kdata.FS
    xc = kdata.xc
    PV = kdata.Pvec
    pdata = kdata.pdata
    ptv = EvalPV(PV, v, xc, pdata)
    return ptv
end

function Jvec2(v, kdata)
    F = kdata.f
    FS = kdata.FS
    xc = kdata.xc
    JV = kdata.Jvec
    delta = kdata.delta
    pdata = kdata.pdata
    atv = EvalJV(JV, v, FS, xc, delta, pdata)
    return atv
end

function EvalPV(PV, v, xc, q::Nothing)
    ptv = PV(v, xc)
    return ptv
end

function EvalPV(PV, v, xc, pdata)
    ptv = PV(v, xc, pdata)
    return ptv
end

function EvalJV(JV, v, FS, xc, delta, q::Nothing)
    atv = JV(v, FS, xc)
    ptcmv!(atv, v, delta)
    #    if delta > 0
    #        atv .= atv + (1.0 / delta) * v
    #    end
    return atv
end

function EvalJV(JV, v, FS, xc, delta, pdata)
    atv = JV(v, FS, xc, pdata)
    ptcmv!(atv, v, delta)
    #    if delta > 0
    #        atv .= atv + (1.0 / delta) * v
    #    end
    return atv
end

function dirder(v, kdata)
    pdata = kdata.pdata
    dx = kdata.dx/norm(v)
    F = kdata.f
    FS = kdata.FS
    xc = kdata.xc
    delta = kdata.delta
    delx = copy(xc)
#    delx .= xc + dx * v
#    delx .+= dx*v
    axpy!(dx, v, delx)
    FPP = copy(xc)
    EvalF!(F, FPP, delx, pdata)
    atv = (FPP - FS) / dx
    ptcmv!(atv, v, delta)
    #    if delta > 0
    #        atv .= atv + (1.0 / delta) * v
    #    end
    return atv
end

function ptcmv!(atv, v, delta)
    (delta == 0.0) || (atv .= atv + (1.0 / delta) * v)
    #return atv
end

"""
forcing(itc, residratio, etag, ItRules, tol)

Compute the Eisenstat-Walker forcing term
"""
function forcing(itc, residratio, etag, ItRules, tol, resnorm)
    gamma = 0.9
    etamax = ItRules.eta
    fixedeta = ItRules.fixedeta
    if fixedeta || (itc == 0)
        etag = etamax
    else
        etaRes = gamma * (etag^2)
        etaA = gamma * (residratio^2)
        etaflim = 0.5 * tol / resnorm
        if etaRes <= 0.1
            etasafe = min(etamax, etaA)
        else
            etasafe = min(etamax, max(etaA, etaRes))
        end
        etag = min(etamax, max(etasafe, etaflim))
    end
    return etag
end

"""
kstore(n, lsolver)

Preallocates the vectors a Krylov method uses internally.
"""
function kstore(n, lsolver)
tmp1 = zeros(n)
tmp2 = zeros(n)
tmp3 = zeros(n)
tmp4 = zeros(n)
if lsolver=="gmres"
return (tmp1, tmp2, tmp3, tmp4)
else
tmp5 = zeros(n)
tmp6 = zeros(n)
tmp7 = zeros(n)
return (tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7)
end
end

function kstep_test(FPS, step, lsolver)
solver_ok = (lsolver == "gmres") || (lsolver == "bicgstab")
solver_ok || error(lsolver, " ", "not supported")
# Do a bit of management
(nk, ) = size(FPS)
n = length(step)
n == nk || error("Krylov storage vectors wrong length")
end


text:
using Cxx
include(Pkg.dir("Gallium","examples","call.jl"))
entrypoint = Gallium.entrypoint_for_julia_expression(dbg,"f","""
buf = IOBuffer()
Base.print_specialized_signature(buf,f)
takebuf_string(buf)
""")
lldb_exec(dbg,"f 7")
lldb_exec(dbg,"f 7")
ctx = Gallium.ctx(dbg)
function getFrameDescription(frame)
  li = Gallium.getASTForFrame(frame).ptr
  str = Gallium.call_prepared_entrypoint(dbg, ctx, entrypoint, Ptr{Void}[li])
  ptr = Gallium.target_call(dbg,:jl_bytestring_ptr,[str.ptr])
  size = Gallium.target_call(dbg,:jl_bytestring_length,[str.ptr])
  bytestring(Gallium.target_read(dbg,convert(UInt64,ptr),size))
end
thread = Gallium.current_thread(ctx)
for frame in thread
    Gallium.isJuliaFrame(frame) || continue
    println(getFrameDescription(frame))
end


text:
@active GlobalRef(x) begin
    if x isa GlobalRef
        (x.mod, x.name)
    else
        nothing
    end
end

@active Symbol(x) begin
    if x isa Symbol
        Some(string(x))
    else
        nothing
    end
end


text:
# T, TSSM_HANDLE
println("including tssm_schroedinger_generalized_laguerre.jl for type ", T)
for (METHOD, SUF, COMPLEX_METHOD, DIM) in (
                 (:SchroedingerGeneralizedLaguerre2D, :_schroedinger_gen_laguerre_2d, true, 2 ),
                 (:SchroedingerGeneralizedLaguerreHermite3D, :_schroedinger_gen_laguerre_hermite_3d, true, 3 ),
                 (:SchroedingerGeneralizedLaguerreReal2D, :_schroedinger_gen_laguerre_real_2d, false, 2 ),
                 (:SchroedingerGeneralizedLaguerreHermiteReal3D, :_schroedinger_gen_laguerre_hermite_real_3d, false, 3 ),
                )
println("    ", METHOD)
if T == :Float128
    SUF = Symbol(SUF, "_wf128")
end

if DIM==2
  @eval begin
    function ($METHOD)(T::Type{($T)}, ntheta::Integer, nfr::Integer,
                                      omega_r::Real, Omega::Real;
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_2D,
                       potential_t::Function=none_3D,
                       potential_t_derivative::Function=none_3D,
                       cubic_coupling::Real=0.0)
        with_potential = potential!=none_2D
        with_potential_t = potential_t!=none_3D
        with_potential_t_derivative = potential_t_derivative!=none_3D
        V_c = cfunction_check_return_type(potential, ($T), (($T),($T)))
        V_t_c = cfunction_check_return_type(potential_t, ($T), (($T),($T),($T)))
        V_t_derivative_c = cfunction_check_return_type(potential_t_derivative, ($T), (($T),($T),($T)))
        c = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new",SUF))), Ptr{Nothing},
                   (Int32, Int32, ($T), ($T),
                   ($T), ($T), Ptr{Nothing}, Bool, Ptr{Nothing}, Bool, Ptr{Nothing}, Bool, ($T)),
                   ntheta, nfr, omega_r, Omega,
                   hbar, mass, V_c, with_potential,
                   V_t_c, with_potential_t, V_t_derivative_c, with_potential_t_derivative,
                   cubic_coupling)
        m = ($METHOD){$T}(c)
        finalizer(x -> ccall( Libdl.dlsym(($TSSM_HANDLE),
                       $(string(PRE,"finalize",SUF))), Nothing, (Ptr{Nothing},), x.m), m)
        m
    end
  end # eval
  if T == :Float64
  @eval begin
    function ($METHOD)(ntheta::Integer, nfr::Integer,
                       omega_r::Real, Omega::Real;
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_2D,
                       potential_t::Function=none_3D,
                       potential_t_derivative::Function=none_3D,
                       cubic_coupling::Real=0.0)
        ($METHOD)(($T), ntheta, nfr, omega_r, Omega,
                       hbar=hbar, mass=mass, potential=potential,
                       potential_t=potential_t, potential_t_derivative=potential_t_derivative,
                       cubic_coupling=cubic_coupling)
    end
  end # eval
  end # if

elseif DIM==3
   @eval begin
     function ($METHOD)(T::Type{($T)}, ntheta::Integer, nfr::Integer,
                                      omega_r::Real, Omega::Real,
                                      nz::Integer, omega_z::Real;
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_3D,
                       potential_t::Function=none_4D,
                       potential_t_derivative::Function=none_4D,
                       cubic_coupling::Real=0.0)
        with_potential = potential!=none_3D
        with_potential_t = potential_t!=none_4D
        with_potential_t_derivative = potential_t_derivative!=none_4D
        V_c = cfunction_check_return_type(potential, ($T), (($T),($T),($T)))
        V_t_c = cfunction_check_return_type(potential_t, ($T), (($T),($T),($T),($T)))
        V_t_derivative_c = cfunction_check_return_type(potential_t_derivative, ($T), (($T),($T),($T)))
        c = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new",SUF))), Ptr{Nothing},
                   (Int32, Int32, ($T), ($T), Int32, ($T),
                   ($T), ($T), Ptr{Nothing}, Bool, Ptr{Nothing}, Bool, Ptr{Nothing}, Bool, ($T)),
                   ntheta, nfr, omega_r, Omega, nz, omega_z,
                   hbar, mass, V_c, with_potential,
                   V_t_c, with_potential_t, V_t_derivative_c, with_potential_t_derivative,
                   cubic_coupling)
        m = ($METHOD){$T}(c)
        finalizer(x -> ccall( Libdl.dlsym(($TSSM_HANDLE),
                       $(string(PRE,"finalize",SUF))), Nothing, (Ptr{Nothing},), x.m), m)
        m
     end
   end # eval
   if T == :Float64
   @eval begin
    function ($METHOD)(ntheta::Integer, nfr::Integer,
                       omega_r::Real, Omega::Real,
                       nz::Integer, omega_z::Real;
                       hbar::Real=1.0, mass::Real=1.0, potential::Function=none_3D,
                       potential_t::Function=none_4D,
                       potential_t_derivative::Function=none_4D,
                       cubic_coupling::Real=0.0)
        ($METHOD)(($T), ntheta, nfr, omega_r, Omega, nz, omega_z,
                       hbar=hbar, mass=mass, potential=potential,
                       potential_t=potential_t, potential_t_derivative=potential_t_derivative,
                       cubic_coupling=cubic_coupling)
    end
  end # eval
  end #if

end # if

@eval begin
    function get_ntheta(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_ntheta",SUF))), Int32,
             (Ptr{Nothing}, ), m.m )
    end

    function get_nfr(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nfr",SUF))), Int32,
             (Ptr{Nothing}, ), m.m )
    end

    function get_nr(m::($METHOD){$T})
      ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nr",SUF))), Int32,
            (Ptr{Nothing}, ), m.m )
    end

    function get_omega_r(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_omega_r",SUF))), ($T),
             (Ptr{Nothing}, ), m.m )
    end

    function get_Omega(m::($METHOD){$T})
       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_Omega",SUF))), ($T),
             (Ptr{Nothing}, ), m.m )
    end
end # eval

if DIM>=3
    @eval begin
        function get_nz(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nz",SUF))), Int32,
                 (Ptr{Nothing}, ), m.m )
        end

        function get_omega_z(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_omega_z",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end
    end # eval
end

if DIM==2
    @eval begin

        function get_weights(m::($METHOD){$T})
           dim =Array(Int32, 1)
           np = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_weights_r",SUF))), Ptr{$T},
                 (Ptr{Nothing}, Ptr{Int32}), m.m, dim)
           n = pointer_to_array(np, dim[1], false)
           copy(n)
        end

        function get_transformation_matrices(m::($METHOD){$T}, unsafe_access::Bool=false)
            dims =Array(Int32, 3)
            Lp = ccall( Libdl.dlsym(tssm_handle, $(string(PRE,"get_L",SUF))), Ptr{$T},
                (Ptr{Nothing}, Ptr{Int32}), m.m, dims )
            L = pointer_to_array(Lp, (dims[1], dims[2], dims[3]), false)
            if unsafe_access
                return L
            else
                return copy(L)
            end
        end

    end # eval

elseif DIM==3
    @eval begin

        function get_weights(m::($METHOD){$T})
           dim =Array(Int32, 1)
           np1 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_weights_r",SUF))), Ptr{$T},
                 (Ptr{Nothing}, Ptr{Int32}), m.m, dim)
           n1 = pointer_to_array(n1p, dim[1], false)
           np2 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_weights_z",SUF))), Ptr{$T},
                 (Ptr{Nothing}, Ptr{Int32}), m.m, dim)
           n2 = pointer_to_array(n2p, dim[1], false)
           copy(n1), copy(n2)
        end

        function get_transformation_matrices(m::($METHOD){$T}, unsafe_access::Bool=false)
            dims =Array(Int32, 3)
            Lp = ccall( Libdl.dlsym(tssm_handle, $(string(PRE,"get_L",SUF))), Ptr{$T},
                (Ptr{Nothing}, Ptr{Int32}), m.m, dims )
            H = pointer_to_array(Hp, (dims[1], dims[2], dims[3]), false)
            Hp = ccall( Libdl.dlsym(tssm_handle, $(string(PRE,"get_H_z",SUF))), Ptr{$T},
                (Ptr{Nothing}, Ptr{Int32}), m.m, dims )
            H = pointer_to_array(Hp, (dims[1], dims[2]), false)
            if unsafe_access
                return L, H
            else
                return copy(L), copy(H)
            end
        end

    end # eval
end #if

end # for


text:
# From SeisNoise
# https://github.com/tclements/SeisNoise.jl/blob/e5210fd11c1f5004c14f288f7ad0ac1cd495a2b0/test/test_slicing.jl#L44
@testset "slide daily data to 1-hour long segment" begin
    N = 86400 # 86400 seconds long signal
    fs = 100. # 100 Hz sampling rate
    A = rand(Float64, Int(N * fs))
    segment_len = 3600 # 1-hour long segment
    segment_step = 1800 # 50% overlapping
    starttime = d2u(DateTime(Date(now())))
    endtime = starttime + length(A) / fs - 1 / fs

    out, starts = slide(A, segment_len, segment_step, fs, starttime, endtime)
    # test first column
    @test out[:, 1] == A[1:Int(segment_len * fs)]
    # test overlap
    @test out[Int(segment_step*fs)+1:end, 1] == out[1:Int(segment_step*fs), 2]
    @test length(starts) == size(out,2)
    @test eltype(out) == eltype(A)
    # test number of slice
    @test size(starts, 1) == 47
end

@testset "slide 1-hour long segment to 15-minute long segment" begin
    N = 3600 # 3600 seconds long signal
    fs = 100. # 100 Hz sampling rate
    A = rand(Float64, Int(N * fs))
    segment_len = 900 # 15-minute long segment
    segment_step = 225 # 75% overlapping
    overlapping_area = segment_len / segment_step - 1
    starttime = d2u(DateTime(Date(now())))
    endtime = starttime + length(A) / fs - 1 / fs

    out, starts = slide(A, segment_len, segment_step, fs, starttime, endtime)
    # test first column
    @test out[:, 1] == A[1:Int(segment_len * fs)]
    # test overlap
    @test out[Int(segment_step*fs)+1:end, 1] == out[1:Int(segment_step*fs*overlapping_area), 2]
    @test length(starts) == size(out, 2)
    @test eltype(out) == eltype(A)
end

@testset "slide to 1-hour long ideal times" begin
    fs = 100.
    cc_len = 3600
    cc_step = 1800
    time_len = 86400
    starttime = DateTime(Date(now()))
    endtime = starttime + Second(time_len)

    expect_length = 47

    starts, ends = ideal_start_end(starttime, endtime, fs, cc_len, cc_step)
    # test number of ideal time slides
    @test size(starts, 1) == expect_length
    @test size(ends, 1) == expect_length
    # test starttime and endtime
    @test starts[1] == starttime
    @test ends[end] == endtime
end


text:
@enum InterpreterDataSource::UInt8 begin
    DATA_TAPE # variables, t in homotopy, instruction results
    DATA_PARAMETERS
    DATA_CONSTANTS
end

struct InterpreterArg
    data::InterpreterDataSource
    ind::Int32
end
function Base.show(io::IO, arg::InterpreterArg)
    if arg.data == DATA_TAPE
        print(io, "t[", arg.ind, "]")
    elseif arg.data == DATA_PARAMETERS
        print(io, "p[", arg.ind, "]")
    elseif arg.data == DATA_CONSTANTS
        print(io, "c[", arg.ind, "]")
    end
end

struct InterpreterInstruction
    op::InstructionOp
    args::NTuple{3,InterpreterArg}
    tape_index::Int
end

function Base.show(io::IO, instr::InterpreterInstruction)
    print(io, "t[$(instr.tape_index)] = ", string(Symbol(instr.op))[7:end], "(")
    if instr.op == INSTR_POW
        print(io, instr.args[1], ",", instr.args[2].ind)
    else
        join(io, instr.args[1:arity(instr.op)], ",")
    end
    print(io, ")")
end

struct Interpreter{T}
    instructions::Vector{InterpreterInstruction}
    constants::Vector{T}
    nexpressions::Int
    assignments::Vector{InterpreterArg}
    variables::Vector{Symbol}
    parameters::Vector{Symbol}
end

function Base.show(io::IO, I::Interpreter{T}) where {T}
    print(io, "Interpreter{$T} for $(length(I.instructions)) instructions")
end

"""
    show_instructions([io::IO], I::Interpreter)

Show the instruction executed by the interpreter `I`.
"""
show_instructions(I::Interpreter) = show_instructions(stdout, I)
function show_instructions(io::IO, I::Interpreter)
    print(io, "c = [")
    join(io, I.constants, ",")
    println("]")
    for i = 1:length(I.variables)
        println(io, "t[", i, "] = x[", i, "]")
    end
    for instr in I.instructions
        println(io, instr)
    end
    for (i, a) in enumerate(I.assignments)
        println(io, "out[", i, "] = ", a)
    end
end

mutable struct InterpreterCache{A<:AbstractArray,T}
    tape::A
    t₁::T
    t₂::T
end
Base.length(C::InterpreterCache) = length(C.tape)
cache_min_length(I::Interpreter) = maximum(instr -> instr.tape_index, I.instructions)
InterpreterCache(T::Type, I::Interpreter) = InterpreterCache(T, cache_min_length(I))
InterpreterCache(T::Type, n::Int) = InterpreterCache(zeros(T, n))
InterpreterCache(A::AbstractArray) = InterpreterCache(A, nothing, nothing)
InterpreterCache(A::AcbRefVector) =
    InterpreterCache(A, Acb(prec = precision(A)), Acb(prec = precision(A)))
InterpreterCache(A::ArbRefVector) =
    InterpreterCache(A, Arb(prec = precision(A)), Arb(prec = precision(A)))

InterpreterCache(::Type{Acb}, I::Interpreter; prec::Int) =
    InterpreterCache(Acb, cache_min_length(I), prec = prec)
InterpreterCache(::Type{Acb}, n::Int; prec::Int) =
    InterpreterCache(AcbRefVector(n, prec = prec))
InterpreterCache(::Type{Arb}, I::Interpreter; prec::Int) =
    InterpreterCache(Arb, cache_min_length(I), prec = prec)
InterpreterCache(::Type{Arb}, n::Int; prec::Int) =
    InterpreterCache(ArbRefVector(n, prec = prec))
function set_arb_precision!(A::InterpreterCache, prec::Int)
    A.tape = setprecision(A.tape, prec)
    A.t₁ = setprecision(A.t₁, prec)
    A.t₂ = setprecision(A.t₂, prec)
    A
end

function InterpreterArg(instr_map, var_map, t, param_map, constants, arg)
    if haskey(instr_map, arg)
        return InterpreterArg(DATA_I, instr_map[arg])
    elseif haskey(var_map, arg)
        return InterpreterArg(DATA_X, var_map[arg])
    elseif haskey(param_map, arg)
        return InterpreterArg(DATA_P, param_map[arg])
    elseif t == arg
        return InterpreterArg(DATA_T, 0)
    else # constant
        c = findfirst(c -> c == arg, constants)
        if isnothing(c)
            push!(constants, arg)
            return InterpreterArg(DATA_C, length(constants))
        else
            return InterpreterArg(DATA_C, c)
        end
    end
end

function Interpreter(
    list::InstructionList,
    assignments;
    variables::Vector{Symbol},
    t::Union{Nothing,Symbol} = nothing,
    parameters::Vector{Symbol} = Symbol[],
    nexpressions::Int = length(assignments),
)
    # compress needed storage
    assignment_tape_indices = Set{Int}()
    for a in assignments
        if a isa InstructionRef
            push!(assignment_tape_indices, a.i)
        end
    end
    alloc_alias = allocation_alias(list, assignment_tape_indices)
    # alloc_alias = collect(1:length(list))
    # we shift all allocation alias by the number of variables (plus t if applicable)
    # since we will always put the variables at the start of the tape
    shift = length(variables) + !isnothing(t)
    alloc_alias .= alloc_alias .+ shift

    var_map = Dict(zip(variables, 1:length(variables)))
    if !isnothing(t)
        var_map[t] = shift
    end
    param_map = Dict(zip(parameters, 1:length(parameters)))

    constants_map = Dict{Any,Int}()
    nconstants = 0
    constants_type = Int32
    instructions = map(list) do (instr_ref, instr)
        args = ntuple(Val(3)) do i
            arg = instr.args[i]
            if isnothing(arg)
                # dummy data
                InterpreterArg(DATA_TAPE, Int32(shift + 1))
            elseif instr.op == INSTR_POW && i == 2
                # We store the power in the INDEX of the argument
                InterpreterArg(DATA_TAPE, Int32(arg))
            elseif haskey(var_map, arg)
                InterpreterArg(DATA_TAPE, Int32(var_map[arg]))
            elseif haskey(param_map, arg)
                InterpreterArg(DATA_PARAMETERS, Int32(param_map[arg]))
            elseif arg isa InstructionRef
                InterpreterArg(DATA_TAPE, Int32(alloc_alias[arg.i]))
            else
                if haskey(constants_map, arg)
                    InterpreterArg(DATA_CONSTANTS, constants_map[arg])
                else
                    nconstants += 1
                    push!(constants_map, arg => nconstants)
                    constants_type = promote_type(typeof(arg), constants_type)
                    InterpreterArg(DATA_CONSTANTS, nconstants)
                end
            end
        end
        InterpreterInstruction(instr.op, args, alloc_alias[instr_ref.i])
    end
    constants = Vector{constants_type}(undef, nconstants)
    for (k, ind) in constants_map
        constants[ind] = k
    end
    assignment_instructions = InterpreterArg[]
    for a in assignments
        if haskey(var_map, a)
            push!(assignment_instructions, InterpreterArg(DATA_TAPE, var_map[a]))
        elseif haskey(param_map, a)
            push!(assignment_instructions, InterpreterArg(DATA_PARAMETERS, param_map[a]))
        elseif haskey(constants_map, a)
            push!(assignment_instructions, InterpreterArg(DATA_CONSTANTS, constants_map[a]))
        elseif a isa InstructionRef
            push!(assignment_instructions, InterpreterArg(DATA_TAPE, alloc_alias[a.i]))
        elseif a isa Number
            nconstants += 1
            push!(constants_map, a => nconstants)
            push!(constants, a)
            push!(assignment_instructions, InterpreterArg(DATA_CONSTANTS, nconstants))
        else
            error(string("Unhandled assignment: ", sprint(show, a)))
        end
    end
    Interpreter(
        instructions,
        constants,
        nexpressions,
        assignment_instructions,
        variables,
        parameters,
    )
end

function interpreter(F::Union{System,Homotopy})
    vars = Symbol.(variables(F))
    params = Symbol.(parameters(F))
    t = F isa Homotopy ? Symbol(F.t) : nothing
    list, out = instruction_list(expressions(F))
    ModelKit.Interpreter(list, out; variables = vars, parameters = params, t = t)
end
function jacobian_interpreter(F::Union{System,Homotopy})
    vars = Symbol.(variables(F))
    params = Symbol.(parameters(F))
    t = F isa Homotopy ? Symbol(F.t) : nothing
    list, out = instruction_list(expressions(F))
    nexpressions = length(out)
    list, out, jac_out = gradient(list, out, vars)
    out′ = [out; vec(jac_out)]
    ModelKit.Interpreter(
        list,
        out′;
        variables = vars,
        parameters = params,
        t = t,
        nexpressions = nexpressions,
    )
end

function load_data(f, var_name, arg, has_parameters::Bool)
    v1, v2, v3 = gensym(var_name), gensym(var_name), gensym(var_name)
    if has_parameters
        quote
            if $arg.data == DATA_TAPE
                $v1 = tape[$arg.ind]
                $(f(v1))
            elseif $arg.data == DATA_PARAMETERS
                $v2 = parameters[$arg.ind]
                $(f(v2))
            else #if $arg.data == DATA_CONSTANTS
                $v3 = constants[$arg.ind]
                $(f(v3))
            end
        end
    else
        quote
            if $arg.data == DATA_TAPE
                $v1 = tape[$arg.ind]
                $(f(v1))
            else #if $arg.data == DATA_CONSTANTS
                $v3 = constants[$arg.ind]
                $(f(v3))
            end
        end
    end
end

Base.@propagate_inbounds neg_instr!(z, i, x) = (z[i] = neg_instr(x))
Base.@propagate_inbounds sqr_instr!(z, i, x) = (z[i] = sqr_instr(x))
Base.@propagate_inbounds sin_instr!(z, i, x) = (z[i] = sin(x))
Base.@propagate_inbounds cos_instr!(z, i, x) = (z[i] = cos(x))
Base.@propagate_inbounds pow_instr!(z, i, x, k::Integer) = (z[i] = pow_instr(x, k))

Base.@propagate_inbounds add_instr!(z, i, x, y, _t₁) = (z[i] = add_instr(x, y))
Base.@propagate_inbounds sub_instr!(z, i, x, y, _t₁) = (z[i] = sub_instr(x, y))
Base.@propagate_inbounds mul_instr!(z, i, x, y, _t₁) = (z[i] = mul_instr(x, y))
Base.@propagate_inbounds div_instr!(z, i, x, y, _t₁) = (z[i] = div_instr(x, y))

Base.@propagate_inbounds muladd_instr!(z, i, u, v, w, _t₁, _t₂) =
    (z[i] = muladd_instr(u, v, w))
Base.@propagate_inbounds mulsub_instr!(z, i, u, v, w, _t₁, _t₂) =
    (z[i] = mulsub_instr(u, v, w))
Base.@propagate_inbounds submul_instr!(z, i, u, v, w, _t₁, _t₂) =
    (z[i] = submul_instr(u, v, w))

# Arb versions
opt_set!(z::Union{Acb,AcbRef}, x::Union{Acb,AcbRef}) = x
opt_set!(z::Union{Acb,AcbRef}, x::Number) = (z[] = x; z)

neg_instr!(z::AcbRefVector, i, x) = Arblib.neg!(z[i], opt_set!(z[i], x))
sqr_instr!(z::AcbRefVector, i, x) = Arblib.sqr!(z[i], opt_set!(z[i], x))
sin_instr!(z::AcbRefVector, i, x) = Arblib.sin!(z[i], opt_set!(z[i], x))
cos_instr!(z::AcbRefVector, i, x) = Arblib.cos!(z[i], opt_set!(z[i], x))
pow_instr!(z::AcbRefVector, i, x, k::Integer) = Arblib.pow!(z[i], opt_set!(z[i], x), k)

add_instr!(z::AcbRefVector, i, x, y, t₁) =
    Arblib.add!(z[i], opt_set!(z[i], x), opt_set!(t₁, y))
sub_instr!(z::AcbRefVector, i, x, y, t₁) =
    Arblib.sub!(z[i], opt_set!(z[i], x), opt_set!(t₁, y))
mul_instr!(z::AcbRefVector, i, x, y, t₁) =
    Arblib.mul!(z[i], opt_set!(z[i], x), opt_set!(t₁, y))
div_instr!(z::AcbRefVector, i, x, y, t₁) =
    Arblib.div!(z[i], opt_set!(z[i], x), opt_set!(t₁, y))

function muladd_instr!(z::AcbRefVector, i, u, v, w, t₁, t₂)
    z[i][] = w
    Arblib.addmul!(z[i], opt_set!(t₁, u), opt_set!(t₂, v))
end
function mulsub_instr!(z::AcbRefVector, i, u, v, w, t₁, t₂)
    Arblib.neg!(z[i], opt_set!(t₁, w))
    Arblib.addmul!(z[i], opt_set!(t₁, u), opt_set!(t₂, v))
end
function submul_instr!(z::AcbRefVector, i, u, v, w, t₁, t₂)
    z[i][] = w
    Arblib.submul!(z[i], opt_set!(t₁, u), opt_set!(t₂, v))
end

function evaluate_block(; has_parameters::Bool = false)
    # This should get compiled to a jump table by LLVM, so order doesn't matter
    load_data(:v₁, :arg₁, has_parameters) do v₁
        quote
            if op == INSTR_NEG
                neg_instr!(tape, i, $v₁)
            elseif op == INSTR_SQR
                sqr_instr!(tape, i, $v₁)
            elseif op == INSTR_SIN
                sin_instr!(tape, i, $v₁)
            elseif op == INSTR_COS
                cos_instr!(tape, i, $v₁)
            elseif op == INSTR_POW
                pow_instr!(tape, i, $v₁, arg₂.ind)
            else
                $(
                    load_data(:v₂, :arg₂, has_parameters) do v₂
                        quote
                            if op == INSTR_ADD
                                add_instr!(tape, i, $v₁, $v₂, t₁)
                            elseif op == INSTR_SUB
                                sub_instr!(tape, i, $v₁, $v₂, t₁)
                            elseif op == INSTR_MUL
                                mul_instr!(tape, i, $v₁, $v₂, t₁)
                            elseif op == INSTR_DIV
                                div_instr!(tape, i, $v₁, $v₂, t₁)
                            else
                                $(
                                    load_data(:v₃, :arg₃, has_parameters) do v₃
                                        quote
                                            if op == INSTR_MULADD
                                                muladd_instr!(
                                                    tape,
                                                    i,
                                                    $v₁,
                                                    $v₂,
                                                    $v₃,
                                                    t₁,
                                                    t₂,
                                                )
                                            elseif op == INSTR_MULSUB
                                                mulsub_instr!(
                                                    tape,
                                                    i,
                                                    $v₁,
                                                    $v₂,
                                                    $v₃,
                                                    t₁,
                                                    t₂,
                                                )
                                            elseif op == INSTR_SUBMUL
                                                submul_instr!(
                                                    tape,
                                                    i,
                                                    $v₁,
                                                    $v₂,
                                                    $v₃,
                                                    t₁,
                                                    t₂,
                                                )
                                            else
                                                throw(
                                                    error(
                                                        string(
                                                            "Unexpected fall-through. instr: ",
                                                            Symbol(op),
                                                        ),
                                                    ),
                                                )
                                            end
                                        end
                                    end
                                )
                            end
                        end
                    end
                )
            end
        end
    end
end

function taylor_block(; has_parameters::Bool = false)
    # This should get compiled to a jump table by LLVM, so order doesn't matter
    three_arg =
        (v₁, v₂) -> load_data(:v₃, :arg₃, has_parameters) do v₃
            quote
                if op == INSTR_MULADD
                    tape[i] = taylor_muladd(Order, $v₁, $v₂, $v₃)
                elseif op == INSTR_MULSUB
                    tape[i] = taylor_mulsub(Order, $v₁, $v₂, $v₃)
                elseif op == INSTR_SUBMUL
                    tape[i] = taylor_submul(Order, $v₁, $v₂, $v₃)
                else
                    throw(error(string("Unexpected fall-through. instr: ", Symbol(op))))
                end
            end
        end

    two_arg = v₁ -> load_data(:v₂, :arg₂, has_parameters) do v₂
        quote
            if op == INSTR_ADD
                tape[i] = taylor_add(Order, $v₁, $v₂)
            elseif op == INSTR_SUB
                tape[i] = taylor_sub(Order, $v₁, $v₂)
            elseif op == INSTR_MUL
                tape[i] = taylor_mul(Order, $v₁, $v₂)
            elseif op == INSTR_DIV
                tape[i] = taylor_div(Order, $v₁, $v₂)
            else
                $(three_arg(v₁, v₂))
            end
        end
    end
    load_data(:v₁, :arg₁, has_parameters) do v₁
        quote
            if op == INSTR_NEG
                tape[i] = taylor_neg(Order, $v₁)
            elseif op == INSTR_SQR
                tape[i] = taylor_sqr(Order, $v₁)
            elseif op == INSTR_SIN
                tape[i] = taylor_sin(Order, $v₁)
            elseif op == INSTR_COS
                tape[i] = taylor_cos(Order, $v₁)
            elseif op == INSTR_POW
                tape[i] = taylor_pow(Order, $v₁, arg₂.ind)
            else
                $(two_arg(v₁))
            end
        end
    end
end

for has_parameters in [true, false], has_t in [true, false]
    @eval function execute!(
        u::AbstractArray,
        I::Interpreter,
        x::AbstractArray,
        $((has_t ? (:t,) : ())...),
        $(has_parameters ? :(parameters::AbstractArray) : :(parameters::Nothing)),
        cache::InterpreterCache,
    )
        isnothing(parameters) || checkbounds(parameters, 1:length(I.parameters))

        constants = I.constants
        tape = cache.tape
        D = length(I.variables)
        for i = 1:D
            tape[i] = x[i]
        end
        $((has_t ? (:(tape[D+1] = t),) : ())...)
        t₁ = cache.t₁
        t₂ = cache.t₂
        @inbounds for instr in I.instructions
            i = instr.tape_index
            op = instr.op
            arg₁, arg₂, arg₃ = instr.args
            $(evaluate_block(has_parameters = has_parameters))
        end
        for (k, assign_arg) in enumerate(I.assignments)
            $(
                load_data(:v, :assign_arg, has_parameters) do v
                    :(u[k] = $v)
                end
            )
        end
        u
    end

    @eval function execute!(
        u::Union{Nothing,AbstractArray},
        U::AbstractArray,
        I::Interpreter,
        x::AbstractArray,
        $((has_t ? (:t,) : ())...),
        $(has_parameters ? :(parameters::AbstractArray) : :(parameters::Nothing)),
        cache::InterpreterCache,
    )
        isnothing(parameters) || checkbounds(parameters, 1:length(I.parameters))

        constants = I.constants
        tape = cache.tape
        D = length(I.variables)
        for i = 1:D
            tape[i] = x[i]
        end
        $((has_t ? (:(tape[D+1] = t),) : ())...)
        t₁ = cache.t₁
        t₂ = cache.t₂
        @inbounds for instr in I.instructions
            i = instr.tape_index
            op = instr.op
            arg₁, arg₂, arg₃ = instr.args
            $(evaluate_block(has_parameters = has_parameters))
        end
        if !(u isa Nothing)
            for k = 1:I.nexpressions
                a = I.assignments[k]
                $(
                    load_data(:v, :a, has_parameters) do v
                        :(u[k] = $v)
                    end
                )
            end
        end
        i = I.nexpressions
        for l = 1:D, k = 1:I.nexpressions
            i += 1
            a = I.assignments[i]
            $(
                load_data(:v, :a, has_parameters) do v
                    :(U[k, l] = $v)
                end
            )
        end
        nothing
    end

    @eval function execute!(
        u::A,
        Order::Val{K},
        I::Interpreter,
        x::AbstractArray,
        $((has_t ? (:t,) : ())...),
        $(has_parameters ? :(parameters::AbstractArray) : :(parameters::Nothing)),
        cache::InterpreterCache,
    ) where {A<:AbstractArray,K}
        isnothing(parameters) || checkbounds(parameters, 1:length(I.parameters))

        constants = I.constants
        tape = cache.tape
        D = length(I.variables)
        for i = 1:D
            tape[i] = x[i]
        end
        $((has_t ? (:(tape[D+1] = (t, one(t))),) : ())...)
        @inbounds for instr in I.instructions
            i = instr.tape_index
            op = instr.op
            arg₁, arg₂, arg₃ = instr.args
            $(taylor_block(has_parameters = has_parameters))
        end
        if u isa TaylorVector
            for (k, assign_arg) in enumerate(I.assignments)
                $(
                    load_data(:v, :assign_arg, has_parameters) do v
                        :(u[k] = $v)
                    end
                )
            end
        else
            for (k, assign_arg) in enumerate(I.assignments)
                $(
                    load_data(:v, :assign_arg, has_parameters) do v
                        :(u[k] = last($v))
                    end
                )
            end
        end
        u
    end
end


text:
# Graph to matrices

using Graphs
using Test
using SparseArrays

# data

global eds = [Edge(1,1,2), Edge(2,1,3), Edge(3,2,4), Edge(4,3,4)]

global a0 = fill(false, 4, 4)
global a0[1, 2] = a0[1, 3] = a0[2, 4] = a0[3, 4] = true

global a0u = copy(a0)
global a0u[2, 1] = a0u[3, 1] = a0u[4, 2] = a0u[4, 3] = true

global a0_s = sparse(a0)
global a0u_s = sparse(a0u)

# graphs

global gd_elst = simple_edgelist(4, eds)
global gu_elst = simple_edgelist(4, eds; is_directed=false)

global gd_ilst = simple_inclist(4)
global gu_ilst = simple_inclist(4, is_directed=false)

let
for ed in eds
  add_edge!(gd_ilst, ed.source, ed.target)
  add_edge!(gu_ilst, ed.source, ed.target)
end
end

# adjacency matrix

@test adjacency_matrix(gd_elst) == a0
@test adjacency_matrix(gu_elst) == a0u

@test adjacency_matrix_sparse(gd_elst) == a0_s
@test adjacency_matrix_sparse(gu_elst) == a0u_s

@test adjacency_matrix(gd_ilst) == a0
@test adjacency_matrix(gu_ilst) == a0u

@test adjacency_matrix_sparse(gd_ilst) == a0_s
@test adjacency_matrix_sparse(gu_ilst) == a0u_s

# weight matrix

global eweights = [1., 2., 3., 4.]
global wm0  = [0. 1. 2. 0.; 0. 0. 0. 3.; 0. 0. 0. 4.; 0. 0. 0. 0.]
global wm0u = [0. 1. 2. 0.; 1. 0. 0. 3.; 2. 0. 0. 4.; 0. 3. 4. 0.]
global wm0_s = sparse(wm0)
global wm0u_s = sparse(wm0u)

@test weight_matrix(gd_elst, eweights) == wm0
@test weight_matrix(gu_elst, eweights) == wm0u

@test weight_matrix_sparse(gd_elst, eweights) == wm0_s
@test weight_matrix_sparse(gu_elst, eweights) == wm0u_s

@test weight_matrix(gd_ilst, eweights) == wm0
@test weight_matrix(gu_ilst, eweights) == wm0u

@test weight_matrix_sparse(gd_ilst, eweights) == wm0_s
@test weight_matrix_sparse(gu_ilst, eweights) == wm0u_s

# distance matrix

global dm0  = [0. 1. 2. Inf; Inf 0. Inf 3.; Inf Inf 0. 4.; Inf Inf Inf 0.]
global dm0u = [0. 1. 2. Inf; 1. 0. Inf 3.; 2. Inf 0. 4.; Inf 3. 4. 0.]

@test distance_matrix(gd_elst, eweights) == dm0
@test distance_matrix(gu_elst, eweights) == dm0u

@test distance_matrix(gd_ilst, eweights) == dm0
@test distance_matrix(gu_ilst, eweights) == dm0u

# Laplacian matrix

global L0 = [2. -1. -1. 0.; -1. 2. 0. -1.; -1. 0. 2. -1.; 0. -1. -1. 2.]
global L0_s = sparse(L0)

@test laplacian_matrix(gu_elst) == L0
@test laplacian_matrix(gu_ilst) == L0

@test laplacian_matrix_sparse(gu_elst) == L0_s
@test laplacian_matrix_sparse(gu_ilst) == L0_s

global Lw = [3. -1. -2. 0.; -1. 4. 0. -3.; -2. 0. 6. -4.; 0. -3. -4. 7.]
global Lw_s = sparse(Lw)

@test laplacian_matrix(gu_elst, eweights) == Lw
@test laplacian_matrix(gu_ilst, eweights) == Lw

@test laplacian_matrix_sparse(gu_elst, eweights) == Lw_s
@test laplacian_matrix_sparse(gu_ilst, eweights) == Lw_s


text:
using Documenter, SuperResolution

makedocs(;
    modules=[SuperResolution],
    sitename="SuperResolution.jl",
    authors="Patrick Schneider",
    #repo="https://github.com/PGSch/SuperResolution.jl/blob/{commit}{path}#L{line}",
    format=Documenter.HTML(prettyurls=false),#;
    #     prettyurls=get(ENV, "CI", "false") == "true",
    #     canonical="https://PGSch.github.io/SuperResolution.jl",
    #     assets=String[],
    # ),
    pages=[
        "Home" => "index.md",
       #  "Algorithms" => [
       #     "algorithms.md",
       #     "init.md",
       #     "kmeans.md",
       #     "kmedoids.md",
       #     "hclust.md",
       #     "mcl.md",
       #     "affprop.md",
       #     "dbscan.md",
       #     "fuzzycmeans.md",
       # ],
       # "validate.md",
    ],
)

deploydocs(;
    repo="github.com/PGSch/SuperResolution.jl",
)

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
# Autogenerated wrapper script for QPALM_jll for x86_64-linux-musl
export libladel, libqpalm

using OpenBLAS32_jll
using SuiteSparse_jll
JLLWrappers.@generate_wrapper_header("QPALM")
JLLWrappers.@declare_library_product(libladel, "libladel.so")
JLLWrappers.@declare_library_product(libqpalm, "libqpalm.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, SuiteSparse_jll)
    JLLWrappers.@init_library_product(
        libladel,
        "lib/libladel.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqpalm,
        "lib/libqpalm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Clang.cindex

base_classes(c) = search(c, isa(c, CursorKind.CXXBaseSpecifier))

function find_sym(name,liblist)
    libs = {dlopen(l)=>l for l in liblist}
    found = false
    dl = C_NULL
    for dl in keys(libs)
        try
            Libdl.dlsym_e(dl, name)
        catch
            continue
        end
        found = true
    end
    if found
        println("symbol $name found in ", libs[dl])
        return splitdir(libs[dl])[end]
    else
        return Union{}
        println("NOT FOUND")
    end
end

###############################################################################
# Extra helpers
#	These functions require libwrapclang to be compiled with
#	-DUSE_CLANG_CPP. They provide some extra information directly
#	from the Clang C++ API.
###############################################################################

function method_vt_index(cursor::cindex.CXCursor)
    !(cindex.CXXMethod_isVirtual(cursor) == 1) && return -1
    ccall(("wci_getCXXMethodVTableIndex", :libwrapclang), Int32,
		  (Ptr{UInt8},), cursor.data)
end

function method_mangled_name(cursor::cindex.CXCursor)
    bufr = zeros(UInt8, 1024)
    ccall(("wci_getCXXMethodMangledName", :libwrapclang), Int32,
		 (Ptr{UInt8},Ptr{UInt8}), cursor.data, bufr)
    unsafe_string(convert(Ptr{UInt8},bufr))
end


text:
struct DetectionCallback
    detection_times::Vector{OptTimePoint}
    detection_types::Vector{UInt8}

    tracing_times::Vector{OptTimePoint}
    tracing_sources::Vector{UInt32}
    tracing_types::Vector{UInt8}

    max_num_infected::UInt32
    time_limit::MocosSim.TimePoint
end

DetectionCallback(sz::Integer, max_num_infected::Integer=10^8, time_limit::MocosSim.TimePoint=typemax(MocosSim.TimePoint)) = DetectionCallback(
    Vector{OptTimePoint}(missing, sz),
    fill(UInt8(0), sz),
    Vector{OptTimePoint}(missing, sz),
    fill(UInt32(0), sz),
    fill(UInt8(0), sz),
    max_num_infected,
    time_limit
)

function reset!(cb::DetectionCallback)
  fill!(cb.detection_times, missing)
  fill!(cb.detection_types, 0)
  fill!(cb.tracing_sources, 0)
  fill!(cb.tracing_types, 0)
end

function (cb::DetectionCallback)(event::MocosSim.Event, state::MocosSim.SimState, params::MocosSim.SimParams)
  eventkind = MocosSim.kind(event)
  contactkind = MocosSim.contactkind(event)
  subject = MocosSim.subject(event)
  if MocosSim.isdetection(eventkind)
    cb.detection_times[subject] = MocosSim.time(event)
    cb.detection_types[subject] = MocosSim.detectionkind(event) |> UInt8
  elseif MocosSim.istracing(eventkind)
    cb.tracing_times[subject] = MocosSim.time(event)
    cb.tracing_sources[subject] = MocosSim.source(event)
    cb.tracing_types[subject] = MocosSim.tracingkind(event) |> UInt8
  end
  return MocosSim.numinfected(state.stats) < cb.max_num_infected && MocosSim.time(event) < cb.time_limit
end

function saveparams(dict, cb::DetectionCallback, prefix::AbstractString="")
  dict[prefix*"detection_times"] = optreal2float32.(cb.detection_times)
  dict[prefix*"detection_types"] = cb.detection_types

  dict[prefix*"tracing_times"] = optreal2float32.(cb.tracing_times)
  dict[prefix*"tracing_sources"] = cb.tracing_sources
  dict[prefix*"tracing_types"] = cb.tracing_types
end

function save_infections_and_detections(path::AbstractString, simstate::MocosSim.SimState, callback::DetectionCallback)
  f = jldopen(path, "w", compress=true)
  try
    MocosSim.saveparams(f, simstate)
    saveparams(f, callback)
  finally
    close(f)
  end
  nothing
end


text:
import Base.parse

# Part 1
struct Policy
    min::Int32
    max::Int32
    char::Char
end

function parse(::Type{Policy}, x)
    min, max, char = split(x, ('-', ' '))
    Policy(parse(Int32, min), parse(Int32, max), only(char))
end

# Load input
parse_line(policy, password) = (parse(Policy, policy), password)
policies_passwords = open("day02_input.txt") do f
    [parse_line(strip.(split(line, ':'))...) for line in eachline(f)]
end

is_valid(policy, password) = policy.min <= count(==(policy.char), password) <= policy.max
println("There are $(sum([is_valid(pp...) for pp in policies_passwords])) valid passwords")

# Part 2
new_is_valid(policy, password) = xor(password[policy.min] == policy.char, password[policy.max] == policy.char)
println("There are now $(sum([new_is_valid(pp...) for pp in policies_passwords])) valid passwords")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
using StaticArrays
using ArraysOfArrays

function split_vec_of_arrays(u)
    vec.(u) |>
    x ->
    VectorOfSimilarVectors(x).data |>
    transpose |>
    VectorOfSimilarVectors
end


text:
function plotModel(m,plotting::Bool=false,includeMeshInfo::Bool=false,Msh = [],cutPad::Int64 = 0,limits = [],filename="")

if limits!=[]
	vmin = limits[1];
	vmax = limits[2];
else
	vmin = minimum(m);
	vmax = maximum(m);
end

limits = tuple([vmax,vmin]...);
if cutPad >= 0
	m,Msh = cutAbsorbingLayer(m,Msh,cutPad);
end
if plotting
	if length(size(m))==2
		T = m';
		imshow(T, clim = limits); colorbar();
		if includeMeshInfo
			Omega = Msh.domain;
			tics = 0:2:floor(Int64,Omega[2])
			xticks(linspace(0,size(m,1)*tics[end]/Omega[2],length(tics)),tics);
			xlabel("Lateral distance (km)");
			tics = 0:floor(Int64,Omega[4])
			# println(tics)
			# println(linspace(0,size(m,2)*tics[end]/Omega[4],length(tics)))
			yticks(linspace(0,size(m,2)*tics[end]/Omega[4],length(tics)),tics);
			ylabel("Depth (km)")
		end
		if filename != ""
			# title(filename[1:end-4],fontsize = 11);
			savefig(string(filename[1:end-4],".png"));
		end
	elseif length(size(m))==3
		lin = zeros(Int64,16);
		v = m;
		for k=1:16
			lin[k] = convert(Int64,round(k*(size(m,2)/16)));
		end
		for k=1:16
			subplot(4,4,k)
			pic = reshape(v[:,lin[k],:],size(m,1),size(m,3))';
			imshow(pic,clim = limits); title(string("frame",lin[k]));colorbar()
		end
		if filename != ""
			# title(filename[1:end-4],fontsize = 11);
			savefig(string(filename[1:end-4],".png"));
		end
	end
end
end


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

name = "GAP"
version = v"4.11.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/gap-system/gap.git", "a20c8f40883c8656a240317d732229dbe7c3b5ab"),
#    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(version)/gap-$(version)-core.tar.bz2",
#                  "6637f66409bc91af21eaa38368153270b71b13b55b75cc1550ed867c629901d1"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gap*

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/configure.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/Makefile.patch

# run autogen.sh if compiling from it source and/or if configure was patched
./autogen.sh

# compile GAP
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-gmp=${prefix} \
    --with-readline=${prefix} \
    --with-zlib=${prefix} \
    --with-gc=julia \
    --with-julia
make -j${nproc}

# install GAP binaries
make install-bin install-headers install-libgap

# get rid of the wrapper shell script, which is useless for us
mv ${WORKSPACE}/destdir/bin/gap.real ${WORKSPACE}/destdir/bin/gap

# install sysinfo.gap but patch out
mkdir -p ${WORKSPACE}/destdir/share/gap/
sed -e 's;$PWD;@GAPROOT@;g' sysinfo.gap > ${WORKSPACE}/destdir/share/gap/sysinfo.gap

# We deliberately do NOT install the GAP library, documentation, etc. because
# they are identical across all platforms; instead, we use another platform
# independent artifact to ship them to the user.
"""

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

# we only care about 64bit builds
filter!(p -> nbits(p) == 64, platforms)

# Windows is not supported
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("gap", :gap),
    LibraryProduct("libgap", :libgap),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("Readline_jll"),
    Dependency("Zlib_jll"),

    # GAP tries hard to produce a binary that works in all Julia versions,
    # regardless of which version of Julia it was compiled again; so the
    # version restriction below could be dropped or changed if necessary
    BuildDependency(PackageSpec(name="libjulia_jll", version=v"1.4.2")),
]

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


text:
"""
Set new world coordinates and calculate transformation data
"""
function _world!(p,xmin,xmax,ymin,ymax)
    p.xmin=xmin
    p.xmax=xmax
    p.ymin=ymin
    p.ymax=ymax
    p.ax= p.w/(p.xmax-p.xmin);
    p.ay=-p.h/(p.ymax-p.ymin);
    p.bx=0   - p.ax *p.xmin;
    p.by=p.h - p.ay *p.ymin;
    nothing
end

"""
Transform a pair of coordinates from world to canvas
"""
_tran2d(p,x,y)=(x*p.ax+p.bx,y*p.ay+p.by)

"""
    command!(p::CanvasPlot,cmd)

Enter new command named `cmd`.

The idea is to pass one single Dict
with all plot data to javascript.

For this purpose we need a some "language".

Plot elements are described by commands executed
one after another. For this purpose, we use
integers converted to strings as dict keys and pass the command
name as the corresponding entry.

The entry `cmdcount` keeps track of the number of
commands. So we can parse all commands from 1 to jsdict[:cmdcount]
in javascript in the same sequence as they have been entered.

Parameters are named and entered into the dictionary with the
command number as prefix.

E.g. for a polyline as command number 5, we create the entres

```
"5" => "polyline"
"5_x" => Vector of x coordinates in canvas coordinate system
"5_y" => Vector of y coordinates in canvas coordinate system
"""
function command!(p,cmd)
    p.jsdict["cmdcount"]=p.jsdict["cmdcount"]+1
    pfx=string(p.jsdict["cmdcount"])
    p.jsdict[pfx]=cmd
    pfx
end

"""
Pass pair of coordinate arrays for `lines!`,`polyline!`,`polygon!`
"""
function _poly!(p::CanvasPlot,cmd,x,y)
    pfx=command!(p,cmd)

    tx=Vector{Float32}(undef,length(x))
    ty=Vector{Float32}(undef,length(y))
    for i=1:length(x)
        tx[i],ty[i]=_tran2d(p,x[i],y[i])
    end
    p.jsdict[pfx*"_x"]=tx
    p.jsdict[pfx*"_y"]=ty
    p

end


text:
struct Rows{Row, Dimensions, Columns, Names} <: AbstractArray{Row, Dimensions}
    columns::Columns
    names::Names
end
@inline function compare_axes(reference_axes, item)
    axes(item) == reference_axes
end
@inline function same_axes(first_column, rest...)
    reduce_unrolled(&, true, partial_map(compare_axes, axes(first_column), rest)...)
end
@inline function same_axes()
    true
end

@inline function Rows{Row, Dimension}(columns::Columns, the_names::Names) where {Row, Dimension, Columns, Names}
    @boundscheck if !same_axes(columns...)
        throw(DimensionMismatch("All arguments to `Rows` must have the same axes"))
    end
    Rows{Row, Dimension, Columns, Names}(columns, the_names)
end

@inline function get_model(columns)
    first(columns)
end
@inline function get_model(::Tuple{})
    1:0
end

@inline function parent(rows::Rows)
    get_model(rows.columns)
end

@inline function name_eltype(::Name, ::Column) where {Name, Column}
    Tuple{Name, eltype(Column)}
end

@inline function Rows(columns, the_names)
    Rows{
        Tuple{map_unrolled(name_eltype, the_names, columns)...},
        ndims(get_model(columns))
    }(columns, the_names)
end

export Rows

"""
    Rows(named_columns)

Iterator over `rows` of a table. Always lazy. Use [`Peek`](@ref) to view.

```jldoctest Rows
julia> using LightQuery

julia> using Test: @inferred

julia> lazy = @name @inferred Rows((a = [1, 2], b = [1.0, 2.0]))
2-element Rows{Tuple{Tuple{Name{:a},Int64},Tuple{Name{:b},Float64}},1,Tuple{Array{Int64,1},Array{Float64,1}},Tuple{Name{:a},Name{:b}}}:
 ((`a`, 1), (`b`, 1.0))
 ((`a`, 2), (`b`, 2.0))

julia> @inferred collect(lazy)
2-element Array{Tuple{Tuple{Name{:a},Int64},Tuple{Name{:b},Float64}},1}:
 ((`a`, 1), (`b`, 1.0))
 ((`a`, 2), (`b`, 2.0))

julia> @name @inferred Rows((a = [1, 2],))
2-element Rows{Tuple{Tuple{Name{:a},Int64}},1,Tuple{Array{Int64,1}},Tuple{Name{:a}}}:
 ((`a`, 1),)
 ((`a`, 2),)
```

All arguments to Rows must have the same axes. Use `@inbounds` to override the
check.

```jldoctest Rows
julia> result = @name Rows((a = 1:2, b = 1:3))
ERROR: DimensionMismatch("All arguments to `Rows` must have the same axes")
[...]
```
"""
@inline function Rows(named_columns)
    Rows(map_unrolled(value, named_columns), map_unrolled(key, named_columns))
end

@inline function to_columns(rows::Rows)
    map_unrolled(tuple, rows.names, rows.columns)
end

@inline function axes(rows::Rows, dimensions...)
    axes(parent(rows), dimensions...)
end
@inline function size(rows::Rows, dimensions...)
    size(parent(rows), dimensions...)
end

@inline function column_getindex((columns, an_index), name)
    name, name(columns)[an_index...]
end
@inline function getindex(rows::Rows, an_index::Int...)
    partial_map(
        column_getindex,
        (to_columns(rows), an_index),
        rows.names
    )
end

@inline function column_setindex!((columns, an_index), (name, value))
    name(columns)[an_index...] = value
    nothing
end
@inline function setindex!(rows::Rows, row, an_index::Int...)
    partial_map(column_setindex!, (to_columns(rows), an_index), row)
    nothing
end

@inline function column_push!(columns, (name, value))
    push!(name(columns), value)
    nothing
end
@inline function push!(rows::Rows, row)
    partial_map(column_push!, to_columns(rows), row)
    nothing
end

@inline function val_fieldtypes(something)
    ()
end
@pure function val_fieldtypes(a_type::DataType)
    if a_type.abstract || (a_type.name == Tuple.name && isvatuple(a_type))
        ()
    else
        map_unrolled(Val, (a_type.types...,))
    end
end

@inline function decompose_named_type(::Val{Tuple{Name{name}, Value}}) where {name, Value}
    Name{name}(), Val{Value}()
end
@inline function decompose_named_type(type)
    missing
end

@inline function decompose_row_type(Row)
    filter_unrolled(
        !ismissing,
        map_unrolled(decompose_named_type, val_fieldtypes(Row))...
    )
end

@inline function similar_val(model, ::Val{Value}, dimensions) where {Value}
    similar(model, Value, dimensions)
end
@inline function similar_column((model, dimensions), (name, val_Value))
    name, similar_val(model, val_Value, dimensions)
end
@inline function similar(rows::Rows, ::Type{ARow}, dimensions::Dims) where {ARow}
    @inbounds Rows(partial_map(
        similar_column,
        (parent(rows), dimensions),
        decompose_row_type(ARow)
    ))
end

@inline function empty(column::Rows{OldRow}, ::Type{NewRow} = OldRow) where {OldRow, NewRow}
    similar(column, NewRow)
end

@inline function widen_column(::HasLength, new_length, an_index, name, column::AbstractArray{Element}, item::Item) where {Element, Item <: Element}
    @inbounds column[an_index] = item
    name, column
end
@inline function widen_column(::HasLength, new_length, an_index, name, column::AbstractArray, item)
    name, setindex_widen_up_to(column, item, an_index)
end
@inline function widen_column(::SizeUnknown, new_length, an_index, name, column::AbstractArray{Element}, item::Item) where {Element, Item <: Element}
    push!(column, item)
    name, column
end
@inline function widen_column(::SizeUnknown, new_length, an_index, name, column::AbstractArray, item)
    name, push_widen(column, item)
end
@inline function widen_column(iterator_size, new_length, an_index, name, ::Missing, item::Item) where {Item}
    new_column = Array{Union{Missing, Item}}(missing, new_length)
    @inbounds new_column[an_index] = item
    name, new_column
end
@inline function widen_column(iterator_size, new_length, an_index, name, ::Missing, ::Missing)
    name, Array{Missing}(missing, new_length)
end
@inline function widen_column_clumped(fixeds, variables)
    widen_column(fixeds..., variables...)
end

@inline function get_new_length(::SizeUnknown, rows, an_index)
    an_index
end
@inline function get_new_length(::HasLength, rows, an_index)
    length(LinearIndices(rows))
end

@inline function lone_column((name, column))
    name, column, missing
end
@inline function lone_value((name, value))
    name, missing, value
end
@inline function column_value((column_name, column), (value_name, value))
    column_name, column, value
end

@inline function widen_named(iterator_size, rows, row, an_index = length(rows) + 1)
    named_columns = to_columns(rows)
    column_names = map_unrolled(key, named_columns)
    value_names = map_unrolled(key, row)
    just_column_names = diff_unrolled(value_names, column_names...)
    in_both_names = diff_unrolled(just_column_names, column_names...)
    @inbounds Rows(partial_map(
        widen_column_clumped,
        (
            iterator_size,
            get_new_length(iterator_size, rows, an_index),
            an_index
        ),
        (
            map_unrolled(lone_column, just_column_names(named_columns))...,
            map_unrolled(
                column_value,
                in_both_names(named_columns),
                in_both_names(row)
            )...,
            map_unrolled(
                lone_value,
                diff_unrolled(column_names, value_names...)(row)
            )...
        )
    ))
end

@inline function push_widen(rows::Rows, row)
    widen_named(SizeUnknown(), rows, row)
end
@inline function setindex_widen_up_to(rows::Rows, row, an_index)
    widen_named(HasLength(), rows, row, an_index)
end

"""
    make_columns(rows)

Collect into columns.

```jldoctest make_columns
julia> using LightQuery

julia> using Test: @inferred

julia> stable(x) = @name (a = x, b = x + 0.0, c = x, d = x + 0.0, e = x, f = x + 0.0);

julia> @inferred make_columns(over(1:4, stable))
((`a`, [1, 2, 3, 4]), (`b`, [1.0, 2.0, 3.0, 4.0]), (`c`, [1, 2, 3, 4]), (`d`, [1.0, 2.0, 3.0, 4.0]), (`e`, [1, 2, 3, 4]), (`f`, [1.0, 2.0, 3.0, 4.0]))

julia> unstable(x) =
            @name if x <= 2
                (a = missing, b = string(x), d = string(x))
            else
                (a = x, b = missing, c = x)
            end;

julia> make_columns(over(1:4, unstable))
((`d`, Union{Missing, String}["1", "2", missing, missing]), (`a`, Union{Missing, Int64}[missing, missing, 3, 4]), (`b`, Union{Missing, String}["1", "2", missing, missing]), (`c`, Union{Missing, Int64}[missing, missing, 3, 4]))

julia> make_columns(when(over(1:4, unstable), row -> true))
((`d`, Union{Missing, String}["1", "2", missing, missing]), (`a`, Union{Missing, Int64}[missing, missing, 3, 4]), (`b`, Union{Missing, String}["1", "2", missing, missing]), (`c`, Union{Missing, Int64}[missing, missing, 3, 4]))
```
"""
function make_columns(rows)
    to_columns(_collect(
        (@inbounds Rows(())),
        rows,
        IteratorEltype(rows),
        IteratorSize(rows)
    ))
end

export make_columns


text:
abstract type LinearModel end

"""
    ESNtrain(esn::AbstractReservoirComputer, beta::Float64[, train_data])

Return the trained output layer using Ridge Regression.
"""
function ESNtrain(esn::AbstractReservoirComputer, beta::Float64; train_data::AbstractArray{Float64} = esn.train_data)

    states_new = nla(esn.nla_type, esn.states)
    W_out = (train_data*states_new')*inv(add_reg(states_new*states_new', beta))

    return W_out
end

function add_reg(X::AbstractArray{Float64}, beta::Float64)
    n = size(X, 1)
    for i=1:n
        X[i,i] += beta
    end
    return X
end

#MLJ Ridge
"""
    Ridge(lambda, solver::MLJLinearModels.Solver)

Return a LinearModel object for the training of the model using Ridge Regression with a MLJLinearModels method.
"""
struct Ridge{T<: AbstractFloat} <: LinearModel
    lambda::T
    solver::MLJLinearModels.Solver
end
#Ridge(lambda, solver) = Ridge{Float64}(lambda, solver)
"""
    ESNtrain(lm::LinearModel, esn::AbstractReservoirComputer[, train_data])

Return the trained output layer using an MLJLinearModels method built into a LinearSolver struct.
"""
ESNtrain(ridge::Ridge{T}, esn::AbstractReservoirComputer; train_data::AbstractArray{Float64} = esn.train_data) where T<: AbstractFloat = _ridge(esn, ridge; train_data = esn.train_data)

function _ridge(esn::AbstractReservoirComputer, ridge::Ridge; train_data::AbstractArray{Float64} = esn.train_data)

    states_new = nla(esn.nla_type, esn.states)
    W_out = zeros(Float64, size(train_data, 1), size(states_new, 1))
    for i=1:size(train_data, 1)
        r = RidgeRegression(lambda = ridge.lambda, fit_intercept = false)
        W_out[i,:] = MLJLinearModels.fit(r, states_new', train_data[i,:], solver = ridge.solver)
    end

    return W_out
end

#MLJ Lasso
"""
    Lasso(lambda, solver::MLJLinearModels.Solver)

Return a LinearModel object for the training of the model using Lasso with a MLJLinearModels method.
"""
struct Lasso{T<: AbstractFloat} <: LinearModel
    lambda::T
    solver::MLJLinearModels.Solver
end
#Lasso(lambda, solver) = Lasso{Float64}(lambda, solver)
ESNtrain(lasso::Lasso{T}, esn::AbstractReservoirComputer; train_data::AbstractArray{Float64} = esn.train_data) where T<: AbstractFloat = _lasso(esn, lasso; train_data = esn.train_data)

function _lasso(esn::AbstractReservoirComputer, lasso::Lasso; train_data::AbstractArray{Float64} = esn.train_data)

    states_new = nla(esn.nla_type, esn.states)
    W_out = zeros(Float64, size(train_data, 1), size(states_new, 1))
    for i=1:size(train_data, 1)
        l = LassoRegression(lambda = lasso.lambda, fit_intercept = false)
        W_out[i,:] = MLJLinearModels.fit(l, states_new', train_data[i,:], solver = lasso.solver)
    end

    return W_out
end

#MLJ ElastNet
"""
    ElastNet(lambda, gamma, solver::MLJLinearModels.Solver)

Return a LinearModel object for the training of the model using Elastic Net with a MLJLinearModels method.
"""
struct ElastNet{T<: AbstractFloat} <: LinearModel
    lambda::T
    gamma::T
    solver::MLJLinearModels.Solver
end
#ElastNet(lambda, gamma, solver) = ElastNet{Float64}(lambda, gamma, solver)
ESNtrain(elastnet::ElastNet{T}, esn::AbstractReservoirComputer; train_data::AbstractArray{Float64} = esn.train_data) where T<: AbstractFloat = _elastnet(esn, elastnet; train_data = esn.train_data)

function _elastnet(esn::AbstractReservoirComputer, elastnet::ElastNet; train_data::AbstractArray{Float64} = esn.train_data)

    states_new = nla(esn.nla_type, esn.states)
    W_out = zeros(Float64, size(train_data, 1), size(states_new, 1))
    for i=1:size(train_data, 1)
        en = ElasticNetRegression(lambda = elastnet.lambda, gamma = elastnet.gamma, fit_intercept = false)
        W_out[i,:] = MLJLinearModels.fit(en, states_new', train_data[i,:], solver = elastnet.solver)
    end

    return W_out
end

#MLJ Huber -> RobustHuber avoids conflict
"""
    ElastNet(delta, lambda, gamma, solver::MLJLinearModels.Solver)

Return a LinearModel object for the training of the model using the Huber function with a MLJLinearModels method.
"""
struct RobustHuber{T<: AbstractFloat} <: LinearModel
    delta::T
    lambda::T
    gamma::T
    solver::MLJLinearModels.Solver
end
#RobustHuber(delta, lambda, gamma, solver) = RobustHuber{Float64}(delta, lambda, gamma, solver)
ESNtrain(huber::RobustHuber{T}, esn::AbstractReservoirComputer; train_data::AbstractArray{Float64} = esn.train_data) where T<: AbstractFloat = _huber(esn, huber; train_data = esn.train_data)

function _huber(esn::AbstractReservoirComputer, huber::RobustHuber; train_data::AbstractArray{Float64} = esn.train_data)

    states_new = nla(esn.nla_type, esn.states)
    W_out = zeros(Float64, size(train_data, 1), size(states_new, 1))
    for i=1:size(train_data, 1)
        h = HuberRegression(delta = huber.delta,
            lambda = huber.lambda,
            gamma = huber.gamma,
            fit_intercept = false)
        W_out[i,:] = MLJLinearModels.fit(h, states_new', train_data[i,:], solver = huber.solver)
    end

    return W_out
end

"""
    HESNtrain(ridge::Ridge{T}, hesn::AbstractReservoirComputer; train_data::AbstractArray{Float64} = hesn.train_data)

Return the trained HESN output layer using an MLJLinearModels method built into a LinearSolver struct.
"""
HESNtrain(ridge::Ridge{T}, hesn::AbstractReservoirComputer; train_data::AbstractArray{Float64} = hesn.train_data) where T<: AbstractFloat = _PIridge(hesn, ridge; train_data = hesn.train_data)

function _PIridge(hesn::AbstractReservoirComputer, ridge::Ridge; train_data::AbstractArray{Float64} = hesn.train_data)

    states_new = nla(hesn.nla_type, hesn.states)
    states_new = vcat(states_new, hesn.physics_model_data[:, 2:end])
    W_out = zeros(Float64, size(train_data, 1), size(states_new, 1))
    for i=1:size(train_data, 1)
        r = RidgeRegression(lambda = ridge.lambda, fit_intercept = false)
        W_out[i,:] = MLJLinearModels.fit(r, states_new', train_data[i,:], solver = ridge.solver)
    end

    return W_out
end


text:
# This file is part of GenericSchur.jl, released under the MIT "Expat" license
# Portions derived from LAPACK, see below.

# invariant subspace computations using complex Schur decompositions

import LinearAlgebra: _ordschur!, _ordschur

function LinearAlgebra._ordschur(T::StridedMatrix{Ty},
                    Z::StridedMatrix{Ty},
                   select::Union{Vector{Bool},BitVector}) where Ty <: Complex
    _ordschur!(copy(T), copy(Z), select)
end

function LinearAlgebra._ordschur!(T::StridedMatrix{Ty},
                    Z::StridedMatrix{Ty},
                    select::Union{Vector{Bool},BitVector}) where Ty <: Complex
    # suppress most checks since this is an internal function expecting
    # components of a Schur object
    n = size(T,1)
    ks = 0
    for k=1:n
        if select[k]
            ks += 1
            if k != ks
                _trexchange!(T,Z,k,ks)
            end
        end
    end
    triu!(T)
    T,Z,diag(T)
end

"""
reorder `T` by a unitary similarity transformation so that `T[iold,iold]`
is moved to `T[inew,inew]`. Also right-multiply `Z` by the same transformation.
"""
function _trexchange!(T,Z,iold,inew)
    # this is Algorithm 7.6.1 in Golub & van Loan 4th ed.
    n = size(T,1)
    if (n < 1) || (iold == inew)
        return
    end
    if iold < inew
        krange = iold:inew-1
    else
        krange = iold-1:-1:inew
    end
    for k in krange
        Tkk = T[k,k]
        Tpp = T[k+1,k+1]
        G,_ = givens(T[k,k+1], Tpp - Tkk, k, k+1)
        if k+1 <= n
            lmul!(G,T)
        end
        rmul!(T,adjoint(G))
        isempty(Z) || rmul!(Z,adjoint(G))
    end
end

# eigvalscond and subspacesep are derived from LAPACK's ztrsen.
# LAPACK is released under a BSD license, and is
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.

"""
    eigvalscond(S::Schur,nsub::Integer) => Real

Estimate the reciprocal of the condition number of the `nsub` leading
eigenvalues of `S`. (Use `ordschur` to move a subspace of interest
to the front of `S`.)

See the LAPACK User's Guide for details of interpretation.
"""
function eigvalscond(S::Schur{Ty},nsub) where Ty
    n = size(S.T,1)
    # solve T₁₁ R - R T₂₂ = σ T₁₂
    R = S.T[1:nsub,nsub+1:n] # copy
    R, scale = trsylvester!(view(S.T,1:nsub,1:nsub),view(S.T,nsub+1:n,nsub+1:n),
                    R)
    rnorm = norm(R) # Frobenius norm, as desired
    s = (rnorm > 0) ?
        scale / (sqrt(scale*scale / rnorm + rnorm ) * sqrt(rnorm)) :
        one(real(Ty))
    s
end

"""
    subspacesep(S::Schur,nsub::Integer) => Real

Estimate the reciprocal condition of the separation angle for the
invariant subspace corresponding
to the leading block of size `nsub` of a Schur decomposition.
(Use `ordschur` to move a subspace of interest
to the front of `S`.)

See the LAPACK User's Guide for details of interpretation.
"""
function subspacesep(S::Schur{Ty}, nsub) where Ty
    n = size(S.T,1)
    scale = one(real(Ty))
    function f(X)
        # solve T₁₁ R - R T₂₂ = σ X
        R, s = trsylvester!(view(S.T,1:nsub,1:nsub),
                            view(S.T,nsub+1:n,nsub+1:n),
                            reshape(X,nsub,n-nsub))
        scale = s
        R
    end
    function fH(X)
        # solve T₁₁ᴴ R - R T₂₂ᴴ = σ X
        R, s = adjtrsylvester!(view(S.T,1:nsub,1:nsub),
                               view(S.T,nsub+1:n,nsub+1:n),
                               reshape(X,nsub,n-nsub))
        scale = s
        R
    end
    est = norm1est!(f,fH,zeros(Ty,nsub*(n-nsub)))
    return scale / est
end

function LinearAlgebra._ordschur(S::StridedMatrix{Ty}, T::StridedMatrix{Ty},
                                 Q::StridedMatrix{Ty}, Z::StridedMatrix{Ty},
                                 select::Union{Vector{Bool},BitVector}
                                 ) where Ty <: Complex{Tr} where Tr <: AbstractFloat
    _ordschur!(copy(S), copy(T), copy(Q), copy(Z), select)
end

function LinearAlgebra._ordschur!(S::StridedMatrix{Ty}, T::StridedMatrix{Ty},
                    Q::StridedMatrix{Ty}, Z::StridedMatrix{Ty},
                    select::Union{Vector{Bool},BitVector}
                    ) where Ty <: Complex{Tr} where Tr <: AbstractFloat
    n = size(S,1)
    ks = 0
    for k=1:n
        if select[k]
            ks += 1
            if k != ks
                _trexchange!(S,T,Q,Z,k,ks)
            end
        end
    end
    triu!(S)
    triu!(T)
    S,T,diag(S),diag(T),Q,Z
end

"""
    IllConditionException

Exception thrown when argument matrix or matrices are too ill-conditioned
for the requested operation. The `index` field may indicate the block
where near-singularity was detected.
"""
struct IllConditionException <: Exception
    index::Integer
end

function _trexchange!(A,B,Q,Z,iold,inew)
    if iold < inew
        icur = iold
        while icur < ilast
            _trexch2!(A,B,Q,Z,icur)
            icur += 1
        end
        icur -= 1
    else
        icur =  iold-1
        while icur >= inew
            _trexch2!(A,B,Q,Z,icur)
            icur -= 1
        end
        icur += 1
    end
    nothing
end

# swap adjacent 1x1 blocks in UT pair (A,B) via unitary equivalence transform
# apply same transform to associated Q,Z
function _trexch2!(A::AbstractMatrix{Ty},B,Q,Z,j1) where Ty
    S = A[j1:j1+1,j1:j1+1]
    T = B[j1:j1+1,j1:j1+1]
    # compute threshold for acceptance
    Tr = real(Ty)
    smlnum = safemin(Tr) / eps(Tr)
    scale = zero(Tr)
    sumsq = one(Tr)
    W = hcat(S,T)
    scale, sumsq = _ssq(vec(W), scale, sumsq)
    sa = scale * sqrt(sumsq)
    thresh = max(Tr(20)*eps(Tr)*sa, smlnum)

    # compute Givens rotations which would swap blocks
    # and tentatively do it
    f = S[2,2] * T[1,1] - T[2,2] * S[1,1]
    g = S[2,2] * T[1,2] - T[2,2] * S[1,2]
    sa = abs(S[2,2])
    sb = abs(T[2,2])
    cz, sz, _ = givensAlgorithm(g,f)
    Gz2 = Givens(1,2,Ty(cz),-sz)
    rmul!(S,adjoint(Gz2))
    rmul!(T,adjoint(Gz2))
    if sa >= sb
        cq, sq, _ = givensAlgorithm(S[1,1],S[2,1])
    else
        cq, sq, _ = givensAlgorithm(T[1,1],T[2,1])
    end
    Gq2 = Givens(1,2,Ty(cq),sq)
    lmul!(Gq2,S)
    lmul!(Gq2,T)
    # weak stability test: subdiags <= O( ϵ norm((S,T),"F"))
    ws = abs(S[2,1]) + abs(T[2,1])
    weak_ok = ws <= thresh

    if !weak_ok
        throw(IllConditionException(j1))
    end

    if false
        # FIXME: not yet checked, needs example

        # Strong stability test
        # Supposedly equiv. to
        # |[A-Qᴴ*S*Z, B-Qᴴ*T*Z]| <= O(1) ϵ |[A,B]| in F-norm
        # which appears rather strange
        Ss = copy(S)
        Ts = copy(T)
        Gz2 = Givens(1,2,Ty(cz),sz)
        rmul!(Ss,adjoint(Gz2))
        rmul!(Ts,adjoint(Gz2))
        Gq2 = Givens(1,2,Ty(cq),-sq)
        lmul!(Gq2,Ss)
        lmul!(Gq2,Ts)
        Ss .-= A[j1:j1+1,j1:j1+1]
        Ts .-= B[j1:j1+1,j1:j1+1]
        scale = zero(Tr)
        sumsq = one(Tr)
        W = hcat(Ss,Ts)
        scale, sumsq = _ssq(vec(W), scale, sumsq)
        ss = scale * sqrt(sumsq)
        strong_ok = ss <= thresh
        if !strong_ok
            throw(IllConditionException(j1))
        end
    end

    Gz = Givens(j1,j1+1,Ty(cz),-sz)
    rmul!(A,adjoint(Gz))
    rmul!(B,adjoint(Gz))
    Gq = Givens(j1,j1+1,Ty(cq),sq)
    lmul!(Gq,A)
    lmul!(Gq,B)
    A[j1+1,j1] = zero(Ty)
    B[j1+1,j1] = zero(Ty)
    rmul!(Z,adjoint(Gz))
    rmul!(Q,adjoint(Gq))
end


text:
# Autogenerated wrapper script for libsndfile_jll for aarch64-apple-darwin
export libsndfile, sndfile_cmp, sndfile_concat, sndfile_convert, sndfile_deinterleave, sndfile_info, sndfile_interleave, sndfile_metadata_get, sndfile_metadata_set, sndfile_play, sndfile_salvage

using FLAC_jll
using Ogg_jll
using libvorbis_jll
using alsa_jll
JLLWrappers.@generate_wrapper_header("libsndfile")
JLLWrappers.@declare_library_product(libsndfile, "@rpath/libsndfile.1.dylib")
JLLWrappers.@declare_executable_product(sndfile_cmp)
JLLWrappers.@declare_executable_product(sndfile_concat)
JLLWrappers.@declare_executable_product(sndfile_convert)
JLLWrappers.@declare_executable_product(sndfile_deinterleave)
JLLWrappers.@declare_executable_product(sndfile_info)
JLLWrappers.@declare_executable_product(sndfile_interleave)
JLLWrappers.@declare_executable_product(sndfile_metadata_get)
JLLWrappers.@declare_executable_product(sndfile_metadata_set)
JLLWrappers.@declare_executable_product(sndfile_play)
JLLWrappers.@declare_executable_product(sndfile_salvage)
function __init__()
    JLLWrappers.@generate_init_header(FLAC_jll, Ogg_jll, libvorbis_jll, alsa_jll)
    JLLWrappers.@init_library_product(
        libsndfile,
        "lib/libsndfile.1.0.31.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        sndfile_cmp,
        "bin/sndfile-cmp",
    )

    JLLWrappers.@init_executable_product(
        sndfile_concat,
        "bin/sndfile-concat",
    )

    JLLWrappers.@init_executable_product(
        sndfile_convert,
        "bin/sndfile-convert",
    )

    JLLWrappers.@init_executable_product(
        sndfile_deinterleave,
        "bin/sndfile-deinterleave",
    )

    JLLWrappers.@init_executable_product(
        sndfile_info,
        "bin/sndfile-info",
    )

    JLLWrappers.@init_executable_product(
        sndfile_interleave,
        "bin/sndfile-interleave",
    )

    JLLWrappers.@init_executable_product(
        sndfile_metadata_get,
        "bin/sndfile-metadata-get",
    )

    JLLWrappers.@init_executable_product(
        sndfile_metadata_set,
        "bin/sndfile-metadata-set",
    )

    JLLWrappers.@init_executable_product(
        sndfile_play,
        "bin/sndfile-play",
    )

    JLLWrappers.@init_executable_product(
        sndfile_salvage,
        "bin/sndfile-salvage",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# this file is the Julia code for the single-step approach in Christensenet al.(2021) when some individuals have no omics data
# the R code can be found in Christensenet al.(2021)
using DelimitedFiles,LinearAlgebra

pct=0.99  #also: 0.1,0.3,0.5,0.7,0.8,0.9,0.95,0.99
repAll=20 #rep1,...,rep20
accu_all=zeros(repAll)
for rep in 1:repAll
    #read missing index
    missing_index=Int.(vec(readdlm("/Users/tianjing/Box/omics/omics_missing_pattern_legarra/missing_index/missingIndex.pct$pct.rep$rep.txt")))
    nInd_noomics = length(missing_index)   #ind with no omics
    all_index=collect(1:1055)
    no_missing_index=setdiff(all_index,missing_index)
    new_order=[missing_index;no_missing_index]
    println("missing pct $pct, rep $rep, #missing ind $nInd_noomics")

    ###read omics
    M_all=readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data1/omics.txt")
    #reorder
    M_all=M_all[new_order,:]

    ###read y
    yy = vec(readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data1/y_no_ar.txt"))  #phenotype
    #reorder
    yy=yy[new_order]

    ### read GRM
    GRM=readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data1/GRM_geno.txt")
    #reorder
    GRM=GRM[new_order,new_order]

    nfeat = 1200          # omics feature
    nInd=1055            # total ind
    #############
    #set missing patter (the last 10 testing inds must have omics data)
    #first nInd_noomics ind without omics data

    #############
    nInd_omics = nInd-nInd_noomics  # ind with omics
    nInd_train = 955       # first 90 ind for training
    nInd_test  =  nInd-nInd_train # last 10 ind for testing
    println("#train ind $nInd_train, #test ind $nInd_test")

    M = M_all[(nInd_noomics+1):end,:]  #remove omics data for first nInd_noomics indviduals
    size(M)

    #standardization
    sdM=vec(std(M,dims=1))
    Mmean = M-ones(size(M,1))*mean(M,dims=1)
    for i in 1:nfeat
        Mmean[:,i]=Mmean[:,i]/sdM[i]/sqrt(nfeat)
    end
    M=Mmean
    size(M)
    sum(round.(mean(M,dims=1),digits=2))
    sum(round.(var(M,dims=1),digits=10))
    M[1:5,1:5]

    MMpsinv = inv(M*M'+I*0.0001)
    MMpsinv[1:5,1:5]

    σ2_am = 2
    h2_m = 0.61
    h2 = 0.47
    h2_ar = 0.2

    c2m = (h2-h2_ar)/h2_m
    s2y = h2_m*σ2_am/(h2-h2_ar)
    σ2_ar = h2_ar*s2y
    σ2_ϵ =(1-c2m-h2_ar)*s2y

    eta1 = σ2_ϵ/(c2m*s2y)
    ## eta1 also equals σ2_ϵ/σ2_am

    zeta = (1-h2_m)/h2_m

    GGzeta = h2_m*GRM

    GGzeta = GGzeta + I*(1-h2_m)
    ## note: M is standardised such that sum over omics features of mean.diag equals 1.
    ## Therfore a priori expectation of MM' is h2_m*GG+(1-h2)*Identity .

    Omega_inv = inv(GGzeta)
    GGzeta22_inv = inv(GGzeta[(nInd_noomics+1):end,(nInd_noomics+1):end])
    Omega_inv[(nInd_noomics+1):end,(nInd_noomics+1):end] = Omega_inv[(nInd_noomics+1):end,(nInd_noomics+1):end] + MMpsinv - GGzeta22_inv

    Ginv=inv(GRM+I*0.0001)
    ## X contains separate intercept for those with/without omics

    RHS = [sum(yy[1:nInd_noomics]) ; sum(yy[(nInd_noomics+1):nInd_train]) ; yy[1:nInd_train] ; zeros(nInd_test)]

    LHS11 = Diagonal([nInd_noomics,nInd_train-nInd_noomics])

    LHS12 = [ [ones(nInd_noomics);zeros(nInd_omics)]'
              [zeros(nInd_noomics) ; ones(nInd_train-nInd_noomics) ; zeros(nInd_test)]' ]

    LHS22 = [     [I
                   zeros(nInd_test,nInd_train)]  zeros(nInd,nInd_test)] + eta1*Omega_inv

    LHS = Matrix([LHS11  LHS12
                  LHS12' LHS22])

    SOL1 = inv(LHS)*RHS   #u=M*alpha, pluggin u into next MME
    # write(SOL1, file="SOL1_mm")

    ## SOL1 is phenotype prediction from omics M%*%alpha

    ## now we aim for the breeding values

    RHS = [sum(SOL1[3:nInd+2]) ; SOL1[3:nInd+2]]

    LHS11 = nInd
    LHS12 = ones(nInd)'
    LHS22 = I + zeta*Ginv

    LHS = [ LHS11  LHS12
            LHS12' LHS22]

    SOL2 = inv(LHS)*RHS
    SOL2[1:5]

    hat_a_m = SOL2[2:(nInd+1)]

    # write.table(hat.a.m, file="EBVss2", quote=FALSE, row.names=FALSE)

    ##prediction accuracy
    tbv=vec(readdlm("/Users/tianjing/Box/Omicsdata/andrew1055/data1/tbv_m.txt"))
    tbv=tbv[new_order]
    accuruacy_all  = cor(hat_a_m,tbv)
    accuruacy_test  = cor(hat_a_m[956:1055],tbv[956:1055])

    println("----------")
    @show pct,rep
    @show round(accuruacy_all,digits=3)
    @show round(accuruacy_test,digits=3)
    accu_all[rep]=accuruacy_test
end

accu_all
round.(accu_all,digits=3)

writedlm("/Users/tianjing/Box/omics/omics_missing_pattern_legarra/legarra_res/res_legarra_ptc$pct.txt",accu_all)

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
using Test
using CUDA
using KernelAbstractions

using ExaPF
import ExaPF: AutoDiff
const PS = ExaPF.PowerSystem
const LS = ExaPF.LinearSolvers

# Test quickstart guide in docs/src/quickstart.md
# If one test is broken, please update the documentation.

@testset "Documentation: quickstart" begin
    pglib_instance = "case1354.m"
    datafile = joinpath(dirname(@__FILE__), "..", "data", pglib_instance)

    # Short version
    polar = ExaPF.PolarForm(datafile, CPU())
    pf_algo = NewtonRaphson(; verbose=0, tol=1e-10)
    convergence = ExaPF.powerflow(polar, pf_algo)
    @test convergence.has_converged
    @test convergence.n_iterations == 5
    @test convergence.norm_residuals <= pf_algo.tol

    # Long version
    pf = PS.PowerNetwork(datafile)
    nbus = PS.get(pf, PS.NumberOfBuses())
    @test nbus == 1354

    pv_indexes = PS.get(pf, PS.PVIndexes())
    # Test only first PV index.
    @test pv_indexes[1] == 17

    # Build-up PolarForm object
    polar = ExaPF.PolarForm(pf, CPU())
    physical_state = get(polar, ExaPF.PhysicalState())
    ExaPF.init_buffer!(polar, physical_state)
    jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())

    linear_solver = LS.DirectSolver()
    convergence = ExaPF.powerflow(
        polar, jx, physical_state, pf_algo;
        linear_solver=linear_solver
    )

    @test convergence.has_converged
    @test convergence.n_iterations == 5
    @test convergence.norm_residuals <= pf_algo.tol

    # Reinit buffer
    ExaPF.init_buffer!(polar, physical_state)
    npartitions = 8
    jac = jx.J
    precond = LS.BlockJacobiPreconditioner(jac, npartitions, CPU())
    iterative_linear_solver = ExaPF.KrylovBICGSTAB(precond)
    @test isa(iterative_linear_solver, LS.AbstractIterativeLinearSolver)
    # Test default tolerance
    @test iterative_linear_solver.atol == 1e-10
    # Build powerflow algorithm
    pf_algo = NewtonRaphson(; verbose=0, tol=1e-7)

    convergence = ExaPF.powerflow(
        polar, jx, physical_state, pf_algo;
        linear_solver=iterative_linear_solver
    )

    @test convergence.has_converged
    @test convergence.n_iterations == 5
    @test convergence.norm_residuals <= pf_algo.tol

    if CUDA.has_cuda_gpu()
        polar_gpu = ExaPF.PolarForm(pf, CUDADevice())
        jx_gpu = AutoDiff.Jacobian(polar_gpu, ExaPF.power_balance, State())
        physical_state_gpu = get(polar_gpu, ExaPF.PhysicalState())
        ExaPF.init_buffer!(polar_gpu, physical_state_gpu)
        linear_solver = LS.DirectSolver()
        convergence = ExaPF.powerflow(
            polar_gpu, jx_gpu, physical_state_gpu, pf_algo;
            linear_solver=linear_solver
        )

        @test convergence.has_converged
        @test convergence.n_iterations == 5
        @test convergence.norm_residuals <= pf_algo.tol

        npartitions = 8
        jac = jx.J # we need to take the Jacobian on the CPU for partitioning!
        precond = LS.BlockJacobiPreconditioner(jac, npartitions, CUDADevice())

        # Reinit buffer
        ExaPF.init_buffer!(polar_gpu, physical_state_gpu)

        linear_solver = ExaPF.KrylovBICGSTAB(precond)
        pf_algo = NewtonRaphson(; verbose=0, tol=1e-7)
        convergence = ExaPF.powerflow(
            polar_gpu, jx_gpu, physical_state_gpu, pf_algo;
            linear_solver=linear_solver
        )

        @test convergence.has_converged
        @test convergence.n_iterations == 5
        @test convergence.norm_residuals <= pf_algo.tol
    end
end


text:
using Libdl
using Printf

const root_path = joinpath(@__DIR__, ".")
# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, ".", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Libxc not installed properly, run Pkg.build(\"Libxc\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

d = Dict{String, String}()
open(joinpath(root_path, "src", "functionals.jl"), "w") do f
    @show f
    for i in 1:1000
        out = ""
        try
            out = readchomp(`$xcinfo $i`)
        catch
            continue
        end
        ss = split(out, r"\n|\t")

        # Dict_keys are (name: family: kink: func_id: comment:)
        for s in ss[1:5]
            kv = split(s, ":")
            k, v = strip.(kv)
            d[k] = v
        end

        #
        out_str = @sprintf "const  %-26s = %6s # %s\n" uppercase(d["name"]) d["func_id"] d["comment"]
        @show
        write(f, out_str)
    end
end


text:
fname = :ARCSpectral

c = ARCTR.Combi(ARCTR.hessian_dense, ARCTR.PDataSpectral{Nothing}, ARCTR.solve_modelARCDiag, ARCTR.preprocessSpectral, ARCTR.decreaseFact, ARCTR.Tparam{Nothing}())

include("Template.jl")


text:
function test_diophantine_solve_diophantine_quadratic()
    soltype = Array{AbstractDiophantineSolutions{DiophantineSolutionXY{Int}},1}

    # linear
    sol = solve(diophantine_equation_quadratic_xy(cx=3, c0=-9))
    @test Set(sol) == Set([diophantine_onex_anyy(3)])
    @test length(sol) == 1
    @test typeof(sol) == soltype

    # simple hyperbolic
    sol = solve(diophantine_equation_quadratic_xy(cxy=4, cx=2, cy=8, c0=4))
    @test Set(sol) == Set([diophantine_onex_anyy(-2)])
    @test length(sol) == 1
    @test typeof(sol) == soltype

    sol = solve(diophantine_equation_quadratic_xy(cxy=4, cx=12, cy=8, c0=24))
    @test Set(sol) == Set([diophantine_onex_anyy(-2), diophantine_anyx_oney(-3)])
    @test length(sol) == 2
    @test typeof(sol) == soltype

    sol = solve(diophantine_equation_quadratic_xy(cxy=2, cx=5, cy=56, c0=7))
    @test Set(sol) == Set([diophantine_solutions((-27, 64), (-29, -69), (-21, 7), (-35, -12), (-9, 1), (-47, -6), (105, -2), (-161, -3))])
    @test length(sol) == 1
    @test typeof(sol) == soltype

    # elliptical
    sol = solve(diophantine_equation_quadratic_xy(cx²=42, cxy=8, cy²=15, cx=23, cy=17, c0=-4915))
    @test Set(sol) == Set([diophantine_solutions((-11, -1))])
    @test length(sol) == 1
    @test typeof(sol) == soltype

    # parabolic
    sol = solve(diophantine_equation_quadratic_xy(cx²=1, cxy=2, cy²=1, cx=2, cy=2, c0=1))
    @test Set(sol) == Set([diophantine_linearx_lineary(1, 0, -1, -1)])
    @test length(sol) == 1
    @test typeof(sol) == soltype

    sol = solve(diophantine_equation_quadratic_xy(cx²=1, cxy=2, cy²=1, cx=6, cy=6, c0=5))
    @test Set(sol) == Set([diophantine_linearx_lineary(1, 0, -1, -1), diophantine_linearx_lineary(1, 0, -1, -5)])
    @test length(sol) == 2
    @test typeof(sol) == soltype

    sol = solve(diophantine_equation_quadratic_xy(cx²=8, cxy=-24, cy²=18, cx=5, cy=7, c0=16))
    @test Set(sol) == Set([diophantine_quadraticx_quadraticy(-174, 17, -2, -116, 21, -2), diophantine_quadraticx_quadraticy(-174, 41, -4, -116, 37, -4)])
    @test length(sol) == 2
    @test typeof(sol) == soltype

    # hyperbolic (WIP)
    sol = solve(diophantine_equation_quadratic_xy(cx²=1, cxy=3, cy²=2))
    @test Set(sol) == Set([diophantine_solutions((0, 0)), diophantine_linearx_lineary(2, 0, -1, 0), diophantine_linearx_lineary(1, 0, -1, 0)])
    @test length(sol) == 3
    @test typeof(sol) == soltype

    sol = solve(diophantine_equation_quadratic_xy(cx²=1, cxy=3, cy²=2, c0=10))
    @test Set(sol) == Set([diophantine_solutions((-21, 11), (21, -11), (-12, 7), (12, -7), (-9, 7), (9, -7), (-12, 11), (12, -11))])
    @test length(sol) == 1
    @test typeof(sol) == soltype
end


text:
cool(x) = x + 1
cool(x, y) = x + y + 1

_second(t) = Base.tuple_type_head(Base.tuple_type_tail(t))

# a rule we define so we can test rules
dummy_identity(x) = x
@scalar_rule(dummy_identity(x), One())

@testset "rules" begin
    @testset "frule and rrule" begin
        @test frule(cool, 1) === nothing
        @test frule(cool, 1; iscool=true) === nothing
        @test rrule(cool, 1) === nothing
        @test rrule(cool, 1; iscool=true) === nothing

        ChainRulesCore.@scalar_rule(Main.cool(x), one(x))
        @test hasmethod(rrule, Tuple{typeof(cool),Number})
        ChainRulesCore.@scalar_rule(Main.cool(x::String), "wow such dfdx")
        @test hasmethod(rrule, Tuple{typeof(cool),String})
        # Ensure those are the *only* methods that have been defined
        cool_methods = Set(m.sig for m in methods(rrule) if _second(m.sig) == typeof(cool))
        only_methods = Set([Tuple{typeof(rrule),typeof(cool),Number},
                            Tuple{typeof(rrule),typeof(cool),String}])
        @test cool_methods == only_methods

        frx, fr = frule(cool, 1)
        @test frx == 2
        @test fr(1) == 1
        rrx, rr = rrule(cool, 1)
        @test rrx == 2
        @test rr(1) == 1
    end
    @testset "iterating and indexing rules" begin
        _, rule = frule(dummy_identity, 1)
        i = 0
        for r in rule
            @test r === rule
            i += 1
        end
        @test i == 1  # rules only iterate once, yielding themselves
        @test rule[1] == rule
        @test_throws BoundsError rule[2]
    end
end


text:
include("../src/sir.jl")
using PyPlot

function no_of_particles()
	K = 50
	σ_o = 0.1
	σ_d = 0.0
	Δ = 5
	Nth = 500
	obsfun = obs

	Npts = 10
	Np = 100

	trj_len = K*Δ
	mean_rmse = zeros(Npts)
	Np_arr = zeros(Int64, Npts)

	for i = 1:Npts
	   @show Np
	   mean_rmse[i] = assimilate(K, Np,
							σ_o, σ_d,
							Δ, Nth,
							obsfun)

	   Np_arr[i] = Np
	   Np = Np + 500
	end
	save("../data/mean_rmse_vs_Np_full.jld", "Np", Np_arr, "mean_rmse", mean_rmse)
end


text:
export ZabbixHostgroup, get_hostgroups

struct ZabbixHostgroup
  id::Int
  name::String
end
ZabbixHostgroup(dict) =
  ZabbixHostgroup(
    parse(Int, dict["groupid"]),
    dict["name"]
  )

showres(h::ZabbixHostgroup; offset=0) =
  println(repeat(" ", offset), "ZabbixHostgroup: ", h.id, " => ", h.name)

function get_hostgroups(z)
  _hostgroups = @checkerr make_request(z, "hostgroup.get")
  hostgroups = ZabbixHostgroup[]
  for h in _hostgroups
    push!(hostgroups, ZabbixHostgroup(h))
  end
  hostgroups
end


text:
using Documenter, FwiFlow
makedocs(sitename="FwiFlow", modules=[FwiFlow],
pages = Any[
    "index.md",
    "api.md",
    "Tutorial" => ["tutorials/fwi.md","tutorials/flow.md", "tutorials/timefrac.md"]

],
authors = "Dongzhuo Li and Kailai Xu")

deploydocs(
    repo = "github.com/lidongzh/FwiFlow.jl.git",
)