text:
## modeling of TCR-MHC-CD4 kinetics from adhesion frequency experiments
##
using Revise
src_dir = "C:/Users/stravaglino3/Downloads/julia_code/src/"
# includet(src_dir*"model_fit_global.jl")
includet(src_dir*"model_fit_global_v2.jl")
includet(src_dir*"models.jl")
includet(src_dir*"helper_functions.jl")
includet(src_dir*"global_plotting.jl")
using Plots
plot(0:10)

# experimental conditions
##
# Site densities (from flow cytometry)
densities_T =   [25.0, 38.0, 1.0]       # TCR, pMHC, Ac
densities_TC1 = [15.0, 38.0, 10.0, 1.0] # TCR, pMHC, CD4, Ac
densities_TC2 = [8.0,  38.0, 17.0, 1.0]
densities_TC3 = [3.0,  38.0, 22.0, 1.0]
densities_C =   [1700., 850., 1.0]       # CD4, pMHC, A

# Importing experimental data
##
using CSV
data_dir = "C:/Users/stravaglino3/Downloads/julia_code/TCR_CD4_project/data/"
tcr_data = CSV.read(data_dir*"tcr_only_data.csv", DataFrame)
tcr_cd4_1 = CSV.read(data_dir*"tcr_cd4_1.csv", DataFrame)
tcr_cd4_2 = CSV.read(data_dir*"tcr_cd4_2.csv", DataFrame)
tcr_cd4_3 = CSV.read(data_dir*"tcr_cd4_3.csv", DataFrame)
cd4_data = CSV.read(data_dir*"cd4_data_from_muaz.csv", DataFrame)
t_t,   Pa_t,   n_se_t   = clean_data(tcr_data)
t_tc1, Pa_tc1, n_se_tc1 = clean_data(tcr_cd4_1)
t_tc2, Pa_tc2, n_se_tc2 = clean_data(tcr_cd4_2)
t_tc3, Pa_tc3, n_se_tc3 = clean_data(tcr_cd4_3)
t_c,   Pa_c,   n_se_c   = clean_data(cd4_data)

n_t   = Pa2n(Pa_t)
n_tc1 = Pa2n(Pa_tc1)
n_tc2 = Pa2n(Pa_tc2)
n_tc3 = Pa2n(Pa_tc3)
n_c   = Pa2n(Pa_c)

scatter(n_se_t[:,1],n_se_t[:,2],yerror=n_se_t[:,3],markercolor=:black,markersize=5)
scatter!(n_se_tc1[:,1],n_se_tc1[:,2],yerror=n_se_tc1[:,3],markercolor=:blue,markersize=5)
scatter!(n_se_tc2[:,1],n_se_tc2[:,2],yerror=n_se_tc2[:,3],markercolor=:red,markersize=5)
scatter!(n_se_tc3[:,1],n_se_tc3[:,2],yerror=n_se_tc3[:,3],markercolor=:cyan,markersize=5,legend=false)
# scatter!(n_se_c[:,1],n_se_c[:,2],yerror=n_se_c[:,3],markercolor=:orange,markersize=5)

# TCR-MHC fit (bimolecular)
##
using BlackBoxOptim, Measurements

tspan = (0,maximum(t_t)+1)
u₀ = [0.0]
bounds = (1.0e-6,100.0)
pinit = [-4.0,0.0]

tcr = Data([n_t],[t_t])
M_tcr = M_global(bimolecular!,u₀,tspan,[densities_T],[[]],0.0,Rodas5)

rates_tcr,hes_tcr,l_tcr = optimization(M_tcr,tcr,bounds)
σ_tcr = hessian2σ(hes_tcr,l_tcr,tcr.n[1],rates_tcr)
@show k_t = rates_tcr .± σ_tcr

f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_T, [])
Color = :black; plot_sum(rates_tcr,:black)

# CD4-MHC fit (bimolecular)
##
tspan = (0,maximum(t_c)+1)
u₀ = [0.0]
bounds = (1.0e-9,100.0)
pinit = [-4.0,0.0]

cd4 = Data([n_c],[t_c])
M_cd4 = M_global(bimolecular!,u₀,tspan,[densities_C],[[]],0.0,Rodas5)

rates_cd4,hes_cd4,l_cd4 = optimization(M_cd4,cd4,bounds)
σ_cd4 = hessian2σ(hes_cd4,l_cd4,cd4.n[1],rates_cd4)
@show k_c = rates_cd4 .± σ_cd4

# TCR-MHC-CD4 fit (1 path)
##
u₀ = [0.0, 0.0]
cons_tcr = rates_tcr
cons_tc = vcat(cons_tcr, rates_cd4)
tspan = (0,maximum(t_tc1)+1)
pinit = [0.0,0.0]

tc1 = Data([n_tc1],[t_tc1])
tc2 = Data([n_tc2],[t_tc2])
tc3 = Data([n_tc3],[t_tc3])
M_tc1 = M_global(trimolecular!,u₀,tspan,[densities_TC1],[cons_tcr],0.0,Rodas5)
M_tc2 = M_global(trimolecular!,u₀,tspan,[densities_TC2],[cons_tcr],0.0,Rodas5)
M_tc3 = M_global(trimolecular!,u₀,tspan,[densities_TC3],[cons_tcr],0.0,Rodas5)

rates_tc1,hes_tc1,l_tc1 = optimization(M_tc1,tc1,bounds)
σ_tc1 = hessian2σ(hes_tc1,l_tc1,tc1.n[1],rates_tc1)
rates_tc2,hes_tc2,l_tc2 = optimization(M_tc2,tc2,bounds)
σ_tc2 = hessian2σ(hes_tc2,l_tc2,tc2.n[1],rates_tc2)
rates_tc3,hes_tc3,l_tc3 = optimization(M_tc3,tc3,(1e-6,1e3))
σ_tc3 = hessian2σ(hes_tc3,l_tc3,tc3.n[1],rates_tc3)

@show k_tc1 = rates_tc1 .± σ_tc1
@show k_tc2 = rates_tc2 .± σ_tc2
@show k_tc3 = rates_tc3 .± σ_tc3

# plot 1 path trimolecular fits
##
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_TC1, rates_tcr)
Color = :blue; plot_sum(rates_tc1,:blue)
# plot_sum([0., 0.])
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_TC2, rates_tcr)
Color = :red; plot_sum(rates_tc2,:red)
# plot_sum([0., 0.])
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_TC3, rates_tcr)
Color = :cyan; plot_sum(rates_tc3,:cyan)

# plot_sum([0., 0.])
#
# TCR-MHC-CD4 fit (1 path, with diffusion?)
##
# u₀ = [0., 0., 0.]
# f = (du,u,p,t) -> trimolecular_1path_diff!(du,u,p,t,densities_TC1, rates_tcr)
# plot_sum([1e-3,1e-3,1,1e-3])
# pinit = [0.0,0.0,0.0,0.0]
# M_tc1_diff = M_global(trimolecular_1path_diff!,u₀,tspan,[densities_TC1],[cons_tcr],0.0,Rodas5)
# rates_tc1_diff,hes_tc1_diff,l_tc1_diff = optimization(M_tc1_diff,tc1,(1e-8,1e3))
# f = (du,u,p,t) -> trimolecular_1path_diff!(du,u,p,t,densities_TC1, rates_tcr)
# plot_sum(rates_tc1_diff)


text:
module HttpHelpers

export state_handler, exception_handling_handler, compress_handler, Route, Router, add_route!

import HTTP, CodecZlib

include("handlers.jl")
include("router.jl")

end

ext: jl


text:
@testset "Sector Performance" begin
    data = sector_performance()
    @test typeof(data) === Dict{String, Any}
    @test length(data) === 11
    sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
end


text:
#!/usr/bin/env julia

input = join(sort(split(ARGS[1], "")))

upper = IOBuffer()
lower = IOBuffer()
even = IOBuffer()
odd = IOBuffer()

for c in input
    if 'a' <= c && c <= 'z'
        print(lower, c)
    elseif 'A' <= c && c <= 'Z'
        print(upper, c)
    elseif occursin(c, "02468")
        print(even, c)
    else
        print(odd, c)
    end
end

print(String(take!(lower)))
print(String(take!(upper)))
print(String(take!(even)))
print(String(take!(odd)))
println()


text:
using StingerGraphs
using Base.Test

function generate!(edges, numvertices::Integer, numedges::Integer)
    for i in 1:numedges
        src, dst = rand(0:numvertices-1), rand(0:numvertices-1)
        edges[:,i] = [0, src, dst, 1, 0]
    end
end

function bench(numvertices, numedges)
    info("Allocating a STINGER")
    s = Stinger()
    edges = zeros(Int64, 5, numedges)
    info("Generating the edges")
    @time generate!(edges, numvertices, numedges)
    info("Inserting $(size(edges,2)) edges")
    @time insert_edges!(s, edges, numedges)
    @test consistency_check(s, numvertices) == true
    removals = edges[:,randperm(numedges)]
    removals = unique(removals, 2)
    @show size(removals)
    info("Removing 1/8 of the edges")
    @time remove_edges!(s, removals, ceil(Int, numedges/8))
    @test consistency_check(s, numvertices) == true
    return s, edges, removals
end

info("Setting STINGER Parameters")
ENV["STINGER_MAX_MEMSIZE"] = 20000000
info("small run to precompile")
bench(100, 100)
info("real benchmark: nv=1000,ne=4000")
bench(1000, 4000)


text:
js_files=["https://www.ishares.com/blk-one-assets/include/bundles/minified-78128b72721fb029b708c5000d861e38.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-78c7fb04771253db3744d40b138f80c5.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-5e5b6f139bb356f365ab24198b7c33b6.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-730b1fd240cc6f60184b43b3d1c766b7.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-6dbcdbeaa354b3ff264700d69e4fc333.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-87cd6f2ff1171cdc6b51e64d29474a30.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-87dacaf861cd7d4f0231f3634f690fd4.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-d8bad6548037e489c44baaf473932769.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-be85168882f9d6a14c02162feb5e5116.js",
          "https://www.ishares.com/blk-one-assets/include/bundles/minified-10612e8f01c8c191917dfaed1fbc78bd.js",
          "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/runtime.458556a34b891ea32398.js",
          "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/polyfills-es5.91147ef724d245a1fc29.js",
          "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/polyfills.921e22471ff1fb602a28.js" ,
          "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/scripts.69c39fe5fecacc5138f1.js" ,
          "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/main.d653d910ff1c4a5b9a4e.js" ]

css_files=["https://www.ishares.com/blk-one-assets/include/bundles/minified-67bdb0a17c5be3039d1755c7ecd2e2e6.css",
           "https://www.ishares.com/blk-one-assets/include/bundles/minified-383ef49ce5f513404820bffb3aff14e8.css",
           "https://www.ishares.com/blk-one-assets/include/bundles/minified-e6ac0be7d09297f70e79618c7732b56a.css",
           "https://www.ishares.com/blk-one-assets/include/bundles/minified-2ce8052867cc56d3708795446ac12ece.css",
           "https://www.ishares.com/blk-one-assets/include/bundles/minified-264c22a9ceb75519e48d564234a2a285.css",
           "https://www.ishares.com/blk-one-assets/include/bundles/minified-d8e35de366c069e4a084c01fc7a4abbd.css",
           "https://www.ishares.com/blk-one-assets/include/bundles/minified-ae9f873241d94ebb3e9288de1c239db2.css",
           "https://www.ishares.com/blk-one-assets/include/bundles/minified-bcc17d66ca66b2742f8c240497925528.css",
           "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/styles.bac52fec26df0e643e60.css"]

using RemoteFiles
@RemoteFileSet JAVASCRIPTS "JAVASCRIPTS" begin

    js1 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-78128b72721fb029b708c5000d861e38.js"
    js2 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-78c7fb04771253db3744d40b138f80c5.js"
    js3 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-5e5b6f139bb356f365ab24198b7c33b6.js"
    js4 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-730b1fd240cc6f60184b43b3d1c766b7.js"
    js5 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-6dbcdbeaa354b3ff264700d69e4fc333.js"
    js6 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-87cd6f2ff1171cdc6b51e64d29474a30.js"
    js7 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-87dacaf861cd7d4f0231f3634f690fd4.js"
    js8 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-d8bad6548037e489c44baaf473932769.js"
    js9 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-be85168882f9d6a14c02162feb5e5116.js"
    js10 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-10612e8f01c8c191917dfaed1fbc78bd.js"
    js11 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/runtime.458556a34b891ea32398.js"
    js12 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/polyfills-es5.91147ef724d245a1fc29.js"
    js13 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/polyfills.921e22471ff1fb602a28.js"
    js14 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/scripts.69c39fe5fecacc5138f1.js"
    js15 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/main.d653d910ff1c4a5b9a4e.js"

end

@RemoteFileSet CSS "CSS" begin

    css1 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-67bdb0a17c5be3039d1755c7ecd2e2e6.css"
    css2 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-383ef49ce5f513404820bffb3aff14e8.css"
    css3 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-e6ac0be7d09297f70e79618c7732b56a.css"
    css4 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-2ce8052867cc56d3708795446ac12ece.css"
    css5 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-264c22a9ceb75519e48d564234a2a285.css"
    css6 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-d8e35de366c069e4a084c01fc7a4abbd.css"
    css7 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-ae9f873241d94ebb3e9288de1c239db2.css"
    css8 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/bundles/minified-bcc17d66ca66b2742f8c240497925528.css"
    css9 = @RemoteFile "https://www.ishares.com/blk-one-assets/include/app-container/ipx-compare-17795c3cc48/styles.bac52fec26df0e643e60.css"

end

download(CSS)


text:
module mean_functions

using ..gpflow
import ..gpflow: instantiate!, predict_f, predict_f_samples
export
    Additive,
    Constant,
    Identity,
    Linear,
    MeanFunction,
    instantiate!

mutable struct Additive{T1,T2} <: AbstractMeanFunction
    first_part::T1
    second_part::T2
    o::Union{PyObject, Nothing}
end

function Additive(first_part, second_part)
    out = Additive(first_part, second_part, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Additive, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.mean_functions.Additive(o.first_part, o.second_part)
    return o.o
end

function (o::Additive)(X)
    if  !(typeof(o.o)<:PyObject)
        instantiate!(o)
    end
    return o.o(X)
end

"""
    y_i = c,,
"""
mutable struct Constant{T1} <: AbstractMeanFunction
    c::T1
    o::Union{PyObject, Nothing}
end

function Constant(;c=nothing)
    out = Constant(c, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Constant, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.mean_functions.Constant(;c=o.c)
    return o.o
end

function (o::Constant)(X)
    if  !(typeof(o.o)<:PyObject)
        instantiate!(o)
    end
    return o.o(X)
end

abstract type AbstractLinear <: AbstractMeanFunction end

"""
    y_i = x_i
"""
mutable struct Identity{T1} <: AbstractLinear
    input_dim::T1
    o::Union{PyObject, Nothing}
end

function Identity(;input_dim=nothing)
    out = Identity(input_dim, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Identity, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.mean_functions.Identity(;input_dim=o.input_dim)
    return o.o
end

function (o::Identity)(X)
    if  !(typeof(o.o)<:PyObject)
        instantiate!(o)
    end
    return o.o(X)
end

"""
    y_i = A x_i + b
"""
mutable struct Linear{T1,T2} <: AbstractLinear
    A::T1
    b::T2
    o::Union{PyObject, Nothing}
end

function Linear(;A=nothing, b=nothing)
    out = Linear(A, b, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{Linear, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.mean_functions.Linear(;A=o.A, b=o.b)
    return o.o
end

function (o::Linear)(X)
    if  !(typeof(o.o)<:PyObject)
        instantiate!(o)
    end
    return o.o(X)
end

"""
    The base mean function class.
    To implement a mean function, write the __call__ method. This takes a
    tensor X and returns a tensor m(X). In accordance with the GPflow
    standard, each row of X represents one datum, and each row of Y is computed
    independently for each row of X.

    MeanFunction classes can have parameters, see the Linear class for an
    example.
"""
mutable struct MeanFunction{T1} <: AbstractMeanFunction
    name::T1
    o::Union{PyObject, Nothing}
end

function MeanFunction(;name=nothing)
    out = MeanFunction(name, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{MeanFunction, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.mean_functions.MeanFunction(;name=o.name)
    return o.o
end

function (o::MeanFunction)(X)
    if  !(typeof(o.o)<:PyObject)
        instantiate!(o)
    end
    return o.o(X)
end

mutable struct SwitchedMeanFunction{T1} <: AbstractMeanFunction
    meanfunction_list::T1
    o::Union{PyObject, Nothing}
end

function SwitchedMeanFunction(meanfunction_list)
    out = SwitchedMeanFunction(meanfunction_list, nothing)
    instantiate!(out)
    out
end

function instantiate!(o::Union{SwitchedMeanFunction, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.mean_functions.SwitchedMeanFunction([mf.o for mf in o.meanfunction_list])
    return o.o
end

function (o::SwitchedMeanFunction)(X)
    if  !(typeof(o.o)<:PyObject)
        instantiate!(o)
    end
    return o.o(X)
end
end # module


text:
module Ising

export IsingModel, update!, simulate!

"""
Ising model with periodic boundary conditions and arbitrary number of dimensions.

To simplify the calculations, the units are set in way that J=1 and kB=1
(the boltzman constant). So, the energy is in units of J and the temperature
is in units of J/kB. β can be thought of as J/(kB*T) or 1/T in this unit system.

`dims` in the constructor takes the number of particles in each dimension,
so that the number of elements of `dims` is equal to the dimensionality of the system.

The simulation here utilizes the Metropolis Monte Carlo method. Also, a checkerboard
coloring pattern is used for flipping many spins at the same time
"""
mutable struct IsingModel{N} <: Integer
    spins::Array{Int8, N}
    energy::Float64
    magnetization::Float64

    β::Float64
    mass::Float64
    dims::NTuple{N, Integer}

    # pre-calculated values for possible state changes
    _probs::Vector{Float64}

    # spins broken down into a checkerboard pattern. spins of each color
    # do not interact with each other, so they can be flipped simultatiously
    _black::SubArray{Int8, 1, Vector{Int8}, Tuple{Vector{Int64}}, false} # even index sum
    _white::SubArray{Int8, 1, Vector{Int8}, Tuple{Vector{Int64}}, false} # odd index sum

    _blackneighbors::SubArray{Int8, 2, Vector{Int8}, Tuple{Matrix{Int64}}, false}
    _whiteneighbors::SubArray{Int8, 2, Vector{Int8}, Tuple{Matrix{Int64}}, false}

    function IsingModel(β, dims::Integer...)
        spins = rand([Int8(-1), Int8(1)], dims...)

        mass = prod(dims)
        indices = CartesianIndices(spins)
        linindices = LinearIndices(spins)
        neighbors = _get_neighbors(indices, linindices, dims, mass)

        energy = -sum(reshape(spins, mass, 1) .* sum(spins[neighbors], dims=2)) / 2
        magnetization = sum(spins) / mass
        probs = exp.(-4β .* collect(-length(dims):length(dims)))

        isblack = sum.(Tuple.(indices)) .% 2 .== 0
        blackindices = linindices[findall(isblack)]
        whiteindices = linindices[findall(map(!, isblack))]

        black = @view spins[blackindices]
        white = @view spins[whiteindices]
        blackneighbors = @view spins[neighbors[blackindices, :]]
        whiteneighbors = @view spins[neighbors[whiteindices, :]]

        return new{length(dims)}(spins, energy, magnetization, β, mass, dims,
            probs, black, white, blackneighbors, whiteneighbors)
    end
end

"""
get the indices of the neighbors of every point in `indices`
"""
@inline function _get_neighbors(indices::CartesianIndices, linindices::LinearIndices,
        dims::Tuple{Vararg{Integer}}, mass::Integer)
    neighborindices = Matrix{Int}(undef, mass, 0)
    for dim in 1:length(dims)
        basevec = zeros(Integer, length(dims))
        basevec[dim] = 1
        step = CartesianIndex(basevec...)
        basevec[dim] = dims[dim]
        modulus = CartesianIndex(basevec...)

        # collecting the CartesianIndices into a vector to enable setindex
        nextneighbors = collect(indices .+ step)
        prevneighbors = collect(indices .- step)
        # referencing fixes broadcasting issues
        selectdim(prevneighbors, dim, 1) .+= Ref(modulus)
        selectdim(nextneighbors, dim, dims[dim]) .-= Ref(modulus)
        neighborindices = hcat(
            neighborindices,
            reshape(linindices[prevneighbors], mass, 1),
            reshape(linindices[nextneighbors], mass, 1))
    end
    return neighborindices
end

"""
Update an `IsingModel` for one time step, where all spins are updated once.
The operations are vectorized on the two "black" and "white" indices.
"""
function update!(ising::IsingModel)
    blackΔEs = vec(2 * ising._black .* sum(ising._blackneighbors, dims=2))
    blackprobs = ising._probs[blackΔEs .÷ 4 .+ (length(ising.dims) + 1)]
    flipblack = rand(length(blackprobs)) .<= blackprobs
    ising._black[flipblack] *= -1

    whiteΔEs = vec(2 * ising._white .* sum(ising._whiteneighbors, dims=2))
    whiteprobs = ising._probs[whiteΔEs .÷ 4 .+ (length(ising.dims) + 1)]
    flipwhite = rand(length(whiteprobs)) .<= whiteprobs
    ising._white[flipwhite] *= -1

    ising.magnetization += 2 / ising.mass * (
        sum(ising._black[flipblack]) + sum(ising._white[flipwhite]))
    ising.energy += sum(blackΔEs[flipblack]) + sum(whiteΔEs[flipwhite])
end

"""
Update an `IsingModel` for `time` steps. If `fullstats` is `false`, calculate and return
the mean and mean squared energy and magnetization (in the order ⟨E⟩, ⟨m⟩, ⟨E²⟩, ⟨m²⟩)
over the simulation time period. If `fullstats` is `true`, record and return the energy
and magnetization for every time step of the simulation (in the order Es, ms)

note: m is averaged over the absolute value
"""
function simulate!(isingmodel::IsingModel, time::Integer; fullstats::Bool=false)
    if fullstats
        energy = Vector{Float64}(undef, time)
        magnetization = Vector{Float64}(undef, time)

        for t in 1:time
            update!(isingmodel)
            energy[t] = isingmodel.energy
            magnetization[t] = isingmodel.magnetization
        end

        return energy, magnetization
    else
        Ē, m̄, Ē², m̄² = zeros(4)

        for _ in 1:time
            update!(isingmodel)
            Ē += isingmodel.energy / time
            m̄ += abs.(isingmodel.magnetization) / time
            Ē² += isingmodel.energy^2 / time
            m̄² += isingmodel.magnetization^2 / time
        end

        return Ē, m̄, Ē², m̄²
    end
end

using Statistics: mean, var

"""
    Calculate the correlation function of an `IsingModel` along the dimension `dim`,
    taking into account the periodic boundary conditions.
"""
function correlation(ising::IsingModel; dim=1)
    variance = var(ising.spins, mean=ising.magnetization)
    if variance == 0
        return zeros(ceil(Integer, size(ising.spins, dim) / 2))
    end

    correlations = Vector{Float64}(undef, ceil(Integer, size(ising.spins, dim) / 2))
    shiftvec = zeros(length(size(ising.spins)))
    for shift in 1:length(correlations)
        shiftvec[dim] = shift
        correlations[shift] = mean(ising.spins .* circshift(ising.spins, shiftvec))
    end

    # demean and normalize
    correlations .-= ising.magnetization ^ 2
    correlations ./= variance

    return correlations
end

"""
    Calculate the correlation length of an `IsingModel`, taking into account the periodic
    boundary conditions. Correlation data smaller than `thresh` is discarded (this is
    considered as "noise"). To counteract cases where the correlation length in one
    direction diverges, in these cases the correlation length is calculated for two
    directions and the minimum of the two correlation lengths is returned.
"""
function corlen(ising::IsingModel; thresh::Float64 = 0.1)
    correlations = correlation(ising, dim=1)
    threshindex = findfirst(x -> abs(x) <= thresh, correlations)

    if isnothing(threshindex)
        threshindex = length(correlations)

        # newcorr = correlation(ising, dim=2)
        # newthresh = findfirst(x -> abs(x) <= thresh, newcorr)
        # if isnothing(newthresh)
        #     newthresh = length(newcorr)
        # elseif newthresh < 3
        #     return 0.0
        # end

        # negativeξ⁻¹ = ([collect(1:threshindex) ones(threshindex)]
        #     \ log.(abs.(correlations[1:threshindex])))[1]
        # ξ = isnan(negativeξ⁻¹) ? 0.0 : 1/negativeξ⁻¹
        # if ξ < 0.0
        #     return 0.0
        # end

        # newnegativeξ⁻¹ = ([collect(1:newthresh) ones(newthresh)]
        #     \ log.(abs.(newcorr[1:newthresh])))[1]
        # newξ = isnan(newnegativeξ⁻¹) ? 0.0 : 1/newnegativeξ⁻¹
        # if newξ < 0.0
        #     return 0.0
        # end

        # return min(ξ, newξ)
    elseif threshindex < 3
        return 0.0
    end

    negativeξ⁻¹ = ([collect(1:threshindex) ones(threshindex)]
        \ log.(abs.(correlations[1:threshindex])))[1]
    if isnan(negativeξ⁻¹)
        return 0.0
    end

    ξ = -1/negativeξ⁻¹
    if ξ < 0
        return 0.0
    end

    return ξ
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
# Ion species
@enum IonSpecies::UInt128 begin
    NonIonic    = 1 << 0
    Sodium      = 1 << 1
    Potassium   = 1 << 2
    Chloride    = 1 << 3
    Calcium     = 1 << 4
    Cationic    = 1 << 5
    Anionic     = 1 << 6
end

const Ca = Calcium
const Na = Sodium
const K = Potassium
const Cl = Chloride
const Mixed = const Leak = NonIonic

const PERIODIC_SYMBOL = IdDict(Na => :Na, K  => :K, Cl => :Cl, Ca => :Ca, Leak => :l)

struct IonConcentration
    ion::IonSpecies
    loc::Location
end

const Concentration = IonConcentration

function IonConcentration(
    ion::IonSpecies,
    val = nothing;
    location::Location = Inside,
    dynamic::Bool = false,
    name::Symbol = PERIODIC_SYMBOL[ion]
)

    sym = Symbol(name,(location == Inside ? "ᵢ" : "ₒ"))
    var = dynamic ? only(@variables $sym(t)) : only(@parameters $sym)
    var = setmetadata(var,  IonConcentration, IonConcentration(ion, location))
    if !isnothing(val)
        if val isa Molarity
            var = setmetadata(var, ConductorUnits, unit(val))
            # FIXME: use proper unit checking
            raw_val = ustrip(Float64, µM, val)
            var = setdefault(var, raw_val)
            return var
        else
            var = setdefault(var, val)
            return var
        end
    end
    return var
end

isconc(x) = hasmetadata(value(x), IonConcentration)
getconc(x) = isconc(x) ? getmetadata(value(x), IonConcentration) : nothing

struct IonCurrent
    ion::IonSpecies
    agg::Bool
end

function IonCurrent(
    ion::IonSpecies,
    val = nothing;
    aggregate::Bool = false,
    dynamic::Bool = true,
    name::Symbol = Symbol("I", PERIODIC_SYMBOL[ion])
)
    var = dynamic ? only(@variables $name(t)) : only(@parameters $name)
    var = setmetadata(var, IonCurrent, IonCurrent(ion, aggregate))
    if !isnothing(val)
        if val isa Current
            var = setmetadata(var, ConductorUnits, unit(val))
            # FIXME: use proper unit checking
            raw_val = ustrip(Float64, µA, val)
            var = setdefault(var, raw_val)
            return var
        else
            var = setdefault(var, val)
            return var
        end
    end
    return var
end

iscurrent(x) = hasmetadata(value(x), IonCurrent)
iscurrent(x::IonCurrent) = true
getcurrent(x) = iscurrent(x) ? getmetadata(value(x), IonCurrent) : nothing
getcurrent(x::IonCurrent) = x
getion(x::IonCurrent) = getfield(getcurrent(x), :ion)
isaggregate(x) = iscurrent(x) ? getfield(getcurrent(x), :agg) : false

struct EquilibriumPotential
    ion::IonSpecies
end

const Equilibrium = EquilibriumPotential

function EquilibriumPotential(ion::IonSpecies, val; dynamic = false, name::Symbol = PERIODIC_SYMBOL[ion])
    sym = Symbol("E", name)
    var = dynamic ? only(@variables $sym(t)) : only(@parameters $sym)
    var = setmetadata(var, EquilibriumPotential, EquilibriumPotential(ion))
    if !isnothing(val)
        if val isa Voltage
            var = setmetadata(var, ConductorUnits, unit(val))
            # FIXME: do proper unit checking
            raw_val = ustrip(Float64, mV, val)
            var = setdefault(var, raw_val)
            return var
        else
            var = setdefault(var, val)
            return var
        end
    end
    return var
end

isreversal(x) = hasmetadata(value(x), EquilibriumPotential)
getreversal(x) = isreversal(x) ? getmetadata(value(x), EquilibriumPotential) : nothing
getion(x::EquilibriumPotential) = getfield(x, :ion)

function getion(x)
    iscurrent(x) && return getion(getcurrent(x))
    isreversal(x) && return getion(getreversal(x))
    return nothing
end

# Alternate constructor; this needs to be better...
function Equilibria(equil::Vector)
    out = Num[]
    for x in equil
        x.first isa IonSpecies || throw("Equilibrium potential must be associated with an ion type.")
        if x.second isa Tuple
            tup = x.second
            tup[2] isa Symbol || throw("Second tuple argument for $(x.first) must be a symbol.")
            push!(out, Equilibrium(x.first, tup[1], dynamic = tup[1] isa Voltage ? false : true, name = tup[2]))
        else
            push!(out, Equilibrium(x.first, x.second, dynamic = x.second isa Voltage ? false : true))
        end
    end
    return out
end


text:
using Parameters, UUIDs, JSON

JSON.lower(x::UUID) = string(x)

@with_kw struct JSXElement
    hash::UUID = uuid4()
    tag::String = "div"
    props::Union{Dict{String, Any}, Nothing} = nothing
    children::Union{Vector{JSXElement}, Nothing} = nothing
    data_type::Union{DataType, Nothing} = nothing
    component_type::String # Input? Output? Static?
end

# BindSet -> Input
@with_kw struct UpdateRule
    # Should be the hash of `BindSet`
    bind_set::UUID = uuid4()
    input::Vector{UUID}
    output::UUID
end

function parse_tree(x::InputComponent)
    local update_rules = Vector{UpdateRule}()
    local components = Vector{Component}()

    local props = isnothing(x.parameters) ? x.parameters : Dict{String, Any}()

    if !isnothing(x.default)
        props["value"] = x.default
    end

    append!(components, [x])

    local jsx_tree = JSXElement(
        hash = x.hash,
        tag = x.component,
        props = props,
        children = nothing,
        component_type = "input",
        data_type = typeof(x).parameters[1]
    )

    return ParsedJSXTree(
        jsx_tree = jsx_tree,
        update_rules = update_rules,
        components = components
    )
end

function parse_tree(x::OutputComponent)
    local update_rules = Vector{UpdateRule}()
    local components = Vector{Component}()

    append!(components, [x])

    local jsx_tree = JSXElement(
        hash = x.hash,
        tag = x.component,
        props = x.parameters,
        children = nothing,
        component_type = "output",
        data_type = typeof(x).parameters[1]
    )

    local inputs = Vector{UUID}()

    for bind in x.callback.bind_set.binds
        append!(inputs, [bind.component.hash])
    end

    append!(update_rules, [UpdateRule(
        bind_set = x.callback.bind_set.hash,
        input = inputs,
        output = x.hash
    )])

    return ParsedJSXTree(
        jsx_tree = jsx_tree,
        update_rules = update_rules,
        components = components
    )
end

@with_kw struct ParsedJSXTree
    jsx_tree::JSXElement
    update_rules::Vector{UpdateRule}
    components::Vector{Component}
end

function parse_tree(x::StaticComponent)
    local update_rules = Vector{UpdateRule}()
    local components = Vector{Component}()

    append!(components, [x])

    if !isnothing(x.children)
        local children = Vector{JSXElement}()

        for i in x.children
            local parsed_tree = parse_tree(i)

            append!(update_rules, parsed_tree.update_rules)
            append!(children, [parsed_tree.jsx_tree])
            append!(components, parsed_tree.components)
        end
    else
        local children = nothing
    end

    local jsx_tree = JSXElement(
        hash = x.hash,
        tag = x.component,
        props = x.parameters,
        children = children,
        component_type = "static",
        data_type = nothing
    )

    return ParsedJSXTree(
        jsx_tree = jsx_tree,
        update_rules = update_rules,
        components = components
    )
end

@with_kw struct MATTApp
    hash::UUID = uuid4()
    jsx_tree::JSXElement
    binds::Dict{UUID, Bind}
    bind_sets::Dict{UUID, BindSet}
    input_components::Dict{UUID, InputComponent}
    output_components::Dict{UUID, OutputComponent}
    update_rules::Dict{UUID, UpdateRule}
    input_bind::Dict{UUID, Vector{UUID}}
    bind_input::Dict{UUID, Vector{UUID}}
    bind_output::Dict{UUID, Vector{UUID}}
    serialized_app::String
end

function setup_app(x::ParsedJSXTree)
    setup_app(x.jsx_tree, x.update_rules, x.components)
end

function setup_app(x::StaticComponent)
    setup_app(parse_tree(x))
end

function setup_app(
    jsx_tree::JSXElement,
    update_rules::Vector{UpdateRule},
    components::Vector{T} where T <: Component)

    local input_components = Dict{UUID, InputComponent}()
    local output_components = Dict{UUID, OutputComponent}()
    local binds = Dict{UUID, Bind}()
    local bind_sets = Dict{UUID, BindSet}()
    local input_bind = Dict{UUID, Vector{UUID}}()
    local bind_input = Dict{UUID, Vector{UUID}}()
    local bind_output = Dict{UUID, Vector{UUID}}()
    local update_rules_dict = Dict{UUID, UpdateRule}()

    for rule in update_rules
        bind_input[rule.bind_set] = rule.input
        update_rules_dict[rule.bind_set] = rule

        if haskey(bind_output, rule.bind_set)
            append!(bind_output[rule.bind_set], [rule.output])
        else
            bind_output[rule.bind_set] = [rule.output]
        end

        for input_hash in rule.input
            if haskey(input_bind, input_hash)
                append!(input_bind[input_hash], [rule.bind_set])
            else
                input_bind[input_hash] = [rule.bind_set]
            end
        end
    end

    local app_ui_def = Dict(
        "jsx_tree" => jsx_tree,
        "input_bind" => input_bind,
        "bind_input" => bind_input,
    )

    for component in components
        if typeof(component) <: InputComponent
            input_components[component.hash] = component
        elseif typeof(component) <: OutputComponent
            output_components[component.hash] = component

            local bind_set_hash = component.callback.bind_set.hash
            if !haskey(bind_sets, bind_set_hash)
                bind_sets[bind_set_hash] = component.callback.bind_set
            end
        end
    end

    for bind_set in values(bind_sets)
        for bind in bind_set.binds
            !haskey(binds, bind.hash) && (binds[bind.hash] = bind)
        end
    end

    MATTApp(
        jsx_tree = jsx_tree,
        binds = binds,
        bind_sets = bind_sets,
        input_components = input_components,
        output_components = output_components,
        update_rules = update_rules_dict,
        input_bind = input_bind,
        bind_input = bind_input,
        bind_output = bind_output,
        serialized_app = json(app_ui_def),
    )
end


text:
using Base.Threads, ProgressMeter, Plots, LaTeXStrings
using FilmDataStructures, FilmUtilities, ResponseModels, OttawaSlabCell
# Open output file.
fileStream = open("photonProfile_test2.txt","w")
# Double check thread initialization.
write(fileStream, "Julia initialized with "*string(nthreads())*" threads.\n\n")
# Radial frequency to electron volts.
const hBEv = 6.582e-16
## External program settings, internal setting contained in uOttawaSlabStruct.jl
enrRng = (0.06, 1.0)
## Begin program execution.
write(fileStream, "Execution:\n")
# Builds slab structure, generating lVar and lPairs.
# Arguments are: temperature of the emitter, background temperature, divisions of NCell
# (Resolution of photon generation is thckInAs / divNCell), division of PCell
stats = @timed (lVar, lPairs) = uOttawaSlabs(800.0, 300.0, 30, 30)
write(fileStream, "Layer structure constructed in " * string(stats.time) * " s.\n")
# Storage for photon number computations.
htPairs = Array{Float64,1}(undef, size(lPairs)[2])
# Compute number of generated photons using heat transfer function.
# Switch from flux values achieved through prefactors declared in uOttawaSlabStruct.jl.
stats = @timed heatTfr!(lVar, lPairs, enrRng, htPairs)
write(fileStream, "Generated excitation profile computed " * string(stats.time) * " s.\n\n")
# Write results to file.
write(fileStream, "Depth (microns)	"*" Excitation Density (cm-2 microns-1)\n\n")
# Location of cell boundary within the slab structure.
cellSurfLoc = lVar.bdrLoc[lPairs[2,1]]
# Output calculation results.
for ind = 1 : length(htPairs)

	write(fileStream, string(round(lVar.bdrLoc[lPairs[2,ind]] - cellSurfLoc, sigdigits = 4)) * "				" * string(round(htPairs[ind], sigdigits = 4)) * "\n")
end
# Flush file stream.
close(fileStream)


text:
# update integrator when u is modified by callbacks
function handle_callback_modifiers!(integrator::DDEIntegrator)
    integrator.reeval_fsal = true # recalculate fsalfirst after applying step

    # update heap of discontinuities
    # discontinuity is assumed to be of order 0, i.e. solution x is discontinuous
    push!(integrator.opts.d_discontinuities, Discontinuity(integrator.t, 0))
end

"""
    reeval_internals_due_to_modification!(integrator::DDEIntegrator)

Recalculate interpolation data and update ODE integrator after changes by callbacks.
"""
function DiffEqBase.reeval_internals_due_to_modification!(integrator::DDEIntegrator,
            x::Type{Val{not_initialization}} = Val{true}) where not_initialization

    # update interpolation data of DDE integrator using old interpolation data
    # of ODE integrator in evaluation of history function that was calculated in
    # `perform_step!`
    if not_initialization
        addsteps!(integrator, integrator.f, true, true, true)
        # copy interpolation data to ODE integrator
        recursivecopy!(integrator.integrator.k, integrator.k)
    end

    # move ODE integrator to new time interval of DDE integrator
    integrator.integrator.t = integrator.t
    integrator.integrator.dt = integrator.dt
    if isinplace(integrator.sol.prob)
        recursivecopy!(integrator.integrator.u, integrator.u)
    else
        integrator.integrator.u = integrator.u
    end

    integrator.u_modified = false
end

# track propagated discontinuities for dependent delays

struct DiscontinuityCallback{F,D<:Discontinuity,A,R,I} <: AbstractContinuousCallback
    lags::F
    discontinuities::Vector{D}
    interp_points::Int
    abstol::A
    reltol::R
    initialize::I
    idxs::Nothing
end

"""
    DiscontinuityCallback(lags, discontinuities::Vector{<:Discontinuity};
                          [interp_points::Int=10, abstol=1e-12, reltol=0])

Callback that tracks `discontinuities` that are propagated by dependent `lags` of the form
`(u,p,t) -> lag(u,p,t)`.

Hereby a number `interp_points` of interpolation points are used to first check for
different signs of functions ``f(t) = T + lag(u(t),p,t) - t``, where ``T`` is time point of a
previous discontinuity and ``t`` is contained in the current time interval. This shows that
the current time interval contains propagated discontinuities of which the exact time point
is then determined by a root finding algorithm. The sign at the lower bound of the time
interval, i.e. at `tprev`, is set to 0 with absolute tolerance `abstol` and relative
tolerance `reltol`.
"""
function DiscontinuityCallback(lags, discontinuities::Vector{<:Discontinuity};
                               interp_points::Int=10, abstol=1e-12, reltol=0)
    DiscontinuityCallback(lags, discontinuities, interp_points, abstol,
                          reltol, initialize!, nothing)
end

# do not initialize discontinuity callback
initialize!(c::DiscontinuityCallback, u, t, integrator::DEIntegrator) = (integrator.u_modified=false)

# find time of first discontinuity in the current time interval (if existent)
function DiffEqBase.find_callback_time(integrator::DDEIntegrator, callback::DiscontinuityCallback, counter)
    # initialize time and order of first discontinuity in the current time interval
    tmin = zero(integrator.t)
    order = 0
    contains_discontinuity = false

    # update interpolation data and calculate interpolation points
    if callback.interp_points != 0
      addsteps!(integrator)
    end
    Θs = range(0, stop=one(integrator.t), length=callback.interp_points)

    for lag in callback.lags
        # define function f that calculates T + lag(u(t),p,t) - t, where t = tprev + Θ*dt is
        # a time point in the time interval of the current step and T ≤ tprev is the time
        # point of a discontinuity
        # hence roots t of this function for fixed T are propagated discontinuities
        function f(Θ, T)
            if typeof(integrator.cache) <: OrdinaryDiffEq.OrdinaryDiffEqMutableCache
                # use temporary array of mutable caches in calculation of interpolants
                tmp = integrator.cache.tmp
                OrdinaryDiffEq.ode_interpolant!(tmp, Θ, integrator, nothing, Val{0})
            else
                tmp = OrdinaryDiffEq.ode_interpolant(Θ, integrator, nothing, Val{0})
            end
            t = integrator.tprev + Θ*integrator.dt
            T + lag(tmp,integrator.p,t) - t
        end

        for d in callback.discontinuities
            # fix time of discontinuity
            T = d.t
            g(Θ) = f(Θ, T)

            # use start and end point of last time interval to check for discontinuities
            previous_condition = T + lag(integrator.uprev, integrator.p, integrator.tprev) -
                integrator.tprev
            if isapprox(previous_condition, 0, rtol=callback.reltol, atol=callback.abstol)
                prev_sign = 0
            else
                prev_sign = cmp(previous_condition, 0)
            end
            next_sign = cmp(T + lag(integrator.u, integrator.p, integrator.t) - integrator.t, 0)

            # find time of discontinuity if one exists
            t,contains_discontinuity = find_discontinuity_time(integrator, callback, prev_sign, next_sign, Θs, g)

            # update time and order of first discontinuity in the current time interval
            if contains_discontinuity && (integrator.tdir * t < integrator.tdir * tmin || iszero(tmin))
                tmin = t
                order = d.order + 1
            end
        end
    end

    tmin, order, contains_discontinuity
end

"""
    find_discontinuity_time(integrator::DDEIntegrator, callback::DiscontinuityCallback,
                            prev_sign::Int, next_sign::Int, Θs, f)

Find time of propagated discontinuity for a certain dependent delay and previous
discontinuity, which is root of the function `f`, in the current time interval of
`integrator` with interpolation points `Θs`. Hereby `f` shows signs `prev_sign` and
`next_sign` at both ends of the time interval.
"""
function find_discontinuity_time(integrator::DDEIntegrator, callback::DiscontinuityCallback,
                                 prev_sign::Int, next_sign::Int, Θs, f)
    if prev_sign * next_sign < 0
        # if signs are different we already know that a discontinuity exists
        contains_discontinuity = true
        interp_index = callback.interp_points
        prev_sign_index = 1
    elseif callback.interp_points != 0
        # recheck interpolation intervals if no discontinuity found yet
        contains_discontinuity, interp_index, prev_sign_index =
            determine_discontinuity_existence(prev_sign, next_sign, Θs, f)
    end

    if contains_discontinuity
        if interp_index != 0
            top_Θ = Θs[interp_index] # Top at the smallest
            bottom_θ = Θs[prev_sign_index]
        else
            # without interpolation points
            top_Θ = one(integrator.t)
            bottom_θ = typeof(integrator.t)(0)
        end

        if f(top_Θ) == 0
          Θ = top_Θ
        else
          Θ = prevfloat(find_zero(f,(bottom_θ,top_Θ),Roots.AlefeldPotraShi(),atol = callback.abstol/100))
        end
        #Θ = prevfloat(...)
        # prevfloat guerentees that the new time is either 1 floating point
        # numbers just before the event or directly at zero, but not after.
        # If there's a barrier which is never supposed to be crossed,
        # then this will ensure that
        # The item never leaves the domain. Otherwise Roots.jl can return
        # a float which is slightly after, making it out of the domain, causing
        # havoc.
        new_t = integrator.dt*Θ
    else
        # no discontinuity found for given combination of lag and previous discontinuity
        new_t = zero(integrator.t)
    end

    new_t,contains_discontinuity
end

"""
    determine_discontinuity_existence(prev_sign::Int, next_sign::Int, Θs, f)

Determine whether function `f` has a root in the interval [0, 1] by checking signs of `f`
at 0 and 1 (`prev_sign` and `next_sign`, respectively) and at interpolation points in `Θs`.

This corresponds to the existence of a propagated discontinuity.
"""
function determine_discontinuity_existence(prev_sign::Int, next_sign::Int, Θs, f)
    contains_discontinuity = false
    interp_index = 0
    prev_sign_index = 1

    # use interpolation intervals to check for discontinuities
    for i in 2:length(Θs)-1
        new_sign = cmp(f(Θs[i]), 0)

        if prev_sign == 0
            prev_sign = new_sign
            prev_sign_index = i
        elseif prev_sign * new_sign < 0
            contains_discontinuity = true
            interp_index = i
            break
        end
    end

    contains_discontinuity, interp_index, prev_sign_index
end

"""
    apply_callback!(integrator::DDEIntegrator, callback::DiscontinuityCallback, cb_time,
                    order)

Handle discontinuity of order `order` at time `integrator.tprev + cb_time` in the current
time interval of `integrator` that was found by `callback`. Cause the current step to fail,
and add the found discontinuity to both the heap of discontinuities and of time stops.
"""
function DiffEqBase.apply_callback!(integrator::DDEIntegrator, callback::DiscontinuityCallback,
                         cb_time, order)
    # do not accept current step
    integrator.t = integrator.tprev
    integrator.force_stepfail = true
    integrator.accept_step = false

    # add new discontinuity of correct order at time of callback
    d = Discontinuity(integrator.t + cb_time, order)
    push!(integrator.opts.d_discontinuities, d)
    push!(integrator.opts.tstops, d.t)

    # u is not modified and step not saved
    # nevertheless callback set `save_cb` to `true` to avoid undesired saving in
    # `handle_callbacks!`
    false, true
end


text:
using .HDF5

export SMAPstack, SMAPseries, SMAPdimMetadata, SMAParrayMetadata, SMAPstackMetadata

const SMAPMISSING = -9999.0
const SMAPGEODATA = "Geophysical_Data"
const SMAPCRS = ProjString("+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

"""
    SMAPdimMetadata(val::Union{Dict,NamedTuple})
    SMAPdimMetadata(pairs::Pair...) => SMAPdimMetadata{Dict}
    SMAPdimMetadata(; kw...) => SMAPdimMetadata{NamedTuple}

`DimMetadata` wrapper for `SMAParray` dimensions.
"""
struct SMAPdimMetadata{T} <: AbstractDimMetadata{T}
    val::T
end

"""
    SMAParrayMetadata(val::Union{Dict,NamedTuple})
    SMAParrayMetadata(pairs::Pair...) => SMAParrayMetadata{Dict}
    SMAParrayMetadata(; kw...) => SMAParrayMetadata{NamedTuple}

`ArrayMetadata` wrapper for `SMAParray`.
"""
struct SMAParrayMetadata{T} <: AbstractArrayMetadata{T}
    val::T
end

"""
    SMAPstackMetadata(val::Union{Dict,NamedTuple})
    SMAPstackMetadata(pairs::Pair...) => SMAPstackMetadata{Dict}
    SMAPstackMetadata(; kw...) => SMAPstackMetadata{NamedTuple}

`StackMetadata` wrapper for `SMAPstack`.
"""
struct SMAPstackMetadata{T} <: AbstractStackMetadata{T}
    val::T
end

# Stack ########################################################################

"""
    SMAPstack(filename::String;
              dims=nothing,
              refdims=nothing,
              window=())

`AbstractGeoStack` for [SMAP](https://smap.jpl.nasa.gov/) datasets.

The simplicity of the format means `dims` and `metadata` are the same for all stack layers,
so we store them as stack fields.

## Arguments

- `filename`: `String` path to a SMAP .h5 file.

## Keyword arguments

- `dims`: Dimensions held on the stack as all layers have identical `Dimension`s.
  These are loaded from the HDF5 by default, but can be passed in to improve performance,
  as is done by [`SMAPseries`](@ref),
- `refdims`: As for `dims`. Often the position time `Dimension` from the `SMAPseries`.
- `metadata`: [`SMAPstackMetadata`](@ref) object. As for `dims`.
- `window`: Like `view` but lazy, for disk based data. Can be a tuple of Dimensions,
  selectors or regular indices. These will be applied when the data is loaded or indexed into.
"""
struct SMAPstack{T,D,R,W,M} <: DiskGeoStack{T}
    filename::T
    dims::D
    refdims::R
    window::W
    metadata::M
end
SMAPstack(filename::String;
          dims=smapread(smapdims, filename),
          refdims=(smap_timedim(smap_timefrompath(filename)),),
          window=(),
          metadata=smapread(smapmetadata, filename),
         ) =
    SMAPstack(filename, dims, refdims, window, metadata)

# Dummy SMAParray type. Actually generates a GeoArray
struct SMAParray end

# SMAP has fixed dims for all layers, so we store them on the stack.
dims(stack::SMAPstack, key::Key...) = stack.dims
dims(stack::SMAPstack, dim) = dims(dims(stack), dim)
refdims(stack::SMAPstack) = stack.refdims
metadata(stack::SMAPstack) = stack.metadata
missingval(stack::SMAPstack, key::Key...) = SMAPMISSING
childtype(stack::SMAPstack) = SMAParray
childkwargs(stack::SMAPstack) = ()

withsource(f, ::Type{SMAParray}, path::AbstractString, key...) =
    smapread(f, path)
withsourcedata(f, ::Type{SMAParray}, path::AbstractString, key) =
    smapread(path) do d
        f(d[smappath(string(key))])
    end

# Base methods

# Override getindex as we already have `dims` - they are
# fixed for the whole stack
Base.getindex(s::SMAPstack, key::Key, i1::Integer, I::Integer...) =
    smapread(filename(s)) do file
        dataset = file[smappath(key)]
        _window = maybewindow2indices(dataset, _dims, window(s))
        readwindowed(dataset, _window, I...)
    end
Base.getindex(s::SMAPstack, key::Key, I::Union{Colon,Integer,AbstractArray}...) =
    smapread(filename(s)) do file
        dataset = file[smappath(key)]
        _dims = dims(s)
        _window = maybewindow2indices(dataset, _dims, window(s))
        _dims, _refdims = slicedims(slicedims(_dims, refdims(s), _window)..., I)
        A = readwindowed(dataset, _window, I...)
        GeoArray(A, _dims, _refdims, Symbol(key), metadata(s), missingval(s))
    end

# HDF5 uses `names` instead of `keys` so we have to special-case it
Base.keys(stack::SMAPstack) =
    smapread(filename(stack)) do dataset
        cleankeys(names(dataset[SMAPGEODATA]))
    end

# Series #######################################################################

"""
    SMAPseries(filenames::AbstractString; kwargs...)
    SMAPseries(filenames::Vector{<:AbstractString}, dims=nothing; kwargs...)

[`GeoSeries`](@ref) loader for SMAP files and whole folders of files,
organised along the time dimension. Returns a [`GeoSeries`](@ref).

## Arguments

- `filenames`: A `String` path to a directory of SMAP files,
  or a vector of `String` paths to specific files.
- `dims`: `Tuple` containing `Ti` dimension for the series.
  Automatically generated form `filenames` unless passed in.

## Keyword Arguments

- `kwargs`: Passed to `GeoSeries`.
"""
SMAPseries(dir::AbstractString; kwargs...) =
    SMAPseries(joinpath.(dir, filter_ext(dir, ".h5")); kwargs...)
SMAPseries(filenames::Vector{<:AbstractString}, dims=nothing; kwargs...) = begin
    if dims isa Nothing
        usedpaths = String[]
        timeseries = []
        errors = []
        for filename in filenames
            println("Loading SMAP file: $filename")
            try
                t = smap_timefrompath(filename)
                push!(timeseries, t)
                push!(usedpaths, filename)
            catch e
                push!(errors, e)
            end
        end
        # Use the first files time dim as a template, but join vals into an array of times.
        timedim = smap_timedim(timeseries)
    else
        usedpaths = filenames
    end
    # Show errors after all files load, or you can't see them:
    if length(errors) > 0
        println("Some errors thrown during file load: ")
        println.(errors)
    end
    # Get the dims once for the whole series
    childkwargs = (
        dims=smapread(smapdims, first(filenames)),
        metadata=smapread(smapmetadata, first(filenames)),
    )
    GeoSeries(usedpaths, (timedim,); childtype=SMAPstack, childkwargs=childkwargs, kwargs...)
end

Base.:*(hrs::Int, ::Type{T}) where T<:Period = T(hrs)

# Utils ########################################################################

smapread(f, filepath::AbstractString) = h5open(f, filepath)

readwindowed(A::HDF5Dataset, window::Tuple{}) = HDF5.read(A)

smappath(key::Key) = SMAPGEODATA * "/" * string(key)

smap_timefrompath(path::String) = begin
    dateformat = DateFormat("yyyymmddTHHMMSS")
    dateregex = r"SMAP_L4_SM_gph_(\d+T\d+)_"
    datematch = match(dateregex, path)
    if !(datematch === nothing)
        DateTime(datematch.captures[1], dateformat)
    else
        error("Date/time not correctly formatted in path: $path")
    end
end

smap_timedim(t::DateTime) = smap_timedim(t:Hour(3):t)
smap_timedim(times::AbstractVector) =
    Ti(times, mode=Sampled(Ordered(), Regular(Hour(3)), Intervals(Start())))

# TODO actually add metadata to the dict
smapmetadata(dataset::HDF5.HDF5File) = SMAPstackMetadata(Dict())

smapdims(dataset::HDF5.HDF5File) = begin
    proj = read(attrs(root(dataset)["EASE2_global_projection"]), "grid_mapping_name")
    if proj == "lambert_cylindrical_equal_area"
        # There are matrices for lookup but all rows/colums are identical.
        # For performance and simplicity we just take a vector slice for each dim.
        extent = attrs(root(dataset)["Metadata/Extent"])
        lonbounds = read(extent["westBoundLongitude"]), read(extent["eastBoundLongitude"])
        latbounds = read(extent["southBoundLatitude"]), read(extent["northBoundLatitude"])
        latvec = read(root(dataset)["cell_lat"])[1, :]
        lonvec = read(root(dataset)["cell_lon"])[:, 1]
        lonmode = Mapped(
           order=Ordered(),
           span=Irregular(lonbounds),
           sampling=Intervals(Center()),
           crs=SMAPCRS,
           mappedcrs=EPSG(4326)
        )
        latmode = Mapped(
            order=Ordered(ReverseIndex(), ReverseArray(), ForwardRelation()),
            span=Irregular(latbounds),
            sampling=Intervals(Center()),
            crs=SMAPCRS,
            mappedcrs=EPSG(4326),
        )
        (Lon(lonvec; mode=lonmode), Lat(latvec; mode=latmode))
    else
        error("projection $proj not supported")
    end
end


text:
# Adapted from:
# SOSDEMO5 --- Upper bound for the structured singular value mu
# Section 3.5 of SOSTOOLS User's Manual

using SparseArrays

function sosdemo5_test(optimizer, config::MOIT.TestConfig, feasible::Bool, γ)
    @polyvar x[1:8]

    # The matrix under consideration
    α = 3 + sqrt(3);
    β = sqrt(3) - 1;
    a = sqrt(2/α);
    b = 1/sqrt(α);
    c = b;
    d = -sqrt(β/α);
    f = (1 + im)*sqrt(1/(α*β));
    U = [a 0; b b; c im*c; d f];
    V = [0 a; b -b; c -im*c; -im*f -d];
    M = U*V';

    Z = monomials(x, 1)

    function build_A(i)
        H = M[i,:]*M[i,:]' - (γ^2)*sparse([i],[i],[1],4,4)
        H = [real(H) -imag(H); imag(H) real(H)]
        dot(Z, H*Z)
    end
    A = build_A.(1:4)

    model = _model(optimizer)
    PolyJuMP.setpolymodule!(model, SumOfSquares)

    # -- Q(x)'s -- : sums of squares
    # Monomial vector: [x1; ... x8]
    Q = Vector{GramMatrix{JuMP.VariableRef, monomialtype(x[1]),
                          monovectype(x[1])}}(undef, 4)
    @variable model Q[1:4] SOSPoly(Z)

    # -- r's -- : constant sum of squares
    Z = monomials(x, 0)
    #r = Matrix{GramMatrix{JuMP.VariableRef}}(4,4) # FIXME doesn't work with 1x1 SDP matrix :(
    @variable model r[i=1:4,j=(i+1):4] >= 0

    # Constraint : -sum(Qi(x)*Ai(x)) - sum(rij*Ai(x)*Aj(x)) + I(x) >= 0
    expr = 0
    # Adding term
    for i = 1:4
        expr -= A[i]*Q[i]
    end
    for i = 1:4
        for j = (i+1):4
            expr -= A[i]*A[j]*r[i,j]
        end
    end
    # Constant term: I(x) = -(x1^4 + ... + x8^4)
    I = -sum(x.^4)
    expr = expr + I

    @constraint model expr >= 0

    JuMP.optimize!(model)

    # Program is feasible, thus 0.8724 is an upper bound for mu.
    if feasible
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    else
        @test JuMP.dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE
    end
end
const γ_opt = 0.8723
sosdemo5_infeasible_test(optimizer, config, ε = 1e-2) = sosdemo5_test(optimizer, config, false, γ_opt - ε)
sosdemo5_feasible_test(optimizer, config, ε = 1e-2) = sosdemo5_test(optimizer, config, true, γ_opt + ε)
sd_tests["sosdemo5_infeasible"] = sosdemo5_infeasible_test
sd_tests["sosdemo5_feasible"] = sosdemo5_feasible_test


text:
############################################################
## joAbstractLinearOperatorInplace - extra functions

# elements(jo)

# hasinverse(jo)

# issquare(jo)

# istall(jo)

# iswide(jo)

# iscomplex(jo)

# islinear(jo)

# isadjoint(jo)


text:
const SpatialDim = Union{XDim,YDim,ZDim}

"""
    Band <: Dimension

    Band(val=:)

Band [`Dimension`]($DDdimdocs) for multi-band rasters.

## Example:
```julia
banddim = Band(10:10:100)
# Or
val = A[Band(1)]
# Or
mean(A; dims=Band)
```
"""
@dim Band

"""
    mappedbounds(x)

Get the bounds converted to the [`mappedcrs`](@ref) value.

Whithout ArchGDAL loaded, this is just the regular bounds.
"""
function mappedbounds end

mappedbounds(dims::Tuple) = map(mappedbounds, dims)
mappedbounds(dim::Dimension) = bounds(dim)
mappedbounds(dim::Union{Y,X}) = mappedbounds(mode(dim), dim)
mappedbounds(::Mapped, dim) = bounds(dim)
@noinline mappedbounds(mode::IndexMode, dim) =
    if mode isa Projected
        error("Load ArchGDAL to convert Projected mode bounds to mapped")
    else
        error("cannot get mapped bounds of a $(nameof(typeof(mode))) mode dim")
    end

projectedbounds(dims::Tuple) = map(projectedbounds, dims)
projectedbounds(dim::Dimension) = bounds(dim)
projectedbounds(dim::Union{Y,X}) = projectedbounds(mode(dim), dim)
projectedbounds(::Projected, dim) = bounds(dim)
@noinline projectedbounds(mode::IndexMode, dim) =
    if mode isa Mapped
        error("Load ArchGDAL to convert Mapped mode dim to projected")
    else
        error("cannot get projected bounds of a $(nameof(typeof(mode))) mode dim")
    end

"""
    mappedindex(x)

Get the index value of a dimension converted to the `mappedcrs` value.

Whithout ArchGDAL loaded, this is just the regular dim value.
"""
function mappedindex end

mappedindex(dims::Tuple) = map(mappedindex, dims)
mappedindex(dim::Dimension) = index(dim)
mappedindex(dim::Union{Y,X}) = mappedindex(mode(dim), dim)
mappedindex(::Mapped, dim) = index(dim)
@noinline mappedindex(mode::IndexMode, dim) =
    if mode isa Projected
        error("Load ArchGDAL to convert Projected mode index to mapped")
    else
        error("cannot get mapped index of a $(nameof(typeof(mode))) mode dim")
    end

projectedindex(dims::Tuple) = map(projectedindex, dims)
projectedindex(dim::Dimension) = index(dim)
projectedindex(dim::Union{Y,X}) = projectedindex(mode(dim), dim)
projectedindex(::Projected, dim) = index(dim)
@noinline projectedindex(mode::IndexMode, dim) =
    if mode isa Mapped
        error("Load ArchGDAL to convert Mapped mode index to projected")
    else
        error("cannot get projected index of a $(nameof(typeof(mode))) mode dim")
    end


text:
using DelimitedFiles
include("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\Materials.jl")
include("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\read_ace.jl")

# Set up grid index mapping
E_min = 1.0E-10
logE_min = log10(E_min)
E_max = 1.0E+02
logE_max = log10(E_max)
log_span = logE_max - logE_min

n_grid = 10001 # Number of grid mapping spaces to use

function grid_index(E::Number)
    #=
    Returns index in grid index map for the given energy.

    :param E: Energy of interest
    =#
    f = (log10(E) - logE_min) / log_span # log fraction through the map
    return convert(Int64, floor(f * n_grid)) + 1
end

function grid_energy(idx::Int)
    #=
    Returns energy at grid index of interest.

    :param idx: Index of interest
    =#
    f = (idx - 1) / n_grid # log fraction through the map
    return exp10(logE_min + f*log_span)
end

# Read in table of cross section directories.
xsdir = readdlm("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\Lib80x\\xsdir.txt", ',')

# Define materials
#=
water = define_Material(
        density = 1.0,
        nuclides = ["O-nat", "H-nat"],
        atomic_fracs = [1.0, 2.0])

UO2 = define_Material(
        density = 10.0,
        nuclides = ["O-nat", "U-234"],
        atomic_fracs = [2.0, 1.0])

=#

B4C = define_Material(
        density = 2.52,
        nuclides = ["B-nat", "C-nat"],
        atomic_fracs = [4.0, 1.0])

Gd2O3 = define_Material(
        density = 7.40,
        nuclides = ["Gd-nat", "O-nat"],
        atomic_fracs = [2.0, 3.0])

#materials = Material[water, UO2]
materials = Material[B4C, Gd2O3]

# Create a list of incident neutron data needed
xs_nuclides = Tuple{Int64, Float64}[] # Pairs of ZA and T.
for i = 1:length(materials) # loop over materials
    for j = 1:length(materials[i].nuclides) # loop over nulcides in material
        ZA_T = (materials[i].nuclides[j], materials[i].temperature)
        # Add the nulcide identifier if it's not in the list yet
        if !(ZA_T in xs_nuclides)
            push!(xs_nuclides, ZA_T)
        end
    end
end

# Extract the paths to the data files needed. Keep nuclide ordering consistent
xs_nuclides_new = Tuple{Int64, Float64}[] # new ordering of xs_nuclides
xs_paths = String[] # paths to data files
for i = 1:length(xsdir[:,1])
    ZA = xsdir[i, 1]
    T = xsdir[i, 2]
    ZA_T = (ZA, T)
    if ZA_T in xs_nuclides
        push!(xs_nuclides_new, ZA_T)
        push!(xs_paths, xsdir[i, 3])
    end
end
xs_nuclides = xs_nuclides_new

# Load the nuclear data and populate the energy grid map
xs_data = AceData[] # list containing nuclear data
grid_map = zeros(Int64, (length(xs_nuclides), n_grid)) # energy grid map
for i = 1:length(xs_nuclides)
    acedata = read_ace(xs_paths[i])
    push!(xs_data, acedata)
    println(xs_nuclides[i])
    nuclide_E_grid = E_grid(acedata)
    for j = 1:n_grid
        grid_E = grid_energy(j)
        grid_map[i, j] = findlast(E -> E < grid_E, nuclide_E_grid)
    end

end

function Σ_T_array(material::Material, E::Number)
    nuclide_idxs = findall(ZA_T -> (ZA_T[1] in material.nuclides) & (ZA_T[2] == material.temperature), xs_nuclides)
    map_lo = grid_index(E)
    idxs_lo = grid_map[nuclide_idxs, map_lo]
    idxs_hi = grid_map[nuclide_idxs, map_lo + 1] .+ 1
    idxs = find_idx.(xs_data[nuclide_idxs], E, idxs_lo, idxs_hi)
    interp_fracs = interp_frac_log.(xs_data[nuclide_idxs], E, idxs)
    σ_T0 = σ_T_idx.(xs_data[nuclide_idxs], idxs)
    σ_T1 = σ_T_idx.(xs_data[nuclide_idxs], idxs.+1)
    Σ_Ts = material.atomic_dens.*exp.(log.(σ_T0) + log.(σ_T1./σ_T0).*interp_fracs) * 1.0E-24
    return Σ_Ts
end

function Σ_a_array(material::Material, E::Number)
    nuclide_idxs = findall(ZA_T -> (ZA_T[1] in material.nuclides) & (ZA_T[2] == material.temperature), xs_nuclides)
    map_lo = grid_index(E)
    idxs_lo = grid_map[nuclide_idxs, map_lo]
    idxs_hi = grid_map[nuclide_idxs, map_lo + 1] .+ 1
    idxs = find_idx.(xs_data[nuclide_idxs], E, idxs_lo, idxs_hi)
    interp_fracs = interp_frac_log.(xs_data[nuclide_idxs], E, idxs)
    σ_a0 = σ_a_idx.(xs_data[nuclide_idxs], idxs)
    σ_a1 = σ_a_idx.(xs_data[nuclide_idxs], idxs.+1)
    Σ_as = material.atomic_dens.*exp.(log.(σ_a0) + log.(σ_a1./σ_a0).*interp_fracs) * 1.0E-24
    return Σ_as
end

#=
Energy = 10 .^ collect(-10:0.001:1)
Σ_T = sum.(Σ_T_array(water, Energy[i]) for i = 1:length(Energy))
using PyPlot
loglog(Energy, Σ_T)
=#

using PyPlot
Energy = 10 .^ collect(-10:0.001:1)
Σ_a = sum.(Σ_a_array(B4C, Energy[i]) for i = 1:length(Energy))
loglog(Energy, Σ_a)
Σ_a = sum.(Σ_a_array(Gd2O3, Energy[i]) for i = 1:length(Energy))
loglog(Energy, Σ_a)

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
"""
A uniform mesh in one space dimension of `N` cells
"""
struct Uniform1DFVMesh{N,T}
  Δx :: T
  cell_centers :: Vector{T}
  cell_faces :: Vector{T}
end

getncells(mesh::Uniform1DFVMesh{N}) where {N} = N
getnnodes(mesh::Uniform1DFVMesh{N}) where {N} = N+1
getnodecoords(mesh::Uniform1DFVMesh, node_idx::Int) = mesh.cell_faces[node_idx]
getdimension(mesh::Uniform1DFVMesh) = 1

Base.summary(::Uniform1DFVMesh{N}) where {N} = string("FV 1D Mesh of ",N," cells")

function Base.show(io::IO, A::Uniform1DFVMesh)
  println(io,summary(A))
  print(io,"cells: ")
  show(io,getncells(A))
  println(io)
  print(io,"Δx: ")
  show(io, A.Δx)
  println(io)
  println(io,"Interval: [",A.cell_faces[1],",",A.cell_faces[end],"]")
end

TreeViews.hastreeview(x::Uniform1DFVMesh) = true
function TreeViews.treelabel(io::IO,x::Uniform1DFVMesh,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Text(Base.summary(x)))
end

function Uniform1DFVMesh(N::Int,bounds::Vector{T}) where {T}
    xend = bounds[end]
    xinit = bounds[1]
    L = xend - xinit
    dx = L/N
    xx = [i*dx+dx/2+xinit for i in 0:(N-1)]
    faces = [xinit + dx*i for i in 0:N]
    Uniform1DFVMesh{N,T}(dx,xx,faces)
end

@inline cell_volume(mesh::Uniform1DFVMesh) = mesh.Δx
@inline cell_volume(mesh::Uniform1DFVMesh,cell::Int) = mesh.Δx
cell_volumes(mesh::Uniform1DFVMesh) = mesh.Δx * ones(mesh.N)
get_nodes_matrix(mesh::Uniform1DFVMesh) = mesh.cell_faces
cell_centers(mesh::Uniform1DFVMesh) = mesh.cell_centers
cell_facets(mesh::Uniform1DFVMesh) = mesh.cell_faces


text:
abstract type AbstractCircleOctant <: AbstractShape end

struct CircleOctant{I <: Integer} <: AbstractCircleOctant
    center::Point{I}
    radius::I
end

struct ThickCircleOctant{I <: Integer} <: AbstractCircleOctant
    center::Point{I}
    radius::I
    thickness::I
end

abstract type AbstractCircle <: AbstractShape end

struct OddCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
end

struct EvenCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
end

struct Circle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
end

struct OddFilledCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
end

struct EvenFilledCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
end

struct FilledCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
end

struct OddThickCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
    thickness::I
end

struct EvenThickCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
    thickness::I
end

struct ThickCircle{I <: Integer} <: AbstractCircle
    position::Point{I}
    diameter::I
    thickness::I
end

#####
##### AbstractCircleOctant
#####

is_valid(shape::AbstractCircleOctant) = shape.radius >= zero(shape.radius)

get_i_min(shape::AbstractCircleOctant) = shape.center.i
get_i_max(shape::AbstractCircleOctant) = shape.center.i + shape.radius

get_j_min(shape::AbstractCircleOctant) = shape.center.j
get_j_max(shape::AbstractCircleOctant) = shape.center.j + shape.radius

get_drawing_optimization_style(::AbstractCircleOctant) = CHECK_BOUNDS

#####
##### CircleOctant
#####

function draw!(f::F, image::AbstractMatrix, shape::CircleOctant, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    center = shape.center
    radius = shape.radius

    I = typeof(radius)

    i_center = center.i
    j_center = center.j

    i = radius
    j = zero(I)

    f(image, i_center + i, j_center + j, color)

    constant = convert(I, 3) - convert(I, 2) * radius * radius

    while i > j + one(I)
        d = convert(I, 2) * i * i + convert(I, 2) * j * j + convert(I, 4) * j - convert(I, 2) * i + constant

        j += one(I)

        if d > zero(I)
            i -= one(I)
        end

        f(image, i_center + i, j_center + j, color)
    end

    return nothing
end

#####
##### ThickCircleOctant
#####

is_valid(shape::ThickCircleOctant) = (shape.radius >= zero(shape.radius)) && (shape.thickness > zero(shape.thickness) && (shape.thickness <= shape.radius + one(shape.radius)))

function draw!(image::AbstractMatrix, shape::ThickCircleOctant, color)
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    center = shape.center
    radius = shape.radius
    thickness = shape.thickness

    if is_outbounds(image, shape)
        return nothing
    end

    if is_inbounds(image, shape)
        draw!(image, ThickCircleOctant(center, radius, thickness), color) do image, i1, j1, i2, j2, color
            draw!(put_pixel_inbounds!, image, VerticalLine(i1, i2, j1), color)
        end
    else
        draw!(image, ThickCircleOctant(center, radius, thickness), color) do image, i1, j1, i2, j2, color
            draw!(put_pixel!, image, VerticalLine(i1, i2, j1), color)
        end
    end

    return nothing
end

function draw!(f::F, image::AbstractMatrix, shape::ThickCircleOctant, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    center = shape.center
    radius = shape.radius
    thickness = shape.thickness

    I = typeof(radius)

    i_center = center.i
    j_center = center.j

    radius_inner = radius - thickness + one(I)
    radius_outer = radius

    i_inner = radius_inner
    j_inner = zero(I)

    i_outer = radius_outer
    j_outer = zero(I)

    constant_inner = convert(I, 3) - convert(I, 2) * radius_inner * radius_inner
    constant_outer = convert(I, 3) - convert(I, 2) * radius_outer * radius_outer

    f(image, i_center + i_inner, j_center + j_inner, i_center + i_outer, j_center + j_outer, color)

    while i_inner > j_inner + one(I)
        d_inner = convert(I, 2) * i_inner * i_inner + convert(I, 2) * j_inner * j_inner + convert(I, 4) * j_inner - convert(I, 2) * i_inner + constant_inner
        d_outer = convert(I, 2) * i_outer * i_outer + convert(I, 2) * j_outer * j_outer + convert(I, 4) * j_outer - convert(I, 2) * i_outer + constant_outer

        j_inner += one(I)
        j_outer += one(I)

        if d_inner > zero(I)
            i_inner -= one(I)
        end

        if d_outer > zero(I)
            i_outer -= one(I)
        end

        f(image, i_center + i_inner, j_center + j_inner, i_center + i_outer, j_center + j_outer, color)
    end

    while i_outer > j_outer + one(I)
        d_outer = convert(I, 2) * i_outer * i_outer + convert(I, 2) * j_outer * j_outer + convert(I, 4) * j_outer - convert(I, 2) * i_outer + constant_outer

        j_outer += one(I)

        if d_outer > zero(I)
            i_outer -= one(I)
        end

        f(image, i_center + j_outer, j_center + j_outer, i_center + i_outer, j_center + j_outer, color)
    end

    return nothing
end

#####
##### AbstractCircle
#####

get_radius(shape::AbstractCircle) = shape.diameter ÷ oftype(shape.diameter, 2)

function get_center(shape::AbstractCircle, radius)
    i_position = shape.position.i
    j_position = shape.position.j
    return Point(i_position + radius, j_position + radius)
end

get_center(shape::AbstractCircle) = get_center(shape, get_radius(shape))

function get_center_radius(shape::AbstractCircle)
    radius = get_radius(shape)
    center = get_center(shape, radius)
    return (center, radius)
end

is_valid(shape::AbstractCircle) = shape.diameter > zero(shape.diameter)

get_i_min(shape::AbstractCircle) = shape.position.i
get_i_max(shape::AbstractCircle) = shape.position.i + shape.diameter - one(shape.diameter)

get_j_min(shape::AbstractCircle) = shape.position.j
get_j_max(shape::AbstractCircle) = shape.position.j + shape.diameter - one(shape.diameter)

get_drawing_optimization_style(::AbstractCircle) = CHECK_BOUNDS

#####
##### OddCircle
#####

is_valid(shape::OddCircle) = isodd(shape.diameter) && shape.diameter > zero(shape.diameter)

function draw!(f::F, image::AbstractMatrix, shape::OddCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    center, radius = get_center_radius(shape)

    draw!(image, CircleOctant(center, radius), color) do image, i, j, color
        draw!(f, image, OddOctantSymmetricShape(center, Point(i, j)), color)
    end

    return nothing
end

#####
##### EvenCircle
#####

is_valid(shape::EvenCircle) = iseven(shape.diameter) && shape.diameter > zero(shape.diameter)

function draw!(f::F, image::AbstractMatrix, shape::EvenCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    center, radius = get_center_radius(shape)

    draw!(image, CircleOctant(center, radius), color) do image, i, j, color
        draw!(f, image, EvenOctantSymmetricShape(center, Point(i, j)), color)
    end

    return nothing
end

#####
##### Circle
#####

function draw!(image::AbstractMatrix, shape::Circle, color)
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    diameter = shape.diameter

    if iseven(diameter)
        draw!(image, EvenCircle(position, diameter), color)
    else
        draw!(image, OddCircle(position, diameter), color)
    end

    return nothing
end

function draw!(f::F, image::AbstractMatrix, shape::Circle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    diameter = shape.diameter

    if iseven(diameter)
        draw!(f, image, EvenCircle(position, diameter), color)
    else
        draw!(f, image, OddCircle(position, diameter), color)
    end

    return nothing
end

#####
##### OddFilledCircle
#####

is_valid(shape::OddFilledCircle) = is_valid(OddCircle(shape.position, shape.diameter))

function draw!(f::F, image::AbstractMatrix, shape::OddFilledCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    center, radius = get_center_radius(shape)

    draw!(image, CircleOctant(center, radius), color) do image, i, j, color
        draw!(f, image, OddSymmetricVerticalLines4(center, Point(i, j)), color)
    end

    return nothing
end

#####
##### EvenFilledCircle
#####

is_valid(shape::EvenFilledCircle) = is_valid(EvenCircle(shape.position, shape.diameter))

function draw!(f::F, image::AbstractMatrix, shape::EvenFilledCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    diameter = shape.diameter

    I = typeof(diameter)

    i_position = position.i
    j_position = position.j

    center, radius = get_center_radius(shape)

    if diameter == convert(I, 2)
        draw!(f, image, FilledRectangle(position, convert(I, 2), convert(I, 2)), color)
    end

    draw!(image, CircleOctant(center, radius), color) do image, i, j, color
        if i > center.i && j > center.j
            draw!(f, image, EvenSymmetricVerticalLines4(center, Point(i, j)), color)
        else
            return nothing
        end
    end

    return nothing
end

#####
##### FilledCircle
#####

function draw!(image::AbstractMatrix, shape::FilledCircle, color)
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    diameter = shape.diameter

    if iseven(diameter)
        draw!(image, EvenFilledCircle(position, diameter), color)
    else
        draw!(image, OddFilledCircle(position, diameter), color)
    end

    return nothing
end

function draw!(f::F, image::AbstractMatrix, shape::FilledCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    diameter = shape.diameter

    if iseven(diameter)
        draw!(f, image, EvenFilledCircle(position, diameter), color)
    else
        draw!(f, image, OddFilledCircle(position, diameter), color)
    end

    return nothing
end

#####
##### OddThickCircle
#####

function is_valid(shape::OddThickCircle)
    diameter = shape.diameter
    thickness = shape.thickness

    I = typeof(shape.diameter)

    return isodd(diameter) && diameter > zero(I) && thickness > zero(I) && convert(I, 2) * thickness <= diameter + one(I)
end

function draw!(f::F, image::AbstractMatrix, shape::OddThickCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    thickness = shape.thickness

    center, radius = get_center_radius(shape)

    draw!(image, ThickCircleOctant(center, radius, thickness), color) do image, i1, j1, i2, j2, color
        draw!(f, image, OddOctantSymmetricShape(center, VerticalLine(i1, i2, j1)), color)
    end

    return nothing
end

#####
##### EvenThickCircle
#####

function is_valid(shape::EvenThickCircle)
    diameter = shape.diameter
    thickness = shape.thickness

    I = typeof(shape.diameter)

    return iseven(diameter) && diameter > zero(I) && thickness > zero(I) && convert(I, 2) * thickness <= diameter
end

function draw!(f::F, image::AbstractMatrix, shape::EvenThickCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    thickness = shape.thickness

    center, radius = get_center_radius(shape)

    draw!(image, ThickCircleOctant(center, radius, thickness), color) do image, i1, j1, i2, j2, color
        draw!(f, image, EvenOctantSymmetricShape(center, VerticalLine(i1, i2, j1)), color)
    end

    return nothing
end

#####
##### ThickCircle
#####

function is_valid(shape::ThickCircle)
    position = shape.position
    diameter = shape.diameter
    thickness = shape.thickness

    if iseven(diameter)
        return is_valid(EvenThickCircle(position, diameter, thickness))
    else
        return is_valid(OddThickCircle(position, diameter, thickness))
    end
end

function draw!(image::AbstractMatrix, shape::ThickCircle, color)
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    diameter = shape.diameter
    thickness = shape.thickness

    if iseven(diameter)
        draw!(image, EvenThickCircle(position, diameter, thickness), color)
    else
        draw!(image, OddThickCircle(position, diameter, thickness), color)
    end

    return nothing
end

function draw!(f::F, image::AbstractMatrix, shape::ThickCircle, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    diameter = shape.diameter
    thickness = shape.thickness

    if iseven(diameter)
        draw!(f, image, EvenThickCircle(position, diameter, thickness), color)
    else
        draw!(f, image, OddThickCircle(position, diameter, thickness), color)
    end

    return nothing
end


text:
module Day01

using Chain
using ShiftedArrays

import ..data_dir # from parent module

export part1, part2

const data_file = joinpath(data_dir, "day01")
const input = readlines(data_file) # read into an array, each line is an element

# PART 1
function part1(input = input)
    @chain input begin
        parse.(Int32, _)  # convert strings to integers
        _ .- lag(_) # compute difference with previous value
        _[2:end] # remove first element, which will be missing
        filter(x -> x > 0, _) # filter to the increases in value
        length # now count
    end
end

# PART 2
function part2(input = input)
    @chain input begin
        parse.(Int32, _) # convert strings to integers
        _ + lag(_) + lag(_, 2) # compute the sum of each value with its previous 2
        _ - lag(_) # now compute the difference with the previous element
        _[4:end] # remove first three elements, since we're summing triples
        filter(x -> x > 0, _) # filter to positive values, which are the increases in value
        length # now count
    end
end

end # module


text:
import Pkg; Pkg.activate("../../")

using StatsFuns
using BenchmarkTools
import Random
using Flux
using Zygote
include(joinpath(@__DIR__, "BnpUtil_old.jl"))

function log_sum_exp(logx::T) where T <: AbstractArray
  mx = maximum(logx)
  return mx + log(sum(exp.(logx .- mx)))
end

Random.seed!(0);
x = randn(10)
@assert StatsFuns.logsumexp(x) == BnpUtil.logsumexpdd(x)[1]

# Logsumexp benchmarks
@benchmark StatsFuns.logsumexp(x)
@benchmark BnpUtil.logsumexp(x)[1]
@benchmark log_sum_exp(x)

### Gmm benchmarks ###
Random.seed!(0);
N, K = 101, 10
X = randn(N, K)
m = randn(K)
s = rand(K)
w = let
    _w = rand(K)
    _w / sum(_w)
end

lpdf_gmm_sf(x, m, s, w) = sum(logsumexp(normlogpdf.(m, s, x) .+ log.(w), dims=2))
@assert sum(BnpUtil.lpdf_gmm(X, m', s', w', dims=2)) == lpdf_gmm_sf(X, m', s', w')
@btime lpdf_gmm_sf(X, m', s', w')
@btime sum(BnpUtil.lpdf_gmm(X, m', s', w', dims=2))

# NOTE:
# logsumexp and lpdf_gmm seems to be comparable. In fact, for forward
# computaton, StatsFuns seems to be a little faster for logsumexp.

### BACKWARD: ###
# BnpUtil
@benchmark gs = Flux.gradient(Flux.params(m)) do
    sum(BnpUtil.lpdf_gmm(X, m', s', w', dims=2))
end

# StatsFuns (error?!)
@benchmark gs = Flux.gradient(Flux.params(m)) do
    lpdf_gmm_sf(X, m', s', w')
end

# Diagnose problem
x, loc, scale = 0.0, 2.0, 1.0
normlogpdf(loc, scale, x)
Flux.gradient(mu -> normlogpdf(mu, scale, x), loc)

function my_normlogpdf(loc, scale, x)
    z = (x - loc) / scale
    return -z * z * 0.5 - 0.5 * log(2 * pi * scale*scale)
end
Flux.gradient(mu -> my_normlogpdf(mu, scale, x), loc)

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
# elimat(m) returns the elimination matrix Lm
# The elimination matrix Lm is for any matrix F, Vech(F)=Lm Vec(F)
function elimat(m::Int64)
    A = eye(m^2)
    L = A[1:m,:]
    for n in 2:m
        S = A[m*(n-1)+1:n*m,:]
        S = S[n:end,:]
        L = [L;S]
    end
    return L
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockWeighter(sides)
    BlockWeighter(side₁, side₂, ...)

A weighting method that assigns weights to points in spatial object
based on blocks of given `sides`. The number `n` of points inside a
block determines the weights `1/n` of these points.
"""
struct BlockWeighter{T,N} <: AbstractWeighter
  sides::SVector{N,T}
end

BlockWeighter(sides::NTuple{N,T}) where {N,T} =
  BlockWeighter{T,N}(sides)

BlockWeighter(sides::Vararg{T,N}) where {N,T} =
  BlockWeighter(sides)

function weight(sdomain::AbstractDomain, weighter::BlockWeighter)
  p = partition(sdomain, BlockPartitioner(weighter.sides))

  weights = Vector{Float64}(undef, npoints(sdomain))
  for s in subsets(p)
    n = length(s)
    weights[s] .= 1/n
  end

  SpatialWeights(sdomain, weights)
end

weight(sdata::AbstractData, weighter::BlockWeighter) =
  weight(domain(sdata), weighter)


text:
using Statistics

function sim(pomdp::POMDP, policy::ModelCheckingPolicy, i::Int64)
    up = DiscreteUpdater(pomdp)
    b0 = initialize_belief(up, initialstate_distribution(pomdp))
    hr = HistoryRecorder(max_steps=50, rng=MersenneTwister(i))
    hist = simulate(hr, pomdp, policy, up, b0)
    return hist
end

function many_sims(pomdp::POMDP, policy::ModelCheckingPolicy, n_sims::Int64 = 100)
    successes = zeros(Int64, n_sims)
    mu = zeros(n_sims - 1)
    sig = zeros(n_sims - 1)
    successes = @showprogress pmap(1:n_sims) do x
                            hist = sim(pomdp, policy, x)
                            return undiscounted_reward(hist) > 0.
                end

    # @showprogress for i=1:n_sims
    #     hist = sim(pomdp, policy)
    #     successes[i] = undiscounted_reward(hist) > 0.
    #     if i > 1
    #         mu[i-1], sig[i-1] = leave_one_out_mean_estimate(successes[1:i])
    #         println("mean: ", mu[i-1], " sig ", sig[i-1])
    #     end
    # end
    mu, sig = running_stats(successes)
    return successes, mu, sig
end

function running_stats(vec)
    mu = zeros(length(vec) - 1)
    sig = zeros(length(vec) - 1)
    for i=2:length(vec)
        mu[i-1], sig[i-1] = leave_one_out_mean_estimate(vec[1:i])
    end
    return mu, sig
end

function leave_one_out_mean_estimate(vec)
    means = zeros(length(vec))
    for i=1:length(vec)
        means[i] = mean(vec[j] for j=1:length(vec) if j != i)
    end
    return mean(means), std(means)/sqrt(length(vec))
end


text:
using Pkg.BinaryPlatforms, JSON
using VersionsJSONUtil
import VersionsJSONUtil: PortableWindows
using Test

const download_urls = Dict(
    v"1.7.0-beta3" => Dict(
        MacOS(:aarch64) =>              "https://julialang-s3.julialang.org/bin/mac/aarch64/1.7/julia-1.7.0-beta3-macaarch64.dmg",
    ),
    v"1.6.2" => Dict(
        Linux(:x86_64) =>               "https://julialang-s3.julialang.org/bin/linux/x64/1.6/julia-1.6.2-linux-x86_64.tar.gz",
        Linux(:i686) =>                 "https://julialang-s3.julialang.org/bin/linux/x86/1.6/julia-1.6.2-linux-i686.tar.gz",
        Linux(:aarch64) =>              "https://julialang-s3.julialang.org/bin/linux/aarch64/1.6/julia-1.6.2-linux-aarch64.tar.gz",
        Linux(:armv7l) =>               "https://julialang-s3.julialang.org/bin/linux/armv7l/1.6/julia-1.6.2-linux-armv7l.tar.gz",
        Linux(:powerpc64le) =>          "https://julialang-s3.julialang.org/bin/linux/ppc64le/1.6/julia-1.6.2-linux-ppc64le.tar.gz",
        Linux(:x86_64, libc = :musl) => "https://julialang-s3.julialang.org/bin/musl/x64/1.6/julia-1.6.2-musl-x86_64.tar.gz",
        MacOS(:x86_64) =>               "https://julialang-s3.julialang.org/bin/mac/x64/1.6/julia-1.6.2-mac64.dmg",
        Windows(:x86_64) =>             "https://julialang-s3.julialang.org/bin/winnt/x64/1.6/julia-1.6.2-win64.exe",
        Windows(:i686) =>               "https://julialang-s3.julialang.org/bin/winnt/x86/1.6/julia-1.6.2-win32.exe",
        PortableWindows(:x86_64) =>     "https://julialang-s3.julialang.org/bin/winnt/x64/1.6/julia-1.6.2-win64.zip",
        PortableWindows(:i686) =>       "https://julialang-s3.julialang.org/bin/winnt/x86/1.6/julia-1.6.2-win32.zip",
        FreeBSD(:x86_64) =>             "https://julialang-s3.julialang.org/bin/freebsd/x64/1.6/julia-1.6.2-freebsd-x86_64.tar.gz",
    ),
)

@testset "VersionsJSONUtil.jl" begin
    @testset "Download URLs for $v" for v in keys(download_urls)
        for (p, url) in download_urls[v]
            @test VersionsJSONUtil.download_url(v, p) == url
        end
    end
end


text:
# Functions to implement goal functions and testing for goal functions
using Printf
export Goal, GoalList, randgoal, randgoallist, rand_env_goallist, env_goal, goal_count, goal_check, my_test_goal
export randgoal_filtered, randgoallist_filtered

# Moved these to aliases.jl
#const Goal =  Vector{MyInt}
#const GoalList = Vector{Goal}

#=
goal1 = MyInt[ 0xb, 0xe ]
goal2 = MyInt[ 0x6, 0x9 ]
goal3 = MyInt[ 0x1, 0xe ]
goal4 = MyInt[ 0x1, 0xe ]
goal5 = MyInt[ 0x3, 0x9 ]
goallist2 = [goal1, goal2, goal3 ]

goallist3 = [0x50, 0x05, 0xf5, 0xa0, 0xe0, 0x0e, 0xf1, 0x80, 0x08, 0xfa]
=#

function component_type( numinputs::Int64 )
  if numinputs <= 2
    comp_type = collect(0x00:0x0f)
  elseif numinputs == 3
    comp_type = collect(0x00:0xff)
  elseif numinputs == 4
    comp_type = UInt16
  elseif numinputs == 5
    comp_type = UInt32
  elseif numinputs == 6
    comp_type = UInt64
  elseif numinputs == 7
    comp_type = UInt128
  elseif numinputs > 7
    error("randgoal doesn't work for numinputs > 7")
  end
  comp_type
end

# generate a random goal assuming that MyInt is set correctly for numinputs
# Creates an inexacterror if conversion to MyInt doesn't work
function randgoal(numinputs::Int64, numoutputs::Int64 )
  map(x->convert( MyInt, x), rand( component_type( numinputs ), numoutputs ) )
end

function randgoal(p::Parameters)
  randgoal(p.numinputs,p.numoutputs)
end

# generate a random goallist of length ngoals
function randgoallist(ngoals::Int64, numinputs::Int64, numoutputs::Int64; repetitions::Int64=1)
  result = Vector{MyInt}[]
  if ngoals % repetitions != 0
    error("ngoals=",ngoals," must be a multiple of repetitions=",repetitions," in randgoallist")
  end
  #println("  collect(1:repetitions.ngoals): ",collect(1:repetitions:ngoals))
  for i = 1:repetitions:ngoals
    goal = rand( component_type(numinputs), numoutputs )
    for j = 1:repetitions
      #println("(i,j): ",(i,j)," i+j-1: ",i+j-1)
      #result[i+j-1] = goal
      push!(result,goal)
    end
  end
  result
end

function randgoallist(ngoals::Int64, p::Parameters; repetitions::Int64=1)
  randgoallist( ngoals, p.numinputs, p.numoutputs, repetitions=repetitions )
end

# Returns a goal list of ngoals which is obtained by generating div(ngoals,repetitions)
#   goals, and then perturing these base goals repetitions-1 times to generate
#   the additional goals.
# If use_env_goals==true, then repeat goals are obtained by flipping specific bits
#    of the base goals
# If use_env_goals==false, then repeat goals are obtained by flipping random bits
#   the additional goals.
function rand_env_goallist( ngoals::Int64, numinputs::Int64, numoutputs::Int64,
      repetitions::Int64, num_flipped_bits::Int64; perturb_goal::Bool=true )
  #println("env goallist: perturb_goal: ",perturb_goal)
  result = Vector{MyInt}[]
  if numoutputs % repetitions != 0
    error("numoutputs=",numoutputs," must be a multiple of repetitions=",repetitions," in rand_env_goallist")
  end
  #println("f_bit_list: ",f_bits_list)
  for i = 1:ngoals
    #g = rand( component_type(numinputs), div( numoutputs, repetitions ) )
    g = randgoal( numinputs, div(numoutputs,repetitions) )
    if perturb_goal
      goal = perturb_bits_goal( g, num_flipped_bits, repetitions, numinputs )
    else
      f_bits_list = flipped_bits_list( num_flipped_bits, repetitions-1, numinputs )
      goal = env_goal( g, f_bits_list )
    end
    Base.push!(result,goal)
  end
  result
end

# return a random goal with a minimum count of min_count in the dataframe column cdf[count_field].
# The goal is filtered by count from the dataframe column cdf[count_field].
function randgoal_filtered( cdf::DataFrame, count_field::Symbol, min_count::Int64, numinputs::Int64, numoutputs::Int64 )
  max_iters = 20    # prevent an infinite loop
  goal = randgoal( numinputs, numoutputs )
  iter = 0
  count = cdf[cdf.goal.==@sprintf("0x%x",goal[1]),count_field][1]
  while iter < max_iters && count < min_count
    goal = randgoal( numinputs, numoutputs )
    println("goal: ",goal)
    count_array = cdf[cdf.goal.==@sprintf("0x%x",goal[1]),count_field]
    println("count_array: ",count_array)
    count = count_array[1]
    iter += 1
  end
  if iter < max_iters
    goal
  else
    error("too many iterations in randgoal filtered")
  end
end

# generate a random goallist of length ngoals
# The goals are filtered by count from the dataframe cdf field count_field.
function randgoallist_filtered( cdf::DataFrame, count_field::Symbol, min_count::Int64, ngoals::Int64, numinputs::Int64,
     numoutputs::Int64; repetitions::Int64=1)
  result = Vector{MyInt}[]
  if ngoals % repetitions != 0
    error("ngoals=",ngoals," must be a multiple of repetitions=",repetitions," in randgoallist")
  end
  #println("  collect(1:repetitions.ngoals): ",collect(1:repetitions:ngoals))
  for i = 1:repetitions:ngoals
    goal = randgoal_filtered( cdf, count_field, min_count, numinputs, numoutputs )
    for j = 1:repetitions
      #println("(i,j): ",(i,j)," i+j-1: ",i+j-1)
      #result[i+j-1] = goal
      push!(result,goal)
    end
  end
  result
end

# The bits that will be flipped to produce the additional goals in function env_goal()
# Example (where MyInt==UInt16):
# julia> flipped_bits_list( 2,3,3)
# 3-element Array{UInt16,1}:
#  0x00c0
#  0x0030
#  0x000c
function flipped_bits_list( num_flipped_bits::Int64, repetitions::Int64, numinputs::Int64 )
  one = MyInt(1)
  if num_flipped_bits > 0
    ones = one
    for i = 1:(num_flipped_bits-1)
      ones = ones<<1 | one
    end
  else
    ones = MyInt(0)   # No flips if num_flipped_bit == 0
  end
  fb = ones << (2^numinputs-num_flipped_bits)
  #Printf.@printf("flipped_bits_list: ones:  0x%04x  fb: 0x%04x\n",ones,fb)
  [ fb >> (num_flipped_bits*i) for i = 0:(repetitions-1) ]
end

# returns a new goal with each component of goal replicated length(flipped_bits_list) times,
#   where each replication has nbits_to_perturb bits of the original component perturbed
function env_goal( goal::Goal, flipped_bits_list::Vector{MyInt} )
  #println("env goal")
  repetitions = length(flipped_bits_list)
  new_goal = Goal()
  for c in goal
    Base.push!( new_goal, c )
    for i = 1:repetitions
      Base.push!( new_goal, xor( c, flipped_bits_list[i] ))
    end
  end
  #println("new_goal: ",new_goal)
  new_goal
end

# returns an new goal with each component of goal replicated repetitions times,
#   where each replication has nbits_to_perturb randomly chosen bits
#   of the original component perturbed
function perturb_bits_goal( goal::Goal, nbits_to_perturb::Int64, repetitions::Int64, numinputs::Int64 )
  #println("perturb goal")
  new_goal = Goal()
  for c in goal
    push!( new_goal, c )
    for i = 1:repetitions-1
      rp = rand_perturbation( 2^numinputs, nbits_to_perturb )
      #Printf.@printf("rp: 0x%04x\n",rp)
      #push!( new_goal, xor( c, rand_perturbation( 2^numinputs, nbits_to_perturb )))
      push!( new_goal, xor( c, rp))
   end
  end
  #println("ptb: Goal: ",goal)
  new_goal
end

# Returns a MyInt with nnbits_to_perturb 1 bits.
# When perturbing a goal, maxlen should be 2^numinputs
function rand_perturbation( maxlength::Int64, nbits_to_perturb::Int64 )
  bits_to_perturb = rand_set( maxlength, nbits_to_perturb )
  #println("bits_to_perturb: ",bits_to_perturb)
  result = MyInt(0)
  for i in bits_to_perturb
    result |= MyInt(1) << i
  end
  result
end

# returns a random set of n_elements intgers from the integer range [0,maxlen-1]
function rand_set( maxlen::Int64, n_elements::Int64 )
  @assert n_elements < maxlen
  S = Int64[]
  while length(unique!(S)) < n_elements
    push!(S,rand(0:(maxlen-1)))
  end
  S
end

# Returns true if testgoal is equal to any of the goals in goallist defined above
function goal_check( testgoal::Goal, goallist::GoalList )
  length( filter( x -> testgoal==x, goallist )) > 0
end

# Returns the number of components of testGoal that are equal to the corresponding component of goal.
# Never used
#function my_test_goal( testGoal::Goal, goal::Goal )
#  length( filter( x -> x == 0, testGoal .⊻ goal ))
#end

function my_test_goal( testGoal::Goal, goallist::GoalList )
  achieved = 0
  for g in goallist
    ga = length( filter( x -> x == 0, testGoal .⊻ g ))
    if ga > achieved
      achieved = ga
    end
    println("my_test_goal: achieved: ",achieved)
  end
  achieved
end

# Returns a 2-tuple of:
#    1.  the  maximum number of matching subgoals of testgoal to components of goallist
#    2.  the  index of the goal of goallist that gave this maximum
# If testgoal equals all subgoals of a goal of goallist, then returns numouputs and the index of the goal in goallist
# If no component of testgoal matches the corresponding component of any goal in goallist, return 0
function goal_count( testgoal::Goal, goallist::GoalList )
   findmax( map( x -> my_test_goal( x, testgoal ), goallist ))[1]
end

# Test if at least numsubgoals components of testgoal match a component of goal (in any order)
#=
function my_test_goal( testgoal::Goal, goal::Goal, numsubgoals::Integer )
  filter( x -> (Base.findfirst( y -> y == x, testgoal ) != nothing ), goal)
end
=#


text:
using BinDeps

@BinDeps.setup

pythia = library_dependency("pythia",aliases=["libpythia","libpythia8"])
fastjet = library_dependency("fastjet",aliases=["libfastjet"])

provides(Sources,URI("http://fastjet.fr/repo/fastjet-3.1.2.tar.gz"),fastjet)
sources = provides(Sources,URI("http://home.thep.lu.se/~torbjorn/pythia8/pythia8205.tgz"),pythia)
builddir = joinpath(BinDeps.builddir(pythia),pythia.name)
bindir = joinpath(builddir,"bin")
steps = BinDeps.generate_steps(pythia, Autotools(configure_options=["--enable-shared"]),Dict())
options = ["--enable-shared","--enable-debug"]
@osx_only push!(options,"--cxx=clang")
provides(BuildProcess,(@build_steps begin
        GetSources(pythia)
        CreateDirectory(builddir)
        # Copy the entire source directory here, since Pythia does not support out of tree builds
        `cp -rpf $(BinDeps.srcdir(pythia,BinDeps.gethelper(pythia,Sources)[1],Dict()))/ $builddir/`
        @build_steps begin
            ChangeDirectory(builddir)
            # Not an actual configure
            `./configure $options --prefix=$(BinDeps.usrdir(pythia))`
            MakeTargets()
            MakeTargets(["install"])
            `chmod +x $(BinDeps.bindir(pythia))/pythia8-config`
        end
    end),pythia)

provides(BuildProcess,Autotools(libtarget = "src/.libs/libfastjet.la"),fastjet)

@BinDeps.install Dict( :pythia => :libpythia, :fastjet => :libfastjet )