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

a = parse(BigInt,"123456789012345678901234567890")
b = parse(BigInt,"123456789012345678901234567891")

@test typeof(a+1) == BigInt
@test a+1 == b
@test isequal(a+1, b)
@test b == a+1
@test !(b == a)
@test b > a
@test b >= a
@test !(b < a)
@test !(b <= a)

c = parse(BigInt,"246913578024691357802469135780")
@test typeof(a * 2) == BigInt
@test a*2 == c
@test c-a == a
@test c == a + a
@test c+1 == a+b

d = parse(BigInt,"-246913578024691357802469135780")
@test typeof(d) == BigInt
@test d == -c

ee = typemax(Int64)
@test typeof(BigInt(ee)) == BigInt
@test BigInt(ee)+1 == parse(BigInt,"9223372036854775808")

#Multiple calls for sanity check, since we're doing direct memory manipulation
@test string(a) == "123456789012345678901234567890"
@test string(b) == "123456789012345678901234567891"
@test string(c) == "246913578024691357802469135780"
@test string(d) == "-246913578024691357802469135780"
@test string(a) == "123456789012345678901234567890"

for i = -10:10, j = [-10:-1; 1:10]
    @test div(BigInt(i), BigInt(j)) == div(i,j)
    @test fld(BigInt(i), BigInt(j)) == fld(i,j)
    @test mod(BigInt(i), BigInt(j)) == mod(i,j)
    @test rem(BigInt(i), BigInt(j)) == rem(i,j)
end

@test typeof(BigInt(typemax(Int8))) == BigInt
@test typeof(BigInt(typemax(Int16))) == BigInt
@test typeof(BigInt(typemax(Int32))) == BigInt
@test typeof(BigInt(typemax(Int64))) == BigInt
@test typeof(BigInt(typemax(Int128))) == BigInt

@test typeof(BigInt(true)) == BigInt
@test typeof(BigInt(typemax(UInt8))) == BigInt
@test typeof(BigInt(typemax(UInt16))) == BigInt
@test typeof(BigInt(typemax(UInt32))) == BigInt
@test typeof(BigInt(typemax(UInt64))) == BigInt
@test typeof(BigInt(typemax(UInt128))) == BigInt

@test typeof(BigInt(BigInt(1))) == BigInt

for x in (Int16(0), 1, 3//4, big(5//6), big(9))
    @test big(typeof(x)) == typeof(big(x))
    @test big(typeof(complex(x, x))) == typeof(big(complex(x, x)))
end

# Signed addition
@test a+Int8(1) == b
@test a+Int16(1) == b
@test a+Int32(1) == b
@test a+Int64(1) == b
@test a+Int128(1) == b
@test Int8(1)+ a == b
@test Int16(1)+a == b
@test Int32(1)+a == b
@test Int64(1)+a == b
@test b+Int8(-1) == a
@test b+Int16(-1) == a
@test b+Int32(-1) == a
@test b+Int64(-1) == a
@test Int8(-1)+ b == a
@test Int16(-1)+b == a
@test Int32(-1)+b == a
@test Int64(-1)+b == a

# Unsigned addition
@test a+true == b
@test a+UInt8(1) == b
@test a+UInt16(1) == b
@test a+UInt32(1) == b
@test a+UInt64(1) == b
@test a+UInt128(1) == b
@test true+a == b
@test UInt8(1)+ a == b
@test UInt16(1)+a == b
@test UInt32(1)+a == b
@test UInt64(1)+a == b

# Signed subtraction
@test b-Int8(1) == a
@test b-Int16(1) == a
@test b-Int32(1) == a
@test b-Int64(1) == a
@test Int8(1)- b == -a
@test Int16(1)-b == -a
@test Int32(1)-b == -a
@test Int64(1)-b == -a
@test a-Int8(-1) == b
@test a-Int16(-1) == b
@test a-Int32(-1) == b
@test a-Int64(-1) == b
@test Int8(-1)- a == -b
@test Int16(-1)-a == -b
@test Int32(-1)-a == -b
@test Int64(-1)-a == -b

# Unsigned subtraction
@test b-true == a
@test b-UInt8(1) == a
@test b-UInt16(1) == a
@test b-UInt32(1) == a
@test b-UInt64(1) == a
@test true-b == -a
@test UInt8(1)- b == -a
@test UInt16(1)-b == -a
@test UInt32(1)-b == -a
@test UInt64(1)-b == -a

# Signed multiplication
@test a*Int8(1) == a
@test a*Int16(1) == a
@test a*Int32(1) == a
@test a*Int64(1) == a
@test Int8(1)* a == a
@test Int16(1)*a == a
@test Int32(1)*a == a
@test Int64(1)*a == a
@test a*Int8(-1) == -a
@test a*Int16(-1) == -a
@test a*Int32(-1) == -a
@test a*Int64(-1) == -a
@test Int8(-1)* a == -a
@test Int16(-1)*a == -a
@test Int32(-1)*a == -a
@test Int64(-1)*a == -a

# Unsigned multiplication
@test a*true == a
@test a*UInt8(1) == a
@test a*UInt16(1) == a
@test a*UInt32(1) == a
@test a*UInt64(1) == a
@test true*a == a
@test UInt8(1)* a == a
@test UInt16(1)*a == a
@test UInt32(1)*a == a
@test UInt64(1)*a == a

@test a+BigInt(1) == b

@test BigInt(5) << 3 == 40
@test BigInt(5) >> 1 == 2
@test BigInt(-5) << 3 == -40
@test BigInt(-5) >> 1 == -3
@test BigInt(5) >> -3 == 40
@test BigInt(5) << -1 == 2
@test BigInt(-5) >> -3 == -40
@test BigInt(-5) << -1 == -3

@test ~BigInt(123) == -124
@test BigInt(123) & BigInt(234) == 106
@test BigInt(123) | BigInt(234) == 251
@test BigInt(123) ⊻ BigInt(234) == 145

@test gcd(BigInt(48), BigInt(180)) == 12
@test lcm(BigInt(48), BigInt(180)) == 720

@test factorial(BigInt(40)) == parse(BigInt,"815915283247897734345611269596115894272000000000")
@test binomial(BigInt(1), -1) == BigInt(0)
@test binomial(BigInt(1), 2)  == BigInt(0)
@test binomial(BigInt(-53), 42) == parse(BigInt,"959509335087854414441273718")
@test binomial(BigInt(113), BigInt(42)) == parse(BigInt,"18672199984318438125634054194360")

a = rand(1:100, 10000)
b = map(BigInt, a)
@test sum(a) == sum(b)

# Iterated arithmetic
a = parse(BigInt,"315135")
b = parse(BigInt,"12412")
c = parse(BigInt,"3426495623485904783478347")
d = parse(BigInt,"-1398984130")
f = parse(BigInt,"2413804710837418037418307081437315263635345357386985747464")
g = parse(BigInt,"-1")

@test +(a, b) == parse(BigInt,"327547")
@test +(a, b, c) == parse(BigInt,"3426495623485904783805894")
@test +(a, b, c, d) == parse(BigInt,"3426495623485903384821764")
@test +(a, b, c, d, f) == parse(BigInt,"2413804710837418037418307081437318690130968843290370569228")
@test +(a, b, c, d, f, g) == parse(BigInt,"2413804710837418037418307081437318690130968843290370569227")

@test *(a, b) == parse(BigInt,"3911455620")
@test *(a, b, c) == parse(BigInt,"13402585563389346256121263521460140")
@test *(a, b, c, d) == parse(BigInt,"-18750004504148804423388563022070650287578200")
@test *(a, b, c, d, f) == parse(BigInt,"-45258849200337190631492857400003938881995610529251881450243326128168934937055005474972396281351684800")
@test *(a, b, c, d, f, g) == parse(BigInt,"45258849200337190631492857400003938881995610529251881450243326128168934937055005474972396281351684800")

@test xor(a, b) == parse(BigInt,"327299")
@test xor(a, b, c) == parse(BigInt,"3426495623485904783798472")
@test xor(a, b, c, d) == parse(BigInt,"-3426495623485906178489610")
@test xor(a, b, c, d, f) == parse(BigInt,"-2413804710837418037418307081437316711364709261074607933698")
@test xor(a, b, c, d, f, g) == parse(BigInt,"2413804710837418037418307081437316711364709261074607933697")

@test (&)(a, b) == parse(BigInt,"124")
@test (&)(a, b, c) == parse(BigInt,"72")
@test (&)(a, b, c, d) == parse(BigInt,"8")
@test (&)(a, b, c, d, f) == parse(BigInt,"8")
@test (&)(a, b, c, d, f, g) == parse(BigInt,"8")

@test (|)(a, b) == parse(BigInt,"327423")
@test (|)(a, b, c) == parse(BigInt,"3426495623485904783802111")
@test (|)(a, b, c, d) == parse(BigInt,"-1396834561")
@test (|)(a, b, c, d, f) == parse(BigInt,"-1358954753")
@test (|)(a, b, c, d, f, g) == parse(BigInt,"-1")

@test trailing_ones(a) == 8
@test trailing_zeros(b) == 2
@test count_ones(a) == 14

# Large Fibonacci to exercise BigInt
# from Bill Hart, https://groups.google.com/group/julia-dev/browse_frm/thread/798e2d1322daf633
function mul(a::Vector{BigInt}, b::Vector{BigInt})
   x = a[2]*b[2]
   c = Array{BigInt,1}(3)
   c[1] = a[1]*b[1] + x
   c[2] = a[1]*b[2] + a[2]*b[3]
   c[3] = x + a[3]*b[3]
   return c
end

function bigfib(n)
    n == 0 && return BigInt(0)
    n -= 1
    r = [BigInt(1), BigInt(1), BigInt(0)]
    s = [BigInt(1), BigInt(0), BigInt(1)]
    while true
       (n & 1) == 1 && (s = mul(s,r))
       (n >>= 1) == 0 && return s[1]
       r = mul(r,r)
    end
end
@test [bigfib(n) for n=0:10] == [0,1,1,2,3,5,8,13,21,34,55]

let s, n = bigfib(1000001)
    @test ndigits(n) == 208988
    @test mod(n,big(10)^15) == 359244926937501
    @test div(n,big(10)^208973) == 316047687386689

    s = string(n)
    @test length(s) == 208988
    @test endswith(s, "359244926937501")
    @test startswith(s, "316047687386689")
end

# serialization (#5133)
let n = parse(BigInt, "359334085968622831041960188598043661065388726959079837"),
    b = IOBuffer()
    serialize(b, n)
    seek(b, 0)
    @test deserialize(b) == n
end

# issue #5873
@test ndigits(big(90)) == 2
@test ndigits(big(99)) == 2
ndigits_mismatch(n) = ndigits(n) != ndigits(BigInt(n))
@test !any(ndigits_mismatch, 8:9)
@test !any(ndigits_mismatch, 64:99)
@test !any(ndigits_mismatch, 512:999)
@test !any(ndigits_mismatch, 8192:9999)

# The following should not crash (#16579)
ndigits(big(rand(Int)), rand(63:typemax(Int)))
ndigits(big(rand(Int)), big(2)^rand(2:999))

for x in big.([-20:20; rand(Int)])
    for _base in -1:1
        @test_throws DomainError ndigits(x, _base)
    end
end

@test Base.ndigits0zpb(big(0), big(rand(2:100))) == 0

# digits with BigInt bases (#16844)
@test digits(big(2)^256, big(2)^128) == [0, 0, 1]

# conversion from float
@test BigInt(2.0) == BigInt(2.0f0) == BigInt(big(2.0)) == 2
@test_throws InexactError convert(BigInt, 2.1)
@test_throws InexactError convert(BigInt, big(2.1))

# issue #13367
@test trunc(BigInt,2.1) == 2
@test round(BigInt,2.1) == 2
@test floor(BigInt,2.1) == 2
@test ceil(BigInt,2.1) == 3

@test trunc(BigInt,2.1f0) == 2
@test round(BigInt,2.1f0) == 2
@test floor(BigInt,2.1f0) == 2
@test ceil(BigInt,2.1f0) == 3

@test_throws InexactError trunc(BigInt,Inf)
@test_throws InexactError round(BigInt,Inf)
@test_throws InexactError floor(BigInt,Inf)
@test_throws InexactError ceil(BigInt,Inf)

@test bin(big(3)) == "11"
@test oct(big(9)) == "11"
@test oct(-big(9)) == "-11"
@test hex(big(12)) == "c"

# Issue #18849: bin, oct, dec, hex should not call sizeof on BigInts
# when padding is desired
let padding = 4, low = big(4), high = big(2^20)
    @test bin(low, padding) == "0100"
    @test oct(low, padding) == "0004"
    @test dec(low, padding) == "0004"
    @test hex(low, padding) == "0004"

    @test bin(high, padding) == "100000000000000000000"
    @test oct(high, padding) == "4000000"
    @test dec(high, padding) == "1048576"
    @test hex(high, padding) == "100000"

    @test bin(-low, padding) == "-0100" # handle negative numbers correctly
    @test oct(-low, padding) == "-0004"
    @test dec(-low, padding) == "-0004"
    @test hex(-low, padding) == "-0004"

    @test bin(-high, padding) == "-100000000000000000000"
    @test oct(-high, padding) == "-4000000"
    @test dec(-high, padding) == "-1048576"
    @test hex(-high, padding) == "-100000"
end

# respect 0-padding on big(0)
for f in (bin, oct, dec, hex)
    local f
    @test f(big(0), 0) == ""
end
@test base(rand(2:62), big(0), 0) == ""

@test isqrt(big(4)) == 2
@test isqrt(big(5)) == 2

@test big(5)^true == big(5)
@test big(5)^false == one(BigInt)

# operations that when applied to Int64 give Float64, should give BigFloat
@test typeof(exp(a)) == BigFloat
@test typeof(exp2(a)) == BigFloat
@test typeof(exp10(a)) == BigFloat
@test typeof(expm1(a)) == BigFloat
@test typeof(cosh(a)) == BigFloat
@test typeof(sinh(a)) == BigFloat
@test typeof(tanh(a)) == BigFloat
@test typeof(sech(a)) == BigFloat
@test typeof(csch(a)) == BigFloat
@test typeof(coth(a)) == BigFloat
@test typeof(cbrt(a)) == BigFloat
@test typeof(tan(a)) == BigFloat
@test typeof(cos(a)) == BigFloat
@test typeof(sin(a)) == BigFloat

@test BigInt <: Signed
@test big(1) isa Signed

let x = big(1)
    @test signed(x) === x
    @test convert(Signed, x) === x
    @test Signed(x) === x
    @test_throws MethodError convert(Unsigned, x) # could change in the future
end

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
include("01_read_data.jl")
include("02_set_options.jl")
include("03_parallel_run.jl")

ext: jl


text:
using BSON: @save, @load
using GeometryBasics
using CairoMakie  # GLMakie
using SLAM
using ProgressMeter
using FileIO

include("kitty.jl")

function main(n_frames)
    base_dir = "/mnt/zy_data/kitty-dataset/"
    sequence = "00"
    stereo = true

    dataset = KittyDataset(base_dir, sequence; stereo)
    println(dataset)

    save_dir = joinpath("/mnt/zy_data/kitty-dataset/projects", "kitty-$sequence")
    isdir(save_dir) || mkdir(save_dir)
    @info "Save directory: $save_dir"

    fx, fy = dataset.K[1, 1], dataset.K[2, 2]
    cx, cy = dataset.K[1:2, 3]
    height, width = 376, 1241
    # height, width = 370, 1226

    camera = SLAM.Camera(fx, fy, cx, cy, 0, 0, 0, 0, height, width)
    right_camera = SLAM.Camera(fx, fy, cx, cy, 0, 0, 0, 0, height, width; Ti0=dataset.Ti0)

    params = Params(;
        stereo,
        window_size=9, max_distance=35, pyramid_levels=3,
        max_nb_keypoints=1000, max_reprojection_error=3.0,
        do_local_bundle_adjustment=false, map_filtering=true)

    saver = ReplaySaver()
    visualizer = nothing
    # visualizer = Visualizer((900, 600))
    # display(visualizer)

    slam_manager = SlamManager(params, camera; right_camera, visualizer=saver)
    slam_manager_thread = Threads.@spawn run!(slam_manager)

    t1 = time()
    @showprogress for i in 1:n_frames
        timestamp = dataset.timestamps[i]
        left_frame, right_frame = dataset[i]
        left_frame = Gray{Float64}.(left_frame)

        if params.stereo
            right_frame = Gray{Float64}.(right_frame)
            add_stereo_image!(slam_manager, left_frame, right_frame, timestamp)
        else
            add_image!(slam_manager, left_frame, timestamp)
        end

        if visualizer ≢ nothing
            SLAM.set_image!(visualizer, rotr90(left_frame))
            process_frame_wc!(visualizer)
        end

        q_size = get_queue_size(slam_manager)
        f_size = length(slam_manager.mapper.estimator.frame_queue)
        m_size = length(slam_manager.mapper.keyframe_queue)
        while q_size > 0 || f_size > 0 || m_size > 0
            sleep(1e-2)
            q_size = get_queue_size(slam_manager)
            f_size = length(slam_manager.mapper.estimator.frame_queue)
            m_size = length(slam_manager.mapper.keyframe_queue)
        end

        sleep(1e-2)
    end

    slam_manager.exit_required = true
    wait(slam_manager_thread)

    t2 = time()
    @info "SLAM took: $(t2 - t1) seconds."

    SLAM.save(saver, save_dir)
    slam_manager, visualizer
end

function replay(n_frames)
    base_dir = "/mnt/zy_data/kitty-dataset/"
    sequence = "00"
    stereo = true

    dataset = KittyDataset(base_dir, sequence; stereo)
    println(dataset)

    save_dir = joinpath("/mnt/zy_data/kitty-dataset/projects", "kitty-$sequence")
    isdir(save_dir) || mkdir(save_dir)
    @info "Save directory: $save_dir"

    saver = ReplaySaver()
    SLAM.load!(saver, save_dir)
    @assert length(saver.positions) == n_frames - 1

    resolution = (900, 600)
    image_resolution = (1241, 376)
    # image_resolution = (1226, 370)
    visualizer = Visualizer(;resolution, image_resolution)
    display(visualizer)

    t1 = time()
    for i in 2:n_frames
        left_frame, right_frame = dataset[i]
        left_frame = Gray{Float64}.(left_frame)

        position = saver.positions[i - 1]
        set_image!(visualizer, rotr90(left_frame))
        set_position!(visualizer, position)

        sleep(1 / 60)
    end
    t2 = time()
    @info "Visualization took: $(t2 - t1) seconds."

    visualizer
    save(joinpath(save_dir, "out_1.jpg"), visualizer.image)  # image  figure
    display(visualizer)
end

# main(4541)
replay(4541)

# 可正常运行, 就是plot和save_plot时没成功.
# http://www.cvlibs.net/datasets/kitti/eval_odometry.php   vslam, kitti


text:
#---------------------------------------------------------------------#
#This function computes the norm of 2 vectors
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function compute_norm(qn,qe,Npoin,DFloat)
    norm2=norm(qn-qe,2)/norm(qe,2)
    return (norm2)
end


text:
using CompressedStacks
using Base.Test

# write your own tests here
@test 1 == 1

# IO test
include("temp.jl")

println("Test IO")
io_test()


text:
module AtomDescribe

using LinearAlgebra
using AtomBase
using XYZFiles

abstract type AbstractDescribe end
abstract type MoleculeDescribe <: AbstractDescribe end

export AbstractDescribe, MoleculeDescribe

function getCharge(symbol::String)
    if symbol == "H"
        return 1.0
    elseif symbol == "C"
        return 12.0
    elseif symbol == "N"
        return 14.0
    elseif symbol == "O"
        return 16.0
    end
    #TODO: Throw unknown chemical symbol.
end

function dist(atom1::T, atom2::T) where {T<:AbstractAtom}
    distanceSquare = 0.0
    for i in 1:3
        distanceSquare += (atom1.position[i]-atom2.position[i])^2
    end
    return sqrt(distanceSquare)
end

include("fileFormat.jl")
include("coulomb.jl")

function formatRead(fileFormat::xyzFormat)
    return XYZFiles.read(fileFormat.fileName, ["energy",])
end

end # module


text:
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport GLS_Apr_weekday_PM_ext

xi_list = GLS_Apr_weekday_PM_ext.xi_list
P = GLS_Apr_weekday_PM_ext.P
L = GLS_Apr_weekday_PM_ext.L  # dimension of xi

println("dimensions of P are: ")
println(size(P, 1), size(P, 2))
println("--------------------")

# println("xi_list is: ")
# println(xi_list)
println("sum of xi's is: ")
println(sum(xi_list))
println("max of xi's is: ")
println(maximum(xi_list))
println("--------------------")

using JuMP

mGLSJulia = Model()

@defVar(mGLSJulia, lam[1:size(P,1)] >= 0)

@defVar(mGLSJulia, p[1:size(P,1), 1:size(P,2)] >= 0)

for i = 1:size(P,1)
    for j = 1:size(P,2)
        if P[i,j] == 0
            @addConstraint(mGLSJulia, p[i,j] == 0)
        end
    end
end

for i = 1:size(P,1)
    @addNLConstraint(mGLSJulia, sum{p[i,j], j = 1:size(P,2)} == 1)
end

for l = 1:L
    @addNLConstraint(mGLSJulia, sum{p[i,l] * lam[i], i = 1:size(P,1)} == xi_list[l])
end

@setNLObjective(mGLSJulia, Min, sum{p[1,j], j = 1:size(P,2)})  # play no actual role, but could not use zero objective

solve(mGLSJulia)

println("lam is: ")
println(getValue(lam))
println("--------------------")

println("objective value is: ")
println(getObjectiveValue(mGLSJulia))
println("--------------------")
println("End")

GLS_Apr_weekday_PM_ext.saveDemandVec(getValue(lam))

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
#### Functions to simulate the Rosenzweig-MacArthur models with red noise
#### "Slow organisms exhibit sudden population disappearances in a reddened world" by Greyson-Gaito, Gellner, & McCann.

# Note - Increase the number of threads by going to Julia Extension settings in VS Code or setting the environment variable in bash
include("packages.jl")
include("slowfast_commoncode.jl")
include("slowfast_canardfinder.jl")

function prop_canard_rednoise(model, ep, effR0, r, reps)
    count_canard = 0
    count_axial = 0
    count_nothing = 0
    tsend = 6000.0
    for i in 1:reps
        pattern = cf_returnmap(model, ep, effR0, 1, r, i, tsend, 0.0:1.0:tsend)
        if pattern == "canard"
            count_canard += 1
        elseif pattern == "axial"
            count_axial += 1
        else
            count_nothing += 1
        end
    end
    return vcat(count_canard / reps, count_axial / reps, count_nothing / reps)
end #function to run reps number of simulations and count the number of simulations exhibiting a quasicanard, axial solution, or nothing.

function prop_canard_rednoise_data(model, ep, effR0, reps)
    r_range = 0.0:0.01:0.9
    data = Vector{Vector{Float64}}(undef,length(r_range))
    @threads for i in eachindex(r_range)
        @inbounds data[i] =  prop_canard_rednoise(model, ep, effR0, r_range[i], reps)
    end

    return prep_data(data, r_range)
end #function to parallelize quasicanard simulations over multiple values of noise correlation

### Rosenzweig-MacArthur model
begin
    rn_ep0079eff05_RozMac = prop_canard_rednoise_data("RozMac", 0.079, 0.5, 1000)
    CSV.write(joinpath(abpath(), "data/rn_ep0079eff05_RozMac.csv"), rn_ep0079eff05_RozMac)
    rn_ep0079eff065_RozMac = prop_canard_rednoise_data("RozMac",0.079, 0.65, 1000)
    CSV.write(joinpath(abpath(),"data/rn_ep0079eff065_RozMac.csv"), rn_ep0079eff065_RozMac)
    rn_ep0079eff07_RozMac = prop_canard_rednoise_data("RozMac",0.079, 0.7, 1000)
    CSV.write(joinpath(abpath(), "data/rn_ep0079eff07_RozMac.csv"), rn_ep0079eff07_RozMac)
end

begin
    rn_ep0004eff05_RozMac = prop_canard_rednoise_data("RozMac",0.004, 0.5, 1000)
    CSV.write(joinpath(abpath(), "data/rn_ep0004eff05_RozMac.csv"), rn_ep0004eff05_RozMac)
    rn_ep0004eff065_RozMac = prop_canard_rednoise_data("RozMac",0.004, 0.65, 1000)
    CSV.write(joinpath(abpath(), "data/rn_ep0004eff065_RozMac.csv"), rn_ep0004eff065_RozMac)
    rn_ep0004eff07_RozMac = prop_canard_rednoise_data("RozMac",0.004, 0.7, 1000)
    CSV.write(joinpath(abpath(), "data/rn_ep0004eff07_RozMac.csv"), rn_ep0004eff07_RozMac)
end


text:
export run_ode_test, run_ode_testsuite

@doc "Runs one ODE TestRun" ->
function run_ode_test{Name}(tr::TestRun{Name}; verbose=true)
    if !isapplicable(tr.tc, tr.solver)
        println("\nTest case $Name not compatible with solver $(tr.solver.solverfn)")
        return nothing
    end

    out, walltime, mem, gc_time = (1,2,3,4) # try makes a new scope?!
    gc()
    try
        out, walltime, mem, gc_time = @timed tr.solver.wrapper(tr)
    catch e
        if verbose
            println("Error occurred: $e")
        end
        return simple_eval(e, tr)
    end
    tend, yend, stats = out
    return simple_eval(tend, yend, stats, walltime, mem, gc_time, tr)
end

function run_ode_test_throwerror{Name}(tr::TestRun{Name})
    if !isapplicable(tr.tc, tr.solver)
        println("\nTest case $Name not compatible with solver $(tr.solver.solverfn)")
        return nothing
    end
    gc()
    out, walltime, mem, gc_time = @timed tr.solver.wrapper(tr)
    tend, yend, stats = out
    return simple_eval(tend, yend, stats, walltime, mem, gc_time, tr)
end

@doc "Runs a suite of ODE TestRun's" ->
function run_ode_testsuite{Name}(suite::TestSuite{Name}; verbose=true, warmup=true, throwerror=false)
    run_fn = throwerror ? run_ode_test_throwerror : run_ode_test
    if !isapplicable(suite.tc, suite.solver)
        println("\nTest case $Name not compatible with solver $(suite.solver.solverfn)")
        return nothing
    end
    if verbose
        println("\nRunning test case $Name with solver $(suite.solver.solverfn)")
    end
    if warmup
        # warmup compiling the function to get accurate timings
        run_fn(first(suite); verbose=false)
    end

    out = TestResults[]
    tot = length(suite)
    for (i,tr) in enumerate(suite)
        if verbose
            print("Running test $i of $tot:")
        end
        push!(out, run_fn(tr; verbose=verbose))
        if verbose
            sd = out[end].scd
            wt = out[end].walltime
            mem = out[end].mem
            if mem<1e6
                println(" sig. digits= $sd, walltime= $(wt)s, memory= $(mem)bytes")
            else
                mem = round(mem/1e6)
                println(" sig. digits= $sd, walltime= $(wt)s, memory= $(mem)MB")
            end
        end
    end
    out
end


text:
# memory operations
# TODO: Properly use dispatch capabilities for these functions

# TODO: cmdbuf cleanup - Was running into errors
function sync_gpu_to_cpu!(dev::MtlDevice, buf::MtlBuffer{T}) where T
    cmdbuf = MtlCommandBuffer(global_queue(dev))
    MtlBlitCommandEncoder(cmdbuf) do enc
        MTL.append_sync!(enc, buf)
    end
    commit!(cmdbuf)
    wait_completed(cmdbuf)
end

# GPU -> GPU
Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, src::MtlBuffer{T}, N::Integer) where T =
    unsafe_copyto!(dev, dst, 1, src, 1, n)
function Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, doff::Integer,
                             src::MtlBuffer{T}, soff::Integer, N::Integer) where T
    cmdbuf = MtlCommandBuffer(global_queue(dev))
    MtlBlitCommandEncoder(cmdbuf) do enc
        MTL.append_copy!(enc, dst, (doff-1) * sizeof(T), src, (soff - 1) * sizeof(T), N * sizeof(T))
    end
    commit!(cmdbuf)
    wait_completed(cmdbuf)
end

# GPU -> CPU
Base.unsafe_copyto!(dev::MtlDevice, dst::Ptr{T}, src::MtlBuffer{T}, N::Integer) where T =
    unsafe_copyto!(dev, dst, src, 1, n)
function Base.unsafe_copyto!(dev::MtlDevice, dst::Ptr{T},
                             src::MtlBuffer{T}, soff::Integer, N::Integer) where T
    storage_type = src.storageMode
    if storage_type ==  MTL.MtStorageModePrivate
        tmp_buf = alloc(T, dev, N, storage=Shared)
        Base.unsafe_copyto!(dev, tmp_buf, 1, src, soff, N)
        Base.unsafe_copyto!(dst, content(tmp_buf), N)
        free(tmp_buf)
    elseif storage_type ==  MTL.MtStorageModeShared
        Base.unsafe_copyto!(dst, content(src, soff), N)
    elseif storage_type ==  MTL.MtStorageModeManaged
        sync_gpu_to_cpu!(dev, src)
        Base.unsafe_copyto!(dst, content(src, soff), N)
    end
    return dst
end

# CPU -> GPU
Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, src::Ptr{T}, N::Integer) where T =
    unsafe_copyto!(dev, dst, 1, src, n)
function Base.unsafe_copyto!(dev::MtlDevice, dst::MtlBuffer{T}, doff::Integer,
                             src::Ptr{T}, N::Integer) where T
    storage_type = dst.storageMode
    if storage_type == MTL.MtStorageModePrivate
        tmp_buf = alloc(T, dev, N, src, storage=Shared)
        Base.unsafe_copyto!(dev, tmp_buf, src, N)
        Base.unsafe_copyto!(dev, dst, doff, tmp_buf, 1, N)
        free(tmp_buf)
    elseif storage_type == MTL.MtStorageModeShared
        Base.unsafe_copyto!(content(dst, doff), src, N)
    elseif storage_type == MTL.MtStorageModeManaged
        Base.unsafe_copyto!(content(dst, doff), src, N)
        MTL.DidModifyRange!(dst, 1:N)
    end
    return dst
end

function unsafe_fill!(dev::MtlDevice, ptr::MtlBuffer{T}, value::Union{UInt8,Int8}, N::Integer) where T
    cmdbuf = MtlCommandBuffer(global_queue(dev))
    MtlBlitCommandEncoder(cmdbuf) do enc
        MTL.append_fillbuffer!(enc, ptr, value, N * sizeof(T))
    end
    commit!(cmdbuf)
    wait_completed(cmdbuf)
end

# TODO: Implement generic fill since mtBlitCommandEncoderFillBuffer is limiting


text:
#=
This is a port of H5Zlz4.c to Julia
https://github.com/HDFGroup/hdf5_plugins/blob/master/LZ4/src/H5Zlz4.c
https://github.com/nexusformat/HDF5-External-Filter-Plugins/blob/master/LZ4/src/H5Zlz4.c
https://github.com/silx-kit/hdf5plugin/blob/main/src/LZ4/H5Zlz4.c

H5Zlz4 is originally a copyright of HDF Group. License: licenses/H5Zlz4_LICENSE.txt

The following license applies to the Julia port.
Copyright (c) 2021 Mark Kittisopikul and Howard Hughes Medical Institute. License MIT, see LICENSE.txt
=#
module H5Zlz4

using CodecLz4
using HDF5.API
import HDF5.Filters: Filter, filterid, register_filter, filtername, filter_func, filter_cfunc

export H5Z_FILTER_LZ4, H5Z_filter_lz4, Lz4Filter

const H5Z_FILTER_LZ4 = API.H5Z_filter_t(32004)

const DEFAULT_BLOCK_SIZE = 1 << 30;
const lz4_name = "HDF5 lz4 filter; see http://www.hdfgroup.org/services/contributions.html"

const LZ4_AGGRESSION = Ref(1)

# flags H5Z_FLAG_REVERSE or H5Z_FLAG_OPTIONAL
# cd_nelmts number of elements in cd_values (0 or 1)
# cd_values the first optional element must be the blockSize
# nbytes - number of valid bytes of data
# buf_size - total size of buffer
# buf - pointer to pointer of data
function H5Z_filter_lz4(flags::Cuint, cd_nelmts::Csize_t,
                        cd_values::Ptr{Cuint}, nbytes::Csize_t,
                        buf_size::Ptr{Csize_t}, buf::Ptr{Ptr{Cvoid}})::Csize_t

    outBuf = C_NULL
    ret_value = Csize_t(0)

    try

    if (flags & API.H5Z_FLAG_REVERSE) != 0 # reverse filter, decompressing
        #i32Buf = Ref{UInt32}()
        blockSize = UInt32(0)
        roBuf = Ref{UInt8}()
        rpos = Ptr{UInt8}(unsafe_load(buf))
        #i64Buf = Ptr{UInt64}(rpos)
        # Load the first 8 bytes from buffer as a big endian UInt64
        # This is the original size of the buffer
        origSize = ntoh(unsafe_load(Ptr{UInt64}(rpos)))
        rpos += 8 # advance the pointer

        # Next read the next four bytes from the buffer as a big endian UInt32
        # This is the blocksize
        #i32Buf[] = rpos
        blockSize = ntoh(unsafe_load(Ptr{UInt32}(rpos)))
        rpos += 4
        if blockSize > origSize
            blockSize = origSize
        end

        # malloc a byte buffer of origSize
        # outBuf = Vector{UInt8}(undef, origSize)
        @debug "OrigSize" origSize
        outBuf = Libc.malloc(origSize)
        # Julia should throw an error if it cannot allocate this
        roBuf = Ptr{UInt8}(outBuf)
        decompSize = 0
        # Start with the first blockSize
        while decompSize < origSize
            # compressedBlockSize = UInt32(0)
            if origSize - decompSize < blockSize # the last block can be smaller than block size
                blockSize = origSize - decompSize
            end

            #i32Buf[] = rpos
            compressedBlockSize = ntoh(unsafe_load(Ptr{UInt32}(rpos)))
            rpos += 4

            if compressedBlockSize == blockSize
                # There was no compression
                # memcpy(roBuf, rpos, blockSize)
                unsafe_copyto!(roBuf, rpos, blockSize)
                decompressedBytes = blockSize
            else
                # do the compression
                # LZ4_decompress_fast, version number 10300 ?
                @debug "decompress_safe" rpos roBuf compressedBlockSize (origSize - decompSize)
                decompressedBytes = CodecLz4.LZ4_decompress_safe(rpos, roBuf, compressedBlockSize, origSize -decompSize)
                @debug "decompressedBytes" decompressedBytes
            end

            rpos += compressedBlockSize
            roBuf += blockSize
            decompSize += decompressedBytes
        end
        Libc.free(unsafe_load(buf))
        unsafe_store!(buf, outBuf)
        outBuf = C_NULL
        ret_value = Csize_t(origSize)
    else
        # forward filter
        # compressing
        #i64Buf = Ref{UInt64}()
        #i32Buf = Ref{UInt32}()

        if nbytes > typemax(Int32)
            error("Can only compress chunks up to 2GB")
        end
        blockSize = unsafe_load(cd_values)
        if cd_nelmts > 0 && blockSize > 0
        else
            blockSize = DEFAULT_BLOCK_SIZE
        end
        if blockSize > nbytes
            blockSize = nbytes
        end
        nBlocks = (nbytes-1) ÷ blockSize + 1
        maxDestSize = nBlocks * CodecLz4.LZ4_compressBound(blockSize) +  4 + 8 + nBlocks*4
        outBuf = Libc.malloc(maxDestSize)

        rpos = Ptr{UInt8}(unsafe_load(buf))
        roBuf = Ptr{UInt8}(outBuf)

        # Header
        unsafe_store!(Ptr{UInt64}(roBuf), hton(UInt64(nbytes)))
        roBuf += 8

        unsafe_store!(Ptr{UInt32}(roBuf), hton(UInt32(blockSize)))
        roBuf += 4

        outSize = 12

        for block = 0:nBlocks-1
            # compBlockSize::UInt32
            origWritten = Csize_t(block*blockSize)
            if nbytes - origWritten < blockSize # the last block may be < blockSize
                blockSize = nbytes - origWritten
            end

            # aggression = 1 is the same LZ4_compress_default
            @debug "LZ4_compress_fast args" rpos outBuf roBuf roBuf+4 blockSize nBlocks CodecLz4.LZ4_compressBound(blockSize)
            compBlockSize = UInt32(CodecLz4.LZ4_compress_fast(rpos, roBuf+4, blockSize, CodecLz4.LZ4_compressBound(blockSize), LZ4_AGGRESSION[]))
            @debug "Compressed block size" compBlockSize

            if compBlockSize == 0
                error("Could not compress block $block")
            end

            if compBlockSize >= blockSize # compression did not save any space, do a memcpy instead
                compBlockSize = blockSize
                unsafe_copyto!(roBuf+4, rpos, blockSize)
            end

            unsafe_store!(Ptr{UInt32}(roBuf), hton(UInt32(compBlockSize))) # write blocksize
            roBuf += 4

            rpos += blockSize
            roBuf += compBlockSize
            outSize += compBlockSize + 4
        end

        Libc.free(unsafe_load(buf))
        unsafe_store!(buf, outBuf)
        unsafe_store!(buf_size, outSize)
        outBuf = C_NULL
        ret_value = Csize_t(outSize)
    end # (flags & API.H5Z_FLAG_REVERSE) != 0

    catch err
        #  "In the case of failure, the return value is 0 (zero) and all pointer arguments are left unchanged."
        ret_value = Csize_t(0)
        @error "H5Zlz4.jl Non-Fatal ERROR: " err
        display(stacktrace(catch_backtrace()))
    finally
        if outBuf != C_NULL
            Libc.free(outBuf)
        end
    end
    return Csize_t(ret_value)
end

# Filters Module

"""
    Lz4Filter(blockSize)

Apply LZ4 compression. `blockSize` is the main argument. The filter id is $H5Z_FILTER_LZ4.

# External Links
* [LZ4 HDF5 Filter ID 32004](https://portal.hdfgroup.org/display/support/Filters#Filters-32004)
* [LZ4 HDF5 Plugin Repository (C code)](https://github.com/nexusformat/HDF5-External-Filter-Plugins/tree/master/LZ4)
"""
struct Lz4Filter <: Filter
    blockSize::Cuint
end
Lz4Filter() = Lz4Filter(DEFAULT_BLOCK_SIZE)

filterid(::Type{Lz4Filter}) = H5Z_FILTER_LZ4
filtername(::Type{Lz4Filter}) = lz4_name
filter_func(::Type{Lz4Filter}) = H5Z_filter_lz4
filter_cfunc(::Type{Lz4Filter}) = @cfunction(H5Z_filter_lz4, Csize_t,
                                             (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                                             Ptr{Csize_t}, Ptr{Ptr{Cvoid}}))

function __init__()
    register_filter(Lz4Filter)
end

end


text:
push!(LOAD_PATH, string(@__DIR__, "/.."))
using GraphMatFun
using Documenter

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

makedocs(;
    modules = [GraphMatFun],
    authors = "Elias Jarlebring <jarlebring@gmail.com> and contributors",
    repo = "https://github.com/matrixfunctions/GraphMatFun.jl/" *
           "blob/{commit}{path}#{line}",
    sitename = "GraphMatFun.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://matrixfunctions.github.io/GraphMatFun.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Manual" => ["graph.md", "generators.md", "optim.md", "code_gen.md"],
        "Index" => "docidx.md",
    ],
)

deploydocs(; repo = "github.com/matrixfunctions/GraphMatFun.jl", devbranch = "main")


text:
solve(i::MSetInstance, ::GreedyAlgorithm; kwargs...) = msets_greedy(i; kwargs...)

function msets_greedy(instance::MSetInstance)
  # Algorithm: sort the weights, take the m largest ones, this is the optimum m-set solution.
  # Implementation: no need for sorting, partialsortperm returns the largest items.
  items = collect(partialsortperm(instance.values, 1:instance.m, rev=true))
  return MSetSolution(instance, items)
end


text:
using Test
using ManifoldsBase
@testset "ManifoldsBase" begin
    bound = 0
    if VERSION >= v"1.1"
        num_ambiguities = length(Test.detect_ambiguities(ManifoldsBase))
        #num_ambiguities > 0 && @warn "The number of ambiguities in ManifoldsBase is $(num_ambiguities)."
        if VERSION >= v"1.8-DEV"
            @test num_ambiguities <= bound + 5
        elseif VERSION >= v"1.6-DEV"
            # At the time of writing there seem to be two ambiguities regarding `getindex`,
            # one with a method from SparseArrays and one from VSCode's JSON processing
            # that's automatically loaded when running code in VSCode.
            @test num_ambiguities <= bound + 3
        else
            @test num_ambiguities == bound + 1
        end
    end

    include("decorator_traits.jl")
    include("passthrough_decorator.jl")
    include("allocation.jl")
    include("numbers.jl")
    include("bases.jl")
    include("manifold_fallbacks.jl")
    include("empty_manifold.jl")
    include("errors.jl")
    include("default_manifold.jl")
    include("complex_manifold.jl")
    include("validation_manifold.jl")
    include("embedded_manifold.jl")
    include("power.jl")
    include("domain_errors.jl")
    include("vector_transport.jl")
end


text:
#///////////////////////////////////////
#// File Name: type_definition_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/15
#// Description: Test code for src/type_definition.jl
#///////////////////////////////////////

using Distributions
using LinearAlgebra
using RobotOS

@testset "Type Definition Test" begin

# WayPoint & Trajectory
p1 = [3., 4.];
p2 = [2., 1.];
wp1 = WayPoint2D(p1, Time(0.0));
wp2 = WayPoint2D(p2, Time(1.0));
traj = Trajectory2D([wp1, wp2]);
@test get_position(traj, Time(0.0)) == p1;
@test get_position(traj, Time(0.1)) ≈ [2.9, 3.7];
@test get_position(traj, Time(1.0)) == p2;

# Robot State
e1 = RobotState([p1; 0.; 0.]);
@test e1.x == [3., 4., 0., 0.];
@test get_position(e1) == [3., 4.];
@test get_velocity(e1) == [0., 0.];
@test e1 == RobotState([p1; 0.; 0;], Time(0.0));
@test e1 ≈ e1;

# World State
ap_dict = Dict("Pedestrian/1" => p2);
w = WorldState(e1, ap_dict);
@test w.t == Time(0.0);
@test isapprox(w, w);

end


text:
function solve()
    m = 1000000007
    s = parse(Int, readline())
    dp = zeros(Int, s)
    for i=3:s
        dp[i] = 1
        for j=3:i-3
            dp[i] = (dp[i] + dp[j]) % m
        end
    end
    dp[s]
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module BasisChange

export rotation_matrix

using LinearAlgebra

rotation_matrix(B, B′) = B*transpose(B′)

end


text:
# moving average functions
@testset "Moving Averages" begin
    Random.seed!(SEED)
    @testset "Array" begin
        x = cumsum(randn(N))
        X = cumsum(randn(N, 2), dims=1)
        tmp = sma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = mama(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
        @test sum(isnan.(tmp)) != N
        tmp = ema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = wma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = hma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = trima(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = mma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = tema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = dema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = swma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = kama(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = alma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = zlema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = vwma(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = vwap(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = hama(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
    end
    @testset "Temporal" begin
        x = TS(cumsum(randn(N)))
        X = TS(cumsum(randn(N, 2), dims=1))
        # moving average functions
        tmp = sma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = mama(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 2
        tmp = ema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = wma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = hma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = trima(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = mma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = tema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = dema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = swma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = kama(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = alma(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = zlema(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = vwma(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = vwap(X)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        tmp = hama(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
    end
end


text:
file_path = TEST_OUTPUTS

function test_reports(file_path::String; backend_pkg::String = "gr")
    if backend_pkg == "gr"
        backend = Plots.gr()
    elseif backend_pkg == "plotlyjs"
        backend = Plots.plotlyjs()
    else
        throw(error("$backend_pkg backend_pkg not supported"))
    end
    cleanup = true

    @testset "testing $backend_pkg report production" begin
        out_path = joinpath(file_path, backend_pkg * "_reports")
        !isdir(out_path) && mkpath(out_path)
        report_out_path = joinpath(out_path, "test_report.html")
        (results_uc, results_ed) = run_test_sim(TEST_RESULT_DIR)

        PG.report(
            results_uc,
            report_out_path,
            generic_template;
            doctype = "md2html",
            backend = backend,
        )
        @test isfile(report_out_path)

        @info("removing test files")
        cleanup && rm(out_path, recursive = true)
    end
end

try
    test_reports(file_path, backend_pkg = "gr")
    test_reports(file_path, backend_pkg = "plotlyjs")
finally
    nothing
end

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

module Sparso

package_path = joinpath(dirname(@__FILE__), "..", "deps")
push!(LOAD_PATH, package_path)

using CompilerTools
using CompilerTools.CFGs
using CompilerTools.LambdaHandling
using CompilerTools.OptFramework
using CompilerTools.LivenessAnalysis
using CompilerTools.AstWalker
using CompilerTools.Loops

# Options user can set
export SA_ENABLE, SA_VERBOSE, SA_USE_JULIA, SA_USE_MKL, SA_USE_SPMP,
       SA_REPLACE_CALLS, SA_USE_SPLITTING_PATTERNS, SA_REORDER, SA_CONTEXT,
       SA_CONTEXT_FUNC, SA_DISABLE_COLLECTIVE_STRUCTURE_PREDICTION

export @acc

# Function interface
export set_options, matrix_market_read

# Data structures user can extend
export FunctionDescription

# Explicitly specify matrix property
export  SA_CONST_VALUED, SA_CONST_SIZED, SA_MAXIMAL_STRUCTURED,
        SA_SYMM_VALUED, SA_SYMM_STRUCTURED,
        SA_STRUCTURE_ONLY, SA_HAS_DEDICATED_MEMORY

export SA_LOWER_OF, SA_UPPER_OF, SA_TRANSPOSE_OF

export set_matrix_property

include("globals.jl")
include("exceptions.jl")
include("function-descriptions.jl")
include("region-formation.jl")
include("call-replacement.jl")
include("code-transformation.jl")
include("constant.jl")
include("lib-interface.jl")
include("general-reordering.jl")
include("context.jl")
include("show.jl")
include("structure-analysis.jl")

end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
# incomplete

EVEN_DIRS = [[0, -1], [0, 1], [-1, -1]]
ODD_DIRS = [[0, -1], [0, 1], [1, 1]]

validcoords(S, r, p) = 1 <= r <= S && 1 <= p <= 2S - 1

getcoords(R, P) = map(x -> [R, P] .+ x, iseven(P) ? EVEN_DIRS : ODD_DIRS)

function simulate(S, RA, PA, RB, PB, construction, painted)
    coordfilter(r, p) = validcoords(S, r, p) && !((r, p) in construction) && !((r, p) in painted[])

    Achecked = filter(x -> coordfilter(x...), getcoords(RA, PA))
    Bchecked = filter(x -> coordfilter(x...), getcoords(RB, PB))

    length(Achecked) == length(Bchecked) == 0 && return 0

    maxscore = 0
    if length(Achecked) > 0 && length(Bchecked) > 0
        for (ra, pa) in Achecked
            for (rb, pb) in Bchecked
                push!(painted[], (ra, pa))
                push!(painted[], (rb, pb))
                maxscore = max(maxscore, simulate(S, ra, pa, rb, pb, construction, painted))
                delete!(painted[], (ra, pa))
                delete!(painted[], (rb, pb))
            end
        end
    elseif length(Achecked) > 0
        for (ra, pa) in Achecked
            push!(painted[], (ra, pa))
            maxscore = max(maxscore, simulate(S, ra, pa, RB, PB, construction, painted)) + 1
            delete!(painted[], (ra, pa))
        end
    elseif length(Bchecked) > 0
        for (rb, pb) in Bchecked
            push!(painted[], (rb, pb))
            maxscore = max(maxscore, simulate(S, RA, PA, rb, pb, construction, painted)) - 1
            delete!(painted[], (rb, pb))
        end
    end

    maxscore
end

function solve(case, S, RA, PA, RB, PB, C, construction)
    score = simulate(S, RA, PA, RB, PB, construction, Ref(Set{Tuple{Int,Int}}()))
    println("Case $case: $score")
end

function main()
    T = parse(Int, readline(stdin))

    for case in 1:T
        S, RA, PA, RB, PB, C = parse.(Int, stdin |> readline |> split)

        construction = Set{Tuple{Int,Int}}()
        for _ in 1:C
            RI, PI = parse.(Int, stdin |> readline |> split)
            push!(construction, (RI, PI))
        end

        solve(case, S, RA, PA, RB, PB, C, construction)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
include("simple.jl")
include("labels.jl")
include("ticks.jl")
include("scale.jl")


text:
module TestGraphvizGraphs
using Test
import JSON

using Catlab.Graphs, Catlab.Graphics.GraphvizGraphs
import Catlab.Graphics: Graphviz

function stmts(graph::Graphviz.Graph, type::Type)
  [ stmt for stmt in graph.stmts if stmt isa type ]
end
function stmts(graph::Graphviz.Graph, type::Type, attr::Symbol)
  [ stmt.attrs[attr] for stmt in graph.stmts
    if stmt isa type && haskey(stmt.attrs, attr) ]
end

# Property graphs
#################

# Parsing
#--------

data_path(name::String) = joinpath(@__DIR__, "data", name)

# Undirected simple graph.
doc = open(JSON.parse, data_path("graphviz_graph.json"), "r")
parsed = parse_graphviz(doc)
@test parsed isa SymmetricPropertyGraph
@test nv(parsed) == 10
@test ne(parsed) ÷ 2 == 13

# Directed simple graph
doc = open(JSON.parse, data_path("graphviz_digraph.json"), "r")
parsed = parse_graphviz(doc)
@test parsed isa PropertyGraph
@test nv(parsed) == 5
@test src(parsed) == [1,1,1,2,3,4]
@test tgt(parsed) == [2,3,4,5,5,5]
@test get_vprop.([parsed], 1:nv(parsed), :name) == ["A", "B", "C", "D", "F"]
@test get_vprop(parsed, 1, :position) == [99, 162]
@test get_vprop(parsed, 1, :size) == [54, 36]
@test get_eprop(parsed, 1, :spline) isa Vector

# Visualizion
#------------

# Symmetric property graph.
g = SymmetricPropertyGraph{String}()
add_vertex!(g, label="v")
add_vertices!(g, 2)
add_edge!(g, 1, 2, xlabel="e")
gv = to_graphviz(g)
@test !gv.directed
@test length(stmts(gv, Graphviz.Node)) == 3
@test length(stmts(gv, Graphviz.Edge)) == 1
@test stmts(gv, Graphviz.Node, :label) == ["v"]
@test stmts(gv, Graphviz.Edge, :xlabel) == ["e"]

# Property graph.
g = PropertyGraph{String}()
add_vertices!(g, 3); add_edge!(g, 1, 2); add_edge!(g, 2, 3)
gv = to_graphviz(g)
@test gv.directed
@test length(stmts(gv, Graphviz.Node)) == 3
@test length(stmts(gv, Graphviz.Edge)) == 2

# Property graph with multiple edges.
g = PropertyGraph{String}()
add_vertices!(g, 2)
add_edge!(g, 1, 2, label="e1")
add_edge!(g, 1, 2, label="e2")
gv = to_graphviz(g)
@test gv.directed
@test length(stmts(gv, Graphviz.Node)) == 2
@test length(stmts(gv, Graphviz.Edge)) == 2
@test stmts(gv, Graphviz.Edge, :label) == ["e1", "e2"]

# Graphs
########

g = path_graph(Graph, 3)
gv = to_graphviz(g)
@test gv.directed
@test stmts(gv, Graphviz.Node, :label) == fill("", 3)

gv = to_graphviz(g, node_labels=true, edge_labels=true)
@test stmts(gv, Graphviz.Node, :label) == ["1", "2", "3"]
@test stmts(gv, Graphviz.Edge, :label) == ["1", "2"]

# Symmetric graphs
##################

g = path_graph(SymmetricGraph, 3)
gv = to_graphviz(g, edge_labels=true)
@test !gv.directed
@test stmts(gv, Graphviz.Node, :label) == fill("", 3)
@test stmts(gv, Graphviz.Edge, :label) == ["(1,3)", "(2,4)"]

# Reflexive graphs
##################

g = path_graph(ReflexiveGraph, 3)
gv = to_graphviz(g, show_reflexive=false)
@test gv.directed
@test length(stmts(gv, Graphviz.Node)) == 3
@test length(stmts(gv, Graphviz.Edge)) == 2

gv = to_graphviz(g, edge_labels=true, show_reflexive=true)
@test length(stmts(gv, Graphviz.Edge)) == 5
@test stmts(gv, Graphviz.Edge, :label) == ["1", "2", "3", "4", "5"]
@test stmts(gv, Graphviz.Edge, :style) == fill("dashed", 3)

# Symmetric reflexive graphs
############################

g = path_graph(SymmetricReflexiveGraph, 3)
gv = to_graphviz(g, show_reflexive=false)
@test !gv.directed
@test length(stmts(gv, Graphviz.Node)) == 3
@test length(stmts(gv, Graphviz.Edge)) == 2

gv = to_graphviz(g, edge_labels=true, show_reflexive=true)
@test length(stmts(gv, Graphviz.Edge)) == 5
@test stmts(gv, Graphviz.Edge, :label) == ["1", "2", "3", "(4,6)", "(5,7)"]
@test stmts(gv, Graphviz.Edge, :style) == fill("dashed", 3)

# Half-edge graphs
##################

g = path_graph(HalfEdgeGraph, 3)
gv = to_graphviz(g, edge_labels=true)
@test !gv.directed
@test stmts(gv, Graphviz.Node, :label) == fill("", 3)
@test stmts(gv, Graphviz.Edge, :label) == ["(1,3)", "(2,4)"]

g = HalfEdgeGraph(2)
add_edge!(g, 1, 2)
add_dangling_edges!(g, [1,2])
gv = to_graphviz(g, node_labels=true, edge_labels=true)
@test stmts(gv, Graphviz.Node, :label) == ["1", "2", "", ""]
@test stmts(gv, Graphviz.Edge, :label) == ["(1,2)", "3", "4"]

end


text:
using CurrenciesBase
using FixedPointDecimals

using Test
using Dates

# Get currencies for tests
@usingcurrencies USD, CAD, EUR, GBP, JPY, AUD, INR
@usingcurrencies CNY  # test one-currency version

# Data tests
include("data.jl")

# Basic functionality tests
include("monetary.jl")
# include("basket.jl")
# include("mixed.jl")

# # Currencies
# include("data-access.jl")

# # Display tests
# include("display.jl")

# # Custom currencies
# include("custom.jl")

# # Computations tests
# include("valuation.jl")

# README & Doc examples
include("examples.jl")


text:
using AxisKeys
using CSV
using DataFrames
using StatsBase
using NamedDims
using Turing

@testset "LooCompare" begin
    using Random
    Random.seed!(129111)

    df = CSV.read("data/WaffleDivorce.csv", DataFrame)
    df.D = zscore(df.Divorce)
    df.M = zscore(df.Marriage)
    df.A = zscore(df.MedianAgeMarriage)
    data = (D=df.D, A=df.A)

    function lin(a, b, c, x...)
        result = @. a + b * c
        for i in 1:2:length(x)
            @. result += x[i] * x[i + 1]
        end
        return result
    end

    @model function m5_1t(A, D)
        a ~ Normal(0, 0.2)
        bA ~ Normal(0, 0.5)
        σ ~ Exponential(1)
        for i in eachindex(D)
            μ = lin(a, A[i], bA)
            D[i] ~ Normal(μ, σ)
        end
    end

    chn5_1t = sample(m5_1t(df.A, df.D), NUTS(1000, 0.9), MCMCThreads(), 1000, 12)

    @model function m5_2t(M, D)
        a ~ Normal(0, 0.2)
        bM ~ Normal(0, 0.5)
        σ ~ Exponential(1)
        for i in eachindex(D)
            μ = lin(a, M[i], bM)
            D[i] ~ Normal(μ, σ)
        end
    end

    chn5_2t = sample(m5_2t(df.M, df.D), NUTS(1000, 0.9), MCMCThreads(), 1000, 12)

    @model function m5_3t(A, M, D)
        a ~ Normal(0, 0.2)
        bA ~ Normal(0, 0.5)
        bM ~ Normal(0, 0.5)
        σ ~ Exponential(1)
        for i in eachindex(D)
            μ = a + M[i] * bM + A[i] * bA
            D[i] ~ Normal(μ, σ)
        end
    end

    chn5_3t = sample(m5_3t(df.A, df.M, df.D), NUTS(1000, 0.9), MCMCThreads(), 1000, 12)
    pw_lls5_1t = pointwise_log_likelihoods(m5_1t(df.A, df.D), chn5_1t)
    psis_loo_output5_1t = psis_loo(m5_1t(df.A, df.D), chn5_1t)
    psis_loo_output5_1t |> display

    pw_lls5_2t = pointwise_log_likelihoods(m5_2t(df.M, df.D), chn5_2t)
    psis_loo_output5_2t = psis_loo(m5_2t(df.M, df.D), chn5_2t)
    psis_loo_output5_2t |> display

    pw_lls5_3t = pointwise_log_likelihoods(m5_3t(df.A, df.M, df.D), chn5_3t)
    psis_loo_output5_3t = psis_loo(m5_3t(df.A, df.M, df.D), chn5_3t)
    psis_loo_output5_3t |> display

    n_tuple = (
        m5_1t=psis_loo_output5_1t,
        m5_2t=psis_loo_output5_2t,
        m5_3t=psis_loo_output5_3t,
    )
    comps = loo_compare(n_tuple)
    comps |> display

    tol = 2 * sqrt(sum([x.mcse ^ 2 for x in n_tuple]))
    @test comps.estimates(:m5_1t, :cv_elpd) ≈ 0.00 atol = .001
    @test comps.std_err[:m5_1t] ≈ 0.00 atol = .001
    @test comps.estimates(:m5_1t, :weight) ≈ .7075 rtol = tol
    @test comps.estimates(:m5_2t, :cv_elpd) ≈ -6.62779 atol = tol
    @test comps.std_err[:m5_2t] ≈ 4.71968 rtol = tol
    @test comps.estimates(:m5_2t, :weight) ≈ 0.00 atol = .001
    @test comps.estimates(:m5_3t, :cv_elpd) ≈ -0.8865456492718806 atol = tol
    @test comps.std_err[:m5_3t] ≈ 0.3638491385753804 rtol = tol
    @test comps.estimates(:m5_3t, :weight) ≈ 0.2915500227946632 atol = tol
    @test sum(comps.estimates(:, :weight, :)) ≈ 1
    total = NamedDims.unname(sum(comps.pointwise(:, :cv_elpd, :); dims=:data))
    @test reshape(total, 3) ≈ comps.estimates(:, :cv_elpd) atol=.001

    gmpd_1 = getproperty.(values(n_tuple), :gmpd)
    gmpd_1 = sort([i for i in gmpd_1])
    comp_gmpd = sort([i for i in values(comps.gmpd)])
    @test comp_gmpd ≈ gmpd_1 atol = .001

end


text:
module TestCopy

using Test
using Transducers

@testset "copy [base types]" begin
    @test copy(Map(x -> x^2), Vector, 1:3) == (1:3) .^ 2
    @test copy(Map(x -> x^2), [1:3;]) == (1:3) .^ 2
    @test copy(Map(x -> x^2), Set, 1:3) == Set((1:3) .^ 2)
    @test copy(Map(x -> x => x^2), Dict, 2:2) == Dict(2 => 4)
end

@testset "Set" begin
    @test Set(eduction(y for x in 1:1000 for y in 1:x)) == Set(1:1000)
end

@testset "Dict" begin
    @test Dict(eduction(x => y for x in 1:1000 for y in 1:x)) == Dict(1:1000 .=> 1:1000)
end

end  # module


text:
struct ZIGammaExpert{T<:Real} <: ZIContinuousExpert
    p::T
    k::T
    θ::T
    ZIGammaExpert{T}(p::T, k::T, θ::T) where {T<:Real} = new{T}(p, k, θ)
end

function ZIGammaExpert(p::T, k::T, θ::T; check_args=true) where {T <: Real}
    check_args && @check_args(ZIGammaExpert, 0 <= p <= 1 && k >= zero(k) && θ > zero(θ))
    return ZIGammaExpert{T}(p, k, θ)
end

#### Outer constructors
ZIGammaExpert(p::Real, k::Real, θ::Real) = ZIGammaExpert(promote(p, k, θ)...)
ZIGammaExpert(p::Integer, k::Integer, θ::Integer) = ZIGammaExpert(float(p), float(k), float(θ))
ZIGammaExpert() = ZIGammaExpert(0.50, 1.0, 1.0)

## Conversion
function convert(::Type{ZIGammaExpert{T}}, p::S, k::S, θ::S) where {T <: Real, S <: Real}
    ZIGammaExpert(T(p), T(k), T(θ))
end
function convert(::Type{ZIGammaExpert{T}}, d::ZIGammaExpert{S}) where {T <: Real, S <: Real}
    ZIGammaExpert(T(d.p), T(d.k), T(d.θ), check_args=false)
end
copy(d::ZIGammaExpert) = ZIGammaExpert(d.p, d.k, d.θ, check_args=false)

## Loglikelihood of Expoert
logpdf(d::ZIGammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? -Inf : Distributions.logpdf.(Distributions.Gamma(d.k, d.θ), x...)
pdf(d::ZIGammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? 0.0 : Distributions.pdf.(Distributions.Gamma(d.k, d.θ), x...)
logcdf(d::ZIGammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? -Inf : Distributions.logcdf.(Distributions.Gamma(d.k, d.θ), x...)
cdf(d::ZIGammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? 0.0 : Distributions.cdf.(Distributions.Gamma(d.k, d.θ), x...)

## expert_ll, etc
expert_ll_exact(d::ZIGammaExpert, x::Real) = (x == 0.) ? log(p_zero(d)) :  log(1-p_zero(d)) + LRMoE.logpdf(d, x)
function expert_ll(d::ZIGammaExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll_pos = LRMoE.expert_ll(LRMoE.GammaExpert(d.k, d.θ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_ll = (yl == 0.) ? log.(p0 + (1-p0)*exp.(expert_ll_pos)) : log.(0.0 + (1-p0)*exp.(expert_ll_pos))
    expert_ll = (tu == 0.) ? log.(p0) : expert_ll
    return expert_ll
end
function expert_tn(d::ZIGammaExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_pos = LRMoE.expert_tn(LRMoE.GammaExpert(d.k, d.θ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_tn = (tl == 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_pos))
    expert_tn = (tu == 0.) ? log.(p0) : expert_tn
    return expert_tn
end
function expert_tn_bar(d::ZIGammaExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar_pos = LRMoE.expert_tn_bar(LRMoE.GammaExpert(d.k, d.θ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_tn_bar = (tl > 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_bar_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_bar_pos))
    return expert_tn_bar
end

exposurize_expert(d::ZIGammaExpert; exposure = 1) = d

## Parameters
params(d::ZIGammaExpert) = (d.p, d.k, d.θ)
p_zero(d::ZIGammaExpert) = d.p
function params_init(y, d::ZIGammaExpert)
    p_init = sum(y .== 0.0) / sum(y .>= 0.0)
    pos_idx = (y .> 0.0)
    μ, σ2 = mean(y[pos_idx]), var(y[pos_idx])
    θ_init = σ2/μ
    k_init = μ/θ_init
    if isnan(θ_init) || isnan(k_init)
        return ZIGammaExpert()
    else
        return ZIGammaExpert(p_init, k_init, θ_init)
    end
end

## KS stats for parameter initialization
function ks_distance(y, d::ZIGammaExpert)
    p_zero = sum(y .== 0.0) / sum(y .>= 0.0)
    return max(abs(p_zero-d.p), (1-d.p)*HypothesisTests.ksstats(y[y .> 0.0], Distributions.Gamma(d.k, d.θ))[2])
end

## Simululation
sim_expert(d::ZIGammaExpert) = (1 .- Distributions.rand(Distributions.Bernoulli(d.p), 1)[1]) .* Distributions.rand(Distributions.Gamma(d.k, d.θ), 1)[1]

## penalty
penalty_init(d::ZIGammaExpert) = [2.0 10.0 2.0 10.0]
no_penalty_init(d::ZIGammaExpert) = [1.0 Inf 1.0 Inf]
penalize(d::ZIGammaExpert, p) = (p[1]-1)*log(d.k) - d.k/p[2] + (p[3]-1)*log(d.θ) - d.θ/p[4]

## statistics
mean(d::ZIGammaExpert) = (1-d.p)*mean(Distributions.Gamma(d.k, d.θ))
var(d::ZIGammaExpert) = (1-d.p)*var(Distributions.Gamma(d.k, d.θ)) + d.p*(1-d.p)*(mean(Distributions.Gamma(d.k, d.θ)))^2
quantile(d::ZIGammaExpert, p) = p <= d.p ? 0.0 : quantile(Distributions.Gamma(d.k, d.θ), p-d.p)
lev(d::ZIGammaExpert, u) = (1-d.p)*lev(GammaExpert(d.k, d.θ), u)
excess(d::ZIGammaExpert, u) = mean(d) - lev(d, u)

## EM: M-Step
function EM_M_expert(d::ZIGammaExpert,
                     tl, yl, yu, tu,
                     exposure,
                    #  expert_ll_pos,
                    #  expert_tn_pos,
                    #  expert_tn_bar_pos,
                     z_e_obs, z_e_lat, k_e;
                     penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Update zero probability
    expert_ll_pos = expert_ll.(LRMoE.GammaExpert(d.k, d.θ), tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(LRMoE.GammaExpert(d.k, d.θ), tl, yl, yu, tu)

    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(yl, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    z_zero_e_lat = z_e_lat .* EM_E_z_zero_lat(tl, p_old, expert_tn_bar_pos)
    z_pos_e_lat = z_e_lat .- z_zero_e_lat
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, z_zero_e_lat, z_pos_e_lat, k_e)

    # Update parameters: call its positive part
    tmp_exp = GammaExpert(d.k, d.θ)
    tmp_update = EM_M_expert(tmp_exp,
                            tl, yl, yu, tu,
                            exposure,
                            # expert_ll_pos,
                            # expert_tn_pos,
                            # expert_tn_bar_pos,
                            # z_e_obs, z_e_lat, k_e,
                            z_pos_e_obs, z_pos_e_lat, k_e,
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIGammaExpert(p_new, tmp_update.k, tmp_update.θ)
end

## EM: M-Step, exact observations
function EM_M_expert_exact(d::ZIGammaExpert,
                     ye, exposure,
                     z_e_obs;
                     penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Update zero probability
    expert_ll_pos = expert_ll_exact.(LRMoE.GammaExpert(d.k, d.θ), ye)

    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(ye, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, 0.0, 0.0, 0.0)

    # Update parameters: call its positive part
    tmp_exp = GammaExpert(d.k, d.θ)
    tmp_update = EM_M_expert_exact(tmp_exp,
                            ye, exposure,
                            z_pos_e_obs;
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIGammaExpert(p_new, tmp_update.k, tmp_update.θ)
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
import Base.+, Base.-, Base./, Base.*, Base.^
import SchumakerSpline.evaluate
function evaluate(f::Undefined_Function, point::Union{Real,Date,DateTime,DatePeriod})
    missing
end
function (s::Undefined_Function)(x::Union{Real,Date,DateTime,DatePeriod})
    return missing
end

function derivative(f::Undefined_Function)
    return f
end
function indefinite_integral(f::Undefined_Function)
    return f
end

function +(f::Undefined_Function, number::Real)
    return f
end
function -(f::Undefined_Function, number::Real)
    return f
end
function *(f::Undefined_Function, number::Real)
    return f
end
function /(f::Undefined_Function, number::Real)
    return f
end
function ^(f::Undefined_Function, number::Integer)
    return f
end

function +(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function +(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function +(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function +(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function +(f1::PE_Function, f2::Undefined_Function)
    return +(f2,f1)
end
function +(f1::Sum_Of_Functions, f2::Undefined_Function)
    return +(f2,f1)
end
function +(f1::Piecewise_Function, f2::Undefined_Function)
    return +(f2,f1)
end

function -(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function -(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function -(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function -(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function -(f1::PE_Function, f2::Undefined_Function)
    return -(f2,f1)
end
function -(f1::Sum_Of_Functions, f2::Undefined_Function)
    return -(f2,f1)
end
function -(f1::Piecewise_Function, f2::Undefined_Function)
    return -(f2,f1)
end

function *(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function *(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function *(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function *(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function *(f1::PE_Function, f2::Undefined_Function)
    return *(f2,f1)
end
function *(f1::Sum_Of_Functions, f2::Undefined_Function)
    return *(f2,f1)
end
function *(f1::Piecewise_Function, f2::Undefined_Function)
    return *(f2,f1)
end

function /(f1::Undefined_Function, f2::Undefined_Function)
    return f1
end
function /(f1::Undefined_Function, f2::PE_Function)
    return f1
end
function /(f1::Undefined_Function, f2::Sum_Of_Functions)
    return f1
end
function /(f1::Undefined_Function, f2::Piecewise_Function)
    return f1
end

function /(f1::PE_Function, f2::Undefined_Function)
    return /(f2,f1)
end
function /(f1::Sum_Of_Functions, f2::Undefined_Function)
    return /(f2,f1)
end
function /(f1::Piecewise_Function, f2::Undefined_Function)
    return /(f2,f1)
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ 73094af4-3595-11eb-269d-e3bf33913758
using Pkg, DrWatson

# ╔═╡ 73323464-3595-11eb-3be0-b9a9f1ea2f8f
begin
    @quickactivate "StatisticalRethinkingTuring"
    using Turing
    using StatisticalRethinking
end

# ╔═╡ 5550278a-3595-11eb-0e3e-f12ab620d448
md"## Coinflip.jl"

# ╔═╡ 7332c5fa-3595-11eb-288b-dd932e0953bf
struct MyDist{T} <: ContinuousUnivariateDistribution
    θ::T
end

# ╔═╡ 734cff6a-359a-11eb-2cde-4756a76d07a7
md"##### I suspect that the for loop in coinflip1 performs a separate evaluation on each iteration with overhead. Perhaps there is someway to improve performance on the Turing side. As a workaround, you can create a new distribution type and place the for loop inside the logpdf function. I use this in cases where
1. I am defining a custom distribution,
2. I am using a non-standard data type,
3. or a for loop is a more natural way to define the model."

# ╔═╡ cf9e35ae-359a-11eb-02ec-0550744aed96
md"##### Below some examples for comparison."

# ╔═╡ 73435406-3595-11eb-3e69-9946981f886d
begin
	import Distributions: logpdf
	function logpdf(d::MyDist, data::Array{Bool,1})
		LL = 0.0
		for i in 1:length(data)
			LL += logpdf(Bernoulli(d.θ), data[i])
		end
		return LL
	end
end

# ╔═╡ 73450f44-3595-11eb-2707-65c5f3ff8451
begin
	import Distributions: loglikelihood
	loglikelihood(d::MyDist, data::Array{Bool,1}) = logpdf(d, data)
end

# ╔═╡ 7350c640-3595-11eb-1f8c-41159f4490cf
begin
    p_true = 0.5
    n = 1000
    Random.seed!(12)
    data = rand(Bernoulli(p_true), n)
end;

# ╔═╡ 73515718-3595-11eb-0c71-b3c8fae4e95b
begin
    n_samples = 2000
    delta = .65
    n_adapt = 1000
    config = NUTS(n_adapt, delta)
end;

# ╔═╡ 735effce-3595-11eb-06c4-9b8202849af3
@model function coinflip1(y)
    p ~ Beta(1.,1.)
    for i in 1:length(y)
        y[i] ~ Bernoulli(p)
    end
end

# ╔═╡ 735fa414-3595-11eb-3c38-8d1ba76fa28e
md"##### Is coinflip2 code correct (no documentation, but it works)?"

# ╔═╡ 736e5b4c-3595-11eb-1ff2-5b2ca39213ac
@model function coinflip2(y)
    p ~ Beta(1.,1.)
    y ~ Bernoulli(p)
end

# ╔═╡ 737a61da-3595-11eb-3f65-ab61047f1fbc
md"##### This model is from 'DynamicPPL: Stan-like Speed for Dynamic Probabilistic Models'"

# ╔═╡ 737d45c6-3595-11eb-0fec-3b044ca5e64c
@model function coinflip3(y)
    p ~ Beta(1.,1.)
    y .~ Bernoulli(p)
end

# ╔═╡ 7b19ca74-3599-11eb-1370-1f773dc7ae9d
md"##### This code uses above defined definition."

# ╔═╡ 738d7f04-3595-11eb-1b38-0315d5ae6083
@model function coinflip4(y)
    p ~ Beta(1.,1.)
    y ~ MyDist(p)
end

# ╔═╡ 4db69812-3598-11eb-31ea-a152bc58aca9
md"#### Note: Output of `@time ...` to Terminal (or in  right lower cell corner)."

# ╔═╡ 73944a32-3595-11eb-21a4-4f76a53efac2
@time chain1 = sample(coinflip1(data), config, n_samples);

# ╔═╡ 49bdbc14-359c-11eb-3e77-fb604ac0bf4f
@time chain2 = sample(coinflip2(data), config, n_samples);

# ╔═╡ 49bde4d2-359c-11eb-3440-0b244f8f784d
@time chain3 = sample(coinflip3(data), config, n_samples);

# ╔═╡ 49be84fa-359c-11eb-0f24-37ff390dc0b4
@time chain4 = sample(coinflip4(data), config, n_samples);

# ╔═╡ 90e7082e-3598-11eb-2ecf-f78e174b3dc8
begin
	chns = chainscat(chain1, chain2, chain3, chain4)
	CHNS(chns)
end

# ╔═╡ 73bbc712-3595-11eb-1eab-5fb16c871e73
md"## End of coinflip.jl"

# ╔═╡ Cell order:
# ╟─5550278a-3595-11eb-0e3e-f12ab620d448
# ╠═73094af4-3595-11eb-269d-e3bf33913758
# ╠═73323464-3595-11eb-3be0-b9a9f1ea2f8f
# ╠═7332c5fa-3595-11eb-288b-dd932e0953bf
# ╟─734cff6a-359a-11eb-2cde-4756a76d07a7
# ╟─cf9e35ae-359a-11eb-02ec-0550744aed96
# ╠═73435406-3595-11eb-3e69-9946981f886d
# ╠═73450f44-3595-11eb-2707-65c5f3ff8451
# ╠═7350c640-3595-11eb-1f8c-41159f4490cf
# ╠═73515718-3595-11eb-0c71-b3c8fae4e95b
# ╠═735effce-3595-11eb-06c4-9b8202849af3
# ╟─735fa414-3595-11eb-3c38-8d1ba76fa28e
# ╠═736e5b4c-3595-11eb-1ff2-5b2ca39213ac
# ╟─737a61da-3595-11eb-3f65-ab61047f1fbc
# ╠═737d45c6-3595-11eb-0fec-3b044ca5e64c
# ╟─7b19ca74-3599-11eb-1370-1f773dc7ae9d
# ╠═738d7f04-3595-11eb-1b38-0315d5ae6083
# ╟─4db69812-3598-11eb-31ea-a152bc58aca9
# ╠═73944a32-3595-11eb-21a4-4f76a53efac2
# ╠═49bdbc14-359c-11eb-3e77-fb604ac0bf4f
# ╠═49bde4d2-359c-11eb-3440-0b244f8f784d
# ╠═49be84fa-359c-11eb-0f24-37ff390dc0b4
# ╠═90e7082e-3598-11eb-2ecf-f78e174b3dc8
# ╟─73bbc712-3595-11eb-1eab-5fb16c871e73


text:
occa_github = "https://github.com/tcew/OCCA2.git";

occapkgdir = Pkg.dir() * "/OCCA";

if !isfile(occapkgdir * "/src/occabuiltwith.jl")
    f=open(occapkgdir * "/src/occabuiltwith.jl","w");
    write(f,"OCCA_USE_OPENMP   = false;\n");
    write(f,"OCCA_USE_PTHREADS = false;\n");
    write(f,"OCCA_USE_OPENCL   = false;\n");
    write(f,"OCCA_USE_CUDA     = false;\n");
    close(f);
end

f = open(occapkgdir * "/src/occapaths.jl","w");
write(f, "const libocca=\"" * occapkgdir * "/deps/OCCA2/lib/libocca.so\"");
close(f);

using OCCA;

currentdir = pwd();
cd(occapkgdir);

occadir = occapkgdir * "/deps/OCCA2";

# If OCCA2 doesn't exist, download source from git.
if !isdir(occadir);
    run(`git submodule add $occa_github deps/OCCA2`);
else
    run(`git submodule foreach git pull`);
end

cd(currentdir);

# Set necessary environment variables
ENV["OCCA_DIR"] = occadir;

if OCCA.USE_OPENMP
    ENV["OCCA_OPENMP_ENABLED"] = 1;
end
if OCCA.USE_PTHREADS
    ENV["OCCA_PTHREADS_ENABLED"] = 1;
end
if OCCA.USE_CUDA
    ENV["OCCA_CUDA_ENABLED"] = 1;
end
if OCCA.USE_OPENCL
    ENV["OCCA_OPENCL_ENABLED"] = 1;
end

# Run the main build command.
run(`make -f $occadir/makefile clean`);
run(`make -j -f $occadir/makefile COMPILING_WITH_JULIA=1`);


text:
using Test
using CommonSubexpressions

@testset "basic usage" begin
    f_count = 0
    function f(x)
        f_count += 1
        1
    end

    g_count = 0
    function g(x)
        g_count += 1
        2
    end

    @test f_count == 0
    @test g_count == 0

    x = 1
    @cse [f(g(x)) for i in 1:10]
    @test f_count == 1
    @test g_count == 1

    f_count = 0
    g_count = 0
    [f(g(x)) for i in 1:10]
    @test f_count == 10
    @test g_count == 10

    h_count = 0
    function h(x, y)
        h_count += 1
        3
    end

    @cse function foo(x)
        [h(h(f(g(x)), g(x)), i) == i for i in 1:5]
    end
    f_count = 0
    g_count = 0
    h_count = 0
    @test foo(1) == [false, false, true, false, false]
    @test f_count == 1
    @test g_count == 1
    @test h_count == 6
end

@testset "function definition" begin
    f_count = 0
    function f(x)
        f_count += 1
        1
    end

    @cse function foo(x)
        [f(x) == i for i in 1:5]
    end

    @test foo(1) == [true, false, false, false, false]
    @test f_count == 1
end

@testset "algebra" begin
    H = [1 2; 3 4]
    W = [2 3; 4 5]
    G = [4 5; 6 7]

    @test (@cse inv(H) * G + W * (inv(H) + W)) == begin
        invH = inv(H)
        invH * G + W * (inv(H) + W)
    end
end

@testset "dict" begin
    @test (@cse begin
        x = Dict("foo" => sin(pi), "bar" => 2 + 2)
        x["bar"] = 100
        x["foo"], x["bar"]
    end) == begin
        x = Dict("foo" => sin(pi), "bar" => 2 + 2)
        x["bar"] = 100
        x["foo"], x["bar"]
    end
end

@testset "int" begin
    @test (@cse 1) == 1
end

@testset "nested usage" begin
    f_count = 0
    function f(x)
        f_count += 1
        x
    end

    g_count = 0
    function g(x)
        g_count += 1
        2x
    end

    @test f_count == 0
    @test g_count == 0

    result = @cse for i in 1:5
        x = f(2)
        @cse for j in 1:5
            y = g(i)
        end
    end
    @test f_count == 1
    @test g_count == 5

    @test result == for i in 1:5
        x = f(2)
        for j in 1:5
            y = g(i)
        end
    end
end

@testset "cse function" begin
    expr = cse(:((x + 1) * (x + 1)))
    @test expr.head == :block
    @test length(expr.args) == 3
    @test expr.args[1].head == :(=)
    @test expr.args[2].head == :(=)
    @test expr.args[1].args[2].head == :call
    @test expr.args[1].args[2].args == [:(+), :x, 1]
    @test expr.args[2].args[2].head == :call
    @test expr.args[2].args[2].args[1] == :(*)
end

@testset "inplace" begin

    x = 1
    @cse begin
        y = 2 + x
        x += y * x
    end
    @test y == 3
    @test x == 4
end


text:
#
# Represents a tensor order that could be set to any order.
#

struct EmptyOrder end

#
# Represents a number that can be set to any type.
#

struct EmptyNumber <: Number end

# This is a backup definition to make:
# A = ITensor(i, j)
# complex!(A)
# work. It acts as if the "default" type is `Float64`
complex(::Type{EmptyNumber}) = ComplexF64

function similartype(::Type{StoreT}, ::Type{ElT}) where {StoreT<:Dense{EmptyNumber},ElT}
  return Dense{ElT,similartype(datatype(StoreT), ElT)}
end

function similartype(
  ::Type{StoreT}, ::Type{ElT}
) where {StoreT<:BlockSparse{EmptyNumber},ElT}
  return BlockSparse{ElT,similartype(datatype(StoreT), ElT),ndims(StoreT)}
end

#
# Empty storage
#

struct EmptyStorage{ElT,StoreT<:TensorStorage} <: TensorStorage{ElT} end

# Get the EmptyStorage version of the TensorStorage
function emptytype(::Type{StoreT}) where {StoreT}
  return EmptyStorage{eltype(StoreT),StoreT}
end

empty(::Type{StoreT}) where {StoreT} = emptytype(StoreT)()

# Defaults to Dense
function EmptyStorage(::Type{ElT}) where {ElT}
  return emptytype(Dense{ElT,Vector{ElT}})()
end

EmptyStorage() = EmptyStorage(Float64)

similar(S::EmptyStorage) = S
similar(S::EmptyStorage, ::Type{ElT}) where {ElT} = empty(similartype(fulltype(S), ElT))

copy(S::EmptyStorage) = S

isempty(::EmptyStorage) = true

nnzblocks(::EmptyStorage) = 0

nnz(::EmptyStorage) = 0

function Base.real(::Type{<:EmptyStorage{ElT,StoreT}}) where {ElT,StoreT}
  return EmptyStorage{real(ElT),real(StoreT)}
end

Base.real(S::EmptyStorage) = real(typeof(S))()

function complex(::Type{<:EmptyStorage{ElT,StoreT}}) where {ElT,StoreT}
  return EmptyStorage{complex(ElT),complex(StoreT)}
end

complex(S::EmptyStorage) = complex(typeof(S))()

#size(::EmptyStorage) = 0

function show(io::IO, mime::MIME"text/plain", S::EmptyStorage)
  return println(io, typeof(S))
end

#
# EmptyTensor (Tensor using EmptyStorage storage)
#

const EmptyTensor{ElT,N,StoreT,IndsT} =
  Tensor{ElT,N,StoreT,IndsT} where {StoreT<:EmptyStorage}

function emptytype(::Type{TensorT}) where {TensorT<:Tensor}
  return Tensor{
    eltype(TensorT),ndims(TensorT),emptytype(storagetype(TensorT)),indstype(TensorT)
  }
end

# XXX TODO: add bounds checking
getindex(T::EmptyTensor, I::Integer...) = zero(eltype(T))
getindex(T::EmptyTensor{EmptyNumber}, I::Integer...) = zero(Float64)

similar(T::EmptyTensor, inds::Tuple) = setinds(T, inds)
function similar(T::EmptyTensor, ::Type{ElT}) where {ElT<:Number}
  return tensor(similar(storage(T), ElT), inds(T))
end

function randn!!(T::EmptyTensor)
  Tf = similar(fulltype(T), inds(T))
  randn!(Tf)
  return Tf
end

# Default to Float64
function randn!!(T::EmptyTensor{EmptyNumber})
  return randn!!(similar(T, Float64))
end

function _fill!!(::Type{ElT}, T::EmptyTensor, α::Number) where {ElT}
  Tf = similar(fulltype(T), ElT, inds(T))
  fill!(Tf, α)
  return Tf
end

fill!!(T::EmptyTensor, α::Number) = _fill!!(eltype(T), T, α)

# Determine the element type from the number you are filling with
fill!!(T::EmptyTensor{EmptyNumber}, α::Number) = _fill!!(eltype(α), T, α)

isempty(::EmptyTensor) = true

function EmptyTensor(::Type{ElT}, inds) where {ElT<:Number}
  return tensor(EmptyStorage(ElT), inds)
end

function EmptyTensor(::Type{StoreT}, inds) where {StoreT<:TensorStorage}
  return tensor(empty(StoreT), inds)
end

function EmptyBlockSparseTensor(::Type{ElT}, inds) where {ElT<:Number}
  StoreT = BlockSparse{ElT,Vector{ElT},length(inds)}
  return EmptyTensor(StoreT, inds)
end

fulltype(::Type{EmptyStorage{ElT,StoreT}}) where {ElT,StoreT} = StoreT
fulltype(T::EmptyStorage) = fulltype(typeof(T))

fulltype(T::Tensor) = fulltype(typeof(T))

# From an EmptyTensor, return the closest Tensor type
function fulltype(::Type{TensorT}) where {TensorT<:Tensor}
  return Tensor{
    eltype(TensorT),ndims(TensorT),fulltype(storetype(TensorT)),indstype(TensorT)
  }
end

function fulltype(
  ::Type{ElR}, ::Type{<:Tensor{ElT,N,EStoreT,IndsT}}
) where {ElR,ElT<:Number,N,EStoreT<:EmptyStorage{ElT,StoreT},IndsT} where {StoreT}
  return Tensor{ElR,N,similartype(StoreT, ElR),IndsT}
end

function zeros(T::TensorT) where {TensorT<:EmptyTensor}
  TensorR = fulltype(TensorT)
  return zeros(TensorR, inds(T))
end

function zeros(::Type{ElT}, T::TensorT) where {ElT,TensorT<:EmptyTensor}
  TensorR = fulltype(ElT, TensorT)
  return zeros(TensorR, inds(T))
end

function insertblock(T::EmptyTensor{<:Number,N}, block) where {N}
  R = zeros(T)
  insertblock!(R, Block(block))
  return R
end

insertblock!!(T::EmptyTensor{<:Number,N}, block) where {N} = insertblock(T, block)

# Special case with element type of EmptyNumber: storage takes the type
# of the input
@propagate_inbounds function _setindex(T::EmptyTensor{EmptyNumber}, x, I...)
  R = zeros(typeof(x), T)
  R[I...] = x
  return R
end

@propagate_inbounds function _setindex(T::EmptyTensor, x, I...)
  R = zeros(T)
  R[I...] = x
  return R
end

@propagate_inbounds function setindex(T::EmptyTensor, x, I...)
  return _setindex(T, x, I...)
end

# This is needed to fix an ambiguity error with ArrayInterface.jl
# https://github.com/ITensor/NDTensors.jl/issues/62
@propagate_inbounds function setindex(T::EmptyTensor, x, I::Int...)
  return _setindex(T, x, I...)
end

setindex!!(T::EmptyTensor, x, I...) = setindex(T, x, I...)

# Version of contraction where output storage is empty
function contract!!(R::EmptyTensor, labelsR, T1::Tensor, labelsT1, T2::Tensor, labelsT2)
  RR = contract(T1, labelsT1, T2, labelsT2, labelsR)
  return RR
end

# When one of the tensors is empty, return an empty
# tensor.
# XXX: make sure `R` is actually correct!
function contract!!(
  R::EmptyTensor, labelsR, T1::EmptyTensor, labelsT1, T2::Tensor, labelsT2
)
  return R
end

# When one of the tensors is empty, return an empty
# tensor.
# XXX: make sure `R` is actually correct!
function contract!!(
  R::EmptyTensor, labelsR, T1::Tensor, labelsT1, T2::EmptyTensor, labelsT2
)
  return R
end

function contract!!(
  R::EmptyTensor, labelsR, T1::EmptyTensor, labelsT1, T2::EmptyTensor, labelsT2
)
  return R
end

# For ambiguity with versions in combiner.jl
function contract!!(
  R::EmptyTensor, labelsR, T1::CombinerTensor, labelsT1, T2::Tensor, labelsT2
)
  RR = contract(T1, labelsT1, T2, labelsT2, labelsR)
  return RR
end

# For ambiguity with versions in combiner.jl
function contract!!(
  R::EmptyTensor, labelsR, T1::Tensor, labelsT1, T2::CombinerTensor, labelsT2
)
  RR = contract(T1, labelsT1, T2, labelsT2, labelsR)
  return RR
end

# For ambiguity with versions in combiner.jl
function contract!!(
  R::EmptyTensor, labelsR, T1::EmptyTensor, labelsT1, T2::CombinerTensor, labelsT2
)
  RR = contraction_output(T1, labelsT1, T2, labelsT2, labelsR)
  return RR
end

promote_rule(::Type{EmptyNumber}, ::Type{T}) where {T<:Number} = T

function promote_rule(
  ::Type{T1}, ::Type{T2}
) where {T1<:EmptyStorage{EmptyNumber},T2<:TensorStorage}
  return T2
end
function promote_rule(::Type{T1}, ::Type{T2}) where {T1<:EmptyStorage,T2<:TensorStorage}
  return promote_type(similartype(T2, eltype(T1)), T2)
end

function contraction_output(T1::EmptyTensor, T2::EmptyTensor, is)
  fulltypeR = contraction_output_type(fulltype(T1), fulltype(T2), typeof(is))
  storagetypeR = storagetype(fulltypeR)
  emptystoragetypeR = emptytype(storagetypeR)
  return Tensor(emptystoragetypeR(), is)
end

function contraction_output(T1::Tensor, T2::EmptyTensor, is)
  fulltypeR = contraction_output_type(typeof(T1), fulltype(T2), typeof(is))
  storagetypeR = storagetype(fulltypeR)
  emptystoragetypeR = emptytype(storagetypeR)
  return Tensor(emptystoragetypeR(), is)
end

function contraction_output(T1::EmptyTensor, T2::Tensor, is)
  fulltypeR = contraction_output_type(fulltype(T1), typeof(T2), typeof(is))
  storagetypeR = storagetype(fulltypeR)
  emptystoragetypeR = emptytype(storagetypeR)
  return Tensor(emptystoragetypeR(), is)
end

function permutedims!!(R::Tensor, T::EmptyTensor, perm::Tuple, f::Function=(r, t) -> t)
  RR = convert(promote_type(typeof(R), typeof(T)), R)
  RR = permutedims!!(RR, RR, ntuple(identity, Val(ndims(R))), (r, t) -> f(r, false))
  return RR
end

function permutedims!!(R::EmptyTensor, T::Tensor, perm::Tuple, f::Function=(r, t) -> t)
  RR = similar(promote_type(typeof(R), typeof(T)), inds(R))
  RR = permutedims!!(RR, T, perm, (r, t) -> f(false, t))
  return RR
end

function permutedims!!(R::EmptyTensor, T::EmptyTensor, perm::Tuple, f::Function=(r, t) -> t)
  RR = convert(promote_type(typeof(R), typeof(T)), R)
  return RR
end

function show(io::IO, mime::MIME"text/plain", T::EmptyTensor)
  summary(io, T)
  return println(io)
end

# XXX: this seems a bit strange and fragile?
# Takes the type very literally.
function HDF5.read(
  parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{StoreT}
) where {StoreT<:EmptyStorage}
  g = open_group(parent, name)
  typestr = string(StoreT)
  if read(attributes(g)["type"]) != typestr
    error("HDF5 group or file does not contain $typestr data")
  end
  return StoreT()
end

function HDF5.write(
  parent::Union{HDF5.File,HDF5.Group}, name::String, ::StoreT
) where {StoreT<:EmptyStorage}
  g = create_group(parent, name)
  attributes(g)["type"] = string(StoreT)
  return attributes(g)["version"] = 1
end


text:
@inline function add_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yhi, ylo = y
    yhi, ylo = two_sum(x, yhi, ylo)
    return yhi, ylo
end

@inline function sub_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yhi, ylo = y
    yhi, ylo = two_diff(x, yhi, ylo)
    return yhi, ylo
end

#= not quite right?
# Algorithm 12 from Tight and rigourous error bounds.  relative error <= 5u²
@inline function mul_fpdd_dd(x::T, y::Tuple{T,T}) where T<:IEEEFloat
    yhi, ylo = y
    hi, lo = two_prod(yhi, x)
    t = ylo * x
    t = t + lo
    hi, lo = two_hilo_sum(yhi, t)
    return hi, lo
end
=#
@inline function mul_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    hi, lo = y
    hihi, hilo = two_prod(x, hi)
    lohi, lolo = two_prod(x, lo)
    hi, lo = two_sum(hihi, hilo, lohi, lolo)
    return hi, lo
end

@inline function dvi_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yinv = inv_dd_dd(y)
    zhi, zlo = mul_fpdd_dd(x, yinv)
    return zhi, zlo
end

@inline function dvi_fpdd_dd_fast(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yinv = inv_dd_dd_fast(y)
    zhi, zlo = mul_fpdd_dd(x, yinv)
    return zhi, zlo
end


text:
# check property, kD blocks, intervals, given delta, dense time
using Reachability

check_kD_dense(model::String, delta::Float64) =
    check_kD_dense([model], [delta])

function check_kD_dense(models::Vector{String}, deltas::Vector{Float64})
    # load models
    models_loaded = false
    for model in models
        if !isdefined(Main, Symbol(model))
            if !models_loaded
                model_library = include(@relpath "model_library.jl")
            end
            path = (@relpath "../") * model_library[model]
            include(path)
            models_loaded = true
        end
    end
    if models_loaded
        println("models have been loaded; please run again for analysis")
        return
    end

    println("-- benchmark suite 'check_kD_dense' --")

    @assert length(models) == length(deltas) "the input lengths must be equal"
    for i in 1:length(models)
        model = models[i]
        delta = deltas[i]

        println("- analyzing model '$model' -")

        # raw model
        func = getfield(Main, Symbol(model))
        S, options_raw = func(:mode => "check")

        # options
        n = MathematicalSystems.statedim(S)
        dict_raw = options_raw.dict
        dict_raw[:verbosity] = "info"
        dict_raw[:δ] = delta

        for i in 1:2
            dict = copy(dict_raw)
            if i == 1
                # warm-up run
                dict[:N] = 3
            else
                # benchmark settings
                dict[:T] = 20.
                dict[:logfile] = "$model-check-kD-dense-givenstep.txt"
            end
            result = solve(S, Options(dict))
        end
        println()
    end
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
make_entry(;id=1) = (id=id,
                     IntField=1, FloatField=0.01,
                     StringField="a", StringField2="b",
                     RandFloat=0.0, RandString="c")

nt_change_key(nt, key, val) = begin
    props = propertynames(nt)
    idxkey = findall(isequal(key), props)[1]  # it is certain, always one
    vals = [v for v in nt];
    vals[idxkey] = val
    return NamedTuple{props, Tuple{(typeof(v) for v in vals)...}}(Tuple(vals))
end

@testset "Environment: creation" begin
    # Basic
    for gfunc in [generate_sample_config_1, generate_sample_config_2]
        cfg = mktemp() do path, io
            write(io, gfunc())
            flush(io)
            parse_configuration(path)
        end
        @test build_search_env(cfg) isa SearchEnv
    end
end

@testset "Environment: push!/pop! primitives" begin
    cfg = mktemp() do path, io
        write(io, generate_sample_config_2())
        flush(io)
        parse_configuration(path)
    end
    for index_type in [:naive, :brutetree, :kdtree, :hnsw, :ivfadc]
        for i in eachindex(cfg.searcher_configs)
            cfg.searcher_configs[i] = nt_change_key(cfg.searcher_configs[i], :search_index, index_type)
            cfg.searcher_configs[i] = nt_change_key(cfg.searcher_configs[i], :search_index_arguments, [])
            cfg.searcher_configs[i] = nt_change_key(cfg.searcher_configs[i], :search_index_kwarguments, Dict{Symbol,Any}())
        end
        env = build_search_env(cfg)
        initial_length = length(env.dbdata)
        index_length(x) = length(x.index)
        get_vec_type(env::SearchEnv{T}) where {T} = T

        @test env isa SearchEnv
        if index_type in [:naive, :brutetree, :ivfadc]  # push!/pop! supported
            # push!
            ientry = make_entry(id=initial_length+1)
            push!(env, ientry)
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length + 1

            # pop!
            oentry, embedded = pop!(env)
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length
            index_type !== :ivfadc && @test ientry == oentry
            @test typeof(embedded) <: AbstractVector{<:AbstractVector{get_vec_type(env)}}

            # pushfirst!
            ientry = make_entry(id=1)
            pushfirst!(env, ientry)
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length + 1

            # popfirst!
            oentry, embedded = popfirst!(env)
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length
            index_type !== :ivfadc && @test ientry == oentry
            @test typeof(embedded) <: AbstractVector{<:AbstractVector{get_vec_type(env)}}

            # deleteat!
            deleteat!(env, [1])
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length - 1
        else
            # push!
            ientry = make_entry(id=initial_length+1)
            push!(env, ientry)
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length

            # pop!
            oentry = pop!(env)  # returns a nothing
            @test oentry === nothing
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length

            # pushfirst!
            ientry = make_entry(id=1)
            pushfirst!(env, ientry)
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length

            # popfirst!
            oentry = popfirst!(env)  # returns a nothing
            @test oentry === nothing
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length

            # deleteat!
            deleteat!(env, [1])
            @test all(length(env.dbdata) .== map(index_length, env.searchers))
            @test length(env.dbdata) == initial_length
        end
    end
end


text:
"""
	$(SIGNATURES)

Initialize an empty, preallocated `WorkHistories` object.
"""
function WorkHistories(F1 :: DataType, nSim :: Integer)
    return WorkHistories(
        zeros(UInt8, nSim),
        zeros(TimeInt, nSim),
        zeros(F1, nSim), zeros(F1, nSim),
        zeros(F1, nSim), zeros(F1, nSim),
        zeros(F1, nSim), zeros(F1, nSim)
    )
end

"""
	$(SIGNATURES)

Create `WorkHistories` with work start endowments filled in.
"""
function WorkHistories(iSchoolV, workStartAgeV,
    hWorkStartV :: Vector{F1}, kWorkStartV :: Vector{F1}) where F1 <: AbstractFloat

    nSim = length(iSchoolV);
    return WorkHistories(
        UInt8.(iSchoolV), TimeInt.(workStartAgeV), hWorkStartV, kWorkStartV,
        zeros(F1, nSim), zeros(F1, nSim),
        zeros(F1, nSim), zeros(F1, nSim)
    )
end

Base.show(io :: IO, wh :: WorkHistories{F1}) where F1 =
    print(io, "WorkHistories of length $(length(wh)).")

Base.length(wh :: WorkHistories) = length(wh.workStartAgeV);
school_index(wh :: WorkHistories) = wh.iSchoolV;
work_start_ages(wh :: WorkHistories) = wh.workStartAgeV;
h_work_start(wh :: WorkHistories{F1}) where F1 = wh.hWorkStartV;
k_work_start(wh :: WorkHistories{F1}) where F1 = wh.kWorkStartV;
n_school(wh :: WorkHistories) = maximum(school_index(wh));

function validate_wh(wh :: WorkHistories{F1}) where F1
    isValid = true;
    if any_less(wh.workStartAgeV, 1)  ||  any_greater(wh.workStartAgeV, 9)
        isValid = false;
        @warn "Invalid work start ages"
    end
    if any_at_most(wh.earnWorkStartV, zero(F1))
        isValid = false;
        @warn "Earnings at work start should be positive"
    end
    if any_at_most(wh.ltEarnV, zero(F1))
        isValid = false;
        @warn "Lifetime earnings should be positive"
    end
    if any_at_most(wh.cWorkStartV, zero(F1))
        isValid = false
        @warn "Negative consumption at work start"
    end
    if any_at_most(wh.muWealthV, zero(F1))
        isValid = false;
        @warn "Negative marginal utility of wealth"
    end
    return isValid
end

"""
	$(SIGNATURES)

Simulate work histories. Input are `WorkHistories` that are initialized with workers' initial conditions: work start age, schooling, initial h and k.
"""
function simulate_workers!(wh :: WorkHistories{F1},
    workerV :: Vector{Worker{F1}}) where F1

    for iSchool = 1 : n_school(wh)
        sIdxV = findall(s -> s == iSchool,  school_index(wh));
        if !isempty(sIdxV)
            simulate_one_worker!(wh, workerV[iSchool], sIdxV);
        end
    end

    @assert validate_wh(wh)
    return nothing
end

"""
	$(SIGNATURES)

Simulate workers for one school level. Workers in rows `sIdxV`.
"""
function simulate_one_worker!(wh :: WorkHistories{F1}, wk :: Worker{F1},
    sIdxV :: AbstractVector{I1}) where {F1 <: AbstractFloat, I1 <: Integer}

    startAgeV = unique(wh.workStartAgeV[sIdxV]);
    for workStartAge in startAgeV
        idxV = findall(a -> a == workStartAge, work_start_ages(wh)[sIdxV]);
        simulate_one_work_start!(wh, wk, workStartAge, sIdxV[idxV]);
    end
    return nothing
end

# Work start age must be common to all
function simulate_one_work_start!(wh :: WorkHistories{F1}, wk :: Worker{F1},
    workStartAge :: Integer, idxV :: AbstractVector{I1}) where {F1 <: AbstractFloat, I1 <: Integer}

    # Earnings at work start
    # Earnings profile, up to `h` factor
    earnWorkStart = earn_profile(wk, workStartAge, 1.0; experV = 1);
    wh.earnWorkStartV[idxV] .= h_work_start(wh)[idxV] .* earnWorkStart;

    wh.ltEarnV[idxV] .= lifetime_earnings(wk, workStartAge,
        h_work_start(wh)[idxV]);
    # Bound below in case we have bad parameter values during testing
    ltIncomeV = max.(0.1, wh.ltEarnV[idxV] .+ k_work_start(wh)[idxV]);
    # Consumption at work start
    wh.cWorkStartV[idxV] .= cons_age1(wk, workStartAge, ltIncomeV);
    # Marginal utility of wealth as of work start
    wh.muWealthV[idxV] .= mu_wealth(wk, workStartAge, ltIncomeV);
    return nothing
end

function make_test_work_histories(nSchool :: Integer, nSim :: Integer)
    rng = MersenneTwister(43);
    iSchoolV = rand(rng, UInt8.(1 : nSchool), nSim);
    workStartAgeV = rand(rng, TimeInt.(1:5), nSim);
    hWorkStartV = 1.0 .+ 2.0 .* rand(rng, nSim);
    kWorkStartV = rand(rng, nSim) .- 0.1;
    earnWorkStartV = 9.0 .+ 3.0 .* rand(rng, nSim);
    ltEarnV = 10.5 .* earnWorkStartV;
    muWealthV = 5.0 ./ earnWorkStartV;
    cWorkStartV = 0.7 .* earnWorkStartV;

    wh = WorkHistories(iSchoolV, workStartAgeV, hWorkStartV, kWorkStartV,
        earnWorkStartV, ltEarnV, muWealthV, cWorkStartV);
    @assert validate_wh(wh)  "Invalid Test work history"
    return wh
end

# ----------------


text:
function affinehull{T<:HRep}(p::T)
  T(eqs(p))
end

# Always type of first arg
function Base.intersect{T1<:HRep, T2<:HRep}(p1::T1, p2::T2)
  if eltype(T1) != eltype(T2)
    error("Cannot take the intersection of polyhedra of different element type")
  end
  if fulldim(T1) != fulldim(T2)
    error("Cannot take the intersection of polyhedra of different dimension")
  end
  T1(HRepIterator([p1, p2]))
end

# Always type of first arg
function (+){T1<:VRep, T2<:VRep}(p1::T1, p2::T2)
  if eltype(T1) != eltype(T2)
    error("Cannot take the intersection of polyhedra of different element type")
  end
  if fulldim(T1) != fulldim(T2)
    error("Cannot take the intersection of polyhedra of different dimension")
  end
  T1(VRepIterator([p1, p2]))
end

# p1 has priority
function usehrep(p1::Polyhedron, p2::Polyhedron)
  hrepiscomputed(p1) && (!vrepiscomputed(p1) || hrepiscomputed(p2))
end

# Always type of first arg
@generated function (*){T1<:Rep, T2<:Rep}(p1::T1, p2::T2)
  if eltype(T1) != eltype(T2)
    error("Cannot take the cartesian product between polyhedra of different element type")
  end
  T = eltype(T1)
  N1 = fulldim(T1)
  N2 = fulldim(T2)
  hashrep = T1 <: HRepresentation || T2 <: HRepresentation
  hasvrep = T1 <: VRepresentation || T2 <: VRepresentation
  Tout = changefulldim(T1, fulldim(T1)+fulldim(T2))
  f = (i, x) -> zeropad(x, i == 1 ? N2 : -N1)
  if hashrep && hasvrep
    error("Cannot take the cartesian product between a H-Representation and a V-Representation")
  elseif hashrep || (!hasvrep && usehrep(p1, p2))
    # TODO fastdecompose
    :(Tout(HRepIterator([p1, p2], f)))
  else
    # TODO fastdecompose
    :(Tout(VRepIterator([p1, p2], f)))
  end
end

function (*){RepT<:HRep}(hrep::RepT, P::AbstractMatrix)
  if size(P, 1) != fulldim(T)
    error("The number of rows of P must match the dimension of the H-representation")
  end
  Tout = mypromote_type(eltype(RepT), eltype(P))
  RepTout = changeboth(RepT, size(P, 2), Tout)
  f = (i, h) -> h * P
  if decomposedhfast(hrep)
    eqs = eqs(hrep, f)
    ineqs = ineqs(hrep, f)
    if RepT <: HRepresentation
      Tout(ineqs=ineqs, eqs=eqs)
    else
      polyhedron(ineqs, eqs, getlibraryfor(p, Tout))
    end
  else
    hreps = hrep(hrep, f)
    if RepT <: HRepresentation
      Tout(hreps)
    else
      polyhedron(hreps, getlibraryfor(p, Tout))
    end
  end
end
function (*){RepT<:VRep}(P::AbstractMatrix, vrep::RepT)
  if size(P, 2) != fulldim(T)
    error("The number of rows of P must match the dimension of the H-representation")
  end
  Tout = mypromote_type(eltype(RepT), eltype(P))
  RepTout = changeboth(RepT, size(P, 1), Tout)
  f = (i, v) -> P * v
  if decomposedvfast(vrep)
    points = points(vrep, f)
    rays = rays(vrep, f)
    if RepT <: VRepresentation
      Tout(points=points, rays=rays)
    else
      polyhedron(points, rays, getlibraryfor(p, Tout))
    end
  else
    vreps = vrep(vrep, f)
    if RepT <: VRepresentation
      Tout(vreps)
    else
      polyhedron(vreps, getlibraryfor(p, Tout))
    end
  end
end

function Base.round{N,T<:AbstractFloat}(rep::HRepresentation{N,T})
  f = (i, h) -> round(h)
  if decomposedfast(rep)
    typeof(rep)(eqs=eqs(rep, f), ineqs=ineqs(rep, f))
  else
    typeof(rep)(hrep(rep, f))
  end
end
function Base.round{N,T<:AbstractFloat}(rep::VRepresentation{N,T})
  f = (i, v) -> round(v)
  if decomposedfast(rep)
    typeof(rep)(eqs=eqs(rep, f), ineqs=ineqs(rep, f))
  else
    typeof(rep)(vrep(rep, f))
  end
end

function gethredundantindices(hrep::HRep)
  red = IntSet([])
  for i in 1:nhreps(hrep)
    if ishredundant(hrep, i)
      push!(red, i)
    end
  end
  red
end
function gethstronglyredundantindices(hrep::HRep)
  red = IntSet([])
  for i in 1:nhreps(hrep)
    if ishstronglyredundant(p, i)[1]
      push!(red, i)
    end
  end
  red
end
function getvredundantindices(vrep::VRep)
  red = IntSet([])
  for i in 1:nvreps(vrep)
    if isvredundant(p, i)
      push!(red, i)
    end
  end
  red
end
function getvstronglyredundantindices(hrep::HRep)
  red = IntSet([])
  for i in 1:nvreps(p)
    if isvstronglyredundant(p, i)
      push!(red, i)
    end
  end
  red
end