text:
using Test
using XLSXasJSON
using JSONPointer
using OrderedCollections
using JSON

data_path = joinpath(@__DIR__, "data")

# testdata
@testset "Adobe Spry Examples" begin
    # source: https://opensource.adobe.com/Spry/samples/data_region/JSONDataSetSample.html
    f = joinpath(data_path, "example.xlsx")

    #example1
    jws = JSONWorksheet(f, "Sheet1")
    @test jws[1]["array_any"] == split("100;200;300;400", ";")
    @test jws[1]["array_int"] == [100,200,300,400]
    @test jws[1]["array_float"] == [0.1,0.2,0.3,0.4]

    @test jws[2]["array_any"] == split("500;600;700;800", ";")
    @test jws[2]["array_int"] == [500,600,700,800]
    @test jws[2]["array_float"] == [0.5,0.6,0.7,0.8]

    @test jws[3]["array_any"] == [900]
    @test jws[3]["array_string"] == ["900"]
    @test jws[3]["array_int"] == [900]
    @test jws[3]["array_float"] == [900.0]

    #example2
    jws = JSONWorksheet(f, "Sheet2")
    @test JSON.json(jws) == replace("""[{"color":"red","value":"#f00"},{"color":"green","value":"#0f0"},{"color":"blue","value":"#00f"},{"color":"cyan","value":"#0ff"},{"color":"magenta","value":"#f0f"},{"color":"yellow","value":"#ff0"},{"color":"black","value":"#000"}]""", "\n"=>"")

    #example5
    jws = JSONWorksheet(f, "Sheet3")
    @test isa(jws[1]["batters"]["batter"], Array)
    @test isa(jws[2]["batters"]["batter"], Array)
    @test isa(jws[3]["batters"]["batter"], Array)

    @test isa(jws[1]["topping"], Array)
    @test isa(jws[2]["topping"], Array)
    @test isa(jws[3]["topping"], Array)

    @test jws[1]["batters"]["batter"][1] == OrderedDict("id"=>1001, "type"=>"Regular")
    @test jws[1]["batters"]["batter"][4] == OrderedDict("id"=>1004, "type"=>"Devil's Food")

    @test jws[1]["topping"][1] == OrderedDict("id"=>5001, "type"=>"None")
    @test jws[2]["topping"][2] == OrderedDict("id"=>5002, "type"=>"Glazed")
    @test jws[3]["topping"][3] == OrderedDict("id"=>5003, "type"=>"Chocolate")

    #example6 - xf_coloriented
    f = joinpath(data_path, "example_coloriented.xlsx")
    jws = JSONWorksheet(f, "Sheet1"; row_oriented = false)
    @test jws[1]["id"] == 1
    @test jws[1]["type"] == "donut"
    @test jws[1]["name"] == "Cake"
    @test jws[1]["image"]["url"] == "images/0001.jpg"
    @test jws[1]["image"]["width"] == 200
    @test jws[1]["image"]["height"] == 2500
    @test jws[1]["thumbnail"]["url"] == "images/thumbnails/0001.jpg"
    @test jws[1]["thumbnail"]["width"] == 32
    @test jws[1]["thumbnail"]["height"] == 32
end

@testset "JSONWorkbook - deleteat!" begin
    xf = joinpath(data_path, "othercase.xlsx")
    jwb = JSONWorkbook(xf)
    @test length(jwb) == 6
    deleteat!(jwb, 1)
    @test length(jwb) == 5

    deleteat!(jwb, :promotion)
    @test length(jwb) == 4
    @test_throws ArgumentError jwb[:promotion]
end

@testset "JSONWorkbook - etc" begin
    xf = joinpath(data_path, "othercase.xlsx")
    jwb = JSONWorkbook(xf)

    @test hassheet(jwb, "Missing")
    @test hassheet(jwb, :Sheet1)
    @test !hassheet(jwb, "Sheet2")
    @test !hassheet(jwb, :Sheet3)

    @test XLSXasJSON.getsheet(jwb, :Sheet1) == jwb["Sheet1"]

    @test jwb[5] == jwb["mergeB"]
    @test jwb[end] == jwb["mergeC"]

    @test jwb[3:5] == [jwb[3],jwb[4],jwb[5]]

    names = sheetnames(jwb)
    for (i, jws) in enumerate(jwb)
        @test sheetnames(jws) == names[i]
    end

end

@testset "JSONWorkbook - writer" begin
    f = joinpath(data_path, "example.xlsx")
    jwb = JSONWorkbook(f)

    XLSXasJSON.write(data_path, jwb)

    prefix = split(basename(f), ".")[1]
    for s in sheetnames(jwb)
        @test isfile(joinpath(data_path, "$(prefix)_$(s).json"))
    end
    @test JSON.parse(JSON.json(jwb[1], 2)) == JSON.parse(JSON.json(jwb[1]))
    @test JSON.parse(JSON.json(jwb[2], 2)) == JSON.parse(JSON.json(jwb[2]))

    # wirte to XLSX
    f2 = joinpath(data_path, "example2.xlsx")
    XLSXasJSON.write_xlsx(f2, jwb)
    jwb2 = JSONWorkbook(f2)

    @test sheetnames(jwb) == sheetnames(jwb2)
    # TODO: Reconstruct pointer string from JSONPointer.Pointer
    # Type info in not preserved
    # for i in 1:length(jwb)
    #     @test jwb[i].pointer == jwb2[i].pointer
    # end
end

@testset "JSONWorksheet - merge" begin
    xf = joinpath(data_path, "othercase.xlsx")
    jwb = JSONWorkbook(xf)

    ws1 = jwb[:mergeA]
    ws2 = jwb[:mergeB]

    @test_throws AssertionError merge(ws1, ws2, "/Something")

    new_sheet = merge(ws1, ws2, "/Key")
    @test collect(keys(new_sheet[1])) == ["Key", "Address", "Name", "Property"]
    @test_throws KeyError jwb[:mergeA][1]["Property"][1]["A"]

    jwb[:mergeA] = new_sheet
    @test keys(jwb[:mergeA][1]) == keys(new_sheet[1])

    @test jwb[:mergeA][1]["Address"]["State"] == "Some"
    @test jwb[:mergeA][1]["Address"]["TEL"] == [555,1111,2222]
    @test jwb[:mergeA][1]["Property"][1]["A"] == "Out"
    @test jwb[:mergeA][1]["Property"][2]["A"] == "think"

    ws3 = jwb[:mergeC]
    new_sheet = merge(ws1, ws3, "/Key")

    @test length(new_sheet) == 6
    @test ws1[1]["Address"]["TEL"] == [555,1111,2222]
    @test new_sheet[1]["Key"] == ws1[1]["Key"]
    @test new_sheet[2]["Address"] == ws1[2]["Address"]

    @test new_sheet[1]["Address"]["TEL"] == ws3[2]["Address"]["TEL"]
    @test new_sheet[3]["Address"]["TEL"] == ws3[1]["Address"]["TEL"]

end

@testset "JSONWorksheet - append!" begin
    xf = joinpath(data_path, "append.xlsx")
    jwb = JSONWorkbook(xf)

    ws1 = jwb["Sheet1"]
    ws2 = jwb["Sheet2"]
    ws3 = jwb["Sheet3"]

    @test_throws AssertionError append!(ws1, ws3)

    @test length(ws1) == 1
    @test length(ws2) == 2

    append!(ws1, ws2)
    @test length(ws1) == 3
    @test length(ws2) == 2

    @test ws1[2] == ws2[1]
    @test ws1[3] == ws2[2]
end

@testset "JSONWorksheet - squeeze" begin
    col1 = rand(100)
    col2 = map(i -> join(rand(20), ";"), 1:100)

    data = ["/a/b/1" "/a/c::array{number}"; col1 col2]

    jws = JSONWorksheet("foo.xlsx", "Sheet1", data,; squeeze = true)
    @test length(jws) == 1
    @test jws[1]["a"]["b"][1] == data[2:end, 1]
    @test length(jws[1]["a"]["c"]) == 100
end

@testset "JSONWorksheet - getindex with index" begin
    data = ["/a" "/b" "/c::array" "/d/1/5/b";
            1     "a"  "A;100;B"  "new"
            2     "b"  "C;200;D"  "test"]

    jws = JSONWorksheet("foo.xlsx", "Sheet1", data)

    @test jws[1, 1] == 1
    @test jws[2, 1] == 2
    @test jws[1, 2] == "a"
    @test jws[2, 2] == "b"
    @test jws[1, 3] == ["A", "100", "B"]
    @test jws[2, 3] == ["C", "200", "D"]

    @test jws[1:2, 1] == [1, 2]
    @test jws[1:end, 1] == [1, 2]

    @test jws[1, 1:2] == [1 "a"]
    @test jws[1, 1:3] == permutedims([1,  "a",  ["A", "100", "B"]])
    @test jws[1, 1:end] == jws[1, 1:4]
    @test jws[:, :] == jws[1:2, 1:4] == jws[1:end, 1:end]
    @test size(jws[:]) == size(jws.data)

    @test jws[1:2, 1:2] == data[2:3, 1:2]

    @test_throws BoundsError jws[3, 1]
    @test_throws BoundsError jws[1, 5]
end

@testset "JSONWorksheet - haskey with a pointer" begin
    data = ["/a/b" "/a/c/1" "/a/d/f" "/a/c/2::array";
    1     "a"      4       "A;100;B"
    2     "b"     "k"      "C;200;D"]

    jws = JSONWorksheet("foo.xlsx", "Sheet1", data)

    @test haskey(jws, j"/a")
    @test haskey(jws, j"/a/b")
    @test haskey(jws, j"/a/c/1")
    @test haskey(jws, j"/a/d")
    @test haskey(jws, j"/a/d/f")
    @test haskey(jws, j"/a/c/2")

    @test !haskey(jws, j"/x")
    @test !haskey(jws, j"/a/1")
    @test !haskey(jws, j"/a/c/5")
    @test !haskey(jws, j"/a/d/f/k")
end

@testset "JSONWorksheet - getindex with a pointer" begin

    data = ["/a/b" "/a/c/1" "/a/d/f" "/a/c/2::array";
                1     "a"      4       "A;100;B"
                2     "b"     "k"      "C;200;D"]

    jws = JSONWorksheet("foo.xlsx", "Sheet1", data)
    x = jws[1, j"/a"]
    @test x["b"] == 1
    @test x["c"] == ["a", ["A", "100", "B"]]
    @test x["d"] == OrderedDict("f" => 4)

    @test jws[1, j"/a/c"] == ["a", ["A", "100", "B"]]

    @test jws[:, j"/a/b"] == [1, 2]

    @test haskey(jws, j"/a/b")
    @test haskey(jws, j"/a")
    @test haskey(jws, j"/a/c/2")

    @test haskey(jws, j"/a/1") == false
    @test haskey(jws, j"/aa") == false
    @test haskey(jws, j"/a/c/3") == false
end

@testset "JSONWorksheet - setindex!" begin
    data = ["/a" "/b" "/c::array";
            1     "a"  "A;100;B"
            2     "b"  "C;200;D"]
    jws = JSONWorksheet("foo.xlsx", "Sheet1", data)

    @test_throws ArgumentError jws[j"/d"] = [1]

    @test !haskey(jws, j"/d")
    jws[j"/d"] = [100, 200]
    @test jws[:, j"/d"] == [100, 200]
    @test haskey(jws, j"/d")

    jws[j"/a"] = ["a", "b"]
    @test jws[:, j"/a"] == ["a", "b"]

    jws[2, j"/a"] = ["change", "world"]
    @test jws[2, j"/a"] == ["change", "world"]

    jws[end, j"/b"] = "hooray"
    @test jws[end, j"/b"] == "hooray"

    @test_throws ErrorException jws[j"/c::array"] = [1, 2]
end

@testset "JSONWorksheet - etc" begin
    f = joinpath(data_path, "example.xlsx")

    #example1
    jws = JSONWorksheet(f, "Sheet1")

    @test xlsxpath(jws) == f
    @test firstindex(jws) == 1
    @test first(jws) == jws[1]
    @test last(jws) == jws[end]

    jws = JSONWorksheet(f, "Sheet2")
    sort!(jws, j"/color")
    @test jws[1][j"/color"] == "black"

end

@testset "Asserts" begin
    xf = joinpath(data_path, "assert.xlsx")
    @test_throws AssertionError JSONWorksheet(xf, "dup")
    @test_throws AssertionError JSONWorksheet(xf, "dup2")
    @test_throws AssertionError JSONWorksheet(xf, "dup3")

    @test_throws Exception JSONWorksheet(xf, "dict_array")
    @test_throws Exception JSONWorksheet(xf, "array_dict")

    @test_throws AssertionError JSONWorksheet(xf, "start_line")
    @test JSONWorksheet(xf, "start_line";start_line=2) isa JSONWorksheet
    @test_throws AssertionError JSONWorksheet(xf, "empty")
end

@testset "missingdata" begin
    xf = joinpath(data_path, "othercase.xlsx")
    jws = JSONWorksheet(xf, "Missing")

    @test size(jws) == (5, 4)
    @test ismissing(jws[4]["Data"]["A"])
    @test all(broadcast(el -> ismissing(el["AllNull"]), jws))
    @test collect(keys(jws[1])) == ["Key", "Data", "AllNull"]
end

@testset "Static Type" begin
    xf = joinpath(data_path, "othercase.xlsx")
    jws = JSONWorksheet(xf, "promotion")
    @test isa(jws[1]["t1"]["A"], Integer)
    @test isa(jws[1]["t1"]["B"], Bool)

    @test isa(jws[1]["t2"]["A"], Integer)
    @test isa(jws[1]["t2"]["B"], Float64)

    @test isa(jws[1]["t3"]["A"], Integer)
    @test isa(jws[1]["t3"]["B"], Bool)
    @test isa(jws[1]["t3"]["C"], Float64)

    data = ["/a::number" "/b::number" "/c::array{number}";
            1.     10   "1.5;2.2"
            2.     20   "3;55"]
    jws = JSONWorksheet("foo.xlsx", "Sheet1", data)

    @test jws[:, j"/a"] == [1, 2]
    @test jws[:, j"/b"] == [10, 20]
    @test jws[:, j"/c"] == [[1.5,2.2], [3.0,55.0]]
end

@testset "Deliminator for a Array in a cell" begin
    data = ["/a::array{number}" "/b::array{integer}" "/c::array";
            "1;2;3"     "4;5;6"   "abc;가나다"
            "1,2,3"     "4,5,6"   "abc,가나다,"]

    jws = JSONWorksheet("foo.xlsx", "Sheet1", data; delim = r";|,")

    @test  jws[1]["a"] == [1.0, 2.0, 3.0]
    @test  jws[2]["a"] == [1.0, 2.0, 3.0]
    @test  jws[1]["b"] == [4, 5, 6]
    @test  jws[2]["b"] == [4, 5, 6]
    @test  jws[1]["c"] == ["abc", "가나다"]
    @test  jws[2]["c"] == ["abc", "가나다"]

end

max_issues_repo_name: devsisters/XLSXasJSON

max_forks_repo_name: devsisters/XLSXasJSON


text:
using AtomicData
using Test
using Unitful

@testset "Units" begin
    @test isa(u"Ha", Unitful.FreeUnits)
    @test isa(u"Ry", Unitful.FreeUnits)
end

@testset "Simple queries" begin
    for unit in [u"cm^-1", u"eV", u"Ry", u"Ha"]
        for el in ["He I", "He II", "W I", "Xe I", "Xe II"]
            df = get_nist_data(el, unit)
            # Slightly bogus checks
            @test size(df,1) > 0
            @test size(df,2) > 0
        end
    end
end

ext: jl


text:
const CONTAINER_ID_BITMASK = 0xffffffffff
const CONTAINER_ID_SPLITTER = '_'
const CONTAINER_PREFIX = "container"
const EPOCH_PREFIX = "e"

parse_container_id() = parse_container_id(ENV["CONTAINER_ID"])

function parse_container_id(cidstr::String)
    parts = split(cidstr, CONTAINER_ID_SPLITTER)
    partidx = 1
    (parts[partidx] == CONTAINER_PREFIX) || error("Invalid ContainerId prefix: $cidstr")
    partidx += 1
    epoch_or_ts = parts[partidx]
    epoch = Int64(0)
    appid = Int32(0)
    clusterts = Int64(0)
    if startswith(epoch_or_ts, EPOCH_PREFIX)
        epoch = parse(Int64, epoch_or_ts[2:end])
        partidx += 1
        clusterts = parse(Int64, parts[partidx])
        partidx += 1
        appid = parse(Int32, parts[partidx])
    else
        clusterts = parse(Int64, epoch_or_ts)
        partidx += 1
        appid = parse(Int32, parts[partidx])
    end
    partidx += 1
    attemptid = parse(Int32, parts[partidx])
    partidx += 1
    cid = parse(Int64, parts[partidx])
    cid = (epoch << 40) | cid

    appid_proto = ApplicationIdProto(; id=appid, cluster_timestamp=clusterts)
    attemptid_proto = ApplicationAttemptIdProto(; application_id=appid_proto, attemptId=attemptid)
    ContainerIdProto(; app_id=appid_proto, app_attempt_id=attemptid_proto, id=cid)
end

function container_id_string(cid::ContainerIdProto)
    app_id = hasproperty(cid, :app_id) ? cid.app_id : cid.app_attempt_id.application_id
    attempt_id = hasproperty(cid, :app_attempt_id) ? cid.app_attempt_id.attemptId : 0
    id = cid.id
    epoch = id >> 40
    id = CONTAINER_ID_BITMASK & id

    parts = [CONTAINER_PREFIX]
    (epoch > 0) && push!(parts, EPOCH_PREFIX * lpad(epoch, 2, "0"))
    push!(parts, string(app_id.cluster_timestamp))
    push!(parts, lpad(app_id.id, 4, "0"))
    push!(parts, lpad(attempt_id, 2, "0"))
    push!(parts, lpad(id, 6, "0"))
    join(parts, CONTAINER_ID_SPLITTER)
end


text:
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
parr = [
    ["Q", "u[1]", 12, 12, 0, 12],
    ["Q", "u[2]", 0, 0, 0, 0],
    ["Q", :η, 12, 12, 0, 12],
    ["Q", :θ, 15, 15, 15, 15],
    ["s_aux", :y, 15, 15, 15, 15],
    ["s_aux", :w, 12, 12, 0, 12],
    ["s_aux", :pkin, 15, 15, 15, 15],
    ["s_aux", :wz0, 12, 12, 0, 12],
    ["aux", "uᵈ[1]", 15, 15, 15, 15],
    ["aux", "uᵈ[2]", 15, 15, 15, 15],
    ["aux", "ΔGᵘ[1]", 15, 15, 15, 15],
    ["aux", "ΔGᵘ[2]", 15, 15, 15, 15],
]

### fully explicit
explicit = [
    [
        "state",
        "u[1]",
        -9.58544066049463849843e-01,
        9.58544066049465071089e-01,
        -6.13908923696726568442e-17,
        4.45400263687296238402e-01,
    ],
    [
        "state",
        "u[2]",
        -4.33260855197780914020e-14,
        1.99397359064900580713e-14,
        -1.30343101521973575132e-16,
        2.95525689323845820606e-15,
    ],
    [
        "state",
        :η,
        -8.52732886154656810618e-01,
        8.52845586939211197652e-01,
        2.20052243093959998331e-14,
        6.02992088522925295813e-01,
    ],
    [
        "state",
        :θ,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        :y,
        0.00000000000000000000e+00,
        1.00000000000000011642e+06,
        5.00000000000000000000e+05,
        2.92775877460665535182e+05,
    ],
    [
        "aux",
        :w,
        -4.04553460063758398447e-04,
        4.04714358463272711169e-04,
        4.75730566051879549438e-19,
        1.63958655681888576441e-04,
    ],
    [
        "aux",
        :pkin,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        :wz0,
        -2.01164684799271339293e-04,
        2.01041968159484089494e-04,
        -2.10942374678779754294e-20,
        1.42228420244455277133e-04,
    ],
    [
        "aux",
        "uᵈ[1]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        "uᵈ[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        "ΔGᵘ[1]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        "ΔGᵘ[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

imex = [
    [
        "state",
        "u[1]",
        -9.57705359685807500192e-01,
        9.57705359685806500991e-01,
        4.54747350886464140750e-17,
        4.45323727947873004851e-01,
    ],
    [
        "state",
        "u[2]",
        -4.07596731389788922865e-14,
        2.53992382675900717845e-14,
        9.88087105439151860723e-18,
        2.63742927063789745855e-15,
    ],
    [
        "state",
        :η,
        -8.55941716778505390373e-01,
        8.56014908798837681481e-01,
        2.16732587432488803528e-14,
        6.05260814296588400829e-01,
    ],
    [
        "state",
        :θ,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        :y,
        0.00000000000000000000e+00,
        1.00000000000000011642e+06,
        5.00000000000000000000e+05,
        2.92775877460665535182e+05,
    ],
    [
        "aux",
        :w,
        -4.03419301824637895407e-04,
        4.03626959596159180614e-04,
        -3.33066907387546970565e-21,
        1.63816237485524198066e-04,
    ],
    [
        "aux",
        :pkin,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        :wz0,
        -1.97016630922629948884e-04,
        1.97066086934914155761e-04,
        -7.49400541621980656312e-19,
        1.39393318574472925390e-04,
    ],
    [
        "aux",
        "uᵈ[1]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        "uᵈ[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        "ΔGᵘ[1]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "aux",
        "ΔGᵘ[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

refVals = (explicit = (explicit, parr), imex = (imex, parr))

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using coolFFT
using Test

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


text:
using GridDensities
using Test

function sum_pdf(d::GridDensity)
	# Take a sample from each bin and ensure that they sum to one
	dims = length(d.bins)
	points = length(d.linear)
	centers = [[d.lo[i] + ((j-1)+0.5)*(d.hi[i]-d.lo[i])/d.bins[i] for j=1:d.bins[i]] for i = 1:dims]
	probs = [pdf(d, [centers[j][d.cart[i][j]] for j = 1:dims]) for i = 1:points]
	return sum(probs)*d.volume
end

function checkbounds(lo, hi, point)
	return all(lo .≤ point .≤ hi)
end

function test_uniform(lo, hi, bins)
	d = GridDensity(ones(prod(bins)), lo, hi, bins)
	point = (hi .+ lo) ./ 2
	prob = 1/(prod(bins)*d.volume)
	return abs(pdf(d, point) - prob) < 1e-8
end

@testset "GridDensities.jl" begin
    # Check for valid pdfs
    d = GridDensity(rand(8), [0.0, 0.0], [2.0, 4.0], [2, 4])
    @test abs(1.0 - sum_pdf(d)) < 1e-8
    d = GridDensity(rand(20*31*52), [0.1, 0.2, 12.0], [1.2, 1.6, 22.5], [20, 31, 52])
    @test abs(1.0 - sum_pdf(d)) < 1e-8

    # Check for samples in ranges
    d = GridDensity(rand(8), [0.0, 0.0], [2.0, 4.0], [2, 4])
    @test all([checkbounds([0.0, 0.0], [2.0, 4.0], rand(d)) for i = 1:50])
    d = GridDensity(rand(20*31*52), [0.1, 0.2, 12.0], [1.2, 1.6, 22.5], [20, 31, 52])
    @test all([checkbounds([0.1, 0.2, 12.0], [1.2, 1.6, 22.5], rand(d)) for i = 1:50])

    # Check pdf for uniform distributions
    @test test_uniform([0.0, 0.0], [2.0, 4.0], [2, 4])
    @test test_uniform([0.1, 0.2, 12.0], [1.2, 1.6, 22.5], [20, 31, 52])
end


text:
using Test
using EllZeroTrendFiltering

include("auxilliary_test_fcns.jl")

# Consider a subset of tests that are relatively small
for problem_fcn in ["discontinuous1",
                    "discontinuous2",
                    "super_exponential1",
                    "super_exponential3",
                    "white_noise",
                    "exponential"]

    include(joinpath(dirname(@__FILE__),"problems", problem_fcn * ".jl"))
    local g, I_sols  # TODO, can remove in julia 1.0?
    g, ζ_vec, I_sols, f_sols = @eval $(Symbol(problem_fcn))()

    M = length(I_sols)
    I_vec, _, f_vec = brute_force_multi(g, M)

    @testset "Data set: $problem_fcn, brute_force_search m=$m" for m in 1:length(f_sols)
        @test f_vec[m] ≈ f_sols[m]  atol=1e-8 # FIXME: The old value 1e-10 didn't work anymore
        if !isempty(I_sols[m])
            @test I_vec[m] == I_sols[m]
        end
    end
end


text:
using Base.Test
import Utils

for rep = 1:100
    @test 40. ≤ Utils.randin(40., 50.) ≤ 50.
    @test sum(Utils.randsimplex(10)) == 1
    @test all(Utils.randsimplex(10) .≥ 0)
end


text:
using PackageCompiler
create_sysimage([:DataFrames,:Distributions,:CairoMakie,:StatsFuns,:Weave],
 sysimage_path="JuliaSysimage.dll", precompile_execution_file="to_compile.jl")


text:
using LinearAlgebra

# função para fazer a fazer a fatoração QR
function fatoracao_QR(A)
    m, n = size(A)
    if m < n
        error("m deve ser ≥ n")
    end
    V = zeros(m,n)
    Q = zeros(m,n)
    R = zeros(n,n)
    # ortogonalizando todos os vetores e armazenando em V
    # v₁ = a₁
    V[:,1] = A[:,1]
    for j in 2:n
        V[:,j] = A[:,j]
        for i in 1:j-1
            # calculando o coeficiente
            αᵢⱼ = dot(V[:,i], A[:,j]) / dot(V[:,i], V[:,i])
            # ortogonalizando o vetor A[:,j]
            V[:,j] = V[:,j] - αᵢⱼ * V[:,i]
        end
    end
    # normalizando os vetores e armazenando em Q
    for j in 1:n
        norma = sqrt(dot(V[:,j], V[:,j]))
        Q[:,j] = V[:,j] / norma
        # preenchendo R
        R[j,j] = norma
        for i in 1:j-1
            αᵢⱼ = dot(V[:,i], A[:,j]) / dot(V[:,i], V[:,i])
            R[i,j] = αᵢⱼ * sqrt(dot(V[:,i], V[:,i]))
        end
    end
    return Q, R
end

max_stars_repo_name: EleonoraAvello/GrupoBetaPacote

max_issues_repo_name: EleonoraAvello/GrupoBetaPacote

max_forks_repo_name: EleonoraAvello/GrupoBetaPacote


text:
@testset "Simulation Unit Tests:" begin

    s = Simulation(test_model, test_driver)

    @test isa(trace(s), BlockOptTrace)

    @test isa(backend(s), BlockOptBackend)

    @testset "inhereted methods" begin

        @test isa(trs_timer(s), EvaluationTimer)

        @test Δt(trs_timer(s)) ≡ 0.0

        @test isa(trs_counter(s), EvaluationCounter)

        @test evaluations(trs_counter(s)) ≡ 0

        @test isa(ghs_timer(s), EvaluationTimer)

        @test Δt(ghs_timer(s)) ≡ 0.0

        @test isa(ghs_counter(s), EvaluationCounter)

        @test evaluations(ghs_counter(s)) ≡ 0

        @test weave!(s, :f_vals, 0.0) ≈ [0.0]

        @test weave!(s, :∇f_norms, 1.0) ≈ [1.0]

        @test weave!(s, :Δ_vals, 2.0) ≈ [2.0]

        @test weave!(s, :p_norms, 3.0) ≈ [3.0]

        @test weave!(s, :ρ_vals, 4.0) ≈ [4.0]

        @test pop!(f_vals(s)) ≈ 0.0

        @test pop!(∇f_norms(s)) ≈ 1.0

        @test pop!(Δ_vals(s)) ≈ 2.0

        @test pop!(p_norms(s)) ≈ 3.0

        @test pop!(ρ_vals(s)) ≈ 4.0

        @test weave!(s, f_vals, 0.0) ≈ [0.0]

        @test weave!(s, ∇f_norms, 1.0) ≈ [1.0]

        @test weave!(s, Δ_vals, 2.0) ≈ [2.0]

        @test weave!(s, p_norms, 3.0) ≈ [3.0]

        @test weave!(s, ρ_vals, 4.0) ≈ [4.0]

        @test weave_level(s) ≡ weave_level(test_driver)

        @test fₖ(s) ≈ ∇fₖ_norm(s) ≈ pₖ_norm(s) ≈ Δₖ(s) ≈ ρ(s) ≈ 0.0

        @test log_level(s) ≡ log_level(test_driver)

        @test isfile(io(s))

        info!(s, "This is an info! message.")

        debug!(s, "This is an debug! message.")

        warn!(s, "This is an warn! message.")

        error!(s, "This is an error! message.")

        # Currently, the statements below print to console. Fix.

        # info!(s,  "Simulation initilaized by:\\n", test_model, test_driver)

        # info!(s,  "Simulation initilized trace and backend:\\n", trace(s), backend(s))

        restricted_type_test(s)
    end

    @testset "edge case testing" begin

    end

    nothing
end

@testset "Simulation Behavior Tests:" begin

    max_iterations!(test_driver, 4)

    s = Simulation(test_model, test_driver)

    @testset "initialize subroutine" begin

        @test evaluations(ghs_counter(s)) ≡ 0

        initialize(s)

        @test evaluations(ghs_counter(s)) ≡ 1
    end

    i = 0

    while !terminal(s)

        build_trs(s)

        @testset "solve_trs subroutine" begin

            k = evaluations(trs_counter(s))

            elapsed = Δt(trs_timer(s))

            solve_trs(s)

            @test evaluations(trs_counter(s)) ≡ k + 1

            @test Δt(trs_timer(s)) > elapsed
        end

        build_trial(s)

        successful_trial = accept_trial(s)

        @testset "accept_trial subroutine" begin

            if successful_trial

                @test last(Δ_vals(s)) ≈ Δₖ(s)

                @test last(ρ_vals(s)) ≈ ρ(s)

                @test last(p_norms(s)) ≈ pₖ_norm(s)
            end

            @test last(f_vals(s)) ≈ fₖ(s) # always holds

            @test last(∇f_norms(s)) ≈ ∇fₖ_norm(s) # always holds
        end

        successful_trial && pflag(s) && secantQN(s)

        successful_trial && update_Sₖ(s)

        successful_trial && @testset "gHS subroutine" begin

            cost = evaluations(ghs_counter(s))

            elapsed = Δt(ghs_timer(s))

            ∇f_norm = ∇fₖ_norm(s)

            gHS(s)

            @test !(last(∇f_norms(s)) ≈ ∇f_norm)

            @test evaluations(ghs_counter(s)) ≡ cost + 1

            @test Δt(ghs_timer(s)) > elapsed
        end

        successful_trial && blockQN(s)

        update_Δₖ(s)

        i += 1
    end

    @test last(∇f_norms(s)) ≈ ∇fₖ_norm(s)

    @test last(f_vals(s)) ≈ fₖ(s)

    if max_iterations(backend(s)) > evaluations(trs_counter(s)) # implying convergence

        @test last(Δ_vals(s)) ≈ Δₖ(s)

        @test last(ρ_vals(s)) ≈ ρ(s)

        @test last(p_norms(s)) ≈ pₖ_norm(s)

        @test ∇fₖ_norm(s) < ϵ_tol(backend(s))
    end

    return nothing
end


text:
abstract type TermColorDepth end
struct TermColor256   <: TermColorDepth end
struct TermColor24bit <: TermColorDepth end

"""
    colorant2ansi(color::Colorant) -> Int

Converts the given colorant into an integer index that corresponds
to the closest 256-colors ANSI code.

```julia
julia> colorant2ansi(RGB(1., 1., 0.))
226
```

This function also tries to make good use of the additional number
of available shades of gray (ANSI codes 232 to 255).

```julia
julia> colorant2ansi(RGB(.5, .5, .5))
244

julia> colorant2ansi(Gray(.5))
244
```
"""
colorant2ansi(color) = _colorant2ansi(color, TermColor256())

# Fallback for non-rgb and transparent colors (convert to rgb)
_colorant2ansi(gr::Color, colordepth::TermColorDepth) =
    _colorant2ansi(convert(RGB, gr), colordepth)
_colorant2ansi(gr::TransparentColor, colordepth::TermColorDepth) =
    _colorant2ansi(color(gr), colordepth)

# 256 colors
function _colorant2ansi(col::AbstractRGB, ::TermColor256)
    r, g, b = clamp01nan(red(col)), clamp01nan(green(col)), clamp01nan(blue(col))
    r24, g24, b24 = map(c->round(Int, c * 23), (r, g, b))
    if r24 == g24 == b24
        # Use grayscales because of higher resultion
        # This way even grayscale RGB images look good.
        232 + r24
    else
        r6, g6, b6 = map(c->round(Int, c * 5), (r, g, b))
        16 + 36 * r6 + 6 * g6 + b6
    end
end

_colorant2ansi(gr::Color{<:Any,1}, ::TermColor256) = round(Int, 232 + clamp01nan(real(gr)) * 23)

# 24 bit colors
function _colorant2ansi(col::AbstractRGB, ::TermColor24bit)
    r, g, b = clamp01nan(red(col)), clamp01nan(green(col)), clamp01nan(blue(col))
    map(c->round(Int, c * 255), (r, g, b))
end

function _colorant2ansi(gr::Color{<:Any,1}, ::TermColor24bit)
    r = round(Int, clamp01nan(real(gr)) * 255)
    r, r, r
end


text:
using DirichletMixtureModels

importall DirichletMixtureModels

import Distributions:
        UnivariateDistribution,
        MultivariateDistribution,
        Exponential,
        Normal,
        NormalCanon,
        MvNormal,
        NormalStats,
        rand,
        suffstats,
        pdf

import ConjugatePriors:
    NormalGamma,
    NormalWishart,
    posterior_canon,
    rand

function example_pdf(y::Float64,μ,λ,μ0,n0,α0,β0)
  pdf(NormalCanon(λ*μ, λ), y)
end
function example_pdf(y::Array{Float64},μ,λ,μ0,n0,α0,β0)
  pdf.(NormalCanon(λ*μ, λ), y)
end
function example_post(Y::Array{Float64,1}, μ0,n0,α0,β0)
  p=posterior_canon(NormalGamma(μ0,n0,α0,β0),suffstats(Normal,Y))
  rand(p)
end
function example_post(y::Float64, μ0,n0,α0,β0)
  p=posterior_canon(NormalGamma(μ0,n0,α0,β0),suffstats(Normal,[y]))
  rand(p)
end
function example_pri(μ0,n0,α0,β0)
  rand(NormalGamma(μ0,n0,α0,β0))
end
function example_marg(y::Float64,μ0,n0,α0,β0)
  gamma(α0+1./2)/gamma(α0) * sqrt(n0/(n0+1)) * 1/sqrt(2*π) * β0^α0 /
    (β0+n0/2/(n0+1)*(y-μ0)^2)^(α0+1./2)
end

max_stars_repo_name: nsdumont/jDirichletMixture

max_issues_repo_name: nsdumont/jDirichletMixture

max_forks_repo_name: nsdumont/jDirichletMixture


text:
using JSON

function parse_notebook(document_body)
    nb = JSON.parse(document_body)
    code_no = 0
    doc_no = 0

    # TODO: handle some of options ?
    options = Dict{Symbol,Any}()
    opt_string = ""

    chunks = map(nb["cells"]) do cell
        text = string('\n', join(cell["source"]), '\n')
        return if cell["cell_type"] == "code"
            CodeChunk(text, code_no += 1, 0, opt_string, options)
        else
            DocChunk(text, doc_no += 1, 0; notebook = true)
        end
    end

    return Dict(), chunks
end


text:
# Gold heater on a substrate
# SI units throughout the script

# Top layer down to substrate information arrays

# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

function main()
    # Half-width of heater line [m]
    b = (12.5/2)*1e-6
    # Length of heater line [m]
    l = 1.0e-3
    # Range of frequencies [Hz]
    f = exp10.(LinRange(0, 9, 1500))
    # Power [W]
    p = 0.030^2*22.11
    # Heater thermal resistance
    ρh = [0. 0.]
    # Interface thermal resistances
    thresistances = [0.]
    # Wrap them into structures
    layers = [
        LayerTOM(310.0, 1.0, 0.2e-6, 2.441e6), # heater
        LayerTOM(160.0, 1.0, 525.0e-6, 2320*700.), # substrate
    ]
    hgeometry = HeaterGeometry(b, l, ρh)
    source = Source(p, f)
    # Solve
    sol = three_omega(layers, hgeometry, source, thresistances; int_limit=5.0e6)
    return sol
end

sol = main()

plot(TOMPlot(), sol)
gui()


text:
##
cd(@__DIR__())
using IterTools: groupby
using Combinatorics: powerset
jolts = map(eachline("input.txt")) do line
    parse(Int, line)
end;

##
function find_chain(jolts)
    res = [0]
    joltsc = copy(jolts)
    while !isempty(joltsc)
        x = last(res)
        ms = findall(joltsc) do y
            abs(y-x) <= 3
        end
        if isnothing(ms) || isempty(ms)
            @info "Not found $x"
            break
        else
            m = argmin([joltsc[m] for m ∈ ms])
        end
        append!(res, joltsc[ms[m]])
        deleteat!(joltsc, ms[m])
    end
    append!(res, last(res) + 3)
    res, diff(res), [filter(==(x), diff(res)) |> length for x ∈ (1,2,3)]
end

##
chain, d, h = find_chain(jolts)
@info h[1] * h[3]

## PARTE 2
function check_chain(chain)
    maximum(diff(chain)) <= 3
end
function find_pivot(chain)
    Iterators.filter(eachindex(chain)) do i
        i == 1 && return false
        i == lastindex(chain) && return false
        chain[i+1] - chain[i-1] <= 3
    end |> collect
end
function count_chains(chain, ns)
    gs = []
    cur = Int[ns[1]]
    for n ∈ ns[2:end]
        if n - last(cur) == 1
            push!(cur, n)
        else
            push!(gs, cur)
            cur = [n]
        end
    end
    push!(gs, cur)
    prod(gs) do g
        length(g) == 1 && return 2
        ss = powerset(g, 2)
        1 + length(g) + count(ss) do s
            c′ = copy(chain)
            check_chain(deleteat!(c′, s))
        end
    end
end;

##
ns = find_pivot(chain)
count_chains(chain, ns)


text:
@enum Building Smelter Constructor Assembler Manufacturer Refinery Miner Foundry WaterExtractor OilExtractor Packager

struct Recipe
    name::String
    out::Vector{Tuple{DataType, Float64}}
    in::Vector{Tuple{DataType, Float64}}
    building::Building
    duration::Float64
end

dataPath = joinpath(@__DIR__, "data.json")

function readData()
    data = JSON3.read(read(dataPath))

    _buildings = [data.buildings[b] for b in first.(filter(!isempty, unique(r.producedIn for r in values(data.recipes))))]
    _recipes = filter(x -> !isempty(x.producedIn), collect(values(data.recipes)))
    _products = unique(union((p.item for r in _recipes for p in r.products), (p.item for r in _recipes for p in r.ingredients)))

    classNameToProduct = Dict{String, Type{<:Product}}(data.items[p].className => nameToType(data.items[p].name) for p in _products)
    classNameToBuilding = Dict{String, Building}(b.className => nameToType(b.name) for b in _buildings)

    empty!(allRecipes)
    empty!(dictProductDependantRecipes)
    empty!(dictProductRecipes)
    empty!(baseRecipes)

    for r in _recipes
        out = [(classNameToProduct[p.item], p.amount * 60 / r.time) for p in r.products]
        in = [(classNameToProduct[p.item], p.amount * 60 / r.time) for p in r.ingredients]
        building = classNameToBuilding[only(r.producedIn)]
        push!(allRecipes, Recipe(r.name, out, in, building, r.time))
    end

    minerMK1 = data.miners.Build_MinerMk1_C
    for p in minerMK1.allowedResources
        product = classNameToProduct[p]
        qty = minerMK1.itemsPerCycle / minerMK1.extractCycleTime * 60
        push!(allRecipes, Recipe("$product", [(product, qty)], [], Miner, 60))
    end
    oilPump = data.miners.Build_OilPump_C
    for p in oilPump.allowedResources
        product = classNameToProduct[p]
        qty = oilPump.itemsPerCycle / oilPump.extractCycleTime * 60 / 1000
        push!(allRecipes, Recipe("$product", [(product, qty)], [], OilExtractor, 60))
    end
    waterPump = data.miners.Build_WaterPump_C
    for p in waterPump.allowedResources
        product = classNameToProduct[p]
        qty = waterPump.itemsPerCycle / waterPump.extractCycleTime * 60 / 1000
        push!(allRecipes, Recipe("$product", [(product, qty)], [], WaterExtractor, 60))
    end

    for r in allRecipes
        for (p, qty) in r.out
            push!(get!(dictProductRecipes, p, Set()), (r, qty))
        end
    end
    for r in allRecipes
        for (p, qty) in r.in
            push!(get!(dictProductDependantRecipes, p, Set()), (r, qty))
        end
    end

    for p in subtypes(Product)
        if p in harvestedProducts
            dictProductRecipes[p] = Set()
        end
        if !haskey(dictProductDependantRecipes, p)
            dictProductDependantRecipes[p] = Set()
        end
        if !haskey(dictProductRecipes, p)
            @warn "No recipe found to make $p"
        end
    end

    append!(baseRecipes, filter(x -> !occursin("Alternate", x.name), allRecipes))
    nothing
end

nameToType(s::AbstractString) = getfield(Satisfactory, Symbol(replace(s, r" |-|\." => ""))) # remove spaces, dashes and dots

recipes(::Type{T}) where T <: Product = dictProductRecipes[T]
recipes(::T) where T <: Product = recipes(T)

dependantRecipes(::Type{T}) where T <: Product = dictProductDependantRecipes[T]
dependantRecipes(::T) where T <: Product = dictProductDependantRecipes(T)

import Base.show, Base.==, Base.hash
show(io::IO, r::Recipe) = print(io, "$(r.name)")
hash(r::Recipe, h::UInt) = hash(r.name, h)
==(r1::Recipe, r2::Recipe) = r1.name == r2.name


text:
using Oxygen, Test

@testset "All Tests" begin
    include("periodic.jl")
    include("electrons.jl")
    include("atom.jl")
    include("molecules.jl")
    include("smiles.jl")
    include("ops.jl")
end


text:
include("dynamics.jl")
using ReferenceFrameRotations
#dynamics.jlのテスト
#クォータニオンとRK法による回転運動のシミュレーションのテスト(振り子)

l=10 #振り子の長さ
m=10 #重りの重さ
d=1 #重りの直径

g=9.8*[0,0,-1] #重力加速度ベクトル
r0=l*[1,0,0] #振り子の初期位置ベクトル
#慣性テンソル
Ix=0.4m*(d^2)
Iy=0.4m*(d^2)+m*(l^2)
Iz=0.4m*(d^2)+m*(l^2)
I=
[Ix 00 00
;00 Iy 00
;00 00 Iz]

#トルクを求める関数
function torque((q,ω)::Tuple{SatelliteToolbox.Quaternion,Vector},(r0,m,g)::Tuple{Vector,Number,Vector})
    r=vect(q*r0*conj(q))
    v=cross(ω,r)
    f=m*g+m*cross(v,ω)
    T= conj(q)*cross(r,f)*q
    return [T.q1,T.q2,T.q3]
end

#力学的エネルギーを求める関数
function energy(q::SatelliteToolbox.Quaternion,ω::Vector,(r0,m,g)::Tuple{Vector,Number,Vector})
    r=vect(q*r0*conj(q))
    ω=vect(conj(q)*ω*q)
    v=cross(ω,r)
    fg=m*g
    E=0.5m*dot(v,v)+0.5dot(ω,[Ix 0 0;0 Ix 0;0 0 Ix]*ω),-dot(r,fg)
    return E
end

#振り子のシミュレーション
q=SatelliteToolbox.Quaternion(1.0,0.0,0.0,0.0)
ω=[0.0,-0.005pi,0.05pi]
E=energy(q,ω,(r0,m,g))
result_all=[(q,ω,E)]

len=500000
dt=0.01
for i in 2:len
    q,ω=(result_all[i-1][1],result_all[i-1][2])
    T=torque((q,ω),(r0,m,g))
    q,ω=dynamics(q,ω,T,I,dt)
    E=energy(q,ω,(r0,m,g))
    push!(result_all,(q,ω,E))
end

#結果から一列取り出す関数
function handle_result(result,len,k)
    b=[result[1][k]]
    for i in 2:len
        push!(b,result[i][k])
    end
    return b
end

q=handle_result(result_all,len,1)
ω=handle_result(result_all,len,2)
Ek=handle_result(handle_result(result_all,len,3),len,1) #運動エネルギー
Ep=handle_result(handle_result(result_all,len,3),len,2) #位置エネルギー
;

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
import Base.Broadcast: instantiate, materialize, materialize!

@inline layout(x, ::NOnes{N}) where {N} = _layout(x, buffer(x), CoordinateTransform{N}(x))
@inline function _layout(x::X, y::Y, t1::CoordinateTransform{NI,NO}) where {X,Y,NI,NO}
    b, t2 = layout(y, NOnes(static(NO)))
    return b, t2(t1)
end
# end recursion b/c no new buffer
_layout(x::X, y::X, t::ComposedTransform) where {X} = x, t
# no new buffer and unkown index transformation, s
_layout(x::X, y::X, ::UnkownTransform{N}) where {X,N} = x, IdentityTransform{N}()
# new buffer, but don't know how to transform indices properly
_layout(x::X, y::Y, ::UnkownTransform{N}) where {X,Y,N} = x, IdentityTransform{N}()

# TODO conversion to pointers
#=
tryptr(x::X) where {X} = tryptr(device(X), x)
tryptr(::CPUPointer, x) = LazyPreserve(x)
tryptr(::AbstractDevice, x) = PseudoPtr(x)  # TODO handle cartesian indexers
=#

# TODO layout
struct Layouted{B,I,T}
    buffer::B
    indices::I
    tform::T
end

@inline function instantiate(lyt::Layouted{B,I}) where {B,I}
    b, t = layout(getfield(lyt, :buffer), NOnes(index_dimsum(I)))
    return Layouted(tryptr(b), getfield(lyt, :indices), t)
end

@inline materialize(lyt::Layouted{B,I,Nothing}) where {B,I} = materialize(instantiate(lyt))
function materialize(lyt::Layouted{B,I,T}) where {B,I,T} end


text:
function compute_reynolds_number(ω, p)

     # --> Compute the velocity.
     vx, vy = compute_velocity(ω, p)

     # --> Mean-squared amplitude.
     E = mean(vx.^2 + vy.^2)

     # --> Reynolds number.
     Re = √E * p.Lx / p.ν

    return Re
end

function compute_velocity(ω, p)

    # --> Extract parameters.
    @unpack Lx, Ly, nx, ny = p

    # --> Get the wavenumbers.
    α = fftfreq(nx, nx/Lx) * 2π
    β = fftfreq(ny, ny/Ly) * 2π

    # --> Initialize arrays.
    vx, vy = zero(ω), zero(ω)

    # --> Compute the velocity component in spectral space.
    for j = 1:nx, i = 1:ny

        # --> Streamfunction.
        k² = α[j]^2 + β[i]^2
        ψ = i == j == 1 ? zero(ω[i, j]) : ω[i, j] / k²

        # --> Velocity.
        vx[i, j] = im * β[i] * ψ
        vy[i, j] = im * α[j] * ψ

    end

    # --> From spectral space to physical space.
    ifft!(vx), ifft!(vy)

    return real(vx), real(vy)
end

function compute_taylor_scale(ω, p)

    # --> Compute the velocity field.
    vx, vy = compute_velocity(ω, p)

    # --> Compute the amplitude of the velocity fluctuation.
    E = mean(vx.^2 + vy.^2)

    # --> Compute the amplitude of the vorticity fluctuation.
    Q = mean(real(ifft(ω)).^2)

    # --> Compute the Taylor micro-scale.
    λ = √(E/Q)

    return λ
end

function mesh(p)

    # --> Unpack the paramters.
    @unpack Lx, Ly, nx, ny = p

    # --> Construct the one-dimensional meshes.
    x = LinRange(-Lx/2, Lx/2, nx+1)[1:end-1]
    y = LinRange(-Ly/2, Ly/2, ny+1)[1:end-1]

    return x, y
end


text:
using Test

@testset "/src/LocalPolytope.jl" begin

using BellScenarip: LocalPolytope

@testset "permutation encodings are not necessary for communication protocols" begin

    # this test verifies that it is okay to use the communication protocols method
    # instead of the strategies method to compute communication protocols.

    # testing for all 1-bit cases with 3 in/out or fewer
    @testset "bob input $y" for y in 1:3
        @testset "bob output $b" for b in 1:3
            @testset "alice input = $x" for x in 1:3
                @testset "protocol out = $num_out" for num_out in 2:2
                    I = Math.id(y)

                    ρ_perms = map(
                        protocol -> kron(I,protocol),
                        LocalPolytope.strategies(x,num_out)
                    )
                    ρ_test = map(
                        protocol -> kron(I,protocol),
                        LocalPolytope.communication_protocols(x,num_out)
                    )

                    β_joint_set = LocalPolytope.strategies(num_out*y,b)

                    test_set = []
                    perm_set = []
                    for ρ in ρ_test
                        test_set = cat(test_set, map(
                            (β) -> β*ρ,
                            β_joint_set
                        ), dims = 1)
                    end

                    for ρ in ρ_perms
                        perm_set = cat(perm_set, map(
                            (β) -> β*ρ,
                            β_joint_set
                        ), dims = 1)
                    end

                    @test issetequal(unique(test_set), unique(perm_set))
                end
            end
        end
    end
end

end


text:
function run_thermal_bubble_regression_test(arch, grid_type)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    if grid_type == :regular
        grid = RectilinearGrid(arch, size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz), halo=(1, 1, 1))
    elseif grid_type == :vertically_unstretched
        zF = range(-Lz, 0, length=Nz+1)
        grid = RectilinearGrid(arch, size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=zF, halo=(1, 1, 1))
    end

    closure = ScalarDiffusivity(ν=4e-2, κ=4e-2)
    model = NonhydrostaticModel(grid=grid, closure=closure, coriolis=FPlane(f=1e-4),
                                buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
    simulation = Simulation(model, Δt=6, stop_iteration=10)

    model.tracers.T.data.parent .= 9.85
    model.tracers.S.data.parent .= 35.0

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    CUDA.@allowscalar model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    datadep_path = "regression_test_data/thermal_bubble_regression.nc"
    regression_data_filepath = @datadep_str datadep_path

    ####
    #### Uncomment the block below to generate regression data.
    ####

    #=
    @warn ("You are generating new data for the thermal bubble regression test.")

    outputs = Dict("v" => model.velocities.v,
                   "u" => model.velocities.u,
                   "w" => model.velocities.w,
                   "T" => model.tracers.T,
                   "S" => model.tracers.S)

    nc_writer = NetCDFOutputWriter(model, outputs, filename=regression_data_filepath, schedule=IterationInterval(10))
    push!(simulation.output_writers, nc_writer)
    =#

    ####
    #### Regression test
    ####

    run!(simulation)

    ds = Dataset(regression_data_filepath, "r")

    test_fields = CUDA.@allowscalar (u = Array(interior(model.velocities.u)),
                                     v = Array(interior(model.velocities.v)),
                                     w = Array(interior(model.velocities.w)),
                                     T = Array(interior(model.tracers.T)),
                                     S = Array(interior(model.tracers.S)))

    correct_fields = CUDA.@allowscalar (u = ds["u"][:, :, :, end],
                                        v = ds["v"][:, :, :, end],
                                        w = ds["w"][:, :, :, end],
                                        T = ds["T"][:, :, :, end],
                                        S = ds["S"][:, :, :, end])

    summarize_regression_test(test_fields, correct_fields)

    CUDA.allowscalar(true)
    @test all(test_fields.u .≈ correct_fields.u)
    @test all(test_fields.v .≈ correct_fields.v)
    @test all(test_fields.w .≈ correct_fields.w)
    @test all(test_fields.T .≈ correct_fields.T)
    @test all(test_fields.S .≈ correct_fields.S)
    CUDA.allowscalar(false)

    return nothing
end


text:
using DataDeps

register(DataDep(
  "online-news-popularity",
  "https://archive.ics.uci.edu/ml/datasets/Online+News+Popularity",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/00332/OnlineNewsPopularity.zip",
  "dba2ae526f62ccef6f2f8efb53b0268319f4d0d2719bb98d88958cfe31d58a22",
  post_fetch_method=(path -> begin
    unpack(path)
    new_path = joinpath(dirname(path), "OnlineNewsPopularity/OnlineNewsPopularity.csv")
    mv(new_path, path)

    UCIData.process_dataset(path,
        target_index=61,
        feature_indices=2:60,
        id_indices=1,
        skipto=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
module v
using JLD
using FFTW
using Trapz
using Printf
using PyPlot
using Random
using Statistics
using Parameters
using SparseArrays
using ProgressMeter
using Distributions
using LinearAlgebra
using DelimitedFiles

@with_kw mutable struct Params
   prefix::String = "D"
   record::Bool = true        # record trajectories if true
   compACF::Bool = true       # compute autocorrelation function if true

   N::Int = 1000              # number of neurons
   K::Int = 1000              # in-degree
   T::Int = 2000              # simulation time
   warmUp::Int = 10           # warm up time
   ncorr::Int = 40            # autocorrelation time
   dt::Float64 = 0.02         # time step
   τ::Float64 = 1.0           # time constant
   rcNsamp::Int = 1000;       @assert rcNsamp <= N;   # recorded sample
   rcDt::Float64 = 0.02;      @assert rcDt >= dt;     # recorded time step
   Φname::String = "rectTanh" # response function name
   λ::Float64 = 1.0           # powerlaw exponent
   Jdist::String = "gamma"    # connectivity distribution name
   ρ_e::Float64 = 0.          # fraction of excitatory neurons

   J_ee::Float64 = 0.         # Je→e (sets the mean coupling)
   J_ei::Float64 = 0.         # Ji→e
   J_ie::Float64 = 0.         # Je→i
   J_ii::Float64 = 2.0        # Ji→i
   scaleJ::Float64 = -0.5     # scaling for mean coupling

   g_ee::Float64 = 0.         # ge→e (sets the std of coupling)
   g_ei::Float64 = 0.         # gi→e
   g_ie::Float64 = 0.         # ge→i
   g_ii::Float64 = 2.0        # gi→i
   scaleg::Float64 = -0.25    # scaling for std coupling

   I_e::Float64 = 0.          # input current to e
   I_i::Float64 = 1.0         # input current to i
   scaleI::Float64 = 0.0      # scaling for input current

   seedJ::Int = 0             # network connectivity seed
   seedIC::Int = 0            # initial condition seed
   progress::Bool = true      # show progress if true

end
dp = Params()                 # default parameters

#################################################
function runsim(;p = dp)
   startTime = time()
   function rectTanh(x)
      if x > 0
         return tanh(x)
      else
         return 0.0
      end
   end
   function rectLin(x)
      if x > 0
         return x
      else
         return 0.0
      end
   end
   function Heaviside(x)
      if x > 0
         return 1.0
      else
         return 0.0
      end
   end
   function rectPower(x)
      if x > 0
         return x^p.λ
      else
         return 0.0
      end
   end
   diluteJ = (p.N != p.K)

   N_i = round(Int, p.N - p.N * p.ρ_e)
   N_e = p.N - N_i
   K_e = round(Int, p.K * p.ρ_e)
   K_i = round(Int, p.K - p.K * p.ρ_e)

   I_e = p.I_e * (K_e ^ p.scaleI)
   I_i = p.I_i * (K_i ^ p.scaleI)

   tstep = floor(Integer, p.T / p.dt)
   tstep_warmUp = floor(Integer, p.warmUp * p.τ / p.dt)
   tstep_tot = tstep + tstep_warmUp
   tcorr_step = floor(Integer, p.ncorr * p.τ / p.rcDt)
   tcorr = 0 : p.rcDt : p.ncorr * p.τ - p.rcDt
   τ⁻¹dt = p.dt / p.τ

   ΦDict = Dict("rectTanh" => rectTanh, "rectLin" => rectLin,
                "rectPower" => rectPower, "Heaviside" => Heaviside)
   Φ = ΦDict[p.Φname]

   rcEvery = floor(Int, p.rcDt/p.dt)
   rctstep = floor(Int, tstep/rcEvery)

   if (p.ρ_e != 0)
      rcN_e        = Int(p.ρ_e * p.rcNsamp)
      rcNsampE     = 1:rcN_e
      rcNsampI     = (rcN_e + 1): p.rcNsamp
      rcNsampRange = vcat(rcNsampE,rcNsampI)
   else
      rcNsampRange = 1:p.rcNsamp
   end

   D = Dict()
   F = Dict()

   if (diluteJ)
      if !iszero(p.ρ_e)
         printstyled("partially-connected EI","\n";color=:green)
         J = randJ_EI_PC(p = p)
      else
         printstyled("partially-connected pure-I","\n";color=:green)
         J = randJ_I_PC(p = p)
      end
   else
      if !iszero(p.ρ_e)
         printstyled("fully-connected EI","\n";color=:green)
         J = randJ_EI_FC(p = p)
      else
         printstyled("fully-connected pure-I","\n";color=:green)
         J = randJ_I_FC(p = p)
      end
   end
   rcCount = 0

   x       = p.g_ii * randn(MersenneTwister(p.seedIC), p.N)
   ϕ       = zeros(Float64,p.N)
   I       = zeros(Float64,p.N)
   percTon = zeros(Int,p.N)

   m   = zeros(Float64,rctstep)
   x̄   = zeros(Float64,rctstep)
   η̄   = zeros(Float64,rctstep)
   n   = zeros(Int,rctstep)
   μ   = zeros(Float64,rctstep)

   if (p.record || p.compACF)
      xs = zeros(Float64,rctstep, p.rcNsamp)
      ϕs = zeros(Float64,rctstep, p.rcNsamp)
      ηs = zeros(Float64,rctstep, p.rcNsamp)
   end

   if !iszero(p.ρ_e)
      mE = zeros(Float64,rctstep)
      mI = zeros(Float64,rctstep)
      nE = zeros(Int,rctstep)
      nI = zeros(Int,rctstep)
   end

   tmx  = zeros(Float64,p.N)
   tmη  = zeros(Float64,p.N)
   tmx2 = zeros(Float64,p.N)
   tmη2 = zeros(Float64,p.N)

   I[1:N_e] .= I_e
   I[(N_e + 1):p.N] .= I_i
   ϕ = Φ.(x)
   iOn = findall(ϕ .> 0)

   ##########################
   t = 0.0
   printParams(J, I_i, p = p)
   prog = Progress(tstep_tot, 0.5)
   for ts = 1:tstep_tot
      p.progress && ProgressMeter.update!(prog,ts);
      t₀ = ts - tstep_warmUp
      #######################
      η = J * ϕ
      x = x + τ⁻¹dt * (-x + η + I)
      ϕ = Φ.(x)
      #######################
      if (t₀ > 0) && (t₀ % rcEvery == 0)
         iOn = findall(x .> 0)
         percTon[iOn] .+= 1
         rcCount += 1

         n[rcCount]   = length(iOn)
         η̄[rcCount]   = mean(η)
         x̄[rcCount]   = mean(x)
         m[rcCount]   = mean(ϕ)
         μ[rcCount]   = mean(ϕ[iOn])
         tmx  .+= x
         tmη  .+= η
         tmx2 .+= x.^2
         tmη2 .+= η.^2

         if !iszero(p.ρ_e)
            mE[rcCount] = mean(ϕ[1:N_e])
            mI[rcCount] = mean(ϕ[N_e+1:end])
            nE[rcCount] = count(iOn .<= N_e)
            nI[rcCount] = count(iOn .> N_e)
         end

         if (p.record || p.compACF)
            xs[rcCount,:] = x[rcNsampRange]
            ϕs[rcCount,:] = ϕ[rcNsampRange]
            ηs[rcCount,:] = η[rcNsampRange]
         end
      end
   end

   #####################
   D[:m] = mean(m)                     # [⟨ϕ⟩]
   D[:x̄] = mean(x̄)                     # [⟨x⟩]
   D[:η̄] = mean(η̄)                     # [⟨η⟩]

   ####################
   D[:Mtm2x] = mean((tmx/rcCount).^2)  # [⟨x⟩²]
   D[:Mtm2η] = mean((tmη/rcCount).^2)  # [⟨η⟩²]

   D[:Mtmx2] = mean((tmx2/rcCount))    # [⟨x²⟩]
   D[:Mtmη2] = mean((tmη2/rcCount))    # [⟨η²⟩]

   D[:tvx] = D[:Mtmx2] - D[:Mtm2x]     # temporal var of x
   D[:tvη] = D[:Mtmη2] - D[:Mtm2η]     # temporal var of η

   D[:qvx] = D[:Mtm2x] - D[:x̄]^2       # quenched var of x
   D[:qvη] = D[:Mtm2η] - D[:η̄]^2       # quenched var of η

   #####################
   D[:f̄], D[:Δf²] = mean(n/p.N), var(n/p.N)
   D[:μ̄], D[:Δμ²] = mean(μ), var(μ)
   D[:percTon] = percTon/rctstep

   #####################
   D[:α] = mean(m) * sqrt(p.K) * (p.J_ii / p.g_ii)^2

   #####################
   if !iszero(p.ρ_e)
      D[:mE] = mE
      D[:mI] = mI
      D[:fE] = nE / N_e
      D[:fI] = nI / N_i
   end

   #####################
   if (p.compACF)
      Rδη = comp_ACF(ηs .- mean(ηs,dims=1), tcorr_step);D[:Rδη]=Rδη
      Rδϕ = comp_ACF(ϕs .- mean(ϕs,dims=1), tcorr_step);D[:Rδϕ]=Rδϕ
      Rδx = comp_ACF(xs .- mean(xs,dims=1), tcorr_step);D[:Rδx]=Rδx

      D[:β] = comp_β(tcorr, Rδη; p = p)
   end

   if (p.record)
      F[:ηs] = ηs
      F[:xs] = xs
      F[:ϕs] = ϕs
   else
      F[:xs] = xs[:,1:10]
      F[:ηs] = ηs[:,1:10]
   end

   D[:tcorr] = tcorr
   D[:δtime] = time() - startTime
   D[:rcNsampRange] = rcNsampRange
   !iszero(p.ρ_e) && (D[:rcNsampE]=rcNsampE; D[:rcNsampI]=rcNsampI)

   return D, F
end

#################################################
function randJ_generic(EorI, J̄, σJ, n, m; p = dp)
   if p.Jdist == "normal"
      J = EorI * J̄ .+ σJ * randn(n, m)
   elseif p.Jdist == "gamma"
      J = EorI * genGamma(J̄, σJ, n, m)
   elseif p.Jdist == "lognormal"
      J = EorI * genLognormal(J̄, σJ, n, m)
   elseif p.Jdist == "bernoulli"
      J = EorI * genBernoulli(J̄, n, m)
   end
   return J
end

#################################################
function randJ_I_FC(;p = dp)
   Random.seed!(p.seedJ)
   J̄ = p.J_ii * (p.K ^ p.scaleJ)
   σJ = p.g_ii * (p.K ^ p.scaleg)
   J = randJ_generic(-1, J̄, σJ, p.N, p.N, p = p)
   return J
end

#################################################
function randJ_I_PC(;p = dp)
   Random.seed!(p.seedJ)
   spMask = sprand(Bool, p.N, p.N, p.K/p.N)
   J_FC = randJ_I_FC(p = p)
   J = spMask .* J_FC
   return J
end

#################################################
function randJ_EI_FC(;p = dp)
   Random.seed!(p.seedJ)
   Ki = round(Int, p.K - p.K * p.ρ_e)
   Ke = p.K - Ki
   Ni = round(Int, p.N - p.N * p.ρ_e)
   Ne = p.N - Ni

   J̄ee = p.J_ee * (Ke ^ p.scaleJ)
   σJee = p.g_ee * (Ke ^ p.scaleg)
   Jee = randJ_generic(1,J̄ee,σJee,Ne,Ne,p=p)

   J̄ei = p.J_ei * (Ki ^ p.scaleJ)
   σJei = p.g_ei * (Ki ^ p.scaleg)
   Jei = randJ_generic(-1,J̄ei,σJei,Ne,Ni,p=p)

   J̄ie = p.J_ie * (Ke ^ p.scaleJ)
   σJie = p.g_ie * (Ke ^ p.scaleg)
   Jie = randJ_generic(1,J̄ie,σJie,Ni,Ne,p=p)

   J̄ii = p.J_ii * (Ki ^ p.scaleJ)
   σJii = p.g_ii * (Ki ^ p.scaleg)
   Jii = randJ_generic(-1,J̄ii,σJii,Ni,Ni,p=p)

   JfromE = vcat(Jee,Jie)
   JfromI = vcat(Jei,Jii)
   J = hcat(JfromE, JfromI)
   return J
end

#################################################
function randJ_EI_PC(;p = dp)
   Random.seed!(p.seedJ)
   spMask = sprand(Bool, p.N, p.N, p.K/p.N)
   J_FC = randJ_EI_FC(p = p)
   J = spMask .* J_FC
   return J
end

#################################################
function genBernoulli(μ, dim...)
   d = Bernoulli(μ)
   return rand(d, dim...)
end

#################################################
function genLognormal(μ, σ, dim...)
   M = log(μ^2 / sqrt(μ^2 + σ^2))
   Σ = sqrt(log((σ^2 / μ^2) + 1))
   d = LogNormal(M, Σ)
   return rand(d, dim...)
end

#################################################
function genGamma(μ, σ, dim...)
   θ = σ^2 / μ
   k = (μ / σ)^2
   d = Gamma(k, θ)
   return rand(d, dim...)
end

#################################################
uniform(a,b,n) = rand(n)*(b - a) .+ a
unroll(x) = reshape(x, length(x))
nanmean(x) = mean(filter(!isnan,x))
nanmean(x,dim) = dropdims(mapslices(nanmean,x,dims=dim),dims=dim)
nanstd(x) = std(filter(!isnan,x))
nanstd(x,dim) = dropdims(mapslices(nanstd,x,dims=dim),dims=dim)

#################################################
function myautocorr(x, maxLag; demean=true, normed=true)
   L = length(x)
   x = demean ? x .-= mean(x) : x
   yApn = zeros(2*L)
   yApn[1:L] .= x
   corr = irfft(abs.(rfft(yApn)).^2, length(yApn))
   B = floor(Int, length(corr)/2)
   corr = corr[1:B]
   lengths = B:-1:1
   corr = (corr ./ lengths)[1:maxLag]
   if var(x) == 0.0
      printstyled("Traj is flat.\n",color=:red)
      return ones(length(corr))
   else
      return (normed) ? corr / corr[1] : corr
   end
end

#################################################
function comp_ACF(y, maxLag)
   N = size(y, 2)
   nonzeroN = findall(!iszero, dropdims(var(y,dims=1),dims=1))
   ac = zeros(maxLag)
   for i in nonzeroN
      ac .+= myautocorr(y[:,i], maxLag, demean=false, normed=false)
   end
   ac /= length(nonzeroN)
   return ac
end

#################################################
function comp_β(τ, Rδη; p = dp)
   β = Rδη[1] ./ trapz(τ, Rδη)
   return β
end

#################################################
function printParams(J, I_i;p = dp)
   runparams = string(
            "N = ",       @sprintf("%d",  p.N),
         " | K = ",       @sprintf("%d",  p.K),
         " | J_ii = ",    @sprintf("%.2f",p.J_ii),
         " | g_ii = ",    @sprintf("%.2f",p.g_ii),
         " | ν = ",       @sprintf("%.2f",p.scaleg),
         " | I_i = ",     @sprintf("%.2f",p.I_i),
         " | I = ",       @sprintf("%.2f",I_i),
         " | mean(J) = ", @sprintf("%.4f",mean(J)),
         " | std(J) = ",  @sprintf("%.4f",std(J)),
         " | Φ = ",       @sprintf("%s",  p.Φname),
         " | λ = ",       @sprintf("%.2f",p.λ),
         " | dist = ",    @sprintf("%s",  p.Jdist),
         " | ρ_e = ",     @sprintf("%.2f",p.ρ_e),
         " | dt = ",      @sprintf("%.4f",p.dt),
         " | T = ",       @sprintf("%d",  p.T),
         " | s = ",       @sprintf("%d",  p.seedJ), "\n")
   printstyled(runparams, color=:green)
end

#################################################
function genfilename(p)
   filename = string("./data/",p.prefix,
                   @sprintf("_N%d",     p.N),
                   @sprintf("_K%d",     p.K),
                   @sprintf("_T%d",     p.T),
                   @sprintf("_d%.4f",   p.dt),
                   @sprintf("_rcDt%.4f",p.rcDt),
                   @sprintf("_tau%.4f", p.τ),
                   @sprintf("_Φ%s",     p.Φname),
                   @sprintf("_λ%.4f",   p.λ),
                   @sprintf("_Λ%s",     p.Jdist),
                   @sprintf("_II%.4f",  p.J_ii),
                   @sprintf("_scJ%.2f", p.scaleJ),
                   @sprintf("_gII%.4f", p.g_ii),
                   @sprintf("_scg%.4f", p.scaleg),
                   @sprintf("_bI%.4f",  p.I_i),
                   @sprintf("_scI%.2f", p.scaleI),
                   @sprintf("_sIC%d",   p.seedIC),
                   @sprintf("_sJ%d",    p.seedJ),
                   ".jld2")
   return filename
end

#################################################
function writefile(results, p)
   filename = genfilename(p)
   save(filename, "data", results)
end

#################################################
function readfile(p)
   filename = genfilename(p)
   return load(filename)["data"]
end

end

max_stars_repo_name: raminkhajeh/sparse-balance

max_issues_repo_name: raminkhajeh/sparse-balance

max_forks_repo_name: raminkhajeh/sparse-balance


text:
function ly_format_file(path::AbstractString; kwargs...)
    format_file(path; LY_JULIAFORMATTER..., kwargs...)
end
ly_format(paths...; kwargs...) = format(paths...; LY_JULIAFORMATTER..., kwargs...)
ly_format(pkgs::Module...; kwargs...) = ly_format(map(pkgdir, pkgs)...; kwargs...)

with_tempdir(f) = mktempdir(dir -> cd(() -> f(dir), dir))

# copied from Documenter.jl
function genkeys(; user = "\$USER", repo = "\$REPO", comment = "lyceumdevs@gmail.com")
    # Error checking. Do the required programs exist?
    success(`which which`) || error("'which' not found.")
    success(`which ssh-keygen`) || error("'ssh-keygen' not found.")

    directory = pwd()
    filename = tempname()

    isfile(filename) && error("temporary file '$(filename)' already exists in working directory")
    isfile("$(filename).pub") &&
    error("temporary file '$(filename).pub' already exists in working directory")

    # Generate the ssh key pair.
    success(`ssh-keygen -N "" -C $comment -f $filename`) ||
    error("failed to generate a SSH key pair.")

    # Prompt user to add public key to github then remove the public key.
    let url = "https://github.com/$user/$repo/settings/keys"
        @info("add the public key below to $url with read/write access:")
        println("\n", read("$filename.pub", String))
        rm("$filename.pub")
    end

    # Base64 encode the private key and prompt user to add it to travis. The key is
    # *not* encoded for the sake of security, but instead to make it easier to
    # copy/paste it over to travis without having to worry about whitespace.
    let travis_url = "https://travis-ci.com/$user/$repo/settings",
        github_url = "https://github.com/$user/$repo/settings/secrets"

        @info(
            "add a secure environment variable named to " *
            "$(travis_url) (if you deploy using Travis CI) or " *
            "$(github_url) (if you deploy using GitHub Actions) with value:"
        )
        println("\n", base64encode(read(filename, String)), "\n")
        rm(filename)
    end
end

hasmanifest(pkgdir::AbstractString) = manifestfile_path(pkgdir, strict = true) !== nothing
hasproject(pkgdir::AbstractString) = projectfile_path(pkgdir, strict = true) !== nothing

function with_sandbox_env(
    fn::Function,
    project::AbstractString = Base.active_project();
    copyall::Bool = false,
    extra_load_paths::Vector{<:AbstractString} = AbstractString[],
    default_load_path::Bool = false,
)
    if isdir(project)
        project = projectfile_path(project, strict = true)
        isnothing(project) && error("No project file found in $project")
    end

    old_dir = pwd()
    old_load_path = deepcopy(LOAD_PATH)
    old_project = Base.active_project()

    mktempdir() do sandbox
        if copyall
            cpinto(dirname(project), sandbox)
        else
            manifest = manifestfile_path(dirname(project), strict = true)
            manifest !== nothing && cp(manifest, joinpath(sandbox, basename(manifest)))
            cp(project, joinpath(sandbox, basename(project)))
        end

        try
            cd(sandbox)
            if default_load_path
                empty!(LOAD_PATH)
                Base.init_load_path()
            end
            Pkg.activate(pwd())
            append!(LOAD_PATH, map(normpath, extra_load_paths))
            return fn()
        finally
            cd(old_dir)
            append!(empty!(LOAD_PATH), old_load_path)
            Pkg.activate(old_project)
        end
    end
end

function cpinto(srcdir::AbstractString, dstdir::AbstractString)
    isdir(srcdir) || throw(ArgumentError("Expected a directory for `srcdir`, got $srcdir"))
    isdir(dstdir) || throw(ArgumentError("Expected a directory for `dstdir`, got $dstdir"))
    for file_or_dir in readdir(srcdir)
        cp(joinpath(srcdir, file_or_dir), joinpath(dstdir, file_or_dir))
    end
end

function flattendir(
    dir::AbstractString = pwd();
    join::Bool = false,
    sort::Bool = true,
    dirs::Bool = true,
    files::Bool = true,
)
    paths = String[]
    for (root, _dirs, _files) in walkdir(dir)
        if dirs
            for dir in _dirs
                push!(paths, joinpath(root, dir))
            end
        end
        if files
            for file in _files
                push!(paths, joinpath(root, file))
            end
        end
    end
    !join && (paths = map(path -> relpath(path, dir), paths))
    sort && sort!(paths)
    return paths
end

# check if child is a subpath of parent e.g.
# issubpath("/tmp/foo", "/tmp/foo/bar/baz.txt") == true
function issubpath(parent::AbstractString, child::AbstractString)
    startswith(realpath(child), realpath(parent))
end


text:
using Test
using GRASP
import GRASP.Configurations: CSFOrbital, CSFDefinition, CSFDefinitionList
import AtomicLevels: @ro_str, @c_str

@testset "csfs.jl" begin

@testset "specname" begin
    @test GRASP.specname(0) == "s"
    @test GRASP.specname(1) == "p"
    @test GRASP.specname(2) == "d"
    @test GRASP.specname(3) == "f"
end

@testset "CSFOrbital" begin
    @test CSFOrbital(3,0) == CSFOrbital(3,0)
    @test CSFOrbital(3,0) != CSFOrbital(3,1)
    @test CSFOrbital(3,0) != CSFOrbital(2,0)
    @test CSFOrbital(3,0) != CSFOrbital(2,1)

    # ordering
    orbitals = [CSFOrbital(n,l) for n=1:4 for l=0:(n-1)]
    for orb in orbitals
        @test orb == orb
    end

    for i=1:length(orbitals), j=(i+1):length(orbitals)
        a, b = orbitals[i], orbitals[j]
        @test a != b
        @test a < b
        @test a <= b
        @test b > a
        @test b >= a
    end

    @test CSFOrbital(1,0) < CSFOrbital(2,0) <  CSFOrbital(2,1)
    @test CSFOrbital(1,0) < CSFOrbital(2,0) <= CSFOrbital(2,0)

    @test maxelectrons(CSFOrbital(1,0)) == 2
    @test maxelectrons(CSFOrbital(2,0)) == 2
    @test maxelectrons(CSFOrbital(5,1)) == 6
    @test maxelectrons(CSFOrbital(5,2)) == 10
    @test maxelectrons(CSFOrbital(5,3)) == 14
    @test maxelectrons(CSFOrbital(5,4)) == 18
end

@testset "CSFDefinition" begin
    cdef1 = let cdef = CSFDefinition()
        for n=1:3, l=0:(n-1)
             push!(cdef, CSFOrbital(n, l), 2*(2l+1), 0)
        end
        cdef
    end

    cdef2 = let cdef = CSFDefinition()
        for n=4:5, l=0:min(3, n-1)
             push!(cdef, CSFOrbital(n, l), 2*(2l+1), 0)
        end
        cdef
    end

    @test string(cdef1) == "1s(2,i)2s(2,i)2p(6,i)3s(2,i)3p(6,i)3d(10,i)"
    @test nelectrons(cdef1) == 28

    cdef_combined = cdef1 * cdef2
    @test string(cdef_combined) == "1s(2,i)2s(2,i)2p(6,i)3s(2,i)3p(6,i)3d(10,i)4s(2,i)4p(6,i)4d(10,i)4f(14,i)5s(2,i)5p(6,i)5d(10,i)5f(14,i)"
    @test nelectrons(cdef_combined) == 28 + 64

    # Test `nexcitations(::CSFDefinition, ::CSFDefinition)`
    cdef3 = CSFDefinition()
    push!(cdef3, CSFOrbital(1, 0), 2, 0)
    push!(cdef3, CSFOrbital(2, 0), 2, 0)
    cdef4 = CSFDefinition()
    push!(cdef4, CSFOrbital(1, 0), 2, 0)
    push!(cdef4, CSFOrbital(2, 1), 2, 0)

    @test nexcitations(cdef3, cdef3) == 0
    @test nexcitations(cdef4, cdef4) == 0
    @test nexcitations(cdef3, cdef4) == 2

    # TODO: Test error for *?
end

import GRASP.Configurations: csfdefinition
import GRASP: CSF, Symmetries, AngularMomentum
@testset "csfdefinition()" begin
    csf = CSF([ro"1s", ro"3d-"], [2, 4], AngularMomentum[0, 0], AngularMomentum[0, 0], Symmetries.even)

    csfdef = csfdefinition(csf)
    @test nelectrons(csf) == nelectrons(csfdef)
end

@testset "parse(::CSFDefinition)" begin
    csfdef = parse(CSFDefinition, "2s(2,2) 3p(3)")
    @test isa(csfdef, CSFDefinition)
    @test length(csfdef.orbitals) == 2
    @test csfdef.orbitals[1] == CSFOrbital(2, 0)
    @test csfdef.orbitals[2] == CSFOrbital(3, 1)
    @test csfdef.nelectrons == [2, 3]
    @test csfdef.nexcitations == [2, 0]

    @test_throws ArgumentError parse(CSFDefinition, "2s(2)3p() 4p()")
    @test_throws ArgumentError parse(CSFDefinition, "2s(2)3p() 4pÖ()")
    @test_throws ArgumentError parse(CSFDefinition, "2s(2)3p() 4p a()")
    @test_throws ArgumentError parse(CSFDefinition, "2s 3pÖ()()")
    @test_throws ArgumentError parse(CSFDefinition, "2s 3pÖÖ(")
    @test_throws ArgumentError parse(CSFDefinition, "2s 3pÖÖ(())")
end

import GRASP.Configurations: parse_l
@testset "parse_l" begin
    @test parse_l("s")   == 0
    @test parse_l("s  ") == 0
    @test parse_l("p")   == 1
    @test parse_l("h")   == 5
end

import GRASP.Configurations: parse_j
@testset "parse_j" begin
    @test parse_j("s")  == -1
    @test parse_j("p-") ==  1
    @test parse_j("p")  == -2
    @test parse_j("d-") ==  2
    @test parse_j("d")  == -3
    @test parse_j("f-") ==  3
end

end # @testset "csfs.jl"


text:
export PowFunction

struct PowFunction{T, N} <: AbstractMathFun{T, N}
    pow::AbstractMathFun
    func::AbstractMathFun

    PowFunction{T, N}(func::AbstractMathFun, pow::T; force_emit=true) where {T, N} = new{T, N}(ConstFunction{T}(pow, force_emit = force_emit), func)
    PowFunction{T, N}(func::AbstractMathFun, pow::AbstractMathFun) where {T, N} = new{T, N}(pow, func)

    (f::PowFunction{T, N})(call::MFCall) where {T, N} = call(f.func(call) ^ f.pow(call))

    function Functions.string(f::PowFunction, emitter::FunctionEmitter)
        return "(" * string(f.func, emitter) * ")^" * string(f.pow, emitter)
    end

    function Functions.emit(f::PowFunction{T, N}, emitter::FunctionEmitter{T}) where {T, N}
        return Expr(:call, :(^), emit(f.func, emitter), emit(f.pow, emitter))
    end

    function Functions.simplify(f::PowFunction{T, N}, emit_constants::Bool) where {T, N}
        pows = simplify(f.pow, emit_constants)
        funcs = simplify(f.func, emit_constants)
        (is_constant(pows, emit_constants) && is_constant(funcs, emit_constants)) && (return ConstFunction{T}(functs.value ^ pows.value))
        (is_constant(pows, emit_constants) && pows.value == 0) && (return ConstFunction{T}(1))
        (is_constant(pows, emit_constants) && pows.value == 1) && (return funcs)
        return PowFunction{T, N}(funcs, pows)
    end
end


text:
"""
    module Faiss

An interface to the Faiss library for similarity-search of vectors (i.e. nearest-neighbour searching).

For basic usage, see [`Index`](@ref), [`add`](@ref), [`search`](@ref), [`add_with_ids`](@ref), [`local_rank`](@ref).
"""
module Faiss

using PythonCall
export np, Index, add, search, add_with_ids, remove_with_ids, add_search, local_rank, add_search_with_ids

const faiss = PythonCall.pynew()
const np = PythonCall.pynew()

function __init__()
    PythonCall.pycopy!(faiss, pyimport("faiss"))
    PythonCall.pycopy!(np, pyimport("numpy"))
end

"""
    Index(dim, spec="Flat")

Create a Faiss index of the given dimension.

The `spec` is an index factory string describing the type of index to construct.
"""
struct Index
    py::Py
end

# cpu index
Index(dim::Integer, str::AbstractString="Flat", metric::Integer=1) =
    Index(faiss.index_factory(convert(Int, dim), convert(String, str), metric))

function Index(dim::Integer; str::AbstractString="Flat", metric::String="L2", gpus::String="")
    # feat数据存储在这里面. 数据量巨大时,容易爆显存
    ENV["CUDA_VISIBLE_DEVICES"] = gpus
    if metric == "L2"
        metric_flag = faiss.METRIC_L2   # 1
    elseif metric == "IP"
        metric_flag = faiss.METRIC_INNER_PRODUCT  # 2
    end
    metric_flag = pyconvert(Integer, metric_flag)
    str_list = split(str, ",")
    if "IDMap2" in str_list
        str = join(str_list[str_list .!="IDMap2"] , ",")  # py faiss not support IDMap2.
        cpu_index = Index(dim, str, metric_flag)
        cpu_index = faiss.IndexIDMap2(cpu_index.py)
    else
        cpu_index = Index(dim, str, metric_flag)
        cpu_index = cpu_index.py
    end

    # println("faiss:", faiss.__version__, " gpus:", ENV["CUDA_VISIBLE_DEVICES"])
    if gpus == ""
        ngpus = 0
    else
        ngpus = length(split(ENV["CUDA_VISIBLE_DEVICES"], ","))
    end

    if ngpus == 0
        index = cpu_index
    elseif ngpus == 1
        res = faiss.StandardGpuResources()
        index = faiss.index_cpu_to_gpu(res, 0, cpu_index)  # use one gpu. make it into a gpu index
    else
        index = faiss.index_cpu_to_all_gpus(cpu_index)  # use all gpus
    end
    return Index(index)
end

"""
    size(idx::Index)

"""
Base.size(idx::Index) = (size(idx, 1), size(idx, 2))

Base.size(idx::Index, i::Integer) = i == 1 ? pyconvert(Int, idx.py.d) : i == 2 ? pyconvert(Int, idx.py.ntotal) : error()

"""
    show(io::IO, ::MIME"text/plain", idx::Index)

"""
function Base.show(io::IO, ::MIME"text/plain", idx::Index)
    metric_dict = Dict(1=>"METRIC_L2", 2=>"METRIC_INNER_PRODUCT")
    println(io, typeof(idx), " of ", size(idx, 2), " vectors of dimension ", size(idx, 1),
    ", metric_type:", metric_dict[pyconvert(Int64, idx.py.metric_type)])
end

"""
    add(idx::Index, vs::AbstractMatrix)

Add the columns of `vs` to the index.
"""
function add(idx::Index, vs::AbstractMatrix)
    size(vs, 2) == size(idx, 1) || error("expecting $(size(idx, 1)) rows")
    # vs_ = Py(convert(AbstractMatrix{Float32}, vs)').__array__()
    vs_ = convert(AbstractMatrix{Float32}, vs)
    vs_ = np.array(pyrowlist(vs_), dtype=np.float32)
    idx.py.add(vs_)
    return nothing
end

"""
    add_with_ids(idx::Index, vs::AbstractMatrix, ids::Array{Int64})

"""
function add_with_ids(idx::Index, vs::AbstractMatrix, ids::Array{Int64})
    size(vs, 2) == size(idx, 1) || error("expecting $(size(idx, 1)) rows")
    size(vs, 1) == size(ids, 1) || error("expecting $(size(vs, 1)) rows")
    # vs_ = Py(convert(AbstractMatrix{Float32}, vs)').__array__()
    vs_ = convert(AbstractMatrix{Float32}, vs)
    vs_ = np.array(pyrowlist(vs_), dtype=np.float32)
    ids_ = np.array(pyrowlist(ids), dtype=np.int64)
    idx.py.add_with_ids(vs_, ids_)
    return nothing
end

"""
    remove_with_ids(idx::Index, ids::Array{Int64})

"""
function remove_with_ids(idx::Index, ids::Array{Int64})
    ids_ = np.array(pyrowlist(ids), dtype=np.int64)
    idx.py.remove_ids(ids_)
end

"""
    search(idx::Index, vs::AbstractMatrix, k::Integer; threshold::Real=0.0)

Search the index for the `k` nearest neighbours of each column of `vs`.

Return `(D, I)` where `I` is a matrix where each column gives the ids of the `k` nearest
neighbours of the corresponding column of `vs` and `D` is the corresponding matrix of
distances.
"""
function search(idx::Index, vs::AbstractMatrix, k::Integer; threshold::Real=0.0)
    size(vs, 2) == size(idx, 1) || error("expecting $(size(idx, 1)) rows")
    # vs_ = Py(convert(AbstractMatrix{Float32}, vs)').__array__()
    vs_ = convert(AbstractMatrix{Float32}, vs)
    vs_ = np.array(pyrowlist(vs_), dtype=np.float32)
    k_ = convert(Int, k)
    th_ = convert(Float32, threshold)

    if pyconvert(Int64, idx.py.ntotal) == 0
        size_1 = size(vs, 1)
        D = zeros(Float32, (size_1, k))
        I = zeros(Int32, (size_1, k))
        if pyconvert(Int64, idx.py.metric_type) == 1
            D = D .+ 2
        end
    else
        D_, I_ = idx.py.search(vs_, k_, threshold=th_)
        D = pyconvert(Array{Float32, 2}, D_)
        I = pyconvert(Array{Int32, 2}, I_)
    end

    return (D, I)
end

"""
    add_search(idx::Index, vs_query::AbstractMatrix, vs_gallery::AbstractMatrix;
                k::Integer=100, threshold::Real=0.0, flag::Bool=true)

Add `vs_gallery` to idx and Search the index for the `k` nearest neighbours of each column of `vs_query`.

Return `(D, I)` where `I` is a matrix where each column gives the ids of the `k` nearest
neighbours of the corresponding column of `vs` and `D` is the corresponding matrix of distances.
"""
function add_search(idx::Index, vs_query::AbstractMatrix, vs_gallery::AbstractMatrix;
                    k::Integer=100, threshold::Real=0.0, flag::Bool=true)
    if flag
        add(idx, vs_gallery)
    end
    D, I = search(idx, vs_query, k; threshold=threshold)

    return (D, I)
end

"""
    add_search_with_ids(idx::Index, vs_query::AbstractMatrix, vs_gallery::AbstractMatrix, ids::Array{Int64};
                        k::Integer=100, threshold::Real=0.0, flag::Bool=true)

Add `vs_gallery` with `ids` to idx and Search the index for the `k` nearest neighbours of each column of `vs_query`.

"""
function add_search_with_ids(idx::Index, vs_query::AbstractMatrix, vs_gallery::AbstractMatrix, ids::Array{Int64};
    k::Integer=100, threshold::Real=0.0, flag::Bool=true)
    if flag
        add_with_ids(idx, vs_gallery, ids)
    end
    D, I = search(idx, vs_query, k; threshold=threshold)

    return (D, I)
end

"""
    local_rank(vs_query::AbstractMatrix, vs_gallery::AbstractMatrix; k::Integer=10,
                str::String="Flat", metric::String="L2", threshold::Real=0.0, gpus::String="")

Create Index and Add `vs_gallery` to idx and Search the index for the `k` nearest neighbours of each column of `vs_query`.

Return `(D, I)` where `I` is a matrix where each column gives the ids of the `k` nearest
neighbours of the corresponding column of `vs` and `D` is the corresponding matrix of distances.
"""
function local_rank(vs_query::AbstractMatrix, vs_gallery::AbstractMatrix; k::Integer=10,
                    str::String="Flat", metric::String="L2", threshold::Real=0.0, gpus::String="")
    feat_dim = size(vs_query, 2)
    idx = Index(feat_dim; str=str, metric=metric, gpus=gpus)
    D, I = add_search(idx, vs_query, vs_gallery; k=k, threshold=threshold)
    # PythonCall.pydel!(idx.py)
    return (D, I)
end

"""
    downcast(idx::Index)

Return the same index downcasted to its most specific type.
"""
downcast(idx::Index) = Index(faiss.downcast_index(idx.py))

end # module


text:
module JordanWigner

using ..AbstractOps
using ..AbstractPaulis: Iop, Xop, Yop, Zop
import ..Paulis, ..FermiTerm, ..FermiSum
using ..FermiOps
import ..PauliDefault
import ..OpSum, ..sort_and_sum_duplicates!

export jordan_wigner, jordan_wigner_fermi

## Jordan-Wigner helper
"""
    fill_pauli(pad, op_ind, fill_op::PauliT, end_op) where PauliT

Return a `Vector{PauliT}` of length `pad` with indices `1:op_ind-1` filled with `fill_op`
and index `op_ind` set to `end_op` and the remaining indices filled with `one(Pauli)`
"""
function fill_pauli(pad, op_ind, fill_op::PauliT, end_op) where PauliT
    str = Vector{PauliT}(undef, pad)
    if pad <  op_ind
        throw(DimensionMismatch("pad is less than op_ind."))
    end
    for i in 1:op_ind-1
       @inbounds str[i] = fill_op
    end
    str[op_ind] = end_op
    for i in op_ind+1:pad
        @inbounds str[i] = PauliT(Iop)
    end
    return str
end

# 1. Future compiler improvements may make this optimization worthless.
# 2. Keno's proposal for immutable Arrays would be very useful here.
const _c1 = [complex(-1/2), complex(-im/2)]
const _c2 = [complex(-1/2), complex(im/2)]
const _c3 = [complex(1/2), complex(-1/2)]
const _c4 = [complex(1/2), complex(1/2)]

## Be careful to copy these when combining. Otherwise _c1, etc. will be mutated
function jordan_wigner(op::FermiOp, op_ind::Integer, pad::Integer, ::Type{PauliT}=PauliDefault) where PauliT
    _fill_pauli(fill_op_ind, end_op_ind) = fill_pauli(pad, op_ind, PauliT(fill_op_ind), PauliT(end_op_ind))
    if op === Lower
        strx = _fill_pauli(Zop, Xop)
        stry = _fill_pauli(Zop, Yop)
        coeffs = _c1
    elseif op === Raise
        strx = _fill_pauli(Zop, Xop)
        stry = _fill_pauli(Zop, Yop)
        coeffs = _c2
    elseif op === NumberOp
        strx = fill(PauliT(Iop), pad)
        stry = _fill_pauli(Iop, Zop)
        coeffs = _c3
    elseif op === Empty
        strx = fill(PauliT(Iop), pad)
        stry = _fill_pauli(Iop, Zop)
        coeffs = _c4
    else
        throw(DomainError(op))
    end
    return OpSum{PauliT}([strx, stry], coeffs; already_sorted=true)
end

function jordan_wigner(term::FermiTerm, ::Type{PauliT}=PauliDefault) where PauliT
    pad = length(term)
    facs = []
    for (i, op) in enumerate(term)
        if op !== I #  op === Raise || op === Lower || op === NumberOp || op === Empty
            push!(facs, jordan_wigner(op, i, pad, PauliT))
        end
    end
    if isempty(facs)  # String is all I
        return(OpSum{PauliT}([fill(one(PauliT), length(term))], [complex(term.coeff)]))
    end
    psum = length(facs) == 1 ? only(facs) : reduce(*, facs) # TODO: better performance
    return term.coeff * psum # NOTE: This performs a copy, which is essential.
end

"""
    jordan_wigner(fsum::FermiSum, ::Type{PauliT}=PauliDefault) where PauliT

Using the Jordan-Wigner transform, convert the electronic Hamiltonian `fsum` to a qubit operator
of type `OpSum{PauliT}`.

    jordan_wigner(term::FermiTerm, ::Type{PauliT}=PauliDefault) where PauliT

Convert `term` to a qubit operator of type `OpSum{PauliT}`.
"""
function jordan_wigner(fsum::FermiSum, ::Type{PauliT}=PauliDefault) where PauliT
    psum = jordan_wigner(fsum[1], PauliT) # could use already sorted flag
    @inbounds for i in 2:length(fsum)
        append!(psum, jordan_wigner(fsum[i], PauliT))
    end
    return sort_and_sum_duplicates!(psum)
end

####
#### Jordan-Wigner using Fermi operators augmented by Z = N - E
####

## These are experimental, and so far do not seem very useful

function fill_fermi(pad, op_ind, fill_op, end_op)
    str = Vector{FermiOps.FermiOp}(undef, pad)
    if pad <  op_ind
        throw(DimensionMismatch("pad is less than op_ind."))
    end
    @inbounds for i in 1:op_ind-1
        str[i] = fill_op
    end
    str[op_ind] = end_op
    @inbounds for i in op_ind+1:pad
        str[i] = FermiOps.I
    end
    return str
end

## Jordan-Wigner
function jordan_wigner_fermi(op::FermiOp, op_ind::Integer, pad::Integer)
    if op === Lower
        str = fill_fermi(pad, op_ind, FermiOps.Z, FermiOps.Lower)
        return FermiTerm(str, complex(1.0))
    elseif op === Raise
        str = fill_fermi(pad, op_ind, FermiOps.Z, FermiOps.Raise)
        return FermiTerm(str, complex(1.0))
    elseif op === NumberOp
        str = fill_fermi(pad, op_ind, FermiOps.I, FermiOps.NumberOp)
        return FermiTerm(str, complex(1.0))
    elseif op === Empty
        str = fill_fermi(pad, op_ind, FermiOps.I, FermiOps.Empty)
        return FermiTerm(str, complex(1.0))
    elseif op === I
        str = fill(FermiOps.I, pad)
        return FermiTerm(str, complex(1.0))
    else
        throw(DomainError(op))
    end
end

function jordan_wigner_fermi(term::FermiTerm)
    pad = length(term)
    facs = [] # FermiTerm{FermiOp, Vector{FermiOp}, ComplexF64}[]
    for (i, op) in enumerate(term)
        if op === Raise || op === Lower || op === NumberOp
            push!(facs, jordan_wigner_fermi(op, i, pad))
        end
    end
    if isempty(facs)
        return (1.0 + 0.0im) * term
    end
    return term.coeff * reduce(*, facs)  # TODO: performance
end

function jordan_wigner_fermi(fsum::FermiSum)
    terms = FermiTerm{FermiOp, Vector{FermiOp}, ComplexF64}[]
    for i in eachindex(fsum) # 1:length(fsum)
        push!(terms, jordan_wigner_fermi(fsum[i]))
    end
    nterms = [x for x in terms]  # nterms is never used
    return FermiSum(terms)
end

end # module JordanWigner


text:
@testset "Simple Reaction System" begin
    rn = @reaction_network ABC begin
        1, A + B --> C
    end

    bg_rn = BondGraph(rn)

    @test bg_rn.name == :ABC
    @test nv(bg_rn) == 5
    @test ne(bg_rn) == 4

    @test any(n -> n.name == :R1, bg_rn.nodes)
    @test any(n -> n isa EqualFlow && name(n) == Symbol("1"), bg_rn.nodes)

    @test length(getnodes(bg_rn, :Ce)) == 3
    @test length(getnodes(bg_rn, EqualFlow)) == 1
    @test length(getnodes(bg_rn, :Re)) == 1

    @test Graphs.degree(bg_rn) == [2, 3, 1, 1, 1]
end

@testset "Reversible MM" begin
    rn = @reaction_network MM_reversible begin
        (1, 1), E + S <--> C
        (1, 1), C <--> E + P
    end

    bg_rn = BondGraph(rn; chemostats=["S", "P"])

    @test bg_rn.name == :MM_reversible
    @test nv(bg_rn) == 10
    @test ne(bg_rn) == 10

    @test length(getnodes(bg_rn, :Ce)) == 2
    @test length(getnodes(bg_rn, :Se)) == 2
    @test length(getnodes(bg_rn, EqualEffort)) == 2
    @test length(getnodes(bg_rn, EqualFlow)) == 2
    @test length(getnodes(bg_rn, :Re)) == 2

    @test Graphs.degree(bg_rn) == [2, 3, 1, 1, 1, 2, 3, 1, 3, 3]
end

@testset "Stoichiometry Test" begin
    rn = @reaction_network Stoichiometry begin
        1, 3A + 2B --> 5C
    end

    bg_rn = BondGraph(rn)

    @test nv(bg_rn) == 8
    @test ne(bg_rn) == 7

    tfs = getnodes(bg_rn, :TF)
    @test length(tfs) == 3
    @test repr.(tfs) == ["TF:3", "TF:2", "TF:5"]
end

@testset "SERCA" begin
    rn = @reaction_network SERCA begin
        (1, 1), P1 + MgATP <--> P2
        (1, 1), P2 + H <--> P2a
        (1, 1), P2 + 2Cai <--> P4
        (1, 1), P4 <--> P5 + 2H
        (1, 1), P5 <--> P6 + MgADP
        (1, 1), P6 <--> P8 + 2Casr
        (1, 1), P8 + 2H <--> P9
        (1, 1), P9 <--> P10 + H
        (1, 1), P10 <--> P1 + Pi
    end

    chemostats = ["MgATP", "MgADP", "Pi", "H", "Cai", "Casr"]
    bg_rn = BondGraph(rn; chemostats)

    @test nv(bg_rn) == 46
    @test ne(bg_rn) == 49
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
_codegen(a::AbstractArray) = a

function _codegen(::Type{T}, f::Function, args::Array{Sym}) where {T}
    ts = _codegen.(args)
    @q begin
        $f($(ts...))
    end
end

Base.getindex(a::Sym{A}, inds...) where {T, A <: AbstractArray{T}} = term(getindex, a, inds...; type=T)

SymbolicUtils.promote_symtype(::typeof(getindex), ::Type{A}, ::Type{<:Integer}...) where {T, A <: AbstractArray{T}} = T


text:
# Normalized is basically wrapper like Symmetric

"""
   Normalized(dist::Union{StringSemiMetric, StringMetric})

Creates a normalized distance. The normalized distance always return a Float64 between 0.0 and 1.0 (or a missing if one of the argument is missing).
A Normalized Distance has a keyword argument `max_dist` that defaults to 1.0. It returns 1.0 if the true distance is higher than `max_dist`.

### Examples
```julia-repl
julia> s1 = "New York Mets vs Atlanta"
julia> s2 = "Atlanta Braves vs New York Mets"
julia> Levenshtein()(s1, s2)
25
julia> StringDistances.Normalized(Levenshtein())(s1, s2)
0.8064
```
"""
struct Normalized{T <: Union{StringSemiMetric, StringMetric}} <: StringSemiMetric
    dist::T
end
Normalized(dist::Normalized) = dist

# Consider all distances to be normalized by default
function (dist::Normalized)(s1, s2; max_dist = 1.0)
    out = dist.dist(s1, s2; max_dist = max_dist)
    max_dist !== nothing && out > max_dist && return 1.0
    return out
end

function (dist::Normalized{<:Union{Hamming, DamerauLevenshtein}})(s1, s2; max_dist = 1.0)
    (s1 === missing) | (s2 === missing) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len2 == 0 && return 0.0
    out = dist.dist(s1, s2) / len2
    max_dist !== nothing && out > max_dist && return 1.0
    return out
end

function (dist::Normalized{<:Union{Levenshtein, OptimalStringAlignement}})(s1, s2; max_dist = 1.0)
    (s1 === missing) | (s2 === missing) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len2 == 0 && return 0.0
    if max_dist == 1.0
        d = dist.dist(s1, s2)
    else
        d = dist.dist(s1, s2; max_dist = ceil(Int, len2 * max_dist))
    end
    out = d / len2
    max_dist !== nothing && out > max_dist && return 1.0
    return out
end

function (dist::Normalized{<:AbstractQGramDistance})(s1, s2; max_dist = 1.0)
    (s1 === missing) | (s2 === missing) && return missing
    # When string length < q for qgram distance, returns s1 == s2
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 <= dist.dist.q - 1 && return Float64(s1 != s2)
    if dist.dist isa QGram
        out = dist.dist(s1, s2) / (len1 + len2 - 2 * dist.dist.q + 2)
    else
        out = dist.dist(s1, s2)
    end
    max_dist !== nothing && out > max_dist && return 1.0
    return out
end


text:
# Use
#
#     DOCUMENTER_DEBUG=true julia --color=yes make.jl local [fixdoctests]
#
# for local builds.

using Documenter
using Pkg
ENV["PYTHON"]=""
Pkg.build("PyCall")
using Plots
pyplot(fmt=:svg)
using RadiationSpectra

makedocs(
    sitename = "RadiationSpectra.jl",
    modules = [RadiationSpectra],
    pages=[
        "Home" => "index.md",
        "Manual" => Any[
            "Peak Finding" => "man/spectrum_deconvolution.md",
            "Fitting" => "man/fitting.md",
            "Calibration" => "man/spectrum_calibration.md",
        ],
        "API" => "api.md",
        "LICENSE" => "LICENSE.md",
    ],
    doctest = ("fixdoctests" in ARGS) ? :fix : true,
    format = Documenter.HTML(   canonical = "https://JuliaHEP.github.io/RadiationSpectra.jl/stable/",
                                prettyurls = !("local" in ARGS)
    )
)

deploydocs(
    repo = "github.com/JuliaHEP/RadiationSpectra.jl.git",
)


text:
#
# Layout combinators, tiling, and rendering.
#

# Layout tree.

abstract type AbstractBlock end

struct Layout
    blk::AbstractBlock
    args::Vector{Layout}
end

struct LiteralBlock <: AbstractBlock
    atom::Any
    len::Int
end

struct HorizontalBlock <: AbstractBlock
end

struct VerticalBlock <: AbstractBlock
end

struct ChoiceBlock <: AbstractBlock
end

struct PenaltyBlock <: AbstractBlock
    cost::Int
    break_factor::Int
    spill_factor::Int

    PenaltyBlock(; cost::Int=0, break_factor::Int=1, spill_factor::Int=1) =
        new(cost, break_factor, spill_factor)
end

let NO_ARGS = Layout[],
    INDENTS = Layout[Layout(LiteralBlock(" " ^ i, i), NO_ARGS) for i = 1:DEFAULT_LINE_WIDTH]

global literal, indent

literal(str::String) =
    Layout(LiteralBlock(str, textwidth(str)), NO_ARGS)

literal(str::Union{AbstractString,Symbol}) =
    literal(string(str))

literal(atom, len::Int) =
    Layout(LiteralBlock(atom, len), NO_ARGS)

indent(i::Int) =
    1 <= i <= length(INDENTS) ? INDENTS[i] : Layout(LiteralBlock(" " ^ i, i), NO_ARGS)

end

const ZERO = literal("")

# Visualization.

show(io::IO, lt::Layout) =
    pprint(io, tile(lt))

# Layout combinators.

function (*)(lts::Vector{Layout})
    lts = filter(!isequal(ZERO), lts)
    if isempty(lts)
        ZERO
    elseif length(lts) == 1
        lts[1]
    else
        Layout(HorizontalBlock(), lts)
    end
end

(/)(lts::Vector{Layout}) =
    Layout(VerticalBlock(), lts)

(|)(lts::Vector{Layout}) =
    Layout(ChoiceBlock(), lts)

(*)(lt1::Layout, lts::Layout...) =
    (*)(Layout[lt1, lts...])

(/)(lt1::Layout, lts::Layout...) =
    (/)(Layout[lt1, lts...])

(|)(lt1::Layout, lts::Layout...) =
    (|)(Layout[lt1, lts...])

penalize(lt::Layout; cost::Int=0, break_factor::Int=1, spill_factor::Int=1) =
    Layout(PenaltyBlock(cost=cost, break_factor=break_factor, spill_factor=spill_factor), Layout[lt])

nobreak(lt::Layout) =
    let lt′ = nobreak(lt.blk, lt.args)
        lt′ isa Layout && lt′.blk == lt.blk && lt′.args == lt.args ? lt : lt′
    end

nobreak(blk::LiteralBlock, args::Vector{Layout}) =
    Layout(blk, args)

function nobreak(blk::HorizontalBlock, args::Vector{Layout})
    args′ = Layout[]
    for arg in args
        arg′ = nobreak(arg)
        arg′ !== nothing || return nothing
        push!(args′, arg′)
    end
    Layout(blk, args′)
end

nobreak(::VerticalBlock, ::Vector{Layout}) = nothing

function nobreak(blk::Union{ChoiceBlock,PenaltyBlock}, args::Vector{Layout})
    args′ = Layout[]
    for arg in args
        arg′ = nobreak(arg)
        arg′ !== nothing || continue
        push!(args′, arg′)
    end
    !isempty(args′) || return nothing
    Layout(blk, args′)
end

function list_layout(items::Vector{Layout};
                     prefix::Union{String,Symbol,Layout}="",
                     par::Tuple{String,String}=("(", ")"),
                     sep::String=", ",
                     sep_brk=:end,  # :start, :end, :both, or :none
                     tab::Int=4,
                     nobrk::Int=10)
    prefix_lt = prefix isa Layout ? prefix : literal(prefix)
    head_lt = prefix_lt * literal(par[1])
    tail_lt = literal(par[2])
    !isempty(items) || return head_lt * tail_lt
    sepc_lt = literal(sep)
    sepl_lt = literal(sep_brk == :start || sep_brk == :both ? lstrip(sep) : "")
    sepr_lt = literal(sep_brk == :end || sep_brk == :both ? rstrip(sep) : "")
    tab_lt = indent(tab)
    vlt = items[1]
    #hlt = penalize(items[1], break_factor=nobrk)
    hlt = items[1]
    for item in items[2:end]
        vlt = (vlt * sepr_lt) / (sepl_lt * item)
        hlt = hlt * sepc_lt * item
    end
    vlt = (head_lt | (head_lt / tab_lt)) * vlt * tail_lt
    #hlt = penalize(head_lt * hlt * tail_lt, break_factor=nobrk)
    #vlt | hlt
    hlt = nobreak(head_lt * hlt * tail_lt)
    (hlt !== nothing ? vlt | hlt : vlt)::Layout
end

function pair_layout(fst::Layout,
                     snd::Layout;
                     sep::String=" => ",
                     sep_brk=:end,  # :start, :end, :both, or :none
                     tab::Int=4)
    nbfst = nobreak(fst)
    sepc_lt = literal(sep)
    sepl_lt = literal(sep_brk == :start || sep_brk == :both ? lstrip(sep) : "")
    sepr_lt = literal(sep_brk == :end || sep_brk == :both ? rstrip(sep) : "")
    tab_lt = indent(tab)
    nbfst !== nothing ?
        ((nobreak(fst) * sepc_lt) | ((fst * sepr_lt) / (tab_lt * sepl_lt))) * snd :
        ((fst * sepr_lt) / (tab_lt * sepl_lt)) * snd
end

# Fallback layout.

tile_repr(obj) =
    literal(repr(obj))

tile(obj) = tile_expr_or_repr(quoteof(obj))

# Layouts for standard types.

function tile(s::String)
    lt = literal(repr(s))
    if '"' in s
        if endswith(s, '"')
            es = escape_string(s[1:end - 1], "") * "\\\""
        else
            es = escape_string(s, "")
        end
        lt = lt | literal("\"\"\"" * replace(es, "\"\"\"" => "\"\"\\\"") * "\"\"\"")
    end
    if '\n' in s
        lines = split(s, '\n')
        n = length(lines)
        for k = 1:n
            line = lines[k]
            if k == n && endswith(line, '"')
                line = escape_string(line[1:end - 1], "") * "\\\""
            else
                line = escape_string(line, "")
            end
            if '"' in line
                line = replace(line, "\"\"\"" => "\"\"\\\"")
            end
            lines[k] = line
        end
        lines[end] *= "\"\"\""
        indented = true
        for line in lines
            if !isempty(line) && !(first(line) in (' ', '\t'))
                indented = false
                break
            end
        end
        if indented
            line = lines[end]
            lines[end] = "\\x" * string(UInt32(first(line)), base = 16, pad = 2) * line[2:end]
        end
        mlt = literal("\"\"\"")
        for line in lines
            mlt = mlt / literal(line)
        end
        lt = lt | mlt
    end
    lt
end

tile(p::Pair) =
    pair_layout(tile(p.first), tile(p.second))

tile(t::Tuple) =
    if length(t) == 1
        literal("(") * tile(t[1]) * literal(",)")
    else
        list_layout(Layout[tile(x) for x in t])
    end

tile(v::Vector) =
    list_layout(Layout[tile(x) for x in v], par=("[", "]"))

tile(s::Set) =
    if isempty(s)
        literal("Set()")
    else
        list_layout(Layout[list_layout([tile(x) for x in s], par=("[", "]"))],
                    prefix=:Set)
    end

tile(t::NamedTuple) =
    if length(t) == 1
        ((key, val),) = pairs(t)
        literal("(") * pair_layout(literal(key), tile(val), sep=" = ") * literal(",)")
    else
        list_layout(Layout[pair_layout(literal(key), tile(val), sep=" = ")
                           for (key, val) in pairs(t)])
    end

tile(d::Dict) =
    list_layout(Layout[pair_layout(tile(key), tile(val)) for (key, val) in d],
                prefix=:Dict)

# Layout of the layout tree.

tile(lt::Layout) =
    tile_layout(lt)

tile_layout(lt::Layout, precedence::Int=0) =
    tile_layout(lt.blk, lt.args, precedence)

function tile_layout(blk::LiteralBlock, ::Vector{Layout}, ::Int)
    if blk.atom isa String && blk.len == textwidth(blk.atom)
        if blk.len > 0 && all(ch == ' ' for ch in blk.atom)
            literal("indent($(blk.len))")
        else
            literal("literal($(repr(blk.atom)))")
        end
    else
        literal("literal($(repr(blk.atom)), $(blk.len))")
    end
end

function tile_layout(blk::PenaltyBlock, args::Vector{Layout}, ::Int)
    items = Layout[tile_layout(arg) for arg in args]
    if blk.cost != 0
        push!(items, literal("cost=$(blk.cost)"))
    end
    if blk.break_factor != 1
        push!(items, literal("break_factor=$(blk.break_factor)"))
    end
    if blk.spill_factor != 1
        push!(items, literal("spill_factor=$(blk.spill_factor)"))
    end
    list_layout(items, prefix=:penalize)
end

function tile_layout(blk::Union{HorizontalBlock, VerticalBlock, ChoiceBlock}, args::Vector{Layout}, precedence::Int)
    precedence′ = blk isa ChoiceBlock ? 0 : 1
    items = Layout[]
    first = true
    for arg in args
        item = tile_layout(arg, first ? precedence′ : precedence′ + 1)
        push!(items, item)
        first = false
    end
    par = precedence > precedence′ ? ("(", ")") : ("", "")
    sep = blk isa HorizontalBlock ? " * " :
          blk isa VerticalBlock ? " / " :
          blk isa ChoiceBlock ? " | " : ""
    return list_layout(items, par=par, sep=sep)
end

# Rendering

function render(io::IO, lt::Layout, col::Int=0, nl::Vector{UInt8}=UInt8[0x0A, 0x20, 0x20, 0x20, 0x20])
    blk = lt.blk
    args = lt.args
    if blk isa LiteralBlock
        col = render(io, blk, args, col, nl)
    elseif blk isa HorizontalBlock
        col = render(io, blk, args, col, nl)
    elseif blk isa VerticalBlock
        col = render(io, blk, args, col, nl)
    else
        col = render(io, blk, args, col, nl)::Int
    end
    col
end

function render(io::IO, blk::LiteralBlock, ::Vector{Layout}, col::Int, ::Vector{UInt8})
    atom = blk.atom
    if atom isa String
        print(io, atom)
    else
        print(io, atom)
    end
    col + blk.len
end

function render(io::IO, ::HorizontalBlock, args::Vector{Layout}, col::Int, nl::Vector{UInt8})
    for arg in args
        col = render(io, arg, col, nl)
    end
    col
end

function render(io::IO, ::VerticalBlock, args::Vector{Layout}, col::Int, nl::Vector{UInt8})
    l = length(nl)
    if col + 1 > l
        resize!(nl, col + 1)
        for k = l:col
            @inbounds nl[k + 1] = 0x20
        end
    end
    col′ = col
    first = true
    for arg in args
        if !first
            GC.@preserve nl unsafe_write(io, pointer(nl), col + 1)
        end
        col′ = render(io, arg, col, nl)
        first = false
    end
    col′
end

render(io::IO, ::Union{ChoiceBlock, PenaltyBlock}, args::Vector{Layout}, col::Int, nl::Vector{UInt8}) =
    !isempty(args) ? render(io, args[1], col, nl) : col