text:
function difference(x::VersionNumber, y::VersionNumber)
    if y.major > x.major
        return VersionNumber(y.major - x.major, y.minor, y.patch)
    elseif y.minor > x.minor
        return VersionNumber(y.major - x.major, y.minor - x.minor, y.patch)
    elseif y.patch > x.patch
        return VersionNumber(y.major - x.major, y.minor - x.minor, y.patch - x.patch)
    else
        throw(ArgumentError("first argument must be strictly less than the second argument"))
    end
end

function latest_version(pkg::String, registry_path::String)
    all_versions = VersionNumber.(keys(Pkg.TOML.parsefile(joinpath(registry_path, uppercase(pkg[1:1]), pkg, "Versions.toml"))))
    return maximum(all_versions)
end

function julia_compat(pkg::String, version::VersionNumber, registry_path::String)
    all_compat_entries_for_julia = Vector{Pkg.Types.VersionRange}(undef, 0)
    compat = Pkg.TOML.parsefile(joinpath(registry_path, uppercase(pkg[1:1]), pkg, "Compat.toml"))
    for version_range in keys(compat)
        if version in Pkg.Types.VersionRange(version_range)
            for compat_entry in compat[version_range]
                name = compat_entry[1]
                if strip(lowercase(name)) == strip(lowercase("julia"))
                    value = compat_entry[2]
                    if value isa Vector
                        for x in value
                            x_range = Pkg.Types.VersionRange(x)
                            push!(all_compat_entries_for_julia, x_range)
                        end
                    else
                        value_range = Pkg.Types.VersionRange(value)
                        push!(all_compat_entries_for_julia, value_range)
                    end
                end
            end
        end
    end
    if length(all_compat_entries_for_julia) < 1
        return Pkg.Types.VersionRange[Pkg.Types.VersionRange("* - *")]
    else
        return all_compat_entries_for_julia
    end
end

function _has_upper_bound(r::Pkg.Types.VersionRange)
    return r.upper != Pkg.Types.VersionBound("*")
end

function range_did_not_narrow(r1::Pkg.Types.VersionRange, r2::Pkg.Types.VersionRange)
    result = !Pkg.Types.isless_ll(r1.lower, r2.lower) && !Pkg.Types.isless_uu(r2.upper, r1.upper)
    return result
end

function range_did_not_narrow(v1::Vector{Pkg.Types.VersionRange}, v2::Vector{Pkg.Types.VersionRange})
    @debug("", v1, v2, repr(v1), repr(v2))
    if isempty(v1) || isempty(v2)
        return false
    else
        n_1 = length(v1)
        n_2 = length(v2)
        results = falses(n_1, n_2) # v1 along the rows, v2 along the columns
        for i = 1:n_1
            for j = 1:n_2
                results[i, j] = range_did_not_narrow(v1[i], v2[j])
            end
        end
        @debug("", results, repr(results))
        return all(results)
    end
end


text:
# a. basic function
function print_me(input="Nothing given :(")
    println("> ", input)
end
print_me("cool")
print_me()
#: > cool
#: > Noting given :(

ext: jl


text:
using TravelingSalesmanExact, GLPK, Test
set_default_optimizer!(GLPK.Optimizer)

function test_valid_tour(tour, L)
    @test length(tour) == L
    @test isperm(tour)
    @test tour isa Vector{Int}
end

function tour_cost(tour, cost::AbstractMatrix)
    n = length(tour)
    inc(a) = a == n ? one(a) : a + 1
    s = zero(eltype(cost))
    for (index, city) in enumerate(tour)
        nextcity = tour[inc(index)]
        s += cost[city, nextcity]
    end
    return s
end

function tour_cost(
    tour,
    cities::AbstractVector;
    distance = TravelingSalesmanExact.euclidean_distance,
)
    n = length(cities)
    cost = [distance(cities[i], cities[j]) for i = 1:n, j = 1:n]
    return tour_cost(tour, cost)
end

function test_tour(input, opt = TravelingSalesmanExact.get_default_optimizer(); kwargs...)
    t, c = get_optimal_tour(input, opt; kwargs...)
    test_valid_tour(t, size(input, 1))
    if :distance in keys(kwargs)
        @test c ≈ tour_cost(t, input; distance = kwargs[:distance])
    else
        @test c ≈ tour_cost(t, input)
    end
    return t, c
end

@testset "Small random asymmetric" begin
    # chosen via rand(5,5)
    cost = 10 *
           [
        0.474886 0.350983 0.262651 0.138455 0.904042;
        0.683586 0.922968 0.278874 0.408406 0.0224372;
        0.513651 0.778167 0.140392 0.981211 0.891122;
        0.20529 0.976361 0.784706 0.98504 0.385203;
        0.489131 0.783738 0.538762 0.998821 0.0324331
    ]
    t1, c1 = test_tour(cost; verbose = true)
    t2, c2 = test_tour(cost; verbose = false)
    t3, c3 = test_tour(cost; verbose = false, lazy_constraints = true)
    @test c1 ≈ c2
    @test c2 ≈ c3

    # incorrect `symmetric` should give the wrong answers
    t4, c4 = get_optimal_tour(cost; verbose = false, symmetric = true)
    test_valid_tour(t4, 5)
    @test !(c2 ≈ c4)
    @test !(c4 ≈ tour_cost(t4, cost))
end

@testset "Medium random asymmetric" begin
    cost = 10 * rand(15, 15)
    t1, c1 = test_tour(cost; verbose = false)
    t2, c2 = test_tour(cost; verbose = true, lazy_constraints = true)
    @test c1 ≈ c2
end

@testset "Small random symmetric" begin
    cost = 10 *
           [
        0.474886 0.350983 0.262651 0.138455 0.904042;
        0.683586 0.922968 0.278874 0.408406 0.0224372;
        0.513651 0.778167 0.140392 0.981211 0.891122;
        0.20529 0.976361 0.784706 0.98504 0.385203;
        0.489131 0.783738 0.538762 0.998821 0.0324331
    ]
    cost_sym = cost + transpose(cost)
    t4, c4 = test_tour(cost_sym; verbose = false)
    t5, c5 = test_tour(cost_sym; symmetric = true, verbose = false)
    t6, c6 = test_tour(cost_sym; symmetric = false, verbose = false)
    t7, c7 = test_tour(cost_sym; symmetric = false, verbose = true, lazy_constraints = true)
    @test c4 ≈ c5
    @test c5 ≈ c6
    @test c6 ≈ c7
end

@testset "Small random cities" begin
    # cities = [ 100*rand(2) for _ in 1:5]
    cities = Array{Float64,1}[
        [48.8885, 41.0517],
        [35.6635, 12.1844],
        [95.6122, 15.9847],
        [67.5772, 9.54407],
        [16.6325, 51.9001],
    ]
    t7, c7 = test_tour(cities; verbose = true)
    t8, c8 = @inferred test_tour(cities; verbose = false, symmetric = false)
    cost = [TravelingSalesmanExact.euclidean_distance(c1, c2) for c1 in cities, c2 in cities]
    t9, c9 = test_tour(cost; verbose = false, symmetric = false)
    t10, c10 = test_tour(cost; verbose = false, symmetric = true)
    @test c7 ≈ c8
    @test c8 ≈ c9
    @test c9 ≈ c10
    @test c8 isa Float64
end

@testset "Medium random cities" begin
    # cities = [ 100*rand(2) for _ in 1:15]
    cities = Array{Float64,1}[
        [53.80996646037832, 19.19022341362766],
        [66.31428678667224, 61.80858818555463],
        [36.353148448445616, 41.255152847278474],
        [90.27061139767463, 6.978822544781482],
        [53.764677734646504, 29.81193167526588],
        [99.44803233011115, 66.55866205554817],
        [38.928236331570986, 18.22266760966855],
        [21.194033394576728, 56.14389699364914],
        [6.430489908372161, 92.08608470530129],
        [67.01950679275805, 53.82196541323463],
        [57.04078593130859, 88.74875672836768],
        [40.91801899608445, 23.47673907020007],
        [8.889246634653713, 44.183494522157886],
        [82.70468776751652, 74.73476674725042],
        [8.637366704185245, 87.23819026270408]
    ]
    t7, c7 = test_tour(cities; verbose = true, lazy_constraints = true)
    t8, c8 = @inferred test_tour(cities; verbose = false, symmetric = false)
    cost = [TravelingSalesmanExact.euclidean_distance(c1, c2) for c1 in cities, c2 in cities]
    t9, c9 = test_tour(cost; verbose = false, symmetric = false)
    t10, c10 = test_tour(cost; verbose = false, symmetric = true)
    @test c7 ≈ c8
    @test c8 ≈ c9
    @test c9 ≈ c10
    @test c8 isa Float64
end

@testset "Exceptions" begin
    cost = rand(5, 4)
    @test_throws ArgumentError test_tour(cost)

    cost = rand(5, 5)
    TravelingSalesmanExact.reset_default_optimizer!()
    @test_throws ArgumentError test_tour(cost)
    set_default_optimizer!(GLPK.Optimizer)
end

@testset "att48.tsp" begin
    cities = simple_parse_tsp(joinpath(@__DIR__, "..", "data", "att48.tsp"))

    @test cities == TravelingSalesmanExact.get_ATT48_cities()

    sym_tour, sym_cost = test_tour(
        cities;
        distance = TravelingSalesmanExact.ATT,
        verbose = true,
    )
    @test sym_cost ≈ 10628

    asym_tour, asym_cost = test_tour(
        cities;
        distance = TravelingSalesmanExact.ATT,
        symmetric = false,
    )
    @test asym_cost ≈ 10628
end

@testset "Time Windows" begin
    cost = zeros(4, 4)
    cost[1,2] = 1
    cost[2,1] = 0
    cost[1,4] = -1

    times = ones(4, 4)
    tw = [
        0 4
        0 1
        0 3
        0 2
    ]
    route, objective = get_optimal_tour(cost, times, tw)
    @test objective ≈ 1
    start = findfirst(route .== 1)
    route = circshift(route, 1-start)
    @test route == [1,2,4,3]

    tw[2,:] .= [0,2]
    route, objective = get_optimal_tour(cost, times, tw)
    @test objective ≈ -1
    times = [
        0 1 0 0
        10 0 10 10
        10 10 0 10
        10 10 10 0
    ]
    route, objective = get_optimal_tour(times, tw)
    @test objective ≈ Inf
    tw = [
        0 100
        0 1
        0 100
        0 100
    ]
    route, objective = get_optimal_tour(times, tw)
    @test objective ≈ 31
end

@testset "Dropped" begin
    cost = ones(10, 10)
    penalties = Dict(
        1 => 3,
        3 => 0,
    )
    route, objective = get_optimal_tour(cost, penalties)
    @test length(route) == 9
    @test objective ≈ 9

    times = [
        0 1 0 0
        10 0 10 10
        10 10 0 10
        10 10 10 0
    ]
    tw = [
        0 100
        0 1
        0 100
        0 100
    ]
    penalties = Dict(2 => 5)
    route, objective = get_optimal_tour(times, tw, penalties)
    @test objective ≈ 25
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to wrong arguments.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Errors" begin
    @test_throws Exception pretty_table(data, vlines = :nothing)
    @test_throws Exception pretty_table(data, hlines = :nothing)
end


text:
## make options for the stationary multichannel lengthscale GP
nminlog10λ, nmaxlog10λ = 2, 200
pnorm = 2.
Klog10λ = transD_GP.GP.Mat32()
λx,λy = 1, 1
x = 0:(0.005λx):λx
y = 0:(0.005λy):0.5λy
λlog10λ = [0.05maximum(x), 0.05maximum(y)]
demean = false
sdev_poslog10λ = [0.01maximum(y), 0.01maximum(x)]
log10bounds = [-1 -0.69; -1 -0.69]
δlog10λ = 0.1
sdev_proplog10λ = [0.1, 0.1]
xall = zeros(2,length(x)*length(y))
for i in 1:size(xall,2)
    xid, yid = Tuple(CartesianIndices((length(x),length(y)))[i])
    xall[:,i] = [x[xid]; y[yid]]
end
xbounds = zeros(Float64,size(xall, 1), 2)
for dim in 1:size(xall, 1)
    xbounds[dim,:] = [minimum(xall[dim,:]), maximum(xall[dim,:])]
end
## Initialize a model using these options
Random.seed!(12)
optlog10λ = transD_GP.OptionsStat(nmin = nminlog10λ,
                        nmax = nmaxlog10λ,
                        xbounds = xbounds,
                        fbounds = log10bounds,
                        xall = xall,
                        λ = λlog10λ,
                        δ = δlog10λ,
                        demean = demean,
                        sdev_prop = sdev_proplog10λ,
                        sdev_pos = sdev_poslog10λ,
                        pnorm = pnorm,
                        quasimultid = false,
                        K = Klog10λ,
                        needλ²fromlog = true,
                        updatenonstat = true,
                        timesλ = 3.6
                        )
@time  log10λ = transD_GP.init(optlog10λ)
## make options for the nonstationary GP -- changes to these should not
# affect the stationary GP
nmin, nmax = 2, 200
fbounds = [-2. 2]
δ = 0.1
sdev_prop = [0.1]
sdev_pos = [0.05, 0.05]
K = transD_GP.GP.Mat32()
demean_ns = false
sampledc_ns = true
## Initialize model for the nonstationary GP
Random.seed!(13)
opt = transD_GP.OptionsNonstat(optlog10λ,
                        nmin = nmin,
                        nmax = nmax,
                        fbounds = fbounds,
                        δ = δ,
                        demean = demean_ns,
                        sampledc = sampledc_ns,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        K = K
                        )
@time m = transD_GP.init(opt, log10λ)
##
fracthresh = 0.05
lsc = transD_GP.testupdate(optlog10λ, log10λ)
ftest = transD_GP.testupdate(opt, log10λ, m)
@testset "Nonstationary GP model changes" begin
    @testset "init test" begin
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        for i = 1:150
            transD_GP.birth!(log10λ, optlog10λ, m, opt)
        end
        lsc[:] = transD_GP.testupdate(optlog10λ, log10λ)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
    end
    @testset "birth tests" begin
        for i = 1:100
            transD_GP.birth!(m, opt, log10λ)
        end
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
    end
    @testset "death tests" begin
        for i = 1:100
            transD_GP.death!(m, opt)
        end
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
    end
    @testset "undo birth" begin
        mold = deepcopy(m)
        transD_GP.birth!(m, opt, log10λ)
        transD_GP.undo_birth!(m)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
        @test norm(mean(mold.fstar - m.fstar)) < 1e-12
    end
    @testset "undo multiple birth" begin
        transD_GP.birth!(m, opt, log10λ )
        transD_GP.birth!(m, opt, log10λ)
        mold = deepcopy(m)
        transD_GP.birth!(m, opt, log10λ)
        transD_GP.undo_birth!(m)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
        @test norm(mean(mold.fstar - m.fstar)) < 1e-12
    end
    @testset "undo death" begin
        mold = deepcopy(m)
        transD_GP.death!(m, opt)
        transD_GP.undo_death!(m, opt)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
        @test norm(mean(mold.fstar - m.fstar)) < 1e-12
    end
    @testset "undo multiple death" begin
        transD_GP.birth!(m, opt, log10λ)
        transD_GP.birth!(m, opt, log10λ)
        transD_GP.birth!(m, opt, log10λ)
        mold = deepcopy(m)
        transD_GP.death!(m, opt)
        transD_GP.undo_death!(m, opt)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
        @test norm(mean(mold.fstar - m.fstar)) < 1e-12
    end
    @testset "property change" begin
        transD_GP.property_change!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
    end
    @testset "undo property change" begin
        mold = deepcopy(m)
        transD_GP.property_change!(m, opt)
        transD_GP.undo_property_change!(m)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
        @test norm(mean(mold.fstar - m.fstar)) < 1e-12
    end
    @testset "position change" begin
        transD_GP.position_change!(m, opt, log10λ)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
    end
    @testset "undo position change" begin
        mold = deepcopy(m)
        transD_GP.position_change!(m, opt, log10λ)
        transD_GP.undo_position_change!(m, opt, log10λ)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
        @test norm(mean(mold.fstar - m.fstar)) < 1e-12
    end
    @testset "dc change" begin
        transD_GP.dc_change!(m, opt)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
    end
    @testset "undo dc change" begin
        mold = deepcopy(m)
        transD_GP.dc_change!(m, opt)
        transD_GP.undo_dc_change!(m)
        transD_GP.sync_model!(m, opt)
        ftest = transD_GP.testupdate(opt, log10λ, m)
        @test norm(mean(lsc' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test mean(abs.((ftest - m.fstar)./ftest)) < fracthresh
        @test norm(mean(mold.fstar - m.fstar)) < 1e-12
    end
end
# ## timing for ns birth in ns model
# log10λ = transD_GP.init(optlog10λ)
# m = transD_GP.init(opt, log10λ)
# for i = 1:148
#     transD_GP.birth!(log10λ, optlog10λ, m, opt)
# end
# @time for i = 1:148
#     transD_GP.birth!(m, opt, log10λ)
# end
# ## time for ntimes birth death
# NTIMES = 20
# ntimes = 150
# T = zeros(NTIMES)
# for I = 1:NTIMES
#     T[I] = time()
#     for i = 1:ntimes
#         transD_GP.birth!(m, opt, log10λ)
#         transD_GP.death!(m, opt)
#     end
#     T[I] = (time() - T[I])/2ntimes
# end
# @info "time for $ntimes birth/death is $(mean(T)) +- $(std(T)/sqrt(NTIMES))"
# ## plot
# l = log10.(log10λ.fstar.^0.5)
# fig,ax = plt.subplots(1,3, sharex=true, sharey=true)
# vmin, vmax = minimum(l[1,:]), maximum(l[1,:])
# im1 = ax[1].imshow(reshape(l[1,:],length(x), length(y)), extent=[y[1],y[end],x[end],x[1]],
#     vmin=vmin, vmax=vmax)
# colorbar(im1, ax=ax[1])
# ax[1].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], marker="+",c="r")
# ax[1].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], c=log10λ.ftrain[1,1:log10λ.n], alpha=0.8)
# vmin, vmax = minimum(l[2,:]), maximum(l[2,:])
# im2 = ax[2].imshow(reshape(l[2,:],length(x), length(y)), extent=[y[1],y[end],x[end],x[1]],
#     vmin=vmin, vmax=vmax)
# colorbar(im2, ax=ax[2])
# ax[2].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], marker="+",c="r")
# ax[2].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], c=log10λ.ftrain[2,1:log10λ.n], alpha=0.8,
#     vmin=vmin, vmax=vmax)
# vmin, vmax = minimum(m.fstar), maximum(m.fstar)
# im3 = ax[3].imshow(reshape(m.fstar,length(x), length(y)), extent=[y[1],y[end],x[end],x[1]],
#     vmin=vmin, vmax=vmax)
# colorbar(im3, ax=ax[3])
# ax[3].scatter(m.xtrain[2,1:m.n], m.xtrain[1,1:m.n], marker="+",c="r")
# ax[3].scatter(m.xtrain[2,1:m.n], m.xtrain[1,1:m.n], c=m.ftrain[1:m.n], alpha=0.8,
#         vmin=vmin, vmax=vmax)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
"""Sample"""
module stats
abstract type Distribution end
abstract type LatheObject end
export Distribution
export LatheObject
#==
Base
==#
include("stats/statbase.jl")
export mean, median, mode, variance, confints, ste, std, q1
export q3, getranks, fact, is_prime, Σ
#==
Distributions
==#
include("stats/distributions.jl")
export NormalDist, T_Dist
#==
Inferential
==#
include("stats/inferential.jl")
export correlationcoeff, TwoTailed, OneTailed
#==
Bayesian
==#
include("stats/bayesian.jl")
export bay_ther, cond_prob
#==
Validate
==#
include("stats/validate.jl")
export mae, mse, r2, catacc
#==
Samples
==#
include("stats/samples.jl")
export sample, sample!
#---------------------------
end


text:
using TensorKit
using TensorOperations
using TensorKitAD
using Zygote
using KrylovKit
using LinearAlgebra
using Plots
using Optim

using Revise
using statmech_tm_solver

T = cmpo_xz()
ψ = cmps(T.Q, T.R)

beta = 16
Fexact = -0.32036477338588226
conv_meas = 1e-3
fidel=-1

log_ovlp(ψ, ψ, beta)

for chi in [36]
    global conv_meas, fidel, conv_meas, ψ
    @show chi
    for ix in 1:100
        Tψ = act(T, ψ)
        status, fidel, ψ1 = compress(Tψ, chi, beta; tol=conv_meas)
        conv_meas1 = convergence_measure(T, ψ1, beta)
        @show fidel, conv_meas
        if status
            ψ = ψ1
            conv_meas = conv_meas1
            f = -(log_ovlp(ψ, act(T, ψ), beta) - log_ovlp(ψ, ψ, beta)) / beta |> real
            @show f, (f - Fexact) / Fexact, conv_meas
        else
            break
        end
    end
    f = -(log_ovlp(ψ, act(T, ψ), beta) - log_ovlp(ψ, ψ, beta)) / beta |> real
    @show chi, get_chi(ψ), f, (f - Fexact) / Fexact
    @show conv_meas

    @show "self-iterating..."
    for ix in 1:50
        print("$(ix) \r")
        Tψ = act(T, ψ)
        status, fidel, ψ1 = compress(Tψ, chi, beta; tol=conv_meas, init=ψ)
        conv_meas1 = convergence_measure(T, ψ1, beta)
        if status && fidel > -conv_meas1 * 0.99
            ψ = ψ1
            conv_meas = conv_meas1
            f = -(log_ovlp(ψ, act(T, ψ), beta) - log_ovlp(ψ, ψ, beta)) / beta |> real
            @show f, (f - Fexact) / Fexact, fidel, conv_meas
        else
            break
        end
        f = -(log_ovlp(ψ, act(T, ψ), beta) - log_ovlp(ψ, ψ, beta)) / beta |> real
    end

    quicksave("ckpt_cMPO_xy_beta$(beta)_chi$(chi)", convert_to_tensormap(ψ))
end


text:
function test_pseudohnf()

  set_assert_level(:NfMaxOrdQuoRing, 1)

  QQ = FlintRationalField()
  Qx, x = PolynomialRing(QQ, "x")

  # Compute a pseudo-hnf of a matrix over Z and check result against the HNF

  K, a = NumberField(x - 1, "a")
  O = maximal_order(K)

  A =
  [ 154830247 724968499 442290149 700116640 313234502;
    384015369 193254623 792471850 452111534 717477652;
    816720700 46054709 339475042 389090910 103149313;
    104945689 475264799 821863415 806964485 676397437;
    281047486 709449227 590950934 18224679 295696061]

  Ahnf =
  [ 1 0 0 0 34095268312495604433062164807036095571620397;
    0 1 0 0 18160914505139254049490851505162505507397915;
    0 0 1 0 37265321283634252189185025806030875371477390;
    0 0 0 1 2276874339297612770861218322365243729516503
    0 0 0 0 37684868701591492337245802520684209569420259]

  de = fmpz(37684868701591492337245802520684209569420259)
  AoverO = MatrixSpace(O, 5, 5)(map(z -> O(z), A))

  Apm = Hecke.PseudoMatrix( AoverO, [(O(1)*O)::Hecke.NfMaxOrdIdl for i in 1:5])

  d = num(det(Apm))

  Apseudohnf = Hecke.pseudo_hnf(Apm, d)

  z = Apseudohnf.matrix

  for i in 1:rows(z)
    Hecke.mul_row!(z, i, K(norm(Apseudohnf.coeffs[i])))
  end

  zinZ = MatrixSpace(FlintZZ, 5, 5)(map(zz -> num(elem_in_basis(O(zz))[1]), z.entries))
  c = parent(zinZ)(Ahnf) - zinZ

  @test all([ mod(c[i,j], de) == 0 for i in 1:5, j in 1:5])

  # Construct random pseudo-matrices over different fields and check if they
  # span the same module

  for i in 2:5
    K, a = NumberField(x^i - 10)
    O = maximal_order(K)
    println("  Testing over field $(x^i - 10)")

    for j in 1:1
      l = rand(10:20) - i + 1
      ll = rand(1:20)
      z = rand(MatrixSpace(O, l, l), fmpz(2)^ll)
      println("    $l x $l matrix with $ll bits")
      cc = [ (O(1)*O)::Hecke.NfMaxOrdIdl for i in 1:l]
      pm = Hecke.PseudoMatrix(z, cc)
      d = det(pm)
      ppm = Hecke.pseudo_hnf(pm, num(d))
      @test Hecke._spans_subset_of_pseudohnf(pm ,ppm)
      @test d == det(ppm)
    end
  end
end


text:
doc"""
  unlike the other files, generate_posit_class_h creates a separate file for
  each header file.
"""
function generate_posit_class_h(io, n, es)
  write(io, """
  #ifndef POSIT_$(n)_$(es)_H
  #define POSIT_$(n)_$(es)_H

  #include "posit.h"

  class $(c(n,es)){
    public:
      uint$(n)_t data;

      //various public constructors
      $(c(n,es))();             //defaults to zero
      $(c(n,es))(const float);        //conversion constructor from IEEE32
      $(c(n,es))(const double);       //conversion constructor from IEEE64
      $(c(n,es))(const $(c(n,es)) &); //copy constructor
      $(c(n,es))(const $(p(n,es)));   //bridge constructor from c functionality.

      //public operators
      $(c(n,es)) operator -() const;
      $(c(n,es)) operator +(const $(c(n,es)) rhs) const;
      $(c(n,es)) operator -(const $(c(n,es)) rhs) const;
      $(c(n,es)) operator *(const $(c(n,es)) rhs) const;
      $(c(n,es)) &operator +=(const $(c(n,es)) rhs);
      $(c(n,es)) &operator -=(const $(c(n,es)) rhs);
      $(c(n,es)) &operator *=(const $(c(n,es)) rhs);
      bool operator <(const $(c(n,es)) rhs) const;
      bool operator <=(const $(c(n,es)) rhs) const;
      bool operator >(const $(c(n,es)) rhs) const;
      bool operator >=(const $(c(n,es)) rhs) const;
      bool operator ==(const $(c(n,es)) rhs) const;
      $(c(n,es)) operator /(const $(c(n,es)) rhs) const;
      $(c(n,es)) &operator /=(const $(c(n,es)) rhs);

      operator float() const;
      operator double() const;
      operator $(p(n,es))() const;
  };

  bool isequal(const $(c(n,es)) a, const $(c(n,es)) b);
  $(c(n,es)) mulinv(const $(c(n,es)) arg);
  $(c(n,es)) log2(const $(c(n,es)) arg);
  $(c(n,es)) exp2(const $(c(n,es)) arg);
  $(c(n,es)) fma(const $(c(n,es)) a, const $(c(n,es)) b, const $(c(n,es)) c);
  $(c(n,es)) fms(const $(c(n,es)) a, const $(c(n,es)) b, const $(c(n,es)) c);
  $(c(n,es)) nfma(const $(c(n,es)) a, const $(c(n,es)) b, const $(c(n,es)) c);
  $(c(n,es)) nfms(const $(c(n,es)) a, const $(c(n,es)) b, const $(c(n,es)) c);
  $(c(n,es)) sqrt(const $(c(n,es)) arg);
  $(c(n,es)) log1p(const $(c(n,es)) arg);
  $(c(n,es)) log(const $(c(n,es)) arg);
  $(c(n,es)) log10(const $(c(n,es)) arg);
  $(c(n,es)) exp(const $(c(n,es)) arg);
  $(c(n,es)) pow(const $(c(n,es)) lhs, const $(c(n,es)) rhs);
  $(c(n,es)) sin(const $(c(n,es)) arg);
  $(c(n,es)) cos(const $(c(n,es)) arg);
  $(c(n,es)) atan(const $(c(n,es)) arg);
  $(c(n,es)) atan2(const $(c(n,es)) y, const $(c(n,es)) x);

  #endif

  """)
end


text:
function distance_matrix_float(x, y)
    @assert length(x) == length(y)
    n = length(x)
    M = zeros(n, n)
    for i in 1:n
        for j in i+1:n
            if i != j
                M[i, j] = sqrt( (x[i]-x[j])^2 + (y[i]-y[j])^2 )
                M[j, i] = M[i, j]
            end
        end
    end
    return M
end

function distance_matrix(x, y)
    M = distance_matrix_float(x, y)
    return Int.(round.(M))
end

function tour_length(tour, M::Matrix{T}) where T
    sum = zero(T)
    for i in 1:length(tour)
        if i < length(tour)
            sum += M[tour[i], tour[i+1]]
        else
            sum += M[tour[i], tour[1]]
        end
    end
    return sum
end

# function tour_length(tour, x, y)
#     M = distance_matrix(x, y)
#     return tour_length(tour, M)
# end

function solve_TSP_Concorde(x, y; norm="EUC_2D",)
    Concorde = pyimport("concorde.tsp")
    solver = Concorde.TSPSolver.from_data(x, y, norm)
    tour_data = solver.solve()
    tour_concorde = tour_data[1] .+ 1  # 0-index -> 1-index
    # length_concorde = tour_length(tour_concorde, x, y)
    length_concorde = Int(tour_data[2])
    return tour_concorde, length_concorde
end

function solve_TSP_LKH(M::Matrix{T}) where T
    Elkai = pyimport("elkai")
    if T == Int
        tour_LKH = Elkai.solve_int_matrix(M) .+ 1 #adding 1 for indexing
    else
        tour_LKH = Elkai.solve_float_matrix(M) .+ 1 #adding 1 for indexing
    end
    length_LKH = tour_length(tour_LKH, M)
    return tour_LKH, length_LKH
end

function solve_TSP_LKH(x, y)
    M = distance_matrix(x, y)
    return solve_TSP_LKH(M)
end


text:
include("potentially_intersecting_simplices.jl")
include("reshape_simplices.jl")
include("function_barriers.jl")

import ..AbstractTransferOperator
import ..TransferOperatorTriangulationApprox

import StateSpaceReconstruction:
    DelaunayTriangulation,
    Simplex,
	point_contained
import StateSpaceReconstruction: subsample_coeffs
import StaticArrays: Size, MMatrix
import Statistics: det
import InplaceOps: @!

#####################################
# Exports
#####################################
export transferoperator_triangulation_approx

function mdet_fm_in_x_view!(M, v, n, starts, stops)
    @inbounds for i = 0:(n - 1)
        M[starts[i+1]:stops[i+1]] = view(v, starts[i+1]:stops[i+1])
    end
    @fastmath det(M)
end

function contained_optim!(signs, s_arr, sx, point, dim, temp_arr, starts, stops)
    # Redefine the temporary simplex. This is in-place, so we don't allocate
    # memory. We could also have re-initialised `signs`, but since we're never
    # comparing more than two consecutive signs, this is not necessary.
    rezero!(s_arr)
    rezero!(signs)
    fill_into!(s_arr, sx)

    #Replace first vertex with the point
    fill_at_inds!(s_arr, point, 1:dim)

    # Signed volume
    signs[1] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

    rezero!(s_arr)
    fill_into!(s_arr, sx) #reset

    for κ = 2:dim # Check remaining signs and stop if sign changes
        # Replace the ith vertex with the point we're cheking (leaving the
        # 1 appended to Vi intact.)
        idxs = ((dim + 1)*(κ - 1)+1):((dim + 1)*(κ - 1)+ 1 + dim - 1)
        fill_at_inds!(s_arr, point, idxs) # ith change

        signs[κ] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

        if !(signs[κ-1] == signs[κ])
            return false
        end

        rezero!(s_arr)
        fill_into!(s_arr, sx)
    end

    # Last the last vertex with the point in question
    idxs = ((dim + 1)*(dim)+1):((dim+1)^2-1)
    fill_at_inds!(s_arr, point, idxs)

    signs[end] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

    if !(signs[end-1] == signs[end])
       return false
    else
        return true
    end
end

function innerloop_optim!(inds::Vector{Int}, signs, s_arr, Sj, pt, dim::Int, M, i::Int, temp_arr, starts, stops)
    for j in 1:length(inds)
        if contained_optim!(signs, s_arr, Sj[j], pt, dim, temp_arr, starts, stops)
            M[inds[j], i] += 1.0
        end
    end
end

"""
    transferoperator_triangulation_approx(invariant_pts,
        DT::DelaunayTriangulation, n_sample_pts::Int = 200, sample_randomly::Bool = false)

Approximate the transfer operator over a triangulation of `invariant_pts`. Does the same
as the single-argument version of the function, but requires that the triangulation
is precomputed on `invariant_pts[1:(end - 1)]` (excluding the last point).

## Example

```julia
# Create some random points and make sure they are invariant
pts = [rand(3) for i = 1:30]
invariant_pts = invariantize(pts)

# Triangulate all but the last point (exclude the last, so that an image simplex
# always exists)
DT = DelaunayTriangulation(invariant_pts[1:(end - 1)])

# Compute the transfer operator
to = transferoperator_triangulation_approx(invariant_pts, DT)
```

## References
The idea of using a triangulation estimator to compute the transfer operator was
introduced in [1], where it is used to compute physical invariant measures. To compute
a physical invariant measure from a transfer operator `to`, simply call
`invariantmeasure(to)`.

1. Froyland, Gary. "Estimating physical invariant measures and space averages of dynamical systems indicators." Bulletin of the Australian Mathematical Society 56.1 (1997): 157-159.

"""
function transferoperator_triangulation_approx(inv_pts,
        DT::DelaunayTriangulation; n_sample_pts::Int = 200, sample_randomly::Bool = false)

    # Some constants used throughout the function
    n_simplices = length(DT)
    dim = length(DT[1]) - 1

    #=
    # Prepare memory-efficient representations of the simplices
    =#
    simplices = reshape_simplices(inv_pts[1:(end- 1)], DT)
    image_simplices = reshape_simplices(inv_pts[2:end], DT)

    #=
    Compute the convex expansions coefficients needed to generate points
    within a simplex. Also, update number of points if a shape-preserving
    simplex subdivision algorithm was used (in that case, we get a
    *minimum* of `n_sample_pts` set of coefficients.
    =#
    convex_coeffs = subsample_coeffs(dim, n_sample_pts, sample_randomly)
    n_coeffs::Int = size(convex_coeffs, 2)
    convex_coeffs = [convex_coeffs[:, i] for i = 1:n_coeffs]

    # Pre-allocated arrays (SizedArrays, for efficiency)
    pt          = Size(1, dim)(zeros(Float64, 1, dim))
    s_arr       = Size((dim+1)^2)(zeros(Float64, (dim+1)^2))
    signs       = Size(dim + 1)(zeros(Float64, dim + 1))

    # Re-arrange simplices so that look-up is a bit more efficient
    simplex_arrs = Vector{Array{Float64, 2}}(undef, n_simplices)
    imsimplex_arrs = Vector{Array{Float64, 2}}(undef, n_simplices)

    for i in 1:n_simplices
        simplex_arrs[i] = hcat(inv_pts[DT[i]]...,)
        imsimplex_arrs[i] = hcat(inv_pts[DT[i] .+ 1]...,)
    end

    # Pre-allocate some information to avoid memory allocation when
    # computing the fractional overlaps. These are a temporary matrix
    # which we use to compute determinants (to check if a point lies
    # inside a simplex), and the indices needed to access its values
    # at different stages of the checking stage (of whether a point
    # lies inside a simplex). Pre-allocating this information and
    # passing it on reduces runtime and memory allocation
    # by half an order of magnitude.

    temp_arr = rand(MMatrix{dim + 1, dim + 1})
    starts = [((dim + 1)*i)+1 for i = 0:dim]
    stops = starts .+ dim

    # The Markov matrix
    M = zeros(Float64, n_simplices, n_simplices)

    for i in 1:n_simplices
        # Find the indices of the simplices who potentially intersect
        # with the image of simplex #i
        idxs = idxs_potentially_intersecting_simplices(inv_pts, DT, i)
        Sj = Vector{AbstractArray}(undef, length(idxs))
        get_simplices_at_inds!(Sj, idxs, simplices)

        # Generate points within the image of simplex #i and check
        # what fraction of those points fall into the simplices.
        # We approximate the intersecting volume between the
        # image simplex and the individual simplices by those
        # fractions.
        @views is = imsimplex_arrs[i]
        for k in 1:n_coeffs
            @! pt = transpose(convex_coeffs[k]) * transpose(is)
            innerloop_optim!(idxs, signs, s_arr, Sj, pt, dim, M, i, temp_arr, starts, stops)
        end
    end
    # Need to normalise, because all we have up until now is counts
    # of how many points inside the image simplex falls into
    # the simplices.
    return TransferOperatorTriangulationApprox(transpose(M) ./ n_coeffs)
end

"""
    transferoperator_triangulation_approx(invariant_pts; n_sample_pts::Int = 200,
        sample_randomly::Bool = false) -> TransferOperatorTriangulationApprox

Compute an approximation of the transfer operator over a triangulated partition of the
provided points (`invariant_pts`). The estimator assumes that the points are invariant
under the forward projection of the simplices in the triangulation, which means that
the last point must be contained within the convex hull of all previous points.

*Tip: The `invariantize` function from `StateSpaceReconstruction` can be used to move the
last point toward the center of the embedding until invariance is achieved.*

This function computes approximate simplex intersections and is slow for very large
triangulations. It should only be used on small to moderately sized datasets (i.e.
in the order of hundreds of points). For larger datasets, use a rectangular estimator.

Returns a  instance.

## Example
```julia
# Create some random points and make sure they are invariant
pts = [rand(3) for i = 1:30]
invariant_pts = invariantize(pts)

# Compute the transfer operator
to = transferoperator_triangulation_approx(invariant_pts)
```

## References
The idea of using a triangulation estimator to compute the transfer operator was
introduced in [1], where it is used to compute physical invariant measures. To compute
a physical invariant measure from a transfer operator `to`, simply call
`invariantmeasure(to)`.

1. Froyland, Gary. "Estimating physical invariant measures and space averages of dynamical systems indicators." Bulletin of the Australian Mathematical Society 56.1 (1997): 157-159.
"""
function transferoperator_triangulation_approx(invariant_pts;
        n_sample_pts::Int = 200, sample_randomly::Bool = false)
    # Triangulate all points but the last (we need the last point
    # for the forward projection of the simplices)
    DT = DelaunayTriangulation(invariant_pts[1:(end - 1)])
    transferoperator_triangulation_approx(invariant_pts, DT)
end


text:
function max_product(num_string)
    greatest = 0
    for i in 1:length(num_string)-12
        num = prod(parse(Int, num_string[j]) for j in i:i+12)
        greatest = max(greatest, num)
    end
    println(greatest)
end

long_num = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"
max_product(long_num)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf
using ImPlot

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
window = GLFW.CreateWindow(1280, 720, "ImPlot Demo")
@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()
pctx = ImPlot.CreateContext()
ImPlot.SetImGuiContext(ctx)
# setup Dear ImGui style
CImGui.StyleColorsDark()

# load Fonts
fonts_dir = joinpath(@__DIR__, "..", "fonts")
fonts = CImGui.GetIO().Fonts
CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

try
    show_demo_window = true
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()

        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        # show the big demo window
        show_demo_window && @c ImPlot.LibCImPlot.ShowDemoWindow(&show_demo_window)

        # show a simple window that we create ourselves.
        # we use a Begin/End pair to created a named window.
        @cstatic f = Cfloat(0.0) counter = Cint(0) begin
            CImGui.Begin("Hello, world!")
            @c CImGui.Checkbox("Show ImPlot Demo", &show_demo_window)
            CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)",
                                 1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))

            CImGui.End()
        end

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end
catch e
    @error "Error in renderloop!" exception = e
    Base.show_backtrace(stderr, catch_backtrace())
finally
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    ImPlot.DestroyContext(pctx)
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end


text:
# This file implements Tables interface and provide compatibility
# to the Queryverse ecosystem.

# -----------------------------------------------------------------------------
# Tables.jl implementation

Tables.istable(::Type{<:ResultSet}) = true

# AbstractColumns interface
Tables.columnaccess(::Type{<:ResultSet}) = true
Tables.columns(rs::ResultSet) = rs
Tables.getcolumn(rs::ResultSet, i::Int) = rs[names(rs)[i]]

# AbstractRow interface
Tables.rowaccess(::Type{<:ResultSet}) = true
Tables.rows(rs::ResultSet) = rs

# Schema definition
Tables.schema(rs::ResultSet) = Tables.Schema(names(rs), eltype.(columns(rs)))

# -----------------------------------------------------------------------------
# Queryverse compatibility

IteratorInterfaceExtensions.getiterator(rs::ResultSet) = Tables.datavaluerows(rs)
TableTraits.isiterabletable(x::ResultSet) = true


text:
using JLD2
using Statistics

function lg_start!(self::E, env, agent) where {E <: ErrorRecorder}
    Q_est = hcat([get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"], self.eval_set["actions"])]...)
    true_values = self.get_true_values(env, self.eval_set["ests"])
    err = mean((Q_est - true_values) .^ 2, dims=2)
    self.error[:,1] = sqrt.(err)
end

function lg_step!(self::E, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total) where {E <: ErrorRecorder}
    if (rem(cur_step_total, self.log_interval) == 0)
        ind = fld(cur_step_total, self.log_interval) + 1
        Q_est = hcat([get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"], self.eval_set["actions"])]...)

        true_values = self.get_true_values(env, self.eval_set["ests"])
        err = mean((Q_est - true_values) .^ 2, dims=2)

        self.error[:,ind] = sqrt.(err)
    end
end

function lg_episode_end!(self::E, cur_step_in_episode, cur_step_total) where {E<:ErrorRecorder}
end

function save_log(self::E, save_dict::Dict) where {E<:ErrorRecorder}
    save_dict[self.save_key] = self.error
end

####
# Error Recording for 1D Tmaze
####
mutable struct OneDTMazeError_dmu <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function OneDTMazeError_dmu(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/OneDTMazeEvalSet_dmu.jld2") OneDTMazeEvalSetDMU
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4,num_logged_steps), OneDTMazeEvalSetDMU, logger_init_info[LoggerInitKey.INTERVAL], :oned_tmaze_dmu_error, OneDTMazeUtils.get_true_values)
    end
end

mutable struct OneDTMazeError <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function OneDTMazeError(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/OneDTMazeEvalSet.jld2") OneDTMazeEvalSet
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4,num_logged_steps), OneDTMazeEvalSet, logger_init_info[LoggerInitKey.INTERVAL], :oned_tmaze_old_error, OneDTMazeUtils.get_true_values)
    end
end

mutable struct OneDTMazeError_Uniform <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function OneDTMazeError_Uniform(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/OneDTMazeEvalSet_Uniform.jld2") OneDTMazeEvalSetUniform
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4,num_logged_steps), OneDTMazeEvalSetUniform, logger_init_info[LoggerInitKey.INTERVAL], :oned_tmaze_uniform_error, OneDTMazeUtils.get_true_values)
    end
end

mutable struct OneDTMazeError_DPI <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function OneDTMazeError_DPI(logger_init_info)
        @load string(pwd(),"/src/data/OneDTMazeEvalSet_d_pi.jld2") eval_set
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4, num_logged_steps), eval_set, logger_init_info[LoggerInitKey.INTERVAL], :oned_tmaze_dpi_error, OneDTMazeUtils.get_true_values)
    end
end

function lg_start!(self::OneDTMazeError_DPI, env, agent)

    for gvf_idx ∈ 1:4
        Q_est = [get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"][gvf_idx], self.eval_set["actions"][gvf_idx])]
        true_values = self.get_true_values(env, self.eval_set["ests"][gvf_idx], gvf_idx)
        err = mean((getindex.(Q_est, gvf_idx) - true_values) .^ 2)
        self.error[gvf_idx, 1] = sqrt.(err)
    end

end

function lg_step!(self::OneDTMazeError_DPI, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if (rem(cur_step_total, self.log_interval) == 0)

        ind = fld(cur_step_total, self.log_interval) + 1
        for gvf_idx ∈ 1:4
            Q_est = [get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"][gvf_idx], self.eval_set["actions"][gvf_idx])]
            true_values = self.get_true_values(env, self.eval_set["ests"][gvf_idx], gvf_idx)
            err = mean((getindex.(Q_est, gvf_idx) - true_values) .^ 2)
            self.error[gvf_idx, ind] = sqrt.(err)
        end
    end
end

####
# Error Recording for Tabular Tmaze
####

mutable struct TTMazeError <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function TTMazeError(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/TTMazeEvalSet.jld2") TTMazeEvalSet
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4,num_logged_steps), TTMazeEvalSet, logger_init_info[LoggerInitKey.INTERVAL],:ttmaze_error, TabularTMazeUtils.get_true_values)
    end
end

mutable struct TTMazeUniformError <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function TTMazeUniformError(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/TTMazeUniformEvalSet.jld2") TTMazeUniformEvalSet
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1

        new(zeros(4,num_logged_steps), TTMazeUniformEvalSet, logger_init_info[LoggerInitKey.INTERVAL], :ttmaze_uniform_error, TabularTMazeUtils.get_true_values)
    end
end

mutable struct TTMazeRoundRobinError <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function TTMazeRoundRobinError(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/TTMazeRoundRobinEvalSet.jld2") TTMazeRoundRobinEvalSet
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4,num_logged_steps), TTMazeRoundRobinEvalSet, logger_init_info[LoggerInitKey.INTERVAL],:ttmaze_round_robin_error, TabularTMazeUtils.get_true_values)
    end
end

####
# Error Recording for TwoD
####

mutable struct TwoDGridWorldError <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function TwoDGridWorldError(logger_init_info)
        @load string(pwd(),"/src/data/TwoDGridWorldSet.jld2") eval_set
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4, num_logged_steps), eval_set, logger_init_info[LoggerInitKey.INTERVAL], :twod_grid_world_error, TwoDGridWorldUtils.get_true_values)
    end
end

mutable struct TwoDGridWorldErrorCenterDPI <: ErrorRecorder
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int
    save_key::Symbol
    get_true_values::Function
    function TwoDGridWorldErrorCenterDPI(logger_init_info)
        @load string(pwd(),"/src/data/TwoDGridWorldSet_dpi.jld2") eval_set
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(4, num_logged_steps), eval_set, logger_init_info[LoggerInitKey.INTERVAL], :twod_grid_world_error_center_dpi, TwoDGridWorldUtils.get_true_values)
    end
end

function lg_start!(self::TwoDGridWorldErrorCenterDPI, env, agent)

    for gvf_idx ∈ 1:4
        Q_est = [get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"][gvf_idx], self.eval_set["actions"][gvf_idx])]
        true_values = self.get_true_values(env, self.eval_set["ests"][gvf_idx], gvf_idx)
        err = mean((getindex.(Q_est, gvf_idx) - true_values) .^ 2)
        self.error[gvf_idx, 1] = sqrt.(err)
    end

end

function lg_step!(self::TwoDGridWorldErrorCenterDPI, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if (rem(cur_step_total, self.log_interval) == 0)

        ind = fld(cur_step_total, self.log_interval) + 1
        for gvf_idx ∈ 1:4
            Q_est = [get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"][gvf_idx], self.eval_set["actions"][gvf_idx])]
            true_values = self.get_true_values(env, self.eval_set["ests"][gvf_idx], gvf_idx)
            err = mean((getindex.(Q_est, gvf_idx) - true_values) .^ 2)
            self.error[gvf_idx, ind] = sqrt.(err)
        end
    end
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function hipStreamSynchronize(stream::hipStream_t)
    @check ccall((:hipStreamSynchronize, "libhip_hcc.so"), hipError_t, (hipStream_t,), stream)
end

function hipStreamQuery(stream::hipStream_t)
    @check ccall((:hipStreamQuery, "libhip_hcc.so"), hipError_t, (hipStream_t,), stream)
end

function hipDeviceSynchronize()
    @check ccall((:hipDeviceSynchronize, "libhip_hcc.so"), hipError_t, ())
end


text:
using NeidArchive
using Documenter

DocMeta.setdocmeta!(NeidArchive, :DocTestSetup, :(using NeidArchive); recursive=true)

makedocs(;
    modules=[NeidArchive],
    authors="Eric Ford",
    repo="https://github.com/RvSpectML/NeidArchive.jl/blob/{commit}{path}#{line}",
    sitename="NeidArchive.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://RvSpectML.github.io/NeidArchive.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/RvSpectML/NeidArchive.jl",
    devbranch="main",
)


text:
struct SparseArray{T,N} <: AbstractArray{T,N}
    data::Dict{NTuple{N,Int64}, T}
    dims::NTuple{N,Int64}
    function SparseArray{T,N}(::UndefInitializer, dims::NTuple{N,Int}) where {T,N}
        data = Dict{NTuple{N,Int64}, T}()
        return new{T,N}(data, dims)
    end
    function SparseArray{T}(::UndefInitializer, dims::NTuple{N,Int}) where {T,N}
        data = Dict{NTuple{N,Int64}, T}()
        return new{T,N}(data, dims)
    end
    function SparseArray(A::SparseArray{T,N}) where {T,N}
        new{T,N}(copy(A.data), A.dims)
    end
end

memsize(A::SparseArray) = memsize(A.data)
@inline function Base.getindex(A::SparseArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    return get(A.data, I, zero(T))
end
@inline function Base.setindex!(A::SparseArray{T,N}, v, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    if v != zero(v)
        A.data[I] = v
    else
        delete!(A.data, I) # does not do anything if there was no key I
    end
    return v
end

Array(a::SparseArray{T,N}) where {T,N} = Array{T,N}(a)
function Array{T,N}(a::SparseArray) where {T,N}
    d = fill(zero(T), size(a))
    for (k,v) in a.data
        d[k...] = v
    end
    d
end

SparseArray(a::AbstractArray{T,N}) where {T,N} = SparseArray{T,N}(a)
function SparseArray{T,N}(a::AbstractArray) where {T,N}
    d = SparseArray{T}(undef, size(a))
    for I in CartesianIndices(a)
        a[I] == zero(T) && continue
        d[I] = a[I]
    end
    return d
end

Base.copy(A::SparseArray) = SparseArray(A)

Base.size(A::SparseArray) = A.dims

Base.similar(A::SparseArray, ::Type{S}, dims::Dims{N}) where {S,N} =
    SparseArray{S,N}(undef, dims)

# TODO: Basic arithmitic

# Vector space functions
#------------------------
function LinearAlgebra.lmul!(a::Number, d::SparseArray)
    lmul!(a, d.data.vals)
    # typical occupation in a dict is about 30% from experimental testing
    # the benefits of scaling all values (e.g. SIMD) largely outweight the extra work
    return d
end
function LinearAlgebra.rmul!(d::SparseArray, a::Number)
    rmul!(d.data.vals, a)
    return d
end
function LinearAlgebra.axpby!(α, x::SparseArray, β, y::SparseArray)
    lmul!(y, β)
    for (k, v) in x
        y[k] += α*v
    end
    return y
end
function LinearAlgebra.axpy!(α, x::SparseArray, y::SparseArray)
    for (k, v) in x
        y[k] += α*v
    end
    return y
end

function LinearAlgebra.norm(x::SparseArray, p::Real = 2)
    norm(Base.Generator(last, A.data), p)
end

function LinearAlgebra.dot(x::SparseArray, y::SparseArray)
    size(x) == size(y) || throw(DimensionMismatch("dot arguments have different size"))
    s = dot(zero(eltype(x)), zero(eltype(y)))
    if length(x.data) >= length(y.data)
        iter = keys(x.data)
    else
        iter = keys(y.data)
    end
    @inbounds for I in iter
        s += dot(x[I...], y[I...])
    end
    return s
end

# TensorOperations compatiblity
#-------------------------------
function add!(α, A::SparseArray{<:Any, N}, CA::Symbol,
                β, C::SparseArray{<:Any, N}, indCinA) where {N}

    (N == length(indCinA) && TupleTools.isperm(indCinA)) ||
        throw(IndexError("Invalid permutation of length $N: $indCinA"))
    size(C) == TupleTools.getindices(size(A), indCinA) ||
        throw(DimensionMismatch("non-matching sizes while adding arrays"))

    β == one(β) || LinearAlgebra.lmul!(β, C);
    for (kA, vA) in A.data
        kC = TupleTools.getindices(kA, indCinA)
        C[kC...] += α* (conjA == :C ? conj(vA) : vA)
    end
    C
end

function trace!(α, A::SparseArray{<:Any, NA}, CA::Symbol, β, C::SparseArray{<:Any, NC},
        indCinA, cindA1, cindA2) where {NA,NC}

    NC == length(indCinA) ||
        throw(IndexError("Invalid selection of $NC out of $NA: $indCinA"))
    NA-NC == 2*length(cindA1) == 2*length(cindA2) ||
        throw(IndexError("invalid number of trace dimension"))
    pA = (indCinA..., cindA1..., cindA2...)
    TupleTools.isperm(pA) ||
        throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))

    sizeA = size(A)
    sizeC = size(C)

    TupleTools.getindices(sizeA, cindA1) == TupleTools.getindices(sizeA, cindA2) ||
        throw(DimensionMismatch("non-matching trace sizes"))
    sizeC == TupleTools.getindices(sizeA, indCinA) ||
        throw(DimensionMismatch("non-matching sizes"))

    β == one(β) || LinearAlgebra.lmul!(β, C);

    for (kA, v) in A.data
        kAc1 = TupleTools.getindices(kA, cindA1)
        kAc2 = TupleTools.getindices(kA, cindA2)
        kAc1 == kAc2 || continue

        kC = TupleTools.getindices(kC, indCinA)
        C[kC...] += α * (conjA == :C ? conj(v) : v)
    end
    return C
end

function contract!(α, A::SparseArray, CA::Symbol, B::SparseArray, CB::Symbol,
        β, C::SparseArray,
        oindA::IndexTuple, cindA::IndexTuple, oindB::IndexTuple, cindB::IndexTuple,
        indCinoAB::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    pA = (oindA...,cindA...)
    (length(pA) == ndims(A) && TupleTools.isperm(pA)) ||
        throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))
    pB = (oindB...,cindB...)
    (length(pB) == ndims(B) && TupleTools.isperm(pB)) ||
        throw(IndexError("invalid permutation of length $(ndims(B)): $pB"))
    (length(oindA) + length(oindB) == ndims(C)) ||
        throw(IndexError("non-matching output indices in contraction"))
    (ndims(C) == length(indCinoAB) && isperm(indCinoAB)) ||
        throw(IndexError("invalid permutation of length $(ndims(C)): $indCinoAB"))

    sizeA = size(A)
    sizeB = size(B)
    sizeC = size(C)

    csizeA = TupleTools.getindices(sizeA, cindA)
    csizeB = TupleTools.getindices(sizeB, cindB)
    osizeA = TupleTools.getindices(sizeA, oindA)
    osizeB = TupleTools.getindices(sizeB, oindB)

    csizeA == csizeB ||
        throw(DimensionMismatch("non-matching sizes in contracted dimensions"))
    TupleTools.getindices((osizeA..., osizeB...), indCinoAB) == size(C) ||
        throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))

    β == one(β) || LinearAlgebra.lmul!(β, C);

    for (kA, vA) in A.data
        kAc = TupleTools.getindices(kA, cindA)
        kAo = TupleTools.getindices(kA, oindA)
        for (kB, vB) in B.data
            kBc = TupleTools.getindices(kB, cindB)
            kAc == kBc || continue

            kBo = TupleTools.getindices(kB, oindB)

            kABo = (kAo..., kBo...)

            kC = TupleTools.getindices(kABo, indCinoAB)

            C[kC...] += α * (CA == :C ? conj(vA) : vA) * (CB == :C ? conj(vB) : vB)
        end
    end
    C
end


text:
"""
This file is for Gym cartpole environment https://github.com/openai/gym/tree/master/gym/envs/classic_control.
These are based on https://github.com/WilsonWangTHU/mbbl
"""

using LyceumBase
using Shapes
using Base: @propagate_inbounds

struct CartpoleGym{D,SP,OP,AP,RP} <: AbstractEnvironment
    state::D
    obs::D
    action::D
    time::D
    nextobs::D
    ctrlrange::D

    sp::SP
    op::OP
    ap::AP
    rp::RP
    function CartpoleGym{T}() where T<:AbstractFloat
        nobs = 4
        naction = 1

        ctrlrange = zeros(T, 2)
        ctrlrange[1] = -1.
        ctrlrange[2] = 1.

        sp = VectorShape(T, nobs)
        op = VectorShape(T, nobs)
        ap = VectorShape(T, naction)
        rp = ScalarShape(T)
        s = zeros(T, nobs)
        o = zeros(T, nobs)
        a = zeros(T, naction)
        env = new{typeof(s),
        typeof(sp), typeof(op), typeof(ap), typeof(rp)}(s,o,a,zeros(T,1),
                                                        zeros(T, nobs),
                                                        ctrlrange,
                                                        sp,op,ap,rp)
        reset!(env)
    end
end

function LyceumBase.tconstruct(::Type{CartpoleGym}, n::Integer)
    Tuple(CartpoleGym{Float64}() for _ in 1:n)
end
CartpoleGym() = first(tconstruct(CartpoleGym, 1))

@inline LyceumBase.statespace(env::CartpoleGym) = env.sp
@inline LyceumBase.obsspace(env::CartpoleGym) = env.op
@inline LyceumBase.actionspace(env::CartpoleGym) = env.ap
@inline LyceumBase.rewardspace(env::CartpoleGym) = env.rp

@inline function LyceumBase.getstate!(state, env::CartpoleGym)
    #state .= env.state
    state[1] = env.state[1] # unroll for speed
    state[2] = env.state[2]
    state[3] = env.state[3]
    state[4] = env.state[4]
end
@inline function LyceumBase.setstate!(env::CartpoleGym, state)
    #env.state .= state
    env.state[1] = state[1] # unroll for speed
    env.state[2] = state[2]
    env.state[3] = state[3]
    env.state[4] = state[4]
end

@inline LyceumBase.getaction!(action, env::CartpoleGym) = action .= env.action
@inline LyceumBase.setaction!(env::CartpoleGym, action) = env.action .= action

@propagate_inbounds function LyceumBase.isdone(state, action, obs, env::CartpoleGym)
    false
end
@propagate_inbounds function LyceumBase.getobs!(obs, env::CartpoleGym)
    getstate!(env.obs, env) # abuse of functions but valid
    getstate!(obs, env)
end
@propagate_inbounds function LyceumBase.getreward(s, a, o, env::CartpoleGym)
    x = o[1]
    theta = o[3]
    up_reward = cos(theta)
    distance_penalty_reward = -0.01 * (x^2)
    return up_reward + distance_penalty_reward
end
@propagate_inbounds function LyceumBase.geteval(s, a, o, env::CartpoleGym)
    x = o[1]
    theta = o[3]
    up_reward = cos(theta)
    distance_penalty_reward = -0.01 * (x^2)
    return up_reward + distance_penalty_reward
end
@propagate_inbounds function LyceumBase.reset!(env::CartpoleGym)
    randreset!(env)
end
@propagate_inbounds function LyceumBase.randreset!(env::CartpoleGym)
    env.state .= 0.
    env.state .= rand(Uniform(-0.05, 0.05), length(env.state))
    getobs!(env.obs, env)
    env.time[1] = 0.0
    env
end

@propagate_inbounds function LyceumBase.step!(env::CartpoleGym)

    state, act, ctrlrange = env.state, env.action, env.ctrlrange
    x, x_dot, theta, theta_dot = state

    # In Benchmark paper  action set to 1 if >0 0 otherwise
    # In Gym it becomes self.force_mag = 10.0 if 1 -10. otherwise
    force = 0.
    if act[1] > 0.; force = 10.; else; force = -10.; end;
    sintheta, costheta = sincos(theta)

    # self.polemass_length = 0.1 * 0.5 = 0.05 self.total_mass=1.0+0.1=1.1
    # self.gravity = 9.8
    temp     = (force + 0.05 * theta_dot^2 * sintheta) / 1.1
    thetaacc = (9.8 * sintheta - costheta * temp) / (0.5 * (4.0 / 3.0 - 0.1 * costheta^2 / 1.1))
    xacc     = temp - 0.05 * thetaacc * costheta / 1.1

    # self.tau = 0.02   EULER default
    x     = x + 0.02 * x_dot;          x_dot     = x_dot + 0.02 * xacc
    theta = theta + 0.02 * theta_dot;  theta_dot = theta_dot + 0.02 * thetaacc

    env.state[1] = x;  env.state[2] = x_dot;  env.state[3] = theta;  env.state[4] = theta_dot
    env.time[1] += timestep(env)
    env
end

@inline LyceumBase.timestep(env::CartpoleGym) = 1.0
@inline Base.time(env::CartpoleGym) = env.time[1]

#199.85387503411363


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

function write_J(fn,J)
  println("Writing to $(fn)")
  fid = open(fn,"w")
  @printf fid "%i %i\n" size(J,1) size(J,2)
  for i in 1:size(J,1)
    for j in 1:size(J,2)
      @printf fid "%f " J[i,j]
    end
    @printf fid "\n"
  end
  close(fid)
end

function write_times(fn,tf,tJ)
  println("Writing to $(fn)")
  fid = open(fn,"w")
  @printf fid "%f %f\r\n" tf tJ
  @printf fid "tf tJ\r\n"
  close(fid)
end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
using KeywordDispatch
using Base.Test

# write your own tests here

add = KeywordDispatchFunction()
# Empty: throws an error "not found"
@test_throws KeywordDispathFunctionNotFound add(x = 2, y = 3)

# Add specializer for x & y
@defkwspec add [:x :y] (args) begin
	args[:x] + args[:y]
end

# Throws unimplemented for x = 2
@test_throws KeywordDispatchFunctionNotImplemented add(x = 2)

@test add(x = 2, y = 3) == 5

@defkwspec add [:x] (args) begin
	args[:x]
end

@test add(x = 2) == 2
@test add(x = 2, z = 7) == 2
@test add(x = 2, y = 3) == 5

remove_specializer!(add, [:x, :y])

@test add(x = 2) == 2
@test_throws KeywordDispatchFunctionNotImplemented add(x = 2, y = 3)
@test add(x = 2, z = 3) == 2


text:
using HDF5,JLD,KUparser,KUnet,Base.Test
@date @load "acl11.trn.jld4"
@show net=newnet("bp10acl11eager13.jld")
ft = Flist.acl11eager
pt = ArcEager13
ncpu = 20
ndeps = length(deprel)
nbatch = 128
nbeam = 64
@date pxy1=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch, ncpu; xy=true)
@date pxy2=bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch; xy=true)
@test @show isequal(pxy1, pxy2)

# Space: Need ncpu x pbatch x nbeam < 7680 to fit in K20 memory

# * iui-k40-memory: Each K40 has 11.5GB for a total of 23GB.  The
# ZN11 ArcEager13 setup takes in MB GPU memory:
#   280*ncpu+0.1*ncpu*nbeam*nbatch
# For ncpu=20, nbeam=64, nbatch=128 we hit the upper limit.

# With that configuration 20cpu parse of acl11.trn takes 2336s.


text:
using Econometrics, Test, Random
function main()
    # dstats
    @testset "dstats" begin
        @test dstats()
    end
    # vech
    @testset "vech" begin
        @test vech()
    end
    # OLS
    x = [ones(10,1) (1:10)]
    Random.seed!(1)
    y = rand(10,1)
    b = ols(y,x,silent=true)[1]
    @testset "ols" begin
        @test b[1,1] ≈ 0.07655377367377458
    end
    # ML
    @testset "mle" begin
        b = mle()
        @test b[1,1] ≈ 0.49466773000112035
    end
    # GMM
    @testset "gmm" begin
        b = gmm()
        @test b[1,1] ≈ 0.4946717833359474
    end
    # samin
    xopt, fopt = samin(1)
    @testset "samin" begin
    @test xopt[1] ≈ 0.0 atol=1e-5
    @test fopt ≈ 2.0 atol=1e-5
    end
    #=
    # this is ill-advised, need to test with non-random data,
    # so it doesn't break when julia changes
    # npreg
    @testset "npreg" begin
        yhat = npreg(1);
        @test yhat[1,1] ≈ 0.02080984791302587
    end
    =#
end
main()

max_issues_repo_name: DawievLill/Econometrics

max_forks_repo_name: DawievLill/Econometrics


text:
module BlockArraysCooTests

using Test
using Gridap.Arrays
using BlockArrays
using LinearAlgebra

blocks = Matrix{Float64}[ [1 2; 3 4], [5 6; 7 8; 9 10] ]
blockids = [(1,1),(2,1)]
ax = (blockedrange([2,3]), blockedrange([2,4]))

a = BlockArrayCoo(blocks,blockids,ax)
@test a[Block(1),Block(1)] === blocks[1]
@test a[Block(1,1)] === blocks[1]
@test a[Block(1,2)] === a.zero_blocks[1]
@test a[BlockIndex((1,1),(2,1))] === blocks[1][2,1]
@test a[BlockIndex(1,2),BlockIndex(1,1)] === blocks[1][2,1]
@test a[2,1] === blocks[1][2,1]
@test a[3,2] === blocks[2][1,2]
@test axes(a) === ax
@test size(a) == (5,6)

@test is_zero_block(a,2,2) == true
@test is_zero_block(a,2,1) == false
@test is_nonzero_block(a,2,1) == true
@test is_zero_block(a,Block(1,2)) == true
@test is_zero_block(a,Block(1),Block(2)) == true

for (i,b) in enumerateblocks(a)
  @test a[i] === b
end

b21 = zeros(3,2)
getblock!(b21,a,Block(2,1))
@test b21 == a[Block(2,1)]

b12 = ones(2,4)
getblock!(b12,a,Block(1,2))
@test b12 == a[Block(1,2)]

blocks = Vector{Float64}[ [1,2,3] ]
blockids = [(2,)]
ax = (blockedrange([2,3,4]),)
a = BlockArrayCoo(blocks,blockids,ax)

@test a[Block(1)] === a.zero_blocks[1]
@test a[Block(2)] === blocks[1]
@test a[Block(3)] === a.zero_blocks[2]
@test a[BlockIndex(2,3)] === blocks[1][3]

blocks = Matrix{Float64}[ [1 2; 3 4], [5 6 7 8; 9 10 11 12; 13 14 15 16], [1 2 3 4; 5 6 7 8], [1 2 3; 4 5 6; 7 8 9] ]
blockids = [(1,1),(2,2),(1,2),(3,3)]
ax = (blockedrange([2,3,3]), blockedrange([2,4,3]))
a = BlockArrayCoo(blocks,blockids,ax)

blocks = Vector{Float64}[ 10*[1,2], 20*[1,2,3] ]
blockids = [(1,),(3,)]
ax = (blockedrange([2,4,3]),)
b = BlockArrayCoo(blocks,blockids,ax)

c = a*b
@test axes(c,1) === axes(a,1)
@test blocksize(c) == (3,)
@test Array(a)*Array(b) == c

mul!(c,a,b)
@test axes(c,1) === axes(a,1)
@test blocksize(c) == (3,)
@test Array(a)*Array(b) == c

d = copy(c)
mul!(d,a,b,2,3)
@test axes(d,1) === axes(a,1)
@test blocksize(d) == (3,)
@test 2*Array(a)*Array(b) + 3*Array(c) == d

b = transpose(a)
@test isa(b,BlockArrayCoo)
c = a*b
@test axes(c,1) === axes(a,1)
@test axes(c,2) === axes(a,1)
@test blocksize(c) == (3,3)
@test Array(a)*Array(b) == c

mul!(c,a,b)
@test axes(c,1) === axes(a,1)
@test axes(c,2) === axes(a,1)
@test blocksize(c) == (3,3)
@test Array(a)*Array(b) == c

cc = CachedArray(c)

axs = (blockedrange([2,3,3]), blockedrange([2,3,3]))
setaxes!(cc,axs)
@test map(blocklasts,axes(cc.array)) == map(blocklasts,axs)
@test cc.array === c

axs = (blockedrange([4,5,3]), blockedrange([2,4,3]))
setaxes!(cc,axs)
@test map(blocklasts,axes(cc.array)) == map(blocklasts,axs)
fill!(cc.array,0)

axs = (blockedrange([5,4,3]), blockedrange([4,2,3]))
setaxes!(cc,axs)
@test map(blocklasts,axes(cc.array)) == map(blocklasts,axs)

axs1 = (blockedrange([5,4,3]), blockedrange([4,2,3]))
axs2 = (blockedrange([4,5,3]), blockedrange([2,4,3]))
axs3 = (blockedrange([4,5,3]), blockedrange([2,4,3]))
@test Arrays._same_axes(axs1,axs2) == false
@test Arrays._same_axes(axs2,axs2)
@test Arrays._same_axes(axs2,axs3)

blocks = [ 10*[1,2], 20*[1,2,3] ]
blockids = [(1,),(3,)]
axs = (blockedrange([2,4,3]),)
b = BlockArrayCoo(blocks,blockids,axs)

cb = CachedArray(b)

setaxes!(cb,axs)
@test map(blocklasts,axes(cb.array)) == map(blocklasts,axs)
@test cb.array === b

axs = (blockedrange([3,2,3]),)
setaxes!(cb,axs)
@test map(blocklasts,axes(cb.array)) == map(blocklasts,axs)
@test size(cb) == (8,)

c = copy(a)
@test isa(c,BlockArrayCoo)
@test c == a
@test axes(c) == axes(a)

fill!(c,0)
copyto!(c,a)
@test c == a
@test axes(c) == axes(a)

c = 2*a
@test isa(c,BlockArrayCoo)
@test 2*Array(a) == c

c = a*2
@test isa(c,BlockArrayCoo)
@test 2*Array(a) == c

c = a + a
@test isa(c,BlockArrayCoo)
@test 2*Array(a) == c

blocks = [ [1 2; 3 4], [5 6 7 8; 9 10 11 12; 13 14 15 16], [1 2 3 4; 5 6 7 8], [1 2 3; 4 5 6; 7 8 9] ]
blockids = [(1,1),(2,2),(1,2),(3,3)]
ax = (blockedrange([2,3,3]), blockedrange([2,4,3]))
a = BlockArrayCoo(blocks,blockids,ax)

blocks = [ [1 2; 3 4], [1 2; 4 5; 8 9], [5 6 7 8; 9 10 11 12; 13 14 15 16] ]
blockids = [(1,1),(2,1),(3,2)]
b = BlockArrayCoo(blocks,blockids,ax)

c = a + b
@test isa(c,BlockArrayCoo)
@test Array(a)+Array(b) == c

c = a - 2*b
@test isa(c,BlockArrayCoo)
@test Array(a)-2*Array(b) == c

@test a[1,2,1] == a[1,2]

blocks = [ 10*[1,2], 20*[1,2,3] ]
blockids = [(1,),(3,)]
axs = (blockedrange([2,4,3]),)
b = BlockArrayCoo(blocks,blockids,axs)

@test b[2] == b[2,1]

axs = (blockedrange([2,4,3]),)

# Blocks of Blocks

bL = BlockArrayCoo([[1,2]],[(1,)],axs)
bR = BlockArrayCoo([[2,3],[4,5,6]],[(1,),(3,)],axs)
axs = (blockedrange([axs[1],axs[1]]),)

bS = BlockArrayCoo([bL,bR],[(1,),(2,)],axs)

zS = zeros_like(bS)
@test isa(zS,BlockArrayCoo)
@test isa(zS[Block(1)],BlockArrayCoo)
@test isa(zS[Block(2)],BlockArrayCoo)
@test length(zS.blocks) == 0
@test length(zS[Block(1)].blocks) == 0
@test length(zS[Block(2)].blocks) == 0

cS = 2*bS
@test isa(cS,BlockArrayCoo)
@test isa(cS[Block(1)],BlockArrayCoo)

bZ = zeros_like(bR)

dS = BlockArrayCoo([bZ,bR],axs,[false,true])

r = cS - dS
@test isa(r,BlockArrayCoo)
@test isa(r[Block(1)],BlockArrayCoo)

r = cS + dS
@test isa(r,BlockArrayCoo)
@test isa(r[Block(1)],BlockArrayCoo)

ax = (blockedrange([2,3]), blockedrange([2,4]))
aLL = BlockArrayCoo([[1 2; 3 4],[5 6; 7 8; 9 10] ],[(1,1),(2,1)],ax)
aLR = BlockArrayCoo([[1 2 5 6; 3 4 1 2; 1 2 3 4],[5 6; 7 8; 9 10] ],[(2,2),(2,1)],ax)
aRL = zeros_like(aLL)
aRR = BlockArrayCoo([[1 2; 3 4],[5 6; 7 8; 9 10] ],[(1,1),(2,1)],ax)

allblocks = Matrix{typeof(aLL)}(undef,2,2)
allblocks[1,1] = aLL
allblocks[1,2] = aLR
allblocks[2,1] = aRL
allblocks[2,2] = aRR

mask = [true true; false true]

axs = ( blockedrange([ax[1],ax[1]]) , blockedrange([ax[2],ax[2]]))

aS = BlockArrayCoo(allblocks,axs,mask)
@test is_zero_block(aS,2,1)
@test is_nonzero_block(aS,2,2)

aS2 = copy(aS)
@test isa(aS2,BlockArrayCoo)
@test isa(aS2[Block(2),Block(2)],BlockArrayCoo)

aSt = transpose(aS)
@test Array(aSt) == transpose(Array(aS))
@test isa(aSt,BlockArrayCoo)
@test isa(aSt[Block(2),Block(2)],BlockArrayCoo)

rS = aS*aSt
@test rS == Array(aS)*Array(aSt)
@test isa(rS,BlockArrayCoo)
@test isa(rS[Block(2),Block(2)],BlockArrayCoo)

mul!(rS,aS,aSt)
@test rS == Array(aS)*Array(aSt)
@test isa(rS,BlockArrayCoo)
@test isa(rS[Block(2),Block(2)],BlockArrayCoo)

axs =(blockedrange([2,4]),)
bL = BlockArrayCoo([[1,2]],[(1,)],axs)
bR = BlockArrayCoo([[2,3],[4,5,6,8]],[(1,),(2,)],axs)
bS = BlockArrayCoo([bL,bR],[(1,),(2,)],(blockedrange([6,6]),))

rS = aS*bS
@test rS == Array(aS)*Array(bS)
@test isa(rS,BlockArrayCoo)
@test isa(rS[Block(2)],BlockArrayCoo)

mul!(rS,aS,bS)
@test rS == Array(aS)*Array(bS)
@test isa(rS,BlockArrayCoo)
@test isa(rS[Block(2)],BlockArrayCoo)

cS = copy(rS)
mul!(rS,aS,bS,3,2)
@test rS == 3*Array(aS)*Array(bS) + 2*cS
@test isa(rS,BlockArrayCoo)
@test isa(rS[Block(2)],BlockArrayCoo)

axs1 = (blockedrange([5,4,3]), blockedrange([4,2,3]))
axs2 = (blockedrange([4,5,3]), blockedrange([2,4,3]))
axsA = (blockedrange([axs1[1],axs1[1]]),blockedrange([axs2[2],axs2[2]]))
axsB = (blockedrange([axs2[1],axs2[1]]),blockedrange([axs1[2],axs1[2]]))
axsC = (blockedrange([axs2[1],axs2[1]]),blockedrange([axs1[2],axs1[2]]))
@test Arrays._same_axes(axsA,axsB) == false
@test Arrays._same_axes(axsA,axsA)
@test Arrays._same_axes(axsB,axsC)

