text:
using FixedPointNumbers

@testset "Parse" begin
    r8(x) = reinterpret(N0f8, x)

    # Color parsing
    # named-color
    redN0f8 = parse(Colorant, "red")
    @test colorant"red" === redN0f8
    @test redN0f8 === RGB{N0f8}(1,0,0)
    @test parse(RGB{Float64}, "red") === RGB{Float64}(1,0,0)
    @test isa(parse(HSV, "blue"), HSV)
    @test_throws ErrorException parse(Colorant, "p ink")
    @test parse(Colorant, "transparent") === RGBA{N0f8}(0,0,0,0)
    @test parse(Colorant, "\nSeaGreen ") === RGB{N0f8}(r8(0x2E),r8(0x8B),r8(0x57))
    @test parse(Colorant, "sea GREEN") === colorant"seagreen"

    # hex-color
    @test parse(Colorant, "#D0FF58") === RGB(r8(0xD0),r8(0xFF),r8(0x58))
    @test parse(Colorant, "0xd0ff58") === RGB(r8(0xD0),r8(0xFF),r8(0x58))
    @test parse(Colorant, "#FB0") === RGB(r8(0xFF),r8(0xBB),r8(0x00))
    @test parse(Colorant, "#FB0A") === RGBA(r8(0xFF),r8(0xBB),r8(0x00),r8(0xAA))
    @test parse(Colorant, "0xFB0A") === ARGB(r8(0xBB),r8(0x00),r8(0xAA),r8(0xFF))
    @test parse(Colorant, "#FFBB00AA") === RGBA(r8(0xFF),r8(0xBB),r8(0x00),r8(0xAA))
    @test parse(Colorant, "0xFFBB00AA") === ARGB(r8(0xBB),r8(0x00),r8(0xAA),r8(0xFF))
    @test_throws ErrorException parse(Colorant, "#BAD05")
    @test_throws ErrorException parse(Colorant, "#BAD0007")
    @test_throws ErrorException parse(Colorant, "#BAD000009")

    # rgb()
    @test parse(Colorant, "rgb(55,217,127)")      === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x7f))
    @test colorant"    rgb(    55, 217, 127 )   " === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x7f))
    @test parse(Colorant, "rgb(22%,85%,50%)")     === RGB{N0f8}(r8(0x38),r8(0xd9),r8(0x80))
    @test parse(Colorant, "rgba(55,217,127,0.5)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5)
    @test parse(Colorant, "rgb( 55,217,127,50%)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5) # CSS Color Module Level 4
    @test parse(Colorant, "rgb( 55 217 127 /.5)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5) # CSS Color Module Level 4
    @test parse(Colorant, "rgb(55, 85%, 50%)")    === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x80)) # this is invalid according to CSS spec.
    @test_throws ErrorException parse(Colorant, "rgb(21.6%,85%,50%)") # this is valid but not supported

    # hsl()
    @test parse(Colorant, "hsl(120, 100%, 50%)")  === HSL{Float32}(120,1.0,.5)
    @test colorant"    hsl(    120, 100%, 50% ) " === HSL{Float32}(120,1.0,.5)
    @test parse(RGB{N0f8},"hsl(120, 100%, 50%)")  === convert(RGB{N0f8}, HSL{Float32}(120,1.0,.5))
    @test_throws ErrorException  parse(Colorant, "hsl(120, 100, 50)")
    @test_throws ErrorException  parse(Colorant, "hsl(120%,100%,50%)")
    @test parse(Colorant, "hsla(120,50%,7%, .6)") === HSLA{Float32}(120,.5,.07,.6)
    @test parse(Colorant, "hsl( 120,50%,7%,60%)") === HSLA{Float32}(120,.5,.07,.6) # CSS Color Module Level 4
    @test parse(Colorant, "hsl( 120 50% 7% / 1)") === HSLA{Float32}(120,.5,.07, 1) # CSS Color Module Level 4

    @test parse(Colorant, :red) === colorant"red"
    @test_deprecated parse(Colorant, colorant"red") === colorant"red"

end


text:
module Bennu

using Adapt
using ArrayInterface
using CUDA
using CUDAKernels
using FillArrays
using GPUArrays
using KernelAbstractions
using LazyArrays
using LinearAlgebra
using LoopVectorization
using SparseArrays
using StaticArrays
using StaticArrays: tuple_prod, tuple_length, size_to_tuple
using StructArrays
using Tullio
using WriteVTK

export LobattoCell, NodalGrid, CartesianOrdering, HilbertOrdering

export spectralderivative, spectralinterpolation, legendregauss,
       legendregausslobatto, partition, quantize,
       hilbertcode, hilbertindices, hilbertperm,
       floattype, arraytype, points, mass, facemass, derivatives,
       derivatives_1d, degrees, celltype, referencecell, vertices,
       connectivity, points, brickgrid, toequallyspaced,
       cubedspheregrid, min_node_distance

export fieldarray

export components

export metrics, facemetrics

export faces, faceindices, boundaryfaces, number_of_faces, faceviews

include("permutations.jl")

include("arrays.jl")
include("cells.jl")
include("gridgenerators.jl")
include("grids.jl")
include("kernelabstractions.jl")
include("kroneckeroperators.jl")
include("operators.jl")
include("partitions.jl")
include("quadratures.jl")
include("sparsearrays.jl")
include("structarrays.jl")
include("tuples.jl")

end

ext: jl


text:
function test_set_met()
    model = MetNets.toy_model()

    met = MetNets.Met("new_met"; b = rand())
    for rxn in model.rxns
        push!(met.rxns, rxn)
        push!(met.S, rand([0,-1, 1]))
    end

    meti = 1
    MetNets.set_met!(model, meti, met)

    # Tests
    @test model.mets[meti] == met.id
    @test model.b[meti] == met.b

    for (rxn, s) in zip(met.rxns, met.S)
        rxni = MetNets.rxnindex(model, rxn)
        @test model.S[meti, rxni] == s
    end

end
test_set_met()


text:
function split_inds(inds::Vector{Int}, percentage::Float64)

    if !(0.0 < percentage < 1.0)
        throw("percentage must be in the range (0,1)")
    end

    N = length(inds)
    splitindex = round(Int, percentage*N)
    copy_inds = copy(inds)
    shuffle!(copy_inds)
    return copy_inds[splitindex+1:N], copy_inds[1:splitindex]
end

function _convertNH(node::treecausation.NodeMeta{S}, indX::Vector{Int}) where {S}
    if node.is_leaf
        return LeafCausalNH(indX[node.region])
    else
        left = _convertNH(node.l, indX)
        right = _convertNH(node.r, indX)
        return NodeCausalNH{S}(node.feature, node.threshold, left, right)
    end
end

function _convertH(node::treecausation.NodeMeta{S}, indX::Vector{Int}) where {S}
    if node.is_leaf
        return LeafCausalH(indX[node.region], Vector{Int}())
    else
        left = _convertH(node.l, indX)
        right = _convertH(node.r, indX)
        return NodeCausalH{S}(node.feature, node.threshold, left, right)
    end
end

function _fill!(
        node     :: LeafOrNodeCausalH{S},
        ind      :: Int,
        X        :: Vector{S}) where {S}

    if is_leaf(node)
        append!(node.inds_pred, ind)
    else
        if X[node.featid] < node.featval
            _fill!(node.left, ind, X)
        else
            _fill!(node.right, ind, X)
        end
    end
end

function fill_tree(
        node     :: treecausation.NodeMeta{S},
        indspred :: AbstractVector{Int},
        indX     :: Vector{Int},
        X        :: AbstractMatrix{S}) where {S}

    nodeH = _convertH(node, indX)
    n_samples = length(indspred)
    for i in 1:n_samples
        ind = indspred[i]
        X_obs = X[ind, :]
        _fill!(nodeH, ind, X_obs)
    end
    return nodeH
end

function build_tree(
        honest             :: Bool,
        indsbuild          :: AbstractVector{Int},
        indspred           :: Union{Nothing, AbstractVector{Int}},
        labels             :: AbstractVector{T},
        treatment          :: AbstractVector{Int},
        features           :: AbstractMatrix{S},
        const_mtry         :: Bool,
        m_pois              = -1,
        max_depth           = -1,
        min_samples_leaf    = 5,
        min_samples_split   = 2;
        rng                 = Random.GLOBAL_RNG) where {S, T <: Float64}

    if max_depth == -1
        max_depth = typemax(Int)
    end

    if m_pois == -1 && !const_mtry
        p = size(features, 2)
        m_pois = floor(Int, min(sqrt(p)+20, p))
    end

    if m_pois == -1 && const_mtry
        p = size(features, 2)
        m_pois = floor(Int, sqrt(p))
    end

    rng = mk_rng(rng)::Random.AbstractRNG

    t = treecausation.fit(
        X                   = features,
        Y                   = labels,
        W                   = treatment,
        indX                = indsbuild,
        const_mtry          = const_mtry,
        m_pois              = Int(m_pois),
        max_depth           = Int(max_depth),
        min_samples_leaf    = Int(min_samples_leaf),
        min_samples_split   = Int(min_samples_split),
        rng                 = rng)

    if honest
        return TreeCausalH{S}(fill_tree(t.root, indspred, t.inds, features), indsbuild, indspred, setdiff(collect(1:length(labels)), union(indsbuild, indspred)))
    else
        return TreeCausalNH{S}(_convertNH(t.root, t.inds), indsbuild, setdiff(collect(1:length(labels)), indsbuild))
    end
end

"""
Build a causal forest.

- if `centering=True` Y and W are centered else they stay unchanged
- if `bootstrap=True` we sample for each tree via bootstrap else we use subsampling
- if `honest=True` we use 2 samples one too build splits and the other one to fill leaves
    otherwise we use the whole sample for the two steps
- if `const_mtry=True` we use a constant mtry otherwise we use a random mtry following
    `min(max(Poisson(m_pois),1),number_of_features)`
"""
function build_forest(
    centering          :: Bool,
    bootstrap          :: Bool,
    honest             :: Bool,
    labels             :: AbstractVector{T},
    treatment          :: AbstractVector{Int},
    features           :: AbstractMatrix{S},
    const_mtry         :: Bool,
    m_pois              = -1,
    n_trees             = 10,
    partial_sampling    = 0.7,
    honest_proportion   = 0.5,
    max_depth           = -1,
    min_samples_leaf    = 5,
    min_samples_split   = 10,
    n_trees_centering   = 100;
    rng                 = Random.GLOBAL_RNG) where {S, T <: Float64}

    if n_trees < 1
        throw("the number of trees must be >= 1")
    end
    if !(0.0 < partial_sampling <= 1.0)&&!bootstrap
        throw("partial_sampling must be in the range (0,1]")
    end

    n_tot_samples = length(labels)

    if bootstrap
        n_samples = length(labels)
    else
        n_samples = floor(Int, partial_sampling * length(labels))
    end

    if honest
        forest = Vector{TreeCausalH{S}}(undef, n_trees)
    else
        forest = Vector{TreeCausalNH{S}}(undef, n_trees)
    end

    if centering
        model_Y = build_forest_oob(labels, features, -1, n_trees_centering)
        model_T = build_forest_oob(treatment, features, -1, n_trees_centering)
        Y_center = labels - apply_forest_oob(model_Y)
        T_center = treatment - apply_forest_oob(model_T)
        Y_vec = Y_center
        T_vec = T_center
    else
        model_Y = nothing
        model_T = nothing
        Y_center = nothing
        T_center = nothing
        Y_vec = labels
        T_vec = treatment
    end

    if rng isa Random.AbstractRNG
        if bootstrap
            Threads.@threads for i in 1:n_trees
                if honest
                    inds = collect(1:n_tot_samples)
                    inds1, inds2 = split_inds(inds, honest_proportion)
                    indsbuild = rand(rng, inds1, length(inds1))
                    indspred = rand(rng, inds2, length(inds2))
                else
                    indsbuild = rand(rng, 1:n_tot_samples, n_samples)
                    indspred = nothing
                end
                forest[i] = build_tree(
                    honest,
                    indsbuild,
                    indspred,
                    Y_vec,
                    T_vec,
                    features,
                    const_mtry,
                    m_pois,
                    max_depth,
                    min_samples_leaf,
                    min_samples_split,
                    rng = rng)
            end
        else
            Threads.@threads for i in 1:n_trees
                if honest
                    inds = StatsBase.sample(rng, 1:n_tot_samples, n_samples; replace=false)
                    indsbuild, indspred = split_inds(inds, honest_proportion)
                else
                    indsbuild = StatsBase.sample(rng, 1:n_tot_samples, n_samples; replace=false)
                    indspred = nothing
                end
                forest[i] = build_tree(
                    honest,
                    indsbuild,
                    indspred,
                    Y_vec,
                    T_vec,
                    features,
                    const_mtry,
                    m_pois,
                    max_depth,
                    min_samples_leaf,
                    min_samples_split,
                    rng = rng)
            end
        end
    elseif rng isa Integer # each thread gets its own seeded rng
        if bootstrap
            Threads.@threads for i in 1:n_trees
                Random.seed!(rng + i)
                if honest
                    inds = collect(1:n_tot_samples)
                    inds1, inds2 = split_inds(inds, honest_proportion)
                    indsbuild = rand(inds1, length(inds1))
                    indspred = rand(inds2, length(inds2))
                else
                    indsbuild = rand(1:n_tot_samples, n_samples)
                    indspred = nothing
                end
                forest[i] = build_tree(
                    honest,
                    indsbuild,
                    indspred,
                    Y_vec,
                    T_vec,
                    features,
                    const_mtry,
                    m_pois,
                    max_depth,
                    min_samples_leaf,
                    min_samples_split)
            end
        else
            Threads.@threads for i in 1:n_trees
                Random.seed!(rng + i)
                if honest
                    inds = StatsBase.sample(1:n_tot_samples, n_samples; replace=false)
                    indsbuild, indspred = split_inds(inds, honest_proportion)
                else
                    indsbuild = StatsBase.sample(1:n_tot_samples, n_samples; replace=false)
                    indspred = nothing
                end
                forest[i] = build_tree(
                    honest,
                    indsbuild,
                    indspred,
                    Y_vec,
                    T_vec,
                    features,
                    const_mtry,
                    m_pois,
                    max_depth,
                    min_samples_leaf,
                    min_samples_split)
            end
        end
    else
        throw("rng must of be type Integer or Random.AbstractRNG")
    end

    return EnsembleCausal{S}(forest, centering, bootstrap, honest, features, labels, treatment, model_Y, model_T, Y_center, T_center)
end

apply_treeH(leaf :: LeafCausalH, x :: AbstractVector{S}) where {S} = leaf.inds_pred
apply_treeNH(leaf :: LeafCausalNH, x :: AbstractVector{S}) where {S} = leaf.inds_build

function apply_treeH(tree :: NodeCausalH{S}, x :: AbstractVector{S}) where {S}
    if tree.featid == 0
        return apply_treeH(tree.left, x)
    elseif x[tree.featid] < tree.featval
        return apply_treeH(tree.left, x)
    else
        return apply_treeH(tree.right, x)
    end
end

function apply_treeNH(tree :: NodeCausalNH{S}, x :: AbstractVector{S}) where {S}
    if tree.featid == 0
        return apply_treeNH(tree.left, x)
    elseif x[tree.featid] < tree.featval
        return apply_treeNH(tree.left, x)
    else
        return apply_treeNH(tree.right, x)
    end
end

function neighbours(
    forest :: EnsembleCausal{S},
    b      :: Int,
    x      :: AbstractVector{S}
    ) where {S}
    if forest.honest
        return apply_treeH(forest.trees[b].tree, x)
    else
        return apply_treeNH(forest.trees[b].tree, x)
    end
end

function apply_forest(
    forest :: EnsembleCausal{S},
    x      :: AbstractVector{S}
    ) where {S}

    centering = forest.centering
    n_samples = length(forest.Y)
    alpha = zeros(n_samples)
    n_trees = length(forest)
    for b in 1:n_trees
        N = neighbours(forest, b, x)
        l = length(N)
        for e in N
            alpha[e] += 1/l
        end
    end
    alpha = alpha/n_trees
    T_alpha = 0
    Y_alpha = 0
    if centering
        Y_vec =  forest.Y_center
        T_vec =  forest.T_center
    else
        Y_vec = forest.Y
        T_vec = forest.T
    end
    for i in 1:n_samples
        T_alpha += alpha[i]*T_vec[i]
        Y_alpha += alpha[i]*Y_vec[i]
    end
    num = 0
    denom = 0
    for i in 1:n_samples
        diffT = T_vec[i]-T_alpha
        num += alpha[i]*diffT*(Y_vec[i]-Y_alpha)
        denom += alpha[i]*diffT^2
    end
    return num/denom
end

"""
Get the causal effect for each row in x given a causal forest
"""
function apply_forest(
    forest :: EnsembleCausal{S},
    x      :: AbstractMatrix{S}
    ) where {S}

    N = size(x, 1)
    predictions = Array{Float64}(undef, N)
    for i in 1:N
        predictions[i] = apply_forest(forest, x[i, :])
    end
    return predictions
end


text:
# On vérifie que le constructeur par défaut fonctionne
ls_at_t = LSAtT(0.0)

@test ls_at_t.x   == 0.0
@test ls_at_t.ht  == nothing
@test ls_at_t.gt  == nothing
@test ls_at_t.h₀  == nothing
@test ls_at_t.g₀  == nothing

@test ls_at_t.current_time == nothing
@test ls_at_t.current_score == nothing

# On test la fonction update!(...)
update!(ls_at_t, x = 1.0, ht = 1.0, gt = 1.0, h₀ = 1.0)
update!(ls_at_t, g₀ = 1.0, current_time = 0.0, current_score = 0.0)

@test ls_at_t.x == 1.0
@test ls_at_t.ht == 1.0
@test ls_at_t.gt == 1.0
@test ls_at_t.h₀ == 1.0
@test ls_at_t.g₀ == 1.0
@test ls_at_t.current_time == 0.0
@test ls_at_t.current_score == 0.0

# on vérifie que la fonction copy fonctionne
ls_at_t_2 = copy(ls_at_t)

@test ls_at_t_2.x == 1.0
@test ls_at_t_2.ht == 1.0
@test ls_at_t_2.gt == 1.0
@test ls_at_t_2.h₀ == 1.0
@test ls_at_t_2.g₀ == 1.0
@test ls_at_t_2.current_time == 0.0
@test ls_at_t_2.current_score == 0.0

ls_64 = LSAtT(0.0)
update!(ls_64, x = 1.0, ht = 1.0, gt = 1.0, h₀ = 1.0)

# ls_32 = convert_ls(Float32, ls_64)
# @test typeof(ls_32.x) == Float32

reinit!(ls_64)
@test ls_64.x == 1.0
@test ls_64.ht == nothing
@test ls_64.current_time == nothing


text:
using ZChop
using Base.Test

z1 = complex(1,1e-15)
zz1 = complex(1,0.0)
z2 = complex(1e-15,1)
zz2 = complex(0.0,1)
za = [z1,z2]
zza = [zz1,zz2]

@test zchop(1.0) == 1.0
@test zchop(1e-15) == 0.0
@test zchop(z1) == zz1
@test zchop(z2) == zz2
@test zchop(za) == zza
t = (za,za)
@test zchop( t ) == (zza,zza)
@test zchop( [za za] ) == [zza zza]
@test zchop(BigFloat[1e-15,1]) == BigFloat[0,1]
@test zchop(1//10^15) == 0//1
@test zchop([1e-16,"cat"]) == [ 0.0, "cat"]
r = r"cat"
@test zchop(r) == r
@test zchop(1) == 1
@test zchop(1,3) == 0
@test zchop(Int) == Int
@test zchop([true,false]) == [true,false]
@test zchop('x') == 'x'
@test zchop(STDOUT) == STDOUT
@test zchop(Base.pi) == float(Base.pi)
@test zchop(Base.pi,4) == 0.0
@test zchop(BigInt(1)) == BigInt(1)
@test zchop(BigInt(1),3) == BigInt(0)
m = [[1.0,2.0,1e-16] [3.0,4.0,5.0]]
@test zchop(m) == [[1.0,2.0,0.0] [3.0,4.0,5.0]]
@test zchop(:a) == :a
@test zchop(:(1+1)) == :(1+1)
@test zchop(:(1+1e-16)) == :(1+0.0)
@test zchop(:( "cat" + [1e-16,2.0] )) == :( "cat" + [0.0,2.0] )

let b = collect(1:100) * 1e-15, c = zchop(b)
    zchop!(b)
    @test c == b
end

# Tuples
@test zchop((1e-15,1,1)) == (0.0,1,1)

@test_throws MethodError zchop()


text:
""" Deserialize abstract wiring diagram from yFiles.

Reads a wiring diagram from the GraphML dialect used by yEd and yFiles. Unlike
the GraphML spec, the yEd data model does not explicitly include ports:

- https://yed.yworks.com/support/qa/102/
- https://yed.yworks.com/support/qa/2531/

We infer the ports of boxes and their order from the geometry of the diagram.
Thus, this module has the nature of a hack. While it may be useful for
interactive and exploratory work, it should not be used in a production system.
"""
module YFilesWiringDiagrams
export read_yfiles_diagram, parse_yfiles_diagram

using Compat
using LightXML
using LightGraphs, MetaGraphs

using ..WiringDiagramCore
using ..GraphMLWiringDiagrams: convert_from_graphml_data, parse_graphml_metagraph
import ..GraphMLWiringDiagrams: parse_graphml_data_value

# Data types
############

struct BoxLayout
  box::Int
  input_coord_map::Dict  # Map from coordinates to input ports
  output_coord_map::Dict # Map from coordinates to output ports
end

# Deserialization
#################

""" Read a wiring diagram from a GraphML file created by yEd and yFiles.
"""
function read_yfiles_diagram(BoxValue::Type, WireValue::Type, filename::String; kw...)
  parse_yfiles_diagram(BoxValue, WireValue, LightXML.parse_file(filename); kw...)
end

""" Parse a wiring diagram from a GraphML string or XML doc created by yFiles.
"""
function parse_yfiles_diagram(BoxValue::Type, WireValue::Type, s::AbstractString; kw...)
  parse_yfiles_diagram(BoxValue, WireValue, LightXML.parse_string(s); kw...)
end
function parse_yfiles_diagram(BoxValue::Type, WireValue::Type, xdoc::XMLDocument;
    direction::Symbol=:vertical, keep_labels::Bool=true)::WiringDiagram
  @assert direction in (:horizontal, :vertical)

  # Clean up GraphML keys before reading.
  xroot = root(xdoc)
  for xkey in xroot["key"]
    yfiles_type = attribute(xkey, "yfiles.type")
    if !isnothing(yfiles_type)
      if !has_attribute(xkey, "attr.name")
        set_attribute(xkey, "attr.name", yfiles_type)
      end
      if !has_attribute(xkey, "attr.type")
        set_attribute(xkey, "attr.type", "yfiles_$yfiles_type")
      end
    end
  end

  # Read the diagram's underlying graph as an attributed graph.
  graph = parse_graphml_metagraph(xdoc, directed=true, multigraph=true)

  # Extract needed information from yFiles' "nodegraphics" and "edgegraphics"
  # and discard the rest. Keep custom data properties, except the blank
  # "description" property inserted by yEd.
  for v in 1:nv(graph)
    v_data = props(graph, v)
    if haskey(v_data, :description) && isempty(v_data[:description])
      delete!(v_data, :description)
    end
    node_graphics = pop!(v_data, :nodegraphics)
    if keep_labels & haskey(node_graphics, :label)
      v_data[:label] = node_graphics[:label]
    end
  end
  for edge in edges(graph)
    for wire_data in get_prop(graph, edge, :edges)
      if haskey(wire_data, :description) && isempty(wire_data[:description])
        delete!(wire_data, :description)
      end
      edge_graphics = pop!(wire_data, :edgegraphics)
      wire_data[:source_coord] = round(Int,
        edge_graphics[direction == :vertical ? :source_x : :source_y])
      wire_data[:target_coord] = round(Int,
        edge_graphics[direction == :vertical ? :target_x : :target_y])
      if keep_labels & haskey(edge_graphics, :label)
        wire_data[:label] = edge_graphics[:label]
      end
    end
  end

  # Add boxes and their ports to diagram, including the outer box.
  diagram = WiringDiagram([], [])
  boxes = BoxLayout[]
  for v in 1:nv(graph)
    box_layout = if pop!(props(graph, v), :input, false)
      # Special case: diagram inputs.
      ports, coord_map = infer_output_ports(graph, v)
      diagram.input_ports = ports
      BoxLayout(input_id(diagram), Dict(), coord_map)
    elseif pop!(props(graph, v), :output, false)
      # Special case: diagram outputs.
      ports, coord_map = infer_input_ports(graph, v)
      diagram.output_ports = ports
      BoxLayout(output_id(diagram), coord_map, Dict())
    else
      # Generic case: a box.
      input_ports, input_coord_map = infer_input_ports(graph, v)
      output_ports, output_coord_map = infer_output_ports(graph, v)
      value = convert_from_graphml_data(BoxValue, props(graph, v))
      box_id = add_box!(diagram, Box(value, input_ports, output_ports))
      BoxLayout(box_id, input_coord_map, output_coord_map)
    end
    push!(boxes, box_layout)
  end

  # Add wires to diagram.
  for edge in edges(graph)
    source, target = boxes[src(edge)], boxes[dst(edge)]
    for wire_data in get_prop(graph, edge, :edges)
      source_port = source.output_coord_map[pop!(wire_data, :source_coord)]
      target_port = target.input_coord_map[pop!(wire_data, :target_coord)]
      value = convert_from_graphml_data(WireValue, wire_data)
      add_wire!(diagram, Wire(value,
        (source.box, source_port) => (target.box, target_port)))
    end
  end

  diagram
end

function infer_input_ports(graph::MetaDiGraph, v::Int)
  in_edges = (get_prop(graph, u, v, :edges) for u in inneighbors(graph, v))
  in_edges = reduce(vcat, in_edges; init=[])
  in_coords = [ edge[:target_coord] for edge in in_edges ]
  infer_ports_from_coordinates(in_coords)
end
function infer_output_ports(graph::MetaDiGraph, v::Int)
  out_edges = (get_prop(graph, v, u, :edges) for u in outneighbors(graph, v))
  out_edges = reduce(vcat, out_edges; init=[])
  out_coords = [ edge[:source_coord] for edge in out_edges ]
  infer_ports_from_coordinates(out_coords)
end

function infer_ports_from_coordinates(coords::Vector{T}) where T
  unique_coords = sort(unique(coords))
  ports = repeat([nothing], length(unique_coords))
  coord_map = Dict{T,Int}(x => i for (i, x) in enumerate(unique_coords))
  (ports, coord_map)
end

function parse_graphml_data_value(::Type{Val{:yfiles_nodegraphics}}, xdata::XMLElement)
  ynode = first(child_elements(xdata)) # e.g., ShapeNode
  ygeom = find_element(ynode, "Geometry")
  data = Dict{Symbol,Any}(
    :x => float_attribute(ygeom, "x"),
    :y => float_attribute(ygeom, "y"),
    :width => float_attribute(ygeom, "width"),
    :height => float_attribute(ygeom, "height"),
  )
  ylabel = find_element(ynode, "NodeLabel")
  if !isnothing(ylabel)
    data[:label] = content(ylabel)
  end
  data
end

function parse_graphml_data_value(::Type{Val{:yfiles_edgegraphics}}, xdata::XMLElement)
  yedge = first(child_elements(xdata)) # e.g., PolyLineEdge
  ypath = find_element(yedge, "Path")
  data = Dict{Symbol,Any}(
    :source_x => float_attribute(ypath, "sx"),
    :source_y => float_attribute(ypath, "sy"),
    :target_x => float_attribute(ypath, "tx"),
    :target_y => float_attribute(ypath, "ty"),
  )
  ylabel = find_element(yedge, "EdgeLabel")
  if !isnothing(ylabel)
    data[:label] = content(ylabel)
  end
  data
end

float_attribute(xelem::XMLElement, name::AbstractString) =
  parse(Float64, attribute(xelem, name, required=true))

end


text:
@testset "Ideals for orders over function fields" begin

  k = GF(7)
  kx, x = RationalFunctionField(k, "x")
  kt = parent(numerator(x))
  ky, y = PolynomialRing(kx, "y")
  F, a = function_field(y^2+x)
  O = integral_closure(kt, F)

  h = O.R(x^2+1)

  f1 = y+5*x+2
  f2 = y+2*x+5

  I = FfOrdIdl(O,h, O(y+5*x+2))
  J = FfOrdIdl(O,h, O(y+2*x+5))
  K2 = FfOrdIdl(O,h)
  @test K2 == I*J

  A = I^3*J
  L = @inferred factor(A)
  G = [(fac,e) for (fac,e) in L]
  @test (I,3) in G
  @test (J,1) in G
  @test length(G)==2

  k = QQ
  kx, x = RationalFunctionField(k, "x")
  kt = parent(numerator(x))
  ky, y = PolynomialRing(kx, "y")
  F, a = function_field(y^2+x*y+x^3+y^3)
  O = integral_closure(kt, F)

  @test @inferred index(O) == x^2 - 1//3*x
  h = O.R(x)
  L = prime_decomposition(O, h)
  @test prod([f[1]^f[2] for f in L]) == FfOrdIdl(O, h)
end


text:
module NormalSplines

#### Inteface deinition
export interpolate, evaluate
export RK_W1, RK_W2, RK_W3, RK_H1, RK_H2, RK_H3
export get_cond, get_epsilon
####

using LinearAlgebra

abstract type ReproducingKernel end
abstract type ReproducingKernel_1 <: ReproducingKernel   end
abstract type ReproducingKernel_2 <: ReproducingKernel_1 end
abstract type ReproducingKernel_3 <: ReproducingKernel_2 end

include("./ReproducingKernels.jl")
include("./GramMatrix.jl")
include("./Interpolate.jl")

"""
`interpolate(knots::Vector{T}, values::Vector{T}, kernel :: RK, factorize = true) where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Create an interpolating normal spline by values of function in knots.

# Arguments
- `knots::Vector{T}`: locations of given function values.
- `values::Vector{T}`: function values in `knots`.
- `kernel::RK`: reproducing kernel of Hilbert space the normal spline is constructed in.
                It must be a struct object of the following type:
                  `RK_W1` or `RK_H1` if the spline is constructing as a continuous function,
                  `RK_W2` or `RK_H2` if the spline is constructing as a differentiable function,
                  `RK_W3` or `RK_H3` if the spline is constructing as a twice differentiable function.
- `factorize::Bool`: defines if it is necessary to compute the Gram matrix factorization.
                  Must be set to `true` if the spline is constructing first time.
                  Can be set to `false` if the spline is constructing with the same knots
                  as it was done first time but with different function values in knots.

Returns : nothing.
"""
function interpolate(knots::Vector{T},
                     values::Vector{T},
                     kernel :: RK,
                     factorize = true) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
    if factorize
        _prepare_spline(kernel, knots)
    end
    _construct_spline(values)
    return
end

"""
`interpolate(knots::Vector{T}, values::Vector{T},
             d_knots::Vector{T}, d_values::Vector{T},
             kernel :: RK, derivative = 1, factorize = true) where {T <: AbstractFloat, RK <: ReproducingKernel_2}`

Create an interpolating normal spline by values of function
and it's first or second derivatives in knots.

# Arguments
- `knots::Vector{T}`: locations of given function values.
- `values::Vector{T}`: function values in `knots`.
- `d_knots::Vector{T}`: locations of given function first derivative values (if `derivative`=`1`)
                        or second derivative values (if `derivative`=`2`).
- `d_values::Vector{T}`: values of function first derivative (if derivative=`1`) or
                         values of function second derivative (if derivative=`2`) in `d_knots`.
- `kernel::RK`: reproducing kernel of Hilbert space the normal spline is constructed in.
                In a case when `derivative`=`1`, `kernel` must be a struct object of the following type:
                     `RK_W2` or `RK_H2` if the spline is constructing as a differentiable function,
                     `RK_W3` or `RK_H3` if the spline is constructing as a twice differentiable function.
                In a case when `derivative`=`2`, `kernel` must be a struct object of the following type:
                     `RK_W3` or `RK_H3` - the spline is constructing as a twice differentiable function,
- `derivative::Int`: `1` - `d_knots` and `d_values` provide data of function first derivative,
                     `2` - `d_knots` and `d_values` provide data of function second derivative.
- `factorize::Bool`: defines if it is necessary to compute the Gram matrix factorization.
                     Must be set to `true` if the spline is constructing first time.
                     Can be set to `false` if the spline is constructing with the same knots
                     as it was done first time but with different function values in knots.
Returns : nothing.
"""
function interpolate(knots::Vector{T},
                     values::Vector{T},
                     d_knots::Vector{T},
                     d_values::Vector{T},
                     kernel :: RK,
                     derivative = 1,
                     factorize = true) where {T <: AbstractFloat, RK <: ReproducingKernel_2}
    if derivative == 1
        if factorize
            _prepare_spline_d1(kernel, knots, d_knots)
        end
        _construct_spline_d1(values, d_values)
    elseif derivative == 2
        if isa(kernel, RK_H2) || isa(kernel, RK_W2)
           Throw(ArgumentError("kernel: Incorrect parameter type (must be RK_H3 or RK_W3)."))
        end
        if factorize
            _prepare_spline_d2(kernel, knots, d_knots)
        end
        _construct_spline_d2(values, d_values)
    else
        Throw(DomainError(derivative, "'derivative' parameter value must be 1 or 2."))
    end

    return
end

"""
`interpolate(knots::Vector{T}, values::Vector{T},
             d1_knots::Vector{T}, d1_values::Vector{T},
             d2_knots::Vector{T}, d2_values::Vector{T},
             kernel :: RK, factorize = true) where {T <: AbstractFloat, RK <: ReproducingKernel_3}`

Create an interpolating normal spline by values of function
and it's first and second derivatives in knots.

# Arguments
- `knots::Vector{T}`: locations of given function values.
- `values::Vector{T}`: function values in `knots`.
- `d1_knots::Vector{T}`: locations of given function derivative values.
- `d1_values::Vector{T}`: values of function derivative in `d1_knots`.
- `d2_knots::Vector{T}`: locations of given function second derivative values.
- `d2_values::Vector{T}`: values of function second derivative in `d2_knots`.
- `kernel::RK`: reproducing kernel of Hilbert space the normal spline is constructed in.
                It must be a struct object of the following type:
                  `RK_W3` or `RK_H3` - the spline is constructing as a twice differentiable function.
- `factorize::Bool`: defines if it is necessary to compute the Gram matrix factorization.
                  Must be set to `true` if the spline is constructing first time.
                  Can be set to `false` if the spline is constructing with the same knots
                  as it was done first time but with different function values in knots.

Returns : nothing.
"""
function interpolate(knots::Vector{T},
                     values::Vector{T},
                     d1_knots::Vector{T},
                     d1_values::Vector{T},
                     d2_knots::Vector{T},
                     d2_values::Vector{T},
                     kernel :: RK,
                     factorize = true) where {T <: AbstractFloat, RK <: ReproducingKernel_3}
    if factorize
        _prepare_spline_d12(kernel, knots, d1_knots, d2_knots)
    end
    _construct_spline_d12(values, d1_values, d2_values)
    return
end

"""
`evaluate(points::Vector{T}, derivative::Int = 0) where T <: AbstractFloat`

Evaluate normal spline and it's derivatives.

# Arguments
- `points::Vector{T}`: locations at which spline or its derivative are evaluating.
- `derivative::Int`: `0` - spline evaluation.
                     `1` - first derivative of spline evaluation.
                     `2` - second derivative of spline evaluation.

Returns: `Vector{T}` of spline values or its derivatives at the locations
          defined in `points`.
"""
function evaluate(points::Vector{T}, derivative::Int = 0) where T <: AbstractFloat
    if derivative == 0
        return _evaluate(points)
    elseif derivative == 1
        return _evaluate_d1(points)
    elseif derivative == 2
        return _evaluate_d2(points)
    else
        Throw(DomainError(derivative, "'derivative' parameter value must be 0, 1 or 2."))
    end
end

"""
`evaluate(point::T, derivative::Int = 0) where T <: AbstractFloat`

Evaluate normal spline and it's derivatives.

# Arguments
- `point::T`: location at which spline or its derivative are evaluating.
- `derivative::Int`: `0` - spline evaluation.
                     `1` - first derivative of spline evaluation.
                     `2` - second derivative of spline evaluation.

Returns: spline value or its derivative at the `point` location.
"""
function evaluate(point::T, derivative::Int = 0) where T <: AbstractFloat
    v_points = Vector{T}(undef, 1)
    v_points[1] = point
    return evaluate(v_points, derivative)[1]
end

"""
`get_cond()`

Estimate condition number of Gram matrix.
"""
function get_cond()
    return _get_cond()
end

"""
`get_epsilon()`

Get value of parameter `ε` of struct object `RK_H1`,
`RK_H2` or `RK_H3`.
"""
function get_epsilon()
    return _get_epsilon()
end

end # module


text:
using LCPSim
using Base.Test

@testset "LCPSim" begin
    include("planar_vs_dummy_link.jl")
    include("joint_limits.jl")
    include("contact_lqr.jl")
    include("linearized_mechanism.jl")
    include("warmstart.jl")
    include("inclined_brick.jl")
    include("acrobot.jl")
    include("../examples/box.jl")
    include("../examples/rimless_wheel.jl")
end


text:
"""
    center_gradient(rbm, ∂, λv, λh)

Given the gradient `∂` of `rbm`, returns the gradient of the equivalent centered RBM
with offsets `λv` and `λh`.
"""
function center_gradient(rbm::RBM, ∂::NamedTuple, λv::AbstractArray, λh::AbstractArray)
    ∂ = oftype(∂, center_gradient_v(rbm, ∂, λv))
    ∂ = oftype(∂, center_gradient_h(rbm, ∂, λh))
    return ∂
end

function center_gradient_v(rbm::RBM, ∂::NamedTuple, λv::AbstractArray)
    @assert size(rbm.visible) == size(λv)
    @assert size(∂.w) == size(rbm.w)
    Δh = grad2ave(rbm.hidden, ∂.hidden)
    ∂w = reshape(∂.w, length(rbm.visible), length(rbm.hidden)) + vec(λv) * vec(Δh)'
    return oftype(∂, (; ∂..., w = reshape(∂w, size(∂.w))))
end

function center_gradient_h(rbm::RBM, ∂::NamedTuple, λh::AbstractArray)
    @assert size(rbm.hidden) == size(λh)
    @assert size(∂.w) == size(rbm.w)
    Δv = grad2ave(rbm.visible, ∂.visible)
    ∂w = reshape(∂.w, length(rbm.visible), length(rbm.hidden)) + vec(Δv) * vec(λh)'
    return oftype(∂, (; ∂..., w = reshape(∂w, size(∂.w))))
end

"""
    uncenter_step(rbm, ∂, λv, λh)

Given parameter update step `∂` of a centered `rbm` with offsets `λv, λh`, returns
the corresponding gradient of the equivalent uncentered RBM.
"""
function uncenter_step(rbm::RBM, ∂::NamedTuple, λv::AbstractArray, λh::AbstractArray)
    @assert size(rbm.visible) == size(λv)
    @assert size(rbm.hidden) == size(λh)
    @assert size(∂.w) == size(rbm.w)
    ∂w = reshape(∂.w, length(rbm.visible), length(rbm.hidden))
    shift_v = reshape(∂w  * vec(λh), size(rbm.visible))
    shift_h = reshape(∂w' * vec(λv), size(rbm.hidden))
    ∂v = uncenter_step(rbm.visible, ∂.visible, shift_v)
    ∂h = uncenter_step(rbm.hidden,  ∂.hidden,  shift_h)
    return (visible = oftype(∂.visible, ∂v), hidden = oftype(∂.hidden, ∂h), w = ∂.w)
end

function centered_gradient(rbm::RBM, ∂::NamedTuple, λv::AbstractArray, λh::AbstractArray)
    ∂c = center_gradient(rbm, ∂, λv, λh)
    return uncenter_step(rbm, ∂c, λv, λh)
end

function uncenter_step(
    layer::Union{Binary,Spin,Potts,Gaussian,ReLU,pReLU,xReLU},
    ∂::NamedTuple, shift::AbstractArray
)
    @assert size(layer) == size(∂.θ) == size(shift)
    return (; ∂..., θ = ∂.θ - shift,)
end

function uncenter_step(layer::dReLU, ∂::NamedTuple, shift::AbstractArray)
    @assert size(layer) == size(∂.θp) == size(∂.θn) == size(shift)
    return (θp = ∂.θp - shift, θn = ∂.θn - shift, γp = ∂.γp, γn = ∂.γn)
end

# get moments from layer gradients, e.g. <v> = -derivative w.r.t. θ
grad2ave(::Union{Binary,Spin,Potts,Gaussian,ReLU,pReLU,xReLU}, ∂::NamedTuple) = -∂.θ
grad2ave(::dReLU, ∂::NamedTuple) = -(∂.θp + ∂.θn)

grad2var(::Union{Binary,Potts}, ∂::NamedTuple) = -∂.θ .* (1 .+ ∂.θ)
grad2var(::Spin, ∂::NamedTuple) = (1 .- ∂.θ) .* (1 .+ ∂.θ)
grad2var(l::Union{Gaussian,ReLU}, ∂::NamedTuple) = 2∂.γ .* sign.(l.γ) - ∂.θ.^2
grad2var(l::dReLU, ∂::NamedTuple) = 2 * (∂.γp .* sign.(l.γp) + ∂.γn .* sign.(l.γn)) - (∂.θp + ∂.θn).^2

function grad2var(l::pReLU, ∂::NamedTuple)
    abs_γ = abs.(l.γ)
    ∂absγ = ∂.γ .* sign.(l.γ)
    @. 2l.η/abs_γ * ((2l.Δ * ∂.Δ + l.η * ∂.η) * l.η - ∂.η - l.Δ * ∂.θ) + 2∂absγ * (1 + l.η^2) - ∂.θ^2
end

function grad2var(l::xReLU, ∂::NamedTuple)
    abs_γ = abs.(l.γ)
    ∂absγ = ∂.γ .* sign.(l.γ)
    ν = @. 2∂absγ - ∂.θ^2
    @. (ν * abs_γ - 2 * (∂.ξ + ∂.θ * l.Δ) * l.ξ + ((ν + 2∂absγ) * abs_γ + 4 * ∂.Δ * l.Δ) * l.ξ^2 - 4∂.ξ * l.ξ^3 + 2abs(l.ξ) * (ν * abs_γ - 3∂.ξ * l.ξ - ∂.θ * l.Δ * l.ξ)) / (abs_γ * (1 + abs(l.ξ))^2)
end


text:
using ECOS, SCS, Gurobi
using JuMP

mysolver = ECOSSolver()
#mysolver = SCSSolver()
#mysolver = GurobiSolver()

items  = [:Gold, :Silver, :Bronze]
values = Dict(:Gold => 5.0,  :Silver => 3.0,  :Bronze => 1.0)
weight = Dict(:Gold => 2.0,  :Silver => 1.5,  :Bronze => 0.3)

m = Model(solver=mysolver)
@defVar(m, 0 <= take_var[items] <= 1)  # Define a variable for each item
@setObjective(m, Max, sum{ values[item] * take_var[item], item in items})
@addConstraint(m, sum{ weight[item] * take_var[item], item in items} <= 3)
@time JuMP.solve(m)

println("Optimal objective: ",getObjectiveValue(m))
println("Optimal vars: ",getValue(take_var))

# take
# [  Gold] = 0.9999999680446406
# [Silver] = 0.46666670881026834
# [Bronze] = 0.9999999633898735

max_stars_repo_name: ToyotaResearchInstitute/rad-robot

max_issues_repo_name: ToyotaResearchInstitute/rad-robot

max_forks_repo_name: ToyotaResearchInstitute/rad-robot


text:
####
#### FermiOp
####

module FermiOps

import Random
import IsApprox

export FermiOp, FermiDefault
export I, NumberOp, Empty, Raise, Lower, Zero, count_raise_lower, is_raise_lower

import ..AbstractOps: _AbstractOp, op_symbols, AbstractOp, _show_op_plain, op_index
import ..AbstractFermiOp

struct FermiOp <: AbstractFermiOp
    ind::Int  # Int is often faster than UInt8
#    ind::UInt8
end

const FermiDefault = FermiOp

op_index(op::FermiOp) = op.ind

Base.one(op::FermiOp) = I
Base.one(::Type{FermiOp}) = I

Base.zero(op::FermiOp) = Zero
Base.zero(::Type{FermiOp}) = Zero

"""
    isless(op1::FermiOp, op2::FermiOp)

Canonical (lexical) order of `FermiOp` is I < N < (I - N) < Raise < Lower < Zero < (Z = N - E)
"""
Base.isless(op1::FermiOp, op2::FermiOp) = isless(op_index(op1), op_index(op2))

const _fermi_chars =
     ('I',   'N',    'E',   '+',   '-',  '0', 'Z', 'X')
#      0      1       2      3      4     5    6    7
const (I, NumberOp, Empty, Raise, Lower, Zero, Z,  NoOp) = FermiOp.(0:7)

const _fermi_string_chars = string.(_fermi_chars)
const _fermi_symbol_chars = Symbol.(_fermi_chars)

op_symbols(::Type{<:AbstractFermiOp}, ::Type{<:AbstractChar}) = _fermi_chars
op_symbols(::Type{<:AbstractFermiOp}, ::Type{<:AbstractString}) = _fermi_string_chars
op_symbols(::Type{<:AbstractFermiOp}, ::Type{Symbol}) = _fermi_symbol_chars

FermiOp(s::Union{AbstractChar, AbstractString, Symbol}) = _AbstractOp(FermiOp, s)

function _show_op_plain(io::IO, op::FermiOp)
    i = op_index(op) + 1
    if i < 1 || i > 8
        i = 8  # uninitialized data is coerced to 8 which means NoOp
    end
    print(io, _fermi_chars[i])
end

function Base.show(io::IO, op::FermiOp)
    print(io, typeof(op), ": ")
    _show_op_plain(io, op)
end

"""
    ferm_op_mult

Multiplication for single-fermion operators. We include `Z=N-E`. Multiplication does not
keep track of phase incurred when `Z` is an operand
"""
const ferm_op_mult =
    #   I       N          E        +         -         0       Z
    ((I, NumberOp, Empty, Raise, Lower, Zero, Z),      # I
     (NumberOp, NumberOp, Zero, Raise, Zero, Zero, Z),    # N
     (Empty, Zero, Empty, Zero, Lower, Zero, Empty),  # E
     (Raise, Zero, Raise, Zero, NumberOp, Zero, Raise), # +
     (Lower, Lower, Zero, Empty, Zero, Zero, Lower),  # -
     (Zero, Zero, Zero, Zero, Zero, Zero, Zero),      # 0
     (Z, NumberOp, Empty, Raise, Lower, Zero, I))      # Z

# const ferm_op_mult =
#     #   I       N          E        +         -         0
#     ((I, NumberOp, Empty, Raise, Lower, Zero),   # I
#      (NumberOp, NumberOp, Zero, Raise, Zero, Zero), # N
#      (Empty, Zero, Empty, Zero, Lower, Zero),   # E
#      (Raise, Zero, Raise, Zero, NumberOp, Zero),  # +
#      (Lower, Lower, Zero, Empty, Zero, Zero),   # -
#      (Zero, Zero, Zero, Zero, Zero, Zero))      # 0

Base.:*(op1::FermiOp, op2::FermiOp) = ferm_op_mult[op_index(op1)+1][op_index(op2)+1]

function Base.inv(op::FermiOp)
    if op === I
        return I
    else
        throw(DomainError(op, "Operator has no inverse"))
    end
end

function Base.:^(op::FermiOp, n::Integer)
    n < 0 && throw(DomainError(n))
    n == 0 && return I
    if op === I || op === NumberOp || op === Empty
        return op
    end
    return Zero
end

function Base.adjoint(op::FermiOp)
    if op === I || op === Empty || op === NumberOp || op === Zero
        return op
    end
    if op === Raise
        return Lower
    end
    if op === Lower
        return Raise
    end
    throw(DomainError(op))
end

"""
    is_raise_lower(x)

Return `true` if `x` is a raising or lowering operator.
"""
is_raise_lower(x) = x === FermiOps.Raise || x === FermiOps.Lower

IsApprox.ishermitian(op::FermiOp) = !is_raise_lower(op)

"""
    count_raise_lower(op::FermiOp)

Return the sum of the numbers of raising and lowering operators in `op`.
Number and complement of number operator each count as two.
"""
function count_raise_lower(op::FermiOp)
    if is_raise_lower(op)
        return 1
    elseif op === NumberOp || op === Empty
        return 2
    else
        return 0
    end
end

const _id_mat = [1. 0.; 0. 1.]
const _zero_mat = [0. 0.; 0. 0.]
const _number_mat = [0. 0.; 0. 1.]
const _empty_mat = [1. 0.; 0. 0.]
const _raise_mat = [0. 1.; 0. 0.]
const _lower_mat = [0. 0.; 1. 0.]

# const _fermi_mats = (_id_mat, _number_mat, _empty_mat, _raise_mat, _lower_mat, _zero_mat)

function Base.Matrix(op::FermiOp)
    _fermi_mats = (_id_mat, _number_mat, _empty_mat, _raise_mat, _lower_mat, _zero_mat)
    _fermi_mats[op_index(op) + 1]
end

end # module FermiOps


text:
"""
Various functions for different calculations
"""

"""
surface_quader:
Calculates the surface of a quader.
"""
function surface_quader(W::Real, L::Real, T::Real)

    return (2*L*W+2*T*L+2*T*W)
end

"""
volume_quader:
Calculates the surface of a quader.
"""
function volume_quader(W::Real, L::Real, T::Real)

    return (W*L*T)
end

"""
ct_and_cA_quader:

Returns the characteristic thickness of a quader.
The characteristic thickness is the greates dimension normal to the greates
surface plane (characteristic cross sectional area see bellow) of the quader
"""
function char_thickness_quader(dims::Vector{T}) where T <: Real

    if length(dims) != 3
        error("char_thickness_quader expects an array length 3")
    end

    A = [dims[1]*dims[3], dims[1]*dims[2], dims[2]*dims[3]]
    n = [dims[2], dims[3], dims[1]]

    return maximum(n[A .== maximum(A)])
end

"""
char_cross_section_area_quader:

Returns the characteristic cross sectional area of a quader.
"""
function char_cross_section_area_quader(dims::Vector{T}) where T <: Real

    if length(dims) != 3
        error("char_cross_section_area_quader expects an array length 3")
    end

    A = [dims[1]*dims[3], dims[1]*dims[2], dims[2]*dims[3]]
    n = [dims[2], dims[3], dims[1]]

    return maximum(A)
end

"""
 surface_area_without_contact():

Calculates the maximal free surface area of pieces with surface
"surface_of_piece" for a given bulk density "bulk_density"
"""
function surface_area_without_contact(bulk_density::Real, weight_of_piece::Real, surface_of_piece::Real )

    no_pieces = (bulk_density/scrap_weight_of_piece)
    return (no_pieces * surface_of_piece)
end

"""
hf_radiation_inf(eps::Real, A::Real,T1::Real, T2::Real)

returns the heat flux from T1 to T2 and this formula is valid if A2 is much
bigger than A1 if epsilon2 >> 0.
"""
function hf_radiation_inf(eps::Real, A1::Real, T1::Real, T2::Real)

    return (eps * A1 * sigma * (T1^4 - T2^4))
end

"""
hf_radiation_par_surf(eps1::Real, eps2::Real, T1::Real, T2::Real)

returns the heat flux (W/m² K) between to infinite parallel surfaces from 1 to 2
"""
function hf_radiation_par_surf(eps1::Real, eps2::Real, T1::Real, T2::Real)

    return (sigma * (T1^4 - T2^4) / (1/eps1 + 1/eps2 - 1))
end

"""
hf_radiation_par_surf_b_weighted(eps::Real, b::Real, T1::Real, T2::Real)

returns the heat flux (W/m² K) between to equally sized parallel surfaces from 1
to 2 which is weighet by a factor of b
"""
function hf_radiation_par_surf_b_weighted(eps1::Real, eps2::Real, b::Real, T1::Real, T2::Real)

    return ( b * hf_radiation_par_surf(eps1, eps2, T1, T2) )
end

function H_of_T(cp_true::MaterialPropertyOfX, T::Real)::Real

    if T > Tmax - eps(typeof(T)) || T < 0 - eps(typeof(T))
        error("value " * string(T), " out of bounds for H_of_T!")
    end

    return integrate_pX(cp_true, 0.0, T)
end

"""
T_of_H:
Caluculates the temperature T out of the enthalpy curve H(T), using bin search
so thats it should be relativly fast ...

If we are using T(H) curve to estimate T is it so that H has to have no
minima or maxima? - it should, as cp is deterimed by heating curves ...
"""
function T_of_H(H_T::Vector{<:Real}, T_H::Vector{<:Real}, H::Real)

    i_low = 1
    i_high =  length(H_T)
    while true
        i = ceil(Int64, (i_low + i_high - 1)/2)

        if H < H_T[i]
            i_high = i
        else
            i_low = i
        end

        if (i_high - i_low) == 1
            break;
        end
    end

    return  ((T_H[i_high] - T_H[i_low])/(H_T[i_high] - H_T[i_low])
                    *(H - H_T[i_low]) + T_H[i_low])
end

function T_of_H(H_T::Vector{<:Real}, T_H::Vector{<:Real}, H::Vector{<:Real})::Vector{<:Real}

    T = zeros(length(H))

    for i = 1:length(H)
        T[i] =  T_of_H(H_T, T_H, H[i])
    end

    return T
end

"""
check_monolithic_growing(p::Vector{<:Real})::Nothing
Checks if values in Vector p are monolithicly growing, used for checking
enthalpy curves, see above.
"""
function check_monolithic_growing(p::Vector{<:Real})::Nothing

    for i = 2:1:length(p)
        ((p[i] < p[i-1]) ? warning("Non growing values detected for H") :
                continue)
    end
end

"""
harm_mean(x::Real...)::Real
Simple (naive) implementation of harmonic mean average
"""
function harm_mean(x::Real...)::Real
    n = length(x)
    denom = 0
    for xi in x
        denom += 1/xi
    end

    return (n/denom)
end

"""
harm_mean(x1::T,x2::T)
optimization if harm mean only gets to values - does not need any
memory allocation ...
"""
function harm_mean(x1::T,x2::T)::T where T <: Real
    denom = 1/x1 + 1/x2

    return (2/denom)
end

"""
print_struct(T::Type, strct)
Function to print a struct and its properties
"""
function print_struct(T::Type, strct)
    for fn in fieldnames(T)
        fval = getfield(strct, fn)

        println("VAR:\n", string(fn), "\nTYPE:\n{", typeof(fval),"}\nVAL:\n", fval,"\n\n")
    end
end


text:
"Abstract Lyapunov function type"
abstract type LyapunovFunction <: CertificateFunction end

"""
    ControlLyapunovFunction

Control Lyapunov function V for a control affine system
# Fields
- `V`: function V(x) that represents the CLF
- `∇V`: function ∇V(x) defining the gradient of the CLF
- `α`: class K function α(V(x))
"""
struct ControlLyapunovFunction <: LyapunovFunction
    V
    ∇V
    α
end

"""
    ControlLyapunovFunction(V, α)

Construct a CLF given V(x) and class K function α(V(x)).
"""
function ControlLyapunovFunction(V, α)
    ∇V(x) = gradient(V, x)[1]'
    return ControlLyapunovFunction(V, ∇V, α)
end

"""
	(V::ControlLyapunovFunction)(x)

Evaluate CLF at state x.
"""
function (V::ControlLyapunovFunction)(x)
    return V.V(x)
end


text:
"""
    getindex(dict, k)

Returns the value `dict[k]`, or throws a `KeyError` if key `k` is not present.
"""
function Base.getindex(dict::MuttsDict{K,V}, k::K) where {K,V}
    hashcode = hash(k)
    node = dict.root

    # INodes
    height = length(dict.config.inode_capacity)
    while isa(node, INode{K,V})
        entries = node.entries
        mask = reinterpret(UInt64, length(entries)-1)
        shift = reinterpret(UInt64, height+1)*UInt64(8)
        slot = 1 + reinterpret(Int64, ((hashcode >> shift) & mask))
        node = @inbounds entries[slot]
        height = height-1
    end

    # Leaf node
    # Note: from the exit condition of the above while loop, the compiler
    # knows ¬isa(node,INode{K,V}) and infers node::Leaf{K,V}.
    num_slots = length(node.entries)
    max_probe2 = min(max_probe-1, num_slots-1)
    for i=0:max_probe2
        index = _hashindex(hashcode, i, num_slots)
        if _slot_occupied(node, index, num_slots)
            (ki,vi) = @inbounds node.entries[index]
            if (k === ki) || isequal(k, ki)
                return vi
            end
        else
            break
        end
    end
    throw(KeyError(k))
end

# d[k1, k2, ks...] is syntactic sugar for d[(k1, k2, ks...)].
getindex(d::MuttsDict, k1, k2, ks...) = getindex(d, tuple(k1, k2, ks...))

function Base.in(_, dict::MuttsDict{K,V}) where {K,V}
    error("""MuttsDict collections only contain Pairs; either look for e.g. A=>B
             instead, or use haskey(k,dict).""")
end

"""
    in((k=>v),dict)

Returns true just when dict contains the key-value pair (k,v).
"""
function Base.in(kv::Pair{K,V}, dict::MuttsDict{K,V}) where {K,V}
    (k,v) = kv
    return Base.in((k,v),dict)
end

"""
    in(k,dict::MuttsDict{K,Nothing})

Useful when using MuttsDict to represent a Set{K}.
"""
function Base.in(k::K, dict::MuttsDict{K,Nothing}) where {K}
    return in((k,nothing), dict)
end

"""
    in((k,v),dict)

Returns true just when dict contains the key-value pair (k,v).
"""
function Base.in((k,v)::Tuple{K,V}, dict::MuttsDict{K,V}) where {K,V}
    hashcode = hash(k)
    node = dict.root

    # Note: there is obvious massive overlap between the implementation of
    # getindex() and in(). This should be ultra-low churn code so I am fine
    # with that.

    # INodes
    height = length(dict.config.inode_capacity)
    while isa(node, INode{K,V})
        entries = node.entries
        mask = reinterpret(UInt64, length(entries)-1)
        shift = reinterpret(UInt64, height+1)*UInt64(8)
        slot = 1 + reinterpret(Int64, ((hashcode >> shift) & mask))
        node = @inbounds entries[slot]
        height = height-1
    end

    # Leaf node
    num_slots = length(node.entries)
    max_probe2 = min(max_probe-1, num_slots-1)
    for i=0:max_probe2
        index = _hashindex(hashcode, i, num_slots)
        if _slot_occupied(node, index, num_slots)
            (ki,vi) = @inbounds node.entries[index]
            if (k === ki) || isequal(k,ki)
                return (v === vi) || isequal(v, vi)
            end
        else
            break
        end
    end
    return false
end

"""
    get(dict, k, default)

Returns the value `dict[k]`, or `default` if key `k` is not present.
"""
function Base.get(dict::MuttsDict{K,V}, k::K, default) where {K,V}
    hashcode = hash(k)
    node = dict.root

    # INodes
    height = length(dict.config.inode_capacity)
    while isa(node, INode{K,V})
        entries = node.entries
        mask = reinterpret(UInt64, length(entries)-1)
        shift = reinterpret(UInt64, height+1)*UInt64(8)
        slot = 1 + reinterpret(Int64, ((hashcode >> shift) & mask))
        node = @inbounds entries[slot]
        height = height-1
    end

    # Leaf node
    # Note: from the exit condition of the above while loop, the compiler
    # knows ¬isa(node,INode{K,V}) and infers node::Leaf{K,V}.
    num_slots = length(node.entries)
    max_probe2 = min(max_probe-1, num_slots-1)
    for i=0:max_probe2
        index = _hashindex(hashcode, i, num_slots)
        if _slot_occupied(node, index, num_slots)
            (ki,vi) = @inbounds node.entries[index]
            if (k === ki) || isequal(k,ki)
                return vi
            end
        else
            break
        end
    end
    return default
end

"""
    haskey(dict,k)

Returns true just when dict contains the key k.
"""
function Base.haskey(dict::MuttsDict{K,V}, k::K) where {K,V}
    return get(dict, k, Base.secret_table_token) != Base.secret_table_token
end


text:
# Population-based evolution to try to approximately replicate Labar & Adami 2016.
# The Julia type PredType is defined in Chromosome.jl
using Main.CGP
using Statistics
using DataFramesMeta
using HypothesisTests

# Run both one population and numpops subpopulations to eithe return a dataframe or
#   write the dataframe to csvfile if it is given as a keyword argument
# mutrate is the probability that a chromosome (circuit) will be mutated in a generation
function run_both_pop_evolve( nreps::Int64, p::Parameters, popsize::Int64, goallistlist::Vector{GoalList},
    ngens::Int64, numpops::Int64, mutrate::CGP.FloatRange=1.0;
    csvfile::String="", uniform_start::Bool=false, prdebug::Bool=false )
  df=DataFrame()
  df.goallist = GoalList[]
  df.numinputs = Int64[]
  df.numoutputs = Int64[]
  df.numgates = Int64[]
  df.levelsback = Int64[]
  df.ngens = Int64[]
  df.numpops = Int64[]
  df.popsize = Int64[]
  df.mutrate = Float64[]
  df.evolvability = Int64[]
  df.maxfit = Float64[]
  df.maxfit_gen = Int64[]
  df.fit_decreases = Int64[]

  mr_goallists_list = [ (mutr,gl) for mutr = mutrate for gl in goallistlist ]
  #df_list = pmap(mr_goallist_pair->run_both_one_goallist(nreps,p,popsize,ngens,numpops,mr_goallist_pair,uniform_start=uniform_start),mr_goallists_list)
  df_list = map(mr_goallist_pair->run_both_one_goallist(nreps,p,popsize,ngens,numpops,mr_goallist_pair,uniform_start=uniform_start),mr_goallists_list)
  df = vcat( df_list... )

  #=
  for goallist in goallistlist
    println("goallist: ",goallist)
    df = run_pop_evolve(nreps,p,popsize,goallist,ngens,df=df,uniform_start=uniform_start)
    mpopsize = Int64(round(popsize/numpops))   # popsize should be divisible by numpops
    df = run_multiple_pop_evolve(nreps,p,mpopsize,goallist,ngens,numpops,mutrate,df=df,uniform_start=uniform_start)
  end
  =#
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# ngens: ",ngens)
      println(f,"# nreps: ",nreps)
      println(f,"# uniform start: ",uniform_start)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end

# Helper function for run_both_pop_evolve().  Used in pmap()
function run_both_one_goallist( nreps::Int64, p::Parameters, popsize::Int64, ngens::Int64, numpops::Int64,
    mr_goal_pair; uniform_start::Bool=false, prdebug::Bool=false )
  (mutrate,goal) = mr_goal_pair
  df = DataFrame()
  df.goallist = GoalList[]
  df.numinputs = Int64[]
  df.numoutputs = Int64[]
  df.numgates = Int64[]
  df.levelsback = Int64[]
  df.ngens = Int64[]
  df.numpops = Int64[]
  df.popsize = Int64[]
  df.mutrate = Float64[]
  df.evolvability = Int64[]
  df.maxfit = Float64[]
  df.maxfit_gen = Int64[]
  df.fit_decreases = Int64[]
  println("goal: ",goal)
  df = run_pop_evolve(nreps,p,popsize,goal,ngens,mutrate,df=df,uniform_start=uniform_start)
  #println("run_both_one_goal A size(df): ",size(df))
  mpopsize = Int64(round(popsize/numpops))   # popsize should be divisible by numpops
  df = run_multiple_pop_evolve(nreps,p,mpopsize,goal,ngens,numpops,mutrate,df=df,uniform_start=uniform_start)
  #println("run_both_one_goal B size(df): ",size(df))
  df
end

function run_multiple_pop_evolve( nreps::Int64, p::Parameters, popsize::Int64, gl::GoalList, ngens::Int64, numpops::Int64,
    mutrate::Float64=1.0;
    df::DataFrame=DataFrame(), csvfile::String="", uniform_start::Bool=false, prdebug::Bool=false )
  println("run_multiple_pop_evolve  mutrate: ",mutrate,"   goallist: ",gl)
  if size(df)[1] == 0
    df.goallist = GoalList[]
    df.numinputs = Int64[]
    df.numoutputs = Int64[]
    df.numgates = Int64[]
    df.levelsback = Int64[]
    df.ngens = Int64[]
    df.numpops = Int64[]
    df.popsize = Int64[]
    df.mutrate = Float64[]
    df.evolvability = Int64[]
    df.maxfit = Float64[]
    df.maxfit_gen = Int64[]
    df.fit_decreases = Int64[]
  end
  for rep = 1:nreps
    (mutrate,mphenotypes,mpredecessors,mfitnesses,mcomplexities,mnumgates)=multiple_pop_evolve(p,popsize,gl,ngens,numpops,mutrate,uniform_start=uniform_start)
    #println("run_multiple_pop_evolve size(mphenotypes): ",size(mphenotypes),"  size(mphenotypes[1]): ",size(mphenotypes[1]))
    push!(df.goallist,gl)
    push!(df.numinputs,p.numinputs)
    push!(df.numoutputs,p.numoutputs)
    push!(df.numgates,p.numinteriors)
    push!(df.levelsback,p.numlevelsback)
    push!(df.ngens,ngens)
    push!(df.numpops,numpops)
    push!(df.popsize,popsize)
    push!(df.mutrate,mutrate)
    push!(df.evolvability,final_evolvability(mphenotypes))
    findmax_mfh=findmax(max_fitness_history(mfitnesses))
    push!(df.maxfit,findmax_mfh[1][1])
    push!(df.maxfit_gen,findmax_mfh[2])
    push!(df.fit_decreases,preds_fit_decreases( mpredecessors, mphenotypes, mfitnesses )[1] )
  end
  df
end

function run_pop_evolve( nreps::Int64, p::Parameters, popsize::Int64, gl::GoalList, ngens::Int64, mutrate::Float64=1.0;
    df::DataFrame=DataFrame(), csvfile::String="", uniform_start::Bool=false, prdebug::Bool=false )
  println("run_pop_evolve  mutrate: ",mutrate,"   goallist: ",gl)
  if size(df)[1] == 0
    df.goallist = GoalList[]
    df.numinputs = Int64[]
    df.numoutputs = Int64[]
    df.numgates = Int64[]
    df.levelsback = Int64[]
    df.ngens = Int64[]
    df.numpops = Int64[]
    df.popsize = Int64[]
    df.mutrate = Float64[]
    df.evolvability = Int64[]
    df.maxfit = Float64[]
    df.maxfit_gen = Int64[]
    df.fit_decreases = Int64[]
  end
  for rep = 1:nreps
    (pop,phenotypes,predecessors,fitnesses,complexities,numgates)=
        pop_evolve(p,popsize,gl,ngens,mutrate,uniform_start=uniform_start)
    #println("run_pop_evolve size(phenotypes): ",size(phenotypes))
    push!(df.goallist,gl)
    push!(df.numinputs,p.numinputs)
    push!(df.numoutputs,p.numoutputs)
    push!(df.numgates,p.numinteriors)
    push!(df.levelsback,p.numlevelsback)
    push!(df.ngens,ngens)
    push!(df.numpops,1)
    push!(df.popsize,popsize)
    push!(df.mutrate,mutrate)
    push!(df.evolvability,final_evolvability(phenotypes))
    findmax_mfh=findmax(max_fitness_history(fitnesses))
    push!(df.maxfit,findmax_mfh[1])
    push!(df.maxfit_gen,findmax_mfh[2])
    push!(df.fit_decreases,preds_fit_decreases( predecessors, phenotypes, fitnesses )[1] )
    if nreps==1
      return (df,pop,phenotypes,predecessors,fitnesses,complexities,numgates)
    end
  end
  df
end

function multiple_pop_evolve( p::Parameters, popsize::Int64, gl::GoalList, ngens::Int64, numpops::Int64, mutrate::Float64=1.0;
    uniform_start::Bool=false, prdebug::Bool=false )
  mphenotypes = fill(zeros(MyInt,0,0),numpops)
  mpredecessors = fill(zeros(PredType,0,0),numpops)
  mfitnesses = fill(zeros(Float64,0,0),numpops)
  mcomplexities = fill(zeros(Float64,0,0),numpops)
  mnumgates = fill(zeros(Int64,0,0),numpops)
  for pind = 1:numpops
    #println("pind: ",pind)
    (pop,mphenotypes[pind],mpredecessors[pind],mfitnesses[pind],mcomplexities[pind],mnumgates[pind]) =
        pop_evolve(p,popsize,gl,ngens,mutrate,uniform_start=false)
  end
  (mutrate,mphenotypes,mpredecessors,mfitnesses,mcomplexities,mnumgates)
end

function pop_evolve( p::Parameters, popsize::Int64, gl::GoalList, ngens::Int64, mutrate::Float64=1.0;
    uniform_start::Bool=false, prdebug::Bool=false )
  #println("pop_evolve: goals: ",gl,"  psize: ",popsize,"  ngens: ",ngens,"  mutrate: ",mutrate)
  count_mutations = 0
  funcs = default_funcs( p.numinputs )
  phenotypes = zeros(MyInt,ngens,popsize)
  predecessors = zeros(PredType,ngens,popsize)
  fitnesses = zeros(Float64,ngens,popsize)
  complexities = zeros(Float64,ngens,popsize)
  numgates = zeros(Int64,ngens,popsize)
  if uniform_start
    c = random_chromosome(p,funcs,ident=PredType(1))
    pop = fill( c, popsize )
  else
    pop = [ random_chromosome(p,funcs,ident=PredType(i)) for i = 1: popsize ]
  end
  for i in 1:popsize
    c = pop[i]
    c.robustness = i
    c.fitness = fitness_funct( p, c, gl )
    prdebug ? print(c.robustness,"  ") : nothing
    prdebug ? print_circuit(c,include_fitness=true,include_robustness=true,include_pheno=true) : nothing
  end
  for gen = 1:ngens
    fitness_vector = rescale_fitnesses([ pop[i].fitness for i = 1:popsize ])
    #fitness_vector = ([ pop[i].fitness for i = 1:popsize ])
    prdebug ? println("fit_vect: ",fitness_vector) : nothing
    prdebug ? println("gen: ",gen,"  max fit: ",maximum(fitness_vector)) : nothing
    propsel!( pop, fitness_vector, maxfit=findmax(fitness_vector)[1] )
    prdebug ? println("after propsel gen: ",gen) : nothing
    for c in pop
      #c.fitness = hamming_distance( gl, output_values(c), p.numinputs )
      c.fitness = fitness_funct( p, c, gl )
      prdebug ? print(c.robustness,"  ") : nothing
      prdebug ? print_circuit(c,include_fitness=true,include_robustness=true,include_pheno=true) : nothing
    end
    #=
    for i in 1:length(pop)
      phenotypes[gen,i] = output_values(pop[i])[1]
      predecessors[gen,i] = pop[i].robustness
      fitnesses[gen,i] = pop[i].fitness
      complexities[gen,i] = complexity5(pop[i])
      numgates[gen,i] = pop[i].params.numinteriors
    end
    if prdebug
      println("phenos:  ",phenotypes[gen,:])
      println("preds:   ",predecessors[gen,:])
      println("fitness: ",fitnesses[gen,:])
      println("complex: ",complexities[gen,:])
      #println("numgates:",numgates[gen,:])
    end
    =#
    sav_pop = deepcopy( pop )
    for i in 1:popsize
      c = pop[i]
      if rand() <= mutrate
        #print("befor mut: ")
        #print_circuit(c)
        pop[i] = c = mutate_chromosome!(deepcopy(c),funcs)[1]
        #print("after mut: ")
        #print_circuit(c)
        #print("newc:      ")
        #print_circuit(new_c)
        #c.fitness = 1.0-hamming_distance( g, output_values(c), p.numinputs )
        c.fitness = fitness_funct( p, c, gl )
        #println("mutation count: ",count_mutations,"  phenotype: ",@sprintf("0x%x",(output_values(c)[1])),"  fitness: ",c.fitness)
        count_mutations += 1
      end
      c.robustness = i  # Robustness is not robustness; instead it is an identifier of the index of circut in the population
    end

    for i in 1:length(pop)
      phenotypes[gen,i] = output_values(pop[i])[1]
      predecessors[gen,i] = pop[i].robustness
      fitnesses[gen,i] = pop[i].fitness = fitness_funct( p, c, gl )
      @assert pop[i].fitness == fitness_funct( p, c, gl )
      complexities[gen,i] = complexity5(pop[i])
      numgates[gen,i] = pop[i].params.numinteriors
    end
    if prdebug
      println("phenos:  ",phenotypes[gen,:])
      println("preds:   ",predecessors[gen,:])
      println("fitness: ",fitnesses[gen,:])
      println("complex: ",complexities[gen,:])
      #println("numgates:",numgates[gen,:])
    end

    prdebug ? println("after mutate gen: ",gen) : nothing
    #=
    for i in 1:popsize
      c = pop[i]
      c.robustness = i  # Robustness is not robustness; instead it is an identifier of the index of circut in the population
      #c.fitness = 1.0-hamming_distance( g, output_values(c), p.numinputs )
      c.fitness = fitness_funct( p, c, gl )
      prdebug ? print(c.robustness,"  ") : nothing
      prdebug ? print_circuit(c,include_fitness=true,include_robustness=true,include_pheno=true) : nothing
    end
    =#
  end
  #println("end pop evolve: len(pop): ",length(pop),"  size: ",size(phenotypes))
  (pop,phenotypes,predecessors,fitnesses,complexities,numgates)
end

function fitness_funct( p::Parameters, c::Chromosome, gl::GoalList )
  maximum( 1.0-hamming_distance( g, output_values(c), p.numinputs ) for g in gl )
  #complexity5(c)
end

function rescale_fitnesses( fit_vect::Vector{Float64} )
  fit_min = minimum( fit_vect )
  #fit_min = quantile( fit_vect, 0.20 )
  fit_max = maximum( fit_vect )
  frange = fit_max - fit_min
  if frange > 0.0
    return [ (f-fit_min)/frange for f in fit_vect ]
  else
    return fit_vect
  end
end

function randgoallistlist( ngoals::Int64, ngoallists::Int64, p::Parameters )
  [ randgoallist( ngoals, p ) for _ = 1:ngoallists ]
end

function pheno_history( phenotypes::Array{MyInt,2} )
  (num_gens, popsize) = size(phenotypes)
  pheno_history = fill(Set(MyInt[]),num_gens)
  pheno_set = Set(MyInt[])
  for g = 1:num_gens
    union!(pheno_set,Set(phenotypes[g,:]))
    pheno_history[g] = deepcopy(pheno_set)
  end
  pheno_history
end

# Result should agree with evolvability()
function pheno_history( mphenotypes::Array{Array{MyInt,2},1} )
  npops = length(mphenotypes)
  (num_gens, popsize) = size(mphenotypes[1])
  #phenohistories = [pheno_history( mphenotypes[i] ) for i = npops ]  # fails due to a bug in Julia
  phenohistories =  Array{Set{UInt16},1}[]
  for i = 1:npops
    push!(phenohistories,pheno_history( mphenotypes[i]))
  end
  phenoHistory = fill(Set(MyInt[]),num_gens)
  for g = 1:num_gens
    phenoHistory[g] = phenohistories[1][g]
    for i = 1:npops
      phenoHistory[g] = union!(phenoHistory[g],phenohistories[i][g])
    end
  end
  pheno_set = Set(MyInt[])
  for g = 1:num_gens
    union!(pheno_set,Set(phenoHistory[g]))
    phenoHistory[g] = deepcopy(pheno_set)
  end
  phenoHistory
end

function final_pheno_history( mphenotypes::Array{Array{MyInt,2},1} )
  length(pheno_history(mphenotypes)[end])
end

function evolvability( phenotypes::Array{MyInt,2} )
  (num_gens, popsize) = size(phenotypes)
  evo_history = zeros(Int64,num_gens)
  pheno_set = Set(MyInt[])
  for i = 1:num_gens
    union!(pheno_set,Set(phenotypes[i,:]))
    evo_history[i] = length(pheno_set)
  end
  pheno_set
  #evo_history
end

function final_evolvability( phenotypes::Array{MyInt,2} )
  length(evolvability(phenotypes))
end

# Result should agree with pheno_history()
function evolvability( mphenotypes::Array{Array{MyInt,2},1} )
  npops = length(mphenotypes)
  (num_gens, popsize) = size(mphenotypes[1])
  # Take the union of the phenotypes of each population for each generation g
  ph_gen = [union( map( union, [ mphenotypes[j][g,:] for j = 1:npops] )... ) for g = 1:num_gens]
  pheno_history = fill(Set(MyInt[]),num_gens)
  pheno_set = Set(MyInt[])
  for g = 1:num_gens
    union!(pheno_set,Set(ph_gen[g]))
    pheno_history[g] = deepcopy(pheno_set)
  end
  pheno_history
end

function final_evolvability( mphenotypes::Array{Array{MyInt,2},1} )
  length(evolvability(mphenotypes)[end])
end

function max_fitness_history( fitnesses::Array{Float64,2} )
  (num_gens, popsize) = size(fitnesses)
  history = [maximum( fitnesses[g,:]) for g = 1:num_gens ]
end

function max_fitness_history( mfitnesses::Array{Array{Float64,2},1} )
  numpops = length(mfitnesses)
  (num_gens, popsize) = size(mfitnesses[1])
  histories = [[ maximum(mfitnesses[i][g,:]) for g = 1:num_gens ] for i = 1:numpops ]
  history = [findmax( [histories[i][g] for i = 1:numpops] ) for g = 1:num_gens ]
end

function mean_fitness_history( fitnesses::Array{Float64,2} )
  (num_gens, popsize) = size(fitnesses)
  history = [mean( fitnesses[g,:]) for g = 1:num_gens ]
end

function max_complexity_history( complexities::Array{Float64,2} )
  (num_gens, popsize) = size(complexities)
  history = [maximum( complexities[g,:]) for g = 1:num_gens ]
end

function mean_complexity_history( complexities::Array{Float64,2} )
  (num_gens, popsize) = size(complexities)
  history = [mean( complexities[g,:]) for g = 1:num_gens ]
end

# Return the list of unique predecessors of element index of populations 1 to gen.
# Usually one will call with gen = maxgenerations
function preds( predecessors::Array{Int64,2}, gen::Int64, index::Int64 )
  indices = zeros(Int64,gen)
  while gen > 0
    indices[gen] = index = predecessors[gen,index]
    #println("gen: ",gen,"  index: ",index)
    gen -= 1
  end
  indices
end

# Returns a list of phenotype-fitness pairs of the predecessors of a random maximum fitness element
#   of the first population to reach the maximum fitness.
function preds_phenos_fits( predecessors::Array{Int64,2}, phenotypes::Array{MyInt,2}, fitnesses::Array{Float64,2} )
  #println("preds_phenos_fits: size(phenotypes): ",size(phenotypes))
  (maxfit,maxfit_gen) = findmax_mfh=findmax(max_fitness_history(fitnesses))
  column = rand(findmaxall( fitnesses[maxfit_gen,:])[2])  # A random columm with maximum fitness
  prds = preds( predecessors, maxfit_gen, column )
  [(phenotypes[i,prds[i]],fitnesses[i,prds[i]]) for i = 1:maxfit_gen ]
end

function preds_phenos_fits( mpredecessors::Array{Array{Int64,2},1}, mphenotypes::Array{Array{MyInt,2},1},
    mfitnesses::Array{Array{Float64,2},1} )
  #println("preds_phenos_fits  size(mphenotypes[1]): ",size(mphenotypes[1]))
  ((maxfit,maxfit_pop),maxfit_gen) = findmax_mfh=findmax(max_fitness_history(mfitnesses))
  #println("preds_phenos_fits  size(mfitnesses[maxfit_pop]): ",size(mfitnesses[maxfit_pop]))
  column = rand(findmaxall( mfitnesses[maxfit_pop][maxfit_gen,:])[2])  # A random columm with maximum fitness
  prds = preds( mpredecessors[maxfit_pop], maxfit_gen, column )
  [(mphenotypes[maxfit_pop][i,prds[i]],mfitnesses[maxfit_pop][i,prds[i]]) for i = 1:maxfit_gen ]
end

# Returns the number of fitness decreases in the fitness history of a random maximum fitness
#   element of the first population with maximum fitness.
# Also returns the maximum fitness over all generations and the first generation to reach this fitness.
function preds_fit_decreases( predecessors::Array{Int64,2}, phenotypes::Array{MyInt,2}, fitnesses::Array{Float64,2} )
  #println("preds_fit_decreases: size(phenotypes): ",size(phenotypes))
  ppf = preds_phenos_fits( predecessors, phenotypes, fitnesses )
  maxfit_gen = length(ppf)
  fits = map( x->x[2], ppf )
  decreases = 0
  for i = 1:(maxfit_gen-1)
    if fits[i+1] < fits[i]
      decreases += 1
    end
  end
  (decreases, fits[maxfit_gen], maxfit_gen)
end

function preds_fit_decreases( mpredecessors::Array{Array{Int64,2},1}, mphenotypes::Array{Array{MyInt,2},1},
    mfitnesses::Array{Array{Float64,2},1} )
  #println("preds_fit_decreases: size(mphenotypes[1]): ",size(mphenotypes[1]))
  ppf = preds_phenos_fits( mpredecessors, mphenotypes, mfitnesses )
  maxfit_gen = length(ppf)
  fits = map( x->x[2], ppf )
  decreases = 0
  for i = 1:(maxfit_gen-1)
    if fits[i+1] < fits[i]
      decreases += 1
    end
  end
  (decreases, fits[maxfit_gen], maxfit_gen)
end

# Return the list of successors of element index of population gen to num_gens
# Usually one will call with gen = 1
function successors( predecessors::Array{Int64,2}, gen::Int64, index::Int64 )
  (num_gens, popsize) = size(predecessors)
  succs = Vector{Int64}[]
  succ = [index]
  for i = gen:num_gens
    succ = findall(x->x in succ,predecessors[i,:])
    push!(succs,succ)
  end
  succs
end

# Check that preds(successors(i) == i (approximately stated)
# Check that for each s in successors(i), preds(s) == i
function test_preds_successors( predecessors::Array{Int64,2} )
  (num_gens, num_indices) = size(predecessors)
  for i = 1:num_indices
    succs = successors( predecessors, 1, i )[num_gens]
    for s in succs
      if i != (prd=preds( predecessors, num_gens, s )[1] )
        println("i: ",i,"  succs: ",succs,"  s: ",s,"  prd: ",prd)
      end
    end
  end
end

# Check that for each c=1:popsize, s in successors(preds(c))
function test_successors_preds( predecessors::Array{Int64,2} )
  (num_gens, popsize) = size(predecessors)
  for c = 1:num_indices
    p = preds(predecessors, num_gens, c )[1]
    succs = successors( predecessors, 1, p )[num_gens]
    if !( c in succs )
      println("p:  ",p,"  c: ",c,"  succs: ",succs)
    end
  end
end

function compare_fract_decreases( df::DataFrame, csvfile::String )
  cdf=compare_fract_decreases( df )
  write_dataframe_with_comments(cdf,csvfile,"$(csvfile[1:(end-4)])fract_dec.csv")
  cdf
end

function compare_fract_decreases( csvfile::String )
  df = read_dataframe( csvfile )
  cdf=compare_fract_decreases( df )
  write_dataframe_with_comments(cdf,csvfile,"$(csvfile[1:(end-4)])fract_dec.csv")
  cdf
end

# Returns a triple of the large population mean, the small population mean, and the p-value of the Mann Whitney test
function compare_fract_decreases( df::DataFrame )
  df.fract_dec = df.fit_decreases./df.maxfit_gen
  num_small_pops = unique(df.numpops)[2]    # The number of pops in the small pop case
  cdf=DataFrame()
  if ("numinputs" in names(df)) || (:numinputs in names(df))
    cdf.numinputs = Int64[]
    cdf.numoutputs = Int64[]
    cdf.numgates = Int64[]
    cdf.levelsback = Int64[]
  end
  cdf.ngens=Int64[]
  cdf.lg_popsize = Int64[]
  cdf.sm_popsize = Int64[]
  cdf.mutrate = Float64[]
  cdf.lge_pop_fdecr = Float64[]
  cdf.sml_pop_fdecr = Float64[]
  cdf.MannWhitPval = Float64[]
  for mr in unique(df.mutrate)
    ldf = @where(df,(:maxfit.>=1.0).&(:numpops.==1).&(:mutrate.==mr))
    large_pop_fract_dec = mean(ldf.fract_dec)
    sdf = @where(df,(:maxfit.>=1.0).&(:numpops.==num_small_pops).&(:mutrate.==mr))
    small_pops_fract_dec = mean(sdf.fract_dec)
    println("mr: ",mr,"  length(ldf.fract_dec): ",length(ldf.fract_dec),"  length(sdf.fract_dec): ",length(sdf.fract_dec))
    if length( sdf.fract_dec ) == 0 || length( ldf.fract_dec ) == 0
      println("warning: length(sdf) == 0  or length(ldf) == 0 for mutrate: ",mr)
    else
      MWT = pvalue(MannWhitneyUTest( ldf.fract_dec, sdf.fract_dec ))
      if ("numinputs" in names(df)) || (:numinputs in names(df))
        push!(cdf,(sdf.numinputs[1],sdf.numoutputs[1],sdf.numgates[1],sdf.levelsback[1],
            sdf.ngens[1],ldf.popsize[1],sdf.popsize[1],mr,large_pop_fract_dec,small_pops_fract_dec,MWT))
      else
        push!(cdf,(sdf.ngens[1],ldf.popsize[1],sdf.popsize[1],mr,large_pop_fract_dec,small_pops_fract_dec,MWT))
      end
    end
  end
  cdf
end


text:
# Autogenerated wrapper script for PCRE2_jll for x86_64-w64-mingw32
export libpcre2_16, libpcre2_32, libpcre2_8

JLLWrappers.@generate_wrapper_header("PCRE2")
JLLWrappers.@declare_library_product(libpcre2_16, "libpcre2-16-0.dll")
JLLWrappers.@declare_library_product(libpcre2_32, "libpcre2-32-0.dll")
JLLWrappers.@declare_library_product(libpcre2_8, "libpcre2-8-0.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpcre2_16,
        "bin\\libpcre2-16-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcre2_32,
        "bin\\libpcre2-32-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcre2_8,
        "bin\\libpcre2-8-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using OneDmin

using NLPModels
#using CUTEst
using OptimizationProblems
using NLPModelsJuMP
using SolverTools   # Pour avoir les utilitaires d'affichage log_header et log_row
using SolverCore

using Stopping
using Logging

#include("params.jl")
γ = 1.0e-5
a = 0.0
b = Inf
scale = 1.0e-5

#nlp = CUTEstModel("BEALE")
nlp = MathOptNLPModel(beale())
#nlp = MathProgNLPModel(beale())

dir = -grad(nlp, nlp.meta.x0) * scale
h = LineModel(nlp, nlp.meta.x0, dir);
dϕt0 = grad(h,a)

ϵ = max(abs(dϕt0),1.0)
α = -γ*ϵ
β =  γ*ϵ

τ₁ = 0.99
τ₀ = 0.001
g₀ = dϕt0

#logger = Logging.ConsoleLogger(stderr,Logging.Warn)
#logger = Logging.ConsoleLogger(stderr,Logging.Info)
logger = Logging.ConsoleLogger(stderr,Logging.Debug)
#logger = Logging.NullLogger()

stp = NLPStopping(h, OneDAtX(0.0,zeros(2)), tol_check=(a,b,c) -> ones(2))
stp.meta.max_iter = 100
stp.meta.atol = γ
stp.meta.rtol = γ

stp = bracket(h, a=0.0, b=Inf, α=α, β =β, pick_in = pick_ins, stp=stp)
@show stp.current_state.x, stp.current_state.fx, stp.current_state.gx, status(stp)
@show nlp.counters
@show h.counters

reset!(h)
reset!(nlp)
reinit!(stp)

stp2 = bracket_s(h, a=0.0, b=Inf, α=α, β=β, stp=stp)
@show stp2.current_state.x, stp2.current_state.fx, stp2.current_state.gx, status(stp2)
@show nlp.counters
@show h.counters

@test stp2.current_state.x == stp.current_state.x

@info "Executing bracket(h)"
stp3 = bracket(h)

@info "Executing TR1D(h)"
stp4 = TR1D(h)


text:
#######################################
#
# Read SSR
#
#######################################

"""
## Read SSR from file or IO stream
Read a file or IO stream and store the data in an `SSR` type.

Matching .mat files are read and modulation frequency information extracted.
Failing that, user passed arguments are used or the modulation frequency is extracted from the file name.

#### Arguments

* `fname`: Name of the file to be read
* `min_epoch_length`: Minimum epoch length in samples. Shorter epochs will be removed (0)
* `max_epoch_length`: Maximum epoch length in samples. Longer epochs will be removed (0 = all)
* `valid_triggers`: Triggers that are considered valid, others are removed ([1,2])
* `stimulation_amplitude`: Amplitude of stimulation (NaN)
* `modulationrate`: Modulation frequency of SSR stimulation (NaN)
* `carrier_frequency`: Carrier frequency (NaN)
* `participant_name`: Name of participant ("")
* `remove_first`: Number of epochs to be removed from start of recording (0)
* `max_epochs`: Maximum number of epochs to retain (0 = all)
* `env` (nothing)
* `bkt` ("")

#### Supported file formats

* BIOSEMI (.bdf)
"""
function read_SSR(
    fname::String;
    stimulation_amplitude::Number = NaN,   # User can set these
    modulationrate::Number = NaN,          # values, but if not
    carrier_frequency::Number = NaN,       # then attempt to read
    stimulation_side::AbstractString = "", # from file name or mat
    participant_name::AbstractString = "",
    valid_triggers::Array{Int} = [1, 2],
    min_epoch_length::Int = 0,
    max_epoch_length::Int = 0,
    remove_first::Int = 0,
    max_epochs::Int = 0,
    env = nothing,
    bkt = "",
    kwargs...,
)

    @info("Importing SSR from file: $fname")

    file_path, file_name, ext = fileparts(fname)

    if env != nothing

        @debug("File type is S3")
        fname = S3.get_object(env, bkt, fname).obj
    end

    #
    # Extract meta data
    #

    # Extract frequency from the file name if not set manually
    if occursin("Hz", file_name) && isnan(modulationrate)
        a = match(r"[-_](\d+[_.]?[\d+]?)Hz|Hz(\d+[_.]?[\d+]?)[-_]", file_name).captures
        modulationrate = parse(Float64, a[[i !== nothing for i in a]][1]) * 1.0u"Hz"
        @debug("Extracted modulation frequency from file name: $modulationrate")
    end

    # Or even better if there is a mat file read it
    mat_path = string(file_path, file_name, ".mat")
    if !isfile(mat_path)
        mat_path = string(file_path, file_name, "-properties.mat")
    end

    if isfile(mat_path)
        modulationrate,
        stimulation_side,
        participant_name,
        stimulation_amplitude,
        carrier_frequency = read_rba_mat(mat_path)
        valid_mat = true
    else
        valid_mat = false
    end

    #
    # Read file data
    #

    # Import raw data
    if ext == "bdf"
        data, triggers, system_codes, samplingrate, reference_channel, header =
            import_biosemi(fname; kwargs...)
    else
        warn("File type $ext is unknown")
    end

    # Create electrodes
    elecs = Electrode[]
    for e in header["chanLabels"]
        push!(elecs, Electrode(e, Talairach(NaN * u"m", NaN * u"m", NaN * u"m"), Dict()))
    end

    # Create SSR type
    if unit(modulationrate) == unit(1.0u"Hz")
        #nothing
    else
        modulationrate = modulationrate * 1.0u"Hz"
    end
    if unit(samplingrate) == unit(1.0u"Hz")
        #nothing
    else
        samplingrate = samplingrate * 1.0u"Hz"
    end

    a = SSR(
        data,
        elecs,
        triggers,
        system_codes,
        samplingrate,
        modulationrate,
        [reference_channel],
        file_path,
        file_name,
        Dict(),
        header,
    )

    # If a valid mat file was found then store that information with the raw header
    if valid_mat
        a.header["rbadata"] = matread(mat_path)
    end

    #
    # Store meta data in processing dictionary
    #

    if stimulation_side != ""
        a.processing["Side"] = stimulation_side
    end
    if participant_name != ""
        a.processing["Name"] = participant_name
    end
    if !isnan(stimulation_amplitude)
        a.processing["Amplitude"] = stimulation_amplitude
    end
    if !isnan(carrier_frequency)
        a.processing["Carrier_Frequency"] = carrier_frequency
    end

    #
    # Clean up
    #

    # Remove status channel information
    remove_channel!(a, "Status")

    # Clean epoch index
    a.triggers = clean_triggers(
        a.triggers,
        valid_triggers,
        min_epoch_length,
        max_epoch_length,
        remove_first,
        max_epochs,
    )

    return a
end

#######################################
#
# Read event files
#
#######################################

function read_evt(a::SSR, fname::AbstractString; kwargs...)
    d = read_evt(fname, a.samplingrate; kwargs...)
    validate_triggers(d)
    a.triggers = d
    return a
end

#######################################
#
# File IO
#
#######################################

function write_SSR(
    a::SSR,
    fname::S;
    chanLabels = channelnames(a),
    subjID = a.header["subjID"],
    startDate = a.header["startDate"],
    startTime = a.header["startTime"],
    kwargs...,
) where {S<:AbstractString}

    fname = convert(String, fname)

    @info("Saving $(size(a.data)[end]) channels to $fname")

    writeBDF(
        fname,
        a.data',
        vec(trigger_channel(a)),
        vec(system_code_channel(a)),
        samplingrate(Int, a),
        chanLabels = convert(Array{S,1}, chanLabels),
        startDate = startDate,
        startTime = startTime,
        subjID = subjID,
    )

end


text:
"""
# module PiGPIOMEM

Julia interface for `/dev/gpiomem` on Raspberry Pi [1].
`/dev/gpiomem` maps the BCM2835 GPIO registers at `0x7E200000` [2, p89].

This module can:
 * Query and set GPIO pins.
 * Configure input/output mode and pullups/pulldowns.

## Interface Summary

Create a `GPIOPin` object:
 * `pin = GPIOPin(7)`

Query pin state:
 * `x = pin[]` (`Bool`)

Set pin state:
 * `pin[] = true` -- logic high
 * `pin[] = false` -- logic low
 * `set(pin)` -- logic high
 * `clear(pin)` -- logic low

Configure input/output:
 * `set_input_mode(pin)` (default)
 * `set_output_mode(pin)`
 * `is_input(pin)::Bool`
 * `is_output(pin)::Bool`

Configure pullup/pulldown:
 * `set_highz(pin)` (default)
 * `set_pulldown(pin)`
 * `set_pullup(pin)`

Set/Get multiple pins:
 * `PiGPIOMEM.set(::GPIOPin...)`
 * `PiGPIOMEM.clear(::GPIOPin...)`
 * `PiGPIOMEM.level(::GPIOPin...)`

## Raspberry Pi GPIO Pinout

```
              Pin 1 Pin2
           +3V3 [ ] [ ] +5V
  SDA / GPIO  2 [ ] [ ] +5V
  SCL / GPIO  3 [ ] [ ] GND
        GPIO  4 [ ] [ ] GPIO 14 / TX
            GND [ ] [ ] GPIO 15 / RX
        GPIO 17 [ ] [ ] GPIO 18
        GPIO 27 [ ] [ ] GND
        GPIO 22 [ ] [ ] GPIO 23
           +3V3 [ ] [ ] GPIO 24
 MOSI / GPIO 10 [ ] [ ] GND
 MISO / GPIO  9 [ ] [ ] GPIO 25
 SCLK / GPIO 11 [ ] [ ] GPIO  8
            GND [ ] [ ] GPIO  7
ID_SD / GPIO  0 [ ] [ ] GPIO  1 / ID_SC
        GPIO  5 [ ] [ ] GND
        GPIO  6 [ ] [ ] GPIO 12
        GPIO 13 [ ] [ ] GND
        GPIO 19 [ ] [ ] GPIO 16
        GPIO 26 [ ] [ ] GPIO 20
            GND [ ] [ ] GPIO 21
             Pin 39 Pin 40
```

## References

 [1] BCM2835-ARM-Peripherals.pdf
     https://www.raspberrypi.org/app/uploads/2012/02/BCM2835-ARM-Peripherals.pdf

 [2] BCM2835 GPIO memory device driver
     https://github.com/raspberrypi/linux/blob/rpi-5.4.y/
                        drivers/char/broadcom/bcm2835-gpiomem.c
"""
module PiGPIOMEM

export GPIOPin
export set, clear
export set_input_mode, set_output_mode, is_input, is_output
export set_highz, set_pullup, set_pulldown

using Mmap

# BCM2835 GPIO Registers.

"Memory map of the BCM2835 GPIO registers (`0x7E200000`)"
const gpiomem = Ref{Vector{UInt32}}()
const gpiomem_length = 40                                        # [1, 6.1, p91]

"""
    Register(offset)

Wraps a pointer to one of the BCM2835 GPIO registers.
Provides efficient access to register content using direct load/store
instructions via the `getindex/setindex!` interface.
The offset is bounds checked once at construction time.
The `/dev/gpiomem` map is loaded on demand.
"""
struct Register
    address::Ptr{UInt32}
    function Register(offset)
        @assert 1 + offset in 1:gpiomem_length
        if !isassigned(gpiomem)
            if ispath("/dev/gpiomem")
                gpiomem[] = Mmap.mmap("/dev/gpiomem",
                                      Vector{UInt32}, gpiomem_length;
                                      grow=false)
            else
                @warn "/dev/gpiomem not found! Using dummy gpiomem."
                gpiomem[] = zeros(UInt32, gpiomem_length)
            end
        end
        new(pointer(gpiomem[], 1 + offset))
    end
end

Base.getindex(r::Register) = unsafe_load(r.address)
Base.setindex!(r::Register, v::UInt32) = unsafe_store!(r.address, v)

gpfsel(n)   = Register(n÷10)                        # 0x 7E20 0000 [1, 6.1, p90]
gpset0()    = Register(0x1C÷4)                      # 0x 7E20 001C [1, 6.1, p90]
gpclr0()    = Register(0x28÷4)                      # 0x 7E20 0028 [1, 6.1, p90]
gplev0()    = Register(0x34÷4)                      # 0x 7E20 0034 [1, 6.1, p90]
gppud()     = Register(0x94÷4)                      # 0x 7E20 0094 [1, 6.1, p91]
gppudclk0() = Register(0x98÷4)                      # 0x 7E20 0098 [1, 6.1, p91]

"""
    GPIOPin(n; [output=false])

Raspberry Pi (BCM2835) General Purpose Input Output pin `n`.

See `help?>` [PiGPIOMEM](@ref) for Interface Summary.

### Implementation Notes

The implementation aims to minimise pin access overhead.
The pin number is bounds checked only once at construction time and
register addresses and bit masks are pre-computed.

The GPIOPin object holds direct pointers to the BCM2835 Registers that
control a specific pin. Rather than storing the pin number, the GPIOPin object
stores bit masks for the register bits that control the pin. This allows
the `set` and `clear` functions to be as simple as possible
(`p.gpset[] = p.pin_bit` and `p.gpclr[] = p.pin_bit` respectively).

`sel_bit` is the Input/output mode bit for this pin in the GPFSELn register.

`pin_bit` is the bit that controls this pin in the GPSETn, GPCLRn and GPLEVn
egisters."
"""
struct GPIOPin

    sel_bit::UInt32                                       # [1, Table 6-2,  p92]
    gpfsel::Register                                      # [1,       6.1,  p90]

    pin_bit::UInt32
    gpset::Register                                       # [1, Table 6-9,  p95]
    gpclr::Register                                       # [1, Table 6-10, p95]
    gplev::Register                                       # [1, Table 6-12, p96]

    function GPIOPin(n; output=false)
        @assert n in 0:27
        pin = new(1 << (3*(n%10)), gpfsel(n),
                  1 << n, gpset0(), gpclr0(), gplev0())
        reset_mode(pin)
        if output
            set_output_mode(pin)
        end
        pin
    end
end

sel_index(p::GPIOPin) = trailing_zeros(p.sel_bit)
pin_index(p::GPIOPin) = trailing_zeros(p.pin_bit)

# Set/Get pin state.

Base.setindex!(p::GPIOPin, v::Bool) = v ? set(p) : clear(p)
Base.setindex!(p::GPIOPin, v) = setindex!(p, !iszero(v))
Base.getindex(p::GPIOPin) = !iszero(level(p))

  set(p::GPIOPin) = (p.gpset[] = p.pin_bit; nothing)      # [1, Table 6-9,  p95]
clear(p::GPIOPin) = (p.gpclr[] = p.pin_bit; nothing)      # [1, Table 6-10, p95]
level(p::GPIOPin) =  p.gplev[] & p.pin_bit                # [1, Table 6-12, p96]

# Input/Output mode.                                         [1, Table 6-2, p92]

reset_mode(p::GPIOPin) = p.gpfsel[] &= ~(UInt32(0b111) << sel_index(p))

set_input_mode(p::GPIOPin)  = (p.gpfsel[] &= ~p.sel_bit; nothing)
set_output_mode(p::GPIOPin) = (p.gpfsel[] |=  p.sel_bit; nothing)
is_input(p::GPIOPin)  = iszero(p.gpfsel[]  &  p.sel_bit)
is_output(p::GPIOPin) = !is_input(p::GPIOPin)

# Operations on multiple pins.

"Combine the `pin_bit`s for multiple `pins` into a single bit mask."
bits(pins...) = reduce(|, p.pin_bit for p in pins)

  set(pins::GPIOPin...) = (first(pins).gpset[] = bits(pins); nothing)
clear(pins::GPIOPin...) = (first(pins).gpclr[] = bits(pins); nothing)
level(pins::GPIOPin...) =  first(pins).gplev[] & bits(pins)

# Pull-up/down.                                            [1, Table 6-28, p101]

set_highz(p::GPIOPin)    = set_pud(p, 0)
set_pulldown(p::GPIOPin) = set_pud(p, 1)
set_pullup(p::GPIOPin)   = set_pud(p, 2)

using LLVM
@noinline nop() = LLVM.Interop.@asmcall("nop")
spin(n) = for i in 1:n nop() end                  # spin(50) ~= 1us on Pi Zero W

function set_pud(p::GPIOPin, mode)
    gppud()[] = mode
    spin(10)                                         # [1, p101] Wait 150 cycles
    gppudclk0()[] = p.pin_bit
    spin(10)                                         # [1, p101] Wait 150 cycles
    gppud()[] = 0
    gppudclk0()[] = 0
    nothing
end

# Pretty printing.

Base.show(io::IO, p::GPIOPin) =
    print(io, "GPIOPin(", pin_index(p), is_output(p) ? ", output=true)" : ")")

Base.show(io::IO, r::Register) =
    print(io, "BCM2835 Register: 0x7E20_",
              uppercase(string(r.address - pointer(gpiomem[]), base=16, pad=4)))

end # module


text:
#=  Abdulkadiroğlu et al. (2015).       =#

using DeferredAcceptance
using Random

function argsort(vec::AbstractArray{<:Real, 1})::AbstractArray{<:Real, 1}
    return invperm(sortperm(vec))
end

# We assume student volume equals total school capacity.
function target_doer(n, m, samp)

    CADA_disutil = zeros(0)
    STB_disutil = zeros(0)
    strat_disutil = zeros(0)

    for _ in 1:samp
        # All students share a master preference order, but they have utility perturbations
        # within it, so they "target" a school that they have less disutility with.
        students = repeat(randperm(m), 1, n)
        students_dist = rand(n)
        students_dist /= sum(students_dist)

        # Students' target schools
        targets = rand(1:m, n)

        # Each student has slightly better utility than expected at their target school.
        students_disutility = 0. .+ students
        for (s, c) in enumerate(targets)
            students_disutility[c, s] -= rand()
        end

        # Compare outcomes if one group of students simply lies and lists
        # their target school first.
        students_strat = copy(students)
        students_strat[targets[1], 1] = 0
        students_strat[:, 1] = argsort(students_strat[:, 1])

        # Schools begin with neutral preferences
        schools_in = ones(Int, n, m)
        capacities = rand(m)
        capacities /= sum(capacities)

        # Break ties
        schools_CADA = choiceaugmented(schools_in, targets)
        schools_STB = singletiebreaking(schools_in)

        out_CADA = nonatomicdeferredacceptance(students, students_dist, schools_CADA, capacities)
        out_STB = nonatomicdeferredacceptance(students, students_dist, schools_STB, capacities)
        out_strat = nonatomicdeferredacceptance(students_strat, students_dist, schools_STB, capacities)

        append!(CADA_disutil, sum(out_CADA[1][1:m, 1] .* students_disutility[:, 1]))
        append!(STB_disutil, sum(out_STB[1][1:m, 1] .* students_disutility[:, 1]))
        append!(strat_disutil, sum(out_strat[1][1:m, 1] .* students_disutility[:, 1]))
    end
    return CADA_disutil, STB_disutil, strat_disutil
end

function target_doer_discrete(m, cap, samp)
    n = m * cap

    CADA_disutil = zeros(0)
    STB_disutil = zeros(0)
    strat_disutil = zeros(0)

    for _ in 1:samp
        # All students share a master preference order, but they have utility perturbations
        # within it, so they "target" a school that they have less disutility with.
        students = repeat(randperm(m), 1, n)

        # Students' target schools
        targets = rand(1:m, n)

        # Each student has slightly better utility than expected at their target school.
        students_disutility = 0. .+ students
        for (s, c) in enumerate(targets)
            students_disutility[c, s] -= rand()
        end

        # Compare outcomes if student simply lies and lists
        # their target school first.
        students_strat = copy(students)
        students_strat[targets[1], 1] = 0
        students_strat[:, 1] = argsort(students_strat[:, 1])

        # Schools begin with neutral preferences
        schools_in = ones(Int, n, m)
        capacities = ones(Int, m) .* cap

        # Break ties
        schools_CADA = choiceaugmented(schools_in, targets)
        schools_STB = singletiebreaking(schools_in)

        out_CADA = deferredacceptance(students, schools_CADA, capacities)
        out_STB = deferredacceptance(students, schools_STB, capacities)
        out_strat = deferredacceptance(students_strat, schools_STB, capacities)

        append!(CADA_disutil, students_disutility[out_CADA[1][1]])
        append!(STB_disutil, students_disutility[out_STB[1][1]])
        append!(strat_disutil, students_disutility[out_strat[1][1]])
    end
    return CADA_disutil, STB_disutil, strat_disutil
end

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
"""
This module is used to define several core, abstract types, in order to keep
compilation simple, and to avoid any circular dependencies between modules.
"""
module core
export  IndividualCollection, Operator, OperatorDefinition, indexed_fitnesses,
        individuals_from_ids, write_to_individual_collection

"""
Used to hold a collection of individuals.
"""
type IndividualCollection{F}
  """
  The fitness values for each individual within this collection, indexed by
  their internal ID.
  """
  fitnesses::Vector{F}

  """
  The developmental stages for each individual within this collection. Indexed
  by the name of the stage, then by the individual's internal ID.
  """
  stages::Dict{AbstractString, Any}

  """
  Constructs a new, empty collection of individuals.
  """
  IndividualCollection() = new([], Dict())

  IndividualCollection(f::Vector{F}, s::Dict{AbstractString, Any}) = new(f, s)
end

"""
Returns a vector containing the fitness of each individual within a given
collection along with their internal ID, in the form of a tuple.
"""
indexed_fitnesses{F}(ic::IndividualCollection{F}) =
  collect(enumerate(ic.fitnesses))

"""
Composes a new collection of individuals by selecting the multi-set of
individuals in a given collection specified by a provided set of internal IDs.
"""
function individuals_from_ids{F}(from::IndividualCollection{F}, ids::Vector{Int})
  to = IndividualCollection{F}()
  to.fitnesses = [from.fitnesses[id] for id in ids]
  for stage in keys(from.stages)
    to.stages[stage] = eltype(from.stages[stage])[from.stages[stage][id] for id in ids]
  end
  to
end

function write_to_individual_collection{F}(to::IndividualCollection{F}, from::IndividualCollection{F}, ids::Vector{Int})
  # Copy fitness.
  for (i, id) in enumerate(ids)
    to.fitnesses[i] = from.fitnesses[id]
  end

  # Copy chromosomes.
  for stage in keys(from.stages)
    for (i, id) in enumerate(ids)
      to.stages[stage][i] = from.stages[stage][id]
    end
  end

  # Return modified buffer.
  to
end

"""
The base type used by all search operators.
"""
abstract Operator

"""
The base type used by all search operator definitions.
"""
abstract OperatorDefinition
end