text:
"""
    ImplicitForm{T<:Number}

Type defining an implicit Bernstein form.

### Fields

- `array` -- the vector of vectors holding the Bernstein coefficients
- `dim`   -- an integer representing the ambient dimension
"""
struct ImplicitForm{N<:Number}
    array::Vector{Vector{N}}
    dim::Int64
    ncoeffs::Int64
end

include("monomials.jl")

ext: jl


text:
function myclamp(x)
	return x < 0 ? 0 : x > 1 ? 1 : x
end

myclamp(0.33)

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
struct LazyHelp
    o::Any
    keys::Tuple{Vararg{String}}
end

LazyHelp(o) = LazyHelp(o, ())
LazyHelp(o, k::AbstractString) = LazyHelp(o, (k,))
LazyHelp(o, k1::AbstractString, k2::AbstractString) = LazyHelp(o, (k1, k2))
LazyHelp(o, s::Symbol) = LazyHelp(o, (string(s),))

function show(io::IO, st, h::LazyHelp)
    py"""
    import inspect
    def get_signature(f):
        try:
            return str(inspect.signature(f))
        except ValueError:
            return ""
    """

    o = h.o
    for k in h.keys
        o = o[k]
    end

    fname = hasproperty(o, "__name__") ? o.__name__ : ""
    sig = hasproperty(o, "__call__") ? py"get_signature"(o) : ""
    fdoc = hasproperty(o, "__doc__") ? o.__doc__ : ""

    docstr = """
        $(fname)$(sig)
    $(fdoc)
    """

    print(io, docstr)
end

Base.show(io::IO, h::LazyHelp) = show(io, "text/plain", h)

function Base.Docs.catdoc(hs::LazyHelp...)
    Base.Docs.Text() do io
        for h in hs
            show(io, MIME"text/plain"(), h)
        end
    end
end

macro pyfunc(_pymod, _func)
	pymod = _pymod |> esc
	func = _func |> esc
	quote
		@doc LazyHelp($(pymod), nameof($(func)))
			function $(func)(args...; kwargs...)
				getproperty($(pymod), nameof($(func)))(args...; kwargs...)
			end
	end
end

macro pyclass(_pymod, _class)
	pymod = _pymod |> esc
	class = _class |> esc
	quote
		@doc LazyHelp($(pymod), nameof($(class)))
			struct $(class)
			    pyobj::$(PyObject)
			    $(class)(args...; kwargs...) = new(getproperty($(pymod), nameof($class))(args...; kwargs...))
			end

		PyCall.PyObject(t::$(class)) = t.pyobj

		function Base.propertynames(t::$(class))
		    propertynames(getfield(t, :pyobj))
		end

		function Base.getproperty(t::$(class), s::Symbol)
		    if s ∈ fieldnames($(class))
		        return getfield(t, s)
		    else
		        return getproperty(getfield(t, :pyobj), s)
		    end
		end
	end
end


text:
function force(x::Float64, potential::Potential)
    f(y) = -derivative(potential.f,y)
    f(x)
end

function friction(z::Float64, thermo::Thermostat)
    g = y -> derivative(thermo.distribution,y)
    g(z)/thermo.distribution(z)
end

function forcederivative(potential::Potential)
    force(x) = -derivative(potential.f,x)
    fprime(x) = derivative(force,x)
end

"""
Equations of motion fot a general potential coupled to any thermostat
"""
function DDfield(r::Vector{Float64}, potential::Potential, beta::Float64, thermo::Thermostat)

    (q, p, z) = r

    f = x-> force(x, potential)

    dq_dt = p
    dp_dt = f(q) + friction(z,thermo)*p/beta
    dz_dt = p^2. - 1.0/beta

    [dq_dt, dp_dt, dz_dt]

end

"""
    jacobian(r, potential, beta, thermo)
Returns the jacobian of the DDfield evaluated at the point `r`
"""
function jacobian(r::Vector{Float64}, potential::Potential, beta::Float64, thermo::Thermostat)

    (q,p,z) = r

    fprime = forcederivative(potential)

    distprime(x) = derivative(thermo.distribution,x)
    distprime2(x) = derivative(distprime,x)

    J = [0. 1. 0.; fprime(q)  friction(z, thermo)/beta p/beta*(distprime2(z)/thermo.distribution(z) - friction(z, thermo)^2.); 0. 2.0*p 0.]

end

"""
Evaluates the variational field associated to the DDfield at the extended vector `r` (q,p,z  and 9 more entries representing the coordinates of three tangent vectors)
"""
function variationalDDfield(r_and_phi::Vector{Float64}, potential::Potential, beta::Float64, thermo::Thermostat)

    v = DDfield(r_and_phi[1:3], potential, beta, thermo)
    J = jacobian(r_and_phi[1:3], potential, beta, thermo)
    rmatrix = reshape(r_and_phi[4:end],3,3)
    DPhi = J*rmatrix'

    return append!(v, DPhi'[:])
end

"""
Given a 3x3 matrix  orthogonalize its columns with the Gram-Schmidt procedure
"""
function gramschmidt(u::Matrix{Float64})
    w = eye(3)
    w[:,1] = u[:,1];
    v1 = w[:,1]/norm(w[:,1])
    w[:,2] = u[:,2] - dot(u[:,2],v1)*v1;
    v2 = w[:,2]/norm(w[:,2]);
    w[:,3] = (u[:,3] - dot(u[:,3],v2)*v2 - dot(u[:,3],v1)*v1)

    return w
end

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
Pkg.update()
using PkgDev
PkgDev.register("UnitRoot")
PkgDev.tag("UnitRoot")
PkgDev.publish()


text:
using MetroMapRuby
using Test
using JuMP
using GLPKMathProgInterface
using MathProgBase

coord1 = EuclideanCoordinate(1, 2)
coord2 = EuclideanCoordinate(1, 3)
station1 = EuclideanStation(1, coord1, "test", false)
station2 = EuclideanStation(2, coord2, "test", false)
station3 = EuclideanStation(3, coord2, "test", false)
station4 = EuclideanStation(4, coord2, "test", false)
station5 = EuclideanStation(5, coord2, "test", false)
station6 = EuclideanStation(6, coord2, "test", false)
station7 = EuclideanStation(7, coord2, "test", false)
station8 = EuclideanStation(8, coord2, "test", false)
station9 = EuclideanStation(9, coord2, "test", false)
line1 = Line("S1")
line2 = Line("S2")
line3 = Line("S3")
line4 = Line("S4")
edge = ProcessedEdge(station1, station2, line1, 0, 1, false)
edge2 = ProcessedEdge(station2, station3, line1, 0, 1, false)
edge3 = ProcessedEdge(station2, station4, line2, 2, 1, false)
edge4 = ProcessedEdge(station5, station2, line2, 2, 1, false)
edge5 = ProcessedEdge(station6, station2, line3, 1, 1, false)
edge6 = ProcessedEdge(station2, station7, line3, 1, 1, false)
edge7 = ProcessedEdge(station8, station2, line4, 3, 1, false)
edge8 = ProcessedEdge(station2, station9, line4, 3, 1, false)
transit_map = InputGraph([station1, station2, station3, station4, station5,
    station6, station7, station8, station9],
    [edge, edge2, edge3, edge4, edge5, edge6, edge7, edge8],
    [line1, line2, line3, line4])

# 9     4    7
#   x   ^   x
#     x | x
# 1 ->  2 -> 3
#     x ^ x
#   x   |   x
# 6     5    8

result = MetroMapRuby.optimize(GLPKSolverMIP(msg_lev = 3), transit_map, false)

@test result.stations[1].coordinate.y ≈ result.stations[2].coordinate.y
@test result.stations[2].coordinate.y ≈ result.stations[3].coordinate.y
@test result.stations[2].coordinate.x ≈ result.stations[4].coordinate.x
@test result.stations[2].coordinate.x ≈ result.stations[5].coordinate.x
@test result.stations[4].coordinate.y > result.stations[5].coordinate.y

@test result.stations[6].coordinate.x ≈ result.stations[9].coordinate.x
@test result.stations[6].coordinate.x ≈ result.stations[1].coordinate.x
@test result.stations[7].coordinate.x ≈ result.stations[8].coordinate.x
@test result.stations[7].coordinate.x ≈ result.stations[3].coordinate.x
@test result.stations[4].coordinate.y ≈ result.stations[9].coordinate.y
@test result.stations[9].coordinate.y ≈ result.stations[7].coordinate.y
@test result.stations[6].coordinate.y ≈ result.stations[8].coordinate.y
@test result.stations[5].coordinate.y ≈ result.stations[6].coordinate.y

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
module AutoDiff

using CUDA
using CUDA.CUSPARSE
using ForwardDiff
using KernelAbstractions
using SparseArrays
using TimerOutputs
using SparseDiffTools
using ..ExaPF: Spmat, xzeros

import Base: show

"""
    AbstractJacobian

Automatic differentiation for the compressed Jacobians of the
constraints `g(x,u)` with respect to the state `x` and the control `u`
(here called design).

TODO: Use dispatch to unify the code of the state and control Jacobian.
This is currently not done because the abstraction of the indexing is not yet resolved.

"""
abstract type AbstractJacobian end

function _init_seed!(t1sseeds, coloring, ncolor, nmap)
    t1sseedvec = zeros(Float64, ncolor)
    @inbounds for i in 1:nmap
        for j in 1:ncolor
            if coloring[i] == j
                t1sseedvec[j] = 1.0
            end
        end
        t1sseeds[i] = ForwardDiff.Partials{ncolor, Float64}(NTuple{ncolor, Float64}(t1sseedvec))
        t1sseedvec .= 0
    end
end

"""
    StateJacobian

Creates an object for the state Jacobian

* `J::SMT`: Sparse uncompressed Jacobian to be used by linear solver. This is either of type `SparseMatrixCSC` or `CuSparseMatrixCSR`.
* `compressedJ::MT`: Dense compressed Jacobian used for updating values through AD either of type `Matrix` or `CuMatrix`.
* `coloring::VI`: Row coloring of the Jacobian.
* `t1sseeds::VP`: The seeding vector for AD built based on the coloring.
* `t1sF::VD`: Output array of active (AD) type.
* `x::VT`: Input array of passive type. This includes both state and control.
* `t1sx::VD`: Input array of active type.
* `map::VI`: State and control mapping to array `x`
* `varx::SubT`: View of `map` on `x`
* `t1svarx::SubD`: Active (AD) view of `map` on `x`
"""
struct StateJacobian{VI, VT, MT, SMT, VP, VD, SubT, SubD} <: AbstractJacobian
    J::SMT
    compressedJ::MT
    coloring::VI
    t1sseeds::VP
    t1sF::VD
    x::VT
    t1sx::VD
    map::VI
    # Cache views on x and its dual vector to avoid reallocating on the GPU
    varx::SubT
    t1svarx::SubD
    function StateJacobian(F, v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus)
        nv_m = size(v_m, 1)
        nv_a = size(v_a, 1)
        if F isa Array
            VI = Vector{Int}
            VT = Vector{Float64}
            MT = Matrix{Float64}
            SMT = SparseMatrixCSC
            A = Vector
        elseif F isa CuArray
            VI = CuVector{Int}
            VT = CuVector{Float64}
            MT = CuMatrix{Float64}
            SMT = CuSparseMatrixCSR
            A = CuVector
        else
            error("Wrong array type ", typeof(F))
        end

        mappv = [i + nv_m for i in pv]
        mappq = [i + nv_m for i in pq]
        # Ordering for x is (θ_pv, θ_pq, v_pq)
        map = VI(vcat(mappv, mappq, pq))
        nmap = size(map,1)

        # Used for sparsity detection with randomized inputs
        function residual_jacobian(V, Ybus, pv, pq)
            n = size(V, 1)
            Ibus = Ybus*V
            diagV       = sparse(1:n, 1:n, V, n, n)
            diagIbus    = sparse(1:n, 1:n, Ibus, n, n)
            diagVnorm   = sparse(1:n, 1:n, V./abs.(V), n, n)

            dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm
            dSbus_dVa = 1im * diagV * conj(diagIbus - Ybus * diagV)

            j11 = real(dSbus_dVa[[pv; pq], [pv; pq]])
            j12 = real(dSbus_dVm[[pv; pq], pq])
            j21 = imag(dSbus_dVa[pq, [pv; pq]])
            j22 = imag(dSbus_dVm[pq, pq])

            J = [j11 j12; j21 j22]
        end

        # Need a host arrays for the sparsity detection below
        spmap = Vector(map)
        hybus_re = Spmat{Vector{Int}, Vector{Float64}}(ybus_re)
        hybus_im = Spmat{Vector{Int}, Vector{Float64}}(ybus_im)
        n = nv_a
        Yre = SparseMatrixCSC{Float64,Int64}(n, n, hybus_re.colptr, hybus_re.rowval, hybus_re.nzval)
        Yim = SparseMatrixCSC{Float64,Int64}(n, n, hybus_im.colptr, hybus_im.rowval, hybus_im.nzval)
        Y = Yre .+ 1im .* Yim
        # Randomized inputs
        Vre = Float64.([i for i in 1:n])
        Vim = Float64.([i for i in n+1:2*n])
        V = Vre .+ 1im .* Vim
        J = residual_jacobian(V, Y, pv, pq)
        coloring = VI(matrix_colors(J))
        ncolor = size(unique(coloring),1)
        if F isa CuArray
            J = CuSparseMatrixCSR(J)
        end
        t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N
        x = VT(zeros(Float64, nv_m + nv_a))
        t1sx = A{t1s{ncolor}}(x)
        t1sF = A{t1s{ncolor}}(zeros(Float64, nmap))

        t1sseeds = A{ForwardDiff.Partials{ncolor,Float64}}(undef, nmap)
        _init_seed!(t1sseeds, coloring, ncolor, nmap)

        compressedJ = MT(zeros(Float64, ncolor, nmap))
        # Views
        varx = view(x, map)
        t1svarx = view(t1sx, map)
        VP = typeof(t1sseeds)
        VD = typeof(t1sx)
        return new{VI, VT, MT, SMT, VP, VD, typeof(varx), typeof(t1svarx)}(
            J, compressedJ, coloring, t1sseeds, t1sF, x, t1sx, map, varx, t1svarx
        )
    end
end

"""
    ControlJacobian

Creates an object for the control Jacobian.

* `J::SMT`: Sparse uncompressed Jacobian to be used by linear solver. This is either of type `SparseMatrixCSC` or `CuSparseMatrixCSR`.
* `compressedJ::MT`: Dense compressed Jacobian used for updating values through AD either of type `Matrix` or `CuMatrix`.
* `coloring::VI`: Row coloring of the Jacobian.
* `t1sseeds::VP`: The seeding vector for AD built based on the coloring.
* `t1sF::VD`: Output array of active (AD) type.
* `x::VT`: Input array of passive type. This includes both state and control.
* `t1sx::VD`: Input array of active type.
* `map::VI`: State and control mapping to array `x`
* `varx::SubT`: View of `map` on `x`
* `t1svarx::SubD`: Active (AD) view of `map` on `x`
"""
struct ControlJacobian{VI, VT, MT, SMT, VP, VD, SubT, SubD} <: AbstractJacobian
    J::SMT
    compressedJ::MT
    coloring::VI
    t1sseeds::VP
    t1sF::VD
    x::VT
    t1sx::VD
    map::VI
    # Cache views on x and its dual vector to avoid reallocating on the GPU
    varx::SubT
    t1svarx::SubD
    function ControlJacobian(F, v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus)
        nv_m = size(v_m, 1)
        nv_a = size(v_a, 1)
        npbus = size(pinj, 1)
        nref = size(ref, 1)
        if F isa Array
            VI = Vector{Int}
            VT = Vector{Float64}
            MT = Matrix{Float64}
            SMT = SparseMatrixCSC
            A = Vector
        elseif F isa CuArray
            VI = CuVector{Int}
            VT = CuVector{Float64}
            MT = CuMatrix{Float64}
            SMT = CuSparseMatrixCSR
            A = CuVector
        else
            error("Wrong array type ", typeof(F))
        end

        mappv =  [i + nv_a for i in pv]
        map = VI(vcat(ref, mappv, pv))
        nmap = size(map,1)

        # Used for sparsity detection with randomized inputs
        function residual_jacobian(V, Ybus, pinj, qinj, ref, pv, pq)
            n = size(V, 1)
            Ibus = Ybus*V
            diagV       = sparse(1:n, 1:n, V, n, n)
            diagIbus    = sparse(1:n, 1:n, Ibus, n, n)
            diagVnorm   = sparse(1:n, 1:n, V./abs.(V), n, n)

            dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm
            dSbus_dpbus = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm

            j11 = real(dSbus_dVm[[pv; pq], [ref; pv; pv]])
            j21 = imag(dSbus_dVm[pq, [ref; pv; pv]])
            J = [j11; j21]
        end

        # Need a host arrays for the sparsity detection below
        spmap = Vector(map)
        hybus_re = Spmat{Vector{Int}, Vector{Float64}}(ybus_re)
        hybus_im = Spmat{Vector{Int}, Vector{Float64}}(ybus_im)
        n = nv_a
        Yre = SparseMatrixCSC{Float64,Int64}(n, n, hybus_re.colptr, hybus_re.rowval, hybus_re.nzval)
        Yim = SparseMatrixCSC{Float64,Int64}(n, n, hybus_im.colptr, hybus_im.rowval, hybus_im.nzval)
        Y = Yre .+ 1im .* Yim
        # Randomized inputs
        Vre = Float64.([i for i in 1:n])
        Vim = Float64.([i for i in n+1:2*n])
        V = Vre .+ 1im .* Vim
        J = residual_jacobian(V, Y, pinj, qinj, ref, pv, pq)
        coloring = VI(matrix_colors(J))
        ncolor = size(unique(coloring),1)
        if F isa CuArray
            J = CuSparseMatrixCSR(J)
        end
        t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N
        x = xzeros(VT, npbus + nv_a)
        t1sx = A{t1s{ncolor}}(x)
        t1sF = A{t1s{ncolor}}(zeros(Float64, length(F)))
        t1sseeds = A{ForwardDiff.Partials{ncolor,Float64}}(undef, nmap)
        _init_seed!(t1sseeds, coloring, ncolor, nmap)

        compressedJ = MT(zeros(Float64, ncolor, length(F)))
        # Views
        varx = view(x, map)
        t1svarx = view(t1sx, map)
        VP = typeof(t1sseeds)
        VD = typeof(t1sx)
        return new{VI, VT, MT, SMT, VP, VD, typeof(varx), typeof(t1svarx)}(
            J, compressedJ, coloring, t1sseeds, t1sF, x, t1sx, map, varx, t1svarx
        )
    end
end

"""
    seed_kernel_cpu!

Seeding on the CPU, not parallelized.

"""
function seed_kernel_cpu!(
    duals::AbstractArray{ForwardDiff.Dual{T,V,N}}, x,
    seeds::AbstractArray{ForwardDiff.Partials{N,V}}
) where {T,V,N}
    for i in 1:size(duals,1)
        duals[i] = ForwardDiff.Dual{T,V,N}(x[i], seeds[i])
    end
end

"""
    seed_kernel_gpu!

Seeding on GPU parallelized over the `ncolor` number of duals

"""
function seed_kernel_gpu!(
    duals::AbstractArray{ForwardDiff.Dual{T,V,N}}, x,
    seeds::AbstractArray{ForwardDiff.Partials{N,V}}
) where {T,V,N}
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:size(duals,1)
        duals[i] = ForwardDiff.Dual{T,V,N}(x[i], seeds[i])
    end
end

"""
    seed_kernel!

Calling the GPU seeding kernel

"""
function seed_kernel!(t1sseeds::CuVector{ForwardDiff.Partials{N,V}}, varx, t1svarx, nbus) where {N, V}
    nthreads = 256
    nblocks = div(nbus, nthreads, RoundUp)
    CUDA.@sync begin
        @cuda threads=nthreads blocks=nblocks seed_kernel_gpu!(
            t1svarx,
            varx,
            t1sseeds,
        )
    end
end

"""
    seed_kernel!(t1sseeds::Vector{ForwardDiff.Partials{N,V}}, varx, t1svarx, nbus) where {N, V}

Calling the CPU seeding kernel

"""
function seed_kernel!(t1sseeds::Vector{ForwardDiff.Partials{N,V}}, varx, t1svarx, nbus) where {N, V}
    seed_kernel_cpu!(t1svarx, varx, t1sseeds)
end

"""
    getpartials_kernel_cpu!(compressedJ, t1sF)

Extract the partials from the AutoDiff dual type on the CPU and put it in the
compressed Jacobian

"""
function getpartials_kernel_cpu!(compressedJ, t1sF)
    for i in 1:size(t1sF,1) # Go over outputs
        compressedJ[:, i] .= ForwardDiff.partials.(t1sF[i]).values
    end
end

"""
    getpartials_kernel_gpu!(compressedJ, t1sF)

Extract the partials from the AutoDiff dual type on the GPU and put it in the
compressed Jacobian

"""
function getpartials_kernel_gpu!(compressedJ, t1sF)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:size(t1sF, 1) # Go over outputs
        for j in eachindex(ForwardDiff.partials.(t1sF[i]).values)
            @inbounds compressedJ[j, i] = ForwardDiff.partials.(t1sF[i]).values[j]
        end
    end
end

"""
    getpartials_kernel!(compressedJ::CuArray{T, 2}, t1sF, nbus) where T

Calling the GPU partial extraction kernel

"""
function getpartials_kernel!(compressedJ::CuArray{T, 2}, t1sF, nbus) where T
    nthreads = 256
    nblocks = div(nbus, nthreads, RoundUp)
    CUDA.@sync begin
        @cuda threads=nthreads blocks=nblocks getpartials_kernel_gpu!(
            compressedJ,
            t1sF
        )
    end
end

"""
    getpartials_kernel!(compressedJ::Array{T, 2}, t1sF, nbus) where T

Calling the CPU partial extraction kernel

"""
function getpartials_kernel!(compressedJ::Array{T, 2}, t1sF, nbus) where T
    getpartials_kernel_cpu!(compressedJ, t1sF)
end

"""
    uncompress_kernel_gpu!(J_nzVal, J_rowPtr, J_colVal, compressedJ, coloring, nmap)

Uncompress the compressed Jacobian matrix from `compressedJ` to sparse CSR on
the GPU. Only bitarguments are allowed for the kernel.
(for GPU only) TODO: should convert to @kernel
"""
function uncompress_kernel_gpu!(J_nzVal, J_rowPtr, J_colVal, compressedJ, coloring, nmap)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:nmap
        for j in J_rowPtr[i]:J_rowPtr[i+1]-1
            @inbounds J_nzVal[j] = compressedJ[coloring[J_colVal[j]], i]
        end
    end
end

"""
    uncompress_kernel!(J::SparseArrays.SparseMatrixCSC, compressedJ, coloring)

Uncompress the compressed Jacobian matrix from `compressedJ` to sparse CSC on
the CPU.
"""
function uncompress_kernel!(J::SparseArrays.SparseMatrixCSC, compressedJ, coloring)
    # CSC is column oriented: nmap is equal to number of columns
    nmap = size(J, 2)
    @assert(maximum(coloring) == size(compressedJ,1))
    for i in 1:nmap
        for j in J.colptr[i]:J.colptr[i+1]-1
            @inbounds J.nzval[j] = compressedJ[coloring[i], J.rowval[j]]
        end
    end
end

"""
    uncompress_kernel!(J::CUDA.CUSPARSE.CuSparseMatrixCSR, compressedJ, coloring)

Uncompress the compressed Jacobian matrix from `compressedJ` to sparse CSC on
the GPU by calling the kernel [`uncompress_kernel_gpu!`](@ref).
"""
function uncompress_kernel!(J::CUDA.CUSPARSE.CuSparseMatrixCSR, compressedJ, coloring)
    # CSR is row oriented: nmap is equal to number of rows
    nmap = size(J, 1)
    nthreads = 256
    nblocks = div(nmap, nthreads, RoundUp)
    CUDA.@sync begin
        @cuda threads=nthreads blocks=nblocks uncompress_kernel_gpu!(
                J.nzVal,
                J.rowPtr,
                J.colVal,
                compressedJ,
                coloring, nmap
        )
    end
end

"""
    residual_jacobian!(arrays::StateJacobian,
                        residual_polar!,
                        v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus,
                        timer = nothing)

Update the sparse Jacobian entries using AutoDiff. No allocations are taking place in this function.

* `arrays::StateJacobian`: Factory created Jacobian object to update
* `residual_polar`: Primal function
* `v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus`: Inputs both
  active and passive parameters. Active inputs are mapped to `x` via the preallocated views.

"""
function residual_jacobian!(arrays::StateJacobian,
                             residual_polar!,
                             v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus,
                             timer = nothing)
    @timeit timer "Before" begin
        @timeit timer "Setup" begin
            nv_m = size(v_m, 1)
            nv_a = size(v_a, 1)
            nmap = size(arrays.map, 1)
            n = nv_m + nv_a
        end
        @timeit timer "Arrays" begin
            arrays.x[1:nv_m] .= v_m
            arrays.x[nv_m+1:nv_m+nv_a] .= v_a
            arrays.t1sx .= arrays.x
            arrays.t1sF .= 0.0
        end
    end
    @timeit timer "Seeding" begin
        seed_kernel!(arrays.t1sseeds, arrays.varx, arrays.t1svarx, nbus)
    end

    @timeit timer "Function" begin
        residual_polar!(
            arrays.t1sF,
            view(arrays.t1sx, 1:nv_m),
            view(arrays.t1sx, nv_m+1:nv_m+nv_a),
            ybus_re, ybus_im,
            pinj, qinj,
            pv, pq, nbus
        )
    end

    @timeit timer "Get partials" begin
        getpartials_kernel!(arrays.compressedJ, arrays.t1sF, nbus)
    end
    @timeit timer "Uncompress" begin
        uncompress_kernel!(arrays.J, arrays.compressedJ, arrays.coloring)
    end
    return nothing
end

"""
    residual_jacobian!(arrays::ControlJacobian,
                        residual_polar!,
                        v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus,
                        timer = nothing)

Update the sparse Jacobian entries using AutoDiff. No allocations are taking place in this function.

* `arrays::ControlJacobian`: Factory created Jacobian object to update
* `residual_polar`: Primal function
* `v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus`: Inputs both active and passive parameters. Active inputs are mapped to `x` via the preallocated views.

"""
function residual_jacobian!(arrays::ControlJacobian,
                             residual_polar!,
                             v_m, v_a, ybus_re, ybus_im, pinj, qinj, pv, pq, ref, nbus,
                             timer = nothing)

    @timeit timer "Before" begin
        @timeit timer "Setup" begin
            npinj = size(pinj , 1)
            nv_m = size(v_m, 1)
            nmap = size(arrays.map, 1)
            n = npinj + nv_m
        end
        @timeit timer "Arrays" begin
            arrays.x[1:nv_m] .= v_m
            arrays.x[nv_m+1:nv_m+npinj] .= pinj
            arrays.t1sx .= arrays.x
            arrays.t1sF .= 0.0
        end
    end
    @timeit timer "Seeding" begin
        seed_kernel!(arrays.t1sseeds, arrays.varx, arrays.t1svarx, nbus)
    end
    @timeit timer "Function" begin
        residual_polar!(
            arrays.t1sF,
            view(arrays.t1sx, 1:nv_m),
            v_a,
            ybus_re, ybus_im,
            view(arrays.t1sx, nv_m+1:nv_m + npinj), qinj,
            pv, pq, nbus
        )
    end

    @timeit timer "Get partials" begin
        getpartials_kernel!(arrays.compressedJ, arrays.t1sF, nbus)
    end
    @timeit timer "Uncompress" begin
        # Uncompress matrix. Sparse matrix elements have different names with CUDA
        uncompress_kernel!(arrays.J, arrays.compressedJ, arrays.coloring)
    end
end

function Base.show(io::IO, jacobian::AbstractJacobian)
    ncolor = size(unique(jacobian.coloring), 1)
    print(io, "Number of Jacobian colors: ", ncolor)
end

end


text:
export Sdd, SddMgr,
    SddLeafNode, SddInnerNode, SddLiteralNode, SddConstantNode,
    Sdd⋀Node, Sdd⋁Node,
    sdd_mgr_for,
    (¬), (∨), (∧), (⇒), (⇐), (⇔)

#############
# Trimmed Sdds
#############

"Root of the trimmed Sdd circuit node hierarchy"
abstract type Sdd <: StructLogicCircuit end

#############
# Elements and XY Partitions
#############

"Represents elements that are not yet compiled into conjunctions"
struct Element
    prime::Sdd
    sub::Sdd
end

"Represent an XY-partition that has not yet been compiled into a disjunction"
const XYPartition = Vector{Element}

"Unique nodes cache for decision nodes"
const Unique⋁Cache = Dict{XYPartition,Sdd}

"Broader definition of an XYPartition, also including vectors of nodes"

Base.hash(x::XYPartition) = mapreduce(e -> hash(e.prime, hash(e.sub)), ⊻, x)

function Base.isequal(x::XYPartition, y::XYPartition)
    length(x) != length(y) && return false
    (x == y) && return true
    # note: the rest of this function could be removed if the XYPartitions were sorted, but that is not advantageous: the sorting is too slow
    return all(x) do e1
        any(y) do e2
            e1.prime === e2.prime && e1.sub === e2.sub
        end
    end
end

#############
# Apply cache
#############

"Representation of the arguments of an Apply call"
struct ApplyArgs
    a1::Sdd
    a2::Sdd
end

Base.hash(aa::ApplyArgs) = hash(aa.a1) ⊻ hash(aa.a2)
Base.isequal(x::ApplyArgs,y::ApplyArgs) =
    (x.a1 === y.a1 && x.a2 === y.a2 ) || (x.a1 === y.a2 && x.a2 === y.a1)

"Apply cache for the result of conjunctions and disjunctions"
const ApplyCache = Dict{ApplyArgs,Sdd}

#############
# Trimmed Sdd managers
#############

"Root of the SDD manager node hierarchy"
abstract type SddMgr <: Vtree end

"SDD manager inner vtree node for trimmed SDD nodes"
mutable struct SddMgrInnerNode <: SddMgr

    left::SddMgr
    right::SddMgr

    parent::Union{SddMgrInnerNode, Nothing}

    variables::BitSet
    unique⋁cache::Unique⋁Cache

    conjoin_cache::ApplyCache

    SddMgrInnerNode(left::SddMgr, right::SddMgr) where D = begin
        @assert isdisjoint(variables(left), variables(right))
        this = new(left, right, nothing,
            union(variables(left), variables(right)),
            Unique⋁Cache(), ApplyCache())
        @assert left.parent === nothing
        left.parent = this
        @assert right.parent === nothing
        right.parent = this
        this
    end

end

"SDD manager leaf vtree node for trimmed SDD nodes"
mutable struct SddMgrLeafNode <: SddMgr

    var::Var
    parent::Union{SddMgrInnerNode, Nothing}

    positive_literal::Sdd
    negative_literal::Sdd

    SddMgrLeafNode(v::Var) = begin
        this = new(v, nothing)
        this.positive_literal = construct_new_sdd_literal(this, var2lit(v))
        this.negative_literal = construct_new_sdd_literal(this, -var2lit(v))
        this
    end

end

@inline NodeType(::SddMgrLeafNode) = Leaf()
@inline NodeType(::SddMgrInnerNode) = Inner()

#############
# Sdd nodes
#############

"A SDD logical leaf node"
abstract type SddLeafNode <: Sdd end

"A SDD logical inner node"
abstract type SddInnerNode <: Sdd end

"A SDD logical literal leaf node, representing the positive or negative literal of its variable"
mutable struct SddLiteralNode <: SddLeafNode
    literal::Lit
    vtree::SddMgrLeafNode
end

"""
A SDD logical constant leaf node, representing true or false.
These are the only structured nodes that don't have an associated vtree node (cf. SDD file format)
"""
struct SddConstantNode <: SddLeafNode
    constant::Bool
end

"Canonical true Sdd node"
const true_sdd = SddConstantNode(true)

"Canonical false Sdd node"
const false_sdd = SddConstantNode(false)

"A SDD logical conjunction node"
mutable struct Sdd⋀Node <: SddInnerNode
    prime::Sdd
    sub::Sdd
    vtree::SddMgrInnerNode
end

"A SDD logical disjunction node"
mutable struct Sdd⋁Node <: SddInnerNode
    children::Vector{Sdd⋀Node}
    vtree::SddMgrInnerNode
    negation::Sdd⋁Node
    Sdd⋁Node(ch,v) = new(ch, v) # leave negation uninitialized
    Sdd⋁Node(ch,v,neg) = new(ch, v, neg)
end

#####################
# traits
#####################

@inline GateType(::Type{<:SddLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:SddConstantNode}) = ConstantGate()
@inline GateType(::Type{<:Sdd⋀Node}) = ⋀Gate()
@inline GateType(::Type{<:Sdd⋁Node}) = ⋁Gate()

#####################
# Constructor
#####################

SddMgr(v::Var) = SddMgrLeafNode(v)
SddMgr(left::SddMgr, right::SddMgr) = SddMgrInnerNode(left, right)

"Obtain an SDD manager that can support compiling the given circuit"
sdd_mgr_for(c::Sdd) = mgr(c)
sdd_mgr_for(c::StructLogicCircuit) =
    (vtree(c) isa SddMgr) ? vtree(c) : SddMgr(vtree(c))
sdd_mgr_for(c::LogicCircuit) =
    SddMgr(SddMgr.(Var.(variables(c))), :balanced) # this could be "smarter" finding a good vtree

"Helper function to construct new SDD literal objects in SddMgrLeafNode constructor"
function construct_new_sdd_literal(n::SddMgrLeafNode, l::Lit)::SddLiteralNode
    @assert n.var == lit2var(l) "Cannot compile literal $l respecting vtree leaf for variable $(n.var)"
    SddLiteralNode(l,n)
end

# logical syntactic sugar

"Logical negation"
(¬)(x::Sdd) = !x
"Logical disjunction"
(∨)(x::Sdd, y::Sdd) = x | y
"Logical conjunction"
(∧)(x::Sdd, y::Sdd) = x & y
"Material logical implication"
(⇒)(x::Sdd, y::Sdd) = ¬x ∨ y
"Material logical implication (reverse)"
(⇐)(x::Sdd, y::Sdd) = (y ⇒ x)
"Bidirectional logical implication"
(⇔)(x::Sdd, y::Sdd) = (x ⇒ y) ∧ (y ⇒ x)

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using DFWannier
using Optim
using DataStructures
function optimize_l()
    T = Float64

    job = ldj("/home/ponet/HfO2/NSOC/")
    wanjob = add_wan_data(job)

    k_points = qe_read_ks_from_bands_file("/home/ponet/HfO2/SOC/bands.out")[2]
    test_eig, , , = calc_observables(wanjob.structure, kpoints[1])

    bands = qe_read_bands_file("/home/ponet/HfO2/SOC/bands.out", T)[b:e]

    at_ids = getfield.(unique(wanjob.structure.atoms), :id)

    lsocs = OrderedDict{Symbol, T}()
    for id in at_ids
        lsocs[id] = zero(T)
    end
    function f(lambdas)
        for (key, l) in zip(keys(lsocs), lambdas)
            lsocs[key] = l
        end
        set_soc!(wanjob.structure, lsocs...)

        energy, , , = calc_observables(wanjob.structure, k_points)

        out_diff = 0.0
        for i=1:length(bands)
            calcbands = [val[i] for val in energy]
            for (calc, exact) in zip(calcbands[80:120], bands[i].eigvals[80:120])
                out_diff += abs(calc-exact)
            end
        end
        return out_diff
    end
    return optimize(f, [lambdaHf, lambdaO], show_trace = true, show_every = 1, extended_trace=true)
end
results = optimize_l()
println(results)


text:
module MICE

using DataFrames
using GLM: fit, LinearModel, predict
using Statistics: mean
using StatsModels: term, schema, modelcols, apply_schema
import CSV
using StatsBase: sample

export mice

function mice(data::DataFrame; max_iterations = 5)
    data_missings = ismissing.(data)

    # Do naive imputation using means
    data_imputed = naive_imputation(data, mean)

    # Do better imputation using LMs
    for iteration in 1:max_iterations

        # TODO: Add early stopping when abs(max(X_t - X_{t-1}))/abs(max(X[known_vals])) < tol

        data_imputed_temp = copy(data_imputed)
        for response_variable in names(data)
            # The explanatory variables are all the other variables
            explanatory_variables = filter(n -> n != response_variable, names(data))

            # Fit a linear model with an intercept and each explanatory variable
            f = term(response_variable) ~ term(1) + foldl(+, term.(explanatory_variables))

            predicted_values = predictive_mean_matching(data_imputed, data_missings, f)

            # # Fit the model
            # # @info data_train
            # # @info f
            # # CSV.write("df.csv", data_train)
            # model = fit(LinearModel, f, data_train)
            # # Predict
            # predicted_values::Vector = predict(model, data_imputed)

            # Replace missing values with the new imputed values
            data_imputed_temp[:, response_variable] = predicted_values
        end
        data_imputed = copy(data_imputed_temp)
    end
    return data_imputed
end

"""
Perform a naïve imputation of the missing values in `data` using the function `f`.
"""
function naive_imputation(data, f = mean)
    data_imputed = copy(data)
    for n in names(data)
        # Calculate the average of the non-missing values
        imputed_value = f(skipmissing(data_imputed[:, n]))
        # Replace the missing values with the imputed value
        data_imputed[:, n] = coalesce.(data_imputed[:, n], imputed_value)
    end
    return data_imputed
end

function predictive_mean_matching(data, data_missings, formula; n = 5)
    formula_data = apply_schema(formula, schema(formula, data))
    formula_missings = apply_schema(formula, schema(formula, data_missings))

    response_variable = modelcols(formula_data.lhs, data)
    response_missings = modelcols(formula_missings.lhs, data_missings)

    # Only fit the model to rows where the response variable is not missing
    data_train = data[.!response_missings, :]

    @info formula
    CSV.write("df.csv", data_train)
    model = fit(LinearModel, formula_data, data_train)
    predicted_values = predict(model, data)

    output = zeros(length(response_variable))
    for (i, is_missing) in enumerate(response_missings)
        if is_missing
            # Calculate the distance between the predicted value being imputed and the predicted values of the other points
            donor_candidates_distances = abs.(predicted_values[i] .- predicted_values)
            # Filter to only the candiates where the true value is known
            donor_candidates = filter((j, v) -> !response_missings[j], enumerate(donor_candidates_distances))
            donor_candidates = sort(donor_candidates, by = (j, v) -> v)[2:(n + 1)]
            # Choose one of the candidates to be the donor
            donor_index, _ = sample(donor_candidates)
            # Set the value being imputed to the donor value
            output[i] = response_column[donor_index]
        else
            output[i] = response_column[i]
        end
    end
    return output
end

end # module


text:
using Test

using ScottishTaxBenefitModel
using .ModelHousehold: Household, Person, People_Dict, default_bu_allocation
using .Pensions: impute_employer_pension!
using .Definitions

@testset "Employer's Pension Imputation" begin

    hh = ExampleHouseholdGetter.get_household( "mel_c2" )
    cpl= get_example( cpl_w_2_children_hh )
    head = get_head( cpl )

    @test head.employment_status == Full_time_Employee
    @test get(head.income,pension_contributions_employer,0.0) == 0.0
    head.income[wages] = 1_000;
    impute_employer_pension!( head )
    @test head.income[wages]*0.3 > head.income[pension_contributions_employer]

end

@testset "Pension correction test on full dataset" begin
    nhhs,npeople = init_data( reset = true, Settings() )
    itsys_scot :: IncomeTaxSys = get_tax( scotland = true )

    for hno in 1:nhhs

    end
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
const LimbicSystem = BilateralStructure(:limbic_system)

const LimbicLobe = BilateralStructure(:limbic_lobe)

const Parahippocampus = BilateralStructure(:parahippocampus)

const SeptalNuclei = BilateralStructure(:septal_nuclei)

const NucleusAccumbens = BilateralStructure(:nucleus_accumbens)

const PapezCircuit = BilateralStructure(:papez_circuit)


text:
export InterpretedSystem, InterpretedHomotopy

"""
    InterpretedSystem <: AbstractSystem

An [`AbstractSystem`](@ref) which automatically generates a program for the
fast evaluation of `F` and its Jacobian. The program is however, not compiled
but rather interpreted. See also [`CompiledSystem`](@ref).

    InterpretedSystem(F::System; optimizations = true)

Construct an `InterpretedSystem` from the given [`System`](@ref) `F`.
If `optimizations = true` then [`optimize`](@ref) is called on `F` before compiling.
"""
mutable struct InterpretedSystem{T₁,T₂} <: AbstractSystem
    system::System
    eval_interpreter::Interpreter{T₁}
    eval_interpreter_cache::InterpreterCache{Vector{ComplexF64},Nothing}
    eval_interpreter_cache_ext::InterpreterCache{Vector{ComplexDF64},Nothing}
    jac_interpreter::Interpreter{T₂}
    jac_interpreter_cache::InterpreterCache{Vector{ComplexF64},Nothing}
    taylor_caches::Tuple{
        InterpreterCache{TaylorVector{2,ComplexF64},Nothing},
        InterpreterCache{TaylorVector{3,ComplexF64},Nothing},
        InterpreterCache{TaylorVector{4,ComplexF64},Nothing},
    }
    arb_eval_interpreter_cache::Union{Nothing,InterpreterCache{ArbRefVector,Arb}}
    arb_jac_interpreter_cache::Union{Nothing,InterpreterCache{ArbRefVector,Arb}}
    acb_eval_interpreter_cache::Union{Nothing,InterpreterCache{AcbRefVector,Acb}}
    acb_jac_interpreter_cache::Union{Nothing,InterpreterCache{AcbRefVector,Acb}}
end

function InterpretedSystem(F::System; optimizations::Bool = true)
    F = optimizations ? optimize(F) : F
    eval_interpreter = interpreter(F)
    jac_interpreter = jacobian_interpreter(F)
    N = cache_min_length(eval_interpreter)
    taylor_caches = (
        InterpreterCache(TaylorVector{2}(ComplexF64, N)),
        InterpreterCache(TaylorVector{3}(ComplexF64, N)),
        InterpreterCache(TaylorVector{4}(ComplexF64, N)),
    )

    InterpretedSystem(
        F,
        eval_interpreter,
        InterpreterCache(ComplexF64, eval_interpreter),
        InterpreterCache(ComplexDF64, eval_interpreter),
        jac_interpreter,
        InterpreterCache(ComplexF64, jac_interpreter),
        taylor_caches,
        nothing,
        nothing,
        nothing,
        nothing,
    )
end

Base.size(F::InterpretedSystem) = size(F.system)
variables(F::InterpretedSystem) = variables(F.system)
parameters(F::InterpretedSystem) = parameters(F.system)
variable_groups(F::InterpretedSystem) = variable_groups(F.system)
Base.:(==)(F::InterpretedSystem, G::InterpretedSystem) = F.system == G.system

function Base.show(io::IO, F::InterpretedSystem)
    print(io, "Interpreted: ")
    show(io, F.system)
end

(F::InterpretedSystem)(x, p = nothing) = F.system(x, p)
function evaluate!(u, F::InterpretedSystem, x, p = nothing)
    execute!(u, F.eval_interpreter, x, p, F.eval_interpreter_cache)
end
function evaluate!(u, F::InterpretedSystem, x::AbstractVector{ComplexDF64}, p = nothing)
    execute!(u, F.eval_interpreter, x, p, F.eval_interpreter_cache_ext)
end
function evaluate_and_jacobian!(u, U, F::InterpretedSystem, x, p = nothing)
    execute!(u, U, F.jac_interpreter, x, p, F.jac_interpreter_cache)
    nothing
end
function jacobian!(U, F::InterpretedSystem, x, p = nothing)
    evaluate_and_jacobian!(nothing, U, F, x, p)
    U
end

for M = 1:3
    @eval function taylor!(
        u::AbstractVector,
        Order::Val{$M},
        F::InterpretedSystem,
        x,
        p = nothing,
    )
        execute!(u, Order, F.eval_interpreter, x, p, F.taylor_caches[$M])
        u
    end
end

# Arb
function evaluate!(
    u,
    F::InterpretedSystem,
    x::AbstractArray{<:Union{Arb,ArbRef}},
    p = nothing,
)
    if isnothing(F.arb_eval_interpreter_cache)
        F.arb_eval_interpreter_cache = InterpreterCache(Arb, F.eval_interpreter; prec = 128)
    end
    if u isa AbstractArray{<:Union{Arb,ArbRef}}
        set_arb_precision!(
            F.arb_eval_interpreter_cache,
            max(precision(first(u)), precision(first(x))),
        )
    else
        set_arb_precision!(F.arb_eval_interpreter_cache, precision(first(x)))
    end
    execute!(u, F.eval_interpreter, x, p, F.arb_eval_interpreter_cache)
end
function evaluate!(
    u,
    F::InterpretedSystem,
    x::AbstractArray{<:Union{Acb,AcbRef}},
    p = nothing,
)
    if isnothing(F.acb_eval_interpreter_cache)
        F.acb_eval_interpreter_cache = InterpreterCache(Acb, F.eval_interpreter; prec = 128)
    end
    if u isa AbstractArray{<:Union{Acb,AcbRef}}
        set_arb_precision!(
            F.acb_eval_interpreter_cache,
            max(precision(first(u)), precision(first(x))),
        )
    else
        set_arb_precision!(F.acb_eval_interpreter_cache, precision(first(x)))
    end
    execute!(u, F.eval_interpreter, x, p, F.acb_eval_interpreter_cache)
end

function evaluate_and_jacobian!(
    u,
    U,
    F::InterpretedSystem,
    x::AbstractArray{<:Union{Arb,ArbRef}},
    p = nothing,
)
    if isnothing(F.arb_eval_interpreter_cache)
        F.arb_jac_interpreter_cache = InterpreterCache(Arb, F.jac_interpreter; prec = 128)
    end
    if U isa AbstractArray{<:Union{Arb,ArbRef}}
        set_arb_precision!(
            F.arb_jac_interpreter_cache,
            max(precision(first(U)), precision(first(x))),
        )
    else
        set_arb_precision!(F.arb_jac_interpreter_cache, precision(first(x)))
    end
    execute!(u, U, F.jac_interpreter, x, p, F.arb_jac_interpreter_cache)
end

function evaluate_and_jacobian!(
    u,
    U,
    F::InterpretedSystem,
    x::AbstractArray{<:Union{Acb,AcbRef}},
    p = nothing,
)
    if isnothing(F.acb_eval_interpreter_cache)
        F.acb_jac_interpreter_cache = InterpreterCache(Acb, F.jac_interpreter; prec = 128)
    end
    if U isa AbstractArray{<:Union{Acb,AcbRef}}
        set_arb_precision!(
            F.acb_jac_interpreter_cache,
            max(precision(first(U)), precision(first(x))),
        )
    else
        set_arb_precision!(F.acb_jac_interpreter_cache, precision(first(x)))
    end
    execute!(u, U, F.jac_interpreter, x, p, F.acb_jac_interpreter_cache)
end

function (F::InterpretedSystem)(x::AbstractArray{<:Union{Arb,ArbRef}}, p = nothing)
    u = ArbVector(first(size(F)); prec = precision(first(x)))
    evaluate!(u, F, x, p)
    u
end
function (F::InterpretedSystem)(x::AbstractArray{<:Union{Acb,AcbRef}}, p = nothing)
    u = AcbVector(first(size(F)); prec = precision(first(x)))
    evaluate!(u, F, x, p)
    u
end
(F::System)(x::AbstractVector{<:Union{Arb,ArbRef,Acb,AcbRef}}, p::Nothing = nothing) =
    InterpretedSystem(F)(x, p)
(F::System)(x::AbstractVector{<:Union{Arb,ArbRef,Acb,AcbRef}}, p::AbstractArray) =
    InterpretedSystem(F)(x, p)
(F::System)(x::AbstractMatrix{<:Union{Arb,ArbRef,Acb,AcbRef}}, p = nothing) =
    InterpretedSystem(F)(x, p)

"""
    InterpretedHomotopy <: AbstractHomotopy

An [`AbstractHomotopy`](@ref) which automatically generates a program for the
fast evaluation of `H` and its Jacobian. The program is however, not compiled
but rather interpreted. See also [`CompiledHomotopy`](@ref).

    InterpretedHomotopy(H::Homotopy; optimizations = true)

Construct an `InterpretedHomotopy` from the given [`Homotopy`](@ref) `H`.
If `optimizations = true` then [`optimize`](@ref) is called on `H` before compiling.
"""
struct InterpretedHomotopy{T₁,T₂} <: AbstractHomotopy
    homotopy::Homotopy
    eval_interpreter::Interpreter{T₁}
    eval_interpreter_cache::InterpreterCache{Vector{ComplexF64}}
    eval_interpreter_cache_ext::InterpreterCache{Vector{ComplexDF64}}
    jac_interpreter::Interpreter{T₂}
    jac_interpreter_cache::InterpreterCache{Vector{ComplexF64}}
    taylor_caches::Tuple{
        InterpreterCache{TaylorVector{2,ComplexF64}},
        InterpreterCache{TaylorVector{3,ComplexF64}},
        InterpreterCache{TaylorVector{4,ComplexF64}},
    }
end

function InterpretedHomotopy(H::Homotopy; optimizations::Bool = true)
    H = optimizations ? optimize(H) : H
    eval_interpreter = interpreter(H)
    jac_interpreter = jacobian_interpreter(H)
    N = cache_min_length(eval_interpreter)
    taylor_caches = (
        InterpreterCache(TaylorVector{2}(ComplexF64, N)),
        InterpreterCache(TaylorVector{3}(ComplexF64, N)),
        InterpreterCache(TaylorVector{4}(ComplexF64, N)),
    )

    InterpretedHomotopy(
        H,
        eval_interpreter,
        InterpreterCache(ComplexF64, eval_interpreter),
        InterpreterCache(ComplexDF64, eval_interpreter),
        jac_interpreter,
        InterpreterCache(ComplexF64, jac_interpreter),
        taylor_caches,
    )
end

Base.size(H::InterpretedHomotopy) = size(H.homotopy)
variables(H::InterpretedHomotopy) = variables(H.homotopy)
parameters(H::InterpretedHomotopy) = parameters(H.homotopy)
Homotopy(H::InterpretedHomotopy) = H.homotopy

function Base.show(io::IO, H::InterpretedHomotopy)
    print(io, "Interpreted: ")
    show(io, H.homotopy)
end

(H::InterpretedHomotopy)(x, t, p = nothing) = H.homotopy(x, t, p)
function evaluate!(
    u,
    H::InterpretedHomotopy,
    x,
    t,
    p = nothing,
    cache = H.eval_interpreter_cache,
)
    execute!(u, H.eval_interpreter, x, t, p, cache)
end
function evaluate!(u, H::InterpretedHomotopy, x::AbstractArray{ComplexDF64}, t, p = nothing)
    execute!(u, H.eval_interpreter, x, t, p, H.eval_interpreter_cache_ext)
end
function evaluate_and_jacobian!(
    u,
    U,
    H::InterpretedHomotopy,
    x,
    t,
    p = nothing,
    cache = H.jac_interpreter_cache,
)
    execute!(u, U, H.jac_interpreter, x, t, p, cache)
    nothing
end

function jacobian!(
    U,
    H::InterpretedHomotopy,
    x,
    t,
    p = nothing,
    cache = H.jac_interpreter_cache,
)
    execute!(nothing, U, H.jac_interpreter, x, t, p, cache)
    nothing
end

for M = 1:3
    @eval function taylor!(
        u::AbstractVecOrMat,
        Order::Val{$M},
        H::InterpretedHomotopy,
        x,
        t,
        p::Union{AbstractVector,Nothing} = nothing,
    )
        execute!(u, Order, H.eval_interpreter, x, t, p, H.taylor_caches[$M])
        u
    end
end


text:
using LinearAlgebra;
using Printf;

function linspace(start::Number, stop::Number, n)
	linx = zeros(Float64, n);
	delta = (stop - start) / (n - 1);
	for i = 1:n
		linx[i] = start + (i - 1) * delta;
	end
	return linx;
end

function tdma(n, a, b, c, r)
	ans = zeros(Float64, n);
	g = zeros(Float64, n);
	g[1] = b[1];
	s = zeros(Float64, n);
	s[1] = r[1];

	ans[1] = s[1] / g[1];

	for i = 2:n-1
		g[i] = b[i] - a[i] * c[i - 1] / g[i - 1];
		s[i] = r[i] - a[i] * s[i - 1] / g[i - 1];
		ans[i] = s[i] / g[i];
	end
	for i = n-2:1
		ans[i] = (s[i] - c[i] * ans[i + 1]) / g[i];
	end
	return ans
end

function init()
	global nx = 5;
	global ny = 5;
	global lx = Float64(100.0);
	global ly = Float64(100.0);
	global x = linspace(0, lx, nx);
	global y = linspace(0, ly, ny);
	global dx = lx / (nx - 1);
	global dy = ly / (ny - 1);
	global t = 0.0;
	global tlims = 300.0;;
	global dt = 10.0;
	global kappa = 0.835;
	global beta_x = kappa * dt / 2.0dx^2;
	global beta_y = kappa * dt / 2.0dy^2;
	global u = zeros(Float64, (ny, nx));
	for j in 1:ny
		u[j, 1] = 75.0;
		u[j, end] = 50.0;
	end
	for i in 1:nx
		u[1, i] = 0.0;
		u[end, i] = 100.0;
	end
	u[begin, begin] = (75.0 + 0.0) / 2.0;
	u[end, begin] = (75.0 + 100.0) / 2.0;
	u[begin, end] = (0.0 + 50.0) / 2.0;
	u[end, end] = (100.0 + 50.0) / 2.0;
end

function output(x::Array, y::Array, u::Array, t::Float64)
	filename = @sprintf("adi_cn2d_t=%03d.dat", t);
	open(filename, "w") do io
		for i in 1:ny
			for j in 1:nx
				write(io, @sprintf("%.5f %.5f %.5f\n", x[j], y[i], u[i, j]));
			end
		end
	end
end

function cn_direction_x(u::Array)
	a = zeros(Float64, nx - 2);
	b = zeros(Float64, nx - 2);
	c = zeros(Float64, nx - 2);
	r = zeros(Float64, nx - 2);
	for j = 1:ny-2
		for i = 1:nx-2
			a[i] = -beta_y;
			b[i] = 1.0 + 2.0beta_y;
			c[i] = -beta_y;
			r[i] = u[j+1, i+1] + beta_x*(u[j+2, i+1] - 2u[j+1, i+1] + u[j, i+1]);
		end
		r[begin] = r[begin] - a[begin] * u[j+1, begin];
		r[end] = r[end] - c[end] * u[j+1, end];
		u[j+1, begin+1:end-1] = tdma(nx - 2, a, b, c, r)
	end
end

function cn_direction_y(u::Array)
	a = zeros(Float64, nx - 2);
	b = zeros(Float64, nx - 2);
	c = zeros(Float64, nx - 2);
	r = zeros(Float64, nx - 2);
	for i in 1:nx-2
		for j in 1:ny-2
			a[j] = -beta_x;
			b[j] = 1.0 + 2.0beta_x;
			c[j] = -beta_x;
			r[j] = u[j+1, i+1] + beta_y * (u[j+1, i+2] - 2.0u[j+1, i+1] + u[j+1, i])
		end
		r[begin] = r[begin] - a[begin]*u[begin, i+1];
		r[end] = r[end] - c[begin]*u[end, i+1];
		u[begin+1:end-1, i] = tdma(ny-2, a, b, c, r);
	end
end

function main()
	init();
	@printf("nx = %d, ny = %d, dt = %.4f, dx = %.4f, dy = %.4f\n", nx, ny, dt, dx, dy);
	loop_count = 1;
	while t <= tlims
		@printf("t = %.5f\n", t);
		## 奇数回目はx方向を陰的にする.
		if loop_count % 2 == 1
			## x方向にcrank nicolson
			cn_direction_x(u);
			## y方向にcrank nicolson
			cn_direction_y(u);
			output(x, y, u, t);
		## 偶数回目はy方向を陰的にする.
		else
			## y方向にcrank nicolson
			cn_direction_y(u);
			## x方向にcrank nicolson
			cn_direction_x(u);
			output(x, y, u, t);
		end
		loop_count += 1;
		global t += dt;
	end
end

main()


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Types

using UUIDs
using Random
using Dates
import LibGit2
import REPL
import Base.string
using REPL.TerminalMenus

using TOML
import ..Pkg, ..Registry
import ..Pkg: GitTools, depots, depots1, logdir, set_readonly, safe_realpath, pkg_server, stdlib_dir, stdlib_path, isurl, stderr_f
import Base.BinaryPlatforms: Platform
using ..Pkg.Versions

import Base: SHA1
using SHA

export UUID, SHA1, VersionRange, VersionSpec,
    PackageSpec, PackageEntry, EnvCache, Context, GitRepo, Context!, Manifest, Project, err_rep,
    PkgError, pkgerror, has_name, has_uuid, is_stdlib, stdlib_version, is_unregistered_stdlib, stdlibs, write_env, write_env_usage, parse_toml, find_registered!,
    project_resolve!, project_deps_resolve!, manifest_resolve!, registry_resolve!, stdlib_resolve!, handle_repos_develop!, handle_repos_add!, ensure_resolved,
    registered_name,
    manifest_info,
    read_project, read_package, read_manifest,
    PackageMode, PKGMODE_MANIFEST, PKGMODE_PROJECT, PKGMODE_COMBINED,
    UpgradeLevel, UPLEVEL_FIXED, UPLEVEL_PATCH, UPLEVEL_MINOR, UPLEVEL_MAJOR,
    PreserveLevel, PRESERVE_ALL, PRESERVE_DIRECT, PRESERVE_SEMVER, PRESERVE_TIERED, PRESERVE_NONE,
    projectfile_path, manifestfile_path

# Load in data about historical stdlibs
include("HistoricalStdlibs.jl")

deepcopy_toml(x) = x
function deepcopy_toml(@nospecialize(x::Vector))
    d = similar(x)
    for (i, v) in enumerate(x)
        d[i] = deepcopy_toml(v)
    end
    return d
end
function deepcopy_toml(x::Dict{String, Any})
    d = Dict{String, Any}()
    sizehint!(d, length(x))
    for (k, v) in x
        d[k] = deepcopy_toml(v)
    end
    return d
end

# See loading.jl
const TOML_CACHE = Base.TOMLCache(TOML.Parser(), Dict{String, Dict{String, Any}}())
const TOML_LOCK = ReentrantLock()
# Some functions mutate the returning Dict so return a copy of the cached value here
parse_toml(toml_file::AbstractString) =
    Base.invokelatest(deepcopy_toml, Base.parsed_toml(toml_file, TOML_CACHE, TOML_LOCK))::Dict{String, Any}

#################
# Pkg Error #
#################
struct PkgError <: Exception
    msg::String
end
pkgerror(msg::String...) = throw(PkgError(join(msg)))
Base.showerror(io::IO, err::PkgError) = print(io, err.msg)

###############
# PackageSpec #
###############
@enum(UpgradeLevel, UPLEVEL_FIXED, UPLEVEL_PATCH, UPLEVEL_MINOR, UPLEVEL_MAJOR)
@enum(PreserveLevel, PRESERVE_ALL, PRESERVE_DIRECT, PRESERVE_SEMVER, PRESERVE_TIERED, PRESERVE_NONE)
@enum(PackageMode, PKGMODE_PROJECT, PKGMODE_MANIFEST, PKGMODE_COMBINED)

const VersionTypes = Union{VersionNumber,VersionSpec,UpgradeLevel}

Base.@kwdef mutable struct GitRepo
    source::Union{Nothing,String} = nothing
    rev::Union{Nothing,String} = nothing
    subdir::Union{String, Nothing} = nothing
end

Base.:(==)(r1::GitRepo, r2::GitRepo) =
    r1.source == r2.source && r1.rev == r2.rev && r1.subdir == r2.subdir

mutable struct PackageSpec
    name::Union{Nothing,String}
    uuid::Union{Nothing,UUID}
    version::Union{Nothing,VersionTypes,String}
    tree_hash::Union{Nothing,SHA1}
    repo::GitRepo
    path::Union{Nothing,String}
    pinned::Bool
    # used for input only
    url
    rev
    subdir

end
function PackageSpec(; name::Union{Nothing,AbstractString} = nothing,
                       uuid::Union{Nothing,UUID,AbstractString} = nothing,
                       version::Union{Nothing,VersionTypes,AbstractString} = VersionSpec(),
                       tree_hash::Union{Nothing,SHA1} = nothing,
                       repo::GitRepo = GitRepo(),
                       path::Union{Nothing,AbstractString} = nothing,
                       pinned::Bool = false,
                       url = nothing,
                       rev = nothing,
                       subdir = nothing)
    uuid = uuid === nothing ? nothing : UUID(uuid)
    return PackageSpec(name, uuid, version, tree_hash, repo, path, pinned, url, rev, subdir)
end
PackageSpec(name::AbstractString) = PackageSpec(;name=name)
PackageSpec(name::AbstractString, uuid::UUID) = PackageSpec(;name=name, uuid=uuid)
PackageSpec(name::AbstractString, version::VersionTypes) = PackageSpec(;name=name, version=version)
PackageSpec(n::AbstractString, u::UUID, v::VersionTypes) = PackageSpec(;name=n, uuid=u, version=v)

function Base.:(==)(a::PackageSpec, b::PackageSpec)
    return a.name == b.name && a.uuid == b.uuid && a.version == b.version &&
    a.tree_hash == b.tree_hash && a.repo == b.repo && a.path == b.path &&
    a.pinned == b.pinned
end

function err_rep(pkg::PackageSpec)
    x = pkg.name !== nothing && pkg.uuid !== nothing ? x = "$(pkg.name) [$(string(pkg.uuid)[1:8])]" :
        pkg.name !== nothing ? pkg.name :
        pkg.uuid !== nothing ? string(pkg.uuid)[1:8] :
        pkg.repo.source
    return "`$x`"
end

has_name(pkg::PackageSpec) = pkg.name !== nothing
has_uuid(pkg::PackageSpec) = pkg.uuid !== nothing
isresolved(pkg::PackageSpec) = pkg.uuid !== nothing && pkg.name !== nothing

function Base.show(io::IO, pkg::PackageSpec)
    vstr = repr(pkg.version)
    f = []
    pkg.name !== nothing && push!(f, "name" => pkg.name)
    pkg.uuid !== nothing && push!(f, "uuid" => pkg.uuid)
    pkg.tree_hash !== nothing && push!(f, "tree_hash" => pkg.tree_hash)
    pkg.path !== nothing && push!(f, "path" => pkg.path)
    pkg.url !== nothing && push!(f, "url" => pkg.url)
    pkg.rev !== nothing && push!(f, "rev" => pkg.rev)
    pkg.subdir !== nothing && push!(f, "subdir" => pkg.subdir)
    pkg.pinned && push!(f, "pinned" => pkg.pinned)
    push!(f, "version" => (vstr == "VersionSpec(\"*\")" ? "*" : vstr))
    if pkg.repo.source !== nothing
        push!(f, "repo/source" => string("\"", pkg.repo.source, "\""))
    end
    if pkg.repo.rev !== nothing
        push!(f, "repo/rev" => pkg.repo.rev)
    end
    if pkg.repo.subdir !== nothing
        push!(f, "repo/subdir" => pkg.repo.subdir)
    end
    print(io, "PackageSpec(\n")
    for (field, value) in f
        print(io, "  ", field, " = ", value, "\n")
    end
    print(io, ")")
end

############
# EnvCache #
############

function projectfile_path(env_path::String; strict=false)
    for name in Base.project_names
        maybe_file = joinpath(env_path, name)
        isfile(maybe_file) && return maybe_file
    end
    return strict ? nothing : joinpath(env_path, "Project.toml")
end

function manifestfile_path(env_path::String; strict=false)
    for name in Base.manifest_names
        maybe_file = joinpath(env_path, name)
        isfile(maybe_file) && return maybe_file
    end
    if strict
        return nothing
    else
        project = basename(projectfile_path(env_path)::String)
        idx = findfirst(x -> x == project, Base.project_names)
        @assert idx !== nothing
        return joinpath(env_path, Base.manifest_names[idx])
    end
end

function find_project_file(env::Union{Nothing,String}=nothing)
    project_file = nothing
    if env isa Nothing
        project_file = Base.active_project()
        project_file === nothing && pkgerror("no active project")
    elseif startswith(env, '@')
        project_file = Base.load_path_expand(env)
        project_file === nothing && pkgerror("package environment does not exist: $env")
    elseif env isa String
        if isdir(env)
            isempty(readdir(env)) || pkgerror("environment is a package directory: $env")
            project_file = joinpath(env, Base.project_names[end])
        else
            project_file = endswith(env, ".toml") ? abspath(env) :
                abspath(env, Base.project_names[end])
        end
    end
    @assert project_file isa String &&
        (isfile(project_file) || !ispath(project_file) ||
         isdir(project_file) && isempty(readdir(project_file)))
    return Pkg.safe_realpath(project_file)
end

Base.@kwdef mutable struct Compat
    val::VersionSpec
    str::String
end

Base.@kwdef mutable struct Project
    other::Dict{String,Any} = Dict{String,Any}()
    # Fields
    name::Union{String, Nothing} = nothing
    uuid::Union{UUID, Nothing} = nothing
    version::Union{VersionTypes, Nothing} = nothing
    manifest::Union{String, Nothing} = nothing
    # Sections
    deps::Dict{String,UUID} = Dict{String,UUID}()
    extras::Dict{String,UUID} = Dict{String,UUID}()
    targets::Dict{String,Vector{String}} = Dict{String,Vector{String}}()
    compat::Dict{String,Compat} = Dict{String,Compat}()
end
Base.:(==)(t1::Project, t2::Project) = all(x -> (getfield(t1, x) == getfield(t2, x))::Bool, fieldnames(Project))
Base.hash(x::Project, h::UInt) = foldr(hash, [getfield(t, x) for x in fieldnames(Project)], init=h)

Base.@kwdef mutable struct PackageEntry
    name::Union{String,Nothing} = nothing
    version::Union{VersionNumber,Nothing} = nothing
    path::Union{String,Nothing} = nothing
    pinned::Bool = false
    repo::GitRepo = GitRepo()
    tree_hash::Union{Nothing,SHA1} = nothing
    deps::Dict{String,UUID} = Dict{String,UUID}()
    uuid::Union{Nothing, UUID} = nothing
    other::Union{Dict,Nothing} = nothing
end
Base.:(==)(t1::PackageEntry, t2::PackageEntry) = t1.name == t2.name &&
    t1.version == t2.version &&
    t1.path == t2.path &&
    t1.pinned == t2.pinned &&
    t1.repo == t2.repo &&
    t1.tree_hash == t2.tree_hash &&
    t1.deps == t2.deps &&
    t1.uuid == t2.uuid
    # omits `other`
Base.hash(x::PackageEntry, h::UInt) = foldr(hash, [x.name, x.version, x.path, x.pinned, x.repo, x.tree_hash, x.deps, x.uuid], init=h)  # omits `other`

Base.@kwdef mutable struct Manifest
    julia_version::Union{Nothing,VersionNumber} = nothing # only set to VERSION when resolving
    manifest_format::VersionNumber = v"2.0.0"
    deps::Dict{UUID,PackageEntry} = Dict{UUID,PackageEntry}()
    other::Dict{String,Any} = Dict{String,Any}()
end
Base.:(==)(t1::Manifest, t2::Manifest) = all(x -> (getfield(t1, x) == getfield(t2, x))::Bool, fieldnames(Manifest))
Base.hash(m::Manifest, h::UInt) = foldr(hash, [getfield(m, x) for x in fieldnames(Manifest)], init=h)
Base.getindex(m::Manifest, i_or_key) = getindex(m.deps, i_or_key)
Base.get(m::Manifest, key, default) = get(m.deps, key, default)
Base.setindex!(m::Manifest, i_or_key, value) = setindex!(m.deps, i_or_key, value)
Base.iterate(m::Manifest) = iterate(m.deps)
Base.iterate(m::Manifest, i::Int) = iterate(m.deps, i)
Base.length(m::Manifest) = length(m.deps)
Base.empty!(m::Manifest) = empty!(m.deps)
Base.values(m::Manifest) = values(m.deps)
Base.keys(m::Manifest) = keys(m.deps)
Base.haskey(m::Manifest, key) = haskey(m.deps, key)

function Base.show(io::IO, pkg::PackageEntry)
    f = []
    pkg.name        !== nothing && push!(f, "name"      => pkg.name)
    pkg.version     !== nothing && push!(f, "version"   => pkg.version)
    pkg.tree_hash   !== nothing && push!(f, "tree_hash" => pkg.tree_hash)
    pkg.path        !== nothing && push!(f, "dev/path"  => pkg.path)
    pkg.pinned                  && push!(f, "pinned"    => pkg.pinned)
    pkg.repo.source !== nothing && push!(f, "url/path"  => "`$(pkg.repo.source)`")
    pkg.repo.rev    !== nothing && push!(f, "rev"       => pkg.repo.rev)
    pkg.repo.subdir !== nothing && push!(f, "subdir"    => pkg.repo.subdir)
    print(io, "PackageEntry(\n")
    for (field, value) in f
        print(io, "  ", field, " = ", value, "\n")
    end
    print(io, ")")
end

mutable struct EnvCache
    # environment info:
    env::Union{Nothing,String}
    # paths for files:
    project_file::String
    manifest_file::String
    # name / uuid of the project
    pkg::Union{PackageSpec, Nothing}
    # cache of metadata:
    project::Project
    manifest::Manifest
    # What these where at creation of the EnvCache
    original_project::Project
    original_manifest::Manifest
end

function EnvCache(env::Union{Nothing,String}=nothing)
    project_file = find_project_file(env)
    project_dir = dirname(project_file)
    # read project file
    project = read_project(project_file)
    # initialize project package
    if project.name !== nothing && project.uuid !== nothing
        project_package = PackageSpec(
            name = project.name,
            uuid = project.uuid,
            version = something(project.version, VersionNumber("0.0")),
            path = project_dir,
        )
    else
        project_package = nothing
    end
    # determine manifest file
    dir = abspath(project_dir)
    manifest_file = project.manifest
    manifest_file = manifest_file !== nothing ?
        abspath(manifest_file) : manifestfile_path(dir)::String
    write_env_usage(manifest_file, "manifest_usage.toml")
    manifest = read_manifest(manifest_file)

    env′ = EnvCache(env,
        project_file,
        manifest_file,
        project_package,
        project,
        manifest,
        deepcopy(project),
        deepcopy(manifest),
        )

    return env′
end

include("project.jl")
include("manifest.jl")

# ENV variables to set some of these defaults?
Base.@kwdef mutable struct Context
    env::EnvCache = EnvCache()
    io::IO = stderr_f()
    use_git_for_all_downloads::Bool = false
    use_only_tarballs_for_downloads::Bool = false
    num_concurrent_downloads::Int = 8

    # Registris
    registries::Vector{Registry.RegistryInstance} = Registry.reachable_registries()

    # The Julia Version to resolve with respect to
    julia_version::Union{VersionNumber,Nothing} = VERSION
end

project_uuid(env::EnvCache) = env.pkg === nothing ? nothing : env.pkg.uuid
collides_with_project(env::EnvCache, pkg::PackageSpec) =
    is_project_name(env, pkg.name) || is_project_uuid(env, pkg.uuid)
is_project(env::EnvCache, pkg::PackageSpec) = is_project_uuid(env, pkg.uuid)
is_project_name(env::EnvCache, name::String) =
    env.pkg !== nothing && env.pkg.name == name
is_project_name(env::EnvCache, name::Nothing) = false
is_project_uuid(env::EnvCache, uuid::UUID) = project_uuid(env) == uuid

###########
# Context #
###########

const STDLIB = Ref{Dict{UUID,String}}()
function load_stdlib()
    stdlib = Dict{UUID,String}()
    for name in readdir(stdlib_dir())
        projfile = projectfile_path(stdlib_path(name); strict=true)
        nothing === projfile && continue
        project = parse_toml(projfile)
        uuid = get(project, "uuid", nothing)
        nothing === uuid && continue
        stdlib[UUID(uuid)] = name
    end
    return stdlib
end

function stdlibs()
    if !isassigned(STDLIB)
        STDLIB[] = load_stdlib()
    end
    return STDLIB[]
end
is_stdlib(uuid::UUID) = uuid in keys(stdlibs())

# Find the entry in `STDLIBS_BY_VERSION`
# that corresponds to the requested version, and use that.
# If we can't find one, defaults to `UNREGISTERED_STDLIBS`
function get_last_stdlibs(julia_version::VersionNumber)
    last_stdlibs = UNREGISTERED_STDLIBS
    for (version, stdlibs) in STDLIBS_BY_VERSION
        if VersionNumber(julia_version.major, julia_version.minor, julia_version.patch) < version
            break
        end
        last_stdlibs = stdlibs
    end
    return last_stdlibs
end
# If `julia_version` is set to `nothing`, that means (essentially) treat all registered
# stdlibs as normal packages so that we get the latest versions of everything, ignoring
# julia compat.  So we set the list of stdlibs to that of only the unregistered stdlibs.
get_last_stdlibs(::Nothing) = UNREGISTERED_STDLIBS

# Allow asking if something is an stdlib for a particular version of Julia
function is_stdlib(uuid::UUID, julia_version::Union{VersionNumber, Nothing})
    # Only use the cache if we are asking for stdlibs in a custom Julia version
    if julia_version == VERSION
        return is_stdlib(uuid)
    end

    last_stdlibs = get_last_stdlibs(julia_version)
    # Note that if the user asks for something like `julia_version = 0.7.0`, we'll
    # fall through with an empty `last_stdlibs`, which will always return `false`.
    return uuid in keys(last_stdlibs)
end

# Return the version of a stdlib with respect to a particular Julia version, or
# `nothing` if that stdlib is not versioned.  We only store version numbers for
# stdlibs that are external and thus could be installed from their repositories,
# e.g. things like `GMP_jll`, `Tar`, etc...
function stdlib_version(uuid::UUID, julia_version::Union{VersionNumber,Nothing})
    last_stdlibs = get_last_stdlibs(julia_version)
    if !(uuid in keys(last_stdlibs))
        return nothing
    end
    return last_stdlibs[uuid][2]
end

is_unregistered_stdlib(uuid::UUID) = haskey(UNREGISTERED_STDLIBS, uuid)

Context!(kw_context::Vector{Pair{Symbol,Any}})::Context =
    Context!(Context(); kw_context...)
function Context!(ctx::Context; kwargs...)
    for (k, v) in kwargs
        setfield!(ctx, k, v)
    end
    return ctx
end

function write_env_usage(source_file::AbstractString, usage_filepath::AbstractString)
    # Don't record ghost usage
    !isfile(source_file) && return

    # Ensure that log dir exists
    !ispath(logdir()) && mkpath(logdir())

    usage_file = joinpath(logdir(), usage_filepath)
    timestamp = now()

    ## Atomically write usage file
    while true
        # read existing usage file
        usage = if isfile(usage_file)
            TOML.parsefile(usage_file)
        else
            Dict{String, Any}()
        end

        # record new usage
        usage[source_file] = [Dict("time" => timestamp)]

        # keep only latest usage info
        for k in keys(usage)
            times = map(d -> Dates.DateTime(d["time"]), usage[k])
            usage[k] = [Dict("time" => maximum(times))]
        end

        # Write to a temp file in the same directory as the destination
        temp_usage_file = tempname(logdir())
        open(temp_usage_file, "w") do io
            TOML.print(io, usage, sorted=true)
        end

        # Move the temp file into place, replacing the original
        mv(temp_usage_file, usage_file, force = true)

        # Check that the new file has what we want in it
        new_usage = if isfile(usage_file)
            TOML.parsefile(usage_file)
        else
            Dict{String, Any}()
        end
        if haskey(new_usage, source_file)
            for e in new_usage[source_file]
                if Dates.DateTime(e["time"]) >= timestamp
                    return
                end
            end
        end
        # If not, try again
    end
end

function read_package(path::String)
    project = read_project(path)
    if project.name === nothing
        pkgerror("expected a `name` entry in project file at `$(abspath(path))`")
    end
    if project.uuid === nothing
        pkgerror("expected a `uuid` entry in project file at `$(abspath(path))`")
    end
    name = project.name
    if !isfile(joinpath(dirname(path), "src", "$name.jl"))
        pkgerror("expected the file `src/$name.jl` to exist for package `$name` at `$(dirname(path))`")
    end
    return project
end

const refspecs = ["+refs/*:refs/remotes/cache/*"]

function relative_project_path(project_file::String, path::String)
    # compute path relative the project
    # realpath needed to expand symlinks before taking the relative path
    return relpath(Pkg.safe_realpath(abspath(path)),
                   Pkg.safe_realpath(dirname(project_file)))
end

function devpath(env::EnvCache, name::AbstractString, shared::Bool)
    @assert name != ""
    dev_dir = shared ? abspath(Pkg.devdir()) : joinpath(dirname(env.project_file), "dev")
    return joinpath(dev_dir, name)
end

function handle_repo_develop!(ctx::Context, pkg::PackageSpec, shared::Bool)
    # First, check if we can compute the path easily (which requires a given local path or name)
    is_local_path = pkg.repo.source !== nothing && !isurl(pkg.repo.source)
    if is_local_path || pkg.name !== nothing
        dev_path = is_local_path ? pkg.repo.source : devpath(ctx.env, pkg.name, shared)
        if pkg.repo.subdir !== nothing
            dev_path = joinpath(dev_path, pkg.repo.subdir)
        end
        # If given an explicit local path, that needs to exist
        if is_local_path && !isdir(dev_path)
            if isfile(dev_path)
                pkgerror("Dev path `$(dev_path)` is a file, but a directory is required.")
            else
                pkgerror("Dev path `$(dev_path)` does not exist.")
            end
        end
        if isdir(dev_path)
            resolve_projectfile!(ctx.env, pkg, dev_path)
            if is_local_path
                pkg.path = isabspath(dev_path) ? dev_path : relative_project_path(ctx.env.project_file, dev_path)
            else
                pkg.path = shared ? dev_path : relative_project_path(ctx.env.project_file, dev_path)
            end
            return false
        end
    end
    # If we dev by name and it is in the Project + tracking a repo in the source we can get the repo from the Manifest
    if pkg.name !== nothing && pkg.uuid === nothing
        uuid = get(ctx.env.project.deps, pkg.name, nothing)
        if uuid !== nothing
            entry = manifest_info(ctx.env.manifest, uuid)
            if entry !== nothing
                pkg.repo.source = entry.repo.source
            end
        end
    end

    # Still haven't found the source, try get it from the registry
    if pkg.repo.source === nothing
        set_repo_source_from_registry!(ctx, pkg)
    end
    @assert pkg.repo.source !== nothing

    repo_path = tempname()
    cloned = false
    package_path = pkg.repo.subdir === nothing ? repo_path : joinpath(repo_path, pkg.repo.subdir)
    if !has_name(pkg)
        LibGit2.close(GitTools.ensure_clone(ctx.io, repo_path, pkg.repo.source))
        cloned = true
        resolve_projectfile!(ctx.env, pkg, package_path)
    end
    if pkg.repo.subdir !== nothing
        repo_name = split(pkg.repo.source, '/', keepempty=false)[end]
        # Make the develop path prettier.
        if endswith(repo_name, ".git")
            repo_name = chop(repo_name, tail=4)
        end
        if endswith(repo_name, ".jl")
            repo_name = chop(repo_name, tail=3)
        end
        dev_path = devpath(ctx.env, repo_name, shared)
    else
        dev_path = devpath(ctx.env, pkg.name, shared)
    end
    if isdir(dev_path)
        println(ctx.io, "Path `$(dev_path)` exists and looks like the correct repo. Using existing path.")
        new = false
    else
        mkpath(dirname(dev_path))
        if !cloned
            LibGit2.close(GitTools.ensure_clone(ctx.io, dev_path, pkg.repo.source))
        else
            mv(repo_path, dev_path)
        end
        new = true
    end
    if !has_uuid(pkg)
        resolve_projectfile!(ctx.env, pkg, dev_path)
    end
    pkg.path = shared ? dev_path : relative_project_path(ctx.env.project_file, dev_path)
    if pkg.repo.subdir !== nothing
        pkg.path = joinpath(pkg.path, pkg.repo.subdir)
    end

    return new
end

function handle_repos_develop!(ctx::Context, pkgs::AbstractVector{PackageSpec}, shared::Bool)
    new_uuids = Set{UUID}()
    for pkg in pkgs
        new = handle_repo_develop!(ctx, pkg, shared)
        new && push!(new_uuids, pkg.uuid)
        @assert pkg.path !== nothing
        @assert has_uuid(pkg)
        pkg.repo = GitRepo() # clear repo field, no longer needed
    end
    return new_uuids
end

add_repo_cache_path(url::String) = joinpath(depots1(), "clones", string(hash(url)))

function set_repo_source_from_registry!(ctx, pkg)
    registry_resolve!(ctx.registries, pkg)
    # Didn't find the package in the registry, but maybe it exists in the updated registry
    if !isresolved(pkg)
        Pkg.Operations.update_registries(ctx; force=false)
        registry_resolve!(ctx.registries, pkg)
    end
    ensure_resolved(ctx.env.manifest, [pkg]; registry=true)
    # We might have been given a name / uuid combo that does not have an entry in the registry
    for reg in ctx.registries
        regpkg = get(reg, pkg.uuid, nothing)
        regpkg === nothing && continue
        info = Pkg.Registry.registry_info(regpkg)
        url = info.repo
        url === nothing && continue
        pkg.repo.source = url
        if info.subdir !== nothing
            pkg.repo.subdir = info.subdir
        end
        return
    end
    pkgerror("Repository for package with UUID `$(pkg.uuid)` could not be found in a registry.")
end

function handle_repo_add!(ctx::Context, pkg::PackageSpec)
    # The first goal is to populate pkg.repo.source if that wasn't given explicitly
    if pkg.repo.source === nothing
        @assert pkg.repo.rev !== nothing
        # First, we try resolving against the manifest and current registry to avoid updating registries if at all possible.
        # This also handles the case where we _only_ wish to switch the tracking branch for a package.
        manifest_resolve!(ctx.env.manifest, [pkg]; force=true)
        if isresolved(pkg)
            entry = manifest_info(ctx.env.manifest, pkg.uuid)
            if entry !== nothing && entry.repo.source !== nothing # reuse source in manifest
                pkg.repo.source = entry.repo.source
            end
        end
        if pkg.repo.source === nothing
            set_repo_source_from_registry!(ctx, pkg)
        end
    end
    @assert pkg.repo.source !== nothing

    # We now have the source of the package repo, check if it is a local path and if that exists
    repo_source = pkg.repo.source
    if !isurl(pkg.repo.source)
        if isdir(pkg.repo.source)
            if !isdir(joinpath(pkg.repo.source, ".git"))
                msg = "Did not find a git repository at `$(pkg.repo.source)`"
                if isfile(joinpath(pkg.repo.source, "Project.toml")) || isfile(joinpath(pkg.repo.source, "JuliaProject.toml"))
                    msg *= ", perhaps you meant `Pkg.develop`?"
                end
                pkgerror(msg)
            end
            LibGit2.with(GitTools.check_valid_HEAD, LibGit2.GitRepo(pkg.repo.source)) # check for valid git HEAD
            pkg.repo.source = isabspath(pkg.repo.source) ? safe_realpath(pkg.repo.source) : relative_project_path(ctx.env.project_file, pkg.repo.source)
            repo_source = normpath(joinpath(dirname(ctx.env.project_file), pkg.repo.source))
        else
            pkgerror("Path `$(pkg.repo.source)` does not exist.")
        end
    end

    let repo_source = repo_source
        # The type-assertions below are necessary presumably due to julia#36454
        LibGit2.with(GitTools.ensure_clone(ctx.io, add_repo_cache_path(repo_source::Union{Nothing,String}), repo_source::Union{Nothing,String}; isbare=true)) do repo
            repo_source_typed = repo_source::Union{Nothing,String}
            GitTools.check_valid_HEAD(repo)

            # If the user didn't specify rev, assume they want the default (master) branch if on a branch, otherwise the current commit
            if pkg.repo.rev === nothing
                pkg.repo.rev = LibGit2.isattached(repo) ? LibGit2.branch(repo) : string(LibGit2.GitHash(LibGit2.head(repo)))
            end

            obj_branch = get_object_or_branch(repo, pkg.repo.rev)
            fetched = false
            if obj_branch === nothing
                fetched = true
                GitTools.fetch(ctx.io, repo, repo_source_typed; refspecs=refspecs)
                obj_branch = get_object_or_branch(repo, pkg.repo.rev)
                if obj_branch === nothing
                    pkgerror("Did not find rev $(pkg.repo.rev) in repository")
                end
            end
            gitobject, isbranch = obj_branch

            # If we are tracking a branch and are not pinned we want to update the repo if we haven't done that yet
            innerentry = manifest_info(ctx.env.manifest, pkg.uuid)
            ispinned = innerentry !== nothing && innerentry.pinned
            if isbranch && !fetched && !ispinned
                GitTools.fetch(ctx.io, repo, repo_source_typed; refspecs=refspecs)
                gitobject, isbranch = get_object_or_branch(repo, pkg.repo.rev)
            end

            # Now we have the gitobject for our ref, time to find the tree hash for it
            tree_hash_object = LibGit2.peel(LibGit2.GitTree, gitobject)
            if pkg.repo.subdir !== nothing
                try
                    tree_hash_object = tree_hash_object[pkg.repo.subdir]
                catch e
                    e isa KeyError || rethrow()
                    pkgerror("Did not find subdirectory `$(pkg.repo.subdir)`")
                end
            end
            pkg.tree_hash = SHA1(string(LibGit2.GitHash(tree_hash_object)))

            # If we already resolved a uuid, we can bail early if this package is already installed at the current tree_hash
            if has_uuid(pkg)
                version_path = Pkg.Operations.source_path(ctx.env.project_file, pkg, ctx.julia_version)
                isdir(version_path) && return false
            end

            temp_path = mktempdir()
            GitTools.checkout_tree_to_path(repo, tree_hash_object, temp_path)
            package = resolve_projectfile!(ctx.env, pkg, temp_path)

            # Now that we are fully resolved (name, UUID, tree_hash, repo.source, repo.rev), we can finally
            # check to see if the package exists at its canonical path.
            version_path = Pkg.Operations.source_path(ctx.env.project_file, pkg, ctx.julia_version)
            isdir(version_path) && return false

            # Otherwise, move the temporary path into its correct place and set read only
            mkpath(version_path)
            mv(temp_path, version_path; force=true)
            set_readonly(version_path)
            return true
        end
    end
end

function handle_repos_add!(ctx::Context, pkgs::AbstractVector{PackageSpec})
    new_uuids = Set{UUID}()
    for pkg in pkgs
        handle_repo_add!(ctx, pkg) && push!(new_uuids, pkg.uuid)
        @assert pkg.name !== nothing && pkg.uuid !== nothing && pkg.tree_hash !== nothing
    end
    return new_uuids
end

function resolve_projectfile!(env::EnvCache, pkg, project_path)
    project_file = projectfile_path(project_path; strict=true)
    project_file === nothing && pkgerror(string("could not find project file in package at `",
                    something(pkg.repo.source, pkg.path, project_path), "` maybe `subdir` needs to be specified"))
    project_data = read_package(project_file)
    if pkg.uuid === nothing || pkg.uuid == project_data.uuid
        pkg.uuid = project_data.uuid
    else
        pkgerror("UUID `$(project_data.uuid)` given by project file `$project_file` does not match given UUID `$(pkg.uuid)`")
    end
    if pkg.name === nothing || pkg.name == project_data.name
        pkg.name = project_data.name
    else
        pkgerror("name `$(project_data.name)` given by project file `$project_file` does not match given name `$(pkg.name)`")
    end
end

get_object_or_branch(repo, rev::SHA1) =
    get_object_or_branch(repo, string(rev))

# Returns nothing if rev could not be found in repo
function get_object_or_branch(repo, rev)
    try
        gitobject = LibGit2.GitObject(repo, "remotes/cache/heads/" * rev)
        return gitobject, true
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow()
    end
    try
        gitobject = LibGit2.GitObject(repo, "remotes/origin/" * rev)
        return gitobject, true
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow()
    end
    try
        gitobject = LibGit2.GitObject(repo, rev)
        return gitobject, false
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow()
    end
    return nothing
end

########################################
# Resolving packages from name or uuid #
########################################

function project_resolve!(env::EnvCache, pkgs::AbstractVector{PackageSpec})
    for pkg in pkgs
        if has_uuid(pkg) && !has_name(pkg) && Types.is_project_uuid(env, pkg.uuid)
            pkg.name = env.pkg.name
        end
        if has_name(pkg) && !has_uuid(pkg) && Types.is_project_name(env, pkg.name)
            pkg.uuid = env.pkg.uuid
        end
    end
end

# Disambiguate name/uuid package specifications using project info.
function project_deps_resolve!(env::EnvCache, pkgs::AbstractVector{PackageSpec})
    uuids = env.project.deps
    names = Dict(uuid => name for (name, uuid) in uuids)
    for pkg in pkgs
        if has_name(pkg) && !has_uuid(pkg) && pkg.name in keys(uuids)
            pkg.uuid = uuids[pkg.name]
        end
        if has_uuid(pkg) && !has_name(pkg) && pkg.uuid in keys(names)
            pkg.name = names[pkg.uuid]
        end
    end
end

# Disambiguate name/uuid package specifications using manifest info.
function manifest_resolve!(manifest::Manifest, pkgs::AbstractVector{PackageSpec}; force=false)
    uuids = Dict{String,Vector{UUID}}()
    names = Dict{UUID,String}()
    for (uuid, entry) in manifest
        push!(get!(uuids, entry.name, UUID[]), uuid)
        names[uuid] = entry.name # can be duplicate but doesn't matter
    end
    for pkg in pkgs
        if has_name(pkg) && !has_uuid(pkg) && pkg.name in keys(uuids)
            length(uuids[pkg.name]) == 1 && (pkg.uuid = uuids[pkg.name][1])
        end
        if has_uuid(pkg) && !has_name(pkg) && pkg.uuid in keys(names)
            pkg.name = names[pkg.uuid]
        end
    end
end

# Disambiguate name/uuid package specifications using registry info.
registry_resolve!(registries::Vector{Registry.RegistryInstance}, pkg::PackageSpec) = registry_resolve!(registries, [pkg])
function registry_resolve!(registries::Vector{Registry.RegistryInstance}, pkgs::AbstractVector{PackageSpec})
    # if there are no half-specified packages, return early
    any(pkg -> has_name(pkg) ⊻ has_uuid(pkg), pkgs) || return
    # collect all names and uuids since we're looking anyway
    names = [pkg.name::String for pkg in pkgs if has_name(pkg)]
    uuids = [pkg.uuid::UUID for pkg in pkgs if has_uuid(pkg)]
    for pkg in pkgs
        @assert has_name(pkg) || has_uuid(pkg)
        if has_name(pkg) && !has_uuid(pkg)
            pkg.uuid = registered_uuid(registries, pkg.name)
        end
        if has_uuid(pkg) && !has_name(pkg)
            pkg.name = registered_name(registries, pkg.uuid)
        end
    end
    return pkgs
end

function stdlib_resolve!(pkgs::AbstractVector{PackageSpec})
    for pkg in pkgs
        @assert has_name(pkg) || has_uuid(pkg)
        if has_name(pkg) && !has_uuid(pkg)
            for (uuid, name) in stdlibs()
                name == pkg.name && (pkg.uuid = uuid)
            end
        end
        if !has_name(pkg) && has_uuid(pkg)
            name = get(stdlibs(), pkg.uuid, nothing)
            nothing !== name && (pkg.name = name)
        end
    end
end

# Ensure that all packages are fully resolved
function ensure_resolved(manifest::Manifest,
        pkgs::AbstractVector{PackageSpec};
        registry::Bool=false,)::Nothing
    unresolved_uuids = Dict{String,Vector{UUID}}()
    for pkg in pkgs
        has_uuid(pkg) && continue
        uuids = [uuid for (uuid, entry) in manifest if entry.name == pkg.name]
        sort!(uuids, by=uuid -> uuid.value)
        unresolved_uuids[pkg.name] = uuids
    end
    unresolved_names = UUID[]
    for pkg in pkgs
        has_name(pkg) && continue
        push!(unresolved_names, pkg.uuid)
    end
    isempty(unresolved_uuids) && isempty(unresolved_names) && return
    msg = sprint() do io
        if !isempty(unresolved_uuids)
            println(io, "The following package names could not be resolved:")
            for (name, uuids) in sort!(collect(unresolved_uuids), by=lowercase ∘ first)
                print(io, " * $name (")
                if length(uuids) == 0
                    what = ["project", "manifest"]
                    registry && push!(what, "registry")
                    print(io, "not found in ")
                    join(io, what, ", ", " or ")
                else
                    join(io, uuids, ", ", " or ")
                    print(io, " in manifest but not in project")
                end
                println(io, ")")
            end
        end
        if !isempty(unresolved_names)
            println(io, "The following package uuids could not be resolved:")
            for uuid in unresolved_names
                println(io, " * $uuid")
            end
        end
    end
    pkgerror(msg)
end

function registered_uuids(registries::Vector{Registry.RegistryInstance}, name::String)
    uuids = Set{UUID}()
    for reg in registries
        union!(uuids, Registry.uuids_from_name(reg, name))
    end
    return uuids
end
# Determine a single UUID for a given name, prompting if needed
function registered_uuid(registries::Vector{Registry.RegistryInstance}, name::String)::Union{Nothing,UUID}
    uuids = registered_uuids(registries, name)
    length(uuids) == 0 && return nothing
    length(uuids) == 1 && return first(uuids)
    repo_infos = Tuple{String, String, UUID}[]
    for uuid in uuids
        for reg in registries
            pkg = get(reg, uuid, nothing)
            pkg === nothing && continue
            info = Pkg.Registry.registry_info(pkg)
            info.repo === nothing && continue
            push!(repo_infos, (reg.name, info.repo, uuid))
        end
    end
    unique!(repo_infos)
    if isinteractive()
        # prompt for which UUID was intended:
        menu = RadioMenu(String["Registry: $(value[1]) - Repo: $(value[2]) - UUID: $(value[3])" for value in repo_infos])
        choice = request("There are multiple registered `$name` packages, choose one:", menu)
        choice == -1 && return nothing
        return repo_infos[choice][3]
    else
        pkgerror("there are multiple registered `$name` packages, explicitly set the uuid")
    end
end

# Determine current name for a given package UUID

function registered_name(registries::Vector{Registry.RegistryInstance}, uuid::UUID)::Union{Nothing,String}
    name = nothing
    for reg in registries
        regpkg = get(reg, uuid, nothing)
        regpkg === nothing && continue
        name′ = regpkg.name
        if name !== nothing
            name′ == name || pkgerror("package `$uuid` has multiple registered name values: $name, $name′")
        end
        name = name′
    end
    return name
end

# Find package by UUID in the manifest file
manifest_info(::Manifest, uuid::Nothing) = nothing
function manifest_info(manifest::Manifest, uuid::UUID)::Union{PackageEntry,Nothing}
    return get(manifest, uuid, nothing)
end
function write_env(env::EnvCache; update_undo=true,
                   skip_writing_project::Bool=false)
    if (env.project != env.original_project) && (!skip_writing_project)
        write_project(env)
    end
    if env.manifest != env.original_manifest
        write_manifest(env)
    end
    update_undo && Pkg.API.add_snapshot_to_undo(env)
end

end # module


text:
################
# NShaskey
################
struct NShaskey{T <: AbstNS} <: Function ___NS_ns::T end

Base.getproperty(x::NShaskey, atr::Symbol) =
    hasfield(typeof(x), atr) ? Base.getfield(x, atr) : (atr ∈ x.___NS_ns._keys)

(x::NShaskey)(atr::Symbol) = atr ∈ x.___NS_ns._keys
(x::NShaskey)(atr::Symbol...) =
    (keys = x.___NS_ns._keys; collect(a ∈ keys for a ∈ atr))

################
# NSdel
################

struct NSdel{T <: AbstNS} <: Function ___NS_ns::T end

Base.getproperty(x::NSdel, atr::Symbol) =
    begin
        hasfield(typeof(x), atr) && (return Base.getfield(x, atr))

        x.___NS_ns._fixed && error("this NS is fixed!")
        x.___NS_ns._lcked && error("this NS is locked!")

        delete!(x.___NS_ns.__dict, atr)
        x.___NS_ns
    end

(x::NSdel)(atr::Symbol...) =
    begin
        if length(atr) == 0
            x.___NS_ns._fixed && error("this NS is fixed!")
            x.___NS_ns._lcked && error("this NS is locked!")
            empty!(x.___NS_ns.__dict)
        else
            for a ∈ atr; Base.getproperty(x, a); end
        end
        x.___NS_ns
    end

################
# NScstize, NSdecstize
################

struct NScstize{T <: AbstNS} <: Function ___NS_ns::T end
struct NSdecstize{T <: AbstNS} <: Function ___NS_ns::T end

Base.getproperty(x::NScstize, atr::Symbol) =
    begin
        hasfield(typeof(x), atr) && (return Base.getfield(x, atr))

        x.___NS_ns._fixed && error("this NS is fixed!")

        d = x.___NS_ns.__dict

        haskey(d, atr) ||
            error("""this NS does not have a property named '$(atr)'.""")

        isa(d[atr], NSnoncst_item) && (d[atr] = NScst_item(d[atr].obj))

        x.___NS_ns
    end

Base.getproperty(x::NSdecstize, atr::Symbol) =
    begin
        hasfield(typeof(x), atr) && (return Base.getfield(x, atr))

        x.___NS_ns._fixed && error("this NS is fixed!")

        d = x.___NS_ns.__dict

        haskey(d, atr) ||
            error("""this NS does not have a property named '$(atr)'.""")

        isa(x.___NS_ns.__dict[atr], NScst_item) &&
            (x.___NS_ns.__dict[atr] = NSnoncst_item(x.___NS_ns.__dict[atr].obj))

        x.___NS_ns
    end

(x::Union{NScstize, NSdecstize})(atr::Symbol...) =
    begin
        length(atr) == 0 && (return x(x.___NS_ns._keys...))

        for a ∈ atr; Base.getproperty(x, a); end
        x.___NS_ns
    end

################
# NSTag
################

struct NSTag{T, CST} ___NS_ns::AbstNS end

Base.getproperty(x::NSTag{TAG}, atr::Symbol) where TAG =
    begin
        Base.hasfield(typeof(x), atr) && (return Base.getfield(x, atr))

        ns = x.___NS_ns

        !ns.haskey(atr) &&
            error("""this NS does not have a property named '$(atr)'.""")

        o = ns.__dict[atr].obj
        to = typeof(o)
        to <: NSTagFunc && to.parameters[1] == TAG && (return o.fnc)

        error("""'$(atr)' is not a $(string(typeof(x).parameters[1])).""")
    end

Base.setproperty!(x::NSTag{TAG, TF}, atr::Symbol, f) where TAG where TF =
    begin
        Base.hasfield(typeof(x), atr) && (Base.setproperty!(x, atr, f); return)

        Base.setproperty!(x.___NS_ns, atr,
                          (TF ? NScst_item : NSnoncst_item)(NSTagFunc{TAG}(f)))
        #=
        Base.setproperty!(x.___NS_ns, atr,
                          typeof(x.___NS_ns) == __NSX_CodeMode
                          ? (TF ? NScst_item : NSnoncst_item)(NSTagFunc{TAG}(f))
                          : _MakeItem(x, f))
        =#
    end

Base.setproperty!(x::NSTag{:prp, TF}, atr::Symbol, f) where TF =
    begin
        Base.hasfield(typeof(x), atr) && (Base.setproperty!(x, atr, f); return)

        typeof(x.___NS_ns) != __NSX_CodeMode && x.___NS_ns.haskey(atr) &&
            x.___NS_ns.del(atr)
        Base.setproperty!(x.___NS_ns, atr,
                          (TF ? NScst_item : NSnoncst_item)(NSTagFunc{:prp}(f)))

        #=
        if typeof(x.___NS_ns) == __NSX_CodeMode
            Base.setproperty!(x.___NS_ns, atr,
                              (TF ? NScst_item : NSnoncst_item)(NSTagFunc{:prp}(f)))
        else
            x.___NS_ns.haskey(atr) && x.___NS_ns.del(atr)
            Base.setproperty!(x.___NS_ns, atr, _MakeItem(x, f))
        end
        =#
    end

################
# NScst
################

struct NScst{T <: AbstNS} ___NS_ns::T end

Base.getproperty(x::NScst, atr::Symbol) =
    (Base.hasfield(typeof(x), atr)
     ? Base.getfield(x, atr)
     : NSTag{atr, true}(x.___NS_ns))

Base.setproperty!(x::NScst, atr::Symbol, o) =
        (Base.hasfield(typeof(x), atr)
         ? Base.setproperty!(x, atr, f)
         : Base.setproperty!(x.___NS_ns, atr, NScst_item(o)))
#=
_MakeItem(x::NSTag{T, false}, f) where T = NSnoncst_item(NSTagFunc{T}(f))
_MakeItem(x::NSTag{T, true},  f) where T = NScst_item(NSTagFunc{T}(f))
_MakeItem(x::NSTag{:dfn, false}, f) = NSnoncst_item(f(x.___NS_ns))
_MakeItem(x::NSTag{:dfn, true},  f) = NScst_item(f(x.___NS_ns))
=#

###############################
# NSTagFunc
###############################

struct NSTagFunc{T} <: Function fnc end

NSDfn = NSTagFunc{:dfn}
NSReq = NSTagFunc{:req}
NSMth = NSTagFunc{:mth}
NSFnc = NSTagFunc{:fnc}
NSPrp = NSTagFunc{:prp}

(f::NSTagFunc)(a...; ka...) = f.fnc(a...; ka...)
(mth::NSMth)(self) = (a...; ka...)->mth.fnc(self, a...; ka...)
# (dfn::NSDfn)(self) = dfn.fnc(self)
# (req::NSReq)(self) = req.fnc(self)
# (mth::NSMth)(a...; ka...) = mth.fnc(a...; ka...)
# (fnc::NSFnc)(a...; ka...) = fnc.fnc(a...; ka...)
# (prp::NSPrp)(a...) = prp.fnc(a...)

macro nsdfn(ex) esc(Expr(:call, :(Wild.NSDfn), ex)) end
macro nsprp(ex) esc(Expr(:call, :(Wild.NSPrp), ex)) end
macro nsreq(ex) esc(Expr(:call, :(Wild.NSReq), ex)) end
macro nsmth(ex) esc(Expr(:call, :(Wild.NSMth), ex)) end
macro nsfnc(ex) esc(Expr(:call, :(Wild.NSFnc), ex)) end
macro nscst(ex) esc(Expr(:call, :(Wild.NScst_item), ex)) end


text:
export rmkey, update, mergef

"""
Remove a `key` from named tuple `nt` has been removed

```jldoctest
x = (a = 3, b = 4, c =12)
rmkey(tag, Val{:c})
(a = 3, b = 4)
```
"""
@generated function rmkey(nt::NamedTuple{K, V}, key::Type{Val{T}}) where {K, V, T}
  T isa Symbol || throw(ArgumentError("Usage: `rmkey(x, Val{:X})`"))
  # args = [:($k = nt.$k) for k in K if k != T]
  ks = tuple((k for k in K if k != T)...)
  vs = [:(nt.$k) for k in K if k != T]
  # FIXME: This can be made more efficient, using TAGS
  # Expr(:tuple, args...)
  :(NamedTuple{$ks}(($(vs...),)))
end

@pre rmkey(nt, key) = key in keys(nt) "Key `key` not in `nt`"
@post rmkey(nt, key) = keys(__ret__) == setdiff(keys(nt), key) "Keys of returned nt same as nt minus key"
@post rmkey(nt, key) = all([__ret__[k] == nt[k] for k in __ret__]) "Values of returned nt same as nt minus key"

"""
Update a named tuple

```
update((x = 3, y = 2, z = 1), Val{:x}, 7)
(x = 7, y = 2, z = 1)
```

"""
@generated function update(nt::NamedTuple{K, V}, key::Type{Val{T}}, val) where {K, V, T}
  args = []
  for k in K
    if k == T
      push!(args, :($k = val))
    else
      push!(args, :($k = nt.$k))
    end
  end
  # FIXME, return named tuple
  Expr(:tuple, args...)
end
@post update(nt, key) = keys(__ret__) == keys(nt) "Keys of returned nt same as nt"
@post update(nt, key) = all([__ret__[k] == nt[k] for k in setdiff(keys(nt), key)])
@post update(nt, key) = __ret__[key] == val "Value of key updated in result"

"""
Merge `nt1` with `nt2`.
For each key :key_i which clashes (is in both `n1` and `n2`), returns
`f(:Val{:key_i}, nt1.key_i, nt2.key_i)`

```
nt1 = (x = 2, y = 4)
nt2 = (x = 4, a = 2)
mergef(f, nt1, nt2)
```
"""
@generated function mergef(f, nt1::NamedTuple{K1, V1}, nt2::NamedTuple{K2, V2}) where {K1, K2, V1, V2}
  if isempty(K1 ∩ K2)
    :(merge(nt1, nt2))
  else
    :(f(nt1, nt2))
  end
end

@post mergef(f, nt1, nt2) = keys(__ret__) == keys(nt1) ∪ keys(nt2)
@post mergef(f, nt1, nt2) = all((__ret__[k] == nt1[k] for k in keys if k in nt1 && k ∉ nt2))
@post mergef(f, nt1, nt2) = all((__ret__[k] == nt2[k] for k in keys if k in nt2 && k ∉ nt1))
@post mergef(f, nt1, nt2) = all((__ret__[k] == f(nt1[k], nt2[k]) for k in keys if k in nt1 && k ∉ nt2))

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect