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

@doc raw"""IngressSpec describes the Ingress the user wishes to exist.

    IoK8sApiNetworkingV1beta1IngressSpec(;
        backend=nothing,
        rules=nothing,
        tls=nothing,
    )

    - backend::IoK8sApiNetworkingV1beta1IngressBackend : A default backend capable of servicing requests that don't match any rule. At least one of 'backend' or 'rules' must be specified. This field is optional to allow the loadbalancer controller or defaulting logic to specify a global default.
    - rules::Vector{IoK8sApiNetworkingV1beta1IngressRule} : A list of host rules used to configure the Ingress. If unspecified, or no rule matches, all traffic is sent to the default backend.
    - tls::Vector{IoK8sApiNetworkingV1beta1IngressTLS} : TLS configuration. Currently the Ingress only supports a single TLS port, 443. If multiple members of this list specify different hosts, they will be multiplexed on the same port according to the hostname specified through the SNI TLS extension, if the ingress controller fulfilling the ingress supports SNI.
"""
mutable struct IoK8sApiNetworkingV1beta1IngressSpec <: SwaggerModel
    backend::Any # spec type: Union{ Nothing, IoK8sApiNetworkingV1beta1IngressBackend } # spec name: backend
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1beta1IngressRule} } # spec name: rules
    tls::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1beta1IngressTLS} } # spec name: tls

    function IoK8sApiNetworkingV1beta1IngressSpec(;backend=nothing, rules=nothing, tls=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1beta1IngressSpec, Symbol("backend"), backend)
        setfield!(o, Symbol("backend"), backend)
        validate_property(IoK8sApiNetworkingV1beta1IngressSpec, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        validate_property(IoK8sApiNetworkingV1beta1IngressSpec, Symbol("tls"), tls)
        setfield!(o, Symbol("tls"), tls)
        o
    end
end # type IoK8sApiNetworkingV1beta1IngressSpec

const _property_map_IoK8sApiNetworkingV1beta1IngressSpec = Dict{Symbol,Symbol}(Symbol("backend")=>Symbol("backend"), Symbol("rules")=>Symbol("rules"), Symbol("tls")=>Symbol("tls"))
const _property_types_IoK8sApiNetworkingV1beta1IngressSpec = Dict{Symbol,String}(Symbol("backend")=>"IoK8sApiNetworkingV1beta1IngressBackend", Symbol("rules")=>"Vector{IoK8sApiNetworkingV1beta1IngressRule}", Symbol("tls")=>"Vector{IoK8sApiNetworkingV1beta1IngressTLS}")
Base.propertynames(::Type{ IoK8sApiNetworkingV1beta1IngressSpec }) = collect(keys(_property_map_IoK8sApiNetworkingV1beta1IngressSpec))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1beta1IngressSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1beta1IngressSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1beta1IngressSpec }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1beta1IngressSpec[property_name]

function check_required(o::IoK8sApiNetworkingV1beta1IngressSpec)
    true
end

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


text:
"""
`PersistentArray{T,N}` wraps an Array{Node, N} to
provide a versioned array interface. The Node type
maintains a history of each element, so the complexity
of each operation on the PersistentArray is dependent on
the corresponding Node operation in relation to the number
of stored versions.
"""
type PersistentArray{T,N} <: AbstractArray{T,N}
    data::AbstractArray
    version
end

function PersistentArray(array::AbstractArray; version=version)
    PersistentArray{Node,size(array)}(array, version)
end

"""
`PersistentArray{V,T}(dimss::Int...; default_value::T=0.0, default_version::V=1)`
builds a persistent array where each Node contains 1 version
(default_version, default_value).
"""
function PersistentArray{V,T}(dims::Int...; default_value::T=0.0, default_version::V=1)
    a = Array(Node{V,T}, dims...)
    map(
        i -> a[i] = Node(default_version, default_value),
        eachindex(a)
    )
    PersistentArray{Node{V,T},length(dims)}(a, default_version)
end

"""
`Base.show(array::PersistentArray)` prints out the array version and data.

NOTE: this doesn't get used by the REPL for some reason.
"""
function Base.show(array::PersistentArray)
    print("PersistentArray(version=$(array.version), data=")
    show(array.data)
    print(")")
end

"`version(array::PersistentArray)` prints the current array version."
version(array::PersistentArray) = array.version

"`commit!(array::PersistentArray` increments the array version."
commit!(array::PersistentArray) = array.version = increment(array.version)

"`increment(x::Int)` provides the increment method for ints (default version type)"
increment(x::Int) = x + 1

"""
`update!{V,T}(array::PersistentArray{Node{V,T}}, x::T, version::V, idx...)`
updates or adds the version with the specifid value for each element specified
by `idx`. Runs in O(mlog(m)) time for each element selected where m is the
number of versions.
"""
function update!{V,T}(array::PersistentArray{Node{V,T}}, x::T, version::V, idx...)
    map(
        i -> update!(i, x, version),
        sub(array.data, idx...)
    )
end

"""
`lookup{V,T}(array::PersistentArray{Node{V,T}}, version::V, idx...)` looks up the
appropriate element selected at the given `version`. Runs in O(log(m)) for each element
selected where m is the number of versions.
"""
function lookup{V,T}(array::PersistentArray{Node{V,T}}, version::V, idx...)
    map(
        i -> lookup(i, version),
        sub(array.data, idx...)
    )
end

"""
`getindex{V,T}(array::PersistentArray{Node{V,T}}, idx...)` returns the
current version of all elements selected. Runs in O(1) for each element.
"""
function Base.getindex{V,T}(array::PersistentArray{Node{V,T}}, idx...)
    result = map(
        i -> getlatest(i),
        sub(array.data, idx...)
    )
    if length(result) == 1
        return result[1]
    else
        return result
    end
end

"""
`setindex!{V,T}(array::PersistentArray{Node{V,T}}, x::T, idx...)`
sets the current version of every node selected. If the node doesn't have a
current version it is added. Runs in O(m) time for each element where
m is the number of versions.
"""
function Base.setindex!{V,T}(array::PersistentArray{Node{V,T}}, x::T, idx...)
    map(
        i -> setlatest!(i, x, array.version),
        sub(array.data, idx...)
    )
end

# Required for subtype of AbstractArray
Base.size(array::PersistentArray, idx::Int...) = size(array.data, idx...)
Base.linearindexing{T<:PersistentArray}(::Type{T}) = Base.LinearFast()
Base.ndims(array::PersistentArray) = ndims(array.data)

ext: jl


text:
@testset "437.path-sum-iii.jl" begin
    @test tree_path_sum(convert(TreeNode{Int},
                                [10, 5, -3, 3, 2, nothing, 11, 3, -2, nothing, 1]), 8) == 3
end


text:
using Turf

geo_data = GeoJSON.parsefile(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/points-inside-polygon/points-inside-polygon.geojson")

polygon = geo_data.features[1].geometry

for i in 2:length(geo_data.features)
    feature = geo_data.features[i]
    point = feature.geometry

    if within(point, polygon)
        feature.properties = Dict("marker-color" => "#ff0000")
    end
end

result = geojson(geo_data)

open(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/points-inside-polygon/points-inside-polygon.result.geojson", "w") do file
    write(file, result)
end


text:
using NPZ
using JuMP
using Ipopt
using LinearAlgebra
using WoodburyMatrices
using .Threads
using ProgressMeter

X1 = npzread("Kernels/G1.npy")
X2 = npzread("Kernels/G2.npy")
X3 = npzread("Kernels/G4.npy") #try on restricted grid
nfactor = npzread("Testing/nfactor.npy")

ϵ = 1e-2
lnorm = -log(ϵ)*(2*size(X1,2))

function d_oed(X1::Matrix{T}, X2::Matrix{T}, X3::Matrix{T}, nfactor, B, ϵ, λ, lnorm; maxiter=1000) where T <: Number
    X = [X3 zeros(T, size(X3));
         zeros(T, size(X3)) X3;
         nfactor*X1 nfactor*X2]
    Xn = maximum(sqrt.(sum(X.^2, dims=2)))
    X ./= Xn

    n = size(X3, 1)
    ndas = size(X1, 1)

    function Ml(μ, X, ϵ)
        M = ϵ*I + X'*Diagonal(μ)*X
        Mc = cholesky(Hermitian(M))
        l = μ.^2 .*diag(X*(Mc\X'))
        return (M, Mc, l)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function f(μ, λ, lnorm, n)
        M, Mc, l = Ml(μ)
        l2 = l[1:2*n] # coherence of point sensors only
        return -logdet(M) + λ*maximum(l2) - lnorm
    end

    f(μ) = f(μ, λ, lnorm, n)

    grad_f(y) = f(vcat([y[i] for i in 1:length(y)], [y[i] for i in 1:length(y)], ones(T, ndas)))

    jump_wrapper_f(x...) = f(vcat([x[i] for i in 1:length(x)], [x[i] for i in 1:length(x)], ones(T, ndas)))

    function gi(i, Mc, μ, u, lmind, X, λ, T)
        ld_term = @views -X[i,:]'*(Mc\X[i,:])
        tmp_A = SymWoodbury(Mc, X[i,:], -μ[i])
        tmp_b = @views tmp_A \ X[i,:]
        tmp_1 = (u'*tmp_b)^2
        tmp_d = @views X[i,:]'*tmp_b
        tmp_den = (1+μ[i]*tmp_d)^2
        lev_term_1 = -μ[lmind]^2*((1+μ[i]*tmp_d)*tmp_1 - μ[i]*tmp_1*tmp_d) / tmp_den
        lev_term_2 = (i==lmind) ? 2*μ[lmind]*(u'*(Mc\u)) : zero(T)
        return ld_term + λ*(lev_term_1 + lev_term_2)
    end

    gi(i, Mc, μ, u, lmind) = gi(i, Mc, μ, u, lmind, X, λ, T)

    function g!(gvec, μ, X, n)
        M, Mc, l = Ml(μ)
        l2 = l[1:n] .+ l[(n+1):2*n] # coherence of point sensors only; sum across channels
        lm, lmind = findmax(l2)
        u1 = @views X[lmind, :]
        u2 = @views X[lmind, :]
        @threads for i = 1:n
            g1 = gi(i, Mc, μ, u1, lmind) # in respect to x component
            g2 = gi(i+n, Mc, μ, u2, lmind) # in respect to y component
            gvec[i] = g1+g2
        end
    end

    g!(gvec, μ) = g!(gvec, μ, X, n)

    jump_wrapper_g!(gvec, x...) = g!(gvec, vcat([x[i] for i in 1:length(x)], [x[i] for i in 1:length(x)], ones(T, ndas)))

    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter"=>maxiter, "tol"=>1e-3))
    JuMP.register(model, :f, n, jump_wrapper_f, jump_wrapper_g!)

    @variable(model, 0.0<=x[1:n]<= 1.0)

    @NLconstraint(model, +(x...) <= B)

    @NLobjective(model, Min, f(x...))

    optimize!(model)

    function ff(μ, X, ϵ, λ, n)
        M = ϵ*I + X'*Diagonal(μ)*X
        l = (μ.^2).*diag(X*(M\X'))
        l2 = l[1:n] .+ l[(n+1):2*n] # coherence of point sensors only; sum across channels
        return -logdet(M) + λ*maximum(l2) - lnorm
    end

    ff(μ) = ff(μ, X, ϵ, λ, n)

    μr = clamp.(value.(x), 0, 1) # because it probably hasn't completely met the constraints, make sure it is feasible.
    μ = vcat(μr, μr, ones(T, ndas))

    function pround(μ, μr, n, B)
        p = Progress(n, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)
        idx = 1:length(μr)
        rμr = copy(μr)
        v1 = copy(μr)
        v2 = copy(μr)
        μ1 = copy(μ)
        μ2 = copy(μ)
        not_zero_or_one(x) = !(isapprox(x,0) || isapprox(x,1))
        feasible(v) = (all(v.>=0) && all(v.<=1) && (sum(v)<=B))
        while length(idx[not_zero_or_one.(rμr)]) > 1
            # j, k = rand(idx[not_zero_or_one.(rμr)], 2)
            idr = idx[not_zero_or_one.(rμr)]
            j = idr[argmax(rμr[idr])]
            k = idr[argmin(rμr[idr])]
            δ = min(1-rμr[j], rμr[k])
            τ = min(1-rμr[k], rμr[j])
            v1[j] += δ
            v1[k] -= δ
            v2[j] -= τ
            v2[k] += τ
            sv1 = sum(v1)
            sv2 = sum(v2)
            if !feasible(v1) && !feasible(v2)
                error("No Feasible solution")
            elseif !feasible(v1)
                rμr .= v2
                v1 .= v2
            elseif !feasible(v2)
                rμr .= v1
                v2 .= v1
            else
                μ1[1:n] .= v1
                μ1[(n+1):(2*n)] .= v1
                f1 = ff(μ1)
                μ2[1:n] .= v2
                μ2[(n+1):(2*n)] .= v2
                f2 = ff(μ2)
                if f1 < f2
                    rμr .= v1
                else
                    rμr .= v2
                end
            end
            next!(p)
        end
        μ1[1:n] .= rμr
        μ1[(n+1):(2*n)] .= rμr
        fr = ff(μ1)
        return (rμr, fr)
    end

    #rμr, fr = pround(μ, μr, n, B)

    npzwrite("IncoherenceResults/res_inc_$B.npy", value.(x))
    #npzwrite("IncoherenceResults/res_inc_$(B)_round.npy", rμr)

    return (model, x)#, rμr, fr)
end

for i in [4,8,16,32,64]
    # model, x, rμr, fr = d_oed(X1, X2, X3, nfactor, i, ϵ, 1/ϵ, lnorm, maxiter=500)
#     model, x = d_oed(X1, X2, X3, nfactor, i, ϵ, 1/ϵ, lnorm, maxiter=500)
    model, x = d_oed(X1, X2, X3, nfactor, i, ϵ, 1.0, lnorm, maxiter=50)
end

# model, x, rμr, fr = d_oed(X1, X2, X3, nfactor, 32, ϵ, 1/ϵ, lnorm, maxiter=2)
# println("Model Objective = $(objective_value(model))")
# # npzwrite("res_mixed.npy", value.(x))
# # ;

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
export SR, sr_none, sr_multiplicative, sr_shift

@enum SRPreconditionType sr_none=1 sr_shift=2 sr_multiplicative=3
################################################################################
######   Preconditioning algorithm definition (structure holding params)  ######
################################################################################
"""
    SR([use_iterative=true, ϵ=0.001, λ0=100, b=0.95, λmin=1e-4, precondition_type=sr_shift)

Stochastic Reconfiguration preconditioner which corrects the gradient according
to the natural gradient computed as S^-1 ∇C. Using this algorithm will lead to
the computation of the S matrix together with the gradient of the cost function
∇C. To compute the natural gradient S^-1∇C an iterative scheme (Minres-QLP) or
a direct inversion is used.

If `use_iterative=true` the inverse matrix `S^-1` is not computed, and an iterative
MINRES-QLP algorithm is used to compute the product S^-1*F

If `precondition_type=sr_shift` then a diagonal uniform shift is added to S
S --> S+ϵ*identity

If `precondition_type=sr_multiplicative` then a diagonal multiplicative shift is added to S
S --> S + max(λ0*b^n,λmin)*Diagonal(diag(S)) where n is the number of the iteration.
"""
struct SR{T1,T2,T3,TP} <: Algorithm
    sr_diag_shift::T1  # cutoff
    sr_diag_mult::T2
    sr_precision::T3   #cutoff
    precondition_type::SRPreconditionType
    use_iterative::Bool
    λ0::TP
    b::TP
    λmin::TP
end

SR(T::Type=STD_REAL_PREC; ϵ=0.001, use_iterative=true, precision=10e-5,
   precondition_type=sr_shift,
   λ0=100.0, b=0.95, λmin=1e-4) = SR(T(ϵ), T(1.0), T(precision),
                                          precondition_type, use_iterative,
                                          T(λ0), T(b), T(λmin))

# -------------- Base.show extension for nice printing -------------- #
Base.show(io::IO, mm::MIME"text/plain", sr::SR) = begin
    print(io, "SR{$(typeof(sr.sr_diag_shift)), $(typeof(sr.λ0))}\n")
    print(io, "\tActive preconditioner: $(sr.precondition_type)\n")
    if sr.precondition_type == sr_multiplicative
        print(io, "\t - λ0   : $(sr.λ0)\n")
        print(io, "\t - b    : $(sr.b)\n")
        print(io, "\t - λmin : $(sr.λmin)")
    else
        print(io, "\t - ϵ : $(sr.sr_diag_shift)")
    end
end

#Base.show(io::IO, sr::SR) = print(io, "(",bs.σ_row,", ", bs.σ_col, ")")

################################################################################
######  Structure holding information computed at the end of a sampling   ######
################################################################################
"""
The SREvaluation structure holds the evaluation of the Network, L, the
vector of generalized forces F acting on it, and the change of basis matrix
S. See S.Sorella et al.
"""
mutable struct SREvaluation{TL,TF,TS} <: EvaluatedAlgorithm
    L::TL
    F::TF
    S::TS

    # Individual values to compute statistical correlators
    LVals::Vector
end

function SREvaluation(net::NeuralNetwork)
    wt = grad_cache(net) # this should be weight_type
    WT = weight_type(net)
    T = out_type(net)

    F = Tuple([zeros(WT,size(w)) for w=wt.tuple_all_weights])
    S = Tuple([zeros(WT,size(w*w')) for w=wt.tuple_all_weights])

    SREvaluation(zero(T),
                 F,
                 S,
                 Vector{T}())
end

EvaluatedNetwork(alg::SR, net) =
    SREvaluation(weights(net))

# Utility method utilised to accumulate results on a single variable
function add!(acc::SREvaluation, o::SREvaluation)
   acc.L  += o.L
   acc.F .+= o.F
   acc.S .+= o.S

   append!(acc.LVals, o.LVals)
end

function precondition!(∇x, params::SR, data::SREvaluation, iter_n)
    ϵ = params.sr_diag_shift

    success = true

    for (Δw, S, F) = zip(∇x, data.S, data.F)
        # new matrix
        if params.precondition_type == sr_none
            Sprecond = S
        elseif params.precondition_type == sr_shift
            Sprecond = S + convert(eltype(S), ϵ)*I
        elseif params.precondition_type == sr_multiplicative
            λ0 = params.λ0; b = params.b; λmin = params.λmin;
            #vv = eigvals(S)
            #println("regol $(max(λ0*b^iter_n, λmin)) --> $(minimum(vv)), a $(maximum(vv))")
            λ = convert(eltype(S),max(λ0*b^iter_n, λmin))
            Sprecond = S + λ*Diagonal(diag(S))

        end

        if !params.use_iterative
            try
                Δw .= pinv(Sprecond)*F
            catch err
                println("Could not invert: $err")
                Δw .= 0.0
                success = false
            end
        else
            x, hist = minresqlp(Sprecond, F, maxiter=size(S,2)*10, log=true, verbose=false, tol=params.sr_precision)
            #x, hist = cg(S.+ ϵ*I, F, maxiter=size(S,2)*10, log=true, verbose=true, tol=10e-10)
            Δw .= x
            add_iters = 1
            while !hist.isconverged
                println("minresqlp not conerged. Additional $(size(S,2)*10) iters for the $add_iters time.")
                x, hist = minresqlp(ΔW, Sprecond, F, maxiter=size(S,2)*10, log=true)
                Δw .= x
                add_iters += 1
                add_iters > 5 && break
            end
            if add_iters > 5 && !hist.isconverged
                success = false
                Δw .= 0.0
            end
        end
    end
    return success
end


text:
using ADFPEPS
using ADFPEPS: parity_conserving
using CUDA
using LinearAlgebra
using Random
using Test
using VUMPS:num_grad
using ADFPEPS:HamiltonianModel
using Zygote
CUDA.allowscalar(false)

@testset "parity_conserving" for atype in [Array], dtype in [ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D = 2
    T = atype(rand(dtype,D,D,D))
    parity_conserving(T)
end


text:
## Whittle.jl : implementation of Whittle covariance function

## Whittle ##
struct Whittle{T} <: IsotropicCovarianceStructure{T}
    λ::T
    σ::T
    p::T

    function Whittle{T}(λ::T, σ::T, p::T) where T
        λ > 0 || throw(DomainError(λ, "correlation length λ of Whittle covariance cannot be negative or zero"))
        σ > 0 || throw(DomainError(σ, "marginal standard deviation σ of Whittle covariance cannot be negative or zero"))
        p >= 1 || throw(DomainError(p, "in p-norm, p must be greater than or equal to 1"))
        isinf(p) && throw(DomainError(p, "in p-norm, p cannot be infinity"))

        new{T}(λ, σ, p)
    end
end

"""
    Whittle(λ, [σ = 1], [p = 2])

Whittle covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm, defined as

``C(x, y) = σ \\displaystyle\\frac{ρ}{λ} K₁\\left(\\frac{ρ}{λ}\\right)``

with ``ρ = ||x-y||_p``.

# Examples
```jldoctest
julia> Whittle(0.1)
Whittle (λ=0.1, σ=1.0, p=2.0)

julia> Whittle(1.0, σ=2)
Whittle (λ=1.0, σ=2.0, p=2.0)

```
See also: [`Exponential`](@ref), [`Linear`](@ref), [`Spherical`](@ref), [`Gaussian`](@ref), [`SquaredExponential`](@ref), [`Matern`](@ref)
"""
Whittle(λ::Real; σ::Real=1.0, p::Real=2) = Whittle{promote_type(typeof(λ),typeof(σ),typeof(p))}(promote(λ, σ, p)...)

# evaluate Whittle covariance
function apply(w::Whittle, x::Real)
    if iszero(x)
        float(one(x))
    else
        ρ = x / w.λ
        ρ * besselk(1, ρ)
    end
end

# short name
shortname(::Whittle) = "Whittle"


text:
# # Example 4: Pulse Parametrization

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`simple_state_to_state.ipynb`](@__NBVIEWER_ROOT_URL__/examples/state_to_state_parametrization.ipynb).

#md # ``\gdef\op#1{\hat{#1}}``
#md # ``\gdef\init{\text{init}}``
#md # ``\gdef\tgt{\text{tgt}}``

#nb # $
#nb # \newcommand{tr}[0]{\operatorname{tr}}
#nb # \newcommand{diag}[0]{\operatorname{diag}}
#nb # \newcommand{abs}[0]{\operatorname{abs}}
#nb # \newcommand{pop}[0]{\operatorname{pop}}
#nb # \newcommand{aux}[0]{\text{aux}}
#nb # \newcommand{opt}[0]{\text{opt}}
#nb # \newcommand{tgt}[0]{\text{tgt}}
#nb # \newcommand{init}[0]{\text{init}}
#nb # \newcommand{lab}[0]{\text{lab}}
#nb # \newcommand{rwa}[0]{\text{rwa}}
#nb # \newcommand{bra}[1]{\langle#1\vert}
#nb # \newcommand{ket}[1]{\vert#1\rangle}
#nb # \newcommand{Bra}[1]{\left\langle#1\right\vert}
#nb # \newcommand{Ket}[1]{\left\vert#1\right\rangle}
#nb # \newcommand{Braket}[2]{\left\langle #1\vphantom{#2}\mid{#2}\vphantom{#1}\right\rangle}
#nb # \newcommand{op}[1]{\hat{#1}}
#nb # \newcommand{Op}[1]{\hat{#1}}
#nb # \newcommand{dd}[0]{\,\text{d}}
#nb # \newcommand{Liouville}[0]{\mathcal{L}}
#nb # \newcommand{DynMap}[0]{\mathcal{E}}
#nb # \newcommand{identity}[0]{\mathbf{1}}
#nb # \newcommand{Norm}[1]{\lVert#1\rVert}
#nb # \newcommand{Abs}[1]{\left\vert#1\right\vert}
#nb # \newcommand{avg}[1]{\langle#1\rangle}
#nb # \newcommand{Avg}[1]{\left\langle#1\right\rangle}
#nb # \newcommand{AbsSq}[1]{\left\vert#1\right\vert^2}
#nb # \newcommand{Re}[0]{\operatorname{Re}}
#nb # \newcommand{Im}[0]{\operatorname{Im}}
#nb # $

# This example illustrates the parametrization of control pulses as a
# form of amplitude constraint.

using QuantumControl
using QuantumControl.Shapes: flattop
using Krotov: SquareParametrization, TanhParametrization, TanhSqParametrization, LogisticParametrization, LogisticSqParametrization
using LinearAlgebra

#jl using Test

using PyPlot
matplotlib.use("Agg")

# ## Parametrizations

# ## Symmetric Bounded Controls

function plot_symmetric_parametrization_comparison()
    fig, axs = matplotlib.pyplot.subplots(figsize=(16, 3), ncols=3)

    u_vals = collect(range(-3, 3, length=101))
    ϵ_vals = collect(range(-1, 1, length=101))

    ϵ_min = -1.0
    ϵ_max = 1.0

    axs[1].plot(u_vals, u_vals, "--", color="black")
    axs[1].plot(u_vals, TanhParametrization(ϵ_min, ϵ_max).epsilon_of_u.(u_vals), label="Tanh")
    axs[1].plot(u_vals, LogisticParametrization(ϵ_min, ϵ_max).epsilon_of_u.(u_vals), label="Logistic(k=1)")
    axs[1].plot(u_vals, LogisticParametrization(ϵ_min, ϵ_max, k=4).epsilon_of_u.(u_vals), label="Logistic(k=4)")
    axs[1].set_ylim(-1.2, 1.2)
    axs[1].set_xlabel("u")
    axs[1].set_ylabel("ϵ")

    axs[2].plot(ϵ_vals, ϵ_vals, "--", color="black")
    axs[2].plot(ϵ_vals, TanhParametrization(ϵ_min, ϵ_max).u_of_epsilon.(ϵ_vals), label="Tanh")
    axs[2].plot(ϵ_vals, LogisticParametrization(ϵ_min, ϵ_max).u_of_epsilon.(ϵ_vals), label="Logistic(k=1)")
    axs[2].plot(ϵ_vals, LogisticParametrization(ϵ_min, ϵ_max, k=4).u_of_epsilon.(ϵ_vals), label="Logistic(k=4)")
    axs[2].set_ylim(-3, 3)
    axs[2].set_ylabel("u")
    axs[2].set_xlabel("ϵ")
    axs[2].legend()

    axs[3].plot(u_vals, [1.0 for _ in u_vals], "--", color="black")
    axs[3].plot(u_vals, TanhParametrization(ϵ_min, ϵ_max).de_du_derivative.(u_vals), label="Tanh")
    axs[3].plot(u_vals, LogisticParametrization(ϵ_min, ϵ_max).de_du_derivative.(u_vals), label="Logistic(k=1)")
    axs[3].plot(u_vals, LogisticParametrization(ϵ_min, ϵ_max, k=4).de_du_derivative.(u_vals), label="Logistic(k=4)")
    axs[3].set_ylim(0, 2)
    axs[3].set_xlabel("u")
    axs[3].set_ylabel("∂ϵ/∂u")

    return fig
end

#!jl plot_symmetric_parametrization_comparison()

# ## Positive (Bounded) Controls

function plot_positive_parametrization_comparison()
    fig, axs = matplotlib.pyplot.subplots(figsize=(16, 3), ncols=3)

    u_vals = collect(range(-3, 3, length=101))
    ϵ_vals = collect(range(0, 1, length=101))
    ϵ_max = 1.0

    axs[1].plot(u_vals, abs.(u_vals), "--", color="black")
    axs[1].plot(u_vals, TanhSqParametrization(ϵ_max).epsilon_of_u.(u_vals), label="TanhSq")
    axs[1].plot(u_vals, LogisticSqParametrization(ϵ_max).epsilon_of_u.(u_vals), label="LogisticSq(k=1)")
    axs[1].plot(u_vals, LogisticSqParametrization(ϵ_max, k=4.0).epsilon_of_u.(u_vals), label="LogisticSq(k=4)")
    axs[1].plot(u_vals, SquareParametrization().epsilon_of_u.(u_vals), label="Square")
    axs[1].set_ylim(0, 1.2)
    axs[1].set_xlabel("u")
    axs[1].set_ylabel("ϵ")

    axs[2].plot(ϵ_vals, ϵ_vals, "--", color="black")
    axs[2].plot(ϵ_vals, TanhSqParametrization(ϵ_max).u_of_epsilon.(ϵ_vals), label="TanhSq")
    axs[2].plot(ϵ_vals, LogisticSqParametrization(ϵ_max).u_of_epsilon.(ϵ_vals), label="LogisticSq(k=1)")
    axs[2].plot(ϵ_vals, LogisticSqParametrization(ϵ_max, k=4.0).u_of_epsilon.(ϵ_vals), label="LogisticSq(k=4)")
    axs[2].plot(ϵ_vals, SquareParametrization().u_of_epsilon.(ϵ_vals), label="Square")
    axs[2].set_ylim(0, 3)
    axs[2].set_ylabel("u")
    axs[2].set_xlabel("ϵ")
    axs[2].legend()

    axs[3].plot(u_vals, sign.(u_vals), "--", color="black")
    axs[3].plot(u_vals, TanhSqParametrization(ϵ_max).de_du_derivative.(u_vals), label="TanhSq")
    axs[3].plot(u_vals, LogisticSqParametrization(ϵ_max).de_du_derivative.(u_vals), label="LogisticSq(k=1)")
    axs[3].plot(u_vals, LogisticSqParametrization(ϵ_max, k=4.0).de_du_derivative.(u_vals), label="LogisticSq(k=4)")
    axs[3].plot(u_vals, SquareParametrization().de_du_derivative.(u_vals), label="Square")
    axs[3].set_ylim(-2, 2)
    axs[3].set_xlabel("u")
    axs[3].set_ylabel("∂ϵ/∂u")

    return fig
end

#!jl plot_positive_parametrization_comparison()

# ## Two-level Hamiltonian

# We consider the Hamiltonian $\op{H}_{0} = - \frac{\omega}{2} \op{\sigma}_{z}$, representing
# a simple qubit with energy level splitting $\omega$ in the basis
# $\{\ket{0},\ket{1}\}$. The control field $\epsilon(t)$ is assumed to couple via
# the Hamiltonian $\op{H}_{1}(t) = \epsilon(t) \op{\sigma}_{x}$ to the qubit,
# i.e., the control field effectively drives transitions between both qubit
# states.
#
# We we will use

ϵ(t) = 0.2 * flattop(t, T=5, t_rise=0.3, func=:blackman);

#-
"""Two-level-system Hamiltonian."""
function hamiltonian(Ω=1.0, ϵ=ϵ)
    σ̂_z = ComplexF64[1 0; 0 -1];
    σ̂_x = ComplexF64[0 1; 1  0];
    Ĥ₀ = -0.5 * Ω * σ̂_z
    Ĥ₁ = σ̂_x
    return (Ĥ₀, (Ĥ₁, ϵ))
end
;
#-

H = hamiltonian();
#jl @test length(H) == 2

# The control field here switches on from zero at $t=0$ to it's maximum amplitude
# 0.2 within the time period 0.3 (the switch-on shape is half a [Blackman pulse](https://en.wikipedia.org/wiki/Window_function#Blackman_window)).
# It switches off again in the time period 0.3 before the
# final time $T=5$). We use a time grid with 500 time steps between 0 and $T$:

tlist = collect(range(0, 5, length=500));

#-
function plot_control(pulse::Vector, tlist)
    fig, ax = matplotlib.pyplot.subplots(figsize=(6, 3))
    ax.plot(tlist, pulse)
    ax.set_xlabel("time")
    ax.set_ylabel("amplitude")
    return fig
end

plot_control(ϵ::T, tlist) where T<:Function =
    plot_control([ϵ(t) for t in tlist], tlist)

#!jl plot_control(H[2][2], tlist)

# ## Optimization target

# The `krotov` package requires the goal of the optimization to be described by a
# list of `Objective` instances. In this example, there is only a single
# objective: the state-to-state transfer from initial state $\ket{\Psi_{\init}} =
# \ket{0}$ to the target state $\ket{\Psi_{\tgt}} = \ket{1}$, under the dynamics
# of the Hamiltonian $\op{H}(t)$:

function ket(label)
    result = Dict(
        "0" => Vector{ComplexF64}([1, 0]),
        "1" => Vector{ComplexF64}([0, 1]),
    )
    return result[string(label)]
end
;

#-
#jl @test dot(ket(0), ket(1)) ≈ 0
#-

objectives = [
    Objective(initial_state=ket(0), generator=H, target_state=ket(1))
]

#-
#jl @test length(objectives) == 1
#-

# ## Square-parametrization for positive pulses

problem = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ  => Dict(
            :lambda_a => 5,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => SquareParametrization(),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    chi=QuantumControl.Functionals.chi_ss!,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence= res -> begin (
            (res.J_T < 1e-3)
            && (res.converged = true)
            && (res.message="J_T < 10⁻³")
        ) end
);
#-
opt_result_positive = optimize(problem, method=:krotov);
#-
opt_result_positive

# We can plot the optimized field:

#-
#!jl plot_control(opt_result_positive.optimized_controls[1], tlist)
#-

#-
#jl @test minimum(opt_result_positive.optimized_controls[1]) ≥ 0.0
#jl @test minimum(opt_result_positive.optimized_controls[1]) < 1e-16
#jl @test maximum(opt_result_positive.optimized_controls[1]) > 0.0
#-

# ## Tanh-Square-Parametrization for positive amplitude-constrained pulses

problem_tanhsq = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ  => Dict(
            :lambda_a => 10,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => TanhSqParametrization(3),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    chi=QuantumControl.Functionals.chi_ss!,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence= res -> begin (
            (res.J_T < 1e-3)
            && (res.converged = true)
            && (res.message="J_T < 10⁻³")
        ) end
);
#-
opt_result_tanhsq = optimize(problem_tanhsq, method=:krotov);
#-
opt_result_tanhsq

# We can plot the optimized field:

#-
#!jl plot_control(opt_result_tanhsq.optimized_controls[1], tlist)
#-

#-
#jl @test minimum(opt_result_tanhsq.optimized_controls[1]) ≥ 0.0
#jl @test minimum(opt_result_tanhsq.optimized_controls[1]) < 1e-16
#jl @test maximum(opt_result_tanhsq.optimized_controls[1]) > 0.0
#jl @test maximum(opt_result_tanhsq.optimized_controls[1]) < 3.0
#-

# ## Logistic-Square-Parametrization for positive amplitude-constrained pulses

problem_logisticsq = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ  => Dict(
            :lambda_a => 1,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => LogisticSqParametrization(3, k=1.0),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    chi=QuantumControl.Functionals.chi_ss!,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence= res -> begin (
            (res.J_T < 1e-3)
            && (res.converged = true)
            && (res.message="J_T < 10⁻³")
        ) end
);
#-
opt_result_logisticsq = optimize(problem_logisticsq, method=:krotov);
# We can plot the optimized field:

#-
#!jl plot_control(opt_result_logisticsq.optimized_controls[1], tlist)
#-

#-
#jl @test minimum(opt_result_logisticsq.optimized_controls[1]) ≥ 0.0
#jl @test minimum(opt_result_logisticsq.optimized_controls[1]) < 1e-16
#jl @test maximum(opt_result_logisticsq.optimized_controls[1]) > 0.0
#jl @test maximum(opt_result_logisticsq.optimized_controls[1]) < 3.0
#-
# ## Tanh-parametrization for amplitude-constrained pulses

problem_tanh = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ  => Dict(
            :lambda_a => 1,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => TanhParametrization(-0.5, 0.5),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    chi=QuantumControl.Functionals.chi_ss!,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence= res -> begin (
            (res.J_T < 1e-3)
            && (res.converged = true)
            && (res.message="J_T < 10⁻³")
        ) end
);
#-
opt_result_tanh = optimize(problem_tanh, method=:krotov);
#-
#!jl plot_control(opt_result_tanh.optimized_controls[1], tlist)
#-
#jl @test minimum(opt_result_tanh.optimized_controls[1]) > -0.5
#jl @test maximum(opt_result_tanh.optimized_controls[1]) < 0.5
#-


text:
using OpenQuantumTools, Test, Random
using OrdinaryDiffEq

H = DenseHamiltonian([(s) -> 1.0], -[σz], unit=:ħ)
u0 = PauliVec[1][1]
coupling = ConstantCouplings(["Z"], unit=:ħ)
bath = EnsembleFluctuator([1.0, 2.0], [2.0, 1.0])
annealing = Annealing(H, u0, coupling=coupling, bath=bath)
tf = 2.0

prob = build_ensembles(annealing, tf, :stochastic)
Random.seed!(1234)
sol1 = solve(prob, Tsit5(), EnsembleSerial(), trajectories=1, save_everystep=false)

Random.seed!(1234)
sol2 = solve(prob, Tsit5(), EnsembleSerial(), trajectories=1, save_everystep=false)

@test sol1[1][end] ≈ sol2[1][end] atol = 1e-6 rtol = 1e-6

sol3 = solve(prob, Tsit5(), EnsembleSerial(), trajectories=1, save_everystep=false)

@test !(isapprox(sol1[1][end], sol3[1][end], atol=1e-6, rtol=1e-6))

max_stars_repo_name: huangjingzheng/OpenQuantumTools

max_issues_repo_name: huangjingzheng/OpenQuantumTools

max_forks_repo_name: huangjingzheng/OpenQuantumTools


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function set_par!(θ_bound, θ_unb, par, opt_transf, MIN, MAX, par_ind, par_size, prior_opt, apriori_rejection)

# ----------------------------------------------------------------------------------------------------------------------
# Update the parameters of the state-space model
# ----------------------------------------------------------------------------------------------------------------------

     if sum(ismissing.(θ_bound)) > 0 || sum(isinf.(θ_bound)) > 0 || sum(isnan.(θ_bound)) > 0 || sum(.~isreal.(θ_bound)) > 0
          apriori_rejection[1] = 1;

     else

          # ------------------------------------------------------------------------------------------------------------
          # Set state-space parameters and par.logprior
          # ------------------------------------------------------------------------------------------------------------

          iend         = 0;
          par.logprior = 0;

          # Observation equations

          if par_size.R > 0
               par.R[par_ind.R .== true] = θ_bound[1:par_size.R];
               par.logprior              = par.logprior + sum(logpdf.(prior_opt.IG, par.R[par_ind.R .== true]));
               iend                      = par_size.R;
          end

          if par_size.d > 0
               par.d[par_ind.d .== true] = θ_bound[iend+1:iend+par_size.d];
               par.logprior              = par.logprior + sum(logpdf.(prior_opt.N, par.d[par_ind.d .== true]));
               iend                      = iend+par_size.d;
          end

          if par_size.Z > 0
               par.Z[par_ind.Z .== true] = θ_bound[iend+1:iend+par_size.Z];
               par.logprior              = par.logprior + sum(logpdf.(prior_opt.N, par.Z[par_ind.Z .== true]));
               iend                      = iend+par_size.Z;
          end

          # Transition equations

          if par_size.Q > 0
               par.Q[par_ind.Q .== true] = θ_bound[iend+1:iend+par_size.Q];
               par.logprior              = par.logprior + sum(logpdf.(prior_opt.IG, par.Q[par_ind.Q .== true]));
               iend                      = iend+par_size.Q;
          end

          if par_size.c > 0
               par.c[par_ind.c .== true] = θ_bound[iend+1:iend+par_size.c];
               par.logprior              = par.logprior + sum(logpdf.(prior_opt.N, par.c[par_ind.c .== true]));
               iend                      = iend+par_size.c;
          end

          if par_size.T > 0 || par_size.λ > 0 || par_size.ρ > 0

               # Set T
               par.T[par_ind.T .== true] = θ_bound[iend+1:iend+par_size.T];
               par.logprior              = par.logprior + sum(logpdf.(prior_opt.N, par.T[par_ind.T .== true]));
               iend                      = iend+par_size.T;

               # Trigonometric states: update T, λ and ρ. Adjust Q and P¹
               if par_size.λ > 0 || par_size.ρ > 0
                    par.λ[par_ind.λ .== true] = θ_bound[iend+1:iend+par_size.λ];
                    par.ρ[par_ind.ρ .== true] = θ_bound[iend+par_size.λ+1:end];
                    par.logprior              = par.logprior + prior_opt.λ + prior_opt.ρ;

                    bool_trig = ((par_ind.λ .== true) .+ (par_ind.ρ .== true)) .> 0;
                    find_trig = findall(bool_trig);

                    for i=1:sum(bool_trig)
                         j = find_trig[i];

                         par.T[j:j+1, j:j+1]  = par.ρ[j]*[[cos(par.λ[j]) sin(par.λ[j])]; [-sin(par.λ[j]) cos(par.λ[j])]];
                         par.Q[j+1, j+1]      = par.Q[j, j];
                         par.P¹[j:j+1, j:j+1] = (par.Q[j, j] ./ (1-par.ρ[j]^2))*Array{Float64}(Matrix(I,2,2));
                    end
               end
          end

          # ------------------------------------------------------------------------------------------------------------
          # par.loglikelihood and par.logposterior
          # ------------------------------------------------------------------------------------------------------------

          kalman_diffuse!(par, 1, 0, 0); # estimate loglikelihood
          par.logprior     = par.logprior + sum(get_logjacobian(θ_unb, MIN, MAX, opt_transf));
          par.logposterior = par.loglik + par.logprior;
     end
end


text:
function depth_image_to_point_cloud(
        depth_map::Matrix, intrinsics::CameraIntrinsics;
        stride_x=1, stride_y=1, flatten=true, cut_far_plane=true)

    height = intrinsics.height
    width = intrinsics.width
    cx = intrinsics.cx
    cy = intrinsics.cy
    fx = intrinsics.fx
    fy = intrinsics.fy

    if size(depth_map) != (height, width)
        error("expected depth map that is (height x width) matrix")
    end

    x = depth_map[1:stride_y:end, 1:stride_x:end] .* (collect(1:stride_x:width) .- cx .- 0.5)' ./ fx
    y = depth_map[1:stride_y:end, 1:stride_x:end] .* (collect(1:stride_y:height) .- cy .- 0.5) ./ fy
    z = depth_map[1:stride_y:end, 1:stride_x:end]
    new_height, new_width = size(x)
    point_cloud = Array{Float64}(undef, new_height, new_width, 3)
    point_cloud[:,:,1] = x
    point_cloud[:,:,2] = y
    point_cloud[:,:,3] = z

    if !flatten && cut_far_plane
        println("Warning!!! cut_far_plane is being ignored because flatten is false.")
    end

    if flatten
        point_cloud = flatten_point_cloud(point_cloud)
        if cut_far_plane
            point_cloud = point_cloud[:, point_cloud[3,:] .< (intrinsics.far - 0.1)]
            if size(point_cloud)[2] == 0
                point_cloud = nothing
            end
        end
    end

    return point_cloud
end

function point_cloud_to_pixel_coordinates(point_cloud::Matrix{Float64}, intrinsics::CameraIntrinsics)
    if size(point_cloud)[1] != 3
        error("expected an 3 x n matrix")
    end

    point_cloud_normalized =  point_cloud ./ point_cloud[3,:]'

    temp = point_cloud_normalized[1:2, :] .* [intrinsics.fx, intrinsics.fy]
    temp = temp .+ [intrinsics.cx + 0.5, intrinsics.cy + 0.5]
    pixel_coords = round.(Int, temp)
    return pixel_coords
end

function flatten_point_cloud(point_cloud::Array{Float64,3})
    (height, width, k) = size(point_cloud)
    if k != 3
        error("expected (height x width x 3) array")
    end
    flat_point_cloud = reshape(permutedims(point_cloud, (3, 1, 2)), (3, width * height))
    @assert size(flat_point_cloud) == (3, width * height)
    return flat_point_cloud
end

function get_points_in_frame_b(points_in_frame_a::Matrix{<:Real}, b_relative_to_a::Pose)
    if size(points_in_frame_a)[1] != 3
        error("expected an 3 x n matrix")
    end
    n = size(points_in_frame_a[2])
    return b_relative_to_a.orientation' * (points_in_frame_a .- b_relative_to_a.pos)
end

function move_points_to_frame_b(points_in_frame_a::Matrix{<:Real}, b_relative_to_a::Pose)
    if size(points_in_frame_a)[1] != 3
        error("expected an 3 x n matrix")
    end
    n = size(points_in_frame_a[2])
    return b_relative_to_a.orientation * points_in_frame_a .+ b_relative_to_a.pos
end

export flatten_point_cloud, depth_image_to_point_cloud, get_points_in_frame_b, move_points_to_frame_b


text:
using Mux
using Mustache
using JSON

abstract type AbstractGrid end
update(grid::AbstractGrid) = nothing

mutable struct D3Grid{T <: AbstractGrid}
    h::Int
    w::Int
    f::T
end

process_main(gr::D3Grid) = render(
    read(joinpath(@__DIR__, "../templates/grid.html"), String),
    Dict("width" => gr.w, "height" => gr.h, "c_width" => 1000 ÷ gr.w, "c_height" => 1000 ÷ gr.h))

function run(gr::D3Grid)
    @app d3 = (
        Mux.defaults,
        page("/", req -> Dict(
            :body => process_main(gr))),
        page("/data", req -> Dict(
                :headers => ["Content-Type" => "application/json"],
                :body => JSON.json(update(gr.f))
            )),
        Mux.notfound()
    )

    serve(d3)
end


text:
# using Revise
# using Pkg
# cd("..")
# Pkg.activate(".")

using Test
using WeightedOnlineStats
import OnlineStatsBase: eachcol, eachrow
using StatsBase
import OnlineStats: Extrema
using Statistics
using Random
using MultivariateStats

Random.seed!(123)
l = 1000
x = rand(l);
xmis = convert(Array{Union{Float64, Missing}}, x);
x2 = rand(l, 5)
x2mis = convert(Array{Union{Float64, Missing}}, x2)
x2mis[end, 1] = missing
x2mis[end - 1, 1] = missing
w = rand(l);
wmis = convert(Array{Union{Float64, Missing}}, w);
wmis[end] = missing
wmis[end - 2] = missing

Base.isapprox(x::Tuple, y::Tuple) =
    reduce(&, map((a, b) -> a ≈ b, x, y))

missing_to_nan(x::Array{Union{Missing, T}}) where T = map(y -> y === missing ? T(NaN) : y, x)

function test_fit(T::Type{<: WeightedOnlineStats.WeightedOnlineStat{S}},
                  data, weights, unpack_fun, jfun) where S
    l = length(data)
    @assert l == length(weights)

    o = T()
    for (xi, wi) in zip(data, weights)
        fit!(o, xi, wi)
    end

    @test unpack_fun(o) ≈ jfun(data, weights)
    @test eltype(unpack_fun(o)) == eltype(o)
end

include("test_hist.jl")
include("test_sum.jl")
include("test_mean.jl")
include("test_var.jl")
include("test_cov.jl")
include("test_pca.jl")


text:
# TODO optimize ChoiceAtTrace using type parameters

struct ChoiceAtTrace <: Trace
    gen_fn::GenerativeFunction # the ChoiceAtCombinator (not the kernel)
    value::Any
    key::Any
    kernel_args::Tuple
    score::Float64
end

get_args(trace::ChoiceAtTrace) = (trace.kernel_args..., trace.key)
get_retval(trace::ChoiceAtTrace) = trace.value
get_score(trace::ChoiceAtTrace) = trace.score
get_gen_fn(trace::ChoiceAtTrace) = trace.gen_fn

struct ChoiceAtChoiceMap{T,K} <: ChoiceMap
    key::K
    value::T
end

get_choices(trace::ChoiceAtTrace) = ChoiceAtChoiceMap(trace.key, trace.value)
Base.isempty(::ChoiceAtChoiceMap) = false
function get_address_schema(::Type{T}) where {T<:ChoiceAtChoiceMap}
    SingleDynamicKeyAddressSchema()
end
get_value(choices::ChoiceAtChoiceMap, addr::Pair) = _get_value(choices, addr)
has_value(choices::ChoiceAtChoiceMap, addr::Pair) = _has_value(choices, addr)
function get_value(choices::ChoiceAtChoiceMap{T,K}, addr::K) where {T,K}
    choices.key == addr ? choices.value : throw(KeyError(choices, addr))
end
get_submaps_shallow(choices::ChoiceAtChoiceMap) = ()
get_values_shallow(choices::ChoiceAtChoiceMap) = ((choices.key, choices.value),)

struct ChoiceAtCombinator{T,K} <: GenerativeFunction{T, ChoiceAtTrace}
    dist::Distribution{T}
end

accepts_output_grad(gen_fn::ChoiceAtCombinator) = has_output_grad(gen_fn.dist)

# TODO
# accepts_output_grad is true if the return value is dependent on the 'gradient source elements'
# if the random choice itself is not a 'gradient source element' then it is independent (false)
# if the random choice is a 'gradient source element', then the return value is dependent (true)
# we will consider the random choice as a gradient source element if the
# distribution has has_output_grad = true)

function choice_at(dist::Distribution{T}, ::Type{K}) where {T,K}
    ChoiceAtCombinator{T,K}(dist)
end

unpack_choice_at_args(args) = (args[end], args[1:end-1])

function assess(gen_fn::ChoiceAtCombinator{T,K}, args::Tuple, choices::ChoiceMap) where {T,K}
    local key::K
    local value::T
    (key, kernel_args) = unpack_choice_at_args(args)
    value = get_value(choices, key)
    weight = logpdf(gen_fn.dist, value, kernel_args...)
    (weight, value)
end

function propose(gen_fn::ChoiceAtCombinator{T,K}, args::Tuple) where {T,K}
    local key::K
    local value::T
    (key, kernel_args) = unpack_choice_at_args(args)
    value = random(gen_fn.dist, kernel_args...)
    score = logpdf(gen_fn.dist, value, kernel_args...)
    choices = ChoiceAtChoiceMap(key, value)
    (choices, score, value)
end

function simulate(gen_fn::ChoiceAtCombinator, args::Tuple)
    (key, kernel_args) = unpack_choice_at_args(args)
    value = random(gen_fn.dist, kernel_args...)
    score = logpdf(gen_fn.dist, value, kernel_args...)
    ChoiceAtTrace(gen_fn, value, key, kernel_args, score)
end

function generate(gen_fn::ChoiceAtCombinator{T,K}, args::Tuple, choices::ChoiceMap) where {T,K}
    local key::K
    local value::T
    (key, kernel_args) = unpack_choice_at_args(args)
    constrained = has_value(choices, key)
    value = constrained ? get_value(choices, key) : random(gen_fn.dist, kernel_args...)
    score = logpdf(gen_fn.dist, value, kernel_args...)
    trace = ChoiceAtTrace(gen_fn, value, key, kernel_args, score)
    weight = constrained ? score : 0.
    (trace, weight)
end

function project(trace::ChoiceAtTrace, selection::Selection)
    (trace.key in selection) ? trace.score : 0.
end

function update(trace::ChoiceAtTrace, args::Tuple, argdiffs::Tuple,
                choices::ChoiceMap)
    (key, kernel_args) = unpack_choice_at_args(args)
    key_changed = (key != trace.key)
    constrained = has_value(choices, key)
    if key_changed && constrained
        new_value = get_value(choices, key)
        discard = ChoiceAtChoiceMap(trace.key, trace.value)
    elseif !key_changed && constrained
        new_value = get_value(choices, key)
        discard = ChoiceAtChoiceMap(key, trace.value)
    elseif !key_changed && !constrained
        new_value = trace.value
        discard = EmptyChoiceMap()
    else
        error("New address $key not constrained in update")
    end
    new_score = logpdf(trace.gen_fn.dist, new_value, kernel_args...)
    new_trace = ChoiceAtTrace(trace.gen_fn, new_value, key, kernel_args, new_score)
    weight = new_score - trace.score
    (new_trace, weight, UnknownChange(), discard)
end

function regenerate(trace::ChoiceAtTrace, args::Tuple, argdiffs::Tuple,
                    selection::Selection)
    (key, kernel_args) = unpack_choice_at_args(args)
    key_changed = (key != trace.key)
    selected = key in selection
    if !key_changed && selected
        new_value = random(trace.gen_fn.dist, kernel_args...)
    elseif !key_changed && !selected
        new_value = trace.value
    elseif key_changed && !selected
        new_value = random(trace.gen_fn.dist, kernel_args...)
    else
        error("Cannot select new address $key in regenerate")
    end
    new_score = logpdf(trace.gen_fn.dist, new_value, kernel_args...)
    if !key_changed && selected
        weight = 0.
    elseif !key_changed && !selected
        weight = new_score - trace.score
    elseif key_changed && !selected
        weight = 0.
    end
    new_trace = ChoiceAtTrace(trace.gen_fn, new_value, key, kernel_args, new_score)
    (new_trace, weight, UnknownChange())
end

function choice_gradients(trace::ChoiceAtTrace, selection::Selection, retval_grad)
    if retval_grad != nothing && !has_output_grad(trace.gen_fn.dist)
        error("return value gradient not accepted but one was provided")
    end
    kernel_arg_grads = logpdf_grad(trace.gen_fn.dist, trace.value, trace.kernel_args...)
    if trace.key in selection
        value_choices = ChoiceAtChoiceMap(trace.key, trace.value)
        choice_grad = kernel_arg_grads[1]
        if choice_grad == nothing
            error("gradient not available for selected choice")
        end
        if retval_grad != nothing
            choice_grad += retval_grad
        end
        gradient_choices = ChoiceAtChoiceMap(trace.key, choice_grad)
    else
        value_choices = EmptyChoiceMap()
        gradient_choices = EmptyChoiceMap()
    end
    input_grads = (kernel_arg_grads[2:end]..., nothing)
    (input_grads, value_choices, gradient_choices)
end

function accumulate_param_gradients!(trace::ChoiceAtTrace, retval_grad)
    if retval_grad != nothing && !has_output_grad(trace.gen_fn.dist)
        error("return value gradient not accepted but one was provided")
    end
    kernel_arg_grads = logpdf_grad(trace.gen_fn.dist, trace.value, trace.kernel_args...)
    (kernel_arg_grads[2:end]..., nothing)
end

export choice_at


text:
##

using IncrementalInference
using Test
using TensorCast
# using Statistics
# using Manifolds # should be done within regular exports

# using Pkg
# Pkg.resolve()

##

@testset "test Mixture sampling" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)

mp = Mixture(Prior, (Normal(), Normal(10,1)),(1/2,1/2) )
addFactor!(fg, [:x0], mp)

##

pts_ = approxConv(fg, :x0f1, :x0)
@cast pts[i,j] := pts_[j][i]

N = size(pts,2)
@test 0.2*N < sum( -5 .< pts .< 5 )
@test 0.2*N < sum( 5 .< pts .< 15 )
@test sum( 15 .< pts  ) < 0.1*N
@test sum( pts .< -5 ) < 0.1*N
@test sum( 3 .< pts .< 7 ) < 0.1*N

# using KernelDensityEstimatePlotting, Gadfly
# Gadfly.set_default_plot_size(25cm,20cm)

# plotKDE(kde!(pts))

##

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg, [:x0], Prior(Normal()), graphinit=false)
initManual!(fg, :x0, [zeros(1) for _ in 1:100])

mlr = Mixture(LinearRelative, (Normal(), Normal(10,1)),(1/2,1/2) )
addFactor!(fg, [:x0;:x1], mlr, graphinit=false)

##

pts_ = approxConv(fg, :x0x1f1, :x1)
@cast pts[i,j] := pts_[j][i]

# plotKDE(kde!(pts))

##

N = size(pts,2)
@test 0.2*N < sum( -5 .< pts .< 5 )
@test 0.2*N < sum( 5 .< pts .< 15 )
@test sum( 15 .< pts  ) < 0.1*N
@test sum( pts .< -5 ) < 0.1*N
@test sum( 3 .< pts .< 7 ) < 0.1*N

##

end

@testset "test packing of Mixture" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

mp = Mixture(Prior, [Normal(); Normal(10,1)], [0.5;0.5])
f0 = addFactor!(fg, [:x0;], mp)

mr = Mixture(LinearRelative, (fancy=manikde!(ContinuousEuclid(1), [randn(1) for _ in 1:75]), naive=Normal(0,10)), [0.4;0.6])
f1 = addFactor!(fg, [:x0;:x1], mr)

##

pf0 = DFG.packFactor(fg, f0)
pf1 = DFG.packFactor(fg, f1)

# now test unpacking
fg_ = initfg();
addVariable!(fg_, :x0, ContinuousScalar)
addVariable!(fg_, :x1, ContinuousScalar)

##

f0_ = DFG.unpackFactor(fg_, pf0)
f1_ = DFG.unpackFactor(fg_, pf1)

##

# ENV["JULIA_DEBUG"] = "DistributedFactorGraphs"
@warn("Skipping pack/unpack compareFactor test for `timezone` and `zone`")
@show typeof(f1)
@show typeof(f1_)

@show  typeof(getSolverData(f1).fnc.params);
@show typeof(getSolverData(f1_).fnc.params);

@test DFG.compareFactor(f1, f1_, skip=[:components;:labels;:timezone;:zone;:vartypes;])

@test IIF._getCCW(f1).usrfnc!.components.naive == IIF._getCCW(f1).usrfnc!.components.naive

# already ManifoldKernelDensity
A = IIF._getCCW(f1).usrfnc!.components.fancy
B = IIF._getCCW(f1_).usrfnc!.components.fancy

# A = ManifoldBelief(Euclid, IIF._getCCW(f1).usrfnc!.components.fancy )
# B = ManifoldBelief(Euclid, IIF._getCCW(f1_).usrfnc!.components.fancy )

@test mmd(A,B) < 1e-6

##

end

@testset "test simple Mixture" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0], Prior(Normal(0.0,0.1)))

##

# require default ::UnitScaling constructor from all factors
mlr = Mixture(LinearRelative(I), [Normal(-1.0, 0.1), Normal(1.0, 0.1)], Categorical([0.5; 0.5]))

# test serialization while we are here
pmlr = convert(PackedMixture, mlr)
umlr = convert(Mixture, pmlr)

@test mlr.mechanics == umlr.mechanics
@test mlr.components == umlr.components
@test mlr.diversity == umlr.diversity

mlr = Mixture(LinearRelative, [Normal(-1.0, 0.1), Normal(1.0, 0.1)], Categorical([0.5; 0.5]))

addFactor!(fg, [:x0,:x1], mlr)

# To look at your factor graph
# if false
# using GraphPlot
# using DistributedFactorGraphs
# plotDFG(fg)
# end

tree = solveTree!(fg)

##

btd = getBelief(getVariable(fg, :x0))
@test isapprox(mean(getKDEfit(btd,distribution=Normal)), 0.0; atol=0.1)

@test isapprox(std(getKDEfit(btd,distribution=Normal)), 0.1; atol=0.05)

btd = getBelief(getVariable(fg, :x1))
pts_ = getPoints(btd)
@cast pts[i,j] := pts_[j][i]
pts_p = pts[pts .>= 0]
pts_n = pts[pts .< 0]

nfit_p = fit(Normal, pts_p)
@test isapprox(mean(nfit_p), 1.0; atol=0.1)
@test isapprox(std(nfit_p), 0.14; atol=0.05) #TODO confirm the correct value and tolerance

nfit_n = fit(Normal, pts_n)
@test isapprox(mean(nfit_n), -1.0; atol=0.1)
@test isapprox(std(nfit_n), 0.14; atol=0.05) #TODO confirm the correct value and tolerance

# To look at your results
# if false
#   using RoMEPlotting
#   plotKDE(fg, ls(fg))
# end

##

end


text:
@testset "GramMatrix tests" begin
    @testset "GramMatrix" begin
        @polyvar x y
        P = GramMatrix{Int}((i,j) -> i + j, [x^2, x*y, y^2])
        @test nvariables(P) == 2
        @test variables(P) == [x, y]
        Pp1 = 2x^4 + 6x^3*y + 12x^2*y^2 + 10x*y^3 + 6y^4 + 1
        @test P + 1   == Pp1
        @test P + 1.0 == Pp1
        zP = zero(typeof(P))
        @test isempty(zP.Q)
        @test zP == 0
        p = polynomial(P)
        @test coefficients(p) == [2, 6, 12, 10, 6]
        @test monomialtype(p) == typeof(x*y)
        @test monomialtype(typeof(p)) == typeof(x*y)
        @test monomials(p) == [x^4, x^3*y, x^2*y^2, x*y^3, y^4]
        for i in 1:3
            for j in 1:3
                @test P[i, j] == i + j
            end
        end
        for P in (GramMatrix((i,j) -> i * j, [y, x]),
                  GramMatrix((i,j) -> (3-i) * (3-j), monovec([y, x])),
                  GramMatrix([1 2; 2 4], [y, x]),
                  GramMatrix([4 2; 2 1], monovec([y, x])))
            @test P.Q.Q == [4, 2, 1]
            @test P.x[1] == x
            @test P.x[2] == y
        end
        P = GramMatrix((i,j) -> ((i,j) == (1,1) ? 2 : 0), [x*y, x^2, y^2])
        Q = GramMatrix([0 1; 1 0], [x^2, y^2])
        @test P == Q
        p = GramMatrix([2 3; 3 2], [x, y])
        @test polynomial(p) isa AbstractPolynomial{Int}
        @test polynomial(p, Float64) isa AbstractPolynomial{Float64}

        @test -2*x + dot(-x - x^2, 0) + GramMatrix{Int}((i,j)->1, [1,x]) == -(-x^2 - 1)
        P = GramMatrix{Int}((i,j) -> i + j, [x^2, x*y, y^2])
        p = polynomial(P)
        @test !iszero(P)
        @test iszero(P-P)
        @test iszero(P-p)
        @test iszero(p-P)
        @test P + P == P + p
        @test x * P == x * p
        @test 2 * P == 2 * p
        @test P * (2x) == (2x) * p

        @test differentiate(GramMatrix{Int}((i,j)->1, [x]), y) == 0
        @test differentiate(GramMatrix{Int}((i,j)->1, [x, y]), y, 1) == 2x + 2y
        #@inferred differentiate(GramMatrix{Int}((i,j)->1, [x, y]), y, 0) # FIXME failing at the moment

        @test GramMatrix([2 3; 3 2], [x, y]) == 2x^2 + 2y^2 + 6x*y

        mp = GramMatrix{Int}((i,j) -> i+j, [x])
        @test norm(mp) == norm(mp, 2) == 2.0

        @polyvar v[1:3]
        P = [1 2 3; 2 4 5; 3 5 6]
        p = GramMatrix(P, v)
        #@inferred p(x => ones(3))
        @test p(v => ones(3)) == 31
        #@inferred subs(p, x => ones(3))
        @test subs(p, v => ones(3)) == 31
    end
    @testset "SOSDecomposition" begin
        @polyvar x y
#       @test isempty(SOSDecomposition(typeof(x)[]))
#       ps = [1, x + y, x^2, x*y, 1 + x + x^2]
#       P = GramMatrix(SOSDecomposition(ps))
#       P.Q == [2 0 1 0 1; 0 1 0 0 0; 1 0 2 1 1; 0 0 1 1 0; 1 0 1 0 2]
#       P.x == [x^2, x*y, x, y, 1]
#       @test P == P
#       @test isapprox(GramMatrix(SOSDecomposition(P)), P)
        @test sprint(show, SOSDecomposition([x+y, x-y])) == "(x + y)^2 + (x - y)^2"
        @testset "SOSDecomposition equality" begin
            @polyvar x y
            @test !isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x+y]))
            @test !isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x+y, x+y]))
            @test isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x+y, x-y]))
            @test isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x-y, x+y+1e-8]), ztol=1e-7)
        end
    end
end


text:
##

using IncrementalInference
using Test

##

@testset "test Mixture sampling" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)

mp = Mixture(Prior, (Normal(), Normal(10,1)),(1/2,1/2) )
addFactor!(fg, [:x0], mp)

##

pts = approxConv(fg, :x0f1, :x0)

N = size(pts,2)
@test 0.2*N < sum( -5 .< pts .< 5 )
@test 0.2*N < sum( 5 .< pts .< 15 )
@test sum( 15 .< pts  ) < 0.1*N
@test sum( pts .< -5 ) < 0.1*N
@test sum( 3 .< pts .< 7 ) < 0.1*N

# using KernelDensityEstimatePlotting, Gadfly
# Gadfly.set_default_plot_size(25cm,20cm)

# plotKDE(kde!(pts))

##

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg, [:x0], Prior(Normal()), graphinit=false)
initManual!(fg, :x0, zeros(1,100))

mlr = Mixture(LinearRelative, (Normal(), Normal(10,1)),(1/2,1/2) )
addFactor!(fg, [:x0;:x1], mlr, graphinit=false)

##

pts = approxConv(fg, :x0x1f1, :x1)

# plotKDE(kde!(pts))

##

N = size(pts,2)
@test 0.2*N < sum( -5 .< pts .< 5 )
@test 0.2*N < sum( 5 .< pts .< 15 )
@test sum( 15 .< pts  ) < 0.1*N
@test sum( pts .< -5 ) < 0.1*N
@test sum( 3 .< pts .< 7 ) < 0.1*N

##

end

@testset "test packing of Mixture" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

mp = Mixture(Prior, [Normal(); Normal(10,1)], [0.5;0.5])
f0 = addFactor!(fg, [:x0;], mp)

mr = Mixture(LinearRelative, (fancy=manikde!(randn(1,75), ContinuousEuclid(1)), naive=Normal(0,10)), [0.4;0.6])
f1 = addFactor!(fg, [:x0;:x1], mr)

pf0 = DFG.packFactor(fg, f0)
pf1 = DFG.packFactor(fg, f1)

# now test unpacking
fg_ = initfg();
addVariable!(fg_, :x0, ContinuousScalar)
addVariable!(fg_, :x1, ContinuousScalar)

f0_ = DFG.unpackFactor(fg_, pf0)
f1_ = DFG.unpackFactor(fg_, pf1)

# ENV["JULIA_DEBUG"] = "DistributedFactorGraphs"
@warn("Skipping pack/unpack compareFactor test for `timezone` and `zone`")
@test DFG.compareFactor(f1, f1_, skip=[:components;:labels;:timezone;:zone])

@test IIF._getCCW(f1).usrfnc!.components.naive == IIF._getCCW(f1).usrfnc!.components.naive

A = ManifoldBelief(Euclid, IIF._getCCW(f1).usrfnc!.components.fancy )
B = ManifoldBelief(Euclid, IIF._getCCW(f1_).usrfnc!.components.fancy )

@test mmd(A,B) < 1e-6

##

end

@testset "test simple Mixture" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0], Prior(Normal(0.0,0.1)))

# require default ::UnitScaling constructor from all factors
mlr = Mixture(LinearRelative(I), [Normal(-1.0, 0.1), Normal(1.0, 0.1)], Categorical([0.5; 0.5]))

# test serialization while we are here
pmlr = convert(PackedMixture, mlr)
umlr = convert(Mixture, pmlr)

@test mlr.mechanics == umlr.mechanics
@test mlr.components == umlr.components
@test mlr.diversity == umlr.diversity

mlr = Mixture(LinearRelative, [Normal(-1.0, 0.1), Normal(1.0, 0.1)], Categorical([0.5; 0.5]))

addFactor!(fg, [:x0,:x1], mlr)

# To look at your factor graph
# if false
# using GraphPlot
# using DistributedFactorGraphs
# dfgplot(fg)
# end

tree, smt, hist = solveTree!(fg)

##

btd = getBelief(getVariable(fg, :x0))
@test isapprox(mean(getKDEfit(btd,distribution=Normal)), 0.0; atol=0.1)

@test isapprox(std(getKDEfit(btd,distribution=Normal)), 0.1; atol=0.05)

btd = getBelief(getVariable(fg, :x1))
pts = getPoints(btd)
pts_p = pts[pts .>= 0]
pts_n = pts[pts .< 0]

nfit_p = fit(Normal, pts_p)
@test isapprox(mean(nfit_p), 1.0; atol=0.1)
@test isapprox(std(nfit_p), 0.14; atol=0.05) #TODO confirm the correct value and tolerance

nfit_n = fit(Normal, pts_n)
@test isapprox(mean(nfit_n), -1.0; atol=0.1)
@test isapprox(std(nfit_n), 0.14; atol=0.05) #TODO confirm the correct value and tolerance

# To look at your results
# if false
#   using RoMEPlotting
#   plotKDE(fg, ls(fg))
# end

##

end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sid_opt_rho2_7(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sid_opt_rho2_7!(A_copy)
end

@inline function exp_sid_opt_rho2_7!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=0.045435112444573475
    coeff2=1.0223246692300103
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.045435112444573475
    coeff2=1.0223246692300103
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.013779624734115105
    coeff2=1.0123617774427427
    coeff3=0.02096699267528694
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.0025872950992219526
    coeff2=0.013669635797719902
    coeff3=1.0123048626491173
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.0024179040551112
    coeff2=0.9995075315498182
    coeff3=0.0417319729524033
    coeff4=0.0010158820545039565
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.00023252325808871646
    coeff2=4.612679859250955e-5
    coeff3=0.010778337993824194
    coeff4=0.9999831553043256
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=7.788069195751279e-6
    coeff2=8.668814613800501e-5
    coeff3=0.0005851748274979261
    coeff4=0.0003724118165732193
    coeff5=1.000002894879292
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=4.367862801332661e-5
    coeff2=4.771005078802816e-6
    coeff3=3.057699543578747e-7
    coeff4=1.4369102577037629e-8
    coeff5=2.8715019361588243e-10
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.0014340288130826606
    coeff2=3.0214170368224798
    coeff3=0.28010545658135216
    coeff4=0.017678291021512645
    coeff5=0.0006690728448219196
    coeff6=0.7423612498678709
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.0009787233688326564
    coeff2=-0.096167679660017
    coeff3=-0.010643200832947241
    coeff4=-0.00024824310516588245
    coeff5=0.00012514140616432543
    coeff6=1.228607675289351
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.0003158514159932758
    coeff2=8.290611574892452
    coeff3=1.7421134180734874
    coeff4=0.30262885192060646
    coeff5=0.02034868019711193
    coeff6=249.89631597753709
    coeff7=1.227713124647877
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.0026162370361472948
    coeff2=0.030992248965416677
    coeff3=0.009171582602004843
    coeff4=0.0014654665589331566
    coeff5=8.173851230520504e-5
    coeff6=1.119040499126371
    coeff7=1.4271746834360811e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9996354364827498
    coeff2=1.0075667227584877
    coeff3=0.181554796011574
    coeff4=0.017961596596304457
    coeff5=0.0006512172198766862
    coeff6=6.761523475937711e-6
    coeff7=-0.036931956487199616
    coeff8=1.1186848822520254
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# Params
struct Params <: AbstractDict{Symbol, String}
    dict::Dict{Symbol, String}
end

Params() = Params(Dict{Symbol, String}())
Params(kv::Iterators.Zip) = Params(Dict(kv))
Params(kv::Pair{Symbol, String}...) = Params(Dict(kv))

function createparams(request::HTTP.Request, path::String)
    paths = URIs.splitpath(path)
    p_names = filter(str -> startswith(str, ':'), paths)
    p_idxs = findall(x -> x ∈ p_names, paths)
    k = map(str -> Symbol(str[2:end]), p_names)
    v = URIs.splitpath(request.target)[p_idxs]
    return Params(zip(k, v))
end

# Query
struct Query <: AbstractDict{Symbol, String}
    dict::Dict{Symbol, String}
end

Query() = Query(Dict{Symbol, String}())
Query(kv::Iterators.Zip) = Query(Dict(kv))
Query(kv::Pair{Symbol, String}...) = Query(Dict(kv))

function createquery(request::HTTP.Request)
    dict = queryparams(URI(request.target))
    if isempty(dict)
        return Query()
    else
        k = Symbol.(keys(dict))
        v = values(dict)
        return Query(zip(k, v))
    end
end

const ParamsOrQuery = Union{Params, Query}

# AbstractDict Interface
Base.Dict(pq::ParamsOrQuery) = copy(getfield(pq, :dict))
Base.length(pq::ParamsOrQuery) = length(getfield(pq, :dict))
Base.get(pq::ParamsOrQuery, key::Symbol, default) = get(getfield(pq, :dict), key, default)
Base.get(pq::ParamsOrQuery, key::AbstractString, default) = get(getfield(pq, :dict), Symbol(key), default)
Base.get(f::Base.Callable, pq::ParamsOrQuery, key::Symbol) = get(f, getfield(pq, :dict), key)
Base.get(f::Base.Callable, pq::ParamsOrQuery, key::AbstractString) = get(f, getfield(pq, :dict), Symbol(key))
Base.iterate(pq::ParamsOrQuery, args...) = iterate(getfield(pq, :dict), args...)
# Calling faster haskey and getindex methods defined for Dict type
Base.haskey(pq::ParamsOrQuery, key::Symbol) = haskey(getfield(pq, :dict), key)
Base.haskey(pq::ParamsOrQuery, key::AbstractString) = haskey(getfield(pq, :dict), Symbol(key))
Base.getindex(pq::ParamsOrQuery, idx::Symbol) = getindex(getfield(pq, :dict), idx)
Base.getindex(pq::ParamsOrQuery, idx::AbstractString) = getindex(getfield(pq, :dict), Symbol(idx))
# params.key and query.key syntax
Base.getproperty(pq::ParamsOrQuery, prop::Symbol) = getindex(getfield(pq, :dict), prop)
Base.propertynames(pq::ParamsOrQuery) = collect(keys(pq))

# Req
const BodyTypes = Union{AbstractString, JSON3.Object, JSON3.Array}

struct Req{B<:BodyTypes}
    method::String
    target::String
    params::Params
    query::Query
    body::B
end


text:
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

const fontface = cv_fontface("DejaVu Sans")

function setup_axis(setup::CV_SceneSetupChain)
    label_style = cv_black → fontface → cv_fontsize(20)
    ticks1  = "%.0f" ⇒ -2:2
    ticks1h =     "" ⇒ -2.5:2.5
    ticks2  = "%.0f" ⇒ -3:3
    ticks3 = ("-π" ⇒ -pi, "-π/2" ⇒ -pi/2, "0" ⇒ 0.0, "π/2" ⇒ pi/2, "π" ⇒ pi,)

    app_l = CV_TickLabelAppearance(; label_style, tick_length=10)
    app_s = CV_TickLabelAppearance(; label_style, tick_length=6)
    return cv_setup_lr_axis(setup,
        (app_l ↦ ticks1, app_s ↦ ticks1h,), (app_l ↦ ticks1, app_s ↦ ticks1h,),
        (app_l ↦ ticks2,), (app_l ↦ ticks3,))
end

function setup_star_arc_painters(setup::CV_SceneSetupChain, cut_test=nothing)
    layout = setup.layout
    multiply_pos = CV_MultiplyByFactor(ComplexF64)

    trafo = cv_get_trafo(layout)
    trafo_domain = multiply_pos
    trafo_codomain = w -> trafo(multiply_pos(w))

    common_style = cv_op_over  → cv_linewidth(3) → cv_antialias_best
    star_style = common_style → cv_color(1,1,1,0.8)
    arc_style = common_style → cv_color(0,0,0,0.8)
    star_lines, arc_lines = cv_star_arc_lines(
        tuple(LinRange(0.2, 1.5, 6)...), tuple(LinRange(-π/4, π/4, 6)...))
    CLP = CV_LinePainter
    star_painter_domain   = star_style ↦ CLP(trafo_domain,   star_lines)
    arc_painter_domain    = arc_style  ↦ CLP(trafo_domain,   arc_lines )
    star_painter_codomain = star_style ↦ CLP(trafo_codomain, star_lines, false,
        trafo_domain, cut_test)
    arc_painter_codomain  = arc_style  ↦ CLP(trafo_codomain, arc_lines, false,
        trafo_domain, cut_test)

    state_counter = cv_get_state_counter(layout)
    cc_can_domain = cv_get_cc_can_domain(layout)
    cc_can_codomain = cv_get_cc_can_codomain(layout)

    update_painter_func = z -> begin
        multiply_pos.factor = z
        state = state_counter.value

        if state == 2
            cv_paint(cc_can_domain, arc_painter_domain)
            cv_paint(cc_can_codomain, arc_painter_codomain)
            cv_paint(cc_can_domain, star_painter_domain)
            cv_paint(cc_can_codomain, star_painter_codomain)
        end
        return nothing
    end

    redraw_func = layout -> begin
        update_painter_func(multiply_pos.factor)
        return nothing
    end
    return cv_combine(setup; update_painter_func, redraw_func)
end

function get_slider_rulers()  # {{{
    major_values = Float64.(10:10:50)
    minor_values = Float64.(15:5:45)
    mini_values = setdiff(setdiff(Float64.(10:1:50), minor_values), major_values)

    intro_values = Float64.(1:1:9)

    major_ticks = (cv_format_ticks("%.0f", major_values...)...,
        CV_TickLabel(55.0, "∞"))
    minor_ticks = cv_format_ticks("", minor_values...)
    mini_ticks = cv_format_ticks("", mini_values...)
    intro_ticks = cv_format_ticks("%.0f", intro_values...)

    major_ruler = CV_Ruler(major_ticks,
        CV_TickLabelAppearance(; tick_length=10,
            label_style=cv_black → fontface → cv_fontsize(20)))
    minor_ruler = CV_Ruler(minor_ticks, CV_TickLabelAppearance(; tick_length=7))
    mini_ruler = CV_Ruler(mini_ticks, CV_TickLabelAppearance(; tick_length=4))
    intro_ruler = CV_Ruler(intro_ticks,
        CV_TickLabelAppearance(; tick_length=4,
            label_style=cv_black → fontface → cv_fontsize(10)))
    return (major_ruler, minor_ruler, mini_ruler, intro_ruler)
end # }}}

function create_slider(setup, slider_pos)
    layout = setup.layout
    can_domain_l = cv_get_can_domain_l(layout)
    can_codomain_l = cv_get_can_codomain_l(layout)
    width = can_codomain_l.rectangle.right - can_domain_l.rectangle.left
    mid = cv_half(can_codomain_l.rectangle.right + can_domain_l.rectangle.left)
    bottom = min(can_domain_l.rectangle.bottom, can_codomain_l.rectangle.bottom)
    height = 15

    slider_data = cv_create_hslider(
        width, height, 0.5, 55.5, get_slider_rulers();
        decoraction_with_layout_and_position_callback=(inner_layout, pos) ->
            cv_border(inner_layout, pos, 1))

    cont_slider = slider_data.slider_container
    cc_cont_slider = slider_data.container_context
    cont_slider_l = cv_add_canvas!(layout,
        cont_slider,
        cv_anchor(cont_slider, :slider_south, cont_slider, :south),
        (mid, bottom - Int32(20)))

    set_slider_value = z -> begin
        slider_pos.value = (z > 53) ? 55.0 : max(1.0, min(z, 50.0))
        return CV_Response(;redraw_flag=true)
    end

    bg_painter = cv_color(.8,.8,.8) ↦ CV_FillPainter()
    mark_painter = (cv_op_source → cv_color(0,0,1) → cv_linewidth(2)) ↦
        CV_ValueMarkPainter(slider_pos,
            0.0, imag(cont_slider.can_slider.corner_ul), false)

    setup = cv_setup_hslider(setup, slider_data, cont_slider_l,
        bg_painter → mark_painter, set_slider_value)
    return setup
end

function create_scene(trafo,
        domain, codomain, slider_pos; cut_test=nothing, gap=80,
        axis_label_style=cv_black → fontface → cv_fontsize(20), padding=30) # {{{
    layout = CV_StateLayout(CV_2DLayout(), CV_CyclicValue(2))
    layout = cv_do_lr_layout(cv_add(layout, trafo, domain, codomain), gap)

    setup = cv_setup_cycle_state(CV_LRSetupChain(layout))
    setup = create_slider(setup, slider_pos)
    setup = cv_setup_lr_painters(setup; cut_test,
        parallel_lines_painter_domain=nothing,
        parallel_lines_painter_codomain=nothing)
    setup = setup_star_arc_painters(setup, cut_test)
    setup = setup_axis(setup)
    setup = cv_setup_lr_border(setup)
    padding > 0 && cv_add_padding!(setup.layout, padding)
    setup = cv_setup_domain_codomain_scene(setup)
    cv_scene_lr_start(setup)
    return setup.layout
end # }}}

slider_pos = CV_TranslateByOffset(Float64)
slider_pos.value = 1.2

cut_func = cv_create_angle_cross_test(π, 0.0, Inf; δ=1e-2)

cut_test = (z, w) -> (slider_pos.value == 1.0) ? false : cut_func(z, w)

trafo = z -> (slider_pos.value == 55.0) ? log(z) :
    slider_pos.value * expm1(log(z)/slider_pos.value)
domain   = CV_Math2DCanvas(-2.5 + 2.5im, 2.5 - 2.5im, 100)
codomain = CV_Math2DCanvas(-3.5 + (pi+0.1)im, 3.5 - (pi+0.1)im, 100)

scene = create_scene(trafo, domain, codomain, slider_pos; cut_test)

cv_save_image(cv_get_can_layout(scene), "./Log01.png")

if CVG !== nothing
    handler = CVG.cvg_visualize(scene)
    CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
@testset "Momentum Gradient Descent" begin
    run_optim_tests(MomentumGradientDescent(),
                    skip = ("Powell", "Rosenbrock"),
                    convergence_exceptions = (("Large Polynomial",1),
                     ("Himmelblau",1), ("Powell", 1)),
                    minimizer_exceptions = (("Powell", 2),))
end


text:
# This file is a part of Julia. Except for the *_kernel functions (see below),
# license is MIT: https://julialang.org/license

struct DoubleFloat64
    hi::Float64
    lo::Float64
end
struct DoubleFloat32
    hi::Float64
end

# sin_kernel and cos_kernel functions are only valid for |x| < pi/4 = 0.7854
# translated from openlibm code: k_sin.c, k_cos.c, k_sinf.c, k_cosf.c.
# asin functions are based on openlibm code: e_asin.c, e_asinf.c. The above
# functions are made available under the following licence:

## Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
##
## Developed at SunPro, a Sun Microsystems, Inc. business.
## Permission to use, copy, modify, and distribute this
## software is freely granted, provided that this notice
## is preserved.

function sin_kernel(x::DoubleFloat64)
    S1 = -1.66666666666666324348e-01
    S2 =  8.33333333332248946124e-03
    S3 = -1.98412698298579493134e-04
    S4 =  2.75573137070700676789e-06
    S5 = -2.50507602534068634195e-08
    S6 =  1.58969099521155010221e-10

    z = x.hi*x.hi
    w = z*z
    r = S2+z*(S3+z*S4) + z*w*(S5+z*S6)
    v = z*x.hi
    x.hi-((z*(0.5*x.lo-v*r)-x.lo)-v*S1)
end

function cos_kernel(x::DoubleFloat64)
    C1 =  4.16666666666666019037e-02
    C2 = -1.38888888888741095749e-03
    C3 =  2.48015872894767294178e-05
    C4 = -2.75573143513906633035e-07
    C5 =  2.08757232129817482790e-09
    C6 = -1.13596475577881948265e-11

    z = x.hi*x.hi
    w = z*z
    r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6))
    hz = 0.5*z
    w = 1.0-hz
    w + (((1.0-w)-hz) + (z*r-x.hi*x.lo))
end

function sin_kernel(x::DoubleFloat32)
    S1 = -0.16666666641626524
    S2 = 0.008333329385889463
    S3 = -0.00019839334836096632
    S4 = 2.718311493989822e-6

    z = x.hi*x.hi
    w = z*z
    r = S3+z*S4
    s = z*x.hi
    Float32((x.hi + s*(S1+z*S2)) + s*w*r)
end

function cos_kernel(x::DoubleFloat32)
    C0 = -0.499999997251031
    C1 = 0.04166662332373906
    C2 = -0.001388676377460993
    C3 = 2.439044879627741e-5

    z = x.hi*x.hi
    w = z*z
    r = C2+z*C3
    Float32(((1.0+z*C0) + w*C1) + (w*z)*r)
end

# fallback methods
sin_kernel(x::Real) = sin(x)
cos_kernel(x::Real) = cos(x)

# Inverse trigonometric functions

# asin methods
ASIN_X_MIN_THRESHOLD(::Type{Float32}) = 2.0f0^-12
ASIN_X_MIN_THRESHOLD(::Type{Float64}) = sqrt(eps(Float64))

arc_p(t::Float64) =
    t*@horner(t,
    1.66666666666666657415e-01,
    -3.25565818622400915405e-01,
    2.01212532134862925881e-01,
    -4.00555345006794114027e-02,
    7.91534994289814532176e-04,
    3.47933107596021167570e-05)

arc_q(t::Float64) =
    @horner(t,
    1.0,
    -2.40339491173441421878e+00,
    2.02094576023350569471e+00,
    -6.88283971605453293030e-01,
    7.70381505559019352791e-02)

arc_p(t::Float32) =
    t*@horner(t,
    1.6666586697f-01,
    -4.2743422091f-02,
    -8.6563630030f-03)

arc_q(t::Float32) = @horner(t, 1.0f0, -7.0662963390f-01)

@inline function asin_kernel(t::Float64, x::Float64)
    pio2_lo = 6.12323399573676603587e-17
    s = sqrt_llvm(t)
    p = arc_p(t) # numerator polynomial
    q = arc_q(t) # denominator polynomial
    if abs(x) >= 0.975 # |x| > 0.975
        Rx = p/q
        return flipsign(pi/2 - (2.0*(s + s*Rx) - pio2_lo), x)
    else
        s0 = reinterpret(Float64, (reinterpret(UInt64, s) >> 32) << 32)
        c = (t - s0*s0)/(s + s0)
        Rx = p/q
        p = 2.0*s*Rx - (pio2_lo - 2.0*c)
        q = pi/4 - 2.0*s0
        return flipsign(pi/4 - (p-q), x)
    end
end
@inline function asin_kernel(t::Float32, x::Float32)
    s = sqrt_llvm(Float64(t))
    p = arc_p(t) # numerator polynomial
    q = arc_q(t) # denominator polynomial
    Rx = p/q # rational approximation
    flipsign(Float32(pi/2 - 2*(s + s*Rx)), x)
end

@noinline asin_domain_error(x) = throw(DomainError(x, "asin(x) is not defined for |x|>1."))
function asin(x::T) where T<:Union{Float32, Float64}
    # Method :
    #    Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
    #    we approximate asin(x) on [0,0.5] by
    #        asin(x) = x + x*x^2*R(x^2)
    #    where
    #        R(x^2) is a rational approximation of (asin(x)-x)/x^3
    #    and its remez error is bounded by
    #        |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
    #
    #    For x in [0.5,1]
    #        asin(x) = pi/2-2*asin(sqrt((1-x)/2))
    #    Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
    #    then for x>0.98
    #        asin(x) = pi/2 - 2*(s+s*z*R(z))
    #            = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
    #    For x<=0.98, let pio4_hi = pio2_hi/2, then
    #        f = hi part of s;
    #        c = sqrt(z) - f = (z-f*f)/(s+f)     ...f+c=sqrt(z)
    #    and
    #        asin(x) = pi/2 - 2*(s+s*z*R(z))
    #            = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
    #            = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
    #
    # Special cases:
    #    if |x|>1, throw DomainError
    absx = abs(x)
    if absx >= T(1.0) # |x|>= 1
        if absx == T(1.0)
            return flipsign(T(pi)/2, x)
        end
        asin_domain_error(x)
    elseif absx < T(1.0)/2
        # if |x| sufficiently small, |x| is a good approximation
        if absx < ASIN_X_MIN_THRESHOLD(T)
            return x
        end
        # else if |x|<0.5 we use a rational approximation R(x)=p(x)/q(x) such that
        # tan(x) ≈ x+x*R(x)
        x² = x*x
        Rx = arc_p(x²)/arc_q(x²) # rational approximation
        return muladd(x, Rx, x)
    end
    # else 1/2 <= |x| < 1
    t = (T(1.0) - absx)/2
    return asin_kernel(t, x)
end

# multiply in extended precision
function mulpi_ext(x::Float64)
    m = 3.141592653589793
    m_hi = 3.1415926218032837
    m_lo = 3.178650954705639e-8

    x_hi = reinterpret(Float64, reinterpret(UInt64,x) & 0xffff_ffff_f800_0000)
    x_lo = x-x_hi

    y_hi = m*x
    y_lo = x_hi * m_lo + (x_lo* m_hi + ((x_hi*m_hi-y_hi) + x_lo*m_lo))

    DoubleFloat64(y_hi,y_lo)
end
mulpi_ext(x::Float32) = DoubleFloat32(pi*Float64(x))
mulpi_ext(x::Rational) = mulpi_ext(float(x))
mulpi_ext(x::Real) = pi*x # Fallback

"""
    sinpi(x)

Compute ``\\sin(\\pi x)`` more accurately than `sin(pi*x)`, especially for large `x`.
"""
function sinpi(x::T) where T<:AbstractFloat
    if !isfinite(x)
        isnan(x) && return x
        throw(DomainError(x, "`x` cannot be infinite."))
    end

    ax = abs(x)
    s = maxintfloat(T)/2
    ax >= s && return copysign(zero(T),x) # integer-valued

    # reduce to interval [-1,1]
    # assumes RoundNearest rounding mode
    t = 3*s
    rx = x-((x+t)-t) # zeros may be incorrectly signed
    arx = abs(rx)

    if (arx == 0) | (arx == 1)
        copysign(zero(T),x)
    elseif arx < 0.25
        sin_kernel(mulpi_ext(rx))
    elseif arx < 0.75
        y = mulpi_ext(T(0.5) - arx)
        copysign(cos_kernel(y),rx)
    else
        y = mulpi_ext(copysign(one(T),rx) - rx)
        sin_kernel(y)
    end
end

# Integers and Rationals
function sinpi(x::T) where T<:Union{Integer,Rational}
    Tf = float(T)
    if !isfinite(x)
        throw(DomainError(x, "`x` must be finite."))
    end

    # until we get an IEEE remainder function (#9283)
    rx = rem(x,2)
    if rx > 1
        rx -= 2
    elseif rx < -1
        rx += 2
    end
    arx = abs(rx)

    if (arx == 0) | (arx == 1)
        copysign(zero(Tf),x)
    elseif arx < 0.25
        sin_kernel(mulpi_ext(rx))
    elseif arx < 0.75
        y = mulpi_ext(T(0.5) - arx)
        copysign(cos_kernel(y),rx)
    else
        y = mulpi_ext(copysign(one(T),rx) - rx)
        sin_kernel(y)
    end
end

"""
    cospi(x)

Compute ``\\cos(\\pi x)`` more accurately than `cos(pi*x)`, especially for large `x`.
"""
function cospi(x::T) where T<:AbstractFloat
    if !isfinite(x)
        isnan(x) && return x
        throw(DomainError(x, "`x` cannot be infinite."))
    end

    ax = abs(x)
    s = maxintfloat(T)
    ax >= s && return one(T) # even integer-valued

    # reduce to interval [-1,1], then [0,1]
    # assumes RoundNearest rounding mode
    rx = abs(ax-((ax+s)-s))

    if rx <= 0.25
        cos_kernel(mulpi_ext(rx))
    elseif rx < 0.75
        y = mulpi_ext(T(0.5) - rx)
        sin_kernel(y)
    else
        y = mulpi_ext(one(T) - rx)
        -cos_kernel(y)
    end
end

# Integers and Rationals
function cospi(x::T) where T<:Union{Integer,Rational}
    if !isfinite(x)
        throw(DomainError(x, "`x` must be finite."))
    end

    ax = abs(x)
    # until we get an IEEE remainder function (#9283)
    rx = rem(ax,2)
    if rx > 1
        rx = 2-rx
    end

    if rx <= 0.25
        cos_kernel(mulpi_ext(rx))
    elseif rx < 0.75
        y = mulpi_ext(T(0.5) - rx)
        sin_kernel(y)
    else
        y = mulpi_ext(one(T) - rx)
        -cos_kernel(y)
    end
end

sinpi(x::Integer) = x >= 0 ? zero(float(x)) : -zero(float(x))
cospi(x::Integer) = isodd(x) ? -one(float(x)) : one(float(x))
sinpi(x::Real) = sinpi(float(x))
cospi(x::Real) = cospi(float(x))

function sinpi(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if isinteger(zr)
        # zr = ...,-2,-1,0,1,2,...
        # sin(pi*zr) == ±0
        # cos(pi*zr) == ±1
        # cosh(pi*zi) > 0
        s = copysign(zero(F),zr)
        c_pos = isa(zr,Integer) ? iseven(zr) : isinteger(zr/2)
        sh = sinh(pi*zi)
        Complex(s, c_pos ? sh : -sh)
    elseif isinteger(2*zr)
        # zr = ...,-1.5,-0.5,0.5,1.5,2.5,...
        # sin(pi*zr) == ±1
        # cos(pi*zr) == +0
        # sign(sinh(pi*zi)) == sign(zi)
        s_pos = isinteger((2*zr-1)/4)
        ch = cosh(pi*zi)
        Complex(s_pos ? ch : -ch, isnan(zi) ? zero(F) : copysign(zero(F),zi))
    elseif !isfinite(zr)
        if zi == 0 || isinf(zi)
            Complex(F(NaN), F(zi))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        pizi = pi*zi
        Complex(sinpi(zr)*cosh(pizi), cospi(zr)*sinh(pizi))
    end
end

function cospi(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if isinteger(zr)
        # zr = ...,-2,-1,0,1,2,...
        # sin(pi*zr) == ±0
        # cos(pi*zr) == ±1
        # sign(sinh(pi*zi)) == sign(zi)
        # cosh(pi*zi) > 0
        s = copysign(zero(F),zr)
        c_pos = isa(zr,Integer) ? iseven(zr) : isinteger(zr/2)
        ch = cosh(pi*zi)
        Complex(c_pos ? ch : -ch, isnan(zi) ? s : -flipsign(s,zi))
    elseif isinteger(2*zr)
        # zr = ...,-1.5,-0.5,0.5,1.5,2.5,...
        # sin(pi*zr) == ±1
        # cos(pi*zr) == +0
        # sign(sinh(pi*zi)) == sign(zi)
        s_pos = isinteger((2*zr-1)/4)
        sh = sinh(pi*zi)
        Complex(zero(F), s_pos ? -sh : sh)
    elseif !isfinite(zr)
        if zi == 0
            Complex(F(NaN), isnan(zr) ? zero(F) : -flipsign(F(zi),zr))
        elseif isinf(zi)
            Complex(F(Inf), F(NaN))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        pizi = pi*zi
        Complex(cospi(zr)*cosh(pizi), -sinpi(zr)*sinh(pizi))
    end
end

"""
    sinc(x)

Compute ``\\sin(\\pi x) / (\\pi x)`` if ``x \\neq 0``, and ``1`` if ``x = 0``.
"""
sinc(x::Number) = x==0 ? one(x)  : oftype(x,sinpi(x)/(pi*x))
sinc(x::Integer) = x==0 ? one(x) : zero(x)
sinc(x::Complex{<:AbstractFloat}) = x==0 ? one(x) : oftype(x, sinpi(x)/(pi*x))
sinc(x::Complex) = sinc(float(x))
sinc(x::Real) = x==0 ? one(x) : isinf(x) ? zero(x) : sinpi(x)/(pi*x)

"""
    cosc(x)

Compute ``\\cos(\\pi x) / x - \\sin(\\pi x) / (\\pi x^2)`` if ``x \\neq 0``, and ``0`` if
``x = 0``. This is the derivative of `sinc(x)`.
"""
cosc(x::Number) = x==0 ? zero(x) : oftype(x,(cospi(x)-sinpi(x)/(pi*x))/x)
cosc(x::Integer) = cosc(float(x))
cosc(x::Complex{<:AbstractFloat}) = x==0 ? zero(x) : oftype(x,(cospi(x)-sinpi(x)/(pi*x))/x)
cosc(x::Complex) = cosc(float(x))
cosc(x::Real) = x==0 || isinf(x) ? zero(x) : (cospi(x)-sinpi(x)/(pi*x))/x

for (finv, f, finvh, fh, finvd, fd, fn) in ((:sec, :cos, :sech, :cosh, :secd, :cosd, "secant"),
                                            (:csc, :sin, :csch, :sinh, :cscd, :sind, "cosecant"),
                                            (:cot, :tan, :coth, :tanh, :cotd, :tand, "cotangent"))
    name = string(finv)
    hname = string(finvh)
    dname = string(finvd)
    @eval begin
        @doc """
            $($name)(x)

        Compute the $($fn) of `x`, where `x` is in radians.
        """ ($finv)(z::T) where {T<:Number} = one(T) / (($f)(z))
        @doc """
            $($hname)(x)

        Compute the hyperbolic $($fn) of `x`.
        """ ($finvh)(z::T) where {T<:Number} = one(T) / (($fh)(z))
        @doc """
            $($dname)(x)

        Compute the $($fn) of `x`, where `x` is in degrees.
        """ ($finvd)(z::T) where {T<:Number} = one(T) / (($fd)(z))
    end
end

for (tfa, tfainv, hfa, hfainv, fn) in ((:asec, :acos, :asech, :acosh, "secant"),
                                       (:acsc, :asin, :acsch, :asinh, "cosecant"),
                                       (:acot, :atan, :acoth, :atanh, "cotangent"))
    tname = string(tfa)
    hname = string(hfa)
    @eval begin
        @doc """
            $($tname)(x)
        Compute the inverse $($fn) of `x`, where the output is in radians. """ ($tfa)(y::T) where {T<:Number} = ($tfainv)(one(T) / y)
        @doc """
            $($hname)(x)
        Compute the inverse hyperbolic $($fn) of `x`. """ ($hfa)(y::T) where {T<:Number} = ($hfainv)(one(T) / y)
    end
end

# multiply in extended precision
function deg2rad_ext(x::Float64)
    m = 0.017453292519943295
    m_hi = 0.01745329238474369
    m_lo = 1.3519960527851425e-10

    u = 134217729.0*x # 0x1p27 + 1
    x_hi = u-(u-x)
    x_lo = x-x_hi

    y_hi = m*x
    y_lo = x_hi * m_lo + (x_lo* m_hi + ((x_hi*m_hi-y_hi) + x_lo*m_lo))

    DoubleFloat64(y_hi,y_lo)
end
deg2rad_ext(x::Float32) = DoubleFloat32(deg2rad(Float64(x)))
deg2rad_ext(x::Real) = deg2rad(x) # Fallback

function sind(x::Real)
    if isinf(x)
        return throw(DomainError(x, "`x` cannot be infinite."))
    elseif isnan(x)
        return oftype(x,NaN)
    end

    rx = copysign(float(rem(x,360)),x)
    arx = abs(rx)

    if rx == zero(rx)
        return rx
    elseif arx < oftype(rx,45)
        return sin_kernel(deg2rad_ext(rx))
    elseif arx <= oftype(rx,135)
        y = deg2rad_ext(oftype(rx,90) - arx)
        return copysign(cos_kernel(y),rx)
    elseif arx == oftype(rx,180)
        return copysign(zero(rx),rx)
    elseif arx < oftype(rx,225)
        y = deg2rad_ext((oftype(rx,180) - arx)*sign(rx))
        return sin_kernel(y)
    elseif arx <= oftype(rx,315)
        y = deg2rad_ext(oftype(rx,270) - arx)
        return -copysign(cos_kernel(y),rx)
    else
        y = deg2rad_ext(rx - copysign(oftype(rx,360),rx))
        return sin_kernel(y)
    end
end

function cosd(x::Real)
    if isinf(x)
        return throw(DomainError(x, "`x` cannot be infinite."))
    elseif isnan(x)
        return oftype(x,NaN)
    end

    rx = abs(float(rem(x,360)))

    if rx <= oftype(rx,45)
        return cos_kernel(deg2rad_ext(rx))
    elseif rx < oftype(rx,135)
        y = deg2rad_ext(oftype(rx,90) - rx)
        return sin_kernel(y)
    elseif rx <= oftype(rx,225)
        y = deg2rad_ext(oftype(rx,180) - rx)
        return -cos_kernel(y)
    elseif rx < oftype(rx,315)
        y = deg2rad_ext(rx - oftype(rx,270))
        return sin_kernel(y)
    else
        y = deg2rad_ext(oftype(rx,360) - rx)
        return cos_kernel(y)
    end
end

tand(x::Real) = sind(x) / cosd(x)

for (fd, f, fn) in ((:sind, :sin, "sine"), (:cosd, :cos, "cosine"), (:tand, :tan, "tangent"))
    name = string(fd)
    @eval begin
        @doc """
            $($name)(x)
        Compute $($fn) of `x`, where `x` is in degrees. """ ($fd)(z) = ($f)(deg2rad(z))
    end
end

for (fd, f, fn) in ((:asind, :asin, "sine"), (:acosd, :acos, "cosine"), (:atand, :atan, "tangent"),
                    (:asecd, :asec, "secant"), (:acscd, :acsc, "cosecant"), (:acotd, :acot, "cotangent"))
    name = string(fd)
    @eval begin
        @doc """
            $($name)(x)

        Compute the inverse $($fn) of `x`, where the output is in degrees. """ ($fd)(y) = rad2deg(($f)(y))
    end
end

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
using Printf

let A = Point(0.0, 0.0),
    B = Point(0.0, 10.0),
    C = Point(10.0, 10.0),
    D = Point(10.0, 0.0),
    E = Point(2.5, 2.5),
    F = Point(2.5, 7.5),
    G = Point(7.5, 7.5),
    H = Point(7.5, 2.5),
    I = Point(3.0, 0.0),
    J = Point(7.0, 0.0),
    K = Point(10.0, 5.0),
    L = Point(7.0, 10.0),
    M = Point(3.0, 10.0),
    N = Point(0.0, 5.0),
    testpts = (Point(5.0, 5.0), Point(5.0, 8.0), Point(-10.0, 5.0), Point(0.0, 5.0),
        Point(10.0, 5.0), Point(8.0, 5.0), Point(10.0, 10.0))

    square = RayCastings.connect(A, B, C, D)
    square_withhole = vcat(square, RayCastings.connect(E, F, G, H))
    strange = RayCastings.connect(A, E, B, F, G, C, D, E)
    exagon = RayCastings.connect(I, J, K, L, M, N)

    println("\n# TESTING WHETHER POINTS ARE WITHIN POLYGONS")
    for poly in (square, square_withhole, strange, exagon)
        println("\nEdges: \n - ", join(poly, "\n - "))
        println("Inside/outside:")
        for p in testpts
            @printf(" - %-12s is %s\n", p, RayCastings.isinside(poly, p) ? "inside" : "outside")
        end
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
####################################
# barplot of cell counts per sample
####################################
plotCounts(fcsRaw, md)

####################################################################
# PCA plot
####################################################################

plotPCA(daf, md)
=#


text:
invoke_map!{F, N}(f::F, dest, As::Vararg{DataValueArray, N}) =
    invoke(map!, Tuple{F, AbstractArray, Vararg{AbstractArray, N}}, f, dest, As...)

"""
    map(f, As::DataValueArray...)

Call `map` with DataValue lifting semantics and return a `DataValueArray`.
Lifting means calling function `f` on the the values wrapped inside `DataValue` entries
of the input arrays, and returning null if any entry is missing.

Note that this method's signature specifies the source `As` arrays as all
`DataValueArray`s. Thus, calling `map` on arguments consisting
of both `Array`s and `DataValueArray`s will fall back to the standard implementation
of `map` (i.e. without lifting).
"""
function Base.map{F}(f::F, As::DataValueArray...)
    T = nullable_broadcast_eltype(f, As...)
    dest = similar(DataValueArray{T}, size(As[1]))
    invoke_map!(f, dest, As...)
end

"""
    map!(f, dest::DataValueArray, As::DataValueArray...)

Call `map!` with DataValue lifting semantics.
Lifting means calling function `f` on the the values wrapped inside `DataValue` entries
of the input arrays, and returning null if any entry is missing.

Note that this method's signature specifies the destination `dest` array as well as the
source `As` arrays as all `DataValueArray`s. Thus, calling `map!` on a arguments
consisting of both `Array`s and `DataValueArray`s will fall back to the standard implementation
of `map!` (i.e. without lifting).
"""
function Base.map!{F}(f::F, dest::DataValueArray, As::DataValueArray...)
    invoke_map!(f, dest, As...)
end

# This definition is needed to avoid dispatch loops going back to the above one
function Base.map!{F}(f::F, dest::DataValueArray)
    invoke_map!(f, dest, dest)
end


text:
# SPDXLicense-Identifier: MIT

# Stable implementation of the log(exp(x) + exp(y)).
function logaddexp(x::T, y::T) where T
    diff = zero(T)
    if x < y
        diff = x - y
        x = y
    else y < x
        diff = y - x
    end

    if diff >= log(eps(T))
        return x + log1p(exp(diff))
    end
    return x
end

"""
    struct LogSemifield{T<:AbstractFloat} <: Semfield
        val::T
    end

Log-semifield is defined as :
  * ``x \\oplus y \\triangleq \\ln( e^x + e^y)``
  * ``x \\otimes y \\triangleq x + y``
  * ``x \\oslash y \\triangleq x - y``
``\\forall x, y \\in \\mathbb{R}``.
"""
struct LogSemifield{T<:AbstractFloat} <: Semifield
    val::T
end

Base.:+(x::LogSemifield{T}, y::LogSemifield{T}) where T =
    LogSemifield{T}(logaddexp(x.val, y.val))
Base.:*(x::LogSemifield, y::LogSemifield) = LogSemifield(x.val + y.val)
Base.:/(x::LogSemifield, y::LogSemifield) = LogSemifield(x.val - y.val)
Base.zero(::Type{LogSemifield{T}}) where T = LogSemifield{T}(T(-Inf))
Base.one(::Type{LogSemifield{T}}) where T = LogSemifield{T}(T(0))
Base.isless(x::LogSemifield, y::LogSemifield) = isless(x.val, y.val)
Base.typemin(x::Type{LogSemifield{T}}) where T = LogSemifield{T}(typemin(T))
Base.typemax(x::Type{LogSemifield{T}}) where T = LogSemifield{T}(typemax(T))
Base.conj(x::LogSemifield{T}) where T = LogSemifield{T}(conj(x.val))

max_stars_repo_name: MartinKocour/HiddenMarkovModel

max_issues_repo_name: MartinKocour/HiddenMarkovModel

max_forks_repo_name: MartinKocour/HiddenMarkovModel