text:
using DataStructures: SortedSet

"""
    WordleGame

A struct tracking the state of a game of Wordle.
"""
struct WordleGame
    target::String
    number::Union{Int, Nothing}
    guesses::Vector{WordleGuess}
    hard::Bool

    function WordleGame(target::String, number::Union{Int, Nothing} = nothing,
               guesses::Vector{WordleGuess} = WordleGuess[]; hard = false)

        target = lowercase(target)

        if target ∉ VALID_WORD_LIST
            error("""Target word "$target" is not a valid word""")
        end

        if !isnothing(number) && !(0 < number <= LATEST_WORDLE_NUMBER)
            error("Given Wordle number is too large ($number).")
        end

        if length(guesses) > 6
            error("There are too many guesses ($(length(guesses))")
        end

        new(target, number, guesses, hard)
    end
end

WordleGame(number::Int; hard = false) = WordleGame(WORDLE_LIST[number], number; hard = hard)
WordleGame(; hard = false) = WordleGame(LATEST_WORDLE_NUMBER; hard = hard)

nguess(game::WordleGame) = length(game.guesses)
target(game::WordleGame) = game.target

function show(io::IO, game::WordleGame)
    header = ["Wordle"]

    if !isnothing(game.number)
        push!(header, string(game.number))
    end

    if game.hard
        # add * for hard mode
        push!(header, "$(nguess(game))/6*")
    else
        push!(header, "$(nguess(game))/6")
    end

    print(io, join(header, ' '))

    if nguess(game) > 0
        println(io)
        println(io)

        for (i, guess) in enumerate(game.guesses)
            if i < nguess(game)
                println(guess)
            else
                print(guess)
            end
        end
    end
end

function (∈)(word::String, game::WordleGame)
    word in getproperty.(game.guesses, :guess)
end

function guess(game::WordleGame, word::String)
    word = lowercase(word)

    if length(word) > 5 || word ∉ VALID_WORD_LIST
        error("$word is not a valid guess")
    end

    if target(game) ∈ game
        error("The game is already over!")
    end

    if game.hard
        letter_outcomes = available_letters(game)

        # might need to account for multiple occurrences of a letter
        required = Set(filter('a':'z') do l
            index = l - 'a' + 1
            letter_outcomes[index] ∈ (CORRECT, PRESENT)
        end)

        if !all(occursin.(required, word))
            error("Invalid guess for hard mode.")
        end
    end

    # initialize a results array to be all incorrect
    results = fill(INCORRECT, 5)

    # we only need to consider letters that are present in both target and guess
    for m in intersect(target(game), word)
        # find the sets of positions for each matched letter
        target_positions = SortedSet(findall(m, target(game)))
        guess_positions = SortedSet(findall(m, word))

        # for each exactly matching position, mark the letter as correct and
        # remove the position from the position sets
        for correct_position in intersect(target_positions, guess_positions)
            results[correct_position] = CORRECT
            delete!(target_positions, correct_position)
            delete!(guess_positions, correct_position)
        end

        # pair off the remainining positions and mark them as present
        # zip will stop after the last pair.
        for (_, guess_position) in zip(target_positions, guess_positions)
            results[guess_position] = PRESENT
        end
    end

	WordleGuess(word, results)
end

function guess!(game::WordleGame, word::String)
    if nguess(game) == 6
        error("The game is already over!")
    end

    push!(game.guesses, guess(game, word))
    game
end

function available_letters(game::WordleGame)
    outcomes = [UNGUESSED for _ in 'a':'z']

    for guess in game.guesses
        for (outcome, letter) in zip(guess.result, guess.guess)
            letter_index = letter - 'a' + 1

            if outcome == CORRECT
                outcomes[letter_index] = CORRECT
            elseif outcome == PRESENT && outcomes[letter_index] != CORRECT
                # if the outcome for a letter in a given guess is present but it
                # was previously correct, then it should still be marked as
                # correct in this list
                outcomes[letter_index] = PRESENT
            else # incorrect
                if outcomes[letter_index] ∉ (CORRECT, PRESENT)
                    outcomes[letter_index] = INCORRECT
                end
            end
        end
    end

    outcomes
end

function print_available_letters(game::WordleGame)
    print("""
    $(join('a':'z', ' '))
    $(available_letters(game) |> join)
    """)
end


text:
using Base.Test
using JavaCall
using Compat

versioninfo();

# JavaCall.init(["-Djava.class.path=$(joinpath(Pkg.dir(), "JavaCall", "test"))"])
JavaCall.init(["-verbose:gc","-Djava.class.path=$(joinpath(Pkg.dir(), "JavaCall", "test"))"])

a=JString("how are you")
@test a.ptr != C_NULL
@test 11==ccall(JavaCall.jnifunc.GetStringUTFLength, jint, (Ptr{JavaCall.JNIEnv}, Ptr{Void}), JavaCall.penv, a.ptr)
b=ccall(JavaCall.jnifunc.GetStringUTFChars, Ptr{UInt8}, (Ptr{JavaCall.JNIEnv}, Ptr{Void}, Ptr{Void}), JavaCall.penv, a.ptr, C_NULL)
@test unsafe_string(b) == "how are you"

# Test parameter passing
T = @jimport Test
@test 10 == jcall(T, "testShort", jshort, (jshort,), 10)
@test 10 == jcall(T, "testInt", jint, (jint,), 10)
@test 10 == jcall(T, "testLong", jlong, (jlong,), 10)
@test typemax(jint) == jcall(T, "testInt", jint, (jint,), typemax(jint))
@test typemax(jlong) == jcall(T, "testLong", jlong, (jlong,), typemax(jlong))
@test "Hello Java"==jcall(T, "testString", JString, (JString,), "Hello Java")
@test @compat Float64(10.02) == jcall(T, "testDouble", jdouble, (jdouble,), 10.02) #Comparing exact float representations hence ==
@test @compat Float32(10.02) == jcall(T, "testFloat", jfloat, (jfloat,), 10.02)
@test realmax(jdouble) == jcall(T, "testDouble", jdouble, (jdouble,), realmax(jdouble))
@test realmax(jfloat) == jcall(T, "testFloat", jfloat, (jfloat,), realmax(jfloat))

c=JString(C_NULL)
@test isnull(c)
@test "" == jcall(T, "testString", JString, (JString,), c)

# Test calling static methods
jlm = @jimport "java.lang.Math"
@test_approx_eq 1.0 jcall(jlm, "sin", jdouble, (jdouble,), pi/2)
@test_approx_eq 1.0 jcall(jlm, "min", jdouble, (jdouble,jdouble), 1,2)
@test 1 == jcall(jlm, "abs", jint, (jint,), -1)

#Test instance creation
jnu = @jimport java.net.URL
gurl = jnu((JString,), "http://www.google.com")
@test "www.google.com"==jcall(gurl, "getHost", JString,())

#Test instance methods
jni=@jimport java.net.URI
guri=jcall(gurl, "toURI", jni,())
@test typeof(guri)==jni

h=jcall(guri, "hashCode", jint,())
typeof(h)==jint

#Test NULL
H=@jimport java.util.HashMap
a=jcall(T, "testNull", H, ())
@test_throws ErrorException jcall(a, "toString", JString, ())

# Arrays
j_u_arrays = @jimport java.util.Arrays
@test 3 == jcall(j_u_arrays, "binarySearch", jint, (Array{jint,1}, jint), [10,20,30,40,50,60], 40)
@test 2 == jcall(j_u_arrays, "binarySearch", jint, (Array{JObject,1}, JObject), ["123","abc","uvw","xyz"], "uvw")

a=jcall(j_u_arrays, "copyOf", Array{jint, 1}, (Array{jint, 1}, jint), [1,2,3], 3)
@test typeof(a) == Array{jint, 1}
@test a[1] == @compat Int32(1)
@test a[2] == @compat Int32(2)
@test a[3] == @compat Int32(3)

a=jcall(j_u_arrays, "copyOf", Array{JObject, 1}, (Array{JObject, 1}, jint), ["a","b","c"], 3)
@test 3==length(a)
@test "a"==unsafe_string(convert(JString, a[1]))
@test "b"==unsafe_string(convert(JString, a[2]))
@test "c"==unsafe_string(convert(JString, a[3]))

#Test for Dates

jd = @jimport(java.util.Date)(())
jcal = @jimport(java.util.GregorianCalendar)(())
jsd =  @jimport(java.sql.Date)((jlong,),round(jlong, time()))

@assert typeof(convert(Dates.DateTime, jd)) == Dates.DateTime
@assert typeof(convert(Dates.DateTime, jcal)) == Dates.DateTime
@assert typeof(convert(Dates.DateTime, jsd)) == Dates.DateTime
nulldate = @jimport(java.util.Date)(C_NULL)
@assert Dates.year(convert(Dates.DateTime, nulldate)) == 1970
nullcal = @jimport(java.util.GregorianCalendar)(C_NULL)
@assert Dates.year(convert(Dates.DateTime, nullcal)) == 1970

Dates.year(convert(Dates.DateTime, nullcal)) == 1970
#Test for Map conversion

JHashMap = @jimport(java.util.HashMap)
p = JHashMap(())
a=@compat Dict("a"=>"A", "b"=>"B")
b=convert(@jimport(java.util.Map), JString, JString, a)
@assert jcall(b, "size", jint, ()) == 2

#Inner Classes
TestInner = @jimport(Test$TestInner)
Test = @jimport(Test)
t=Test(())
inner = TestInner((Test,), t)
@assert jcall(inner, "innerString", JString, ()) == "from inner"

#Static Field Access
@assert jfield(@jimport(java.lang.Math), "E", jdouble) == 2.718281828459045
@assert jfield(@jimport(java.lang.Math), "PI", jdouble) == 3.141592653589793
@assert jfield(@jimport(java.text.NumberFormat), "INTEGER_FIELD", jint) == 0
Locale = @jimport java.util.Locale
lc = jfield(@jimport(java.util.Locale), "CANADA", @jimport(java.util.Locale))
#Instance field access
#Disabled for now. Need to verify stability
@assert jfield(@jimport(java.util.logging.Logger), "GLOBAL_LOGGER_NAME", JString ) == "global"
@assert jcall(lc, "getCountry", JString, ()) == "CA"
@assert jfield(t, "integerField", jint) == 100
@assert jfield(t, "stringField", JString) == "A STRING"

# Test Memory allocation and de-allocatios
# the following loop fails with an OutOfMemoryException in the absence of de-allocation
# However, since Java and Julia memory are not linked, and manual gc() is required.
gc()
for i in 1:100000
	a=JString("A"^10000); #deleteref(a);
	if (i%10000 == 0); gc(); end
end

#Test for Issue #8
@test_throws ErrorException jcall(jlm, "sinx", jdouble, (jdouble,), 1.0)
@test_throws ErrorException jcall(jlm, "sinx", jdouble, (jdouble,), 1.0)

@test length(listmethods(JString("test"))) >= 72
@test length(listmethods(JString("test"), "indexOf")) >= 3
m = listmethods(JString("test"), "indexOf")[1]
@test getname(getreturntype(m)) == "int"
@test [getname(typ) for typ in getparametertypes(m)] == ["java.lang.String", "int"]

#Test for double free bug, #20
#Fix in #28. The following lines will segfault without the fix
JHashtable = @jimport java.util.Hashtable
JProperties = @jimport java.util.Properties
ta_20=Any[]
for i=1:100; push!(ta_20, convert(JHashtable, JProperties((),))); end
gc(); gc()
for i=1:100; @test jcall(ta_20[i], "size", jint, ()) == 0; end

# Test array conversions
jobj = jcall(T, "testArrayAsObject", JObject, ())
arr = convert(Array{Array{UInt8, 1}, 1}, jobj)
@test ["Hello", "World"] == map(Compat.String, arr)

# At the end, unload the JVM before exiting
JavaCall.destroy()

ext: jl


text:
# ============================================================
# Structs of Interface Fluxes
# compatible with solution structs
# ============================================================

struct Flux{T1,T2,ND} <: AbstractFlux
    n::T1
    w::T2
    fw::T2
end

struct Flux1F{T1,T2,T3,ND} <: AbstractFlux
    n::T1
    w::T2
    fw::T2
    ff::T3
end

struct Flux2F{T1,T2,T3,ND} <: AbstractFlux
    n::T1
    w::T2
    fw::T2
    fh::T3
    fb::T3
end

function Flux1D(fw::AA)
    n = ones(axes(fw)[end])
    w = deepcopy(fw)
    if w[1] isa Number
        w .= 0.0
    else
        for e in w
            e .= 0.0
        end
    end

    return Flux{typeof(n),typeof(w),1}(n, w, fw)
end

function Flux1D(fw::AA, ff::AA)
    n = ones(axes(fw)[end])
    w = deepcopy(fw)
    if w[1] isa Number
        w .= 0.0
    else
        for e in w
            e .= 0.0
        end
    end

    return Flux1F{typeof(n),typeof(w),typeof(ff),1}(n, w, fw, ff)
end

function Flux1D(fw::AA, fh::AA, fb::AA)
    n = ones(axes(fw)[end])
    w = deepcopy(fw)
    if w[1] isa Number
        w .= 0.0
    else
        for e in w
            e .= 0.0
        end
    end

    return Flux2F{typeof(n),typeof(w),typeof(fh),1}(n, w, fw, fh, fb)
end

function Flux2D(n::AA, fw::AA)
    w = deepcopy(fw)
    if w[1] isa Number
        w .= 0.0
    else
        for e in w
            e .= 0.0
        end
    end

    return Flux{typeof(n),typeof(w),2}(n, w, fw)
end

function Flux2D(n::AA, fw::AA, ff::AA)
    w = deepcopy(fw)
    if w[1] isa Number
        w .= 0.0
    else
        for e in w
            e .= 0.0
        end
    end

    return Flux1F{typeof(n),typeof(w),typeof(ff),2}(n, w, fw, ff)
end

function Flux2D(n::AA, fw::AA, fh::AA, fb::AA)
    w = deepcopy(fw)
    if w[1] isa Number
        w .= 0.0
    else
        for e in w
            e .= 0.0
        end
    end

    return Flux2F{typeof(n),typeof(w),typeof(fh),2}(n, w, fw, fh, fb)
end


text:
include("MultiArmBanditsEnv.jl")
include("RandomWalk1D.jl")
include("TigerProblemEnv.jl")
include("MontyHallEnv.jl")
include("RockPaperScissorsEnv.jl")
include("TicTacToeEnv.jl")
include("TinyHanabiEnv.jl")
include("PigEnv.jl")
include("KuhnPokerEnv.jl")
include("CartPoleEnv.jl")
include("MountainCarEnv.jl")
include("PendulumEnv.jl")
include("BitFlippingEnv.jl")
include("StockTradingEnv.jl")
include("GraphShortestPathEnv.jl")
include("SpeakerListenerEnv.jl")


text:
module BerkeleyABC
using ABC_jll

export start_abc, stop_abc, restart_abc, parse_timing

function start_abc(;load_aliases=true)
    ccall((:Abc_Start, libabc), Cvoid, ())
    pAbc = ccall((:Abc_FrameGetGlobalFrame, libabc), Ptr{UInt}, ())
    abc_cmd(cmd) = ccall((:Cmd_CommandExecute, libabc), Int, (Ptr{Int}, Cstring), pAbc, cmd)
    if load_aliases
        abc_cmd("source $(abc_rc)")
    end
    function cmd(args)
        status = 0
        res = ""
        mktemp() do _, io
            redirect_stdout(io) do
                status = abc_cmd(args)
                Base.Libc.flush_cstdio()
                seek(io, 0)
                res = readlines(io)
            end
        end
        return status, res
    end
end

function stop_abc()
    ccall((:Abc_Stop, libabc), Cvoid, ())
end

function restart_abc()
    stop_abc()
    start_abc()
end

# So yeah, kinda brittle parsing of the abc stime command
function parse_timing((st, res))
    if st != 0
        println("Error parsing timing: bad ABC result")
        return
    end
    # first find the correct line to parse
    line = ""
    for outer line in res
        all(occursin.(["Gates", "Area", "Delay"], (line,))) && break
    end

    toks = filter(!=(""), split(line, r"[ =]"))
    g, a, d = findfirst.(contains.(["Gates", "Area", "Delay"]), (toks,))

    gates = 0
    area = delay = 0.0

    try
        gates = parse(Int, split(toks[g+1], '\e')[1])
        area  = parse(Float64, split(toks[a+1], '\e')[1])
        delay = parse(Float64, split(toks[d+1], '\e')[1])
    catch
        println("Parsing area and delay failed on the following line:")
        println(res[1])
        println(toks)
    end

    (gates=gates, area=area, delay=delay)
end

end


text:
function preprocess(PData::PDataLDLt, H, g, n1, n2)
    # printstyled("dans preprocessLDLt ⤈ \n", color = :red)
    T = eltype(H)
    # @show T
    # @show H[1, 1]
    # @show H[1, 2]
    # @show H[2, 1]
    # @show H[2, 2]
    # @show g[1]
    # @show g[2]
    n, = size(g)
    # @show n
    global L = Matrix{T}(undef, n, n)
    global D = Matrix{T}(undef, n, n)
    global pp = Vector{Int}(undef, n)
    global ρ = nothing
    global ncomp = nothing

    # @show H

    try
        (L, D, pp, ρ, ncomp) = ldlt_symm(H, 'r')
    catch
        println("*******   Problem in LDLt")
        res = PDataLDLt()
        res.OK = false
        return res
    end

    # printstyled("dans preprocessLDLt apres ldlt_symm eltype(L) = $(eltype(L)) \n", color = :cyan)
    # @show L
    # @show D
    # @show L[1, 1]
    # @show L[2, 1]
    # @show L[1, 2]
    # @show L[2, 2]
    # @show D[1, 1]
    # @show D[2, 1]
    # @show D[1, 2]
    # @show D[2, 2]
    # @show pp
    # @show ρ
    # @show ncomp

    if true in isnan.(D)
        println("*******   Problem in D from LDLt: NaN")
        println(" cond (H) = $(cond(H))")
        res = PDataLDLt()
        res.OK = false
        return res
    end

    # @show L
    # @show D
    # @show pp

    # D = Hermitian(D)
    # @show D
    X = eigen(SymTridiagonal(D))
    Δ = X.values
    # @show Δ[1]
    # @show Δ[2]
    Q = X.vectors
    # @show Q[1, 1]
    # @show Q[2, 1]
    # @show Q[1, 2]
    # @show Q[2, 2]
    l_m, = findmin(Δ)
    # @show l_m
    ĝ = L \ (g[pp])
    # @show ĝ[1]
    # @show ĝ[2]
    # @show eltype(ĝ)
    # @show ĝ
    g̃ = Q' * ĝ
    # @show g̃[1]
    # @show g̃[2]
    # @show eltype(g̃)
    # @show g̃

    n_g = norm(g)
    # @show n_g
    λ = max(-l_m, 0.0)
    # @show eltype(λ)
    # @show eltype(λ)
    # @show λ
    # printstyled("avant de sortir de PreProcessLDLt ↑ \n", color = :red)

    PData.L = L
    PData.D = D
    PData.pp = pp
    PData.Δ = Δ
    PData.Q = Q
    PData.g̃ = g̃
    PData.λ = λ
    PData.success = true
    PData.OK = true

    return PData
end

function AInv(X::PDataLDLt, d̃::Array{T,1}) where {T}
    d̂ = X.Q * d̃
    u = X.L' \ d̂
    # @show u[invperm(X.pp)][1]
    # @show u[invperm(X.pp)][2]
    # printstyled("avant la sortie de AInv ⟰ \n", color = :red)
    return u[invperm(X.pp)]
end


text:
function stacktraces(d,h;normalize=true)
	h_out = [h[1]]
	h_out[1].trid = size(d[:,:],2)
	if (normalize == true)
		val = sum(d[:,:],2)/size(d[:,:],2)
	else
		val = sum(d[:,:],2)
	end
	return val, h_out
end

function SeisStack(in::String,out::String;key=["imx" "imy"],normalize=true)

	@compat parameters = Dict(:normalize=>true)
	SeisProcess(in,out,[stacktraces],[parameters],key=key)
	DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
	filename_data_out = join([DATAPATH out "@data@"])
	filename_headers_out = join([DATAPATH out "@headers@"])
	extent = ReadTextHeader(in)

	stream_h = open(filename_headers_out)
	@compat extent.n2 = round(Int,filesize(stream_h)/(4*length(fieldnames(Header))))
	close(stream_h)
	extent.d2 = 1.; extent.o2 = 0.; extent.label2 = ""; extent.unit2 = "";
	extent.n3 = 1; extent.d3 = 1.; extent.o3 = 0.; extent.label3 = ""; extent.unit3 = "";
	extent.n4 = 1; extent.d4 = 1.; extent.o4 = 0.; extent.label4 = ""; extent.unit4 = "";
	extent.n5 = 1; extent.d5 = 1.; extent.o5 = 0.; extent.label5 = ""; extent.unit5 = "";
	WriteTextHeader(out,extent,"native_float",4,filename_data_out,filename_headers_out)

end


text:
Pkg.clone(pwd())
Pkg.add("BinDeps")
using BinDeps
println(BinDeps.debug("Gtk"))
Pkg.build("GtkBuilderAid")
Pkg.test("GtkBuilderAid"; coverage = true)


text:
#=##############################################################################
# # DESCRIPTION
#     Examples of VTK formatting and usage of VTKtools
# # AUTHORSHIP
#   * Author    : Eduardo J Alvarez
#   * Email     : Edo.AlvarezR@gmail.com
#   * Created   : Nov 2017
#   * License   : MIT License
=###############################################################################
# ------------ MODULE IMPORTS ------------------------------------------------
import Statistics
using PyPlot

module_path,_ = splitdir(@__FILE__);      # Path to this module
data_path = module_path*"/../docs/data/"       # Path to data folder

import GeometricTools
gt = GeometricTools

# ------------ HEADERS -------------------------------------------------------
for header_name in ["simple", "airfoil", "wing", "blade", "taylor_wing"]
include("example_"*header_name*".jl")
end


text:
module postwidderTests

using Test
using InverseLaplace: postwid, _PWcoeffs

# Compare against results reported in Algorithm 368: Numerical inversion of Laplace transforms [D5].

#### F(s) = 1/sqrt(s), f(t) = 1/sqrt(π*t)

# test for single time point
setprecision(53)
t = 1.0
@test isapprox(postwid(s -> 1/sqrt(s), t, v = _PWcoeffs(18)), 1/sqrt(pi*t), rtol = 1e-4)
t = 5.0
@test isapprox(postwid(s -> 1/sqrt(s), t, v = _PWcoeffs(18)), 1/sqrt(pi*t), rtol = 1e-4)

setprecision(128)
t = 1.0
@test isapprox(postwid(s -> 1/sqrt(s), t, v = _PWcoeffs(36)), 1/sqrt(pi*t), rtol = 1e-12)
t = 5.0
@test isapprox(postwid(s -> 1/sqrt(s), t, v = _PWcoeffs(36)), 1/sqrt(pi*t), rtol = 1e-12)

# test across t array
t = 1.0:10.0

setprecision(53)
@test isapprox(postwid(s -> 1/sqrt(s), t, v = _PWcoeffs(18)), 1 ./ sqrt.(pi.*t), rtol = 1e-4)

setprecision(128)
@test isapprox(postwid(s -> 1/sqrt(s), t, v = _PWcoeffs(36)), 1 ./ sqrt.(pi.*t), rtol = 1e-12)

#### F(s) = 1/s^4, f(t) = t^3 / 6

# test for single time point
setprecision(53)
t = 1.0
@test isapprox(postwid(s -> 1/s^4, t, v = _PWcoeffs(18)), t^3/6, rtol = 1e-4)
t = 5.0
@test isapprox(postwid(s -> 1/s^4, t, v = _PWcoeffs(18)), t^3/6, rtol = 1e-4)

setprecision(128)
t = 1.0
@test isapprox(postwid(s -> 1/s^4, t, v = _PWcoeffs(36)), t^3/6, rtol = 1e-12)
t = 5.0
@test isapprox(postwid(s -> 1/s^4, t, v = _PWcoeffs(36)), t^3/6, rtol = 1e-12)

# test across t array
t = 1.0:10.0

setprecision(53)
@test isapprox(postwid(s -> 1/s^4, t, v = _PWcoeffs(18)), t.^3 ./6, rtol = 1e-4)

setprecision(128)
@test isapprox(postwid(s -> 1/s^4, t, v = _PWcoeffs(36)), t.^3 ./6, rtol = 1e-12)

#### F(s) = 1/(s+1), f(t) = exp(-t)
## for larger t you need a higher amount of coeffs depending on accuracy

# test for single time point
setprecision(53)
t = 1.0
@test isapprox(postwid(s -> 1/(s+1), t, v = _PWcoeffs(18)), exp(-t), rtol = 1e-4)
t = 5.0
@test isapprox(postwid(s -> 1/(s+1), t, v = _PWcoeffs(18)), exp(-t), rtol = 1e-3)

setprecision(128)
t = 1.0
@test isapprox(postwid(s -> 1/(s+1), t, v = _PWcoeffs(36)), exp(-t), rtol = 1e-12)
t = 5.0
@test isapprox(postwid(s -> 1/(s+1), t, v = _PWcoeffs(36)), exp(-t), rtol = 1e-8)

# test across t array
t = 1.0:10.0

setprecision(53)
@test isapprox(postwid(s -> 1/(s+1), t, v = _PWcoeffs(18)), exp.(-t), rtol = 1e-3)

setprecision(128)
@test isapprox(postwid(s -> 1/(s+1), t, v = _PWcoeffs(36)), exp.(-t), rtol = 1e-8)

end # module


text:
abstract type  AbsTransmission end
#初级传动
struct PrimTransmission <: AbsTransmission
    drivingGear::Int
    drivenGear::Int
    name::String
end
#主轴
struct DrivingTransmission <: AbsTransmission
    firstGear::Int
    secondGear::Int
    thirdGear::Int
    fourthGear::Int
    fifthGear::Int
    sixthGear::Int
    name::String
end
#副轴
struct DrivenTransmission <: AbsTransmission
    firstGear::Int
    secondGear::Int
    thirdGear::Int
    fourthGear::Int
    fifthGear::Int
    sixthGear::Int
    name::String
end
struct EngRatio
    ratio::Tuple
end
#大小飞轮
struct Flywheel <: AbsTransmission
    s_wheel::Int
    b_wheel::Int
    name::String
end
#发动机传动
struct Transmission <: AbsTransmission
    crankshaft::PrimTransmission
    drivingshaft::DrivingTransmission
    drivenshaft::DrivenTransmission
    flywheel::Flywheel

end

#轮周长
struct TireCircumference
    circumfeience::Float64
    name::String
end

#车速
struct VehicleSpeed
    engineName::String
    prim_ratio::Float64
    second_ratio::Vector{Float64}
    third_ratio::Float64
    circumfeience::Float64
    ca_speed::Float64
    speed::Vector{Float64}

    function VehicleSpeed(engineName,prim_ratio,second_ratio,third_ratio,circumfeience,ca_speed)
        n = length(second_ratio)
        sp = Vector{Float64}(undef,n)
        for i in 1:n
            sp[i] = ca_speed*circumfeience*60/((prim_ratio * second_ratio[i] * third_ratio)*1000)
        end
        new(engineName,prim_ratio,second_ratio,third_ratio,circumfeience,ca_speed,sp)
    end
end

max_stars_repo_name: powerzf/VehicleCalculate

max_issues_repo_name: powerzf/VehicleCalculate

max_forks_repo_name: powerzf/VehicleCalculate


text:
using Pkg
Pkg.activate(".")

using PlotlyJS
using Dash

 include("dash_utils.jl")
 include("ode.jl")

mutable struct ParametricModel
    name::String # Name for the model (e.g linear, quadratic, periodic, etc)
    formula::Function # Formula for the curvature function (depending on parameters)
    params::Vector{String} # Parameters in formula
end

mutable struct Example
    name::String
    model::ParametricModel
    param_values::Vector{Float64}
    smin::Float64
    smax::Float64
    resolution::Int64
end

# CurveData constructor for and Example
function CurveData(e::Example)
    κ_formula = e.model.formula(e.param_values...) # Evaluate model using parameter values
    CurveData(κ_formula, e.smin, e.smax)
end

# Parametric models for curvature function
linear = ParametricModel("Linear Model", (a,b)->"$a*s+$b", ["a","b"])
quadratic = ParametricModel("Quadratic Model", (a,b,c)->"$a*s^2+$b*s+$c", ["a","b","c"])
cubic = ParametricModel("Cubic Model", (a,b,c,d)->"$a*s^3+$b*s^2+$c*s+$d", ["a","b","c","d"])
periodic = ParametricModel("Periodic Model", (a,b,ω)->"$a*($b-cos($ω*s))", ["a","b","ω"])

# Examples
mpi = 3.1416

examples_list = [
    Example("Circle", linear, [0.0,1.0], -mpi, mpi, 250),
    Example("Double Spiral", linear, [1.0,0.0], -5,5, 300),
    Example("Hair", quadratic, [1.0,0.0,1.0],-5,5,500),
    Example("Periodic", periodic,[1.4142,1.0,1.0],-5*mpi,5*mpi,1000 ),
    Example("Cross", periodic,[7/4,1.0,1.0],-5*mpi,5*mpi,1000 ),
    Example("Snail", cubic,[1.0,0.0,0.0,2.0],-2.5,2.5,500 )
]
examples = Dict(e.name => e for e in examples_list)
default_example = examples["Circle"]
cd_default = CurveData(default_example)
sol = curve_from_curvature(CurveData(default_example),n = default_example.resolution)

plot_width = 500
plot_height = 500

function plot_solution(sol)
    x = [q[1] for q in sol]
    y = [q[2] for q in sol]
    min_x = minimum(x)
    max_x = maximum(x)
    min_y = minimum(y)
    max_y = maximum(y)
    Δx = max_x - min_x
    Δy = max_y - min_y
    mx = (min_x + max_x)/2
    my = (min_y + max_y)/2
    if Δy < Δx
        min_y = my - Δx/2
        max_y = my + Δx/2
    else
        min_x = mx - Δy/2
        max_x = mx + Δy/2
    end

    ϵ = 0.25
    layout = Layout(xaxis_title="x", yaxis_title="y",
    xaxis_range=[min_x-ϵ,max_x+ϵ],yaxis_range=[min_y-ϵ,max_y+ϵ], title="Curve",
    width=plot_width, height=plot_height)
    plot_data =[scatter(;x=x, y=y, mode="lines", name="Solution")]
    plot(plot_data, layout)
end

function plot_κ(cd::CurveData; sres=300)
    ss = LinRange(cd.smin,cd.smax, sres)
    κs = kappa.(ss,cd.κ_formula)
    layout = Layout(xaxis_title="arc-length", yaxis_title="curvature", title="Curvature function",
    xaxis_range=[cd.smin,cd.smax],
    width=plot_width, height=plot_height)
    plot_data =[scatter(;x=ss, y=κs, mode="lines")]
    plot(plot_data, layout)
end

mathjax = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js?config=TeX-MML-AM_CHTML"
app = dash(external_scripts=[mathjax])
app.title = "DashCurvature"

basics = load_markdown("basics.md")
details = load_markdown("details.md")

info = html_div([
    dcc_tabs(id="tabs", value="tab-1", children=[
        dcc_tab(label="Basics", value="tab-1"),
        dcc_tab(label="Details", value="tab-2"),
    ]),
    html_div(id="tabs-content")
])

sol_plot = plot_solution(sol)
κ_plot = plot_κ(CurveData(default_example))

curve_graph = dcc_loading([dcc_graph(id="curve_graph",figure=sol_plot)])
κ_graph = dcc_loading([dcc_graph(id="κ_graph",figure=κ_plot)])
main_graph = row_div([κ_graph,curve_graph])

inputs =html_div([
    html_center(html_h4("Curvature input")),
    html_break(1),
    html_center(html_div([
    input_field("kappa",raw"$\kappa=\kappa(s)$", value = cd_default.κ_formula),
    input_field("smin", "length min", value = string(cd_default.smin)),
    input_field("smax", "length max", value = string(cd_default.smax)),
    input_field("res", "resolution", value = string(300)),
    html_button("Solve", id="solve-button", n_clicks=0),
    html_break(1)
    ]))])

dropdown = html_div(style=Dict("width"=>"250px"), dcc_dropdown(id="examples", options=[(label = e.name, value = e.name) for e in values(examples)], value=default_example.name))

function to_info(e::Example)
    str = "Parameter's values: "
    for (i, p) in enumerate(e.model.params)
        val = e.param_values[i]
        str = str*p*"= $val, "
    end
    str = str[1:end-2]
    ["Parametric model: s->"*e.model.formula(e.model.params...),html_br(), str]
end

examples_control = html_div([
    dropdown,
    html_break(1),
    html_div(id="model_info", to_info(default_example)),
    html_br(),
    ])

controls = row_div([inputs, examples_control])

# App Layout
app.layout = html_div(
    [
        html_center(html_h1("Curve from Curvature")),
        html_h5(),
        basics,
        html_break(1),
        main_graph,
        html_break(1),
        controls,
        html_break(1),
        details
    ]
)

callback!(
    app,
    Output("model_info","children"),
    Output("kappa", "value"),
    Output("smin","value"),
    Output("smax","value"),
    Output("res","value"),
    Output("solve-button","n_clicks"),
    Input("examples", "value")) do value
    e = examples[value]
    formula = e.model.formula(e.param_values...)
    cd = CurveData(formula,e.smin,e.smax)
    κ_plot = plot_κ(cd)
    return (to_info(e),formula, string(e.smin), string(e.smax), string(e.resolution),2)
end

callback!(
    app,
    Output("curve_graph", "figure"),
    Output("κ_graph", "figure"),
    Input("solve-button", "n_clicks"),
    State("kappa", "value"),
    State("smin", "value"),
    State("smax", "value"),
    State("res", "value")
) do clicks, kappa, smin, smax, res
    sminf = parse(Float64,smin)
    smaxf = parse(Float64,smax)
    resf = parse(Int64, res)
    cd = CurveData(kappa,sminf,smaxf)
    sol = curve_from_curvature(cd,n=resf)
    sol_plot = plot_solution(sol)
    κ_plot = plot_κ(cd)
    return sol_plot, κ_plot
end

run_server(app, "0.0.0.0", debug=true)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://twitter.com/yujitach/status/1424030835771023363

# %%
VERSION

# %%
]st

# %%
Threads.nthreads()

# %%
"""
    Original

* https://gist.github.com/yujitach/c30d7a174bbc3d3d3e40a3c0f9f9d47f
* TAB を "    " で置換
"""
module Original

using LinearAlgebra,LinearMaps
import Arpack

const L=20

diag_ = zeros(Float64,2^L)

function prepareDiag(diag)
    for state = 1 : 2^L
        for i = 1 : L
            j = i==L ? 1 : i+1
            diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
end

function Hfunc!(C,B,diag)
    for state = 1 : 2^L
        C[state] = diag[state] * B[state]
    end
    for state = 1 : 2^L
        for i = 1 : L
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            if newstate==0
                newstate = 2^L
            end
            C[newstate] -= B[state]
        end
    end
end

println("preparing...")
prepareDiag(diag_)

println("computing the lowest eigenvalue...")
H=LinearMap((C,B)->Hfunc!(C,B,diag_),2^L,ismutating=true,issymmetric=true,isposdef=false)
@time e,v = Arpack.eigs(H,nev=1,which=:SR)
@time e,v = Arpack.eigs(H,nev=1,which=:SR)

println("obtained:")
println(e[1])

println("theoretical:")
println(-2sum([ abs(sin((n-1/2) * pi/L)) for n in 1 : L]))

end;

# %%
"""
    Rev0

* This revision is almost equivalent to the original.
* Stop using constants.
* Always pass global variables to functions as arguments.
* Swap the order of the for loop.
* Add @inbounds macro.
* Revise sum([f(x) for x in X]) to sum(f(x) for x in X).
"""
module Rev0

using LinearAlgebra, LinearMaps
import Arpack

function prepareDiag(L)
    diag = zeros(2^L)
    for state = 1:2^L
        for i = 1:L
            j = i==L ? 1 : i+1
            @inbounds diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
    diag
end

function Hfunc!(C, B, diag, L)
    for state = 1:2^L
        @inbounds C[state] = diag[state] * B[state]
    end
    for i = 1:L
        for state = 1:2^L
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            if newstate == 0
                newstate = 2^L
            end
            @inbounds C[newstate] -= B[state]
        end
    end
end
prepareHfunc!(diag, L) = (C, B) -> Hfunc!(C, B, diag, L)

L = 20

println("preparing...")
diag_ = prepareDiag(L)

println("computing the lowest eigenvalue...")
H = LinearMap(prepareHfunc!(diag_, L), 2^L, ismutating=true, issymmetric=true, isposdef=false)
@time e, v = Arpack.eigs(H, nev=1, which=:SR)
@time e, v = Arpack.eigs(H, nev=1, which=:SR)

println("obtained:")
println(e[1])

println("theoretical:")
println(-2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L))

end;

# %%
"""
    Rev1

* Use Threads.@threads macro.
"""
module Rev1

using LinearAlgebra, LinearMaps
import Arpack

function prepareDiag(L)
    diag = zeros(2^L)
    for state = 1:2^L
        for i = 1:L
            j = i==L ? 1 : i+1
            @inbounds diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
    diag
end

function Hfunc!(C, B, diag, L)
    Threads.@threads for state = 1:2^L
        @inbounds C[state] = diag[state] * B[state]
    end
    for i = 1:L
        Threads.@threads for state = 1:2^L
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            if newstate == 0
                newstate = 2^L
            end
            @inbounds C[newstate] -= B[state]
        end
    end
end
prepareHfunc!(diag, L) = (C, B) -> Hfunc!(C, B, diag, L)

L = 20

println("preparing...")
diag_ = prepareDiag(L)

println("computing the lowest eigenvalue...")
H = LinearMap(prepareHfunc!(diag_, L), 2^L, ismutating=true, issymmetric=true, isposdef=false)
@time e, v = Arpack.eigs(H, nev=1, which=:SR)
@time e, v = Arpack.eigs(H, nev=1, which=:SR)

println("obtained:")
println(e[1])

println("theoretical:")
println(-2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L))

end;

# %%
"""
    Rev2

* Use LoopVectorization.@tturbo macro.
"""
module Rev2

using LinearAlgebra, LinearMaps
import Arpack
using LoopVectorization

function prepareDiag(L)
    diag = zeros(2^L)
    for state = 1:2^L
        for i = 1:L
            j = i==L ? 1 : i+1
            @inbounds diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
    diag
end

function Hfunc!(C, B, diag, L)
    N = length(diag)
    @tturbo for state = 1:N
        C[state] = diag[state] * B[state]
    end
    for i = 1:L
        @tturbo for state = 1:N
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            c = newstate == 0
            newstate = !c*newstate + c*N # remove if statement
            C[newstate] -= B[state]
        end
    end
end
prepareHfunc!(diag, L) = (C, B) -> Hfunc!(C, B, diag, L)

L = 20

println("preparing...")
diag_ = prepareDiag(L)

println("computing the lowest eigenvalue...")
H = LinearMap(prepareHfunc!(diag_, L), 2^L, ismutating=true, issymmetric=true, isposdef=false)
@time e, v = Arpack.eigs(H, nev=1, which=:SR)
@time e, v = Arpack.eigs(H, nev=1, which=:SR)

println("obtained:")
println(e[1])

println("theoretical:")
println(-2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L))

end;

# %%
using BenchmarkTools
using LinearAlgebra
using Arpack: Arpack

H = Original.H
H0 = Rev0.H
H1 = Rev1.H
H2 = Rev2.H

B, C = similar(Original.diag_), similar(Original.diag_)
println("Hamiltonian bench")
print("  Original:                ")
@btime mul!($C, $H, $B)
print("  Rev0 (almost original):  ")
@btime mul!($C, $H0, $B)
print("  Rev1 (Threads.@threads): ")
@btime mul!($C, $H1, $B)
print("  Rev2 (LoopVectorization):")
@btime mul!($C, $H2, $B);

# %%
println("Arpack.eigs bench")
print("  Original:                ")
@btime e, v = Arpack.eigs($H, nev=1, which=:SR)
print("  Rev0 (almost original):  ")
@btime e, v = Arpack.eigs($H0, nev=1, which=:SR)
print("  Rev1 (Threads.@threads): ")
@btime e, v = Arpack.eigs($H1, nev=1, which=:SR)
print("  Rev2 (LoopVectorization):")
@btime e, v = Arpack.eigs($H2, nev=1, which=:SR);

# %%
@show a = 2
@show foo = x -> a*x
@show foo(3)
@show a = 10
@show foo(3);

# %%
@show a = 2
@show makebar(a) = x -> a*x
@show bar = makebar(a)
@show bar(3)
@show a = 10
@show bar(3);

# %%

max_stars_repo_path: 0015/Arpack/1D Ising example - Arpack, LinearMaps, Threads, LoopVectorization.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# graph represented by a light-weight edge list

type GenericEdgeList{V,E,VList,EList} <: AbstractGraph{V,E}
    is_directed::Bool
    vertices::VList
    edges::EList
end

@graph_implements GenericEdgeList vertex_list edge_list vertex_map edge_map

typealias SimpleEdgeList{E} GenericEdgeList{Int,E,UnitRange{Int},Vector{E}}
typealias EdgeList{V,E} GenericEdgeList{V,E,Vector{V},Vector{E}}

# construction

simple_edgelist{E}(nv::Integer, edges::Vector{E}; is_directed::Bool=true) =
    SimpleEdgeList{E}(is_directed, intrange(nv), edges)

edgelist{V,E}(vertices::Vector{V}, edges::Vector{E}; is_directed::Bool=true) =
    EdgeList{V,E}(is_directed, vertices, edges)

# required interface

is_directed(g::GenericEdgeList) = g.is_directed

num_vertices(g::GenericEdgeList) = length(g.vertices)
vertices(g::GenericEdgeList) = g.vertices

num_edges(g::GenericEdgeList) = length(g.edges)
edges(g::GenericEdgeList) = g.edges
edge_index(e, g::GenericEdgeList) = edge_index(e)

# mutation

add_vertex!{V}(g::GenericEdgeList{V}, v::V) = (push!(g.vertices, v); v)
add_vertex!{V}(g::GenericEdgeList{V}, x) = add_vertex!(g, make_vertex(g, x))

add_edge!{V,E}(g::GenericEdgeList{V,E}, e::E) = (push!(g.edges, e); e)
add_edge!{V,E}(g::GenericEdgeList{V,E}, u::V, v::V) = add_edge!(g, make_edge(g, u, v))


text:
import .libffmpeg: AVDictionary
import Base: getindex, setindex!, iterate, length, empty!

mutable struct AVDict <: AbstractDict{String, String}
    ref_ptr_dict::Ref{Ptr{AVDictionary}}
end

function AVDict()
    d = AVDict(Ref{Ptr{AVDictionary}}(C_NULL))
    finalizer(empty!, d)
end

function AVDict(ps::Pair...)
    d = AVDict()
    for (key, value) in ps
        d[key] = value
    end
    return d
end

function AVDict(o::AbstractDict)
    d = AVDict()
    for (key, value) in o
        d[key] = value
    end
    return d
end

Base.empty!(d::AVDict) = libffmpeg.av_dict_free(d.ref_ptr_dict)

Base.convert(::Type{Ref{Ptr{AVDictionary}}}, d::AVDict) = d.ref_ptr_dict

function setindex!(d::AVDict, value, key)
    libffmpeg.av_dict_set(d.ref_ptr_dict, string(key), string(value), 0)
    nothing
end

function getindex(d::AVDict, key::AbstractString)
    pItem = libffmpeg.av_dict_get(d.ref_ptr_dict[], key, C_NULL, 0)
    item = unsafe_load(pItem)
    value = unsafe_string(item.value)
    return value
end

function iterate(d::AVDict, state=C_NULL)
    pItem = libffmpeg.av_dict_get(d.ref_ptr_dict[], "", state, libffmpeg.AV_DICT_IGNORE_SUFFIX)
    pItem == C_NULL && return nothing

    item = unsafe_load(pItem)
    key = unsafe_string(item.key)
    value = unsafe_string(item.value)

    return (key => value, pItem)
end

function length(d::AVDict)
    return libffmpeg.av_dict_count(d.ref_ptr_dict[])
end


text:
# using Pkg; for p in ("Knet","ArgParse"); haskey(Pkg.installed(),p) || Pkg.add(p); end

"""

julia vgg.jl image-file-or-url

This example implements the VGG model from `Very Deep Convolutional
Networks for Large-Scale Image Recognition', Karen Simonyan and Andrew
Zisserman, arXiv technical report 1409.1556, 2014. This example works
for D and E models currently. VGG-D is the default model if you do not
specify any model.

* Paper url: https://arxiv.org/abs/1409.1556
* Project page: http://www.robots.ox.ac.uk/~vgg/research/very_deep
* MatConvNet weights used here: http://www.vlfeat.org/matconvnet/pretrained

"""
module VGG
using Knet,CUDA,ArgParse
include(Knet.dir("data","imagenet.jl"))

const imgurl = "https://github.com/BVLC/caffe/raw/master/examples/images/cat.jpg"
const vggurl = "http://www.vlfeat.org/matconvnet/models/imagenet-vgg-verydeep-16.mat"
const LAYER_TYPES = ["conv", "relu", "pool", "fc", "prob"]

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="vgg.jl (c) Deniz Yuret, İlker Kesen, 2016. Classifying images with the VGG model from http://www.robots.ox.ac.uk/~vgg/research/very_deep."
    # s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("image"; default=imgurl; help="Image file or URL.")
        ("--model"; default="imagenet-vgg-verydeep-16"; help="Model name")
        ("--top"; default=5; arg_type=Int; help="Display the top N classes")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array and float type to use")
    end
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end
    println(s.description)
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    atype = eval(Meta.parse(o[:atype]))

    global _vggcache
    if !@isdefined(_vggcache); _vggcache=Dict(); end
    cachekey = (o[:model],o[:atype])
    if !haskey(_vggcache,cachekey)
        vgg = matconvnet(o[:model])
        params = get_params(vgg, atype)
        convnet = get_convnet(params...)
        description = vgg["meta"]["classes"]["description"]
        averageImage = convert(Array{Float32},vgg["meta"]["normalization"]["averageImage"])
        _vggcache[cachekey] = vgg, params, convnet, description, averageImage
    else
        vgg, params, convnet, description, averageImage = _vggcache[cachekey]
    end

    image = imgdata(o[:image], averageImage)
    image = convert(atype, image)
    @info("Classifying")
    @time y1 = convnet(image)
    z1 = vec(Array(y1))
    s1 = sortperm(z1,rev=true)
    p1 = exp.(logp(z1))
    display(hcat(p1[s1[1:o[:top]]], description[s1[1:o[:top]]]))
    println()
end

# This procedure makes pretrained MatConvNet VGG parameters convenient for Knet
# Also, if you want to extract features, specify the last layer you want to use
function get_params(CNN, atype; last_layer="prob")
    layers = CNN["layers"]
    weights, operations, derivatives = [], [], []

    for l in layers
        get_layer_type(x) = startswith(l["name"], x)
        operation = filter(x -> get_layer_type(x), LAYER_TYPES)[1]
        push!(operations, operation)
        push!(derivatives, haskey(l, "weights") && length(l["weights"]) != 0)

        if derivatives[end]
            w = copy(l["weights"])
            if operation == "conv"
                w[2] = reshape(w[2], (1,1,length(w[2]),1))
            elseif operation == "fc"
                w[1] = transpose(mat(w[1]))
            end
            push!(weights, w)
        end

        last_layer != nothing && get_layer_type(last_layer) && break
    end

    map(w -> map(wi->convert(atype,wi), w), weights), operations, derivatives
end

# get convolutional network by interpreting parameters
function get_convnet(weights, operations, derivatives)
    function convnet(xs)
        i, j = 1, 1
        num_weights, num_operations = length(weights), length(operations)
        while i <= num_operations && j <= num_weights
            if derivatives[i]
                xs = forw(xs, operations[i], weights[j])
                j += 1
            else
                xs = forw(xs, operations[i])
            end

            i += 1
        end
        convert(Array{Float32}, xs)
    end
end

# convolutional network operations
convx(x,w) = conv4(w[1], x; padding=1, mode=1) .+ w[2]
relux(x) = relu.(x)
poolx = pool
probx(x) = x
fcx(x,w) = w[1] * mat(x) .+ w[2]
tofunc(op) = eval(Meta.parse(string(op, "x")))
forw(x,op) = tofunc(op)(x)
forw(x,op,w) = tofunc(op)(x,w)

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia vgg.jl cat.jpg
# julia> VGG.main("cat.jpg")
PROGRAM_FILE=="vgg.jl" && main(ARGS)

end # module


text:
@testset "update_y_l" begin
#test update_y_l.jl
Random.seed!(123)

TF = Float64
TI = Int64

#without BLAS
x=randn(100)
p=2
i=10
Blas_active = false

y=Vector{Vector{Float64}}(undef, p)
y[1]=randn(51)
y[2]=randn(100)
y_old=Vector{Vector{Float64}}(undef, p)
y_old[1]=randn(51)
y_old[2]=randn(100)
l_old=Vector{Vector{Float64}}(undef, p)
l_old[1]=randn(51)
l_old[2]=randn(100)
l=Vector{Vector{Float64}}(undef, p)
l[1]=randn(51)
l[2]=randn(100)

rho=[1.234;10.23432]
gamma=Vector{Float64}(undef, p)
gamma[1]=1.0
gamma[2]=1.345

prox=Vector{Any}(undef, p)
prox[1] = x -> 1.0.*x
m=randn(100)
prox[2] = x -> prox_l2s!(x,rho[2],m)
TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(undef, p)
TD_OP[1] = sparse(1.0*I, 51,100)*2.0
TD_OP[2] = sparse(1.0*I, 100, 100)
maxit=39
log_PSDMM = log_type_PARSDMM(zeros(maxit,p-1),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit,p),
zeros(maxit,p),zeros(maxit),zeros(maxit),TimerOutput());

P_sub = Vector{Any}(undef, 1)
P_sub[1] = prox[1]
counter=12
x_hat=Vector{Vector{Float64}}(undef, p)
x_hat[1]=randn(51)
x_hat[2]=randn(100)
r_pri=Vector{Vector{Float64}}(undef, p)
r_pri[1]=randn(51)
r_pri[2]=randn(100)
s=Vector{Vector{Float64}}(undef, p)
s[1]=randn(51)
s[2]=randn(100)

#construct rerefrence solution
x_hat2=deepcopy(x_hat)
y2=deepcopy(y)
l2=deepcopy(l)
x2=deepcopy(x)
for i=1:p
  x_hat2[i]=gamma[i]*TD_OP[i]*x2 + (1-gamma[i])*y2[i]
  y2[i]=prox[i](x_hat2[i]-l2[i]./rho[i])
  l2[i].=l2[i].+rho[i].*(y2[i].-x_hat2[i])
end

(y,l,r_pri,s,log_PSDMM,counter,y_old,l_old)=update_y_l(x,p,i,Blas_active,
  y,y_old,l,l_old,rho,gamma,prox,TD_OP,log_PSDMM,P_sub,counter,x_hat,r_pri,s)

for i=1:p #don't test relative norm here, because l may be zeros/eps() valued
  @test norm(y[i]-y2[i]) <=1e-14
  @test norm(l[i]-l2[i]) <=1e-14
  @test norm(s[i]-TD_OP[i]*x) <= 1e-14
  @test norm(r_pri[i]-(-TD_OP[i]*x+y[i])) <=1e-14
end

#check that gamma and rho do not change
@test gamma[2]==1.345
@test gamma[1]==1.0
@test rho[1]==1.234
@test rho[2]==10.23432

#verify that y_old, l_old are not the same as current y and l
for i=1:p #don't test relative norm here, because l may be zeros/eps() valued
  @test norm(y[i]-y_old[i]) >10*eps()
  @test norm(l[i]-l_old[i]) >10*eps()
end

#now with BLAS
x=randn(100)
p=2
i=10
Blas_active = true

y=Vector{Vector{Float64}}(undef, p)
y[1]=randn(51)
y[2]=randn(100)
y_old=Vector{Vector{Float64}}(undef, p)
y_old[1]=randn(51)
y_old[2]=randn(100)
l_old=Vector{Vector{Float64}}(undef, p)
l_old[1]=randn(51)
l_old[2]=randn(100)
l=Vector{Vector{Float64}}(undef, p)
l[1]=randn(51)
l[2]=randn(100)

rho=[1.234;10.23432]
gamma=Vector{Float64}(undef, p)
gamma[1]=1.0
gamma[2]=1.345

prox=Vector{Any}(undef, p)
prox[1] = x -> 1.0.*x
m=randn(100)
prox[2] = x -> prox_l2s!(x,rho[2],m)
TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(undef, p)

TD_OP[1] = sparse(1.0*I, 51,100)*2.0
TD_OP[2] = sparse(1.0*I, 100, 100)
maxit=39
log_PSDMM = log_type_PARSDMM(zeros(maxit,p-1),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),TimerOutput());
P_sub = Vector{Any}(undef, 1)
P_sub[1] = prox[1]
counter=12
x_hat=Vector{Vector{Float64}}(undef, p)
x_hat[1]=randn(51)
x_hat[2]=randn(100)
r_pri=Vector{Vector{Float64}}(undef, p)
r_pri[1]=randn(51)
r_pri[2]=randn(100)
s=Vector{Vector{Float64}}(undef, p)
s[1]=randn(51)
s[2]=randn(100)

#construct rerefrence solution
x_hat2=deepcopy(x_hat)
y2=deepcopy(y)
l2=deepcopy(l)
x2=deepcopy(x)
for i=1:p
  x_hat2[i]=gamma[i]*TD_OP[i]*x2 + (1-gamma[i])*y2[i]
  y2[i]=prox[i](x_hat2[i]-l2[i]./rho[i])
  l2[i].=l2[i].+rho[i].*(y2[i].-x_hat2[i])
end

(y,l,r_pri,s,log_PSDMM,counter,y_old,l_old)=update_y_l(x,p,i,Blas_active,
y,y_old,l,l_old,rho,gamma,prox,TD_OP,log_PSDMM,P_sub,counter,x_hat,r_pri,s)

for i=1:p #don't test relative norm here, because l may be zeros/eps() valued
  @test norm(y[i]-y2[i]) <=1e-14
  @test norm(l[i]-l2[i]) <=1e-14
  @test norm(s[i]-TD_OP[i]*x) <= 1e-14
  @test norm(r_pri[i]-(-TD_OP[i]*x+y[i])) <=1e-14
end

#check that gamma and rho do not change
@test gamma[2]==1.345
@test gamma[1]==1.0
@test rho[1]==1.234
@test rho[2]==10.23432

#verify that y_old, l_old are not the same as current y and l
for i=1:p #don't test relative norm here, because l may be zeros/eps() valued
  @test norm(y[i]-y_old[i]) >10*eps()
  @test norm(l[i]-l_old[i]) >10*eps()
end

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
"""
    geometric_rng(p, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Geometric(`p`) distribution. Optionally you can set a specific seed.

# Notes

The Geometric distribution is given:

``f(x,p) = (1-p)^{x-1}p) \\quad x = 1,2,3,…``

where ``0 ≤ p ≤ 1``

# Examples

```julia-repl
julia> geometric_rng(.8)
1-element Vector{Int64}:
 1

julia> geometric_rng(.8, 5)
5-element Vector{Int64}:
 2
 3
 1
 1
 1

julia> geometric_rng(.8, (2,2), seed=45)
2×2 Matrix{Int64}:
 1  1
 1  1

```

# References

D.P. Kroese, T. Taimre, Z.I. Botev. Handbook of Monte Carlo Methods.
  Wiley Series in Probability and Statistics, John Wiley & Sons, New York, 2011.

C. Alexopoulos, D. Goldsman. Random variate generation. 2020.
"""
function geometric_rng(p::Real, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    check_p(p)
    c = log(1 - p)
    U = get_std_uniform(shape, seed=seed)
    X = zeros(Int, shape)
    X = floor.(Int, log.(U) ./ c)  # here sub U for 1-U is fine
    return X
end


text:
"""
SPARSETOCSR
-----------
    convert sparse matrix into compressed row storage arrays
    and returns the row pointer (rp), column index (ci) and value index (ai) arrays
    of a compressed sparse representation of the matrix A (or its triplet format storage)

Functions
---------
- (rp,ci,ai,m) = sparse_to_csr{T}(A::SparseMatrixCSC{T,Int64})
- (rp,ci,ai,m) = sparse_to_csr{T}(nzi::Array{Int64,1},nzj::Array{Int64,1},nzv::Array{T,1})

Example
-------
~~~
i = [1;2;3]
j = [3;4;4]
v = [8;9;10]
(rp,ci,ai,m) = sparse_to_csr(i,j,v)
~~~
"""
function sparse_to_csr(A::SparseMatrixCSC{T,Int64}) where T
    At = copy(A');
    return (At.colptr,At.rowval,At.nzval,At.m);
end

function sparse_to_csr(nzi::Array{Int64,1},nzj::Array{Int64,1},
                         nzv::Array{T,1}) where T
    At = sparse(nzj,nzi,nzv)
    return (At.colptr,At.rowval,At.nzval,At.m)
end

max_stars_repo_name: nassarhuda/MatrixNetworks

max_issues_repo_name: nassarhuda/MatrixNetworks

max_forks_repo_name: nassarhuda/MatrixNetworks


text:
using Merlin
using HDF5

const wordembeds_file = "wordembeds_nyt100.h5"

function train()
    words = h5read(wordembeds_file, "s")
    worddict = Dict(words[i] => i for i=1:length(words))
    chardict = Dict{Char,Int}(' ' => 1)
    tagdict = Dict{String,Int}()
    #traindoc = UD_English.traindata()
    #testdoc = UD_English.testdata()
    train_w, train_c, train_t = readdata(".data/wsj_00-18.conll", worddict, chardict, tagdict)
    test_w, test_c, test_t = readdata(".data/wsj_22-24.conll", worddict, chardict, tagdict)
    info("# Training sentences:\t$(length(train_w))")
    info("# Testing sentences:\t$(length(test_w))")
    info("# Words:\t$(length(worddict))")
    info("# Chars:\t$(length(chardict))")
    info("# Tags:\t$(length(tagdict))")

    wordembeds = h5read(wordembeds_file, "v")
    T = eltype(wordembeds)
    charembeds = rand(T, 10, length(chardict))
    charembeds = charembeds * T(0.002) - T(0.001)
    nn = Model(wordembeds, charembeds, length(tagdict))
    opt = SGD(momentum=0.95, nesterov=true)

    batchsize = 8
    train_data = makebatch(batchsize, train_w, train_c, train_t)
    train_data = collect(zip(train_data...))
    test_data = makebatch(100, test_w, test_c, test_t)
    test_data = collect(zip(test_data...))

    for epoch = 1:10
        println("Epoch: $epoch")
        opt.rate = 0.001 / batchsize / epoch
        #opt.rate = 0.0075 / (1 + 0.05*(epoch-1))

        loss = 0
        for i in randperm(length(train_data))
            w, c, t = train_data[i]
            y = nn(w, c)
            y = softmax_crossentropy(t, y)
            loss += sum(y.data)
            gradient!(y)
            update!(nn.fw, opt)
            update!(nn.fc, opt)
            update!(nn.fs, opt)
        end
        loss /= batchsize * length(train_data)
        # loss = minimize!(train_f, opt, train_data)
        println("Loss: $loss")

        # test
        ys = Int[]
        zs = Int[]
        for (w,c,t) in test_data
            append!(ys, t.data)
            y = nn(w, c)
            z = vec(argmax(y.data,1))
            append!(zs, z)
        end
        length(ys) == length(zs) || throw("Length mismatch.")

        acc = mean(i -> ys[i] == zs[i] ? 1.0 : 0.0, 1:length(ys))
        acc = round(acc, 5)
        println("Test acc.: $acc")
        println()
    end
    # Merlin.save("postagger.h5", nn)
end

include("data.jl")
include("model.jl")
train()

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Shapes (Modia3D/Solids/_module.jl)
#

# This algorithm accords with [D. Eberly (2002), Polyhedral Mass Properties (Revisited)]
# https://www.geometrictools.com/Documentation/PolyhedralMassProperties.pdf
# it returns volume, centroid and inertia of a polyhedral shape.
# it input values are:
#   vertices          ... all nodes of an object
#   triangle_indices  ... indices of faces, which span the surface.
#                         All faces should be specified in right-handed/counter-clockwise order!

neps = sqrt(eps())

function helpingFunc(w0,w1,w2)
    tmp0 = w0 + w1
    f1   = tmp0 + w2
    tmp1 = w0.*w0
    tmp2 = tmp1 + w1.*tmp0
    f2   = tmp2 + w2.*f1
    f3   = w0.*tmp1 + w1.*tmp2 + w2.*f2
    g0   = f2 + w0.*(f1+w0)
    g1   = f2 + w1.*(f1+w1)
    g2   = f2 + w2.*(f1+w2)
    return (f1,f2,f3,g0,g1,g2)
end

function computeMassProperties(vertices, triangle_indices; bodyCoords::Bool=false)
    if isempty(triangle_indices)
        # (volume, centroid, inertia)
        return (0.0, zeros(3), zeros(9))
    else
    weights = [1/6, 1/24, 1/24, 1/24, 1/60, 1/60, 1/60, 1/120, 1/120, 1/120]
    integral = zeros(10)
    volume     = 0.0
    centroid       = zeros(3)
    inertia  = SMatrix{3,3,Float64,9}(zeros(9))
    for i=eachindex(triangle_indices)
        # get vertices of i-th triangle
        v1 = vertices[triangle_indices[i][1]]
        v2 = vertices[triangle_indices[i][2]]
        v3 = vertices[triangle_indices[i][3]]

        # get edges and cross products of edges
        v1v2 = v2 - v1
        v1v3 = v3 - v1
        n    = cross(v1v2, v1v3) # )

        (f1,f2,f3,g0,g1,g2) = helpingFunc(v1, v2, v3)

        integral[1] += n[1] * f1[1]
        integral[2] += n[1] * f2[1]
        integral[3] += n[2] * f2[2]
        integral[4] += n[3] * f2[3]
        integral[5] += n[1] * f3[1]
        integral[6] += n[2] * f3[2]
        integral[7] += n[3] * f3[3]
        integral[8] += n[1] * (v1[2]*g0[1] + v2[2]*g1[1] + v3[2]*g2[1])
        integral[9] += n[2] * (v1[3]*g0[2] + v2[3]*g1[2] + v3[3]*g2[2])
        integral[10] += n[3] * (v1[1]*g0[3] + v2[1]*g1[3] + v3[1]*g2[3])
    end

    for i=1:10
        integral[i] *= weights[i]
    end

    volume = integral[1]
    if volume > neps
        # center of volume
        centroid = [integral[2], integral[3], integral[4]]/volume

        # inertia tensor
        inertia_xx = integral[6] + integral[7]
        inertia_yy = integral[5] + integral[7]
        inertia_zz = integral[5] + integral[6]
        inertia_xy = -integral[8]
        inertia_yz = -integral[9]
        inertia_xz = -integral[10]

        if bodyCoords
            inertia_xx -= volume*(centroid[2]^2 + centroid[3]^2)
            inertia_yy -= volume*(centroid[3]^2 + centroid[1]^2)
            inertia_zz -= volume*(centroid[1]^2 + centroid[2]^2)
            inertia_xy += volume*centroid[1]*centroid[2]
            inertia_yz += volume*centroid[2]*centroid[3]
            inertia_xz += volume*centroid[3]*centroid[1]
        end
        inertia = SMatrix{3,3,Float64,9}([inertia_xx inertia_xy inertia_xz; inertia_xy inertia_yy inertia_yz; inertia_xz inertia_yz inertia_zz])
    end
    return (volume, centroid, inertia)
end; end


text:
using Documenter, TensorMatrices_lemon

makedocs()
# makedocs(modules=[TensorMatrices_lemon],
# 		 doctest=true)
#
 deploydocs(
			deps = Deps.pip("mkdocs", "python-markdown-math"),
 		  repo = "github.com/yoyoyoju/TensorMatrices_lemon.jl.git",
 		  julia = "0.5",
 		  osname = "osx"
		  )


text:
"""
Stores the results from running a judgement, see [`judge`](@ref).

The following (unexported) methods are defined on a `BenchmarkJudgement` (written below as `judgement`):

* `target_result(judgement)::BenchmarkResults` - the [`BenchmarkResults`](@ref) of the `target`.
* `baseline_result(judgement)::BenchmarkResults` -  the [`BenchmarkResults`](@ref) of the `baseline`.
* `benchmarkgroup(judgement)::BenchmarkGroup` - a [`BenchmarkGroup`](https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/doc/manual.md#the-benchmarkgroup-type)
   contaning the estimated results

A `BenchmarkJudgement` can be exported to markdown using the function [`export_markdown`](@ref).

See also [`BenchmarkResults`](@ref)
"""
struct BenchmarkJudgement
    target_results::BenchmarkResults
    baseline_results::BenchmarkResults
    benchmarkgroup::BenchmarkGroup
end

target_result(judgement::BenchmarkJudgement) = judgement.target_results
baseline_result(judgement::BenchmarkJudgement) = judgement.baseline_results
benchmarkgroup(judgement::BenchmarkJudgement) = judgement.benchmarkgroup

BenchmarkTools.isinvariant(f, judgement::BenchmarkJudgement) = BenchmarkTools.isinvariant(f, benchmarkgroup(judgement))
BenchmarkTools.isinvariant(judgement::BenchmarkJudgement) = BenchmarkTools.isinvariant(benchmarkgroup(judgement))

BenchmarkTools.isregression(f, judgement::BenchmarkJudgement) = BenchmarkTools.isregression(f, benchmarkgroup(judgement))
BenchmarkTools.isregression(judgement::BenchmarkJudgement) = BenchmarkTools.isregression(benchmarkgroup(judgement))

BenchmarkTools.isimprovement(f, judgement::BenchmarkJudgement) = BenchmarkTools.isimprovement(f, benchmarkgroup(judgement))
BenchmarkTools.isimprovement(judgement::BenchmarkJudgement) = BenchmarkTools.isimprovement(benchmarkgroup(judgement))

BenchmarkTools.invariants(f, judgement::BenchmarkJudgement) = BenchmarkTools.invariants(f, benchmarkgroup(judgement))
BenchmarkTools.invariants(judgement::BenchmarkJudgement) = BenchmarkTools.invariants(benchmarkgroup(judgement))

BenchmarkTools.regressions(f, judgement::BenchmarkJudgement) = BenchmarkTools.regressions(f, benchmarkgroup(judgement))
BenchmarkTools.regressions(judgement::BenchmarkJudgement) = BenchmarkTools.regressions(benchmarkgroup(judgement))

BenchmarkTools.improvements(f, judgement::BenchmarkJudgement) = BenchmarkTools.improvements(f, benchmarkgroup(judgement))
BenchmarkTools.improvements(judgement::BenchmarkJudgement) = BenchmarkTools.improvements(benchmarkgroup(judgement))

function Base.show(io::IO, judgement::BenchmarkJudgement)
    target, base = judgement.target_results, judgement.baseline_results
    print(io, "Benchmarkjudgement (target / baseline):\n")
    println(io, "    Package: ", target.name)
    println(io, "    Dates: ", Dates.format(target.date,  "d u Y - H:M"), " / ",
                               Dates.format(base.date, "d u Y - H:M"))
    println(io, "    Package commits: ", target.commit[1:min(length(target.commit), 6)], " / ",
                                        base.commit[1:min(length(base.commit), 6)])
    println(io, "    Julia commits: ", target.julia_commit[1:6], " / ",
                                       base.julia_commit[1:6])
end

function export_markdown(file::String, results::BenchmarkJudgement; kwargs...)
    open(file, "w") do f
        export_markdown(f, results; kwargs...)
    end
end

function export_markdown(io::IO, judgement::BenchmarkJudgement; export_invariants::Bool = false)
    target, baseline = judgement.target_results, judgement.baseline_results
    function env_strs(res)
        return if isempty(benchmarkconfig(res).env)
            "None"
        else
            join(String[string("`", k, " => ", v, "`") for (k, v) in benchmarkconfig(res).env], " ")
        end
    end

    function jlstr(res)
        jlcmd = benchmarkconfig(res).juliacmd
        flags = length(jlcmd) <= 1 ? [] : jlcmd[2:end]
        return if isempty(flags)
            "None"
        else
            """`$(join(flags, ","))`"""
        end
    end

    println(io, """
                # Benchmark Report for *$(name(target))*

                ## Job Properties
                * Time of benchmarks:
                    - Target: $(Dates.format(date(target), "d u Y - HH:MM"))
                    - Baseline: $(Dates.format(date(baseline), "d u Y - HH:MM"))
                * Package commits:
                    - Target: $(commit(target)[1:min(6, length(commit(target)))])
                    - Baseline: $(commit(baseline)[1:min(6, length(commit(baseline)))])
                * Julia commits:
                    - Target: $(juliacommit(target)[1:min(6, length(juliacommit(target)))])
                    - Baseline: $(juliacommit(baseline)[1:min(6, length(juliacommit(baseline)))])
                * Julia command flags:
                    - Target: $(jlstr(target))
                    - Baseline: $(jlstr(baseline))
                * Environment variables:
                    - Target: $(env_strs(target))
                    - Baseline: $(env_strs(baseline))
                """)

    entries = BenchmarkTools.leaves(benchmarkgroup(judgement))
    entries = entries[sortperm(map(x -> string(first(x)), entries))]

    cw = [2, 10, 12]
    for (ids, t) in entries
        _update_col_widths!(cw, ids, t)
    end

    if export_invariants
        print(io, """
                    ## Results
                    A ratio greater than `1.0` denotes a possible regression (marked with $(_REGRESS_MARK)), while a ratio less
                    than `1.0` denotes a possible improvement (marked with $(_IMPROVE_MARK)). All results are shown below.

                    | ID$(" "^(cw[1]-2)) | time ratio$(" "^(cw[2]-10)) | memory ratio$(" "^(cw[3]-12)) |
                    |---$("-"^(cw[1]-2))-|-----------$("-"^(cw[2]-10))-|-------------$("-"^(cw[3]-12))-|
                    """)

        for (ids, t) in entries
            println(io, _resultrow(ids, t, cw))
        end
    else
        print(io, """
                    ## Results
                    A ratio greater than `1.0` denotes a possible regression (marked with $(_REGRESS_MARK)), while a ratio less
                    than `1.0` denotes a possible improvement (marked with $(_IMPROVE_MARK)). Only significant results - results
                    that indicate possible regressions or improvements - are shown below (thus, an empty table means that all
                    benchmark results remained invariant between builds).

                    | ID$(" "^(cw[1]-2)) | time ratio$(" "^(cw[2]-10)) | memory ratio$(" "^(cw[3]-12)) |
                    |---$("-"^(cw[1]-2))-|-----------$("-"^(cw[2]-10))-|-------------$("-"^(cw[3]-12))-|
                    """)

        for (ids, t) in entries
            if BenchmarkTools.isregression(t) || BenchmarkTools.isimprovement(t)
                println(io, _resultrow(ids, t, cw))
            end
        end
    end

    println(io)
    println(io, """
                ## Benchmark Group List
                Here's a list of all the benchmark groups executed by this job:
                """)

    for id in unique(map(pair -> pair[1][1:end-1], entries))
        println(io, "- `", _idrepr(id), "`")
    end

    println(io)
    println(io, "## Julia versioninfo")

    println(io, "\n### Target")
    print(io, "```\n", versioninfo(target), "```")

    println(io, "\n\n### Baseline")
    print(io, "```\n", versioninfo(baseline), "```")

    return nothing
end


text:
module TestOutputSize

include("preamble.jl")
using Base: HasLength, HasShape, IteratorSize, SizeUnknown
using Transducers: OutputSize, SizeStable, SizeChanging

@testset "SizeStable" begin
    @testset for xf in [
        Map(sin),
        MapSplat(*),
        Replace((a = 1,)),
        FlagFirst(),
        Iterated(identity, nothing),
        Count(),
        GetIndex(1:3),
        SetIndex([]),
        Enumerate(),
        Broadcasting(),
        opcompose(MapSplat(*), Map(sin)),
        opcompose(MapSplat(*), Map(sin), FlagFirst()),
    ]
        @test OutputSize(typeof(xf)) isa SizeStable

        @testset for xs in [1:3, reshape(1:6, 2, 3), 1]
            @test IteratorSize(xs |> xf) isa HasShape
            @test size(xs |> xf) == size(xs)
            @test length(xs |> xf) == length(xs)
        end

        @testset for xs in [(1, 2), (a = 1, b = 2)]
            @test IteratorSize(xs |> xf) isa HasLength
            @test length(xs |> xf) == length(xs)
            @test_error size(xs |> xf)
        end

        @testset """eachline(IOBuffer("a"))""" begin
            xs = eachline(IOBuffer("a"))
            @test IteratorSize(xs |> xf) isa SizeUnknown
            @test_error size(xs |> xf)
            @test_error length(xs |> xf)
        end
    end
end

@testset "SizeChanging" begin
    @testset for xf in [
        Cat(),
        opcompose(Cat(), MapSplat(*)),
        opcompose(Cat(), MapSplat(*), Map(sin)),
    ]
        @test OutputSize(typeof(xf)) isa SizeChanging

        @testset for xs in [
            1:3,
            reshape(1:6, 2, 3),
            1,
            (1, 2),
            (a = 1, b = 2),
            eachline(IOBuffer("a")),
        ]
            @test IteratorSize(xs |> xf) isa SizeUnknown
            @test_error size(xs |> xf)
            @test_error length(xs |> xf)
        end
    end
end

end  # module


text:
# Julia wrapper for header: /usr/include/libavdevice/avdevice.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avdevice_version,
    avdevice_configuration,
    avdevice_license,
    avdevice_register_all,
    av_input_audio_device_next,
    av_input_video_device_next,
    av_output_audio_device_next,
    av_output_video_device_next,
    avdevice_app_to_dev_control_message,
    avdevice_dev_to_app_control_message,
    avdevice_capabilities_create,
    avdevice_capabilities_free,
    avdevice_list_devices,
    avdevice_free_list_devices,
    avdevice_list_input_sources,
    avdevice_list_output_sinks

function avdevice_version()
    ccall((:avdevice_version, libavdevice), UInt32, ())
end

function avdevice_configuration()
    ccall((:avdevice_configuration, libavdevice), Cstring, ())
end

function avdevice_license()
    ccall((:avdevice_license, libavdevice), Cstring, ())
end

function avdevice_register_all()
    ccall((:avdevice_register_all, libavdevice), Cvoid, ())
end

function av_input_audio_device_next(d)
    ccall((:av_input_audio_device_next, libavdevice), Ptr{AVInputFormat}, (Ptr{AVInputFormat},), d)
end

function av_input_video_device_next(d)
    ccall((:av_input_video_device_next, libavdevice), Ptr{AVInputFormat}, (Ptr{AVInputFormat},), d)
end

function av_output_audio_device_next(d)
    ccall((:av_output_audio_device_next, libavdevice), Ptr{AVOutputFormat}, (Ptr{AVOutputFormat},), d)
end

function av_output_video_device_next(d)
    ccall((:av_output_video_device_next, libavdevice), Ptr{AVOutputFormat}, (Ptr{AVOutputFormat},), d)
end

function avdevice_app_to_dev_control_message(s, _type::AVAppToDevMessageType, data, data_size::Csize_t)
    ccall((:avdevice_app_to_dev_control_message, libavdevice), Cint, (Ptr{AVFormatContext}, AVAppToDevMessageType, Ptr{Cvoid}, Csize_t), s, _type, data, data_size)
end

function avdevice_dev_to_app_control_message(s, _type::AVDevToAppMessageType, data, data_size::Csize_t)
    ccall((:avdevice_dev_to_app_control_message, libavdevice), Cint, (Ptr{AVFormatContext}, AVDevToAppMessageType, Ptr{Cvoid}, Csize_t), s, _type, data, data_size)
end

function avdevice_capabilities_create(caps, s, device_options)
    ccall((:avdevice_capabilities_create, libavdevice), Cint, (Ptr{Ptr{AVDeviceCapabilitiesQuery}}, Ptr{AVFormatContext}, Ref{Ptr{AVDictionary}}), caps, s, device_options)
end

function avdevice_capabilities_free(caps, s)
    ccall((:avdevice_capabilities_free, libavdevice), Cvoid, (Ptr{Ptr{AVDeviceCapabilitiesQuery}}, Ptr{AVFormatContext}), caps, s)
end

function avdevice_list_devices(s, device_list)
    ccall((:avdevice_list_devices, libavdevice), Cint, (Ptr{AVFormatContext}, Ptr{Ptr{AVDeviceInfoList}}), s, device_list)
end

function avdevice_free_list_devices(device_list)
    ccall((:avdevice_free_list_devices, libavdevice), Cvoid, (Ptr{Ptr{AVDeviceInfoList}},), device_list)
end

function avdevice_list_input_sources(device, device_name, device_options, device_list)
    ccall((:avdevice_list_input_sources, libavdevice), Cint, (Ptr{AVInputFormat}, Cstring, Ptr{AVDictionary}, Ptr{Ptr{AVDeviceInfoList}}), device, device_name, device_options, device_list)
end

function avdevice_list_output_sinks(device, device_name, device_options, device_list)
    ccall((:avdevice_list_output_sinks, libavdevice), Cint, (Ptr{AVOutputFormat}, Cstring, Ptr{AVDictionary}, Ptr{Ptr{AVDeviceInfoList}}), device, device_name, device_options, device_list)
end


text:
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end
parsed_args = parse_args(ARGS, s)

function parseLine(l)
    [parse(UInt8, d) for d in l]
end

test_suffix = parsed_args["test"] ? "_test" : ""
depths = [parseLine(l) for l in readlines("input/day_9$test_suffix.txt")]
mx = length(depths)
my = length(depths[1])
function getValue(x, y)
    if x <= 0 || x > mx || y <= 0 || y > my
        return typemax(UInt8)
    end
    depths[x][y]
end

sinks = []
for x = 1:mx
    for y = 1:my
        v = getValue(x, y)
        if getValue(x - 1, y) > v && getValue(x + 1, y) > v && getValue(x, y - 1) > v && getValue(x, y + 1) > v
            push!(sinks, (x = x, y = y, v = v))
        end
    end
end

function part1()
    sum(s -> s.v + 1, sinks)
end

function basinsize(sink)
    basin = Set()
    function addtobasin(x, y)
        if (x, y) ∉ basin && getValue(x, y) < 9
            push!(basin, (x, y))
            addtobasin(x - 1, y)
            addtobasin(x + 1, y)
            addtobasin(x, y - 1)
            addtobasin(x, y + 1)
        end
    end
    addtobasin(sink.x, sink.y)
    length(basin)
end

function part2()
    sizes = sort(map(basinsize, sinks))
    s1 = pop!(sizes)
    s2 = pop!(sizes)
    s3 = pop!(sizes)
    s1 * s2 * s3
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
module IrisClustering

using Logging, LoggingExtras

function main()
  Base.eval(Main, :(const UserApp = IrisClustering))

  include(joinpath("..", "genie.jl"))
end; main()

end


text:
e1 = parse("x * sqrt(9) - 5")
e2 = :(x * sqrt(9) - 5)
e3 = Expr(:call, :-, Expr(:call, :*, :x, Expr(:call, :sqrt, 9)), 5)

@assert e1 == e2 && e2 == e3 && typeof(e1) == Expr


text:
em = ErgodicManagerR2("double gaussian", K=5, bins=100)

# Trajectory manager
x0 = [0.89,0.01]
N = 400
h = 0.1
ci = ConstantInitializer([-0.01,0.01])
tm = TrajectoryManager(x0, h, N, ci)

xd, ud = smc_trajectory(em, tm, umax=.15)

# now test with two agents...

# Set up different domains with different discretizations
d = Domain([1,1], 100)
num_agents = 2

# Set up distribution and ergodic manager
K = 5
means = [[.3,.7], [.7,.3]]
Sigmas = [.025*eye(2), .025*eye(2)]
phi = gaussian(d, means, Sigmas)
em = ErgodicManagerR2(d, phi, K)

# Set up first trajectory manager
x0 = [0.49,0.01]
N = 50
h = 0.6
ci = ConstantInitializer([0.0, 0.0])
tm1 = TrajectoryManager(x0, h, N, ci)
dynamics!(tm1, SingleIntegrator(2,h))

# second tm is like the first, but different starting point
tm2 = deepcopy(tm1)
tm2.x0 = [.79,.99]

# array of trajectory managers
vtm = [tm1, tm2]

# Generate the trajectories
ddc = 1e-4
xd, ud = smc_trajectory(em, vtm, umax=.15)


text:
"""
$(TYPEDEF)

# Fields

Accessing fields directly is part of the API. The additional property `cap` is supported.

$(FIELDS)
"""
struct HomogeneousBasis{N,U}
    "The univariate family, eg `Chebyshev()`."
    univariate::U
    "Indexes by dimensions."
    I::NTuple{N,Int}
    "The level (`≥ 0`) of the Smolyak approximation. `cap = level + N`."
    level::Int
end

function Base.propertynames(basis::HomogeneousBasis, private = false)
    (fieldnames(HomogeneousBasis)..., :cap)
end

function Base.getproperty(basis::HomogeneousBasis{N}, name::Symbol) where {N}
    if name ≡ :cap
        N + basis.level
    else
        getfield(basis, name)
    end
end

function degrees_of_freedom(basis::HomogeneousBasis)
    @unpack univariate, I, cap = basis
    C = mapfoldl((C1, C2) -> count_combinations(cap, C1, C2), I) do i
        indexes = 1:min(i, cap)
        Counts(indexes, set_length.(univariate, indexes))
    end
    sum(C)
end

"""
$(SIGNATURES)

Calculate the matrix `A[R[i[end]], R[j[end]] ⊗ … ⊗ A[R[i[2]], R[j[2]] ⊗ A[R[i[1]], R[j[1]]`,
where

- `A` is the basis matrix for all points and basis functions,

- `R` contains the ranges for each index,

- `i` and `j` are iterables of indexes in `R`.

Internal helper function used in building basis matrices.
"""
function _basis_block(A, R, i, j)
    mapreduce(((i, j), ) -> A[R[i], R[j]], (A1, A2) -> kron(A2, A1), zip(i, j))
end

"""
$(SIGNATURES)

Calculate all combinations of `x[R[i[1]]], x[R[i[2]]], …`, with the first index changing the
fastest (column major).

Internal helper function used in building coordinates.
"""
function _coordinate_block(x, R, i)
    mapfoldl(i -> x[R[i]],
             (y, x) -> vec(((x, y) -> (x, y...)).(x, permutedims(y))),
             reverse(i), init = [()])
end

"""
$(SIGNATURES)

Return `A, v` where `A` is a ``d × d`` matrix and `v` is a `d`-element vector, where

- `d = degrees_of_freedom(basis)`,

- `v` is a vector of *coordinates*,

- `A \\ f.(x)` would calculate the *coefficients* of approximating `f`.
"""
function basis_matrix_and_coordinates(basis::HomogeneousBasis)
    @unpack univariate, I, cap = basis
    iter = CappedCartesianIndices(cap, I)
    K = maximum(largest_indices(iter))

    # univariate building blocks
    R = set_range.(univariate, 1:K)
    x0 = all_gridpoints(univariate, K)
    A0 = evaluate(last(last(R)), univariate, x0)

    # combinations
    # FIXME these could be made much more efficient by pre-allocating, and making
    # `_basis_block` and `_coordinate_block` mutate a pre-allocated matrix
    A = mapreduce(j -> mapreduce(i -> _basis_block(A0, R, i, j), vcat, iter), hcat, iter)
    x = mapreduce(i -> _coordinate_block(x0, R, i), vcat, iter)
    A, SVector.(x)
end

"""
$(SIGNATURES)

Helper function for iterating when interpolating `basis` at `x`.

Return a `NamedTuple` of

- `iter`: an iterator for the set indices,
- `R`: a vector of set ranges
- `A0s`: basis matrices for each coordinate, to be indexed by elements of `R`
"""
function _interpolation_helper(basis::HomogeneousBasis{N}, x::AbstractVector) where {N}
    @argcheck length(x) == N
    @unpack univariate, I, cap = basis
    iter = CappedCartesianIndices(cap, I)
    Ks = largest_indices(iter)
    R = set_range.(univariate, 1:maximum(Ks))
    A0s = map((K, x) -> evaluate(last(R[K]), univariate, x), Ks, x)
    (iter = iter, R = R, A0s = A0s)
end

"""
$(SIGNATURES)

Interpolate `basis` with coefficients `c` at `x`.
"""
function interpolate(basis::HomogeneousBasis, c::AbstractVector{T1},
                     x::AbstractVector{T2}) where {T1, T2}
    @unpack iter, R, A0s = _interpolation_helper(basis, x)
    S = float(promote_type(T1, T2))
    v = zero(S)
    p = 1
    for i in iter
        for I in CartesianIndices(map(i -> R[i], i))
            v += prod(map(getindex, A0s, Tuple(I))) * c[p]
            p += 1
        end
    end
    v
end

"""
$(SIGNATURES)

In-place version of [`interpolated_basis`](@ref). Overwrites the first argument,
which should have length `degrees_of_freedom(basis)`.
"""
function interpolated_basis!(a::AbstractVector, basis::HomogeneousBasis, x::AbstractVector)
    @unpack iter, R, A0s = _interpolation_helper(basis, x)
    p = 1
    for i in iter
        for I in CartesianIndices(map(i -> R[i], i))
            a[p] = prod(map(getindex, A0s, Tuple(I)))
            p += 1
        end
    end
    a
end

"""
$(SIGNATURES)

Return a vector `a` such that

```julia
dot(a, c) ≈ interpolate(basis, c, x)
```

This is useful when repeatedly interpolating at the same set of coordinates `x` (`a` is
specific to `x`).
"""
function interpolated_basis(basis::HomogeneousBasis, x::AbstractVector{T}) where{T}
    interpolated_basis!(Vector{float(T)}(undef, degrees_of_freedom(basis)), basis, x)
end


text:
using Flux

using Flux: params, sigmoid
using Flux.Losses: binarycrossentropy
using Zygote: @nograd
using AccuracyAtTop
using AdversarialPrediction

# -------------------------------------------------------------------------------
# Models
# -------------------------------------------------------------------------------
abstract type Model end

sqsum(x) = sum(abs2, x)

# Baseline model
abstract type BaseLine <: Model end
abstract type BaseLineSimple <: Model end

weights(y, s) = oftype(s, y ./ sum(y .== 1) .+ (1 .- y) ./ sum(y .== 0))

@nograd weights

function build_loss(::Type{<:BaseLine}, arg, surrogate, reg)

    function loss(x, y, model, pars)
        return loss(y, model(x), pars)
    end

    function loss(y, s, pars)
        w = weights(y, s)
        return binarycrossentropy(sigmoid.(s), y; agg = x -> sum(w .* x)) + eltype(s)(reg) * sum(sqsum, pars)
    end
    return loss
end

function build_loss(::Type{<:BaseLineSimple}, arg, surrogate, reg)

    function loss(x, y, model, pars)
        return loss(y, model(x), pars)
    end

    function loss(y, s, pars)
        return binarycrossentropy(sigmoid.(s), y) + eltype(s)(reg) * sum(sqsum, pars)
    end
    return loss
end

# Threshold models
abstract type DeepTopPush <: Model end
abstract type PatMat <: Model end
abstract type PatMatNP <: Model end

build_threshold(::Type{DeepTopPush}, arg) = Maximum(; samples = NegSamples)
build_threshold(::Type{PatMat}, τ::Real) = Quantile(τ; samples = AllSamples, rev = true)
build_threshold(::Type{PatMatNP}, τ::Real) = FPRate(τ)

function build_loss(type, arg, surrogate, reg)
    thres = build_threshold(type, arg)

    function loss(x, y, model, pars)
        return loss(y, model(x), pars)
    end

    function loss(y, s, pars)
        t = threshold(thres, y, s)
        return fnr(y, s, t, surrogate) + eltype(s)(reg) * sum(sqsum, pars)
    end
    return loss
end

# Ap-Perf model
abstract type APPerf <: Model end

@metric FnrAtFpr th
function define(::Type{FnrAtFpr}, C::AdversarialPrediction.ConfusionMatrix, th)
    return C.fn / C.pp
end

function constraint(::Type{FnrAtFpr}, C::AdversarialPrediction.ConfusionMatrix, th)
    return C.fp / C.nn <= th
end

function build_loss(::Type{<:APPerf}, arg, surrogate, reg; admm_iter = 100)

    metric = FnrAtFpr(arg)
    special_case_positive!(metric)
    cs_special_case_positive!(metric, false)

    function loss(x, y, model, pars)
        return loss(y, model(x), pars)
    end

    function loss(y, s, pars)
        return ap_objective(vec(s), vec(y), metric; max_iter = admm_iter) + eltype(s)(reg) * sum(sqsum, pars)
    end
    return loss
end

# TFCO model
abstract type TFCO <: Model end


text:
module Stan

using StanBase
using StanSample
using StanOptimize
using StanDiagnose
using StanVariational
using StanQuap

end # module


text:
#!/usr/bin/env julia

# Start Test Script
using NewtonMethods
using Base.Test

include("ChebyQuad.jl")

# Run tests
tic()
println("Testing NLEQ1")
@time @test include("testNLEQ1.jl")

println("Testing NLEQ2")
@time @test include("testNLEQ2.jl")
toc()


text:
"""
	Estimates the integral of f between a and b by the trapezoidal rule

    trapz(f, a, b, n)

	 Parameters
	 ----------
     `f` is the function to integrate
     `a` lower bound for the interval
     `b` upper bound for the interval
     `n` number of partition of [a,b]
"""

function trapz(f, a, b, n)
    h = (b - a)/n
    integral = f(a) + f(b)
    i = 1
    while i < n
        x1 = a + i*h
        integral += 2*f(x1)
        i += 1
    end
    integral *= h/2
    return integral
end


text:
## LocalRule
export NagyRule

struct NagyRule{T<:Vector} <: MCMCRule
    adjacency_list::T
end

function NagyRule(H::AbsLinearOperator)
    N_sites = nsites(basis(H))

    all_sites = sites(H)
    if !(eltype(all_sites) <: Vector)
        all_sites = [all_sites]
    end

    couplings = Vector{Tuple{Int,Int}}()
    for coupling=all_sites
        length(coupling) == 1 && continue
        if length(coupling) > 2
            @warn "Can't exchange between 3-site couplings. This coupling is ignored"
            continue
        end
        push!(couplings, (coupling[1], coupling[2]))
    end

    return NagyRule{typeof(couplings)}(couplings)
end

struct NagyRuleCache
end

function RuleSamplerCache(rule::NagyRule, s, v, net, par_cache)
    return NagyRuleCache()
end

function propose_step!(σp::ADoubleState, s::MetropolisSampler{<:NagyRule},
                       net::Union{MatrixNet,KetNet}, c, rc)
    adj_list = s.rule.adjacency_list
    hilb     = physical(c.hilb)

    # determine move type
    move_type = rand(c.rng, 1:8)#1:n_moves)
    sites_to_flip = 1:nsites(hilb)

    # 1- Flip a spin
    if move_type == 1 || move_type == 2                 # hopping in σ
        flipat_r1 = rand(c.rng, sites_to_flip)
        flipat_r2 = rand(c.rng, adj_list[flipat_r1])

        @inbounds flipat!(c.rng, row(σp), hilb, flipat_r1)
        @inbounds flipat!(c.rng, row(σp), hilb, flipat_r2)

    elseif move_type == 3 || move_type == 4             # hopping in η
        flipat_c1 = rand(c.rng, sites_to_flip)
        flipat_c2 = rand(c.rng, adj_list[flipat_c1])

        @inbounds flipat!(c.rng, col(σp), hilb, flipat_c1)
        @inbounds flipat!(c.rng, col(σp), hilb, flipat_c2)

    elseif move_type == 5                               # exctitation in σ
        flipat_r1 = rand(c.rng, sites_to_flip)

        @inbounds flipat!(c.rng, row(σp), hilb, flipat_r1)

    elseif move_type == 6                               # excitations in η
        flipat_r1 = rand(c.rng, sites_to_flip)

        @inbounds flipat!(c.rng, col(σp), hilb, flipat_r1)

    elseif move_type == 7                                # Dissipator
        flipat_r  = rand(c.rng, sites_to_flip)
        old_val_r = row(σp)[flipat_r]

        if old_val_r == 0
            if rand(c.rng, 1:10) == 1 # excite with 10% chance
                @inbounds flipat!(c.rng, row(σp), hilb, flipat_r)
            end
        else # always dissipate
            @inbounds flipat!(c.rng, row(σp), hilb, flipat_r)
        end

        flipat_c  = flipat_r #rand(c.rng, sites_to_flip)
        old_val_c = col(σp)[flipat_c]
        if old_val_c == 0
            if rand(c.rng, 1:10) == 1 # excite with 10% chance
                @inbounds flipat!(c.rng, col(σp), hilb, flipat_c)
            end
        else # always dissipate
            @inbounds flipat!(c.rng, col(σp), hilb, flipat_c)
        end

    elseif move_type == 8                                # jumper
        flipat_r = rand(c.rng, sites_to_flip)
        flipat_c = rand(c.rng, sites_to_flip)

        @inbounds old_val_r, _ =flipat!(c.rng, row(σp), hilb, flipat_r)
        @inbounds old_val_c, _ =flipat!(c.rng, col(σp), hilb, flipat_c)
    end
#    elseif move_type == 9
#        if rand(c.rng, 1:100) == 1
#            zero!(σ)
#        end
        #prob = exp(2*real(net(σ)) - logψ_σ)
        #rv = rand(c.rng)
        #if !(rv < prob)
        #    setat!(row(σ), flipat_r, old_val_r)
        #    setat!(col(σ), flipat_c, old_val_c)
        #else
        #end
#    end

    return true
end


text:
# Test that a hitran fixed-width file is correctly parsed
@testset "read_hitran" begin

    println("Testing read_hitran...")

    ####
    #### These tests check that the correct values & types are extracted from the file,
    #### especially that molecule/isotope match and ν between ν_min and ν_max)
    ####

    CO2_test_file = "test_profiles/testCO2.data"

    test_ht = Absorption.read_hitran(CO2_test_file, mol=2, iso=1, ν_min=6000, ν_max=6400)
    @test test_ht.mol == [2, 2, 2, 2] && test_ht.mol[1] isa Int64
    @test test_ht.iso == [1, 1, 1, 1] && test_ht.iso[1] isa Int64
    @test test_ht.νᵢ == [6000.542970, 6286.403343, 6317.417493, 6380.824116] && test_ht.νᵢ[1] isa Float64
    @test test_ht.Sᵢ == [ 1.098E-28, 9.843E-30, 5.613E-27, 1.809E-30] && test_ht.Sᵢ[1] isa Float64
    @test test_ht.Aᵢ == [ 9.993e-08, 1.179e-08, 1.324e-05, 1.601e-02] && test_ht.Aᵢ[1] isa Float64
    @test test_ht.γ_air == [.0880, .0687, .0682, .0671] && test_ht.γ_air[1] isa Float64
    @test test_ht.γ_self == [0.118, 0.087, 0.081, 0.073] && test_ht.γ_self[1] isa Float64
    @test test_ht.E″ == [7.8043, 464.1717, 639.6004, 3798.2095] && test_ht.E″[1] isa Float64
    @test test_ht.n_air == [0.77, 0.76, 0.76, 0.73] && test_ht.n_air[1] isa Float64
    @test test_ht.δ_air == [-.004342, -.007362, -.007443, -.007669] && test_ht.δ_air[1] isa Float64
    @test test_ht.global_upper_quanta == ["       4 1 1 03", "       2 2 2 12", "       2 2 2 12", "       4 2 2 12"] && test_ht.global_upper_quanta[1] isa String
    @test test_ht.global_lower_quanta == ["       0 0 0 01", "       0 0 0 01", "       0 0 0 01", "       1 2 2 01"] && test_ht.global_lower_quanta[1] isa String
    @test test_ht.local_upper_quanta == ["               ", "               ", "               ", "               "] && test_ht.local_upper_quanta[1] isa String
    @test test_ht.local_lower_quanta == ["     Q  4e     ", "     Q 34e     ", "     R 40e     ", "     R 51f     "] && test_ht.local_lower_quanta[1] isa String
    @test test_ht.ierr == ["367774", "367764", "367764", "367774"] && test_ht.ierr[1] isa String
    @test test_ht.iref == ["2029 5 4 5 7", "2029 5 4 5 7", "2029 5 4 5 7", "2029 5 4 5 7"] && test_ht.iref[1] isa String
    @test test_ht.line_mixing_flag == [" ", " ", " ", " "] && test_ht.line_mixing_flag[1] isa String
    @test test_ht.g′ == [9.0, 69.0, 83.0, 105.0] && test_ht.g′[1] isa Float64
    @test test_ht.g″ == [9.0, 69.0, 81.0, 103.0] && test_ht.g″[1] isa Float64

    ####
    #### These tests check that the optional params work correctly
    ####

    # Not specifying the molecule #
    test_ht = Absorption.read_hitran(CO2_test_file, iso=1, ν_min=6000, ν_max=6400)
    @test test_ht.mol == [1, 2, 2, 2, 2]
    @test test_ht.iso == [1, 1, 1, 1, 1]
    @test test_ht.νᵢ == [6286.403343, 6000.542970, 6286.403343, 6317.417493, 6380.824116]
    @test test_ht.g″ == [69.0, 9.0, 69.0, 81.0, 103.0]

    # Not specifying the isotope #
    test_ht = Absorption.read_hitran(CO2_test_file, mol=2, ν_min=6000, ν_max=6400)
    @test test_ht.mol == [2, 2, 2, 2, 2]
    @test test_ht.iso == [2, 1, 1, 1, 1]
    @test test_ht.νᵢ == [6000.542970, 6000.542970, 6286.403343, 6317.417493, 6380.824116]
    @test test_ht.g″ == [9.0, 9.0, 69.0, 81.0, 103.0]

    # Not specifying the molecule # OR isotope #
    test_ht = Absorption.read_hitran(CO2_test_file, ν_min=6000, ν_max=6400)
    @test test_ht.mol == [1, 2, 2, 2, 2, 2]
    @test test_ht.iso == [1, 2, 1, 1, 1, 1]
    @test test_ht.νᵢ == [6286.403343, 6000.542970, 6000.542970, 6286.403343, 6317.417493, 6380.824116]
    @test test_ht.g″ == [69.0, 9.0, 9.0, 69.0, 81.0, 103.0]

    # Not specifying ν_min
    test_ht = Absorption.read_hitran(CO2_test_file, mol=2, iso=1, ν_max=6400)
    @test length(test_ht.mol) == 9

    # Not specifying ν_max
    test_ht = Absorption.read_hitran(CO2_test_file, mol=2, iso=1, ν_min=6000)
    @test length(test_ht.mol) == 7

    # Not specifying ν_min OR ν_max
    test_ht = Absorption.read_hitran(CO2_test_file, mol=2, iso=1)
    @test length(test_ht.mol) == 12

end

# Test that absorption cross sections are calculated correctly
# (Not using pre-saved interpolator)

@testset "absorption_cross_section_hitran" begin

    println("Testing absorption_cross_section_hitran...")

    ####
    #### These tests check that for one molecule (CO2), over a temperature/
    #### pressure grid and in a specific band (6000<ν<6400), our Voigt
    #### broadening implementation closely matches the HAPI output under
    #### the same conditions
    ####

    # Temperature and pressure grids.
    # Note that these are pre-defined for testing -- the values are defined in generateHapiTests.py
    # To test against a different range, you must change the grids in that
    # file and rerun it. Then, you can change it here.

    pressures = 250:250:1250
    temperatures = 100:75:400

    # Get the test data
    CO2_file = artifact("CO2")
    test_ht = Absorption.read_hitran(CO2_file, mol=2, iso=1, ν_min=6000, ν_max=6400)

    grid = 6000:0.01:6400;

    # Threshold -- our value must be within ϵ of the HAPI value
    ϵ = 3.6e-27

    # Create a HitranModel
    model = make_hitran_model(test_ht, Voigt(), CEF=HumlicekWeidemann32SDErrorFunction())

    # Loop over every temperature/pressure combo and test that the results match HAPI
    @showprogress 1 "Testing HAPI equivalence (On CO2 Band)..." for temp in temperatures
        for pres in pressures
            jl_cs = absorption_cross_section(model, grid, pres, temp)
            py_cs = array_type(default_architecture)(readdlm("test_profiles/Voigt_CO2_T" * string(temp) * "_P" * string(pres) * ".csv"))
            Δcs = abs.(jl_cs - py_cs)
            @test maximum(Δcs) < ϵ
        end
    end

    ####
    #### Now test HAPI equivalence with other molecules
    ####

    names = ["H2O", "CO2", "O3", "N2O", "CO"]
    molecules = [1, 2, 3, 4, 5]
    isotopes = [1, 1, 1, 1, 1]

    # Loop over every temperature/pressure combo and test that the results match HAPI
    # (Doing this for other molecules)
    @showprogress 1 "Testing HAPI equivalence (On Other Molecules)..." for name in names

        pres = 1000
        temp = 250

        # Get the test data
        test_ht = Absorption.read_hitran(artifact(name), iso=1, ν_min=6000, ν_max=6400)
        # Create a HitranModel
        model = make_hitran_model(test_ht, Voigt(), CEF=HumlicekWeidemann32SDErrorFunction())

        jl_cs = absorption_cross_section(model, grid, pres, temp)
        py_cs = array_type(default_architecture)(readdlm("test_profiles/Voigt_" * name * "_T250_P1000.csv"))
        Δcs = abs.(jl_cs - py_cs)
        @test maximum(Δcs) < ϵ
    end
end

# Test that absorption cross sections if calculated using wavenumbers

@testset "absorption_cross_section_wavelengths" begin

    println("Testing absorption_cross_section_wavelengths...")

    hitran_data = read_hitran(artifact("CO2"), mol=2, iso=1, ν_min=6000, ν_max=6400)
    model = make_hitran_model(hitran_data, Voigt(), wing_cutoff = 40, CEF=HumlicekWeidemann32SDErrorFunction(), architecture=CPU())
    grid = 6000:0.01:6400
    @test absorption_cross_section(model, grid, 1000.1, 296.1) ≈ reverse(absorption_cross_section(model, reverse(1e7 ./ collect(grid)), 1000.1, 296.1, wavelength_flag=true));
end

# Test that absorption cross sections are calculated correctly
# using a new interpolator

@testset "absorption_cross_section_interpolator" begin

    println("Testing absorption_cross_section_interpolator...")

    # Get the test data
    test_ht = Absorption.read_hitran(artifact("CO2"), iso=1, ν_min=6000, ν_max=6400)

    # Pressure and temperature grids
    pressures = 250:250:1250
    temperatures = 100:75:400

    # Wavelength grid
    ν_grid = 6000:0.01:6400

    # Create the Interpolation Model
    interp_model = make_interpolation_model(test_ht, Voigt(), ν_grid, pressures,
                                            temperatures, CEF=HumlicekWeidemann32SDErrorFunction())

    # Threshold -- our value must be within ϵ of the HAPI value
    ϵ = 3.6e-27

    # Loop over every temperature/pressure combo and test that the results match HAPI
    for temp in temperatures
        for pres in pressures
            jl_cs = absorption_cross_section(interp_model, ν_grid, pres, temp)
            py_cs = readdlm("test_profiles/Voigt_CO2_T" * string(temp) * "_P" * string(pres) * ".csv")
            Δcs = abs.(jl_cs - py_cs)
            @test maximum(Δcs) < ϵ
        end
    end

end

# Test that checks whether the cross-section auto-differentiation works

@testset "absorption_cross_section_autodiff" begin

    println("Testing absorption_cross_section_autodiff...")

    # Load HITRAN data
    hitran_data = read_hitran(artifact("CO2"), mol=2, iso=1, ν_min=6000, ν_max=6400)

    # Create the model with parameters
    model = make_hitran_model(hitran_data, Voigt(), wing_cutoff = 40, CEF=HumlicekWeidemann32SDErrorFunction(), architecture=CPU())

    # Compute the cross-section with autodifferentiation
    value, derivs = absorption_cross_section(model, 6000:0.01:6400, 1000.1, 296.1, autodiff=true);

end


text:
using CmdStan, Turing, ContinuousBenchmarks

# Model taken from https://github.com/goedman/Stan.jl/blob/master/Examples/Mamba/EightSchools/schools8.jl

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "school8-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "school8-stan.model.jl"))

global stanmodel, rc, sim
# stanmodel = Stanmodel(name="schools8", model=eightschools);
stanmodel = Stanmodel(Sample(algorithm=CmdStan.Hmc(CmdStan.Static(0.75*5),CmdStan.diag_e(),0.75,0.0),
  save_warmup=true,adapt=CmdStan.Adapt(engaged=false)),
  num_samples=2000, num_warmup=0, thin=1,
  name="schools8", model=eightschools, nchains=1);

rc, sim = stan(stanmodel, schools8data, CmdStanDir=ContinuousBenchmarks.CMDSTAN_HOME, summary=false)

stan_d = Dict()

for i = 1:8
  stan_d["eta[$i]"] = sim[:, ["eta.$i"], :].value[:]
  stan_d["theta[$i]"] = sim[:, ["theta.$i"], :].value[:]
end

stan_d["mu"] = sim[:, ["mu"], :].value[:]
stan_d["tau"] =sim[:, ["tau"], :].value[:]

for k = keys(stan_d)
  stan_d[k] = mean(stan_d[k])
end

# println("CmdStan time: $stan_time")

# describe(sim)


text:
module DensifyInnerMostBlockLevelMapsTests
   using Test
   using Gridap
   using Gridap.Fields
   using Gridap.Arrays

   m=DensifyInnerMostBlockLevelMap()

   # VectorBlock{Vector} test
   y=Vector{Vector{Float64}}(undef,3)
   touched=Vector{Bool}(undef,3)
   touched .= true
   x=rand(3)
   y[1]=x
   y[2]=x.+1
   y[3]=x.+2
   by=ArrayBlock(y,touched)
   cache=return_cache(m,by)
   @test vcat(y...)==evaluate!(cache,m,by)

   # MatrixBlock{Vector} test
   x=rand(3)
   y=Array{Vector{Float64},2}(undef,(1,4))
   y[1,1]=x
   y[1,3]=x.+2
   y[1,4]=x.+3
   touched  = Array{Bool,2}(undef,(1,4))
   touched .= true
   touched[1,2]=false
   by=ArrayBlock(y,touched)
   cache=return_cache(m,by)
   @test hcat(y[1,1],zeros(3),y[1,3],y[1,4]) == evaluate!(cache,m,by)

   # VectorBlock{Matrix} test
   x=rand(1,3)
   y=Array{Matrix{Float64}}(undef,(4,))
   y[1]=x
   y[2]=x.+1
   y[3]=x.+2
   y[4]=x.+3
   touched  = Array{Bool,1}(undef,(4,))
   touched .= true
   by=ArrayBlock(y,touched)
   cache=return_cache(m,by)
   @test vcat(y...)==evaluate!(cache,m,by)

   # MatrixBlock{Matrix} test
   x1=rand(2,3)
   x2=rand(3,3)
   x3=rand(4,3)
   y=Array{Matrix{Float64}}(undef,(3,3))
   y[1,1]=x
   y[2,1]=x2
   y[3,1]=x3
   y[1,2]=x.+3
   y[1,3]=x.+5
   touched  = Array{Bool,2}(undef,(3,3))
   touched .= true
   touched[2:3,2:3].=false
   by=ArrayBlock(y,touched)
   cache=return_cache(m,ArrayBlock(y,touched))
   @test vcat(hcat(y[1,1:3]...),
              hcat(y[2,1],zeros(3,3),zeros(3,3)),
              hcat(y[3,1],zeros(4,3),zeros(4,3)))==evaluate!(cache,m,by)

   # Nested blocking test
   touched_parent      = Array{Bool,2}(undef,(3,3))
   touched_parent     .= false
   touched_parent[2,3] = true
   y_parent            = Array{ArrayBlock{Matrix{Float64}},2}(undef,(3,3))
   y_parent[2,3]       = by
   parent              = ArrayBlock(y_parent,touched_parent)
   cache=return_cache(m,parent)

   result_touched  = Array{Bool,2}(undef,(3,3))
   result_array    = Array{Matrix{Float64}}(undef,(3,3))
   result_touched .= touched_parent
   result_array[2,3] = vcat(hcat(y[1,1:3]...),
                            hcat(y[2,1],zeros(3,3),zeros(3,3)),
                            hcat(y[3,1],zeros(4,3),zeros(4,3)))

   @test result_array[2,3] == evaluate!(cache,m,parent).array[2,3]
end