text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    let fbody = try __lookup_kwbody__(which(ITensors.addtags, (ITensor{3},String,Vararg{Any,N} where N,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},typeof(addtags),ITensor{3},String,Vararg{Any,N} where N,))
        end
    end
    let fbody = try __lookup_kwbody__(which(ITensors.op, (Index{Int64},String,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},typeof(op),Index{Int64},String,))
        end
    end
    let fbody = try __lookup_kwbody__(which(ITensors.prime, (ITensor{3},Int64,Vararg{Any,N} where N,))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},typeof(prime),ITensor{3},Int64,Vararg{Any,N} where N,))
        end
    end
    let fbody = try __lookup_kwbody__(which(sortperm!, (Array{Int64,1},Array{ITensors.SiteOp,1},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Base.Sort.InsertionSortAlg,Function,Function,Nothing,Base.Order.ForwardOrdering,Bool,typeof(sortperm!),Array{Int64,1},Array{ITensors.SiteOp,1},))
        end
    end
    precompile(Tuple{Core.kwftype(typeof(ITensors.NDTensors.eigen)),NamedTuple{(:ishermitian, :which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho),Tuple{Bool,String,TagSet,Int64,Int64,Float64,Nothing,String}},typeof(eigen),LinearAlgebra.Hermitian{Float64,ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}}}})
    precompile(Tuple{Core.kwftype(typeof(ITensors.factorize)),NamedTuple{(:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho),Tuple{String,TagSet,Int64,Int64,Float64,Nothing,String}},typeof(factorize),ITensor{3},IndexSet{2,Index{Int64}}})
    precompile(Tuple{Core.kwftype(typeof(ITensors.replacebond!)),NamedTuple{(:maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :which_decomp),Tuple{Int64,Int64,Float64,Nothing,String,String}},typeof(replacebond!),MPS,Int64,ITensor{3}})
    precompile(Tuple{Core.kwftype(typeof(KrylovKit.eigsolve)),NamedTuple{(:ishermitian, :tol, :krylovdim, :maxiter),Tuple{Bool,Float64,Int64,Int64}},typeof(KrylovKit.eigsolve),ProjMPO,ITensor{3},Int64,Symbol})
    precompile(Tuple{Core.kwftype(typeof(KrylovKit.expand!)),NamedTuple{(:verbosity,),Tuple{Int64}},typeof(KrylovKit.expand!),KrylovKit.LanczosIterator{ProjMPO,ITensor{3},KrylovKit.ModifiedGramSchmidt2},KrylovKit.LanczosFactorization{ITensor{3},Float64}})
    precompile(Tuple{ProjMPO,ITensor{3}})
    precompile(Tuple{Type{MPO},AutoMPO,Array{Index{Int64},1}})
    precompile(Tuple{Type{TagSet},String})
    precompile(Tuple{typeof(*),Float64,ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}}})
    precompile(Tuple{typeof(*),ITensor{2},ITensor{3}})
    precompile(Tuple{typeof(*),ITensor{3},ITensor{1}})
    precompile(Tuple{typeof(*),ITensor{3},ITensor{2}})
    precompile(Tuple{typeof(*),ITensor{3},ITensor{3}})
    precompile(Tuple{typeof(*),ITensor{4},ITensor{1}})
    precompile(Tuple{typeof(*),ITensor{4},ITensor{4}})
    precompile(Tuple{typeof(==),Array{Index{Int64},1},Array{Index{Int64},1}})
    precompile(Tuple{typeof(==),IndexSet{1,Index{Int64}},IndexSet{1,Index{Int64}}})
    precompile(Tuple{typeof(==),IndexSet{2,Index{Int64}},IndexSet{2,Index{Int64}}})
    precompile(Tuple{typeof(==),IndexSet{3,Index{Int64}},IndexSet{3,Index{Int64}}})
    precompile(Tuple{typeof(==),IndexVal{Index{Int64}},Index{Int64}})
    precompile(Tuple{typeof(==),TagSet,TagSet})
    precompile(Tuple{typeof(Base.Broadcast.materialize!),ITensor{4},Base.Broadcast.Broadcasted{ITensors.ITensorStyle,Nothing,typeof(+),Tuple{ITensor{4},ITensor{4}}}})
    precompile(Tuple{typeof(Base.afoldl),typeof(*),ITensor{4},ITensor{3}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,0,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{0,Union{}}},Tuple{},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,1,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64}}},Tuple{Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,1,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64}}},Tuple{Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contract!!),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64}})
    precompile(Tuple{typeof(ITensors.NDTensors.contraction_output),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},IndexSet{2,Index{Int64}}})
    precompile(Tuple{typeof(ITensors.NDTensors.getperm),IndexSet{4,Index{Int64}},IndexSet{4,Index{Int64}}})
    precompile(Tuple{typeof(ITensors.NDTensors.permutedims!!),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},Function})
    precompile(Tuple{typeof(ITensors.NDTensors.permutedims!!),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},Function})
    precompile(Tuple{typeof(ITensors.compute_contraction_labels),IndexSet{2,Index{Int64}},IndexSet{3,Index{Int64}}})
    precompile(Tuple{typeof(ITensors.compute_contraction_labels),IndexSet{4,Index{Int64}},IndexSet{3,Index{Int64}}})
    precompile(Tuple{typeof(KrylovKit.orthogonalize!),ITensor{3},ITensor{3},KrylovKit.ModifiedGramSchmidt2})
    precompile(Tuple{typeof(Random.randn!),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}}})
    precompile(Tuple{typeof(axpy!),Float64,ITensor{3},ITensor{3}})
    precompile(Tuple{typeof(combiner),Index{Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Number,2,ITensors.NDTensors.Combiner,IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Number,2,ITensors.NDTensors.Combiner,IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},NTuple{4,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Number,2,ITensors.NDTensors.Combiner,IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Number,2,ITensors.NDTensors.Combiner,IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,1,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64}}},Tuple{Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,1,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64}}},Tuple{Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},NTuple{4,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},Tuple{}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,1,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64}}},Tuple{Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,1,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64}}},Tuple{Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},NTuple{4,Int64}})
    precompile(Tuple{typeof(contract),ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64},ITensors.NDTensors.Tensor{Float64,4,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64}}},NTuple{4,Int64}})
    precompile(Tuple{typeof(dmrg),MPO,MPS,Sweeps})
    precompile(Tuple{typeof(hassameinds),ITensor{3},Tuple{Index{Int64},Index{Int64},Index{Int64}}})
    precompile(Tuple{typeof(hastags),TagSet,String})
    precompile(Tuple{typeof(itensor),Array{Float64,3},Index{Int64},Index{Int64},Vararg{Index{Int64},N} where N})
    precompile(Tuple{typeof(mul!),ITensor{3},ITensor{3},Float64})
    precompile(Tuple{typeof(prime),ITensor{3},ITensors.Not{TagSet}})
    precompile(Tuple{typeof(randomITensor),Index{Int64},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(randomMPS),Array{Index{Int64},1},Int64})
    precompile(Tuple{typeof(replaceind!),ITensor{2},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(replaceind),IndexSet{2,Index{Int64}},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(setdiff),IndexSet{2,Index{Int64}},IndexSet{1,Index{Int64}}})
    precompile(Tuple{typeof(setindex!),Dict{Array{ITensors.SiteOp,1},Array{Float64,2}},Array{Float64,2},Array{ITensors.SiteOp,1}})
    precompile(Tuple{typeof(show),Base.GenericIOBuffer{Array{UInt8,1}},TagSet})
    precompile(Tuple{typeof(svd),ITensors.NDTensors.Tensor{Complex{Float64},2,ITensors.NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64}}}})
    precompile(Tuple{typeof(svd),ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64}}}})
end


text:
"""
For an initial guess from the index of a Fold bifurcation point located in ContResult.bifpoint, returns a point which will be refined using `newtonFold`.
"""
function FoldPoint(br::ContResult, index::Int64)
	@assert br.foldpoint[index].type == :fold "The provided index does not refer to a Fold point"
	bifpoint = br.foldpoint[index]
	return BorderedArray(_copy(bifpoint.u), bifpoint.param)
end

####################################################################################################
# Method using Minimally Augmented formulation
@with_kw struct FoldProblemMinimallyAugmented{TF, TJ, TJa, Td2f, vectype, S <: AbstractLinearSolver, Sa <: AbstractLinearSolver, Sbd <: AbstractBorderedLinearSolver}
	F::TF								# Function F(x, p) = 0
	J::TJ								# Jacobian of F wrt x
	Jadjoint::TJa						# Adjoint of the Jacobian of F
	d2F::Td2f = nothing					# Hessian of F
	a::vectype							# close to null vector of J^T
	b::vectype							# close to null vector of J
	linsolver::S						# linear solver
	linsolverAdjoint::Sa = linsolver	# linear solver for the jacobian adjoint
	linbdsolver::Sbd					# bordered linear solver
end

hasHessian(pb::FoldProblemMinimallyAugmented{TF, TJ, TJa, Td2f, vectype, S, Sa, Sbd}) where {TF, TJ, TJa, Td2f, vectype, S, Sa, Sbd} = Td2f != Nothing

FoldProblemMinimallyAugmented(F, J, Ja, a, b, linsolve) = FoldProblemMinimallyAugmented(F, J, Ja, nothing, a, b, linsolve, linsolve, BorderingBLS(linsolve))

FoldProblemMinimallyAugmented(F, J, Ja, d2F, a, b, linsolve) = FoldProblemMinimallyAugmented(F, J, Ja, d2F, a, b, linsolve, linsolve, BorderingBLS(linsolve))

function (fp::FoldProblemMinimallyAugmented)(x::vectype, p::T) where {vectype, T}
	# These are the equations of the minimally augmented (MA) formulation of the Fold bifurcation point
	# input:
	# - x guess for the point at which the jacobian is singular
	# - p guess for the parameter at which the jacobian is singular
	# The jacobian of the MA problem is solved with a bordering method
	a = fp.a
	b = fp.b

	# we solve Jv + a σ1 = 0 with <b, v> = n
	# the solution is v = -σ1 J\a with σ1 = -n/<b, J^{-1}a>
	n = T(1)
	v, _, _ = fp.linsolver(fp.J(x, p), a)
	bv = dot(b, v)
	bv == T(0) && @error "Error with Method using Minimally Augmented formulation for Fold bifurcation"

	σ1 = -n / bv
	return fp.F(x, p), σ1
end

function (foldpb::FoldProblemMinimallyAugmented)(x::BorderedArray)
	res = foldpb(x.u, x.p)
	return BorderedArray(res[1], res[2])
end

# Struct to invert the jacobian of the fold MA problem.
struct FoldLinearSolveMinAug <: AbstractLinearSolver; end

function foldMALinearSolver(x, p::T, pbMA::FoldProblemMinimallyAugmented,
							duu, dup;
							debug_::Bool = false) where T
	# The jacobian should be passed as a tuple as Jac_fold_MA(u0, pb::FoldProblemMinimallyAugmented) = (return (u0, pb, d2F::Bool))
	# The Jacobian J of the vector field is expressed at (x, p)
	# We solve here Jfold⋅res = du := [duu, dup]
	# The Jacobian expression of the Fold problem is Jfold = [J dpF ; σx σp] where σx := ∂_xσ
	# We recall the expression of σx = -< w, d2F(x,p)[v, x2]> where (w, σ2) is solution of J'w + b σ2 = 0 with <a, w> = n
	############### Extraction of function names #################
	# N = length(du) - 1

	F = pbMA.F
	J = pbMA.J
	Jadjoint = pbMA.Jadjoint
	d2F = pbMA.d2F
	a = pbMA.a
	b = pbMA.b

	# we define the following jacobian. It is used at least 3 times below. This avoid doing 3 times the possibly costly building of J(x, p)
	J_at_xp = J(x, p)

	n = T(1)
	# we solve Jv + a σ1 = 0 with <b, v> = n
	# the solution is v = -σ1 J\a with σ1 = -n/<b, J\a>
	v = pbMA.linsolver(J_at_xp, a)[1]
	σ1 = -n / dot(b, v)
	rmul!(v, -σ1)

	# we solve J'w + b σ2 = 0 with <a, w> = n
	# the solution is w = -σ2 J'\b with σ2 = -n/<a, J'\b>
	w = pbMA.linsolverAdjoint(Jadjoint(x, p), b)[1]
	σ2 = -n / dot(a, w)
	rmul!(w, -σ2)

	δ = T(1e-8)
	ϵ1, ϵ2, ϵ3 = T(δ), T(δ), T(δ)
	################### computation of σx σp ####################
	dpF = minus(F(x, p + ϵ1), F(x, p - ϵ1)); rmul!(dpF, T(1) / T(2ϵ1))
	dJvdp = minus(apply(J(x, p + ϵ3), v), apply(J(x, p - ϵ3), v)); rmul!(dJvdp, T(1) / T(2ϵ3))
	σp = -dot(w, dJvdp) / n

	if hasHessian(pbMA) == false
		# We invert the jacobian of the Fold problem when the Hessian of x -> F(x, p) is not known analytically. We thus rely on finite differences which can be slow for large dimensions
		prod(size(x)) > 1e4 && @warn "You might want to pass the Hessian, you are using finite differences with $(prod(size(x))) unknowns to compute a gradient"
		# we first need to compute the value of ∂_xσ written σx
		σx = zero(x)
		e  = zero(x)

		# this is the part which does not work if x is not AbstractArray. We use CartesianIndices to support AbstractArray as a type for the solution we are looking for
		for ii in CartesianIndices(x)
			e[ii] = T(1)
			# d2Fve := d2F(x,p)[v,e]
			d2Fve = (apply(J(x + ϵ2 * e, p), v) - apply(J(x - ϵ2 * e, p), v)) / T(2ϵ2)
			σx[ii] = -dot(w, d2Fve) / n
			e[ii] = T(0)
		end

		########## Resolution of the bordered linear system ########
		dX, dsig, flag, it = pbMA.linbdsolver(J_at_xp, dpF, σx, σp, duu, dup)

	else
		# We invert the jacobian of the Fold problem when the Hessian of x -> F(x, p) is known analytically. Much faster than the previous case

		# we solve it here instead of calling linearBorderedSolver because this removes the need to pass the linear form associated to σx
		x1, x2, _, it = pbMA.linsolver(J_at_xp, duu, dpF)

		d2Fv = d2F(x, p, x1, v)
		σx1 = -dot(w, d2Fv ) / n

		copyto!(d2Fv, d2F(x, p, x2, v))
		σx2 = -dot(w, d2Fv ) / n

		dsig = (dup - σx1) / (σp - σx2)

		# dX = x1 .- dsig .* x2
		dX = copyto!(similar(x1), x1); axpy!(-dsig, x2, dX)
	end

	if debug_
		return dX, dsig, true, sum(it), 0., [J(x, p) dpF ; σx' σp]
	else
		return dX, dsig, true, sum(it), 0.
	end
end

function (foldl::FoldLinearSolveMinAug)(Jfold, du::BorderedArray{vectype, T}, debug_::Bool = false) where {vectype, T}
	out =  foldMALinearSolver(Jfold[1].u,
				 Jfold[1].p,
				 Jfold[2],
				 du.u, du.p;
				 debug_ = debug_)

	if debug_ == false
		return BorderedArray(out[1], out[2]), out[3], out[4], out[5]
	else
		return BorderedArray(out[1], out[2]), out[3], out[4], out[5], out[6]
	end
end

################################################################################################### Newton / Continuation functions
"""
This function turns an initial guess for a Fold point into a solution to the Fold problem based on a Minimally Augmented formulation. The arguments are as follows
- `F   = (x, p) -> F(x, p)` where `p` is the parameter associated to the Fold point
- `dF  = (x, p) -> d_xF(x, p)` associated jacobian
- `dFt = (x, p) -> transpose(d_xF(x, p))` associated jacobian transpose, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method.
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x,p)[v1,v2]`.
- `foldpointguess` initial guess (x_0, p_0) for the Fold point. It should be a `BorderedArray` as given by the function FoldPoint
- `eigenvec` guess for the 0 eigenvector
- `options::NewtonPar`
"""
function newtonFold(F, J, Jt, d2F, foldpointguess::BorderedArray{vectype, T}, eigenvec, options::NewtonPar; normN = norm) where {T, vectype}
	foldvariable = FoldProblemMinimallyAugmented(
		(x, p) ->  F(x, p),
		(x, p) ->  J(x, p),
		(x, p) -> Jt(x, p),
		d2F,
		_copy(eigenvec), #copy(eigenvec),
		_copy(eigenvec), #copy(eigenvec),
		options.linsolver)

	foldPb = u -> foldvariable(u)

	# Jacobian for the Fold problem
	Jac_fold_MA(u0, pb::FoldProblemMinimallyAugmented) = (return (u0, pb))

	opt_fold = @set options.linsolver = FoldLinearSolveMinAug()

	# solve the Fold equations
	return newton(x ->  foldPb(x),
				x -> Jac_fold_MA(x, foldvariable),
				foldpointguess,
				opt_fold, normN = normN)
end

# calls when hessian is unknown, finite differences are then used
newtonFold(F, J, Jt, foldpointguess::BorderedArray{vectype, T}, eigenvec, options::NewtonPar; normN = norm) where {T,vectype} = newtonFold(F, J, Jt, nothing, foldpointguess, eigenvec, options; normN = normN)

newtonFold(F, J, foldpointguess::BorderedArray{vectype, T}, eigenvec, options::NewtonPar; normN = norm) where {T,vectype} = newtonFold(F, J, (x, p) -> transpose(J(x, p)), foldpointguess, eigenvec, options; normN = normN)

"""
Simplified call to refine an initial guess for a Fold point. More precisely, the call is as follows

	`newtonFold(F, J, Jt, br::ContResult, index::Int64, options::NewtonPar)`

or

	`newtonFold(F, J, Jt, d2F, br::ContResult, index::Int64, options::NewtonPar)`

whether the Hessian d2F is known analytically or not. The parameters / options are as usual except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine.
"""
function newtonFold(F, J, Jt, d2F, br::ContResult, ind_fold::Int64, options::NewtonPar;kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.foldpoint[ind_fold]
	eigenvec = bifpt.tau

	# solve the Fold equations
	outfold, hist, flag =  newtonFold(F, J, Jt, d2F, foldpointguess, eigenvec, options; kwargs...)

	return outfold, hist, flag
end

function newtonFold(F, J, Jt, br::ContResult, ind_fold::Int64, options::NewtonPar; kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.foldpoint[ind_fold]
	eigenvec = bifpt.tau

	# solve the Fold equations
	outfold, hist, flag =  newtonFold(F, J, Jt, foldpointguess, eigenvec, options; kwargs...)

	return outfold, hist, flag
end

newtonFold(F, J, br::ContResult, ind_fold::Int64, options::NewtonPar; kwargs...) = newtonFold(F, J, (x, p) -> transpose(J(x, p)), br, ind_fold, options; kwargs...)

newtonFold(F, br::ContResult, ind_fold::Int64, options::NewtonPar; kwargs...) = newtonFold(F, (x0, p) -> finiteDifferences(x -> F(x, p), x0), br, ind_fold, options; kwargs...)

"""
Codim 2 continuation of Fold points. This function turns an initial guess for a Fold point into a curve of Fold points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p1, p2) ->	F(x, p1, p2)` where `p` is the parameter associated to the Fold point
- `J = (x, p1, p2) -> d_xF(x, p1, p2)` associated jacobian
- `Jt = (x, p1, p2) -> transpose(d_xF(x, p1, p2))` associated jacobian transpose
- `d2F = p2 -> ((x, p1, v1, v2) -> d2F(x, p1, p2, v1, v2))` this is the hessian of `F` computed at `(x, p1, p2)` and evaluated at `(v1, v2)`.
- `foldpointguess` initial guess (x_0, p1_0) for the Fold point. It should be a `BorderedArray` as given by the function FoldPoint
- `p2` parameter p2 for which `foldpointguess` is a good guess
- `eigenvec` guess for the 0 eigenvector at p1_0
- `options::NewtonPar`
"""
function continuationFold(F, J, Jt, d2F, foldpointguess::BorderedArray{vectype, T}, p2_0::T, eigenvec, options_cont::ContinuationPar ; kwargs...) where {T,vectype}
	# Bad way it creates a struct for every p2?
	# Jacobian for the Fold problem
	Jac_fold_MA(u0, pb) = (return (u0, pb))

	# options for the Newton Solver inheritated from the ones the user provided
	options_newton = options_cont.newtonOptions

	foldvariable = p2 -> FoldProblemMinimallyAugmented(
		(x, p1) ->  F(x, p1, p2),
		(x, p1) ->  J(x, p1, p2),
		(x, p1) -> Jt(x, p1, p2),
		d2F(p2),
		_copy(eigenvec), #copy(eigenvec),
		_copy(eigenvec), #copy(eigenvec),
		options_newton.linsolver)
	foldPb = (u, p2) -> foldvariable(p2)(u)

	opt_fold_cont = @set options_cont.newtonOptions.linsolver = FoldLinearSolveMinAug()

	# solve the Fold equations
	return continuation(
		(x, p2) -> foldPb(x, p2),
		(x, p2) -> Jac_fold_MA(x, foldvariable(p2)),
		foldpointguess, p2_0,
		opt_fold_cont,
		printSolution = (u, p) -> u.p,
		plotSolution = (x; kwargs...) -> (xlabel!("p2", subplot=1); ylabel!("p1", subplot=1)  ); kwargs...)
end

"""
codim 2 continuation of Fold points. This function turns an initial guess for a Fold point into a curve of Fold points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p1, p2) -> F(x, p1, p2)` where `p` is the parameter associated to the Fold point
- `J = (x, p1, p2) -> d_xF(x, p1, p2)` associated jacobian
- `foldpointguess` initial guess (x_0, p1_0) for the Fold point. It should be a `BorderedArray` as given by the function FoldPoint
- `p2` parameter p2 for which foldpointguess is a good guess
- `eigenvec` guess for the 0 eigenvector at p1_0
- `options::NewtonPar`

!!! warning "Hessian"
    The hessian of `F` in this case is computed with Finite differences. This can be slow for many variables, e.g. ~1e6
"""
continuationFold(F, J, Jt, foldpointguess::BorderedArray{vectype, T}, p2_0::T, eigenvec, options_cont::ContinuationPar ;kwargs...) where {T,vectype} = continuationFold(F, J, Jt, p2 -> nothing, foldpointguess, p2_0, eigenvec, options_cont ; kwargs...)

continuationFold(F, J, foldpointguess::BorderedArray{vectype, T}, p2_0::Real, eigenvec, options_cont::ContinuationPar ; kwargs...)  where {T,vectype} = continuationFold(F, J, (x, p1, p2) -> transpose(J(x, p1, p2)), foldpointguess, p2_0, eigenvec, options_cont ; kwargs...)

function continuationFold(F, foldpointguess::BorderedArray{vectype, T}, p2_0::Real, eigenvec, options::ContinuationPar ; kwargs...)  where {T,vectype}
	return continuationFold(F,
		(x0, p) -> finiteDifferences(x -> F(x, p), x0),
		foldpointguess, p2_0,
		eigenvec,
		options ; kwargs...)
end

"""
Simplified call for continuation of Fold point. More precisely, the call is as follows `continuationFold(F, J, Jt, d2F, br::ContResult, index::Int64, options)` where the parameters are as for `continuationFold` except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine.

Simplified calls are also provided but at the cost of using finite differences.
"""
function continuationFold(F, J, Jt, d2F, br::ContResult, ind_fold::Int64, p2_0::Real, options_cont::ContinuationPar ; kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.foldpoint[ind_fold]
	eigenvec = bifpt.tau
	return continuationFold(F, J, Jt, d2F, foldpointguess, p2_0, eigenvec, options_cont ;kwargs...)
end

function continuationFold(F, J, Jt, br::ContResult, ind_fold::Int64, p2_0::Real, options_cont::ContinuationPar ; kwargs...)
	foldpointguess = FoldPoint(br, ind_fold)
	bifpt = br.foldpoint[ind_fold]
	eigenvec = bifpt.tau
	return continuationFold(F, J, Jt, foldpointguess, p2_0, eigenvec, options_cont; kwargs...)
end

continuationFold(F, J, br::ContResult, ind_fold::Int64, p2_0::Real, options_cont::ContinuationPar ; kwargs...) = continuationFold(F, J, (x, p1, p2) -> transpose(J(x, p1, p2)), br, ind_fold, p2_0, options_cont::ContinuationPar; kwargs...)

continuationFold(F, br::ContResult, ind_fold::Int64, p2_0::Real, options_cont::ContinuationPar ; kwargs...) = continuationFold(F, (x0, p1, p2) -> finiteDifferences(x -> F(x, p1, p2), x0), br, ind_fold, p2_0, options_cont::ContinuationPar; kwargs...)

ext: jl


text:
function solve()
    n = parse(Int, readline())
    ps = []
    for _=1:n
        p = [parse(Int, x) for x in split(readline())]
        push!(ps, p)
    end
    for i=1:n-2, j=i+1:n-1, k=j+1:n
        l = (ps[k][2] - ps[i][2]) * (ps[j][1] - ps[i][1])
        r = (ps[k][1] - ps[i][1]) * (ps[j][2] - ps[i][2])
        if l == r
            return "Yes"
        end
    end
    "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
struct Box <: Widget
    center::Point{2,Float64}
    dims::Point{2,Float64}
    color::RGB
end

GeometryExperiments.PointSet(b::Box) = (Translation(b.center) ∘ Scaling(b.dims ./ 2))(PointSet(HyperCube, Point2f))

vertex_data_type(::Type{Box}) = PosColor{Point2f,RGB{Float16}}
mesh_encoding_type(T::Type{Box}) = MeshVertexEncoding{TriangleStrip, vertex_data_type(T)}

function GeometryExperiments.MeshVertexEncoding(b::Box)
    pos = (Translation(-1.0, -1.0) ∘ inv(Scaling(1920 / 2, 1080 / 2)))(PointSet(b))
    MeshVertexEncoding(collect(vertex_data_type(b).(pos.points, Ref(b.color))), Triangle)
end

nvertices(::Type{Box}) = 4
resource_types(::Type{Box}) = ()


text:
## Funciones para renormalización y cómputo de MAI

## Algoritmos de renormalización de MAI-G

# Función para propósitos ilustrativos. Se utiliza en el cuaderno introductorio de la medida de inflación subyacente MAI. Para realizar el cómputo en InflationCoreMai se utiliza renorm_g_glp_perf.

# Renormalización de distribución g con distribución glp utilizando n segmentos
# Esta función devuelve la distribución glpₜ (tiene menor desempeño)
function renorm_g_glp(G, GLP, glp, n)

    # Percentiles con n segmentos y precisión
    p = (0:n) ./ n
    ε = step(glp.vspace)

    # Obtener los percentiles de los n segmentos
    q_g = quantile(G, p)
    q_glp = quantile(GLP, p)

    # Obtener lista de segmentos para renormalizar
    segments = get_segments_list(q_g, q_glp, n)
    # @debug "Percentiles y segmentos" q_g q_glp ε segments

    # Crear una copia de la distribución glp
    glpₜ = deepcopy(glp)

    # Renormalizar cada segmento
    for i in 2:length(segments)
        k = segments[i]
        t = segments[i-1]

        # Renormalizar el segmento
        v1 = i == 2 ? min(q_g[t], q_glp[t]) : q_g[t]
        vl = i == length(segments) ? max(q_g[k], q_glp[k]) : q_g[k]
        norm_c = (G(vl) - G(v1)) / (GLP(vl) - GLP(v1))

        isnan(norm_c) && continue
        # @debug "Renormalizando segmento $i" t k q_g[t] q_g[k] norm_c

        # Si es el primer segmento, incluir el límite inferior.
        # De lo contrario, renormalizar a partir de la siguiente posición en la
        # grilla
        if i == 2
            renormalize!(glpₜ, v1, vl, norm_c)
        else
            renormalize!(glpₜ, v1 + ε, vl, norm_c)
        end
    end

    # Devolver la distribución renormalizada
    glpₜ
end

# Renormalización de distribución g con distribución glp utilizando n segmentos.
# Esta función no devuelve la distribución glpₜ y aplica algunas optimizaciones
# de memoria
function renorm_g_glp_perf(G::AccumulatedDistr{T}, GLP, glp, q_g, q_glp, n) where T

    # Percentiles con n segmentos y precisión
    ε = step(glp.vspace)

    # Obtener números de segmentos para renormalizar
    r̲, r̄ = get_segments(q_g, q_glp)

    e_glpt::T = zero(T)

    # Renormalizar cada segmento
    @inbounds for i in 2:n+1
        # Obtener segmentos a renormalizar. Saltar en el segmento especial
        r̲ < i < r̄ && continue
        k = i
        t = i == r̄ ? r̲ : i-1

        # Renormalizar el segmento
        if t == 1
            v1 = min(q_g[t], q_glp[t])
        else
            v1 = q_g[t]
        end

        if k == n+1
            vl = max(q_g[k], q_glp[k])
        else
            vl = q_g[k]
        end

        # Constante de normalización
        norm_c = (G(vl) - G(v1)) / (GLP(vl) - GLP(v1))

        (isnan(norm_c) || isinf(norm_c)) && continue

        # Si es el primer segmento, incluir el límite inferior. De lo contrario,
        # renormalizar a partir de la siguiente posición en la grilla
        if t == 1
            e_glpt += renorm_sum(glp, v1, vl, norm_c)
        else
            e_glpt += renorm_sum(glp, v1 + ε, vl, norm_c)
        end
    end

    # Devolver el promedio ponderado de la distribución renormalizada
    e_glpt
end

## Algoritmos de renormalización de variante MAI-F

# Función para propósitos ilustrativos. Se utiliza en el cuaderno introductorio de la medida de inflación subyacente MAI. Para realizar el cómputo en InflationCoreMai se utiliza renorm_g_flp_perf.

# Renormalización de distribución f con distribución flp utilizando n segmentos.
# Esta función devuelve la distribución flpₜ y por lo tanto, tiene menor desempeño.
function renorm_f_flp(F, FLP, GLP, glp, n)

    # Percentiles con n segmentos y precisión
    p = (0:n) ./ n
    ε = step(glp.vspace)

    # Obtener distribuciones acumuladas
    # Obtener los percentiles de los n segmentos
    q_f = quantile(F, p)
    q_flp = quantile(FLP, p)

    # Obtener lista de segmentos para renormalizar
    segments = get_segments_list(q_f, q_flp, n)
    # @debug "Percentiles y segmentos" q_g q_flp ε segments

    # Crear una copia de la distribución glp
    flpₜ = deepcopy(glp)

    # Renormalizar cada segmento
    S = length(segments)

    for i in 2:S
        k = segments[i]
        t = segments[i-1]

        # Renormalizar el segmento
        v1_num = i == 2 ? min(q_f[t], q_flp[t]) : q_flp[t]
        vl_num = i == S ? max(q_f[k], q_flp[k]) : q_flp[k]

        v1 = i == 2 ? min(q_f[t], q_flp[t]) : q_f[t]
        vl = i == S ? max(q_f[k], q_flp[k]) : q_f[k]
        norm_c = (GLP(vl_num) - GLP(v1_num)) / (GLP(vl) - GLP(v1))

        isnan(norm_c) && continue
        # @debug "Renormalizando segmento $i" t k q_g[t] q_g[k] norm_c

        # Si es el primer segmento, incluir el límite inferior.
        # De lo contrario, renormalizar a partir de la siguiente posición en la
        # grilla
        if i == 2
            renormalize!(flpₜ, v1, vl, norm_c)
        else
            renormalize!(flpₜ, v1 + ε, vl, norm_c)
        end
    end

    # Devolver la distribución renormalizada
    flpₜ
end

# Renormalización de distribución f con distribución flp utilizando n segmentos.
# Esta función no devuelve la distribución flpₜ y aplica algunas optimizaciones
# de memoria
function renorm_f_flp_perf(F::AccumulatedDistr{T}, GLP, glp, q_f, q_flp, n) where T

    # Percentiles con n segmentos y precisión
    ε = step(glp.vspace)

    # Obtener números de segmentos para renormalizar
    r̲, r̄ = get_segments(q_f, q_flp)

    e_flpt::T = zero(T)

    # Renormalizar cada segmento
    for i in 2:n+1
        # Obtener segmentos a renormalizar. Saltar en el segmento especial
        r̲ < i < r̄ && continue
        k = i
        t = i == r̄ ? r̲ : i-1

        # Renormalizar el segmento
        if t == 1
            v1 = v1_num = min(q_f[t], q_flp[t])
        else
            v1 = q_f[t]
            v1_num = q_flp[t]
        end

        if k == n+1
            vl = vl_num = max(q_f[k], q_flp[k])
        else
            vl = q_f[k]
            vl_num = q_flp[k]
        end

        # Constante de normalización
        norm_c = (GLP(vl_num) - GLP(v1_num)) / (GLP(vl) - GLP(v1))

        (isnan(norm_c) || isinf(norm_c)) && continue

        # Si es el primer segmento, incluir el límite inferior. De lo contrario,
        # renormalizar a partir de la siguiente posición en la grilla
        if t == 1
            e_flpt += renorm_sum(glp, v1, vl, norm_c)
        else
            e_flpt += renorm_sum(glp, v1 + ε, vl, norm_c)
        end
    end

    # Devolver el promedio ponderado de la distribución renormalizada
    e_flpt
end

## Funciones auxiliares

# Función para renormalizar segmento dado entre dos variaciones intermensuales a y b
# por el valor k
function renormalize!(tdistr::TransversalDistr, a, b, k)
    # Obtener índices de variaciones a y b en vspace
    ia = vposition(a, tdistr.vspace)
    ib = vposition(b, tdistr.vspace)

    # Renormalizar el vector de distribución subyacente
    @views tdistr.distr[ia:ib] .*= k
    nothing
end

# Función para computar la suma ponderada de renormalización del segmento dado
# entre dos variaciones intermensuales a y b por el valor k
function renorm_sum(tdistr::TransversalDistr{T}, a, b, k) where T
    # Obtener índices de variaciones a y b en vspace
    ia = vposition(a, tdistr.vspace)
    ib = vposition(b, tdistr.vspace)

    # Renormalizar el vector de distribución subyacente y obtener la suma ponderada
    s = zero(T)
    @inbounds for j in ia:ib
        s += tdistr.distr[j] * tdistr.vspace[j]
    end

    s * k
end

# Función para obtener segmento especial de renormalización
function get_segments(q_cp, q_lp)
    # Obtener número de percentiles que conforman segmento especial en la
    # disribución de largo plazo
    k̄ = findfirst(q_lp .> 0)
    k̲ = findlast(q_lp .< 0)

    # Obtener número de percentiles que conforman segmento especial en la
    # disribución del mes o ventana
    s̄ = findfirst(q_cp .> 0)
    s̲ = findlast(q_cp .< 0)

    # Obtener los números comunes
    r̲ = min(k̲, s̲)
    r̄ = max(k̄, s̄)

    r̲, r̄
end

function get_segments_list(q_cp, q_lp, n)
    # Devolver lista de segmentos para renormalizar
    r̲, r̄ = get_segments(q_cp, q_lp)
    union(1:r̲, r̄:n+1)
end


text:
addqueen(queens::Array{Vector{Int}}, queen::Vector{Int}) = push!(copy(queens), queen)

hitsany(queen::Vector{Int}, queens::Array{Vector{Int}}) = any(map(x->hits(queen, x), queens))
hits(a::Array{Int}, b::Array{Int}) = any(a .== b) || abs(a-b)[1] == abs(a-b)[2]

function solve(x, y, n, d=Array(Vector{Int}, 0))
  if n == 0
    return d
  end
  for px = 1:x
    for py = 1:y
      if !hitsany([px, py], d)
        s = solve(x, y, n-1, addqueen(d, [px, py]))
        s != nothing && return s
      end
    end
  end
  return nothing
end

max_stars_repo_name: alanedelman/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
@testset "testing copper plate network construction" begin
    ps_model = PSI._canonical_model_init(bus_numbers, GLPK_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PSI.CopperPlatePowerModel, sys5b, time_range);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PSI.CopperPlatePowerModel, sys5b, time_range);
    PSI.construct_network!(ps_model, PSI.CopperPlatePowerModel, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 120
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 24

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)
    @test termination_status(ps_model.JuMPmodel) == MOI.OPTIMAL

    ps_model = PSI._canonical_model_init(bus_numbers, GLPK_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods; parameters = false)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PSI.CopperPlatePowerModel, sys5b, time_range; parameters = false);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PSI.CopperPlatePowerModel, sys5b, time_range; parameters = false);
    PSI.construct_network!(ps_model, PSI.CopperPlatePowerModel, sys5b, time_range; parameters = false);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 120
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 24

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)
    @test termination_status(ps_model.JuMPmodel) == MOI.OPTIMAL
end

@testset "testing DC-PF with PTDF formulation" begin
    PTDF, A = PowerSystems.buildptdf(branches5, nodes5)
    ps_model = PSI._canonical_model_init(bus_numbers, GLPK_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PSI.StandardPTDFForm, sys5b, time_range);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PSI.StandardPTDFForm, sys5b, time_range);
    PSI.construct_network!(ps_model, PSI.StandardPTDFForm, sys5b, time_range; PTDF = PTDF)
    PSI.construct_device!(ps_model, PSY.Branch, PSI.SeriesLine, PSI.StandardPTDFForm, sys5b, time_range)
    @test JuMP.num_variables(ps_model.JuMPmodel) == 264
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 264
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 264

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)

    @test termination_status(ps_model.JuMPmodel) == MOI.OPTIMAL

    ps_model = PSI._canonical_model_init(bus_numbers, GLPK_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods; parameters = false)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PSI.StandardPTDFForm, sys5b, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PSI.StandardPTDFForm, sys5b, time_range; parameters = false)
    PSI.construct_network!(ps_model, PSI.StandardPTDFForm, sys5b, time_range; PTDF = PTDF, parameters = false)
    PSI.construct_device!(ps_model, PSY.Branch, PSI.SeriesLine, PSI.StandardPTDFForm, sys5b, time_range; parameters = false)
    @test JuMP.num_variables(ps_model.JuMPmodel) == 264
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.Interval{Float64}) == 264
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 264

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)

    @test termination_status(ps_model.JuMPmodel) == MOI.OPTIMAL
end

 @testset "PTDF ArgumentError" begin
    ps_model = PSI._canonical_model_init(bus_numbers, GLPK_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PSI.StandardPTDFForm, sys5b, time_range);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PSI.StandardPTDFForm, sys5b, time_range);
    @test_throws ArgumentError PSI.construct_network!(ps_model, PSI.StandardPTDFForm, sys5b, time_range)
end

@testset "testing DC-PF network construction" begin
    ps_model = PSI._canonical_model_init(bus_numbers, GLPK_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PM.DCPlosslessForm, sys5b, time_range);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.DCPlosslessForm, sys5b, time_range);
    PSI.construct_network!(ps_model, PM.DCPlosslessForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 384
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 288

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)

    @test termination_status(ps_model.JuMPmodel) == MOI.OPTIMAL

    ps_model = PSI._canonical_model_init(bus_numbers, GLPK_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods; parameters = false)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PM.DCPlosslessForm, sys5b, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.DCPlosslessForm, sys5b, time_range; parameters = false)
    PSI.construct_network!(ps_model, PM.DCPlosslessForm, sys5b, time_range; parameters = false)
    @test JuMP.num_variables(ps_model.JuMPmodel) == 384
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 288

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)

    @test termination_status(ps_model.JuMPmodel) == MOI.OPTIMAL
end

@testset  "testing AC-PF network construction" begin
    ps_model = PSI._canonical_model_init(bus_numbers, ipopt_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PM.StandardACPForm, sys5b, time_range);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.StandardACPForm, sys5b, time_range);
    PSI.construct_network!(ps_model, PM.StandardACPForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 1056
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 264

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)

    @test termination_status(ps_model.JuMPmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    ps_model = PSI._canonical_model_init(bus_numbers, ipopt_optimizer, PM.AbstractPowerFormulation, 1:sys5b.time_periods; parameters = false)
    PSI.construct_device!(ps_model, PSY.ThermalGen, PSI.ThermalDispatch, PM.StandardACPForm, sys5b, time_range; parameters = false)
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.StandardACPForm, sys5b, time_range; parameters = false)
    PSI.construct_network!(ps_model, PM.StandardACPForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 1056
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 144
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 264

    JuMP.@objective(ps_model.JuMPmodel, Min, AffExpr(0))
    JuMP.optimize!(ps_model.JuMPmodel)

    @test termination_status(ps_model.JuMPmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
end


text:
using DelimitedFiles
using PyPlot
using Statistics
include("compute_elements.jl")

n = 8
mass_grid = readdlm("../../../data/T1_masses_libration_min_max.txt")
fnames = ["T1_xv_libration_min_1.txt","T1_xv_libration_max_2.txt","T1_positions_velocities0001.txt"]

ifile = 1
fname = string("../../../data/",fnames[ifile])
t,eloft1 = compute_elements(fname,[1.0;mass_grid[ifile,:]],n)
ifile = 2
fname = string("../../../data/",fnames[ifile])
t,eloft2 = compute_elements(fname,[1.0;mass_grid[ifile,:]],n)
ifile = 3
fname = string("../../../data/",fnames[ifile])
t,eloft3 = compute_elements(fname,[1.0;mass_grid[ifile,:]],n)
#mass_grid = [ 4.5651630946375855e-5, 4.509775588066495e-5, 1.2793574298379903e-5, 2.2636788066718058e-5, 3.4259988884100445e-5, 4.390338504975531e-5, 9.898177265368541e-6]
#t,eloft = compute_elements("Hyak/Fabrycky/T1_positions_velocities0001.txt",[1.0;mass_grid],n)

function normalize(phi)
while phi < 0
  phi += 2pi
end
while phi > 2pi
  phi -= 2pi
end
return phi
end

function offset(phi,nt)
for i=1:nt-1
  if phi[i+1] > (phi[i]+3.0)
    phi[i+1] -= pi
  end
end
return phi
end

nskip = 100
rtd = 180/pi
cp = ["C0","C1","C2","C3","C4","C5","C6","C7"]
nt = length(t)
phi11 = normalize.(2eloft1[:,2,1] .- 5eloft1[:,2,2] + 3eloft1[:,2,3])
phi11 = offset(phi11,nt)
phi12 = normalize.(2eloft2[:,2,1] .- 5eloft2[:,2,2] + 3eloft2[:,2,3])
phi12 = offset(phi12,nt)
phi13 = normalize.(2eloft3[:,2,1] .- 5eloft3[:,2,2] + 3eloft3[:,2,3])
phi13 = offset(phi13,nt)

clf();
figure(figsize=(8,8))
#plot(t[1:nskip:end],360 .- phi1[1:nskip:end].*rtd,label=L"360-$\phi_{bcd}$",color=cp[1],alpha=0.5);
plot(t[1:nskip:end], phi11[1:nskip:end].*rtd,label=L"$\phi_{bcd}$",color=cp[1],alpha=0.5);
plot(t[1:nskip:end], phi12[1:nskip:end].*rtd,color=cp[1],alpha=0.5);
plot(t[1:nskip:end], phi13[1:nskip:end].*rtd,color=cp[1],alpha=0.5);
plot([t[1],t[end]],[0.0,0.0] .+ (360 - 196.5),color=cp[1],linestyle="--");
phi21 = normalize.(eloft1[:,2,2] .- 3eloft1[:,2,3] + 2eloft1[:,2,4] .+8pi);
phi21 = offset(phi21,nt);
phi22 = normalize.(eloft2[:,2,2] .- 3eloft2[:,2,3] + 2eloft2[:,2,4] .+8pi);
phi22 = offset(phi22,nt);
phi23 = normalize.(eloft3[:,2,2] .- 3eloft3[:,2,3] + 2eloft3[:,2,4] .+8pi);
phi23 = offset(phi23,nt);
#plot(t[1:nskip:end],360 .- phi2[1:nskip:end].*rtd,label=L"360-$\phi_{cde}$",color=cp[2],alpha=0.5);
plot(t[1:nskip:end],phi21[1:nskip:end].*rtd,label=L"$\phi_{cde}$",color=cp[2],alpha=0.5);
plot(t[1:nskip:end],phi22[1:nskip:end].*rtd,color=cp[2],alpha=0.5);
plot(t[1:nskip:end],phi23[1:nskip:end].*rtd,color=cp[2],alpha=0.5);
plot([t[1],t[end]],[0.0,0.0] .+ (360-308.8),color=cp[2],linestyle="--");
phi31 = normalize.(2eloft1[:,2,3] .- 5eloft1[:,2,4] + 3eloft1[:,2,5] .-4pi);
phi31 = offset(phi31,nt);
phi32 = normalize.(2eloft2[:,2,3] .- 5eloft2[:,2,4] + 3eloft2[:,2,5] .-4pi);
phi32 = offset(phi32,nt);
phi33 = normalize.(2eloft3[:,2,3] .- 5eloft3[:,2,4] + 3eloft3[:,2,5] .-4pi);
phi33 = offset(phi33,nt);
plot(t[1:nskip:end],phi31[1:nskip:end].*rtd,label=L"$\phi_{def}$",color=cp[3],alpha=0.5);
plot(t[1:nskip:end],phi32[1:nskip:end].*rtd,color=cp[3],alpha=0.5);
plot(t[1:nskip:end],phi33[1:nskip:end].*rtd,color=cp[3],alpha=0.5);
plot([t[1],t[end]],[0.0,0.0] .+ 219.7,color=cp[3],linestyle="--");
phi41 = normalize.(eloft1[:,2,4] .- 3eloft1[:,2,5] + 2eloft1[:,2,6] .+2pi);
phi41 = offset(phi41,nt);
phi42 = normalize.(eloft2[:,2,4] .- 3eloft2[:,2,5] + 2eloft2[:,2,6] .+2pi);
phi42 = offset(phi42,nt);
phi43 = normalize.(eloft3[:,2,4] .- 3eloft3[:,2,5] + 2eloft3[:,2,6] .+2pi);
phi43 = offset(phi43,nt);
plot(t[1:nskip:end],phi41[1:nskip:end].*rtd,label=L"$\phi_{efg}$",color=cp[4],alpha=0.5);
plot(t[1:nskip:end],phi42[1:nskip:end].*rtd,color=cp[4],alpha=0.5);
plot(t[1:nskip:end],phi43[1:nskip:end].*rtd,color=cp[4],alpha=0.5);
plot([t[1],t[end]],[0.0,0.0].+ 282.3,color=cp[4],linestyle="--");
phi51 = normalize.(eloft1[:,2,5] .- 2eloft1[:,2,6] +  eloft1[:,2,7] .+2pi);
phi51 = offset(phi51,nt);
phi52 = normalize.(eloft2[:,2,5] .- 2eloft2[:,2,6] +  eloft2[:,2,7] .+2pi);
phi52 = offset(phi52,nt);
phi53 = normalize.(eloft3[:,2,5] .- 2eloft3[:,2,6] +  eloft3[:,2,7] .+2pi);
phi53 = offset(phi53,nt);
plot(t[1:nskip:end],phi51[1:nskip:end].*rtd,label=L"$\phi_{fgh}$",color=cp[5],alpha=0.5);
plot(t[1:nskip:end],phi52[1:nskip:end].*rtd,color=cp[5],alpha=0.5);
plot(t[1:nskip:end],phi53[1:nskip:end].*rtd,color=cp[5],alpha=0.5);
plot([t[1],t[end]],[0.0,0.0] .+180.4,color=cp[5],linestyle="--");
legend(ncol=3,fontsize=20);
xlabel("Time [days]",fontsize=15);
ylabel(L"$\phi$ [$^\circ$]",fontsize=15);
;
println("bcd: ",mean(phi11)*rtd," ",mean(phi12)*rtd," ",mean(phi13)*rtd);
println("cde: ",mean(phi21)*rtd," ",mean(phi22)*rtd," ",mean(phi23)*rtd);
println("def: ",mean(phi31)*rtd," ",mean(phi32)*rtd," ",mean(phi33)*rtd);
println("efg: ",mean(phi41)*rtd," ",mean(phi42)*rtd," ",mean(phi43)*rtd);
println("fgh: ",mean(phi51)*rtd," ",mean(phi52)*rtd," ",mean(phi53)*rtd);
axis([t[1],t[end],0,360]);
tight_layout()
savefig("../Laplace_angle.pdf",bbox_inches="tight")


text:
using Luxor

"""
    drawlogo(logo_width::Int)

Draw a `logo_width × logo_width`-sized square image logo for JuliaConstraints.
"""
function drawlogo(logo_width::Int)

    base_width = 1024
    s = logo_width/base_width

    origin() # Set origin to center of image
    scale(s) # Call to scale() must come after origin()

    triangle_side = 0.73*base_width
    triangle_height = triangle_side*√3/2

    ball_centers = [Point(-triangle_side/2,triangle_height/2), Point(0,-triangle_height/2), Point(triangle_side/2,triangle_height/2)]
    ball_radius = 0.18*triangle_side
    ball_colors = ["brown3", "forestgreen", "mediumorchid3"]

    link_num = 7 # number of chain links in each direction
    link_small_width = triangle_side/48
    link_large_width = triangle_side/24
    link_thickness = 2*link_small_width*s # line thickness is not scaled http://juliagraphics.github.io/Luxor.jl/stable/transforms/#Scaling-of-lines
    link_length = triangle_side/link_num/2 + link_small_width/2

    link_color = "royalblue"

    @layer begin # Draw rounded rectangles for the "large chain links"
        setcolor(link_color)
        angles = (-π/3, π/3, 0)
        for i in 1:3
            pt1 = ball_centers[mod1(i, 3)]
            pt2 = ball_centers[mod1(i + 1, 3)]
            for j in 1:link_num
                if mod(j,2) == 1
                    link_shape = squircle(O, link_length, link_small_width, rt=0.5, vertices=true) # no drawing, just build array with shape
                    style = :fill
                else
                    link_shape = squircle(O, link_length, link_large_width, rt=0.5, vertices=true) # no drawing, just build array with shape
                    style = :stroke
                end
                link_large = poly(link_shape, :transparent) # transparent => no drawing
                polyrotate!(link_large, angles[i]) # rotate transparent polygonal, so no drawing
                polymove!(link_large, O, between(pt1, pt2, j/(link_num+1))) # move transparent polygonal, still no drawing
                setline(link_thickness)
                poly(link_large, style) # finally draws it, in the right place
            end
        end
    end

    @layer begin # Draw balls
        for (ball_center, ball_color) in zip(ball_centers, ball_colors)
            setcolor(ball_color)
            circle(ball_center, ball_radius, :fill)
        end
    end
end

for logo_width in (1024, 512, 128, 32) # large, medium, small and favicon sizes
    Drawing(logo_width, logo_width, joinpath("logo_images", "jc-logo_"*string(logo_width)*"x"*string(logo_width)*".png"))
    drawlogo(logo_width)
    finish()
end

max_stars_repo_name: JuliaConstraints/jc-logo

max_issues_repo_name: JuliaConstraints/jc-logo

max_forks_repo_name: JuliaConstraints/jc-logo


text:
test = """
    ecl:gry pid:860033327 eyr:2020 hcl:#fffffd
    byr:1937 iyr:2017 cid:147 hgt:183cm

    iyr:2013 ecl:amb cid:350 eyr:2023 pid:028048884
    hcl:#cfa07d byr:1929

    hcl:#ae17e1 iyr:2013
    eyr:2024
    ecl:brn pid:760753108 byr:1931
    hgt:179cm

    hcl:#cfa07d eyr:2025 pid:166559648
    iyr:2011 ecl:brn hgt:59in
    """

function read_passport(ppt::AbstractString)::Dict{Symbol, <: AbstractString}
    return Dict(let (key, value) = split(field, ':')
                    Symbol(key) => value
                end
                for field in split(ppt, r"\s+"; keepempty = false))
end

function read_passports(io)::Vector
    return read_passport.(split(read(io, String), "\n\n"))
end

const REQ_FIELDS = [:byr, :iyr, :eyr, :hgt, :hcl, :ecl, :pid]

q1(ppt::AbstractDict)::Bool = REQ_FIELDS ⊆ keys(ppt)

byr(v) = 1920 <= parse(Int, v) <= 2002
iyr(v) = 2010 <= parse(Int, v) <= 2020
eyr(v) = 2020 <= parse(Int, v) <= 2030
hcl(v) = length(v) == 7 && v[1] == '#' && all(∈(['0':'9'; 'a':'f']), v[2:7])
ecl(v) = Symbol(v) ∈ [:amb, :blu, :brn, :gry, :grn, :hzl, :oth]
pid(v) = length(v) == 9 && tryparse(Int, v) !== nothing
cid(v) = true

function hgt(s)::Bool
    try
        v, u = parse(Int, s[1:end - 2]), Symbol(s[end - 1:end])
        return u == :in && 59 <= v <= 76 ||
            u == :cm && 150 <= v <= 193
    catch
        return false
    end
end

function q2(ppt::AbstractDict)::Bool
    return q1(ppt) && all(ppt) do (key, value)
        return getfield(@__MODULE__, key)(value)
    end
end

let passports = read_passports(IOBuffer(test))
    @assert count(q1, passports) == 2
    @assert count(q2, passports) == 2
end

passports = read_passports("day04.in")
println("Q1: ", count(q1, passports))
println("Q2: ", count(q2, passports))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using MyFunRep
using Test
using Plots

@testset "MyFunRep.jl" begin
    # Write your tests here.
end


text:
export sigmoid

"""
    Sigmoid <: Func
"""
@func mutable struct Sigmoid end

forward(f::Sigmoid, x) = @. tanh(x * 0.5) * 0.5 + 0.5

function backward(self, gy)
    y = self.outputs[1].value
    return gy .* y .* (1 - y)
end

sigmoid(x) = Sigmoid()(x)


text:
# Type
"""
	(type) Segment{T<:AnyComplex} in the complex plane

Each `Segment` type is parameterized according to the common type of its complex input arguments.
"""
struct Segment{T<:AnyComplex} <: AbstractCurve
	za::T
	zb::T
	function Segment{T}(a,b) where T<:AnyComplex
		@assert isfinite(a) && isfinite(b)
		new(a,b)
	end
end

# Untyped constructor
"""
	Segment(a,b)

Consruct a segment that starts at value `a` and ends at `b`.
"""
function Segment(a::Number,b::Number)
	a,b = promote(complex(float(a)),b)
	Segment{typeof(a)}(a,b)
end

# Complex type converters
for ctype in [:Spherical,:Polar,:Complex]
	@eval begin
		function $ctype(S::Segment{T}) where T<:AnyComplex
			Segment($ctype(S.za),$ctype(S.zb))
		end
	end
end

# Required methods
arclength(S::Segment) = abs(S.zb-S.za)
point(S::Segment,t::Real) = (1-t)*S.za + t*S.zb
(C::Segment)(t::Real) = point(C,t)
"""
	arg(S::Segment,z)

Find the parameter argument `t` such that `S(t)==z` is true.

This gives undefined results if `z` is not actually on the segment.
"""
arg(S::Segment,z::Number) = (real(z) - real(S.za)) / (real(S.zb) - real(S.za))
tangent(S::Segment,t::Real) = S.zb - S.za
unittangent(S::Segment,t::Real=0) = sign(S.zb-S.za)

+(S::Segment,z::Number) = Segment(S.za+z,S.zb+z)
-(S::Segment) = Segment(-S.za,-S.zb)
*(S::Segment,z::Number) = Segment(S.za*z,S.zb*z)

"""
	inv(S)
Invert the segment `S` through the origin. In general the inverse is an `Arc`, though the result is a `Segment` if `S` would pass through the origin when extended.
"""
function inv(S::Segment)
	w = 1 ./ point(S,[0,0.5,1])
	Arc(w...)
end

# Other methods
isfinite(::Segment) = true
conj(S::Segment) = Segment(conj(S.za),conj(S.zb))
reverse(S::Segment) = Segment(S.zb,S.za)

sign(S::Segment) = unittangent(S)

"""
	isapprox(S1::Segment,S2::Segment; tol=<default>)
	S1 ≈ S2

Determine if `S1` and `S2` represent the same segment, irrespective of the type or values of its parameters. Identity is determined by agreement within `tol`, which is interpreted as the weaker of absolute and relative differences.
"""
function isapprox(S1::Segment,S2::Segment;tol=DEFAULT[:tol])
	return isapprox(S1.za,S2.za,rtol=tol,atol=tol) &&
		isapprox(S1.zb,S2.zb,rtol=tol,atol=tol)
end

"""
	isleft(z,S::Segment)

Determine whether the number `z` lies "to the left" of segment `S`. This means that the angle it makes with `tangent(S)` is in the interval (0,π).

Note that `isleft` and `isright` are *not* logical opposites; a point on the (extended) segment should give `false` in both cases.
"""
function isleft(z::Number,S::Segment)
	a,b = S.za,S.zb
	(real(b)-real(a)) * (imag(z)-imag(a)) > (real(z)-real(a)) * (imag(b)-imag(a))
end
"""
	isright(z,S::Segment)

Determine whether the number `z` lies "to the right" of segment `S`. This means that the angle it makes with `tangent(S)` is in the interval (-π,0).

Note that `isleft` and `isright` are *not* logical opposites; a point on the (extended) segment should give `false` in both cases.
"""
function isright(z::Number,S::Segment)
	a,b = S.za,S.zb
	(real(b)-real(a)) * (imag(z)-imag(a)) < (real(z)-real(a)) * (imag(b)-imag(a))
end

"""
	dist(z,S::Segment)

Compute the distance from number `z` to the segment `S`.
"""
dist(z::Number,S::Segment) = abs(z - closest(z,S))
"""
	closest(z,S::Segment)

Find the point on segment `S` that lies closest to `z`.
"""
function closest(z::Number,S::Segment)
	# translate and rotate segment to positive Re axis
	d = S.zb-S.za
	s = sign(d)
	ζ = (z-S.za)/s
	S.za + s*min( max(real(ζ),0), abs(d) )
end
"""
	reflect(z,S::Segment)

Reflect the value `z` across the extension of segment `S` to a line. (For reflection of a segment through a point, use translation and negation.)
"""
reflect(z::Number,S::Segment) = reflect(z,Line(S.za,S.zb))

# Display methods
function show(io::IO,S::Segment{T}) where {T}
	print(IOContext(io,:compact=>true),"Segment(",point(S,0),",",point(S,1),")")
end

function show(io::IO,::MIME"text/plain",S::Segment{T}) where {T}
	print(io,"Segment{$T} in the complex plane:\n   from (",point(S,0),") to (",point(S,1),")")
end

plotdata(S::Segment{T}) where T<:Union{Complex,Polar} = [S.za,S.zb]


text:
using Parquet
using Test

function get_parcompat(parcompat=joinpath(dirname(@__FILE__), "parquet-compatibility"))
    # look for parquet-compatibility in test folder, clone to tempdir if not found
    if !isdir(parcompat)
        parcompat = joinpath(dirname(@__FILE__), "parquet-compatibility")
        run(`git clone https://github.com/Parquet/parquet-compatibility.git $parcompat`)
    end
end

get_parcompat()


text:
module MathOptInterfaceIpopt

# todo
# single get/set
# fix chg coeff
# inplace getters
# updates!

export IpoptSolverInstance

const IPOPTMOI = MathOptInterfaceIpopt

import Base.show, Base.copy

using Ipopt
using MathOptInterface
const MOI = MathOptInterface

const SUPPORTED_OBJECTIVES = [
]
const SUPPORTED_CONSTRAINTS = [
]
const SUPPORTED_VARIABLE = [
]

end # module


text:
"""
    TwosArray{N, T, L}

Represents a static array with N dimensions whose size is exactly 2 along each dimension.
This makes templating on the number of dimensions easier than with a regular SArray.
"""
struct TwosArray{N, T, L} <: StaticArray{NTuple{N, 2}, T, N}
    data::NTuple{L, T}
end

function TwosArray(x::NTuple{L, T}) where {L, T}
    @assert log2(L) == round(log2(L)) "L must equal 2^N for integer N"
    N = Int(log2(L))
    TwosArray{N, T, L}(x)
end

getindex(b::TwosArray, i::Int) = b.data[i]

similar_type(::Type{TwosArray{N, T, L}}, ::Type{T2}) where {N, T, L, T2} = TwosArray{N, T2, L}


text:
pwd() ∉ LOAD_PATH && push!(LOAD_PATH, pwd())
"C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src" ∉ LOAD_PATH && push!(LOAD_PATH, "C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src")
"C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src\\3DFDTD" ∉ LOAD_PATH && push!(LOAD_PATH, "C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\src\\3DFDTD")

"C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\tests" ∉ LOAD_PATH && push!(LOAD_PATH, "C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\tests")
#"/home/zander/Code/julia/RCS-Invert/src" ∉ LOAD_PATH && push!(LOAD_PATH, "/home/zander/Code/julia/RCS-Invert/src")
#"/home/zander/Code/julia/RCS-Invert/src" ∉ LOAD_PATH && push!(LOAD_PATH, "/home/zander/Code/julia/RCS-Invert/src")
using Test
using Plots
using Yee

dB(x::Number) = 20*log10(abs(x))
include("C:\\Users\\Jones\\Desktop\\home\\Code\\julia\\RCS-Invert\\tests\\fdtd_test.jl")

detectorindex = (1, 2, [:, :, :]) # detectorindex follows the form (Field, which component to read, [location of detector])
typeof(detectorindex)
field, field_extended = PML_Test(detectorindex)

plotly()
heatmap(dB.(field[100][1][2][:,:,div(54,2)] .-
        field_extended[100][1][2][div(107,2)-div(54,2)+1:div(107,2)+div(54,2),div(107,2)-div(54,2)+1:div(107,2)+div(54,2),div(107,2)]))

heatmap(dB.(field_extended[100][1][2][:,div(107,2),:]))

axislen = size(FDTD_out[1][1][1], 1)
timelen = size(FDTD_out, 1)

lightcone = zeros(axislen, timelen-50)
for i in 1:timelen-50
    lightcone[:,i] .= FDTD_out[i][1][2][div(54,2),div(54,2),:]
end
heatmap(dB.(lightcone))
#plot((wierd))
#@test !any(isnan.(FDTD_out))
#@test !any(isinf.(FDTD_out))
#=
pmllightcone = zeros(size(FDTD_out[1].E.min[1][2],1), timelen)
for i in 1:timelen
    pmllightcone[:,i] = FDTD_out[i].E.max[1][3][:,div(54,2),div(54,2)]
end
heatmap(dB.(pmllightcone))
=#


text:
"""
    SGLD(n_iters::Int, epsilon::Float64)

 Stochastic Gradient Langevin Dynamics sampler.

Usage:

```julia
SGLD(1000, 0.5)
```

Arguments:

- `n_iters::Int` : Number of samples to pull.
- `epsilon::Float64` : The scaling factor for the learing rate.

Example:

```julia
@model example begin
  ...
end

sample(example, SGLD(1000, 0.5))
```
"""
mutable struct SGLD{AD, T} <: StaticHamiltonian{AD}
    n_iters :: Int       # number of samples
    epsilon :: Float64   # constant scale factor of learning rate
    space   :: Set{T}    # sampling space, emtpy means all
    gid     :: Int
end
SGLD(args...; kwargs...) = SGLD{ADBackend()}(args...; kwargs...)
function SGLD{AD}(epsilon::Float64, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    SGLD{AD, eltype(_space)}(1, epsilon, _space, 0)
end
function SGLD{AD}(n_iters, epsilon) where AD
    SGLD{AD, Any}(n_iters, epsilon, Set(), 0)
end
function SGLD{AD}(n_iters, epsilon, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return SGLD{AD, eltype(_space)}(n_iters, epsilon, _space, 0)
end
function SGLD{AD1}(alg::SGLD{AD2, T}, new_gid::Int) where {AD1, AD2, T}
    SGLD{AD1, T}(alg.n_iters, alg.epsilon, alg.space, new_gid)
end
function SGLD{AD, T}(alg::SGLD, new_gid::Int) where {AD, T}
    SGLD{AD, T}(alg.n_iters, alg.epsilon, alg.space, new_gid)
end

function step(model, spl::Sampler{<:SGLD}, vi::VarInfo, is_first::Val{true})
    spl.alg.gid != 0 && link!(vi, spl)

    spl.info[:wum] = NaiveCompAdapter(UnitPreConditioner(), ManualSSAdapter(MSSState(spl.alg.epsilon)))

    # Initialize iteration counter
    spl.info[:t] = 0

    spl.alg.gid != 0 && invlink!(vi, spl)
    return vi, true
end

function step(model, spl::Sampler{<:SGLD}, vi::VarInfo, is_first::Val{false})
    # Update iteration counter
    spl.info[:t] += 1

    @debug "compute current step size..."
    γ = .35
    ϵ_t = spl.alg.epsilon / spl.info[:t]^γ # NOTE: Choose γ=.55 in paper
    mssa = spl.info[:wum].ssa
    mssa.state.ϵ = ϵ_t

    @debug "X-> R..."
    if spl.alg.gid != 0
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    @debug "recording old variables..."
    θ = vi[spl]
    _, grad = gradient(θ, vi, model, spl)
    verifygrad(grad)

    @debug "update latent variables..."
    θ .-= ϵ_t .* grad ./ 2 .+ rand.(Normal.(zeros(length(θ)), sqrt(ϵ_t)))

    @debug "always accept..."
    vi[spl] = θ

    @debug "R -> X..."
    spl.alg.gid != 0 && invlink!(vi, spl)

    return vi, true
end


text:
using BARON
BARON_SILENT = with_optimizer(BARON.Optimizer, OutputFlag = 0)

# Loading BARON examples (that haven't been loaded yet)
for i in ["nlp1.jl", "nlp2.jl", "nlp3.jl"]
    include(OCTHaGOn.BARON_DIR * i)
end

BARON_PROBLEMS = [minlp, pool1, nlp1, nlp2, nlp3]

bms = [p(false) for p in BARON_PROBLEMS]
globalsolve_and_time!.(bms)

gms = [p(true) for p in BARON_PROBLEMS]
globalsolve_and_time!.(gms)

# save_fit.(gms)
# load_fit.(gms)


text:
function projection_onto_bounds_ipopt_style( nlp::Class_CUTEst, pars::Class_parameters, x::Array{Float64,1} )
    if pars.output_level >= 4
        println("projecting onto bounds ...")
    end
    ifree = _i_not_fixed(nlp.nlp)
    uvar = deepcopy(nlp.nlp.meta.uvar[ifree])
    lvar =  deepcopy(nlp.nlp.meta.lvar[ifree])

    x = deepcopy(x)
    b_L = zeros(length(x))
    b_U = zeros(length(x))

    κ_1 = 1e-2;
    κ_2 = 1e-2;

    threshold_too_close = 1e-8

    for i = 1:length(x)
      if (uvar[i] < lvar[i] + threshold_too_close)
        @show uvar[i] - lvar[i]
        println("ERROR: bounds too close!")
        error("bounds too close!")
      end

      p_L = min(κ_1 * max(1, abs(lvar[i])), κ_2 * (uvar[i] - lvar[i]) )
      p_U = min(κ_1 * max(1, abs(uvar[i])), κ_2 * (uvar[i] - lvar[i]) )

      if lvar[i] > -Inf
        b_L[i] = lvar[i] + p_L
      else
        b_L[i] = -Inf
      end

      if uvar[i] < Inf
        b_U[i] = uvar[i] - p_U
      else
        b_U[i] = Inf
      end

      x[i] = x[i] #+ rand() * 100.0

      if b_L[i] < b_U[i]
        x[i] = min(max(b_L[i], x[i]), b_U[i])
      else
        @show b_L[i], b_U[i]
        error("this shouldn't happen!")
        #x[i] = (uvar[i] + lvar[i]) / 2.0
      end

      if(x[i] <= lvar[i] || x[i] >= uvar[i])
         @show (b_L[i], b_U[i]), (lvar[i], uvar[i]) , x[i]
         error("Error in projection_onto_bounds")
      end
    end

    return x
end


text:
# issue 288
function not_a_reduction!(A, B)
  @turbo for j ∈ 1:size(A.re, 2)
    jre = B.re[j]
    jim = B.im[j]
    for i ∈ 1:size(A.re, 1)
      ire = B.re[i]
      iim = B.im[i]
      cisim = iim * jre - ire * jim
      cisre = ire * jre + iim * jim
      ρre_i = A.re[i,j]
      ρim_i = A.im[i,j]
      re_out = ρre_i * cisre - ρim_i * cisim
      im_out = ρre_i * cisim + ρim_i * cisre
      A.re[i,j] = re_out
      A.im[i,j] = im_out
    end
  end
  return nothing
end
function not_a_reduction_noturbo!(A, B)
  @turbo for j ∈ 1:size(A.re, 2)
    jre = B.re[j]
    jim = B.im[j]
    for i ∈ 1:size(A.re, 1)
      ire = B.re[i]
      iim = B.im[i]
      cisim = iim * jre - ire * jim
      cisre = ire * jre + iim * jim
      ρre_i = A.re[i,j]
      ρim_i = A.im[i,j]
      re_out = ρre_i * cisre - ρim_i * cisim
      im_out = ρre_i * cisim + ρim_i * cisre
      A.re[i,j] = re_out
      A.im[i,j] = im_out
    end
  end
  return nothing
end

@testset "Untangle reductions" begin
  @show @__LINE__
  N = 11
  A1 = (re = rand(N,N), im = rand(N,N))
  A2 = deepcopy(A1)
  B = (re = rand(N), im = rand(N))
  not_a_reduction!(A1, B)
  not_a_reduction_noturbo!(A2, B)
  @test A1.re ≈ A2.re
  @test A1.im ≈ A2.im
end


text:
const GAS_CONSTANT_SI = 8.31446261815324 # Exact since May20, 2019; https://en.wikipedia.org/wiki/Gas_constant
const PA_IN_BAR = 1e5               # Сколько паскалей в баре
const BOLTZMANN_SI = 1.380649e-23 # Exact since May 20, 2019; https://en.wikipedia.org/wiki/Boltzmann_constant
const C3_TC = 369.85 # propane critical temperature (K)
const C3_VC = 2e-4 # propane critical volume (m³ / mol)
const C3_ZC = 0.27626 # propane critical Z-factor (dimensionless)
const C3_OMEGA = 0.1524 # propane acentric factor


text:
using Random
using LinearAlgebra
using Plots
using CSV, DataFrames
using RandomizedLasso
const RL = RandomizedLasso

## Load data
file = CSV.read("/Users/theodiamandis/Downloads/file7b5323e77330.csv", DataFrame)
M = Matrix(file)
b = M[:, 1]
A = M[:, 2:end]
m, n = size(A)
λmax = maximum(1/m*A'*b)
λ = 0.01*λmax

## Solve ADMM
prob = RL.LassoADMMProblem(A, b, λ)

RL.reset!(prob)
res_npc = RL.solve!(
    prob;
    print_iter=10,
    precondition=false,
    max_iters=200,
    tol=1e-5,
    logging=true
)

RL.reset!(prob)
res = RL.solve!(
    prob;
    print_iter=10,
    precondition=true,
    max_iters=200,
    tol=1e-5,
    logging=true
)

## Plotting
log_npc = res_npc.log
log = res.log

function make_plot(y1, y2, yaxis_label; yaxis=:log)
    plt = plot(
        log_npc.iter_time,
        y1,
        dpi=300,
        lw=3,
        label="No Preconditioner",
        ylabel=yaxis_label,
        xlabel="Wall clock time (s)",
        title="$yaxis_label Convergence",
        legend=:topright,
        yaxis=yaxis
    )
    plot!(plt, log.iter_time, y2, label="Nystrom", lw=3)
    return plt
end

plt_rp = make_plot(log_npc.rp, log.rp, "Primal Residual")
plt_rd = make_plot(log_npc.rd, log.rd, "Dual Residual")
plt_dual_gap = make_plot(log_npc.dual_gap, log.dual_gap, "Duality Gap")
plt_obj_val =make_plot(log_npc.obj_val, log.obj_val, "Objective Value"; yaxis=:none)
savefig(plt_rp, joinpath(@__DIR__, "figs/rp.pdf"))
savefig(plt_rd, joinpath(@__DIR__, "figs/rd.pdf"))
savefig(plt_dual_gap, joinpath(@__DIR__, "figs/dg.pdf"))
savefig(plt_obj_val, joinpath(@__DIR__, "figs/obj.pdf"))


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

# module Vehicles

# export plot_vehicles!

using Agents
using InteractiveDynamics
using GLMakie
using Observables

include("Environment.jl")
include("Parameters.jl")
include("DrawEnvironment.jl")

import .Parameters as P

mutable struct VehicleAgent <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
end

function initialize()
    space2d = ContinuousSpace((P.EXTENT_WIDTH, P.EXTENT_HEIGHT))
    properties = Dict()
    properties[:tick] = 0
    properties[:steptext] = Observable("Step: " * string(properties[:tick]))
    intersectingRoads = TwoWayIntersectingRoads(2000, 0, 4000, 2000, 4000, 0)
    properties[:env] = intersectingRoads
    properties[:spawn_rate] = 1400
    model = ABM(VehicleAgent, space2d, scheduler=Schedulers.randomly; properties=properties)
    add_vehicle!((4000.0, 1924.0), model, (-0.1, 0.0))
    return model
end

function add_vehicle!(spawn_pos, model, initial_vel=(P.VEHICLE_INITIAL_SPEED, 0.0))
    add_agent!(
        spawn_pos,
        model,
        initial_vel
    )
end

# Generalize
function vehicle_poly()
    hw = P.VEHICLE_WIDTH / 2
    pt = Point2f[(0, -hw), (P.VEHICLE_LENGTH, -hw), (P.VEHICLE_LENGTH, hw), (0, hw)]
    return Polygon(pt)
end

function vehicle_marker(v::VehicleAgent)
    φ = atan(v.vel[2], v.vel[1])
    rotate2D(vehicle_poly(), φ)
end

function vehicle_step!(agent, model)
    move_agent!(agent, model)
end

function plot_environment!(model)
    custom_setup_environment!(model.env)
    drawRoad!(model.env)
end

function custom_setup_environment!(twir::TwoWayIntersectingRoads)
    # Make the signals on the edge inactive
    twir.leftRoadSeg.R2Lroad.signal.active = false
    twir.rightRoadSeg.L2Rroad.signal.active = false
    twir.topRoadSeg.B2Troad.signal.active = false
    twir.bottomRoadSeg.T2Broad.signal.active = false
    # Alternate L2R and T2B road signals
    setSignalState!(twir.leftRoadSeg.L2Rroad.signal, :green)
    setSignalState!(twir.rightRoadSeg.R2Lroad.signal, :green)
    setSignalState!(twir.topRoadSeg.T2Broad.signal, :red)
    setSignalState!(twir.bottomRoadSeg.B2Troad.signal, :red)
    foreach(draw_cross_path!, cross_paths(twir))
end

function cross_paths(twir)
    cps = []
    xpos = top_boundary(twir.leftRoadSeg)[2][1] - P.CROSS_PATH_WIDTH - P.CROSS_PATH_MARGIN
    push!(cps, HcrossPath(twir.leftRoadSeg, xpos))
    xpos = top_boundary(twir.rightRoadSeg)[1][1] + P.CROSS_PATH_WIDTH + P.CROSS_PATH_MARGIN
    push!(cps, HcrossPath(twir.rightRoadSeg, xpos))
    ypos = right_boundary(twir.topRoadSeg)[2][2] + P.CROSS_PATH_WIDTH + P.CROSS_PATH_MARGIN
    push!(cps, VcrossPath(twir.topRoadSeg, ypos))
    ypos = right_boundary(twir.bottomRoadSeg)[1][2] - P.CROSS_PATH_WIDTH - P.CROSS_PATH_MARGIN
    push!(cps, VcrossPath(twir.bottomRoadSeg, ypos))
    return cps
end

function mark_spawn_positions(twir::TwoWayIntersectingRoads)
    for sp in spawnPos(twir)
        spawnPt = Point2f(sp.pos...)
        scatter!(spawnPt, color=:green, markersize=P.SIGNAL_MS)
    end
end

function model_step!(model)
    model.tick += 1
    model.steptext[] = "Step: " * string(model.tick)
    if (model.tick % model.spawn_rate == 0)
        rnd_spawn_pos = rand(spawnPos(model.env))
        @show rnd_spawn_pos
        spawn_vel = rnd_spawn_pos.orient .* P.VEHICLE_INITIAL_SPEED
        @show spawn_vel
        add_vehicle!(rnd_spawn_pos.pos, model, spawn_vel)
    end
    draw_signal!(model.env)
    @show model.tick
end

function plot_vehicles!()
    axiskwargs = (title=P.PLOT_TITLE, titlealign=P.PLOT_TITLE_ALIGN)
    model = initialize()
    params = Dict(
        :spawn_rate => P.VSR_MIN:P.VSR_INC:P.VSR_MAX,
    )
    fig, _ = abmplot(
        model;
        agent_step! = vehicle_step!,
        model_step! = model_step!,
        am=vehicle_marker,
        params=params,
        ac=:green,
        axis=axiskwargs
    )
    text!(model.steptext, position = (1000, 1000), textsize = 10)
    plot_environment!(model)
    return fig
end

#intersectingRoads = TwoWayIntersectingRoads(2000,0,4000,2000,4000,0)

#a macro for easing development only. will be removed later
macro pv()
    quote
        plot_vehicles!()
    end
end

# end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.answer;
#
# import typedb.client.api.concept.Concept;
#
# import javax.annotation.CheckReturnValue;
#
# public interface NumericGroup {
#
#     @CheckReturnValue
#     Concept owner();
#
#     @CheckReturnValue
#     Numeric numeric();
# }


text:
########################################################################
# Mapped atoms.
########################################################################

"Mapped atom."
struct MappedAtom{T1<:Int64, T2<:AbstractLinearOp, T3<:AbstractAtom} <: AbstractAtom
    n::T1
    map::T2
    atom::T3
    function MappedAtom(map::AbstractLinearOp, atom::AbstractAtom)
        size(map, 2) == length(atom) || error("size(M) incompatible with atom")
        n = size(map, 1)
        new{Int64, AbstractLinearOp, AbstractAtom}(n, map, atom)
    end
end

Base.length(a::MappedAtom) = a.n
Base.vec(a::MappedAtom) = a.map * a.atom
Base.:(*)(M::AbstractLinearOp, a::MappedAtom) = M*vec(a)
LinearAlgebra.mul!(Ma::AbstractVector{Float64}, M::LinearOp, a::MappedAtom) = mul!(Ma, M, vec(a))

########################################################################
# Mapped atomic set.
########################################################################

"""
    MappedAtomicSet(n::Integer, A::AbstractAtomicSet, M::LinearMap)

Linear transformation another atomic set, including scaling.
"""
struct MappedAtomicSet{T1<:Int64, T2<:AbstractLinearOp, T3<:AbstractAtomicSet} <: AbstractAtomicSet
    n::T1
    map::T2
    set::T3
    function MappedAtomicSet(map::AbstractLinearOp, set::AbstractAtomicSet)
        size(map, 2) == length(set) || error("size(M) incompatible with atomic set")
        n = size(map, 1)
        return new{Int64, AbstractLinearOp, AbstractAtomicSet}(n, map, set)
    end
end

"""
    *(M, A::AbstractAtomicSet)

Apply the linear operator `M` on the atomic set `A`. Results in
another atomic set (MA) with the defining property

    MA = { Ma | a ∈ A }.
"""
function Base.:(*)(M::AbstractLinearOp, A::AbstractAtomicSet)
    return MappedAtomicSet(M, A)
end

function Base.:(*)(λ::Float64, A::AbstractAtomicSet)
    return MappedAtomicSet( LinearMaps.UniformScalingMap(λ, length(A)), A)
end

"""
    support(A::MappedAtomicSet, z::Vector)

Support value of a gauge with respect to a scaled atomic set.
"""
support(A::MappedAtomicSet, z::Vector{Float64}) = support(A.set, A.map'*z)

"""
    expose(A::MappedAtomicSet, z::Vector)

Obtain an atom in the face of the face of the scaled atomic set exposed by
the vector `z`.
"""
function expose!(A::MappedAtomicSet, z::Vector{Float64}, a::MappedAtom)
    expose!(A.set, A.map'*z, a.atom)
    return nothing
end

function expose(A::MappedAtomicSet, z::Vector{Float64})
    atom = expose(A.set, A.map'*z)
    return MappedAtom(A.map, atom)
end

Base.length(A::MappedAtomicSet) = A.n
atom_name(A::MappedAtomicSet) = "Map of $(atom_name(A.set))"
atom_description(A::MappedAtomicSet) = "M ⋅ $(atom_description(A.set))"
atom_parameters(A::MappedAtomicSet) = "n = $(length(A))"

########################################################################
# Face of the Mapped Atomic Set
########################################################################

struct MappedFace{T1<:Int64, T2<:AbstractLinearOp, T3<:AbstractFace} <: AbstractFace
    n::T1
    k::T1
    map::T2
    fa::T3
    function MappedFace(map::AbstractLinearOp, fa::AbstractFace)
        size(map, 2) == length(fa) || error("size(M) incompatible with face")
        n = size(map, 1)
        k = rank(fa)
        new{Int64, AbstractLinearOp, AbstractFace}(n, k, map, fa)
    end
end

"""
    face(A, z)

Return a face of the atomic set `A` exposed by the vector `z`.
"""
function face!(A::MappedAtomicSet, z::Vector{Float64}, F::MappedFace)
    face!(A.set, A.map'*z, F.fa)
    return nothing
end

function face(A::MappedAtomicSet, z::Vector{Float64})
    fa = face(A.set, A.map'*z)
    return MappedFace(A.map, fa)
end

"""
Given a face and a set of weights, reveal the corresponding
point on the face.
"""
Base.:(*)(F::MappedFace, c::Vector) = (F.map*F.fa)*c
Base.:(*)(M::AbstractLinearOp, F::MappedFace) = M*F.map*F.fa
Base.:(*)(λ::Real, F::MappedFace) = λ*F.map*F.fa
Base.length(F::MappedFace) = F.n
rank(F::MappedFace) = F.k
vec(F::MappedFace) = F.map*F.fa
face_name(F::MappedFace) = "Mapped Atomic Set"
face_parameters(F::MappedFace) = "rank = $(rank(F)); n = $(length(F))"


text:
using Test

import HiGHS

function test_model()
    cc = [1.0, -2.0]
    cl = [0.0, 0.0]
    cu = [10.0, 10.0]
    ru = [2.0, 1.0]
    rl = [0.0, 0.0]
    astart = [0, 2, 4]
    aindex = [0, 1, 0, 1]
    avalue = [1.0, 2.0, 1.0, 3.0]
    return (cc, cl, cu, rl, ru, astart, aindex, avalue)
end

@testset "Direct C call" begin
    (colcost, collower, colupper, rowlower, rowupper, astart, aindex, avalue) = test_model()
    n_col = convert(Cint, size(colcost, 1))
    n_row = convert(Cint, size(rowlower, 1))
    n_nz = convert(Cint, size(aindex, 1))

    colcost = convert(Array{Cdouble}, colcost)
    collower = convert(Array{Cdouble}, collower)
    colupper = convert(Array{Cdouble}, colupper)

    rowlower = convert(Array{Cdouble}, rowlower)
    rowupper = convert(Array{Cdouble}, rowupper)
    matstart = convert(Array{Cint}, astart)
    matindex = convert(Array{Cint}, aindex)
    matvalue = convert(Array{Cdouble}, avalue)

    # solution = HighsSolution(Array{Cdouble, 1}(undef, n_col), Array{Cdouble, 1}(undef, n_col), Array{Cdouble, 1}(undef, n_row),  Array{Cdouble, 1}(undef, n_row))
    # basis = HighsBasis(Array{Cint, 1}(undef, n_col), Array{Cint, 1}(undef, n_row))
    colvalue, coldual = (Array{Cdouble, 1}(undef, n_col), Array{Cdouble, 1}(undef, n_col))
    rowvalue, rowdual = (Array{Cdouble, 1}(undef, n_row),  Array{Cdouble, 1}(undef, n_row))
    colbasisstatus, rowbasisstatus = (Array{Cint, 1}(undef, n_col), Array{Cint, 1}(undef, n_row))

    modelstatus = Ref{Cint}(42)
    status = HiGHS.Highs_call(n_col, n_row, n_nz, colcost, collower, colupper, rowlower, rowupper, matstart, matindex, matvalue, colvalue, coldual, rowvalue, rowdual, colbasisstatus, rowbasisstatus, modelstatus)
    @test status == 0
    @test modelstatus[] == 9 # optimal
end

include("MOI_wrapper.jl")