text:
using Flux, Zygote, LinearAlgebra, Statistics
println("Heat han")
using Test, StochasticDiffEq

using NeuralPDE

using Random
Random.seed!(100)

opt = Flux.ADAM(0.005)  #optimizer
# high-dimensional heat equation
d = 30 # number of dimensions
x0 = fill(8.0f0,d)
tspan = (0.0f0,2.0f0)
dt = 0.5
time_steps = div(tspan[2]-tspan[1],dt)
m = 30 # number of trajectories (batch size)

g(X) = sum(X.^2)
f(X,u,σᵀ∇u,p,t) = 0.0
μ_f(X,p,t) = 0.0
σ_f(X,p,t) = 1.0
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = 10 + d #hidden layer size
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = [Flux.Chain(Dense(d,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d)) for i in 1:time_steps]

alg = NNPDEHan(u0, σᵀ∇u, opt = opt)

ans = solve(prob, alg, verbose = true, abstol=1e-8, maxiters = 1000, dt=dt, trajectories=m)

u_analytical(x,t) = sum(x.^2) .+ d*t
analytical_ans = u_analytical(x0, tspan[end])
error_l2 = sqrt((ans - analytical_ans)^2/ans^2)

println("high-dimensional heat equation")
# println("numerical = ", ans)
# println("analytical = " ,analytical_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 1.0


text:
using TreatmentPanels
using Test
using DataFrames

@testset "Input tests" begin
    test_df = DataFrame(id = ["a", "a", "b", "b"], period = [1, 2, 1, 2], value = 1.0:4.0)

    # Basic functionality
    @test BalancedPanel(test_df, ["a" => 2];
        id_var = :id, t_var = :period, outcome_var = :value) isa BalancedPanel{SingleUnitTreatment, ContinuousTreatment}
    @test BalancedPanel(test_df, "a" => 2;
        id_var = :id, t_var = :period, outcome_var = :value) isa BalancedPanel{SingleUnitTreatment, ContinuousTreatment}

    # Column not present in data
    @test_throws "" BalancedPanel(test_df, "a" => 2; id_var = :wrong_var, t_var = :period, outcome_var = :value)
    @test_throws "" BalancedPanel(test_df, "a" => 2; id_var = :id, t_var = :wrong_var, outcome_var = :value)
    @test_throws "" BalancedPanel(test_df, "a" => 2; id_var = :id, t_var = :period, outcome_var = :wrong_var)

    # Treatment assignment has nonexistent units/time periods
    @test_throws "" BalancedPanel(test_df, "c" => 2; id_var = :id, t_var = :period, outcome_var = :value)
end

ext: jl


text:
# TODO: add keyword `shrink` to give a consistant result on Base
#       when this is done, then we can propose this change to upstream Base
"""
    fdiff(A::AbstractArray; dims::Int, rev=false, boundary=:periodic)

A one-dimension finite difference operator on array `A`. Unlike `Base.diff`, this function doesn't
shrink the array size.

Take vector as an example, it computes `(A[2]-A[1], A[3]-A[2], ..., A[1]-A[end])`.

# Keywords

- `rev::Bool`
  If `rev==true`, then it computes the backward difference
  `(A[end]-A[1], A[1]-A[2], ..., A[end-1]-A[end])`.
- `boundary`
  By default it computes periodically in the boundary, i.e., `:periodic`.
  In some cases, one can fill zero values with `boundary=:zero`.

# Examples

```jldoctest; setup=:(using ImageBase: fdiff)
julia> A = [2 4 8; 3 9 27; 4 16 64]
3×3 $(Matrix{Int}):
 2   4   8
 3   9  27
 4  16  64

julia> diff(A, dims=2) # this function exists in Base
3×2 $(Matrix{Int}):
  2   4
  6  18
 12  48

julia> fdiff(A, dims=2)
3×3 $(Matrix{Int}):
  2   4   -6
  6  18  -24
 12  48  -60

julia> fdiff(A, dims=2, rev=true) # reverse diff
3×3 $(Matrix{Int}):
  -6   2   4
 -24   6  18
 -60  12  48

julia> fdiff(A, dims=2, boundary=:zero) # fill boundary with zeros
3×3 $(Matrix{Int}):
  2   4  0
  6  18  0
 12  48  0
```

See also [`fdiff!`](@ref) for the in-place version.
"""
fdiff(A::AbstractArray; kwargs...) = fdiff!(similar(A), A; kwargs...)

"""
    fdiff!(dst::AbstractArray, src::AbstractArray; dims::Int, rev=false, boundary=:periodic)

The in-place version of [`ImageBase.fdiff`](@ref)
"""
function fdiff!(dst::AbstractArray, src::AbstractArray;
        dims=_fdiff_default_dims(src),
        rev=false,
        boundary::Symbol=:periodic)
    isnothing(dims) && throw(UndefKeywordError(:dims))
    axes(dst) == axes(src) || throw(ArgumentError("axes of all input arrays should be equal. Instead they are $(axes(dst)) and $(axes(src))."))
    N = ndims(src)
    1 <= dims <= N || throw(ArgumentError("dimension $dims out of range (1:$N)"))

    r = axes(src)
    r0 = ntuple(i -> i == dims ? UnitRange(first(r[i]), last(r[i]) - 1) : UnitRange(r[i]), N)
    r1 = ntuple(i -> i == dims ? UnitRange(first(r[i])+1, last(r[i])) : UnitRange(r[i]), N)

    d0 = ntuple(i -> i == dims ? UnitRange(last(r[i]), last(r[i])) : UnitRange(r[i]), N)
    d1 = ntuple(i -> i == dims ? UnitRange(first(r[i]), first(r[i])) : UnitRange(r[i]), N)

    if rev
        dst[r1...] .= view(src, r1...) .- view(src, r0...)
        if boundary == :periodic
            dst[d1...] .= view(src, d1...) .- view(src, d0...)
        elseif boundary == :zero
            dst[d1...] .= zero(eltype(dst))
        else
            throw(ArgumentError("Wrong boundary condition $boundary"))
        end
    else
        dst[r0...] .= view(src, r1...) .- view(src, r0...)
        if boundary == :periodic
            dst[d0...] .= view(src, d1...) .- view(src, d0...)
        elseif boundary == :zero
            dst[d0...] .= zero(eltype(dst))
        else
            throw(ArgumentError("Wrong boundary condition $boundary"))
        end
    end

    return dst
end

_fdiff_default_dims(A) = nothing
_fdiff_default_dims(A::AbstractVector) = 1


text:
struct MonteCarlo <: AbstractMonteCarlo
    n::Integer
    MonteCarlo(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

struct SobolSampling <: AbstractQuasiMonteCarlo
    n::Integer
    SobolSampling(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

struct HaltonSampling <: AbstractQuasiMonteCarlo
    n::Integer
    HaltonSampling(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

struct LatinHypercubeSampling <: AbstractQuasiMonteCarlo
    n::Integer
    LatinHypercubeSampling(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

function sample(inputs::Array{<:UQInput}, sim::MonteCarlo)
    return sample(inputs, sim.n)
end

function sample(inputs::Array{<:UQInput}, sim::AbstractQuasiMonteCarlo)
    random_inputs = filter(i -> isa(i, RandomUQInput), inputs)
    deterministic_inputs = filter(i -> isa(i, DeterministicUQInput), inputs)

    n_rv = count_rvs(random_inputs)

    u = qmc_samples(sim, n_rv)

    samples = quantile.(Normal(), u)
    samples = DataFrame(names(random_inputs) .=> eachcol(samples))

    if !isempty(deterministic_inputs)
        samples = hcat(samples, sample(deterministic_inputs, sim.n))
    end

    to_physical_space!(inputs, samples)

    return samples
end

function qmc_samples(sim::SobolSampling, rvs::Integer)
    return transpose(QuasiMonteCarlo.sample(sim.n, zeros(rvs), ones(rvs), SobolSample()))
end

function qmc_samples(sim::HaltonSampling, rvs::Integer)
    h = HaltonPoint(rvs; length=sim.n)

    u = zeros(sim.n, rvs)

    for (i, hp) in enumerate(h)
        u[i, :] = hp
    end

    return u
end

function qmc_samples(sim::LatinHypercubeSampling, rvs::Integer)
    return transpose(
        QuasiMonteCarlo.sample(sim.n, zeros(rvs), ones(rvs), LatinHypercubeSample())
    )
end


text:
struct SQLParserNoParseTreeException <: Exception end

function Base.showerror(io::IO, ::SQLParserNoParseTreeException)
    Base.print(io, "SQLParserNoParseTreeException")
end

struct SQLParserException <: Exception
    inner::Ptr{PgQueryError}
end

function convert_cstring_to_str(c_str)
    c_str == C_NULL ? "" : unsafe_string(c_str)
end

function Base.showerror(io::IO, pe::SQLParserException)
    err = unsafe_load(pe.inner)
    Base.print(io, "SQLParserException: " * convert_cstring_to_str(err.message) *
        " at location " * string(err.cursorpos) * " (" * convert_cstring_to_str(err.funcname) *
        " in " * convert_cstring_to_str(err.filename) * ":" * string(err.lineno) * ")\n")
end


text:
using SurrealNumbers
out_dir = "Data/"

# addition table
x = convert.(SurrealFinite, [0, 1/2, 1, 2])
A = zeros(SurrealFinite, length(x), length(x))
B = zeros(Int, length(x), length(x))
for i=1:length(x)
    for j=1:length(x)
        A[i,j] = x[i] + x[j]
        B[i,j] = generation(x[i]) + generation(x[j])
    end
end
float.(A)
A
generation.(A)
A .== canonicalise.(A)
iscanonical.(A)
B

a10 = convert(SurrealFinite, 3/2)
a11 = convert(SurrealFinite, 1/2) + convert(SurrealFinite, 1)
file = "$(out_dir)addition_ex_10.dot"
FID = open(file, "w")
surreal2dag(FID, a10)
close(FID)
run(`dot -Tpdf -O $file`)
file = "$(out_dir)addition_ex_11.dot"
FID = open(file, "w")
surreal2dag(FID, a11)
close(FID)
run(`dot -Tpdf -O $file`)

# subtraction table
x = convert.(SurrealFinite, [0, 1/2, 1, 2])
A = zeros(SurrealFinite, length(x), length(x))
B = zeros(Int, length(x), length(x))
for i=1:length(x)
    for j=1:length(x)
        A[i,j] = x[i] - x[j]
        B[i,j] = generation(x[i]) + generation(x[j])
    end
end
float.(A)
A
generation.(A)
A .== canonicalise.(A)
iscanonical.(A)
B

n = 2
x = convert.(SurrealFinite, (0:2^n)/2^n)
A = zeros(SurrealFinite, 2^n+1, 2^n+1)
for i=1:2^n +1
    for j=1:2^n +1
        A[i,j] = x[i] + x[j]
    end
end
float.(A)
A
generation.(A)
A .== canonicalise.(A)
iscanonical.(A)

n = 1
x = convert.(SurrealFinite, 0:n)
A = zeros(SurrealFinite, n+1, n+1)
for i=1:n+1
    for j=1:n+1
        A[i,j] = x[i] * x[j]
    end
end
float.(A)
generation.(A)
A .== canonicalise.(A)
iscanonical.(A)

# analyse.(A)

# function analyse(s::SurrealFinite)
#     println("depth = ", generation(s) )

# end


text:
module NetcdfIO

using DataFrames: DataFrame
using DocStringExtensions: METHODLIST
using NCDatasets: Dataset, defDim, defVar

# export public functions
export add_nc_dim!, append_nc!, create_nc!, dimname_nc, grow_nc!, read_nc, save_nc!, size_nc, varname_nc

# constants
const ATTR_LAT   = Dict("description" => "Latitude", "unit" => "°");
const ATTR_LON   = Dict("description" => "Longitude", "unit" => "°");
const ATTR_CYC   = Dict("description" => "Cycle index", "unit" => "-");
const ATTR_ABOUT = Dict("about" => "This is a file generated using NetcdfIO.jl",
                        "notes" => "NetcdfIO.jl uses NCDatasets.jl to read and write NC files");

# include the files
include("append.jl")
include("create.jl")
include("grow.jl"  )
include("info.jl"  )
include("read.jl"  )
include("save.jl"  )

end # module


text:
# ===============================
# Test minimal (default) options
# ===============================

# linear ODE: x' = Ax
A = [ 0.0509836  0.168159  0.95246   0.33644
      0.42377    0.67972   0.129232  0.126662
      0.518654   0.981313  0.489854  0.588326
      0.38318    0.616014  0.518412  0.778765]

# initial set
X0 = BallInf(ones(4), 0.1)

# default options (computes all variables)
s = solve(IVP(LCS(A), X0), :T=>0.1)
@test dim(s.Xk[1].X) == 4
@test s.options isa Options

# two variables and custom partition
s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4]))
@test dim(s.Xk[1].X) == 4

# the default partition is used.
s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1),
          op=BFFPSV18(:vars=>[1,3]))
@test dim(s.Xk[1].X) == 2

# ===============================
# Test projection
# ===============================
s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1, :project_reachset=>false),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4]))
ps = project(s)

# ===============
# Test check mode
# ===============

A = [ 0.0509836  0.168159  0.95246   0.33644
      0.42377    0.67972   0.129232  0.126662
      0.518654   0.981313  0.489854  0.588326
      0.38318    0.616014  0.518412  0.778765]

# check that x1 + x2 <= 2 doesn't hold
# (this computes only block 1, needed for the property)
s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1, :mode=>"check",
          :property=>SafeStatesProperty(LinearConstraint([1., 1., 0., 0.], 2.))),
          op=BFFPSV18(:vars=>[1,2], :partition=>[1:2, 3:4]))
@test s.violation == 1

# same but computing the two blocks, even if the second block is not needed
# for the property
s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1, :mode=>"check",
          :property=>SafeStatesProperty(LinearConstraint([1., 1., 0., 0.], 2.))),
          op=BFFPSV18(:vars=>[1,2,3], :partition=>[1:2, 3:4]))
@test s.violation == 1

# check that x1 - x2 <= 2 holds
s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1, :mode=>"check",
          :property=>SafeStatesProperty(LinearConstraint([1., -1., 0., 0.], 2.))),
          op=BFFPSV18(:vars=>[1,2,3], :partition=>[1:2, 3:4]))
@test s.violation == -1

# ===============================
# Test reachability options
# ===============================

# template directions
s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1, :ε_proj=>1e-5, :set_type_proj=>HPolygon),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:4],
                      :block_options => LazySets.Approximations.OctDirections))

s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :sih_method=>"lazy"))

s = solve(IVP(LCS(A), X0),
          Options(:T=>0.1),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :sih_method=>"concrete"))

s = solve(IVP(LCS(sparse(A)), X0),
          Options(:T=>0.1),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :exp_method=>"lazy",
                      :assume_sparse=>false))

s = solve(IVP(LCS(sparse(A)), X0),
          Options(:T=>0.1),
          op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4], :exp_method=>"lazy",
                      :assume_sparse=>true))

# uses Interval for set_type_init and set_type_iter but Hyperrectangle for
# set_type_proj
s = solve(IVP(LCS(sparse(A)), X0),
          Options(:T=>0.1),
          op=BFFPSV18(:vars=>[1,3], :partition=>[[i] for i in 1:4],
                      :exp_method=>"lazy",
                      :block_options=>Interval))

# ===============================
# System with an odd dimension
# ===============================
A = randn(5, 5); X0 = BallInf(ones(5), 0.1)
s = solve(IVP(LCS(sparse(A)), X0), Options(:T=>0.1),
    op=BFFPSV18(:vars=>[1,3], :partition=>[1:2, 3:4, [5]], :block_options=>
                [HPolygon, HPolygon, Interval]))

A = [1 2 1.; 0 0. 1; -2 1 4]
X0 = BallInf(ones(3), 0.1)
system = IVP(LCS(sparse(A)), X0)
s = solve(system, Options(:T=>0.1),
          op=BFFPSV18(:vars=>1:3, :partition=>[1:2, 3:3]))

# =======================================================
# Affine ODE with a nondeterministic input, x' = Ax + Bu
# =======================================================
# linear ODE: x' = Ax + Bu, u ∈ U
A = [0.0509836  0.168159  0.95246   0.33644
     0.42377    0.67972   0.129232  0.126662
     0.518654   0.981313  0.489854  0.588326
     0.38318    0.616014  0.518412  0.778765]

B = [0.866688  0.771231
     0.065935  0.349839
     0.109643  0.904222
     0.292474  0.227857]

# non-deterministic inputs
U = Interval(0.99, 1.01) × Interval(0.99, 1.01)

# initial set
X0 = BallInf(ones(4), 0.1)

# inputs
U1 = ConstantInput(U)
U2 = U  # use internal wrapping

for inputs in [U1, U2]
    # default options (computes all variables)
    s = solve(IVP(CLCCS(A, B, nothing, U), X0), :T=>0.1)
end

# ===============================
# Test GLGM06
# ===============================

# linear ODE: x' = Ax
A = [ 0.0509836  0.168159  0.95246   0.33644
      0.42377    0.67972   0.129232  0.126662
      0.518654   0.981313  0.489854  0.588326
      0.38318    0.616014  0.518412  0.778765]

# initial set
X0 = BallInf(ones(4), 0.1)

# default options
s = solve(IVP(LCS(A), X0), Options(:T=>0.1), op=GLGM06())

# specify maximum order of zonotopes
s = solve(IVP(LCS(A), X0), Options(:T=>0.1), op=GLGM06(:max_order=>15))

# affine system, x' = Ax + Bu
s = solve(IVP(CLCCS(A, B, nothing, U), X0), Options(:T=>0.1), op=GLGM06())

# ======================
# Unbounded-time setting
# ======================
X = LinearConstraint([1., 1., 0., 0.], 9.5)
property = SafeStatesProperty(LinearConstraint([1., 1., 0., 0.], 10.))

# default options (computes all variables)
s = Reachability.solve(IVP(CLCCS(A, B, X, U), X0), Options(:T=>Inf))
s = solve(IVP(LCS(A), X0), :T=>Inf, :mode=>"check", :property=>property)
s = solve(IVP(CLCCS(A, B, X, U), X0),
          :T=>Inf, :mode=>"check", :property=>property)

# =================================================================
# Affine ODE with a nondeterministic input, x' = Ax + b, no inputs
# =================================================================
A = [ 0.0509836  0.168159  0.95246   0.33644
      0.42377    0.67972   0.129232  0.126662
      0.518654   0.981313  0.489854  0.588326
      0.38318    0.616014  0.518412  0.778765]
b =  [1.0, -1, 0, 1.]
X = HalfSpace([-1.0, 0.0, 0.0, 0.0], 0.0)
X0 = BallInf(zeros(4), 0.5)
p = IVP(ConstrainedAffineContinuousSystem(A, b, X), X0)
sol = solve(p, :T=>0.1)

# ==================================
# Test TMJets with Van der Pol model
# ==================================
using TaylorIntegration
using TaylorModels: @taylorize
using Reachability: solve

@taylorize function vanderPol!(t, x, dx)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end

𝑆 = BlackBoxContinuousSystem(vanderPol!, 2)
X0 = Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45])
𝑃 = InitialValueProblem(𝑆, X0)

# reach mode
𝑂 = Options(:T=>7.0, :mode=>"reach")
solve(𝑃, 𝑂, op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2));

# check mode
property=(t,x)->x[2] <= 2.75
𝑂 = Options(:T=>7.0, :mode=>"check", :property=>property)
solve(𝑃, 𝑂, op=TMJets(:abs_tol=>1e-10, :orderT=>10, :orderQ=>2))


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Variables #
# ========================== #
include("variable.jl")
# Objective #
# ========================== #
include("objectives/objective.jl")
# Constraint #
# ========================== #
include("constraints/constraints.jl")
# SingleDecision #
# ========================== #
# SingleDecision support is achieved by adding a few methods to MathOptInterface
# as well as two specialized functionize bridges which behave as their SingleVariable
# counterparts
include("functionize.jl")

function MOIB.Variable.function_for(map::MOIB.Variable.Map, ci::MOI.ConstraintIndex{SingleDecision})
    return SingleDecision(MOI.VariableIndex(ci.value))
end

function MOIB.Variable.function_for(map::MOIB.Variable.Map, ci::MOI.ConstraintIndex{VectorOfDecisions})
    decisions = MOI.VariableIndex[]
    for i in ci.value:-1:-length(map.bridges)
        vi = MOI.VariableIndex(i)
        if map.index_in_vector[-vi.value] == -1
            continue
        elseif bridge_index(map, vi) == -ci.value
            push!(decisions, vi)
        else
            break
        end
    end
    return VectorOfDecisions(decisions)
end

function MOIB.is_variable_bridged(
    b::MOIB.AbstractBridgeOptimizer,
    ci::MOI.ConstraintIndex{<:Union{SingleDecision, VectorOfDecisions}})
    return ci.value < 0 && !haskey(MOIB.Constraint.bridges(b), ci)
end

function MOIB.bridged_function(b::MOIB.AbstractBridgeOptimizer, f::SingleDecision)
    return f
end

function MOIB.unbridged_function(bridge::MOIB.AbstractBridgeOptimizer,
                                 f::Union{SingleDecision, VectorOfDecisions})
    return f
end

function MOI.set(b::MOIB.AbstractBridgeOptimizer,
                 attr::MOI.ObjectiveFunction,
                 f::SingleDecision)
    if MOIB.Variable.has_bridges(MOIB.Variable.bridges(b))
        if MOIB.is_bridged(b, f.decision)
            BridgeType = MOIB.Objective.concrete_bridge_type(
                FunctionizeDecisionObjectiveBridge{Float64}, typeof(f))
            MOIB._bridge_objective(b, BridgeType, f)
        end
    end
end

function MOI.add_constraint(b::MOIB.AbstractBridgeOptimizer,
                            f::SingleDecision,
                            s::MOI.AbstractSet)
    if MOIB.is_bridged(b, f.decision)
        if MOI.is_valid(b, MOI.ConstraintIndex{MOI.SingleVariable, typeof(s)}(f.decision.value))
            # The other constraint could have been through a variable bridge.
            error("Cannot add two `SingleDecision`-in-`$(typeof(s))`",
                  " on the same decision $(f.decision).")
        end
        BridgeType = MOIB.Constraint.concrete_bridge_type(
            SingleDecisionConstraintBridge{Float64}, typeof(f), typeof(s))
        return MOIB.add_bridged_constraint(b, BridgeType, f, s)
    end
    error("`SingleDecision`-in-`$(typeof(s))` is only supported through variable bridging.")
end

function MOI.add_constraint(b::MOIB.AbstractBridgeOptimizer,
                            f::VectorOfDecisions,
                            s::MOI.AbstractSet)
    if any(vi -> MOIB.is_bridged(b, vi), f.decisions)
        if MOI.is_valid(b, MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(s)}(first(f.decisions).value))
            # The other constraint could have been through a variable bridge.
            error("Cannot add two `VectorOfDecisions`-in-`$(typeof(s))`",
                  " on the same first decision $(first(f.decisions)).")
        end
        if !MOIB.is_bridged(b, first(f.decisions)) && !MOIB.is_bridged(b, typeof(f), typeof(s))
            # The index of the contraint will have positive value hence
            # it would clash with the index space of `b.model` since
            # the constraint type is normally not bridged.
            error("Cannot `VectorOfDecisions`-in-`$(typeof(s))` for",
                  " which some decisions are bridged but not the",
                  " first one `$(first(f.decisions))`.")
        end
        BridgeType = MOIB.Constraint.concrete_bridge_type(
            VectorDecisionConstraintBridge{Float64}, typeof(f), typeof(s))
        return MOIB.add_bridged_constraint(b, BridgeType, f, s)
    end
    error("`VectorOfDecisions`-in-`$(typeof(s))` is only supported through variable bridging.")
end

function MOI.add_constraints(b::MOIB.AbstractBridgeOptimizer,
                             f::Vector{SingleDecision},
                             s::Vector{S}) where S <: MOI.AbstractSet
    MOI.add_constraint.(b, f, s)
end

function MOI.add_constraints(b::MOIB.AbstractBridgeOptimizer,
                             f::Vector{VectorOfDecisions},
                             s::Vector{S}) where S <: MOI.AbstractSet
    MOI.add_constraint.(b, f, s)
end

function MOI.get(b::MOIB.AbstractBridgeOptimizer, attr::MOI.ConstraintSet,
                 ci::MOI.ConstraintIndex{SingleDecision})
    return if MOIB.is_bridged(b, ci)
        MOI.throw_if_not_valid(b, ci)
        MOIB.call_in_context(b, ci, bridge -> MOI.get(b, attr, bridge))
    else
        MOI.get(b.model, attr, ci)
    end
end

function MOIB.bridged_constraint_function(
    b::MOIB.AbstractBridgeOptimizer, f::SingleDecision,
    set::MOI.AbstractScalarSet)
    return MOIB.bridged_function(b, f), set
end

function MOIB.unbridged_constraint_function(
    b::MOIB.AbstractBridgeOptimizer, f::SingleDecision)
    return MOIB.unbridged_function(b, f)
end

# Bridge addition #
# ========================== #
function add_decision_bridges!(model::JuMP.Model)
    add_bridge(model, DecisionBridge)
    add_bridge(model, DecisionsBridge)
    add_bridge(model, AffineDecisionObjectiveBridge)
    add_bridge(model, QuadraticDecisionObjectiveBridge)
    add_bridge(model, AffineDecisionConstraintBridge)
    add_bridge(model, QuadraticDecisionConstraintBridge)
    add_bridge(model, VectorAffineDecisionConstraintBridge)
    add_bridge(model, FunctionizeDecisionObjectiveBridge)
end


text:
# using PyPlot
using ConvDiffMIPDECO
using jInv.Mesh
# using jInvVis
using jInv.LinearSolvers
using MUMPSjInv

domain      = [0 2 0 1]
nc		    = [128 64]*4
pNoise      = .05
noBoreHoles = 30
Mc          = getRegularMesh(domain,nc)
P,mf        = getRandomBoreholes2D(noBoreHoles,nc,nc) 		# Randomly selecting borehols
pForFEM     = getConvDiffFEMParam(Mc,sig=0.01,Ainv=getMUMPSsolver());
v = kron([1;0.0],ones(prod(Mc.n+2)))
pFor  = getConvDiffParam(Mc,v,sig=pForFEM.sig,Ainv=getMUMPSsolver(),bc=(:dir,:neu,:neu,:neu));
w = zeros(Mc.n[1],Mc.n[2]);
w[30:40,30:40] = 1;
viewImage2D(w,Mc)
title("source term")
colorbar()
pFor.sig

dobsFEM, = getData(vec(w),pForFEM);
dobs,    = getData(vec(w),pFor);

subplot(3,1,1)
viewImage2D(pForFEM.Fields,getPaddedMesh(Mc))
ylabel("FEM")
colorbar()
title("Fields")

subplot(3,1,2)
viewImage2D(pFor.Fields,Mc)
ylabel("FD")
colorbar()

subplot(3,1,3)
viewImage2D(getNodalAverageMatrix(Mc)*pForFEM.Fields-pFor.Fields,Mc)
colorbar()
ylabel("Difference")

println("relative error (PDE solve FEM vs FD): $(norm(vec(getNodalAverageMatrix(Mc)*pForFEM.Fields-pFor.Fields))/norm(vec(pFor.Fields)))")

"""
fdAnisoTV(W,Mc)

prototype implementation of regularizer using for loops
"""
function fdAnisoTV(W,Mc)
    W = reshape(W,Mc.n[1], Mc.n[2])
    L = Mc.h

    Reg = 0.0;
    # compute first term of regularizer (derivatives in first dimension)
    for i=2:Mc.n[1]
        for j=1:Mc.n[2]
            Reg += Mc.h[2]*abs(W[i,j] - W[i-1,j])
        end
    end
    # compute first term of regularizer (derivatives in second dimension)
    for i=1:Mc.n[1]
        for j=2:Mc.n[2]
            Reg += Mc.h[1]*abs(W[i,j] - W[i,j-1])
        end
    end
    return Reg
end

Rfd = fdAnisoTV(w,Mc)

using jInv.InverseSolve
Rjinv = anisoTVReg(vec(w),0*vec(w),Mc,eps=0)[1]
println("Relative error (regularizer for loop): $(abs(Rfd-Rjinv)/abs(Rjinv))")

"""
fdPDE(U,Mc)

prototype implementation of PDE using for loops

here, assume that v = [1;0]
"""
function fdPDE(U,Mc;c=0.01)
    Up = zeros(Mc.n[1]+2, Mc.n[2]+2)
    Up[2:end-1,2:end-1] = reshape(U,Mc.n[1], Mc.n[2])

    # left boundary: 0-Dirichlet
    Up[1,2:end-1] = -U[1,:]
    # right boundary: 0-Neumann
    Up[end,2:end-1] = U[end,:]
    # top boundary: 0-Neumann
    Up[2:end-1,end] = U[:,end]
    # bottom boundary: 0-Neumann
    Up[2:end-1,1] = U[:,1]

    L = Mc.h

    # allocate output (with ghostpoints for simpler indexing)
    V = zeros(Mc.n[1]+2, Mc.n[2]+2)
    for i=2:Mc.n[1]+1
        for j=2:Mc.n[2]+1
        V[i,j] = (c/(L[1]^2))*(2*Up[i,j]-Up[i-1,j]-Up[i+1,j]) +
                 (c/(L[2]^2))*(2*Up[i,j]-Up[i,j-1]-Up[i,j+1]) +
                 (Up[i,j]-Up[i-1,j])/L[1]
        end
    end
    V = V[2:Mc.n[1]+1,2:Mc.n[2]+1]
    return V
end
U = randn(Mc.n[1],Mc.n[2])
# U = copy(w)
Vfd = fdPDE(copy(U),Mc,c=pFor.sig);

Vjinv = reshape(pFor.A*vec(U),Mc.n[1],Mc.n[2]);

subplot(3,1,1);
viewImage2D(Vfd,Mc)
colorbar()

subplot(3,1,2);
viewImage2D(Vjinv,Mc)
colorbar()

subplot(3,1,3);
viewImage2D(abs.(Vfd-Vjinv),Mc)
colorbar()

println("relative error (conv diff for loop): $(norm(vec(Vjinv-Vfd))/norm(vec(Vjinv)))")

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
#-----------------------------------------------------------------------# LinReg
"""
    LinReg()

Linear regression, optionally with element-wise ridge regularization.

# Example

    x = randn(100, 5)
    y = x * (1:5) + randn(100)
    o = fit!(LinReg(), zip(eachrow(x),y))
    coef(o)
    coef(o, .1)
    coef(o, [0,0,0,0,Inf])
"""
mutable struct LinReg{W} <: OnlineStat{XY}
    β::Vector{Float64}
    A::Matrix{Float64}
    weight::W
    n::Int
end
LinReg(p=0;weight=EqualWeight()) = LinReg(zeros(p), zeros(1, 1), weight, 0)

function _fit!(o::LinReg, xy)
    γ = o.weight(o.n += 1)
    x, y = xy
    if o.n == 1
        p = length(x)
        o.β = zeros(p)
        o.A = zeros(p + 1, p + 1)
    end
    for j in 1:(size(o.A, 2)-1)
        xj = x[j]
        o.A[j, end] = smooth(o.A[j, end], xj * y, γ)    # x'y
        for i in 1:j
            o.A[i, j] = smooth(o.A[i,j], x[i] * xj, γ)  # x'x
        end
    end
    o.A[end] = smooth(o.A[end], y * y, γ)  # y'y
end
value(o::LinReg, args...) = coef(o, args...)
function coef(o::LinReg)
    o.β[:] = Symmetric(o.A[1:(end-1), 1:(end-1)]) \ o.A[1:(end-1), end]
end
function coef(o::LinReg, λ::Real)
    o.β[:] = Symmetric(o.A[1:(end-1), 1:(end-1)] + λ*I) \ o.A[1:(end-1), end]
end
function coef(o::LinReg, λ::AbstractVector{<:Real})
    o.β[:] = Symmetric(o.A[1:(end-1), 1:(end-1)] + Diagonal(λ)) \ o.A[1:(end-1), end]
end

predict(o::LinReg, x::AbstractVector) = x'o.β
predict(o::LinReg, x::AbstractMatrix) = x * o.β

function _merge!(o::LinReg, o2::LinReg)
    o.n += o2.n
    smooth!(o.A, o2.A, nobs(o2) / nobs(o))
end

#-----------------------------------------------------------------------# LinRegBuilder
"""
    LinRegBuilder(p)

Create an object from which any variable can be regressed on any other set of variables,
optionally with element-wise ridge regularization.  The main function to use with
`LinRegBuilder` is `coef`:

    coef(o::LinRegBuilder, λ = 0; y=1, x=[2,3,...], bias=true, verbose=false)

Return the coefficients of a regressing column `y` on columns `x` with ridge (`abs2` penalty)
parameter `λ`.  An intercept (`bias`) term is added by default.

# Examples

    x = randn(1000, 10)
    o = fit!(LinRegBuilder(), eachrow(x))

    coef(o; y=3, verbose=true)

    coef(o; y=7, x=[2,5,4])
"""
mutable struct LinRegBuilder{W} <: OnlineStat{VectorOb}
    A::Matrix{Float64}  #  x'x, pretend that x = [x, 1]
    weight::W
    n::Int
end
function LinRegBuilder(p=0; weight = EqualWeight())
    LinRegBuilder(Matrix{Float64}(undef, p + 1, p + 1), weight, 0)
end
function Base.show(io::IO, o::LinRegBuilder)
    print(io, "LinRegBuilder of $(nvars(o)) variables")
end
nvars(o::LinRegBuilder) = size(o.A, 1) - 1

function _fit!(o::LinRegBuilder, x)
    o.n += 1
    if o.n == 1
        o.A = zeros(length(x) + 1, length(x) + 1)
        o.A[end] = 1.0
    end
    smooth_syr!(o.A, BiasVec(x), o.weight(o.n))
end
value(o::LinRegBuilder, args...) = coef(o, args...)

xs(o::LinRegBuilder, y) = setdiff(1:size(o.A, 1) - 1, y)
function add_diag!(S, λ::Number)
    for i in 1:(size(S, 1) - 1)
        S[i, i] += λ
    end
end
function add_diag!(S, λ::AbstractVector)
    for i in 1:(size(S, 1) - 1)
        S[i, i] += λ[i]
    end
end

function coef(o::LinRegBuilder, λ=0.0; y=1, x=xs(o,y), bias=true, verbose=false)
    verbose && @info("Regress $y on $x $(bias ? "with bias" : "")")
    inds = collect(x)
    bias && push!(inds, size(o.A, 1))
    push!(inds, y)
    S = Symmetric(o.A)[inds, inds]
    add_diag!(S, λ)
    p = length(inds) - 1
    return S[1:p, 1:p] \ S[1:p, end]
end

function _merge!(o::LinRegBuilder, o2::LinRegBuilder)
    o.n += o2.n
    smooth!(o.A, o2.A, nobs(o2) / nobs(o))
end


text:
function compute_s(n)
    sn = 0
    for i = 1:n
        sn += 1/i
    end
    sn
end

compute_s(100)

for x = 1:2:20      # Steps of 2 - only odd values x
    print(x, " ")
end

for x = 1:0.2:5     # Non-integers are ok
    print(x, " ")
end

for x = 10:-1:-5    # Use a negative step to count down
    print(x, " ")
end

for x = 10:2:5      # No repetitions since start > stop and step is positive
    print(x, " ")
end

function my_factorial(n)
    y = 1
    for i = 2:n
        y *= i
    end
    y
end

my_factorial(5)

my_factorial(20)     # This is OK

my_factorial(30)     # Too large - overflow

my_factorial(Int128(30))

factorial(30)

function taylor_cos_bad(x,n)
    y = 0
    for k = 0:n
        y += (-1)^k / factorial(2k) * x^2k
    end
    y
end

# x = 0.25, excellent approximation of degree 4
println(taylor_cos_bad(0.25, 2))  # Taylor approximation
println(cos(0.25))                # true value

# x = 10, very poor approximation of degree 4
println(taylor_cos_bad(10, 2))  # Taylor approximation
println(cos(10))                # true value

# x = 10, try approximation of degree 30
println(taylor_cos_bad(10, 15)) # Taylor approximation
println(cos(10))                # true value

function taylor_cos(x,n)
    term = 1
    y = 1
    for k = 1:n
        term *= -x^2 / ((2k-1) * 2k)
        y += term
    end
    y
end

# x = 10, try approximation of degree 30
println(taylor_cos(10, 15)) # Taylor approximation
println(cos(10))            # true value

# x = 10, try approximation of degree 100
println(taylor_cos(10, 50)) # Taylor approximation
println(cos(10))            # true value

x = 10
y = 10

for i = 1:5
    z = i        # Local scope, only visible inside for-loop
    x = z        # Using x from parent scope
    local y = z  # Local scope, only visible inside for-loop (not using y from parent scope)
end

println(x)    #  = 5, since for loop modifies parent variable x
println(y)    #  = 10, since for loop uses local variable y
println(z)    # Error: z only defined in the local scope of the for-loop


text:
import SimplePosets
using SimpleGraphs: SimpleGraph, AbstractSimpleGraph, SimpleDigraph, eltype, NV, elist, in_neighbors, add_edges!, vertex2idx

export toposort
function toposort(m::Model)
    names = toposort(poset(m).D)
    setdiff(names, arguments(m))
end

# Thanks to @CameronBieganek for this implementation:
# https://discourse.julialang.org/t/lightgraphs-jl-transition/69526/52?u=cscherrer
function toposort(g::SimpleDigraph{T}) where {T}
    g = deepcopy(g)
    order = T[]
    s = Symbol[v for v in vlist(g) if SimpleGraphs.in_deg(g, v) == 0]

    while !isempty(s)
        u = pop!(s)
        push!(order, u)

        for v in out_neighbors(g, u)
            SimpleGraphs.delete!(g, u, v)
            if SimpleGraphs.in_deg(g, v) == 0
                push!(s, v)
            end
        end
    end

    if SimpleGraphs.NE(g) > 0
        error("Graph contains cycles.")
    end

    order
end


text:
@testset "Spy" begin

    @mockable caller(x, y) = double(x) + double(y)
    @mockable double(x) = 2x
    @mockable triple(x) = 3x

    # basic statistics
    spy() do
        @test caller(1,2) == 6
        @test called_exactly_once(double, 1)
        @test called_exactly_once(double, 2)
        @test was_not_called(triple, 2)
    end

    # stats are gathered with `apply` as well
    apply() do
        @test caller(1,2) == 6
        @test called_exactly_once(double, 1)
        @test called_exactly_once(double, 2)
    end

    # stats for functions with kwargs
    @mockable div2(x; n = 2) = x / n
    spy() do
        @test div2(10) == 5
        @test div2(10; n = 5) == 2
        @test called_exactly_once(div2, 10; n = 2)
        @test called_exactly_once(div2, 10; n = 5)
    end
end


text:
module AutodetectDependabot

function _get_possible_branch_names()
    possible_branch_names = [
        get(ENV, "GITHUB_BASE_REF", ""),
        get(ENV, "GITHUB_HEAD_REF", ""),
        get(ENV, "GITHUB_REF", ""),
    ]
    return possible_branch_names
end

function _chop_refs_head(branch_name::AbstractString)
    replace(branch_name, r"^(refs\/heads\/)" => "")
end

function _is_dependabot_branch(branch_name::AbstractString)
    return startswith(branch_name, "dependabot/julia") || startswith(branch_name, "compathelper/")
end

function is_dependabot_job()
    possible_branch_names = _get_possible_branch_names()
    return any(_is_dependabot_branch.(_chop_refs_head.(possible_branch_names)))
end

end # module

max_stars_repo_name: fonsp/julia-runtest

max_issues_repo_name: snu-development/julia-runtest

max_forks_repo_name: snu-development/julia-runtest


text:
"""
This  package, which depends only  on the packages `LaurentPolynomials` and
`ModuleElts`,  implements Puiseux polynomials,  that is linear combinations
of  monomials of the  type `x₁^{a₁}… xₙ^{aₙ}`  where `xᵢ` are variables and
`aᵢ`  are exponents which can be arbitrary rational numbers.

Recall  that  on  an  algebraically  closed coefficient field, the quotient
field  of the Puiseux polynomials is  the algebraic closure of the quotient
field  of  the  multivariate  polynomials.  We use Puiseux polynomials with
cyclotomic coefficients as splitting fields of cyclotomic Hecke algebras.

This  package is in particular  a perfectly usable (and  quite good I hope)
implementation   of  multivariate  polynomials  and  multivariate  rational
fractions if you are only interested in that.

Some  functions described below work  only with polynomials where variables
are  raised to integral powers;  we will refer to  such objects as "Laurent
polynomials"; some functions require further that variables are raised only
to  positive powers: we refer then to "true polynomials".

Puiseux  polynomials have the  parametric type `Mvp{C,E}`  where `C` is the
type  of the coefficients and `E` is the type of the exponents (`E=Int` for
Laurent   polynomials;   `E=Rational{Int}`   for   more   general   Puiseux
polynomials).  When printing only `C` is  printed if `E==Int`. The rational
fractions  have type  `Frac{Mvp{C,Int}}`(and in  addition the numerator and
denominator should be true polynomials).

We first look at how to make Puiseux polynomials.

`@Mvp x₁,…,xₙ`

assigns  to each  Julia name  `xᵢ` an  `Mvp` representing  an indeterminate
suitable   to  build   multivariate  polynomials   or  rational  fractions.
`Mvp(:x₁)` creates the same `Mvp` without assigning it to variable `x₁`.

```julia-repl
julia> @Mvp x,y

julia> (x+y^-1)^3
Mvp{Int64}: x³+3x²y⁻¹+3xy⁻²+y⁻³

julia> x+Mvp(:z)
Mvp{Int64}: x+z

julia> x^(1//2)
Mvp{Int64,Rational{Int64}}: x½

julia> Mvp(3)
Mvp{Int64}: 3
```

`Mvp(x::Number)`  returns the multivariate polynomial  with only a constant
term, equal to `x`.

It  is convenient to create `Mvp`s using  variables such as `x,y` above. To
create  them more  directly, `Monomial(:x=>1,:y=>-2)`  creates the monomial
`xy⁻²`, and then `Mvp(Monomial(:x=>1,:y=>-2)=>3,Monomial()=>4)` creates the
`Mvp`  `3xy⁻²+4`. This is the way `Mvp` are printed in another context than
the repl, IJulia or Pluto (where they display nicely as shown above).

```julia-rep1
julia> print(3x*y^-2+4)
Mvp(Monomial(:x,:y => -2) => 3,Monomial() => 4) # :x is shorthand for :x=>1

julia> print(x^(1//2))
Mvp(Monomial(:x => 1//2) => 1)
```

Only  monomials and one-term `Mvp`s can  be raised to a non-integral power;
the  `Mvp` with one term constant `c`  times the monomial `m` can be raised
to  a fractional  power of  denominator `d`  if and  only if `root(c,d)` is
defined (this is equivalent to `c^{1//d}` for floats);

```julia-repl
julia> (4x)^(1//2)
Mvp{Int64,Rational{Int64}}: 2x½

julia> (2.0x)^(1//2)
Mvp{Float64,Rational{Int64}}: 1.4142135623730951x½

julia> root(2.0x)
Mvp{Float64,Rational{Int64}}: 1.4142135623730951x½
```

One  may  want  to  define  `root`  differently;  for instance, in my other
package   `CylotomicNumbers`  I   define  square   roots  of  rationals  as
cyclotomics;  I also have implemented in `CylotomicNumbers` arbitrary roots
of roots of unity). Using `CylotomicNumbers`:

```julia-rep1
julia> (2x)^(1//2)
Mvp{Cyc{Int64},Rational{Int64}}: √2x½

julia> (E(3)*x)^(2//3)
Mvp{Cyc{Int64},Rational{Int64}}: ζ₉²x⅔
```

There  are various ways to take an  `Mvp` apart. Below are the most direct;
look   also  at  the   functions  `coefficient`,  `coefficients`,  `pairs`,
`monomials`, `variables` and `powers`.

```julia-repl
julia> p=3x*y^-2+4
Mvp{Int64}: 3xy⁻²+4

julia> term(p,1)
xy⁻² => 3

julia> term(p,2)
 => 4

julia> length(p) # the number of terms
2

julia> term.(p,1:length(p)) # same as pairs(p)
2-element Vector{Pair{Monomial{Int64}, Int64}}:
 xy⁻² => 3
      => 4

julia> m=first(term(p,1)) # same as first(monomials(p))
Monomial{Int64}:xy⁻²

julia> length(m) # how many variables in m
2

julia> m[:x] # power of x in m
1

julia> m[:y] # power of y in m
-2

julia> map((x,y)->x=>y,variables(m),powers(m)) # same as pairs(m)
2-element Vector{Pair{Symbol, Int64}}:
 :x => 1
 :y => -2
```

The valuation and degree of an Mvp can be inspected globally or variable by
variable.

```julia-repl
julia> p
Mvp{Int64}: 3xy⁻²+4

julia> variables(p)
2-element Vector{Symbol}:
 :x
 :y

julia> degree(p),degree(p,:x),degree(p,:y)
(0, 1, 0)

julia> valuation(p),valuation(p,:x),valuation(p,:y)
(-1, 0, -2)
```

Terms  are totally ordered in an `Mvp`  by a monomial ordering (that is, an
ordering  on  monomials  so  that  `x<y`  implies `xz<yz` for any monomials
`x,y,z`).  By default, the  ordering is `lex`.  The terms are in decreasing
order,  so that the  first term is  the highest. The  orderings `grlex` and
`grevlex` are also implemented (see `grobner_basis` for how to use them).

An  `Mvp` is a *scalar*  if the valuation and  degree are `0`. The function
`scalar`  returns the  constant coefficient  if the  `Mvp` is a scalar, and
`nothing` otherwise.

Usual  arithmetic (`+`, `-`,  `*`, `^`, `/`,  `//`, `one`, `isone`, `zero`,
`iszero`,  `==`)  works.  Elements  of  type  `<:Number`  are considered as
scalars for scalar operations on the coefficients.

```julia-repl
julia> p
Mvp{Int64}: 3xy⁻²+4

julia> p^2
Mvp{Int64}: 9x²y⁻⁴+24xy⁻²+16

julia> p/2
Mvp{Float64}: 1.5xy⁻²+2.0

julia> p//2
Mvp{Rational{Int64}}: (3//2)xy⁻²+2//1
```

When  converting an `Mvp` one needs to specify the two type parameters (the
type of the coefficients and the type of the exponents).

```julia-repl
julia> Mvp{Float64,Rational{Int}}(p)
Mvp{Float64,Rational{Int64}}: 3.0xy⁻²+4.0
```

One can evaluate an `Mvp` when setting the value of some variables by using
the  function call syntax (actually, the keyword syntax for the object used
as a function). There is also a more explicit `value` function.

```julia-repl
julia> p=x+y
Mvp{Int64}: x+y

julia> value(p,:x=>2)
Mvp{Int64}: y+2

julia> p(x=2)
Mvp{Int64}: y+2

julia> p(x=2,y=x)
Mvp{Int64}: x+2

julia> value(p,:x=>2,:y=>x)
Mvp{Int64}: x+2
```

Note  that  an  `Mvp`  always  evaluates  to an `Mvp`, for consistency. You
should  use `scalar`  on the  result of  evaluating all  variables to get a
number.

```julia-repl
julia> p(x=1,y=2)
Mvp{Int64}: 3

julia> scalar(p(x=1,y=2))
3

julia> v=(x^(1//2))(x=2.0)
Mvp{Float64,Rational{Int64}}: 1.4142135623730951

julia> scalar(v)
1.4142135623730951
```

One  can divide an `Mvp` by another when the division is exact, compute the
`gcd` and `lcm` of two `Mvp`.

```julia-repl
julia> exactdiv(x^2-y^2,x-y)
Mvp{Int64}: x+y

julia> (x+y)/(2x^2)   # or by a monomial
Mvp{Float64}: 0.5x⁻¹+0.5x⁻²y

julia> (x+y)//(2x^2)
Mvp{Rational{Int64}}: (1//2)x⁻¹+(1//2)x⁻²y

julia> (x+y)/(x-y)    # otherwise one gets a rational fraction
Frac{Mvp{Int64, Int64}}: (x+y)/(x-y)
```

Raising  a non-monomial  Laurent polynomial  to a  negative power returns a
rational   fraction.  Rational  fractions  are  normalized  such  that  the
numerator  and denominators are true polynomials  prime to each other. They
have  the  arithmetic  operations  `+`,  `-`  , `*`, `/`, `//`, `^`, `inv`,
`one`,  `isone`, `zero`, `iszero` (these  operations can operate between an
`Mvp` or a `Number` and a `Frac{<:Mvp}`).

```julia-repl
julia> (x+1)^-2
Frac{Mvp{Int64, Int64}}: 1/(x²+2x+1)

julia> x+1/(y+1)
Frac{Mvp{Int64, Int64}}: (xy+x+1)/(y+1)

julia> 1-1/(y+1)
Frac{Mvp{Int64, Int64}}: y/(y+1)
```

One  can evaluate a `Frac` when setting  the value of some
variables by using the function call syntax or the `value` function:

```julia-repl
julia> ((x+y)/(x-y))(x=y+1)
Mvp{Float64}: 2.0y+1.0

julia> value((x+y)/(x-y),:x=>y+1;Rational=true) # use // for division
Mvp{Int64}: 2y+1
```

`Frac`  can  be  dissected  using  `numerator`  and `denominator`. They are
scalars  for  broadcasting  and  can  be  sorted  (have  `cmp` and `isless`
methods).

```julia-repl
julia> m=[x+y x-y;x+1 y+1]
2×2 Matrix{Mvp{Int64, Int64}}:
 x+y  x-y
 x+1  y+1

julia> n=inv(Frac.(m))
2×2 Matrix{Frac{Mvp{Int64, Int64}}}:
 (-y-1)/(x²-2xy-y²-2y)  (x-y)/(x²-2xy-y²-2y)
 (x+1)/(x²-2xy-y²-2y)   (-x-y)/(x²-2xy-y²-2y)

julia> lcm(denominator.(n))
Mvp{Int64}: x²-2xy-y²-2y
```

Finally,   `Mvp`s  have   methods  `conj`,   `adjoint`  which   operate  on
coefficients,   a   `derivative`   method,   and  methods  `positive_part`,
`negative_part` and `bar` (useful for Kazhdan-Lusztig theory).

```julia_repl
julia> @Mvp z
Mvp{Int64}: z

julia> hessian(p,vars)=[derivative(derivative(p,x),y) for x in vars, y in vars]
hessian (generic function with 1 method)

julia> hessian(x^2*y^2*z^2,[:x,:y,:z])
3×3 Matrix{Mvp{Int64, Int64}}:
 2y²z²  4xyz²  4xy²z
 4xyz²  2x²z²  4x²yz
 4xy²z  4x²yz  2x²y²

julia> jacobian(pols,vars)=[derivative(p,v) for p in pols, v in vars]
jacobian (generic function with 1 method)

julia> jacobian([x,y,z],[:x,:y,:z])
3×3 Matrix{Mvp{Int64, Int64}}:
 1  0  0
 0  1  0
 0  0  1

julia> p=(x+y^-1)^4
Mvp{Int64}: x⁴+4x³y⁻¹+6x²y⁻²+4xy⁻³+y⁻⁴

julia> positive_part(p)
Mvp{Int64}: x⁴

julia> negative_part(p)
Mvp{Int64}: y⁻⁴

julia> bar(p)
Mvp{Int64}: y⁴+4x⁻¹y³+6x⁻²y²+4x⁻³y+x⁻⁴
```

Despite  the degree of generality of our  polynomials, the speed is not too
shabby. For the Fateman test f(f+1) where f=(1+x+y+z+t)^15, we take 3sec.
According to the Nemo paper, Sagemath takes 10sec and Nemo takes 1.6sec.
"""
module PuiseuxPolynomials
using ModuleElts
using LaurentPolynomials
export coefficient, monomials, powers
export Mvp, Monomial, @Mvp, variables, value, laurent_denominator, term,
       lex, grlex, grevlex, grobner_basis, rename_variables
#------------------ Monomials ---------------------------------------------
"""
`Monomials` are implemented as a list of pairs `:variable=>degree` sorted
in the alphabetic order of vraiables.
"""
struct Monomial{T} # T is Int or Rational{Int}
  d::ModuleElt{Symbol,T}
end

Monomial(a::Union{Pair,Symbol}...)=Monomial(ModuleElt(map(x->x isa Symbol ?
          x=>1 : x,a)...)) # to be avoided when preformance needed
Monomial()=one(Monomial{Int})

Base.convert(::Type{Monomial{T}},v::Symbol) where T=Monomial(v=>T(1))
Base.convert(::Type{Monomial{T}},m::Monomial{N}) where {T,N}=
  Monomial(convert(ModuleElt{Symbol,T},m.d))

# we need a promote rule to mix fractional momonials with others in an Mvp
Base.promote_rule(::Type{Monomial{N}},::Type{Monomial{M}}) where {N,M}=
  Monomial{promote_type(N,M)}

Base.:*(a::Monomial, b::Monomial)=Monomial(a.d+b.d)
Base.isone(a::Monomial)=iszero(a.d)
#Base.iszero(a::Monomial)=false
Base.one(::Type{Monomial{T}}) where T=Monomial(zero(ModuleElt{Symbol,T}))
Base.one(m::Monomial)=Monomial(zero(m.d))
Base.inv(a::Monomial)=Monomial(-a.d)
LaurentPolynomials.exactdiv(a::Monomial, b::Monomial)=a*inv(b)
Base.:/(a::Monomial, b::Monomial)=a*inv(b)
Base.://(a::Monomial, b::Monomial)=a*inv(b)
Base.:^(x::Monomial,p)=Monomial(x.d*p)
Base.getindex(a::Monomial,k)=getindex(a.d,k)
Base.length(a::Monomial)=length(a.d)
"`variables(a::Monomial)` iterator on the variables of `a` (a sorted list of `Symbol`s)"
variables(a::Monomial)=keys(a.d)
"`powers(a::Monomial)` iterator on the powers of variables in `a`"
powers(a::Monomial)=values(a.d)
"""
`pairs(a::Monomial)`

returns the pairs `:variable=>power` in `a`.
"""
Base.pairs(a::Monomial)=pairs(a.d)
ispositive(a::Monomial)=all(>=(0),powers(a))

const unicodeFrac=Dict((1,2)=>'½',(1,3)=>'⅓',(2,3)=>'⅔',
  (1,4)=>'¼',(3,4)=>'¾',(1,5)=>'⅕',(2,5)=>'⅖',(3,5)=>'⅗',
  (4,5)=>'⅘',(1,6)=>'⅙',(5,6)=>'⅚',(1,8)=>'⅛',(3,8)=>'⅜',
  (5,8)=>'⅝',(7,8)=>'⅞',(1,9)=>'⅑',(1,10)=>'⅒',(1,7)=>'⅐')

const subvec=['₀','₁','₂','₃','₄','₅','₆','₇','₈','₉']

function LaurentPolynomials.stringexp(io::IO,r::Rational{<:Integer})
  d=denominator(r); n=numerator(r)
  if isone(d) return LaurentPolynomials.stringexp(io,n) end
  if get(io,:TeX,false) return "^{\\frac{$n}{$d}}" end
  res=Char[]
  if n<0 push!(res,'⁻'); n=-n end
  if haskey(unicodeFrac,(n,d)) push!(res,unicodeFrac[(n,d)])
  else
   if isone(n) push!(res,'\U215F')
   else append!(res,map(x->LaurentPolynomials.supvec[x+1],reverse(digits(n))))
     push!(res,'⁄')
   end
   append!(res,map(x->subvec[x+1],reverse(digits(d))))
  end
  String(res)
end

function Base.show(io::IO,m::Monomial)
  replorTeX=get(io,:TeX,false) || get(io,:limit,false)
  if !replorTeX
    print(io,"Monomial(")
    join(io,map(m.d)do (s,c)
      if c==1 repr(s)
      elseif c==2 join([repr(s),repr(s)],",")
      else repr(s=>c)
      end
    end,",")
    print(io,")")
    return
  end
  if isone(m) return end
  start=true
  for (v,d) in m.d
    if !(start || replorTeX) print(io,"*") end
    print(io,string(v))
    if !isone(d)
      if isone(denominator(d)) d=numerator(d) end
      if replorTeX print(io,LaurentPolynomials.stringexp(io,d))
      elseif d isa Integer print(io,"^$d")
      else print(io,"^($d)")
      end
    end
    start=false
  end
end

function Base.show(io::IO, ::MIME"text/plain", m::Monomial)
  if !haskey(io,:typeinfo)
    print(io,typeof(m),":")
    io=IOContext(io,:typeinfo=>typeof(m))
  end
  show(io,m)
end

"""
`lex(a::Monomial, b::Monomial)`
The  "lex" ordering,  where `a<b`  if the  first variable  in `a/b`
occurs to a positive power.
"""
function lex(a::Monomial, b::Monomial)
  for ((va,pa),(vb,pb)) in zip(a.d,b.d)
    if va!=vb return va<vb ? pa>0 : pb<0 end
    if pa!=pb return isless(pb,pa) end
  end
  la=length(a.d)
  lb=length(b.d)
  @inbounds la>lb ? last(a.d.d[lb+1])>0 : lb>la ? last(b.d.d[la+1])<0 : false
end

"""
`grlex(a::Monomial, b::Monomial)`
The "grlex" ordering, where `a<b̀` if `degree(a)>degree(b)` or the degrees
are equal but `lex(a,b)`.
"""
function grlex(a::Monomial, b::Monomial)
  da=degree(a);db=degree(b)
  if da!=db return isless(db,da) end
  lex(a,b)
end

"""
`grevlex(a::Monomial, b::Monomial)`
The "grevlex" ordering, where `a<b̀` if `degree(a)>degree(b)` or the degrees
are equal but the last variable in `a/b` occurs to a negative power.
"""
function grevlex(a::Monomial, b::Monomial)
  da=degree(a);db=degree(b)
  if da!=db return isless(db,da) end
  lex(Monomial(ModuleElt(reverse(b.d.d))),
      Monomial(ModuleElt(reverse(a.d.d))))
end

"""
`isless(a::Monomial,b::Monomial)`

For  our implementation of `Mvp`s to  work, `isless` must define a monomial
order (that is, for monomials `m,a,b` we have `a<b => a*m<b*m`). By default
we  use the  "lex" ordering.
"""
@inline Base.isless(a::Monomial, b::Monomial)=lex(a,b)

Base.:(==)(a::Monomial, b::Monomial)=a.d==b.d

#monomial d of greatest degree in each variable such that (m/d,n/d) positive
Base.gcd(m::Monomial,n::Monomial)=Monomial(ModuleElts.merge2(min,m.d,n.d))
Base.gcd(v::AbstractArray{<:Monomial})=reduce(gcd,v)

#monomial l of smallest degree in each variable such that (l/m,l/n) positive
Base.lcm(m::Monomial,n::Monomial)=Monomial(ModuleElts.merge2(max,m.d,n.d))
Base.lcm(v::AbstractArray{<:Monomial})=reduce(lcm,v)

Base.hash(a::Monomial, h::UInt)=hash(a.d,h)

LaurentPolynomials.degree(m::Monomial)=sum(powers(m);init=0)
LaurentPolynomials.degree(m::Monomial,var::Symbol)=m.d[var]

function LaurentPolynomials.root(m::Monomial,n::Integer=2)
  if all(x->iszero(x%n),powers(m))
       Monomial(ModuleElt(k=>div(v,n) for (k,v) in m.d))
  else Monomial(ModuleElt(k=>v//n for (k,v) in m.d))
  end
end

#------------------------------------------------------------------------------
"""
`Mvp`s are implemented as a list of pairs `monomial=>coefficient` sorted
by the monomial order `lex`.
"""
struct Mvp{T,N} # N=type of exponents T=type of coeffs
  d::ModuleElt{Monomial{N},T}
end

Mvp(a::Pair...;c...)=Mvp(ModuleElt(a...;c...))
Mvp(;c...)=zero(Mvp{Int,Int}) # for some calls to map() to work

"""
 `@Mvp x,y`

 is  equivalent to `x=Mvp(:x);y=Mvp(:y)`  excepted it creates  `x,y` in the
 global scope of the current module, since it uses `eval`.
"""
macro Mvp(t)
  if t isa Expr
    for v in t.args
      Base.eval(Main,:($v=Mvp($(Core.QuoteNode(Symbol(v))))))
    end
  elseif t isa Symbol
    Base.eval(Main,:($t=Mvp($(Core.QuoteNode(t)))))
  end
end

Base.broadcastable(p::Mvp)=Ref(p)
Base.cmp(a::Mvp,b::Mvp)=cmp(a.d,b.d)
Base.isless(a::Mvp,b::Mvp)=cmp(a,b)==-1
Base.hash(a::Mvp,h::UInt)=hash(a.d,h)

function Base.show(io::IO, ::MIME"text/html", a::Mvp)
  print(IOContext(io,:TeX=>true),"\$",a,"\$")
end

function Base.show(io::IO, ::MIME"text/plain", a::Mvp{T,N}) where{T,N}
  if !haskey(io,:typeinfo)
    print(io,N==Int ? "Mvp{$T}" : "Mvp{$T,$N}",": ")
    io=IOContext(io,:typeinfo=>typeof(a))
  end
  show(io,a)
end

function Base.show(io::IO, x::Mvp)
  if get(io,:limit,false) || get(io,:TeX,false)
    show(IOContext(io,:showbasis=>nothing),x.d)
    # :showbasis=>nothing necessary if called when already showing a ModuleElt
  else
    print(io,"Mvp(")
    join(io,repr.(pairs(x)),",")
    print(io,")")
  end
end

Base.length(x::Mvp)=length(x.d)
term(x::Mvp,i)=x.d.d[i]
ismonomial(x::Mvp)=length(x.d)==1 # x is a non-zero monomial
Base.zero(p::Mvp)=Mvp(zero(p.d))
Base.zero(::Type{Mvp{T,N}}) where {T,N}=Mvp(zero(ModuleElt{Monomial{N},T}))
Base.one(::Type{Mvp{T,N}}) where {T,N}=Mvp(one(Monomial{N})=>one(T);check=false)
Base.one(::Type{Mvp{T}}) where T=one(Mvp{T,Int})
Base.one(::Type{Mvp})=Mvp(1)
Base.one(p::Mvp{T,N}) where {T,N}=iszero(p) ? one(Mvp{T,N}) : Mvp(one(Monomial{N})=>one(first(coefficients(p)));check=false)
Base.isone(x::Mvp)=ismonomial(x) && isone(first(term(x,1)))&&isone(last(term(x,1)))
Base.copy(p::Mvp)=Mvp(p.d)
Base.iszero(p::Mvp)=iszero(p.d)
Base.convert(::Type{Mvp},a::Number)=convert(Mvp{typeof(a),Int},a)
Base.convert(::Type{Mvp{T,N}},a::Number) where {T,N}=iszero(a) ?
 zero(Mvp{T,N}) : Mvp(one(Monomial{N})=>convert(T,a);check=false)
(::Type{Mvp{T,N}})(a::Number) where {T,N}=convert(Mvp{T,N},a)
(::Type{Mvp{T,N}})(a::Mvp) where {T,N}=convert(Mvp{T,N},a)
Base.convert(::Type{Mvp{T,N}},a::Mvp{T1,N1}) where {T,T1,N,N1}=
  Mvp(convert(ModuleElt{Monomial{N},T},a.d))
Base.convert(::Type{Mvp},x::Mvp)=x
Base.convert(::Type{Mvp},v::Symbol)=Mvp(Monomial(v)=>1;check=false)
"""
`Mvp(x::Symbol)` creates the `Mvp` with one term of degree one and coefficient
1 with variable `x`
"""
Mvp(x::Symbol)=convert(Mvp,x)
Mvp(x::Number)=convert(Mvp,x)
Mvp(x::Mvp)=x
# stupid stuff to make LU work
Base.adjoint(a::Mvp)=conj(a)
Base.transpose(a::Mvp)=a
Base.abs(a::Mvp)=a

Base.:(==)(a::Mvp, b::Mvp)=a.d==b.d
Base.:(==)(a::Mvp,x::Number)=(iszero(a) && iszero(x)) || (!isnothing(x) && x==scalar(a))
Base.:(==)(x::Number,a::Mvp)= a==x

function Base.convert(::Type{T},a::Mvp) where T<:Number
  if iszero(a) return zero(T) end
  x=scalar(a)
  if isnothing(x) throw(InexactError(:convert,T,a)) end
  T(x)
end
(::Type{T})(a::Mvp) where T<: Number=convert(T,a)

Base.isinteger(p::Mvp)=scalar(p)!==nothing && isinteger(scalar(p))

# we need a promote rule to handle Vectors of Mvps of different types
Base.promote_rule(::Type{Mvp{T1,N1}},::Type{Mvp{T2,N2}}) where {T1,T2,N1,N2} =
  Mvp{promote_type(T1,T2),promote_type(N1,N2)}

# and this rule to handle Vectors mixing Mvps with numbers
Base.promote_rule(::Type{Mvp{T1,N}},::Type{T2}) where {T1,N,T2<:Number} =
  Mvp{promote_type(T1,T2),N}

Base.:+(a::Mvp, b::Mvp)=Mvp(a.d+b.d) # ModuleElts.+ takes care of promotion
Base.:+(a::Number, b::Mvp)=+(promote(a,b)...)
Base.:+(a::Mvp, b::Number)=+(promote(a,b)...)

Base.:-(a::Mvp)=Mvp(-a.d)
Base.:-(a::Mvp, b::Mvp)=Mvp(a.d-b.d)
Base.:-(a::Mvp, b::Number)=-(promote(a,b)...)
Base.:-(b::Number, a::Mvp)=-(promote(b,a)...)

Base.:*(a::Number, b::Mvp)=Mvp(b.d*a)
Base.:*(b::Mvp, a::Number)=a*b
# we have a monomial order so there is no ordering check in next line
Base.:*(a::Monomial, b::Mvp)=Mvp(ModuleElt(m*a=>c for (m,c) in pairs(b);check=false))
Base.:*(b::Mvp,a::Monomial)=a*b

function Base.:*(a::Mvp, b::Mvp)
  if length(a)>length(b) a,b=(b,a) end
  a,b=promote(a,b)
  if iszero(a) return a
  elseif iszero(b) return b
  elseif isone(a) return b
  elseif isone(b) return a
  end
  let b=b # needed !!!!
    sum(b*m*c for (m,c) in pairs(a))
  end
end

Base.:/(p::Mvp,q::Number)=Mvp(p.d/q)
Base.://(p::Mvp,q::Number)=Mvp(p.d//q)
Base.div(a::Mvp,b::Number)=Mvp(merge(div,a.d,b))
LaurentPolynomials.exactdiv(a::Mvp,b::Number)=Mvp(merge(exactdiv,a.d,b;check=false))

"""
`conj(p::Mvp)` acts on the coefficients of `p`

```julia-repl
julia> @Mvp x;conj(im*x+1)
Mvp{Complex{Int64}}: (0 - 1im)x+1 + 0im
```
"""
Base.conj(a::Mvp)=Mvp(merge(conj,a.d;check=false))

function Base.:^(x::Mvp, p::Union{Integer,Rational})
  if isinteger(p) p=Int(p) end
  if iszero(p) return one(x)
  elseif iszero(x) || isone(p) return x
  elseif p isa Rational return root(x,denominator(p))^numerator(p)
  elseif ismonomial(x)
    (m,c)=term(x,1);return isone(c) ? Mvp(m^p=>c;check=false) : Mvp(m^p=>c^p;check=false)
  elseif p>=0 return Base.power_by_squaring(x,p)
  else return Base.power_by_squaring(inv(x),-p)
  end
end

"""
`degree(m::Mvp[,v::Symbol])`

The `degree` of a monomial is the sum of  the exponents of the variables.
The `degree` of an `Mvp` is the largest degree of a monomial.

With  second argument a  variable name, `degree`  returns the degree of the
polynomial in that variable.

```julia-repl
julia> a=x^2+x*y
Mvp{Int64}: x²+xy

julia> degree(a), degree(a,:y), degree(a,:x)
(2, 1, 2)
```
"""
LaurentPolynomials.degree(p::Mvp)=iszero(p) ? 0 : maximum(degree,monomials(p))
LaurentPolynomials.degree(p::Mvp,v::Symbol)=iszero(p) ? 0 : maximum(x->degree(x,v),monomials(p))

"""
`valuation(m::Mvp[,v::Symbol])`

The `valuation` of an `Mvp` is the minimal degree of a monomial.

With  second argument a variable name, `valuation` returns the valuation of
the polynomial in that variable.

```julia-repl
julia> @Mvp x,y; a=x^2+x*y
Mvp{Int64}: x²+xy

julia> valuation(a), valuation(a,:y), valuation(a,:x)
(2, 0, 1)
```

"""
LaurentPolynomials.valuation(p::Mvp)=iszero(p) ? 0 : minimum(degree,monomials(p))
LaurentPolynomials.valuation(p::Mvp,v::Symbol)=iszero(p) ? 0 : minimum(x->degree(x,v),monomials(p))

"""
`coefficient(p::Mvp,m::Monomial)`

The coefficient of the polynomial `p` on the monomial `m`.

```julia-repl
julia> @Mvp x,y; p=(x-y)^3
Mvp{Int64}: x³-3x²y+3xy²-y³

julia> coefficient(p,Monomial(:x,:x,:y)) # coefficient on x²y
-3

julia> coefficient(p,Monomial()) # constant coefficient
0
```
"""
coefficient(p::Mvp,m::Monomial)=p.d[m]

"""
`pairs(p::Mvp)`

returns the pairs monomial=>coefficient in `p`.
"""
Base.pairs(p::Mvp)=pairs(p.d)

"""
`monomials(p::Mvp)`

is an efficient iterator over the monomials of `p`.
"""
monomials(p::Mvp)=keys(p.d)

"""
`coefficients(p::Mvp)`

is an efficient iterator over the coefficients of the monomials in `p`.
"""
LaurentPolynomials.coefficients(p::Mvp)=values(p.d)

"""
  `coefficients(p::Mvp, var::Symbol)`

returns  a Dict with keys the degree  in `var` and values the corresponding
coefficient of `p` with respect to `var`.

```julia-repl
julia> p=(x+y+inv(y))^4
Mvp{Int64}: x⁴+4x³y+4x³y⁻¹+6x²y²+12x²+6x²y⁻²+4xy³+12xy+12xy⁻¹+4xy⁻³+y⁴+4y²+6+4y⁻²+y⁻⁴

julia> coefficients(p,:x)
Dict{Int64, Mvp{Int64, Int64}} with 5 entries:
  0 => y⁴+4y²+6+4y⁻²+y⁻⁴
  4 => 1
  2 => 6y²+12+6y⁻²
  3 => 4y+4y⁻¹
  1 => 4y³+12y+12y⁻¹+4y⁻³

julia> coefficients(p,:y)
Dict{Int64, Mvp{Int64, Int64}} with 9 entries:
  0  => x⁴+12x²+6
  4  => 1
  -1 => 4x³+12x
  2  => 6x²+4
  -3 => 4x
  -2 => 6x²+4
  -4 => 1
  3  => 4x
  1  => 4x³+12x
```

The  same  caveat  is  applicable  to  `coefficients` as to evaluating: the
values  are always `Mvp`s. To get a list of scalars for the coefficients of
a  univariate polynomial represented as a `Mvp`, one should use `scalar` on
the values of `coefficients`.
"""
function LaurentPolynomials.coefficients(p::Mvp{T,N},v::Symbol)where {T,N}
  if iszero(p) return Dict{Int,typeof(p)}() end
  d=Dict{N,typeof(p.d.d)}()
  for (m,c) in pairs(p)
    found=false
    for (i,(v1,deg)) in enumerate(m.d)
      if v1==v
        found=true
        d[deg]=push!(get!(d,deg,empty(p.d.d)),
                Monomial(ModuleElt(deleteat!(copy(m.d.d),i);check=false))=>c)
      end
    end
    if !found  d[0]=push!(get!(d,0,empty(p.d.d)),m=>c) end
  end
  Dict(dg=>Mvp(c...;check=false) for (dg,c) in d) # c... is sorted by defn of monomial order
end

"""
  `coefficient(p::Mvp, var::Symbol, d)`

returns  the coefficient of degree `d` in the variable `var` in the `Mvp` `p`.

```julia-repl
julia> @Mvp x,y; p=(x+y^(1//2)+1)^3
Mvp{Int64,Rational{Int64}}: x³+3x²y½+3x²+3xy+6xy½+3x+y³⁄₂+3y+3y½+1

julia> coefficient(p,:y,1//2)
Mvp{Int64,Rational{Int64}}: 3x²+6x+3

julia> coefficient(p,:x,1)
Mvp{Int64,Rational{Int64}}: 3y+6y½+3
```
"""
function coefficient(p::Mvp,var::Symbol,d)
  res=empty(p.d.d)
  for (m,c) in pairs(p)
    found=false
    for (i,(v1,deg)) in enumerate(m.d)
      if v1==var && deg==d
        found=true
        push!(res,Monomial(ModuleElt(deleteat!(copy(m.d.d),i);check=false))=>c)
        break
      elseif v1>var break
      end
    end
    if d==0 && !found push!(res,m=>c) end
  end
  Mvp(ModuleElt(res;check=false))
end

"""
  `Pol(p::Mvp{T}) where T`

  converts the one-variable `Mvp{T}` `p` to a `Pol{T}`. It is an error if `p`
  has more than one variable.

```julia-repl
julia> @Mvp x; @Pol q; Pol(x^2+x)
Pol{Int64}: q²+q
```
"""
function LaurentPolynomials.Pol(x::Mvp{T})where T
  l=variables(x)
  if isempty(l) return Pol(scalar(x)) end
  if length(l)>1 error("cannot convert $(length(l))-variate Mvp to Pol") end
  val=Int(valuation(x))
  p=zeros(T,Int(degree(x))-val+1)
  for (mon,coeff) in pairs(x) p[Int(degree(mon))-val+1]=coeff end
  Pol(p,val)
end

"""
  `Pol(p::Mvp,v::Symbol)`

returns  a polynomial whose coefficients are  the coefficients of the `Mvp`
`p`  with respect to the variable `v`.  The variable `v` should appear only
with integral powers in `p`.

```julia-repl
julia> p=(x+y^(1//2))^3
Mvp{Int64,Rational{Int64}}: x³+3x²y½+3xy+y³⁄₂

julia> Pol(:q); Pol(p,:x)
Pol{Mvp{Int64, Rational{Int64}}}: q³+3y½q²+3yq+y³⁄₂
```
"""
function LaurentPolynomials.Pol(p::Mvp{T,N},var::Symbol)where{T,N}
  v=Int(valuation(p,var))
  res=[Pair{Monomial{N},T}[] for i in v:Int(degree(p,var))]
  for (m,c) in pairs(p)
    found=false
    for (i,(v1,deg)) in enumerate(m.d)
      if v1==var
        found=true
        push!(res[Int(deg)-v+1],
             Monomial(ModuleElt(deleteat!(copy(m.d.d),i);check=false))=>c)
        break
      elseif v1>var break
      end
    end
    if !found push!(res[-v+1],m=>c) end
  end
  Pol(map(x->Mvp(ModuleElt(x;check=false)),res),v)
end

"""
`Mvp(p::Pol)` converts `p` to  an  `Mvp` (with the same variable name).

```julia-repl
julia> @Pol q
Pol{Int64}: q

julia> Mvp(q^2+q)
Mvp{Int64}: q²+q
```
"""
Mvp(x::Pol)=convert(Mvp,x)

Base.convert(::Type{Mvp{T,N}},p::Pol) where{T,N}=
                     p(Mvp(convert(Monomial{N},LaurentPolynomials.varname[])=>one(T)))
Base.convert(::Type{Mvp},p::Pol)=p(Mvp(LaurentPolynomials.varname[]))

"""
`variables(p::Mvp)`

`variables(v::Array{Mvp})`

returns  the list of  variables of `p`  (resp. all `p`  in `v`) as a sorted
list of `Symbol`s.

```julia-repl
julia> @Mvp x,y,z

julia> variables([x+y+1,z])
3-element Vector{Symbol}:
 :x
 :y
 :z
```
"""
variables(pp::AbstractArray{<:Mvp})=sort(unique(k1 for p in pp for k in monomials(p) for k1 in variables(k)))

variables(p::Mvp)=sort(unique(k1 for k in monomials(p) for k1 in variables(k)))

"""
`scalar(p::Mvp)`

If `p` is a scalar, return that scalar, otherwise return `nothing`.

```julia-repl
julia> p=Mvp(:x)+1
Mvp{Int64}: x+1

julia> w=p(x=4)
Mvp{Int64}: 5

julia> scalar(w)
5

julia> typeof(scalar(w))
Int64
```
if  `p` is an array, then apply `scalar` to its elements and return `nothing`
if it contains any `Mvp` which is not a scalar.
"""
function LaurentPolynomials.scalar(p::Mvp{T})where T
  if iszero(p) return zero(T) end
  if ismonomial(p)
    (m,c)=term(p,1)
    if isone(m) return c end
  end
end

function LaurentPolynomials.scalar(m::AbstractArray{<:Mvp})
  p=scalar.(m)
  if !any(isnothing,p) return p end
end

"""
    `value(p::Mvp,:x₁=>v₁,:x₂=>v₂,...)`
     ̀(p::Mvp)(x₁=v₁,…,xₙ=vₙ)`

returns  the value of  `p` when doing  the simultaneous substitution of the
variable `:x1` by `v1`, of `x2` by `v2`, …

```julia-repl
julia> p=-2+7x^5*inv(y)
Mvp{Int64}: 7x⁵y⁻¹-2

julia> p(x=2)
Mvp{Int64}: -2+224y⁻¹

julia> p(y=1)
Mvp{Int64}: 7x⁵-2

julia> p(x=2,y=1)
Mvp{Int64}: 222
```

One should pay attention to the fact that the last value is not an integer,
but  a constant `Mvp` (for consistency).  See the function `scalar` for how
to convert such constants to their base ring.

```julia-repl
julia> p(x=y)
Mvp{Int64}: 7y⁴-2

julia> p(x=y,y=x)
Mvp{Int64}: -2+7x⁻¹y⁵
```
Evaluating an `Mvp` which is a Puiseux polynomial may cause calls to `root`

```julia-repl
julia> p=x^(1//2)*y^(1//3)
Mvp{Int64,Rational{Int64}}: x½y⅓

julia> p(;x=y)
Mvp{Int64,Rational{Int64}}: y⅚

julia> p(;x=4)
Mvp{Int64,Rational{Int64}}: 2y⅓

julia> p(;y=2.0)
Mvp{Float64,Rational{Int64}}: 1.2599210498948732x½
```
"""
function value(p::Mvp,k::Pair...)
  res1=badi=nothing
  if length(k)==1 (key,val)=k[1]
  # special code when length(k)==1 and !(key∈ variables(p)) is 10x faster
  elseif isempty(k) return p
  else vv=Dict(k)
  end
  for (i,(m,c1)) in enumerate(pairs(p))
    res=badj=nothing
    for (j,(v,c)) in enumerate(m.d)
      if length(k)==1 && key!=v continue end
      if length(k)>1
        if !haskey(vv,v) continue
        else val=vv[v]
        end
      end
      if badj==nothing
        badj=Int[]
        res=Mvp(c1)
      end
      if isinteger(c) res*=val^Int(c)
      else res*=root(val,denominator(c))^numerator(c)
      end
      push!(badj,j)
    end
    if badj!==nothing
      res*=Monomial(ModuleElt(deleteat!(copy(m.d.d),badj);check=false))
      badj=nothing
      if badi===nothing res1=res
        badi=Int[]
      else res1+=res
      end
      push!(badi,i)
    end
 #  println("badi=$badi m=$m c=$c res1=$res1")
  end
  if badi===nothing return p end
  Mvp(deleteat!(copy(pairs(p)),badi)...;check=false)+res1
end

(p::Mvp)(;arg...)=value(p,arg...)

function LaurentPolynomials.root(a::Integer,n=2)
  if n!=2 return end
  r=isqrt(a)
  if r^2==a return r end
end

function LaurentPolynomials.root(p::Mvp,n=2)
  if iszero(p) return p end
  n=Int(n)
  if !ismonomial(p)
    throw(DomainError("root($p,$n) non-monomial not implemented"))
  end
  (m,c)=term(p,1)
  isone(c) ? Mvp(root(m,n)=>c;check=false) : Mvp(root(m,n)=>root(c,n))
end

"""
`Base.:^(p,m;vars=variables(p))`

Implements  the action of  a matrix on  `Mvp`s. `vars` should  be a list of
symbols   representing  variables.   The  polynomial   `p`  is  changed  by
simultaneous  substitution in it of  `vᵢ` by `(v×m)ᵢ` where  `v` is the row
vector  of  the  `Mvp(vᵢ)`.  If  `vars`  is  omitted,  it  is  taken  to be
`variables(p)`.

```julia-repl
julia> @Mvp x,y

julia> (x+y)^[1 2;3 1]
Mvp{Int64}: 3x+4y
```
"""
Base.:^(p::Mvp,m::AbstractMatrix;vars=variables(p))=
  p(;map(Pair,vars,permutedims(Mvp.(vars))*m)...)

LaurentPolynomials.positive_part(p::Mvp)=
  Mvp(ModuleElt(m=>c for (m,c) in pairs(p) if ispositive(m);check=false))

LaurentPolynomials.negative_part(p::Mvp)=
  Mvp(ModuleElt(m=>c for (m,c) in pairs(p) if all(<(0),powers(m));check=false))

LaurentPolynomials.bar(p::Mvp)=Mvp(ModuleElt(inv(m)=>c for (m,c) in pairs(p)))

"""
The  function 'derivative(p,v₁,…,vₙ)' returns the  derivative of 'p' with
respect to  the variable given by the symbol 'v₁', then `v₂`, ...

```julia-repl
julia> @Mvp x,y;p=7x^5*y^-1-2
Mvp{Int64}: 7x⁵y⁻¹-2

julia> derivative(p,:x)
Mvp{Int64}: 35x⁴y⁻¹

julia> derivative(p,:y)
Mvp{Int64}: -7x⁵y⁻²

julia> derivative(p,:x,:y)
Mvp{Int64}: -35x⁴y⁻²

julia> p=x^(1//2)*y^(1//3)
Mvp{Int64,Rational{Int64}}: x½y⅓

julia> derivative(p,:x)
Mvp{Rational{Int64},Rational{Int64}}: (1//2)x⁻½y⅓

julia> derivative(p,:y)
Mvp{Rational{Int64},Rational{Int64}}: (1//3)x½y⁻⅔

julia> derivative(p,:z)
Mvp{Rational{Int64},Rational{Int64}}: 0
```
"""
function LaurentPolynomials.derivative(p::Mvp,vv...)
  # check needed because 0 could appear in coeffs
  for v in vv
    p=Mvp(ModuleElt(m*Monomial(v=>-1)=>c*degree(m,v) for (m,c) in pairs(p)))
  end
  p
end

# returns p/q when the division is exact, nothing otherwise
# Arguments must be true polynomials
function LaurentPolynomials.exactdiv(p::Mvp,q::Mvp)
  if isone(q) return p end
  if iszero(q) error("cannot divide by 0")
  elseif iszero(p) || isone(q) return p
  elseif ismonomial(q)
    m,c=term(q,1)
    return Mvp(ModuleElt(inv(m)*m1=>exactdiv(c1,c) for (m1,c1) in pairs(p);check=false))
   elseif ismonomial(p) error(q," does not exactly divide ",p)
  end
  var=first(variables(first(monomials(p))))
  res=zero(p)
  mq=degree(q,var)
  cq=coefficient(q,var,mq)
  while !iszero(p)
#   if length(p)<length(q) return nothing end
    mp=degree(p,var)
    t=exactdiv(coefficient(p,var,mp),cq)
    if mp!=mq t=Monomial(ModuleElt(var=>mp-mq))*t end
    res+=t
    p-=t*q
  end
  res
end

Base.gcd(a::AbstractFloat,b::AbstractFloat)=one(a)*one(b)

"""
`gcd(p::Mvp,  q::Mvp)`  computes  the  `gcd`  of  the  'Mvp' arguments. The
arguments must be true polynomials.

```julia-repl
julia> gcd(x^2-y^2,(x+y)^2)
Mvp{Int64}: -x-y
```
"""
function Base.gcd(a::Mvp,b::Mvp)
  if iszero(a) return b
  elseif iszero(b) return a
  elseif ismonomial(a)
    (m,c)=term(a,1)
    return Mvp(gcd(m,reduce(gcd,monomials(b)))=>gcd(c,reduce(gcd,coefficients(b)));check=false)
  elseif ismonomial(b) return gcd(b,a)
  end
  va=variables(a)
  vb=variables(b)
  vars=intersect(va,vb)
  if isempty(vars)
    if iszero(a) return b
    elseif iszero(b) return a
    else
      return Mvp(gcd(reduce(gcd,coefficients(a)),reduce(gcd,coefficients(b))))
    end
  end
  v=first(vars)
  if length(union(va,vb))==1 srgcd(Pol(a),Pol(b))(Mvp(v))
  else              srgcd(Pol(a,v),Pol(b,v))(Mvp(v))
  end
end

Base.gcd(v::AbstractArray{<:Mvp})=reduce(gcd,v;init=Mvp(0))

"""
`laurent_denominator(p1,p2,…)`

returns  the unique true monomial  `m` of minimal degree  such that for all
the Laurent polynomials `p1,p2,…` the product `m*pᵢ` is a true polynomial.

```julia-repl
julia> laurent_denominator(x^-1,y^-2+x^4)
Monomial{Int64}:xy²
```
"""
laurent_denominator(pp::Mvp...)=inv(gcd([m for p in pp for (m,c) in pairs(p)]))

"""
`lcm(p1,p2,...)`

Returns  the Lcm  of the  `Mvp` arguments.  The arguments  must be  true
polynomials.

```julia-repl
julia> lcm(x^2-y^2,(x+y)^2)
Mvp{Int64}: -x³-x²y+xy²+y³
```
"""
Base.lcm(a::Mvp,b::Mvp)=exactdiv(a*b,gcd(a,b))
Base.lcm(a::AbstractArray{<:Mvp})=reduce(lcm,a)

Base.eltype(p::Mvp{T,N}) where{T,N} =T

Base.denominator(p::Mvp)=lcm(denominator.(coefficients(p)))
Base.numerator(p::Mvp{<:Rational{T},N}) where{T,N} =convert(Mvp{T,N},p*denominator(p))
#----------------------------- Frac{Mvp{T,Int}} -------------------------------
# make both pols positive without common monomial factor
function make_positive(a::Mvp,b::Mvp)
  d=laurent_denominator(a,b)
  isone(d) ? (a,b) : (a*d,b*d)
end

LaurentPolynomials.Frac(a::Mvp,b::Mvp;k...)=Frac(promote(a,b)...;k...)

"""
`Frac(a::Mvp,b::Mvp;pol=false,prime=false)`

`Mvp`s  `a` and `b` are promoted to  same coefficient type, and checked for
being  true polynomials  without common  monomial factor (unless `pol=true`
asserts  that this  is already  the case)  and unless `prime=true` they are
made prime to each other by dividing by their gcd.
"""
function LaurentPolynomials.Frac(a::T,b::T;pol=false,prime=false)::Frac{T} where T<:Mvp
  if iszero(a) return LaurentPolynomials.Frac_(a,one(a))
  elseif iszero(b) error("division by 0")
  end
  if !pol
    (a,b)=make_positive(a,b)
  end
  if !prime && !ismonomial(b) && !ismonomial(a)
    d=gcd(a,b)
    a,b=exactdiv(a,d),exactdiv(b,d)
  end
  if scalar(b)==-1 a,b=(-a,-b) end
  return LaurentPolynomials.Frac_(a,b)
end

LaurentPolynomials.Frac(a::Mvp)=Frac(a,one(a);prime=true)

function LaurentPolynomials.Frac(a::Mvp{<:Rational,Int},b::Mvp{<:Rational,Int};k...)
  Frac(numerator(a)*denominator(b),numerator(b)*denominator(a);k...)
end

function Mvp(p::Frac{<:Mvp})
  if isone(p.den) return p.num
  elseif ismonomial(p.den)
    (m,c)=term(p.den,1)
    return p.num*inv(m)*inv(c)
  end
  error("cannot convert ",p," to Mvp")
end

Base.convert(::Type{Mvp{T,Int}},p::Frac{<:Mvp}) where T=convert(Mvp{T,Int},Mvp(p))

function Base.convert(::Type{Frac{T}},p::Mvp) where T<:Mvp
  q=convert(T,p)
  Frac(q,one(q);prime=true)
end

function Base.convert(::Type{Frac{T}},p::Number) where T<:Mvp
  Frac(convert(T,p),convert(T,1);pol=true,prime=true)
end

function Base.promote_rule(a::Type{T1},b::Type{Frac{T2}})where {T1<:Mvp,T2<:Mvp}
  Frac{promote_type(T1,T2)}
end

function Base.promote_rule(a::Type{Frac{T1}},b::Type{Frac{T2}})where {T1<:Mvp,T2<:Mvp}
  Frac{promote_type(T1,T2)}
end

function Base.promote_rule(a::Type{T1},b::Type{Frac{T2}})where {T1<:Number,T2<:Mvp}
  Frac{promote_type(Mvp{T1,Int},T2)}
end

function Base.inv(p::Mvp)
  if ismonomial(p)
    (m,c)=term(p,1)
    return Mvp(inv(m)=>c^2==1 ? c : 1/c)
  end
  Frac(Mvp(1),p;prime=true)
end

function Base.://(a::Mvp,b::Mvp)
  if iszero(a) return a end
  if ismonomial(b)
    (m,c)=term(b,1)
    return Mvp(ModuleElt(m1/m=>c^2==1 ? c1*c : c1//c for (m1,c1) in pairs(a);check=false))
  end
  Frac(a,b)
end

function Base.:/(a::Mvp,b::Mvp)
  if iszero(a) return a end
  if ismonomial(b)
    (m,c)=term(b,1)
    return Mvp(ModuleElt(m1/m=>c1/c for (m1,c1) in pairs(a);check=false))
  end
  Frac(a,b)
end

Base.://(a::Frac{<:Mvp},b::Number)=Frac(a.num,a.den*b;pol=true,prime=true)
Base.://(a::Frac{<:Mvp},b::Mvp)=Frac(a.num,a.den*b)
Base.:/(a::Frac{<:Mvp},b::Union{Frac{<:Mvp},Number,Mvp})=a//b
Base.:/(a::Union{Number,Mvp},b::Frac{<:Mvp})=a*inv(b)
Base.://(p::Number,q::Mvp)=Mvp(p)//q
Base.:/(p::Number,q::Mvp)=Mvp(p)/q

Base.:*(a::Frac{<:Mvp},b::Number)=Frac(a.num*b,a.den;pol=true,prime=true)
Base.:*(b::Number,a::Frac{<:Mvp})=a*b

value(p::Frac{<:Mvp},k::Pair...;Rational=false)=Rational ?
  value(p.num,k...)//value(p.den,k...) : value(p.num,k...)/value(p.den,k...)

(p::Frac{<:Mvp})(;arg...)=value(p,arg...)

#--------------------  Grobner bases -------------------------------------
# the reference is Cox, Little, O'Shea chapter 2

# minimum term and its index for monomial order lt
# findmin does not have the keywords lt and by so I must spin my own
function fmin(l;lt=lex,by=first)
  res=(first(l),1)
  for i in 2:length(l)
   if lt(by(l[i]),by(first(res))) res=(l[i],i) end
  end
  res
end

# Leading term for monomial order lt
LT(p;lt=lex)=lt==lex ? first(pairs(p)) : fmin(pairs(p);lt)[1]

# drop from p leading term for monomial order lt
function dropLT(p;lt=lex)
  if lt==lex return Mvp(ModuleElt(pairs(p)[2:end];check=false)) end
  Mvp(ModuleElt(deleteat!(pairs(p),fmin(pairs(p);lt)[2]);check=false))
end

# quotient of leading terms
function quotientLT(p,q;lt=lex)
  pm,pc=LT(p;lt)
  qm,qc=LT(q;lt)
  t=pm/qm
  if ispositive(t) Mvp(t=>exactdiv(pc,qc)) end
end

# remainder on division of p by list F
# Cox-Little-O'Shea Th. 3 §3 chap.2
function remainder(p,F;lt=lex)
  q=zero(F)
  r=zero(p)
  while !iszero(p)
    gotquotient=false
    for i in eachindex(F)
      t=quotientLT(p,F[i];lt)
      if t!==nothing
        q[i]+=t
#       xprintln("p=",p," F[i]=",F[i]," t=",t)
        p-=t*F[i]
        gotquotient=true
        break
      end
    end
    if !gotquotient
      r+=Mvp(LT(p;lt))
      p=dropLT(p;lt)
    end
  end
  (q,r)
end

# Cox-Little-O'Shea def. 4.(ii) §6 chap.2
function S_polynomial(p,q;lt=lex)
  pm,pc=LT(p;lt)
  qm,qc=LT(q;lt)
  c=lcm(pm,qm)
  (c/pm)*exactdiv(p,pc)-(c/qm)*exactdiv(q,qc)
end

function reduce_basis(F;lt=lex)
# F=sort(F,by=length)
  F=copy(F)
  i=1
  while i<=length(F)
    if any(j->j!=i && quotientLT(F[i],F[j];lt)!==nothing,eachindex(F))
      deleteat!(F,i)
    else i+=1
    end
  end
  F
end

"""
`grobner_basis(F;lt=lex)`

computes  a Gröbner basis  of the polynomial  ideal generated by the `Mvp`s
given  by the vector `F`. The  keyword `lt` describes the monomial ordering
to use.
```julia-repl
julia> @Mvp x,y,z; F=[x^2+y^2+z^2-1,x^2-y+z^2,x-z]
3-element Vector{Mvp{Int64, Int64}}:
 x²+y²+z²-1
 x²-y+z²
 x-z

julia> grobner_basis(F)
3-element Vector{Mvp{Int64, Int64}}:
 x-z
 -y+2z²
 4z⁴+2z²-1

julia> grobner_basis(F;lt=grlex)
3-element Vector{Mvp{Int64, Int64}}:
 x-z
 y²+y-1
 -y+2z²

julia> grobner_basis(F;lt=grevlex)
3-element Vector{Mvp{Int64, Int64}}:
 x-z
 y²+y-1
 2x²-y
```

There is no keyword to change the ordering of the variables. We suggest
to use `rename_variables` for this purpose.
"""
function grobner_basis(F;lt=lex)
# Cox-Little-O'Shea Th. 9 §10 chap.2
  F=copy(F)
  B=[(i,j) for j in 1:length(F) for i in 1:j-1]
  t=length(F)
  while !isempty(B)
    i,j=popfirst!(B)
    li,_=LT(F[i];lt)
    lj,_=LT(F[j];lt)
    lij=lcm(li,lj)
    if lij==li*lj continue end
    ll=filter(l->l!=i && l!=j && !((i,l) in B) && !((j,l) in B)
               && !((l,i) in B) && !((l,j) in B),1:length(F))
    if any(l->ispositive(lij/first(LT(F[l];lt))),ll) continue end
    s=S_polynomial(F[i],F[j];lt)
    r=remainder(s,F;lt)[2]
    if !iszero(r)
      t+=1
      push!(F,r)
      append!(B,tuple.(1:t-1,t))
    end
  end
  reduce_basis(F;lt)
end

"""
`rename_variables(p,v)` renames `variables(p)` to `v`

`rename_variables(p,l::Pair{Symbol,Symbol}...)` renames the
variables in `p` as indicated by the pairs in `l`.

```julia-repl
julia> @Mvp x,y,z; p=x+y+z
Mvp{Int64}: x+y+z

julia> rename_variables(p,Symbol.('A':'Z'))
Mvp{Int64}: A+B+C

julia> rename_variables(p,[:U,:V])
Mvp{Int64}: U+V+z

julia> rename_variables(p,:x=>:U,:z=>:V) # faster than p(;x=Mvp(:U),z=Mvp(:V))
Mvp{Int64}: U+V+y
```
"""
function rename_variables(p::Mvp,l::Pair{Symbol,Symbol}...)
  d=Dict(l)
  Mvp(ModuleElt(map(pairs(p)) do (m,c)
    Monomial(ModuleElt(map(((v,i),)->(haskey(d,v) ? d[v] : v)=>i,m.d.d)))=>c
  end))
end

rename_variables(p::Mvp,l::AbstractVector{Symbol})=
     rename_variables(p,map((x,y)->x=>y,variables(p),l)...)
#--------------------  benchmarks -------------------------------------
#julia> @btime (x+y+z)^3;
#  2.292 μs (49 allocations: 5.42 KiB)

#julia1.6.3> @btime PuiseuxPolynomials.fateman(15)
# 4.040 s (15219390 allocations: 5.10 GiB)
function fateman(n)
  f=(1+sum(Mvp.((:x,:y,:z,:t))))*one(n)
  f=f^n
  length(f*(f+1))
end

# @btime inv(Frac.([x+y x-y;x+1 y+1])) setup=(x=Mvp(:x);y=Mvp(:y))
# 195.739 μs (4483 allocations: 296.94 KiB)
end


text:
@testset "testing TICR, Dirichlet distribution" begin
# previously in PhyloNetworks

df = CSV.read(joinpath(dirname(Base.find_package("PhyloNetworks")),"..","examples","buckyCF.csv"));
d = readTableCF(df);

@testset "ticr! on data frame, tree" begin
net1 = readTopology("((((D:0.4,C:0.4):4.8,((A:0.8,B:0.8):2.2)#H1:2.2::0.7):4.0,(#H1:0::0.3,E:3.0):6.2):2.0,O:11.2);");
setGamma!(net1.edge[9],0.0);
net2 = deepcopy(net1);
# without optimizing branch lengthes
result1 = ticr!(net1,df,false);
@test result1[2] ≈ 1.480872194397731   # z statistic, from R: prop.test(2,15, p=0.05, alternative="greater", correct=F)
@test result1[1] ≈ 0.06932031690660927 # p-value, from R
@test result1[3] == Dict("[0.0, 0.01)"=>2, "[0.01, 0.05)"=>0, "[0.05, 0.1)"=>2, "[0.1, 1.0]"=>11)
@test result1[4][2] ≈ 48.152697007372566418 # pseudo log-lik obtained from R
@test result1[4][1] ≈ 10.576940922426542713 atol=1e-5 # alpha obtained from R
#with optimizing branch lengthes
result2 = ticr!(net2,d,true);
net2_1 = readTopology("(D,C,(((A,B):1.297)#H1:1.412::1.0,((#H1:0.0::0.0,E):6.877,O):1.354):9.352);");
result3 = ticr!(net2_1,d,false);
@test result3[2] ≈ 1.480872194397731   # z statistic, from R, same as above
@test result3[1] ≈ 0.06932031690660927 # p-value, from R
@test result3[3] == Dict("[0.0, 0.01)"=>2, "[0.01, 0.05)"=>0, "[0.05, 0.1)"=>2, "[0.1, 1.0]"=>11)
@test result3[4][2] ≈ 54.241562916699216 # pseudo log-lik obtained from R
@test result3[4][1] ≈ 20.128258663235194 # alpha obtained from R
result3_1 = ticr!(net2_1,d,false; test=:goodness);
@test result3_1[2] ≈ 25.962962962962965463 # chi-squared statistic obtained from R
@test result3_1[1] ≈ 9.7092282251534852702e-06 # p-value obtained from R
@test result3_1[3] == Dict("[0.0, 0.01)"=>2, "[0.01, 0.05)"=>0, "[0.05, 0.1)"=>2, "[0.1, 1.0]"=>11)
@test result3[4][2] ≈ 54.241562916699216 atol=1e-3 # pseudo log-lik obtained from R
@test result3[4][1] ≈ 20.128258663235194 atol=1e-4 # alpha obtained from R
end

@testset "ticr! on data frame, network, p value from maximum observed values, dirichlet distribution" begin
net3 = readTopology("((((D:0.4,C:0.4):4.8,((A:0.8,B:0.8):2.2)#H1:2.2::0.7):4.0,(#H1:0::0.3,E:3.0):6.2):2.0,O:11.2);");
# without optimizing branch lengths
netresult1 = ticr!(net3,d,false);
@test netresult1[2] ≈ -0.8885233166386386  # z stat, from R: prop.test(0,15, p=0.05, alternative="greater", correct=F)
@test netresult1[1] ≈ 0.8128703403598878   # p-value, from R
@test netresult1[3] == Dict("[0.0, 0.01)"=>0, "[0.01, 0.05)"=>0, "[0.05, 0.1)"=>2, "[0.1, 1.0]"=>13)
@test netresult1[4][2] ≈ 68.03708830981597 # pseudo log-lik
@test netresult1[4][1] ≈ 29.34808731515701 atol=1e-5 # alpha
end

@testset "ticr! on data frame, network, minimum pvalue, dirichlet distribution" begin
net3 = readTopology("((((D:0.4,C:0.4):4.8,((A:0.8,B:0.8):2.2)#H1:2.2::0.7):4.0,(#H1:0::0.3,E:3.0):6.2):2.0,O:11.2);");
# without optimizing branch lengths
netresult1 = ticr!(net3,d,false; quartetstat=:minpval);
@test netresult1[2] ≈ 8.589058727506838  # z stat
@test netresult1[1] ≈ 4.384234705965304e-18   # p-value
@test netresult1[3] == Dict("[0.05, 0.1)" => 1, "[0.0, 0.01)"  => 4, "[0.01, 0.05)" => 4, "[0.1, 1.0]"  => 6)
@test netresult1[4][2] ≈ 68.03708830981597 # pseudo log-lik
@test netresult1[4][1] ≈ 29.34808731515701 atol=1e-5 # alpha
end

end


text:
export TracyWidom

"""
Tracy-Widom distribution

The probability distribution of the normalized largest eigenvalue of a random
matrix with iid Gaussian matrix elements of variance 1/2 and mean 0.

The cdf of Tracy-Widom is given by

``
F_beta (s) = lim_{n→∞} Pr(√2 n^{1/6} (λ_max - √(2n) ≤ s)
``

where beta = 1, 2, or 4 for the orthogonal, unitary, or symplectic ensembles.

References:

1. doi:10.1016/0370-2693(93)91114-3
2. doi:10.1007/BF02100489
3. doi.org/10.1007/BF02099545

Numerical routines adapted from Alan Edelman's course notes for MIT 18.338,
Random Matrix Theory, 2016.
"""
struct TracyWidom <: ContinuousUnivariateDistribution end

"""
Cumulative density function of the Tracy-Widom distribution.

Computes the Tracy-Widom distribution by Bornemann's method of evaluating
a finite dimensional approximation to the Fredholm determinant using quadrature.

doi.org/10.1090/S0025-5718-09-02280-7

# Arguments
* `d::TracyWidom` or `Type{TracyWidom}`: an instance of `TracyWidom` or the type itself
* `s::Real`: The point at which to evaluate the cdf
* `beta::Integer = 2`: The Dyson index defining the distribution. Takes values 1, 2, or 4
* `num_points::Integer = 25`: The number of points in the quadrature
"""
function cdf(d::TracyWidom, s::T; beta::Integer=2, num_points::Integer=25) where {T<:Real}
    beta ∈ (1,2,4) || throw(ArgumentError("Beta must be 1, 2, or 4"))
    quad = gausslegendre(num_points)
    _TWcdf(s, beta, quad)
end

function cdf(d::Type{TracyWidom}, s::T; beta::Integer=2, num_points::Integer=25) where {T<: Real}
    cdf(d(), s, beta=beta, num_points=num_points)
end

function cdf(d::TracyWidom, s_vals::AbstractArray{T}; beta::Integer=2, num_points::Integer=25) where {T<:Real}
    beta ∈ (1,2,4) || throw(ArgumentError("Beta must be 1, 2, or 4"))
    quad = gausslegendre(num_points)
    [_TWcdf(s, beta, quad) for s in s_vals]
end

function cdf(d::Type{TracyWidom}, s_vals::AbstractArray{T}; beta::Integer=2, num_points::Integer=25) where {T<:Real}
    cdf(d(), s_vals, beta=beta, num_points=num_points)
end

function _TWcdf(s::T, beta::Integer, quad::Tuple{Array{Float64,1},Array{Float64,1}}) where {T<:Real}
    if beta == 2
        kernel = ((ξ,η) -> _K2tilde(ξ,η,s))
        return _fredholm_det(kernel, quad)
    elseif beta == 1
        kernel = ((ξ,η) -> _K1tilde(ξ,η,s))
        return _fredholm_det(kernel, quad)
    elseif beta == 4
        kernel2 = ((ξ,η) -> _K2tilde(ξ,η,s*sqrt(2)))
        kernel1 = ((ξ,η) -> _K1tilde(ξ,η,s*sqrt(2)))
        F2 = _fredholm_det(kernel2, quad)
        F1 = _fredholm_det(kernel1, quad)
        return (F1 + F2/F1) / 2
    end
end

function _fredholm_det(kernel::Function, quad::Tuple{Array{T,1},Array{T,1}}) where {T<:Real}
    nodes, weights = quad
    N = length(nodes)
    sqrt_weights = sqrt.(weights)
    weights_matrix = kron(transpose(sqrt_weights),sqrt_weights)
    K_matrix = [kernel(ξ,η) for ξ in nodes, η in nodes]
    det(I - weights_matrix .* K_matrix)
end

_ϕ(ξ, s) =  s + 10*tan(π*(ξ+1)/4)
_ϕprime(ξ) = (5π/2)*(sec(π*(ξ+1)/4))^2

# For beta = 2
function _airy_kernel(x, y)
    if x==y
        return (airyaiprime(x))^2 - x * (airyai(x))^2
    else
        return (airyai(x) * airyaiprime(y) - airyai(y) * airyaiprime(x)) / (x - y)
    end
end

_K2tilde(ξ,η,s) = sqrt(_ϕprime(ξ) * _ϕprime(η)) * _airy_kernel(_ϕ(ξ,s), _ϕ(η,s))

# For beta = 1
_A_kernel(x,y) = airyai((x+y)/2) / 2
_K1tilde(ξ,η,s) = sqrt(_ϕprime(ξ) * _ϕprime(η)) * _A_kernel(_ϕ(ξ,s), _ϕ(η,s))

"""
Samples the largest eigenvalue of the n × n GUE matrix
"""
function rand(d::TracyWidom, n::Int)
    n > 1 || error("Cannot construct $n × $n matrix")
    if n < 100
        k = n
    else #Exploit the fact that the eigenvector is concentrated in the top left corner
	k = round(Int, n-10*n^(1/3)-1)
    end
    a=randn(n-k+1)
    b=[χ(i) for i=(n-1):-1:k]
    v=eigmax(SymTridiagonal(a, b))
end
rand(d::Type{TracyWidom}, t::Integer) = rand(d(), t)


text:
using Test, JuLIP
using IPFitting
using IPFitting: Dat
using IPFitting.Data: configtype
using JuLIP: write_dict, read_dict
using JuLIP.Testing: print_tf
##

@info("Testing Dat:")
println("  * Dat constructor")
println("  * vec, devec")
println("  * energy(dat), forces(dat), virial(dat)")
println("  * (de-)dictionisation, ==")

at = rattle!(bulk(:Fe) * 3, 0.1)
lj = LennardJones(r0 = rnn(:Fe)) * C2Shift(2.5*rnn(:Fe))
dat = Dat(at, "test"; E = energy(lj, at),
                      F = forces(lj, at),
                      V = virial(lj, at) )

print_tf(@test(configtype(dat) == "test"))
print_tf(@test(energy(dat) == energy(lj, at)))
print_tf(@test(forces(dat) == forces(lj, at)))
print_tf(@test(virial(dat) ≈ virial(lj, at)))
println()

##

@info("Test (de-)serialisation of det")
print_tf(@test all(JuLIP.Testing.test_fio(dat)))
println()


text:
@testset "Writing" begin
    server = InfluxServer("http://localhost:8086")

    create_database(server, "write_test")
    dbs, = list_databases(server)
    @test "write_test" in dbs[!, :name]

    t = 2000
    measurements = [
        Measurement(
            "performance",
            Dict("cpu_load" => 0.1, "mem_free" => 123456, "stringint" => "1"),
            Dict("host" => "nureha"),
            t-1.5,
        ),
        Measurement(
            "performance",
            Dict("cpu_load" => 100.0, "mem_free" => 0, "stringint" => "1"),
            Dict("host" => "sadboi"),
            t-1.5,
        ),
        Measurement(
            "performance",
            Dict("cpu_load" => 0.2, "mem_free" => 234567, "stringint" => "1"),
            Dict("host" => "nureha"),
            t,
        ),
    ]

    write(server, "write_test", measurements)
    cpu_load, = query(server, "write_test", SELECT(;measurements=["performance"], condition=@querify(:host != "sadboi")))
    @test size(cpu_load) == (2,5)
    @test cpu_load[1, :cpu_load] == 0.1
    @test cpu_load[2, :cpu_load] == 0.2
    @test cpu_load[1, :mem_free] == 123456
    @test cpu_load[2, :mem_free] == 234567
    @test cpu_load[1, :host] == "nureha"
    @test cpu_load[2, :host] == "nureha"
    @test cpu_load[1, :stringint] == "1"
end