text:
module Paisley

export bfbfb,
       catalog,
       cgdae,
       codex,
       dragon,
       eadgbe,
       erase!,
       fkbjdn,
       place!,
       sentinel,
       tacit

include("shelves.jl")

using .Shelves

function bounds(span::Int)
  35 < span < 61
end

function tacit()
  local mute = repeat("__ ", 12)
  bounds(strwidth(mute)) ? mute : "Check tacit"
end

function sentinel(yarn::String)
  local rex = r"^([jknz]+\d+)+([lm]\d+)?[o-z]*$"
  ismatch(rex, lowercase(yarn))
end

function erase!(yarn::String)
  local seal = symbol(lowercase(yarn))
  if haskey(codex, seal)
    delete!(codex, seal)
    println("\n\t$seal deleted")
  else
    println("\n\t$seal not found")
  end
end

function place!(yarn::String, wire::String)
  if sentinel(yarn) && is_valid_ascii(wire)
    local seal = symbol(lowercase(yarn))
    local span = strwidth(wire)
    if bounds(span)
      codex[seal] = wire
      return nothing
    elseif isless(span, 36)
      println("\n\tSize: $span? $wire")
    else
      println("\n\tSize: $span?")
    end
  elseif sentinel(yarn)
    println("\n\tASCII? $wire")
  else
    println("\n\t$yarn?")
  end
end

place!("z0", tacit()) # default

function catalog()
  local art = collect(keys(codex))
  sort!(art)
  local columns = 7
  let nth = 1
    println()
    for (nth, its) in enumerate(art)
      print("\t", its)
      if ==(nth % columns, 0)
        println()
      end
    end
    if !=(nth % columns, 0)
      println()
    end
  end
end

function transmute!(cord::String)
  let line = copy(cord)
    if trust
      for (old, new) in duets
        line = replace(line, old, new)
      end
    end
    return line
  end
end

function pitch(seal::Symbol, nth::Int)
  try
    local wire = get(codex, seal, tacit())
    if bounds(strwidth(wire))
      local cord = string(wire[nth:end], wire[1:nth + 1])
      transmute!(cord)
    else
      tacit()
    end
  catch anomaly
    println("\nCause $anomaly")
  end
end

function hermit(seal)
  pitch(seal, 1)
end

function jovian(seal)
  pitch(seal, 7)
end

function copper(seal)
  pitch(seal, 13)
end

function saturn(seal)
  pitch(seal, 16)
end

function gemini(seal)
  pitch(seal, 19)
end

function helios(seal)
  pitch(seal, 22)
end

function silver(seal)
  pitch(seal, 28)
end

function aquari(seal)
  pitch(seal, 31)
end

function vulcan(seal)
  pitch(seal, 34)
end

function bfbfb(seal)
  local fn = saturn(seal)
  local bn = vulcan(seal)
  (bn, fn, bn, fn, bn)
end

function cgdae(seal)
  copper(seal),
  silver(seal),
  jovian(seal),
  helios(seal),
  hermit(seal)
end

function eadgbe(seal)
  local ingot = copper(seal)
  ingot,
  vulcan(seal),
  helios(seal),
  jovian(seal),
  silver(seal),
  ingot
end

function fkbjdn(seal)
  local dn = jovian(seal)
  local fk = gemini(seal)
  local bj = aquari(seal)
  local rows = (fk, bj, dn, fk, bj, dn)
  reverse(rows)
end

function dragon(seal)
  aquari(seal),
  saturn(seal),
  hermit(seal),
  helios(seal),
  jovian(seal),
  silver(seal),
  copper(seal),
  vulcan(seal),
  gemini(seal)
end

end # module

max_stars_repo_name: reidiiius/masquerade

max_issues_repo_name: reidiiius/masquerade

max_forks_repo_name: reidiiius/masquerade


text:
#TODO(Corneliu): documentation
function translate(imodel,
                   omodel,
                   word::AbstractString;
                   n=10,
                   threshold=0.5,
                   vocab=nothing,
                   ofilter=x->true)
    if !in(word, imodel.vocab)
        return [word], false
    end
	pos, scores = cosine_vec(omodel, get_vector(imodel, word), n, vocab=vocab)
    pos = pos[scores .>= threshold]
    if isempty(pos)
        return [word], false
    else
        candidates = filter(ofilter, omodel.vocab[pos])
        if isempty(candidates)
            return [word], false
        else
            return candidates, true
        end
    end
end

function translate(imodel,
                   omodel,
                   words::Vector{T};
                   ngram_complexity=1,
                   separator="_",
                   kwargs...) where {T}
    if ngram_complexity == 1
        _words = words
    else
        _words = T[]
        pos = ngram_search(imodel.vocab, words; max_length=ngram_complexity)
        fpos = Iterators.flatten(pos)
        not_found_positions = setdiff(1:length(words), fpos)
        for i in eachindex(words)
            if i in not_found_positions
                # word does not exist in input languge
                push!(_words, words[i])
            else
                # Attempt to translate n-gram
                # if the word is the last part of it
                for p in pos
                    if i == maximum(p)
                        joined = join(words[p], separator, separator)
                        push!(_words, joined)
                    end
                end
            end
        end

    end
    return [translate(imodel, omodel, word; kwargs...) for word in _words]
end

function ngram_search(vocab, tokens; separator="_", max_length=3)
    # Initializations
    found = Vector{UnitRange{Int}}()
    n = length(tokens)
    i = 1
    j = n
    while i <= n
        if j-i+1 <= max_length
            token = join(tokens[i:j], separator, separator)
            is_match = token in vocab
            if is_match
                push!(found, i:j)
                i = j + 1
                j = n
                continue
            end
        end
        if i == j
            j = n
            i+= 1
        else
            j-= 1
        end
    end
    return found
end

ext: jl

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
using MetacommunityDynamics
using DynamicGrids
using Plots
using Distributions
using Dispersal: Moore, DispersalKernel
using EcologicalNetworks: mpnmodel, trophic_level, UnipartiteNetwork, degree, richness
using EcologicalNetworksPlots
using ColorSchemes

number_of_species = 300
connectance = 0.05
forbiddenlinkprob = 0.5

foodweb = mpnmodel(number_of_species, connectance, forbiddenlinkprob)
heatmap(Matrix(foodweb.edges), legend=:none, aspectratio=1)

function foodwebplot(N)
    I = initial(FoodwebInitialLayout, N)
    L = SpringElectric(1.2; gravity=0.05)
    L.move = (true, false)
    for step in 1:(100richness(N))
      position!(L, I, N)
    end
    plot(I, N, la=0.9, frame=:box, dpi=300)
    scatter!(I, N, nodefill=trophic_level(N), nodesize=degree(N), c=:YlGn)
end
foodwebplot(foodweb)

savefig("foodweb.png")

speciespool = DiscreteUnipartiteSpeciesPool(Symbol.(foodweb.S), Matrix(foodweb.edges)) # move these type changes to a method
trophicdict = trophic_level(foodweb)  # returns a dictionary

resource = filter(s -> trophicdict[String(s)] == 1, species(speciespool))
consumers = filter(s -> trophicdict[String(s)] > 1, species(speciespool))

consumermodel =
    FoodWebEating(consumers, resource, LotkaVolterra(0.2), metaweb(speciespool)) +
    AdjacentBernoulliDispersal(consumers, DispersalKernel(radius=1), 0.1) +
    RandomExtinction(consumers, probability=0.1) +
    LinearMortality(consumers, 0.01);

plantmodel =
    LogisticGrowth(resource) +
    AdjacentBernoulliDispersal(resource, DispersalKernel(radius=2), 0.1) ;

fullmodel = consumermodel + plantmodel;

dim = (50,50)
ntimesteps = 300

init = NamedTuple(merge(
    rand(Biomass, resource, Exponential(10), dim...),
    rand(Biomass, consumers, Exponential(10), dim...)));

arrayout = ArrayOutput(init, tspan=1:ntimesteps)
@time sim!(arrayout, fullmodel);

"""
    plots :
"""

mat = zeros(number_of_species, ntimesteps)

spnames = []
for t = 1:ntimesteps
  ind = 1

  for (k,v) in zip(keys(arrayout[t]), values(arrayout[t]))
    mat[ind, t] = sum(v)
    push!(spnames, k)
    ind += 1
  end
end

plt = plot(legend=:none)
xaxis!(plt, xticks=[i for i in 0:20:ntimesteps],"timestep")
yaxis!(plt, "log2(biomass)")
for s in 1:number_of_species
    speciesname = spnames[s]
    thiscol = get(ColorSchemes.thermal, trophicdict[String(speciesname)]/max(values(trophicdict)...))

    @show thiscol, speciesname, trophicdict[String(speciesname)]
    plot!(plt, 1:ntimesteps, mat[s,:], la=0.7, lw=3, c=thiscol, colorbar_title="trophic level")
end
plt

savefig(plt, "manyspeciestimeseries.png")

# ----
sz = dim[1]
ntimesteps = ntimesteps

biomass =  [zeros(sz,sz) for i in 1:ntimesteps]
speciesrichness = [zeros(sz,sz) for i in 1:ntimesteps]
for t in 1:ntimesteps
    now = arrayout[t]

    for (sp, abundmatrix) in zip(keys(now), values(now))
        for x in 1:sz, y in 1:sz
            if abundmatrix[x,y] > 0.0001
                speciesrichness[t][x,y] += 1
                biomass[t][x,y] += abundmatrix[x,y]
            end
        end
    end
end


text:
using POMDPs, POMDPGym, Test, Crux, Flux, Random, BSON
function test_solver(𝒮fn, mdp, π...)
    # run it once
    Random.seed!(0)
    S1 = 𝒮fn(deepcopy.(π)...)
    π1 = solve(S1, deepcopy(mdp))

    # run it again on the cpu with the same rng
    Random.seed!(0)
    S2 = 𝒮fn(deepcopy.(π)...)
    π2 = solve(S2, deepcopy(mdp))

    # Run it on the gpu
    Random.seed!(0)
    S3 = 𝒮fn(gpu.(deepcopy.(π))...)
    π3 = solve(S3, deepcopy(mdp))

    # compare the results
    s = rand(Crux.dim(state_space(mdp))...)
    try
        value(π[1], s)
        @test all(value(π1, s) .≈ value(π2, s))
        @test all(abs.(value(π2, s) .- value(π3, s)) .< 1e-2)
    catch
        if action(π1, s)[1] isa Symbol
            @test all(action(π1, s) .== action(π2, s))
            @test all(action(π2, s) .== action(π3, s))
        else
            @test all(action(π1, s) .≈ action(π2, s))
            @test all(abs.(action(π2, s) .- action(π3, s)) .< 1e-2)
        end
    end
end

## Training params
N = 100
ΔN = 50

## discrete RL
discrete_mdp = GridWorldMDP()
S = state_space(discrete_mdp)
A() = DiscreteNetwork(Chain(Dense(2, 32, relu), Dense(32, 4)), actions(discrete_mdp))
V() = ContinuousNetwork(Chain(Dense(2, 32, relu), Dense(32, 1)))
AC() = ActorCritic(A(), V())

test_solver((π) -> REINFORCE(π=π, S=S, N=N, ΔN=ΔN), discrete_mdp, A())
test_solver((π) -> A2C(π=π, S=S, N=N, ΔN=ΔN), discrete_mdp, AC())
test_solver((π) -> PPO(π=π, S=S, N=N, ΔN=ΔN), discrete_mdp, AC())
test_solver((π) -> DQN(π=π, S=S, N=N), discrete_mdp, A())

## Continuous RL
continuous_mdp = PendulumPOMDP()
S = state_space(continuous_mdp)
QSA() = ContinuousNetwork(Chain(Dense(3, 32, tanh), Dense(32, 1)))
V() = ContinuousNetwork(Chain( Dense(2, 32, relu), Dense(32, 1)))
A() = ContinuousNetwork(Chain(Dense(2, 32, relu), Dense(32, 1, tanh)), 1)
G() = GaussianPolicy(A(), zeros(Float32, 1))

test_solver((π) -> REINFORCE(π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, G())
test_solver((π) -> A2C(π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(G(), V()))
test_solver((π) -> PPO(π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(G(), V()))
test_solver((π) -> DDPG(π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(A(), QSA()))
test_solver((π) -> TD3(π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(A(), DoubleNetwork(QSA(), QSA())))
test_solver((π) -> SAC(π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(G(), DoubleNetwork(QSA(), QSA())))

# Continuous IL
𝒟_demo = expert_trajectories = BSON.load("examples/il/expert_data/pendulum.bson")[:data]
D(output=1) = ContinuousNetwork(Chain(DenseSN(3, 12, relu), DenseSN(12, output)))

test_solver((π, D) -> OnPolicyGAIL(D=D, 𝒟_demo=𝒟_demo, π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(G(), V()), D())
test_solver((π, D) -> OffPolicyGAIL(D=D, 𝒟_demo=𝒟_demo, π=π, S=S, N=50, ΔN=ΔN), continuous_mdp, ActorCritic(G(), DoubleNetwork(QSA(), QSA())), D(2))
test_solver(π -> BC(π=π, 𝒟_demo=𝒟_demo, S=S, opt=(epochs=1,), log=(period=50,)), continuous_mdp, A())
# NOTE: gradient penalty on the gpu only plays nicely with tanh, not relus in the discriminator?
test_solver((π) -> AdVIL(𝒟_demo=𝒟_demo, π=π, S=S, a_opt=(epochs=1,), log=(period=50,)), continuous_mdp, ActorCritic(A(), QSA()))
test_solver((π) -> SQIL(𝒟_demo=𝒟_demo, π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(G(), DoubleNetwork(QSA(), QSA())))
test_solver((π) -> AdRIL(𝒟_demo=𝒟_demo, π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, ActorCritic(G(), DoubleNetwork(QSA(), QSA())))
test_solver((π) -> ASAF(𝒟_demo=𝒟_demo, π=π, S=S, N=N, ΔN=ΔN), continuous_mdp, G())

# Batch RL
test_solver(π -> BatchSAC(π=π, 𝒟_train=𝒟_demo, S=S, a_opt=(epochs=1,)), continuous_mdp, ActorCritic(G(), DoubleNetwork(QSA(), QSA())))
test_solver(π -> CQL(π=π, 𝒟_train=𝒟_demo, S=S, a_opt=(epochs=1,)), continuous_mdp, ActorCritic(G(), DoubleNetwork(QSA(), QSA())))


text:
## Julia 1.1.0

# Initialize
using Plots
plotlyjs() # Backend PlotlyJS

x = range(-4 * pi, stop = 4 * pi, length = 1000)

plot(x, sin.(x))
savefig("./png/sin.png")

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
"""A face of Box{3}"""
struct Rectangle
    # Origin
    O::Point{3}
    # Two perp vector defining the plane ...
    t1::Point{3}
    t2::Point{3}
    # ... with normal
    normal::Point{3}
    # Lengths of sides
    l1::Float64
    l2::Float64
    # A
    # |
    # O--B
    function Rectangle(O::Point{3}, t1::Point{3}, t2::Point{3}, n::Point{3})
        l1 = norm(t1)
        t1 /= l1

        l2 = norm(t2)
        t2 /= l2

        n0 = cross(t1, t2)
        @assert n == normalize(n0)
        new(O, t1, t2, n, l1, l2)
    end
end

"""Construct from enumerating"""
function Rectangle(b::Box{3}, facet::Int)
    const dX = b.Y - b.X

    facet == 1 && return Rectangle(b.X, Point(dX[1], 0, 0), Point(0, dX[2], 0), Point(0, 0, 1))
    facet == 2 && return Rectangle(b.X, Point(0, dX[2], 0), Point(0, 0, dX[3]), Point(1, 0, 0))
    facet == 3 && return Rectangle(b.X, Point(0, 0, dX[3]), Point(dX[1], 0, 0), Point(0, 1, 0))

    facet == 4 && return Rectangle(b.Y, Point(-dX[1], 0, 0), Point(0, -dX[2], 0), Point(0, 0, 1))
    facet == 5 && return Rectangle(b.Y, Point(0, -dX[2], 0), Point(0, 0, -dX[3]), Point(1, 0, 0))
    facet == 6 && return Rectangle(b.Y, Point(0, 0, -dX[3]), Point(-dX[1], 0, 0), Point(0, 1, 0))

end

# Eval the curve using [0, 1]^2 parametrization
(r::Rectangle)(s::Number, t::Number) = r.O + s*r.l1*r.t1 + t*r.l2*r.t2

show(io::IO, b::Rectangle) = print(io, "$(b.O)-$(b(1, 1))")

"""
    Quad tree like division

    Bottom part    Top part
    |---|---|      |---|---|
    | 2 | 4 |      | 6 | 8 |
    |---|---|      |---|---|
    | 1 | 3 |      | 5 | 7 |
    |---|---|      |---|---|

"""
function split(box::Box{3})
    X, Y = box.X, box.Y
    M = 0.5(X + Y)
    [Box(Point(X[1], X[2], X[2]), Point(M[1], M[2], M[2])),
     Box(Point(X[1], M[2], X[2]), Point(M[1], Y[2], M[2])),
     Box(Point(M[1], X[2], X[2]), Point(Y[1], M[2], M[2])),
     Box(Point(M[1], M[2], X[2]), Point(Y[1], Y[2], M[2])),
     Box(Point(X[1], X[2], M[2]), Point(M[1], M[2], Y[2])),
     Box(Point(X[1], M[2], M[2]), Point(M[1], Y[2], Y[2])),
     Box(Point(M[1], X[2], M[2]), Point(Y[1], M[2], Y[2])),
     Box(Point(M[1], M[2], M[2]), Point(Y[1], Y[2], Y[2]))]
end

"""6 faces of the prism"""
function surface(box::Box{3})
    map(i -> Rectangle(box, i), 1:6)
end

"""Does a face(rectangle) collide with segment"""
function collides(line::Segment{3}, box::Rectangle)
    D = line.O - box.O
    D1 = dot(D, box.t1)
    D2 = dot(D, box.t2)
    D3 = dot(D, line.t)

    c1 = dot(box.t1, line.t)
    c2 = dot(box.t2, line.t)
    C = c1^2 + c2^2 - 1

    # @show [1 0 -c1; 0 1 -c2; c1 c2 -1]\[D1, D2, D3]
    # We solve the 3x3 system
    α = (-c1*c1/C + 1)*D1 + (-c1*c2/C)*D2 + (c1/C)*D3
    # @show α
    !(-EPS < α < box.l1 + EPS) && return false

    β = (-c1*c2/C)*D1 + (-c2*c2/C + 1)*D2 + (c2/C)*D3
    # @show β
    !(-EPS < β < box.l2 + EPS) && return false

    γ = (-c1/C)*D1 + (-c2/C)*D2 + (1/C)*D3
    # @show γ
    !(-EPS < γ < line.length + EPS) && return false

    true
end

"""See if segment collides with a box"""
function collides(segment::Segment{3}, box::Box{3})
    !bbox_collides(segment.box, box) && return false

    (segment.O ∈ box && segment.O ∈ box) && return true

    # Check collisions with cases --[---]--
    for line in surface(box)
        collides(segment, line) && return true
    end

    false
end

max_stars_repo_name: MiroK/fractional-dim


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.core.test.behaviour;
#
# import typedb.common.test.server.TypeDBRunner;
# import typedb.common.test.server.TypeDBSingleton;
# import org.junit.AfterClass;
#
# public abstract class BehaviourTest {
#     // The following code is for running the TypeDB Core distribution imported as an artifact.
#     // If you wish to debug locally against an instance of TypeDB that is already running in
#     // the background, comment out all the code in this file that references 'runner'
#     // and update ConnectionSteps to connect to TypeDBClient.DEFAULT_URI.
#
#     @AfterClass
#     public static void afterAll() {
#         TypeDBRunner server = TypeDBSingleton.getTypeDBRunner();
#         assert server != null;
#         server.stop();
#     }
# }


text:
using OrdinaryDiffEq, Test, LinearAlgebra

function f(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
function foop(u,p,t)
  dx = p[1]*u[1] - p[2]*u[1]*u[2]
  dy = -p[3]*u[2] + p[4]*u[1]*u[2]
  [dx,dy]
end

p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)
proboop = ODEProblem(foop,[1.0;1.0],(0.0,10.0),p)

SPECIAL_INTERPS = [Tsit5(),DP5(),SSPRK22(),OwrenZen3(),OwrenZen4(),OwrenZen5(),
                   BS5(),Vern6(),Vern7(),Vern8(),Vern9(),DP8(),Rosenbrock23(),
                   Rodas4(),Rodas5()]

y1 = zeros(2); y2 = zeros(2)
for alg in SPECIAL_INTERPS
  @show alg
  sol = solve(prob,alg,dt=0.0033,abstol=1e-14,reltol=1e-14)
  soloop = solve(proboop,alg,adaptive=false,tstops=sol.t,abstol=1e-14,reltol=1e-14)
  @test maximum(norm(soloop(t) - sol(t)) for t in 0:0.001:10) < 1e-10
  @test maximum(norm(soloop(y1,t) - sol(y2,t)) for t in 0:0.001:10) < 1e-10
end


text:
using DashBootstrapComponents

popovers = html_div([
    dbc_button("Hidden Arrow", id = "hide-arrow-target", className = "me-1", n_clicks = 0),
    dbc_popover(
        "I am a popover without an arrow!",
        target = "hide-arrow-target",
        trigger = "legacy",
        hide_arrow = true,
    ),
    dbc_button("Offset Popover", id = "offset-target", n_clicks = 0),
    dbc_popover(
        "I am a popover that's been offset!",
        target = "offset-target",
        trigger = "legacy",
        hide_arrow = true,
        offset = "50,20",
    ),
])

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
function _git_cmd(str::AbstractString;
                  adjust_PATH::Bool = true,
                  adjust_LIBPATH::Bool = true)
    git_path, env_mapping = _env_mapping(; adjust_PATH = adjust_PATH,
                                           adjust_LIBPATH = adjust_LIBPATH)
    new_env = copy(ENV)
    for p in env_mapping
        new_env[p[1]] = p[2]
    end
    new_cmd = Cmd(`$(git_path) $(split(str))`; env = new_env)
    return new_cmd
end


text:
include("binarytrees.jl")
include("buildtree.jl")
include("buildnodetree.jl")
include("comparetrees.jl")


text:
using ProjectedGPE, ApproxFun, FastGaussQuadrature, PyPlot, Combinatorics
M=350
tol = 1e-14
#c=Vector{BigFloat}(randn(M))+im*Vector{BigFloat}(randn(M))
#c=randn(M)+im*randn(M)
c=zeros(M)
c[end]=1.
N1=sum(abs.(c).^2)

#evaluate modes using BigFloats
#can then use ApproxFun.jl for polynomials, applying correct normalization.
#Combinatorics.jl needed to evaluate factorial of Vector.

n=0:BigFloat(M)-1;
Hnorm = sqrt(BigFloat(π))*BigFloat(2).^n.*factorial.(n)

#Inserting norm.
f=Fun(Hermite(),c./sqrt.(Hnorm))
f2=Fun(GaussWeight(Hermite()),c./sqrt.(Hnorm))
x=linspace(-1.5sqrt(2M),1.5sqrt(2M),1000)
plot(x,abs.(f.(x)).^2.*exp.(-x.^2))

x,w=gausshermite(M)

N2=sum(w.*abs.(f.(x)).^2)
N3=sum(w.*abs.(f2.(x).*exp.(x.^2)).^2)

Float64(N2)

abs((N1-N2)/N1)

3*M*tol

#block to create T matrix
M=10
x,w=gausshermite(M)
T=zeros(x*ones(1,M))
n=0:BigFloat(M)-1;
Hnorm = sqrt.(BigFloat(π))*BigFloat(2).^n.*factorial.(n)

for j = 0:M-1
  c=zeros(M)
  c[j+1]=1.
  #Hn=Fun(Hermite(),c./sqrt(Hnorm));
  #T[:,j+1]=Hn.(x).*exp(-x.^2/2)
  Hn=Fun(GaussWeight(Hermite()),c)
  T[:,j+1]=Hn.(x).*exp.(x.^2/2)
end

x1,y1,T1=nfieldtrans("Hermite",M,2,2)

mind = M
maximum(abs.(T1[:,mind]-T[:,mind]))

T2=eigmat("Hermite",M,x,2)
maximum(abs.(T1[:,mind]-T2[:,mind]))

#-------------------------------
#Make and test Laguerre transform and quadrature
M=27
α = 0.
ω = 1.
x,w=gausslaguerre(M,α)
T=zeros(x*ones(1,M))
n=0:M-1;
Lnorm = exp.(lgamma.(n+α+1)-lgamma.(n+1))/ω

for j = 0:M-1
  c=zeros(M)
  c[j+1]=1.
  Lna=Fun(Laguerre(α),c./sqrt(Lnorm));
  T[:,j+1]=Lna.(sqrt(ω)*x).*exp.(-sqrt(ω)*x/2).*(sqrt(ω)*x).^(α/2)
end

c=zeros(M);
c[end]=BigFloat(1.)

#Gauss Laguerre quadrature can absorb the factor x^α/2 using associated Laguerre
#c=randn(M)+im*rand(M)
N1 = sum(abs.(c).^2)
L=Fun(Laguerre(α),c./sqrt(Lnorm))

L2=Fun(WeightedLaguerre(α),c./sqrt(Lnorm))
N2=sum(w.*abs(L.(sqrt(ω)*x)).^2)

N3=sum(w.*abs(L2.(sqrt(ω)*x).*exp.(sqrt(ω)*x)./(sqrt(ω)x).^(α/2)).^2)
x=linspace(0,20sqrt(M),1000)
plot(x,abs.(L.(x)).^2.*exp.(-x).*x.^α)

Lnorm

#------------------------------------
#test Laguerre
M=26
α = 0.
ω = 1.
#x,w=gausslaguerre(M,α)

#T = eigmat("Laguerre",M,x/ω,ω,α)

x,w,T=nfieldtrans(M,2;basis="Laguerre")
#norm test

c=zeros(M);
c[end]=1.
N1 = sum(abs(c).^2)

Ln = T*c

N2=sum(w.*abs(Ln).^2)
f1 = Fun(GaussWeight(Hermite()),[1])

#------------------------------------
#test Hermite
M=21

ω = 10
#x,w=gausslaguerre(M,α)

#T = eigmat("Laguerre",M,x/ω,ω,α)

x,w,T=nfieldtrans("Hermite",M,2,ω)
#norm test

c=zeros(M);
c[end]=1.
N1 = sum(abs(c).^2)

Hn = T*c

N2=sum(w.*abs(Hn).^2)


text:
import PyPlot
const plt = PyPlot
printKeys(x) = [print("$k  ") for k in keys(x)]
using Distributions
#=
plt.ioff() # turn off interactive graphics. turn on with plt.ion()
plt.ion() # turn on interactive graphics. turn on with plt.ion()
=#

function pyRange(n::Int)
  return collect(0:(n-1))
end

function saveimg(path::String)
  plt.savefig(path)
  plt.close()
end

function plotZ(Z::Matrix{T}; xticks=true, kw...) where {T <: Number}
  J, K = size(Z)
  cm = plt.cm_get_cmap(:Greys)
  img = plt.imshow(Z, aspect="auto", vmin=0, vmax=1, cmap=cm; kw...);
  [ plt.axhline(y=i+.5, color="grey", linewidth=.5) for i in pyRange(J)];
  [ plt.axvline(x=i+.5, color="grey", linewidth=.5) for i in pyRange(K)];
  plt.yticks(pyRange(J), 1:J);
  if xticks
    plt.xticks(pyRange(K), 1:K, rotation=:vertical);
  end
  ax = plt.gca()
  ax[:tick_params](length=0)
  return ax
end

function plotY(Y::Matrix{T}, colorbar=true; kw...) where {T  <: Number}
  cm = plt.cm_get_cmap(:bwr)
  cm[:set_under](color=:blue)
  cm[:set_over](color=:red)
  cm[:set_bad](color=:black)
  img = plt.imshow(Y, aspect="auto", cmap=cm; kw...) # vmin, vmax
  if colorbar
    plt.plt[:colorbar]();
  end
  ax = plt.gca()
  ax[:tick_params](length=0)
  return img
end

### IBP ###
J = 32
K = 8
Z = Int.(rand(J, K) .> .6);

plotZ(Matrix(Z'))
ax = plt.gca()
ax[:set_xlabel]("Features", fontsize=12)
ax[:set_ylabel]("Observations", fontsize=12)
plt.tight_layout()
saveimg("img/Z.pdf")

### Heatmap ###
Y = randn(30000, 32) * 2;
Y[:, 3] .= NaN;

plotY(Y, vmin=-3, vmax=3)
plt.tight_layout()
saveimg("img/heatmap.pdf")

# yZ plot
#plt.rc_context(Dict("axes.edgecolor" => "grey", "xtick.color" => "grey", "ytick.color" => "grey"))

# Y
plt.subplot2grid((10, 1), (0, 0), rowspan=7)
plotY(Matrix(Y), false, vmin=-3, vmax=3);
plt.yticks(fontsize=5)
plt.plt[:colorbar]()
plt.xticks(pyRange(J), 1:J, rotation=:vertical);

# Z
plt.subplot2grid((10, 1), (8, 0), rowspan=3)
ax = plotZ(Matrix(Z'))
plt.xticks([])
plt.yticks(fontsize=5)
# Add ticks to other axis
ax = plt.gca()
ax2 = ax[:twinx]()
plt.plt[:colorbar](aspect=5)
# Change default ticks! Comment this out to see the difference.
ax2[:set_yticks](pyRange(K))
plt.yticks((K-1) / K * pyRange(K) .+ .5, (1:K) / 10, fontsize=5)
ax2[:tick_params](length=0)

plt.savefig("img/yZ.pdf", bbox_inches="tight")
plt.close()

#plt.rcdefaults()

# kde
function kde(x::Vector{T}; from::T=minimum(x), to::T=maximum(x), numPoints::Int=1000,
             bw::Float64=0.0, kernel::Function=z->pdf(Normal(), z)) where {T <: Number}
  n = length(x)
  if bw == 0.0
    bw = std(x) * n ^ (-0.2) * 1.06
  end
  x_grid = collect(range(from, stop=to, length=numPoints))
  return Dict(:dx => [sum(kernel.((xi .- x) / bw)) / (n * bw) for xi in x_grid], :x => x_grid)
end

# Density
dist = Gamma(3, 200)
@time x = rand(dist, 10000);
@time d = kde(x, numPoints=100, bw=100.);
plt.plot(d[:x], pdf.(dist, d[:x]), label=:truth, color=:blue)
plt.plot(d[:x], d[:dx], color=:orange, label=:kde)
x_ci = quantile(x, [.025, .975])
plt.fill_between(d[:x][x_ci[1] .< d[:x] .< x_ci[2]], d[:dx][x_ci[1] .< d[:x] .< x_ci[2]],
                color=:orange)
plt.tight_layout()
plt.legend()
saveimg("img/den.pdf")

# Histogram vs KDE
x = sort([rand(Normal(10, 1), 3000); rand(Normal(-10, 1), 3000)])
true_den = (pdf.(Normal(10, 1), x) .+ pdf.(Normal(-10, 1), x)) / 2
@time d = kde(x, numPoints=1000, from=-20., to=20., bw=.1)
@time d_auto = kde(x, numPoints=1000, from=-20., to=20.)
plt.plot(d[:x], d[:dx], color=:orange, label="kde: bw fixed at .1")
plt.plot(d_auto[:x], d_auto[:dx], color=:red, label="kde: auto-bw")
plt.plt[:hist](x, density=true, label=:hist, bins=50)
plt.plot(x, true_den, label="true density", c=:black)
plt.tight_layout()
plt.legend()
#plt.show()
saveimg("img/kde_vs_den.pdf")

# Axis everywhere
# https://matplotlib.org/examples/api/two_scales.html
fig, ax1 = plt.subplots()
ax1[:plot](randn(100))
ax1[:set_xlabel]("x")
# Make the y-axis label, ticks and tick labels match the line color.
ax1[:set_ylabel](:leftLabel, color=:blue)
ax1[:tick_params]("y", colors=:blue)

ax2 = ax1[:twinx]()
ax2[:set_ylabel](:blabla)
ax2[:tick_params](:blable)
saveimg("img/axisGalore.pdf")

# TODO: Plot in plot
# https://stackoverflow.com/questions/17458580/embedding-small-plots-inside-subplots-in-matplotlib
# https://matplotlib.org/api/_as_gen/matplotlib.pyplot.subplots.html
function infig_position(ax, rect)
  # Get current figure
  fig = plt.gcf()

  box = ax[:get_position]()
  inax_position = ax[:transAxes][:transform](rect[1:2])
  transFigure = fig[:transFigure][:inverted]()
  infig_pos = transFigure[:transform](inax_position)
  width = box[:width] * rect[3]
  height = box[:height] * rect[4]
  xpos = infig_pos[1]
  ypos = infig_pos[2]

  return xpos, ypos, width, height
end

function plotPost(y, ax::PyPlot.PyObject; rect=[.7, .7, .3, .2], showTrace::Bool=true,
                  showAcc::Bool=true, digits_ci=2, numPoints=100, c=:royalblue,
                  alpha_level=.05, useHistogram=false, normalizeHist=true, histBins=0,
                  accFontsize=6, annfs=6, hide_yaxis=true)
  # Get current figure
  fig = plt.gcf()

  if showTrace
    # Get box position / dimensions
    xpos, ypos, width, height = infig_position(ax, rect)

    subax = fig[:add_axes]([xpos, ypos, width, height])
    subax[:axis](:off)
    subax[:plot](1:length(y), y, c=:grey, linewidth=.5)
    accRate = length(unique(y)) / length(y)
    if showAcc
      subax[:set_title]("acc: $(accRate*100)%", fontsize=accFontsize)
    end
  end

  # plot histogram / density
  if useHistogram
    # Not really supported
    if histBins == 0
      ax[:hist](y, density=normalizeHist)
    else
      ax[:hist](y, density=normalizeHist, bins=histBins)
    end
  else
    if length(unique(y)) == 1
      ax[:annotate]("mean: $(round(mean(y), digits=digits_ci))", color=:red,
                    xy=(.8, .6), fontsize=annfs, xycoords="axes fraction")
      ax[:annotate]("SD: $(round(std(y), digits=digits_ci))",
                    xy=(.8, .4), fontsize=annfs, xycoords="axes fraction")
    else
      # Confidence Interval
      a_lower = alpha_level / 2
      a_upper = 1 - a_lower
      y_ci = quantile(y, [a_lower, a_upper])
      d_ci = kde(y, numPoints=numPoints, from=y_ci[1], to=y_ci[2])
      d = kde(y, numPoints=numPoints)
      ax[:plot](d_ci[:x], d_ci[:dx], c=c)
      ax[:plot](d[:x], d[:dx], c=c, alpha=.5)
      ax[:fill_between](d[:x][y_ci[1] .< d[:x] .< y_ci[2]],
                        d[:dx][y_ci[1] .< d[:x] .< y_ci[2]], color=c)
      ax[:fill_between](d[:x][minimum(y) .< d[:x] .< maximum(y)],
                        d[:dx][minimum(y) .< d[:x] .< maximum(y)], color=c, alpha=.3)
      ax[:set_xticks](round.(y_ci, digits=digits_ci))
      y_mean = mean(y)
      ax[:vlines](y_mean, ymin=0, ymax=d[:dx][argmin(abs.(d[:x] .- y_mean))],
                  color=:red, linewidth=1)

      # Plot Annotations
      ax[:annotate]("mean: $(round(y_mean, digits=digits_ci))", color=:red,
                    xy=(.8, .6), fontsize=annfs, xycoords="axes fraction")
      ax[:annotate]("$(Int(round((1-alpha_level) * 100)))% CI", color=c,
                    xy=(.8, .5), fontsize=annfs, xycoords="axes fraction")
      ax[:annotate]("SD: $(round(std(y), digits=digits_ci))",
                    xy=(.8, .4), fontsize=annfs, xycoords="axes fraction")
    end
  end

  ax[:spines]["right"][:set_visible](false)
  ax[:spines]["top"][:set_visible](false)

  if hide_yaxis
    ax[:axes][:yaxis][:set_ticklabels]([])
    ax[:spines]["left"][:set_visible](false)
    ax[:tick_params](length=0, axis="y")
  end

  ax[:spines]["bottom"][:set_color]("grey")
end

function plotPost(y; kw...)
  fig, ax = plt.subplots()
  plotPost(y, ax; kw...)
  return fig, ax
end

#plotPost(Y[:,1], useHistogram=true, histBins=30)
#plotPost(Y[:,1])

function plotPosts(Y; use_tight_layout=true,
                   fig_size_inches=missing,
                   fig_subplots_adjust=missing,
                   kw...)
  J = size(Y, 2)
  fig, ax = plt.subplots(J,J)
  if use_tight_layout
    fig[:tight_layout]()
  end

  if !ismissing(fig_size_inches)
    fig[:set_size_inches](fig_size_inches[1], fig_size_inches[2])
  end

  if !ismissing(fig_subplots_adjust)
    fig[:subplots_adjust](wspace=fig_subplots_adjust[1],
                          hspace=fig_subplots_adjust[2])
  end

  plotPosts(Y, fig, ax; kw...)
end

function plotPosts(Y, fig, ax; details=false, digits_cor=2, scatter_alpha=.7, scatter_s=.5,
                   warn=true, kw...)
  J = size(Y, 2)
  counter = 0
  for c in 1:J
    for r in 1:J
      counter += 1
      if r == c
        plotPost(Y[:, r], ax[counter]; kw...)
      elseif r > c
        ax[counter][:axis](:off)
        #xpos, ypos, w, h = infig_position(ax[counter], [.5, .5, .5, .5])
        cor_ycr = round(cor(Y[:, c], Y[:, r]), digits=digits_cor)
        if !isnan(cor_ycr)
          ax[counter][:text](0, .4, "r = $cor_ycr", fontsize=abs(cor_ycr)*10+10)
        else
          if warn
            @warn "Correlation was NaN."
          end
          ax[counter][:text](0, .4, "r = ???", fontsize=10)
        end
      else
        #ax[counter][:plot](Y[:, c], Y[:, r], c=:grey, linewidth=.5)
        ax[counter][:scatter](Y[:, c], Y[:, r], c=:grey, alpha=scatter_alpha, s=scatter_s)
        if !details
          ax[counter][:axis](:off)
        end
      end
    end
  end
  return fig, ax
end

#Y = randn(300, 5)
#Y[:,1] .= Y[:, 5] .+ 10
s = [[-1., 0., -3., 2.] [0., .5, .8, .2] [-3., .8, .6, .2] [2, .2, .2, -1.]]
S = s * s'
Y = rand(MvNormal([1,2,3,4], S), 1000)'
#Y[:, 2] .= 2.0

fig, ax = plotPosts(Y, digits_ci=1, fig_subplots_adjust=(.1, .1), fig_size_inches=(7, 7), annfs=8);
saveimg("img/plotPosts.pdf")

#= Or for more control
J = size(Y, 2)
fig, ax = plt.subplots(J,J);
fig[:tight_layout]()
fig[:set_size_inches](8, 8)
fig[:subplots_adjust](wspace=.1, hspace=.1)
plotPosts(Y, fig, ax, digits_ci=1);
saveimg("img/plotPosts.pdf")
=#

# Heatmap with annotations
W = rand(3, 10)
function plotW(W; digits_W=2, pad_colorbar=.05, shrink_colorbar=1, show_colorbar=true, kw...)
  I, K = size(W)
  img = plt.imshow(W, aspect=:auto; kw...)
  for i in pyRange(I)
    for k in pyRange(K)
      img[:axes][:text](k, i, round(W[i+1, k+1], digits=digits_W),
                        ha=:center, va=:center, color=:white)
    end
  end
  plt.xticks(pyRange(K), 1:K)
  plt.yticks(pyRange(I), 1:I)
  if show_colorbar
    plt.plt[:colorbar](pad=pad_colorbar, shrink=shrink_colorbar);
  end
  return img
end

img = plotW(W; pad_colorbar=.01, vmin=0., vmax=1.)
ax = plt.gca()
ax[:set_xlabel]("Features", fontsize=12)
ax[:set_ylabel]("Samples", fontsize=12)
plt.tight_layout()
saveimg("img/W.pdf")


text:
using ConservationLawsParticles
using RecursiveArrayTools, DifferentialEquations, Plots

# model
V(t, x) = 2 + 0.5sin(2x) + 0.3sin(2sqrt(2)*x) + 0.2cos(2sqrt(7)*x)
W(t, x) = x^2/2 - abs(x)
W′(t, x) = x - sign(x)
# W(t, x) = (a=abs(x); a*log(a+1/4))
# W′(t, x) = (a=abs(x); sign(x) * (log(a+1/4) + a/(a+1/4)))
mob(ρ) = max((ρ-1)^2 * (2ρ+1), 0)
smodel = SampledModel((V,), ((W′,),), (mob,))
imodel = IntegratedModel((V,), ((W,),), (mob,))

# initial condition
x0 = ArrayPartition(gaussian_particles(2, 401))

# time span
tspan = (0., 15.)

# ODE system for the particles
sprob = ODEProblem(velocities_gen!, x0, tspan, smodel)
iprob = ODEProblem(velocities_gen!, x0, tspan, imodel)

abstol = reltol = 2e-8

@time ssol = solve(sprob, BS5(); abstol=abstol, reltol=reltol);
@time isol = solve(iprob, BS5(); abstol=abstol, reltol=reltol);

# plot the particle trajectories
plot(legend=false, title="Trajectories", xlabel="time", ylabel="position")
plot!(ssol, vars=1:20:401; color=:blue)
plot!(isol, vars=1:20:401; color=:red)
savefig("traffic2.png")

# plot an animation of the density
anim = @animate for t in range(tspan...; step=1/24)
    p = plot(legend=false, xlims=(-2,20), ylims=(0,1),
        title="Density", xlabel="position", ylabel="density")
    x = ssol(t).x[1]
    min_dens = minimum(pwc_densities(x)[1][1,:,:], dims=1)[:];
    plot!(p, repeat(x, inner=(1,2))[1:20:401,:]',
        hcat(zero(min_dens), min_dens)[1:20:401,:]',
        color=:gray, opacity=.25)
    plot_density!(p, ssol(t); color=:blue)
    plot_density!(p, isol(t); color=:red)
end
gif(anim, "traffic2.gif")


text:
module TestUtils

using Test
using TMLEEpistasis
using MLJBase
using TOML
using BGEN
using DataFrames
using Serialization
using TMLE
using MLJLinearModels

include("helper_fns.jl")

@testset "Test parse_queries" begin
    build_query_file()
    queries = TMLEEpistasis.parse_queries(queryfile)
    expected_queries = [
        Query(case=(RSID_10="AG", RSID_100="AG"), control=(RSID_10="GG", RSID_100="GG"), name="QUERY_1"),
        Query(case=(RSID_10="AG", RSID_100="AA"), control=(RSID_10="GG", RSID_100="GG"), name="QUERY_2"),
        Query(case=(RSID_10="AG", RSID_100="AA"), control=(RSID_10="TT", RSID_100="GG"), name="QUERY_3"),
    ]
    test_queries(queries, expected_queries)

    rm(queryfile)
end

@testset "Test phenotypes parsing" begin
    # Fallback when no list is specified
    @test nothing === TMLEEpistasis.phenotypesnames(nothing)
    # Test with a list of phenotypes
    @test ["SAMPLE_ID", "CONTINUOUS_1"] == TMLEEpistasis.phenotypesnames(phenotypelist_file)
    # phenotypes loading
    phenotypes = TMLEEpistasis.load_phenotypes(binary_phenotypefile, nothing)
    @test size(phenotypes) == (490, 3)
    @test names(phenotypes) == ["SAMPLE_ID", "BINARY_1", "BINARY_2"]

    phenotypes = TMLEEpistasis.load_phenotypes(continuous_phenotypefile, phenotypelist_file)
    @test size(phenotypes) == (490, 2)
    @test names(phenotypes) == ["SAMPLE_ID", "CONTINUOUS_1"]
end

@testset "Test AdaptiveCV" begin
    # Continuous target
    cv = TMLEEpistasis.AdaptiveCV(CV())
    n_sample_to_nfolds = ((10, 10), (200, 20), (1000, 10), (5000, 5), (20000, 3))
    for (n, expected_nfolds) in n_sample_to_nfolds
        y = rand(n)
        ttp = MLJBase.train_test_pairs(cv, 1:n, y)
        @test length(ttp) == expected_nfolds
        @test ttp == MLJBase.train_test_pairs(CV(nfolds=expected_nfolds), 1:n, y)
    end
    # Categorical target
    cv = TMLEEpistasis.AdaptiveCV(StratifiedCV(nfolds=2))
    y = categorical(["a", "a", "a", "b", "b", "c", "c"])
    @test TMLEEpistasis.countuniques(y) == [3, 2, 2]
    ## neff < 30 => nfolds = 5*neff = 7
    ttp = MLJBase.train_test_pairs(cv, 1:7, y)
    @test length(ttp) == 7
    @test ttp == MLJBase.train_test_pairs(StratifiedCV(nfolds=7), 1:7, y)

    ## neff = 2500 => 10
    n = 50_500
    y = categorical(vcat(repeat([true], 50_000), repeat([false], 500)))
    @test TMLEEpistasis.countuniques(y) == [50_000, 500]
    ttp = MLJBase.train_test_pairs(cv, 1:n, y)
    @test length(ttp)== 10
    @test ttp == MLJBase.train_test_pairs(StratifiedCV(nfolds=10), 1:n, y)

end

end;

true


text:
module DataIndexer

export dataindex

STOPS_START = [".", "_", "-", "id_", "label ", "sdk_", "gcc-", "git-", "glew-", "bullet3-", "jdk1", "roi-manually-",
    "openssl-"]
STOPS_FULL = ["__pycache__", "node_modules", "site-packages", "bower_components", "src", "include", "runs", "man3",
    "lib", "bin", "doc", "docs", "build", "terminfo", "debug", "tool", "tools", "gcc", "op", "ops", "env", "envs",
    "3rdparty", "crops", "experiment", "experiments", "install", "json", "mmcv", "mmdet", "openssl", "dep", "blender",
    "pkg", "pkgs", "deps", "gentoo", "awscli", "intel", "software", "softwares", "gmake", "dependencies", "playground",
    "scripts", "curl", "tmp", "dmlc"]
STOPS_IN = ["linux", "finetune", "ffmpeg", "baseline", "inception", "opencv", "snapshots", "ccache", "epoch", "intel64",
    "caffe", "conda", "tensorrt", "cuda", "cudnn", "x86_64", "cmake", "tfevents", "egg-info", "backup", "_recovery_",
    "torch", "tensorflow", "python", "matlab", "tfrecords", "openmpi", "julia/stdlib/v"]
STOPS_END = ["frame", "frames", "model", "models", "config", "configs", "module", "modules", "checkpoint", "checkpoints",
    "logs", "result", "results", "workdir", "workdirs", "work_dir", "work_dirs", "snapshot", "snapshots", "cache",
    "cached", "_bk"]
SKIP_START = [".", "train-", "validation-", "json_", "model_", "pre_", "extract_feature_", "log_extract_",
    "eval_record_", "train_record_"]
SKIP_FULL = ["label", "label~", "_success", "checkpoint", "copying", "readme", "license", "dockerfile", "makefile",
    "copyright", "authors", "contributors", "thumbnail", "keyshot", "labels", "id_rsa", "log", "notice", "install",
    "sequence", "dummy", "changes", "md5file", "changelog", "n_frames"]
SKIP_IN = ["tfevents", "model.ckpt", "model.pth", ".ckpt.", ".txt.", "~backup", ".txt_", ".sh~", "model.data", ".json_"]
PASS_EXT = [".jpg", ".png", ".avi", ".jpeg", ".mkv", ".py", ".webm", ".pth.tar", ".pth", ".tif", ".tiff", ".mp4", ".mp3",
    ".bmp", ".npy", ".gif", ".txt", ".mdb", ".flv", ".json", ".json5", ".yuv", ".h264", ".list", ".sh", ".md", ".whl",
    ".egg", ".conf", ".yaml", ".rviz", ".launch", ".xml", ".ipynb", ".c", ".cpp", ".h", ".cu", ".tsv", ".o", ".iso",
    ".pub", ".zsh", ".bash", ".yml", ".pickle", ".meta", ".pkl", ".ply", ".bin", ".eps", ".gp", ".skv", ".mat", ".cuh",
    ".key", ".csv", ".proto", ".pcap", ".timestamps", ".m", ".nrrd", ".bak", ".pyc", ".hpp", ".pk", ".mov", ".svs",
    ".prototxt", "xlsx", ".lua", ".a", ".odgt", ".log", ".pyx", ".lock", ".rb", ".html", ".ini", ".caffemodel", ".dng",
    ".out", ".cmake", ".java", ".pdm", ".index", ".pdf", ".t7", ".pcd", ".dat", ".idx", ".prof", ".rst", ".js", ".in",
    ".cc", ".tfrecords", ".pbtxt", ".scsv", ".dict", ".proper", ".txt2", ".tag", ".value", ".mexa64", ".tpl", ".pptx",
    ".exe", ".pxd", ".binaryproto", ".calibration", ".ig", ".tx", ".tfrecord", ".filter", ".solverstate", ".mhd",
    ".zraw", ".raw", ".pb", ".pgm", ".ppm", ".uff", ".pyt", ".wts", ".cloc", ".runlog", ".pfm", ".model", ".lst", ".wav",
    ".2runlog", ".gt", "nii", ".sdf", ".ann", ".pm", ".mo", ".sample", ".la", ".so", ".so.1", ".so.2", ".so.3", ".so.4",
    ".spec", ".nii.gz", ".cfg", ".ckpt", ".arw", ".go", ".m4", ".plist", ".sln", ".vcproj", ".params", ".states",
    ".rsrc", ".sip", ".qm", ".onnx", ".kep", ".score", ".input", ".supp", ".natvis", ".lds", ".npz", ".exr", ".svg",
    ".rawplain16", ".nv21", ".p", ".mexw32", ".mexw64", ".mexglx", ".mexmaci64", ".mexmaci", ".mexmaci64", ".bag",
    ".ttf", ".acv", ".crdownload", ".pl", ".args", ".config", ".aif", ".ape", ".mpc", ".wmv", ".aa", ".bpg", ".btf",
    ".cr3", ".crw", ".dr4", ".vrd", ".eip", ".dcm", ".dpx", ".dv", ".djvu", ".dylib", ".elf", ".macho", ".flac", ".ogg",
    ".flif", ".fpf", ".swf", ".ppt", ".afm", ".dfont", ".pfa", ".pfb", ".raf", ".xcf", ".gpx", ".igc", ".kml", ".txt~",
    ".icc", ".itc", ".indd", ".j2c", ".jp2", ".lnk", ".lfp", ".mts", ".mie", ".miff", ".moi", ".mxf", ".mrw", ".nef",
    ".docx", ".dss", ".ods", ".opus", ".pgf", ".pict", ".xmp", ".psp", ".mobi", ".rw2", ".iiq", ".psd", ".m4a", ".webp",
    ".rtf", ".hdr", ".ra", ".ram", ".rm", ".r3d", ".x3f", ".pmp", ".torrent", ".ics", ".vcf", ".wtv", ".inx", ".numbers",
    ".cr2", ".cntk", ".feature", ".runner", ".toml", ".mk", ".stamp", ".hh", ".f", ".sub", ".am", ".ac", ".vim", ".el",
    ".vbs", ".xls", ".lml", ".db", ".avro", ".flo", ".dump", ".lic", ".filtered", ".path", ".rockspec", ".dav", ".hik",
    ".nef", ".template", ".bat", ".interp", ".vtp", ".xsl", ".ebuild", ".patch", ".mask", ".defaults", ".inc", ".sens",
    ".py~", ".md5", ".info", ".trans", ".vg", ".cdf", ".h5", ".ta.gz", ".txz", ".old", ".np", ".inv", ".pem", ".obj",
    ".heic", ".nv12", ".eff", ".dll", ".res", ".lck", ".tmp", ".cache", ".bullet", ".plan", ".filters", ".xcscheme",
    ".diff", ".cav", ".install", ".properties", ".f90", ".dist", ".rc", ".stl", ".xacro", ".urdf", ".dae", ".mtl",
    ".inl", ".ino", ".glsl", ".cl", ".mm", ".command", "bsp", ".d", ".make", ".pov", ".osl", ".manifest", ".nsi", ".ico",
    ".po", ".blend", ".desktop", ".spi1d", ".spi3d", ".spimtx", ".icns", ".blender", ".pt", ".pak", ".css", ".dtd",
    ".gir", ".its", ".loc", ".url", ".sed", ".sin", ".header", ".charset", ".tcl", ".alias", ".qml", ".qmltypes", ".prf",
    ".qph", ".jar", ".gradle", ".mpg", ".rmvb", ".part", ".okl", ".rec", ".graph", ".run", ".bk", ".autopkg", ".repo",
    ".nvvp", ".vaml", ".record", ".txt1", ".txt3", ".asv", ".j", ".jp", ".x-ms-bmp", ".py.1", ".py.2", ".faissindex",
    ".ui", ".qrc", ".doxy", ".vtk", ".pclzf", ".thrift", ".so.3.0", ".port", ".json1", ".jsonl", ".label", ".so.3.0.2",
    ".pot", ".gmo", ".netrc", ".aria2", ".success", ".ids", ".aar", ".pc", ".mmp", ".pkg", ".iby", ".inf", ".mak",
    ".def", ".com", ".opt", ".msg", ".sdl", ".cnf", ".crt", ".csr", ".cacert", ".der", ".dhp", ".prm", ".tex", ".thm",
    ".report", ".mf", ".landmark", ".mapfile", ".ymlst", ".dummy", ".env", ".guess", ".lib", ".xyz", ".rank0", ".rank1",
    ".rank2", ".rank3", ".so.4.4.5", ".rle", ".py.1", ".py.2", ".output", ".lrprev", ".lrcat", ".orig", ".map", ".woff",
    ".woff2", ".eot", ".timelog", ".vmd", ".htm", ".sfx", ".pts", ".ts", ".doc", ".jfif", ".part1", ".part2", ".part3",
    ".part4", ".part5", ".part6", ".part7", ".part8", ".feat", ".gemspec", ".rake", ".encrypt", ".m~", ".rpm", ".nv",
    ".nv2", ".fmt", ".t", ".textgrid", ".dot", ".matlab", ".sd", ".chk", ".xib", ".pddm", ".swift", ".php", ".vcxproj",
    ".props", ".targets", ".asm", ".hex", ".errwarn", ".class", ".rda", ".r", ".rmd", ".wt", ".app", ".xib", ".cnf",
    ".c1", ".c2", ".c3", ".ncf", ".stdout", ".brew", ".port", ".status"]

FILE_EXT = [".zip", ".tar", ".tar.gz", ".rar", ".tgz", ".tar.bz2", ".7z", ".tar.xz", ".zip.1", ".rar.1", ".zip.001",
    ".z01"]
NUM_SHOW_DIR = 222
NUM_SUPPRESS = 11111
PATTERN_DETECT_THRESHOLD = 33
MIN_DEPTH = 1
DEBUG = [false]

iskeepdir(d) = try
    return !(match(r"^\d+$", d) !== nothing ||
            any(x->endswith(lowercase(d), x), STOPS_END) ||
            any(x->startswith(lowercase(d), x), STOPS_START) ||
            any(x->lowercase(d) == x, STOPS_FULL) ||
            any(x->occursin(x, lowercase(d)), STOPS_IN))
catch
    return false
end

isfileheurstic(f) = try
    return any(x->endswith(lowercase(f), x), FILE_EXT)
catch
    return false
end

iskeepfile(f) = try
    if any(x->endswith(lowercase(f), x), FILE_EXT)
        ('0' <= f[1] <= '9') || return true
        return tryparse(Int, rsplit(f, "."; limit = 2)[1]) === nothing
    end
    return false
catch
    return false
end

ispasscheck(f, p) = try
    any(x->startswith(lowercase(f), x), SKIP_START) && return true
    any(x->occursin(x, lowercase(f)), SKIP_IN) && return true
    any(x->endswith(lowercase(f), x), PASS_EXT) && return true
    any(x->lowercase(f) == x, SKIP_FULL) && return true
    if startswith(f, p)
        s = split(f, r"[._-]", limit = 2, keepempty = false)
        length(s) == 2 && s[1] == p && tryparse(Int, s[2]) !== nothing && return true
    end
    return false
catch
    return true
end

haspatternquick(list) = try
    l = length(list[1])
    if l > 0 && !occursin("batch", list[1]) && all(length(f) == l for f in list)
        l > 1 && length(list) > PATTERN_DETECT_THRESHOLD * 3 && return true
        l > 2 && length(list) > PATTERN_DETECT_THRESHOLD && return true
        l > 6 && length(list) > PATTERN_DETECT_THRESHOLD / 2 && return true
        l > 10 && length(list) > PATTERN_DETECT_THRESHOLD / 3 && return true
        l > 15 && length(list) > PATTERN_DETECT_THRESHOLD / 4 && return true
        l > 30 && length(list) > PATTERN_DETECT_THRESHOLD / 5 && return true
        l > 40 && length(list) > PATTERN_DETECT_THRESHOLD / 6 && return true
        l > 50 && length(list) > PATTERN_DETECT_THRESHOLD / 9 && return true
    end
    return false
catch
    return false
end

haspattern(list, p) = try
    list = [rsplit(f, "."; limit = 2)[1] for f in list]
    if length(list) > PATTERN_DETECT_THRESHOLD / 6 && !occursin("batch", list[1]) && (m = match(r"^(.*?[^\d])\d+$", list[1])) !== nothing
        length(m.captures[1]) > 1 && all(startswith.(list, m.captures[1])) && return true
    end
    length(list) < PATTERN_DETECT_THRESHOLD / 4 && return false
    all(occursin.("batch", list)) && return false
    nu = length(unique(length.(list)))
    nu == 1 && length(list[1]) == 2 && length(list) > PATTERN_DETECT_THRESHOLD * 2 && return true
    nu == 1 && length(list[1]) > 2 && length(list) > PATTERN_DETECT_THRESHOLD && return true
    nu == 1 && length(list[1]) > 15 && length(list) > PATTERN_DETECT_THRESHOLD / 3 && return true
    nu == 2 && length(list) > PATTERN_DETECT_THRESHOLD * 3 && return true
    nu == 1 && !occursin("batch", list[1]) && all(startswith.(list, p)) && return true
    if length(list) > PATTERN_DETECT_THRESHOLD / 2 && match(r"^(.*?)[\d_-]*$", list[1]).captures[1] == ""
        all(match(r"^(.*?)[\d_-]*$", f).captures[1] == "" for f in list) && return true
    end
    if length(list) > PATTERN_DETECT_THRESHOLD * 3 && match(r"^(.*?[^\d])\d+$", list[1]) !== nothing
        nreduce = length(unique([match(r"^(.*?)[\d_-]*$", f).captures[1] for f in list]))
        nreduce / length(list) < 0.1 && return true
    end
    if length(list) > PATTERN_DETECT_THRESHOLD * 3 && nu < 10
        lens = length.(list)
        unilens = unique(lens)
        any(sum(lens .== l) > length(list) - 5 for l in unilens) && return true
        lens = collect(setdiff(Set(unilens), Set([24, 32])))
        length(lens) <= 5 && sum(lens .== lens[1]) < max(10, length(list) / 50) && return true
    end
    if length(list) > PATTERN_DETECT_THRESHOLD * 4
        pattern_starts = ["giphy-", "bing-", "yt-", "vine-", "getty-", "flickr-"]
        any(sum(startswith.(list, site)) > PATTERN_DETECT_THRESHOLD for site in pattern_starts) && return true
    end
    return false
catch
    return false
end

function walk(root, depth = 0, parentname = "^")
    content = nothing
    try
        content = readdir(root)
    catch
        c = Channel(0)
        close(c)
        return c
    end

    dirs = Vector{eltype(content)}()
    files = Vector{eltype(content)}()

    if length(content) < NUM_SUPPRESS * 2
        if !haspatternquick(content)
            for name in content
                ispasscheck(name, parentname) && continue
                try
                    if depth == 0 || length(dirs) < NUM_SHOW_DIR
                        DEBUG[] && !isfileheurstic(name) && !isdir(joinpath(root, name)) && println(joinpath(root, name))
                        if (!isfileheurstic(name)) && isdir(joinpath(root, name))
                            iskeepdir(name) && push!(dirs, name)
                        else
                            iskeepfile(name) && push!(files, name)
                        end
                    else
                        push!(files, name)
                    end
                catch
                end
            end
        end
    end

    # if > NUM_SHOW_DIR, do not go deeper, transform them to files
    if depth > MIN_DEPTH && length(dirs) >= NUM_SHOW_DIR
        prepend!(files, dirs)
        empty!(dirs)
    end

    # if > NUM_SUPPRESS, do not show files
    if length(files) > NUM_SUPPRESS
        empty!(files)
    end

    # detect pattern
    haspattern(dirs, parentname) && empty!(dirs)
    haspattern(files, parentname) && empty!(files)

    function _it(chnl)
        put!(chnl, (root, dirs, files))
        for dir in dirs
            path = joinpath(root, dir)
            if !islink(path)
                for (root_l, dirs_l, files_l) in walk(path, depth + 1, dir)
                    put!(chnl, (root_l, dirs_l, files_l))
                end
            end
        end
    end

    return Channel(_it)
end

function dataindex(ROOTDIR)
    for (root, dirs, files) in walk(ROOTDIR)
        DEBUG[] || println(root)
        sleep(0.001)
        for f in files
            DEBUG[] || println(joinpath(root, f))
        end
    end
end

end # module


text:
######################
# Logit and Logistic #
######################
using StatsFuns: logit, logistic

struct Logit{N, T<:Real} <: Bijector{N}
    a::T
    b::T
end
function Logit(a, b)
    T = promote_type(typeof(a), typeof(b))
    Logit{0, T}(a, b)
end
up1(b::Logit{N, T}) where {N, T} = Logit{N + 1, T}(b.a, b.b)
# For equality of Logit with Float64 fields to one with Duals
Base.:(==)(b1::Logit, b2::Logit) = b1.a == b2.a && b1.b == b2.b

(b::Logit{0})(x::Real) = _logit(x, b.a, b.b)
(b::Logit{0})(x) = _logit.(x, b.a, b.b)
(b::Logit{1})(x::AbstractVector) = _logit.(x, b.a, b.b)
(b::Logit{1})(x::AbstractMatrix) = _logit.(x, b.a, b.b)
(b::Logit{2})(x::AbstractMatrix) = _logit.(x, b.a, b.b)
(b::Logit{2})(x::AbstractArray{<:AbstractMatrix}) = map(b, x)
_logit(x, a, b) = logit((x - a) / (b - a))

(ib::Inverse{<:Logit{0}})(y::Real) = _ilogit(y, ib.orig.a, ib.orig.b)
(ib::Inverse{<:Logit{0}})(y) = _ilogit.(y, ib.orig.a, ib.orig.b)
(ib::Inverse{<:Logit{1}})(x::AbstractVecOrMat) = _ilogit.(x, ib.orig.a, ib.orig.b)
(ib::Inverse{<:Logit{2}})(x::AbstractMatrix) = _ilogit.(x, ib.orig.a, ib.orig.b)
(ib::Inverse{<:Logit{2}})(x::AbstractArray{<:AbstractMatrix}) = map(ib, x)
_ilogit(y, a, b) = (b - a) * logistic(y) + a

logabsdetjac(b::Logit{0}, x::Real) = logit_logabsdetjac(x, b.a, b.b)
logabsdetjac(b::Logit{0}, x) = logit_logabsdetjac.(x, b.a, b.b)
logabsdetjac(b::Logit{1}, x::AbstractVector) = sum(logit_logabsdetjac.(x, b.a, b.b))
logabsdetjac(b::Logit{1}, x::AbstractMatrix) = vec(sum(logit_logabsdetjac.(x, b.a, b.b), dims = 1))
logabsdetjac(b::Logit{2}, x::AbstractMatrix) = sum(logit_logabsdetjac.(x, b.a, b.b))
logabsdetjac(b::Logit{2}, x::AbstractArray{<:AbstractMatrix}) = map(x) do x
    logabsdetjac(b, x)
end
logit_logabsdetjac(x, a, b) = -log((x - a) * (b - x) / (b - a))


text:
include(joinpath("deps", "Compute.jl"))

using GridDX
using GridShape

function main(argd)
    basis_matrix = eye(3)
    box_size = tuple([convert(Int, x) for x in argd["size"]]...)
    origin = convert(Array{Float64, 1}, argd["origin"])
    for output_file in argd["outputs"]
        fn_compute = Compute.get_fn(argd["shape"], argd["params"])
        @time grid = GridDX.make(
            fn_compute,
            box_size,
            basis_matrix;
            origin=origin
        )

        sizes = [float(x) for x in box_size]
        corner = -sizes/2.0 + origin

        @time GridDX.save(
            output_file,
            GridDX.clone(grid; origin=corner)
        )
    end
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
using CartesianFDM

n = (6, 6)
bc = ntuple(dirichlet, length(n))

ops = fdmoperators(bc, n)

### symbolic
ε = only(@variables(ε))

# metrics
V = replace!(mask(ops, scalar(:V, n))) do el
    iszero(el) ? ε : el
end
A = mask(ops, vector(:A, n))

# primary
U = vector(:U, n)
W = vector(:W, n)

# secondary
E = strainrate(ops, A, V, U, W)
L = divergence2(ops, A, E)

# remove ε
L = substitute.(L, Ref(Dict([ε => 0])))

###
using SparseArrays
using Base.Iterators

J = map(product(L, U)) do (Li, Uj)
    sparse(Symbolics.jacobian(Li, Uj))
end

### numerical
V̂ = map(mask(ops, ones(Float64, prod(n)))) do x
    iszero(x) ? eps(x) : x
end
Â = mask(ops, ntuple(length(n)) do i
                  ones(Bool, prod(n))
              end)
Ŵ = zero.(Â)

Ê = strainrate(ops, Â, V̂, U, Ŵ)
L̂ = divergence2(ops, Â, Ê)

Ĵ = map(product(L̂, U)) do (Li, Uj)
    sparse(Symbolics.jacobian(Li, Uj))
end

### vector heat conduction
@variables t p
(rhs, rhs!) = build_function(vcat(L̂...), vcat(U...), p, t)

using DifferentialEquations

tspan = (0.0, 100.0)
Û = ones(length(n) * prod(n))
prob = ODEProblem(eval(rhs!), Û, tspan)
sol = solve(prob)

using Plots

heatmap(first(reshape(sol(0.1), n..., 2)))


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
struct graph
    pb::problem
    mono_pb::mono_problem
    layer::Vector{vertex} #last layer of the graph
end
Base.isempty(g::graph) = isempty(g.layer)

function build_graph(Δ::Triangle, output::Bool)
    xr = Δ.xr
    xs = Δ.xs
    pb = xr.pb
    λ1, λ2 = Δ.λ
    mono_pb = mono_problem(pb, λ1, λ2)
    # @assert obj(solve_mono(mono_pb)) == Δ.ub "$xr => $xs isn't a supporting edge of Y ($(solution(pb,solve_mono(mono_pb))))"
    reduce!(mono_pb, Δ, output)
    size(mono_pb) == 0 && return graph(pb, mono_pb, vertex[])

    s = source(mono_pb)

    current_layer = [s]
    for l = 1:size(mono_pb)

        # @assert mono_pb.variables[l] == current_layer[1].i
        # @assert l == current_layer[1].layer

        previous_layer = current_layer
        current_layer = vertex[]

        #Adds all vertices in the next layer obtained without picking the object
        for v in previous_layer
            if relax(v, mono_pb) >= Δ.lb #&& relax_z1(v, sortperm_z1) >= obj_1(xr) || relax_z2(v, sortperm_z2) >= obj_2(xs)
                push!(current_layer, vertex_skip(v, mono_pb))
            end
        end

        #Adds all vertices in the next layer obtained by picking the object (if possible)
        for v in previous_layer

            let wplus1 = v.w + pb.w[v.i] #Calculate the weight of the solution when we pick item i
                wplus1 > mono_pb.c && continue #If the solution is too heavy for the knapsack, skip it.
                child = vertex_keep(v, mono_pb) #else, create the vertex

                # Check if a node with that weight already exists so we can merge them.
                # The nodes in the layer are sorted by their weight so we can use a dichotomic search
                range = searchsorted(current_layer, wplus1, lt = weight_lt) # (weight_lt compares a vertex to an Integer weight)
                if isempty(range)
                    # ∄ node with that weight, we can add it (at the right place) to the current layer
                    insert!(current_layer, range.start , child)
                else
                    # a node with this weight already exists, we merge them
                    merge!(current_layer[range.start], child)
                end
            end
        end
    end
    return graph(pb, mono_pb, current_layer)
end

#Calculates the relaxation for the combined problem
function relax(v::vertex, mpb::mono_problem)

    #we assume the items are already sorted by utility
    #we only consider picking items which come after v.i

    vars = view(variables(mpb), v.layer+1:size(mpb))::SubArray{Int}
    zλ= v.zλ
    w = v.w
    c = mpb.c
    i = 1

    while i <= length(vars) && w + mpb.w[vars[i]] <= c
        var = vars[i]
        w += mpb.w[var]
        zλ += mpb.p[var]
        i += 1
    end

    if i <= length(vars) && length(vars) > 1
        cleft = mpb.c - w
        if i == 1
            varplus1 = vars[i+1]
            zλ += floor(Int, cleft * mpb.p[varplus1]/mpb.w[varplus1])
        elseif i == length(vars)
            varminus1 = vars[i-1]
            zλ += floor(Int, mpb.p[var] - (mpb.w[var] - cleft)*(mpb.p[varminus1]/mpb.w[varminus1]))
        else
            var = vars[i]
            varplus1 = vars[i+1]
            varminus1 = vars[i-1]
            U0 = zλ + floor(Int, cleft * mpb.p[varplus1]/mpb.w[varplus1])
            U1 = zλ + floor(Int, mpb.p[var] - (mpb.w[var] - cleft)*(mpb.p[varminus1]/mpb.w[varminus1]))
            zλ = max(U0,U1)
        end
    end

    return zλ
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
module TZLookup

using TimeZoneLookup
using TimeZoneLookup: buildnode, faceof, E, V
using Test

@testset "faceof" begin
    @test faceof(EData(E(V(0, 0), V(0, 10)), -1, 1), V(5, 5)) == 1
    @test faceof(EData(E(V(0, 0), V(0, 10)), -1, 1), V(-5, 5)) == -1
    @test faceof(EData(E(V(0, 0), V(10, 0)), -1, 1), V(5, 5)) == -1
    @test faceof(EData(E(V(0, 0), V(10, 0)), -1, 1), V(5, -5)) == 1
    @test faceof(EData(E(V(0, 0), V(10, 10)), -1, 1), V(0, 10)) == -1
    @test faceof(EData(E(V(0, 0), V(10, 10)), -1, 1), V(10, 0)) == 1

    @test faceof(EData(E(V(100, 100), V(100, 110)), -1, 1), V(105, 1055)) == 1
    @test faceof(EData(E(V(100, 100), V(100, 110)), -1, 1), V(95, 105)) == -1
    @test faceof(EData(E(V(100, 100), V(110, 100)), -1, 1), V(105, 105)) == -1
    @test faceof(EData(E(V(100, 100), V(110, 100)), -1, 1), V(105, 95)) == 1
    @test faceof(EData(E(V(100, 100), V(110, 110)), -1, 1), V(105, 110)) == -1
    @test faceof(EData(E(V(100, 100), V(110, 110)), -1, 1), V(110, 100)) == 1
end

@testset "update!" begin
    S = [EData(E(V(10, 10), V(30, 10)), 1, -1),
         EData(E(V(10, 10), V(20, 20)), -1, 1),
         EData(E(V(30, 10), V(40, 20)), 1, -1),
         EData(E(V(20, 20), V(40, 20)), -1, 1)
        ]
    bb = buildnode(V(0, 0), V(50, 50), E(V(0, 50), V(50, 50)), E(V(0, 0), V(50, 0)), -1)
    root = bb.node

    update!(root, S[2])

    @test query(root, V(5, 15)).face == -1
    @test query(root, V(25, 15)).face == -1
    @test query(root, V(15, 5)).face == 1
    @test query(root, V(15, 25)).face == -1

    update!(root, S[3])

    @test query(root, V(25, 15)).face == -1
    @test query(root, V(45, 15)).face == -1
    @test query(root, V(35, 10)).face == -1
    @test query(root, V(35, 20)).face == 1
end

end # module


text:
struct Null
end

isnull(x) = issingletontype(typeof(x))

function merge_merge(fields::Vector{Vector{Symbol}}, options::PGFPlotsX.Options...;uniq::Bool=true)
    # fields is a vector of subvector of symbols
    # fields in a subvector, they should be deduped together
    # for field in fields, append them, overwrite otherwise
    mentioned_keys = vcat(fields...)
    @assert allunique(mentioned_keys) "elem appears multiple times in multiple subvectors"
    all_keys = unique(mapreduce(o->collect(keys(o.dict)), vcat ,options))
    left_keys = setdiff(all_keys,mentioned_keys)
    merged = reduce(options; init=PGFPlotsX.Options()) do x,y
        dict = OrderedDict()
        for fs in fields
            x_f = zip([vcat(get(x.dict, f, [])) for f in fs]...)
            y_f = zip([vcat(get(y.dict, f, [])) for f in fs]...)
            new_f = vcat(x_f...,y_f...)
            if length(new_f) == 0
                continue
            end
            if uniq
                unique!(new_f)
            end
            merge!(dict,OrderedDict(fk=>collect(fv) for (fk,fv) in zip(fs, zip(new_f...))))
        end
        merge!(dict,OrderedDict(fk=>get(merge(y.dict,x.dict), fk, nothing) for fk in left_keys))
        return PGFPlotsX.Options(dict, x.print_empty|y.print_empty)
    end
    return merged
end

function PGFPlotsX.print_tex(io::IO, n::Null)
end

struct LineWithEnds <: PGFPlotsX.TikzElement
    end1::Tuple
    end2::Tuple
    options::PGFPlotsX.Options
end

function PGFPlotsX.print_tex(io::IO, line::LineWithEnds)
    PGFPlotsX.print_indent(io, "%%%%%%%%% plot line of the coefplot %%%%%%%%%%")
    print(io, "\\draw")
    PGFPlotsX.print_options(io, line.options; newline = false)
    println(io, "(axis cs:$(join(line.end1,","))) -- (axis cs:$(join(line.end2,",")));")
end

struct Dot <: PGFPlotsX.TikzElement
    loc::Tuple
    options::PGFPlotsX.Options
end

default_dot_options() = merge(PGFPlotsX.Options(:circle => nothing, :"inner sep" => "0pt",:"minimum size" => "4pt"), color_as_fill_option(:navy))

function PGFPlotsX.print_tex(io::IO, dot::Dot)
    PGFPlotsX.print_indent(io, "%%%%%%%%% plot dot of the coefplot %%%%%%%%%%")
    print(io, "\\filldraw")
    print(io,"(axis cs:$(join(dot.loc,","))) node")
    PGFPlotsX.print_options(io, dot.options; newline = false)
    println(io, "{}; ")
end

struct AbstractCaption <: PGFPlotsX.TikzElement
    caption::String
    options::PGFPlotsX.Options
end

default_note_options() = PGFPlotsX.Options(:name => "note", :anchor => "north west", :font => "{\\fontsize{4.5}{4.5}\\selectfont}",
                                           :"text width" => "0.7\\textwidth", :at => "(current axis.outer south west)")

function PGFPlotsX.print_tex(io::IO, caption::AbstractCaption)
    PGFPlotsX.print_indent(io, "%%%%%%%%% plot caption %%%%%%%%%%")
    print(io, "\\draw node")
    PGFPlotsX.print_options(io, caption.options; newline = false)
    println(io, "{"*caption.caption*"};")
end

mutable struct Legend <: PGFPlotsX.TikzElement
    l::String
    line_options::PGFPlotsX.Options
    dot_options::PGFPlotsX.Options
end

function PGFPlotsX.print_tex(io::IO, lgd::Legend)
    PGFPlotsX.print_indent(io, "%%%%%%%%% plot legends %%%%%%%%%%")
    print(io, "\\addlegendimage{")
    PGFPlotsX.print_opt(io, lgd.line_options)
    PGFPlotsX.print_opt(io, ", mark = *, mark options = {")
    PGFPlotsX.print_opt(io, lgd.dot_options)
    print(io, "}};")
end

Legend(l::String,representitive_scoptions::SinglecoefplotOption) = Legend(l, get_line_options(representitive_scoptions),
                                                                          merge(get_dot_options(representitive_scoptions),
                                                                          color_as_draw_option(changeopacity(representitive_scoptions.dotcolor,0)))) # need to set draw opacity = 0

mutable struct Plotable <: PGFPlotsX.TikzElement
    options::PGFPlotsX.Options
    note::Union{AbstractCaption,Null}
    contents::Vector{T} where T
    function Plotable(options::PGFPlotsX.Options, note::Union{AbstractCaption,Null}, contents...)
        new(options, note, collect(contents))
    end
end

function PGFPlotsX.print_tex(io::IO, p::Plotable)
    print(io, "\\begin{axis}")
    PGFPlotsX.print_options(io, p.options)
    PGFPlotsX.print_indent(io) do io
        for elt in reverse(p.contents)
            PGFPlotsX.print_tex(io, elt)
        end
    end
    println(io, "\\end{axis}")
    PGFPlotsX.print_tex(io, p.note)
end


text:
function shownodes(node::Node, k::Int = 0)
    if has_children(node)
        for nd in node.children
            print("k = $(k + 1): ")
            println(nd.item_ids)
        end
        for nd in node.children
            shownodes(nd, k+1)
        end
    end
end

function randstr(n::Int, len::Int = 16)
    vals = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m",
            "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"]
    upper = map(uppercase, vals)
    append!(vals, upper)
    append!(vals, map(string, 0:9))
    res = Array{String,1}(undef, n)
    for i = 1:n
        res[i] = join(rand(vals, len))
    end
    res
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: forecastquery
using AWS.Compat
using AWS.UUIDs

"""
    QueryForecast()

Retrieves a forecast for a single item, filtered by the supplied criteria. The criteria is a key-value pair. The key is either item_id (or the equivalent non-timestamp, non-target field) from the TARGET_TIME_SERIES dataset, or one of the forecast dimensions specified as part of the FeaturizationConfig object. By default, QueryForecast returns the complete date range for the filtered forecast. You can request a specific date range. To get the full forecast, use the CreateForecastExportJob operation.  The forecasts generated by Amazon Forecast are in the same timezone as the dataset that was used to create the predictor.

# Required Parameters
- `Filters`: The filtering criteria to apply when retrieving the forecast. For example, to get the forecast for client_21 in the electricity usage dataset, specify the following:  {\"item_id\" : \"client_21\"}  To get the full forecast, use the CreateForecastExportJob operation.
- `ForecastArn`: The Amazon Resource Name (ARN) of the forecast to query.

# Optional Parameters
- `EndDate`: The end date for the forecast. Specify the date using this format: yyyy-MM-dd'T'HH:mm:ss (ISO 8601 format). For example, 2015-01-01T20:00:00.
- `NextToken`: If the result of the previous request was truncated, the response includes a NextToken. To retrieve the next set of results, use the token in the next request. Tokens expire after 24 hours.
- `StartDate`: The start date for the forecast. Specify the date using this format: yyyy-MM-dd'T'HH:mm:ss (ISO 8601 format). For example, 2015-01-01T08:00:00.
"""
query_forecast(Filters, ForecastArn; aws_config::AWSConfig=global_aws_config()) = forecastquery("QueryForecast", Dict{String, Any}("Filters"=>Filters, "ForecastArn"=>ForecastArn); aws_config=aws_config)
query_forecast(Filters, ForecastArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = forecastquery("QueryForecast", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Filters"=>Filters, "ForecastArn"=>ForecastArn), args)); aws_config=aws_config)


text:
using Cassette
using Documenter

makedocs(;
    modules=[Cassette],
    authors="Jarrett Revels",
    repo="https://github.com/jrevels/Cassette.jl/blob/{commit}{path}#L{line}",
    sitename="Cassette",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jrevels.github.io/Cassette.jl",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "Disclaimers" => "disclaimers.md",
        "Why Cassette?" => "whycassette.md",
        "The Overdubbing Mechanism" => "overdub.md",
        "Contextual Dispatch" => "contextualdispatch.md",
        "Contextual Compiler Pass Injection" => "contextualpass.md",
        "Contextual Tagging of Values" => "contextualtagging.md",
        "Cassette API Documentation" => "api.md",
        "Related Work" => "relatedwork.md",
    ],
    strict=false,
)

deploydocs(;
    repo="github.com/jrevels/Cassette.jl",
)


text:
using PyPlot, JLD

data = JLD.load("./182602_23-03_qpScalingTest.jld")
resCost = data["resCost"]
resIter = data["resIter"]

cc =  ["#1f77b4" "#ff7f0e" "#2ca02c" "#d62728" "#9467bd" "#8c564b"]

PyPlot.figure(1,facecolor="white",figsize=(12,5))
PyPlot.xticks(fontsize=15)
PyPlot.yticks(fontsize=15)
PyPlot.plot(1:1:nn,resIter[1:nn,1],label="OSSDP(off)","--",color=cc[1])
PyPlot.plot(1:1:nn,resIter[1:nn,2],label="OSSDP(on)",color=cc[1])
PyPlot.plot(1:1:nn,resIter[1:nn,7],label="OSQP1(off)","--",color=cc[2])
PyPlot.plot(1:1:nn,resIter[1:nn,8],label="OSQP1(on)",color=cc[2])
# PyPlot.plot(1:1:nn,resIter[:,5],label="OSQP2(off)","--",color=cc[3])
# PyPlot.plot(1:1:nn,resIter[:,6],label="OSQP2(on)",color=cc[3])
PyPlot.grid(true)
PyPlot.xticks(1:1:nn)
PyPlot.xlabel("Problem number",fontsize=15)
PyPlot.ylabel("Iterations to Convergence",fontsize=15)
PyPlot.legend(ncol = 6,bbox_to_anchor=(0., 1.02, 1., .102),loc=3)
PyPlot.savefig("/Users/Micha/Dropbox/Research/OSSDP/Notes/Scaling_experiments/figs/qp_compare.eps")

PyPlot.figure(2,facecolor="white",figsize=(12,5))
 PyPlot.xticks(fontsize=15)
      PyPlot.yticks(fontsize=15)
PyPlot.plot(1:1:nn,resIter[1:nn,1],label=L"\rho=1.0 (off)","--",color=cc[1])
PyPlot.plot(1:1:nn,resIter[1:nn,2],label=L"\rho=1.0 (on)",color=cc[1])
PyPlot.plot(1:1:nn,resIter[1:nn,3],label=L"\rho=100.0 (off)","--",color=cc[2])
PyPlot.plot(1:1:nn,resIter[1:nn,4],label=L"\rho=100.0 (on)",color=cc[2])
PyPlot.plot(1:1:nn,resIter[1:nn,5],label=L"\rho=1000.0 (off)","--",color=cc[3])
PyPlot.plot(1:1:nn,resIter[1:nn,6],label=L"\rho=1000.0 (on)",color=cc[3])
PyPlot.grid(true)
PyPlot.xticks(1:1:nn)
PyPlot.xlabel("Problem number",fontsize=15)
PyPlot.ylabel("Iterations to Convergence",fontsize=15)
PyPlot.legend(ncol = 6,bbox_to_anchor=(0., 1.02, 1., .102),loc=3)

PyPlot.savefig("/Users/Micha/Dropbox/Research/OSSDP/Notes/Scaling_experiments/figs/qp_rho.eps")


text:
bitstring(x::Posit8) = bitstring(reinterpret(UInt8,x))
bitstring(x::Posit16) = bitstring(reinterpret(UInt16,x))
bitstring(x::Posit32) = bitstring(reinterpret(UInt32,x))

bitstring(x::Posit8_1) = bitstring(UInt8(reinterpret(UInt32,x) >> 24))
bitstring(x::Posit8_2) = bitstring(UInt8(reinterpret(UInt32,x) >> 24))

bitstring(x::Posit16_1) = bitstring(UInt16(reinterpret(UInt32,x) >> 16))
bitstring(x::Posit16_2) = bitstring(UInt16(reinterpret(UInt32,x) >> 16))

# cut off the string as UInt24 does not exist
bitstring(x::Posit24_1) = bitstring(reinterpret(UInt32,x))[1:24]
bitstring(x::Posit24_2) = bitstring(reinterpret(UInt32,x))[1:24]

"""Bitstring split into sign, regime, exponent and fraction bits.
# Example
julia> bitstring(Posit16(123)," ")
"0 11110 0 111011000"
"""
function bitstring(x::AbstractPosit,mode::String)
    if mode == " "
        # number of exponent bits
        ne = if x isa PositEx1 1 elseif x isa PositEx2 2 else 0 end
        s = bitstring(x)
        n = n_regimebits(s)

        # concatenate the string
        if n >= length(s)-2     # only exponent bits (possibly including terminating bit)
            return "$(s[1]) $(s[2:end])"
        elseif n > length(s)-ne-3
            return "$(s[1]) $(s[2:2+n]) $(s[2+n+1:end])"
        else
            if ne == 0
                return "$(s[1]) $(s[2:2+n]) $(s[2+n+1:end])"
            else
                return "$(s[1]) $(s[2:2+n]) $(s[2+n+1:2+n+ne]) $(s[2+n+ne+1:end])"
            end
        end
    else
        return bitstring(x)
    end
end

"""Number of regime bits for a posit bitstring, excluding the terminating bit.
# Example
julia> n_regimebits("0111100111011000")
4"""
function n_regimebits(s::String)
    n = 2   # skip the sign bit
    slength = length(s)
    while n < slength && s[n+1] == s[n]
        n += 1
    end
    n-1     # subtract the sign bit
end


text:
@testset "Querying" begin
    server = InfluxServer("http://localhost:8086")
    dbname = "NOAA_water_database"

    # Simple query:
    data, = query(server, dbname, SELECT(measurements=["h2o_feet"]))
    @test size(data) == (15258, 4)

    # Simple query with limits:
    for limit in (1, 5, 100)
        data, = query(server, dbname, SELECT(measurements=["h2o_pH"], limit=limit))
        @test size(data) == (limit, 3)
    end

    # Double-measurement query, showing that the columns get merged
    feet_data, pH_data = query(server, dbname, SELECT(measurements=["h2o_feet", "h2o_pH"], limit=10))
    @test size(feet_data) == (10, 5)
    @test size(pH_data) == (10, 5)

    h2o_feet_data, h2o_pH_data = query(server, dbname,
        SELECT(;
            measurements=["h2o_feet", "h2o_pH"],
            limit=10
        ),
    )
    @test size(h2o_feet_data) == (10, 5)
    @test size(h2o_pH_data) == (10, 5)

    # Test WHERE clauses, fields/tags, and sub-limit results
    h2o_feet_data = query(server, dbname,
        SELECT(;
            fields=["water_level"],
            tags=["location"],
            measurements=["h2o_feet"],
            limit=20,
            condition=@querify(:location = "santa_monica" && :water_level > 7.0),
        ),
    )[1]
    @test size(h2o_feet_data) == (13, 3)
    @test mean(h2o_feet_data[!, :water_level]) > 7.0
end


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

using BAT
using Test

using Random, Distributions, StatsBase

@testset "importance_samplers" begin
    function test_moments(dist::AnyDensityLike, algo::BAT.AbstractSamplingAlgorithm; rtol::Real=0.01)
        # ToDo: Wrap in @inferred when type stable
        samples = bat_sample(dist, algo).result

        @test isapprox(@inferred(mean(samples)), @inferred(mean(dist)), rtol=rtol)
        @test isapprox(@inferred(var(samples)), @inferred(var(dist)), rtol=rtol)
    end

    @testset "sobol_sampler" begin
        dist = MvNormal([0.4, 0.6], [2.0 1.2; 1.2 3.0])
        algo = SobolSampler(nsamples=10^5)

        test_moments(dist, algo)
    end

    @testset "grid_sampler" begin
        dist = MvNormal([0.4, 0.6], [2.0 1.2; 1.2 3.0])
        algo = GridSampler(ppa=500)

        test_moments(dist, algo, rtol=0.1)
    end
end


text:
module Godunov

include(joinpath(dirname(@__FILE__),"Limiters.jl"))

function calc_limiter(limiter_func, f, i, i⁻, i⁺)
    if f[i] - f[i⁻] ≈ 0.0
        r = 0.0
    elseif f[i⁺] - f[i] ≈ 0.0
        r = 2.0
    else
        r = (f[i] - f[i⁻])/(f[i⁺] - f[i])
    end
    return limiter_func(r)
end

function Φ_constant(f, i, i⁻, i⁻², c)
    return abs(c)*f[i⁻]
end

function generate_Φ_linear(limiter_func)
    function Φ(f, i, i⁻, i⁻², c)
        return abs(c)*(f[i⁻] + calc_limiter(limiter_func, f, i⁻, i⁻², i)*0.5*(f[i]-f[i⁻]))
    end
    return Φ
end

function generate_Φ(Riemann_solver_type, limiter_func)
    if Riemann_solver_type == :Riemann_constant
        return Φ_constant
    elseif Riemann_solver_type == :Riemann_linear
        return generate_Φ_linear(limiter_func)
    end
end

function generate_solve_pm(f₀, f, Riemann_solver_type; flux_limiter)
    if flux_limiter == :None
        limiter_func = x -> 1.0
    else
        limiter_func = Limiters.generate_limiter(flux_limiter)
    end
    Φ = generate_Φ(Riemann_solver_type, limiter_func)

    function godunov⁺!(g, f, i, i⁻, i⁺, i⁻², c)
        g[i] = f[i] + Φ(f, i, i⁻, i⁻², c) - Φ(f, i⁺, i, i⁻, c)
    end

    function godunov⁻!(g, f, i, i⁻, i⁺, i⁺², c)
        return godunov⁺!(g, f, i, i⁺, i⁻, i⁺², c)
    end

    function solve⁺!(c)
        godunov⁺!(f, f₀, 1, length(f), 2, length(f)-1, c)
        godunov⁺!(f, f₀, 2, 1, 3, length(f), c)
        @inbounds @fastmath @simd for i = 3:length(f)-1
            godunov⁺!(f, f₀, i, i-1, i+1, i-2, c)
        end
        godunov⁺!(f, f₀, length(f), length(f)-1, 1, length(f)-2, c)
    end

    function solve⁻!(c)
        godunov⁻!(f, f₀, 1, length(f), 2, 3, c)
        @inbounds @fastmath @simd for i = 2:length(f)-2
            godunov⁻!(f, f₀, i, i-1, i+1, i+2, c)
        end
        godunov⁻!(f, f₀, length(f)-1, length(f)-2, length(f), 1, c)
        godunov⁻!(f, f₀, length(f), length(f)-1, 1, 2, c)
    end

    return solve⁻!, solve⁺!
end

function generate_solver(f₀, f, Riemann_solver_type; flux_limiter=:None)
    solve⁻!, solve⁺! = generate_solve_pm(f₀, f, Riemann_solver_type; flux_limiter=flux_limiter)

    function solve!(c)
        if c > 0
            solve⁺!(c)
        else
            solve⁻!(c)
        end
    end

    return solve!
end

function generate_solver(f₀, f, c, Riemann_solver_type; flux_limiter=:None)
    solve⁻!, solve⁺! = generate_solve_pm(f₀, f, Riemann_solver_type; flux_limiter=flux_limiter)

    if c > 0
        return () -> solve⁺!(c)
    else
        return () -> solve⁻!(c)
    end
end

end # module


text:
import SpecialFunctions: gamma

abstract type AbstractWindResourceModel end

"""
    DiscritizedWindResource(wind_directions, wind_speeds, wind_probabilities, measurement_heights, air_density, ti_model, wind_shear_model)

    Struct defining a wind resource

# Arguments
- `wind_directions::Array{Float}(Nstates)`: an array of wind directions corresponding to each wind farm state
- `wind_speeds::Array{Float}(Nstates)`: an array of wind speeds corresponding to each wind farm state
- `wind_probabilities::Array{Float}(Nstates)`: an array of probabilities corresponding to each wind farm state
- `measurement_heights::Array{Float}(Nstates)`: an array of measurement heights corresponding to each wind farm state
- `air_density::Float`: the air density
- `ambient_ti::Array{Float}`: an array of the ambient turbulence intensity for each wind direction
- `wind_shear_model::Array{AbstractWindShearModel}(1)`: contains a struct defining the desired turbulence intensity model
"""
struct DiscretizedWindResource{AF, TF, ASM} <: AbstractWindResourceModel

    wind_directions::AF
    wind_speeds::AF
    wind_probabilities::AF
    measurement_heights::AF
    air_density::TF
    ambient_tis::AF
    wind_shear_model::ASM

end

function _dist_weibull(x,L;k=2.0)
    if L < 0.0001
        L = 0.0001
    end
    return (k/L)*(x/L)^(k-1)*exp(-(x/L)^k)
end

function setup_weibull_distribution(wind_directions,wind_probabilities,wind_speeds,nspeeds)
        if nspeeds == 1
            return wind_directions, wind_probabilities, wind_speeds
        end

        ndirections = length(wind_directions)
        dirs = zeros(ndirections*nspeeds)
        freqs = zeros(ndirections*nspeeds)
        speeds = zeros(ndirections*nspeeds)

        #direction loops
        k = 1
        for i = 1:ndirections
            for j = 1:nspeeds
                dirs[k] = wind_directions[i]
                k += 1
            end
        end

        #speed and frequency loops
        m = 1
        for i = 1:ndirections
            avg_speed = wind_speeds[i]
            speed_dist = range(25.0/(2.0*nspeeds)+0.001,stop=25.0-25.0/(2.0*nspeeds),length=nspeeds)
            dspeed = speed_dist[2]-speed_dist[1]
            num_int = 1000

            for j = 1:nspeeds
                speed_int = range(speed_dist[j]-dspeed/2.0,stop=speed_dist[j]+dspeed/2.0,length=num_int)
                k = 2.0
                scale = avg_speed/(gamma(1.0+1.0/k))
                freq_int = _dist_weibull.(speed_int,scale)

                speed_freq = trapz(speed_int,freq_int)
                speeds[m] = speed_dist[j]
                freqs[m] = speed_freq*wind_probabilities[i]
                m += 1
            end
        end

        return dirs, freqs, speeds
end


text:
# Autogenerated wrapper script for SCIP_jll for x86_64-apple-darwin
export libscip

using bliss_jll
using GMP_jll
using CompilerSupportLibraries_jll
using Ipopt_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("SCIP")
JLLWrappers.@declare_library_product(libscip, "@rpath/libscip.7.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(bliss_jll, GMP_jll, CompilerSupportLibraries_jll, Ipopt_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libscip,
        "lib/libscip.7.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()