text:
@test mysignif(SpecificG(3.0, 0.391_798_509E3), TestDigits) ≈ mysignif(SpecificG_Ph(3.0, 500.0), TestDigits)
@test mysignif(SpecificG(3.0u"MPa", 0.391_798_509E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(3.0u"MPa", 500.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificG(80.0, 0.378_108_626E3), TestDigits) ≈ mysignif(SpecificG_Ph(80.0, 500.0), TestDigits)
@test mysignif(SpecificG(80.0u"MPa", 0.378_108_626E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(80.0u"MPa", 500.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificG(80.0, 0.611_041_229E3), TestDigits) ≈ mysignif(SpecificG_Ph(80.0, 1500.0), TestDigits)
@test mysignif(SpecificG(80.0u"MPa", 0.611_041_229E3u"K"), TestDigits) ≈ mysignif(SpecificG_Ph(80.0u"MPa", 1500.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(3.0, 0.391_798_509E3), TestDigits) ≈ mysignif(SpecificF_Ph(3.0, 500.0), TestDigits)
@test mysignif(SpecificF(3.0u"MPa", 0.391_798_509E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(3.0u"MPa", 500.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(80.0, 0.378_108_626E3), TestDigits) ≈ mysignif(SpecificF_Ph(80.0, 500.0), TestDigits)
@test mysignif(SpecificF(80.0u"MPa", 0.378_108_626E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(80.0u"MPa", 500.0u"kJ/kg"), TestDigits)

@test mysignif(SpecificF(80.0, 0.611_041_229E3), TestDigits) ≈ mysignif(SpecificF_Ph(80.0, 1500.0), TestDigits)
@test mysignif(SpecificF(80.0u"MPa", 0.611_041_229E3u"K"), TestDigits) ≈ mysignif(SpecificF_Ph(80.0u"MPa", 1500.0u"kJ/kg"), TestDigits)

@test SpecificV(3.0, 0.391_798_509E3) ≈ SpecificV_Ph(3.0, 500.0)
@test SpecificV(3.0u"MPa", 0.391_798_509E3u"K") ≈ SpecificV_Ph(3.0u"MPa", 500.0u"kJ/kg")

@test SpecificV(80.0, 0.378_108_626E3) ≈ SpecificV_Ph(80.0, 500.0)
@test SpecificV(80.0u"MPa", 0.378_108_626E3u"K") ≈ SpecificV_Ph(80.0u"MPa", 500.0u"kJ/kg")

@test SpecificV(80.0, 0.611_041_229E3) ≈ SpecificV_Ph(80.0, 1500.0)
@test SpecificV(80.0u"MPa", 0.611_041_229E3u"K") ≈ SpecificV_Ph(80.0u"MPa", 1500.0u"kJ/kg")

@test SpecificU(3.0, 0.391_798_509E3) ≈ SpecificU_Ph(3.0, 500.0)
@test SpecificU(3.0u"MPa", 0.391_798_509E3u"K") ≈ SpecificU_Ph(3.0u"MPa", 500.0u"kJ/kg")

@test SpecificU(80.0, 0.378_108_626E3) ≈ SpecificU_Ph(80.0, 500.0)
@test SpecificU(80.0u"MPa", 0.378_108_626E3u"K") ≈ SpecificU_Ph(80.0u"MPa", 500.0u"kJ/kg")

@test SpecificU(80.0, 0.611_041_229E3) ≈ SpecificU_Ph(80.0, 1500.0)
@test SpecificU(80.0u"MPa", 0.611_041_229E3u"K") ≈ SpecificU_Ph(80.0u"MPa", 1500.0u"kJ/kg")

@test SpecificS(3.0, 0.391_798_509E3) ≈ SpecificS_Ph(3.0, 500.0)
@test SpecificS(3.0u"MPa", 0.391_798_509E3u"K") ≈ SpecificS_Ph(3.0u"MPa", 500.0u"kJ/kg")

@test SpecificS(80.0, 0.378_108_626E3) ≈ SpecificS_Ph(80.0, 500.0)
@test SpecificS(80.0u"MPa", 0.378_108_626E3u"K") ≈ SpecificS_Ph(80.0u"MPa", 500.0u"kJ/kg")

@test SpecificS(80.0, 0.611_041_229E3) ≈ SpecificS_Ph(80.0, 1500.0)
@test SpecificS(80.0u"MPa", 0.611_041_229E3u"K") ≈ SpecificS_Ph(80.0u"MPa", 1500.0u"kJ/kg")

@test SpecificCP(3.0, 0.391_798_509E3) ≈ SpecificCP_Ph(3.0, 500.0)
@test SpecificCP(3.0u"MPa", 0.391_798_509E3u"K") ≈ SpecificCP_Ph(3.0u"MPa", 500.0u"kJ/kg")

@test SpecificCP(80.0, 0.378_108_626E3) ≈ SpecificCP_Ph(80.0, 500.0)
@test SpecificCP(80.0u"MPa", 0.378_108_626E3u"K") ≈ SpecificCP_Ph(80.0u"MPa", 500.0u"kJ/kg")

@test SpecificCP(80.0, 0.611_041_229E3) ≈ SpecificCP_Ph(80.0, 1500.0)
@test SpecificCP(80.0u"MPa", 0.611_041_229E3u"K") ≈ SpecificCP_Ph(80.0u"MPa", 1500.0u"kJ/kg")

@test SpecificCV(3.0, 0.391_798_509E3) ≈ SpecificCV_Ph(3.0, 500.0)
@test SpecificCV(3.0u"MPa", 0.391_798_509E3u"K") ≈ SpecificCV_Ph(3.0u"MPa", 500.0u"kJ/kg")

@test SpecificCV(80.0, 0.378_108_626E3) ≈ SpecificCV_Ph(80.0, 500.0)
@test SpecificCV(80.0u"MPa", 0.378_108_626E3u"K") ≈ SpecificCV_Ph(80.0u"MPa", 500.0u"kJ/kg")

@test SpecificCV(80.0, 0.611_041_229E3) ≈ SpecificCV_Ph(80.0, 1500.0)
@test SpecificCV(80.0u"MPa", 0.611_041_229E3u"K") ≈ SpecificCV_Ph(80.0u"MPa", 1500.0u"kJ/kg")

@test SpeedOfSound(3.0, 0.391_798_509E3) ≈ SpeedOfSound_Ph(3.0, 500.0)
@test SpeedOfSound(3.0u"MPa", 0.391_798_509E3u"K") ≈ SpeedOfSound_Ph(3.0u"MPa", 500.0u"kJ/kg")

@test SpeedOfSound(80.0, 0.378_108_626E3) ≈ SpeedOfSound_Ph(80.0, 500.0)
@test SpeedOfSound(80.0u"MPa", 0.378_108_626E3u"K") ≈ SpeedOfSound_Ph(80.0u"MPa", 500.0u"kJ/kg")

@test SpeedOfSound(80.0, 0.611_041_229E3) ≈ SpeedOfSound_Ph(80.0, 1500.0)
@test SpeedOfSound(80.0u"MPa", 0.611_041_229E3u"K") ≈ SpeedOfSound_Ph(80.0u"MPa", 1500.0u"kJ/kg")


text:
using ClearStacktrace
using Test

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

module ModuleC
    func_c() = error("Error")
end

module ModuleB
    import ..ModuleC
    func_b(::String, ::Rational) = ModuleC.func_c()
end

module ModuleA
    import ..ModuleB
    func_a(::Int) = ModuleB.func_b("a", 1//3)
end

ModuleA.func_a(1)

ext: jl


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
struct combo_item
    p::Clonglong
    w::Clonglong
    x::Cint
    i::Cint
end

function solve_mono(pb::mono_problem,lb=0,ub=0)

    if isempty(variables(pb))
        return mono_solution(pb, BitArray(0), 0, 0)
    end

    #If all items fit in the KP (which happens sometimes when doing variable reductions),
    #combo throws a ReadOnlyMemoryError()
    if sum(x -> pb.w[x], variables(pb)) + pb.ω <= pb.c
        return mono_solution(pb, trues(size(pb)), sum(x -> pb.p[x], variables(pb)), sum(x -> pb.w[x], variables(pb)))
    end

    #@assert all(x -> x>=0, pb.p) "$(pb.p)"
    #@assert all(x -> x>0, pb.w)

    items = [combo_item(pb.p[j], pb.w[j], 0, i) for (i,j) in enumerate(variables(pb))]

    z = ccall((:solve,libcomboPath),
    Clonglong,
    (Ref{combo_item}, Cint, Clonglong, Clonglong, Clonglong),
    items, size(pb), pb.c - pb.ω, lb, ub)

    # #@assert lb <= z <= ub

    vars = falses(size(pb))
    for it in items
        it.x==1 && (vars[it.i]=true)
    end

    #If combo can't find a solution better than the lower bound,
    #it returns z = lb with an empty solution...
    #in that case, we recalculate z
    if z == lb
        z = dot(pb.p, vars, pb.variables)
    end

    #@assert z == dot(pb.p, vars, pb.variables) "$z != $(dot(pb.p, vars, pb.variables)) $lb $ub"

    return mono_solution(pb, vars, z, dot(pb.w, vars, pb.variables))

end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
using Pkg
using Test

testcode = raw"""
module MyPkg
using Cassette
Cassette.@context Ctx
const mypass = Cassette.@pass (ctx, ref) -> ref.code_info
end # module
"""

mktempdir() do dir
# for debugging use:
# let dir = mktempdir()
# @show dir
    cd(dir) do
        Pkg.generate("MyPkg")
        open(joinpath("MyPkg", "src", "MyPkg.jl"), "w") do io
            write(io, testcode)
        end
        Pkg.activate("MyPkg")
        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, ".."))) # add Cassette

        run(pipeline(`$(Base.julia_cmd()) --project=./MyPkg -e "using MyPkg"`, stderr="errs.log"))
        errs = read("errs.log", String)
        @testset "Precompilation" begin
            @test !occursin("WARNING: Method definition overdub", errs)
        end
    end
end


text:
"""Compute A' * Filter( A( LoG * f ) ) from Projections

where
LoG : Laplacian of Gaussian filter
f : object
A : Radon transform

Ref:

[1] 1992_Srinivasa,Ramakrishnan,Rajgopal_Detection_of_edges_from_projections IEEE_Transactions_on_Medical_Imaging
"""

using TomoForward
using XfromProjections
using SparseArrays
using PyPlot

include("./edge_from_proj.jl")
# img = imread("test_data/shepplogan512.png")[:,:,1]

img = zeros(100, 100)
img[20:50, 30:40] .= 1

nangles = 100
detcount = Int(floor(size(img,1)*1.4))
proj_geom = ProjGeom(1.0, detcount, LinRange(0,pi,nangles+1)[1:nangles])

A = fp_op_parallel2d_strip(proj_geom, size(img, 1), size(img, 2))
@time p = A * vec(img);
p = reshape(Array(p), (:, detcount));
q = filter_proj(p)

p_log = radon_log(q, 1.0)

@time LoG = A' * vec(p_log)
fbp = A' * vec(q)

LoG_img = reshape(LoG, size(img))
fbp_img = reshape(fbp, size(img))

ax00 = plt.subplot2grid((2,3), (0,0))
ax01 = plt.subplot2grid((2,3), (0,1))
ax02 = plt.subplot2grid((2,3), (0,2))
ax10 = plt.subplot2grid((2,3), (1,0))
ax11 = plt.subplot2grid((2,3), (1,1))
ax12 = plt.subplot2grid((2,3), (1,2))
ax00.imshow(p'); ax00.set_title("sinogram p")
ax01.imshow(fbp_img); ax01.set_title("fbp")
# ax02.imshow(bp_img_strip); ax02.set_title("bp_strip")
ax10.imshow(p_log); ax10.set_title("p * R(LoG)")
ax11.imshow(LoG_img, cmap="coolwarm"); ax11.set_title("LoG")
# ax12.imshow(fbp_img_strip); ax12.set_title("fbp_strip")
plt.show()

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
#!/usr/bin/julia

# Trizen
# Date: 27 August 2016
# https://github.com/trizen

# Generalized implementation of Knuth's up-arrow hyperoperation (modulo some n).

# See also:
#   https://en.wikipedia.org/wiki/Knuth%27s_up-arrow_notation

using Memoize
using Printf

const modulo = 10^3

@memoize function hyper1(n::Int64, k::Int64)
    powermod(n, k, modulo)
end

@memoize function hyper2(n::Int64, k::Int64)
    k <= 1 && return n
    hyper1(n, hyper2(n, k - 1))
end

@memoize function hyper3(n::Int64, k::Int64)
    k <= 1 && return n
    hyper2(n, hyper3(n, k - 1))
end

@memoize function hyper4(n::Int64, k::Int64)
    k <= 1 && return n
    hyper3(n, hyper4(n, k - 1))
end

@memoize function knuth(k::Int64, n::Int64, g::Int64)

    k %= modulo
    g %= modulo

    n >= 1 && g == 0 && return 1

    n == 0 && return ((k * g) % modulo)
    n == 1 && return (hyper1(k, g))
    n == 2 && return (hyper2(k, g))
    n == 3 && return (hyper3(k, g))
    n == 4 && return (hyper4(k, g))

    knuth(k, n - 1, knuth(k, n, g - 1))
end

for i in 0:6
    x = 1 + trunc(Int64, rand()*100)
    y = 1 + trunc(Int64, rand()*100)

    n = knuth(x, i, y)
    @printf("%5s %10s %5s = %5s   (mod %s)\n", x, repeat("^", i), y, n, modulo)
end

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
struct LineData{D, C, LP, SP, DP}
    line_params::LP
    swept_params::SP
    data::D
    data_pms::DP
    c::C
end

Base.show(io::IO, ld::LineData) =
    print(io, "LineData(", ld.line_params, ", ", ld.swept_params, ", ", ld.c, ")")

struct DataCollection{LD}
    data::Vector{LD}
end

Base.keys(dc::DataCollection) = keys(dc.data)

Base.show(io::IO, ld::DataCollection) =
    print(io, "DataCollection(LineParams: $(keys(ld[1].line_params)), SweptParams: $(keys(ld[1].swept_params)))")

Base.getindex(dc::DataCollection, idx) = if idx isa AbstractArray
    DataCollection(dc.data[idx])
else
    dc.data[idx]
end

Base.eltype(dc::DataCollection) = eltype(dc.data)

function load_runs(ic, get_data)
    tmp = get_data(FileIO.load(joinpath(ic[1].folder_str, "results.jld2")))
    d = typeof(tmp)[]

    diff_dict = diff(ic)
    # pms_template = NamedTuple(Symbol(k)=>diff_dict[k][1] for k ∈ keys(diff_dict))
    pms_template = @static if VERSION <= v"1.5"
        (; Symbol(k)=>diff_dict[k][1] for k ∈ keys(diff_dict))
    else
        NamedTuple(Symbol(k)=>diff_dict[k][1] for k ∈ keys(diff_dict))
    end
    pms = typeof(pms_template)[]

    for (idx, item) ∈ enumerate(ic)
        push!(d, get_data(FileIO.load(joinpath(item.folder_str, "results.jld2"))))
        @static if VERSION <= v"1.5"
            push!(pms, (; Symbol(k)=>item.parsed_args[k] for k ∈ keys(diff_dict)))
        else
            push!(pms, NamedTuple(Symbol(k)=>item.parsed_args[k] for k ∈ keys(diff_dict)))
        end
    end
    d, pms
end

function get_comp_function(comp::Symbol)
    if comp == :max
        findmax
    elseif comp == :min
        findmin
    end
end

function get_data_function(get_data_sym::Symbol)
end

function find_best_params(ic,
	                  param_keys,
		          comp,
		          get_comp_data,
		          get_data=get_comp_data)

    if comp isa Symbol
        comp = get_comp_function(comp)
    end

    ic_diff = diff(ic)
    params = if param_keys isa String
        ic_diff[param_keys]
    else
        collect(Iterators.product([ic_diff[k] for k ∈ param_keys]...))
    end

    s = zeros(length(params))
    for (p_idx, p) ∈ enumerate(params)
        sub_ic = if param_keys isa String
            search(ic, Dict(param_keys=>p))
        else
            search(ic, Dict(param_keys[i]=>p[i] for i ∈ 1:length(p)))
        end

        s[p_idx] = mean(load_runs(sub_ic, get_comp_data)[1])
    end

    v, idx = comp(s)

    best_ic = if param_keys isa String
	search(
	    ic,
	    Dict(param_keys=>params[idx]))
    else
	search(
	    ic,
	    Dict(param_keys[i]=>params[idx][i] for i ∈ 1:length(params[idx])))

    end

    data, data_pms = load_runs(best_ic, get_data)

    data, data_pms, v, params[idx]
end

function get_line_data_for(
    ic::ItemCollection,
    line_keys,
    param_keys;
    comp,
    get_comp_data,
    get_data)

    ic_diff = diff(ic)
    params = if line_keys isa String
        ic_diff[line_keys]
    else
        collect(Iterators.product([ic_diff[k] for k ∈ line_keys]...))
    end

    strg = LineData[]
    lk = ReentrantLock()

    Threads.@threads for p_idx ∈ 1:length(params)
	p = params[p_idx]

	sub_ic = if line_keys isa String
            search(ic, Dict(line_keys=>p))
        else
            search(ic, Dict(line_keys[i]=>p[i] for i ∈ 1:length(p)))
        end

        if !isempty(sub_ic)
	    d, dp, c, ps = find_best_params(
	        sub_ic,
	        param_keys,
                comp,
	        get_comp_data,
	        get_data)

            ld = if line_keys isa String
                LineData(Dict(line_keys=>params[p_idx]), ps, d, dp, c)
            else
                LineData(Dict(line_keys[i]=>params[p_idx][i] for i in 1:length(line_keys)), ps, d, dp, c)
            end

            lock(lk) do
                push!(strg, ld)
            end

        end
    end
    DataCollection(convert(Vector{typeof(strg[1])}, strg))
end

# function get_data_frame_for(
#     ic::ItemCollection,
#     line_keys,
#     param_keys;
#     comp,
#     get_comp_data,
#     get_data)

#     ic_diff = diff(ic)
#     params = if line_keys isa String
#         ic_diff[line_keys]
#     else
#         collect(Iterators.product([ic_diff[k] for k ∈ line_keys]...))
#     end

#     strg = LineData[]

#     Threads.@threads for p_idx ∈ 1:length(params)
# 	p = params[p_idx]

# 	sub_ic = if line_keys isa String
#             search(ic, Dict(line_keys=>p))
#         else
#             search(ic, Dict(line_keys[i]=>p[i] for i ∈ 1:length(p)))
#         end

#         if !isempty(sub_ic)
# 	    d, dp, c, ps = find_best_params(
# 	        sub_ic,
# 	        param_keys,
#                 comp,
# 	        get_comp_data,
# 	        get_data)

#             if line_keys isa String
#                 push!(strg, LineData(Dict(line_keys=>params[p_idx]), ps, d, dp, c))
#             else
#                 push!(strg, LineData(Dict(line_keys[i]=>params[p_idx][i] for i in 1:length(line_keys)), ps, d, dp, c))
#             end

#         end
#     end
#     DataCollection(convert(Vector{typeof(strg[1])}, strg))
# end


text:
function unit_test_client(verbose=true)
    bits_tab = Vector{UInt8}(make_bitvector(UT_BITS_TAB, UT_BITS_NB))
    input_bits_tab = Vector{UInt8}(make_bitvector(UT_INPUT_BITS_TAB, UT_INPUT_BITS_NB))

    errno = Libc.errno
    ctx_cl = TcpContext("127.0.0.1", 1502)
    modbus_set_debug(ctx_cl, verbose)
    modbus_flush(ctx_cl)
    #modbus_set_error_recovery(ctx_cl, Int(MODBUS_ERROR_RECOVERY_LINK) |
    #    Int(MODBUS_ERROR_RECOVERY_PROTOCOL))
   @testset "LibModbus.jl tests" begin
    rc,old_response_to_sec,old_response_to_usec = modbus_get_response_timeout(ctx_cl)
    modbus_connect(ctx_cl) >= 0 || error("connection failed")
    @testset "1/1 No response timeout modification on connect" begin
        rc,new_response_to_sec,new_response_to_usec = modbus_get_response_timeout(ctx_cl)
        @test ((old_response_to_sec,old_response_to_usec)
               == (new_response_to_sec,new_response_to_usec))
    end
        @testset "Single coil bits" begin
            verbose && println("1/2 modbus_write_bit: ")
            rc = modbus_write_bit(ctx_cl, UT_BITS_ADDRESS, 0x1)
            @test rc == 1;modbus_flush(ctx_cl)
            verbose && println("2/2 modbus_read_bits: ")
            rc,rp_bits = modbus_read_bits(ctx_cl, UT_BITS_ADDRESS, 0x1)
            @test rc == 1
            @test length(rp_bits) > 0 && rp_bits[1] == 1
        end

        @testset "Multiple coil bits" begin
            verbose && println("1/2 modbus_write_bits: ")
            rc = modbus_write_bits(ctx_cl, UT_BITS_ADDRESS, bits_tab)
            @test rc == UT_BITS_NB
            verbose && println("2/2 modbus_read_bits: ")
            rc,rp_bits = modbus_read_bits(ctx_cl, UT_BITS_ADDRESS, UT_BITS_NB)
            @test rc == UT_BITS_NB
            @test rp_bits == bits_tab
        end

        @testset "Discrete inputs" begin
            verbose && println("1/1 modbus_read_input_bits: ")
            rc,rp_bits = modbus_read_input_bits(ctx_cl, UT_INPUT_BITS_ADDRESS,
                                                UT_INPUT_BITS_NB)
            @test rc == UT_INPUT_BITS_NB
            @test rp_bits == input_bits_tab
        end

        @testset "Holding registers" begin
            #  Single register
            verbose && println("1/2 modbus_write_register: ")
            rc = modbus_write_register(ctx_cl, UT_REGISTERS_ADDRESS, 0x1234)
            @test rc == 1
            verbose && println("2/2 modbus_read_registers: ")
            rc,rp_regs = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS, 1)
            @test rc == 1
            @test rp_regs[] == 0x1234
            # Many registers
            verbose && println("1/5 modbus_write_registers: ")
            rc = modbus_write_registers(ctx_cl, UT_REGISTERS_ADDRESS, UT_REGISTERS_TAB)
            @test rc == UT_REGISTERS_NB
            verbose && println("2/5 modbus_read_registers: ")
            rc,rp_regs = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS,
                                               UT_REGISTERS_NB)
            @test rc == UT_REGISTERS_NB
            @test rp_regs == UT_REGISTERS_TAB
            verbose && println("3/5 modbus_read_registers (0): ")
            rc,_ = @test_logs((:warn, r"Illegal data value"),
                              modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS, 0))
            errn = errno()
            @test rc == -1 && errn == EMBXILVAL
            rp_regs = zeros(UInt16, UT_REGISTERS_NB - 1)
            verbose && println("4/5 modbus_write_and_read_registers: ")
            rc,rp_regs = modbus_write_and_read_registers(ctx_cl,
                                                         UT_REGISTERS_ADDRESS + 1,
                                                         rp_regs,
                                                         UT_REGISTERS_ADDRESS,
                                                         UT_REGISTERS_NB)
            @test rc == UT_REGISTERS_NB
            @test rp_regs[1] == UT_REGISTERS_TAB[1]
        end

        @testset "Input registers" begin
            verbose && println("1/1 modbus_read_input_registers: ")
            rc,rp_reg = modbus_read_input_registers(ctx_cl,
                                                    UT_INPUT_REGISTERS_ADDRESS,
                                                    UT_INPUT_REGISTERS_NB)
            @test rc == UT_INPUT_REGISTERS_NB

            verbose && println("1/1 Write mask: ")
            rc = modbus_write_register(ctx_cl, UT_REGISTERS_ADDRESS, 0x12)
            rc = modbus_mask_write_register(ctx_cl, UT_REGISTERS_ADDRESS, 0x00F2, 0x0025)
            @test rc != -1
            rc,rp_regs = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS, 1)
            @test rp_regs[1] == 0x17
        end

        @testset "Illegal data address" begin
            ## The mapping begins at the defined addresses and ends at address +
            ## nb_points so these addresses are not valid.
            verbose && println("* modbus_read_bits (0): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_read_bits(ctx_cl, 0, 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_read_bits (max): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_read_bits(ctx_cl, UT_BITS_ADDRESS, UT_BITS_NB + 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_read_input_bits (0): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_read_input_bits(ctx_cl, 0, 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_read_input_bits (max): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_read_input_bits(ctx_cl, UT_INPUT_BITS_ADDRESS,
                                                      UT_INPUT_BITS_NB + 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_read_registers (0): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_read_registers(ctx_cl, 0, 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_read_registers (max): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS,
                                                     UT_REGISTERS_NB_MAX + 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_read_input_registers (0): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_read_input_registers(ctx_cl, 0, 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_read_input_registers (max): ")
            rc,_ =
                @test_logs((:warn, r"Illegal data address"),
                           modbus_read_input_registers(ctx_cl, UT_INPUT_REGISTERS_ADDRESS,
                                                       UT_INPUT_REGISTERS_NB + 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_bit (0): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_bit(ctx_cl, 0, 0x01))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_bit (max): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_bit(ctx_cl, UT_BITS_ADDRESS + UT_BITS_NB, 0x01))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_coils (0): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_bits(ctx_cl, 0, [0x01]))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_coils (max): ")
            rp_bits = zeros(UInt8, UT_BITS_NB)
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_bits(ctx_cl, UT_BITS_ADDRESS + UT_BITS_NB, rp_bits))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_register (0): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_register(ctx_cl, 0x0, 0x01))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_register (max): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_register(ctx_cl,
                                                  UT_REGISTERS_ADDRESS + UT_REGISTERS_NB_MAX,
                                                  0x01))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_registers (0): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_registers(ctx_cl, 0, UInt16[1]))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_registers (max): ")
            rp_regs = zeros(UInt16, UT_REGISTERS_NB)
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_write_registers(ctx_cl,
                                                   UT_REGISTERS_ADDRESS +
                                                       UT_REGISTERS_NB_MAX,
                                                   rp_regs))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_mask_write_registers (0): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_mask_write_register(ctx_cl, 0, 0x00F2, 0x0025))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_mask_write_registers (max): ")
            rc = @test_logs((:warn, r"Illegal data address"),
                            modbus_mask_write_register(ctx_cl,
                                                       UT_REGISTERS_ADDRESS +
                                                           UT_REGISTERS_NB_MAX,
                                                       0x00F2, 0x0025))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_and_read_registers (0): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_write_and_read_registers(ctx_cl, 0, rp_regs[1:1],
                                                               0, 1))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD

            verbose && println("* modbus_write_and_read_registers (max): ")
            rc,_ =  @test_logs((:warn, r"Illegal data address"),
                               modbus_write_and_read_registers(ctx_cl,
                                                               UT_REGISTERS_ADDRESS +
                                                                   UT_REGISTERS_NB_MAX,
                                                               rp_regs,
                                                               UT_REGISTERS_ADDRESS +
                                                                   UT_REGISTERS_NB_MAX,
                                                               UT_REGISTERS_NB))
            errn = errno()
            @test rc == -1 && errn == EMBXILADD
        end

        @testset "Too many data" begin
            verbose && println("* modbus_read_bits: ")
            rc, _ = @test_logs((:warn, r"Too many data"),
                               modbus_read_bits(ctx_cl, UT_BITS_ADDRESS,
                                                MODBUS_MAX_READ_BITS + 1))
            errn = errno()
            @test rc == -1 && errn == EMBMDATA

            verbose && println("* modbus_read_input_bits: ")
            rc,_ = @test_logs((:warn, r"Too many data"),
                              modbus_read_input_bits(ctx_cl, UT_INPUT_BITS_ADDRESS,
                                                     MODBUS_MAX_READ_BITS + 1))
            errn = errno()
            @test rc == -1 && errn == EMBMDATA

            verbose && println("* modbus_read_registers: ")
            rc,_ = @test_logs((:warn, r"Too many data"),
                              modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS,
                                                    MODBUS_MAX_READ_REGISTERS + 1))
            errn = errno()
            @test rc == -1 && errn == EMBMDATA

            verbose && println("* modbus_read_input_registers: ")
            rc,_ = @test_logs((:warn, r"Too many data"),
                              modbus_read_input_registers(ctx_cl, UT_INPUT_REGISTERS_ADDRESS,
                                                          MODBUS_MAX_READ_REGISTERS + 1))
            errn = errno()
            @test rc == -1 && errn == EMBMDATA

            verbose && println("* modbus_write_bits: ")
            rp_bits = zeros(UInt8,  MODBUS_MAX_WRITE_BITS + 1)
            rc = @test_logs((:warn, r"Too many data"),
                            modbus_write_bits(ctx_cl, UT_BITS_ADDRESS, rp_bits))
            errn = errno()
            @test rc == -1 && errn == EMBMDATA

            verbose && println("* modbus_write_registers: ")
            rp_regs = zeros(UInt16, MODBUS_MAX_WRITE_REGISTERS + 1)
            rc = @test_logs((:warn, r"Too many data"),
                            modbus_write_registers(ctx_cl, UT_REGISTERS_ADDRESS, rp_regs))
            errn = errno()
            @test rc == -1 && errn == EMBMDATA
        end

        @testset "Slave reply" begin

            old_slave = modbus_get_slave(ctx_cl)
            verbose && println("1/3 Response from slave $(INVALID_SERVER_ID): ")
            modbus_set_slave(ctx_cl, INVALID_SERVER_ID)
            rc,_ = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS, UT_REGISTERS_NB)
            @test rc == UT_REGISTERS_NB

            rc = modbus_set_slave(ctx_cl, MODBUS_BROADCAST_ADDRESS)
            @test rc != -1

            verbose && println("2/3 Reply after a query with unit id == 0: ")
            rc,_ = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS, UT_REGISTERS_NB)
            @test rc == UT_REGISTERS_NB

            # Restore slave
            modbus_set_slave(ctx_cl, old_slave)
            verbose && println("3/3 Response with an invalid TID or slave: ")
            rc,_= @test_logs((:warn, r"Invalid data"),
                             modbus_read_registers(ctx_cl,
                                                   UT_REGISTERS_ADDRESS_INVALID_TID_OR_SLAVE,
                                                   1))
            errn = errno()
            @test rc == -1  && errn == EMBBADDATA
            modbus_flush(ctx_cl)
            verbose && println("2/2 Report slave ID:")
            #
            rc,rp_bits = modbus_report_slave_id(ctx_cl, NB_REPORT_SLAVE_ID)
            # Run status indicator is ON
            @test rc > 1 && rp_bits[2] == 0xFF

            # Save original timeout
            rc,old_response_to_sec,old_response_to_usec =
                modbus_get_response_timeout(ctx_cl)
            rc,old_byte_to_sec,old_byte_to_usec = modbus_get_byte_timeout(ctx_cl)
            einval_reg = Regex("$(Libc.strerror(Libc.EINVAL))")
            verbose && println("1/6 Invalid response timeout (zero): ")
            rc =  @test_logs((:warn, einval_reg),
                             modbus_set_response_timeout(ctx_cl, 0, 0))
            errn = errno()
            @test rc == -1 && errn == Libc.EINVAL

            verbose && println("2/6 Invalid response timeout (too large us): ")
            rc = @test_logs((:warn, einval_reg),
                            modbus_set_response_timeout(ctx_cl, 0, 1000000))
            errn = errno()
            @test rc == -1 && errn == Libc.EINVAL

            verbose && println("3/6 Invalid byte timeout (too large us): ")
            rc = @test_logs((:warn, einval_reg),
                            modbus_set_byte_timeout(ctx_cl, 0, 1000000))
            errn = errno()
            @test rc == -1 && errn == Libc.EINVAL

            # # verbose && println("4/6 1us response timeout: ")
            # etimedout_reg = Regex("$(Libc.strerror(Libc.ETIMEDOUT))")
            # modbus_set_response_timeout(ctx_cl, 0, 1)
            # rc,_ = @test_logs((:warn, etimedout_reg),
            #                   modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS,
            #                                         UT_REGISTERS_NB))
            # errn = errno()
            # @test rc == -1 && errn == Libc.ETIMEDOUT

            # # # A wait and flush operation is done by the error recovery code of
            # # # libmodbus but after a sleep of current response timeout
            # # # so 0 can be too short!
            # sleep(old_response_to_sec + 1e-6*old_response_to_usec)
            # modbus_flush(ctx_cl)

            # # # Trigger a special behaviour on server to wait for 0.5 second before
            # # # replying whereas allowed timeout is 0.2 second
            # verbose && println("5/6 Too short response timeout (0.2s < 0.5s): ")
            # modbus_set_response_timeout(ctx_cl, 0, 200000)
            # rc,_ =  @test_logs((:warn, etimedout_reg),
            #                    modbus_read_registers(ctx_cl,
            #                                          UT_REGISTERS_ADDRESS_SLEEP_500_MS, 1))
            # errn = errno()
            # @test rc == -1 && errn == Libc.ETIMEDOUT

            # # Wait for reply (0.2 + 0.4 > 0.5 s) and flush before continue
            # sleep(0.4)
            # modbus_flush(ctx_cl)

            # verbose && println("6/6 Adequate response timeout (0.6s > 0.5s): ")
            # modbus_set_response_timeout(ctx_cl, 0, 600000)
            # rc,_ = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS_SLEEP_500_MS, 1)
            # @test rc == 1
            # modbus_flush(ctx_cl)
            # # Disable the byte timeout.
            #   The full response must be available in the 600ms interval
            # verbose && println("7/7 Disable byte timeout: ")
            # rc=modbus_set_byte_timeout(ctx_cl, 0, 0);println("byte timeout disable $(rc)")
            # rc,_ = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS_SLEEP_500_MS, 1)
            # @test rc == 1
            # # # Restore original response and byte timeout
            # modbus_set_response_timeout(ctx_cl, old_response_to_sec, old_response_to_usec)
            # modbus_set_byte_timeout(ctx_cl, old_byte_to_sec, old_byte_to_usec)
            # sleep(0.4)
            # modbus_flush(ctx_cl)
            # # The test server is only able to test byte timeouts with the TCP
            # # backend
            # # Timeout of 3ms between bytes
            # verbose && println("1/2 Too small byte timeout (3ms < 5ms): ")
            # modbus_set_byte_timeout(ctx_cl, 0, 3000)
            # rc,_ = @test_logs((:warn, etimedout_reg),
            #                   modbus_read_registers(ctx_cl,
            #                                         UT_REGISTERS_ADDRESS_BYTE_SLEEP_5_MS,
            #                                         1))
            # errn = errno()
            # @test rc == -1 && errn == Libc.ETIMEDOUT

            # # Wait remaing bytes before flushing
            # sleep(11 * 0.005)
            # modbus_flush(ctx_cl)

            # # # Timeout of 7ms between bytes
            # modbus_set_byte_timeout(ctx_cl, 0, 7000)
            # verbose && println("2/2 Adapted byte timeout (7ms > 5ms): ")
            # rc,_ = modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS_BYTE_SLEEP_5_MS, 1)
            # @test rc == 1

            # # Restore original timeouts
            modbus_set_response_timeout(ctx_cl, old_response_to_sec, old_response_to_usec)
            modbus_set_byte_timeout(ctx_cl, old_byte_to_sec, old_byte_to_usec)
            sleep(0.4)
            modbus_flush(ctx_cl)
        end
        @testset "Bad response" begin
            verbose && println("* modbus_read_registers with invalid data: ")
            rc,_= @test_logs((:warn, r"Invalid data"),
                             modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS,
                                                   UT_REGISTERS_NB_SPECIAL))
            errn = errno()
            @test rc == -1 && errn == EMBBADDATA
            verbose && println("* modbus_read_registers at special address: ")
            rc,_ =  @test_logs((:warn, r"Slave device or server is busy"),
                               modbus_read_registers(ctx_cl, UT_REGISTERS_ADDRESS_SPECIAL,
                                                     UT_REGISTERS_NB))
            errn = errno()
            @test rc == -1 && errn == EMBXSBUSY
        end
        modbus_close(ctx_cl)
        modbus_free!(ctx_cl)
        @testset "Invalid initialization" begin
            einval_reg = Regex("$(Libc.strerror(Libc.EINVAL))")
            ctx_cl = @test_logs((:warn, einval_reg), RtuContext("", 1, :none, 0, 0))
            errn = errno()
            @test !modbus_context_valid(ctx_cl) && errn == Libc.EINVAL

            ctx = @test_logs((:warn, einval_reg), RtuContext("/dev/dummy", 0, :none, 0, 0))
            errn = errno()
            @test !modbus_context_valid(ctx_cl) && errn == Libc.EINVAL
        end
    end
    nothing
end


text:
"""
## glob_to_loc!

This subroutine transforms the global end reactions and
moments into the local system (2- or 3-d). Called from hinge!().

### Function
```julia
glob_to_loc!(loc, glob, gamma, coord)
```

### Arguments
```julia
* loc::Vector{Float64}       : Local force and momemts (Updated)
* glob::Vector{Float64}      : Globale forces and moments
* gamma::Float64             : Element orientation angle (3D)
* coord::Matrix{Float64}     : Nodal coordinates
```
"""
function glob_to_loc!(loc::Vector{Float64}, glob::Vector{Float64},
  gamma::Float64, coord::LinearAlgebra.Adjoint{Float64,Array{Float64,2}})
  #
  # This subroutine transforms the global end reactions and
  # moments into the local system (2- or 3-d).
  #
  # Updates loc[:]
  #

  ndim = size(coord, 2)
  if ndim == 2
    x1 = coord[1, 1]
    y1 = coord[1, 2]
    x2 = coord[2, 1]
    y2 = coord[2, 2]
    ell = sqrt((x2-x1)^2+(y2-y1)^2)
    cg = (x2-x1)/ell
    sg = (y2-y1)/ell
    loc[1] = cg*glob[1]+sg*glob[2]
    loc[2] = cg*glob[2]-sg*glob[1]
    loc[3] = glob[3]
    loc[4] = cg*glob[4]+sg*glob[5]
    loc[5] = cg*glob[5]-sg*glob[4]
    loc[6] = glob[6]
  elseif ndim == 3
    t = zeros(12,12)
    r0 = zeros(3, 3)
    x1 = coord[1, 1]
    y1 = coord[1, 2]
    z1 = coord[1, 3]
    x2 = coord[2, 1]
    y2 = coord[2, 2]
    z2 = coord[2, 3]
    xl = x2 - x1
    yl = y2 - y1
    zl = z2 - z1
    ell = sqrt(xl*xl+yl*yl+zl*zl)
    gamrad = gamma*pi/180.0
    cg = cos(gamrad)
    sg = sin(gamrad)
    den = ell*sqrt(xl*xl+zl*zl)
    if den !== 0.0
      r0[1,1] = xl/ell
      r0[1,2] = yl/ell
      r0[1,3] = zl/ell
      r0[2,1] = (-xl*yl*cg-ell*zl*sg)/den
      r0[2,2] = den*cg/(ell*ell)
      r0[2,3] = (-yl*zl*cg+ell*xl*sg)/den
      r0[3,1] = (xl*yl*sg-ell*zl*cg)/den
      r0[3,2] = -den*sg/(ell*ell)
      r0[3,3] = (yl*zl*sg+ell*xl*cg)/den
    else
      r0[1,1] = 0.0
      r0[1,3] = 0.0
      r0[2,2] = 0.0
      r0[3,2] = 0.0
      r0[1,2] = 1.0
      r0[2,1] = -cg
      r0[3,3] = cg
      r0[2,3] = sg
      r0[3,1] = sg
    end
    for i in 1:3
      for j in 1:3
        x = r0[i,j]
        for  k in 0:3:9
          t[i+k,j+k] = x
        end
      end
    end
    for i in 1:12
      sum = 0.0
      for j in 1:12
        sum = sum+t[i,j]*glob[j]
      end
      loc[i] = sum
    end
  end
end


text:
using MimiFUND
using Mimi
using Test
using DataFrames
using CSVFiles
using Query

@testset "fund" begin

#------------------------------------------------------------------------------
#   1. Run the whole model
#------------------------------------------------------------------------------

@testset "fund-model" begin

#default model exported by fund module
default_nsteps = 1050
m = MimiFUND.get_model()
run(m)
@test Mimi.time_labels(m) == collect(1950:1:1950+default_nsteps)

#default model created by MimiFUND.get_model()
m1 = MimiFUND.get_model()
run(m1)
@test Mimi.time_labels(m1) == collect(1950:1:1950+default_nsteps)

# use optional `nsteps` arg for MimiFUND.get_model()
@test_throws ErrorException m2 = MimiFUND.get_model(nsteps = 2000) # should error because it's longer than the default without providing different parameters
new_nsteps = 350
m2 = MimiFUND.get_model(nsteps = new_nsteps)
@test length(Mimi.dimension(m2, :time)) == new_nsteps + 1
run(m2)
@test length(m2[:climatedynamics, :temp]) == new_nsteps + 1

end #fund-model testset

#------------------------------------------------------------------------------
#   2. Run tests to make sure integration version (Mimi v0.5.0)
#   values match Mimi 0.4.0 values
#------------------------------------------------------------------------------
@testset "test-integration" begin

m = MimiFUND.get_model()
run(m)

missingvalue = -999.999
err_number = 1.0e-9
err_array = 0.0

# for c in map(name, Mimi.compdefs(m)), v in Mimi.variable_names(m, c)
for c in Mimi.compdefs(m), v in Mimi.variable_names(m, nameof(c))

    # load data for comparison
    orig_name = c.comp_id.comp_name
    filename = joinpath(@__DIR__, "../contrib/validation_data_v040/$orig_name-$v.csv")
    results = m[nameof(c), v]

    df = load(filename) |> DataFrame
    if typeof(results) <: Number
        validation_results = df[1,1]
        @test results ≈ validation_results atol = err_number #slight imprecision with these values due to rounding

    else
        validation_results = Matrix(df)

        #replace missings with missingvalue so they can be compared
        results[ismissing.(results)] .= missingvalue
        validation_results[ismissing.(validation_results)] .= missingvalue

        #match dimensions
        if size(validation_results,1) == 1
            validation_results = validation_results'
        end

        @test results ≈ validation_results atol = err_array

    end
end

end #fund-integration testset

#------------------------------------------------------------------------------
# 3. Test marginal damages functions (test that each function does not error)
#------------------------------------------------------------------------------

@testset "test-compute_sc" begin

# Test functions from file "new_marginaldamages.jl"

# Test the compute_scco2 function with various keyword arguments
scc0 = MimiFUND.compute_sc(year = 2020, gas = :CO2)
scc1 = MimiFUND.compute_scco2(year = 2020)
@test scc1 isa Float64   # test that it's not missing or a NaN
@test scc0 == scc1
scc2 = MimiFUND.compute_scco2(year = 2020, last_year=2300)
@test scc2 < scc1  # test that shorter horizon made it smaller
scc3 = MimiFUND.compute_scco2(year = 2020, last_year=2300, equity_weights=true)
@test scc3 > scc2  # test that equity weights made itbigger
scc4 = MimiFUND.compute_scco2(year = 2020, last_year=2300, equity_weights=true, eta=.8, prtp=0.01)
@test scc4 > scc3   # test that lower eta and prtp make scc higher
scch4 = MimiFUND.compute_scch4(year = 2020)
@test scch4 > scc1   # test social cost of methane is higher

# Test with a modified model
m_high_cs = MimiFUND.get_model()
set_param!(m_high_cs, :climatesensitivity, 5)
scc6 = MimiFUND.compute_scco2(m_high_cs, year=2020, last_year=2300)
@test scc6 > scc2 # test that it's higher than the default because of a higher climate sensitivity

# Test get_marginal_model
mm = MimiFUND.get_marginal_model(year=2020, gas=:CH4)
run(mm)

# Test return_mm = true
result = MimiFUND.compute_sc(year=2050, return_mm = true)
@test result.sc isa Float64
@test result.mm isa Mimi.MarginalModel

# Test other gases
scch4 = MimiFUND.compute_scch4(year = 2020)
scn2o = MimiFUND.compute_scn2o(year = 2020)
scsf6 = MimiFUND.compute_scsf6(year = 2020)
@test scch4 < scn2o < scsf6

# Test that modifying the pulse_size keyword changes the values, but not by much
scch4_2 = MimiFUND.compute_scch4(year = 2020, pulse_size = 1e3)
scn2o_2 = MimiFUND.compute_scn2o(year = 2020, pulse_size = 1e3)
scsf6_2 = MimiFUND.compute_scsf6(year = 2020, pulse_size = 1e3)
@test scch4 != scch4_2
@test scch4 ≈ scch4_2 rtol = 1e-3
@test scn2o != scn2o_2
@test scn2o ≈ scn2o_2 rtol = 1e-3
@test scsf6 != scsf6_2
@test scsf6 ≈ scsf6_2 rtol = 1e-1 # TODO Lisa Rennels - ok to increase tolerance here?

# Test monte carlo simulation
scco2_values = MimiFUND.compute_sc(year = 2020, gas = :CO2, n = 10)
@test all(!isnan, scco2_values)

# Test the seed functionality
results1 = MimiFUND.compute_sc(gas = :CH4, year = 2020, n = 10, seed = 350)
results2 = MimiFUND.compute_sc(gas = :CH4, year = 2020, n = 10, seed = 350, trials_output_filename = "tmp_trials.csv")  # test that saving trial values does not affec the seed/sampling order
rm("tmp_trials.csv")    # remove the save trials data
@test results1 == results2
results_shorter = MimiFUND.compute_sc(gas = :CH4, year = 2020, last_year = 2300, n = 10, seed = 350)
@test results_shorter != results1

# Test old exported marginal damage function
md = MimiFUND.getmarginaldamages()

end #marginaldamages testset

#------------------------------------------------------------------------------
# 4. Run basic test of Marginal Damages and MCS functionality
#------------------------------------------------------------------------------

@testset "test-mcs" begin

# mcs
MimiFUND.run_fund_mcs(10)        # Run 10 trials of basic FUND MCS
MimiFUND.run_fund_scc_mcs(10)    # Run 10 trials of FUND MCS SCC calculations

end #test-mcs testset

#------------------------------------------------------------------------------
# 5. Validation tests for the exact values of social cost calculations
#------------------------------------------------------------------------------

@testset "sc-validation" begin

datadir = joinpath(@__DIR__, "SC validation data")
atol = 1e-2

# Test a subset of all validation configurations against the pre-saved values from MimiFUND v3.13.0

validation_results = load(joinpath(datadir, "deterministic_sc_values_v3-13-0.csv")) |> DataFrame

for spec in [
    (gas = :CO2, year = 2020, eta = 0., prtp = 0.03, equity_weights = true, equity_weights_normalization_region = 1, last_year = 3000, pulse_size = 1.),
    (gas = :CH4, year = 2055, eta = 1.45, prtp = 0.015, equity_weights = false, equity_weights_normalization_region = 0, last_year = 2300, pulse_size = 1e7),
    (gas = :N2O, year = 2020, eta = 1.45, prtp = 0.03, equity_weights = false, equity_weights_normalization_region = 0, last_year = 3000, pulse_size = 1e7),
    (gas = :SF6, year = 2055, eta = 0., prtp = 0.015, equity_weights = true, equity_weights_normalization_region = 10, last_year = 2300, pulse_size = 1.)
]
    # Compute the value for this set of keyword arguments
    sc = MimiFUND.compute_sc(gas = spec.gas, year = spec.year, eta = spec.eta, prtp = spec.prtp, equity_weights = spec.equity_weights, equity_weights_normalization_region = spec.equity_weights_normalization_region, last_year = spec.last_year, pulse_size = spec.pulse_size)

    # Get the pre-saved value for this set of keyword arguments
    validation_value = validation_results |>
        @filter(_.gas == string(spec.gas) && _.year == spec.year && _.eta == spec.eta && _.prtp == spec.prtp && _.equity_weights == string(spec.equity_weights) && _.equity_weights_normalization_region == spec.equity_weights_normalization_region && _.last_year == spec.last_year && _.pulse_size == spec.pulse_size) |>
        DataFrame

    @test sc ≈ validation_value[1, :SC] atol = atol
end

# Test Monte Carlo results with the same configuration and seed - note here we use
# FUND 3.11.7 since 3.13.0 only updated SC-N2O and SC-SF6
sc_mcs = MimiFUND.compute_sc(gas = :CO2, year = 2020, eta = 1.45, prtp = 0.015, n = 25, seed = 350)
validation_mcs = load(joinpath(datadir, "mcs_sc_values_v3-11-7.csv")) |> DataFrame
@test all(isapprox.(sc_mcs, validation_mcs[!, :SCCO2], atol = atol))

end

end #fund testset

nothing


text:
include(joinpath(@__DIR__, "../../rand_skewnormal.jl"))

function simulatedata(; Z, W, N::Vector{Int},
                      mus::Dict{Int, Vector{Float64}},
                      sig2=Vector{Float64},
                      skew::Float64,
                      sortLambda=false, seed=nothing,
                      propmissingscale=0.3,
                      expressed_not_missing=false,
                      beta=[-9.2, -2.3],  # linear missing mechanism
                      eps_mus_dist=Uniform(-.3, .3))
  if seed != nothing
    Random.seed!(seed)
  end

  J, K = size(Z)
  I = length(N)

  # Create dict of L
  L = Dict(0 => length(mus[0]),
           1 => length(mus[1]))

  @assert size(W) == (I, K)
  @assert length(sig2) == I
  @assert all(sig2 .> 0)
  @assert all(length(mus[key]) == L[key] for key in keys(L))
  @assert all(N .> 0)

  # Simulate cell phenotypes (lambda)
  lam = [begin
           rand(Categorical(W[i, :]), N[i])
         end for i in 1:I]

  if sortLambda
    lam = [sort(lami) for lami in lam]
  end

  # Generate y
  y = [zeros(Ni, J) for Ni in N]
  y_complete = deepcopy(y)

  z(i::Int, n::Int, j::Int) = Z[j, lam[i][n]]

  eps_mus = let
    if eps_mus_dist != nothing
      Dict(0 => rand(eps_mus_dist, I, J),
           1 => rand(eps_mus_dist, I, J))
    else
      Dict(0 => zeros(I, J),
           1 => zeros(I, J))
    end
  end

  function mu(i::Int, n::Int, j::Int)
    @assert L[0] == L[1] == 1
    z_inj = z(i, n, j)
    return mus[z_inj][1] + eps_mus[z_inj][i, j]
  end

  for i in 1:I
    sig_i = sqrt(sig2[i])
    for j in 1:J
      for n in 1:N[i]
        # y_complete[i][n, j] = rand(Normal(mu(i, n, j), sig_i))
        y_complete[i][n, j] = rand_skewnormal(mu(i, n, j), sig_i, skew)
      end
      # Set some y to be missing
      # Using a linear missing mechanism to rate missingness
      p_miss = [CytofRepFAM.Model.prob_miss(y_complete[i][n, j], beta)
                for n in 1:N[i]]

      # If I want to truly expressed cells to not be missing:
      if expressed_not_missing
        # If Z_inj == 1, then p_miss = 0
        # Otherwise, p_miss = p_miss
        p_miss .*= [1 - z(i, n, j) for n in 1:N[i]]
      end

      # If this is not used, then many positive markers may be missing too.
      prop_not_expressed = sum(W[i,:] .* (1 .- Z[j,:]))
      # If propmissingscale is 1, about all the non-expressed observations
      # will be missing. Best to set propmissingscale to about 0.3.
      prop_missing = propmissingscale * prop_not_expressed

      num_missing = Int(round(N[i] * prop_missing))
      idx_missing = Distributions.wsample(1:N[i], p_miss, num_missing,
                                          replace=false)
      y[i][:, j] .= y_complete[i][:, j] .+ 0

      if propmissingscale > 0
        y[i][idx_missing, j] .= NaN
      end
    end
  end

  return Dict(:Z => Z, :N => N, :L => L, :mus => mus, :W => W,
              :seed => seed, :lam => lam, :sig2 => sig2, :y => y,
              :y_complete => y_complete, :beta => beta,
              :eps_mus => eps_mus, :skew => skew)
end


text:
# developed with Julia 1.1.1
#
# functions for simulating micro grid control

function compute_stage_cost(battery::Battery, prices::Prices, t::Int64,
    control::Float64, net_energy_demand::Float64)

    control = control*battery.power*0.25
    imported_energy = (control + net_energy_demand)
    return (prices.buy[t]*max(0.,imported_energy)
        - prices.sell[t]*max(0.,-imported_energy))

end

function compute_stage_dynamics(battery::Battery, state::Float64, control::Float64)

    scale_factor = battery.power*0.25/battery.capacity
    soc = state + (battery.charge_efficiency*max(0.,control)
        - max(0.,-control)/battery.discharge_efficiency)*scale_factor
    soc = max(0., min(1., soc))
    return soc

end

function save_simulations(site::Site, simulations::Array{Simulation})
    path_to_score = joinpath(site.path_to_save_folder, site.id*".jld2")
    save(path_to_score, "simulations", simulations)
end

function group_all_simulations(sites::Array{Site})

    scores = Dict()

    for site in sites

        path_to_score = joinpath(site.path_to_save_folder, site.id*".jld2")
        simulation = load(path_to_score, "simulations")
        scores[site.id] = simulation

    end

    save(joinpath(sites[1].path_to_save_folder, "score.jld2"), scores)

    if isfile(joinpath(sites[1].path_to_save_folder, "score.jld2"))
        for site in sites
            rm(joinpath(site.path_to_save_folder, site.id*".jld2"))
        end
    end

end

### hackable functions

function initialize_site_controller(controller::AbstractController, site::Site, prices::Prices)
    """
    hackable function to update site dependent data and parameters
    and initialize the controller
    """
    return DummyController()
end

function compute_control(controller::AbstractController, information::Information)
    """
    hackable function to implement the control technique of the controller
    """
    return 0.
end


text:
using CSFML.LibCSFML
using Test

@testset "version" begin
    @test LibCSFML.CSFML_VERSION_MAJOR == 2
    @test LibCSFML.CSFML_VERSION_MINOR == 5
    @test_nowarn mode = sfVideoMode(1280, 720, 32)
end


text:
using MCMCChains
using FFTW

using Random
using Statistics
using Test

@testset "ESS per second" begin
    c1 = Chains(randn(100,5, 1), info = (start_time=time(), stop_time = time()+1))
    c2 = Chains(randn(100,5, 1), info = (start_time=time()+1, stop_time = time()+2))
    c = chainscat(c1, c2)

    wall = MCMCChains.wall_duration(c)
    compute = MCMCChains.compute_duration(c)

    @test round(wall, digits=1) ≈ round(c2.info.stop_time - c1.info.start_time, digits=1)
    @test compute ≈ (MCMCChains.compute_duration(c1) + MCMCChains.compute_duration(c2))

    s = ess_rhat(c)
    @test length(s[:,:ess_per_sec]) == 5
    @test all(map(!ismissing, s[:,:ess_per_sec]))
end

@testset "ESS and R̂ (chains)" begin
    x = rand(10_000, 40, 10)
    chain = Chains(x)

    for method in (ESSMethod(), FFTESSMethod(), BDAESSMethod())
        # analyze chain
        ess_df = ess_rhat(chain; method = method)

        # analyze array
        ess_array, rhat_array = ess_rhat(x; method = method)

        @test ess_df[:,2] == ess_array
        @test ess_df[:,3] == rhat_array
    end
end

@testset "ESS and R̂ (single sample)" begin # check that issue #137 is fixed
    val = rand(1, 5, 3)
    chain = Chains(val, ["a", "b", "c", "d", "e"])

    for method in (ESSMethod(), FFTESSMethod(), BDAESSMethod())
        # analyze chain
        ess_df = ess_rhat(chain; method = method)

        # analyze array
        ess_array, rhat_array = ess_rhat(val; method = method)

        @test ismissing(ess_df[:,2][1]) # since min(maxlag, niter - 1) = 0
        @test ismissing(ess_df[:,3][1])
    end
end


text:
using Printf

include("alias.jl")
include("utils.jl")
include("PGBF.jl")
include("CGBF.jl")
include("overlap.jl")

function test_CGBF()

    c1 = init_CGBF( 0.0,0.0,0.0 )  # specifying center
    push!(c1, 1.0, 1.0)
    println(c1)

    @printf("Test CGBF evaluate: %f\n", abs(evaluate(c1, 0.0, 0.0, 0.0)-0.71270547))

    c2 = init_CGBF( 0.0, 0.0, 0.0 )
    push!( c2, 1.0, 0.2 )
    push!( c2, 0.5, 0.2 )
    @printf("Test CGBF overlap, <c2|c2>: %f\n", abs(overlap(c2,c2)-1.0))

    @printf("Test CGBF overlap, <c1|c2>: %f\n", overlap(c1,c2))
end

test_CGBF()


text:
import Base: copy, isapprox
using Distributions: Beta, cdf
using LinearAlgebra
using Optim
using Statistics: mean

function wXageF0(ages)
    1.5 .* ages ./ (2 .+ ages)
end

function bandSparse(N, k::Integer)
    m = zeros(N, N)
    if k > 0
        for i in 1:(N-k)
            m[i, i + k] = 1.0
        end
    else
        for i in 1:(N + k)
            m[i - k, i] = 1.0
        end
    end
    m
end

function bandSparse(N, ks::AbstractArray)
    m = zeros(N, N)
    for k in ks
        if k > 0
            for i in 1:(N-k)
                m[i, i + k] = 1.0
            end
        else
            for i in 1:(N + k)
                m[i - k, i] = 1.0
            end
        end
    end
    m
end

function shiftO(N)
  bandSparse(N, 1)
end

function stageO(N, p)
  bandSparse(N, 0) * (1 - p) .+ bandSparse(N, 1) * p
end

function bdiag(A, B)
  calO = zeros(size(A)[1] + size(B)[1], size(A)[2] + size(B)[2])
  calO[1:size(A)[1], 1:size(A)[2]] .= A
  calO[(size(A)[1] + 1):end, (size(A)[2] + 1):end] .= B
  calO
end

function twoBlockO(A, B, p)
  nA = size(A)[1]
  calO = bdiag(A, B)
  calO[nA, nA + 1] = p
  calO
end

function ages_1066(unit = 365)
  x1 = collect(1:18) * 10
  x2 = maximum(x1) .+ collect(1:18) * (30 + 1 / 1.8)
  x3 = maximum(x2) .+ collect(1:18) * 365
  x4 = maximum(x3) .+ collect(1:12) * 365 * 5
  vcat(x1, x2, x3, x4) / unit
end

function calO_1066()
  A = shiftO(18)
  B = stageO(18, 1 / 3)
  C = stageO(18, 1 / 36.5)
  D = stageO(12, 1 / 36.5 / 5)
  calO = twoBlockO(A, B, 1 / 3)
  calO = twoBlockO(calO, C, 1 / 36.5)
  calO = twoBlockO(calO, D, 1 / 36.5 / 5)
  return(calO)
end

function calP_SoI(N, tau, nu, r)
  s = r .* (1 .- nu)
  v = r .* nu
  # Remain in state with probability s, go to next state probability v
  PSoI = bandSparse(N, 0) .* s .+ transpose(bandSparse(N, 1) .* v)
  PSoI = vcat(transpose(1 .- r), PSoI) # rbind
  PSoI = hcat(zeros(10), PSoI) # cbind
  # Uninfected remain uninfected.
  PSoI[1, 1] = 1
  PSoI
end

function SuSblock_Adam(N)
  SuS = bandSparse(N, 0:(N - 1))
  vcat(transpose(zeros(N + 1)), hcat(ones(N), SuS))
end

function ageXimmune(ageInDays)
    ageInDays .> (8 * 365)
end

"""
The original R:
xiF_0 <- function(x, p.xi = 2, N = 9, C = 10, D = 0) {
  diff(pbeta(c(0:N) / N, (1 - exp(-p.xi * x)) * C, exp(-p.xi * x) * C + D))
}
The beta distribution has the same density in both languages. I checked.
"""
function xiF_inner(x; p_xi = 2, N = 9, C = 10, D = 0)
    if (x > 0)
        α = (1 - exp(-p_xi * x)) * C
        β = exp(-p_xi * x) * C + D
        dist = Beta(α, β)
        diff(cdf.(dist, collect(0:N) / N))
    else
        bins = zeros(N)
        bins[1] = 1
        bins
    end
end

# I can't tell whether V is an array or a single float when it gets used.
# If it's an array, use pdf.() above, and it will work out.
function xiF_h(pit, V ; p_xi1 = 5, p_xi2 = 1.5, N = 9, eps = 0.001)
  if V > 0
      v = 1 - exp(-(pit[2] + p_xi2 * pit[3]) / V)
  else
      v = 1
  end
  xi1 = xiF_inner(v, p_xi = p_xi1)
  xi = xi1 .+ eps
  xi / sum(xi)
end

function make_parameters_adam()
    tau = 10
    N = 9
    nuVals = [1/2, 1/6, 1/7, 1/8, 1/9, 1/9, 1/9, 1/12, 0]
    rVals = vcat([0, 1/1000], ones(6) / 250, [1/400])
    calDtrue_SoI =vcat([0], ones(9))
    calDlm_SoI = vcat([0, 0.99, .95], 1 .- collect(3:7).^3 / 8.6^3, [.05, .001])
    calDimm_SoI = vcat([0, 0.95, .9, .6, .3], zeros(5))
    calK_SoI = vcat([0, 0], collect(8:-1:1) / 10)
    calBsevere_SoI = vcat([0, .3], zeros(8))
    calBfever_SoI = vcat([0, 1, .7], exp.(-collect(1:6)), [0])
    rho = 0.2
    etaXstage = vcat([0, 0.99, .95], 1 .- collect(3:7).^3 / 8.6^3, [.05, .001]).^2
    rhoXstage = [0, 1, 1, 1, .8, .8, .6, .6, .4, .4]
    ageForImmune = 8
    delta = 1 / 730
    mu1 = 0.001
    mu2 = 0.01
    mu3 = 0.001
    wF = wXageF0
    agesF = ages_1066
    calOF = calO_1066
    # Probability of taking drugs, given _any_ state. It's a background rate.
    Dd = 1 - exp(-delta / tau)
    r = exp.(-rVals * tau)

    PSoI = calP_SoI(N, tau, nuVals, r)
    SuS = SuSblock_Adam(N)
    ages = agesF()
    aXi = (ages .> 8)
    demog = vcat(diff(ages), [15])
    demog = demog / sum(demog)
    demog29 = (ages .>= 2) .& (ages .< 10)
    demog29 = demog29 / sum(demog29)
    U5 = copy(demog)
    O5 = copy(demog)
    O5[ages .< 5] .= 0
    U5[ages .>= 5] .= 0
    calO = calOF()
    w = wF(ages)
    calBfever = vcat(calBfever_SoI, calBfever_SoI, [0, 0])
    trt = rho .* rhoXstage .* etaXstage
    calBtrtdmal = vcat(trt .* calBfever_SoI, trt .* calBfever_SoI, [0, 0])
    calBuntrmal = vcat((1 .- trt) .* calBfever_SoI, (1 .- trt) .* calBfever_SoI, [0, 0])
    calBsevere = vcat(calBsevere_SoI, calBsevere_SoI, [0, 0])
    calBtrtdsev = vcat(trt .* calBsevere_SoI, trt .* calBsevere_SoI, [0, 0])
    calBuntrsev = vcat((1 .- trt) .* calBsevere_SoI, (1 .- trt) .* calBsevere_SoI, [0, 0])
    calK = vcat(calK_SoI, calK_SoI, [0, 0])
    calK[12] = calK[3]
    Dict{Symbol,Any}(
      :tau => tau,
      :N => N,
      :calDtrue => vcat(calDtrue_SoI, ones(10), [0, 0]),
      :calDlm => vcat(calDlm_SoI, calDlm_SoI, [0, 0]),
      :calDimm => vcat(calDimm_SoI, calDimm_SoI, [0, 0]),
      :calBfever => calBfever,
      :calBtrtdmal => calBtrtdmal,
      :calBuntrmal => calBuntrmal,
      :calBsevere => calBsevere,
      :calBtrtdsev => calBtrtdsev,
      :calBuntrsev => calBuntrsev,
      :calK => calK,
      :nu => nuVals,
      :r => r,
      :Dd => Dd,
      :rho => rho,
      :etaXstage => etaXstage,
      :rhoXstage => rhoXstage,
      :PSoI => PSoI,
      :SuS => SuS,
      :wF => wF,
      :w => w,
      :ages => ages,
      :ageXimmune => ageXimmune,
      :aXi => aXi,
      :demog => demog,
      :demog29 => demog29,
      :U5 => U5,
      :O5 => O5,
      :calO => calO,
      :mu1 => mu1,
      :mu2 => mu2,
      :mu3 => mu3,
      :xiF => xiF_h
    )
end

### Now we make the Adam model from those parameters

"""
A single group. The code may use cohort/population in the opposite sense.
I mean here that this is 22 compartments for one group of people.
"""
mutable struct CohortState
    X::Array{Float64,2}
    pit::Array{Float64,1}
    V::Float64
    ageInDays::Float64
    alpha1::Float64
    alpha2::Float64
end

function emptyX_Adam()
    X = zeros(22, 1)
    X[1] = 1
    CohortState(X, zeros(3), 0, 0, 0, 0)
end

function isapprox(a::CohortState, b::CohortState)
    (a.alpha1 ≈ b.alpha1 && a.alpha2 ≈ b.alpha2 && a.V ≈ b.V) || return(false)
    a.X ≈ b.X || return(false)
    a.pit ≈ b.pit || return(false)
    a.ageInDays ≈ b.ageInDays || return(false)
    true
end

"""
This tells me where they differ.
"""
function differences(a::CohortState, b::CohortState)
    msg = String[]
    a.alpha1 ≈ b.alpha1 || push!(msg, "alpha1")
    a.alpha2 ≈ b.alpha2 || push!(msg, "alpha2")
    a.V ≈ b.V || push!(msg, "V")
    a.X ≈ b.X || push!(msg, "X")
    a.pit ≈ b.pit || push!(msg, "pit")
    a.ageInDays ≈ b.ageInDays || push!(msg, "ageInDays")
    msg
end

"""
Multiple groups. There are 22 compartments for each of the
groups. The code may use cohort to mean this and population to mean
the CohortState, which I get, because this tracks multiple cohorts.
"""
mutable struct PopulationState
    X::Array{Float64,2}
    pit::Array{Float64,2}
    V::Float64
    ageInDays::Array{Float64,2}
    alpha1::Float64
    alpha2::Float64
end

function emptyXX_Adam(; L = 66)
    X = zeros(22, L)
    X[1, :] .= 1
    pit = zeros(3, L)
    ageInDays = zeros(1, L)
    PopulationState(X, pit, 0, ageInDays, 0, 0)
end

function isapprox(a::PopulationState, b::PopulationState)
    (a.alpha1 ≈ b.alpha1 && a.alpha2 ≈ b.alpha2 && a.V ≈ b.V) || return(false)
    a.X ≈ b.X || return(false)
    a.pit ≈ b.pit || return(false)
    a.ageInDays ≈ b.ageInDays || return(false)
    true
end

"""
This tells me where they differ.
"""
function differences(a::PopulationState, b::PopulationState)
    msg = String[]
    a.alpha1 ≈ b.alpha1 || push!(msg, "alpha1")
    a.alpha2 ≈ b.alpha2 || push!(msg, "alpha2")
    a.V ≈ b.V || push!(msg, "V")
    a.X ≈ b.X || push!(msg, "X")
    a.pit ≈ b.pit || push!(msg, "pit")
    a.ageInDays ≈ b.ageInDays || push!(msg, "ageInDays")
    msg
end

"""
Acts on members of the cohort state.
"""
function calP_Adam(alpha, pit, V, params)
    SuS = params[:SuS]
    PSoI = params[:PSoI]
    alpha = min(alpha, 1)
    xi = params[:xiF](pit, V)
    E1 = SuS .* vcat(0, xi)  # multiply each column
    e1 = sum(E1; dims=(1,))
    e1[isapprox.(e1, 0)] .= 1
    hatxi = E1 * ((1 ./ e1) .* I(10))
    Deta = min.(params[:rho] * params[:etaXstage] .* params[:rhoXstage], 1)
    vecD = Deta .+ (1 .- Deta) .* params[:Dd]
    notTreat = Diagonal(1 .- vecD)
    treat = vcat(vecD, vecD)
    B1 = PSoI * notTreat
    @assert size(B1) == (10, 10)
    B2 = hatxi * notTreat
    @assert size(B2) == (10, 10)
    chi = vcat([1], cumsum(xi))
    calP = vcat(
        hcat((1 .- alpha .* chi) .* B1, (1 -alpha) .* B2),
        hcat(alpha * chi .* B1, alpha * B2),
        transpose(treat),
        transpose(0 * treat)
    )
    calP = hcat(calP, zeros(size(calP)[1], 2))
    calP[22, 21] = 1
    calP[1, 22] = 1
    @assert size(calP) == (22, 22)
    calP
end

function PxX_Adam!(vecX::CohortState, alpha, params)
    vecX.ageInDays += params[:tau]
    w = params[:wF](vecX.ageInDays / 365)
    Xt = calP_Adam(w * alpha, vecX.pit, vecX.V, params) * vecX.X

    vecX.V = alpha + vecX.V * (1 - params[:mu1])
    vecX.pit[1] = vecX.pit[1] * (1 - params[:mu1]) + vecX.alpha2
    vecX.pit[2] = (vecX.pit[2] * (1 - params[:mu2]) +
            dot(params[:calDimm],  vecX.X))
    vecX.pit[3] = (vecX.pit[3] * (1 - params[:mu3]) +
            vecX.alpha2 * params[:ageXimmune](vecX.ageInDays))
    vecX.alpha2 = vecX.alpha1
    vecX.alpha1 = w * alpha
    vecX.X = Xt
end

"""
Simulates a cohort from birth over mx time steps and returns it
as a PopulationState.
"""
function cohortXX_Adam(alpha, params; mx = 2920)
    # Grow a cohort so you can assign its stages to a population.
    vecX = emptyX_Adam()
    XX = emptyXX_Adam(L = mx)
    vecX.V = alpha / params[:mu1]
    for i in 1:mx
        # An array column is size (22,) but an array with one column is (22,1)
        XX.X[:, i] .= vecX.X[:, 1]
        XX.pit[:, i] = vecX.pit[:, 1]
        XX.ageInDays[1, i] = vecX.ageInDays
        PxX_Adam!(vecX, alpha, params)
    end
    XX.alpha2 = vecX.alpha2
    XX.alpha1 = vecX.alpha1
    XX.V = vecX.V
    XX
end

function cohort2ages_Adam(cXX::PopulationState, params)
    ages = params[:ages]
    # ageInDays is a row vector. Easier to work with column here.
    ageInYears = vec(cXX.ageInDays / 365)
    XX = emptyXX_Adam()
    XX.V = cXX.V
    XX.ageInDays = zeros(1, length(ages))
    no_ages = []
    lower = 0.0
    for i in 1:length(ages)
        upper = ages[i]
        ix = (ageInYears .>= lower) .& (ageInYears .< upper)
        age_cnt = sum(ix)
        if age_cnt > 0
            XX.X[:, i] = sum(cXX.X[:, ix], dims = 2) / age_cnt
            XX.pit[:, i] = sum(cXX.pit[:, ix], dims = 2) / age_cnt
            XX.ageInDays[1, i] = mean(cXX.ageInDays[ix])
        else
            push!(no_ages, ages[i])
        end
        lower = upper
    end
    if length(no_ages) > 0
        println("no ages for $no_ages")
        @assert length(no_ages) == 0
    end
    XX
end

function test_cohort2ages_Adam()
    params = make_parameters_adam()
    XX = cohortXX_Adam(0.3, params)
    cohort2ages_Adam(XX, params)
end

function PtimesX_Adam!(XX::PopulationState, alpha, params)
    ageInDays = XX.ageInDays .+ params[:tau] * params[:calO]
    Xt = copy(XX.X)
    for a in 1:66
        Xt[:, a] .= calP_Adam(
            alpha * params[:w][a],
            XX.pit[:, a],
            XX.V,
            params
            ) * XX.X[:, a]
    end
    Xt *= params[:calO]
    Xt[1, 1] = 1

    XX.V = alpha + XX.V * (1 - params[:mu1])
    XX.pit[1, :] .= XX.pit[1, :] * (1 - params[:mu1]) .+ alpha
    XX.pit[2, :] .= (XX.pit[2, :] * (1 - params[:mu2]) .+
            transpose(XX.X) * params[:calDimm])
    XX.pit[3, :] .= XX.pit[3, :] * (1 - params[:mu3]) .+ XX.alpha2 * params[:aXi]

    XX.pit *= params[:calO]
    XX.ageInDays = XX.ageInDays * params[:calO]
    XX.alpha2 = XX.alpha1
    XX.alpha1 = alpha
    XX.X = Xt
end

function XX2pr29_Adam(XX, params)
    transpose(params[:calDlm]) * XX.X * params[:demog29]
end

function ar2stableWave_Adam(alpha, params; tol = 0.01)
    cXX = cohortXX_Adam(alpha, params)
    XX = cohort2ages_Adam(cXX, params)
    for i in 1:100
        PtimesX_Adam!(XX, alpha, params)
    end
    x = XX2pr29_Adam(XX, params)
    df = 1
    while df > tol
        original = copy(XX.X)
        PtimesX_Adam!(XX, alpha, params)
        x = XX2pr29_Adam(XX, params)
        df = sum(abs.(original .- XX.X))
    end
    XX
end

function ar2pr_Adam(alpha, params)
    XX = ar2stableWave_Adam(alpha, params)
    XX2pr29_Adam(XX, params)
end

"""
    pr2arSS_Adam(x, params)

Given a PfPR value and model parameters, return the attack rate
that would produce this PfPR for a population. This is an optimization
step.
"""
function pr2arSS_Adam(x, params)
    objective = x -> (x - ar2pr_Adam(alpha, params))^2
    lower = 0.0
    upper = 1.0
    res = optimize(objective, lower, upper)
    Optim.minimizer(res), Optim.iterations(res), Optim.converged(res)
end

function pr2ar_Adam(xx, params; alpha0 = nothing, XX0 = nothing, tol = 0.01)
    if alpha0 === nothing
        alpha = pr2arSS_Adam(xx[1], params)
    else
        alpha = alpha0
    end
end


text:
function compute(n::Int)::Int
    limit_powers, limit_numbers = trunc(Int, log2(n)), trunc(Int, √n)
    power_counts = repeat([0], limit_powers)
    for i ∈ 2:limit_numbers
        power = trunc(Int, log(i, n))
        power_counts[power] += 1
        j = i * i
        while j ≤ limit_numbers
            power_counts[trunc(Int, log(j, n))] -= 1
            j *= i
        end
    end
    sieve = repeat([false], limit_powers * n)
    duplicate_counts = repeat([0], limit_powers)
    count = 0
    for i ∈ 1:limit_powers
        for j ∈ 2:n
            if sieve[i * j]
                count += 1
            else
                sieve[i * j] = true
            end
        end
        duplicate_counts[i] = count
    end
    return (n - 1) ^ 2 - sum(power_counts .* duplicate_counts)
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# These tests are not currently run from runtests.jl

# Tests may currently fail 'randomly' because of the limitations
# of running everything in a single thread, using libuv.
# For example, trying to connect while the websocket is
# busy in a blocking operation can cause timeouts somewhere.
#
# The test set outside of this is now complete, and /benchmark
# is a more effective context for testing the efficiency.
#
# TODO
# 'functions_open_browsers.jl' is an older version
#     of the same file name from /benchmark. Using that requires
#     some restructuring of this file.
# 'functions_log_test.jl' is an older version of
#     /logutils/...  Use the new methods.

# Launches available browsers on browsertest.html.
# Websockets are intiated by the browser / web page.
# Handler_functions_websockets respond and store allocations & etc.
# After a successfull exchange, the browsers navigate to a second html page
# for binary performance tests. In the end, all websockets are closed
# and a summary is output.

cd(Pkg.dir("WebSockets", "test", "optional"))
using Compat
using HttpServer
using WebSockets
using Base.Test
using JSON
const WEBSOCKETS = Dict{String, WebSockets.WebSocket}()
const RECEIVED_WS_MSGS_TIME = Dict{String, Vector{Float64}}()
const RECEIVED_WS_MSGS_ALLOCATED = Dict{String, Vector{Int64}}()
const RECEIVED_WS_MSGS_LENGTH = Dict{String, Vector{Int64}}()
const WEBSOCKETS_SUBPROTOCOL = Dict{String, WebSockets.WebSocket}()
const WEBSOCKETS_BINARY = Dict{String, WebSockets.WebSocket}()
# Note that only text messages are stored. Binary messages are discarded.
const RECEIVED_WS_MSGS = Dict{String, Vector{String}}()

# Logs from travis indicate it takes 7 s to fire up Safari and
# have the first websocket open. This can easily increase if more
# browsers become available.
const FIRSTWAIT = Base.Dates.Second(2)
const MAXWAIT = Base.Dates.Second(60*15)
global n_responders = 0
include("functions_server.jl")
closeall()
server = start_ws_server_async()
include("functions_open_browsers.jl")
info("This OS is $(string(Sys.KERNEL))\n")
n_browsers = 0
#n_browsers += open_testpage("phantom")
n_browsers += open_all_browsers()#

# Control flow passes to async handler functions

if n_browsers > 0
    info("Sleeping main thread for an initial minimum of $FIRSTWAIT\n")
    t0 = now()
    while now()-t0 < FIRSTWAIT
        sleep(1)
    end
    contwait = true
    while contwait && now()-t0 < MAXWAIT
        # All websockets close briefly when navigating from browsertest.html to browsertest2.html
        # So we require two consecutive checks for websockets before exiting, with 2 seconds between.
        contwait = count_open_websockets() > 0
        if !contwait
            sleep(2)
            passedtime = now() - t0
            exitlatestin = MAXWAIT + t0- now()
            info("$(div(passedtime.value, 1000)) s passed. Max remaining test time $(div(exitlatestin.value, 1000)) s\n")
        end
        sleep(.5)
        contwait = contwait && count_open_websockets() > 0
    end
end
n_opensockets = count_open_websockets()

closeall()
info(Dates.format(now(), "HH:MM:SS"), "\tClosed web sockets and servers.")
# sum up
allocs = Vector{Int64}()
times = Vector{Float64}()
lengths = Vector{Int64}()
n_msgs = 0
n_text = 0
n_ws = 0
for ke in keys(WEBSOCKETS)
    n_ws += 1
    for msgno = 1:length(RECEIVED_WS_MSGS_LENGTH[ke])
        if RECEIVED_WS_MSGS_LENGTH[ke][msgno] > 0
            n_msgs += 1
            push!(allocs, RECEIVED_WS_MSGS_ALLOCATED[ke][msgno])
            push!(times, RECEIVED_WS_MSGS_TIME[ke][msgno])
            push!(lengths, RECEIVED_WS_MSGS_LENGTH[ke][msgno])
            if haskey(RECEIVED_WS_MSGS, ke)
                n_text += 1
            end
        end
    end
end
n_binary = n_msgs - n_text

# print summary
info("Spawned $n_browsers browsers. $n_browsers made requests. Opened $n_ws sockets, $n_opensockets did not close as intended.")
info("Received $n_msgs messages, $n_text text and $n_binary binary, $(round(sum(lengths)/ 1000 / 1000,3)) Mb, sent a similar amount.")

if n_msgs > 0
    maxlength = maximum( lengths )
    minlength = minimum( va -> va > 0.0 ? va:typemax(va), lengths )
    avglength = sum(lengths) / n_msgs

    maxaloc = maximum(  va -> va < Inf ? va:0.0, allocs ./ lengths )
    minaloc = minimum( va -> va > 0.0 ? va:typemax(va), allocs ./ lengths )
    avgaloc = sum(allocs) / sum(lengths)

    maxtime = maximum(va -> va < Inf ? va:0.0, times)
    mintime = minimum(va -> va > 0.0 ? va:typemax(va), times)
    avgtime = sum(times) / n_msgs

    maxspeed = maximum(va -> va < Inf ? va:0.0, lengths ./ times)
    minspeed = minimum(va -> va > 0.0 ? va:typemax(va), lengths ./ times)
    avgspeed = sum(lengths) / sum(times)

    info("Length of messages received\n",
        "\t\tAverage length:\t\t", round(avglength / 1000, 3), " kB\n",
        "\t\t\Minimum :\t\t", minlength , " b\n",
        "\t\t\Maximum :\t\t", round(maxlength / 1000 / 1000 , 3), " Mb\n")

    info("Time spent reading and waiting for received messages\n",
        "\t\tAverage time:\t\t", round(avgtime, 4), " s\n",
        "\t\t\Minimum :\t\t", mintime , " s\n",
        "\t\t\Maximum :\t\t", maxtime, " s\n")

    info("Reading speed (strongly affected by the active browsers)\n",
        "\t\tAverage speed:\t", round(avgspeed / 1000 / 1000, 3), " Mb/s\n",
        "\t\t\Minimum :\t", round(minspeed / 1000 / 1000 , 10), " Mb/s\n",
        "\t\t\Maximum :\t", round(maxspeed / 1000 / 1000, 3), " Mb/s\n")

    info("Allocations for reading, bytes allocated per byte received\n",
        "\t\tAverage allocation:\t", round(avgaloc, 3), "\n",
        "\t\t\Minimum :\t\t", round(minaloc, 3), "\n",
        "\t\t\Maximum :\t\t", round(maxaloc, 3), "\n",)

    info("Text messages received per websocket:")
    for m in RECEIVED_WS_MSGS
        display(m)
    end

else
    info("Failure on speed / allocation test.")
end

@test n_opensockets == 0
@test n_msgs == n_responders * 37


text:
push!(LOAD_PATH,"../src/")
using Documenter
using Plots
using ParticleInCell

ENV["GKSwstype"] = "100"

println("using Pkg; Pkg.develop(PackageSpec(path=pwd()))")

makedocs(modules=[ParticleInCell],
         sitename = "ParticleInCell.jl",
         authors="Pierre Navaro",
         format=Documenter.HTML(;
         prettyurls=get(ENV, "CI", "false") == "true",
         canonical="https://juliavlasov.github.io/ParticleInCell.jl",
         assets=String[],
         ),
         doctest = false,
         pages = ["Home"           => "index.md",
                  "Maxwell solver" => "maxwell.md",
                  "Landau damping" => "landau_damping.md",
                  "Two-stream instability" => "tsi.md",
                  "Vlasov-Maxwell 2D" => "vlasov-maxwell.md",
                  "Vlasov-Poisson 2D" => "vlasov-poisson.md",
                  "Types"     => "types.md",
                  "Functions" => "functions.md",
                  "Contents"  => "contents.md"])

deploydocs(;
    branch = "gh-pages",
    devbranch = "master",
    repo   = "github.com/JuliaVlasov/ParticleInCell.jl"
)


text:
using BinaryBuilder

include("../../../fancy_toys.jl")

name = "CUDA_full"
version = v"10.2.89"

sources_linux = [
    FileSource("http://developer.download.nvidia.com/compute/cuda/10.2/Prod/local_installers/cuda_10.2.89_440.33.01_linux.run",
               "560d07fdcf4a46717f2242948cd4f92c5f9b6fc7eae10dd996614da913d5ca11", "installer.run")
]
sources_macos = [
    FileSource("http://developer.download.nvidia.com/compute/cuda/10.2/Prod/local_installers/cuda_10.2.89_mac.dmg",
               "51193fff427aad0a3a15223b1a202a6c6f0964fcc6fb0e6c77ca7cd5b6944d20", "installer.dmg")
]
sources_windows = [
    FileSource("http://developer.download.nvidia.com/compute/cuda/10.2/Prod/local_installers/cuda_10.2.89_441.22_win10.exe",
               "b538271c4d9ffce1a8520bf992d9bd23854f0f29cee67f48c6139e4cf301e253", "installer.exe")
]

sources_linux_aarch64 = [
    # Expected install order for cuda-toolkit-10-2, i.e. output from `apt-get -s install --no-install-recommends cuda-toolkit-10-2` (using source: "deb https://repo.download.nvidia.com/jetson/common r32.5 main")
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-license-10-2_10.2.89-1_arm64.deb", "8a862acbff5b33904bfe7ec3e92a553a8312da1db9f651b6cfe14db137a139ce"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvdisasm-10-2_10.2.89-1_arm64.deb", "9523033692cca5a2b29cd69942d69036deedacb8eb3273395662f6024b2c27f9"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cuobjdump-10-2_10.2.89-1_arm64.deb", "719b32f039cd8ed6123e1f9e3fa9badcf4c6ba5ac4a0b24a97d2db88e0764e1e"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-gdb-10-2_10.2.89-1_arm64.deb", "d3fef316b7d5b215cf11ff0190c69d9e8a2652b9dd37f454c7350110954e3496"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvprof-10-2_10.2.89-1_arm64.deb", "54d96b9cdba5a53da92b2cdaada27ec7a886b3a6a29e7d33b5fdf429ad788681"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-memcheck-10-2_10.2.89-1_arm64.deb", "f8f498108e1e4e3fdc40a1cb80cec44bb3da987b2160a3128d3bcf85a8533bbf"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cudart-10-2_10.2.89-1_arm64.deb", "2a718596cf1162bf0076d4ec6db52a5f7c3617b7ab7cd243887376e841a99915"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-driver-dev-10-2_10.2.89-1_arm64.deb", "c1c55ba59d8a28a7d56800504c65683a4d392893f35d08ec7bddf6b45efda468"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cudart-dev-10-2_10.2.89-1_arm64.deb", "5aa2bf1e8e9d467dacbff778b0d2d4a7bd31077a443b7a49711cc798562ea37d"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-misc-headers-10-2_10.2.89-1_arm64.deb", "e92834f576241295c74393b478cb7121d9e20cb3454aed26c464c3523eaeadde"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvcc-10-2_10.2.89-1_arm64.deb", "1a0ea57d4c1b1d9394d7e4f6ab94baa2aa49883f4ba2d59a60b750bb88d0fdeb"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cupti-10-2_10.2.89-1_arm64.deb", "319771f42db1d9a4a273bc5ff753148247ece9cf7544d2008d57d9061e6964f9"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cupti-dev-10-2_10.2.89-1_arm64.deb", "3bd27507b8eef4ae9d5faf8671e8843ef7bd4cede82dc76a8e6519390e22a5dc"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvtx-10-2_10.2.89-1_arm64.deb", "b4f9dedf0c21dbc75daadbabfde6aa17a73b6bdad0d5a18fae3da6e9717bfc99"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-command-line-tools-10-2_10.2.89-1_arm64.deb", "7425afa9751b073b709f969149babb65411e2ab96bbe1744af0d89a517a6f1a2"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvprune-10-2_10.2.89-1_arm64.deb", "1e04820f53fb96b737c6b5d92ae2c1c54414e84a87ae9a55a00fc78c05e4e33f"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-compiler-10-2_10.2.89-1_arm64.deb", "cf15fd18669b88dfd9e49be6e8189e03936ce993be47dfc88b99d1b3d86d3a6d"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cufft-10-2_10.2.89-1_arm64.deb", "7d06a0bdd16b16c2384de79e488cde5bb580b3bc9ee46bd1e23d17ce6d260e01"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cufft-dev-10-2_10.2.89-1_arm64.deb", "239060c552e98511b8a752127c57bd86cd1261af5062b99f85d94a12ef61d876"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-curand-10-2_10.2.89-1_arm64.deb", "1cf024814383a2ec21e562e1e89809050dbae36bceed6ee276affe03185ca266"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-curand-dev-10-2_10.2.89-1_arm64.deb", "f6594cd55b89a4e45ca3a382cef99c6bd3d2efb7ab222a82f261a79cc9e066ce"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cusolver-10-2_10.2.89-1_arm64.deb", "ede4930b8a7d8098590b3691785cf225fdbb08e3cb7853ca8f0ffc7313e89b3f"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cusolver-dev-10-2_10.2.89-1_arm64.deb", "e3d786efc6d92e24c2527a2a76685de4cfc92e609a3b6ab1258cc820a5c867e4"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cusparse-10-2_10.2.89-1_arm64.deb", "459f9e41ac458eef3c94ba0a5d367182af83c070454d9e38b3565252c039a827"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-cusparse-dev-10-2_10.2.89-1_arm64.deb", "3a9ccd17e6916fda5cede14c206a8199fdba27f1118947a1f31cbf33132241df"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvrtc-10-2_10.2.89-1_arm64.deb", "fe9cc7cbdab29035371e5f77d190aecd97bbf31ca66b3f62612bf62580417b16"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvrtc-dev-10-2_10.2.89-1_arm64.deb", "5c3bd9bc84170ec0c3465444409e7912046e2523e6a18d331b53901c0e57d229"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cublas/libcublas10_10.2.2.89-1_arm64.deb", "d0299b139a163136432dfb2c028769944b6c5636ad9238614860c196a1c91aea"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cublas/libcublas-dev_10.2.2.89-1_arm64.deb", "5fa7e3e8fe266fdea7e91778610b7e8d3d85d8950875a4915ce3626c9e564365"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-npp-10-2_10.2.89-1_arm64.deb", "a01a204e6afefb0072424817c7133c5c18a9f3996fbf48f0b31297cb7c07e1ac"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-npp-dev-10-2_10.2.89-1_arm64.deb", "4d572a2472564f5008c3cfe9b24dbef765a1531a73c43436c66a146b4d16cace"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvgraph-10-2_10.2.89-1_arm64.deb", "f800974ac6cb3fb6595a08e613fc0f376d7b91f43eed1b5306b0496f9588e441"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvgraph-dev-10-2_10.2.89-1_arm64.deb", "cc2d9897c54f27a20f90bea5df391875b3f8163ec69745fb7546c1ce57e1d718"),
    # FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda-samples/cuda-samples-10-2_10.2.89-1_arm64.deb", "121e273d8586bde904ceeab72a603a86d781f3bac6d3a21732703ca2ca9ec528"),
    # FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-documentation-10-2_10.2.89-1_arm64.deb", "746f675b08f29cdf7d9fe39edb325a960ae29106da805369b8dc509e72ec2329"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-libraries-10-2_10.2.89-1_arm64.deb", "71309637bbc86b3d7f014c5ba5da709b88f69dfc0710aa8caa661b3aa40a51b4"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-libraries-dev-10-2_10.2.89-1_arm64.deb", "9d3e1e4cf097f6e04c7f881b4af1a54c52552fb2828f73f7ae8fa5b9efe32d29"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-nvml-dev-10-2_10.2.89-1_arm64.deb", "c8743e69c84a432c5e6dea6edfcacf1bb6b09b028bee61c8aece7a41d0447265"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-tools-10-2_10.2.89-1_arm64.deb", "eeaae6a103cc5d950bf4e774a277d90286492df0e0815e161f7c12291c9aa5ba"),
    FileSource("https://repo.download.nvidia.com/jetson/common/pool/main/c/cuda/cuda-toolkit-10-2_10.2.89-1_arm64.deb", "47a1c8c8bde2c763396a68c38cc901217e1606bb536ec5880d4c072c4eee9073"),
]
sources_linux_aarch64_manifest = join([src.filename for src in sources_linux_aarch64], "\n")

script = """
if [[ \${target} == aarch64-linux-gnu ]]; then
    echo \"$(sources_linux_aarch64_manifest)\" >> \${WORKSPACE}/srcdir/aarch64-linux-gnu_install.manifest
fi
""" * raw"""
cd ${WORKSPACE}/srcdir

# use a temporary directory to avoid running out of tmpfs in srcdir on Travis
temp=${WORKSPACE}/tmpdir
mkdir ${temp}

apk add p7zip

mkdir ${prefix}/cuda
if [[ ${target} == x86_64-linux-gnu ]]; then
    sh installer.run --tmpdir="${temp}" --target "${temp}" --noexec
    cd ${temp}/builds/cuda-toolkit
    find .

    # clean-up
    rm -r libnsight libnvvp nsightee_plugins nsight-compute-2019.5.0 nsight-systems-2019.5.2 doc

    mv * ${prefix}/cuda
elif [[ ${target} == x86_64-apple-darwin* ]]; then
    7z x installer.dmg 5.hfs -o${temp}
    cd ${temp}
    7z x 5.hfs
    tar -zxf CUDAMacOSXInstaller/CUDAMacOSXInstaller.app/Contents/Resources/payload/cuda_mac_installer_tk.tar.gz
    cd Developer/NVIDIA/CUDA-*/
    find .

    # clean-up
    rm -r libnsight libnvvp nsightee_plugins nsight-compute-2019.5.0 NsightSystems-2019.5.2.16 doc

    mv * ${prefix}/cuda
elif [[ ${target} == x86_64-w64-mingw32 ]]; then
    7z x installer.exe -o${temp}
    cd ${temp}
    find .

    mv EULA.txt ${prefix}/cuda

    for project in cuobjdump memcheck nvcc cupti nvdisasm curand cusparse npp cufft \
                   cublas cudart cusolver nvrtc nvgraph nvprof nvprune; do
        [[ -d ${project} ]] || { echo "${project} does not exist!"; exit 1; }
        cp -a ${project}/* ${prefix}/cuda
    done

    # NVIDIA Tools Extension Library
    7z x "CUDAVisualStudioIntegration/NVIDIA NVTX Installer.x86_64".*.msi -o${temp}/nvtx_installer
    find nvtx_installer
    for file in nvtx_installer/*.*_*; do
        mv $file $(echo $file | sed 's/\.\(\w*\)_.*/.\1/')
    done
    mv nvtx_installer/*.dll ${prefix}/cuda/bin
    mv nvtx_installer/*64_*.lib ${prefix}/cuda/lib/x64
    mv nvtx_installer/*32_*.lib ${prefix}/cuda/lib/Win32
    mv nvtx_installer/*.h ${prefix}/cuda/include

    # fixup
    chmod +x ${prefix}/cuda/bin/*.{exe,dll}

    # clean-up
    rm ${prefix}/cuda/*.nvi
elif [[ ${target} == aarch64-linux-gnu ]]; then
    apk add dpkg
    cat ${target}_install.manifest | xargs -t -Ideb_file dpkg-deb -x deb_file $temp
    find $temp
    mv $temp/usr/local/cuda-10.2/* ${prefix}/cuda

    rsync -aOv --remove-source-files $temp/usr/include/ ${prefix}/cuda/targets/aarch64-linux/include/
    rsync -aOv --remove-source-files $temp/usr/lib/$target/ ${prefix}/cuda/targets/aarch64-linux/lib/

    # Move CUPTI into same paths as CUPTI on x86_64-linux-gnu
    mkdir $prefix/cuda/extras/CUPTI/{include,lib64}
    CUDA_TARGETS_REGEX='./usr/local/cuda-10.2/(targets/[^ ]+).*$'
    CUPTI_PKGS="cuda-cupti-10-2 cuda-cupti-dev-10-2"
    CUPTI_PATHS=""
    for CUPTI_PKG in $CUPTI_PKGS; do
        CUPTI_PATHS+=" "`dpkg-deb -c $CUPTI_PKG* | grep -E "$CUDA_TARGETS_REGEX" | sed -E -e "s#.*$CUDA_TARGETS_REGEX#\1#"`
    done
    for CUPTI_PATH in $CUPTI_PATHS; do
        CUPTI_PATH=$prefix/cuda/$CUPTI_PATH
        [ -f $CUPTI_PATH ] || [ -L $CUPTI_PATH ] && [[ $CUPTI_PATH == */include* ]] && mv -nv $CUPTI_PATH $prefix/cuda/extras/CUPTI/include
        [ -f $CUPTI_PATH ] || [ -L $CUPTI_PATH ] && [[ $CUPTI_PATH == */lib* ]] && mv -nv $CUPTI_PATH $prefix/cuda/extras/CUPTI/lib64
    done

    mv ${prefix}/cuda/doc/EULA.txt ${prefix}/cuda
    rmdir ${prefix}/cuda/doc
fi

cd ${prefix}/cuda
install_license EULA.txt
"""

products = Product[
    # this JLL isn't meant for use by Julia packages, but only as build dependency
]

dependencies = []

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

if should_build_platform("x86_64-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux, script,
                   [Platform("x86_64", "linux")], products, dependencies;
                   skip_audit=true)
end

if should_build_platform("x86_64-apple-darwin14")
    build_tarballs(non_reg_ARGS, name, version, sources_macos, script,
                   [Platform("x86_64", "macos")], products, dependencies;
                   skip_audit=true)
end

if should_build_platform("x86_64-w64-mingw32")
    build_tarballs(ARGS, name, version, sources_windows, script,
                   [Platform("x86_64", "windows")], products, dependencies;
                   skip_audit=true)
end

if should_build_platform("aarch64-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux_aarch64, script,
                   [Platform("aarch64", "linux")], products, dependencies;
                   skip_audit=true)
end

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
struct Simplex{T<:Number, U<:Complex}
  vertices::AbstractVector{Vertex{T,U}}
  function Simplex{T,U}(vertices::AbstractVector{Vertex{T,U}}
      ) where {T<:Number, U<:Complex}
    sort!(vertices, by=v->angle(value(v)))
    return new{T,U}(vertices)
  end
end

import Base: length, iterate, push!, iterate, getindex
import Base: eachindex, sort!, setindex!, hash, isequal
Base.length(s::Simplex) = length(s.vertices)
Base.iterate(s::Simplex) = iterate(s.vertices)
Base.iterate(s::Simplex, counter) = iterate(s.vertices, counter)
Base.getindex(s::Simplex, index) = s.vertices[index]
Base.sort!(s::Simplex; kwargs...) = sort!(s.vertices; kwargs...)
function Base.setindex!(s::Simplex, entry, index)
  s.vertices[index] = entry
  sortbyangle!(s)
  return nothing
end
Base.hash(s::Simplex) = hash(s, hash(:Simplex))
Base.hash(s::Simplex, h::UInt64) = hash(hash.(s), h)
Base.isequal(a::Simplex, b::Simplex) = all(isequal.(a, b))

sortbyangle!(s::Simplex) = sort!(s, by=v->angle(value(v)))
issortedbyangle(s::Simplex) = issorted(s, by=v->angle(value(v)))

dimensionality(s::Simplex) = length(s) - 1

function getvertex(s::Simplex, i::Int)
  @assert 1 <= i <= length(s)
  @assert issortedbyangle(s)
  return s[i]
end

function findabsvaluevertex(s::Simplex, op::T) where {T<:Function}
  _, index = op(map(v -> abs(value(v)), s))
  return s[index]
end

minabsvaluevertex(s::Simplex) = findabsvaluevertex(s, findmin)
maxabsvaluevertex(s::Simplex) = findabsvaluevertex(s, findmax)

function centroidignorevertex(f::T, s::Simplex, vertextoignore::Vertex
    ) where {T<:Function, U<:Function}
  g(v) = isequal(v, vertextoignore) ? zero(position(v)) : position(v)
  x = mapreduce(g, +, s) ./ (length(s) - 1)
  return Vertex(x, f(x))
end

centroid(s::Simplex) = mapreduce(position, +, s) ./ length(s)
function closestomiddlevertex(s::Simplex)
  mid = centroid(s)
  _, index = findmin(map(v->sum((position(v) - mid).^2), s))
  return s[index]
end

function swap!(s::Simplex, this::Vertex, forthat::Vertex)
  index = findfirst(x->isequal(x, this), s.vertices)
  s[index] = forthat
  @assert issortedbyangle(s)
  return nothing
end

function assessconvergence(simplex, config::NamedTuple)
  if abs(value(minabsvaluevertex(simplex))) <= config[:stopvalroot]
    return true, :STOPVAL_ROOT_REACHED
  elseif abs(value(maxabsvaluevertex(simplex))) >= config[:stopvalpole]
    return true, :STOPVAL_POLE_REACHED
  end

  toprocess = Set{Int}(1)
  processed = Set{Int}()
  while !isempty(toprocess)
    vi = pop!(toprocess)
    v = getvertex(simplex, vi)
    connectedto = Set{Int}()
    for (qi, q) ∈ enumerate(simplex)
      thisxtol = true
      for (i, (pv, pq)) ∈ enumerate(zip(position(v), position(q)))
        thisxtol &= isapprox(pv, pq, rtol=config[:xtol_rel][i],
                                     atol=config[:xtol_abs][i])
      end
      thisxtol && push!(connectedto, qi)
      thisxtol && for i in connectedto if i ∉ processed push!(toprocess, i) end end
    end
    push!(processed, vi)
  end
  allxtol = all(i ∈ processed for i ∈ 1:length(simplex))
  allxtol && return true, :XTOL_REACHED

  allftol = true
  for (vi, v) ∈ enumerate(simplex)
    for qi ∈ vi+1:length(simplex)
      q = getvertex(simplex, qi)
      allftol &= all(isapprox(value(v), value(q), rtol=config[:ftol_rel],
                              atol=config[:ftol_abs]))
    end
  end
  allftol && return true, :FTOL_REACHED

  return false, :CONTINUE
end

function _πtoπ(ϕ)
  ϕ < -π && return _πtoπ(ϕ + 2π)
  ϕ > π && return _πtoπ(ϕ - 2π)
  return ϕ
end

function windingangle(s::Simplex)
  return sum(_πtoπ.(angle.(value.(circshift(s.vertices, -1))) .-
                    angle.(value.(s.vertices))))
end

function windingnumber(s::Simplex)
  radians = windingangle(s)
  return isfinite(radians) ? Int64(round(radians / 2π)) : Int64(0)
end


text:
module BenchUnique

using BenchmarkTools
using Random: shuffle
using ThreadsX

n = 1_000_000
datasets = [
    # (label, setup)
    ("rand(1:10, $n)", rand(1:10, n)),
    ("rand(1:1000, $n)", rand(1:1000, n)),
]

suite = BenchmarkGroup()

for (label, xs) in datasets
    s = suite[label] = BenchmarkGroup()
    s["tx"] = @benchmarkable(ThreadsX.unique($xs))
    s["base"] = @benchmarkable(unique($xs))
end

end  # module
BenchUnique.suite


text:
module TestValidation

using CovidModels
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEProblem
using Petri
using Test
using SemanticModels.PetriModels

mkpath("outputs")

# states are [S, I, R]
Psir = PetriModel(Petri.Model(1:3, [
  ([1, 2], [2, 2]), # exposure
  ([2], [3]),     # recovery
]))

m = Psir.model

u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1

# probability of infection, recovery
β = [10 / sum(u0), 1 / 5]

tspan = (0, 100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg = Tsit5())

savedata(sol, "S,I,R", "outputs/sirdata.csv")

# states are [S, I, R, D]
Psird = PetriModel(Petri.Model(
  1:4,
  [
    ([1, 2], [2, 2]), # exposure
    ([2], [3]),     # recovery
    ([2], [4]),     # death
  ]
))

m = Psird.model

u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1

# probability of infection, recovery, death
β = [10 / sum(u0), 1 / 5, 1 / 10]

tspan = (0, 100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg = Tsit5())

savedata(sol, "S,I,R,D", "outputs/sirddata.csv")

# states are [S, I, E, R]
Pseir = PetriModel(Petri.Model(
  1:5,
  [
    ([1, 2], [3, 2]), # exposure
    ([3], [2]),     # onset
    ([2], [4]),     # recovery
  ],
))

m = Pseir.model

u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1

# probability of exposure, progression, recovery
β = [10 / sum(u0), 1 / 2, 1 / 5]

tspan = (0, 100.0)
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob, alg = Tsit5())

savedata(sol, "S,I,E,R", "outputs/seirdata.csv")

# states are [S, I, E, R, D]
Pseird = PetriModel(Petri.Model(
  1:5,
  [
    ([1, 2], [3, 2]), # exposure
    ([3], [2]),     # onset
    ([2], [4]),     # recovery
    ([2], [5]),     # death
  ],
))

m = Pseird.model

u0 = zeros(Float64, length(m.S))
u0[1] = 10000
u0[2] = 1

seirdparams() = begin
  # probability of exposure, progression, recovery, death
  βseird = [10 / sum(u0), 1 / 2, 1 / 5, 1 / 16]
end

tspan = (0, 100.0)
prob = ODEProblem(fluxes(m), u0, tspan, seirdparams())
sol = OrdinaryDiffEq.solve(prob, alg = Tsit5())

savedata(sol, "S,I,E,R,D", "outputs/seirddata.csv")
end


text:
# @doc "a DReal typed Expression" ->
immutable Ex{T}
  e::Ptr{Void}
  vars::Set{ASCIIString}
end

no_vars() = Set{ASCIIString}()

## TODO
## Check that lb <= ub
function safe_add!(name::ASCIIString, ctx::Context)
  name ∈ ctx.vars && error("variable $name already declared within this context")
  add_vars!(ctx,name)
end

# FIXME: Way too much duplication here

# @doc "Create a variable" ->
Var(ctx::Context, T::Type{Float64}, name::ASCIIString, lb::Float64, ub::Float64) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_real_var(ctx.ctx, name, lb, ub), Set([name])))
Var(ctx::Context, T::Type{Int64}, name::ASCIIString, lb::Int32, ub::Int32) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_int_var(ctx.ctx, name, Int64(lb), Int64(ub)), Set([name])))
Var(ctx::Context, T::Type{Int64}, name::ASCIIString, lb::Int64, ub::Int64) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_int_var(ctx.ctx, name, lb, ub), Set([name])))
Var(ctx::Context, T::Type{Int64}, name::ASCIIString) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_int_var(ctx.ctx, name, typemin(Int64), typemax(Int64)),Set([name])))
Var(ctx::Context,T::Type{Bool}, name::ASCIIString) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_bool_var(ctx.ctx, name),Set([name])))

Var(ctx::Context, T::Type{Float64}, name::ASCIIString) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_unbounded_real_var(ctx.ctx, name),Set([name])))

# Global Context
Var(T, name::ASCIIString, lb, ub) = Var(global_context(), T, name, lb, ub)
Var(T, name::ASCIIString) = Var(global_context(), T, name)

# Default Name
Var{R<:Real}(ctx::Context, T::Type, lb::R, ub::R) = Var(ctx, T, genvar(), lb, ub)
Var(ctx::Context, T::Type) = Var(ctx, T, genvar())
Var{R<:Real}(T::Type, lb::R, ub::R) = Var(T, genvar(), lb, ub)
Var(T::Type) = Var(T, genvar())

ForallVar(T::Type{Float64}, lb::Float64, ub::Float64; name::ASCIIString = genvar(), ctx::Context = global_context()) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_forall_real_var(ctx.ctx, name, lb, ub), Set([name])))

ForallVar(T::Type{Float64}; name::ASCIIString = genvar(), ctx::Context = global_context()) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_forall_unbounded_real_var(ctx.ctx, name), Set([name])))

ForallVar(T::Type{Int64}, lb::Int64, ub::Int64; name::ASCIIString = genvar(), ctx::Context = global_context()) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_forall_int_var(ctx.ctx, name, lb, ub), Set([name])))

ForallVar(T::Type{Int64}; name::ASCIIString = genvar(), ctx::Context = global_context()) =
  (safe_add!(name, ctx); Ex{T}(opensmt_mk_forall_int_var(ctx.ctx, name, typemin(Int), typemax(Int)), Set([name])))
## Printing
## ========
print(io::IO, e::Ex) = opensmt_print_expr(e.e)