text:
module FlatArrays

# Data structures are not exported because their constructors should not be
# directly called.
export
    flatten,
    isflat

using Base:
    @propagate_inbounds,
    OneTo

import Base:
    IndexStyle,
    IndexLinear,
    axes,
    eltype,
    length,
    ndims,
    size,
    firstindex,
    lastindex,
    getindex,
    setindex!,
    iterate,
    has_offset_axes,
    convert,
    unsafe_convert

macro _inline_meta(); Expr(:meta, :inline); end

"""

All flat array types are sub-types of `AbstractFlatArray{T,N}`.

"""
abstract type AbstractFlatArray{T,N} <: DenseArray{T,N} end

"""

`FlatVector` type stores a *flat* vector.  Never directly call the
constructor but rather `flatten(A,n)` to build an instance of `FlatVector`.

"""
struct FlatVector{T,A<:AbstractArray{T}} <: AbstractFlatArray{T,1}
    parent::A
    nelem::Int
end

"""

`FlatMatrix` type stores a *flat* matrix.  Never directly call the
constructor but rather `flatten(A,m,n)` to build an instance of
`FlatMatrix`.

"""
struct FlatMatrix{T,A<:AbstractArray{T}} <: AbstractFlatArray{T,2}
    parent::A
    nelem::Int
    dims::NTuple{2,Int}
end

# Extend base methods so that flat arrays behave like arrays.

eltype(::AbstractFlatArray{T,N}) where {T,N} = T
ndims(::AbstractFlatArray{T,N}) where {T,N} = N
length(A::AbstractFlatArray) = A.nelem
has_offset_axes(::AbstractFlatArray) = false
firstindex(::AbstractFlatArray) = 1
firstindex(::AbstractFlatArray{T,N}, d::Integer) where {T,N} =
    (@_inline_meta; (d % UInt) - 1 < N ? 1 : error("dimension out of range"))
lastindex(A::AbstractFlatArray) = length(A)
lastindex(A::AbstractFlatArray, d) = size(A, d)
IndexStyle(::Type{<:AbstractFlatArray}) = IndexLinear()

unsafe_convert(::Type{Ptr{Cvoid}}, A::AbstractFlatArray) =
    unsafe_convert(Ptr{Cvoid}, A.parent)
unsafe_convert(::Type{Ptr{T}}, A::AbstractFlatArray{T}) where {T} =
    unsafe_convert(Ptr{T}, A.parent)

# The following is not needed:
# Base.IteratorSize(::Type{<:AbstractFlatArray{<:Any,N}}) where {N} =
#     Base.HasShape{N}()
# Base.IteratorEltype(::Type{<:AbstractFlatArray{<:Any,N}}) where {N} =
#     Base.HasEltype()

iterate(A::AbstractFlatArray, i=1) =
    (@_inline_meta; (i % UInt) - 1 < length(A) ?
     (@inbounds A[i], i + 1) : nothing)

size(A::FlatVector) = (A.nelem,)
size(A::FlatMatrix) = A.dims
size(A::FlatVector, d::Integer) = (d == 1 ? A.nelem :
                                   d > 1 ? 1 :
                                   error("dimension out of range"))
size(A::FlatMatrix, d::Integer) = (d > 2 ? 1 :
                                   d > 0 ? A.dims[d] :
                                   error("dimension out of range"))

axes(A::FlatVector) = (OneTo(A.nelem),)
axes(A::FlatMatrix) = map(OneTo, A.dims)
axes(A::FlatVector, d::Integer) = (d == 1 ? OneTo(A.nelem) :
                                   d > 1 ? OneTo(1) :
                                   error("dimension out of range"))
axes(A::FlatMatrix, d::Integer) = (d > 2 ? OneTo(1) :
                                   d > 0 ? OneTo(A.dims[d]) :
                                   error("dimension out of range"))

strides(A::FlatVector) = (1,)
strides(A::FlatMatrix) = (1, size(A, 1))
stride(A::FlatVector, d::Integer) = (d == 1 ? 1 :
                                     d >= 2 ? length(A) :
                                     error("dimension out of range"))
stride(A::FlatMatrix, d::Integer) = (d == 1 ? 1 :
                                     d == 2 ? size(A, 1) :
                                     d >= 3 ? length(A) :
                                     error("dimension out of range"))

# Indexation.

@inline @propagate_inbounds getindex(A::AbstractFlatArray, i) =
    getindex(A, _index1(A, i)...)

@inline @propagate_inbounds getindex(A::AbstractFlatArray, i, j) =
    getindex(A, _index1(A, i, j)...)

@inline function getindex(A::AbstractFlatArray, i::Int)
    @boundscheck checkbounds(A.parent, i)
    @inbounds val = A.parent[i]
    return val
end

@inline function getindex(A::AbstractFlatArray, i::Int, j::Int)
    @boundscheck checkbounds(A.parent, i, j)
    @inbounds val = A.parent[_index2(A, i, j)]
    return val
end

@inline @propagate_inbounds setindex!(A::AbstractFlatArray, val, i) =
    setindex!(A, val, _index1(A, i)...)

@inline @propagate_inbounds setindex!(A::AbstractFlatArray, val, i, j) =
    setindex!(A, val, _index1(A, i, j)...)

@inline function setindex!(A::AbstractFlatArray, val, i::Int)
    @boundscheck checkbounds(A.parent, i)
    @inbounds A.parent[i] = val
    return val
end

@inline function setindex!(A::AbstractFlatArray, val, i::Int, j::Int)
    @boundscheck checkbounds(A.parent, i, j)
    @inbounds A.parent[_index2(A, i, j)] = val
    return val
end

# First pass on indices to make them integers (must return a tuple).
@inline _index1(A::AbstractFlatArray, i::Integer) = (Int(i),)
@inline _index1(A::AbstractFlatArray, i::Integer, j::Integer) = (Int(i), Int(j))
@inline _index1(A::AbstractFlatArray{<:Any,1}, I::CartesianIndex{1}) = (I[1],)
@inline _index1(A::AbstractFlatArray{<:Any,2}, I::CartesianIndex{2}) = (I[1], I[2])

# Second pass on indices to make them a single index.
@inline _index2(A::FlatVector, i::Int, j::Int) = i + A.nelem*(j - 1)
@inline _index2(A::FlatMatrix, i::Int, j::Int) = i + A.dims[1]*(j - 1)

"""

```julia
flatten(A, n) -> V
```

yields a *flat* vector `V` of length `n` whose elements are given by array `A`.

```julia
flatten(A, m, n) -> M
```

yields a *flat* matrix `M` with `m` rows and `n` columns whose elements are
given by array `A`.

A *flat* array has 1-based linear indexing, column-major storage order and all
its elements are contiguous.  This may be very convenient for calling efficient
code which expects such kind of arrays.

If the argument `A` of `flatten` is already a *flat* array, the returned array
will share the contents of `A`; otherwise an independent copy of the elements
is made.

If the contents of `A` can be shared, calling `flatten` is faster (by a factor
~ 7 for a vector, ~ 5 for a matrix) and uses less memory than calling
[`reshape`](@ref).  This may be significant for very small arrays.  Also most
base methods like `firstindex`, `lastindex`, etc. applied to an
`AbstractFlatArray` can be optimized out at compilation time.  Finally, the
array returned by `flatten` is guaranteed to have 1-based linear indexing and
column-major storage order.  However, contrary to `resize`, the array returned
by `flatten` may not share its contents with `A`.  In fact, if `isflat(A)` is
true, then `flatten(A,...)` will return an array sharing its contents with `A`;
otherwise, the contents is not shared.

!!! warning
    Many fast operations on flat arrays rely on the fact that the parent array
    backing the storage of the flat array is never resized (*e.g.*, with the
    [`resize!`](@ref) method).

See also: [`isflat`](@ref), [`resize`](@ref).

""" flatten

flatten(arr::AbstractArray, dims::NTuple{1,Integer}) =
    flatten(arr, dims[1])
flatten(arr::AbstractArray, nelem::Integer) =
    flatten(arr, Int(nelem))
@inline flatten(arr::A, nelem::Int) where {T,A<:AbstractArray{T}} =
    FlatVector(_flatten(arr, nelem), nelem)

flatten(arr::AbstractArray, dims::NTuple{2,Integer}) =
    flatten(arr, dims[1], dims[2])
flatten(arr::AbstractArray, nrows::Integer, ncols::Integer) =
    flatten(arr, Int(nrows), Int(ncols))
@inline function flatten(arr::A,
                         nrows::Int, ncols::Int) where {T,A<:AbstractArray{T}}
    nrows > 0 && ncols > 0 || error("invalid dimensions")
    nelem = nrows*ncols
    return FlatMatrix(_flatten(arr, nelem), nelem, (nrows, ncols))
end

@inline function _flatten(A::Array, nelem::Int)
    (nelem % UInt) - 1 < length(A) || error("incompatible length")
    return A
end

@inline function _flatten(A::AbstractFlatArray, nelem::Int)
    (nelem % UInt) - 1 < length(A) || error("incompatible length")
    return A.parent
end

@inline function _flatten(A::AbstractArray, nelem::Int)
    (nelem % UInt) - 1 < length(A) || error("incompatible length")
    return (isflat(A) ? A : _flatcopy(A, nelem))
end

function _flatcopy(A::AbstractArray{T}, nelem::Int) :: Vector{T} where {T}
    F = Array{T,1}(undef, nelem)
    F[1], state = iterate(A)
    @inbounds @simd for j in 2:nelem
        F[j], state = iterate(A, state)
    end
    return F
end

"""

`isflat(A)` yields whether array `A` is a *flat* array, that is an array
with 1-based indices and whose elements are all continuous.

See also: [`flatten`](@ref).

"""
isflat(::AbstractFlatArray) = true
isflat(::Array) = true
isflat(A::T) where {T<:AbstractArray} = _isflat(IndexStyle(T), A)

# The test below could also be `has_offset_axes(A) == false`
_isflat(::IndexLinear, A::AbstractArray) = (firstindex(A) == 1)

function _isflat(::IndexCartesian, A::AbstractArray{T,N}) where {T,N}
    #
    # We could have checked the strides with the following code:
    #
    #     n = 1
    #     @inbounds for d in 1:N
    #         stride(A, d) == n || return false
    #         n *= size(A, d)
    #     end
    #     return true
    #
    # However, (i) not all array types implement the stride method, (ii)
    # non-linear indexing style indicates that the array is certainly
    # non-flat (otherwise indexation is not optimized for that kind of array).
    # Therefore we simply return false.
    return false
end

end # module

max_stars_repo_name: emmt/FlatArrays

max_issues_repo_name: emmt/FlatArrays

max_forks_repo_name: emmt/FlatArrays


text:
using StructJuMP, JuMP
using StructJuMPSolverInterface

include("select_solver.jl")

#############
# A sample model
#############
scen = 2
m = StructuredModel(num_scenarios=scen)
@variable(m, x[1:2])
@variable(m, z[1:2])
@NLobjective(m, Min, x[1]^2+x[2]^2+x[1]*z[2] + x[1]*x[2])
@NLconstraint(m, x[1] + x[2]<= 100)

for i in 1:scen
    bl = StructuredModel(parent=m)
    @variable(bl, y[1:2])
    @NLconstraint(bl, -100<= x[2]*x[1] + x[1]*y[1]*y[2] <= 100)
    @NLconstraint(bl, -100<= x[2]*x[1]*y[1]*y[2] <= 100)
    @NLobjective(bl, Min, (y[1]+y[2])^2)
end

structJuMPSolve(m)

getVarValue(m)

#verification
# x[1]^2+x[2]^2+x[1]*x[4] + x[1]*x[2] + (y1[1]+y1[2])^2 + (y2[1]+y2[2])^2

# x[1] + x[2]
# x[2]*x[1] + x[1]*y1[1]*y1[2]
# x[2]*x[1]*y1[1]*y1[2]
# x[2]*x[1] + x[1]*y2[1]*y2[2]
# x[2]*x[1]*y2[1]*y2[2]

ext: jl


text:
# test of deleteHybridEdge!, functions to extract displayed trees/subnetworks,
#      used to compare networks with the hardwired cluster distance.
# Cecile March 2016

if !isdefined(:doalltests) doalltests = false; end

#----------------------------------------------------------#
#   testing functions to delete edges and nodes            #
#----------------------------------------------------------#
@testset "testing deleteHybridEdge!" begin

if doalltests
# example of network with one hybrid edge connected to the root:
net = readTopology("((Adif:1.0,(Aech:0.122,#H6:10.0::0.047):10.0):1.614,Aten:1.0,((Asub:1.0,Agem:1.0):0.0)#H6:5.062::0.953);");
# plot(net, showEdgeNumber=true, showNodeNumber=true)
deleteHybridEdge!(net, net.edge[10]);
writeTopologyLevel1(net) == "(Adif:1.0,(Aech:0.122,(Asub:1.0,Agem:1.0):10.0):10.0,Aten:2.614);" ||
 error("deleteHybridEdge! didn't work on 10th edge")
net = readTopology("((Adif:1.0,(Aech:0.122,#H6:10.0::0.047):10.0):1.614,Aten:1.0,((Asub:1.0,Agem:1.0):0.0)#H6:5.062::0.953);");
deleteHybridEdge!(net, net.edge[3]);
writeTopologyLevel1(net) == "((Adif:1.0,Aech:10.122):1.614,Aten:1.0,(Asub:1.0,Agem:1.0):5.062);" ||
 error("deleteHybridEdge! didn't work on 3rd edge")
# plot(net, showEdgeNumber=true, showNodeNumber=true)

net=readTopology("(4,((1,(2)#H7:::0.864):2.069,(6,5):3.423):0.265,(3,#H7:::0.1361111):10.0);");
# plot(net, showEdgeNumber=true, showNodeNumber=true)
deleteHybridEdge!(net, net.edge[11]);
writeTopologyLevel1(net) == "(4,((1,2):2.069,(6,5):3.423):0.265,3);" ||
 error("deleteHybridEdge! didn't work on 11th edge")
net=readTopology("(4,((1,(2)#H7:::0.864):2.069,(6,5):3.423):0.265,(3,#H7:::0.1361111):10.0);");
deleteHybridEdge!(net, net.edge[4]);
writeTopologyLevel1(net) == "(4,((6,5):3.423,1):0.265,(3,2):10.0);" ||
 error("deleteHybridEdge! didn't work on 4th edge")
end

# example with wrong attributed inChild1
net=readTopology("(4,((1,(2)#H7:::0.864):2.069,(6,5):3.423):0.265,(3,#H7:::0.1361111):10.0);");
net.edge[5].isChild1 = false;
@test_nowarn deleteHybridEdge!(net, net.edge[4]);
@test writeTopologyLevel1(net) == "(4,((6,5):3.423,1):0.265,(3,2):10.0);"
# or: deleteHybridEdge! didn't work on 4th edge when isChild1 was outdated

if doalltests
net = readTopology("((Adif:1.0,(Aech:0.122,#H6:10.0::0.047):10.0):1.614,Aten:1.0,((Asub:1.0,Agem:1.0):0.0)#H6:5.062::0.953);");
net.edge[5].isChild1 = false # edge 5 from -1 to -2
deleteHybridEdge!(net, net.edge[10]);
println("a warning is expected: \"node -1 being the root is contradicted by isChild1 of its edges.\"")
writeTopologyLevel1(net) == "(Adif:1.0,(Aech:0.122,(Asub:1.0,Agem:1.0):10.0):10.0,Aten:2.614);" ||
 error("deleteHybridEdge! didn't work on 10th edge after isChild1 was changed")
# plot(net, showEdgeNumber=true, showNodeNumber=true)
end

end # of testing deleteHybridEdge!

@testset "testing deleteleaf! and hardwiredClusterDistance" begin

cui2str = "(Xgordoni,Xmeyeri,(Xcouchianus,(Xvariatus,(Xevelynae,((Xxiphidium,#H25:9.992::0.167):1.383,(Xmilleri,(Xandersi,(Xmaculatus,((((Xhellerii,(Xalvarezi,Xmayae):0.327):0.259,Xsignum):1.866,(Xclemenciae_F2,Xmonticolus):1.461):0.786,((((Xmontezumae,(Xnezahuacoyotl)#H26:0.247::0.807):0.372,((Xbirchmanni_GARC,Xmalinche_CHIC2):1.003,Xcortezi):0.454):0.63,((Xcontinens,Xpygmaeus):1.927,((Xnigrensis,Xmultilineatus):1.304,#H26:0.0::0.193):0.059):2.492):2.034)#H25:0.707::0.833):1.029):0.654):0.469):0.295):0.41):0.646):3.509):0.263);"
cui3str = "(Xmayae,((Xhellerii,(((Xclemenciae_F2,Xmonticolus):1.458,(((((Xmontezumae,(Xnezahuacoyotl)#H26:0.247::0.804):0.375,((Xbirchmanni_GARC,Xmalinche_CHIC2):0.997,Xcortezi):0.455):0.63,(#H26:0.0::0.196,((Xcontinens,Xpygmaeus):1.932,(Xnigrensis,Xmultilineatus):1.401):0.042):2.439):2.0)#H7:0.787::0.835,(Xmaculatus,(Xandersi,(Xmilleri,((Xxiphidium,#H7:9.563::0.165):1.409,(Xevelynae,(Xvariatus,(Xcouchianus,(Xgordoni,Xmeyeri):0.263):3.532):0.642):0.411):0.295):0.468):0.654):1.022):0.788):1.917)#H27:0.149::0.572):0.668,Xalvarezi):0.257,(Xsignum,#H27:1.381::0.428):4.669);"

net3  = readTopology(cui3str);
net2  = readTopology(cui2str);

if doalltests
# major tree, root with outgroup then delete 2 leaves:
tree2 = majorTree(net2);        tree3 = majorTree(net3);
rootatnode!(tree2,"Xmayae");    rootatnode!(tree3,"Xmayae");
deleteleaf!(tree2,"Xhellerii"); deleteleaf!(tree3,"Xhellerii");
deleteleaf!(tree2,"Xsignum");   deleteleaf!(tree3,"Xsignum");
hardwiredClusterDistance(tree2, tree3, false) == 0 ||
  error("HWDist not 0, major tree, root then prune");
# major tree, delete 2 leaves then root with outgroup:
tree2 = majorTree(net2);        tree3 = majorTree(net3);
deleteleaf!(tree2,"Xhellerii"); deleteleaf!(tree3,"Xhellerii");
deleteleaf!(tree2,"Xsignum");   deleteleaf!(tree3,"Xsignum");
hardwiredClusterDistance(tree2, tree3, false) == 0 || error("HWD not 0, major tree - 2 taxa");
hardwiredClusterDistance(tree2, tree3, true) == 20 || error("rooted RF dist not 20");
rootatnode!(tree3,"Xmaculatus");
hardwiredClusterDistance(tree2, tree3, true) == 15 || error("rooted RF dist not 15");
rootatnode!(tree2,"Xgordoni");
hardwiredClusterDistance(tree2, tree3, true) == 16 || error("rooted RF dist not 16");
hardwiredClusterDistance(tree2, tree3, false) == 0 || error("HWD not 0, major tree - 2 taxa");
end

# network: delete 2 leaves
net2  = readTopology(cui2str);
@test_nowarn deleteleaf!(net2,"Xhellerii");
@test_nowarn deleteleaf!(net2,"Xsignum");
@test_warn """node 13 is a leaf. Will create a new node if needed, to set taxon "Xmayae" as outgroup.""" rootatnode!(net2,"Xmayae");
net3  = readTopology(cui3str);
@test_nowarn deleteleaf!(net3,"Xhellerii");
@test_nowarn deleteleaf!(net3,"Xsignum");
@test_warn """node 13 is a leaf. Will create a new node if needed, to set taxon "Xmayae" as outgroup.""" rootatnode!(net2,"Xmayae");
@test hardwiredClusterDistance(net2, net3, true) == 3
@test_nowarn deleteleaf!(net3,"Xmayae");    #plot(net3);
@test net3.numHybrids == 2
# using simplify=false in deleteleaf!
net3  = readTopology(cui3str);
deleteleaf!(net3,"Xhellerii"); deleteleaf!(net3,"Xsignum");
deleteleaf!(net3,"Xmayae", simplify=false);
@test net3.numHybrids==3 # or: deleteleaf wrong on mayae with simplify=false
# plot(net3); # looks weird though: k=2 cycle at the root. 3 root edges:
# one to a leaf, 1 major & 1 minor hybrid edge to the same child.

end # of testset for deleteleaf! and hardwiredClusterDistance

#----------------------------------------------------------#
#   testing functions to display trees / subnetworks       #
#----------------------------------------------------------#

@testset "testing deleteHybridThreshold!" begin

if doalltests
net21 = readTopology("(A,((B,#H1),(C,(D)#H1)));");
# manual bug fix to get good gamma's (0.5 not 1.0) and major/minor
net21.edge[3].gamma = 0.5;
net21.edge[7].gamma = 0.5;
deleteHybridThreshold!(net21,0.2);
writeTopologyLevel1(net21) == "(A,((B,#H1:::0.5),(C,(D)#H1:::0.5)));" ||
 error("deleteHybridThreshold! didn't work on net21, gamma=0.2")
deleteHybridThreshold!(net21,0.5);
writeTopologyLevel1(net21) == "(A,((C,D),B));" ||
 error("deleteHybridThreshold! didn't work on net21, gamma=0.5")

net22 = readTopology("(A,((B,#H1:::0.2),(C,(D)#H1:::0.8)));");
deleteHybridThreshold!(net22,0.3);
writeTopologyLevel1(net22) == "(A,((C,D),B));" ||
 error("deleteHybridThreshold! didn't work on net22, gamma=0.3")

net31 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(C:0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
net32 = deepcopy(net31); for e in net32.edge e.length=-1.0; end
# plot(net31)
deleteHybridThreshold!(net31,0.3);
writeTopologyLevel1(net31) == "(A:1.0,((C:0.9,D:1.1):1.3,B:2.3):0.7);" ||
 error("deleteHybridThreshold! didn't work on net31, gamma=0.3")
deleteHybridThreshold!(net32,0.3)
writeTopologyLevel1(net32) == "(A,((C,D),B));" ||
 error("deleteHybridThreshold! didn't work on net32, gamma=0.3")

net42 = readTopology("(A:1.0,((B:1.1,#H1:0.2):1.2,(C:0.9,(D:0.8)#H1:0.3):1.3):0.7):0.1;");
deleteHybridThreshold!(net42,0.5);
writeTopologyLevel1(net42) == "(A:1.0,((C:0.9,D:1.1):1.3,B:2.3):0.7);" ||
 error("deleteHybridThreshold! didn't work on net42, gamma=0.5")
end

net5 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
# plot(net5)
@test_nowarn deleteHybridThreshold!(net5,0.5);  # both H1 and H2 eliminated
@test writeTopologyLevel1(net5) == "(A:1.0,((((C:0.52,E:0.52):0.6,F:1.5):0.9,D:1.1):1.3,B:2.3):0.7);"
# or: deleteHybridThreshold! didn't work on net5, gamma=0.5
net5 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
@test_nowarn deleteHybridThreshold!(net5,0.3);  # H2 remains
@test writeTopologyLevel1(net5) == "(A:1.0,((((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,D:1.1):1.3,B:2.3):0.7);"
# or: deleteHybridThreshold! didn't work on net5, gamma=0.3

end # of testset, deleteHybridThreshold!

@testset "testing displayedNetworks! and displayedTrees" begin

net3 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(C:0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
net31 = displayedNetworks!(net3, net3.node[6]); #H1 = 6th node
@test writeTopologyLevel1(net31) == "(A:1.0,((B:1.1,D:1.0):1.2,C:2.2):0.7);"
# or: displayedNetworks! didn't work on net3, minor at 6th node
@test writeTopologyLevel1(net3)  == "(A:1.0,((C:0.9,D:1.1):1.3,B:2.3):0.7);"
# or: displayedNetworks! didn't work on net3, major at 6th node

if doalltests
net3 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(C:0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
a = displayedTrees(net3, 0.2);
length(a) == 2 ||
 error("displayedTrees didn't work on net3, gamma=0.2: output not of length 2")
writeTopologyLevel1(a[1]) == "(A:1.0,((C:0.9,D:1.1):1.3,B:2.3):0.7);" ||
 error("displayedTrees didn't work on net3, gamma=0.2: 1st tree wrong")
writeTopologyLevel1(a[2]) == "(A:1.0,((B:1.1,D:1.0):1.2,C:2.2):0.7);" ||
 error("displayedTrees didn't work on net3, gamma=0.2: 2nd tree wrong")
end

net5 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
a = displayedTrees(net5, 0.5);
@test length(a) == 1 # or: displayedTrees didn't work on net5, gamma=0.5
@test writeTopologyLevel1(a[1]) == "(A:1.0,((((C:0.52,E:0.52):0.6,F:1.5):0.9,D:1.1):1.3,B:2.3):0.7);"
# or: displayedTrees didn't work on net5, gamma=0.5
a = displayedTrees(net5, 0.1);
@test length(a) == 4 # or: displayedTrees didn't work on net5, gamma=0.1
@test writeTopologyLevel1(a[1]) == "(A:1.0,((((C:0.52,E:0.52):0.6,F:1.5):0.9,D:1.1):1.3,B:2.3):0.7);"
@test writeTopologyLevel1(a[2]) == "(A:1.0,((B:1.1,D:1.0):1.2,((C:0.52,E:0.52):0.6,F:1.5):2.2):0.7);"
@test writeTopologyLevel1(a[3]) == "(A:1.0,((((F:0.7,E:0.51):0.8,C:1.12):0.9,D:1.1):1.3,B:2.3):0.7);"
@test writeTopologyLevel1(a[4]) == "(A:1.0,((B:1.1,D:1.0):1.2,((F:0.7,E:0.51):0.8,C:1.12):2.2):0.7);"

end # of testset, displayedNetworks! & displayedTrees

@testset "testing majorTree and displayedNetworkAt!" begin

net5 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
@test writeTopology(majorTree(net5)) == "(A:1.0,((((C:0.52,E:0.52):0.6,F:1.5):0.9,D:1.1):1.3,B:2.3):0.7);"
@test_nowarn displayedNetworkAt!(net5, net5.hybrid[1]);
@test writeTopology(net5) == "(A:1.0,((((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,D:1.1):1.3,B:2.3):0.7);"

end # of testset, majorTree & displayedNetworkAt!

#----------------------------------------------------------#
#   testing functions to compare trees                     #
#----------------------------------------------------------#

@testset "testing tree2Matrix" begin

if doalltests
net5 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
tree = displayedTrees(net5, 0.0);
taxa = tipLabels(net5);
M1 = tree2Matrix(tree[1], taxa, rooted=false);
M2 = tree2Matrix(tree[2], taxa, rooted=false);
M1 ==
[15 0 0 1 1 1 1;
 12 0 0 1 1 1 0;
  8 0 0 1 1 0 0] || error("bad M1, from net5")
M2 ==
[ 4 0 1 0 0 0 1;
 12 0 0 1 1 1 0;
  8 0 0 1 1 0 0] || error("bad M2, from net5")
hardwiredClusterDistance(tree[1], tree[2], true) == 2 || error("wrong dist between rooted trees 1 and 2");
hardwiredClusterDistance(tree[1], tree[2], false)== 2 || error("wrong dist between unrooted trees 1 and 2");
hardwiredClusters(net5, taxa) ==
[16 0 1 1 1 1 1 10;
  4 0 1 0 0 0 1 10;
  3 0 0 0 0 0 1 11;
 15 0 0 1 1 1 1 10;
 12 0 0 1 1 1 0 10;
  8 0 0 1 1 0 0 10;
  7 0 0 0 1 0 0 11;
 11 0 0 0 1 1 0 10] || error("wrong matrix of hardwired clusters for net5");
hardwiredClusterDistance(net5, tree[2], true) == 4 || error("wrong dist between net5 and tree 2");
end

## check things with R
# library(ape); set.seed(15); phy <- rmtree(10,8); write.tree(phy,"tenRtrees.tre")
# library(phangorn); RF.dist(phy,rooted=F)
#     1  2  3  4  5  6  7  8  9
# 2   8
# 3  10 10
# 4  10 10 10
# 5  10  8  8 10
# 6  10 10  8 10  8
# 7  10 10 10 10 10 10
# 8  10  8  8 10 10  8 10
# 9  10 10  8  8 10 10 10 10
# 10  8  6 10  8  8 10 10 10  8
# i=10; j=1; RF.dist(phy[[i]],phy[[j]],rooted=T) # 10
# i=10; j=2; RF.dist(phy[[i]],phy[[j]],rooted=T) #  8
## now checking if we get the same results in julia
# phy = readInputTrees("tenRtrees.tre")
# or
phy1 = readTopology("((t8:0.8623136566,(((t6:0.141187073,t2:0.7767125128):0.9646669542,t4:0.8037273993):0.447443719,t5:0.7933459524):0.8417851452):0.7066285675,(t1:0.0580010619,(t7:0.6590069213,t3:0.1069735419):0.5657461432):0.3575631182);");
phy2 = readTopology("((t3:0.9152618761,t4:0.4574306419):0.7603277895,(((t1:0.4291725352,t8:0.3302786439):0.3437780738,(t5:0.8438980761,(t6:0.6667000714,t2:0.7141199473):0.01087239943):0.752832541):0.2591188031,t7:0.7685037958):0.9210739341);");
phy3 = readTopology("(((t7:0.3309174306,t6:0.8330178803):0.7741786113,(((t2:0.4048132468,t8:0.6809111023):0.6810255498,(t4:0.6540613638,t5:0.2610215396):0.8490990005):0.6802781771,t3:0.2325445588):0.911911567):0.94644987,t1:0.09404937108);");
phy10= readTopology("((t4:0.1083955287,((t1:0.8376079942,t8:0.1745392387):0.6178579947,((t6:0.3196466176,t2:0.9228881211):0.3112748025,t7:0.05162345758):0.7137957355):0.5162231021):0.06693460606,(t5:0.005652675638,t3:0.2584615161):0.7333540542);");

@test hardwiredClusterDistance(phy1, phy10, false)== 8
@test hardwiredClusterDistance(phy2, phy10, false)== 6
@test hardwiredClusterDistance(phy3, phy10, false)==10
@test hardwiredClusterDistance(phy1, phy2,  false)== 8
@test hardwiredClusterDistance(phy1, phy3,  false)==10
@test hardwiredClusterDistance(phy2, phy3,  false)==10
@test hardwiredClusterDistance(phy1, phy10, true) ==10
@test hardwiredClusterDistance(phy2, phy10, true) == 8
# or: wrong RF distance between some of the trees

end # of testset, tree2Matrix

#----------------------------------------------------------#
#   testing function to compare networks                   #
#   with hardwired clusters                                #
#   used for detection of given hybridization event        #
#----------------------------------------------------------#

@testset "test displayedTrees, hardwiredClusters, hardwiredClusterDistance, displayedNetworkAt!" begin

estnet = readTopology("(6,((5,#H7:0.0::0.402):8.735,((1,2):6.107,((3,4):1.069)#H7:9.509::0.598):6.029):0.752);")
# originally from "../msSNaQ/simulations/estimatedNetworks/baseline/nloci10/1_julia.out"
trunet = readTopology("((((1,2),((3,4))#H1),(#H1,5)),6);");
@test hardwiredClusterDistance(majorTree(trunet), majorTree(estnet),false) == 0 # false: unrooted
truminor = minorTreeAt(trunet, 1); # (1:1.0,2:1.0,((5:1.0,(3:1.0,4:1.0):2.0):1.0,6:1.0):2.0);
estminor = minorTreeAt(estnet, 1); # (5:1.0,(3:1.0,4:1.0):1.069,(6:1.0,(1:1.0,2:1.0):10.0):8.735);
@test writeTopology(truminor) == "((((1,2),(3,4)),5),6);"
@test writeTopology(estminor) == "(6,(((1,2):6.107,(3,4):10.578):6.029,5):0.752);"
@test hardwiredClusterDistance(truminor, estminor, false) == 0 # false: unrooted
# so the hybrid edge was estimated correctly!!

net5 = readTopology("(A,((B,#H1:::0.2),(((C,(E)#H2:::0.7),(#H2:::0.3,F)),(D)#H1:::0.8)));");
tree = displayedTrees(net5, 0.0);
taxa = tipLabels(net5);
@test hardwiredClusters(tree[1], taxa) ==
[16 0 1 1 1 1 1 10;
 15 0 0 1 1 1 1 10;
 12 0 0 1 1 1 0 10;
  8 0 0 1 1 0 0 10]
@test hardwiredClusters(tree[2], taxa) ==
[16 0 1 1 1 1 1 10;
  4 0 1 0 0 0 1 10;
 12 0 0 1 1 1 0 10;
  8 0 0 1 1 0 0 10]
@test hardwiredClusters(net5, taxa) ==
[16 0 1 1 1 1 1 10;
  4 0 1 0 0 0 1 10;
  3 0 0 0 0 0 1 11;
 15 0 0 1 1 1 1 10;
 12 0 0 1 1 1 0 10;
  8 0 0 1 1 0 0 10;
  7 0 0 0 1 0 0 11;
 11 0 0 0 1 1 0 10]

if doalltests
trunet = readTopology("(((1,2),((3,4))#H1),(#H1,5),6);"); # unrooted
taxa = tipLabels(trunet);
hardwiredClusters(trunet, taxa) ==
[8 1 1 1 1 0 0 10
 3 1 1 0 0 0 0 10
 7 0 0 1 1 0 0 11
 6 0 0 1 1 0 0 10
11 0 0 1 1 1 0 10] || error("wrong hardwired cluster matrix for unrooted trunet");
trunet = readTopology("((((1,2),((3,4))#H1),(#H1,5)),6);"); # rooted: good!
hardwiredClusters(trunet, taxa) ==
[12 1 1 1 1 1 0 10;
  8 1 1 1 1 0 0 10;
  3 1 1 0 0 0 0 10;
  7 0 0 1 1 0 0 11;
  6 0 0 1 1 0 0 10;
 11 0 0 1 1 1 0 10] || error("wrong hardwired cluster matrix for trunet");
hardwiredClusters(estnet, taxa) ==
[13 1 1 1 1 1 0 10
  4 0 0 1 1 1 0 10
  3 0 0 1 1 0 0 11
 10 0 0 1 1 0 0 10
 12 1 1 1 1 0 0 10
  7 1 1 0 0 0 0 10] || error("wrong hardwired cluster matrix for estnet")
hardwiredClusterDistance(trunet,estnet,true) == 0 ||
 error("trunet and estnet should be found to be at HWDist 0");

net51 = readTopologyLevel1("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;")
directEdges!(net51); # doing this avoids a warning on the next line:
displayedNetworkAt!(net51, net51.hybrid[1]) # H2
# "WARNING: node -3 being the root is contradicted by isChild1 of its edges."
rootatnode!(net51, "A");
writeTopologyLevel1(net51) == "(A:0.5,((((C:0.52,(E:0.5)#H2:0.02):0.6,(#H2:0.01,F:0.7):0.8):0.9,D:1.1):1.3,B:2.3):0.5);" ||
  error("wrong net51 after displayedNetworkAt!");
net52 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7):0.1;");
displayedNetworkAt!(net52, net52.hybrid[2])
writeTopologyLevel1(net52) == "(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,E:0.52):0.6,F:1.5):0.9,(D:0.8)#H1:0.3::0.8):1.3):0.7);" ||
  error("wrong net52 after displayedNetworkAt!");
taxa = tipLabels(net52); # order: A B C E F D
hardwiredClusters(net51, taxa) ==
[16 0 1 1 1 1 1 10;
 15 0 0 1 1 1 1 10;
 12 0 0 1 1 1 0 10;
  8 0 0 1 1 0 0 10;
  7 0 0 0 1 0 0 11;
 11 0 0 0 1 1 0 10] || error("wrong hardwired clusters for net51");
hardwiredClusters(net52, taxa) ==
[16 0 1 1 1 1 1 10;
  4 0 1 0 0 0 1 10;
  3 0 0 0 0 0 1 11;
 15 0 0 1 1 1 1 10;
 12 0 0 1 1 1 0 10;
  8 0 0 1 1 0 0 10] || error("wrong hardwired clusters for net52");
hardwiredClusterDistance(net51,net52,true) == 4 ||
 error("wrong HWDist between net51 and net52");
end

end # of testset: displayedTrees, hardwiredClusters, hardwiredClusterDistance, displayedNetworkAt!

@testset "testing hardwiredCluster! on single nodes" begin

net5 = "(A,((B,#H1),(((C,(E)#H2),(#H2,F)),(D)#H1)));" |> readTopology |> directEdges! ;
taxa = net5 |> tipLabels # ABC EF D
m = hcat([true,false,false,false,false,false],
[false,true,false,false,false,false],
[false,false,false,false,false,true],
[false,true,false,false,false,true],
[false,false,true,false,false,false],
[false,false,false,true,false,false],
[false,false,false,true,false,false],
[false,false,true,true,false,false],
[false,false,false,true,false,false],
[false,false,false,false,true,false],
[false,false,false,true,true,false],
[false,false,true,true,true,false],
[false,false,false,false,false,true],
[false,false,false,false,false,true],
[false,false,true,true,true,true],
[false,true,true,true,true,true])
for i = 1:16
  @test hardwiredCluster(net5.edge[i], taxa) == m[:,i]
end
end # of testset, hardwiredCluster! on single nodes


text:
let
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4], 1, 4) == [2,3,4,1]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4], 2, 4) == [3,4,1,2]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4,5,6,7], 1, 4) == [2,3,4,1,5,6,7]
    @test AutomotiveDrivingModels.cyclic_shift_left!([1,2,3,4,5,6,7], 2, 4) == [3,4,1,2,5,6,7]

    td = get_test_trajdata()
    scene = Scene()
#    print(`scene=`)
#    println(scene)

    g = get!(scene, td, 1)
#    print(`g=`)
#    println(g)
    col = get_first_collision(g)
#    print(`col=`)
#    println(col)
    @test col.A == 1
    @test col.B == 2
    @test col.is_colliding == true

    @test get_first_collision(get!(scene, td, 2), CPAMemory()).is_colliding == true
end


text:
import DelimitedFiles

data = DelimitedFiles.readdlm("data/day_one.txt", ',', Int)
complements = 2020 .- data

println(prod(intersect(complements, data)))

for j in data
    complements = 2020 - j .- data
    matches = intersect(complements, data)
    if length(matches) > 0
        println(j * prod(matches))
        break
    end
end

max_stars_repo_name: ColCarroll/adventofcode

max_issues_repo_name: ColCarroll/adventofcode

max_forks_repo_name: ColCarroll/adventofcode


text:
"""
    μ ± σ

Creates $DEFAUL_NUM_PARTICLES `Particles` with mean `μ` and std `σ`.
If `μ` is a vector, the constructor `MvNormal` is used, and `σ` is thus treated as std if it's a scalar, and variances if it's a matrix or vector.
"""
±(μ::Real,σ) = μ + σ*Particles(DEFAUL_NUM_PARTICLES)
±(μ::AbstractVector,σ) = Particles(DEFAUL_NUM_PARTICLES, MvNormal(μ, σ))
∓(μ::Real,σ) = μ + σ*StaticParticles(DEFAUL_STATIC_NUM_PARTICLES)
∓(μ::AbstractVector,σ) = StaticParticles(DEFAUL_STATIC_NUM_PARTICLES, MvNormal(μ, σ))

(..)(a,b) = Particles(DEFAUL_NUM_PARTICLES, Uniform(a,b))

"""
    ⊗(μ,σ) = outer_product(Normal.(μ,σ))

See also `outer_product`
"""
⊗(μ,σ) = outer_product(Normal.(μ,σ))

"""
    p = outer_product(dists::Vector{<:Distribution}, N=100_000)

Creates a multivariate systematic sample where each dimension is sampled according to the corresponding univariate distribution in `dists`. Returns `p::Vector{Particles}` where each Particles has a length approximately equal to `N`.
The particles form the outer product between `d` systematically sampled vectors with length given by the d:th root of N, where `d` is the length of `dists`, All particles will be independent and have marginal distributions given by `dists`.

See also `MonteCarloMeasurements.⊗`
"""
function outer_product(dists::AbstractVector{<:Distribution}, N=100_000)
    d = length(dists)
    N = floor(Int,N^(1/d))
    dims = map(dists) do dist
        v = systematic_sample(N,dist; permute=true)
    end
    cart_prod = vec(collect(Iterators.product(dims...)))
    p = map(1:d) do i
        Particles(getindex.(cart_prod,i))
    end
end

# StaticParticles(N::Integer = DEFAUL_NUM_PARTICLES; permute=true) = StaticParticles{Float64,N}(SVector{N,Float64}(systematic_sample(N, permute=permute)))

function print_functions_to_extend()
    excluded_functions = [fill, |>, <, display, show, promote, promote_rule, promote_type, size, length, ndims, convert, isapprox, ≈, <, (<=), (==), zeros, zero, eltype, getproperty, fieldtype, rand, randn]
    functions_to_extend = setdiff(names(Base), Symbol.(excluded_functions))
    for fs in functions_to_extend
        ff = @eval $fs
        ff isa Function || continue
        isempty(methods(ff)) && continue # Sort out intrinsics and builtins
        f = nameof(ff)
        if !isempty(methods(ff, (Real,Real)))
            println(f, ",")
        end
    end
end
"""
    shortform(p::AbstractParticles)
Return a short string describing the type
"""
shortform(p::Particles) = "Part"
shortform(p::StaticParticles) = "SPart"
shortform(p::WeightedParticles) = "WPart"
function to_num_str(p, d=3)
    s = std(p)
    if s < eps(p)
        string(round(mean(p), sigdigits=d))
    else
        string(round(mean(p), sigdigits=d), " ± ", round(s, sigdigits=d-1))
    end
end
function Base.show(io::IO, ::MIME"text/plain", p::AbstractParticles{T,N}) where {T,N}
    sPT = shortform(p)
    print(io, "$(sPT)$N(", to_num_str(p, 4),")")
end

function Base.show(io::IO, p::AbstractParticles{T,N}) where {T,N}
    print(io, to_num_str(p, 3))
end
# function Base.show(io::IO, p::MvParticles)
#     sPT = shortform(p)
#     print(io, "(", N, " $sPT with mean ", round.(mean(p), sigdigits=3), " and std ", round.(sqrt.(diag(cov(p))), sigdigits=3),")")
# end
for mime in (MIME"text/x-tex", MIME"text/x-latex")
    @eval function Base.show(io::IO, ::$mime, p::AbstractParticles)
        print(io, "\$"); show(io, p); print("\$")
    end
end

# Two-argument functions
foreach(register_primitive_multi, [+,-,*,/,//,^,
max,min,mod,mod1,atan,atand,add_sum,hypot])
# One-argument functions
foreach(register_primitive_single, [*,+,-,/,
exp,exp2,exp10,expm1,
log,log10,log2,log1p,
sin,cos,tan,sind,cosd,tand,sinh,cosh,tanh,
asin,acos,atan,asind,acosd,atand,asinh,acosh,atanh,
zero,sign,abs,sqrt,rad2deg,deg2rad])

for PT in (:Particles, :StaticParticles)
    # Constructors
    @eval begin
        $PT(v::Vector) = $PT{eltype(v),length(v)}(v)
        function $PT{T,N}(n::Real) where {T,N} # This constructor is potentially dangerous, replace with convert?
            v = fill(n,N)
            $PT{T,N}(v)
        end

        """
            ℝⁿ2ℝⁿ_function(f::Function, p::AbstractArray{T})
        Applies  `f : ℝⁿ → ℝⁿ` to an array of particles.
        """
        function ℝⁿ2ℝⁿ_function(f::F, p::AbstractArray{$PT{T,N}}) where {F,T,N}
            individuals = map(1:length(p[1])) do i
                f(getindex.(p,i))
            end
            RT = promote_type(T,eltype(eltype(individuals)))
            PRT = $PT{RT,N}
            out = similar(p, PRT)
            for i = 1:length(p)
                out[i] = PRT(getindex.(individuals,i))
            end
            reshape(out, size(p))
        end

        function ℝⁿ2ℝⁿ_function(f::F, p::AbstractArray{$PT{T,N}}, p2::AbstractArray{$PT{T,N}}) where {F,T,N}
            individuals = map(1:length(p[1])) do i
                f(getindex.(p,i), getindex.(p2,i))
            end
            RT = promote_type(T,eltype(eltype(individuals)))
            PRT = $PT{RT,N}
            out = similar(p, PRT)
            for i = 1:length(p)
                out[i] = PRT(getindex.(individuals,i))
            end
            reshape(out, size(p))
        end
    end
    for ff in (var, std)
        f = nameof(ff)
        @eval function (Statistics.$f)(p::$PT{T,N},args...;kwargs...) where {T,N}
            N == 1 && (return zero(T))
            $f(p.particles, args...;kwargs...)
        end
    end
    @forward @eval($PT).particles Statistics.mean, Statistics.cov, Statistics.median, Statistics.quantile, Statistics.middle
end

for PT in (:WeightedParticles,)
    # Constructors
    @eval begin
        $PT(v::Vector, w=fill(-log(length(v)), length(v))) = $PT{eltype(v),length(v)}(v,w)
        function $PT{T,N}(n::Real) where {T,N} # This constructor is potentially dangerous, replace with convert?
            v = fill(n,N)
            w = fill(-log(N),N)
            $PT{T,N}(v,w)
        end
    end
    # Two-argument functions
    for ff in (+,-,*,/,//,^, max,min,mod,mod1,atan,add_sum)
        f = nameof(ff)
        @eval begin
            function (Base.$f)(p::$PT{T,N},a::Real...) where {T,N}
                # $PT{T,N}(map(x->$f(x,a...), p.particles),p.logweights)
                $PT{T,N}($f.(p.particles, maybe_particles.(a)...), .+(p.logweights, maybe_logweights.(a)...))
            end
            function (Base.$f)(a::Real,p::$PT{T,N}) where {T,N}
                $PT{T,N}(map(x->$f(a,x), p.particles),p.logweights)
            end
            function (Base.$f)(p1::$PT{T,N},p2::$PT{T,N}) where {T,N} # TODO: add logweights for multiplication but not for addition
                sf = string($f)
                $f == Base.:(*) || @warn("p1 $sf p2 not yet fully supported for WeightedParticles")
                $PT{T,N}(map($f, p1.particles, p2.particles),p1.logweights+p2.logweights)
            end
            function (Base.$f)(p1::$PT{T,N},p2::$PT{S,N}) where {T,S,N} # Needed for particles of different float types :/
                $PT{promote_type(T,S),N}(map($f, p1.particles, p2.particles),p1.logweights+p2.logweights)
            end
        end
    end
    # One-argument functions
    for ff in [*,+,-,/,
        exp,exp2,exp10,expm1,
        log,log10,log2,log1p,
        sin,cos,tan,sind,cosd,tand,sinh,cosh,tanh,
        asin,acos,atan,asind,acosd,atand,asinh,acosh,atanh,
        zero,sign,abs,sqrt,rad2deg,deg2rad]
        f = nameof(ff)
        @eval function (Base.$f)(p::$PT)
            $PT(map($f, p.particles), p.logweights)
        end
    end
    for ff in [mean, median, quantile]
        f = nameof(ff)
        @eval function (Statistics.$f)(p::$PT,args...;kwargs...)
            $f(p.particles, Weights(exp.(p.logweights)), args...;kwargs...)
        end
    end
    for ff in [var, std]
        f = nameof(ff)
        @eval function (Statistics.$f)(p::$PT{T,N},args...;kwargs...) where {T,N}
            N == 1 && (return zero(T))
            $f(p.particles, AnalyticWeights(exp.(p.logweights)), args...;kwargs...)
        end
    end
    @eval Statistics.cov(p::$PT,args...;kwargs...) = var(p,args...;kwargs...)
end

function Particles(d::Distribution;kwargs...)
    Particles(DEFAUL_NUM_PARTICLES, d; kwargs...)
end

function StaticParticles(d::Distribution;kwargs...)
    StaticParticles(DEFAUL_STATIC_NUM_PARTICLES, d; kwargs...)
end

for PT in (:Particles, :StaticParticles, :WeightedParticles)
    @forward @eval($PT).particles Base.iterate, Base.extrema, Base.minimum, Base.maximum

    @eval begin
        $PT{T,N}(p::$PT{T,N}) where {T,N} = p

        function $PT(m::Matrix)
            map(1:size(m,2)) do i
                $PT{eltype(m),size(m,1)}(@view(m[:,i]))
            end
        end

        function $PT(N::Integer=DEFAUL_NUM_PARTICLES, d::Distribution=Normal(0,1); permute=true, systematic=true)
            if systematic
                v = systematic_sample(N,d; permute=permute)
            else
                v = rand(d, N)
            end
            $PT{eltype(v),N}(v)
        end

        function $PT(N::Integer, d::MultivariateDistribution)
            v = rand(d,N)' |> copy # For cache locality
            $PT(v)
        end

        nakedtypeof(p::$PT{T,N}) where {T,N} = $PT
        nakedtypeof(::Type{$PT{T,N}}) where {T,N} = $PT
    end
    # @eval begin

    # end
    @eval begin
        Base.length(::Type{$PT{T,N}}) where {T,N} = N
        Base.eltype(::Type{$PT{T,N}}) where {T,N} = $PT{T,N}
        Base.promote_rule(::Type{S}, ::Type{$PT{T,N}}) where {S,T,N} = $PT{promote_type(S,T),N} # This is hard to hit due to method for real 3 lines down
        Base.promote_rule(::Type{Bool}, ::Type{$PT{T,N}}) where {T,N} = $PT{promote_type(Bool,T),N} # Needed since above is not specific enough.

        Base.convert(::Type{StaticParticles{T,N}}, p::$PT{T,N}) where {T,N} = StaticParticles(p.particles)
        Base.convert(::Type{$PT{T,N}}, f::Real) where {T,N} = $PT{T,N}(fill(T(f),N))
        Base.convert(::Type{$PT{T,N}}, f::$PT{S,N}) where {T,N,S} = $PT{promote_type(T,S),N}(promote_type(T,S).(f.particles))
        function Base.convert(::Type{S}, p::$PT{T,N}) where {S<:ConcreteFloat,T,N}
            N == 1 && (return S(p[1]))
            std(p) < eps(S) || throw(ArgumentError("Cannot convert a particle distribution to a float if not all particles are the same."))
            return S(p[1])
        end
        function Base.convert(::Type{S}, p::$PT{T,N}) where {S<:ConcreteInt,T,N}
            isinteger(p) || throw(ArgumentError("Cannot convert a particle distribution to an int if not all particles are the same."))
            return S(p[1])
        end
        Base.zeros(::Type{$PT{T,N}}, dim::Integer) where {T,N} = [$PT(zeros(eltype(T),N)) for d = 1:dim]
        Base.zero(::Type{$PT{T,N}}) where {T,N} = $PT(zeros(eltype(T),N))
        Base.isfinite(p::$PT{T,N}) where {T,N} = isfinite(mean(p))
        Base.round(p::$PT{T,N}, r::RoundingMode, args...; kwargs...) where {T,N} = round(mean(p), r, args...; kwargs...)
        Base.round(::Type{S}, p::$PT{T,N}, args...; kwargs...) where {S,T,N} = round(S, mean(p), args...; kwargs...)
        function Base.AbstractFloat(p::$PT{T,N}) where {T,N}
            N == 1 && (return p[1])
            std(p) < eps(T) || throw(ArgumentError("Cannot convert a particle distribution to a number if not all particles are the same."))
            return p[1]
        end

        """
        union(p1::AbstractParticles, p2::AbstractParticles)

        A `Particles` containing all particles from both `p1` and `p2`. Note, this will be twice as long as `p1` or `p2` and thus of a different type.
        `pu = Particles([p1.particles; p2.particles])`
        """
        function Base.union(p1::$PT{T,NT},p2::$PT{T,NS}) where {T,NT,NS}
            $PT([p1.particles; p2.particles])
        end

        """
        intersect(p1::AbstractParticles, p2::AbstractParticles)

        A `Particles` containing all particles from the common support of `p1` and `p2`. Note, this will be of undetermined length and thus undetermined type.
        """
        function Base.intersect(p1::$PT,p2::$PT)
            mi = max(minimum(p1),minimum(p2))
            ma = min(maximum(p1),maximum(p2))
            f = x-> mi <= x <= ma
            $PT([filter(f, p1.particles); filter(f, p2.particles)])
        end

        Base.:^(p::$PT, i::Integer) = $PT(p.particles.^i) # Resolves ambiguity
        Base.:\(p::Vector{<:$PT}, p2::Vector{<:$PT}) = Matrix(p)\Matrix(p2) # Must be here to be most specific
    end

end

Base.length(p::AbstractParticles{T,N}) where {T,N} = N
Base.ndims(p::AbstractParticles{T,N}) where {T,N} = ndims(T)
Base.:\(H::MvParticles,p::AbstractParticles) = Matrix(H)\p.particles
# Base.:\(p::AbstractParticles, H) = p.particles\H
# Base.:\(p::MvParticles, H) = Matrix(p)\H
# Base.:\(H,p::MvParticles) = H\Matrix(p)

Base.Broadcast.broadcastable(p::AbstractParticles) = Ref(p)
Base.setindex!(p::AbstractParticles, val, i::Integer) = setindex!(p.particles, val, i)
Base.getindex(p::AbstractParticles, i::Integer) = getindex(p.particles, i)
# Base.getindex(v::MvParticles, i::Int, j::Int) = v[j][i] # Defining this methods screws with show(::MvParticles)

Base.Matrix(v::MvParticles) = reduce(hcat, getfield.(v,:particles))
# function Statistics.var(v::MvParticles,args...;kwargs...) # Not sure if it's a good idea to define this. Is needed for when var(v::AbstractArray) is used
#     s2 = map(1:length(v[1])) do i
#         var(getindex.(v,i))
#     end
#     eltype(v)(s2)
# end

Statistics.mean(v::MvParticles) = mean.(v)
Statistics.cov(v::MvParticles,args...;kwargs...) = cov(Matrix(v), args...; kwargs...)
Distributions.fit(d::Type{<:MultivariateDistribution}, p::MvParticles) = fit(d,Matrix(p)')
Distributions.fit(d::Type{<:Distribution}, p::AbstractParticles) = fit(d,p.particles)
check_similar_weights(v) = all(v[i].logweights ≈ v[1].logweights for i in 2:length(v))
function Statistics.cov(v::MvWParticles,args...;kwargs...)
    check_similar_weights(v)
    cov(Matrix(v), AnalyticWeights(v[1].logweights); kwargs...)
end
Distributions.fit(d::Type{<:MultivariateDistribution}, p::MvWParticles) = error("Not implemented for weighted particles yet")
Distributions.fit(d::Type{<:Distribution}, p::WeightedParticles) = error("Not implemented for weighted particles yet")
Distributions.fit(d::Type{<:MvNormal}, p::MvWParticles) = MvNormal(p)
Distributions.fit(d::Type{<:Normal}, p::WeightedParticles) = Normal(p)
Distributions.Normal(p::AbstractParticles) = Normal(mean(p), std(p))
Distributions.MvNormal(p::AbstractParticles) = MvNormal(mean(p), cov(p))
Distributions.MvNormal(p::MvParticles) = MvNormal(mean(p), cov(p))

meanstd(p::AbstractParticles) = std(p)/sqrt(length(p))
meanvar(p::AbstractParticles) = var(p)/length(p)

Base.:(==)(p1::AbstractParticles{T,N},p2::AbstractParticles{T,N}) where {T,N} = p1.particles == p2.particles
Base.:(!=)(p1::AbstractParticles{T,N},p2::AbstractParticles{T,N}) where {T,N} = p1.particles != p2.particles

function _comparison_operator(p)
    length(p) == 1 && return
    USE_UNSAFE_COMPARIONS[] || error("Comparison operators are not well defined for uncertain values and are currently turned off. Call `unsafe_comparisons(true)` to enable comparison operators for particles using the current reduction function $(COMPARISON_FUNCTION[]). Change this function using `set_comparison_function(f)`.")
end

function Base.:<(a::Real,p::AbstractParticles)
    _comparison_operator(p)
    a < COMPARISON_FUNCTION[](p)
end
function Base.:<(p::AbstractParticles,a::Real)
    _comparison_operator(p)
    COMPARISON_FUNCTION[](p) < a
end
function Base.:<(p::AbstractParticles, a::AbstractParticles)
    _comparison_operator(p)
    COMPARISON_FUNCTION[](p) < COMPARISON_FUNCTION[](a)
end
function Base.:(<=)(p::AbstractParticles{T,N}, a::AbstractParticles{T,N}) where {T,N}
    _comparison_operator(p)
    COMPARISON_FUNCTION[](p) <= COMPARISON_FUNCTION[](a)
end

"""
    p1 ≈ p2

Determine if two particles are not significantly different
"""
Base.:≈(p::AbstractParticles, a::AbstractParticles, lim=2) = abs(mean(p)-mean(a))/(2sqrt(std(p)^2 + std(a)^2)) < lim
Base.:≈(a::Real,p::AbstractParticles, lim=2) = abs(mean(p)-a)/std(p) < lim
Base.:≈(p::AbstractParticles, a::Real, lim=2) = abs(mean(p)-a)/std(p) < lim
Base.:≈(p::MvParticles, a::AbstractVector) = all(a ≈ b for (a,b) in zip(a,p))
Base.:≈(a::AbstractVector, p::MvParticles) = all(a ≈ b for (a,b) in zip(a,p))
Base.:≈(a::MvParticles, p::MvParticles) = all(a ≈ b for (a,b) in zip(a,p))
Base.:≉(a,b::AbstractParticles,lim=2) = !(≈(a,b,lim))
Base.:≉(a::AbstractParticles,b,lim=2) = !(≈(a,b,lim))
"""
    p1 ≉ p2

Determine if two particles are significantly different
"""
Base.:≉(a::AbstractParticles,b::AbstractParticles,lim=2) = !(≈(a,b,lim))

Base.sincos(x::AbstractParticles) = sin(x),cos(x)
Base.minmax(x::AbstractParticles,y::AbstractParticles) = (min(x,y), max(x,y))

Base.:!(p::AbstractParticles) = all(p.particles .== 0)

Base.isinteger(p::AbstractParticles) = all(isinteger, p.particles)
Base.iszero(p::AbstractParticles) = all(iszero, p.particles)
Base.iszero(p::AbstractParticles, tol) = abs(mean(p.particles)) < tol

≲(a::Real,p::AbstractParticles,lim=2) = (mean(p)-a)/std(p) > lim
≲(p::AbstractParticles,a::Real,lim=2) = (a-mean(p))/std(p) > lim
≲(p::AbstractParticles,a::AbstractParticles,lim=2) = (mean(p)-mean(a))/(2sqrt(std(p)^2 + std(a)^2)) > lim
≳(a::Real,p::AbstractParticles,lim=2) = ≲(p,a,lim)
≳(p::AbstractParticles,a::Real,lim=2) = ≲(a,p,lim)
≳(p::AbstractParticles,a::AbstractParticles,lim=2) = ≲(a,p,lim)
Base.eps(p::Type{<:AbstractParticles{T,N}}) where {T,N} = eps(T)
Base.eps(p::AbstractParticles{T,N}) where {T,N} = eps(T)
Base.eps(p::AbstractParticles{<:Complex{T},N}) where {T,N} = eps(T)

"""
    norm(x::AbstractParticles, p=2)

if p == 2: return abs(mean(x))
elseif p == Inf: return max(extrema(x)...)
"""
function LinearAlgebra.norm(x::AbstractParticles, p::Union{AbstractFloat, Integer}=2)
    if p == 2
        return abs(mean(x))
    elseif p == Inf
        return max(extrema(x)...)
    end
    throw(ArgumentError("Cannot take $(p)-norm of particles"))
end

Base.exp(p::AbstractMatrix{<:AbstractParticles}) = ℝⁿ2ℝⁿ_function(exp, p)
LinearAlgebra.lyap(p1::Matrix{<:AbstractParticles}, p2::Matrix{<:AbstractParticles}) = ℝⁿ2ℝⁿ_function(lyap, p1, p2)

# OBS: defining this was a very bad idea, eigvals jump around and get confused with each other etc.
# function LinearAlgebra.eigvals(p::Matrix{$PT{T,N}}) where {T,N} # Special case to propte types differently
#     individuals = map(1:length(p[1])) do i
#         eigvals(getindex.(p,i))
#     end
#
#     PRT = Complex{$PT{T,N}}
#     out = Vector{PRT}(undef, length(individuals[1]))
#     for i = eachindex(out)
#         c = getindex.(individuals,i)
#         out[i] = complex($PT{T,N}(real(c)),$PT{T,N}(imag(c)))
#     end
#     out
# end


text:
# by Lazaro Alonso
# original post
# https://nextjournal.com/lazarus/extracting-dominant-colours-from-pictures
using CairoMakie, Clustering, Colors, Random
using Images: load, RGB, channelview, permutedims
using Downloads: download
CairoMakie.activate!()
#let
    """
    getImg(; url = "", file = "")
    """
    function getImg(; url = "", file = "")
        try
            out = file != "" ? RGB.(load(file)) : RGB.(load(download(url)))
            return out
        catch e
            out = [RGB{Float64}(rand(3)...) for i in 1:600, j in 1:400]
            return out
        end
    end

    function getColorPalette(image; ncolors = 5)
        img_CHWa = channelview(image)
        #img_CHW = permutedims(img_CHWa, (1,3,2))
        testmat = reshape(img_CHWa, (3, size(image)[1]*size(image)[2])) #input shape
        sol = kmeans(testmat, ncolors)
        csize = counts(sol) # get the cluster sizes
        colores = sol.centers # get the cluster centers, dominat colors
        indxc = sortperm(csize)
        colores, indxc
    end
    #https://commons.wikimedia.org/wiki/Category:Ceyx_erithaca?uselang=de
    #url = "https://upload.wikimedia.org/wikipedia/commons/b/b8/Oriental_dwarf_kingfisher.jpg"
    url = "https://upload.wikimedia.org/wikipedia/commons/thumb/1/12/Oriental_dwarf_kingfisher_a_Jewel_of_India%27s_Rainforests_07.jpg/1280px-Oriental_dwarf_kingfisher_a_Jewel_of_India%27s_Rainforests_07.jpg?uselang=de"
    #path = "/path/file/file.jpg"
    img = getImg(; url = url) # ;url= url, file = file
    #img = img[150:end-40,250:end-40]
    #img = img[550:end-560,1200:end-1400]
    img = img[150:end-100,300:end-350]
    # colores, indxc = getColorPalette(img; ncolors = ncolors)
    # color = [RGB(colores[:,indxc[i]]...) for i in 1:ncolors]

    function plotPalette(;img=img, marker = :diamond, ncolors = 8)
        Random.seed!(123)
        yshift = 0.04*size(img)[1]
        colores, indxc = getColorPalette(img; ncolors = ncolors)
        color = [RGB(colores[:,indxc[i]]...) for i in 1:ncolors]
        hexNames = [hex(color[i]) for i in 1:ncolors]

        xcolor = fill(1.1*size(img)[2], ncolors)
        ycolor = LinRange(yshift,size(img)[1]-yshift,ncolors)
        fig, ax, _= image(rotr90(img), figure = (resolution = (1.45size(img)[2], size(img)[1]),),
            axis = (aspect = DataAspect(),))
        scatterlines!(ax, xcolor, ycolor, color = color, markercolor = color, marker = marker, #'■',
            markersize = 0.1*size(img)[1], linewidth = 2)
        for i in 1:ncolors
            text!(ax, "\"#"*hexNames[i]*"\"", position = (1.07xcolor[i], ycolor[i]-5),
                textsize = 0.025*size(img)[1])
        end
        xlims!(nothing, 1.45*size(img)[2])
        hidedecorations!(ax)
        hidespines!(ax)
        fig
    end
    fig = plotPalette()
    # fig = with_theme(plotPalette, theme_black())
    save(joinpath(@__DIR__, "output", "getColorPalette.png"), fig, px_per_unit = 2.0) # HIDE
#end
# output image under License
# CC BY-SA 4.0
# https://creativecommons.org/licenses/by-sa/4.0/
# due to:
# https://commons.wikimedia.org/wiki/Category:Ceyx_erithaca?uselang=de#/media/File:Oriental_dwarf_kingfisher_a_Jewel_of_India's_Rainforests_07.jpg

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Images

"""
    expand dimensions, zeros, randn
"""
expand_dim(dim...) = (dim..., 3, 1)
zeros_like(T::Type, dims...) = fill!(similar(T, dims...), 0f0)
zeros_like(xs::AbstractArray, dims...) = fill!(similar(xs, dims...), 0f0)
randn_like(T::Type, dims...) = randn!(similar(T, dims...))
randn_like(xs::AbstractArray, dims...) = randn!(similar(xs, dims...))

"""
    size
"""
function size_pyramid(scale, min_size, image_size)
    current_size = min_size
    pyramid = Vector{Tuple{Int64,Int64}}()
    for i in 1:100
        push!(pyramid, current_size)
        current_size == image_size && break
        current_size = @. floor(Int64, min_size * scale^i)
        current_size = min.(current_size, image_size)
    end
    return pyramid
end

channel_pyramid(n_stage) = min.(map(s->32 * 2^(floor(Int64, s / 4)), 1:n_stage), 128)

"""
    resize
"""
function resize_and_padding(x::Array{Float32,4},
            image_shape::Tuple{Int64,Int64}, padded_shape::Tuple{Int64,Int64})
    # println(size(x), image_shape, padded_shape)
    x_large = imresize(view(x, :, :, :, 1), image_shape...)
    xx = zeros(Float32, expand_dim(padded_shape...))
    pad1, pad2 = (@. (padded_shape - image_shape) ÷ 2)[1:2]
    xx[1 + pad1:image_shape[1] + pad1, 1 + pad2:image_shape[2] + pad2, : , 1] = x_large
    return xx
end

function resize_and_padding(x::CuArray{Float32,4},
            image_shape::Tuple{Int64,Int64}, padded_shape::Tuple{Int64,Int64})
    return cu(resize_and_padding(adapt(Array{Float32}, x), image_shape, padded_shape))
end

"""
    image pyramid
"""
function build_image_pyramid(img::AbstractArray, image_shapes::Vector{Tuple{Int64,Int64}}, noise_shapes::Vector{Tuple{Int64,Int64}})
    return map((is, ns)->resize_and_padding(img, is, ns), image_shapes, noise_shapes)
end

function build_zero_pyramid(xs::AbstractArray, shapes::Vector{Tuple{Int64,Int64}})
    return map(s->zeros_like(xs, expand_dim(s...)), shapes)
end

"""
    noise pyramid
"""
function build_noise_pyramid(xs::AbstractArray, shapes::Vector{Tuple{Int64,Int64}}, amplifiers::Vector{Float32})
    return map((s, a)->a * randn_like(xs, expand_dim(s...)), shapes, amplifiers)
end

function build_rec_pyramid(xs::AbstractArray, shapes::Vector{Tuple{Int64,Int64}}, amplifier::Float32)
    v = build_zero_pyramid(xs, shapes)
    randn!(v[1])
    v[1] *= amplifier
    return v
end

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
# validate PROWIM cl-alpha sweep, props off:
alphas = range(-5, stop=12, length = 18) .* pi/180
wing_b = PROWIMData.wing["span"]
wing_TR = PROWIMData.wing["TR"] # hershey bar
wing_c = PROWIMData.wing["chord"]
wing_AR = wing_b / wing_c
wing_θroot = PROWIMData.wing["theta_root"]
wing_θtip = PROWIMData.wing["theta_tip"]
wing_le_sweep = PROWIMData.wing["le_sweep"]
wing_ϕ = PROWIMData.wing["dihedral"]
plot_base_name = "PROWIM_props_off"
data_PROWIM_CL_alpha = AS.cl_alpha_sweep_template(alphas, wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip, wing_le_sweep, wing_ϕ;
        plot_directory = joinpath(AS.topdirectory, "data","plots",AS.TODAY),
        plot_base_name = plot_base_name,
        plot_extension = ".pdf"
    )
objective = AS.runsimulation!(data_PROWIM_CL_alpha...)

# Fetch data.
CL_data = PROWIMData.props_off["CL_alpha"]
CD_data_balance = PROWIMData.props_off["CD_alpha_balance"]
CD_data_wake_survey = PROWIMData.props_off["CD_alpha_wake_survey"]

# Update plots.
fig_props_off_cl_alpha = plt.figure(plot_base_name * "_cl_alpha_sweep")
axs = fig_props_off_cl_alpha.get_axes()
axs[3].scatter(CL_data[:,1], CL_data[:,2], marker = "x", label="Veldhuis")
axs[3].legend(loc="upper left", bbox_to_anchor=(1.01,1))
axs[1].scatter(CD_data_balance[:,1], CD_data_balance[:,2], marker = "+", label = "balance")
axs[1].scatter(CD_data_wake_survey[:,1], CD_data_wake_survey[:,2], marker = "x", label = "wake survey")
axs[1].legend(loc="upper left", bbox_to_anchor=(1.01,1))

fig_props_off_cl_alpha.tight_layout()
if save_figs
    fig_props_off_cl_alpha.savefig(joinpath(plot_directory,"PROWIM_props_off_cl_alpha_sweep.pdf"), bbox_inches="tight")
end

# @test objective == 0

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using IOLA
using IOLA.Utils
using MDCT
using DSP
using WAV
using Gnuplot
using Statistics

codec = "wma"

params = Codec.getparams(Codec.WMA)
seg_hop = params.hop
transform = Codec.gettransform(params)

bitrates = [128, 196, 256, 320]

diffs = Vector{Array{Float64, 2}}()

(y, fs, nbits, opt) = wavread("data/audio.wav")
seg_len = 4*params.length

push!(diffs, analyze(y[1:Int(10fs),1], transform, seg_len, seg_hop))

for bitr in bitrates
    (y, fs, nbits, opt) = wavread("data/$codec/audio-$bitr.wav")
    push!(diffs, analyze(y[1:Int(10fs),1], transform, seg_len, seg_hop))
end

xs = 3:seg_hop
periodf = zeros(length(xs), 5)
for i = 1:length(xs), j = 1:5
    per = mod.(diffs[j][:,2],xs[i])
    periodf[i, j] = std(per; corrected=false)/mean(per)
end

@gp(:GP2,
    xs, periodf[:,2], "w l tit '128k'",
    xs, periodf[:,3], "w l tit '196k'",
    xs, periodf[:,4], "w l tit '256k'",
    xs, periodf[:,5], "w l tit '320k'",
    xs, periodf[:,1], "w l tit 'none'"
   )

period = zeros(5)
phis = zeros(length(diffs[1][:,2]), 5)
for i = 1:5
    period[i] = findperiod(diffs[i][:,2], xs)
    phis[:, i] = 2π*mod.(diffs[i][:,2], period[i])./period[i]
end

@gp(:GP1, "set angles radians",
    "set polar",
    "set grid polar 15. lt -1 dt 0 lw 0.5",
    "unset border", "unset param",
    "unset xtics",
    phis[:,2], diffs[2][:,1], "pt 7 ps 2 t '128k'",
    phis[:,3], diffs[3][:,1], "pt 7 ps 2 t '196k'",
    phis[:,4], diffs[4][:,1], "pt 7 ps 2 t '256k'",
    phis[:,5], diffs[5][:,1], "pt 7 ps 2 t '320k'",
    phis[:,1], diffs[1][:,1], "pt 7 ps 2 t 'none'"
   )

roms = zeros(5)
for i = 1:5
    roms[i] = radiusofmean(diffs[i][:,1], phis[:,i])
end

max_stars_repo_name: Ashymad/IOLA

max_issues_repo_name: Ashymad/IOLA

max_forks_repo_name: Ashymad/IOLA


text:
# old goldstandard
const goldstandard_file = "C:\\Users\\Daniel\\Documents\\Bitbucket\\musical-grammars\\data\\annotations\\Goldstandard.txt"
const chroma = Dict("C"=>0,"D"=>2,"E"=>4,"F"=>5,"G"=>7,"A"=>9,"B"=>11)

function pitch_class(note_name::AbstractString)
    m = match(r"([A-G])(b*|#*)", note_name)
    if contains(m[2], "#")
        chroma[m[1]] + length(m[2])
    else
        chroma[m[1]] - length(m[2])
    end
end

function translate_goldstandard_notation1(chord::AbstractString)
    m = match(r"([A-G]b*#*)(.*)", chord)
    string(pitch_class(m[1]), m[2])
end

function translate_goldstandard_notation1(sequence::AbstractVector)
    map(translate_goldstandard_notation1, sequence)
end

const chord_form_dict = Dict(""=>"maj", "-"=>"min", "%"=>"hdim", "o"=>"dim")

function translate_goldstandard_notation2(chord::AbstractString)
    # println(chord)
    m = match(r"([A-G]b*#*)(.*)", chord)
    # println(m)
    string(pitch_class(m[1]), chord_form_dict[m[2]])
end

function translate_goldstandard_notation2(sequence::AbstractVector)
    map(translate_goldstandard_notation2, sequence)
end

function produce_goldstandard_sequences()
    lines = map(strip, readlines(goldstandard_file))
    for line in lines
        if !isempty(line) && line[1] == '\$' && line[end] == '\$'
            if contains(line[2:end-1], "maj")
                produce(translate_goldstandard_notation1(split(line[2:end-1])))
            else
                produce(translate_goldstandard_notation2(split(line[2:end-1])))
            end
        end
    end
end

function read_goldstandard()
    collect(@task(produce_goldstandard_sequences()))
end

#new treebank
const treebank_chord_form_dict = Dict(""=>"maj", "7"=>"maj", "^7"=>"maj", "6"=>"maj",
                                      "m"=>"min", "m7"=>"min", "m6"=>"min",
                                      "%7"=>"hdim", "o"=>"dim", "o7"=>"dim", "+"=>"aug")

function translate_treebank_notation(chord::AbstractString)
    m = match(r"([A-G]b*#*)(.*)", chord)
    string(pitch_class(m[1]), treebank_chord_form_dict[m[2]])
end

function translate_treebank_notation(chords::AbstractVector)
    map(translate_treebank_notation, chords)
end

function right_branching_tree(num_leafs::Int)
    if num_leafs == 1
        TreeNode("x")
    else
        n = TreeNode("x")
        insert_child!(n, right_branching_tree(num_leafs - 1))
        insert_child!(n, TreeNode("x"))
        n
    end
end

function read_tree(treebank_file::AbstractString)
    chomp(readlines(treebank_file)[4])[3:end]
end

function read_treebank_trees()
    dir = joinpath(dirname(@__FILE__), "..", "..", "trees")
    files = ["afternoon_in_paris.tsv", "all_of_me.tsv", "alone_together.tsv", "anthropology.tsv", "a_night_in_tunisia.tsv"]
    # files = ["afternoon_in_paris.tsv", "all_of_me.tsv"]
    trees = TreeNode{String}[]
    for file in files
        t = tree(read_tree(joinpath(dir, file)))
        for n in t
            n.data = translate_treebank_notation(n.data)
        end
        push!(trees, t)
    end
    trees
end

function read_goldstandard_trees_1()
    trees = map(tree,
    ["[Cmaj[Cmaj[Cmaj[Cmaj[Cmaj][Cmaj[Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Dmin][Gmaj]]][Cmaj]]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Dmin][Gmaj]]][Cmaj]]]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Dmin][Gmaj]]][Cmaj[Cmaj][Cmaj]]]]][Cmaj[Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Gmaj[Dmin[Amaj][Dmin]][Gmaj]][Cmaj]]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Dmin][Gmaj]]][Cmaj]]]]"
    "[Cmaj[Cmaj[Cmaj][Gmaj[Dmin[Amaj[Emaj][Amaj]][Dmin]][Gmaj[Dmaj[Amin[Emaj[Bmin][Emaj]][Amin][Dmaj]][Gmaj[Dmin][Gmaj]]]]]][Cmaj[Cmaj[Cmaj][Cmaj[Gmaj[Dmin[Amaj[Emaj][Amaj]][Dmin]][Gmaj[Fmaj][Gmaj[Fmin][Gmaj[Dmin[Amaj[Emin][Amaj]][Dmin]][Gmaj]]]]][Cmaj]]][Gmaj[Dmin[Ebdim][Dmin]][Gmaj]]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Gmaj[Gmaj[Dmin[Amin[Emaj[Bhdim[Fmaj][Bhdim]][Emaj]][Amin]][Dmin[Ebmaj][Dmin]]][Gmaj]][Gmaj[Dmin[Amin[Emin][Amin]][Dmin]][Gmaj]]]][Cmaj[Cmaj[Gmaj[Dmaj[Amaj[Emin][Amaj]][Dmaj[Amin][Dmaj]]][Gmaj[Dmin][Gmaj]]][Cmaj]][Cmaj[Gmaj[Gmaj[Dmin[Amin[Emaj[Bhdim[Fmaj][Bhdim]][Emaj]][Amin]][Dmin[Ebmaj][Dmin]]][Gmaj]][Gmaj[Dmin[Amin[Emin][Amin]][Dmin]][Gmaj]]][Cmaj]]]][Cmaj[Gmaj[Amin[Gmaj[Dhdim[Amin][Dhdim]][Gmaj]][Amin[Emin][Amin]]][Gmaj[Fmaj[Cmaj[Gmaj[Dhdim][Gmaj]][Cmaj]][Fmaj]][Gmaj[Dmin[Dmin[Amaj[Emin[Bmaj[F#min][Bmaj]][Emin]][Amaj]][Dmin]][Dmin[Dmin[Amaj][Dmin]][Dmin[Amaj][Dmin]]]][Gmaj[Abmaj][Gmaj]]]]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Gmaj[Gmaj[Dmin[Amin[Emaj[Bhdim[Fmaj][Bhdim]][Emaj]][Amin]][Dmin[Ebmaj][Dmin]]][Gmaj]][Gmaj[Dmin[Amin[Emin][Amin]][Dmin]][Gmaj]]][Cmaj]]]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj[Cmaj[Cmaj][Cmaj[Fmaj][Cmaj]]][Gmaj[Dmaj[Amin][Dmaj]][Gmaj[Gmaj][Gmaj]]]][Cmaj[Cmaj[Cmaj][Cmaj[Fmaj][Cmaj]]][Cmaj[Gmaj][Cmaj]]]][Cmaj[Cmaj[Cmaj[Cmaj][Cmaj[Fmaj][Cmaj]]][Gmaj[Dmaj[Amin][Dmaj]][Gmaj[Gmaj][Gmaj]]]][Cmaj[Cmaj[Cmaj][Cmaj[Fmaj][Cmaj]]][Gmaj[Gmaj][Gmaj[Abmaj][Gmaj]]]]]][Cmaj[Cmaj[Cmaj][Cmaj[Fmaj][Cmaj]]][Cmaj[Gmaj[Dmaj[Amin][Dmaj]][Gmaj]][Cmaj]]]]"
    "[Cmin[Cmin[Cmin[Cmin][Cmin[Gmaj][Cmin]]][Gmaj[Abmaj][Gmaj]]][Cmin[Cmin[Cmin][Cmin[Gmaj][Cmin]]][Cmin[Gmaj[Abmaj][Gmaj]][Cmin]]]]"
    "[Cmaj[Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Fmaj][Cmaj]]]][Cmaj[Cmaj[Fmaj][Cmaj]][Cmaj[Gmaj[Dmin[Amaj[Emin][Amaj]][Dmin]][Gmaj]][Cmaj]]]]"
    "[Cmaj[Cmaj[Cmaj][Cmaj[Bbmaj][Cmaj]]][Cmaj[Dbmaj[Dmaj[Ebmaj[Bbmaj][Ebmaj]][Dmaj]][Dbmaj]][Cmaj]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj][Cmaj[Gmaj[Dmaj][Gmaj[Dmin][Gmaj]]][Cmaj]]][Cmaj[Cmaj][Cmaj[Gmaj[Dmaj][Gmaj[Dmin][Gmaj]]][Cmaj]]]][Cmaj[Cmaj[Gmaj[Fmaj][Gmaj[Dmaj][Gmaj[Dmin][Gmaj]]]][Cmaj[Dbmaj][Cmaj]]][Cmaj[Gmaj[Dmaj][Gmaj[Dmin][Gmaj]]][Cmaj[Cmaj][Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]]]]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj][Amin[Emaj[Bhdim][Emaj]][Amin]]][Cmaj[Cmaj[Bbmaj[Fmaj[Cmaj[Gmin][Cmaj]][Fmaj]][Bbmaj]][Cmaj]][Gmaj[Dmaj[Amin][Dmaj]][Gmaj[Dmin][Gmaj]]]]][Cmaj[Cmaj[Cmaj[Cmaj][Amin[Emaj[Bhdim][Emaj]][Amin]]][Cmaj[Bbmaj[Fmaj[Cmaj[Gmin][Cmaj]][Fmaj]][Bbmaj]][Cmaj]]][Cmaj[Gmaj[Dmaj][Gmaj[Dmin[Amaj[Emin[Emin[D#dim][Emin]][Emin[Fmaj][Emin]]][Amaj]][Dmin]][Gmaj]]][Cmaj[Cmaj][Cmaj[Gmaj][Cmaj]]]]]]"
    "[Cmaj[Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Fmin][Gmaj[Dmin][Gmaj[Abmaj[Ebmaj][Abmaj]][Gmaj]]]]]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Cmaj[Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Fmin][Gmaj[Dmin][Gmaj]]]][Cmaj]]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj[Dbmaj][Cmaj[Gaug][Cmaj[Cmaj][Cmaj]]]][Cmaj[Galt[Fmin[Cmaj[Ghdim[Dmaj[Cmaj][Dmaj]][Ghdim]][Cmaj]][Fmin[Fmin][Fmin]]][Galt[Dhdim][Galt]]][Cmaj[Cmaj][Cmaj]]]][Cmaj[Galt[Fmin[Cmaj[Ghdim[Dmaj[Cmaj][Dmaj]][Ghdim]][Cmaj]][Fmin[Fmin][Fmin]]][Galt[Dhdim][Galt]]][Cmaj[Cmaj][Cmaj[Cmaj][Cmaj]]]]][Cmaj[Gmaj[Bbmaj[Fmaj[Cmin[Cmin][Cmin]][Fmaj]][Bbmaj[Bbmaj][Bbmaj]]][Gmaj[Abmaj[Abmaj][Abmaj]][Gmaj[Dmin][Gmaj]]]][Cmaj[Gaug[Fmin[Cmaj[Ghdim][Cmaj]][Fmin[Fmin][Fmin]]][Gaug[Dbmaj][Gaug]]][Cmaj[Cmaj[Cmaj][Cmaj[Cmaj][Cmaj]]][Cmaj[Dmaj][Cmaj]]]]]]"
    "[Cmaj[Cmaj[Cmaj][Cmaj[Cmin[Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Dmin][Gmaj]]][Cmin]][Cmaj[Gmaj[Dhdim][Gmaj]][Cmaj]]]][Cmaj[Cmaj[Gmaj[Gmaj[Dmin][Gmaj]][Gmaj[Dmin[Ebmaj[Emin][Ebmaj]][Dmin]][Gmaj]]][Cmaj]][Cmaj[Cmaj[Gmaj[Abmaj[Ebmaj[Bbmaj[Fmaj[Cmin][Fmaj]][Bbmaj]][Ebmaj[Bbmin][Ebmaj]]][Abmaj]][Gmaj[Dmin][Gmaj]]][Cmaj]][Cmaj[Cmaj[Gmaj[Gmaj[Dmin][Gmaj]][Gmaj[Dmin[Ebmaj[Emin][Ebmaj]][Dmin]][Gmaj]]][Cmaj]][Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]]]]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj][Cmaj[Fmaj][Cmaj]]][Cmaj[Fmaj][Cmaj]]][Cmaj[Gmaj][Cmaj[Fmaj][Cmaj]]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj][Amin[Emaj[Bhdim][Emaj]][Amin]]][Gmaj[Dmaj[Amaj[Bbmaj[Fmaj[Cmaj[Gmin][Cmaj]][Fmaj]][Bbmaj]][Amaj[Emin][Amaj]]][Dmaj]][Gmaj[Dmin][Gmaj]]]][Cmaj[Cmaj[Cmaj][Amin[Emaj[Bhdim][Emaj]][Amin]]][Gmaj[Abmaj[Ebmin[Bbmaj[Fmaj[Cmaj[Gmin][Cmaj]][Fmaj]][Bbmaj]][Ebmin[Emin][Ebmin]]][Abmaj]][Gmaj[Dmin][Gmaj]]]]]"
    "[Cmaj[Cmaj[Cmaj][Cmaj[Bbmaj[Fmin][Bbmaj]][Cmaj]]][Cmaj[Cmaj[Gmaj[Abmaj[Ebmaj[Bbmin][Ebmaj]][Abmaj]][Gmaj[Dmaj[Amin][Dmaj]][Gmaj[Dmin][Gmaj]]]][Cmaj]][Dbmaj[Abmaj[Ebmaj][Abmaj]][Dbmaj]]]]"
    "[Cmaj[Cmaj[Cmaj[Cmaj[Cmaj][Cmaj[Gmaj[Dmin[Amin][Dmin]][Gmaj]][Cmaj]]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Gmaj[Gmaj[Dmin][Gmaj]][Gmaj[Dmin[Amaj[Emin][Amaj]][Dmin]][Gmaj]]]]][Gmaj[Gmaj[Dmaj[Ebmaj[Ebmaj][Ebmaj[Bbmaj[Fmin][Bbmaj]][Ebmaj[Ebmaj][Ebmaj]]]][Dmaj]][Gmaj]][Gmaj[Dmin[Ebmaj[Emin][Ebmaj]][Dmin]][Gmaj]]]][Cmaj[Cmaj[Cmaj[Cmaj][Cmaj[Gmaj[Dmin[Amin][Dmin]][Gmaj]][Cmaj]]][Cmaj[Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]][Gmaj[Gmaj[Dmin][Gmaj]][Gmaj[Dmin[Amaj[Emin][Amaj]][Dmin]][Gmaj]]]]][Cmaj[Cmaj[Dmin[Fmaj[Cmaj[Gmin][Cmaj]][Fmaj]][Dmin]][Cmaj]][Cmaj[Cmaj[Gmaj[Dmin[Amin][Dmin]][Gmaj]][Cmaj]][Cmaj[Gmaj[Dmin][Gmaj]][Cmaj]]]]]]"])
    for i in 1:length(trees)
        trees[i] = map(translate_goldstandard_notation1, trees[i])
    end
    trees
end

function read_parsable_goldstandard_trees_1()
    map(tree, ["[0maj[0maj[0maj[7maj[2min][7maj]][0maj]][0maj[7maj[2min[9maj][2min]][7maj]][0maj]]][0maj[0maj[7maj[2min][7maj]][0maj]][0maj[7maj[8maj[3maj[10maj[5maj[0min][5maj]][10maj]][3maj[10min][3maj]]][8maj]][7maj[2min][7maj]]][0maj]]]]",
    "[0min[0min[0min[0min][0min[7maj][0min]]][7maj[8maj][7maj]]][0min[0min[0min][0min[7maj][0min]]][0min[7maj[8maj][7maj]][0min]]]]",
    "[0maj[0maj[0maj[7maj[2min][7maj]][0maj]][0maj[0maj[7maj[2min][7maj]][0maj]][0maj[5maj][0maj]]]][0maj[0maj[5maj][0maj]][0maj[7maj[2min[9maj[4min][9maj]][2min]][7maj]][0maj]]]]",
    "[0maj[0maj[0maj][0maj[10maj][0maj]]][0maj[1maj[2maj[3maj[10maj][3maj]][2maj]][1maj]][0maj]]]",
    "[0maj[0maj[0maj[0maj][0maj[7maj[2maj][7maj[2min][7maj]]][0maj]]][0maj[0maj][0maj[7maj[2maj][7maj[2min][7maj]]][0maj]]]][0maj[0maj[7maj[5maj][7maj[2maj][7maj[2min][7maj]]]][0maj[1maj][0maj]]][0maj[7maj[2maj][7maj[2min][7maj]]][0maj[0maj][0maj[7maj[2min][7maj]][0maj]]]]]]",
    "[0maj[0maj[0maj][0maj[0min[7maj[8maj[3maj[10maj[5maj[0min][5maj]][10maj]][3maj[10min][3maj]]][8maj]][7maj[2min][7maj]]][0min]][0maj[7maj[2hdim][7maj]][0maj]]]][0maj[0maj[7maj[7maj[2min][7maj]][7maj[2min[3maj[4min][3maj]][2min]][7maj]]][0maj]][0maj[0maj[7maj[8maj[3maj[10maj[5maj[0min][5maj]][10maj]][3maj[10min][3maj]]][8maj]][7maj[2min][7maj]]][0maj]][0maj[0maj[7maj[7maj[2min][7maj]][7maj[2min[3maj[4min][3maj]][2min]][7maj]]][0maj]][0maj[7maj[2min][7maj]][0maj]]]]]]",
    "[0maj[0maj[0maj[0maj][0maj[5maj][0maj]]][0maj[5maj][0maj]]][0maj[7maj][0maj[5maj][0maj]]]]"])
end

function read_goldstandard_trees_2()
    trees = map(tree,
    # ["[C[C[C[C[C[C][A-]][C[Db[Ab][Db]]C]][C[Db[D[Eb][D]]Db][C[C][A-]]]][C[Db[Ab][Db]]C]][C[C[G[F[F[C[G[B][G[Ab][G]]][C]][F[F][D-]]][F[Gb[Db][Gb]][F]]][G[Ab][G]]][C[Db][C[C][A-]]]][C[C[Db[Ab][Db]][C]][C[Db][C[G[Ab][G]][C]]]]]]",
    ["[C-[C-[C-[C-[G[D%][G]][C-]][Eb[Bb[F-[C][F-]][Bb]][Eb]]][C-[C-[G[D%][G]][C-]][C-[G[F-][G[Ab][G]]][C-]]]][C-[G[D%[A][D%]][G]][C-[C-[C-[C-[G[D%][G]][C-]][Eb[Bb[F-[C][F-]][Bb]][Eb]]][C-[C-[G[D%][G]][C-]][C-[G[F-][G[Ab][G]]][C-]]]] [C-[G[Ab[A][Ab]][G]][C-]]]]]",
    "[C-[C-[C-[Eb[Bb[F-][Bb]][Eb]][C-[G[D%[Ab][D%]][G]][C-]]][C-[Eb[Bb[F-[C][F-]][Bb]][Eb]][C-[G[D%[Ab][D%]][G]][C-]]]][C-[C-[C-[G[D%][G]][C-]][Eb[Bb[F-][Bb]][Eb]]][C-[C-[G[D%][G]][C-]][C-[G[D%[Ab[Bbb[Bb-[Cb][Bb-]][Bbb]][Ab]][D%]][G]][C-]]]]]",
    "[C-[C-[C-[C-[C-[C-][C-[G[D%][G]][C-]]][C-[C-[G[D%][G]][C-]][Eb[Bb[F-][Bb]][Eb]]]][C-[Eb[Bb[F-[Eo][F-[C][F-]]][Bb]][Eb]][C-[G[D%[Ab][D%]][G]][C-]]]][C-[C-[C-[G[D%][G]][C-]][C-[G[D%][G]][C-]]][C-[G[F-[C[G[D%][G]][C[G%][C]]][F-]][G[D%[Ab[C-[G[F-][G[D%][G]]][C-]][Ab]][D%]][G]]][C-]]]][C-[C-[C-[F-][C-]][C-[F-][C-]]][C-[G-[F-][G-]][C-]]]]",
    "[C-[C-[C-[C-][C-[G[F-][G[D%][G]]][C-]]][C-[Db[Ab[Eb-][Ab]][Db]][C-[G[D%][G]][C-]]]][C-[G[D%][G]][C-]]]",
    "[C[C[C][C[C[G[F[C[G-[D[A-[E[B%][E]][A-]][D]][G-]][C]][F]][G[Ab[A[Bb[F-][Bb]][A[E-][A]]][Ab[Eb-][Ab]]][G[D-][G]]]][C]][A-]]][C[G[D-][G]][C]]]",
    # "[C[C[C[C][C[Bb[D-[A[Bb][A]][D-]][Bb[F-][Bb]]][C[G[G[D-[A-[E-][A-]][D-]][G]][G[D[A-[E[B%][E]][A-]]D]][G[D-][G]]][C]]]][C[Bb[D-[A[Bb][A]][D-]][Bb[F-][Bb]]][C[G[D-[A-[A-[E-][A-]][A-[E-[F[F#%][F]][E-]][A-]]][D-]][G]][C[C][C[F][C]]]]]][C[G[G[G[D-[A-[E-][A-]][D-]][G]][G[D-[A-[E-][A-]][D-]][G]]][G[D-[A-[E-][A-]][D-]][G]]][C[C][C[F][C]]]]]",
    "[C-[C-[C-][C-[F-][C-]]][C-[Db[D][Db]][C-]]]",
    "[C[C[C[C[G[D-[A-][D-]][G]][C]][A-[E[B%[F][B%]][E]][A-]]][C[C[G[D-[A][D-]][G]][C]][C[G[D-][G]][C]]]][C[C[C[G[D-[A-[E[B-][E]][A-]][D-]][G]][C]][A-[E[B%[F][B%]][E]][A-]]][C[C[G[D-[A][D-]][G]][C]][C[G[F[C][F]][G]][C]]]]]",])
    for i in 1:length(trees)
        # println(i)
        trees[i] = map(translate_goldstandard_notation2, trees[i])
    end
    trees
end

function read_parsable_goldstandard_trees_2()
    map(tree, [
    "[0min[0min[0min[0min[7maj[2hdim][7maj]][0min]][3maj[10maj[5min[0maj][5min]][10maj]][3maj]]][0min[0min[7maj[2hdim][7maj]][0min]][0min[7maj[5min][7maj[8maj][7maj]]][0min]]]][0min[7maj[2hdim[9maj][2hdim]][7maj]][0min[0min[0min[0min[7maj[2hdim][7maj]][0min]][3maj[10maj[5min[0maj][5min]][10maj]][3maj]]][0min[0min[7maj[2hdim][7maj]][0min]][0min[7maj[5min][7maj[8maj][7maj]]][0min]]]][0min[7maj[8maj[9maj][8maj]][7maj]][0min]]]]]",
    "[0min[0min[0min[0min][0min[7maj[5min][7maj[2hdim][7maj]]][0min]]][0min[1maj[8maj[3min][8maj]][1maj]][0min[7maj[2hdim][7maj]][0min]]]][0min[7maj[2hdim][7maj]][0min]]]",
    "[0maj[0maj[0maj][0maj[0maj[7maj[5maj[0maj[7min[2maj[9min[4maj[11hdim][4maj]][9min]][2maj]][7min]][0maj]][5maj]][7maj[8maj[9maj[10maj[5min][10maj]][9maj[4min][9maj]]][8maj[3min][8maj]]][7maj[2min][7maj]]]][0maj]][9min]]][0maj[7maj[2min][7maj]][0maj]]]",
    "[0min[0min[0min][0min[5min][0min]]][0min[1maj[2maj][1maj]][0min]]]",
    "[0maj[0maj[0maj[0maj[7maj[2min[9min][2min]][7maj]][0maj]][9min[4maj[11hdim[5maj][11hdim]][4maj]][9min]]][0maj[0maj[7maj[2min[9maj][2min]][7maj]][0maj]][0maj[7maj[2min][7maj]][0maj]]]][0maj[0maj[0maj[7maj[2min[9min[4maj[11min][4maj]][9min]][2min]][7maj]][0maj]][9min[4maj[11hdim[5maj][11hdim]][4maj]][9min]]][0maj[0maj[7maj[2min[9maj][2min]][7maj]][0maj]][0maj[7maj[5maj[0maj][5maj]][7maj]][0maj]]]]]",
    ])
end

# read_goldstandard_trees_2()

function train!(grammar, trees)
    for tree in trees
        forest = run_chartparser(tree, grammar)
        add_obs!(grammar, best_tree(forest))
    end
    grammar
end

function single_fold_cross_validation(trees, make_grammar=make_jazz_grammar, with_training=true)
    precision_vals = zeros(Float64, length(trees))
    for (i,test_tree) in enumerate(trees)
        grammar = make_grammar()
        if with_training
            train!(grammar, [trees[1:i-1];trees[i+1:end]])
        end
        @time forest  = run_chartparser(leaf_data(test_tree), grammar)
        precision_vals[i] = precision(best_tree(forest), test_tree)
    end
    precision_vals
end

function single_fold_cross_validation_baseline(trees)
    precision_vals = zeros(Float64, length(trees))
    for (i,test_tree) in enumerate(trees)
        baseline_tree = right_branching_tree(length(leafs(test_tree)))
        precision_vals[i] = precision(baseline_tree, test_tree)
    end
    precision_vals
end

# function single_fold_cross_validation(trees)
#     precision_vals = Float64[]
#     recall_vals = Float64[]
#     precision_base = Float64[]
#     recall_base = Float64[]
#
#     for (i,test_tree) in enumerate(trees)
#         # grammar = make_jazz_grammar()
#         grammar = make_plain_PCFG_jazz_grammar()
#
#         # training
#         for tree in [trees[1:i-1];trees[i+1:end]]
#             @time forest = run_chartparser(tree, grammar)
#             add_obs!(grammar, best_tree(forest))
#         end
#
#         # evaluating
#         @time forest  = run_chartparser(leaf_data(test_tree), grammar)
#         parsed_tree   = best_tree(forest)
#         baseline_tree = right_branching_tree(length(leafs(test_tree)))
#         push!(precision_vals, precision(parsed_tree,   test_tree))
#         push!(recall_vals,       recall(parsed_tree,   test_tree))
#         push!(precision_base, precision(baseline_tree, test_tree))
#         push!(recall_base,       recall(baseline_tree, test_tree))
#     end
#     precision_vals, recall_vals, precision_base, recall_base
# end

function parse_leaf_data(trees)
    grammar = make_jazz_grammar(Int)
    # grammar = make_plain_PCFG_jazz_grammar(BigInt)
    for k in 1:2
        println("id\tlength\ttime_in_seconds\t number_of_parses")
        for (i,t) in enumerate(trees)
            chords = leaf_data(t)
            tic()
            forest = run_chartparser(chords, grammar)
            time = toq()
            @assert is_complete(forest)
            println(i, "\t", length(chords), "\t", time, "\t", score(forest))
        end
    end
end

function make_parsable_trees()
    parsable_trees = [
        read_treebank_trees();
        read_parsable_goldstandard_trees_1();
        read_parsable_goldstandard_trees_2()
    ]
    [parsable_trees[1:6]; parsable_trees[8:end]]
end

# parse_leaf_data(parsable_trees)

# p,r,pb,rb = single_fold_cross_validation(parsable_trees)

# parse_trees(parsable_trees)


text:
using Arrow, CSV

csvs = filter(readdir(".")) do f
    return last(splitext(f)) == ".csv"
end

for f in csvs
    Arrow.write("$(first(splitext(f))).arrow", CSV.File(f); compress=:zstd)
end


text:
# patented simplex noise algorithm

# will remove this soon

struct Grad
  x::Int64
  y::Int64
  z::Int64
end

const grad3 = [Grad(1, 1, 0), Grad(-1, 1, 0), Grad(1, -1, 0), Grad(-1, -1, 0), Grad(1, 0, 1),
  Grad(-1, 0, 1), Grad(1, 0, -1), Grad(-1, 0, -1), Grad(0, 1, 1), Grad(0, -1, 1),
  Grad(0, 1, -1), Grad(0, -1, -1)]

p = [151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225,
140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120,
234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 88,
237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139,
48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105,
92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73,
209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86,
164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38,
147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189,
28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153,
101, 155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224,
232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144,
12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214,
31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254,
138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215,
61, 156, 180]

# To remove the need for index wrapping, double the _permutation table length

const _perm =      Array{Int64}(undef, 512)
const _permMod12 = Array{Int64}(undef, 512)

function initpatentednoise()
    for i in 1:512
          _perm[i] = p[max(1, i & 255)] # julia indexing
          _permMod12[i] = mod1(_perm[i], 12) # julia indexing
    end
end

initpatentednoise()

function quickfloor(x)
    return convert(Int, floor(x))
end

function dotp(g::Grad, x, y)
    return g.x * x + g.y * y
end

function dotp(g::Grad, x, y, z)
    return g.x * x + g.y * y + g.z * z
end

# based on code by Stefan Gustavson, Linköping University
function _noise3(xin, yin, zin)
    # Skewing and unskewing factors
    _G3 = 1.0/6.0
    # n0, n1, n2, n3 are noise contributions from the four corners
    # skew the input space to determine which simplex cell we're in
    s = (xin + yin + zin) * 1.0/3.0 # simple skew factor for 3D
    i = quickfloor(xin + s)
    j = quickfloor(yin + s)
    k = quickfloor(zin + s)
    t = (i + j + k) * _G3
    X0 = i - t # Unskew the cell origin back to (x, y, z) space
    Y0 = j - t
    Z0 = k - t
    x0 = xin - X0 # the x, y, z distances from the cell origin
    y0 = yin - Y0
    z0 = zin - Z0
    # For the 3D case, the simplex shape is a slightly irregular tetrahedron.
    # Determine which simplex we are in.
    # i1, j1, k1 are offsets for second corner of simplex in (i, j, k) coords
    # i2, j2, k2 are offsets for third corner of simplex in (i, j, k) coords
    if x0 >= y0
        if y0 >= z0
            i1=1
            j1=0
            k1=0
            i2=1
            j2=1
            k2=0 # X Y Z order
        elseif x0 >= z0
            i1=1
            j1=0
            k1=0
            i2=1
            j2=0
            k2=1 # X Z Y order
        else
            i1=0
            j1=0
            k1=1
            i2=1
            j2=0
            k2=1 # Z X Y order
        end
    else  # x0 < y0
        if y0 < z0
            i1=0
            j1=0
            k1=1
            i2=0
            j2=1
            k2=1 # Z Y X order
        elseif x0 < z0
            i1=0
            j1=1
            k1=0
            i2=0
            j2=1
            k2=1 # Y Z X order
        else
            i1=0
            j1=1
            k1=0
            i2=1
            j2=1
            k2=0 # Y X Z order
        end
    end
    # A step of (1, 0, 0) in (i, j, k) means a step of (1-c, -c, -c) in (x, y, z),
    # a step of (0, 1, 0) in (i, j, k) means a step of (-c, 1-c, -c) in (x, y, z), and
    # a step of (0, 0, 1) in (i, j, k) means a step of (-c, -c, 1-c) in (x, y, z), where
    # c = 1/6.
    x1 = x0 - i1 + _G3 # Offsets for second corner in (x, y, z) coords
    y1 = y0 - j1 + _G3
    z1 = z0 - k1 + _G3
    x2 = x0 - i2 + 2.0 * _G3 # Offsets for third corner in (x, y, z) coords
    y2 = y0 - j2 + 2.0 * _G3
    z2 = z0 - k2 + 2.0 * _G3
    x3 = x0 - 1.0 + 3.0 * _G3 # Offsets for last corner in (x, y, z) coords
    y3 = y0 - 1.0 + 3.0 * _G3
    z3 = z0 - 1.0 + 3.0 * _G3

    # Work out the hashed gradient indices of the four simplex corners
    ii = mod1(i & 255, 256)
    jj = mod1(j & 255, 256)
    kk = mod1(k & 255, 256)
    gi0 = _permMod12[mod1(ii + _perm[jj + _perm[kk]], 256)]
    gi1 = _permMod12[mod1(ii + i1 + _perm[jj + j1 + _perm[kk + k1]], 256)]
    gi2 = _permMod12[mod1(ii + i2 + _perm[jj + j2 + _perm[kk + k2]], 256)]
    gi3 = _permMod12[mod1(ii + 1 +  _perm[jj + 1 + _perm[kk + 1]], 256)]
    # Calculate the contribution from the four corners
    t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0
    if t0 < 0
        n0 = 0.0
    else
        t0 *= t0
        n0 = t0 * t0 * dotp(grad3[gi0], x0, y0, z0)
    end
    t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1
    if t1 < 0
        n1 = 0.0
    else
        t1 *= t1
        n1 = t1 * t1 * dotp(grad3[gi1], x1, y1, z1)
    end
    t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2
    if t2 < 0
        n2 = 0.0
    else
        t2 *= t2
        n2 = t2 * t2 * dotp(grad3[gi2], x2, y2, z2)
    end
    t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3
    if t3 < 0
        n3 = 0.0
    else
        t3 *= t3
        n3 = t3 * t3 * dotp(grad3[gi3], x3, y3, z3)
    end
    # Add contributions from each corner to get the final noise value.
    # The result is scaled to stay just inside [-1, 1]
    res = 32.0 * (n0 + n1 + n2 + n3)
    # then to [0, 1]
    return (res + 1) / 2
end

function _octaves(x, y, z, octaves::Int, persistence)
    total     = 0.0
    frequency = 1.0
    amplitude = 1.0
    maxval    = 0.0
    for i in 1:octaves
        total += _noise3(x * frequency, y * frequency, z * frequency) * amplitude
        maxval += amplitude
        amplitude *= persistence
        frequency *= 2
    end
    return total / maxval
end

function seedpatentednoise(a::Array{Int64})
    _permMod12[:] = a
end

"""
    patentednoise(x, y = 1, z = 1;
        detail::Int64 = 1,
        persistence   = 0.0)

Generate a noise value between 0.0 and 1.0 corresponding to the `x`, `y`, and
`z` values.

The `detail` value is an integer specifying how many octaves of noise you want.

The `persistence` value, typically between 0.0 and 1.0, controls how quickly the
amplitude diminishes for each successive octave for values of `detail` greater
than 1.
"""
function patentednoise(x, y=1, z=1;
    detail      = 1,
    persistence = 0.0)
    return _octaves(x, y, z, detail, persistence)
end


text:
abstract type CopperPlatePowerModel <: PM.AbstractActivePowerFormulation  end

abstract type StandardPTDFForm <: PM.DCPlosslessForm end

#= This code is from PowerModel's network definitions. Added here just for reference.

export
    # exact non-convex models
    ACPPowerModel, StandardACPForm,
    ACRPowerModel, StandardACRForm,
    ACTPowerModel, StandardACTForm,

    # linear approximations
    DCPPowerModel, DCPlosslessForm,
    NFAPowerModel, NFAForm,

    # quadratic approximations
    DCPLLPowerModel, StandardDCPLLForm,
    LPACCPowerModel, AbstractLPACCForm,

    # quadratic relaxations
    SOCWRPowerModel, SOCWRForm,
    QCWRPowerModel, QCWRForm,
    SOCWRConicPowerModel, SOCWRConicForm,
    QCWRTriPowerModel, QCWRTriForm,
    SOCBFPowerModel, SOCBFForm,
    SOCBFConicPowerModel, SOCBFConicForm,

    # sdp relaxations
    SDPWRMPowerModel, SDPWRMForm,
    SparseSDPWRMPowerModel, SparseSDPWRMForm

##### Top Level Abstract Types #####

"active power only models"
abstract type AbstractActivePowerFormulation <: AbstractPowerFormulation end

"variants that target conic solvers"
abstract type AbstractConicPowerFormulation <: AbstractPowerFormulation end

"for branch flow models"
abstract type AbstractBFForm <: AbstractPowerFormulation end

"for variants of branch flow models that target QP or NLP solvers"
abstract type AbstractBFQPForm <: AbstractBFForm end

"for variants of branch flow models that target conic solvers"
abstract type AbstractBFConicForm <: AbstractBFForm end

##### Exact Non-Convex Models #####

""
abstract type AbstractACPForm <: AbstractPowerFormulation end

""
abstract type StandardACPForm <: AbstractACPForm end

"""
AC power flow formulation with polar bus voltage variables.

The seminal reference of AC OPF:
```
@article{carpentier1962contribution,
  title={Contribution to the economic dispatch problem},
  author={Carpentier, J},
  journal={Bulletin de la Societe Francoise des Electriciens},
  volume={3},
  number={8},
  pages={431--447},
  year={1962}
}
```

History and discussion:
```
@techreport{Cain2012,
  author = {Cain, Mary B and {O' Neill}, Richard P and Castillo, Anya},
  title = {{History of optimal power flow and formulations}},
  year = {2012}
  pages = {1--36},
  url = {https://www.ferc.gov/industries/electric/indus-act/market-planning/opf-papers/acopf-1-history-formulation-testing.pdf}
}
```
"""
const ACPPowerModel = GenericPowerModel{StandardACPForm}

""
ACPPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, StandardACPForm; kwargs...)

""
abstract type AbstractACRForm <: AbstractPowerFormulation end

""
abstract type StandardACRForm <: AbstractACRForm end

"""
AC power flow formulation with rectangular bus voltage variables.

```
@techreport{Cain2012,
  author = {Cain, Mary B and {O' Neill}, Richard P and Castillo, Anya},
  pages = {1--36},
  title = {{History of optimal power flow and formulations}},
  url = {https://www.ferc.gov/industries/electric/indus-act/market-planning/opf-papers/acopf-1-history-formulation-testing.pdf}
  year = {2012}
}
```
"""
const ACRPowerModel = GenericPowerModel{StandardACRForm}

"default rectangular AC constructor"
ACRPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, StandardACRForm; kwargs...)

""
abstract type AbstractACTForm <: AbstractPowerFormulation end

""
abstract type StandardACTForm <: AbstractACTForm end

"""
AC power flow formulation (nonconvex) with variables for voltage angle, voltage magnitude squared, and real and imaginary part of voltage crossproducts. A tangens constraint is added to represent meshed networks in an exact manner.
```
@ARTICLE{4349090,
  author={R. A. Jabr},
  title={A Conic Quadratic Format for the Load Flow Equations of Meshed Networks},
  journal={IEEE Transactions on Power Systems},
  year={2007},
  month={Nov},
  volume={22},
  number={4},
  pages={2285-2286},
  doi={10.1109/TPWRS.2007.907590},
  ISSN={0885-8950}
}
```
"""
const ACTPowerModel = GenericPowerModel{StandardACTForm}

"default AC constructor"
ACTPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, StandardACTForm; kwargs...)

##### Linear Approximations #####

""
abstract type AbstractDCPForm <: AbstractActivePowerFormulation end

"active power only formulations where p[(i,j)] = -p[(j,i)]"
abstract type DCPlosslessForm <: AbstractDCPForm end

"""
Linearized 'DC' power flow formulation with polar voltage variables.

```
@ARTICLE{4956966,
  author={B. Stott and J. Jardim and O. Alsac},
  journal={IEEE Transactions on Power Systems},
  title={DC Power Flow Revisited},
  year={2009},
  month={Aug},
  volume={24},
  number={3},
  pages={1290-1300},
  doi={10.1109/TPWRS.2009.2021235},
  ISSN={0885-8950}
}
```
"""
const DCPPowerModel = GenericPowerModel{DCPlosslessForm}

"default DC constructor"
DCPPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, DCPlosslessForm; kwargs...)

abstract type NFAForm <: DCPlosslessForm end

"""
The an active power only network flow approximation, also known as the transportation model.
"""
const NFAPowerModel = GenericPowerModel{NFAForm}

"default DC constructor"
NFAPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, NFAForm; kwargs...)

##### Quadratic Approximations #####

""
abstract type AbstractDCPLLForm <: AbstractDCPForm end

""
abstract type StandardDCPLLForm <: AbstractDCPLLForm end

""
const DCPLLPowerModel = GenericPowerModel{StandardDCPLLForm}

"default DC constructor"
DCPLLPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, StandardDCPLLForm; kwargs...)

""
abstract type AbstractLPACForm <: AbstractPowerFormulation end

abstract type AbstractLPACCForm <: AbstractLPACForm end

"""
The LPAC Cold-Start AC Power Flow Approximation.

Note that the LPAC Cold-Start model requires the least amount of information
but is also the least accurate variant of the LPAC formulations.  If a
nominal AC operating point is available, the LPAC Warm-Start model will provide
improved accuracy.

The original publication suggests to use polyhedral outer approximations for
the cosine and line thermal lit constraints.  Given the recent improvements in
MIQCQP solvers, this implementation uses quadratic functions for those
constraints.

```
@article{doi:10.1287/ijoc.2014.0594,
  author = {Coffrin, Carleton and Van Hentenryck, Pascal},
  title = {A Linear-Programming Approximation of AC Power Flows},
  journal = {INFORMS Journal on Computing},
  volume = {26},
  number = {4},
  pages = {718-734},
  year = {2014},
  doi = {10.1287/ijoc.2014.0594},
  eprint = {https://doi.org/10.1287/ijoc.2014.0594}
}
```
"""
const LPACCPowerModel = GenericPowerModel{AbstractLPACCForm}

"default LPACC constructor"
LPACCPowerModel(data::Dict{String,Any}; kwargs...) =
    GenericPowerModel(data, AbstractLPACCForm; kwargs...)

##### Quadratic Relaxations #####

""
abstract type AbstractWRForm <: AbstractPowerFormulation end

""
abstract type AbstractWRConicForm <: AbstractConicPowerFormulation end

""
abstract type SOCWRConicForm <: AbstractWRConicForm end

""
abstract type SOCWRForm <: AbstractWRForm end

"""
Second-order cone relaxation of bus injection model of AC OPF.

The implementation casts this as a convex quadratically constrained problem.
```
@article{1664986,
  author={R. A. Jabr},
  title={Radial distribution load flow using conic programming},
  journal={IEEE Transactions on Power Systems},
  year={2006},
  month={Aug},
  volume={21},
  number={3},
  pages={1458-1459},
  doi={10.1109/TPWRS.2006.879234},
  ISSN={0885-8950}
}
```
"""
const SOCWRPowerModel = GenericPowerModel{SOCWRForm}

"""
Second-order cone relaxation of bus injection model of AC OPF.

This implementation casts the problem as a convex conic problem.
"""
const SOCWRConicPowerModel = GenericPowerModel{SOCWRConicForm}

"default SOC constructor"
SOCWRPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, SOCWRForm; kwargs...)

SOCWRConicPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, SOCWRConicForm; kwargs...)

""
abstract type QCWRForm <: AbstractWRForm end

"""
"Quadratic-Convex" relaxation of AC OPF
```
@Article{Hijazi2017,
  author="Hijazi, Hassan and Coffrin, Carleton and Hentenryck, Pascal Van",
  title="Convex quadratic relaxations for mixed-integer nonlinear programs in power systems",
  journal="Mathematical Programming Computation",
  year="2017",
  month="Sep",
  volume="9",
  number="3",
  pages="321--367",
  issn="1867-2957",
  doi="10.1007/s12532-016-0112-z",
  url="https://doi.org/10.1007/s12532-016-0112-z"
}
```
"""
const QCWRPowerModel = GenericPowerModel{QCWRForm}

"default QC constructor"
QCWRPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, QCWRForm; kwargs...)

""
abstract type QCWRTriForm <: QCWRForm end

"""
"Quadratic-Convex" relaxation of AC OPF with convex hull of triple product
```
@Article{Hijazi2017,
  author="Hijazi, Hassan and Coffrin, Carleton and Hentenryck, Pascal Van",
  title="Convex quadratic relaxations for mixed-integer nonlinear programs in power systems",
  journal="Mathematical Programming Computation",
  year="2017",
  month="Sep",
  volume="9",
  number="3",
  pages="321--367",
  issn="1867-2957",
  doi="10.1007/s12532-016-0112-z",
  url="https://doi.org/10.1007/s12532-016-0112-z"
}
```
"""
const QCWRTriPowerModel = GenericPowerModel{QCWRTriForm}

"default QC trilinear model constructor"
QCWRTriPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, QCWRTriForm; kwargs...)

""
abstract type SOCBFForm <: AbstractBFQPForm end

"""
Second-order cone relaxation of branch flow model

The implementation casts this as a convex quadratically constrained problem.
```
@INPROCEEDINGS{6425870,
  author={M. Farivar and S. H. Low},
  title={Branch flow model: Relaxations and convexification},
  booktitle={2012 IEEE 51st IEEE Conference on Decision and Control (CDC)},
  year={2012},
  month={Dec},
  pages={3672-3679},
  doi={10.1109/CDC.2012.6425870},
  ISSN={0191-2216}
}
```
Extended as discussed in:
```
@misc{1506.04773,
  author = {Carleton Coffrin and Hassan L. Hijazi and Pascal Van Hentenryck},
  title = {DistFlow Extensions for AC Transmission Systems},
  year = {2018},
  eprint = {arXiv:1506.04773},
  url = {https://arxiv.org/abs/1506.04773}
}
```
"""
const SOCBFPowerModel = GenericPowerModel{SOCBFForm}

"default SOC constructor"
SOCBFPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, SOCBFForm; kwargs...)

""
abstract type SOCBFConicForm <: AbstractBFConicForm end

""
const SOCBFConicPowerModel = GenericPowerModel{SOCBFConicForm}

"default SOC constructor"
SOCBFConicPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, SOCBFConicForm; kwargs...)

###### SDP Relaxations ######

""
abstract type AbstractWRMForm <: AbstractConicPowerFormulation end

""
abstract type SDPWRMForm <: AbstractWRMForm end

"""
Semi-definite relaxation of AC OPF

Originally proposed by:
```
@article{BAI2008383,
  author = "Xiaoqing Bai and Hua Wei and Katsuki Fujisawa and Yong Wang",
  title = "Semidefinite programming for optimal power flow problems",
  journal = "International Journal of Electrical Power & Energy Systems",
  volume = "30",
  number = "6",
  pages = "383 - 392",
  year = "2008",
  issn = "0142-0615",
  doi = "https://doi.org/10.1016/j.ijepes.2007.12.003",
  url = "http://www.sciencedirect.com/science/article/pii/S0142061507001378",
}
```
First paper to use "W" variables in the BIM of AC OPF:
```
@INPROCEEDINGS{6345272,
  author={S. Sojoudi and J. Lavaei},
  title={Physics of power networks makes hard optimization problems easy to solve},
  booktitle={2012 IEEE Power and Energy Society General Meeting},
  year={2012},
  month={July},
  pages={1-8},
  doi={10.1109/PESGM.2012.6345272},
  ISSN={1932-5517}
}
```
"""
const SDPWRMPowerModel = GenericPowerModel{SDPWRMForm}

""
SDPWRMPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, SDPWRMForm; kwargs...)

abstract type SparseSDPWRMForm <: SDPWRMForm end

"""
Sparsity-exploiting semidefinite relaxation of AC OPF

Proposed in:
```
@article{doi:10.1137/S1052623400366218,
  author = {Fukuda, M. and Kojima, M. and Murota, K. and Nakata, K.},
  title = {Exploiting Sparsity in Semidefinite Programming via Matrix Completion I: General Framework},
  journal = {SIAM Journal on Optimization},
  volume = {11},
  number = {3},
  pages = {647-674},
  year = {2001},
  doi = {10.1137/S1052623400366218},
  URL = {https://doi.org/10.1137/S1052623400366218},
  eprint = {https://doi.org/10.1137/S1052623400366218}
}
```
Original application to OPF by:
```
@ARTICLE{6064917,
  author={R. A. Jabr},
  title={Exploiting Sparsity in SDP Relaxations of the OPF Problem},
  journal={IEEE Transactions on Power Systems},
  volume={27},
  number={2},
  pages={1138-1139},
  year={2012},
  month={May},
  doi={10.1109/TPWRS.2011.2170772},
  ISSN={0885-8950}
}
```
"""
const SparseSDPWRMPowerModel = GenericPowerModel{SparseSDPWRMForm}

""
SparseSDPWRMPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, SparseSDPWRMForm; kwargs...)

##### Union Types #####
#
# These types should not be exported because they exist only to prevent code
# replication
#
# Note that Union types are discouraged in Julia,
# https://docs.julialang.org/en/v1/manual/style-guide/#Avoid-strange-type-Unions-1
# and should be used with discretion.
#
# If you are about to add a union type, first double check if refactoring the
# type hierarchy can resolve the issue instead.
#

AbstractWRForms = Union{AbstractACTForm, AbstractWRForm, AbstractWRConicForm, AbstractWRMForm}
AbstractWForms = Union{AbstractWRForms, AbstractBFForm}
AbstractPForms = Union{AbstractACPForm, AbstractACTForm, AbstractDCPForm, AbstractLPACForm}

"union of all conic form branches"
AbstractConicForms = Union{AbstractConicPowerFormulation, AbstractBFConicForm}

=#


text:
__precompile__()

module TestModule1

export func1, func2, test_interp, print_nprocs

using Distributed
using SharedArrays
using ProgressMeter
using Interpolations

include("SubModule.jl")
using .SubMod1

function func1()
    n = 200000
    arr = SharedArray{Float64}(n)
    @sync @distributed for i = 1:n
        arr[i] = i^2
    end
    res = sum(arr)
    return res
end

function func2()
    n = 100
    arr = SharedArray{Float64}(n)
    @sync @distributed for i = 1:n
        arr[i] = SubMod1.test_me(i)
    end
    res = sum(arr)
    return res
end

function sine_func(x)
    sleep(22*1e-3)
    return sin(x)
end

function create_interp()
    @info "create_interp"
    x = range(0, stop = 10, length = 1000)

    @info "Computing Interpolation..."
    y = @showprogress pmap(sine_func, x, batch_size=10)

    knots = (x,)
    itp = interpolate(knots, y, Gridded(Linear()))

    return itp
end

const test_interp = create_interp()

function print_nprocs()
    return nworkers()
end

end # module

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples