text:
const LookupArraySetters = Union{AllMetadata,LookupArray,LookupArrayTrait,Nothing}
set(lookup::LookupArray, x::LookupArraySetters) = _set(lookup, x)

# _set(lookup::LookupArray, newlookup::LookupArray) = lookup
_set(lookup::AbstractCategorical, newlookup::AutoLookup) = begin
    lookup = _set(lookup, parent(newlookup))
    o = _set(order(lookup), order(newlookup))
    md = _set(metadata(lookup), metadata(newlookup))
    rebuild(lookup; order=o, metadata=md)
end
_set(lookup::LookupArray, newlookup::AbstractCategorical) = begin
    lookup = _set(lookup, parent(newlookup))
    o = _set(order(lookup), order(newlookup))
    md = _set(metadata(lookup), metadata(newlookup))
    rebuild(newlookup; data=parent(lookup), order=o, metadata=md)
end
_set(lookup::AbstractSampled, newlookup::AutoLookup) = begin
    # Update index
    lookup = _set(lookup, parent(newlookup))
    o = _set(order(lookup), order(newlookup))
    sa = _set(sampling(lookup), sampling(newlookup))
    sp = _set(span(lookup), span(newlookup))
    md = _set(metadata(lookup), metadata(newlookup))
    rebuild(lookup; data=parent(lookup), order=o, span=sp, sampling=sa, metadata=md)
end
_set(lookup::LookupArray, newlookup::AbstractSampled) = begin
    # Update each field separately. The old lookup may not have these fields, or may have
    # a subset with the rest being traits. The new lookup may have some auto fields.
    lookup = _set(lookup, parent(newlookup))
    o = _set(order(lookup), order(newlookup))
    sp = _set(span(lookup), span(newlookup))
    sa = _set(sampling(lookup), sampling(newlookup))
    md = _set(metadata(lookup), metadata(newlookup))
    # Rebuild the new lookup with the merged fields
    rebuild(newlookup; data=parent(lookup), order=o, span=sp, sampling=sa, metadata=md)
end
_set(lookup::LookupArray, newlookup::NoLookup{<:AutoIndex}) = NoLookup(axes(lookup, 1))
_set(lookup::LookupArray, newlookup::NoLookup) = newlookup

# Set the index
_set(lookup::LookupArray, index::Val) = rebuild(lookup; data=index)
_set(lookup::LookupArray, index::Colon) = lookup
_set(lookup::LookupArray, index::AbstractArray) = rebuild(lookup; data=index)
_set(lookup::LookupArray, index::AutoIndex) = lookup
_set(lookup::LookupArray, index::AbstractRange) =
    rebuild(lookup; data=_set(parent(lookup), index), order=orderof(index))
# Update the Sampling lookup of Sampled dims - it must match the range.
_set(lookup::AbstractSampled, index::AbstractRange) = begin
    i = _set(parent(lookup), index)
    o = orderof(index)
    sp = Regular(step(index))
    rebuild(lookup; data=i, span=sp, order=o)
end

# Order
_set(lookup::LookupArray, neworder::Order) = rebuild(lookup; order=_set(order(lookup), neworder))
_set(lookup::NoLookup, neworder::Order) = lookup
_set(order::Order, neworder::Order) = neworder
_set(order::Order, neworder::AutoOrder) = order

# Span
_set(lookup::AbstractSampled, span::Span) = rebuild(lookup; span=span)
_set(lookup::AbstractSampled, span::AutoSpan) = lookup
_set(span::Span, newspan::Span) = newspan
_set(span::Span, newspan::AutoSpan) = span

# Sampling
_set(lookup::AbstractSampled, newsampling::Sampling) =
    rebuild(lookup; sampling=_set(sampling(lookup), newsampling))
_set(lookup::AbstractSampled, sampling::AutoSampling) = lookup
_set(sampling::Sampling, newsampling::Sampling) = newsampling
_set(sampling::Sampling, newsampling::AutoSampling) = sampling

# Locus
_set(lookup::AbstractSampled, locus::Locus) =
    rebuild(lookup; sampling=_set(sampling(lookup), locus))
_set(sampling::Points, locus::Union{AutoLocus,Center}) = Points()
_set(sampling::Points, locus::Locus) = _locuserror()
_set(sampling::Intervals, locus::Locus) = Intervals(locus)
_set(sampling::Intervals, locus::AutoLocus) = sampling

# Metadata
_set(lookup::LookupArray, newmetadata::AllMetadata) = rebuild(lookup; metadata=newmetadata)
_set(metadata::AllMetadata, newmetadata::AllMetadata) = newmetadata

# Index
_set(index::AbstractArray, newindex::AbstractArray) = newindex
_set(index::AbstractArray, newindex::AutoLookup) = index
_set(index::Colon, newindex::AbstractArray) = newindex

_set(A, x) = _cantseterror(A, x)

@noinline _locuserror() = throw(ArgumentError("Can't set a locus for `Points` sampling other than `Center` - the index values are the exact points"))
@noinline _cantseterror(a, b) = throw(ArgumentError("Can not set any fields of $(typeof(a)) to $(typeof(b))"))


text:
using PosDefManifold, PosDefManifoldML, Statistics, Test

# Give two symmetric positive definite matrices A and B,
# 80 marices will be created randomly as points on the the geodesic
# relying A to B and 80 marices as points on the the geodesic relying B to A,
# for 80 arc-lengths given by function rand(Float64)*exp(-i/(160))),
# for 1=1:80. This makes 80 two class problems with 160 observations each
# that must be progressively better separable as i increases.
# For each of three default machine learning models (MDM, ENLR and SVM)
# the average accuracy is found by a 10-fold cross-validation for each i
# dataset. The Pearson correlation of the accuracy with the vector [1 2...80]
# is computed and a statistical test on the probability that the correlation is
# equal to zero is performed, requesting the p-value be <1e-8
# (i.e., the correlation with 80-2 degrees of freedom must be > 0.6)

# This is not a formal test, but a general test on the functioning of all
# main functions of the packages and on the expected behavior of ML models.

A= ℍ([1.98 0.97 0.17; 0.96 1.45 1.37; 0.17 1.37 6.33])
B= ℍ([1.21 2.82 -1.48; 2.82  7.09 -3.69; -1.48 -3.69  2.00])
n=80

g(α::Float64)=geodesic(Fisher, A, B, α)

gVec(n, i, P, Q)=
   ℍVector([geodesic(Fisher, P, Q, rand(Float64)*exp(-i/(n*2))) for j=1:n])

PTr=ℍVector₂([vcat(gVec(n, i, A, B), gVec(n, i, B, A)) for i=1:n])
PTe=ℍVector₂([vcat(gVec(n, i, A, B), gVec(n, i, B, A)) for i=1:n])
y=vcat(fill(1, n), fill(2, n))

println("\nTesting MDM model...")
cv=[cvAcc(MDM(), PTr[i], y; verbose=false, scoring=:b).avgAcc for i=1:n]
@test cor(cv, [i for i=1:n]) > 0.6

println("\nTesting ENLR model...")
cv=[cvAcc(ENLR(), PTr[i], y; verbose=false).avgAcc for i=1:n]
@test cor(cv, [i for i=1:n]) > 0.6

println("\nTesting SVM model...")
cv=[cvAcc(SVM(), PTr[i], y; verbose=false).avgAcc for i=1:n]
@test cor(cv, [i for i=1:n]) > 0.6

ext: jl


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PlaneElasticity.jl/blob/master/LICENSE

using Documenter, PlaneElasticity

makedocs(modules=[PlaneElasticity],
         format = :html,
         checkdocs = :all,
         sitename = "PlaneElasticity.jl",
         pages = ["index.md"]
        )


text:
reload("SAMC")
reload("MCBN")
using Base.Test

include("entropy.jl")
include("nets.jl")
include("cod.jl")
# include("accuracy.jl")
# include("zero.jl")


text:
module TestPkgAPI

using DocumentationOverview: DocumentationOverview, API
using Test

module ExampleModule

"""
    ExampleModule.f(x)

Example docstring.
"""
f(x) = nothing

nodoc(x) = nothing

end  # module ExampleModule

function test_api_from_str()
    api = API("ExampleModule.f", namespace = TestPkgAPI)
    @test api.module === ExampleModule
    @test api.name === :f
    @test api.value === ExampleModule.f
    @test api.maybevalue === Some(ExampleModule.f)

    nonexisting = API("ExampleModule.nonexisting", namespace = TestPkgAPI)
    @test nonexisting.maybevalue === nothing
end

end  # module


text:
""" Algorithms operating on wiring diagrams.
"""
module WiringDiagramAlgorithms
export Junction, add_junctions!, rem_junctions!,
  normalize_cartesian!, normalize_copy!, normalize_delete!,
  topological_sort, crossing_minimization_by_sort

using AutoHashEquals
import LightGraphs
using UnionFind
using Statistics: mean

using ..WiringDiagramCore, ..WiringLayers
import ..WiringDiagramCore: input_ports, output_ports, set_box

# Traversal
###########

""" Topological sort of boxes in wiring diagram.

Returns a list of box IDs, excluding the outer box's input and output IDs.
"""
function topological_sort(d::WiringDiagram)::Vector{Int}
  vs = LightGraphs.topological_sort_by_dfs(graph(d))
  outer_ids = (input_id(d), output_id(d))
  filter(v -> v ∉ outer_ids, vs)
end

# Diagonals and codiagonals
###########################

""" Junction node in a wiring diagram.

Used to explicitly represent copies, merges, deletions, and creations.
"""
@auto_hash_equals struct Junction{Value} <: AbstractBox
  value::Value
  ninputs::Int
  noutputs::Int
end
input_ports(junction::Junction) = repeat([junction.value], junction.ninputs)
output_ports(junction::Junction) = repeat([junction.value], junction.noutputs)

""" Add junction nodes to wiring diagram.

Transforms from implicit to explicit representation of (co)diagonals.
"""
function add_junctions!(d::WiringDiagram)
  add_output_junctions!(d, input_id(d))
  add_input_junctions!(d, output_id(d))
  for v in box_ids(d)
    add_input_junctions!(d, v)
    add_output_junctions!(d, v)
  end
  return d
end
function add_input_junctions!(d::WiringDiagram, v::Int)
  for (port, port_value) in enumerate(input_ports(d, v))
    wires = in_wires(d, v, port)
    nwires = length(wires)
    if nwires != 1
      rem_wires!(d, wires)
      jv = add_box!(d, Junction(port_value, nwires, 1))
      add_wire!(d, Port(jv, OutputPort, 1) => Port(v, InputPort, port))
      add_wires!(d, [ wire.source => Port(jv, InputPort, i)
                      for (i, wire) in enumerate(wires) ])
    end
  end
end
function add_output_junctions!(d::WiringDiagram, v::Int)
  for (port, port_value) in enumerate(output_ports(d, v))
    wires = out_wires(d, v, port)
    nwires = length(wires)
    if nwires != 1
      rem_wires!(d, wires)
      jv = add_box!(d, Junction(port_value, 1, nwires))
      add_wire!(d, Port(v, OutputPort, port) => Port(jv, InputPort, 1))
      add_wires!(d, [ Port(jv, OutputPort, i) => wire.target
                      for (i, wire) in enumerate(wires) ])
    end
  end
end

""" Remove junction nodes from wiring diagram.

Transforms from explicit to implicit representation of (co)diagonals.
"""
function rem_junctions!(d::WiringDiagram)
  junction_ids = filter(v -> box(d,v) isa Junction, box_ids(d))
  junction_diagrams = map(junction_ids) do v
    junction = box(d,v)::Junction
    layer = complete_layer(junction.ninputs, junction.noutputs)
    to_wiring_diagram(layer, input_ports(junction), output_ports(junction))
  end
  substitute!(d, junction_ids, junction_diagrams)
  return d
end

""" Put a wiring diagram for a cartesian category into normal form.

This function puts a wiring diagram representing a morphism in a free cartesian
category into normal form. Copies and deletions are simplified as much as
possible.
"""
function normalize_cartesian!(d::WiringDiagram)
  # TODO: By duality, the same algorithm will normalize a cocartesian diagram.
  # We should support cartesian, cocartesian, and bicartesian diagrams.
  normalize_copy!(normalize_delete!(d))
end

""" Normalize copies in a wiring diagram.

This function maximizes sharing of intermediate computations in a wiring
diagram where copies are natural.

This algorithm is basically the same as the congruence closure algorithm on
term graphs, in the special case of the empty relation R = ∅
(Baader & Nipkow, 1998, Term Rewriting and All That, Sec. 4.4).
The main difference is the possibility of zero or many function outputs.
"""
function normalize_copy!(d::WiringDiagram)
  # Compute equivalence classes of boxes (without modifying the diagram).
  uf = UnionFinder(nboxes(d)+2)
  initial = filter(box_ids(d)) do v
    all(u == input_id(d) for u in inneighbors(d,v))
  end
  for v1 in initial
    for v2 in initial
      merge_if_congruent!(d, uf, v1, v2)
    end
  end

  # Keep only the designated representative of each equivalence class.
  extra = Int[]
  for v in box_ids(d)
    rep = find!(uf, v)
    if v != rep
      for wire in out_wires(d, v)
        add_wire!(d, Wire(set_box(wire.source, rep), wire.target))
      end
      push!(extra, v)
    end
  end
  rem_boxes!(d, extra)
  d
end

function merge_if_congruent!(d::WiringDiagram, uf::UnionFinder, v1::Int, v2::Int)
  if v1 == v2 || (find!(uf, v1) != find!(uf, v2) && is_congruent(d, uf, v1, v2))
    union!(uf, v1, v2)
    for out1 in filter(v -> v != output_id(d), outneighbors(d, v1))
      for out2 in filter(v -> v != output_id(d), outneighbors(d, v2))
        merge_if_congruent!(d, uf, out1, out2)
      end
    end
  end
end

function is_congruent(d::WiringDiagram, uf::UnionFinder, v1::Int, v2::Int)::Bool
  box(d, v1) == box(d, v2) && all(eachindex(input_ports(box(d,v1)))) do port
    wires1, wires2 = in_wires(d,v1,port), in_wires(d,v2,port)
    n1, n2 = length(wires1), length(wires2)
    @assert n1 <= 1 && n2 <= 1 # TODO: Handle merges?
    n1 == n2 && all(zip(wires1, wires2)) do pair
      src1, src2 = pair[1].source, pair[2].source
      set_box(src1, find!(uf, src1.box)) == set_box(src2, find!(uf, src2.box))
    end
  end
end

""" Normalize deletions in a wiring diagram.

This function removes all unused intermediate computations in a wiring diagram
where deletion is natural.
"""
function normalize_delete!(d::WiringDiagram)
  unused = Set{Int}()
  for v in reverse(topological_sort(d))
    if all(wire.target.box in unused for wire in out_wires(d, v))
      push!(unused, v)
    end
  end
  rem_boxes!(d, unused)
  d
end

# Layout
########

""" Crossing minimization by sorting a univariate statistic.

The boxes in `sources` and/or `targets` are fixed and the boxes in `vs` are
permuted. A permutation `σ` of the latter is returned, such that `vs[σ]` are the
sorted box IDs. Both one-sided and two-sided crossing minimization are
supported, depending on whether just one, or both, of `sources` and `targets`
are given.

In this simple but popular heuristic algorithm, the boxes are permuted by
sorting a univariate statistic of the positions of incoming and/or outgoing
wires. Typical choices are:

- `mean`: the sample mean, yielding the "barycenter method"
- `median`: the sample median

In both cases, this algorithm has the property that if there is a permutation
with no crossings, it will find it.
"""
function crossing_minimization_by_sort(d::WiringDiagram, vs::Vector{Int};
    sources::Vector{Int}=Int[], targets::Vector{Int}=Int[],
    statistic::Function=mean)::Vector{Int}
  @assert allunique(vs) && allunique(sources) && allunique(targets)
  if isempty(sources) && isempty(targets)
    # Degenerate case: nothing to sort, so preserve original order.
    return collect(eachindex(vs))
  end

  source_coord = port_coords(d, sources, OutputPort)
  target_coord = port_coords(d, targets, InputPort)
  stats = map(vs) do v
    source_coords = mapreduce(vcat, sources; init=[]) do source
      [ source_coord(wire.source) for wire in wires(d, source, v) ]
    end
    target_coords = mapreduce(vcat, targets; init=[]) do target
      [ target_coord(wire.target) for wire in wires(d, v, target) ]
    end
    statistic(vcat(source_coords, target_coords))
  end
  sortperm(stats)
end

""" Make function mapping ports to logical coordinates.
"""
function port_coords(d::WiringDiagram, vs::Vector{Int}, kind::PortKind)
  get_ports = kind == InputPort ? input_ports : output_ports
  index = Dict(vs[i] => i for i in eachindex(vs))
  sizes = [ length(get_ports(d,v)) for v in vs ]
  offsets = cumsum(vcat([0], sizes))
  (port::Port -> offsets[index[port.box]] + port.port)
end

end


text:
t = Template(; user=me)
pkg_dir = joinpath(t.dir, test_pkg)

@testset "CITATION" begin
    @testset "Plugin creation" begin
        p = Citation()
        @test isempty(p.gitignore)
        @test p.dest == "CITATION.bib"
        @test isempty(p.badges)
        @test isempty(p.view)
        @test !p.readme_section
        p = Citation(; readme_section=true)
        @test p.readme_section
    end

    @testset "File generation" begin
        p = Citation()
        @test gen_plugin(p, t, test_pkg) == ["CITATION.bib"]
        @test isfile(joinpath(pkg_dir, "CITATION.bib"))
        citation = read(joinpath(pkg_dir, "CITATION.bib"), String)

        @test occursin("@misc", citation)
        @test occursin("$(t.authors)", citation)
        @test occursin("v0.1.0", citation)
    end

    @testset "Readme untouched" begin
        p = Citation(; readme_section=false)
        t.plugins[Citation] = p
        isdir(pkg_dir) && rm(pkg_dir; recursive=true)
        generate(test_pkg, t, git=false)
        readme = read(joinpath(pkg_dir, "README.md"), String)
        @test !occursin("## Citing", readme)
        @test !occursin("CITATION.bib", readme)
    end

    @testset "Readme modification" begin
        p = Citation(; readme_section=true)
        t.plugins[Citation] = p
        isdir(pkg_dir) && rm(pkg_dir; recursive=true)
        generate(test_pkg, t, git=false)
        readme = read(joinpath(pkg_dir, "README.md"), String)
        @test occursin("## Citing", readme)
        @test occursin("CITATION.bib", readme)
    end
end

rm(pkg_dir; recursive=true)


text:
export eachneighborindex, nnodes

const IndexSetNetwork = Vector{Dict{Int,IndexSet}}

IndexSetNetwork(n::Integer) = [eltype(IndexSetNetwork)() for _ in 1:n]

struct ITensorNetwork
  itensors::Vector{ITensor}
  indsnetwork::IndexSetNetwork
end

Base.length(tn::ITensorNetwork) = length(itensors(tn))
Base.copy(tn::ITensorNetwork) = ITensorNetwork(copy(itensors(tn)), copy.(inds(tn)))

# Call this with tn[1:end]
itensors(tn::ITensorNetwork) = tn.itensors

ITensors.inds(tn::ITensorNetwork) = tn.indsnetwork

# TODO: should this return an ITensorNetwork with 1 node?
Base.getindex(tn::ITensorNetwork, i::Integer) = itensors(tn)[i]

# TODO: should this update reverse(ij)?
function Base.setindex!(isn::IndexSetNetwork, is::IndexSet, ij::Pair{<:Integer,<:Integer})
  i, j = ij
  isn[i][j] = is
  #isn[j][i] = dag(is)
  return isn
end

# TODO: also update indsnetwork with the new indices
Base.setindex!(tn::ITensorNetwork, t::ITensor, i::Integer) = (itensors(tn)[i] = t)

function setinds!(tn::ITensorNetwork, is::IndexSet, ij::Pair{<:Integer,<:Integer})
  return inds(tn)[ij] = is
end

# Iteration iterates through the ITensors of the network
Base.iterate(tn::ITensorNetwork, args...) = iterate(itensors(tn), args...)

# Iterate over the nodes that are neighbors of the specified node
# For example, for an ITensorNetwork `tn = [A₁, A₂, A₃, A₄]`:
#
# --A₁--A₂--A₃--A₄--
#
# `collect(eachneighborindex(tn, 2)) == [1, 3]`
eachneighborindex(isn::IndexSetNetwork, i::Integer) = Iterators.filter(==(i), keys(isn[i]))

# Iterator over each IndexSet connected to the node `i`
eachlinkinds(isn::IndexSetNetwork, i::Integer) = values(isn[i])

eachneighborindex(tn::ITensorNetwork, i::Integer) = eachneighborindex(inds(tn))

function ITensorNetwork(tensors::Vector{<:ITensor})
  N = length(tensors)
  indsnetwork = [eltype(IndexSetNetwork)() for n in 1:N]
  # Get the "link" indices
  # Determine the neighborhood of each node/tensor
  # through common indices
  for i in 1:N
    for j in (i + 1):N
      linkindsᵢⱼ = commoninds(tensors[i], tensors[j])
      if !isempty(linkindsᵢⱼ)
        indsnetwork[i][j] = linkindsᵢⱼ
        indsnetwork[j][i] = dag(indsnetwork[i][j])
      end
    end
  end
  # Get the "site" indices
  for i in 1:N
    indsnetwork[i][i] = uniqueinds(tensors[i], eachlinkinds(indsnetwork, i)...)
  end
  return ITensorNetwork(tensors, indsnetwork)
end

function Base.getindex(tn::IndexSetNetwork, ij::Pair{<:Integer,<:Integer})
  i, j = ij
  @assert j ∈ 1:nnodes(tn)
  return get(tn[i], j, IndexSet())
end

ITensors.inds(tn::ITensorNetwork, ij::Pair{<:Integer,<:Integer}) = inds(tn)[ij]

function ITensors.linkinds(tn::ITensorNetwork, ij::Pair{<:Integer,<:Integer})
  i, j = ij
  @assert i ≠ j
  return inds(tn, ij)
end

# Get the self edges/site indices of node `i`
ITensors.siteinds(tn::ITensorNetwork, i::Integer) = inds(tn, i => i)

# Get all of the self edges/site indices of the network, returning `s` such
# that `s[i] == siteinds(tn, i)`
ITensors.siteinds(tn::ITensorNetwork) = IndexSet[siteinds(tn, i) for i in 1:nnodes(tn)]

# Number of nodes/vertices/tensors of the graph/network
nnodes(isn::IndexSetNetwork) = length(isn)
function nnodes(tn::ITensorNetwork)
  n = nnodes(inds(tn))
  @assert n == length(itensors(tn))
  return n
end

⊗(Ts::ITensor...) = ITensorNetwork([Ts...])

# TODO: build this by appending T to the network instead
# of rebuilding the entire network
⊗(Ts::ITensorNetwork, T::ITensor) = ITensorNetwork([Ts..., T])

⊗(TN1::ITensorNetwork, TN2::ITensorNetwork) = ITensorNetwork([TN1..., TN2...])

#
# ITensors.jl extensions
#

ITensors.unioninds(is::Vector{<:IndexSet}) = unioninds(is...)

Base.:*(tn::ITensorNetwork) = *(tn...)
Base.:*(t::ITensor, tn::ITensorNetwork) = *(t, tn...)
Base.:*(tn::ITensorNetwork, t::ITensor) = *(tn..., t)
Base.:*(tn1::ITensorNetwork, tn2::ITensorNetwork) = *(tn1..., tn2...)

function ITensors.dag(tn::ITensorNetwork)
  tn = copy(tn)
  N = nnodes(tn)
  for i in 1:N
    tn[i] = dag(tn[i])
    for j in keys(inds(tn)[i])
      setinds!(tn, dag(inds(tn, i => j)), i => j)
    end
  end
  return tn
end

#
# Contraction sequence optimization
#

function trivial_contraction_sequence(tn::ITensorNetwork)
  N = nnodes(tn)
  # N-1 pairwise contractions
  t = 1
  for n in 2:N
    t = [n, t]
  end
  return t
end

function contract_pair!(T⃗, ab)
  a, b = ab
  Tᵃ, Tᵇ = T⃗[a], T⃗[b]
  # Contract a and b
  Tᵈ = Any[Tᵃ, Tᵇ]
  # Remove tensors Tᵃ and Tᵇ from the list of tensors
  deleteat!(T⃗, [a, b])
  # Append Tᵈ to the list of tensors
  push!(T⃗, Tᵈ)
  return T⃗
end

# Get all tensor pairs for a given number of tensors
function tensor_pairs(n::Integer)
  pairs = Pair{Int,Int}[]
  for i in 1:(n - 1)
    for j in (i + 1):n
      push!(pairs, i => j)
    end
  end
  return pairs
end

# Depth-first constructive approach
# TODO: this isn't pruning repeats
function each_contraction_sequence(n::Integer)
  Ts = Any[]
  each_each_ab = Iterators.product(
    (ITensorInfiniteMPS.tensor_pairs(n) for n in reverse(2:n))...
  )
  for each_ab in each_each_ab
    T = Any[1:n...]
    for ab in each_ab
      contract_pair!(T, ab)
    end
    push!(Ts, T[])
  end
  return Ts
end

each_contraction_sequence(tn::ITensorNetwork) = each_contraction_sequence(nnodes(tn))
function each_contraction_sequence(isn::Vector{<:Vector{<:Index}})
  return each_contraction_sequence(length(isn))
end
each_contraction_sequence(isn::IndexSetNetwork) = each_contraction_sequence(length(isn))

#
# Contraction cost using Vector{<: Vector{<: IndexSet}}
#

contraction_cost(T1::ITensor, T2::ITensor) = contraction_cost(inds(T1), inds(T2))

#
# TODO: use the following for symdiff
# Index[setdiff(inds(TN[1]), inds(TN[2]))..., setdiff(inds(TN[2]), inds(TN[1]))...]
#

function _intersect(A::Vector{IndexT}, B::Vector{IndexT}) where {IndexT<:Index}
  R = IndexT[]
  for a in A
    a ∈ B && push!(R, a)
  end
  return R
end

function _setdiff(A::Vector{IndexT}, B::Vector{IndexT}) where {IndexT<:Index}
  R = IndexT[]
  for a in A
    a ∉ B && push!(R, a)
  end
  return R
end

# Faster symdiff
function _symdiff(is1::Vector{IndexT}, is2::Vector{IndexT}) where {IndexT<:Index}
  setdiff12 = _setdiff(is1, is2)
  setdiff21 = _setdiff(is2, is1)
  display(setdiff12)
  display(setdiff21)
  return IndexT[_setdiff(is1, is2)..., _setdiff(is2, is1)...]
end

# Return the noncommon indices and the cost of contraction
function _contract_inds(is1::Vector{IndexT}, is2::Vector{IndexT}) where {IndexT<:Index}
  N1 = length(is1)
  N2 = length(is2)
  Ncommon = 0
  dim_common = 1
  is1_contracted = fill(false, N1)
  is2_contracted = fill(false, N2)
  # Determine the contracted indices
  for (n1, i1) in pairs(is1)
    n2 = findfirst(==(i1), is2)
    if !isnothing(n2)
      Ncommon += 1
      dim_common *= dim(i1)
      is1_contracted[n1] = true
      is2_contracted[n2] = true
    end
  end
  cost = _dim(is1) * _dim(is2) ÷ dim_common
  Nnoncommon = N1 + N2 - 2 * Ncommon
  isR = Vector{IndexT}(undef, Nnoncommon)
  n = 1
  for n1 in 1:N1
    if !is1_contracted[n1]
      isR[n] = is1[n1]
      n += 1
    end
  end
  for n2 in 1:N2
    if !is2_contracted[n2]
      isR[n] = is2[n2]
      n += 1
    end
  end
  return isR, cost
end

function _dim(is::Vector{<:Index})
  isempty(is) && return 1
  return mapreduce(dim, *, is)
end

function contraction_cost(is1::Vector{IndexT}, is2::Vector{IndexT}) where {IndexT<:Index}
  # 1.169417 s for N = 7 network
  #common_is = _intersect(is1, is2)
  #noncommon_is = _symdiff(is1, is2)
  #return noncommon_is, _dim(common_is) * _dim(noncommon_is)

  # 0.635962 s for N = 7 network
  return _contract_inds(is1, is2)
end

function contraction_cost!(
  cost::Ref{Int}, is1::Vector{IndexT}, is2::Vector{IndexT}
) where {IndexT<:Index}
  contracted_is, current_cost = contraction_cost(is1, is2)
  cost[] += current_cost
  return contracted_is
end

contraction_cost!(cost::Ref{Int}, ::Vector{Union{}}, ::Vector{Union{}}) = (cost[] += 1)

function contraction_cost!(cost::Ref{Int}, is::Vector{<:Vector{<:Index}}, n::Integer)
  return is[n]
end

function contraction_cost!(cost::Ref{Int}, is::Vector{<:Vector{<:Index}}, sequence)
  return contraction_cost!(
    cost, contraction_cost!(cost, is, sequence[1]), contraction_cost!(cost, is, sequence[2])
  )
end

function contraction_cost(isn::Vector{<:Vector{<:Index}}, sequence)
  # TODO: use the network itself
  cost = Ref(0)
  contraction_cost!(cost, isn, sequence)
  return cost[]
end

function contraction_cost(tn::ITensorNetwork, sequence)
  return contraction_cost(collect.(inds.(tn)), sequence)
end

# Return the optimal contraction sequence
function optimal_contraction_sequence(isn::Vector{<:Vector{<:Index}})
  sequences = each_contraction_sequence(isn)
  contraction_costs = map(sequence -> contraction_cost(isn, sequence), sequences)
  mincost, minsequenceindex = findmin(contraction_costs)
  return sequences[minsequenceindex], mincost
end

# Return the optimal contraction sequence
function optimal_contraction_sequence(tn::ITensorNetwork)
  if length(tn) == 1
    return Any[1], 0
  end
  #return optimal_contraction_sequence(inds(tn))
  #return optimal_contraction_sequence(inds.(itensors(tn)))
  return optimal_contraction_sequence(collect.(inds.(itensors(tn))))
end


text:
using DataStructures: nlargest

@doc raw"""
    GlobalPool(aggr)

Global pooling layer for graph neural networks.
Takes a graph and feature nodes as inputs
and performs the operation

```math
\mathbf{u}_V = \square_{i \in V} \mathbf{x}_i
```
where ``V`` is the set of nodes of the input graph and
the type of aggregation represented by ``\square`` is selected by the `aggr` argument.
Commonly used aggregations are `mean`, `max`, and `+`.

```julia
using Flux, GraphNeuralNetworks, LightGraphs

pool = GlobalPool(mean)

g = GNNGraph(erdos_renyi(10, 4))
X = rand(32, 10)
pool(g, X) # => 32x1 matrix

g = Flux.batch([GNNGraph(erdos_renyi(10, 4)) for _ in 1:5])
X = rand(32, 50)
pool(g, X) # => 32x5 matrix
```
"""
struct GlobalPool{F} <: GNNLayer
    aggr::F
end

function (l::GlobalPool)(g::GNNGraph, X::AbstractArray)
    if isnothing(g.graph_indicator)
        # assume only one graph
        indexes = fill!(similar(X, Int, g.num_nodes), 1)
    else
        indexes = g.graph_indicator
    end
    return NNlib.scatter(l.aggr, X, indexes)
end

"""
    TopKPool(adj, k, in_channel)

Top-k pooling layer.

# Arguments

- `adj`: Adjacency matrix  of a graph.
- `k`: Top-k nodes are selected to pool together.
- `in_channel`: The dimension of input channel.
"""
struct TopKPool{T,S}
    A::AbstractMatrix{T}
    k::Int
    p::AbstractVector{S}
    Ã::AbstractMatrix{T}
end

function TopKPool(adj::AbstractMatrix, k::Int, in_channel::Int; init=glorot_uniform)
    TopKPool(adj, k, init(in_channel), similar(adj, k, k))
end

function (t::TopKPool)(X::AbstractArray)
    y = t.p' * X / norm(t.p)
    idx = topk_index(y, t.k)
    t.Ã .= view(t.A, idx, idx)
    X_ = view(X, :, idx) .* σ.(view(y, idx)')
    return X_
end

function topk_index(y::AbstractVector, k::Int)
    v = nlargest(k, y)
    return collect(1:length(y))[y .>= v[end]]
end

topk_index(y::Adjoint, k::Int) = topk_index(y', k)


text:
using Luxor, Test

function rotate_points_test(fname)
    Drawing(300, 300, fname)
    origin()
    background("black")
    pts = ngon(O, 120, 5, vertices=true)

    sethue("red")
    circle.(pts, 10, :stroke)

    # rotate each point by 2π
    pts1 = rotate_point_around_point.(pts, O, 2π)

    sethue("blue")
    circle.(pts1, 20, :stroke)

    # they should be equal
    for pr in zip(pts, pts1)
        @test first(pr) ≈ last(pr)
    end

    # now rotate each point by π
    pts2 = rotate_point_around_point.(pts, O, π)
    sethue("green")
    circle.(pts2, 30, :stroke)

    # they should not be equal
    for pr in zip(pts, pts2)
        @test !isapprox(first(pr), last(pr))
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
end

function tick_test(fname)

    Drawing(400, 400, fname)
    origin()
    background("antiquewhite")
    sethue("red")
    majorticks, minorticks = tickline(O - (100, 0), O + (100, 0), minor=1)

    @test length(majorticks) == 3
    @test length(minorticks) == 5

    @test majorticks[1] == Point(-100, 0)
    @test majorticks[2] == Point(0, 0)
    @test majorticks[3] == Point(100, 0)

    # all minor ticks should be 5.0 apart
    @test all(pt -> isapprox(pt, Point(50.0, 0.0)), diff(minorticks))

    @test finish() == true
    println("...finished test: output in $(fname)")
end

rotate_points_test("rotatepoints.png")
tick_test("ticktest.png")


text:
using Test
using SparseIR

using LinearAlgebra: I
using QuadGK: gauss

function validategauss(rule)
    @test rule.a ≤ rule.b
    @test all(≤(rule.b), rule.x)
    @test all(≥(rule.a), rule.x)
    @test issorted(rule.x)
    @test length(rule.x) == length(rule.w)
    @test rule.x_forward ≈ rule.x .- rule.a
    @test rule.x_backward ≈ rule.b .- rule.x
end

@testset "gauss.jl" begin
    @testset "collocate" begin
        r = SparseIR.legendre(20)
        cmat = SparseIR.legendre_collocation(r)
        emat = SparseIR.legvander(r.x, length(r.x) - 1)
        @test emat * cmat ≈ I(20)
    end

    @testset "gauss legendre" begin
        rule = SparseIR.legendre(200)
        validategauss(rule)
        x, w = gauss(200)
        @test rule.x ≈ x
        @test rule.w ≈ w
    end

    @testset "piecewise" begin
        edges = [-4, -1, 1, 3]
        rule = SparseIR.piecewise(SparseIR.legendre(20), edges)
        validategauss(rule)
    end

    @testset "integrals" begin
        r = SparseIR.legendre(40)
        f(x) = cos(5x) + cos(50x)
        @test SparseIR.quadrature(r, f) ≈ 2sin(5) / 5 + 2sin(50) / 50
        g(x) = sin(x^3)
        @test SparseIR.quadrature(SparseIR.reseat(r, 1, 4), g) ≈ 0.2042835844987353218
    end
end


text:
using Distance
using Test

#   Testes das distancias
@testset "Distances test" begin
    @test GetDistance(Euclidean(), [2,2], [2,4]) == 2
    @test GetDistance(CityBlock(), [2,2], [2,4]) == 2
    @test GetDistance(TotalVariation(), [2,2], [2,4]) == 1
    @test GetDistance(Chebyshev(), [2,2], [2,4]) == 2
    @test GetDistance(Jaccard(), [2,2], [2,4]) == 1 - 2/3
    @test GetDistance(BrayCurtis(), [2,2], [2,4]) == 1/5
    @test GetDistance(CosineDist(), [2,2], [2,4]) == 1 - 12/(4 * sqrt(10))
    @test GetDistance(SpanNormDist(), [2,2], [2,4]) == 2
    @test GetDistance(CityBlock(), [2,2], [2,4]) == 2
end

#   Testes de casos de erro
@testset "Errors test" begin
    @test_throws DimensionMismatch GetDistance(Euclidean(), [], [])
    @test_throws DimensionMismatch GetDistance(CityBlock(), [1, 3, 4], [1, 2])
end


text:
function gschur!(H::Hessenberg{Tq, S}, Z=Matrix(H.Q), alg=QRIteration(); kwargs...
                 ) where {S <: SymTridiagonal{T}} where {Tq, T <: Real}
    _gschur!(H.H, alg, Z; kwargs...)
    # aliasing might be a cruel trap, so don't.
    v = copy(H.H.dv)
    if Z === nothing
        return Schur(zeros(Tq,0,0),zeros(Tq,0,0),v)
    end
    # Schur form must have same type as Z
    n = length(v)
    Tschur = similar(Z,n,n)
    fill!(Tschur, zero(eltype(Z)))
    @inbounds for j=1:n
        Tschur[j,j] = v[j]
    end
    return Schur(Tschur, Z, v)
end

function geigen!(A::SymTridiagonal{T}, alg=QRIteration()) where T <: AbstractFloat
    n = length(A.dv)
    V = Matrix{T}(I,n,n)
    _gschur!(A, alg, V)
    λ = copy(A.dv)
    LinearAlgebra.Eigen(sorteig!(λ, V, eigsortby)...)
end

function geigvals!(A::SymTridiagonal{T}, alg=QRIteration()) where T <: AbstractFloat
    _gschur!(A, alg, nothing)
    λ = copy(A.dv)
    return λ
end

# the rest are internal implementation methods

function _gschur!(A::SymTridiagonal{T},
                  alg::Algorithm,
                  Z::Union{Nothing, AbstractArray} = nothing;
                  maxiter=30*size(A,1)) where {T}
    throw(ArgumentError("Unsupported algorithm $alg"))
end

#translated from LAPACK::zsteqr
# LAPACK Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
"""
     _gschur!(A::SymTridiagonal, ::QRIteration, Z)

    overwrites `A.dv` with eigenvalues, optionally applying orthogonal transforms to `Z`.
"""
function _gschur!(A::SymTridiagonal{T},
                  alg::QRIteration,
                  Z::Union{Nothing, AbstractArray} = nothing;
                  maxiter=30*size(A,1)) where {T}
    d = A.dv
    e = A.ev
    n = length(d)
    wantZ = (Z !== nothing) && (size(Z,1) > 0)
    if wantZ
        nz = size(Z,2)
        (nz == n) || throw(DimensionMismatch("second dimension of Z must match tridiagonal A"))
    end

    epsT = eps(T)
    eps2 = epsT^2
    safmin = floatmin(T)
    safmax = 1 / safmin
    ssfmax = sqrt(safmax) / 3
    ssfmin = sqrt(safmin) / eps2

    niter = 0
    gvec = Vector{Givens}(undef, wantZ ? n : 0)

    # Determine where matrix splits and choose QL or QR for each block,
    # according to whether top or bottom diagonal element is smaller.

    l1 = 1
    while l1 <= n
        if l1 > 1
            e[l1-1] = zero(T)
        end
        msplit = _st_findsplit_fwd!(d,e,l1,n)
        l = l1
        lsave = l
        lend = msplit
        lendsave = lend
        l1 = msplit + 1
        if lend == l
            @mydebug println("empty at $l, continuing")
            continue
        end
        # scale A[l:lend,l:lend]
        anorm = opnorm(SymTridiagonal(view(d, l:lend), view(e, l:lend-1)), Inf)
        iscale = 0
        if anorm == 0
            @mydebug println("null at $l:$lend, continuing")
            continue
        end
        if anorm > ssfmax
            iscale = 1
            safescale!(view(d, l:lend), anorm, ssfmax)
            safescale!(view(e, l:lend-1), anorm, ssfmax)
        elseif anorm < ssfmin
            iscale = 2
            safescale!(view(d, l:lend), anorm, ssfmin)
            safescale!(view(e, l:lend-1), anorm, ssfmin)
        end
        if abs(d[lend]) < abs(d[l])
            lend = lsave
            l = lendsave
        end

        if lend > l
            l, nit = _st_QLIter!(d,e,Z,l,lend,gvec)
        else
            l, nit = _st_QRIter!(d,e,Z,l,lend,gvec)
        end
        niter +- nit
        # undo scaling if necessary
        if iscale == 1
            safescale!(view(d, lsave:lendsave), ssfmax, anorm)
            safescale!(view(e, lsave:lendsave-1), ssfmax, anorm)
        elseif iscale == 2
            safescale!(view(d, lsave:lendsave), ssfmin, anorm)
            safescale!(view(e, lsave:lendsave-1), ssfmin, anorm)
        end

        # check for iteration limit
        if niter >= maxiter
            nbad = count(x -> !iszero(x), e)
            @warn "$nbad eigenvalues failed to converge in $maxiter iterations"
            return
        end
        @mydebug begin
            if l1 <= n
                println("iscale was $iscale; next: l1=$l1")
            else
                println("iscale was $iscale; done.")
            end
        end
    end # outer loop
    nothing
end

# QL iteration
function _st_QLIter!(d::AbstractVector{T},e,Z,l,lend,gvec) where {T}
    @mydebug println("QL $l:$lend")
    wantZ = Z !== nothing
    safmin = floatmin(T)
    eps2 = eps(T)^2
    niter = 0
    while l <= lend
        # find small subdiagonal
        msplit = lend
        for m=l:lend-1
            if abs2(e[m]) <= (eps2 * abs(d[m])) * abs(d[m+1]) + safmin
                msplit = m
                break
            end
        end
        if msplit < lend
            e[msplit] = zero(T)
        end
        p = d[l]
        if msplit != l
            # If remaining matrix is 2x2, diagonalize it
            if msplit == l+1
                @mydebug println("deflating 2 at $l")
                if wantZ
                    rt1, rt2, cs, sn = _dgz_sym2x2(d[l],e[l],d[l+1], true)
                    G = Givens(l,l+1,cs,sn)
                    gvec[l] = G
                    rmul!(Z, G')
                else
                    rt1, rt2 = _dgz_sym2x2(d[l],e[l],d[l+1], false)
                end
                d[l] = rt1
                d[l+1] = rt2
                e[l] = zero(T)
                l += 2
                if l <= lend
                    continue # to top of QL loop
                end
                break
            end
            niter += 1
            # compute shift
            g = (d[l+1] - p) / (2 * e[l])
            r = hypot(g, one(T))
            g = d[msplit] - p + (e[l] / (g + copysign(r,g)))
            s = one(T)
            c = one(T)
            p = zero(T)
            # inner loop
            for i = msplit-1:-1:l
                f = s * e[i]
                b = c * e[i]
                c,s,r = givensAlgorithm(g,f)
                if i < msplit-1
                    e[i+1] = r
                end
                g = d[i+1] - p
                r = (d[i] - g) * s + 2 * c * b
                p = s * r
                d[i+1] = g+p
                g = c*r - b
                if wantZ
                    gvec[i] = Givens(i,i+1,c, -s)
                end
            end
            if wantZ
                # xlasr RVB
                for i=msplit-1:-1:l
                    rmul!(Z, gvec[i]')
                end
            end
            d[l] -= p
            e[l] = g
        else
            d[l] = p
            l += 1
        end
    end # l loop for QL
    @mydebug println("QL end")
    return l, niter
end

# QR iteration
function _st_QRIter!(d::AbstractVector{T},e,Z,l,lend,gvec) where {T}
    @mydebug println("QR $lend:$l")
    wantZ = Z !== nothing
    safmin = floatmin(T)
    eps2 = eps(T)^2
    niter = 0
    while l >= lend
        msplit = lend
        for m = l:-1:lend+1
            if abs2(e[m-1]) <= (eps2 * abs(d[m])) * abs(d[m-1]) + safmin
                msplit = m
                break
            end
        end
        if msplit > lend
            e[msplit-1] = zero(T)
        end
        p = d[l]
        if msplit != l
            # If remaining matrix is 2x2, diagonalize it
            if msplit == l-1
                @mydebug println("deflating 2 at $msplit")
                if wantZ
                    rt1, rt2, cs,sn = _dgz_sym2x2(d[l-1],e[l-1],d[l], true)
                    G = Givens(l-1,l,cs,sn)
                    gvec[msplit] = G
                    # xlasr RVF
                    rmul!(Z, G')
                else
                    rt1, rt2, _, _ = _dgz_sym2x2(d[l-1],e[l-1],d[l], false)
                end
                d[l-1] = rt1
                d[l] = rt2
                e[l-1] = zero(T)
                l -= 2
                if l >= lend
                    continue # to top of QR loop
                end
                break
            end
            niter += 1
            # compute shift
            g = (d[l-1] - p) / (2 * e[l-1])
            r = hypot(g, one(T))
            g = d[msplit] - p + (e[l-1] / (g + copysign(r,g)))
            s = one(T)
            c = one(T)
            p = zero(T)
            # inner loop
            for i = msplit:l-1
                f = s * e[i]
                b = c * e[i]
                c,s,r = givensAlgorithm(g,f)
                if i != msplit
                    e[i-1] = r
                end
                g = d[i] - p
                r = (d[i+1] - g) * s + 2 * c * b
                p = s * r
                d[i] = g+p
                g = c*r - b
                if wantZ
                    gvec[i] = Givens(i,i+1,c, s)
                end
            end
            if wantZ
                # xlasr RVF
                for i=msplit:l-1
                    rmul!(Z, gvec[i]')
                end
            end
            d[l] -= p
            e[l-1] = g
        else
            d[l] = p
            l -= 1
        end
    end # l loop for QR
    @mydebug println("QR end")
    return l, niter
end

@inline function _st_findsplit_fwd!(d::AbstractVector{T},e,l1,n) where T
    msplit = n
    if l1 <= n-1
        for m=l1:n-1
            tst = abs(e[m])
            if iszero(tst)
                msplit = m
                break
            end
            if tst <= sqrt(abs(d[m])) * sqrt(abs(d[m+1])) * eps(T)
                e[m] = zero(T)
                msplit = m
                break
            end
        end
    end
    return msplit
end

function _dgz_sym2x2(a::T, b, c, wantvec) where T
    sm = a+c
    df = a-c
    adf = abs(df)
    tb = 2b
    ab = abs(tb)
    (acmx, acmn) = (abs(a) > abs(c)) ? (a,c) : (c,a)
    rt = hypot(ab,adf)
    if sm < zero(T)
        rt1 = (sm - rt) / 2
        sgn1 = -1
        # order of execution important
        # for full accuracy, next line needs higher precision
        rt2 = (acmx / rt1) * acmn - (b / rt1) * b
    elseif sm > zero(T)
        rt1 = (sm + rt) / 2
        sgn1 = 1
        # order of execution important
        # for full accuracy, next line needs higher precision
        rt2 = (acmx / rt1) * acmn - (b / rt1) * b
    else
        rt1, rt2 = rt/2, -rt/2
        sgn1 = 1
    end
    if wantvec
        # compute eigenvector
        if df > 0
            cs = df + rt
            sgn2 = 1
        else
            cs = df - rt
            sgn2 = -1
        end
        acs = abs(cs)
        if acs > ab
            ct = -tb / cs
            sn1 = one(T) / sqrt(one(T) + ct^2)
            cs1 = ct * sn1
        else
            if iszero(ab)
                cs1 = one(T)
                sn1 = zero(T)
            else
                tn = -cs / tb
                cs1 = one(T) / sqrt(one(T) + tn^2)
                sn1 = tn * cs1
            end
        end
        if sgn1 == sgn2
            tn = cs1
            cs1 = -sn1
            sn1 = tn
        end
        return rt1, rt2, cs1, sn1
    else
        return rt1, rt2, zero(T), zero(T)
    end
end


text:
dv = DataArray([1, 2, 3], bitpack([false, false, false]))

dv = DataArray([1, 2, 3], [false, false, false])

failNA(dv)
removeNA(dv)
replaceNA(dv, 3)
for v in each_failNA(dv)
	println(v)
end
for v in each_removeNA(dv)
	println(v)
end
for v in each_replaceNA(dv, 3)
	println(v)
end

dv[1] = NA

failNA(dv)
removeNA(dv)
replaceNA(dv, 3)
for v in each_failNA(dv)
	println(v)
end
for v in each_removeNA(dv)
	println(v)
end
for v in each_replaceNA(dv, 3)
	println(v)
end

dv = DataArray(ComplexPair, 5)

type MyType
	a::Int64
	b::Int64
end

dv = DataArray(MyType, 5)


text:
using AlphaStructures
using ViewerGL
GL =  ViewerGL

"""
Return points from file.
"""
function load_points(filename::String)::Matrix
    io = open(filename, "r")
    point = readlines(io)
    close(io)

    b = [tryparse.(Float64,split(point[i], " ")) for i in 1:length(point)]
    V = hcat(b...)
    return V
end

V = load_points("examples/examples3D/model/griffin.txt")

# view points
GL.VIEW([
	GL.GLPoints(permutedims(V))
])

# points by columns
filtration = AlphaStructures.alphaFilter(V);
VV, EV, FV, TV = AlphaStructures.alphaSimplex(V, filtration, 0.35)

#view all simplices
GL.VIEW(
	[
		GL.GLGrid(V, EV, GL.COLORS[1], 0.6) # White
		GL.GLGrid(V, FV, GL.COLORS[2], 0.8) # Red
		GL.GLGrid(V, TV, GL.COLORS[3], 0.6) # Green
	]
);


text:
using Random
using Flux: normalise
using StatsBase: countmap

function get_data(filename::String; split=true, dlm=",", delete_first_line=true)
    file = open(filename, "r")
    number_lines = countlines(file)
    close(file)
    lines = Vector{String}(undef, number_lines)
    i = 1
    for l in eachline(filename)
        lines[i] = l
        i += 1
    end
    if delete_first_line lines = lines[2:end] end
    if split lines = split_columns(lines; dlm=dlm) end
    lines
end

function split_columns(data::Vector{String}; dlm=",")::Vector{Vector{String}}
    data_length = length(data)
    columns_data = Vector{Vector{String}}(undef, data_length)
    for line = 1:data_length
        columns_data[line] = dlm == " " ? split(data[line]) : split(data[line], dlm)
    end
    columns_data
end

function separate_data(data::Vector{T} where T <: Any, test_percentage::Float64; shuffle=false)
    if shuffle shuffle!(data) end
    bound = Int(floor((1-test_percentage)*length(data)))
    train_set = data[1:bound]
    test_set = data[(bound+1):end]
    train_set, test_set
end

function convert_values(data::Vector{Vector{T}}, column::Int, conversion_table::Dict{String, <: Number}) where T <: Any
    data_length = length(data)
    if data_length >= 1 && length(data[1]) < column
        error("Column $column is out of data column bound.")
    end
    converted_data = Vector{Vector{Number}}(undef, data_length)
    num_columns  = length(data[1])
    for l = 1:data_length
        conv_dt = Vector{Number}(undef, num_columns)
        for c = 1:num_columns
            conv_dt[c] = c == column ? conversion_table[data[l][c]] : data[l][c]
        end
        converted_data[l] = conv_dt
    end
    converted_data
end

function convert_values(data::Vector{Vector{T}}, column::Int, conversion::Type{<: Number}) where T <: Any
    data_length = length(data)
    if data_length >= 1 && length(data[1]) < column
        error("Column $column is out of data column bound.")
    end
    converted_data = Vector{Vector{Any}}(undef, data_length)
    num_columns  = length(data[1])
    for l = 1:data_length
        conv_dt = Vector{Any}(undef, num_columns)
        for c = 1:num_columns
            conv_dt[c] = c == column ? parse(conversion, data[l][c]) : data[l][c]
        end
        converted_data[l] = conv_dt
    end
    converted_data
end

# Vector{Union{Type{<: Number}, Dict{String, <: Number}}}
function convert_values(data::Vector{Vector{String}}, conversion_vector::Vector{S} where S <: Any)
    data_length = length(data)
    converted_data = Vector{Vector{Float64}}(undef, data_length)
    num_columns  = length(data[1])
    if num_columns != length(conversion_vector)
        error("There is no conversion function to all columns.")
    end
    for l = 1:data_length
        conv_dt = Vector{Any}(undef, num_columns)
        for c = 1:num_columns
            conv_dt[c] = conversion(data[l][c], conversion_vector[c])
        end
        converted_data[l] = conv_dt
    end
    converted_data
end

conversion(data::String, table::Dict{String, <: Number}) = table[data]
conversion(data::String, type::Type{T}) where T <: Number = parse(type, data)

function separate_input_output(data::Vector{Vector{T}};output_columns::Int=1, normalised::Bool=false) where T <: Any
    data_length = length(data)
    num_columns = length(data[1])
    if output_columns >= data_length
        error("The quantity of output columns is equal or greater than the number of data columns.")
    end
    separator = num_columns-output_columns+1
    input_data = Vector{Vector{Float64}}(undef, data_length)
    output_data = Vector{Vector{Float64}}(undef, data_length)
    for i = 1:data_length
        input_data[i] = data[i][1:(separator-1)]
        output_data[i] = data[i][separator:end]
    end
    if normalised input_data = normalise.(input_data, dims=1) end
    input_data, output_data
end

#=
function separate_dataset(data::Vector{Vector{T}} where T <: Number; output_columns::Int=1,
                          test_percentage::Float64=0.2, shuffle::Bool=false,
                          normalise::Bool=false, separate_test_data::Bool=true,
                          balance_classes::Bool=false)
    #if normalise normalise_inputs!(data, output_columns) end
    if separate_test_data
        train_data, test_data = separate_data(data, test_percentage; shuffle=shuffle)
        if balance_classes balance_classes!(train_data) end
        input_train_data, output_train_data = separate_input_output(train_data, normalised=normalise)
        input_test_data, output_test_data = separate_input_output(test_data, normalised=normalise)
        return input_train_data, output_train_data, input_test_data, output_test_data
    else
        if balance_classes balance_classes!(data) end
        input_data, output_data = separate_input_output(data, normalised=normalise)
    end
end
=#

function separate_dataset(data::Vector{Vector{T}} where T <: Number; output_columns::Int=1,
                          test_percentage::Float64=0.2, shuffle::Bool=false,
                          normalise::Bool=false, separate_test_data::Bool=true,
                          balance_classes::Bool=false)
    if normalise normalise_inputs!(data, output_columns) end
    if separate_test_data
        train_data, test_data = separate_data(data, test_percentage; shuffle=shuffle)
        if balance_classes balance_classes!(train_data) end
        input_train_data, output_train_data = separate_input_output(train_data, normalised=false)
        input_test_data, output_test_data = separate_input_output(test_data, normalised=false)
        return input_train_data, output_train_data, input_test_data, output_test_data
    else
        if balance_classes balance_classes!(data) end
        input_data, output_data = separate_input_output(data, normalised=false)
    end
end

#!!!
function normalise_inputs!(data::Vector{Vector{T}} where T <: Number, output_columns::Int=1)
    number_features = length(data[1])
    normalised_features = Vector{Vector{Float64}}(undef, number_features)

    for feature in 1:number_features
        if feature <= (number_features - output_columns)
            normalised_features[feature] = Flux.normalise(getindex.(data, feature))
        else
            normalised_features[feature] = getindex.(data, feature)
        end
    end

    for i in 1:length(data)
        data[i] = getindex.(normalised_features, i)
    end
end

function create_numerical_map(data::Vector{T})::Dict{T, <: Any}  where {T <: Any}
    set_values = sort(collect(Set(data)))
    indexes = 1.0:length(set_values)
    Dict(collect(zip(set_values, indexes)))
end

function balance_classes!(data::Vector{Vector{T}} where T <: Number)
    # The class (output) is the last element in the vector
    outputs = map(v -> v[end], data)
    distribution = countmap(outputs)
    biggest_class = maximum(values(distribution))
    if any(c -> c[2] != biggest_class, distribution)
        for class in keys(distribution)
            if distribution[class] !=  biggest_class
                diff = biggest_class - distribution[class]
                class_inputs = filter(v -> v[end] == class, data)
                class_len = length(class_inputs)
                number_copies = floor(Int, diff/class_len)
                if number_copies != 0
                    append!(data, repeat(class_inputs, number_copies))
                end
                fill = diff - number_copies*class_len
                if fill > 0
                    append!(data, class_inputs[1:fill])
                end
            end
        end
    end
    data
end

abstract type DatasetSetup end

struct Dataset <: DatasetSetup
    input_train::Vector{Vector{Float64}}
    output_train::Vector{Vector{Float64}}
    labels::Vector{Float64}
    n_features::Int
    dataset_id::Int
end

struct SeparatedDataset <: DatasetSetup
    input_train::Vector{Vector{Float64}}
    output_train::Vector{Vector{Float64}}
    input_test::Vector{Vector{Float64}}
    output_test::Vector{Vector{Float64}}
    labels::Vector{Float64}
    n_features::Int
    dataset_id::Int
end

function Base.show(io::IO, dt::DatasetSetup)
    print(io, "Dataset -> id:$(dt.dataset_id), features:$(dt.n_features), labels:$(dt.labels)")
end

include("dataset_plot.jl")

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
module BinaryProvider

import Pkg3: iswindows, isapple, islinux

# Include our subprocess running funtionality
include("OutputCollector.jl")
# External utilities such as downloading/decompressing tarballs
include("PlatformEngines.jl")
# Platform naming
include("PlatformNames.jl")
# Everything related to file/path management
include("Prefix.jl")
# Abstraction of "needing" a file, that would trigger an install
include("Products.jl")
# Abstraction of bundled binary package
include("BinaryPackage.jl")

# BinDeps support, disabled for now because I don't particularly want to force
# users to install BinDeps to install this package.  That seems counter-productive
#include("BinDepsIntegration.jl")

function __init__()
    #global global_prefix

    # Initialize our global_prefix
    # global_prefix = Prefix(joinpath(dirname(@__FILE__), "../", "global_prefix"))
    # activate(global_prefix)

    # Find the right download/compression engines for this platform
    probe_platform_engines!()

    # If we're on a julia that's too old, then fixup the color mappings
    # if !haskey(Base.text_colors, :default)
    #     Base.text_colors[:default] = Base.color_normal
    # end
end

end # module


text:
function run_multiple_params(
      folder::String,
      objective::String,
      x_axis::Array,
      y_axis::Array,
      params,
      gen_areas=Int64[],
      area_transfer=Int64[],
      gen_zones=Int64[],
      zone_transfer=Int64[],
      vert_line=([],""),
      horiz_line=([],""),
      xlimits=[],
      ylimits=[],
      output_plot_label=("",""),
      scopf=false,
      ctg=0,
      points_params=(),
      suffix="";
      no_converter_loss::Bool=false,
      output_location_base="",
      output_results_folder=""
   )

   mn_data = read_sn_data(folder, no_converter_loss=no_converter_loss)

   folder = abspath(folder)
   if length(output_location_base) > 0
      output_folder = output_location_base
   else
      folder_split = splitpath(folder)
      toplevels = folder_split[1:end-3]
      output_folder = joinpath(toplevels...,"results/$(folder_split[end-1])/$(folder_split[end])$suffix")
   end
   output_folder = joinpath(output_folder, output_results_folder)
   if !isdir(output_folder)
      mkpath(output_folder)
   end
   println("output location: $output_folder")

   (results_dict, output_plot) = run_multiple_params(
         mn_data,
         output_folder,
         objective,
         x_axis,
         y_axis,
         params,
         gen_areas,
         area_transfer,
         gen_zones,
         zone_transfer,
         vert_line,
         horiz_line,
         xlimits,
         ylimits,
         output_plot_label,
         scopf,
         ctg,
         points_params;
         output_results_folder
      )

   return results_dict, output_plot
end

function run_multiple_params(
      mn_data::Dict,
      output_folder::String,
      objective::String,
      x_axis::Array,
      y_axis::Array,
      params,
      gen_areas=Int64[],
      area_transfer=Int64[],
      gen_zones=Int64[],
      zone_transfer=Int64[],
      vert_line=([],""),
      horiz_line=([],""),
      xlimits=[],
      ylimits=[],
      output_plot_label=("",""),
      scopf=false,
      ctg=0,
      points_params=();
      output_results_folder=""
   )

   output_folder = joinpath(output_folder, output_results_folder)

   (results_dict, subnet_arr,
   summary_df,
   binding_cnstr_fulldict,
   constraintlog,
   x_label_arr,
   constraint_transitions) = apply_params(
      mn_data,
      output_folder,
      objective,
      x_axis,
      y_axis,
      params,
      gen_areas,
      area_transfer,
      gen_zones,
      zone_transfer,
      vert_line,
      horiz_line,
      xlimits,
      ylimits,
      output_plot_label,
      scopf,
      ctg,
      false
   )
   if length(points_params) > 0
      (points_dict, _,
      pts_summary_df,
      _,_,_,_
      ) = apply_params(
         mn_data,
         output_folder,
         objective,
         x_axis,
         y_axis,
         points_params,
         gen_areas,
         area_transfer,
         gen_zones,
         zone_transfer,
         vert_line,
         horiz_line,
         xlimits,
         ylimits,
         output_plot_label,
         scopf,
         ctg,
         true
      )
      # concatenate summary dataframes
      summary_df = vcat(summary_df, pts_summary_df)
   end

   # println("final results_dict: $results_dict")
   results_dict_string = JSON.json(results_dict)
   println("Saving results at $(output_folder)/all_results_dict.json")
   open("$output_folder/all_results_dict.json", "w") do f
      write(f, results_dict_string)
   end

   println("Saving summary at $(output_folder)/all_param_summary.csv")
   CSV.write("$(output_folder)/all_param_summary.csv", summary_df)

   binding_cnstr_dict_string = JSON.json(binding_cnstr_fulldict)
   open("$output_folder/all_binding_constraints.json", "w") do f
      write(f, binding_cnstr_dict_string)
   end

   open("$output_folder/constraint_log.txt", "w") do f
      write(f, lstrip(constraintlog, '\n'))
   end

   x_label_arr=unique(x_label_arr)
   # println("x_label_arr: $x_label_arr")
   output_plot = nothing
   for x_label in x_label_arr
      # println("x_label: $x_label")
      xvals = []
      println("x_label: $x_label")
      if x_label[1] == "frequency (Hz)"
         xvals = constraint_transitions
         # println("xvals: $xvals")
      end
      if !(length(points_params) > 0)
         output_tmp = plot_results_dict_line(
            results_dict,
            subnet_arr,
            x_label,
            output_folder,
            vert_line,
            horiz_line,
            xlimits,
            ylimits,
            output_plot_label,
            xvals,
            plot_infeasible_boundaries=true
         )
      else
         output_tmp = plot_results_dict_line(
            results_dict,
            subnet_arr,
            x_label,
            output_folder,
            vert_line,
            horiz_line,
            xlimits,
            ylimits,
            output_plot_label,
            xvals,
            points_dict,
            plot_infeasible_boundaries=true
         )
      end

      if output_tmp != nothing
         output_plot = deepcopy(output_tmp)
      end
   end
   return results_dict, output_plot
end

function apply_params(
   folder::String,
   output_folder::String,
   objective::String,
   x_axis::Array,
   y_axis::Array,
   params,
   gen_areas,
   area_transfer,
   gen_zones,
   zone_transfer,
   vert_line,
   horiz_line,
   xlimits,
   ylimits,
   output_plot_label,
   scopf,
   ctg,
   pts_param=false;
   no_converter_loss::Bool=false
)

   mn_data = read_sn_data(folder, no_converter_loss=no_converter_loss)

   (
   results_dict,
   subnet_arr,
   summary_df,
   binding_cnstr_fulldict,
   constraintlog,
   x_label_arr,
   constraint_transitions
   ) = apply_params(
      mn_data,
      output_folder,
      objective,
      x_axis,
      y_axis,
      params,
      gen_areas,
      area_transfer,
      gen_zones,
      zone_transfer,
      vert_line,
      horiz_line,
      xlimits,
      ylimits,
      output_plot_label,
      scopf,
      ctg,
      pts_param
   )
   return results_dict, subnet_arr, summary_df, binding_cnstr_fulldict, constraintlog, x_label_arr, constraint_transitions
end

function apply_params(
   mn_data::Dict,
   output_folder::String,
   objective::String,
   x_axis::Array,
   y_axis::Array,
   params,
   gen_areas,
   area_transfer,
   gen_zones,
   zone_transfer,
   vert_line,
   horiz_line,
   xlimits,
   ylimits,
   output_plot_label,
   scopf,
   ctg,
   pts_param=false
)
   (p_keys,p_values)=params
   override_param = Dict{Any, Any}()
   results_dict = Dict{String, Dict{Any, Any}}()
   summary_df = DataFrame()
   subnet_arr = Array{Int64,1}()
   x_label_arr = []
   binding_cnstr_fulldict = Dict{String,Any}()
   prev_keystring = ""
   constraintlog = ""
   constraint_transitions = []
   for i in 1:length(p_values[1])
      # println("i=$i")
      for (val_idx,value) in enumerate(p_values)
         # println(val_idx)
         # println("keys: $(string.(p_keys[val_idx]))")
         # println("val: $(value[i])")
         set_nested!(override_param, string.(p_keys[val_idx]), value[i])
      end
      # println("running $(folder)")
      # println("override_param: $override_param")
      print_results = false
      fix_f_override = false
      direct_pq = true
      master_subnet = 1
      start_vals=Dict{String, Dict}("sn"=>Dict())
      uniform_gen_scaling = false
      unbounded_pg = false
      output_to_files = true
      (result, res_summary, solution_pm, binding_cnstr) = multifrequency_opf(
         mn_data,
         output_folder,
         objective,
         gen_areas,
         area_transfer,
         gen_zones,
         zone_transfer,
         print_results,
         override_param,
         fix_f_override,
         direct_pq,
         master_subnet,
         start_vals,
         uniform_gen_scaling,
         unbounded_pg,
         output_to_files
      )

      if (length(binding_cnstr_fulldict) == 0)
         # keystring = rstrip(join([p_keys[val_idx][end]*"=$(value[i]); " for (val_idx,value) in enumerate(p_values)]), ';')
         keystring = p_keys[1][end]*"=$(p_values[1][i])"
         binding_cnstr_fulldict[keystring] = binding_cnstr
         prev_keystring = keystring
      else
         if (res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]) && (binding_cnstr_fulldict[prev_keystring] != binding_cnstr)
            println("Active constraint set has changed. Saving set.")
            # keystring = rstrip(join([p_keys[val_idx][end]*"=$(value[i]); " for (val_idx,value) in enumerate(p_values)]), ';')
            keystring = p_keys[1][end]*"=$(p_values[1][i])"
            binding_cnstr_fulldict[keystring] = binding_cnstr
            constraintlog *= "\n$keystring\n"
            for (idx, class) in binding_cnstr
               for (subnet, constr_set) in class
                  prev_set = []
                  try
                     prev_set = binding_cnstr_fulldict[prev_keystring][idx][subnet]
                  catch KeyError
                     prev_set = []
                  end
                  constr_added = setdiff(constr_set, prev_set)
                  if length(constr_added) > 0
                     println("Added constraints in $subnet:")
                     println(constr_added)
                     constraintlog *= "Added constraints in $subnet:\n$constr_added\n"
                  end
               end
            end
            for (idx, class) in binding_cnstr_fulldict[prev_keystring]
               for (subnet, constr_set) in class
                  new_set = []
                  try
                     new_set = binding_cnstr[idx][subnet]
                  catch KeyError
                     new_set = []
                  end
                  constr_removed = setdiff(constr_set, new_set)
                  if length(constr_removed) > 0
                     println("Removed constraints in $subnet:")
                     println(constr_removed)
                     constraintlog *= "Removed constraints in $subnet:\n$constr_removed\n"
                  end
               end
            end
            # Add to contraint transition points if this is a frequency sweep
            key_idx = indexin(["f_min", "f_max"], [p_keys[val_idx][end] for (val_idx,value) in enumerate(p_values)])
            if sum(key_idx .!= nothing) == 2
               if abs(p_values[key_idx[1]][i] - p_values[key_idx[2]][i]) <= 1e-6
                  append!(constraint_transitions, p_values[key_idx[1]][i])
               end
            end
            if !pts_param
               # println("constraint_transitions: $constraint_transitions")
               open("$output_folder/constraint_log.txt", "w") do f
                  write(f, lstrip(constraintlog, '\n'))
               end
            end
            prev_keystring = keystring
         end
      end
      # println("Result keys: $(keys(result))")
      if i == 1
         subnet_arr = result["subnet"]
      end
      # println(override_param[:network_wide][:f_max][subnet])
      println("trial $i of $(length(p_values[1]))")
      # println(res_summary)
      res_summary[:,:network] .= "trial $i"
      # res_summary.network = "trial $i"
      for (ax_idx,axis) in enumerate([x_axis, y_axis])
         # ax_res_key = []
         # println("axis: $axis")
         for ax_key in axis
            if !(ax_key[1] in keys(result))
               if ax_key[1] == "sn"
                  for sn_key in ax_key[2]
                     elem = ax_key[3]
                     for elem_key in ax_key[4]
                        if isa(ax_key[5], Array)
                           for quantity in ax_key[5]
                              new_key = ("sn", "$sn_key", elem, "$elem_key", quantity)
                              label = join(new_key, " ")
                              if (ax_idx == 1) && !(new_key in x_label_arr)
                                 push!(x_label_arr,new_key)
                              end
                              value = get_nested(solution_pm, new_key)
                              if value == nothing
                                 # println(new_key)
                                 # println(keys(solution_pm))
                                 # println(solution_pm["sn"])
                                 # throw(KeyError("The value $new_key for axis $ax_idx is not found."))
                                 # println("The value $new_key for axis $ax_idx is not found.")
                                 value = NaN
                              end
                              if !(label in keys(results_dict))
                                 results_dict[label] = Dict{Any, Any}()
                              end
                              # Comment the '|| true' to prevent saving outputs when not solved to optimum
                              if !pts_param
                                 if res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED] || true
                                    results_dict[label][i] = value
                                 end
                              else
                                 if res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
                                    lbl = join(["$(p_keys[j][end])=$(round(p_values[j][i],digits=3)), " for j in 1:length(p_keys)])
                                    lbl = rstrip(lbl, ' ')
                                    lbl = rstrip(lbl, ',')
                                    println("lbl=$lbl")
                                    results_dict[label][lbl] = value
                                 end
                              end
                           end
                        else
                           quantity = ax_key[5]
                           new_key = ("sn", "$sn_key", elem, "$elem_key", quantity)
                           label = join(new_key, " ")
                           if (ax_idx == 1) && !(new_key in x_label_arr)
                              push!(x_label_arr,new_key)
                           end
                           value = get_nested(solution_pm, new_key)
                           if value == nothing
                              # println(new_key)
                              # println(keys(solution_pm))
                              throw(KeyError("The value $new_key for axis $ax_idx is not found."))
                           end
                           if !(label in keys(results_dict))
                              results_dict[label] = Dict{Float64, Any}()
                           end
                           # Comment '|| true' to prevent saving outputs when not solved to optimum
                           if !pts_param
                              if res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED] || true
                                 results_dict[label][i] = value
                              end
                           else
                              if res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
                                 lbl = join(["$(p_keys[j][end])=$(round(p_values[j][i],digits=3)), " for j in 1:length(p_keys)])
                                 lbl = rstrip(lbl, ' ')
                                 lbl = rstrip(lbl, ',')
                                 println("lbl=$lbl")
                                 results_dict[label][lbl] = value
                              end
                           end
                        end
                     end
                  end
               elseif (ax_idx == 1) && !(ax_key in x_label_arr) # ax_key[1] is not "sn": custom axis values and label
                  push!(x_label_arr,ax_key)
               end
            elseif (ax_idx == 1) && !(ax_key in x_label_arr) # ax_key is in the results dict and should be in the x axis array
               push!(x_label_arr,ax_key)
            end
         end
      end
      for (label,value) in result
         if !(label in keys(results_dict))
            results_dict[label] = Dict{Float64, Any}()
         end
         if !pts_param
            if res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED] || true
               results_dict[label][i] = value
            end
         else
            if res_summary.status[1] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
               lbl = join(["$(p_keys[j][end])=$(round(p_values[j][i],digits=3)), " for j in 1:length(p_keys)])
               lbl = rstrip(lbl, ' ')
               lbl = rstrip(lbl, ',')
               println("lbl=$lbl")
               results_dict[label][lbl] = value
            end
         end
      end

      if size(summary_df,1) > 0
         summary_df = vcat(summary_df, res_summary)
      else
         summary_df = res_summary
      end
   end
   return results_dict, subnet_arr, summary_df, binding_cnstr_fulldict, constraintlog, x_label_arr, constraint_transitions
end


text:
## There are many Poly class methods
## Instances of Poly can be created via:
## @syms x y z
## p = poly(x^2 -x - 1)
##

const SymPoly = Sym

Base.divrem(p::Sym, q::Sym) = sympy_meth(:div, p, q)
Base.div(p::Sym, q::Sym) = divrem(p,q)[1]
Base.rem(p::Sym, q::Sym) = divrem(p,q)[2]

## rename these, their use is special to polynomials, so we prefix
## Renamed polynomial methods.

## The functions `div`, `rem, `divrem` had slightly different meaning in
## Julia than SymPy, where they are polynomial methods. Here we rename
## them by prefixing with "poly". We do the same for `roots`, so as to
## not conflict with the `roots` function from `Polynomials.jl`.

## [This isn't really needed. We do keep the `polyroots`, as `roots` does not live in base.]

"""

Polynomial division. Renamed from `div` in SymPy to avoid confusion with Julia's `div`. Also, same as `divrem`.

"""
polydiv(ex::Sym, args...; kwargs...) = sympy_meth(:div, ex, args...; kwargs...)

"""

Polynomial division remainer.

"""

polyrem(ex::Sym, args...; kwargs...) = rem(ex, args...; kwargs...)

"""

Polynomial division with remainder. Renamed from `divrem` in SymPy to avoid confusion with Julia's `divrem`

"""
polydivrem(ex::Sym, args...; kwargs...) = sympy_meth(:div, ex, args...; kwargs...)

## deprecate these
Base.@deprecate polydivrem(p::Sym, q::Sym) divrem(p::Sym, q::Sym)
Base.@deprecate polydiv(p::Sym, q::Sym) divrem(p::Sym, q::Sym)
Base.@deprecate polyrem(p::Sym, q::Sym) rem(p::Sym, q::Sym)

"""

Find roots of a polynomial. Renamed from `roots` in
SymPy to avoid confusion with the `roots` function of `Polynomials`

"""
polyroots(ex::Sym, args...; kwargs...) = sympy_meth(:roots, ex, args...; kwargs...)
export polydiv, polyrem, polydivrem, polyroots

polynomial_sympy_methods = (
                            :sqf_list,
                            :groebner,
                            :solve_poly_system,
                            :resultant,
                            :cancel,
                            :apart, :together,
                            :poly,
                            :poly_from_expr,
                            :degree, :degree_list,
                            :LC, :LM, :LT,
                            :pdiv, :prem, :pquo, :pexquo,
                            :quo, :exquo,
                            :half_gcdex, :gcdex,
                            :invert,
                            :subresultants,
                            :discriminant,
                            :terms_gcd,
                            :cofactors,
                            :gcd_list,
                            :lcm_list,
                            :monic,
                            :content,
                            :compose, :decompose,
                            :sturm,
                            :gff_list, :gff,
                            :sqf_norm, :sqf_part, :sqf_list, :sqf,
                            :factor_list,
                            :intervals,
                            :refine_root,
                            :count_roots,
                            :real_roots,
                            :nroots,
                            :ground_roots,
                            :nth_power_roots_poly,
                            :reduced,
                            :is_zero_dimensional,
                            :symmetrize,
                            :horner,
                            :interpolate,
                            :viete,
                            :construct_domain,
                            :minimal_polynomial,
                            :minpoly, :primitive_element,
                            :field_isomorphism, :to_number_field
#                            :numer#, :denom
#                            :roots ## conflict with Roots.roots and functionality provided by solve
                            )

denom(x::Sym) = sympy_meth(:denom, x)
numer(x::Sym) = sympy_meth(:numer, x)
export denom, numer

polynomial_sympy_methods_import = (:expand,
                                   :factor #,
                                   #:trunc
                                   )

for meth in polynomial_sympy_methods_import
    meth_name = string(meth)
    @eval ($meth)(ex::Sym, args...; kwargs...) = sympy_meth($meth_name, ex, args...; kwargs...)
end

## SymPy Poly class

"""

Constructor for polynomials. This allows certain polynomial-specifec
methods to be called, such as `coeffs`. There is also `poly` that does
the same thing. (We keep `Poly`, as it is named like a constructor.)

Examples:
```
p = Poly(x^2 + 2x + 1, x)  #  a poly in x over Z
coeffs(p)  # [1, 2, 1]
p = Poly(a*x^2 + b*x + c, x)  # must specify variable that polynomial is over.
coeffs(p)  ## [a,b,c]
```

"""
Poly(x::Sym) = sympy[:Poly](x)
Poly(x::Sym, args...; kwargs...) = sympy_meth(:Poly, x, args...; kwargs...)
export Poly

## Poly class methods that aren't sympy methods
## e.g., must be called as obj.method(...), but not method(obj, ...)
## not coeffs(x^2 -1), rather poly(x^2 - 1)  |> coeffs
## or poly(a*x^2 + b*x + c, x) |> coeffs to specify a specific variable
polynomial_instance_methods = (:EC, :ET, :LC, :LM, :LT, :TC, ## no .abs()
                              :add_ground,
                              :all_coeffs,
                              :all_monoms, :all_roots, :all_terms,
                              :as_dict, :as_expr, :as_list,
                              :clear_denoms,
                              :coeffs,
                              :deflate,
                              :exclude,
                              :l1_norm,
                              :lift,
                              :ltrim,
                              :max_norm,
                              :monoms,
                              :mul_ground,
                              :neg,
                              :nth, :nth_power_roots_poly,
                              :per,
                              :rat_clear_denoms,
                              :retract,
                              :root,
                              :set_domain, :set_modulus,
                              :sqr, :sub_ground,
                              :to_exact, :to_field, :to_ring, :total_degree,
                              :unify
                              )

polynomial_predicates = (
                         :is_cyclotomic,
                         :is_ground,
                         :is_homogeneous,
                         :is_irreducible,
                         :is_linear,
                         :is_monic,
                         :is_monomial,
                         :is_multivariate,
                         :is_one,
                         :is_primitive,
                         :is_quadratic,
                         :is_sqf,
                         :is_univariate,
                         :is_zero)

## special cases due to output

## {Sym => val} dictionary
#roots(p::Sym, args...; kwargs...) = sympy.roots(project(p), args...; kwargs...)


text:
export ModelKitSystem

"""
    ModelKitSystem(system::ModelKit.System) <: AbstractSystem
"""
struct ModelKitSystem{T} <: AbstractSystem
    system::ModelKit.CompiledSystem{T}
    degrees::Vector{Int}
end
#
function ModelKitSystem(system::ModelKit.System)
    csystem = ModelKit.compile(system)
    degrees = maxdegrees(system)
    ModelKitSystem(csystem, degrees)
end

Base.size(F::ModelKitSystem) = size(F.system)

cache(F::ModelKitSystem, x, p = nothing) = SystemNullCache()

@propagate_inbounds evaluate!(u, F::ModelKitSystem, x, ::SystemNullCache) =
    ModelKit.evaluate!(u, F.system, x)
evaluate(F::ModelKitSystem, x, ::SystemNullCache) = ModelKit.evaluate(F.system, x)
(F::ModelKitSystem)(x...) = evaluate(F, x...)

@propagate_inbounds jacobian!(U, F::ModelKitSystem, x, ::SystemNullCache) =
    ModelKit.jacobian!(U, F.system, x)
jacobian(F::ModelKitSystem, x, ::SystemNullCache) = ModelKit.jacobian(F.system, x)

@propagate_inbounds evaluate_and_jacobian!(u, U, F::ModelKitSystem, x, ::SystemNullCache) =
    ModelKit.evaluate_and_jacobian!(u, U, F.system, x)

degrees(F::ModelKitSystem) = F.degrees

# # parameter version
# @propagate_inbounds evaluate!(u, F::ModelKitSystem, x, p, ::SystemNullCache) =
#     ModelKit.evaluate!(u, F.system, x, p)
# evaluate(F::ModelKitSystem, x, p, ::SystemNullCache) = ModelKit.evaluate(F.system, x, p)
# @propagate_inbounds jacobian!(U, F::ModelKitSystem, x, p, ::SystemNullCache) =
#     ModelKit.jacobian!(U, F.system, x, p)
# jacobian(F::ModelKitSystem, x, p, ::SystemNullCache) = ModelKit.jacobian(F.system, x, p)
# @propagate_inbounds evaluate_jacobian!(u, U, F::ModelKitSystem, x, p, ::SystemNullCache) =
#     ModelKit.evaluate_jacobian!(u, U, F.system, x, p)
# evaluate_jacobian(F::ModelKitSystem, x, p, ::SystemNullCache) =
#     ModelKit.evaluate_jacobian(F.system, x, p)


text:
type TimematrFast{T}
    vals::DataFrame
    idx::Array{T, 1}

    function TimematrFast(vals::Array{Float64, 2}, idx::Array{T, 1})
        ## chkIdx(idx)
        ## chkNum(vals)
        if(size(vals, 1) != length(idx))
            if (length(idx) == 0) | (size(vals, 1) == 0)
                return new(DataFrame([]), Array{T, 1}[])
            end
            error(length(idx), " idx entries, but ", size(vals, 1), " rows of data")
        end
        return new(DataFrame(vals), idx)
    end
end

function TimematrFast{T}(vals::Array{Float64, 2}, idx::Array{T, 1})
    return TimematrFast{T}(vals, idx)
end


text:
module Traceur

using MacroTools
# using Vinyl: @primitive, overdub
# using ASTInterpreter2: linearize!

import Core.MethodInstance

export @trace, @trace_static

include("util.jl")
include("analysis.jl")
include("trace.jl")
# include("trace_static.jl")

end # module


text:
import Base.*

"""
    Kernel

An abstract kernel.
"""
abstract type Kernel end

"""
    CovarianceKernel

A covariance kernel, assumed symmetric in its arguments.
"""
abstract type CovarianceKernel <: Kernel end

"""
    CrossCovarianceKernel

A cross-covariance kernel, NOT necessarily symmetric in its arguments.
"""
abstract type CrossCovarianceKernel <: Kernel end

"""
    pairwise_column_difference(x::AbstractMatrix, y::AbstractMatrix)

Computes the 3-dimensional distance array, where out dimensions
are `(input_dimension, x_data_point_dimension, y_data_point_dimension)`.
"""
function pairwise_column_difference(x::AbstractMatrix, y::AbstractMatrix)
  (d1,m) = size(x)
  (d2,n) = size(y)
  reshape(x,(d1,m,1)) .- reshape(y,(d2,1,n))
end

"""
    EuclideanKernel

An abstract covariance kernel defined over a Euclidean space.
"""
abstract type EuclideanKernel <: CovarianceKernel end

include("kernels/operators.jl")
include("kernels/euclidean.jl")
include("kernels/circular.jl")

"""
    (k::Kernel)(x1::AbstractArray{<:Any,3}, x2::AbstractMatrix)

Computes a kernel matrix in batched form.
"""
function (k::Kernel)(x1::AbstractArray{<:Any,3}, x2::AbstractMatrix)
  (d,n1,n2) = size(x1)
  x1r = reshape(x1, (d,n1*n2))
  k(x1r, x2)
end


text:
using Test
using Bonsai
using Bonsai: register!, Node, match_middleware
using URIs
using HTTP

@testset "match" begin
    app = App()
    app.get("/{x}") do stream
    end
    req = Request()
    req.url = URI("http://locahost:4040/pet")
    match(app, req)
    @test_skip haskey(req.context[:params], "x")
end

@testset "register!" begin
    r = Node("*")

    # matches
    register!(r, GET, "/fish/{id}", x -> 1)
    register!(r, GET, "/fish/super", x -> 2)
    register!(r, GET, "/fish/**", x -> 3)
    # wrong method
    register!(r, PUT, "/fish/**", x -> 4)
    # wrong path
    register!(r, PUT, "/turtle/super", x -> 5)
    p = "/fish/super"
    params = Dict()
    segments = split(p, "/"; keepempty=false)
    ms = match_middleware(r, params, "GET", segments, 1)
    @test sort(map(x -> x(nothing), ms)) == [1, 2, 3]
end

@testset "app" begin

    app = App()
    app.get("**") do stream, next

    end
    app.get("**") do stream

    end

    req = HTTP.Request()
    req.method = "GET"
    req.target = "/"
    handler, middleware = match(app, req)
    @test !isnothing(handler)
    @test !isempty(middleware)
end

@testset "/" begin
    app = App()

    app.get("/") do stream
    end

    req = HTTP.Request()
    req.method = "GET"
    req.target = "/"
    handler, middleware = match(app, req)
    @test !isnothing(handler)
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
libtorch = pyimport_conda("torch", "pytorch")
py"""
import torch
import random

# test
# 0 - 4 array-dim
# strides
# reshape
# offset
# slice
# mutation
#
# combine

def random_slice(x):
  shape = x.shape
  def randrange(s):
    s1 = random.randint(0, s-1)
    s2 = random.randint(s1+1, s)
    return (s1, s2)
  slices = list(map(randrange, shape))
  if len(shape) == 1:
    s1, s2 = slices[0]
    y = x[s1:s2]
  elif len(shape) == 2:
    s11, s12 = slices[0]
    s21, s22 = slices[1]
    y = x[s11:s12, s21:s22]
  elif len(shape) == 3:
    s11, s12 = slices[0]
    s21, s22 = slices[1]
    s31, s32 = slices[2]
    y = x[s11:s12, s21:s22, s31:s32]
  return y

def random_stride(x):
  shape = x.shape
  def randrange(s):
    s1 = random.randint(0, s//2-1)
    s2 = random.randint(s//2+2, s-1)
    si = random.randint(2, 3)
    return (s1, si, s2)
  strides = list(map(randrange, shape))
  if len(shape) == 1:
    s1, si, s2 = strides[0]
    y = x[s1:s2:si]
  elif len(shape) == 2:
    s11, si1, s12 = strides[0]
    s21, si2, s22 = strides[1]
    y = x[s11:s12:si1, s21:s22:si2]
  elif len(shape) == 3:
    s11, si1, s12 = strides[0]
    s21, si2, s22 = strides[1]
    s31, si3, s32 = strides[2]
    y = x[s11:s12:si1, s21:s22:si2, s31:s32:si3]
  return y

def pyslice(x, i):
  return x[i:]

"""
thload = libtorch.load
thsave = libtorch.save
random_slice = py"random_slice"
random_stride = py"random_stride"
pyslice = py"pyslice"


text:
using Reexport

using MLStyle
using NestedTuples
import NestedTuples
using TransformVariables

function NestedTuples.schema(::Type{TransformVariables.TransformTuple{T}}) where {T}
    schema(T)
end

# In Bijectors.jl,
# logdensity_with_trans(dist, x, true) == logdensity(transformed(dist), link(dist, x))

export xform

xform(m::ConditionalModel{A, B}, _data::NamedTuple) where {A,B} = xform(m | _data)

function xform(m::ConditionalModel{A, B}) where {A,B}
    return _xform(getmoduletypencoding(m), Model(m), argvals(m), observations(m))
end

# function xform(m::Model{EmptyNTtype, B}) where {B}
#     return xform(m,NamedTuple())
# end

export sourceXform

sourceXform(m::Model) = sourceXform()(m)

function sourceXform(_data=NamedTuple())
    function(_m::Model)

        _datakeys = getntkeys(_data)
        proc(_m, st::Assign)        = :($(st.x) = $(st.rhs))
        proc(_m, st::Return)     = nothing
        proc(_m, st::LineNumber) = nothing

        function proc(_m, st::Sample)
            x = st.x
            xname = QuoteNode(x)
            rhs = st.rhs

            thecode = @q begin
                _t = xform($rhs, get(_data, $xname, NamedTuple()))
                if !isnothing(_t)
                    _result = merge(_result, ($x=_t,))
                end
            end

            # Non-leaves might be referenced later, so we need to be sure they
            # have a value
            isleaf(_m, st.x) || pushfirst!(thecode.args, :($x = testvalue($rhs)))

            return thecode
        end

        wrap(kernel) = @q begin
            _result = NamedTuple()
            $kernel
            $as(_result)
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten

    end
end

function xform(d, _data::NamedTuple)
    if hasmethod(support, (typeof(d),))
        return asTransform(support(d))
    end

    error("Not implemented:\nxform($d)")
end

using TransformVariables: ShiftedExp, ScaledShiftedLogistic

function asTransform(supp:: Dists.RealInterval)
    (lb, ub) = (supp.lb, supp.ub)

    (lb, ub) == (-Inf, Inf) && (return asℝ)
    isinf(ub) && return ShiftedExp(true,lb)
    isinf(lb) && return ShiftedExp(false,lb)
    return ScaledShiftedLogistic(ub-lb, lb)
end

# export xform
# xform(::Normal)       = asℝ
# xform(::Cauchy)       = asℝ
# xform(::Flat)         = asℝ

# xform(::HalfCauchy)   = asℝ₊
# xform(::HalfNormal)   = asℝ₊
# xform(::HalfFlat)     = asℝ₊
# xform(::InverseGamma) = asℝ₊
# xform(::Gamma)        = asℝ₊
# xform(::Exponential)  = asℝ₊

# xform(::Beta)         = as𝕀
# xform(::Uniform)      = as𝕀

xform(d, _data) = nothing

# TODO: Convert this to use `ProductMeasure`
# function xform(d::For{T,NTuple{N,Int}}, _data)  where {N,T}
#     xf1 = xform(d.f(getindex.(d.θ, 1)...), _data)
#     return as(Array, xf1, d.θ...)

#     # TODO: Implement case of unequal supports
# end

# TODO: Convert this to use `ProductMeasure`
# function xform(d::For{T,NTuple{N,UnitRange{Int}}}, _data::NamedTuple)  where {N,T}
#     xf1 = xform(d.f(getindex.(d.θ, 1)...), _data)
#     return as(Array, xf1, length.(d.θ)...)

#     # TODO: Implement case of unequal supports
# end

# TODO: Convert this to use `ProductMeasure`
# function xform(d::iid, _data::NamedTuple)
#     as(Array, xform(d.dist, _data), d.size...)
# end

# xform(d::MvNormal, _data::NamedTuple=NamedTuple()) = as(Array, size(d))

function xform(μ::AbstractMeasure,  _data::NamedTuple=NamedTuple())
    xform(representative(μ))
end

xform(μ::ProductMeasure) = as(Array, xform(first(μ.data)), size(μ.data)...)

using MeasureTheory

xform(::Lebesgue{ℝ}) = asℝ

xform(::Lebesgue{𝕀}) = as𝕀

xform(::Lebesgue{ℝ₊}) = asℝ₊

@gg function _xform(M::Type{<:TypeLevel}, _m::Model{Asub,B}, _args::A, _data) where {Asub,A,B}
    body = type2model(_m) |> sourceXform(_data) |> loadvals(_args, _data)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end