text:
# Tests for Radial DEA Models
@testset "DirectionalDEAModel" begin

    ## Test Directional DF DEA Models with FLS Book data
    # Test against results in R
    X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17]
    Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12]

    # Observed CRS
    deaddfobs = deaddf(X, Y, Gx = X, Gy = Y, rts = :CRS)

    @test typeof(deaddfobs) == DirectionalDEAModel

    @test nobs(deaddfobs) == 11
    @test ninputs(deaddfobs) == 2
    @test noutputs(deaddfobs) == 1
    @test efficiency(deaddfobs) ≈ [0.00000000000;
                               0.23282544774;
                               0.09898790422;
                               0.00000000000;
                               0.52628538417;
                               0.28571428571;
                               0.00000000000;
                               0.13787061050;
                               0.09883720930;
                               0.34177215190;
                               0.00000000000]
    @test slacks(deaddfobs, :X) ≈ [0.000000000   0;
                               0.000000000   0;
                               0.000000000   0;
                               0.000000000   0;
                               0.000000000   0;
                               5.714285714   0;
                               0.000000000   0;
                               0.000000000   0;
                               1.802325581   0;
                               0.000000000   0;
                               0.000000000   4]
    @test slacks(deaddfobs, :Y) ≈ zeros(11)

    @test efficiency(deaddf(X, Y, Gx = :Observed, Gy = :Observed, rts = :CRS)) == efficiency(deaddfobs)

    @test efficiency(deaddf(targets(deaddfobs, :X), targets(deaddfobs, :Y), Gx = :Observed, Gy = :Observed, rts = :CRS, slack = false)) ≈ zeros(11) atol=1e-15
    @test efficiency(deaadd(targets(deaddfobs, :X), targets(deaddfobs, :Y))) ≈ zeros(11) atol=1e-13

    # Observed VRS
    deaddfobsvrs = deaddf(X, Y, Gx = X, Gy = Y, rts = :VRS)

    @test nobs(deaddfobsvrs) == 11
    @test ninputs(deaddfobsvrs) == 2
    @test noutputs(deaddfobsvrs) == 1
    @test efficiency(deaddfobsvrs) ≈ [0.00000000000;
                              0.1076130509;
                              0.0000000000;
                              0.0000000000;
                              0.2883597884;
                              0.0000000000;
                              0.00000000000;
                              0.0000000000;
                              0.0000000000;
                              0.1821192053;
                              0.00000000000]
    @test slacks(deaddfobsvrs, :X) ≈ [0.000000000   0;
                              0.000000000   0;
                              0.000000000   0;
                              0.000000000   0;
                              0.000000000   0;
                              0   0;
                              0.000000000   0;
                              0.000000000   0;
                              0   0;
                              0.000000000   4.32781457;
                              0.000000000   4]
    @test slacks(deaddfobsvrs, :Y) ≈ [0.000000000;
                             0.000000000;
                             0.000000000;
                             0.000000000;
                             0.3915343915;
                             0.000000000;
                             0.000000000;
                             0.000000000;
                             0.000000000;
                             0.000000000;
                             0.000000000]

    @test efficiency(deaddf(X, Y, Gx = :Observed, Gy = :Observed, rts = :VRS)) == efficiency(deaddfobsvrs)

    @test efficiency(deaddf(targets(deaddfobsvrs, :X), targets(deaddfobsvrs, :Y), Gx = :Observed, Gy = :Observed, rts = :VRS, slack = false)) ≈ zeros(11) atol=1e-15
    @test efficiency(deaadd(targets(deaddfobsvrs, :X), targets(deaddfobsvrs, :Y))) ≈ zeros(11) atol=1e-12

    # Ones CRS
    deaddfones = deaddf(X, Y, Gx = ones(size(X)), Gy = ones(size(Y)), rts = :CRS)

    @test nobs(deaddfones) == 11
    @test ninputs(deaddfones) == 2
    @test noutputs(deaddfones) == 1
    @test efficiency(deaddfones) ≈ [0.00000000000;
                             3.219963031;
                             2.121693122;
                             0.00000000000;
                             6.735674677;
                             1.945945946;
                             0.00000000000;
                             3.635859519;
                             1.837837838;
                             10.231053604;
                             0.00000000000]
    @test slacks(deaddfones, :X) ≈ [0.000000000   0;
                             0.000000000   0;
                             0.000000000   0;
                             0.000000000   0;
                             0.000000000   0;
                             10.918918919   0;
                             0.000000000   0;
                             0.000000000   0;
                             4.756756757    0;
                             0.000000000   0;
                             0.000000000   4]
    @test slacks(deaddfones, :Y) ≈ zeros(11)

    @test efficiency(deaddf(X, Y, Gx = :Ones, Gy = :Ones, rts = :CRS)) == efficiency(deaddfones)

    # Ones VRS
    deaddfonesvrs = deaddf(X, Y, Gx = ones(size(X)), Gy = ones(size(Y)), rts = :VRS)

    @test nobs(deaddfonesvrs) == 11
    @test ninputs(deaddfonesvrs) == 2
    @test noutputs(deaddfonesvrs) == 1
    @test efficiency(deaddfonesvrs) ≈ [0.000000000;
                            1.418867925;
                            0.000000000;
                            0.000000000;
                            4.067924528;
                            0.000000000;
                            0.000000000;
                            0.000000000;
                            0.000000000;
                            5.000000000;
                            0.000000000]
    @test slacks(deaddfonesvrs, :X) ≈ [0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0   0;
                            0.000000000   0;
                            0.000000000   0;
                            0   0;
                            0.000000000   6;
                            0.000000000   4]
    @test slacks(deaddfonesvrs, :Y) ≈ zeros(11) atol = 1e-14

    @test efficiency(deaddf(X, Y, Gx = :Ones, Gy = :Ones, rts = :VRS)) == efficiency(deaddfonesvrs)

    # Only X CRS
    deaddfonlyX = deaddf(X, Y, Gx =  X, Gy = zeros(size(Y)), rts = :CRS)

    @test nobs(deaddfonlyX) == 11
    @test ninputs(deaddfonlyX) == 2
    @test noutputs(deaddfonlyX) == 1
    @test efficiency(deaddfonlyX) ≈ [0.0000000000;
                            0.3777103209;
                            0.1801437556;
                            0.0000000000;
                            0.6896290689;
                            0.4444444444;
                            0.0000000000;
                            0.2423309104;
                            0.1798941799;
                            0.5094339623;
                            0.0000000000]
    @test slacks(deaddfonlyX, :X) ≈ [0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            4.444444444   0;
                            0.000000000   0;
                            0.000000000   0;
                            1.640211640   0;
                            0.000000000   0;
                            0.000000000   4]
    @test slacks(deaddfonlyX, :Y) ≈ zeros(11)

    @test efficiency(deaddf(X, Y, Gx = :Observed, Gy = :Zeros, rts = :CRS)) == efficiency(deaddfonlyX)

    # Only X VRS
    deaddfonlyXvrs = deaddf(X, Y, Gx = X, Gy = zeros(size(Y)), rts = :VRS)

    @test nobs(deaddfonlyXvrs) == 11
    @test ninputs(deaddfonlyXvrs) == 2
    @test noutputs(deaddfonlyXvrs) == 1
    @test efficiency(deaddfonlyXvrs) ≈ [0.0000000000;
                            0.1300138313;
                            0.0000000000;
                            0.0000000000;
                            0.2883597884;
                            0.0000000000;
                            0.0000000000;
                            0.0000000000;
                            0.0000000000;
                            0.5068790731;
                            0.0000000000]
    @test slacks(deaddfonlyXvrs, :X) ≈ [0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0   0;
                            0.000000000   0;
                            0.000000000   0;
                            0   0;
                            0.000000000   0;
                            0.000000000   4]
    @test slacks(deaddfonlyXvrs, :Y) ≈ [0.000000000;
                            0.000000000;
                            0.000000000;
                            0.000000000;
                            2.698412698;
                            0.000000000;
                            0.000000000;
                            0.000000000;
                            0.000000000;
                            0.000000000;
                            0.000000000]

    @test efficiency(deaddf(X, Y, Gx = :Observed, Gy = :Zeros, rts = :VRS)) == efficiency(deaddfonlyXvrs)

    # Only Y CRS
    deaddfonlyY = deaddf(X, Y, Gx = zeros(size(X)), Gy = Y, rts = :CRS)

    @test nobs(deaddfonlyY) == 11
    @test ninputs(deaddfonlyY) == 2
    @test noutputs(deaddfonlyY) == 1
    @test efficiency(deaddfonlyY) ≈ [0.0000000000;
                            0.6069686411;
                            0.2197260274;
                            0.0000000000;
                            2.2219512195;
                            0.8000000000;
                            0.0000000000;
                            0.3198373984;
                            0.2193548387;
                            1.0384615385;
                            0.0000000000]
    @test slacks(deaddfonlyY, :X) ≈ [0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            8   0;
                            0.000000000   0;
                            0.000000000   0;
                            2   0;
                            0.000000000   0;
                            0.000000000   4]
    @test slacks(deaddfonlyY, :Y) ≈ zeros(11)

    @test efficiency(deaddf(X, Y, Gx = :Zeros, Gy = :Observed, rts = :CRS)) == efficiency(deaddfonlyY)

    # Only Y VRS
    deaddfonlyYvrs = deaddf(X, Y, Gx = zeros(size(X)), Gy = Y, rts = :VRS)

    @test nobs(deaddfonlyYvrs) == 11
    @test ninputs(deaddfonlyYvrs) == 2
    @test noutputs(deaddfonlyYvrs) == 1
    @test efficiency(deaddfonlyYvrs) ≈ [0.0000000000;
                            0.5075187970;
                            0.0000000000;
                            0.0000000000;
                            2.2039473684;
                            0.0000000000;
                            0.0000000000;
                            0.0000000000;
                            0.0000000000;
                            0.1923076923;
                            0.0000000000]
    @test slacks(deaddfonlyYvrs, :X) ≈ [0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0.000000000   0;
                            0   0;
                            0.000000000   0;
                            0.000000000   0;
                            0   0;
                            5   11;
                            0.000000000   4]
    @test slacks(deaddfonlyYvrs, :Y) ≈ zeros(11) atol = 1e-14

    @test efficiency(deaddf(X, Y, Gx = :Zeros, Gy = :Observed, rts = :VRS)) == efficiency(deaddfonlyYvrs)

    # Mean CRS
    deaddfmean = deaddf(X, Y, Gx = :Mean, Gy = :Mean, rts = :CRS)

    @test efficiency(deaddfmean) ≈ [0; 0.1713509018; 0.1082533683; 0; 0.3584401184; 0.1148158887; 0; 0.1934829069; 0.1084372282; 0.5444473258; 0]

    # Mean VRS
    deaddfmeanvrs = deaddf(X, Y, Gx = :Mean, Gy = :Mean, rts = :VRS)

    @test efficiency(deaddfmeanvrs) ≈ [0; 0.08056567388; 0; 0;0.23098350118; 0; 0; 0; 0; 0.24886877828; 0]

    # Test no slacks
    deaddfnoslack = deaddf(X, Y, Gx = X, Gy = Y, slack = false)
    @test efficiency(deaddfnoslack) == efficiency(deaddfobs)
    @test isempty(slacks(deaddfnoslack, :X)) == 1
    @test isempty(slacks(deaddfnoslack, :Y)) == 1

    @test efficiency(dea(targets(deaddfnoslack, :X), targets(deaddfnoslack, :Y), slack = false)) ≈ ones(11)
    @test efficiency(deaadd(targets(deaddfnoslack, :X), targets(deaddfnoslack, :Y))) != zeros(11) # Different as there is no slacks in first model

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deaddfobs_ref_eff = zeros(size(X, 1))

    deaddfobsvs_ref_eff = zeros(size(X, 1))

    deaddfvrs_ref_slackX = zeros(size(X))
    deaddfvrs_ref_slackY = zeros(size(Y))

    Gx = X[:,:]
    Gy = Y[:,:]

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]
        Gxeval = Gx[i:i,:]
        Gxeval = Gxeval[:,:]
        Gyeval = Gy[i:i,:]
        Gyeval = Gyeval[:,:]

        deaddfobs_ref_eff[i] = efficiency(deaddf(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :CRS, Xref = Xref, Yref = Yref))[1]

        deaddfobsvs_ref_eff[i] = efficiency(deaddf(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deaddfvrs_ref_slackX[i,:] = slacks(deaddf(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :VRS, Xref = Xref, Yref = Yref), :X)
        deaddfvrs_ref_slackY[i,:] = slacks(deaddf(Xeval, Yeval, Gx = Gxeval, Gy = Gyeval, rts = :VRS, Xref = Xref, Yref = Yref), :Y)
    end

    @test deaddfobs_ref_eff ≈ efficiency(deaddfobs)

    @test deaddfobsvs_ref_eff ≈ efficiency(deaddfobsvrs)

    @test deaddfvrs_ref_slackX ≈ slacks(deaddfobsvrs, :X) atol=1e-14
    @test deaddfvrs_ref_slackY ≈ slacks(deaddfobsvrs, :Y) atol=1e-15

    # Print
    show(IOBuffer(), deaddfobs)
    show(IOBuffer(), deaddfnoslack)

    # Test errors
    @test_throws DimensionMismatch deaddf([1; 2 ; 3], [4 ; 5], Gx = [1; 2 ; 3], Gy = [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch deaddf([1; 2], [4 ; 5], Gx = [1; 2], Gy = [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deaddf([1 1; 2 2], [4 4; 5 5], Gx = [1 1; 2 2], Gy = [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deaddf([1 1; 2 2], [4 4; 5 5], Gx = [1 1; 2 2], Gy = [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deaddf([1; 2; 3], [4; 5; 6], Gx = [1; 2; 3], Gy = [4; 5; 6], rts = :Error) # Invalid returns to scale
    @test_throws DimensionMismatch deaddf([1 1; 2 2; 3 3], [4; 5; 6], Gx = [1 1 1; 2 2 2; 3 3 3], Gy = [4; 5; 6]) # Different size of inputs direction
    @test_throws DimensionMismatch deaddf([1; 2; 3], [4 4; 5 5; 6 6], Gx = [1; 2; 3], Gy = [4 4 4; 5 5 5; 6 6 6]) # Different size of inputs direction
    @test_throws ArgumentError deaddf([1; 2; 3], [1; 2; 3], Gx = :Error, Gy = :Ones) # Invalid inuts direction
    @test_throws ArgumentError deaddf([1; 2; 3], [1; 2; 3], Gx = :Ones, Gy = :Error) # Invalid outputs direction

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------
    # Tests against results in R

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deaddf(X, Y, Gx = X, Gy = Y)) ≈ [0; 0; 0; 0.25; 0.4285714286; 0; 0.2; 0.2307692308]

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    @test efficiency(deaddf(X, Y, Gx = X, Gy = Y)) ≈ [0; 0; 0; 0.2173913043; 0.4; 0; 0.12; 0.2307692308]

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    @test efficiency(deaddf(X, Y, Gx = X, Gy = Y)) ≈ [0.4285714286; 0; 0.1111111111; 0.25; 0.4285714286; 0.4285714286; 0.3207547170; 0.6666666667]

    # ------------------
    # Test model with zero Directions
    # ------------------
    deamddfzeros = deaddf(X, Y, Gx = zeros(size(X)), Gy = zeros(size(Y)), rts = :VRS)
    @test efficiency(deamddfzeros) == zeros(8)

end


text:
# global PyObject constants that get initialized at runtime.  We
# initialize them here (rather than via "global foo = ..." in __init__)
# so that their type is known at compile-time.

const PY_THERMO = PyNULL()
const PY_CHEM = PyNULL()
const PY_COPY = PyNULL()

function __init__()
    copy!(PY_THERMO, pyimport_conda("thermo", "thermo", "conda-forge"))
    copy!(PY_CHEM, pyimport("thermo.chemical"))
    copy!(PY_COPY, pyimport("copy"))
end

ext: jl


text:
using Base: prompt, getpass, SecretBuffer, SecretBuffer!, shred!
using ..Common

mutable struct PKISecHandler <: SecHandler
    subtype::CosName
    version::Tuple{Int, Int}
    length::Int
    r::Int
    o::Vector{UInt8}
    u::Vector{UInt8}
    oe::Vector{UInt8}
    ue::Vector{UInt8}
    p::UInt32
    perms::Vector{UInt8}
    id::Vector{UInt8}
    cf::Union{CosDict, CosNullType}
    stmf::CosName
    strf::CosName
    eff::CosName
    access::Function
    skey_path::String
    iv_path::String
    keys::Dict{CosName, Tuple{UInt32, Vector{UInt8}}}
end

function PKISecHandler(doc::CosDoc, subtype::CosName, access::Function)
    enc = doc.encrypt
    length = get(get(enc, cn"Length", CosInt(128)))
    vrsn   = doc.version
    subtype in [cn"adbe.pkcs7.s3", cn"adbe.pkcs7.s4", cn"adbe.pkcs7.s5"] ||
        error("Unsupported release of security handler")

    r = subtype === cn"adbe.pkcs7.s3" ? 2 :
        subtype === cn"adbe.pkcs7.s4" ? 3 : 4
    o      = Vector{UInt8}(get(enc, cn"O"))
    u      = Vector{UInt8}(get(enc, cn"U"))
    oe     = Vector{UInt8}(get(enc, cn"OE"))
    ue     = Vector{UInt8}(get(enc, cn"UE"))
    p      = reinterpret(UInt32, Int32(get(get(enc, cn"P", CosInt(0)))))
    perms  = Vector{UInt8}(get(enc, cn"Perms"))
    ids = cosDocGetID(doc)
    id = Vector{UInt8}(ids[1])

    cf   = get(enc, cn"CF", CosDict())
    if r < 4
        d = CosDict()
        set!(cf, cn"DefaultCryptFilter", d)
        set!(d, cn"Type", cn"CryptFilter")
        set!(d, cn"CFM", cn"V2")
        set!(d, cn"AuthEvent", cn"DocOpen")
        set!(d, cn"Recipients", get(enc, cn"Recipients"))
        set!(d, cn"EncryptMetadata", get(enc, cn"EncryptMetadata", CosTrue))
        set!(d, cn"Length", get(enc, cn"Length", CosInt(128)))
        stmf = get(enc, cn"StmF", cn"DefaultCryptFilter")
        strf = get(enc, cn"StrF", cn"DefaultCryptFilter")
    else
        stmf = get(enc, cn"StmF", cn"Identity")
        strf = get(enc, cn"StrF", cn"Identity")
    end
    eff  = get(enc, cn"EFF", stmf)
    eff !== stmf &&
        Base.@warn("Embedded file streams without Crypt filters may not decrypt")
    access === identity && (access = get_digital_id)
    skey_path, io = get_tempfilepath()
    # Ensure the files have enough bytes in them.
    shred!(SecretBuffer!(crypto_random(rand(1000:10000)))) do s
        write(io, s)
        close(io)
    end
    iv_path, io = get_tempfilepath()
    # Ensure the files have enough bytes in them.
    shred!(SecretBuffer!(crypto_random(rand(1000:10000)))) do s
        write(io, s)
        close(io)
    end
    keys = Dict(cn"Identity" => (0xffffffff, Vector{UInt8}()))
    return PKISecHandler(subtype, vrsn, length, r, o, u, oe, ue,
                         p, perms, id, cf, stmf, strf, eff, access,
                         skey_path, iv_path, keys)
end

function get_digital_id()
    p12file = ""
    while !isfile(p12file)
        p12file =
            prompt("Select the PCKS#12 (.p12) certificate for the recepient")
    end
    p12pass = getpass("Enter the password to open the PKCS#12 (.p12) file")
    return shred!(x->read_pkcs12(p12file, x), p12pass)
end

function get_key(h::PKISecHandler, cfn::CosName, access::Function)
    cert, pkey = access()
    cf = get(h.cf, cfn)
    rs = get(cf, cn"Recipients")
    data = nothing
    flags = Cint(CMS_BINARY)
    digbuf, rbuf, perm = nothing, IOBuffer(), UInt32(0)
    try
        for r in get(rs)
            rb = Vector{UInt8}(r)
            write(rbuf, rb)
            data !== nothing && continue
            ci = CMSContentInfo(rb)
            data = decrypt(ci, pkey, cert, Vector{UInt8}(undef, 0), flags)
            if data !== nothing
                if length(data) > 20
                    bperm = data[21:24]
                    ENDIAN_BOM == 0x04030201 && reverse!(bperm)
                    perm = reinterpret(UInt32, bperm)[1]
                end
                digbuf = SecretBuffer!(data)
                digbuf.size = 20
                seekend(digbuf)
            end
        end
        data === nothing &&
            error("Unable to decrypt with recepient certificate")
        write(digbuf, take!(rbuf))
        encMetadata = get(get(cf, cn"EncryptMetadata", CosTrue))
        !encMetadata && write(digbuf, fill(0xff, 4))
        cfm = get_cfm(h, cfn)
        algo = cfm === cn"AESV3" ? "sha256" : "sha1"
        mdctx = DigestContext(algo)
        update!(mdctx, (@view digbuf.data[1:digbuf.size]))
        key = SecretBuffer!(close(mdctx))
        h.length = get(get(cf, cn"Length"))
        key.size = h.length / 8
        return store_key(h, cfn, (perm, key))
    finally
        digbuf !== nothing && shred!(digbuf)
    end
end


text:
dupdate = df[nonunique(df[:,[:Date]]),:][:Date]
println("The following rows have duplicate DATESTAMP:")
println(df[df[:Date] .== dupdate,:])
println("All values good in these rows:")
println(df[df[:ValidValues] .== 24,:])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# Applcation wide configuration. It contains default values for
# different settings. You may change these settings using val()
# function. Also may be used like intermodule data sharing mechanism.
# One module may set a value and other module may read it later.
# Every configuration value may be obtained by unique key. Some
# settings are used like data containers in sense that, they
# just store the values, which may be changed many times during
# application life cicle.
# TODO: update usage
# Usage:
#     import Config
#     ...
#     cfg = Config.create()                       # creates new config instance
#     cfg.SECTION_ID_KEY_ID                       # returns value or nothing
#     cfg.SECTION_ID_KEY_ID = newVal              # sets new value - newVal
#     Config.save(cfg, "config.data")             # saves all to file
#     Config.load("config.data")                  # loads all from file
#
# @singleton
# @author DeadbraiN
#
module Config
  import Helper
  import CommandLine

  export ORGANISM_MAX_MUTATION_PERIOD

  export if_debug
  export if_test
  export if_not_test
  export if_profile
  export if_not_profile

  export create
  export save
  export load
  export isEmpty
  export format
  export merge

  export ConfigData

  export MODE_QUIET_ALL
  export MODE_QUIET_IMPORTANT
  export MODE_QUIET_NO
  #
  # Show all messages
  #
  const MODE_QUIET_ALL = 0
  #
  # Show only important messages
  #
  const MODE_QUIET_IMPORTANT = 1
  #
  # Supress showing messages
  #
  const MODE_QUIET_NO = 2
  #
  # Maximum period for mutations. Related to orgRainMutationPeriod config
  #
  const ORGANISM_MAX_MUTATION_PERIOD = 1000
  #
  # Data type for storing configuration data. Is used in pair with GData
  # type. For accessing use Gonfig.val(man.cfg, :SYMBOL[, value]) function
  #
  type ConfigData
    #
    # Comman line arguments returned by CommanLine.create() function
    #
    cmdLineArgs::Dict{String, String}
    #
    # Probabilities with which mutator decides what to do: add,
    # change, delete character of the code; change amount of
    # mutations or change mutations period... Depending on these
    # values, organism may have different strategies of living.
    # For example: if add value is bigger then del and change,
    # then code size will be grow up all the time. If del value is
    # bigger then other, then it will be decreased to zero lines
    # of code and will die.
    # Format: [
    #     add          - Probability of adding of new character to the code
    #     change       - Probability of changing existing character in a code
    #     small-change - Probability of "small change" - change of expression part
    #     delete       - Probability of deleting of a character in a code
    #     clone        - Probability for amount of mutations on clone
    #     period       - Probability of period of organism mutations
    #     amount       - Probability of amount of mutations per period
    # ]
    #
    orgMutationProbs::Array{Int, 1}
    #
    # Max value, which we may use in orgMutationProbs array.
    #
    orgMutationProbsMaxValue::Int
    #
    # Percent of mutations from code size, which will be applied to
    # organism after clonning. Should be <= 1.0
    #
    orgCloneMutation::Float64
    #
    # Amount of iterations before clonning process
    #
    orgClonePeriod::Int
    #
    # Amount of iterations within organism's life loop, after that we
    # do mutations according to orgRainMutationPercent config. If 0, then
    # mutations will be disabled. Should be less then ORGANISM_MAX_MUTATION_PERIOD
    #
    orgRainMutationPeriod::Int
    #
    # Value, which will be used like amount of mutations per
    # orgRainMutationPeriod iterations. 0 is a possible value if
    # we want to disable mutations. Should be less then 100
    #
    orgRainMutationPercent::Float64
    #
    # Amount of organisms we have to create on program start
    #
    orgStartAmount::Int
    #
    # Amount of energy for first organisms. They are like Adam and
    # Eve. It means that these empty (without code) organism were created
    # by operator and not by evolution.
    #
    orgStartEnergy::Int
    #
    # Begin color of "empty" organism (organism without code).
    #
    orgStartColor::Int
    #
    # Only after this amount of mutations organism should update it's color
    #
    orgColorPeriod::Int
    #
    # Amount of iterations within organism's life loop, after that we decrease
    # some amount of energy. If 0, then energy decreasing will be disabled.
    #
    orgEnergySpendPeriod::Int
    #
    # Amount of energy, which will be decreased in case of organis's
    # code error or exception
    #
    orgEnergySpendOnError::Int
    #
    # Amount of iterations when organism is alive. It will die after
    # this period. If 0, then will not be used.
    #
    orgAlivePeriod::Int
    #
    # This value means the period between organism codeSizes, which
    # affects energy grabbing by the system. For example: we have two
    # organisms: org1.energy = 10, org2.energy = 10, org1.codeSize = 6,
    # org2.codeSize = 9, Config.orgGarbagePeriod = 5. It means that
    # during energy grabbing by the system org1 and org2 will spend the
    # same amount of energy - 2 units. This is because the period goes
    # from 1..5, 6..10,... and both organisms are in the same period.
    #
    orgGarbagePeriod::Int
    #
    # Amount of errors in organisms codes in current population
    #
    orgErrors::Int
    #
    # Amount of eval calls for generatin organisms code
    #
    orgEvals::Int
    #
    # Maximum amount of arguments in custom functions. Minimum 1. Maximum
    # <= amount of default variables.
    #
    codeFuncParamAmount::Int
    #
    # Amount of iterations in a loop (for operator)
    #
    codeLoopAmount::Int
    #
    # Amount of organisms entire code runs per some period of time
    #
    codeRuns::Int
    #
    # If organism reach this limit of amount of code lines, then codeSizeCoef
    # will be used during it's energy grabbing by system. We use this approach,
    # because our CPU's are slow and organisms with big codes are very slow. But
    # it's possible for organisms to go outside the limit by inventing new
    # effective mechanisms of energy obtaining.
    #
    codeMaxSize::Int
    #
    # This coefficiend is used for calculating of amount of energy,
    # which grabbed from each organism depending on his codeSize.
    # This coefficient affects entire code size of population and
    # entire system speed. It depends on CPU speed also. So, for
    # different PC's it may be different.
    # Formula is the following: grabEnergy = cfg.codeSizeCoef * org.codeSize
    # See Config.codeMaxSize for details. This config will be turn on only if
    # organism reaches code size limit Config.codeMaxSize
    #
    codeSizeCoef::Int
    #
    # World width
    #
    worldWidth::Int
    #
    # World height
    #
    worldHeight::Int
    #
    # Turns on ciclic world mode. It means that organisms may go outside
    # it's border, but still be inside. For example, if the world has 10x10
    # size and the organism has 10x5 position in it, one step right will move
    # this organism at the position 1x5. The same scenario regarding Y
    # coordinate (height).
    worldCyclical::Bool
    #
    # Maximum amount of organisms in a world. If some organisms will
    # try to clone itself, when entire amount of organisms are equal
    # this value, then it(clonning) will not happen.
    #
    worldMaxOrgs::Int
    #
    # Minimum amount of orgaisms in a world. If this value riached,
    # then remove minimum energetic organisms mechanism will be disabled
    # until total amount will be more then this value.
    #
    worldMinOrgs::Int
    #
    # Amount of energy blocks in a world. Blocks will be placed in a
    # random way...
    #
    worldStartEnergyDots::Int
    #
    # Amount of energy in every block. See worldStartEnergyDots
    # config for details.
    #
    worldStartEnergyInDot::UInt16
    #
    # Minimum percent of energy in current world. Under percent i mean
    # percent from entire world area (100%). If the energy will be less
    # or equal then this percent, then new random energy should be added.
    # Should be less then 100.0 and more and equal to 0.0. 0.17 is a
    # normal percent for this system.
    #
    worldEnergyCheckPercent::Float64
    #
    # An amount of iteration, after which we have to check world energy
    # amount. Works in pair with worldEnergyCheckPercent. May be 0 if
    # you want to disable it
    #
    worldEnergyCheckPeriod::Int
    #
    # World scaling. On todays monitors pixel are so small, so we have
    # to zoom them with a coefficient.
    #
    worldZoom::Int
    #
    # Period of seconds, which is user for checking IPS value. It's
    # possible to increase it to reduce amount of requests and additional
    # code in main loop
    #
    worldIpsPeriod::Float64
    #
    # Period of making automatic backup of application. In seconds
    #
    backupPeriod::Float64
    #
    # Amount of backup files stored on HDD. Old files will be removed
    #
    backupAmount::Int
    #
    # The period of time between yield() calls in "stand by" mode.
    # In this mode manager waits for data in sockets and new connections.
    # In this mode yield() is called only once in a period, because
    # it eats CPU cicles. In case of data in sockets or new connections
    # yield() will be called more often.
    #
    conYieldPeriod::Float64
    #
    # Percent of energy, which will be minused from organism after
    # stepping from one instance to another.
    #
    conStepEnergySpendPercent::Int
    #
    # Starting number for TCP/IP listening
    #
    conServerPort::Int
    #
    # Works in pair with conServerPort. An IP of current
    # server/instance.
    # TODO: IPv6?
    conServerIp::IPv4
    #
    # Port number for "fast" mode. It uses, for example, for pooling
    #
    conFastServerPort::Int
    #
    # Left side server's (instance) port we want connect to. May be
    # zero (0) if no left side server available.
    #
    conLeftServerPort::Int
    #
    # Left server(instance) IP address. Works in pair with
    # conLeftServerPort
    #
    conLeftServerIp::IPv4
    #
    # Right side server's (instance) port we want connect to. May be
    # zero (0) if no right side server available.
    #
    conRightServerPort::Int
    #
    # Right server(instance) IP address. Works in pair with
    # conRightServerPort
    #
    conRightServerIp::IPv4
    #
    # Left up server's (instance) port we want connect to. May be
    # zero (0) if no up side server available.
    #
    conUpServerPort::Int
    #
    # Up server(instance) IP address. Works in pair with
    # conUpServerPort
    #
    conUpServerIp::IPv4
    #
    # Left down server's (instance) port we want connect to. May be
    # zero (0) if no down side server available.
    #
    conDownServerPort::Int
    #
    # Down server(instance) IP address. Works in pair with
    # conDownServerPort
    #
    conDownServerIp::IPv4
    #
    # Array of included plugins
    #
    plugIncluded::Array{String, 1}
    #
    # Array of excluded plugins. Affects plugIncluded list
    #
    plugExcluded::Array{String, 1}
    #
    # Debug mode. This mode means, that all debug messages
    # will be posted to the terminal
    #
    modeDebug::Bool
    #
    # Testing mode. In this mode user may run jevo step by step
    # and test it'sinternal parts. For example, during unit tests
    #
    modeTest::Bool
    #
    # Is used for profiling the application with ProfileView
    # package. See run-profiling.sh for details
    #
    modeProfile::Bool
    #
    # Amount of iterations in profile mode after which ProfileView
    # package will draw performance flame chart
    #
    modeProfilePeriod::Int
    #
    # Amount of seconds for status showing in terminal
    #
    modeStatusPeriod::Int
    #
    # In this mode status report will be short or full
    #
    modeStatusFull::Bool
    #
    # Mode for showing/supressing of messages. Possible values:
    #   0 - all messages
    #   1 - only important messages
    #   2 - no messages
    #
    modeQuiet::Int
  end
  #
  # Returns one configuration setting from command line parameters
  # or from ConfigData instance. Default values are set in create() method
  # @param args Command line arguments
  # @param name Name of the setting
  # @return {String} Value from command line or default value
  #
  function _setting(data::ConfigData, name::Symbol)
    local sName::String = string(name)
    local val::Any

    if CommandLine.has(data.cmdLineArgs, sName)
      val = CommandLine.val(data.cmdLineArgs, sName, true)
      #
      # In case of Int type we have to figure out exact type
      # Int8, UInt16 and so on...
      #
      if typeof(val) === Int return convert(typeof(getfield(data, name)), val) end
      return val
    end

    getfield(data, name)
  end
  #
  # Merges default configuration obtained by calling create() method
  # and command line arguments. Command line arguments have higher
  # priority, then default. For example, if use runs jevo like this:
  #
  #   julia src/AppManager.jl worldWidth:1024
  #
  # it means, that worldWidth config will be 1024 instead of default value
  # @param data Configuration data object
  # @return {ConfigData} Updated configuration object
  #
  function merge(data::ConfigData)
    local name::Symbol

    for name in fieldnames(data)
      #
      # We have to skip command line arguments object
      #
      if name === :cmdLineArgs continue end
      setfield!(data, name, _setting(data, name))
    end

    data
  end
  #
  # Creates configuration object. It will be used in all config
  # related functions like val(), save(), load(),... We don't need
  # to use global objects, because they are very slow!
  # @param empty Flag to create empty configuration
  # @return {ConfigData} Data type instance
  #
  function create(isMerge::Bool = true, empty::Bool = false)
    local cfg::ConfigData = ConfigData(
      CommandLine.create(),                    # cmdLineArgs
      empty ? [] : [50,100,1,0,1,1,1,1,1],     # orgMutationProbs (add,change,del,small-change,clone,period,amount,probs,cloneEnergy)
      100,                                     # orgMutationProbsMaxValue
      0.02,                                    # orgCloneMutation
      20,                                      # orgClonePeriod
      0,                                       # orgRainMutationPeriod
      0.01,                                    # orgRainMutationPercent
      500,                                     # orgStartAmount
      10000,                                   # orgStartEnergy
      100,                                     # orgStartColor
      50,                                      # orgColorPeriod
      500,                                     # orgEnergySpendPeriod
      1000,                                    # orgEnergySpendOnError
      10000,                                   # orgAlivePeriod (amountOfSeconds * averageIPSperSecond)
      7,                                       # orgGarbagePeriod
      0,                                       # orgErrors
      0,                                       # orgEvals

      2,                                       # codeFuncParamAmount
      8,                                       # codeLoopAmount
      0,                                       # codeRuns
      21,                                      # codeMaxSize
      20,                                      # codeSizeCoef

      1900,                                    # worldWidth
      940,                                     # worldHeight
      true,                                    # worldCyclical
      1000,                                    # worldMaxOrgs
      0,                                       # worldMinOrgs
      1000,                                    # worldStartEnergyDots
      UInt16(0x0001F4),                        # worldStartEnergyInDot
      0.3,                                     # worldEnergyCheckPercent
      20000,                                   # worldEnergyCheckPeriod
      1,                                       # worldZoom
      1.0,                                     # worldIpsPeriod

      300.0,                                   # backupPeriod
      10,                                      # backupAmount

      0.01,                                    # conYieldPeriod
      20,                                      # conStepEnergySpendPercent
      2010,                                    # conServerPort (current server port)
      ip"127.0.0.1",                           # conServerIp
      2011,                                    # conFastServerPort (current server "fast" mode port)
      0,                                       # conLeftServerPort
      ip"127.0.0.1",                           # conLeftServerIp
      0,                                       # conRightServerPort
      ip"127.0.0.1",                           # conRightServerIp
      0,                                       # conUpServerPort
      ip"127.0.0.1",                           # conUpServerIp
      0,                                       # conDownServerPort
      ip"127.0.0.1",                           # conDownServerIp

      [#="Phylogen", "Status"=#],              # plugIncluded
      [],                                      # plugExcluded

      false,                                   # modeDebug
      false,                                   # modeTest
      false,                                   # modeProfile
      2000,                                    # modeProfilePeriod
      10.0,                                    # modeStatusPeriod
      false,                                   # modeStatusFull
      0                                        # modeQuiet
    )

    isMerge ? merge(cfg) : cfg
  end
  #
  # This macro turns on DEBUG code if modeDebug. For example, it
  # shows debug messages in terminal. Don't change this macro position
  #
  macro if_debug(ex) @static if _setting(create(), :modeDebug) esc(ex) end end
  #
  # This macro turns on special test mode, which is called "iterational".
  # In this mode we may run Manager for special amount of iterations like a task.
  #
  macro if_test(ex) @static if _setting(create(), :modeTest) esc(ex) end end
  #
  # This macro turns on if special test mode is off.
  #
  macro if_not_test(ex) @static if !_setting(create(), :modeTest) esc(ex) end end
  #
  # This macro turns on specified code in profiling mode. You may use
  # run-profiling.sh for this
  #
  macro if_profile(ex) @static if _setting(create(), :modeProfile) esc(ex) end end
  #
  # This macro turns on specified code in not profiling mode. You may use
  # run-server.bat for this.
  #
  macro if_not_profile(ex) @static if !_setting(create(), :modeProfile) esc(ex) end end
  #
  # Saves all data into the file. If file exists, it will be overriden
  # @param data Config Data type
  # @param file File name
  # @return {Bool} saving result
  #
  function save(data::ConfigData, file::String = "config.data")
    Helper.save(data, file)
  end
  #
  # Loads all data from the file. If something went wrong new,
  # default configuration data type will be created.
  # @param file File name
  # @return {Tuple{ConfigData, Bool}} loading result
  #
  function load(file::String = "config.data")
    local ret = nothing
    try
      ret = Helper.load(file)
    catch e
      Helper.warn("Config.load(): $e")
      @if_debug showerror(STDOUT, e, catch_backtrace())
    end

    typeof(ret) !== ConfigData ? create(true, true) : ret
  end
  #
  # Checks if config data type is empty
  # @param data Config Data type
  # @return {Bool}
  #
  function isEmpty(data::ConfigData)
    length(data.orgMutationProbs) == 0
  end
  #
  # Formats configuration for usable user's view
  # @return {String}
  # TODO: do we need this?
  function format(data::ConfigData)
    local fields::Array{Symbol, 1} = fieldnames(ConfigData)
    local i::Symbol
    local arr::Array{String, 1} = []

    for i in fields
      push!(arr, string(i) * ": " * string(getfield(data, i)))
    end

    arr
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
type PCA
    solver::String
    n_components::Integer
    components::Matrix
    mean_::Vector
end

function PCA(;
             solver::String = "svd",
             n_components::Integer = 2,
             components::Matrix = zeros(4,4),
             mean_::Vector = zeros(4))
    return PCA(solver, n_components, components, mean_)
end

function train!(model::PCA, X::Matrix)
    model.mean_ = vec(mean(X,1))
    n_sample = size(X, 1)
    X_de_mean = X - repmat(model.mean_', n_sample, 1)

    if model.solver == "svd"
        U,S,V = svd(X_de_mean)
    elseif model.solver == "eig"
        cov_ = cov(X_de_mean)
        D,V = eigs(cov_, nev = model.n_components)
    end
    model.components = V[:, 1:model.n_components]
end

function transform(model::PCA,
                 x::Matrix)
    n = size(x,1)
    res = zeros(n, model.n_components)
    for i = 1:n
        res[i, :] = transform(model, x[i,:])
    end
    return res
end

function transform(model::PCA,
                 x::Vector)
    x = x - model.mean_
    x = vec(x' * model.components)
    return x
end

function plot_in_2d(model::PCA, X::Matrix, y::Vector, title::String)
    X = transform(model, X)
    x1 = X[:, 1]
    x2 = X[:, 2]
    df = DataFrame(x = x1, y = x2, clu = y)
    println("Computing finished")
    println("Drawing the plot.....Please Wait(Actually Gadfly is quite slow in drawing the first plot)")
    Gadfly.plot(df, x = "x", y = "y", color = "clu", Geom.point, Guide.title(title))
end

function test_PCA()
    X_train, X_test, y_train, y_test = make_digits()
    model = PCA()
    train!(model,X_train)
    plot_in_2d(model, X_train, y_train, "PCA")
end


text:
@diffusion_process FitzHughDiffusionAltAuxLin begin
    :dimensions
    process --> 2
    wiener --> 1

    :parameters
    (ϵ, s, γ, β, σ, t, T) --> Float64
    (u, v) --> ℝ{2}

    :additional
    constdiff --> true
    linear --> true
end

function B(t, P::FitzHughDiffusionAltAuxLin)
    @SMatrix [
        0.0  1.0;
        (1.0-P.γ-3.0*P.v[1]^2-6*P.v[1]*P.v[2])/P.ϵ (1.0-P.ϵ-3.0*P.v[1]^2)/P.ϵ
    ]
end

function β(t, P::FitzHughDiffusionAltAuxLin)
    ℝ{2}(
        0.0,
        (2*P.v[1]^3+P.s-P.β+6*P.v[1]^2*P.v[2])/P.ϵ
    )
end

σ(t, P::FitzHughDiffusionAltAuxLin) = ℝ{2}(0.0, P.σ/P.ϵ)


text:
using CuArrays
include("../../src/core.jl")

using MLDataUtils
using MLDatasets

using CUDNN2

include("model.jl")

@diffrule CUDAnative.log(x::Real) x ds / x
@diffrule CUDAnative.log(x::AbstractArray) x ds ./ x

function main()
    data, labels = MNIST.traindata()
    X_train = reshape(data, 28, 28, 1, size(data, 3))
    y_train = one_hot(maximum(labels) + 1, labels .+ 1)
    m = C2(randn(5, 5, 1, 3), randn(5, 5, 3, 6), randn(10, 4*4*6), randn(10))
    m = fit!(m, X_train, y_train; cuda=true, n_epochs=700)
    # test accuracy on a single batch
    X = to_cuda(X_train[:,:,:,101:200])
    y = to_cuda(y_train[:,101:200])
    accuracy(predict(m, X), y)
end


text:
module Coordinates

export AbstractCoordinates
export Cartesian, Polar

abstract type AbstractCoordinates end

struct Cartesian <: AbstractCoordinates end
struct Polar <: AbstractCoordinates end

end # module


text:
module MadNLPGraph

import MadNLP:
    optimize!, Optimizer, InteriorPointSolver,
    MOIU, MOI, get_nnz_hess, get_nnz_jac, set_x!, set_g!, INITIAL, is_jac_hess_constant,
    jacobian_structure, hessian_lagrangian_structure, eval_objective,
    eval_objective_gradient, eval_function, eval_constraint, eval_hessian_lagrangian, eval_constraint_jacobian,
    @kwdef, Logger, @debug, @warn, @error, @sprintf,
    AbstractOptions, AbstractLinearSolver, EmptyLinearSolver, set_options!,
    SparseMatrixCSC, SubVector, StrideOneVector,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia,
    default_linear_solver, default_dense_solver, get_csc_view, get_cscsy_view, nnz, mul!,
    set_blas_num_threads, blas_num_threads, @blas_safe_threads,
    MOIModel, NLPModelsCounters, MadNLPExecutionStats, AbstractNLPModel, NLPModelMeta,
    get_lcon, get_ucon, jac_structure!, hess_structure!, obj, grad!, cons!, jac_coord!, hess_coord!
import Metis: partition
import LightGraphs: Graph, Edge, add_edge!, edges, src, dst, neighbors, nv
import Plasmo: OptiGraph, OptiNode, OptiEdge, all_nodes, all_edges, all_variables, num_all_nodes, link_constraints, getnode, num_variables, num_constraints
import JuMP: _create_nlp_block_data, set_optimizer, GenericAffExpr

include("plasmo_interface.jl")
include("graphtools.jl")
include("schur.jl")
include("schwarz.jl")

export MadNLPSchur, MadNLPSchwarz

end # module


text:
import PredictMDAPI
import Test


text:
@testset "complex_powers" begin
    complex_powers_comparison(z, m, T=Float64) = (
        complex_powers(Complex{T}(z), m),
        z.^collect(0:m),
        eps(T)*m
    )

    for T in [Float64, Float32, Float16]
        fudge = one(T) + 2 * sqrt(eps(T))
        for k in 0:25
            z = exp(k*big(π)im/10)
            for m in [0, 1, 2, 3, 4, 1_000]
                mine, theirs, ϵ = complex_powers_comparison(z, m, T)
                @test mine ≈ theirs rtol=2ϵ
                @test_throws DomainError complex_powers(Complex{T}(z*fudge), m)
                @test_throws DomainError complex_powers(Complex{T}(z/fudge), m)
            end
        end
    end
end


text:
import Base.Libc: TmStruct
import SeisIO: buf_to_uint,
  check_for_gap!,
  dtconst,
  endtime,
  is_u8_digit,
  mk_t,
  split_id,
  stream_float,
  sμ,
  μs


text:
using Test
using Logging
# using Revise
using LarSurf
# using Plasm
# Logging.configure(level==Logging.Debug)

# include("../src/LarSurf.jl")
# include("../src/block.jl")


text:
#############################################################################
# weights_of_graphs.jl
# handles the weights of AbstractGraph, complementary of LightGraphs.
#############################################################################

export WeightedGraph
export add_edge!, vertice_weights, edge_weights
export nv, ne, edges, induced_subgraph, connected_components

struct WeightedGraph
  vertice_weights::Dict{Int64, Number}
  edge_weights::Dict{Tuple{Int64,Int64},Number}
  graph::AbstractGraph

  function WeightedGraph(graph::AbstractGraph, vertice_weights::Dict{Int64, Number} = Dict{Int64, Number}())
    for i = 1:nv(graph)
      if in(i, keys(vertice_weights)) == false
        push!(vertice_weights, i => 1)
      end
    end
    return new(vertice_weights, Dict{Tuple{Int64,Int64}, Number}(), graph)
  end
end

WeightedGraph(n::Int) = WeightedGraph(Graph(n))

WeightedGraph(n::Int, vertice_weights::Dict{Int64, Number}) = WeightedGraph(Graph(n), vertice_weights)

function add_edge!(wg::WeightedGraph, src::Int, dst::Int, weight::Number = 1)
  add_edge!(wg.graph, src, dst)
  push!(wg.edge_weights, (src, dst) => weight)
end

function vertice_weights(wg::WeightedGraph)
  return wg.vertice_weights
end

function edge_weights(wg::WeightedGraph)
  return wg.edge_weights
end

nv(wg::WeightedGraph) = nv(wg.graph)

ne(wg::WeightedGraph) = ne(wg.graph)

edges(wg::WeightedGraph) = edges(wg.graph)

induced_subgraph(wg::WeightedGraph, elist) = induced_subgraph(wg.graph, elist)

connected_components(wg::WeightedGraph) = connected_components(wg.graph)


text:
# Binary search for zero of continuous function f
# with one zero-crossing in [lo, hi]
function binary_search(f, lo, hi; ϵ = 1e-10, maxiter = 100)
    if f(lo) > 0
        @warn "f($(lo))=$(f(lo)) should be negative at low end";
        return lo;
    end
    if f(hi) < 0
        @warn "f($(hi))=$(f(hi)) should be positive at high end";
        return hi;
    end

    # @assert f(lo) < +ϵ "f($(lo))=$(f(lo)) should be negative at low end";
    # @assert f(hi) > -ϵ "f($(hi))=$(f(hi)) should be positive at high end";

    for i in 1:maxiter
        mid = (lo+hi)/2;
        if mid == lo || mid == hi
            # not going to get any better ...
            #@warn "binary_search collapsed after $i iterations.\nf($(lo)) = $(f(lo))\nf($(hi)) = $(f(hi))";
            return mid;
        end
        fmid = f(mid);
        if fmid < -ϵ
            lo = mid
        elseif fmid > ϵ
            hi = mid;
        else
            return mid;
        end
    end
    @warn "binary_search did not reach tolerance $ϵ in $maxiter iterations.\nf($(lo)) = $(f(lo))\nf($(hi)) = $(f(hi)),\nmid would be $((lo+hi)/2)";
    return (lo+hi)/2;
end


text:
# This pkey is just a randomly generated key so we have something to use. This isn't used
# for authenticating anywhere.
const TEST_PKEY = """
-----BEGIN OPENSSH PRIVATE KEY-----
b3BlbnNzaC1rZXktdjEAAAAABG5vbmUAAAAEbm9uZQAAAAAAAAABAAAAlwAAAAdzc2gtcn
NhAAAAAwEAAQAAAIEA3nPemZYYdXxCpkn9z/50Pk2LO5+8m4FkmQBKUY0arKlsymrRy88O
dYu9Y9fiu3WplgG2PLpqj2Vd1SJq7DzdUOEmZ9Vp47rO2Tx/bkRn1AfqLGS7AYwMMejqdd
rC03o/6m1mEQdI4AblLiK2U/gWOe3dZdiq1AuhFD6Bmp/+7DkAAAIQ+Ce3m/gnt5sAAAAH
c3NoLXJzYQAAAIEA3nPemZYYdXxCpkn9z/50Pk2LO5+8m4FkmQBKUY0arKlsymrRy88OdY
u9Y9fiu3WplgG2PLpqj2Vd1SJq7DzdUOEmZ9Vp47rO2Tx/bkRn1AfqLGS7AYwMMejqddrC
03o/6m1mEQdI4AblLiK2U/gWOe3dZdiq1AuhFD6Bmp/+7DkAAAADAQABAAAAgQCSs69Vcl
rm/++kYp90D8bxX4o24/0qQIbrL/nTFk9FFuacBx4cXoyWkHVx5umr3sjcGHzqR7YGoz7i
VDRXXzxD47SRdYJ0IHvUB8dCpZHJt6qDxndMUvp5z/eLccAATm77L1i3ve3D4Eferj+mvR
JV+ZLzLT/e3MBzjdkgcXZDQQAAAEANLoS7phEEVMxZKiC03q0NNmTUFmx85P3TiDSiXtYe
7ZlAK9fns0873twevFWeVBYOQ30QhEJc1mKanjEvXcc1AAAAQQDwNqt5yA6KvzdV3dH1ZX
rTwWx1CO5ERpuMSOAQhk3y6hfmeLa5HnbSCP5Z3uVqizEDnGXkvRfhSyfEUu/jJTU1AAAA
QQDtEmQvWdgz+HtIuTG1ySJ9FYO6LeCEXHtQX78aOfNaj2jqLTXHdqrMr0V5exJcNV4XSc
8e9dZXl8OX+9Ub0Y91AAAAGmZjaG9ybmV5QEZlcm5hbmRvcy1NQlAubGFu
-----END OPENSSH PRIVATE KEY-----
"""

@testset "git_push" begin
    function create_local_remote(dir::AbstractString)
        remote_path = joinpath(dir, "localremote.git")
        run(`git init --bare $remote_path`)

        return remote_path
    end

    pushed_str = "HEAD -> master, origin/master"

    @testset "basic push" begin
        mktempdir() do local_remote_dir
            mktempdir() do f
                local_remote_path = create_local_remote(local_remote_dir)

                cd(f) do
                    run(`git init`)
                    run(`git remote add origin $local_remote_path`)

                    run(`touch foobar.txt`)
                    CompatHelper.git_add()
                    CompatHelper.git_commit("Message")

                    output = read(`git log --decorate`, String)
                    @test !occursin(pushed_str, output)

                    CompatHelper.git_push("origin", "master")

                    output = read(`git log --decorate`, String)
                    @test occursin(pushed_str, output)
                end
            end
        end
    end

    @testset "force push" begin
        mktempdir() do local_remote_dir
            mktempdir() do f
                local_remote_path = create_local_remote(local_remote_dir)

                cd(f) do
                    run(`git init`)
                    run(`git remote add origin $local_remote_path`)

                    run(`touch foobar.txt`)
                    CompatHelper.git_add()
                    CompatHelper.git_commit("Message 1")

                    output = read(`git log --decorate`, String)
                    @test !occursin(pushed_str, output)

                    CompatHelper.git_push("origin", "master")

                    output = read(`git log --decorate`, String)
                    @test occursin(pushed_str, output)
                    @test occursin("Message 1", output)

                    run(`touch baz.txt`)
                    CompatHelper.git_add()
                    run(
                        `git -c user.name=user -c user.email=email commit --amend --no-edit -m "Message 2"`,
                    )

                    CompatHelper.git_push("origin", "master"; force=true)
                    output = read(`git log --decorate`, String)
                    @test occursin(pushed_str, output)
                    @test occursin("Message 2", output)
                end
            end
        end
    end

    @testset "SSH push" begin
        mktempdir() do local_remote_dir
            mktempdir() do f
                local_path = joinpath(f, CompatHelper.LOCAL_REPO_NAME)
                pkey = joinpath(f, "privatekey")

                open(pkey, "w") do io
                    print(io, TEST_PKEY)
                end

                local_remote_path = create_local_remote(local_remote_dir)

                cd(f) do
                    run(`git init`)
                    run(`git remote add origin $local_remote_path`)

                    run(`touch foobar.txt`)
                    CompatHelper.git_add()
                    CompatHelper.git_commit("Message")

                    output = read(`git log --decorate`, String)
                    @test !occursin(pushed_str, output)

                    CompatHelper.git_push("origin", "master", pkey)

                    output = read(`git log --decorate`, String)
                    @test occursin(pushed_str, output)
                end
            end
        end
    end
end

@testset "git_reset" begin
    mktempdir() do f
        cd(f) do
            run(`git init`)

            @test !isfile("foobar.txt")

            run(`touch foobar.txt`)
            CompatHelper.git_add()
            @test CompatHelper.git_commit("Message")

            @test !isfile("bazbar.txt")

            run(`touch bazbar.txt`)
            CompatHelper.git_add()
            @test CompatHelper.git_commit("Message2")

            @test isfile("foobar.txt")
            @test isfile("bazbar.txt")

            hash = read(`git rev-parse HEAD`, String)
            CompatHelper.git_reset("HEAD~1"; soft=true)

            @test isfile("foobar.txt")
            @test isfile("bazbar.txt")

            sleep(1) # To make sure we get a new timestamp
            @test CompatHelper.git_commit("Message2")

            new_hash = read(`git rev-parse HEAD`, String)
            @test hash != new_hash
        end
    end
end

@testset "git_commit" begin
    @testset "success" begin
        mktempdir() do f
            cd(f) do
                run(`git init`)
                run(`touch foobar.txt`)
                CompatHelper.git_add()

                @test CompatHelper.git_commit("Message")
            end
        end
    end

    @testset "failure" begin
        mktempdir() do f
            cd(f) do
                run(`git init`)
                run(`touch foobar.txt`)
                CompatHelper.git_add()

                # Manually create an index lock file, so that the commit fails
                run(`touch .git/index.lock`)

                @test !CompatHelper.git_commit("Message")
            end
        end
    end
end

@testset "git_branch" begin
    branch = "foobar"

    @testset "no checkout" begin
        mktempdir() do f
            cd(f) do
                run(`git init`)
                run(`touch foobar.txt`)
                CompatHelper.git_add()
                CompatHelper.git_commit("Message")

                output = strip(read(`git branch`, String))
                @test output == "* master"

                CompatHelper.git_branch(branch)

                output = strip(read(`git branch`, String))
                @test output == "$branch\n* master"

                CompatHelper.git_checkout(branch)

                output = strip(read(`git branch`, String))
                @test output == "* $branch\n  master"
            end
        end
    end

    @testset "with checkout" begin
        mktempdir() do f
            cd(f) do
                run(`git init`)
                run(`touch foobar.txt`)
                CompatHelper.git_add()
                CompatHelper.git_commit("Message")

                output = strip(read(`git branch`, String))
                @test output == "* master"

                CompatHelper.git_branch(branch; checkout=true)

                output = strip(read(`git branch`, String))
                @test output == "* $branch\n  master"
            end
        end
    end
end

@testset "git_add" begin
    untracked_str = "Untracked files"
    committed_str = "Changes to be committed"

    mktempdir() do f
        cd(f) do
            run(`git init`)
            run(`touch foo.txt`)
            run(`touch bar.txt`)

            output = read(`git status`, String)
            @test occursin(untracked_str, output)

            CompatHelper.git_add()

            output = read(`git status`, String)
            @test !occursin(untracked_str, output)
            @test occursin(committed_str, output)
        end
    end
end

@testset "get_git_name_and_email" begin
    @testset "using defaults" begin
        expected = (
            CompatHelper.COMPATHELPER_GIT_COMMITTER_NAME,
            CompatHelper.COMPATHELPER_GIT_COMMITTER_EMAIL,
        )

        results = CompatHelper.get_git_name_and_email()
        @test results == expected
    end

    @testset "setting env vars" begin
        expected = ("User", "Email")

        withenv("GIT_COMMITTER_NAME" => "User", "GIT_COMMITTER_EMAIL" => "Email") do
            results = CompatHelper.get_git_name_and_email()
            @test results == expected
        end
    end
end

@testset "git_clone" begin
    @testset "HTTPS" begin
        mktempdir() do f
            cd(f) do
                local_path = joinpath(f, CompatHelper.LOCAL_REPO_NAME)
                CompatHelper.git_clone(
                    "https://github.com/JuliaRegistries/CompatHelper.jl/", local_path
                )

                @test !isempty(readdir(local_path))
            end
        end
    end

    @testset "SSH" begin
        mktempdir() do f
            cd(f) do
                local_path = joinpath(f, CompatHelper.LOCAL_REPO_NAME)

                apply(make_ssh_clone_patch(local_path)) do
                    CompatHelper.git_clone(
                        "git@github.com:JuliaRegistries/CompatHelper.jl.git", local_path
                    )
                end

                @test !isempty(readdir(local_path))
            end
        end
    end
end

@testset "git_checkout" begin
    master = "master"

    mktempdir() do f
        cd(f) do
            run(`git init`)
            # Need to create a commit before hand, see below
            # https://stackoverflow.com/a/63480330/1327636
            run(`touch foobar.txt`)
            CompatHelper.git_add()
            CompatHelper.git_commit("Message")

            CompatHelper.git_checkout(master)
            result = String(read((`git branch --show-current`)))

            @test contains(result, master)
        end
    end
end

@testset "git_get_master_branch" begin
    @testset "default branch" begin
        branch = "master"
        mktempdir() do f
            cd(f) do
                run(`git init`)
                run(`touch foobar.txt`)
                CompatHelper.git_add()
                CompatHelper.git_commit("Message")
                CompatHelper.git_checkout(branch)

                @test CompatHelper.git_get_master_branch(CompatHelper.DefaultBranch()) ==
                    branch
            end
        end
    end

    @testset "master branch" begin
        branch = "main"
        mktempdir() do f
            cd(f) do
                run(`git init`)
                run(`git branch -m $branch`)
                run(`touch foobar.txt`)
                CompatHelper.git_add()
                CompatHelper.git_commit("Message")
                CompatHelper.git_checkout(branch)

                @test CompatHelper.git_get_master_branch(branch) == branch
            end
        end
    end
end


text:
@testset "evaluate.jl" begin
    @testset "evaluate erroring" begin
        @test_throws DexError evaluate("(1 : Int) + 2.0")
    end

    @testset "evaluate  show" begin
        @test repr(evaluate("1")) == repr("1")
        @test repr(evaluate("1.5")) == repr("1.5")
        @test repr(evaluate("[1, 2]")) == repr("[1, 2]")
        @test repr(evaluate("1+3")) == repr("4")
        @test repr(evaluate("for i. [1, 2].i + 1")) == repr("[2, 3]")

        # This seems weird: why is it doubly quoted? 😕
        @test repr(evaluate("IToW8 65")) === repr(repr('A'))
    end

    @testset "evaluate juliaize" begin
        @test juliaize(evaluate("1")) === Int32(1)
        @test juliaize(evaluate("1.5")) === 1.5f0
        @test juliaize(evaluate("IToW8 65")) === UInt8(65)
    end

    @testset "juliaize-dexize round-trip" begin
        @test juliaize(dexize(Int64(3))) === Int64(3)
        @test juliaize(dexize(Int32(3))) === Int32(3)
        @test juliaize(dexize(UInt8(3))) === UInt8(3)
        @test juliaize(dexize(Float64(3))) === Float64(3)
        @test juliaize(dexize(Float32(3))) === Float32(3)
        @test juliaize(dexize(UInt64(3))) === UInt64(3)
        @test juliaize(dexize(UInt32(3))) === UInt32(3)
    end

    @testset "Atom function call" begin
        m = DexModule("""
        def addOne (x: Float) : Float = x + 1.0
        """)
        x = evaluate("2.5")
        y = evaluate("[2, 3, 4]")
        @test repr(m.addOne(x)) == repr("3.5")

        # This is a function that is in `m` from dex's prelude
        @test repr(m.sum(y)) == repr("9")
    end

    @testset "convert Atom" begin
        atom = convert(Atom, 1f0)
        @test convert(Number, atom) === 1f0
        @test convert(Real, atom) === 1f0
        @test convert(Float64, atom) === 1.0

        atom = convert(Atom, Int32(2))
        @test convert(Number, atom) === Int32(2)
        @test convert(Real, atom) === Int32(2)
        @test convert(Float64, atom) === 2.0
    end

    @testset "DexModule" begin
        m = DexModule("""
        x = 2.5
        y = [2, 3, 4]
        """)
        @test repr(m.x) == repr("2.5")
        @test repr(m.y) == repr("[2, 3, 4]")
    end
end

max_stars_repo_name: justinchiu/dex-lang

max_issues_repo_name: justinchiu/dex-lang

max_forks_repo_name: justinchiu/dex-lang


text:
using DiffEqBase, DiffEqJump, DiffEqBiological
using DiffEqJumpExtensions
using BenchmarkTools
using BioSimulator

# number of save times to look at
num_save_vec = 25000 #10.^(3:6)

# diffusion model
N = 256
diffnetwork = "@reaction_network difftype begin\n"
for i in 1:(N-1)
    diffnetwork *= "\t K, X$(i) --> X$(i+1)\n"
    diffnetwork *= "\t K, X$(i+1) --> X$(i)\n"
end
diffnetwork *= "end K"
#println(diffnetwork)
rs = eval( parse(diffnetwork) )
tf = 200.
params = (1.,)
prob   = DiscreteProblem(10*ones(Int64,N), (0.0, tf), params)

# Biosimlator model of the same network
model = Network("Diffusion")
for i in 1:N
    model <= Species("X$(i)", 10)
end
for i in 1:(N-1)
    model <= Reaction("diffusion$(i)a", params[1], "X$(i) --> X$(i+1)")
    model <= Reaction("diffusion$(i)b", params[1], "X$(i+1) --> X$(i)")
end

# mass action version
scaled_rates   = ones(Float64,2*(N-1))
reactant_stoch = Vector{Vector{Pair{Int64,Int64}}}();
net_stoch      = Vector{Vector{Pair{Int64,Int64}}}();
for i = 1:(N-1)
    push!(reactant_stoch, [i => 1])
    push!(reactant_stoch, [(i+1) => 1])
    push!(net_stoch, [i => -1, (i+1) => 1])
    push!(net_stoch, [(i+1) => -1, i => 1])
end
# reactant_stoch = Vector{Array{Int64,2}}();
# net_stoch      = Vector{Array{Int64,2}}();
# for i = 1:(N-1)
#     push!(reactant_stoch, [i,1]')
#     push!(reactant_stoch, [(i+1),1]')
#     push!(net_stoch, [i (i+1); -1 1])
#     push!(net_stoch, [(i+1) i; -1  1])
# end

# DiffEqJump methods to benchmark
#methods = (Direct2(), DirectVEC(), DiffEqJumpExtensions.FRM(), FRMVEC())
#methods = (DirectVEC(), FRMVEC())
methods = (DirectMA(scaled_rates, reactant_stoch, net_stoch), DirectVEC())

# exact methods
# println("Full solution paths saving:\n")
# for method in methods
#     println("Solving with method: ", typeof(method), ", using SSAStepper")
#     jump_prob = JumpProblem(prob, method, rs)
#     @btime solve($jump_prob, SSAStepper())
# end
# println()

# methods with different numbers of points
println("Saving at fixed times:\n")
for num_save in num_save_vec
    dtsave = tf / num_save

    for method in methods
        println("Solving with method: ", typeof(method), ", using SSAStepper and ", num_save, " points")
        jump_prob = JumpProblem(prob, method, rs, save_positions=(false,false))
        sol = solve(jump_prob, SSAStepper(), saveat=dtsave);
        @btime solve($jump_prob, SSAStepper(), saveat=$dtsave)
    end

    # BioSimulator:
    println("BioSimulator, Direct, saving at ", num_save, " points:")
    result = simulate(model, BioSimulator.SSA, time=tf, epochs=num_save, trials=1);
    @btime simulate($model, BioSimulator.SSA, time=$tf, epochs=$num_save, trials=1);

    # println("BioSimulator, FRM, saving at ", num_save, " points:")
    # @btime simulate($model, BioSimulator.FRM, time=$tf, epochs=$num_save, trials=1);

    println()
end

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
using Test

@testset "atexit.jl" begin
    function _atexit_tests_gen_cmd_eval(expr::String)
        cmd_eval = ```
        $(Base.julia_cmd()) -e $(expr)
        ```
        return cmd_eval
    end
    function _atexit_tests_gen_cmd_script(temp_dir::String, expr::String)
        script, io = mktemp(temp_dir)
        println(io, expr)
        close(io)
        cmd_script = ```
        $(Base.julia_cmd()) $(script)
        ```
        return cmd_script
    end
    atexit_temp_dir = mktempdir()
    atexit(() -> rm(atexit_temp_dir; force = true, recursive = true))
    @testset "these should exit with exit code 0" begin
        julia_expr_list = Dict(
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(0))
            exit(22)
            """ => 0,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            )
        for julia_expr in keys(julia_expr_list)
            cmd_eval = _atexit_tests_gen_cmd_eval(julia_expr)
            cmd_script = _atexit_tests_gen_cmd_script(atexit_temp_dir, julia_expr)
            expected_exit_code = julia_expr_list[julia_expr]
            @test success(cmd_eval)
            @test success(cmd_script)
            p_eval = run(cmd_eval; wait = false)
            p_script = run(cmd_script; wait = false)
            wait(p_eval)
            wait(p_script)
            @test p_eval.exitcode == expected_exit_code
            @test p_script.exitcode == expected_exit_code
        end
    end
    @testset "these should exit with exit code 1" begin
        julia_expr_list = Dict(
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(1))
            exit(22)
            """ => 1,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> println("No error"))
            atexit(() -> exit(1))
            exit(22)
            """ => 1,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(1))
            atexit(() -> exit(1))
            exit(22)
            """ => 1,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            )
        for julia_expr in keys(julia_expr_list)
            cmd_eval = _atexit_tests_gen_cmd_eval(julia_expr)
            cmd_script = _atexit_tests_gen_cmd_script(atexit_temp_dir, julia_expr)
            expected_exit_code = julia_expr_list[julia_expr]
            @test_throws ProcessFailedException run(cmd_eval)
            @test_throws ProcessFailedException run(cmd_script)
            p_eval = run(cmd_eval; wait = false)
            p_script = run(cmd_script; wait = false)
            wait(p_eval)
            wait(p_script)
            @test p_eval.exitcode == expected_exit_code
            @test p_script.exitcode == expected_exit_code
        end
    end
    @testset "test exit codes other than 0 or 1" begin
        julia_expr_list = Dict(
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(13))
            exit(22)
            """ => 13,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> println("No error"))
            atexit(() -> exit(5))
            exit(22)
            """ => 5,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(33))
            atexit(() -> exit(33))
            exit(22)
            """ => 33,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(21))
            atexit(() -> exit(21))
            exit(22)
            """ => 21,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            )
        for julia_expr in keys(julia_expr_list)
            cmd_eval = _atexit_tests_gen_cmd_eval(julia_expr)
            cmd_script = _atexit_tests_gen_cmd_script(atexit_temp_dir, julia_expr)
            expected_exit_code = julia_expr_list[julia_expr]
            @test_throws(ProcessFailedException, run(cmd_eval))
            @test_throws(ProcessFailedException, run(cmd_script))
            p_eval = run(cmd_eval; wait = false)
            p_script = run(cmd_script; wait = false)
            wait(p_eval)
            wait(p_script)
            @test p_eval.exitcode == expected_exit_code
            @test p_script.exitcode == expected_exit_code
        end
    end
    @testset "test what happens if multiple places call exit(n)" begin
        julia_expr_list = Dict(
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(22))
            atexit(() -> exit(11))
            atexit(() -> exit(33))
            """ => 22,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            """
            atexit(() -> exit(4))
            atexit(() -> exit(16))
            atexit(() -> exit(7))
            exit(22)
            """ => 4,
            # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            )
        for julia_expr in keys(julia_expr_list)
            cmd_eval = _atexit_tests_gen_cmd_eval(julia_expr)
            cmd_script = _atexit_tests_gen_cmd_script(atexit_temp_dir, julia_expr)
            expected_exit_code = julia_expr_list[julia_expr]
            @test_throws(ProcessFailedException, run(cmd_eval))
            @test_throws(ProcessFailedException, run(cmd_script))
            p_eval = run(cmd_eval; wait = false)
            p_script = run(cmd_script; wait = false)
            wait(p_eval)
            wait(p_script)
            @test p_eval.exitcode == expected_exit_code
            @test p_script.exitcode == expected_exit_code
        end
    end
    rm(atexit_temp_dir; force = true, recursive = true)
end

max_stars_repo_name: rezarastak/julia

max_issues_repo_name: rezarastak/julia

max_forks_repo_name: rezarastak/julia


text:
using TimingsInference
using Base.Test

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

t1 = GeoTrip(0.0f0, 1.0f0, 2.0f0, 3.0f0, DateTime(2016, 2, 9, 12, 35, 00), 60.0f0)
t2 = GeoTrip(0.0f0, 1.0f0, 2.0f0, 3.0f0, DateTime(2016, 2, 7, 13, 00, 00), 60.0f0)

@test inTimeWindow(t1, 12, 13) == true
@test inTimeWindow(t2, 12, 13) == false
@test onlyWeekdays(t1) == true
@test onlyWeekdays(t2) == false


text:
module MyExample

# Write your package code here.
greet() = print("Hello Torres!")

include("extra_file.jl")
export my_f

end #module


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

function read_mat(filename::String; model_id = nothing)
    mat_dict = MAT.matread(filename)
    if isnothing(model_id)
        model_id, model_mat = first(mat_dict)
    else
        model_mat = mat_dict[model_id]
    end
    return MetNet(model_mat; reshape = true)
end


text:
using LinearAlgebra

"""
    A = getcalmatrix(F::Array, H::Array, s::Vector)

In:
* `F` [nvoxel nshim] Acquired field maps
* `H` [nvoxel nbasis] Spatial basis functions
* `s` [nshim] applied shim currents (e.g., pairwise differences) used to acquire `F`
(excluding center frequency)
Out:
* `A` [nbasis nshim+1] (first column corresponds to center frequency)
"""
function getcalmatrix(F::AbstractMatrix{<:Real}, H::AbstractMatrix{<:Real}, s::AbstractVector)

	nb = size(H,2)      # number of basis functions. Includes center frequency.
	ns = length(s)      # number of shim terms (excluding center frequency)

	A = zeros(nb, ns+1)

	# center frequency (identity by definition)
	A[1,1] = 1.0

	# fit each acquired fieldmap to basis
	for ii = 1:ns
		A[:,ii+1] = s[ii] * H \ F[:,ii]
	end

	# A = inv(H'*H)*H'*F*inv(S);   [nb nb] Includes cf offset term.

	return A
end


text:
#=
# Stopwords list from sources
# * https://www.ranks.nl/stopwords
=#

default_english_stopwords = Set([
  "a", "able", "about", "above", "abst", "accordance", "according", "accordingly",
  "across", "act", "actually", "added", "adj", "affected", "affecting", "affects",
  "after", "afterwards", "again", "against", "ah", "all", "almost", "alone", "along",
  "already", "also", "although", "always", "am", "among", "amongst", "an", "and",
  "announce", "another", "any", "anybody", "anyhow", "anymore", "anyone", "anything",
  "anyway", "anyways", "anywhere", "apparently", "approximately", "are", "aren",
  "arent", "arise", "around", "as", "aside", "ask", "asking", "at", "auth",
  "available", "away", "awfully",
  "b", "back", "be", "became", "because", "become", "becomes", "becoming", "been",
  "before", "beforehand", "begin", "beginning", "beginnings", "begins", "behind",
  "being", "believe", "below", "beside", "besides", "between", "beyond", "biol",
  "both", "brief", "briefly", "but", "by",
  "c", "ca", "came", "can", "cannot", "can't", "cause", "causes", "certain",
  "certainly", "co", "com", "come", "comes", "contain", "containing", "contains",
  "could", "couldnt",
  "d", "date", "did", "didn't", "different", "do", "does", "doesn't", "doing",
  "done", "don't", "down", "downwards", "due", "during",
  "e", "each", "ed", "edu", "effect", "eg", "eight", "eighty", "either", "else",
  "elsewhere", "end", "ending", "enough", "especially", "et", "et-al", "etc",
  "even", "ever", "every", "everybody", "everyone", "everything", "everywhere",
  "ex", "except",
  "f", "far", "few", "ff", "fifth", "first", "five", "fix", "followed",
  "following", "follows", "for", "former", "formerly", "forth", "found",
  "four", "from", "further", "furthermore",
  "g", "gave", "get", "gets", "getting", "give", "given", "gives", "giving",
  "go", "goes", "gone", "got", "gotten",
  "h", "had", "happens", "hardly", "has", "hasn't", "have", "haven't",
  "having", "he", "hed", "hence", "her", "here", "hereafter", "hereby",
  "herein", "heres", "hereupon", "hers", "herself", "hes", "hi", "hid",
  "him", "himself", "his", "hither", "home", "how", "howbeit", "however", "hundred",
  "i", "id", "ie", "if", "i'll", "im", "immediate", "immediately", "importance",
  "important", "in", "inc", "indeed", "index", "information", "instead", "into",
  "invention", "inward", "is", "isn't", "it", "itd", "it'll", "its", "itself", "i've",
  "j", "just",
  "k", "keep","keeps", "kept", "kg", "km", "know", "known", "knows",
  "l", "largely", "last", "lately", "later", "latter", "latterly", "least",
  "less", "lest", "let", "lets", "like", "liked", "likely", "line", "little",
  "'ll", "look", "looking", "looks", "ltd",
  "m", "made", "mainly", "make", "makes", "many", "may", "maybe", "me",
  "mean", "means", "meantime", "meanwhile", "merely", "mg", "might", "million",
  "miss", "ml", "more", "moreover", "most", "mostly", "mr", "mrs", "much",
  "mug", "must", "my", "myself",
  "n", "na", "name", "namely", "nay", "nd", "near", "nearly", "necessarily",
  "necessary", "need", "needs", "neither", "never", "nevertheless", "new",
  "next", "nine", "ninety", "no", "nobody", "non", "none", "nonetheless",
  "noone", "nor", "normally", "nos", "not", "noted", "nothing", "now", "nowhere",
  "o", "obtain", "obtained", "obviously", "of", "off", "often", "oh", "ok",
  "okay", "old", "omitted", "on", "once", "one", "ones", "only", "onto",
  "or", "ord", "other", "others", "otherwise", "ought", "our", "ours", "ourselves",
  "out", "outside", "over", "overall", "owing", "own",
  "p", "page", "pages", "part", "particular", "particularly", "past", "per",
  "perhaps", "placed", "please", "plus", "poorly", "possible", "possibly",
  "potentially", "pp", "predominantly", "present", "previously", "primarily",
  "probably", "promptly", "proud", "provides", "put",
  "q", "que", "quickly", "quite", "qv", "r", "ran", "rather", "rd", "re",
  "readily", "really", "recent", "recently", "ref", "refs", "regarding",
  "regardless", "regards", "related", "relatively", "research", "respectively",
  "resulted", "resulting", "results", "right", "run",
  "s", "said", "same", "saw", "say", "saying", "says", "sec", "section", "see",
  "seeing", "seem", "seemed", "seeming", "seems", "seen", "self", "selves",
  "sent", "seven", "several", "shall", "she", "shed", "she'll", "shes",
  "should", "shouldn't", "show", "showed", "shown", "showns", "shows",
  "significant", "significantly", "similar", "similarly", "since", "six",
  "slightly", "so", "some", "somebody", "somehow", "someone", "somethan",
  "something", "sometime", "sometimes", "somewhat", "somewhere", "soon",
  "sorry", "specifically", "specified", "specify", "specifying", "still",
  "stop", "strongly", "sub", "substantially", "successfully", "such",
  "sufficiently", "suggest", "sup", "sure",
  "t", "take", "taken", "taking", "tell", "tends", "th", "than", "thank",
  "thanks", "thanx", "that", "that'll", "thats", "that've", "the", "their",
  "theirs", "them", "themselves", "then", "thence", "there", "thereafter",
  "thereby", "thered", "therefore", "therein", "there'll", "thereof", "therere",
  "theres", "thereto", "thereupon", "there've", "these", "they", "theyd",
  "they'll", "theyre", "they've", "think", "this", "those", "thou", "though",
  "thoughh", "thousand", "throug", "through", "throughout", "thru", "thus",
  "til", "tip", "to", "together", "too", "took", "toward", "towards", "tried",
  "tries", "truly", "try", "trying", "ts", "twice", "two",
  "u", "un", "under", "unfortunately", "unless", "unlike", "unlikely",
  "until", "unto", "up", "upon", "ups", "us", "use", "used", "useful",
  "usefully", "usefulness", "uses", "using", "usually",
  "v", "value", "various", "'ve", "very", "via", "viz", "vol", "vols", "vs",
  "w", "want", "wants", "was", "wasnt", "way", "we", "wed", "welcome",
  "we'll", "went", "were", "werent", "we've", "what", "whatever", "what'll",
  "whats", "when", "whence", "whenever", "where", "whereafter", "whereas",
  "whereby", "wherein", "wheres", "whereupon", "wherever", "whether",
  "which", "while", "whim", "whither", "who", "whod", "whoever", "whole",
  "who'll", "whom", "whomever", "whos", "whose", "why", "widely", "willing",
  "wish", "with", "within", "without", "wont", "words", "world", "would",
  "wouldnt", "www",
  "x",
  "y", "yes", "yet", "you", "youd", "you'll", "your", "youre",
  "yours", "yourself", "yourselves", "you've",
  "z", "zero"
])


text:
include("trackingcell.jl");
include("lineage.jl");
include("tiffxml.jl");
include("segmentation3d.jl")

data_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF";
ret_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_hamamatsu-60X-TIRF/20200316_result";

paths_all = []
for s in 1
	img = load(File(format"TIFF","$data_dir/20200316/HE7-11-1-80uw-PWM_1_s$s.ome.tiff"));
	img_edge, mask_edge, mask_clear = split_cell_LoG(img, size(img)[3]÷20);
	save("$ret_dir/img_edge_s$s.tiff", img_edge);

	nucleus_pos = [ component_centroids(label_components(mask_clear[:,:,i])) for i in 1:size(mask_clear)[3] ];

	paths = generate_path(nucleus_pos);
	push!(paths_all, paths)

	for i in 1:length(paths)
		save("$ret_dir/s$s-c$i.tiff", N0f16.(export_cell(paths[i], img, mask_edge)));
		embedxml(700, 700, 20, length(paths[i]), "$ret_dir/s$s-c$i.tiff");
	end
end


text:
"""
    EasyPlot

Masked APIs of plotting; based on PyPlot
"""
module EasyPlot
   import PyPlot  # based on PyPlot

# -----------------------------------------------------------------------------
# a public Dict of the attributes of fonts; modify a copy of it in functions
public_fontdict = Dict(
   "family" => "Times New Roman",
   "weight" => "normal",
   "color" => "black",
   "size" => 16
)

# ----------
"""
    Plot_SteadyState( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ; outpdf::Union{Nothing,String} = nothing )

plot 4 sub-plots: capital distribution, consumption distribution, labor distribution, asset & UEBMI-Indi distribution;
the keyword parameter outpdf= indicates whether to output to a pdf figure & the file directory,
it is nothing ("not output") by default;
"""
function Plot_SteadyState( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ;
      outpdf::Union{Nothing,String} = nothing, picsize::Tuple = (12,8) )
   # define an age-based x-axis
   local Xval = env.START_AGE : env.START_AGE + env.S - 1
   # plotting
   PyPlot.figure( figsize = picsize )
   PyPlot.subplot(2,2,1)
      PyPlot.plot(Xval, Dst[:𝒜][t,:] .* Ps[:N][t,:] )
      PyPlot.xlabel("Age"); PyPlot.ylabel("k Distribution among ages"); PyPlot.grid(true)
   PyPlot.subplot(2,2,2)
      PyPlot.plot(Xval, Dst[:c][t,:] .* Ps[:N][t,:] )
      PyPlot.xlabel("Age"); PyPlot.ylabel("Consumption Distribution"); PyPlot.grid(true)
   PyPlot.subplot(2,2,3)
      PyPlot.plot(Xval[1:env.Sr], Dst[:Lab][t,:] .* Ps[:N][t,1:env.Sr] )
      PyPlot.xlabel("Age"); PyPlot.ylabel("Labor Distribution"); PyPlot.grid(true)
   PyPlot.subplot(2,2,4)
      PyPlot.plot(Xval, Dst[:a][t,:] .* Ps[:N][t,:] )
      PyPlot.plot(Xval, Dst[:Φ][t,:] .* Ps[:N][t,:] )
      PyPlot.legend(["Asset Distri","UEBMI-Indi Distri"]); PyPlot.grid(true)
   PyPlot.tight_layout()

   if isa(outpdf, String)
   PyPlot.savefig( outpdf, format = "pdf" )
   end

   return nothing
end
# ----------
"""
   Plot_PerformProfile( PerfLog::Dict ; startiter::Int = 1 )

Plots for the Dict of performance/convergence returned by EasySearch.Transition!(ReturnLog=true);
startiter= receives an integer indicating which round to start to plot;
returns nothing
"""
function Plot_PerformProfile( PerfLog::Dict ; startiter::Int = 1 )
   local XVAL = startiter:length(PerfLog[:K])
   PyPlot.figure( figsize = (8,4) )
   PyPlot.subplot(1,2,1)
      PyPlot.plot( XVAL, PerfLog[:K][XVAL], linewidth = 1.0 )
      PyPlot.title("Capital (K)")
      PyPlot.xlabel("Iteration Round")
      PyPlot.ylabel("Relative Error (%)")
      PyPlot.grid(true)
   PyPlot.subplot(1,2,2)
      PyPlot.plot( XVAL, PerfLog[:L][XVAL], linewidth = 1.0 )
      PyPlot.title("Labor (L)")
      PyPlot.xlabel("Iteration Round")
      PyPlot.grid(true)
   PyPlot.suptitle("Convergence Profile")
   # nominal returns
   return nothing
end
# ----------
"""
   Plot_Transition( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ; outpdf::Union{Nothing,String} = nothing, picsize::Tuple = (12,8) )

the keyword parameter outpdf= indicates whether to output to a pdf figure & the file directory,
it is nothing ("not output") by default;

Plots:
1. Y,   ΔY/Y,   LI/Y,   Δ(LI/Y)/(LI/Y)
2. r,   w̄,      K,      L
3. C,   I,      G,      SubstituionRate(Pension)
"""
function Plot_Transition( Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ;
   YearRange::Tuple = ( env.START_YEAR, env.START_YEAR + env.T - 1 ),  # the range of years to plot
   LineWidth::Float64 = 1.0,  # the width of lines to plot
   outpdf::Union{Nothing,String} = nothing, picsize::Tuple = (12,8) )
   # -----------------
   # assertion: years
   @assert( env.START_YEAR <= YearRange[1] < YearRange[2] < (env.START_YEAR + env.T - 1)  ,  "out of the range of years! Maximum T-1 years" )
   # prepare the values of X axis
   local tmpDisplaceYear = 1 - env.START_YEAR  # the shift to convert real years to index
   local XYEAR = YearRange[1]:YearRange[2]  # the years to plot on x-axis
   local XLOC = XYEAR .+ tmpDisplaceYear  # the index to slice data

   # prepare the growth rate of LI/Y
   local LI2Ygrowth = (Dt[:LI2Y][2:env.T] ./ Dt[:LI2Y][1:env.T-1] .- 1.0) .* 100

   # maximum layout
   local tmpLayout = (3,4)
   # Plotting
   PyPlot.figure( figsize = picsize )
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 1)
      PyPlot.plot( XYEAR, Dt[:Y][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("GDP")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 2)
      PyPlot.plot( XYEAR, Dt[:GDPgrowth][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("GDP Growth Rate (%)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 3)
      PyPlot.plot( XYEAR, Dt[:LI2Y][XLOC] .* 100, linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("UEBMI-POOL Gaps / GDP (%)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 4)
      PyPlot.plot( XYEAR, LI2Ygrowth[XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Growth of UEBMI-POOL Gaps / GDP (%)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 5)
      PyPlot.plot( XYEAR, Dt[:r][XLOC] .* 100, linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Net Interest Rate (%)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 6)
      PyPlot.plot( XYEAR, Dt[:w̄][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Average Wage Level")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 7)
      PyPlot.plot( XYEAR, Dt[:K][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Aggregated Capital (K)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 8)
      PyPlot.plot( XYEAR, Dt[:L][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Aggregated Labor (L)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 9)
      PyPlot.plot( XYEAR, Dt[:C][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Aggregated Consumption (C)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 10)
      PyPlot.plot( XYEAR, Dt[:I][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Aggregated Investment (I)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 11)
      PyPlot.plot( XYEAR, Dt[:G][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Government Purchase (G)")
      # ------------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 12)
      PyPlot.plot( XYEAR, Dt[:SubstiRat][XLOC] .* 100, linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Substitution Rate of Pension (%)")

      # tight layout
      PyPlot.tight_layout()

      # output
      if isa(outpdf, String)
      PyPlot.savefig( outpdf, format = "pdf" )
      end

   return nothing
end
# --------------
"""
   Plot_Calibrate( Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ;
      YearRange::Tuple = ( env.START_YEAR, env.START_YEAR + env.T - 1 ),  # the range of years to plot
      LineWidth::Float64 = 1.0,  # the width of lines to plot
      outpdf::Union{Nothing,String} = nothing, picsize::Tuple = (12,8) )

plots key indicators in calibration:
   1. LI/GDP
   2. LI/pool account expenditure
   3. LI/pool account incomes
"""
function Plot_Calibrate( Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ;
   YearRange::Tuple = ( env.START_YEAR, env.START_YEAR + env.T - 1 ),  # the range of years to plot
   LineWidth::Float64 = 1.0,  # the width of lines to plot
   outpdf::Union{Nothing,String} = nothing, picsize::Tuple = (12,8), tmpLayout::Tuple = (1,2),
   GlobalFontSize::Int = 12 )
   #  ---------------
   # assertion: years
   @assert( env.START_YEAR <= YearRange[1] < YearRange[2] < (env.START_YEAR + env.T - 1)  ,  "out of the range of years! Maximum T-1 years" )
   # prepare the values of X axis
   local tmpDisplaceYear = 1 - env.START_YEAR  # the shift to convert real years to index
   local XYEAR = YearRange[1]:YearRange[2]  # the years to plot on x-axis
   local XLOC = XYEAR .+ tmpDisplaceYear  # the index to slice data

   # prepare the growth rate of LI/Y
   local LI2Ygrowth = (Dt[:LI2Y][2:env.T] ./ Dt[:LI2Y][1:env.T-1] .- 1.0) .* 100

   # maximum layout
   # local tmpLayout = (1,2)
   # Plotting
   PyPlot.figure( figsize = picsize )
      # PyPlot.subplot(tmpLayout[1],tmpLayout[2], 1)
      # PyPlot.plot( XYEAR, 100.0 .* Dt[:LI][XLOC] ./ Dt[:Y][XLOC], linewidth = LineWidth )
      #    PyPlot.xlabel("Year"); PyPlot.grid(true)
      #    PyPlot.ylabel("LI/Y (%)")
      # -----
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 1)
      PyPlot.plot( XYEAR, 100.0 .* Dt[:LI][XLOC] ./ Dt[:AggPoolExp][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Percentage (%)")
      # ------
      PyPlot.subplot(tmpLayout[1],tmpLayout[2], 2)
      PyPlot.plot( XYEAR, 100.0 .* Dt[:LI][XLOC] ./ Dt[:AggPoolIn][XLOC], linewidth = LineWidth )
         PyPlot.xlabel("Year"); PyPlot.grid(true)
         PyPlot.ylabel("Percentage (%)")

      # tight layout
      PyPlot.tight_layout()

      # output
      if isa(outpdf, String)
      PyPlot.savefig( outpdf, format = "pdf" )
      end

end  # function ends
#  -----------

end  # module ends
#


text:
module Histogram

# simple histogram functions copied from Julia v0.4.7 to avoid having to
# pull in all of Statistics.jl

export hist, hist!

function histrange(v::AbstractArray{T}, n::Integer) where {T<:AbstractFloat}
    nv = length(v)
    if nv == 0 && n < 0
        throw(ArgumentError("number of bins must be ≥ 0 for an empty array, got $n"))
    elseif nv > 0 && n < 1
        throw(ArgumentError("number of bins must be ≥ 1 for a non-empty array, got $n"))
    end
    if nv == 0
        return 0.0:1.0:0.0
    end
    lo, hi = extrema(v)
    if hi == lo
        step = 1.0
    else
        bw = (hi - lo) / n
        e = 10.0^floor(log10(bw))
        r = bw / e
        if r <= 2
            step = 2*e
        elseif r <= 5
            step = 5*e
        else
            step = 10*e
        end
    end
    start = step*(ceil(lo/step)-1)
    nm1 = ceil(Int,(hi - start)/step)
    return start:step:(start + nm1*step)
end

function histrange(v::AbstractArray{T}, n::Integer) where {T<:Integer}
    nv = length(v)
    if nv == 0 && n < 0
        throw(ArgumentError("number of bins must be ≥ 0 for an empty array, got $n"))
    elseif nv > 0 && n < 1
        throw(ArgumentError("number of bins must be ≥ 1 for a non-empty array, got $n"))
    end
    if nv == 0
        return 0:1:0
    end
    if n <= 0
        throw(ArgumentError("number of bins n = $n must be positive"))
    end
    lo, hi = extrema(v)
    if hi == lo
        step = 1
    else
        bw = (hi - lo) / n
        e = 10^max(0,floor(Int,log10(bw)))
        r = bw / e
        if r <= 1
            step = e
        elseif r <= 2
            step = 2*e
        elseif r <= 5
            step = 5*e
        else
            step = 10*e
        end
    end
    start = step*(ceil(lo/step)-1)
    nm1 = ceil(Int,(hi - start)/step)
    return start:step:(start + nm1*step)
end

# Sturges' formula
function sturges(n::Integer)
    n==0 && return one(n)
    return ceil(Int, log2(n)) + 1
end

function hist!(h::AbstractArray{HT}, v::AbstractVector, edg::AbstractVector) where HT
    n = length(edg) - 1
    length(h) == n || throw(DimensionMismatch("length(histogram) must equal length(edges) - 1"))
    fill!(h, zero(HT))
    for x in v
        i = searchsortedfirst(edg, x)-1
        if 1 <= i <= n
            h[i] += 1
        end
    end
    return edg, h
end

hist(v::AbstractVector, edg::AbstractVector) = hist!(Array{Int}(undef, length(edg)-1), v, edg)
hist(v::AbstractVector, n::Integer) = hist(v,histrange(v,n))
hist(v::AbstractVector) = hist(v,sturges(length(v)))

end


text:
# # Parsing MATPOWER Files

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# An example of how to parse MATPOWER files and create a `System` using [PowerSystems.jl](github.com/NREL-SIIP/PowerSystems.jl)

# ### Environment
# This notebook depends on the SIIPExamples.jl environment which is loaded by default

using Pkg
Pkg.status()

# ### Dependencies
using SIIPExamples
using PowerSystems
using TimeSeries

# ### Fetch Data
# PowerSystems.jl links to some test data that is suitable for this example.
# Let's download the test data
base_dir = PowerSystems.download(PowerSystems.TestData; branch = "master");

# ### Create a `System`
sys = System(PowerModelsData(joinpath(base_dir, "matpower", "case5_re.m")))
sys


text:
include("hyperelasticity.jl")
#
ts = ExplicitSolverTime(Δt, NT)
ubd, abd = compute_boundary_info(domain, globaldata, ts)
Fext = compute_external_force(domain, globaldata, ts)

d0 = zeros(2domain.nnodes)
v0 = zeros(2domain.nnodes)
a0 = zeros(2domain.nnodes)

mode = "consistent_tangent"
if length(ARGS)>=1
  global mode = ARGS[1]
end

if mode=="linear"
  global Dc, nn_law
  Dc = Variable(rand(3,3))
  Dc = spd(Dc)
  function nn_law(strain)
    strain_tensor = strain_voigt_to_tensor(strain)
    stress = strain*Dc
    stress
  end

elseif mode=="consistent_tangent"
  global Dc, nn_law
  Dc = Variable(rand(3,3))
  Dc = spd(Dc)
  function nn_law(strain)
    coef = ae(strain, [20,20,20,6])
    coef = [coef constant(ones(size(coef,1), 1))]
    H = consistent_tangent_matrix(coef, Dc)
    stress = batch_matmul(H, strain)
    stress
  end

elseif mode=="nn"
  global Dc, nn_law
  # general neural network
  function nn_law(ε)
    ae(ε, [20,20,20,3])
  end
elseif mode=="free_energy"
  global Dc, nn_law
  # free energy
  function nn_law(ε)
    φ = squeeze(ae(ε, [20,20,20,1]))
    tf.gradients(φ, ε)[1]
  end
else
  error("$mode not valid")
end

d, v, a= ExplicitSolver(globaldata, domain, d0, v0, a0, Δt, NT, nn_law, Fext, ubd, abd; strain_type="finite")

# test  = open("test.txt", "w")
sess = Session(); init(sess)
# run(sess, d)
idx = @. (0:n)*m + 1
idx = [idx; (@. idx + domain.nnodes)]
loss = sum((d[:, idx] - d_[:, idx])^2)
sess = Session(); init(sess)
@info run(sess, loss)
# error()

!isdir(mode) && mkdir(mode)

for i = 1:100
  loss_ = BFGS!(sess, loss, 1000)
  d0 = run(sess, d)

  # visualize
  close("all")
  p = visualize_displacement(d0, domain)
  saveanim(p, "$mode/$i.gif")
  close("all")
  plot(d_[:,1], "-", color="C1")
  plot(d0[:,1], "--", color="C1")
  plot(d_[:,1+domain.nnodes], color="C2")
  plot(d0[:,1+domain.nnodes], "--", color="C2")
  savefig("$mode/$i.png")

  open("$mode/loss.txt", "a") do io
    writedlm(io, loss_)
  end
  ADCME.save(sess, "$mode/$i.mat")
  if length(loss_)<1000
    break
  end
end


text:
export binary_search

function binary_search(array::Vector{T}, target::T) where T
    left = 1
    right = length(array)

    while left <= right
        mid = trunc(Int, left+(right-left)/2)
        if array[mid] > target
            right = mid-1
        elseif array[mid] < target
            left = mid+1
        else
            return mid
        end
    end

    return -1
end


text:
using HDF5
using Test

@testset "mmap" begin

# Create a new file
fn = tempname()
f = h5open(fn, "w")
@test isopen(f)

# Create two datasets, one with late allocation (the default for contiguous
# datasets) and the other with explicit early allocation.
hdf5_A = d_create(f, "A", datatype(Int64), dataspace(3,3))
hdf5_B = d_create(f, "B", datatype(Float64), dataspace(3,3);
                  alloc_time = HDF5.H5D_ALLOC_TIME_EARLY)
# The late case cannot be mapped yet.
@test_throws ErrorException("Error getting offset") readmmap(f["A"])
# Then write and fill dataset A, making it mappable. B was filled with 0.0 at
# creation.
A = rand(Int64,3,3)
hdf5_A[:,:] = A
flush(f)
close(f)
# Read HDF5 file & MMAP
f = h5open(fn,"r")
A_mmaped = readmmap(f["A"])
@test all(A .== A_mmaped)
@test all(iszero, readmmap(f["B"]))
# Check that it is read only
@test_throws ReadOnlyMemoryError A_mmaped[1,1] = 33
close(f)
# Now check if we can write
f = h5open(fn,"r+")
A_mmaped = readmmap(f["A"])
A_mmaped[1,1] = 33
close(f)

end