text:
"""
    _diamond_lattice_code(code_matrix::Matrix{TensorNetworkCode})

Takes a matrix of `TensorNetworkCode`s and returns a big `TensorNetworkCode` with all
constituent codes contracted in a diamond lattice.
"""
function _diamond_lattice_code(code_matrix::Matrix{TensorNetworkCode})

    L,W = size(code_matrix)

    # Coordinate shifts for each code in code_matrix
    coord_shifts =
    [(isodd(l) ? [0,0] : [6,0]) + [12*(w-1),6*(l-1)]
        for l in 1:L, w in 1:W]

    # Starting coords for each node (first is virtual & the rest are
    # physical qubits) in a seed code
    coords = [[0,0],[3,-3],[3,3],[0,0.5],[-3,-3],[-3,3]]

    output_code = code_matrix[1,1]
    set_coords!(output_code,coords)

    # contract on the rest of the codes in code_matrix
    for α in 1:L,β in 1:W
        if (α == 1 && β == 1) || (iseven(α) && β == W)
            continue
        end

        next_seed_code = code_matrix[α,β]
        set_coords!(next_seed_code,coords)
        shift = coord_shifts[α,β]
        shift_coords!(next_seed_code,shift)

        output_code = contract_by_coords(output_code,next_seed_code)

        if num_qubits(output_code) == 0
            return output_code
        end
    end

    return output_code
end

_diamond_lattice_code(code_matrix::Matrix{SimpleCode}) =
_diamond_lattice_code(TensorNetworkCode.(code_matrix))

"""
    surface_code_bulk(L::Int64)

Generates the bulk `TensorNetworkCode` for the surface code.
"""
function _surface_code_bulk(L::Int64)

    # Surface code with 5 qubits
    small_surface = five_qubit_surface_code()

    # Other constituent codes are made by choosing certain logical
    # operators to be stabilizers, resulting in stabilizers states
    x_small_surface = gauge(small_surface, 1, 1)
    z_small_surface = gauge(small_surface, 1, 3)

    # this one is rotated:
    z_small_surface_rotated = permute(z_small_surface,[2,5,3,1,4])

    code_matrix = [isodd(i) ? x_small_surface : z_small_surface_rotated
        for i in 1:2*L-1, j in 1:L]

    code_matrix[1,1] = small_surface # this is the only one with a logical
    for i in 2:2*L-1
        if isodd(i)
            code_matrix[i,1] = z_small_surface
        end
    end

    return _diamond_lattice_code(code_matrix)
end

"""
    boundary_fusion(
        L::Int64,
        bulk_code_function=surface_code_bulk)

To create the surface code, we need to get the right boundary
conditions, which this code does.
"""
function _boundary_fusion(L::Int64,bulk_code::TensorNetworkCode)

    # Repitition codes for boundary contraction
    x_repitition_zero = SimpleCode(
        "X repitition code",
        [[3,3,0],[0,3,3],[1,1,1]],
        [],
        [[1,0,0],[0,0,1],[3,0,0]])
    z_repitition_zero = SimpleCode(
        "Z repitition code",
        [[1,1,0],[0,1,1],[3,3,3]],
        [],
        [[3,0,0],[0,0,3],[1,0,0]])
    z_rep = TensorNetworkCode(z_repitition_zero)
    x_rep = TensorNetworkCode(x_repitition_zero)

    # depending on the side (top, left, etc), we choose different coordinates
    left_z_rep = deepcopy(z_rep)
    coords = [[0,0],[-3,0],[-3,-3],[-3,3]]
    set_coords!(left_z_rep,coords)

    right_z_rep = z_rep
    coords = [[0,0],[3,0],[3,-3],[3,3]]
    set_coords!(right_z_rep,coords)

    top_x_rep = deepcopy(x_rep)
    coords = [[0,0],[0,-3],[-3,-3],[3,-3]]
    set_coords!(top_x_rep,coords)

    bottom_x_rep = x_rep
    coords = [[0,0],[0,3],[-3,3],[3,3]]
    set_coords!(bottom_x_rep,coords)

    # Now we can put the code together
    output = bulk_code
    if num_qubits(output) == 0
        return output
    end

    # Let's start fusing boundary codes, starting with the top row
    code = top_x_rep
    for l in 1:L-1
        l==1 ? shift = [6,0] : shift = [12,0]
        shift_coords!(code,shift)
        output = contract_by_coords(output,code)
    end

    code = bottom_x_rep
    for l in 1:L-1
        l == 1 ? shift = [6,12*(L-1)] : shift = [12,0]
        shift_coords!(code, shift)
        output = contract_by_coords(output,code)
    end

    code = left_z_rep
    for l in 1:L-1
        l==1 ? shift = [0,6] : shift = [0,12]
        shift_coords!(code, shift)
        output = contract_by_coords(output,code)
    end

    code = right_z_rep
    for l in 1:L-1
        l==1 ? shift = [12*(L-1),6] : shift = [0,12]
        shift_coords!(code, shift)
        output = contract_by_coords(output,code)
    end

    return output
end

"""
    surface_code(L::Int64)

Returns an (L+1) x L surface code `TensorNetworkCode`.

# Examples
```jldoctest
julia> surf = surface_code(2);

julia> dist = find_distance_logicals(surf)[1] # code distance
3
```
"""
function surface_code(L::Int64)

    bulk_code = _surface_code_bulk(L)
    return _boundary_fusion(L,bulk_code)
end

"""
    almost_rotated_surface_code(
    L::Int64,
    input_seed_code::SimpleCode,
    input_coords::Array{Int64}) -> TensorNetworkCode

Generates a `TensorNetworkCode` that is almost the rotated surface
code.  Requires the `input_seed_code` to have five qubits and to be
in the bulk (`input_coords` with components not equal to 1 or L).

# Examples
```jldoctest
julia> code = almost_rotated_surface_code(3,five_qubit_code(),[2,2]);

julia> verify(code) # try this out and check we get a real code
true
```
"""
function almost_rotated_surface_code(
    L::Int64,
    input_seed_code::SimpleCode,
    input_coords::Array{Int64})

"""
The codes and labels below refer to seed codes laid out in the
following pattern.

3x3:
A  B  C
D  S  F
G  H  I

5x5:
A  B  Br B  C
D  E  Er E  F
Dr Er S  Er Fr
D  E  Er E  F
G  H  Hr H  I

S is a five-qubit surface code with a single logical.
E is a five qubit code with a logical as a stabilizer.
Er is the same as E with Hadamards applied to every qubit.
The rest are detailed below.
"""

if iseven(L)
    error("This only works for odd L!")
end

code_dictionary = Dict{String,TensorNetworkCode}()

# A
stabilizers = [[1,1,1],[0,3,3],[3,0,3]]
logicals = []
code = SimpleCode("A",stabilizers,logicals)
A = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[1,0],[0.3,0.3]]
set_coords!(A,coordinates)
code_dictionary["A"] = A

# B
stabilizers = [[1,1,0,1],[3,0,0,3],[3,0,3,3],[0,3,0,3]]
logicals = []
code = SimpleCode("B",stabilizers,logicals)
B = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[1,0],[0.3,0.3]]
set_coords!(B,coordinates)
code_dictionary["B"] = B

# Br
stabilizers = [[3,3,0,3],[0,1,1,1],[1,0,1,1],[0,0,3,3]]
logicals = []
code = SimpleCode("Br",stabilizers,logicals)
Br = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[1,0],[0.3,0.3]]
set_coords!(Br,coordinates)
code_dictionary["Br"] = Br

# C
stabilizers = [[3,3,3],[0,1,1],[1,0,1]]
logicals = []
code = SimpleCode("C",stabilizers,logicals)
C = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[0.3,0.3]]
set_coords!(C,coordinates)
code_dictionary["C"] = C

# D
stabilizers = [[0,0,1,1],[0,1,1,1],[3,0,3,3],[1,0,0,1]]
logicals = []
code = SimpleCode("D",stabilizers,logicals)
D = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[0,-1],[1,0],[0.3,0.3]]
set_coords!(D,coordinates)
code_dictionary["D"] = D

# Dr
stabilizers = [[3,3,0,3],[0,3,3,3],[1,0,1,1],[0,1,0,1]]
logicals = []
code = SimpleCode("Dr",stabilizers,logicals)
Dr = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[0,-1],[1,0],[0.3,0.3]]
set_coords!(Dr,coordinates)
code_dictionary["Dr"] = Dr

# E
stabilizers = [[0,1,1,0,1],[1,0,0,1,1],[3,3,0,0,3],[0,0,3,3,3],[1,1,0,0,0]]
logicals = []
code = SimpleCode("E",stabilizers,logicals)
E = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[0,-1],[1,0],[0.3,0.3]]
set_coords!(E,coordinates)
code_dictionary["E"] = E

# Er
stabilizers = [[0,3,3,0,3],[3,0,0,3,3],[1,1,0,0,1],[0,0,1,1,1],[3,3,0,0,0]]
logicals = []
code = SimpleCode("Er",stabilizers,logicals)
Er = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[0,-1],[1,0],[0.3,0.3]]
set_coords!(Er,coordinates)
code_dictionary["Er"] = Er

# F
stabilizers = [[0,3,3,3],[1,1,0,1],[0,1,0,1],[0,0,1,1]]
logicals = []
code = SimpleCode("F",stabilizers,logicals)
F = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[0,-1],[0.3,0.3]]
set_coords!(F,coordinates)
code_dictionary["F"] = F

# Fr
stabilizers = [[0,1,1,1],[3,3,0,3],[3,0,3,3],[1,0,0,1]]
logicals = []
code = SimpleCode("Fr",stabilizers,logicals)
Fr = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[0,-1],[0.3,0.3]]
set_coords!(Fr,coordinates)
code_dictionary["Fr"] = Fr

# G
stabilizers = [[3,3,3],[1,0,1],[0,1,1]]
logicals = []
code = SimpleCode("G",stabilizers,logicals)
G = TensorNetworkCode(code)
coordinates = [[0,0],[0,-1],[1,0],[0.3,0.3]]
set_coords!(G,coordinates)
code_dictionary["G"] = G

# H
stabilizers = [[3,3,0,3],[0,1,1,1],[0,3,0,3],[0,0,3,3]]
logicals = []
code = SimpleCode("H",stabilizers,logicals)
H = TensorNetworkCode(code)
coordinates = [[0,0],[-1,0],[0,-1],[1,0],[0.3,0.3]]
set_coords!(H,coordinates)
code_dictionary["H"] = H

# Hr
stabilizers = [[1,1,0,1],[0,3,3,3],[1,0,1,1],[3,0,0,3]]
logicals = []
code = SimpleCode("Hr",stabilizers,logicals)
Hr = TensorNetworkCode(code)
coordinates = [[0,0],[-1,0],[0,-1],[1,0],[0.3,0.3]]
set_coords!(Hr,coordinates)
code_dictionary["Hr"] = Hr

# I
stabilizers = [[1,1,1],[3,0,3],[0,3,3]]
logicals = []
code = SimpleCode("I",stabilizers,logicals)
I = TensorNetworkCode(code)
coordinates = [[0,0],[-1,0],[0,-1],[0.3,0.3]]
set_coords!(I,coordinates)
code_dictionary["I"] = I

# Assign code tensor labels to a matrix
code_label_matrix = fill("", L, L)
for i in 1:L, j in 1:L
    if iseven(i+j)
        code_label_matrix[i,j] = "E"
    else
        code_label_matrix[i,j] = "Er"
    end
end

for j in 1:L
    if iseven(j)
        code_label_matrix[1,j] = "B"
    else
        code_label_matrix[1,j] = "Br"
    end
end

for j in 1:L
    if iseven(j)
        code_label_matrix[L,j] = "H"
    else
        code_label_matrix[L,j] = "Hr"
    end
end

for i in 1:L
    if iseven(i)
        code_label_matrix[i,1] = "D"
    else
        code_label_matrix[i,1] = "Dr"
    end
end

for i in 1:L
    if iseven(i)
        code_label_matrix[i,L] = "F"
    else
        code_label_matrix[i,L] = "Fr"
    end
end

# Take care of the corners
code_label_matrix[1,1] = "A"
code_label_matrix[1,L] = "C"
code_label_matrix[L,1] = "G"
code_label_matrix[L,L] = "I"

# Centre code
stabilizers = [[0,1,1,0,1],[1,0,0,1,1],[3,3,0,0,3],[0,0,3,3,3]]
logicals = [[3,0,3,0,3],[0,1,0,1,1]]
code = SimpleCode("central",stabilizers,logicals)
central = TensorNetworkCode(code)
coordinates = [[0,0],[0,1],[-1,0],[0,-1],[1,0],[0.3,0.3]]
set_coords!(central,coordinates)
code_dictionary["central"] = central
# Add to the code label matrix
centre = ceil(Int,L/2)
code_label_matrix[centre,centre] = "central"

# Include input code
if num_qubits(input_seed_code) != 5
    error("input_seed_code does not have five physical qubits!")
elseif 1 in input_coords || L in input_coords
    error("location of input_seed_code needs to be in the bulk!")
end

input_tn_code = TensorNetworkCode(input_seed_code)
coordinates = [[0,0],[0,1],[-1,0],[0,-1],[1,0],[0.3,0.3]]
set_coords!(input_tn_code,coordinates)
code_dictionary["input"] = input_tn_code

# Add to the code label matrix
code_label_matrix[input_coords...] = "input"

output_code = deepcopy(code_dictionary[code_label_matrix[1,1]])
for i in 1:L
    for j in 1:L
        if i==j && i==1
            continue
        end
        code = deepcopy(code_dictionary[code_label_matrix[i,j]])
        shift = [Float64(j-1)*2,Float64(i-1)*2]
        shift_coords!(code,shift)
        output_code = contract_by_coords(output_code,code)
    end
end

return output_code
end

"""
    rotated_surface_code(L::Int64) -> TensorNetworkCode

Generates the LxL rotated surface code as a `TensorNetworkCode`.

# Examples
```jldoctest
julia> rot = rotated_surface_code(3);

julia> pauli_weight(rot.stabilizers[2]) # weight of a plaquette stabilizer
4
```
"""
function rotated_surface_code(L::Int64)

stabilizers = [[0,1,1,0,1],[1,0,0,1,1],[3,3,0,0,3],[0,0,3,3,3]]
logicals = [[3,0,3,0,3],[0,1,0,1,1]]
code = SimpleCode("central",stabilizers,logicals)

# Add to the code label matrix
centre = ceil(Int,L/2)
input_coords = [centre,centre]

return almost_rotated_surface_code(L,code,input_coords)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule lrslib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("lrslib")
JLLWrappers.@generate_main_file("lrslib", UUID("3873f7d0-7b7c-52c3-bdf4-8ab39b8f337a"))
end  # module lrslib_jll

ext: jl


text:
function state_space_validation(A,B,C,D)
    nx = size(A, 1)
    nu = size(B, 2)
    ny = C == I ? nx : size(C, 1)
    colC = C == I ? nx : size(C, 2)

    if size(A, 2) != nx && nx != 0
        error("A has dimentions $(size(A)), but must be square")
    elseif size(B, 1) != nx
        error("The number of rows of A ($(size(A,1))) and B ($(size(B,1))) are not equal")
    elseif colC != nx
        error("The number of columns of A ($(size(A,2))) and C ($(colC)) are not equal")
    elseif nu != size(D, 2)
        error("The number of columns of B ($(size(B,2))) and D ($(size(D,2))) are not equal")
    elseif ny != size(D, 1)
        error("The number of rows of C ($(ny)) and D ($(size(D,1))) are not equal")
    end

    nx,nu,ny
end

mutable struct LdtiSubsystem{T<:Real} <:LinearSubsystem{T}
    A::AbstractMatrix{T}
    B::AbstractMatrix{T}
    C::OutputMatrix{T}
    D::AbstractMatrix{T}
    E::AbstractMatrix{T}
    index::UInt
    neighbours::Array{LdtiSubsystem{T}}
    Aij::Array{AbstractMatrix{T}}
    x::Vector{T}
    xnext::Vector{T}
    function LdtiSubsystem(i, A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::OutputMatrix{T}, D::AbstractMatrix{T}) where {T<:Real}
        (nx, nu, ny) = state_space_validation(A, B, C, D)
        new{T}(A, B, C, D, zeros(T, nx, 0), i, [], [], zeros(T, nx), zeros(T, nx))
    end
    function LdtiSubsystem(i, A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::OutputMatrix{T}) where {T<:Real}
        ny = C == I ? size(A,1) : size(C,1)
        D = zeros(T, ny, size(B, 2))
        LdtiSubsystem(i, A, B, C, D)
    end
end

function Base.zero(s::LdtiSubsystem{T}) where {T}
    LdtiSubsystem(0, zeros(T, s.nx,s.nx), zeros(T, s.nx, 1), I)
end

function setInitialState(s::LdtiSubsystem{T}, x::Vector{T}) where {T}
    setfield!(s, :x, x)
    setfield!(s, :xnext, x)
end

function updateState(s::LdtiSubsystem{T}, u::Union{T, Vector{T}}, w::Vector{T}) where {T}
    if isempty(s.neighbours)
        setfield!(s, :xnext, s.A * s.x + (s.B * u)[:] .+ w)
    else
        barx = getNeighbourStates(s)
        setfield!(s, :xnext, s.A * s.x + s.B * u + s.E * barx + w)
    end
    s.xnext
end
updateState(s::LdtiSubsystem{T}, u::Union{T, Vector{T}}) where {T} = updateState(s, u, zeros(T, s.nx))

function outputMap(s::LinearSubsystem{T}, u::Union{T, Vector{T}}, v::Vector{T}) where {T}
    setfield!(s, :x, s.xnext)
    s.C * s.x + (s.D * u)[:] + v
end
outputMap(s::LinearSubsystem{T}, u::Union{T, Vector{T}}) where {T} = outputMap(s, u, zeros(T, s.ny))

getNeighbourStates(s::AbstractSubsystem) = isempty(s.neighbours) ? zero(s.x) : foldl(vcat, [s.x for s in s.neighbours])

function addNeighbour(s::LinearSubsystem{T}, neighbour::LinearSubsystem{T}, weight::AbstractMatrix{T}, conservation::Bool = false) where {T}
    if neighbour.index == s.index
        error("Self loops are not allowed")
    end
    push!(s.neighbours, neighbour)
    push!(s.Aij, weight)
    updateE(s)
    if conservation
       s.A -= weight
    end
    map(x-> convert(Int, x), getNeighbourIndices(s))
end

function removeNeighbour(s::AbstractSubsystem, neighbour::Integer, conservation::Bool = false)
    idx = findfirst(x -> x.index == neighbour, s.neighbours)
    deleteat!(s.neighbours, idx)
    weight = splice!(s.Aij, idx)
    if conservation
        s.A += weight
    end
    updateE(s)
end

function setNeighbour(s::LinearSubsystem, neighbour::LinearSubsystem, neighbourIdx::Integer, weight::AbstractMatrix{T}, conservation::Bool = false) where {T}
    j = findfirst(x -> x.index == neighbourIdx, s.neighbours)
    s.neighbours[j] = neighbour
    if conservation
        s.A = s.A + s.Aij[j] - weight
    end
    s.Aij[j] = weight
    updateE(s)
end

function getNeighbourIndices(s::AbstractSubsystem)
    Tuple([s.index for s in s.neighbours])
end

function updateE(s::LinearSubsystem{T}) where {T}
    setfield!(s, :E, isempty(s.Aij) ? zeros(T, s.nx, 0) : foldl(hcat, s.Aij))
end

function getOutboundWeights(s::LinearSubsystem{T}) where {T}
    outWeights = Vector{AbstractMatrix{T}}()
    for sj in s.neighbours
        ji = findfirst(x -> s.index == x, getNeighbourIndices(sj))
        push!(outWeights, sj.Aij[ji])
    end
    outWeights
end

Base.eltype(::Type{S}) where {S<:LinearSubsystem} = S

function Base.getproperty(sys::LinearSubsystem{T}, s::Symbol) where {T}
    if s === :nx
        return size(sys.A, 1)
    elseif s === :nu
        return size(sys.B, 2)
    elseif s === :ny
        return size(sys.D, 1)
    elseif s === :nx̄
        return size(sys.E,2)
    elseif s === :Cmat
        return sys.C == I ? AbstractMatrix{T}(I, sys.ny, sys.ny) : sys.C
    else
        return getfield(sys, s)
    end
end

function Base.setproperty!(sys::LinearSubsystem, s::Symbol, val)
    if s ∉ [:A, :B, :C, :D, :index]
        error("Protected property")
    else
        if size(val) != size(getfield(sys, s))
            error("Cannot change size")
        else
            return setfield!(sys, s, val)
        end
    end
end

function Base.show(io::IO, s::LinearSubsystem)
    println(io, "Subsystem $(s.index):")
    if size(s.A, 1) > 0
        println(io, "A = \n", _string_mat_with_headers(s.A))
        println(io, "B = \n", _string_mat_with_headers(s.B))
        println(io, "C = ", s.C == I ? "I" : "\n"*_string_mat_with_headers(s.C))
    end
    println(io, "D = \n", _string_mat_with_headers(s.D), "\n")
    if !isempty(s.neighbours)
        println(io, "Neighbours = $(convert.(Int64, getNeighbourIndices(s)))")
        println(io, "E = \n", _string_mat_with_headers(s.E))
    end
end


text:
# Tests that take more time to run

## FIXME: bitrot, returns unsigned int, must be converted to int
## test passes, but we should not have to to conversion
@test convert(Int,ntcountprimes("10^19", "10^19+10^6")) == 23069

# fixed overflow bugs
@test countprimes(Int128(10)^19, Int128(10)^19+10^3) == 28
@test countprimes(Int128(10)^19+10^9) == Int128(234057667299198865)
## FIXME: bitrot, sign error
##@test countprimes("10^19 + 10^9") == 234057667299198865


text:
# # Example 6 : Climatological Mean
#
# An example of one tile with a climatology time axis from the global ocean domain available in https://github.com/gaelforget/nctiles-testcases
#

using NCTiles,NCDatasets

# File Paths
inputs=NCTiles.NCTILES_TESTCASES
NCTiles.ensure_testcases_installed()

outputs = joinpath(tempdir(),"NCTILES_TESTCASES_OUTPUT/")
if ~ispath(outputs); mkpath(outputs); end

savedir = joinpath(outputs,"ex6")
if ~ispath(savedir); mkpath(savedir); end

field_name = "FeT"
README = [field_name*" -- Source: Gael Forget; version: alpha."];

# ### One tile example

ncvars,ncdims,fileatts = readncfile(joinpath(inputs,"diags_nctiles/FeT.0062.nc"))
rm(joinpath(savedir,"ex6a.nc"),force=true)
write(ncvars,joinpath(savedir,"ex6a.nc"),README=README)

# ### Define Dimensions
#
# For reference:
#
#
# ```
# - i_c: NCvar("i_c", "1", 30, 1.0:30.0, Dict("units" => "1","long_name" => "Cartesian coordinate 1"), NCDatasets)
# - j_c: NCvar("j_c", "1", 30, 1.0:30.0, Dict("units" => "1","long_name" => "Cartesian coordinate 2"), NCDatasets)
# - k_c: NCvar("k_c", "1", 50, 1.0:50.0, Dict("units" => "1","long_name" => "Cartesian coordinate 3"), NCDatasets)
# - t:   NCvar("t", "1", 12, 1.0:12.0, Dict("units" => "1","long_name" => "Time coordinate"), NCDatasets)
# - tcb:  NCvar("tcb", "", 2, Any[], Dict{Any,Any}(), NCDatasets)
# ```
#
#
# Note : tcb is unitles, only has a dimension, no values

FeT_dims = [ncdims["i_c"],
                ncdims["j_c"],
                ncdims["k_c"],
                ncdims["t"]]
clim_dims = [ncdims["tcb"],
                ncdims["t"]]

# ### Define Variables
#
# Note : lat, lon, dep, tim, thic, area, land defined same as non-climatology example

FeT = NCvar("FeT", # name
                "mmol Fe", # units
                FeT_dims, # dimensions
                NCData(joinpath(inputs,"diags_nctiles/FeT.0062.nc"), "FeT", NCDatasets, Float32), # values- to be read from file
                Dict("coordinates" => "lon lat dep tim","long_name" => "FeT concentration"), # attributes
                NCDatasets) # backend

climatology_bounds = NCvar("climatology_bounds", # name
                                "days since 1992-1-1 0:0:0", # units
                                clim_dims, # dimensions
                                NCData(joinpath(inputs,"diags_nctiles/FeT.0062.nc"), "climatology_bounds", NCDatasets, Float32), # values- to be read from file
                                Dict("long_name" => "climatology_bounds"), # attributes
                                ncvars["climatology_bounds"].backend) # backend

# Write to file

writevars = Dict(["FeT" => FeT,
                "lon" => ncvars["lon"],
                "lat" => ncvars["lat"],
                "dep" => ncvars["dep"],
                "tim" => ncvars["tim"],
                "thic" => ncvars["thic"],
                "area" => ncvars["area"],
                "land" => ncvars["land"],
                "climatology_bounds" => climatology_bounds])
rm(joinpath(savedir,"ex6b.nc"),force=true)
write(writevars,joinpath(savedir,"ex6b.nc"),README=README)


text:
mas2rad(mas) = mas * pi / 648000000

function gauss(uv, inc, pa, fwhm, fluxg, fluxp)
    bl = baseline(uv, inc, pa)
    Θ = mas2rad(fwhm)
    return fluxg * exp.(-(π * Θ * bl).^2 / (4 * log(2))) .+ fluxp
end

function disk(uv, inc, pa, diam, fluxd, fluxp)
    bl = baseline(uv, inc, pa)
    Θ = mas2rad(diam)
    output = fill(1.0, size(bl))
    blmask = .!iszero.(bl)
    output[blmask] = 2 * besselj1.(π * Θ * bl[blmask]) ./ (π * Θ * bl[blmask])
    return fluxd * output .+ fluxp
end

function annulus(uv, inc, pa, diam, fluxa, fluxp)
    bl = baseline(uv, inc, pa)
    rho_0 = mas2rad(diam) / 2
    return fluxa * besselj0.(2 * π * rho_0 * bl) .+ fluxp
end

function pseudolorentzian(uv, inc, pa, hfdiam, fluxl, fluxp)
    bl = baseline(uv, inc, pa)
    theta = mas2rad(hfdiam)
    return fluxl * exp.(-π * theta * bl / sqrt(3)) .+ fluxp
end


text:
using JuMP, EAGO, Gurobi

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[5], 0.0)
JuMP.set_lower_bound(x[4], 0.0)
JuMP.set_lower_bound(x[2], 0.0)
JuMP.set_lower_bound(x[3], 0.0)
JuMP.set_lower_bound(x[1], 10.0)
JuMP.set_upper_bound(x[1], 2000.0)
JuMP.set_upper_bound(x[2], 16000.0)
JuMP.set_upper_bound(x[3], 120.0)
JuMP.set_upper_bound(x[4], 5000.0)
JuMP.set_upper_bound(x[5], 2000.0)
JuMP.set_lower_bound(x[6], 85.0)
JuMP.set_upper_bound(x[6], 93.0)
JuMP.set_lower_bound(x[7], 90.0)
JuMP.set_upper_bound(x[7], 95.0)
JuMP.set_lower_bound(x[8], 3.0)
JuMP.set_upper_bound(x[8], 12.0)
JuMP.set_lower_bound(x[9], 1.2)
JuMP.set_upper_bound(x[9], 4.0)
JuMP.set_lower_bound(x[10], 145.0)
JuMP.set_upper_bound(x[10], 162.0)

# ----- Constraints ----- #
@NLconstraint(m, e1, -x[1]*(1.12+0.13167*x[8]-0.00667* (x[8])^2)+x[4] == 0.0)
@constraint(m, e2, -x[1]+1.22*x[4]-x[5] == 0.0)
@NLconstraint(m, e3, -0.001*x[4]*x[9]*x[6]/(98-x[6])+x[3] == 0.0)
@NLconstraint(m, e4, -(1.098*x[8]-0.038* (x[8])^2)-0.325*x[6]+x[7] == 57.425)
@NLconstraint(m, e5, -(x[2]+x[5])/x[1]+x[8] == 0.0)
@constraint(m, e6, x[9]+0.222*x[10] == 35.82)
@constraint(m, e7, -3*x[7]+x[10] == -133.0)

# ----- Objective ----- #
@NLobjective(m, Max, 0.063*x[4]*x[7] - 5.04*x[1] - 0.035*x[2] - 10*x[3] - 3.36*x[5]) # CAN BE QUADRATIC

JuMP.optimize!(m)

run_time = backend(m).optimizer.model.optimizer._run_time
#objective_value =
#objective_bound =
#println("objective between: ")
println("run time: $run_time")


text:
push!(LOAD_PATH, "../src/")
using Documenter, Neuroimaging

makedocs(
    modules = [Neuroimaging],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Neuroimaging.jl",
    authors = "Robert Luke",
    pages = [
        "Home" => "index.md",
        "Types" => "types.md",
        "Usage/Implementation Details" =>
            Any["Accessing data"=>"usage/access.md", "Filtering"=>"usage/filter.md"],
        "General EEG Processing" => Any[
            "Overview"=>"eeg/eeg.md",
            "Example"=>"eeg/examples.md",
            "API"=>"eeg/functions.md",
        ],
        "Steady State Responses (EEG)" => Any[
            "Overview"=>"assr/assr.md",
            "Example"=>"assr/examples.md",
            "API"=>"assr/functions.md",
        ],
        "Transient Responses (EEG)" =>
            Any["Overview"=>"tr/transient.md", "Example"=>"tr/examples.md"],
        "Coordinate Systems" => Any["Example"=>"coord/examples.md",],
        "Source Modelling" => Any["Example"=>"source/examples.md",],
        "Input/Output Support" => "IO.md",
        "Low-Level API" => "api.md",
    ],
)

deploydocs(
    repo = "github.com/rob-luke/Neuroimaging.jl.git",
    push_preview = true,
    devbranch = "main",
)


text:
using LinearAlgebra: diagind

"""
    OuyangModelOne(A=0.01, B=1.6, Γ=1e-4, N=10, ΔE=1.6e-2, D=1.0)

Model #1 from [Ouyang and Subotnik](http://aip.scitation.org/doi/10.1063/1.4908032).
See also [Ouyang's thesis](https://repository.upenn.edu/edissertations/2508/).
"""
Parameters.@with_kw struct OuyangModelOne <: DiabaticFrictionModel
    A::Float64 = 0.01
    B::Float64 = 1.6
    Γ::Float64 = 1e-4
    N::Int = 10
    ΔE::Float64 = 1.6e-2
    ρ::Float64 = N / ΔE
    C::Float64 = sqrt(Γ / 2π*ρ) / N
    D::Float64 = 1.0

    n_states::Int = N + 1
end

NQCModels.nstates(model::OuyangModelOne) = model.n_states
NQCModels.ndofs(::OuyangModelOne) = 1

function NQCModels.potential!(model::OuyangModelOne, V::Hermitian, r::Real)
    Parameters.@unpack A, B, C, D, ΔE, N = model

    Vsys(x) = A*tanh(B*x)
    Vbath(x, Δ) = -Vsys(x) + Δ
    Vsb(x) = C*exp(-D*x^2)

    V[1,1] = Vsys(r)
    Δs = range(-ΔE/2, ΔE/2, length=N)
    V.data[1,2:N+1] .= Vsb(r)
    V[diagind(N+1,N+1)[2:end]] .= Vbath.(r, Δs)

    return V
end

function NQCModels.derivative!(model::OuyangModelOne, derivative::Hermitian, r::Real)
    Parameters.@unpack A, B, C, D, N = model

    Dsys(x) = B*A*sech(B*x)^2
    Dbath(x) = -Dsys(x)
    Dsb(x) = -2D*x*C*exp(-D*x^2)

    derivative[1,1] = Dsys(r)
    derivative.data[1,2:N+1] .= Dsb(r)
    derivative[diagind(N+1,N+1)[2:end]] .= Dbath(r)

    return derivative
end


text:
module TestSpheroidal

export test_m_symmetry, test_s_symmetry, test_conj_symmetry, compare_to_qnm

include("../src/CustomTypes.jl")
include("../src/Spheroidal.jl")
include("ReadQNM.jl")
using .CustomTypes
using .ReadQNM
import .Spheroidal as SH

import Test: @test

const tolerance = 1e-8 ## tolerance we compare to

##============================================================
"""
Test la_{s,l,m,n}(-c) = la_{s,l,-m,n}(c)
"""
function test_m_symmetry(
      nl::myI,
      neig::myI,
      s::myI,
      m::myI,
      a::myF,
      om::myC
   )
   la_pm, ph_pm = SH.eig_vals_vecs(nl, neig, s,  m, -a*om)
   la_nm, ph_nm = SH.eig_vals_vecs(nl, neig, s, -m,  a*om)

   for (i,_) in enumerate(la_pm)
      @test abs(la_pm[i] - la_nm[i]) < tolerance
   end

end
"""
Test la_{-s,l,m,n}(c) = 2s + la_{s,l,m,n}(c)
"""
function test_s_symmetry(
      nl::myI,
      neig::myI,
      s::myI,
      m::myI,
      a::myF,
      om::myC
   )
   la_ps, ph_ps = SH.eig_vals_vecs(nl, neig,  s, m, a*om)
   la_ns, ph_ns = SH.eig_vals_vecs(nl, neig, -s, m, a*om)

   for (i,_) in enumerate(la_ps)
      @test abs(la_ps[i] + 2.0*s - la_ns[i]) < tolerance
   end
end
"""
Test (la_{s,l,m,n}(c))* = 2s + la_{s,l,m,n}(c*)
"""
function test_conj_symmetry(
      nl::myI,
      neig::myI,
      s::myI,
      m::myI,
      a::myF,
      om::myC
   )
   la,   ph   = SH.eig_vals_vecs(nl, neig, s, m, a*     om )
   la_c, ph_c = SH.eig_vals_vecs(nl, neig, s, m, a*conj(om))

   for (i,_) in enumerate(la)
      @test abs(conj(la[i]) - la_c[i]) < tolerance
   end
end

"""
Compare against values computed by Leo Stein's qnm code,

J.Open Source Softw. 4 (2019) 42, 1683

This is written in python so we need to bind to the library
"""
function compare_to_qnm(
      nl::myI,
      neig::myI,
      s::myI,
      n::myI,
      l::myI,
      m::myI,
      avals::Vector{myF}
   )

   lmin = SH.compute_l_min(s,m)

   println("Comparison to qnm library: s=$s, n=$n, l=$l, m=$m")

   for a in avals
      println("testing a=$a")
      om, la = qnm(n,s,m,l,a)
      ls, vs = SH.eig_vals_vecs(nl, neig, s, m, a*om)

      @test abs(ls[l-lmin+1]-la)<tolerance
   end
end

end


text:
"""
    get_fhir_access_token(smart_result) -> AbstractString
"""
function get_fhir_access_token end

"""
    has_fhir_patient_id(smart_result) -> Bool
"""
function has_fhir_patient_id end

has_fhir_patient_id(smart_result) = false

"""
    get_fhir_patient_id(smart_result) -> AbstractString
"""
function get_fhir_patient_id end

"""
    has_fhir_encounter_id(smart_result) -> Bool
"""
function has_fhir_encounter_id end

has_fhir_encounter_id(smart_result) = false

"""
    get_fhir_encounter_id(smart_result) -> AbstractString
"""
function get_fhir_encounter_id end


text:
"""
    sum(N::AbstractEcologicalNetwork; dims=nothing)

This function will return the sum of all interactions in the network. For
quantitative networks, this is the sum of interaction strengths. For binary
networks, this is the number of interactions. For probabilistic networks, this
is the expected number of realized interactions.

Optionally, one can give the argument dims, simular to the native `sum`, which
computes the sum of the interactions for the lower (`dims=2`) or higher
(`dims=1`) level.

#### References

- Delmas, E., Besson, M., Brice, M.-H., Burkle, L.A., Dalla Riva, G.V., Fortin,
  M.-J., Gravel, D., Guimarães, P.R., Hembry, D.H., Newman, E.A., Olesen, J.M.,
  Pires, M.M., Yeakel, J.D., Poisot, T., 2018. Analysing ecological networks of
  species interactions. Biological Reviews 112540.
  https://doi.org/10.1111/brv.12433

- Dunne, J.A., 2006. The Network Structure of Food Webs, in: Dunne, J.A.,
  Pascual, M. (Eds.), Ecological Networks: Linking Structure and Dynamics.
  Oxford University Press, pp. 27–86.

- Martinez, N.D., 1992. Constant Connectance in Community Food Webs. The
  American Naturalist 139, 1208–1218.
"""
function Base.sum(N::AbstractEcologicalNetwork; dims::Union{Nothing,Int}=nothing)
    @assert isnothing(dims) || dims == 1 || dims == 2
    isnothing(dims) && return sum(N.edges)
    return sum(N.edges, dims=dims)
end

"""
    links(N::BinaryNetwork)

Number of non-zero interactions in a deterministic network.

#### References

- Delmas, E., Besson, M., Brice, M.-H., Burkle, L.A., Dalla Riva, G.V., Fortin,
  M.-J., Gravel, D., Guimarães, P.R., Hembry, D.H., Newman, E.A., Olesen, J.M.,
  Pires, M.M., Yeakel, J.D., Poisot, T., 2018. Analysing ecological networks of
  species interactions. Biological Reviews 112540.
  https://doi.org/10.1111/brv.12433

- Dunne, J.A., 2006. The Network Structure of Food Webs, in: Dunne, J.A.,
  Pascual, M. (Eds.), Ecological Networks: Linking Structure and Dynamics.
  Oxford University Press, pp. 27–86.

- Martinez, N.D., 1992. Constant Connectance in Community Food Webs. The
  American Naturalist 139, 1208–1218.
"""
function links(N::BinaryNetwork)
    return count(!iszero, N.edges)
end

"""
    links(N::QuantitativeNetwork)

Number of non-zero interactions in a quantitative network (use `sum` to get the
sum of interaction strengths).

#### References

- Delmas, E., Besson, M., Brice, M.-H., Burkle, L.A., Dalla Riva, G.V., Fortin,
  M.-J., Gravel, D., Guimarães, P.R., Hembry, D.H., Newman, E.A., Olesen, J.M.,
  Pires, M.M., Yeakel, J.D., Poisot, T., 2018. Analysing ecological networks of
  species interactions. Biological Reviews 112540.
  https://doi.org/10.1111/brv.12433

- Dunne, J.A., 2006. The Network Structure of Food Webs, in: Dunne, J.A.,
  Pascual, M. (Eds.), Ecological Networks: Linking Structure and Dynamics.
  Oxford University Press, pp. 27–86.

- Martinez, N.D., 1992. Constant Connectance in Community Food Webs. The
  American Naturalist 139, 1208–1218.
"""
function links(N::QuantitativeNetwork)
    return count(!iszero, N.edges)
end

"""
    links(N::ProbabilisticNetwork)

Expected number of interactions in a probabilistic network. To get the number of
interactions that have a non-zero probability, use *e.g.* `links(N>0.0)`.

#### References

- Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
  D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
  Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function links(N::ProbabilisticNetwork)
    return sum(N)
end

"""
    links_var(N::ProbabilisticNetwork)

Expected variance of the number of links for a probabilistic network.

#### References

- Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
  D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
  Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function links_var(N::ProbabilisticNetwork)
   return sum(N.edges .* (1 .- N.edges))
end

"""
    connectance(N::AbstractEcologicalNetwork)

Number of links divided by the number of possible interactions. In unipartite
networks, this is ``L/S^2``. In bipartite networks, this is ``L/(T × B)``. It is
worth noting that while the maximal connectance is always 1 (i.e. the graph is
complete), the minimum value (assuming that the network is not degenerate) is
*not* 0. Instead, the minimum number of interactions in a unipartite network is
`S-1`, and in a bipartite network it is `max(T,B)`.

Connectance can therefore be transformed between 0 and 1, using the following
approach: let `m` be the minimum number of interactions, and Co be the measured
connectance, then the corrected value is `(Co-m)/(1-m)`. To our best knowledge,
this is not standard practice, and therefore is not suggested as a function in
the package.

#### References

- Delmas, E., Besson, M., Brice, M.-H., Burkle, L.A., Dalla Riva, G.V., Fortin,
  M.-J., Gravel, D., Guimarães, P.R., Hembry, D.H., Newman, E.A., Olesen, J.M.,
  Pires, M.M., Yeakel, J.D., Poisot, T., 2018. Analysing ecological networks of
  species interactions. Biological Reviews 112540.
  https://doi.org/10.1111/brv.12433

- Dunne, J.A., 2006. The Network Structure of Food Webs, in: Dunne, J.A.,
  Pascual, M. (Eds.), Ecological Networks: Linking Structure and Dynamics.
  Oxford University Press, pp. 27–86.

- Martinez, N.D., 1992. Constant Connectance in Community Food Webs. The
  American Naturalist 139, 1208–1218.
"""
function connectance(N::AbstractEcologicalNetwork)
    return links(N) / (richness(N; dims=1)*richness(N; dims=2))
end

"""
    linkage_density(N::AbstractEcologicalNetwork)

Number of links divided by species richness.

#### References

- Delmas, E., Besson, M., Brice, M.-H., Burkle, L.A., Dalla Riva, G.V., Fortin,
  M.-J., Gravel, D., Guimarães, P.R., Hembry, D.H., Newman, E.A., Olesen, J.M.,
  Pires, M.M., Yeakel, J.D., Poisot, T., 2018. Analysing ecological networks of
  species interactions. Biological Reviews 112540.
  https://doi.org/10.1111/brv.12433

- Dunne, J.A., 2006. The Network Structure of Food Webs, in: Dunne, J.A.,
  Pascual, M. (Eds.), Ecological Networks: Linking Structure and Dynamics.
  Oxford University Press, pp. 27–86.
"""
function linkage_density(N::AbstractEcologicalNetwork)
    return links(N) / richness(N)
end

"""
    connectance_var(N::ProbabilisticNetwork)

Expected variance of the connectance for a probabilistic matrix, measured as the
variance of the number of links divided by the squared size of the matrix.

#### References

- Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
  D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
  Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function connectance_var(N::ProbabilisticNetwork)
   return links_var(N) / (prod(size(N))^2)
end


text:
using Plots, Printf, LinearAlgebra

# enable plotting by default
if !@isdefined do_visu; do_visu = true end

# finite-difference support functions
@views av_xi(A) = 0.5*(A[1:end-1,2:end-1].+A[2:end,2:end-1])
@views av_yi(A) = 0.5*(A[2:end-1,1:end-1].+A[2:end-1,2:end])
@views   inn(A) = A[2:end-1,2:end-1]

@views function diffusion_2D_expl(; do_visu=true, save_fig=false)
    # Physics
    lx, ly = 10.0, 10.0   # domain size
    npow   = 3            # power-law exponent
    ttot   = 1.0          # total simulation time
    # Numerics
    nx, ny = 128, 128     # numerical grid resolution
    # Derived numerics
    dx, dy = lx/nx, ly/ny # grid size
    xc, yc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    qHx    = zeros(nx-1, ny-2) # on staggered grid
    qHy    = zeros(nx-2, ny-1) # on staggered grid
    dHdt   = zeros(nx-2, ny-2) # normal grid, without boundary points
    # Initial condition
    H      = exp.(.-(xc.-lx/2).^2 .-(yc.-ly/2)'.^2)
    dt     = minimum(min(dx, dy)^2 ./inn(H).^npow./4.1)  # time step (obeys ~CFL condition)
    t = 0.0; it = 0
    # Physical time loop
    while t<ttot
        qHx    .= -av_xi(H).^npow.*diff(H[:,2:end-1], dims=1)/dx  # flux
        qHy    .= -av_yi(H).^npow.*diff(H[2:end-1,:], dims=2)/dy  # flux
        dHdt   .= -diff(qHx, dims=1)/dx -diff(qHy, dims=2)/dy     # rate of change
        H[2:end-1,2:end-1] .= inn(H) .+ dt.*dHdt      # update rule, sets the BC as H[1]=H[end]=0
        t += dt; it += 1
    end
    @printf("Total time = %1.2f, time steps = %d \n", round(ttot, sigdigits=2), it)
    # Visualize
    if do_visu
        fontsize = 12
        opts = (aspect_ratio=1, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"),
                ticks=nothing, framestyle=:box, titlefontsize=fontsize, titlefont="Courier", colorbar_title="",
                xlabel="Lx", ylabel="Ly", xlims=(xc[1],xc[end]), ylims=(yc[1],yc[end]), clims=(0.,1.))
        display(heatmap(xc, yc, H'; c=:davos, title="explicit diffusion (nt=$it)", opts...))
        if save_fig savefig("diff2D_expl.png") end
    end
    return xc, yc, H
end

diffusion_2D_expl(; do_visu=do_visu);


text:
# this file is just a bunch of ugly code to make sure that the Julia Array{UInt8,1} becomes a JS Uint8Array() instead of a normal array. This improves performance of the client.
# ignore it if you are not interested in that kind of stuff

import UUIDs: UUID
import MsgPack

# MsgPack.jl doesn't define a serialization method for MIME and UUID objects, so we these ourselves:
MsgPack.msgpack_type(m::Type{<:MIME}) = MsgPack.StringType()
MsgPack.msgpack_type(u::Type{UUID}) = MsgPack.StringType()
MsgPack.to_msgpack(::MsgPack.StringType, m::MIME) = string(m)
MsgPack.to_msgpack(::MsgPack.StringType, u::UUID) = string(u)

# We want typed integer arrays to arrive as JS typed integer arrays:
const JSTypedIntSupport = [Int8, UInt8, Int16, UInt16, Int32, UInt32, Float32, Float64]
JSTypedInt = Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Float32,Float64}

MsgPack.msgpack_type(m::Type{Vector{T}}) where T <: JSTypedInt = MsgPack.ExtensionType()
MsgPack.to_msgpack(::MsgPack.ExtensionType, x::Vector{T}) where T <: JSTypedInt = let
    type = findfirst(isequal(T), JSTypedIntSupport) + 0x10
    MsgPack.Extension(type, reinterpret(UInt8, x))
end

# The other side does the same (/frontend/common/MsgPack.js), and we decode it here:
function decode_extension_and_addbits(x::MsgPack.Extension)
    julia_type = JSTypedIntSupport[x.type - 0x10]
    reinterpret(julia_type, x.data)
end

function decode_extension_and_addbits(x::Dict)
    # we mutate the dictionary, that's fine in our use case and saves memory?
    for (k, v) in x
        x[k] = decode_extension_and_addbits(v)
    end
    x
end
decode_extension_and_addbits(x::Array) = map(decode_extension_and_addbits, x)

# We also convert everything (except the JS typed arrays) to 64 bit numbers, just to make it easier to work with.
decode_extension_and_addbits(x::T) where T <: Union{Signed,Unsigned} = Int64(x)
decode_extension_and_addbits(x::T) where T <: AbstractFloat = Float64(x)
decode_extension_and_addbits(x::Any) = x

function pack(args...)
    MsgPack.pack(args...)
end

function unpack(args...)
    MsgPack.unpack(args...) |> decode_extension_and_addbits
end


text:
# Overloadable method, e.g to insert OffsetPrecalc's precalculated stride multiples
@inline tdot(ptr::AbstractStridedPointer{T}, a, b, c) where {T} = tdot(T, a, b, c)

@inline tdot(::Type{T}, a::Tuple{A}, b::Tuple{B,Vararg}, ::Tuple{Val{false},Vararg}) where {T,A,B} = lazymul(first(a), first(b))
@inline tdot(::Type{T}, a::Tuple{A}, b::Tuple{B,Vararg}, ::Tuple{Val{true},Vararg}) where {T,A,B} = lazymul_no_promote(T, first(a), first(b))

@inline function tdot(::Type{T}, a::Tuple{A1,A2,Vararg}, b::Tuple{B1,B2,Vararg}, c::Tuple{Val{false},Vararg}) where {T,A1,A2,B1,B2}
    vadd(lazymul(first(a), first(b)), tdot(T,Base.tail(a), Base.tail(b), Base.tail(c)))
end
@inline function tdot(::Type{T}, a::Tuple{A1,A2,Vararg}, b::Tuple{B1,B2,Vararg}, c::Tuple{Val{true},Vararg}) where {T,A1,A2,B1,B2}
    vadd(lazymul_no_promote(T,first(a), first(b)), tdot(T,Base.tail(a), Base.tail(b), Base.tail(c)))
end

@inline function tdot(::Type{T}, a::Tuple{A}, b::Tuple{B,Vararg}, c::Tuple{C,Vararg}, ::Tuple{Val{false},Vararg}) where {T,A,B,C}
    lazymul(first(a), first(b), first(c))
end
# @inline tdot(::Type{T}, a::Tuple{A}, b::Tuple{B,Vararg}, c::Tuple{C,Vararg}, ::Tuple{Val{true},Vararg}) where {T,A,B,C} = lazymul_no_promote(T, first(a), first(b))
@inline function tdot(::Type{T}, a::Tuple{A}, b::Tuple{B,Vararg}, c::Tuple{C,Vararg}, ::Tuple{Val{true},Vararg}) where {T,A,B,C}
    lazymul_no_promote(T,first(a), first(b), first(c))
end

@inline function tdot(::Type{T}, a::Tuple{A1,A2,Vararg}, b::Tuple{B1,B2,Vararg}, c::Tuple{C1,C2,Vararg}, d::Tuple{Val{false},Vararg}) where {T,A1,A2,B1,B2,C1,C2}
    vadd(lazymul(first(a), first(b), first(c)), tdot(T,Base.tail(a), Base.tail(b), Base.tail(c), Base.tail(d)))
end
# @inline function tdot(::Type{T}, a::Tuple{A1,A2,Vararg}, b::Tuple{B1,B2,Vararg}, c::Tuple{C1,C2,Vararg}, d::Tuple{Val{true},Vararg}) where {T,A1,A2,B1,B2,C1,C2}
#     vadd(lazymul_no_promote(T,first(a), first(b)), tdot(T,Base.tail(a), Base.tail(b), Base.tail(c), Base.tail(d)))
# end
@inline function tdot(::Type{T}, a::Tuple{A1,A2,Vararg}, b::Tuple{B1,B2,Vararg}, c::Tuple{C1,C2,Vararg}, d::Tuple{Val{true},Vararg}) where {T,A1,A2,B1,B2,C1,C2}
    vadd(lazymul_no_promote(T,first(a), first(b), first(c)), tdot(T,Base.tail(a), Base.tail(b), Base.tail(c), Base.tail(d)))
end


text:
"""
    AbstractCache

Abstract supertype used for information storage object the directed acyclic graph.
"""
abstract type AbstractCache end

"""
    AbstractCacheAttribute

Abstract supertype used for attributes stored in a cache.
"""
abstract type AbstractCacheAttribute end

Base.@kwdef mutable struct VariableValues{T<:Real}
    x0::Vector{T}                          = T[]
    x::Vector{T}                           = T[]
    lower_variable_bounds::Vector{T}       = T[]
    upper_variable_bounds::Vector{T}       = T[]
    node_to_variable_map::Vector{Int}      = Int[]
    variable_to_node_map::Vector{Int}      = Int[]
    variable_types::Vector{VariableType}   = VariableType[]
end

@inline val(b::VariableValues{T}, i::Int) where T = @inbounds b.x[i]
@inline lbd(b::VariableValues{T}, i::Int) where T = @inbounds b.lower_variable_bounds[i]
@inline ubd(b::VariableValues{T}, i::Int) where T = @inbounds b.upper_variable_bounds[i]
val(b::VariableValues{T}) where T = b.x
lbd(b::VariableValues{T}) where T = b.lower_variable_bounds
ubd(b::VariableValues{T}) where T = b.upper_variable_bounds

function _get_x!(::Type{BranchVar}, out::Vector{T}, v::VariableValues{T}) where T<:Real
    @inbounds for i = 1:length(v.node_to_variable_map)
        out[i] = v.x[v.node_to_variable_map[i]]
    end
    return nothing
end

function _initialize_or_copy!(y::VariableValues{T}, x::VariableValues{T}, s::Symbol) where {T<:Real}
    isempty(getfield(y, s)) ? setfield!(y, s, copy(getfield(x, s))) : copy!(getfield(y, s), getfield(x, s))
end
function update_box_and_pnt!(y::VariableValues{T}, x::VariableValues{T}, update_box::Bool) where {T<:Real}
    if update_box
        _initialize_or_copy!(y, x, :x0)
        _initialize_or_copy!(y, x, :lower_variable_bounds)
        _initialize_or_copy!(y, x, :upper_variable_bounds)
    end
    _initialize_or_copy!(y, x, :x)
end

"""
    DirectedTree

A tree graph with a single sink node.
"""
Base.@kwdef mutable struct DirectedTree <: AbstractDirectedAcyclicGraph
    "List of nodes"
    nodes::Vector{Node}                       = Node[]
    "List of index of variables in this tree"
    variables::Dict{Int,Int}                  = Dict{Int,Int}()
    "Information on all variables..."
    v::VariableValues{Float64}                = VariableValues{Float64}()
    "List of constant values"
    constant_values::Vector{Float64}          = Float64[]
    "List of constant values"
    parameter_values::Vector{Float64}         = Float64[]
    "Number of nodes"
    node_count::Int                           = 0
    "Number of variables"
    variable_count::Int                       = 0
    "Number of constants"
    constant_count::Int                       = 0
    sink_bnd::Interval{Float64}               = Interval{Float64}(-Inf,Inf)
    ""
    sparsity::Vector{Int}                     = Int[]
    ""
    rev_sparsity::Dict{Int,Int}               = Dict{Int,Int}()
    dependent_variable_count::Int             = 0
    dep_subexpr_count::Int                    = 0
    dependent_subexpressions::Vector{Int}     = Int[]
    linearity::Linearity                      = LIN_CONSTANT
    user_operators::OperatorRegistry          = OperatorRegistry()
    #children::SpraseMatrixCSC{Bool,Int}      = spzeros(Bool,1,1)
    #parents::SparseMatrixCSC{Bool,Int}       = spzeros(Bool,1,1)
end
#DirectedTree(n::Int) = DirectedTree(children = spzeros(Bool,n,n), parents=spzeros(Bool,n,n))
const DAT = DirectedTree

# node property access functions that can be defined at abstract type
node(g::DAT, i)            = g.nodes[i]
nodes(g::DAT)              = g.nodes

variable(g::DAT, i)        = g.variables[i]
variables(g::DAT)          = g.variables

constant_value(g::DAT, i)  = g.constant_values[i]
constant_values(g::DAT)    = g.constant_values
parameter_value(g::DAT, i) = g.parameter_values[i]
parameter_values(g::DAT)   = g.parameter_values

node_class(g::DAT, i)      = node_class(node(g, i))
ex_type(g::DAT, i)         = ex_type(node(g, i))
first_index(g::DAT, i)     = first_index(node(g, i))
secondary_index(g::DAT, i) = secondary_index(node(g, i))
arity(g::DAT, i)           = arity(node(g, i))
children(g::DAT, i)        = children(node(g, i))
child(g::DAT, i, j)        = child(node(g, j), i)

is_binary(g::DAT, i) = arity(g, i) == 2

node_count(g::DAT)     = g.node_count
variable_count(g::DAT) = g.variable_count
constant_count(g::DAT) = g.constant_count

dep_subexpr_count(g::DAT)            = length(g.dependent_subexpressions)
sparsity(g::DAT, i)                  = g.sparsity
rev_sparsity(g::DAT, i::Int, k::Int) = g.rev_sparsity[i]

user_univariate_operator(g::DAT, i) = g.user_operators.univariate_operator_f[i]
user_multivariate_operator(g::DAT, i) = g.user_operators.multivariate_operator_evaluator[i]

function DirectedTree(aux_info, d, op::OperatorRegistry, sub_sparsity::Dict{Int,Vector{Int}}, subexpr_linearity, parameter_values, is_sub, subexpr_indx)

    nd = copy(d.nd)
    adj = copy(d.adj)
    const_values = copy(d.const_values)

    sparsity, dependent_subexpressions = _compute_sparsity(d, sub_sparsity, is_sub, subexpr_indx)
    rev_sparsity = Dict{Int,Int}()
    for (i,s) in enumerate(sparsity)
        rev_sparsity[s] = i
    end

    nodes = _convert_node_list(aux_info, d.nd, op)
    lin = linearity(nd, adj, subexpr_linearity)
    DirectedTree(nodes = nodes,
                    variables = rev_sparsity,
                    constant_values = const_values,
                    parameter_values = parameter_values,
                    node_count = length(nodes),
                    variable_count = length(sparsity),
                    constant_count = length(const_values),
                    sparsity = sparsity,
                    rev_sparsity = rev_sparsity,
                    dependent_variable_count = length(sparsity),
                    dep_subexpr_count = length(dependent_subexpressions),
                    dependent_subexpressions = copy(dependent_subexpressions),
                    linearity = lin[1],
                    user_operators = op
                    )
end

forward_uni = [i for i in instances(AtomType)]
setdiff!(forward_uni, [VAR_ATOM; PARAM_ATOM; CONST_ATOM; SELECT_ATOM; SUBEXPR])
f_switch = binary_switch(forward_uni, is_forward = true)
@eval function fprop!(t::T, ex::Expression, g::DAT, c::AbstractCache , k::Int) where T<:AbstractCacheAttribute
    id = ex_type(g, k)
    $f_switch
    error("fprop! for ex_type = $id not defined.")
    return
end

reverse_uni = [i for i in instances(AtomType)]
setdiff!(reverse_uni, [VAR_ATOM; PARAM_ATOM; CONST_ATOM; SELECT_ATOM; SUBEXPR])
r_switch = binary_switch(reverse_uni, is_forward = false)
@eval function rprop!(t::T, ex::Expression, g::DAT, c::AbstractCache, k::Int) where T<:AbstractCacheAttribute
    id = ex_type(g, k)
    $r_switch
    error("rprop! for ex_type = $id not defined.")
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# --- Utilities

ieltype(T) =
    if Base.IteratorEltype(T) isa Base.HasEltype
        eltype(T)
    else
        Any
    end

const DenseSubVector{T} =
    SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true}

identityof(::typeof(+), e) = zero(e)
identityof(::typeof(*), e) = one(e)
identityof(::typeof(min), e) = typemax(e)
identityof(::typeof(max), e) = typemin(e)
identityof(::typeof(append!), e) = empty(e)
ridentityof(::typeof(append!), e) = ()

_cljapiurl(name) =
    "https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/$name"
_cljref(name) =
    "[`$name` in Clojure]($(_cljapiurl(name)))"
_thx_clj(name) =
    "This API is modeled after $(_cljref(name))."


text:
function objective_min_redispatch_cost(pm::GenericPowerModel)
    PowerModels.check_cost_models(pm)

    gen_cost = Dict()
    for (n, nw_ref) in nws(pm)
        for (i,gen) in nw_ref[:gen]
            pg_delta = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
            if length(gen["cost"]) == 1
                gen_cost[(n,i)] = gen["cost"][1]
            elseif length(gen["cost"]) == 2
                gen_cost[(n,i)] = gen["cost"][1]*pg_delta + gen["cost"][2]
            elseif length(gen["cost"]) == 3
                gen_cost[(n,i)] = gen["cost"][1]*pg_delta^2 + gen["cost"][2]*pg_delta + gen["cost"][3]
            else
                gen_cost[(n,i)] = 0.0
            end
        end
    end
    return @objective(pm.model, Min,
        sum(
            sum( gen_cost[(n,i)] for (i,gen) in nw_ref[:gen] ) +
            sum(
                load["voll"]*sum( PowerModels.var(pm, n, c, :pl_delta, i) for c in PowerModels.conductor_ids(pm, n)) for (i,load) in PowerModels.ref(pm, :load))
        for (n, nw_ref) in PowerModels.nws(pm))
    )
end

function objective_total_risk(pm::GenericPowerModel, first_stage_network_id, second_stage_network_ids)
    constraint_first_stage_cost(pm::GenericPowerModel, first_stage_network_id)
    constraint_second_stage_risk(pm::GenericPowerModel, second_stage_network_ids)
    for (n, network) in pm.ref[:nw]
        constraint_dispatch_cost(pm::GenericPowerModel, n)
        constraint_redispatch_cost(pm::GenericPowerModel, n)
        constraint_loadshedding_cost(pm::GenericPowerModel, n)
    end
    first_stage_cost = PowerModels.var(pm, :first_stage_cost)
    second_stage_risk = PowerModels.var(pm, :second_stage_risk)

    return @objective(pm.model, Min, first_stage_cost + second_stage_risk)
end

function constraint_first_stage_cost(pm::GenericPowerModel, first_stage_network_id; cnd::Int = pm.ccnd)
    pg_delta = PowerModels.var(pm, first_stage_network_id, cnd, :pg_delta)
    pl_delta = PowerModels.var(pm, first_stage_network_id, cnd, :pl_delta)
    first_stage_cost = PowerModels.var(pm, :first_stage_cost)

    gen_cost_fs = Dict()
    n = first_stage_network_id
    for (i,gen) in pm.ref[:nw][n][:gen]
        pg_delta_fs = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
        if length(gen["cost"]) == 1
            gen_cost_fs[(n,i)] = 0
        elseif length(gen["cost"]) == 2
            gen_cost_fs[(n,i)] = gen["cost"][1]*pg_delta_fs + gen["cost"][2]
        elseif length(gen["cost"]) == 3
            gen_cost_fs[(n,i)] = gen["cost"][1]*pg_delta_fs^2 + gen["cost"][2]*pg_delta_fs + gen["cost"][3]
        else
            gen_cost_fs[(n,i)] = 0.0
        end
    end

    return  @constraint(pm.model, first_stage_cost ==
        sum(gen_cost_fs[(n,i)] for (i,gen) in PowerModels.ref(pm, :gen)) +
        sum(load["voll"]*pl_delta[i] for (i,load) in PowerModels.ref(pm, :load)))
end

function constraint_first_stage_cost(pm::GenericPowerModel{T}, first_stage_network_id; cnd::Int = pm.ccnd)  where T <: PowerModels.AbstractDCPForm
    pg_delta = PowerModels.var(pm, first_stage_network_id, cnd, :pg_delta)
    pl_delta = PowerModels.var(pm, first_stage_network_id, cnd, :pl_delta)
    first_stage_cost = PowerModels.var(pm, :first_stage_cost)

    gen_cost_fs = Dict()
    n = first_stage_network_id
    for (i,gen) in pm.ref[:nw][n][:gen]
        pg_delta_fs = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
        if length(gen["cost"]) == 1
            gen_cost_fs[(n,i)] = 0
        elseif length(gen["cost"]) == 2
            gen_cost_fs[(n,i)] = gen["cost"][1]*pg_delta_fs + gen["cost"][2]
        elseif length(gen["cost"]) == 3
            gen_cost_fs[(n,i)] = gen["cost"][2]*pg_delta_fs + gen["cost"][3]
        else
            gen_cost_fs[(n,i)] = 0.0
        end
    end

    return  @constraint(pm.model, first_stage_cost ==
        sum(gen_cost_fs[(n,i)] for (i,gen) in PowerModels.ref(pm, :gen)) +
        sum(load["voll"]*pl_delta[i] for (i,load) in PowerModels.ref(pm, :load)))
end

function constraint_second_stage_risk(pm::GenericPowerModel, second_stage_network_ids; cnd::Int = pm.ccnd)
    pg_delta = Dict(n => PowerModels.var(pm, n, cnd, :pg_delta) for (n, contingency_id) in second_stage_network_ids)
    pl_delta = Dict(n => PowerModels.var(pm, n, cnd, :pl_delta) for (n, contingency_id) in second_stage_network_ids)
    second_stage_risk = PowerModels.var(pm, :second_stage_risk)

    gen_cost_ss = Dict()
    for (cont_id, n)  in second_stage_network_ids
        for (i, gen) in pm.ref[:nw][n][:gen]
            pg_delta_ss = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
            if length(gen["cost"]) == 1
                gen_cost_ss[(n,i)] = 0
            elseif length(gen["cost"]) == 2
                gen_cost_ss[(n,i)] = gen["cost"][1]*pg_delta_ss + gen["cost"][2]
            elseif length(gen["cost"]) == 3
                gen_cost_ss[(n,i)] = gen["cost"][1]*pg_delta_ss^2 + gen["cost"][2]*pg_delta_ss + gen["cost"][3]
            else
                gen_cost_ss[(n,i)] = 0.0
            end
        end
    end
    return @constraint(pm.model, second_stage_risk ==
    sum(    PowerModels.ref(pm, n, :contingencies, cont_id)["prob"] *
            (
            sum(gen_cost_ss[(n,i)] for (i,gen) in PowerModels.ref(pm, n, :gen))  +
            sum(load["voll"]*pl_delta[n][i] for (i,load) in PowerModels.ref(pm, n, :load))
            )
            for (cont_id, n)  in second_stage_network_ids
    ))
end

function constraint_second_stage_risk(pm::GenericPowerModel{T}, second_stage_network_ids; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    pg_delta = Dict(n => PowerModels.var(pm, n, cnd, :pg_delta) for (n, contingency_id) in second_stage_network_ids)
    pl_delta = Dict(n => PowerModels.var(pm, n, cnd, :pl_delta) for (n, contingency_id) in second_stage_network_ids)
    second_stage_risk = PowerModels.var(pm, :second_stage_risk)

    gen_cost_ss = Dict()
    for (cont_id, n)  in second_stage_network_ids
        for (i, gen) in pm.ref[:nw][n][:gen]
            pg_delta_ss = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
            if length(gen["cost"]) == 1
                gen_cost_ss[(n,i)] = 0
            elseif length(gen["cost"]) == 2
                gen_cost_ss[(n,i)] = gen["cost"][1]*pg_delta_ss + gen["cost"][2]
            elseif length(gen["cost"]) == 3
                gen_cost_ss[(n,i)] = gen["cost"][2]*pg_delta_ss + gen["cost"][3]
            else
                gen_cost_ss[(n,i)] = 0.0
            end
        end
    end
    return @constraint(pm.model, second_stage_risk ==
    sum(    PowerModels.ref(pm, n, :contingencies, cont_id)["prob"] *
            (
            sum(gen_cost_ss[(n,i)] for (i,gen) in PowerModels.ref(pm, n, :gen))  +
            sum(load["voll"]*pl_delta[n][i] for (i,load) in PowerModels.ref(pm, n, :load))
            )
            for (cont_id, n)  in second_stage_network_ids
    ))
end
""
function constraint_dispatch_cost(pm::GenericPowerModel, n; cnd::Int = pm.ccnd)
    gen_cost = Dict()
    dispatch_cost = PowerModels.var(pm, :dispatch_cost)
    for (i,gen) in pm.ref[:nw][n][:gen]
        pg = sum( var(pm, n, c, :pg, i) for c in conductor_ids(pm, n) )
        if length(gen["cost"]) == 1
            gen_cost[(i)] = gen["cost"][1]
        elseif length(gen["cost"]) == 2
            gen_cost[(i)] = gen["cost"][1]*pg + gen["cost"][2]
        elseif length(gen["cost"]) == 3
            gen_cost[(i)] = gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
        else
            gen_cost[(i)] = 0.0
        end
    end
    return @constraint(pm.model, dispatch_cost[n] ==
            sum( gen_cost[(i)] for (i,gen) in pm.ref[:nw][n][:gen] )
    )
end

function constraint_dispatch_cost(pm::GenericPowerModel{T}, n; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    gen_cost = Dict()
    dispatch_cost = PowerModels.var(pm, :dispatch_cost)
    for (i,gen) in pm.ref[:nw][n][:gen]
        pg = sum( var(pm, n, c, :pg, i) for c in conductor_ids(pm, n) )
        if length(gen["cost"]) == 1
            gen_cost[(i)] = gen["cost"][1]
        elseif length(gen["cost"]) == 2
            gen_cost[(i)] = gen["cost"][1]*pg + gen["cost"][2]
        elseif length(gen["cost"]) == 3
            gen_cost[(i)] = gen["cost"][2]*pg + gen["cost"][3]
        else
            gen_cost[(i)] = 0.0
        end
    end
    return @constraint(pm.model, dispatch_cost[n] ==
            sum( gen_cost[(i)] for (i,gen) in pm.ref[:nw][n][:gen] )
    )
end

function constraint_redispatch_cost(pm::GenericPowerModel, n; cnd::Int = pm.ccnd)
    redispatch_cost = PowerModels.var(pm, :redispatch_cost)
    gen_cost = Dict()

    for (i,gen) in pm.ref[:nw][n][:gen]
        pg_delta = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
        if length(gen["cost"]) == 1
            gen_cost[(i)] = 0
        elseif length(gen["cost"]) == 2
            gen_cost[(i)] = gen["cost"][1]*pg_delta
        elseif length(gen["cost"]) == 3
            gen_cost[(i)] = gen["cost"][1]*pg_delta^2 + gen["cost"][2]*pg_delta
        else
            gen_cost[(i)] = 0.0
        end
    end

    return @constraint(pm.model, redispatch_cost[n] ==
            sum( gen_cost[(i)] for (i,gen) in pm.ref[:nw][n][:gen] )
    )
end

function constraint_redispatch_cost(pm::GenericPowerModel{T}, n; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    redispatch_cost = PowerModels.var(pm, :redispatch_cost)
    gen_cost = Dict()

    for (i,gen) in pm.ref[:nw][n][:gen]
        pg_delta = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
        if length(gen["cost"]) == 1
            gen_cost[(i)] = 0
        elseif length(gen["cost"]) == 2
            gen_cost[(i)] = gen["cost"][1]*pg_delta
        elseif length(gen["cost"]) == 3
            gen_cost[(i)] = gen["cost"][1]*pg_delta^2 + gen["cost"][2]*pg_delta
        else
            gen_cost[(i)] = 0.0
        end
    end

    return @constraint(pm.model, redispatch_cost[n] ==
            sum( gen_cost[(i)] for (i,gen) in pm.ref[:nw][n][:gen] )
    )
end

function constraint_loadshedding_cost(pm::GenericPowerModel, n; cnd::Int = pm.ccnd)
    pl_delta = PowerModels.var(pm, n, cnd, :pl_delta)
    loadshedding_cost = PowerModels.var(pm, :loadshedding_cost)

    return @constraint(pm.model, loadshedding_cost[n] ==
            sum(load["voll"]*pl_delta[i] for (i,load) in PowerModels.ref(pm, n,:load))
    )
end

# function objective_total_risk(pm::GenericPowerModel{T}, first_stage_network_id, second_stage_network_ids; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
#
#     for (n, network) in pm.ref[:nw]
#         constraint_loadshedding_cost(pm::GenericPowerModel, n)
#     end
#     pg_delta = PowerModels.var(pm, first_stage_network_id, cnd, :pg_delta)
#     pl_delta = PowerModels.var(pm, first_stage_network_id, cnd, :pl_delta)
#
#     pg_delta_ss = Dict(n => PowerModels.var(pm, n, cnd, :pg_delta) for (n, contingency_id) in second_stage_network_ids)
#     pl_delta_ss = Dict(n => PowerModels.var(pm, n, cnd, :pl_delta) for (n, contingency_id) in second_stage_network_ids)
#
#     gen_cost_fs = Dict()
#     gen_cost_ss = Dict()
#     for n in first_stage_network_id
#         for (i,gen) in pm.ref[:nw][n][:gen]
#             pg_fs = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
#             if length(gen["cost"]) == 1
#                 gen_cost_fs[(i)] = 0
#             elseif length(gen["cost"]) == 2
#                 gen_cost_fs[(i)] = gen["cost"][1]*pg_fs
#             elseif length(gen["cost"]) == 3
#                 gen_cost_fs[(i)] = gen["cost"][2]*pg_fs
#             else
#                 gen_cost_fs[(i)] = 0.0
#             end
#         end
#     end
#     for (cont_id, n) in second_stage_network_ids
#         for (i,gen) in pm.ref[:nw][n][:gen]
#             pg_ss = sum( var(pm, n, c, :pg_delta, i) for c in conductor_ids(pm, n) )
#             if length(gen["cost"]) == 1
#                 gen_cost_ss[(n,i)] = 0
#             elseif length(gen["cost"]) == 2
#                 gen_cost_ss[(n,i)] = gen["cost"][1]*pg_ss
#             elseif length(gen["cost"]) == 3
#                 gen_cost_ss[(n,i)] = gen["cost"][2]*pg_ss
#             else
#                 gen_cost_ss[(n,i)] = 0.0
#             end
#         end
#     end
#
#     return @objective(pm.model, Min,
#             sum(gen_cost_fs[(i)] for (i,gen) in PowerModels.ref(pm, :gen)) +
#             sum(load["voll"]*pl_delta[i] for (i,load) in PowerModels.ref(pm, :load)) +
#             sum(    pm.ref[:nw][n][:contingencies][contingency_id]["prob"] *
#             (
#             sum(gen_cost_ss[(n,i)] for (i,gen) in pm.ref[:nw][n][:gen]) +
#             sum(load["voll"]*pl_delta_ss[n][i] for (i,load) in pm.ref[:nw][n][:load])
#             )
#             for (contingency_id, n) in second_stage_network_ids
#     ))
# end

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
module BeaData

import HTTP
import JSON
using Dates
using DataFrames, DataStructures

import Base.show

export
    # Types
    Bea,
    BeaTable,
    # Methods
    get_bea_datasets,
    get_bea_parameterlist,
    get_nipa_table

const DEFAULT_API_URL = "https://APPS.BEA.GOV/api/data"
const API_KEY_LENGTH  = 36

"""

A connection to the U.S. Bureau of Economic Analysis (BEA) Data API.

Constructors
------------
* `Bea()`
* `Bea(key::AbstractString)`

Arguments
---------
* `key`: Registration key provided by the BEA.

A valid registration key is required to retrieve data from the BEA's API.  A key can be obtained by registering at the BEA website.

A default API key can be specified in a ~/.beadatarc file.

"""
mutable struct Bea
    url::AbstractString
    key::AbstractString
end

Bea(key) = Bea(DEFAULT_API_URL, key)

function Bea()
    key = ""
    if "BEA_KEY" in keys(ENV)
        key = ENV["BEA_KEY"]
    elseif isfile(joinpath(homedir(), ".beadatarc"))
        open(joinpath(homedir(),".beadatarc"), "r") do f
            key = read(f, String)
        end
        key = rstrip(key)
    else
        error("No API key found, connection not initialized")
    end

    println("API key loaded.")

    # Key validation
    if length(key) > API_KEY_LENGTH || length(key) < API_KEY_LENGTH
        error("Invalid key length (≠ ", API_KEY_LENGTH, " chars), connection not initialized")
    end

    return Bea(key)
end

function Base.show(io::IO, b::Bea)
    println(io, "BEA API Connection")
    println(io, "url: $(b.url)")
    println(io, "key: $(b.key)")
end

"""
A NIPA table with data and metadata returned from a [`get_nipa_table`](@ref) call.

Fields
---

* tablenum - NIPA table number
* tablename - API TableName
* tabledesc - The table title (e.g., 'Real Gross Domestic Product, Chained Dollars' for Table 1.1.6)
* linedesc - `OrderedDict` of descriptions for each line of the table
* tablenotes - Table notes, if any
* frequency
* startyear
* endyear
* df - `DataFrame` containing the data values from the table; column names are the line numbers from the table, the first column contains the date for each observation in Julia `Date` format

"""
mutable struct BeaTable
    tablenum::AbstractString
    tablename::AbstractString
    tabledesc::AbstractString
    linedesc::OrderedDict
    tablenotes::Any
    frequency::AbstractString
    startyear::Int
    endyear::Int
    df::DataFrame
end

function Base.show(io::IO, b::BeaTable)
    println(io, "BEA NIPA Table")
    println(io, "Table: $(b.tablenum)")
    println(io, "TableName: $(b.tablename)")
    println(io, "Description: $(b.tabledesc)")
    println(io, "Coverage: $(b.frequency), from $(b.startyear) to $(b.endyear)")
end

include("get_bea_datasets.jl")
include("get_bea_parameterlist.jl")
include("get_nipa_table.jl")

end # module


text:
using Dubins

q0 = [0.0, 0.0, 0.0]
q1 = [3.0, 3.0, 0.0]
ρ = 1.0

errcode, path = dubins_shortest_path(q0, q1, ρ)
path_length = dubins_path_length(path)

sample_times = linspace(0, path_length, 100)

f = open("./sample.tex", "w")
write(f,"\\documentclass[10pt]{standalone}\n\\usepackage[T1]{fontenc}\n\\usepackage{pgfplots}\n\\pgfplotsset{compat=newest}
\\begin{document}
\\begin{tikzpicture}
\\begin{axis}[axis equal,grid=both,title={Dubins' Path}]
\\addplot[only marks, mark=*,color=red] plot coordinates {
($(q0[1]), $(q0[2])) ($(q1[1]), $(q1[2]))};
\\addplot[smooth, mark=none,color=red] plot coordinates {\n ($(q0[1]), $(q0[2]))  ")

for t in sample_times
    q = dubins_path_sample(path, t)
    write(f,"($(q[1]), $(q[2])) ")
end
write(f,"($(q1[1]), $(q1[2]))")
write(f,"};\n\\end{axis}\n\\end{tikzpicture}\n\\end{document}")


text:
# Ordinal Models
@testset "Ordinal Models" begin
    data =
        dataset("Ecdat", "Kakadu") |>
        (data -> select(data, [:RecParks, :Sex, :Age, :Schooling]))
    data.RecParks = convert(Vector{Int}, data.RecParks)
    data.RecParks = levels!(categorical(data.RecParks, ordered = true, compress = true), collect(1:5))
    model = fit(EconometricModel, @formula(RecParks ~ Age + Sex + Schooling), data)
    @test sprint(show, model) ==
        "Probability Model for Ordinal Response\nCategories: 1 < 2 < 3 < 4 < 5\nNumber of observations: 1827\nNull Loglikelihood: -2677.60\nLoglikelihood: -2657.40\nR-squared: 0.0075\nLR Test: 40.42 ∼ χ²(3) ⟹ Pr > χ² = 0.0000\nFormula: RecParks ~ Age + Sex + Schooling\n──────────────────────────────────────────────────────────────────────────────────────────\n                          PE          SE        t-value  Pr > |t|        2.50%      97.50%\n──────────────────────────────────────────────────────────────────────────────────────────\nAge                  0.00943647  0.00254031    3.71469     0.0002   0.00445424   0.0144187\nSex: male           -0.0151659   0.0846365    -0.179188    0.8578  -0.181161     0.150829\nSchooling           -0.103902    0.0248742    -4.17711     <1e-04  -0.152687    -0.0551174\n(Intercept): 1 | 2  -2.92405     0.191927    -15.2352      <1e-48  -3.30047     -2.54763\n(Intercept): 2 | 3  -1.54922     0.171444     -9.03632     <1e-18  -1.88547     -1.21297\n(Intercept): 3 | 4  -0.298938    0.166904     -1.79108     0.0734  -0.626281     0.0284051\n(Intercept): 4 | 5   0.669835    0.167622      3.9961      <1e-04   0.341083     0.998587\n──────────────────────────────────────────────────────────────────────────────────────────"
    β, V, σ = coef(model), vcov(model), stderror(model)
    @test β ≈ [
        0.009437926,
        -0.015143049,
        -0.103911316,
        -2.92391240,
        -1.549266200,
        -0.298963900,
        0.6698249,
        ] rtol = 1e-4
    @test V ≈ [
        0.00000646 -0.0000007000 0.000013250 0.00031473 0.0003161358 0.000321276 0.0003270734
        -0.00000070 0.0071633800 -0.000184930 0.00294759 0.0029647442 0.002995366 0.0029904485
        0.00001325 -0.0001849300 0.000618730 0.00288365 0.0028548209 0.002791361 0.0027318928
        0.00031473 0.0029475900 0.002883650 0.03684071 0.0289239290 0.026986237 0.026498765
        0.00031614 0.0029647400 0.002854820 0.02892393 0.0293971651 0.027160433 0.0266006675
        0.000321276 0.0029953660 0.002791361 0.02698623 0.0271604331 0.027860143 0.0270590982
        0.0003270734 0.002990448 0.002731893 0.02649876 0.0266006675 0.027059098 0.0281002038
        ] rtol = 1e-3
    @test σ ≈ [
        0.002540789,
        0.084636770,
        0.024874381,
        0.191939349,
        0.171456015,
        0.166913580,
        0.167631154,
        ] rtol = 1e-4
    data = data[data.RecParks .≠ 5, :]
    model = fit(EconometricModel, @formula(RecParks ~ Age + Sex + Schooling), data)
    @test sprint(show, model) ==
        "Probability Model for Ordinal Response\nCategories: 1 < 2 < 3 < 4\nNumber of observations: 1200\nNull Loglikelihood: -1502.61\nLoglikelihood: -1500.73\nR-squared: 0.0013\nLR Test: 3.76 ∼ χ²(3) ⟹ Pr > χ² = 0.2886\nFormula: RecParks ~ Age + Sex + Schooling\n───────────────────────────────────────────────────────────────────────────────────────────\n                          PE          SE        t-value  Pr > |t|        2.50%       97.50%\n───────────────────────────────────────────────────────────────────────────────────────────\nAge                  0.00063565  0.00329173    0.193105    0.8469  -0.00582255   0.00709385\nSex: male            0.0896101   0.106068      0.84484     0.3984  -0.118489     0.297709\nSchooling           -0.0509934   0.0300878    -1.69482     0.0904  -0.110024     0.00803731\n(Intercept): 1 | 2  -2.58123     0.228989    -11.2723      <1e-27  -3.0305      -2.13197\n(Intercept): 2 | 3  -1.12053     0.212192     -5.28071     <1e-06  -1.53684     -0.704217\n(Intercept): 3 | 4   0.486702    0.210541      2.31167     0.0210   0.0736316    0.899772\n───────────────────────────────────────────────────────────────────────────────────────────"
    β, V, σ = coef(model), vcov(model), stderror(model)
    @test β ≈ [0.00063565, 0.08961012, -0.05099339, -2.5812334, -1.1205278, 0.48670172] rtol = 1e-4
    @test V ≈ [
        0.00001084 -5.386e-06 0.00002186 0.00052578 0.0005259 0.00052967
        -5.386e-06 0.01125032 -0.00021416 0.00465658 0.00470423 0.00484966
        0.00002186 -0.00021416 0.00090528 0.00435314 0.00431928 0.00424373
        0.00052578 0.00465658 0.00435314 0.05243606 0.04424607 0.04221716
        0.0005259 0.00470423 0.00431928 0.04424607 0.04502565 0.0424849
        0.00052967 0.00484966 0.00424373 0.04221716 0.0424849 0.04432745
        ] rtol = 1e-3
    @test σ ≈ [0.0032917, 0.1060675, 0.0300878, 0.2289892, 0.2121925, 0.2105408] rtol = 1e-4
    data = CSV.read(joinpath(pkgdir(Econometrics), "data", "auto.csv"), DataFrame) |>
        (data -> select(data, [:rep77, :foreign, :length, :mpg])) |>
        dropmissing
    data.foreign = categorical(data.foreign, compress = true)
    data.rep77 = levels!(categorical(data.rep77; ordered = true, compress = true),
        ["Poor", "Fair", "Average", "Good", "Excellent"],
        )
    model = fit(EconometricModel, @formula(rep77 ~ foreign + length + mpg), data)
    @test sprint(show, model) ==
        "Probability Model for Ordinal Response\nCategories: Poor < Fair < Average < Good < Excellent\nNumber of observations: 66\nNull Loglikelihood: -89.90\nLoglikelihood: -78.25\nR-squared: 0.1295\nLR Test: 23.29 ∼ χ²(3) ⟹ Pr > χ² = 0.0000\nFormula: rep77 ~ foreign + length + mpg\n─────────────────────────────────────────────────────────────────────────────────────────────\n                                    PE         SE     t-value  Pr > |t|       2.50%    97.50%\n─────────────────────────────────────────────────────────────────────────────────────────────\nforeign: Foreign                2.89681    0.790641   3.66387    0.0005   1.31684     4.47678\nlength                          0.0828275  0.02272    3.64558    0.0005   0.0374253   0.12823\nmpg                             0.230768   0.0704548  3.2754     0.0017   0.0899749   0.37156\n(Intercept): Poor | Fair       17.9275     5.55119    3.22948    0.0020   6.83431    29.0206\n(Intercept): Fair | Average    19.8651     5.59648    3.54956    0.0007   8.68139    31.0487\n(Intercept): Average | Good    22.1033     5.70894    3.8717     0.0003  10.6949     33.5117\n(Intercept): Good | Excellent  24.6921     5.89075    4.19168    <1e-04  12.9204     36.4639\n─────────────────────────────────────────────────────────────────────────────────────────────"
    β, V, σ = coef(model), vcov(model), stderror(model)
    @test β ≈ [2.89679875, 0.08282676, 0.23076532, 17.92728, 19.86486, 22.10311, 24.69193] atol = 1e-3
    @test V ≈ [
        0.62578335 0.0111044038 0.012321411 2.4421745 2.4836507 2.5750058 2.6957587
        0.0111044 0.0005186744 0.001222385 0.1243671 0.1258979 0.1284727 0.1319772
        0.01232141 0.0012223847 0.004974678 0.3328169 0.3368057 0.3430257 0.3549711
        2.44217451 0.1243670884 0.332816913 30.9544666 31.0279194 31.6005528 32.5099187
        2.48365067 0.1258979463 0.336805731 31.0279194 31.4396158 31.9986912 32.9170873
        2.57500582 0.1284726666 0.343025653 31.6005528 31.9986912 32.7044311 33.6189243
        2.69575871 0.1319772302 0.354971081 32.5099187 32.9170873 33.6189243 34.8533724
        ] rtol = 1e-2
    @test σ ≈ [
        0.79106469,
        0.02277442,
        0.07053140,
        5.56367384,
        5.60710405,
        5.71877881,
        5.90367448,
        ] rtol = 1e-2
end


text:
#!/usr/bin/env julia
using Test
using TreeViewWidget
using WebIO

@testset "TreeViewWidget.jl" begin

    function test_text_node_html(tree_node::TreeViewWidget.AbstractTreeViewNode, html::Node)
        # tree_node.label
        @test html.instanceof.tag == :span
        @test html.props[:id] == tree_node.id
        @test html.children[1] == tree_node.label
    end

    function test_node_html(tree_node::TreeViewWidget.AbstractTreeViewNode, html::Node)
        #

        #   tree_node.label
        #

        @test html.instanceof.tag == :li
        test_text_node_html(tree_node, html.children[1])
    end

    function test_nested_node_html(tree_node::TreeViewWidget.AbstractTreeViewNode, html::Node)
        #

        #
        #   label
        #
children

        #

        @test html.instanceof.tag == :li
        @test length(html.children) == 3
        test_text_node_html(tree_node, html.children[2])
        nested_list = html.children[3].children
        tree_node_children = tree_node.children
        @test length(nested_list) == length(tree_node_children)
        for (child_node, html_element) in zip(tree_node_children, nested_list)
            if isempty(child_node.children)
                test_node_html(child_node, html_element)
            else
                test_nested_node_html(child_node, html_element)
            end
        end
    end

    function test_example_root_html(root::TreeViewWidget.TreeViewRoot, html::Node)
        #

        #
...beverages with children... coffee with children...

        #
...meat with children...

        #
...vegetables...

        #

        @test html isa Node
        @test html.instanceof.tag == :ul
        @test length(html.children) == 3
        @test html.children[1].instanceof.tag == :li
        test_nested_node_html(root.nodes[1], html.children[1])
        test_nested_node_html(root.nodes[2], html.children[2])
        test_node_html(root.nodes[3], html.children[3])
    end

    @testset "TreeViewNode" begin
        tree_node = TreeViewNode("foo")
        @test tree_node.label == "foo"
        @test !isempty(tree_node.id)

        tree_child = TreeViewNode("bar")
        push!(tree_node, tree_child)
        @test first(tree_node.children) == tree_child
    end

    @testset "TreeViewRoot" begin
        root = TreeViewRoot()
        push!(root, TreeViewNode("foo"))
        push!(root, TreeViewNode("bar"))
        @test length(root.nodes) == 2
    end

    @testset "ChildlessNodeToHTML" begin
        tree_node = TreeViewNode("foo")
        html = TreeViewWidget.tohtml(tree_node)
        @test html isa Node
        @test length(html.children) == 1
        test_node_html(tree_node, html)
    end

    @testset "NestedNodeToHTML" begin
        tree_node = TreeViewNode("parent_node")
        push!(tree_node, TreeViewNode("child1"))
        push!(tree_node, TreeViewNode("child2"))
        nr_tree_children = length(tree_node.children)

        html = TreeViewWidget.tohtml(tree_node)
        @test html isa Node
        test_nested_node_html(tree_node, html)
    end

    @testset "DoubleNestedNodeToHTML" begin
        tree_node = TreeViewNode("parent_node")
        push!(tree_node, TreeViewNode("child1"))
        child2 = TreeViewNode("child2")
        push!(tree_node, child2)
        push!(child2, TreeViewNode("grandchild1"))
        push!(child2, TreeViewNode("grandchild2"))
        nr_tree_children = length(tree_node.children)

        html = TreeViewWidget.tohtml(tree_node)
        @test html isa Node
        test_nested_node_html(tree_node, html)

        # explicitly test the 2nd level again to be sure, recursiveness can be difficult
        test_nested_node_html(tree_node.children[2], html.children[3].children[2])
    end

    @testset "TreeviewRootToHTML" begin
        root = TreeViewWidget.example_treeview_root()
        html = TreeViewWidget.tohtml(root)
        test_example_root_html(root, html)
    end

    @testset "TreeView" begin
        example_root = TreeViewWidget.example_treeview_root()
        treeview = TreeView(example_root)
        test_example_root_html(example_root, treeview.scope.dom.children[end])

        # Observerable behavior as used in Interact.jl etc
        @test treeview[] == treeview.selection
    end

    @testset "NodeIterator" begin
        x = TreeViewNode("x")
        @test collect(eachnode(x)) == [x]

        x2 = TreeViewNode("x2"); push!(x, x2)
        x3 = TreeViewNode("x3"); push!(x, x3)
        @test collect(eachnode(x)) == [x, x2, x3]

        x4 = TreeViewNode("x4"); push!(x2, x4)
        x5 = TreeViewNode("x5"); push!(x3, x5)
        @test collect(eachnode(x)) == [x, x2, x4, x3, x5]

        tree = TreeViewRoot([x])
        @test collect(eachnode(tree)) == [x, x2, x4, x3, x5]

        treeview = TreeView(tree)
        @test collect(eachnode(treeview)) == [x, x2, x4, x3, x5]
    end

    @testset "FindById" begin
        tree = TreeViewWidget.example_treeview_root()
        espresso_node = tree.nodes[1].children[3].children[1]
        @test TreeViewWidget.find_by_id(tree, espresso_node.id) == espresso_node
        @test tree[espresso_node.id] == espresso_node

        treeview = TreeView(tree)
        @test treeview[espresso_node.id] == espresso_node

        d = Dict{String, String}()
        d["id"] = espresso_node.id
        @test treeview[d] == espresso_node
    end

end


text:
# Illustrative examples from Flux.jl

using Solaris
using Solaris.Flux, Solaris.DiffEqFlux

"""
gradient: scalar -> scalar
"""
gradient(x -> 3x^2 + 2x + 1, 5) # (32, )
gradient((a, b) -> a * b, 2, 3) # (3, 2)

"""
gradient: scalar -> vector

    f = (x₁ - y₁)² + (x₂ - y₂)²
    ∂f/∂x₁ = 2x₁ - 2y₁
    ⋯
"""
f(x, y) = sum((x .- y) .^ 2)
gradient(f, [2, 1], [2, 0])

"""
params: a syntax sugar for machine learning with thousands of parameters
"""
x = [2, 1]
y = [2, 0]
gs = gradient(params(x, y)) do
    f(x, y)
end
gs[x]
gs[y]

"""
toy model
"""
W = rand(2, 5)
b = rand(2)

predict(x) = W * x .+ b
function loss(x, y)
    ŷ = predict(x)
    sum((y .- ŷ) .^ 2)
end

X, Y = rand(5), rand(2)
L0 = loss(X, Y)

gs = gradient(() -> loss(X, Y), params(W, b))
W̄ = gs[W]
W .-= 0.1 .* W̄
L1 = loss(X, Y)

@assert L1 < L0

"""
affine chain
"""
m = Chain(
    Solaris.Affine(28^2, 128, relu),
    Solaris.Affine(128, 32, relu),
    Solaris.Affine(32, 10),
)

# regularization
sqnorm(x) = sum(abs2, x)
loss(x, y) = logitcrossentropy(m(x), y) + sum(sqnorm, Flux.params(m))

loss(rand(Float32, 28^2), rand(Float32, 10))


text:
# norm1est.jl is adapted from
# https://github.com/JuliaLang/julia/blob/master/base/sparse/linalg.jl
# and is licensed under the MIT License.
#
# Copyright (c) 2009-2015: Jeff Bezanson, Stefan Karpinski, Viral
# B. Shah, and other contributors:
#
# https://github.com/JuliaLang/julia/contributors
#
#   Permission is hereby granted, free of charge, to any person
#   obtaining a copy of this software and associated documentation
#   files (the "Software"), to deal in the Software without
#   restriction, including without limitation the rights to use, copy,
#   modify, merge, publish, distribute, sublicense, and/or sell copies
#   of the Software, and to permit persons to whom the Software is
#   furnished to do so, subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function A_pow_n_B!(res, A, n::Integer, v)
    tmp = similar(v)
    mul!(res, A, v)
    for i in 1:n-1
        mul!(tmp, A, res)
        copyto!(res, tmp)
    end
end

function At_pow_n_B!(res, A, n::Integer, v)
    tmp = similar(v)
    mul!(res, adjoint(A), v)
    for i in 1:n-1
        mul!(tmp, adjoint(A), res)
        copyto!(res, tmp)
    end
end

function norm1est(A, m::Integer, t::Integer = min(2,maximum(size(A))))::Real
    # Effectively implements Algorithm 2.4 of Higham, Tisseur, SIAM J. Mat. Anal. Appl. 21, 1185 (2000)
    # The first argument is the power to which A is raised.

    maxiter = 5
    # Check the input
    @assert size(A,1) == size(A,2) "Matrix exponential only defined for square matrices"

    n = size(A,1)

    if t <= 0
        throw(ArgumentError("number of blocks must be a positive integer"))
    end
    if t > n
        throw(ArgumentError("number of blocks must not be greater than $n"))
    end

    ind = Array{Integer}(undef, n)
    ind_hist = Array{Integer}(undef, maxiter * t)

    #Ti = typeof(float(zero(T)))
    Ti = eltype(A)
    S = zeros(Ti, n, t)

    function _rand_pm1!(v)
        for i in eachindex(v)
            v[i] = rand() < 0.5 ? 1 : -1
        end
    end

    function _any_abs_eq(v, n::Integer)
        for i in eachindex(v)
            if abs(v[i]) == n
                return true
            end
        end
        return false
    end

    # Generate the block matrix
    X = Array{Ti}(undef, n, t)
    Y = similar(X)
    Z = similar(X)

    X[1:n,1] .= 1
    for j = 2:t
        while true
            _rand_pm1!(view(X,1:n,j)) # Chose random
            yaux = X[1:n,j]' * X[1:n,1:j-1]
            if !_any_abs_eq(yaux, n)
                break
            end
        end
    end
    rmul!(X, 1. /n)

    iter = 0
    local est
    local est_old
    est_ind = 0

    while iter < maxiter
        iter += 1
        A_pow_n_B!(Y, A, m, X)

        est = zero(real(eltype(Y)))

        est_ind = 0

        # Find the best estimate by looking at the maximum 1-norm of the columns of Y
        # and save its position
        for i = 1:t
            y = norm(Y[1:n, i], 1)
            if y > est
                est = y
                est_ind = i
            end
        end

        if iter == 1
            est_old = est
        end

        if est > est_old || iter == 2
            ind_best = est_ind
        end

        # If, after the second iteration, the estimate is smaller than before,
        # break, cause we have found the best lower bound
        if iter >= 2 && est <= est_old
            est = est_old
            break
        end

        est_old = est
        S_old = copy(S)

        for j = 1:t
            for i = 1:n
                S[i,j] = Y[i,j] == 0 ? one(Y[i,j]) : sign(Y[i,j])
            end
        end

        # TODO: Check if every column of S is parallel to a column of S_old:
        #       you should break in this case

        if Ti <: Real
            # Check wether cols of S are parallel to cols of S or S_old
            for j = 1:t
                while true
                    repeated = false
                    if j > 1
                        saux = S[1:n,j]' * S[1:n,1:j-1]
                        if _any_abs_eq(saux,n)
                            repeated = true
                        end
                    end
                    if !repeated
                        saux2 = S[1:n,j]' * S_old[1:n,1:t]
                        if _any_abs_eq(saux2,n)
                            repeated = true
                        end
                    end
                    if repeated
                        _rand_pm1!(view(S,1:n,j))
                    else
                        break
                    end
                end
            end
        end

        # Use the conjugate transpose
        At_pow_n_B!(Z, A, m, S)
        h_max = zero(real(eltype(Z)))
        h = zeros(real(eltype(Z)), n)
        h_ind = 0
        for i = 1:n
            h[i] = norm(Z[i,1:t], Inf)
            if h[i] > h_max
                h_max = h[i]
                h_ind = i
            end
            ind[i] = i
        end
        if iter >=2 && ind_best == h_ind
            break
        end
        p = sortperm(h, rev=true)
        h = h[p]
        permute!(ind, p)
        if t > 1
            addcounter = t
            elemcounter = 0
            while addcounter > 0 && elemcounter < n
                elemcounter = elemcounter + 1
                current_element = ind[elemcounter]
                found = false
                for i = 1:t * (iter - 1)
                    if current_element == ind_hist[i]
                        found = true
                        break
                    end
                end
                if !found
                    addcounter = addcounter - 1
                    for i = 1:current_element - 1
                        X[i,t-addcounter] = 0
                    end
                    X[current_element,t-addcounter] = 1
                    for i = current_element + 1:n
                        X[i,t-addcounter] = 0
                    end
                    ind_hist[iter * t - addcounter] = current_element
                else
                    if elemcounter == t && addcounter == t
                        break
                    end
                end
            end
        else
            ind_hist[1:t] = ind[1:t]
            for j = 1:t
                for i = 1:ind[j] - 1
                    X[i,j] = 0
                end
                X[ind[j],j] = 1
                for i = ind[j] + 1:n
                    X[i,j] = 0
                end
            end
        end
    end
    return est
end


text:
tests = ["test-pval-adjustment",
         "test-pval-pi0-adjustment",
         "test-pi0-estimators",
         "test-grenander",
         "test-utils",
         "test-model",
         "test-doc"
         ]

for t in tests
    test_file = "$(t).jl"
    println(" * $(t) ...")
    include(test_file)
end


text:
"""
Base type for all problem descriptors. It is expected to has the following fields:
- `n`: the size of the problem.
- `saneComputationLimit`: maximum size of the problem for which computing results for
    adversarial lower bound, recoverable lower bound, selection lower bound or
    lagrangian lower bound makes sense.
- `equalCardinalityProperty`: specifies if problem possess equal cardinality property.
- `cardinality`: cardinality of the problem if any.
"""
abstract type ProblemDescriptor end

"""
    getProblemSize(pd::ProblemDescriptor)

Returns the size of the problem.

# Arguments
- `pd`: an instance of [`ProblemDescriptor`](@ref)
"""
getProblemSize(pd::ProblemDescriptor) = pd.n

"""
    getSaneComputationLimit(pd::ProblemDescriptor)

Returns maximum size of the problem for which computing results for adversarial
lower bound, recoverable lower bound, selection lower bound or lagrangian lower bound makes sense.

# Arguments
- `pd`: an instance of [`ProblemDescriptor`](@ref)
"""
getSaneComputationLimit(pd::ProblemDescriptor) = pd.saneComputationLimit

"""
    hasEqualCardinalityProperty(pd::ProblemDescriptor)

Returns whether the problem possess equal cardinality property.

# Arguments
- `pd`: an instance of [`ProblemDescriptor`](@ref)
"""
hasEqualCardinalityProperty(pd::ProblemDescriptor) = pd.equalCardinalityProperty

"""
    getCardinality(pd::ProblemDescriptor)

Returns cardinality of the problem if any.

# Arguments
- `pd`: an instance of [`ProblemDescriptor`](@ref)
"""
getCardinality(pd::ProblemDescriptor) = pd.cardinality

"""
`KnapsackProblemDescriptor` is an implementation [`ProblemDescriptor`](@ref) for
minimum knapsack problem.
"""
struct KnapsackProblemDescriptor <: ProblemDescriptor
    equalCardinalityProperty::Bool
    cardinality::Int
    n::Int
    saneComputationLimit::Int

    KnapsackProblemDescriptor(n::Integer) = new(false, -1, n, 400)
end

"""
`AssignmentProblemDescriptor` is an implementation [`ProblemDescriptor`](@ref) for
minimum assignment problem.
"""
struct AssignmentProblemDescriptor <: ProblemDescriptor
    equalCardinalityProperty::Bool
    cardinality::Int
    n::Int
    saneComputationLimit::Int

    AssignmentProblemDescriptor(n::Integer) = new(true, n, n, 25)
end


text:
"""
    count_nucleotides(strand)

The frequency of each nucleotide within `strand` as a dictionary.

Invalid strands raise a `DomainError`.

"""
function count_nucleotides(strand)
    count = Dict(
        'A' => 0, 'G' => 0,
        'C' => 0, 'T' => 0
    )
    for x in strand
        try
            count[x] += 1
        catch e
            throw(DomainError(x, "Invalid nucleotide"))
        end
    end

    return count

end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions