text:
precompile(Tuple{typeof(Base.stat), Base.Libc.RawFD})
precompile(Tuple{typeof(Base.isvalid_file_crc), Base.IOStream})
precompile(Tuple{typeof(Base.stale_cachefile), String, String})
precompile(Tuple{typeof(Base._include_from_serialized), String, Array{Any, 1}})
precompile(Tuple{typeof(Base._tryrequire_from_serialized), Base.PkgId, UInt64, String})
precompile(Tuple{typeof(Base._require_search_from_serialized), Base.PkgId, String})
precompile(Tuple{typeof(Base.create_expr_cache), String, String, Array{Base.Pair{Base.PkgId, UInt64}, 1}, Nothing})
precompile(Tuple{typeof(Base.compilecache), Base.PkgId, String})
precompile(Tuple{typeof(Base._tryrequire_from_serialized), Base.PkgId, UInt64, Nothing})
precompile(Tuple{typeof(Base._require_from_serialized), String})
precompile(Tuple{typeof(Base._require), Base.PkgId})
precompile(Tuple{typeof(Base.require), Base.PkgId})
precompile(Tuple{typeof(Base.require), Module, Symbol})
precompile(Tuple{typeof(Base.isassigned), Core.SimpleVector, Int64})
precompile(Tuple{typeof(Compat.Sys.__init__)})
precompile(Tuple{typeof(Blosc.check_deps)})
precompile(Tuple{typeof(Blosc.__init__)})
precompile(Tuple{typeof(HDF5.check_deps)})
precompile(Tuple{typeof(HDF5.blosc_set_local), Int64, Int64, Int64})
precompile(Tuple{typeof(HDF5.blosc_filter), UInt32, UInt64, Ptr{UInt32}, UInt64, Ptr{UInt64}, Ptr{Ptr{Nothing}}})
precompile(Tuple{typeof(HDF5.__init__)})
precompile(Tuple{typeof(Base.rehash!), Base.Dict{Tuple{Int64, Any, UInt64}, DataType}, Int64})
precompile(Tuple{typeof(Base.hash), Tuple{Nothing, UInt64}, UInt64})
precompile(Tuple{typeof(Base.hash), Tuple{Int64, Nothing, UInt64}, UInt64})
precompile(Tuple{typeof(Base.hashindex), Tuple{Int64, Nothing, UInt64}, Int64})
precompile(Tuple{typeof(Base.hash), Tuple{Int32, UInt64}, UInt64})
precompile(Tuple{typeof(Base.hash), Tuple{Int64, Int32, UInt64}, UInt64})
precompile(Tuple{typeof(Base.hashindex), Tuple{Int64, Int32, UInt64}, Int64})
precompile(Tuple{typeof(Base.rehash!), Base.Dict{String, Tuple{Any, Any, Int64}}, Int64})
precompile(Tuple{typeof(Base.rehash!), Base.Dict{Int64, Tuple{Function, Int64, Vararg{Int64, N} where N}}, Int64})
precompile(Tuple{typeof(Base.rehash!), Base.Dict{String, Type}, Int64})
precompile(Tuple{typeof(Base.rehash!), Base.Dict{Type, Nothing}, Int64})
precompile(Tuple{typeof(Base.active_project), Bool})
precompile(Tuple{typeof(Base.load_path_expand), String})
precompile(Tuple{typeof(Base.load_path)})
precompile(Tuple{typeof(Base.identify_package), String})
precompile(Tuple{typeof(Base.identify_package), Base.PkgId, String})
precompile(Tuple{typeof(Base.locate_package), Base.PkgId})
precompile(Tuple{typeof(Base.stale_cachefile), String, String})
precompile(Tuple{typeof(Base._tryrequire_from_serialized), Base.PkgId, UInt64, String})
precompile(Tuple{typeof(Base._require_search_from_serialized), Base.PkgId, String})
precompile(Tuple{typeof(Base.compilecache_path), Base.PkgId})
precompile(Tuple{typeof(Base.load_path_setup_code), Bool})
precompile(Tuple{typeof(Base.create_expr_cache), String, String, Array{Base.Pair{Base.PkgId, UInt64}, 1}, Nothing})
precompile(Tuple{typeof(Base.compilecache), Base.PkgId, String})
precompile(Tuple{typeof(Base._tryrequire_from_serialized), Base.PkgId, UInt64, Nothing})
precompile(Tuple{typeof(Base._require_from_serialized), String})
precompile(Tuple{typeof(Base._require), Base.PkgId})
precompile(Tuple{typeof(Base.require), Base.PkgId})
precompile(Tuple{typeof(Base.require), Module, Symbol})
precompile(Tuple{typeof(Requires.__init__)})
precompile(Tuple{typeof(Base.rehash!), Base.Dict{Base.PkgId, Array{Function, 1}}, Int64})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(Base.copyto!), Base.IndexLinear, Array{Function, 1}, Base.IndexLinear, Array{Any, 1}})
precompile(Tuple{typeof(Base.get!), getfield(Requires, Symbol("#1#2")), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(Base.throw_boundserror), Array{Array{Function, 1}, 1}, Tuple{Int64}})
precompile(Tuple{typeof(Base._delete!), Base.Dict{Base.PkgId, Array{Function, 1}}, Int64})
precompile(Tuple{typeof(Base._collect), Array{Function, 1}, Base.Generator{Array{Function, 1}, getfield(Requires, Symbol("#3#4"))}, Base.EltypeUnknown, Base.HasShape{1}})
precompile(Tuple{typeof(Requires.loadpkg), Base.PkgId})
precompile(Tuple{typeof(MPI.__init__)})
precompile(Tuple{typeof(Libdl.dlpath), String})
precompile(Tuple{typeof(Requires.listenpkg), getfield(MPI, Symbol("#41#44")), Base.PkgId})
precompile(Tuple{typeof(Base.Docs.docm), LineNumberNode, Module, Any, Any, Bool})
precompile(Tuple{typeof(Base.Docs.docm), LineNumberNode, Module, Any, Any})
precompile(Tuple{getfield(Base, Symbol("#@views")), LineNumberNode, Module, Any})
precompile(Tuple{typeof(Base.replace_ref_end_!), Float64, Nothing})
precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Int64})
precompile(Tuple{typeof(Base.getindex), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Int64})
precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!), Base.BitArray{1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Base.OneTo{Int64}, Int64, Int64})
precompile(Tuple{typeof(Base.copy), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Array{Any, 1}}}})
precompile(Tuple{typeof(Base.Broadcast.materialize), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(Base._views), Tuple{Array{Any, 1}}}})
precompile(Tuple{typeof(Base._views), Expr})
precompile(Tuple{Type{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Axes, F, Args} where Args<:Tuple where F where Axes}, typeof(Base._views), Tuple{Array{Any, 1}}})
precompile(Tuple{typeof(Base.Broadcast.instantiate), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(Base._views), Tuple{Array{Any, 1}}}})
precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!), Array{Symbol, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Base.OneTo{Int64}, Int64, Int64})
precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf), typeof(Base._views), Expr})
precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!), Array{Any, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Base.OneTo{Int64}, Int64, Int64})
precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!), Array{Expr, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, Base.OneTo{Int64}, Int64, Int64})
precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf), typeof(Base._views), Float64})
precompile(Tuple{typeof(Base.Docs.splitexpr), Expr})
precompile(Tuple{typeof(Base.Docs.bindingexpr), Any})
precompile(Tuple{typeof(Base.Docs.objectdoc), Any, Any, Any, Any, Any, Any})
precompile(Tuple{typeof(Base.Docs.docstr), Any, Any})
precompile(Tuple{typeof(Base.Docs._docstr), Core.SimpleVector, Any})
precompile(Tuple{typeof(Base.replace_ref_end_!), LineNumberNode, Nothing})
precompile(Tuple{typeof(Base.Broadcast._broadcast_getindex_evalf), typeof(Base._views), LineNumberNode})
precompile(Tuple{typeof(Base.Broadcast.restart_copyto_nonleaf!), Array{Any, 1}, Array{Symbol, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(Base._views), Tuple{Base.Broadcast.Extruded{Array{Any, 1}, Tuple{Bool}, Tuple{Int64}}}}, LineNumberNode, Int64, Base.OneTo{Int64}, Int64, Int64})
precompile(Tuple{typeof(Parsers.__init__)})
precompile(Tuple{typeof(Base.GMP.MPZ.set), Base.GMP.BigInt})
precompile(Tuple{typeof(Base.deepcopy), Base.GMP.BigInt})
precompile(Tuple{typeof(Base.Threads.resize_nthreads!), Array{Base.GMP.BigInt, 1}, Base.GMP.BigInt})
precompile(Tuple{typeof(Base.foreach), getfield(Parsers, Symbol("#33#38")), Array{Base.GMP.BigInt, 1}})
precompile(Tuple{typeof(Base.foreach), getfield(Parsers, Symbol("#34#39")), Array{Base.GMP.BigInt, 1}})
precompile(Tuple{typeof(Base.foreach), getfield(Parsers, Symbol("#35#40")), Array{Base.GMP.BigInt, 1}})
precompile(Tuple{typeof(Base.foreach), getfield(Parsers, Symbol("#36#41")), Array{Base.GMP.BigInt, 1}})
precompile(Tuple{typeof(Base.foreach), getfield(Parsers, Symbol("#37#42")), Array{Base.GMP.BigInt, 1}})
precompile(Tuple{getfield(Blosc, Symbol("#1#2"))})


text:
mutable struct UpdatableQR{T} <: Factorization{T}
    """
    Gives the qr factorization an (n, m) matrix as Q1*R1
    Q2 is such that Q := [Q1 Q2] is orthogonal and R is an (n, n) matrix where R1 "views into".
    """
    Q::Matrix{T}
    R::Matrix{T}
    n::Int
    m::Int

    Q1::SubArray{T, 2, Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}}, true}
    Q2::SubArray{T, 2, Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}}, true}
    R1::UpperTriangular{T, SubArray{T, 2, Matrix{T}, Tuple{UnitRange{Int},UnitRange{Int}}, false}}

    function UpdatableQR(A::AbstractMatrix{T}) where {T}
        n, m = size(A)
        @assert(m <= n, "Too many columns in the matrix.")

        F = qr(A)
        Q = F.Q*Matrix(I, n, n)
        R = zeros(T, n, n)
        R[1:m, 1:m] .= F.R

        new{T}(Q, R, n, m,
            view(Q, :, 1:m), view(Q, :, m+1:n),
            UpperTriangular(view(R, 1:m, 1:m)))
    end
end

function add_column!(F::UpdatableQR{T}, a::AbstractVector{T}) where {T}
    a1 = F.Q1'*a;
    a2 = F.Q2'*a;

    x = copy(a2)
    for i = length(x):-1:2
        G, r = givens(x[i-1], x[i], i-1, i)
        lmul!(G, x)
        lmul!(G, F.Q2')
    end

    F.R[1:F.m, F.m+1] .= a1
    F.R[F.m+1, F.m+1] = x[1]

    F.m += 1; update_views!(F)

    return a2
end

function add_column_householder!(F::UpdatableQR{T}, a::AbstractVector{T}) where {T}
    a1 = F.Q1'*a;
    a2 = F.Q2'*a;

    Z = qr(a2)
    LAPACK.gemqrt!('R','N', Z.factors, Z.T, F.Q2) # Q2 .= Q2*F.Q
    F.R[1:F.m, F.m+1] .= a1
    F.R[F.m+1, F.m+1] = Z.factors[1, 1]
    F.m += 1; update_views!(F)

    return Z
end

function remove_column!(F::UpdatableQR{T}, idx::Int) where {T}
    Q12 = view(F.Q, :, idx:F.m)
    R12 = view(F.R, idx:F.m, idx+1:F.m)

    for i in 1:size(R12, 1)-1
        G, r = givens(R12[i, i], R12[i + 1, i], i, i+1)
        lmul!(G, R12)
        rmul!(Q12, G')
    end

    for i in 1:F.m, j in idx:F.m-1
        F.R[i, j] = F.R[i, j+1]
    end
    F.R[:, F.m] .= zero(T)

    F.m -= 1; update_views!(F)

    return nothing
end

function update_views!(F::UpdatableQR{T}) where {T}
    F.R1 = UpperTriangular(view(F.R, 1:F.m, 1:F.m))
    F.Q1 = view(F.Q, :, 1:F.m)
    F.Q2 = view(F.Q, :, F.m+1:F.n)
end

mutable struct NullspaceHessianLDL{T}
    # Struct for the LDL factorization of the reduced Hessian
    # on the nullspace Z of the working constraints
    # U'*D*U = Z[:, end:-1:1]'*P*Z[:, end:-1:1]
    # The projected hessian can include at most one negative eigenvector

    n::Int # Dimension of the original space
    m::Int # Dimension of the nullspace
    artificial_constraints::Int

    P::Matrix{T}     # The full hessian
    # Z is the nullspace, i.e. QR.Q2 where QR is defined below
    Z::SubArray{T, 2, Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}}, true}
    U::UpperTriangular{T, SubArray{T, 2, Matrix{T}, Tuple{UnitRange{Int},UnitRange{Int}}, false}}
    D::Diagonal{T, SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}}
    QR::UpdatableQR{T}
    data::Matrix{T}  # That's where U is viewing into
    d::Vector{T}     # That's where D.diag views into
    indefinite_tolerance::T

    function NullspaceHessianLDL(P::Matrix{T}, A::AbstractMatrix{T}) where {T}
        @assert(size(A, 1) == size(P, 1) == size(P, 2), "Matrix dimensions do not match.")

        F = UpdatableQR(A)
        n = F.n
        m = F.n - F.m

        data = zeros(T, n, n)
        F.Q2 .= F.Q2[:, m:-1:1]
        WPW = F.Q2'*P*F.Q2; WPW .= (WPW .+ WPW')./2
        indefinite_tolerance = 1e-12

        C = cholesky(WPW, Val(true); check=false)
        if all(diag(C.U) .> indefinite_tolerance)
            F.Q2 .= F.Q2[:, C.p]
            artificial_constraints = 0
        else
            idx = findfirst(diag(C.U) .<= indefinite_tolerance)
            F.Q2 .= [F.Q2[:, C.p[idx:end]] F.Q2[:, C.p[1:idx-1]]]
            artificial_constraints = m + 1 - idx
        end
        m -= artificial_constraints
        U = view(data, 1:m, 1:m)

        F.Q2[:, artificial_constraints+1:end].= F.Q2[:, end:-1:artificial_constraints+1]
        Z = view(F.Q2, :, artificial_constraints+1:size(F.Q2, 2))
        U .= view(C.U, 1:m, 1:m)
        d = ones(T, n)
        D = Diagonal(view(d, 1:m))

        new{T}(n, m, artificial_constraints, P, Z, UpperTriangular(U), D, F, data, d, indefinite_tolerance)
    end

end

function add_constraint!(F::NullspaceHessianLDL{T}, a::AbstractVector{T}) where {T}
    a2 = add_column!(F.QR, a)
    if F.m == 1 # Nothing to write
        F.m -= 1; update_views!(F)
        return nothing
    end

    l = length(a2)
    for i = l:-1:2
        if l-i+2 <= F.m
            G, _ = givens(a2[i-1], a2[i], l-i+1, l-i+2)
            rmul!(F.U.data, G)
        end
        G, _ = givens(a2[i-1], a2[i], i-1, i)
        lmul!(G, a2)
    end
    F.m -= 1; update_views!(F)
    lmul!(sqrt.(F.D), F.U.data)
    hessenberg_to_triangular!(F.U.data)

    z = view(F.Z, :, 1)
    Pz = F.P*z
    if F.d[F.m + 1] <= -10
        # Recompute the last column of U
        l = reverse!(F.Z'*Pz)
        F.U[:, end] = F.U'\l
    end
    # Correct last element of U
    u1 = view(F.U, 1:F.m-1, F.m)
    d_new = dot(z, F.P*z) - dot(u1, u1)
    # Prevent F.U having zero columns
    F.U[end] = max(sqrt(abs(d_new)), F.indefinite_tolerance)

    # Scale matrices so that diag(F.U) = ones(..)
    F.D.diag .= one(T)
    F.D.diag[end] *= sign(d_new)

    return nothing
end

function remove_constraint!(F::NullspaceHessianLDL{T}, idx::Int) where{T}
    @assert(F.m == 0 || F.D[end] > F.indefinite_tolerance,
        "Constraints can be removed only when the reduced Hessian was already Positive Semidefinite.")

    if F.artificial_constraints == 0
        remove_column!(F.QR, idx)
    else
        idx != 0 && @warn("Ignoring non-zero index: Removing an artificial constraint.")
        F.artificial_constraints -= 1
    end

    z = view(F.QR.Q2, :, F.artificial_constraints + 1)
    Pz = F.P*z
    u = F.D\(F.U'\reverse!(F.Z'*Pz))
    d_new = dot(z, Pz) - dot(u, F.D*u)

    F.m += 1; update_views!(F)

    F.U[1:end-1, end] .= u
    F.U[end, end] = one(T)
    F.D[end] = d_new

    return nothing
end

function update_views!(F::NullspaceHessianLDL{T}) where {T}
    if size(F.U, 1) > F.m
        F.U.data[F.m+1, :] .= zero(T)
        F.U.data[:, F.m+1] .= zero(T)
    end
    F.U = UpperTriangular(view(F.data, 1:F.m, 1:F.m))
    F.D = Diagonal(view(F.d, 1:F.m))
    F.Z = view(F.QR.Q2, :, F.artificial_constraints+1:size(F.QR.Q2, 2))
end

function hessenberg_to_triangular!(A::AbstractMatrix{T}) where {T}
    n = size(A, 1)
    for i in 1:n-1
        G, _ = givens(A[i, i], A[i + 1, i], i, i+1)
        lmul!(G, A)
    end
    return A
end

mutable struct NullspaceHessian{T}
    n::Int # Dimension of the original space
    m::Int # Dimension of the nullspace

    P::Matrix{T}     # The full hessian
    # Z is the nullspace, i.e. QR.Q2 where QR is defined below
    Z::SubArray{T, 2, Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}}, true}
    ZPZ::SubArray{T, 2, Matrix{T}, Tuple{UnitRange{Int},UnitRange{Int}}, false}  # equal to Z'*P*Z
    QR::UpdatableQR{T}
    data::Matrix{T}  # That's where ZPZ is viewing into

    function NullspaceHessian{T}(P::Matrix{T}, A::Matrix{T}) where {T}
        @assert(size(A, 1) == size(P, 1) == size(P, 2), "Matrix dimensions do not match.")

        F = UpdatableQR(A)
        n = F.n
        m = F.n - F.m

        data = zeros(T, n, n)
        ZPZ = view(data, n-m+1:n, n-m+1:n)
        ZPZ .= F.Q2'*P*F.Q2; ZPZ .= (ZPZ .+ ZPZ')./2

        new{T}(n, m, P, F.Q2, ZPZ, F, data)
    end

    function NullspaceHessian{T}(P::Matrix{T}, F::UpdatableQR{T}) where {T}
        n = F.n
        m = F.n - F.m
        @assert(n == size(P, 1) == size(P,2), "Dimensions do not match.")

        data = zeros(T, n, n)
        ZPZ = view(data, n-m+1:n, n-m+1:n)
        ZPZ .= F.Q2'*P*F.Q2; ZPZ .= (ZPZ .+ ZPZ')./2

        new{T}(n, m, P, F.Q2, ZPZ, F, data)
    end
end

function add_constraint!(H::NullspaceHessian{T}, a::Vector{T}) where {T}
    Z = add_column_householder!(H.QR, a)

    LAPACK.gemqrt!('L','T',Z.factors,Z.T,H.ZPZ)
    LAPACK.gemqrt!('R','N',Z.factors,Z.T,H.ZPZ)
    # ToDo: Force symmetry? (i.e. H.ZPZ .= (H.ZPZ .+ H.ZPZ')./2)
    H.m -= 1; update_views!(H)
    # H.ZPZ .= (H.ZPZ .+ H.ZPZ')./2

    return nothing
end

function remove_constraint!(H::NullspaceHessian{T}, idx::Int) where{T}
    remove_column!(H.QR, idx)
    H.m += 1; update_views!(H)

    Pz = H.P*view(H.Z, :, 1)  # ToDo: avoid memory allocation
    mul!(view(H.ZPZ, 1, :), H.Z', Pz)
    for i = 2:H.m
        H.ZPZ[i, 1] = H.ZPZ[1, i]
    end

    return nothing
end

function update_views!(H::NullspaceHessian{T}) where {T}
    range = H.n-H.m+1:H.n
    H.ZPZ = view(H.data, range, range)
    H.Z = H.QR.Q2
end

ext: jl


text:
Base.display(doc::Document) = println("Doc: $(basename(doc._uri)) ", doc.code.index)
function Base.display(server::LanguageServerInstance)
    println("Root: ", server.workspaceFolders)
    for (uri, d) in server.documents
        display(d)
    end
end


text:
using HungarianAlg: hungarian, display

msize = 4
matrix = rand(0:200, (msize, msize))
println(string(join([string(matrix[a, :]) for a in 1:msize], "\n")))
println()

result = hungarian(matrix)
println(display(result))
println()

println(result.match)
println(result.revenues)
println(result.row_weights)
println(result.col_weights)
println(result.revenue_sum)

max_stars_repo_name: jbrightuniverse/hungarianalg-julia

max_issues_repo_name: jbrightuniverse/hungarianalg-julia

max_forks_repo_name: jbrightuniverse/hungarianalg-julia


text:
using Makie
using Observables

function cleanaxes()
    scene, layout = layoutscene()
    ax = layout[1, 1] = Axis(scene)
    axbox = pixelarea(ax.scene)[]
    lim = ax.finallimits[]
    e = events(ax.scene)
    return ax, axbox, lim, e
end

@testset "zoom Axis" begin
    ax, axbox, lim, e = cleanaxes()
    # Put the mouse in the center
    e.mouseposition[] = Tuple(axbox.origin + axbox.widths/2)
    # zoom in
    e.scroll[] = (0.0, -1.0)
    newlim = ax.finallimits[]
    @test newlim.widths ≈ 0.9 * lim.widths
    @test newlim.origin ≈ lim.origin + (lim.widths - newlim.widths)/2
    # zoom out restores original position
    e.scroll[] = (0.0, 1.0)
    newlim = ax.finallimits[]
    @test newlim.widths ≈ lim.widths
    @test all(abs.(newlim.origin - lim.origin) .< 1e-7*lim.widths)
    ax.finallimits[] = lim
    # Put mouse in corner
    e.mouseposition[] = Tuple(axbox.origin)
    # zoom in
    e.scroll[] = (0.0, -1.0)
    newlim = ax.finallimits[]
    @test newlim.widths ≈ 0.9 * lim.widths
    @test all(abs.(newlim.origin - lim.origin) .< 1e-7*lim.widths)
    # zoom out
    e.scroll[] = (0.0, 1.0)
    newlim = ax.finallimits[]
    @test newlim.widths ≈ lim.widths
    @test all(abs.(newlim.origin - lim.origin) .< 1e-7*lim.widths)
    ax.finallimits[] = lim

    # Zoom only x or y
    for (lockname, idx, zoomkey) in ((:xzoomlock, 1, :yzoomkey), (:yzoomlock, 2, :xzoomkey))
        ax, axbox, lim, e = cleanaxes()
        lock = getproperty(ax, lockname)
        @test !lock[]
        lock[] = true
        e.mouseposition[] = Tuple(axbox.origin + axbox.widths/2)
        e.scroll[] = (0.0, -1.0)
        newlim = ax.finallimits[]
        @test newlim.widths[idx] == lim.widths[idx]
        @test newlim.widths[3-idx] ≈ 0.9 * lim.widths[3-idx]
        @test newlim.origin[idx] == lim.origin[idx]
        @test newlim.origin[3-idx] ≈ lim.origin[3-idx] + (lim.widths[3-idx] - newlim.widths[3-idx])/2
        e.scroll[] = (0.0, 1.0)
        newlim = ax.finallimits[]
        @test newlim.widths ≈ lim.widths
        @test all(abs.(newlim.origin - lim.origin) .< 1e-7*lim.widths)
        ax.finallimits[] = lim
        lock[] = false
        # Simulate pressing the keys
        key = getproperty(ax, zoomkey)[]
        keypresses = getfield(e, Makie.button_key(key))[]
        @test isempty(keypresses)
        push!(keypresses, key)
        e.scroll[] = (0.0, -1.0)
        newlim = ax.finallimits[]
        @test newlim.widths[idx] == lim.widths[idx]
        @test newlim.widths[3-idx] ≈ 0.9 * lim.widths[3-idx]
        @test newlim.origin[idx] == lim.origin[idx]
        @test newlim.origin[3-idx] ≈ lim.origin[3-idx] + (lim.widths[3-idx] - newlim.widths[3-idx])/2
        e.scroll[] = (0.0, 1.0)
        newlim = ax.finallimits[]
        @test newlim.widths ≈ lim.widths
        @test all(abs.(newlim.origin - lim.origin) .< 1e-7*lim.widths)
    end

    # Rubber band selection
    scene, layout = layoutscene()
    ax = layout[1, 1] = Axis(scene)
    plot!(ax, [10, 15, 20])
    axbox = pixelarea(ax.scene)[]
    lim = ax.finallimits[]
    e = events(ax.scene)

    e.mouseposition[] = Tuple(axbox.origin)
    e.mousebuttons[] = Set([Mouse.left])
    e.mousedrag[] = Mouse.down
    e.mouseposition[] = Tuple(axbox.origin + axbox.widths ./ Vec2(2, 3))
    e.mousedrag[] = Mouse.pressed
    e.mousebuttons[] = Set{typeof(Mouse.left)}()
    e.mousedrag[] = Mouse.up
    newlim = ax.finallimits[]
    @test newlim.origin ≈ lim.origin
    @test newlim.widths ≈ lim.widths ./ Vec2(2, 3)
    # Ctrl-click to restore
    key = Makie.Keyboard.left_control
    keypresses = e.keyboardbuttons[]
    @test isempty(keypresses)
    push!(keypresses, key)
    buttons = e.mousebuttons[]
    @test isempty(buttons)
    e.mousebuttons[] = Set([Mouse.left])
    empty!(e.mousebuttons[])
    empty!(keypresses)
    newlim = ax.finallimits[]
    @test all(lim.origin .>= newlim.origin) && all(lim.origin+lim.widths .<= newlim.origin+newlim.widths)
end

@testset "pan Axis" begin
    ax, axbox, lim, e = cleanaxes()
    e.mouseposition[] = Tuple(axbox.origin + axbox.widths/2)
    e.scroll[] = (0.0, -1.0)
    newlim = ax.finallimits[]
    e.mouseposition[] = Tuple(axbox.origin)
    panbtn = ax.panbutton[]
    e.mousebuttons[] = Set([panbtn])
    e.mousedrag[] = Mouse.down
    e.mouseposition[] = Tuple(axbox.origin + axbox.widths/10)
    e.mousedrag[] = Mouse.pressed
    e.mousebuttons[] = Set{typeof(panbtn)}()
    e.mousedrag[] = Mouse.up
    panlim = ax.finallimits[]
    @test panlim.widths == newlim.widths
    @test (5/4)*panlim.origin ≈ -newlim.origin
end


text:
function load_and_init(
    vires::ViresConnection,
    #scenario::Union(String, ETree),
    scenario::Union{String, ETree},
    mode::Symbol;
    timeout::Float64=TIMEOUT_LONG_DEFAULT,
    ntries::Int=3
    )

    completed = false
    while ntries > 0 && !completed
        ntries -= 1
        write(vires.SCP, SCPMessage(get_xml_simctrl_load_init(scenario, mode)))
        completed = wait_for_packet_with_element(vires.SCP, "SimCtrl", "InitDone", timeout=timeout)
    end
    if !completed
        error("failed to init")
    end
end

function step_with_udp_checking_frameid(vires::ViresConnection, nsteps::Int, frame_start::Int=0)

    #=
    Make sure that UDP messages are received from frame_start to frame_start+nsteps-1
    =#

    # tic()
    stepcount = 0
    framecount = frame_start-1
    while stepcount < nsteps
        stepcount += 1
        pkg = read(vires.UDP, RDB_Package)
        if !isempty(pkg.entry_headers)
            if pkg.entry_headers[1].pkgId == RDB_PKG_ID_END_OF_FRAME
                framecount = pkg.header.frameNo
                # println("found frame: ", framecount, " / ", nsteps)
            end
        end
    end
    # println("done!", toc())
end
function step_with_udp_explicit_stepwidth(
    vires::ViresConnection,
    nsteps::Int;
    framestep::Real=0.04, # [s]
    frame_start::Int=0,
    simtime_start::Float64=0.0
    )

    #=
    Make sure that UDP messages are received from frame_start to frame_start+nsteps-1
    =#

    stepcount = 0
    Δt = convert(Cfloat, framestep)
    simTime = simtime_start
    frameNo = convert(UInt32, frame_start)

   # playerId:        1
   # steeringWheel:   +0.0000000000000000e+00
   # steeringSpeed:   +6.2826038338243961e-03
   # throttlePedal:   +0.0000000000000000e+00
   # brakePedal:      +0.0000000000000000e+00
   # clutchPedal:     +0.0000000000000000e+00
   # accelTgt:        +3.4999719355255365e-04
   # steeringTgt:     -8.7401109340135008e-05
   # curvatureTgt:    +0.0000000000000000e+00
   # steeringTorque:  +0.0000000000000000e+00
   # engineTorqueTgt: +0.0000000000000000e+00
   # speedTgt:        +0.0000000000000000e+00
   # gear:            4
   # sourceId:        0
   # spare0:          (0, 0)
   # validityFlags:   40e2 <flags, gear, tgt accel, tgt steer, steer speed>
   # flags:           0
   # mockupInput0:    0
   # mockupInput1:    0
   # mockupInput2:    0
   # spare:           0

    trigger = RDB_TRIGGER_t(Δt, frameNo, 0x00000000)

    # myDriver->playerId       = mOwnPlayerId;
    # myDriver->steeringWheel  = 9.0 * sin( 0.8 * simTime );
    # myDriver->throttlePedal  = 0.5 * ( 1 + sin( 0.8 * simTime ) );
    # myDriver->brakePedal     = 0.0;
    # myDriver->clutchPedal    = 0.0f;
    # myDriver->gear           = RDB_GEAR_BOX_POS_D;
    # myDriver->validityFlags  = RDB_DRIVER_INPUT_VALIDITY_STEERING_WHEEL | RDB_DRIVER_INPUT_VALIDITY_THROTTLE | RDB_DRIVER_INPUT_VALIDITY_BRAKE;

    # ctrl = RDB_DRIVER_CTRL_t(
        # convert(Uint32, 1), # playerId
        # convert(Cfloat, 4.0), # steeringWheel
        # convert(Cfloat, 0.0), # steeringSpeed::Cfloat
        # convert(Cfloat, 1.0), # throttlePedal::Cfloat
        # convert(Cfloat, 0.0), # brakePedal::Cfloat
        # convert(Cfloat, 0.0), # clutchPedal::Cfloat
        # convert(Cfloat, 0.0), # accelTgt::Cfloat
        # convert(Cfloat, 0.0), # steeringTgt::Cfloat
        # convert(Cdouble, 0.0), # curvatureTgt::Cdouble
        # convert(Cfloat, 0.0), # steeringTorque::Cfloat
        # convert(Cfloat, 0.0), # engineTorqueTgt::Cfloat
        # convert(Cfloat, 0.0), # speedTgt::Cfloat
        # convert(Uint8, RDB_GEAR_BOX_POS_D), # gear::Uint8
        # convert(Uint8, RDB_DRIVER_SOURCE_UNKNOWN), # sourceId::Uint8
        # (convert(Uint8, 0), convert(Uint8, 0)), # spare0::(Uint8, Uint8)
        # convert(Uint32, RDB_DRIVER_INPUT_VALIDITY_STEERING_WHEEL | RDB_DRIVER_INPUT_VALIDITY_THROTTLE |
        #                 RDB_DRIVER_INPUT_VALIDITY_BRAKE), # validityFlags::Uint32
        # convert(Uint32, RDB_DRIVER_FLAG_NONE), # flags::Uint32
        # convert(Uint32, 0), # mockupInput0::Uint32
        # convert(Uint32, 0), # mockupInput1::Uint32
        # convert(Uint32, 0), # mockupInput2::Uint32
        # convert(Uint32, 0), # spare::Uint32
        # )
    # start_of_frame = RDB_START_OF_FRAME_t()
    # end_of_frame = RDB_END_OF_FRAME_t()

    println("listening")

    while stepcount < nsteps
        stepcount += 1

        finished = false
        while !finished
            pkg = read(vires.UDP, RDB_Package)
            println(pkg)
            if !isempty(pkg.entry_headers)
                if pkg.entry_headers[1].pkgId == RDB_PKG_ID_START_OF_FRAME
                    frameNo = pkg.header.frameNo
                    simTime = pkg.header.simTime
                elseif pkg.entry_headers[1].pkgId == RDB_PKG_ID_END_OF_FRAME
                    finished = true
                # elseif pkg.entry_headers[1].pkgId == RDB_PKG_ID_DRIVER_CTRL
                #     println(STDOUT, bytes2hex(pkg.data_blocks[1]))
                    # buf = IOBuffer(pkg.data_blocks[1])
                    # ctrl2 = read(buf, RDB_DRIVER_CTRL_t)
                    # println(ctrl2)
                end
            end
        end

        # println("frameNo: ", frameNo)
        trigger.frameNo = frameNo + 1
        # write_udp_packet(vires.UDP, start_of_frame, frameNo, simTime)
        # write_bytehex_to_stdout((buf)->write_udp_packet(buf, ctrl, frameNo, simTime))
        # write_udp_packet(vires.UDP, ctrl, frameNo, simTime)
        write_udp_packet(vires.UDP, trigger, frameNo, simTime)
        # write_udp_packet(vires.UDP, end_of_frame, frameNo, simTime)
        # sleep(2.0)

        # ctrl.steeringWheel  = 9.0 * sin( 0.8 * simTime )
        # ctrl.throttlePedal  = 0.5 * ( 1 + sin( 0.8 * simTime ) )
        # ctrl.brakePedal     = 0.0
        # ctrl.clutchPedal    = 0.0
        # ctrl.gear           = RDB_GEAR_BOX_POS_D
        # ctrl.validityFlags  = RDB_DRIVER_INPUT_VALIDITY_STEERING_WHEEL | RDB_DRIVER_INPUT_VALIDITY_THROTTLE | RDB_DRIVER_INPUT_VALIDITY_BRAKE

        # write_udp_packet(vires.UDP, start_of_frame, frameNo, simTime)
        # write_udp_packet(vires.UDP, ctrl, convert(Uint32, frameNo+1), convert(Cfloat, simTime))
        # write_udp_packet(vires.UDP, trigger, frameNo, simTime)
        # write_udp_packet(vires.UDP, end_of_frame, frameNo, simTime)
    end
end
function step_with_udp_authorization(vires::ViresConnection, nsteps::Int)

    udp_sync = create_udp_entry(RDB_SYNC_t, vires.sync_mask)
    println("MASK: ", vires.sync_mask)
    println("UDP SYNC: ", udp_sync)

    tic()
    stepcount = 0
    while stepcount < nsteps
        stepcount += 1

        pkg = read(vires.UDP, RDB_Package)

        if !isempty(pkg.entry_headers)
            if pkg.entry_headers[1].pkgId == RDB_PKG_ID_START_OF_FRAME
                framecount = pkg.header.frameNo
                println("found frame: ", framecount, " / ", nsteps)

                println("received: ")
                println(pkg)

                write_udp_packet(vires.UDP, udp_sync, pkg.header.frameNo, pkg.header.simTime)

                # buf = IOBuffer()
                # write_udp_packet(buf, udp_sync, pkg.header.frameNo, pkg.header.simTime)
                # arr = takebuf_array(buf)
                # println(STDOUT, bytes2hex(arr))

                # write_udp_packet(buf, udp_sync, pkg.header.frameNo, pkg.header.simTime)
                # buf2 = IOBuffer(takebuf_string(buf))
                # pkg2 = read_special(buf2, RDB_Package)
                # println(STDOUT, pkg2)
            end
        end
    end
    toc()
end

function record_scenario_run(
    vires::ViresConnection,
    nsteps::Int,
    destpath::String,
    timeout::Float64=TIMEOUT_DEFAULT;
    remove_dat_files=true
    )

    dir, file = splitdir(destpath)
    dir = isempty(dir) ? "." : dir
    file_base = splitext(file)[1]
    file_dat = file_base * ".dat"

    write_and_wait_for_mirrored_message(vires.SCP,
                                        SCPMessage(get_xml_record_overwrite(dir, file_dat)),
                                        "record timed out",
                                        timeout=timeout)

    # create_camera_relative_to_ego(vires, dx=18.0, dz=100.0)

    # write_and_wait_for_mirrored_message(vires.SCP,
    #     SCPMessage("<Player name=\"Ego\" visible=\"true\"/>"),
    #     "make player visible timed out")
    # write_and_wait_for_mirrored_message(vires.SCP,
    #     SCPMessage("<Camera><PosEyepoint player=\"Ego\"/><ViewRelative dh=\"0.0\"  dp=\"0.0\" dr=\"0.0\"/></Camera>"),
    #     "set camera timed out")

    write_and_wait_for_mirrored_message(vires.SCP,
                                        SCPMessage("<SimCtrl><Start/></SimCtrl>"),
                                        "did not step",
                                        timeout=timeout)

    # write_and_wait_for_mirrored_message(vires.SCP,
    #                                     SCPMessage("<Player name=\"Ego\"><Driver ctrlLatLong=\"extern\"/></Player>"),
    #                                     "record timed out",
    #                                     timeout=timeout)

    # if vires.registered_to_authorize_framestep
    #     step_with_udp_authorization(vires, nsteps)
    # else
    #     step_with_udp_checking_frameid(vires, nsteps)
    # end
    step_with_udp_explicit_stepwidth(vires, nsteps)

    write_and_wait_for_mirrored_message(vires.SCP,
                                        SCPMessage(get_xml_record_stop()),
                                        "did not stop record",
                                        timeout=timeout)

    convert_recording_to_csv(dir, file_dat, vires, timeout=timeout)

    if remove_dat_files
        rm(joinpath(dir, file_dat)) # remove the .dat files
    end

    # idle_and_print_messages(vires.SCP, 3.0)
end
function set_driver_behavior(
    vires::ViresConnection,
    def::DriverBehavior;
    #name::Union(Nothing,String)=nothing,
    name::Union{Void,String}=nothing,
    #id::Union(Nothing,Int)=nothing,
    id::Union{Void,Int}=nothing,
    #visible::Union(Nothing,Bool)=nothing,
    visible::Union{Void,Bool}=nothing,
    timeout::Float64=TIMEOUT_DEFAULT
    )

    write_and_wait_for_mirrored_message(vires.SCP,
        SCPMessage("SCP","TaskControl", get_xml_player_driver_behavior_normalized(def, name=name, id=id, visible=visible)),
        "set driver behavior timed out",
        timeout=timeout)
end
function convert_recording_to_csv(
    dir::String,
    file_dat::String,
    vires::ViresConnection;
    timeout::Float64=TIMEOUT_LONG_DEFAULT
    )

    write(vires.SCP, SCPMessage(get_xml_replay_convert(dir, file_dat)))

    finished = false
    start_time = time()
    while !finished && time() - start_time < timeout
        root = message_payload_to_etree(read(vires.SCP, SCPMessage))
        if root.name == "Info"
            index = findfirst(elem->isa(elem, ETree) && elem.name == "Message", root.elements)
            if index != 0
                element = root.elements[index]
                finished = get(element.attr, "text", "") == "Data export complete."
            end
        end
    end
end

function create_camera_directly_above_ego(
    vires::ViresConnection;
    distance::Float64=50.0
    )

    write_and_wait_for_mirrored_message(vires.SCP,
        SCPMessage("<Player name=\"Ego\" visible=\"true\"/>"),
        "make player visible timed out")
    write_and_wait_for_mirrored_message(vires.SCP,
        SCPMessage("<Camera><PosTether player=\"Ego\" distance=\"" * string(distance) *
                             "\" azimuth=\"3.14159\" elevation=\"1.5708\" slew=\"1\" stayLevel=\"true\"/><ViewPlayer player=\"Ego\"/></Camera>"),
        "set camera timed out")
end
function create_camera_relative_to_ego(
    vires::ViresConnection;
    dx::Float64=0.0,
    dy::Float64=0.0,
    dz::Float64=50.0,
    )

    write_and_wait_for_mirrored_message(vires.SCP,
        SCPMessage("<Player name=\"Ego\" visible=\"true\"/>"),
        "make player visible timed out")
    write_and_wait_for_mirrored_message(vires.SCP,
        SCPMessage("<Camera><PosRelative player=\"Bot03\" dx=\"" * string(dx) * "\" " *
                                                       "dy=\"" * string(dy) * "\" " *
                                                       "dz=\"" * string(dz) * "\"/>" *
                            "<ViewRelative dh=\"0.0\" dp=\"1.5708\" dr=\"0.0\"/></Camera>"),
        "set camera timed out")
end

function register_to_authorize_framesteps!(
    vires::ViresConnection, timeout::Float64=TIMEOUT_DEFAULT;
    source="USER")

    mask = -1

    println(vires.registered_to_authorize_framestep)

    if !vires.registered_to_authorize_framestep

        write(vires.SCP, SCPMessage("SCP","TaskControl", get_xml_query_sync(source=source)))

        finished = false
        start_time = time()
        while !finished && time() - start_time < timeout
            root = message_payload_to_etree(read(vires.SCP, SCPMessage))
            if root.name == "Reply" && get(root.attr, "entity", "") == "RDB"
                index = findfirst(elem->isa(elem, ETree) && elem.name == "Sync", root.elements)
                if index != 0
                    element = root.elements[index]
                    mask = int(get(element.attr, "mask", "-1"))
                    finished = get(element.attr, "source", "") == source

                    if get(element.attr, "mode", "") != "frame"
                        mask = -1
                    end
                end
            end
        end

        if mask ≥ 0
            vires.registered_to_authorize_framestep = true
            vires.sync_mask = convert(Uint32, mask)
        end
    end
end
function unregister_to_authorize_framesteps!(
    vires::ViresConnection,
    timeout::Float64=TIMEOUT_DEFAULT;
    source="USER")

    freed = true

    if vires.registered_to_authorize_framestep

        write(vires.SCP, SCPMessage("SCP","TaskControl", get_xml_query_sync(source=source, delete=true)))

        finished = false
        start_time = time()
        while !finished && time() - start_time < timeout
            root = message_payload_to_etree(read(vires.SCP, SCPMessage))
            if root.name == "Reply" && get(root.attr, "entity", "") == "RDB"
                index = findfirst(elem->isa(elem, ETree) && elem.name == "Sync", root.elements)
                if index != 0
                    element = root.elements[index]
                    finished = get(element.attr, "source", "") == source
                    freed = get(element.attr, "mode", "") == "free"
                end
            end
        end
    end

    freed
end

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using ParallelAccelerator
using DocOpt

include("image-notypes.jl")
using .Image

#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)
ParallelAccelerator.ParallelIR.PIRNumThreadsMode(2)

@acc begin
# Block Jacobi preconditioner
# Used in singleScaleOpticalFlow below
# Needs to be toplevel function for use with ParallelAccelerator
@inline function blockJacobiPreconditioner(Ix, Iy, ru, rv, lam)
    a = Ix.*Ix.+4.0f0*lam
    b = Ix.*Iy
    c = Iy.*Iy.+4.0f0*lam
    invdet = 1.0f0./(a.*c-b.*b)
    rru = (c.*ru-b.*rv).*invdet
    rrv = (a.*rv-b.*ru).*invdet
    rru, rrv
end

# Solve optical flow problem at one scale
# Flow problem is formulated as the solution of a linear system where matrix is sparse
# Use an iterative method to
# Pre: size(i1)==size(i2)
function singleScaleOpticalFlow(i1, i2, lam, ni)
# Returns two Matrix{Float32} size=size(i1)
  (w, h) = size(i1)
  # Compute the gradient of i1
  # This is basically two stencils, one for x and one for y
  # The stencil is [1 8 0 -8 -1]/12 and leave zeros in the boundary
  Ix = Array(Float32, w, h)
  runStencil(Ix, i1, 1, :oob_dst_zero) do b, a
    b[0,0] = (a[-2,0] + 8.0f0 * a[-1,0] - 8.0f0 * a[1,0] - a[2,0])/12.0f0
  end
  Iy = Array(Float32, w, h)
  runStencil(Iy, i1, :oob_dst_zero) do b, a
    b[0,0] = (a[0,-2] + 8.0f0 * a[0,-1] - 8.0f0 * a[0,1] - a[0,2])/12.0f0
  end
  # Compute the time partial derivative
  It = i1-i2

  # Initialise for iterative solve
  xu = zeros(Float32, w, h)
  xv = zeros(Float32, w, h)
  ru = -Ix.*It
  rv = -Iy.*It
  zu, zv = blockJacobiPreconditioner(Ix, Iy, ru, rv, lam)
  pu, pv = zu, zv
  rsold = sum(ru.*zu)+sum(rv.*zv)

  # Iterate to solution
  for i in 1:ni
    # The matrix vector product is equivalent to a stencil
    # Want:
    #   -lam for (-1,0), (+1,0), (0,-1), and (0,+1)
    #   Ix*Iy for (0,0) of the other array
    #   Ix^2+4*lam and Iy^2+4*lam for (0,0) respectively (Ix for u, Iy for v)
    # At boundaries just include terms that exist
    Apu = Array(Float32, w, h)
    Apv = Array(Float32, w, h)
    runStencil(Apu, Apv, pu, pv, Ix, Iy, :oob_src_zero) do Apu, Apv, pu, pv, Ix, Iy
      ix = Ix[0,0]
      iy = Iy[0,0]
      Apu[0,0] = ix * (ix*pu[0,0] + iy*pv[0,0]) + lam*(4.0f0*pu[0,0]-(pu[-1,0]+pu[1,0]+pu[0,-1]+pu[0,1]))
      Apv[0,0] = iy * (ix*pu[0,0] + iy*pv[0,0]) + lam*(4.0f0*pv[0,0]-(pv[-1,0]+pv[1,0]+pv[0,-1]+pv[0,1]))
    end
    pTAp = sum(pu.*Apu)+sum(pv.*Apv)
    alpha = rsold/pTAp
    xu += alpha*pu
    xv += alpha*pv
    ru -= alpha*Apu
    rv -= alpha*Apv
    zu, zv = blockJacobiPreconditioner(Ix, Iy, ru, rv, lam)
    rsnew = sum(ru.*zu)+sum(rv.*zv)
    beta = rsnew/rsold
    pu = zu + beta*pu
    pv = zv + beta*pv
    rsold = rsnew
  end

  return xu, xv
end

# Solve optical flow problem at all scales
# Pre: size(i1)==size(i2)
@fastmath function multiScaleOpticalFlow(i1, i2, lam, ni, ns)
# Returns two Matrix{Float32} size=size(i1)
  (w, h) = size(i1)
  if ns==1
    return singleScaleOpticalFlow(i1, i2, lam, ni)
  end
  scale = Float32((Float64(w)/50.0)^(-1.0/Float32(ns)))
  local u, v # initialised in the first iteration of the loop
  i = ns
  while i >= 0
    cw = round(Int, Float32(w)*scale^Float32(i))
    ch = round(Int, Float32(h)*scale^Float32(i))
    # println("cw,ch=",cw,",",ch)
    # println("Doing scale: ", cw, "x", ch)
    # Downsample the images
    si1 = Image.downSample(i1, cw, ch)
    si2 = Image.downSample(i2, cw, ch)
    # Interpolate the current flow and warp image 2 by current flow
    if i==ns
      wi = si2
      u = zeros(Float32, cw, ch)
      v = zeros(Float32, cw, ch)
    else
      # Coarsest scale, no flow yet, use zeros, warp is image 2
      u, v = Image.interpolateFlow(u, v, cw, ch)
      wi = Image.warpMotion(si2, u, v, si1)
    end
    # Compute the flow from image 1 to warp
    # println("checksum si1 wi: ", sum(si1), " ", sum(wi))
    du, dv = singleScaleOpticalFlow(si1, wi, lam, (i==0 ? ni*5 : ni))
    # Update the flow
    u += du
    v += dv
    i -= 1
  end
  return u, v
end

function multiFrameOpticalFlow(frames, nframes)
  cartesianarray(Tuple{Matrix{Float32}, Matrix{Float32}}, (nframes - 1,)) do i
    multiScaleOpticalFlow(frames[i], frames[i+1], Float32(0.025), 100, 44)
  end
end

end

function main()
    doc = """opt-flow.jl

Horn-Schunck multi-frame optical flow estimator.

Usage:
  opt-flow.jl -h | --help
  opt-flow.jl [--image-name-prefix=<image-name-prefix>] [--image-name-suffix=<image-name-suffix>] [--num-frames=<num-frames>]

Options:
  -h --help                     Show this screen.
  --image-name-prefix=<image-name-prefix>  Specify an image filename prefix [default: small_0].
  --image-name-suffix=<image-name-suffix>  Specify an image filename prefix [default: .dat].
  --num-frames=<num-frames>                Specify a number of frames (at least 2) [default: 2].

Assumes that image files are in the same directory as opt-flow.jl.
"""
    arguments = docopt(doc)

    fname_prefix = arguments["--image-name-prefix"]
    fname_suffix = arguments["--image-name-suffix"]
    nframes = parse(Int, arguments["--num-frames"])

    images_dir = dirname(@__FILE__)

    filenames = [ string(fname_prefix, string(i < 10 ? "0" : "", i), fname_suffix) for i = 1:nframes ]
    filepaths = [ joinpath(images_dir, i) for i in filenames ]

    println("nframes = ", nframes)
    assert(nframes > 1)
    frames = Array(Matrix{Float32}, nframes)
    for i = 1:nframes
        frames[i] = Image.readImage(filepaths[i])
        if i > 1 && size(frames[i])!=size(frames[i-1])
            write(STDERR, "images are different sizes")
            exit(-1)
        end
    end

    println("filenames = ", filenames)
    println("checksums = ", map(sum, frames))
    println("Image size: ", size(frames[1], 1), "x", size(frames[1], 2))

    tic()
    multiFrameOpticalFlow(frames, 1)
    println("SELFPRIMED ", toq())

    tic()
    u,v = multiFrameOpticalFlow(frames, nframes)
    selftimed = toq()
    Image.writeFlo(u[1], v[1], "out.flo")
    println("checksum: ", sum(u[1]), " ", sum(v[1]))
    println("SELFTIMED ", selftimed)
end

main()


text:
using Crux
using Test
using POMDPModels
using Flux
using LinearAlgebra
using CUDA
using Distributions

# Constant Layer
c1 = ConstantLayer(ones(10))
@test Crux.device(c1) == cpu
@test c1(rand(100)) == c1.vec
@test Flux.params(c1)[1] == c1.vec

c2 = c1 |> gpu
@test Crux.device(c2) == gpu
@test c2(rand(100)) == c2.vec

# Distribution stuff
objs = [:up, :down]
o = ObjectCategorical(objs)
@test o.objs == objs
@test o.cat.p == Categorical(2).p

@test rand(o) in objs
@test size(rand(o, 10)) == (10,)

@test logpdf(o, :up) == logpdf(o, :down)
@test size(logpdf(o, rand(o,10))) == (1,10)

## Flux Stuff
W = rand(2, 5)
b = rand(2)

predict(x) = (W * x) .+ b
loss(x, y) = sum((predict(x) .- y).^2)

x, y = rand(5), rand(2) # Dummy data
l = loss(x, y) # ~ 3

θ = Flux.params(W, b)
grads = Flux.gradient(() -> loss(x, y), θ)

@test norm(grads) > 1

##  MultitaskDecay Schedule
m = MultitaskDecaySchedule(10, [1,2,3])
l = Crux.LinearDecaySchedule(1.0, 0.1, 10)

for i=1:10
    @test m(i) == l(i)
end

for i=11:20
    @test m(i) == l(i-10)
end

for i=21:30
    @test m(i) == l(i-20)
end

m = MultitaskDecaySchedule(10, [1,2,1])

for i=1:10
    @test m(i) == l(i)
end

for i=11:20
    @test m(i) == l(i-10)
end

for i=21:30
    @test m(i) == l(i-10)
end

@test m(31) == 0.1
@test m(0) == 1


text:
#!/usr/bin/env julia

using GaussDCA
using FastaIO
using DelimitedFiles
using ArgParse

# Some needed scripts
# PATH_SCRIPT = "/home/jrodrig5/landscapes_coevolution/src/giancarlo"
# include(joinpath(PATH_SCRIPT,"read_script.jl"))
# include(joinpath(PATH_SCRIPT,"compute_fi_fij.jl"))

# PATH_SCRIPT = "/home/fenix/projects/sars_cov2/src/energies"
include(joinpath(@__DIR__,"tools_energy_analysis.jl"))
include(joinpath(@__DIR__,"compute_fi_fij.jl"))
# @__DIR__

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--refseq", "-r"
            help = "Path to the reference sequence in FASTA format"
            required = true
        "--ali", "-i"
            help = "Path to the alignment in FASTA format"
            required = true
        "--root_parameters", "-p"
            help = "Root path to the gzip compressed files with the h and J parameters. They should be at root.h.gz and root.J.gz"
            required = true
        "--outfile", "-o"
            help = "Path for the output file"
            required = true
    end

    return parse_args(s)
end

let alphabet = [ 1,21, 2, 3, 4, 5, 6, 7, 8,21, 9,10,11,12,21,13,14,15,16,17,21,18,19,21,20]
    # A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y
    global letter2num
    function letter2num(c::Union{Char,UInt8})
        i = UInt8(c) - 0x40
        1 <= i <= 25 && return alphabet[i]
        return 21
    end
end

aas = Dict(1 => "A",
    2 => "C",
    3 => "D",
    4 => "E",
    5 => "F",
    6 => "G",
    7 => "H",
    8 => "I",
    9 => "K",
    10 => "L",
    11 => "M",
    12 => "N",
    13 => "P",
    14 => "Q",
    15 => "R",
    16 => "S",
    17 => "T",
    18 => "V",
    19 => "W",
    20 => "Y",
    21 => "-")

function read_par(path_par::AbstractString)

   # read fields
   h_name = path_par*".h"
   h = readdlm(h_name)

   q, N = size(h)

   J_name = path_par*".J"
   J_m = readdlm(J_name)
   J_m += J_m'

   J = zeros(q,q,N,N);

   for i=1:N
      for j = (i+1):N
           J[:,:,i,j] = J_m[(1+q*(i-1)):(q+q*(i-1)),(1+q*(j-1)):(q+q*(j-1))]
      J[:,:,j,i] = (J_m[(1+q*(i-1)):(q+q*(i-1)),(1+q*(j-1)):(q+q*(j-1))])'
       end
   end

   return h,J

end

function fasta2matrix(filename::AbstractString, max_gap_fraction::Real)

    f = FastaReader(filename)

    max_gap_fraction = Float64(max_gap_fraction)

    # pass 1

    seqs = Int[]
    inds = Int[]
    fseqlen = 0

    for (name, seq) in f
        ngaps = 0
        if f.num_parsed == 1
            ls = length(seq)
            resize!(inds, ls)
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    fseqlen += 1
                    inds[fseqlen] = i
                    c == '-' && (ngaps += 1)
                end
            end
        else
            ls = length(seq)
            ls == length(inds) || error("inputs are not aligned")
            tstfseqlen = 0
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    tstfseqlen += 1
                    inds[tstfseqlen] == i || error("inconsistent inputs")
                    c == '-' && (ngaps += 1)
                end
            end
            tstfseqlen == fseqlen || error("inconsistent inputs")
        end
        ngaps / fseqlen <= max_gap_fraction && push!(seqs, f.num_parsed)
    end

    length(seqs) > 0 || error("Out of $(f.num_parsed) sequences, none passed the filter (max_gap_fraction=$max_gap_fraction)")

    # pass 2
    # Z = Array{Int8}(fseqlen, length(seqs))
    Z = zeros(Int8, (fseqlen, length(seqs)))

    seqid = 1
    for (name, seq) in f
        seqs[end] < f.num_parsed && break
        seqs[seqid] == f.num_parsed || continue
        for i = 1:fseqlen
            c = seq[inds[i]]
            Z[i, seqid] = letter2num(c)
        end
        seqid += 1
    end
    @assert seqid == length(seqs) + 1

    close(f)

    return Z'
end

function main()
    parsed_args = parse_commandline()

    reference_sequence_file = parsed_args["refseq"]
    input_ali = parsed_args["ali"]
    rootname_parameters = parsed_args["root_parameters"]
    output_file = parsed_args["outfile"]

    # Load info for IND model
    Z = GaussDCA.read_fasta_alignment(input_ali,0.9)
    theta = 0.2
    N, M = size(Z)
    q = Int(maximum(Z))
    fi_tens, _, _, _, Meff, _ = compute_new_fi_fij(Z,q,theta)
    fi_tens_pseudo = add_pseudocounts(fi_tens, 0.01)

    # Load info for DCA model
    h,J = read_par(rootname_parameters)
    S = fasta2matrix(reference_sequence_file,0.9)

    # To vector
    ref_seq = S[1,:]

    # Compute energies
    energy_ref_dca = compute_energy_single_sequence(h,J,ref_seq)
    energy_ref_ind = compute_energy_independent_model(fi_tens_pseudo, ref_seq)

    file_out = open(output_file,"w")
    tam = length(ref_seq)
    # effects = Array{Float64,2}(tam,19)
    effects = zeros(Float64,tam,19)
    println(file_out, "position\taa_wildtype\taa_mutation\tenergy_IND\tscore_IND\tenergy_DCA\tscore_DCA")
    for pos in 1:tam
        original = ref_seq[pos]
        for mut in 1:21
            mut_seq = copy(ref_seq)
            mut_seq[pos] = mut
            energy_mut_dca = compute_energy_single_sequence(h,J,mut_seq)
            energy_mut_ind = compute_energy_independent_model(fi_tens_pseudo, mut_seq)

            diff_energy_dca = -(energy_mut_dca - energy_ref_dca)
            diff_energy_ind = -(energy_mut_ind - energy_ref_ind)
            println(file_out, pos,"\t", aas[original],"\t",aas[mut],"\t",round(energy_mut_ind, digits=4),"\t",round(diff_energy_ind, digits=4),"\t",round(energy_mut_dca, digits=4),"\t",
            round(diff_energy_dca, digits=4))
        end
#         println(file_out, pos,"\t", aas[original],"\t","*","\t","NA","\t","NA")
    end
end

main()

max_stars_repo_name: juan-rodriguez-rivas/covmut

max_issues_repo_name: juan-rodriguez-rivas/covmut

max_forks_repo_name: juan-rodriguez-rivas/covmut


text:
"""
    UnitCell{T<:AbstractFloat}

A type defining a unit cell.

# Fields
$(TYPEDFIELDS)
"""
struct UnitCell{T<:AbstractFloat}

    "Number of spatial dimensions."
    D::Int

    "Orbitals per unit cell."
    n::Int

    "Matrix where the columns are the lattice vectors."
    lattice_vecs::Matrix{T}

    "Matrix where the columns are the reciprocal lattice vectors."
    reciprocal_vecs::Matrix{T}

    "Matrix where the columns are then basis vectors."
    basis_vecs::Matrix{T}
end

"""
    UnitCell(lattice_vecs::Matrix{T}, basis_vecs::Matrix{T}) where {T<:AbstractFloat}

Constrcuts a [`UnitCell`](@ref).
"""
function UnitCell(lattice_vecs::Matrix{T}, basis_vecs::Matrix{T}) where {T<:AbstractFloat}

    @assert 1 <= size(lattice_vecs,1) <= 3
    @assert size(basis_vecs,2) >= 1
    @assert size(lattice_vecs,1)==size(basis_vecs,1)
    @assert size(lattice_vecs,1)==size(lattice_vecs,2)

    # dimension of unit cell
    D = size(lattice_vecs,1)

    # number of orbitals in unit cell
    n = size(basis_vecs,2)

    # calculating reciprocal lattice vectors corresponding to lattice vectors
    reciprocal_vecs = 2π * inv(lattice_vecs)'

    return UnitCell{T}(D, n, lattice_vecs, reciprocal_vecs, basis_vecs)
end

"""
    UnitCell(lattice_vecs::AbstractVector{Vector{T}}, basis_vecs::AbstractVector{Vector{T}})
        where {T<:AbstractFloat}

Constrcuts a [`UnitCell`](@ref).
"""
function UnitCell(lattice_vecs::AbstractVector{Vector{T}}, basis_vecs::AbstractVector{Vector{T}}) where {T<:AbstractFloat}

    return UnitCell(hcat(lattice_vecs...), hcat(basis_vecs...))
end

UnitCell(; lattice_vecs, basis_vecs) = UnitCell(lattice_vecs, basis_vecs)

"""
    Base.show(io::IO, uc::UnitCell{T}) where {T}
    Base.show(io::IO, ::MIME"text/plain", uc::UnitCell{T}) where {T}

Show unit cell.
"""
Base.show(io::IO, uc::UnitCell{T}) where {T} = print(io,"UnitCell{$T}(D=$(uc.D), n=$(uc.n))")
function Base.show(io::IO, ::MIME"text/plain", uc::UnitCell{T}) where {T}

    (; D, n, lattice_vecs, reciprocal_vecs, basis_vecs) = uc
    print(io, "UnitCell{$T}:\n")
    print(io," • D = $D\n")
    print(io," • n = $n\n")
    print(io," • lattice_vecs =\n")
    show(io,"text/plain",lattice_vecs)
    print(io,"\n")
    print(io," • reciprocal_vecs =\n")
    show(io,"text/plain",reciprocal_vecs)
    print(io,"\n")
    print(io," • basis_vecs =\n")
    show(io,"text/plain",basis_vecs)
    return nothing
end

"""
    loc_to_pos!(r::AbstractVector{T}, l::AbstractVector{Int},
        unit_cell::UnitCell{T}) where {T}

Calculate the position `r` of a unit cell at location `l`.
"""
function loc_to_pos!(r::AbstractVector{T}, l::AbstractVector{Int}, unit_cell::UnitCell{T}) where {T}

    (; D, lattice_vecs) = unit_cell
    @assert length(r) == length(l) == D

    fill!(r,0.0)
    for d in 1:D
        @views @. r += l[d] * lattice_vecs[:,d]
    end

    return nothing
end

loc_to_pos!(; r, l, unit_cell) = loc_to_pos!(r, l, unit_cell)

"""
    loc_to_pos!(r::AbstractVector{T}, l::AbstractVector{Int}, s::Int,
        unit_cell::UnitCell{T}) where {T}

Calculate the position `r` of a orbital `o` at location `l`.
"""
function loc_to_pos!(r::AbstractVector{T}, l::AbstractVector{Int}, o::Int, unit_cell::UnitCell{T}) where {T}

    loc_to_pos!(r,l,unit_cell)
    @views @. r += unit_cell.basis_vecs[:,o]

    return nothing
end

loc_to_pos!(; r, l, o, unit_cell) = loc_to_pos!(r, l, o, unit_cell)

"""
    loc_to_pos(l::AbstractVector{Int}, unit_cell::UnitCell{T})
        where {T}

Return the position `r` of a unit cell at location `l`.
"""
function loc_to_pos(l::AbstractVector{Int}, unit_cell::UnitCell{T}) where {T}

    r = zeros(T,unit_cell.D)
    loc_to_pos!(r,l,unit_cell)

    return r
end

loc_to_pos(; l, unit_cell) = loc_to_pos(l, unit_cell)

"""
    loc_to_pos(l::AbstractVector{Int}, s::Int, unit_cell::UnitCell{T})
        where {T}

Return the position `r` of a orbital `o` at location `l`.
"""
function loc_to_pos(l::AbstractVector{Int}, o::Int, unit_cell::UnitCell{T}) where {T}

    r = zeros(T,unit_cell.D)
    loc_to_pos!(r,l,o,unit_cell)

    return r
end

loc_to_pos(; l, o, unit_cell) = loc_to_pos(l, o, unit_cell)

"""
    displacement_to_vec!(Δr::AbstractVector{T}, Δl::AbstractVector{Int},
        o₁::Int, o₂::Int, unit_cell::UnitCell{T}) where {T}

Computes the position space displacement vector `Δr` corresponding to a displacement definition given by initial and final
orbitals `o₁` and `o₂` in the unit cell respectively, along with a displacement in unit cells `Δl`.

# Arguments
- `Δr::AbstractVector{T}`: displacement vector in position space.
- `Δl::AbstractVector{Int}`: displacement in unit cells.
- `o₁::Int`: initial orbital in unit cell.
- `o₂::Int`: final orbital in unit cell.
- `unit_cell::UnitCell{T}`: unit cell.
"""
function displacement_to_vec!(Δr::AbstractVector{T}, Δl::AbstractVector{Int}, o₁::Int, o₂::Int, unit_cell::UnitCell{T}) where {T}

    (; D, n, lattice_vecs, basis_vecs) = unit_cell
    @assert length(Δr) == length(Δl) == D
    @assert 1 <= o₁ <= n
    @assert 1 <= o₂ <= n

    fill!(Δr,0.0)
    for d in 1:D
        @views @. Δr += Δl[d] * lattice_vecs[:,d]
    end
    @views @. Δr += basis_vecs[:,o₂] - basis_vecs[:,o₁]

    return nothing
end

displacement_to_vec!(; Δr, Δl, o₁, o₂, unit_cell) = displacement_to_vec!(Δr, Δl, o₁, o₂, unit_cell)

"""
    displacement_to_vec(Δl::AbstractVector{Int}, o₁::Int, o₂::Int,
        unit_cell::UnitCell{T}) where {T}

Returns the position space displacement vector `Δr` corresponding to a displacement definition given by initial and final
orbitals `o₁` and `o₂` in the unit cell respectively, along with a displacement in unit cells `Δl`.
"""
function displacement_to_vec(Δl::AbstractVector{Int}, o₁::Int, o₂::Int, unit_cell::UnitCell{T}) where {T}

    Δr = zeros(T,unit_cell.D)
    displacement_to_vec!(Δr, Δl, o₁, o₂, unit_cell)

    return Δr
end

displacement_to_vec(; Δl, o₁, o₂, unit_cell) = displacement_to_vec(Δl, o₁, o₂, unit_cell)


text:
module TestBioArgParse

using Base.Test, BioArgParse, ArgParse
using Bio: Seq

@testset "Parsing sequences" begin

    s = ArgParseSettings()
    s.add_help = true

    @add_arg_table s begin
        "--dna1", "-d"
            help = "A DNA Sequence"
            required = false
            arg_type = BioSequence{DNAAlphabet{4}}
        "--dna2", "-D"
            help = "A DNA Sequence, stored in 2-bit"
            required = false
            arg_type = BioSequence{DNAAlphabet{2}}
        "--rna1", "-r"
            help = "An RNA Sequence"
            required = false
            arg_type = BioSequence{RNAAlphabet{4}}
        "--rna2", "-R"
            help = "An RNA Sequence, stored in 2-bit"
            required = false
            arg_type = BioSequence{RNAAlphabet{2}}
        "--prot", "-P"
            help = "An amino acid sequence"
            required = false
            arg_type = BioSequence{AminoAcidAlphabet}
        "--nuc", "-n"
            help = "A DNA nucleotide"
            required = false
            arg_type = DNANucleotide
        "--rnuc", "-q"
            help = "An RNA nucleotide"
            required = false
            arg_type = RNANucleotide
    end

    runtest(args) = parse_args(args, s)

    @test runtest([]) == Dict{AbstractString, Any}("dna1"=>nothing,
                                                   "dna2"=>nothing,
                                                   "rna1"=>nothing,
                                                   "rna2"=>nothing,
                                                   "prot"=>nothing,
                                                   "nuc"=>nothing,
                                                   "rnuc"=>nothing)

    answer1 = Dict{AbstractString, Any}("dna1"=>BioSequence{DNAAlphabet{4}}("atcgatcg"),
                                        "dna2"=>BioSequence{DNAAlphabet{2}}("aaatttcccggg"),
                                        "rna1"=>BioSequence{RNAAlphabet{4}}("aucgaucg"),
                                        "rna2"=>BioSequence{RNAAlphabet{2}}("aaauuucccggg"),
                                        "prot"=>BioSequence{AminoAcidAlphabet}("andeqhm"),
                                        "nuc"=>DNA_A,
                                        "rnuc"=>RNA_U
                                       )
    @test runtest(["--dna1",
                   "atcgatcg",
                   "--dna2",
                   "aaatttcccggg",
                   "--rna1",
                   "aucgaucg",
                   "--rna2",
                   "aaauuucccggg",
                   "--prot",
                   "andeqhm",
                   "--nuc",
                   "A",
                   "--rnuc",
                   "U"]) == answer1

    @test runtest(["-d",
                   "atcgatcg",
                   "-D",
                   "aaatttcccggg",
                   "-r",
                   "aucgaucg",
                   "-R",
                   "aaauuucccggg",
                   "-P",
                   "andeqhm",
                   "-n",
                   "a",
                   "-q",
                   "u"]) == answer1

end

end


text:
####################

using Soss
using SampleChainsDynamicHMC

p = @model k,x begin
    σ ~ Exponential()
    α ~ Normal()
    β ~ Normal() |> iid(k)
    yhat = α .+ x * β
    y ~ For(eachindex(yhat)) do j
        Normal(μ = yhat[j], σ=σ)
    end
    return y
end;

m = @model k,x begin
    σ ~ Exponential()
    α ~ Normal()
    β ~ Normal() |> iid(k)
    yhat = α .+ x * β
    y ~ For(eachindex(yhat)) do j
        Normal(μ = yhat[j], σ=σ)
    end
    return y
end;

k=1

x = randn(1000,k)

y = rand(p(k=k,x=x))

post = m(k=k, x=x) | (;y)

ℓ(nt) = logdensity(post, nt) ;

t = xform(post);

chains = initialize!(DynamicHMCChain, ℓ, t);
drawsamples!(chains, 199)

drawsamples!(chains, 800)

pred = [predict(Soss.withmeasures(m)(k=k, x=x), c) for c in chains];
yhat = hcat(getproperty.(pred, :yhat)...)
d = getproperty.(pred, :_y_dist)
ppc = hcat((cdf.(dj.data, y) for dj in d)...)

using StatsBase
using AverageShiftedHistograms

function grassfire(v,m)
    ordv = ordinalrank(v)
    n = length(v)
    (lo,hi) = extrema(m)
    Δ = (hi-lo)/20
    lo = lo - Δ
    hi = hi + Δ
    a = ash(ordv, m[:,1]; rngx=range(1-n/20,n+n/20,length=800), rngy=range(lo,hi,length=600), kernelx=AverageShiftedHistograms.Kernels.uniform, kernely=AverageShiftedHistograms.Kernels.uniform)
    for j in 2:size(m,2)
        ash!(a,ordv, m[:,j])
    end
    ash!(a; mx=3,my=10)

    # [ash(ordv, m[:,j]; rngx=range(1-n/20,n+n/20,length=800), rngy=range(lo,hi,length=600)) for j in 1:size(m,2)]
end

# a = grassfire(y,ppc)

a = grassfire(vec(mean(yhat, dims=2)),ppc)
using Plots
plot(a)

b = grassfire(vec(x),y .-yhat)
plot(b)

using Plots

scatter(ordinalrank(y), ppc)

using NestedTuples: rmap
rmap(length, chains)

# chain1 = initialize!(DynamicHMCChain, ℓ, t);
# drawsamples!(chain1,1000)

# chain2 = initialize!(DynamicHMCChain, ℓ, t);
# drawsamples!(chain2,1000)

# chain3 = initialize!(DynamicHMCChain, ℓ, t);
# drawsamples!(chain3,1000)

# chains = MultiChain(chain1, chain2, chain3)


text:
# Belief Propagation                #
# --------------------------------- #
@rule NormalMeanPrecision(:out, Marginalisation) (m_μ::PointMass, m_τ::PointMass) = NormalMeanPrecision(mean(m_μ), mean(m_τ))

@rule NormalMeanPrecision(:out, Marginalisation) (m_μ::UnivariateNormalDistributionsFamily, m_τ::PointMass) = begin
    m_μ_mean, m_μ_cov = mean_cov(m_μ)
    return NormalMeanPrecision(m_μ_mean, inv(m_μ_cov + inv(mean(m_τ))))
end

# Variational                       #
# --------------------------------- #
@rule NormalMeanPrecision(:out, Marginalisation) (q_μ::Any, q_τ::Any) = NormalMeanPrecision(mean(q_μ), mean(q_τ))

@rule NormalMeanPrecision(:out, Marginalisation) (m_μ::NormalMeanPrecision, q_τ::Any) = NormalMeanPrecision(mean(m_μ), cholinv( cov(m_μ) + cholinv(mean(q_τ)) ))


text:
###############################################################################
#
# Capacitance buffer flow
#
###############################################################################
"""
    buffer_rate(pv::AbstractCapacity{FT}) where {FT<:AbstractFloat}

Return the buffer rate, given
- `pv` [`AbstractCapacity`](@ref) type struct

Note that only symplastic water can be used as capacitance
"""
function buffer_rate(
            pv::PVCurveLinear{FT}
) where {FT<:AbstractFloat}
    return pv.k_refill
end

function buffer_rate(
            pv::PVCurveSegmented{FT}
) where {FT<:AbstractFloat}
    return pv.k_refill * (1 - pv.RWC_apo)
end

###############################################################################
#
# Calculate leaf-level critical_flow
#
###############################################################################
"""
    critical_flow(
                hs::LeafHydraulics{FT},
                ini::FT = FT(0.5)
    ) where {FT<:AbstractFloat}
    critical_flow(
                tree::TreeSimple{FT},
                ini::FT = FT(0.5)
    ) where {FT<:AbstractFloat}

Calculate the critical flow rate (K ≈ 0), given
- `hs` [`LeafHydraulics`](@ref) type struct
- `ini` Initial guess
- `tree` [`TreeSimple`](@ref) type struct

Note, for the safety of no NaN, update critical_flow when ΔP >= -0.01
"""
function critical_flow(
            hs::LeafHydraulics{FT},
            ini::FT = FT(0.5)
) where {FT<:AbstractFloat}
    # calculate maximal flow
    _fh     = (hs.p_ups - hs.p_crt) * hs.k_sla / hs.f_vis;
    _fl     = FT(0);
    _fx     = min((_fh+_fl)/2, ini);
    _ms     = NewtonBisectionMethod{FT}(_fl, _fh, _fx);
    _rt     = SolutionTolerance{FT}(1e-5, 50);
    @inline f(x) = end_pressure(hs, x) - hs.p_crt;
    _solut  = find_zero(f, _ms, _rt);

    if isnan(_solut)
        @warn "E_crit is NaN, please check the settings...";
        @show hs.p_ups;
    end

    return _solut
end

function critical_flow(
            tree::TreeSimple{FT},
            ini::FT = FT(0.5)
) where {FT<:AbstractFloat}
    # calculate maximal flow for whole tree, remember to times leaf area
    _kr = (tree.root).k_max;
    _ks = (tree.stem).k_max;
    _kl = (tree.leaf).k_sla / (tree.leaf).f_vis * (tree.leaf).area;
    _kt = 1 / (1 / _kr + 1 / _ks + 1/ _kl);
    _fh = -1 * (tree.leaf).p_crt * _kt;
    _fl = FT(0);
    _fx = min((_fh+_fl)/2, ini);
    _ms = NewtonBisectionMethod{FT}(_fl, _fh, _fx);
    _rt = SolutionTolerance{FT}(1e-3, 50);
    @inline f(x) = end_pressure(tree, x) - (tree.leaf).p_crt;
    _solut  = find_zero(f, _ms, _rt);

    return _solut
end

###############################################################################
#
# Calculate root flow rates and pressure from total flow rate
#
###############################################################################
"""
    roots_flow!(roots::Array{RootHydraulics{FT},1},
                ks::Array{FT,1},
                ps::Array{FT,1},
                qs::Array{FT,1},
                flow::FT,
                recalculate::Bool
    ) where {FT<:AbstractFloat}
    roots_flow!(roots::Array{RootHydraulics{FT},1},
                ks::Array{FT,1},
                ps::Array{FT,1},
                qs::Array{FT,1},
                flow::FT
    ) where {FT<:AbstractFloat}
    roots_flow!(plant::Union{GrassLikeOrganism{FT},
                             PalmLikeOrganism{FT},
                             TreeLikeOrganism{FT}},
                flow::FT
    ) where {FT<:AbstractFloat}

Recalculate the flow rates in the root from the pressure and conductance
    profiles in each root at non-steady state, given
- `roots` Array of [`RootHydraulics`](@ref) structs
- `ks` Container for conductance in each root layer
- `ps` Container for end xylem pressure in each layer
- `qs` Container for flow rate out of each layer
- `flow` Total flow rate out of the roots
- `recalculate` A paceholder indicator of recalculating root flow (useless)
- `plant` [`AbstractPlantOrganism`](@ref) type struct
"""
function roots_flow!(
            roots::Array{RootHydraulics{FT},1},
            ks::Array{FT,1},
            ps::Array{FT,1},
            qs::Array{FT,1},
            flow::FT,
            recalculate::Bool
) where {FT<:AbstractFloat}
    N = length(roots);

    # flush the values to ks, ps, and qs
    for i in 1:N
        _p,_k = root_pk(roots[i], qs[i]);
        ps[i] = _p;
        ks[i] = _k;
    end

    # use ps and ks to compute the Δq to adjust
    pm = mean(ps);
    for i in 1:N
        qs[i] -= (pm - ps[i]) * ks[i];
    end

    # adjust the qs so that sum(qs) = flow
    q_diff = sum(qs) - flow;
    k_sum  = sum(ks);
    for i in 1:N
        qs[i] -= q_diff * ks[1] / k_sum;
    end

    return nothing
end

function roots_flow!(
            roots::Array{RootHydraulics{FT},1},
            ks::Array{FT,1},
            ps::Array{FT,1},
            qs::Array{FT,1},
            flow::FT
) where {FT<:AbstractFloat}
    count = 0
    while true
        roots_flow!(roots, ks, ps, qs, flow, true);
        count += 1;

        if maximum(ps) - minimum(ps) < 1e-4
            break
        end

        if count>20
            break
        end
    end

    return nothing
end

function roots_flow!(
            plant::Union{GrassLikeOrganism{FT},
                         PalmLikeOrganism{FT},
                         TreeLikeOrganism{FT}},
            flow::FT
) where {FT<:AbstractFloat}
    roots_flow!(plant.roots,
                plant.cache_k,
                plant.cache_p,
                plant.cache_q,
                flow);

    return nothing
end

###############################################################################
#
# Calculate root flow rates and pressure from total flow rate
#
###############################################################################
"""
    xylem_flow(root::RootHydraulics{FT},
               pressure::FT,
               ini::FT = FT(1)
    ) where {FT<:AbstractFloat}

Calculate the flow rate from a given tree base pressure, given
- `root` [`RootHydraulics`](@ref) type struct
- `pressure` Given tree base pressure in `[MPa]`
- `ini` Initial guess
"""
function xylem_flow(
            root::RootHydraulics{FT},
            pressure::FT,
            ini::FT = FT(1)
) where {FT<:AbstractFloat}
    _fh    = (root.p_ups + root.p_osm * root.T_sap / K_25(FT) - pressure) *
             root.k_max / root.f_vis;
    _fl    = -_fh;
    _fx    = min(_fh, ini);
    _ms    = NewtonBisectionMethod{FT}(_fl, _fh, _fx);
    _st    = SolutionTolerance{FT}(1e-4, 50);
    @inline f(x) = end_pressure(root, x) - pressure;
    _solut = find_zero(f, _ms, _st);

    return _solut
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using FinEtools
using FinEtools.MeshQuadrilateralModule
using FinEtoolsFlexStructures.FESetShellQ4SRIModule: FESetShellQ4SRI
using FinEtoolsFlexStructures.FESetShellDSG3Module: FESetShellDSG3
using VisualStructures: plot_nodes, plot_midsurface, render, plot_space_box, plot_solid, space_aspectratio, default_layout_3d, plot_from_json
using PlotlyJS
using JSON

L = 42.0
H = 21.0

fens, fes = Q4block(L, H, 3, 5)
fens.xyz = xyz3(fens)
fes = FESetShellQ4SRI(connasarray(fes))
plots = cat(plot_nodes(fens),
    plot_midsurface(fens, fes; facecolor = "rgb(15, 15, 225)", lwidth = 4),
    dims = 1)
layout = default_layout_3d(autosize=true)
# layout = default_layout_3d(width=400, height=400)
layout[:showLegend] = true

config  = PlotConfig(plotlyServerURL="https://chart-studio.plotly.com", showLink=true)
pl = plot(plots, layout; config = config)

display(pl)

fens, fes = T3block(L, H, 7, 5)
fens.xyz = xyz3(fens)
fes = FESetShellDSG3(connasarray(fes))
plots = cat(plot_nodes(fens),
    plot_midsurface(fens, fes; facecolor = "rgb(15, 155, 225)", lwidth = 4),
    dims = 1)
layout = default_layout_3d(autosize=true)
# layout = default_layout_3d(width=400, height=400)
layout[:showLegend] = true

config  = PlotConfig(plotlyServerURL="https://chart-studio.plotly.com", showLink=true)
pl = plot(plots, layout; config = config)

display(pl)

true


text:
"""
    one-moment bulk Microphysics scheme

Microphysics parameterization based on the ideas of Kessler_1995:
  - condensation/evaporation as relaxation to equilibrium
  - autoconversion
  - accretion
  - rain evaporation
  - rain terminal velocity
"""
module Microphysics

using ClimateMachine.MoistThermodynamics
using CLIMAParameters
using CLIMAParameters.Planet: ρ_cloud_liq, R_v, grav
using CLIMAParameters.Atmos.Microphysics:
    C_drag,
    MP_n_0,
    K_therm,
    τ_cond_evap,
    τ_sub_resub,
    τ_acnv,
    q_liq_threshold,
    E_col,
    D_vapor,
    ν_air,
    N_Sc,
    a_vent,
    b_vent

const APS = AbstractParameterSet

# rain fall speed
export terminal_velocity_single_drop_coeff
export terminal_velocity

# rates of conversion between microphysics categories
export conv_q_vap_to_q_liq
export conv_q_vap_to_q_ice
export conv_q_liq_to_q_rai_acnv
export conv_q_liq_to_q_rai_accr
export conv_q_rai_to_q_vap

"""
    terminal_velocity_single_drop_coeff(ρ)

where:
  - `ρ` is the density of air

Returns the proportionality coefficient between terminal velocity of an
individual water drop and the square root of its radius * g.
"""
function terminal_velocity_single_drop_coeff(
    param_set::APS,
    ρ::FT,
) where {FT <: Real}
    _ρ_cloud_liq::FT = ρ_cloud_liq(param_set)
    _C_drag::FT = C_drag(param_set)

    # terminal_vel_of_individual_drop = v_drop_coeff * (g * drop_radius)^(1/2)
    return sqrt(FT(8 / 3) / _C_drag * (_ρ_cloud_liq / ρ - FT(1)))
end

"""
    terminal_velocity(q_rai, ρ)

where:
  - `q_rai` - rain water specific humidity
  - `ρ`     - density of air

Returns the mass weighted average rain terminal velocity assuming
Marshall Palmer 1948 distribution of rain drops.
"""
function terminal_velocity(param_set::APS, q_rai::FT, ρ::FT) where {FT <: Real}

    rain_w = FT(0)

    if q_rai > FT(0)
        v_c = terminal_velocity_single_drop_coeff(param_set, ρ)
        # gamma(9/2)
        gamma_9_2 = FT(11.631728396567448)
        _grav::FT = grav(param_set)
        _ρ_cloud_liq::FT = ρ_cloud_liq(param_set)
        _MP_n_0::FT = MP_n_0(param_set)

        lambda::FT = (FT(8) * π * _ρ_cloud_liq * _MP_n_0 / ρ / q_rai)^FT(1 / 4)
        rain_w = gamma_9_2 * v_c / FT(6) * sqrt(_grav / lambda)
    end

    return rain_w
end

"""
    conv_q_vap_to_q_liq(q_sat, q)

where:
- `q_sat` - PhasePartition at equilibrium
- `q`     - current PhasePartition

Returns the q_liq tendency due to condensation/evaporation.
The tendency is obtained assuming a relaxation to equilibrium with
constant timescale.
"""
function conv_q_vap_to_q_liq(
    param_set::APS,
    q_sat::PhasePartition{FT},
    q::PhasePartition{FT},
) where {FT <: Real}
    _τ_cond_evap::FT = τ_cond_evap(param_set)
    return (q_sat.liq - q.liq) / _τ_cond_evap
end

"""
    conv_q_vap_to_q_ice(q_sat, q)

where:
- `q_sat` - PhasePartition at equilibrium
- `q`     - current PhasePartition

Returns the q_ice tendency due to condensation/evaporation.
The tendency is obtained assuming a relaxation to equilibrium with
constant timescale.
"""
function conv_q_vap_to_q_ice(
    param_set::APS,
    q_sat::PhasePartition{FT},
    q::PhasePartition{FT},
) where {FT <: Real}

    _τ_sub_resub::FT = τ_sub_resub(param_set)
    return (q_sat.ice - q.ice) / _τ_sub_resub
end

"""
    conv_q_liq_to_q_rai_acnv(q_liq)

where:
- `q_liq` - is the liquid water specific humidity

Returns the q_rai tendency due to collisions between cloud droplets
(autoconversion) parametrized following Kessler 1995.
"""
function conv_q_liq_to_q_rai_acnv(param_set::APS, q_liq::FT) where {FT <: Real}

    _τ_acnv::FT = τ_acnv(param_set)
    _q_liq_threshold::FT = q_liq_threshold(param_set)
    return max(FT(0), q_liq - _q_liq_threshold) / _τ_acnv
end

"""
    conv_q_liq_to_q_rai_accr(param_set, q_liq, q_rai, ρ)

where:
- `param_set` - is an `AbstractParameterSet`
- `q_liq` - is the liquid water specific humidity
- `q_rai` - is the rain water specific humidity
- `ρ` - is the density of air

Returns the q_rai tendency due to collisions between cloud droplets
and rain drops (accretion) parametrized following Kessler 1995.
"""
function conv_q_liq_to_q_rai_accr(
    param_set::APS,
    q_liq::FT,
    q_rai::FT,
    ρ::FT,
) where {FT <: Real}

    accr_rate = FT(0)
    _MP_n_0::FT = MP_n_0(param_set)
    _E_col::FT = E_col(param_set)
    _ρ_cloud_liq::FT = ρ_cloud_liq(param_set)
    _grav::FT = grav(param_set)

    if (q_rai > FT(0) && q_liq > FT(0))
        # terminal_vel_of_individual_drop = v_drop_coeff * drop_radius^(1/2)
        v_c = terminal_velocity_single_drop_coeff(param_set, ρ)

        #gamma(7/2)
        gamma_7_2 = FT(3.3233509704478426)

        accr_coeff::FT =
            gamma_7_2 *
            FT(8)^FT(-7 / 8) *
            π^FT(1 / 8) *
            v_c *
            _E_col *
            (ρ / _ρ_cloud_liq)^FT(7 / 8)

        accr_rate =
            accr_coeff *
            _MP_n_0^FT(1 / 8) *
            sqrt(_grav) *
            q_liq *
            q_rai^FT(7 / 8)
    end
    return accr_rate
end

"""
    conv_q_rai_to_q_vap(q_rai, q, T, p, ρ)

where:
 - q_rai - rain water specific humidity
 - q - current PhasePartition
 - T - temperature
 - p - pressure
 - ρ - air density

Returns the q_rai tendency due to rain evaporation. Parameterized following
Smolarkiewicz and Grabowski 1996.
"""
function conv_q_rai_to_q_vap(
    param_set::APS,
    qr::FT,
    q::PhasePartition{FT},
    T::FT,
    p::FT,
    ρ::FT,
) where {FT <: Real}

    evap_rate = FT(0)
    _R_v = R_v(param_set)
    _D_vapor = D_vapor(param_set)
    _ν_air = ν_air(param_set)
    _grav = grav(param_set)
    _ρ_cloud_liq = ρ_cloud_liq(param_set)
    _N_Sc::FT = N_Sc(param_set)
    _b_vent::FT = b_vent(param_set)
    _a_vent::FT = a_vent(param_set)
    _K_therm::FT = K_therm(param_set)
    _MP_n_0::FT = MP_n_0(param_set)

    if qr > FT(0)
        qv_sat = q_vap_saturation(param_set, T, ρ, q)
        q_v = q.tot - q.liq - q.ice
        S = q_v / qv_sat - FT(1)

        L = latent_heat_vapor(param_set, T)
        p_vs = saturation_vapor_pressure(param_set, T, Liquid())
        G::FT =
            FT(1) / (
                L / _K_therm / T * (L / _R_v / T - FT(1)) +
                _R_v * T / _D_vapor / p_vs
            )

        # gamma(11/4)
        gamma_11_4 = FT(1.6083594219855457)
        v_c = terminal_velocity_single_drop_coeff(param_set, ρ)

        av::FT = sqrt(2 * π) * _a_vent * sqrt(ρ / _ρ_cloud_liq)
        bv::FT =
            FT(2)^FT(7 / 16) *
            gamma_11_4 *
            π^FT(5 / 16) *
            _b_vent *
            (_N_Sc)^FT(1 / 3) *
            sqrt(v_c) *
            (ρ / _ρ_cloud_liq)^FT(11 / 16)

        F::FT =
            av * sqrt(qr) +
            bv * _grav^FT(1 / 4) / (_MP_n_0)^FT(3 / 16) / sqrt(_ν_air) *
            qr^FT(11 / 16)

        evap_rate = min(FT(0), S * F * G * sqrt(_MP_n_0) / ρ)
    end
    return evap_rate
end
end #module Microphysics.jl


text:
using OrdinaryDiffEq, RecursiveArrayTools, Unitful
using LinearAlgebra, Test

@testset "Algorithms" begin
algs = [Euler(),Midpoint(),Heun(),Ralston(),RK4(),SSPRK104(),SSPRK22(),SSPRK33(),
        SSPRK43(),SSPRK432(),BS3(),BS5(),DP5(),DP8(),Feagin10(),Feagin12(),
        Feagin14(),TanYam7(),Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9()]

@testset "Scalar units" begin
    f(y,p,t) = 0.5*y / 3.0u"s"
    u0 = 1.0u"N"
    prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))

    for alg in algs
        @show alg
        sol = solve(prob,alg,dt=1u"s"/10)
    end

    sol = solve(prob,ExplicitRK())
end

@testset "2D units" begin
    f(dy,y,p,t) = (dy .= 0.5.*y ./ 3.0u"s")
    u0 = [1.0u"N" 2.0u"N"
          3.0u"N" 1.0u"N"]
    prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))

    for alg in algs
        @show alg
        sol = solve(prob,alg,dt=1u"s"/10)
    end

    sol = solve(prob,ExplicitRK())
end
end

@testset "Mixed units" begin
@testset "With ArrayPartition" begin
    r0 = [1131.340, -2282.343, 6672.423]u"km"
    v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
    Δt = 86400.0*365u"s"
    μ = 398600.4418u"km^3/s^2"
    rv0 = ArrayPartition(r0,v0)

    function f(dy, y, μ, t)
        r = norm(y.x[1])
        dy.x[1] .= y.x[2]
        dy.x[2] .= -μ .* y.x[1] / r^3
    end

    prob = ODEProblem(f,rv0,(0.0u"s",1.0u"s"),μ)
    for alg in [Tsit5()]
        @show alg
        sol = solve(prob,alg)
    end

    for alg in [    AutoVern6(Rodas5(autodiff=false)),
                    AutoVern7(Rodas5(autodiff=false)),
                    AutoVern8(Rodas5(autodiff=false)),
                    AutoVern9(Rodas5(autodiff=false))]
        @show alg
        @test_broken sol = solve(prob,alg)
    end
end

@testset "Without ArrayPartition" begin
    # coordinate: u = [position, momentum]
    # parameters: p = [mass, force constanst]
    function f_harmonic!(du,u,p,t)
        du[1] = u[2]/p[1]
        du[2] = -p[2]*u[1]
    end

    mass = 1.0u"kg"
    k = 1.0u"N/m"
    p = [mass, k]

    u0 = [1.0u"m", 0.0u"kg*m/s"] # initial values (position, momentum)
    tspan = (0.0u"s", 10.0u"s")
    prob = ODEProblem(f_harmonic!, u0, tspan, p)
    @test_broken solve(prob, Tsit5()).retcode == :Success
end

@testset "Unitful time with unitless state" begin
  using OrdinaryDiffEq, Unitful
  u0 = 30.0
  tspan = (0.0u"s",10.0u"s")
  prob1 = ODEProblem((du,u,t,p)->(du[1] = -0.2u"1/s"*u[1]),[u0],tspan)
  prob2 = ODEProblem((u,t,p)->(-0.2u"1/s"*u[1]),u0,tspan)
  prob3 = ODEProblem((u,t,p)->[-0.2u"1/s"*u[1]],[u0],tspan)
  for prob in [prob1, prob2, prob3]
    solve(prob, Tsit5()).retcode === :Success
  end
end

end


text:
function try_candidate(candidate::Candidate, truth::Candidate)
    codes = fill(Wordle.PRESENT_NOWHERE, 5)

    truth_letter_counts = LetterCounts(truth)

    for (i, char) in enumerate(candidate.chars)
        if char == truth.chars[i]
            @assert truth_letter_counts[char] > 0
            codes[i] = PRESENT_HERE
            truth_letter_counts[char] -= 1
        end
    end

    for (i, char) in enumerate(candidate.chars)
        if truth_letter_counts[char] > 0 && char != truth.chars[i]
            codes[i] = PRESENT_ELSEWHERE
            truth_letter_counts[char] -= 1
        end
    end

    Trial(candidate, Result(tuple(codes...)))
end


text:
using COSMO, JuMP, LinearAlgebra
using Statistics
using JuMP
import Plots as Plots
import ProgressMeter as Pm
MOI = JuMP.MathOptInterface

include("utils.jl")

function MakeLassoOptimizationProblem(A::Matrix, y::Matrix, λ::Float64)
    """
        Phrase the quadratic programming problem for Lasso regularization
        problem.

    """
    m, n = size(A)
    @assert size(y, 1) == m && size(y, 2) == 1 "the label and"*
    " the data matrix doesn't have the maching dimension X:" *
    string(m, "×",  n) * string(";y: ", size(y))
    @assert λ ≤ 1 && λ ≥ 0 "Regularization λ should be between (0, 1)"

    Pb = Pm.ProgressUnknown("Building Optimization Model: ", spinner=true)

    model = Model(with_optimizer(COSMO.Optimizer); bridge_constraints = false)

    set_optimizer_attribute(model, MOI.Silent(), true)
    set_optimizer_attribute(model, "max_iter", 5000*10)

    @variable(model, x[1:n]); Pm.next!(Pb)
    setvalue.(x, A\y); Pm.next!(Pb)
    @variable(model, η[1:n]); Pm.next!(Pb)
    @constraint(model, -η .<= x); Pm.next!(Pb)
    @constraint(model, x .<= η); Pm.next!(Pb)
    @objective(model, Min, λ*sum(η) + sum((A*x - y).^2)); Pm.finish!(Pb)

    return model

end

function MakeQuantileOptimizationProblem(A::Matrix, y::Matrix, λ::Float64)
    """
        Phrase the Quantile regression problem into a Quadratic Programming
        problem

    """
end

mutable struct LassoSCOP
    ## It's just a collection of data.

    A::Matrix # Feature matrix
    y::Matrix # label vector
    μ::Matrix # feature mean
    u::Number # label mean
    Z::Matrix # Starndardized Matrix
    l::Matrix # Zero mean label

    OptModel::Model  # The JuMP model for getting it right.
    λ::Float64       # the regularization parameter.

    LassoPath::Union{Matrix, Nothing} # going along a fixed row is the fixing
    # the feature while varying the lambda quantity.
    λs::Union{Vector, Nothing} # the lambda values.

    function LassoSCOP(A::Matrix, y::Union{Matrix, Vector}, λ::Float64=0.0)
        A = copy(A)
        m, _ = size(A)
        @assert size(y, 1) == m ""*
        "The rows of X should match of the columns of y, but the size of"*
        string("X, Y is: ", size(A), " ", size(y))

        y = copy(y)
        if ndims(y) == 1
            y = reshape(y, (length(y), 1))
        end
        μ = mean(A, dims=1)::Matrix
        u = mean(y)
        Z = A .- μ
        l = y .- u
        OptModel = MakeLassoOptimizationProblem(Z, l, λ)

        return new(A, y, μ, u, Z, l, OptModel, λ)
    end