#using BenchmarkTools
#@btime Arrays._same_axes($axsA,$axsA)

end # module


text:
"""
    LocalizationGroup

A set of [`Localization`](@ref)s determined to be associated together and treated as a single entity.
"""
mutable struct LocalizationGroup <: DataEntity
    index::Int            # unique identifier within this dataset
    channel::String
    x::Float64
    y::Float64
    z::Float64
    localizations::Vector{Localization}

    LocalizationGroup(localization::Localization) =
        new(localization.index, localization.channel, localization.x, localization.y, localization.z, [localization])
end

import Base.push!
function push!(group::LocalizationGroup, newlocalization::Localization...)
    push!(group.localizations, newlocalization...)
    calcposition!(group)
end

import Base.append!
function append!(group::LocalizationGroup, newlocalizations::Vector{Localization})
    append!(group.localizations, newlocalizations)
    calcposition!(group)
end

append!(group::LocalizationGroup, othergroup::LocalizationGroup) =
    append!(group, othergroup.localizations)

function calcposition!(group::LocalizationGroup)
    count = length(group.localizations)
    group.x = sum(l -> l.x, group.localizations) / count
    group.y = sum(l -> l.y, group.localizations) / count
    group.z = sum(l -> l.z, group.localizations) / count
    return group
end


text:
using FFTW, LinearAlgebra

import Splittings:UniformMesh, compute_rho, compute_e
import Splittings:Landau, distribution

mesh1 = UniformMesh( 0., 4π,  64, endpoint=false)
mesh2 = UniformMesh(-6., 6., 128, endpoint=false)

ϵ, kx = 0.001, 0.5

f = zeros(Complex{Float64},(mesh1.length,mesh2.length))

f .= distribution( mesh1, mesh2, Landau(ϵ, kx))

@testset "Poisson 1D using FFT" begin

ρ = compute_rho(mesh2, f)
@test real(ρ) ≈ ϵ*cos.(kx*mesh1.points)
e = compute_e(mesh1, ρ)
@test real(e) ≈ ϵ*sin.(kx*mesh1.points)/kx

import Splittings:compute_rho!, compute_e!

ex  = zeros(Complex{Float64},mesh1.length)
rho = similar(ex)

compute_rho!(rho, mesh2, f)
@test real(rho) ≈ ϵ*cos.(kx*mesh1.points)
compute_e!(ex, mesh1, rho)
@test real(ex) ≈ ϵ*sin.(kx*mesh1.points)/kx

end


text:
module Teuchos
using CxxWrap, MPI
import ..libjltrilinos
import ..CxxUnion
import Base.get

export ParameterList

abstract type PLAssociative <: CxxWrap.CppAssociative{String, Any} end

struct ArrayView{T}
  array::Ptr{T}
  size::Int
end

ArrayView(a::AbstractArray{T,1}) where {T} = ArrayView(pointer(a),length(a))
ArrayView(a::AbstractArray{T,1}, length::Integer) where {T} = ArrayView(pointer(a),length)
Base.convert(::Type{ArrayView{T}}, a::AbstractArray{T,1}) where {T} = ArrayView(a, length(a))

@readmodule(libjltrilinos, :register_teuchos)
@wraptypes

CxxWrap.argument_overloads(t::Type{ArrayView{T}}) where {T} = [AbstractArray{T,1}]

@wrapfunctions

const ParUnion = Union{CxxWrap.SmartPointer{ParameterList}, ParameterList}

get(::Type{Nothing}, ::ParUnion, ::String) = "nullptr"

# High-level interface for ParameterList
Base.length(pl::ParUnion) = numParams(pl)
function Base.getindex(pl::ParUnion, key)
  if !isParameter(pl, key)
    throw(KeyError(key))
  end
  if isSublist(pl, key)
    return sublist(pl, key)
  end
  return get(get_type(pl, key), pl, key)
end
Base.setindex!(pl::ParUnion, v, key) = set(pl, key, v)
Base.setindex!(pl::ParUnion, v::Int64, key) = set(pl, key, Int32(v))
Base.setindex!(pl::ParUnion, v::CppEnum, key) = set(pl, key, convert(Int32,v))
Base.keys(pl::ParUnion) = keys(pl)

function Base.iterate(pl::ParUnion, state=(1, keys(pl)))
  idx, keys = state

  if idx >= length(keys)
    return nothing
  end

  return ((keys[idx], pl[keys[idx]]), (idx+1,keys))
end

Base.haskey(pl::ParUnion, key) = isParameter(pl, key)
Base.isempty(pl::ParUnion) = (numParams(pl) == 0)

get(pl::ParUnion, key, default_value) = isParameter(pl, key) ? pl[key] : default_value

"""
A pair of ParameterLists, one for input and one for storing the actually used parameters and any added default values
"""
struct ParameterListPair
  input
  output::ParUnion
end

function sublist(pl::ParameterListPair, name)
  param_val = get(pl.input, name, Teuchos.ParameterList(name))
  return ParameterListPair(param_val, Teuchos.sublist(pl.output, name))
end

function get(pl::ParameterListPair, key, default_value)
  result = get(pl.input, key, default_value)
  pl.output[key] = result
  return result
end

# Interface to Array
Teuchos.Array(nelems::Integer, elem::T) where {T} = Teuchos.Array{T}(nelems, elem)
Teuchos.Array(nelems::Integer, elem::T) where {T <: AbstractString} = Teuchos.Array{AbstractString}(nelems, elem)
Base.IndexStyle(::Teuchos.Array) = IndexLinear()
Base.size(arr::Teuchos.Array) = (size(arr),)
Base.getindex(arr::Teuchos.Array, i::Integer) = at(arr,i-1)

# Convert Teuchos Comm to MPI.Comm
MPI.Comm(comm::CxxUnion{Teuchos.MpiComm}) = MPI.Comm(getRawMpiComm(comm))

end


text:
#volume_solver(model::HelmholtzModel,args...;kwargs...) = volume_solver(model_type(model),volume_solver_type(model),model,args...;kwargs...)

function fugacity_coeff_impl(mt::SingleVT,model::HelmholtzModel,v,t)
    rho = one(v)/v
    ar(_rho) = αR_impl(mt,model,_rho,t)
    (ar,dadrho) = f_dfdx(ar,rho)
    p = rho*(one(rho)+rho*dadrho)*RGAS*t
    Z = p*v/(RGAS*t)
    return ar + (Z-1) - log(Z)
end

function flash_impl(mt::SingleSatP,model::HelmholtzModel, _p)
    _pt = QuickStates.pt()
    p = float(normalize_units(_p))
    TTT = typeof(p)

    if (Pc = pressure(model,CriticalPoint())) ≈ p
        vc = mol_volume(model,CriticalPoint())
        tc=temperature(model,CriticalPoint())
        return state(mol_v=vc,t=tc)
    elseif p > Pc
        throw(error("the phase is supercritical at pressure = $p Pa"))
    else
    end
    Tc = temperature(model,CriticalPoint())
    model_pred = single_sat_aprox(model)
    t_pred = temperature_impl(mt,model_pred,p) #prediction temperature
    vv = v_zeros(_pt,model,p,t_pred)

    v1 = first(vv)
    #(v1 == v2) && return v1
    v2 = last(vv)
    if v1 ≈ v2
        throw(error("cannot converge at pressure = $p Pa"))
    else
        px = p
        _A(z, t) = mol_helmholtz_impl(QuickStates.vt(),model, z, t)
        v1old = zero(TTT)
        v2old = zero(TTT)
        Told = TTT(-100.0)
        Tx = t_pred
        for i = 1:20
            if i > 1
                v1old = v1
                v2old = v2
            end
            A = z -> _A(z, Tx)
            _v1 = Threads.@spawn v_zero($_pt,model,$p,$Tx,0.95*$v1;phase=:liquid)
            _v2 = Threads.@spawn v_zero($_pt,model,$p,$Tx,1.05*$v2;phase=:gas)
            v1 = fetch(_v1)
            v2 = fetch(_v2)

            if abs(v1 - v1old) / v1 < 1e-15 && i > 1
                #println("v1 condition")
                break
            elseif abs(v2 - v2old) / v2 < 1e-15 && i > 1
                #println("v2 condition")
                break
            end

            _px(T) = (_A(v1, T) - _A(v2, T)) / (v2 - v1) - p
            #_px(T) = exp(_A(v1,T)-_A(v2,T)) - one(Tx)
            Told = Tx
            Tx = Roots.find_zero(_px, Tx)
            if abs(Told - Tx) < 1e-10 * Tx
                #println("P condition")
                break
            end
        end
    end

    phase1 = state(mol_v = v1, t=Tx ,phase=:liquid)
    phase2 = state(mol_v = v2, t=Tx, phase=:gas)
    return (phase1, phase2)
end

function flash_impl(mt::SingleSatT,model::HelmholtzModel, _t)
    _pt = QuickStates.pt()
    t = float(normalize_units(_t))
    tc = ThermoState.temperature(model,CriticalPoint())
    if tc ≈ t
        return state(mol_v=mol_volume(model,CriticalPoint()),t=tc)
    elseif t > tc
        throw(error("the phase is supercritical at temperature = $T K"))
    end
    model_pred = single_sat_aprox(model)
    p_pred = pressure_impl(mt,model_pred,t)
    p = p_pred
    _p(z) = pressure_impl(QuickStates.vt(),model, z, t)

    vv = v_zeros(_pt,model,p,t)

    v1 = first(vv)
    #(v1 == v2) && return v1
    v2 = last(vv)
    if v1 ≈ v2
        throw(error("cannot converge at temperature = $T K"))
    else

        p1 = _p(v1)
        p2 = _p(v2)
        px = p_pred
        ϕ(z) = fugacity_coeff_impl(QuickStates.vt(),model, z, t)
        v1old = 0
        v2old = 0
        for i = 1:20
            if i > 1
                v1old = v1
                v2old = v2
            end

            _v1 = Threads.@spawn v_zero($_pt,model,$px,$t,0.9*v1;phase=:liquid)
            _v2 = Threads.@spawn v_zero($_pt,model,$px,$t,1.05*v2;phase=:gas)

            v1 = fetch(_v1)
            v2 = fetch(_v2)
            if abs(v1 - v1old) / v1 < 1e-15 && i > 1
                #println("v1 condition")
                break
            elseif abs(v2 - v2old) / v2 < 1e-15 && i > 1
                #println("v2 condition")
                break
            end
            pold = px
            px = px*exp(ϕ(v1)-ϕ(v2))
            #px = (A(v1) - A(v2)) / (v2 - v1)
            #@show px
            if abs(pold - px) < 1e-10 * px
                #println("P condition")
                break
            end

        end
        #exp(logphi)*P= f
        phase1 = state(mol_v = v1, t= t,phase=:liquid)
        phase2 = state(mol_v = v2, t= t,phase=:gas)
        return (phase1, phase2)
    end
end

function equilibria(model::HelmholtzModel,st::ThermodynamicState)
    return equilibria(state_type(st),model,st)
end

function equilibria(mt::SingleSatP,model::HelmholtzModel,st::ThermodynamicState)
    p = pressure(FromState(),st)
    return flash_impl(mt,model,p)
end

function equilibria(mt::SingleΦP,model::HelmholtzModel,st::ThermodynamicState)
    p = pressure(FromState(),st)
    ϕ = quality(FromState(),st)
    m = mass(FromState()st,u"kg",molecular_weight(model))
    mv = ϕ*m
    ml = (1-ϕ)*m
    st_l, st_g = flash_impl(mt,model,p)
    Teq = temperature(FromState(),st_l)
    vl = mol_volume(FromState(),st_l)
    vg = mol_volume(FromState(),st_g)
    tagl = :liquid
    tagv = :gas
    st1 = state(mass=ml,t=Teq,mol_v=vl,phase=tagl)
    st2 = state(mass=mv,t=Teq,mol_v=vv,phase=tagv)
    return (st1,st2)
end

function equilibria(mt::SingleSatT,model::HelmholtzModel,st::ThermodynamicState)
    t = temperature(FromState(),st)
    return flash_impl(mt,model,t)
end

function equilibria(mt::SingleΦT,model::HelmholtzModel,st::ThermodynamicState)
    t = temperature(FromState(),st)
    ϕ = quality(FromState(),st)
    m = mass(FromState()st,u"kg",molecular_weight(model))
    mv = ϕ*m
    ml = (1-ϕ)*m
    st_l, st_g = flash_impl(mt,model,t)
    Peq = pressure(FromState(),st_l)
    vl = mol_volume(FromState(),st_l)
    vg = mol_volume(FromState(),st_g)
    tagl = :liquid
    tagv = :gas
    st1 = state(mass=ml,p=Peq,mol_v=vl,phase=tagl)
    st2 = state(mass=mv,p=Peq,mol_v=vv,phase=tagv)
    return (st1,st2)
end

function pressure(mt::SingleSatT,model::HelmholtzModel,st::ThermodynamicState,unit)
    t = temperature(FromState(),st)
    st_l, st_g = flash_impl(mt,model,t)
    p1 = pressure(model,st_l)
    p2 = pressure(model,st_g)
    Peq = 0.5*(p1+p2)
    return convert_unit(u"Pa",unit,Peq)
end

function temperature(mt::SingleSatP,model::HelmholtzModel,st::ThermodynamicState,unit)
    p = pressure(FromState(),st)
    st_l, st_g = flash_impl(mt,model,p)
    Teq = temperature(FromState(),st_l)
    return convert_unit(u"K",unit,Teq)
end

function pressure(mt::SingleΦT,model::HelmholtzModel,st::ThermodynamicState,unit)
    return pressure(QuickStates.sat_t(),model,st,unit)
end

function temperature(mt::SingleΦP,model::HelmholtzModel,st::ThermodynamicState,unit)
    return temperature(QuickStates.sat_p(),model,st,unit)
end


text: using CuYao: cu


text:
"""
    all_libs()
"""
function all_libs()
    libs = Symbol[
        :CUBLAS,
        :GemmKernels,
        :Tullio,
    ]
    return libs
end

function getbenchfunc(lib::Symbol)
    if lib === :CUBLAS
        return benchmark_gemm_cublas!
    elseif lib === :GemmKernels
        return benchmark_gemm_gemmkernels!
    elseif lib === :Tullio
        return benchmark_gemm_tullio!
    end
    throw(ArgumentError("Unknown library \"$(lib)\""))
end

function gemm_cublas!(C, A, B)
    CUDA.CUBLAS.cublasSetMathMode(CUDA.CUBLAS.handle(), CUDA.CUBLAS.CUBLAS_TENSOR_OP_MATH)
    CUDA.CUBLAS.gemmEx!('N', 'N', 1, A, B, 0, C)
    return C
end

function benchmark_gemm_cublas!(C, A, B)
    CUDA.CUBLAS.cublasSetMathMode(CUDA.CUBLAS.handle(), CUDA.CUBLAS.CUBLAS_TENSOR_OP_MATH)
    return BenchmarkTools.@benchmark CUDA.@sync begin
        CUDA.CUBLAS.gemmEx!('N', 'N', 1, $A, $B, 0, $C)
    end
end

function gemm_gemmkernels!(C, A, B)
    GemmKernels.BLAS.gemmEx!('N', 'N', 1, A, B, 0, C)
    return C
end

function benchmark_gemm_gemmkernels!(C, A, B)
    return BenchmarkTools.@benchmark CUDA.@sync begin
        GemmKernels.BLAS.gemmEx!('N', 'N', 1, $A, $B, 0, $C)
    end
end

function gemm_tullio!(C, A, B)
    Tullio.@tullio C[i, j] = A[i, k] * B[k, j] cuda=256 grad=false
    return C
end

function benchmark_gemm_tullio!(C, A, B)
    return BenchmarkTools.@benchmark CUDA.@sync begin
        gemm_tullio!($C, $A, $B)
    end
end


text:
using RoadRunner
using Test

ant_str = """
    k1 := sin (time)
    k2 := k1 + 3.14

    S1 -> S2; v
    S1 = 1; S2 = 0; v = 0
"""

rr = RoadRunner.loada(ant_str)

@testset "compartment" begin
    @test RoadRunner.getNumberOfCompartments(rr) == 1
end

@testset "reaction" begin
    @test RoadRunner.getNumberOfReactions(rr) == 1
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase, Test

@testset "get integration points for Element{Seg2}" begin
    element = Element(Seg2, (1, 2))
    ips = get_integration_points(element)
    @test length(ips) == 2
end


text:
function validateGraphVizInstalled()
    # Check if GraphViz is installed
    try
        (read(`dot -'?'`, String)[1:10] == "Usage: dot") || error()
    catch
        error("GraphViz is not installed correctly. Make sure GraphViz is installed. If you are on Windows, manually add the path to GraphViz to your path variable. You should be able to run 'dot' from the command line.")
    end
end

function dot2png(dot_graph::AbstractString)
    # Generate PNG image from DOT graph
    validateGraphVizInstalled()
    proc = open(`dot -Tpng`, "r+")
    write(proc.in, dot_graph)
    close(proc.in)
    return read(proc.out, String)
end
function dot2jpg(dot_graph::AbstractString)
    # Generate PNG image from DOT graph
    validateGraphVizInstalled()
    proc = open(`dot -Tjpg`, "r+")
    write(proc.in, dot_graph)
    close(proc.in)
    return read(proc.out, String)
end

export dot2png,dot2jpg

colors = Dict(
    :failure =>"#ff9a20",
    :running => "#00bbd3",
    :success => "#49b156",
)

function toDotContent(tree, parent_id)
    name = BehaviorTree.format(tree)
    node_id = string(parent_id, replace(name, "!"=>""))
    if length(children(tree)) == 0
        shape = if startswith(name, "is") "ellipse" else "box" end
        return string([
            """$node_id:n\n""",
            """$node_id [shape=$shape, label="$name"]\n""",
            """$parent_id -> $node_id""",
        ]...)
    end
    shape = "box"
    label = if(typeof(tree) == Selector) "?" else "->" end
    out = string(
        """$node_id:n\n""",
        """$node_id [shape=$shape, label="$label"]\n""",
        ["""$(toDotContent(c, node_id))\n""" for c in children(tree)]...)
    if parent_id != ""
        out = string(
            out,
            """$parent_id -> $node_id""",
        )
    end
    out
end

function toStatusDotContent(tree, parent_id)
    name = BehaviorTree.format(tree.tree.x)
    node_id = string(parent_id, replace(name, "!"=>""))
    if length(children(tree)) == 0
        shape = if startswith(name, "is") "ellipse" else "box" end
        color = colors[tree.shadow.x]
        return string([
            """$node_id:n\n""",
            """$node_id [shape=$shape, label="$name", style=filled,color="$color"]\n""",
            """$parent_id -> $node_id""",
        ]...)
    end
    shape = "box"
    label = if(typeof(tree.tree.x) == Selector) "?" else "->" end
    #TODO: more efficient implementation?
    lastchild = tree.shadow.x
    while length(children(lastchild)) > 0
        lastchild = last(children(lastchild))
    end
    color = colors[lastchild]
    out = string(
        """$node_id:n\n""",
        """$node_id [shape=$shape, label="$label", style=filled, color="$color"]\n""",
        ["""$(toStatusDotContent(c, node_id))\n""" for c in children(tree)]...)
    if parent_id != ""
        out = string(
            out,
            """$parent_id -> $node_id""",
        )
    end
    out
end
function toDot(tree::BT, results)
    st = ShadowTree(tree, results)
    content = toStatusDotContent(st, "")
    return """digraph tree {
    $(content)
    }"""
end
function toDot(tree::BT)
    content = toDotContent(tree, "")
    return """digraph tree {
    $(content)
    }"""
end

export toDot


text:
__precompile__(true)
module Optionals

import Base: convert, promote_rule, show, repr, isless, ==
import Missings: Missing, MissingException, ismissing, missing, coalesce
using Nullables

export Missing, ismissing, missing, coalesce, Optional

@static if VERSION < v"0.7.0-alpha"
    struct Optional{T}
        nullable::Nullable{T}
        Optional{T}() where T = new(Nullable{T}())
        Optional{T}(x::Any) where T = new(Nullable{T}(x))
        Optional{T}(::Missing) where T = new(Nullable{T}())
    end

    ismissing(x::Optional) = isnull(x.nullable)

    coalesce(x::Optional) =
        ismissing(x) ? missing : unsafe_get(x.nullable)

    _unwrap(x::Optional{T}) where T =
        ismissing(x) ? throw(MissingException()) : unsafe_get(x.nullable)

    repr(x::Optional) = ismissing(x) ? repr(missing) : repr(_unwrap(x))
    repr(mime::Union{MIME, AbstractString}, x::Optional) =
        ismissing(x) ? repr(mime, missing) : repr(mime, _unwrap(x))
else

    struct Optional{T}
        value::Union{T, Missing}
        Optional{T}() where T = new(missing)
        Optional{T}(x) where T = new(x)
        Optional{T}(::Missing) where T = new(missing)
    end

    ismissing(x::Optional) = x.value === missing

    coalesce(x::Optional) = x.value

    _unwrap(x::Optional{T}) where T = x.value::T

    repr(x::Optional, context=nothing) =
        ismissing(x) ?
          repr(missing, context=context) :
          repr(_unwrap(x), context=context)

    repr(mime::Union{MIME, AbstractString}, x::Optional, context=nothing) =
        ismissing(x) ?
          repr(mime, missing, context=context) :
          repr(mime, _unwrap(x), context=context)
end

Optional(x) = Optional{typeof(x)}(x)

convert(::Type{Optional{T}}, x::Optional{T}) where T = x
convert(::Type{Optional{T}}, x::Optional) where T =
    ismissing(x) ? Optional{T}() : Optional(convert(T, unwrap(x)))
convert(::Type{Optional{T}}, ::Missing) where T = Optional{T}()
convert(::Type{Optional{T}}, x::Any) where T = Optional{T}(convert(T, x))
convert(::Type{Optional}, x::Optional) = x
convert(::Type{Optional}, x::Any) = Optional(x)

promote_rule(::Type{Optional{T}}, ::Type{Optional{U}}) where {T, U} =
    Optional{promote_type(T, U)}
promote_rule(::Type{Optional{T}}, ::Type{Missing}) where T = Optional{T}
promote_rule(::Type{Optional{T}}, U::Type) where T =
    Optional{promote_type(T, U)}

coalesce(x::Optional, ys...) = coalesce(coalesce(x), ys...)
coelesce(x::Optional, y) = ismissing(x) ? y : _unwrap(x)

show(io::IO, x::Optional) =
    ismissing(x) ? show(io, missing) : show(io, _unwrap(x))
show(io::IO, mime, x::Optional) =
    ismissing(x) ? show(io, mime, missing) : show(io, mime, _unwrap(x))

function isless(x::Optional, y::Optional)
    ismissing(x) && return missing
    ismissing(y) && return missing
    return isless(_unwrap(x), _unwrap(y))
end

function isless(x, y::Optional)
    ismissing(y) && return missing
    return isless(x, _unwrap(y))
end

function isless(x::Optional, y)
    ismissing(x) && return missing
    return isless(_unwrap(x), y)
end

function ==(x::Optional, y::Optional)
    ismissing(x) && return missing
    ismissing(y) && return missing
    return _unwrap(x) == _unwrap(y)
end

function ==(x, y::Optional)
    ismissing(y) && return missing
    return x == _unwrap(y)
end

function ==(x::Optional, y)
    ismissing(x) && return missing
    return _unwrap(x) == y
end

end # module


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2020 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

# system packages - these are required to be installed to solve the modeling problem
# check: - are these are installed?
# y => all is good with the world.
# n => Ooops, we need to install them
using LinearAlgebra # pre-installed w/Julia
using Statistics    # pre-installed w/Julia
using Pkg           # pre-installed w/Julia
installed_package_set = keys(Pkg.installed())

# Do we have GLPK?
if (in("GLPK",installed_package_set) == false)
    Pkg.add("GLPK")
end

# Do we have JSON?
if (in("TOML",installed_package_set) == false)
    Pkg.add(PackageSpec(url="https://github.com/JuliaLang/TOML.jl.git"))
end

# Do we have DelimitedFiles?
if (in("DelimitedFiles",installed_package_set) == false)
    Pkg.add("DelimitedFiles")
end

# Do we have CSV?
if (in("CSV",installed_package_set) == false)
    Pkg.add("CSV")
end

# Do we have DataFrames?
if (in("DataFrames",installed_package_set) == false)
    Pkg.add("DataFrames")
end

# Do we have Logging?
if (in("Logging",installed_package_set) == false)
    Pkg.add("Logging")
end

# load the required system packages -
using GLPK
using DelimitedFiles
using TOML
using CSV
using DataFrames
using Logging

# load my model files -
include("./src/Checks.jl")
include("./src/Types.jl")
include("./src/Data.jl")
include("./src/Kinetics.jl")
include("./src/Control.jl")
include("./src/Solver.jl")
include("./src/Utility.jl")
include("./src/Flux.jl")
include("./src/Constraints.jl")

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG