text:
using ClassifierMetrics
using Random
Random.seed!(123)

# Generate synthetic data
labels = rand(Bool, 500);

# function to add noise
add_noise(label::Bool, λ=0.0) = label ? 1 - λ*rand() : λ*rand()
# simulate good classifier
good = add_noise.(labels, 0.6)
# simulate bad classifier
bad  = add_noise.(labels, 0.98)

roc_good = ROC(labels, good)
roc_bad  = ROC(labels,  bad)

pr_bad  = PR(labels,  bad)
pr_good  = PR(labels,  good)

det_bad  = DET(labels,  bad)
det_good  = DET(labels,  good)

println("ROCAUC (good detector): $(auc(roc_good)) ")
println("ROCAUC (bad  detector): $(auc(roc_bad )) ")

println("PRAUC (good detector): $(auc(pr_good)) ")
println("PRAUC (bad  detector): $(auc(pr_bad )) ")

using Plots
gr()

roc_plot = plot(roc_good, label="good");
plot!(roc_plot,roc_bad, label="bad", legend=:bottomright)

pr_plot = plot(pr_good);
plot!(pr_plot,pr_bad, legend=false)

det_plot = plot(det_good);
plot!(det_plot,det_bad, legend=false)

plot(roc_plot, pr_plot, det_plot)


text:
# kspace_spiral.jl

using MIRT: mri_kspace_spiral

using Plots; default(markerstrokecolor=:auto, markersize=1, markershape=:circle)

N = 64

k0, o0, g0 = mri_kspace_spiral() # default 22,-1
for fov in (20,21)
	mri_kspace_spiral( ; fov=fov, Nt=-1, warn_nk=false)
end
k5l,_,g5l = mri_kspace_spiral(nl = 5) # interleaves

plot(xlabel="kx", ylabel="ky", aspect_ratio=1)
p1 = scatter!(k0[:,1], k0[:,2], label = "1-shot spiral")

p2 = plot()
scatter!(g0[:,1], label="gx")
scatter!(g0[:,2], label="gy")

p4 = plot(g5l[:,1,:], label="")
plot!(g5l[:,2,:], label="")

p3 = plot(xlabel="kx", ylabel="ky", aspect_ratio=1, title="5-shot spiral")
for ii=1:5
	scatter!(k5l[:,1,ii], k5l[:,2,ii], label="")
end

plot(p1, p2, p3, p4)

ext: jl


text:
using Mocking

get_bot_1_patch = @patch function get_bot(a...; b...)
    return Dict{String, Any}("name" => "DocOrderPizzaBot", "idleSessionTTLInSeconds" => 300, "checksum" => "20172ee3-fa06-49b2-bbc5-667c090303e9", "status" => "NOT_BUILT", "childDirected" => true, "lastUpdatedDate" => 1.494360160133e9, "abortStatement" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "I don't understand. Can you try again?", "contentType" => "PlainText"), Dict{String, Any}("content" => "I'm sorry, I don't understand.", "contentType" => "PlainText")]), "clarificationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "I'm sorry, I didn't hear that. Can you repeate what you just said?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Can you say that again?", "contentType" => "PlainText")], "maxAttempts" => 1), "version" => "\$LATEST", "createdDate" => 1.494360160133e9, "description" => "Orders a pizza from a local pizzeria.", "locale" => "en-US", "intents" => Any[Dict{String, Any}("intentVersion" => "\$LATEST", "intentName" => "DocOrderPizza")])
end

get_bots_1_patch = @patch function get_bots(a...; b...)
    return Dict{String, Any}("bots" => Any[Dict{String, Any}("name" => "DocOrderPizzaBot", "status" => "NOT_BUILT", "lastUpdatedDate" => 1.494360160133e9, "version" => "\$LATEST", "createdDate" => 1.494360160133e9, "description" => "Orders a pizza from a local pizzeria.")])
end

get_intent_1_patch = @patch function get_intent(a...; b...)
    return Dict{String, Any}("rejectionStatement" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "Ok, I'll cancel your order.", "contentType" => "PlainText"), Dict{String, Any}("content" => "I cancelled your order.", "contentType" => "PlainText")]), "name" => "DocOrderPizza", "checksum" => "ca9bc13d-afc8-4706-bbaf-091f7a5935d6", "sampleUtterances" => Any["Order me a pizza.", "Order me a {Type} pizza.", "I want a {Crust} crust {Type} pizza", "I want a {Crust} crust {Type} pizza with {Sauce} sauce."], "lastUpdatedDate" => 1.494359783453e9, "confirmationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "Should I order  your {Crust} crust {Type} pizza with {Sauce} sauce?", "contentType" => "PlainText")], "maxAttempts" => 1), "description" => "Order a pizza from a local pizzeria.", "version" => "\$LATEST", "createdDate" => 1.494359783453e9, "conclusionStatement" => Dict{String, Any}("responseCard" => "foo", "messages" => Any[Dict{String, Any}("content" => "All right, I ordered  you a {Crust} crust {Type} pizza with {Sauce} sauce.", "contentType" => "PlainText"), Dict{String, Any}("content" => "OK, your {Crust} crust {Type} pizza with {Sauce} sauce is on the way.", "contentType" => "PlainText")]), "fulfillmentActivity" => Dict{String, Any}("type" => "ReturnIntent"), "slots" => Any[Dict{String, Any}("name" => "Type", "priority" => 1, "valueElicitationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "What type of pizza would you like?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Vegie or cheese pizza?", "contentType" => "PlainText"), Dict{String, Any}("content" => "I can get you a vegie or a cheese pizza.", "contentType" => "PlainText")], "maxAttempts" => 1), "slotTypeVersion" => "\$LATEST", "sampleUtterances" => Any["Get me a {Type} pizza.", "A {Type} pizza please.", "I'd like a {Type} pizza."], "slotConstraint" => "Required", "description" => "The type of pizza to order.", "slotType" => "DocPizzaType"), Dict{String, Any}("name" => "Crust", "priority" => 2, "valueElicitationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "What type of crust would you like?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Thick or thin crust?", "contentType" => "PlainText")], "maxAttempts" => 1), "slotTypeVersion" => "\$LATEST", "sampleUtterances" => Any["Make it a {Crust} crust.", "I'd like a {Crust} crust."], "slotConstraint" => "Required", "description" => "The type of pizza crust to order.", "slotType" => "DocPizzaCrustType"), Dict{String, Any}("name" => "Sauce", "priority" => 3, "valueElicitationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "White or red sauce?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Garlic or tomato sauce?", "contentType" => "PlainText")], "maxAttempts" => 1), "slotTypeVersion" => "\$LATEST", "sampleUtterances" => Any["Make it {Sauce} sauce.", "I'd like {Sauce} sauce."], "slotConstraint" => "Required", "description" => "The type of sauce to use on the pizza.", "slotType" => "DocPizzaSauceType")])
end

get_intents_1_patch = @patch function get_intents(a...; b...)
    return Dict{String, Any}("intents" => Any[Dict{String, Any}("name" => "DocOrderPizza", "lastUpdatedDate" => 1.494359783453e9, "version" => "\$LATEST", "createdDate" => 1.494359783453e9, "description" => "Order a pizza from a local pizzeria.")])
end

get_slot_type_1_patch = @patch function get_slot_type(a...; b...)
    return Dict{String, Any}("name" => "DocPizzaCrustType", "checksum" => "210b3d5a-90a3-4b22-ac7e-f50c2c71095f", "lastUpdatedDate" => 1.494359274403e9, "enumerationValues" => Any[Dict{String, Any}("value" => "thick"), Dict{String, Any}("value" => "thin")], "version" => "\$LATEST", "createdDate" => 1.494359274403e9, "description" => "Available crust types")
end

get_slot_types_1_patch = @patch function get_slot_types(a...; b...)
    return Dict{String, Any}("slotTypes" => Any[Dict{String, Any}("name" => "DocPizzaCrustType", "lastUpdatedDate" => 1.494359274403e9, "version" => "\$LATEST", "createdDate" => 1.494359274403e9, "description" => "Available crust types"), Dict{String, Any}("name" => "DocPizzaSauceType", "lastUpdatedDate" => 1.49435644223e9, "version" => "\$LATEST", "createdDate" => 1.49435644223e9, "description" => "Available pizza sauces"), Dict{String, Any}("name" => "DocPizzaType", "lastUpdatedDate" => 1.494359198656e9, "version" => "\$LATEST", "createdDate" => 1.494359198656e9, "description" => "Available pizzas")])
end

put_bot_1_patch = @patch function put_bot(a...; b...)
    return Dict{String, Any}("name" => "DocOrderPizzaBot", "idleSessionTTLInSeconds" => 300, "checksum" => "20172ee3-fa06-49b2-bbc5-667c090303e9", "status" => "NOT_BUILT", "childDirected" => true, "lastUpdatedDate" => 1.494360160133e9, "abortStatement" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "I don't understand. Can you try again?", "contentType" => "PlainText"), Dict{String, Any}("content" => "I'm sorry, I don't understand.", "contentType" => "PlainText")]), "clarificationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "I'm sorry, I didn't hear that. Can you repeate what you just said?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Can you say that again?", "contentType" => "PlainText")], "maxAttempts" => 1), "version" => "\$LATEST", "createdDate" => 1.494360160133e9, "description" => "Orders a pizza from a local pizzeria.", "locale" => "en-US", "intents" => Any[Dict{String, Any}("intentVersion" => "\$LATEST", "intentName" => "DocOrderPizza")])
end

put_intent_1_patch = @patch function put_intent(a...; b...)
    return Dict{String, Any}("rejectionStatement" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "Ok, I'll cancel your order.", "contentType" => "PlainText"), Dict{String, Any}("content" => "I cancelled your order.", "contentType" => "PlainText")]), "name" => "DocOrderPizza", "checksum" => "ca9bc13d-afc8-4706-bbaf-091f7a5935d6", "sampleUtterances" => Any["Order me a pizza.", "Order me a {Type} pizza.", "I want a {Crust} crust {Type} pizza", "I want a {Crust} crust {Type} pizza with {Sauce} sauce."], "lastUpdatedDate" => 1.494359783453e9, "confirmationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "Should I order  your {Crust} crust {Type} pizza with {Sauce} sauce?", "contentType" => "PlainText")], "maxAttempts" => 1), "description" => "Order a pizza from a local pizzeria.", "version" => "\$LATEST", "createdDate" => 1.494359783453e9, "conclusionStatement" => Dict{String, Any}("responseCard" => "foo", "messages" => Any[Dict{String, Any}("content" => "All right, I ordered  you a {Crust} crust {Type} pizza with {Sauce} sauce.", "contentType" => "PlainText"), Dict{String, Any}("content" => "OK, your {Crust} crust {Type} pizza with {Sauce} sauce is on the way.", "contentType" => "PlainText")]), "fulfillmentActivity" => Dict{String, Any}("type" => "ReturnIntent"), "slots" => Any[Dict{String, Any}("name" => "Sauce", "priority" => 3, "valueElicitationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "White or red sauce?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Garlic or tomato sauce?", "contentType" => "PlainText")], "maxAttempts" => 1), "slotTypeVersion" => "\$LATEST", "sampleUtterances" => Any["Make it {Sauce} sauce.", "I'd like {Sauce} sauce."], "slotConstraint" => "Required", "description" => "The type of sauce to use on the pizza.", "slotType" => "DocPizzaSauceType"), Dict{String, Any}("name" => "Type", "priority" => 1, "valueElicitationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "What type of pizza would you like?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Vegie or cheese pizza?", "contentType" => "PlainText"), Dict{String, Any}("content" => "I can get you a vegie or a cheese pizza.", "contentType" => "PlainText")], "maxAttempts" => 1), "slotTypeVersion" => "\$LATEST", "sampleUtterances" => Any["Get me a {Type} pizza.", "A {Type} pizza please.", "I'd like a {Type} pizza."], "slotConstraint" => "Required", "description" => "The type of pizza to order.", "slotType" => "DocPizzaType"), Dict{String, Any}("name" => "Crust", "priority" => 2, "valueElicitationPrompt" => Dict{String, Any}("messages" => Any[Dict{String, Any}("content" => "What type of crust would you like?", "contentType" => "PlainText"), Dict{String, Any}("content" => "Thick or thin crust?", "contentType" => "PlainText")], "maxAttempts" => 1), "slotTypeVersion" => "\$LATEST", "sampleUtterances" => Any["Make it a {Crust} crust.", "I'd like a {Crust} crust."], "slotConstraint" => "Required", "description" => "The type of pizza crust to order.", "slotType" => "DocPizzaCrustType")])
end

put_slot_type_1_patch = @patch function put_slot_type(a...; b...)
    return Dict{String, Any}("name" => "DocPizzaSauceType", "checksum" => "cfd00ed1-775d-4357-947c-aca7e73b44ba", "lastUpdatedDate" => 1.49435644223e9, "enumerationValues" => Any[Dict{String, Any}("value" => "red"), Dict{String, Any}("value" => "white")], "version" => "\$LATEST", "createdDate" => 1.49435644223e9, "description" => "Available pizza sauces")
end


text:
using SacredDataStore
using Test

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


text:
using Random
using Test
using SSHT

bitsign(b::Bool) = b ? -1 : 1
bitsign(i::Integer) = bitsign(isodd(i))

chop(x) = abs2(x) < 100eps(x) ? zero(x) : x
chop(x::Complex) = Complex(chop(real(x)), chop(imag(x)))

function integrate(f::AbstractArray{T,2}, g::AbstractArray{U,2}, L::Integer) where {T,U}
    nphi = SSHT.sampling_dh_nphi(L)
    ntheta = SSHT.sampling_dh_ntheta(L)
    @assert size(f) == size(g) == (nphi, ntheta)

    s = zero(T) * zero(U) * zero(Float64)
    for p in 1:nphi, t in 1:ntheta
        phi = SSHT.sampling_dh_p2phi(p, L)
        theta = SSHT.sampling_dh_t2theta(t, L)
        dtheta = SSHT.sampling_weight_dh(theta, L)
        dphi = 2π / nphi
        s += conj(f[p, t]) * g[p, t] * dtheta * dphi
    end

    return s
end

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

# Half angle formulae:
#     sin(θ/2)^2 = (1-cos(θ))/2
#     cos(θ/2)^2 = (1+cos(θ))/2

# These sYlm are taken from Wikipedia and black-holes.org:
# <https://en.wikipedia.org/wiki/Spherical_harmonics> and
# <https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics>
# <https://www.black-holes.org/SpinWeightedSphericalHarmonics.nb>
function sYlm(s::Integer, l::Integer, m::Integer, θ::Real, ϕ::Real)
    @assert abs(s) ≤ l
    @assert -l ≤ m ≤ l

    # Parity:
    s < 0 && return bitsign(s + m) * conj(sYlm(-s, l, -m, θ, ϕ))

    (s, l, m) == (0, 0, 0) && return sqrt(1 / 4π)
    (s, l, m) == (0, 1, -1) && return sqrt(3 / 8π) * sin(θ) * cis(-ϕ)
    (s, l, m) == (0, 1, 0) && return sqrt(3 / 4π) * cos(θ)
    (s, l, m) == (0, 1, +1) && return -sqrt(3 / 8π) * sin(θ) * cis(ϕ)
    (s, l, m) == (0, 2, -2) && return sqrt(15 / 2π) * cos(θ / 2)^2 * sin(θ / 2)^2 * cis(-2ϕ)
    (s, l, m) == (0, 2, -1) && return -sqrt(15 / 2π) * cos(θ / 2) * sin(θ / 2) * (-cos(θ / 2)^2 + sin(θ / 2)^2) * cis(-ϕ)
    (s, l, m) == (0, 2, 0) && return sqrt(5 / 4π) * (cos(θ / 2)^4 - 4 * cos(θ / 2)^2 * sin(θ / 2)^2 + sin(θ / 2)^4)
    (s, l, m) == (0, 2, +1) && return -sqrt(15 / 2π) * cos(θ / 2) * sin(θ / 2) * (cos(θ / 2)^2 - sin(θ / 2)^2) * cis(ϕ)
    (s, l, m) == (0, 2, +2) && return sqrt(15 / 2π) * cos(θ / 2)^2 * sin(θ / 2)^2 * cis(2ϕ)

    (s, l, m) == (+1, 1, -1) && return -sqrt(3 / 16π) * (1 + cos(θ)) * cis(-ϕ)
    (s, l, m) == (+1, 1, 0) && return sqrt(3 / 8π) * sin(θ)
    (s, l, m) == (+1, 1, +1) && return -sqrt(3 / 16π) * (1 - cos(θ)) * cis(ϕ)
    (s, l, m) == (+1, 2, -2) && return -sqrt(5 / π) * cos(θ / 2)^3 * sin(θ / 2) * cis(-2ϕ)
    (s, l, m) == (+1, 2, -1) && return -sqrt(5 / 4π) * cos(θ / 2)^2 * (cos(θ / 2)^2 - 3 * sin(θ / 2)^2) * cis(-ϕ)
    (s, l, m) == (+1, 2, 0) && return sqrt(15 / 2π) * cos(θ / 2) * sin(θ / 2) * (cos(θ / 2)^2 - sin(θ / 2)^2)
    (s, l, m) == (+1, 2, +1) && return sqrt(5 / 4π) * sin(θ / 2)^2 * (-3 * cos(θ / 2)^2 + sin(θ / 2)^2) * cis(ϕ)
    (s, l, m) == (+1, 2, +2) && return sqrt(5 / π) * cos(θ / 2) * sin(θ / 2)^3 * cis(2ϕ)

    (s, l, m) == (+2, 2, -2) && return sqrt(5 / 4π) * cos(θ / 2)^4 * cis(-2ϕ)
    (s, l, m) == (+2, 2, -1) && return -sqrt(5 / π) * cos(θ / 2)^3 * sin(θ / 2) * cis(-ϕ)
    (s, l, m) == (+2, 2, 0) && return sqrt(15 / 2π) * cos(θ / 2)^2 * sin(θ / 2)^2
    (s, l, m) == (+2, 2, +1) && return -sqrt(5 / π) * cos(θ / 2) * sin(θ / 2)^3 * cis(ϕ)
    (s, l, m) == (+2, 2, +2) && return sqrt(5 / 4π) * sin(θ / 2)^4 * cis(2ϕ)

    @assert false
end

# ð F = - (sin θ)^s (∂_θ + i / sin(θ) ∂_ϕ) (sin θ)^-s F
# (Calculated manually from sYlm above)
function ðsYlm(s::Integer, l::Integer, m::Integer, θ::Real, ϕ::Real)
    @assert abs(s) ≤ l
    @assert -l ≤ m ≤ l

    (s, l, m) == (0, 0, 0) && return 0
    (s, l, m) == (0, 1, -1) && return -sqrt(3 / 8π) * (1 + cos(θ)) * cis(-ϕ)
    (s, l, m) == (0, 1, 0) && return sqrt(3 / 4π) * sin(θ)
    (s, l, m) == (0, 1, +1) && return -sqrt(3 / 8π) * (1 - cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 2, -2) && return -sqrt(15 / 8π) * (1 + cos(θ)) * sin(θ) * cis(-2ϕ)
    (s, l, m) == (0, 2, -1) && return -sqrt(15 / 2π) * cos(θ / 2)^2 * (-1 + 2 * cos(θ)) * cis(-ϕ)
    (s, l, m) == (0, 2, 0) && return sqrt(45 / 16π) * sin(2θ)
    (s, l, m) == (0, 2, +1) && return -sqrt(15 / 2π) * (1 + 2 * cos(θ)) * sin(θ / 2)^2 * cis(ϕ)
    (s, l, m) == (0, 2, +2) && return -sqrt(15 / 8π) * (-1 + cos(θ)) * sin(θ) * cis(2ϕ)

    (s, l, m) == (+1, 1, -1) && return 0
    (s, l, m) == (+1, 1, 0) && return 0
    (s, l, m) == (+1, 1, +1) && return 0
    (s, l, m) == (+1, 2, -2) && return sqrt(5 / π) * cos(θ / 2)^4 * cis(-2ϕ)
    (s, l, m) == (+1, 2, -1) && return -sqrt(5 / π) * cos(θ / 2)^2 * sin(θ) * cis(-ϕ)
    (s, l, m) == (+1, 2, 0) && return sqrt(15 / 8π) * sin(θ)^2
    (s, l, m) == (+1, 2, +1) && return -sqrt(5 / π) * sin(θ / 2)^2 * sin(θ) * cis(ϕ)
    (s, l, m) == (+1, 2, +2) && return sqrt(5 / π) * sin(θ / 2)^4 * cis(2ϕ)

    (s, l, m) == (+2, 2, -2) && return 0
    (s, l, m) == (+2, 2, -1) && return 0
    (s, l, m) == (+2, 2, 0) && return 0
    (s, l, m) == (+2, 2, +1) && return 0
    (s, l, m) == (+2, 2, +2) && return 0

    @assert false
end

# ð̄ F = - (sin θ)^-s (∂_θ - i / sin(θ) ∂_ϕ) (sin θ)^s F
# (Calculated manually from sYlm above)
function ð̄sYlm(s::Integer, l::Integer, m::Integer, θ::Real, ϕ::Real)
    @assert abs(s) ≤ l
    @assert -l ≤ m ≤ l

    (s, l, m) == (0, 0, 0) && return 0
    (s, l, m) == (0, 1, -1) && return sqrt(3 / 8π) * (1 - cos(θ)) * cis(-ϕ)
    (s, l, m) == (0, 1, 0) && return sqrt(3 / 4π) * sin(θ)
    (s, l, m) == (0, 1, +1) && return sqrt(3 / 8π) * (1 + cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 1, +1) && return -sqrt(3 / 8π) * (1 - cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 2, -2) && return -sqrt(15 / 8π) * (-1 + cos(θ)) * sin(θ) * cis(-2ϕ)
    (s, l, m) == (0, 2, -1) && return sqrt(15 / 2π) * (1 + 2 * cos(θ)) * sin(θ / 2)^2 * cis(-ϕ)
    (s, l, m) == (0, 2, 0) && return sqrt(45 / 16π) * sin(2θ)
    (s, l, m) == (0, 2, +1) && return sqrt(15 / 2π) * cos(θ / 2)^2 * (-1 + 2 * cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 2, +2) && return -sqrt(15 / 8π) * (1 + cos(θ)) * sin(θ) * cis(2ϕ)

    (s, l, m) == (+1, 1, -1) && return -sqrt(3 / 4π) * sin(θ) * cis(-ϕ)
    (s, l, m) == (+1, 1, 0) && return -sqrt(3 / 2π) * cos(θ)
    (s, l, m) == (+1, 1, +1) && return sqrt(3 / 4π) * sin(θ) * cis(ϕ)
    (s, l, m) == (+1, 2, -2) && return -sqrt(45 / 16π) * sin(θ)^2 * cis(-2ϕ)
    (s, l, m) == (+1, 2, -1) && return -sqrt(45 / 16π) * sin(2θ) * cis(-ϕ)
    (s, l, m) == (+1, 2, 0) && return -sqrt(15 / 32π) * (1 + 3 * cos(2θ))
    (s, l, m) == (+1, 2, +1) && return sqrt(45 / 16π) * sin(2θ) * cis(ϕ)
    (s, l, m) == (+1, 2, +2) && return -sqrt(45 / 16π) * sin(θ)^2 * cis(2ϕ)

    (s, l, m) == (+2, 2, -2) && return sqrt(5 / 4π) * (1 + cos(θ)) * sin(θ) * cis(-2ϕ)
    (s, l, m) == (+2, 2, -1) && return sqrt(5 / 4π) * (cos(θ) + cos(2θ)) * cis(-ϕ)
    (s, l, m) == (+2, 2, 0) && return -sqrt(15 / 2π) * cos(θ) * sin(θ)
    (s, l, m) == (+2, 2, +1) && return sqrt(5 / 4π) * (cos(θ) - cos(2θ)) * cis(ϕ)
    (s, l, m) == (+2, 2, +2) && return -sqrt(5 / π) * sin(θ / 2)^2 * sin(θ) * cis(2ϕ)

    @assert false
end

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

@testset "Grid point counts (Driscoll & Healy)" begin
    for L in 1:10
        n = SSHT.sampling_dh_n(L)
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)
        @test nphi == 2 * L - 1
        @test ntheta == 2 * L
        @test n == nphi * ntheta
    end
end

@testset "Coordinates (Driscoll & Healy)" begin
    for L in 1:100
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)
        for p in 1:nphi
            phi = SSHT.sampling_dh_p2phi(p, L)
            @test phi ≈ 2π * (p - 1) / nphi
        end
        for t in 1:ntheta
            theta = SSHT.sampling_dh_t2theta(t, L)
            @test theta ≈ π * (t - 1 / 2) / ntheta
        end
    end
end

@testset "Grid point counts (McEwen & Wiaux)" begin
    for L in 1:2:10
        n = SSHT.sampling_mw_n(L)
        nphi = SSHT.sampling_mw_nphi(L)
        ntheta = SSHT.sampling_mw_ntheta(L)
        @test nphi == 2 * L - 1
        @test ntheta == L
        # The south pole is covered only by a single point
        @test n == nphi * (ntheta - 1) + 1
    end
end

@testset "Coordinates (McEwen & Wiaux)" begin
    for L in 1:100
        nphi = SSHT.sampling_mw_nphi(L)
        ntheta = SSHT.sampling_mw_ntheta(L)
        for p in 1:nphi
            phi = SSHT.sampling_mw_p2phi(p, L)
            @test phi ≈ 2π * (p - 1) / nphi
        end
        for t in 1:ntheta
            theta = SSHT.sampling_mw_t2theta(t, L)
            @test theta ≈ π * (2 * t - 1) / (2 * ntheta - 1)
        end
    end
end

@testset "Mode layout and counts" begin
    for L in 1:10
        mode_seen = falses(L^2)
        for l in 0:(L - 1), m in (-l):l
            ind = SSHT.sampling_elm2ind(l, m)
            @test 1 ≤ ind ≤ L^2
            @test SSHT.sampling_ind2elm(ind) == (l, m)
            @test !mode_seen[ind]
            mode_seen[ind] = true
        end
        @test all(mode_seen)
    end
end

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

# TODO: Test phase and parity of sYlm
#
# {\displaystyle Y_{\ell }^{m}(\theta ,\phi )\to Y_{\ell }^{m}(\pi -\theta ,\pi +\phi )=(-1)^{\ell }Y_{\ell }^{m}(\theta ,\phi )}
#
# <math display="block">Y_\ell^m(\theta,\phi) \to Y_\ell^m(\pi-\theta,\pi+\phi) = (-1)^\ell Y_\ell^m(\theta,\phi)</math>
#
# {}_s\bar Y_{l m} &= \left(-1\right)^{s+m}{}_{-s}Y_{l(-m)}\\
#
# {}_sY_{l m}(\pi-\theta,\phi+\pi) &= \left(-1\right)^l {}_{-s}Y_{l m}(\theta,\phi).

Random.seed!(100)
modes = [(name="(l=0,m=0)", fun=(θ, ϕ) -> sYlm(0, 0, 0, θ, ϕ), modes=[1, 0, 0, 0]),
         (name="(l=1,m=-1)", fun=(θ, ϕ) -> im * sYlm(0, 1, +1, θ, ϕ) + im * sYlm(0, 1, -1, θ, ϕ), modes=[0, im, 0, im]),
         (name="(l=1,m=0)", fun=(θ, ϕ) -> sYlm(0, 1, 0, θ, ϕ), modes=[0, 0, 1, 0]),
         (name="(l=1,m=+1)", fun=(θ, ϕ) -> sYlm(0, 1, +1, θ, ϕ) - sYlm(0, 1, -1, θ, ϕ), modes=[0, -1, 0, 1])]
@testset "Simple real transforms: $(mode.name)" for mode in modes
    for L in 2:20
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)

        f = Array{Float64}(undef, nphi, ntheta)
        for p in 1:nphi, t in 1:ntheta
            ϕ = SSHT.sampling_dh_p2phi(p, L)
            θ = SSHT.sampling_dh_t2theta(t, L)
            f[p, t] = mode.fun(θ, ϕ)
        end

        flm = SSHT.core_dh_forward_sov_real(f, L)
        @assert L ≥ 2
        @test flm[1:4] ≈ mode.modes
        if L > 2
            @test all(isapprox(0; atol=100eps()), flm[5:end])
        end

        f′ = SSHT.core_dh_inverse_sov_real(flm, L)
        @test isapprox(f′, f; atol=1000eps())
    end
end

Random.seed!(100)
modes = [(name="(s=$s,l=$l,m=$m)", spin=s, el=l, fun=(θ, ϕ) -> sYlm(s, l, m, θ, ϕ),
          modes=L -> [l == l′ && m == m′ for l′ in 0:(L - 1) for m′ in (-l′):l′]) for s in -2:+2 for l in abs(s):2 for m in (-l):l]
@testset "Simple complex transforms: $(mode.name)" for mode in modes
    for L in (mode.el + 1):20
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)

        f = Array{Complex{Float64}}(undef, nphi, ntheta)
        for p in 1:nphi, t in 1:ntheta
            ϕ = SSHT.sampling_dh_p2phi(p, L)
            θ = SSHT.sampling_dh_t2theta(t, L)
            f[p, t] = mode.fun(θ, ϕ)
        end

        flm = SSHT.core_dh_forward_sov(f, L, mode.spin)
        @test flm ≈ mode.modes(L)

        f′ = SSHT.core_dh_inverse_sov(flm, L, mode.spin)
        @test isapprox(f′, f; atol=1000eps())
    end
end

Random.seed!(100)
@testset "Linearity of real transforms" begin
    for iter in 1:100
        L = rand(1:100)
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)

        f = randn(Float64, nphi, ntheta)
        g = randn(Float64, nphi, ntheta)
        α = randn(Float64)
        h = f + α * g

        flm = SSHT.core_dh_forward_sov_real(f, L)
        glm = SSHT.core_dh_forward_sov_real(g, L)
        hlm = SSHT.core_dh_forward_sov_real(h, L)

        @test flm + α * glm ≈ hlm
    end
end

Random.seed!(100)
@testset "Linearity of complex transforms" begin
    for iter in 1:100
        L = rand(1:100)
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)

        spin = rand(-4:4)

        f = randn(Complex{Float64}, nphi, ntheta)
        g = randn(Complex{Float64}, nphi, ntheta)
        α = randn(Complex{Float64})
        h = f + α * g

        flm = SSHT.core_dh_forward_sov(f, L, spin)
        glm = SSHT.core_dh_forward_sov(g, L, spin)
        hlm = SSHT.core_dh_forward_sov(h, L, spin)

        @test flm + α * glm ≈ hlm

        flm = randn(Complex{Float64}, L^2)
        glm = randn(Complex{Float64}, L^2)
        hlm = flm + α * glm

        f = SSHT.core_dh_inverse_sov(flm, L, spin)
        g = SSHT.core_dh_inverse_sov(glm, L, spin)
        h = SSHT.core_dh_inverse_sov(hlm, L, spin)

        @test f + α * g ≈ h
    end
end

Random.seed!(100)
@testset "Orthonormality of complex transforms" begin
    for iter in 1:100
        L = rand(1:100)
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)

        smax = min(4, L - 1)
        spin = rand((-smax):smax)

        flm = zeros(Complex{Float64}, L^2)
        glm = zeros(Complex{Float64}, L^2)

        lf = rand(abs(spin):(L - 1))
        mf = rand((-lf):lf)
        lg = rand(abs(spin):(L - 1))
        mg = rand((-lg):lg)

        flm[SSHT.sampling_elm2ind(lf, mf)] = 1
        glm[SSHT.sampling_elm2ind(lg, mg)] = 1

        f = SSHT.core_dh_inverse_sov(flm, L, spin)
        g = SSHT.core_dh_inverse_sov(glm, L, spin)

        @test isapprox(integrate(f, f, L), 1; atol=1 / L^2)
        @test isapprox(integrate(f, g, L), (lf == lg) * (mf == mg); atol=1 / L^2)

        h = conj(f) .* f
        hlm = SSHT.core_dh_forward_sov(h, L, 0)
        @test isapprox(hlm[SSHT.sampling_elm2ind(0, 0)], sqrt(1 / 4π); atol=sqrt(eps()))

        h = conj(f) .* g
        hlm = SSHT.core_dh_forward_sov(h, L, 0)
        @test isapprox(hlm[SSHT.sampling_elm2ind(0, 0)], (lf == lg) * (mf == mg) * sqrt(1 / 4π); atol=sqrt(eps()))
    end
end

Random.seed!(100)
modes = [(name="(s=$s,l=$l,m=$m)", spin=s, el=l, fun=(θ, ϕ) -> sYlm(s, l, m, θ, ϕ), ðfun=(θ, ϕ) -> ðsYlm(s, l, m, θ, ϕ),
          ð̄fun=(θ, ϕ) -> ð̄sYlm(s, l, m, θ, ϕ)) for s in 0:+2 for l in abs(s):2 for m in (-l):l]
@testset "Simple derivatives (eth, eth-bar): $(mode.name)" for mode in modes
    for L in (mode.el + 1):20
        nphi = SSHT.sampling_dh_nphi(L)
        ntheta = SSHT.sampling_dh_ntheta(L)

        f = Array{Complex{Float64}}(undef, nphi, ntheta)
        ðf₀ = Array{Complex{Float64}}(undef, nphi, ntheta)
        ð̄f₀ = Array{Complex{Float64}}(undef, nphi, ntheta)
        for p in 1:nphi, t in 1:ntheta
            ϕ = SSHT.sampling_dh_p2phi(p, L)
            θ = SSHT.sampling_dh_t2theta(t, L)
            f[p, t] = mode.fun(θ, ϕ)
            ðf₀[p, t] = mode.ðfun(θ, ϕ)
            ð̄f₀[p, t] = mode.ð̄fun(θ, ϕ)
        end

        flm = SSHT.core_dh_forward_sov(f, L, mode.spin)

        ðflm = SSHT.eth(flm, L, mode.spin)
        ðf = SSHT.core_dh_inverse_sov(ðflm, L, mode.spin + 1)
        @test isapprox(ðf, ðf₀; atol=10000eps())

        ð̄flm = SSHT.ethbar(flm, L, mode.spin)
        ð̄f = SSHT.core_dh_inverse_sov(ð̄flm, L, mode.spin - 1)
        @test isapprox(ð̄f, ð̄f₀; atol=10000eps())
    end
end


text:
import FirstOrderLp
import CSV

include("../src/PrimalDualMethods.jl")
include("../src/shared_variables.jl")

function find_best_primal_weight(
  params::PrimalDualOptimizerParameters,
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  primal_weight_trial_list::Vector{Float64},
  primal_weight_for_evaluation::Float64,
  step_size_for_evaluation::Float64,
)
  best_kkt_error = Inf
  best_primal_weight = NaN
  for primal_weight in primal_weight_trial_list
    println("run for primal weight = $primal_weight")
    params.primal_weight = primal_weight
    solver_output = optimize(params, problem)

    kkt_error = min(
      solver_output.iteration_stats.kkt_error_average_iterate[end],
      solver_output.iteration_stats.kkt_error_current_iterate[end],
    )

    println("kkt_error = $kkt_error")
    if kkt_error < best_kkt_error
      best_kkt_error = kkt_error
      best_primal_weight = primal_weight
    end
  end
  return best_primal_weight
end

mutable struct ResultsCsvPaths
  no_restarts::String
  adaptive_restarts::String
  dynamic_adaptive_restarts::String
  restart_lengths::Dict{Int64,String}
end

function get_csv_paths(
  results_directory::String,
  restart_length_list::Vector{Int64},
)
  restart_length_csv_paths = Dict{Int64,String}()
  for restart_length in restart_length_list
    restart_length_csv_paths[restart_length] = "$results_directory/restart_length$(restart_length).csv"
  end
  return ResultsCsvPaths(
    "$results_directory/no_restarts.csv",
    "$results_directory/adaptive_restarts.csv",
    "$results_directory/dynamic_adaptive_restarts.csv",
    restart_length_csv_paths,
  )
end

function solve_problem_with_multiple_methods(;
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  restart_length_list::Vector{Int64},
  results_directory::String,
  iteration_limit_for_finding_best_primal_weight::Int64,
  iteration_limit::Int64,
  kkt_tolerance::Float64,
  method::PrimalDualMethod,
)
  println("solving problem with: ")
  print("rows = ", size(problem.constraint_matrix, 1), ", ")
  print("cols = ", size(problem.constraint_matrix, 2), ", ")
  println("nnz = ", length(problem.constraint_matrix.nzval), ".")

  saved_stdout = stdout
  log_file = open("$results_directory/best_primal_weight.txt", "w")
  redirect_stdout(log_file)
  params = PrimalDualOptimizerParameters(
    method,
    nothing, # step_size (forces the solver to use a provably correct step size)
    1.0, # primal_weight
    30, # record every
    100, # print every
    true, # verbose
    iteration_limit_for_finding_best_primal_weight, # iteration limit
    NoRestarts(), # restart scheme
    nothing, # don't record information on subspace
    kkt_tolerance,
  )

  println("FIND BEST PRIMAL WEIGHT")
  best_primal_weight =
    find_best_primal_weight(params, problem, 4.0 .^ collect(-5:5), 1.0, 1.0)
  println("best_primal_weight = $best_primal_weight")
  close(log_file)

  csv_paths = get_csv_paths(results_directory, restart_length_list)

  log_file = open("$results_directory/no_restarts.txt", "w")
  redirect_stdout(log_file)
  params.iteration_limit = iteration_limit
  params.primal_weight = best_primal_weight
  println("Run no restarts")
  solver_output = optimize(params, problem)
  close(log_file)
  CSV.write(csv_paths.no_restarts, solver_output.iteration_stats)

  log_file = open("$results_directory/adaptive_restarts.txt", "w")
  redirect_stdout(log_file)
  println("Run adaptive restarts")
  params.restart_scheme = AdaptiveRestarts(exp(-1), true)
  solver_output = optimize(params, problem)
  close(log_file)
  CSV.write(csv_paths.adaptive_restarts, solver_output.iteration_stats)

  log_file = open("$results_directory/dynamic_adaptive_restarts.txt", "w")
  redirect_stdout(log_file)
  println("Run extra adaptive restarts")
  params.restart_scheme = AdaptiveRestarts(exp(-1), false)
  solver_output = optimize(params, problem)
  close(log_file)
  CSV.write(csv_paths.dynamic_adaptive_restarts, solver_output.iteration_stats)

  for restart_length in restart_length_list
    log_file =
      open("$results_directory/restart_length_$(restart_length).txt", "w")
    redirect_stdout(log_file)
    println("restart_length = $restart_length")
    params.restart_scheme = FixedFrequencyRestarts(restart_length, true)
    solver_output = optimize(params, problem)
    close(log_file)
    CSV.write(
      csv_paths.restart_lengths[restart_length],
      solver_output.iteration_stats,
    )
  end
  redirect_stdout(saved_stdout)
end

@assert length(ARGS) == 4
test_problem_folder = ARGS[1]
results_directory = ARGS[2]
method = RecoverPrimalDualMethodFromString(ARGS[3])
problem_name_arg = ARGS[4]

function main()
  ################
  # run problems #
  ################
  iteration_limit_for_finding_best_primal_weight = 5000
  kkt_tolerance = 1e-9

  if problem_name_arg == "all"
    problem_name_list = ALL_PROBLEM_NAMES
  elseif problem_name_arg in ALL_PROBLEM_NAMES
    problem_name_list = [problem_name_arg]
  else
    error("unknown problem name arguement supplied")
  end

  for problem_name in problem_name_list
    instance_path = joinpath(test_problem_folder, "$(problem_name).mps.gz")
    lp = FirstOrderLp.qps_reader_to_standard_form(instance_path)
    lp.variable_upper_bound .= Inf

    results_subdirectory = joinpath(results_directory, problem_name)
    mkpath(results_subdirectory)

    solve_problem_with_multiple_methods(;
      problem = lp,
      restart_length_list = RESTART_LENGTHS_DICT[problem_name],
      results_directory = results_subdirectory,
      iteration_limit_for_finding_best_primal_weight = iteration_limit_for_finding_best_primal_weight,
      iteration_limit = ITERATION_LIMIT,
      kkt_tolerance = kkt_tolerance,
      method,
    )
  end
end

main()

max_stars_repo_name: pedersor/google-research

max_issues_repo_name: pedersor/google-research

max_forks_repo_name: pedersor/google-research


text:
using Documenter, StatsBase, Statistics, Random

makedocs(
    format = :html,
    sitename = "StatsBase.jl",
    modules = [StatsBase],
    pages = ["index.md",
             "weights.md",
             "means.md",
             "scalarstats.md",
             "robust.md",
             "deviation.md",
             "cov.md",
             "counts.md",
             "ranking.md",
             "sampling.md",
             "empirical.md",
             "signalcorr.md",
             "multivariate.md",
             "misc.md",
             "statmodels.md",
             "transformations.md"]
)

deploydocs(
    repo = "github.com/JuliaStats/StatsBase.jl.git",
    target = "build",
    julia  = "1.0",
    deps   = nothing,
    make   = nothing
)


text:
"""
Subroutine dgpstep

This subroutine computes the gradient projection step

  s = P[x + alpha*w] - x,

where P is the projection on the n-dimensional interval [xl,xu].
"""
function dgpstep(n,x,xl,xu,alpha,w,s)
    # TODO
    # This computation of the gradient projection step avoids
    # rouding errors for the components that are feasible.

    @inbounds for i=1:n
        if x[i] + alpha*w[i] < xl[i]
            s[i] = xl[i] - x[i]
        elseif x[i] + alpha*w[i] > xu[i]
            s[i] = xu[i] - x[i]
        else
            s[i] = alpha*w[i]
        end
    end
    return
end

@inline function dgpstep(n::Int,x::CuDeviceArray{Float64,1},xl::CuDeviceArray{Float64,1},
                         xu::CuDeviceArray{Float64,1},alpha,w::CuDeviceArray{Float64,1},
                         s::CuDeviceArray{Float64,1})
    tx = threadIdx().x
    ty = threadIdx().y

    if tx <= n && ty == 1
        @inbounds begin
            # It might be better to process this using just a single thread,
            # rather than diverging between multiple threads.

            if x[tx] + alpha*w[tx] < xl[tx]
                s[tx] = xl[tx] - x[tx]
            elseif x[tx] + alpha*w[tx] > xu[tx]
                s[tx] = xu[tx] - x[tx]
            else
                s[tx] = alpha*w[tx]
            end
        end
    end
    CUDA.sync_threads()

    return
end


text:
__precompile__()
module gauss_quadrature
using Dierckx
export  do_quad, delta_measure,Gauss_Hermite, Gauss_Hermite_per,QUAD, Golub_Welsch!,periodic_simplify!,simplify!,sortquad!,quad_to_pdf
##########################
function init(fcn, opt)
  #=
  =#
end
#######################
type QUAD
  length
  X
  w
end
#######################
function do_quad(test,quad::QUAD)
  # evaluate test fn wrt to quad
  #
  sum=0.0
  for i=1:quad.length
    sum += test(quad.X[i]) * quad.w[i]
  end
  return sum
end
##############################################
function lanczos(N,xw)
  # Adapted from Gautschi's OPQ Matlab routines
  Ncap=size(xw,1)
  if N<=0 || N>Ncap
      print("N:",N,"; Ncap:",Ncap)
     error("N out of range")
  end
  p0=xw[:,1]; p1=zeros(Ncap,1); p1[1]=xw[1,2]
  for n=1:Ncap-1
    pn=xw[n+1,2]; gam=1; sig=0; t=0; xlam=xw[n+1,1]
    for k=1:n+1
      rho=p1[k]+pn; tmp=gam*rho; tsig=sig

      if (rho)<=0
        gam=1; sig=0
      else
        gam=p1[k]/rho
        sig=pn/rho
      end
      tk=sig*(p0[k]-xlam)-gam*t
      p0[k]=p0[k]-(tk-t); t=tk
      if (sig)<=0
        pn=tsig*p1[k]
      else
        pn=(t^2)/sig
      end
      tsig=sig; p1[k]=tmp
    end
  end
  ab=[p0[1:N] p1[1:N]]
  return ab
end
######################################
function gauss(N,ab)
  # Adapted from Gautschi's OPQ Matlab routines
  N0=size(ab,1)
  if N0<N
    error("input array ab too short")
  end
  J=zeros(N,N)
  for n=1:N
     J[n,n]=ab[n,1]
  end
  for n=2:N
    J[n,n-1]=sqrt(ab[n,2])
    J[n-1,n]=J[n,n-1]
  end
  D,V=eig(J)
  #V=V(:,I);
  xw=[D ab[1,2]*V[1,:]'.^2]
  return xw
end
########################################
function Golub_Welsch!(quad::QUAD, m)
  # Compute m-point GQ
  # No change if input quad has too few points
  if quad.length> m
    m=round(Int16,m)
    ab=lanczos(m, [quad.X quad.w])
    xw=gauss(m,ab)
    quad.length=m
    quad.X=xw[:,1]
    quad.w=xw[:,2]
  end
end
###########################
function sortquad!(quad::QUAD)
  # sort quad points
  perm=sortperm(quad.X)
  quad.X=quad.X[perm]
  quad.w=quad.w[perm]
end
###########################
function simplify!(quad::QUAD, m, cutoff)
# Return reduced quadrature
# Group points outside interval [cutoff[1],cutoff[2]]
# Apply m-point GQ on interior
# No action if input quad less than m
  if quad.length > m
     ind=find(x->(x<=cutoff[1]), quad.X) # small nodes
     m1=length(ind); w1=sum(quad.w[ind])
     ind=find(x->(x>=cutoff[2]), quad.X) # big nodes
     m2=length(ind); w2=sum(quad.w[ind])
     if (quad.length-m1-m2>m)
       ind=find(y->(y>cutoff[1] && y<cutoff[2]), quad.X) # interior nodes
       ab=lanczos(m, [quad.X[ind] quad.w[ind]])
       xw=gauss(m,ab)
       quad.length=m+2
       quad.X=[cutoff[1]; xw[:,1]; cutoff[2]]
       quad.w=[w1;        xw[:,2]; w2] # update quad
     end
  end
end
#################################################
function partition_simplify!(quad::QUAD, m, d, cutoff)
# Apply periodic conditions and find m-point GQ
# on d sub-intervals.
  if quad.length > m
   ind=find(x->(x<cutoff[1]), quad.X) # small nodes
   m1=length(ind); w1=sum(quad.w[ind])
   ind=find(x->(x>cutoff[2]), quad.X) # big nodes
   m2=length(ind); w2=sum(quad.w[ind])
   if (quad.length-m1-m2>m)
    partition=collect(linspace(cutoff[1], cutoff[2], d+1))
    x=[]; w=[];
    for i=1:d
      inds=find(y->(y>partition[i] && y<=partition[i+1]), quad.X)
      quad1=QUAD(length(inds),quad.X[inds],quad.w[inds])
      Golub_Welsch!(quad1, max(5,m*min(1,10*sum(quad1.w) ) ) )
      x=[x; quad1.X]
      w=[w; quad1.w]
    end
    x=[cutoff[1]; x; cutoff[2]]
    w=[w1; w; w2]
    quad.length=length(x)
    quad.X=x
    quad.w=w
  end
 end
end
#####
function l1(q1::QUAD, q2::QUAD)
  mysum=0.
  i1=2; i2=2; c1=cumsum(q1.w); c2=cumsum(q2.w)
  if (q1.X[i1] < q2.X[i2]) # find smallest point
    b=q1.X[i1];  i1=i1+1
  else
    b=q2.X[i2];  i2=i2+1
  end
  while (i1<=q1.length && i2<=q2.length)
    if (q1.X[i1]<q2.X[i2]) # q1 is next smallest
      delta=q1.X[i1]-b; # work out delta X
      # interpolate on q2
      if (q2.X[i2]>q2.X[i2-1])
        lambda=(q1.X[i1]-q2.X[i2-1])/(q2.X[i2]-q2.X[i2-1])
        cinterp=lambda*c2[i2]+(1-lambda)*c2[i2-1]#
        c_delta=abs(cinterp-c1[i1]);
      else
        c_delta=0
      end
      # record current x
      b=q1.X[i1]; i1+=1
    else
      delta=q2.X[i2]-b; # work out delta X
      # interpolate on q2
      if q1.X[i1-1]<q1.X[i1]
        lambda=(q2.X[i2]-q1.X[i1-1])/(q1.X[i1]-q1.X[i1-1])
        cinterp=lambda*c1[i1]+(1-lambda)*c1[i1-1]
        c_delta=abs(cinterp-c2[i2])
      else
        c_delta=0
      end
      # record current x
      b=q2.X[i2]; i2+=1
    end
    mysum+=c_delta*delta
  end
  if (i1<q1.length)
    mysum+=sum(abs(c1[i1:end-1]).*(q1.w[i1+1:end]-q1.w[i1:end-1]))
  end
  if (i2<q2.length)
    mysum+=sum(abs(c2[i2:end-1]).*(q2.w[i2+1:end]-q2.w[i2:end-1]))
  end
  return mysum
end
#

#####
function l1cdf(q1::QUAD, cdf)
  mysum=0.
  c1=cumsum(q1.w)
  for i=2:q1.length,
    delta1=abs(c1[i-1]-cdf(q1.X[i-1] ))
    delta2=abs(c1[i]  -cdf(q1.X[i]) )
    mysum+= (q1.X[i]-q1.X[i-1])*(delta1+delta2)/2
  end
  return mysum
end
#################################################
function periodic_simplify!(quad::QUAD, m, length_scale)
# Apply periodic conditions and find m-point GQ
  quad.X=mod(quad.X, length_scale)
  Golub_Welsch!(quad, m)
end
#################################################
function partition_periodic_simplify!(quad::QUAD, m, d, length_scale)
# Apply periodic conditions and find m-point GQ
# on d sub-intervals.
  theta=mod(quad.X, length_scale)
  partition=collect(linspace(-1e-6, length_scale, d+1))
  x=[]; w=[];
  for i=1:d
    inds=find(y->(y>partition[i] && y<=partition[i+1]), theta)
    quad1=QUAD(length(inds),theta[inds],quad.w[inds])
    Golub_Welsch!(quad1, m)
    x=[x; quad1.X]
    w=[w; quad1.w]
  end
  quad.length=length(x)
  quad.X=x
  quad.w=w
end
###########################
function Gauss_Hermite(mu_sig_sq, n)
  # Returns m-point Gauss-Hermite quadrature rule
  # for N(m,s2), for mu_sig_sq=[m,s2]
  #
  # using Hermite recursion rule for w(x)=e^(-x^2)
  # and change of variable
  #
  println("GH: [mean,variance]=",mu_sig_sq)
  v =sqrt((1:n-1)/2.0);
  J=diagm(v,1)+diagm(v,-1)
  D,V=eig(J)
  xw=[D sqrt(pi)*V[1,:]'.^2]
  # change of variable
  f=x->mu_sig_sq[1]+x*sqrt(2*mu_sig_sq[2])
  weight_scale=sqrt(pi)
  quad=QUAD(n,f(xw[:,1]),xw[:,2]/weight_scale)
  #
  #
  println(quad.length,"-pt Gauss-Hermite rule.")
  sortquad!(quad)
  return quad
end
###########################
function Gauss_Hermite_per(mu_sig_sq, m, length_scale)
  # Returns m-point Gauss-Hermite quadrature rule
  # for N(m,s2), for mu_sig_sq=[m,s2]
  #
  # using Hermite recursion rule for w(x)=e^(-x^2)
  # and change of variable
  #
  # Apply periodic conditions on domain with length_scale length
  println("GH: [mean,variance]=",mu_sig_sq)
  nn=10*m
  v =sqrt((1:nn-1)/2.0);
  J=diagm(v,1)+diagm(v,-1)
  D,V=eig(J)
  xw=[D sqrt(pi)*V[1,:]'.^2]
  # change of variable
  f=x->mu_sig_sq[1]+x*sqrt(2*mu_sig_sq[2])
  weight_scale=sqrt(pi)
  quad=QUAD(nn,f(xw[:,1]),xw[:,2]/weight_scale)
  #
  periodic_simplify!(quad, m, length_scale)
  #
  println(quad.length,"-pt Gauss-Hermite rule.")
  sortquad!(quad)
  return quad
end
################################
function Gauss_Hermite_sech(mu_nu, n)
  # Returns m-point Gauss-Hermite quadrature rule
  # for sech^2( (x-mu)/(4 nu))/(8 nu)
  # for mu_nu=[mu,nu]
  # via Gauss-Hermite
  #
  nu=mu_nu[2]
  mu=mu_nu[1]
  s=nu*sqrt(8)
  mu_sig_sq=[mu, s^2]
  quad=Gauss_Hermite(mu_sig_sq, n)
  c=4*sqrt(pi*nu^2);
  helper=z-> c/(4*nu* (exp(-0.25*z/nu)+exp(0.25*z/nu))*exp(-(0.25*z/nu)^2) )
  for i=1:quad.length
    quad.w[i]=quad.w[i]*helper(quad.X[i]-mu)
  end
  quad.w=quad.w/sum(quad.w)
  return quad
end
#############################
function delta_measure(x)
  return  QUAD(1,x,1)
end
######################
function quad_to_pdf(quad::QUAD,n,s_in)
  sortquad!(quad)
  cdf_all=cumsum(quad.w)
  kk=3
  cdf_all=[zeros(kk);cdf_all;ones(kk)]
  eps=1e-2
  X=[quad.X[1]+linspace(-kk*eps,-eps,kk);quad.X;quad.X[end]+linspace(eps,kk*eps,kk)]
  n=min(quad.length,n)
  k=3
  range=collect(linspace(quad.X[1],quad.X[end],n))
  #print(length(X)," ", length(cdf_all), " ", length(range))
  #spl=Spline1D(X,cdf_all,range;w=ones(length(X)), k=k, bc="nearest")
  spl=Spline1D(X,cdf_all;w=ones(length(X)), k=k, bc="nearest",s=s_in)
  pdf=derivative(spl,range)
  cdf=spl(range)
  return range,pdf,cdf,X,cdf_all
end

######################
function quad_to_pdf(quad::QUAD)
  sortquad!(quad)
  cdf_all=cumsum(quad.w)
  kk=3
  cdf_all=[zeros(kk);cdf_all;ones(kk)]
  eps=1e-2
  X=[quad.X[1]+linspace(-kk*eps,-eps,kk);quad.X;quad.X[end]+linspace(eps,kk*eps,kk)]
  n=min(quad.length,n)
  k=3
  range=collect(linspace(quad.X[1],quad.X[end],n))
  #print(length(X)," ", length(cdf_all), " ", length(range))
  #spl=Spline1D(X,cdf_all,range;w=ones(length(X)), k=k, bc="nearest")
  spl=Spline1D(X,cdf_all)
  pdf=derivative(spl,range)
  cdf=spl(range)
  return range,pdf,cdf,X,cdf_all
end
#
end # module defn


text:
using Documenter, RipQP

makedocs(
  modules = [RipQP],
  doctest = true,
  linkcheck = true,
  strict = true,
  format = Documenter.HTML(assets = ["assets/style.css"], prettyurls = get(ENV, "CI", nothing) == "true"),
  sitename = "RipQP.jl",
  pages = ["Home" => "index.md",
           "Tutorial" => "tutorial.md",
           "API" => "api.md",
           "Reference" => "reference.md",
          ]
)

deploydocs(deps = nothing, make = nothing,
  repo = "github.com/JuliaSmoothOptimizers/RipQP.jl.git",
  target = "build",
  devbranch = "master"
)


text:
###
  Name     : c7_21_JuMP02.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 3/15/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using JuMP
using Gurobi
#master = Model()
#master = Model(solver = GLPKSolverLP(method=:Exact))
master = Model(solver =GurobiSolver())
@variable(master, x[1:2] >= 0)
w=[14 31 36 45]
A=[1 0; 1 0; 0 2; 1 0]
b=[211; 395; 610; 97]
@defConstrRef myCons[1:4]
for i=1:4
    myCons[i] = @addConstraint(master, dot(x, vec(A[i,:]))>=b[i])
end
@setObjective(master, Min, sum(x))
master
status=solve(master)
getValue(x)
println("\nOptimal Solution is:\n")s
println("width: ", w)
epsilon=1e-6
for i=1:size(A,2)
    if getValue(x[i])>epsilon
        println("Cutting Pattern: ", A[:,i], ", Number of Paper Rolls Cut Using this Pattern: ", getValue(x[i]))
    end
end

# http://www.juliaopt.org/notebooks/Chiwei%20Yan%20-%20Cutting%20Stock.html


text:
using Hecke
using Base.Test

#include("NfMaxOrd-test.jl")
#include("NfMaxOrd.jl")
include("NfOrd.jl")

#test_pseudohnf()

set_verbose_level(:UnitGroup, 1)
set_verbose_level(:ClassGroup, 0)

test_NfOrd()
test_NfOrd_Elem()
test_NfOrd_Idl()
test_NfOrd_FracIdl()

Qx, x = PolynomialRing(QQ, "x")

K, a = NumberField(x^2 - 2)
print("Testing $(x^2 - 2)\n")
O = maximal_order(K)
#
# x^5 + 514944*x^2 + 123904 test prime decomposition with this (2 is index divisor and only one prime ideal over 2)
#
print("Tentative class and unit group computation ... \n")
c, U, b = Hecke._class_unit_group(O);
print("Saturating the tentative unit group ... \n")
Hecke._refine_with_saturation(c, U)

K, a = NumberField(x^3 - 2)
print("Testing $(x^3 - 2)\n")
O = maximal_order(K)
#
# x^5 + 514944*x^2 + 123904 test prime decomposition with this (2 is index divisor and only one prime ideal over 2)
#
print("Tentative class and unit group computation ... \n")
c, U, b = Hecke._class_unit_group(O);
print("Saturating the tentative unit group ... \n")
Hecke._refine_with_saturation(c, U)

K, a = NumberField(x^5 - 11^2 * 7)
print("Testing $(x^5 - 11^2 * 7)\n")
O = maximal_order(K)
#
# x^5 + 514944*x^2 + 123904 test prime decomposition with this (2 is index divisor and only one prime ideal over 2)
@test discriminant(O) == 109853253125
#
print("Tentative class and unit group computation ... \n")
c, U, b = Hecke._class_unit_group(O);
print("Saturating the tentative unit group ... \n")
Hecke._refine_with_saturation(c, U)
##
f = x^18 + 18*x^16 + 135*x^14 + 192*x^12 - 2961*x^10 - 17334*x^8+ 20361*x^6 +  315108*x^4 + 514944*x^2 + 123904
#
print("Testing $f\n")
println(signature(f))
K, a = NumberField(f)
O = maximal_order(K)
#
print("Tentative class and unit group computation ... \n")
c, U, b = Hecke._class_unit_group(O);
print("Saturating the tentative unit group ... \n")
Hecke._refine_with_saturation(c, U)
#
@test c.h == 36
@test Hecke._validate_class_unit_group(c, U) == 1
##
#f = x^5 + 514944*x^2 + 123904
##
#print("Testing $f\n")
#println(signature(f))
#K, a = NumberField(f)
#O = maximal_order(K)
##
#print("Tentative class and unit group computation ... \n")
#c, U, b = Hecke._class_unit_group(O);
#print("Saturating the tentative unit group ... \n")
#Hecke._refine_with_saturation(c, U)
#
#@test c.h == 2


text:
export
    isregular, mirror, upsample,
    downsample, coefficientcount, partition,
    combine, copyrange, copyrange!

#####
##### Collection utilities
#####

"""
    isregular(signal::AbstractArray)

Return `true` if all dimensions of a `signal` are of equal size.
"""
isregular(signal::AbstractArray) = isone(length(unique(size(signal))))

"""
    mirror(signal::AbstractVector)

Flip the sign of each odd index in `signal` and return the resulting vector.
Roughly ten times faster than the original implementation.
"""
mirror(signal::AbstractVector) = [flipsign(sample, -iseven(index)) for (index, sample) in enumerate(signal)]

"""
    mirror(signal::StaticVector)

Flip the sign of each odd index in `signal` and return the resulting vector.
About 1/3 faster than when dispatching on `Base` arrays and non-allocating.
"""
mirror(signal::StaticVector) = typeof(signal)(_mirror(signal))

@generated function _mirror(signal::StaticVector{N,T}) where {N,T}
    return Expr(:tuple, (iseven(i) ? :(-signal[$i]) : :(signal[$i]) for i in 1:N)...)
end

"""
    upsample(signal; skipfirst=true)

Doubles the number of samples in `signal` by inserting a zero of `eltype(signal)`
on one side of each sample. If `skipfirst` is `true`, these zeros are inserted on
the left side of each element; otherwise, the zeros are inserted on the right side.
"""
function upsample(signal; skipfirst=false)
    dest = zeros(eltype(signal), length(signal) << 1)
    _upsample!(dest, signal; skipfirst)
    return dest
end

_upsample!(dest, src; skipfirst=false) = copyto!(@view(dest[(begin + skipfirst):2:end]), src)

"""
    downsample(signal; skipfirst=true)

Selects every other sample in `signal` and returns those samples as a vector.
If `skipfirst` is `true`, the sample selection process begins with the first
sample; otherwise, the process begins with the second sample in the signal.
"""
downsample(signal; skipfirst=false) = _downsample!(similar(signal, length(signal) >> 1), signal; skipfirst)

_downsample!(dest, src; skipfirst=false) = copyto!(dest, @view(src[(begin + skipfirst):2:end]))

"""
    coefficientcount(signal; threshold=zero(eltype(signal)), level::Integer=1)

Returns the number of coefficients in `signal` greater than or equal to
`threshold`, excluding coefficients in levels below `level`.
"""
function coefficientcount(signal; threshold=zero(eltype(signal)), level::Integer=1)
    level > 0 || throw(ArgumentError("transform level must be greater than zero"))
    start = level > 1 ? 1 + (1 << (level - 2)) : 1
    return count(>=(threshold) ∘ abs, @view(signal[start:end]))
end

"""
    partition(signal; upto=lastindex(signal))

Split the odd- and even-indexed elements in `signal` from the start
to the index `upto`, storing the odd and even elements in-place, respectively,
then the remaining elements from `signal`, returning the result.
"""
partition(signal; upto=lastindex(signal)) = _partition!(copy(signal), signal; upto)

function _partition!(dest, src; upto=lastindex(src))
    odd_range, even_range = 1:2:upto, 2:2:upto
    odd_size, even_size = length(odd_range), length(even_range)
    odd_dest, even_dest = 1:odd_size, (odd_size + 1):upto
    @views begin
        copyto!(dest[odd_dest], src[odd_range])
        copyto!(dest[even_dest], src[even_range])
        if upto != lastindex(src)
            remainder = (upto + 1):lastindex(src)
            copyto!(dest[remainder], src[remainder])
        end
    end
    return dest
end

"""
    partition(signal, range::AbstractRange)

Split the odd- and even-indexed elements in `signal` in the range
given by `range`, storing the odd and even elements in-place, respectively,
then the remaining elements from `signal`, returning the result.
"""
partition(signal, range::AbstractRange) = copyto!(copy(signal), partition(view(signal, range)))

"""
    combine(signal; upto=lastindex(signal))

Return the inverse operation of `partition(signal; upto=lastindex(signal))`.
"""
combine(signal; upto=lastindex(signal)) = _combine!(copy(signal), signal; upto)

function _combine!(dest, src; upto=lastindex(src))
    pivot = ceil(Int, upto / 2)
    @views begin
        copyto!(dest[1:2:upto], src[begin:pivot])
        copyto!(dest[2:2:upto], src[(pivot + 1):upto])
        copyto!(dest[upto:end], src[upto:end])
    end
    return dest
end

"""
    combine(signal; range::AbstractRange)

Return the inverse operation of `partition(signal, range::AbstractRange)`.
"""
function combine(signal, range::AbstractRange)
    dest = copy(signal)
    src = combine(signal; upto=length(range))
    @views copyto!(dest[range], src[begin:length(range)])
    return dest
end

"""
    copyrange!(dest, src, range)

Copy a `range` of elements from `src` into `dest` and return the resulting vector.
Adapted from `Wavelets.jl`'s `stridedcopy!`.
"""
copyrange!(dest, src, range) = copyto!(dest, view(src, range))

# Notable omission: `circshift!`, which now has an implementation in Base Julia.


text:
module DOM

using Lazy

# Basic Types

type Node{tag}
  attrs::Dict{Any, Any}
  children::Vector{Any}
end

tag{T}(node::Node{T}) = T
attrs(node::Node) = node.attrs
children(node::Node) = node.children

Node(tag::Symbol, attrs::Dict = Dict(), content::AbstractVector = c(); kws...) =
  Node{tag}(isempty(kws) ? attrs : merge(attrs, Dict([kws...])),
            content)

Node(tag::Symbol, attrs::Dict, content; kws...) =
  Node(tag, attrs, c(content); kws...)

Node(tag::Symbol, content; kws...) =
  Node(tag, Dict(), content; kws...)

Node(tag::Symbol, content::Node...) = Node(tag, collect(content))

# CSS selector parsing

function cssparse(s)
  attrs = Dict()
  id = match(r"#[A-Za-z0-9]+", s)
  id == nothing || (attrs[:id] = id.match[2:end])
  classes = matchall(r"\.[A-Za-z0-9]+", s)
  isempty(classes) || (attrs[:class] = map(s->s[2:end], classes))
  return attrs
end

Node(tag::Symbol, selector::AbstractString, props::Dict, args...; kws...) =
  Node(tag, merge!(cssparse(selector), props), args...; kws...)

Node(tag::Symbol, selector::AbstractString, content, args...; kws...) =
  Node(tag, cssparse(selector), content, args...; kws...)

# Rendering

export htmlescape

attrstring(xs::Vector) = join(xs, " ")
attrstring(x) = string(x)
attrstring(d::Dict) = @as _ d map(t->"$(t[1])=\"$(attrstring(t[2]))\"", _) join(_, " ")

htmlescape(s::AbstractString) =
    @> s replace(r"&(?!(\w+|\#\d+);)", "&") replace("<", "<") replace(">", ">") replace("\"", """)

render(io::IO, s::AbstractString) = print(io, htmlescape(s))

function render(io::IO, node::Node)
  print(io, "<", tag(node))
  isempty(attrs(node)) || print(io, " ", attrstring(attrs(node)))
  print(io, ">")
  render(io, children(node))

  print(io, "</", tag(node), ">")
end

function render(io::IO, xs::Vector)
  for x in xs
    render(io, x)
  end
end

render(io::IO, x) = writemime(io, MIME"text/html"(), x)

Base.writemime(io::IO, ::MIME"text/html", node::Node) = render(io, node)

Base.show(io::IO, node::Node) = render(io, node)

Node(tag::Symbol, io::IO, args...) = render(io, Node(tag, args...))

# Specific elements

for t in :[div span a h1 h2 h3 b table tr td em strong].args
  eval(Expr(:export, t))
  @eval ($t)(args...; kws...) = Node($(Expr(:quote, t)), args...; kws...)
end

end


text:
using PWS
using Plots

system = PWS.chemotaxis_system(
    mean_L = 20,
    num_receptors = 10_000,
    Y_tot = 5_000,
    LR_timescale = 1e-2,
    Y_timescale = 1e-1,
    Y_ratio = 1/8,
    LR_ratio = 0.5,
    q=0,
    dtimes = 0:0.04:2
)

@time sol = PWS._solve(system)

p = plot(sol)

@time conf = PWS.generate_configuration(system)

cond_ens = PWS.ConditionalEnsemble(system)

alg = PWS.SMCEstimate(256)

@profview result = PWS.simulate(alg, conf, cond_ens).samples

result

using Distributions

for x in eachrow(result)
    @show mean(x) / var(x)
end

θ = 11.5
p = plot(xlabel="Path Action", ylabel="Probability")
fitted_dists = []
for (c, i) in enumerate(11:10:51)
    x = result[i,:]
    t = dtimes[i]
    x = (mean(x) + var(x) / θ) .- x
    plot!(p, x, seriestype=:stephist, xlim=(-10,400), fill=true, normalize=true, fillalpha=0.1, label="T=$t", color=c)

    fitted_gamma = @show fit(Gamma{Float64}, x)
    push!(fitted_dists, fitted_gamma)
    plot!(p, 0:1:400, y->pdf(fitted_gamma, y), color=c, label="Gamma distribution fit")
end
p
savefig(plot(p, dpi=300), "~/Downloads/histograms.png")

p = plot(result[end,:], seriestype=:stephist, fill=true, legend=false, yscale=:log10, ylim=(1e-6,1e-2), normalize=true)
savefig(p, "~/Downloads/histogram.png")

using Statistics

function bootstrap(f, data::AbstractVector; bsamples=100)
    N = length(data)
    map(1:bsamples) do i
        sample_indices = rand(1:N, N)
        f(@view data[sample_indices])
    end
end

PWS.logmeanexp(result[end, :])

bdata = bootstrap(PWS.logmeanexp, result[end,:], bsamples=10^4)
h = histogram(bdata, title="Bootstrapping Histogram", xlabel="Path Action", ylabel="Count")
savefig(h, "~/Downloads/bootstrap.png")

using Distributions

p=plot()
for (i, shape_param) in enumerate([2, 5, 10, 20])
    dist = Gamma(shape_param, 1.0)
    plot!(p, 0.0:0.01:30, x->pdf(dist, x), color=i, linestyle=:solid, label="Gamma α=$shape_param")
    plot!(p, 0.0:0.01:30, x->pdf(dist, x)*exp(-x)*(scale(dist) + 1)^shape(dist), color=i, linestyle=:dot, fill=true, fillalpha=0.1, label="exp(-x)*Gamma(x)")
end
p
savefig(plot(p, dpi=300), "~/Downloads/gamma.png")

p=plot()
for (i, dist) in enumerate(fitted_dists)
    plot!(p, 0.0:1:400, x->pdf(dist, x), color=i, linestyle=:solid)
    plot!(p, 0.0:1:400, x->pdf(dist, x)*exp(-x)*(scale(dist) + 1)^shape(dist), color=i, linestyle=:dot, fill=true, fillalpha=0.1)
end
p


text:
module bfgaFibonacci

    import Base.isless

    include("Types.jl")
    using .Types

    #using GeneticAlgorithms
    #include("GeneticAlgorithms.jl")
    #using .GeneticAlgorithms

    include("bfga.jl")
    using .bfga

    include("Bf.jl")
    using .BfInterpreter

    export _trainingExamples, fitness, nb_predictions

    _trainingExamples = [ UInt8[1,2] , UInt8[3,5] , UInt8[0,1] , UInt8[8,13] ]
    nb_predictions = 2

    function fitness(ent, instructionsSet)
        #println(" $(join(genesToBfInstructions(ent.dna),"")) ")
        score = 0
        @sync for i in 1:4
            @async score += fitness_aux(ent, i, instructionsSet)
        end
        ent.fitness = score
        score
    end

    function fitness_aux(ent, instructionsSet)
        n = rand(1:4)
        fitness_aux(ent, n, instructionsSet)
    end

    function fitness_aux(ent, num, instructionsSet)
        input = _trainingExamples[num]
        goal =  Int(sum(input))
        last_value = Int(input[2])
        score =0
        #target_length = length(goal)
        #target_score = target_length*256
        try
            output, m_Ticks = execute(ent.program, input, instructionsSet)
            n= length(output)

            #score = 256 - abs(Int(output[1]) - goal)
            for i in 1:min(n,nb_predictions)
                #@show Int(output[i]), goal, last_value
                #@show Int(256 - abs(Int(output[i]) - goal))
                score += 256 - abs(Int(output[i]) - goal)
                goal, last_value = last_value + goal, goal
            end

            ent.bonus += (2000 - m_Ticks)
            abs(score) # - target_score)
        catch y
            0
        end
    end

    function getBfCode(ent)
        join( ent.program , "")
    end

    function simulate_entity(ent, instructionsSet)
        #bft = bfType(ent.dna)
        _res = "code : $(ent.program ) "
        #println(_res)
        input = _trainingExamples[rand(1:4)]
        goal =  sum(input)
        last_value = input[2]
        mem = []
        for i in 1:nb_predictions
            push!(mem, goal)
            goal, last_value = last_value + goal, goal
        end

        try
            output, _ = execute(ent.program, input, instructionsSet)
            n= length(output)

            if length(output) > nb_predictions
                _res = _res* "\n $input -> $(map(e -> Int(e),output[1:nb_predictions])) "
            else
                _res = _res* "\n $input -> $(map(e -> Int(e),output)) "
            end
        catch y
            _res = _res* "\n BEST raises Errors \n Error : $y "

        end

        return _res
    end

    function getTargetFitness()
        256 * 4 * nb_predictions
    end

    function getParams()
        #=
        populationSize ::  Int
        generations ::  Int
        genomeSize ::  Int
        maxGenomeSize :: Int
        crossoverRate :: Float64
        mutationRate :: Float64
        elitism  :: Bool # Keep previous generation's fittest individual in place of worst in current
        historyPath :: String  # Path to save log of fitness history at each generation. Can be used to plot on Excel chart etc.

        totalFitness :: Float64
        targetFitness :: Float64
        targetFitnessCount ::  Int
        currentGeneration ::  Int

        #thisGeneration :: Array
        #nextGeneration :: Array
        #fitnessTable :: Array{Float64,1}
        =#
        logfile = open("../Results/logFibonacci.log", "w")

        tgFitness =  getTargetFitness()
        println("targetFitness = $tgFitness ")
        write(logfile, "targetFitness = $tgFitness \n")
        return Main.GeneticAlgorithms.Types.GAParams(136, 1000000 , 120, 150, 0.7, 0.01, true, logfile ,  0.0 , tgFitness, 0.0 , 0 )
    end

end

include("../src/GeneticAlgorithms.jl")
using .GeneticAlgorithms
include("bfga.jl")
using .bfga

using Distributed
using Pkg

function test_serial()
    GeneticAlgorithms.runga(bfga, bfgaFibonacci) #, initial_pop_size = 136)
end


text:
###################################################################################################################
#Auxiliar functions
###################################################################################################################

function gridShape(dims::Int,box::Array{<:Real,2},radius::Union{<:Real,Array{<:Real,1}})

    #Check inconsistencies in shape
    if (dims,2) != size(box)
        error("Community.simulationBox has to be of dimensions of (dims,2), where dims are the dimensions of the model.")
    end

    if typeof(radius) <: Array{<:Real,1} && length(radius) != dims
        error("Community.radiusInteraction has to be a Real Number or an array with the dimensions of the model.")
    end

    if any([box[i,2]<= box[i,1] for i in 1:dims]) && dims > 0
        error("Community.simulationBox must have all the minimum values bigger than the maximum.")
    end

    if dims == 0
        error("For models with dimension 0, it does not make sense to use grid neighbours computation. Compile with other method instead.")
    end

    #Compute additional parameters for the structure
    if dims > 0
        if typeof(radius) <: Array{<:Real,1}
            axisSize = [ceil(INT,(box[i,2]-box[i,1])/radius[i])+2 for i in 1:size(box)[1]]
        elseif typeof(radius) <: Real
            axisSize = [ceil(INT,(box[i,2]-box[i,1])/radius)+2 for i in 1:size(box)[1]]
        end
        # cum = cumprod(axisSize)
        # cumSize = [1;cum[1:end-1]]
        # n = cum[end]
    else
        axisSize = []
    end

    return axisSize
end

function position2gridVectorPosition_(x,boxX,dX,nX)

    pos = min(max(ceil(Int,(x-boxX)/dX+1),1),nX)

    return pos
end

function position2gridVectorPosition_(x,boxX,dX,nX,y,boxY,dY,nY)

    posX = min(max(ceil(Int,(x-boxX)/dX+1),1),nX)
    posY = min(max(ceil(Int,(y-boxY)/dY+1),1),nY)

    return posX + nX * (posY-1)
end

function position2gridVectorPosition_(x,boxX,dX,nX,y,boxY,dY,nY,z,boxZ,dZ,nZ)

    posX = min(max(ceil(Int,(x-boxX)/dX+1),1),nX)
    posY = min(max(ceil(Int,(y-boxY)/dY+1),1),nY)
    posZ = min(max(ceil(Int,(z-boxZ)/dZ+1),1),nZ)

    return posX + nX * (posY-1) + nX * nY * (posZ-1)
end

function gridVectorPositionNeighbour_(p,neighbour,nX,periodicX)
    posX = p+(neighbour-2)

    if periodicX == true
        if posX == 1
            return nX - 1
        elseif posX == nX
            return 2
        else
            return posX
        end
    else
        if posX == 0
            return -1
        elseif posX == nX + 1
            return -1
        else
            return posX
        end
    end
end

function gridVectorPositionNeighbour_(p,neighbour,nX,periodicX,nY,periodicY)

    posY = floor(Int,p/nX-1/nX)
    p -= nX*posY
    posX = p

    addY = floor(Int,neighbour/3-1/3)
    neighbour -= addY*3
    addX = neighbour - 2

    posY += addY
    posX += addX

    abort = false

    if periodicX == true
        if posX == 1
            posX =  nX - 1
        elseif posX == nX
            posX = 2
        end
    else
        if posX == 0
            abort = true
        elseif posX == nX + 1
            abort = true
        end
    end

    if periodicY == true
        if posY == 1
            posY = nY - 1
        elseif posY == nY
            posY = 2
        end
    else
        if posY == 0
            abort= true
        elseif posY == nY + 1
            abort = true
        end
    end

    if abort
        return -1
    else
        return posX + nX*(posY-1)
    end
end

function gridVectorPositionNeighbour_(p,neighbour,nX,periodicX,nY,periodicY,nZ,periodicZ)

    nXY = nX*nY
    posZ = floor(Int,p/nXY-1/nXY)
    p -= nXY*posZ
    posY = floor(Int,p/nX-1/nX)
    p -= nX*posY
    posX = p

    addZ = floor(Int,neighbour/9-1/9)
    neighbour -= addZ*9
    addY = floor(Int,neighbour/3-1/3)
    neighbour -= addY*3
    addX = neighbour - 2

    posZ += addZ
    posY += addY
    posX += addX

    abort = false

    if periodicX == true
        if posX == 1
            posX =  nX - 1
        elseif posX == nX
            posX = 2
        end
    else
        if posX == 0
            abort = true
        elseif posX == nX + 1
            abort = true
        end
    end

    if periodicY == true
        if posY == 1
            posY = nY - 1
        elseif posY == nY
            posY = 2
        end
    else
        if posY == 0
            abort= true
        elseif posY == nY + 1
            abort = true
        end
    end

    if periodicZ == true
        if posZ == 1
            posZ = nZ - 1
        elseif posZ == nZ
            posZ = 2
        end
    else
        if posZ == 0
            abort= true
        elseif posZ == nZ + 1
            abort = true
        end
    end

    if abort
        return -1
    else
        return posX + nX*(posY-1) + nXY*(posZ-1)
    end
end

#########################################################
#cpu#####################################################
#########################################################
function gridInsertCounts1D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)

    lockadd_ = Threads.SpinLock()
    Threads.@threads for ic1_ = 1:N

        nnPosIdCell_[ic1_] = AgentBasedModels.position2gridVectorPosition_(localV[ic1_,1],box[1,1],r[1],n[1])
        lock(lockadd_)
            nnGridCounts_[nnPosIdCell_[ic1_]]+=1
        unlock(lockadd_)

    end

    return nothing

end

function gridInsertCounts2D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)

    lockadd_ = Threads.SpinLock()
    Threads.@threads for ic1_ = 1:N

        nnPosIdCell_[ic1_] = AgentBasedModels.position2gridVectorPosition_(localV[ic1_,1],box[1,1],r[1],n[1],
                                                                            localV[ic1_,2],box[2,1],r[2],n[2])
        lock(lockadd_)
            nnGridCounts_[nnPosIdCell_[ic1_]]+=1
        unlock(lockadd_)

    end

    return nothing

end

function gridInsertCounts3D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)

    lockadd_ = Threads.SpinLock()
    Threads.@threads for ic1_ = 1:N

        nnPosIdCell_[ic1_] = AgentBasedModels.position2gridVectorPosition_(localV[ic1_,1],box[1,1],r[1],n[1],
                                                                            localV[ic1_,2],box[2,1],r[2],n[2],
                                                                            localV[ic1_,3],box[3,1],r[3],n[3])
        lock(lockadd_)
            nnGridCounts_[nnPosIdCell_[ic1_]]+=1
        unlock(lockadd_)

    end

    return nothing

end

function gridCountingSort_cpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    lockadd_ = Threads.SpinLock()
    Threads.@threads for ic1_ = 1:N

        id_ = nnPosIdCell_[ic1_]
        if id_ == 1
            posInit_ = 0
        else
            posInit_ = nnGridCountsCum_[id_-1]
        end
        lock(lockadd_)
            posCell_ = nnGridCountsAux_[id_]
            nnGridCountsAux_[id_]+=1
        unlock(lockadd_)
        nnCellIdPos_[posInit_+posCell_] = ic1_

    end

    return nothing

end

function gridComputeNN1D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    nnGridCounts_ .= 0
    nnGridCountsAux_ .= 1
    gridInsertCounts1D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)
    nnGridCountsCum_ .= cumsum(nnGridCounts_)
    gridCountingSort_cpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    return nothing
end

function gridComputeNN2D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    nnGridCounts_ .= 0
    nnGridCountsAux_ .= 1
    gridInsertCounts2D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)
    nnGridCountsCum_ .= cumsum(nnGridCounts_)
    gridCountingSort_cpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    return nothing
end

function gridComputeNN3D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    nnGridCounts_ .= 0
    nnGridCountsAux_ .= 1
    gridInsertCounts3D_cpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)
    nnGridCountsCum_ .= cumsum(nnGridCounts_)
    gridCountingSort_cpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    return nothing
end

#########################################################
#gpu#####################################################
#########################################################
function gridInsertCounts1D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:N
        nnPosIdCell_[ic1_] = AgentBasedModels.position2gridVectorPosition_(localV[ic1_,1],box[1,1],r[1],n[1])
        CUDA.atomic_add!(CUDA.pointer(nnGridCounts_,nnPosIdCell_[ic1_]),Int32(1))
    end

    return nothing

end

function gridInsertCounts2D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:N
        nnPosIdCell_[ic1_] = AgentBasedModels.position2gridVectorPosition_(localV[ic1_,1],box[1,1],r[1],n[1],
                                                                            localV[ic1_,2],box[2,1],r[2],n[2])
        CUDA.atomic_add!(CUDA.pointer(nnGridCounts_,nnPosIdCell_[ic1_]),Int32(1))
    end

    return nothing

end

function gridInsertCounts3D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:N
        nnPosIdCell_[ic1_] = AgentBasedModels.position2gridVectorPosition_(localV[ic1_,1],box[1,1],r[1],n[1],
                                                                            localV[ic1_,2],box[2,1],r[2],n[2],
                                                                            localV[ic1_,3],box[3,1],r[3],n[3])
        CUDA.atomic_add!(CUDA.pointer(nnGridCounts_,nnPosIdCell_[ic1_]),Int32(1))
    end

    return nothing

end

function gridCountingSort_gpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
    stride_ = blockDim().x * gridDim().x

    for ic1_ in index_:stride_:N
        id_ = nnPosIdCell_[ic1_]
        if id_ == 1
            posInit_ = 0
        else
            posInit_ = nnGridCountsCum_[id_-1]
        end
        posCell_ = CUDA.atomic_add!(CUDA.pointer(nnGridCountsAux_,id_),Int32(1))
        nnCellIdPos_[posInit_+posCell_] = ic1_
    end

    return nothing

end

function gridComputeNN1D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    nnGridCounts_ .= 0
    nnGridCountsAux_ .= 1
    #Insert counts
    kernel_ = @cuda launch = false gridInsertCounts1D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_);
    prop_ = AgentBasedModels.configurator_(kernel_,N);
    kernel_(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_;threads=prop_[1],blocks=prop_[2])
    #Cumsum
    nnGridCountsCum_ .= cumsum(nnGridCounts_)
    #Insert counts Sort
    kernel_ = @cuda launch = false gridCountingSort_gpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_);
    prop_ = AgentBasedModels.configurator_(kernel_,N);
    kernel_(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_;threads=prop_[1],blocks=prop_[2])

    return nothing
end

function gridComputeNN2D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    nnGridCounts_ .= 0
    nnGridCountsAux_ .= 1
    #Insert counts
    kernel_ = @cuda launch = false gridInsertCounts2D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_);
    prop_ = AgentBasedModels.configurator_(kernel_,N);
    kernel_(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_;threads=prop_[1],blocks=prop_[2])
    #Cumsum
    nnGridCountsCum_ .= cumsum(nnGridCounts_)
    #Insert counts Sort
    kernel_ = @cuda launch = false gridCountingSort_gpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_);
    prop_ = AgentBasedModels.configurator_(kernel_,N);
    kernel_(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_;threads=prop_[1],blocks=prop_[2])

    return nothing
end

function gridComputeNN3D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)

    nnGridCounts_ .= 0
    nnGridCountsAux_ .= 1
    #Insert counts
    kernel_ = @cuda launch = false gridInsertCounts3D_gpu!(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_);
    prop_ = AgentBasedModels.configurator_(kernel_,N);
    kernel_(N,localV,box,r,n,nnPosIdCell_,nnGridCounts_;threads=prop_[1],blocks=prop_[2])
    #Cumsum
    nnGridCountsCum_ .= cumsum(nnGridCounts_)
    #Insert counts Sort
    kernel_ = @cuda launch = false gridCountingSort_gpu!(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_);
    prop_ = AgentBasedModels.configurator_(kernel_,N);
    kernel_(N,nnPosIdCell_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_;threads=prop_[1],blocks=prop_[2])

    return nothing
end

###################################################################################################################
#Call functions
###################################################################################################################

function argumentsGrid_!(program::Program_, platform::String)

    append!(program.declareVar.args,
    (quote
        n_ = AgentBasedModels.gridShape($(program.agent.dims),com.simulationBox,com.radiusInteraction)
        nnPosIdCell_ = zeros($INT,nMax)
        nnGridCounts_ = zeros($INT,prod(n_))
        nnGridCountsAux_ = zeros($INT,prod(n_))
        nnGridCountsCum_ = zeros($INT,prod(n_))
        nnCellIdPos_ = zeros($INT,nMax)
        if typeof(com.radiusInteraction) <: Real
            radiusInteraction = com.radiusInteraction.*ones($INT,$(program.agent.dims))
        else
            radiusInteraction = Array(com.radiusInteraction)
        end
        n_ = Array(n_)
    end).args
    )

    #Create kernels for the algorithm
    if platform != "cpu" && platform != "gpu"
        error("Platform should be or cpu or gpu. ", platform, " was given.")
    end

    append!(program.args, [:nnPosIdCell_,:nnGridCounts_,:nnGridCountsAux_,:nnGridCountsCum_,:nnCellIdPos_,:n_])

    aux = Meta.parse(string("AgentBasedModels.gridComputeNN",program.agent.dims,"D_",platform,"!(N,localV,simulationBox,radiusInteraction,n_,nnPosIdCell_,nnGridCounts_,nnGridCountsCum_,nnGridCountsAux_,nnCellIdPos_)"))

    push!(program.execInit.args, :($aux))
    push!(program.execInloop.args, :($aux))

    return nothing
end

function loopGrid_(program::Program_, code::Expr, platform::String)

    code = vectorize_(program.agent, code, program, interaction=true)

    #Compute the cell id for each agent in grid and linearized format for any box possible box dimensions.
    if program.agent.dims == 1
        periodicX = typeof(program.agent.boundary.boundaries[1]) == Periodic
        loop=:(begin
            for iAux1_ in 1:3
                pos_ = AgentBasedModels.gridVectorPositionNeighbour_(nnPosIdCell_[ic1_],iAux1_,n_[1],$periodicX)
                if pos_ != -1
                    max_ = nnGridCountsCum_[pos_]
                    box_ = max_ - nnGridCounts_[pos_] + 1
                    for ic2_ in box_:max_
                        ic2_ = nnCellIdPos_[ic2_]
                        $code
                    end
                end
            end
        end)
    elseif program.agent.dims == 2
        periodicX = typeof(program.agent.boundary.boundaries[1]) == Periodic
        periodicY = typeof(program.agent.boundary.boundaries[2]) == Periodic
        loop=:(begin
            for iAux1_ in 1:9
                pos_ = AgentBasedModels.gridVectorPositionNeighbour_(nnPosIdCell_[ic1_],iAux1_,n_[1],$periodicX,n_[2],$periodicY)
                if pos_ != -1
                    max_ = nnGridCountsCum_[pos_]
                    box_ = max_ - nnGridCounts_[pos_] + 1
                    for ic2_ in box_:max_
                        ic2_ = nnCellIdPos_[ic2_]
                        $code
                    end
                end
            end
        end)
    elseif program.agent.dims == 3
        periodicX = typeof(program.agent.boundary.boundaries[1]) == Periodic
        periodicY = typeof(program.agent.boundary.boundaries[2]) == Periodic
        periodicZ = typeof(program.agent.boundary.boundaries[3]) == Periodic
        loop=:(begin
            for iAux1_ in 1:27
                pos_ = AgentBasedModels.gridVectorPositionNeighbour_(nnPosIdCell_[ic1_],iAux1_,n_[1],$periodicX,n_[2],$periodicY,n_[3],$periodicZ)
                if pos_ != -1
                    max_ = nnGridCountsCum_[pos_]
                    box_ = max_ - nnGridCounts_[pos_] + 1
                    for ic2_ in box_:max_
                        ic2_ = nnCellIdPos_[ic2_]
                        $code
                    end
                end
            end
        end)
    end

    loop = subs_(loop,:nnic2_,:(ic2_))
#    loop = postwalk(x->@capture(x,nnic2_) ? :(nnId_[ic2_]) : x, loop)
    loop = simpleFirstLoop_(platform, loop)

    return loop
end


text:
module Solvers
    using Tsp.Model
    include("nn.jl")
    include("Tabu.jl")
    include("Anneal.jl")
    include("Genetic.jl")
    @enum SolverType tabu=1 annealing=2 genetic=3

    export solver, nearest_neighbour

    for s in instances(SolverType)
        @eval export $(Symbol(s))
    end

    """
    Obtain a solver distance
    """
    function solver(kind::SolverType)::Function
        if kind == tabu
            return Tabu.find_optimal
        elseif kind == annealing
            return Annealing.find_optimal
        end
    end

end


text:
equalto(x) = y -> y == x

array = rand(Bool, (300, 400, 500))
indices = CartesianIndices(array)

for fn in [equalto(true), equalto(false)]
    tracker = EulerTracker(array, fn)

    for i in 1:2000
        idx = rand(indices)
        tracker[idx] = true - tracker[idx]
    end

    @test euler_characteristic(tracker) == BitArray(tracker) .|> fn |> euler_characteristic
end


text:
module BioEnergeticFoodWebs

using Distributions
using OrdinaryDiffEq, DiffEqCallbacks
using JSON
using JLD2
using StatsBase
using Statistics
using LinearAlgebra

export trophic_rank,
  model_parameters,
  #rewire_parameters,
  simulate,
  nichemodel,
  save,
  population_stability,
  total_biomass,
  population_biomass,
  foodweb_evenness,
  species_richness,
  species_persistence,
  producer_growth,
  nutrient_intake,
  consumer_intake,
  metabolism,
  adbm_model,
  temperature_size_rule,
  NoEffectTemperature,
  ExtendedEppley,
  ExponentialBA,
  ExtendedBA,
  Gaussian

# Includes
include(joinpath(".", "trophic_rank.jl"))
include(joinpath(".", "checks.jl"))
include(joinpath(".", "dBdt.jl"))
include(joinpath(".", "make_parameters.jl"))
include(joinpath(".", "simulate.jl"))
include(joinpath(".", "random.jl"))
include(joinpath(".", "measures.jl"))

include(joinpath(".", "rewiring/ADBM.jl"))
include(joinpath(".", "rewiring/GilljamRewire.jl"))
include(joinpath(".", "rewiring/StaniczenkoRewire.jl"))
include(joinpath(".", "rewiring/parameters/checkParameters.jl"))
include(joinpath(".", "rewiring/parameters/updateParameters.jl"))

include(joinpath(".", "temperature_dependence_functions.jl"))
include(joinpath(".", "temperature_size_rule_function.jl"))
include(joinpath(".", "biological_rates.jl"))

end


text:
using SparseArrays
Base.getproperty(S::SparseMatrixCSC, ::Symbol) = error("use accessor function")
Base.getproperty(S::SparseVector, ::Symbol) = error("use accessor function")

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: milesfrain/julia

max_forks_repo_name: milesfrain/julia


text:
# module Bukdu

# in filter.jl
export plugins, before, after

# in controller.jl
export ApplicationController
export get, post, delete, patch, put # verbs
export index, edit, new, show, create, update, delete # actions

# in layout.jl
export ApplicationLayout, Layout, layout, /

# in octo.jl
export Assoc
export validates

# in view.jl
export ApplicationView, View

# in renderers.jl
export render
export Tag # html

# in logger.jl
export Logger

# in router.jl
export ApplicationRouter, Router, NoRouteError, reset
export scope, resources

# in endpoint.jl
export ApplicationEndpoint, Endpoint, reload

# in plug.jl
export Plug, plug

# in server.jl
export start, stop


text:
#md# The positions of words can be initialized with pre-trained word vectors so that similar words will appear near each other.
#md# ### Words
using WordCloud
stwords = ["us"];
words_weights = processtext(open(pkgdir(WordCloud)*"/res/Barack Obama's First Inaugural Address.txt"), stopwords=WordCloud.stopwords_en ∪ stwords)
words_weights = Dict(zip(words_weights...))
#md# ### Embeddings
using Embeddings
using TSne
const embtable = load_embeddings(GloVe{:en})
const get_word_index = Dict(word=>ii for (ii,word) in enumerate(embtable.vocab))
function get_embedding(word)
    ind = get_word_index[word]
    emb = embtable.embeddings[:,ind]
    return emb
end
wordvec = Dict()
for k in keys(words_weights)
    if k in keys(get_word_index)
        wordvec[k] = get_embedding(k)
    elseif lowercase(k) in keys(get_word_index)
        wordvec[k] = get_embedding(lowercase(k))
    else
        pop!(words_weights, k)
        println("remove ", k)
    end
end
embedded = tsne(hcat(values(wordvec)...)', 2)
#md# ### WordCloud
wc = wordcloud(
    words_weights,
    maskshape = box,
    masksize = (1000, 1000, 0),
    density=0.3,
    run = initwords!
)

pos = embedded
mean = sum(pos, dims=1) / size(pos, 1)
r = maximum(sqrt.(pos[:,1].^2 + pos[:,2].^2 ))
pos = (pos .- mean) ./ 2r
sz = collect(size(wc.mask))'
pos = round.(Int, pos .* sz .+ sz ./ 2)

setpositions!(wc, keys(wordvec)|>collect, eachrow(pos), type=setcenter!)
setstate!(wc, :placewords!)
generate!(wc, teleporting=false)
println("results are saved to embedding.png")
paint(wc, "embedding.png")
wc
#eval# runexample(:embedding)
#md# ![](embedding.png)


text:
function test_crop_array()
    x = rand(16,16,32,8)
    y = crop_array(x,(:,:,30))

    return size(y) == (16,16,30,8)
end

function test_init0()
    x = init0(16,8,4)

    return size(x) == (16,8,4) && sum(x) ≈ 0
end

function test_convertKA()
    x = rand(8,8)
    x = convert2KnetArray(x)

    if CUDA.functional()
        return x isa Knet.KnetArray{Float32}
    else
        return x isa Array{Float32}
    end
end

function test_blowup()
    x = rand(4,4)
    y = blowup_array(x, 3)
    return size(y) == (4,4,3)
end

function test_recycle()
    x = rand(4,4)
    y = recycle_array(x, 4)
    return size(y) == (4,16)
end

function test_de_embed()
    x = rand(8,4,12)
    y = de_embed(x)
    return size(y) == (4,12)
end


text:
module AES
	include("aes-modes.jl")
	export AESEncrypt, AESDecrypt, AESParameters
	export AESECB, AESCBC, AESCFB, AESOFB, AESCTR
end


text:
#!/usr/bin/env julia

pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "ext"))

include("loadmeta.jl")
include("utils.jl")
include("gitmeta.jl")
include("genstdlib.jl")
include("generate.jl")

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
import Statistics
include("quantile.jl")


text:
# File to specify the project exports
# It is not a module as it is just included by the JavaCall module
# to facilitate the exports

export
# Init Options
    defaultopts, fromcurrentvm, forjavahome, setfromcurrentvm!, unsetfromcurrentvm!,
    setjavahome!, unsetjavahome!, pushclasspath!, pushoptions!,
# JNI
    # Types.jl
    jint, jlong, jbyte, jboolean, jchar, jshort, jfloat, jdouble, jsize,
    jvoid, jobject, jclass, jthrowable, jweak, jmethodID, jfieldID, jstring, jarray,
    JNINativeMethod, jobjectArray, jbooleanArray, jbyteArray, jshortArray, jintArray,
    jlongArray, jfloatArray, jdoubleArray, jcharArray, jvalue, jobjectRefType,
    # Constants.jl
    JNI_FALSE, JNI_TRUE,
# Java VM
    init, destroy,
# Jimport
    @jimport


text:
"""
The overall HHPC framework solver object.

Attributes:
    - `graph_tracker::GraphTracker{D,C,AD}`
    - `modal_policies::Dict{D,Policy}` A map from the mode to the modal policies (in and out-horizon)
    - `replan_time_threshold::Int64` The Delta T parameter for periodic replanning (discrete time-steps)
    - `heuristic::Function` The heuristic method for the global layer
    - `curr_state::CMSSPState{D,C}` The current state of the agent, tracked by the solver.
    - `max_steps::Int64` The maximum length of a problem episode (may not be applicable)
    - `start_state::CMSSPState{C,AC}`
    - `goal_modes::Vector{D}`
"""
mutable struct HHPCSolver{D,C,AD,M,B,RNG <: AbstractRNG} <: Solver
    graph_tracker::GraphTracker{D,C,AD,M,B,RNG}
    modal_policies::Dict
    replan_time_threshold::Int64
    goal_modes::Vector{D}
    curr_state::CMSSPState{D,C}
    rng::RNG
    heuristic::Function
    max_steps::Int64
end

function HHPCSolver{D,C,AD,M,B,RNG}(num_samples::Int64, modal_policies::Dict,
                                    deltaT::Int64, goal_modes::Vector{D},
                                    start_state::CMSSPState{D,C},
                                    bookkeeping::B, rng::RNG=Random.GLOBAL_RNG,
                                    heuristic::Function = n->0, max_steps::Int64=1000) where {D,C,AD,M,B,RNG <: AbstractRNG}
    return HHPCSolver(GraphTracker{D,C,AD,M,B,RNG}(num_samples, bookkeeping, rng),
                      modal_policies, deltaT, goal_modes,
                      start_state, rng, heuristic, max_steps)
end

function set_start_state!(solver::HHPCSolver, start_state::CMSSPState{D,C}) where {D,C}
    solver.curr_state = start_state
end

function set_open_loop_samples!(solver::HHPCSolver, num_samples::Int64)
    solver.graph_tracker.num_samples = num_samples
end

"""
Creates the weight function for the edges, used by the global layer graph search.

Arguments:
    - `solver::HHPCSolver{D,C,AD,AC}` The HHPC solver instance
    - `cmssp::CMSSP{D,C,AD,AC}` The CMSSP instance
    - `u::OpenLoopVertex{D,C,AD}` The parent vertex
    - `v::OpenLoopVertex{D,C,AD}` The child vertex
"""
function edge_weight_fn(solver::HHPCSolver, cmssp::CMSSP,
                        u::OpenLoopVertex, v::OpenLoopVertex)

    # u's state to v's pre-state should be in a single region
    @assert u.state.mode == v.pre_bridge_state.mode

    # We care about relative tp_dist between u and v
    # If inf, have modal policy overload it!
    mode = u.state.mode

    # Branch on finite horizon and infinite horizon
    if is_inf_hor(v.tp)

        policy = solver.modal_policies[mode].inf_hor_policy
        pol_value = inf_hor_value(policy, u.state.continuous, v.pre_bridge_state.continuous)

    else
        policy = solver.modal_policies[mode].fin_hor_policy
        mdp = get_mdp(policy)
        temp_curr_time = convert(Int64, round(mean(u.tp)))

        pol_value, weighted_minvalue = horizon_weighted_value(policy,
                                                temp_curr_time,
                                                v.tp,
                                                u.state.continuous,
                                                v.pre_bridge_state.continuous)

        # Filter out edges that nominally would trigger an interrupt
        if pol_value <= mdp.beta_threshold*weighted_minvalue
            return Inf
        end
    end

    edge_weight = -1.0*pol_value

    # Also add the cost due to the assumed mode switch
    if v.pre_bridge_state.mode != v.state.mode
        ac_idx = mode_actionindex(cmssp,v.bridging_action)
        mode_idx = mode_index(cmssp,v.pre_bridge_state.mode)
        edge_weight += -1.0*cmssp.modeswitch_mdp.R[mode_idx,ac_idx]
    end

    return edge_weight
end

# TODO : Fill this in
@POMDP_require solve(solver::HHPCSolver, cmssp::CMSSP) begin

    P = typeof(cmssp)
    S = statetype(P)
    A = actiontype(P)
    D = modetype(cmssp)
    C = continuoustype(cmssp)
    AD = modeactiontype(cmssp)
    AC = controlactiontype(cmssp)
    CS = typeof(cmssp.curr_context_set)
    PR = typeof(cmssp.params)
    M = metadatatype(solver.graph_tracker)
    B = bookkeepingtype(solver.graph_tracker)
    MDPType = Union{ModalFinHorMDP{D,C,AC,PR},ModalInfHorMDP{D,C,AC,PR}}

    @req isterminal(::P, ::S)
    @req generate_sr(::MDPType, ::C, ::AC, ::typeof(solver.rng)) # OR transition + reward?
    @req isterminal(::MDPType, ::C)

    # Global layer requirements
    @req update_vertices_with_context!(::P, ::typeof(solver.graph_tracker), ::Int64)
    @req generate_goal_vertex_set!(::P, ::OpenLoopVertex{D, C, AD, M}, ::typeof(solver.graph_tracker), ::RNG where {RNG <: AbstractRNG})
    @req generate_next_valid_modes(::P, ::D)
    @req generate_bridge_vertex_set!(::P, ::OpenLoopVertex{D,C,AD,M}, ::Tuple{D,D}, ::typeof(solver.graph_tracker), ::AD, ::RNG where {RNG <: AbstractRNG})
    @req update_next_target!(::P, ::typeof(solver), ::Int64)
    @req zero(::Type{M})

    # Local layer requirements
    @req get_relative_state(::MDPType, ::C, ::C)
    @req expected_reward(::MDPType, ::C, ::AC)
    @req convert_s(::Type{Vector{Float64}}, ::C, ::ModalFinHorMDP{D,C,AC,PR})
    @req convert_s(::Type{C}, ::AbstractVector{Float64}, ::ModalFinHorMDP{D,C,AC,PR})

    # HHPC requirements
    @req simulate_cmssp!(::P, ::S, ::Union{Nothing, A}, ::Int64, ::RNG where {RNG <: AbstractRNG})
    @req update_context_set!(::P, ::RNG where {RNG <: AbstractRNG})
    @req get_bridging_action(::OpenLoopVertex{D, C, AD, M})
    @req display_context_future(::CS, ::Int64)

end

# TODO - Update context before or after simulate step???
"""
Executes the top-level behavior of HHPC. Utilizes both global and local layer logic, and interleaving.
"""
function POMDPs.solve(solver::HHPCSolver, cmssp::CMSSP{D,C,AD,AC,P,CS}) where {D,C,AD,AC,P,CS}

    @warn_requirements solve(solver,cmssp)

    # Initialize Global flags
    t = 0
    plan = true
    last_plan_time = 0

    # Defs needed for open-loop-plan
    edge_weight(u::OpenLoopVertex, v::OpenLoopVertex) = edge_weight_fn(solver, cmssp, u, v)

    # Diagnostic info
    total_cost = 0.0
    successful = false
    mode_switches = 0

    start_metadata = zero(metadatatype(solver.graph_tracker))

    next_target = nothing

    while t <= solver.max_steps

        @debug "Curr state - ",solver.curr_state
        @debug "Curr time - ", t

        if plan == true
            open_loop_plan!(cmssp, solver.curr_state, t, edge_weight, solver.heuristic,
                                  solver.goal_modes, solver.graph_tracker, start_metadata)
            last_plan_time = t
            plan = false
            next_target = solver.graph_tracker.curr_graph.vertices[solver.graph_tracker.curr_soln_path_idxs[2]]
        else
            if is_inf_hor(next_target.tp) == false
                valid_update = update_next_target!(cmssp, solver, t)

                if valid_update == false
                    plan = false
                    continue
                end
                next_target = solver.graph_tracker.curr_graph.vertices[solver.graph_tracker.curr_soln_path_idxs[2]]
            end
        end

        # now do work for macro-actions
         # Second vertex in full plan
        # start_metadata = next_target.metadata
        @debug "Next Target - ", next_target

        # Check if 'terminal' state as per modal MDP
        @assert solver.curr_state.mode == next_target.pre_bridge_state.mode

        relative_time = mean(next_target.tp) - t
        @debug "Time to next switch - ".relative_time

        if relative_time < 2

            curr_action = get_bridging_action(next_target)

        else
            if is_inf_hor(next_target.tp)
                policy = solver.modal_policies[solver.curr_state.mode].inf_hor_policy
            else
                policy = solver.modal_policies[solver.curr_state.mode].fin_hor_policy
            end
            curr_action = get_best_intramodal_action(policy,
                                                     t,
                                                     next_target.tp,
                                                     solver.curr_state.continuous,
                                                     next_target.pre_bridge_state.continuous)
            if curr_action != nothing
                curr_action = curr_action.action
            end
        end

        @debug "Current Action - ", curr_action

        if curr_action != nothing
            temp_full_action = CMSSPAction{AD,AC}(curr_action, 1)
        else
            temp_full_action = nothing
        end

        # If curr_action = NOTHING, means interrupt
        # So simulator should handle NOTHING
        # Will typically be bound to other environment variables
        # update_context_set!(cmssp, solver.rng)
        (new_state, reward, failed_mode_switch, timeout) = simulate_cmssp!(cmssp, solver.curr_state, temp_full_action, t, solver.rng)
        t = t+1
        total_cost += -1.0*reward

        if curr_action == nothing || failed_mode_switch || t - last_plan_time > solver.replan_time_threshold ||
            new_state.mode != solver.curr_state.mode
            plan = true
        end

        if new_state.mode != solver.curr_state.mode
            start_metadata = next_target.metadata
            mode_switches += 1
        end

        # Update current state and context
        solver.curr_state = new_state
        update_context_set!(cmssp, solver.rng)

        @debug new_state

        # Check terminal conditions
        if POMDPs.isterminal(cmssp, solver.curr_state)
            successful = true
            break
        end

        if timeout
            successful = false
            break
        end
    end

    # TODO : Diagnostics
    return total_cost, t, successful, mode_switches
end


text:
using SlurmAddAllocatedProcs
using Test

@testset "SlurmAddAllocatedProcs.jl" begin
    @test_throws ErrorException addprocs_slurm_allocated()
    @testset "parse_slurm_tasks_per_node" begin
        @test parse_slurm_tasks_per_node("28(x3)") == 28*3
        @test parse_slurm_tasks_per_node("28(x3),1") == 28*3+1
        @test parse_slurm_tasks_per_node("2,28(x2),1") == 2+28*2+1
        @test parse_slurm_tasks_per_node("2,28(x2),1,28(x4)") == 2+28*2+1+28*4
    end
end

# tests from SlurmClusterManager .jl
# https://github.com/kleinhenz/SlurmClusterManager.jl/blob/master/test/runtests.jl
JULIAPATH = joinpath(Sys.BINDIR, Base.julia_exename())

function testslurmjob(; ntasks = nothing, nnodes = nothing, cpupertask = nothing)
    project_path = abspath(joinpath(@__DIR__, ".."))
    # project should point to top level dir so that SlurmAddAllocatedProcs is available to script.jl
    script_path = joinpath(@__DIR__, "script.jl")
    juliacmd = "$JULIAPATH --startup=no --project=$project_path $script_path"
    # submit job
    mktempdir(@__DIR__) do dir
        cd(dir)
        @show pwd()
        open("jobscript.slurm", "w") do f
            println(f, "#!/bin/bash")
            println(f, juliacmd)
        end
        tasks = ntasks === nothing ? `` : `-n $ntasks`
        nodes = nnodes === nothing ? `` : `-N $nnodes`
        cpuspernode = cpupertask === nothing ? `` : `--cpus-per-task $cpupertask`
        outfile = joinpath(dir, "test.out")
        sbatchcmd = `sbatch --export=ALL --parsable $tasks $nodes $cpuspernode -o $outfile jobscript.slurm`
        jobid = read(sbatchcmd, String)
        println("jobid = $jobid")

        # get job state from jobid
        getjobstate = jobid -> read(`sacct -j $jobid --format=state --noheader`, String)

        # wait for job to complete
        status = timedwait(300.0, pollint=10.0) do
          state = getjobstate(jobid)
          state == "" && return false
          state = first(split(state)) # don't care about jobsteps
          println("jobstate = $state")
          return state == "COMPLETED" || state == "FAILED"
        end

        # check that job finished running within timelimit (either completed or failed)
        @test status == :ok

        # print job output
        output = read(outfile, String)
        println("script output:")
        println(output)

        state = getjobstate(jobid) |> split

        # check that everything exited without errors
        @test all(state .== "COMPLETED")
    end
end

# test job submission if slurm is available
if Sys.which("sinfo") !== nothing
    testslurmjob(ntasks = 4) # 4 tasks on 1 node
    testslurmjob(nnodes = 2) # all workers on 2 nodes
    testslurmjob(ntasks = 4, nnodes = 1) # 4 tasks on 1 nodes
    testslurmjob(ntasks = 4, nnodes = 2) # 2 tasks each on 2 nodes
    testslurmjob(ntasks = 2, cpupertask = 2) # 2 tasks on 1 node with 2 CPUs per task
    testslurmjob(ntasks = 2, nnodes = 2, cpupertask = 2) # 1 task each on 2 nodes with 2 CPUs per task
end


text:
using AbstractMCMC
using Distributions
using LinearAlgebra
using MCMCChains
using NestedSamplers
using StableRNGs
using StatsBase
using StatsFuns
using Test

rng = StableRNG(8425)
AbstractMCMC.setprogress!(get(ENV, "CI", "false") == "false")

include("utils.jl")

@testset "Bounds" begin include("bounds/bounds.jl") end
@testset "Proposals" begin include("proposals/proposals.jl") end
@testset "Sampler" begin include("sampler.jl") end
@testset "Sampling" begin include("sampling.jl") end
@testset "Models" begin include("models.jl") end