text:
module AVFormat
  include(joinpath(dirname(@__FILE__),"..","..","..","init.jl"))
  w(f) = joinpath(avformat_dir, f)

  using ..AVUtil
  using ..AVCodecs

  include(w("LIBAVFORMAT.jl"))
end


text:
using MolecularBoxes

using Test
using StaticArrays

#a 3D vector for testing
const Vec = SVector{3,Float64}

const BoxV =  Box{Vec,3,(true,true,true)}

@testset "Test Box constructor" begin
    vectors  = ( Vec(3,0,0),
                 Vec(0,4,0),
                 Vec(0,0,5))
    avectors = ([3.,0.,0.],
                [0.,4.,0.],
                [0.,0.,5.])
    ppp = (true,true,true)
    pp = (true,true)
    @test Box{Vec,3,ppp}(vectors) != nothing
    box = BoxV(vectors)
    @test_throws MethodError Box{Vec,2,ppp}(vectors)
    @test_throws ErrorException Box{Vec,3,pp}(vectors)
    @test_throws ErrorException Box{Vec,3,ppp}((
        Vec(1,1,0),
        Vec(0,2,0),
        Vec(0,0,3),
    ))

    @test BoxV(vectors) == Box(Vec(3,4,5))

end

@testset "Test getters" begin
    p = (true,false,true)
    vectors  = (
        Vec(3,0,0),
        Vec(0,4,0),
        Vec(0,0,5),
    )
    box = Box{Vec,3,p}(vectors)
    @test isperiodic(box) == p
    @test box.vectors == vectors
    @test box.lengths == convert(Vec, collect(vectors[i][i] for i in 1:3) )
end

lengths = Vec(3,4,5)
hi = lengths

box = Box(lengths, periodic=(true,true,true))
boxpfp = Box(lengths, periodic=(true,false,true))
boxppf = Box(lengths, periodic=(true,true,false))

v = Vec(-1,3,8)
v1 = Vec(0.5,1,1)
v2 = lengths+Vec(-0.5,-1,-1)

@testset "test dimensions/coordinates access functions" begin
    @test box.lengths == lengths
end

@testset "Test separation" begin
    @test separation(v1,v2,box) == Vec(1,-2,2)
    @test separation(v1,v2,boxppf) == Vec(1,-2,-3)
    @test separation(v2,v1,boxppf) == Vec(-1,-2, 3)
    #@pending separation(need,more,tests) --> Vec(x,x,x)
    # especially for boundary cases
end

@testset "Test centre-of-mass calculation" begin
    @test center_of_mass([v1,v2], box) == Vec(3,2,5)
    @test center_of_mass([v1,v2], box, weights = [2.0, 2.0]) == Vec(3,2,5)
    # could user a better test of weights
    @test_throws ErrorException center_of_mass([v1,v2], boxpfp)
end

ext: jl


text:
module FunctionalStateMachine

using Dates
using Requires
using DocStringExtensions
using ProgressMeter

export
  StateMachine,
  emptyState,
  exitStateMachine,
  sandboxStateMachineStep,
  getStateLabel,
  histStateMachineTransitions

include("StateMachine.jl")

function __init__()
  @require Graphs="86223c79-3864-5bf0-83f7-82e725a168b6" begin
    @info "Adding graph visualization tools."

    export
      histGraphStateMachineTransitions,
      drawStateTransitionStep,
      drawStateMachineHistory,
      animateStateMachineHistoryByTime,
      animateStateMachineHistoryByTimeCompound

    include("StateMachineAnimation.jl")
  end
end

end


text:
# Contains script to run SIRSDE model using command line.

using LinearAlgebra
using StochasticDiffEq
using Base: @kwdef

using Random
using NPZ
using ArgParse

# Notation
# s: Susceptible
# i: Infected
# r: Removed/recovered/dead
# β: Infection rate
# γ: recovery rate
# R₀: infection_rate/recovery rate  (i.e. β/γ)
# σ: Controls volitility of R₀
# η: Controls mean reversion strength of R₀

"""
SIR model described using stochastic differential equations.
Ref: Inspired by https://julia.quantecon.org/continuous_time/covid_sde.html
"""
@kwdef struct SIRSDETask
    "Volatitility of R₀."
    σ::Float64 = 0.05
    "Mean reversion strength of R₀."
    η::Float64 = 0.05
    "Number of time steps."
    T::Int = 365
    "Number of individuals"
    N::Int = 100000
    "Initial proportions of [s,i,r]."
    sir_0::Vector{Float64} = [0.999, 0.001, 0]
    "Deterministic portion of SDE."
    f::Function = _sir_f!
    "Stochastic portion of SDE."
    g::Function = _sir_g!
    "Parameter names."
    θ_names::Vector{String} = ["β", "γ"]
    name::String = "SIRSDE"
end

"""
Deterministic portion of SDE (s, i, r are described using proportions).
"""
function _sir_f!(du, u, p, t)
    s, i, r, R₀ = u
    (; γ, R̄₀, η, σ) = p
    β = γ * R₀
    du[1] = -β * s * i        # ds/dt
    du[2] = β * s * i - γ * i   # di/dt
    du[3] = γ * i           # dr/dt
    du[4] = η * (R̄₀(t, p) - R₀) # dR₀/dt
    return nothing
end

"""
Stochastic portion of SDE (R₀ volatitility)
"""
function _sir_g!(du, u, p, t)
    s, i, r, R₀ = u
    (; γ, R̄₀, η, σ) = p
    du[1:3] .= 0
    du[4] = σ * sqrt(abs(R₀))
    return nothing
end

"""
Batched simulations using SIR model. θ has columns [β γ].
"""
function simulate(
    task::SIRSDETask,
    θ::AbstractMatrix{Float64})
    @assert all(θ .>= 0)
    (; σ, η, T, sir_0, f, g, N) = task

    x = Array{Float64}(undef, size(θ, 1), T)  # To store infection data
    for (θᵢ, xᵢ) in zip(eachrow(θ), eachrow(x))
        β, γ = θᵢ
        p = (
            R̄₀_ref=β / γ,  # Fixed reference (R̄₀ must be a function)
            R̄₀=(t, p) -> p.R̄₀_ref,
            γ=γ,
            σ=σ,
            η=η,
        )
        u_0 = [sir_0; β / γ]
        prob = SDEProblem(f, g, u_0, (0, T), p)
        sol = solve(prob, SOSRI())
        xᵢ .= [sol(t)[2] for t in 0:(T-1)]
    end
    return x .* N
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--seed"
        help = "Seed for simulations"
        arg_type = Int
        required = true

        "--theta_path"
        help = "Path to prior samples (npz file) to use for simulations"
        arg_type = String
        required = true

        "--output_path"
        help = "Path to output raw simulations"
        arg_type = String
        required = true
    end
    return parse_args(s)
end

function main()
    args = parse_commandline()
    Random.seed!(args["seed"])
    θ = npzread(args["theta_path"])["theta"]
    θ = convert(Matrix{Float64}, θ)
    task = SIRSDETask()
    x = simulate(task, θ)
    npzwrite(args["output_path"], x=x)
end

main()


text:
"""
	adjMat(network::AbstractNetwork, sparse::Bool=false)

Return the adjacency matrix of the network.
If `sparse` is true, return a sparse version of the matrix.

See also: `adjVet`
"""
function adjMat(network::AbstractNetwork; sparse::Bool=false)
	isnothing(network._adjMat) && calcAdjMat!(network)
	return sparse ? copy(network._adjMat) : BitArray(network._adjMat)
end

"""
	adjVet(network::AbstractNetwork)

Return the adjacency vector of the network

See also: `adjMat`
"""
function adjVet(network::AbstractNetwork)
	isnothing(network._adjMat) && calcAdjMat!(network)
	return adjMatToVet(network._adjMat)
end

"""
	numnodes(network::AbstractNetwork)

Return the number of nodes on the network.
"""
numnodes(network::AbstractNetwork) = network.N

"""
	numconnections(network::AbstractNetwork)

Return the number of connections present on the network.
"""
function numconnections(network::AbstractNetwork)
	isnothing(network._props.numConnections) && calcNumConnections!(network)
	return network._props.numConnections
end

"""
	maxconnections(network::AbstractNetwork)

Return the maximum amount of connections the network can support.
"""
maxconnections(network::AbstractNetwork) =
	Int(network.N * (network.N - 1) / (isdirected(network) ? 1 : 2))

"""
	isdirected(network::AbstractNetwork)

Return whether or not the network is directed.
"""
isdirected(network::AbstractNetwork) = network._props.directed

"""
	connectivity(network::AbstractNetwork, idx_node::Integer; dir_behaviour::Symbol=:total)

Return the connectivity of the specified node.
"""
function connectivity(
		network::AbstractNetwork,
		idx_node::Integer;
		dir_behaviour::Union{Symbol, String}=:total
	)
	hasnodeOrError(network, idx_node)

	dir_behaviour = Symbol(dir_behaviour)
	if isdirected(network) && ! (dir_behaviour in [:total, :in, :out, :mean, :both, :bi])
		throw(ArgumentError(
			"dir_behaviour parameter must be one of the follwing: :total, :in, :out, :mean, :both, :bi"
		))
	end

	return calcConnectivity(network, idx_node; dir_behaviour=dir_behaviour)
end

"""
	connectivities(network::AbstractNetwork; dir_behaviour::Symbol=:total)

Return the connectivities of the all node of the network.

If the network is directed, the `dir_behaviour` parameter controls which connectivity is returned,
`:in` for in-degree, `out` for out-degree, `:total` for the sum of in- and out-degree,
`:mean` for the mean of in- and out-degree, `:both` for both in- and out-degree, and
`:bi` for the number of neighbors for which both in and out connections are present.

See also: `connectivity`
"""
connectivities(network::AbstractNetwork; dir_behaviour::Symbol=:total) =
	[connectivity(network, i; dir_behaviour=dir_behaviour) for i in 1:network.N]

"""
	connectivity(network::AbstractNetwork; dir_behaviour::Symbol=:total)

Return the mean connectivity of the network.

If the network is directed, the `degree` parameter controls which connectivity is returned,
`:in` for in-degree, `out` for out-degree, `:total` for the sum of in- and out-degree,
`:mean` for the mean of in- and out-degree, `:both` for both in- and out-degree, and
`:bi` for the number of neighbors for which both in and out connections are present.

See also: `connectivities`
"""
connectivity(network::AbstractNetwork; dir_behaviour::Symbol=:total) =
	mean(connectivities(network; dir_behaviour=dir_behaviour))

"""
	clusteringcoefficient(network::AbstractNetwork, idx_node::Integer)

Return the clustering coefficient of the specified node.
"""
function clusteringcoefficient(network::AbstractNetwork, idx_node::Integer)
	hasnodeOrError(network, idx_node)

	return isnothing(network._props.clusteringCoefficients) ?
		calcClusteringCoefficient(network, idx_node) :
		network._props.clusteringCoefficients[idx_node]
end

"""
	clusteringcoefficient(network::AbstractNetwork)

Return the average clustering coefficient of the network.

See also: `clusteringcoefficients`
"""
function clusteringcoefficient(network::AbstractNetwork)
	isnothing(network._props.clusteringCoefficients) && calcClusteringCoefficients!(network)
	return mean(network._props.clusteringCoefficients)
end

"""
	clusteringcoefficients(network::AbstractNetwork)

Return the clustering coefficients of all nodes of the network.

See also: `clusteringcoefficient`
"""
function clusteringcoefficients(network::AbstractNetwork)
	isnothing(network._props.clusteringCoefficients) && calcClusteringCoefficients!(network)
	return network._props.clusteringCoefficients
end

"""
	transitivity(network::AbstractNetwork)

Return the transitivity of the network.
"""
function transitivity(network::AbstractNetwork)
	if isdirected(network)
		println("Transitivity not supported for directed networks")
		return
	end

	isnothing(network._props.transitivity) && calcTransitivity!(network)
	return network._props.transitivity
end

"""
	shortestpath(network::AbstractNetwork, source::Integer,	target::Integer=nothing)

Return the length of the shortest paths from the source node the every other node.
If `target` is specified, only the path from the source to the target is returned.
"""
function shortestpath(
	network::AbstractNetwork,
	source::Integer,
	target::Union{Integer, Nothing}=nothing
)
	hasnodeOrError(network, source)
	isnothing(target) || hasnodeOrError(network, target)

	isnothing(network._props.shortestPaths[source]) && calcShortestPath!(network, source)
	return isnothing(target) ?
		network._props.shortestPaths[source] :
		network._props.shortestPaths[source][target]
end

"""
	shortestpath(network::AbstractNetwork)

Return the average shortest path length of the network.
If the network is disconnected, return the average path length of
the largest connected component.
"""
function shortestpath(network::AbstractNetwork)
	# Ensures all shortest paths are calculated
	for i in 1:network.N
		shortestpath(network, i)
	end

	shortestPaths = @view network._props.shortestPaths[:]
	totalLength = 0

	# If the network is disconnected, find the largest connected component
	if any([any(isinf, sp_neuron_i) for sp_neuron_i in shortestPaths])
		components = findComponents(shortestPaths)
		largestComp = components[ findmax(length.(components))[2] ]

		(length(largestComp) == 1) && (return Inf)
	else
		# If the network is not disconnected, the largest component is the whole net
		largestComp = 1:network.N
	end

	for i in largestComp, j in largestComp
		(i == j) && continue
		totalLength += shortestPaths[i][j]
	end

	return totalLength / (length(largestComp) * (length(largestComp) - 1))
end

"""
	averagepathlength(network::AbstractNetwork)

Return the average shortest path length of the network.
"""
averagepathlength(network::AbstractNetwork) = shortestpath(network)

"""
	sigma(network::AbstractNetwork)

Compare the network with an equivalent random network and return
the small-world-ness measure ``σ = (C / C_ran) / (L / L_ran)``

See also: `omega`, `smallworldness`
"""
function sigma(network::AbstractNetwork)
	isnothing(network._props.σ) && smallworldness(network, verbose=false)
	return network._props.σ
end

"""
	omega(network::AbstractNetwork)

Compare the network with an equivalent random network and return
the small-world-ness measure ``ω = (L_ran / L) - (C / C_reg)``

See also: `sigma`, `smallworldness`
"""
function omega(network::AbstractNetwork)
	isnothing(network._props.ω) && smallworldness(network, verbose=false)
	return network._props.ω
end

"""
	smallworldness(network::AbstractNetwork; verbose::Bool=true)

Return the small-world-ness measures σ and ω.
If `verbose` is true, also print an analysis.

See also: `sigma`, `omega`
"""
function smallworldness(network::AbstractNetwork; verbose::Bool=true)
	equivRan = equivalentRandomNetwork(network, seed=network.seed)
	equivLat = equivalentLatticeNetwork(network, seed=network.seed)

	net_cc     = clusteringcoefficient(network)
	random_cc  = clusteringcoefficient(equivRan)
	lattice_cc = clusteringcoefficient(equivLat)

	net_apl    = averagepathlength(network)
	random_apl = averagepathlength(equivRan)
	verbose && (lattice_apl = averagepathlength(equivLat))

	network._props.σ = (net_cc / random_cc) / (net_apl / random_apl)
	network._props.ω = random_apl / net_apl - net_cc / lattice_cc

	if verbose
		println("Small-world-ness analysis")
		println("-------------------------")
		println("Clutering Coefficient:")
		println("- Analysed network:           $(net_cc)")
		println("- Equivalent random network:  $(random_cc)")
		println("- Equivalent lattice network: $(lattice_cc)")
		println("")
		println("Average Path Length")
		println("- Analysed network:           $(net_apl)")
		println("- Equivalent random network:  $(random_apl)")
		println("- Equivalent lattice network: $(lattice_apl)")
		println("")
		println("σ = $(network._props.σ)")
		println("ω = $(network._props.ω)")
		println("")
		println("References")
		println("σ > 1  =>  Small-world")
		println("ω -> -1  => Regular / Lattice")
		println("ω ->  0  => Small-world")
		println("ω ->  1  => Random")
	end

	return network._props.σ, network._props.ω
end

#---------- Calculators ----------
function calcNumConnections!(network::AbstractNetwork)
	network._props.numConnections =
		Int( sum(adjMat(network, sparse=true)) / (isdirected(network) ? 1 : 2) )
end

function calcConnectivity(network::AbstractNetwork, idx_node::Integer; dir_behaviour::Symbol)
	isdirected(network) || (return sum(adjMat(network, sparse=true)[:,idx_node]))

	mat = adjMat(network, sparse=true)
	inDegree = sum(mat[idx_node,:])
	outDegree = sum(mat[:,idx_node])

	(dir_behaviour == :total) && (return inDegree + outDegree)
	(dir_behaviour == :in)    && (return inDegree)
	(dir_behaviour == :out)   && (return outDegree)
	(dir_behaviour == :mean)  && (return (inDegree + outDegree) / 2)
	(dir_behaviour == :both)  && (return (inDegree, outDegree))
	if dir_behaviour == :bi
		biConnecs = Bool.(mat[idx_node, :]) .& Bool.(mat[:, idx_node])
		(isa(biConnecs, SparseArrays.AbstractSparseVector) && isempty(biConnecs.nzind)) &&
			(return 0)
		return sum(biConnecs)
	end
end

function calcAdjMat!(network::AbstractNetwork)
	error("No adjacency matrix calculators were found for a network of type $(typeof(network))")
end

function calcClusteringCoefficient(network::AbstractNetwork, idx_node::Integer)
	# As in DOI: 10.1103/PhysRevE.76.026107
	if isdirected(network)
		nb = neighbors(network, idx_node; directed_behaviour=:any)
		totalDegree = connectivity(network, idx_node, dir_behaviour=:total)
		biDegree = connectivity(network, idx_node, dir_behaviour=:bi)

		numPossibleTriangles = totalDegree * (totalDegree - 1) - 2 * biDegree
		numExistingTriangles = 0

		for n1 in nb, n2 in nb
			(n1 == n2) && continue

			hasconnection(network, idx_node, n1) && hasconnection(network, n1, n2) &&
				hasconnection(network, n2, idx_node) && (numExistingTriangles += 1)

			hasconnection(network, idx_node, n1) && hasconnection(network, n1, n2) &&
				hasconnection(network, idx_node, n2) && (numExistingTriangles += 1)

			hasconnection(network, idx_node, n1) && hasconnection(network, n2, n1) &&
				hasconnection(network, n2, idx_node) && (numExistingTriangles += 1)

			hasconnection(network, n1, idx_node) && hasconnection(network, n1, n2) &&
				hasconnection(network, n2, idx_node) && (numExistingTriangles += 1)
		end
	else
		nb = neighbors(network, idx_node)
		numNeighbors = length(nb)

		numPossibleTriangles = numNeighbors * (numNeighbors - 1) / 2
		numExistingTriangles = 0

		for i in 1:numNeighbors for j in i+1:numNeighbors
			hasconnection(network, nb[i], nb[j]) && (numExistingTriangles += 1)
		end end
	end

	return numExistingTriangles == 0 ? 0 : numExistingTriangles / numPossibleTriangles
end

function calcClusteringCoefficients!(network::AbstractNetwork)
	network._props.clusteringCoefficients =
		[calcClusteringCoefficient(network, i) for i in 1:network.N]
end

function calcTransitivity!(network::AbstractNetwork)
	numExistingTriangles = 0
	numPossibleTriangles = 0

	for idx_node in 1:network.N
		nb = neighbors(network, idx_node)
		numNeighbors = length(nb)

		numPossibleTriangles += numNeighbors * (numNeighbors - 1) / 2

		for i in 1:numNeighbors for j in i+1:numNeighbors
			hasconnection(network, nb[i], nb[j]) && (numExistingTriangles += 1)
		end end
	end

	network._props.transitivity =
		numExistingTriangles == 0 ? 0 : numExistingTriangles / numPossibleTriangles
end

function calcShortestPath!(
	network::AbstractNetwork,
	source::Integer,
	target::Union{Integer,Nothing}=nothing
)
	hasnodeOrError(network, source)
	isnothing(target) || hasnodeOrError(network, target)

	network._props.shortestPaths[source] = dijkstra(network, source, target)
end

function dijkstra(
	network::AbstractNetwork,
	source::Integer,
	target::Union{Integer,Nothing}=nothing
)
	unvisited = Vector(1:network.N)
	distance = repeat([Inf], network.N)
	distance[source] = 0

	while !isempty(unvisited) && !all(x -> isinf(x), unvisited)
		currentNode = unvisited[ findmin(distance[unvisited])[2] ]
		unvisited = filter(x -> x != currentNode, unvisited)

		if !isnothing(target) && currentNode == target
			return distance[currentNode]
		end

		for nb in neighbors(network, currentNode, directed_behaviour=:destination)
			testVal = distance[currentNode] + 1
			(testVal < distance[nb]) && (distance[nb] = testVal)
		end
	end

	return distance
end

function findComponents(shortestPaths)
	N = length(shortestPaths)
	components = []

	while isempty(components) || sum(length.(components)) < N
		component = Vector(1:N)
		for found_comp in components
			component = filter(x -> !(x in found_comp), component)
		end
		out = []

		for i in component
			if isinf(shortestPaths[i][component[1]])
				push!(out, i)
			end
		end

		component = filter(x -> !(x in out), component)
		out = []

		for i in component
			for j in component
				if isinf(shortestPaths[i][j])
					push!(out, j)
				end
			end
		end

		component = filter(x -> !(x in out), component)
		push!(components, component)
	end

	return components
end


text:
using AnsiColor
using Base.Test

println(black("black"))
@test black("Foo") == "\e[0;30;49mFoo\e[0m"
@test black("Foo") == colorize(:black, "Foo")

println(red("red"))
@test red("Foo") == "\e[0;31;49mFoo\e[0m"
@test red("Foo") == colorize(:red, "Foo")

println(green("green"))
@test green("Foo") == "\e[0;32;49mFoo\e[0m"
@test green("Foo") == colorize(:green, "Foo")

println(yellow("yellow"))
@test yellow("Foo") == "\e[0;33;49mFoo\e[0m"
@test yellow("Foo") == colorize(:yellow, "Foo")

println(blue("blue"))
@test blue("Foo") == "\e[0;34;49mFoo\e[0m"
@test blue("Foo") == colorize(:blue, "Foo")

println(magenta("magenta"))
@test magenta("Foo") == "\e[0;35;49mFoo\e[0m"
@test magenta("Foo") == colorize(:magenta, "Foo")

println(cyan("cyan"))
@test cyan("Foo") == "\e[0;36;49mFoo\e[0m"
@test cyan("Foo") == colorize(:cyan, "Foo")

println(white("white"))
@test white("Foo") == "\e[0;37;49mFoo\e[0m"
@test white("Foo") == colorize(:white, "Foo")

println(uncolor("default"))
@test uncolor("Foo") == "\e[0;39;49mFoo\e[0m"
@test uncolor("Foo") == colorize(:default, "Foo")


text:
using Plots, DifferentialEquations
dt = 0.001
σnoise = 0.007
τAMPA = 2e-3
RNG = MersenneTwister(10)
RNG2 = MersenneTwister(10)

function model1(du, u, p, t, W)
    x, I = u
    du[1] = 0.0
    du[2] = (W[1]*σnoise*√τAMPA - I) / τAMPA
end

prob = RODEProblem(model1, [0.0, 0.0], (0.0, 1.0), noise=WienerProcess(0.0, 0.0, rng=RNG, reseed=false))
sol = solve(prob, saveat=dt, dt=dt, reltol=1e-15, abstol=1e-15)
display(plot(sol))

function model2(du, u, p, t, W)
    du[1] = 2*u[1]*sin(W[1])
end

W = OrnsteinUhlenbeckProcess(τAMPA, 0.0, σnoise*√τAMPA, 0.0, [0.0], rng=RNG2, reseed=false)
prob2 = RODEProblem(model2, [0.0], (0.0, 1.0), noise=W)
sol2 = solve(prob2, saveat=dt, dt=dt, reltol=1e-15, abstol=1e-15)
display(plot(sol2.W.t, reduce(vcat, sol2.W.u)))


text:
using StatsPlots
using Test
using StableRNGs
using NaNMath

@testset "Grouped histogram" begin
    rng = StableRNG(1337)
    gpl = groupedhist( rand(rng, 1000), yscale=:log10, ylims=(1e-2, 1e4), bar_position = :stack)
    @test NaNMath.minimum(gpl[1][1][:y]) <= 1e-2
    @test NaNMath.minimum(gpl[1][1][:y]) > 0
    rng = StableRNG(1337)
    gpl = groupedhist( rand(rng, 1000), yscale=:log10, ylims=(1e-2, 1e4), bar_position = :dodge)
    @test NaNMath.minimum(gpl[1][1][:y]) <= 1e-2
    @test NaNMath.minimum(gpl[1][1][:y]) > 0

    data = [1,1,1,1,2,1]
    mask = (collect(1:6) .< 5)
    gpl1 = groupedhist(data[mask], group = mask[mask], color = 1)
    gpl2 = groupedhist(data[.!mask], group = mask[.!mask], color = 2)
    gpl12 = groupedhist(data, group = mask, nbins = 5, bar_position = :stack)
    @test NaNMath.maximum(gpl12[1][2][:y]) == NaNMath.maximum(gpl1[1][1][:y])
    data=[10 12; 1 1; 0.25 0.25]
    gplr = groupedbar(data)
    @test NaNMath.maximum(gplr[1][1][:y]) == 10
    @test NaNMath.maximum(gplr[1][2][:y]) == 12
    gplr = groupedbar(data, bar_position = :stack)
    @test NaNMath.maximum(gplr[1][1][:y]) == 22
    @test NaNMath.maximum(gplr[1][2][:y]) == 12
end # testset


text:
MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.inner.variable_info)
MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices) = VI.(1 : length(model.inner.variable_info))

function MOIU.allocate_variables(model::Optimizer, nvars::Integer)
    variable_info = model.inner.variable_info
    resize!(variable_info, nvars)
    @inbounds for i = 1 : nvars
        variable_info[i] = VariableInfo()
    end
    return VI.(1 : nvars)
end

MOIU.load_variables(model::Optimizer, nvars::Integer) = nothing

MOI.supports(model::Optimizer, ::MOI.VariableName, ::Type{VI}) = true

function MOIU.load(model::Optimizer, attr::MOI.VariableName, vi::VI, value)
    check_variable_indices(model, vi)
    model.inner.variable_info[vi.value].name = value
end

MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{VI}) = true

function MOIU.load(model::Optimizer, ::MOI.VariablePrimalStart, vi::VI, value::Union{Real, Nothing})
    check_variable_indices(model, vi)
    model.inner.variable_info[vi.value].start = value
    return
end


text:
# unit test the bGrad

include("../src/DMD_Util.jl");
include("../src/DMD_Type.jl");
include("../src/DMD_Func.jl");

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
k = 3;      # number of modes
T = Float64;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
#------------------------------------------------------------------------------
params = DMDParams(k, X, t, lossFunc, lossGrad);
#
# Initialize with true parameter
#------------------------------------------------------------------------------
copyto!(params.a, at);
copyto!(params.B, Bt);
#
# Evaluate bGrad
#------------------------------------------------------------------------------
update_P!(params);
id = rand(1:n);
function func(br, params, id)
	copyto!(params.br[id], br);
	return bFunc(params, id)
end
func(br) = func(br, params, id);
#
function grad(gbr, br, params, id)
	copyto!(params.br[id], br);
	bGrad(gbr, params, id);
end
grad(gbr, br) = grad(gbr, br, params, id);
#
br = randn(k<<1);
gbr = zeros(k<<1);
grad(gbr, br);
#
# finite difference method verification
ϵ = 1e-8;
br_fd = zeros(k<<1);
gbr_fd = zeros(k<<1);
#
for i in eachindex(gbr_fd)
	copyto!(br_fd, br);
	br_fd[i] += ϵ;
	gbr_fd[i] = (func(br_fd) - func(br))/ϵ;
end
#
err = sum(abs2, gbr_fd - gbr);
if err < 1e-6
	println("bGrad: OK");
else
	println("bGrad: Wrong, err: $err");
end


text:
# essa é uma análise exploratória (EDA) da base pública de CNPJs (QSA) da Receita Federal do Brasil (RFB)

## carrega bibliotecas relevantes
using DataFrames # para manipulação de data frames
#using CSV # para ler CSVs
using SQLite, Tables # para ler SQLites
using StatsBase, Statistics # para estatística básica
using Queryverse # metapacote para data science
using Gadfly

## lendo arquivo de entrada
db = SQLite.DB("bd_dados_qsa_cnpj.db") # carrega o arquivo SQLite com as tabelas

## cria índices com colunas úteis para consultas
# SQLite.createindex!(db, "cnpj_dados_cadastrais_pj", "idx_situacao_cadastral", "situacao_cadastral",
#     unique=false, ifnotexists=true)

# SQLite.createindex!(db, "cnpj_dados_cadastrais_pj", "idx_identificador_matriz_filial", "identificador_matriz_filial",
#     unique=false, ifnotexists=true)

# SQLite.createindex!(db, "cnpj_dados_cadastrais_pj", "idx_porte_empresa", "porte_empresa",
#     unique=false, ifnotexists=true)

# SQLite.createindex!(db, "cnpj_dados_cadastrais_pj", "idx_opcao_pelo_simples", "opcao_pelo_simples",
#     unique=false, ifnotexists=true)

# SQLite.createindex!(db, "cnpj_dados_cadastrais_pj", "idx_opcao_pelo_mei", "opcao_pelo_mei",
#     unique=false, ifnotexists=true)

# SQLite.createindex!(db, "cnpj_dados_cadastrais_pj", "idx_cnpj", "cnpj",
#     unique=true, ifnotexists=true)

## executa a consulta base e retorna um result set, que será depois usado
# para trazer as empresas de interesse
# IMPORTANTE: como a base é muito grande, para reduzir o uso massivo de memória RAM
# é crucial ter uma consulta base que já filtre as empresas de principal interesse -
# no nosso caso só estamos interessados nas empresas ativas (situacao_cadastral='02'),
# matrizes (identificador_matriz_filial = '1'), e sem considerar MEIs (opcao_pelo_mei = 'N')
# - e restrinja os campos de interesse;
# como teremos de fazer já a consulta, vamos aproveitar para fazer joins e transformações
# interessantes
dados_pj_rs =
    DBInterface.execute(db,
"
    SELECT
        pj.cnpj,
        pj.nome_fantasia,
        pj.situacao_cadastral as situacao,
        case pj.porte_empresa
            when '00' then '?'
            when '01' then 'MICRO'
            when '03' then 'PEQUENA'
            when '05' then 'DEMAIS'
        end as porte,
        case pj.opcao_pelo_simples
            when '' then 'N'
            when '0' then 'N'
            when '5' then 'S'
            when '7' then 'S'
            when '6' then 'S - Excluido'
            when '8' then 'S - Excluido'
        end as simples,
        case pj.opcao_pelo_mei
            when 'S' then 'S'
            else 'N'
        end as MEI,
        nat_jur.nm_natureza_juridica as nat_juridica,
        pj.data_inicio_atividade as dat_ini_atividade,
        pj.cnae_fiscal as cod_cnae_principal,
        cnae.nm_cnae as cnae_principal,
        cnae.nm_divisao as divisao_cnae_principal,
        cnae.nm_classe as classe_cnae_principal,
        cnae.nm_grupo as grupo_cnae_principal,
        pj.bairro, pj.codigo_municipio as cod_municipio, pj.municipio, pj.uf
    FROM `cnpj_dados_cadastrais_pj` pj
        INNER JOIN `tab_cnae` cnae on cnae.cod_cnae = pj.cnae_fiscal
        INNER JOIN `tab_natureza_juridica` nat_jur on nat_jur.cod_subclass_natureza_juridica = pj.codigo_natureza_juridica
    WHERE
        situacao_cadastral='02' -- somente ativas
        and identificador_matriz_filial='1' -- somente matrizes
        and pj.opcao_pelo_mei = 'N' -- não considerar MEIs (muito numerosos, não são empresas, mas empresários individuais)
    ;
") # consulta base

# OBS: o result set somente suporta ser consumido uma única vez (single pass, forware only).
# Para usá-lo de novo, é preciso resetá-lo (rodar a query de novo) com: SQLite.reset!(dados_pj_rs)

## tendo o result set inicial, podemos executar as análises de interesse

### configurações gerais de plot
include("estilo_plot.jl")
Gadfly.push_theme(estilo_plot())
set_default_plot_size(21cm, 10cm)

## Análise 1: quantos cartórios temos no Brasil?
SQLite.reset!(dados_pj_rs) # para recarregar o result set, se necessário
dados_cartorios = dados_pj_rs |>
    @filter(_.cnae_principal=="Cartórios") |>
    DataFrame

print("O Brasil possui $(nrow(dados_cartorios)) cartórios!")

# mas onde eles estão localizados? são muitos porque temos muitas cidades?
dados_cartorios_uf = dados_cartorios |>
    @groupby(_.uf) |> @map({uf=key(_), n_cartorios=length(_), n_cartorios_str=string(length(_)) }) |> @orderby_descending(_.n_cartorios)

dados_cartorios_uf_municipios = dados_cartorios |>
    @select(:uf, :municipio) |> @unique() |>
    @groupby(_.uf) |> @map({uf=key(_), n_municipios=length(_), n_municipios_str=string(length(_))}) |> @orderby_descending(_.n_municipios)

comparacao_cartorios_uf = dados_cartorios_uf |>
    @join(dados_cartorios_uf_municipios, _.uf, _.uf, {_.uf, _.n_cartorios, _.n_cartorios_str, __.n_municipios, __.n_municipios_str}) |>
    @orderby_descending(_.n_cartorios) |> DataFrame

print(comparacao_cartorios_uf)

set_default_plot_size(30cm, 18cm)

p1 = plot(comparacao_cartorios_uf |> @take(10) |> DataFrame,
    y=:uf, x=:n_cartorios, label=:n_cartorios_str,
    Geom.bar(orientation=:horizontal), Geom.label(position=:right),
    Scale.x_continuous(format=:plain),
    Guide.title("Top 10 UFs em nr cartórios"), Guide.xlabel("Cartórios"), Guide.ylabel("UF"))

p2 = plot(comparacao_cartorios_uf |> @orderby_descending(_.n_municipios) |> @take(10) |> DataFrame,
    y=:uf, x=:n_municipios, label=:n_municipios_str,
    Geom.bar(orientation=:horizontal), Geom.label(position=:right),
    Scale.x_continuous(format=:plain),
    Guide.title("Top 10 UFs em nr municípios"),  Guide.xlabel("Municípios"), Guide.ylabel("UF"))

vstack(p1, p2)

max_stars_repo_name: llrt/eda-qsa-rfb-jl

max_issues_repo_name: llrt/eda-qsa-rfb-jl

max_forks_repo_name: llrt/eda-qsa-rfb-jl


text:
#!/usr/bin/env julia
using FASTX
using BioAlignments

# Read sequence dictionary for HLA
seqFile = "data/S2_callseqs.fa"
name2seq = Dict()
open(FASTA.Reader,seqFile) do reader
	for record in reader
		name2seq[FASTA.identifier(record)]=FASTA.sequence(record)
	end
end

#target = "HLA:HLA00602" # >HLA:HLA00602 DQA1*01:02:01:01 6484 bp
#chrStart=32628178 # DQA1
#chrEnd=32637531
target = "HLA:HLA00010" # >HLA:HLA00010 A*02:05:01:01 3517 bp
chrStart = 29941259
chrEnd = 29945884

hlaseq = name2seq[target]

# Read chromosome 6
chr6File= "data/chr6_hg38.fa"
chr6record = FASTA.Record()
reader = open(FASTA.Reader,chr6File)
read!(reader,chr6record)
close(reader)

# Pull out region of chromosome6
chr6seq = FASTA.sequence(chr6record)
hlaref = chr6seq[chrStart:chrEnd]

# Align that region to big HLA sequence

for j in 0:40
	scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=-1);
	res = pairalign(LocalAlignment(),hlaref,hlaseq,scoremodel)
end

#print(res)

# HLA-A
# ./cutandalign.jl  2.69s user 0.54s system 104% cpu 3.080 total

# 10X loop:
# cutandalign2.jl  3.48s

# 20X loop:
# cutandalign2.jl  4.54s

# 40X ./cutandalign2.jl  6.65s
# So 2sec = 20X extra so 10X/second steady state.

# HLA-DQA1
# ./cutandalign.jl  2.58s user 0.86s system 101% cpu 3.398 total

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
using CSV
using TimeSeries
using DataFrames

###### native data loader #################

function _load_data(source::Val{:native}, region::AbstractString, path::AbstractString)
    # fn = _build_filename(region)
    fn = region
    fp = rstrip(path,'/')*"/"*fn*".csv"
    if fn != region
        @warn("loading data from `"*fp*"'")
    end
    header = ["Date", "Confirmed", "Deaths"];
    df = CSV.read(fp,header=header,datarow=2,dateformat="yyyy-mm-dd",ignoreemptylines=true)
    TimeArray(df.Date,Array{Int64}(df[!,2:3]),Symbol.(header[2:3]))
end


text:
import NLsolve: nlsolve
import ReverseDiff
import LinearAlgebra: tr

export OrderedBasis, TransportFamily, TransportMap, inverse, sample_cost

"""
    abstract type OrderedBasis end

Abstract type corresponding to an ordered basis of scalar-functions
``(ψ_j)_{j\\in\\mathcal J_i}``.  An instance `tuple` of this type will act as a map, so that
`tuple(θ)` will correspond to the vector ``(ψ_j(θ))_{j \\in \\mathcal J_i}``.

To implement a subtype `B <: OrderedBasis`, one must implement:

```julia
(tuple::B)(θ::AbstractArray{<:Real})::AbstractArray{<:Real}
Base.length(tuple::B)::Int64
```
"""
abstract type OrderedBasis end

function Base.length(tuple::OrderedBasis)
  error("must implement Base.length(::$(typeof(tuple)))")
end

"""
    T = TransportFamily(bases::AbstractArray{<:OrderedBasis})

Provide an array `bases` of `OrderedBasis` instances where each element `bases[i]`
corresponds to some ordered basis ``(ψ_j)_{j\\in\\mathcal J_i}`` and receive the induced
family ``\\{ \\tilde T(\\cdot; γ) \\}_{γ \\in \\mathcal A}`` (`T`).

With an instance `T` of a `TransportFamily`, one may apply `T(γ)` to a parameter `γ` and get a
`TransportMap` instance corresponding to the map ``T(\\cdot; γ)``.
Recall, this is defined so that the ``i``-th component ``T_i(\\cdot; γ_i)`` is as follows.

```math
T_i(θ; γ_i) = \\sum_{j\\in\\mathcal J_i} γ_{i,j} ψ_j(θ)
```
"""
struct TransportFamily{B <: OrderedBasis}
  bases::AbstractArray{B}
end

function Base.display(T::TransportFamily)
  println(
    """
    TransportFamily on R^$(length(T.bases))
    $(join(map(p -> "  J_$(p[1]) = $(p[2])", enumerate(T.bases)), "\n"))"""
  )
end

"""
    map = TransportMap(dim, parameter, f)

Wrap a provided function `f` with a `TransportMap` type for the type checking.
In addition, we keep track of the dimension (`dim`) of the vector space on which `f` is
defined and the parameter (`parameter`) of the specific map.
"""
struct TransportMap{B <: OrderedBasis}
  dim::Int
  parameter::AbstractArray{<:AbstractArray{<:Real}}
  f::Function
end

function (F::TransportMap)(theta)::AbstractArray{<:Real}
  if length(theta) !== F.dim
    throw(ArgumentError("$theta must have length as $(F.dim)"))
  end
  r = F.f(theta)
  if length(r) !== F.dim
    println("Warning: the codomain does not have the appropriate dimension")
  end
  r
end

function Base.display(F::TransportMap{B}) where B <: OrderedBasis
  println(
    """
    TransportMap{$B}
    $(join(map(p -> "  γ_$(p[1]) = $(p[2])", enumerate(F.parameter)), "\n"))"""
  )
end

function (T::TransportFamily{B})(gamma::AbstractArray{<:AbstractArray{<:Real}})::TransportMap{B} where B
  # make sure there is a parameter for each component function
  if (length(gamma) !== length(T.bases))
    throw(ArgumentError("$(T.bases) must have same length as $gamma"))
  end
  d = length(gamma)

  # span the components
  TransportMap{B}(d, gamma, theta -> begin
    map(pair -> begin
      # test if the i-th component basis has the same length as the matching parameter
      basis, gam_i = pair
      tuple = basis(theta)
      if length(gam_i) !== length(tuple)
        throw(ArgumentError(
          "The ordered basis $(basis) must have the same amount of components as $(gam_i)"
        ))
      end

      # take the linear combination
      reduce((s, p) -> s + p[1] * p[2], zip(gam_i, tuple); init=0.0)
    end, zip(T.bases, gamma))
  end)
end

"""
    theta = inverse(tmap::TransportMap, r)

Return ``θ`` such that ``T(θ; γ) = r``, where ``T(\\cdot; γ)`` is our transport map
(`tmap`).
"""
function inverse(tmap::TransportMap, r::AbstractArray{<:Real})::AbstractArray{<:Real}
  if length(r) !== tmap.dim
    throw(ArgumentError("$r must have length $(tmap.dim)"))
  end
  reduce(
    (theta, pair) -> begin
      # unpack the index and image component
      i, r_i = pair

      # preallocate the theta components we have solved and those that are unused
      theta_solved = theta[1:(i-1)]
      theta_unused = theta[(i+1):end]

      # create the partially-evaluated component function, ready for NLsolve
      func! = (F, x) -> begin
        F[1] = tmap([theta_solved..., x[1], theta_unused...])[i] - r_i
      end

      # solve the variable and return the update
      theta_i = nlsolve(func!, [0.0]).zero[1]
      [theta_solved..., theta_i, theta_unused...]
    end,
    enumerate(r);
    init=zeros(tmap.dim)
  )
end

"""
    cost = sample_cost(tmap::TransportMap, samples::AbstractArray{<:AbstractArray{<:Real}})

Compute the cost ``C(\\tilde T(\\cdot; γ))`` of a transport map
``\\tilde T(\\cdot; γ)`` based off of samples ``θ_1, \\ldots, θ_K``.

```math
C\\big(\\tilde T(\\cdot; γ)\\big) = \\sum_{i=1}^d\\sum_{k=1}^K \\Big(
  \\frac{1}{2} \\tilde T_i^2(θ^{(k)}; γ_i) - \\log\\frac{∂\\tilde T_i(\\cdot; γ_i)}{∂θ_i}(θ^{(k)})
\\Big)
```
"""
function sample_cost(tmap::TransportMap, samples::AbstractArray{<:AbstractArray{<:Real}})
  sum(
    map(
      sample -> 0.5 * sum(tmap(sample).^2) - tr(log.(ReverseDiff.jacobian(tmap, sample))),
      samples,
    )
  )
end

#---------------------------------------------------------------------------------------------
# specific implementations
#---------------------------------------------------------------------------------------------
include("polynomials.jl")

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
function simulate(nsamples :: Int, params :: MCFMParams, factor_func, error_func)
    if nsamples <= 0
        throw(ArgumentError("Number of samples must be positive"))
    end
    factorlayout = params.factorlayout
    channellayout = params.channellayout
    factors = factor_func(nsamples, factorlayout) :: MultiChannelFactors
    #TODO : errors should take Sigma from params
    errors  = error_func(nsamples, channellayout) :: MultiChannelData
    exactvalues  = params.loading * factors.factors
    exactdata  = unstack(exactvalues, channellayout)
    data = exactdata + errors
    return data
end

function periodic_factors(nsamples, factorlayout::MultiChannelFactorLayout)
    freqs = 2*π .* (1:factorlayout.nfactors)
    factors = transpose(hcat([sin.((1:nsamples)/nsamples * f) for f in freqs ]...))
    fcommon = factors[1:factorlayout.nchannelshared, :]
    funique = factors[(factorlayout.nchannelshared+1):end, :]
    out = MultiChannelFactors(factorlayout, fcommon, funique)
    return out
end

function classical_factors(nsamples, factorlayout::MultiChannelFactorLayout)
    # Independent N(0, 1) for all samples is classical
    dist = Normal(0, 1)
    factors = rand(dist, (factorlayout.nfactors, nsamples))
    fcommon = factors[1:factorlayout.nchannelshared, :]
    funique = factors[(factorlayout.nchannelshared+1):end, :]
    out = MultiChannelFactors(factorlayout, fcommon, funique)
    return out
end

function indep_gauss_errors(nsamples, channellayout::MultiChannelLayout; sd=0.1)
    dist = Normal(0, sd)
    errors = unstack(rand(dist, (channellayout.nobs, nsamples)), channellayout)
    return errors
end

function simChannel(H, G, fc, fu, edist::Distribution)
    if length(edist) == 1
        # Univariate dist, sample iid
        errs = rand(edist, size(H)[1])
    elseif length(edist) == size(H)[1]
        # Multivariate dist, sample once
        errs = rand(edist, 1)
        errs = errs[:, 1]
    else
        throw(ArgumentError("Distribution dim invalid"))
    end
    x = H * fc + G * fu +  errs
    return(x)
end

function simFactor(ncommon, nuniques)
    dist = Normal(0, 1) # Classical FA: factors iid N(0, 1)
    fc = rand(dist, ncommon)
    fus = [rand(dist, n) for n in nuniques]
    return fc, fus
end

function simACData(N, Hs, Gs, edists; ρ=0.9)
    if length(Hs) != length(Gs) || length(Hs) != length(edists)
        throw(ArgumentError("Length mismatch"))
    end
    nc = size(Hs[1])[2]
    nus = [size(G)[2] for G in Gs ]
    dat = Array{Array{Array{Float64, 1}, 1}, 1}(undef, N)
    for n in 1:N
        fc, fus = simFactor(nc, nus) # Almost certain we want a new factor draw for each sample
        xs = Array{Array{Float64, 1}, 1}(undef, length(Hs))
        for i in 1:length(Hs)
            x = simChannel(Hs[i], Gs[i], fc, fus[i], edists[i])
            xs[i] = x
        end
        dat[n] = xs
    end
    return(dat)
end

function simData(N, Hs, Gs, edists; factor_func = simFactor)
    if length(Hs) != length(Gs) || length(Hs) != length(edists)
        throw(ArgumentError("Length mismatch"))
    end
    nc = size(Hs[1])[2]
    nus = [size(G)[2] for G in Gs ]
    dat = Array{Array{Array{Float64, 1}, 1}, 1}(undef, N)
    for n in 1:N
        fc, fus = simFactor(nc, nus) # Almost certain we want a new factor draw for each sample
        xs = Array{Array{Float64, 1}, 1}(undef, length(Hs))
        for i in 1:length(Hs)
            x = simChannel(Hs[i], Gs[i], fc, fus[i], edists[i])
            xs[i] = x
        end
        dat[n] = xs
    end
    return(dat)
end

function randomLoadingNorm(L, p; target_tr=1)
    top = rand(Normal(0, 1), p, p)
    Q1 = LowerTriangular(top)
    Q2 = rand(Normal(0, 1), (L-p), p)
    Q = vcat(Q1, Q2)
    trace = tr(Q * transpose(Q))
    Q = sqrt(target_tr/trace) * Q
    diag_signs = sign.(diag(Q))
    Q= Q * diagm(diag_signs)
    return Q
end

function randomLoadingUnif(L, p)
    dist = DiscreteUniform(-2, 2)
    Q = rand(dist, (L, p))
    return(Q)
end

function randomLoadings(Ls, nc, nus, loadingfunc)
    nchan = length(Ls)
    Hs = [loadingfunc(L, nc) for L in Ls]
    Gs = [loadingfunc(L, nu) for (L, nu) in zip(Ls, nus)]
    return Hs, Gs
end

function LyuouLoadings1(Ls, nc, nus)
    N = sum(Ls)
    Hs = [ transpose(hcat([[0.8^f * (-1)^(i+f) for f in 1:nc] for i in 1:L]...)) for L in Ls]
    Gs = [ones(L, nu) for (L, nu) in zip(Ls, nus)]
    Gs = [(length(Ls) - i + 1) .* Gs[i] for i in 1:length(Gs)]
    Gs[1] = Gs[1] .+ [i/(2*Ls[1]) for i in 1:Ls[1]]
    Σs = [diagm(repeat([1], L)) for L in Ls ]
    return Hs, Gs, Σs
end

function powerStructLoadings(Ls, nc, nus;commonpow=0.3, uniquepow=0.3, sigmapow=0.4)
    Hs = extract_blocks(randomLoadingNorm(sum(Ls), nc), Ls, nc; diag=false)
    Gs = [randomLoadingNorm(L, nu) for (L, nu) in zip(Ls, nus)]
    Σs = [diagm(repeat([4/(3*L)], L)) for L in Ls ]

    return Hs, Gs, Σs
end


text:
"Type-safe wrapper for `Int64`, references a variable."
struct VarRef
    val::Int64
end

"Type-safe wrapper for `Int64`, references a constraint."
struct ConsRef
    val::Int64
end

"ManagedSCIP holds pointers to SCIP data and takes care of memory management."
mutable struct ManagedSCIP
    scip::Ref{Ptr{SCIP_}}
    vars::Dict{VarRef, Ref{Ptr{SCIP_VAR}}}
    conss::Dict{ConsRef, Ref{Ptr{SCIP_CONS}}}
    var_count::Int64
    cons_count::Int64

    # Map from user-defined Julia types (keys are <: AbstractConstraintHandler
    # or <: AbstractConstraint, respectively) to the corresponding SCIP objects.
    # The reverse mapping is handled by SCIP itself.
    # This also serves to prevent premature GC.
    conshdlrs::Dict{Any, Ptr{SCIP_CONSHDLR}}
    conshdlrconss::Dict{Any, Ptr{SCIP_CONS}}

    # Map from user-defined types (keys are <: AbstractSeparator) to the
    # corresponding SCIP objects.
    sepas::Dict{Any, Ptr{SCIP_SEPA}}

    function ManagedSCIP()
        scip = Ref{Ptr{SCIP_}}(C_NULL)
        @SC SCIPcreate(scip)
        @assert scip[] != C_NULL
        @SC SCIPincludeDefaultPlugins(scip[])
        @SC SCIP.SCIPcreateProbBasic(scip[], "")

        mscip = new(scip, Dict(), Dict(), 0, 0, Dict(), Dict(), Dict())
        finalizer(free_scip, mscip)
    end
end

# Protect ManagedSCIP from GC for ccall with Ptr{SCIP_} argument.
Base.unsafe_convert(::Type{Ptr{SCIP_}}, mscip::ManagedSCIP) = mscip.scip[]

"Release references and free memory."
function free_scip(mscip::ManagedSCIP)
    # Avoid double-free (SCIP will set the pointers to NULL).
    if mscip.scip[] != C_NULL
        for c in values(mscip.conss)
            @SC SCIPreleaseCons(mscip, c)
        end
        for v in values(mscip.vars)
            @SC SCIPreleaseVar(mscip, v)
        end
        # only mscip.scip is GC-protected during ccall!
        GC.@preserve mscip begin
            @SC SCIPfree(mscip.scip)
        end
    end
    @assert mscip.scip[] == C_NULL
end

"Set a parameter's current value."
function get_parameter(mscip::ManagedSCIP, name::AbstractString)
    param = SCIPgetParam(mscip, name)
    if param == C_NULL
        error("Unrecognized parameter: $name")
    end
    paramtype = SCIPparamGetType(param)
    if paramtype === SCIP_PARAMTYPE_BOOL
        value = Ref{SCIP_Bool}()
        @SC SCIPgetBoolParam(mscip, name, value)
        if value[] == TRUE
            return true
        elseif value[] == FALSE
            return false
        else
            error("encountered invalid value for a boolean: $(value[])")
        end
    elseif paramtype === SCIP_PARAMTYPE_INT
        value = Ref{Cint}()
        @SC SCIPgetIntParam(mscip, name, value)
        return value[]
    elseif paramtype === SCIP_PARAMTYPE_LONGINT
        value = Ref{Clonglong}()
        @SC SCIPgetLongintParam(mscip, name, value)
        return value[]
    elseif paramtype === SCIP_PARAMTYPE_REAL
        value = Ref{Cdouble}()
        @SC SCIPgetRealParam(mscip, name, value)
        return value[]
    elseif paramtype === SCIP_PARAMTYPE_CHAR
        value = Ref{Cchar}()
        @SC SCIPgetCharParam(mscip, name, value)
        return Char(value[])
    elseif paramtype === SCIP_PARAMTYPE_STRING
        value = Ref{Cstring}()
        @SC SCIPgetStringParam(mscip, name, value)
        return unsafe_string(value[])
    else
        error("Unexpected parameter type: $paramtype")
    end
end

"Set a parameter."
function set_parameter(mscip::ManagedSCIP, name::AbstractString, value)
    param = SCIPgetParam(mscip, name)
    if param == C_NULL
        error("Unrecognized parameter: $name")
    end
    paramtype = SCIPparamGetType(param)
    if paramtype === SCIP_PARAMTYPE_BOOL
        @SC SCIPsetBoolParam(mscip, name, value)
    elseif paramtype === SCIP_PARAMTYPE_INT
        @SC SCIPsetIntParam(mscip, name, value)
    elseif paramtype === SCIP_PARAMTYPE_LONGINT
        @SC SCIPsetLongintParam(mscip, name, value)
    elseif paramtype === SCIP_PARAMTYPE_REAL
        @SC SCIPsetRealParam(mscip, name, value)
    elseif paramtype === SCIP_PARAMTYPE_CHAR
        @SC SCIPsetCharParam(mscip, name, value)
    elseif paramtype === SCIP_PARAMTYPE_STRING
        @SC SCIPsetStringParam(mscip, name, value)
    else
        error("Unexpected parameter type: $paramtype")
    end
    return nothing
end

"Return pointer to SCIP variable."
function var(mscip::ManagedSCIP, vr::VarRef)::Ptr{SCIP_VAR}
    return mscip.vars[vr][]
end

"Return pointer to SCIP constraint."
function cons(mscip::ManagedSCIP, cr::ConsRef)::Ptr{SCIP_CONS}
    return mscip.conss[cr][]
end

"Store reference to variable, return VarRef"
function store_var!(mscip::ManagedSCIP, var__::Ref{Ptr{SCIP_VAR}})
    mscip.var_count += 1
    vr = VarRef(mscip.var_count)
    mscip.vars[vr] = var__
    return vr
end

"Store reference to constraint, return ConsRef"
function store_cons!(mscip::ManagedSCIP, cons__::Ref{Ptr{SCIP_CONS}})
    mscip.cons_count += 1
    cr = ConsRef(mscip.cons_count)
    mscip.conss[cr] = cons__
    return cr
end

"Add variable to problem (continuous, no bounds), return var ref."
function add_variable(mscip::ManagedSCIP)
    var__ = Ref{Ptr{SCIP_VAR}}(C_NULL)
    @SC SCIPcreateVarBasic(mscip, var__, "", -SCIPinfinity(mscip), SCIPinfinity(mscip),
                           0.0, SCIP_VARTYPE_CONTINUOUS)
    @SC SCIPaddVar(mscip, var__[])
    return store_var!(mscip, var__)
end

"Delete variable from problem."
function delete(mscip::ManagedSCIP, vr::VarRef)
    # delete variable from SCIP problem
    deleted = Ref{SCIP_Bool}()
    @SC SCIPdelVar(mscip, var(mscip, vr), deleted)
    deleted[] == TRUE || error("Variable at $(vr.val) could not be deleted!")

    # release memory and remove reference
    @SC SCIPreleaseVar(mscip, mscip.vars[vr])
    delete!(mscip.vars, vr)
    return nothing
end

"Delete constraint from problem."
function delete(mscip::ManagedSCIP, cr::ConsRef)
    @SC SCIPdelCons(mscip, cons(mscip, cr))

    # release memory and remove reference
    @SC SCIPreleaseCons(mscip, mscip.conss[cr])
    delete!(mscip.conss, cr)
    return nothing
end

"""
Add (ranged) linear constraint to problem, return cons ref.

# Arguments
- `varrefs::AbstractArray{VarRef}`: variable references for affine terms.
- `coefs::AbstractArray{Float64}`: coefficients for affine terms.
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
"""
function add_linear_constraint(mscip::ManagedSCIP, varrefs, coefs, lhs, rhs)
    @assert length(varrefs) == length(coefs)
    vars = [var(mscip, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateConsBasicLinear(
        mscip, cons__, "", length(vars), vars, coefs, lhs, rhs)
    @SC SCIPaddCons(mscip, cons__[])
    return store_cons!(mscip, cons__)
end

"""
Add (ranged) quadratic constraint to problem, return cons ref.

# Arguments
- `linrefs::AbstractArray{VarRef}`: variable references for affine terms.
- `lincoefs::AbstractArray{Float64}`: coefficients for affine terms.
- `quadrefs1::AbstractArray{VarRef}`: first variable references for quadratic terms.
- `quadrefs2::AbstractArray{VarRef}`: second variable references for quadratic terms.
- `quadcoefs::AbstractArray{Float64}`: coefficients for quadratic terms.
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
"""
function add_quadratic_constraint(mscip::ManagedSCIP, linrefs, lincoefs,
                                  quadrefs1, quadrefs2, quadcoefs, lhs, rhs)
    @assert length(linrefs) == length(lincoefs)
    @assert length(quadrefs1) == length(quadrefs2)
    @assert length(quadrefs1) == length(quadcoefs)

    linvars = [var(mscip, vr) for vr in linrefs]
    quadvars1 = [var(mscip, vr) for vr in quadrefs1]
    quadvars2 = [var(mscip, vr) for vr in quadrefs2]

    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateConsBasicQuadratic(
        mscip, cons__, "", length(linvars), linvars, lincoefs,
        length(quadvars1), quadvars1, quadvars2, quadcoefs, lhs, rhs)
    @SC SCIPaddCons(mscip, cons__[])
    return store_cons!(mscip, cons__)
end

"""
Add second-order-cone constraint to problem, return cons ref.

Does not support the full generality of SCIP's constraint (offsets and
coefficients). The first entry in `varrefs` is used for the special variable on
the right-hand side.
"""
function add_second_order_cone_constraint(mscip::ManagedSCIP, varrefs)
    vars = [var(mscip, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateConsBasicSOC(mscip, cons__, "", length(vars) - 1,
                               vars[2:end], C_NULL, C_NULL, 0.0, vars[1], 1.0, 0.0)
    @SC SCIPaddCons(mscip, cons__[])
    return store_cons!(mscip, cons__)
end

"""
Add special-ordered-set of type 1 to problem, return cons ref.

# Arguments
- `varrefs::AbstractArray{VarRef}`: variable references
- `weights::AbstractArray{Float64}`: numeric weights
"""
function add_special_ordered_set_type1(mscip::ManagedSCIP, varrefs, weights)
    @assert length(varrefs) == length(weights)
    vars = [var(mscip, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateConsBasicSOS1(mscip, cons__, "", length(vars), vars, weights)
    @SC SCIPaddCons(mscip, cons__[])
    return store_cons!(mscip, cons__)
end

"""
Add special-ordered-set of type 2 to problem, return cons ref.

# Arguments
- `varrefs::AbstractArray{VarRef}`: variable references
- `weights::AbstractArray{Float64}`: numeric weights
"""
function add_special_ordered_set_type2(mscip::ManagedSCIP, varrefs, weights)
    @assert length(varrefs) == length(weights)
    vars = [var(mscip, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateConsBasicSOS2(mscip, cons__, "", length(vars), vars, weights)
    @SC SCIPaddCons(mscip, cons__[])
    return store_cons!(mscip, cons__)
end

"""
Add abspower constraint to problem, return cons ref.

    lhs ≤ sign(x + a) * abs(x + a)^n + c*z ≤ rhs

# Arguments
- `x::VarRef`: reference for power variable
- `a::Float64`: offset for power variable
- `n::Float64`: exponent for power variable, n >= 1
- `z::VarRef`: reference for linear variable
- `c::Float64`: coefficient for linear variable
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
"""
function add_abspower_constraint(mscip::ManagedSCIP, x, a, n, z, c, lhs, rhs)
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateConsBasicAbspower(
        mscip, cons__, "", var(mscip, x), var(mscip, z), n, a, c, lhs, rhs)
    @SC SCIPaddCons(mscip, cons__[])
    return store_cons!(mscip, cons__)
end

"""
Add indicator constraint to problem, return cons ref.

    y = 1 ==> a^T x ≤ rhs

y has to be a binary variable, or SCIP will error.

# Arguments
- `y::VarRef`: reference for binary indicator variable
- `x::Vector{VarRef}`: reference vector for variables
- `a::Float64`: coefficients for x variable
- `rhs::Float64`: right-hand side for linear constraint
"""
function add_indicator_constraint(mscip::ManagedSCIP, y, x, a, rhs)
    SCIPvarIsBinary(var(mscip, y)) > 0 || error("indicator variable must be binary.")
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    xref = [var(mscip, x[i]) for i in eachindex(x)]
    @SC SCIPcreateConsBasicIndicator(
        mscip, cons__, "", var(mscip, y), length(x), xref, a, rhs
    )
    @SC SCIPaddCons(mscip, cons__[])
    return store_cons!(mscip, cons__)
end

# Transform SCIP C function name as follows:
# 1. Remove leading "SCIP" part (drop the first four characters).
# 2. Convert camel case to snake case.
# For example, `SCIPprintStatusStatistics` becomes `print_status_statistics`.
const STATISTICS_FUNCS = map(x -> Symbol(camel_case_to_snake_case(string(x)[5 : end])), SCIP_STATISTICS_FUNCS)

for (scip_statistics_func, statistics_func) in zip(SCIP_STATISTICS_FUNCS, STATISTICS_FUNCS)
    @eval begin
        """
            $($statistics_func)(mscip::ManagedSCIP)

        Print statistics (calls `$($scip_statistics_func)`) to standard output.
        """
        function $statistics_func end

        function $statistics_func(mscip::ManagedSCIP)
            ret = $scip_statistics_func(mscip, C_NULL)
            ret !== nothing && @assert ret == SCIP_OKAY
            return nothing
        end
    end
end


text:
################################################################################
# Test Suite: Run all test suites
################################################################################

import PyPlot
const Py = PyPlot
using IVPTestSuite

################################################################################
## Select test cases to use
################################################################################
totest = IVPTestSuite.tc_all
# totest = similar(totest)
# totest[:plei] = IVPTestSuite.tc_all[:plei]

################################################################################
## Set Tolerancs and step sizes for tests
################################################################################
# For adaptive solvers
abstols = 10.0.^(-5:-1:-10)
reltols = abstols

# For fixed step solvers
ntsteps = vcat(collect(10.^(1:5)), 500_000)

################################################################################
## Exclude specific solvers
################################################################################
## The best way to exclude specific solvers or include solvers not included
## in release version of a given package, is to comment it out or include it in
## the dictionary of allsolvers (see src/solvers)
################################################################################
## Include suite to be run
################################################################################
include("suite_Sundials.jl")
include("suite_DASSL.jl")
include("suite_ODE_adaptive.jl")
include("suite_ODE_fixedstep.jl")

################################################################################
## Plot results
################################################################################
include("plot_suites.jl")


text:
""
function constraint_power_balance_shunt_dispatch(pm::AbstractActivePowerModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    p    = get(var(pm, n),    :p, Dict()); PowerModels._check_var_keys(p, bus_arcs, "active power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); PowerModels._check_var_keys(pg, bus_gens, "active power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); PowerModels._check_var_keys(ps, bus_storage, "active power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); PowerModels._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); PowerModels._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")

    cstr_p = JuMP.@constraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*1.0)

    if report_duals(pm)
        sol(pm, n, :bus, i)[:lam_kcl_r] = cstr_p
        sol(pm, n, :bus, i)[:lam_kcl_i] = NaN
    end
end

""
function constraint_thermal_limit_from_soft(pm::AbstractActivePowerModel, n::Int, f_idx, rate_a)
    l,i,j = f_idx
    p_fr = var(pm, n, :p, f_idx)
    sm_slack = var(pm, :sm_slack, l)

    JuMP.@constraint(pm.model, p_fr <= rate_a + sm_slack)
end

""
function constraint_thermal_limit_to_soft(pm::AbstractActivePowerModel, n::Int, t_idx, rate_a)
    l,i,j = t_idx
    p_to = var(pm, n, :p, t_idx)
    sm_slack = var(pm, :sm_slack, l)

    JuMP.@constraint(pm.model, p_to <= rate_a + sm_slack)
end

""
function expression_bus_generation(pm::AbstractActivePowerModel, n::Int, i::Int, bus_gens)
    pg = get(var(pm, n), :pg, Dict()); PowerModels._check_var_keys(pg, bus_gens, "active power", "generator")

    pg_total = 0.0
    if length(bus_gens) > 0
        pg_total = sum(pg[g] for g in bus_gens)
    end

    var(pm, n, :bus_pg)[i] = pg_total
end

""
function expression_bus_withdrawal(pm::AbstractActivePowerModel, n::Int, i::Int, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
    ps = get(var(pm, n), :ps, Dict()); PowerModels._check_var_keys(ps, bus_storage, "active power", "storage")

    ps_total = 0.0
    if length(bus_storage) > 0
        ps_total = sum(ps[s] for s in bus_storage)
    end

    pd_total = 0.0
    if length(bus_pd) > 0
        pd_total = sum(pd for pd in values(bus_pd))
    end

    gs_total = 0.0
    if length(bus_gs) > 0
        gs_total = sum(gs for gs in values(bus_gs))*1.0^2
    end

    var(pm, n, :bus_wdp)[i] = ps_total + pd_total + gs_total
end


text:
module roessler
using ReachabilityAnalysis
@taylorize function roessler!(du, u, p, t)
    local a, b, c = 0.2, 0.2, 5.7

    du[1] = -u[2] - u[3]
    du[2] = u[1] + (a * u[2])
    du[3] = b + (u[3]*(u[1] - c))

    return du
end

function model(X0)
    S = @system(x' = roessler!(x), dim:3)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
struct arf_struct
    exponent::UInt      # fmpz
    size::UInt          # mp_size_t
    mantissa1::UInt     # mantissa_struct of length 128
    mantissa2::UInt
end

struct arb_struct
                        # ┌ arf_struct (midpoint)
    exponent::UInt      # │ fmpz
    size::UInt          # │ mp_size_t
    mantissa1::UInt     # │ mantissa_struct of length 128
    mantissa2::UInt     # │
                        # └
                        # ┌ mag_struct (radius)
    exponent_mag::UInt  # │ fmpz
    mantissa_mag::UInt  # │ mp_limb_t
                        # └
end

struct acb_struct
                          # ┌ arb_struct (real)
    exponent_r::UInt      # │ fmpz
    size_r::UInt          # │ mp_size_t
    mantissa1_r::UInt     # │ mantissa_struct of length 128
    mantissa2_r::UInt     # │
    exp_mag_r::Int        # │ fmpz
    mantissa_mag_r::UInt  # │ mp_limb_t
                          # └
                          # ┌ arb_struct (imag)
    exponent_i::Int       # │ fmpz
    size_i::UInt          # │ mp_size_t
    mantissa1_i::UInt     # │ mantissa_struct of length 128
    mantissa2_i::UInt     # │
    exp_mag_i::Int        # │ fmpz
    mantissa_mag_i::UInt  # │ mp_limb_t
                          # └
end

struct mag_struct
    exponent::UInt       # fmpz
    mantissa::UInt       # mp_limb_t
end

for prefix in (:arf, :arb, :acb, :mag)
    arbstruct = Symbol(prefix, :_struct)
    arb_init = Symbol(prefix, :_init)
    arb_clear = Symbol(prefix, :_clear)
    @eval begin
        init!(t::$arbstruct) =
            ccall(@libarb($arb_init), Cvoid, (Ref{$arbstruct},), t)
        clear!(t::$arbstruct) =
            ccall(@libarb($arb_clear), Cvoid, (Ref{$arbstruct},), t)
    end
end