text:
# using Routing
include("../src/Routing_include.jl")
using Test

# For testing purpose
using Random

for i in 1:100
rnd_seed = rand(0:1000000)
Random.seed!(rnd_seed)
# Random.seed!(3232) # old bug instance
solomon_dataset_name = "R104_050"
solomon = load_solomon(solomon_dataset_name)
n_customers = length(solomon.nodes) - 1
dual_var = rand(0:20, n_customers)

# Random.seed!(0)

# For Debugging
include("debugging.jl")

# solomon_dataset_name = "C102_050"
# solomon = load_solomon(solomon_dataset_name)
# n_customers = length(solomon.nodes) - 1
# dual_var = rand(0:20, n_customers)
# dual_var = ones(n_customers) .* 10

# dual_var = [3.8237089201877836,5.935211267605691,3.5788732394366782,7.671830985915506,6.599999999999994,4.678873239436655,11.300000000000072,16.204225352112584,5.099999999999973,12.976291079812267,6.776291079812243,12.121126760563307,6.208450704225527,12.900000000000034,18.743661971831017,0.6999999999999886,18.164788732394427,12.098356807511898,10.623708920187784,14.301291079812213,15.335211267605604,8.099999999999994,8.43521126760563,8.878873239436668,23.042253521126618,14.214084507042223,0.5999999999999943,0.7999999999999545,26.52887323943669,18.499999999999915,4.376291079812155,14.477582159624472,7.247417840375427,17.850000000000094,29.799999999999997,22.476291079812142,4.199999999999996,20.352112676056265,17.021126760563277,6.28591549295777,8.099999999999994,12.75633802816887,31.16478873239444,5.300000000000068,13.80000000000009,14.600000000000009,1.5237089201877154,16.66478873239444,39.23521126760568,14.007042253521107,15.048122065727835,10.67629107981211,8.70000000000001,15.235915492957815,6.135915492957793,11.999999999999936,6.1563380281688715,4.756338028168925,4.3352112676056045,5.400000000000006,12.200000000000017,6.9525821596244235,26.847417840375446,19.40000000000003,37.94553990610343,15.521126760563313,30.90000000000004,7.721126760563319,10.123708920187838,7.798708920187778,20.40704225352111,6.578873239436646,9.799999999999972,16.33591549295781,5.300000000000011,11.928873239436705,5.300000000000026,5.599999999999959,6.77887323943666,3.70000000000001,2.7230046948355415,9.499999999999911,11.20000000000001,12.956338028168986,2.752112676056356,20.487323943661863,16.48732394366207,6.399999999999917,2.3999999999999773,21.001291079812308,6.299999999999997,1.700000000000017,11.700000000000017,9.142253521126836,1.0000000000000036,8.435211267605553,4.900000000000006,11.26478873239446,6.587323943661858,3.299999999999983] # R201_100
# dual_var = [2.200000000000003,5.599999999999952,1.1000000000000085,8.423809523809622,1.8000000000000114,9.380555555555656,3.5,10.158333333332985,6.0999999999999375,13.735317460317347,6.599999999999966,0.29999999999996163,6.400000000000006,10.251984126984016,16.325000000000088,1.2527777777780447,11.099999999999994,20.30000000000009,11.000000000000043,15.299999999999983,1.200000000000017,2.0726190476186197,35.595634920635035,9.463888888889038,21.800793650793302,0.16865079365103952,0.09999999999996587,3.731349206349016,25.20000000000003,17.400000000000034,3.5,15.299999999999997,7.595238095237917,8.295238095237977,29.904761904761997,13.052777777777813,1.600000000000028,23.500000000000085,11.709126984126762,10.447619047618819,15.636904761905036,11.584126984126765,28.427777777777656,3.7027777777785484,22.113888888888617,20.600000000000023,14.816666666666872,3.3305555555553386,37.372222222222206,3.304761904762046,7.8333333333336554,3.799999999999997,8.8,15.828571428571502,7.363888888889036,11.51428571428557,5.372222222222355,4.9523809523811835,2.60000000000007,8.730555555555316,17.116666666666205,5.336111111111094,25.50000000000002,20.361111111111263,47.14087301587317,12.59523809523806,13.200000000000003,7.499999999999993,11.636111111111113,13.09523809523806,8.100000000000001,0.5999999999999943,11.461904761904886,16.928571428571495,0.5857142857144382,15.000000000000014,4.300000000000033,16.64007936507921,4.399999999999977,2.500000000000121,3.490476190476002,3.8999999999999773,0.7000000000000028,16.01666666666718,1.8999999999999773,19.80555555555539,13.227777777777698,10.125000000000128,1.300000000000078,16.833333333333478,3.0527777777777754,2.8000000000000003,4.846825396825898,0.30000000000000426,1.2999999999999727,0.7999999999999687,7.178174603174295,1.6999999999999957,7.727777777777462,11.753174603174118] #R202_100

# dual_var = [4.599999999999994,6.6107556618580645,6.507937812990541,5.033395620001583,12.031500859556067,5.393557067045491,11.484098213618402,7.100000000000023,7.067975932431393,5.599999999999994,4.599999999999968,6.799999999999983,4.099999999999947,10.735395022049524,14.400000000000006,1.3999999999999897,17.599999999999994,9.777191867852551,7.381822258763812,5.540791538978972,2.8646049779504636,5.763771582330683,1.5000000000000178,11.235395022049524,9.564604977950488,11.299999999999983,2.0558225577396882,5.700000000000017,20.81116301666787,13.400000000000013,5.370509754092353,13.71841692204201,11.052612302862656,11.419519396068452,14.5,14.040896180581768,3.4000000000000363,18.956734434561607,15.30615516854769,14.776186561028444,8.099999999999994,10.648471485163338,24.419624037670967,8.700000000000019,13.494622169070894,18.38678899768276,13.882678077583972,26.37821959787736,18.49460348307046,12.600000000000023,9.491232528589592,9.983866507212749,3.1099222662382644,6.335395022049545,10.499999999999996,12.883922565213926,7.263177367516152,3.9000000000000465,2.599999999999966,6.325229837805605,3.2000000000000015,6.162183272292367,13.874512295388278,19.667478884819456,33.23208386277004,17.700000000000017,32.668790642051135,1.0826818147844754,8.279254802302148,11.168248748038039,11.176683608640403,7.711163016667868,4.399405785185504,8.733395620001625,5.39999999999997,8.11116301666788,9.900000000000006,4.856667164959989,3.656667164959971,7.317318185215502,6.597604454742546,6.564604977950532,5.584632633231141,11.890126317363125,5.052784214066731,15.199999999999994,8.900000000000006,6.780701098736721,2.6000000000000085,9.24079153897899,6.300000000000023,2.194622169070948,8.407324912175962,7.7253606398087316,7.116133492787222,4.004140817699197,6.211465729875209,3.694622169070905,5.8692615292621175,4.107324912175948] # R205_100

@show length(dual_var)

pg = solomon_to_espprc(solomon, dual_var)
print_pct_neg_arcs(pg)

##################################################

@info("ESPPRC $(solomon_dataset_name) testing...")
print("Pulse     : "); @time pulse = solveESPPRC(pg, method="pulse")
print("Mono      : "); @time mono0 = solveESPPRC(pg, method="monodirectional")
print("Mono DSSR : "); @time mono1 = solveESPPRC(pg, method="monodirectional", DSSR=true)
print("Bi        : "); @time bidi0 = solveESPPRC(pg, method="bidirectional")
print("Bi   DSSR : "); @time bidi1 = solveESPPRC(pg, method="bidirectional", DSSR=true)

@show pulse.cost, pulse.load, pulse.time
@show mono0.cost, mono0.load, mono0.time
@show mono1.cost, mono1.load, mono1.time
@show bidi0.cost, bidi0.load, bidi0.time
@show bidi1.cost, bidi1.load, bidi1.time
@show pulse.path
@show mono0.path
@show mono1.path
@show bidi0.path
@show bidi1.path

# show_details(bidi0.path, pg)

@testset "ESPPRC $(solomon_dataset_name) Test" begin
    @test isapprox(pulse.cost, mono0.cost, atol=1e-7)
    @test isapprox(mono0.cost, mono1.cost, atol=1e-7)
    @test isapprox(pulse.cost, bidi0.cost, atol=1e-7)
    @test isapprox(bidi0.cost, bidi1.cost, atol=1e-7)
end
@show rnd_seed
println("done.")

end

############################################################
# max_neg = 20
# @info("ESPPRC $(solomon_dataset_name) testing with max_neg_routes=$(max_neg)...")

# @time sol, neg_sols = solveESPPRC_vrp(pg, method="pulse", max_neg_routes=max_neg)
# @time lab1, neg_labs1 = solveESPPRC_vrp(pg, method="monodirectional", max_neg_routes=max_neg)
# @time lab1d, neg_labs1d = solveESPPRC_vrp(pg, method="monodirectional", max_neg_routes=max_neg, DSSR=true)

# # @time lab2, neg_labs2 = solveESPPRC(pg, method="bidirectional", max_neg_routes=max_neg)

# @testset "ESPPRC $(solomon_dataset_name) Test with max_neg_routes" begin
#     @test length(neg_sols) <= max_neg
#     @test length(neg_labs1) <= max_neg
#     @test length(neg_labs1d) <= max_neg
#     # @test length(neg_labs2) <= max_neg
#     for l in neg_sols
#         @test l.cost < 0.0
#     end
#     for l in neg_labs1
#         @test l.cost < 0.0
#     end
#     for l in neg_labs1d
#         @test l.cost < 0.0
#     end
#     # for l in neg_labs2
#     #     @test l.cost < 0.0
#     # end
# end


text:
using NearestNeighbors
using MNIST

X_train, y_train = traindata()
X_test, y_test = testdata()

@elapsed t1 = NaiveNeighborTree(X_train)
@elapsed t2 = BallTree(X_train) # VERY SLOW!!

v = X_train[:, 1]

@elapsed k_nearest(v, t1, 30)
@elapsed k_nearest(v, t2, 30)

errors = 0
n = length(y_test)

for i in 1:25
	is, ds = k_nearest(X_test[:, i], t2, 3)
	if mode(y_train[is]) != y_test[i]
		errors += 1
	end
end

errors / 25

ext: jl


text:
using Test
using HTTP
using HTTP.IOExtras, HTTP.Sockets

@testset "websockets.jl" begin
    p = 8085 # rand(8000:8999)
    socket_type = ["wss", "ws"]

    function listen_localhost()
        @async HTTP.listen(Sockets.localhost, p) do http
            if HTTP.WebSockets.is_upgrade(http.message)
                HTTP.WebSockets.upgrade(http) do ws
                    while !eof(ws)
                        data = readavailable(ws)
                        write(ws, data)
                    end
                end
            end
        end
    end

    @testset "External Host - $s" for s in socket_type
        HTTP.WebSockets.open("$s://echo.websocket.org") do io
            write(io, "Foo")
            @test !eof(io)
            @test String(readavailable(io)) == "Foo"

            write(io, "Hello")
            write(io, " There")
            write(io, " World", "!")
            closewrite(io)

            buf = IOBuffer()
            write(buf, io)
            @test String(take!(buf)) == "Hello There World!"
        end
    end

    @testset "Localhost" begin
       listen_localhost()

        HTTP.WebSockets.open("ws://127.0.0.1:$(p)") do ws
            write(ws, "Foo")
            @test String(readavailable(ws)) == "Foo"

            write(ws, "Bar")
            @test String(readavailable(ws)) == "Bar"
        end
    end
end


text:
# by LJR
using Plots
using LaTeXStrings

fileName = ARGS

function plotFile(file::String; method::Symbol=:no, order=5)
    fit     = true
    if method == :linear
        M = 1
    elseif method == :poly
        M = 2
    elseif method == :exp
        M = 3
    elseif method == :log
        M = 4
    elseif method == :no
        fit = false
    else
        error("Method not implemented")
    end
    # Get the filename without extention name like .csv
    fileRegex = match(r"(.+?)(\.[^.]*$|$)", file)
    csv = readdlm(file, ',', Float64)
    y = csv[:, 1]; x = csv[:, 2]
    # Extent plot domain by 20%
    x_max = maximum(x)
    x_min = 0
    @show    xlen = x_max-x_min
    x_plot = 0:xlen/100:(x_max+xlen/5)
    # Plot data
    scatter(x, y, fmt = :pdf, ylabel = "Voltage (V)", xlabel = "Current (A)", label="Data")
    # Plot the best fit line
    if fit
        plot!(x_plot, dataFit(x, y, x_plot, Val{M}, order), label="Best fit line")
        savefig(fileRegex[1]*"_$method"*".pdf")
    end
    savefig(fileRegex[1]*".pdf")
end

# Dispatch into the right argument
dataFit(x, y, x_plot, ::Type{Val{1}}, kwarg...) = dataFit(x, y, x_plot, Val{1})
dataFit(x, y, x_plot, ::Type{Val{3}}, kwarg...) = dataFit(x, y, x_plot, Val{3})
dataFit(x, y, x_plot, ::Type{Val{4}}, kwarg...) = dataFit(x, y, x_plot, Val{4})

# Linear fit by using linreg
function dataFit(x, y, x_plot, ::Type{Val{1}})
    b, a = linreg(x, y)
    map(x -> a*x + b, x_plot)
end

# Polynomial fit by QR factorization
function dataFit(x, y, x_plot, ::Type{Val{2}}, order)
    function horner(coef,x)
        sum = coef[end]
        for k = length(coef)-1:-1:1
            sum = coef[k] + x*sum
        end
        sum
    end
  A = [ float(x[i])^pow for i = 1:length(x), pow = 0:order ]
  co = A \ y
  map(x -> horner(co, x), x_plot)
end

function dataFit(x, y, x_plot, ::Type{Val{3}})
    fit = linreg(x, log(y))
    a1 = exp(fit[1]); a2 = fit[2]
    map(x -> a1*exp(x*a2), x_plot)
end

function dataFit(x, y, x_plot, ::Type{Val{4}})
    a1, a2 = linreg(log(x), y)
    map(x -> a1 + a2*log(x), x_plot)
end

plotFile.(fileName)

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
function make_rules_from_one_line(rule_string::AbstractString)
    lhs, rhs = map(strip, split(rule_string, "-->"))
    rhss = map(rhs->map(strip, rhs), map(split, map(strip, split(rhs, "|"))))

    rules = map(rhs->make_cfg_rule(lhs, rhs...), rhss)

    (binary_rules   = [r.rule for r in rules if !r.isterminal],
     terminal_rules = [r.rule for r in rules if  r.isterminal])
end

function make_rules(rule_string::AbstractString)
    rules = map(make_rules_from_one_line, split(strip(rule_string), "\n"))
    (binary_rules   = vcat(map(rs -> rs.binary_rules,   rules)...),
     terminal_rules = vcat(map(rs -> rs.terminal_rules, rules)...))
end

function make_cfg_rule(lhs, rhs1, rhs2)
    r = NamedFunction(
        "$lhs --> $rhs1 $rhs2",
        @closure x -> x == lhs ? (rhs1, rhs2) : nothing)

    (rule = r, isterminal = false)
end

function make_cfg_rule(lhs, terminal)
    r = NamedFunction(
        "$lhs --> $terminal",
        @closure x -> x == lhs ? terminal : nothing)
    (rule = r, isterminal = true)
end

function grammar_from_string(str::AbstractString)
    categories = map(s -> string(strip(split(s, "-->")[1])), split(strip(str), "\n"))
    start_categories = [categories[1]]

    binary_rules, terminal_rules = make_rules(str)

    Grammar(categories, start_categories, binary_rules, terminal_rules, identity)
end


text:
@testset "Occu GOF" begin

  #Test separating cohorts by NA pattern
  Random.seed!(123)
  umd1 = simulate(Occu, @formula(ψ~elev+forest), @formula(p~wind+precip),
                  (50,3), [0,0,0,0,0,0])
  fit1 = occu(@formula(ψ~elev), @formula(p~wind), umd1)

  gc = Unmarked.get_cohorts(fit1)
  @test isequal(length(gc), 1)
  @test isequal(gc[1].y, umd1.y)

  yna = Array{Union{Int,Missing}}(deepcopy(umd1.y))
  yna[1,:] = fill(missing,3)
  yna[2,1] = missing
  yna[3,[1,3]] = fill(missing,2)
  umd2 = UmData(yna, umd1.site_covs, umd1.obs_covs)
  fit2 = occu(@formula(ψ~elev), @formula(p~wind), umd2)

  gc = Unmarked.get_cohorts(fit2)
  @test isequal(length(gc), 3)
  @test isequal(gc[1].y, umd2[4:50].y)
  @test isequal(gc[2].y, umd2[3].y)
  @test isequal(gc[3].y, umd2[2].y)

  #Test encounter history probability
  psi = 0.5
  pval = [0.2, 0.8, 0.3]
  eh1 = "0 1 1"
  ehp1 = Unmarked.get_prob_eh(Occu, eh1, pval, psi)
  @test isequal(ehp1, 0.5*prod(pval .* [0,1,1]
                               + (1 .- pval) .* (1 .- [0,1,1])))
  eh2 = "0 missing 1"
  ehp2 = Unmarked.get_prob_eh(Occu, eh2, pval, psi)
  @test isequal(ehp2, 0.5*prod([0.2,0.3] .* [0,1]
                               + (1 .- [0.2,0.3]) .* (1 .- [0,1])))

  #Test getting obs and expected counts
  oc = Unmarked.get_obs_counts(fit1, umd1)
  @test isequal(length(oc), 8)
  @test isequal(sum(collect(values(oc))), 50)
  @test isequal(oc["0 0 0"], 35)
  ec = Unmarked.get_exp_counts(fit1, umd1)
  @test isequal(keys(oc), keys(ec))
  @test isequal(round(ec["1 1 0"],digits=4), 2.2077)

  #Test chi-square calc
  @test isequal(round(Unmarked.mb_chisq(fit1),digits=4), 2.095)
  #with missing values
  @test isequal(round(Unmarked.mb_chisq(fit2),digits=4), 14.9232)

  #Test parametric bootstrap
  Random.seed!(123)
  gf = gof(fit1)
  @test isequal(length(gf.tstar), 50)
  @test isequal(mean(gf.t0 .< gf.tstar), 0.86)
  gf2 = gof(fit2)
  @test isequal(mean(gf2.t0 .< gf2.tstar), 0.12)

end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
# Dynamics Tests
using Schrodinger
using Test, LinearAlgebra, SparseArrays
println("Testing Dynamics...")

@testset "Time Dynamics" begin
N = 5 # cavity levels for tests

@testset "Simple Qubit" begin
println("  - Simple qubit (1/3)")
ω = 0.2*2π
H = 0.5ω*σx
ψ₀ = basis(2,0)
res = @inferred sesolve(H,ψ₀,(0.0,10.0),(σz,σy),saveat=0.1)
@test real.(res.evals) ≈ [cos.(ω.*res.times) sin.(π.+ω.*res.times)]
@test SchrodingerProp(H,(0,10))(ψ₀) ≈ res.states[end]
end

@testset "Jaynes-Cummings" begin
println("  - Jaynes-Cummings (2/3)")
ωc = 1.0*2π
ωa = 1.0*2π
g = 0.05*2π
κ = 0.1
a = destroy(N) ⊗ qeye(2)
sm = qeye(N) ⊗ destroy(2)
Hj = ωc*a'*a + ωa*sm'*sm + g*(a'*sm + a*sm')
O = (a'*a, sm'*sm)
ψ₀ = normalize!(basis(N,1)+0.5*basis(N,3)+1im*basis(N,4)) ⊗ basis(2,0)
res = @inferred sesolve(Hj, ψ₀, (0.0,25.0), qeye(N) ⊗ Operator([0 0;0 1]))
f(t,g) = 0.5*(1 - (4/9*cos(2g*t) + 1/9*cos(√3*2g*t) + 4/9*cos(√4*2g*t)))
@test real.(res.evals) ≈ f.(res.times,g)
@test SchrodingerProp(Hj,(0.0,25.0))(ψ₀) ≈ res.states[end]

ψ₀ = basis(N,0) ⊗ basis(2,1)
c_ops = (√(κ)*a, √(κ)*sm) # use same rate to solve eq. analytically
res = @inferred mesolve(Hj,c_ops,Operator(ψ₀),(0.0,25.0),O)
res = mesolve(Hj,c_ops,ψ₀,(0.0,25.0),O)
f(t,g,κ,ϕ) = 0.5*(1+cos(t*2g+ϕ))*exp(-0.5κ*t)*exp(-0.5κ*t)
@test real.(res.evals) ≈ [f.(res.times,g,κ,-π) f.(res.times,g,κ,0)]
@test LindbladProp(Hj,c_ops,(0.0,25.0))(Operator(ψ₀)) ≈ res.states[end]
L = @inferred LindbladEvo(Hj,c_ops)
@test LindbladProp(L,(0.0,25.0))(Operator(ψ₀)) ≈ res.states[end]
end

@testset "Three-Level Atom and Cavity" begin
println("  - Three-Level Atom and Cavity (3/3)")
c = 4.4 # coupling strength
κ = 0.1
g,e,u = basis(3,0), basis(3,1), basis(3,2)
σ_ge = qeye(N) ⊗ (g*e') # |g><e|
σ_ue = qeye(N) ⊗ (u*e') # |u><e|
a = destroy(N) ⊗ qeye(3)
n = numberop(N) ⊗ qeye(3)
ψ₀ = normalize!(basis(N,0) + 0.5*basis(N,1)) ⊗ u # Define initial state
σ_gg = Operator(basis(N,1) ⊗ g) # Define states onto which to project
σ_uu = Operator(basis(N,0) ⊗ u)
H0 = n - c * (σ_ge' * a + a' * σ_ge) # time-independent term
H1 = (σ_ue' + σ_ue);  # time-dependent term
H1_coeff(t,p) = 9*exp(-(t/5)^2)
S = @inferred SchrodingerEvo(H0,(H1,H1_coeff))
res = @inferred lsolve(S, ψ₀, (-15.0,15.0), (n, σ_uu, σ_gg), Schrodinger.Tsit5())
U1 = @inferred SchrodingerProp(H0,(H1,H1_coeff),(-15.0,15.0),2000)
U2 = @inferred SchrodingerProp(S,(-15.0,15.0))
@test U1(ψ₀) ≈ res.states[end] rtol=5E-6 # okay showing here...
@test U2(ψ₀) ≈ res.states[end] # this one is much better, need to figure out why
@test U2 ≈ U1 rtol=4E-5

L = @inferred LindbladEvo((H0,(H1,H1_coeff,[1,2,3])),√(κ)*a)
res = @inferred lsolve(L, Operator(ψ₀), (-15.0,15.0), (n, σ_uu, σ_gg), Schrodinger.Tsit5())
U1 = @inferred LindbladProp(H0,(H1,H1_coeff),√(κ)*a,(-15.0,15.0),500)
U2 = @inferred LindbladProp(L,(-15.0,15.0))
@test U1(Operator(ψ₀)) ≈ res.states[end] rtol=4E-4 # very poor showing here...
@test U2(Operator(ψ₀)) ≈ res.states[end] # this one is much better, need to figure out why
@test U2 ≈ U1 rtol=3E-3 #...
end

end


text:
"""
    cartesian_to_lat_lon(x, y, z)

Convert 3D coordinates (x, y, z) from the unit sphere to (lat, lon). Assumes "lat" is positive
with "z", equatorial plane falls at z=0  and "lon" is measured anti-clockwise (eastward)
from x-axis (y=0) about z-axis.

This is a Julia translation of MATLAB code from MITgcm [1].

[1]: http://wwwcvs.mitgcm.org/viewvc/MITgcm/MITgcm_contrib/high_res_cube/matlab-grid-generator/map_xyz2lonlat.m?view=markup
"""
cartesian_to_lat_lon(x, y, z) = cartesian_to_latitude(x, y, z), cartesian_to_longitude(x, y, z)

cartesian_to_latitude(x, y, z) = atand(z, hypot(x, y))

cartesian_to_longitude(x, y, z) = atand(y, x)


text:
# ---
# title: 475. Heaters
# id: problem475
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search
# link: <https://leetcode.com/problems/heaters/description/>
# hidden: true
# ---
#
# Winter is coming! During the contest, your first job is to design a standard
# heater with a fixed warm radius to warm all the houses.
#
# Every house can be warmed, as long as the house is within the heater's warm
# radius range.
#
# Given the positions of `houses` and `heaters` on a horizontal line, return
# _the minimum radius standard of heaters  so that those heaters could cover all
# houses._
#
# **Notice** that all the `heaters` follow your radius standard, and the warm
# radius will the same.
#
#
#
# **Example 1:**
#
#
#
#     Input: houses = [1,2,3], heaters = [2]
#     Output: 1
#     Explanation: The only heater was placed in the position 2, and if we use the radius 1 standard, then all the houses can be warmed.
#
#
# **Example 2:**
#
#
#
#     Input: houses = [1,2,3,4], heaters = [1,4]
#     Output: 1
#     Explanation: The two heater was placed in the position 1 and 4. We need to use radius 1 standard, then all the houses can be warmed.
#
#
# **Example 3:**
#
#
#
#     Input: houses = [1,5], heaters = [2]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= houses.length, heaters.length <= 3 * 104`
#   * `1 <= houses[i], heaters[i] <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
function HPDE_∂F∂x(func_F,u,Δx,dim)
    Uₓᴸ, Uₓᴿ = GeneralFluxLimiter(u,Δx)
    n = size(u,dim)
    ∂f∂x = copy(u)

    for i=1:n
        A⁺ , A⁻ = LRJacobian(func_F,u[:,i])
        ∂f∂x[:,i] = (A⁻ * Uₓᴸ[:,i] .+ A⁺ * Uₓᴿ[:,i])
    end
    return ∂f∂x
end

function HPDE_∂u∂t(prob::HPDEProblem,u)

    ∂F∂x = HPDE_∂F∂x(prob.func_F,u,prob.Grid.Δx,2)

    # future extentions
    #∂F∂x = HPDE_∂F∂x(prob.func_F,u,prob.Grid.Δx,1)
    #∂G∂y = HPDE_∂F∂x(prob.func_G,u,prob.Grid.Δy,2)
    #∂H∂z = HPDE_∂F∂x(prob.func_H,u,prob.Grid.Δz,3)
    #S = prob.func_S(u)

    ∂G∂y = 0
    ∂H∂z = 0
    #S = func_S(u)
    S = 0

    ∂u∂t = S .- (∂F∂x .+ ∂G∂y .+ ∂H∂z)

    return ∂u∂t
end

function LRJacobian(func_F,u)
    A = ForwardDiff.jacobian(func_F, u)
    #A = (abs.(A) .> 1e-10) .* A
    #println(A)
    Λ = diagm(eigvals(A))
    Λ⁺ = 0.5 .* (Λ .+ abs.(Λ))
    Λ⁻ = 0.5 .* (Λ .- abs.(Λ))
    P = eigvecs(A)
    P⁻¹ = inv(P)
    A⁺ = P * Λ⁺ * P⁻¹
    A⁻ = P * Λ⁻ * P⁻¹
    return A⁺, A⁻
end

function GeneralFluxLimiter(u,Δx)
    ϵ = 1.0e-8
    ∂u∂x(U⁺,U⁻) = (U⁺ .- U⁻) ./ Δx
    func_Ψ(r) = 0
    #func_Ψ(r) = κ_Scheme_new(r,1)
    func_Ψ(r) = SuperBee(r)
    #func_Ψ(r) = Koren(r)
    #func_Ψ(r) = 0
    func_r(Uₓ₁,Uₓ₂) = (Uₓ₁ .+ ϵ) ./ (Uₓ₂ .+ ϵ)
    func_U(Δxᵢ,r,uₓ) = (Δxᵢ / 2) .* func_Ψ(r) .* (uₓ .+ ϵ)
    Uᵢ₋₂ = [u[:,1] u[:,1] u[:,1:end-2]]
    Uᵢ₋₁ = [u[:,1] u[:,1:end-1]]
    Uᵢ   = copy(u)
    Uᵢ₊₁ = [u[:,2:end] u[:,end]]
    Uᵢ₊₂ = [u[:,3:end] u[:,end] u[:,end]]

    Uₓi₊¾ = ∂u∂x(Uᵢ₊₂,Uᵢ₊₁)
    Uₓi₊½ = ∂u∂x(Uᵢ₊₁,Uᵢ)
    Uₓi₋½ = ∂u∂x(Uᵢ,Uᵢ₋₁)
    Uₓi₋¾ = ∂u∂x(Uᵢ₋₁,Uᵢ₋₂)

    Uᴸi₊½ = Uᵢ₊₁ .- func_U(Δx,func_r(Uₓi₊½,Uₓi₊¾),Uₓi₊¾)
    Uᴸi₋½ = Uᵢ   .- func_U(Δx,func_r(Uₓi₋½,Uₓi₊½),Uₓi₊½)
    Uᴿi₊½ = Uᵢ   .+ func_U(Δx,func_r(Uₓi₊½,Uₓi₋½),Uₓi₋½)
    Uᴿi₋½ = Uᵢ₋₁ .+ func_U(Δx,func_r(Uₓi₋½,Uₓi₋¾),Uₓi₋¾)

    Uₓᴸ = (Uᴸi₊½ .- Uᴸi₋½) ./ Δx
    Uₓᴿ = (Uᴿi₊½ .- Uᴿi₋½) ./ Δx
    #return Uᴿi₊½, Uᴿi₋½, Uᴸi₊½ , Uᴸi₋½
    #return Uₓi₊¾, Uₓi₊½, Uₓi₋½, Uₓi₋¾
    return Uₓᴸ, Uₓᴿ
end

################## METHODS ####################

function Koren(r)
    return max.(0,min.(2 .* r,min.(1/3 .+ 2 /3 .* r,2)))
end

function UpWind_FO(r)
    return 0
end

function UpWind_SO(r)
    return κ_Scheme_new(r,-1)
end

function SCD(r)
    return κ_Scheme_new(r,1)
end

function UpWind_Q(r)
    return κ_Scheme_new(r,0.5)
end

function UpWind_C(r)
    return κ_Scheme_new(r,1/3.0)
end

function SuperBee(r)
    return max.(0,min.(r,1))
end


text:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
# Given some model parameters, synthesize and plot the channel maps and
# integrated spectrum

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    # default = 0
    "--norad"
    help = "Use the image alread here."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    required = true
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using constants
using image
using emodel
using HDF5

# read the wavelengths for all 23 channels
fid = h5open("../" * config["data_file"], "r")
lams = read(fid["lams"]) # [μm]
close(fid)

pp = config["parameters"]
params = ["M_star", "a_c", "T_10", "q", "gamma", "logM_CO", "ksi", "dpc", "incl", "PA", "e", "w", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
starting_param = Array(Float64, nparam)

for i=1:nparam
    starting_param[i] = pp[params[i]][1]
end

# Convert logM_CO to M_CO
starting_param[6] = 10^starting_param[6]

pars = Parameters(starting_param...)

vel = pars.vel # [km/s]
# RADMC conventions for inclination and PA
incl = pars.incl # [deg]
PA = pars.PA # [deg] Position angle runs counter clockwise
npix = config["npix"] # number of pixels

# Doppler shift the dataset wavelength to rest-frame wavelength
beta = vel/c_kms # relativistic Doppler formula
shift_lams =  lams .* sqrt((1. - beta) / (1. + beta)) # [microns]

grd = config["grid"]
grid = Grid(grd["nr"], grd["ntheta"], grd["nphi"], grd["r_in"], grd["r_out"], grd["na"], true)

write_grid("", grid)
write_model(pars, "", grid)
write_lambda(shift_lams, "")

println("Write model")
tic()
write_model(pars, "", grid)
toc()

println("Write lams")
tic()
write_lambda(shift_lams, "")
toc()

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
@inline nstatic(::Val{N}) where {N} = ntuple(StaticInt, Val(N))

invariant_permutation(::Any, ::Any) = False()
function invariant_permutation(x::T, y::T) where {N,T<:Tuple{Vararg{StaticInt,N}}}
    if x === nstatic(Val(N))
        return True()
    else
        return False()
    end
end

permute(x::Tuple, perm::Val) = permute(x, static(perm))
permute(x::Tuple{Vararg{Any}}, perm::Tuple{Vararg{StaticInt}}) = eachop(getindex, perm, x)
function permute(x::Tuple{Vararg{Any,K}}, perm::Tuple{Vararg{StaticInt,K}}) where {K}
    if invariant_permutation(perm, perm) === False()
        return eachop(getindex, perm, x)
    else
        return x
    end
end

"""
    eachop(op, args...; iterator::Tuple{Vararg{StaticInt}}) -> Tuple

Produces a tuple of `(op(args..., iterator[1]), op(args..., iterator[2]),...)`.
"""
@inline function eachop(op::F, itr::Tuple{I1,I2,Vararg}, args::Vararg{Any,K}) where {F,I1,I2,K}
    return (op(args..., first(itr)), eachop(op, Base.tail(itr), args...)...)
end
@inline function eachop(op::F, itr::Tuple{I}, args::Vararg{Any,K}) where {F,I,K}
    return (op(args..., first(itr)),)
end

"""
    eachop_tuple(op, arg, args...; iterator::Tuple{Vararg{StaticInt}}) -> Type{Tuple}

Produces a tuple type of `Tuple{op(arg, args..., iterator[1]), op(arg, args..., iterator[2]),...}`.
Note that if one of the arguments passed to `op` is a `Tuple` type then it should be the first argument
instead of one of the trailing arguments, ensuring type inference of each element of the tuple.
"""
eachop_tuple(op, itr, arg, args...) = _eachop_tuple(op, itr, arg, args)
@generated function _eachop_tuple(op, ::I, arg, args::A) where {A,I}
    t = Expr(:curly, Tuple)
    narg = length(A.parameters)
    for p in I.parameters
        call_expr = Expr(:call, :op, :arg)
        if narg > 0
            for i in 1:narg
                push!(call_expr.args, :(getfield(args, $i)))
            end
        end
        push!(call_expr.args, :(StaticInt{$(p.parameters[1])}()))
        push!(t.args, call_expr)
    end
    Expr(:block, Expr(:meta, :inline), t)
end

#=
    find_first_eq(x, collection::Tuple)

Finds the position in the tuple `collection` that is exactly equal (i.e. `===`) to `x`.
If `x` and `collection` are static (`is_static`) and `x` is in `collection` then the return
value is a `StaticInt`.
=#
@generated function find_first_eq(x::X, itr::I) where {X,N,I<:Tuple{Vararg{Any,N}}}
    if (is_static(X) & is_static(I)) === True()
        return Expr(:block, Expr(:meta, :inline),
            :(Base.Cartesian.@nif $(N + 1) d->(x === getfield(itr, d)) d->(static(d)) d->(nothing)))
    else
        return Expr(:block, Expr(:meta, :inline),
            :(Base.Cartesian.@nif $(N + 1) d->(x === getfield(itr, d)) d->(d) d->(nothing)))
    end
end


text:
"""
	Cart-pole
"""

struct Cartpole{I, T} <: Model{I, T}
	n::Int
    m::Int
	d::Int

    mc     # mass of the cart in kg
    mp     # mass of the pole (point mass at the end) in kg
    l      # length of the pole in m
    g      # gravity m/s^2
end

function f(model::Cartpole, x, u, w)
    H = @SMatrix [model.mc + model.mp model.mp * model.l * cos(x[2]);
				  model.mp * model.l * cos(x[2]) model.mp * model.l^2.0]
    C = @SMatrix [0.0 -1.0 * model.mp * x[4] * model.l * sin(x[2]);
	 			  0.0 0.0]
    G = @SVector [0.0,
				  model.mp * model.g * model.l * sin(x[2])]
    B = @SVector [1.0,
				  0.0]
    qdd = SVector{2}(-H \ (C * view(x, 3:4) + G - B * u[1]))

    return @SVector [x[3],
					 x[4],
					 qdd[1],
					 qdd[2]]
end

n, m, d = 4, 1, 4
model = Cartpole{Midpoint, FixedTime}(n, m, d, 1.0, 0.2, 0.5, 9.81)

"""
	Cart-pole with Coulomb friction
"""
struct CartpoleFriction{I, T} <: Model{I, T}
	n::Int
    m::Int
	d::Int

    mc     # mass of the cart in kg
    mp     # mass of the pole (point mass at the end) in kg
    l      # length of the pole in m
    g      # gravity m/s^2
	μ      # coefficient of friction
end

function f(model::CartpoleFriction, x, u, w)
    H = @SMatrix [model.mc + model.mp model.mp * model.l * cos(x[2]);
				  model.mp * model.l * cos(x[2]) model.mp * model.l^2.0]
    C = @SMatrix [0.0 -1.0 * model.mp * x[4] * model.l * sin(x[2]);
	 			  0.0 0.0]
    G = @SVector [0.0,
				  model.mp * model.g * model.l * sin(x[2])]
    B = @SVector [1.0,
				  0.0]
    qdd = SVector{2}(-H \ (C * view(x, 3:4) + G - B * (u[1] + u[2] - u[3])))

    return @SVector [x[3],
					 x[4],
					 qdd[1],
					 qdd[2]]
end

n, m, d = 4, 7, 4
model_friction = CartpoleFriction{Midpoint, FixedTime}(n, m, d, 1.0, 0.2, 0.5, 9.81, 0.1)

function visualize!(vis, model, q;
       Δt = 0.1,
	   color = RGBA(1,0,0,1.0))

	default_background!(vis)

	l2 = Cylinder(Point3f0(-model.l * 2.0, 0.0, 0.0),
		Point3f0(model.l * 2.0, 0.0, 0.0),
		convert(Float32, 0.025))

	setobject!(vis["slider"], l2, MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    l1 = Cylinder(Point3f0(0.0, 0.0, 0.0),
		Point3f0(0.0, 0.0, model.l),
		convert(Float32, 0.025))

    setobject!(vis["arm"], l1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    setobject!(vis["base"], HyperSphere(Point3f0(0.0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = color))

    setobject!(vis["ee"], HyperSphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = color))

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	for t = 1:length(q)
	    MeshCat.atframe(anim,t) do
			x = q[t]
			px = x[1] + model.l * sin(x[2])
			pz = -model.l * cos(x[2])
            settransform!(vis["arm"], cable_transform([x[1]; 0;0], [px; 0.0; pz]))
            settransform!(vis["base"], Translation([x[1]; 0.0; 0.0]))
            settransform!(vis["ee"], Translation([px; 0.0; pz]))
	    end
	end

	settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(- pi / 2))))

	MeshCat.setanimation!(vis,anim)
end


text:
# Example 3.35 and 3.38 of
# Blekherman, G., Parrilo, P. A., & Thomas, R. R. (Eds.).
# Semidefinite optimization and convex algebraic geometry SIAM 2013
@testset "Example 3.35 with $(factory.optimizer_constructor)" for factory in sdp_factories
    @polyvar x
    m = SOSModel(factory)
    @constraint m x^4+4x^3+6x^2+4x+5 in SOSCone() # equivalent to >= 0
    JuMP.optimize!(m)
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    # p(x) = (x^2+2x)^2 + 2(1+x)^2 + 3
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
module ParticleDA

using Random
using Distributions, Statistics, MPI, Base.Threads, YAML, HDF5
using TimerOutputs
import Future
using EllipsisNotation

export run_particle_filter, BootstrapFilter, OptimalFilter

include("params.jl")
include("io.jl")
include("OptimalFilter.jl")

using .Default_params

# Functions to extend in the model

"""
    ParticleDA.get_grid_size(model_data) -> NTuple{N, Int} where N

Return a tuple with the dimensions (number of nodes) of the grid.
"""
function get_grid_size end
"""
    ParticleDA.get_grid_domain_size(model_data) -> NTuple{N, float} where N

Return a tuple with the dimensions (metres) of the grid domain.
"""
function get_grid_domain_size end
"""
    ParticleDA.get_grid_cell_size(model_data) -> NTuple{N, float} where N

Return a tuple with the dimensions (metres) of a grid cell.
"""
function get_grid_cell_size end

"""
    ParticleDA.get_n_state_var(model_data) -> Int

Return the number of state variables.
"""
function get_n_state_var end

"""
    ParticleDa.get_obs_noise_std(model_data) -> Float

Return standard deviation of observation noise. Required for optimal filter only.
"""
function get_obs_noise_std end

"""
    ParticleDa.get_model_noise(model_data) -> NamedTuple({:sigma, :lambda, :nu},Tuple{Float, Float, Float})

Return parameters of the noise covariance function of the first state variable (height). Required for optimal filter only.
"""
function get_model_noise_params end

"""
    ParticleDA.get_particle(model_data) -> particles

Return the vector of particles.  This method is intended to be extended by the
user with the above signature, specifying the type of `model_data`.  Note: this
function should return the vector of particles itself and not a copy, because it
will be modified in-place.
"""
function get_particles end

"""
   ParticleDA.set_particles!(model_data, particles)

Overwrite particle state in model_data with the vector particles. This method is intended to be extended by the
user with the above signature, specifying the type of `model_data`.
"""
function set_particles! end

"""
    ParticleDA.get_truth(model_data) -> truth_observations

Return the vector of true observations.  This method is intended to be extended
by the user with the above signature, specifying the type of `model_data`.
"""
function get_truth end

"""
    ParticleDA.get_stations(model_data) -> NamedTuple({:nst,:ist,:jst},Tuple{Int, Array{Float,1}, Array{Float,1}})

Return a named tuple with number of stations and their coordinates (nst,ist,jst) of the
points of observation. This method is intended to be extended by the user with the above signature,
specifying the type of `model_data`.
Required for optimal filter only.
"""
function get_stations end

"""
    ParticleDA.update_truth!(model_data, nprt_per_rank::Int) -> truth_observations

Update the true observations using the dynamic of the model and return the
vector of the true observations.  `nprt_per_rank` is the number of particles per
each MPI rank.  This method is intended to be extended by the user with the
above signature, specifying the type of `model_data`.
"""
function update_truth! end

"""
    ParticleDA.update_particle_dynamics!(model_data, nprt_per_rank::Int)

Update the particles using the dynamic of the model.  `nprt_per_rank` is the
number of particles per each MPI rank.  This method is intended to be extended
by the user with the above signature, specifying the type of `model_data`.
"""
function update_particle_dynamics! end

"""
    ParticleDA.update_particle_noise!(model_data, nprt_per_rank::Int)

Update the particles using the noise of the model and return the vector of the
particles.  `nprt_per_rank` is the number of particles per each MPI rank.  This
method is intended to be extended by the user with the above signature,
specifying the type of `model_data`.
"""
function update_particle_noise! end

"""
    ParticleDA.get_particle_observations!(model_data, nprt_per_rank::Int) -> particles_observations

Return the vector of the particles observations.  `nprt_per_rank` is the number
of particles per each MPI rank.  This method is intended to be extended by the
user with the above signature, specifying the type of `model_data`.
"""
function get_particle_observations! end

"""
    ParticleDA.write_snapshot(output_filename, model_data, avg_arr, var_arr, weights, it)

Write a snapshot of the data after an update of the particles to the HDF5 file
`output_filename`.  `avg_arr` is the array of the mean of the particles,
`var_arr` is the array of the standard deviation of the particles, `weights` is
the array of the weigths of the particles, `it` is the index of the time step
(`it==0` is the initial state, before moving forward the model for the first
time).  This method is intended to be extended by the user with the above
signature, specifying the type of `model_data`.
"""
function write_snapshot end

"""
    ParticleFilter

Abstract type for the particle filter to use.  Currently used subtypes are:
* [`BootstrapFilter`](@ref)
* [`OptimalFilter`](@ref)
"""
abstract type ParticleFilter end
"""
    BootstrapFilter()

Instantiate the singleton type `BootstrapFilter`.  This can be used as argument
of [`run_particle_filter`](@ref) to select the bootstrap filter.
"""
struct BootstrapFilter <: ParticleFilter end
struct OptimalFilter <: ParticleFilter end

# Get weights for particles by evaluating the probability of the observations predicted by the model
# from independent normal pdfs for each observation.
function get_log_weights!(weight::AbstractVector{T},
                          obs::AbstractVector{T},
                          obs_model::AbstractMatrix{T},
                          weight_std::T) where T

    nobs = size(obs_model,1)
    @assert size(obs,1) == nobs

    weight .= 1.0

    for iobs = 1:nobs
        weight .+= logpdf.(Normal(obs[iobs], weight_std), @view(obs_model[iobs,:]))
    end

end

# Get weights for particles by evaluating the probability of the observations predicted by the model
# from a multivariate normal pdf with mean equal to real observations and covariance equal to cov_obs
function get_log_weights!(weight::AbstractVector{T},
                          obs::AbstractVector{T},
                          obs_model::AbstractMatrix{T},
                          cov_obs::AbstractMatrix{T}) where T

    weight .= Distributions.logpdf(Distributions.MvNormal(obs, cov_obs), obs_model)

end

#
function normalized_exp!(weight::AbstractVector)

    weight .-= maximum(weight)
    @. weight = exp(weight)
    weight ./= sum(weight)

end

# Resample particles from given weights using Stochastic Universal Sampling
function resample!(resampled_indices::AbstractVector{Int}, weight::AbstractVector{T}, rng::Random.AbstractRNG=Random.default_rng()) where T

    nprt = length(weight)
    nprt_inv = 1.0 / nprt
    k = 1

    #TODO: Do we need to sort state by weight here?

    weight_cdf = cumsum(weight)
    u0 = nprt_inv * rand(rng, T)

    # Note: To parallelise this loop, updates to k and u have to be atomic.
    # TODO: search for better parallel implementations
    for ip in 1:nprt

        u = u0 + (ip - 1) * nprt_inv

        while(u > weight_cdf[k])
            k += 1
        end

        resampled_indices[ip] = k

    end

end

struct SummaryStat{T}
    avg::T
    var::T
    n::Int
end

function SummaryStat(X::AbstractVector)
    m = mean(X)
    v = varm(X,m, corrected=true)
    n = length(X)
    SummaryStat(m,v,n)
end

function stats_reduction(S1::SummaryStat, S2::SummaryStat)

    n = S1.n + S2.n
    m = (S1.avg*S1.n + S2.avg*S2.n) / n

    # Calculate pooled unbiased sample variance of two groups. From https://stats.stackexchange.com/q/384951
    # Can be found in https://www.tandfonline.com/doi/abs/10.1080/00031305.2014.966589
    # To get the uncorrected variance, use
    # v = (S1.n * (S1.var + S1.avg * (S1.avg-m)) + S2.n * (S2.var + S2.avg * (S2.avg-m)))/n
    v = ((S1.n-1) * S1.var + (S2.n-1) * S2.var + S1.n*S2.n/n * (S2.avg - S1.avg)^2 )/(n-1)

    SummaryStat(m, v, n)

end

function get_mean_and_var!(statistics::Array{SummaryStat{T}},
                           particles::AbstractArray{T},
                           master_rank::Int) where T

    ndims(particles) != ndims(statistics) + 1 &&
        error("particles must have one dimension more than statistics")
    Threads.@threads for idx in CartesianIndices(statistics)
        statistics[idx] = SummaryStat(@view(particles[idx,:]))
    end

    MPI.Reduce!(statistics, stats_reduction, master_rank, MPI.COMM_WORLD)

end

function unpack_statistics!(avg::AbstractArray{T}, var::AbstractArray{T}, statistics::AbstractArray{SummaryStat{T}}) where T

    for idx in CartesianIndices(statistics)
        avg[idx] = statistics[idx].avg
        var[idx] = statistics[idx].var
    end
end

function copy_states!(particles::AbstractArray{T},
                      buffer::AbstractArray{T},
                      resampling_indices::Vector{Int},
                      my_rank::Int,
                      nprt_per_rank::Int) where T

    # These are the particle indices stored on this rank
    particles_have = my_rank * nprt_per_rank + 1:(my_rank + 1) * nprt_per_rank

    # These are the particle indices this rank should have after resampling
    particles_want = resampling_indices[particles_have]

    # These are the ranks that have the particles this rank should have
    rank_has = floor.(Int, (particles_want .- 1) / nprt_per_rank)

    # We could work out how many sends and receives we have to do and allocate
    # this appropriately but, lazy
    reqs = Vector{MPI.Request}(undef, 0)

    # Send particles to processes that want them
    for (k,id) in enumerate(resampling_indices)
        rank_wants = floor(Int, (k - 1) / nprt_per_rank)
        if id in particles_have && rank_wants != my_rank
            local_id = id - my_rank * nprt_per_rank
            req = MPI.Isend(@view(particles[..,local_id]), rank_wants, id, MPI.COMM_WORLD)
            push!(reqs, req)
        end
    end

    # Receive particles this rank wants from ranks that have them
    # If I already have them, just do a local copy
    # Receive into a buffer so we dont accidentally overwrite stuff
    for (k,proc,id) in zip(1:nprt_per_rank, rank_has, particles_want)
        if proc == my_rank
            local_id = id - my_rank * nprt_per_rank
            @view(buffer[..,k]) .= @view(particles[..,local_id])
        else
            req = MPI.Irecv!(@view(buffer[..,k]), proc, id, MPI.COMM_WORLD)
            push!(reqs,req)
        end
    end

    # Wait for all comms to complete
    MPI.Waitall!(reqs)

    particles .= buffer

end

# Initialize arrays used by the filter
function init_filter(filter_params::FilterParameters, model_data, nprt_per_rank::Int, ::Vector{<:Random.AbstractRNG}, T::Type, ::BootstrapFilter)

    if MPI.Comm_rank(MPI.COMM_WORLD) == filter_params.master_rank
        weights = Vector{T}(undef, filter_params.nprt)
    else
        weights = Vector{T}(undef, nprt_per_rank)
    end

    resampling_indices = Vector{Int}(undef, filter_params.nprt)

    size = get_grid_size(model_data)
    n_state_var = get_n_state_var(model_data)

    statistics = Array{SummaryStat{T}, length(size) + 1}(undef, size..., n_state_var)
    avg_arr = Array{T, length(size) + 1}(undef, size..., n_state_var)
    var_arr = Array{T, length(size) + 1}(undef, size..., n_state_var)

    # Memory buffer used during copy of the states
    copy_buffer = Array{T, length(size) + 2}(undef, size..., n_state_var, nprt_per_rank)

    return (;weights, resampling_indices, statistics, avg_arr, var_arr, copy_buffer)
end

# Initialize arrays used by the filter
function init_filter(filter_params::FilterParameters, model_data, nprt_per_rank::Int, rng::Vector{<:Random.AbstractRNG}, T::Type, ::OptimalFilter)

    filter_data = init_filter(filter_params, model_data, nprt_per_rank, rng, T, BootstrapFilter())

    stations = get_stations(model_data)
    size = get_grid_size(model_data)
    domain_size = get_grid_domain_size(model_data)
    cell_size = get_grid_cell_size(model_data)

    grid = Grid(size...,cell_size...,domain_size...)
    grid_ext = Grid((grid.nx-1)*2, (grid.ny-1)*2, grid.dx, grid.dy, (grid.x_length-grid.dx)*2, (grid.y_length-grid.dy)*2)

    model_noise_params = get_model_noise_params(model_data)
    obs_noise_std = get_obs_noise_std(model_data)
    # Precompute two FFT plans, one in-place and the other out-of-place
    C = complex(T)
    tmp_array = Matrix{C}(undef, grid_ext.nx, grid_ext.ny)
    fft_plan, fft_plan! = FFTW.plan_fft(tmp_array), FFTW.plan_fft!(tmp_array)

    offline_matrices = init_offline_matrices(grid, grid_ext, stations, model_noise_params, obs_noise_std, fft_plan, fft_plan!, T)
    online_matrices = init_online_matrices(grid, grid_ext, stations, nprt_per_rank, T)

    return (; filter_data..., offline_matrices, online_matrices, stations, grid, grid_ext, rng, obs_noise_std, fft_plan, fft_plan!)
end

function update_particle_proposal!(model_data, filter_data, truth_observations, nprt_per_rank, filter_type::BootstrapFilter)

    update_particle_noise!(model_data, nprt_per_rank)

end

function update_particle_proposal!(model_data, filter_data, truth_observations, nprt_per_rank, filter_type::OptimalFilter)

        # Optimal Filter: After updating the particle dynamics, we apply the "optimal proposal" in
        #                 sample_height_proposal!() to the first state variable (height). We apply
        #                 a sample from the gaussian random field in update_particle_noise!() to the other
        #                 state variables (velocity).

        particles = get_particles(model_data)

        # Apply optimal proposal, the result will be in offline_matrices.samples
        sample_height_proposal!(@view(particles[:,:,1,:]),
                                filter_data.offline_matrices,
                                filter_data.online_matrices,
                                truth_observations,
                                filter_data.stations,
                                filter_data.grid,
                                filter_data.grid_ext,
                                filter_data.fft_plan,
                                filter_data.fft_plan!,
                                nprt_per_rank,
                                filter_data.rng[threadid()],
                                filter_data.obs_noise_std)

        # Add noise from the standard gaussian random field to all state variables in model_data
        update_particle_noise!(model_data, nprt_per_rank)

        # Overwrite the height state variable with the samples of the optimal proposal
        particles[:,:,1,:] .= filter_data.online_matrices.samples
        set_particles!(model_data, particles)

end

function run_particle_filter(init, filter_params::FilterParameters, model_params_dict::Dict, filter_type; rng::Union{Random.AbstractRNG,Nothing}=nothing)

    MPI.Init()

    my_rank = MPI.Comm_rank(MPI.COMM_WORLD)
    my_size = MPI.Comm_size(MPI.COMM_WORLD)

    # For now, assume that the particles can be evenly divided between ranks
    @assert mod(filter_params.nprt, my_size) == 0

    nprt_per_rank = Int(filter_params.nprt / my_size)

    if filter_params.enable_timers
        TimerOutputs.enable_debug_timings(ParticleDA)
    end
    timer = TimerOutput()

    nprt_per_rank = Int(filter_params.nprt / MPI.Comm_size(MPI.COMM_WORLD))

    _rng = let
        m = if isnothing(rng)
            Random.MersenneTwister(filter_params.random_seed + my_rank)
        else
            rng
        end
        [m; accumulate(Future.randjump, fill(big(10)^20, nthreads()-1), init=m)]
    end

    # Do memory allocations
    @timeit_debug timer "Model initialization" model_data = init(model_params_dict, nprt_per_rank, my_rank, _rng)

    # TODO: put the body of this block in a function
    @timeit_debug timer "Filter initialization" filter_data = init_filter(filter_params, model_data, nprt_per_rank, _rng, Float64, filter_type)

    @timeit_debug timer "get_particles" particles = get_particles(model_data)
    @timeit_debug timer "Mean and Var" get_mean_and_var!(filter_data.statistics, particles, filter_params.master_rank)

    # Write initial state (time = 0) + metadata
    if(filter_params.verbose && my_rank == filter_params.master_rank)
        @timeit_debug timer "IO" begin
            unpack_statistics!(filter_data.avg_arr, filter_data.var_arr, filter_data.statistics)
            write_snapshot(filter_params.output_filename,
                           model_data,
                           filter_data.avg_arr,
                           filter_data.var_arr,
                           filter_data.weights,
                           0)
        end
    end

    for it in 1:filter_params.n_time_step

        # integrate true synthetic wavefield
        @timeit_debug timer "True State Update and Process Noise" truth_observations = update_truth!(model_data, nprt_per_rank)

        # Forecast: Update tsunami forecast and get observations from it
        # Parallelised with threads.

        @timeit_debug timer "Particle Dynamics" update_particle_dynamics!(model_data, nprt_per_rank);
        @timeit_debug timer "Particle Proposal" update_particle_proposal!(model_data, filter_data, truth_observations, nprt_per_rank, filter_type)
        @timeit_debug timer "Particle Observations" model_observations = get_particle_observations!(model_data, nprt_per_rank)

        @timeit_debug timer "Particle Weights" get_log_weights!(@view(filter_data.weights[1:nprt_per_rank]),
                                                       truth_observations,
                                                       model_observations,
                                                       filter_params.weight_std)

        # Gather weights to master rank and resample particles.
        # Doing MPI collectives in place to save memory allocations.
        # This style with if statmeents is recommended instead of MPI.Gather_in_place! which is a bit strange.
        # Note that only master_rank allocates memory for all particles. Other ranks only allocate
        # for their chunk of state.
        if my_rank == filter_params.master_rank
            @timeit_debug timer "MPI Gather" MPI.Gather!(MPI.IN_PLACE,
                                                         UBuffer(filter_data.weights, nprt_per_rank),
                                                         filter_params.master_rank,
                                                         MPI.COMM_WORLD)
            @timeit_debug timer "Weights" normalized_exp!(filter_data.weights)
            @timeit_debug timer "Resample" resample!(filter_data.resampling_indices, filter_data.weights, _rng[threadid()])

        else
            @timeit_debug timer "MPI Gather" MPI.Gather!(filter_data.weights,
                                                         nothing,
                                                         filter_params.master_rank,
                                                         MPI.COMM_WORLD)
        end

        # Broadcast resampled particle indices to all ranks
        MPI.Bcast!(filter_data.resampling_indices, filter_params.master_rank, MPI.COMM_WORLD)

        @timeit_debug timer "get_particles" particles = get_particles(model_data)
        @timeit_debug timer "State Copy" copy_states!(particles,
                                                      filter_data.copy_buffer,
                                                      filter_data.resampling_indices,
                                                      my_rank,
                                                      nprt_per_rank)

        @timeit_debug timer "Mean and Var" get_mean_and_var!(filter_data.statistics, particles, filter_params.master_rank)

        if my_rank == filter_params.master_rank && filter_params.verbose

            @timeit_debug timer "IO" begin
                unpack_statistics!(filter_data.avg_arr, filter_data.var_arr, filter_data.statistics)
                write_snapshot(filter_params.output_filename, model_data, filter_data.avg_arr, filter_data.var_arr, filter_data.weights, it)
            end

        end

    end

    if filter_params.enable_timers

        if my_rank == filter_params.master_rank
            print_timer(timer)
        end

        if filter_params.verbose
            # Gather string representations of timers from all ranks and write them on master
            str_timer = string(timer)

            timer_lengths = MPI.Gather(sizeof(str_timer), filter_params.master_rank, MPI.COMM_WORLD)

            if my_rank == filter_params.master_rank
                timer_chars = MPI.Gatherv!(str_timer,
                                           MPI.VBuffer(Vector{UInt8}(undef, sum(timer_lengths)), timer_lengths),
                                           filter_params.master_rank,
                                           MPI.COMM_WORLD)
                @timeit_debug timer "IO" write_timers(timer_lengths, my_size, timer_chars, filter_params)
            else
                MPI.Gatherv!(str_timer, nothing, filter_params.master_rank, MPI.COMM_WORLD)
            end
        end
    end

    unpack_statistics!(filter_data.avg_arr, filter_data.var_arr, filter_data.statistics)

    return get_truth(model_data), filter_data.avg_arr, filter_data.var_arr
end

# Initialise params struct with user-defined dict of values.
function get_params(T, user_input_dict::Dict)

    user_input = (; (Symbol(k) => v for (k,v) in user_input_dict)...)
    params = T(;user_input...)

end

get_params(user_input_dict::Dict) = get_params(FilterParameters, user_input_dict)

# Initialise params struct with default values
get_params() = FilterParameters()

function read_input_file(path_to_input_file::String)

    # Read input provided in a yaml file. Overwrite default input parameters with the values provided.
    if isfile(path_to_input_file)
        user_input_dict = YAML.load_file(path_to_input_file)
    else
        @warn "Input file " * path_to_input_file * " not found, using default parameters"
        user_input_dict = Dict()
    end
    return user_input_dict

end

"""
    run_particle_filter(init, path_to_input_file::String, filter_type::ParticleFilter)

Run the particle filter.  `init` is the function which initialise the model,
`path_to_input_file` is the path to the YAML file with the input parameters.
`filter_type` is the particle filter to use.  See [`ParticleFilter`](@ref) for
the possible values.
"""
function run_particle_filter(init, path_to_input_file::String, filter_type::ParticleFilter; rng::Union{Random.AbstractRNG,Nothing}=nothing)

    MPI.Init()

    # Do I/O on rank 0 only and then broadcast params
    if MPI.Comm_rank(MPI.COMM_WORLD) == 0

        user_input_dict = read_input_file(path_to_input_file)

    else

        user_input_dict = nothing

    end

    user_input_dict = MPI.bcast(user_input_dict, 0, MPI.COMM_WORLD)

    return run_particle_filter(init, user_input_dict, filter_type; rng)

end

"""
    run_particle_filter(init, user_input_dict::Dict, filter_type::ParticleFilter)

Run the particle filter.  `init` is the function which initialise the model,
`user_input_dict` is the list of input parameters, as a `Dict`.  `filter_type`
is the particle filter to use.  See [`ParticleFilter`](@ref) for the possible
values.
"""
function run_particle_filter(init, user_input_dict::Dict, filter_type::ParticleFilter; rng::Union{Random.AbstractRNG,Nothing}=nothing)

    filter_params = get_params(FilterParameters, get(user_input_dict, "filter", Dict()))
    model_params_dict = get(user_input_dict, "model", Dict())

    return run_particle_filter(init, filter_params, model_params_dict, filter_type; rng)

end

end # module


text:
function GlobalIonosphereMapModel(arg0::JString)
    return GlobalIonosphereMapModel((JString,), arg0)
end

function GlobalIonosphereMapModel(arg0::JString, arg1::DataProvidersManager, arg2::TimeScale)
    return GlobalIonosphereMapModel((JString, DataProvidersManager, TimeScale), arg0, arg1, arg2)
end

function get_parameters_drivers(obj::GlobalIonosphereMapModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function get_tec(obj::GlobalIonosphereMapModel, arg0::AbsoluteDate, arg1::GeodeticPoint)
    return jcall(obj, "getTEC", jdouble, (AbsoluteDate, GeodeticPoint), arg0, arg1)
end

function get_tec(obj::GlobalIonosphereMapModel, arg0::FieldAbsoluteDate, arg1::GeodeticPoint)
    return jcall(obj, "getTEC", RealFieldElement, (FieldAbsoluteDate, GeodeticPoint), arg0, arg1)
end

function path_delay(obj::GlobalIonosphereMapModel, arg0::AbsoluteDate, arg1::GeodeticPoint, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "pathDelay", jdouble, (AbsoluteDate, GeodeticPoint, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function path_delay(obj::GlobalIonosphereMapModel, arg0::FieldAbsoluteDate, arg1::GeodeticPoint, arg2::RealFieldElement, arg3::jdouble)
    return jcall(obj, "pathDelay", RealFieldElement, (FieldAbsoluteDate, GeodeticPoint, RealFieldElement, jdouble), arg0, arg1, arg2, arg3)
end

function path_delay(obj::GlobalIonosphereMapModel, arg0::FieldSpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{RealFieldElement})
    return jcall(obj, "pathDelay", RealFieldElement, (FieldSpacecraftState, TopocentricFrame, jdouble, Vector{RealFieldElement}), arg0, arg1, arg2, arg3)
end

function path_delay(obj::GlobalIonosphereMapModel, arg0::SpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{jdouble})
    return jcall(obj, "pathDelay", jdouble, (SpacecraftState, TopocentricFrame, jdouble, Vector{jdouble}), arg0, arg1, arg2, arg3)
end


text:
using PyCall, Test

function test_venv_has_python(path)
    newpython = PyCall.python_cmd(venv=path).exec[1]
    if !isfile(newpython)
        @info """
        Python executable $newpython does not exists.
        This directory contains only the following files:
        $(join(readdir(dirname(newpython)), '\n'))
        """
    end
    @test isfile(newpython)
end

function test_venv_activation(path)
    newpython = PyCall.python_cmd(venv=path).exec[1]

    # Run a fresh Julia process with new Python environment
    code = """
    $(Base.load_path_setup_code())
    using PyCall
    println(PyCall.pyimport("sys").executable)
    println(PyCall.pyimport("sys").exec_prefix)
    println(PyCall.pyimport("pip").__file__)
    """
    # Note that `pip` is just some arbitrary non-standard
    # library.  Using standard library like `os` does not work
    # because those files are not created.
    env = copy(ENV)
    env["PYCALL_JL_RUNTIME_PYTHON"] = newpython
    jlcmd = setenv(`$(Base.julia_cmd()) --startup-file=no -e $code`, env)
    if Sys.iswindows()
        # Marking the test broken in Windows.  It seems that
        # venv copies .dll on Windows and libpython check in
        # PyCall.__init__ detects that.
        @test_broken begin
            output = read(jlcmd, String)
            sys_executable, exec_prefix, mod_file = split(output, "\n")
            newpython == sys_executable
        end
    else
        output = read(jlcmd, String)
        sys_executable, exec_prefix, mod_file = split(output, "\n")
        @test newpython == sys_executable
        @test startswith(exec_prefix, path)
        @test startswith(mod_file, path)
    end
end

@testset "virtualenv activation" begin
    pyname = "python$(pyversion.major).$(pyversion.minor)"
    if Sys.which("virtualenv") === nothing
        @info "No virtualenv command. Skipping the test..."
    elseif Sys.which(pyname) === nothing
        @info "No $pyname command. Skipping the test..."
    else
        mktempdir() do tmppath
            if PyCall.pyversion.major == 2
                path = joinpath(tmppath, "kind")
            else
                path = joinpath(tmppath, "ϵνιℓ")
            end
            run(`virtualenv --python=$pyname $path`)
            test_venv_has_python(path)

            newpython = PyCall.python_cmd(venv=path).exec[1]
            venv_libpython = PyCall.find_libpython(newpython)
            if venv_libpython != PyCall.libpython
                @info """
                virtualenv created an environment with incompatible libpython:
                    $venv_libpython
                """
                return
            end

            test_venv_activation(path)
        end
    end
end

@testset "venv activation" begin
    # In case PyCall is built with a Python executable created by
    # `virtualenv`, let's try to find the original Python executable.
    # Otherwise, `venv` does not work with this Python executable:
    # https://bugs.python.org/issue30811
    sys = PyCall.pyimport("sys")
    if hasproperty(sys, :real_prefix)
        # sys.real_prefix is set by virtualenv and does not exist in
        # standard Python:
        # https://github.com/pypa/virtualenv/blob/16.0.0/virtualenv_embedded/site.py#L554
        candidates = [
            PyCall.venv_python(sys.real_prefix, "$(pyversion.major).$(pyversion.minor)"),
            PyCall.venv_python(sys.real_prefix, "$(pyversion.major)"),
            PyCall.venv_python(sys.real_prefix),
            PyCall.pyprogramname,  # must exists
        ]
        python = candidates[findfirst(isfile, candidates)]
    else
        python = PyCall.pyprogramname
    end

    if PyCall.conda
        @info "Skip venv test with conda."
    elseif !success(PyCall.python_cmd(`-c "import venv"`, python=python))
        @info "Skip venv test since venv package is missing."
    else
        mktempdir() do tmppath
            if PyCall.pyversion.major == 2
                path = joinpath(tmppath, "kind")
            else
                path = joinpath(tmppath, "ϵνιℓ")
            end
            # Create a new virtual environment
            run(PyCall.python_cmd(`-m venv $path`, python=python))
            test_venv_has_python(path)
            test_venv_activation(path)
        end
    end
end


text:
function GetBatch(minibatchstart,BatchSize,Ntrain)

    if minibatchstart>Ntrain-BatchSize+1
        minibatchstart=Ntrain-BatchSize+1
    end
    minibatch=minibatchstart:minibatchstart+BatchSize-1
    minibatchstart=minibatchstart+BatchSize
    if minibatchstart>Ntrain
        minibatchstart=1
    end
    return minibatchstart,minibatch
end
export GetBatch

function GetNewBatch(OldMiniBatch,Ntrain)
    MiniBatch=rem(OldMiniBatch[end]+1:OldMiniBatch[end]+length(OldMiniBatch),Ntrain)
    MiniBatch[MiniBatch.==0]=Ntrain
    return MiniBatch
end
export GetNewBatch

# Standard Gradient Descent:

function GradientDescentUpdate!(x,grad,LearningRate)
    # xnew=x-LearningRate*grad
    axpy!(-LearningRate,grad,x)
end
export GradientDescentUpdate!

# Gradient Descent with Momentum:

function GradientDescentMomentumUpdate!(x,grad,avgrad,Momentum,LearningRate)
    # avgrad_new=alpha*avgrad_old +(1-alpha)*avgrad
    # xnew=x-LearningRate*avgrad
    # where alpha=Momentum
    scale!(Momentum,avgrad)
    axpy!(1-Momentum,grad,avgrad)
    axpy!(-LearningRate,avgrad,x)
end
export GradientDescentMomentumUpdate!

function GradientDescentMomentumInit(net)
    avgrad=Array(Any,length(net.node))
    for par in Parameters(net)
        avgrad[par]=cArray(net.gpu,zeros(size(net.value[par]))) # initial average gradient
    end
    return avgrad
end
export GradientDescentMomentumInit

# Nesterov Accelerated Gradient Descent:

function NesterovInit(net)
    velo=Array(Any,length(net.node)) # needed for Nesterov
    for par in Parameters(net)
#        velo[par]=cArray(size(net.value[par])); fill!(velo[par],0.0) # Nesterov velocity
        velo[par]=cArray(net.gpu,zeros(net.eltype,size(net.value[par]))) # Nesterov velocity
    end
    return velo
end
export NesterovInit

function NesterovGradientDescentUpdate!(thetaP,gradP,v,LearningRate,t)
    # A standard Nesterov Update (see for example appendix A.1 in mlr.org/proceedings/papers/v28/sutskever13-supp.pdf) uses
    # v_new = mu*v_old-epsilon*gradient(theta_old+mu*v_old)
    # theta_new=theta_old+v_new
    # However, we can reparameterise this using
    # thetaP = theta + v
    # and write an update for thetaP
    # v_new = mu*v_old-epsilon*gradient(thetaP_old)
    # thetaP_new = thetaP_old-mu_old*v_old+v_new*(1+mu_new)
    # Formally, we would then need to retransform back to get the value for theta.
    # However, as we converge, mu goes to 1 and the gradient goes to zero.
    # This means that thetaP tends to theta anyway as we converge to the optimum.

    mu_new=min(0.99,1-3/(t+5))
    mu_old=min(0.99,1-3/(t+4))
    scale!(mu_old,v)
    axpy!(mu_new,v,thetaP)
    axpy!(-LearningRate,gradP,v)
    axpy!(-(1.0+mu_new)*LearningRate,gradP,thetaP)

end
export NesterovGradientDescentUpdate!

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# Autogenerated wrapper script for Goalign_jll for i686-linux-musl
export goalign

JLLWrappers.@generate_wrapper_header("Goalign")
JLLWrappers.@declare_executable_product(goalign)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        goalign,
        "bin/goalign",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module BioMakie

export  bondshape,
        resletterdict,
        viewmsa,
		atomcolors,
        atomcoords,
        atomradii,
		viewstruc,
        resatoms,
        resbonds,
        backbonebonds,
        defaultatom,
        defaultresidue,
        covrad,
        vdwrad

using BioStructures, MolecularGraph, MIToS
using Lazy, SplitApplyCombine, TensorCast
using DataStructures, DelimitedFiles, JLD2, JSON3, HTTP
using Distances, Distributions, GeometryBasics, Colors
using FileIO, FastaIO, OrderedCollections
using GraphMakie, JSServe, Meshes
using Makie, GLMakie, WGLMakie
GLMakie.activate!()

include("../src/utils.jl")
include("../src/chemdata.jl")
include("../src/bonds.jl")
include("../src/structure.jl")
include("../src/msa.jl")

end # BioMakie


text:
# Comment here
const bar = '\n'

"""
    test_fun(a::Int)
For test only
"""
function test_fun(a::Int, b::T) where T <: Number
    println(a)
    println("foo $(bar)")
end

@enum Unicode α=1 β=2

res = [√i for i in 1:10]

#= Dummy function =#
test_fun²(:sym, true, raw"test", `echo 1`)

# function to calculate the volume of a sphere
function sphere_vol(r)
    # julia allows Unicode names (in UTF-8 encoding)
    # so either "pi" or the symbol π can be used
    return 4/3*pi*r^3
end

# functions can also be defined more succinctly
quadratic(a, sqr_term, b) = (-b + sqr_term) / 2a

# calculates x for 0 = a*x^2+b*x+c, arguments types can be defined in function definitions
function quadratic2(a::Float64, b::Float64, c::Float64)
    # unlike other languages 2a is equivalent to 2*a
    # a^2 is used instead of a**2 or pow(a,2)
    sqr_term = sqrt(b^2-4a*c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)
    # multiple values can be returned from a function using tuples
    # if the return keyword is omitted, the last term is returned
    r1, r2
end

vol = sphere_vol(3)
# @printf allows number formatting but does not automatically append the \n to statements, see below
using Printf
@printf "volume = %0.3f\n" vol
#> volume = 113.097

quad1, quad2 = quadratic2(2.0, -2.0, -12.0)
println("result 1: ", quad1)
#> result 1: 3.0
println("result 2: ", quad2)
#> result 2: -2.0

s1 = "The quick brown fox jumps over the lazy dog α,β,γ"

# search returns the first index of a char
i = findfirst(isequal('b'), s1)
println(i)
#> 11
# the second argument is equivalent to the second argument of split, see below

# or a range if called with another string
r = findfirst("brown", s1)
println(r)
#> 11:15

# string replace is done thus:
r = replace(s1, "brown" => "red")
show(r); println()
#> "The quick red fox jumps over the lazy dog α,β,γ"

# search and replace can also take a regular expressions by preceding the string with 'r'
r = findfirst(r"b[\w]*n", s1)
println(r)
#> 11:15

# again with a regular expression
r = replace(s1, r"b[\w]*n" => "red")
show(r); println()
#> "The quick red fox jumps over the lazy dog α,β,γ"

# there are also functions for regular expressions that return RegexMatch types
# match scans left to right for the first match (specified starting index optional)
r = match(r"b[\w]*n", s1)
println(r)
#> RegexMatch("brown")

# RegexMatch types have a property match that holds the matched string
show(r.match); println()
#> "brown"

# eachmatch returns an iterator over all the matches
r = eachmatch(r"[\w]{4,}", s1)
for i in r print("\"$(i.match)\" ") end
#> "quick" "brown" "jumps" "over" "lazy"
println()

r = collect(m.match for m = eachmatch(r"[\w]{4,}", s1))
println(r)
#> SubString{String}["quick", "brown", "jumps", "over", "lazy"]

# a string can be repeated using the repeat function,
# or more succinctly with the ^ syntax:
r = "hello "^3
show(r); println() #> "hello hello hello "

# the strip function works the same as python:
# e.g., with one argument it strips the outer whitespace
r = strip("hello ")
show(r); println() #> "hello"
# or with a second argument of an array of chars it strips any of them;
r = strip("hello ", ['h', ' '])
show(r); println() #> "ello"
# (note the array is of chars and not strings)

# similarly split works in basically the same way as python:
r = split("hello, there,bob", ',')
show(r); println() #> SubString{String}["hello", " there", "bob"]
r = split("hello, there,bob", ", ")
show(r); println() #> SubString{String}["hello", "there,bob"]
r = split("hello, there,bob", [',', ' '], limit=0, keepempty=false)
show(r); println() #> SubString{String}["hello", "there", "bob"]
# (the last two arguements are limit and include_empty, see docs)

# the opposite of split: join is simply
r = join(collect(1:10), ", ")
println(r) #> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

function printsum(a)
    # summary generates a summary of an object
    println(summary(a), ": ", repr(a))
end

# arrays can be initialised directly:
a1 = [1,2,3]
printsum(a1)
#> 3-element Array{Int64,1}: [1, 2, 3]

# or initialised empty:
a2 = []
printsum(a2)
#> 0-element Array{Any,1}: Any[]

# since this array has no type, functions like push! (see below) don't work
# instead arrays can be initialised with a type:
a3 = Int64[]
printsum(a3)
#> 0-element Array{Int64,1}: Int64[]

# ranges are different from arrays:
a4 = 1:20
printsum(a4)
#> 20-element UnitRange{Int64}: 1:20

# however they can be used to create arrays thus:
a4 = collect(1:20)
printsum(a4)
#> 20-element Array{Int64,1}: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
#>  15, 16, 17, 18, 19, 20]

# arrays can also be generated from comprehensions:
a5 = [2^i for i = 1:10]
printsum(a5)
#> 10-element Array{Int64,1}: [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]

# arrays can be any type, so arrays of arrays can be created:
a6 = (Array{Int64, 1})[]
printsum(a6)
#> 0-element Array{Array{Int64,1},1}: Array{Int64,1}[]
# (note this is a "jagged array" (i.e., an array of arrays), not a multidimensional array,
# these are not covered here)

# Julia provided a number of "Dequeue" functions, the most common
# for appending to the end of arrays is push!
# ! at the end of a function name indicates that the first argument is updated.

push!(a1, 4)
printsum(a1)
#> 4-element Array{Int64,1}: [1, 2, 3, 4]

# push!(a2, 1) would cause error:

push!(a3, 1)
printsum(a3) #> 1-element Array{Int64,1}: [1]
#> 1-element Array{Int64,1}: [1]

push!(a6, [1,2,3])
printsum(a6)
#> 1-element Array{Array{Int64,1},1}: Array{Int64,1}[[1, 2, 3]]

# using repeat() to create arrays
# you must use the keywords "inner" and "outer"
# all arguments must be arrays (not ranges)
a7 = repeat(a1,inner=[2],outer=[1])
printsum(a7)
#> 8-element Array{Int64,1}: [1, 1, 2, 2, 3, 3, 4, 4]
a8 = repeat(collect(4:-1:1),inner=[1],outer=[2])
printsum(a8)
#> 8-element Array{Int64,1}: [4, 3, 2, 1, 4, 3, 2, 1]

# try, catch can be used to deal with errors as with many other languages
try
    push!(a,1)
catch err
    showerror(stdout, err, backtrace());println()
end
#> UndefVarError: a not defined
#> Stacktrace:
#>  [1] top-level scope at C:\JuliaByExample\src\error_handling.jl:5
#>  [2] include at .\boot.jl:317 [inlined]
#>  [3] include_relative(::Module, ::String) at .\loading.jl:1038
#>  [4] include(::Module, ::String) at .\sysimg.jl:29
#>  [5] exec_options(::Base.JLOptions) at .\client.jl:229
#>  [6] _start() at .\client.jl:421
println("Continuing after error")
#> Continuing after error

# repeat can be useful to expand a grid
# as in R's expand.grid() function:

m1 = hcat(repeat([1,2],inner=[1],outer=[3*2]),
          repeat([1,2,3],inner=[2],outer=[2]),
          repeat([1,2,3,4],inner=[3],outer=[1]))
printsum(m1)
#> 12×3 Array{Int64,2}: [1 1 1; 2 1 1; 1 2 1; 2 2 2; 1 3 2; 2 3 2; 1 1 3; 2 1 3;
#>   1 2 3; 2 2 4; 1 3 4; 2 3 4]

# for simple repetitions of arrays,
# use repeat
m2 = repeat(m1,1,2)     # replicate a9 once into dim1 and twice into dim2
println("size: ", size(m2))
#> size: (12, 6)

m3 = repeat(m1,2,1)     # replicate a9 twice into dim1 and once into dim2
println("size: ", size(m3))
#> size: (24, 3)

# Julia comprehensions are another way to easily create
# multidimensional arrays

m4 = [i+j+k for i=1:2, j=1:3, k=1:2] # creates a 2x3x2 array of Int64
m5 = ["Hi Im # $(i+2*(j-1 + 3*(k-1)))" for i=1:2, j=1:3, k=1:2]
# expressions are very flexible
# you can specify the type of the array by just
# placing it in front of the expression

import LegacyStrings
m5 = LegacyStrings.ASCIIString["Hi Im element # $(i+2*(j-1 + 3*(k-1)))" for i=1:2, j=1:3, k=1:2]
printsum(m5)
#> 2×3×2 Array{LegacyStrings.ASCIIString,3}: LegacyStrings.ASCIIString[
#>   "Hi Im element # 1" "Hi Im element # 3" "Hi Im element # 5";
#>   "Hi Im element # 2" "Hi Im element # 4" "Hi Im element # 6"]
#>
#> LegacyStrings.ASCIIString["Hi Im element # 7" "Hi Im element # 9"
#>   "Hi Im element # 11"; "Hi Im element # 8" "Hi Im element # 10" "Hi Im element # 12"]

# Array reductions
# many functions in Julia have an array method
# to be applied to specific dimensions of an array:

sum(m4, dims=3)        # takes the sum over the third dimension
sum(m4, dims=(1,3))  # sum over first and third dim

maximum(m4, dims=2)    # find the max elt along dim 2
findmax(m4, dims=3)    # find the max elt and its index along dim 3
                    # (available only in very recent Julia versions)

# Broadcasting
# when you combine arrays of different sizes in an operation,
# an attempt is made to "spread" or "broadcast" the smaller array
# so that the sizes match up. broadcast operators are preceded by a dot:

m4 .+ 3       # add 3 to all elements
m4 .+ [1,2]      # adds vector [1,2] to all elements along first dim

# slices and views
m4=m4[:,:,1] # holds dim 3 fixed
m4[:,2,:]  # that's a 2x1x2 array. not very intuititive to look at

# get rid of dimensions with size 1:
dropdims(m4[:,2,:], dims=2) # that's better

# assign new values to a certain view
m4[:,:,1] = rand(1:6,2,3)
printsum(m4)
#> 2×3 Array{Int64,2}: [3 5 3; 1 3 5]

# (for more examples of try, catch see Error Handling above)
try
    # this will cause an error, you have to assign the correct type
    m4[:,:,1] = rand(2,3)
catch err
    println(err)
end
#> InexactError(:Int64, Int64, 0.7603891754678744)

try
    # this will cause an error, you have to assign the right shape
    m4[:,:,1] = rand(1:6,3,2)
catch err
    println(err)
end
#> DimensionMismatch("tried to assign 3×2 array to 2×3×1 destination")

# dicts can be initialised directly:
a1 = Dict(1=>"one", 2=>"two")
printsum(a1)
#> Dict{Int64,String} with 2 entries: Dict(2=>"two",1=>"one")

# then added to:
a1[3]="three"
printsum(a1)
#> Dict{Int64,String} with 3 entries: Dict(2=>"two",3=>"three",1=>"one")
# (note dicts cannot be assumed to keep their original order)

# dicts may also be created with the type explicitly set
a2 = Dict{Int64, AbstractString}()
a2[0]="zero"
printsum(a2)
#> Dict{Int64,AbstractString} with 1 entry: Dict{Int64,AbstractString}(0=>"zero")

# dicts, like arrays, may also be created from comprehensions
using Printf
a3 = Dict([i => @sprintf("%d", i) for i = 1:10])
printsum(a3)
#> Dict{Int64,String} with 10 entries: Dict(7=>"7",4=>"4",9=>"9",10=>"10",
#>  2=>"2",3=>"3",5=>"5",8=>"8",6=>"6",1=>"1")

# as you would expect, Julia comes with all the normal helper functions
# for dicts, e.g., haskey
println(haskey(a1,1)) #> true

# which is equivalent to
println(1 in keys(a1)) #> true
# where keys creates an iterator over the keys of the dictionary

# similar to keys, values get iterators over the dict's values:
printsum(values(a1))
#> Base.ValueIterator for a Dict{Int64,String} with 3 entries: ["two", "three", "one"]

# use collect to get an array:
printsum(collect(values(a1)))
#> 3-element Array{String,1}: ["two", "three", "one"]

if true
    println("It's true!")
else
    println("It's false!")
end
#> It's true!

if false
   println("It's true!")
else
   println("It's false!")
end
#> It's false!

# Numbers can be compared with opperators like <, >, ==, !=

1 == 1.
#> true

1 > 2
#> false

"foo" != "bar"
#> true

# and many functions return boolean values

occursin("that", "this and that")
#> true

# More complex logical statments can be achieved with `elseif`

function checktype(x)
   if x isa Int
      println("Look! An Int!")
   elseif x isa AbstractFloat
      println("Look! A Float!")
   elseif x isa Complex
      println("Whoa, that's complex!")
   else
      println("I have no idea what that is")
   end
end

checktype(2)
#> Look! An Int!

checktype(√2)
#> Look! A Float!

checktype(√Complex(-2))
#> Whoa, that's complex!

checktype("who am I?")
#> I have no idea what that is

# For simple logical statements, one can be more terse using the "ternary operator",
# which takes the form `predicate ? do_if_true : do_if_false`

1 > 2 ? println("that's true!") : println("that's false!")
#> that's false!

noisy_sqrt(x) = x ≥ 0 ? sqrt(x) : "That's negative!"
noisy_sqrt(4)
#> 2.0
noisy_sqrt(-4)
#> That's negative!

# "Short-circuit evaluation" offers another option for conditional statements.
# The opperators `&&` for AND and `||` for OR only evaluate the right-hand
# statement if necessary based on the predicate.
# Logically, if I want to know if `42 == 0 AND x < y`,
# it doesn't matter what `x` and `y` are, since the first statement is false.
# This can be exploited to only evaluate a statement if something is true -
# the second statement doesn't even have to be boolean!

everything = 42
everything < 100 && println("that's true!")
#> "that's true!"
everything == 0 && println("that's true!")
#> false

√everything > 0 || println("that's false!")
#> true
√everything == everything || println("that's false!")
#> that's false!


text:
include("train.jl")

stack = @stack MetaTac stack_input=true begin
  Conv(1024, relu, stride = 3, padding = 0)
  Dense(512, relu)
end

model = NeuralModel(MetaTac, stack) |> to_gpu |> Async
player = MCTSPlayer(model, power = 250)

mcts_opponents = [MCTSPlayer(power = p) for p in [10, 50, 100, 200, 500]]
self_opponents = [MCTSPlayer(model, power = p, name = "current$p") for p in [10, 50]]

loss = Loss(policy = 0.5, reg = 1e-5)

with_contest( train_self!
            , player
            , loss = loss
            , playings = 50
            , batchsize = 100
            , branching = 0.02
            , iterations = 10
            , epochs = 100
            , augment = true
            , opponents = [mcts_opponents; self_opponents]
            , length = 1000
            , interval = 5
            )

max_stars_repo_name: roSievers/jtac


text:
module TestDoctest

import ThreadsX
using Documenter: doctest
using Test

@testset "doctest" begin
    doctest(ThreadsX)
end

end  # module


text:
#=

    pool_data.jl

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

    Original Source: https://github.com/ap-browning/Spheroids/blob/main/Module/Inference/pool_data.jl

=#

"""
    pooled_data(data::DataFrame, groups, vars)

Pool the data in variables `vars` (i.e., subtract the group mean) by `groups`.
"""
function pooled_data(data::DataFrame, groups::Array{Symbol,1}, vars::Array{Symbol,1})

    group_vals  = [unique(data[:,group]) for group ∈ groups]
    filter_vals = collect(Base.product(group_vals...))[:]

    pooled_data = Array{Array{Float64,2},1}(undef,length(filter_vals))

    # Loop through groups
    for i = 1:length(filter_vals)

        ingroup  = .*([data[:,groups[j]] .== filter_vals[i][j] for j = 1:length(groups)]...)
        data_i   = Array(data[ingroup,vars])
        pooled_data[i] = data_i .- mean(data_i,dims=1)
    end

    vcat(pooled_data...)

end
function pooled_data(data::DataFrame, groups, vars)
    res = pooled_data(data,
        typeof(groups) <: Array ? groups : [groups],
        typeof(vars)   <: Array ? vars   : [vars]
    )
    if typeof(vars) == Symbol
        return res[1]
    else
        return res
    end
end

pooled_std(data::DataFrame, groups, vars) = std(pooled_data(data,groups,vars),dims=1)[:]
pooled_cov(data::DataFrame, groups, vars) = cov(pooled_data(data,groups,vars))
pooled_cor(data::DataFrame, groups, vars) = cor(pooled_data(data,groups,vars))

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
using DataFrames
using StanSample, Test

ProjDir = @__DIR__
cd(ProjDir)

gq = "
  data {
    int<lower=0> N;
    int<lower=0> y[N];
  }
  parameters {
    real<lower=0.00001> Theta;
  }
  model {
    Theta ~ gamma(4, 1000);
    for (n in 1:N)
      y[n] ~ exponential(Theta);
  }
  generated quantities{
    real y_pred;
    y_pred = exponential_rng(Theta);
  }
";

data = Dict(
  "N" => 3,
  "y" => [100, 950, 450]
);

sm = SampleModel("Generate_quantities", gq, joinpath(@__DIR__, "tmp"));

rc = stan_sample(sm; data)

if success(rc)
  df = stan_generate_quantities(sm, 1, "1")
  df |> display
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
using GAIO
using Interpolations
using QuadGK
using LinearAlgebra
using StaticArrays

# a vector field exhibiting a knotted trajectory
function knotted_v(x, tol)
    A = 2
    B = 1
    L = 0.9
    h = 2
    r = 2
    p = 3
    q = 4
    T = 2*π*A

    γ = t -> (
        cos(t) - (L*A)/(A+B) * cos((A+B)/A * t),
        sin(t) - (L*A)/(A+B) * sin((A+B)/A * t),
       -sin(B/A * t) + h * (1-2*t/T)^p
    )

    γ_t = t -> (
       -sin(t) + (L*A)/A * sin((A+B)/A * t),
        cos(t) - (L*A)/A * cos((A+B)/A * t),
       -cos(B/A * t) * B/A - h * p * (1-2*t/T)^(p-1) * 2/T
    )

    d = (x, y) -> inv(sum(abs.(x .- y) .^ q))

    σ = x -> (x[1]/(1-x[3]), x[2]/(1-x[3]))

    σ_inv = x -> (
        2*x[1]/(1 + x[1]^2 + x[2]^2),
        2*x[2]/(1 + x[1]^2 + x[2]^2),
         1 - 2/(1 + x[1]^2 + x[2]^2)
    )

    # approximate the integral with adaptive Gauss quadrature for given tolerance
    int, = quadgk(t -> begin
        γt = γ(t)
        γdt = γ_t(t)
        nγdt = norm(γdt)
        dγt = d(x, γt)
        σγ = σ(γdt ./ nγdt)
        val = nγdt .* dγt .* σγ
        return SVector(val[1], val[2], nγdt * dγt)
    end, 0, T, rtol=tol)

    nom = (int[1], int[2])
    denom1 = int[3]
    y = (r * x[1] / norm((x[1], x[2])), r * x[2] / norm((x[1], x[2])), x[3])
    denom2 = d(x, y)

    w = nom ./ (denom1 + denom2)

    return σ_inv(w)
end
# interpolate the function v on the mesh [a,b]^3 with n equidistant nodes
# using linear B-Splines
function interpolate_v(v, a, b, n, tol)
    val = Array{SVector{3,Float64},3}(undef, n, n, n)
    x = range(a, b, length = n)
    y = range(a, b, length = n)
    z = range(a, b, length = n)

    @Threads.threads for i in 1:n
        for j in 1:n
            for k in 1:n
                val[i,j,k] = SVector(v((x[i], y[j], z[k]), tol))
            end
        end
    end

    itp = interpolate(val, BSpline(Linear()))
    sitp = scale(itp, x, y, z)

    return t -> sitp(t[1], t[2], t[3])
end

vf = interpolate_v(knotted_v, -3.5, 3.5, 30, 1e-12)
f(x) = rk4_flow_map(vf, x, step_size = 0.075)

center, radius = (0,0,0), (2,2,2)
P = BoxPartition(Box(center, radius))
F = BoxMap(f, P, no_of_points = 200)

C = chain_recurrent_set(F, P[:], steps = 18)

plot(C)


text:
function count_fixed_points(array)
    result = 0
    for (index, value) in enumerate(array)
        if index == value
            result += 1
        end
    end
    return result
end

function count_cycles(array)
    result = 0
    visited = trues(size(array, 1))
    for (index, value) in enumerate(array)
        if visited[index]
            j = index
            while true
                visited[j] = false
                j = array[j]
                if j == index
                    break
                end
            end
            result += 1
        end
    end
    return result
end

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
using Colors

export intersperse,
       class,
       memoize

# Dict macro
macro d(xs...)
  if VERSION < v"0.4-"
    Expr(:dict, map(esc, xs)...)
  else
    :(Dict($(map(esc, xs)...)))
  end
end

convert(::Type{Color}, s::AbstractString) =
    parse(Colorant, s)

render_color(c) = string("#" * hex(c))

# Utility functions for Elem

boolattr(a) = a ? true : nothing # I know, this is insane

@doc """
Intersperse a value in between elements in a vector

Optionally you can tell it to enclose the result in the seperator element.

e.g.
```
intersperse(0, [1, 2, 3])
# => [1, 0, 2, 0, 3]

intersperse(0, [1, 2, 3], true)
# => [0, 1, 0, 2, 0, 3, 0]
```
""" ->
intersperse(x, xs, enclose=false) = begin
    if length(xs) > 1
        res = foldl((acc, nxt) -> vcat(acc, x, nxt),
                    Any[xs[1]], xs[2:end])
    else
        res = xs
    end
    enclose ? [x, res, x] : res
end

getproperty(el::Elem, prop, default) =
    hasproperties(el) ? get(properties(el), prop, default) : default

wrapmany(t::TileList, wrap, state) =
    length(t.tiles) == 1 ?
        render(t.tiles[1], state) :
        Elem(wrap, map(t -> render(t, state), t.tiles))

@api class => (Class <: Tile) begin
    doc("Add a HTML class.")
    arg(class::AbstractString, doc="Space separated classes.")
    curry(content::TileList, doc="A tile or a vector of tiles.")
    kwarg(
        forcewrap::Bool=false,
        doc="""If set to true, contents will be put in a containing tag and the
               classes are set on the container, even if there is only one tile.
               """
    )
    kwarg(wrap::Symbol=:div, doc="The tag to use for the container.")
end

maybestring(s::AbstractString, state) = s
maybestring(s::TileList, state) =
    length(s.tiles) == 1 ? maybestring(s.tiles[1], state) : render(s.tiles, state)
maybestring(s, state) = render(s, state)

addclasses(t, cs) =
    t & @d(:className => cs * " " * getproperty(t, :className, ""))

render(c::Class, state) =
    addclasses(c.forcewrap ? Elem(c.wrap, maybestring(c.content, state)) :
                           wrapmany(c.content, c.wrap, state),
               c.class)

@doc """
given a sentinal, vector of parts, prefix, suffix and a value,

if the vector of parts is referentially equal to the sentinal, then returns
    [prefix * suffix => value]
otherwise returns
    [prefix * name(part) * suffix => value for part in parts]

Can be used while rendering border/padding for different sides, or border
radius for different corners etc.
""" ->
mapparts(sentinal, parts, prefix, suffix, value) =
   parts === sentinal ?
       @d(prefix * suffix => value) :
       Dict(Pair[prefix * name(part) * suffix => value for part in parts])

style(x) = @d(:style => x)

teeprint(x, fn=println) = begin
    fn(x)
    x
end

# @doc """
# memoize anything's rendered output and state
# """ ->
#
# @api memoize => (Memoized<:Tile) begin
#     arg(tile::Any)
#     arg(store::WeakKeyDict=WeakKeyDict())
# end

deepmerge!(a::Associative, b::Associative) = begin
    for (k, v) in b
        @show k, v
        if isa(v, Associative) && haskey(a, k) && isa(a[k], Associative)
            @show k, "AAA"
            a[k] = deepmerge!(a[k], v)
        else
            a[k] = b[k]
        end
    end
    a
end

# Todo: allow specifying hash function,
# allow mirroring to JLD optionally
# render(m::Memoized, state) = begin
#     if haskey(m.store, m.tile)
#         deepmerge!(state, st) # This.
#         m.store[m.tile]
#     else
#         st = Dict("embedded_signals"=>Dict())
#         elem = render(convert(Tile, m.tile), st)
#         deepmerge!(state, st)
#       m.store[m.tile] = (elem, st)
#       elem
#   end
#end


text:
"""
# DuplexIOs.jl

Combine two one-direction `IO` objects into a single bi-direction `IO` object.
"""
module DuplexIOs

export DuplexIO

struct DuplexIO{IN <: IO, OUT <: IO} <:IO
    in::IN
    out::OUT
end

# Duplex methods.

Base.isopen(d::DuplexIO) = (Base.isopen(d.in) && Base.isopen(d.out))
Base.close(d::DuplexIO) = (Base.close(d.in); Base.close(d.out))
Base.lock(d::DuplexIO) = (Base.lock(d.in); Base.lock(d.out))
Base.unlock(d::DuplexIO) = (Base.unlock(d.in); Base.unlock(d.out))

Base.skipchars(p, d::DuplexIO; k...) = Base.skipchars(p, d.in; k...)

function Base.stat(d::DuplexIO)
    s1 = stat(d.in)
    s2 = stat(d.out)
    s1 == s2 || throw(ArgumentError("""
        Can't `stat` a DuplexIO unless `.in` and `.out` refer to the same file.
        Try `stat(io.in)` or `stat(io.out)`.
        """))
    return s1
end

@static if isdefined(Base, :shutdown)
    Base.shutdown(d::DuplexIO) = close(d.out)
end

"""
Types of arg 2 for methods of f.
"""
arg2_types(f) = (m.nargs < 3 || m.sig isa UnionAll ? missing :
                 m.sig.types[3] == Vararg          ? missing :
                                                     m.sig.types[3]
                 for m in methods(f)
                ) |> skipmissing |> unique

"""
    @wrap [x] f

Wrap method `f` for DuplexIO field `x` (`in` or `out`):

    Base.f(d::DuplexIO, a...) = Base.f(d.x, a...)

Or if `x` is omitted, wrap for both fields:

    Base.f(d::DuplexIO, a...) = (Base.f(d.in, a...); Base.f(d.out, a...))
"""
macro wrap(a, b=nothing)

    x, f = (b == nothing) ? (nothing, a) : (a, b)

    if x != nothing
        Expr(:block,

             esc(:(($f(d::DuplexIO; k...) = $f(d.$x; k...)))),

            (esc(:(($f(d::DuplexIO, a::$T, aa...; k...) =
                    $f(d.$x, a, aa...; k...))))
             for T in arg2_types(Main.eval(:($f))))...
        )
    else
         esc(:(($f(d::DuplexIO, a...; k...) = ($f(d.in, a...; k...),
                                              ($f(d.out, a...; k...))))))
    end
end

# Output methods.

@wrap out Base.iswritable
@wrap out Base.buffer_writes
@wrap out Base.flush
@wrap out Base.write
@wrap out Base.unsafe_write

# Input methods.

@wrap in Base.isreadable
@wrap in Base.eof
@wrap in Base.bytesavailable
@wrap in Base.read
@wrap in Base.read!
#@wrap in Base.readbytes! FIXME "WARNING: Method definition readbytes...." ???
@wrap in Base.readuntil
@wrap in Base.readline
@wrap in Base.countlines
@wrap in Base.eachline
@static if isdefined(Base, :readeach)
    @wrap in Base.readeach
end
@wrap in Base.unsafe_read
@wrap in Base.peek
@wrap in Base.readavailable
@wrap in Base.mark
@wrap in Base.ismarked
@wrap in Base.unmark
@wrap in Base.reset

end #module DuplexIOs


text:
struct GappedTopKBatchQs <: ExtendingBatchQs
    sequentialQs::SequentialPQs
    k::Int
    m::Int

    function GappedTopKBatchQs(sequentialQs::SequentialPQs; k::Int, m::Int=0)::RandomBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))
        (m < 0) && throw(ArgumentError("Invalid batch increase m=$(m)."))
        return new(sequentialQs, k, m)
    end
end

"""
Use sequential strategy, randomly select batch_size observations from batch_size + m
observations with highest usefulness.
"""
function select_batch(qs::GappedTopKBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    sequential_strategy = qs.sequentialQs
    candidate_scores = qs_score(sequential_strategy, x, labels)[candidate_indices]
    descending_indices = sortperm(candidate_scores; rev=true)
    # find best k+m elements
    best_k_m = candidate_indices[descending_indices[1:min(qs.k + qs.m, end)]]
    # take random sample of size k from best elements
    best_sampled = sample(best_k_m, qs.k, replace=false)
    return best_sampled
end


text:
__precompile__()

"""
    Phylo package

The `Phylo` package provides some simple phylogenetics types (e.g. NamedTree)
to interface to the `Diversity` package for measuring phylogenetic diversity.
It also provides an interface to `R` for copying trees to and from that
language and can read newick and nexus tree files (including `TreeSet`s that
contain multiple trees).

Finally it also provides a standard abstract interface to
phylogenetic trees, by defining `AbstractNode`, `AbstractBranch` and
`AbstractTree` supertypes, and methods to interface to them. It also
provides (through the `Phylo.API` submodule) methods to
(re)define to write your own phylogenetic type in a way that will
interact cleanly with other phylogenetic packages.
"""
module Phylo

import Base: Pair, Tuple, show, eltype, length, getindex
import Compat: IteratorSize, IteratorEltype

abstract type AbstractNode end
abstract type AbstractBranch end
abstract type AbstractTree{NodeLabel, BranchLabel} end
abstract type AbstractBranchTree{NL, BL} <: AbstractTree{NL, BL} end
export AbstractNode, AbstractBranch, AbstractTree

"""
    Phylo.API submodule

The `Phylo.API` submodule should be `import`ed if you want to
create a new phylogeny, node or branch subtype. Otherwise it can be
ignored.
"""
module API
include("API.jl")
# AbstractTree methods
export _ntrees, _addbranch!, _deletebranch!, _branch!, _setbranch!
export _addnode!, _addnodes!, _deletenode!, _setnode!
export _getnodenames, _hasnode, _getnode, _getnodes
export _getbranchnames, _hasbranch, _getbranch, _getbranches
export _hasrootheight, _getrootheight, _setrootheight!, _clearrootheight!
export _nodetype, _branchtype
export _extractnode, _extractbranch
export _extractnodename, _extractbranchname
export _getleafinfo, _setleafinfo!, _leafinfotype
export _getnoderecord, _setnoderecord!
export _hasheight, _getheight, _setheight!
export _hasparent, _getparent, _getancestors
export _haschildren, _getchildren, _getdescendants
export _validate
export _getleafnames, _resetleaves!, _nleaves

# AbstractNode methods
export _isleaf, _isroot, _isinternal, _isunattached
export _indegree, _hasinbound, _getinbound, _setinbound!, _deleteoutbound!
export _outdegree, _getoutbounds, _addoutbound!, _deleteoutbound!
export _hasoutboundspace, _hasinboundspace

# AbstractBranch methods
export _src, _dst, _getlength
export _setsrc!, _setdst!

# Label names
export _newnodelabel, _newbranchlabel

end

include("Interface.jl")
# AbstractTree methods
export ntrees, nodetype, branchtype, nodenametype, branchnametype
export addbranch!, deletebranch!, branch!
export addnode!, addnodes!, deletenode!
export getnodenames, hasnode, getnode, getnodes
export getbranchnames, hasbranch, getbranch, getbranches
export hasrootheight, getrootheight, setrootheight!
export hasparent, getparent, getancestors
export haschildren, getchildren, getdescendants
export validate

# AbstractTree / AbstractNode methods
export isleaf, isroot, isinternal, isunattached
export indegree, outdegree, hasinbound, getinbound, getoutbounds
export hasoutboundspace, hasinboundspace
export getleafnames, resetleaves, nleaves
export getleafinfo, setleafinfo!, leafinfotype
export getnoderecord, setnoderecord!
export hasheight, getheight, setheight!

# AbstractTree / AbstractBranch methods
export src, dst, getlength
export changesrc!, changedst!

include("Info.jl")
export LeafInfo

include("Branch.jl")
export Branch

include("Node.jl")
export BinaryNode, Node

include("Tree.jl")
export BinaryTree, NamedBinaryTree, NamedTree
export PolytomousTree, NamedPolytomousTree
export hasrootheight, getrootheight, setrootheight!, clearrootheight!
export getrootdistance

include("routes.jl")
export branchhistory, branchroute, nodehistory, noderoute
export distance, distances, heighttoroot, heightstoroot

# Iterator methods expanded
include("Iterators.jl")
export nodeiter, nodefilter, nodenameiter, nodenamefilter,
    branchiter, branchfilter, branchnameiter, branchnamefilter

# A set of multiple trees
include("TreeSet.jl")
export TreeSet, treeiter, treenameiter, treeinfoiter

# Random tree generator
include("rand.jl")
export Nonultrametric, Ultrametric

# Read Newick Tree
include("newick.jl")
export parsenewick, parsenexus

# Display methods expanded
include("show.jl")

# Method for trimming trees
include("trim.jl")
export droptips!, keeptips!

# Path into package
path(path...; dir::String = "test") = joinpath(@__DIR__, "..", dir, path...)

using Requires
@static if VERSION < v"0.7.0-"
    @require RCall begin
        println("Creating Phylo RCall interface...")
        include("rcall.jl")
    end
else
    function __init__()
        @require RCall="6f49c342-dc21-5d91-9882-a32aef131414" begin
            println("Creating Phylo RCall interface...")
            include("rcall.jl")
        end
    end
end

end # module


text:
using Interpolations
using Test, LinearAlgebra

@testset "Scaling" begin
    # Model linear interpolation of y = -3 + .5x by interpolating y=x
    # and then scaling to the new x range

    itp = interpolate(1:1.0:10, BSpline(Linear()))

    sitp = @inferred(scale(itp, -3:.5:1.5))
    @test typeof(sitp) <: Interpolations.ScaledInterpolation
    @test parent(sitp) === itp

    for (x,y) in zip(-3:.05:1.5, 1:.1:10)
        @test sitp(x) ≈ y
    end

    # Verify that it works in >1D, with different types of ranges

    gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
    testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)

    xs = -5:.5:5
    ys = -4:.2:4
    zs = Float64[testfunction(x,y) for x in xs, y in ys]

    itp2 = interpolate(zs, BSpline(Quadratic(Flat(OnGrid()))))
    sitp2 = @inferred scale(itp2, xs, ys)

    for x in xs, y in ys
        @test testfunction(x,y) ≈ sitp2(x,y)
    end

    # Test gradients of scaled grids
    xs = -pi:.1:pi
    ys = map(sin, xs)
    itp = interpolate(ys, BSpline(Linear()))
    sitp = @inferred scale(itp, xs)

    for x in -pi:.1:pi
        g = @inferred(Interpolations.gradient(sitp, x))[1]
        @test ≈(cos(x),g,atol=0.05)
    end

    # Verify that return types are reasonable
    @inferred(sitp2(-3.4, 1.2))
    @inferred(sitp2(-3, 1))
    @inferred(sitp2(-3.4, 1))

    sitp32 = @inferred scale(interpolate(Float32[testfunction(x,y) for x in -5:.5:5, y in -4:.2:4], BSpline(Quadratic(Flat(OnGrid())))), -5f0:.5f0:5f0, -4f0:.2f0:4f0)
    @test typeof(@inferred(sitp32(-3.4f0, 1.2f0))) == Float32

    # Iteration
    itp = interpolate(rand(3,3,3), BSpline(Quadratic(Flat(OnCell()))))
    knots = map(d->1:10:21, 1:3)
    sitp = @inferred scale(itp, knots...)

    iter = @inferred(eachvalue(sitp))

    function foo!(dest, sitp)
        i = 0
        for s in eachvalue(sitp)
            dest[i+=1] = s
        end
        dest
    end
    function bar!(dest, sitp)
        for I in CartesianIndices(size(dest))
            dest[I] = sitp(I)
        end
        dest
    end
    rfoo = Array{Float64}(undef, Interpolations.ssize(sitp))
    rbar = similar(rfoo)
    foo!(rfoo, sitp)
    bar!(rbar, sitp)
    @test rfoo ≈ rbar
end