end

function LassoPath(this::LassoSCOP, tol::Float64=1e-8)
    """
        Analyze the Lasso Problem by drawing a lasso path. It will start with
        a parameter that will make all predictors zero and then solve it
        iterative by chopping the regularization λ by half each iteration.

        **this**:
            An instance of the LassoSCOP
        **tol**:
            If the infinity norm of vector of the change in the weights is
            less than this quantity, then it stops and return all the results.

    """

    u = this.u
    A = this.Z
    y = this.l
    Results = Vector{Vector}()

    function λMax(A, y)
        ToMax = A'*(y .- u)
        ToMax *= 2
        return maximum(abs.(ToMax))
    end

    λ = λMax(A, y)
    λs = Vector{Float64}()
    Changeλ(this, λ)
    x = SolveForx(this)
    dx = Inf
    push!(Results, x)
    push!(λs, λ)
    MaxItr = 100
    pb = Pm.ProgressThresh(tol, "inf norm of δx: ")
    while dx >= tol && MaxItr >= 0
        push!(λs, λ)
        λ /= 2
        MaxItr -= 1
        setvalue.(this.OptModel[:x], x) # warm start!
        Changeλ(this, λ)
        x = SolveForx(this)
        push!(Results, x)
        dx = norm(Results[end - 1] - Results[end], Inf)
        Pm.update!(pb, dx)
    end

    ResultsMatrix = zeros(length(x), length(Results))
    for II ∈ 1:length(Results)
        ResultsMatrix[:, II] = Results[II]
    end

    # Store it.
    this.LassoPath = ResultsMatrix
    this.λs = λs
    return ResultsMatrix, λs
end

function VisualizeLassoPath(
                            this::LassoSCOP,
                            fname::Union{String, Nothing}=nothing,
                            title::Union{String, Nothing}=nothing
                            )
    """
        Make a plots for the lasso path and save it in the pwd.

    """
    @assert isdefined(this, :LassoPath) "Lasso Path not defined for the object"*
    "yet".
    error("Haven't implemented it yet.")
    λs = this.λs
    Paths = this.LassoPath
    Loggedλ = log2.(λs)
    p1 = Plots.heatmap(
        Paths[end:-1:begin, :],
        title=title===nothing ? "Lasso Path" : title
        )
    p2 = plot(Loggedλ, Paths', label=nothing)
    Plots.xlabel!(p2, "log_2(λ)")
    p = Plots.plot(p1, p2, layout=(2, 1))
    Plots.plot!(size=(600, 800))
    Plots.plot!(dpi=400)
    Plots.savefig(p, fname===nothing ? "plot.png" : fname)
    return
end

function CaptureImportantWeights(
                                this::LassoSCOP,
                                top_k::Union{Float64, Int64} = 0.5,
                                threshold::Float64=1e-10
                                )
    """
        Caputre the indices for the most important predictors from the
        regression. Returns the indices of the important weights.

        ---
        **this::LassoSCOP**:
            An instance of the type LassoSCOP

        ---
        **top_k::Union{Float64, Int64} = 0.5**
            When it's in between 0, 1, it's specifying the portion of weights
            that are non-zero, when it's not in in this range and it's positive,
            then it's specifying the number of top weights, measured by the
            absolute value of the weights.

        ---
        **threshold::Float64=1e-10**
            Weigthts less than threshold is viewed as unimportant and it's
            discarded as zero.

    """
    # TODO: test this
    @assert isdefined(this, :LassoPath) "Lasso Path not defined for this"*
    "Object yet. "
    @assert top_k >= 0 "this parameters, should be a positive number"
    top_k = min(top_k, size(this.LassoPath, 1))
    @assert threshold >= 0 "This parameters shouold be a positive number"

    Paths = this.LassoPath
    top_k::Int64 = top_k <= 1 ? ceil(size(Paths, 1)*0.5) : round(top_k)
    for JJ in size(Paths, 2)
        Col = view(Paths, :, JJ)
        NonNegative = sum(abs.(Col) .>= threshold)
        # rank them by abs and returns the indices for top k weights
        if NonNegative >= top_k
            Indices = sortperm(abs.(Col), rev=true)
            return Indices[begin:top_k]
        end
    end
    return Vector{Int64}()
end

function Changeλ(this::LassoSCOP, λ)
    """
        Change the Lasso regularizer of the current model
    """
    model = this.OptModel
    x = model[:x]  # objects can be indexed with symbols!
    η = model[:η]
    y = this.l
    A = this.Z
    @objective(model, Min, λ*sum(η) + sum((A*x - y).^2))

end

function SolveForx(this::LassoSCOP)::Vector{Float64}
    """
        Solve for the weights of the current model,
        given the current configuration of the model.
    """
    optimize!(this.OptModel)
    TernimationStatus = termination_status(this.OptModel)
    # @assert TernimationStatus == MOI.OPTIMAL "Terminated with non-optimal value when solving for x. "*
    # string("The status is: ", TernimationStatus)*"\n this is the results \n $(OptResults)"

    if !(TernimationStatus == MOI.OPTIMAL)
        Warn("\nWarning: convergence status for solver: $(TernimationStatus)")
        Warn("Current Value λ: $(this.λ)")
        PrintTitle("Here is the summary for the solution: ")
        display(solution_summary(this.OptModel))
    end
    return value.(this.OptModel[:x])
end

function Getαβ(this::LassoSCOP, lambda::Float64)
    """
        Get the weights and biases, for the original model
        (The model trained is normalized),
        for a particular regularization value.
    """
    error("Not yet implemented. ")

end

# TODO: Override Base.show for this LASSOPath TYPE.

# TODO: Iterate, Use Base.Iterate, override it.
function Base.iterate(this::LassoSCOP)
    """
        Iterate through the parameters for the Lasso program:
            * λ
            * solutions
    """
    u = this.u
    A = this.Z
    y = this.l

    function λMax(A, y)
        ToMax = A'*(y .- u)
        ToMax *= 2
        return maximum(abs.(ToMax))
    end
    λ = λMax(A, y)
    x = SolveForx(this)
    return λ, x

end

function Base.Iterate(this::LassoSCOP, state::Int64)

end


text:
using Setfield

export MarkovChain

"""
    MarkovChain

`MarkovChain(pars, step)` defines a Markov Chain with global parameters `pars` and transition kernel `step`. Here, `pars` is a named tuple, and `step` is a Soss model that takes arguments `(pars, state)` and returns a `next` value containing the new `pars` and `state`.

NOTE: This is experimental, and may change in the near future.

```jldoctest
mstep = @model pars,state begin
    σ = pars.σ
    x0 = state.x
    x ~ Normal(x0, σ)
    next = (pars=pars, state=(x=x,))
end;

m = @model s0 begin
    σ ~ Exponential()
    pars = (σ=σ,)
    chain ~ MarkovChain(pars, mstep(pars=pars, state=s0))
end;

r = rand(m(s0=(x=2,),));

for s in Iterators.take(r.chain,3)
    println(s)
end

# output

(x = -6.596883394256064,)
(x = 0.48200039561318864,)
(x = -2.838556784903994,)
```
"""
struct MarkovChain{P,D}
    pars :: P
    step :: D
end

struct MarkovChainRand{P,D}
    mc::MarkovChain{P,D}
end

Base.IteratorSize(::MarkovChainRand) = Base.IsInfinite()

Base.eltype(mc::MarkovChainRand{P,D}) where {P,D} = typeof(mc.dist.args.state)

export next
function next(mc::MarkovChain{P,D}, state) where {P,D}
    @set mc.step.args.state = state
end

function Distributions.logpdf(mc::MarkovChain{P,D}, x::AbstractVector{X}) where {P,D,X}
    ℓ = 0.0
    for xj in x
        ℓ += logpdf(mc.step,xj)
        @set! mc.step.args.state = xj
    end
    return ℓ
end

function Base.iterate(r::MarkovChainRand{P,D}) where {P,D}
    state = rand(r.mc.step).next.state
    return (state, state)
end

function Base.iterate(r::MarkovChainRand{P,D}, state::NamedTuple) where {P,D}
    step = next(r.mc,state).step
    newstate = rand(step).next.state
    return (newstate, newstate)
end

Base.rand(mc::MarkovChain) = MarkovChainRand(mc)

# EXAMPLE

# pars = (σ = 1.2,)
# s0 = (x=0.0,)

# mstep = @model pars,state begin
#     σ = pars.σ
#     x0 = state.x
#     x ~ Normal(x0, σ)
#     next = (pars=pars, state=(x=x,))
# end

# m = @model s0 begin
#     σ ~ Exponential()
#     pars = (σ=σ,)
#     x ~ MarkovChain(pars, mstep(pars=pars, state=s0))
# end

# r = rand(m(s0=s0))

# for (n,s) in enumerate(r.x)
#     n > 10 && break
#     println(s)
# end


text:
"""
# Description

Rewrite an expression to escape all symbols except those in a specified ignore
list.

See docs/escaping.md for details.

# Arguments

1. `e::Any`: An expression.
2. `ignore::Tuple`:

# Return Values

1. `e::Any`:

# Examples

```
julia> escape_symbols(:(a + b), (:a, :b))
:(($(Expr(:escape, :+)))(a, b))
```
"""
function escape_symbols(
    @nospecialize(e::Any),
    ignore::Tuple=(:uses_default_lifting, ),
)::Any
    if isa(e, Expr)
        Expr(
            e.head,
            map(a -> escape_symbols(a, ignore), e.args)...,
        )
    elseif isa(e, Symbol)
        if e in ignore
            e
        else
            esc(e)
        end
    else
        e
    end
end


text:
using ChemistryFeaturization.Codec

@testset "Codecs" begin
    codec_tests = ["OneHotOneCold_tests", "SimpleCodec_tests", "DirectCodec_tests"]
    for t in codec_tests
        tp = abspath(testdir, "codecs", "$(t).jl")
        include(tp)
    end
end


text:
type Context
    request::Request
    response::Response
    data::Dict{Any, Any}
    Context(req::Request) = new(req, Response(), Dict{Any, Any}())
end

"""
# Response(::Context)

converts the context to a response
"""
Response(ctx::Context) = ctx.response