text:
export is_homogeneous,
    homogenize,
    uniquevar,
    Composition,
    expand,
    compose,
    validate,
    precondition,
    normalize_coefficients,
    linear_system

const MPPoly{T} = MP.AbstractPolynomialLike{T}
const MPPolys = AbstractVector{<:MP.AbstractPolynomialLike}
const WeightedVariable = Tuple{<:MP.AbstractVariable,Int}

##############
# COMPOSITION
##############

abstract type AbstractComposition end

"""
    Composition

A `Composition` is a composition of polynomial systems. This is the result of [`compose`](@ref).
"""
mutable struct Composition{T<:MP.AbstractPolynomialLike} <: AbstractComposition
    polys::Vector{Vector{T}} # polys = [f1, f2, f3] -> f1 ∘ f2 ∘ f3
end

Base.length(C::Composition) = length(C.polys[1])
Base.:(==)(C::Composition, D::Composition) = C.polys == D.polys

"""
    compose(g, f)::Composition

Compose the polynomial systems `g` and `f`.
You can also use the infix operator `∘` (written by \\circ).

```julia-repl
julia> @polyvar a b c x y z;
julia> g = [a * b * c];
julia> f = [x+y, y + z, x + z];
julia> expand(compose(g, f))
1-element Array{DynamicPolynomials.Polynomial{true,Int64},1}:
 x²y + x²z + xy² + 2xyz + xz² + y²z + yz²
```
"""
compose(g::Composition, f::Composition) = Composition([g.polys..., f.polys...])
compose(g::Composition, f::Composition, h::Composition...) =
    compose(Composition([g.polys..., f.polys...]), h...)
compose(g::Composition, fs::MPPolys...) = Composition([g.polys..., fs...])
compose(g::MPPolys, f::Composition) = Composition([g, f.polys...])
compose(g::MPPolys, f::Composition, h::Composition...) =
    compose(Composition([g, f.polys...]), h...)
compose(g::MPPolys, fs::MPPolys...) = Composition([g, fs...])

import Base: ∘
∘(g::Union{Composition,<:MPPolys}, f::MPPolys) = compose(g, f)
∘(g::Union{Composition,<:MPPolys}, f::MPPolys...) = compose(g, f...)
∘(g::Union{Composition,<:MPPolys}, f::Composition...) = compose(g, f...)

"""
    scale(C::Composition, λ)

Scale a composition by λ.
"""
scale(C::Composition, λ) = Composition([[C.polys[1] .* λ]; C.polys[2:end]])

function MP.subs(C::Composition, args::MP.AbstractSubstitution...)
    p = map(pᵢ -> MP.subs(pᵢ, args...), C.polys[end])
    Composition([C.polys[1:end-1]; [p]])
end

"""
    expand(C::Composition; parameters=nothing)

Expand the composition to the polynomial system is originally represents.

## Example
```julia-repl
julia> @polyvar a b c x y z;
julia> f = [a * b * c];
julia> g = [x+y, y + z, x + z];
julia> expand(f ∘ g)
1-element Array{DynamicPolynomials.Polynomial{true,Int64},1}:
 x²y + x²z + xy² + 2xyz + xz² + y²z + yz²
```
"""
function expand(C::Composition; parameters = nothing)
    g = nothing
    for f in reverse(C.polys)
        if g === nothing
            g = f
            continue
        end
        vars = variables(f, parameters)
        g = map(fᵢ -> MP.subs(fᵢ, vars => g), f)
    end
    g
end

"""
    validate(C::Composition; parameters=nothing)

Validates that the composition is well defined.
"""
function validate(C::Composition; parameters = nothing)
    g = nothing
    for f in reverse(C.polys)
        if g === nothing
            g = f
            continue
        end
        nvars = length(variables(f, parameters))
        if nvars !== length(g)
            return false
        end
        g = f
    end
    true
end

#############################
# HomogenizationInformation
#############################

"""
	HomogenizationInformation(;homvar=nothing, homvars=nothing, variable_groups=nothing)

This captures the information about the desired homogenization of the problem.
`homvar(s) = nothing` indicates that the problem should be homogenized.
`variable_groups` indicate the desired variable grouping for multi-homogeneous homotopies.

## Fields
* `homvars`
* `vargroups`
"""
struct HomogenizationInformation{N,T}
# one of the two fields is always not Nothing
    homvars::Union{Nothing,NTuple{N,T}}
    vargroups::Union{Nothing,NTuple{N,Vector{T}}}
end

function HomogenizationInformation(;
    homvar = nothing,
    homvars = nothing,
    variable_groups = nothing,
)
    if homvars === nothing && homvar !== nothing
        homvars = (homvar,)
        if variable_groups !== nothing && length(variable_groups) > 1
            error("You provided more than 1 variables group but only declared `homvar=$homvar`. Use instead `homvars=...`.")
        end
    end

    if variable_groups === nothing && homvars === nothing
        return nothing
    elseif variable_groups === nothing && length(homvars) > 1
        error("Currently we cannot find variable groups just from given `homvars`. Please provide explicit variable groups using `variable_groups=...`.")
    elseif variable_groups === nothing
        HomogenizationInformation(homvars, nothing)
    else
        HomogenizationInformation(homvars, tuple(collect.(variable_groups)...))
    end
end

function add_variable_groups(
    HI::HomogenizationInformation,
    F::MPPolys;
    parameters = nothing,
)
    if HI.vargroups === nothing
        if length(HI.homvars) == 1
            vargroups = (variables(F, parameters),)
            return HomogenizationInformation(HI.homvars, vargroups)
        else
            error("Cannot add variable groups")
        end
    end
    HI
end
add_variable_groups(::Nothing, F; parameters = nothing) = nothing

homvars(H::HomogenizationInformation) = H.homvars
homvars(::Nothing) = nothing

"""
    VariableGroups{N, Homvars, V}

A `VariableGroups` stores an `NTuple` of indices mapping to the indices of the original variables
of type `V`. `Homvars` is true if there have been homogenization variables declared.
"""
struct VariableGroups{N,Homvars,V<:Union{MP.AbstractVariable,Int}}
    variables::Vector{V}
    groups::NTuple{N,Vector{Int}}
    pull_back_mapping::Vector{Tuple{Int,Int}}
end

function VariableGroups(
    variables,
    groups::NTuple{N,Vector{Int}},
    dedicated_homvars::Bool,
) where {N}
    hom_vars = last.(groups)
    elements = vcat(groups...)
    pull_back_mapping = Tuple{Int,Int}[]
    proj_homvar_indices = cumsum([length(g) for g in groups])
    for i = 1:sum(length, groups)
        i ∉ hom_vars || continue
        l = 0
        for (gᵢ, group) in enumerate(groups), k in group
            l += 1
            if i == k
                push!(pull_back_mapping, (l, proj_homvar_indices[gᵢ]))
                break
            end
        end
    end
    VariableGroups{N,dedicated_homvars,eltype(variables)}(
        variables,
        groups,
        pull_back_mapping,
    )
end

"""
    VariableGroups(all_variables, variable_groups, homvars)

# Example
```julia
all_vars = @polyvar a b c x y z
variable_groups = ([c, b, a], [x, y, z])
homvars = (b, y)

vargroups = VariableGroups(all_vars, variable_groups, homvars)
vargroups.groups == ([3, 1, 2], [4, 6, 5])
vargroups.homvars == true
```
"""
function VariableGroups(
    all_variables,
    variable_groups::NTuple{N,<:Vector{<:Union{Int,MP.AbstractVariable}}},
    homvars::Union{Nothing,NTuple{N,<:Union{Int,MP.AbstractVariable}}} = nothing,
) where {N}
    groups = ntuple(N) do k
        variable_group = variable_groups[k]
        h = homvars === nothing ? 0 : homvars[k]
        group = Int[]
        homvar = 0
        for v in variable_group
            i = findfirst(w -> w == v, all_variables)
            if v == h
                homvar = i
            else
                push!(group, i)
            end
        end
        if homvar != 0
            push!(group, homvar)
        end
        group
    end
    VariableGroups(all_variables, groups, homvars !== nothing)
end
function VariableGroups(all_variables, homvar_info::HomogenizationInformation)
    if homvar_info.vargroups === nothing
        VariableGroups(all_variables, (collect(all_variables),), homvar_info.homvars)
    else
        VariableGroups(all_variables, homvar_info.vargroups, homvar_info.homvars)
    end
end

function VariableGroups(all_variables, homvar::MP.AbstractVariable)
    VariableGroups(all_variables, (collect(all_variables),), (homvar,))
end
function VariableGroups(all_variables, homvar::Nothing = nothing)
    VariableGroups(all_variables, (collect(all_variables),), homvar)
end
function VariableGroups(nvariables::Int, homvar::Nothing = nothing)
    VariableGroups(collect(1:nvariables), (collect(1:nvariables),), nothing)
end
function VariableGroups(nvariables::Int, homvar::Int)
    VariableGroups(collect(1:nvariables), (collect(1:nvariables),), (homvar,))
end
function VariableGroups(nvariables::Int, hominfo::HomogenizationInformation)
    if hominfo.vargroups === nothing
        VariableGroups(collect(1:nvariables), (collect(1:nvariables),), hominfo.homvars)
    else
        VariableGroups(collect(1:nvariables), hominfo.vargroups, hominfo.homvars)
    end
end

"""
	has_dedicated_homvars(variable_groups)

Returns `true` if there have been homogenization variables declared.
"""
has_dedicated_homvars(groups::VariableGroups{N,Homvars}) where {N,Homvars} = Homvars

"""
	projective_dims(variable_groups)

Returns the projective dimension of each variable group.
"""
function projective_dims(groups::VariableGroups)
    length.(groups.groups) .- 1
end

nvariables(G::VariableGroups) = sum(length.(G.groups))
Base.length(::VariableGroups{N}) where {N} = N

"""
	variable_groups(VG::VariableGroups)

Returns the variable groups.
"""
variable_groups(VG::VariableGroups) = map(g -> variables(VG)[g], VG.groups)

"""
	flattened_groups(VG::VariableGroups)

Returns the variable groups.
"""
flattened_variable_groups(VG::VariableGroups) = vcat(variable_groups(VG)...)

"""
	variables(VG::VariableGroups)

Returns all variables in their original order.
"""
variables(VG::VariableGroups) = VG.variables

"""
	ngroups(VG::VariableGroups)

Group the given variables in their corresponding groups.
"""
ngroups(VG::VariableGroups{M}) where {M} = M

function embed_projective(VG::VariableGroups{M}, v::AbstractVector) where {M}
    dims = projective_dims(VG)
    x = PVector(zeros(eltype(v), sum(dims) + M), dims)
    embed_projective!(x, VG, v)
end

function embed_projective!(
    x::PVector{<:Number,M},
    VG::VariableGroups{M},
    v::PVector{<:Number,M},
) where {M}
    x.data .= v.data
    x
end
function embed_projective!(
    x::PVector{<:Number,M},
    VG::VariableGroups{M},
    v::AbstractVector,
) where {M}
    if sum(projective_dims(VG)) == length(v)
        k = 1
        for group in VG.groups
            for i = 1:(length(group)-1)
                x[k] = v[group[i]]
                k += 1
            end
            x[k] = one(eltype(x))
            k += 1
        end
    else
        k = 1
        for group in VG.groups
            for i in group
                x[k] = v[i]
                k += 1
            end
        end
    end
    x
end

##############
# POLYNOMIALS
#############
function variables(
    polys::Union{MPPoly,MPPolys},
    params::AbstractVector{<:MP.AbstractVariable},
)
    setdiff!(variables(polys), params)
end
variables(C::Composition) = variables(C.polys[end])
variables(C::Composition, p::AbstractVector{<:MP.AbstractVariable}) =
    variables(C.polys[end], p)
variables(p::Union{MPPoly,MPPolys,Composition}, ::Nothing) = variables(p)

"""
    hasparameters(F, parameters)

Returns `true` if the parameters occur in F.
"""
function hasparameters(polys::Union{MPPoly,MPPolys}, parameters = nothing)
    isnothing(parameters) && return false
    vars = variables(polys)
    any(p -> p in vars, parameters)
end

"""
    nvariables(polys, parameters)

Returns the number of variables occuring in `polys`.
"""
nvariables(F::MPPolys, p::AbstractVector{<:MP.AbstractVariable}) = length(variables(F, p))
nvariables(C::Composition, p = nothing) = length(variables(C, p))

"""
    is_homogeneous(f::MP.AbstractPolynomialLike)

Checks whether `f` is homogeneous.

    is_homogeneous(f::MP.AbstractPolynomialLike, vars)

Checks whether `f` is homogeneous in the variables `vars` with possible weights.
"""
is_homogeneous(f::MP.AbstractPolynomialLike) = MP.mindegree(f) == MP.maxdegree(f)
function is_homogeneous(f::MP.AbstractPolynomialLike, variables)
    d_min, d_max = minmaxdegree(f, variables)
    d_min == d_max
end

"""
    is_homogeneous(F::Vector{MP.AbstractPolynomialLike}, variables)

Checks whether each polynomial in `F` is homogeneous in the variables `variables`.
"""
function is_homogeneous(F::MPPolys, variables)
    all(f -> is_homogeneous(f, variables), F)
end
function is_homogeneous(F::MPPolys, ::Nothing = nothing; parameters = nothing)
    if parameters !== nothing
        is_homogeneous(F, variables(F, parameters))
    else
        all(is_homogeneous, F)
    end
end
function is_homogeneous(
    F::MPPolys,
    hominfo::HomogenizationInformation;
    parameters = nothing,
)
    if hominfo.vargroups !== nothing
        all(vars -> is_homogeneous(F, vars), hominfo.vargroups)
    elseif parameters !== nothing
        is_homogeneous(F, variables(F, parameters))
    else
        all(is_homogeneous, F)
    end
end
function is_homogeneous(C::Composition, hominfo = nothing; kwargs...)
    homogeneous_degrees_helper(C; kwargs...) !== nothing
end
function is_homogeneous(C::Composition, hominfo::HomogenizationInformation; kwargs...)
    if hominfo.vargroups && length(hominfo.vargroups) > 1
        error("`variable_groups` and compositions are currently not supported.")
    end

    homogeneous_degrees_helper(C; kwargs...) !== nothing
end

function homogeneous_degrees_helper(C::Composition; parameters = nothing, weights = nothing)
    for f in reverse(C.polys)
        weights = homogeneous_degrees_helper(f, parameters = parameters, weights = weights)
        weights === nothing && return nothing
    end
    weights
end
function homogeneous_degrees_helper(F::MPPolys; parameters = nothing, weights = nothing)
    vars = variables(F, parameters)
    if !isnothing(weights)
        vars = zip(vars, weights)
    end
    degrees = Int[]
    for f in F
        mindeg, maxdeg = minmaxdegree(f, vars)
        if mindeg != maxdeg
            return nothing
        else
            push!(degrees, mindeg)
        end
    end
    degrees
end

"""
    degree(term::MP.AbstractTermLike, vars)

Compute the (weighted) degree of `f` in the variables `vars`.
"""
function degree(term::MP.AbstractTermLike, variables::Vector{<:MP.AbstractVariable})
    sum(MP.degree(term, v) for v in variables)
end
function degree(
    term::MP.AbstractTermLike,
    variables::NTuple{N,<:MP.AbstractVariable},
) where {N}
    sum(MP.degree(term, v) for v in variables)
end
function degree(term::MP.AbstractTermLike, weighted_variables)
    deg = 0
    for (v, w) in weighted_variables
        deg += w * MP.degree(term, v)
    end
    deg
end

"""
    minmaxdegree(f::MP.AbstractPolynomialLike, variables)

Compute the minimum and maximum (weighted total) degree of `f` with respect to the given variables.
"""
function minmaxdegree(f::MP.AbstractPolynomialLike, variables)
    d_min, d_max = typemax(Int), 0
    for term in MP.terms(f)
        d = degree(term, variables)
        d_min, d_max = min(d, d_min), max(d, d_max)
    end
    d_min, d_max
end

"""
    maxdegrees(F, parameters=nothing)

Computes the degrees of the polynomials of `F`.
"""
function maxdegrees(F::MPPolys; parameters = nothing, weights = nothing)
    vars = variables(F, parameters)
    if !isnothing(weights)
        vars = zip(vars, weights)
    end
    last.(minmaxdegree.(F, Ref(vars)))
end
function maxdegrees(C::Composition; parameters = nothing)
    weights = nothing
    for f in reverse(C.polys)
        weights = maxdegrees(f, parameters = parameters, weights = weights)
    end
    weights
end

"""
	multidegrees(F, variable_groups)

Computes the multi-degrees of the polynomials of `F` with respect to the given
variable groups. Returns an integer matrix, each column is the degree with respect
to the corresponding variable group and each row is one of the polynomials in the system.

## Example
```
julia> @polyvar x y;
julia> multidegrees([x*y-2, x^2-4], ([x], [y]))
2×2 Array{Int64,2}:
 1  2
 1  0
 ```
"""
multidegrees(F, vargroups) = [minmaxdegree(f, vars)[2] for vars in vargroups, f in F]
multidegrees(F, VG::VariableGroups) = multidegrees(F, variable_groups(VG))

"""
	support_coefficients(F::MPPolys, vars=variables(F))

Returns a tuple `(A, c)` where `A` is the vector of supports of `F` and `c` is the vector
of corresponding coefficients.
"""
function support_coefficients(F::MPPolys, vars = variables(F))
    terms₁ = MP.terms(F[1])
    c₁ = float.(MP.coefficient.(terms₁))
    A₁ = [Int32(MP.degree(term, var)) for var in vars, term in terms₁]
    A = [A₁]
    c = [c₁]
    for i = 2:length(F)
        termsᵢ = MP.terms(F[i])
        cᵢ = float.(MP.coefficient.(termsᵢ))
        Aᵢ = [Int32(MP.degree(term, var)) for var in vars, term in termsᵢ]
        push!(A, Aᵢ)
        push!(c, cᵢ)
    end
    A, c
end

"""
    compute_numerically_degrees(F::AbstractSystem)

This tries to compute numerically the degrees of `F`. If successfull
this returns a vector containing the degrees. If not `nothign` is returned.
"""
function compute_numerically_degrees(F)
    n, N = size(F)
    # The number of variables match, but it still cannot be homogeneous.
    # We evaluate the system with y:=rand(N) and 2y. If homogeneous then the output
    # scales accordingly to the degrees which we can obtain by taking logarithms.
    x = PVector(rand(ComplexF64, N))
    system_cache = cache(F, x)
    y = evaluate(F, x, system_cache)
    LinearAlgebra.rmul!(x, 2)
    y2 = evaluate(F, x, system_cache)

    degrees = Int[]
    for (y2ᵢ, yᵢ) in zip(y2, y)
        # y2ᵢ = 2^dᵢ yᵢ
        float_dᵢ = log2(abs(y2ᵢ / yᵢ))
        dᵢ = round(Int, float_dᵢ)
        if abs(dᵢ - float_dᵢ) > 1e-10
            return nothing
        end
        push!(degrees, dᵢ)
    end
    degrees
end

"""
    uniquevar(f::MP.AbstractPolynomialLike, tag=:x0)
    uniquevar(F::MPPolys, tag=:x0)

Creates a unique variable.
"""
uniquevar(f::MP.AbstractPolynomialLike, tag = :x0) = MP.similarvariable(f, gensym(tag))
uniquevar(F::MPPolys, tag = :x0) = uniquevar(F[1], tag)
uniquevar(C::Composition, tag = :x0) = uniquevar(C.polys[1][1], tag)

"""
    homogenize(f::MP.AbstractPolynomial, variable=uniquevar(f))

Homogenize the polynomial `f` by using the given variable `variable`.

    homogenize(F::Vector{<:MP.AbstractPolynomial}, variable=uniquevar(F))

Homogenize each polynomial in `F` by using the given variable `variable`.

    homogenize(f::MP.AbstractPolynomial, v::Vector{<:MP.AbstractVariable}, variable=uniquevar(f))

Homogenize the variables `v` in the polynomial `f` by using the given variable `variable`.

    homogenize(F::Vector{<:MP.AbstractPolynomial}, v::Vector{<:MP.AbstractVariable}, variable=uniquevar(F))

Homogenize the variables `v` in each polynomial in `F` by using the given variable `variable`.
"""
function homogenize(f::MP.AbstractPolynomialLike, var = uniquevar(f); parameters = nothing)
    homogenize(f, variables(f, parameters), var)
end
function homogenize(f::MP.AbstractPolynomialLike, vars::Vector, var = uniquevar(f))
    _, d_max = minmaxdegree(f, vars)
    MP.polynomial(map(f) do t
        d = degree(t, vars)
        var^(d_max - d) * t
    end)
end
function homogenize(f::MP.AbstractPolynomialLike, hominfo::HomogenizationInformation)
    g = f
    for (vargroup, homvar) in zip(hominfo.vargroups, hominfo.homvars)
        g = homogenize(g, vargroup, homvar)
    end
    g
end

function homogenize(F::MPPolys, var = uniquevar(F); parameters = nothing)
    homogenize(F, variables(F, parameters), var)
end
function homogenize(F::MPPolys, variables::Vector, var::MP.AbstractVariable = uniquevar(F))
    map(f -> homogenize(f, variables, var), F)
end
function homogenize(F::MPPolys, hominfo::HomogenizationInformation)
    homogenize.(F, Ref(hominfo))
end
function homogenize(
    C::Composition,
    var = uniquevar(C.polys[1]);
    parameters = nothing,
    weights = nothing,
)
    polys = map(length(C.polys):-1:1) do k
        f̄, weights = homogenize_degrees(
            C.polys[k],
            var;
            parameters = parameters, weights = weights,
        )
        if k > 1
            push!(f̄, var)
            push!(weights, 1)
        end
        f̄
    end
    Composition(reverse!(polys))
end

"""
    homogenize_degree(f::MPPoly, variables, var)

Homogenize `f` by using the variable `v`. Returns the homogenized polynomial
and its degree.
"""
function homogenize_degree(f::MPPoly, variables, var::MP.AbstractVariable)
    _, d = minmaxdegree(f, variables)
    MP.polynomial(map(t -> t * var^(d - degree(t, variables)), f)), d
end

"""
    homogenize_degrees(F::MPPolys, var; parameters=nothing, weights=nothing)
    homogenize_degrees(F::MPPolys, variables, var)

Homogenize the system `F` using the variable `var`. Returns the homogenized system
and the degrees of each polynomial.
"""
function homogenize_degrees(F::MPPolys, variables, var::MP.AbstractVariable)
    F̄ = similar(F)
    degrees = Int[]
    for (i, f) in enumerate(F)
        f̄, d = homogenize_degree(f, variables, var)
        F̄[i] = f̄
        push!(degrees, d)
    end

    F̄, degrees
end
function homogenize_degrees(
    F::MPPolys,
    var::MP.AbstractVariable;
    parameters = nothing,
    weights = nothing,
)
    allvars = variables(F, parameters)
    if weights !== nothing
        homogenize_degrees(F, zip(allvars, weights), var)
    else
        homogenize_degrees(F, allvars, var)
    end
end

"""
    remove_zeros!(F::Vector{<:MP.AbstractPolynomialLike})

Remove zero polynomials from the given system `F`.
"""
function remove_zeros!(F::MPPolys)
    filter!(!iszero, F)
    F
end
function remove_zeros!(C::Composition)
    filter!(!iszero, C.polys[1])
    C
end

"""
    has_constant_polynomial(F::Vector{<:MP.AbstractPolynomialLike})

Returns true if the system `F` contains a constant polynomial.
"""
function has_constant_polynomial(F::MPPolys)
    for f in F
        MP.nterms(f) == 1 && MP.isconstant(first(MP.terms(f))) && return true
    end
    false
end
has_constant_polynomial(C::Composition) = has_constant_polynomial(C.polys[1])

permute(F::MPPolys, perm) = F[perm]
permute(C::Composition, perm) = Composition([[C.polys[1][perm]]; C.polys[2:end]])

"""
    check_zero_dimensional(F::Vector{<:MP.AbstractPolynomial})

Check that the given polynomial system can have zero dimensional components.
"""
function check_zero_dimensional(F::Union{MPPolys,Composition})
    N = nvariables(F)
    n = length(F)

    if n ≥ N || (n == N - 1 && is_homogeneous(F))
        return nothing
    end
    error("The input system will not result in a finite number of solutions.")
end

"""
    homogenize_if_necessary(F::MPPolys)

Homogenizes the system `F` if necessary and returns the (new) system `F` its variables
and a subtype of [`AbstractHomogenization`] indicating whether it was homegenized.
If it was homogenized and no then the new variable is the **first one**.
"""
function homogenize_if_necessary(
    F::Union{MPPolys,Composition},
    hominfo::Union{Nothing,HomogenizationInformation};
    vars = nothing,
    parameters = nothing,
)
    if vars === nothing
        vars = variables(F, parameters)
    end

    # This fills in the simple variable group (allvars,)
    hominfo = add_variable_groups(hominfo, F; parameters = parameters)

    if is_homogeneous(F, hominfo; parameters = parameters)
        vargroups = VariableGroups(vars, hominfo)
        F, vargroups, homvars(hominfo)
    elseif hominfo === nothing
        homvar = uniquevar(F)
        push!(vars, homvar)

        F′ = homogenize(F, homvar; parameters = parameters)
        vargroups = VariableGroups(vars, homvar)

        F′, vargroups, (homvar,)
    else
        if hominfo.homvars !== nothing
            error("Input system is not homogeneous although `homvars=$(hominfo.homvars)` was passed.")
        end

        # We create a new variable for each variable group to homogenize the system
        new_homvars = map(_ -> uniquevar(F), hominfo.vargroups)
        push!(vars, new_homvars...)
        new_vargroups = map(hominfo.vargroups, new_homvars) do group, v
            vcat(group, v)
        end
        hominfo = HomogenizationInformation(new_homvars, new_vargroups)

        F′ = homogenize(F, hominfo)
        vargroups = VariableGroups(vars, hominfo)

        F′, vargroups, new_homvars
    end
end

"""
    classify_system(F, vargroups::VariableGroups; affine_tracking=false)

Returns a symbol indicating whether `F` is `:square`, `:overdetermined` or `:underdetermined`.
"""
function classify_system(F, vargroups::VariableGroups; affine_tracking = false)
    n = length(F) - nvariables(vargroups)
    if !affine_tracking
        n += ngroups(vargroups)
    end

    if n == 0
        :square
    elseif n > 0
        :overdetermined
    else
        :underdetermined
    end
end

exponent(term::MP.AbstractTermLike, vars) = [MP.degree(term, v) for v in vars]

"""
    weyldot(f::Polynomial, g::Polynomial)

Compute the [Bombieri-Weyl dot product](https://en.wikipedia.org/wiki/Bombieri_norm).
Note that this is only properly defined if `f` and `g` are homogeneous.

    weyldot(f::Vector{Polynomial}, g::Vector{Polynomial})

Compute the dot product for vectors of polynomials.
"""
function weyldot(
    f::MP.AbstractPolynomialLike{T},
    g::MP.AbstractPolynomialLike{S},
    vars = variables([f, g]),
) where {T,S}
    if f === g
        return sum(f) do term
            abs2(MP.coefficient(term)) / multinomial(exponent(term, vars))
        end
    end
    result = zero(promote_type(T, S))
    for term_f in f
        c_f = MP.coefficient(term_f)
        exp_f = exponent(term_f, vars)
        for term_g in g
            c_g = MP.coefficient(term_g)
            exp_g = exponent(term_g, vars)
            if exp_f == exp_g
                result += (c_f * conj(c_g)) / multinomial(exp_f)
                break
            end
        end
    end
    result
end
function weyldot(F::MPPolys, G::MPPolys, vars = variables(F))
    sum(fg -> weyldot(fg..., vars), zip(F, G))
end

"""
    weylnorm(F, vars=variables(F))

Compute the [Bombieri-Weyl norm](https://en.wikipedia.org/wiki/Bombieri_norm).
Note that this is only properly defined if `f` is homogeneous.
"""
weylnorm(F::MPPolys, vars = variables(F)) = √(sum(f -> weyldot(f, f, vars), F))
weylnorm(f::MPPoly, vars = variables(f)) = √(weyldot(f, f, vars))

"Computes the multinomial coefficient (|k| \\over k)"
function multinomial(k::Vector{Int})
    s = 0
    result = 1
    for i in k
        s += i
        result *= binomial(s, i)
    end
    result
end

"""
	precondition(f::Union{MPPolys,Composition}, variables=variables(f))

Rescale the equations and variables of `f` following the heuristic described in [1].
Returns the scaled system and a vector of scaling factors of the variables.

[1]: HOM4PS-2.0: a software package for solving polynomial systems by the polyhedral homotopy continuation method,
	Lee, T.L., Li, T.Y. & Tsai, C.H. Computing (2008) 83: 109.
	https://doi.org/10.1007/s00607-008-0015-6
"""
function precondition(f::Composition, variables = variables(f))
    vars, equations = preconditioning_factors(expand(f), variables)
    y = vars .* variables
    scale(subs(f, variables => y), equations), vars
end

function precondition(f::MPPolys, variables = variables(f))
    vars, equations = preconditioning_factors(f, variables)
    y = vars .* variables
    map(equations, f) do c, fᵢ
        c * MP.subs(fᵢ, variables => y)
    end, vars
end
function preconditioning_factors(f::MPPolys, vars)
    nvars = length(vars)
    m = sum(nᵢ -> nᵢ + div(nᵢ * (nᵢ - 1), 2), MP.nterms.(f))
    n = nvars + length(f)
    A, b = zeros(Int, m, n), zeros(Float64, m)

    row = 1
    for (i, fᵢ) in enumerate(f)
        termsᵢ = MP.terms(fᵢ)
        ntermsᵢ = length(termsᵢ)
        blockᵢ_start = row - 1
        for t in termsᵢ
            for (j, vⱼ) in enumerate(vars)
                A[row, j] = MP.degree(t, vⱼ)
            end
            A[row, nvars+i] = 1
            b[row] = -log10(convert(Float64, abs(MP.coefficient(t))))
            row += 1
        end

        for k = 1:ntermsᵢ, l = (k+1):ntermsᵢ
            for j = 1:nvars
                A[row, j] = A[blockᵢ_start+k, j] - A[blockᵢ_start+l, j]
            end
            b[row] = b[blockᵢ_start+k] - b[blockᵢ_start+l]
            row += 1
        end
    end
    c = ones(n)
    try
        c = LinearAlgebra.qr(A, Val(true)) \ b
    catch e
        if !isa(e, LinearAlgebra.LAPACKException)
            rethrow(e)
        end
    end
    (variables = exp10.(c[1:nvars]), equations = exp10.(c[nvars+1:end]))
end

"""
	normalize_coefficients(f::Union{MPPolys,Composition})

Rescale the equations of `f` such that each coefficient has absolute value between 0 and 1.
"""
function normalize_coefficients(f::Composition)
    Λ = maximum.(abs, MP.coefficients.(expand(f)))
    scale(f, Λ)
end
function normalize_coefficients(f::MPPolys)
    map(f) do fᵢ
        cmax = maximum(abs, MP.coefficients(fᵢ))
        map(MP.terms(fᵢ)) do t
            c = MP.coefficient(t)
            m = MP.monomial(t)
            c / cmax * m
        end |> MP.polynomial
    end
end

"""
    linear_system(f::Vector{<:MP.AbstractPolynomialLike})

Given a polynomial system which represents a linear system ``Ax=b`` return
`A` and `b`. If `f ` is not a linear system `nothing` is returned.

### Example
```
julia> @polyvar x y z;
julia> f = [2x+3y+z+5, -1x+2y+4z-2];
julia> A, b = linear_system(f);
julia> A
2×3 Array{Int64,2}:
  2  3  1
 -1  2  4

julia> b
2-element Array{Int64,1}:
 -5
  2
```
"""
function linear_system(f::Vector{<:MP.AbstractPolynomialLike}, vars = variables(f))
    n = length(vars)
    A = zeros(MP.coefficienttype(f[1]), length(f), n)
    b = zeros(eltype(A), length(f))
    for (i, fᵢ) in enumerate(f)
        for t in MP.terms(fᵢ)
            constant = true
            for (j, v) in enumerate(vars)
                d = MP.degree(t, v)
                d ≤ 1 || return nothing

                if d == 1
                    A[i, j] = MP.coefficient(t)
                    constant = false
                    break
                end
            end
            if constant
                b[i] = -MP.coefficient(t)
            end
        end
    end
    A, b
end

"""
    monomials(vars::Vector{<:PolyVar}, d; homogeneous::Bool = false)

Create all monomials of a given degree.

```
julia> @polyvar x y
(x, y)

julia> monomials([x,y], 2)
6-element Array{Monomial{true},1}:
 x²
 xy
 y²
 x
 y
 1

julia> monomials([x,y], 2; homogeneous = true)
3-element Array{Monomial{true},1}:
 x²
 xy
 y²
```
"""
function monomials(
    vars::AbstractVector{<:DynamicPolynomials.PolyVar},
    d::Int;
    homogeneous::Bool = false,
)
    n = length(vars)
    if homogeneous
        pred = x -> sum(x) == d
    else
        pred = x -> sum(x) ≤ d
    end
    exps = collect(Iterators.filter(pred, Iterators.product(Iterators.repeated(0:d, n)...)))
    sort!(exps, lt = td_order, rev = true)
    map(exps) do exp
        prod(i -> vars[i]^exp[i], 1:n)
    end
end
function td_order(x, y)
    sx = sum(x)
    sy = sum(y)
    sx == sy ? x < y : sx < sy
end

################
## ModelKit ##
################

polyvar(v::ModelKit.Variable) = DynamicPolynomials.PolyVar{true}(string(v.name))
function maxdegrees(F::ModelKit.System; parameters = F.parameters)
    variables = setdiff!([F.variables; F.parameters], parameters)
    vars = polyvar.(variables)
    params = polyvar.(parameters)
    P = evaluate(F.expressions, variables => vars, parameters => params)
    maxdegrees(P; parameters = params)
end


text:
function assemble(op::Identity,
        testnfs::AbstractTimeBasisFunction,
        trialfns::AbstractTimeBasisFunction)

    tbf = convolve(testnfs, trialfns)
    T = scalartype(tbf)
    z = zeros(T, numintervals(tbf)-1)
    Δt = timestep(tbf)
    for i in eachindex(z)
        p = tbf.polys[i]
        t = (i-1)*Δt
        z[i] = evaluate(p,t)
    end

    return z
end

type TensorOperator <: Operator
    spatial_factor
    temporal_factor
end

⊗(A::AbstractOperator, B::AbstractOperator) = TensorOperator(A, B)
function scalartype(A::TensorOperator)
    promote_type(
        scalartype(A.spatial_factor),
        scalartype(A.temporal_factor))
end

function assemble(operator::TensorOperator, testfns, trialfns)

    test_spatial_basis  = spatialbasis(testfns)
    trial_spatial_basis = spatialbasis(trialfns)

    M = numfunctions(test_spatial_basis)
    N = numfunctions(trial_spatial_basis)

    tbf = convolve(
        temporalbasis(testfns),
        temporalbasis(trialfns))

    kmax = numintervals(tbf) - 1
    k0 = fill(1,    (M,N))
    k1 = fill(kmax, (M,N))

    T = promote_type(scalartype(operator), scalartype(testfns), scalartype(trialfns))
    data = zeros(T, M, N, kmax)
    Z = SparseND.Banded3D(k0, k1, data)

    store(v, m, n, k) = (Z[m,n,k] += v)
    assemble!(operator, testfns, trialfns, store)
    return Z

end

function assemble!(operator::TensorOperator, testfns, trialfns, store)

    space_operator = operator.spatial_factor
    time_operator  = operator.temporal_factor

    space_testfns = spatialbasis(testfns)
    space_trialfns = spatialbasis(trialfns)

    time_testfns = temporalbasis(testfns)
    time_trialfns = temporalbasis(trialfns)

    zt = assemble(time_operator, time_testfns, time_trialfns)
    function store1(v,m,n)
        for (k,w) in enumerate(zt)
            store(w*v,m,n,k)
        end
    end
    assemble!(space_operator, space_testfns, space_trialfns, store1)

end

type TemporalDifferentiation
    operator
end

function assemble!(operator::TemporalDifferentiation, testfns, trialfns, store)

    trial_time_fns  = temporalbasis(trialfns)
    trial_space_fns = spatialbasis(trialfns)

    trialfns = SpaceTimeBasis(
        trial_space_fns,
        derive(trial_time_fns)
    )

    assemble!(operator.operator, testfns, trialfns, store)

end

ext: jl


text:
using Base.Test
using QuantumOptics

mutable struct test_operators <: Operator
  basis_l::Basis
  basis_r::Basis
  data::Matrix{Complex128}
  test_operators(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new(b1, b2, data) : throw(DimensionMismatch())
end

@testset "operators" begin

srand(0)

b1 = GenericBasis(5)
b2 = GenericBasis(3)
b = b1 ⊗ b2
op1 = randoperator(b1)
op = randoperator(b, b)
op_test = test_operators(b, b, op.data)
op_test2 = test_operators(b1, b, randoperator(b1, b).data)
op_test3 = test_operators(b1 ⊗ b2, b2 ⊗ b1, randoperator(b, b).data)
ψ = randstate(b)
ρ = randoperator(b)

@test basis(op1) == b1
@test length(op1) == length(op1.data) == 25

@test_throws ArgumentError op_test*op_test
@test_throws ArgumentError -op_test

@test_throws ArgumentError 1 + op_test
@test_throws ArgumentError op_test + 1
@test_throws ArgumentError 1 - op_test
@test_throws ArgumentError op_test - 1

@test_throws ArgumentError dagger(op_test)
@test_throws ArgumentError identityoperator(test_operators, b, b)
@test_throws ArgumentError trace(op_test)
@test_throws ArgumentError ptrace(op_test, [1])
@test_throws ArgumentError ishermitian(op_test)
@test_throws ArgumentError full(op_test)
@test_throws ArgumentError sparse(op_test)

@test expect(1, op1, ρ) ≈ expect(embed(b, 1, op1), ρ)
@test expect(1, op1, ψ) ≈ expect(embed(b, 1, op1), ψ)
@test expect(op, [ρ, ρ]) == [expect(op, ρ) for i=1:2]
@test expect(1, op1, [ρ, ψ]) == [expect(1, op1, ρ), expect(1, op1, ψ)]

@test variance(1, op1, ρ) ≈ variance(embed(b, 1, op1), ρ)
@test variance(1, op1, ψ) ≈ variance(embed(b, 1, op1), ψ)
@test variance(op, [ρ, ρ]) == [variance(op, ρ) for i=1:2]
@test variance(1, op1, [ρ, ψ]) == [variance(1, op1, ρ), variance(1, op1, ψ)]

@test tensor(op_test) === op_test
@test_throws ArgumentError tensor(op_test, op_test)
@test_throws ArgumentError permutesystems(op_test, [1, 2])

@test embed(b, b, 1, op) == embed(b, 1, op)
@test embed(b, Dict{Vector{Int}, SparseOperator}()) == identityoperator(b)

@test_throws ErrorException QuantumOptics.operators.gemm!()
@test_throws ErrorException QuantumOptics.operators.gemv!()

@test_throws ArgumentError expm(sparse(op1))

@test one(b1).data == diagm(ones(b1.shape[1]))
@test one(op1).data == diagm(ones(b1.shape[1]))

end # testset


text:
#=
1D advection-diffusion using FV
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVaddiff1d");

useLog("FVaddiff1dlog", level=3)

# Configuration setup
domain(1)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 40 # number of elements
mesh(LINEMESH, elsperdim=n, bids=2)

# Variables and BCs
u = variable("u", location=CELL)
v = variable("v", location=CELL)
boundary(u, 1, FLUX, 0)
boundary(u, 2, NO_BC)
boundary(v, 1, FLUX, 0)
boundary(v, 2, NO_BC)

# Time interval and initial condition
T = 0.5;
timeInterval(T)
initial(u, "x<0.1 ? 1 : 0")
initial(v, "x<0.1 ? 1 : 0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F.n ds)
coefficient("a", 1) # advection velocity
# The "upwind" function applies upwinding to the term (a.n)*u with flow velocity a.
# The optional third parameter is for tuning. Default upwind = 0, central = 1. Choose something between these.
coefficient("D", 0.01) # Diffusion rate
flux([u, v], ["upwind(a,u) - D*dot(grad(u),normal())", "upwind(a,v)"])
# Note that there is no source() for this problem.

#@exportCode("fvad1dcode") # uncomment to export generated code to a file
#@importCode("fvad1dcode") # uncomment to import code from a file

solve([u,v])

finalize_finch()

##### Uncomment below to plot

# # The initial condition
# init = zeros(n);
# x = Finch.fv_info.cellCenters[:]
# for i=1:n
#     if x[i] < 0.1
#         init[i] = 1
#     end
# end

# using Plots
# pyplot();
# display(plot([x x x], [init u.values[:] v.values[:]], markershape=:circle, label=["initial" "upwind advection, central diffusion" "advection only"]))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
#=
maximize likelihood of d observations at discrete points appearing with random
frequencies, subject to probability vector not being too far from a uniform
prior
=#

struct DiscreteMaxLikelihood{T <: Real} <: ExampleInstanceJuMP{T}
    d::Int
    use_EF::Bool
end

function build(inst::DiscreteMaxLikelihood{T}) where {T <: Float64}
    d = inst.d
    @assert d >= 2
    freq = rand(1:(2 * d), d)
    freq /= sum(freq)

    model = JuMP.Model()
    JuMP.@variable(model, p[1:d])
    JuMP.@variable(model, hypo)
    JuMP.@objective(model, Max, hypo)
    JuMP.@constraint(model, sum(p) == 1)

    JuMP.@constraint(model, vcat(hypo, p) in Hypatia.HypoPowerMeanCone{T}(freq))

    ext = (inst.use_EF ? VecNegEntropyEF() : VecNegEntropy())
    add_spectral(ext, d, vcat(inv(d), inv(d), p), model)

    # save for use in tests
    model.ext[:freq] = freq
    model.ext[:p_var] = p

    return model
end

function test_extra(inst::DiscreteMaxLikelihood{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective and constraints
    tol = eps(T)^0.2
    freq = model.ext[:freq]
    p_opt = JuMP.value.(model.ext[:p_var])
    @test sum(p_opt) ≈ 1 atol=tol rtol=tol
    @test minimum(p_opt) >= -tol
    p_opt = pos_only(p_opt)
    obj_result = exp(sum(f_i * log(p_i) for (f_i, p_i) in zip(freq, p_opt)))
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    return
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VerificationIPFlowResult <: SwaggerModel
    access::Any # spec type: Union{ Nothing, Access } # spec name: access
    ruleName::Any # spec type: Union{ Nothing, String } # spec name: ruleName

    function VerificationIPFlowResult(;access=nothing, ruleName=nothing)
        o = new()
        validate_property(VerificationIPFlowResult, Symbol("access"), access)
        setfield!(o, Symbol("access"), access)
        validate_property(VerificationIPFlowResult, Symbol("ruleName"), ruleName)
        setfield!(o, Symbol("ruleName"), ruleName)
        o
    end
end # type VerificationIPFlowResult

const _property_map_VerificationIPFlowResult = Dict{Symbol,Symbol}(Symbol("access")=>Symbol("access"), Symbol("ruleName")=>Symbol("ruleName"))
const _property_types_VerificationIPFlowResult = Dict{Symbol,String}(Symbol("access")=>"Access", Symbol("ruleName")=>"String")
Base.propertynames(::Type{ VerificationIPFlowResult }) = collect(keys(_property_map_VerificationIPFlowResult))
Swagger.property_type(::Type{ VerificationIPFlowResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VerificationIPFlowResult[name]))}
Swagger.field_name(::Type{ VerificationIPFlowResult }, property_name::Symbol) =  _property_map_VerificationIPFlowResult[property_name]

function check_required(o::VerificationIPFlowResult)
    true
end

function validate_property(::Type{ VerificationIPFlowResult }, name::Symbol, val)
end


text:
chr2ind(s::String, n::Integer) = Base.chr2ind(s, n)
function chr2ind(s::String, range::Range)
    if length(s) == 0 || length(range) == 0
        return 0:-1
    else
        ec = chr2ind(s, range.stop)
        sc = chr2ind(s, range.start)
        if step(range)>0
            sc:step(range):(ec+sizeof(s[ec:ec])-1)
        else
            sc+sizeof(s[sc:sc])-1:step(range):ec
        end
    end
end
"""
unirange(s::String, range::Range)

Like s[range] but for unicode.
"""
function unirange(s::String, range::Range)
    if length(s) == 0 || length(range) == 0
        ""
    elseif step(range) == 1
        r = chr2ind(s, range)
        s[r.start:r.stop]
    elseif step(range) == -1
        s = reverse(s)
        r = chr2ind(s, length(s)-range+1)
        s[r.start:r.stop]
    else
        error("Only single step range permitted")
    end
end

"""
replace!(c::Array{String, 1}, pattern, repl)

Replace pattern with repl in each member of c.
"""
function replace!(c::Array{String, 1}, pattern, repl)
    map!(s->replace(s, pattern, repl), c)
end


text:
begin
#Linear Search

# Julia code to linearly search x in arr[].
# If x is present then return its location,
# otherwise return -1
import Base:length,show,AbstractArray
arr = [2, 3, 4, 10, 40]
x = 10
n = length(arr)

function linearsearch(arr::AbstractArray{},n,x)
	for i in 1:n
		if arr[i] == x
		 	return i
		else
				continue
		end
	end
	return -1
end
global result=-1
# Function call
# for i in 1:n
# 	if arr[i] == x
# 		global result=i
# 		break
# 	else
# 		continue
# 	end
# end
result=linearsearch(arr::AbstractArray{},n,x)

println("Result after function execution:",result)
if result == -1
    println("Element:",x," is not present in array",arr)
else
    println("Element:",x," is present at index:", result," in array",arr)
end
end

max_stars_repo_name: Learning-Julia/Algorithms

max_issues_repo_name: Learning-Julia/Algorithms

max_forks_repo_name: Learning-Julia/Algorithms


text:
#=
    NEMO: Next Energy Modeling system for Optimization.
    https://github.com/sei-international/NemoMod.jl

    Copyright © 2018: Stockholm Environment Institute U.S.

    File description: NemoMod module.
=#

module NemoMod
export calculatescenario, createnemodb, dropdefaultviews, dropresulttables, logmsg, setparamdefault

#= List of module global variables.
    • csdbpath - dbpath argument in last invocation of calculatescenario()
    • csquiet - quiet argument in last invocation of calculatescenario()
    • csjumpmodel - jumpmodel argument in last invocation of calculatescenario() (only set if a customconstraints include is performed)
=#

# BEGIN: Access other modules and code files.
using JuMP, SQLite, DataFrames, Distributed, Dates, ConfParser, MathOptInterface
using GLPK, Cbc  # Open-source solvers

# Proprietary solvers - enclosed in try blocks for users who aren't using NEMO installer
try
    using CPLEX
catch
    # Just continue
end

try
    using Gurobi
catch
    # Just continue
end

try
    using Mosek
catch
    # Just continue
end

# Do not load Xpress since it throws an unrecoverable error when loaded in a package on a machine without an Xpress license

include("db_structure.jl")  # Functions for manipulating structure of scenario databases
include("other_functions.jl")  # Core NEMO functions
include("scenario_calculation.jl")  # Functions for calculating a scenario with NEMO
# END: Access other modules and code files.

end  # module NemoMod


text:
using CSV

function ccf_orders(order_list_timeseres::AbstractChunkListTimeseries,  line_list_df::DataFrame, pipeline::PipelinePlan; verbose::Bool = false, recalc::Bool = false,
   range_no_mask_change::Real=RvSpectMLBase.max_bc, ccf_mid_velocity::Real=0.0, mask_scale_factor::Real=1, mask_type::Symbol = :tophat )

   if need_to(pipeline, :ccf_orders)  # Compute order CCF's & measure RVs
      if mask_type == :tophat
        mask_shape = CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor*mask_scale_factor)
      else
        @error("Requested mask shape (" * string(mask_type) * " not avaliable.")
      end

      line_list = CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
      ccf_plan = CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=range_no_mask_change)
      v_grid = CCF.calc_ccf_v_grid(ccf_plan)
      tstart = now()    # Compute CCFs for each order
      order_ccfs = CCF.calc_order_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
      if verbose   println("# Order CCFs runtime: ", now()-tstart)  end

      set_cache!(pipeline,:ccf_orders, (order_ccfs = order_ccfs, v_grid=v_grid) )
      if save_data(pipeline, :ccf_orders)
         for (i, order) in orders_to_use
            if !(sum(order_ccfs[:,i,:]) > 0)   continue    end
            local t = Tables.table( order_ccfs[:,i,:]', header=Symbol.(v_grid) )
            CSV.write(joinpath(output_dir,target_subdir * "_ccf_order=" * string(order) * ".csv"),t)
         end
      end
      dont_need_to!(pipeline, :ccf_orders);
   end
   if has_cache(pipeline,:ccf_orders) return read_cache(pipeline,:ccf_orders)
   else   @error("Invalid pipeline state.")          end
end


text:
using DifferentialEquations
using Plots
plotlyjs()

# constants
l = 1.0                             # length [m]
m = 1.0                             # mass[Kg]
g = 9.81                            # gravitational acceleration [m/s²]
b = 0.05                           # damping coefficient [s^-1]

#function to compute the ODEs
function pendulum!(du,u,p,t)
    du[1] = u[2]                    # θ'(t) = ω(t)
    du[2] = -(b*u[2])-(g\l)*sin(u[1]) # ω'(t) = -(g/l) sin θ(t) + 3/(ml^2)M(t)
end

#initial conditions
θ₀ = π/6                         # initial angular deflection [rad]
ω₀ = 0                            # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀]                       # initial state vector
tspan = (0.0,200.0)                  # time interval

#defining the problem and solving
prob = ODEProblem(pendulum!,u₀,tspan)
sol = solve(prob,saveat=0.1)

#plotting the solution
display(plot(sol,linewidth=2,
    xaxis="t (s)",color=["green" "blue"],
    label=["θ [rad]" "ω [rad/s]"],layout=(2,1)))
savefig("./plots/simple damped pendulum-plot.png")

#trajectory
t=sol.t
U=sol[1:end,:]
(θ,ω)=[U[x,:] for x in 1:size(U,1)]
x=l.*sin.(π.-θ)
y=zeros(size(t))
z=l.*cos.(π.-θ)
display(plot(x,y,z,
    linewidth=2,label="trajectory",xaxis="x (m)",
    yaxis="y (m)",zaxis="z (m)",grid=(:on,:black)))
savefig("./plots/simple damped pendulum-trajectory.png")

max_stars_repo_path: Modeling the Dragon Module as a Pendulum/Scripts/pendulum damped.jl


text:
#=
  # Follows Ziller et al 2013.
  # There is a straight port of get_beta_parameters from the R code
  # This is called get_beta_parametersPort
=#

using Distributions
using DataFrames
using Docile
using Lumberjack

type CpGMethEstimate
    alpha::Float64
    beta::Float64
    methcounts::Array{Int64}
    unmethcounts::Array{Int64}
end

"""
 CpGMethCount

 Records the methylation counts in each replicate for a CpG position

 methcounts:  array of methylated counts in replicates at CpG site
 coverages:   array of coverage in replicates at CpG site
"""
type CpGMethCount
    methcounts::Array{Int64}
    coverages::Array{Int64}
end

function issame(array)
  for i=1:length(array)
     if array[i] != array[1]
       return false
     end
  end
  return true
end

"""
 computeBDI
 monte carlo approximation of beta difference distribution

 #### arguments
  a1: alpha of condition1
  b1: beta of condition1
  a2: alpha of condition2
  b2: beta of condition2
  nSamp: Number of samples to use
  alpha: Report credible interval 1-alpha

 #### returns
  array with values corresponding in order to the following:
   s1: one sided p value
   qt1: 1-α credibility interval left boundary
   qt2: 1-α credibiliy interval right boundary
   μ : mean differences
   v : stdev difference
"""

function computeBDI(a1,b1,a2,b2,nSamp=10000,alpha=0.1)
    if (alpha < 0) | (alpha > 1)
        throw(ArgumentError("alpha out of range"))
    end
    #monte carlo approximation
    if(a1 <= 0 ) | ( b1 <= 0) | (a2 <= 0) | (b2 <= 0)
       throw(ArgumentError("a or b out of range ($a1,$b1, $a2, $b2)"))
    end

    beta1  = 0.0
    beta2  = 0.0

    try
        beta1  = Beta(a1,b1)
        beta2  = Beta(a2,b2)
    catch e
        println("a or b out of range ($a1,$b1, $a2, $b2)")
        error(e)
    end

    theta1 = rand(beta1,nSamp)
    theta2 = rand(beta2,nSamp)
    #calculate differences between two distributions
    di=theta1-theta2
    qt=quantile(di,[alpha/2,1-alpha/2])
    # s1 one sided p value:
    # i.e. if they were the same distribution, we would have ~50% > or <
    #         min(mean( [1,0,1,0,1,1,1,0,01,0,1,0,1]), mean( [1,0,1,0,1,1,1,0,01,0,1,0,1]))
    #      if different, then we will bias towards one distribution being different something like
    #         min(mean([1,0,0,0,0,1,0,0,0,1,0,0,1,0,0] ) mean([0,1,1,1,1,0,1,1,1,0,1,1,0,1,1]))
    s1=min( mean(theta1 .>=theta2), mean(theta1 .<=theta2))
    mu=mean(di)
    v=std(di)
    # note that qt gets splatted in - it has two values, so returning size 5 array
    qs= [ s1, qt, mu, v ]
    return(qs)
end

"""
 getBDI

"""
function getBDI(alpha1,beta1,navg1,alpha2,beta2,navg2;alpha=0.1,nSamp=10000)
  na_found = all( x-> typeof(x) == nothing, [alpha1,beta1,navg1,alpha2,beta2,navg2])
  nan_found = any( x-> Base.isnan(x), [alpha1,beta1,navg1,alpha2,beta2,navg2])
  if na_found | nan_found
      #return( [nan, nan, nan, nan, nan] )
      return([NaN,NaN,NaN,NaN,NaN])
  else
      a1=(alpha1+1)/(alpha1-1)-1
      b1=(beta1+1)/(beta1-1)-1
      if (navg1 < 30) & (navg2 < 30)
          x1=alpha1+1
          x2=navg1+1
          # monte carlo approximation
          #println("monte carlo estimated BDI")
          return  computeBDI(alpha1,beta1,alpha2,beta2,nSamp,alpha )
      else
           #use standard normal to approximate BDI for n>=30
           #println("Normal estimated BDI")
           norm=Normal()
           mu=alpha1/(alpha1+beta1) - alpha2 / (alpha2+beta2)
           sigma= alpha1*beta1 / ( (alpha1+beta1)^2 * (alpha1+beta1+1) )
           sigma=sqrt( sigma+alpha2*beta2 / ( (alpha2+beta2)^2 * (alpha2+beta2 +1) ) )
           np=quantile(norm, (1 - alpha/2) )
           p=min( cdf(norm, mu/sigma ) ,  1- cdf( norm, mu/sigma )  ) #as.numeric(mu/sigma>=np*(1-alpha/2))-as.numeric(mu/sigma<=(-1*(np*(1-alpha/2))))
           return [p, mu-(np*sigma), mu+np*sigma, mu, sigma]
      end
  end
end

function getDCpGStatistic(cpg_meth_est1::CpGMethEstimate,cpg_meth_est2::CpGMethEstimate;alpha=0.1,nSamp=10000)
    meth_level1 = cpg_meth_est1.alpha /  ( cpg_meth_est1.alpha + cpg_meth_est1.beta )
    meth_level2 = cpg_meth_est2.alpha /  ( cpg_meth_est2.alpha + cpg_meth_est2.beta )
    ave_cov1 = ( sum( cpg_meth_est1.methcounts) + sum( cpg_meth_est1.unmethcounts) ) /2
    ave_cov2 = ( sum( cpg_meth_est2.methcounts) + sum( cpg_meth_est2.unmethcounts) ) /2
    if meth_level1 <= meth_level2
        difference = getBDI(cpg_meth_est1.alpha,cpg_meth_est1.beta,ave_cov1,
                                     cpg_meth_est2.alpha,cpg_meth_est2.beta,ave_cov2,alpha,nSamp)
    else
        difference = getBDI(cpg_meth_est2.alpha,cpg_meth_est2.beta,ave_cov2,
                                     cpg_meth_est1.alpha,cpg_meth_est1.beta,ave_cov1,alpha,nSamp)
    end
    return difference
end

function _getDCpGStatistic(cpg_meth_est1::Array{CpGMethEstimate},cpg_meth_est2::Array{CpGMethEstimate}, sample_names)
    num_cpgs = length(cpg_meth_est1)
    differences = Any[]
    for i=1:num_cpgs
        meth_level1 = cpg_meth_est1[i].alpha /  ( cpg_meth_est1[i].alpha + cpg_meth_est1[i].beta )
        meth_level2 = cpg_meth_est2[i].alpha /  ( cpg_meth_est2[i].alpha + cpg_meth_est2[i].beta )
        ave_cov1 = ( sum( cpg_meth_est1[i].methcounts) + sum( cpg_meth_est1[i].unmethcounts) ) /2
        ave_cov2 = ( sum( cpg_meth_est2[i].methcounts) + sum( cpg_meth_est2[i].unmethcounts) ) /2
        if meth_level1 <= meth_level2
           push!(differences, getBDI(cpg_meth_est1[i].alpha,cpg_meth_est1[i].beta,ave_cov1,
                                     cpg_meth_est2[i].alpha,cpg_meth_est2[i].beta,ave_cov2))
        else
           push!(differences, getBDI(cpg_meth_est2[i].alpha,cpg_meth_est2[i].beta,ave_cov2,
                                     cpg_meth_est1[i].alpha,cpg_meth_est1[i].beta,ave_cov1) )
        end
    end
    return differences
end

"""
  getDCpGStatistic

  Get Statistics on Differential Methylation at a CpG.

  #### Arguments
      * cpg_meth_est1: Array of CpGMethEstimates condition 1
      * cpg_meth_est2: Array of CpGMethEstimates condition 2
      * sample_names:  Currently unused

  #### Returns
      * DataFrame with columns:
           significance              Float64[],
           credibilityIntervalMin    Float64[],
           credibilityIntervalMax    Float64[],
           estimatedDifferenceMean   Float64[],
           estimatedVariance         Float64[]

"""
function getDCpGStatistic(cpg_meth_est1::Array{CpGMethEstimate},cpg_meth_est2::Array{CpGMethEstimate}, sample_names;verbose=false)
    num_cpgs = length(cpg_meth_est1)
    differences = Any[]
    for i=1:num_cpgs
        meth_level1 = cpg_meth_est1[i].alpha /  ( cpg_meth_est1[i].alpha + cpg_meth_est1[i].beta )
        meth_level2 = cpg_meth_est2[i].alpha /  ( cpg_meth_est2[i].alpha + cpg_meth_est2[i].beta )
        ave_cov1 = ( sum( cpg_meth_est1[i].methcounts) + sum( cpg_meth_est1[i].unmethcounts) ) /2
        ave_cov2 = ( sum( cpg_meth_est2[i].methcounts) + sum( cpg_meth_est2[i].unmethcounts) ) /2
        if meth_level1 <= meth_level2
           push!(differences, getBDI(cpg_meth_est1[i].alpha,cpg_meth_est1[i].beta,ave_cov1,
                                     cpg_meth_est2[i].alpha,cpg_meth_est2[i].beta,ave_cov2))
        else
           push!(differences, getBDI(cpg_meth_est2[i].alpha,cpg_meth_est2[i].beta,ave_cov2,
                                     cpg_meth_est1[i].alpha,cpg_meth_est1[i].beta,ave_cov1) )
        end

        if verbose
            if i % 100 == 0
                Lumberjack.info("processed $i CpGs")
            end
        end

    end
    outdf = DataFrame(
           significance              =Float64[],
           credibilityIntervalMin    =Float64[],
           credibilityIntervalMax    =Float64[],
           estimatedDifferenceMean   =Float64[],
           estimatedVariance         =Float64[]
    )
    for row in differences
        push!(outdf, row)
    end
    return outdf
end

"""
  getDCpGStatistics

  Arguments:
   * cond1_meth_counts: array of CpGMethCount objects
   * cond2_meth_counts: array of CpGMethCount objects
  Returns
    DataFrame of statistics on differences of CpGs with Columns:
    * significance
    * credibilityIntervalMin
    * credibilityIntervalMax
    * estimatedDifferenceMean
    * estimatedVariance
"""
function getDCpGStatistic( cond1_meth_counts::Array{CpGMethCount}, cond2_meth_counts::Array{CpGMethCount})
    meth_ests_cond1 = Any[]
    meth_ests_cond2 = Any[]
    for i=1:length(cond1_meth_counts)
         cpg_meth_est_cond1=get_beta_parameters(cond1_meth_counts[i].methcounts, cond1_meth_counts[i].coverages)
         cpg_meth_est_cond2=get_beta_parameters(cond1_meth_counts[i].methcounts, cond1_meth_counts[i].coverages)
         push!(meth_ests_cond1, cpg_meth_est_cond1)
         push!(meth_ests_cond2, cpg_meth_est_cond2)
    end
    diff_cpg_stats = getDCpGStatistic(meth_ests_cond1,meth_ests_cond2,["A","B"])
    return diff_cpg_stats
end

#expects vector size n=4, meth1,total1,meth2,total2
function get_beta_parametersPort(dat)
    p=[dat[1]/dat[2],dat[3]/dat[4]]
    if(p[1] == p[2])  #case one, use uninformative prior
        (max,i)=findmax(dat[ [2,4] ])
        m=dat[[1,3]][i]
        n=dat[[2,4]][i]
        return [m+1,n-m+1]
    else
        mu=(dat[1]+dat[3])/(dat[2]+dat[4])
        #weighted variance
        w=[dat[2]/(dat[2]+dat[4]),dat[4]/(dat[2]+dat[4])]#basic weights
        s2=sum(w[1]*(dat[1]/dat[2]-mu)^2+w[2]*(dat[3]/dat[4]-mu)^2)
        nm=(dat[2]+dat[4])/2#mean sample coverage
        #estimate systematic variance between replicates
        V=(s2-1/nm*(mu*(1-mu)))/(1-1/nm)
        #M=(mu*(1-mu)-s2)/(s2-1/nm*(mu*(1-mu)))
        th=(1/maximum(dat[[2,4]]))^2
    end

    #println("V = $V s2 = $s2")

    if (V<th)
        M=0
    else #get M parameter
        M=mu*(1-mu)/V -1
    end
    params=[ mu*M, M*(1-mu)]
    #special treatment for overdispersion case a,b negative
    if ( (params[1] <=0 ) | (params[2] <= 0 )  )
        (max,i)=findmax(dat[ [2,4] ] )
        return [ dat[ [1,3]  ][i]+1, dat[  [2,4] ][i]-dat[  [1,3] ][i]+1 ]
    else
        return(params)
    end
end

"""
 Estimate Alpha and Beta parameters of Beta Distribution given replicates

 ### Arguments:
     * `methcounts`:  Array  of replicate counts corresponding to methylated reads at CpG/Cytosine Position
     * `totalcounts`: Array  of replicate counts corresponding to total coverage at CpG/Cytosine Position

 ### Returns
    * Floating Point Array [alpha,beta]

"""

function _get_beta_parameters(methcounts,totalcounts)
    num_reps = length(methcounts)

    # if num_reps == 1
    # then return...
    if num_reps==1 #one replicate case
        return [methcounts[1]+1, totalcounts[1]-methcounts[1]+1 ]
    end
    p=fill(0.0,num_reps)
    for r=1:num_reps
        p[r] = methcounts[r] / totalcounts[r]
    end

    if issame(p)
         # case one - use uninformative prior
         # where the probability estimates are
         # the same for each replicate
        (max_cov,max_cov_idx)=findmax(totalcounts)
         m = methcounts[max_cov_idx]
         n = totalcounts[max_cov_idx]
         return [m+1,n-m+1]
    else
    # case two
         μ = sum(methcounts) / sum(totalcounts)

         # weighted variance
         s2 = 0.0
         for r=1:num_reps
             w  = totalcounts[r] / sum(totalcounts) #basic weights
             s2 += (  w * (p[r] - μ)^2  )
         end
         # mean sample coverage
         nm= sum(totalcounts) / num_reps
         # estimate systematic variance between replicates
         V=( s2 - 1/nm * ( μ * (1 - μ )  ) ) /  ( 1 - 1/nm )
         θ = ( 1 / maximum( totalcounts  ) )^2
    end

    #println("V = $V s2 = $s2")

    if(V < θ )
      M = 0
    else
      M = μ * (1-μ)/V -1
    end

    alpha = μ * M
    beta  = M*(1-μ )

    if ( (alpha <= 0) | (beta <= 0 ) )
        (max_cov,max_cov_idx)=findmax( totalcounts )
        return [ methcounts[max_cov_idx]+1, totalcounts[max_cov_idx] - methcounts[max_cov_idx] + 1 ]
    else
        return( [alpha,beta]  )
    end
end

"""
 Estimate Alpha and Beta parameters of Beta Distribution given replicates

 ### Arguments:
     * `methcounts`:  Array  of replicate counts corresponding to methylated reads at CpG/Cytosine Position
     * `totalcounts`: Array  of replicate counts corresponding to total coverage at CpG/Cytosine Position

 ### Returns
    * CpGMethEstimate Object

"""

function get_beta_parameters(methcounts,totalcounts)
     beta_params  = _get_beta_parameters(methcounts,totalcounts)
     unmethcounts = totalcounts - methcounts
     return CpGMethEstimate( beta_params[1], beta_params[2], methcounts, unmethcounts)
end

# this differs only in the way we calculate M, which follows the formula in the
# paper - I have not done the algebra so this is just to show that we get equivalent
# values
function get_beta_parameters3(methcounts,totalcounts)
    num_reps = length(methcounts)
    s2 = 0.0
    # if num_reps == 1
    # then return...

    p=fill(0.0,num_reps)
    for r=1:num_reps
        p[r] = methcounts[r] / totalcounts[r]
    end

    if isequal(p,p)
    # case one - use uninformative prior
         # where the probability estimates are
         # the same for each replicate
        (max_cov,max_cov_idx)=findmax(totalcounts)
         m = methcounts[max_cov_idx]
         n = totalcounts[max_cov_idx]
         return [m+1,n-m+1]
    else
    # case two
         μ = sum(methcounts) / sum(totalcounts)

         # weighted variance
         # in the text we have this expressed as:
         #    sum( n_ij ) * (e_ij - mu_i)^2 / sum n_ij
         # which could be expressed
         #     w * (e_ij - mu_i)^2
         #   where w = n_ij / sum(n_ij)
         # So we can first just calculate the "weights"

        for r=1:num_reps
             w  = totalcounts[r] / sum(totalcounts) #basic weights
             s2  += (  w * (p[r] - μ)^2  )
         end

         # mean sample coverage
         nm= sum(totalcounts) / num_reps
         # estimate systematic variance between replicates
         V=( s2 - 1/nm * ( μ * (1 - μ )  ) ) /  ( 1 - 1/nm )
         θ = ( 1 / maximum( totalcounts  ) )^2
    end

    if(V < θ )
      M = 0
    else
      #M = μ * (1-μ)/V
      M = ( μ*(1-μ) - s^2 ) / ( s^2 - 1/nm * (μ)(1-μ) )
    end

    alpha = μ * M
    beta  = M*(1-μ )
    if ( alpha <= 0 | beta <= 0  )
        max_cov_idx=findmax( totalcounts )
        return [ methcounts[max_cov_idx]+1, totalcounts[max_cov_idx] - methcounts[max_cov_idx] + 1 ]
    else
        return( [alpha,beta]  )
    end
end


text:
import TimeseriesSurrogates.randomshuffle
import TimeseriesSurrogates.randomphases
import TimeseriesSurrogates.randomamplitudes
import TimeseriesSurrogates.aaft
import TimeseriesSurrogates.iaaft
import DynamicalSystems.Dataset
import StateSpaceReconstruction: Embeddings

############################
# Multivariate time series
############################

"""
    randomshuffle(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))

Column-wise random shuffle surrogate of an array,
where each column is a scalar-valued time series. `cols` controls which
variables of the embedding are shuffled.
"""
function randomshuffle(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = randomshuffle(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    d_shuffled
end

"""
    randomphases(a::AbstractArray{Number, 2}, cols = 1:size(d, 2))

Column-wise random phases Fourier surrogate of an array,
where each column is a scalar-valued time series. `cols` controls
which variables of the embedding are shuffled.
"""
function randomphases(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = randomphases(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    d_shuffled
end

"""
    randomamplitudes(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))

Column-wise random amplitude Fourier surrogate of an array,
where each column is a scalar-valued time series. `cols` controls
which variables of the embedding are shuffled.
"""
function randomamplitudes(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = randomamplitudes(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    d_shuffled
end

"""
    aaft(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))

Column-wise amplitude-adjusted Fourier transform (AAFT) surrogate of an array,
where each column is a scalar-valued time series.
`cols` controls which variables of the embedding are shuffled.
"""
function aaft(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = aaft(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    d_shuffled
end

"""
    iaaft(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))

Iterated amplitude-adjusted Fourier transform (IAAFT) surrogate of an array,
where each column is a scalar-valued time series.
`cols` controls which variables of the embedding are shuffled.
"""
function iaaft(a::AbstractArray{Number, 2}; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = iaaft(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    Dataset(d_shuffled)
end

##############
# Embeddings
##############
"""
    randomshuffle(E::Embeddings.AbstractEmbedding;
                    cols = 1:size(E.points, 1))

Column-wise random shuffle surrogate of an embedding. `cols` controls which
variables of the embedding are shuffled.
"""
function randomshuffle(E::Embeddings.AbstractEmbedding;
                        cols = 1:size(E.points, 1))
    n_variables = size(E.points, 1)
    E_shuffled = similar(E.points)
    for i = 1:n_variables
        if i ∈ cols
            E_shuffled[i, :] = randomshuffle(E.points[i, :])
        else
            E_shuffled[i, :] = E.points[i, :]
        end
    end

    embed(transpose(E_shuffled))
end

"""
    randomphases(E::Embeddings.AbstractEmbedding; cols = 1:size(E.points, 1))

Column-wise random phases Fourier surrogate of an embedding.
`cols` controls which variables of the embedding are shuffled.
"""
function randomphases(E::Embeddings.AbstractEmbedding;
                        cols = 1:size(E.points, 1))
    n_variables = size(E.points, 1)
    E_shuffled = similar(E.points)
    for i = 1:n_variables
        if i ∈ cols
            E_shuffled[i, :] = randomphases(E.points[i, :])
        else
            E_shuffled[i, :] = E.points[i, :]
        end
    end

    embed(transpose(E_shuffled))
end

"""
    randomamplitudes(E::Embeddings.AbstractEmbedding;
                        cols = 1:size(E.points, 1))

Column-wise random amplitude Fourier surrogate of an embedding. `cols`
controls which variables of the embedding are shuffled.
"""
function randomamplitudes(E::Embeddings.AbstractEmbedding;
                            cols = 1:size(E.points, 1))
    n_variables = size(E.points, 1)
    E_shuffled = similar(E.points)
    for i = 1:n_variables
        if i ∈ cols
            E_shuffled[i, :] = randomamplitudes(E.points[i, :])
        else
            E_shuffled[i, :] = E.points[i, :]
        end
    end

    embed(transpose(E_shuffled))
end

"""
    aaft(E::Embeddings.AbstractEmbedding; cols = 1:size(E.points, 1))

Column-wise amplitude-adjusted Fourier transform (AAFT) surrogate of an
embedding. `cols` controls which variables of the embedding are shuffled.
"""
function aaft(E::Embeddings.AbstractEmbedding; cols = 1:size(E.points, 1))
    n_variables = size(E.points, 1)
    E_shuffled = similar(E.points)
    for i = 1:n_variables
        if i ∈ cols
            E_shuffled[i, :] = aaft(E.points[i, :])
        else
            E_shuffled[i, :] = E.points[i, :]
        end
    end

    embed(transpose(E_shuffled))
end

"""
    iaaft(E::Embeddings.AbstractEmbedding; cols = 1:size(E.points, 1))

Column-wise iterated amplitude-adjusted Fourier transform (IAAFT) surrogate
of an embedding. `cols` controls which variables of the embedding are shuffled.
"""
function iaaft(E::Embeddings.AbstractEmbedding; cols = 1:size(E.points, 1),
				n_maxiter = 200, tol = 1e-6, n_windows = 50)
    n_variables = size(E.points, 1)
    E_shuffled = similar(E.points)
    for i = 1:n_variables
        if i ∈ cols
            E_shuffled[i, :] = iaaft(E.points[i, :],
								n_maxiter = n_maxiter,
								tol = tol,
								n_windows = n_windows)
        else
            E_shuffled[i, :] = E.points[i, :]
        end
    end

    embed(transpose(E_shuffled))
end

##############
# Datasets
##############

"""
    randomshuffle(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))

Column-wise random shuffle surrogate of a Dataset. `cols` controls which
variables of the embedding are shuffled.
"""
function randomshuffle(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = randomshuffle(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    Dataset(d_shuffled)
end

"""
    randomphases(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))

Column-wise random phases Fourier surrogate of a Dataset. `cols` controls
which variables of the embedding are shuffled.
"""
function randomphases(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = randomphases(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    Dataset(d_shuffled)
end

"""
    randomamplitudes(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))

Column-wise random amplitude Fourier surrogate of a Dataset. `cols` controls
which variables of the embedding are shuffled.
"""
function randomamplitudes(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = randomamplitudes(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    Dataset(d_shuffled)
end

"""
    aaft(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))

Column-wise amplitude-adjusted Fourier transform (AAFT) surrogate of a Dataset.
`cols` controls which variables of the embedding are shuffled.
"""
function aaft(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = aaft(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    Dataset(d_shuffled)
end

"""
    iaaft(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))

Iterated amplitude-adjusted Fourier transform (IAAFT) surrogate of a Dataset.
`cols` controls which variables of the embedding are shuffled.
"""
function iaaft(d::DynamicalSystemsBase.Dataset; cols = 1:size(d, 2))
    n_variables = size(d, 2)
    d_shuffled = zeros(eltype(d), size(d))
    for i = 1:n_variables
        if i ∈ cols
            d_shuffled[:, i] = iaaft(d[:, i])
        else
            d_shuffled[:, i] = d[:, i]
        end
    end

    Dataset(d_shuffled)
end

##########################
# CustomReconstructions
##########################

randomshuffle(d::CustomReconstruction; cols = 1:size(d, 2)) =
    randomshuffle(d.reconstructed_pts; cols = cols)

randomphases(d::CustomReconstruction; cols = 1:size(d, 2)) =
    randomphases(d.reconstructed_pts; cols = cols)

randomamplitudes(d::CustomReconstruction; cols = 1:size(d, 2)) =
    randomamplitudes(d.reconstructed_pts; cols = cols)

aaft(d::CustomReconstruction; cols = 1:size(d, 2)) =
    aaft(d.reconstructed_pts; cols = cols)

iaaft(d::CustomReconstruction; cols = 1:size(d, 2)) =
    iaaft(d.reconstructed_pts; cols = cols)

export randomshuffle, randomphases, randomamplitudes, aaft, iaaft


text:
using OrdinaryDiffEq, Test, Statistics
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

@testset "Scalar Discrete Problem" begin
  prob = DiscreteProblem(0.5,(0.0,1.0))
  sol =solve(prob,FunctionMap())
  @test sol[1] == sol[end]
  @test sol(0.5:0.1:0.7) == [sol[1],sol[1],sol[1]]

  sol =solve(prob_ode_linear,FunctionMap())
  @test sol[end] ≈ .505
  sol =solve(prob_ode_linear,FunctionMap(scale_by_time=true),dt=1/4)
  sol2 =solve(prob_ode_linear,Euler(),dt=1/4)
  @test sol[end] == sol2[end]
  @test sol(0.53) != sol2(0.53)
end

@testset "Array Discrete Problem" begin
  prob2 = DiscreteProblem(rand(4,2),(0.0,1.0))
  sol =solve(prob2,FunctionMap())
  @test sol[1] == sol[end]
  @test sol(0.5) == sol[1]

  @test_nowarn sol =solve(prob_ode_2Dlinear,FunctionMap())
  @test_nowarn sol2 =solve(prob_ode_2Dlinear,Euler(),dt=1)
  sol =solve(prob_ode_2Dlinear,FunctionMap(scale_by_time=true),dt=1/4)
  sol2 =solve(prob_ode_2Dlinear,Euler(),dt=1/4)
  @test sol[end] == sol2[end]
  @test sol(0.35) != sol2(0.53)
end

@testset "Discrete Problem (Henon Map)" begin
  function henon_map!(u_next, u, _p, t)
    u_next[1] = 1 + u[2] - 1.4 * u[1]^2
    u_next[2] = 0.3 * u[1]
  end

  prob = DiscreteProblem(henon_map!, [0.5, 0.5], (0, 100)) # Integer time
  sol = solve(prob, FunctionMap())
  @test eltype(sol.t) <: Integer
end


text:
using NeuralStyle
using Test

@info "Testing NeuralStyle.jl ..."

@testset "NeuralStyle" begin

  @testset "Utils" begin
    include("utils.jl")
  end

end


text:
#### Entrainment-Detrainment kernels

"""
    entr_detr(
        m::AtmosModel{FT},
        entr::EntrainmentDetrainment,
        state::Vars,
        aux::Vars,
        t::Real,
        ts,
        env,
        i,
    ) where {FT}
Returns the dynamic entrainment and detrainment rates,
as well as the turbulent entrainment rate, following
Cohen et al. (JAMES, 2020), given:
 - `m`, an `AtmosModel`
 - `entr`, an `EntrainmentDetrainment` model
 - `state`, state variables
 - `aux`, auxiliary variables
 - `t`, the time
 - `ts`, NamedTuple of thermodynamic states
 - `env`, NamedTuple of environment variables
 - `i`, index of the updraft
"""
function entr_detr(
    m::AtmosModel{FT},
    entr::EntrainmentDetrainment,
    state::Vars,
    aux::Vars,
    t::Real,
    ts,
    env,
    i,
) where {FT}

    # Alias convention:
    gm = state
    en = state.turbconv.environment
    up = state.turbconv.updraft
    en_aux = aux.turbconv.environment
    up_aux = aux.turbconv.updraft

    N_up = n_updrafts(m.turbconv)
    ρ_inv = 1 / gm.ρ
    a_up_i = up[i].ρa * ρ_inv

    lim_ϵ = entr.lim_ϵ
    lim_amp = entr.lim_amp
    w_min = entr.w_min
    # precompute vars
    w_up_i = up[i].ρaw / up[i].ρa
    sqrt_tke = sqrt(max(en.ρatke, 0) * ρ_inv / env.a)
    # ensure far from zero
    Δw = filter_w(w_up_i - env.w, w_min)
    w_up_i = filter_w(w_up_i, w_min)

    Δb = up_aux[i].buoyancy - en_aux.buoyancy

    D_ε, D_δ, M_δ, M_ε =
        nondimensional_exchange_functions(m, entr, state, aux, t, ts, env, i)

    # I am commenting this out for now, to make sure there is no slowdown here
    Λ_w = abs(Δb / Δw)
    Λ_tke = entr.c_λ * abs(Δb / (max(en.ρatke * ρ_inv, 0) + w_min))
    λ = lamb_smooth_minimum(
        SVector(Λ_w, Λ_tke),
        m.turbconv.mix_len.smin_ub,
        m.turbconv.mix_len.smin_rm,
    )

    # compute limiters
    εt_lim = εt_limiter(w_up_i, lim_ϵ, lim_amp)
    ε_lim = ε_limiter(a_up_i, lim_ϵ, lim_amp)
    δ_lim = δ_limiter(a_up_i, lim_ϵ, lim_amp)
    # compute entrainment/detrainment components
    ε_trb =
        2 * a_up_i * entr.c_t * sqrt_tke /
        max((w_up_i * a_up_i * m.turbconv.pressure.H_up), entr.εt_min)
    ε_dyn = λ / w_up_i * (D_ε + M_ε) * ε_lim
    δ_dyn = λ / w_up_i * (D_δ + M_δ) * δ_lim

    ε_dyn = min(max(ε_dyn, FT(0)), FT(1))
    δ_dyn = min(max(δ_dyn, FT(0)), FT(1))
    ε_trb = min(max(ε_trb, FT(0)), FT(1))

    return ε_dyn, δ_dyn, ε_trb
end;

"""
    ε_limiter(a_up::FT, ϵ::FT) where {FT}
Returns the asymptotic value of entrainment
needed to ensure boundedness of the updraft
area fraction between 0 and 1, given
 - `a_up`, the updraft area fraction
 - `ϵ`, a minimum threshold value
"""
ε_limiter(a_up::FT, ϵ::FT, lim_amp::FT) where {FT} =
    FT(1) + lim_amp * exp(-a_up^2 / (2 * ϵ)) - exp(-(FT(1) - a_up)^2 / (2 * ϵ))

"""
    δ_limiter(a_up::FT, ϵ::FT) where {FT}
Returns the asymptotic value of detrainment
needed to ensure boundedness of the updraft
area fraction between 0 and 1, given
 - `a_up`, the updraft area fraction
 - `ϵ`, a minimum threshold value
"""
δ_limiter(a_up::FT, ϵ::FT, lim_amp::FT) where {FT} =
    FT(1) - exp(-a_up^2 / (2 * ϵ)) + lim_amp * exp(-(FT(1) - a_up)^2 / (2 * ϵ))

"""
    εt_limiter(w_up_i::FT, ϵ::FT) where {FT}
Returns the asymptotic value of turbulent
entrainment needed to ensure positiveness
of the updraft vertical velocity, given
 - `w_up_i`, the updraft vertical velocity
 - `ϵ`, a minimum threshold value
"""
εt_limiter(w_up_i::FT, ϵ::FT, lim_amp::FT) where {FT} =
    FT(1) + lim_amp * exp(-w_up_i^2 / (2 * ϵ))

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
"""
Subroutine daxpy

This subroutine computes constant times a vector plus a vector.
It uses unrolled loops for increments equal to one.
Jack Dongarra, LINPACK, 3/11/78.
"""
function tron_daxpy(n,da,dx,incx,dy,incy)
    if n <= 0
        return
    end
    if da == 0.0
        return
    end
    if incx != 1 || incy != 1
        ix = 1
        iy = 1
        if incx < 0
            ix = (-n+1)*incx + 1
        end
        if incy < 0
            iy = (-n+1)*incy + 1
        end
        for i=1:n
            dy[iy] += da*dx[ix]
            ix += incx
            iy += incy
        end
        return
    else

        # Code for both increments equal to 1

        m = mod(n, 4)
        if m != 0
            for i=1:m
                dy[i] += da*dx[i]
            end
        end
        if n < 4
            return
        end

        mp1 = m + 1
        for i=mp1:4:n
            dy[i] += da*dx[i]
            dy[i + 1] += da*dx[i + 1]
            dy[i + 2] += da*dx[i + 2]
            dy[i + 3] += da*dx[i + 3]
        end
        return
    end
end

if isequal(BLAS_LIBRARY, :OpenBlas)
    daxpy(n,da,dx,incx,dy,incy) = BLAS.axpy!(da,dx,dy)
    daxpy(n,da,dx::CuArray,incx,dy::CuArray,incy) = CUBLAS.axpy!(n, da, dx, dy)
else
    daxpy(n,da,dx,incx,dy,incy) = tron_daxpy(n,da,dx,incx,dy,incy)
end

@inline function daxpy(n::Int,da::Float64,
                       dx::CuDeviceArray{Float64,1},incx::Int,
                       dy::CuDeviceArray{Float64,1},incy::Int)
    tx = threadIdx().x
    ty = threadIdx().y

    if tx <= n && ty == 1
        @inbounds dy[tx] = dy[tx] + da*dx[tx]
    end
    CUDA.sync_threads()

    return
end


text:
using Test, Colors
using FixedPointNumbers

@testset "Parse" begin
    r8(x) = reinterpret(N0f8, x)

    # Color parsing
    # named-color
    redN0f8 = parse(Colorant, "red")
    @test colorant"red" === redN0f8
    @test redN0f8 === RGB{N0f8}(1,0,0)
    @test parse(RGB{Float64}, "red") === RGB{Float64}(1,0,0)
    @test isa(parse(HSV, "blue"), HSV)
    @test_throws ArgumentError parse(Colorant, "p ink")
    @test parse(Colorant, "transparent") === RGBA{N0f8}(0,0,0,0)
    @test parse(Colorant, "\nSeaGreen ") === RGB{N0f8}(r8(0x2E),r8(0x8B),r8(0x57))
    @test parse(Colorant, "sea GREEN") === colorant"seagreen"

    # hex-color
    @test parse(Colorant, "#D0FF58") === RGB(r8(0xD0),r8(0xFF),r8(0x58))
    @test parse(Colorant, "0xd0ff58") === RGB(r8(0xD0),r8(0xFF),r8(0x58))
    @test parse(Colorant, "#FB0") === RGB(r8(0xFF),r8(0xBB),r8(0x00))
    @test parse(Colorant, "#FB0A") === RGBA(r8(0xFF),r8(0xBB),r8(0x00),r8(0xAA))
    @test parse(Colorant, "0xFB0A") === ARGB(r8(0xBB),r8(0x00),r8(0xAA),r8(0xFF))
    @test parse(Colorant, "#FFBB00AA") === RGBA(r8(0xFF),r8(0xBB),r8(0x00),r8(0xAA))
    @test parse(Colorant, "0xFFBB00AA") === ARGB(r8(0xBB),r8(0x00),r8(0xAA),r8(0xFF))
    @test_throws ArgumentError parse(Colorant, "#BAD05")
    @test_throws ArgumentError parse(Colorant, "#BAD0007")
    @test_throws ArgumentError parse(Colorant, "#BAD000009")

    # rgb()
    @test parse(Colorant, "rgb(55,217,127)")      === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x7f))
    @test colorant"    Rgb(  55.1, 217, 127 )   " === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x7f))
    @test parse(Colorant, "rgb(22%,85%,50%)")     === RGB{N0f8}(r8(0x38),r8(0xd9),r8(0x80))
    @test parse(Colorant, "rgba(55,217,127,0.5)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5)
    @test parse(Colorant, "rgb( 55,217,127,50%)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5) # CSS Color Module Level 4
    @test parse(Colorant, "rgb( 55 217 127 /.5)") === RGBA{N0f8}(r8(0x37),r8(0xd9),r8(0x7f),0.5) # CSS Color Module Level 4
    @test_throws ArgumentError parse(Colorant, "rgb(55, 85%, 50%)") # this is invalid according to CSS spec.
    @test parse(Colorant, "rgb(21.6%,85%,50%)") === RGB{N0f8}(r8(0x37),r8(0xd9),r8(0x80))

    # hsl()
    @test parse(Colorant, "hsl(120, 100%, 50%)")  === HSL{Float32}(120,1.0,.5)
    @test colorant"    Hsl(    120, 100%, 50% ) " === HSL{Float32}(120,1.0,.5)
    @test parse(RGB{N0f8},"hsl(120, 100%, 50%)")  === convert(RGB{N0f8}, HSL{Float32}(120,1.0,.5))
    @test_throws ArgumentError parse(Colorant, "hsl(120, 100, 50)")
    @test_throws ArgumentError parse(Colorant, "hsl(120%,100%,50%)")
    @test parse(Colorant, "hsla(120,50%,7%, .6)") === HSLA{Float32}(120,.5,.07,.6)
    @test parse(Colorant, "hsl( 120,50%,7%,60%)") === HSLA{Float32}(120,.5,.07,.6) # CSS Color Module Level 4
    @test parse(Colorant, "hsl( 120 50% 7% / 1)") === HSLA{Float32}(120,.5,.07, 1) # CSS Color Module Level 4
    @test parse(Colorant, "hsl(         90.0, 100%,   0%)") === HSL{Float32}(90, 1, 0)
    @test parse(Colorant, "hsl(        90Deg, 100%,   0%)") === HSL{Float32}(90, 1, 0)
    @test parse(Colorant, "hsl(     0.25turn, 120%, 0.0%)") === HSL{Float32}(90, 1, 0)
    @test parse(Colorant, "hsl(1.57079633RAD, 100%, -10%)") === HSL{Float32}(90, 1, 0)
    @test parse(Colorant, "hsl(      100grad, 100%, 0e2%)") === HSL{Float32}(90, 1, 0)

    @test parse(Colorant, :red) === colorant"red"
    @test_deprecated parse(Colorant, colorant"red") === colorant"red"

end


text:
import Base: abs, *, /, +, -, convert
import LinearAlgebra: dot, cross, transpose

struct Point{T}
    x::T
    y::T
    z::T

    function Point(x::Number, y::Number, z::Number)
        x_p, y_p, z_p = promote(x, y, z)
        return new{typeof(x_p)}(x_p, y_p, z_p)
    end
end

convert(::Type{Point{T}}, p::Point) where T = Point(convert(T, p.x), convert(T, p.y), convert(T, p.z))

Point(x::SVector) = Point(x[1], x[2], x[3])

svector(p::Point) = SVector((p.x, p.y, p.z))

cross(a, b) = Point(cross(svector(a), svector(b)))

+(a::Point, b::Point) = Point(a.x + b.x, a.y + b.y, a.z + b.z)
-(a::Point, b::Point) = Point(a.x - b.x, a.y - b.y, a.z - b.z)

-(p::Point) = Point(-p.x, -p.y, -p.z)

abs(p::Point) = norm(svector(p))

*(s::Number, p::Point) = Point(s * p.x, s * p.y, s * p.z)

/(p::Point, s::Number) = Point(p.x / s, p.y / s, p.z / s)

@inline *(A::SMatrix{3,3}, p::Point) = Point(A * svector(p))
@inline *(A::SDiagonal{3}, p::Point) = Point(A * svector(p))

conjugate(p::Point) = Point(p.x', p.y', p.z')

dot(a::Point, b::Point) = a.x * b.x + a.y * b.y + a.z * b.z

normalize_abs(p::Point) = inv(abs(p)) * p

normalize_z(p::Point) = inv(p.z) * p

distance(a::Point, b::Point) = abs(normalize_z(a) - normalize_z(b))

function cross_operator(p::Point)
    return @SMatrix [ 0   -p.z  p.y;
                      p.z  0   -p.x;
                     -p.y  p.x  0]
end

const CircleMatrix{T} = SMatrix{3, 3, T, 9}

center_of_circle(c::CircleMatrix) = Point(c[3,1], c[3,2], -c[1,1])

function circle_from_mp(m::Point, p::Point)
    x = m.x
    y = m.y

    mz = -m.z

    tang = @SMatrix [
        mz 0  x;
        0  mz y;
        x  y  0
    ]

    fundmat = SDiagonal(SVector(0, 0, 1))

    mu = dot(tang * p, p)

    la = dot(fundmat * p, p)

    m1 = mu * fundmat
    m2 = la * tang

    return m1 - m2
end

function circle_from_mr(m::Point, r::Real)
    sr = m.z * r
    rr = sr * sr

    return @SMatrix [
        m.z*m.z   0       -m.z*m.x;
        0         m.z*m.z -m.z*m.y;
        -m.z*m.x -m.z*m.y  m.x*m.x + m.y*m.y - rr
    ]
end

function intersect_line_circle(l::Point, c::CircleMatrix)
    l1 = cross_operator(l)

    s = transpose(l1) * (c * l1)

    maxidx = argmax((abs(l.x), abs(l.y), abs(l.z)))

    if maxidx == 1
        a11 = s[2,2]
        a12 = s[2,3]
        a21 = s[3,2]
        a22 = s[3,3]
        b = l.x
    elseif maxidx == 2
        a11 = s[1,1]
        a12 = s[1,3]
        a21 = s[3,1]
        a22 = s[3,3]
        b = l.y
    else
        a11 = s[1,1]
        a12 = s[1,2]
        a21 = s[2,1]
        a22 = s[2,2]
        b = l.z
    end

    alp = sqrt(a12 * a21 - a11 * a22) / b

    erg = s + alp * l1

    erg1 = Point(erg[argmax((norm(erg[1,:]), norm(erg[2,:]), norm(erg[3,:]))),:])
    erg2 = Point(erg[:,argmax((norm(erg[:,1]), norm(erg[:,2]), norm(erg[:,3])))])

    return erg1 / abs(erg1), erg2 / abs(erg2)
end

function intersect_circle_circle(c1::CircleMatrix, c2::CircleMatrix)
    ct1 = c2[1,1]
    line1 = ct1 * c1[3,:]

    ct2 = c1[1,1]
    line2 = ct2 * c2[3,:]

    ll = line1 - line2

    return intersect_line_circle(Point(ll[1], ll[2], ll[3] / 2), c1)
end

function projective_dist_min_scal(a::Point, b::Point)
    sa = abs(a)
    sb = abs(b)

    if sa == 0 || sb == 0
        return zero(sa)
    end

    cb = conjugate(b)

    p = dot(a, cb)

    if abs(p) < 1e-13
        #error("p too small")
        p = 1.0
    end

    np = p / abs(p)

    na = a / sa
    nb = np * (b / sb)

    d1 = abs(na + nb)
    d2 = abs(na - nb)

    return min(d1, d2)
end

function tracing2(old::Tuple{Point{T},Point{T}}, new::Tuple{Point{T},Point{T}}) where T
    o1, o2 = old
    n1, n2 = new

    safety = 3

    do1n1 = projective_dist_min_scal(o1, n1)
    do1n2 = projective_dist_min_scal(o1, n2)
    do2n1 = projective_dist_min_scal(o2, n1)
    do2n2 = projective_dist_min_scal(o2, n2)
    do1o2 = projective_dist_min_scal(o1, o2)
    dn1n2 = projective_dist_min_scal(n1, n2)
    cost1 = do1n1 + do2n2
    cost2 = do1n2 + do2n1

    # Always sort output: we don't know yet whether it's correct, but
    # it's our best bet.

    if cost1 > cost2
        res = (n2, n1)
        cost = cost2
    else
        res = (n1, n2)
        cost = cost1
    end

    # is nan test?

    if do1o2 > cost * safety && dn1n2 > cost * safety
        # Distance within matching considerably smaller than distance
        # across matching, so we could probably match correctly.
        return res
    elseif dn1n2 < 1e-5
        error("New points too close: we presumably are inside a singularity.")
    elseif do1o2 < 1e-5
        error("Moved out of singularity.")
    end

    # refine
    return nothing
end

function project_point_to_line(point::Point, line::Point)
    return cross(cross(Point(line.x, line.y, 0), point), line)
end

function point_reflection(center::Point, point::Point)
    return 2 * point.z * center - center.z * point
end


text:
include("6.3.10 (c+d x)^m (a+b tanh)^n.jl")
include("6.3.11 (e x)^m (a+b tanh(c+d x^n))^p.jl")
include("6.3.12 (d+e x)^m tanh(a+b x+c x^2)^n.jl")

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.3 Hyperbolic tangent/.jl


text:
# ==============================================================================
# Trop{T} types

# Fake templates to make difference between Min-Plus and Max-Plus numbers
struct Min end
struct Max end
const MM = Union{Min, Max}

# Base class for Max-Plus and Min-Plus structures
struct Trop{T <: MM} <: Number
    v::Float64
end

# Max-Plus structure
const MP = Trop{Max}

# Min-Plus structure
const MI = Trop{Min}

# Sparse matrix (classic algebra) shorter name
const SparseMatrix{T,U} = SparseMatrixCSC{T,U}

# Sparse matrix (max-plus algebra) shorter name
const SpaTrop{T,U} = SparseMatrix{Trop{T},U}
const SpaMP{U} = SparseMatrix{MP,U}
const SpaMI{U} = SparseMatrix{MI,U}

# Sparse vector (max-plus algebra) shorter name
const SpvTrop{T,U} = SparseVector{Trop{T},U}
const SpvMP{U} = SparseVector{MP,U}
const SpvMI{U} = SparseVector{MI,U}

# Dense matrix (max-plus algebra) shorter name
const ArrTrop{T,N} = Array{Trop{T},N}
const ArrMP{N} = Array{MP,N}
const ArrMI{N} = Array{MI,N}
const AbsMatMP = AbstractMatrix{MP}

# Dense vector (max-plus algebra) shorter name
const VecTrop{T} = Vector{Trop{T}}
const VecMP = Vector{MP}
const VecMI = Vector{MI}


text:
using Tracker, ForwardDiff, Zygote, DistributionsAD, Bijectors, ReverseDiff
using Random, LinearAlgebra, Combinatorics, Test
using DistributionsAD: TuringUniform, TuringMvNormal, TuringMvLogNormal,
                        TuringPoissonBinomial

Random.seed!(123456)

function get_stage()
    if get(ENV, "TRAVIS", "") == "true" || get(ENV, "GITHUB_ACTIONS", "") == "true"
        if "STAGE" in keys(ENV)
            return ENV["STAGE"]
        else
            return "nonAD"
        end
    end

    return "all"
end

stg = get_stage()
if stg in ("nonAD", "all")
    @testset "Interface" begin
        include("interface.jl")
    end
    include("transform.jl")
    include("norm_flows.jl")
    include("bijectors/permute.jl")
end
if stg != "nonAD"
    include("ad/ad_test_utils.jl")
    include("ad/distributions.jl")
end


text:
struct MaximumState
    iteration::Int
    numintervals::Int
    maxenclosure::arb
    metadata::Dict
end

function MaximumState(i, numintervals, maxenclosure)
    MaximumState(Int(i), numintervals, maxenclosure, Dict())
end

struct MaximumTrace
    states::Vector{MaximumState}
end

MaximumTrace() = MaximumTrace(Array{MaximumState}(undef, 0))

function Base.show(io::IO, s::MaximumState)
    @printf io "%6d   %11d    %s\n" s.iteration s.numintervals string(s.maxenclosure)
    return
end

function Base.show(io::IO, t::MaximumTrace)
    @printf io "Iter     Numintervals   Enclosure \n"
    @printf io "------   ------------   --------------\n"
    for state in t.states
        show(io, state)
    end
    return
end

Base.push!(t::MaximumTrace, s::MaximumState) = push!(t.states, s)

function update!(tr::MaximumTrace,
                 st::MaximumState,
                 store_trace::Bool,
                 show_trace::Bool,
                 callback = nothing)
    if store_trace
        push!(tr, st)
    end
    if show_trace
        show(st)
    end
    if !isnothing(callback)
        callback(st)
    end
    return
end

@recipe function f(st::MaximumState;
                   level = nothing)
    xs = [Float64(interval[j]) for j in 1:2, interval in st.metadata["intervals"]]
    if isnothing(level)
        ys = [[Float64(eval[j]) for eval in st.metadata["evaluations"]] for j in 1:2]
    else
        ys = [[level + j for _ in st.metadata["evaluations"]] for j in 1:2]
    end
    seriescolor --> :blue
    label := ""
    fillrange := ys[2]'
    (xs, ys[1]')
end

@userplot TracePlot

@recipe function f(h::TracePlot; cname = "oranges")
    if typeof(h.args[1]) == MaximumTrace
        states = h.args[1].states
    elseif typeof(h.args[1]) == Vector{MaximumState}
        states = h.args[1]
    else
        error("Trace plots should be given a trace or vector of states.  Got: $(typeof(h.args[1]))")
    end

    N = length(states)

    coloroffset = 1
    cm = colormap(cname, coloroffset + N)

    for i in 1:N
        @series begin
            seriescolor := cm[i + coloroffset]
            level := i
            states[i]
        end
    end
end