text:
@with_kw struct Cut
    arcs::Vector{Pair{Int64,Int64}}
    imbalance::Float64
    expansion::Float64
    size::Int64 = length(arcs)
end

"""
    forward_grow!(set, g, flow_matrix, capacity_matrix, reverse)
Start a breadth first search in a graph g from nodes in set following only non
satured arcs. This function update the set and add to it every nodes that can be
reached from the nodes in set. In the context of undirected graphs this is the same
as backward growing.

# Arguments
- `set::BitVector` the set of nodes to consider.
- `g::SimpleGraph` the graph to explore.
- `capacity_matrix::SparseMatrixCSC{Int64, Int64}` capacities of arcs.
- `flow_matrix::SparseMatrixCSC{Int64,Int64}` current flow of arcs.

# Return
- `nothing` (only set is modified)
"""
function forward_grow!(
    set::BitVector,
    g::SimpleGraph,
    flow_matrix::SparseMatrixCSC{Int64,Int64},
    capacity_matrix::SparseMatrixCSC{Int64,Int64},
)::Nothing
    # init the breadth first search algorithm
    q = Queue{Int64}()
    enqueue!.([q], findall(x -> x > 0, set))
    set .= false
    @inbounds while !isempty(q)
        cur = dequeue!(q)
        set[cur] = true
        @inbounds for nei in neighbors(g, cur)
            if (!set[nei]) &&
               abs(flow_matrix[cur, nei]) < capacity_matrix[cur, nei]
                set[nei] = true
                enqueue!(q, nei)
            end
        end
    end
    nothing
end

"""
	augment_flow!(flow_matrix, capacity_matrix, g, source, target)
Find an augmenting path and augment the flow of 1 unit along it. The flow_matrix
is modified for every arc (u,v) in the path it does flow_matrix[u, v] += 1 and
flow_matrix[v,u] += -1.

# Arguments
-`flow_matrix::SparseMatrixCSC{Int64,Int64}` the value of the current flow which may be modified.
-`capacity_matrix::SparseMatrixCSC{Int64, Int64}` the value of the arcs capacities.
-`g::SimpleGraph` the graph to consider.
-`source::Int64` index of the source node.
-`target::Int64` index of the target node.

# Return
- `augment::Int64` the augmentation of the flow so 1 if path is found, 0 otherwise.
"""
function augment_flow!(
    flow_matrix::SparseMatrixCSC{Int64,Int64},
    capacity_matrix::SparseMatrixCSC{Int64,Int64},
    g::SimpleGraph,
    source::Int64,
    target::Int64,
)::Int64
    # first find an augmented path from source to target with a DFS
    stack = Stack{Int64}()
    visited = falses(nv(g))
    push!(stack, source)
    visited[source] = true
    prevs = Array{Int64,1}(undef, nv(g))
    cur = source
    while (cur != target) && (!isempty(stack))
        cur = pop!(stack)
        @inbounds for nei in neighbors(g, cur)
            if (!visited[nei]) &&
               abs(flow_matrix[cur, nei]) < capacity_matrix[cur, nei]
                visited[nei] = true
                push!(stack, nei)
                prevs[nei] = cur
            end
        end
    end

    return if cur == target
        # augmenting path exists so augment flow along it
        @inbounds while cur != source
            next = cur
            cur = prevs[cur]
            flow_matrix[cur, next] += 1
            flow_matrix[next, cur] += -1
        end
        1
    else
        0
    end
end

"""
	piercing_node(cut, to_increase, to_avoid, increase_node, avoid_node, dist)
Compute which node will become a new source or target to balance the current cut.

# Arguments
-`cut::Vector{Pair{Int64, Int64}}` the current cut set of edges.
-`to_increase::BitVector` the side of the cut to increase.
-`to_avoid::BitVector` the side to the cut to not increase.
-`increase_node::Int64` the original source or target of to_increase.
-`avoid_node::Int64` the original source or target of to_avoid.
-`dist::Array{Int64, 2}` a matrix of distances (shortest path 1 per edge) between all nodes of the graph.

# Return
-`node::Int64` the new piercing node to balance the cut.
"""
function piercing_node(
    cut::Vector{Pair{Int64,Int64}},
    to_increase::BitVector,
    to_avoid::BitVector,
    source_increase::Bool,
    dist::Array{Int64,2},
)::Int64
    # first heuristic
    nodes = map(arc -> to_increase[arc.first] ? arc.second : arc.first, cut)
    @inbounds res = findfirst(p -> ~(to_increase[p] | to_avoid[p]), nodes)
    if isnothing(res)
        # second heuristic
        increase_node, avoid_node = source_increase ? (1, 2) : (2, 1)
        res = findmax(map(p -> dist[p, avoid_node] - dist[p, increase_node], nodes))[2]
    end
    return nodes[res]
end

"""
	flowcutter(graph, source, target)
Computes multiple cuts more and more balanced in a graph g.

# Arguments
-`graph::SimpleGraph{Int64}` the graph to consider.
-`source::Int64` index of the source node.
-`target::Int64` index of the target node.

# Return
-`cuts::Vector{Cut}` all the cuts computed by flowcutter.
"""
function flowcutter(
    graph::SimpleGraph{Int64},
    source::Int64,
    target::Int64,
)::Vector{Cut}
    g = copy(graph)
    add_vertex!(g)
    add_vertex!(g)
    n = nv(g)

    super_s = n - 1
    super_t = n

    dist = [dijkstra_shortest_paths(g, source).dists dijkstra_shortest_paths(g, target).dists]

    add_edge!(g, super_s, source)
    add_edge!(g, target, super_t)

    flow_matrix = spzeros(Int64, n, n)
    capacity_matrix = SparseMatrixCSC{Int64,Int64}(adjacency_matrix(g))
    @inbounds capacity_matrix[super_s, source] = typemax(Int64)
    @inbounds capacity_matrix[target, super_t] = typemax(Int64)

    S = falses(n)
    @inbounds S[[super_s, source]] .= 1
    T = falses(n)
    @inbounds T[[super_t, target]] .= 1

    S_reachable = copy(S)
    T_reachable = copy(T)

    cuts = Vector{Cut}()

    forward_grow!(S_reachable, g, flow_matrix, capacity_matrix)
    forward_grow!(T_reachable, g, flow_matrix, capacity_matrix)

    while (!any(S .& T)) && (sum(S .| T) < n)
        if any(S_reachable .& T_reachable)
            augment_flow!(flow_matrix, capacity_matrix, g, super_s, super_t)
            S_reachable = copy(S)
            T_reachable = copy(T)
            forward_grow!(S_reachable, g, flow_matrix, capacity_matrix)
            forward_grow!(T_reachable, g, flow_matrix, capacity_matrix)
        else
            cut_arcs = Vector{Pair{Int64,Int64}}()
            size_SR = sum(S_reachable)
            size_TR = sum(T_reachable)
            if size_SR <= size_TR
                forward_grow!(S, g, flow_matrix, capacity_matrix)
                # output source side cut edges
                for e in edges(g)
                    # keep arcs such as in u => v u is always from the larger side of the cut
                    @inbounds if S_reachable[e.src] && (!S_reachable[e.dst])
                        push!(cut_arcs, e.dst => e.src)
                    elseif (!S_reachable[e.src]) && S_reachable[e.dst]
                        push!(cut_arcs, e.src => e.dst)
                    end
                end
                push!(
                    cuts,
                    Cut(
                        arcs = cut_arcs,
                        imbalance = 2 * (n - size_SR - 1) / (n - 2) - 1,
                        expansion = length(cut_arcs) / (size_SR - 1)
                    ),
                )

                x = piercing_node(
                    cut_arcs,
                    S_reachable,
                    T_reachable,
                    true,
                    dist,
                )
                @inbounds S[x] = 1
                @inbounds S_reachable[x] = 1
                add_edge!(g, super_s, x)
                @inbounds capacity_matrix[super_s, x] = typemax(Int64)

                forward_grow!(S_reachable, g, flow_matrix, capacity_matrix)
            else
                forward_grow!(T, g, flow_matrix, capacity_matrix)
                # output target side cut edges

                for e in edges(g)
                    # keep arcs such as in u => v u is always from the larger side of the cut
                    @inbounds if T_reachable[e.src] && (!T_reachable[e.dst])
                        push!(cut_arcs, e.dst => e.src)
                    elseif (!T_reachable[e.src]) && T_reachable[e.dst]
                        push!(cut_arcs, e.src => e.dst)
                    end
                end
                push!(
                    cuts,
                    Cut(
                        arcs = cut_arcs,
                        imbalance =  2 * (n - size_TR - 1) / (n - 2) - 1,
                        expansion = length(cut_arcs) / (size_TR - 1)
                    )
                )

                x = piercing_node(
                    cut_arcs,
                    T_reachable,
                    S_reachable,
                    false,
                    dist,
                )
                @inbounds T[x] = 1
                @inbounds T_reachable[x] = 1
                add_edge!(g, x, super_t)
                @inbounds capacity_matrix[x, super_t] = typemax(Int64)

                forward_grow!(
                    T_reachable,
                    g,
                    flow_matrix,
                    capacity_matrix,
                )
            end
        end
    end
    rem_vertex!(g, super_t)
    rem_vertex!(g, super_s)
    return cuts
end


text:
# Translation from our variables to Nocedal and Wright's
# JMW's dx <=> NW's s
# JMW's dg <=> NW' y

macro bfgstrace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(g)
                dt["~inv(H)"] = copy(invH)
            end
            g_norm = vecnorm(g, Inf)
            update!(tr,
                    iteration,
                    f_x,
                    g_norm,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

immutable BFGS <: Optimizer
    linesearch!::Function
end

BFGS(; linesearch!::Function = hz_linesearch!) =
  BFGS(linesearch!)

function optimize{T}(d::DifferentiableFunction,
                     initial_x::Vector{T},
                     mo::BFGS,
                     o::OptimizationOptions;
                     initial_invH::Matrix = eye(length(initial_x)))
    # Print header if show_trace is set
    print_header(o)

    # Maintain current state in x and previous state in x_previous
    x, x_previous = copy(initial_x), copy(initial_x)

    # Count the total number of iterations
    iteration = 0

    # Track calls to function and gradient
    f_calls, g_calls = 0, 0

    # Count number of parameters
    n = length(x)

    # Maintain current gradient in g and previous gradient in g_previous
    g, g_previous = Array(T, n), Array(T, n)

    # Store the approximate inverse Hessian in invH
    invH = copy(initial_invH)

    # The current search direction
    s = Array(T, n)

    # Intermediate storage
    u = Array(T, n)

    # Buffers for use in line search
    x_ls, g_ls = Array(T, n), Array(T, n)

    # Store f(x) in f_x
    f_x_previous, f_x = NaN, d.fg!(x, g)
    f_calls, g_calls = f_calls + 1, g_calls + 1
    copy!(g_previous, g)

    # Keep track of step-sizes
    alpha = alphainit(one(T), x, g, f_x)

    # TODO: How should this flag be set?
    mayterminate = false

    # Maintain a cache for line search results
    lsr = LineSearchResults(T)

    # Maintain record of changes in position and gradient
    dx, dg = Array(T, n), Array(T, n)

    # Maintain a cached copy of the identity matrix
    I = eye(size(invH)...)

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.store_trace || o.show_trace || o.extended_trace || o.callback != nothing
    @bfgstrace

    # Assess multiple types of convergence
    x_converged, f_converged = false, false
    g_converged = vecnorm(g, Inf) < o.g_tol

    # Iterate until convergence
    converged = g_converged
    while !converged && iteration < o.iterations
        # Increment the number of steps we've had to perform
        iteration += 1

        # Set the search direction
        # Search direction is the negative gradient divided by the approximate Hessian
        A_mul_B!(s, invH, g)
        scale!(s, -1)

        # Refresh the line search cache
        dphi0 = vecdot(g, s)
        # If invH is not positive definite, reset it to I
        if dphi0 > 0.0
            copy!(invH, I)
            @simd for i in 1:n
                @inbounds s[i] = -g[i]
            end
            dphi0 = vecdot(g, s)
        end
        clear!(lsr)
        push!(lsr, zero(T), f_x, dphi0)

        # Determine the distance of movement along the search line
        alpha, f_update, g_update =
          mo.linesearch!(d, x, s, x_ls, g_ls, lsr, alpha, mayterminate)
        f_calls, g_calls = f_calls + f_update, g_calls + g_update

        # Maintain a record of previous position
        copy!(x_previous, x)

        # Update current position
        @simd for i in 1:n
            @inbounds dx[i] = alpha * s[i]
            @inbounds x[i] = x[i] + dx[i]
        end

        # Maintain a record of the previous gradient
        copy!(g_previous, g)

        # Update the function value and gradient
        f_x_previous, f_x = f_x, d.fg!(x, g)
        f_calls, g_calls = f_calls + 1, g_calls + 1

        x_converged,
        f_converged,
        g_converged,
        converged = assess_convergence(x,
                                       x_previous,
                                       f_x,
                                       f_x_previous,
                                       g,
                                       o.x_tol,
                                       o.f_tol,
                                       o.g_tol)

        # Measure the change in the gradient
        @simd for i in 1:n
            @inbounds dg[i] = g[i] - g_previous[i]
        end

        # Update the inverse Hessian approximation using Sherman-Morrison
        dx_dg = vecdot(dx, dg)
        if dx_dg == 0.0
            break
        end
        A_mul_B!(u, invH, dg)

        c1 = (dx_dg + vecdot(dg, u)) / (dx_dg * dx_dg)
        c2 = 1 / dx_dg

        # invH = invH + c1 * (s * s') - c2 * (u * s' + s * u')
        for i in 1:n
            @simd for j in 1:n
                @inbounds invH[i, j] += c1 * dx[i] * dx[j] - c2 * (u[i] * dx[j] + u[j] * dx[i])
            end
        end

        @bfgstrace
    end

    return MultivariateOptimizationResults("BFGS",
                                           initial_x,
                                           x,
                                           Float64(f_x),
                                           iteration,
                                           iteration == o.iterations,
                                           x_converged,
                                           o.x_tol,
                                           f_converged,
                                           o.f_tol,
                                           g_converged,
                                           o.g_tol,
                                           tr,
                                           f_calls,
                                           g_calls)
end

ext: jl


text:
# Autogenerated wrapper script for DGGRID7_jll for armv7l-linux-musleabihf-cxx03
export dggrid

using GDAL_jll
JLLWrappers.@generate_wrapper_header("DGGRID7")
JLLWrappers.@declare_executable_product(dggrid)
function __init__()
    JLLWrappers.@generate_init_header(GDAL_jll)
    JLLWrappers.@init_executable_product(
        dggrid,
        "bin/dggrid",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
This example uses the Housing dataset from the UCI Machine Learning
Repository to demonstrate a linear regression model. The dataset has
housing related information for 506 neighborhoods in Boston from
1978. Each neighborhood has 14 attributes, the goal is to use the
first 13, such as average number of rooms per house, or distance to
employment centers, to predict the 14’th attribute: median dollar
value of the houses.

To run the demo, simply `include("housing.jl")` and run `Housing.train()`.
The dataset will be automatically downloaded.  You can provide the
initial weights as an optional argument, which should be a pair of
1x13 weight matrix and a scalar bias.  `train` also accepts the
following keyword arguments: `lr` specifies the learning rate,
`epochs` gives number of epochs, and `seed` specifies the random
number seed.  The quadratic loss for the train and test sets will be
printed at every epoch and optimized parameters will be returned.
"""
module Housing
using AutoGrad
using Main

function loss(wb, x=xtrn, ygold=ytrn)
    (w,b) = wb
    ypred = w*x .+ b
    sum(abs2(ypred - ygold)) / size(ygold,2)
end

function train(w=Any[0.1*randn(1,13), 0.0]; lr=.1, epochs=20, seed=nothing)
    seed==nothing || (srand(seed); loaddata())
    isdefined(Housing,:xtrn) || loaddata()
    gradfun = grad(loss)
    println((0, :trnloss, loss(w,xtrn,ytrn), :tstloss, loss(w,xtst,ytst)))
    for epoch=1:epochs
        g = gradfun(w, xtrn, ytrn)
        for i in 1:length(w)
            w[i] -= lr * g[i]
        end
        println((epoch, :trnloss, loss(w,xtrn,ytrn), :tstloss, loss(w,xtst,ytst)))
    end
    return w
end

function loaddata()
    global xtrn, ytrn, xtst, ytst
    url = "https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data"
    file=joinpath(AutoGrad.datapath, "housing.data")
    if !isfile(file)
        info("Downloading $url to $file")
        download(url, file)
    end
    data = readdlm(file)'
    @show size(data) # (14,506)
    x = data[1:13,:]
    y = data[14:14,:]
    x = (x .- mean(x,2)) ./ std(x,2) # Data normalization
    r = randperm(size(x,2))          # trn/tst split
    xtrn=x[:,r[1:400]]
    ytrn=y[:,r[1:400]]
    xtst=x[:,r[401:end]]
    ytst=y[:,r[401:end]]
end

end # module Housing


text:
module collision_2_boxes

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of SolidFileMesh
vmat2 = deepcopy(vmat1)                                           # material of convex decomposition of SolidFileMesh
vmat2.transparency = 0.7

cmat = "Steel"

@assembly ThreeDFiles begin
  world = Modia3D.Object3D(visualizeFrame=true)

  boxMoving     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(0.5,0.5,0.5) , "Aluminium", vmat1; contactMaterial = cmat); r=[0.0, 0.0, 5.0], fixed=false) # ,R=ModiaMath.rot2(-pi/3),

  box     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(3.0,2.0,5.0) , "Aluminium", vmat1; contactMaterial = cmat); r=[-3.0, 0.0, 2.5], fixed=true) # R=ModiaMath.rot2(-pi/3),
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[-1,0,0])
threeD = ThreeDFiles(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=true))

# Modia3D.visualizeAssembly!( threeD )

model = Modia3D.SimulationModel( threeD )
result = ModiaMath.simulate!(model; stopTime=0.71, tolerance=1e-6,interval=0.001, log=false)

println("... success of collision_2_boxes.jl!")

end


text:
using SIMD

StaticNumbers.@glue_to SIMD

@testset "SIMD - Vec generator" begin
    @test Vec(i^2 for i ∈ static(1:3)) === Vec((1, 4, 9))
end


text:
"""
A script checking if the oscilatory assemblies are a chaotic attractor.

We see the exponential increase in error even though the initial condition
is only one neuron different.

We see about the same distribution of visitations to each neuron.

Works without making assemblies in the areas! It is just about the assembly
input creating an attractor.

When the assemblies are projected, the attractor actually oscillated through the
assemblies, but this is probably because some of the neurons incorperated into
the projected assemblies already belonged to the attractor and this is why they
were incorperated into the assemblies in the first place.

We are just watching an attractor move through space.

Without an input stimulus, the orbit appears to wander, but with a stimulus,
this wandering is constrained to visit the same nodes with the same frequencies.

Amazing!

Questions:
1. Can we predict the attractor frequencies via undriven network frequencies and
   the stimulus vector? All the stimulus vector does is increase the probability
   that the given neurons will win.
2. If used in combination with the network frequencies will the stimulus vector
   give us the probabilities that the neuron ends up in the attractor?
"""

using AssemblyCalculus
using AssemblyCalculus: overlap, freeze_synapses!, synapse_noise!
using Plots
using StatsBase

K = 20
N = K^2
MEAN_DEGREE = K
P = MEAN_DEGREE / N

MAX_ITERS = 1000
TOL = floor(Int, 0.95 * K)

brain = Brain(num_areas=3, n=N, k=K, β=0.001, p=P)
synapse_noise!(brain)

assems, sp, ds = simulate!([rand_stim(brain[3])], max_iters=MAX_ITERS)
a0 = assems[1]

# Project alternating pattern across areas 1 and 2
# assems, sp, ds = simulate!([rand_stim(brain[1])], [a0], tol = TOL, max_iters=MAX_ITERS)
# a1 = assems[1]
# assems, sp, ds = simulate!([rand_stim(brain[2])], [a0], tol = TOL, max_iters=MAX_ITERS)
# a2 = assems[1]
# assems, sp, ds = simulate!([rand_stim(brain[1])], [a0], tol = TOL, max_iters=MAX_ITERS)
# a3 = assems[1]
# assems, sp, ds = simulate!([rand_stim(brain[2])], [a0], tol = TOL, max_iters=MAX_ITERS)
# a4 = assems[1]

# println("overlap a1, a3: $(overlap(a1, a3))")
# println("overlap a2, a4: $(overlap(a2, a4))")

freeze_synapses!(brain)
map(random_firing!, brain.areas[1:2])
init_fire = deepcopy([brain[1].firing, brain[2].firing])

timesteps = 500

assems, sp, ds = simulate!(
   [zero_stim(brain[1]), zero_stim(brain[2])],
   [a0],
   max_iters=timesteps,
   tol=0,
   random_initial=false,
)

delta_init = init_fire
delta_init[1][1] = rand(1:N)
brain[1].firing = delta_init[1]
brain[2].firing = delta_init[2]

assems, sp_delta, ds_delta = simulate!(
   [zero_stim(brain[1]), zero_stim(brain[2])],
   [a0],
   max_iters=timesteps,
   tol=0,
   random_initial=false,
)

function spike_to_vec(firing::Array{Int, 1}, n::Int)
   x = zeros(n)
   x[firing] .= 1.0
   return x
end

function spike_to_vec(spikes::Array{Array{Int, 2}, N}, n::Int) where N
   mapreduce(A -> mapslices(x -> spike_to_vec(x, n), A, dims=1), vcat, spikes)
end

num_diff(x) = sum(x[x .> 0])
V = spike_to_vec(sp, N)
Vdelta = spike_to_vec(sp_delta, N)
diff = abs.(V .- Vdelta)
shared = V .* Vdelta

err = reshape(mapslices(num_diff, V - Vdelta, dims=1), :)

p1 = plot()
plot!(
   p1,
   err,
   linewidth=3,
   color=:purple,
   label="Orbit error",
   legend=:bottomright
)
xlabel!("Timestep")
ylabel!("Difference in Assembly")

p2 = plot()
histogram!(p2, reshape(sp[1], :), alpha=0.5, bins=N)
histogram!(p2, reshape(sp_delta[1], :), alpha=0.5, bins=N)

p3 = plot()
histogram!(p3, reshape(sp[2], :), alpha=0.5, bins=N)
histogram!(p3, reshape(sp_delta[2], :), alpha=0.5, bins=N)

map(x -> x.plasticity = 0.05, brain.areas)
_, sp_plastic, ds_plastic = simulate!(
   [zero_stim(brain[1]), zero_stim(brain[2])],
   [a0],
   max_iters=timesteps,
   tol=0,
   random_initial=false,
)

V_plastic = spike_to_vec(sp_plastic, N);

plot(p1, p2, p3, size=(1200, 400), layout=@layout [a{0.25w} [b; c]])


text:
module biclusterseval

export get_biclusters, evaluate_fitness, compress_chromes, initialize_input_on_gpus

using CUDA
using Base.Threads: @threads, nthreads, threadid
using CSV: File
using DataFrames: DataFrame, groupby, Not, select

include("constants.jl")

function evaluate_fitness(
    fitness,
    input_data,
    rows_number::UInt32,
    cchromes, # compressed chromes
    cids, # compressed chromes indices
)::Nothing
    idx_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x # bicluster/chromo number
    idx_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y # row number

    trend_check = @cuStaticSharedMem(Int32, BLOCK_SIZE)
    trend_check[threadIdx().y] = 0

    idx_y > rows_number && return nothing

    evaluate_trends(trend_check, input_data, cchromes, cids)

    if BLOCK_SIZE == 1024
        if threadIdx().y <= 512 && idx_y + 512 <= rows_number
            trend_check[threadIdx().y] += trend_check[threadIdx().y + 512]
        end
        sync_threads()
    end

    if BLOCK_SIZE >= 512
        if threadIdx().y <= 256 && idx_y + 256 <= rows_number
            trend_check[threadIdx().y] += trend_check[threadIdx().y + 256]
        end
        sync_threads()
    end

    if BLOCK_SIZE >= 256
        if threadIdx().y <= 128 && idx_y + 128 <= rows_number
            trend_check[threadIdx().y] += trend_check[threadIdx().y + 128]
        end
        sync_threads()
    end

    if BLOCK_SIZE >= 128
        if threadIdx().y <= 64 && idx_y + 64 <= rows_number
            trend_check[threadIdx().y] += trend_check[threadIdx().y + 64]
        end
        sync_threads()
    end

    if BLOCK_SIZE >= 64
        if threadIdx().y <= 32 && idx_y + 32 <= rows_number
            trend_check[threadIdx().y] += trend_check[threadIdx().y + 32]
        end
        sync_threads()
    end

    if threadIdx().y <= 32
        trend_check[threadIdx().y] += trend_check[threadIdx().y + 16]
        sync_warp()
        trend_check[threadIdx().y] += trend_check[threadIdx().y + 8]
        sync_warp()
        trend_check[threadIdx().y] += trend_check[threadIdx().y + 4]
        sync_warp()
        trend_check[threadIdx().y] += trend_check[threadIdx().y + 2]
        sync_warp()
        trend_check[threadIdx().y] += trend_check[threadIdx().y + 1]

        if threadIdx().y == 1
            @atomic fitness[idx_x] += trend_check[1]
        end
    end

    return nothing
end

function get_biclusters(
    d_input_data::Vector{CuArray{Float32,2}},
    population::Population,
    gpus_num::Int,
    negative_trends,
    approx_trends_ratio,
)
    compressed_chromes, chromes_ids = compress_chromes(population)

    matrices = Vector(undef, nthreads())
    @threads for (dev, d_data_subset) in collect(zip(devices(), d_input_data))
        device!(dev)

        rows_number = size(d_data_subset, 1)

        d_matrix = CUDA.zeros(Int32, (rows_number, length(population)))

        blocks_per_chromo = ceil(Int, rows_number / BLOCK_SIZE)

        d_compressed_chromes = CuArray{Int32}(undef, length(compressed_chromes))
        copyto!(d_compressed_chromes, compressed_chromes)
        d_chromes_ids = CuArray{Int32}(undef, length(chromes_ids))
        copyto!(d_chromes_ids, chromes_ids)

        @cuda blocks = (length(population), blocks_per_chromo) threads = (1, BLOCK_SIZE) get_biclusters_rows(
            d_matrix,
            d_data_subset,
            rows_number,
            d_compressed_chromes,
            d_chromes_ids,
            negative_trends,
            approx_trends_ratio,
        )

        synchronize()

        matrices[threadid()] = Array(d_matrix)

        CUDA.unsafe_free!(d_matrix)
        CUDA.unsafe_free!(d_data_subset)
        CUDA.unsafe_free!(d_compressed_chromes)
        CUDA.unsafe_free!(d_chromes_ids)
    end

    matrix = vcat(matrices...)

    return [
        Dict("cols" => chromo, "rows" => findall(isone, matrix[:, i]))
        for (i, chromo) in enumerate(population)
    ]
end

function get_biclusters_rows(
    rows_matrix,
    input_data,
    rows_number::Int,
    cchromes, # compressed chromes
    cids, # compressed chromes indices
    negative_trends,
    approx_trends_ratio,
)::Nothing
    idx_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x # bicluster/chromo number
    idx_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y # row number

    trend_check = @cuStaticSharedMem(Int32, BLOCK_SIZE)
    trend_check[threadIdx().y] = 0

    idx_y > rows_number && return nothing

    evaluate_trends(
        trend_check,
        input_data,
        cchromes,
        cids,
        approx_trends_ratio = approx_trends_ratio,
        trend_sign = 1,
    )

    if negative_trends
        evaluate_trends(
            trend_check,
            input_data,
            cchromes,
            cids,
            approx_trends_ratio = approx_trends_ratio,
            trend_sign = -1,
        )
    end

    rows_matrix[idx_y, idx_x] = trend_check[threadIdx().y]

    return nothing
end

function evaluate_trends(
    trend_check,
    input_data,
    cchromes,
    cids;
    approx_trends_ratio = 1,
    trend_sign = 1,
)::Nothing
    idx_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x # bicluster/chromo number
    idx_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y # row number

    prev_value = input_data[idx_y, cchromes[cids[idx_x]]]

    trend_count::Int32 = 0
    for i = (cids[idx_x] + 1):(cids[idx_x + 1] - 1)
        next_value = input_data[idx_y, cchromes[i]]

        trend_count +=
            trend_sign * (next_value - prev_value + EPSILON) >= 0 &&
            prev_value != typemax(Float32)

        prev_value = next_value
    end

    chromo_len = cids[idx_x + 1] - cids[idx_x]
    trend_check[threadIdx().y] += trend_count + 1 >= chromo_len * approx_trends_ratio

    sync_threads()
end

function initialize_input_on_gpus(
    input_path::String,
    gpus_num::Int,
)::Vector{CuArray{Float32,2}}
    length(devices()) < gpus_num &&
        error("Not enough GPUs available: $(length(devices())) < $(gpus_num).")
    nthreads() == gpus_num ||
        error("The number of GPUs and threads must be equal, use '-t' option.")

    data = DataFrame(File(input_path))
    data = data[!, 2:end]
    data = coalesce.(data, typemax(Float32))

    nrows = size(data, 1)
    data.gpu_no = repeat(1:gpus_num, inner = ceil(Int, nrows / gpus_num))[1:nrows]

    d_input_data = Vector(undef, gpus_num)
    @threads for (dev, data_subset) in collect(zip(devices(), groupby(data, :gpu_no)))
        device!(dev)
        d_input_data[threadid()] = CuArray(Matrix(select(data_subset, Not(:gpu_no))))
    end

    return d_input_data
end

function compress_chromes(population::Population)
    compressed_chromes = Vector{Int32}()
    chromes_ids = Vector{Int32}()

    for chromo in population
        push!(chromes_ids, length(compressed_chromes) + 1)
        append!(compressed_chromes, chromo)
    end
    push!(chromes_ids, length(compressed_chromes) + 1)

    return compressed_chromes, chromes_ids
end

end


text:
include("candlesticks.jl")
include("levels.jl")
include("movingaverages.jl")
include("momentum.jl")
include("volatility.jl")
include("volume.jl")


text:
using Documenter
using OpticSimVis

makedocs(
    sitename = "OpticSimVis.jl",
    format = Documenter.HTML(
        # prettyurls = get(ENV, "CI", nothing) == "true",
        # assets = [asset("assets/GlimmerJulia.png", class = :ico, islocal = true)],
    ),
    modules = [Glimmer],
    pages = [
        "Home" => "index.md",
        "Reference" => "ref.md",
    ],
    expandfirst = [])

deploydocs(
    repo = "github.com/galran/OpticSimVis.jl.git",
    branch = "gh-pages",
    devbranch = "main",
    push_preview = true,
)


text:
@testset "accumulation.jl" begin
    @testset "scalar" begin
        @test 16 == add!!(12, 4)
    end

    @testset "Differentials" begin
        @test 16 == add!!(12, @thunk(2*2))
        @test 16 == add!!(16, Zero())

        @test 16 == add!!(16, DoesNotExist())  # Should this be an error?
    end

    @testset "Array" begin
        @testset "Happy Path" begin
            @testset "RHS Array" begin
                A = [1.0 2.0; 3.0 4.0]
                result = -1.0*ones(2,2)
                ret = add!!(result, A)
                @test ret === result  # must be same object
                @test result == [0.0 1.0; 2.0 3.0]
            end

            @testset "RHS StaticArray" begin
                A = @SMatrix[1.0 2.0; 3.0 4.0]
                result = -1.0*ones(2,2)
                ret = add!!(result, A)
                @test ret === result  # must be same object
                @test result == [0.0 1.0; 2.0 3.0]
            end

            @testset "RHS Diagonal" begin
                A = Diagonal([1.0, 2.0])
                result = -1.0*ones(2,2)
                ret = add!!(result, A)
                @test ret === result  # must be same object
                @test result == [0.0 -1.0; -1.0 1.0]
            end
        end

        @testset "Unhappy Path" begin
            # wrong length
            @test_throws DimensionMismatch add!!(ones(4,4), ones(2,2))
            # wrong shape
            @test_throws DimensionMismatch add!!(ones(4,4), ones(16))
            # wrong type (adding scalar to array)
            @test_throws MethodError add!!(ones(4), 21.0)
        end
    end

    @testset "InplaceableThunk" begin
        A=[1.0 2.0; 3.0 4.0]
        ithunk = InplaceableThunk(
            @thunk(A*B),
            x -> x.+=A
        )

        result = -1.0*ones(2,2)
        ret = add!!(result, ithunk)
        @test ret === result  # must be same object
        @test result == [0.0 1.0; 2.0 3.0]
    end
end


text:
"""
    mutable struct Structure_C  <: Str

This contains information for model constructing of C and its compounds.

# Supported list:
- structure_graphene()
- structure_graphene_ort()
"""
mutable struct Structure_C <: Str
    atom_vec::Matrix{Float64}
    cell_vec::Matrix{Float64}
    atom_type
    atom_name
    atom_charge
    para_mass
    num_atoms
    num_atom_types
    vec_type_id
    res_name
end

"""
    structure_graphene()

Do this will generate a Structure_C type which contains all information needed to build a Graphene moedel.
"""
function structure_graphene()
    atom_vec = [
        1/3 1/3 0
        2/3 2/3 0
    ]
    cell_vec = [
        2.46 0                0
        1.23 2.46 * sin(pi/3) 0
        0    0                3.35
    ]
    atom_type = [1 1]
    atom_charge = 0 .* atom_type
    para_mass = [12]
    atom_name = split("C")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "Graph"
    Structure_C(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1], res_name)
end

"""
    structure_graphene_ort()

Do this will generate a Famliy_C type which contains all information needed to build a Graphene moedel with orthogonal unit cell.
"""
function structure_graphene_ort()
    atom_vec = [
        0   0   0
        0   1/3 0
        1/2 1/2 0
        1/2 5/6 0
    ]
    cell_vec = [
        2.46 0                0
        0    3.69 / cos(pi/6) 0
        0    0                0.35
    ]
    atom_type = [1 1 1 1]
    atom_charge = 0 .* atom_type
    para_mass = [12]
    atom_name = split("C")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "Graph"
    Structure_C(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1], res_name)
end


text:
import Base: length, map
import Base: @propagate_inbounds
import Rocket: tap

## Equality node is a special case and has a special implementation
## It should not be used during model creation but instead is a part of variable node implementation

abstract type EqualityNodeOutboundType end

struct EqualityLeftOutbound <: EqualityNodeOutboundType end
struct EqualityRightOutbound <: EqualityNodeOutboundType end

"""
    EqualityNode

Generic idea of an equality node is to keep track of intermediate `left` and `right` messages
in a convenient manner so they can be reused during `prod` for random variables with large `degree`.

 <- left    -------     right ->
     ----- |   =   | ------
            -------
               |
               | outbound

# Attributes:
- `left`: An observable that indicates that `left` message can be computed, always sends `missing` since we are not interested in actual value. Can be treated as notification.
- `right`: An observable that indicates that `right` message can be computed, always sends `missing` since we are not interested in actual value. Can be treated as notification.
- `cache_left`: Keeps cached value for left outbound message, node itself does not track if cache is valid or not. Cache integrity is tracked by `EqualityChain`.
- `cache_right`: Keeps cached value for right outbound message, node itself does not track if cache is valid or not. Cache integrity is tracked by `EqualityChain`.

See also: [`EqualityChain`](@ref)
"""
mutable struct EqualityNode
    left         :: LazyObservable{Missing}
    right        :: LazyObservable{Missing}
    cache_left   :: Message
    cache_right  :: Message

    EqualityNode() = new(lazy(Missing), lazy(Missing), Message(missing, true, true), Message(missing, true, true))
end

getoutbound(::EqualityLeftOutbound, node::EqualityNode)  = node.left
getoutbound(::EqualityRightOutbound, node::EqualityNode) = node.right

setoutbound!(::EqualityLeftOutbound, node::EqualityNode, left)   = set!(node.left, left)
setoutbound!(::EqualityRightOutbound, node::EqualityNode, right) = set!(node.right, right)

getcache(::EqualityLeftOutbound, node::EqualityNode)  = node.cache_left
getcache(::EqualityRightOutbound, node::EqualityNode) = node.cache_right

setcache!(::EqualityLeftOutbound, node::EqualityNode, cache::Message)  = node.cache_left = cache
setcache!(::EqualityRightOutbound, node::EqualityNode, cache::Message) = node.cache_right = cache

# Equality chain

"""
    EqualityChain
"""
struct EqualityChain{P, F}
    length     :: Int
    nodes      :: Vector{EqualityNode}
    inputmsgs  :: Vector{MessageObservable{AbstractMessage}}
    cacheleft  :: BitVector
    cacheright :: BitVector
    pipeline   :: P
    prod_fn    :: F

    function EqualityChain(inputmsgs::Vector{MessageObservable{AbstractMessage}}, pipeline::P, prod_fn::F) where { P, F }
        n     = length(inputmsgs)
        nodes = map(_ -> EqualityNode(), 1:n)
        return new{P, F}(n, nodes, inputmsgs, falses(n), falses(n), pipeline, prod_fn)
    end
end

Base.length(chain::EqualityChain) = chain.length

prod(chain::EqualityChain, left, right) = chain.prod_fn((left, right))

getpipeline(chain::EqualityChain) = chain.pipeline

@propagate_inbounds getnode(chain::EqualityChain, node_index) = chain.nodes[node_index]

__check_indices(::EqualityLeftOutbound, chain::EqualityChain, node_index)  = 1 < node_index <= length(chain)
__check_indices(::EqualityRightOutbound, chain::EqualityChain, node_index) = 1 <= node_index < length(chain)

@propagate_inbounds getoutbound(type::EqualityNodeOutboundType, chain::EqualityChain, node_index) = __check_indices(type, chain, node_index) ? getoutbound(type, getnode(chain, node_index)) : (of(missing))
@propagate_inbounds getinbound(chain::EqualityChain, node_index)                                  = chain.inputmsgs[node_index]

@propagate_inbounds iscached(::EqualityLeftOutbound, chain::EqualityChain, node_index)  = chain.cacheleft[node_index]
@propagate_inbounds iscached(::EqualityRightOutbound, chain::EqualityChain, node_index) = chain.cacheright[node_index]

@propagate_inbounds setcache!(::EqualityLeftOutbound, chain::EqualityChain, node_index)  = chain.cacheleft[node_index] = true
@propagate_inbounds setcache!(::EqualityRightOutbound, chain::EqualityChain, node_index) = chain.cacheright[node_index] = true

@propagate_inbounds setcache!(::EqualityLeftOutbound, chain::EqualityChain, range::OrdinalRange)  = fill_bitarray!(view(chain.cacheleft, forward_range(range)), true)
@propagate_inbounds setcache!(::EqualityRightOutbound, chain::EqualityChain, range::OrdinalRange) = fill_bitarray!(view(chain.cacheright, forward_range(range)), true)

@propagate_inbounds function getcache(type::EqualityNodeOutboundType, chain::EqualityChain, node_index)
    if __check_indices(type, chain, node_index)
        return getcache(type, getnode(chain, node_index))
    else
        return Message(missing, true, true)
    end
end

nextindex(::EqualityLeftOutbound, node_index)  = node_index + 1
nextindex(::EqualityRightOutbound, node_index) = node_index - 1

@propagate_inbounds first_unmaterialized_index(::EqualityLeftOutbound, chain::EqualityChain, node_index)::Int  = default_if_nothing(findfirst(view(chain.cacheleft, node_index:length(chain))), length(chain))
@propagate_inbounds first_unmaterialized_index(::EqualityRightOutbound, chain::EqualityChain, node_index)::Int = default_if_nothing(findlast(view(chain.cacheright, 1:node_index)), 1)

@propagate_inbounds precompute_range(type::EqualityLeftOutbound, chain::EqualityChain, node_index)  = first_unmaterialized_index(type, chain, node_index):-1:node_index
@propagate_inbounds precompute_range(type::EqualityRightOutbound, chain::EqualityChain, node_index) = first_unmaterialized_index(type, chain, node_index):node_index

@propagate_inbounds function materialize!(type::EqualityNodeOutboundType, chain::EqualityChain, node_index)
    if __check_indices(type, chain, node_index)
        node = getnode(chain, node_index)
        if iscached(type, chain, node_index)
            return getcache(type, node)
        else
            # precompute messages in linear fashion
            range = precompute_range(type, chain, node_index)
            for index in range
                arg1 = as_message(getrecent(getinbound(chain, index)))
                arg2 = as_message(getcache(type, chain, nextindex(type, index)))
                result = prod(chain, arg1, arg2)
                setcache!(type, getnode(chain, index), result)
            end
            setcache!(type, chain, range)
            return materialize!(type, chain, node_index)
        end
    else
        return Message(missing, true, true)
    end
end

##

struct ChainInvalidationCallback
    index      :: Int
    cacheleft  :: BitVector
    cacheright :: BitVector

    function ChainInvalidationCallback(index::Int, chain::EqualityChain)
        return new(index, chain.cacheleft, chain.cacheright)
    end
end

Rocket.tap(callback::ChainInvalidationCallback) = Rocket.TapOperator{ChainInvalidationCallback}(callback)

function (callback::ChainInvalidationCallback)(_)
    fill_bitarray!(view(callback.cacheleft, firstindex(callback.cacheleft):(callback.index)), false)
    fill_bitarray!(view(callback.cacheright, (callback.index):lastindex(callback.cacheright)), false)
end

##

struct ChainOutboundMapping
    index :: Int
    chain :: EqualityChain
end

function (mapping::ChainOutboundMapping)(_)
    from_left  = materialize!(EqualityRightOutbound(), mapping.chain, nextindex(EqualityRightOutbound(), mapping.index))
    from_right = materialize!(EqualityLeftOutbound(), mapping.chain, nextindex(EqualityLeftOutbound(), mapping.index))
    return as_message(prod(mapping.chain, from_left, from_right))
end

Base.map(::Type{Message}, mapping::ChainOutboundMapping) = Rocket.MapOperator{Message, ChainOutboundMapping}(mapping)

function initialize!(chain::EqualityChain, outputmsgs::AbstractVector)
    n = length(chain)

    pipeline = getpipeline(chain)

    Left  = EqualityLeftOutbound()
    Right = EqualityRightOutbound()

    @inbounds for index in 1:n
        node = getnode(chain, index)

        # As soon as we receive new inbound message - we invalidate cache for part of the chain: see ChainInvalidationCallback
        input = getinbound(chain, index) |> tap(ChainInvalidationCallback(index, chain)) |> share_recent()

        left  = combineLatestUpdates((getoutbound(Left, chain, nextindex(Left, index)), input), PushNew()) |> pipeline |> map_to(missing) |> share_recent()
        right = combineLatestUpdates((getoutbound(Right, chain, nextindex(Right, index)), input), PushNew()) |> pipeline |> map_to(missing) |> share_recent()

        setoutbound!(Left, node, left)
        setoutbound!(Right, node, right)

        from_left  = getoutbound(Right, chain, nextindex(Right, index)) # Inbound message comming from left direction  (is a right from `index - 1`)
        from_right = getoutbound(Left, chain, nextindex(Left, index))  # Inbound message comming from right direction (is a left from `index + 1`)

        connect!(outputmsgs[index], combineLatestUpdates((from_left, from_right), PushNew()) |> map(Message, ChainOutboundMapping(index, chain)))
    end

    return nothing
end


text:
function mountaincar_run_agnostic_sysid(
    rock_c::Float64,
    position_sigma::Float64;
    value_aware = true,
)
    mountaincar = MountainCar(rock_c; position_sigma = position_sigma)
    model = MountainCar(0.0)
    horizon = 500
    agent = MountainCarAgnosticSysIDAgent(
        mountaincar,
        model,
        horizon,
    )
    run(agent; debug = true, value_aware = value_aware)
end


text:
module Helpers
    function retroAffect(xs, ys, reward, horizon=10)

        cuttingIndex = max(1, length(ys) - horizon)
        correctedHorizon = min(length(ys), horizon)
        for (i, y) in enumerate(ys[cuttingIndex:end])
            ys[cuttingIndex + (i - 1)][1] = y[1] + (reward / max((correctedHorizon - (i - 1)), 1))
        end

        return (xs, ys)
    end
    export retroAffect
end

max_stars_repo_name: evrardco/neural-snake

max_issues_repo_name: evrardco/neural-snake

max_forks_repo_name: evrardco/neural-snake


text:
using ComplexMixtures
using PDBTools

# Load pdb file of the system
system = readPDB("./simulation/equilibrated.pdb")

# The trajectory file can be downloaded from:
# https://drive.google.com/file/d/12TT5tblkFp1NtFOAQgjjGhmnYaXA8vQi/view?usp=sharing
trajectory_file = "../trajectories/traj_POPC.dcd"

# Select the atoms corresponding to glycerol and water
popc = select(system,"resname POPC")
O_water = select(system,"water and element O")
O_ethanol = select(system,"resname ETOH and element O")

# Set the complete membrane as the solute. We use nmols=1 here such
# that the membrane is considered a single solute in the calculation.
solute = Selection(popc,nmols=1)

# Compute water-POPC distribution and KB integral
solvent = Selection(O_water,natomspermol=1)
traj = Trajectory(trajectory_file,solute,solvent)

# We want to get reasonably converged KB integrals, which usually
# require large solute domains. Distribution functions converge
# rapidly (~10Angs or less), on the other side.
opt = ComplexMixtures.Options(dbulk=20.,n_random_samples=2)

# Compute the MDDF
mddf_O_water_POPC = mddf(traj,opt)

# Save results for later use
save(mddf_O_water_POPC,"./results/mddf_O_water_POPC.json")

# Compute water-POPC distribution and KB integral
solvent = Selection(O_ethanol,natomspermol=1)
traj = Trajectory(trajectory_file,solute,solvent)
mddf_O_ethanol_POPC = mddf(traj,opt)

# Save results for later use
save(mddf_O_ethanol_POPC,"./results/mddf_O_ethanol_POPC.json")

C2_ethanol = select(system,"resname ETOH and name C2")
solvent = Selection(C2_ethanol,natomspermol=1)
traj = Trajectory(trajectory_file,solute,solvent)
mddf_C2_ethanol_POPC = mddf(traj,opt)
save(mddf_C2_ethanol_POPC,"./results/mddf_C2_ethanol_POPC.json")


text:
function compute_inverse_size(network)
    b0 = Base.gc_bytes()
    network_ref = PowerModels.build_ref(network)
    network_ref = network_ref[:nw][0]
    bus_idx2id, bus_id2idx = build_indices(network_ref[:bus])
    b_base = compute_B_sparse(network_ref, bus_idx2id, bus_id2idx)
    b_inv = compute_B_inverse(network_ref, bus_idx2id, bus_id2idx)
    return Base.gc_bytes() - b0
end

"maps dict ids into the range 1-to-n"
function build_indices(d::Dict)
    original_indices = sort(collect(keys(d)));
    num_indices = length(original_indices)
    reverse_indices = Dict(zip(original_indices,1:num_indices))
    return original_indices, reverse_indices
end

"computes the bus admittance matrix B, with indices given by `bus_id2idx::Dict`"
function compute_B_sparse(ref::Dict{Symbol,<:Any}, bus_idx2id::Vector{Int}, bus_id2idx::Dict{Int,Int}; inactive_branches::Set{Int}=Set(Int[]))
    num_buses = length(ref[:bus])
    I = Int64[]
    J = Int64[]
    V = Float64[]

    for (i,branch) in ref[:branch]
        if !(branch["index"] in inactive_branches)
            f_bus = bus_id2idx[ref[:bus][branch["f_bus"]]["index"]]
            t_bus = bus_id2idx[ref[:bus][branch["t_bus"]]["index"]]
            b_val = -branch["br_x"]/(branch["br_x"]^2+branch["br_r"]^2)
            push!(I, f_bus); push!(J, t_bus); push!(V,  b_val)
            push!(I, t_bus); push!(J, f_bus); push!(V,  b_val)
            push!(I, f_bus); push!(J, f_bus); push!(V, -b_val)
            push!(I, t_bus); push!(J, t_bus); push!(V, -b_val)
        end
    end

    return sparse(I,J,V)
end

"computes the bus admittance matrix B inverse, with indices given by `bus_id2idx::Dict`"
function compute_B_inverse(ref::Dict{Symbol,<:Any}, bus_idx2id::Vector{Int}, bus_id2idx::Dict{Int,Int}; inactive_branches::Set{Int}=Set(Int[]))
    num_buses = length(ref[:bus])
    B = compute_B(ref, bus_idx2id, bus_id2idx; inactive_branches=inactive_branches)

    ref_bus_id, ref_bus = first(ref[:ref_buses])
    r = bus_id2idx[ref_bus["index"]]
    nonref_indices = Int64[b for b in 1:num_buses if b != r]
    inv_B = zeros(Float64, num_buses, num_buses)
    inv_B[nonref_indices, nonref_indices] = inv(B[nonref_indices, nonref_indices])

    # matrix is dence, does not seem to help performance
    #inv_B = sparse(inv_B)
    #println(inv_B)

    return inv_B
end

"computes the bus admittance matrix B, with indices given by `bus_id2idx::Dict`"
function compute_B(ref::Dict{Symbol,<:Any}, bus_idx2id::Vector{Int}, bus_id2idx::Dict{Int,Int}; inactive_branches::Set{Int}=Set(Int[]))
    num_buses = length(ref[:bus])
    B = zeros(Float64, num_buses, num_buses)

    for (i,branch) in ref[:branch]
        if !(branch["index"] in inactive_branches)
            f_bus = bus_id2idx[ref[:bus][branch["f_bus"]]["index"]]
            t_bus = bus_id2idx[ref[:bus][branch["t_bus"]]["index"]]
            b_val = -branch["br_x"]/(branch["br_x"]^2+branch["br_r"]^2)
            B[f_bus, t_bus] += b_val
            B[t_bus, f_bus] += b_val
            B[f_bus, f_bus] += -b_val
            B[t_bus, t_bus] += -b_val
        end
    end

    return B
end

"compute bus injections"
function compute_bus_injections(ref::Dict{Symbol,<:Any}, bus_idx2id::Vector{Int})
    @assert length(ref[:bus]) == length(bus_idx2id)
    num_buses = length(bus_idx2id)
    bus_inj = [0.0 for i in 1:num_buses]
    for i in 1:num_buses
        for g in ref[:bus_gens][bus_idx2id[i]]
            bus_inj[i] += ref[:gen][g]["pg"]
        end
        for l in ref[:bus_loads][bus_idx2id[i]]
            bus_inj[i] -= ref[:load][l]["pd"]
        end
        for s in ref[:bus_shunts][bus_idx2id[i]]
            bus_inj[i] -= ref[:shunt][s]["gs"]*1.0^2
        end
    end
    return bus_inj
end

"compute bus injections"
function compute_bus_injections(ref::Dict{Symbol,<:Any}, gen_setpoint::Dict{String,<:Any}, bus_idx2id::Vector{Int})
    @assert length(ref[:bus]) == length(bus_idx2id)
    num_buses = length(bus_idx2id)
    bus_inj = [0.0 for i in 1:num_buses]
    for i in 1:num_buses
        for g in ref[:bus_gens][bus_idx2id[i]]
            bus_inj[i] += gen_setpoint["$(g)"]["pg"]
        end
        for l in ref[:bus_loads][bus_idx2id[i]]
            bus_inj[i] -= ref[:load][l]["pd"]
        end
        for s in ref[:bus_shunts][bus_idx2id[i]]
            bus_inj[i] -= ref[:shunt][s]["gs"]*1.0^2
        end
    end
    return bus_inj
end

"compute bus injections"
function compute_bus_injections(ref::Dict{Symbol,<:Any}, gen_setpoint::Dict{String,<:Any}, load_setpoint::Dict{String,<:Any}, bus_idx2id::Vector{Int}; bus_inj_offset::Float64=0.0)
    @assert length(ref[:bus]) == length(bus_idx2id)
    num_buses = length(bus_idx2id)
    bus_inj = [bus_inj_offset for i in 1:num_buses]
    for i in 1:num_buses
        for g in ref[:bus_gens][bus_idx2id[i]]
            bus_inj[i] += gen_setpoint["$(g)"]["pg"]
        end
        for l in ref[:bus_loads][bus_idx2id[i]]
            bus_inj[i] -= load_setpoint["$(l)"]["pd"]
        end
        for s in ref[:bus_shunts][bus_idx2id[i]]
            bus_inj[i] -= ref[:shunt][s]["gs"]*1.0^2
        end
    end
    return bus_inj
end

function solve_theta(B, bus_injection::Vector{Float64}, ref_bus::Int64)
    theta = B \ bus_injection
    theta = theta .- theta[ref_bus]
    return theta
end

compute_flows(ref::Dict{Symbol,<:Any}, theta, bus_id2idx::Dict{Int,Int}) = Dict(l => (-branch["br_x"]/(branch["br_x"]^2+branch["br_r"]^2))*(theta[bus_id2idx[branch["f_bus"]]] - theta[bus_id2idx[branch["t_bus"]]]) for (l,branch) in ref[:branch])

function compute_branch_ptdf_single(ref::Dict{Symbol,<:Any}, B_inverse::Matrix{Float64}, bus_id2idx::Dict{Int,Int}, branch_id::Int)
    branch_ptdf = Dict{Int,Any}()
    branch = ref[:branch][branch_id]
    row_fr = bus_id2idx[branch["f_bus"]]
    row_to = bus_id2idx[branch["t_bus"]]
    b = (-branch["br_x"]/(branch["br_x"]^2+branch["br_r"]^2))

    branch_ptdf[branch_id] = -b*(B_inverse[row_fr,:] - B_inverse[row_to,:])
    return branch_ptdf
end

function built_flow_cut(cont_label, branch_id, rating_level, branch_ptdf::Dict{Int,<:Any}, bus_idx2id::Vector{Int})
    bus_injection = Dict(bus_idx2id[idx] => value for (idx,value) in enumerate(branch_ptdf[branch_id]))
    cut = (cont_label=cont_label, branch_id=branch_id, rating_level=rating_level, bus_injection=bus_injection)
end

# a global network variable used for iterative computations
network_global = Dict{String,Any}()

# a global contingency list used for iterative computations
contingency_order_global = []

""
function load_network_global(con_file, inl_file, raw_file, rop_file, scenario_id)
    info(_LOGGER, "skipping goc and power models data warnings")
    pm_logger_level = getlevel(getlogger(PowerModels))
    goc_logger_level = getlevel(_LOGGER)

    setlevel!(getlogger(PowerModels), "error")
    setlevel!(_LOGGER, "error")

    goc_data = parse_goc_files(con_file, inl_file, raw_file, rop_file, scenario_id=scenario_id)
    global network_global = build_pm_model(goc_data)
    global contingency_order_global = contingency_order(network_global)

    setlevel!(getlogger(PowerModels), pm_logger_level)
    setlevel!(_LOGGER, goc_logger_level)

    return 0
end

"build a static ordering of all contigencies"
function contingency_order(network)
    gen_cont_order = sort(network["gen_contingencies"], by=(x) -> x.label)
    branch_cont_order = sort(network["branch_contingencies"], by=(x) -> x.label)

    gen_cont_total = length(gen_cont_order)
    branch_cont_total = length(branch_cont_order)

    gen_rate = 1.0
    branch_rate = 1.0
    steps = 1

    if gen_cont_total == 0 && branch_cont_total == 0
        # defaults are good
    elseif gen_cont_total == 0 && branch_cont_total != 0
        steps = branch_cont_total
    elseif gen_cont_total != 0 < branch_cont_total == 0
        steps = gen_cont_total
    elseif gen_cont_total == branch_cont_total
        steps = branch_cont_total
    elseif gen_cont_total < branch_cont_total
        gen_rate = 1.0
        branch_rate = branch_cont_total/gen_cont_total
        steps = gen_cont_total
    elseif gen_cont_total > branch_cont_total
        gen_rate = gen_cont_total/branch_cont_total
        branch_rate = 1.0
        steps = branch_cont_total
    end

    #println(gen_cont_total)
    #println(branch_cont_total)
    #println(steps)

    #println(gen_rate)
    #println(branch_rate)
    #println("")

    cont_order = []
    gen_cont_start = 1
    branch_cont_start = 1
    for s in 1:steps
        gen_cont_end = min(gen_cont_total, trunc(Int,ceil(s*gen_rate)))
        #println(gen_cont_start:gen_cont_end)
        for j in gen_cont_start:gen_cont_end
            push!(cont_order, gen_cont_order[j])
        end
        gen_cont_start = gen_cont_end+1

        branch_cont_end = min(branch_cont_total, trunc(Int,ceil(s*branch_rate)))
        #println("$(s) - $(branch_cont_start:branch_cont_end)")
        for j in branch_cont_start:branch_cont_end
            push!(cont_order, branch_cont_order[j])
        end
        branch_cont_start = branch_cont_end+1
    end

    #=
    for s in 1:steps
        gen_cont_start = trunc(Int, ceil(1+(s-1)*gen_rate))
        gen_cont_end = min(gen_cont_total, trunc(Int,ceil(s*gen_rate)))
        #println(gen_cont_start:gen_cont_end)
        for j in gen_cont_start:gen_cont_end
            push!(cont_order, gen_cont_order[j])
        end

        branch_cont_start = trunc(Int, ceil(1+(s-1)*branch_rate))
        branch_cont_end = min(branch_cont_total, trunc(Int,ceil(s*branch_rate)))
        println("$(s) - $(branch_cont_start:branch_cont_end)")
        for j in branch_cont_start:branch_cont_end
            push!(cont_order, branch_cont_order[j])
        end
    end
    =#

    #println(length(cont_order))
    #println(gen_cont_total + branch_cont_total)

    @assert(length(cont_order) == gen_cont_total + branch_cont_total)

    return cont_order
end

function write_contingencies(network; output_dir="", file_name="contingencies.txt")
    if length(output_dir) > 0
        file_path = joinpath(output_dir, file_name)
    else
        file_path = file_name
    end

    open(file_path, "w") do cont_file
        for c in network["gen_contingencies_active"]
            write(cont_file, "$(c.label)\n")
        end
        for c in network["branch_contingencies_active"]
            write(cont_file, "$(c.label)\n")
        end
    end
end

function write_active_flow_cuts(network; output_dir="", file_name="active_flow_cuts.txt")
    if length(output_dir) > 0
        file_path = joinpath(output_dir, file_name)
    else
        file_path = file_name
    end

    open(file_path, "w") do cont_file
        for cut in network["gen_flow_cuts"]
            write(cont_file, "$(cut.cont_label), gen, $(cut.branch_id), $(cut.rating_level)\n")
        end

        for cut in network["branch_flow_cuts"]
            write(cont_file, "$(cut.cont_label), branch, $(cut.branch_id), $(cut.rating_level)\n")
        end
    end
end

function read_active_flow_cuts(;output_dir="", file_name="active_flow_cuts.txt")
    if length(output_dir) > 0
        file_path = joinpath(output_dir, file_name)
    else
        file_path = file_name
    end

    if isfile(file_path)
        info(_LOGGER, "loading flow cuts file: $(file_path)")
        return parse_flow_cuts_file(file_path)
    else
        info(_LOGGER, "flow cuts file not found: $(file_path)")
        return []
    end
end

function parse_flow_cuts_file(file::String)
    open(file) do io
        return parse_flow_cuts_file(io)
    end
end

function parse_flow_cuts_file(io::IO)
    cuts_list = []

    for line in readlines(io)
        if length(strip(line)) == 0
            warn(_LOGGER, "skipping blank line in cuts file")
            continue
        end
        line_parts = split(line, ",")
        if length(line_parts) != 4
            warn(_LOGGER, "skipping ill formated line\n   $(line)")
            continue
        end

        cont_label = strip(line_parts[1])
        cont_type = strip(line_parts[2])
        branch_id = parse(Int, line_parts[3])
        rating_level = parse(Float64, line_parts[4])

        push!(cuts_list, (cont_label=cont_label, cont_type=cont_type, branch_id=branch_id, rating_level=rating_level))
    end

    return cuts_list
end

""
function check_contingencies_branch_power_bpv_remote(cont_range, output_dir, cut_limit=1, solution_file="solution1.txt")
    if length(network_global) <= 0 || length(contingency_order_global) <= 0
        error(_LOGGER, "check_contingencies_branch_flow_remote called before load_network_global")
    end

    sol = read_solution1(network_global, output_dir=output_dir, state_file=solution_file)
    PowerModels.update_data!(network_global, sol)

    active_cuts = read_active_flow_cuts(output_dir=output_dir)
    gen_flow_cuts = []
    branch_flow_cuts = []
    for cut in active_cuts
        if cut.cont_type == "gen"
            push!(gen_flow_cuts, cut)
        elseif cut.cont_type == "branch"
            push!(branch_flow_cuts, cut)
        else
            warn(_LOGGER, "unknown contingency type in cut $(cut)")
        end
    end

    network = copy(network_global)
    contingencies = contingency_order_global[cont_range]
    network["gen_contingencies"] = [c for c in contingencies if c.type == "gen"]
    network["branch_contingencies"] = [c for c in contingencies if c.type == "branch"]

    cuts = check_contingencies_branch_power_bpv(network, total_cut_limit=cut_limit, gen_flow_cuts=gen_flow_cuts, branch_flow_cuts=branch_flow_cuts)

    return cuts
end

"""
A variant of `check_contingencies_branch_flow_remote`, which ignores the
participation factor based generator response model from the ARPA-e GOC
Challenge 1 specification and instead injects active power equally at all buses
in the network.
"""
function check_contingencies_branch_power_bpv(network;
        gen_flow_cut_limit=10, branch_flow_cut_limit=10, total_cut_limit=typemax(Int64),
        gen_eval_limit=typemax(Int64), branch_eval_limit=typemax(Int64), sm_threshold=0.01,
        gen_flow_cuts=[], branch_flow_cuts=[]
        )

    if InfrastructureModels.ismultinetwork(network)
        error(_LOGGER, "the branch flow cut generator can only be used on single networks")
    end
    time_contingencies_start = time()

    gen_cuts_active = Dict()
    for gen_cut in gen_flow_cuts
        if !haskey(gen_cuts_active, gen_cut.cont_label)
            gen_cuts_active[gen_cut.cont_label] = Set{Int}()
        end
        push!(gen_cuts_active[gen_cut.cont_label], gen_cut.branch_id)
    end

    branch_cuts_active = Dict()
    for branch_cut in branch_flow_cuts
        if !haskey(branch_cuts_active, branch_cut.cont_label)
            branch_cuts_active[branch_cut.cont_label] = Set{Int}()
        end
        push!(branch_cuts_active[branch_cut.cont_label], branch_cut.branch_id)
    end

    network_ref = PowerModels.build_ref(network)
    network_ref = network_ref[:nw][0]

    solution_base = extract_solution(network; branch_flow=true)

    pd_total = sum(load["pd"] for (i,load) in network_ref[:load])
    p_losses = sum(gen["pg"] for (i,gen) in network_ref[:gen]) - pd_total
    load_setpoint = network["load"]
    p_delta = 0.0

    if p_losses > pg_loss_tol
        load_count = length(load_setpoint)
        p_delta = p_losses/load_count
        for (i,load) in load_setpoint
            load["pd"] = load["pd"] + p_delta
        end
        warn(_LOGGER, "active power losses found $(p_losses) increasing loads by $(p_delta)")
    end

    bus_idx2id, bus_id2idx = build_indices(network_ref[:bus])
    b_base = compute_B_sparse(network_ref, bus_idx2id, bus_id2idx)
    b_inv_base = compute_B_inverse(network_ref, bus_idx2id, bus_id2idx)
    #branch_ptdf_base = compute_branch_ptdf(network_ref, b_inv_base, bus_id2idx)
    bus_inj_base = compute_bus_injections(network_ref, solution_base["gen"], load_setpoint, bus_idx2id)

    @assert length(network_ref[:ref_buses]) == 1
    ref_bus_id, ref_bus = first(network_ref[:ref_buses])
    ref_bus_idx = bus_id2idx[ref_bus["index"]]

    gen_cont_total = length(network["gen_contingencies"])
    branch_cont_total = length(network["branch_contingencies"])

    gen_eval_limit = min(gen_eval_limit, gen_cont_total)
    branch_eval_limit = min(branch_eval_limit, branch_cont_total)

    gen_cap = Dict(gen["index"] => sqrt(max(abs(gen["pmin"]), abs(gen["pmax"]))^2 + max(abs(gen["qmin"]), abs(gen["qmax"]))^2) for (i,gen) in network["gen"])
    network["gen_contingencies"] = sort(network["gen_contingencies"], rev=true, by=x -> gen_cap[x.idx])
    gen_contingencies = network["gen_contingencies"][1:gen_eval_limit]

    line_imp_mag = Dict(branch["index"] => branch["rate_a"]*sqrt(branch["br_r"]^2 + branch["br_x"]^2) for (i,branch) in network["branch"])
    network["branch_contingencies"] = sort(network["branch_contingencies"], rev=true, by=x -> line_imp_mag[x.idx])
    branch_contingencies = network["branch_contingencies"][1:branch_eval_limit]

    gen_cuts = []
    for (i,cont) in enumerate(gen_contingencies)
        if length(gen_cuts) >= gen_flow_cut_limit
            info(_LOGGER, "hit gen flow cut limit $(gen_flow_cut_limit)")
            break
        end
        if length(gen_cuts) >= total_cut_limit
            info(_LOGGER, "hit total cut limit $(total_cut_limit)")
            break
        end
        #info(_LOGGER, "working on ($(i)/$(gen_eval_limit)/$(gen_cont_total)): $(cont.label)")

        sol_tmp = deepcopy(solution_base)
        cont_gen = sol_tmp["gen"]["$(cont.idx)"]
        pg_lost = cont_gen["pg"]
        qg_lost = cont_gen["qg"]
        cont_gen["pg"] = 0.0
        cont_gen["qg"] = 0.0

        gen = network_ref[:gen][cont.idx]
        gen_bus = network_ref[:bus][gen["gen_bus"]]

        sol_tmp["delta"] = 0.0

        bus_inj_cont = compute_bus_injections(network_ref, sol_tmp["gen"], load_setpoint, bus_idx2id, bus_inj_offset=pg_lost/length(bus_idx2id))

        va_vector = []
        try
            va_vector = solve_theta(b_base, bus_inj_cont, ref_bus_idx)
        catch exception
            warn(_LOGGER, "linear solve failed on $(cont.label)")
            continue
        end
        p_flows = compute_flows(network_ref, va_vector, bus_id2idx)

        for (l,flow) in p_flows
            sol_branch = sol_tmp["branch"]["$(l)"]
            sol_branch["pf"] =  flow
            sol_branch["pt"] = -flow
            sol_branch["qf"] = 0.0
            sol_branch["qt"] = 0.0
        end

        #PowerModels.print_summary(sol_tmp)

        network["gen"]["$(cont.idx)"]["gen_status"] = 0
        vio = compute_violations_ratec(network, sol_tmp)
        network["gen"]["$(cont.idx)"]["gen_status"] = 1

        #info(_LOGGER, "$(cont.label) violations $(vio)")
        #if vio.vm > vm_threshold || vio.pg > pg_threshold || vio.qg > qg_threshold || vio.sm > sm_threshold
        if vio.sm > sm_threshold
            branch_vio = branch_violations_sorted_ratec(network, sol_tmp)[1]
            if !haskey(gen_cuts_active, cont.label) || !(branch_vio.branch_id in gen_cuts_active[cont.label])
                info(_LOGGER, "adding flow cut on cont $(cont.label) branch $(branch_vio.branch_id) due to constraint flow violations $(branch_vio.sm_vio)")

                branch_ptdf_cont = compute_branch_ptdf_single(network_ref, b_inv_base, bus_id2idx, branch_vio.branch_id)
                cut = built_flow_cut(cont.label, branch_vio.branch_id, 1.0, branch_ptdf_cont, bus_idx2id)
                cut = (gen_id=cont.idx, cont_label=cut.cont_label, branch_id=cut.branch_id, rating_level=cut.rating_level, bus_injection=cut.bus_injection)
                push!(gen_cuts, cut)
                break
            else
                warn(_LOGGER, "skipping active flow cut on cont $(cont.label) branch $(branch_vio.branch_id) with constraint flow violations $(branch_vio.sm_vio)")
            end
        end
    end

    # work around for julia compiler bug
    #num_buses = length(network_ref[:bus])
    #b_cont = zeros(Float64, num_buses, num_buses)

    branch_cuts = []
    for (i,cont) in enumerate(branch_contingencies)
        if length(branch_cuts) >= branch_flow_cut_limit
            info(_LOGGER, "hit branch flow cut limit $(branch_flow_cut_limit)")
            break
        end
        if length(gen_cuts) + length(branch_cuts) >= total_cut_limit
            info(_LOGGER, "hit total cut limit $(total_cut_limit)")
            break
        end

        #info(_LOGGER, "working on ($(i)/$(branch_eval_limit)/$(branch_cont_total)): $(cont.label)")
        sol_tmp = deepcopy(solution_base)
        if haskey(sol_tmp, "branch")
            cont_branch = sol_tmp["branch"]["$(cont.idx)"]
            cont_branch["pf"] = 0.0
            cont_branch["pt"] = 0.0
            cont_branch["qf"] = 0.0
            cont_branch["qt"] = 0.0
        end

        sol_tmp["delta"] = 0.0

        b_cont = compute_B_sparse(network_ref, bus_idx2id, bus_id2idx, inactive_branches=Set([cont.idx]))
        # work around for julia compiler bug
        #update_B!(b_cont, network_ref, bus_idx2id, bus_id2idx, inactive_branches=Set([cont.idx]))

        va_vector = []
        try
            va_vector = solve_theta(b_cont, bus_inj_base, ref_bus_idx)
        catch exception
            warn(_LOGGER, "linear solve failed on $(cont.label)")
            continue
        end
        p_flows = compute_flows(network_ref, va_vector, bus_id2idx)

        for (l,flow) in p_flows
            sol_branch = sol_tmp["branch"]["$(l)"]
            if l != cont.idx
                sol_branch["pf"] =  flow
                sol_branch["pt"] = -flow
            else
                sol_branch["pf"] = 0.0
                sol_branch["pt"] = 0.0
            end

            sol_branch["qf"] = 0.0
            sol_branch["qt"] = 0.0
        end

        #PowerModels.print_summary(sol_tmp)

        network["branch"]["$(cont.idx)"]["br_status"] = 0
        vio = compute_violations_ratec(network, sol_tmp)
        network["branch"]["$(cont.idx)"]["br_status"] = 1

        #info(_LOGGER, "$(cont.label) violations $(vio)")
        #if vio.vm > vm_threshold || vio.pg > pg_threshold || vio.qg > qg_threshold || vio.sm > sm_threshold
        if vio.sm > sm_threshold
            branch_vio = branch_violations_sorted_ratec(network, sol_tmp)[1]

            if !haskey(branch_cuts_active, cont.label) || !(branch_vio.branch_id in branch_cuts_active[cont.label])
                info(_LOGGER, "adding flow cut on cont $(cont.label) branch $(branch_vio.branch_id) due to constraint flow violations $(branch_vio.sm_vio)")

                b_inv_cont = compute_B_inverse(network_ref, bus_idx2id, bus_id2idx, inactive_branches=Set([cont.idx]))
                branch_ptdf_cont = compute_branch_ptdf_single(network_ref, b_inv_cont, bus_id2idx, branch_vio.branch_id)

                cut = built_flow_cut(cont.label, branch_vio.branch_id, 1.0, branch_ptdf_cont, bus_idx2id)
                push!(branch_cuts, cut)
                break
            else
                warn(_LOGGER, "skipping active flow cut on cont $(cont.label) branch $(branch_vio.branch_id) with constraint flow violations $(branch_vio.sm_vio)")
            end
        end
    end

    if p_delta != 0.0
        warn(_LOGGER, "re-adjusting loads by $(-p_delta)")
        for (i,load) in load_setpoint
            load["pd"] = load["pd"] - p_delta
        end
    end

    time_contingencies = time() - time_contingencies_start
    info(_LOGGER, "cont eval time: $(time_contingencies)")

    return (gen_cuts=gen_cuts, branch_cuts=branch_cuts)
end

""
function check_contingencies_branch_power_remote(cont_range, output_dir, cut_limit=1, solution_file="solution1.txt")
    if length(network_global) <= 0 || length(contingency_order_global) <= 0
        error(_LOGGER, "check_contingencies_branch_flow_remote called before load_network_global")
    end

    sol = read_solution1(network_global, output_dir=output_dir, state_file=solution_file)
    PowerModels.update_data!(network_global, sol)

    active_cuts = read_active_flow_cuts(output_dir=output_dir)
    gen_flow_cuts = []
    branch_flow_cuts = []
    for cut in active_cuts
        if cut.cont_type == "gen"
            push!(gen_flow_cuts, cut)
        elseif cut.cont_type == "branch"
            push!(branch_flow_cuts, cut)
        else
            warn(_LOGGER, "unknown contingency type in cut $(cut)")
        end
    end

    network = copy(network_global)
    contingencies = contingency_order_global[cont_range]
    network["gen_contingencies"] = [c for c in contingencies if c.type == "gen"]
    network["branch_contingencies"] = [c for c in contingencies if c.type == "branch"]

    cuts = check_contingencies_branch_power(network, total_cut_limit=cut_limit, gen_flow_cuts=gen_flow_cuts, branch_flow_cuts=branch_flow_cuts)

    return cuts
end

"""
Checks a given operating point against the contingencies to look for branch
flow violations.  The DC Power Flow approximation is used for flow simulation.
If a violation is found, computes a PTDF cut based on bus injections.  Uses the
participation factor based generator response model from the ARPA-e GOC
Challenge 1 specification.
"""
function check_contingencies_branch_power(network;
        gen_flow_cut_limit=10, branch_flow_cut_limit=10, total_cut_limit=typemax(Int64),
        gen_eval_limit=typemax(Int64), branch_eval_limit=typemax(Int64), sm_threshold=0.01,
        gen_flow_cuts=[], branch_flow_cuts=[]
        )

    if InfrastructureModels.ismultinetwork(network)
        error(_LOGGER, "the branch flow cut generator can only be used on single networks")
    end
    time_contingencies_start = time()

    gen_cuts_active = Dict()
    for gen_cut in gen_flow_cuts
        if !haskey(gen_cuts_active, gen_cut.cont_label)
            gen_cuts_active[gen_cut.cont_label] = Set{Int}()
        end
        push!(gen_cuts_active[gen_cut.cont_label], gen_cut.branch_id)
    end

    branch_cuts_active = Dict()
    for branch_cut in branch_flow_cuts
        if !haskey(branch_cuts_active, branch_cut.cont_label)
            branch_cuts_active[branch_cut.cont_label] = Set{Int}()
        end
        push!(branch_cuts_active[branch_cut.cont_label], branch_cut.branch_id)
    end

    network_ref = PowerModels.build_ref(network)
    network_ref = network_ref[:nw][0]

    solution_base = extract_solution(network; branch_flow=true)

    pd_total = sum(load["pd"] for (i,load) in network_ref[:load])
    p_losses = sum(gen["pg"] for (i,gen) in network_ref[:gen]) - pd_total
    load_setpoint = network["load"]
    p_delta = 0.0

    if p_losses > pg_loss_tol
        load_count = length(load_setpoint)
        p_delta = p_losses/load_count
        for (i,load) in load_setpoint
            load["pd"] = load["pd"] + p_delta
        end
        warn(_LOGGER, "active power losses found $(p_losses) increasing loads by $(p_delta)")
    end

    bus_idx2id, bus_id2idx = build_indices(network_ref[:bus])
    b_base = compute_B_sparse(network_ref, bus_idx2id, bus_id2idx)
    b_inv_base = compute_B_inverse(network_ref, bus_idx2id, bus_id2idx)
    #branch_ptdf_base = compute_branch_ptdf(network_ref, b_inv_base, bus_id2idx)
    bus_inj_base = compute_bus_injections(network_ref, solution_base["gen"], load_setpoint, bus_idx2id)

    @assert length(network_ref[:ref_buses]) == 1
    ref_bus_id, ref_bus = first(network_ref[:ref_buses])
    ref_bus_idx = bus_id2idx[ref_bus["index"]]

    gen_cont_total = length(network["gen_contingencies"])
    branch_cont_total = length(network["branch_contingencies"])

    gen_eval_limit = min(gen_eval_limit, gen_cont_total)
    branch_eval_limit = min(branch_eval_limit, branch_cont_total)

    gen_cap = Dict(gen["index"] => sqrt(max(abs(gen["pmin"]), abs(gen["pmax"]))^2 + max(abs(gen["qmin"]), abs(gen["qmax"]))^2) for (i,gen) in network["gen"])
    network["gen_contingencies"] = sort(network["gen_contingencies"], rev=true, by=x -> gen_cap[x.idx])
    gen_contingencies = network["gen_contingencies"][1:gen_eval_limit]

    line_imp_mag = Dict(branch["index"] => branch["rate_a"]*sqrt(branch["br_r"]^2 + branch["br_x"]^2) for (i,branch) in network["branch"])
    network["branch_contingencies"] = sort(network["branch_contingencies"], rev=true, by=x -> line_imp_mag[x.idx])
    branch_contingencies = network["branch_contingencies"][1:branch_eval_limit]

    gen_cuts = []
    for (i,cont) in enumerate(gen_contingencies)
        if length(gen_cuts) >= gen_flow_cut_limit
            info(_LOGGER, "hit gen flow cut limit $(gen_flow_cut_limit)")
            break
        end
        if length(gen_cuts) >= total_cut_limit
            info(_LOGGER, "hit total cut limit $(total_cut_limit)")
            break
        end
        #info(_LOGGER, "working on ($(i)/$(gen_eval_limit)/$(gen_cont_total)): $(cont.label)")

        sol_tmp = deepcopy(solution_base)
        cont_gen = sol_tmp["gen"]["$(cont.idx)"]
        pg_lost = cont_gen["pg"]
        qg_lost = cont_gen["qg"]
        cont_gen["pg"] = 0.0
        cont_gen["qg"] = 0.0

        gen = network_ref[:gen][cont.idx]
        gen_bus = network_ref[:bus][gen["gen_bus"]]
        gen_set = network_ref[:area_gens][gen_bus["area"]]

        alpha_gens = [gen["alpha"] for (i,gen) in network_ref[:gen] if gen["index"] != cont.idx && gen["index"] in gen_set]
        if length(alpha_gens) == 0 || isapprox(sum(alpha_gens), 0.0)
            warn(_LOGGER, "no available active power response in cont $(cont.label), active gens $(length(alpha_gens))")
            continue
        end
        alpha_total = sum(alpha_gens)
        delta = pg_lost/alpha_total
        sol_tmp["delta"] = delta
        #info(_LOGGER, "$(pg_lost) - $(alpha_total) - $(delta)")

        for (i,gen) in network_ref[:gen]
            sol_gen = sol_tmp["gen"]["$(i)"]
            if gen["index"] != cont.idx && gen["index"] in gen_set
                sol_gen["pg"] = sol_gen["pg"] + gen["alpha"]*delta
            end
            sol_gen["qg"] = 0.0
        end

        bus_inj_cont = compute_bus_injections(network_ref, sol_tmp["gen"], load_setpoint, bus_idx2id)
        #p_flows = compute_branch_flows(branch_ptdf_base, bus_inj_cont)
        va_vector = []
        try
            va_vector = solve_theta(b_base, bus_inj_cont, ref_bus_idx)
        catch exception
            warn(_LOGGER, "linear solve failed on $(cont.label)")
            continue
        end
        p_flows = compute_flows(network_ref, va_vector, bus_id2idx)

        for (l,flow) in p_flows
            sol_branch = sol_tmp["branch"]["$(l)"]
            sol_branch["pf"] =  flow
            sol_branch["pt"] = -flow
            sol_branch["qf"] = 0.0
            sol_branch["qt"] = 0.0
        end

        #PowerModels.print_summary(sol_tmp)

        network["gen"]["$(cont.idx)"]["gen_status"] = 0
        vio = compute_violations_ratec(network, sol_tmp)
        network["gen"]["$(cont.idx)"]["gen_status"] = 1

        #info(_LOGGER, "$(cont.label) violations $(vio)")
        #if vio.vm > vm_threshold || vio.pg > pg_threshold || vio.qg > qg_threshold || vio.sm > sm_threshold
        if vio.sm > sm_threshold
            branch_vios = branch_violations_sorted_ratec(network, sol_tmp)
            branch_vio = branch_vios[1]

            if !haskey(gen_cuts_active, cont.label) || !(branch_vio.branch_id in gen_cuts_active[cont.label])
                info(_LOGGER, "adding flow cut on cont $(cont.label) branch $(branch_vio.branch_id) due to constraint flow violations $(branch_vio.sm_vio)")

                branch_ptdf_cont = compute_branch_ptdf_single(network_ref, b_inv_base, bus_id2idx, branch_vio.branch_id)
                cut = built_flow_cut(cont.label, branch_vio.branch_id, 1.0, branch_ptdf_cont, bus_idx2id)
                cut = (gen_id=cont.idx, cont_label=cut.cont_label, branch_id=cut.branch_id, rating_level=cut.rating_level, bus_injection=cut.bus_injection)
                push!(gen_cuts, cut)
            else
                warn(_LOGGER, "skipping active flow cut on cont $(cont.label) branch $(branch_vio.branch_id) with constraint flow violations $(branch_vio.sm_vio)")
            end

        end
    end

    # work around for julia compiler bug
    #num_buses = length(network_ref[:bus])
    #b_cont = zeros(Float64, num_buses, num_buses)

    branch_cuts = []
    for (i,cont) in enumerate(branch_contingencies)
        if length(branch_cuts) >= branch_flow_cut_limit
            info(_LOGGER, "hit branch flow cut limit $(branch_flow_cut_limit)")
            break
        end
        if length(gen_cuts) + length(branch_cuts) >= total_cut_limit
            info(_LOGGER, "hit total cut limit $(total_cut_limit)")
            break
        end

        #info(_LOGGER, "working on ($(i)/$(branch_eval_limit)/$(branch_cont_total)): $(cont.label)")
        sol_tmp = deepcopy(solution_base)
        if haskey(sol_tmp, "branch")
            cont_branch = sol_tmp["branch"]["$(cont.idx)"]
            cont_branch["pf"] = 0.0
            cont_branch["pt"] = 0.0
            cont_branch["qf"] = 0.0
            cont_branch["qt"] = 0.0
        end

        #delta = 0.0
        for (i,gen) in network_ref[:gen]
            sol_gen = sol_tmp["gen"]["$(i)"]
            sol_gen["pg"] = sol_gen["pg"] #+ gen["alpha"]*delta
        end

        b_cont = compute_B_sparse(network_ref, bus_idx2id, bus_id2idx, inactive_branches=Set([cont.idx]))
        # work around for julia compiler bug
        #update_B!(b_cont, network_ref, bus_idx2id, bus_id2idx, inactive_branches=Set([cont.idx]))

        va_vector = []
        try
            va_vector = solve_theta(b_cont, bus_inj_base, ref_bus_idx)
        catch exception
            warn(_LOGGER, "linear solve failed on $(cont.label)")
            continue
        end
        p_flows = compute_flows(network_ref, va_vector, bus_id2idx)

        for (l,flow) in p_flows
            sol_branch = sol_tmp["branch"]["$(l)"]
            if l != cont.idx
                sol_branch["pf"] =  flow
                sol_branch["pt"] = -flow
            else
                sol_branch["pf"] = 0.0
                sol_branch["pt"] = 0.0
            end

            sol_branch["qf"] = 0.0
            sol_branch["qt"] = 0.0
        end

        #PowerModels.print_summary(sol_tmp)

        network["branch"]["$(cont.idx)"]["br_status"] = 0
        vio = compute_violations_ratec(network, sol_tmp)
        network["branch"]["$(cont.idx)"]["br_status"] = 1

        #info(_LOGGER, "$(cont.label) violations $(vio)")
        #if vio.vm > vm_threshold || vio.pg > pg_threshold || vio.qg > qg_threshold || vio.sm > sm_threshold
        if vio.sm > sm_threshold
            branch_vio = branch_violations_sorted_ratec(network, sol_tmp)[1]
            if !haskey(branch_cuts_active, cont.label) || !(branch_vio.branch_id in branch_cuts_active[cont.label])
                info(_LOGGER, "adding flow cut on cont $(cont.label) branch $(branch_vio.branch_id) due to constraint flow violations $(branch_vio.sm_vio)")

                b_inv_cont = compute_B_inverse(network_ref, bus_idx2id, bus_id2idx, inactive_branches=Set([cont.idx]))
                branch_ptdf_cont = compute_branch_ptdf_single(network_ref, b_inv_cont, bus_id2idx, branch_vio.branch_id)

                cut = built_flow_cut(cont.label, branch_vio.branch_id, 1.0, branch_ptdf_cont, bus_idx2id)
                push!(branch_cuts, cut)
            else
                warn(_LOGGER, "skipping active flow cut on cont $(cont.label) branch $(branch_vio.branch_id) with constraint flow violations $(branch_vio.sm_vio)")
            end
        end
    end

    if p_delta != 0.0
        warn(_LOGGER, "re-adjusting loads by $(-p_delta)")
        for (i,load) in load_setpoint
            load["pd"] = load["pd"] - p_delta
        end
    end

    time_contingencies = time() - time_contingencies_start
    info(_LOGGER, "cont eval time: $(time_contingencies)")

    return (gen_cuts=gen_cuts, branch_cuts=branch_cuts)
end


text:
using Test
using YaoCompiler
using YaoLocations
using LinearAlgebra
using YaoArrayRegister
using YaoCompiler.Intrinsics

struct JLEmulationTarget <: NativeJuliaTarget end

for G in [:X, :Y, :Z, :H, :S, :T]
    @eval function Intrinsics.apply(r::ArrayReg, ::Intrinsics.$(Symbol(G, :Gate)), locs::Locations)
        instruct!(r, Val($(QuoteNode(G))), Tuple(locs))
        return
    end

    @eval function Intrinsics.apply(
        r::ArrayReg,
        ::Intrinsics.$(Symbol(G, :Gate)),
        locs::Locations,
        ctrl::CtrlLocations,
    )
        ctrl_locs = Tuple(ctrl.storage)
        ctrl_configs = ntuple(length(ctrl)) do k
            Int(ctrl.flags[k])
        end
        instruct!(r, Val($(QuoteNode(G))), Tuple(locs), ctrl_locs, ctrl_configs)
        return
    end
end

for axis in [:x, :y, :z]
    @eval function Intrinsics.apply(r::ArrayReg, op::Intrinsics.$(Symbol(:R, axis)), locs::Locations)
        instruct!(r, Val($(QuoteNode(Symbol(:R, axis)))), Tuple(locs), op.θ)
        return
    end

    @eval function Intrinsics.apply(
        r::ArrayReg,
        op::Intrinsics.$(Symbol(:R, axis)),
        locs::Locations,
        ctrl::CtrlLocations,
    )
        ctrl_locs = Tuple(ctrl.storage)
        ctrl_configs = ntuple(length(ctrl)) do k
            Int(ctrl.flags[k])
        end
        instruct!(r, Val($(QuoteNode(Symbol(:R, axis)))), Tuple(locs), ctrl_locs, ctrl_configs, op.θ)
        return
    end
end

function Intrinsics.apply(r::ArrayReg, op::Intrinsics.shift, locs::Locations)
    instruct!(r, Diagonal(ComplexF64[1.0, exp(im * op.θ)]), Tuple(locs))
end

@noinline function Intrinsics.measure(r::ArrayReg, locs::Locations)
    result = YaoArrayRegister.measure!(r, Tuple(locs))
    return MeasureResult{typeof(result)}(result)
end

@device function test_intrinsic(theta, phi)
    # syntax sugar
    1 => X
    @gate 2 => Z
    @ctrl 1 4 => Rx(theta)
    @ctrl 2 4 => Ry(phi)
    return
end

@device function test_location_map(theta)
    3:6 => test_intrinsic(theta, 2.5)
    @ctrl 1 3:6 => test_intrinsic(theta, 2.1)
    return
end

r = rand_state(10)
op = test_location_map(1.0)
f = YaoCompiler.compile(JLEmulationTarget(), Intrinsics.apply, Tuple{typeof(r),typeof(op)})

function location_map(r, theta)
    instruct!(r, Val(:X), (3,))
    instruct!(r, Val(:Z), (4,))
    instruct!(r, Val(:Rx), (6,), (3,), (1,), theta)
    instruct!(r, Val(:Ry), (6,), (4,), (1,), 2.5)

    instruct!(r, Val(:X), (3,), (1,), (1,))
    instruct!(r, Val(:Z), (4,), (1,), (1,))
    instruct!(r, Val(:Rx), (6,), (1, 3), (1, 1), theta)
    instruct!(r, Val(:Ry), (6,), (1, 4), (1, 1), 2.1)
    return
end

r1, r2 = copy(r), copy(r)
f(r1, op)
location_map(r2, 1.0)
@test r1 ≈ r2


text:
###
# Configuration
###
include("drawcircle.jl")

function loadConfig(dict)
    for key in keys(dict)
        sym = symbol(key)
        @eval $sym = $(dict[key])
    end
end

function loadConfig(dict, blacklist)
    for key in keys(dict)
        if !(key in blacklist)
            sym = symbol(key)
            @eval $sym = $(dict[key])
        end
    end
end

function saveConfig(keys)
    dict = Dict{String, Any}()
    for key in keys
        sym = symbol(key)
        dict[key] = eval(sym)
    end
    return dict
end

function checkVars(vars, keys)
    for v in vars
        if !(v in keys)
            return false
        end
    end
    return true
end

baseConfig = {
    "fieldSize"             => 2400,
    "fieldRes"              => 40,
    "timeTotal"             => 15000,
    "stepIntegration"       => 0.5,
    "visInterval"           => 5,
    "storeStep"             => 100,
    "epsilon"               => 0.01,
    "deathEpsilon"          => 0.01,
    "stableTime"            => 5000,
    "v_mR"                  => 180, # membraneRadius
    "v_mT"                  => 170,  # membraneThickness
    "v_fD"                   => 20, # food refill distance
    "avgM"                  => 0.8, #initial concentration of M in membrane
    "avgA"                  => 0.6,
    "avgF"                  => 0.18,
    "avgT"                  => 0.0,
    "diffusionF"            => 20, # absolute diffusion rate
    "diffusionA"            => 20, # absolute diffusion rate
    "diffusionM"            => 20, # absolute diffusion rate
    "diffusionW"            => 20, # absolute diffusion rate
    "attractionRate"        => 0.1, # attraction of directionality

    "MW_repulsion"          => 7.0,
    "MA_repulsion"          => 7.0,
    "MM_repulsion"          => 0.5,
    "AA_repulsion"          => 0.5,
    "long_direction"        => 20.0,

    "decayA"                => 0.001, # this is per timestep, not integration step.
    "decayM"                => 0.001,

    # rules for food
    "flowRateF"             => 0.8,

    # reaction 1 : a11*A + f11*F + m11*M + w11*W <-> a12*A + f12*F + m12*M + w12*W
    "a11" => 2, "f11" => 1, "m11" => 0, "w11" => 0, # negative numbers inhibit
    "a12" => 3, "f12" => 0, "m12" => 0, "w12" => 0, # inhibition needs to balanced for there not to be creation

    "kf1" => 0.1, # forward reaction rate
    "kb1" => 0.0, # backwards reaction rate

    # reaction 2 : a21*A + f21*F + m21*M + w21*W <-> a22*A + f22*F + m22*M + w22*W
    "a21" => 2, "f21" => 1, "m21" => 0, "w21" => 0,
    "a22" => 2, "f22" => 0, "m22" => 1, "w22" => 0,

    "kf2" => 0.1, # forward reaction rate
    "kb2" => 0.0, # backwards reaction rate

    #reaction 3 : a31*A + f31*F + m31*M + w31*W <-> a32*A + f32*F + m32*M + w32*W
    "a31" => 0, "f31" => 0, "m31" => 0, "w31" => 0,
    "a32" => 0, "f32" => 0, "m32" => 0, "w32" => 0,

    "kf3" => 0.0, # forward reaction rate
    "kb3" => 0.0 # backwards reaction rate
}

historyVars = ["history_W", "history_A", "history_M", "history_M_pot", "history_F", "history_dir"]
dataVars = ["history_W", "history_A", "history_M", "history_M_pot", "history_F", "history_dir", "Avec", "Fvec", "Mvec", "Wvec", "DAvec", "DMvec", "tx"]

function updateDependentValues()
    @eval begin
        saturationF = avgF
        fieldSizeY = fieldSize
        fieldSizeX = fieldSizeY
        fieldResY = fieldRes
        fieldResX = fieldResY
        _xc = fieldSize/2 # circle center's x position
        _yc = fieldSize/2 # circle center's y position

        # adjust membrane for given resolution
        xc = ceil(_xc*fieldRes/fieldSize)
        yc = ceil(_yc*fieldRes/fieldSize)
        mR = ceil(v_mR*fieldRes/fieldSize) # membraneRadius
        mT = ceil(v_mT*fieldRes/fieldSize) # membraneThickness
        fD = ceil(v_fD*fieldRes/fieldSize) # membraneThickness

        Frefill = zeros(fieldRes, fieldRes)
        for r in 0:0.01:(mR+mT+fD)
            drawcircle!(Frefill, xc, yc, r)
        end
        FrefillBinMask = Frefill .< 0.5
    end
end


text:
@testset "distance" begin
    x = [3, 4]
    y = [0.0, 0.0]

    @test dist(x, y) ≈ 5

    X = [x y]'
    Y = [y x x]'

    D = dist(X, Y)
    @test size(D) == (2, 3)
    @test D[1, 1] ≈ 5
    @test D[1, 3] ≈ 0

    @test size(dist(Y, Y)) == (3, 3)

    @test hamming("banana", "banaan") == 2
end

@testset "graph" begin

    @test [(1, 2), (3, 3)] isa EdgeList
    @test !([(1, 2), (3, 3, 4)] isa EdgeList)
    edges = [(0.1, 3, 2), (2.0, 1, 2)]
    @test edges isa WeightedEdgeList

    adjlist = edges2adjlist(edges)
    @test adjlist isa AdjList{Float64, Int}
    @test (2.0, 1) in adjlist[2]
    @test (2.0, 1) ∉ edges2adjlist(edges, double=false)[2]

    edgesrec = adjlist2edges(adjlist)
    @test (0.1, 3, 2) in edgesrec && (0.1, 2, 3) in edgesrec

    @test nvertices(adjlist) == nvertices(edges) == 3
    @test vertices(adjlist) == vertices(edges) == Set([1, 2, 3])
    @test isconnected(adjlist) == isconnected(edges) == true
    @test !isconnected([(0.1, 1, 2), (0.1, 2, 3), (0.01, 5, 4)])

end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()
##
using Plots
using MCPhylo
using CSV
using DataFrames
using LinearAlgebra

##

function getEqulibrium(rates)
      q = zeros(4, 4)
      indices = [
            (2, 1),
            (3, 1),
            (4, 1),
            (1, 2),
            (3, 2),
            (4, 2),
            (1, 3),
            (2, 3),
            (4, 3),
            (1, 4),
            (2, 4),
            (3, 4),
      ]
      q[[CartesianIndex(x...) for x in indices]] = rates
      dia = sum(q, dims = 2)
      q[diagind(4, 4)] = -dia
      eq = svd(q).U[:, end]
      eq /= sum(eq)
end

##

fpairs = CSV.read("../data/fpairs.txt", DataFrame, header = false)[:, 1]

nrc = findall(occursin.("NRc", fpairs))
##

colnames = ["r_$(i)_$(j)" for j in 1:4 for i in 1:4 if i !=j]

##

df_ = []
for charNum = 1:28
      charName = fpairs[charNum]
      sim = read(
            "modelFitting/output/universal_$(lpad(charNum, 2, "0")).jls",
            ModelChains,
      )
      rates = exp.(vcat([sim[:, :s_rates, :].value[:, :, k] for k in sim.chains]...))
      if charNum ∈ nrc
            samples = samples[:, [4, 6, 5, 1, 3, 2, 11, 10, 12, 8, 7, 9]]
      end

      eqs = mapslices(getEqulibrium, rates, dims=2)

      rateMedians = mapslices(median, rates, dims=1)[1,:]
      rateHPD = permutedims(mapslices(hpd, rates, dims=1))

      rDF = DataFrame(
            :featurePair => charName,
            :variable => colnames,
            :median => rateMedians,
            :HPD_l => rateHPD[:,1],
            :HPD_u => rateHPD[:,2],
      )

      eqMedians = mapslices(median, eqs, dims=1)[1,:]
      eqHPD = permutedims(mapslices(hpd, eqs, dims=1))

      eqDF = DataFrame(
            :featurePair => charName,
            :variable => ["e1", "e2", "e3", "e4"],
            :median => eqMedians,
            :HPD_l => eqHPD[:,1],
            :HPD_u => eqHPD[:,2],
      )
      push!(df_, vcat(rDF, eqDF))
end

df = vcat(df_...)

CSV.write("../data/posteriorStatistics.csv", df)

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
module ImageViewer

using GLMakie
using ColorSchemes

export imview

function imview(imgin; color_map_menu=true)
    # create figure to draw in
    fig = Figure()
    # that will be the image
    ax = Axis(fig[1,2], alignmode=Inside(), xrectzoom=false, yrectzoom=false)
    ax.aspect = DataAspect()

    # for small widgets
    gl = fig[1,1, Left()] = GridLayout()

    # fix image layout
    colsize!(fig.layout, 1, Fixed(100))
    #rowsize!(fig.layout, 0, Fixed(900))

    # specialize for three dimensions
    if ndims(imgin) == 3
        # add some sliders
        img = add_depth_slier(fig, (1,), imgin)
        img = add_γ_slider(gl[1,1], (1,), img)
        img = add_γ_slider(gl[2,1], (1,), img)

        # plot
        h = heatmap!(ax, img, interpolate=false, tellheight=true)

        # add color map menu
        if color_map_menu
            cmap_menu = add_colormap_menu(gl[3,1], (1,))
            cmap = on(cmap_menu.selection) do s
                h.colormap = s
            end
        end
        # add color bar
        cb = Colorbar(fig[1, 3], h, tellheight=true, width = 20, ticklabelspace=30f0)
        fix_sizes_to_img(fig, ax)
    end

    return fig
end

function fix_sizes_to_img(fig, ax)
    rowsize!(fig.layout, 1, ax.scene.px_area[].widths[2])
    colsize!(fig.layout, 2, ax.scene.px_area[].widths[1])
end

function add_γ_slider(gr, pos, img)
    γ_r = γ_range()
    Label(gr[pos..., 1, Top()], "γ:", textsize=20, halign=:left, valign=:top)
    γ = Slider(gr[pos..., 1], range=γ_r, startvalue=1, width=100)
    γ_text = @lift to_str($(γ.value))
    Label(gr[pos..., 1, Top()], γ_text)
    img = @lift($img .^ eltype($img)($(γ.value)))
    return img
end

function add_depth_slier(gr, pos, img)
    depth = Slider(gr[2,2, Top()], range=axes(img, 3))
    rowsize!(gr.layout, 2, Fixed(10))

    imgout = lift(depth.value) do d
            return view(img, :, :, d)
    end

    return imgout
end

function add_colormap_menu(gr, pos;
            cmap_options=["curl", "thermal", "roma", "gist_rainbow", "hsv", "twilight"]
        )

    Label(gr[pos..., 1, Top()], "Color Map")
    cmap = Menu(gr[pos..., 1], options=cmap_options)
    return cmap
end

function apply_gamma(img, γ)
    return img .^ γ
end

function γ_range()
    γ_range = 10.0 .^ (-2:0.05:1)
end

function to_str(x::AbstractFloat; digits=2)
    string(round(x, digits=digits))
end

# precompilation try
let
    while true
        imview(randn(Float32, (5,5,5)))
        imview(randn(Float64, (5,5,5)))
        break
    end
end

precompile(imview, (Array{Float32, 3}, Array{Float64, 3}))

end # module


text:
using Colors, Distributions, Bijectors, Zygote, StatsPlots, MCMCChains
using StatsFuns: logsumexp, logit, logistic
using Revise, AdvancedHMC

########
# Data #
########
const DATA = [RGB(0.5, 0.0, 0.0), RGB(0.6, 0.0, 0.0), RGB(0.4, 0.0, 0.0)]
const DATA_MIXTURE = vcat(DATA, [RGB(0.5, 0.0, 1.0), RGB(0.6, 0.0, 1.0), RGB(0.4, 0.0, 1.0)])

#######################
# Color distributions #
#######################

struct UniformRGB end
Distributions.rand(::UniformRGB) = RGB{Float64}(rand(3)...)
Distributions.logpdf(::UniformRGB, ::RGB{Float64}) = zero(Float64)

# struct BetaRGB
#     α
#     β
# end
# Distributions.rand(b::BetaRGB) =
#     RGB(map(rand, [Beta(b.α[i], b.β[i]) for i = 1:3])...)
# function Distributions.logpdf(b::BetaRGB, x::RGB{T}) where {T<:AbstractFloat}
#     betas = [Beta(b.α[i], b.β[i]) for i = 1:3]
#     xs = [x.r, x.g, x.b]
#     return sum(map(i -> logpdf(betas[i], xs[i]), 1:3))
# end

function Bijectors.logpdf_with_trans(d, c::RGB{Float64})
    θ = [c.r, c.g, c.b]
    return logpdf(d, c) + sum(log.(θ .* (one(Float64) .- θ)))
end

struct NormalRGB
    μ::RGB{Float64}
    σ::Float64
end
Distributions.logpdf(d::NormalRGB, x::RGB{Float64}) = logpdf(Normal(zero(Float64), d.σ), colordiff(d.μ, x)^2)

###################
# Sampling helper #
###################

function sample_colors(θ_init, logjoint, grad; stepsize=1e-3, n_samples=2_000, n_adapts=1_000)
    # Dimension
    dim = length(θ_init)

    # Define metric space, Hamiltonian and sampling method
    metric = DiagEuclideanMetric(dim)
    h = Hamiltonian(metric, logjoint, grad)
    # prop = NUTS(Leapfrog(find_good_eps(h, θ_init)), 5, 1000.0)
    prop = StaticTrajectory(Leapfrog(stepsize), 20)
    # adaptor = StanNUTSAdaptor(n_adapts, Preconditioner(metric), NesterovDualAveraging(0.8, prop.integrator.ϵ))

    # Sampling
    # samples = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; progress=true)
    samples = sample(h, prop, θ_init, n_samples; progress=true)
    return samples
end

############################
# Model 1: single Gaussian #
############################
# Log-joint and it's gradient
function logπ(θ::Vector{Float64}) :: Float64
    c = RGB{Float64}(θ...)
    logprior = logpdf(UniformRGB(), c)
    loglikelihood = sum([logpdf(NormalRGB(t, 5one(Float64)), c) for t in  DATA])
    return logprior + loglikelihood
end
# @code_warntype logπ(rand(3))
∂logπ∂θ(θ::Vector{Float64}) :: Vector{Float64} = vec(collect(gradient(logπ, θ)[1]))
# @code_warntype ∂logπ∂θ(rand(6))
samples = sample_colors(rand(3), logπ, ∂logπ∂θ)
plot(
    Chains(
        reshape(hcat(samples...)', (2_000, 3, 1)),
        ["c.r", "c.g", "c.b"]
    )
)

colors = map(s -> RGB(s...), samples)
mean_samples = mean(samples)
mean_color = RGB(mean_samples...)

###################################
# Model 1: mixture of 2 Gaussians #
###################################
# Log-joint and it's gradient
logπ_mixture(θlink::Vector{Float64}) = logπ_mixture(θlink[1:3], θlink[4:6])
# @Mike: I cannot define `logπ_mixture(θ::Vector{Float64})` directly
#        and construct colors by `c1 = RGB(θ[1:3]...)`, which gives me
#        error when calling `gradient`.
function logπ_mixture(θ1link::Vector{Float64}, θ2link::Vector{Float64}) :: Float64
    p = ones(Float64) / 2
    σ = 20one(Float64)
    c1 = RGB{Float64}(logistic.(θ1link)...)
    c2 = RGB{Float64}(logistic.(θ2link)...)
    logprior = logpdf_with_trans(UniformRGB(), c1) + logpdf_with_trans(UniformRGB(), c2)
    loglikelihood = sum(
        [logsumexp([logpdf(NormalRGB(t, σ), c1) - log(2one(Float64)),
                    logpdf(NormalRGB(t, σ), c2) - log(2one(Float64))]) for t in DATA_MIXTURE]
    )
    return logprior + loglikelihood
end
# @code_warntype logπ(rand(6))
∂logπ∂θ_mixture(θlink::Vector{Float64}) :: Vector{Float64} = ∂logπ∂θ_mixture(θlink[1:3], θlink[4:6])
∂logπ∂θ_mixture(θ1link::Vector{Float64}, θ2link::Vector{Float64}) :: Vector{Float64} =
    vcat(map(collect, gradient(logπ_mixture, θ1link, θ2link))...)
# @code_warntype ∂logπ∂θ_mixture(rand(6))
sampleslink_mixture = sample_colors(randn(6), logπ_mixture, ∂logπ∂θ_mixture; stepsize=1e-2)
samples_mixture = map(s -> logistic.(s), sampleslink_mixture)
plot(
    Chains(
        reshape(hcat(samples_mixture...)', (2_000, 6, 1)),
        ["c1.r", "c1.g", "c1.b", "c2.r", "c2.g", "c2.b"]
    )
)

colors_mixture = map(s -> [RGB(s[1:3]...), RGB(s[4:6]...)], samples_mixture)
map(cs->cs[1], colors_mixture)
map(cs->cs[2], colors_mixture)
mean_samples_mixture = mean(samples_mixture)
mean_color = [RGB(mean_samples_mixture[1:3]...), RGB(mean_samples_mixture[4:6]...)]

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper