text:
## Morikawa and Fujiwara 2013 GMPE tests
@testset "Morikawa and Fujiwara 2013 GMPE PGA" begin
  # init model parameters
  include("../examples/morikawa-fujiwara-2013.conf")
  ## PGA at epicenter on grid M7.0, Dl = 250 (constant)
  @test gmpe_mf2013(eq_7,config_mf2013_crustal_pga,grid_epicenter)[1].pga == 53.28
  ## PGA,PGV,PSA on grid, M6, ASID false, Dl - constant
  S_c = gmpe_mf2013(eq_6,config_mf2013_crustal_pga,grid)
  @test length(S_c) == TEST_GRID_SIZE
  @test typeof(S_c) == Array{GroundMotion.Point_pga_out,1}
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 3.4
  S_c = gmpe_mf2013(eq_6,config_mf2013_crustal_pgv,grid)
  @test typeof(S_c) == Array{GroundMotion.Point_pgv_out,1}
  @test round(sum([S_c[i].pgv for i=1:length(S_c)]),digits=2) == 4.61
  S_c = gmpe_mf2013(eq_6,config_mf2013_crustal_psa_03,grid)
  @test typeof(S_c) == Array{GroundMotion.Point_psa_out,1}
  S_c = gmpe_mf2013(eq_6,config_mf2013_crustal_psa_10,grid)
  @test round(sum([S_c[i].psa for i=1:length(S_c)]),digits=2) == 5.45
  S_c = gmpe_mf2013(eq_6,config_mf2013_crustal_psa_30,grid)
  @test round(sum([S_c[i].psa for i=1:length(S_c)]),digits=2) == 1.43
  ## PGA modeling on grid M8.5, withoit minpga, ASID false, Dl - constant
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_pga,grid)
  @test length(S_c) == TEST_GRID_SIZE
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 37.6
  ## PGA modeling on grid M8.5, withoit minpga, ASID true, Dl - constant
  S_c = gmpe_mf2013(eq_85,config_mf2013_intraplate_pga_asid,grid)
  @test length(S_c) == TEST_GRID_SIZE
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 74.4
  ## run PGA modeling on grid M6, withoit minpga, ASID false, Dl on GRID
  S_c = gmpe_mf2013(eq_6,config_mf2013_crustal_pga,grid_dl)
  @test length(S_c) == TEST_GRID_SIZE
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 506.54
  ## run PGA modeling on grid M8.5, withoit minpga, ASID true, Dl on GRID
  S_c = gmpe_mf2013(eq_85,config_mf2013_intraplate_pga_asid,grid_dl,Xvf=40)
  @test length(S_c) == TEST_GRID_SIZE
  @test round(sum([S_c[i].pga for i=1:length(S_c)]),digits=2) == 4110.09
  ## PGV,PSA modeling on grid M6, withoit minpga, ASID false, Dl on GRID
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_pgv,grid_dl)
  @test length(S_c) == TEST_GRID_SIZE
  @test typeof(S_c) == Array{GroundMotion.Point_pgv_out,1}
  @test round(sum([S_c[i].pgv for i=1:length(S_c)]),digits=2) == 2989.47
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_psa_03,grid_dl)
  @test length(S_c) == TEST_GRID_SIZE
  @test typeof(S_c) == Array{GroundMotion.Point_psa_out,1}
  @test round(sum([S_c[i].psa for i=1:length(S_c)]),digits=2) == 4177.5
  ## PGA,PGV,PSA without grid Dl and VS30 by default, ASID=false
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_pga)
  @test length(S_c) == SIMULATION_ARRAY_SIZE
  @test round(sum(S_c),digits=2) == 5468.53
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_pgv)
  @test round(sum(S_c),digits=2) == 7710.7
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_psa_03)
  @test round(sum(S_c),digits=2) == 13116.13
  ## PGA without grid Dl=500 VS30 by default, ASID=false
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_pga,Dl=500)
  @test round(sum(S_c),digits=2) == 5725.89
  ## PGA without grid Dl by default VS30=500, ASID=false
  S_c = gmpe_mf2013(eq_85,config_mf2013_crustal_pga,VS30=500)
  @test round(sum(S_c),digits=2) == 4790.89
  ## ## PGA without grid Dl and VS30 by default, ASID=true, Xvf=40
  S_c = gmpe_mf2013(eq_85,config_mf2013_intraplate_pga_asid,Xvf=40)
  @test round(sum(S_c),digits=2) == 12845.95
end


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
reorder = 1
periods = [0,1,0]
dims = MPI.Dims_create(nnodes, [0,0,0])
comm_cart = MPI.Cart_create(comm, dims; periodic=periods, reorder=reorder)

rank = MPI.Comm_rank(comm)
ccoords = MPI.Cart_coords(comm_cart, rank)
@test all(ccoords .>= 0)
@test all(ccoords .< dims)

coords = MPI.Cart_coords(comm_cart)
@test all(coords .>= 0)
@test all(coords .< dims)

MPI.free(comm_cart)

MPI.Finalize()
# @test MPI.Finalized()

ext: jl


text:
export IFS

"""
`IFS` is an iterated function system.
"""
immutable IFS
  funcs::Vector{AffineMap}
end

function IFS(args...)
  fun_list = [ f for f in args ]
  IFS(fun_list)
end

# Apply F to a single vector
function (F::IFS)(x::Vector)
  y = [ f(x) for f in F.funcs ]
  return Set(y)
end

# Apply F to a list of vectors
function (F::IFS)(xlist::Set)
  A = Set{Vector{Float64}}()
  for x in xlist
    A = union(A,F(x))
  end
  return A
end

function show(io::IO,F::IFS)
  nf = size(F.funcs)[1]
  for k=1:nf
    println(io,"#$k:\t$(F.funcs[k])")
  end
end

"""
`square_check(F)` checks if all functions `f` in the `IFS` pass
`square_check(f)`.
"""
function square_check(F::IFS)
  return all(square_check(f) for f in F.funcs)
end

"""
`Triangle_IFS()` returns an `IFS` that creates a Triangle triangle.
"""
function Triangle_IFS()
  A = 0.5 * eye(2)
  f1 = AffineMap(A,[0,0])
  f2 = AffineMap(A,[0.5,0])
  f3 = AffineMap(A,[0.25,0.5])
  return IFS(f1,f2,f3)
end

"""
`Carpet_IFS()` returns an `IFS` for Triangle's Carpet.
"""
function Carpet_IFS()
  A = (1/3)*eye(2)
  f1 = AffineMap(A,[0,0])
  f2 = AffineMap(A,[0,1/3])
  f3 = AffineMap(A,[0,2/3])
  g1 = AffineMap(A,[1/3,0])
  g3 = AffineMap(A,[1/3,2/3])
  h1 = AffineMap(A,[2/3,0])
  h2 = AffineMap(A,[2/3,1/3])
  h3 = AffineMap(A,[2/3,2/3])
  return IFS(f1,f2,f3,g1,g3,h1,h2,h3)
end

"""
`visualize(F)` for an `IFS` draws the unit square (dotted)
and the image of the functions in `F` (solid frames).
"""
function visualize(F::IFS)
  clf()
  for f in F.funcs
    if !square_check(f)
      warn("$f does not map unit square to itself.")
    end
    visualize(f)
  end
end


text:
module BinanceHistoricalData
    using Dates, HTTP, JSON3
    export download_klines_spot, download_aggtrades_spot, symbol_info

    const FREQS = ["1m","5m","15m","30m","1h","2h","4h","6h","8h","12h","1d","3d","1w","1mo"]
    const AGGLEVELS = ["monthly","daily"]
    const SYMTYPE = ["spot","futures"]

    const tbl_colnames_klines = ["OpenTime","OpenPrice","HighPrice","LowPrice","ClosePrice","Volume","CloseTime","QuoteVolume","NTrades","TakerBuyBaseVolume","TakerBuyQuoteVolume","IgnoreFlag"]
    const tbl_colnames_aggtrades = ["AggTradeId","Price","Quantity","FirstTradeId","LastTradeId","Timestamp","IsMarketMakerBuyer","IsTradeBestPriceMatch"]

    """
        download_klines_spot(sym,date::Date,write_dir;freq="1m",agg_level="monthly",verbose=true)

    Download historical klines data starting on `date` for symbol `sym` to directory `write_dir`.

    Other argument descriptions:
     - Data frequency `freq`. (See `BinanceHistoricalData.FREQS` for possible values )
     - Data aggregation level `agg_level`. (See `BinanceHistoricalData.AGGLEVELS` for possible values )

    Downloads zipped CSV with column names `BinanceHistoricalData.tbl_colnames_klines`.
    See [this link](https://github.com/binance/binance-public-data) for more details
    on the CSV file's contents.

    """
    function download_klines_spot(sym,date::Date,write_dir;freq="1m",agg_level="monthly",verbose=true)
        # Check arguments
        (agg_level ∈ AGGLEVELS) || error("must have agg_level ∈ $AGGLEVELS")
        (freq ∈ FREQS) || error("must have freq ∈ $FREQS")
        # Format information into strings
        sym_str = uppercase(String(sym))
        date_str = (agg_level == "monthly") ? Dates.format(Date(date),"yyyy-mm") : Dates.format(Date(date),"yyyy-mm-dd")
        # Put together full data url
        sym_type = "spot"
        data_url = "https://data.binance.vision/data/$sym_type/$agg_level/klines/$sym/$freq/$sym_str-$freq-$date_str.zip"
        # Download the data
        verbose && @info "[$(now())] - Fetching binance historical kline data" sym, date, freq, sym_type, agg_level
        return download(data_url,joinpath(write_dir,"kline-$sym_type-$sym_str-$freq-$date_str.zip"))
    end

    """
        download_aggtrades_spot(sym,date::Date,write_dir::AbstractString;sym_type="spot",agg_level="monthly",verbose=true)

    Download historical aggregate trade data starting on `date` for symbol `sym` to directory `write_dir`.

    Other argument descriptions:
     - Data frequency `freq`. (See `BinanceHistoricalData.FREQS` for possible values )
     - Data aggregation level `agg_level`. (See `BinanceHistoricalData.AGGLEVELS` for possible values )

    Downloads zipped CSV with column names `BinanceHistoricalData.tbl_colnames_aggtrades`.
    See [this link](https://github.com/binance/binance-public-data) for more details
    on the CSV file's contents.
    """
    function download_aggtrades_spot(sym,date::Date,write_dir::AbstractString; sym_type="spot",agg_level="monthly",verbose=true)
        # Check arguments
        (agg_level ∈ AGGLEVELS) || error("must have agg_level ∈ $AGGLEVELS")
        # Format information into strings
        sym_str = String(sym)
        date_str = (agg_level == "monthly") ? Dates.format(Date(date),"yyyy-mm") : Dates.format(Date(date),"yyyy-mm-dd")
        # Put together full data url
        data_url = "https://data.binance.vision/data/spot/$agg_level/aggTrades/$sym_str/$sym_str-aggTrades-$date_str.zip"
        # Download the data
        verbose && @info "[$(now())] - Fetching binance historical aggtrade data" sym, date, sym_type, agg_level
        return download(data_url,joinpath(write_dir,"aggTrades-$sym_type-$sym_str-$date_str.zip"))
    end

    """
        symbol_info()

    Get table of all available symbols on the Binance and some of their properties.
    """
    function symbol_info()
        resp = HTTP.get("https://api.binance.com/api/v3/exchangeInfo")
        js_tmp = JSON3.read(resp.body)
        return js_tmp.symbols
    end

end


text:
#
#--------------------------------------#
export ABu
#--------------------------------------#
"""
 Tensor product operator
 (As ⊗ Br) * u
"""
function ABu(As::AbstractArray,Br::AbstractArray,u::AbstractArray)

m,n = size(u);

Bu = u;
if(length(Br)!=0)
    mb,nb = size(Br);
    m     = Int(m*mb/nb);
    Bu    = reshape(Bu,nb,:)
    Bu    = Br * Bu;
    Bu    = reshape(Bu,m,n);
end

ABu = Bu;
if(length(As)!=0)
    ma,na = size(As);
    Ey    = n/na;
    n     = Int(Ey*ma);
    tmp   = zeros(eltype(u),m,n);
    for i=1:Ey
        ii=Int((i-1)*ma+1):Int(i*ma);
        jj=Int((i-1)*na+1):Int(i*na);
        tmp[:,ii] = ABu[:,jj]*As';
    end
    ABu = tmp
end

return ABu
end
#--------------------------------------#
Zygote.@adjoint function ABu(As,Br,u)
return ABu(As,Br,u),dv->(nothing,nothing,ABu(As',Br',dv));
end
#--------------------------------------#
export ABu1
#--------------------------------------#
"""
 (As ⊗ Br) * u
"""
function ABu1(As,Br,u)

    m = size(As)[1]
    n = size(Br)[1]
    E = size(u )[3]

    v = zeros(m,n,E)

    ABu1!(v,As,Br,u,E=E)

return v
end
#--------------------------------------#
export ABu1!
#--------------------------------------#
"""
 v = (As ⊗ Br) * u

 To adopt this (general) implementation of the tensor product operation,
 we need to get the following to work with arrays of shape (nx,ny,ne)

 1. geometry: x,y <= semreshape
 2. mask: <= construct in square domain, and then do sem_reshape
 3. gather scatter ops <= global ⟺   local numbering and do sem_reshape;

 sem_reshape -- should only be done when generating Mesh object

 focus on deformed([-1,1]²) for now. think about general geomtries later

"""
function ABu1!(v,As,Br,u)

    @views for ie=1:size(E)[end]
        v[:,:,ie] .= Br * u[:,:,ie] * As'
    end

return
end
#--------------------------------------#


text:
using Test
using Random
using NeXLCore

@testset "NeXLCore Tests" begin
    include("element.jl")
    include("xray.jl")
    include("material.jl")
    include("matparser.jl")
    include("stoichiometry.jl")
    include("matu.jl")
    include("kratio.jl")
    include("eta.jl")
    include("custommacs.jl")
    include("bethe.jl")
    include("meanionizationpotential.jl")
    include("electron.jl")
    include("scattering.jl")
    include("staging.jl")
    include("standardize.jl")
end

# To test coverage (from NWMR's Win10 box)
# 1. Open the command line
# 2. cd c:\Users\nritchie\.julia\dev\NeXLCore\test
# 3. > C:\Users\nritchie\AppData\Local\Programs\Julia\Julia-1.4.1\bin\julia.exe --code-coverage=user
# 4. julia> using Pkg
# 5. julia> Pkg.test("NeXLCore",coverage=true)
# 6. Look at .cov files in src for test coverage data


text:
using TuringModels
#using StatsPlots

ProjDir = @__DIR__

# Define a simple Normal model with unknown mean and variance.

@model gdemo(x, y) = begin
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))
end

#  Run sampler, collect results

chn = sample(gdemo(1.5, 2), NUTS(0.65), 1000)

# Summarise results (currently requires the master branch from MCMCChains)

chn |> display

# Plot and save results if in ./dev

#p = plot(chn)
#savefig("$(ProjDir)/basic-example-plot.png")


text:
module Format

using TextWrap

export printf,
    printfT, printHeader, printEnd, printLine, formatFrac, toFrac, interval, center

LEFT_WIDTH = 50
RIGHT_WIDTH = 25
LINE_WIDTH = LEFT_WIDTH + RIGHT_WIDTH

# Formats and outputs a line of text
function printf(text...; line = false)
    line && println()
    println_wrapped(text..., width = LINE_WIDTH)
    line && printLine()
end

# Formats and outputs a line of text with a theorem reference
function printfT(theorem, text...)
    wrapped = join(
        map(
            line -> wrap(string(line), width = LEFT_WIDTH, break_on_hyphens = false),
            [text...],
        ),
        "\n",
    )
    lastPos = occursin("\n", wrapped) ? findlast('\n', wrapped) : 0
    println()
    print(rpad(wrapped, lastPos + LEFT_WIDTH))
    println(lpad(theorem, RIGHT_WIDTH))
    printLine()
end

# Formats and outputs text as a section header
function printHeader(text...)
    println("\n")
    println(text...)
    printLine()
end

# Outputs the end of a procedure
function printEnd()
    println()
    println(center("END"))
end

# Outputs a section separator
function printLine()
    println('—'^LINE_WIDTH)
end

# Formats a Rational type variable into its string representation
function formatFrac(frac::Rational{Int64}, den = denominator(frac))
    (n, d) = (numerator(frac), denominator(frac))
    if n % d == 0 && den == d
        return string(Int64(n / d))
    end
    m = den % d == 0 ? Int64(den / d) : 1
    string(m * n, "/", m * d)
end

# Converts a string representation into a Rational type variable
function toFrac(frac::String)
    try
        f = vcat(map(x -> parse(Int64, x), split(frac, "/")), [1])
        f[1] // f[2]
    catch e
        Int64(frac)
    end
end

# Outputs a string interval that concatenates given ranges labeled with given labels
function interval(rngs...; labels = repeat([" "], length(rngs) - 1))
    lines = []
    w = maximum(length(rng[2]) for rng in rngs) + 1
    int = center(rngs[1][1], width = w)
    val = center(rngs[1][2], width = w)
    for i = 2:length(rngs)
        label = labels[i-1]
        rng = rngs[i]
        int *= label
        val *= ' '^length(label) * center(rng[2], width = w)
        if i < length(rngs) && length(int * labels[i]) + 2w > LINE_WIDTH
            int *= center(rng[1][1:1], width = w)
            append!(lines, [center(int)])
            append!(lines, [center(val)])
            int = center(rng[1][length(rng[1]):length(rng[1])], width = w)
            val = center(rng[2], width = w)
        else
            int *= center(rng[1], width = w)
        end
    end
    append!(lines, [center(int)])
    append!(lines, [center(val)])
    join(lines, "\n")
end

# Centers text within width spaces
function center(text...; width::Int64 = LINE_WIDTH)
    out = []
    half = Int64(floor(width / 2))
    for line in text
        mid = Int64(floor(length(line) / 2))
        append!(
            out,
            [lpad(line[1:mid], half) * rpad(line[(mid+1):length(line)], width - half)],
        )
    end
    join(out, "\n")
end

# Helper function for printfT -- Unicode-compatible version of built-in findlast
function findlast(char::Char, str::String)
    o = 0
    e = 1
    i = 1
    while i <= length(str)
        if str[e] == char
            o = i
        end
        e += nextind(string(str[e]), 1) - 1
        i += 1
    end
    o
end

end


text:
module CoverRanges #src
export CompInf, min_covering, solve_min_covering #src

include("access_dict.jl") #src
export AccessDict, total_gets, total_sets #src

# The idea of the algorithm is that we can simply start from the left-most element
# of the range `r` that we wish to cover, and decide whether or not to include it
# in our partial cover. If we do want to include it, the optimal range to use to
# cover it is the one which extends the furthest to the right (out of all the
# ranges which include the element). Then it remains to solve a smaller instance
# of the problem. Likewise, if we don't include the element, we again are left
# with a smaller instance of the problem. In other words, we aim to solve the
# dynamic programming equation
#
# ```julia
# min_number_to_cover(a:b, K) = min(1 + min_number_to_cover(f:b, K - (f-a+1)), min_number_to_cover((a+1):b, K))
# ```
#
# where
# `min_number_to_cover(r, K)` is the minimum number of ranges needed to cover at
# least `K` elements of `r` and `f` is the endpoint furthest to the right of any
# range in `C` that includes `a`.
#
# Why is this the right recursive equation for this problem? Let us say we reach
# node `a`. We can either
#
# * include node `a` in our set of nodes covered, in which case we choose the
#   range that includes `a` and stretches as far as possible to the right; in that
#   case, we have to solve the problem on the remaining set of nodes (this
#   corresponds to `min_number_to_cover(f:b, K - (f-a+1))`).
# * Otherwise, we skip node `a` and need to solve the problem with the remaining
# nodes `a+1:b`,
#   still needing to cover `K` nodes (this corresponds to `min_number_to_cover(
#   (a+1):b, K))`).
#

# Let's begin the implementation.
# First, we need an object to represent an impossible task, which is never chosen in a
# minimum and which is invariant under adding `+1`. Julia's `Inf` does this for
# us, but `Inf` is a floating point number, and we'd prefer to stick to integers
# (with the exception of this infinity).
struct ComparisonInf end
const CompInf = ComparisonInf()
Base.isless(::ComparisonInf, ::ComparisonInf) = false

for T in (Int8, Int16, Int32, Int64, Int128)
    @eval Base.isless(::ComparisonInf, ::$T) = false
    @eval Base.isless(::$T, ::ComparisonInf) = true
    @eval Base.:+(::$T, c::ComparisonInf) = c
    @eval Base.:+(c::ComparisonInf, ::$T) = c
end

# In the above list, we could add other types if needed, such as
# [SaferIntegers.jl](https://github.com/JeffreySarnoff/SaferIntegers.jl)
# numeric types to ensure our algorithm does not overflow.

# Next, we define an struct to hold an instance or subinstance of our problem.
struct LeftToCover{T}
    "Range to partially cover"
    r::UnitRange{T}
    "Number of nodes we are required to cover"
    K::T
end

# And likewise, one to hold a solution (or part of one).
struct PartialCover{T}
    "Minimum number needed to partially cover the range"
    m::Union{T, ComparisonInf}
    "Index of latest range added"
    i::Int
    "Key to point towards next range in the partial cover"
    key::LeftToCover{T}
end

# Next, to find the range with the furthest right endpoint
# including a given node, we just write a simple loop.
"""
    find_furthest_range(a, C) -> index, right_endpoint

Given a node `a`, finds the range `r` in `C` such that `r` contains `a` and has
the largest right endpoint.
"""
function find_furthest_range(a, C)
    ## This could be optimized more by sorting `C` first and terminating early
    best_f_so_far = typemin(a) # right endpoint
    best_i_so_far = nothing # index
    for (i, r) in enumerate(C)
        a ∈ r || continue
        if last(r) > best_f_so_far
            best_f_so_far = last(r)
            best_i_so_far = i
        end
    end
    return best_i_so_far, best_f_so_far
end

# Now comes the heart of the algorithm. We compute
#
# ```julia
# min_number_to_cover(a:b, K) = min(1 + min_number_to_cover(f:b, K - (f-a+1)), min_number_to_cover((a+1):b, K))
# ```
#
# recursively with the help of a memoization cache, which maps from
# (sub)instances of the problem to solutions.
function min_covering(key::LeftToCover{T}, furthest_ranges, cache::AbstractDict{LeftToCover{T}, PartialCover{T}}) where {T}
    get!(cache, key) do
        r = key.r;  K = key.K
        K <= 0 && return PartialCover{T}(zero(T), 0, key) # Nothing left to cover
        K > length(r) && return PartialCover{T}(CompInf, 0, key) # Impossible task
        a, b = first(r), last(r)
        i_a, f_a = furthest_ranges[a] # index and right endpoint of optimal range including `a`

        ## Solve the subproblem in which we skip `a`
        skip_cover = min_covering(LeftToCover{T}(a+1:b, K), furthest_ranges, cache)
        i_a === nothing && return skip_cover # cannot cover `a` means we must skip it

        ## Solve the subproblem in which we include `a` in our partial cover
        add_key = LeftToCover{T}( (f_a+1) : b, K - (f_a - a + 1))
        add_cover = min_covering(add_key, furthest_ranges, cache)

        ## Decide whether or not to include `a` and return the better result
        if add_cover.m + 1 <= skip_cover.m
            return PartialCover{T}(add_cover.m + T(1), i_a, add_key)
        else
            return skip_cover
        end
    end
end

# Note that in the case that we choose to skip `a`, we return `skip_cover`, which
# is the solution for the subproblem that skips `a`. But when we choose to include
# `a` in the partial cover, we make a new `PartialCover` object instead of
# returning `add_cover`. That's because if we include `a` by using the range that
# ends at `f_a`, we really are adding a range to the partial cover represented by
# `add_cover`, so we make a new `PartialCover` object, one whose value of `m` is one
# larger (corresponding to the extra range we added), which holds the index `i_a`
# of the range that we added, and which holds `add_key`, which we can use to
# recover `add_cover` from `cache` later. In all, this lets us hold the sequence
# of ranges used in the final optimal cover efficiently (as a linked list).

# Lastly, we wrap `min_covering` in a helper function that constructs `furthest_ranges`,
# the memoization cache, and unwinds the linked list that holds the solution in order
# to recover the indices of the ranges included in the minimal partial cover we have found.
"""
    solve_min_covering(r::UnitRange{T}, K, C; [cache]) -> Union{Nothing, Vector{Int}}

Returns either `nothing` or a vector of indices of `C` corresponding to a
minimal set of ranges which covers at least `K` elements of `r`. Returns
`nothing` if and only if the task is impossible.

* `r`: range to partially cover
* `K`: number of elements of `r` that must be
covered
* `C`: vector of ranges to use in constructing the partial cover
* `cache`: optionally provide an `AbstractDict{LeftToCover{T}, PartialCover{T}}`
which is used as a memoization cache

!!! warning
    Overflow is not checked.
"""
function solve_min_covering(r::UnitRange{T}, K, C; cache =  Dict{LeftToCover{T},
PartialCover{T}}()) where {T}
    ## `furthest_ranges` is all we need from `C`
    furthest_ranges = Dict(a => find_furthest_range(a, C) for a in r)

    ## Recursive function to compute the covering
    result = min_covering(LeftToCover{T}(r, K), furthest_ranges, cache)

    k = result.m # minimum number of ranges needed
    k == CompInf && return nothing

    ## We store the list of ranges essentially as a linked list.
    ## To recover the optimal list of ranges, we traverse through
    ## the cache, at each step getting the `key` for the next element in the cache.
    range_inds = Int[]
    sizehint!(range_inds, k)
    while k > 1
        push!(range_inds, result.i)
        k = result.m
        result = cache[result.key]
    end
    return range_inds
end

end #src


text:
function getTableauSLRK(s, o, tsym, q, p, q̃, p̃, ω, d=Nothing)
    TableauSLRK(tsym, o, s, q, p, q̃, p̃, ω, d)
end

function TableauSLRKLobattoIII(s, q, p)
    o = 2s-2
    tsym = Symbol("SLRKLobattoIII")
    getTableauSLRK(s, o, tsym, q, p, q, p, get_lobatto_ω_matrix(s), get_lobatto_nullvector(s))
end

function TableauSLRKLobattoIIIAB(s)
    q = TableauLobattoIIIA(s)
    p = TableauLobattoIIIB(s)
    TableauSLRKLobattoIII(s, q, p)
end

function TableauSLRKLobattoIIIBA(s)
    q = TableauLobattoIIIB(s)
    p = TableauLobattoIIIA(s)
    TableauSLRKLobattoIII(s, q, p)
end

function TableauSLRKLobattoIIICC̄(s)
    q = TableauLobattoIIIC(s)
    p = TableauLobattoIIIC̄(s)
    TableauSLRKLobattoIII(s, q, p)
end

function TableauSLRKLobattoIIIC̄C(s)
    q = TableauLobattoIIIC̄(s)
    p = TableauLobattoIIIC(s)
    TableauSLRKLobattoIII(s, q, p)
end

function TableauSLRKLobattoIIID(s)
    l = TableauLobattoIIID(s)
    TableauSLRKLobattoIII(s, l, l)
end

function TableauSLRKLobattoIIIE(s)
    l = TableauLobattoIIIE(s)
    TableauSLRKLobattoIII(s, l, l)
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct GatePipe <: Pipe
        name::String
        arc::Arc
        available::Bool
        diameter::Float64
        length::Float64
        roughness::Float64
        flowlimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        open_status::Union{Nothing,Bool}
        flow::Union{Nothing,Float64}
        headloss::Union{Nothing,Float64}
        internal::InfrastructureSystemsInternal
    end

A bidirectional-flow pipe with a gate (shutoff) valve.

# Arguments
- `name::String`
- `arc::Arc`
- `available::Bool`
- `diameter::Float64`
- `length::Float64`
- `roughness::Float64`
- `flowlimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `open_status::Union{Nothing,Bool}`
- `flow::Union{Nothing,Float64}`
- `headloss::Union{Nothing,Float64}`
- `internal::InfrastructureSystemsInternal`: internal reference, do not modify
"""
mutable struct GatePipe <: Pipe
    name::String
    arc::Arc
    available::Bool
    diameter::Float64
    length::Float64
    roughness::Float64
    flowlimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    open_status::Union{Nothing,Bool}
    flow::Union{Nothing,Float64}
    headloss::Union{Nothing,Float64}
    "internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function GatePipe(name, arc, available, diameter, length, roughness, flowlimits, open_status, flow, headloss, )
    GatePipe(name, arc, available, diameter, length, roughness, flowlimits, open_status, flow, headloss, InfrastructureSystemsInternal(), )
end

function GatePipe(; name, arc, available, diameter, length, roughness, flowlimits, open_status, flow, headloss, )
    GatePipe(name, arc, available, diameter, length, roughness, flowlimits, open_status, flow, headloss, )
end

# Constructor for demo purposes; non-functional.
function GatePipe(::Nothing)
    GatePipe(;
        name="init",
        arc=Arc(nothing),
        available=true,
        diameter=0.0,
        length=0.0,
        roughness=0.0,
        flowlimits=(min=-10.0, max=10.0),
        open_status=nothing,
        flow=nothing,
        headloss=nothing,
    )
end

"""Get GatePipe name."""
get_name(value::GatePipe) = value.name
"""Get GatePipe arc."""
get_arc(value::GatePipe) = value.arc
"""Get GatePipe available."""
get_available(value::GatePipe) = value.available
"""Get GatePipe diameter."""
get_diameter(value::GatePipe) = value.diameter
"""Get GatePipe length."""
get_length(value::GatePipe) = value.length
"""Get GatePipe roughness."""
get_roughness(value::GatePipe) = value.roughness
"""Get GatePipe flowlimits."""
get_flowlimits(value::GatePipe) = value.flowlimits
"""Get GatePipe open_status."""
get_open_status(value::GatePipe) = value.open_status
"""Get GatePipe flow."""
get_flow(value::GatePipe) = value.flow
"""Get GatePipe headloss."""
get_headloss(value::GatePipe) = value.headloss
"""Get GatePipe internal."""
get_internal(value::GatePipe) = value.internal


text:
"""
Maximum Gap Tree
================

Also called **Gap Tree Lasso (GapLas)**.

In the central function `gaplas`, a tree having the largest gap
values (computed by `gap_vec!`) is selected.  The non-tree edge-flows
are forwarded into the input `y`.  Then the tree solver is used and
the tree-edges are updated.

Graph Indexes
-------------

For the different steps, several indexes are necessary
1. [`Dual.gap_vec!`](@ref) needs to access every edge once.

2. Prim's minimum spanning tree: `IncidenceIndex` providing outgoing
   edges for a node (i.e. every edge is included in two node's
   adjacency lists)

3. Once the tree has been determined we need to have a `ChildrenIndex`
   from the `parent` (e.g. for `dfs_walk`); also needed for computing
   the `Queues` layout (parenthesis tree representation) and the
   processing order.
"""
module MGT

import Statistics
import Printf: @sprintf

import GraphIdx
import GraphIdx: Const, Vec, Ones, Weights
import GraphIdx: PrimMstMem, prim_mst_edges, Graph, EdgeGraph, Edge
import GraphIdx: WeightedGraph, enumerate_edges, num_edges, num_nodes
import GraphIdx.Tree: RootedTree

import ..TreeDP: TreeDPMem, tree_dp!
import ..Utils: sum2, primal_objective
import ..Dual: dual!, gap_vec!, primal_from_dual!
import ..TreeFlow: update_tree!, extract_non_tree!

"""
    gaplas(y, graph, λ; learn=1.0, max_iter=3, ...)

Optimize in each iteration along a tree.

The `learn` parameter controls how much of the new tree
solution should be taken (should be between ``0`` and ``1.0``).
"""
function gaplas(
    y::Array{Float64,N},
    graph::EdgeGraph,
    lambda::Weights{Float64};
    root_node::Int = 1,
    mu::Wmu = Ones{Float64}(),
    max_iter::I = 3,
    verbose::Bool = true,
    process::Fu1 = x->nothing,
    dprocess::Fu2 = α->nothing,
    tprocess::Fu3 = (t,w)->nothing,
    learn::Float64 = 1.0,
)::Array{Float64,N} where {N, I<:Integer, Fu1<:Function, Fu2<:Function, Fu3<:Function, Wmu}
    local m = num_edges(graph)
    local n = length(y)
    @assert n == num_nodes(graph)
    local alpha = zeros(m)
    local γ = zeros(m)
    local x = copy(y)
    local z = similar(y)
    local tlam = Vector{Float64}(undef, n)
    local dp_mem = TreeDPMem(n)
    local mst_mem = PrimMstMem(graph)
    local selected = mst_mem.selected
    local parent = mst_mem.parent
    local tree = RootedTree(root_node, parent)
    local wgraph::WeightedGraph = WeightedGraph(mst_mem.neighbors, collect(lambda, m))
    local x_new::Array{Float64,N} = learn >= 1.0 ? x : copy(x)
    local alpha_new::Array{Float64} = learn >= 1.0 ? alpha : copy(alpha)
    if verbose
        local γ_sorted = similar(γ)
        local x2 = similar(x)
    end

    for it in 1:max_iter
        @timed begin
        gap_vec!(γ, x, alpha, wgraph, -1.0)
        @timed if verbose
            if n < 30 && it > 1 && false
                tree_gamma_check(γ, alpha, tlam, selected,
                                 x, z, dp_mem.proc_order, parent)
            end
            local quant = Statistics.quantile!(γ_sorted .= γ .* -1,
                                               [0.90, 0.95, 0.98])
            local gap = -sum(γ)
            local dual_obj = 0.5*sum2(x)
            local prim_obj = primal_objective(x, y, wgraph)
            primal_from_dual!(x2 .= y, alpha, wgraph)
            # @assert x ≈ x2 "‖x - x2‖_∞ = $(maximum(abs.(x - x2)))"
            @assert prim_obj + dual_obj - gap ≈ 0.5sum2(y)

            print(@sprintf("%4d %12.4f %12.4f %12.4f",
                           it,
                           gap,
                           dual_obj,
                           prim_obj,
                           ))

            if false
                print(@sprintf("%8f %8f %8f", quant...))
            end
            println()
        end

        @timed prim_mst_edges(γ, root_node, mst_mem)
        tprocess(γ, parent)
        z .= y
        extract_non_tree!(graph, parent, z, alpha, tlam, lambda)
        x_new .= x
        @timed tree_dp!(x_new, z, tree, Vec(tlam), mu, dp_mem)
        if !(x === x_new)
            x .= (1 - learn) .* x .+ learn .* x_new
        end
        process(x)
        let tree_alpha = tlam   # alpha within the tree (tlam is not needed)
            tree_alpha .= vec(x_new) .- vec(z)
            dual!(tree_alpha, dp_mem.proc_order, parent)
            alpha_new .= alpha
            update_tree!(alpha_new, tree_alpha, selected, graph, parent)
            if !(alpha === alpha_new)
                alpha .= (1 - learn) .* alpha .+ learn .* alpha_new
            end
            if false
                duality_check(alpha, lambda)
            end
        end
        dprocess(alpha)
        end
    end

    return x
end

@deprecate(
    gaplas(y::Array, edges, λ::Array{Float64}; args::Args...) where {Args},
    gaplas(y, edges, GraphIdx.Vec(λ); args...),
)

@deprecate(
    gaplas(y::Array, edges, λ::Float64; args::Args...) where {Args},
    gaplas(y, edges, GraphIdx.Const(λ); args...),
)

gaplas(
    y::Array, edges::Vector{Edge{Int}}, λ::Weights{Float64}; args::Args...
) where {Args} =
    gaplas(y, EdgeGraph(length(y), edges), λ; args...)

"""
    gaplas(...)

Graph has to implement several methods:
- `iter_edges(::Function, ::Graph)`
- `IncidenceIndex(::Graph)`
"""
function gaplas(
    y::Array{F,N},
    g::Graph,
    λ::Weights,
    μ::Weights = Ones{F}();
    args::Args...
)::Array{F,N} where {F, N, Args}
    edge_graph::EdgeGraph = collect(g)
    gaplas(y, edge_graph, λ; mu=μ, args...)
end

"""
    gaplas!(...)

Perform one iteration. Premises:
1. `x = y + D'*α`
2. `α` dually feasible, `abs(α[e]) ≤ λ[e]` for all edges `e`.
"""
function gaplas!(
    x::Array{Float64}, α::Vector{Float64}, γ::Vector{Float64}, y::Array{Float64}, g::Graph
)
    error("notimplemented")
end

"""
Ensure that the alpha is dually feasible, i.e.
componentwise absolutely not greater than lambda.
"""
function duality_check(alpha, lambda)
    for i in 1:length(alpha)
        @assert(abs(alpha[i]) <= (1 + 1e-8)*lambda[i],
                "$i: $u--$v " *
                "$(alpha[i]) < $(lambda[i])?")
    end
end

function tree_gamma_check(γ, alpha, tlam, selected, x, z, proc_order, parent, edges)
    let γ = round.(-γ, digits=2)
        # @show γ
    end
    m = length(alpha)

    @assert selected[1] < 0
    for i in @view selected[2:end]
        if abs(γ[i]) >= 1e-15
            talp = vec(x) - vec(z)
            for i in @view proc_order[1:end-1]
                talp[parent[i]] += talp[i]
            end
            alp = zeros(m)
            for i in @view selected[2:end]
                u, v = edges[i]
                alp[i] = parent[u] == v ? talp[u] : -talp[v]
            end
            for i in 1:m
                println(@sprintf("%7.3f  %7.3f", alp[i], alpha[i]))
            end
            @assert(
                abs(γ[i]) < 1e-15,
                "i=$i: $(edges[i])\n γ=$(γ[i])\n α=$(alpha[i])" *
                "\n $(x[edges[i][1]]) $(x[edges[i][2]])" *
                "\n $(dp_mem.proc_order)" *
                "\n $(y≈z)" *
                "\n $(talp ≈ tlam)" * ""
            )
        end
    end
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

# using SymPy
using Plots

include("linearcode.jl")

# using additive over stabilizer
abstract type AbstractAdditiveCode <: AbstractCode end
abstract type AbstractQuantumCode <: AbstractAdditiveCode end
abstract type AbstractCSSCode <: AbstractQuantumCode end
# can also build a AbstractQuantumLinearCode if the additive code is also linear
# will probably have to tweak all of these later

mutable struct CSSCode <: AbstractCSSCode
    F::FqNmodFiniteField # base field (symplectic)
    E::FqNmodFiniteField # additive field
    n::Integer
    k::Integer
    d::Union{Integer, Missing}
    dx::Union{Integer, Missing}
    dz::Union{Integer, Missing}
    stabs::fq_nmod_mat
    Xstabs::fq_nmod_mat
    Zstabs::fq_nmod_mat
    Xorigcode::Union{LinearCode, Missing}
    ZorigCode::Union{LinearCode, Missing}
    signs::Vector{Int64} # make bools, uint8's?
    Xsigns::Vector{Int64} # make bools, uint8's?
    Zsigns::Vector{Int64} # make bools, uint8's?
    dualgens::fq_nmod_mat
    logspace::Union{fq_nmod_mat, Missing}
    logicals::Union{Vector{Tuple{fq_nmod_mat, fq_nmod_mat}}, Missing}
    charvec::Vector{Int64} # make bools, uint8's?
    dualweightenum::Union{WeightEnumerator, Missing}
    logsweightenum::Union{WeightEnumerator, Missing}
    Pauliweightenum::Union{WeightEnumerator, Missing}
end

mutable struct QuantumCode <: AbstractQuantumCode
    F::FqNmodFiniteField # base field (symplectic)
    E::FqNmodFiniteField # additive field
    n::Integer
    k::Union{Integer, Rational{Int64}}
    d::Union{Integer, Missing}
    stabs::fq_nmod_mat
    dualgens::fq_nmod_mat
    logspace::Union{fq_nmod_mat, Missing}
    logicals::Union{Vector{Tuple{fq_nmod_mat, fq_nmod_mat}}, Missing}
    charvec::Vector{Int64} # make bools, uint8's?
    signs::Vector{Int64}
    dualweightenum::Union{WeightEnumerator, Missing}
    logsweightenum::Union{WeightEnumerator, Missing}
    Pauliweightenum::Union{WeightEnumerator, Missing}
end

field(S::AbstractQuantumCode) = S.F
quadraticfield(S::AbstractQuantumCode) = S.E
length(S::AbstractQuantumCode) = S.n
numqubits(S::AbstractQuantumCode) = S.n
dimension(S::AbstractQuantumCode) = S.k
signs(S::AbstractQuantumCode) = S.signs
Xsigns(S::CSSCode) = S.Xsigns
Zsigns(S::CSSCode) = S.Zsigns
stabilizers(S::AbstractQuantumCode) = S.stabs
Xstabilizers(S::CSSCode) = S.Xstabs
Zstabilizers(S::CSSCode) = S.Zstabs
numXstabs(S::CSSCode) = size(S.Xstabs, 1)
numZstabs(S::CSSCode) = size(S.Zstabs, 1)
normalizermatrix(S::AbstractQuantumCode) = S.dualgens
charactervector(S::AbstractQuantumCode) = S.charvec
minimumdistanceX(S::CSSCode) = S.dx
minimumdistanceZ(S::CSSCode) = S.dz
minimumdistance(S::CSSCode) = S.d

function quadratictosymplectic(M::fq_nmod_mat)
    E = base_ring(M)
    iseven(degree(E)) || error("The base ring of the given matrix is not a quadratic extension.")
    F, _ = FiniteField(Int64(characteristic(E)), div(degree(E), 2), "ω")
    nrows = size(M, 1)
    ncols = size(M, 2)
    Msym = zero_matrix(F, nrows, 2 * ncols)
    for r in 1:nrows
        for c in 1:ncols
            if !iszero(M[r, c])
                Msym[r, c] = F(coeff(M[r, c], 0))
                Msym[r, c + ncols] = F(coeff(M[r, c], 1))
            end
        end
    end
    return Msym
end

function symplectictoquadratic(M::fq_nmod_mat)
    iseven(size(M, 2)) || error("Input to symplectictoquadratic is not of even length.")
    nrows = size(M, 1)
    ncols = div(size(M, 2), 2)
    F = base_ring(M)
    E, ω = FiniteField(Int64(characteristic(F)), 2 * degree(F), "ω")
    ϕ = embed(F, E)
    Mquad = zero_matrix(E, nrows, ncols)
    for r in 1:nrows
        for c in 1:ncols
            Mquad[r, c] = ϕ(M[r, c]) + ϕ(M[r, c + ncols]) * ω
        end
    end
    return Mquad
end

function _Paulistringtosymplectic(str::T) where T <: Union{String, Vector{Char}}
    n = length(str)
    F, _ = FiniteField(2, 1, "ω")
    sym = zero_matrix(F, 1, 2 * n)
    for (i, c) in enumerate(str)
        if c == 'X'
            sym[1, i] = F(1)
        elseif c == 'Z'
            sym[1, i + n] = F(1)
        elseif c == 'Y'
            sym[1, i] = 1
            sym[1, i + n] = F(1)
        elseif c != 'I'
            error("Encountered non-{I, X, Y, Z} character in Pauli string. This function is only defined for binary strings.")
        end
    end
    return sym
end
_Paulistringtosymplectic(A::Vector{T}) where T <: Union{String, Vector{Char}} = vcat([_Paulistringtosymplectic(s) for s in A]...)
_Paulistringstofield(str::T) where T <: Union{String, Vector{Char}} = symplectictoquadratic(_Paulistringtosymplectic(str))
_Paulistringstofield(A::Vector{T}) where T <: Union{String, Vector{Char}} = vcat([_Paulistringstofield(s) for s in A]...)

# need symplectictoPaulistring
# quadratictoPaulistring

symplecticstabilizers(S::AbstractQuantumCode) = quadratictosymplectic(S)

function _getsigns(A::fq_nmod_mat, charvec::Vector{Int64})
    if length(charvec) == 2 * size(A, 2)
        A = quadratictosymplectic(A)
    end
    length(charvec) == size(A, 2) || error("Input to _getsigns is expected to be in symplectic form and of the same length as the characteristic vector.")
    if charvec == ones(size(A, 2))
        return ones(Int64, div(size(A, 2), 2))
    end

    signs = Vector{Int64}()
    for r in 1:size(A, 1)
        parity = 1
        for c = 1:size(A, 2)
            if !iszero(A[r, c]) && charvec[c] == -1
                parity *= -1
            end
        end
        append!(signs, parity)
    end
    return signs
end

function _processstrings(SPauli::Vector{T}, charvec::Union{Vector{Int64}, Vector{Any}}=[]) where T <: Union{String, Vector{Char}}
    Paulisigns = Vector{Int64}()
    SPaulistripped = Vector{String}()
    for (i, s) in enumerate(SPauli)
        if s[1] ∈ ['I', 'X', 'Y', 'Z']
            append!(Paulisigns, 1)
            push!(SPaulistripped, s)
        elseif s[1] == '+'
            append!(Paulisigns, 1)
            push!(SPaulistripped, s[2:end])
        elseif s[1] == '-'
            append!(Paulisigns, -1)
            push!(SPaulistripped, s[2:end])
        else
            error("The first element of Pauli string $i is neither a Pauli character or +/-: $s.")
        end
    end

    for s in SPaulistripped
        for i in s
            if i ∉ ['I', 'X', 'Y', 'Z']
                error("Element of provided Pauli string is not a Pauli character: $s.")
            end
        end
    end

    n = length(SPaulistripped[1])
    for s in SPaulistripped
        length(s) == n || error("Provided Pauli strings are not all of the same length.")
    end

    if isempty(charvec)
        charvec = ones(Int64, 2 * n)
    else
        if length(charvec) != 2 * n
            error("Unexpected length of character vector. Expected: $(2 * n), received: $(length(charvec)).")
        end

        for s in charvec
            (s == 1 || s == -1) || error("Qubit phases must be +/- 1; received: $s")
        end
    end
    return SPaulistripped, charvec
end

function splitsymplecticstabilizers(S::fq_nmod_mat, signs::Vector{Int64})
    Xstabs = []
    Xsigns = Vector{Int64}()
    Zstabs = []
    Zsigns = Vector{Int64}()
    mixedstabs = []
    mixedsigns = Vector{Int64}()

    half = div(size(S, 2), 2)
    for r in 1:size(S, 1)
        # use views?
        s = S[r, :]
        if iszero(s)
            continue
        else
            sx = iszero(s[1, 1:half])
            sz = iszero(s[1, half + 1:end])
            if (sx && !sz)
                push!(Zstabs, s)
                append!(Zsigns, signs[r])
            elseif !sx && sz
                push!(Xstabs, s)
                append!(Xsigns, signs[r])
            elseif !sx && !sz
                push!(mixedstabs, s)
                append!(mixedsigns, signs[r])
            end
        end
    end

    if !isempty(Xstabs)
        Xstabs = vcat(Xstabs...)
    end
    if !isempty(Zstabs)
        Zstabs = vcat(Zstabs...)
    end
    if !isempty(mixedstabs)
        mixedstabs = vcat(mixedstabs...)
    end
    return Xstabs, Xsigns, Zstabs, Zsigns, mixedstabs, mixedsigns
end

function isCSSsymplectic(S::fq_nmod_mat, signs::Vector{Int64}=[], trim::Bool=true)
    Xstabs, Xsigns, Zstabs, Zsigns, mixedstabs, mixedsigns = splitsymplecticstabilizers(S, signs)
    if isempty(mixedstabs)
        if trim
            half = div(size(Xstabs, 2), 2)
            return true, Xstabs[:, 1:half], Xsigns, Zstabs[:, half + 1:end], Zsigns
        else
            return true, Xstabs, Xsigns, Zstabs, Zsigns
        end
    else
        if trim
            half = div(size(Xstabs, 2), 2)
            return false, Xstabs[:, 1:half], Xsigns, Zstabs[:, half + 1:end], Zsigns, mixedstabs, mixedsigns
        else
            return false, Xstabs, Xsigns, Zstabs, Zsigns, mixedstabs, mixedsigns
        end
    end
end

function CSSCode(C1::AbstractLinearCode, C2::AbstractLinearCode, charvec::Union{Vector{Int64}, Vector{Any}}=[])
    length(C1) ==  length(C2) || error("Both codes must have the same length in the CSS construction.")
    field(C1) == field(C2) || error("Both codes must be over the same base field in the CSS construction.")
    C2 ⊆ C1 || error("The second argument must be a subset of the first in the CSS construction.")
    if !isempty(charvec)
        for s in charvec
            (s == 1 || s == -1) || error("X-stabilizer signs must be +/- 1; received: $s")
        end
    end

    # C2 ⊆ C1
    # k = k1 - k2
    # d >= minimum(d1, d2^⟂)
    # X - H(C2^⟂), Z - H(C1)
    D2 = dual(C2)
    S = directsum(paritycheckmatrix(D2), paritycheckmatrix(C1))
    Sq2 = symplectictoquadratic(S)
    if isempty(charvec)
        charvec = ones(Int64, 2 * length(C1))
        signs = ones(Int64, size(S, 1))
        Xsigns = ones(Int64, size(paritycheckmatrix(D2), 1))
        Zsigns = ones(Int64, size(paritycheckmatrix(C1), 1))
    else
        signs = _getsigns(S, charvec)
        Xsigns = signs[1:size(paritycheckmatrix(D2), 1), :]
        Zsigns = signs[size(paritycheckmatrix(D2), 1) + 1:end, :]
    end
    dualgens = directsum(generatormatrix(C1), generatormatrix(D2))

    return CSSCode(field(C1), base_ring(Sq2), length(C1), dimension(C1) - dimension(C2),
        missing, minimumdistance(D2), minimumdistance(C1), Sq2, paritycheckmatrix(D2),
        paritycheckmatrix(C1), C2, C1, signs, Xsigns, Zsigns, dualgens, missing,
        missing, charvec, missing, missing, missing)
end

function CSSCode(C::LinearCode, charvec::Union{Vector{Int64}, Vector{Any}}=[])
    # this should have Xstabs = Zstabs
    D = dual(C)
    C ⊆ D || error("The single code CSS construction requires C ⊆ C^⟂.")
    if !isempty(charvec)
        for s in charvec
            (s == 1 || s == -1) || error("X-stabilizer signs must be +/- 1; received: $s")
        end
    end

    # C2 ⊆ C1
    # k = k1 - k2
    # d >= minimum(d1, d2^⟂)
    # X - H(C2^⟂), Z - H(C1)
    S = directsum(paritycheckmatrix(D), paritycheckmatrix(D))
    Sq2 = symplectictoquadratic(S)
    k = dimension(D) - dimension(C)
    if isempty(charvec)
        charvec = ones(Int64, 2 * length(C))
        signs = ones(Int64, size(S, 1))
        Xsigns = ones(Int64, size(paritycheckmatrix(D), 1))
        Zsigns = ones(Int64, size(paritycheckmatrix(D), 1))
    else
        signs = _getsigns(S, charvec)
        Xsigns = signs[1:size(paritycheckmatrix(D), 1), :]
        Zsigns = signs[size(paritycheckmatrix(D), 1) + 1:end, :]
    end
    dualgens = directsum(generatormatrix(D), generatormatrix(D))

    return CSSCode(field(C), base_ring(Sq2), length(C), k, missing, minimumdistance(D),
        minimumdistance(D), Sq2, paritycheckmatrix(D), paritycheckmatrix(D), C, D,
        signs, Xsigns, Zsigns, dualgens, missing, missing, charvec, missing,
        missing, missing)
end

function CSSCode(Xmatrix::fq_nmod_mat, Zmatrix::fq_nmod_mat, charvec::Union{Vector{Int64}, Vector{Any}}=[])
    # this should have Zstabs * Xstabs' = 0
    # set dual containing if Xstabs == Zstabs, else not
    size(Xmatrix, 2) ==  size(Zmatrix, 2) || error("Both matrices must have the same length in the CSS construction.")
    base_ring(Xmatrix) == base_ring(Zmatrix) || error("Both matrices must be over the same base field in the CSS construction.")
    # TODO can remove rank check now that we have full character vector
    rank(Xmatrix) == size(Xmatrix, 1) || error("Provided X-stabilizer matrix is not full rank.")
    rank(Zmatrix) == size(Zmatrix, 1) || error("Provided Z-stabilizer matrix is not full rank.")
    if !isempty(charvec)
        for s in charvec
            (s == 1 || s == -1) || error("X-stabilizer signs must be +/- 1; received: $s")
        end
    end

    S = directsum(Xmatrix, Zmatrix)
    aresymplecticorthogonal(S, S) || error("The given matrices are not symplectic orthogonal.")
    Sq2 = symplectictoquadratic(S)
    if isempty(charvec)
        charvec = ones(Int64, size(S, 2))
        signs = ones(Int64, size(S, 1))
        Xsigns = ones(Int64, size(Xmatrix, 1))
        Zsigns = ones(Int64, size(Zmatrix, 1))
    else
        signs = _getsigns(S, charvec)
        Xsigns = signs[1:size(Xmatrix, 1), :]
        Zsigns = signs[size(Xmatrix, 1) + 1:end, :]
    end

    G = hcat(S[:, size(Xmatrix, 2) + 1:end], -S[:, 1:size(Xmatrix, 2)])
    _, H = right_kernel(G)
    size(H', 1) == 2 * size(Sq2, 2) - size(Sq2, 1) || error("Normalizer matrix is not size n + k.")
    # note the H here is transpose of the standard definition
    !(!iszero(G * H) || !iszero(H' * G')) || error("Normalizer matrix is not transpose, symplectic orthogonal.")
    dualgens = symplectictoquadratic(H')

    return CSSCode(base_ring(Xmatrix), base_ring(Sq2), size(Xmatrix, 2),
        size(Xmatrix, 2) - size(Xmatrix, 1) - size(Zmatrix, 1),
        missing, missing, missing, Sq2, Xmatrix, Zmatrix, missing, missing, signs,
        Xsigns, Zsigns, dualgens, missing, missing, charvec, missing, missing,
        missing)
end

function CSSCode(SPauli::Vector{T}, charvec::Union{Vector{Int64}, Vector{Any}}=[]) where T <: Union{String, Vector{Char}}
    SPaulistripped, charvec = _processstrings(SPauli, charvec)
    S = _Paulistringtosymplectic(SPaulistripped)
    aresymplecticorthogonal(S, S) || error("The given stabilizers are not symplectic orthogonal.")
    if isempty(charvec)
        charvec = ones(Int64, size(S, 2))
    end
    signs = _getsigns(S, charvec)
    Sq2 = symplectictoquadratic(S)

    del = Vector{Int64}()
    for r in 1:size(Sq2, 1)
        if iszero(Sq2[r, :])
            append!(del, r)
        end
    end

    # not necessary but faster to skip
    if !isempty(del)
        Sq2 = Sq2[setdiff(1:size(Sq2, 1), del), :]
        S = S[setdiff(1:size(Sq2, 1), del), :]
        signs = signs[setdiff(1:size(Sq2, 1), del)]
    end

    ###############
    # need to figure out how to track signs as it drops rank
    # would be easier to reassign given character vector instead of stabilizer signs
    ###############

    # if rank(Sq2) != size(Sq2, 1)
    # rk = rank(G)
    # !iszero(rk) || error("Rank zero matrix passed into LinearCode constructor.")
    # if rk < size(G, 1)
    #     _, G = rref(G)
    # end

    # find dual
    G = hcat(S[:, size(Sq2, 2) + 1:end], -S[:, 1:size(Sq2, 2)])
    _, H = right_kernel(G)
    size(H', 1) == 2 * size(Sq2, 2) - size(Sq2, 1) || error("Normalizer matrix is not size n + k.")
    # note the H here is transpose of the standard definition
    !(!iszero(G * H) || !iszero(H' * G')) || error("Normalizer matrix is not transpose, symplectic orthogonal.")
    dualgens = symplectictoquadratic(H')

    args = isCSSsymplectic(S, signs, true)
    if args[1]
        return CSSCode(base_ring(S), base_ring(Sq2), size(Sq2, 2), size(Sq2, 2) - size(Sq2, 1),
            missing, missing, missing, Sq2, args[2], args[4], missing, missing, signs,
            args[3], args[5], dualgens, missing, missing, charvec, missing, missing,
            missing)
    else
        error("Provided Pauli strings are not CSS.")
    end
end

# entanglement-assisted is not symplectic orthogonal
function QuantumCode(SPauli::Vector{T}, charvec::Union{Vector{Int64}, Vector{Any}}=[]) where T <: Union{String, Vector{Char}}
    SPaulistripped, charvec = _processstrings(SPauli, charvec)
    S = _Paulistringtosymplectic(SPaulistripped)
    aresymplecticorthogonal(S, S) || error("The given stabilizers are not symplectic orthogonal.")
    if isempty(charvec)
        charvec = ones(Int64, size(S, 2))
    end
    signs = _getsigns(S, charvec)
    Sq2 = symplectictoquadratic(S)

    # make del zero rows, cols helper functions
    # which will return a flag if this was ever done
    del = Vector{Int64}()
    for r in 1:size(Sq2, 1)
        if iszero(Sq2[r, :])
            append!(del, r)
        end
    end

    # not necessary but faster to skip
    if !isempty(del)
        Sq2 = Sq2[setdiff(1:size(Sq2, 1), del), :]
        S = S[setdiff(1:size(Sq2, 1), del), :]
        signs = signs[setdiff(1:size(Sq2, 1), del)]
    end

    # need to code an additive only rref but then apply _getsigns
    # if rank(Sq2) != size(Sq2, 1)
    # rk = rank(G)
    # !iszero(rk) || error("Rank zero matrix passed into LinearCode constructor.")
    # if rk < size(G, 1)
    #     _, G = rref(G)
    # end

    # find dual
    G = hcat(S[:, size(Sq2, 2) + 1:end], -S[:, 1:size(Sq2, 2)])
    _, H = right_kernel(G)
    size(H', 1) == 2 * size(Sq2, 2) - size(Sq2, 1) || error("Normalizer matrix is not size n + k.")
    # note the H here is transpose of the standard definition
    !(!iszero(G * H) || !iszero(H' * G')) || error("Normalizer matrix is not transpose, symplectic orthogonal.")
    dualgens = symplectictoquadratic(H')

    # q^n / p^k but rows is n - k
    dimcode = Int64(order(base_ring(S)))^size(Sq2, 2) // Int64(characteristic(base_ring(S)))^(size(Sq2, 1))
    if isinteger(dimcode)
        dimcode = Int64(log(Int64(characteristic(base_ring(S))), Int64(dimcode)))
    end

    args = isCSSsymplectic(S, signs, true)
    if args[1]
        return CSSCode(base_ring(S), base_ring(Sq2), size(Sq2, 2), dimcode,
            missing, missing, missing, Sq2, args[2], args[4], missing, missing, signs,
            args[3], args[5], dualgens, missing, missing, charvec, missing, missing,
            missing)
    else
        return QuantumCode(base_ring(S), base_ring(Sq2), size(Sq2, 2), dimcode,
            missing, Sq2, dualgens, missing, missing, charvec, signs, missing,
            missing, missing)
    end
end

function QuantumCode(Sq2::fq_nmod_mat, symp::Bool=false, charvec::Union{Vector{Int64}, Vector{Any}}=[]) where T <: Union{String, Vector{Char}}
    if symp
        iseven(size(Sq2, 2)) || error("Expected a symplectic input but the input matrix has an odd number of columns.")
        if !isempty(charvec)
            size(Sq2, 2) == length(charvec) || error("The characteristic value is of incorrect length.")
        end
        S = Sq2
        Sq2 = symplectictoquadratic(Sq2)
    else
        E = base_ring(Sq2)
        iseven(degree(E)) || error("The base ring of the given matrix is not a quadratic extension.")
        S = quadratictosymplectic(Sq2)
    end
    if isempty(charvec)
        charvec = ones(Int64, size(S, 2))
    end
    aresymplecticorthogonal(S, S) || error("The given stabilizers are not symplectic orthogonal.")
    signs = _getsigns(S, charvec)
    # make del zero rows, cols helper functions
    # which will return a flag if this was ever done
    del = Vector{Int64}()
    for r in 1:size(Sq2, 1)
        if iszero(Sq2[r, :])
            append!(del, r)
        end
    end

    # not necessary but faster to skip
    if !isempty(del)
        Sq2 = Sq2[setdiff(1:size(Sq2, 1), del), :]
        S = S[setdiff(1:size(Sq2, 1), del), :]
        signs = signs[setdiff(1:size(Sq2, 1), del)]
    end

    # need to code an additive only rref but then apply _getsigns
    # if rank(Sq2) != size(Sq2, 1)
    # rk = rank(G)
    # !iszero(rk) || error("Rank zero matrix passed into LinearCode constructor.")
    # if rk < size(G, 1)
    #     _, G = rref(G)
    # end

    # find dual
    G = hcat(S[:, size(Sq2, 2) + 1:end], -S[:, 1:size(Sq2, 2)])
    _, H = right_kernel(G)
    size(H', 1) == 2 * size(Sq2, 2) - size(Sq2, 1) || error("Normalizer matrix is not size n + k.")
    # note the H here is transpose of the standard definition
    !(!iszero(G * H) || !iszero(H' * G')) || error("Normalizer matrix is not transpose, symplectic orthogonal.")
    dualgens = symplectictoquadratic(H')

    try
        temp = div(Int64(order(base_ring(Sq2))), Int64(characteristic(base_ring(Sq2))))
        dimcode = temp^(size(Sq2, 2) - size(Sq2, 1))
    catch
        # over flows hard here - make rational{BigInt}?
        dimcode = Int64(order(base_ring(Sq2)))^size(Sq2, 2) // Int64(characteristic(base_ring(Sq2)))^(size(Sq2, 1))
    end
    if isinteger(dimcode)
        dimcode = Int64(log(Int64(characteristic(base_ring(S))), Int64(dimcode)))
    end

    args = isCSSsymplectic(S, signs, true)
    if args[1]
        return CSSCode(base_ring(S), base_ring(Sq2), size(Sq2, 2), dimcode,
            missing, missing, missing, Sq2, args[2], args[4], missing, missing, signs,
            args[3], args[5], dualgens, missing, missing, charvec, missing, missing,
            missing)
    else
        return QuantumCode(base_ring(S), base_ring(Sq2), size(Sq2, 2), dimcode,
            missing, Sq2, dualgens, missing, missing, charvec, signs, missing,
            missing, missing)
    end
end

# this returns a basis for the space S^⟂ / S
# work would need to be done to extract appropriate X_i, Z_i pairs from this
# can implement Gottesman algorithm later
function logicalspace(S::AbstractQuantumCode, alg::String="quotient")
    if !ismissing(S.logspace)
        return S.logspace
    end

    alg ∈ ["quotient"] || error("Algorithm '$alg' not specified in logicalspace.")
    if alg == "quotient"
        F = quadraticfield(S)
        G = stabilizers(S)
        Gdual = normalizermatrix(S)
        V = VectorSpace(F, length(S))
        U, UtoV = sub(V, [V(G[i, :]) for i in 1:size(G, 1)])
        W, WtoV = sub(V, [V(Gdual[i, :]) for i in 1:size(Gdual, 1)])
        gensofUinW = [preimage(WtoV, UtoV(g)) for g in gens(U)]
        UinW, UinWtoW = sub(W, gensofUinW)
        Q, WtoQ = quo(W, UinW)
        C2modC1basis = [WtoV(x) for x in [preimage(WtoQ, g) for g in gens(Q)]]
        Fbasis = [[F(C2modC1basis[j][i]) for i in 1:AbstractAlgebra.dim(parent(C2modC1basis[1]))] for j in 1:length(C2modC1basis)]
        G2 = matrix(F, length(Fbasis), length(Fbasis[1]), vcat(Fbasis...))

        # returns linear object, make additive if necessary
        if size(G2, 1) == dimension(S)
            S.logspace = vcat(G2, gen(F) * G2)
            return S.logspace
        elseif size(G2, 1) == 2 * dimension(S)
            S.logspace = G2
            return S.logspace
        else
            error("Logical space produced of incorrect dimension; expected: ", 2 * dimension(S), ", received: ", size(G2, 1))
        end
    end
end

# need a canonical phase to set Weyl pairs to
# does not check which elements in Weyl pair are X or Z
function logicaloperators(Q::AbstractQuantumCode)
    logspace = logicalspace(Q)
    Weylpairs = Vector{Tuple{fq_nmod_mat, fq_nmod_mat}}()
    todo = [quadratictosymplectic(logspace)[i, :] for i in 1:size(logspace, 1)]

    count = 1
    while !isempty(todo)
        curr1 = popfirst!(todo)
        commutes = Vector{fq_nmod_mat}()
        doesnt = Vector{Tuple{fq_nmod_mat, fq_nmod}}()
        for i in 1:length(todo)
            SIP = symplecticinnerproduct(curr1, todo[i])
            if iszero(SIP)
                push!(commutes, todo[i])
            else
                push!(doesnt, (todo[i], SIP))
            end
        end

        if !isempty(doesnt)
            curr2 = popfirst!(doesnt)
            curr2, curr2SIP = curr2[1], curr2[2]
        else
            error("Logical basis element anti-commutes with no elements.")
        end

        todo = Vector{fq_nmod_mat}()
        for i in 1:length(doesnt)
            elm, SIP = doesnt[i]
            elm .-= SIP * inv(curr2SIP) * curr2
            push!(todo, elm)
        end

        push!(Weylpairs, [symplectictoquadratic(curr1), symplectictoquadratic(curr2)])
        count < 50 ||error("Logical operator loop has gone through 50 iterations. Increase if more than this many logical qudits.")
        count += 1
    end

    length(Weylpairs) == dimension(Q) || error("Completed symplectic Gram-Schmidt without generating all logicals.")
    # isempty(todo) || error("Completed symplectic Gram-Schmidt without using all basis elements.")

    Q.logicals = Weylpairs
    return Weylpairs
end

# can't check rank here because a self-dual code will have half rank since function is not additive
# need to think about what this means for an irrational dimension
# should check commutation relations
function setlogicals!(Q::AbstractQuantumCode, L::fq_nmod_mat)
    size(L) == (2 * dimension(Q), length(Q)) || error("Provided matrix is of incorrect size for the logical space.")
    # rank(L) == size(L, 1) || error("Provided matrix is not of full rank.")
    # aresymplecticorthogonal(stabilizers(Q), L) || error("Provided matrix does not commute with the code.")
    # !aresymplecticorthogonal(L, L) || error("Provided matrix should not be symplectic self-orthogonal.")
    Q.logspace = L
end

# update for roots of unity
function changesigns!(Q::AbstractQuantumCode, charvec::Vector{Int64})
    length(charvec) == 2 * length(Q) || error("Characteristic vector is of improper length for the code.")
    for s in charvec
        (s == 1 || s == -1) || error("Qubit phases must be +/- 1; received: $s")
    end
    Q.signs = _getsigns(stabilizers(Q), charvec)
    Q.charvec = charvec
end

# add for minimum distance check
function show(io::IO, S::AbstractQuantumCode)
    if get(io, :compact, false)
        if typeof(S) <: CSSCode
            if typeof(dimension(S)) <: Integer
                println(io, "[[$(length(S)), $(dimension(S))]]_$(order(field(S))) CSS code.")
            else # don't think this can ever be reached
                println(io, "(($(length(S)), $(dimension(S))))_$(order(field(S))) CSS code.")
            end
        else
            if typeof(dimension(S)) <: Integer
                println(io, "[[$(length(S)), $(dimension(S))]]_$(order(field(S))) quantum (additive) code.")
            else
                println(io, "(($(length(S)), $(dimension(S))))_$(order(field(S))) quantum (additive) code.")
            end
        end
    else
        if typeof(S) <: CSSCode
            if typeof(dimension(S)) <: Integer
                println(io, "[[$(length(S)), $(dimension(S))]]_$(order(field(S))) CSS code.")
            else # don't think this can ever be reached
                println(io, "(($(length(S)), $(dimension(S))))_$(order(field(S))) CSS code.")
            end
            println(io, "X-stabilizer matrix: $(numXstabs(S)) × $(length(S))")
            for i in 1:numXstabs(S)
                print(io, "\t")
                for j in 1:length(S)
                    if j != length(S)
                        print(io, "$(S.Xstabs[i, j]) ")
                    elseif j == length(S) && i != length(S) - dimension(S)
                        println(io, "$(S.Xstabs[i, j])")
                    else
                        print(io, "$(S.Xstabs[i, j])")
                    end
                end
            end
            println(io, "Z-stabilizer matrix: $(numZstabs(S)) × $(length(S))")
            for i in 1:numZstabs(S)
                print(io, "\t")
                for j in 1:length(S)
                    if j != length(S)
                        print(io, "$(S.Zstabs[i, j]) ")
                    elseif j == length(S) && i != length(S) - dimension(S)
                        println(io, "$(S.Zstabs[i, j])")
                    else
                        print(io, "$(S.Zstabs[i, j])")
                    end
                end
            end
        else
            if typeof(dimension(S)) <: Integer
                println(io, "[[$(length(S)), $(dimension(S))]]_$(order(field(S))) quantum (additive) code.")
            else
                println(io, "(($(length(S)), $(dimension(S))))_$(order(field(S))) quantum (additive) code.")
            end
            println(io, "Stabilizer matrix: $(size(S.stabs, 1)) × $(length(S))")
            for i in 1:size(S.stabs, 1)
                print(io, "\t")
                for j in 1:length(S)
                    if j != length(S)
                        print(io, "$(S.stabs[i, j]) ")
                    elseif j == length(S) && i != length(S) - dimension(S)
                        println(io, "$(S.stabs[i, j])")
                    else
                        print(io, "$(S.stabs[i, j])")
                    end
                end
            end
        end
    end
end

function Xsyndrome(S::CSSCode, v::fq_nmod_mat)
    n = length(S)
    if length(v) == 2 * n
        v = v[n + 1:end]
    end
    !(size(v) != (n, 1) && size(v) != (1, n)) ||
        error("Vector to be tested is of incorrect dimension; expected length $n, received: $(size(v)).")
    base_ring(v) == field(S) || error("Vector must have the same base ring as the stabilizers.")

    size(v, 1) != 1 || return Xstabilizers(S) * v'
    return Xstabilizers(S) * v
end

function Zsyndrome(S::CSSCode, v::fq_nmod_mat)
    n = length(S)
    if length(v) == 2 * n
        v = v[1:n]
    end
    !(size(v) != (n, 1) && size(v) != (1, n)) ||
        error("Vector to be tested is of incorrect dimension; expected length $n, received: $(size(v)).")
    base_ring(v) == field(S) || error("Vector must have the same base ring as the stabilizers.")

    size(v, 1) != 1 || return Zstabilizers(S) * v'
    return Zstabilizers(S) * v
end

function syndrome(S::AbstractQuantumCode, v::fq_nmod_mat)
    n = length(S)
    !(size(v) != (2 * n, 1) && size(v) != (1, 2 * n)) ||
        error("Vector to be tested is of incorrect dimension; expected length $(2 * n), received: $(size(v)).")
    # base_ring(v) == field(S) || error("Vector must have the same base ring as the stabilizers.")

    size(v, 1) != 1 || return symplecticstabilizers(S) * v'
    return symplecticstabilizers(S) * v
end

# function Singletonbound()
#     n - k ≧ 2(d - 1)
# end

# will have to get signs on them separately
# could make this considerably faster by using combinatorics package and skipping 0's
# using Combinatorics
# iter = combinations(1:size(stabs, 2))
# but this only in base 2
function allstabilizers(Q::AbstractQuantumCode)
    E = quadraticfield(Q)
    all = Vector{fq_nmod_mat}()
    stabs = stabilizers(Q)
    for iter in Base.Iterators.product([0:(Int64(characteristic(field(Q))) - 1) for _ in 1:size(stabs, 1)]...)
        stab = E(iter[1]) * stabs[1, :]
        for r in 2:size(stabs, 1)
            if !iszero(iter[r])
                stab += E(iter[r]) * stabs[r, :]
            end
        end
        push!(all, stab)
    end
    return all
end

#############################
#   Generator Coefficients  #
#############################

# # to run this, need an error model,
# # need wrtV,
# # need clarification on general stabilizer codes
# # need to process Pauli here
# # sum of syndrome and logical?
# function generatorcoefficients(Q::AbstractQuantumCode, θ::Union{Float64, Missing},
#     Paui::Char=' ')
#
#     synlen = size(stabilizers(Q), 1)
#     numsyns = BigInt(characteristic(field(Q)))^synlen
#     # do some checks here for feasibility
#     W, T = Pauliweightenumerator(Q, Pauli, true, true)
#     prevsyn = zeros(Int64, synlen)
#     logs = logicalspace(Q)
#
#     for i in 1:size(logs, 1)
#         γ = logs[i, :]
#         for j in 1:numsyns
#             μ = digits(j, base=2, pad=synlen)
#             shifttrellis!(T, μ .+ γ, wrtV, err_model, chractervector(Q), Pauli)
#             push!(W, _Pauliweightenumerator(T))
#             _reducepoly(W)
#         end
#     end
#
#     if ismissing(θ)
#         # @vars θ
#         n = length(Q)
#         sum = Complex(0.0) # typeof(sum) = Sym
#         for term in W
#             if term[2] < 0 || term[3] < 0 || term[4] < 0
#                 sum -= term[1] * cos(θ / 2)^abs(term[2] + term[3]) * (1im * sin(θ / 2))^abs(term[3] + term[4])
#             else
#                 sum += term[1] * cos(θ / 2)^abs(term[2] + term[3]) * (1im * sin(θ / 2))^abs(term[3] + term[4])
#             end
#         end
#         return sum, W
#     else
#         n = length(Q)
#         sum = Complex(0.0)
#         for term in W
#             if term[2] < 0 || term[3] < 0 || term[4] < 0
#                 sum -= term[1] * cos(θ / 2)^(n - abs(term[2] + term[3])) * (1im * sin(θ / 2))^abs(term[3] + term[4])
#             else
#                 sum += term[1] * cos(θ / 2)^(n - abs(term[2] + term[3])) * (1im * sin(θ / 2))^abs(term[3] + term[4])
#             end
#         end
#         return sum, W
#     end
# end

# function plotgeneratorcoefficients(W, θmin, Qmax)
#
# end
#
# function plotgeneratorcoefficients(Q, θmin, Qmax)
#
# end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
#=

data initializer for tribe model

Ignacio Quintero Mächler

t(-_-t)

June 14 2017

=#

"""
    read_data(tree_file::String, data_file::String; delim::Char = '\t', eol::Char = '\r')

Read a phylogenetic tree using **ape** package in R through
`RCall` and the data file with the trait and biogeographic information.
"""
function read_data_tribe(tree_file::String,
                         data_file::String)

  tree, bts = read_tree(tree_file)

  tip_labels = Dict(i => val for (val,i) = enumerate(tree.tlab))

  data = readdlm(data_file)

  if size(data,1) != (tree.nnod + 1)
    data = readdlm(data_file, '\t', '\r')
  end

  if size(data,1) != (tree.nnod + 1)
    data = readdlm(data_file, '\t', '\n')
  end

  if size(data,1) != (tree.nnod + 1)
    error("Data file cannot be made of the right dimensions.\n Make sure the data file has the same number of rows as tips in the tree")
  end

  data_tlab   = convert(Array{String,1}, data[:,1])
  data_values = convert(Array{Float64,1},data[:,2])
  data_areas  = convert(Array{Int64,2},  data[:,3:end])

  # create dictionaries
  tip_areas = Dict(tip_labels[val] => data_areas[i,:]
                   for (i,val) = enumerate(data_tlab))

  tip_values = Dict(tip_labels[val] => data_values[i]
                    for (i,val) = enumerate(data_tlab))

  return tip_values, tip_areas, tree, bts
end

"""
    initialize_data(tip_values::Dict{Int64,Float64}, tip_areas ::Dict{Int64,Array{Int64,1}}, m::Int64, tree::rtree, bts::Array{Float64,1})

Function to initialize `X` and `Y` matrix given
the tip_values and tip_areas (as Dictionaries).
"""
function initialize_data(tip_values::Dict{Int64,Float64},
                         tip_areas ::Dict{Int64,Array{Int64,1}},
                         min_dt    ::Float64,
                         tree      ::rtree,
                         bts       ::Array{Float64,1})

  br     = branching_times(tree)
  n      = tree.nnod + 1
  nareas = length(tip_areas[1])

  #*
  # make times and δt vector
  #*

  # make sure each branch has nareas + 1 sampling points
  ets = unique(br[:,4])
  for i = sortperm(br[:,3])
    # number of times that cross the branch
    nover = length(findall(x -> br[i,4] > x > br[i,5], ets))
    if nover < nareas
      nets = convert(Array{Float64,1},
        range(br[i,4], stop = br[i,5], length = nareas-nover+2))
      if length(nets) > 2
        append!(ets,nets[2:(end-1)])
      end
    end
  end

  # epoch times
  sort!(ets, rev = true)

  tr_height = ets[1]
  mδt       = min_dt*tr_height

  # incorporate more 'ets' according to min_dt
  new_ets = Float64[]
  for i in eachindex(ets)

    if i == lastindex(ets)
      if ets[i]/tr_height > min_dt
        append!(new_ets, collect(0:mδt:ets[i])[2:end])
      end
    else
      if (ets[i] - ets[i+1])/tr_height > min_dt
        append!(new_ets, collect(ets[i+1]:mδt:ets[i])[2:end])
      end
    end
  end

  # add new_ets
  append!(ets, new_ets)

  # sort epoch times from start to end
  sort!(ets, rev = true)

  # push present
  push!(ets, 0.0)

  #create δt vector
  δt = abs.(diff(ets))

  # initialize data augmentation matrices
  X = fill(NaN, length(ets), n)
  B = copy(X)
  Y = fill(23, length(ets), n, nareas)

  # which rows are branching points
  wch = indexin(bts, ets)

  # coupled nodes (cells that are coupled in array)
  coup = zeros(Int64, tree.nnod, 3)

  bord = sortperm(br[:,5])[(1:tree.nnod-1) .+ n]

  alive = collect(1:n)

  #which column alive
  wca = 1:n

  wts = sort(wch, rev = true) .+ 1

  setindex!(coup, wts .- 1, :,3)

  wrtf = wts[1]:length(ets)
  X[wrtf, wca] .= 1.0
  B[wrtf, wca]  = repeat(alive, inner = length(wrtf))

  for i = Base.OneTo(tree.nnod-1)
    fn  = br[bord[i],2]
    wda = br[findall(isequal(fn), br[:,1]),2]

    cda = indexin(wda,alive)
    coup[i,1:2] = cda

    alive[cda] = [fn,0]
    wrtf = wts[i+1]:(wts[i]-1)

    B[wrtf,:] = repeat(alive, inner = length(wrtf))

    wca = findall(x -> x > 0, alive)

    X[wts[i+1]:length(ets),wca] .= 1.0
  end

  coup[tree.nnod,1:2] = findall(x -> x > 0, alive)

  ncoup = zeros(Int64,tree.nnod,2)

  for j = Base.OneTo(tree.nnod), i=1:2
    ncoup[j,i] = vecind(coup[j,3], coup[j,i], lastindex(ets))
  end

  X[1,1]      = 1.0
  B[1,1]      = n + 1
  B[B .== 0] .= NaN

  # Brownian bridges initialization for X
  si = initialize_X!(tip_values, X, B, ncoup, δt, tree)

  # declare non-23s for Y
  initialize_Y!(tip_areas, Y, B)

  return X, Y, B, ncoup, δt, tree, si
end

"""
    initialize_X!(tip_values::Dict{Int64,Float64}, X::Array{Float64,2}, B::Array{Float64,2}, ncoup::Array{Int64,2}, δt::Array{Float64,1}, tree::rtree)

Create an initial trait matrix, `X`, using Brownian bridges.
"""
function initialize_X!(tip_values::Dict{Int64,Float64},
                       X         ::Array{Float64,2},
                       B         ::Array{Float64,2},
                       ncoup     ::Array{Int64,2},
                       δt        ::Array{Float64,1},
                       tree      ::rtree)

  co1    = ncoup[:,1]
  co2    = ncoup[:,2]
  nr, nc = size(X)
  nbrs   = 2*nc - 1

  # matrix of delta times
  δtM = zeros(lastindex(δt)+1,nc)
  for i = Base.OneTo(nc)
    δtM[2:end,i] = δt
  end

  Xnod = ncoup[size(ncoup,1):-1:1,:]

  sort!(Xnod, dims=1)

  # brownian motion
  bm_ll = make_bm_ll(tip_values, tree)
  op    = optimize(bm_ll, rand(nc), Optim.Options(
                                    g_tol = 1e-6,
                                    iterations = 100_000,
                                    store_trace = false))

  ar = minimizer(op)[2:end]
  si = minimizer(op)[1]

  X[Xnod[:,1]] = ar

  for i in Base.OneTo(nc)
    X[nr,i] = tip_values[i]
  end

  X[co2] = X[co1]

  for i = setdiff(1:nbrs,nc+1)
    wbranch = findall(isequal(i), B)
    l_i = wbranch[1] - CartesianIndex(1,0)
    l_f = wbranch[end]
    X[CartesianIndices((l_i[1]:l_f[1], l_i[2]:l_i[2]))] =
      bb(X[l_i], X[l_f], δtM[wbranch])
  end

  X[co2] .= NaN

  return si
end

"""
    initialize_Y!(tip_areas::Dict{Int64,Array{Int64,1}}, Y::Array{Int64,3}, B::Array{Float64,2})

Simple function to initialize Y with all 1s (the
real biogeographic sampling is done during the MCMC).
"""
function initialize_Y!(tip_areas::Dict{Int64,Array{Int64,1}},
                       Y        ::Array{Int64,3},
                       B        ::Array{Float64,2})

  # index non 23 for Y
  ind = findall(!isnan, B)
  lB  = length(ind)

  # #D indices
  indY = CartesianIndex{3}[]
  for a = 1:size(Y,3), i = ind
    push!(indY, CartesianIndex(i,a))
  end

  Y[indY] .= 1

  for i in Base.OneTo(size(Y,2))
    Y[end,i,:] = tip_areas[i]
  end

end

"""
Maximum likelihood Brownian Motion.
"""
function make_bm_ll(tip_values::Dict{Int64,Float64},
                    tree      ::rtree)

  wt = tree.ed .<= (tree.nnod + 1)
  tips_ind = findall(wt)

  # base with trait values
  ntr  = zeros(size(tree.ed))

  # assign tip values
  for i = eachindex(tip_values)
    ntr[tips_ind[i]] = tip_values[i]
  end

  # internal nodes
  ins  = unique(tree.ed[:,1])
  lins = length(ins)

  # make triads for all internal nodes
  # including the root
  trios = Array{Int64,1}[]
  ndi  = ins[1]
  daus = findall(isequal(ndi), tree.ed[:,1])
  pushfirst!(daus, 0)
  push!(trios, daus)

  # for all internal nodes
  for i = 2:lins
    ndi  = ins[i]
    daus = findall(isequal(ndi), tree.ed[:,1])
    pushfirst!(daus, findall(isequal(ndi), tree.ed[:,2])[1])
    push!(trios, daus)
  end

  el = tree.el

  function f(p::Array{Float64,1})

    @inbounds begin

      σ² = p[1]

      if σ² <= 0
        return Inf
      end

      for i = eachindex(trios)
        pr, d1, d2 = trios[i]
        if pr == 0
          ntr[d1,1] = ntr[d2,1] = p[i+1]
        else
          ntr[pr,2] = ntr[d1,1] = ntr[d2,1] = p[i+1]
        end
      end

      ll = 0.0
      for i in eachindex(el)
        ll -= logdnorm_tc(ntr[i,2], ntr[i,1], el[i]*σ²)
      end
    end

    return ll
  end

  return f
end

"""
    branching_times(tree::rtree)

Function to estimate absolute
branching times, with time 0 at the
present, time going backward.
"""
function branching_times(tree::rtree)

  @inbounds begin

    n    = tree.nnod + 1
    el_t = findall(x -> x <= n, tree.ed[:,2])

    brs = zeros(lastindex(tree.el),5)

    brs[:, 1:2] = tree.ed
    brs[:, 3]   = tree.el

    for j = eachindex(tree.el)
      if brs[j,1] == (n+1)
        brs[j,4] = 0.0
        brs[j,5] = brs[j,4] + brs[j,3]
      else
        brs[j,4] = brs[brs[j,1] .== brs[:,2],5][1]
        brs[j,5] = brs[j,4] + brs[j,3]
      end
    end

     # change time forward order
    @views tree_depth = brs[n .== brs[:,2],5][1]

    for j = eachindex(tree.el)
      brs[j,4] = tree_depth - brs[j,4]
      brs[j,5] = tree_depth - brs[j,5]
    end

    brs[el_t,5] .= 0.0

  end

  return brs
end

"""
    bb(xs::Float64, xf::Float64, δt::Array{Float64,1})

Brownian bridge simulation function for
a vector of times δt.
"""
function bb(xs::Float64, xf::Float64, δt::Array{Float64,1})

  t  = pushfirst!(cumsum(δt),0.0)
  tl = lastindex(t)
  w  = zeros(tl)

  for i = Base.OneTo(tl-1)
    w[i+1] = randn()*sqrt(δt[i])
  end

  cumsum!(w, w)
  wf = w[tl]
  tf = t[tl]

  return @. xs + w - t/tf * (wf - xf + xs)
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
module Fields

import ..CryptoGroups: bitlength, order, <|, tobits, value, modulus

include("abstract_fields.jl")
include("galois_fields.jl")

export Field, BinaryField, FP, F2PB, F2GNB

end


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

using ValidatedNumerics
using FactCheck

setprecision(Interval, Float64)
setrounding(Interval, :narrow)

facts("Numeric tests") do

    a = @interval(0.1, 1.1)
    b = @interval(0.9, 2.0)
    c = @interval(0.25, 4.0)

    ## Basic arithmetic
    @fact @interval(0.1) --> Interval(9.9999999999999992e-02, 1.0000000000000001e-01)
    @fact +a == a --> true
    @fact a+b --> Interval(9.9999999999999989e-01, 3.1000000000000001e+00)
    @fact -a --> Interval(-1.1000000000000001e+00, -9.9999999999999992e-02)
    @fact a-b --> Interval(-1.9000000000000001e+00, 2.0000000000000018e-01)
    @fact Interval(1//4,1//2) + Interval(2//3) --> Interval(11//12, 7//6)
    @fact Interval(1//4,1//2) - Interval(2//3) --> Interval(-5//12, -1//6)

    @fact 10a --> @interval(10a)
    #@fact 10Interval(1//10) --> one(@interval(1//10))
    @fact Interval(-30.0,-15.0) / Interval(-5.0,-3.0) --> Interval(3.0, 10.0)
    @fact @interval(-30,-15) / @interval(-5,-3) --> Interval(3.0, 10.0)
    @fact b/a --> Interval(8.18181818181818e-01, 2.0000000000000004e+01)
    @fact a/c --> Interval(2.4999999999999998e-02, 4.4000000000000004e+00)
    @fact c/4.0 --> Interval(6.25e-02, 1e+00)
    @fact c/zero(c) --> emptyinterval(c)
    @fact Interval(0.0, 1.0)/Interval(0.0,1.0) --> Interval(0.0, Inf)
    @fact Interval(-1.0, 1.0)/Interval(0.0,1.0) --> entireinterval(c)
    @fact Interval(-1.0, 1.0)/Interval(-1.0,1.0) --> entireinterval(c)
    a = @interval(1.e-20)
    @fact a --> Interval(1.0e-20, 1.0000000000000001e-20)
    @fact diam(a) --> eps(1.e-20)
end

facts("Power tests") do
    @fact @interval(0,3) ^ 2 --> Interval(0, 9)
    @fact @interval(2,3) ^ 2 --> Interval(4, 9)
    @fact @interval(-3,0) ^ 2 --> Interval(0, 9)
    @fact @interval(-3,-2) ^ 2 --> Interval(4, 9)
    @fact @interval(-3,2) ^ 2 --> Interval(0, 9)
    @fact @interval(0,3) ^ 3 --> Interval(0, 27)
    @fact @interval(2,3) ^ 3 --> Interval(8, 27)
    @fact @interval(-3,0) ^ 3 --> @interval(-27., 0.)
    @fact @interval(-3,-2) ^ 3 --> @interval(-27, -8)
    @fact @interval(-3,2) ^ 3 --> @interval(-27., 8.)
    @fact Interval(0,3) ^ -2 --> Interval(1/9, Inf)
    @fact Interval(-3,0) ^ -2 --> Interval(1/9, Inf)
    @fact Interval(-3,2) ^ -2 --> Interval(1/9, Inf)
    @fact Interval(2,3) ^ -2 --> Interval(1/9, 1/4)
    @fact Interval(1,2) ^ -3 --> Interval(1/8, 1.0)
    @fact Interval(0,3) ^ -3 --> @interval(1/27, Inf)
    @fact Interval(-1,2) ^ -3 --> entireinterval()
    @fact_throws ArgumentError Interval(-1, -2) ^ -3  # wrong way round
    @fact Interval(-3,2) ^ (3//1) --> Interval(-27, 8)
    @fact Interval(0.0) ^ 1.1 --> Interval(0, 0)
    @fact Interval(0.0) ^ 0.0 --> emptyinterval()
    @fact Interval(0.0) ^ (1//10) --> Interval(0, 0)
    @fact Interval(0.0) ^ (-1//10) --> emptyinterval()
    @fact ∅ ^ 0 --> ∅
    @fact Interval(2.5)^3 --> Interval(15.625, 15.625)
    #@fact Interval(5//2)^3.0 --> Interval(125//8)

    x = @interval(-3,2)
    @fact x^3 --> @interval(-27, 8)

    @fact @interval(-3,4) ^ 0.5 --> @interval(0, 2)
    @fact @interval(-3,4) ^ 0.5 --> @interval(-3,4)^(1//2)
    @fact @interval(-3,2) ^ @interval(2) --> Interval(0.0, 4.0)
    @fact @interval(-3,4) ^ @interval(0.5) --> Interval(0, 2)
    @fact @biginterval(-3,4) ^ 0.5 --> @biginterval(0, 2)

    @fact @interval(1,27)^@interval(1/3) --> roughly(Interval(1., 3.))
    @fact @interval(1,27)^(1/3) --> roughly(Interval(1., 3.))
    @fact Interval(1., 3.) ⊆ @interval(1,27)^(1//3) --> true
    @fact @interval(0.1,0.7)^(1//3) --> Interval(0.46415888336127786, 0.8879040017426008)
    @fact @interval(0.1,0.7)^(1/3)  --> roughly(Interval(0.46415888336127786, 0.8879040017426008))

    setprecision(Interval, 256)
    x = @biginterval(27)
    y = x^(1//3)
    @fact (0 < diam(y) < 1e-76) --> true
    y = x^(1/3)
    @fact (0 < diam(y) < 1e-76) --> true
    @fact x^(1//3) --> x^(1/3)
end

setprecision(Interval, Float64)

facts("Exp and log tests") do
    @fact exp(@biginterval(1//2)) ⊆ exp(@interval(1//2)) --> true
    @fact exp(@interval(1//2)) --> Interval(1.648721270700128, 1.6487212707001282)
    @fact exp(@biginterval(0.1)) ⊆ exp(@interval(0.1)) --> true
    @fact exp(@interval(0.1)) --> Interval(1.1051709180756475e+00, 1.1051709180756477e+00)
    @fact diam(exp(@interval(0.1))) --> eps(exp(0.1))

    @fact log(@biginterval(1//2)) ⊆ log(@interval(1//2)) --> true
    @fact log(@interval(1//2)) --> Interval(-6.931471805599454e-01, -6.9314718055994529e-01)
    @fact log(@biginterval(0.1)) ⊆ log(@interval(0.1)) --> true
    @fact log(@interval(0.1)) --> Interval(-2.3025850929940459e+00, -2.3025850929940455e+00)
    @fact diam(log(@interval(0.1))) --> eps(log(0.1))

    @fact exp2(@biginterval(1//2)) ⊆ exp2(@interval(1//2)) --> true
    @fact exp2(Interval(1024.0)) --> Interval(1.7976931348623157e308, Inf)
    @fact exp10(@biginterval(1//2)) ⊆ exp10(@interval(1//2)) --> true
    @fact exp10(Interval(308.5)) --> Interval(1.7976931348623157e308, Inf)

    @fact log2(@biginterval(1//2)) ⊆ log2(@interval(1//2)) --> true
    @fact log2(@interval(0.25, 0.5)) --> Interval(-2.0, -1.0)
    @fact log10(@biginterval(1//10)) ⊆ log10(@interval(1//10)) --> true
    @fact log10(@interval(0.01, 0.1)) --> @interval(log10(0.01), log10(0.1))
end

facts("Comparison tests") do
    d = @interval(0.1, 2)

    @fact d < 3 --> true
    @fact d <= 2 --> true
    @fact d < 2 --> false
    @fact -1 < d --> true
    @fact !(d < 0.15) --> true

    # abs
    @fact abs(@interval(0.1, 0.2)) --> Interval(9.9999999999999992e-02, 2.0000000000000001e-01)
    @fact abs(@interval(-1, 2)) --> Interval(0, 2)

    # real
    @fact real(@interval(-1, 1)) --> Interval(-1, 1)
end

facts("Rational tests") do

    f = 1 // 3
    g = 1 // 3

    @fact @interval(f*g) --> Interval(1.1111111111111109e-01, 1.1111111111111115e-01)
    @fact big(1.)/9 ∈ @interval(f*g) --> true
    @fact @interval(1)/9 ⊆ @interval(f*g) --> true
    @fact @interval(1)/9 ≠ @interval(f*g) --> true

    h = 1/3
    i = 1/3

    @fact @interval(h*i) --> Interval(1.1111111111111109e-01, 1.1111111111111115e-01)
    @fact big(1.)/9 ∈ @interval(1/9) --> true

    @fact @interval(1/9) == @interval(1//9) --> true
end

facts("Floor etc. tests") do
    a = @interval(0.1)
    b = Interval(0.1, 0.1)
    @fact dist(a,b) <= eps(a) --> true

    @fact floor(@interval(0.1, 1.1)) --> Interval(0, 1)
    @fact round(@interval(0.1, 1.1), RoundDown) --> Interval(0, 1)
    @fact ceil(@interval(0.1, 1.1)) --> Interval(1, 2)
    @fact round(@interval(0.1, 1.1), RoundUp) --> Interval(1, 2)
    @fact sign(@interval(0.1, 1.1)) --> Interval(1.0)
    @fact trunc(@interval(0.1, 1.1)) --> Interval(0.0, 1.0)
    @fact round(@interval(0.1, 1.1), RoundToZero) --> Interval(0.0, 1.0)
    @fact round(@interval(0.1, 1.1)) --> Interval(0.0, 1.0)
    @fact round(@interval(0.1, 1.5)) --> Interval(0.0, 2.0)
    @fact round(@interval(-1.5, 0.1)) --> Interval(-2.0, 0.0)
    @fact round(@interval(-2.5, 0.1)) --> Interval(-2.0, 0.0)
    @fact round(@interval(0.1, 1.1), RoundTiesToEven) --> Interval(0.0, 1.0)
    @fact round(@interval(0.1, 1.5), RoundTiesToEven) --> Interval(0.0, 2.0)
    @fact round(@interval(-1.5, 0.1), RoundTiesToEven) --> Interval(-2.0, 0.0)
    @fact round(@interval(-2.5, 0.1), RoundTiesToEven) --> Interval(-2.0, 0.0)
    @fact round(@interval(0.1, 1.1), RoundTiesToAway) --> Interval(0.0, 1.0)
    @fact round(@interval(0.1, 1.5), RoundTiesToAway) --> Interval(0.0, 2.0)
    @fact round(@interval(-1.5, 0.1), RoundTiesToAway) --> Interval(-2.0, 0.0)
    @fact round(@interval(-2.5, 0.1), RoundTiesToAway) --> Interval(-3.0, 0.0)

    # :wide tests
    setrounding(Interval, :wide)
    setprecision(Interval, Float64)

    a = @interval(-3.0, 2.0)
    @fact a --> Interval(-3.0, 2.0)
    @fact a^3 --> Interval(-27.000000000000004, 8.000000000000002)
    @fact Interval(-3,2)^3 --> Interval(-27.000000000000004, 8.000000000000002)

    @fact Interval(-27.0, 8.0)^(1//3) --> Interval(-5.0e-324, 2.0000000000000004)

    setrounding(Interval, :narrow)
end


text:
using Polynomials
import Base: show, *, +

"""
    PolyExp(p::Vector{T}, a::T, b::T)
    PolyExp(pol::Polynomial{T},a::T,b::T)

On the model of `Polynomial` from package `Polynomials`, construct a function that is a polynome multiply by an exponential
function. The exponential is an exponential of an affine function ``a x + b``.
The polynome is construct from its coefficients `p`, lowest order first.

If ``f = (p_n x^n + \\ldots + p_2 x^2 + p_1 x + p_0)e^{a x + b}``, we construct this through
`PolyExp([p_0, p_1, ..., p_n], a, b)`.
It is also possible to construct it directly from the polynome.

In the sequels some methods with the same name than for Polynomial are implemented
(`derivative`, `integrate`, `strings`, ...) but not all, only the methods needed are developped.

# Arguments :
- `p::Vector{T}` or pol::Polynomial{T} : vector of coefficients of the polynome, or directly the polynome.
- `a::T`, `b::T` : coefficients of affine exponentiated function.

# Examples
```julia
julia> pe=PolyExp([1,2,3],2,1)
PolyExp(Polynomial(1 + 2*x + 3*x^2)*exp(2*x + 1))

julia> pe(0)
2.718281828459045

julia> pe(1)
120.51322153912601
```
"""
struct PolyExp{T}
    p::Polynomial{T}
    a::T
    b::T
    PolyExp(p::Vector{T}, a::T, b::T) where {T<:Number} = new{T}(Polynomial{T}(p), a, b)
    PolyExp(pol::Polynomial{T}, a::T, b::T) where {T<:Number} = new{T}(pol, a, b)
end
function _printNumberPar(x::Number)
    return isreal(x) ? "$(real(x))" : (iszero(real(x)) ? "$(imag(x))im" : "($x)")
end
function _printNumber(x::Number)
    return isreal(x) ? "$(real(x))" : (iszero(real(x)) ? "$(imag(x))im" : "$x")
end
function Base.show(io::IO, pe::PolyExp)
    return print(
        io,
        "PolyExp($(pe.p)*exp($(_printNumberPar(pe.a))*x + $(_printNumber(pe.b))))",
    )
end
"""
    derivative(pe::PolyExp)

Construct the derivative of the `pe` function.

# Examples
```julia
julia> derivative(PolyExp([1, 3, -1],3,1))
PolyExp(Polynomial(6 + 7*x - 3*x^2)*exp(3*x + 1))

julia> derivative(PolyExp([1.0+im, 3im, -1, 4.0], 2.0+1.5im,1.0im))
PolyExp(Polynomial((0.5 + 6.5im) - (6.5 - 6.0im)*x + (10.0 - 1.5im)*x^2 + (8.0 + 6.0im)*x^3)*exp((2.0 + 1.5im)*x + 1.0im))
```
"""
function Polynomials.derivative(pe::PolyExp)
    return PolyExp(pe.p * pe.a + Polynomials.derivative(pe.p), pe.a, pe.b)
end

"""
    integrate(pe::PolyExp)

Construct the integrate function of `pe` which is of `PolyExp` type.
The algorithm used is a recursive integration by parts.

# Examples
```julia
julia> integrate(PolyExp([1.0,2,3],2.0,5.0))
PolyExp(Polynomial(0.75 - 0.5*x + 1.5*x^2)*exp(2.0*x + 5.0))

julia> integrate(PolyExp([1.0+0im,2],2.0im,3.0+0im))
PolyExp(Polynomial((0.5 - 0.5im) - 1.0im*x)*exp(2.0im*x + 3.0))
```
"""
function Polynomials.integrate(pe::PolyExp)
    if (pe.a != 0)
        pol = pe.p / pe.a
        if Polynomials.degree(pe.p) > 0
            pol -=
                Polynomials.integrate(PolyExp(Polynomials.derivative(pe.p), pe.a, pe.b)).p /
                pe.a
        end
    else
        pol = Polynomials.integrate(pol)
    end
    return PolyExp(pol, pe.a, pe.b)
end
(pe::PolyExp)(x) = pe.p(x) * exp(pe.a * x + pe.b)
coeffs(pe::PolyExp) = vcat(Polynomials.coeffs(pe.p), pe.a, pe.b)
function Polynomials.integrate(p::PolyExp, v_begin::Number, v_end::Number)
    pol = Polynomials.integrate(p)
    return pol(v_end) - pol(v_begin)
end
*(p1::PolyExp, p2::PolyExp) = PolyExp(p1.p * p2.p, p1.a + p2.a, p1.b + p2.b)
function +(p1::PolyExp, p2::PolyExp)
    @assert (p1.a == p2.a && p1.b == p2.b) "for adding exponents must be equal"
    return PolyExp(p1.p + p2.p, p1.a, p1.b)
end


text:
"""
scipy.optimize module

- [Optimization and Root Finding (scipy.optimize) Reference Guide](https://docs.scipy.org/doc/scipy/reference/optimize.html)

# Examples

You can optimize the Rosen function:

```julia-repl
julia> x0 = [1.3, 0.7, 0.8, 1.9, 1.2];

julia> res = SciPy.optimize.minimize(SciPy.optimize.rosen, x0, method="Nelder-Mead", tol=1e-6)
Dict{Any,Any} with 8 entries:
  "fun"           => 1.94206e-13
  "nit"           => 295
  "nfev"          => 494
  "status"        => 0
  "success"       => true
  "message"       => "Optimization terminated successfully."
  "x"             => [1.0, 1.0, 1.0, 1.0, 1.0]
  "final_simplex" => ([1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0], [1…

```
"""
module optimize

using PyCall
import PyCall: hasproperty # Base.hasproperty in Julia 1.2

import ..pyoptimize
import .._generate_docstring
import ..LazyHelp

const _ignore_funcs = ["optimize"]

# This list can be generated by print_scipy_api_list.py script
all_properties = ["BFGS", "Bounds", "HessianUpdateStrategy", "LbfgsInvHessProduct", "LinearConstraint", "NonlinearConstraint", "OptimizeResult", "OptimizeWarning", "RootResults", "SR1"
, "anderson", "approx_fprime", "basinhopping", "bisect", "bracket", "brent", "brenth", "brentq", "broyden1", "broyden2", "brute", "check_grad", "cobyla", "curve_fit",
"diagbroyden", "differential_evolution", "dual_annealing", "excitingmixing", "fixed_point", "fmin", "fmin_bfgs", "fmin_cg", "fmin_cobyla", "fmin_l_bfgs_b", "fmin_ncg",
 "fmin_powell", "fmin_slsqp", "fmin_tnc", "fminbound", "fsolve", "golden", "lbfgsb", "least_squares", "leastsq", "line_search", "linear_sum_assignment", "linearmixing"
, "linesearch", "linprog", "linprog_verbose_callback", "lsq_linear", "minimize", "minimize_scalar", "minpack", "minpack2", "moduleTNC", "newton", "newton_krylov", "nnl
s", "nonlin", "optimize", "quadratic_assignment", "ridder", "root", "root_scalar", "rosen", "rosen_der", "rosen_hess", "rosen_hess_prod", "shgo", "show_options", "slsq
p", "tnc", "toms748", "zeros"]

for f in all_properties
    f in _ignore_funcs && continue

    sf = Symbol(f)
    @eval @doc LazyHelp(pyoptimize, $f) $sf(args...; kws...) = pycall(pyoptimize.$f, PyAny, args...; kws...)
end

function __init__()
    copy!(pyoptimize, pyimport_conda("scipy.stats", "scipy"))
end

end # module