text:
import Base: setindex!, getindex, size, IndexStyle, sum, eltype
using Base: IndexLinear

using Transpiler.cli: LocalPointer
using Transpiler.cli: get_local_id, get_global_id, barrier,  CLK_LOCAL_MEM_FENCE
using Transpiler.cli: get_local_size, get_global_size, get_group_id, get_num_groups
import GPUArrays: synchronize, synchronize_threads, device, global_size, linear_index

import GPUArrays: LocalMemory
using GPUArrays: AbstractDeviceArray

"""
Array type on the device
"""
struct DeviceArray{T, N, Ptr} <: AbstractDeviceArray{T, N}
    size::NTuple{N, Cuint}
    ptr::Ptr
end
# shaninagans for uploading CLArrays to OpenCL as a DeviceArray
# (spoiler alert: they can't contain pointers while uploading, but can on the device)
"""
Dummy pointer type for inlining into structs that get uploaded to the GPU
"""
struct HostPtr{T}
    ptr::Int32
    (::Type{HostPtr{T}})() where T = new{T}(Int32(0))
end
eltype(::Type{HostPtr{T}}) where T = T
const PreDeviceArray{T, N} = DeviceArray{T, N, HostPtr{T}} # Pointer free variant for kernel upload
const GlobalArray{T, N} = DeviceArray{T, N, GlobalPointer{T}}
const LocalArray{T, N} = DeviceArray{T, N, LocalPointer{T}}

const OnDeviceArray{T, N} = Union{GlobalArray{T, N}, LocalArray{T, N}} # Variant on the device containing the correct pointer

size(x::DeviceArray) = x.size
size(x::DeviceArray, i::Integer) = x.size[i]

getindex(x::OnDeviceArray, ilin::Integer) =  x.ptr[ilin]
function getindex(x::OnDeviceArray{T, N}, i::Vararg{Integer, N}) where {T, N}
    ilin = gpu_sub2ind(size(x), Cuint.(i))
    return x.ptr[ilin]
end
function setindex!(x::OnDeviceArray{T, N}, val, ilin::Integer) where {T, N}
    x.ptr[ilin] = T(val)
    return
end

function setindex!(x::OnDeviceArray{T, N}, val, i::Vararg{Integer, N}) where {T, N}
    ilin = gpu_sub2ind(size(x), Cuint.(i))
    x.ptr[ilin] = T(val)
    return
end

predevice_type(::Type{T}) where T = T

kernel_convert(A::CLArray{T, N}) where {T, N} = PreDeviceArray{T, N}(A.size, HostPtr{T}())
predevice_type(::Type{GlobalArray{T, N}}) where {T, N} = PreDeviceArray{T, N}
device_type(::CLArray{T, N}) where {T, N} = GlobalArray{T, N}
reconstruct(x::PreDeviceArray{T, N}, ptr::GlobalPointer{T}) where {T, N} = GlobalArray{T, N}(x.size, ptr)

# some converts to inline CLArrays into tuples and refs
kernel_convert(x::RefValue{T}) where T <: CLArray = RefValue(kernel_convert(x[]))
predevice_type(::Type{RefValue{T}}) where T <: GlobalArray = RefValue{predevice_type(T)}
device_type(x::RefValue{T}) where T <: CLArray = RefValue{device_type(x[])}
reconstruct(x::RefValue{T}, ptr::GlobalPointer) where T <: PreDeviceArray = RefValue(reconstruct(x[], ptr))

kernel_convert(x::Tuple) = kernel_convert.(x)
predevice_type(::Type{T}) where T <: Tuple = Tuple{predevice_type.((T.parameters...))...}
device_type(x::T) where T <: Tuple = Tuple{device_type.(x)...}

@generated function reconstruct(x::Tuple, ptrs::GlobalPointer...)
    ptrlist = to_tuple(ptrs)
    tup = Expr(:tuple)
    ptr_idx = 0
    for (xi, T) in enumerate(to_tuple(x))
        hasptr, fields = contains_pointer(T)
        if hasptr
            # consume the n pointers that T contains
            ptr_args = ntuple(i-> :(ptrs[$(i + ptr_idx)]), length(fields))
            ptr_idx += 1
            push!(tup.args, :(reconstruct(x[$xi], $(ptr_args...))))
        else
            push!(tup.args, :(x[$xi]))
        end
    end
    return tup
end

#synchronize
function synchronize(x::CLArray)
    cl.finish(global_queue(x)) # TODO figure out the diverse ways of synchronization
end

immutable KernelState
    empty::Int32
    KernelState() = new(Int32(0))
end

for (i, sym) in enumerate((:x, :y, :z))
    for (f, fcl, isidx) in (
            (:blockidx, get_group_id, true),
            (:blockdim, get_local_size, false),
            (:threadidx, get_local_id, true),
            (:griddim, get_num_groups, false)
        )

        fname = Symbol(string(f, '_', sym))
        if isidx
            @eval GPUArrays.$fname(::KernelState)::Cuint = $fcl($(i-1)) + Cuint(1)
        else
            @eval GPUArrays.$fname(::KernelState)::Cuint = $fcl($(i-1))
        end
    end
end

global_size(state::KernelState) = get_global_size(0)
linear_index(state::KernelState) = get_global_id(0) + Cuint(1)

synchronize_threads(::KernelState) = cli.barrier(CLK_LOCAL_MEM_FENCE)

LocalMemory(state::KernelState, ::Type{T}, ::Val{N}, ::Val{C}) where {T, N, C} = Transpiler.cli.LocalPointer{T}()

function (::Type{AbstractDeviceArray})(ptr::PtrT, shape::Vararg{Integer, N}) where PtrT <: Transpiler.cli.LocalPointer{T} where {T, N}
    DeviceArray{T, N, PtrT}(shape, ptr)
end
function (::Type{AbstractDeviceArray})(ptr::PtrT, shape::NTuple{N, Integer}) where PtrT <: Transpiler.cli.LocalPointer{T} where {T, N}
    DeviceArray{T, N, PtrT}(shape, ptr)
end


text:
using Test, LinearAlgebra
using OrdinaryDiffEq, DiffEqBase

const ALGOS = (
    (SymplecticEuler, true, 1),
    (VelocityVerlet, false, 2),
    (VerletLeapfrog, true, 2),
    (PseudoVerletLeapfrog, true, 2),
    (McAte2, true, 2),
    (Ruth3, true, 3),
    (McAte3, true, 3),
    (CandyRoz4, true, 4),
    (McAte4, true, 4),
    (CalvoSanz4, true, 4),
    # (McAte42, true, 4), known to be broken
    (McAte5, true, 5),
    (Yoshida6, true, 6),
    (KahanLi6, true, 6),
    (McAte8, true, 8),
    (KahanLi8, true, 8),
    (SofSpa10, true, 10),
)

function dp(p, q, pa, t)
    0q .+ pa[2]
end

function dq(p, q, pa, t)
    p .* pa[1]
end

dp(res, p, q, pa, t) = (res .= dp(p, q, pa, t))
dq(res, p, q, pa, t) = (res .= dq(p, q, pa, t))

dynode(iip, dp, dq) = DynamicalODEFunction{iip}(dp, dq)

# [0:1] used in dp, dq; [3:4] start values for p0, q0
const PARAMS = ((1.0, 0.1, 1.0, 0.0), (0.1, 1.0, 1.0, -1.0))
const IIPS = (true, false)
const TSPAN = (0.0, 1.0)

solution(t, w) = (w[2] * t + w[3], (w[2] / 2 * t + w[3]) * w[1] * t + w[4])
apa(iip::Bool, x) = iip ? vcat.(x) : x
errorbound(dt, d, x) = 100 * abs(dt)^d + 1000 * eps(norm(x))
function printerrors(text, calc, solution, pa, t1)
    print(text, ": ")
    print(norm(calc[1] - solution(t1, pa)[1]), " ")
    print(norm(calc[2] - solution(t1, pa)[2]))
    println()
end

@testset "symplectic $alg-$iip-$pa" for (alg, x, d) in ALGOS, iip in IIPS, pa = PARAMS

    dt = 0.01
    tspan = TSPAN
    t0, t1 = tspan
    dynfun = dynode(iip, dp, dq)
    p0, q0 = apa(iip, solution(t0, pa))
    prob = DynamicalODEProblem(dynfun, p0, q0, tspan, pa)

    if x || pa[1] == 1
        sol = solve(prob, alg(); dt = dt)
        calc = sol(t1)
        # printerrors("$alg-$iip-$pa", calc, solution, pa, t1)
        @test calc[1] ≈ solution(t1, pa)[1] rtol = errorbound(dt, d, calc[1])
        @test calc[2] ≈ solution(t1, pa)[2] rtol = errorbound(dt, d, calc[2])
    else
        @test_throws ArgumentError solve(prob, alg(); dt = dt)
    end
end

ext: jl


text:
@testset "Config" begin

    @testset "ArgParse.parse_item" begin
        robot = sb.ArgParse.parse_item(sb.Robot{Float64}, "\"--position 0.2 2.0\"")
        @test sb.get_position(robot) == [0.2, 2.0]
    end;

    @testset "get_config" begin
        config = sb.get_config(["--strategy", "TestStrategy"])
        @test length(config.ws.myPlayers) == 0
        @test typeof(sb.get_strategy(config)) <: Module
        @test config.active_strategy == "TestStrategy"
    end;

end;


text:
#using Pkg; Pkg.activate(first(splitdir(@__DIR__)))

using CUDA
using DLPack
using StaticArrays
using Zygote

#include("butane.jl")

# Advanced Sampling
#include(joinpath("..", "src", "SSAGES.jl"))
#using .SSAGES

# Collective Variables
ξ = SSAGES.DihedralAngle(1, 5, 8, 11)
# Define grid if necessary
grid = SSAGES.Grid(-π, π, 64, true)
# Sampling algorithm
sampling_method = SSAGES.ABF(ξ, grid, 100)

context = hoomd.context.current
sampler = SSAGES.bind(context, sampling_method)

# Run simulation
hoomd.run(100_000)
#hoomd.run(400_000)

# Data analysis
using Plots
using Statistics

θ = let #grid = sampler.grid
    Δ = only(grid.Δl)
    r = range(Δ / 2 + only(grid.r₀); step = Δ, length = only(grid.size))
    collect(r)
end

dA = let sm = sampler.sample_and_bias.update
    F = Array(vec(sm.ΣF) ./ max.(sm.hist, 1))
    dA = (F .- mean(F)) ./ std(F)
end

plot(θ, dA; label = "")
scatter!(θ, dA; label = "")

@inline function withgrad(f::F, ps) where {F}
    ξ, pb = Zygote._pullback(Zygote.Context(), f, ps...)
    ∂ξ = Zygote.tailmemaybe(pb(Zygote.sensitivity(ξ)))
    #ξ, back = Diffractor.∂⃖¹(f, ps)
    #∂ξ = back(ps)
    return ξ, SVector(∂ξ)
end

function test(sm)
    ss = sm.snapshot
    cv = sm.ξ
    rs = ss.positions
    ps = reinterpret(reshape, SVector{3, eltype(rs)}, rs)
    ids = cv.indices
    inds = cv.ξ.indices
    f = SSAGES.CollectiveVariables.cv_function(cv)
    g = ps -> withgrad(f, ps)
    #g = ps -> ForwardDiff.gradient(f, ps)
    #ξ = CUDA.fill(0.0)
    #u = CUDA.fill(zeros(SVector{3}), 4)
    #v = reinterpret(reshape, SVector{4, SVector{3, Float64}}, u)

    #function gputest(ps, ids, inds, u, v, x)
    #    i = threadIdx().x
    #    if i ≤ length(u)
    #        u[i] = ps[ids[inds[i]] + 1]
    #    end
    #    if i == length(u) + 1
    #        sync_threads()
    #        t = g(v[])
    #        x[] = t[1]
    #    end
    #    return nothing
    #end

    #@show g(SVector((SVector((rand() for _ in 1:3)...) for _ in 1:4)...))

    #@cuda threads=(length(u) + 1) gputest(ps, ids, inds, u, v, ξ)
    #CUDA.@allowscalar g(reinterpret(reshape, SVector{4, SVector{3, Float64}}, view(ps, view(ids, inds) .+ 1))[])
    #@device_code_warntype interactive=true gputest(g, ps, ids, inds, u, v, ξ)
    #return x, u
    val_and_grad(g, rs, ids, inds)
end

function val_and_grad(f::F, rs, ids, inds) where {F}
    # size(rs) == (3, 4)
    T = eltype(rs)
    #g = (p₁, p₂, p₃, p₄) -> withgrad(f, p₁, p₂, p₃, p₄)
    ps = reinterpret(reshape, SVector{3, T}, rs)
    x = CUDA.fill(zero(T))
    u = CUDA.fill(zeros(SVector{3}), 4)
    v = reinterpret(reshape, SVector{4, SVector{3, Float64}}, u)

    function kernel(x, u, v, ps, ids, inds)
        i = threadIdx().x
        @inbounds if i ≤ length(u)
            u[i] = ps[ids[inds[i]] + 1]
        end
        @inbounds if i == length(u) + 1
            sync_threads()
            ξ, ∂ξ = f(v[])
            x[] = ξ
            v[] = ∂ξ
        end
        return nothing
    end

    @cuda threads=(length(u) + 1) kernel(x, u, v, ps, ids, inds)

    return x, u
end

#let
#    sm = sampler.sample_and_bias.update
#    ss = sm.snapshot
#    cv = sm.ξ
#    rs = ss.positions
#    cv(rs)
#end

#test(sampler.sample_and_bias.update)


text:
function solution_error(num_pts, u_extended)
    x = range(0,1,length=num_pts)
    y = range(0,1,length=num_pts)'

    u_true = zeros((num_pts, num_pts));

    for ix = 1:num_pts
        for iy = 1:num_pts
            u_true[ix,iy]=u(x[ix],y[iy])
        end
    end

    error = u_true-u_extended;
    error = error[:]
    return x,y, error
end


text:
using JuMP
using Gurobi
using Plasmo

## Model on x
# Min 16x[1] + 10x[2]
# s.t. x[1] + x[2] <= 1OutputFlag=0
#      x ∈ {0,1}
m1 = Model(solver=GurobiSolver(OutputFlag=0))
@variable(m1, xm[i in 1:2],Bin)
@constraint(m1, xm[1] + xm[2] <= 1)
@objective(m1, Max, 16xm[1] + 10xm[2])

## Model on y`
# Max  4y[2]
# s.t. y[1] + y[2] <= 1
#      8x[1] + 2x[2] + y[2] + 4y[2] <= 10
#      x, y ∈ {0,1}

m2 = Model(solver=GurobiSolver(OutputFlag=0))
@variable(m2, xs[i in 1:2],Bin)
@variable(m2, y[i in 1:2], Bin)
@constraint(m2, y[1] + y[2] <= 1)
@constraint(m2, 8xs[1] + 2xs[2] + y[1] + 4y[2] <= 10)
@objective(m2, Max, 4y[2])

## Plasmo Graph
g = ModelGraph()
setsolver(g, GurobiSolver(OutputFlag=0))
n1 = add_node(g)
setmodel(n1,m1)
n2 = add_node(g)
setmodel(n2,m2)
add_edge(g, n1,n2)
## Linking
# m1[x] = m2[x]  ∀i ∈ {1,2}
@linkconstraint(g, [i in 1:2], n1[:xm][i] == n2[:xs][i])


text:
using Test
using Metatheory
using TermInterface

function prove(t, ex, steps = 1, timeout = 10, eclasslimit = 5000)
  params = SaturationParams(
    timeout = timeout,
    eclasslimit = eclasslimit,
    # scheduler=Schedulers.ScoredScheduler, schedulerparams=(1000,5, Schedulers.exprsize))
    scheduler = Schedulers.BackoffScheduler,
    schedulerparams = (6000, 5),
  )

  hist = UInt64[]
  push!(hist, hash(ex))
  for i in 1:steps
    g = EGraph(ex)

    exprs = [true, g[g.root]]
    ids = [addexpr!(g, e)[1].id for e in exprs]

    goal = EqualityGoal(exprs, ids)
    params.goal = goal
    saturate!(g, t, params)
    ex = extract!(g, astsize)
    if !TermInterface.istree(ex)
      return ex
    end
    if hash(ex) ∈ hist
      return ex
    end
    push!(hist, hash(ex))
  end
  return ex
end

@testset "Prop logic" begin
  or_alg = @theory p q r begin
    ((p ∨ q) ∨ r) == (p ∨ (q ∨ r))
    (p ∨ q) == (q ∨ p)
    (p ∨ p) --> p
    (p ∨ true) --> true
    (p ∨ false) --> p
  end

  and_alg = @theory p q r begin
    ((p ∧ q) ∧ r) == (p ∧ (q ∧ r))
    (p ∧ q) == (q ∧ p)
    (p ∧ p) --> p
    (p ∧ true) --> p
    (p ∧ false) --> false
  end

  comb = @theory p q r begin
    # DeMorgan
    ¬(p ∨ q) == (¬p ∧ ¬q)
    ¬(p ∧ q) == (¬p ∨ ¬q)
    # distrib
    (p ∧ (q ∨ r)) == ((p ∧ q) ∨ (p ∧ r))
    (p ∨ (q ∧ r)) == ((p ∨ q) ∧ (p ∨ r))
    # absorb
    (p ∧ (p ∨ q)) --> p
    (p ∨ (p ∧ q)) --> p
    # complement
    (p ∧ (¬p ∨ q)) --> p ∧ q
    (p ∨ (¬p ∧ q)) --> p ∨ q
  end

  negt = @theory p begin
    (p ∧ ¬p) --> false
    (p ∨ ¬(p)) --> true
    ¬(¬p) == p
  end

  impl = @theory p q begin
    (p == ¬p) --> false
    (p == p) --> true
    (p == q) --> (¬p ∨ q) ∧ (¬q ∨ p)
    (p => q) --> (¬p ∨ q)
  end

  fold = @theory p q begin
    (p::Bool == q::Bool) => (p == q)
    (p::Bool ∨ q::Bool)  => (p || q)
    (p::Bool => q::Bool) => ((p || q) == q)
    (p::Bool ∧ q::Bool)  => (p && q)
    ¬(p::Bool)           => (!p)
  end

  t = or_alg ∪ and_alg ∪ comb ∪ negt ∪ impl ∪ fold

  ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
  @test prove(t, ex, 3, 10, 5000)

  @test @areequal t true ((¬p == p) == false)
  @test @areequal t true ((¬p == ¬p) == true)
  @test @areequal t true ((¬p ∨ ¬p) == ¬p) (¬p ∨ p) ¬(¬p ∧ p)
  @test @areequal t true ((p => (p ∨ p)))
  @test @areequal t true ((p => (p ∨ p)) == ((¬(p) ∧ q) => q)) == true

  # Frege's theorem
  @test @areequal t true (p => (q => r)) => ((p => q) => (p => r))

  # Demorgan's
  @test @areequal t true (¬(p ∨ q) == (¬p ∧ ¬q))

  # Consensus theorem
  @test @areequal t ((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)) ((x ∧ y) ∨ (¬x ∧ z))
end

# https://www.cs.cornell.edu/gries/Logic/Axioms.html
# The axioms of calculational propositional logic C are listed in the order in
# which they are usually presented and taught. Note that equivalence comes
# first. Note also that, after the first axiom, we take advantage of
# associativity of equivalence and write sequences of equivalences without
# parentheses. We use == for equivalence, | for disjunction, & for conjunction,
# ~ for negation (not), => for implication, and <= for consequence.
#
# Associativity of ==: ((p == q) == r) == (p == (q == r))
# Symmetry of ==: p == q == q == p
# Identity of ==: true == q == q
#
# Definition of false: false == ~true
# Distributivity of not: ~(p == q) == ~p == q
# Definition of =/=: (p =/= q) == ~(p == q)
#
# Associativity of |: (p | q) & r == p | (q | r)
# Symmetry of |: p | q == q | p
# Idempotency of |: p | p == p
# Distributivity of |: p |(q == r) == p | q == p | r
# Excluded Middle: p | ~p
#
# Golden rule: p & q == p == q == p | q
#
# Implication: p => q == p | q == q
# Consequence: p <= q == q => p

@testset "Calculational Logic" begin
  calc = @theory p q r begin
    ((p == q) == r) == (p == (q == r))
    (p == q) == (q == p)
    (q == q) --> true

    ¬(p == q) == (¬(p) == q)
    (p != q) == ¬(p == q)

    ((p ∨ q) ∨ r) == (p ∨ (q ∨ r))
    (p ∨ q) == (q ∨ p)
    (p ∨ p) --> p
    (p ∨ (q == r)) == (p ∨ q == p ∨ r)
    (p ∨ ¬(p)) --> true

    # DeMorgan
    ¬(p ∨ q) == (¬p ∧ ¬q)
    ¬(p ∧ q) == (¬p ∨ ¬q)

    (p ∧ q) == ((p == q) == p ∨ q)

    (p => q) == ((p ∨ q) == q)
    # (p => q)            ==  (¬p ∨ q)
    # (p <= q)            =>  (q => p)
  end

  fold = @theory p q begin
    (p::Bool == q::Bool) => (p == q)
    (p::Bool ∨ q::Bool)  => (p || q)
    (p::Bool => q::Bool) => ((p || q) == q)
    (p::Bool ∧ q::Bool)  => (p && q)
    ¬(p::Bool)           => (!p)
  end

  # t = or_alg ∪ and_alg ∪ neg_alg ∪ demorgan ∪ and_or_distrib ∪
  #     absorption ∪ calc

  t = calc ∪ fold

  g = EGraph(:(((¬p == p) == false)))
  saturate!(g, t)
  extract!(g, astsize)

  @test @areequal t true ((¬p == p) == false)
  @test @areequal t true ((¬p == ¬p) == true)
  @test @areequal t true ((¬p ∨ ¬p) == ¬p) (¬p ∨ p) ¬(¬p ∧ p)
  @test @areequal t true ((p => (p ∨ p)) == true)
  params = SaturationParams(timeout = 12, eclasslimit = 10000, schedulerparams = (1000, 5))

  @test areequal(t, true, :(((p => (p ∨ p)) == ((¬(p) ∧ q) => q)) == true); params = params)

  # Frege's theorem
  # params = SaturationParams(timeout=12, eclasslimit=15000, scheduler=Schedulers.ScoredScheduler)
  # params = SaturationParams(timeout=12, eclasslimit=15000, schedulerparams=(500, 2))
  @test_skip areequal(t, true, :((p => (q => r)) => ((p => q) => (p => r))); params = params)

  # Demorgan's
  @test @areequal t true (¬(p ∨ q) == (¬p ∧ ¬q))

  # Consensus theorem
  # @test_skip
  areequal(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)), :((x ∧ y) ∨ (¬x ∧ z)); params = params)
end


text:
using DEBmicroTrait
using CSV, DataFrames, Statistics
using HypothesisTests
using JLD

########################################
# I/O
df_isolates     = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates2traits.csv", DataFrame, missingstring="N/A")
########################################

########################################
V_cell          = DEBmicroTrait.genome_size_to_cell_volume(convert(Array{Float64,1}, df_isolates.Genome_size))
Min_gen_time    = df_isolates.Min_gen_time
Gram_stain      = convert(Array{String,1}, df_isolates.gram_stain)
#α               = DEBmicroTrait.constrain_enzyme_allocation(V_cell, Min_gen_time, Gram_stain, df_isolates.z_hydrolases./df_isolates.Genome_size*1e6)
zh              =  df_isolates.z_hydrolases./df_isolates.Genome_size*1e6
α_X             =  1e-2*(df_isolates.z_hydrolases./df_isolates.Genome_size*1e6)./maximum(df_isolates.z_hydrolases./df_isolates.Genome_size*1e6)
########################################

########################################
# statistics
df_isolates_out = DataFrame()
df_isolates_out.response = df_isolates.Rhizosphere_response
df_isolates_out.zh = df_isolates.z_hydrolases./df_isolates.Genome_size*1e6
df_isolates_out.α_X = α_X

df_p_α = filter(x->(x.response.=="positive") , df_isolates_out)
df_n_α = filter(x->(x.response.=="negative") , df_isolates_out)
df_u_α = filter(x->(x.response.=="undefined") , df_isolates_out)
zh_kw   = KruskalWallisTest(df_p_α.zh, df_n_α.zh)
a_kw   = KruskalWallisTest(df_p_α.α_X, df_n_α.α_X)
########################################

########################################
# I/O
save("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates_enzymes.jld", "zh", zh, "alpha", α_X)
########################################

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using LRUCache
using Test
using Random
using Base.Threads

# Test insertion and reinsertion ordering
@testset "Single-threaded insertion and reinsertion" begin
    for cache in [LRU(; maxsize = 100), LRU{Int, Int}(; maxsize = 100)]
        r = 1:100
        for i in reverse(r)
            cache[i] = i
        end
        @test collect(cache) == collect(i=>i for i in r)

        @threads for i = 1:10:100
            @test haskey(cache, i)
            @test !haskey(cache, 100+i)
        end

        # reinsert in random order
        p = randperm(100)
        for i in reverse(p)
            cache[i] = i
        end
        @test collect(cache) == collect(i=>i for i in p)
    end
end

@testset "Multi-threaded insertion and reinsertion" begin
    cache = LRU{Int, Int}(; maxsize = 100)

    r = 1:100
    @threads for i in reverse(r)
        cache[i] = i
    end

    @threads for i = 1:10:100
        @test haskey(cache, i)
        @test !haskey(cache, 100+i)
    end

    @test Set(cache) == Set(i=>i for i=1:100)

    @threads for i in 100:-1:1
        cache[i] = 2*i
    end

    @test Set(cache) == Set(i=>2*i for i=1:100)
end

@testset "Single-threaded getting and setting" begin
    cache = LRU{Int, Int}(; maxsize = 50*sizeof(Int), by = sizeof)

    for i in 100:-1:1
        @test 2*i == get(cache, i, 2*i)
        @test 2*i == get(()->2*i, cache, i)
        @test i == (iseven(i) ? get!(cache, i, i) : get!(()->i, cache, i))
        @test i == get!(cache, i, 2*i)
        @test i == get!(()->error("this should not happen"), cache, i)
        @test i == get(cache, i, 2*i)
        @test i == get(()->2*i, cache, i)
    end
    for i in 1:50
        @test haskey(cache, i)
        @test !haskey(cache, i+50)
        @test_throws KeyError getindex(cache, i+50)
    end

    @test collect(cache) == collect(i=>i for i in 1:50)

    p = randperm(50)
    for i in reverse(p)
        cache[i] = i
    end
    @test collect(cache) == collect(i=>i for i in p)

    p10 = p[1:10]
    resize!(cache; maxsize = 10*sizeof(Int))
    @test collect(cache) == collect(i=>i for i in p10)

    resize!(cache; maxsize = 50*sizeof(Int))
    @test collect(cache) == collect(i=>i for i in p10)

    for i in 50:-1:1
        @test get!(cache, i, 2*i) == (i in p10 ? i : 2*i)
    end

    for i in reverse(p10)
        @test cache[i] == i
    end
    resize!(cache; maxsize = 10*sizeof(Int))
    @test collect(cache) == collect(i=>i for i in p10)

    delete!(cache, p10[1])
    @test !haskey(cache, p10[1])

    @test p10[10] == pop!(cache, p10[10])
    @test !haskey(cache, p10[10])

    @test length(empty!(cache)) == 0
    @test isempty(cache)
    for i in p10
        @test !haskey(cache, p10)
    end
end

@testset "Multi-threaded getting and setting" begin
    cache = LRU{Int, Int}(; maxsize = 50*sizeof(Int), by = sizeof)

    @threads for i in 50:-1:1
        @test 2*i == get(cache, i, 2*i)
        @test 2*i == get(()->2*i, cache, i)
        @test i == (iseven(i) ? get!(cache, i, i) : get!(()->i, cache, i))
        @test i == get!(cache, i, 2*i)
        @test i == get!(()->error("this should not happen"), cache, i)
        @test i == get(cache, i, 2*i)
        @test i == get(()->2*i, cache, i)
    end

    p = randperm(50)
    p10 = p[1:10]
    @threads for i in p10
        @test cache[i] == i
    end
    resize!(cache; maxsize = 10*sizeof(Int))
    @test Set(cache) == Set(i=>i for i in p10)

    resize!(cache; maxsize = 50*sizeof(Int))
    @test Set(cache) == Set(i=>i for i in p10)

    delete!(cache, p10[1])
    @test !haskey(cache, p10[1])
    @test_throws KeyError getindex(cache, p10[1])

    @test p10[10] == pop!(cache, p10[10])
    @test !haskey(cache, p10[10])
    @test_throws KeyError getindex(cache, p10[1])
end

include("originaltests.jl")


text:
"""
$(SIGNATURES)

Calculate the Direction Of Arrival (DOA) with respect to the manifold
`manifold`. The calculated DOA is returned as `Spherical` object, where the
azimuth is counted counter-clockwise from the x axis and the elevation is
counted from the xy-plane towards the zenith.
If the initial azimuth `init_az` and elevation `init_el` are provided, this
will use the hill climbing algorithm.
The hill climbing algorithm should be a lot faster than plain maximum search,
but may only find a local maximum, which might be desired.

# Examples
Using the signal subspace:
```julia-repl
julia> est_doa(manifold, a -> norm(signal_subspace' * a) / norm(a))
```

Using the noise subspace (MUSIC):
```julia-repl
julia> est_doa(manifold, a -> norm(a) / norm(noise_subspace' * a))
```

# Arguments:
- `manifold::AbstractManifold`: Manifold of the antenna array
- `reduction_function`: Function that reduces the manifold vector to a scalar
"""
function est_doa(
    manifold::AbstractManifold,
    reduction_function;
    init_az = nothing,
    init_el = nothing,
    kwargs...
)
    p = Pattern(manifold, reduction_function; kwargs...)
    if isnothing(init_az) || isnothing(init_el)
        _, idx = findmax(p.values)
    else
        init_az_idx = floor(Int, mod2pi(init_az) / Float64(p.azs.step) + 1)
        init_el_idx = floor(Int, (π / 2 - init_el) / Float64(p.els.step) + 1)
        _, idx = hill_climbing(p.values, init_el_idx, init_az_idx)
    end
    az = p.azs[idx.I[2]]
    el = p.els[idx.I[1]]
    return Spherical(1, az, π / 2 - el)
end

"""
$(SIGNATURES)

This is a convenience function of `est_doa`.
See `est_doa` for more details.
"""
function est_doa_by_noise_subspace(
    manifold::AbstractManifold{N},
    noise_subspace;
    init_az = nothing,
    init_el = nothing,
    kwargs...
) where N
    size(noise_subspace, 1) == N || ArgumentError("Number of antenna channels must match")
    est_doa(
        manifold,
        a -> norm(a) / norm(noise_subspace' * a);
        init_az = init_az,
        init_el = init_el,
        kwargs...
    )
end

"""
$(SIGNATURES)

This is a convenience function of `est_doa`.
See `est_doa` for more details.
"""
function est_doa_by_signal_subspace(
    manifold::AbstractManifold{N},
    signal_subspace;
    init_az = nothing,
    init_el = nothing,
    kwargs...
) where N
    size(signal_subspace, 1) == N || ArgumentError("Number of antenna channels must match")
    est_doa(
        manifold,
        a -> norm(signal_subspace' * a) / norm(a);
        init_az = init_az,
        init_el = init_el,
        kwargs...
    )
end

"""
$(SIGNATURES)

This is a convenience function of `est_doa`.
See `est_doa` for more details.
"""
est_doa_by_music(m, ns; init_az = nothing, init_el = nothing, k...) =
    est_doa_by_noise_subspace(m, ns, init_az = init_az, init_el = init_el, k...)

"""
$(SIGNATURES)

Multi-dimensionial hill climbing.

# Arguments:
- `values::AbstractArray`: Some multi-dimensional array
- `inits`: Starting indices
"""
function hill_climbing(values, inits...)
    ndims(values) == length(inits) || ArgumentError("Number of inits does not fit dimensions of values")
    indices = CartesianIndices(values)
    lower_index_bounds, upper_index_bounds = first(indices), last(indices)
    step_size = oneunit(lower_index_bounds)
    current_index = CartesianIndex(inits...)
    current_value = values[current_index]
    while true
        next_index = current_index
        for J in max(lower_index_bounds, current_index - step_size):min(upper_index_bounds, current_index + step_size)
            if values[J] > current_value
                next_index = J
                current_value = values[J]
            end
        end
        if current_index == next_index
            return current_value, current_index
        else
            current_index = next_index
        end
    end
end


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb04md(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    LDA = NMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    SCALE = Array{Float64,1}(undef, NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    MAXRED = parse(Float64, replace(vs[2],'D'=>'E'))
if ( N<=0 || N>NMAX )
else

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    MAXRED, INFO = SLICOT.mb04md!(N, MAXRED, A, SCALE)
    println(io, "MAXRED = $MAXRED")
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io,"A:")
    _nc = N
    _nr = N
    show(io,"text/plain",A[1:_nr,1:_nc])
    println(io,)

# interp output 2
    println(io,"SCALE:")
    _nr = N
    show(io,"text/plain",SCALE[1:_nr])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


text:
mutable struct RunDump <: Output
  path_prefix::String
  num_trajectories::Int
  RunDump(path_prefix::AbstractString, ::Integer) = new(path_prefix, 0)
end

function pushtrajectory!(d::RunDump, state::Simulation.SimState, ::Simulation.SimParams, callback::DetectionCallback)
  d.num_trajectories += 1
  f = jldopen(d.path_prefix*"_$(d.num_trajectories).jld2", "w", compress=true)
  try
    Simulation.saveparams(f, state)
    saveparams(f, callback)
  finally
    close(f)
  end
  nothing
end


text:
mutable struct SimulationRef
    raw::String
    models::String
    results::String
    run_count::Dict{Int64, Dict{Int64, Int64}}
    date_ref::Dict{Int64, Dates.DateTime}
    current_time::Dates.DateTime
    reset::Bool
end

function _initialize_sim_ref(steps::Int64, stages_keys::Base.KeySet)

    count_dict = Dict{Int64, Dict{Int64, Int64}}()

    for s in 1:steps
        count_dict[s] = Dict{Int64, Int64}()
        for st in stages_keys
            count_dict[s][st] = 0
        end
    end

    return sim_ref = SimulationRef("init",
                                   "init",
                                   "init",
                                   count_dict,
                                   Dict{Int64, Dates.DateTime}(),
                                   Dates.now(),
                                   true
                                   )

end

mutable struct Simulation
    steps::Int64
    stages::Vector{_Stage}
    valid_timeseries::Bool
    daterange::NTuple{2, Dates.DateTime} #Inital Time of the first forecast and Inital Time of the last forecast
    ref::SimulationRef

    function Simulation(base_name::String,
                        steps::Int64,
                        stages::Dict{Int64, Stage},
                        simulation_folder::String;
                        verbose::Bool = false, kwargs...)

    sim_ref = _initialize_sim_ref(steps, keys(stages))

    dates, validation, stages_vector = _build_simulation!(sim_ref,
                                                        base_name,
                                                        steps,
                                                        stages,
                                                        simulation_folder;
                                                        verbose = verbose, kwargs...)

    new(steps,
        stages_vector,
        validation,
        dates,
        sim_ref)

    end

end

################# accessor functions ####################

get_steps(s::Simulation) = s.steps
get_daterange(s::Simulation) = s.daterange


text:
function depthfirstsearch!(A::Matrix{Float64}, u::Int64,
			visited::Vector{Bool}, contents::Vector{Int64})
	n = size(A,1)
	for v=1:n
		adj = A[u, v]
		if (abs(adj) > 0.5) && (u != v)
			if visited[v] == false
				visited[v] = true
				push!(contents, v)
				contents = depthfirstsearch!(A, v, visited, contents)
			end
		end
	end
	return contents
end

function connectedcomponents(A::Matrix{Float64})
	n = size(A,1)
	components = 0
	visited = Array{Bool, 1}(undef, n)
	visited[:] .= false
	table = Array{Array{Int64, 1}, 1}([])
	for u=1:n
		if visited[u] == false
			visited[u] = true
			components += 1
			contents = Array{Int64, 1}([])
			push!(contents, u)
			contents = depthfirstsearch!(A, u, visited, contents)
			push!(table, contents)
		end
	end
	return components, table
end

function partitionadjacency(A::Matrix{Float64}, v::Vector{Int64})
	n = size(A,1)
	n1 = size(v,1)
	Atemp = zeros(Float64, n, n1)
	Afinal = zeros(Float64, n1, n1)
	for (i, val) ∈ enumerate(v)
		Atemp[:,i] = A[:,val]
	end
	for (i, val) ∈ enumerate(v)
		Afinal[i,:] = Atemp[val,:]
	end
	return Afinal
end

function avcompontentsize(A::Matrix{Float64})
	components, table = connectedcomponents(A)
	s = 0
	for i=1:components
		s += size(table[i],1)
	end
	return s/components
end

function clustersize(A::Matrix{Float64}, v::Vector{Int64})
	Anew = partitionadjacency(A, v)
	avs = avcompontentsize(Anew)
	return avs
end

function nodevoltages(L::Matrix{Float64}, P::Vector{Float64})
	n = size(L,1)
	Lred = L[2:end,2:end]
	C = zeros(Float64, n, n)
	C[2:end,2:end] = inv(Lred)
	return C*P
end

function edgecurrents(theta::Vector{Float64}, E::Matrix{Float64}, X::Vector{Float64})
	tmp = transpose(E)*theta
	return tmp./X
end

function fracture!(net::Graph, nodeset::Vector{Int64}, P::Vector{Float64}, X::Vector{Float64}, Finit::Vector{Float64}, info::NetworkRecorder, t::Int64, alpha::Float64, fmaxtemp::Float64, seq::Bool)

	tol = 1e-5
	simtime = t + 1

	#Check the status of P, see if has any sources or sinks. If not, halt recursion
	nodecardinality = length(nodeset)
	P1 = zeros(Float64, nodecardinality)
	sourcecounter = 0
	weighedsourcecounter = 0.0
	sinkcounter = 0
	for i=1:nodecardinality
		vi = nodeset[i]
		P1[i] = P[vi]
		if P1[i] > tol
			sourcecounter += 1
			weighedsourcecounter += P1[i]
		end
		if P1[i] < -tol
			sinkcounter += 1
		end
	end
	if (sourcecounter == 0) || (sinkcounter == 0)
		if simtime ∉ info.simulationtime
			push!(info.simulationtime, simtime)
			push!(info.liveedges, 0)
			push!(info.livesources, 0)
			push!(info.weightedlivesources, 0.0)
			push!(info.livesinks, 0)
			push!(info.totalpower, 0.0)
			if seq
				flowvec = zeros(Float64, info.liveedges[1])
				flowvec[:] .= NaN
				push!(info.flowsequence, flowvec)
			end
		end
		return info
	end

	# #Balance P1 - if there is power excess, increase loads
	# delta = sum(P1)
	# individualdelta = delta/Float64(sinkcounter)
	# for i=1:nodecardinality
	# 	if P1[i] < -tol
	# 		P1[i] -= individualdelta
	# 	end
	# end

	#Balance P1 homogeneously
	delta = sum(P1)/2.0
	individualdelta_sink = delta/Float64(sinkcounter)
	individualdelta_source = delta/Float64(sourcecounter)
	for i=1:nodecardinality
		if P1[i] < -tol
			P1[i] -= individualdelta_sink
		end
		if P1[i] > tol
			P1[i] -= individualdelta_source
		end
	end

	# # Balance P1 heterogeneously
	# delta = sum(P1)/2.0 #This is the amount of surplus given to sources (or sinks) to deal with
	# avP = 0.0
	# for i=1:nodecardinality
	# 	if P1[i] > tol
	# 		avP += P1[i]
	# 	end
	# end
	# avP /= sourcecounter
	# for i=1:nodecardinality
	# 	if P1[i] < -tol
	# 		P1[i] -=  abs(P1[i])*delta/(Float64(sinkcounter)*avP)
	# 	end
	# 	if P1[i] > tol
	# 		P1[i] -= abs(P1[i])*delta/(Float64(sourcecounter)*avP)
	# 	end
	# end
	# println("sum", "    ", sum(P1))

	#Find which edges are connected to the node-set
	edgeset = Array{Int64, 1}([])
	for vi ∈ nodeset
		v = net.nodelist[vi]
		for e ∈ v.outedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
		for e ∈ v.inedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
	end
	edgecardinality = length(edgeset)

	#Construct X
	X1 = zeros(Float64, edgecardinality)
	for i=1:edgecardinality
		X1[i] = X[edgeset[i]]
	end

	#Construct E1
	E1 = zeros(Float64, nodecardinality, edgecardinality)
	for i=1:nodecardinality
		vi = nodeset[i]
		v = net.nodelist[vi]
		for j=1:edgecardinality
			if edgeset[j] ∈ v.outedges
				E1[i,j] = 1.0
			end
			if edgeset[j] ∈ v.inedges
				E1[i,j] = -1.0
			end
		end
	end

	# A = adjacencymatrix(E1)
	# comt, tab = connectedcomponents(A)
	# println("com... = ", comt)
	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E1, 1.0./X1)
	theta = nodevoltages(L, P1)
	F = edgecurrents(theta, E1, X1)./fmaxtemp

	#Record the number of edges in a tuple of recursion depth and number of edges
	if simtime ∈ info.simulationtime
		i = findall(y->y==simtime, info.simulationtime)
		info.liveedges[i] .+= edgecardinality
		info.livesources[i] .+= sourcecounter
		info.weightedlivesources .+= weighedsourcecounter
		info.livesinks[i] .+= sinkcounter
		info.totalpower[i] .+= sum(abs.(P1))/2.0
		if seq
			for j=1:edgecardinality
				ej = edgeset[j]
				info.flowsequence[i][ej] = abs(F[j])
			end
		end
	else
		push!(info.simulationtime, simtime)
		push!(info.liveedges, edgecardinality)
		push!(info.livesources, sourcecounter)
		push!(info.weightedlivesources, weighedsourcecounter)
		push!(info.livesinks, sinkcounter)
		push!(info.totalpower, sum(abs.(P1))/2.0)
		if seq
			flowvec = zeros(Float64, info.liveedges[1])
			flowvec[:] .= NaN
			for j=1:edgecardinality
				ej = edgeset[j]
				flowvec[ej] = abs(F[j])
			end
			push!(info.flowsequence, flowvec)
		end
	end

	#Scan through new flow vector and record any overloads
	#If no new overloads, halt recursion
	overloadcounter = 0
	survivors = Array{Int64, 1}([])
	for i=1:edgecardinality
		ei = edgeset[i]
		if abs(F[i]) > alpha #*abs(Finit[ei])
			net.edgelist[ei].active = false
			overloadcounter += 1
		else
			push!(survivors, i)
		end
	end
	if overloadcounter == 0
		return info
	end

	#Delete overloads from E1
	newedgecard = length(survivors)
	E2 = zeros(Float64, nodecardinality, newedgecard)
	for j=1:newedgecard
		E2[:,j] = E1[:,survivors[j]]
	end

	#Check to see if its broken
	Adj = adjacencymatrix(E2)
	components, table = connectedcomponents(Adj)

	#Clear any superfluous data
	E1 = nothing
	E2 = nothing
	Adj = nothing
	L = nothing
	F = nothing

	#Call fracture function on each of the components, passing down info
	for i=1:components
		chunk = table[i]
		nodesubset = zeros(Int64, length(chunk))
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
		end
		info = fracture!(net, nodesubset, P, X, Finit, info, simtime, alpha, fmaxtemp, seq)
	end

	return info
end

function fracture2!(net::Graph, nodeset::Vector{Int64}, P::Vector{Float64}, X::Vector{Float64}, Finit::Vector{Float64}, t::Int64, alpha::Float64, fmaxtemp::Float64)

	tol = 1e-5
	simtime = t + 1

	#Check the status of P, see if has any sources or sinks. If not, halt recursion
	nodecardinality = length(nodeset)
	P1 = zeros(Float64, nodecardinality)
	sourcecounter = 0
	weighedsourcecounter = 0.0
	sinkcounter = 0
	for i=1:nodecardinality
		vi = nodeset[i]
		P1[i] = P[vi]
		if P1[i] > tol
			sourcecounter += 1
			weighedsourcecounter += P1[i]
		end
		if P1[i] < -tol
			sinkcounter += 1
		end
	end
	if (sourcecounter == 0) || (sinkcounter == 0)
		return 0
	end

	# #Balance P1 - if there is power excess, increase loads
	# delta = sum(P1)
	# individualdelta = delta/Float64(sinkcounter)
	# for i=1:nodecardinality
	# 	if P1[i] < -tol
	# 		P1[i] -= individualdelta
	# 	end
	# end

	# #Balance P1 homogeneously
	# delta = sum(P1)/2.0
	# individualdelta_sink = delta/Float64(sinkcounter)
	# individualdelta_source = delta/Float64(sourcecounter)
	# for i=1:nodecardinality
	# 	if P1[i] < -tol
	# 		P1[i] -= individualdelta_sink
	# 	end
	# 	if P1[i] > tol
	# 		P1[i] -= individualdelta_source
	# 	end
	# end

	# Balance P1 heterogeneously
	delta = sum(P1)/2.0 #This is the amount of surplus given to sources (or sinks) to deal with
	avP = 0.0
	for i=1:nodecardinality
		if P1[i] > tol
			avP += P1[i]
		end
	end
	avP /= sourcecounter
	for i=1:nodecardinality
		if P1[i] < -tol
			P1[i] -=  abs(P1[i])*delta/(Float64(sinkcounter)*avP)
		end
		if P1[i] > tol
			P1[i] -= abs(P1[i])*delta/(Float64(sourcecounter)*avP)
		end
	end
	# println("sum", "    ", sum(P1))

	#Find which edges are connected to the node-set
	edgeset = Array{Int64, 1}([])
	for vi ∈ nodeset
		v = net.nodelist[vi]
		for e ∈ v.outedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
		for e ∈ v.inedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
	end
	edgecardinality = length(edgeset)

	#Construct X
	X1 = zeros(Float64, edgecardinality)
	for i=1:edgecardinality
		X1[i] = X[edgeset[i]]
	end

	#Construct E1
	E1 = zeros(Float64, nodecardinality, edgecardinality)
	for i=1:nodecardinality
		vi = nodeset[i]
		v = net.nodelist[vi]
		for j=1:edgecardinality
			if edgeset[j] ∈ v.outedges
				E1[i,j] = 1.0
			end
			if edgeset[j] ∈ v.inedges
				E1[i,j] = -1.0
			end
		end
	end

	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E1, 1.0./X1)
	theta = nodevoltages(L, P1)
	F = edgecurrents(theta, E1, X1)./fmaxtemp

	#Scan through new flow vector and record any overloads
	overloadcounter = 0
	survivors = Array{Int64, 1}([])
	for i=1:edgecardinality
		ei = edgeset[i]
		if abs(F[i]) > alpha #*abs(Finit[ei])
			net.edgelist[ei].active = false
			overloadcounter += 1
		else
			push!(survivors, i)
		end
	end

	#If no new overloads, halt recursion
	if overloadcounter == 0
		return edgecardinality
	end

	#Delete overloads from E1
	newedgecard = length(survivors)
	E2 = zeros(Float64, nodecardinality, newedgecard)
	for j=1:newedgecard
		E2[:,j] = E1[:,survivors[j]]
	end

	#Check to see if its broken
	Adj = adjacencymatrix(E2)
	components, table = connectedcomponents(Adj)

	#Clear any superfluous data
	E1 = nothing
	E2 = nothing
	Adj = nothing
	L = nothing
	F = nothing

	#Call fracture function on each of the components, passing down info
	descendent_survivingedges = 0
	for i=1:components
		chunk = table[i]
		nodesubset = zeros(Int64, length(chunk))
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
		end
		survivingedges = fracture2!(net, nodesubset, P, X, Finit, simtime, alpha, fmaxtemp)
		descendent_survivingedges += survivingedges
	end

	return descendent_survivingedges
end

function fracturerand!(net::Graph, nodeset::Vector{Int64}, P::Vector{Float64}, X::Vector{Float64}, Finit::Vector{Float64}, capnoise::Vector{Float64}, t::Int64, alpha::Float64, fmaxtemp::Float64)

	tol = 1e-5
	simtime = t + 1

	#Check the status of P, see if has any sources or sinks. If not, halt recursion
	nodecardinality = length(nodeset)
	P1 = zeros(Float64, nodecardinality)
	sourcecounter = 0
	weighedsourcecounter = 0.0
	sinkcounter = 0
	for i=1:nodecardinality
		vi = nodeset[i]
		P1[i] = P[vi]
		if P1[i] > tol
			sourcecounter += 1
			weighedsourcecounter += P1[i]
		end
		if P1[i] < -tol
			sinkcounter += 1
		end
	end
	if (sourcecounter == 0) || (sinkcounter == 0)
		return 0
	end

	#Balance P1 homogeneously
	delta = sum(P1)/2.0
	individualdelta_sink = delta/Float64(sinkcounter)
	individualdelta_source = delta/Float64(sourcecounter)
	for i=1:nodecardinality
		if P1[i] < -tol
			P1[i] -= individualdelta_sink
		end
		if P1[i] > tol
			P1[i] -= individualdelta_source
		end
	end

	#Find which edges are connected to the node-set
	edgeset = Array{Int64, 1}([])
	for vi ∈ nodeset
		v = net.nodelist[vi]
		for e ∈ v.outedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
		for e ∈ v.inedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
	end
	edgecardinality = length(edgeset)

	#Construct X
	X1 = zeros(Float64, edgecardinality)
	for i=1:edgecardinality
		X1[i] = X[edgeset[i]]
	end

	#Construct E1
	E1 = zeros(Float64, nodecardinality, edgecardinality)
	for i=1:nodecardinality
		vi = nodeset[i]
		v = net.nodelist[vi]
		for j=1:edgecardinality
			if edgeset[j] ∈ v.outedges
				E1[i,j] = 1.0
			end
			if edgeset[j] ∈ v.inedges
				E1[i,j] = -1.0
			end
		end
	end

	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E1, 1.0./X1)
	theta = nodevoltages(L, P1)
	F = edgecurrents(theta, E1, X1)./fmaxtemp

	#Scan through new flow vector and record any overloads
	overloadcounter = 0
	survivors = Array{Int64, 1}([])
	for i=1:edgecardinality
		ei = edgeset[i]
		if abs(F[i]) > alpha*(1.0 + capnoise[ei])
			net.edgelist[ei].active = false
			overloadcounter += 1
		else
			push!(survivors, i)
		end
	end

	#If no new overloads, halt recursion
	if overloadcounter == 0
		return edgecardinality
	end

	#Delete overloads from E1
	newedgecard = length(survivors)
	E2 = zeros(Float64, nodecardinality, newedgecard)
	for j=1:newedgecard
		E2[:,j] = E1[:,survivors[j]]
	end

	#Check to see if its broken
	Adj = adjacencymatrix(E2)
	components, table = connectedcomponents(Adj)

	#Clear any superfluous data
	E1 = nothing
	E2 = nothing
	Adj = nothing
	L = nothing
	F = nothing

	#Call fracture function on each of the components, passing down info
	descendent_survivingedges = 0
	for i=1:components
		chunk = table[i]
		nodesubset = zeros(Int64, length(chunk))
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
		end
		survivingedges = fracture2!(net, nodesubset, P, X, Finit, simtime, alpha, fmaxtemp)
		descendent_survivingedges += survivingedges
	end

	return descendent_survivingedges
end

function fracture3!(net::Graph, nodeset::Vector{Int64}, P::Vector{Float64}, X::Vector{Float64}, Finit::Vector{Float64}, t::Int64, alpha::Float64, fmaxtemp::Float64)

	tol = 1e-5
	simtime = t + 1

	#Check the status of P, see if has any sources or sinks. If not, halt recursion
	nodecardinality = length(nodeset)
	P1 = zeros(Float64, nodecardinality)
	sourcecounter = 0
	weighedsourcecounter = 0.0
	sinkcounter = 0
	for i=1:nodecardinality
		vi = nodeset[i]
		P1[i] = P[vi]
		if P1[i] > tol
			sourcecounter += 1
			weighedsourcecounter += P1[i]
		end
		if P1[i] < -tol
			sinkcounter += 1
		end
	end
	if (sourcecounter == 0) || (sinkcounter == 0)
		return 1
	end

	#Balance P1 homogeneously
	delta = sum(P1)/2.0
	individualdelta_sink = delta/Float64(sinkcounter)
	individualdelta_source = delta/Float64(sourcecounter)
	for i=1:nodecardinality
		if P1[i] < -tol
			P1[i] -= individualdelta_sink
		end
		if P1[i] > tol
			P1[i] -= individualdelta_source
		end
	end

	#Find which edges are connected to the node-set
	edgeset = Array{Int64, 1}([])
	for vi ∈ nodeset
		v = net.nodelist[vi]
		for e ∈ v.outedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
		for e ∈ v.inedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
	end
	edgecardinality = length(edgeset)

	#Construct X
	X1 = zeros(Float64, edgecardinality)
	for i=1:edgecardinality
		X1[i] = X[edgeset[i]]
	end

	#Construct E1
	E1 = zeros(Float64, nodecardinality, edgecardinality)
	for i=1:nodecardinality
		vi = nodeset[i]
		v = net.nodelist[vi]
		for j=1:edgecardinality
			if edgeset[j] ∈ v.outedges
				E1[i,j] = 1.0
			end
			if edgeset[j] ∈ v.inedges
				E1[i,j] = -1.0
			end
		end
	end

	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E1, 1.0./X1)
	theta = nodevoltages(L, P1)
	F = edgecurrents(theta, E1, X1)./fmaxtemp

	#Scan through new flow vector and record any overloads
	overloadcounter = 0
	survivors = Array{Int64, 1}([])
	for i=1:edgecardinality
		ei = edgeset[i]
		if abs(F[i]) > alpha #*abs(Finit[ei])
			net.edgelist[ei].active = false
			overloadcounter += 1
		else
			push!(survivors, i)
		end
	end

	#If no new overloads, halt recursion
	if overloadcounter == 0
		return 1
	end

	#Delete overloads from E1
	newedgecard = length(survivors)
	E2 = zeros(Float64, nodecardinality, newedgecard)
	for j=1:newedgecard
		E2[:,j] = E1[:,survivors[j]]
	end

	#Check to see if its broken
	Adj = adjacencymatrix(E2)
	components, table = connectedcomponents(Adj)

	#Clear any superfluous data
	E1 = nothing
	E2 = nothing
	Adj = nothing
	L = nothing
	F = nothing

	#Call fracture function on each of the components, passing down info
	descendent_leaves = 0
	for i=1:components
		chunk = table[i]
		nodesubset = zeros(Int64, length(chunk))
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
		end
		leaves = fracture2!(net, nodesubset, P, X, Finit, simtime, alpha, fmaxtemp)
		descendent_leaves += leaves
	end

	return descendent_leaves
end

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
using StructsOfArrays
using Base.Test

regular = complex(randn(10000), randn(10000))
soa = convert(StructOfArrays, regular)
@test regular == soa
@test_approx_eq sum(regular) sum(soa)

soa64 = convert(StructOfArrays{Complex64}, regular)
@test convert(Array{Complex64}, regular) == soa64

sim = similar(soa)
@test typeof(sim) == typeof(soa)
@test size(sim) == size(soa)

regular = complex(randn(10, 5), randn(10, 5))
soa = convert(StructOfArrays, regular)
for i = 1:10, j = 1:5
    @test regular[i, j] == soa[i, j]
end
@test size(soa, 1) == 10
@test size(soa, 2) == 5

immutable OneField
    x::Int
end

small = StructOfArrays(Complex64, 2)
@test typeof(similar(small, Complex)) === Vector{Complex}
@test typeof(similar(small, Int)) === Vector{Int}
@test typeof(similar(small, SubString)) === Vector{SubString}
@test typeof(similar(small, OneField)) === Vector{OneField}
@test typeof(similar(small, Complex128)) <: StructOfArrays


text:
const Match{A,B} = Tuple{A,B}

function fuzziness(matches)
    corrs = [m[2] for m in matches]
    maxcorr = maximum(corrs)
    corrs ./ maxcorr
end

"""Returns cross-correlation local maxima between `signal` and `chunk`
with gaussian attenuation (STD = `σ``) about `shift`."""
function alignments(signal, chunk, σ=0.0, shift=0)
    l = length(signal) # assuming length(signal) = length(chunk)
    xc = xcorr(signal, chunk)
    rng = -shift:(length(signal)-length(chunk)-shift)
    rng = -(l-shift-1):(l-shift-1)
    if σ != 0.0
        xc .*= exp.(-rng.^2/2σ^2) # apply gaussian window function
    end
    filter(x->x[2] > mincorr, local_maxima(xc))
end

function mincheck(matches::AbstractArray{T} where T<:Match, min_threshold)
    filter(m->m[2]>min_threshold, matches)
end

mincheck(matches, min_threshold) = mincheck.(matches, min_threshold)

function hitcheck(matches::AbstractArray{T} where T<:Match, hit_threshold)
    maximum(m[2] for m in region_matches) > hit_threshold
end

function hitcheck(matches, hit_threshold)
    h = hitcheck.(matches, hit_threshold)
    if any(isempty, h)
        Array{Int}[]
    else
        h
    end
end

function fuzzcheck(matches::AbstractArray{T} where T<:Match, max_fuzz)
    fuzz = fuzziness(matches)
    [matches[i] for i in eachindex(matches) if fuzz[i] < maxfuzz]
end
function alignments(s::Spectrum, ref::Spectrum, hitcorr, mincorr, σ=0.0)
    signal = s[:]
    matches = [alignments(s, c, σ, mincorr) for c in intrng_data(ref)]
    if any(isempty(matches))
        Array{Int}[]
    else
        matches
    end
end

function alignments(s::Spectrum, ref::Spectrum, hitcorr, mincorr, max_fuzz, σ=0.0)
    matches = alignments(s, ref, hitcorr, mincorr, σ)
end

function alignments(signal, chunk, σ, hitcorr, mincorr, target_n)
    a = alignments(signal, chunk, shift, σ)

end

function synthesize(ref::Spectrum, alignment)
end


text:
using Documenter: Documenter, makedocs, deploydocs
using PkgTemplates: PkgTemplates

makedocs(;
    modules=[PkgTemplates],
    authors="Chris de Graaf, Invenia Technical Computing Corporation",
    repo="https://github.com/invenia/PkgTemplates.jl/blob/{commit}{path}#L{line}",
    sitename="PkgTemplates.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://invenia.github.io/PkgTemplates.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "User Guide" => "user.md",
        "Developer Guide" => "developer.md",
        "Migrating To PkgTemplates 0.7+" => "migrating.md",
    ],
)

deploydocs(;
    repo="github.com/invenia/PkgTemplates.jl",
)


text:
"""
    DynamicArray

A dynamic aray implementation.
"""
mutable struct DynamicArray{T}
    size::Int
    capacity::Int
    data::AbstractVector{T}
    DynamicArray{T}() where T = new(0, 0, Vector{T}(undef, 0))
end

"""
    push!(a::DynamicArray{T}, x::T)

Add new element x to end of existing array a. Amortised time complexity O(1).
"""
function push!(a::DynamicArray{T}, x::T) where T
    if a.size >= a.capacity
        a.capacity = a.capacity == 0 ? 1 : a.capacity
        resize!(a, a.capacity * 2)
    end

    a.size += 1
    a[a.size] = x

    return nothing
end

"""
    pop!(a::DynamicArray)

Remove last element from array a. Amortised time complexity O(1).
"""
function pop!(a::DynamicArray)
    x = a[a.size]
    a.size -= 1

    if a.size <= a.capacity ÷ 4
        resize!(a, a.capacity ÷ 2)
    end

    return x
end

"""
    getindex(aDynamicArray, i::Int)

Get element from array a at index i.
"""
function getindex(a::DynamicArray, i::Int)
    i <= a.size || throw(BoundsError(a, i))

    return a.data[i]
end

"""
    setindex!(a::DynamicArray{T}, x::T, i::Int)

Set array a at index i to x.
"""
setindex!(a::DynamicArray{T}, x::T, i::Int) where T = (a.data[i] = x)

"""
    size(a::DynamicArray)

Get size of dynamic array.
"""
size(a::DynamicArray) = a.size

"""
    resize!(a::DynamicArray, n::Int)

Resize the array.
"""
resize!(a::DynamicArray, n::Int) = (a.capacity = n; resize!(a.data, n))

"""
    isempty(a::DynamicArray)

Check if array is empty.
"""
isempty(a::DynamicArray) = size(a) == 0