text:
__precompile__()
module CholmodSolve2

import Base.SparseArrays.CHOLMOD: SuiteSparse_long, Factor,
        C_Dense, C_Sparse, C_Factor, @cholmod_name, Cint, CHOLMOD_A, common,
        CHOLMODException, xtyp

const DEBUG = false
# Only allow Float64 for now
const Tv = Float64 # TODO Relax to VTypes

# Julia, dont CG these, or we will have problems!!
const CHOLMOD_ptrs = Vector{Ptr{C_Dense{Tv}}}(4)

# Needed for precompilation
function __init__()
    # Make sure we have 4 NULL pointers allocated that we can reference to
    map!(_-> Ptr{C_Dense{Tv}}(C_NULL), CHOLMOD_ptrs, 1:4)
    # Workspace pointers that CHOLMOD will use. CHOLMOD will take care of ALL allocations
    global const CHOLMOD_BHandle = pointer(CHOLMOD_ptrs,1)
    global const CHOLMOD_XHandle = pointer(CHOLMOD_ptrs,2)
    global const CHOLMOD_YHandle = pointer(CHOLMOD_ptrs,3)
    global const CHOLMOD_EHandle = pointer(CHOLMOD_ptrs,4)
end

function solve2!(sys::Integer, F::Factor{Tv}, B::Ptr{C_Dense{Tv}})
    if size(F,1) != unsafe_load(B).nrow
        throw(DimensionMismatch("B must have as many rows as F in A_ldiv_B!"))
    end

    cmn = common()
    d = ccall((@cholmod_name("solve2", SuiteSparse_long),:libcholmod), Cint,
            (Cint, Ptr{C_Factor{Tv}},       # sys, *L
            Ptr{C_Dense{Tv}},   Ptr{Void},  # *B  *Bset
            Ptr{Void},          Ptr{Void},  # **X **Xset
            Ptr{Void},          Ptr{Void},  # **Y **E
            Ref{UInt8}),                    # *common
            sys, get(F.p),
            B, C_NULL,                      # Bset = C_NULL
            CHOLMOD_XHandle, C_NULL,        # Xset = C_NULL
            CHOLMOD_YHandle, CHOLMOD_EHandle,
            cmn)
    if CHOLMOD_XHandle == C_NULL
        throw(CHOLMODException("Solve failed with unexpected error"))
    elseif unsafe_load(CHOLMOD_XHandle) == C_NULL
        throw(CHOLMODException("Solve failed with NULL output"))
    end
    return d
end

function Base.A_ldiv_B!(C::AbstractVecOrMat{Tv}, F::Factor{Tv}, B::AbstractVecOrMat{Tv}) where Tv
    if size(B,1) != size(F,1) || size(C,1) != size(F,1)
        throw(DimensionMismatch("B and C must have same number of rows as F, got size(F)=$(size(F)), size(B)=$(size(B)), size(C)=$(size(C))."))
    end
    if size(C,2) != size(B,2)
        throw(DimensionMismatch("B and C must have same number of columns, got size(B)=$(size(B)), size(C)=$(size(C))."))
    end
    # Allocate Dense B if needed
    DEBUG && println("Julia ensure B: : $(size(B,1)), $(size(B,2))")
    ensure_dense(CHOLMOD_BHandle, size(B,1), size(B,2), size(B,1))
    # Copy B to CHOLMOD input
    _copy!(unsafe_load(CHOLMOD_BHandle), B)
    #Try to resize Y to avoid unnessesary allocation in solve2
    if CHOLMOD_YHandle != C_NULL && unsafe_load(CHOLMOD_YHandle) != C_NULL
        Yp = unsafe_load(CHOLMOD_YHandle)
        Y = unsafe_load(Yp)
        # This is how many rows solve2 wants in Y
        newrows = max(4, size(B,2))
        # TODO We don't dare to change ncol yet
        if Y.ncol == size(B,1) && Y.nrow == Y.d && Y.nzmax >= Y.ncol*newrows
            # Do change nrow and d to what solve2 wants
            DEBUG && println("Julia ensure Y: $newrows, $(Y.ncol)")
            ensure_dense(CHOLMOD_YHandle, newrows, Y.ncol, newrows)
        end
    end
    # Only needed when dimensions of B actually change, but solve2 won't handle it nicely
    DEBUG && println("Julia ensure X: : $(size(B,1)), $(size(B,2))")
    ensure_dense(CHOLMOD_XHandle, size(B,1), size(B,2), size(B,1))

    status = solve2!(CHOLMOD_A, F, unsafe_load(CHOLMOD_BHandle))
    if status != 1
        throw(DimensionMismatch("Solve failed"))
    end
    # Copy the CHOLMOD output to C
    _copy!(C, unsafe_load(CHOLMOD_XHandle))

    return C
end

function ensure_dense(XHandle::Ptr{Ptr{C_Dense{Tv}}}, nrow, ncol, d)
    if XHandle != C_NULL
        XDensep = unsafe_load(XHandle)
        if XDensep != C_NULL
            XDense = unsafe_load(XDensep)
            # To be safe, only do if nrow=ncol
            if XDense.xtype == xtyp(Tv) && nrow == d && XDense.nzmax >= d*ncol
                #We can just change the dimensions, by owerwringing memory we are not allowed to write to...
                DEBUG && println("Julia changed dimensions")
                unsafe_store!(Base.convert(Ptr{Csize_t}, XDensep), nrow, 1) # nrow
                unsafe_store!(Base.convert(Ptr{Csize_t}, XDensep), ncol, 2)     # ncol
                unsafe_store!(Base.convert(Ptr{Csize_t}, XDensep), nrow, 4) # d
                return
            end
        end
    end
    DEBUG && println("Julia couldn't change dimensions")
    X = ccall((@cholmod_name("ensure_dense", SuiteSparse_long),:libcholmod), Ptr{C_Dense{Tv}},
           (Ptr{Ptr{C_Dense{Tv}}},
            Csize_t, Csize_t, Csize_t, Cint, Ptr{UInt8}),
            XHandle,
            nrow, ncol, d, xtyp(Tv), common())
    if XHandle == C_NULL
        throw(CHOLMODException("Ensure dense failed with NULL pointer, this should not happen"))
    elseif unsafe_load(XHandle) == C_NULL
        throw(CHOLMODException("Ensure dense failed in CHOLMOD"))
    end
    return
end

function _copy!(dest::AbstractArray{Tv}, D::Ptr{C_Dense{Tv}})
    s = unsafe_load(D)
    n = s.nrow*s.ncol
    n <= length(dest) || throw(BoundsError(dest, n))
    if s.d == s.nrow && (isa(dest, Array) || Base.iscontiguous(dest))
        unsafe_copy!(pointer(dest), s.x, s.d*s.ncol)
    else
        k = 0
        for j = 1:s.ncol
            for i = 1:s.nrow
                dest[k+=1] = unsafe_load(s.x, i + (j - 1)*s.d)
            end
        end
    end
    dest
end

# TODO verify correctness when s.d != s.nrow
function _copy!(D::Ptr{C_Dense{Tv}}, orig::AbstractVecOrMat{Tv})
    s = unsafe_load(D)
    n = s.nrow*s.ncol
    length(orig) <= n || throw(BoundsError(n, orig))
    if s.d == s.nrow && (isa(orig, Array) || Base.iscontiguous(orig))
        unsafe_copy!(s.x, pointer(orig), length(orig))
    else
        k = 0
        for j = 1:s.ncol
            for i = 1:s.nrow
                unsafe_store!(s.x, orig[k+=1], i + (j - 1)*s.d)
            end
        end
    end
    D
end

# For debugging memory use
function print_common()
    Base.SparseArrays.CHOLMOD.set_print_level(common(), 5)
    ccall((@cholmod_name("print_common", SuiteSparse_long),:libcholmod), Cint,
                (Ptr{UInt8}, Ref{UInt8}),                                    # *common
                    "Test1", common())
end

end # module


text:
abstract type Animal end
abstract type Dog <: Animal end
abstract type Cat <: Animal end

struct Lab <: Dog end
struct Collie <: Dog end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    get_interaction_matrix(::LinearOptics)
"""
function get_interaction_matrix(problem::LinearOptics)
    @debug "start: get_interaction_matrix"

    G = get_empty_matrix(problem.physic, problem.atoms)
    problem.kernelFunction(problem.atoms, problem.laser, G)

    @debug "end  : get_interaction_matrix"
    return G
end

function get_empty_matrix(physic::Scalar, atoms::Atom{<:ThreeD})
    Array{ComplexF64}(undef, atoms.N, atoms.N)
end
function get_empty_matrix(physic::Vectorial, atoms::Atom{<:ThreeD})
    Array{ComplexF64}(undef, 3atoms.N, 3atoms.N)
end

function get_interaction_matrix(problem::NonLinearOptics{MeanField})
    @debug "start: get_interaction_matrix"

    if haskey(problem.data, :G)
        return problem.data[:G]
    end

    temp_scalar_problem = LinearOptics(Scalar(), problem.atoms, problem.laser)
    G = get_empty_matrix(temp_scalar_problem.physic, temp_scalar_problem.atoms)
    temp_scalar_problem.kernelFunction(temp_scalar_problem.atoms, problem.laser, G)

    temp_scalar_problem = 1; GC.gc()
    problem.data[:G] = G
    @debug "end  : get_interaction_matrix"
    return G
end

"""
    green_scalar!(atoms, laser, G)

Computes:
    @. G = -(Γ/2)*exp(1im*k₀ * R_jk) / (1im * k₀ * R_jk)
    G[diagind(G)] .= 1im * laser.Δ - Γ/2
"""
function green_scalar!(atoms, laser, G)
    @debug "start: green_scalar!"

    G[:] = get_pairwise_matrix(atoms.r) # R_jk

    Threads.@threads for j in eachindex(G) #
        @inbounds G[j] = -(Γ / 2) * cis(k₀ * G[j]) / (1im * k₀ * G[j])
    end
    G[diagind(G)] .= 1im * laser.Δ - Γ / 2

    @debug "end  : green_scalar!"
    return nothing
end

"""
    experimental: green_vectorial!(atoms, laser, G)
"""
function green_vectorial!(atoms, laser, G)
    @debug "start: green_vectorial!"

    δ(x,y) = Int(==(x,y))

    Xt, Yt, Zt = r[:,1], r[:,2] ,r[:,3]

    Xjn = Xt*ones(1,N) - ones(N,1)*Xt'
    Yjn = Yt*ones(1,N) - ones(N,1)*Yt'
    Zjn = Zt*ones(1,N) - ones(N,1)*Zt'

    array_XYZ_jn = [Xjn, Yjn, Zjn]
    Rjn = sqrt.(Xjn.^2 + Zjn.^2 + Yjn.^2)

    P(x) =  1 - 1/x + 1/x^2
    Q(x) = -1 + 3/x - 3/x^2

    α_range = β_range = [-1, +1, 0]

    term2 = (3/2)*exp.(im*Rjn)./Rjn
    term2[findall(isnan.(term2))] .= 0
    term2[findall(isinf.(term2))] .= 0

    P_Rjn = P.(im*Rjn)
    P_Rjn[findall(isnan.(P_Rjn))] .= 0

    Q_Rjn_over_Rjn_squared = Q.(im*Rjn)./Rjn.^2
    Q_Rjn_over_Rjn_squared[findall(isnan.(Q_Rjn_over_Rjn_squared))] .= 0

    A = []
    for (α_idx, α) ∈ enumerate(α_range)
        B = []
        for (β_idx, β) ∈ enumerate(β_range)
            term1 = im*I(N).*δ(α,β)
            # term2 = (3/2)*exp.(im*Rjn)./Rjn  ## Defined outside
            term3 = P_Rjn.*δ(α,β)
            term4 = Q_Rjn_over_Rjn_squared.*(array_XYZ_jn[α_idx].*array_XYZ_jn[β_idx])

            K = term1 + term2.*(term3 .+ term4)
            push!(B, K)
        end
        push!(A, vcat(B[1:length(β_range)]...)  )
    end

    G[:] = hcat(A[1:length(α_range)]...)

    @debug "end  : green_vectorial!"
    return nothing
end


text:
@testset "Warm up point generation" begin
    model = load_model(StandardModel, model_paths["e_coli_core.json"])

    pts, lbs, ubs = warmup_from_variability(
        model,
        Tulip.Optimizer,
        100;
        modifications = [change_constraint("EX_glc__D_e"; lb = -2, ub = 2)],
        workers = W,
    )

    idx = first(indexin(["EX_glc__D_e"], reactions(model)))
    @test size(pts) == (95, 100)
    @test all(pts[idx, :] .>= -2)
    @test all(pts[idx, :] .<= 2)
    @test lbs[idx] == -2
    @test ubs[idx] == 2
end


text:
# This file is used for performing experiments
# Run by changing path in REPL to scripts folder
# Then `include("experiments.jl")`

using AutomotiveSimulator
using AutomotiveInteraction
using PGFPlots
using JLD

include("helpers.jl")

# Scenarios for upper merge
# List of tuples. Elem 1 of tuple is veh id list, elem 2 is list with ts and te
scenarios_upper = [
        ([6,8,13,19,28,29],[1,120]), #19,29 merge lane
        ([34,37,40,42,43,49,50],[88,230]), #34,37,43,50 merge lane
        ([59,60,62,66,67,71,73],[188,390]), #59,62,66,71 merge lane
        ([66,67,71,73,76,82,88,92,95],[310,440]), #71,76,95 merge lane
        ([187,191,193,194],[940,985]), #187,191,194 merge lane

        ([211,213,215,217,221],[1135,1190]), #211,217,221 merge lane
        ([248,249,250,256],[1323,1420]), # 249 merge lane
        ([272,275,276,277,279],[1440,1615]), #275, 277 merge lane
        ([294,295,300,307,308,311],[1576,1780]), #295,307,308 merge lane
        ([324,327,329,333,335,339,341,346,354,357],[1810,1950]) #327,333,339,346,354
]

# Scenarios for lower merge
scenarios_lower = [
        ([2,7,10,18,25],[1,30]), #7,25 merge vehicles
        ([64,68,70,74,77],[198,230]), #70 merge lane
        ([105,106,108,110],[382,430]), # 110 merge lane
        ([176,177,178],[827,870]), #177 merge lane

        ([252,253,255],[1315,1370]), #253 merge lane
        ([361,363,364,366],[1830,1885]) #363 merge lane
]

#***********Run filtering and store resulting models to jld files****
# s = scenarios_lower
# f = FilteringEnvironment(mergeenv=MergingEnvironmentLower())

s = scenarios_upper
f = FilteringEnvironment()
for i in 1:length(s)
        veh_id_list = s[i][1]
        ts = s[i][2][1]
        te = s[i][2][2]
        m,p,md = obtain_driver_models(f,veh_id_list=veh_id_list,num_p=50,ts=ts,te=te)
        #JLD.save("media/lower_$i.jld","m",m,"p",p,"md",md,"veh_id_list",veh_id_list,"ts",ts,"te",te)
        JLD.save("media/upper_$i.jld","m",m,"p",p,"md",md,"veh_id_list",veh_id_list,"ts",ts,"te",te)
end

#************Docstring example code from multiscenarios_pf in helpers.jl*****
# USEFUL to keep around for making tikz plots later on by rerunning this script
# Select a block of code and press shift+enter to run block selectively

# Lower merge
rmse_pos_mat_idm,rmse_vel_mat_idm,coll_mat_idm =
        multiscenarios_idm(mergetype="lower",modelmaker=make_IDM_models);
rmse_pos_mat_cidm,rmse_vel_mat_cidm, coll_mat_cidm =
        multiscenarios_idm(mergetype="lower",modelmaker=make_cidm_models);
rmse_pos_mat_lmidm, rmse_vel_mat_lmidm, coll_mat_lmidm =
        multiscenarios_idm(mergetype="lower",modelmaker=make_lmidm_models);
rmse_pos_mat_pf,rmse_vel_mat_pf, coll_mat_pf = multiscenarios_pf(mergetype="lower");

rmse_pos_list = [rmse_pos_mat_idm,rmse_pos_mat_cidm,rmse_pos_mat_lmidm,rmse_pos_mat_pf];
rmse_plots_modelscompare(rmse_pos_list,filename = "media/rmse_pos_lower.svg");

coll_mat_list = [coll_mat_idm,coll_mat_cidm,coll_mat_lmidm,coll_mat_pf];
coll_barchart(coll_mat_list,filename = "media/coll_barchart_lower.svg");

# Upper merge
rmse_pos_mat_idm,rmse_vel_mat_idm,coll_mat_idm =
        multiscenarios_idm(mergetype="upper",modelmaker=make_IDM_models);
rmse_pos_mat_cidm, rmse_vel_mat_cidm, coll_mat_cidm =
        multiscenarios_idm(mergetype="upper",modelmaker=make_cidm_models);
rmse_pos_mat_lmidm, rmse_vel_mat_lmidm, coll_mat_lmidm =
        multiscenarios_idm(mergetype="upper",modelmaker=make_lmidm_models);
rmse_pos_mat_pf, rmse_vel_mat_pf, coll_mat_pf = multiscenarios_pf(mergetype="upper");

rmse_pos_list = [rmse_pos_mat_idm,rmse_pos_mat_cidm,rmse_pos_mat_lmidm,rmse_pos_mat_pf];
rmse_plots_modelscompare(rmse_pos_list,filename = "media/rmse_pos_upper.svg");

coll_mat_list = [coll_mat_idm,coll_mat_cidm,coll_mat_lmidm,coll_mat_pf];
coll_barchart(coll_mat_list,filename = "media/coll_barchart_upper.svg");

#********************Train upper test lower******************
# We need to show a variability in the generated scenarios
# So we need to combine the particle sets of different vehicles together
# And then for the same set of vehicles in the lower merge i.e. test domain
# We show significantly different driving behavior by sampling from the particle set

#************LsqFit failed experimentation***************
# model(t, p) = p[1] * exp.(-p[2] * t)
# tdata = range(0,stop=10,length=20)
# ydata = model(tdata, [1.0 2.0]) + 0.01*randn(length(tdata))
# p0 = [0.5,0.5]
# fit = curve_fit(model, tdata, ydata, p0)

# function test_lmfit(d,p)
#         print("test_lmfit called\n")
#         print("says:d=$d\n")
#         t = d["t"]
#         return p[1]*exp.(-p[2]*t)
# end

# test_tdata = []
# for t in tdata
#         push!(test_tdata,Dict("t"=>t))
# end

# curve_fit(test_lmfit,test_tdata,ydata,[0.5,0.5])

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
# This file was copied from Transducers.jl
# which is available under an MIT license (see LICENSE).
using PkgBenchmark
using Markdown

function displayresult(result)
    md = sprint(export_markdown, result)
    md = replace(md, ":x:" => "❌")
    md = replace(md, ":white_check_mark:" => "✅")
    return display(Markdown.parse(md))
end

function printnewsection(name)
    println()
    println()
    println()
    printstyled("▃"^displaysize(stdout)[2]; color = :blue)
    println()
    printstyled(name; bold = true)
    println()
    return println()
end


text:
#println("hello world")
#mysuperdupermathvar = 3333 + 2222 + 11111
#println("type of that math var: ", typeof(mysuperdupermathvar))
#println("just wanna know something...")
#firstvar = readline()
#println("something was: ", firstvar)

# and this might be a comment
#print("enter something: ")
#line = readline()
#println("youve entered: ", line)

#using Sockets

#server = listen(IPv4(0), 8080) #0 listens on any ip
#while true
#  conn = accept(server)
#  @async begin
#    try
      #while true
        #line = readline(conn)
        #write(conn,line)
#      println("Hello World\n")
      #end
#    catch err
#      print("connection ended with error $err")
#    end
#  end
#end

#@async begin
#    server = listen(IPv4(0), 8080)
#    while true # always true... thats a never ending loop... old trick...
#        sock = accept(server)
#        println("Hello World\n")
#    end
#end

logfile = ""

samplearray = ["element 1", "element 2", "element 3", "element 4"]

# base julia funct struct
function readIndex()
    return read(open("index.html"), String)
end

# logfile
function openLog()
    logfile = open(string("log-", today(), ".log"), "w")
end

# using date functions
#using Dates

#serving
#using HTTP

#HTTP.listen("0.0.0.0", 8080) do http::HTTP.Stream
    #@show http.message #logs in terminal
    #@show HTTP.header(http, "Content-Type") #logs in terminal
    # more logging below
    #while !eof(http)
    #    println("body data: ", String(readIndex())) #save mode?
    #end
#    HTTP.setstatus(http, 200)
    #HTTP.setheader(http, "Foo-Header" => "bar") #man in the middle... literally
#    HTTP.startwrite(http)
#    write(http, readIndex())
    #write(http, "more response body")
#end

#using HTTP

#HTTP.serve("0.0.0.0", 8080) do request::HTTP.Request
#   @show request
#   @show request.method
#   @show HTTP.header(request, "Content-Type")
#   @show HTTP.payload(request)
#   try
#       return HTTP.Response(readIndex())
#   catch e
#       return HTTP.Response(404, "Error: $e")
#   end
#end

using Genie, Genie.Router
using Genie.Renderer, Genie.Renderer.Html, Genie.Renderer.Json

route("/") do
    html(Genie.Renderer.Html.read_template_file("./partial.html"), name="Oliver Leitner", samplearray=samplearray)
end

route("/hello") do
    #html(readIndex()) # reads cleanly... if you want text, you gotta say it
    html(Genie.Renderer.Html.read_template_file("./index.html"))
end

Genie.up(8080, "0.0.0.0", async=false)

max_stars_repo_name: OliverLeitner/juliademo

max_issues_repo_name: OliverLeitner/juliademo

max_forks_repo_name: OliverLeitner/juliademo


text:
import Dyno

filename = joinpath(rootdir,"examples","rbc.mod")

model_data = Dyno.modfile_parser(filename)

# print(model_data)


text:
## Residual Constrained Alternating Minimization
# Function Definitions

function pAltMinChol{ET<:Number}(b::AbstractArray{ET}, r::Int, K::Int)
# Performs RCAM on a single matricised monochromatic slice
# Calls: minFrobQR

    # Generate L0
    m,n = size(b)
    L = collect(ET, randn(m,r))

    indn = (1:n)::UnitRange{Int64}
    local R

    # Reduce the data down to nonzero elements of each column
    ### change to do in one loop ###
    indnz = [find(x -> x != 0, b[:,i]) for i in indn]
    bsub = [b[indnz[i],i] for i in indn]

    # This needs to be cleaned up for non-square entries
    indnz_tran = [find(x -> x != 0, b[i,:]) for i in indn]
    bsub_tran = [b[i,indnz_tran[i]] for i in indn]

    # Distribute the reduced data
    bsubD = distribute(bsub)
    bsubD_tran = distribute(bsub_tran)

    # Begin Alterations
    for k = 1:K

        # Solve for R
        R = pminFrobChol(L,bsubD,indn,indnz)

        # Solve for L
        L = pminFrobChol(R,bsubD_tran,indn,indnz_tran)

    end

    return L, R
end
export pAltMinChol

function pminFrobChol{ET<:Number}(L::AbstractArray{ET}, bsubD::AbstractArray, indi::UnitRange{Int64},
                    indnz::AbstractArray)
# Solves for each row of input
# Calls: minFrobRsub

    # Init L2 and bsub
    L2 = [L[indnz[i],:] for i in indi]

    # Distribute L2 and bsub
    LD = distribute(L2)

    # Map solving function to distributed inputs
    Rcol = map(pminFrobRsubChol, LD, bsubD)

    # Reform R to matrix
    R = Array{ET}(size(L))
    for i = 1:length(Rcol)
        R[i,:] = Rcol[i]
    end

    R = conj(R)

    return R
end

function pminFrobRsubChol{ET<:Number}(L2::AbstractArray{ET},bsub::AbstractArray{ET})
# Solves a single row, meant to be mapped over a DArray

    bsub = transpose(L2)*bsub

    # Chol fact
    L2 = L2'*L2

    # Faster to do in w/o mem alloc, even though it does chol twice
    v = backsub(chol(L2),fwdsub((chol(L2))',bsub))

    return v

end

function backsub{ET<:Number}(A::AbstractArray{ET}, y::AbstractArray{ET})

    m = size(A,2)
    x = zeros(ET,m,1)

    x[m] = y[m]/A[m,m]

    for i = m-1:-1:1
        x[i] = (y[i]-dot(A[i,i+1:m], x[i+1:m]))/A[i,i]
    end

    return x
end

function fwdsub{ET<:Number}(A::AbstractArray{ET},y::AbstractArray{ET})

    m = size(A,2)
    x = zeros(ET,m,1)
    x[1] = y[1]/A[1,1]

    for i = 2:1:m
        x[i] = (y[i]-dot(A[i,i-1:-1:1], x[i-1:-1:1]))/A[i,i]
    end

    return x
end


text:
using DifferentiableNAS
using Flux
using Flux: throttle, logitcrossentropy, onecold, onehotbatch, Optimiser
using Zygote: @nograd
using StatsBase: mean
using Parameters
using CUDA
using Distributions
using BSON
using Dates

include("../CIFAR10.jl")
include("../training_utils.jl")

@show args = parse_commandline()

train, val = get_processed_data(args["val_split"], args["batchsize"], 0.02f0, args["random_seed"])
@show size(train[1][1])

loss(m, x, y) = sum(m(x))

if args["random_seed"] > -1
    Random.seed!(args["random_seed"])
end
base_folder = prepare_folder("optest", args)

for _ in 1:2
    for prim in PRIMITIVES
        @show beginscript = now()
        global args
        global base_folder
        global cbepoch
        global cbbatch

        m = Chain([MixedOp(1,"1-2",3,1,[prim]) for _ in 1:8*14*8]...) |> gpu

        optimiser_α = Optimiser(WeightDecay(1f-3),ADAM(3f-4,(0.5f0,0.999f0)))
        optimiser_w = Optimiser(WeightDecay(3f-4),CosineAnnealing(args["epochs"]),Momentum(0.025f0, 0.9f0))

        losses = [0f0, 0f0]

        for epoch in 1:1
            @show epoch
            display(Dates.format(convert(DateTime,now()-beginscript), "HH:MM:SS"))
            @show prim
            @time begin
                w = Flux.params(m)
                for (train_batch, val_batch) in zip(TrainCuIterator(train), TrainCuIterator(val))
                    gsw = gradient(w) do
                        train_loss = loss(m, train_batch...)
                        return train_loss
                    end
                end
            end
        end
        display(("done", Dates.format(convert(DateTime,now()-beginscript), "HH:MM:SS")))
    end
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
"""
Download seismic data from server
May 23, 2019
Kurama Okubo
"""

using SeisIO, Printf, PlotlyJS, JLD2, FileIO

finame = "./data/SeisData_BP.jld2"
sr = 10 #downsampling for plotting
netstaID = 1 #to choose station
channelID = 1 #to choose station
#----------------#

# load data
@load finame S

#plotting
layout = Layout(width=1200, height=800, xaxis_range = [0, 60],
				yaxis_range = [-1, 25],
				xaxis=attr(title="Time [min]", dtick=10.0),
				yaxis=attr(title="Hours"),
                font =attr(size=18),
                showlegend=true,
                title =  S[1][netstaID][plotnscID].id)
p = plot([NaN], layout)

for stid = 1:size(S,1)

    requeststr = S[stid][netstaID][plotnscID].id
    sttime = u2d(S[stid][netstaID][plotnscID].t[1,2]*1e-6)

	println("test1")
    dt = 1.0./S[stid][netstaID][plotnscID].fs[1]
    t = LinRange(0.0:dt[1]:(S[stid][netstaID][plotnscID].t[2]-1)*dt[1])
    pt = t[1:sr:end]./60 # plot in [min]
	println("test2")

    #arbitrary amplitude normalization for plotting
    normalized_amp=2.0*maximum(abs.(S[1][netstaID][plotnscID].x))
	if maximum(abs.(S[stid][netstaID][plotnscID].x)) > 2.0*normalized_amp
		normalized_amp=0.2*maximum(abs.(S[stid][netstaID][plotnscID].x))
	end
	println("test3")

    trace1 = scatter(;x=pt, y= S[stid][netstaID][plotnscID].x[1:sr:end]./normalized_amp .+ (stid-1.0), mode="lines", line=attr(dash = false),
    name=@sprintf("%s: %s", requeststr, string(sttime)))
	println("test4")

    addtraces!(p, trace1)
end
deletetraces!(p, 1)

max_stars_repo_name: kura-okubo/SeisDownloadwJulia

max_issues_repo_name: kura-okubo/SeisDownloadwJulia

max_forks_repo_name: kura-okubo/SeisDownloadwJulia


text:
#=
    Support for >= and >
=#
const UnionGT{T} = Union{CPE.Strictly{MOI.GreaterThan{T}, T}, MOI.GreaterThan{T}}

function supports_concreteB(concrete_B)
    added_constraints = MOIB.added_constraint_types(concrete_B)
    length(added_constraints) > 1 && return false
    # The inner constraint should not create any variable (might have unexpected consequences)
    return isempty(MOIB.added_constrained_variable_types(concrete_B))
end

function get_num_vars(fct::SAF)
    return length(fct.terms)
end

function get_num_vars(fct::VAF)
    return length(fct.terms)
end

function get_num_vars(fct::MOI.VectorOfVariables)
    return length(fct.variables)
end

"""
    map_function(bridge, fct, set)

Default for when set is not needed so return `MOIBC.map_function(bridge, fct)`
If the set is needed a specific method should be implemented
"""
function map_function(
    bridge::Type{<:MOIBC.AbstractBridge},
    fct,
    set
)
    MOIBC.map_function(bridge, fct)
end

"""
    added_constraint_types(bridge, set)

Default for when set is not needed so return `MOIBC.added_constraint_types(bridge)`
If the set is needed a specific method should be implemented
"""
function added_constraint_types(
    B::Type{<:MOIBC.AbstractBridge},
    S,
)
    return MOIB.added_constraint_types(B)
end


text:
"""
    ArrayModel{T <: MillFunction} <: AbstractMillModel

A model node for processing [`ArrayNode`](@ref)s. It applies a (sub)model `m` stored in it to data of
the [`ArrayNode`](@ref).

# Examples
```jldoctest array_model
julia> Random.seed!(0);
```

```jldoctest array_model; filter=$(DOCTEST_FILTER)
julia> n = ArrayNode(randn(Float32, 2, 2))
2×2 ArrayNode{Matrix{Float32}, Nothing}:
 0.679...  -0.353...
 0.828...  -0.135...
```

```jldoctest array_model
julia> m = ArrayModel(Dense(2, 2))
ArrayModel(Dense(2, 2))
```

```jldoctest array_model; filter=$(DOCTEST_FILTER)
julia> m(n)
2×2 ArrayNode{Matrix{Float32}, Nothing}:
 0.661...  -0.188...
 0.101...   0.275...
```

See also: [`AbstractMillModel`](@ref), [`IdentityModel`](@ref), [`identity_model`](@ref), [`ArrayNode`](@ref).
"""
struct ArrayModel{T <: MillFunction} <: AbstractMillModel
    m::T
end

Flux.@functor ArrayModel

# (m::ArrayModel)(x::ArrayNode) = ArrayNode(m.m(getfield(x, :data)))

function (m::ArrayModel)(x::ArrayNode)
    c = getfield(m, :m)
    a = getfield(x, :data)
    ArrayNode(c(a))
end

"""
    identity_model()

Returns an [`ArrayModel`](@ref) realising the `identity` transformation.

# Examples
```jldoctest
julia> identity_model()
ArrayModel(identity)
```

See also: [`ArrayModel`](@ref), [`IdentityModel`](@ref).
"""
identity_model() = ArrayModel(identity)

"""
    IdentityModel

Alias for `ArrayModel{typeof(identity)}`.

See also: [`ArrayModel`](@ref), [`identity_model`](@ref).
"""
const IdentityModel = ArrayModel{typeof(identity)}

# Base.hash(m::ArrayModel{T}, h::UInt) where {T} = hash((T, m.m), h)
# (m1::ArrayModel{T} == m2::ArrayModel{T}) where {T} = m1.m == m2.m


text:
"""
    rectangular_meander!(c::Cell{T,S}, length, nlegs::Int,
        conductor_width, conductor_spacing, meta::Meta=GDSMeta(0,0);
        first_length=length/2, last_length=length/2, start_upward=true)
Create a positive trace of rectangular meander specified by its
vertical length (`length`), number of legs (`nlegs::Int`), conductor
width and spacing (`conductor_width` and `conductor_spacing`). The
first and last length are initialized as half the `length`, and can
be specified with keyword arguments `first_length` and `last_length`.
The meander starts upward unless the keyword argument
`start_upward' is enforced to `false`.
"""
function rectangular_meander!(c::Cell{T,S}, length, nlegs::Int,
    conductor_width, conductor_spacing, meta::Meta=GDSMeta(0,0);
    first_length=length/2, last_length=length/2,
    start_upward=true) where {T<:Coordinate, S<:GDSMeta}
    if nlegs < 2
        error("The number of legs must be 2 or larger.")
    end

    # starting direction:
    # +1 if start_upward=true, -1 if start_upward=false
    dir = 2 * start_upward - 1

    p = Path(Point(conductor_width/2, zero(T)), α0=dir*π/2)

    path_style = Paths.Trace(conductor_width)
    corner_style = Paths.SimpleTraceCorner()

    straight!(p, first_length + conductor_width/2, path_style)
    corner!(p, -dir * π/2, corner_style)
    straight!(p, (conductor_spacing + conductor_width)/2, path_style)
    for n in 1:(nlegs - 2)
        straight!(p, (conductor_spacing + conductor_width)/2, path_style)
        corner!(p, dir * (-1)^n * π/2, corner_style)
        straight!(p, conductor_width + length, path_style)
        corner!(p, dir * (-1)^n * -π/2, corner_style)
        straight!(p, (conductor_spacing + conductor_width)/2, path_style)
    end
    straight!(p, (conductor_spacing + conductor_width)/2, path_style)
    corner!(p, dir * (-1)^(nlegs-1) * π/2, corner_style)
    straight!(p, conductor_width/2 + last_length, path_style)

    render!(c, p, meta)
end


text:
export CSMOccuranceType
export parseCSMVerboseLog, calcCSMOccurancesFolders, calcCSMOccuranceMax, printCSMOccuranceMax, reconstructCSMHistoryLogical

# [cliqId][fsmIterNumber][fsmFunctionName] => (nr. call occurances, list global call sequence position, list of status)
const CSMOccuranceType = Dict{Int, Dict{Int, Dict{Symbol, Tuple{Int, Vector{Int}, Vector{String}}}}}

function parseCSMVerboseLog(resultsDir::AbstractString;verboseName::AbstractString="csmVerbose.log")
  #
  fid = open(joinpath(resultsDir, verboseName), "r")
  fsmLines = readlines(fid)
  close(fid)

  # parse lines into usable format
  sfsmL = split.(fsmLines, r" -- ")
  cids = split.(sfsmL .|> x->match(r"cliq\d+", x[1]).match, r"cliq") .|> x->parse(Int,x[end])
  iters = split.(sfsmL .|> x->match(r"iter=\d+", x[1]).match, r"iter=") .|> x->parse(Int,x[end])
  smfnc = sfsmL .|> x->split(x[2], ',')[1] .|> Symbol
  statu = sfsmL .|> x->split(x[2], ',')[2] .|> x->lstrip(rstrip(x))

  return cids, iters, smfnc, statu
end

## Make lookup from all runs
function calcCSMOccurancesFolders(folderList::Vector{<:AbstractString};
                                  verboseName::AbstractString="csmVerbose.log" )
  #
  # lookup for histogram on each step per fsm
  # [cliqId][fsmIterNumber][fsmFunctionName] => (nr. call occurances, list global call sequence position)
  csmCounter = CSMOccuranceType()
  # lookup for transition counts per fsm function
  trxCounter = Dict{Symbol, Dict{Symbol, Int}}()

  prevFnc = Dict{Int, Symbol}()

  for rDir in folderList
    ## load the sequence from each file
    cids, iters, smfnc, statu = parseCSMVerboseLog(rDir, verboseName=verboseName)

    # populate histogram
    for (idx,smfi) in enumerate(smfnc)
      if !haskey(csmCounter, cids[idx])
        csmCounter[cids[idx]] = Dict{Int, Dict{Symbol, Tuple{Int,Vector{Int}}}}()
      end
      if !haskey(csmCounter[cids[idx]], iters[idx])
        # Tuple{Int,Int[]} == (nr. occurances of call, list global call sequence position)
        csmCounter[cids[idx]][iters[idx]] = Dict{Symbol, Tuple{Int,Vector{Int}, Vector{String}}}()
      end
      easyRef = csmCounter[cids[idx]][iters[idx]]
      if !haskey(easyRef,smfi)
        easyRef[smfi] = (0,Int[],String[])
      end
      # add position in call sequence (global per solve)
      globalSeqIdx = easyRef[smfi][2]
      push!(globalSeqIdx, idx)
      statSeq = easyRef[smfi][3]
      push!(statSeq, statu[idx])
      easyRef[smfi] = (easyRef[smfi][1]+1, globalSeqIdx, statSeq)

      ## also track the transitions
      if haskey(prevFnc, cids[idx])
        if !haskey(trxCounter, prevFnc[cids[idx]])
          # add function lookup if not previously seen
          trxCounter[prevFnc[cids[idx]]] = Dict{Symbol, Int}()
        end
        if !haskey(trxCounter[prevFnc[cids[idx]]], smfi)
          # add previously unseen transition
          trxCounter[prevFnc[cids[idx]]][smfi] = 0
        end
        # from previous to next function
        trxCounter[prevFnc[cids[idx]]][smfi] += 1
      end
      # always update prevFnc register
      prevFnc[cids[idx]] = smfi
    end
  end

  return csmCounter, trxCounter
end

"""
    $SIGNATURES

Use maximum occurance from `csmCounter::CSMOccuranceType` to summarize many CSM results.

Notes
- `percentage::Bool=false` shows median global sequence occurance ('m'), or
  - `percentage::Bool=true` of occurance ('%')
"""
function calcCSMOccuranceMax( csmCounter::CSMOccuranceType;
                              percentage::Bool=false)
  #
  ncsm = length(keys(csmCounter))
  maxOccuran = Dict()
  # max steps
  for i in 1:ncsm
    # sequence of functions that occur most often
    maxOccuran[i] = Vector{Tuple{Symbol, String, String}}()
  end

  # pick out the max for each CSM iter
  for (csmID, csmD) in csmCounter, stp in 1:length(keys(csmD))
    maxFnc = :null
    maxCount = 0
    totalCount = 0
    for (fnc, cnt) in csmCounter[csmID][stp]
      totalCount += cnt[1]
      if maxCount < cnt[1]
        maxCount = cnt[1]
        maxFnc = fnc
      end
    end
    # occurance count
    perc = if percentage
      "$(round(Int,(maxCount/totalCount)*100))"
    else
      # get medial position (proxy to most frequent)
      "$(round(Int,Statistics.median(csmCounter[csmID][stp][maxFnc][2])))"
    end
    # get status
      allst = csmCounter[csmID][stp][maxFnc][3]
      qst = unique(allst)
      mqst = qst .|> y->count(x->x==y, allst)
      midx = findfirst(x->x==maximum(mqst),mqst)
    maxStatus = qst[midx]
    push!(maxOccuran[csmID], (maxFnc, perc, maxStatus) ) # position in vector == stp
  end
  maxOccuran
end

"""
    $SIGNATURES

Print the most likely FSM function at each step per state machine, as swim lanes.

Example

```julia
csmCo = calcCSMOccurancesFolders(resultFolder[maskTrue])
maxOcc = calcCSMOccuranceMax(csmCo)
printCSMOccuranceMax(maxOcc)
```
"""
function printCSMOccuranceMax(maxOcc;
                              fid=stdout,
                              percentage::Bool=false )
  #
  ncsm = length(keys(maxOcc))

  # print titles
  titles = Tuple[]
  for cid in 1:ncsm
    tpl = ("","","$cid                   ","    ")
    push!(titles, tpl)
  end
  IIF.printHistoryLane(fid, "", titles)
  print(fid,"----")
  for i in 1:ncsm
    print(fid,"+--------------------")
  end
  println(fid,"")

  maxsteps=0
  for i in 1:ncsm
    maxsteps = maxsteps < length(maxOcc[i]) ? length(maxOcc[i]) : maxsteps
  end

  for stp in 1:maxsteps
    TPL = Tuple[]
    for cid in 1:ncsm
      tpl = ("","","                     ","    ")
      if stp <= length(maxOcc[cid])
        fncName = maxOcc[cid][stp][1]
        # either show percentage or sequence index
        percOrSeq = "$(maxOcc[cid][stp][2])"
        percOrSeq *= percentage ? "%" : "m"
        # get status
        tpl = ("",percOrSeq,fncName,maxOcc[cid][stp][3])
      end
      push!(TPL, tpl)
    end
    IIF.printHistoryLane(fid, stp, TPL)
  end
end

"""
    $SIGNATURES

Use `solveTree!`'s` `verbose` output to reconstruct the swim lanes Logical sequence of CSM function calls.

Notes
- This is a secondary function to primary `printCSMHistoryLogical`.

Related

printCSMHistoryLogical
"""
function reconstructCSMHistoryLogical(resultsDir::AbstractString;
                                      fid::IO=stdout,
                                      verboseName::AbstractString="csmVerbose.log" )
  #
  csmCounter, trxCounter = calcCSMOccurancesFolders([resultsDir], verboseName=verboseName)

  # print with sequence position
  maxOcc = calcCSMOccuranceMax(csmCounter, percentage=false)

  printCSMOccuranceMax(maxOcc, fid=fid)
end

#


text:
module TestGRM
using Test
using TMLEEpistasis

@testset "Test readGRM" begin
    prefix = joinpath("data", "grm", "test.grm")
    GRM, ids = TMLEEpistasis.readGRM(prefix)
    @test size(GRM, 1) == 18915
    @test size(ids, 1) == 194
end

end;

true


text:
using StochasticSearch, FactCheck

facts("[Run]") do
    context("constructors") do
        tuning_run = Run()
        @fact (typeof(tuning_run) <: Run) --> true
    end
end


text:
#######################################################################
# TightbindingModelTerm
#######################################################################
abstract type TightbindingModelTerm end
abstract type AbstractOnsiteTerm <: TightbindingModelTerm end
abstract type AbstractHoppingTerm <: TightbindingModelTerm end

struct OnsiteTerm{F,
                  S<:Union{Missing,Tuple{Vararg{NameType}}},
                  C} <: AbstractOnsiteTerm
    o::F
    sublats::S
    coefficient::C
    forcehermitian::Bool
end

struct HoppingTerm{F,
                   S<:Union{Missing,Tuple{Vararg{Tuple{NameType,NameType}}}},
                   D<:Union{Missing,Tuple{Vararg{SVector{L,Int}}} where L},
                   R<:Union{Missing,AbstractFloat},
                   C} <: AbstractHoppingTerm
    h::F
    sublats::S
    dns::D
    range::R
    coefficient::C
    forcehermitian::Bool
end

(o::OnsiteTerm{<:Function})(r,dr) = o.coefficient * o.o(r)
(o::OnsiteTerm)(r,dr) = o.coefficient * o.o

(h::HoppingTerm{<:Function})(r, dr) = h.coefficient * h.h(r, dr)
(h::HoppingTerm)(r, dr) = h.coefficient * h.h

sanitize_sublats(s::Missing) = missing
sanitize_sublats(s::Integer) = (nametype(s),)
sanitize_sublats(s::NameType) = (s,)
sanitize_sublats(s::Tuple) where {N} = nametype.(s)
sanitize_sublats(s::Tuple{}) = ()
sanitize_sublats(n) = throw(ErrorException(
    "`sublats` for `onsite` must be either `missing`, an `s` or a tuple of `s`s, with `s::$NameType` is a sublattice name"))

sanitize_sublatpairs(s::Missing) = missing
sanitize_sublatpairs((s1, s2)::NTuple{2,Union{Integer,NameType}}) = ((nametype(s1), nametype(s2)),)
sanitize_sublatpairs((s2, s1)::Pair) = (sanitize_sublatpairs(s1, s2),)
sanitize_sublatpairs(s::Union{Integer,NameType}) = sanitize_sublatpairs((s,s))
sanitize_sublatpairs(s::NTuple{N,Any}) where {N} =
    ntuple(n -> first(sanitize_sublatpairs(s[n])), Val(N))
sanitize_sublatpairs(s) = throw(ErrorException(
    "`sublats` for `hopping` must be either `missing`, a tuple `(s₁, s₂)`, or a tuple of such tuples, with `sᵢ::$NameType` a sublattice name"))

sanitize_dn(dn::Missing) = missing
sanitize_dn(dn::Tuple{Vararg{Tuple}}) = SVector.(dn)
sanitize_dn(dn::Tuple{Vararg{Integer}}) = (SVector(dn),)
sanitize_dn(dn::Tuple{}) = ()

sublats(t::OnsiteTerm{<:Any,Missing}, lat::AbstractLattice) = collect(1:nsublats(lat))
function sublats(t::OnsiteTerm{<:Any,<:Tuple}, lat::AbstractLattice)
    names = lat.unitcell.names
    s = Int[]
    for name in t.sublats
        i = findfirst(isequal(name), names)
        i !== nothing && push!(s, i)
    end
    return s
end

sublats(t::HoppingTerm{<:Any,Missing}, lat::AbstractLattice) =
    vec(collect(Iterators.product(1:nsublats(lat), 1:nsublats(lat))))
function sublats(t::HoppingTerm{<:Any,<:Tuple}, lat::AbstractLattice)
    names = lat.unitcell.names
    s = Tuple{Int,Int}[]
    for (n1, n2) in t.sublats
        i1 = findfirst(isequal(n1), names)
        i2 = findfirst(isequal(n2), names)
        i1 !== nothing && i2 !== nothing && push!(s, (i1, i2))
    end
    return s
end

displayparameter(::Type{<:Function}) = "Function"
displayparameter(::Type{T}) where {T} = "$T"

function Base.show(io::IO, o::OnsiteTerm{F}) where {F}
    i = get(io, :indent, "")
    print(io,
"$(i)OnsiteTerm{$(displayparameter(F))}:
$(i)  Sublattices      : $(o.sublats === missing ? "any" : o.sublats)
$(i)  Force hermitian  : $(o.forcehermitian)
$(i)  Coefficient      : $(o.coefficient)")
end

function Base.show(io::IO, h::HoppingTerm{F}) where {F}
    i = get(io, :indent, "")
    print(io,
"$(i)HoppingTerm{$(displayparameter(F))}:
$(i)  Sublattice pairs : $(h.sublats === missing ? "any" : (t -> Pair(reverse(t)...)).(h.sublats))
$(i)  dn cell jumps    : $(h.dns === missing ? "any" : h.dns)
$(i)  Hopping range    : $(round(h.range, digits = 6))
$(i)  Force hermitian  : $(h.forcehermitian)
$(i)  Coefficient      : $(h.coefficient)")
end

# External API #
"""
    onsite(o; sublats = missing, forcehermitian = true)

Create an `TightbindingModelTerm` that applies an onsite energy `o` to a `Lattice` when
creating a `Hamiltonian` with `hamiltonian`.

The onsite energy `o` can be a number, a matrix (preferably `SMatrix`) or a function of the
form `r -> ...` for a position-dependent onsite energy. If `sublats` is specified as a
sublattice name or tuple thereof, `onsite` is only applied to sublattices with said names.
If `forcehermitian` is true, the model will produce an hermitian Hamiltonian.

The dimension of `o::AbstractMatrix` must match the orbital dimension of applicable
sublattices (see also `orbitals` option for `hamiltonian`). If `o::Number` it will be
treated as `o * I` (proportional to identity matrix) when applied to multiorbital
sublattices.

`TightbindingModelTerm`s created with `onsite` or `hopping` can be added or substracted
together to build more complicated `TightbindingModel`s.

# Examples
```
julia> onsite(1, sublats = (1,2)) - hopping(2)
TightbindingModel{2}: model with 2 terms
  OnsiteTerm{Int64}:
    Sublattices      : (1, 2)
    Force hermitian  : true
    Coefficient      : 1
  HoppingTerm{Int64}:
    Sublattice pairs : any
    dn cell jumps    : any
    Hopping range    : 1.0
    Force hermitian  : true
    Coefficient      : -1

julia> hamiltonian(LatticePresets.honeycomb(orbitals = (:a, :b)), onsite(r->@SMatrix[1 2; 3 4]))
Hamiltonian{<:Lattice} : 2D Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 1 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Elements         : 2 × 2 blocks
  Onsites          : 2
  Hoppings         : 0
  Coordination     : 0.0
```

# See also:
    `hopping`
"""
function onsite(o; sublats = missing, forcehermitian::Bool = true)
    return TightbindingModel(OnsiteTerm(o, sanitize_sublats(sublats), 1, forcehermitian))
end

"""
    hopping(h; sublats = missing, range = 1, dn = missing, forcehermitian = true)

Create an `TightbindingModelTerm` that applies a hopping `h` to a `Lattice` when
creating a `Hamiltonian` with `hamiltonian`.

The maximum distance between coupled sites is given by `range::Real`. If a cell distance
`dn::NTuple{L,Int}` or distances `dn::NTuple{M,NTuple{L,Int}}` are given, only unit cells
at that distance will be coupled.

The hopping amplitude `h` can be a number, a matrix (preferably `SMatrix`) or a function
of the form `(r, dr) -> ...` for a position-dependent hopping (`r` is the bond center,
and `dr` the bond vector). If `sublats` is specified as a sublattice name pair, or tuple
thereof, `hopping` is only applied between sublattices with said names. If `forcehermitian`
is true, the model will produce an hermitian Hamiltonian.

The dimension of `h::AbstractMatrix` must match the orbital dimension of applicable
sublattices (see also `orbitals` option for `hamiltonian`). If `h::Number` it will be
treated as `h * I` (proportional to identity matrix) when applied to multiorbital
sublattices.

`TightbindingModelTerm`s created with `onsite` or `hopping` can be added or substracted
together to build more complicated `TightbindingModel`s.

# Examples
```
julia> onsite(1) - hopping(2, dn = ((1,2), (0,0)), sublats = (1,1))
TightbindingModel{2}: model with 2 terms
  OnsiteTerm{Int64}:
    Sublattices      : any
    Force hermitian  : true
    Coefficient      : 1
  HoppingTerm{Int64}:
    Sublattice pairs : (1 => 1,)
    dn cell jumps    : ([1, 2], [0, 0])
    Hopping range    : 1.0
    Force hermitian  : true
    Coefficient      : -1

julia> hamiltonian(LatticePresets.honeycomb(), hopping((r,dr) -> cos(r[1]), sublats = ((1,1), (2,2))))
Hamiltonian{<:Lattice} : 2D Hamiltonian on a 2D Lattice in 2D space
  Bloch harmonics  : 7 (SparseMatrixCSC, sparse)
  Harmonic size    : 2 × 2
  Elements         : scalars
  Onsites          : 0
  Hoppings         : 12
  Coordination     : 6.0
```

# See also:
    `onsite`
"""
function hopping(h; sublats = missing, range::Real = 1, dn = missing, forcehermitian::Bool = true)
    return TightbindingModel(HoppingTerm(h, sanitize_sublatpairs(sublats), sanitize_dn(dn),
        float(range) + sqrt(eps(float(range))), 1, forcehermitian))
end

Base.:*(x, o::OnsiteTerm) =
    OnsiteTerm(o.o, o.sublats, x * o.coefficient, o.forcehermitian)
Base.:*(x, t::HoppingTerm) =
    HoppingTerm(t.h, t.sublats, t.dns, t.range, x * t.coefficient, t.forcehermitian)
Base.:*(t::TightbindingModelTerm, x) = x * t
Base.:-(t::TightbindingModelTerm) = (-1) * t

LinearAlgebra.ishermitian(t::OnsiteTerm) = t.forcehermitian
LinearAlgebra.ishermitian(t::HoppingTerm) = t.forcehermitian

# Base.:+(t1::TightbindingModelTerm, t2::TightbindingModelTerm) = TightbindingModel((t1, t2))
# Base.:-(t1::TightbindingModelTerm, t2::TightbindingModelTerm) = TightbindingModel((t1, -t2))

#######################################################################
# TightbindingModel
#######################################################################
struct TightbindingModel{N,T<:Tuple{Vararg{TightbindingModelTerm,N}}}
    terms::T
end

terms(t::TightbindingModel) = t.terms

TightbindingModel(ts::TightbindingModelTerm...) = TightbindingModel(ts)

(m::TightbindingModel)(r, dr) = sum(t -> t(r, dr), m.terms)

# External API #

Base.:*(x, m::TightbindingModel) = TightbindingModel(x .* m.terms)
Base.:*(m::TightbindingModel, x) = x * m
Base.:-(m::TightbindingModel) = TightbindingModel((-1) .* m.terms)

Base.:+(m::TightbindingModel, t::TightbindingModel) = TightbindingModel((m.terms..., t.terms...))
Base.:-(m::TightbindingModel, t::TightbindingModel) = m + (-t)

function Base.show(io::IO, m::TightbindingModel{N}) where {N}
    ioindent = IOContext(io, :indent => "  ")
    print(io, "TightbindingModel{$N}: model with $N terms", "\n")
    foreach(t -> print(ioindent, t, "\n"), m.terms)
end

LinearAlgebra.ishermitian(m::TightbindingModel) = all(t -> ishermitian(t), m.terms)

#######################################################################
# nondiagonal
#######################################################################

struct NondiagonalTerm{T<:HoppingTerm,S<:Tuple{Vararg{UnitRange{Int}}}} <: AbstractHoppingTerm
    term::T
    sublatranges::S
end

function nondiagonal(m::TightbindingModel{N}, nsublats::Tuple{Vararg{Int}}) where {N}
    sranges = sublatranges(nsublats...)
    return TightbindingModel(ntuple(i -> NondiagonalTerm(m.terms[i], sranges), Val(N)))
end

sublatranges(i::Int, is::Int...) = _sublatranges((1:i,), is...)
_sublatranges(rs::Tuple, i::Int, is...) = _sublatranges((rs..., last(last(rs)) + 1: last(last(rs)) + i), is...)
_sublatranges(rs::Tuple) = rs

function Base.show(io::IO, h::NondiagonalTerm)
    i = get(io, :indent, "")
    show(io, h.term)
    print(io,"\n$(i)  Nondiagonal only : $(h.sublatranges)")
end

Base.:*(x, t::NondiagonalTerm) = NondiagonalTerm(x * t.term)

LinearAlgebra.ishermitian(t::NondiagonalTerm) = ishermitian(t.term)

function sublats(t::NondiagonalTerm, lat::AbstractLattice)
    s = sublats(t.term, lat)
    sr = t.sublatranges
    filter!(spair ->  findblock(first(spair), sr) != findblock(last(spair), sr), s)
    return s
end

findblock(s, sr) = findfirst(r -> s in r, sr)

#######################################################################
# checkmodelorbs - check for inconsistent orbital dimensions
#######################################################################
checkmodelorbs(model::TightbindingModel, orbs, lat) =
    foreach(term -> _checkmodelorbs(term, orbs, lat), model.terms)

function _checkmodelorbs(term::HoppingTerm, orbs, lat)
    for (s1, s2) in sublats(term, lat)
        _checkmodelorbs(term(first(sites(lat, s1)), first(sites(lat, s2))), length(orbs[s1]), length(orbs[s2]))
    end
    return nothing
end

function _checkmodelorbs(term::OnsiteTerm, orbs, lat)
    for s in sublats(term, lat)
        _checkmodelorbs(term(first(sites(lat, s)), first(sites(lat, s))), length(orbs[s]))
    end
    return nothing
end

_checkmodelorbs(nt::NondiagonalTerm, orbs, lat) = _checkmodelorbs(nt.term, orbs, lat)

_checkmodelorbs(s::SMatrix, m, n = m) =
    size(s) == (m, n) || @warn("Possible dimension mismatch between model and Hamiltonian. Did you correctly specify the `orbitals` in hamiltonian?")

_checkmodelorbs(s::Number, m, n = m) = _checkmodelorbs(SMatrix{1,1}(s), m, n)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the transformations between Geodetic and Geocentric frames.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/transformations/geodetic_geocentric.jl
# ==================================================

# Function: ECEFtoGeodetic
# ------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-3: Converting ECEF to Lat Lon [1, p. 173].
#
# According to this example, using:
#
#   r_ecef = 6524.834 i + 6862.875 j + 6448.296 k [km]
#
# one gets:
#
#   Geodetic Latitude  = 34.352496°
#            Longitude = 46.4464°
#            Altitude  = 5085.22 km
#
# Scenario 02
# ===========
#
# At the poles, we have a singularity. We know that, if:
#
#   r_ecef = 0 i + 0 j + Z k
#
# then
#
#   Geodetic Latitude  = 90° for Z > 0 and -90° for Z < 0
#            Longitude =  0°
#            Altitude  = Z - b_wgs84
#
################################################################################

@testset "Function ECEFtoGeodetic" begin
    # Scenario 01
    # ===========

    r = [6524.834e3, 6862.875e3, 6448.296e3]

    ϕ_gd, λ_gd, h = ECEFtoGeodetic(r)

    @test rad2deg(ϕ_gd) ≈ 34.352496 atol = 1e-6
    @test rad2deg(λ_gd) ≈ 46.4464   atol = 1e-4
    @test h/1000        ≈ 5085.22   atol = 1e-2

    # Scenario 02
    # ===========

    aux = rand(0:1000)

    Z = R0 + aux
    ϕ_gd, λ_gd, h = ECEFtoGeodetic([0;0;Z])

    @test rad2deg(ϕ_gd) ≈ 90
    @test rad2deg(λ_gd) ≈ 0
    @test h             ≈ Z - SatelliteToolbox.b_wgs84

    Z = -R0 + aux
    ϕ_gd, λ_gd, h = ECEFtoGeodetic([0;0;Z])

    @test rad2deg(ϕ_gd) ≈ -90
    @test rad2deg(λ_gd) ≈ 0
    @test h             ≈ -Z - SatelliteToolbox.b_wgs84
end


text:
s = "Hello"
s *= ", world!"


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@doc doc"""
    GenericDensity{F<:Function} <: AbstractDensity

Constructors:

    GenericDensity(f)

Turns the density function `f` into a BAT-compatible [`AbstractDensity)(@ref).
`f(v)` must return either a [`LogDVal`](@ref) (recommended) or a
[`LinDVal`](@ref).

It must be safe to execute `f` in parallel on multiple threads and
processes.
"""
struct GenericDensity{F<:Function} <: AbstractDensity
    f::F
end

Base.convert(::Type{GenericDensity}, f::Function) = GenericDensity(f)
Base.convert(::Type{AbstractDensity}, f::Function) = GenericDensity(f)

Base.parent(density::GenericDensity) = density.f

function eval_logval_unchecked(density::GenericDensity, v::Any)
    logvalof(density.f(v))
end


text:
import Mmap

function write_vector(path,filename,head,data)
    open(joinpath(path,filename), "w+") do s
        write(s, head)
        write(s, length(data))
        write(s, data)
    end
end

function read_vector(path,filename)
    A=Nothing
    open(joinpath(path,filename)) do s  # default is read-only
        head=read(s,Int)
        m = read(s, Int)
        A = Mmap.mmap(s, Vector{Float64}, m)
    end
    return A
end

function write_matrix(path,filename,head,data)
    open(joinpath(path,filename), "w+") do s
        write(s, head)
        write(s, size(data,1))
        write(s, size(data,2))
        write(s, data)
    end
end

function read_matrix(path,filename)
    A=Nothing
    open(joinpath(path,filename)) do s  # default is read-only
        head=read(s,Int)
        m = read(s, Int)
        n = read(s, Int)
        A = Mmap.mmap(s, Matrix{Float64}, (m,n))
    end
    return A
end


text:
"""
usage:
$ julia -p auto -L buffon.jl
$ buffon(convert(Int,1e9))
$ buffon_parallel(convert(Int,1e9))
"""

using BenchmarkTools

function buffon(n_trial)
    hit_counter = zero(Bool)
    for i in 1:n_trial
        position = rand()
        φ = (rand() * π) - π / 2
        x_right = position + cos(φ) / 2
        x_left = position - cos(φ) / 2
        is_hit = (x_right >= 1 || x_left <= 0) ? 1 : 0
        hit_counter += is_hit
    end
    approx_π = n_trial/hit_counter * 2
end

using Distributed

function buffon_parallel(n_trial)
    hit_counter = @distributed (+) for i in 1:n_trial
        position = rand()
        φ = (rand() * π) - π / 2
        x_right = position + cos(φ) / 2
        x_left = position - cos(φ) / 2
        (x_right >= 1 || x_left <= 0) ? 1 : 0
    end
    approx_π = n_trial/hit_counter * 2
end

function main()
    @benchmark buffon(convert(Int,1e9))
    @benchmark buffon_parallel(convert(Int,1e9))
end

#main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
import HypothesisTests: OneSampleTTest

# Example system
sys = logistic2_unidir(c_xy = 1.0)

# Run test on time series with 200 observations, treating the first
# variable as the source and the second variable as the target
setup = DiscreteSystemSetup(source = 1, target = 2, n_pts = 200)

# Causality tests
test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
test_cm = CrossMappingTest(n_reps = 10)
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
test_jdd = JointDistanceDistributionTest()
test_jddt = JointDistanceDistributionTTest()
predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
test_pa = PredictiveAsymmetryTest(predictive_test = predtest)

# Run test
@test causality(sys, setup, test_ccm) isa Vector{Vector{T}} where {T <: Real}
@test causality(sys, setup, test_cm) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_vf) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_tog) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_jdd) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_jddt) isa OneSampleTTest


text:
# # Stability

# In this example, we will demonstrate the stability of persistent homology. The stability
# theorem roughly states that a small change in the input data will result in a small
# change in the resulting persistence diagram. In other words, persistent homology is very
# tolerant of noisy data. Also, see the
# [Distances](https://mtsch.github.io/PersistenceDiagrams.jl/dev/generated/distances/)
# example in [PersistenceDiagrams.jl](https://github.com/mtsch/PersistenceDiagrams.jl).

# Again, start with loading some packages.

using Ripserer
using Plots
using Random; # hide
Random.seed!(1337); # hide
gr(); # hide
nothing; # hide

# As in the Basics example, we will look at the persistent homology of a noisy circle.

function noisy_circle(n; r1=1, r2=1, noise=0.1)
    points = NTuple{2,Float64}[]
    for _ in 1:n
        θ = 2π * rand()
        point = (
            r1 * sin(θ) + noise * rand() - noise / 2,
            r2 * cos(θ) + noise * rand() - noise / 2,
        )
        push!(points, point)
    end
    return points
end

# We will look at the first persistent homology group of this space.

# First, let's see what happens if we repeatedly sample 100 random points from a circle.

anim = @animate for _ in 1:200
    points = noisy_circle(100; noise=0)
    result = ripserer(points)

    plt_pts = scatter(
        points;
        legend=false,
        aspect_ratio=1,
        xlim=(-2.2, 2.2),
        ylim=(-2.2, 2.2),
        title="Data",
    )
    plt_diag = plot(result; infinity=3)

    plot(plt_pts, plt_diag; size=(800, 400))
end
gif(anim, "stability_anim_1.gif") # hide

# We notice that an interval in ``H_1`` always stands out and that its death remains
# constant. The only thing that changes is the birth time. The birth time is equal to the
# largest distance between adjacent points in the circle. At the birth time, the circle is
# connected.

# Now, let's add some noise!

anim = @animate for _ in 1:200
    points = noisy_circle(100; noise=0.2)
    result = ripserer(points)

    plt_pts = scatter(
        points;
        legend=false,
        aspect_ratio=1,
        xlim=(-2.2, 2.2),
        ylim=(-2.2, 2.2),
        title="Data",
    )
    plt_diag = plot(result; infinity=3)

    plot(plt_pts, plt_diag; size=(800, 400))
end
gif(anim, "stability_anim_2.gif") # hide

# The interval is jumping around a lot more now, but it hovers around the same general area.
# It's still clearly the most persistent feature of our space.

# Next, let's look at how adding more and more noise affects the diagram.

anim = @animate for noise in vcat(0:0.01:1, 1:-0.01:0)
    points = noisy_circle(100; noise=noise)
    result = ripserer(points)

    plt_pts = scatter(
        points;
        legend=false,
        aspect_ratio=1,
        xlim=(-2.2, 2.2),
        ylim=(-2.2, 2.2),
        title="Data",
    )
    plt_diag = plot(result; infinity=3)

    plot(plt_pts, plt_diag; size=(800, 400))
end
gif(anim, "stability_anim_3.gif") # hide

# We see we have to add quite a bit of noise to destroy the diagram. Notice how the death
# time of the interval decreases as we add noise. This is the result of the diameter of the
# hole in our circle shrinking.

# Finally, let's stretch our circle.

anim = @animate for r in vcat(0.0:0.02:2, 2:-0.02:0.0)
    points = noisy_circle(100; noise=0.1, r1=r)
    result = ripserer(points)

    plt_pts = scatter(
        points;
        legend=false,
        aspect_ratio=1,
        xlim=(-2.2, 2.2),
        ylim=(-2.2, 2.2),
        title="Data",
    )
    plt_diag = plot(result; infinity=3)

    plot(plt_pts, plt_diag; size=(800, 400))
end
gif(anim, "stability_anim_4.gif") # hide

# Again, we see the persistent homology stays stable, as long as the data at least somewhat
# resembles a circle.


text:
using ChaosTools
using Test
using DynamicalSystemsBase, DelayEmbeddings
using StatsBase
standardize = DelayEmbeddings.standardize
test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

ti = time()
@testset "ChaosTools tests" begin

include("basins/basins_tests.jl")
include("basins/uncertainty_tests.jl")
include("basins/tipping_points_tests.jl")

include("orbitdiagrams/orbitdiagram_tests.jl")
include("orbitdiagrams/poincare_tests.jl")

include("chaosdetection/lyapunov_exponents.jl")
#include("chaosdetection/gali_tests.jl") # TODO: make those faster
include("chaosdetection/partially_predictable_tests.jl")
include("chaosdetection/01test.jl")
#include("chaosdetection/expansionentropy_tests.jl")  # TODO: make those faster

include("period_return/periodicity_tests.jl")
include("period_return/period_tests.jl")
include("period_return/yin_tests.jl")
# include("period_return/transit_time_tests.jl")

include("dimensions/dims.jl")
include("dimensions/correlationdim.jl")
include("nlts_tests.jl")
# include("dyca_tests.jl") # TODO: fix DyCa tests once and for
# all by using https://github.com/JuliaDynamics/JuliaDynamics/tree/master/timeseries

end

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits = 33), " seconds or ",
round(ti/60, digits = 3), " minutes")


text:
include("../src/PowerFlowAnalysis.jl")

using JuMP
using Ipopt
using PowerModels
using PowerFlowAnalysis
using MATLAB
using PyPlot

casename = "case9" #Matpower test case
case = PowerModels.parse_file(join(["data/", casename, ".m"]))
pm = ACPPowerModel(case)
sol = run_ac_pf(case, IpoptSolver())
sys = PowerFlowSystem(pm, sol)
lf = LurieFormAd(sys)

# Choose between different choices of the u - constraints
net = Network(pm)
pqload = sys.p⁰[sys.pq]
nzload = find(pqload -> pqload < 0, pqload)
buspJ = nzload[[1,2]]##the bus ordering in this code
busp = [net.buses[k] for k in sys.pq[buspJ]] #the bus ordering in Matpower cases

npv = length(sys.pv)
nb = sys.nb

ratio = 1; #aspect ratio of two loads
ϕ = acot(ratio)
v⁰ = lf.sys.v⁰
R = create_loadability_R(lf, buspJ, [cos(ϕ), sin(ϕ)]) #loading scenario

lc = BaseLurieAdCertificate(lf, R)

Qlim = get_Qlims(lc, pm)
dV = 0.01 #voltage derivation limit
thermalrate = 2.
T = 0.01#angle separation limit
K = 0.01#logarithmic voltage difference limit

limitQQ = 1 #if reactive power constraints are considered
limitIQ = 1 #if reactive thermal constraints are considered

#core optimization problem for construction
mLP = brouwer_certificate_LP(lc, lc.B, lc.C, Qlim, thermalrate, T, K, dV, 100, limitQQ, limitIQ; region=:uniform)
sol_LP = solve(mLP)

Δ₁u = getvalue(getvariable(mLP, :Δ₁u))
Δ₁v = getvalue(getvariable(mLP, :Δ₁v))

#options for Continuation Power Flow
dV = findmax(Δ₁v[1,:])[1]
stepsize = 0.03
adaptive = 0

basep = -sys.p⁰[sys.pq[buspJ]] #the base point
resB = cal_box(R, Δ₁u, Δ₁v, npv, buspJ, basep, v⁰) #inner approx. region
nump = 4 #number of points to plot

polytope_plot(casename,busp,nump,resB,basep,nb,ratio,thermalrate,dV,stepsize,adaptive,limitQQ,limitIQ,Δ₁u)

max_stars_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_issues_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_forks_repo_name: hunghtd/PowerFlowAnalysis-tutorial


text:
using BrainFlow

# specify logging library to use
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)
sleep(5)
BrainFlow.stop_stream(board_shim)
data = BrainFlow.get_current_board_data(256, board_shim)
BrainFlow.release_session(board_shim)


text:
@testset "Utils" begin

    A = [missing 1.0; 2.0 missing]
    ClimateTools.replace_missing!(A)
    @test isnan(A[1,1])
    @test isnan(A[2,2])

end


text:
# We force specialisation on each function to avoid (tiny) allocations.
#
# Note that, for mapreduce, we can assume that the operation is commutative,
# which allows MPI to freely reorder operations.
#
# We also define mapfoldl (and mapfoldr) for completeness, even though the global
# operations are not strictly performed from left to right (or from right to
# left), since each process locally reduces first.
for (func, commutative) in [:mapreduce => true, :mapfoldl => false, :mapfoldr => false]
    @eval function Base.$func(f::F, op::OP, u::PencilArray; kws...) where {F, OP}
        rlocal = $func(f, op, parent(u); kws...)
        op_mpi = MPI.Op(op, typeof(rlocal); iscommutative = $commutative)
        MPI.Allreduce(rlocal, op_mpi, get_comm(u))
    end
end


text:
@testset "▶ types/figure                " begin
    f = Figure()
    io = IOBuffer()
    Base.show(io, MIME("text/plain"), f)
    s = String(take!(io))
    @test all(split(s) .== split(raw"""
    GPlot.Figure{GLE}
        Name:          default ("_fig_")
        Size:          (12.0, 9.0)
        Bg. color:     white
        N. axes:       0
        LaTeX:         false
        Transparent:   false"""))
end

@testset "▶ /figure                     " begin
    # CONSTRUCTORS
    # --> nothing
    f  = Figure()
    @test f.id == "_fig_"
    @test f.g isa GPlot.Backend
    @test f.axes isa Vector
    @test eltype(f.axes) == GPlot.Axes{GPlot.GLE}
    @test f.size == (12., 9.)
    @test f.textstyle.hei == 0.35 # default
    @test f.textstyle.font == "texcmss"
    @test isnothing(f.textstyle.color)
    @test isnothing(f.texlabels)
    @test isnothing(f.texscale)
    @test isnothing(f.texpreamble)
    @test isnothing(f.transparency)
    @test gcf() == f
    @test gca() === nothing

    # --> name
    f2 = Figure("test")
    @test f2.id == "test"
    @test gcf() == f2
    @test gca() === nothing
    @test GPlot.GP_ENV["ALLFIGS"]["test"] == f2

    # --> recuperate via id
    f3 = Figure("test")
    @test f2 === f3

    # add axes and check isempty
    @test isempty(f3.axes)
    GPlot.add_axes2d!()
    @test !isempty(f3.axes)

    GPlot.erase!(f3)
    @test isempty(f3.axes)

    GPlot.destroy(f3)
    @test "test" ∉ keys(GPlot.GP_ENV["ALLFIGS"])

    # Subroutine
    f = Figure()
    m = G.MarkerStyle("circle", 0.5, colorant"blue")
    G.add_sub_marker!(f, m)
    f.subroutines["mk_circle_rgba_0_0_0_0_1_0_1_0_"] == raw"""
    sub _mk_circle_rgba_0_0_0_0_1_0_1_0_ size mdata
        gsave
        set color rgba(0.0,0.0,1.0,1.0)
        mk circle 1
        grestore
    end sub
    define mk mk_circle_rgba_0_0_0_0_1_0_1_0_ _mk_circle_rgba_0_0_0_0_1_0_1_0_
    """
end

@testset "▶ set_prop/figure             " begin
    f = Figure(size=(5, 7), tex=true, texscale="fixed",
               alpha=true, transparency=true,
               preamble=tex"\usepackage{amssymb}",
               font="helvetica", fontsize=11)

    @test f.size == (5., 7.)
    @test f.texlabels == true
    @test f.texscale == "fixed"
    @test f.transparency == true
    @test f.texpreamble == t"\usepackage{amssymb}"
    @test f.textstyle.font == "psh"
    @test f.textstyle.hei  == 11 * GPlot.PT_TO_CM

    # when things go wrong, specifics for figure
    @test_throws MethodError Figure(tex=0)
    @test_throws GPlot.UnknownOptionError Figure(something=5)
    @test_throws GPlot.OptionValueError Figure(texscale="not-fixed")
end

@testset "▶ apply_gle/figure            " begin
    f = Figure(size=(10, 6), bgcol=c"red", bgalpha=0.1)
    s = G.assemble_figure(f; debug=true)

    isin(s, "size 10.0 6.0")
    isin(s, "amove -0.05 -0.05")
    # note the float rounding error, it's fine though
    isin(s, "box 10.1 6.1 fill rgba(1.0,0.0,0.0,0.102) nobox")
    isin(s, "set font texcmss hei 0.35") # default
end


text:
function PREOS(mixture, thermo)

mixing_rule_num = thermo.mixingrule;
activityfun = thermo.activity_model;
phase1 = thermo.phase;
fug_need = thermo.fugacity_switch;
critical_pres = [mixture.components.Pc]; %[Pa]
critical_temp = [mixture.components.Tc]; %[K]
acentric_fact = [mixture.components.acentric_factor]; %[-]
BIP = mixture.bip;
x = mixture.mole_fraction;
p = mixture.pressure; %[Pa]
T = mixture.temperature;
N = length(critical_temp);
fugacity = zeros(1,N);
R=8.314;
s1 = 0.623225;  %Huron Vidal

bi = 0.077796*R*critical_temp./critical_pres;
aci=0.457235*(R*critical_temp).^2 ./critical_pres;
mi = 0.37646+(1.54226-0.26992*acentric_fact).*acentric_fact;
Tr = T./critical_temp;
alfai = 1+mi.*(1-sqrt(Tr));   %//alfai=ai^0.5
alfa = alfai.^2.0;           %//alfa=ai
ai = aci .* alfa;
%Q is the parameter for MHV1 and MHV2 mixing rule
%it depends on the EOS
Q = (mixing_rule_num==3)*[-0.53 0]+(mixing_rule_num==4)*[-0.4347 -0.003654];
[a, b] = mixing_rule(mixture, thermo, ai, bi, s1, Q);

A_coef=a*p/(R*T)^2;
B_coef=b*p/(R*T);
%        poly_coef(4)=-B_coef*(A_coef-B_coef*(1+B_coef));
%        poly_coef(3)=A_coef-B_coef*(2.0+3*B_coef);
%        poly_coef(2)=-1+B_coef;
%        poly_coef(1)=1;
poly_coef = [1 -1+B_coef A_coef-B_coef*(2.0+3*B_coef) ...
    -B_coef*(A_coef-B_coef*(1+B_coef))];

z_root = roots(poly_coef);
if (sum(imag(z_root)~=0)==0)
    liquid_z = min(z_root);
    vapor_z = max(z_root);
else
    liquid_z = z_root(imag(z_root)==0);
    vapor_z = liquid_z;
end
% %{now we should calculate vapor and liquid
% %compressibility factor by the following method :}
% real_part = real(z_root);
% imag_part = imag(z_root);
% % findmaxmin(real_part,img_part,zv,zl);
% max_no=real_part(3);
% min_no=real_part(3);
% for i=1:3
%    if ((real_part(i)>max_no) && (real_part(i)>0) && (imag_part(i)==0))
%        max_no=real_part(i);
%    end
%    if ((real_part(i)<min_no) && (real_part(i)>0) && (imag_part(i)==0))
%          min_no=real_part(i);
%    end
% end
%      liquid_z = min_no;
%      vapor_z = max_no;
if (phase1==1)     %then  //phase 1 is liquid
   zz=liquid_z;
elseif (phase1==2) %then    //phase 2 is vapor
   zz=vapor_z;
end

if (fug_need==1)
    if (mixing_rule_num==1) %simple mixing rule
        part1=bi/b*(zz-1)-log(zz-b*p/(R*T));
        part2=x*(sqrt(ai'*ai).*(1-[BIP.EOScons]-[BIP.EOStdep]*T))';
        part3=A_coef/(2.828*B_coef)*(bi/b-2/a*part2) ...
            *log((zz+2.414*b*p/(R*T))/(zz-0.414*b*p/(R*T)));
        fugacity=exp(part1+part3);
    elseif (mixing_rule_num==2)  %Huron Vidal mixing rule
        [~, gama] = activityfun(T, x, component, BIP);
        part1 = bi/b*(zz-1)-log(zz-b*p/(R*T));
        part3 = -1/(2*sqrt(2))*(ai./bi/R/T - ...
                   log(gama)/0.623225)* ...
                   log((zz+(1+sqrt(2))*B_coef)/(zz+(1-sqrt(2))*B_coef));
        fugacity = exp(part1+part3);
    elseif (mixing_rule_num==3) % MHV1 mixing rule
        [~, gama] = activityfun(T, x, component, BIP);
        q1 = -0.53;  %Michelsen for PR
        logfi = bi/b*(zz-1) - log(zz-B_coef) - 1/(2*sqrt(2))*(ai./(bi*R*T) ...
            + log(gama)/q1+log(b./bi)/q1+(bi/b-1)/q1)*log((zz+(1+sqrt(2))*B_coef)/ ...
            (zz+(1-sqrt(2))*B_coef));
        fugacity = exp(logfi);
    elseif (mixing_rule_num==4) % MHV2 mixing rule
        [~, gama] = activityfun(T, x, component, BIP);
        q1 = -0.4347;
        q2 = -0.003654;
        alphai = ai./(bi*R*T);
        alpha = a/(b*R*T);
        logfi = bi/b*(zz-1) - log(zz-B_coef) - ...
            1/(2*sqrt(2))*(q1*alphai+q2*(alpha^2+alphai.^2)+log(gama) ...
            +log(b./bi)+bi/b-1)/(q1+2*q2*alpha)* ...
            log((zz+(1+sqrt(2))*B_coef)/(zz+(1-sqrt(2))*B_coef));
        fugacity = exp(logfi);
    end
end

%% calculate the residual enthalpy
for i=1:N
   Abar(i)=0;
   for j=1:N
     Abar(i)=Abar(i)+sqrt(ai(j))/(R*T)*x(j);
   end
   Abar(i)=Abar(i)*sqrt(ai(i))/(R*T);
end
 part1=0;
 for i=1:N
   part1=part1+x(i)*Abar(i)*mi(i)*(-1)/sqrt(critical_temp(i))/alfai(i);
 end
 part1=(part1*T^0.5/(A_coef/p)-1)*R*T*A_coef/(B_coef*2*sqrt(2))*log((zz+B_coef*(1+sqrt(2)))/(zz+B_coef*(1-sqrt(2))));
 HR=part1+R*T*(zz-1);
% HR=0;


text:
using MacroTools

walk(x, inner, outer) = outer(x)
walk(x::Expr, inner, outer) = outer(Expr(x.head, map(inner, x.args)...))
postwalk(f, x) = walk(x, x -> postwalk(f, x), f)

function referred(expr::Expr)
    if expr.head == :$
        :($(Expr(:$, :(Ref($(expr.args...)))))[])
    else
        expr
    end
end
referred(x)  = x

"""
    @noelide _bmacro_ expression
where _bmacro_ is one of @btime, @belapsed, @benchmark
Wraps all interpolated code in _expression_ in a __Ref()__ to
stop the compiler from cheating at simple benchmarks. Works
with any macro that accepts interpolation
#Example
    julia> @btime \$a + \$b
      0.024 ns (0 allocations: 0 bytes)
    3
    julia> @noelide @btime \$a + \$b
      1.277 ns (0 allocations: 0 bytes)
    3
"""
macro noelide(expr)
    out = postwalk(referred, expr) |> esc
end


text:
"""
    SeisBinData(in,out; <keyword arguments>)

Sequentially bin seismic data using already binned trace headers (SeisBinHeaders).
Input arguments should be consistent with SeisBinHeaders input arguments.

# Arguments
* `in`: filename of input, irregularly sampled data
* `out`: filename of output, regularly sampled data

# Keyword arguments
* `style="sxsygxgy"`: bin style. Options: "mxmyhxhy","mxmyhaz","sxsyhxhy","gxgyhxhy","sxsyhaz","gxgyhaz"
* `ang=90`: inline direction measured in degrees CC from East
* `gamma=1`: vp/vs ratio for PS Asymptotic Conversion Point gathers (use gamma=1 for PP data)
* `osx=0`,`osy=0`,`ogx=0`,`ogy=0` : origin for source and receiver coordinate system
* `omx=0`,`omy=0`,`ohx=0`,`ohy=0`: origin for midpoint and offset coordinate system
* `oaz=0`,`oh=0` : origin for azimuth and offset coordinate system
* `dsx=1`,`dsy=1`,`dgx=1`,`dgy=1`: source and receiver step-size
* `dmx=1`,`dmy=1`,`dhx=1`,`dhy=1`: midpoint and offset step-size
* `dh=1`,`daz=1`: offset and azimuth step-size
* `min_isx=0`,`max_isx=0`,`min_isy=0`,`max_isy=0`: grid extreme values for sources
* `min_igx=0`,`max_igx=0`,`min_igy=0`,`max_igy=0`: grid extreme values for receivers
* `min_imx=0`,`max_imx=0`,`min_imy=0`,`max_imy=0`: grid extreme values for midpoints
* `min_ihx=0`,`max_ihx=0`,`min_ihy=0`,`max_ihy=0`: grid extreme values for offsets
* `min_ih=0`,`max_ih=0`,`min_iaz=0`,`max_iaz=0`: grid extreme values for azimuth and offset
* `ntrace=10000`: maximum number of traces processed at a time

# Output
In file `out`, the binned data is created.

*Credits: Aaron Stanton, Fernanda Carozzi, 2017*

"""

function SeisBinData(in,out;style="sxsygxgy",ang=90,gamma=1,osx=0,osy=0,ogx=0,ogy=0,omx=0,omy=0,ohx=0,ohy=0,oh=0,oaz=0,dsx=1,dsy=1,dgx=1,dgy=1,dmx=1,dmy=1,dhx=1,dhy=1,dh=1,daz=1,min_isx=0,max_isx=0,min_isy=0,max_isy=0,min_igx=0,max_igx=0,min_igy=0,max_igy=0,min_imx=0,max_imx=0,min_imy=0,max_imy=0,min_ihx=0,max_ihx=0,min_ihy=0,max_ihy=0,min_ih=0,max_ih=0,min_iaz=0,max_iaz=0,ntrace=10000)

	rad2deg = 180/pi;
	deg2rad = pi/180;
	gammainv = 1/gamma;
	if (ang > 90)
		ang2=-deg2rad*(ang-90)
	else
		ang2=deg2rad*(90-ang)
	end

	naz=convert(Int32,360/daz)

	if (style=="sxsygxgy")
		nsx = max_isx - min_isx + 1
		nsy = max_isy - min_isy + 1
		ngx = max_igx - min_igx + 1
		ngy = max_igy - min_igy + 1

		nx1=nsx;nx2=nsy;nx3=ngx;nx4=ngy;

	elseif (style=="mxmyhxhy")
		nmx = max_imx - min_imx + 1
		nmy = max_imy - min_imy + 1
		nhx = max_ihx - min_ihx + 1
		nhy = max_ihy - min_ihy + 1

		nx1=nmx;nx2=nmy;nx3=nhx;nx4=nhy;

	elseif (style=="mxmyhaz")
	  nmx = max_imx - min_imx + 1
		nmy = max_imy - min_imy + 1
		nh  = max_ih  - min_ih  + 1
		naz = max_iaz - min_iaz + 1

		nx1=nmx;nx2=nmy;nx3=nh;nx4=naz;

	elseif (style=="sxsyhxhy")
	  nsx = max_isx - min_isx + 1
	  nsy = max_isy - min_isy + 1
		nhx = max_ihx - min_ihx + 1
		nhy = max_ihy - min_ihy + 1

		nx1=nsx;nx2=nsy;nx3=nhx;nx4=nhy;

	elseif (style=="gxgyhxhy")
		ngx = max_igx - min_igx + 1
		ngy = max_igy - min_igy + 1
		nhx = max_ihx - min_ihx + 1
		nhy = max_ihy - min_ihy + 1

		nx1=ngx;nx2=ngy;nx3=nhx;nx4=nhy;

	elseif (style=="sxsyhaz")
	  nsx = max_isx - min_isx + 1
	 	nsy = max_isy - min_isy + 1
		nh  = max_ih  - min_ih  + 1
		naz = max_iaz - min_iaz + 1

		nx1=nsx;nx2=nsy;nx3=nh;nx4=naz;

	elseif (style=="gxgyhaz")
    ngx = max_igx - min_igx + 1
		ngy = max_igy - min_igy + 1
		nh  = max_ih  - min_ih  + 1
		naz = max_iaz - min_iaz + 1

		nx1=ngx;nx2=ngy;nx3=nh;nx4=naz;

	else
		error("style not recognized.")
	end

	nx_out = nx1*nx2*nx3*nx4

	stream = open(ParseHeaderName(in))

	@compat nx_in = round(Int,filesize(stream)/(4*length(fieldnames(Header))))

	seek(stream, header_count["n1"])
	nt = read(stream,Int32)

	seek(stream, header_count["o1"])
	o1 = read(stream,Float32)

	seek(stream, header_count["d1"])
	dt = read(stream,Float32)

	close(stream)

	d = zeros(Float32,nt,1)

	out_d = ParseDataName(out)
	out_h = ParseHeaderName(out)
	stream_d = open(out_d,"a")

#Create 0 volume
	j = 1
	for ix4 = 1 : nx4
			for ix3 = 1 : nx3
				 for ix2 = 1 : nx2
					 for ix1 = 1 : nx1

		    		position_d = 4*nt*(j - 1)
		    		seek(stream_d,position_d)
						write(stream_d,d[:,1])
						j += 1
		    	end
				end
  	  end
	end
	close(stream_d)

	stream_d = open(out_d,"a")

	if (style=="sxsygxgy")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = ntrace
			for k = 1 : num_traces_in
				itrace = (h[k].igy - min_igy)*nx1*nx2*nx3 + (h[k].igx - min_igx)*nx2*nx1 + (h[k].isy - min_isy)*nx1 + h[k].igx - min_igx + 1
				if (itrace > 0 && itrace <= nx_out)

					if (h[k].isx >= min_isx && h[k].isx <= max_isx && h[k].isy >= min_isy && h[k].isy <= max_isy && h[k].igx >= min_igx && h[k].igx <= max_igx && h[k].igy >= min_igy && h[k].igy <= max_igy)

						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])

				  end
				end
			end
			j += num_traces_in
		end

	elseif (style=="mxmyhxhy")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)

			for k = 1 : num_traces_in
				itrace = (h[k].ihy - min_ihy)*nx3*nx2*nx1 + (h[k].ihx - min_ihx)*nx2*nx1 + (h[k].imy - min_imy)*nx1 + h[k].imx - min_imx + 1
	 		 if (itrace > 0 && itrace <= nx_out)

	 				 if (h[k].imx >= min_imx && h[k].imx <= max_imx && h[k].imy >= min_imy && h[k].imy <= max_imy && h[k].ihx >= min_ihx && h[k].ihx <= max_ihx && h[k].ihy >= min_ihy && h[k].ihy <= max_ihy)

						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])

					end
				end
			end
			j += num_traces_in
		end

	elseif (style=="mxmyhaz")
		j = 1
		while j <= nx_in
			d,h_in,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)

			for k = 1 : num_traces_in

				itrace = (h_in[k].iaz - min_iaz)*nx3*nx2*nx1 + (h_in[k].ih - min_ih)*nx2*nx1 + (h_in[k].imy - min_imy)*nx1 + h_in[k].imx - min_imx + 1
			  if (itrace > 0 && itrace <= nx_out)

			 		 if (h_in[k].imx >= min_imx && h_in[k].imx <= max_imx && h_in[k].imy >= min_imy && h_in[k].imy <= max_imy && h_in[k].ih >= min_ih && h_in[k].ih <= max_ih && h_in[k].iaz >= min_iaz && h_in[k].iaz <= max_iaz)

						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])

					end
				end
			end
			j += num_traces_in
		end

	elseif (style=="sxsyhxhy")
		j = 1
		while j <= nx_in
			d,h_in,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)

			for k = 1 : num_traces_in
				itrace = (h_in[k].ihy - min_ihy)*nx3*nx2*nx1 + (h_in[k].ihx - min_ihx)*nx2*nx1 + (h_in[k].isy - min_isy)*nx1 + h_in[k].isx - min_isx + 1
				if (itrace > 0 && itrace <= nx_out)

					 if (h_in[k].isx >= min_isx && h_in[k].isx <= max_isx && h_in[k].isy >= min_isy && h_in[k].isy <= max_isy && h_in[k].ihx >= min_ihx && h_in[k].ihx <= max_ihx && h_in[k].ihy >= min_ihy && h_in[k].ihy <= max_ihy)

						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
					end
				end
			end
			j += num_traces_in
		end

	elseif (style=="gxgyhxhy")
		j = 1
		while j <= nx_in
			d,h_in,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h_in[k].ihy - min_ihy)*nx3*nx2*nx1 + (h_in[k].ihx - min_ihx)*nx2*nx1 + (h_in[k].igy - min_igy)*nx1 + h_in[k].igx - min_igx + 1
				if (itrace > 0 && itrace <= nx_out)

					 if (h_in[k].igx >= min_igx && h_in[k].igx <= max_igx && h_in[k].igy >= min_igy && h_in[k].igy <= max_igy && h_in[k].ihx >= min_ihx && h_in[k].ihx <= max_ihx && h_in[k].ihy >= min_ihy && h_in[k].ihy <= max_ihy)

						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])

					end
				end
			end
			j += num_traces_in
		end

	elseif (style=="sxsyhaz")
		j = 1
		while j <= nx_in
			d,h_in,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in
				itrace = (h_in[k].iaz - min_iaz)*nx3*nx2*nx1 + (h_in[k].ih - min_ih)*nx2*nx1 + (h_in[k].isy - min_isy)*nx1 + h_in[k].isx - min_isx + 1
				if (itrace > 0 && itrace <= nx_out)

					 if (h_in[k].isx >= min_isx && h_in[k].isx <= max_isx && h_in[k].isy >= min_isy && h_in[k].isy <= max_isy && h_in[k].ih >= min_ih && h_in[k].ih <= max_ih && h_in[k].ias >= min_iaz && h_in[k].iaz <= max_iaz)

						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])
					end
				end
			end
			j += num_traces_in
		end

	elseif (style=="gxgyhaz")
		j = 1
		while j <= nx_in
			d,h,e = SeisRead(in,group="some",itrace=j,ntrace=ntrace)
			num_traces_in = size(d[:,:],2)
			for k = 1 : num_traces_in

				itrace = (h[k].iaz - min_iaz)*nx1*nx2*nx3 + (h[k].ih - min_ih)*nx1*nx2 + (h[k].igy - min_igy)*nx1 + h[k].igx - min_igx + 1
				if (itrace > 0 && itrace <= nx_out)

					if (h[k].igx >= min_igx && h[k].igx <= max_igx && h[k].igy >= min_igy && h[k].igy <= max_igy && h[k].ih >= min_ih && h[k].ih <= max_ih && h[k].iaz >= min_iaz && h[k].iaz <= max_iaz)

						position_d = 4*nt*(itrace - 1)
						seek(stream_d,position_d)
						write(stream_d,d[:,k])

					end
				end
			end
			j += num_traces_in
		end
	end
	close(stream_d)

end


text:
using StochasticSearch
using FactCheck

include("parameter.jl")
include("configuration.jl")
include("run.jl")
include("simulated_annealing.jl")
include("iterative_first_improvement.jl")
include("iterative_probabilistic_improvement.jl")
include("randomized_first_improvement.jl")
include("iterative_greedy_construction.jl")
include("iterated_local_search.jl")
include("unit_value.jl")

FactCheck.exitstatus()


text:
# ---
# title: 1122. Relative Sort Array
# id: problem1122
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Sort
# link: <https://leetcode.com/problems/relative-sort-array/description/>
# hidden: true
# ---
#
# Given two arrays `arr1` and `arr2`, the elements of `arr2` are distinct, and
# all elements in `arr2` are also in `arr1`.
#
# Sort the elements of `arr1` such that the relative ordering of items in `arr1`
# are the same as in `arr2`.  Elements that don't appear in `arr2` should be
# placed at the end of `arr1` in **ascending** order.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr1 = [2,3,1,3,2,4,6,7,9,2,19], arr2 = [2,1,4,3,9,6]
#     Output: [2,2,2,1,4,3,3,9,6,7,19]
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr1.length, arr2.length <= 1000`
#   * `0 <= arr1[i], arr2[i] <= 1000`
#   * All the elements of `arr2` are **distinct**.
#   * Each `arr2[i]` is in `arr1`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ---
# title: 160. Intersection of Two Linked Lists
# id: problem160
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/intersection-of-two-linked-lists/description/>
# hidden: true
# ---
#
# Write a program to find the node at which the intersection of two singly
# linked lists begins.
#
# For example, the following two linked lists:
#
# [![](https://assets.leetcode.com/uploads/2018/12/13/160_statement.png)](https://assets.leetcode.com/uploads/2018/12/13/160_statement.png)
#
# begin to intersect at node c1.
#
#
#
# **Example 1:**
#
# [![](https://assets.leetcode.com/uploads/2020/06/29/160_example_1_1.png)](https://assets.leetcode.com/uploads/2020/06/29/160_example_1_1.png)
#
#
#
#     Input: intersectVal = 8, listA = [4,1,8,4,5], listB = [5,6,1,8,4,5], skipA = 2, skipB = 3
#     Output: Reference of the node with value = 8
#     **Input Explanation:** The intersected node's value is 8 (note that this must not be 0 if the two lists intersect). From the head of A, it reads as [4,1,8,4,5]. From the head of B, it reads as [5,6,1,8,4,5]. There are 2 nodes before the intersected node in A; There are 3 nodes before the intersected node in B.
#
#
#
# **Example 2:**
#
# [![](https://assets.leetcode.com/uploads/2020/06/29/160_example_2.png)](https://assets.leetcode.com/uploads/2020/06/29/160_example_2.png)
#
#
#
#     Input: intersectVal = 2, listA = [1,9,1,2,4], listB = [3,2,4], skipA = 3, skipB = 1
#     Output: Reference of the node with value = 2
#     **Input Explanation:**  The intersected node's value is 2 (note that this must not be 0 if the two lists intersect). From the head of A, it reads as [1,9,1,2,4]. From the head of B, it reads as [3,2,4]. There are 3 nodes before the intersected node in A; There are 1 node before the intersected node in B.
#
#
#
#
# **Example 3:**
#
# [![](https://assets.leetcode.com/uploads/2018/12/13/160_example_3.png)](https://assets.leetcode.com/uploads/2018/12/13/160_example_3.png)
#
#
#
#     Input: intersectVal = 0, listA = [2,6,4], listB = [1,5], skipA = 3, skipB = 2
#     Output: null
#     **Input Explanation:** From the head of A, it reads as [2,6,4]. From the head of B, it reads as [1,5]. Since the two lists do not intersect, intersectVal must be 0, while skipA and skipB can be arbitrary values.
#     Explanation: The two lists do not intersect, so return null.
#
#
#
#
# **Notes:**
#
#   * If the two linked lists have no intersection at all, return `null`.
#   * The linked lists must retain their original structure after the function returns.
#   * You may assume there are no cycles anywhere in the entire linked structure.
#   * Each value on each linked list is in the range `[1, 10^9]`.
#   * Your code should preferably run in O(n) time and use only O(1) memory.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Revise, WaltonForaging
##
init_forage_time = [1,0.5,0.227]
reward = 10
θ = 5
e1 = Exponential(θ)
mean(e1)
rand(e1)
quantile(e1,0.025)
quantile(e1,0.925)

function update_time(rew_num, initial)
    M = initial*(1.3^(rew_num-1))
    EXP = Exponential(M)
    CIlow = quantile(EXP,0.025)
    CIhigh =  quantile(EXP,0.925)
    return M, CIlow, CIhigh
end

function update_time(initial; n_rewards = 15)
    df = DataFrame(Mean = Float64[],CI_L = Float64[], CI_H = Float64[], Rew = Float64[])
    prev = 0
    for i in 1:n_rewards
        m, cil, cih = update_time(i,initial)
        push!(df,(Mean = m+prev,CI_L = cil+prev, CI_H = cih+prev, Rew = i))
        prev +=m
    end
    return df
end
##
poor = update_time(1)
medium = update_time(0.5)
rich = update_time(0.227)

@df poor plot(:Rew,:Mean, ribbon = (:Mean - :CI_L,:CI_H - :Mean), fillalpha = 0.2, label = "poor")
@df medium plot!(:Rew,:Mean, ribbon = (:Mean - :CI_L,:CI_H - :Mean), fillalpha = 0.2, label = "medium")
@df rich plot!(:Rew,:Mean, ribbon = (:Mean - :CI_L,:CI_H - :Mean), fillalpha = 0.2, label = "high", legend = :topleft)
##
@df poor plot(:Mean, :Rew, ribbon = (:CI_L,:CI_H), fillalpha = 0.2,label = "poor")
@df medium plot!(:Mean, :Rew, ribbon = (:CI_L,:CI_H), fillalpha = 0.2, label = "medium")
@df rich plot!(:Mean, :Rew, ribbon = (:CI_L,:CI_H), fillalpha = 0.2, label = "high",
     legend = :topleft, xlims = (0,15), ylims=(1,18),xticks = 0:1:15,yticks = 0:1:18)
##
xaxis = 0:20
low = [0.075,32.5]
high = [0.075,57.5]
plot(xaxis,basic(xaxis,low),label = "low", color = :grey)
plot!(xaxis,basic(xaxis,high),label = "high", color = :black)
low2 = [0.075*1.1,32.5]
high2 = [0.075*1.1,57.5]
plot!(xaxis,basic(xaxis,low2),label = "low2", color = :grey, linestyle = :dash)
plot!(xaxis,basic(xaxis,high2),label = "high2", color = :black, linestyle = :dash)
##
inverse(20,low2) - inverse(20,low)
inverse(20,high2) - inverse(20,high)
vline!([inverse(20,high)], color = :red)
vline!([inverse(20,high2)], color = :red, linestyle = :dash)
vline!([inverse(20,low)], color = :green)
vline!([inverse(20,low2)], color = :green, linestyle = :dash)
hline!([20], color = :blue)
##
xaxis = 0:20
low = [0.075,32.5]
high = [0.075,57.5]
plot(xaxis,basic(xaxis,low),label = "low", color = :grey, legend = :topleft)
plot!(xaxis,basic(xaxis,high),label = "high", color = :black)
hline!([20], color = :blue, label = "threshold")
scaling_factor = 1.1
low2 = [0.075*scaling_factor,32.5]
high2 = [0.075*scaling_factor,57.5]
plot!(xaxis,basic(xaxis,low2),label = "scaled_low", color = :grey, linestyle = :dash)
plot!(xaxis,basic(xaxis,high2),label = "scaled_high", color = :black, linestyle = :dash)
vline!([inverse(20,high)], color = :red, label = "")
vline!([inverse(20,high2)], color = :red, linestyle = :dash, label = "")
vline!([inverse(20,low)], color = :green, label = "")
vline!([inverse(20,low2)], color = :green, linestyle = :dash, label = "")
##
inverse(20,low2) - inverse(20,low)
inverse(20,high2) - inverse(20,high)
##
xaxis = 0:20
low = [0.075,32.5]
high = [0.075,57.5]
plot(xaxis,basic(xaxis,low),label = "low", color = :grey, legend = :topleft)
plot!(xaxis,basic(xaxis,high),label = "high", color = :black)
hline!([20], color = :blue, label = "threshold")
scaling_factor = 1.1
hline!([20*scaling_factor], color = :blue, linestyle = :dash, label = "scaled_treshold")
vline!([inverse(20,high)], color = :red, label = "")
vline!([inverse(20*scaling_factor,high)], color = :red, linestyle = :dash, label = "")
vline!([inverse(20,low)], color = :green, label = "")
vline!([inverse(20*scaling_factor,low)], color = :green, linestyle = :dash, label = "")
##
inverse(20*scaling_factor,low) - inverse(20,low)
inverse(20*scaling_factor,high) - inverse(20,high)


text:
function compare_with_experimets(base_model)
    for stst in RathData.STSTS
        model = deepcopy(base_model)
        println("\nStst: $stst")
        ξ = RathData.val(:ξ, stst)
        println("exp xi: $ξ")
        exp_μ = RathData.val(:μ, stst)
        println("exp growth: $exp_μ")
        base_intake_info = Human1.HumanGEM.stst_base_intake_info(stst)
        Chemostat.apply_bound!(model, ξ, base_intake_info;
            emptyfirst = true, ignore_miss = true
        )
        fbaout_μ = MetLP.objval(MetLP.fba!(model, Human1.HumanGEM.HUMAN_BIOMASS_IDER))
        println("fba growth: $(fbaout_μ)")
        fbaout_μ < exp_μ && @warn("fbaout_μ ($fbaout_μ) > exp_μ ($exp_μ)")
    end
end


text:
# Simple propositional logic

module PropositionalLogic

using ..LogicBase

export PropositionalVariable, BinaryConnective, assignments, ⊢

abstract Formula

immutable PropositionalVariable <: Formula
    var::Symbol
end

immutable Negate <: Formula
    f::Formula
end

immutable BinaryConnective{T} <: Formula
    x::Formula
    y::Formula
end

function Base.writemime(io::IO, ::MIME"text/plain", pv::PropositionalVariable)
    print(io, pv.var)
end

function Base.writemime(io::IO, ::MIME"text/plain", n::Negate)
    print(io, '¬')
    writemime(io, "text/plain", n.f)
end

function Base.writemime{T}(io::IO, ::MIME"text/plain", bc::BinaryConnective{T})
    print(io, '(')
    writemime(io, "text/plain", bc.x)
    print(io, ' ', T, ' ')
    writemime(io, "text/plain", bc.y)
    print(io, ')')
end

immutable Assignment
    table::Dict{PropositionalVariable, Bool}
end

assign(a::Assignment, v::PropositionalVariable) = a.table[v]
assign(a::Assignment, n::Negate) = ¬assign(a, n.f)
assign(a::Assignment, c::BinaryConnective{:∧}) =
    assign(a, c.x) && assign(a, c.y)
assign(a::Assignment, c::BinaryConnective{:∨}) =
    assign(a, c.x) || assign(a, c.y)
assign(a::Assignment, c::BinaryConnective{:⇒}) = assign(a, c.x) ⇒ assign(a, c.y)
assign(a::Assignment, c::BinaryConnective{:⇔}) = assign(a, c.x) ⇔ assign(a, c.y)

immutable AssignmentIterator
    over::Vector{PropositionalVariable}
end

Base.start(ai::AssignmentIterator) = UInt(1) << length(ai.over)
function Base.next(ai::AssignmentIterator, state::UInt)
    state -= 1
    len = length(ai.over)
    assignmentdata = [
        v => state & (UInt(1) << (len - i)) ≠ 0
        for (i, v) in enumerate(ai.over)]
    Assignment(assignmentdata), state
end
Base.done(ai::AssignmentIterator, state::UInt) = state == UInt(0)

vars(v::PropositionalVariable) = Set([v])
vars(n::Negate) = vars(n.f)
vars(bc::BinaryConnective) = vars(bc.x) ∪ vars(bc.y)

assignments(vs::Vector{PropositionalVariable}) = AssignmentIterator(vs)
assignments(vs::Set{PropositionalVariable}) = assignments(collect(vs))
assignments(f::Formula) = assignments(vars(f))
assignments(f::Formula, g::Formula) = assignments(vars(f) ∪ vars(g))
assignments(vs::Vector{Formula}) = assignments(foldl(∪, Set(), vs))

#= Tautology & Entailment =#
⊢(::Tuple{}, f::Formula) = all(a -> assign(a, f), assignments(f))
⊢(f::Formula, g::Formula) =
    all(a -> assign(a, f) ⇒ assign(a, g), assignments(f, g))
⊢(fs, g::Formula) = all(a -> ¬assign(a, g) ∧ ¬⋀(map(f -> assign(a, f), fs)))

#= Convenience Functions =#
LogicBase. ¬(f::Formula) = Negate(f)
LogicBase. ∧(f::Formula, g::Formula) = BinaryConnective{:∧}(f, g)
LogicBase. ∨(f::Formula, g::Formula) = BinaryConnective{:∨}(f, g)
LogicBase. ⊻(f::Formula, g::Formula) = BinaryConnective{:⊻}(f, g)
LogicBase. ⇒(f::Formula, g::Formula) = BinaryConnective{:⇒}(f, g)
LogicBase. ⇔(f::Formula, g::Formula) = BinaryConnective{:⇔}(f, g)

end

max_stars_repo_name: TotalVerb/FirstOrderLogic

max_issues_repo_name: TotalVerb/FirstOrderLogic

max_forks_repo_name: TotalVerb/FirstOrderLogic


text:
abstract type ATruth end

struct Axiom <: ATruth  # its id is the position in the truth table
    name::String
    left
    right
end
struct Definition <: ATruth
    name::String
    left
    right
end

struct Proof
    lines::Array{Tuple{Node,Int64,Array{Int64}}}
end
lines(x::Proof) = x.lines
steps(x::Proof) = [line[1] for line in x.lines]
function english(x::Proof)
    result = ""
    for line in lines(x)
        result *= "$(line[1]) because $(name(truths[line[2]]))\n"
    end
    return result
end
# Base.:(==)(x::Proof, y::Proof) = ((steps(x)[end] in steps(y)) || (steps(y)[end] in steps(x)))
Base.hash(x::Proof, h::UInt) = hash(steps(x)[end], h)
Base.:(==)(x::Proof, y::Proof) = steps(x)[end] == steps(y)[end]
extensions(x::Proof, truth_table=truths) = [Proof(push!(copy(lines(x)), re)) for re in rewrites_wp(steps(x)[end], Int64[], truth_table)]

struct Theorem <: ATruth
    name::String
    left
    right
    proof::Proof
end
name(x::ATruth) = x.name