text:
"""
*S*ingle *P*ole *A*pproximation for extra nonlinear scattering solutions
"""
module SPA

export spa

using ...Common
using LinearAlgebra
using NLsolve
using RecipesBase

import ..HelmholtzNLS
import ..HelmholtzCF
import ..helmholtzeigen
import ..scattering!
import ...Common.SPA_ACCEPTANCE_PHASE

################################################################################
# spa main

"""
    spa(nls; [refine=false, nev=1, verbose=false])

Single Pole Approximation of nonlinear scattering problem `nls` for solutions besides
that given in `nls.ψ`.

`refine` | if `true`, additional refinement through nonlinear CF equation
`nev=1` | number of CF states used in refinement, making it a "few pole approximation".
`verbose` | if `true`, shows traces for CF nonlinear refinement
"""
function spa(nls::HelmholtzNLS; refine::Bool=false, forcerefine::Bool=false, nev::Integer=1, takemeout=false, nev_th::Integer=100, kwargs...)
    sim = nls.simulation

    # generate CF states
        cf = HelmholtzCF(sim)
        if nev > nev_th
            u = ScalarField(sim, 2nev_th + nev + 1)
            N = nev_th
            η1, u1 = helmholtzeigen(cf, nls.ω, [nls.ω], nls.ψ; nev=nev_th)
            η = fill(η1[1],size(u,2))
            for i ∈ 1:N u.values[:,i] = u1[:,i] end
            for i ∈ 1:N η[i] = η1[i] end
            F = Vector(diag(cf.F))
            while nev > N
                η1, u1 = helmholtzeigen(cf, nls.ω, [nls.ω], nls.ψ; η=minimum(real(η)), nev=nev_th)
                η2, u2 = helmholtzeigen(cf, nls.ω, [nls.ω], nls.ψ; η=maximum(real(η)), nev=nev_th)
                dinds1 = Int[]
                dinds2 = Int[]
                for i ∈ 1:N
                    for j ∈ 1:nev_th
                        if abs(η[i] - η1[j]) < 1e-1
                            if abs(1-abs(sum(u.values[:,i].*u1.values[:,j].*F)*sim.dx)) < 1e-1
                                push!(dinds1,j)
                            end
                        end
                    end
                    for j ∈ 1:nev_th
                        if abs(η[i] - η2[j]) < 1e-1
                            if abs(1-abs(sum(u.values[:,i].*u2.values[:,j].*F)*sim.dx)) < 1e-1
                                push!(dinds2,j)
                            end
                        end
                    end
                end
                kinds1 = setdiff(1:nev_th,dinds1)
                kinds2 = setdiff(1:nev_th,dinds2)
                sort!(dinds1)
                sort!(dinds2)
                deleteat!(η1,dinds1)
                deleteat!(η2,dinds2)
                N1 = length(η1)
                N2 = length(η2)
                println("size(η): ",size(η))
                println("N: ",N)
                println("N1: ",N1)
                println("N2: ",N2)
                for i ∈ 1:N1 η[N + i] = η1[i] end
                for i ∈ 1:N2 η[N + N1 + i] = η2[i] end
                for i ∈ 1:N1 u.values[:,N+i] = u1.values[:,kinds1[i]] end
                for i ∈ 1:N2 u.values[:,N+N1+i] = u2.values[:,kinds2[i]] end
                N += N1 + N2
            end
            η = η[1:N]
            perm = sortperm(η,by=abs)
            u0 = u
            u = ScalarField(sim,nev)
            for i ∈ 1:nev u.values[:,i] = u0.values[:,perm[i]] end
            η = η[perm[1:nev]]
        else
            η, u = helmholtzeigen(cf, nls.ω, [nls.ω], nls.ψ; nev=nev)
            F = Vector(diag(cf.F))
        end

    # initialize CFAmplitude structure with CF eigenpairs
        cfa1 = CFAmplitudeVector(nls, η, u)
        cfa2 = CFAmplitudeVector(nls, η, u)
        cfa3 = CFAmplitudeVector(nls, η, u)
        cfa4 = CFAmplitudeVector(nls, η, u)

    # check that nonlinear solution exists
    if nls.converged[]
        # single-cf approximations
            phase = angle(sum(F.* abs2.(u[:,1]) .* u[:,1] .* nls.solution.total[:]))
            τ = cis(-phase)*nls.solution.total[:]

            A = η[1]*cfa1.Γ*sum(F .* abs2.(u[:,1]) .* u[:,1].^2)*sim.dx
            B = (η[1] + cfa1.γ*cfa1.D₀)*cfa1.Γ*sum(F .* abs2.(u[:,1]) .* u[:,1] .* τ)*sim.dx
            C = η[1]*cfa1.Γ*sum(F .* u[:,1].^3 .* conj(τ))*sim.dx
            E = cfa1.γ*cfa1.D₀*cfa1.Γ*sum(F .* abs2.(u[:,1]) .* τ.^2)*sim.dx
            G = -cfa1.γ*cfa1.D₀ + η[1]*sum(sim.F .* u[:,1].^2)*sim.dx # yes, this one supposed to be sim.F

            small_angle_condition = abs(sin(angle(G/A))) ≤ abs(sin(SPA_ACCEPTANCE_PHASE))

            b = real(B/A)
            c = real(C/A)
            e = real(E/A)
            g = real(G/A)

            Σ = (b+c); cfa1.Σ = Σ; cfa2.Σ = Σ
            Π = (e+g); cfa1.Π = Π; cfa2.Π = Π
            Δ = (Σ/2)^2 - Π; cfa1.Δ = Δ; cfa2.Δ = Δ

            cfa1.scf_estimate = -cis(phase)*(Σ/2 - sqrt(complex(Δ)))
            cfa2.scf_estimate = -cis(phase)*(Σ/2 + sqrt(complex(Δ)))

            # can't exect to resolve better than imaginary parts of these, since in theory they should vanish
            discriminant_condition = Δ > 2maximum(abs∘imag,[(B+C)/A,(E+G)/A])

        if discriminant_condition && small_angle_condition
            cfa1.multivalued = true
            cfa2.multivalued = true
        else
            cfa1.multivalued = false
            cfa2.multivalued = false
        end

        # single pole approximations
            # cfa1
            ρ, θ = abs(cfa1.scf_estimate), angle(cfa1.scf_estimate)
            spa1 = SPAsolver(cfa1)
            try # try spa, if works record results, otherwise declare single-valued
                res = nlsolve(spa1,[ρ,θ]; method=:anderson, m=1, iterations=100, ftol=1e-6)
                if converged(res)
                    cfa1.converged = true
                    cfa1.multivalued = true
                    cfa1.spa_estimate = res.zero[1]*cis(res.zero[2])
                else
                    cfa1.converged = false
                    cfa1.multivalued = false
                end
            catch
                cfa1.converged = false
                cfa1.multivalued = false
            end
            if takemeout
            function f!(Fg,x)
                a = Vector{ComplexF64}(undef,length(x)÷2)
                for i ∈ eachindex(a) a[i] = complex(x[2(i-1)+1],x[2(i-1)+2]) end
                fill!(cfa1.a,0)
                copyto!(cfa1.a,a)
                ψ = cfa1(length(a))
                # Fg[1],Fg[2] = reim(A*abs2(a)+cis(phase)*B*conj(a)+cis(-phase)*C*a+cis(2phase)*E*conj(a)/a+G)
                hψ = 1 .+ cfa1.Γ .* abs2.(ψ.values)
                hφ = 1 .+ cfa1.Γ .* abs2.(cfa1.φ.values)
                q = zeros(ComplexF64,size(u,1))
                for i ∈ eachindex(a) q += a[i]*η[i]*cfa1.u[:,i] end
                for i ∈ eachindex(a)
                    A = sum(sim.F .* ψ.values  .* cfa1.u.values[:,i])*sim.dx*cfa1.γ*cfa1.D₀
                    B = sum(F .* cfa1.φ[:] .* cfa1.u.values[:,i] .* hψ)*sim.dx*cfa1.γ*cfa1.D₀
                    H = sum(hψ .* F.* q .* u[:,i])*sim.dx
                    Fg[2(i-1)+1], Fg[2(i-1)+2] = reim((-H+A-B)/norm(a))
                end
            end
            x = Vector{Float64}(undef,2nev)
            for i ∈ 1:nev x[2(i-1)+1], x[2(i-1)+2] = reim(cfa1.a[i]) end
            return nlsolve(f!,x)
        end

            #cfa2
            ρ, θ = abs(cfa2.scf_estimate), angle(cfa2.scf_estimate)
            spa2 = SPAsolver(cfa2)
            try # try spa, if works record results, otherwise declare single-valued
                res = nlsolve(spa2,[ρ,θ]; method=:anderson, m=1, iterations=100, ftol=1e-6)
                if converged(res)
                    cfa2.converged = true
                    cfa2.multivalued = true
                    cfa2.spa_estimate = res.zero[1]*cis(res.zero[2])
                else
                    cfa2.converged = false
                    cfa2.multivalued = false
                end
            catch
                cfa2.converged = false
                cfa2.multivalued = false
            end

        if forcerefine || (refine && cfa1.multivalued)
            #cfa1
                cfa1.refined = true
                x = a_to_x(cfa1,nev)
                res = nlsolve(cfa1, x; kwargs...)
                cfa1.converged = converged(res)
                x_to_a!(cfa1,res.zero)

            converged(res) ? append!(cfa2.a0,cfa1.a) : nothing

            #cfa2
                cfa2.refined = true
                x = a_to_x(cfa2,nev)
                res = nlsolve(cfa2, x; kwargs...)
                cfa2.converged = converged(res)
                x_to_a!(cfa2,res.zero)
        end

        ρ² = (e-g)/(1+c*(c-b)/e)
        if ρ² > 0
            ρ = sqrt(ρ²)
            cosθ = (c-b)*ρ/2/e
            if cosθ^2 ≤ 1
                θ = phase + acos(cosθ)
                cfa3.Σ, cfa3.Π = ρ, θ
                θ = phase - acos(cosθ)
                cfa4.Σ, cfa4.Π = ρ, θ
                cfa3.scf_estimate = cis(real(cfa3.θ))*real(cfa3.ρ)
                cfa4.scf_estimate = cis(real(cfa4.θ))*real(cfa4.ρ)
                cfa3.multivalued = true
                cfa4.multivalued = true
            else
                cfa3.multivalued = false
                cfa4.multivalued = false
            end
        else
            cfa3.multivalued = false
            cfa4.multivalued = false
        end

        if forcerefine || (refine && cfa3.multivalued)
            #cfa3
                cfa3.refined = true
                x = a_to_x(cfa3,nev)
                res = nlsolve(cfa3, x; kwargs...)
                cfa3.converged = converged(res)
                x_to_a!(cfa3,res.zero)

            converged(res) ? append!(cfa4.a0,cfa3.a) : nothing

            #cfa4
                cfa4.refined = true
                x = a_to_x(cfa4,nev)
                res = nlsolve(cfa4, x; kwargs...)
                cfa4.converged = converged(res)
                x_to_a!(cfa4,res.zero)
        end
    else #if get(kwargs,:show_trace,false)
        @warn "must solve NonLinearScatteringProblem first"
    end

    return cfa1, cfa2, cfa3, cfa4
end

################################################################################
# structures

# structure used in nonlinear CF solve
struct CFAmplitudeVector{TF}
    scf_estimate::Base.RefValue{ComplexF64}
    spa_estimate::Base.RefValue{ComplexF64}
    Σ::Base.RefValue{Float64}
    Π::Base.RefValue{Float64}
    Δ::Base.RefValue{Float64}
    a::Vector{ComplexF64}
    a0::Vector{ComplexF64}
    φ::TF
    ζ::TF
    u::TF
    η::Vector{ComplexF64}
    B::Vector{ComplexF64}
    dx::Float64
    D₀::Float64
    γ::ComplexF64
    Γ::Float64
    F::Vector{Float64}
    converged::Base.RefValue{Bool}
    multivalued::Base.RefValue{Bool}
    refined::Base.RefValue{Bool}
end

"""
    CFAmplitudeVector(nls, η, u, [a0]) -> cfa

construct nonlinear CF solver object with CF evals `η` and efuns `u`.
`estimate` is the SPA estimate, `a0` is a previously found solution to be avoided
"""
function CFAmplitudeVector(
            nls::HelmholtzNLS{N},
            η::Vector,
            u::ScalarField{N},
            a0::Vector=ComplexF64[]
            ) where N

    nev = size(u,2)
    sim = nls.simulation
    φ = nls.solution.total

    tls = nls.simulation.dispersive_domains[1].χ
    Γ = Common.Dispersions.TwoLevelSystems.Γ(tls,nls.ω)
    γ = Common.Dispersions.TwoLevelSystems.γ(tls,nls.ω)
    D₀ = tls.D₀

    a = zeros(ComplexF64,nev)

    ζ = ScalarField(sim,1)
    B = Vector{ComplexF64}(undef, nev)
    Fφ = sim.F .* φ[:] ./ (1 .+ Γ .* abs2.(φ[:]) )
    @inbounds for i ∈ eachindex(B) B[i] = γ*D₀*sum(Fφ .* u[:,i]) * sim.dx end
    return CFAmplitudeVector(Ref(complex(0.0)), Ref(complex(0.0)), Ref(0.0), Ref(0.0), Ref(0.0), a, a0, φ, ζ, u, η, B, sim.dx, D₀, γ, Γ, sim.F, Ref(true), Ref(false), Ref(false))
end

# create wavefunction φ + ∑ⱼaⱼuⱼ
"""
    (cfa)([nev]) -> φ

Generate `φ` from background field in `cfa` plus CF corrections
"""
function (cfa::CFAmplitudeVector)(nev::Integer=length(cfa.a))
    copyto!(cfa.ζ.values,cfa.φ.values)
    @inbounds for i ∈ 1:nev
        @inbounds @fastmath @simd for j ∈ eachindex(cfa.ζ)
            cfa.ζ[j] = cfa.ζ[j] + cfa.a[i]*cfa.u[j,i]
        end
    end
    return cfa.ζ
end

# cfa(F,x) used in nlsolve.
function (cfa::CFAmplitudeVector)(F::Vector{Float64},x::Vector{Float64})
    n = length(x)÷2
    # load data
        x_to_a!(cfa,x)
    # create wavefunction associated with a's
        ψ = cfa(n)
    # normalization penalizes being close to 0 or a0
        𝒩⁻¹ = 1/norm(cfa.a) + (isempty(cfa.a0) ? 0.0 : 1/norm(cfa.a-cfa.a0) )
    @inbounds @fastmath for i ∈ 1:n
        A = complex(0.0)
        @inbounds @fastmath @simd ivdep for j ∈ eachindex(cfa.F)
            A += cfa.F[j] * ψ.values[j] * cfa.u.values[j,i] / (1 + cfa.Γ * abs2(ψ.values[j]) )
        end
        A *= cfa.γ*cfa.D₀*cfa.dx
        F[i], F[n+i] = reim((-cfa.η[i]*cfa.a[i] + A - cfa.B[i])*𝒩⁻¹)
    end
    return nothing
end

function Base.getproperty(cfa::CFAmplitudeVector, sym::Symbol)
    if Base.sym_in(sym,(:Σ,:σ,:sigma,:Sigma,:s,:S,:ρ,:rho,:r))
        return getfield(cfa,:Σ)[]
    elseif Base.sym_in(sym,(:Π,:π,:Pi,:pi,:P,:p,:θ,:ϑ,:theta))
        return getfield(cfa,:Π)[]
    elseif Base.sym_in(sym,(:Δ,:δ,:d,:D,:delta,:Delta,:discriminant,:Discriminant,:dis,:Dis))
        return getfield(cfa,:Δ)[]
    elseif Base.sym_in(sym, (:scf_estimate, :spa_estimate, :converged, :multivalued, :refined))
        return getfield(cfa,sym)[]
    else
        return getfield(cfa,sym)
    end
end

function Base.setproperty!(cfa::CFAmplitudeVector, sym::Symbol, x::Number)
    if Base.sym_in(sym,(:scf_estimate, :spa_estimate))
        getfield(cfa,sym)[] = x
        return cfa.a[1] = x
    elseif Base.sym_in(sym, (:Σ, :Π, :Δ, :converged, :multivalued, :refined))
        getfield(cfa,sym)[] = x
    else
        return setfield!(cfa,sym,x)
    end
end

################################################################################
# data loaders

a_to_x(cfa::CFAmplitudeVector,n::Integer) = a_to_x(cfa.a,n)
a_to_x(a::Vector,n::Integer) = a_to_x!(Vector{Float64}(undef,2n),a)

a_to_x!(x::Vector, cfa::CFAmplitudeVector) = a_to_x!(x,cfa.a)
@inline function a_to_x!(x::Vector,a::Vector)
    n = length(x)÷2
    @inbounds @simd ivdep for i ∈ 1:n x[i], x[n+i] = reim(a[i]) end
    return x
end

x_to_a!(cfa::CFAmplitudeVector, x::Vector) = x_to_a!(cfa.a,x)
@inline function x_to_a!(a::Vector,x::Vector)
    n = length(x)÷2
    @inbounds @simd ivdep for i ∈ 1:n a[i] = complex(x[i], x[n+i]) end
    return a
end

################################################################################

struct SPAsolver
    Γ::Float64
    φ::Vector{ComplexF64}
    u::Vector{ComplexF64}
    η::ComplexF64
    F::Vector{Float64}
    dx::Float64
    γ::ComplexF64
    D₀::Float64
    h::Vector{Float64}
end

function SPAsolver(cfa::CFAmplitudeVector)
    F = cfa.F ./ (1 .+ cfa.Γ .* abs2.(cfa.φ[:]) )
    SPAsolver(cfa.Γ, cfa.φ.values[:], cfa.u.values[:,1], cfa.η[1], F, cfa.dx, cfa.γ, cfa.D₀, Vector{Float64}(undef,length(F)))
end

function (spa::SPAsolver)(F::Vector,x::Vector)
    a = x[1]*cis(x[2])
    for i ∈ eachindex(spa.h) spa.h[i] = 1/(1 + spa.Γ*abs2( spa.φ[i] + a*spa.u[i])) end
    I1 = sum(spa.F .* spa.u.^2 .* spa.h)*spa.dx
    I2 = sum(spa.F .* abs2.(spa.u) .* spa.φ.^2 .* spa.h)*spa.dx*spa.Γ
    I3 = sum(spa.F .* abs2.(spa.u) .* spa.u .* spa.φ .* spa.h)*spa.dx*spa.Γ
    b, β = reim((I1 - spa.η/(spa.γ*spa.D₀))/I3)
    c, κ = reim(-I2/I3)
    θ = -atan(β+κ,b-c)
    ρ = (b+c)*cos(θ) - (β-κ)*sin(θ)
    F[1], F[2] = ρ-x[1], θ-x[2]
end

################################################################################
# Pretty Printing

import ..Common.PRINTED_COLOR_LIGHT
import ..Common.PRINTED_COLOR_GOOD
import ..Common.PRINTED_COLOR_NUMBER
import ..Common.PRINTED_COLOR_BAD
import ..Common.PRINTED_COLOR_WARN
import ..Common.PRINTED_COLOR_VARIABLE
import ..Common.PRINTED_COLOR_INSTRUCTION

function Base.show(io::IO,cfa::CFAmplitudeVector)
    printstyled(io,"CF Amplitudes ",color=PRINTED_COLOR_LIGHT)
    if cfa.refined[]
        if cfa.converged[]
            printstyled(io,"converged ", color=PRINTED_COLOR_GOOD)
        else
            printstyled(io,"not converged ", color=PRINTED_COLOR_BAD)
        end
    end
    if cfa.multivalued[]
        printstyled(io,"multi-valued ", color=PRINTED_COLOR_WARN)
        show(IOContext(io,:compact=>true),cfa.a)
    else
        printstyled(io,"single-valued", color=PRINTED_COLOR_INSTRUCTION)
    end
end

function Base.show(io::IO,mime::MIME"text/plain",cfa::CFAmplitudeVector)
    printstyled(io,"CF Amplitudes ",color=PRINTED_COLOR_LIGHT)
    if cfa.refined[]
        if cfa.converged[]
            printstyled(io,"converged ", color=PRINTED_COLOR_GOOD)
        else
            printstyled(io,"not converged ", color=PRINTED_COLOR_BAD)
        end
    end
    if cfa.multivalued[]
        printstyled(io,"multi-valued ", color=PRINTED_COLOR_WARN)
        print(io,"(SPA estimate of ")
        printstyled(io,"a[1]",color=PRINTED_COLOR_VARIABLE)
        print(io,": ")
        printstyled(IOContext(io,:compact=>true), cfa.spa_estimate, color=PRINTED_COLOR_NUMBER)
        print(io,")")
        print(io,"\nVector ")
        printstyled(io,"a",color=PRINTED_COLOR_VARIABLE)
        print(io,": ")
        show(IOContext(io,:compact=>true),mime,cfa.a)
    else
        printstyled(io,"single-valued", color=PRINTED_COLOR_INSTRUCTION)
    end
end

################################################################################
# Plotting

@recipe function f(cfa::CFAmplitudeVector)
    cfa()
end

end # module

using .SPA


text:
using SingleFloats
using Test

include("frombase.jl")

@test Float64(sqrt(Single32(2))) == sqrt(2.0)
@test Float32(sqrt(Single32(2.0))) == Float32(sqrt(Float64(2.0)))

#@test sum(cot.(Single32.(collect(1.0:20)))) == sum(cot.(reverse(Single32.(collect(1.0:20)))))

ext: jl


text:
# DiffusionMap.jl
module DiffusionMap

# dependency
using LinearAlgebra, Arpack, KrylovKit, SparseArrays

export
# kernels.jl
AbstractKernel, InverseDistanceKernel,
GaussianKernel, CosineKernel, CustomKernel,

# laplacians.jl
AbstractLaplacian, RowNormalizedLaplacian,
Adjacency, NormalizedAdjacency,
SymmetricLaplacian, RegularLaplacian,

# eigensolvers.jl
AbstractEigenSolver, FullEigen,
ArpackEigen, KrylovEigen,

# diffusionTypes.jl

Diffusionmap,

# diffusionCalculation.jl
get_laplacian, get_adjacency, thresholding!, solve

include("kernels.jl")
include("laplacians.jl")
include("eigensolvers.jl")
include("diffusionTypes.jl")
include("diffusionCalculation.jl")
# include("utils.jl")

end # module


text:
# This file is a part of GERDAMetadata.jl, licensed under the MIT License (MIT).

export DataLocations

mutable struct DataLocations
    raw::String
    gen::String
    meta::String
end

DataLocations(p::PropDict) = DataLocations(p.raw, p.gen, p.meta)

Base.convert(::Type{DataLocations}, p::PropDict) = DataLocations(p)

export SetupConfig

mutable struct SetupConfig
    data::DataLocations
end

SetupConfig(p::PropDict) = SetupConfig(convert(DataLocations, p.data))

Base.convert(::Type{SetupConfig}, p::PropDict) = SetupConfig(p)

export DataConfig

mutable struct DataConfig
    setups::Dict{Symbol,SetupConfig}
end

function DataConfig(p::PropDict)
    setups = Dict{Symbol,SetupConfig}()
    for (k, v) in p.setups
        setups[k] = SetupConfig(v)
    end
    DataConfig(setups)
end

Base.convert(::Type{DataConfig}, p::PropDict) = DataConfig(p)

function Base.read(::Type{DataConfig}, filename::AbstractString)
    p = read(PropDict, filename, subst_pathvar = true, subst_env = true)
    convert(DataConfig, p)
end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

## Auxiliary function ##

"""
    resize_coeffs1!{T<Number}(coeffs::Array{T,1}, order::Int)

If the length of `coeffs` is smaller than `order+1`, it resizes
`coeffs` appropriately filling it with zeros.
"""
function resize_coeffs1!(coeffs::Array{T,1}, order::Int) where {T<:Number}
    lencoef = length(coeffs)
    resize!(coeffs, order+1)
    if order > lencoef-1
        z = zero(coeffs[1])
        @simd for ord in lencoef+1:order+1
            @inbounds coeffs[ord] = z
        end
    end
    return nothing
end

"""
    resize_coeffsHP!{T<Number}(coeffs::Array{T,1}, order::Int)

If the length of `coeffs` is smaller than the number of coefficients
correspondinf to `order` (given by `size_table[order+1]`), it resizes
`coeffs` appropriately filling it with zeros.
"""
function resize_coeffsHP!(coeffs::Array{T,1}, order::Int) where {T<:Number}
    lencoef = length( coeffs )
    @inbounds num_coeffs = size_table[order+1]
    @assert order ≤ get_order() && lencoef ≤ num_coeffs
    num_coeffs == lencoef && return nothing
    resize!(coeffs, num_coeffs)
    z = zero(coeffs[1])
    @simd for ord in lencoef+1:num_coeffs
        @inbounds coeffs[ord] = z
    end
    return nothing
end

## Minimum order of an HomogeneousPolynomial compatible with the vector's length
function orderH(coeffs::Array{T,1}) where {T<:Number}
    ord = 0
    ll = length(coeffs)
    for i = 1:get_order()+1
        @inbounds num_coeffs = size_table[i]
        ll ≤ num_coeffs && break
        ord += 1
    end
    return ord
end

## Maximum order of a HomogeneousPolynomial vector; used by TaylorN constructor
function maxorderH(v::Array{HomogeneousPolynomial{T},1}) where {T<:Number}
    m = 0
    @inbounds for i in eachindex(v)
        m = max(m, v[i].order)
    end
    return m
end

## getcoeff ##
"""
    getcoeff(a, n)

Return the coefficient of order `n::Int` of a `a::Taylor1` polynomial.
"""
getcoeff(a::Taylor1, n::Int) = (@assert 0 ≤ n ≤ a.order; return a[n])

getindex(a::Taylor1, n::Int) = a.coeffs[n+1]
getindex(a::Taylor1, u::UnitRange{Int}) = view(a.coeffs, u .+ 1 )
getindex(a::Taylor1, c::Colon) = view(a.coeffs, c)
getindex(a::Taylor1{T}, u::StepRange{Int,Int}) where {T<:Number} =
    view(a.coeffs, u[:] .+ 1)

setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:Number} = a.coeffs[n+1] = x
setindex!(a::Taylor1{T}, x::T, u::UnitRange{Int}) where {T<:Number} =
    a.coeffs[u .+ 1] .= x
function setindex!(a::Taylor1{T}, x::Array{T,1}, u::UnitRange{Int}) where {T<:Number}
    @assert length(u) == length(x)
    for ind in eachindex(x)
        a.coeffs[u[ind]+1] = x[ind]
    end
end
setindex!(a::Taylor1{T}, x::T, c::Colon) where {T<:Number} = a.coeffs[c] .= x
setindex!(a::Taylor1{T}, x::Array{T,1}, c::Colon) where {T<:Number} = a.coeffs[c] .= x
setindex!(a::Taylor1{T}, x::T, u::StepRange{Int,Int}) where {T<:Number} =
    a.coeffs[u[:] .+ 1] .= x
function setindex!(a::Taylor1{T}, x::Array{T,1}, u::StepRange{Int,Int}) where {T<:Number}
    @assert length(u) == length(x)
    for ind in eachindex(x)
        a.coeffs[u[ind]+1] = x[ind]
    end
end

"""
    getcoeff(a, v)

Return the coefficient of `a::HomogeneousPolynomial`, specified by `v`,
which is a tuple (or vector) with the indices of the specific
monomial.
"""
function getcoeff(a::HomogeneousPolynomial, v::NTuple{N,Int}) where {N}
    @assert N == get_numvars() && all(v .>= 0)
    kdic = in_base(get_order(),v)
    @inbounds n = pos_table[a.order+1][kdic]
    a[n]
end
getcoeff(a::HomogeneousPolynomial, v::Array{Int,1}) = getcoeff(a, (v...,))

getindex(a::HomogeneousPolynomial, n::Int) = a.coeffs[n]
getindex(a::HomogeneousPolynomial, n::UnitRange{Int}) = view(a.coeffs, n)
getindex(a::HomogeneousPolynomial, c::Colon) = view(a.coeffs, c)
getindex(a::HomogeneousPolynomial, u::StepRange{Int,Int}) where {T<:Number} =
    view(a.coeffs, u[:])

setindex!(a::HomogeneousPolynomial{T}, x::T, n::Int) where {T<:Number} =
    a.coeffs[n] = x
setindex!(a::HomogeneousPolynomial{T}, x::T, n::UnitRange{Int}) where {T<:Number} =
    a.coeffs[n] .= x
setindex!(a::HomogeneousPolynomial{T}, x::Array{T,1}, n::UnitRange{Int}) where {T<:Number} =
    a.coeffs[n] .= x
setindex!(a::HomogeneousPolynomial{T}, x::T, c::Colon) where {T<:Number} =
    a.coeffs[c] .= x
setindex!(a::HomogeneousPolynomial{T}, x::Array{T,1}, c::Colon) where {T<:Number} =
    a.coeffs[c] = x
setindex!(a::HomogeneousPolynomial{T}, x::T, u::StepRange{Int,Int}) where {T<:Number} =
    a.coeffs[u[:]] .= x
setindex!(a::HomogeneousPolynomial{T}, x::Array{T,1}, u::StepRange{Int,Int}) where {T<:Number} =
    a.coeffs[u[:]] .= x[:]

"""
    getcoeff(a, v)

Return the coefficient of `a::TaylorN`, specified by `v`,
which is a tuple (or vector) with the indices of the specific
monomial.
"""
function getcoeff(a::TaylorN, v::NTuple{N,Int}) where {N}
    order = sum(v)
    @assert order ≤ a.order
    getcoeff(a[order], v)
end
getcoeff(a::TaylorN, v::Array{Int,1}) = getcoeff(a, (v...,))

getindex(a::TaylorN, n::Int) = a.coeffs[n+1]
getindex(a::TaylorN, u::UnitRange{Int}) = view(a.coeffs, u .+ 1)
getindex(a::TaylorN, c::Colon) = view(a.coeffs, c)
getindex(a::TaylorN, u::StepRange{Int,Int}) where {T<:Number} =
    view(a.coeffs, u[:] .+ 1)

function setindex!(a::TaylorN{T}, x::HomogeneousPolynomial{T}, n::Int) where
        {T<:Number}
    @assert x.order == n
    a.coeffs[n+1] = x
end
setindex!(a::TaylorN{T}, x::T, n::Int) where {T<:Number} =
    a.coeffs[n+1] = HomogeneousPolynomial(x, n)
function setindex!(a::TaylorN{T}, x::T, u::UnitRange{Int}) where {T<:Number}
    for ind in u
        a[ind] = x
    end
    a[u]
end
function setindex!(a::TaylorN{T}, x::Array{HomogeneousPolynomial{T},1}, u::UnitRange{Int}) where {T<:Number}
    @assert length(u) == length(x)
    for ind in eachindex(x)
        a[u[ind]] = x[ind]
    end
end
function setindex!(a::TaylorN{T}, x::Array{T,1}, u::UnitRange{Int}) where {T<:Number}
    @assert length(u) == length(x)
    for ind in eachindex(x)
        a[u[ind]] = x[ind]
    end
end
setindex!(a::TaylorN{T}, x::T, ::Colon) where {T<:Number} =
    (a[0:end] = x; a[:])
setindex!(a::TaylorN{T}, x::Array{HomogeneousPolynomial{T},1}, ::Colon) where
    {T<:Number} = (a[0:end] = x; a[:])
setindex!(a::TaylorN{T}, x::Array{T,1}, ::Colon) where {T<:Number} =
    (a[0:end] = x; a[:])
function setindex!(a::TaylorN{T}, x::T, u::StepRange{Int,Int}) where {T<:Number}
    for ind in u
        a[ind] = x
    end
    a[u]
end
function setindex!(a::TaylorN{T}, x::Array{HomogeneousPolynomial{T},1}, u::StepRange{Int,Int}) where {T<:Number}
    # a[u[:]] .= x[:]
    @assert length(u) == length(x)
    for ind in eachindex(x)
        a[u[ind]] = x[ind]
    end
end
function setindex!(a::TaylorN{T}, x::Array{T,1}, u::StepRange{Int,Int}) where {T<:Number}
    @assert length(u) == length(x)
    for ind in eachindex(x)
        a[u[ind]] = x[ind]
    end
end

## eltype, length, get_order, etc ##
for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
    @eval begin
        if $T == HomogeneousPolynomial
            @inline iterate(a::$T, state=1) = state > length(a) ? nothing : (a.coeffs[state], state+1)
            # Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1)
            @inline length(a::$T) = size_table[a.order+1]
            @inline firstindex(a::$T) = 1
            @inline lastindex(a::$T) = length(a)
        else
            @inline iterate(a::$T, state=0) = state > a.order ? nothing : (a.coeffs[state+1], state+1)
            # Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1)
            @inline length(a::$T) = length(a.coeffs)
            @inline firstindex(a::$T) = 0
            @inline lastindex(a::$T) = a.order
        end
        @inline eachindex(a::$T) = firstindex(a):lastindex(a)
        @inline eltype(::$T{S}) where {S<:Number} = S
        @inline size(a::$T) = size(a.coeffs)
        @inline get_order(a::$T) = a.order
        @inline axes(a::$T) = ()
    end
end

## fixorder ##
for T in (:Taylor1, :TaylorN)
    @eval begin
        @inline function fixorder(a::$T, b::$T)
            a.order == b.order && return a, b
            minorder, maxorder = minmax(a.order, b.order)
            if minorder ≤ 0
                minorder = maxorder
            end
            return $T(copy(a.coeffs), minorder), $T(copy(b.coeffs), minorder)
        end
    end
end

function fixorder(a::HomogeneousPolynomial, b::HomogeneousPolynomial)
    @assert a.order == b.order
    return a, b
end

# Finds the first non zero entry; extended to Taylor1
function Base.findfirst(a::Taylor1{T}) where {T<:Number}
    first = findfirst(x->!iszero(x), a.coeffs)
    isa(first, Nothing) && return -1
    return first-1
end
# Finds the last non-zero entry; extended to Taylor1
function Base.findlast(a::Taylor1{T}) where {T<:Number}
    last = findlast(x->!iszero(x), a.coeffs)
    isa(last, Nothing) && return -1
    return last-1
end

## copyto! ##
# Inspired from base/abstractarray.jl, line 665
for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
    @eval function copyto!(dst::$T{T}, src::$T{T}) where {T}
        length(dst) < length(src) && throw(ArgumentError(string("Destination has fewer elements than required; no copy performed")))
        destiter = eachindex(dst)
        y = iterate(destiter)
        for x in src
            dst[y[1]] = x
            y = iterate(destiter, y[2])
        end
        return dst
    end
end

"""
    constant_term(a)

Return the constant value (zero order coefficient) for `Taylor1`
and `TaylorN`. The fallback behavior is to return `a` itself if
`a::Number`, or `a[1]` when `a::Vector`.
"""
constant_term(a::Taylor1) = a[0]

constant_term(a::TaylorN) = a[0][1]

constant_term(a::Vector{T}) where {T<:Number} = constant_term.(a)

constant_term(a::Number) = a

"""
    linear_polynomial(a)

Return the linear part of `a` as a polynomial (`Taylor1`
or `TaylorN`). The fallback behavior is to return `a` itself if
`a::Number`, or `a[1]` when `a::Vector`.
"""
linear_polynomial(a::Taylor1) = Taylor1([zero(a[1]), a[1]])

linear_polynomial(a::TaylorN) = TaylorN([a[1]])

linear_polynomial(a::Vector{T}) where {T<:Number} = linear_polynomial.(a)

linear_polynomial(a::Number) = a

#=
Tuple operations taken from ForwardDiff.jl

Copyright (c) 2015: Jarrett Revels, Theodore Papamarkou, Miles Lubin, and other contributors

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

=#
function tupexpr(f, N)
    ex = Expr(:tuple, [f(i) for i=1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $ex
    end
end

@generated function scale_tuple(tup::NTuple{N}, x) where N
    return tupexpr(i -> :(tup[$i] * x), N)
end

@generated function div_tuple_by_scalar(tup::NTuple{N}, x) where N
    return tupexpr(i -> :(tup[$i] / x), N)
end

@generated function add_tuples(a::NTuple{N}, b::NTuple{N})  where N
    return tupexpr(i -> :(a[$i] + b[$i]), N)
end

@generated function sub_tuples(a::NTuple{N}, b::NTuple{N})  where N
    return tupexpr(i -> :(a[$i] - b[$i]), N)
end

@generated function minus_tuple(tup::NTuple{N}) where N
    return tupexpr(i -> :(-tup[$i]), N)
end

@generated function mul_tuples(a::NTuple{N}, b::NTuple{N}, afactor, bfactor) where N
    return tupexpr(i -> :((afactor * a[$i]) + (bfactor * b[$i])), N)
end


text:
## -----------  struct2dict

"""
	$(SIGNATURES)

Struct that contains information about one node in the `Dict` created by `struct2dict`.
"""
struct NodeInfo{T}
    # type :: DataType
    size
    value :: T
end

NodeInfo(s :: Number) = NodeInfo{typeof(s)}(nothing, s);
NodeInfo(s :: AbstractString) = NodeInfo{String}(nothing, String(s));
NodeInfo(s :: Symbol) = NodeInfo{Symbol}(nothing, s);
# NodeInfo(s :: AbstractArray{F1}) where F1 <: AbstractFloat =
#     NodeInfo{typeof(s)}(size(s), s);

node_value(n :: NodeInfo{T}) where T = n.value;
node_type(n :: NodeInfo{T}) where T = T;
node_size(n :: NodeInfo{T}) where T = n.size;

Base.isequal(n1 :: NodeInfo{T}, n2 :: NodeInfo{T}) where T =
    isequal(node_size(n1), node_size(n2))  &&
    isequal(node_value(n1), node_value(n2));
Base.isequal(n1 :: NodeInfo{T1}, n2 :: NodeInfo{T2}) where {T1, T2} = false;

"""
	$(SIGNATURES)

Convert a `struct` to a nested `Dict`. Known leaf types, such as strings, arrays become `NodeInfo` entries in the `Dict`. `struct`s become their own `Dict`s.

The purpose is to convert a `struct` into a `Dict` of simple "built-in" types that can be serialized. Each leaf entry contains information about the `DataType` of the leaf, so that the object can be reconstructed. In case the data are saved to a format that does not preserve size information, the size of the object is also saved.
"""
function struct2dict(s)
    d = Dict{Symbol, Any}();
    # d[:ownType] = typeof(s);
    for pName in propertynames(s)
        d[pName] = struct2dict(getproperty(s, pName));
    end
    return d
end

# For each known type, `struct2dict` produces a `NodeInfo` object.
# The user can extend with new methods.
scalar2dict(s) = NodeInfo{typeof(s)}(nothing, s);
struct2dict(s :: Number) = scalar2dict(s);
struct2dict(s :: AbstractString) = scalar2dict(s);
struct2dict(s :: Symbol) = scalar2dict(s);
struct2dict(s :: Dict) = scalar2dict(s);

array2dict(s :: AbstractArray{T}) where T = NodeInfo{typeof(s)}(size(s), s);
struct2dict(s :: AbstractArray{F1}) where F1 <: Number = array2dict(s);
struct2dict(s :: AbstractArray{F1}) where F1 <: AbstractString = array2dict(s);
struct2dict(s :: AbstractArray{F1}) where F1 <: Symbol = array2dict(s);

"""
	$(SIGNATURES)

The inverse of `struct2dict`. Takes a fully initialized MUTABLE `struct` as input. Copies the values from the `Dict` into the `struct`.
Mismatches between `Dict` keys and `propertynames(s)` are ignored. That means: fields that are missing from the `Dict` are not initialized. Extra fields in the `Dict` are ignored.
"""
function dict2struct!(s, d :: Dict)
    # errorCnt = 0;
    @assert !isimmutable(s)  "Not mutable: $(typeof(s))"
    for pName in propertynames(s)
        if haskey(d, pName)
            if isa(d[pName], Dict)
                dict2struct!(getproperty(s, pName), d[pName]);
            else
                node2field!(s, pName, d[pName]);
            end
        end
    end
    return nothing
end

"""
	$(SIGNATURES)

Set one field in a `struct` from a `NodeInfo`. Catch errors and return `success` as `Bool`.
User can override for specific `NodeInfo` types `T`.
"""
function node2field!(s, pName :: Symbol, n :: NodeInfo{T}) where T
    success = true;
    try
        setproperty!(s, pName, node_value(n));
    catch
        success = false;
        nodeType = node_type(n);
        @warn """
            Cannot set property $pName of $(typeof(s))
              Node of type $nodeType
            """
    end
    return success
end

# ---------------


text:
__precompile__()

module LohmannRouting

include("routing.jl")

export upstream_catchments,
       flow_sequence
       calc_xmask,
       findnearest,
       upsample,
       match_dates,
       make_irf,
       make_river_irf,
       make_grid_uh,
       rout

end # module


text:
function en(an,anm1,Ln)
	return sqrt((1-abs2(an))/((one(an)-abs2(anm1))*(one(an)-abs2(Ln))))
end

function ln(z::Vector{ComplexPair{BigFloat}},a::Vector{ComplexPair{BigFloat}},phinm1::Vector{ComplexPair{BigFloat}})
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(ComplexPair{BigFloat})
        phinm1star = phinm1
        f1 = [((z[i]-anm1)/(one(ComplexPair{BigFloat})-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
        f2 = [((one(ComplexPair{BigFloat})-conj(anm1)*z[i])/(one(ComplexPair{BigFloat})-conj(an)*z[i]))*phinm1star[i] for i=1:length(z)]
        num = sum(f1)
        den = sum(f2)
        r = -num/den
    else
        an = copy(a[end])
        anm1 = copy(a[end-1])
        bProd = length(a)>=3 ? conjBlaschkeProd(z,[a[1:end-2]]) : ones(ComplexPair{BigFloat},length(z))
        f1 = [((z[i]-anm1)/(one(ComplexPair{BigFloat})-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
        f2 = [((z[i]-anm1)/(one(ComplexPair{BigFloat})-conj(an)*z[i]))*phinm1[i]*bProd[i] for i=1:length(z)]
        num = sum(f1)
        den = sum(f2)
        r = -num/den
    end
    return r
end

function phiNext(z::Vector{ComplexPair{BigFloat}},a::Vector{ComplexPair{BigFloat}},phinm1::Vector{ComplexPair{BigFloat}})
    phin = zeros(ComplexPair{BigFloat},length(phinm1))
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(ComplexPair{BigFloat})
        L = ln(z,a,phinm1)
        c = en(an,anm1,L)
        for k=1:length(z)
            q = (z[k]-anm1)
            p = (one(ComplexPair{BigFloat})-conj(anm1)*z[k])
            C = c/(z[k]-an)
            PHI = p*phinm1[k]
            PHISTAR = q*phinm1[k]
            phin[k] = C*(PHI+conj(L)*PHISTAR)
        end
    else
        an = a[end]
        anm1 = a[end-1]
        L = ln(z,a,phinm1)
        c = en(an,anm1,L)
        bProd = length(a)>=3 ? blaschkeProd(z,[a[1:end-2]]) : ones(ComplexPair{BigFloat},length(z))
        for k=1:length(z)
            q = (z[k]-an)
            p = (one(ComplexPair{BigFloat})-conj(anm1)*z[k])
            C = c*p/q
            phin[k] = C*(phinm1[k]+conj(L)*bProd[k]*conj(phinm1[k]))
        end
    end
    return (phin, L)
end

function phiNext(z::Vector{ComplexPair{BigFloat}},a::Vector{ComplexPair{BigFloat}},phinm1::Vector{ComplexPair{BigFloat}},L::ComplexPair{BigFloat})
    phin = zeros(ComplexPair{BigFloat},length(phinm1))
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(ComplexPair{BigFloat})
        c = en(an,anm1,L)
        for k=1:length(z)
            q = (z[k]-anm1)
            p = (one(ComplexPair{BigFloat})-conj(anm1)*z[k])
            C = c/(z[k]-an)
            PHI = p*phinm1[k]
            PHISTAR = q*phinm1[k]
            phin[k] = C*(PHI+conj(L)*PHISTAR)
        end
    else
        an = a[end]
        anm1 = a[end-1]
        c = en(an,anm1,L)
        bProd = length(a)>=3 ? blaschkeProd(z,[a[1:end-2]]) : ones(ComplexPair{BigFloat},length(z))
        for k=1:length(z)
            q = (z[k]-an)
            p = (one(ComplexPair{BigFloat})-conj(anm1)*z[k])
            C = c*p/q
            phin[k] = C*(phinm1[k]+conj(L)*bProd[k]*conj(phinm1[k]))
        end
    end
    return phin
end

function blaschkeProd(z::Vector{ComplexPair{BigFloat}},aj::Vector{ComplexPair{BigFloat}})
    r = ones(ComplexPair{BigFloat},length(z))
    for i=1:length(aj)
        t = [(one(ComplexPair{BigFloat})-conj(aj[i])*z[k])./(z[k]-aj[i]) for k=1:length(z)]
        r = r.*t
    end
    return r
end

function conjBlaschkeProd(z::Vector{ComplexPair{BigFloat}},aj::Vector{ComplexPair{BigFloat}})
    r = ones(ComplexPair{BigFloat},length(z))
    for i=1:length(aj)
        t = [(z[k]-aj[i])/(one(ComplexPair{BigFloat})-conj(aj[i])*z[k]) for k=1:length(z)]
        r = r.*t
    end
    return r
end

function genZArb(N::Int)
	return exp(convert(Vector{ComplexPair{BigFloat}},im*linspace2(0.0,2*pi,N)))
end

function genZArcArb(N::Int,Div::Int)
    return exp(convert(Vector{ComplexPair{BigFloat}},im*linspace(0.0,2*pi/Div,N)))
end

function genZArbUpSample(N::Int,Mult::Int)
    return exp(convert(Vector{ComplexPair{BigFloat}},im*linspace2upsample(0.0,2*pi,N,Mult)))
end

function genBasis(z::Vector,a::Vector{ComplexPair{BigFloat}})
    P = zeros(ComplexPair{BigFloat},(length(z),length(a)+1))
    P[:,1] = ones(ComplexPair{BigFloat},length(z))/sqrt(length(z));
    L = zeros(ComplexPair{BigFloat},length(a))
    for k=1:length(a)
        (P[:,k+1],L[k]) = phiNext(z,[a[1:k]],P[:,k])
    end
    return (P,L)
end

function genBasisGivenL(z::Vector{ComplexPair{BigFloat}},a::Vector{ComplexPair{BigFloat}},L::Vector{ComplexPair{BigFloat}},N::Int)
    P = zeros(ComplexPair{BigFloat},(length(z),length(a)+1))
    P[:,1] = ones(ComplexPair{BigFloat},length(z))/sqrt(N);
    for k=1:length(a)
        P[:,k+1] = phiNext(z,[a[1:k]],P[:,k],L[k])
    end
    return P
end

function unitDiskGrid(n::Int,delta::Float64)
    N = nextpow2(n)
    M = N>>1
    v = linspace(delta,one(Float64)-delta,M)
    Zjk = zeros(ComplexPair{BigFloat},N,N)
    for j=1:(M)
            #bottom, running counter-clockwise
            Zjk[(M+j),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(5*pi/4,7*pi/4,2*j))
            #top, running clockwise
            Zjk[(M-j+1),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(3*pi/4,1*pi/4,2*j))
    end
    for j=1:(M)
            #right, running clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M+j] = v[j]*exp(im*linspace(pi/4,-pi/4,2*j)) #right
            #left, running counter-clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M-j+1] = v[j]*exp(im*linspace(3*pi/4,5*pi/4,2*j)) #left
    end
    return (reshape(Zjk,N^2))
end


text:
# Tags:
Cassette.@context HessianSparsityContext

const TaggedOf{T} = Tagged{A, T} where A

const HTagType = Union{Input, TermCombination}
Cassette.metadatatype(::Type{<:HessianSparsityContext}, ::DataType) = HTagType

istainted(ctx::HessianSparsityContext, x) = ismetatype(x, ctx, TermCombination)

@inline Cassette.overdub(ctx::HessianSparsityContext, f::typeof(istainted), x) =
    istainted(ctx, x)
@inline Cassette.overdub(ctx::HessianSparsityContext, f::typeof(this_here_predicate!)) =
    this_here_predicate!(ctx.metadata)

# getindex on the input
@inline function Cassette.overdub(ctx::HessianSparsityContext,
                                  f::typeof(getindex),
                                  X::Tagged,
                                  idx::Tagged...)
    if any(i->ismetatype(i, ctx, TermCombination) && !isone(metadata(i, ctx)), idx)
        error("getindex call depends on input. Cannot determine Hessian sparsity")
    end
    Cassette.overdub(ctx, f, X, map(i->untag(i, ctx), idx)...)
end

# plugs an ambiguity
@inline function Cassette.overdub(ctx::HessianSparsityContext,
                                  f::typeof(getindex),
                                  X::Tagged)
    Cassette.recurse(ctx, f, X)
end

@inline function Cassette.overdub(ctx::HessianSparsityContext,
                                  f::typeof(getindex),
                                  X::Tagged,
                                  idx::Integer...)
    if ismetatype(X, ctx, Input)
        val = Cassette.fallback(ctx, f, X, idx...)
        i = LinearIndices(untag(X, ctx))[idx...]
        tag(val, ctx, TermCombination(Set([Dict(i=>1)])))
    else
        Cassette.recurse(ctx, f, X, idx...)
    end
end

@inline function Cassette.overdub(ctx::HessianSparsityContext,
                                  f::typeof(Base.unsafe_copyto!),
                                  X::Tagged,
                                  xstart,
                                  Y::Tagged,
                                  ystart,
                                  len)
    if ismetatype(Y, ctx, Input)
        val = Cassette.fallback(ctx, f, X, xstart, Y, ystart, len)
        nometa = Cassette.NoMetaMeta()
        X.meta.meta[xstart:xstart+len-1] .= (i->Cassette.Meta(TermCombination(Set([Dict(i=>1)])), nometa)).(ystart:ystart+len-1)
        val
    else
        Cassette.recurse(ctx, f, X, xstart, Y, ystart, len)
    end
end

@inline function Cassette.overdub(ctx::HessianSparsityContext,
                                  f::typeof(copy),
                                  X::Tagged)
    if ismetatype(X, ctx, Input)
        val = Cassette.fallback(ctx, f, X)
        tag(val, ctx, Input())
    else
        Cassette.recurse(ctx, f, X)
    end
end

combine_terms(::Nothing, terms...) = one(TermCombination)

# 1-arg functions
combine_terms(::Val{true}, term) = term
combine_terms(::Val{false}, term) = term * term

# 2-arg functions
function combine_terms(::Val{linearity}, term1, term2) where linearity

    linear11, linear22, linear12 = linearity
    term = zero(TermCombination)
    if linear11
        if !linear12
            term += term1
        end
    else
        term += term1 * term1
    end

    if linear22
        if !linear12
            term += term2
        end
    else
        term += term2 * term2
    end

    if linear12
        term += term1 + term2
    else
        term += term1 * term2
    end
    term
end

# Hessian overdub
#
function getterms(ctx, x)
    ismetatype(x, ctx, TermCombination) ? metadata(x, ctx) : one(TermCombination)
end

function hessian_overdub(ctx::HessianSparsityContext, f, linearity, args...)
    t = combine_terms(linearity, map(x->getterms(ctx, x), args)...)
    val = Cassette.fallback(ctx, f, args...)
    tag(val, ctx, t)
end

@inline function Cassette.overdub(ctx::HessianSparsityContext,
                                  f::typeof(getproperty),
                                  x::Tagged, prop)
    if ismetatype(x, ctx, TermCombination) && !isone(metadata(x, ctx))
        error("property of a non-constant term accessed")
    else
        Cassette.recurse(ctx, f, x, prop)
    end
end

haslinearity(ctx::HessianSparsityContext, f, nargs) = haslinearity(untag(f, ctx), nargs)
linearity(ctx::HessianSparsityContext, f, nargs) = linearity(untag(f, ctx), nargs)

@inline function Cassette.overdub(ctx::HessianSparsityContext,
                                  f,
                                  args...)
    tainted = any(x->ismetatype(x, ctx, TermCombination), args)
    val = if tainted && haslinearity(ctx, f, Val{nfields(args)}())
        l = linearity(ctx, f, Val{nfields(args)}())
        hessian_overdub(ctx, f, l, args...)
    else
        val = Cassette.recurse(ctx, f, args...)
       #if tainted && !ismetatype(val, ctx, TermCombination)
       #    @warn("Don't know the linearity of function $f")
       #end
        val
    end
    val
end


text:
# Tests automatic differentiation on sn.jl:
include("../src/sn_jacobian.jl")
using PyPlot

@testset "sn_jacobian" begin

function test_sn_jacobian2(l_max)
r0 = [0.01,100.0]
nb = 10
#l_max = 20
#l_max = 5
n_max = l_max^2+2*l_max

# Now, carry out finite-difference derivative computation:
function sn_jac_num(l_max::Int64,r::T,b::T) where {T <: Real}
  dq = big(1e-24)
  n_max = l_max^2+2*l_max
# Allocate an array for s_n:
  sn_big = zeros(BigFloat,n_max+1)
# Make BigFloat versions of r & b:
  r_big = big(r); b_big = big(b)
# Compute s_n to BigFloat precision:
  s_n!(l_max,r_big,b_big,sn_big)
# Now, compute finite differences:
  sn_jac_big= zeros(BigFloat,n_max+1,2)
  sn_plus = zeros(BigFloat,n_max+1)
  s_n!(l_max,r_big+dq,b_big,sn_plus)
  sn_minus = zeros(BigFloat,n_max+1)
  s_n!(l_max,r_big-dq,b_big,sn_minus)
  sn_jac_big[:,1] = (sn_plus-sn_minus)*.5/dq
  s_n!(l_max,r_big,b_big+dq,sn_plus)
  s_n!(l_max,r_big,b_big-dq,sn_minus)
  sn_jac_big[:,2] = (sn_plus-sn_minus)*.5/dq
return convert(Array{Float64,1},sn_big),convert(Array{Float64,2},sn_jac_big)
end

#sn_jacobian
#convert(Array{Float64,2},sn_jac_big)-sn_jacobian

get_cmap("plasma")
epsilon = 1e-10; delta = 1e-3
for i=1:2
  fig,axes = subplots(1,1)
  r=r0[i]
  if r < 1.0
    b = [linearspace(epsilon,delta,nb); linearspace(delta,r-delta,nb);
     -logarithmspace(log10(delta),log10(epsilon),nb) .+r; linearspace(r-epsilon,r+epsilon,nb); logarithmspace(log10(epsilon),log10(delta),nb) .+r;
     linearspace(r+delta,1-r-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+(1-r); linearspace(1-r-epsilon,1-r+epsilon,nb);
     logarithmspace(log10(epsilon),log10(delta),nb).+(1-r); linearspace(1-r+delta,1+r-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+(1+r)]
     nticks = 14
     xticknames=[L"$10^{-15}$",L"$10^{-12}$",L"$10^{-3}$",L"$r-10^{-3}$",L"$r-10^{-12}$",L"$r+10^{-12}$",L"$r+10^{-3}$",
     L"$1-r-10^{-3}$",L"$1-r-10^{-12}$",L"$1-r+10^{-12}$",L"$1-r+10^{-3}$",L"$1+r-10^{-3}$",L"$1+r-10^{-12}$",L"$1+r-10^{-15}$"]
  else
    b = [logarithmspace(log10(epsilon),log10(delta),nb).+(r-1); linearspace(r-1+delta,r-delta,nb);
     -logarithmspace(log10(delta),log10(epsilon),nb) .+r; linearspace(r-epsilon,r+epsilon,nb); logarithmspace(log10(epsilon),log10(delta),nb) .+r;
     linearspace(r+delta,r+1-delta,nb); -logarithmspace(log10(delta),log10(epsilon),nb) .+(r+1)]
     nticks = 8
     xticknames=[L"$r-1+10^{-12}$",L"$r-1+10^{-3}$",L"$r-10^{-3}$",L"$r-10^{-12}$",L"$r+10^{-12}$",L"$r+10^{-3}$",
     L"$r+1-10^{-3}$",L"$r+1-10^{-12}$"]
  end
#  if r < 1.0
#    b = [linearspace(1e-8,epsilon,nb); linearspace(epsilon,delta,nb); linearspace(delta,r-delta,nb);
##     linearspace(r-delta,r-epsilon,nb); linearspace(r-epsilon,r,nb); linearspace(r,r+epsilon,nb); linearspace(r+epsilon,r+delta,nb);
#     linearspace(r-delta,r-epsilon,nb); linearspace(r-epsilon,r+epsilon,nb); linearspace(r+epsilon,r+delta,nb);
#     linearspace(r+delta,1-r-delta,nb); linearspace(1-r-delta,1-r-epsilon,nb); linearspace(1-r-epsilon,1-r+epsilon,nb);
#     linearspace(1-r+epsilon,1-r+delta,nb); linearspace(1-r+delta,1+r-delta,nb); linearspace(1+r-delta,1+r-epsilon,nb);linearspace(1+r-epsilon,1+r-1e-8,nb)]
#  else
#    b = [linearspace(r-1+1e-8,r-1+epsilon,nb); linearspace(r-1+epsilon,r-1+delta,nb); linearspace(r-1+delta,r-delta,nb);
##     linearspace(r-delta,r-epsilon,nb); linearspace(r-epsilon,r,nb); linearspace(r,r+epsilon,nb); linearspace(r+epsilon,r+delta,nb);
#     linearspace(r-delta,r-epsilon,nb); linearspace(r-epsilon,r+epsilon,nb); linearspace(r+epsilon,r+delta,nb);
#     linearspace(r+delta,r+1-delta,nb); linearspace(r+1-delta,r+1-epsilon,nb); linearspace(r+1-epsilon,r+1-1e-8,nb)]
#  end
  k2 = (1.-(r-b).^2)./(4.*b.*r)
  k2H = ((b+1).^2-r.^2)./(4b)
  sn_jac_grid = zeros(length(b),n_max+1,2)
  sn_jac_array= zeros(n_max+1,2)
  sn_grid = zeros(length(b),n_max+1)
  sn_grid_big = zeros(length(b),n_max+1)
  sn_array= zeros(n_max+1)
  sn_jac_grid_num = zeros(length(b),n_max+1,2)
  for j=1:length(b)
#    println("r: ",r," b: ",b[j])
    sn_array,sn_jac_array= sn_jac(l_max,r,b[j])
    s_n!(l_max,r,b[j],sn_array)
    sn_grid[j,:]=sn_array
    sn_jac_grid[j,:,:]=sn_jac_array
    sn_array,sn_jac_array = sn_jac_num(l_max,r,b[j])
    sn_jac_grid_num[j,:,:]=sn_jac_array
    sn_grid_big[j,:]=sn_array
  end
# Now, make plots:
  ax = axes
  m=0;l=0
  for n=0:n_max
    if m==0
#      ax[:plot](abs.(sn_jac_grid[:,n+1,1]-sn_jac_grid_num[:,n+1,1]),color=cmap(l/(l_max+2)),lw=1,label=string("l=",l))
#      ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,1])-sn_jac_grid_num[:,n+1,1]),lw=1,label=string("l=",l))
#      ax[:semilogy](b,abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1,label=string("l=",l))
#      ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1,label=string("l=",l))
      ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1,label=string("l=",l))
#      semilogy(abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1,label=string("l=",l))
    else
#      ax[:plot](abs.(sn_jac_grid[:,n+1,1]-sn_jac_grid_num[:,n+1,1]),color=cmap(l/(l_max+2)),lw=1)
#      ax[:semilogy](abs.(sn_jac_grid[:,n+1,1]-sn_jac_grid_num[:,n+1,1]),lw=1)
      ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1)
#      semilogy(abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1)
#      ax[:semilogy](b,abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1)
    end
#    ax[:plot](abs.(sn_jac_grid[:,n+1,2]-sn_jac_grid_num[:,n+1,2]),color=cmap(l/(l_max+2)),lw=1)
#    ax[:semilogy](abs.(sn_jac_grid[:,n+1,2]-sn_jac_grid_num[:,n+1,2]),lw=1)
#    ax[:semilogy](b,abs.(asinh.(sn_jac_grid[:,n+1,2])-asinh.(sn_jac_grid_num[:,n+1,2])),lw=1)
    ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,2])-asinh.(sn_jac_grid_num[:,n+1,2])),lw=1)
    ax[:semilogy](abs.(asinh.(sn_grid[:,n+1])-asinh.(sn_grid_big[:,n+1])),lw=1)
#    semilogy(abs.(asinh.(sn_jac_grid[:,n+1,2])-asinh.(sn_jac_grid_num[:,n+1,2])),lw=1)
    println("n: ",n," m: ",m," l: ",l," mu: ",l-m," nu: ",l+m," max dsn/dr: ",maximum(abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1]))),
      " maximum dsn/db: ",maximum(abs.(asinh.(sn_jac_grid[:,n+1,2])-asinh.(sn_jac_grid_num[:,n+1,2]))))
    m +=1
    if m > l
      l += 1
      m = -l
    end
  end
  ax[:legend](loc="upper right",fontsize=6)
  ax[:set_xlabel]("b values")
  ax[:set_ylabel]("Derivative Error")
  ax[:axis]([0,length(b),1e-16,1])
  ax[:set_xticks](nb*linearspace(0,nticks-1,nticks))
  ax[:set_xticklabels](xticknames,rotation=45)
  if i==1
    ax[:set_title]("r = 0.01")
  else
    ax[:set_title]("r = 100")
  end
#  ax[:axis]([minimum(b),maximum(b),1e-16,1])
#  ax[:axis]([0,length(b),1e-16,1])
#  read(STDIN,Char)
  clf()
  for n=0:n_max
      plot(b,sn_grid[:,n+1])
  end
#  read(STDIN,Char)
## Loop over n and see where the differences between the finite-difference
## and AutoDiff are greater than the derivative value:
l=0; m=0
for n=0:n_max
  diff = sn_jac_grid[:,n+1,:]-sn_jac_grid_num[:,n+1,:]
  for ib=1:length(b)
    test1 = isapprox(sn_jac_grid[ib,n+1,1],sn_jac_grid_num[ib,n+1,1],atol=1e-8)
    @test test1
    if ~test1
      println("b: ",b[ib]," n: ",n," dsn/dr: ",sn_jac_grid[ib,n+1,1]," dsn/dr_num: ",sn_jac_grid_num[ib,n+1,1]," diff: ",sn_jac_grid[ib,n+1,1]-sn_jac_grid_num[ib,n+1,1])
    end
    test2 = isapprox(sn_jac_grid[ib,n+1,2],sn_jac_grid_num[ib,n+1,2],atol=1e-8)
    @test test2
    if ~test2
      println("b: ",b[ib]," n: ",n," dsn/db: ",sn_jac_grid[ib,n+1,2]," dsn/db_num: ",sn_jac_grid_num[ib,n+1,2]," diff: ",sn_jac_grid[ib,n+1,2]-sn_jac_grid_num[ib,n+1,2])
    end
  end
  mask = (abs.(diff) .> 1e-3*abs.(sn_jac_grid[:,n+1,:])) .& (abs.(sn_jac_grid[:,n+1,:]) .> 1e-5)
  if sum(mask) > 0 || mod(l-m,4) == 0
    println("n: ",n," max dsn/dr: ",maximum(abs.(diff)))
    println("b: ",b[mask[:,1]],b[mask[:,2]]," k^2_H: ",k2H[mask[:,1]],k2H[mask[:,2]])
# end
#end
# Now, plot derivatives for each value of n:
#for n=0:n_max
    clf()
    plot(b,sn_grid[:,n+1])
    plot(b,sn_jac_grid[:,n+1,1])
    plot(b[mask[:,1]],sn_jac_grid[mask[:,1],n+1,1],"o")
    plot(b,sn_jac_grid[:,n+1,2])
    plot(b[mask[:,2]],sn_jac_grid[mask[:,2],n+1,2],"o")
    plot(b,sn_jac_grid_num[:,n+1,1],linestyle="--",linewidth=2)
    plot(b,sn_jac_grid_num[:,n+1,2],linestyle="--",linewidth=2)
    println("n: ",n," l: ",l," m: ",m," mu: ",l-m," nu: ",l+m)
#    read(STDIN,Char)
  end
  m +=1
  if m > l
    l += 1
    m = -l
  end
end
#
#
##end
#
##loglog(abs.(reshape(sn_jac_grid,length(b)*(n_max+1)*2)),abs.(reshape(sn_jac_grid-sn_jac_grid_num,length(b)*(n_max+1)*2)),".")
##loglog(abs.(reshape(sn_jac_grid,length(b)*(n_max+1)*2)),abs.(reshape(sn_jac_grid_num,length(b)*(n_max+1)*2)),".")
end
return
end

test_sn_jacobian2(10)
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using ThArrays
using Test

@testset "Tensor Indexing" begin

    j_vect = rand(10)
    t_vect = Tensor(j_vect)
    j_data = rand(2, 3, 4)
    t_data = Tensor(j_data)

    @testset "Indexing with Int" begin
        for i in 1:length(j_vect)
            @test t_vect[i][] == j_vect[i]
        end
        for i in 1:length(j_data)
            @test t_data[i][] == j_data[i]
        end
    end

    @testset "Indexing with Range" begin
        @test t_data[1:6] == Tensor(j_data[1:6])
        @test t_data[8:12] == Tensor(j_data[8:12])
        @test t_data[18:24] == Tensor(j_data[18:24])
    end

    @testset "Indexing with CartesianIndex" begin
        @test t_data[1, 2, :] == Tensor(j_data[1, 2, :])
        @test t_data[[1], 2, :] == Tensor(j_data[[1], 2, :])
        @test t_data[:, 2, :] == Tensor(j_data[:, 2, :])
        @test t_data[:, [2], :] == Tensor(j_data[:, [2], :])
        @test t_data[1:2, 2:3, 2:4] == Tensor(j_data[1:2, 2:3, 2:4])
    end

    @testset "Set Index with Int" begin
        for i in 1:length(j_data)
            t_data[i] = i
        end
        for i in 1:length(j_data)
            @test t_data[i][] == i
        end
    end

    @testset "Set Index with Range" begin
        t_data[1:4] = collect(1:4)
        for i in 1:4
            @test t_data[i][] == i
        end
    end

    @testset "Set Index with CartesianIndex" begin
        t_data[1, 1:2, 1:2] = zeros(2, 2)
        @test t_data[1, 1:2, 1:2] == Tensor(zeros(2, 2))
    end
end


text:
"""
set_pvals!(model::CS.Optimizer)

Set the possible values for each constraint.
"""
function set_pvals!(model::CS.Optimizer)
    com = model.inner
    for constraint in com.constraints
        set_pvals!(com, constraint)
    end
end

"""
    set_var_in_all_different!(model::CS.Optimizer)

Set the vector `in_all_different` for each variable
"""
function set_var_in_all_different!(model::CS.Optimizer)
    com = model.inner
    n_all_different = com.info.n_constraint_types.alldifferent
    for var in com.search_space
        var.in_all_different = zeros(Bool, n_all_different)
    end
    ci = 1
    for constraint in com.constraints
        if constraint isa AllDifferentConstraint
            for vidx in constraint.indices
                com.search_space[vidx].in_all_different[ci] = true
            end
            ci += 1
        end
    end
end

"""
init_constraints!(com::CS.CoM; constraints=com.constraints)

Initializes all `constraints` which implement the `init_constraint!` method.
It also activates all constraints which implement the `activate_constraint!` method.
Return if feasible after initalization
"""
function init_constraints!(com::CS.CoM; constraints = com.constraints)
    feasible = true
    for constraint in constraints
        constraint.is_deactivated && continue
        if constraint.impl.init
            feasible = init_constraint!(com, constraint, constraint.fct, constraint.set)
            !feasible && break
        end
        constraint.is_initialized = true
        if constraint.impl.activate
            feasible = activate_constraint!(com, constraint, constraint.fct, constraint.set)
            !feasible && break
        end
        constraint.is_activated = true
    end
    return feasible
end

"""
init_constraints!(com::CS.CoM; constraints=com.constraints)

Initializes all constraints of the model as new `constraints` were added.
Return if feasible after the update of the initalization
"""
function update_init_constraints!(com::CS.CoM; constraints = com.constraints)
    feasible = true
    for constraint in com.constraints
        if constraint.impl.update_init
            feasible = update_init_constraint!(
                com,
                constraint,
                constraint.fct,
                constraint.set,
                constraints,
            )
            !feasible && break
        end
    end
    return feasible
end

"""
set_impl_functions!(com, constraint::Constraint)

Set std.impl.[] for each constraint
"""
function set_impl_functions!(com, constraint::Constraint)
    if com.sense != MOI.FEASIBILITY_SENSE
        set_impl_update_best_bound!(constraint)
    end
    set_impl_init!(constraint)
    set_impl_activate!(constraint)
    set_impl_update_init!(constraint)
    set_impl_finished_pruning!(constraint)
    set_impl_restore_pruning!(constraint)
    set_impl_reverse_pruning!(constraint)
end

"""
set_impl_functions!(com::CS.CoM)

Set std.impl.[] for each constraint
"""
function set_impl_functions!(com::CS.CoM; constraints = com.constraints)
    for constraint in constraints
        set_impl_functions!(com, constraint)
    end
end

"""
    set_impl_init!(constraint::Constraint)
Sets `std.impl.init` if the constraint type has a `init_constraint!` method
"""
function set_impl_init!(constraint::Constraint)
    c_type = typeof(constraint)
    c_fct_type = typeof(constraint.fct)
    c_set_type = typeof(constraint.set)
    if hasmethod(init_constraint!, (CS.CoM, c_type, c_fct_type, c_set_type))
        constraint.impl.init = true
    end
end

"""
    set_impl_activate!(constraint::Constraint)
Sets `std.activate.init` if the constraint type has a `activate_constraint!` method
"""
function set_impl_activate!(constraint::Constraint)
    c_type = typeof(constraint)
    c_fct_type = typeof(constraint.fct)
    c_set_type = typeof(constraint.set)
    if hasmethod(activate_constraint!, (CS.CoM, c_type, c_fct_type, c_set_type))
        constraint.impl.activate = true
    end
end

"""
set_impl_update_init!(constraint::Constraint)
Sets `std.impl.update_init` if the constraint type has a `update_init_constraint!` method
"""
function set_impl_update_init!(constraint::Constraint)
    c_type = typeof(constraint)
    c_fct_type = typeof(constraint.fct)
    c_set_type = typeof(constraint.set)
    if hasmethod(
        update_init_constraint!,
        (CS.CoM, c_type, c_fct_type, c_set_type, Vector{<:Constraint}),
    )
        constraint.impl.update_init = true
    end
end

"""
set_impl_update_best_bound!(constraint::Constraint)

Sets `update_best_bound` if the constraint type has a `update_best_bound_constraint!` method
"""
function set_impl_update_best_bound!(constraint::Constraint)
    c_type = typeof(constraint)
    c_fct_type = typeof(constraint.fct)
    c_set_type = typeof(constraint.set)
    if hasmethod(
        update_best_bound_constraint!,
        (CS.CoM, c_type, c_fct_type, c_set_type, Int, Int, Int),
    )
        constraint.impl.update_best_bound = true
    else # just to be sure => set it to false otherwise
        constraint.impl.update_best_bound = false
    end
end

"""
set_impl_reverse_pruning!(constraint::Constraint)
Sets `std.impl.single_reverse_pruning` and `std.impl.reverse_pruning`
if `single_reverse_pruning_constraint!`, `reverse_pruning_constraint!` are implemented for `constraint`.
"""
function set_impl_reverse_pruning!(constraint::Constraint)
    c_type = typeof(constraint)
    c_fct_type = typeof(constraint.fct)
    c_set_type = typeof(constraint.set)
    if hasmethod(
        single_reverse_pruning_constraint!,
        (CS.CoM, c_type, c_fct_type, c_set_type, CS.Variable, Int),
    )
        constraint.impl.single_reverse_pruning = true
    else # just to be sure => set it to false otherwise
        constraint.impl.single_reverse_pruning = false
    end

    if hasmethod(reverse_pruning_constraint!, (CS.CoM, c_type, c_fct_type, c_set_type, Int))
        constraint.impl.reverse_pruning = true
    else # just to be sure => set it to false otherwise
        constraint.impl.reverse_pruning = false
    end
end

"""
set_impl_finished_pruning!(constraint::Constraint)
Sets `std.impl.finished_pruning` if `finished_pruning_constraint!`  is implemented for `constraint`.
"""
function set_impl_finished_pruning!(constraint::Constraint)
    c_type = typeof(constraint)
    c_fct_type = typeof(constraint.fct)
    c_set_type = typeof(constraint.set)
    if hasmethod(finished_pruning_constraint!, (CS.CoM, c_type, c_fct_type, c_set_type))
        constraint.impl.finished_pruning = true
    else # just to be sure => set it to false otherwise
        constraint.impl.finished_pruning = false
    end
end

"""
set_impl_restore_pruning!(constraint::Constraint)
Sets `std.impl.restore_pruning` if `restore_pruning_constraint!`  is implemented for the `constraint`.
"""
function set_impl_restore_pruning!(constraint::Constraint)
    c_type = typeof(constraint)
    c_fct_type = typeof(constraint.fct)
    c_set_type = typeof(constraint.set)
    if hasmethod(
        restore_pruning_constraint!,
        (CS.CoM, c_type, c_fct_type, c_set_type, Union{Int,Vector{Int}}),
    )
        constraint.impl.restore_pruning = true
    else # just to be sure => set it to false otherwise
        constraint.impl.restore_pruning = false
    end
end

"""
call_finished_pruning!(com)

Call `finished_pruning_constraint!` for every constraint which implements that function as saved in `constraint.impl.finished_pruning`
"""
function call_finished_pruning!(com)
    for constraint in com.constraints
        if constraint.impl.finished_pruning
            finished_pruning_constraint!(com, constraint, constraint.fct, constraint.set)
        end
    end
end

"""
call_restore_pruning!(com, prune_steps)

Call `call_restore_pruning!` for every constraint which implements that function as saved in `constraint.impl.restore_pruning`
"""
function call_restore_pruning!(com, prune_steps)
    for constraint in com.constraints
        if constraint.impl.restore_pruning
            restore_pruning_constraint!(
                com,
                constraint,
                constraint.fct,
                constraint.set,
                prune_steps,
            )
        end
    end
end

"""
    count_unfixed(com, constraint::Constraint)

Count number of unfixed variables in the given Constraint
"""
function count_unfixed(com::CS.CoM, constraint::Constraint)
    count(!isfixed, com.search_space[i] for i in constraint.indices)
end

"""
    get_two_unfixed(com, constraint::Constraint)

Get two unfixed indices and local indices of a constraint. One should use
[`count_unfixed`](@ref) to make sure that there are exactly two unfixed indices.
"""
function get_two_unfixed(com::CS.CoM, constraint::Constraint)
    local_vidx_1 = 0
    vidx_1 = 0
    local_vidx_2 = 0
    vidx_2 = 0
    for local_i in 1:length(constraint.indices)
        var = com.search_space[constraint.indices[local_i]]
        if !isfixed(var)
            if local_vidx_1 == 0
                local_vidx_1 = local_i
                vidx_1 = var.idx
            else
                local_vidx_2 = local_i
                vidx_2 = var.idx
                break
            end
        end
    end
    return local_vidx_1, vidx_1, local_vidx_2, vidx_2
end


text:
function plot_lines!(vis::Visualizer, model::Quadruped, q::AbstractVector;
		r=0.0205, offset=0.05, size=10, name::Symbol=:Quadruped, col::Bool=true, α::Real=1.0)
	p_shift = [0.0, 0.0, r]
	orange_mat, blue_mat, black_mat = get_line_material(size, α=α)
	black_mat.color = RGBA(0.0,153/256,153/256,α)
	black_mat.color = RGBA(209/256,0/256,209/256,α)
	black_mat.color = RGBA(136/256,008/256,008/256,α)
	black_mat.color = RGBA(170/256,074/256,068/256,α)

	# Point Traj
	torso_point = Vector{Point{3,Float64}}()
	f1_point = Vector{Point{3,Float64}}()
	f2_point = Vector{Point{3,Float64}}()
	f3_point = Vector{Point{3,Float64}}()
	f4_point = Vector{Point{3,Float64}}()

	for qi in q
		k_torso = kinematics_1(model, qi, body = :torso, mode = :com)
		p_torso = [k_torso[1], -offset, k_torso[2]] + p_shift

		k_calf_1 = kinematics_2(model, qi, body = :calf_1, mode = :ee)
		p_calf_1 = [k_calf_1[1], -offset, k_calf_1[2]] + p_shift

		k_calf_2 = kinematics_2(model, qi, body = :calf_2, mode = :ee)
		p_calf_2 = [k_calf_2[1], -offset, k_calf_2[2]] + p_shift

		k_calf_3 = kinematics_3(model, qi, body = :calf_3, mode = :ee)
		p_calf_3 = [k_calf_3[1], -offset, k_calf_3[2]] + p_shift

		k_calf_4 = kinematics_3(model, qi, body = :calf_4, mode = :ee)
		p_calf_4 = [k_calf_4[1], -offset, k_calf_4[2]] + p_shift

		push!(torso_point, Point(p_torso...))
		push!(f1_point, Point(p_calf_1...))
		push!(f2_point, Point(p_calf_2...))
		push!(f3_point, Point(p_calf_3...))
		push!(f4_point, Point(p_calf_4...))
	end

	# Set lines
	if col
		setobject!(vis[name]["lines/torso"], MeshCat.Line(torso_point, orange_mat))
		setobject!(vis[name]["lines/foot1"], MeshCat.Line(f1_point, blue_mat))
		setobject!(vis[name]["lines/foot2"], MeshCat.Line(f2_point, blue_mat))
		setobject!(vis[name]["lines/foot3"], MeshCat.Line(f3_point, blue_mat))
		setobject!(vis[name]["lines/foot4"], MeshCat.Line(f4_point, blue_mat))
	else
		setobject!(vis[name]["lines/torso"], MeshCat.Line(torso_point, black_mat))
		setobject!(vis[name]["lines/foot1"], MeshCat.Line(f1_point, black_mat))
		setobject!(vis[name]["lines/foot2"], MeshCat.Line(f2_point, black_mat))
		setobject!(vis[name]["lines/foot3"], MeshCat.Line(f3_point, black_mat))
		setobject!(vis[name]["lines/foot4"], MeshCat.Line(f4_point, black_mat))
	end
	return nothing
end

function build_payload!(vis::Visualizer, model::Quadruped; name::Symbol=:Quadruped,
		object=:mesh, r=0.0205, rp=0.20, α=1.0)
	rp = convert(Float32, rp)
	pl_mat = MeshPhongMaterial(color = RGBA(0.4, 0.4, 0.4, α))
	thickness = 0.10
	quad_thickness = 0.267
	default_background!(vis)

	if object == :mesh
		obj = MeshFileObject(joinpath(@__DIR__, "payload_mesh", "Box.obj"))
		setobject!(vis[name][:payload][:object], obj)
		tr = Translation(0.058, quad_thickness/2, 0.18)
		rt = LinearMap(rp*0.5*I(3)*RotZ(π/2)*RotX(pi))
		tform = compose(tr, rt)
		settransform!(vis[name][:payload][:object], tform)
	elseif object == :rectangle
		obj = Rect(Vec(0.058, quad_thickness/2-thickness/2, -rp/2+0.18), Vec(rp, thickness, rp))
		setobject!(vis[name][:payload][:object], obj, pl_mat)
	end
	return nothing
end

function set_payload!(vis::Visualizer, model::Quadruped, q::AbstractVector;
		name::Symbol=:Quadruped, r=0.0205, rp=0.15, offset=0.00)

	r = convert(Float32, r)
	p_shift = [0.0; offset; r]

	p = [q[1]; 0.0; q[2]] + p_shift

	k_torso = kinematics_1(model, q, body = :torso, mode = :ee)
	p_torso = [k_torso[1], 0.0, k_torso[2]] + p_shift
	settransform!(vis[name][:payload], cable_transform(p, p_torso))
	return nothing
end

function animate_payload!(vis::Visualizer, anim::MeshCat.Animation, model::ContactModel,
		q::AbstractVector; name::Symbol=model_name(model))
	for t in 1:length(q)
		MeshCat.atframe(anim, t) do
			set_payload!(vis, model, q[t], name=name)
		end
	end
	setanimation!(vis, anim)
	return nothing
end

function visualize_payload!(vis::Visualizer, model::ContactModel, q::AbstractVector;
		h=0.01, α=1.0, object::Symbol=:mesh,
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=model_name(model))

	build_payload!(vis, model, name=name, object=object, α=α)
	animate_payload!(vis, anim, model, q, name=name)
	return anim
end

function visualize_payload!(vis::Visualizer, model::ContactModel, traj::ContactTraj;
		sample=max(1, Int(floor(traj.H / 100))), h=traj.h*sample,  α=1.0, object::Symbol=:mesh,
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=model_name(model))

	visualize_payload!(vis, model, traj.q[3:sample:end]; anim=anim, name=name, object=object, h=h, α=α)
	return anim
end

function build_robot!(vis::Visualizer, model::Quadruped; name::Symbol=:Quadruped, r=0.0205, α=1.0)
	r = convert(Float32, r)
	body_mat = MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, α))
	contact_mat = MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, α))

	default_background!(vis)

	torso = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_torso),
		convert(Float32, 0.035))
	setobject!(vis[name][:robot]["torso"], torso, body_mat)

	thigh_1 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh1),
		convert(Float32, 0.0175))
	setobject!(vis[name][:robot]["thigh1"], thigh_1, body_mat)

	calf_1 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf1),
		convert(Float32, 0.0125))
	setobject!(vis[name][:robot]["leg1"], calf_1, body_mat)

	thigh_2 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh2),
		convert(Float32, 0.0175))
	setobject!(vis[name][:robot]["thigh2"], thigh_2, body_mat)

	calf_2 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf2),
		convert(Float32, 0.0125))
	setobject!(vis[name][:robot]["leg2"], calf_2, body_mat)

	thigh_3 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh3),
		convert(Float32, 0.0175))
	setobject!(vis[name][:robot]["thigh3"], thigh_3, body_mat)

	calf_3 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf3),
		convert(Float32, 0.0125))
	setobject!(vis[name][:robot]["leg3"], calf_3, body_mat)

	thigh_4 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh4),
		convert(Float32, 0.0175))
	setobject!(vis[name][:robot]["thigh4"], thigh_4, body_mat)

	calf_4 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf4),
		convert(Float32, 0.0125))
	setobject!(vis[name][:robot]["leg4"], calf_4, body_mat)

	hip1 = setobject!(vis[name][:robot]["hip1"],
		GeometryBasics.Sphere(Point3f0(0), convert(Float32, 0.035)), body_mat)
	hip2 = setobject!(vis[name][:robot]["hip2"],
		GeometryBasics.Sphere(Point3f0(0), convert(Float32, 0.035)), body_mat)
	knee1 = setobject!(vis[name][:robot]["knee1"],
		GeometryBasics.Sphere(Point3f0(0), r), body_mat)
	knee2 = setobject!(vis[name][:robot]["knee2"],
		GeometryBasics.Sphere(Point3f0(0), r), body_mat)
	knee3 = setobject!(vis[name][:robot]["knee3"],
		GeometryBasics.Sphere(Point3f0(0), r), body_mat)
	knee4 = setobject!(vis[name][:robot]["knee4"],
		GeometryBasics.Sphere(Point3f0(0), r), body_mat)

	feet1 = setobject!(vis[name][:robot]["feet1"],
		GeometryBasics.Sphere(Point3f0(0), r), contact_mat)
	feet2 = setobject!(vis[name][:robot]["feet2"],
		GeometryBasics.Sphere(Point3f0(0), r), contact_mat)
	feet3 = setobject!(vis[name][:robot]["feet3"],
		GeometryBasics.Sphere(Point3f0(0), r), contact_mat)
	feet4 = setobject!(vis[name][:robot]["feet4"],
		GeometryBasics.Sphere(Point3f0(0), r), contact_mat)
	return nothing
end

function set_robot!(vis::Visualizer, model::Quadruped, q::AbstractVector;
		name::Symbol=:Quadruped, r=0.0205, offset=0.00)

	r = convert(Float32, r)
	p_shift = [0.0; offset; r]

	p = [q[1]; 0.0; q[2]] + p_shift

	k_torso = kinematics_1(model, q, body = :torso, mode = :ee)
	p_torso = [k_torso[1], 0.0, k_torso[2]] + p_shift

	k_thigh_1 = kinematics_1(model, q, body = :thigh_1, mode = :ee)
	p_thigh_1 = [k_thigh_1[1], 0.0, k_thigh_1[2]] + p_shift

	k_calf_1 = kinematics_2(model, q, body = :calf_1, mode = :ee)
	p_calf_1 = [k_calf_1[1], 0.0, k_calf_1[2]] + p_shift

	k_thigh_2 = kinematics_1(model, q, body = :thigh_2, mode = :ee)
	p_thigh_2 = [k_thigh_2[1], 0.0, k_thigh_2[2]] + p_shift

	k_calf_2 = kinematics_2(model, q, body = :calf_2, mode = :ee)
	p_calf_2 = [k_calf_2[1], 0.0, k_calf_2[2]] + p_shift

	k_thigh_3 = kinematics_2(model, q, body = :thigh_3, mode = :ee)
	p_thigh_3 = [k_thigh_3[1], 0.0, k_thigh_3[2]] + p_shift

	k_calf_3 = kinematics_3(model, q, body = :calf_3, mode = :ee)
	p_calf_3 = [k_calf_3[1], 0.0, k_calf_3[2]] + p_shift

	k_thigh_4 = kinematics_2(model, q, body = :thigh_4, mode = :ee)
	p_thigh_4 = [k_thigh_4[1], 0.0, k_thigh_4[2]] + p_shift

	k_calf_4 = kinematics_3(model, q, body = :calf_4, mode = :ee)
	p_calf_4 = [k_calf_4[1], 0.0, k_calf_4[2]] + p_shift

	settransform!(vis[name][:robot]["thigh1"], cable_transform(p, p_thigh_1))
	settransform!(vis[name][:robot]["leg1"], cable_transform(p_thigh_1, p_calf_1))
	settransform!(vis[name][:robot]["thigh2"], cable_transform(p, p_thigh_2))
	settransform!(vis[name][:robot]["leg2"], cable_transform(p_thigh_2, p_calf_2))
	settransform!(vis[name][:robot]["thigh3"], cable_transform(p_torso, p_thigh_3))
	settransform!(vis[name][:robot]["leg3"], cable_transform(p_thigh_3, p_calf_3))
	settransform!(vis[name][:robot]["thigh4"], cable_transform(p_torso, p_thigh_4))
	settransform!(vis[name][:robot]["leg4"], cable_transform(p_thigh_4, p_calf_4))
	settransform!(vis[name][:robot]["torso"], cable_transform(p, p_torso))
	settransform!(vis[name][:robot]["hip1"], MeshCat.Translation(p))
	settransform!(vis[name][:robot]["hip2"], MeshCat.Translation(p_torso))
	settransform!(vis[name][:robot]["knee1"], MeshCat.Translation(p_thigh_1))
	settransform!(vis[name][:robot]["knee2"], MeshCat.Translation(p_thigh_2))
	settransform!(vis[name][:robot]["knee3"], MeshCat.Translation(p_thigh_3))
	settransform!(vis[name][:robot]["knee4"], MeshCat.Translation(p_thigh_4))
	settransform!(vis[name][:robot]["feet1"], MeshCat.Translation(p_calf_1))
	settransform!(vis[name][:robot]["feet2"], MeshCat.Translation(p_calf_2))
	settransform!(vis[name][:robot]["feet3"], MeshCat.Translation(p_calf_3))
	settransform!(vis[name][:robot]["feet4"], MeshCat.Translation(p_calf_4))

	return nothing
end

function contact_point(model::Quadruped, q::AbstractVector)
	k_calf_1 = kinematics_2(model, q, body = :calf_1, mode = :ee)
	p_calf_1 = [k_calf_1[1], 0.0, k_calf_1[2]]

	k_calf_2 = kinematics_2(model, q, body = :calf_2, mode = :ee)
	p_calf_2 = [k_calf_2[1], 0.0, k_calf_2[2]]

	k_calf_3 = kinematics_3(model, q, body = :calf_3, mode = :ee)
	p_calf_3 = [k_calf_3[1], 0.0, k_calf_3[2]]

	k_calf_4 = kinematics_3(model, q, body = :calf_4, mode = :ee)
	p_calf_4 = [k_calf_4[1], 0.0, k_calf_4[2]]

	pc = [p_calf_1, p_calf_2, p_calf_3, p_calf_4]
	return pc
end

function build_meshrobot!(vis::Visualizer, model::Quadruped; name::Symbol=:Quadruped, α=1.0)
	default_background!(vis)

	urdf = joinpath(@__DIR__, "mesh", "a1.urdf")
	package_path = @__DIR__
	build_meshrobot!(vis, model, urdf, package_path; name=name, α=α)
end

function convert_config(model::Quadruped, q::AbstractVector)
    # Quadruped configuration
    # 1.  position long axis +x
    # 2.  position long axis +z
    # 3.  trunk rotation along -y
    # 4.  back  left  shoulder rotation along -y
    # 5.  back  left  elbow    rotation along -y
    # 6.  back  right shoulder rotation along -y
    # 7.  back  right elbow    rotation along -y
    # 8.  front right shoulder rotation along -y
    # 9.  front right elbow    rotation along -y
	# 10. front left  shoulder rotation along -y
	# 11. front left  elbow    rotation along -y

    # URDF configuration
    # 1.  front right clavicle rotation along +x
    # 2.  front left  clavicle rotation along +x
    # 3.  back  right clavicle rotation along +x
    # 4.  back  left  clavicle rotation along +x
    # 5.  front right shoulder rotation along +y
    # 6.  front left  shoulder rotation along +y
    # 7.  back  right shoulder rotation along +y
    # 8.  back  left  shoulder rotation along +y
    # 9.  front right elbow    rotation relative to shoulder along +y
    # 10. front left  elbow    rotation relative to shoulder along +y
    # 11. back  right elbow    rotation relative to shoulder along +y
    # 12. back  left  elbow    rotation relative to shoulder along +y
    q_ = zeros(12)
    x, z, θ = q[1:3]
    q_[5]  = -q[10] + θ - π/2
    q_[6]  = -q[8] + θ - π/2
	q_[7]  = -q[4] + θ - π/2
    q_[8]  = -q[6] + θ - π/2
    q_[9]  = +q[10]-q[11]
    q_[10] = +q[8]-q[9]
	q_[11] = +q[4]-q[5]
    q_[12] = +q[6]-q[7]

	r_foot = 0.02
	trunk_length = 2*0.183
	trunk_width = 2*0.132
    tform = compose(
		Translation(x , trunk_width/2, z+r_foot),
		compose(
			LinearMap(AngleAxis(π/2-θ, 0, 1.0, 0)),
			Translation(trunk_length/2, 0,0)
			)
		)
    return q_, tform
end


text:
#example.jl
function acronym(phrase)
    regex = "[A-Z]+['a-z]*|['a-z]+"
    return join([uppercase(word.match[1]) for word in eachmatch(Regex(regex), phrase)], "")
end


text:
chunkedtest1 = include("./data/chunked/test1.jl")

@testset "read" begin
    io = IOBuffer(chunkedtest1)
    s = ChunkedTransferDecodeStream(io)

    data = read(s)

    expected = "'Twas brillig, and the slithy toves\n" *
        "      Did gyre and gimble in the wabe:\n" *
        "All mimsy were the borogoves,\n" *
        "      And the mome raths outgrabe."

    @test String(data) == expected

    close(s)
end


text:
module SuiteSparseGraphBLAS
__precompile__(true)

using Libdl: dlsym, dlopen, dlclose

# Allow users to specify a non-Artifact shared lib.
using Preferences
include("find_binary.jl")
const libgraphblas_handle = Ref{Ptr{Nothing}}()
@static if artifact_or_path == "default"
    using SSGraphBLAS_jll
    const libgraphblas = SSGraphBLAS_jll.libgraphblas
else
    const libgraphblas = artifact_or_path
end

using SparseArrays
using SparseArrays: nonzeroinds
using MacroTools
using LinearAlgebra
using Random: randsubseq, default_rng, AbstractRNG, GLOBAL_RNG
using SpecialFunctions: lgamma, gamma, erf, erfc
using Base.Broadcast
include("abstracts.jl")
include("libutils.jl")
include("lib/LibGraphBLAS.jl")
using .libgb
include("operators/libgbops.jl")

# Globals

include("types.jl")
include("gbtypes.jl")

include("constants.jl")

include("operators/operatorutils.jl")
include("operators/unaryops.jl")
include("operators/binaryops.jl")
include("operators/monoids.jl")
include("operators/semirings.jl")
include("operators/selectops.jl")
include("descriptors.jl")
using .UnaryOps
using .BinaryOps
using .Monoids
using .Semirings

# Create typed operators
_createunaryops()
_createbinaryops()
_createmonoids()
_createsemirings()

include("operators/oplist.jl")
include("indexutils.jl")

include("operations/extract.jl")
include("scalar.jl")
include("vector.jl")
include("matrix.jl")
include("random.jl")

include("operations/operationutils.jl")
include("operations/transpose.jl")
include("operations/mul.jl")
include("operations/ewise.jl")
include("operations/map.jl")
include("operations/select.jl")
include("operations/reduce.jl")
include("operations/kronecker.jl")
include("operations/concat.jl")
include("operations/resize.jl")

include("print.jl")
include("import.jl")
include("export.jl")
include("pack.jl")
include("unpack.jl")
include("options.jl")
#EXPERIMENTAL
include("operations/argminmax.jl")
include("operations/broadcasts.jl")
include("chainrules/chainruleutils.jl")
include("chainrules/mulrules.jl")
include("chainrules/ewiserules.jl")
include("chainrules/maprules.jl")
include("chainrules/reducerules.jl")
include("chainrules/selectrules.jl")
include("chainrules/constructorrules.jl")
#include("random.jl")
include("misc.jl")
include("asjulia.jl")
include("sparsemat.jl")
export SparseArrayCompat
export libgb
export UnaryOps, BinaryOps, Monoids, Semirings #Submodules
export UnaryOp, BinaryOp, Monoid, Semiring #UDFs
export Descriptor #Types
export gbset, gbget
export xtype, ytype, ztype, validtypes #Determine input/output types of operators
export GBScalar, GBVector, GBMatrix #arrays
export lgamma, gamma, erf, erfc #reexport of SpecialFunctions.

# Function arguments not found elsewhere in Julia
#UnaryOps not found in Julia/stdlibs.
export frexpe, frexpx, positioni, positionj
#BinaryOps not found in Julia/stdlibs.
export second, rminus, pair, ∨, ∧, lxor, fmod, firsti,
    firstj, secondi, secondj
#SelectOps not found in Julia/stdlibs
export offdiag

export clear!, extract, extract!, subassign!, assign!, hvcat! #array functions

#operations
export mul, select, select!, eadd, eadd!, emul, emul!, map, map!, gbtranspose, gbtranspose!,
gbrand, eunion, eunion!
# Reexports from LinAlg
export diag, diagm, mul!, kron, kron!, transpose, reduce, tril, triu

# Reexports from SparseArrays
export nnz, sprand, findnz, nonzeros, nonzeroinds

function __init__()
    @static if artifact_or_path != "default"
        libgraphblas_handle[] = dlopen(libgraphblas)
    else
        #The artifact does dlopen for us.
        libgraphblas_handle[] = SSGraphBLAS_jll.libgraphblas_handle
    end
    _load_globaltypes()
    # We initialize GraphBLAS by giving it Julia's GC wrapped memory management functions.
    # In the future this should hopefully allow us to do no-copy passing of arrays between Julia and SS:GrB.
    # In the meantime it helps Julia respond to memory pressure from SS:GrB and finalize things in a timely fashion.
    libgb.GxB_init(libgb.GrB_NONBLOCKING, cglobal(:jl_malloc), cglobal(:jl_calloc), cglobal(:jl_realloc), cglobal(:jl_free))
    gbset(:nthreads, Sys.CPU_THREADS ÷ 2)
    # Eagerly load selectops constants.
    _loadselectops()
    # Set printing done by SuiteSparse:GraphBLAS to base-1 rather than base-0.
    gbset(BASE1, 1)
    atexit() do
        # Finalize the lib, for now only frees a small internal memory pool.
        libgb.GrB_finalize()
        @static if artifact_or_path != "default"
            dlclose(libgraphblas_handle[])
        end
    end
end

include("operators/ztypes.jl")
end #end of module


text:
using DataFrames: DataFrame
using Latexify
using Test
d = DataFrame(A = 11:13, B = [:X, :Y, :Z])

@test latexify(d; env=:table, side=1:3, latex=false) ==
raw"\begin{tabular}{ccc}
 & A & B\\
1 & 11 & X\\
2 & 12 & Y\\
3 & 13 & Z\\
\end{tabular}
"

arr = ["x/(y-1)", 1.0, 3//2, :(x-y), :symb]

M = vcat(hcat(arr...), hcat(arr...))
head = ["col$i" for i in 1:size(M, 2)]
side = ["row$i" for i in 1:size(M, 1)]

@test latexify(M; env=:table, head=1:2, adjustment=:l, latex=false, transpose=true) ==
raw"\begin{tabular}{ll}
1 & 2\\
x/(y-1) & x/(y-1)\\
1.0 & 1.0\\
3//2 & 3//2\\
x - y & x - y\\
symb & symb\\
\end{tabular}
"

@test latexify(M; env=:table, head=1:2, adjustment=:l, transpose=true) ==
raw"\begin{tabular}{ll}
$1$ & $2$\\
$\frac{x}{y - 1}$ & $\frac{x}{y - 1}$\\
$1.0$ & $1.0$\\
$\frac{3}{2}$ & $\frac{3}{2}$\\
$x - y$ & $x - y$\\
$symb$ & $symb$\\
\end{tabular}
"

@test latexify(M; env=:table, head=1:5) ==
raw"\begin{tabular}{ccccc}
$1$ & $2$ & $3$ & $4$ & $5$\\
$\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
$\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
\end{tabular}
"

@test latexify(M; env=:table, side=[1, 2]) ==
raw"\begin{tabular}{cccccc}
$1$ & $\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
$2$ & $\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
\end{tabular}
"

@test latexify(M; env=:table, booktabs=true) ==
raw"\begin{tabular}{ccccc}
\toprule
$\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
$\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
\bottomrule
\end{tabular}
"

@test latexify(M; env=:table, head=1:5, booktabs=true) ==
raw"\begin{tabular}{ccccc}
\toprule
$1$ & $2$ & $3$ & $4$ & $5$\\
\midrule
$\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
$\frac{x}{y - 1}$ & $1.0$ & $\frac{3}{2}$ & $x - y$ & $symb$\\
\bottomrule
\end{tabular}
"

D = Dict(:a=>"x/(k+x)", :b=>"x - y")
@test latexify(D; env=:tabular) ==
raw"\begin{tabular}{cc}
$a$ & $\frac{x}{k + x}$\\
$b$ & $x - y$\\
\end{tabular}
"

@test latexify(D; env=:tabular, head=["Keys", "Values"]) ==
raw"\begin{tabular}{cc}
$Keys$ & $Values$\\
$a$ & $\frac{x}{k + x}$\\
$b$ & $x - y$\\
\end{tabular}
"


text:
function logprior(real::State{A1, A2, A3},
                  tran::State{A1, A2, A3},
                  mps::StateMP,
                  c::Constants) where {A1, A2, A3}

  lp = 0.0
  for key in fieldnames(State)
    if !(key in (:y_m, :y_log_s))
      if key == :v
        a = tran.alpha[1]
        lpdf = ADVI.compute_lpdf(c.priors.v(a), tran.v)
        labsdj = ADVI.logabsdetJ(mps.v, real.v, tran.v)
      else # key != :v
        p = getfield(c.priors, key)
        t = getfield(tran, key)
        r = getfield(real, key)
        mp = getfield(mps, key)

        lpdf = ADVI.compute_lpdf(p, t)
        labsdj = ADVI.logabsdetJ(mp, r, t)
      end
      lp += sum(lpdf + labsdj)
    end
  end

  return lp
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using GeometryExperiments
using Test

const P2 = Point{2,Float64}
const P3 = Point{3,Float64}

@testset "GeometryExperiments.jl" begin
    @testset "Transforms" begin
        @testset "Translations" begin
            tr = Translation(1., 3.)
            tr_inv = inv(tr)
            @test tr ∘ tr_inv == Translation(0., 0.)
            @test tr(Point(1., 2.)) == Point(2., 5.)
        end

        @testset "Scalings" begin
            sc = Scaling(1., 2.)
            sc_inv = inv(sc)
            @test sc ∘ sc_inv == Scaling(1., 1.)
            @test sc(Point(1., 2.)) == Point(1., 4.)

            us = UniformScaling(2.)
            @test us ∘ inv(us) == UniformScaling(1.)
            @test us(Point(1., 2.)) == Point(2., 4.)
        end

        @testset "Rotations" begin
        end

        @testset "Composition" begin
            @test Translation(2., 3.) ∘ Translation(1., 2.) == Translation(3., 5.)
            @test Scaling(2., 3.) ∘ Scaling(1., 2.) == Scaling(2., 6.)

            tr = Translation(1., 2.) ∘ Scaling(1., 2.) ∘ Translation(-1., -2.)
            @test all(transforms(tr) .== [Translation(1., 2.), Scaling(1., 2.), Translation(-1., -2.)])

            tr = Scaling(2., 3.) ∘ Translation(1., 2.)
            @test tr(Point(0., 0.)) == Point(2., 6.)

            p = Point(1., 2.)
            @test (inv(tr) ∘ tr)(p) == p
        end
    end

    @testset "Geometry" begin
        eval_sph(radius, p::Point) = hypot(p...) - radius

        p = Point(0., 1., 0.)

        hc = HyperCube(0.2)
        @test p ∉ hc
        @test Point(0., 0., 0.) ∈ hc
        @test Translated(hc, Translation(0., 0., 0.))(p) == hc(p) == 0.8
        @test Translated(hc, Translation(0.05, 1., 0.))(p) ≈ -0.15
        @test Transformed(hc, Scaling(1., 2., 3.)) === Box(0.2, Scaling(1., 2., 3.))

        sph = HyperSphere(0.2)
        @test p ∉ sph
        @test Point(0., 0., 0.) ∈ sph
        @test sph(Point(0.3, 0.4, 0.5)) == eval_sph(0.2, Point(0.3, 0.4, 0.5))
        @test sph(Point(0.5, 0.2)) == eval_sph(0.2, Point(0.5, 0.2))
        @test Translated(sph, Translation(0., 0., 0.))(p) == sph(p) == 0.8
        @test Translated(sph, Translation(0.05, 1., 0.))(p) ≈ -0.15
        @test p ∈ Translated(sph, Translation(0.05, 1., 0.))

        elps = Ellipsoid(0.2, Scaling(1., 2., 3.))
        @test Scaled(sph, Scaling(1., 2., 3.)) === elps
        @test elps ≈ Ellipsoid(Point(0.2, 0.4, 0.6))
        @test p ∉ elps

        @testset "Projections" begin
            p = Projection{2}(HyperSphere(eval_sph(0., Point(0.2, 0.2))))
            @test Point(0.1,0.1) ∈ p
            @test Point(0.1,2.) ∉ p
            @test Point(0.1,2.,0.) ∉ p
            @test Point(0.1,0.1,0.) ∈ p
            @test Point(0.1,0.1,0.1) ∉ p
            @test p(Point(0.2,0.2)) == 0.
            @test p(Point(0.2,0.2,0.)) == 0.
            @test p(Point(0.2,0.2,0.5)) == 0.5
            @test p(Point(0.2,0.2,0.5,0.3)) == eval_sph(0., Point(0.5,0.3))
        end
    end

    @testset "Point sets" begin
        set = PointSet([Point(0., 0.), Point(1., 0.), Point(0., 1.), Point(1., 1.)])
        set2 = PointSet([Point(0., 0.5), Point(0.5, 0.), Point(0.5, 1.), Point(1., 0.5)])
        @test boundingelement(set) == Translated(Scaled(HyperCube(1.), Scaling(0.5, 0.5)), Translation(0.5, 0.5))
        @test boundingelement(set) == boundingelement(set2)

        set = PointSet([Point(-1., -1.), Point(1., -1.), Point(-1., 1.), Point(1., 1.)])
        @test boundingelement(set) == Translated(Scaled(HyperCube(1.), Scaling(1., 1.)), Translation(0., 0.))
        @test Scaling(1., 1.)(set) == set
        @test Translation(0., 0.)(set) == set
        (Scaling(1., 1.) ∘ Translation(0., 0.))(set) == set

        set = PointSet([Point(0., 0., 0.), Point(1., 0., 0.), Point(0., 1., 0.), Point(0., 0., 1.), Point(1., 1., 0.), Point(1., 0., 1.), Point(0., 1., 1.), Point(1., 1., 1.)])
        @test boundingelement(set) == Translated(Scaled(HyperCube(1.), Scaling(0.5, 0.5, 0.5)), Translation(0.5, 0.5, 0.5))

        @test PointSet(HyperCube, P2) == PointSet(P2[(-1,-1),(1,-1),(-1,1),(1,1)])
        @test PointSet(HyperCube(0.5), P2) == PointSet(P2[(-0.5,-0.5),(0.5,-0.5),(-0.5,0.5),(0.5,0.5)])
        @test PointSet(HyperCube, P3) == PointSet(P3[(-1,-1,-1),(1,-1,-1),(-1,1,-1),(1,1,-1),(-1,-1,1),(1,-1,1),(-1,1,1),(1,1,1)])

        @testset "Nearest" begin
            set = PointSet([Point(0., 1.), Point(0.5, 0.5), Point(1., 0)])
            @test sort_nearest(set, Point(1., 0.75)) == set.points[[2, 3, 1]]
        end
    end

    @testset "Curves" begin
        @testset "Bezier curves" begin
            T = Float64
            b = BezierCurve()
            points = P2[(0,0),(0.5,1),(1,0)]
            @test b(T(0), points) == P2(0,0)
            @test b(T(1), points) == P2(1,0)
            @test b(T(0.5), points) == P2(0.5,0.5)
            @test b(T(0.5), points) == P2(0.5,0.5)
            @test_throws DomainError(T(-0.1), "b(t) is not defined for t outside [0, 1].") b(T(-0.1), points)
            @test_throws DomainError(T(1.2), "b(t) is not defined for t outside [0, 1].") b(T(1.2), points)
        end

        @testset "Patches" begin
            p = Patch(BezierCurve(), 3)
            points = P2[(0,0),(0.5,1),(1,0), (1.5,-1),(2,0)]
            @test p(0, points) == P2(0,0)
            @test p(0.25, points) == P2(0.5, 0.5)
            @test p(0.5, points) == P2(1,0)
            @test p(0.75, points) == P2(1.5, -0.5)
            @test p(1, points) == P2(2,0)
            @test split(points, p) == [@view(points[1:3]), @view(points[3:5])]
            @test curve_points(p, points) == @view points[1:3]
            @test curve_points(p, points, 1) == @view points[3:5]
        end
    end

    @testset "Mesh encodings" begin
        strip = TriangleStrip(1:5)
        list = TriangleList([(1, 2, 3), (2, 4, 3), (3, 4, 5)])
        @test TriangleList(strip) == list

        fan = TriangleFan(1:5)
        list = TriangleList([(1, 2, 3), (1, 3, 4), (1, 4, 5)])
        @test GeometryExperiments.topology_class(typeof(fan)) == Triangle
        @test TriangleList(fan) == list

        @test MeshVertexEncoding(P2[(1.2,1.4),(0.1,0.2),(0.3,0.4),(0.5,0.2)], Triangle) isa MeshVertexEncoding{<:TriangleStrip,P2}
        @test MeshVertexEncoding(fan, P2[(1.2,1.4),(0.1,0.2),(0.3,0.4),(0.5,0.2)]) isa MeshVertexEncoding{<:TriangleFan,P2}

        strip = LineStrip(1:5)
        @test GeometryExperiments.topology_class(typeof(strip)) == Line
        @test LineList(strip) == LineList([(1, 2), (2, 3), (3, 4), (4, 5)])
    end
end


text:
using Plots

rotation(x,y,theta)= [cos(theta) -sin(theta); sin(theta) cos(theta)]*[x;y]
points=[[0.0;0.0],[1.0;0.0]]
tempo=(points[2]-points[1])/sqrt(2)
point=rotation(tempo[1], tempo[2], pi/4)+points[1]
insert!(points,2, point)

for step in 1:20

    i=1
    while points[i]!= points[end]
        temp1=(points[i+1]-points[i])/sqrt(2)
        point1=rotation(temp1[1], temp1[2], +pi/4)+points[i]
        insert!(points, i+1, point1)
        temp2=(points[i+3]-points[i+2])/sqrt(2)
        point2=rotation(temp2[1], temp2[2], -pi/4)+points[i+2]
        insert!(points,i+3, point2)
        i+=4
    end
end

plot(hcat(points...)[1,1:floor(Int,length(points)/2)] , hcat(points...)[2,1:floor(Int, length(points)/2)], color=:blue, border=:none)
plot!(hcat(points...)[1,floor(Int, length(points)/2):end], hcat(points...)[2,floor(Int, length(points)/2):end], color=:red, border=:none)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter2\\report\\Fig\\DragonFractal.png")


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# We can't guarantee that these are correct, but we can at least check
# that they run
@test length(Base.Sys.cpu_info()) > 0
sprint(Base.Sys.cpu_summary)
@test Base.Sys.uptime() > 0
Base.Sys.loadavg()