text:
# ExactReach, ai2, maxSens

using NeuralVerification, LazySets, Test, LinearAlgebra, GLPKMathProgInterface
import NeuralVerification: ReLU, Id

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

small_nnet_file = "$(@__DIR__)/../examples/networks/small_nnet.nnet"
small_nnet_id_file = "$(@__DIR__)/../examples/networks/small_nnet_id.nnet"

mnist1_file = "$(@__DIR__)/../examples/networks/mnist1.nnet"
mnist2_file = "$(@__DIR__)/../examples/networks/mnist2.nnet"
mnist3_file = "$(@__DIR__)/../examples/networks/mnist3.nnet"
mnist4_file = "$(@__DIR__)/../examples/networks/mnist4.nnet"

small_nnet = mnist1 = read_nnet(small_nnet_file, last_layer_activation = ReLU())
small_nnet_id = mnist1 = read_nnet(small_nnet_id_file, last_layer_activation = Id())

mnist1 = read_nnet(mnist1_file, last_layer_activation = Id())
mnist2 = read_nnet(mnist2_file, last_layer_activation = Id())
mnist3 = read_nnet(mnist3_file, last_layer_activation = Id())
mnist4 = read_nnet(mnist4_file, last_layer_activation = Id())

println("###### Problem type - input:HPolytope, output:HPolytope ######")
println("###### Network: small_nnet, problem: holds              ######")

in_hyper  = Hyperrectangle(low = [-0.9], high = [0.9])
in_hpoly  = convert(HPolytope, in_hyper)
#out_superset    = Hyperrectangle(low = [30.0], high = [80.0])

problem_holds    = Problem(small_nnet, in_hyper, HPolytope([HalfSpace([1.], 100.)]))

solver = FastLin(10, 10.0, 1.0)
println("$(typeof(solver)) - small_nnet")
timed_result =@timed solve(solver, problem_holds)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

solver = FastLip(10, 10.0, 1.0)
println("$(typeof(solver)) - small_nnet")
timed_result =@timed solve(solver, problem_holds)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = MIPVerify(optimizer)
println("$(typeof(solver)) - small_nnet")
timed_result =@timed solve(solver, problem_holds)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = ILP(optimizer, 1)
println("$(typeof(solver)) - small_nnet")
timed_result =@timed solve(solver, problem_holds)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

println("###### Problem type - input:HPolytope, output:HPolytope ######")

# entry 23 in MNIST datset
input_center = [0.0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [-1311.1257826380004,4633.767704436501,-654.0718535670002,-1325.349417307,1175.2361184373997,-1897.8607293569007,-470.3405972940001,830.8337987382,-377.7467076115001,572.3674015264198]

in_epsilon = 1.0 #0-255
out_epsilon = 10.0 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)

A = Matrix(undef, 2, 1)
A = [1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ,0.0]'
b = [0.0]
outputSet = HPolytope(A, b)

problem_mnist1 = Problem(mnist1, inputSet, outputSet)
## MNIST1 ##
print("\n\n\n")
println("###### Network: mnist1                                  ######")
print("\n\n\n")

solver = FastLin(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist1")
timed_result =@timed solve(solver, problem_mnist1)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

solver = FastLip(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist1")
timed_result =@timed solve(solver, problem_mnist1)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = MIPVerify(optimizer)
println("$(typeof(solver)) - mnist1")
timed_result =@timed solve(solver, problem_mnist1)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = ILP(optimizer, 1)
println("$(typeof(solver)) - mnist1")
timed_result =@timed solve(solver, problem_mnist1)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

problem_mnist2 = Problem(mnist2, inputSet, outputSet)
## MNIST1 ##
print("\n\n\n")
println("###### Network: mnist2                                  ######")
print("\n\n\n")

solver = FastLin(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist2")
timed_result =@timed solve(solver, problem_mnist2)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

solver = FastLip(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist2")
timed_result =@timed solve(solver, problem_mnist2)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = MIPVerify(optimizer)
println("$(typeof(solver)) - mnist2")
timed_result =@timed solve(solver, problem_mnist2)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = ILP(optimizer, 1)
println("$(typeof(solver)) - mnist2")
timed_result =@timed solve(solver, problem_mnist2)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

problem_mnist3 = Problem(mnist3, inputSet, outputSet)
## MNIST1 ##
print("\n\n\n")
println("###### Network: mnist3                                  ######")
print("\n\n\n")

solver = FastLin(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist3")
timed_result =@timed solve(solver, problem_mnist3)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

solver = FastLip(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist3")
timed_result =@timed solve(solver, problem_mnist3)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = MIPVerify(optimizer)
println("$(typeof(solver)) - mnist3")
timed_result =@timed solve(solver, problem_mnist3)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

optimizer = GLPKSolverMIP()
solver = ILP(optimizer, 1)
println("$(typeof(solver)) - mnist3")
timed_result =@timed solve(solver, problem_mnist3)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

problem_mnist4 = Problem(mnist4, inputSet, outputSet)
## MNIST1 ##
print("\n\n\n")
println("###### Network: mnist4                                  ######")
print("\n\n\n")

solver = FastLin(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist4")
timed_result =@timed solve(solver, problem_mnist4)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

solver = FastLip(10, 10.0, 1.0)
println("$(typeof(solver)) - mnist4")
timed_result =@timed solve(solver, problem_mnist4)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

#=
optimizer = GLPKSolverMIP()
solver = MIPVerify(optimizer)
println("$(typeof(solver)) - mnist4")
timed_result =@timed solve(solver, problem_mnist4)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")
=#
optimizer = GLPKSolverMIP()
solver = ILP(optimizer, 1)
println("$(typeof(solver)) - mnist4")
timed_result =@timed solve(solver, problem_mnist4)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")

# ACAS
print("\n\n\n")
println("###### Problem type - input:HPolytope, output:HPolytope ######")
print("\n\n\n")
println("###### Network: acas                                    ######")
print("\n\n\n")

#acas_file = "$(@__DIR__)/../examples/networks/ACASXU_run2a_1_1_tiny_4.nnet"
acas_file = "$(@__DIR__)/../examples/networks/ACASXU_run2a_4_5_batch_2000.nnet"
#ACASXU_run2a_4_5_batch_2000.nnet

acas_nnet = read_nnet(acas_file, last_layer_activation = Id())

# ACAS PROPERTY 10 - modified
# Original input range:
# LOWER BOUND: array([[ 0.21466922,  0.11140846, -0.4999999 ,  0.3920202 ,  0.15      ]])
# UPPER BOUND: array([[ 0.58819589,  0.4999999 , -0.49840835,  0.66474747,  0.65      ]])
A0 = Matrix{Float64}(I, 5, 5)
A1 = -Matrix{Float64}(I, 5, 5)
A = vcat(A0, A1)

#b_lower = [0.21466922, 0.11140846, -0.4999999, 0.52838384, 0.4]
#b_upper = [0.58819589, 0.4999999 , -0.4999999, 0.52838384, 0.4]

#b_lower = [ 0.21466922,  0.11140846, -0.4999999 ,  0.3920202 ,  0.15      ]
#b_upper = [ 0.58819589,  0.4999999 , -0.49840835,  0.66474747,  0.65      ]

#b_lower = [ 0.21466922,  0.11140846, -0.4999999 ,  0.3920202 ,  0.20      ]
#b_upper = [ 0.58819589,  0.4999999 , -0.49840835,  0.66474747,  0.65      ]

b_lower = [ 0.21466922,  0.11140846, -0.4999999 ,  0.3920202 ,  0.4      ]
b_upper = [ 0.58819589,  0.4999999 , -0.49840835,  0.66474747,  0.4      ]

in_hyper  = Hyperrectangle(low = b_lower, high = b_upper)
inputSet = convert(HPolytope, in_hyper)

#inputSet = HPolytope(A, b)

#A = Matrix(undef, 2, 1)
#A = [1.0, 0.0, 0.0, 0.0, -1.0]'
#b = [0.0]
#outputSet = HPolytope(A, b)

outputSet = HPolytope([HalfSpace([1.0, 0.0, 0.0, 0.0, -1.0], 0.0)])

problem_polytope_polytope_acas = Problem(acas_nnet, in_hyper, outputSet)

print("\n\n\n")
println("###### Network: acas                                  ######")
print("\n\n\n")

solver = FastLin(10, 10.0, 1.0)
println("$(typeof(solver)) - acas")
timed_result =@timed solve(solver, problem_polytope_polytope_acas)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

solver = FastLip(10, 10.0, 1.0)
println("$(typeof(solver)) - acas")
timed_result =@timed solve(solver, problem_polytope_polytope_acas)
print(" - Time: " * string(timed_result[2]) * " s")
print(" - Output: ")
print(timed_result[1])
println("")

#=
optimizer = GLPKSolverMIP()
solver = MIPVerify(optimizer)
println("$(typeof(solver)) - acas")
timed_result =@timed solve(solver, problem_polytope_polytope_acas)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")
=#
optimizer = GLPKSolverMIP()
solver = ILP(optimizer, 1)
println("$(typeof(solver)) - acas")
timed_result =@timed solve(solver, problem_polytope_polytope_acas)
println(" - Time: " * string(timed_result[2]) * " s")
println(" - Output: ")
println(timed_result[1])
println("")


text:
using LinearAlgebra, Optim
include("GenerateExamples.jl")
include("DampedSystem.jl")
include("CauchyLike.jl")
include("ComplexSymDPR1.jl")

# Generate Small example

ext: jl


text:
precompile(Tuple{typeof(Base.similar), Array{Base.Grisu.Bignums.Bignum, 1}})
precompile(Tuple{typeof(Base.length), Array{Base.Grisu.Bignums.Bignum, 1}})
precompile(Tuple{typeof(Base.deepcopy_internal), Any, Base.IdDict{Any, Any}})
precompile(Tuple{typeof(Base.deepcopy_internal), Array{UInt32, 1}, Base.IdDict{Any, Any}})
precompile(Tuple{typeof(MacroTools.__init__)})
precompile(Tuple{typeof(Random.shuffle!), Random.MersenneTwister, Array{Symbol, 1}})
precompile(Tuple{typeof(Zlib_jll.__init__)})
precompile(Tuple{Core.var"#Type##kw", NamedTuple{(:libc, :compiler_abi), Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.FreeBSD}, Symbol})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.FreeBSD, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.FreeBSD})
precompile(Tuple{typeof(Libdl.dlopen), String, UInt32})
precompile(Tuple{typeof(JLLWrappers.get_julia_libpaths)})
precompile(Tuple{typeof(Base.join), Array{String, 1}, Char})
precompile(Tuple{typeof(Lz4_jll.__init__)})
precompile(Tuple{typeof(Libdl.dlopen), String})
precompile(Tuple{typeof(Zstd_jll.__init__)})
precompile(Tuple{typeof(Blosc_jll.__init__)})
precompile(Tuple{Type{NamedTuple{(:libgfortran_version, :libstdcxx_version, :cxxstring_abi), T} where T<:Tuple}, Tuple{Nothing, Nothing, Symbol}})
precompile(Tuple{Core.var"#Type##kw", NamedTuple{(:libgfortran_version, :libstdcxx_version, :cxxstring_abi), Tuple{Nothing, Nothing, Symbol}}, Type{Pkg.BinaryPlatforms.CompilerABI}})
precompile(Tuple{typeof(Base.foreach), Blosc_jll.var"#7#9", Tuple{Array{String, 1}, Array{String, 1}, Array{String, 1}}})
precompile(Tuple{typeof(Base.foreach), Blosc_jll.var"#8#10", Tuple{Array{String, 1}, Array{String, 1}, Array{String, 1}}})
precompile(Tuple{typeof(Blosc.__init__)})
precompile(Tuple{typeof(Requires.__init__)})
precompile(Tuple{typeof(Requires.loadpkg), Base.PkgId})
precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(nghttp2_jll.__init__)})
precompile(Tuple{typeof(OpenSSL_jll.__init__)})
precompile(Tuple{typeof(MbedTLS_jll.__init__)})
precompile(Tuple{typeof(LibSSH2_jll.__init__)})
precompile(Tuple{typeof(Base.foreach), LibSSH2_jll.var"#7#9", Tuple{Array{String, 1}}})
precompile(Tuple{typeof(Base.foreach), LibSSH2_jll.var"#8#10", Tuple{Array{String, 1}}})
precompile(Tuple{typeof(LibCURL_jll.__init__)})
precompile(Tuple{typeof(Base.foreach), LibCURL_jll.var"#10#12", Tuple{Array{String, 1}, Array{String, 1}, Array{String, 1}, Array{String, 1}}})
precompile(Tuple{typeof(Base.foreach), LibCURL_jll.var"#11#13", Tuple{Array{String, 1}, Array{String, 1}, Array{String, 1}, Array{String, 1}}})
precompile(Tuple{typeof(HDF5_jll.__init__)})
precompile(Tuple{typeof(HDF5.blosc_set_local), Int64, Int64, Int64})
precompile(Tuple{typeof(HDF5.blosc_filter), UInt32, UInt64, Ptr{UInt32}, UInt64, Ptr{UInt64}, Ptr{Ptr{Nothing}}})
precompile(Tuple{typeof(HDF5.__init__)})
precompile(Tuple{typeof(Base.iterate), Int64})
precompile(Tuple{typeof(HDF5.h5p_set_char_encoding), HDF5.Properties, Int64})
precompile(Tuple{typeof(HDF5.h5p_set_create_intermediate_group), HDF5.Properties, Int64})
precompile(Tuple{typeof(Base.get!), Type{Array{Function, 1}}, Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(Base.push!), Array{Function, 1}, Function})
precompile(Tuple{typeof(DocStringExtensions.__init__)})
precompile(Tuple{typeof(MPI.__init__)})
precompile(Tuple{typeof(MPI.Get_library_version)})
precompile(Tuple{typeof(Base.allocatedinline), Type{Array{Function, 1}}})
precompile(Tuple{typeof(Base._delete!), Base.Dict{Base.PkgId, Array{Function, 1}}, Int64})
precompile(Tuple{HDF5.var"#57#60"})
precompile(Tuple{HDF5.var"#58#61"})
precompile(Tuple{HDF5.var"#59#62"})
precompile(Tuple{Requires.var"#@include", LineNumberNode, Module, Any})
precompile(Tuple{typeof(Base.read), String})
precompile(Tuple{typeof(Requires._include_path), String})
precompile(Tuple{typeof(Base.Filesystem.isfile), String})
precompile(Tuple{typeof(HDF5.include), String})
precompile(Tuple{typeof(Libdl.dlopen), Function, String})
precompile(Tuple{Libdl.var"##dlopen#4", Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(Libdl.dlopen), HDF5.var"#82#83", String})
precompile(Tuple{typeof(Base.in), Int64, Tuple{Int64, Int64}})
precompile(Tuple{DocStringExtensions.var"#29#30"{typeof(DocStringExtensions.template_hook)}, LineNumberNode, Vararg{Any, N} where N})
precompile(Tuple{typeof(DocStringExtensions.template_hook), LineNumberNode, Module, String, Expr})
precompile(Tuple{typeof(libcxxwrap_julia_jll.__init__)})
precompile(Tuple{typeof(Pkg.Artifacts.do_artifact_str), String, Base.Dict{String, Any}, String, Module})
precompile(Tuple{Pkg.Artifacts.var"#ensure_artifact_installed##kw", NamedTuple{(:platform,), Tuple{Pkg.BinaryPlatforms.Linux}}, typeof(Pkg.Artifacts.ensure_artifact_installed), String, Base.Dict{String, Any}, String})
precompile(Tuple{Pkg.Artifacts.var"##ensure_artifact_installed#42", Pkg.BinaryPlatforms.Platform, Bool, Bool, typeof(Pkg.Artifacts.ensure_artifact_installed), String, Base.Dict{String, Any}, String})
precompile(Tuple{Pkg.Artifacts.var"##query_override#7", Base.Dict{Symbol, Base.Dict{K, V} where V where K}, typeof(Pkg.Artifacts.query_override), Base.SHA1})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.__init__)})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.checkversion)})
precompile(Tuple{typeof(CxxWrap.StdLib.__init__)})
precompile(Tuple{typeof(Libdl.dlsym), Ptr{Nothing}, Symbol})
precompile(Tuple{typeof(Base.iterate), Array{CxxWrap.CxxWrapCore.CppFunctionInfo, 1}})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore._register_function_pointers), CxxWrap.CxxWrapCore.CppFunctionInfo, Bool})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.methodkey), CxxWrap.CxxWrapCore.CppFunctionInfo})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore._module_name_hash), Module, UInt64})
precompile(Tuple{typeof(Base.iterate), Array{CxxWrap.CxxWrapCore.CppFunctionInfo, 1}, Int64})
precompile(Tuple{typeof(Base.Iterators.enumerate), Array{Tuple{Symbol, Symbol, Symbol, UInt64}, 1}})
precompile(Tuple{typeof(Base.iterate), Base.Iterators.Enumerate{Array{Tuple{Symbol, Symbol, Symbol, UInt64}, 1}}, Tuple{Int64}})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore._get_function_pointer), Tuple{Symbol, Symbol, Symbol, UInt64}})
precompile(Tuple{typeof(Base.setindex!), Array{Tuple{Ptr{Nothing}, Ptr{Nothing}, Bool}, 1}, Tuple{Ptr{Nothing}, Ptr{Nothing}, Bool}, Int64})
precompile(Tuple{typeof(Base.iterate), Base.Iterators.Enumerate{Array{Tuple{Symbol, Symbol, Symbol, UInt64}, 1}}, Tuple{Int64, Int64}})
precompile(Tuple{typeof(Parsers.__init__)})
precompile(Tuple{typeof(Base.Threads.resize_nthreads!), Array{Base.MPFR.BigFloat, 1}, Base.MPFR.BigFloat})
precompile(Tuple{typeof(JuliaChem.JERI.__init__)})
precompile(Tuple{typeof(Base.typeintersect), Any, Any})
precompile(Tuple{typeof(Base.allocatedinline), Type{JuliaChem.JCModules.Atom}})
precompile(Tuple{typeof(Base.allocatedinline), Type{JuliaChem.JCModules.Shell}})
precompile(Tuple{typeof(Base.allocatedinline), Type{Array{Float64, 1}}})
precompile(Tuple{Base.Cartesian.var"#@nloops", LineNumberNode, Module, Any, Any, Any, Vararg{Any, N} where N})
precompile(Tuple{typeof(Base.Cartesian._nloops), Int64, Symbol, Expr, Expr})
precompile(Tuple{typeof(Base.Cartesian.lreplace), Expr, Symbol, Int64})
precompile(Tuple{Base.var"##s91#155", Any, Any, Any, Any, Any})
precompile(Tuple{typeof(Base.merge_types), Tuple{Symbol}, Type{NamedTuple{(:swmr,), Tuple{Bool}}}, Type{NamedTuple{(), Tuple{}}}})
precompile(Tuple{Type{StaticArrays.Size{S} where S}, Type{Tuple{3}}})
precompile(Tuple{typeof(StaticArrays._Length), Int64})
precompile(Tuple{StaticArrays.var"##s173#184", Any, Any, Any, Any})
precompile(Tuple{Type{Base.Generator{I, F} where F where I}, StaticArrays.var"#185#186", Base.UnitRange{Int64}})
precompile(Tuple{typeof(Base.collect), Base.Generator{Base.UnitRange{Int64}, StaticArrays.var"#185#186"}})
precompile(Tuple{StaticArrays.var"##s37#6", Any, Any, Any, Any, Any, Type{T} where T, Type{T} where T, Type{T} where T, Any})
precompile(Tuple{typeof(Base.all), Function, Core.SimpleVector})
precompile(Tuple{typeof(Base._all), StaticArrays.var"#7#8", Core.SimpleVector, Base.Colon})
precompile(Tuple{typeof(StaticArrays.tuple_prod), Type{Tuple{3}}})
precompile(Tuple{typeof(Base.:(*)), Int64})
precompile(Tuple{typeof(StaticArrays.tuple_minimum), Type{Tuple{3}}})
precompile(Tuple{typeof(Base.minimum), Tuple{Int64}})
precompile(Tuple{typeof(StaticArrays.tuple_length), Type{Tuple{3}}})
precompile(Tuple{Base.Cartesian.var"#@nexprs", LineNumberNode, Module, Int64, Expr})
precompile(Tuple{Base.Cartesian.var"#@ncall", LineNumberNode, Module, Int64, Any, Vararg{Any, N} where N})
precompile(Tuple{typeof(Base._compute_eltype), Type{Tuple{Base.StackTraces.StackFrame, Int64}}})
precompile(Tuple{Base.var"##s828#449", Any, Any, Any, Any, Any, Any})
precompile(Tuple{Base.Cartesian.var"#@ntuple", LineNumberNode, Module, Int64, Any})
precompile(Tuple{typeof(Base.Cartesian._nloops), Int64, Symbol, Symbol, Expr, Vararg{Expr, N} where N})
precompile(Tuple{typeof(Base.Cartesian._nloops), Int64, Symbol, Expr, Expr, Vararg{Expr, N} where N})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{Type{Base.Val{x} where x}, Symbol})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{typeof(Base.allocatedinline), Type{Union{Tuple{Int64, Int64, Int64, Int64}, Symbol}}})
precompile(Tuple{Type{StaticArrays.Size{S} where S}, Int64})
precompile(Tuple{StaticArrays.var"##s37#12", Any, Any, Any, Any})
precompile(Tuple{typeof(StaticArrays.tuple_length), Type{Tuple{1}}})
precompile(Tuple{typeof(StaticArrays.tuple_prod), Type{Tuple{1}}})
precompile(Tuple{typeof(StaticArrays.tuple_minimum), Type{Tuple{1}}})
precompile(Tuple{StaticArrays.var"##s37#2", Any, Any, Any, Any, Any})
precompile(Tuple{typeof(Base.ntuple), StaticArrays.var"#3#4", Base.Val{1}})
precompile(Tuple{Type{StaticArrays.Size{S} where S}, Type{Tuple{12}}})
precompile(Tuple{typeof(StaticArrays.tuple_prod), Type{Tuple{12}}})
precompile(Tuple{typeof(StaticArrays.tuple_minimum), Type{Tuple{12}}})
precompile(Tuple{typeof(StaticArrays.tuple_length), Type{Tuple{12}}})
precompile(Tuple{Type{StaticArrays.Size{S} where S}, Type{Tuple{24}}})
precompile(Tuple{typeof(StaticArrays.tuple_prod), Type{Tuple{24}}})
precompile(Tuple{typeof(StaticArrays.tuple_minimum), Type{Tuple{24}}})
precompile(Tuple{typeof(StaticArrays.tuple_length), Type{Tuple{24}}})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{typeof(Base.allocatedinline), Type{MPI.Status}})
precompile(Tuple{typeof(Base.allocatedinline), Type{Array{Float64, 2}}})
precompile(Tuple{typeof(Base._compute_eltype), Type{Tuple{Base.Pair{String, Bool}, Base.Pair{String, Base.Dict{String, String}}}}})
precompile(Tuple{typeof(Base._compute_eltype), Type{Tuple{Base.Pair{String, Float64}, Base.Pair{String, Bool}, Base.Pair{String, Base.Dict{String, Float64}}}}})
precompile(Tuple{typeof(Base._compute_eltype), Type{Tuple{Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Float64}, Base.Pair{String, Bool}}}})
precompile(Tuple{Type{Base.Broadcast.BroadcastStyle}, Base.Broadcast.DefaultArrayStyle{0}, StaticArrays.StaticArrayStyle{1}})
precompile(Tuple{typeof(Base.typename), DataType})
precompile(Tuple{Type{Base.Broadcast.BroadcastStyle}, Base.Broadcast.DefaultArrayStyle{1}, StaticArrays.StaticArrayStyle{1}})
precompile(Tuple{typeof(Main.full_rhf), String})
precompile(Tuple{Distributed.var"#133#134"})
precompile(Tuple{MPI.var"#19#20"{Int32, MPI.Info}})
precompile(Tuple{MPI.var"#22#23"{Int32, MPI.Comm}})
precompile(Tuple{MPI.var"#32#33"{Int32, MPI.Datatype}})
precompile(Tuple{MPI.var"#59#60"{Int32, MPI.Op}})
precompile(Tuple{MPI.var"#69#70"{Int32, MPI.Request}})
precompile(Tuple{MPI.var"#83#84"{Int32, MPI.Win}})
precompile(Tuple{MPI.var"#88#89"{Int32, MPI.FileHandle}})
precompile(Tuple{typeof(Base.println), Base.TTY, String, Vararg{Any, N} where N})
precompile(Tuple{typeof(Base.print), Base.TTY, String, Int64, Vararg{Any, N} where N})
precompile(Tuple{typeof(Base.print), Base.TTY, Int64})
precompile(Tuple{JSON.Parser.var"##parse#1", Type{T} where T, Type{Int64}, Bool, Nothing, typeof(JSON.Parser.parse), String})
precompile(Tuple{typeof(JSON.Parser.parse_value), JSON.Parser.ParserContext{Base.Dict{String, Any}, Int64, true, nothing}, JSON.Parser.MemoryParserState})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Int64, String})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Float64, String})
precompile(Tuple{typeof(JSON.Parser.skip!), JSON.Parser.MemoryParserState, UInt8, UInt8, Vararg{UInt8, N} where N})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, String, Array{Any, 1}})
precompile(Tuple{Type{Base.Dict{K, V} where V where K}, Base.Pair{String, Array{Any, 1}}})
precompile(Tuple{typeof(Base.merge!), Base.Dict{Any, Any}, Base.Dict{String, Array{Any, 1}}})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, String, Int64})
precompile(Tuple{Type{Base.Dict{K, V} where V where K}, Base.Pair{String, Int64}})
precompile(Tuple{typeof(Base.merge!), Base.Dict{Any, Any}, Base.Dict{String, Int64}})
precompile(Tuple{typeof(Base.merge!), Base.Dict{Any, Any}, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.convert), Type{Any}, String})
precompile(Tuple{typeof(Base.convert), Type{Any}, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Any, Any}, Base.Dict{String, Any}, String})
precompile(Tuple{JuliaChem.JCBasis.var"#run##kw", NamedTuple{(:output,), Tuple{String}}, typeof(JuliaChem.JCBasis.run), Base.Dict{Any, Any}, Base.Dict{Any, Any}})
precompile(Tuple{typeof(Base.convert), Type{Array{Float64, 1}}, Array{Any, 1}})
precompile(Tuple{typeof(Base.setindex!), Array{Float64, 1}, Float64, Int64})
precompile(Tuple{HDF5.var"##h5open#10", Bool, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(HDF5.h5open), JuliaChem.JCBasis.var"#4#7"{String, MPI.Comm, Array{String, 1}, String, Base.Dict{String, Int64}, Base.Dict{String, Int64}, JuliaChem.JCModules.Molecule, Array{JuliaChem.JCModules.Shell, 1}, CxxWrap.StdLib.StdVectorAllocated{CxxWrap.StdLib.StdVector{JuliaChem.JERI.Shell}}, Array{Bool, 1}, Array{Array{Float64, 1}, 1}}, String, Vararg{String, N} where N})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.Shell}})
precompile(Tuple{HDF5.var"#h5open##kw", NamedTuple{(:swmr,), Tuple{Bool}}, typeof(HDF5.h5open), String, String})
precompile(Tuple{typeof(HDF5.h5p_set_fclose_degree), HDF5.Properties, Int64})
precompile(Tuple{typeof(Base.push!), CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.Atom}, JuliaChem.JERI.AtomAllocated})
precompile(Tuple{typeof(Base.read), HDF5.Group})
precompile(Tuple{typeof(HDF5.get_mem_compatible_jl_type), HDF5.Datatype})
precompile(Tuple{typeof(Base.read), HDF5.Dataset, Type{Float64}})
precompile(Tuple{Type{Array{Float64, N} where N}, UndefInitializer, UInt64})
precompile(Tuple{typeof(HDF5.dataspace), Array{Float64, 1}})
precompile(Tuple{typeof(HDF5.h5d_read), HDF5.Dataset, HDF5.Datatype, HDF5.Dataspace, HDF5.Dataspace, HDF5.Properties, Array{Float64, 1}})
precompile(Tuple{typeof(Base.close), HDF5.Dataspace})
precompile(Tuple{typeof(Base.read), HDF5.Dataset, Type{HDF5.FixedString{1, 0}}})
precompile(Tuple{Type{Array{HDF5.FixedString{1, 0}, N} where N}, UndefInitializer, Int64})
precompile(Tuple{typeof(HDF5.dataspace), Array{HDF5.FixedString{1, 0}, 1}})
precompile(Tuple{typeof(HDF5.h5d_read), HDF5.Dataset, HDF5.Datatype, HDF5.Dataspace, HDF5.Dataspace, HDF5.Properties, Array{HDF5.FixedString{1, 0}, 1}})
precompile(Tuple{typeof(Base.Broadcast.broadcasted), Function, Array{HDF5.FixedString{1, 0}, 1}})
precompile(Tuple{Type{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Axes, F, Args} where Args<:Tuple where F where Axes}, typeof(HDF5.normalize_types), Tuple{Array{HDF5.FixedString{1, 0}, 1}}})
precompile(Tuple{typeof(Base.Broadcast.instantiate), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(HDF5.normalize_types), Tuple{Array{HDF5.FixedString{1, 0}, 1}}}})
precompile(Tuple{typeof(Base.copy), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(HDF5.normalize_types), Tuple{Array{HDF5.FixedString{1, 0}, 1}}}})
precompile(Tuple{typeof(Base.close), HDF5.Datatype})
precompile(Tuple{typeof(Base.close), HDF5.Properties})
precompile(Tuple{typeof(Base.close), HDF5.Dataset})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), JuliaChem.JERI.AtomAllocated})
precompile(Tuple{typeof(Base.similar), Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(HDF5.normalize_types), Tuple{Base.Broadcast.Extruded{Array{HDF5.FixedString{1, 0}, 1}, Tuple{Bool}, Tuple{Int64}}}}, Type{String}})
precompile(Tuple{typeof(Base.Broadcast.copyto_nonleaf!), Array{String, 1}, Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(HDF5.normalize_types), Tuple{Base.Broadcast.Extruded{Array{HDF5.FixedString{1, 0}, 1}, Tuple{Bool}, Tuple{Int64}}}}, Base.OneTo{Int64}, Int64, Int64})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Array{Float64, 1}, String})
precompile(Tuple{Type{Array{Float64, N} where N}, UndefInitializer, UInt64, UInt64})
precompile(Tuple{typeof(HDF5.dataspace), Array{Float64, 2}})
precompile(Tuple{typeof(HDF5.h5d_read), HDF5.Dataset, HDF5.Datatype, HDF5.Dataspace, HDF5.Dataspace, HDF5.Properties, Array{Float64, 2}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Array{Float64, 2}, String})
precompile(Tuple{typeof(Base.getindex), Base.Dict{String, Int64}, String})
precompile(Tuple{typeof(Base.convert), Type{Array{Float64, 1}}, Array{Float64, 1}})
precompile(Tuple{typeof(Base.convert), Type{Array{Float64, N} where N}, Array{Float64, 1}})
precompile(Tuple{Type{NamedTuple{(:formatters,), T} where T<:Tuple}, Tuple{PrettyTables.var"#31#33"{Array{String, 1}, Array{Int64, 1}}}})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{PrettyTables.var"#pretty_table##kw", NamedTuple{(:formatters,), Tuple{PrettyTables.var"#31#33"{Array{String, 1}, Array{Int64, 1}}}}, typeof(PrettyTables.pretty_table), Array{Float64, 2}, Array{String, 2}})
precompile(Tuple{Type{Base.IOContext{IO_t} where IO_t<:IO}, Base.TTY, Base.Pair{Symbol, Bool}})
precompile(Tuple{PrettyTables.var"#_pretty_table##kw", NamedTuple{(:formatters,), Tuple{PrettyTables.var"#31#33"{Array{String, 1}, Array{Int64, 1}}}}, typeof(PrettyTables._pretty_table), Base.IOContext{Base.TTY}, Array{Float64, 2}, Array{String, 2}})
precompile(Tuple{typeof(Base.close), HDF5.Group})
precompile(Tuple{typeof(Base._array_for), Type{Symbol}, Base.UnitRange{Int64}, Base.HasShape{1}})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Symbol, 1}, Symbol, Base.Generator{Base.UnitRange{Int64}, PrettyTables.var"#59#63"}, Int64})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Symbol, 1}, Symbol, Base.Generator{Base.UnitRange{Int64}, PrettyTables.var"#60#64"}, Int64})
precompile(Tuple{Type{PrettyTables.PrintInfo}, Array{Float64, 2}, Array{String, 2}, Array{Int64, 1}, Array{Int64, 1}, Int64, Int64, Int64, Int64, Int64, Int64, Bool, String, Bool, Nothing, Symbol, String, Array{Symbol, 1}, Tuple{}, Tuple{PrettyTables.var"#31#33"{Array{String, 1}, Array{Int64, 1}}}, Bool, String, Symbol, Array{Symbol, 1}, Tuple{}, Bool, Base.Val{:print}})
precompile(Tuple{Type{Base.IOContext{IO_t} where IO_t<:IO}, Base.GenericIOBuffer{Array{UInt8, 1}}, Base.Pair{Symbol, Bool}})
precompile(Tuple{Type{NamedTuple{(:has_color,), T} where T<:Tuple}, Tuple{Bool}})
precompile(Tuple{Core.var"#Type##kw", NamedTuple{(:has_color,), Tuple{Bool}}, Type{PrettyTables.Display}})
precompile(Tuple{typeof(Base._array_for), Type{Crayons.Crayon}, Base.UnitRange{Int64}, Base.HasShape{1}})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Crayons.Crayon, 1}, Crayons.Crayon, Base.Generator{Base.UnitRange{Int64}, PrettyTables.var"#97#101"}, Int64})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Crayons.Crayon, 1}, Crayons.Crayon, Base.Generator{Base.UnitRange{Int64}, PrettyTables.var"#98#102"}, Int64})
precompile(Tuple{typeof(Base._compute_eltype), Type{Tuple{Bool, Bool, Int64, Bool, Bool, Bool, Base.Val{:print}}}})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{typeof(PrettyTables._fill_matrix_data!), Array{String, 2}, Array{Int64, 2}, Array{Array{String, 1}, 2}, Array{Array{Int64, 1}, 2}, Array{Int64, 1}, Array{Int64, 1}, Array{Int64, 1}, Array{Int64, 1}, Array{Int64, 1}, Array{Int64, 1}, Int64, Array{Int64, 1}, Array{Float64, 2}, Array{String, 2}, Tuple{PrettyTables.var"#31#33"{Array{String, 1}, Array{Int64, 1}}}, PrettyTables.Display, Bool, Bool, Bool, Bool, Array{Bool, 1}, Bool, Array{Int64, 1}, Bool, Base.Val{:print}, Symbol})
precompile(Tuple{Base.var"##sprint#355", Nothing, Int64, typeof(Base.sprint), Function, String, Vararg{String, N} where N})
precompile(Tuple{typeof(Base.escape_string), Base.GenericIOBuffer{Array{UInt8, 1}}, String, String})
precompile(Tuple{PrettyTables.var"#31#33"{Array{String, 1}, Array{Int64, 1}}, Float64, Int64, Int64})
precompile(Tuple{Type{NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :linebreaks, :renderer), T} where T<:Tuple}, Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Base.Val{:print}}})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{PrettyTables.var"#_parse_cell_text##kw", NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Base.Val{:print}}}, typeof(PrettyTables._parse_cell_text), Float64})
precompile(Tuple{Base.var"##sprint#355", Base.Pair{Symbol, Bool}, Int64, typeof(Base.sprint), Function, Float64})
precompile(Tuple{typeof(Base.print), Base.IOContext{Base.GenericIOBuffer{Array{UInt8, 1}}}, Float64})
precompile(Tuple{typeof(Base.indexed_iterate), Tuple{Array{String, 1}, Array{Int64, 1}, Int64}, Int64})
precompile(Tuple{typeof(Base.indexed_iterate), Tuple{Array{String, 1}, Array{Int64, 1}, Int64}, Int64, Int64})
precompile(Tuple{typeof(Base.convert), Type{Array{Int64, 1}}, Array{Int64, 1}})
precompile(Tuple{typeof(Base.haskey), Base.Dict{String, Function}, String})
precompile(Tuple{Printf.var"#@sprintf", LineNumberNode, Module, Vararg{Any, N} where N})
precompile(Tuple{typeof(Printf._printf), String, Expr, String, Tuple{Symbol}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Function}, Function, String})
precompile(Tuple{Formatting.var"#1#4", Float64})
precompile(Tuple{Formatting.var"#11#12", Float64})
precompile(Tuple{PrettyTables.var"#_parse_cell_text##kw", NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Base.Val{:print}}}, typeof(PrettyTables._parse_cell_text), String})
precompile(Tuple{typeof(Base.getindex), Base.Dict{String, Function}, String})
precompile(Tuple{typeof(Base.push!), Base.Set{Any}, Int64})
precompile(Tuple{typeof(Base.in), Int64, Base.Set{Any}})
precompile(Tuple{typeof(Base.isequal), Int64, Int64})
precompile(Tuple{Type{Array{Int64, 1}}, Array{Any, 1}})
precompile(Tuple{typeof(Base.sum), Array{Int64, 1}})
precompile(Tuple{typeof(Base.println), Base.IOContext{Base.GenericIOBuffer{Array{UInt8, 1}}}, Base.SubString{String}})
precompile(Tuple{typeof(PrettyTables._print_table_header!), Base.IOContext{Base.GenericIOBuffer{Array{UInt8, 1}}}, PrettyTables.Display, Array{String, 2}, Array{String, 2}, Array{Int64, 2}, Array{Int64, 1}, Array{Int64, 1}, Int64, Int64, Array{Int64, 1}, Array{Int64, 1}, Array{Symbol, 1}, Array{Symbol, 1}, Tuple{}, Bool, Bool, PrettyTables.TextFormat, Crayons.Crayon, Array{Crayons.Crayon, 1}, Array{Crayons.Crayon, 1}, Crayons.Crayon, Crayons.Crayon})
precompile(Tuple{typeof(Base.print), Base.GenericIOBuffer{Array{UInt8, 1}}, Crayons.Crayon, String, Vararg{Any, N} where N})
precompile(Tuple{typeof(Base.print), Base.GenericIOBuffer{Array{UInt8, 1}}, Crayons.Crayon})
precompile(Tuple{typeof(PrettyTables._print_table_data), Base.IOContext{Base.GenericIOBuffer{Array{UInt8, 1}}}, PrettyTables.Display, Array{Float64, 2}, Array{Array{String, 1}, 2}, Array{Array{Int64, 1}, 2}, Array{Int64, 1}, Array{Int64, 1}, Int64, Array{Int64, 1}, Array{Int64, 1}, Array{Union{Tuple{Int64, Int64, Int64, Int64}, Symbol}, 1}, Array{Union{Int64, Symbol}, 1}, Array{Symbol, 1}, Tuple{Char, Char, Char, Char}, Tuple{}, Symbol, Int64, Tuple{}, Bool, Bool, Bool, PrettyTables.TextFormat, Crayons.Crayon, Crayons.Crayon, Crayons.Crayon})
precompile(Tuple{typeof(Base.reshape), Array{Float64, 1}, Tuple{Int64}})
precompile(Tuple{typeof(Base.copy), Array{Float64, 1}})
precompile(Tuple{Type{JuliaChem.JCModules.Shell}, Int64, Int64, Int64, Array{Float64, 1}, Array{Float64, 1}, Array{Float64, 1}, Int64, Int64, Int64, Bool})
precompile(Tuple{typeof(Base.convert), Type{Array{Float64, N} where N}, Array{Float64, 2}})
precompile(Tuple{typeof(Base.getindex), Array{Float64, 2}, Function, Int64})
precompile(Tuple{Type{JuliaChem.JERI.BasisSet}, CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.Atom}, CxxWrap.StdLib.StdVectorAllocated{CxxWrap.StdLib.StdVector{JuliaChem.JERI.Shell}}})
precompile(Tuple{typeof(JuliaChem.JERI.precompute_shell_pair_data), CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.ShellPair}, JuliaChem.JERI.BasisSetAllocated})
precompile(Tuple{typeof(Base.println), Base.TTY, Int64})
precompile(Tuple{typeof(Base.println), String, String, Vararg{Any, N} where N})
precompile(Tuple{typeof(Base.print), Base.TTY, Float64})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), JuliaChem.JERI.ShellAllocated})
precompile(Tuple{JuliaChem.JCRHF.Energy.var"#run##kw", NamedTuple{(:output,), Tuple{String}}, typeof(JuliaChem.JCRHF.Energy.run), JuliaChem.JCModules.Molecule, JuliaChem.JCModules.Basis, Base.Dict{String, Any}})
precompile(Tuple{Type{JuliaChem.JERI.OEIEngine}, CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.Atom}, JuliaChem.JERI.BasisSetAllocated, Int64})
precompile(Tuple{typeof(Printf.print_fixed), Base.TTY, Int64, Int32, Int32, Bool, Array{UInt8, 1}})
precompile(Tuple{Type{JuliaChem.JERI.TEIEngine}, JuliaChem.JERI.BasisSetAllocated, CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.ShellPair}})
precompile(Tuple{JuliaChem.JCRHF.Energy.var"#10#14"{Int64, Array{Array{Float64, 2}, 1}, Array{Float64, 2}, Array{Float64, 2}, JuliaChem.JCModules.Basis, Array{Float64, 2}, Array{Float64, 2}, Array{Array{Float64, 1}, 1}, Float64, Bool, Int64, Int64, Base.Threads.Atomic{Int64}}})
precompile(Tuple{typeof(Base.task_done_hook), Task})
precompile(Tuple{typeof(Base.unsafe_write), Base.TTY, Ptr{UInt8}, Int32})
precompile(Tuple{typeof(Base.empty), Base.Dict{Any, Any}, Type{String}, Type{Float64}})
precompile(Tuple{typeof(Base.grow_to!), Base.Dict{String, Float64}, Tuple{Base.Pair{String, Float64}, Base.Pair{String, Bool}, Base.Pair{String, Base.Dict{String, Float64}}}, Int64})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), JuliaChem.JERI.TEIEngineAllocated})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), JuliaChem.JERI.OEIEngineAllocated})
precompile(Tuple{typeof(Base.empty), Base.Dict{String, Float64}, Type{String}, Type{Real}})
precompile(Tuple{typeof(Base.merge!), Base.Dict{String, Real}, Base.Dict{String, Float64}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Real}, Bool, String})
precompile(Tuple{typeof(Base.grow_to!), Base.Dict{String, Real}, Tuple{Base.Pair{String, Float64}, Base.Pair{String, Bool}, Base.Pair{String, Base.Dict{String, Float64}}}, Int64})
precompile(Tuple{typeof(Base.empty), Base.Dict{String, Real}, Type{String}, Type{Any}})
precompile(Tuple{typeof(Base.merge!), Base.Dict{String, Any}, Base.Dict{String, Real}})
precompile(Tuple{typeof(Base.convert), Type{Real}, Float64})
precompile(Tuple{typeof(Base.convert), Type{Real}, Bool})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Base.Dict{String, Float64}, String})
precompile(Tuple{typeof(Base.grow_to!), Base.Dict{String, Any}, Tuple{Base.Pair{String, Float64}, Base.Pair{String, Bool}, Base.Pair{String, Base.Dict{String, Float64}}}, Int64})
precompile(Tuple{typeof(Base.convert), Type{Any}, Float64})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Any, Any}, Float64, String})
precompile(Tuple{typeof(Base.convert), Type{Any}, Bool})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Any, Any}, Bool, String})
precompile(Tuple{typeof(Base.convert), Type{Any}, Base.Dict{String, Float64}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Any, Any}, Base.Dict{String, Float64}, String})
precompile(Tuple{typeof(Base.empty), Base.Dict{Any, Any}, Type{String}, Type{Array{Float64, 2}}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Array{Float64, 2}}, Array{Float64, 2}, String})
precompile(Tuple{typeof(Base.grow_to!), Base.Dict{String, Array{Float64, 2}}, Tuple{Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Float64}, Base.Pair{String, Bool}}, Int64})
precompile(Tuple{typeof(Base.empty), Base.Dict{String, Array{Float64, 2}}, Type{String}, Type{Any}})
precompile(Tuple{typeof(Base.merge!), Base.Dict{String, Any}, Base.Dict{String, Array{Float64, 2}}})
precompile(Tuple{typeof(Base.grow_to!), Base.Dict{String, Any}, Tuple{Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Array{Float64, 2}}, Base.Pair{String, Float64}, Base.Pair{String, Bool}}, Int64})
precompile(Tuple{JuliaChem.JCProperties.var"#run##kw", NamedTuple{(:output,), Tuple{String}}, typeof(JuliaChem.JCProperties.run), JuliaChem.JCModules.Molecule, JuliaChem.JCModules.Basis, Base.Dict{String, Any}, Base.Dict{String, Any}})
precompile(Tuple{Type{JuliaChem.JERI.PropEngine}, CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.Atom}, JuliaChem.JERI.BasisSetAllocated})
precompile(Tuple{typeof(Base._atexit)})
precompile(Tuple{typeof(MPI._finalize)})
precompile(Tuple{Blosc.var"#1#2"})
precompile(Tuple{typeof(Base.Filesystem.temp_cleanup_purge)})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), JuliaChem.JERI.PropEngineAllocated})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.ShellPair}})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), JuliaChem.JERI.BasisSetAllocated})
precompile(Tuple{typeof(Base.close), HDF5.File})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), CxxWrap.StdLib.StdVectorAllocated{CxxWrap.StdLib.StdVector{JuliaChem.JERI.Shell}}})
precompile(Tuple{typeof(CxxWrap.CxxWrapCore.delete), CxxWrap.StdLib.StdVectorAllocated{JuliaChem.JERI.Atom}})


text:
import Base: *, sort
using Random, Statistics, DataFrames, CSV
using BioSequences, FASTX, XLSX, CPUTime
import BioSequences: iscompatible

#The following functions are the base set of tools that this program uses

# Overload the concatenation operator to combine a sequence and
# a single base, i.e. dna"AGCGTGC" * DNA_T
function *(seq::LongDNASeq, base::DNA)
    seq * LongDNASeq([base])
end

# Overload sort to order DNA codes by degeneracy:
#   A,G,C,T,M,R,W,S,Y,K,V,H,D,B,N
function sort(seq::LongSequence{DNAAlphabet{4}}; rev::Bool=false)
    sort(convert(Vector{DNA}, seq), by=degeneracy, rev=rev)
end

# Determine whether or not a candidate sequence belongs to any member of the list of sites. We only want seq that belong to sites
"""
    isvalid(seq,sites)

Check if the degeneracy of seq is less than or equal to the number of sites it occurs in.
"""
function isvalid(seq, sites)
    counter=0;
    for site in sites
        if occursin(seq,site)
            counter+=1;
        end
    end
    if degeneracy(seq;uselog2=false)<=counter # If the degeneracy of the randomer is greater than the number of sites it hits, then it must hit sites that aren't in the list
        return true
    else
    return false
    end
end

"""
    find_valid_randomer(prefix, bases, sites)

Randomly select bases to add the prefix, ensuring at least one site
appears in the randomer.
"""
function find_valid_randomer(prefix, bases, sites)
    n = length(bases)
    pre_n = length(prefix)
    randomer = prefix * (dna"-" ^ n)
    for i = 1:n
        candidates = shuffle(bases[i])
        for cand in candidates
            randomer[pre_n+i] = cand
            if isvalid(randomer[1:pre_n+i], sites)
                break
            else
                randomer[pre_n+i] = DNA_Gap
            end
        end
        if randomer[pre_n+i] == DNA_Gap
            break
        end
    end
    return randomer
end

"""
    degeneracy(base::DNA[, uselog2=true])

Calculate the degeneracy of a single base. The degeneracy
is the number of bases in the code, so
    degeneracy(DNA_N) == 4
and
    degeneracy(DNA_A) == 1

If `uselog2`, the degeneracy is log2 transformed. Because degeneracy(DNA_Gap)=0, we approximate
the log2 transformation as 2^-100.
"""
function degeneracy(base::DNA; uselog2=true)
    deg = 0.0
    if base == DNA_N
        deg = 4.0
    elseif base == DNA_A || base == DNA_C || base == DNA_G || base == DNA_T
        deg = 1.0
    elseif base == DNA_B || base == DNA_D || base == DNA_H || base == DNA_V
        deg = 3.0
    elseif (base == DNA_K || base == DNA_Y || base == DNA_S || base == DNA_W
            || base == DNA_R || base == DNA_M)
        deg = 2.0
    else # DNA_Gap
        deg = 0.0
    end

    if uselog2
        if base == DNA_Gap
            deg = 2^-100  #log2(0) is negative infinity, so we use something close to zero
        end

        return log2(deg)

    else
        return deg
    end
end

"""
    degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)

Calculate the degeneracy of a DNA sequence.
"""
function degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)
    if uselog2
        return sum(degeneracy.(seq, uselog2=true))
    else
        return prod(degeneracy.(seq, uselog2=false))
    end
end

"""
    pooled_degeneracy(pool)

Calcualte the sum degeneracy of pool.
"""
function pooled_degeneracy(pool)
    deg=0;
    for randomer in pool
        deg +=degeneracy(randomer;uselog2=false);
    end
    return deg
end

"""
    iscompatible(seq_A::LongSequence{DNAAlphabet{4}},seq_B::LongSequence{DNAAlphabet{4}})

Check if seq_A and seq_B contain at least one common sequence. Seq_A and Seq_B must be the same length.
"""
function iscompatible(seq_A::LongSequence{DNAAlphabet{4}},seq_B::LongSequence{DNAAlphabet{4}})
    if length(seq_A) == length(seq_B)
        return prod(iscompatible.(seq_A,seq_B))
    else # Sequences are different lenghts
        return false
    end
end

## Create the list of available actions. Reffered to as either ns or bases in other funcitons

"""
    define_action_space(sites)

Find the available base codes for an oligo, such that only codes that appear in sites at each postion are included.
"""
function define_action_space(sites)
    randomerlen=length(sites[1]);
    nsites=length(sites)
    bases=[];
    for i=1:randomerlen
        states=dna"-"^nsites;
        nucleotides=dna"----"
        #println("TEST: $nucleotides")
        for j=1:nsites
            states[j]=sites[j][i];
        end
        if !occursin(dna"A",states)
            nucleotides[1]=DNA_A
        else
            nucleotides[1]=DNA_Gap
        end
        if !occursin(dna"T",states)
            nucleotides[2]=DNA_T
        else
            nucleotides[2]=DNA_Gap
        end
        if !occursin(dna"C",states)
            nucleotides[3]=DNA_C
        else
            nucleotides[3]=DNA_Gap
        end
        if !occursin(dna"G",states)
            nucleotides[4]=DNA_G
        else
            nucleotides[4]=DNA_Gap
        end
        #println("TEST: $states")
        push!(bases,omit_codes(nucleotides))
        #println("Test: $nucleotides")

    end
    return bases
end

"""
    omit_codes(bases)

Helper function for define_action_space.
"""
function omit_codes(bases)
    testcodes=dna"AGCTMRWSYKVHDBN"
    testcodes_copy=deepcopy(testcodes)
    indexarray=[]
    for i= 1:length(testcodes)
        count=0;
        for nucleotide in bases
            if iscompatible(nucleotide,testcodes[i])
                count+=1;
            end
        end
        if count >0
            push!(indexarray,i)
        end
    end
    check=sort(indexarray,rev=true)
    for index in check
        deleteat!(testcodes_copy,index)
    end

    return testcodes_copy
end

#Simulator function
function simulate_random_hit_score(prefix,bases,sites; nsims=1000, horizon=length(bases), kwargs...)
    hit_scores=zeros(nsims)
    true_horizon = min(horizon, length(bases))
    Threads.@threads for i = 1:nsims   #This is where multithreading is incorporated
        hit_scores[i]= degeneracy(find_valid_randomer(prefix,bases[1:true_horizon],sites);uselog2=false)
    end
    return mean(hit_scores)
end

### MAIN FUNCTION
"""
    oligo_rollout(bases::Array{LongSequence{DNAAlphabet{4}},1},sites::Array{LongSequence{DNAAlphabet{4}},1}; simulate=simulate_random_hit_score, kwargs...)

Design an oligo that captures as many oligos from sites as possible

# Arguments
- `bases`: An array allowed bases at each postion, usually all 15 degenerate bases. Option is given because some companies restrict which degernate bases are allowed.
- `sites`: An array of non-degenerate oligos to be captured by an oligo.
- `simulate`: The choice of policy for rollout simulations. simulate_random_hit_score will perform a random rollout with reward for capturing sites.
"""
function oligo_rollout(bases::Array{LongSequence{DNAAlphabet{4}},1},sites::Array{LongSequence{DNAAlphabet{4}},1}; simulate=simulate_random_hit_score, kwargs...)
    n=length(bases)
    randomer= dna"-" ^ n
    # Find the longest oligo and use this as the horizon. (they should all be equal)
    max_len = map(length, sites) |> maximum
    horizon = max_len
    for i= 1:n
        best_base= DNA_Gap
        best_score=0;
        start = max(i - max_len, 1) # only need to look back max_len when checking the oligos
        bases=define_action_space(sites[findall(j -> iscompatible(j[start:i-1],randomer[start:i-1]),sites)])
        if length(bases[i])==1
            best_base=bases[i][1];
            #println("Shortcut Triggered at position $i")
        else
            for base in bases[i]
                if ~isvalid(randomer[start:i-1] * base, sites)
                # adding this bases creates sequence not in the oligo pool; skip it
                    continue
                end
                score= simulate(randomer[start:i-1] * base, bases[i+1:end], sites; horizon=horizon,kwargs... )

                # Update use this base if:
                #   1. the mean randomer mRNA coverage is higher than the previous best, OR
                #   2. the randomer mRNA coverage is tied with the previous best, but the new base has
                #      higher individual degeneracy
                if score > best_score || (score == best_score && degeneracy(base) > degeneracy(best_base))
                    best_score = score
                    best_base = base
                end
            end
        end
        if best_base == DNA_Gap
            # The sequence has terminated; there are no bases that can be added without making
            # a restriction site.
            break
        else
            randomer[i] = best_base
        end
    end

    return randomer
end

"""
    remove_hit_oligos(randomer,sites)

Remove oligos from sites that are captured by randomer.
"""
function remove_hit_oligos(randomer,sites)
    n=length(sites);
    hit_indicies=[]
    for i=1:n
        if occursin(sites[i],randomer)
            push!(hit_indicies,i)
        end
    end
    deleteat!(sites,hit_indicies)
    return sites
end

##Read in the oligo pool. Excel file must be just a list of randomers in column a with no header on sheet 1
function read_oligo_pool(file)
    xf=XLSX.readxlsx(file);
    x=xf[XLSX.sheetnames(xf)[1]]
    sites=[]
    n=length(x[:]);
    poolarray=x[:];
 for i in 1:n
        seq=(LongDNASeq(string(poolarray[i])))
    push!(sites,seq)
    end
    return sites
end
##Read in the oligo pool as a data frame. The oligos must have a header called "Sequence"
function read_oligo_pool(data)
    sites=[];
    for i=1:nrow(data)
        seq=(LongDNASeq(string(data[:Sequence][i])))
        push!(sites,seq)
    end
    return sites
end

#Random pool generator for testing
"""
    generate_random_pool(;randomerlen=6,randpoolsize=100,nucleotides=dna"ACGT")

Generate a random pool of oligos with length=randomerlen and pool size=randpoolsize using nucleotides
"""
function generate_random_pool(;randomerlen=6,randpoolsize=100,nucleotides=dna"ACGT")
    randpool=[]
    randomer=dna"-"^randomerlen
    for i = 1:randomerlen
        bases=shuffle(nucleotides)
        randomer[i]=bases[1];
    end
    push!(randpool,randomer)
    while length(randpool)<randpoolsize
        randomer=dna"-"^randomerlen
        for i = 1:randomerlen
            bases=shuffle(nucleotides)
            randomer[i]=bases[1];
        end
        check=true
        for j = 1:length(randpool)
            if randomer==randpool[j]
                check=false
                break
            end
        end
        if check
            push!(randpool,randomer)
        end
        if length(randpool)==length(nucleotides)^randomerlen
            break
        end
    end

    return LongDNASeq.(randpool)
end

"""
    decompress_oligo(comp_oligo)

Expand a degnerate oligo into the non-degnerate oligos it defines.
"""
function decompress_oligo(comp_oligo)
    basedict=Dict(DNA_A => dna"A", DNA_T => dna"T", DNA_C => dna"C", DNA_G => dna"G",DNA_R =>dna"AG",DNA_Y=>dna"CT",DNA_S=>dna"GC",
    DNA_W=>dna"AT",DNA_K=>dna"GT",DNA_M=>dna"AC",DNA_B=>dna"CGT",DNA_D=>dna"AGT",DNA_H=>dna"ACT",DNA_V=>dna"ACG",DNA_N=>dna"ACGT",DNA_Gap=>dna"-")
    nseqs=trunc(Int128,degeneracy(comp_oligo;uselog2=false))
    randomerlen=length(comp_oligo)
    blockarray=zeros(Int128,randomerlen+1)
    blockarray[1]=nseqs
    degsarray=zeros(Int128,randomerlen)
    repsarray=zeros(Int128,randomerlen+1)
    repsarray[1]=1
    uncomp_oligos=[dna"-"^randomerlen for _ in 1:nseqs]
    for i = 1:randomerlen
        nucleotides=basedict[comp_oligo[i]]
        degsarray[i]=length(nucleotides)
        blockarray[i+1]=trunc(Int128,blockarray[i]/degsarray[i]);
        bases=dna"-"^blockarray[i];
        repsarray[i+1]=trunc(Int128,prod(degsarray[1:i]));
            for  k in 1:degsarray[i]
                for l in 1:blockarray[i+1]
                    bases[l+blockarray[i+1]*(k-1)]=nucleotides[k]
                end
            end

        bases=bases^repsarray[i];
        for m in 1:nseqs
            uncomp_oligos[m][i]=bases[m];
        end
    end
    return uncomp_oligos
end

"""
    decompress_pool(compressed_pool)

Expand a pool of degenerate oligos into a single pool of non-degnerate oligos.
"""

function decompress_pool(compressed_pool)
    uncomppool=decompress_oligo(compressed_pool[1])
    for i = 2:length(compressed_pool)
        subset=decompress_oligo(compressed_pool[i])
        for i in subset
            push!(uncomppool,i)
        end
    end
    return uncomppool
end
function generate_poissrnd_pool(;randomerlen=6,randpoolsize=4,desiredsize=16)
    randpool=[];
    randomer=dna"-"^randomerlen;
    nucleotides=[dna"ATGC",dna"TRYKMSW",dna"BDHV",dna"N"]
    lambda=(desiredsize/randpoolsize)^(1/randomerlen);
    for i = 1:randomerlen
        deg,=poissrnd_degeneracy(lambda);
        bases=shuffle(nucleotides[deg])
        randomer[i]=bases[1];
    end
    push!(randpool,randomer)
    while length(randpool)<randpoolsize
        randomer=dna"-"^randomerlen
        for i = 1:randomerlen
            deg, =poissrnd_degeneracy(lambda);
            bases=shuffle(nucleotides[deg])
            randomer[i]=bases[1];
        end
        check=true
        for j = 1:length(randpool)
            if occursin(randomer,randpool[j])
                check=false
                break
            end
        end
        if check
            push!(randpool,randomer)
        end
        if length(randpool)==length(nucleotides[1])^randomerlen
            break
        end
    end

    return LongDNASeq.(randpool)
end

function poissrnd_degeneracy(lambda)
    pdf=zeros(4);
    cdf=zeros(4);
    for i=1:4
        pdf[i]=poisson_pmf(i,lambda)
    end
    pdf=pdf./sum(pdf)
    for i=1:4
        cdf[i]=sum(pdf[1:i])
    end
    r1=rand()
    scale=(cdf.-r1)
    for i=1:4
        if scale[i]>0
            return i , cdf
            continue
        end
    end
end
function poisson_pmf(x,lambda)
    Px=(lambda^x*exp(-lambda))/factorial(x)
end
function norm_degeneracy(lambda)
    pdfun=zeros(4);
    cdfun=zeros(4);
    for i=1:4
        values=cdf(Normal(lambda, .8),[i,i+1])
        pdfun[i]=values[2]-values[1]
    end
    pdfun=pdfun./sum(pdfun)
    for i=1:4
        cdfun[i]=sum(pdfun[1:i])
    end
    r1=rand()
    scale=(cdfun.-r1)
    for i=1:4
        if scale[i]>0
            return i ,cdfun
            continue
        end
    end
end

## Running function ##
"""
    oligo_pool_compressor(uncompressed_pool::Array{LongSequence{DNAAlphabet{4}},1},nucleotides::LongSequence{DNAAlphabet{4}}=dna"AGCTMRWSYKVHDBN";kwargs...)

Compress an uncompresssed_pool of non-degenerate oligos into a smaller compresssed_pool of degenerate oligos.

# Arguments
- `uncompressed_pool`: A set of non-degenerate oligos to be compressed.
- `nucleotides=dna"AGCTMRWSYKVHDBN"`: Base codes available to be used.
# Common Keyword Arguments
- `nsims`: number of rollout simulations per action
"""
function oligo_pool_compressor(uncompressed_pool::Array{LongSequence{DNAAlphabet{4}},1},nucleotides::LongSequence{DNAAlphabet{4}}=dna"AGCTMRWSYKVHDBN";kwargs...)
    sites=unique(uncompressed_pool);
    original_poolsize=length(sites)
    compressed_pool=[]
    randomerlen=length(sites[1]);
    ns= [nucleotides for _ in 1:randomerlen];
    while length(sites)>0

        randomer=oligo_rollout(ns,sites;kwargs...)
        sites=remove_hit_oligos(randomer,sites);
        push!(compressed_pool,randomer);
        #checkpoint=length(compressed_pool);
        #println("Compressed Pool Size: $checkpoint")
    end
    n=length(compressed_pool)
    failed_sites=[]
    for i = 1:length(compressed_pool)
        if degeneracy(compressed_pool[i];uselog2=false)==0
            push!(failed_sites,i)
        end
    end
    deleteat!(compressed_pool,failed_sites);
    new_poolsize=length(compressed_pool)
    println("Compressed Pool contains $new_poolsize oligos. Original Pool contained $original_poolsize oligos.")
    return compressed_pool
end

function benchmark_oligo_compressor(;all_bases=dna"AGCTMRWSYKVHDBN",kwargs...)
    nreps=10;
    testsizes=[100];
    poolsize=repeat(testsizes,nreps)
    nruns=length(poolsize)
    run=1:nruns
    time=zeros(nruns)
    compressedsize=zeros(nruns)
    pools=[]
    for i = 1:nruns
        randpool=generate_random_pool(;randomerlen=6, randpoolsize=poolsize[i])
        strpool=join(convert.(String,randpool),",");
        push!(pools,strpool)
        stats= @timed oligo_pool_compressor(randpool,all_bases;kwargs...)
        compressedsize[i]=length(stats.value)
        time[i]=stats.time
        filename="oligo_compressor_4_15_21.csv"
        data=DataFrame(poolsize=poolsize,compressedsize=compressedsize,time=time)
        #CSV.write(filename,data)
    end
    data=DataFrame(pool=pools,poolsize=poolsize,compressedsize=compressedsize,time=time)
    return data
end
function run_recompression_experiment(;nruns=50,all_bases = dna"AGCTMRWSYKVHDBN",kwargs...)
    sizes=repeat([1,5,10,15,20],10)
    runs=1:nruns
    pools=[]
    poolsize=zeros(nruns)
    decompressedsize=zeros(nruns)
    recompressedsize=zeros(nruns)
    for i = 1:nruns
        randpool=generate_poissrnd_pool(;randpoolsize=10, desiredsize=sizes[i])
        poolsize[i]=length(randpool)
        strpool=join(convert.(String,randpool),",");
        push!(pools,strpool)
        uncomrandpool=decompress_pool(randpool);
        decompressedsize[i]=length(uncomrandpool);
        recomprandpool=oligo_pool_compressor(uncomrandpool,all_bases;kwargs...)
        recompressedsize[i]=length(recomprandpool);
    end
    data=DataFrame(A=runs,B=pools,C=poolsize,D=decompressedsize,E=recompressedsize)
    rename!(data,[:runs,:pools,:poolsize,:decompressedsize,:recompressedsize])
    CSV.write("recompression_experiment_3-19-21.csv",data)
end
#benchmark_oligo_compressor()
function rerun_recompression_experiment(filename;nruns=50,all_bases = dna"AGCTMRWSYKVHDBN",kwargs...)
    sizes=repeat([1,5,10,15,20],10)
    runs=1:nruns

    data=CSV.read("recompression_experiment_3-19-21.csv",DataFrame)
    data[:recompressed100]=0;
    for i = 1:nruns
        randpool=LongDNASeq.(split(data[i,2],","))

        uncomrandpool=decompress_pool(randpool);
        recomprandpool=oligo_pool_compressor(uncomrandpool,all_bases;kwargs...)
        data[i,:recompressed100]=length(recomprandpool);
    end
    CSV.write(filename,data)
end

#rerun_recompression_experiment("recompression_3_19_21_add_nsims100.csv";nsims=100)
#data=benchmark_oligo_compressor(;nsims=100)

function run_compression_experiment()
    orig_pool=LongDNASeq.(CSV.read("./Oligo_Compressor/Experiments/Nature_2018_NSR_Primers.csv",DataFrame)[:Sequence])
    all_bases = dna"AGCTMRWSYKVHDBN"
    new_pool=oligo_pool_compressor(orig_pool,all_bases;nsims=100)
    data=DataFrame(randomers=new_pool)
    orig_len=length(orig_pool)
    new_len=length(new_pool)
    println("Original Length: $orig_len, New Length: $new_len")
    filename="./Oligo_Compressor/Experiments/Nature_2018_NSR_Primers_Compressed.csv"
    CSV.write(filename,data)
end


text:
# Contiguous on-device arrays

import Base: length, size, ndims, eltype, copy!

export
    DevicePtr, CuArray, free, to_host, ndims

type CuArray{T,N}
    ptr::DevicePtr{Void}
    shape::NTuple{N,Int}
    len::Int
end

function CuArray(T::Type, len::Integer)
    n = Int(len)
    p = cualloc(T, n)
    CuArray{T,1}(p, (n,), n)
end

function CuArray(T::Type, len::Integer, value::Cuint)
    n = Int(len)
    p = cualloc(T, n)
    cumemset(p, value, n)
    CuArray{T,1}(p, (n,), n)
end

function CuArray{N}(T::Type, shape::NTuple{N,Int})
    n = prod(shape)
    p = cualloc(T, n)
    CuArray{T,N}(p, shape, n)
end

function CuArray{N}(T::Type, shape::NTuple{N,Int}, value::Cuint)
    n = prod(shape)
    p = cualloc(T, n)
    cumemset(p, value, n)
    CuArray{T,N}(p, shape, n)
end

length(g::CuArray) = g.len
size(g::CuArray) = g.shape
ndims{T,N}(::CuArray{T,N}) = N
eltype{T,N}(::CuArray{T,N}) = T
eltype{T,N}(::Type{CuArray{T,N}}) = T

function size{T,N}(g::CuArray{T,N}, d::Integer)
    d >= 1 ? (d <= N ? g.shape[d] : 1) : error("Invalid index of dimension.")
end

function free(g::CuArray)
    if !isnull(g.ptr)
        free(g.ptr)
        g.ptr = DevicePtr(C_NULL)
    end
end

function copy!{T}(dst::Array{T}, src::CuArray{T})
    if length(dst) != length(src)
        throw(ArgumentError("Inconsistent array length."))
    end
    nbytes = length(src) * sizeof(T)
    @cucall(:cuMemcpyDtoH, (Ptr{Void}, DevicePtr{Void}, Csize_t),
                           pointer(dst), src.ptr, nbytes)
    return dst
end

function copy!{T}(dst::CuArray{T}, src::Array{T})
    if length(dst) != length(src)
        throw(ArgumentError("Inconsistent array length."))
    end
    nbytes = length(src) * sizeof(T)
    @cucall(:cuMemcpyHtoD, (DevicePtr{Void}, Ptr{Void}, Csize_t),
                           dst.ptr, pointer(src), nbytes)
    return dst
end

CuArray{T,N}(a::Array{T,N}) = copy!(CuArray(T, size(a)), a)
to_host{T}(g::CuArray{T}) = copy!(Array(T, size(g)), g)


text:
#=
Compare call option prices to consol price.

=#

using Plots
include("asset_pricing.jl")
plotly()

ap = AssetPriceModel(beta=0.9)
zeta = 1.0
strike_price = 40.0

x = ap.mc.state_values
p = consol_price(ap, zeta)
w = call_option(ap, zeta, strike_price)

plot(x, p, color="blue", lw=2, xlabel="state", label="consol price")
plot!(x, w, color="green", lw=2, label="value of call option")

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
using RunTests
using Base.Test

@testmodule TestModuleWithOneXFailingTest begin

@xfail function test_xfailing()
  @test false
end

end


text:
using Jecco.AdS5_3_1

grid = SpecCartGrid3D(
    x_min            = -5.0,
    x_max            =  5.0,
    x_nodes          =  12,
    y_min            = -5.0,
    y_max            =  5.0,
    y_nodes          =  12,
    u_outer_min      =  0.1,
    u_outer_max      =  1.003,
    u_outer_domains  =  3,
    u_outer_nodes    =  24,
    u_inner_nodes    =  12,
    fd_order         =  4,
    sigma_diss       =  0.2,
)

id   = BlackBrane(
    AH_pos = 1.001,
)

evoleq = AffineNull(
    phi0       = 0.0,
    potential  = ZeroPotential(),
    gaugecondition = ConstantAH(u_AH = 1.0),
)

io = InOut(
    out_boundary_every  = 10,
    out_bulk_every      = 100,
    out_gauge_every     = 10,
    # remove_existing     = true,
)

integration = Integration(
    dt              = 0.001,
    tmax            = 5.0,
    ODE_method      = AdS5_3_1.AB3(),
    filter_poststep = true,
)

run_model(grid, id, evoleq, integration, io)


text:
# copyright Lénaïc Chizat
# the code shows the training dynamics in parameter and predictor space for a 2-layer relu NN on a two class classification task
# we optimize the exponential loss (the behavior is sensibly the same with the logistic loss)
# we use a specific step-size schedule motivated by the theory in this ref https://arxiv.org/abs/2002.04486 (sections 4 & 5)

# HOW TO RUN THE CODE
# in a prompt first run:
# include("exp_output.jl")
# Then for fast prototyping, run for instance the following code (takes 1 min)
# illustration(4, 60, 200, 4, 100000, 10, 0.05);
# These parameters are explained in the header of the function "illustration"
# A good illustration is obtained as follows (takes 20 min):
# Random.seed!(9); # 9 generates a nice random problem
# illustration(4, 60, 400, 1, 3000000, 400, 0.005);

using LinearAlgebra, Random
using PyPlot, ProgressMeter # these packages need to be installed (via "] add NamePackage" )
using3D()

"""
Gradient descent to train a 2-layers ReLU neural network for the exponential loss
We use the step-size schedule from in https://arxiv.org/abs/2002.04486 (sections 4 & 5)
INPUT: X (training input), Y (training output), m (nb neurons), both: training
        both layers or just the output layer
OUTPUT: Ws (training trajectory)
"""
function twonet(X, Y, m, stepsize, niter)

    (n,d) = size(X) # n samples in R^d
    # initialize
    W = randn(m, d+1)
    # input weights initialized uniformly on the unit sphere
    W[:,1:d] = W[:,1:d] ./ sqrt.(sum(W[:,1:d].^2, dims=2))
    # output weights initialized by -1 or 1 with equal proba
    W[1:div(m,2),end] .= 0.0
    W[div(m,2)+1:end,end] .= 0.0

    Ws    = zeros(m, d+1, niter) # store optimization path
    loss  = zeros(niter) # loss is -log of the empirical risk
    margins = zeros(niter)
    betas = zeros(niter)

    @showprogress 1 "Training neural network..." for iter = 1:niter
        Ws[:,:,iter] = W
        act  =  max.( W[:,1:end-1] * X', 0.0) # activations
        out  =  (1/m) * sum( W[:,end] .* act , dims=1) # predictions of the network
        perf = Y .* out[:]
        margin = minimum(perf)
        temp = exp.(margin .- perf) # numerical stabilization of exp
        gradR = temp .* Y ./ sum(temp)' # gradient of the loss
        grad_w1 = zeros(m,d) #(W[:,end] .* float.(act .> 0) * ( X .* gradR  ))  # gradient for input weights
        grad_w2 = act * gradR  # gradient for output weights

        grad = cat(grad_w1, grad_w2, dims=2) # size (m × d+1)
        betas[iter] = max(1, sqrt(sum(W[:,end].^2)/m))
        loss[iter] = margin - log(sum(exp.(margin .- perf))/n)
        margins[iter] = margin/betas[iter] # margin in F_1 norm
        W = W + stepsize * grad * betas[iter]/(sqrt(iter+1))
    end

    return Ws, loss, margins, betas
end

##
"Coordinates of the 2d cluster centers, p is k^2 the number of clusters"
function cluster_center(p,k)
    p1 = mod.(p .- 1,k) .+ 1
    p2 = div.(p .- 1,k) .+ 1
    Delta = 1/(3k-1)
    x1 =  Delta*(1 .+ 3*(p1 .- 1)) .- 1/2
    x2 =  Delta*(1 .+ 3*(p2 .- 1)) .- 1/2
    return x1, x2
end

"""
Plot the classifier for a test case, comparing training both  layer
k: number of clusters per dimension in the data set (choose 3, 4 or 5)
n: number of training samples
m: number of neurons
stepsize, niter: parameters for training
nframes, resolution: parameters for the plots
"""
function illustration(k, n, m, stepsize, niter, nframes, resolution)
    # data distribution (itself random)
    sd = 0 # number of spurious dimensions with random noise
    Delta = 1/(3k-1) # interclass distance
    A = ones(k^2) # cluster affectation
    A[randperm(k^2)[1:div(k^2,2)]] .= -1

    # sample from this data distribution
    P = rand(1:k^2,n) # cluster label
    T = 2π*rand(n)  # shift angle
    R = Delta*rand(n) # shift magnitude
    X = cat(ones(n), cluster_center(P,k)[1] .+ R .* cos.(T),cluster_center(P,k)[2] + R .* sin.(T), (rand(n,sd) .- 1/2), dims=2)
    Y = A[P]

    # plot training set
    X1 = X[(Y .== 1),:]
    X2 = X[(Y .== -1),:]
    fig = figure(figsize=[3,3])
    plot(X1[:,2],X1[:,3],"+k")
    plot(X2[:,2],X2[:,3],"_k")
    axis("equal");axis("off");
    display(fig)

    # train the neural network
    Ws, loss, margins, betas = twonet(X, Y, m, stepsize, niter)

    # define the sequence of time steps ts to be plotted
    a = (niter-1)/(nframes-1)^4
    ts = setdiff(Int.(floor.(a*(0:nframes-1).^4)) .+ 1)
    Ws = Ws[:,:,ts]

    # these are the positions we plot for the particles
    Wproj = Ws[:,1:end-1,:] .* abs.(Ws[:,end:end,:])
    WN   = sqrt.(sum(Wproj.^2, dims = 2))
    Wdir = Wproj ./ WN
    Wlog = tanh.(0.5*WN) .* Wdir

    @showprogress 1 "Plotting images..." for k = 1:length(ts)
        ioff() # turns off interactive plotting
        fig = figure(figsize=[7,4])
        ax1 = subplot(121, projection="3d")
        ax1.set_position([0,0.1,0.5,0.8])

        if k<11
            indt = 1:k
        else
            indt = (k-10):k
        end

        for i = 1:size(Wlog,1)
            plot3D(Wlog[i,2,indt],Wlog[i,3,indt],Wlog[i,1,indt], color="k", linewidth=0.2) # tail
        end
        ind_p = (Ws[:,end,k] .>= 0.0)[:]
        ind_m = (Ws[:,end,k] .< 0.0)[:]
        plot3D(Wlog[ind_p,2,k],Wlog[ind_p,3,k],Wlog[ind_p,1,k],"o",color="C3", markersize=1)
        plot3D(Wlog[ind_m,2,k],Wlog[ind_m,3,k],Wlog[ind_m,1,k],"o",color="C0", markersize=1)

        ax1.set_xlim3d(-1, 1)
        ax1.set_ylim3d(-1, 1)
        ax1.set_zlim3d(-1, 1)
        ax1.set_xticks([-1/2, 0, 1/2])
        ax1.set_yticks([-1/2, 0, 1/2])
        ax1.set_zticks([-1/2, 0, 1/2])
        ax1.view_init(25-20*sinpi(k/length(ts)),45-45*cospi(k/length(ts)))

        ax2 = subplot(122)
        ax2.set_position([0.45,0.25,0.5,0.5])

        f(x1,x2,k) = (1/m) * sum( Ws[:,end,k] .* max.( Ws[:,1:3,k] * [1;x1;x2], 0.0)) # prediction function

        xs = -0.8:resolution:0.8
        tab = [f(xs[i],xs[j],k) for i=1:length(xs), j=1:length(xs)]
        pcolormesh(xs', xs, tanh.(tab'), cmap="coolwarm", shading="gouraud", vmin=-1.0, vmax=1.0, edgecolor="face")

        xs = -0.8:resolution:0.8
        tab = [f(xs[i],xs[j],k) for i=1:length(xs), j=1:length(xs)]
        contour(xs', xs, tanh.(tab'), levels =0, colors="k", antialiased = true, linewidths=2)

        # plot training set
        X1 = X[(Y .== 1),:]
        X2 = X[(Y .== -1),:]
        plot(X1[:,2],X1[:,3],"+k")
        plot(X2[:,2],X2[:,3],"_k")
        axis("equal");axis("off");
        ax2.set_xticks([-1/2, 0, 1/2])
        ax2.set_yticks([-1/2, 0, 1/2])

        savefig("dynamics_output_$(k).png",bbox_inches="tight", dpi=300)
        close(fig)
    end

end

##

max_stars_repo_name: brechetp/ImplicitBias

max_issues_repo_name: brechetp/ImplicitBias

max_forks_repo_name: brechetp/ImplicitBias


text:
"""
Model for ZIP Load model given by:

P_zip = P_power + P_current * (V / V0) + P_impedance * (V / V0)^2
Q_zip = Q_power + Q_current * (V / V0) + Q_impedance * (V / V0)^2

with V = sqrt(V_r^2 + V_i^2) and V0 the voltage magnitude from the power flow solution

The current taken for the load is computed as:
I_zip = (P_zip + j Q_zip)^* / (V_r + j V_i)^*
I_zip = (P_zip - j Q_zip) / (V_r - j V_i)

For constant impedance it is obtained:
Iz_re = (1 / V0)^2 * (V_r * P_impedance + V_i * Q_impedance)
Iz_im = (1 / V0)^2 * (V_i * P_impedance - V_r * Q_impedance)

For constant current it is obtained:
Ii_re = (1 / V0) * ( (V_r * P_current + V_i * Q_current) / V )
Ii_im = (1 / V0) * ( (V_i * P_current - V_r * Q_current) / V )

For constant power it is obtained:
Ip_re =  (V_r * P_power + V_i * Q_power) / V^2
Ip_im =  (V_i * P_power - V_r * Q_power) / V^2

Model for Exponential Load model given by:

P_exp = P0 * (V / V0)^α
Q_exp = Q0 * (V / V0)^β

The current taken for the load is computed as:
I_exp = (P_exp + j Q_exp)^* / (V_r + j V_i)^*
I_exp = (P_exp - j Q_exp) / (V_r - j V_i)

It results:
Ir_exp = V_r * P0 * (V^(α - 2) / V0^α) + V_i * Q0 * (V^(β - 2)/ V0^β)
Ii_im  = V_i * P0 * (V^(α - 2) / V0^α) - V_r * Q0 * (V^(β - 2)/ V0^β)

"""
function mdl_zip_load!(
    voltage_r::T,
    voltage_i::T,
    current_r::AbstractArray{T},
    current_i::AbstractArray{T},
    wrapper::StaticLoadWrapper,
) where {T <: ACCEPTED_REAL_TYPES}
    # Read power flow voltages
    #V0_mag_inv = 1.0 / get_V_ref(wrapper)
    V0_mag_inv = 1.0 / PSY.get_magnitude(PSY.get_bus(wrapper))
    V0_mag_sq_inv = V0_mag_inv^2

    V_mag = sqrt(voltage_r^2 + voltage_i^2)
    V_mag_inv = 1.0 / V_mag
    V_mag_sq_inv = V_mag_inv^2

    # Load device parameters
    P_power = get_P_power(wrapper)
    P_current = get_P_current(wrapper)
    P_impedance = get_P_impedance(wrapper)
    Q_power = get_Q_power(wrapper)
    Q_current = get_Q_current(wrapper)
    Q_impedance = get_Q_impedance(wrapper)
    exp_params = get_exp_params(wrapper)
    Ir_exp = zero(T)
    Ii_exp = zero(T)

    # Compute ZIP currents
    Iz_re = V0_mag_sq_inv * (voltage_r * P_impedance + voltage_i * Q_impedance)
    Iz_im = V0_mag_sq_inv * (voltage_i * P_impedance - voltage_r * Q_impedance)

    Ii_re = V0_mag_inv * V_mag_inv * (voltage_r * P_current + voltage_i * Q_current)
    Ii_im = V0_mag_inv * V_mag_inv * (voltage_i * P_current - voltage_r * Q_current)

    Ip_re = V_mag_sq_inv * (voltage_r * P_power + voltage_i * Q_power)
    Ip_im = V_mag_sq_inv * (voltage_i * P_power - voltage_r * Q_power)

    # Compute Exponential Currents
    if !isempty(exp_params)
        for tuple in exp_params
            P0 = tuple.P_exp
            α = tuple.P_coeff
            Q0 = tuple.Q_exp
            β = tuple.Q_coeff

            # Compute Auxiliary terms
            V_coeff_active = V_mag^(α - 2.0) * V0_mag_inv^α
            V_coeff_reactive = V_mag^(β - 2.0) * V0_mag_inv^β

            # Compute currents
            Ir_exp += voltage_r * P0 * V_coeff_active + voltage_i * Q0 * V_coeff_reactive
            Ii_exp += voltage_i * P0 * V_coeff_active - voltage_r * Q0 * V_coeff_reactive
        end
    end

    # Update current
    current_r[1] += -(Iz_re + Ii_re + Ip_re + Ir_exp) #in system pu flowing out
    current_i[1] += -(Iz_im + Ii_im + Ip_im + Ii_exp) #in system pu flowing out

    return
end


text:
#=============================#
#== put needed modules here ==#
#=============================#
import JuliaChem

#================================#
#== JuliaChem execution script ==#
#================================#
function travis_rhf(input_file)
  try
    #== read in input file ==#
    molecule, driver, model, keywords = JuliaChem.JCInput.run(input_file;
      output="verbose")

    #== generate basis set ==#
    mol, basis = JuliaChem.JCBasis.run(molecule, model; output="verbose")

    #== compute molecular inform ation ==#
    JuliaChem.JCMolecule.run(mol)

    #== perform scf calculation ==#
    rhf_energy = JuliaChem.JCRHF.Energy.run(mol, basis, keywords["scf"];
      output="verbose")

    #== compute molecular properties ==#
    rhf_properties = JuliaChem.JCRHF.Properties.run(mol, basis, rhf_energy, keywords["prop"],
      output="verbose")

    return (Energy = rhf_energy, Properties = rhf_properties)
  catch e
    bt = catch_backtrace()
    msg = sprint(showerror, e, bt)
    println(msg)

    JuliaChem.finalize()
    exit()
  end
end


text:
using HwlocViz
using Test

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

max_stars_repo_name: crstnbr/HwlocViz

max_issues_repo_name: crstnbr/HwlocViz

max_forks_repo_name: crstnbr/HwlocViz


text:
__precompile__()

module Hymod

using DataFrames
using Dates
using Distributions
using Statistics
using Logging
using JSON

export randomparams, hargreaves, simulate, calibrate

"""
    randomparams()

Function to create a random set of parameters as kwargs for Hymod.simulate()
"""
function randomparams()

    param_bounds = Dict{Symbol,Dict}(
        :cmax => Dict{Symbol,Float64}(:lower => 1.0, :upper => 100),
        :bexp => Dict{Symbol,Float64}(:lower => 0.0, :upper => 2.0),
        :alpha => Dict{Symbol,Float64}(:lower => 0.2, :upper => 0.99),
        :ks => Dict{Symbol,Float64}(:lower => 0.01, :upper => 0.5),
        :kq => Dict{Symbol,Float64}(:lower => 0.5, :upper => 1.2),
    )

    out = Dict{Symbol,Float64}()

    for k in keys(param_bounds)
        minV = param_bounds[k][:lower]
        maxV = param_bounds[k][:upper]
        p = rand(Uniform(minV, maxV), 1)
        out[k] = p[1]
    end

    return out
end

function abspower(x::Float64, y::Float64)
    # Needed to capture invalid overflow with netgative values
    x = abs(x)
    return x^y
end

function precipexcess(
    x_loss::Float64,
    cmax::Float64,
    bexp::Float64,
    prc::Float64,
    pet::Float64,
)
    xn_prev = x_loss
    ct_prev =
        cmax *
        (1 - abspower((1 - ((bexp + 1) * (xn_prev) / cmax)), (1 / (bexp + 1))))
    # Calculate Effective rainfall 1
    ER1 = max((prc - cmax + ct_prev), 0.0)
    prc = prc - ER1
    dummy = min(((ct_prev + prc) / cmax), 1)
    xn = (cmax / (bexp + 1)) * (1 - abspower((1 - dummy), (bexp + 1)))

    # Calculate Effective rainfall 2
    ER2 = max(prc - (xn - xn_prev), 0)

    # Alternative approach
    # actual ET is linearly related to the soil moisture state
    evap = (1 - (((cmax / (bexp + 1)) - xn) / (cmax / (bexp + 1)))) * pet
    xn = max(xn - evap, 0)  # update state

    return ER1, ER2, xn
end

"""
    randomparams(xslow::Float64, inflow::Float64, Rs::Float64)

Linear reservoir routing function
"""
function linearreservoir(xslow::Float64, inflow::Float64, Rs::Float64)
    # Linear reservoir routing
    xslow = (1 - Rs) * xslow + (1 - Rs) * inflow
    outflow = (Rs / (1 - Rs)) * xslow
    return xslow, outflow
end

"""
    hargreaves(forcings::DataFrame; tmincol::Symbol = :tmin, tmaxcol::Symbol = :tmax, dtcol::Symbol = :datetime)

Function to calculate PET using Hargreves equation.
"""
function hargreaves(
    forcings::DataFrame;
    tmincol::Symbol = :tmin,
    tmaxcol::Symbol = :tmax,
    dtcol::Symbol = :datetime,
)

    dts = forcings[!, dtcol]
    tmin = forcings[!, tmincol]
    tmax = forcings[!, tmaxcol]

    hargreaves(tmin, tmax, dts)
end

"""
    hargreaves(tmin::AbstractArray, tmax::AbstractArray, dts::AbstractArray)

Function to calculate PET using Hargreves equation.
"""
function hargreaves(
    tmin::AbstractArray,
    tmax::AbstractArray,
    dts::AbstractArray,
)
    hargreaves.(tmin, tmax, dts)
end

"""
    hargreaves(tmin::Real, tmax::Real, dts::Date)

Function to calculate PET using Hargreves equation.
"""
function hargreaves(tmin::Real, tmax::Real, dts::Date)
    Gsc = 367
    lhov = 2.257

    doy = dayofyear(dts)

    tavg = mean((tmin, tmax))

    b = 2 * pi * (Int16(doy) / 365)
    Rav =
        1.00011 +
        0.034221 * cos(b) +
        0.00128 * sin(b) +
        0.000719 * cos(2 * b) +
        0.000077 * sin(2 * b)
    Ho = ((Gsc * Rav) * 86400) / 1e6

    eto = (0.0023 * Ho * (tmax - tmin)^0.5 * (tavg + 17.8))

    return eto
end

"""
    simulate(forcings::DataFrame; precipcol::Symbol = :precip, petcol::Symbol = :pet, kwargs...)

Function to simulate discharge using Hymod rainfall-runoff model.
"""
function simulate(
    forcings::DataFrame;
    precipcol::Symbol = :precip,
    petcol::Symbol = :pet,
    kwargs...,
)
    p = forcings[:, precipcol]
    e = forcings[:, petcol]

    simulate(p, e; kwargs...)

end

"""
    simulate(precip::AbstractArray, pet::AbstractArray; initflow::Bool = true, cmax::Float64 = 1.0, bexp::Float64 = 0.0, alpha::Float64 = 0.2, ks::Float64 = 0.01, kq::Float64 = 0.5)

Function to simulate discharge using Hymod rainfall-runoff model
"""
function simulate(
    precip::AbstractArray,
    pet::AbstractArray;
    initflow::Bool = true,
    cmax::Float64 = 1.0,
    bexp::Float64 = 0.0,
    alpha::Float64 = 0.2,
    ks::Float64 = 0.01,
    kq::Float64 = 0.5,
)

    lt_to_m::Float64 = 0.001

    n_iters = length(precip)

    # HYMOD PROGRAM IS SIMPLE RAINFALL RUNOFF MODEL
    x_loss::Float64 = 0.0
    # Initialize slow tank state
    # value of 0 init flow works ok if calibration data starts with low discharge
    x_slow::Float64 = initflow ? (2.3503 / (ks * 22.5)) : 0.0
    # Initialize state(s) of quick tank(s)
    x_quick = zeros(Float64, 3)
    outflow = zeros(Float64, n_iters)
    output = zeros(Float64, n_iters)
    # START PROGRAMMING LOOP WITH DETERMINING RAINFALL - RUNOFF AMOUNTS

    # prealloc values within loops
    Pval::Float64 = 0.0
    PETval::Float64 = 0.0
    ER1::Float64 = 0.0
    ER2::Float64 = 0.0
    UQ::Float64 = 0.0
    US::Float64 = 0.0
    QS::Float64 = 0.0
    inflow::Float64 = 0.0

    @inbounds for i = 1:n_iters
        Pval = precip[i]
        PETval = pet[i]
        # Compute excess precipitation and evaporation
        ER1, ER2, x_loss = precipexcess(x_loss, cmax, bexp, Pval, PETval)
        # Calculate total effective rainfall
        ET = ER1 + ER2
        #  Now partition ER between quick and slow flow reservoirs
        UQ = alpha * ET
        US = (1 - alpha) * ET
        # Route slow flow component with single linear reservoir
        x_slow, QS = linearreservoir(x_slow, US, ks)
        # Route quick flow component with linear reservoirs
        inflow = UQ

        for j in eachindex(x_quick)
            # Linear reservoir
            x_quick[j], outflow = linearreservoir(x_quick[j], inflow, kq)
            inflow = outflow
        end

        # Compute total flow for timestep
        output[i] = ((QS + outflow) / lt_to_m)
    end

    output[output.<0] .= 0

    return output
end

"""
    calibrate(forcing::DataFrame, paramspace::Dict, nsamples::Int64; precipcol::Symbol = :precip, petcol::Symbol = :pet, obscol::Symbol = :obs, savefinal::Bool = false)

Function to calibrate Hymod rainfall-runoff model
"""
function calibrate(
    forcing::DataFrame,
    paramspace::Dict,
    nsamples::Int64;
    precipcol::Symbol = :precip,
    petcol::Symbol = :pet,
    obscol::Symbol = :obs,
    savefinal::Bool = false,
)
    # implementation of a monte carlo sampling

    obs = forcing[!, obscol]

    params = Dict{Symbol,Array}()

    keyList = collect(keys(paramspace))

    for k in keyList
        p = paramspace[k]
        minV = p[:lower]
        maxV = p[:upper]
        params[k] = collect(rand(Uniform(minV, maxV), nsamples))
    end

    nIter = nsamples
    losses = zeros(nIter)

    @info "Calibrating model with $nIter iterations..."
    for i = 1:nIter
        vals = [params[k][i] for k in keyList]
        pars = Dict(keyList .=> vals)
        q = simulate(forcing, precipcol = :precip, petcol = :pet; pars...)
        # only test after one year of spinup
        losses[i] = nse(q[365:end], obs[365:end])
    end

    # check to make sure that the loss metric is a valid number
    losses = convert(Vector{Union{Missing,Float64}}, losses)
    losses[map(isnan, losses)] .= missing

    finalLoss, idx = findmax(skipmissing(losses))
    vals = [params[k][idx] for k in keyList]
    finalPars = Dict(keyList .=> vals)
    finalQ = simulate(forcing, precipcol = :precip, petcol = :pet; finalPars...)

    if savefinal
        t = now()
        finalPars[:loss] = finalLoss
        jsondata = Dict(String(k) => v for (k, v) in finalPars)
        @debug (jsondata)
        open("hymod_calibration_results_$t.json", "w") do io
            JSON.print(io, jsondata, 4)
        end
    end

    return finalQ, finalPars, finalLoss

end

"""
    nse(y_true::AbstractArray, y_pred::AbstractArray)

Function to calculate Nash-Sutcliffe model efficiency coefficient
"""
function nse(y_true::AbstractArray, y_pred::AbstractArray)
    # sum square difference between observed and simulated
    numerator = sum((y_true - y_pred) .^ 2)
    # sum square difference between observed and observed mean
    denominator = sum((y_true .- mean(y_true)) .^ 2)

    return 1 - (numerator / denominator)
end

end # module