text:
export MPO,
       randomMPO,
       applyMPO,
       multMPO,
       errorMPOprod,
       maxlinkdim,
       orthogonalize!,
       truncate!

mutable struct MPO
  N_::Int
  A_::Vector{ITensor}
  llim_::Int
  rlim_::Int

  MPO() = new(0,Vector{ITensor}(), 0, 0)

  function MPO(N::Int, A::Vector{<:ITensor}, llim::Int=0, rlim::Int=N+1)
    new(N, A, llim, rlim)
  end

  function MPO(sites::Vector{Index})
    N = length(sites)
    v = Vector{ITensor}(undef, N)
    l = [Index(1, "Link,l=$ii") for ii ∈ 1:N-1]
    for ii ∈ eachindex(sites)
      s = sites[ii]
      sp = prime(s)
      if ii == 1
        v[ii] = ITensor(s, sp, l[ii])
      elseif ii == N
        v[ii] = ITensor(l[ii-1], s, sp)
      else
        v[ii] = ITensor(l[ii-1], s, sp, l[ii])
      end
    end
    new(N,v,0,N+1)
  end

end

MPO(N::Int) = MPO(N,Vector{ITensor}(undef,N))

function MPO(sites,
             ops::Vector{String})
  N = length(sites)
  its = Vector{ITensor}(undef, N)
  links = Vector{Index}(undef, N)
  for ii ∈ eachindex(sites)
    si = sites[ii]
    d = dim(si)
    spin_op = op(sites, ops[ii], ii)
    links[ii] = Index(1, "Link,n=$ii")
    local this_it
    if ii == 1
      this_it = ITensor(links[ii], si, si')
      for jj in 1:d, jjp in 1:d
        this_it[links[ii](1), si[jj], si'[jjp]] = spin_op[si[jj], si'[jjp]]
      end
    elseif ii == N
      this_it = ITensor(links[ii-1], si, si')
      for jj in 1:d, jjp in 1:d
        this_it[links[ii-1](1), si[jj], si'[jjp]] = spin_op[si[jj], si'[jjp]]
      end
    else
      this_it = ITensor(links[ii-1], links[ii], si, si')
      for jj in 1:d, jjp in 1:d
        this_it[links[ii-1](1), links[ii](1), si[jj], si'[jjp]] = spin_op[si[jj], si'[jjp]]
      end
    end
    its[ii] = this_it
  end
  MPO(N,its)
end

MPO(sites, ops::String) = MPO(sites, fill(ops, length(sites)))

function randomMPO(sites, m::Int=1)
  M = MPO(sites)
  for i ∈ eachindex(sites)
    randn!(M[i])
    normalize!(M[i])
  end
  m > 1 && throw(ArgumentError("randomMPO: currently only m==1 supported"))
  return M
end

length(m::MPO) = m.N_
tensors(m::MPO) = m.A_
leftlim(m::MPO) = m.llim_
rightlim(m::MPO) = m.rlim_

function set_leftlim!(m::MPO,new_ll::Int)
  m.llim_ = new_ll
end

function set_rightlim!(m::MPO,new_rl::Int)
  m.rlim_ = new_rl
end

getindex(m::MPO, n::Integer) = getindex(tensors(m), n)

function setindex!(M::MPO,T::ITensor,n::Integer)
  (n <= leftlim(M)) && set_leftlim!(M,n-1)
  (n >= rightlim(M)) && set_rightlim!(M,n+1)
  setindex!(tensors(M),T,n)
end

copy(m::MPO) = MPO(m.N_, copy(tensors(m)))
similar(m::MPO) = MPO(m.N_, similar(tensors(m)), 0, m.N_)

function deepcopy(m::T) where {T <: Union{MPO,MPS}}
    res = similar(m)
    # otherwise we will end up modifying the elements of A!
    res.A_ = deepcopy(tensors(m))
    return res
end

eachindex(m::MPO) = 1:length(m)

# TODO: optimize finding the index a little bit
# First do: scom = commonindex(A[j],x[j])
# Then do: uniqueindex(A[j],A[j-1],A[j+1],(scom,))
function siteindex(A::MPO,x::MPS,j::Integer)
  N = length(A)
  if j == 1
    si = uniqueindex(A[j],(A[j+1],x[j]))
  elseif j == N
    si = uniqueindex(A[j],(A[j-1],x[j]))
  else
    si = uniqueindex(A[j],(A[j-1],A[j+1],x[j]))
  end
  return si
end

siteinds(A::MPO,x::MPS) = [siteindex(A,x,j) for j ∈ 1:length(A)]

"""
    dag(m::MPS)
    dag(m::MPO)

Hermitian conjugation of a matrix product state or operator `m`.
"""

function dag(m::T) where {T <: Union{MPS, MPO}}
  N = length(m)
  mdag = T(N)
  @inbounds for i ∈ eachindex(m)
    mdag[i] = dag(m[i])
  end
  return mdag
end

function prime!(M::T,vargs...) where {T <: Union{MPS,MPO}}
  @inbounds for i ∈ eachindex(M)
    prime!(M[i],vargs...)
  end
end

function primelinks!(M::T, plinc::Integer = 1) where {T <: Union{MPS,MPO}}
  @inbounds for i ∈ eachindex(M)[1:end-1]
    l = linkindex(M,i)
    prime!(M[i],plinc,l)
    prime!(M[i+1],plinc,l)
  end
end

function simlinks!(M::T) where {T <: Union{MPS,MPO}}
  @inbounds for i ∈ eachindex(M)[1:end-1]
    l = linkindex(M,i)
    l̃ = sim(l)
    #M[i] *= δ(l,l̃)
    replaceindex!(M[i],l,l̃)
    #M[i+1] *= δ(l,l̃)
    replaceindex!(M[i+1],l,l̃)
  end
end

"""
maxlinkdim(M::MPS)
maxlinkdim(M::MPO)

Get the maximum link dimension of the MPS or MPO.
"""
function maxlinkdim(M::T) where {T <: Union{MPS,MPO}}
  md = 0
  for b ∈ eachindex(M)[1:end-1]
    md = max(md,dim(linkindex(M,b)))
  end
  md
end

function show(io::IO, W::MPO)
  print(io,"MPO")
  (length(W) > 0) && print(io,"\n")
  @inbounds for (i, A) ∈ enumerate(tensors(W))
    println(io,"$i  $(inds(A))")
  end
end

function linkindex(M::MPO,j::Integer)
  N = length(M)
  j ≥ length(M) && error("No link index to the right of site $j (length of MPO is $N)")
  li = commonindex(M[j],M[j+1])
  if isnothing(li)
    error("linkindex: no MPO link index at link $j")
  end
  return li
end

"""
inner(y::MPS, A::MPO, x::MPS)

Compute <y|A|x>
"""
function inner(y::MPS,
               A::MPO,
               x::MPS)::Number
  N = length(A)
  if length(y) != N || length(x) != N
      throw(DimensionMismatch("inner: mismatched lengths $N and $(length(x)) or $(length(y))"))
  end
  ydag = dag(y)
  simlinks!(ydag)
  sAx = siteinds(A,x)
  replacesites!(ydag,sAx)
  O = ydag[1]*A[1]*x[1]
  @inbounds for j ∈ eachindex(y)[2:end]
    O = O*ydag[j]*A[j]*x[j]
  end
  return O[]
end

"""
inner(B::MPO, y::MPS, A::MPO, x::MPS)

Compute <By|A|x>
"""
function inner(B::MPO,
               y::MPS,
               A::MPO,
               x::MPS)::Number
  N = length(B)
  if length(y) != N || length(x) != N || length(A) != N
      throw(DimensionMismatch("inner: mismatched lengths $N and $(length(x)) or $(length(y)) or $(length(A))"))
  end
  ydag = dag(y)
  prime!(ydag, 2)
  Bdag = dag(B)
  prime!(Bdag)
  # Swap prime levels 1 -> 2 and 2 -> 1.
  @inbounds for j ∈ eachindex(Bdag)
    swapprime!(inds(Bdag[j]),2,3)
    swapprime!(inds(Bdag[j]),1,2)
    swapprime!(inds(Bdag[j]),3,1)
  end
  O = ydag[1]*Bdag[1]*A[1]*x[1]
  @inbounds for j ∈ eachindex(y)[2:end]
    O = O*ydag[j]*Bdag[j]*A[j]*x[j]
  end
  return O[]
end

"""
errorMPOprod(y::MPS, A::MPO, x::MPS)

Compute the distance between A|x> and an approximation MPS y:
| |y> - A|x> |/| A|x> | = √(1 + (<y|y> - 2*real(<y|A|x>))/<Ax|A|x>)
"""
function errorMPOprod(y::MPS, A::MPO, x::MPS)
  N = length(A)
  if length(y) != N || length(x) != N
      throw(DimensionMismatch("inner: mismatched lengths $N and $(length(x)) or $(length(y))"))
  end
  return sqrt(abs(1. + (inner(y,y) - 2*real(inner(y,A,x)))/inner(A,x,A,x)))
end

function plussers(::Type{T}, left_ind::Index, right_ind::Index, sum_ind::Index) where {T<:Array}
    #if dir(left_ind) == dir(right_ind) == Neither
        total_dim    = dim(left_ind) + dim(right_ind)
        total_dim    = max(total_dim, 1)
        # TODO: I am not sure if we should be using delta
        # tensors for this purpose? I think we should consider
        # not allowing them to be made with different index sizes
        #left_tensor  = δ(left_ind, sum_ind)
        left_tensor  = diagITensor(1.0,left_ind, sum_ind)
        right_tensor = ITensor(right_ind, sum_ind)
        for i in 1:dim(right_ind)
            right_tensor[right_ind(i), sum_ind(dim(left_ind) + i)] = 1
        end
        return left_tensor, right_tensor
    #else # tensors have QNs
    #    throw(ArgumentError("support for adding MPOs with defined quantum numbers not implemented yet."))
    #end
end

function sum(A::T, B::T; kwargs...) where {T <: Union{MPS, MPO}}
    n = A.N_
    length(B) =! n && throw(DimensionMismatch("lengths of MPOs A ($n) and B ($(length(B))) do not match"))
    orthogonalize!(A, 1; kwargs...)
    orthogonalize!(B, 1; kwargs...)
    C = similar(A)
    rand_plev = 13124
    lAs = [linkindex(A, i) for i in 1:n-1]
    prime!(A, rand_plev, "Link")

    first  = Vector{ITensor{2}}(undef,n-1)
    second = Vector{ITensor{2}}(undef,n-1)
    for i in 1:n-1
        lA = linkindex(A, i)
        lB = linkindex(B, i)
        r  = Index(dim(lA) + dim(lB), tags(lA))
        f, s = plussers(typeof(data(store(A[1]))), lA, lB, r)
        first[i]  = f
        second[i] = s
    end
    C[1] = A[1] * first[1] + B[1] * second[1]
    for i in 2:n-1
        C[i] = dag(first[i-1]) * A[i] * first[i] + dag(second[i-1]) * B[i] * second[i]
    end
    C[n] = dag(first[n-1]) * A[n] + dag(second[n-1]) * B[n]
    prime!(C, -rand_plev, "Link")
    truncate!(C; kwargs...)
    return C
end

function applyMPO(A::MPO, psi::MPS; kwargs...)::MPS
  method = get(kwargs, :method, "DensityMatrix")
  if method == "DensityMatrix" || method == "densitymatrix"
    return densityMatrixApplyMPO(A, psi; kwargs...)
  elseif method == "naive" || method == "Naive"
    return naiveApplyMPO(A, psi; kwargs...)
  end
  throw(ArgumentError("Method $method not supported"))
end

function densityMatrixApplyMPO(A::MPO, psi::MPS; kwargs...)::MPS
    n = length(A)
    n != length(psi) && throw(DimensionMismatch("lengths of MPO ($n) and MPS ($(length(psi))) do not match"))
    psi_out         = similar(psi)
    cutoff::Float64 = get(kwargs, :cutoff, 1e-13)

    maxdim::Int     = get(kwargs,:maxdim,maxlinkdim(psi))
    mindim::Int     = max(get(kwargs,:mindim,1), 1)
    normalize::Bool = get(kwargs, :normalize, false)
    all(x -> x != Index(), [siteindex(A, psi, j) for j in 1:n]) || throw(ErrorException("MPS and MPO have different site indices in applyMPO method 'DensityMatrix'"))

    rand_plev = 14741
    psi_c     = dag(copy(psi))
    A_c       = dag(copy(A))
    prime!(psi_c, rand_plev)
    prime!(A_c, rand_plev)
    for j in 1:n-1
        unique_site_ind = setdiff(findinds(A_c[j], "Site"), findindex(psi_c[j], "Site"))[1]
        pl = id(unique_site_ind) == id(commonindex(A_c[j], psi_c[j])) ? 1 : 0
        A_c[j] = setprime(A_c[j], pl, unique_site_ind)
    end
    E = Vector{ITensor}(undef, n-1)
    E[1] = psi[1]*A[1]*A_c[1]*psi_c[1]
    for j in 2:n-1
        E[j] = E[j-1]*psi[j]*A[j]*A_c[j]*psi_c[j]
    end
    O     = psi[n] * A[n]
    ρ     = E[n-1] * O * dag(prime(O, rand_plev))
    ts    = tags(commonindex(psi[n], psi[n-1]))
    Lis   = commonindex(ρ, A[n])
    Ris   = uniqueinds(ρ, Lis)
    FU, D = eigenHermitian(ρ, Lis, Ris; ispossemidef=true,
                                        tags=ts,
                                        kwargs...)
    psi_out[n] = setprime(dag(FU), 0, "Site")
    O     = O * FU * psi[n-1] * A[n-1]
    O     = prime(O, -1, "Site")
    for j in reverse(2:n-1)
        dO  = prime(dag(O), rand_plev)
        ρ   = E[j-1] * O * dO
        ts  = tags(commonindex(psi[j], psi[j-1]))
        Lis = IndexSet(commonindex(ρ, A[j]), commonindex(ρ, psi_out[j+1]))
        Ris = uniqueinds(ρ, Lis)
        FU, D = eigenHermitian(ρ, Lis, Ris; ispossemidef=true,
                                            tags=ts,
                                            kwargs...)
        psi_out[j] = dag(FU)
        O = O * FU * psi[j-1] * A[j-1]
        O = prime(O, -1, "Site")
    end
    if normalize
        O /= norm(O)
    end
    psi_out[1]    = copy(O)
    psi_out.llim_ = 0
    psi_out.rlim_ = 2
    return psi_out
end

function naiveApplyMPO(A::MPO, psi::MPS; kwargs...)::MPS
  N = length(A)
  if N != length(psi)
    throw(DimensionMismatch("lengths of MPO ($N) and MPS ($(length(psi))) do not match"))
  end

  psi_out = MPS(N)
  for j=1:N
    psi_out[j] = noprime(A[j]*psi[j])
  end

  for b=1:(N-1)
    Al = commonindex(A[b],A[b+1])
    pl = commonindex(psi[b],psi[b+1])
    C,_ = combiner(Al,pl)
    psi_out[b] *= C
    psi_out[b+1] *= C
  end

  truncate!(psi_out;kwargs...)

  return psi_out
end

function multMPO(A::MPO, B::MPO; kwargs...)::MPO
    cutoff::Float64 = get(kwargs, :cutoff, 1e-14)
    resp_degen::Bool = get(kwargs, :respect_degenerate, true)
    maxdim::Int = get(kwargs,:maxdim,maxlinkdim(A)*maxlinkdim(B))
    mindim::Int = max(get(kwargs,:mindim,1), 1)
    N = length(A)
    N != length(B) && throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
    A_ = copy(A)
    orthogonalize!(A_, 1)
    B_ = copy(B)
    orthogonalize!(B_, 1)

    links_A = findinds.(A.A_, "Link")
    links_B = findinds.(B.A_, "Link")

    for i in 1:N
        if length(commoninds(findinds(A_[i], "Site"), findinds(B_[i], "Site"))) == 2
            A_[i] = prime(A_[i], "Site")
        end
    end
    res = deepcopy(A_)
    for i in 1:N-1
        ci = commonindex(res[i], res[i+1])
        new_ci = Index(dim(ci), tags(ci))
        replaceindex!(res[i], ci, new_ci)
        replaceindex!(res[i+1], ci, new_ci)
        @assert commonindex(res[i], res[i+1]) != commonindex(A[i], A[i+1])
    end
    sites_A = Index[]
    sites_B = Index[]
    @inbounds for (AA, BB) in zip(tensors(A_), tensors(B_))
        sda = setdiff(findinds(AA, "Site"), findinds(BB, "Site"))
        sdb = setdiff(findinds(BB, "Site"), findinds(AA, "Site"))
        sda_ind = setprime(sda[1], 0) == sdb[1] ? plev(sda[1]) == 1 ? sda[1] : setprime(sda[1], 1) : setprime(sda[1], 0)
        push!(sites_A, sda_ind)
        push!(sites_B, sdb[1])
    end
    res[1] = ITensor(sites_A[1], sites_B[1], commonindex(res[1], res[2]))
    for i in 1:N-2
        if i == 1
            clust = A_[i] * B_[i]
        else
            clust = nfork * A_[i] * B_[i]
        end
        lA = commonindex(A_[i], A_[i+1])
        lB = commonindex(B_[i], B_[i+1])
        nfork = ITensor(lA, lB, commonindex(res[i], res[i+1]))
        res[i], nfork = factorize(mapprime(clust,2,1), inds(res[i]), dir="fromleft", tags=tags(lA), cutoff=cutoff, maxdim=maxdim, mindim=mindim)
        mid = dag(commonindex(res[i], nfork))
        res[i+1] = ITensor(mid, sites_A[i+1], sites_B[i+1], commonindex(res[i+1], res[i+2]))
    end
    clust = nfork * A_[N-1] * B_[N-1]
    nfork = clust * A_[N] * B_[N]

    # in case we primed A
    A_ind = uniqueindex(findinds(A_[N-1], "Site"), findinds(B_[N-1], "Site"))
    Lis = IndexSet(A_ind, sites_B[N-1], commonindex(res[N-2], res[N-1]))
    U, V = factorize(nfork,Lis,dir="fromright",cutoff=cutoff,which_factorization="svd",tags="Link,n=$(N-1)",maxdim=maxdim,mindim=mindim)
    res[N-1] = U
    res[N] = V
    truncate!(res;kwargs...)
    for i in 1:N
        res[i] = mapprime(res[i], 2, 1)
    end
    return res
end

function orthogonalize!(M::Union{MPS,MPO},
                        j::Int;
                        kwargs...)
  while leftlim(M) < (j-1)
    (leftlim(M) < 0) && set_leftlim!(M,0)
    b = leftlim(M)+1
    linds = uniqueinds(M[b],M[b+1])
    Q,R = qr(M[b], linds)
    M[b] = Q
    M[b+1] *= R
    set_leftlim!(M,b)
    if rightlim(M) < leftlim(M)+2
      set_rightlim!(M,leftlim(M)+2)
    end
  end

  N = length(M)

  while rightlim(M) > (j+1)
    (rightlim(M) > (N+1)) && set_rightlim!(M,N+1)
    b = rightlim(M)-2
    rinds = uniqueinds(M[b+1],M[b])
    Q,R = qr(M[b+1], rinds)
    M[b+1] = Q
    M[b] *= R
    set_rightlim!(M,b+1)
    if leftlim(M) > rightlim(M)-2
      set_leftlim!(M,rightlim(M)-2)
    end
  end
end

function truncate!(M::Union{MPS,MPO}; kwargs...)

  N = length(M)

  # Left-orthogonalize all tensors to make
  # truncations controlled
  orthogonalize!(M,N)

  # Perform truncations in a right-to-left sweep
  for j in reverse(2:N)
    rinds = uniqueinds(M[j],M[j-1])
    U,S,V = svd(M[j],rinds;kwargs...)
    M[j] = U
    M[j-1] *= (S*V)
    set_rightlim!(M,j)
  end

end

@doc """
orthogonalize!(M::MPS, j::Int; kwargs...)

Move the orthogonality center of the MPS
to site j. No observable property of the
MPS will be changed, and no truncation of the
bond indices is performed. Afterward, tensors
1,2,...,j-1 will be left-orthogonal and tensors
j+1,j+2,...,N will be right-orthogonal.

orthogonalize!(W::MPO, j::Int; kwargs...)

Move the orthogonality center of an MPO to site j.
""" orthogonalize!

@doc """
truncate!(M::MPS; kwargs...)

Perform a truncation of all bonds of an MPS,
using the truncation parameters (cutoff,maxdim, etc.)
provided as keyword arguments.

truncate!(M::MPO; kwargs...)

Perform a truncation of all bonds of an MPO,
using the truncation parameters (cutoff,maxdim, etc.)
provided as keyword arguments.
""" truncate!


text:
"""
    sankey(names, source, target, value)

Creates an `EChart` diagram displaying the path and frequencies between nodes.

## Methods
```julia
sankey(names::AbstractVector, source::AbstractVector{<:Union{Missing, Int}},
    target::AbstractVector{<:Union{Missing, Int}}, value::AbstractVector{<:Union{Missing, Real}};
```

## Arguments
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

## Examples
```julia
name = ["Agricultural waste", "Bio-conversion", "Liquid", "Losses", "Solid", "Gas", "Biofuel imports",
    "Biomass imports", "Coal imports", "Coal"]
source = [0, 1, 1, 1, 1, 6, 7, 8, 10]
target = [1, 2, 3, 4, 5, 2, 4, 9, 9]
value = [124.729, 0.597, 26.862, 280.322, 81.144, 35, 35, 11.606, 63.965]
s = sankey(name, source, target, value)
```
"""
function sankey(names::AbstractVector,
                source::AbstractVector{<:Union{Missing, Int}},
                target::AbstractVector{<:Union{Missing, Int}},
                value::AbstractVector{<:Union{Missing, Real}};
                kwargs...)

    ec = newplot(kwargs, ec_charttype = "sankey")
    ec.series = [SankeySeries(name = "Series 1", layout = "none", data = arrayofdicts(name = names),
        links = arrayofdicts(source = source, target = target, value = value))]

    return ec

end

ext: jl


text:
module MyModel

using Mimi

include("region_parameters.jl")
include("gross_economy.jl")
include("emissions.jl")

export construct_MyModel

function construct_MyModel()

	m = Model()

	set_dimension!(m, :time, collect(2015:5:2110))
	set_dimension!(m, :regions, [:Region1, :Region2, :Region3])	 # Note that the regions of your model must be specified here

	add_comp!(m, grosseconomy)
	add_comp!(m, emissions)

	set_param!(m, :grosseconomy, :l, l)
	set_param!(m, :grosseconomy, :tfp, tfp)
	set_param!(m, :grosseconomy, :s, s)
	set_param!(m, :grosseconomy, :depk,depk)
	set_param!(m, :grosseconomy, :k0, k0)
	set_param!(m, :grosseconomy, :share, 0.3)

	# set parameters for emissions component
	set_param!(m, :emissions, :sigma, sigma)
	connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    return m

end

end #module


text:
module PointGroup

using Molecules

include("sea.jl")
include("moit.jl")
end


text:
using EquationsOfState
using Test

@testset "EquationsOfState.jl" begin
    include("Collections.jl")
    include("FiniteStrains.jl")
    include("NonlinearFitting.jl")
    include("LinearFitting.jl")
    # include("Find.jl")
end


text:
uT = u"yr"
uL = u"m"
uM = u"g"
uQ = uM*uL^(-2)  # unit of pools, i.e. quantities: g/m^2

ustrip_num(num::Num) = Unitful.ustrip(ModelingToolkit.get_unit(num), ModelingToolkit.value(num))

#ustrip(::)
# Unitful.ustrip(ModelingToolkit.get_unit(k),v)


text:
module Util
export detailindex, detailrange, scalingrange, detailn, nscales, maxlevel, tl2level, level2tl,
    mirror, upsample, downsample, iscube, isdyadic, wcount, circshift!,
    split!, merge!, stridedcopy!,
    isvalidtree, maketree,
    makewavelet, testfunction

# WAVELET INDEXING AND SIZES

# detail coef at level j location i (i starting at 1) -> vector index
detailindex(j::Integer,i::Integer) = 2^j+i
# the range of detail coefs at level j
detailrange(j::Integer) = (2^j+1):(2^(j+1))
# the range of scaling coefs at level j
scalingrange(j::Integer) = 1:2^j
# number of detail coefs at level j
detailn(j::Integer) = 2^j
# number of scales of dyadic length signal (n=2^J)
nscales(n::Integer) = int(log2(n))
nscales(x::Vector) = nscales(length(x))
nscales(x::Array) = nscales(size(x,1))
# the largest detail level
maxlevel(n::Integer) = nscales(n)-1
maxlevel(x::Vector) = maxlevel(length(x))
# convert number of transformed levels L to minimum level j
tl2level(n::Integer, L::Integer) = nscales(n)-L
tl2level(x::Vector, L::Integer) = tl2level(length(x), L)
# convert maximum level j to number of transformed levels L
level2tl(arg...) = tl2level(arg...)

# UTILITY FUNCTIONS

# mirror of filter
mirror{T<:Number}(f::AbstractVector{T}) = f.*(-1).^(0:length(f)-1)
# upsample
function upsample(x::AbstractVector, sw::Int=0)
    @assert sw==0 || sw==1
    n = length(x)
    y = zeros(eltype(x), n<<1)
    sw -= 1

    for i = 1:n
        @inbounds y[i<<1 + sw] = x[i]
    end
    return y
end
# downsample
function downsample(x::AbstractVector, sw::Int=0)
    @assert sw==0 || sw==1
    n = length(x)
    @assert n%2 == 0
    y = zeros(eltype(x), n>>1)
    sw -= 1

    for i = 1:length(y)
        @inbounds y[i] = x[i<<1 + sw]
    end
    return y
end

# are all dimensions equal length?
function iscube(x::AbstractArray)
    for i = 1:ndims(x)
        size(x,1)!=size(x,i) && return false
    end
    return true
end
# are all dimensions dyadic length?
function isdyadic(x::AbstractArray)
    for i = 1:ndims(x)
        n = size(x,i)
        J = nscales(n)
        n != 2^J && return false
    end
    return true
end

# count coefficients above threshold t (>=), excluding coefficients in levels < level
# where level -1 is the x[1] coefficient
function wcount(x::AbstractVector, t::Real=0; level::Int=-1)
    @assert level >= -1
    c = 0
    si = 1
    level >= 0 && (si = 1 + 2^level)
    @inbounds for i = si:length(x)
        if abs(x[i]) >= t
            c += 1
        end
    end
    return c
end
# count coefficients above threshold t (>=)
function wcount(x::AbstractArray, t::Real=0)
    c = 0
    @inbounds for i = 1:length(x)
        if abs(x[i]) >= t
            c += 1
        end
    end
    return c
end

# inplace circular shift of vector by shift, such that anew[i]=aold[i-shift] (mod length(a))
function circshift!(a::AbstractVector, shift::Integer)
    atype = typeof(a)
    s = length(a)
    shift = mod(shift,s)
    shift == 0 && return a
    shift = ifelse(s>>1 < shift, shift-s, shift)  # shift a smaller distance if possible
    if shift < 0
        tmp = a[1:-shift]
        for i = 1:s+shift
            @inbounds a[i] = a[i-shift]
        end
        a[s+1+shift:s] = tmp
    else
        tmp = a[s+1-shift:s]
        for i = s:-1:1+shift
            @inbounds a[i] = a[i-shift]
        end
        a[1:shift] = tmp
    end
    return a::atype
end
# out of place circular shift of vector by shift, such that b[i]=a[i-shift] (mod length(a))
function circshift!(b::AbstractVector, a::AbstractVector, shift::Integer)
    @assert length(a) == length(b)
    atype = typeof(a)
    s = length(a)
    shift = mod(shift,s)
    shift == 0 && return copy!(b,a)
    shift = ifelse(s>>1 < shift, shift-s, shift)  # shift a smaller distance if possible
    if shift < 0
        for i = 1:s+shift
            @inbounds b[i] = a[i-shift]
        end
        sh = s+shift
        for i = s+shift+1:s
            @inbounds b[i] = a[i-sh]
        end
    else
        for i = s:-1:1+shift
            @inbounds b[i] = a[i-shift]
        end
        sh = -s+shift
        for i = shift:-1:1
            @inbounds b[i] = a[i-sh]
        end
    end
    return b::atype
end

# put odd elements into first half, even into second half
function split!{T<:Number}(a::AbstractVector{T})
    n = length(a)
    nt = n>>2 + (n>>1)%2
    tmp = Array(T, nt)
    split!(a, n, tmp)
    return a
end

# split only the range 1:n
function split!{T<:Number}(a::AbstractVector{T}, n::Integer, tmp::Vector{T})
    @assert n <= length(a)
    @assert n%2 == 0
    n == 2 && return a
    nt = n>>2 + (n>>1)%2
    @assert nt <= length(tmp)

    for i=1:nt # store evens
        @inbounds tmp[i] = a[i<<1]
    end
    for i=1:n>>1  # odds to first part
        @inbounds a[i] = a[(i-1)<<1 + 1]
    end
    for i=0:nt - 1  # evens to end
        @inbounds a[n-i] = a[n - i<<1]
    end
    copy!(a,n>>1+1,tmp,1,nt)
    return a
end

# out of place split from a to b, only the range 1:n
function split!{T<:Number}(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer)
    @assert n <= length(a) && n <= length(b)
    @assert n%2 == 0
    if n == 2
        b[1] = a[1]
        b[2] = a[2]
        return b
    end

    h = n>>1
    for i=1:h  # odds to b
        @inbounds b[i] = a[(i-1)<<1 + 1]
    end
    for i=h+1:n  # evens to b
        @inbounds b[i] = a[(i - h)<<1]
    end

    return b
end
# out of place split from a to b, only the range a[ia:inca:ia+(n-1)*inca] to b[1:n]
function split!{T<:Number}(b::AbstractVector{T}, a::AbstractArray{T}, ia::Integer, inca::Integer, n::Integer)
    @assert ia+(n-1)*inca <= length(a) && n <= length(b)
    @assert n%2 == 0
    if n == 2
        b[1] = a[ia]
        b[2] = a[ia + inca]
        return b
    end

    h = n>>1
    inca2 = inca<<1
    for i=1:h  # odds to b
        @inbounds b[i] = a[ia + (i-1)*inca2]
    end
    iainca = ia + inca
    hp1 = h + 1
    for i=h+1:n  # evens to b
        @inbounds b[i] = a[iainca + (i - hp1)*inca2]
    end

    return b
end

# inverse the operation of split!
function merge!{T<:Number}(a::AbstractVector{T})
    n = length(a)
    nt = n>>2 + (n>>1)%2
    tmp = Array(T,nt)
    merge!(a, n, tmp)
    return a
end

# merge only the range 1:n
function merge!{T<:Number}(a::AbstractVector{T}, n::Integer, tmp::Vector{T})
    @assert n <= length(a)
    @assert n%2 == 0
    n == 2 && return a
    nt = n>>2 + (n>>1)%2
    @assert nt <= length(tmp)

    copy!(tmp,1,a,n>>1+1,nt)
    for i=nt-1:-1:0  # evens from end
        @inbounds a[n - 2*i] = a[n-i]
    end
    for i=n>>1:-1:1  # odds from first part
        @inbounds a[(i-1)<<1 + 1] = a[i]
    end
    for i=nt:-1:1 # retrieve evens
        @inbounds a[i<<1] = tmp[i]
    end
    return a
end

# out of place merge from a to b, only the range 1:n
function merge!{T<:Number}(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer)
    @assert n <= length(a) && n <= length(b)
    @assert n%2 == 0
    if n == 2
        b[1] = a[1]
        b[2] = a[2]
        return b
    end

    h = n>>1
    for i=1:h  # odds to b
        @inbounds b[(i-1)<<1 + 1] = a[i]
    end
    for i=h+1:n  # evens to b
        @inbounds b[(i - h)<<1] = a[i]
    end

    return b
end
# out of place merge from a to b, only the range a[1:n] to b[ib:incb:ib+(n-1)*incb]
function merge!{T<:Number}(b::AbstractArray{T}, ib::Integer, incb::Integer, a::AbstractVector{T}, n::Integer)
    @assert n <= length(a) && ib+(n-1)*incb <= length(b)
    @assert n%2 == 0
    if n == 2
        b[ib] = a[1]
        b[ib + incb] = a[2]
        return b
    end

    h = n>>1
    incb2 = incb<<1
    for i=1:h  # odds to b
        @inbounds b[ib + (i-1)*incb2] = a[i]
    end
    ibincb = ib + incb
    hp1 = h + 1
    for i=h+1:n  # evens to b
        @inbounds b[ibincb + (i - hp1)*incb2] = a[i]
    end

    return b
end

function stridedcopy!{T<:Number}(b::AbstractVector{T}, a::AbstractArray{T}, ia::Integer, inca::Integer, n::Integer)
    @assert ia+(n-1)*inca <= length(a) && n <= length(b)

    @inbounds for i = 1:n
        b[i] = a[ia + (i-1)*inca]
    end
    return b
end
function stridedcopy!{T<:Number}(b::AbstractArray{T}, ib::Integer, incb::Integer, a::AbstractVector{T}, n::Integer)
    @assert ib+(n-1)*incb <= length(b) && n <= length(a)

    @inbounds for i = 1:n
        b[ib + (i-1)*incb] = a[i]
    end
    return b
end

# wavelet packet transforms WPT
# valid if 0 nodes have 0 children and length is correct
function isvalidtree(x::AbstractVector, b::BitVector)
    ns = nscales(x)
    nb = length(b)
    length(b) == 2^(ns)-1 || return false
    @assert (2^(ns-1)-1)<<1+1 <= nb

    for i in 1:2^(ns-1)-1
        @inbounds if !b[i] && (b[i<<1] || b[i<<1+1])
            return false
        end
    end
    return true
end
# return a tree (BitVector)
# s=:full, all nodes for first L levels equal 1, others 0
# s=:dwt, nodes corresponding to a dwt for first L levels equal 1, others 0
function maketree(n::Int, L::Int, s::Symbol=:full)
    ns = nscales(n)
    @assert 0 <= L <= ns
    nb = 2^(ns)-1
    b = BitArray(nb)
    fill!(b, false)
    @assert (2^(ns-1)-1)<<1+1 <= nb

    t = true
    if s == :full
        for i in 1:2^(L)-1
            @inbounds b[i] = t
        end
    elseif s == :dwt
        for i in 1:L
            @inbounds b[2^(i-1)] = t
        end
    else
        error("uknown symbol")
    end
    return b
end

# return scaling and wavelet functions and location vector, made from filter h
# iterated with a cascade algorithm with N steps
makewavelet(h, arg...) = makewavelet(h.qmf, arg...)
function makewavelet(h::AbstractVector, N::Integer=8)
	@assert N>=0
	sc = norm(h)
	h = h*sqrt(2)/sc
	phi = copy(h)
	psi = mirror(reverse(h))

	for i=1:N
		phi = conv(upsample(phi), h)
		psi = conv(upsample(psi), h)
	end
	phi = phi[1:end-2^(N)+1]
	psi = psi[1:end-2^(N)+1]
	return scale!(phi,sc/sqrt(2)), scale!(psi,sc/sqrt(2)), linspace(0,length(h)-1,length(psi))
end

# return a vector of test function values on [0,1), see
#Donoho, D.L.; I.M. Johnstone (1994), "Ideal spatial adaptation by wavelet shrinkage," Biometrika, vol. 81, pp. 425–455.
function testfunction(n::Int, ft::String)
    @assert n >= 1
    f = Array(Float64,n)
    range = 0:1/n:1-eps()
    i = 1
    if ft=="Blocks"
        tj = [0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81]
        hj = [4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,2.1,-4.2]
        for t in range
            f[i] = 0
            for j = 1:11
                f[i] += hj[j]*(1+sign(t-tj[j]))/2
            end
            i += 1
        end
    elseif ft=="Bumps"
        tj = [0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81]
        hj = [4,5,3,4,5,4.2,2.1,4.3,3.1,5.1,4.2]
        wj = [0.005,0.005,0.006,0.01,0.01,0.03,0.01,0.01,0.005,0.008,0.005]
        for t in range
            f[i] = 0
            for j = 1:11
                f[i] += hj[j]/(1+abs((t-tj[j])/wj[j]))^4
            end
            i += 1
        end
    elseif ft=="HeaviSine"
        for t in range
            f[i] = 4*sin(4*pi*t) - sign(t-0.3) - sign(0.72-t)
            i += 1
        end
    elseif ft=="Doppler"
        for t in range
            f[i] = sqrt(t*(1-t))*sin(2*pi*1.05/(t+0.05))
            i += 1
        end
    else
        error("test function not found")
    end
    return f
end

function copygeneral2!(y::AbstractArray, x::AbstractMatrix, t::Range, d::Integer)
    k = 1
    for i in t
        @inbounds y[k] = x[i,d]
        k +=1
    end
    return y
end
function copygeneral2!(y::AbstractArray, x::AbstractMatrix, d::Integer, t::Range)
    k = 1
    for i in t
        @inbounds y[k] = x[d,i]
        k +=1
    end
    return y
end
function copygeneral1!(y::AbstractMatrix, t::Range, d::Integer, x::AbstractArray)
    k = 1
    for i in t
        @inbounds y[i,d] = x[k]
        k +=1
    end
    return y
end
function copygeneral1!(y::AbstractMatrix, d::Integer, t::Range, x::AbstractArray)
    k = 1
    for i in t
        @inbounds y[d,i] = x[k]
        k +=1
    end
    return y
end

end # module


text:
#include("../src/ttv.jl")
#include("/Users/ericagol/Computer/Julia/regress.jl")

#maxabs(x) = maximum(abs.(x))

#using Base.Test

#KEPLER_TOL  = 1e-15
#TRANSIT_TOL = 1e-15

#@testset "findtransit3" begin

#n = 8
n = 3
t0 = 7257.93115525
#h  = 0.12
h  = 0.05
#tmax = 600.0
#tmax = 100.0
tmax = 10.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')
# Increase masses for debugging:
elements[2,1] *= 10.0
elements[3,1] *= 10.0

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
ntt = zeros(Int64,n)
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
tt  = zeros(n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
# Now, compute derivatives (with respect to initial cartesian positions/masses at
# beginning of each step):
dtdq0 = zeros(n,maximum(ntt),7,n)
dtdq0_num = zeros(BigFloat,n,maximum(ntt),7,n)
dlnq = big(1e-10)
# Make radius of star large:
rstar = 1e12
dtdelements_num = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,dtdq0_num,dlnq,rstar)
# Now do computation in BigFloat precision:
dtdq0_big = big.(dtdq0)
dtdelements_num = ttv_elements!(n,big(t0),big(h),big(tmax),big.(elements),big.(tt),count,dtdq0_big,big(rstar))

dtdq0_num = convert(Array{Float64,4},dtdq0_num)
dtdq0_big = convert(Array{Float64,4},dtdq0_big)

mask = zeros(Bool, size(dtdq0))
for i=2:n, j=1:count[i], k=1:5, l=1:n
  mask[i,j,k,l] = true
end
#println("Max diff log(dtdq0): ",maximum(abs.(dtdq0_num[mask]./dtdq0[mask]-1.0)))
println("Max diff asinh(dtdq0): ",maximum(abs.(asinh.(dtdq0_num[mask])-asinh.(dtdq0[mask]))))
println("Max diff     dtdq0 : ",maximum((dtdq0_num[mask]-dtdq0[mask])))
#@test isapprox(dtdq0[mask],convert(Array{Float64,4},dtdq0_num)[mask];norm=maxabs)
@test isapprox(asinh.(dtdq0[mask]),asinh.(convert(Array{Float64,4},dtdq0_num)[mask]);norm=maxabs)
#unit = ones(dtdq0[mask])
#@test isapprox(dtdq0[mask]./convert(Array{Float64,4},dtdq0_num)[mask],unit;norm=maxabs)
#end


text:
################
# Bounding boxes
#################

# A bounding box is an Interval or ProductDomain of intervals that encompasses the domain.

# If the boundingbox is not a product of intervals, something has gone wrong.

boundingbox(a::SVector{1}, b::SVector{1}) = a[1]..b[1]

boundingbox(a::Number, b::Number) = a..b

boundingbox(a, b) = ProductDomain(map((ai,bi)->ClosedInterval(ai,bi), a, b)...)

boundingbox(d::AbstractInterval) = d

boundingbox(::UnitHyperBall{N,T}) where {N,T} = boundingbox(-ones(SVector{N,T}), ones(SVector{N,T}))

boundingbox(d::ProductDomain) = cartesianproduct(map(boundingbox, elements(d))...)

boundingbox(d::DerivedDomain) = boundingbox(source(d))

boundingbox(d::DifferenceDomain) = boundingbox(d.d1)

function boundingbox(d::UnionDomain)
    left = SVector(minimum(hcat(map(infimum,map(boundingbox,elements(d)))...);dims=2)...)
    right = SVector(maximum(hcat(map(supremum,map(boundingbox,elements(d)))...);dims=2)...)
    boundingbox(left,right)
end

function boundingbox(d::IntersectionDomain)
    left = SVector(maximum(hcat(map(infimum,map(boundingbox,elements(d)))...);dims=2)...)
    right = SVector(minimum(hcat(map(supremum,map(boundingbox,elements(d)))...);dims=2)...)
    boundingbox(left,right)
end

DomainSets.superdomain(d::DomainSets.MappedDomain) = DomainSets.source(d)

# Now here is a problem: how do we compute a bounding box, without extra knowledge
# of the map? We can only do this for some maps.
boundingbox(d::DomainSets.MappedDomain) = mapped_boundingbox(boundingbox(source(d)), forward_map(d))

function mapped_boundingbox(box::Interval, fmap)
    l,r = (infimum(box),supremum(box))
    ml = fmap*l
    mr = fmap*r
    boundingbox(min(ml,mr), max(ml,mr))
end

# In general, we can at least map all the corners of the bounding box of the
# underlying domain, and compute a bounding box for those points. This will be
# correct for affine maps.
function mapped_boundingbox(box::ProductDomain, fmap)
    crn = corners(infimum(box),supremum(box))
    mapped_corners = [fmap*crn[:,i] for i in 1:size(crn,2)]
    left = [minimum([mapped_corners[i][j] for i in 1:length(mapped_corners)]) for j in 1:size(crn,1)]
    right = [maximum([mapped_corners[i][j] for i in 1:length(mapped_corners)]) for j in 1:size(crn,1)]
    boundingbox(left, right)
end

# Auxiliary functions to rotate a bounding box when mapping it.
function corners(left::AbstractVector, right::AbstractVector)
    @assert length(left)==length(right)
    N=length(left)
    corners = zeros(N,2^N)
    # All possible permutations of the corners
    for i=1:2^length(left)
        for j=1:N
            corners[j,i] = ((i>>(j-1))%2==0) ? left[j] : right[j]
        end
    end
    corners
end


text:
using Distributed
#w=WorkerPool([2, 3])
procs=3

addprocs(procs) # total: procs + main proc
println("Processes: ",nprocs())

import Base.show

@everywhere mutable struct Worker
  alive::Bool
  orderID::Int
  tick::UInt
  tasks::Channel{Task} #tasks.sz_max
  Worker() = new(true,0,0,Channel{Task}(0))
end

@everywhere function worker_stop() WORKER.alive=false end
@everywhere function worker_alive() WORKER.alive end
@everywhere function worker_tick(t=0); if t>0 WORKER.tick=t; end; WORKER.tick; end
@everywhere worker_ids() = (myid(),WORKER.orderID)
@everywhere worker_task(f::Function) = push!(WORKER.tasks, Task(f))

@everywhere function Base.show(io::IO, this::Worker)
  print(io, "#", this.orderID, " ")
end

@everywhere WORKER = Worker()
#@everywhere JOBS = Channel{Task}(0)

@everywhere function write(i)
  t=rand()
  sleep(t) #Libc.systemsleep(t)
  println((i,myid(),t))
  (i,myid())
end

@everywhere NOTASK=()->nothing

@everywhere function worker(update::Function; start=NOTASK,close=NOTASK,tick=0)
  println(WORKER, "start")
  start()
  running=true
  while running == true running=update() end
  println(WORKER, "close")
  close()
  (worker_ids(),"end")
end

MESSAGES=RemoteChannel(()->Channel{String}(0))
#RemoteChannel(()->Channel(messages))

Printer = ()->begin
  println(take!(MESSAGES))
  sleep(1)
  true
end

Writer = ()->begin
  put!(MESSAGES,"test")
  println("writes")
  sleep(1)
  true
end

WORKERS = Future[]
wp=WorkerPool([2, 3])

print("START")
push!(WORKERS, remotecall(worker, wp, Printer)); print(".");
push!(WORKERS, remotecall(worker, wp, Writer)); print(".");
#for i=1:length(wp.workers) push!(WORKERS, remotecall(()->worker(i), wp)); print("."); end
println()
#=
#remotecall(worker_task(()->println("hi")), wp)
#sleep(3)

println("CLOSE") # start tasks on the workers to process requests in parallel
for w in wp.workers remotecall(worker_stop, wp) end
sleep(3)
=#
println("READ")
for w in WORKERS println(fetch(w)) end
println("END")

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module Day9

export get_inputs, get_solution1, get_solution2

## Input getting
read_input(f_name) = parse.(Int, readlines(joinpath(@__DIR__, f_name)))

function get_inputs()
    test_input1 = test_input2 = (read_input("test_input1.txt"), 5)
    test_output1 = 127
    test_output2 = 62
    data = (read_input("input.txt"), 25)
    return (; test_input1, test_output1, test_input2, test_output2, data)
end

## Solution functions
function is_valid(number, preamble)
    for i in eachindex(preamble), j in i:length(preamble)
        if preamble[i] + preamble[j] == number return true end
    end
    return false
end

# Part 1
function get_solution1((data, preamble_length))
    for i in eachindex(data)
        j = i + preamble_length
        preamble = @view data[i:j-1]
        val = data[j]
        if !is_valid(val, preamble) return val end
    end
    return -1
end

# Part 2
function get_solution2((data, preamble_length))
    sol1 = get_solution1((data, preamble_length))
    for i in eachindex(data)
        val = 0
        for j in i:length(data)
            val += data[j]
            if val==sol1
                this_range = @view data[i:j]
                return sum(extrema(this_range))
            elseif val > sol1
                break
            end
        end
    end
    return -1
end

end


text:
"""
```
sum!(R1::Result, R2::Result)
```

Sum the counts of two Results structures, adding the result to the first structure
as in R1 = R1 + R2.

"""
function sum!(R1::Result, R2::Result)

    @. R1.md_count += R2.md_count
    @. R1.md_count_random += R2.md_count_random

    @. R1.solute_atom += R2.solute_atom
    @. R1.solvent_atom += R2.solvent_atom

    @. R1.rdf_count += R2.rdf_count
    @. R1.rdf_count_random += R2.rdf_count_random

    sum!(R1.density, R2.density)
    sum!(R1.volume, R2.volume)

    return nothing
end

function sum!(D1::Density, D2::Density)
    D1.solute += D2.solute
    D1.solvent += D2.solvent
    D1.solvent_bulk += D2.solvent_bulk
    return nothing
end

function sum!(V1::Volume, V2::Volume)
    V1.total += V2.total
    V1.bulk += V2.bulk
    V1.domain += V2.domain
    @. V1.shell += V2.shell
    return nothing
end


text:
#!/usr/bin/env julia
using SkyCoords
using BenchmarkTools

function myprintln(io, string)
    println(io, string)
    println(string)
end

io = open("julia_times.csv", "w")
println(io, "n,system,time")
for n in [1, 10, 100, 1000, 10_000, 100_000, 1_000_000]
    if n == 1
        c = ICRSCoords(2pi * rand(), pi * (rand() - 0.5))
        t1 = minimum(@benchmark(convert(GalCoords, $c)).times) / 1e9
        myprintln(io, "$n,galactic,$t1")
        t2 = minimum(@benchmark(convert(FK5Coords{2000}, $c)).times) / 1e9
        myprintln(io, "$n,fk5j2000,$t2")
        t3 = minimum(@benchmark(convert(FK5Coords{1975}, $c)).times) / 1e9
        myprintln(io, "$n,fk5j1975,$t3")
    else
        c = [ICRSCoords(2pi * rand(), pi * (rand() - 0.5)) for i = 1:n]
        t1 = minimum(@benchmark(convert(Vector{GalCoords{Float64}}, $c)).times) / 1e9
        myprintln(io, "$n,galactic,$t1")
        t2 = minimum(@benchmark(convert(Vector{FK5Coords{2000,Float64}}, $c)).times) / 1e9
        myprintln(io, "$n,fk5j2000,$t2")
        t3 = minimum(@benchmark(convert(Vector{FK5Coords{1975,Float64}}, $c)).times) / 1e9
        myprintln(io, "$n,fk5j1975,$t3")
    end
    println()
end
close(io)


text:
@testset "Preprocessing" begin

    fname = joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf")
    s = read_SSR(fname)

    @testset "Triggers" begin

        @testset "Validation" begin

            validate_triggers(s.triggers)

            s1 = deepcopy(s)
            delete!(s1.triggers, "Index")
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            delete!(s1.triggers, "Code")
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            delete!(s1.triggers, "Duration")
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            s1.triggers["test"] = [1]
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            s1.triggers["Duration"] = s1.triggers["Duration"][1:4]
            @test_throws KeyError validate_triggers(s1.triggers)

            s1 = deepcopy(s)
            s1.triggers["Code"] = s1.triggers["Code"][1:4]
            @test_throws KeyError validate_triggers(s1.triggers)

        end
    end
end

@testset "Referencing" begin

    signals = [0 1 2] .* ones(5, 3)

    @testset "Remove template" begin

        signals = [0 1 2] .* ones(5, 3)
        template = vec(2 * ones(5))
        @test remove_template(signals, template) == [-2 -1 0] .* ones(5, 3)

    end
    @testset "Reference to channel" begin

        @test rereference(signals, 3) == [-2 -1 0] .* ones(5, 3)
        @test rereference(signals, "C2", ["C1", "C2", "C3"]) == [-1 0 1] .* ones(5, 3)

    end

    @testset "Reference to group of channels" begin

        @test rereference(signals, [1, 2, 3]) == [-1 0 1] .* ones(5, 3)
        @test rereference(signals, ["C2", "C1", "C3"], ["C1", "C2", "C3"]) ==
              [-1 0 1] .* ones(5, 3)
        @test rereference(signals, "car", ["C1", "C2", "C3"]) == [-1 0 1] .* ones(5, 3)
        @test rereference(signals, "average", ["C1", "C2", "C3"]) == [-1 0 1] .* ones(5, 3)

    end

end


text:
## dew pressure solver
function x0_dew_pressure(model::EoSModel,T,y)
    #TODO
    #on sufficiently large temps,
    #the joule-thompson inversion occurs
    #making the virial coeff positive
    #on those cases, use an strategy that supposes pure gas on that side
    #Pbi = inf
    #xi = 0

    #check each T with T_scale, if treshold is over, replace Pi with inf
    pure = split_model(model)
    crit = crit_pure.(pure)

    T_c = [tup[1] for tup in crit]
    V_c = [tup[3] for tup in crit]
    _0 = zero(T+first(y))
    nan = _0/_0
    sat_nan = (nan,nan,nan)
    replaceP = ifelse.(T_c .< T,true,false)

    eachx = eachcol(Diagonal(ones(eltype(y),length(y))))
#     Bi = second_virial_coefficient.(model,T,eachx)
    #using P_B(2B) as a sat aproximation
    #z = 1 + B/v
    #P_B = RT/v(1+B/v)
    #P_B(2B) = -RT/2B(1-B/2B)
    #P_B(2B) = -0.25*RT/B
    sat = [if !replaceP[i] saturation_pressure(pure[i],T) else sat_nan end for i in 1:length(pure)]

    P_sat = [tup[1] for tup in sat]
    V_l_sat = [tup[2] for tup in sat]
    V_v_sat = [tup[3] for tup in sat]
#     P_Bi = @. -0.25*R̄*T/Bi
    #=xP0 = yP
    #dot(x,P0) = P
    P = dot(x,P0)
    =#
    P⁻¹ = zero(T)
    V0_l = zero(T)
    V0_v = zero(T)
    Pi   = zero(y)
    for i in 1:length(y)
        if !replaceP[i]
            Pi[i] = P_sat[i][1]
            P⁻¹+=y[i]/Pi[i]
            V0_v += y[i]*V_v_sat[i]
        else
            Pi[i] = pressure(pure[i],V_c[i],T)
            P⁻¹+=y[i]/Pi[i]
            V0_v += y[i]*V_c[i]*1.2
        end
    end
    #@show P_Bi
    #P = dot(x,P_Bi)
    P = 1/P⁻¹
    x = @. y*P/Pi
    xsum = 1/∑(x)
    x    = x.*xsum

    for i in 1:length(y)
        if !replaceP[i]
            V0_l += x[i]*V_l_sat[i]
        else
            V0_l += x[i]*V_c[i]
        end
    end

    prepend!(x,log10.([V0_l,V0_v]))
    return x
end

function dew_pressure(model::EoSModel, T, y; v0 =nothing)
    TYPE = promote_type(eltype(T),eltype(y))
#     lb_v = lb_volume(model,x)
    ts = T_scales(model,y)
    pmix = p_scale(model,y)
    if v0 === nothing
        v0 = x0_dew_pressure(model,T,y)
    end
    len = length(v0[1:end-1])
    #xcache = zeros(eltype(x0),len)
    Fcache = zeros(eltype(v0[1:end-1]),len)
    f!(F,z) = Obj_dew_pressure(model, F, T, exp10(z[1]), exp10(z[2]), z[3:end],y,ts,pmix)
    r  =Solvers.nlsolve(f!,v0[1:end-1],LineSearch(Newton()))
    sol = Solvers.x_sol(r)
    v_l = exp10(sol[1])
    v_v = exp10(sol[2])
    x = FractionVector(sol[3:end])
    P_sat = pressure(model,v_v,T,y)
    return (P_sat, v_l, v_v, x)
end

function Obj_dew_pressure(model::EoSModel, F, T, v_l, v_v, x, y,ts,ps)
    x   = FractionVector(x) #julia magic, check misc.jl
    μ_l = VT_chemical_potential(model,v_l,T,x)
    μ_v = VT_chemical_potential(model,v_v,T,y)
    p_l = pressure(model,v_l,T,x)
    p_v = pressure(model,v_v,T,y)
    for i in 1:length(x)
        F[i] = (μ_l[i]-μ_v[i])/(R̄*ts[i])
    end
    F[end] = (p_l-p_v)/ps
    return F
end

function dew_temperature(model,p,y)
    f(z) = Obj_dew_temperature(model,z,p,y)
    pure = split_model(model)
    sat = saturation_temperature.(pure,p)
    Ti   = zero(y)
    for i ∈ 1:length(y)
        if isnan(sat[i][1])
            Tc,pc,vc = crit_pure(pure[i])
            g(x) = p-pressure(pure[i],vc,x,[1.])
            Ti[i] = Roots.find_zero(g,(Tc))
        else
            Ti[i] = sat[i][1]
        end
    end
    T = Roots.find_zero(f,sum(Ti)/length(y))
    p,v_l,v_v,x = dew_pressure(model,T,y)
    return T,v_l,v_v,x
end

function Obj_dew_temperature(model,T,p,y)
    p̃,v_l,v_v,y = dew_pressure(model,T,y)
    return p̃-p
end


text:
nload(buses::Vector{<:Bus}) = length([b for b in buses if isa(b, LoadBus)])

ngen(buses::Vector{<:Bus}) = length([b for b in buses if isa(b, GenBus)])

count_bus_type(grid::Grid) = (nload(buses(grid)), ngen(buses(grid)))

#######################
### CHECK FUNCTIONS ###
#######################
# Everything here uses connectivty matrices as input instead of grids because
# all functions are supposed to work both for buses and for substations,
# effectively giving us two representations of the network (fine and coarse)
"""
    laplacian(con_mat::AbstractMatrix{<:Real})

Return the Laplacian matrix for a graph defined by adjacency matrix `con_mat`.
"""
function laplacian(conmat::AbstractMatrix{Bool})
    lap = zeros(Int, size(conmat))
    degree = sum(conmat, dims=1) # this gives the node degrees
    n = length(degree)
    for j in 1:(n - 1), i in (j + 1):n
        if conmat[i,j] > 0
            lap[i,j] = -1
        end
        lap[j,j] = degree[j]
    end
    lap[n,n] = degree[n]
    return Symmetric(lap, :L)
end

function laplacian(con_mat::AbstractMatrix{T}) where T <: Real
    return laplacian(con_mat .> zero(T)) # treating the substation graph as a simple graph.
end

# Shorter code but worse performance.
# function laplacian(con_mat::AbstractMatrix{<:Bool})
#     return sum(con_mat, 1) .* eye(Int, size(con_mat)[1]) - con_mat
# end

"""
    total_links(con_mat::AbstractMatrix{Bool})

Return the total number of links in a system with connectivity matrix `con_mat`.
"""
function total_links(con_mat::AbstractMatrix{Bool})
    return (1 / 2) * sum(con_mat)
end

"""
    total_links(con_mat::AbstractMatrix{Integer})

Return the total number of links in a system with connectivity matrix `con_mat`. The graph
will be treated as a simple graph.
"""
function total_links(con_mat::AbstractMatrix{T}) where T <: Real
    return total_links(con_mat .> zero(T)) # treating the substation graph as a simple graph.
end

"""
    mean_node_deg(con_mat::AbstractMatrix)

Return the average nodal degree of a system with connectivity matrix `con_mat`.
"""
function mean_node_deg(con_mat::AbstractArray{Bool, 2})
    return mean(sum(con_mat, dims=1))
end

function mean_node_deg(con_mat::AbstractArray{Int, 2})
    return mean_node_deg(con_mat .> 0) # treating the substation graph as a simple graph.
end

"""
    cluster_coeff(con_mat::AbstractMatrix)

Return the clustering coefficient of a system with connectivity matrix `con_mat`.
"""
function cluster_coeff(conmat::AbstractMatrix{Bool})
    degree = sum(conmat, 1) # this gives the node degrees
    clust_coeff = Vector{Float64}(length(degree))
    for i in 1:length(degree)
        maxedges = degree[i] * (degree[i] - 1) / 2
        cumul = 0
        for j in 1:length(degree)
            if conmat[i,j]
                cumul += dot(conmat[i,:],conmat[:,j])
            end
        end
        push!(clust_coeff, maxedges > 0 ? cumul/maxedges : 0)
    end
    return mean(clust_coeff)
end

function cluster_coeff(con_mat::AbstractMatrix{Int})
    return cluster_coeff(con_mat .> 0) # treating the substation graph as a simple graph.
end

"""
    cluster_coeff_degw(con_mat::AbstractMatrix{<:Real})

Return the degree-weighted clustering coefficient of a system with connectivity
matrix `con_mat`.
"""
function cluster_coeff_degw(conmat::AbstractMatrix{Bool})
    degree = sum(conmat, 1) # this gives the node degrees
    clust_coeff = Vector{Float64}(length(degree))
    for i in 1:length(degree)
        maxedges = degree[i] * (degree[i] - 1) / 2
        cumul = 0
        for j in 1:length(degree)
            if conmat[i, j]
                cumul += dot(conmat[i,:], conmat[:,j])
            end
        end
        push!(clust_coeff, maxedges > 0 ? degree[i] * cumul / maxedges : 0)
    end
    return mean(clust_coeff) / mean(degree)
end

function cluster_coeff_degw(con_mat::AbstractMatrix{Int})
    return cluster_coeff(con_mat .> 0) # treating the substation graph as a simple graph.
end

"""
    test_connectivity(con_mat::AbstractMatrix{<:Integer}, verb=true)

Return `true` if the system with connectivity matrix `con_mat`is connected and `false`
otherwise. If `verb=true` the result will be printed onscreen together with the Fiedler
eigenvalue.
"""
function test_connectivity(con_mat::AbstractMatrix{<:Integer}, verb=true)
    ZERO_TOL = 1e-13 # tolerance for zero due to numerical noise
    lap = laplacian(con_mat)
    evals = eigvals(lap)
    sort!(evals)
    if length(evals) > 1 && evals[2] > ZERO_TOL
        fiedler = evals[2]
        if verb
            println("Fiedler eigenvalue: $fiedler. System is connected.")
        end
        return true
    else
        if verb
            println("Fiedler eigenvalue equals zero. Network is NOT connected!")
        end
        return false
    end
end

"""
    dijkstra(con_mat::AbstractMatrix{<:Real}, v1)

Compute minimum distances to each node in the network with connectivity matrix
con_mat, starting from node v1, using Dijkstra's algorithm.
"""
function dijkstra(con_mat::AbstractMatrix{<:Real}, v1)
    n = size(con_mat, 1)
    c = v1
    visited = [v1]
    dist = [i == v1 ? 0 : Inf for i in 1:n]
    r_dist = [Inf for d in dist]
    while length(visited) < n
        neigh = [i for i in 1:n if con_mat[c, i] > 0 && !(i in visited)]
        for i in neigh
            if dist[i] > (dist[c] + con_mat[c, i])
                dist[i] = (dist[c] + con_mat[c, i])
                r_dist[i] = dist[i]
            end
        end
        c = indmin(r_dist)
        push!(visited, c)
        r_dist[c] = Inf
    end
    return dist
end

"""
    floyd_warshall(weights::AbstractMatrix{<:Real})

Compute minimum distances to each node in the network using Floyd-Warshall's algorithm. Here
we assume that the weigths matrix is symmetric.
"""
function floyd_warshall(weights::AbstractMatrix{<:Real})
    s = size(weights)
    dist = zeros(Float64, s)
    for j in 2:s[1]
        for i in 1:(j - 1)
            if weights[i, j] > 0
                dist[i, j] = weights[i, j] # we are assuming that weights is symmetric as it
                dist[j, i] = weights[i, j] # should be in our case.
            else
                dist[i, j] = Inf
                dist[j, i] = Inf
            end
        end
    end

    for k in 1:s[1]
        for i in 1:s[1]
            for j in 1:s[1]
                if dist[i, j] > dist[i, k] + dist[k, j]
                    dist[i, j] = dist[i, k] + dist[k, j]
                end
            end
        end
    end
    return dist
end

"""
    mean_shortest_path(conmat::AbstractMatrix{<:Real})

Compute the mean shortest distance between nodes for a network represented by `conmat`.
The `conmat` matrix should specify the distance between each connected node and be equal to
zero whenever two nodes are not connected. Using the adjacency matrix will return the mean
shortest path in hops. `conmat` is assumed to be symmetric. Floyd-Warshall's algorithm is
used to compute the distances.
"""
function mean_shortest_path(conmat::AbstractMatrix{<:Real})
    n = size(conmat, 1)
    return sum(floyd_warshall(conmat))/(n * n)
end

"""
    mean_shortest_path(conmat::AbstractMatrix{<:Integer}, distmat::AbstractMatrix{Float64})

Compute the mean shortest distance between nodes for a network with connectivity matrix
`conmat` and distances between nodes given by `distmat`. Distances are computed using
Floyd-Warshall's algorithm.
"""
function mean_shortest_path(
        conmat::AbstractMatrix{<:Integer},
        distmat::AbstractMatrix{Float64}
    )
    return mean_shortest_path(conmat.*distmat)
end

# """
#     mean_shortest_path(conmat::AbstractMatrix{<:Real})
#
# Compute the mean shortest distance between nodes for a network represented by `conmat`.
# The `conmat` matrix should specify the distance between each connected node and be equal to
# zero whenever two nodes are not connected. Distances are computed using Dijkstra's algorithm.
# """
# function mean_shortest_path(conmat::AbstractMatrix{<:Real})
#     n = size(conmat, 1)
#     mdist = [mean(dijkstra(conmat, i)) for i in 1:n]
#     return mean(mdist)
# end

# """
#     mean_shortest_path(conmat::AbstractMatrix{<:Real}, distmat::AbstractMatrix{<:Real})
#
# Compute the mean shortest distance between nodes for a network with connectivity matrix
# `conmat` and distances between nodes given by `distmat`. Distances are computed using
# Dijkstra's algorithm.
# """
# function mean_shortest_path(
#         conmat::AbstractMatrix{<:Real},
#         distmat::AbstractMatrix{<:Real}
#     )
#     return mean_shortest_path(conmat.*distmat)
# end

"""
    robustness_line(con_mat::AbstractMatrix{<:Integer}, n=1000)

Return the average number of transmission lines that have to be randomly removed from a
system with connectivity matrix `con_mat` before it becomes disconnected. Average is
computed over n iterations.
"""
function robustness_line(con_mat::AbstractMatrix{Bool}, n=1000)
    s = size(con_mat, 1)
    conns = []
    for i in 1:(s - 1), j in (i + 1):s
        con_mat[j, i] > 0 ? push!(conns, (i, j)) : nothing
    end
    total = 0
    for rep in 1:n
        connscopy = copy(conns)
        conmat = copy(con_mat)
        step = 0
        while test_connectivity(conmat, false)
            k = rand(1:length(connscopy))
            ij = connscopy[k]
            conmat[ij[1],ij[2]] = false
            conmat[ij[2],ij[1]] = false
            deleteat!(connscopy, k)
            step += 1
        end
        total += step
    end
    return total / n
end

function robustness_line(con_mat::AbstractMatrix{Int}, n=1000)
    return robustness_line(con_mat .> 0 , n)
end

"""
    robustness_node(con_mat::AbstractMatrix{<:Integer}, n=1000)

Return the average number of nodes that have to be randomly removed from a system with
connectivity matrix `con_mat` before it becomes disconnected. Average is computed over n
iterations.
"""
function robustness_node(con_mat::AbstractMatrix{Bool}, n=1000)
    s = size(con_mat, 1)
    total = 0
    for rep in 1:n
        len = s
        conmat = copy(con_mat)
        step = 0
        while test_connectivity(conmat, false)
            k = rand(1:len)
            mask = ones(Bool, len)
            mask[k] = false
            conmat = conmat[mask, mask]
            step += 1
            len -= 1
        end
        total += step
    end
    return total / n
end

function robustness_node(con_mat::AbstractMatrix{Int}, n=1000)
    return robustness_node(con_mat .> 0, n)
end


text:
using Documenter
using Crystals

makedocs(modules = [Crystals],
         clean = true,
         format = :html,
         sitename = "Crystals.jl",
         authors = "Mayeul d'Avezac",
         analytics = "UA-89508993-1",
         pages = Crystals._doc_pages())

deploydocs(
    repo = "github.com/mdavezac/Crystals.jl.git",
    target = "build",
    julia = "0.6",
    deps = nothing,
    make = nothing,
)


text:
using DSGE
using Base.Test, DataFrames, HDF5

path = dirname(@__FILE__)

@testset "Test various forms of data loading" begin
    # Can we actually test? Require that FRED API key exists
    if haskey(ENV, "FRED_API_KEY") || isfile(joinpath(ENV["HOME"],".freddatarc"))

        # Specify vintage and dates
        custom_settings = Dict{Symbol, Setting}(
                :data_vintage            => Setting(:data_vintage, "160812"),
                :cond_vintage            => Setting(:cond_vintage, "160812"),
                :cond_id                 => Setting(:cond_id, 0),
                :use_population_forecast => Setting(:use_population_forecast, true),
                :date_forecast_start     => Setting(:date_forecast_start, quartertodate("2016-Q3")),
                :date_conditional_end    => Setting(:date_forecast_start, quartertodate("2016-Q3")),
                :n_anticipated_shocks    => Setting(:n_anticipated_shocks, 6))

        m = Model990(custom_settings = custom_settings, testing = true)

        # Read expected results
        exp_data, exp_cond_data, exp_semicond_data =
            h5open("$path/../reference/load_data_out.h5", "r") do h5
                read(h5, "data"), read(h5, "cond_data"), read(h5, "semicond_data")
            end

        # Unconditional data
        println("The following warnings are expected test behavior:")
        df = load_data(m; try_disk=false, verbose=:none)
        data = df_to_matrix(m, df)
        @test @test_matrix_approx_eq exp_data data

        # Conditional data
        cond_df = load_data(m; cond_type=:full, try_disk=false, verbose=:none)
        cond_data = df_to_matrix(m, cond_df; cond_type=:full)
        @test @test_matrix_approx_eq exp_cond_data cond_data

        # Semiconditional data
        semicond_df = load_data(m; cond_type=:semi, try_disk=false, verbose=:none)
        semicond_data = df_to_matrix(m, semicond_df; cond_type=:semi)
        @test @test_matrix_approx_eq exp_semicond_data semicond_data
    else
        warn("Skipping load_data test because FRED_API_KEY not present")
    end
end

nothing


text:
using DrWatson
@quickactivate "AccuracyAtTop_DeepTopPush"

using Tar, ProgressMeter, MAT, CSV, DataFrames, DelimitedFiles

# extraction
const IMAGENET_DIR = "/disk/macha/data_deeptoppush/datasets/ImageNet2012"
const IMAGENET_VERSION = "ILSVRC2012"

function extract_train(file)
    dir = joinpath(IMAGENET_DIR, "train")
    tmp1 = joinpath(IMAGENET_DIR, "tmp")
    mkpath(dir)
    Tar.extract(file, tmp1)

    @showprogress for file in readdir(tmp1; join = true)
        tmp2 = Tar.extract(file)
        datadir = joinpath(dir, basename(file)[1:end-4])
        mkpath(datadir)
        files = readdir(tmp2)
        mv.(joinpath.(tmp2, files), joinpath.(datadir, files))
    end
    rm(tmp1; recursive = true)
    return
end

function extract_test(file)
    dir = joinpath(IMAGENET_DIR, "test")
    mkpath(dir)
    Tar.extract(file, dir)
    return
end

extract_train(IMAGENET_DIR * "/$(IMAGENET_VERSION)_img_train.tar")
extract_test(IMAGENET_DIR * "/$(IMAGENET_VERSION)_img_val.tar")

# label extraction
function create_label_map(meta::String)
    m = matread(meta)["synsets"]
    return Dict(String.(m["WNID"]) .=> Int.(m["$(IMAGENET_VERSION)_ID"]))
end

function data_train(label_map, image_dir = IMAGENET_DIR)
    dir_train = joinpath(image_dir, "train")

    dfs = map(readdir(dir_train)) do dir
        return DataFrame(
            files = readdir(joinpath(dir_train, dir); join = true),
            labels = label_map[dir],
            category = dir,
        )
    end
    df = reduce(vcat, dfs)
    CSV.write(joinpath(image_dir, "data_train.csv"), df)
    return df
end

function data_test(file, label_map, image_dir = IMAGENET_DIR)
    dir_test = joinpath(image_dir, "test")
    rev_map = Dict(values(label_map) .=> keys(label_map))
    labels = readdlm(file, ' ', Int, '\n')[:]

    df = DataFrame(
        files = readdir(dir_test; join = true),
        labels = labels,
        category = [rev_map[label] for label in labels],
    )
    CSV.write(joinpath(image_dir, "data_test.csv"), df)
    return df
end

label_map = create_label_map(IMAGENET_DIR * "/meta.mat")

data_train(label_map)
data_test(IMAGENET_DIR * "/$(IMAGENET_VERSION)_validation_ground_truth.txt", label_map)

# reshape all images
using Images, ImageMagick
using Base.Threads
using PyCall

py"""
from PIL import Image
import numpy

def convert_to_rgb(file):
    img = Image.open(file).convert('RGB')
    img.save(file)
    return img
"""

function load_image(file, resize)
    img = RGB.(ImageMagick.load(file))
    if !isnothing(resize) && size(img) != resize
        img = Images.imresize(img, resize)
    end
    return img
end

df = vcat(
    CSV.read(joinpath(IMAGENET_DIR, "data_test.csv"), DataFrame),
    CSV.read(joinpath(IMAGENET_DIR, "data_train.csv"), DataFrame),
)

@time Threads.@threads for file in files
    file_new = replace(file, "ImageNet2012" => "ImageNet224")
    file_new = replace(file_new, "JPEG" => "png")
    if !isfile(file_new)
        try
            img = load_image(file, (224, 224))
        catch
            py"convert_to_rgb"(file)
            img = load_image(file, (224, 224))
        end
        mkpath(dirname(file_new))
        Images.save(file_new, img)
    end
end

image_dir = "/disk/macha/data_deeptoppush/datasets/ImageNet224"

data_train(label_map, image_dir)
data_test(IMAGENET_DIR * "/$(IMAGENET_VERSION)_validation_ground_truth.txt", label_map, image_dir)


text:
# Command Line Components
import Base: show
export AbstractOption, Option, ShortOption

"""
    AbstractOption{T}
"""
abstract type AbstractOption{T} end

isoption(s::AbstractString) = startswith(s, "-")
islong(s::AbstractString) = startswith(s, "--")
isoption(::Void) = false
islong(::Void) = false

import Base: ==, in, parse

function in(item::Union{AbstractString, Char}, collection::Vector{AbstractOption})
    for each in collection
        if item == each
            return true
        end
    end
    return false
end

struct Option{T} <: AbstractOption{T}
    name::String
    key::Symbol
    doc::String
end

Option(::Type{T}, name::String, key::Symbol, doc::String="") where T =
    Option{T}(name, key, doc)
Option(::Type{T}, key::Symbol, doc::String="") where T =
    Option(T, string(key), key, doc)

Option(name::String, key::Symbol, doc::String="") =
    Option(Bool, name, key, doc)
Option(key::Symbol, doc::String="") = Option(Bool, key, doc)

optionstring(op::Option{T}) where T = "--$(op.name)[=$(T)]"
optionstring(op::Option{Bool}) = "--$(op.name)"

==(s::Char, op::Option) = false
==(s::AbstractString, op::Option) = (s == op.name) || (s == "--"*op.name)

function show(io::IO, op::Option{T}) where T
    print(io, "--$(op.name)[=$(T)]")
    print(io, "\t$(op.doc)")
end

function show(io::IO, op::Option{Bool})
    print(io, "--$(op.name)")
    print(io, "\t$(op.doc)")
end

"""
    ShortOption{T} <: AbstractOption{T}

Option with Shorthand. e.g option (flag) `help`
can be declared by following commands.

```shell
> cmd -h
> cmd --help
```
"""
struct ShortOption{T} <: AbstractOption{T}
    short::Char
    long::String
    key::Symbol
    doc::String
end

ShortOption(::Type{T}, short::Char, long::String, key::Symbol, doc::String="") where T =
    ShortOption{T}(short, long, key, doc)
ShortOption(::Type{T}, long::String, key::Symbol, doc::String="") where T =
    ShortOption(T, long[1], long, key, doc)
ShortOption(::Type{T}, key::Symbol, doc::String="") where T =
    ShortOption(T, string(key), key, doc)

ShortOption(short::Char, long::String, key::Symbol, doc::String="") =
    ShortOption(Bool, short, long, key, doc)
ShortOption(long::String, key::Symbol, doc::String="") =
    ShortOption(long[1], long, key, doc)
ShortOption(key::Symbol, doc::String="") =
    ShortOption(string(key), key, doc)

optionstring(op::ShortOption{T}) where T = "-$(op.short)[$(T)],--$(op.long)[=$(T)]"
optionstring(op::ShortOption{Bool}) = "-$(op.short),--$(op.long)"

==(s::Char, op::ShortOption) = (s == op.short)
==(s::AbstractString, op::ShortOption) =
    length(s) == 1? s[1] == op : ((s == op.long) || (s == "--$(op.long)")) || (s == "-$(op.short)")

function show(io::IO, op::ShortOption{T}) where T
    print(io, "-$(op.short)[$(T)],--$(op.long)[=$(T)]")
    print(io, "\t$(op.doc)")
end

function show(io::IO, op::ShortOption{Bool})
    print(io, "-$(op.short),--$(op.long)")
    print(io, "\t$(op.doc)")
end

# Preserved Options
const help = ShortOption(:help, "help information")
const version = ShortOption(:version, "version information")

# Signature
export Signature

struct Signature
    args::Vector{DataType}
    options::Vector{AbstractOption}
    doc::String
end

Signature(args::Vector{DataType}, options::Vector) =
    Signature(args, [options..., help], "\tNo Documents")

function show(io::IO, op::Signature)
    if length(op.args) > 0
        println("[arguments]:")
        println(io, "$(op.doc)")
    end

    println("[options]:")
    if length(op.options) > 0
        for each in op.options
            print(io, "  ")
            println(io, each)
        end
    end
end

# Parsers
"""
--flag is the shorthands for --flag=true
"""
function parse(val::AbstractString, opt::AbstractOption{Bool})
    val == "" && return opt.key, true
    return opt.key, parse(Bool, val)
end

function parse(val::AbstractString, opt::AbstractOption{T}) where {T <: Union{Integer, AbstractFloat}}
    return opt.key, parse(T, val)
end

function parse(val::AbstractString, opt::AbstractOption{T}) where {T <: AbstractString}
    return opt.key, T(val)
end

function parse(val::AbstractString, opt::AbstractOption{T}) where {T}
    value = try
        eval(parse(val))
    catch e
        throw(CLIError("Invalid Argument: $val"))
    end

    if !isa(value, T)
        throw(CLIArgError(T, typeof(value)))
    end
    return opt.key, value
end


text:
##############################################################################
#
# All Document types share a common metadata profile using DocumentMetadata
#
##############################################################################

mutable struct DocumentMetadata
    language
    title::String
    author::String
    timestamp::String
end

"""
    DocumentMetadata(language, title::String, author::String, timestamp::String)

Stores basic metadata about Document.

...
# Arguments
- `language`: What language is the document in? Defaults to Languages.English(), a Language instance defined by the Languages package.
- `title::String` : What is the title of the document? Defaults to "Untitled Document".
- `author::String` : Who wrote the document? Defaults to "Unknown Author".
- `timestamp::String` : When was the document written? Defaults to "Unknown Time".
...
"""
DocumentMetadata() = DocumentMetadata(
    Languages.English(),
    "Untitled Document",
    "Unknown Author",
    "Unknown Time"
)

##############################################################################
#
# The abstract Document type
#
##############################################################################

abstract type AbstractDocument; end

mutable struct FileDocument <: AbstractDocument
    filename::String
    metadata::DocumentMetadata
end

"""
    FileDocument(pathname::AbstractString)

Represents a document using a plain text file on disk.

# Example
```julia-repl
julia> pathname = "/usr/share/dict/words"
"/usr/share/dict/words"

julia> fd = FileDocument(pathname)
A FileDocument
 * Language: Languages.English()
 * Title: /usr/share/dict/words
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: A A's AMD AMD's AOL AOL's Aachen Aachen's Aaliyah
```
"""
function FileDocument(pathname::AbstractString)
    doc = FileDocument(String(pathname), DocumentMetadata())
    doc.metadata.title = pathname
    return doc
end

mutable struct StringDocument{T<:AbstractString} <: AbstractDocument
    text::T
    metadata::DocumentMetadata
end

"""
    StringDocument(txt::AbstractString)

Represents a document using a UTF8 String stored in RAM.

# Example
```julia-repl
julia> str = "To be or not to be..."
"To be or not to be..."

julia> sd = StringDocument(str)
A StringDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: To be or not to be...
```
"""
StringDocument(txt::AbstractString) = StringDocument(txt, DocumentMetadata())

mutable struct TokenDocument{T<:AbstractString} <: AbstractDocument
    tokens::Vector{T}
    metadata::DocumentMetadata
end

"""
    TokenDocument(txt::AbstractString)
    TokenDocument(txt::AbstractString, dm::DocumentMetadata)
    TokenDocument(tkns::Vector{T}) where T <: AbstractString

Represents a document as a sequence of UTF8 tokens.

# Example
```julia-repl
julia> my_tokens = String["To", "be", "or", "not", "to", "be..."]
6-element Array{String,1}:
    "To"
    "be"
    "or"
    "not"
    "to"
    "be..."

julia> td = TokenDocument(my_tokens)
A TokenDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: ***SAMPLE TEXT NOT AVAILABLE***
```
"""
function TokenDocument(txt::AbstractString, dm::DocumentMetadata)
    TokenDocument(tokenize(dm.language, String(txt)), dm)
end
function TokenDocument(tkns::Vector{T}) where T <: AbstractString
    TokenDocument(tkns, DocumentMetadata())
end
TokenDocument(txt::AbstractString) = TokenDocument(String(txt), DocumentMetadata())

mutable struct NGramDocument{T<:AbstractString} <: AbstractDocument
    ngrams::Dict{T,Int}
    n::Union{Int,Vector{Int}}
    metadata::DocumentMetadata
end

"""
    NGramDocument(txt::AbstractString, n::Integer=1)
    NGramDocument(txt::AbstractString, dm::DocumentMetadata, n::Integer=1)
    NGramDocument(ng::Dict{T, Int}, n::Integer=1) where T <: AbstractString

Represents a document as a bag of n-grams, which are UTF8 n-grams and map to counts.

# Example
```julia-repl
julia> my_ngrams = Dict{String, Int}("To" => 1, "be" => 2,
                                     "or" => 1, "not" => 1,
                                     "to" => 1, "be..." => 1)
Dict{String,Int64} with 6 entries:
  "or"    => 1
  "be..." => 1
  "not"   => 1
  "to"    => 1
  "To"    => 1
  "be"    => 2

julia> ngd = NGramDocument(my_ngrams)
A NGramDocument{AbstractString}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: ***SAMPLE TEXT NOT AVAILABLE***
```
"""
function NGramDocument(txt::AbstractString, dm::DocumentMetadata, n::Integer...=1)
    NGramDocument(ngramize(dm.language, tokenize(dm.language, String(txt)), n...), (length(n) == 1) ? Int(first(n)) : Int[n...], dm)
end
function NGramDocument(txt::AbstractString, n::Integer...=1)
    NGramDocument(txt, DocumentMetadata(), n...)
end
function NGramDocument(ng::Dict{T, Int}, n::Integer...=1) where T <: AbstractString
    NGramDocument(merge(Dict{AbstractString,Int}(), ng), (length(n) == 1) ? Int(first(n)) : Int[n...], DocumentMetadata())
end

##############################################################################
#
# text() / text!(): Access to document text as a string
#
##############################################################################
"""
    text(fd::FileDocument)
    text(sd::StringDocument)
    text(ngd::NGramDocument)

Access the text of Document as a string.

# Example
```julia-repl
julia> sd = StringDocument("To be or not to be...")
A StringDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: To be or not to be...

julia> text(sd)
"To be or not to be..."
```
"""
function text(fd::FileDocument)
    !isfile(fd.filename) && error("Can't find file: $(fd.filename)")
    read(fd.filename, String)
end

text(sd::StringDocument) = sd.text
function text(td::TokenDocument)
    @warn("TokenDocument's can only approximate the original text")
    join(td.tokens, " ")
end
function text(ngd::NGramDocument)
    error("The text of an NGramDocument cannot be reconstructed")
end

text!(sd::StringDocument, new_text::AbstractString) = (sd.text = new_text)
function text!(d::AbstractDocument, new_text::AbstractString)
    error("The text of a $(typeof(d)) cannot be edited")
end

##############################################################################
#
# tokens() / tokens!(): Access to document text as a token array
#
##############################################################################
"""
    tokens(d::TokenDocument)
    tokens(d::(Union{FileDocument, StringDocument}))

Access the document text as a token array.

# Example
```julia-repl
julia> sd = StringDocument("To be or not to be...")
A StringDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: To be or not to be...

julia> tokens(sd)
7-element Array{String,1}:
    "To"
    "be"
    "or"
    "not"
    "to"
    "be.."
    "."
```
"""
tokens(d::(Union{FileDocument, StringDocument})) = tokenize(language(d), text(d))
tokens(d::TokenDocument) = d.tokens
function tokens(d::NGramDocument)
    error("The tokens of an NGramDocument cannot be reconstructed")
end

tokens!(d::TokenDocument, new_tokens::Vector{T}) where {T <: AbstractString} = (d.tokens = new_tokens)
function tokens!(d::AbstractDocument, new_tokens::Vector{T}) where T <: AbstractString
    error("The tokens of a $(typeof(d)) cannot be directly edited")
end

##############################################################################
#
# ngrams() / ngrams!(): Access to document text as n-gram counts
#
##############################################################################
"""
    ngrams(ngd::NGramDocument, n::Integer)
    ngrams(d::AbstractDocument, n::Integer)
    ngrams(d::NGramDocument)
    ngrams(d::AbstractDocument)

Access the document text as n-gram counts.

# Example
```julia-repl
julia> sd = StringDocument("To be or not to be...")
A StringDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: To be or not to be...

julia> ngrams(sd)
 Dict{String,Int64} with 7 entries:
  "or"   => 1
  "not"  => 1
  "to"   => 1
  "To"   => 1
  "be"   => 1
  "be.." => 1
  "."    => 1
```
"""
function ngrams(ngd::NGramDocument, n::Integer)
    error("The n-gram complexity of an NGramDocument cannot be increased")
end
ngrams(d::AbstractDocument, n::Integer...) = ngramize(language(d), tokens(d), n...)
ngrams(d::NGramDocument) = d.ngrams
ngrams(d::AbstractDocument) = ngrams(d, 1)

ngrams!(d::NGramDocument, new_ngrams::Dict{AbstractString, Int}) = (d.ngrams = new_ngrams)
function ngrams!(d::AbstractDocument, new_ngrams::Dict)
    error("The n-grams of $(typeof(d)) cannot be directly edited")
end

##############################################################################
#
# Length describes length of document in characters
#
##############################################################################

function Base.length(d::NGramDocument)
    error("NGramDocument's do not have a well-defined length")
end

Base.length(d::AbstractDocument) = length(text(d))

##############################################################################
#
# Length describes length of document in characters
#
##############################################################################

ngram_complexity(ngd::NGramDocument) = ngd.n

function ngram_complexity(d::AbstractDocument)
    error("$(typeof(d))'s have no n-gram complexity")
end

##############################################################################
#
# Union type that refers to a generic, non-abstract document type
#
##############################################################################

const GenericDocument = Union{
    FileDocument,
    StringDocument,
    TokenDocument,
    NGramDocument
}

##############################################################################
#
# Easier Document() constructor that decides types based on inputs
#
##############################################################################

Document(str::AbstractString) = isfile(str) ? FileDocument(str) : StringDocument(str)
Document(tkns::Vector{T}) where {T <: AbstractString} = TokenDocument(tkns)
Document(ng::Dict{String, Int}) = NGramDocument(ng)

##############################################################################
#
# Conversion rules
#
##############################################################################

function Base.convert(::Type{StringDocument}, d::FileDocument)
    StringDocument(text(d), d.metadata)
end
function Base.convert(::Type{TokenDocument}, d::(Union{FileDocument, StringDocument}))
    TokenDocument(tokens(d), d.metadata)
end
function Base.convert(::Type{NGramDocument},
            d::(Union{FileDocument, StringDocument, TokenDocument}))
    NGramDocument(ngrams(d), 1, d.metadata)
end
Base.convert(::Type{TokenDocument}, d::TokenDocument) = d
Base.convert(::Type{NGramDocument}, d::NGramDocument) = d

##############################################################################
#
# getindex() methods: StringDocument("This is text and that is not")["is"]
#
##############################################################################

Base.getindex(d::AbstractDocument, term::AbstractString) = ngrams(d)[term]


text:
using Test
using SumOfSquares
import MultivariateBases
const MB = MultivariateBases

function xor_complement_test()
    @test Certificate.Sparsity.xor_complement([1], 1) == Int[]
    @test Certificate.Sparsity.xor_complement(Int[], 1) == [1]
    @test Certificate.Sparsity.xor_complement([1], 2) == [2]
    @test Certificate.Sparsity.xor_complement([2], 2) == [1]
    @test Certificate.Sparsity.xor_complement([1, 2], 2) == Int[]
    @test Certificate.Sparsity.xor_complement([1, 3], 2) == Int[]
    @test Certificate.Sparsity.xor_complement(Int[], 2) == [1, 2]
    @test Certificate.Sparsity.xor_complement([7], 3) == [3, 5]
    @test Certificate.Sparsity.xor_complement([5, 6, 3], 3) == [7]
    @test Certificate.Sparsity.xor_complement([3], 3) == [3, 4]
end

set_monos(bases::Vector{<:MB.MonomialBasis}) = Set([basis.monomials for basis in bases])

function Certificate.Sparsity.sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, domain::SemialgebraicSets.BasicSemialgebraicSet, sp::Sparsity.Monomial, maxdegree, degs)
    half_monos = Certificate.maxdegree_gram_basis(MB.MonomialBasis, variables, div(maxdegree, 2))
    P = Set(monos)
end

"""
    wml19()

Examples of [MWL19].

[WML19] Wang, Jie, Victor Magron, and Jean-Bernard Lasserre. "TSSOS: A Moment-SOS hierarchy that exploits term sparsity." arXiv preprint arXiv:1912.08899 (2019).
"""
function wml19()
    certificate = Certificate.Newton(SOSCone(), MB.MonomialBasis, tuple())
    @testset "Example 4.2" begin
        @polyvar x[1:3]
        f = 1 + x[1]^4 + x[2]^4 + x[3]^4 + prod(x) + x[2]
        expected_1_false = Set(monovec.([
            [x[3]^2],
            [x[1] * x[3], x[2]],
            [x[2], 1],
            [x[2]^2],
            [x[2] * x[3], x[1]],
            [x[1] * x[2], x[3]],
            [x[1]^2]
        ]))
        expected_2 = Set(monovec.([
            [x[1]^2, x[2]^2, x[3]^2, 1],
            [x[2], 1],
            [x[2] * x[3], x[1]],
            [x[1] * x[3], x[2]],
            [x[1] * x[2], x[3]]
        ]))
        cluster_expected_1_false = Set(monovec.([
            [1, x[2], x[1]*x[3]],
            [x[1]^2],
            [x[2]^2],
            [x[3]^2],
            [x[1], x[2]*x[3]],
            [x[3], x[1]*x[2]],
        ]))
        cluster_expected_1 = Set(monovec.([
            [1, x[2], x[1]^2, x[2]^2, x[1]*x[3], x[3]^2],
            [x[1], x[2]*x[3]],
            [x[3], x[1]*x[2]],
        ]))
        cluster_expected_2 = Set(monovec.([
            [1, x[2], x[1]^2, x[2]^2, x[1]*x[3], x[3]^2],
            [x[1], x[2]*x[3], x[3], x[1]*x[2]],
        ]))
        @testset "$completion $k $use_all_monomials" for completion in [ClusterCompletion(), ChordalCompletion()], k in 0:2, use_all_monomials in [false, true]
            if completion isa ClusterCompletion
                expected = k == 1 ? (use_all_monomials ? cluster_expected_1 : cluster_expected_1_false) : cluster_expected_2
            else
                expected = (k == 1 && !use_all_monomials) ? expected_1_false : expected_2
            end
            @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(completion, k, use_all_monomials), certificate)) == expected
        end
        expected = Set(monovec.([
            [x[1]^2, x[1] * x[3], x[2]^2, x[3]^2, x[2], 1],
            [x[1] * x[2], x[2] * x[3], x[1], x[3]]
        ]))
        @test set_monos(Certificate.Sparsity.sparsity(f, SignSymmetry(), certificate)) == expected
    end
    @testset "Example 5.4" begin
        preorder_certificate = Certificate.Putinar(Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 4), SOSCone(), MB.MonomialBasis, 4)
        @polyvar x[1:2]
        f = x[1]^4 + x[2]^4 + x[1] * x[2]
        K = @set 1 - 2x[1]^2 - x[2]^2 >= 0
        @testset "$completion $k $use_all_monomials" for completion in [ClusterCompletion(), ChordalCompletion()], k in 0:2, use_all_monomials in [false, true]
            basis, preorder_bases = Certificate.Sparsity.sparsity(f, K, Sparsity.Monomial(completion, k, use_all_monomials), preorder_certificate)
            if k == 1 && (!use_all_monomials || completion isa ChordalCompletion)
                if use_all_monomials
                    @test set_monos(preorder_bases[1]) == Set(monovec.([[x[1], x[2]], [constantmonomial(x[1] * x[2])]]))
                    @test set_monos(basis) == Set(monovec.([[x[1], x[2]], [x[1] * x[2], 1], [x[1]^2, x[2]^2, 1]]))
                else
                    @test set_monos(preorder_bases[1]) == Set(monovec.([[x[1], x[2]]]))
                    @test set_monos(basis) == Set(monovec.([[x[1]^2], [x[1]*x[2], 1], [x[2]^2], [x[1], x[2]]]))
                end
            else
                @test set_monos(preorder_bases[1]) == Set(monovec.([[constantmonomial(x[1] * x[2])], [x[1], x[2]]]))
                @test set_monos(basis) == Set(monovec.([[x[1]^2, x[1]*x[2], x[2]^2, 1], [x[1], x[2]]]))
            end
        end
    end
    @testset "Example 6.7" begin
        @polyvar x[1:2]
        f = 1 + x[1]^2 * x[2]^4 + x[1]^4 * x[2]^2 + x[1]^4 * x[2]^4 - x[1] * x[2]^2 - 3x[1]^2 * x[2]^2
        @testset "$completion $k $use_all_monomials" for completion in [ClusterCompletion(), ChordalCompletion()], k in 0:2, use_all_monomials in [false, true]
            expected = if completion isa ClusterCompletion
                Set(monovec.([[x[1]^2 * x[2]^2, x[1] * x[2]^2, 1], [x[1] * x[2]], [x[1]^2 * x[2]]]))
            else
                Set(monovec.([[x[1] * x[2]^2, 1], [x[1]^2 * x[2]^2, 1], [x[1] * x[2]], [x[1]^2 * x[2]]]))
            end
            @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(completion, k, use_all_monomials), certificate)) == expected
        end
        @test set_monos(Certificate.Sparsity.sparsity(f, SignSymmetry(), certificate)) == Set(monovec.([
            [x[1]^2 * x[2]^2, x[1] * x[2]^2, 1], [x[1]^2 * x[2], x[1] * x[2]]
        ]))
    end
end

"""
    l09()

Examples of [MWL19].

[L09] Lofberg, Johan. "Pre-and post-processing sum-of-squares programs in practice." IEEE transactions on automatic control 54.5 (2009): 1007-1011.
"""
function l09()
    certificate = Certificate.Newton(SOSCone(), MB.MonomialBasis, tuple())
    @testset "Example 1 and 2" begin
        @polyvar x[1:2]
        f = 1 + x[1]^4 * x[2]^2 + x[1]^2 * x[2]^4
        @test Certificate.monomials_half_newton_polytope(monomials(f), tuple()) == [
            x[1]^2 * x[2], x[1] * x[2]^2, 1
        ]
        expected = Set(monovec.([
            [x[1]^2 * x[2]], [x[1] * x[2]^2], [constantmonomial(x[1] * x[2])]
        ]))
        for i in 0:2
            @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(ChordalCompletion(), i), certificate)) == expected
        end
        @test set_monos(Certificate.Sparsity.sparsity(f, SignSymmetry(), certificate)) == expected
    end
    @testset "Example 3 and 4" begin
        @polyvar x[1:3]
        f = 1 + x[1]^4 + x[1] * x[2] + x[2]^4 + x[3]^2
        @testset "$k $use_all_monomials" for k in 0:2, use_all_monomials in [false, true]
            if k == 1
                if use_all_monomials
                    @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), certificate)) == Set(monovec.([
                        [x[1]^2, x[2]^2, 1], [x[1], x[2]], [x[1] * x[2], 1], [x[3]]
                    ]))
                else
                    @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), certificate)) == Set(monovec.([
                        [x[1]^2], [x[2]^2], [x[1], x[2]], [x[1] * x[2], 1], [x[3]]
                    ]))
                end
            else
                @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), certificate)) == Set(monovec.([
                    [x[1], x[2]], [x[3]], [x[1]^2, x[2]^2, 1], [x[1] * x[2], 1]
                ]))
            end
        end
        @test set_monos(Certificate.Sparsity.sparsity(f, SignSymmetry(), certificate)) == Set(monovec.([
            [x[1], x[2]], [x[3]], [x[1]^2, x[1] * x[2], x[2]^2, 1]
        ]))
    end
end
function square_domain()
    d = 6
    preorder_certificate = Certificate.Putinar(Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 6), SOSCone(), MB.MonomialBasis, 6)
    @polyvar x y
    f = x^2*y^4 + x^4*y^2 - 3*x^2*y*2 + 1
    K = @set(1 - x^2 >= 0 && 1 - y^2 >= 0)
    @testset "Square domain $k $use_all_monomials" for k in 0:4, use_all_monomials in [false, true]
        basis, preorder_bases = Certificate.Sparsity.sparsity(f, K, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), preorder_certificate)
        if k == 1
            if use_all_monomials
                @test set_monos(basis) == Set(monovec.([[x^2 * y, y, 1], [x^3, x * y^2, x], [x^2, y, 1], [x^2, y^2, 1], [x^2 * y, y^3, y], [x * y, x]]))
            else
                @test set_monos(basis) == Set(monovec.([[x^2 * y, y^3], [x^2, y], [x^2 * y, 1], [x^3, x * y^2], [x * y, x]]))
            end
        elseif k == 2
            @test set_monos(basis) == Set(monovec.([[x^2 * y, x^2, y^2, 1], [x^3, x * y^2, x * y, x], [x^2 * y, x^2, y, 1], [x^2 * y, y^3, x^2, y]]))
        elseif k == 3
            @test set_monos(basis) == Set(monovec.([[x^3, x * y^2, x * y, x], [x^2 * y, x^2, y^2, y, 1], [x^2 * y, y^3, x^2, y, 1]]))
        else
            @test set_monos(basis) == Set(monovec.([[x^3, x * y^2, x * y, x], [x^2 * y, y^3, x^2, y^2, y, 1]]))
        end
        expected = Set(monovec.([[x^2, y^2, y, 1], [x * y, x]]))
        if k == 1
            if use_all_monomials
                @test set_monos(preorder_bases[1]) == Set(monovec.([[x * y, x], [x^2, y, 1], [x^2, y^2, 1]]))
            else
                @test set_monos(preorder_bases[1]) == Set(monovec.([[y, 1], [x * y, x], [x^2, y], [x^2, y^2]]))
            end
        else
            @test set_monos(preorder_bases[1]) == expected
        end
        if k == 1
            if use_all_monomials
                @test set_monos(preorder_bases[2]) == Set(monovec.([[x * y, x], [x^2, y], [x^2, y^2, 1]]))
            else
                @test set_monos(preorder_bases[2]) == Set(monovec.([[x * y, x], [x^2, y], [x^2, y^2], [constantmonomial(x * y)]]))
            end
        elseif k == 2
            @test set_monos(preorder_bases[2]) == Set(monovec.([[x^2, y^2, 1], [x^2, y, 1], [x * y, x]]))
        else
            @test set_monos(preorder_bases[2]) == expected
        end
    end

end
function sum_square(n)
    @testset "Sum square" begin
        @polyvar x[1:(2n)]
        certificate = Certificate.Newton(SOSCone(), MB.MonomialBasis, tuple())
        f = sum((x[1:2:(2n-1)] .- x[2:2:(2n)]).^2)
        expected = Set(monovec.([monovec([x[(2i - 1)], x[2i], 1]) for i in 1:n]))
        @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Variable(), Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 2))) == expected
        expected = Set(monovec.([[x[(2i - 1)], x[2i]] for i in 1:n]))
        @test set_monos(Certificate.Sparsity.sparsity(f, SignSymmetry(), certificate)) == expected
    end
end
function drop_monomials()
    @testset "Drop monomials" begin
        @polyvar x
        f = polynomial(x^2)
        certificate = Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 2)
        @testset "$k $use_all_monomials" for k in 0:2, use_all_monomials in [false, true]
            # The monomial `1˘ is dropped as it is useless.
            if use_all_monomials
                expected = Set(monovec.([[x], [constantmonomial(x^2)]]))
            else
                expected = Set([monovec([x])])
            end
            @test set_monos(Certificate.Sparsity.sparsity(f, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), certificate)) == expected
        end
        preorder_certificate = Certificate.Putinar(Certificate.MaxDegree(SOSCone(), MB.MonomialBasis, 4), SOSCone(), MB.MonomialBasis, 3)
        f = polynomial(x^3)
        K = @set x >= 0
        @testset "$k $use_all_monomials" for k in 0:3, use_all_monomials in [false, true]
            basis, preorder_bases = Certificate.Sparsity.sparsity(f, K, Sparsity.Monomial(ChordalCompletion(), k, use_all_monomials), preorder_certificate)
            if k == 1 && !use_all_monomials
                @test set_monos(basis) == Set(monovec.([[x^2, x]]))
            elseif (k == 2 && !use_all_monomials) || (k == 1 && use_all_monomials)
                @test set_monos(basis) == Set(monovec.([[x^2, 1], [x^2, x]]))
            else
                @test set_monos(basis) == Set(monovec.([[x^2, x, 1]]))
            end
            if k == 1 && !use_all_monomials
                @test set_monos(preorder_bases[1]) == Set(monovec.([[x]]))
            else
                @test set_monos(preorder_bases[1]) == Set(monovec.([[x, 1]]))
            end
        end
    end
end
@testset "Sparsity" begin
    xor_complement_test()
    wml19()
    l09()
    square_domain()
    sum_square(8)
    @test Certificate.Sparsity.appropriate_type(32) == Int64
    sum_square(32)
    @test Certificate.Sparsity.appropriate_type(64) == Int128
    sum_square(64)
    @test Certificate.Sparsity.appropriate_type(128) == BigInt
    sum_square(128)
    drop_monomials()
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares