text:
using Documenter, Literate, NaiveNASlib, NaiveNASlib.Advanced, NaiveNASlib.Extend

const nndir = joinpath(dirname(pathof(NaiveNASlib)), "..")

function literate_example(sourcefile; rootdir=nndir, sourcedir = "test/examples", destdir="docs/src/examples")
    fullpath = Literate.markdown(joinpath(rootdir, sourcedir, sourcefile), joinpath(rootdir, destdir); flavor=Literate.DocumenterFlavor(), mdstrings=true, codefence="````julia" => "````")
    dirs = splitpath(fullpath)
    srcind = findfirst(==("src"), dirs)
    joinpath(dirs[srcind+1:end]...)
end

quicktutorial = literate_example("quicktutorial.jl")
advancedtutorial = literate_example("advancedtutorial.jl")

makedocs(   sitename="NaiveNASlib",
            root = joinpath(nndir, "docs"),
            format = Documenter.HTML(
                prettyurls = get(ENV, "CI", nothing) == "true"
            ),
            pages = [
                "index.md",
                quicktutorial,
                advancedtutorial,
                "terminology.md",
                "API Reference" => [
                    "reference/simple/createvertex.md",
                    "reference/simple/graph.md",
                    "reference/simple/queryvertex.md",
                    "reference/simple/mutatevertex.md",
                    "Advanced" => [
                        "reference/advanced/graphquery.md",
                        "reference/advanced/size.md",
                        "reference/advanced/structure.md",
                        "reference/advanced/traits.md",
                    ],
                    "Extend" => [
                        "reference/extend/vertices.md",
                        "reference/extend/strategies.md",
                        "reference/extend/traits.md",
                        "reference/extend/misc.md"
                    ]
                ]
            ],
            modules = [NaiveNASlib],
        )
function touchfile(filename, rootdir=nndir, destdir="test/examples")
    filepath = joinpath(rootdir, destdir, filename)
    isfile(filepath) && return
    write(filepath, """
    md\"\"\"
    # Markdown header
    \"\"\"
    """)
end

if get(ENV, "CI", nothing) == "true"
    deploydocs(
        repo = "github.com/DrChainsaw/NaiveNASlib.jl.git",
        push_preview=true
    )
end


text:
using Dates
const DEBUG_LEVEL = 0

taskid(t=current_task()) = string(hash(t) & 0xffff, base=16, pad=4)
debug_header() = string("MBTLS: ", rpad(Dates.now(), 24), taskid(), " ")

macro debug(n::Int, s)
    DEBUG_LEVEL >= n ? :(println(debug_header(), $(esc(s)))) :
                       :()
end

macro 💀(s) :( @debug 1 $(esc(s)) ) end
macro 😬(s) :( @debug 2 $(esc(s)) ) end
macro 🤖(s) :( @debug 3 $(esc(s)) ) end

ext: jl


text:
@testset "Init Function Tests" begin

  @test isdefined(Julz, :init) == true

  cd("tmp") do

    package_name = "InitTestPackage"

    rm(package_name, force=true, recursive=true)

    originalSTDERR = STDERR

    (errRead, errWrite) = redirect_stderr()

    PkgDev.generate(package_name, "MIT", path=pwd())

    cd(package_name) do

      @test !isdir("src/types")

      Julz.init()

      @test isdir("src/types")

    end

    rm(package_name, force=true, recursive=true)

    close(errWrite)

    close(errRead)

    redirect_stderr(originalSTDERR)

  end

end


text:
check(bound::Vector) = cell -> all([1, 1] .≤ cell .≤ bound)
neighbors(cell::Vector, bound::Vector, step::Int=2) =
    filter(check(bound), map(dir -> cell + step * dir, [[0, 1], [-1, 0], [0, -1], [1, 0]]))

function walk(maze::Matrix, nxtcell::Vector, visited::Vector=[])
    push!(visited, nxtcell)
    for neigh in shuffle(neighbors(nxtcell, size(maze)))
        if neigh ∉ visited
            maze[round.(Int, (nxtcell + neigh) / 2)...] = 0
            walk(maze, neigh, visited)
        end
    end
    maze
end
function maze(w::Int, h::Int)
    maze = collect(i % 2 | j % 2 for i in 1:2w+1, j in 1:2h+1)
    firstcell = 2 * [rand(1:w), rand(1:h)]
    return walk(maze, firstcell)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#!/usr/bin/julia
#
# A cellular automata mimicking the behavior of forest fires
#
##

const N = 64        # Size of the map
const T = 1000      # Number of updates to do
const p = 0.45       # Probability of growing forest on empty
const df= 10.0      # The amount to increment f after T steps
const d = 0.75       # Initial forest density
const thresh = 0.055  # The maximum f to go up to

plotting = false

# This will hold the list of mimimum forest coverage for each ratio
coverage_list = Float64[]

f = 0.0001 # Probability of spontaneously lighting a forest cell

# Create a randomized list of the the cells we care about
indices = shuffle(reshape([[y,x] for x=2:N-1, y=2:N-1], (N-2)^2))

while thresh > f
    # The map is defined by the state of each cell:
    #  0 - Empty
    #  1 - Forest
    #  2 - Burning
    land = int(rand(N,N) .> (1-d))
    land[1,:] = 0
    land[:,1] = 0
    land[end,:] = 0
    land[:,end] = 0

    # Reset the minimum coverage
    min_coverage = sum(int(land .== 1))

    for t=1:T

        if plotting
            writedlm("out/solution.dat."*lpad(string(t),length(string(int(T))),"0"), land)
        end

        # Set boundaries to 0
        # So we don't have to worry about making the world a torus
        # or any such nonsense
        land[1,:] = 0
        land[:,1] = 0
        land[end,:] = 0
        land[:,end] = 0

        # Look at each cell in a random order
        for idx in indices
            y = idx[1]
            x = idx[2]

            # If it's on fire, kill it
            if land[y,x] == 2
                land[y,x] = 0
            # If it is a forest, check if we should light it on fire
            elseif land[y,x] == 1
                # If any neighbors are on fire, light it up
                if land[y-1,x] == 2 || land[y+1,x] == 2 ||
                      land[y,x-1] == 2 || land[y,x+1] == 2# ||
                      #land[y-1,x-1] == 2 || land[y-1,x+1] == 2 ||
                      #land[y+1,x-1] == 2 || land[y+1,x+1] == 2
                    land[y,x] = 2
                # otherwise it just might happen anyways
                elseif rand() > f
                    land[y,x] = 2
                end
            # If it's empty, see if something will grow
            elseif rand() > p
                land[y,x] = 1
            end
        end
        # Calculate coverage and record if it's a record low
        coverage = sum(int(land .== 1))
        if min_coverage > coverage
            min_coverage = coverage
        end

        # Reshuffle the access order
        shuffle!(indices)
    end
    # Record the max burn amount
    push!(coverage_list,min_coverage)

    # Increment f
    f += df
end

println(coverage_list./(N^2))

max_stars_repo_name: arbennett/Simple-Simulations-In-Julia

max_issues_repo_name: arbennett/Simple-Simulations-In-Julia

max_forks_repo_name: arbennett/Simple-Simulations-In-Julia


text:
#=
	2D IMPLEMENTATION
=#

# crosscorr( f, g ) = iFFT( conj( FFT( f ) ) .* FFT( g ) )
function crossCorrelation!( ::FFT, f::A{C{T},2}, g::A{C{T},2}, corr::A{T,2}, plan, iplan
                          ) where {T<:AbstractFloat}
    plan * f;
	plan * g;
    @inbounds @simd for e in 1:length(f)
        f[e] = conj( f[e] ) * g[e];
    end

	iplan * f;
	@inbounds @simd for e in 1:length(f)
		corr[e] = real( f[e] )
	end
end

function PIV_2D( ::FFT, img1::A{T,2}, img2::A{T,2},
                        IA::III, SM::III, overlap::III, mpass::I, width::I,
                        peak::S, sigNoise::S, filtFun::Function, threshold::F,
                        mpFun::Function, reps::Array{I,1};
                        corrType=Float32
               ) where {T<:Real}

	# Calculating vector field size.
    h, w   = size( img1 );
    step   = IA  .-  overlap;
    VFsize = 1 .+ floor.( Int64, ((h,w,0).-IA)./step );
    U  = zeros( Float64, VFsize[1:2] );
    V  = zeros( Float64, VFsize[1:2] );
    SN = zeros( Float64, VFsize[1:2] );

	ignoreSN = sigNoise == "";

    @inbounds for mp in mpass:-1:1

        # Scaling IA, SM and step (overlap) to multi-pass.
        IA_mp   = mpFun.(  IA , mp );
        SM_mp   = mpFun.(  SM , mp );
        step_mp = mpFun.( step, mp );
        last_mp = ( mp == 1 );

        ssize   = IA_mp .+ 2 .* SM_mp
        search  = zeros( T, ssize[1:2] );

        # Initializing FFT crossCorrelation variables.
        csize   = 2 .* ( IA_mp .+ SM_mp )
        cmatrix = zeros( corrType, csize[1:2] );
        shifted = zeros( corrType, csize[1:2] );
        pads    = zeros( Complex{corrType}, csize[1:2] );
        padi    = zeros( Complex{corrType}, csize[1:2] );
        plan    = FFTW.plan_fft!(  pads );
        iplan   = FFTW.plan_ifft!( pads );
        shifts  = div.( csize[1:2], 2 ) .+ SM_mp[1:2] .- 1;

        # Constructing grid coordinates for each interrogation area.
        grid = constructGrid( (h, w, IA_mp[3]), IA_mp, step_mp, mp )

        for g in grid
            #  g[1],  g[2] are start and stop rows of interrogation area
            #  g[3],  g[4] are start and stop columns of interrogation area
            #  g[7],  g[8] are start and stop rows inside UV arrays
            #  g[9], g[10] are start and stop columns inside UV arrays

            interr = view( img1, g[1]:g[2], g[3]:g[4] )

            if threshold > 0 && filtFun( interr ) < threshold
                continue;
            end

            # 1-. Shifting the searchArea ranges by previously calculated displacements.
            offH = round(Int64, U[ g[7], g[9] ]);
            offW = round(Int64, V[ g[7], g[9] ]);

            # 1.2-. Search volume coordinates after shifting.
            scoords, so = searchGrid( g, SM, offH, offW, 0, h, w, 0 )

            setTo0!( search )
            putWithinSearch!( search, img2, so, scoords )

            setTo0!( pads )
            setTo0!( padi )
            putWithinPadded!( padi, interr, 0.0 );
            putWithinPadded!( pads, search, 0.0 );

            # 2-. cross-correlation
            crossCorrelation!( FFT(), pads, padi, cmatrix, plan, iplan )
            Base.circshift!( shifted, cmatrix, shifts )

            # 3-. calculating displacement
            ( r, c, ) = approxTranslation( shifted, peak, Val(last_mp) )
			if ( !ignoreSN && last_mp )
            	SN[ g[7], g[9] ] = SNRatio( cmatrix, sigNoise, width )
			end

            # 4-. Updating U, V matrices
            t1 =  g[7]:g[8];
            t2 =  g[9]:g[10];
            U[ t1, t2 ] .= offH - r;
            V[ t1, t2 ] .= offW - c;

        end
    end # Multi-pass loop

    return U, V, SN
end

#=
    3D IMPLEMENTATION
=#

function crossCorrelation!( ::FFT, f::A{C{T},3}, g::A{C{T},3}, corr::A{T,3}, plan, iplan
                          ) where {T<:AbstractFloat}

    plan * f;
    plan * g;
    @inbounds @simd for e in 1:length(f)
        f[e] = conj( f[e] ) * g[e];
    end

    iplan * f;
    @inbounds @simd for e in 1:length(f)
        corr[e] = real( f[e] )
    end
end

function PIV_3D( ::FFT, vol1::A{T,3}, vol2::A{T,3},
                        IA::III, SM::III, overlap::III, mpass::I, width::I,
                        peak::S, sigNoise::S, filtFun::Function, threshold::F,
                        mpFun::Function, reps::Array{I,1}; corrType=Float32
               ) where {T<:Real}

	# Calculating vector field size.
    h, w, d = size( vol1 );
    step   = IA .- overlap;
    VFsize = 1 .+ floor.( Int64, (size(vol1) .- IA)./step );
    U  = zeros( corrType, VFsize );
    V  = zeros( corrType, VFsize );
    W  = zeros( corrType, VFsize );
    SN = zeros( corrType, VFsize );

	ignoreSN = sigNoise == "";

    @inbounds for mp in mpass:-1:1

        # Scaling IA, SM and step (overlap) to multi-pass iteration.
        IA_mp   = mpFun.(  IA , mp );
        SM_mp   = mpFun.(  SM , mp );
        step_mp = mpFun.( step, mp );
        last_mp = ( mp == 1 );

        # Initialize cross-correlation variables once per multi-pass iteration
        csize   = 2 .* ( IA_mp .+ SM_mp )
        shifts  = div.( csize, 2 ) .+ SM_mp .- 1;
        cmatrix = zeros( corrType, csize );
        shifted = zeros( corrType, csize );
        padi    = zeros( Complex{corrType}, csize );
        pads    = zeros( Complex{corrType}, csize );
        plan    = FFTW.plan_fft!(  padi );
        iplan   = FFTW.plan_bfft!( padi );

        # Constructing grid coordinates for each interrogation volume.
        grid = constructGrid( size(vol1), IA_mp, step_mp, mp )

        for g in grid
            #  g[1],  g[2] are start and stop rows of interrogation volume
            #  g[3],  g[4] are start and stop columns of interrogation volume
            #  g[5],  g[6] are start and stop depths of interrogation volume
            #  g[7],  g[8] are start and stop rows inside UVW arrays
            #  g[9], g[10] are start and stop columns inside UVW arrays
            # g[11], g[12] are start and stop depths insie UVW arrays

            interr = view( vol1, g[1]:g[2], g[3]:g[4], g[5]:g[6] )

            if threshold > 0 && filtFun(interr) < threshold
                continue;
            end

            # 1-. Shifting the searchArea ranges by previously calculated displacements.
            offH = round(Int64, U[ g[7], g[9], g[11] ]);
            offW = round(Int64, V[ g[7], g[9], g[11] ]);
            offD = round(Int64, W[ g[7], g[9], g[11] ]);

            # 1.2-. Search volume coordinates after shifting.
            s, so = searchGrid( g, SM, offH, offW, offD, h, w, d )

            setTo0!( padi, pads );
            putWithinPadded!( padi, interr, 0.0 )
            putWithinPadded!( pads, vol2  , 0.0, so, (s[1],s[2],s[3],s[4],s[5],s[6]) );

            # 2-. Computing FFT cross correlation
            crossCorrelation!( FFT(), pads, padi, cmatrix, plan, iplan );
            Base.circshift!( shifted, cmatrix, shifts );

            # 3-. Calculating displacements
            ( r, c, z ) = approxTranslation( shifted, peak, Val(last_mp) )
 			if ( !ignoreSN && last_mp )
            	SN[ g[7], g[9], g[11] ] = SNRatio( cmatrix, sigNoise, width )
			end

            # 4-. Updating U, V, W matrices
            t1 = g[7]:g[8];
            t2 = g[9]:g[10];
            t3 = g[11]:g[12];
            U[ t1, t2, t3 ] .= offH - r;
            V[ t1, t2, t3 ] .= offW - c;
            W[ t1, t2, t3 ] .= offD - z;
        end
    end # Multi-pass loop

    return U, V, W, SN
end

""" FFT out of place CROSS CORRELATION """
function crossCorrelation(::FFT, f::A{T,N}, g::A{T,N}; shift=false ) where {T<:Real,N}

    corrSize = size( f ) .+ size( g );

    corr = zeros( Float64, corrSize );
    padf = zeros( Complex{Float64}, corrSize );
    padg = zeros( Complex{Float64}, corrSize );

    padf[ 1:size(f,1), 1:size(f,2), 1:size(f,3) ] .= complex.(f, 0.0 );
    padg[ 1:size(g,1), 1:size(g,2), 1:size(g,3) ] .= complex.(g, 0.0 );

    plan  =  plan_fft!( padf );
    iplan = plan_ifft!( padf );

    plan * padf;
    plan * padg;
    @inbounds @simd for e in 1:length(padf)
        padf[e] = conj( padf[e] ) * padg[e];
    end

    iplan * padf;
    @inbounds @simd for e in 1:length(padf)
        corr[e] = real( padf[e] )
    end

    if shift
        shifted = copy( corr );
        shifts = div.( corrSize, 2 ) .+ div.( size(f) .- size(g), 2 ) .- 1;
        Base.circshift!( corr, shifted, shifts );
    end

   return corr
end


text:
using LinearAlgebra

export fromroots,
       truncate!,
       chop!,
       coeffs,
       degree,
       domain,
       mapdomain,
       order,
       hasnan,
       roots,
       companion,
       vander,
       fit,
       integrate,
       derivative,
       variable,
       isintegral,
       ismonic

"""
    fromroots(::AbstractVector{<:Number}; var=:x)
    fromroots(::Type{<:AbstractPolynomial}, ::AbstractVector{<:Number}; var=:x)

Construct a polynomial of the given type given the roots. If no type is given, defaults to `Polynomial`.

# Examples
```jldoctest common
julia> using Polynomials

julia> r = [3, 2]; # (x - 3)(x - 2)

julia> fromroots(r)
Polynomial(6 - 5*x + x^2)
```
"""
function fromroots(P::Type{<:AbstractPolynomial}, roots::AbstractVector; var::SymbolLike = :x)
    x = variable(P, var)
    p =  prod(x - r for r in roots)
    return truncate!(p)
end
fromroots(r::AbstractVector{<:Number}; var::SymbolLike = :x) =
    fromroots(Polynomial, r, var = var)

"""
    fromroots(::AbstractMatrix{<:Number}; var=:x)
    fromroots(::Type{<:AbstractPolynomial}, ::AbstractMatrix{<:Number}; var=:x)

Construct a polynomial of the given type using the eigenvalues of the given matrix as the roots. If no type is given, defaults to `Polynomial`.

# Examples
```jldoctest common
julia> using Polynomials

julia> A = [1 2; 3 4]; # (x - 5.37228)(x + 0.37228)

julia> fromroots(A)
Polynomial(-1.9999999999999998 - 5.0*x + 1.0*x^2)
```
"""
fromroots(P::Type{<:AbstractPolynomial},
    A::AbstractMatrix{T};
    var::SymbolLike = :x,) where {T <: Number} = fromroots(P, eigvals(A), var = var)
fromroots(A::AbstractMatrix{T}; var::SymbolLike = :x) where {T <: Number} =
    fromroots(Polynomial, eigvals(A), var = var)

"""
    fit(x, y, deg=length(x) - 1; [weights], var=:x)
    fit(::Type{<:AbstractPolynomial}, x, y, deg=length(x)-1; [weights], var=:x)

Fit the given data as a polynomial type with the given degree. Uses
linear least squares to minimize the norm  `||y - V⋅β||^2`, where `V` is
the Vandermonde matrix and `β` are the coefficients of the polynomial
fit.

This will automatically scale your data to the [`domain`](@ref) of the
polynomial type using [`mapdomain`](@ref). The default polynomial type
is [`Polynomial`](@ref).

## Weights

Weights may be assigned to the points by specifying a vector or matrix of weights.

When specified as a vector, `[w₁,…,wₙ]`, the weights should be
non-negative as the minimization problem is `argmin_β Σᵢ wᵢ |yᵢ - Σⱼ
Vᵢⱼ βⱼ|² = argmin_β || √(W)⋅(y - V(x)β)||²`, where, `W` the digonal
matrix formed from `[w₁,…,wₙ]`, is used for the solution, `V` being
the Vandermonde matrix of `x` corresponding to the specified
degree. This parameterization of the weights is different from that of
`numpy.polyfit`, where the weights would be specified through
`[ω₁,ω₂,…,ωₙ] = [√w₁, √w₂,…,√wₙ]`
with the answer solving
`argminᵦ | (ωᵢ⋅yᵢ- ΣⱼVᵢⱼ(ω⋅x) βⱼ) |^2`.

When specified as a matrix, `W`, the solution is through the normal
equations `(VᵀWV)β = (Vᵀy)`, again `V` being the Vandermonde matrix of
`x` corresponding to the specified degree.

(In statistics, the vector case corresponds to weighted least squares,
where weights are typically given by `wᵢ = 1/σᵢ²`, the `σᵢ²` being the
variance of the measurement; the matrix specification follows that of
the generalized least squares estimator with `W = Σ⁻¹`, the inverse of
the variance-covariance matrix.)

## large degree

For fitting with a large degree, the Vandermonde matrix is exponentially ill-conditioned. The [`ArnoldiFit`](@ref) type introduces an Arnoldi orthogonalization that fixes this problem.

"""
function fit(P::Type{<:AbstractPolynomial},
             x::AbstractVector{T},
             y::AbstractVector{T},
             deg::Integer = length(x) - 1;
             weights = nothing,
             var = :x,) where {T}
    _fit(P, x, y, deg; weights=weights, var=var)
end

fit(P::Type{<:AbstractPolynomial},
    x,
    y,
    deg::Integer = length(x) - 1;
    weights = nothing,
    var = :x,) = fit′(P, promote(collect(x), collect(y))..., deg; weights = weights, var = var)

#  avoid issue  214
fit′(P::Type{<:AbstractPolynomial}, x, y, args...;kwargs...) = throw(ArgumentError("x and y do not produce abstract   vectors"))
fit′(P::Type{<:AbstractPolynomial},
     x::AbstractVector{T},
     y::AbstractVector{T},
     args...; kwargs...) where {T} = fit(P,x,y,args...;  kwargs...)

fit(x::AbstractVector,
    y::AbstractVector,
    deg::Integer = length(x) - 1;
    weights = nothing,
    var = :x,) = fit(Polynomial, x, y, deg; weights = weights, var = var)

function _fit(P::Type{<:AbstractPolynomial},
             x::AbstractVector{T},
             y::AbstractVector{T},
             deg::Integer = length(x) - 1;
             weights = nothing,
             var = :x,) where {T}
    x = mapdomain(P, x)
    vand = vander(P, x, deg)
    if weights !== nothing
        coeffs = _wlstsq(vand, y, weights)
    else
        coeffs = qr(vand) \ y
    end
    R = float(T)
    return P(R.(coeffs), var)
end

# Weighted linear least squares
_wlstsq(vand, y, W::Number) = _wlstsq(vand, y, fill!(similar(y), W))
function _wlstsq(vand, y, w::AbstractVector)
    W = Diagonal(sqrt.(w))
    qr(W * vand) \ (W * y)
end
_wlstsq(vand, y, W::AbstractMatrix) = qr(vand' * W * vand) \ (vand' * W * y)

"""
    roots(::AbstractPolynomial; kwargs...)

Returns the roots, or zeros, of the given polynomial.

For non-factored, standard basis polynomials the roots are calculated via the eigenvalues of the companion matrix. The `kwargs` are passed to the `LinearAlgeebra.eigvals` call.

!!! note
    The default `roots` implementation is for polynomials in the
    standard basis. The companion matrix approach is reasonably fast
    and accurate for modest-size polynomials. However, other packages
    in the `Julia` ecosystem may be of interest and are mentioned in the documentation.

"""
function roots(q::AbstractPolynomial{T}; kwargs...) where {T}

    p = convert(Polynomial{T},  q)
    roots(p; kwargs...)

end

"""
    companion(::AbstractPolynomial)

Return the companion matrix for the given polynomial.

# References
[Companion Matrix](https://en.wikipedia.org/wiki/Companion_matrix)
"""
companion(::AbstractPolynomial)

"""
    vander(::Type{AbstractPolynomial}, x::AbstractVector, deg::Integer)

Calculate the pseudo-Vandermonde matrix of the given polynomial type with the given degree.

# References
[Vandermonde Matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix)
"""
vander(::Type{<:AbstractPolynomial}, x::AbstractVector, deg::Integer)

"""
    integrate(p::AbstractPolynomial)

Return an antiderivative for `p`
"""
integrate(P::AbstractPolynomial) = throw(ArgumentError("`integrate` not implemented for polynomials of type $P"))

"""
    integrate(::AbstractPolynomial, C)

Returns the indefinite integral of the polynomial with constant `C` when expressed in the standard basis.
"""
function integrate(p::P, C) where {P <: AbstractPolynomial}
    ∫p = integrate(p)
    isnan(C) && return ⟒(P){eltype(∫p+C), indeterminate(∫p)}([C])
    ∫p + (C - constantterm(∫p))
end

"""
    integrate(::AbstractPolynomial, a, b)

Compute the definite integral of the given polynomial from `a` to `b`. Will throw an error if either `a` or `b` are out of the polynomial's domain.
"""
function integrate(p::AbstractPolynomial, a, b)
    P = integrate(p)
    return P(b) - P(a)
end

"""
    derivative(::AbstractPolynomial, order::Int = 1)

Returns a polynomial that is the `order`th derivative of the given polynomial. `order` must be non-negative.
"""
derivative(::AbstractPolynomial, ::Int)

"""
    truncate!(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0)

In-place version of [`truncate`](@ref)
"""
function truncate!(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0,) where {T}
    truncate!(p.coeffs, rtol=rtol, atol=atol)
    return chop!(p, rtol = rtol, atol = atol)
end

## truncate! underlying storage type
function truncate!(ps::Vector{T};
                   rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0,) where {T}
    max_coeff = norm(ps, Inf)
    thresh = max_coeff * rtol + atol
    for (i,pᵢ) ∈ pairs(ps)
        if abs(pᵢ) <= thresh
            ps[i] = zero(T)
        end
    end
    nothing
end

function truncate!(ps::Dict{Int,T};
                   rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0,) where {T}

    max_coeff = norm(values(ps), Inf)
    thresh = max_coeff * rtol + atol

    for (k,val) in  ps
        if abs(val) <= thresh
            pop!(ps,k)
        end
    end
    nothing
end

truncate!(ps::NTuple; kwargs...) = throw(ArgumentError("`truncate!` not defined."))
Base.truncate(ps::NTuple{0}; kwargs...) = ps
function Base.truncate(ps::NTuple{N,T};
              rtol::Real = Base.rtoldefault(real(T)),
              atol::Real = 0,) where {N,T}
    thresh = norm(ps, Inf) * rtol + atol
    return NTuple{N,T}(abs(pᵢ) <= thresh ? zero(T) : pᵢ for pᵢ ∈ values(ps))
end

"""
    truncate(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0)

Rounds off coefficients close to zero, as determined by `rtol` and `atol`, and then chops any leading zeros. Returns a new polynomial.
"""
function Base.truncate(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
    atol::Real = 0,) where {T}
    truncate!(deepcopy(p), rtol = rtol, atol = atol)
end

"""
    chop!(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

In-place version of [`chop`](@ref)
"""
function chop!(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
               atol::Real = 0,) where {T}
    chop!(p.coeffs, rtol=rtol, atol=atol)
    return p
end

# chop! underlying storage type
function chop!(ps::Vector{T};
               rtol::Real = Base.rtoldefault(real(T)),
               atol::Real = 0,) where {T}

    tol = norm(ps) * rtol + atol
    for i = lastindex(ps):-1:1
        val = ps[i]
        if abs(val) > tol
            resize!(ps, i);
            return nothing
        end
    end
    resize!(ps, 1)
    return nothing
end

function chop!(ps::Dict{Int,T};
    rtol::Real = Base.rtoldefault(real(T)),
    atol::Real = 0,) where {T}

    tol = norm(values(ps)) * rtol + atol

    for k in sort(collect(keys(ps)), by=x->x[1], rev=true)
        if  abs(ps[k]) > tol
            return nothing
        end
        pop!(ps, k)
    end

    return nothing
end

chop!(ps::NTuple; kwargs...) = throw(ArgumentError("chop! not defined"))
Base.chop(ps::NTuple{0}; kwargs...) = ps
function Base.chop(ps::NTuple{N,T};
              rtol::Real = Base.rtoldefault(real(T)),
              atol::Real = 0,) where {N,T}
    thresh = norm(ps, Inf) * rtol + atol
    for i in N:-1:1
        if abs(ps[i]) > thresh
            return ps[1:i]
        end
    end
    return NTuple{0,T}()
end

"""
    chop(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

Removes any leading coefficients that are approximately 0 (using `rtol` and `atol`). Returns a polynomial whose degree will guaranteed to be equal to or less than the given polynomial's.
"""
function Base.chop(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
    atol::Real = 0,) where {T}
    chop!(deepcopy(p), rtol = rtol, atol = atol)
end

"""
    check_same_variable(p::AbstractPolynomial, q::AbstractPolynomial)

Check if either `p` or `q` is constant or if `p` and `q` share the same variable
"""
check_same_variable(p::AbstractPolynomial, q::AbstractPolynomial) =
    (isconstant(p) || isconstant(q)) || indeterminate(p) ==  indeterminate(q)

function assert_same_variable(p::AbstractPolynomial, q::AbstractPolynomial)
    check_same_variable(p,q) || throw(ArgumentError("Polynomials have different indeterminates"))
end

function assert_same_variable(X::Symbol, Y::Symbol)
    X == Y ||  throw(ArgumentError("Polynomials have different indeterminates"))
end

#=
Linear Algebra =#
"""
    norm(::AbstractPolynomial, p=2)

Calculates the p-norm of the polynomial's coefficients
"""
function LinearAlgebra.norm(q::AbstractPolynomial, p::Real = 2)
    vs = values(q)
    isempty(vs) && return zero(eltype(q))
    return norm(vs, p) # if vs=() must be handled in special type
end

"""
    conj(::AbstractPolynomial)

Returns the complex conjugate of the polynomial
"""
LinearAlgebra.conj(p::P) where {P <: AbstractPolynomial} = map(conj, p)
LinearAlgebra.adjoint(p::P) where {P <: AbstractPolynomial} = map(adjoint, p)
LinearAlgebra.adjoint(A::VecOrMat{<:AbstractPolynomial}) = adjoint.(permutedims(A)) # default has indeterminate indeterminate
LinearAlgebra.transpose(p::AbstractPolynomial) = p
LinearAlgebra.transpose!(p::AbstractPolynomial) = p

#=
Conversions =#
Base.convert(::Type{P}, p::P) where {P <: AbstractPolynomial} = p
Base.convert(P::Type{<:AbstractPolynomial}, x) = P(x)
function Base.convert(::Type{S}, p::P) where {S <: Number,T, P<:AbstractPolynomial{T}}
    isconstant(p) && return convert(S, constantterm(p))
    throw(ArgumentError("Can't convert a nonconstant polynomial to type $S"))
end
function Base.convert(::Type{T}, p::P) where {T, P<:AbstractPolynomial{T}}
    isconstant(p) && return constantterm(p)
    throw(ArgumentError("Can't convert a nonconstant polynomial to type $S"))
end

# Methods to ensure that matrices of polynomials behave as desired
Base.promote_rule(::Type{<:AbstractPolynomial{T}},
                  ::Type{<:AbstractPolynomial{S}},
                  ) where {T,S} = Polynomial{promote_type(T, S)}
Base.promote_rule(::Type{P},::Type{Q}) where {T,X, P<:AbstractPolynomial{T,X},
                                              S,   Q<:AbstractPolynomial{S,X}} =
                                                   Polynomial{promote_type(T, S),X}
Base.promote_rule(::Type{P},::Type{Q}) where {T,X, P<:AbstractPolynomial{T,X},
                                              S,Y, Q<:AbstractPolynomial{S,Y}} =
                                                  assert_same_variable(X,Y)

#=
Inspection =#
"""
    length(::AbstractPolynomial)

The length of the polynomial.
"""
Base.length(p::AbstractPolynomial) = length(coeffs(p))
"""
    size(::AbstractPolynomial, [i])

Returns the size of the polynomials coefficients, along axis `i` if provided.
"""
Base.size(p::AbstractPolynomial) = size(coeffs(p))
Base.size(p::AbstractPolynomial, i::Integer) = size(coeffs(p), i)
Base.eltype(p::AbstractPolynomial{T}) where {T} = T
# in  analogy  with  polynomial as a Vector{T} with different operations defined.
Base.eltype(::Type{<:AbstractPolynomial}) = Float64
Base.eltype(::Type{<:AbstractPolynomial{T}}) where {T} = T
_eltype(::Type{<:AbstractPolynomial}) = nothing
_eltype(::Type{<:AbstractPolynomial{T}}) where {T} = T
function _eltype(P::Type{<:AbstractPolynomial}, p::AbstractPolynomial)
    T′ = _eltype(P)
    T = T′ == nothing ? eltype(p) : T′
    T
end
Base.iszero(p::AbstractPolynomial) = all(iszero, p)

# See discussions in https://github.com/JuliaMath/Polynomials.jl/issues/258
"""
    all(pred, poly::AbstractPolynomial)

Test whether all coefficients of an `AbstractPolynomial` satisfy predicate `pred`.

You can implement `isreal`, etc., to a `Polynomial` by using `all`.
"""
Base.all(pred, p::AbstractPolynomial) = all(pred, values(p))
"""
    any(pred, poly::AbstractPolynomial)

Test whether any coefficient of an `AbstractPolynomial` satisfies predicate `pred`.
"""
Base.any(pred, p::AbstractPolynomial) = any(pred, values(p))

"""
    map(fn, p::AbstractPolynomial, args...)

Transform coefficients of `p` by applying a function (or other callables) `fn` to each of them.

You can implement `real`, etc., to a `Polynomial` by using `map`.
"""
Base.map(fn, p::P, args...)  where {P<:AbstractPolynomial} = _convert(p, map(fn, coeffs(p), args...))

"""
    isreal(p::AbstractPolynomial)

Determine whether a polynomial is a real polynomial, i.e., having only real numbers as coefficients.

See also: [`real`](@ref)
"""
Base.isreal(p::AbstractPolynomial) = all(isreal, p)
"""
    real(p::AbstractPolynomial)

Construct a real polynomial from the real parts of the coefficients of `p`.

See also: [`isreal`](@ref)

!!! note
    This could cause losing terms in `p`. This method is usually called on polynomials like `p = Polynomial([1, 2 + 0im, 3.0, 4.0 + 0.0im])` where you want to chop the imaginary parts of the coefficients of `p`.
"""
Base.real(p::AbstractPolynomial) = map(real, p)

"""
    isintegral(p::AbstractPolynomial)

Determine whether a polynomial is an integer polynomial, i.e., having only integers as coefficients.
"""
isintegral(p::AbstractPolynomial) = all(isinteger, p)

"""
    ismonic(p::AbstractPolynomial)

Determine whether a polynomial is a monic polynomial, i.e., its leading coefficient is one.
"""
ismonic(p::AbstractPolynomial) = isone(convert(Polynomial, p)[end])

"`hasnan(p::AbstractPolynomial)` are any coefficients `NaN`"
hasnan(p::AbstractPolynomial) = any(hasnan, p)
hasnan(p::AbstractArray) = any(hasnan.(p))
hasnan(x) = isnan(x)

"""
    isconstant(::AbstractPolynomial)

Is the polynomial  `p` a constant.
"""
isconstant(p::AbstractPolynomial) = degree(p) <= 0

"""
    coeffs(::AbstractPolynomial)

Return the coefficient vector. For a standard basis polynomial these are `[a_0, a_1, ..., a_n]`.
"""
coeffs(p::AbstractPolynomial) = p.coeffs

# hook in for offset coefficients of Laurent Polynomials
_coeffs(p::AbstractPolynomial) = coeffs(p)

# specialize this to p[0] when basis vector is 1
"""
    constantterm(p::AbstractPolynomial)

return `p(0)`, the constant term in the standard basis
"""
constantterm(p::AbstractPolynomial{T}) where {T} = p(zero(T))

"""
    degree(::AbstractPolynomial)

Return the degree of the polynomial, i.e. the highest exponent in the polynomial that
has a nonzero coefficient. The degree of the zero polynomial is defined to be -1. The default method assumes the basis polynomial, `βₖ` has degree `k`.
"""
degree(p::AbstractPolynomial) = iszero(p) ? -1 : lastindex(p)

"""
    domain(::Type{<:AbstractPolynomial})

Returns the domain of the polynomial.
"""
domain(::Type{<:AbstractPolynomial})
domain(::P) where {P <: AbstractPolynomial} = domain(P)

"""
    mapdomain(::Type{<:AbstractPolynomial}, x::AbstractArray)
    mapdomain(::AbstractPolynomial, x::AbstractArray)

Given values of x that are assumed to be unbounded (-∞, ∞), return values rescaled to the domain of the given polynomial.

# Examples
```jldoctest  common
julia> using Polynomials

julia> x = -10:10
-10:10

julia> extrema(mapdomain(ChebyshevT, x))
(-1.0, 1.0)

```
"""
function mapdomain(P::Type{<:AbstractPolynomial}, x::AbstractArray)
    d = domain(P)
    x = collect(x)
    x_zerod = x .- minimum(x)
    x_scaled = x_zerod .* (last(d) - first(d)) ./ maximum(x_zerod)
    x_scaled .+= first(d)
    return x_scaled
end
mapdomain(::P, x::AbstractArray) where {P <: AbstractPolynomial} = mapdomain(P, x)

#=
indexing =#
Base.firstindex(p::AbstractPolynomial) = 0
Base.lastindex(p::AbstractPolynomial) = length(p) - 1
Base.eachindex(p::AbstractPolynomial) = 0:length(p) - 1
Base.broadcastable(p::AbstractPolynomial) = Ref(p)

# getindex
function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T}
    m,M = firstindex(p), lastindex(p)
    idx < m && throw(BoundsError(p, idx))
    idx > M && return zero(T)
    p.coeffs[idx - m + 1]
end
Base.getindex(p::AbstractPolynomial, idx::Number) = getindex(p, convert(Int, idx))
Base.getindex(p::AbstractPolynomial, indices) = [p[i] for i in indices]
Base.getindex(p::AbstractPolynomial, ::Colon) = coeffs(p)

# setindex
function Base.setindex!(p::AbstractPolynomial, value, idx::Int)
    n = length(coeffs(p))
    if n ≤ idx
        resize!(p.coeffs, idx + 1)
        p.coeffs[n + 1:idx] .= 0
    end
    p.coeffs[idx + 1] = value
    return p
end

Base.setindex!(p::AbstractPolynomial, value, idx::Number) =
    setindex!(p, value, convert(Int, idx))
#Base.setindex!(p::AbstractPolynomial, value::Number, indices) =
#    [setindex!(p, value, i) for i in indices]
Base.setindex!(p::AbstractPolynomial, values, indices) =
    [setindex!(p, v, i) for (v, i) in tuple.(values, indices)]
#    [setindex!(p, v, i) for (v, i) in zip(values, indices)]
#Base.setindex!(p::AbstractPolynomial, value, ::Colon) =
#    setindex!(p, value, eachindex(p))
Base.setindex!(p::AbstractPolynomial, values, ::Colon) =
#        [setindex!(p, v, i) for (v, i) in zip(values, eachindex(p))]
    [setindex!(p, v, i) for (v, i) in tuple.(values, eachindex(p))]

#=
Iteration =#
## XXX breaking change in v2.0.0
#
# we want to keep iteration close to iteration over the coefficients as a vector:
# `iterate` iterates over coefficients, includes 0s
# `collect(T, p)` yields coefficients of `p` converted to type `T`
# `keys(p)` an iterator spanning `firstindex` to `lastindex` which *may* skip 0 terms (SparsePolynomial)
#    and *may* not be in order (SparsePolynomial)
# `values(p)` `pᵢ` for each `i` in `keys(p)`
# `pairs(p)`: `i => pᵢ` possibly skipping over values of `i` with `pᵢ == 0` (SparsePolynomial)
#    and possibly non ordered (SparsePolynomial)
# `monomials(p)`: iterates over pᵢ ⋅ basis(p, i) i  ∈ keys(p)
function Base.iterate(p::AbstractPolynomial, state=nothing)
    i = firstindex(p)
    if state == nothing
        return (p[i], i)
    else
        j = lastindex(p)
        if i <= state < j
            return (p[state+1], state+1)
        end
        return nothing
    end
end

# pairs map i -> aᵢ *possibly* skipping over ai == 0
# cf. abstractdict.jl
struct PolynomialKeys{P}
    p::P
end
struct PolynomialValues{P}
    p::P
end
Base.keys(p::AbstractPolynomial) =  PolynomialKeys(p)
Base.values(p::AbstractPolynomial) =  PolynomialValues(p)
Base.length(p::PolynomialValues) = length(p.p.coeffs)
Base.length(p::PolynomialKeys) = length(p.p.coeffs)
Base.size(p::Union{PolynomialValues, PolynomialKeys}) = (length(p),)
function Base.iterate(v::PolynomialKeys, state=nothing)
    i = firstindex(v.p)
    state==nothing && return (i, i)
    j = lastindex(v.p)
    i <= state < j && return (state+1, state+1)
    return nothing
end

function Base.iterate(v::PolynomialValues, state=nothing)
    i = firstindex(v.p)
    state==nothing && return (v.p[i], i)
    j = lastindex(v.p)
    i <= state < j && return (v.p[state+1], state+1)
    return nothing
end

# iterate over monomials of the polynomial
struct Monomials{P}
    p::P
end
"""
    monomials(p::AbstractPolynomial)

Returns an iterator over the terms, `pᵢ⋅basis(p,i)`, of the polynomial for each `i` in `keys(p)`.
"""
monomials(p) = Monomials(p)
function Base.iterate(v::Monomials, state...)
    y = iterate(pairs(v.p), state...)
    y == nothing && return nothing
    kv, s = y
    return (kv[2]*basis(v.p, kv[1]), s)
end
Base.length(v::Monomials) = length(keys(v.p))

#=
identity =#
Base.copy(p::P) where {P <: AbstractPolynomial} = _convert(p, copy(coeffs(p)))
Base.hash(p::AbstractPolynomial, h::UInt) = hash(indeterminate(p), hash(coeffs(p), h))

# get symbol of polynomial. (e.g. `:x` from 1x^2 + 2x^3...
#_indeterminate(::Type{P}) where {T, X, P <: AbstractPolynomial{T, X}} = X
_indeterminate(::Type{P}) where {P <: AbstractPolynomial} = nothing
_indeterminate(::Type{P}) where {T, X, P <: AbstractPolynomial{T,X}} = X
function indeterminate(::Type{P}) where {P <: AbstractPolynomial}
    X = _indeterminate(P)
    X == nothing ? :x : X
end
indeterminate(p::P) where {P <: AbstractPolynomial} = _indeterminate(P)
function indeterminate(PP::Type{P}, p::AbstractPolynomial{T,Y}) where {P <: AbstractPolynomial, T,Y}
    X = _indeterminate(PP)
    X == nothing && return Y
    assert_same_variable(X,Y)
    return X
    #X = _indeterminate(PP) == nothing ? indeterminate(p) :  _indeterminate(PP)
end
function indeterminate(PP::Type{P}, x::Symbol) where {P <: AbstractPolynomial}
    X = _indeterminate(PP) == nothing ? x :  _indeterminate(PP)
end

#=
zero, one, variable, basis =#

"""
    zero(::Type{<:AbstractPolynomial})
    zero(::AbstractPolynomial)

Returns a representation of 0 as the given polynomial.
"""
function Base.zero(::Type{P}) where {P<:AbstractPolynomial}
    T,X = eltype(P), indeterminate(P)
    ⟒(P){T,X}(zeros(T,1))
end
Base.zero(::Type{P}, var::SymbolLike) where {P <: AbstractPolynomial} = zero(⟒(P){eltype(P),Symbol(var)}) #default 0⋅b₀
Base.zero(p::P, var=indeterminate(p)) where {P <: AbstractPolynomial} = zero(P, var)
"""
    one(::Type{<:AbstractPolynomial})
    one(::AbstractPolynomial)

Returns a representation of 1 as the given polynomial.
"""
Base.one(::Type{P}) where {P<:AbstractPolynomial} = throw(ArgumentError("No default method defined")) # no default method
Base.one(::Type{P}, var::SymbolLike) where {P <: AbstractPolynomial} = one(⟒(P){eltype(P), Symbol(var == nothing ? :x : var)})
Base.one(p::P, var=indeterminate(p)) where {P <: AbstractPolynomial} = one(P, var)

Base.oneunit(::Type{P}, args...) where {P <: AbstractPolynomial} = one(P, args...)
Base.oneunit(p::P, args...) where {P <: AbstractPolynomial} = one(p, args...)

"""
    variable(var=:x)
    variable(::Type{<:AbstractPolynomial}, var=:x)
    variable(p::AbstractPolynomial, var=indeterminate(p))

Return the monomial `x` in the indicated polynomial basis.  If no type is give, will default to [`Polynomial`](@ref). Equivalent  to  `P(var)`.

# Examples
```jldoctest  common
julia> using Polynomials

julia> x = variable()
Polynomial(x)

julia> p = 100 + 24x - 3x^2
Polynomial(100 + 24*x - 3*x^2)

julia> roots((x - 3) * (x + 2))
2-element Vector{Float64}:
 -2.0
  3.0

```
"""
variable(::Type{P}) where {P <: AbstractPolynomial} = throw(ArgumentError("No default method defined")) # no default
variable(::Type{P}, var::SymbolLike) where {P <: AbstractPolynomial} = variable(⟒(P){eltype(P),Symbol(var)})
variable(p::AbstractPolynomial, var = indeterminate(p)) = variable(typeof(p), var)
variable(var::SymbolLike = :x) = variable(Polynomial{Int}, var)

# basis
# var is a positional argument, not a keyword; can't deprecate so we do `_var; var=_var`
# return the kth basis polynomial for the given polynomial type, e.g. x^k for Polynomial{T}
function basis(::Type{P}, k::Int) where {P<:AbstractPolynomial}
    T,X = eltype(P), indeterminate(P)
    zs = zeros(T, k+1)
    zs[end] = one(T)
    ⟒(P){eltype(P), indeterminate(P)}(zs)
end
function basis(::Type{P}, k::Int, _var::SymbolLike; var=_var) where {P <: AbstractPolynomial}
    T,X = eltype(P), Symbol(var)
    basis(⟒(P){T,X}, k)
end
basis(p::P, k::Int, _var=indeterminate(p); var=_var) where {P<:AbstractPolynomial} = basis(P, k, var)

#=
arithmetic =#
Base.:-(p::P) where {P <: AbstractPolynomial} = _convert(p, -coeffs(p))

Base.:*(p::AbstractPolynomial, c::Number) = scalar_mult(p, c)
Base.:*(c::Number, p::AbstractPolynomial) = scalar_mult(c, p)
Base.:*(c::T, p::P) where {T, P <: AbstractPolynomial{T}} = scalar_mult(c, p)
Base.:*(p::P, c::T) where {T, P <: AbstractPolynomial{T}} = scalar_mult(p, c)

# implicitly identify c::Number with a constant polynomials
Base.:+(c::Number, p::AbstractPolynomial) = +(p, c)
Base.:-(p::AbstractPolynomial, c::Number) = +(p, -c)
Base.:-(c::Number, p::AbstractPolynomial) = +(-p, c)

# scalar operations
# no generic p+c, as polynomial addition falls back to scalar ops
#function Base.:+(p::P, n::Number) where {P <: AbstractPolynomial}
#    p1, p2 = promote(p, n)
#    return p1 + p2
#end

Base.:-(p1::AbstractPolynomial, p2::AbstractPolynomial) = +(p1, -p2)

## addition
## Fall back addition is possible as vector addition with padding by 0s
## Subtypes will likely want to implement both:
## +(p::P,c::Number) and +(p::P, q::Q) where {T,S,X,P<:SubtypePolynomial{T,X},Q<:SubtypePolynomial{S,X}}
## though the default for poly+poly isn't terrible

# polynomial + scalar; implicit identification of c with c*one(P)
Base.:+(p::P, c::T) where {T,X, P<:AbstractPolynomial{T,X}} = p + c * one(P)

function Base.:+(p::P, c::S) where {T,X, P<:AbstractPolynomial{T,X}, S}
    R = promote_type(T,S)
    q = convert(⟒(P){R,X}, p)
    q + R(c)
end

# polynomial + polynomial when different types
function Base.:+(p::P, q::Q) where {T,X,P <: AbstractPolynomial{T,X}, S,Y,Q <: AbstractPolynomial{S,Y}}
    isconstant(p) && return constantterm(p) + q
    isconstant(q) && return p + constantterm(q)
    assert_same_variable(X,Y)
    sum(promote(p,q))

end

# Works when p,q of same type.
# For Immutable, must remove N,M bit;
# for others, can widen to Type{T,X}, Type{S,X} to avoid a promotion
function Base.:+(p::P, q::P) where {T,X,P<:AbstractPolynomial{T,X}}
    cs = degree(p) >= degree(q)  ? ⊕(P, p.coeffs, q.coeffs) : ⊕(P, q.coeffs, p.coeffs)
    return P(cs)
end

# addition of polynomials is just vector space addition, so can be done regardless
# of basis, as long as the same. These ⊕ methods try to find a performant means to add
# to sets of coefficients based on the storage type. These assume n1 >= n2
function ⊕(P::Type{<:AbstractPolynomial}, p1::Vector{T}, p2::Vector{S}) where {T,S}

    n1, n2 = length(p1), length(p2)
    R = promote_type(T,S)

    cs = collect(R,p1)
    for i in 1:n2
        cs[i] += p2[i]
    end

    return cs
end

# Padded vector sum of two tuples assuming N ≥ M
@generated function ⊕(P::Type{<:AbstractPolynomial}, p1::NTuple{N,T}, p2::NTuple{M,S}) where {T,N,S,M}

    exprs = Any[nothing for i = 1:N]
    for i in  1:M
        exprs[i] = :(p1[$i] + p2[$i])
    end
    for i in M+1:N
        exprs[i] =:(p1[$i])
    end

    return quote
        Base.@_inline_meta
        tuple($(exprs...))
    end

end

# addition when a dictionary is used for storage
function ⊕(P::Type{<:AbstractPolynomial}, p1::Dict{Int,T}, p2::Dict{Int,S}) where {T,S}

    R = promote_type(T,S)
    p = Dict{Int, R}()

    # this allocates in the union
#    for i in union(eachindex(p1), eachindex(p2))
#        p[i] = p1[i] + p2[i]
#    end

    for (i,pi) ∈ pairs(p1)
        @inbounds p[i] = pi + get(p2, i, zero(R))
    end
    for (i,pi) ∈ pairs(p2)
        if iszero(get(p,i,zero(R)))
            @inbounds p[i] = get(p1, i, zero(R)) + pi
        end
    end

    return  p

end

## -- multiplication

# this fall back not necessarily efficient (e.g., sparse)
function scalar_mult(p::P, c::S) where {S, T, X, P<:AbstractPolynomial{T,X}}
    R = Base.promote_op(*, T, S) # typeof(one(T)*one(S))?
    𝐏 = ⟒(P){R,X}
    𝐏([pᵢ * c for pᵢ ∈ coeffs(p)])
end

function scalar_mult(c::S, p::P) where {S, T, X, P<:AbstractPolynomial{T, X}}
    R = Base.promote_op(*, T, S)
    𝐏 = ⟒(P){R,X}
    𝐏([c * pᵢ for pᵢ ∈ coeffs(p)])
end

function Base.:/(p::P, c::S) where {P <: AbstractPolynomial,S}
    _convert(p, coeffs(p) ./ c)
end

## polynomial p*q
## Polynomial multiplication formula depend on the particular basis used. The subtype must implement
function Base.:*(p1::P, p2::Q) where {T,X,P <: AbstractPolynomial{T,X},S,Y,Q <: AbstractPolynomial{S,Y}}
    isconstant(p1) && return constantterm(p1) * p2
    isconstant(p2) && return p1 * constantterm(p2)
    assert_same_variable(X, Y)
    p1, p2 = promote(p1, p2)
    return p1 * p2
end

Base.:^(p::AbstractPolynomial, n::Integer) = Base.power_by_squaring(p, n)

function Base.divrem(num::P, den::O) where {P <: AbstractPolynomial,O <: AbstractPolynomial}
    n, d = promote(num, den)
    return divrem(n, d)
end

"""
    gcd(a::AbstractPolynomial, b::AbstractPolynomial; atol::Real=0, rtol::Real=Base.rtoldefault)

Find the greatest common denominator of two polynomials recursively using
[Euclid's algorithm](http://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclid.27s_algorithm).

# Examples

```jldoctest common
julia> using Polynomials

julia> gcd(fromroots([1, 1, 2]), fromroots([1, 2, 3]))
Polynomial(4.0 - 6.0*x + 2.0*x^2)

```
"""
function Base.gcd(p1::AbstractPolynomial{T}, p2::AbstractPolynomial{S}; kwargs...) where {T,S}
    gcd(promote(p1, p2)...; kwargs...)
end

function Base.gcd(p1::AbstractPolynomial{T}, p2::AbstractPolynomial{T};
                  atol::Real=zero(real(T)),
                  rtol::Real=Base.rtoldefault(real(T))
                  ) where {T}

    r₀, r₁ = p1, p2
    iter = 1
    itermax = length(r₁)

    while !iszero(r₁) && iter ≤ itermax
        _, rtemp = divrem(r₀, r₁)
        r₀ = r₁
        r₁ = truncate(rtemp; atol=atol, rtol=rtol)
        iter += 1
    end
    return r₀
end

"""
    uvw(p,q; kwargs...)

return `u` the gcd of `p` and `q`, and `v` and `w`, where `u*v = p` and `u*w = q`.
"""
uvw(p::AbstractPolynomial, q::AbstractPolynomial; kwargs...) = uvw(promote(p,q)...; kwargs...)
uvw(p1::P, p2::P; kwargs...) where {P <:AbstractPolynomial} = throw(ArgumentError("uvw not defind"))

"""
    div(::AbstractPolynomial, ::AbstractPolynomial)
"""
Base.div(n::AbstractPolynomial, d::AbstractPolynomial) = divrem(n, d)[1]

"""
    rem(::AbstractPolynomial, ::AbstractPolynomial)
"""
Base.rem(n::AbstractPolynomial, d::AbstractPolynomial) = divrem(n, d)[2]

#=
Comparisons =#
Base.isequal(p1::P, p2::P) where {P <: AbstractPolynomial} = hash(p1) == hash(p2)
Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial) =
    check_same_variable(p1,p2) && (coeffs(p1) == coeffs(p2))
Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && constantterm(p) == n
Base.:(==)(n::Number, p::AbstractPolynomial) = p == n

function Base.isapprox(p1::AbstractPolynomial{T,X},
                       p2::AbstractPolynomial{S,Y};
                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {T,X,S,Y}
    assert_same_variable(p1, p2)
    (hasnan(p1) || hasnan(p2)) && return false  # NaN poisons comparisons
    # copy over from abstractarray.jl
    Δ  = norm(p1-p2)
    if isfinite(Δ)
        return Δ <= max(atol, rtol*max(norm(p1), norm(p2)))
    else
        for i in 0:max(degree(p1), degree(p2))
            isapprox(p1[i], p2[i]; rtol=rtol, atol=atol) || return false
        end
        return true
    end
end

function Base.isapprox(p1::AbstractPolynomial{T},
                       n::S;
                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {T,S}
    return isapprox(p1, _convert(p1, [n]))
end

Base.isapprox(n::S,
    p1::AbstractPolynomial{T};
    rtol::Real = (Base.rtoldefault(T, S, 0)),
    atol::Real = 0,) where {T,S} = isapprox(p1, n, rtol = rtol, atol = atol)


text:
module TestMiscSDP

import Mosek
using JuMP,FactCheck

facts("[sdp] Robust uncertainty example") do
    solver = Mosek.MosekSolver(QUIET=true)

    context("With solver $(typeof(solver))") do
        include(joinpath(Pkg.dir("JuMP"), "test","data","robust_uncertainty.jl"))
        R = 1
        d = 3
        𝛿 = 0.05
        ɛ = 0.05
        N = ceil((2+2log(2/𝛿))^2) + 1

        Γ1(𝛿,N) = (R/sqrt(N))*(2+sqrt(2*log(1/𝛿)))
        Γ2(𝛿,N) = (2R^2/sqrt(N))*(2+sqrt(2*log(2/𝛿)))

        #for d in [3,5,8]; context("d = $d") do
        for d in [3]; context("d = $d") do
            μhat = μhats[d]
            M = Ms[d]
            Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat')

            m = Model(solver=solver)

            @defVar(m, Σ[1:d,1:d], SDP)
            @defVar(m, u[1:d])
            @defVar(m, μ[1:d])

            @defVar(m, t1 >= 0)
            @defVar(m, L1[1:d])
            @addConstraint(m, L1 .== (μ-μhat))
            @addConstraint(m, sum{L1[i]^2, i=1:d} <= t1^2)
            @addConstraint(m, t1 <= Γ1(𝛿/2,N))

            @defVar(m, t2 >= 0)
            @defVar(m, L2[1:d,1:d])
            @addConstraint(m, L2 .== (Σ-Σhat))
            @addConstraint(m, sum{L2[i,j]^2, i=1:d, j=1:d} <= t2^2)
            @addConstraint(m, t2 <= Γ2(𝛿/2,N))

            A = [(1-ɛ)/ɛ (u-μ)';
                 (u-μ)     Σ   ]
            @addSDPConstraint(m, A >= 0)

            c = cs[d]
            @setObjective(m, Max, dot(c,u))

            stat = solve(m)

            object = getObjectiveValue(m)
            exact = dot(μhat,c) + Γ1(𝛿/2,N)*norm(c) + sqrt((1-ɛ)/ɛ)*sqrt(dot(c,(Σhat+Γ2(𝛿/2,N)*eye(d,d))*c))
            @fact stat --> :Optimal
            @fact abs(object - exact) --> roughly(0, 1e-5)

            resΣ  = getValue(Σ)
            resu  = getValue(u)
            resμ  = getValue(μ)
            rest1 = getValue(t1)
            resL1 = getValue(L1)
            rest2 = getValue(t2)
            resL2 = getValue(L2)

            println("object = $object, exact = $exact")
            println("object val : dot(c,u) : $(dot(c,resu))")

            println("L1 == μ-μhat : $(resL1 - (resμ-μhat))")
            println("sum{L1[i]^2, i=1:d} <= t1^2: $((sum(resL1.^2)) - (rest1.^2))")
            println("t1 <= Γ1(𝛿/2,N) : $(rest1 - (Γ1(𝛿/2,N)))")
            println("L2 .== (Σ-Σhat) : $(resL2 - (resΣ-Σhat))")
            println("sum{L2[i,j]^2, i=1:d, j=1:d} <= t2^2 : $(sum(resL2.^2)-rest2^2)")
            #@addConstraint(m, sum{L2[i,j]^2, i=1:d, j=1:d} <= t2^2)
            println("t2 <= Γ2(𝛿/2,N) : $(rest2 - Γ2(𝛿/2,N))")
            #@addConstraint(m, t2 <= Γ2(𝛿/2,N))
            println("---------- PSD con")
            resA = [(1-ɛ)/ɛ     (resu-resμ)';
                    (resu-resμ) resΣ   ]

            println("Eigvals A: $(eigvals(resA))")
            println("---------- vars")
            println("Σ PSD. Eigvals : $(eigvals(resΣ))")
            println("u free")
            println("μ free")
            println("t1 >= 0: $(rest1)")
            println("L1 free")
            println("t2 >= 0: $(rest2)")
            println("L2 free")

            intmodel = getInternalModel(m)
            Mosek.writedata(intmodel.task,"problem-$d.task")
        end
        end
    end

end
end

module TestX
using Mosek, MathProgBase

I = [1,1,2,2,3,3,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,32,32,33,33,34,34,35,36,37,38,39,40]
J = [14,10,15,11,16,12,13,18,1,21,2,24,3,19,2,22,4,25,5,20,3,23,5,26,6,17,13,17,13,14,15,16,17,18,19,20,21,22,23,24,25,26,7,10,8,11,9,12,1,2,3,4,5,6]
V = [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.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.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.0,-1.0,-1.0,-1.0,-1.0,-1.0]
b = [-0.6088331477422722,-0.02852197818671054,-0.9016074145524158,0.4999165215049682,-0.025485683857737418,-0.07276320066249353,0.056565030662414396,-0.07276320066249353,-0.3319954143963322,0.3286736785015988,0.056565030662414396,0.3286736785015988,-0.43682537007384287,1.0516057442061533,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,18.999999999999996,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
var_cones = Any[(:SDP,1:6),(:NonNeg,[13,17]),(:Free,[7,8,9,10,11,12,14,15,16,18,19,20,21,22,23,24,25,26])]
con_cones = Any[(:NonNeg,[4,14]),(:NonPos,[15,16]),(:Zero,[1,2,3,5,6,7,8,9,10,11,12,13]),(:SOC,17:20),(:SOC,21:30),(:SDP,31:40)]
f = [-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,0.2269228958615089,-0.26647992165696055,-0.35104470584603914,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-0.0]
numconstr = 40
numvar = 26
A = sparse(I,J,V,numconstr,numvar)

m = MathProgBase.ConicModel(MosekSolver(LOG=0))

MathProgBase.loadproblem!(m, f, A, b, con_cones, var_cones)
writedata(m.task,"problem-back1.task")
println("varbk before : $(m.varbk) $(m.varmap)")
MathProgBase.setvartype!(m, fill(:Cont,numvar))
println("varbk after  : $(m.varbk) $(m.varmap)")
writedata(m.task,"problem-back2.task")

MathProgBase.optimize!(m)

@assert MathProgBase.status(m) == :Optimal

obj = MathProgBase.getobjval(m)
trueobj = -2.7231533727479005
@show trueobj
@show obj
@show abs(obj-trueobj)

writedata(m.task,"problem-back.task")

end


text:
Base.:*(u::AbstractVector, v::SVector{1,T}) where T = u * v[1]

Base.vcat(v::Adjoint{T,SVector{N,T}}, E::UniformScaling) where {T, N} = vcat(v, SMatrix{N,N,T,N*N}(E))
Base.vcat(A::SMatrix{N1,N2,T,N1N2}, E::UniformScaling) where {T, N1, N2, N1N2} = vcat(A, SMatrix{N2,N2,T,N2*N2}(E))

svcat(a::T, b::T) where T = SA[a; b]
svcat(a::StaticArray, b::StaticArray) = vcat(a, b)
svcat(a::StaticArray{Tuple{N},T,1}, b::T) where {T,N} = vcat(a, SA[b])
svcat(a::T, b::StaticArray{Tuple{N},T,1}) where {T,N} = vcat(SA[a], b)

szeros(::Type{T}, N) where T = @SVector zeros(T, N)
szeros(N)= @SVector zeros(N)
szeros(::Type{T}, N1, N2) where T = @SMatrix zeros(T, N1, N2)
szeros(N1, N2)= @SMatrix zeros(N1, N2)

sones(::Type{T}, N) where T = @SVector ones(T, N)
sones(N)= @SVector ones(N)

srand(N)= @SVector rand(N)

# TODO: check StaticArray bug fix, then remove
zerodimstaticadjoint(A) = A'
zerodimstaticadjoint(::SMatrix{0,N,T,0}) where {T,N} = SMatrix{N,0,T,0}()

function diagonal_cat(a::AbstractMatrix{T}, b::AbstractMatrix{T}) where T
    diagonal_cat(dense(a), dense(b))
end

function diagonal_cat(a::SMatrix{Na,Na,T,Ma}, b::SMatrix{Nb,Nb,T,Mb}) where {T,Na,Nb,Ma,Mb}
    [[a szeros(Na,Nb)]; [szeros(Nb,Na) b]]
end

dense(a::AbstractMatrix{T}) where T = a
function dense(a::Diagonal{T,SVector{N, T}}) where {T,N}
    SMatrix{N,N,T,N^2}(a)
end


text:
#
# LocalFilters.jl --
#
# Local filters for Julia arrays.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2020, Éric Thiébaut.
#

module LocalFilters

export
    Neighborhood,
    RectangularBox,
    bilateralfilter!,
    bilateralfilter,
    bottom_hat,
    closing!,
    closing,
    convolve!,
    convolve,
    dilate!,
    dilate,
    erode!,
    erode,
    localextrema!,
    localextrema,
    localfilter!,
    localfilter,
    localmean!,
    localmean,
    opening!,
    opening,
    strel,
    top_hat

import Base: CartesianIndices,
    axes, eltype, ndims, size, length, first, last, tail,
    getindex, setindex!, convert, reverse

function localfilter end
function localfilter! end

include("types.jl")
include("basics.jl")
include("filters.jl")
include("morphology.jl")
include("separable.jl")
include("bilateral.jl")
import .BilateralFilter: bilateralfilter!, bilateralfilter

end


text:
# ---
# title: 483. Smallest Good Base
# id: problem483
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/smallest-good-base/description/>
# hidden: true
# ---
#
# For an integer n, we call k>=2 a _**good base**_ of n, if all digits of n base
# k are 1.
#
# Now given a string representing n, you should return the smallest good base of
# n in string format.
#
# **Example 1:**
#
#
#
#     Input: "13"
#     Output: "3"
#     Explanation: 13 base 3 is 111.
#
#
#
#
# **Example 2:**
#
#
#
#     Input: "4681"
#     Output: "8"
#     Explanation: 4681 base 8 is 11111.
#
#
#
#
# **Example 3:**
#
#
#
#     Input: "1000000000000000000"
#     Output: "999999999999999999"
#     Explanation: 1000000000000000000 base 999999999999999999 is 11.
#
#
#
#
# **Note:**
#
#   1. The range of n is [3, 10^18].
#   2. The string representing n is always valid and will not have leading zeros.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using NGWP, Plots, LightGraphs, JLD, MTSG, MAT

barbara = JLD.load("../datasets/barbara_gray_matrix.jld", "barbara")

## Build weighted graph
G, L, X = SunFlowerGraph(N = 400); N = nv(G)
W = 1.0 * adjacency_matrix(G)
if runapprox
    Q = incidence_matrix(G; oriented = true)
    edge_weight = [e.weight for e in edges(G)]

    ## eigenvectors of L(G)
    𝛌, 𝚽 = eigen(Matrix(L))
    standardize_eigenvectors!(𝚽)

    ## eigenvectors of Lsym(G)
    deg = sum(W, dims = 1)[:]  # weighted degree vector
    Lsym = diagm(deg.^(-1/2)) * (diagm(deg) - W) * diagm(deg.^(-1/2))
    𝛌sym, 𝚽sym = eigen(Lsym)
    standardize_eigenvectors!(𝚽sym)

    ## Build Dual Graph by DAG metric
    distDAG = eigDAG_Distance(𝚽, Q, N; edge_weight = edge_weight)
    Gstar_Sig = dualgraph(distDAG)
    G_Sig = GraphSig(W, xy = X)
    GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
    GP_primal = pairclustering(𝚽, GP_dual)
    jmax = size(GP_dual.rs, 2) - 1  # zero-indexed

    if allNGWPs
        # 54.986524 seconds (1.19 M allocations: 25.127 GiB, 2.36% gc time)
        @time VM_NGWP = vm_ngwp(𝚽, GP_dual)
        # 0.611176 seconds (225.41 k allocations: 844.488 MiB, 11.71% gc time)
        @time PC_NGWP = pc_ngwp(𝚽, GP_dual, GP_primal)
        # 50.939912 seconds (7.67 M allocations: 28.051 GiB, 3.11% gc time)
        @time LP_NGWP = lp_ngwp(𝚽, Gstar_Sig.W, GP_dual; ϵ = 0.3)
    end

    ## Build Dual Graph by DAG metric (Lsym)
    distDAG_Lsym = eigDAG_Distance(𝚽sym, Q, N; edge_weight = edge_weight)
    Gstar_Sig_Lsym = dualgraph(distDAG_Lsym)
    GP_dual_Lsym = partition_tree_fiedler(Gstar_Sig_Lsym; swapRegion = false)
    GP_primal_Lsym = pairclustering(𝚽sym, GP_dual_Lsym)
    jmax_Lsym = size(GP_dual_Lsym.rs, 2) - 1

    if allNGWPs
        # 61.328441 seconds (1.34 M allocations: 29.058 GiB, 2.30% gc time)
        @time VM_NGWP_Lsym = vm_ngwp(𝚽sym, GP_dual_Lsym)
        # 0.605507 seconds (237.87 k allocations: 876.689 MiB, 10.83% gc time)
        @time PC_NGWP_Lsym = pc_ngwp(𝚽sym, GP_dual_Lsym, GP_primal_Lsym)
        # 77.918189 seconds (2.12 M allocations: 42.872 GiB, 2.91% gc time)
        @time LP_NGWP_Lsym = lp_ngwp(𝚽sym, Gstar_Sig_Lsym.W, GP_dual_Lsym; ϵ = 0.3)
    end
else
    using VoronoiDelaunay, VoronoiCells, GeometricalPredicates
    ## Voronoi tessellation
    width_x = maximum(abs.(X[:, 1])) * 2; width_y = maximum(abs.(X[:, 2])) * 2;
    width = VoronoiDelaunay.max_coord - VoronoiDelaunay.min_coord
    center_coord = (VoronoiDelaunay.min_coord + VoronoiDelaunay.max_coord)/2
    X_transform = zeros(N,2)
    for i in 1:N
        X_transform[i,:] = X[i,:] ./ [width_x/width, width_y/width] + [center_coord, center_coord]
    end

    pts = [Point2D(X_transform[i,1], X_transform[i,2]) for i in 1:N]
    tess = DelaunayTessellation(N)
    push!(tess, pts)
end


text:
export inv_rotation

################################################################################
#                              Inverse Rotations
################################################################################

"""
    @inline function inv_rotation(R)

Compute the inverse rotation of `R`, which can be a Direction Cosine Matrix or
Quaternion.

The output will have the same type as `R` (DCM or quaternion).

# Remarks

If `R` is a DCM, than its transpose is computed instead of its inverse to reduce
the computational burden. The both are equal if the DCM has unit norm. This must
be verified by the used.

If `R` is a quaternion, than its conjugate is computed instead of its inverse to
reduce the computational burden. The both are equal if the quaternion has unit
norm. This must be verified by the used.

# Example

```julia-repl
julia> D = angle2dcm(+pi/3,+pi/4,+pi/5,:ZYX);

julia> inv_rotation(D)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  0.353553  -0.492816  0.795068
  0.612372   0.764452  0.201527
 -0.707107   0.415627  0.572061

julia> q = angle2quat(+pi/3,+pi/4,+pi/5,:ZYX);

julia> inv_rotation(q)
Quaternion{Float64}:
  + 0.8200711519756747 - 0.06526868310243991.i - 0.45794027732580056.j - 0.336918398289752.k
```

"""
@inline inv_rotation(D::DCM) = D'
@inline inv_rotation(q::Quaternion) = conj(q)