text:
function var_info_test(T)
    rng = MersenneTwister(123);
    @testset "VarInfo $T" begin
        vi = make_test_var_info(:x, T, nothing, rng);
        x = make_test_values(vi, rng);
        @test check_var(vi, x; silent = false)
        lb, ub = bounds(vi);
        @test lb == get_option(vi, :lb);
        @test ub == get_option(vi, :ub);

        vi2 = deepcopy(vi);
        @test isequal(vi, vi2)

        vType = eltype(x);
        if isscalar(vi)
            x2 = ub + one(vType);
        else
            x2 = copy(x);
            x2[1] = ub + one(vType);
        end
        @test !check_var(vi, x2; silent = true)

        # Test grandSum +++++

        # Test dimSums +++++

        # x2 = copy(x);
        # x2[2,1] = ubs[2,1] + T(1e-5);
        # @test !var_check(vi, x2; silent = true)

        # @test !var_check(vi, zeros(Float16, sz...))
    end
end

# function frac_info_test(T)
#     @testset "Fractions" begin
#         rng = MersenneTwister(12);
#         sz = (4,3);
#         vi = fraction_info(:x, "x", T, sz);
#         x = rand(rng, T, sz...);
#         @test var_check(vi, x; silent = false);
#         x[3] = T(1.0 + 1e-5);
#         @test !var_check(vi, x)
#     end
# end

# function share_info_test(T)
#     @testset "Shares" begin
#         rng = MersenneTwister(12);
#         sz = (4,3);
#         vi = share_info(:x, "x", T, sz);
#         x = rand(rng, T, sz...);
#         @test !var_check(vi, x);
#         @test var_check(vi, x ./ sum(x))
#     end
# end

@testset "All" begin
    for T in (Float64, Vector{Float32}, Matrix{Int}, UInt8)
        var_info_test(T);
        # frac_info_test(T);
        # share_info_test(T);
    end
    # gs_basics_test();
end

# ----------


text:
using MatchCore

# example assuming * operation is always binary

# ENV["JULIA_DEBUG"] = Metatheory

struct NumberFold <: AbstractAnalysis
    egraph::EGraph
    data::Dict{Int64, Any}
end

# Mandatory for AbstractAnalysis
NumberFold(g::EGraph) = NumberFold(g, Dict{Int64, Any}())

# This should be auto-generated by a macro
function EGraphs.make(an::NumberFold, n)
    n isa Number && return n

    if n isa Expr && Meta.isexpr(n, :call)
        if n.args[1] == :*
            l = n.args[2]
            r = n.args[3]

            if get(an, l, nothing) isa Number && get(an, r, nothing) isa Number
                return an[l] * an[r]
            end
        elseif n.args[1] == :+
            l = n.args[2]
            r = n.args[3]

            if get(an, l, nothing) isa Number && get(an, r, nothing) isa Number
                return an[l] + an[r]
            end
        end
    end
    return nothing
end

function EGraphs.join(analysis::NumberFold, from, to)
    if from isa Number
        if to isa Number
            @assert from == to
        else return from
        end
    end
    return to
end

function EGraphs.modify!(an::NumberFold, id::Int64)
    g = an.egraph
    # !haskey(an, id) && return nothing
    if an[id] isa Number
        newclass = EGraphs.add!(g, an[id])
        merge!(g, newclass.id, id)
    end
end

Base.setindex!(an::NumberFold, value, id::Int64) = setindex!(an.data, value, id)
Base.getindex(an::NumberFold, id::Int64) = an.data[id]
Base.haskey(an::NumberFold, id::Int64) = haskey(an.data, id)
Base.delete!(an::NumberFold, id::Int64) = delete!(an.data, id)

comm_monoid = @theory begin
    a * b => b * a
    a * 1 => a
    a * (b * c) => (a * b) * c
end

G = EGraph(:(3 * 4))
addanalysis!(G, NumberFold)
@testset "Basic Constant Folding Example - Commutative Monoid" begin
    # addanalysis!(G, NumberFold())
    @test (true == @areequalg G comm_monoid 3 * 4 12)

    @test (true == @areequalg G comm_monoid 3 * 4 12 4*3  6*2)
end

@testset "Basic Constant Folding Example 2 - Commutative Monoid" begin
    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    addanalysis!(G, NumberFold)
    display(G.M); println()
    display(G.analyses[1].data); println()
    @test (true == @areequalg G comm_monoid (3 * a) * (4 * b) (12*a)*b ((6*2)*b)*a)
end

@testset "Basic Constant Folding Example - Adding analysis after saturation" begin
    G = EGraph(:(3 * 4))
    # addexpr!(G, 12)
    saturate!(G, comm_monoid)
    addexpr!(G, :(a * 2))
    addanalysis!(G, NumberFold)
    saturate!(G, comm_monoid)

    display(G.M); println()
    println(G.root)
    display(G.analyses[1].data); println()

    @test (true == areequal(G, comm_monoid, :(3 * 4), 12, :(4*3), :(6*2)))

    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    addanalysis!(G, NumberFold)
    @test (true == areequal(G, comm_monoid, :((3 * a) * (4 * b)), :((12*a)*b),
        :(((6*2)*b)*a)))
end

@testset "Infinite Loops analysis" begin
    boson = @theory begin
        1 * x => x
    end

    G = EGraph(Util.cleanast( :(1 * x) ))
    saturate!(G,boson, timeout=100)
    extractor = addanalysis!(G, ExtractionAnalysis, astsize)
    ex = extract!(G, extractor)

    display(G.M); println()
    display(G.analyses[1].data); println()

    println(ex)

    using Metatheory.EGraphs
    boson = @theory begin
        (:c * :cdag) => :cdag * :c + 1
        a * (b + c) => (a * b) + (a * c)
        (b + c) * a => (b * a) + (c * a)
        # 1 * x => x
        (a * b) * c => a * (b * c)
        a * (b * c) => (a * b) * c
    end

    G = EGraph(Util.cleanast( :(c * c * cdag * cdag) ))
    saturate!(G,boson)
    extractor = addanalysis!(G, ExtractionAnalysis, astsize_inv)
    ex = extract!(G, extractor)

    println(ex)
end

ext: jl


text:
sqeuclidean(a,b) = sqrt(sum((a.-b).^2))
quantization_objective(sumBTB::AbstractMatrix, sumBTx::AbstractArray, C::AbstractArray) = dot((sumBTB*C),C)-(2*dot(sumBTx,C))

"""
    function decompose_error(x_i, x_approx_i)

Computes the parallel and orthogonal error of the approximated codeword `̃x_i` w.r.t. the original data point `x_i`.
h_par = x((x - ̃x)⋅x) / (x/x⋅x).\n
h_orthog = x-̃x - h_par
"""
function decompose_error(x_i, x_approx_i)
    parallel = (dot((x_i - x_approx_i), x_i)*x_i)/dot(x_i,x_i)
    orthog = (x_i - x_approx_i) - parallel
    return parallel, orthog
end

"""
    function anisotropic_loss(x_i::AbstractArray, x_approx_i::AbstractArray, η::AnisotropicWeights)

Computes and returns the anisotropic loss between the approximated codeword `̃x_i` and
the original data point `x_i`, based on the input weights for the parallel loss and orthogonal loss `η`.
"""
function anisotropic_loss(x_i::AbstractArray, x_approx_i::AbstractArray, η::AnisotropicWeights)
    parallel, orthog = decompose_error(x_i, x_approx_i)
    return η.h_par * dot(parallel,parallel) + η.h_orthog * dot(orthog,orthog)
end


text:
export CDS_scaled_add!

"""
 Computes A = A + alpha * B for A and B in the compressed diagonal storage format (CDS/DIA)
 TO DO: make this function multi-threaded per column of A and B, or use BLAS functions
"""
function CDS_scaled_add!(
                  A           ::Array{TF,2},
                  B           ::Array{TF,2},
                  A_offsets   ::Vector{TI},
                  B_offsets   ::Vector{TI},
                  alpha       ::TF
                  ) where {TF<:Real,TI<:Integer}

for k=1:length(B_offsets)
  A_update_col = findall((in)(B_offsets[k]),A_offsets)
  if isempty(A_update_col) == true
    error("attempted to update a diagonal in A in CDS storage that does not exist. A and B need to have the same nonzero diagonals")
  end
  A[:,A_update_col] .= A[:,A_update_col] .+ alpha .* B[:,k];
end

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# indicator of the zero cone

export IndZero

"""
**Indicator of the zero cone**

    IndZero()

Returns the indicator function of the set containing the origin, the "zero cone".
"""
struct IndZero <: ProximableFunction end

is_separable(f::IndZero) = true
is_convex(f::IndZero) = true
is_singleton(f::IndZero) = true
is_cone(f::IndZero) = true
is_affine(f::IndZero) = true

function (f::IndZero)(x::AbstractArray{T}) where T <: RealOrComplex
  for k in eachindex(x)
    if x[k] != zero(T)
      return Inf
    end
  end
  return 0.0
end

function prox!(y::AbstractArray{T}, f::IndZero, x::AbstractArray{T}, gamma::Real=1.0) where T <: RealOrComplex
  for k in eachindex(x)
    y[k] = zero(T)
  end
  return 0.0
end

prox!(y::AbstractArray{T}, f::IndZero, x::AbstractArray{T}, gamma::AbstractArray) where {T <: RealOrComplex} = prox!(y, f, x, 1.0)

fun_name(f::IndZero) = "indicator of the zero cone"
fun_dom(f::IndZero) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndZero) = "x ↦ 0 if all(x = 0), +∞ otherwise"
fun_params(f::IndZero) = "none"

function prox_naive(f::IndZero, x::AbstractArray, gamma::Real=1.0)
  return zero(x), 0.0
end

prox_naive(f::IndZero, x::AbstractArray, gamma::AbstractArray) = prox_naive(f, x, 1.0)


text:
"""
Dump all the resort scripts here.  Currently supports:
    - EOS-SG
    - NGL-UNAVCO
"""

function gnssfol(groot::AbstractString,stationname::AbstractString)
    gfol = joinpath(groot,stationname)
    if !isdir(gfol)
        @info "$(Dates.now()) - GNSS Zenith Wet Delay data directory for $(stationname) does not exist."
        @info "$(Dates.now()) - Creating data directory $(gfol)."; mkpath(gfol);
    end
    return gfol
end

function gnssresort(dataID::AbstractString,region,path::AbstractString)

    groot = gnssroot(path,dataID);
    if dataID == "EOS-SG"
        stations = eosload(region); eosresort(stations,groot);
    #elseif dataID == "NGL-UNAVCO"
    #    stations = nglload(region); nglresort(stations,groot);
    else
        error("$(Dates.now()) - $(dataID) is currently not a data resource supported by ClimateGNSS.jl.  Should you be interested, please consider approaching Nathanael Wong (nathanaelwong@fas.harvard.edu) to contribute to ClimateGNSS.jl and create an additional functionality for $(dataID).")
    end

end


text:
"""
    simulation(r, v, m, sld::Dict{Any, Any}, g::Number = NaN, t_step::Number = NaN, t_max::Number = NaN)

Runs the Wisdom-Holfman map of the given N-body problem that defines the following Hamiltonian system.

```math
\\mathcal{H}(\\textbf{q}, \\textbf{p}) = \\frac{\\textbf{p'}^{2}_{1}}{2 m'_{1}} + \\sum_{i=2}^{N} \\frac{\\textbf{p'}^{2}_{i}}{2 m_{i}'} - \\sum_{i=2}^{N} \\frac{Gm_{i}^{'}M_{i}}{\\mid \\textbf{q}_{i}^{'} \\mid} + \\sum_{i=2}^{N} \\frac{Gm_{i}^{'}M_{i}}{\\mid \\textbf{q}_{i}^{'} \\mid} - \\sum_{i=1}^{N} \\sum_{j=i+1}^{N} \\frac{Gm_{i}m_{j}}{\\mid \\textbf{q}_{i} - \\textbf{q}_{j} \\mid}
```

# Args

* `r`: NxD matrix, where N is the number of bodies and D the number of dimensions, it contains the initial positions of each body.
* `v`: NxD matrix. It contains the initial velocities of each body.
* `m`: array of N elements. It contains the mass of each body.
* `sld`: simulation logic dictionary, a dictionary where all posible options of the simulation can be specified. The following keys can be defined:
  * `bool_energy`: if true, energy of the hamiltonian will be calculated every output step.
  * `bool_print_kep`: information about the kepler-flow will be displayed in the console.
  * `ignore_H0`: boolean that if true, the linear movement of the system will be taken into account.
  * `out_step`: the output will be calculated on every 'out_step'-th step.
* `g`: gravitational constant.
* `t_step`: time that will be simulated in each step.
* `t_max`: specifies when the simulation will be stopped, after how much time.

# Returns

* `r_out`: OxNxD matrix, where O is the number of outputs, N the number of bodies and D the number of dimensions, it contains the positions calculated on specified time-steps of each body.
* `v_out`: OxNxD matrix, it contains the velocities calculated on specified time-steps of each body.
* `e_out`: array of O elements, it contains the energies of each output step of the whole N-body system.
"""
function simulation(r, v, m, sld::Dict{Any, Any}, g::Number = NaN, t_step::Number = NaN, t_max::Number = NaN)

    ### Analyze given input and calculate some variables
    exit = analyzeInput!(r, v, m, g, t_step, t_max, sld)
    if exit
        return
    end

    ### Allocate vectors we'll use during the execution of the program
    r_c = copy(r) # cartesian positions
    v_c = copy(v) # cartesian velocities

    r_j = zero(r) # jacobi positions
    v_j = zero(r) # jacobi velocities

    m_j = zero(m) # jacobi masses
    mu = zero(m) # standard gravitational parameters
    mu_sum = zero(m) # standard gravitational parameters of the sum of inner bodies
    m_sum = zero(m) # sum of masses of inner bodies

    Gz = zeros(typeof(r[1,1]),4) # g-functions will be saved here

    # if used values haven't got enough precision, define auxiliary varaibles for Kahan
    if eps(typeof(g)) > 1e-30
        r_kahan = zero(r)
        v_kahan = zero(r)
    end

    # if there are only 2 bodies, don't define the following
    if size(r)[1] > 2
        a_1 = zero(r)
        a_2 = zero(r)

        dp = zero(r[1,:])
        dp_sum = zero(r)

        if eps(typeof(g)) > 1e-30
            kick_kahan = zero(r)
        end
    end

    # auxiliary variables to use when calculating energy
    if sld["bool_energy"]
        p_j = zero(r)
        r_dist = zero(m)
    end

    # iteration count and step variables
    dt = t_step # full time step
    half_dt = dt/2.0 # half time step
    iter = 0 # current iteration
    max_iter = Int64(round(t_max/dt)) # number of iterations that the simulation will last
    steps_until_out = Int64(sld["out_step"]) # how many steps are needed to get an output
    index = 1 # how many outputs have been saved
    n_out = Int64(ceil(max_iter / steps_until_out)) + 1 # number of total outputs

    # output will be saved here
    r_out = zeros(typeof(r[1,1]), n_out, size(r)[1], size(r)[2])
    v_out = zeros(typeof(r[1,1]), n_out, size(r)[1], size(r)[2])
    if sld["bool_energy"]
        e_out = zeros(typeof(r[1,1]), n_out)
    end

    # calculate different mass types and standard gravitational parameter (mu)
    calculateMassMu!(m, g, m_j, mu, mu_sum, m_sum)

    # Calculate jacobi coordinates
    cartesian2jacobi!(r_c, r_j, m, m_sum)
    cartesian2jacobi!(v_c, v_j, m, m_sum)

    # In order to avoid calculating H_0 (or the movement of the whole system)
    if sld["ignore_H0"]
        v_j[1,:] -= v_j[1,:]
    end

    r_out[1,:,:] = r_c
    v_out[1,:,:] = v_c

    if sld["bool_energy"]
        for i in 1:size(r)[1]
            r_dist[i] = norm(r_j[i,:])
            p_j[i,:] = v_j[i,:]*m_j[i]
        end
        e_out[1] = calculateEnergy(r_c, r_j, r_dist, p_j, m, m_j, mu, mu_sum, sld["ignore_H0"])
    end

    ### drift the system under H0 for half a timestep dt
    if !sld["ignore_H0"]
        r_j[1,:] += v_j[1,:]*half_dt
    end

    ### drift all particles under H Kepler for half a timestep dt
    if eps(typeof(g)) <= 1e-30
        keplerStep!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"])
    else
        keplerStepWithKahan!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"], r_kahan, v_kahan)
    end
    total_dt = half_dt

    while iter < max_iter
        iter += 1
        steps_until_out -= 1

        # calculate 1st part of H Interaction in Jacobi coordinates
        if size(r)[1] > 2
            firstInteractionStep!(r_j, mu_sum, a_1)

            # update positions in the inertial frame
            jacobi2cartesian!(r_c, r_j, m, m_sum)

            # calculate 2nd part of H Interaction in inertial frame and convert them into Jacobi coordinates
            secondInteractionStep!(r_c, m, m_sum, m_j, mu, dp, dp_sum, a_2)

            # apply kick from Jacobi accelerations to Jacobi velocities
            if eps(typeof(g)) <= 1e-30
                v_j += dt * (a_1 + a_2)
            else
                kahanSum!(v_j, dt * (a_1 + a_2), kick_kahan)
            end
        end

        if iter != max_iter
            if steps_until_out == 0
                index += 1

                ### drift the system under H0 for half a timestep dt
                if !sld["ignore_H0"]
                    r_j[1,:] += v_j[1,:]*half_dt
                end

                ### drift all particles under H Kepler for half a timestep (half_dt)
                if eps(typeof(g)) <= 1e-30
                    keplerStep!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"])
                else
                    keplerStepWithKahan!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"], r_kahan, v_kahan)
                end
                total_dt += half_dt

                jacobi2cartesian!(r_c, r_j, m, m_sum)
                jacobi2cartesian!(v_c, v_j, m, m_sum)

                r_out[index,:,:] = r_c
                v_out[index,:,:] = v_c

                if sld["bool_energy"]
                    for i in 1:size(r)[1]
                        r_dist[i] = norm(r_j[i,:])
                        p_j[i,:] = v_j[i,:]*m_j[i]
                    end
                    e_out[index] = calculateEnergy(r_c, r_j, r_dist, p_j, m, m_j, mu, mu_sum, sld["ignore_H0"])
                end

                if !sld["ignore_H0"]
                    r_j[1,:] += v_j[1,:]*half_dt
                end

                if eps(typeof(g)) <= 1e-30
                    keplerStep!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"])
                else
                    keplerStepWithKahan!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"], r_kahan, v_kahan)
                end
                total_dt += half_dt

                steps_until_out = Int64(sld["out_step"])

            else
                ### drift the system under H0 for a full timestep dt
                if !sld["ignore_H0"]
                    r_j[1,:] += v_j[1,:]*dt
                end

                ### drift all particles under H Kepler for a full timestep (dt)
                if eps(typeof(g)) <= 1e-30
                    keplerStep!(r_j, v_j, dt, mu_sum, Gz, sld["bool_print"])
                else
                    keplerStepWithKahan!(r_j, v_j, dt, mu_sum, Gz, sld["bool_print"], r_kahan, v_kahan)
                end
                total_dt += dt
            end
        end
    end

    index += 1

    ### drift the system under H0 for half a timestep dt
    if !sld["ignore_H0"]
        r_j[1,:] += v_j[1,:]*half_dt
    end

    if eps(typeof(g)) <= 1e-30
        keplerStep!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"])
    else
        keplerStepWithKahan!(r_j, v_j, half_dt, mu_sum, Gz, sld["bool_print"], r_kahan, v_kahan)
    end
    total_dt += half_dt

    jacobi2cartesian!(r_c, r_j, m, m_sum)
    jacobi2cartesian!(v_c, v_j, m, m_sum)

    r_out[index,:,:] = r_c
    v_out[index,:,:] = v_c

    if sld["bool_energy"]
        for i in 1:size(r)[1]
            r_dist[i] = norm(r_j[i,:])
            p_j[i,:] = v_j[i,:]*m_j[i]
        end
        e_out[index] = calculateEnergy(r_c, r_j, r_dist, p_j, m, m_j, mu, mu_sum, sld["ignore_H0"])
    end

    if sld["bool_energy"]
        (r_out, v_out, e_out)
    else
        (r_out, v_out)
    end

end

max_stars_repo_name: salanueva/KeplerFlow

max_issues_repo_name: salanueva/KeplerFlow

max_forks_repo_name: salanueva/KeplerFlow


text:
using Test
using AtlasRobot
using RigidBodyDynamics
using RigidBodyDynamics.Contact: location

@testset "atlas" begin
    mechanism = AtlasRobot.mechanism()
    @test num_velocities(mechanism) == 36
    meshdir = joinpath(AtlasRobot.packagepath(), "Atlas", "urdf", "meshes")
    @test isdir(meshdir)
    @test count(x -> endswith(x, ".obj"), readdir(meshdir)) == 45
end

@testset "setnominal!" begin
    mechanism = AtlasRobot.mechanism(add_flat_ground=true)
    state = MechanismState(mechanism)
    AtlasRobot.setnominal!(state)
    for sideprefix in ('l', 'r')
        foot = findbody(mechanism, "$(sideprefix)_foot")
        for point in contact_points(foot)
            contact_location_world = transform(state, location(point), root_frame(mechanism))
            @test contact_location_world.v[3] ≈ 0 atol=2e-2
        end
    end
    @test center_of_mass(state).v[3] > 1
end


text:
using Random: AbstractRNG

"""
    hazards(::Function, process, ::AbstractRNG)

A process can pass transition distributions by creating a `hazards`
function that follows a visitor pattern. The first argument is a
callback function that takes three arguments: an identifier for a
clock, a distribution, and whether to enable or disable it.
"""
function hazards(::Function, process, ::AbstractRNG)
end

export hazards


text:
@testset "$TEST $G" begin

if !@isdefined(istree)
    function istree(parents::Vector{Int}, maxdepth)
        flag = true
        for i in 1:n
            s = i
            depth = 0
            while parents[s] > 0 && parents[s] != s
                s = parents[s]
                depth += 1
                if depth > maxdepth
                    return false
                end
            end
        end
        return flag
    end
end

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)

z = bfs_tree(g5, 1)
@test nv(z) == 4 && ne(z) == 3 && !has_edge(z, 2, 3)

g6 = graph(:house, G)
@test gdistances(g6, 2) == [1, 0, 2, 1, 2]
@test gdistances(g6, [1,2]) == [0, 0, 1, 1, 2]
@test !is_bipartite(g6)

g = G(5)
add_edge!(g,1,2); add_edge!(g,1,4)
add_edge!(g,2,3); add_edge!(g,2,5)
add_edge!(g,3,4)
@test is_bipartite(g)

g10 = CompleteGraph(10, G)
@test bipartite_map(g10) == Vector{Int}()

g10 = CompleteBipartiteGraph(10,10, G)
@test bipartite_map(g10) == Vector{Int}([ones(10); 2*ones(10)])

h10 = blockdiag(g10,g10)
@test bipartite_map(h10) == Vector{Int}([ones(10); 2*ones(10); ones(10); 2*ones(10)])

end # testset


text:
"""
N-dimensional ellipsoid. Note that only N=3 dimensional ellipsoids can be rendered
in the visualizer.
"""
struct HyperEllipsoid{N, T} <: GeometryPrimitive{N, T}
    center::Point{N, T}
    radii::Vec{N, T}
end

GeometryBasics.origin(geometry::HyperEllipsoid{N, T}) where {N, T} = geometry.center
radii(geometry::HyperEllipsoid{N, T}) where {N, T} = geometry.radii

"""
A collection of points which will be rendered as a point cloud in the visualiezr.

You can set the color of each point indepenedntly when you construct a PointCloud,
or you can set the overall point color by passing a `PointsMaterial` as the material
when you call `setobject!(vis, point_cloud, material)`
"""
struct PointCloud{T, PointType <: StaticVector{3, T}, C <: Colorant} <: AbstractGeometry{3, T}
    position::Vector{PointType}
    color::Vector{C}
end

function PointCloud(position::AbstractVector{<:AbstractVector{T}},
           color::AbstractVector{C}=RGB{Float32}[]) where {T, C <: Colorant}
    PointCloud{T, Point{3, T}, C}(position, color)
end

"""
A visual representation of the origin of a coordinate system, drawn as three
lines in red, green, and blue along the x, y, and z axes. The `scale` parameter
controls the length of the three lines.
"""
struct Triad <: AbstractGeometry{3, Float64}
    scale::Float64

    Triad(scale=1.0) = new(scale)
end

"""
A 3D cone, represented by the position of the origin of its base, the position of its tip,
and the radius of its base.
"""
struct Cone{N, T} <: AbstractGeometry{N, T}
    origin::Point{N, T}
    apex::Point{N, T}
    r::T
end

GeometryBasics.origin(geometry::Cone) = geometry.origin

center(geometry::HyperEllipsoid) = origin(geometry)
center(geometry::HyperRectangle) = minimum(geometry) + 0.5 * widths(geometry)
center(geometry::HyperSphere) = origin(geometry)
center(geometry::Cylinder) = (origin(geometry) + geometry.extremity) / 2
center(geometry::Cone) = (origin(geometry) + geometry.apex) / 2

"""
$(SIGNATURES)

Different tools disagree about what various geometric primitives mean. For example,
GeometryBasics.jl considers the "origin" of a cube to be its bottom-left corner, where
DrakeVisualizer and MeshCat consider its origin to be the center. The
intrinsic_transform(g) returns the transform from the GeometryBasics origin to the
MeshCat origin.
"""
intrinsic_transform(g) = IdentityTransformation()
intrinsic_transform(g::HyperRectangle) = Translation(center(g)...)
intrinsic_transform(g::HyperSphere) = Translation(center(g)...)
intrinsic_transform(g::HyperEllipsoid) = Translation(center(g)...) ∘ LinearMap(SMatrix{3, 3}(SDiagonal(radii(g)...)))

function intrinsic_transform(g::Cylinder{3})
    # Three.js wants a cylinder to lie along the y axis
    R = rotation_between(SVector(0, 1, 0), g.extremity - origin(g))
    Translation(center(g)) ∘ LinearMap(R)
end

function intrinsic_transform(g::Cone{3})
    # Three.js wants a cone to lie along the y axis
    R = rotation_between(SVector(0, 1, 0), g.apex - g.origin)
    Translation(center(g)) ∘ LinearMap(R)
end


text:
#ubox.jl - Ubox method

@universal function same_exp(x::Unum, y::Unum)
  is_negative(x) && is_positive(y) && return false
  (decode_exp(x) - (x.exponent == z64)) == (decode_exp(y) - (y.exponent == z64))
end

function isterminal{ESS,FSS}(x::Utype{ESS,FSS})
  isterminal(x.val)
end
function isterminal{ESS,FSS}(a::Vector{Utype{ESS,FSS}})
  for element in a
    isterminal(element) || return false
  end
  return true
end
@universal function isterminal(x::Unum)
  is_exact(x) && return true
  (x.fsize == max_fsize(FSS)) && return true
  return false
end
@universal function isterminal(x::Ubound)
  return false
end

upper_avg(x, y) = ((x + y) ÷ 2) + ((x + y < 0) ? 0 : 1) * isodd(x + y)

@universal function udivide(x::Unum)
  isterminal(x) && throw(ArgumentError("can't udivide an terminal unum"))

  if is_positive(x)
    #the low_ulp is the original unum, except with an extended ubit.
    low_ulp = copy(x)
    low_ulp.fsize += 1
    mid_exact = lub(low_ulp)
    high_ulp = copy(mid_exact)
    high_ulp.flags |= UNUM_UBIT_MASK
  else
    high_ulp = copy(x)
    high_ulp.fsize += 1
    mid_exact = glb(high_ulp)
    low_ulp = copy(mid_exact)
    low_ulp.flags |= UNUM_UBIT_MASK
  end

  (low_ulp, mid_exact, high_ulp)
end

function count_same(a::UInt64, b::UInt64)
  UInt16(leading_zeros(a $ b))
end

function count_same{FSS}(a::ArrayNum{FSS}, b::ArrayNum{FSS})
  total::UInt16 = 0
  for idx = 1:__cell_length(FSS)
    read = leading_zeros(a.a[idx] $ b.a[idx])
    total += read
    (read != 64) && return total
  end
end

function maskout(a::UInt64, count::UInt16)
  mask_top(count) & a
end
function maskout!{FSS}(a::ArrayNum{FSS}, count::UInt16)
  t = zero(ArrayNum{FSS})
  mask_top!(t, count)
  for idx = 1:__cell_length(FSS)
    a.a[idx] &= t.a[idx]
  end
end
@fracproc maskout count

function bitin(a::UInt64, count::UInt16)
  bottom_bit(count) | a
end
function bitin!{FSS}(a::ArrayNum{FSS}, count::UInt16)
  t = zero(ArrayNum{FSS})
  bottom_bit!(t, count)
  for idx = 1:__cell_length(FSS)
    a.a[idx] |= t.a[idx]
  end
end
@fracproc bitin count

@universal function udivide(x::Ubound)
  #special case for infinite bounds.
  if is_inf(x.lower)
    if is_inf(x.upper)
      middlebound = Ubound(neg_mmr(U), pos_mmr(U))
      return (x.lower, middlebound, x.upper)
    end

    return (x.lower, Ubound(neg_mmr(U), x.upper), nothing)
  elseif is_inf(x.upper)
    return (Ubound(x.lower, neg_mmr(U)), x.upper, nothing)
  end

  #similar special case for nearly infinite bounds.
  if is_neg_mmr(x.lower)
    if is_pos_mmr(x.upper)
      middlebound = Ubound(neg_big_exact(U), pos_big_exact(U))
      return (x.lower, middlebound, x.upper)
    end

    return (x.lower, Ubound(neg_big_exact(U), x.upper), nothing)
  elseif is_pos_mmr(x.upper)
    return (Ubound(x.lower, pos_big_exact(U)), x.upper, nothing)
  end

  #cases where zero is one of the bounds.
  if is_zero(x.lower)
    return (x.lower, Ubound(pos_sss(U), x.upper), nothing)
  elseif is_zero(x.upper)
    return (Ubound(x.lower, neg_sss(U)), x.upper, nothing)
  end

  #cases where we have exact values.
  if is_exact(x.lower)
    if is_exact(x.upper)
      return (x.lower, Ubound(upper_ulp(x.lower), lower_ulp(x.upper)), x.upper)
    end
    return (x.lower, Ubound(upper_ulp(x.lower), x.upper), nothing)
  elseif is_exact(x.upper)
    return (Ubound(x.lower, lower_ulp(x.upper)), x.upper, nothing)
  end

  #if the ubound straddles zero, then naturally cleave across zero.
  if (is_negative(x.lower) && is_positive(x.upper))
    return (resolve_as_utype!(copy(x.lower), neg_sss(Unum{ESS,FSS})),
            zero(Unum{ESS,FSS}),
            resolve_as_utype!(pos_sss(Unum{ESS,FSS}), copy(x.upper)))
  end

  #if the bounds are not exact...  First check if the exponents are the same
  if same_exp(x.lower, x.upper)
    #each exponential range cleaves over two single ubounds.  First, check to
    #see if they are in the upper half or the lower half.
    count::UInt16 = count_same(x.lower.fraction, x.upper.fraction)

    middle_value = copy(x.lower)
    middle_value.flags = x.lower.flags & (~UNUM_UBIT_MASK)
    frac_maskout!(middle_value, count)
    frac_bitin!(middle_value, count)

    #set the middle_value to look like
    return (resolve_as_utype!(copy(x.lower), lower_ulp(middle_value)), middle_value, resolve_as_utype!(upper_ulp(middle_value), x.upper))
  else
    #if they're not, then do a binary search on the exponents.
    middle_exp = upper_avg(decode_exp(x.lower), decode_exp(x.upper))
    middle_value = zero(U)
    middle_value.flags = x.lower.flags & UNUM_SIGN_MASK
    (middle_value.esize, middle_value.exponent) = encode_exp(middle_exp)

    lhs = Ubound(copy(x.lower), lower_ulp(middle_value))
    rhs = Ubound(upper_ulp(middle_value), copy(x.upper))

    return (lhs, middle_value, rhs)
  end
end

udivide(x::Utype) = map(Utype, udivide(x.val))

doc"""
  ufilter(f::Function, uvalue::Utype)
  ufilter(f::Function, uarray::Vector{Utype,1})

  passes an array of utype boxed objects and does a binary search to whittle down
  the results to a minimal array of solutions.  Function f should be a predicate
  fuction, aka, returns a boolean value.
"""
function ufilter{ESS,FSS, verbose}(f::Function, u::Utype{ESS,FSS}, ::Type{Val{verbose}} = Val{false})
  #first create an anonymous function g which takes an array of Utypes and passes
  #the first member to f.
  g = (a::Vector{Utype{ESS, FSS}}) -> f(a[1])
  #then execute the g function on the array ufilter method.
  ufilter(g, [u], Val{verbose})
end

const max_iters = 2
iters = 0

function ufilter{ESS,FSS, verbose}(f::Function, v::Vector{Utype{ESS,FSS}}, ::Type{Val{verbose}} = Val{false})

  global iters

  l = length(v)
  result = Array{Utype,2}(l, 0)

  #trigger a terminal check.
  isterminal(v) && return f(v) ? v : Matrix{Utype{ESS,FSS}}(l,0)

  #pick a random index that isn't terminal.
  validindices = filter((idx)->(!isterminal(v[idx])), 1:l)
  ridx = rand(validindices)

  (lower_u, middle_u, upper_u) = udivide(v[ridx])

  vl = copy(v)
  vm = copy(v)

  vl[ridx] = lower_u
  vm[ridx] = middle_u

  if f(vl)
    verbose && (print("searching "); describe.(vl))
    iters == max_iters && exit()
    result = hcat(result, ufilter(f, vl, Val{verbose}))
  end

  if f(vm)
    verbose && (print("searching "); describe.(vm))
    iters == max_iters && exit()
    (result = hcat(result, ufilter(f, vm, Val{verbose})))
  end

  if (upper_u != nothing)
    vu = copy(v)
    vu[ridx] = upper_u
    if f(vu)
      verbose && (print("searching "); describe.(vu))
      iters == max_iters && exit()
      (result = hcat(result, ufilter(f, vu, Val{verbose})))
    end
  end

  result
end

export ufilter

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
module CMUDict

export cmudict

using ..Data: deps

const version = "0.7b"
const cache_prefix = "https://cache.julialang.org"

function load()
  suffixes = ["", ".phones", ".symbols"]
  if isdir(deps("cmudict"))
    if all(isfile(deps("cmudict", "cmudict$x")) for x in suffixes)
      return
    end
  end
  info("Downloading CMUDict dataset")
  mkpath(deps("cmudict"))
  for x in suffixes
    download("$cache_prefix/http://svn.code.sf.net/p/cmusphinx/code/trunk/cmudict/cmudict-$version$x",
             deps("cmudict", "cmudict$x"))
  end
end

function phones()
  load()
  Symbol.(first.(split.(split(readstring(deps("cmudict", "cmudict.phones")),
                        "\n", keep = false), "\t")))
end

function symbols()
  load()
  Symbol.(split(readstring(deps("cmudict", "cmudict.symbols")),
                "\n", keep = false))
end

function rawdict()
  load()
  Dict(String(xs[1]) => Symbol.(xs[2:end]) for xs in
       filter(!isempty, split.(split(readstring(deps("cmudict", "cmudict")), "\n"))))
end

validword(s) = ismatch(r"^[\w\-\.]+$", s)

cmudict() = filter((s, ps) -> validword(s), rawdict())

alphabet() = ['A':'Z'..., '0':'9'..., '_', '-', '.']

end


text:
using ChemistryFeaturization, AtomGraphs
using ChemistryFeaturization.ElementFeature
using CSV, DataFrames

# custom lookup table...
test_df = CSV.read(abspath(@__DIR__, "test_data", "lookup_table.csv"), DataFrame)

fnames = ["X", "Block", "Atomic mass"]
fs = ElementFeatureDescriptor.(fnames)

triangle_C = AtomGraph(Float32.([0 1 1; 1 0 1; 1 1 0]), ["C", "C", "C"])

@testset "Encode/Featurize" begin
    using ChemistryFeaturization.ElementFeature

    # make sure both constructors give the same results
    fzn1 = GraphNodeFeaturization(fs)
    fzn2 = GraphNodeFeaturization(fnames)
    fzn3 = GraphNodeFeaturization(fnames, nbins = [6, 4, 10])

    @test output_shape(fzn1) == 24 # 10 (default length for X) + 4 (default for block) + 10 (default for atomic mass)
    @test output_shape(fzn2) == 24
    @test output_shape(fzn3) == 20

    encoded_1, encoded_2 = encode.([triangle_C, triangle_C], Ref(fzn1)) # encode can be broadcasted
    @test encoded_1 == encoded_2

    encoded_1, encoded_2 = encode.([triangle_C, triangle_C], [fzn1, fzn2]) # encode can be broadcasted
    @test encoded_1 == encoded_2

    featurized_1, featurized_2 = featurize.([triangle_C, triangle_C], [fzn1, fzn2])
    @test featurized_1.encoded_features == featurized_2.encoded_features
end

@testset "Decode" begin
    fzn1 = GraphNodeFeaturization(fs)
    fa = featurize(triangle_C, fzn1)
    decoded = decode(fa)

    @test all(map(i -> decoded[i]["Block"] == "p", 1:3))
    @test all(map(i -> decoded[i]["X"][1] <= 2.55 <= decoded[i]["X"][2], 1:3))
end

# encodable_elements
@testset "Encodable Elements" begin
    fzn4 = GraphNodeFeaturization(ElementFeatureDescriptor.(["Boiling point", "6d"]))
    fzn5 = GraphNodeFeaturization(ElementFeatureDescriptor.(["7s", "Valence"]))
    @test encodable_elements(fzn4) == ["Ac", "Th", "U"]
    @test encodable_elements(fzn5) == ["Fr", "Ra", "Ac", "Th"]

    # Custom lookup_table
    feature_1 = ElementFeatureDescriptor("MeaningOfLife", test_df)   # zero-value case - `As` has a value = 0
    feature_2 = ElementFeatureDescriptor("neg_nums", test_df)
    feature_3 = ElementFeatureDescriptor("first_letter", test_df)
    feature_4 = ElementFeatureDescriptor("noarsenic", test_df) # missing value case - `As` has a missing value

    @test encodable_elements(GraphNodeFeaturization([feature_1, feature_2])) ==
          ["C", "As", "Tc"]
    @test encodable_elements(GraphNodeFeaturization([feature_1, feature_4])) == ["C", "Tc"]
    @test encodable_elements(GraphNodeFeaturization([feature_2, feature_3])) ==
          ["C", "As", "Tc"]
    @test encodable_elements(GraphNodeFeaturization([feature_2, feature_4])) == ["C", "Tc"]
end

# chunk_vec helper fcn
@testset "chunk_vec" begin
    vec = [1, 1, 0, 1, 0, 1, 0]
    @test_throws AssertionError chunk_vec(vec, [3, 3])
    @test chunk_vec(vec, [4, 1, 2]) == [[1, 1, 0, 1], [0], [1, 0]]
end


text:
function walk_nodes(ind::CGPInd)
    dot_strs = []
    push!(dot_strs, "digraph cgpgraph {\n")
    function visit(node_num, current_nd)
        if current_nd.active
            fname = String(Symbol(current_nd.f))
            visit(current_nd.x, ind.nodes[current_nd.x])
            if(CGPFunctions.arity[fname] > 1)
               visit(current_nd.y, ind.nodes[current_nd.y])

            end
            @show (node_num, current_nd)

            xdotline = "N$(current_nd.x)-> N$node_num ;"
            push!(dot_strs, "N$node_num [label=\"$fname:$node_num\"];")
            push!(dot_strs, xdotline)
            if(CGPFunctions.arity[fname] > 1)
               ydotline = "N$(current_nd.y)-> N$node_num ;"
               push!(dot_strs, ydotline)
            end
        end
    end
    for current_node_num in ind.outputs
        current_node     = ind.nodes[current_node_num]
        visit(current_node_num, current_node)
    end
    push!(dot_strs, "}")
    dot_strs
end


text:
using FunManifolds
using Test
using ForwardDiff
using StaticArrays

include("utils.jl")

@testset "FunManifolds tests" begin
    include("manifolds/DiscretizedCurves.jl")
    include("manifolds/CurveWarpingSRSFSpace.jl")
    include("manifolds/FunctionCurve.jl")
    include("curve_warping.jl")
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "Load a small PSDD and test methods" begin
   file = zoo_psdd_file("little_4var.psdd")
   prob_circuit = load_prob_circuit(file);
   @test prob_circuit isa ProbΔ

   # Testing number of nodes and parameters
   @test  9 == num_parameters(prob_circuit)
   @test 20 == size(prob_circuit)[1]

   # Testing Read Parameters
   EPS = 1e-7
   @test abs(prob_circuit[13].log_thetas[1] - (-1.6094379124341003)) < EPS
   @test abs(prob_circuit[13].log_thetas[2] - (-1.2039728043259361)) < EPS
   @test abs(prob_circuit[13].log_thetas[3] - (-0.916290731874155)) < EPS
   @test abs(prob_circuit[13].log_thetas[4] - (-2.3025850929940455)) < EPS

   @test abs(prob_circuit[18].log_thetas[1] - (-2.3025850929940455)) < EPS
   @test abs(prob_circuit[18].log_thetas[2] - (-2.3025850929940455)) < EPS
   @test abs(prob_circuit[18].log_thetas[3] - (-2.3025850929940455)) < EPS
   @test abs(prob_circuit[18].log_thetas[4] - (-0.35667494393873245)) < EPS

   @test abs(prob_circuit[20].log_thetas[1] - (0.0)) < EPS
end

psdd_files = ["little_4var.psdd", "msnbc-yitao-a.psdd", "msnbc-yitao-b.psdd", "msnbc-yitao-c.psdd", "msnbc-yitao-d.psdd", "msnbc-yitao-e.psdd", "mnist-antonio.psdd"]

@testset "Test parameter integrity of loaded PSDDs" begin
   for psdd_file in psdd_files
      @test check_parameter_integrity(load_prob_circuit(zoo_psdd_file(psdd_file)))
   end
end

@testset "Test parameter integrity of loaded structured PSDDs" begin
   circuit, vtree = load_struct_prob_circuit(
      zoo_psdd_file("little_4var.psdd"), zoo_vtree_file("little_4var.vtree"))
   @test check_parameter_integrity(circuit)
   @test vtree isa PlainVtree
   # no other combinations of vtree and psdd are in this repo?
   # @test check_parameter_integrity(load_struct_prob_circuit(
   #          "circuits/mnist-antonio.psdd", "circuits/balanced.vtree"))
end


text:
using ImmutableArrays

function normal(a, b)
    dx = b[1] - a[1]
    dy = b[2] - a[2]
    c = sqrt(dx * dx + dy * dy)
    return Vector2(dx/c, dy/c)
end
function dist(a, b)
    dx = b[1] - a[1]
    dy = b[2] - a[2]
    c = sqrt(dx * dx + dy * dy)
    return c
end

typealias Coordinate Vector2{Int16}

immutable PointElement
    x::Uint16
end

immutable Triangle{T <: Integer} # for this purpose uint16 is used
    a::T
    b::T
    c::T
end
immutable LineEdge
    x::Int16
    y::Int16
    ex::Int8
    ey::Int8
    linesofar1::Int8
    linesofar2::Int8
    const extrudeScale = Int8(63)
    function LineEdge(x::Int16, y::Int16, ex::Float32, ey::Float32, tx::Int8, ty::Int8, linesofar::Int32)
        new(
            (x * 2) | tx,
            (y * 2) | ty,
            round(extrudeScale * ex),
            round(extrudeScale * ey),
            linesofar / 128,
            linesofar % 128,
        )
    end
end

function add(buffer::Vector{LineEdge}, x::Int16, y::Int16, ex::Float32, ey::Float32, tx::Int8, ty::Int8, linesofar::Int32 = 0)
    push!(buffer, LineEdge(x, y, ex, ey, tx, ty, linesofar))
    return length(buffer)
end
function add{T}(buffer::Vector{Triangle{T}}, a::T, b::T, c::T)
    push!(buffer, Triangle(a,b,c))
    return length(buffer)
end
function add(buffer::Vector{PointElement}, a::Uint16)
    push!(buffer, PointElement(a))
    return length(buffer)
end

immutable ElementGroup{N, T}
    array::NTuple{N, T}
    vertex_length::Uint32
    elements_length::Uint32

    ElementGroup() = new(0, 0)
    ElementGroup(vertex_length::Uint32, elements_length::Uint32) = new(vertex_length, elements_length)
end
const triangleGroups = Ptr{Triangle{Uint16}}[]
const pointGroups    = Ptr{PointElement}[]

function addGeometry(vertices)
    # TODO: use roundLimit
    # const float roundLimit = geometry.round_limit

    if length(vertices) < 2
        error("a line must have at least two vertices\n")
    end

    firstVertex = first(vertices)
    lastVertex  = last(vertices)
    closed      = firstVertex[1] == lastVertex[1] && firstVertex[2] == lastVertex[2]

    if length(vertices) == 2 && closed
        # fprintf(stderr, "a line may not have coincident points\n")
        return
    end

    beginCap    = properties.cap
    endCap      = closed ? CapType.Butt : properties.cap

    currentJoin = JoinType.Miter

    currentVertex   = Coordinate(0, 0)
    prevVertex      = Coordinate(0, 0)
    nextVertex      = Coordinate(0, 0)
    prevNormal      = Vector(0.0)
    nextNormal      = Vector(0.0)

    EType       = Int32
    e1          = -one(EType)
    e2          = -one(EType)
    e3          = -one(EType)

    flip        = one(Int8)
    distance    = 0.0

    if (closed)
        currentVertex   = vertices[end - 2]
        nextNormal      = normal(currentVertex, lastVertex)
    end

    start_vertex = convert(Int32, length(vertexBuffer))

    triangle_store  = Triangle{Uint16}[]
    point_store     = PointElement[]

    for i=1:length(vertices)

        if nextNormal
            prevNormal = Vec2(-nextNormal[1], -nextNormal[2])
        end
        if currentVertex
            prevVertex = currentVertex
        end

        currentVertex   = vertices[i]
        currentJoin     = properties.join

        if prevVertex
            distance += dist(currentVertex, prevVertex)
        end

        # Find the next vertex.
        if i + 1 < length(vertices)
            nextVertex = vertices[i + 1]
        else
            nextVertex = Coordinate(0, 0)
        end

        # If the line is closed, we treat the last vertex like the first vertex.
        if (!nextVertex && closed)
            nextVertex = vertices[1]
        end

        if (nextVertex)
            # if two consecutive vertices exist, skip one
            (currentVertex[1] == nextVertex[1] && currentVertex[2] == nextVertex[2]) && continue
        end

        # Calculate the normal towards the next vertex in this line. In case
        # there is no next vertex, pretend that the line is continuing straight,
        # meaning that we are just reversing the previous normal
        if (nextVertex)
            nextNormal = normal(currentVertex, nextVertex)
        else
            nextNormal = Vector2(-prevNormal[1], -prevNormal[2])
        end

        # If we still don't have a previous normal, this is the beginning of a
        # non-closed line, so we're doing a straight "join".
        if (!prevNormal)
            prevNormal = Vector2(-nextNormal[1], -nextNormal[2])
        end

        # Determine the normal of the join extrusion. It is the angle bisector
        # of the segments between the previous line and the next line.
        joinNormal = Vector2(
            prevNormal[1] + nextNormal[1],
            prevNormal[2] + nextNormal[2]
        )

        # Cross product yields 0..1 depending on whether they are parallel
        # or perpendicular.
        joinAngularity = nextNormal[1] * joinNormal[2] - nextNormal[2] * joinNormal[1]
        joinNormal[1] /= joinAngularity
        joinNormal[2] /= joinAngularity
        roundness::Float64 = max(abs(joinNormal[1]), abs(joinNormal[2]))

        # Switch to miter joins if the angle is very low.
        if (currentJoin != JoinType.Miter)
            if (abs(joinAngularity) < 0.5 && roundness < properties.miter_limit)
                currentJoin = JoinType.Miter
            end
        end

        # Add offset square begin cap.
        if (!prevVertex && beginCap == CapType.Square)
            # Add first vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   flip * (prevNormal[1] + prevNormal[2]), flip * (-prevNormal[1] + prevNormal[2]), # extrude normal
                                   0, 0, distance) - start_vertex)::Int32  # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))
            e1 = e2
            e2 = e3

            # Add second vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   flip * (prevNormal[1] - prevNormal[2]), flip * (prevNormal[1] + prevNormal[2]), # extrude normal
                                   0, 1, distance) - start_vertex)::Int32  # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))
            e1 = e2
            e2 = e3
        # Add offset square end cap.
        elseif !nextVertex && endCap == CapType.Square
            # Add first vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   nextNormal[1] - flip * nextNormal[2], flip * nextNormal[1] + nextNormal[2], # extrude normal
                                   0, 0, distance) - start_vertex)::Int32 # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))

            e1 = e2::Int32
            e2 = e3::Int32

            # Add second vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   nextNormal[1] + flip * nextNormal[2], -flip * nextNormal[1] + nextNormal[2], # extrude normal
                                   0, 1, distance) - start_vertex)::Int32  # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))
            e1 = e2::Int32
            e2 = e3::Int32

        elseif currentJoin == JoinType.Miter
            # MITER JOIN
            if (abs(joinAngularity) < 0.01)
                # The two normals are almost parallel.
                joinNormal[1] = -nextNormal[2]
                joinNormal[2] = nextNormal[1]
            elseif roundness > properties.miter_limit
                # If the miter grows too large, flip the direction to make a
                # bevel join.
                joinNormal[1] = (prevNormal[1] - nextNormal[1]) / joinAngularity
                joinNormal[2] = (prevNormal[2] - nextNormal[2]) / joinAngularity
            end

            if roundness > properties.miter_limit
                flip = -flip
            end

            # Add first vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   flip * joinNormal[1], flip * joinNormal[2], # extrude normal
                                   0, 0, distance) - start_vertex)::Int32  # texture normal

            e1 >= 0 && e2 >= 0 && e3 >= 0 && push!(triangle_store, Triangle(e1, e2, e3))

            e1 = e2::Int32
            e2 = e3::Int32

            # Add second vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   -flip * joinNormal[1], -flip * joinNormal[2], # extrude normal
                                   0, 1, distance) - start_vertex)::Int32  # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))

            e1 = e2::Int32
            e2 = e3::Int32

            if ((!prevVertex && beginCap == CapType.Round) ||
                    (!nextVertex && endCap == CapType.Round))
                push!(point_store, PointElement(e1))
            end
        else
            # Close up the previous line
            # Add first vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   flip * prevNormal[2], -flip * prevNormal[1], # extrude normal
                                   0, 0, distance) - start_vertex)::Int32  # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))
            e1 = e2::Int32
            e2 = e3::Int32

            # Add second vertex.
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   -flip * prevNormal[2], flip * prevNormal[1], # extrude normal
                                   0, 1, distance) - start_vertex)::Int32  # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))
            e1 = e2::Int32
            e2 = e3::Int32

            prevNormal = Vector2(-nextNormal[1], -nextNormal[2])
            flip = 1

            # begin/end caps
            if ((!prevVertex && beginCap == CapType.Round) ||
                    (!nextVertex && endCap == CapType.Round))
                push!(point_store, PointElement(e1))
            end

            if (currentJoin == JoinType.Round)
                if (prevVertex && nextVertex && (!closed || i > 0))
                    push!(point_store, PointElement(e1))
                end

                # Reset the previous vertices so that we don't accidentally create
                # any triangles.
                e1 = -one(Int32)
                e2 = -one(Int32)
                e3 = -one(Int32)
            end

            # Start the new quad.
            # Add first vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   -flip * nextNormal[2], flip * nextNormal[1], # extrude normal
                                   0, 0, distance) - start_vertex)::Int32 # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))

            e1 = e2::Int32
            e2 = e3::Int32

            # Add second vertex
            e3 = (add(vertexBuffer, currentVertex[1], currentVertex[2], # vertex pos
                                   flip * nextNormal[2], -flip * nextNormal[1], # extrude normal
                                   0, 1, distance) - start_vertex)::Int32 # texture normal

            (e1 >= 0 && e2 >= 0 && e3 >= 0) && push!(triangle_store, Triangle(e1, e2, e3))

            e1 = e2::Int32
            e2 = e3::Int32
        end
    end

    end_vertex      = length(vertexBuffer)
    vertex_count    = end_vertex - start_vertex

    # Store the triangle/line groups.

        if (!length(triangleGroups) || (last(triangleGroups).vertex_length + vertex_count > 65535))
            # Move to a new group because the old one can't hold the geometry.
            triangleGroups.emplace_back()
        end

        group = last(triangleGroups)
        for triangle in triangle_store
            add(triangleElementsBuffer,
                group.vertex_length + triangle.a,
                group.vertex_length + triangle.b,
                group.vertex_length + triangle.c
            )
        end

        group.vertex_length += vertex_count
        group.elements_length += length(triangle_store)
    end

    # Store the line join/cap groups.

        if (!length(pointGroups) || (last(pointGroups).vertex_length + vertex_count > 65535))
            # Move to a new group because the old one can't hold the geometry.
            pointGroups.emplace_back()
        end

        group = last(pointGroups)
        for point in point_store
            add(pointElementsBuffer, group.vertex_length + point)
        end

        group.vertex_length += vertex_count
        group.elements_length += length(point_store)
    end
end


text:
name = "libLLVM"
LLVM_full_version = v"8.0.1+3"

# Include common tools
include("../common.jl")
build_tarballs(ARGS, configure_extraction(ARGS, LLVM_full_version, name)...; skip_audit=true)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "continuous" begin
    prec = 128
    RR = RealField(prec)
    atol = sqrt(eps())
    rtol = sqrt(eps())

    problems = [(sin, 0, 1, sin(RR(1))),
                (sinpi, 0, 1, 1),
                (x -> 2x + 1, 0, 1, 3),
                (x -> -(x + RR("1e-1"))^3 , 0, 1, RR("-1e-3"))]

    for (f, a, b, value) in problems
        a = RR(a)
        b = RR(b)
        value = RR(value)

        for evaltype in (:ball, :taylor)
            enclosure = enclosemaximum(f, a, b, atol = atol, rtol = rtol, evaltype = evaltype)

            @test contains(enclosure, value)
            if !contains(enclosure, value)
                @show enclosure value
            end
            @test radius(value) < atol || radius(value)/abs(value) < rtol
        end
    end
end

@testset "unbounded" begin
    prec = 128
    RR = RealField(prec)

    a = RR(0)
    b = RR(1)

    @test isnan(enclosemaximum(tanpi, a, b))
    @test isnan(enclosemaximum(x -> 1/x, a, b))
end

@testset "discontinuous" begin
    prec = 128
    RR = RealField(prec)
    atol = sqrt(eps())
    rtol = sqrt(eps())

    a = RR(0)
    b = RR(1)

    enclosure = enclosemaximum(a, b, atol = atol, rtol = rtol) do x
        if x < 0.5
            return sin(x)
        elseif x > 0.5
            return cos(x)
        end

        setunion(sin(x), cos(x))
    end
    value = cos(RR(0.5))

    @test contains(enclosure, value)
    @test radius(value) < atol || radius(value)/abs(value) < rtol
end

@testset "taylor" begin
    prec = 128
    RR = RealField(prec)
    atol = sqrt(eps())
    rtol = sqrt(eps())

    a = RR(0)
    b = RR(1)

    enclosure = enclosemaximum(sin, a, b, atol = atol, rtol = rtol, evaltype = :taylor)
    value = sin(b)
    @test contains(enclosure, value)
    @test radius(value) < atol || radius(value)/abs(value) < rtol
end


text:
struct SURE <: AbstractEBayesMethod end

mutable struct FittedEBayes{SS,
                            L,
                            FL,
                            PRED,
                            EBM,
                            FEBM,
                            O}
    ss::SS
    loc::L
    fitted_loc::FL
    predictor::PRED
    method::EBM
    fitted_obj::FEBM
    report::O #Typically nothing or result from Optim.jl
end

MLJBase.fitted_params(feb::FittedEBayes) = (feb.fitted_obj, feb.fitted_loc)

function predict(feb::FittedEBayes)
    prior_μs = predict(feb.loc, feb.fitted_loc, feb.ss)
    predict(feb.predictor, feb.fitted_obj, prior_μs, feb.ss)
end

# TODO: add nicer default interface

function fit(predictor::Normal, method::SURE, loc::EBayesShrinkageLocation, ss::NormalSamples;
             optimizer = IPNewton())
    fitted_loc = fit(loc, ss)
    x0 = optim_start_point(predictor, method, loc, fitted_loc, ss)
    f_diff = optim_objective(predictor, method, loc, fitted_loc, ss, x0)
    dfc = optim_constraints(predictor, method, loc, fitted_loc, ss, x0)
    opt_res = optimize(f_diff, dfc, x0, optimizer)
    fitted_obj = optim_extract(predictor, method, loc, fitted_loc, ss, opt_res)
    FittedEBayes(ss, loc, fitted_loc,
                 predictor, method, fitted_obj,
                 opt_res)
end

function fit(predictor::Normal, method::SURE, ss::NormalSamples; optimizer=IPNewton())
    fit(predictor, method, GrandMeanLocation(), ss; optimizer=optimizer)
end

#------------------------------------------------------------------------------------------
#----------------- Normal prior with with SURE --------------------------------------------
#------------------------------------------------------------------------------------------

function optim_start_point(predictor::Normal,
                           sure::SURE,
                           loc::EBayesShrinkageLocation,
                           fitted_loc,
                           ss::NormalSamples)
    μs = predict(loc, fitted_loc, ss)
    Zs = response(ss)
    σs_squared = var(ss)

    x0  = sqrt(max( mean((Zs .- μs).^2 .- σs_squared), 0.0001))
    x0 = [x0]
    x0
end

function optim_constraints(predictor::Normal,
                           sure::SURE,
                           loc::EBayesShrinkageLocation,
                           fitted_loc,
                           ss::NormalSamples,
                           x0)
    dfc = TwiceDifferentiableConstraints([0.0], [Inf])
    dfc
end

function _f_sure(A, Z, μ, h_ii, σ_squared)
    #(σ^2/(A^2 + σ^2))^2*(Z - μ)^2 + σ^2 - 2*(1 - h_ii)*σ^4/(A^2 + σ^2)
    (σ_squared/(A + σ_squared))^2*(Z - μ)^2 - 2*(1 - h_ii)*σ_squared^2/(A + σ_squared)
end

function optim_objective(predictor::Normal,
                         sure::SURE,
                         loc::EBayesShrinkageLocation,
                         fitted_loc,
                         ss::NormalSamples,
                         x0)
    #f_sure(A, Z, σ) = (σ^2/(A^2 + σ^2))^2*Z^2 + σ^2*(A^2 - σ^2)/(A^2 + σ^2)
                                        #             + σ^2 - 2σ^4/(A^2 + σ^2)
    #f_sure(A, Z, σ) = (σ^2/(A^2 + σ^2))^2*(Z - Zs_bar)^2 + σ^2*(A^2 - σ^2 + 2*σ^2/n)/(A^2 + σ^2)

    Zs = response(ss)
    σs_squared = var(ss)
    μs = predict(loc, fitted_loc, ss)
    hat_diags = leverage(loc, fitted_loc, ss)

    f_diff = TwiceDifferentiable(A -> mean( _f_sure.(first(A), Zs, μs,
                                                     hat_diags, σs_squared)),
                                 x0; autodiff=:forward)
end

function optim_extract(predictor::Normal,
                       sure::SURE,
                       loc::EBayesShrinkageLocation,
                       fitted_loc,
                       ss::NormalSamples,
                       res)
    A = first(Optim.minimizer(res))
    A
end

#function optim_extract(s::NormalSamples,
#                       sure::SURE{D, FixedLocation{T}},
#                       res) where {D<:UnivariateDistribution, T}
#    σ_prior = first(Optim.minimizer(res))
#    D_type = typeof(sure.d)
#    D_type.(sure.locs.μs, σ_prior)
#end

#function optim_extract(s::NormalSamples,
#                sure::SURE{Normal{T}, GrandMeanLocation{S, W}},
#                res) where {T,S,W}

#    σ_prior = first(Optim.minimizer(res))
#    Zs_bar = mean(s.Z)
#    Normal.(Zs_bar .+ sure.locs.μs, σ_prior)
#end

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

#function fit(sure::SURE{Normal{T}, FixedLocation{T}},
#             s::HomoskedasticNormalSamples) where T
#     σ_squared_lik = s.σ^2
#     Zs = s.Zs
#     μ = sure.locs.μs
#     σ_prior_squared = max( mean((Zs .- μ).^2) - σ_squared_lik, zero(T))
#     EmpiricalBayesProblem(s, Normal.(μ, sqrt(σ_prior_squared)), sure)
#end

#----------------------------------------------------------------------
# Semiparametric SURE
#----------------------------------------------------------------------

#(1-b_i)Z_i + b_i \bar{Z}
#    `` \sum_{i=1}^n b_i^2 (Z_i - \bar{Z})^2 + (1 - 2(1-1/n)b_i))\sigma_i^2 ``
#
#  so (Z_i - \bar{Z})^2 in front of b_i^2 and
# -2\sigma^2(1-1/n) in front of b...
# hence could use      \sigma^2(1-1/n)/(Z_i - \bar{Z})^2 as my pseudo-observation.

#function fit(sure::SURE{SemiparametricLinearShrinker,
#                        GrandMeanLocation{T, Nothing}},
#             s::NormalSamples) where T
#    n = length(s)
#    Zs = s.Z
#    σs = s.σ
#    Zbar = mean(Zs)

#    ws = (Zs .- Zbar).^2
#    Y_pseudo = replace(ws, 0.0 => 10.0) #will get 0 weight in optim anyway

#    Y_pseudo = σs.^2 .* (1 - 1/n) ./ Y_pseudo

#    iso_fit = fit(IsotonicRegression, σs, Y_pseudo; wts=ws, lb=0.0, ub=1.0)
#    bs = predict(iso_fit)

#    linshrink = LinearShrinker(Zbar,  1 .- bs)
#    EBayesOptimResult(s, linshrink, sure)
#end

# More or less the same
#   \sum_{i=1}^n b_i^2 (Z_i - \mu_i)^2 + (1 - 2b_i))\sigma_i^2 `


text:
Base.dataids(sa::HybridArray) = Base.dataids(sa.data)

@inline size(sa::HybridArray{S,T,N,N}) where {S,T,N} = size(sa.data)

@inline length(sa::HybridArray) = length(sa.data)

@inline strides(sa::HybridArray{S,T,N,N}) where {S,T,N} = strides(sa.data)

@inline pointer(sa::HybridArray) = pointer(sa.data)

@generated function _sized_abstract_array_axes(::Type{S}, ax::Tuple) where S<:Tuple
    exprs = Any[]
    map(enumerate(S.parameters)) do (i, si)
        if isa(si, Dynamic)
            push!(exprs, :(ax[$i]))
        else
            push!(exprs, SOneTo(si))
        end
    end
    return Expr(:tuple, exprs...)
end

function axes(sa::HybridArray{S}) where S
    ax = axes(sa.data)
    return _sized_abstract_array_axes(S, ax)
end

function promote_rule(::Type{<:HybridArray{S,T,N,M,TDataA}}, ::Type{<:HybridArray{S,U,N,M,TDataB}}) where {S,T,U,N,M,TDataA,TDataB}
    TU = promote_type(T,U)
    HybridArray{S,TU,N,M,promote_type(TDataA, TDataB)::Type{<:AbstractArray{TU}}}
end

@inline copy(a::HybridArray) = typeof(a)(copy(a.data))

similar(a::HA, ::Type{T2}) where {S, HA<:HybridArray{S}, T2} = HybridArray{S, T2}(similar(a.data, T2))

similar(::Type{<:HybridArray{S,T,N,M}},::Type{T2}) where {S,T,N,M,T2} = HybridArray{S,T2,N,M}(undef)
similar(::Type{SA},::Type{T},s::Size{S}) where {SA<:HybridArray,T,S} = hybridarray_similar_type(T,s,StaticArrays.length_val(s))(undef)
hybridarray_similar_type(::Type{T},s::Size{S},::Type{Val{D}}) where {T,S,D} = HybridArray{Tuple{S...},T,D,length(s)}

# for internal use only (used in vcat and hcat)
# TODO: try to make this less hacky
# adding method to similar_type ends up being used in broadcasting for some reason?
_h_similar_type(::Type{A},::Type{T},s::Size{S}) where {A<:HybridArray,T,S} = hybridarray_similar_type(T,s,StaticArrays.length_val(s))

Size(::Type{<:HybridArray{S}}) where {S} = Size(S)

Base.IndexStyle(a::HybridArray) = Base.IndexStyle(a.data)
Base.IndexStyle(::Type{HA}) where {S,T,N,M,TData,HA<:HybridArray{S,T,N,M,TData}} = Base.IndexStyle(TData)

Base.vec(a::HybridArray) = vec(a.data)

StaticArrays.similar_type(::Type{<:SSubArray},::Type{T},s::Size{S}) where {S,T} = StaticArrays.default_similar_type(T,s,StaticArrays.length_val(s))


text:
module Vcov
using LinearAlgebra
using StatsBase
using Combinatorics
using FixedEffects
using DataFrames

##############################################################################
##
## Any DataType must define the following methods:
##
##############################################################################
struct VcovData{T, N} <: RegressionModel
    modelmatrix::Matrix{Float64}       # X
    crossmodelmatrix::T                    # X'X in the simplest case. Can be Matrix but preferably Factorization
    residuals::Array{Float64, N}      # vector or matrix of residuals (matrix in the case of IV, residuals of Xendo on (Z, Xexo))
    dof_residual::Int
end
StatsBase.modelmatrix(x::VcovData) = x.modelmatrix
StatsBase.crossmodelmatrix(x::VcovData) = x.crossmodelmatrix
StatsBase.residuals(x::VcovData) = x.residuals
StatsBase.dof_residual(x::VcovData) = x.dof_residual

##############################################################################
##
## Any type used for standard errors must define the following methods:
##
##############################################################################
materialize(df::AbstractDataFrame, v::CovarianceEstimator) = v
completecases(df::AbstractDataFrame, ::CovarianceEstimator) = trues(size(df, 1))
S_hat(x::RegressionModel, ::CovarianceEstimator) = error("S_hat not defined for this type")
df_FStat(x::RegressionModel, ::CovarianceEstimator, hasintercept::Bool) = dof_residual(x) - hasintercept

include("utils.jl")
include("vcovsimple.jl")
include("vcovrobust.jl")
include("vcovcluster.jl")

end


text:
"""
二维模型的抽象类型
"""

module Fermi

abstract type Abstract2DModel{T} end

using ..Basics

export QuadrateSystem, TriangularSystem

"""
正方晶系
"""
struct QuadrateSystem{T} <: Abstract2DModel{T}
    brillouin :: Basics.AbstractRectangle{:AXISSQUARE}
    dispersion :: Vector{Function}
    bandnum :: Int64
    kadd :: Function
    QuadrateSystem{T}(disps) where T = begin
        new{T}(
            Square(Point2D(0., 0.), 2pi),
            disps,
            length(disps),
            square_kadd
        )
    end
end

"""
正方晶系的动量平移
"""
function square_kadd(pt1::Point2D, pt2::Point2D)
    destx = pt1.x + pt2.x
    desty = pt1.y + pt2.y
    #如果大于pi，就像0靠拢
    while abs(destx) > pi
        destx = destx - sign(destx) * 2pi
    end
    while abs(desty) > pi
        desty = desty - sign(desty) * 2pi
    end
    return Point2D(destx, desty)
end

"""
三角晶系，这里三角晶系的的初基矢量的选择
a1=(-1/2, √3/2); a2=(1/2, √3/2);
b1=π(-2, 2/√3); b2=π(2, 2/√3);
"""
struct TriangularSystem{T} <: Abstract2DModel{T}
    brillouin :: Basics.AbstractHexagon{:EQ}
    dispersion :: Vector{Function}
    bandnum :: Int64
    kadd :: Function
    TriangularSystem{T}(disps) where T = begin
        new{T}(
            EqHexagon(Point2D(0., 0.), 4pi/3.0),
            disps,
            length(disps),
            hexagon_kadd2
        )
    end
end

"""
三角晶系的动量平移
"""
function hexagon_kadd(pt1::Point2D, pt2::Point2D)
    @warn "hexagon_kadd很慢，用hexagon_kadd2代替"
    #找到离目标最近的一个第一布里渊区的中心
    #如果这个点距离中心是最近的，那么他就在第一布里渊区里面了，如果
    #离其他的点近，就平移一次，再检查一遍
    b1x = -2pi
    b1y = 2pi/sqrt(3)
    brlu_cnts = Vector{Point2D}(undef, 7)
    brlu_cnts[1] = Point2D(0., 0.)
    brlu_cnts[2] = Point2D(b1x, b1y)
    brlu_cnts[3] = Point2D(0., 2*b1y)
    brlu_cnts[4] = Point2D(-b1x, b1y)
    brlu_cnts[5] = Point2D(-b1x, -b1y)
    brlu_cnts[6] = Point2D(0., -2*b1y)
    brlu_cnts[7] = Point2D(b1x, -b1y)
    #
    dest = pt1 + pt2
    dis2cents = [dest - cnt for cnt in brlu_cnts]
    dis2cents = [pt.x^2 + pt.y^2 for pt in dis2cents]
    #
    (minv, mina) = findmin(dis2cents)
    while mina != 1
        dest = dest - brlu_cnts[mina]
        dis2cents = [dest - cnt for cnt in brlu_cnts]
        dis2cents = [pt.x^2 + pt.y^2 for pt in dis2cents]
        (minv, mina) = findmin(dis2cents)
    end
    return dest
end

"""
三角晶系的动量平移
b1=π(-2, 2/√3); b2=π(2, 2/√3);
"""
function hexagon_kadd2(pt1::Point2D, pt2::Point2D)
    dest = pt1 + pt2
    destx, desty = dest.x, dest.y
    #用一个长方型的框包起来
    while abs(desty) > π*2/√3
        desty = desty - sign(desty) * π*4/√3
    end
    while abs(destx) > 2*π
        destx = destx - sign(destx) * π*4
    end
    #如果在竖线上，肯定已经在第一布里渊区
    if isapprox(destx, 0.)
        return Point2D(destx, desty)
    end
    #如果在横线上
    if isapprox(desty, 0.)
        if destx > pi*4/3
            return Point2D(destx-π*2, π*2/√3)
        end
        if desty < -pi*4/3
            return Point2D(destx+π*2, π*2/√3)
        end
        return Point2D(destx, desty)
    end
    #如果在四个角落上
    lslope = -√3 * sign(destx) * sign(desty)
    y0 = sign(desty)*π*4/√3
    #println(lslope)
    #println(y0)
    #println(lslope*destx + y0)
    if desty*sign(desty) > (lslope*destx + y0)*sign(desty)
        destx = destx - sign(destx)*π*2
        desty = desty - sign(desty)*π*2/√3
    end
    return Point2D(destx, desty)
end

##"""不比kadd2要快
##三角晶系的动量平移
##b1=π(-2, 2/√3); b2=π(2, 2/√3);
##"""
##function hexagon_kadd3(pt1::Point2D, pt2::Point2D)
##    dest = pt1 + pt2
##    destx, desty = dest.x, dest.y
##    #用一个长方型的框包起来
##    while abs(desty) > π*2/√3
##        desty = desty - sign(desty) * π*4/√3
##    end
##    while abs(destx) > 2*π
##        destx = destx - sign(destx) * π*4
##    end
##    #如果在竖线上，肯定已经在第一布里渊区
##    if isapprox(destx, 0.)
##        return Point2D(destx, desty)
##    end
##    #
##    if sign(destx) * sign(desty) > 0
##        inp = 0.5*(√3)*destx + 0.5*desty
##        if inp > π*2/√3
##            return Point2D(destx - π*2, desty - π*2/√3)
##        elseif inp < -π*2/√3
##            return Point2D(destx + π*2, desty + π*2/√3)
##        end
##        return Point2D(destx, desty)
##    else
##        inp = -0.5*(√3)*destx + 0.5*desty
##        if inp > π*2/√3
##            return Point2D(destx + π*2, desty - π*2/√3)
##        elseif inp < -π*2/√3
##            return Point2D(destx - π*2, desty + π*2/√3)
##        end
##        return Point2D(destx, desty)
##    end
##end

#
#一些内建的模型
#不能在include上面加"""document"""

include("quadrate_lattice.jl")
export common_square_lattice

include("triangular_lattice.jl")
export common_triangle_lattice

include("kagome_lattice.jl")

include("surface.jl")

include("patch.jl")

end


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

module SIS3316Digitizers

using Base.Threads

using ArgCheck
using ArraysOfArrays
using BitOperations
using ElasticArrays
using Observables
using ParallelProcessingTools
using UnsafeArrays
using Sockets
using Tables
using TypedTables

using ..MemRegisters
using ..VMEGateways

include("evtformat.jl")
include("read_data.jl")
include("memory.jl")
include("memregs.jl")
include("memfifo.jl")
include("sis3316_digitizer.jl")

end # module


text:
#############################################################################
#############################################################################
#
# This file defines the polynomial type with several operations
#
#############################################################################
#############################################################################

####################################
# Polynomial and PolynomialModP Type #
####################################

    """
    A Polynomial type - designed to be for polynomials with integer coefficients.
    """
    struct Polynomial
        terms::Array{Term}
        Polynomial() = new(Term[])
        #Inner constructor
        Polynomial(h::Array{Term}, clean::Bool) = clean && new(h)
    end

    """
    Construct a polynomial with a single term.
    """
    function Polynomial(t::Term)
        terms = Term[]
        t.coeff != 0 && push!(terms, t)
        return Polynomial(terms,true)
    end

    """
    Construct a polynomial with a vector of terms.
    """
    function Polynomial(tv::Vector{Term})
        n = length(tv)
       sort!(tv, by = x -> -x.degree) #sorts
        if n == 1
            return Polynomial(tv, true)
        else
            terms = Term[]
            i = 1
            while i < n+1 #removes like terms
                if tv[i].coeff != 0
                    if i == n
                        push!(terms,tv[i])

                    elseif tv[i].degree == tv[i+1].degree
                        tv[i+1] = tv[i]+tv[i+1]
                    else
                        push!(terms,tv[i])
                    end
                end
                i = i + 1
            end
        end

        if length(terms)== 0 || terms[1].coeff == 0
            return Polynomial()
        end

        return Polynomial(terms, true)
        #Once made pretty (no repreated factors and ordered), the vector is sent to the inner constructor to
        #made into a polynomial
    end

    """
    Polymonials mod P
    """
    struct PolynomialModP
        polynomial::Polynomial
        mod::Int
        PolynomialModP(p) = new(Polynomial(Term[]),p)
        PolynomialModP(poly::Polynomial, mod::Int, modded::Bool) = modded && new(poly,mod) #see 109 for the explination of the bool
    end

    function PolynomialModP(t::Term,p)
        terms = Term[]
        modcoeff = mod(t.coeff,p)
        modcoeff != 0 && push!(terms,Term(modcoeff,t.degree))
        return PolynomialModP(Polynomial(terms),p,true)
    end

    function PolynomialModP(tv::Vector{Term},p::Int)
        n = length(tv)
        sort!(tv, by = x -> -x.degree) #it has to be sorted in order, soo this does that
        tv = mod.(tv,p)
        if n == 1
            terms = tv
        else
            terms = Term[]
            i = 1
            while i < n+1 # This adds like terms together
                if tv[i].coeff != 0
                    tv[i].coeff
                    if i == n
                        push!(terms,tv[i])
                    elseif tv[i].degree == tv[i+1].degree
                        tv[i+1] = tv[i]+tv[i+1]
                        tv[i+1].coeff = tv[i+1].coeff%p
                    else
                        push!(terms,tv[i])
                    end
                end
                i = i + 1
            end
        end

        if length(terms)== 0 || terms[1].coeff == 0
            return PolynomialModP(p)
        end

        terms = Polynomial(terms,true) #The true signifies that polynomial is odered and has no repeats,
        #as such can just be turned into a poly straight away
        return PolynomialModP(terms, p, true)
    end

    PolynomialModP(poly::Polynomial, p::Int) = PolynomialModP(mod(poly,p), p, true)

####################################
# Polynomial Construction #
####################################
    """
    Construct a polynomial of the form x^p-x.
    """
        cyclotonic_polynomial(p::Int) = Polynomial([Term(1,p), Term(-1,0)])

    """
    Construct a polynomial of the form x-n.
    """
        linear_monic_polynomial(n::Int) = Polynomial([Term(1,1), Term(-n,0)])

    """
    Construct a polynomial of the form x.
    """
    x_poly() = Polynomial(Term(1,1))
    x_poly(p::Int) = PolynomialModP(Term(1,1),p)

    """
    Creates the zero polynomial.
    """
    zero(::Type{Polynomial})::Polynomial = Polynomial()
    zero(::Type{PolynomialModP},n::Int)::PolynomialModP = PolynomialModP(n)
    """
    Creates the unit polynomial.
    """
    one(::Type{Polynomial})::Polynomial = Polynomial(Term(1,0))
    one(p::Polynomial) = one(typeof(p))

    one(::Type{PolynomialModP},p)::PolynomialModP = PolynomialModP(one(Term),p)
    one(p::PolynomialModP) = one(typeof(p),p.mod)

    """
    Generates a random polynomial.
    """
    function rand(::Type{Polynomial} ;
                    degree::Int = -1,
                    terms::Int = -1,
                    max_coeff::Int = 100,
                    mean_degree::Float64 = 5.0,
                    prob_term::Float64  = 0.7,
                    monic = false,
                    condition = (p)->true)

        while true
            _degree = degree == -1 ? rand(Poisson(mean_degree)) : degree
            _terms = terms == -1 ? rand(Binomial(_degree,prob_term)) : terms
            degrees = vcat(sort(sample(0:_degree-1,_terms,replace = false)),_degree)
            coeffs = rand(1:max_coeff,_terms+1)
            monic && (coeffs[end] = 1)
            p = Polynomial( [Term(coeffs[i],degrees[i]) for i in 1:length(degrees)] )
            condition(p) && return p
        end
    end

    function rand(::Type{PolynomialModP}, mod::Int ;
        degree::Int = -1,
        terms::Int = -1,
        mean_degree::Float64 = 5.0,
        prob_term::Float64  = 0.7,
        monic = false,
        condition = (p)->true)

        while true
            max_coeff = mod
            _degree = degree == -1 ? rand(Poisson(mean_degree)) : degree
            _terms = terms == -1 ? rand(Binomial(_degree,prob_term)) : terms
            degrees = vcat(sort(sample(0:_degree-1,_terms,replace = false)),_degree)
            coeffs = rand(1:max_coeff-1,_terms+1)
            monic && (coeffs[end] = 1)
            p = [Term(coeffs[i],degrees[i]) for i in 1:length(degrees)]
            poly = PolynomialModP(p,mod)
            condition(poly) && return poly
        end
    end

###########
# Display #
###########

    """
    Show a polynomial.
    """
    function show(io::IO, p::Polynomial)
        p = deepcopy(p)

        if iszero(p)
            print(io,"0")
        else
            n = length(p.terms)
            for (i,t) in enumerate(p.terms)

                if t.coeff < 0
                    print(io,i != 1 ? " - " : "-")
                    print(io,-t)
                else
                    print(io,i != 1 ? " + " : "")
                    print(io,t)
                end
            end
        end
    end

    show(io::IO,poly::PolynomialModP) = print(io, poly.polynomial, " (mod ", poly.mod, ")")

##############################################
# Iteration over the terms of the polynomial #
##############################################

    """
    Allows to do iteration over the terms of the polynomial. The iteration is in an arbitrary order.
    """
    iterate(p::Polynomial, state=1) = iterate(p.terms, state)
    iterate(p::PolynomialModP, state=1) = iterate(p.polynomial.terms, state)

    function coeffvector(p::Polynomial, size::Int = -1)::Vector{Int}
        size == -1 && (size = degree(p)+1)
        a = zeros(size)
        if length(p.terms) != size
            for t in p.terms
                a[size - t.degree] = t.coeff
            end
        else
            a = coeffs(p)
        end
        return a
    end

    coeffvector(p::PolynomialModP,size::Int = -1) = coeffvector(p.polynomial, size)

##############################
# Queries about a polynomial #
##############################

    """
    The number of terms of the polynomial.
    """
    length(p::Polynomial) = length(p.terms)
    length(p::PolynomialModP) = length(p.polynomial)

    """
    The leading term of the polynomial.
    """
    leading(p::Polynomial)::Term = isempty(p.terms) ? zero(Term) : first(p.terms)
    leading(p::PolynomialModP)::Term = leading(p.polynomial)

    """
    Returns the coefficients of the polynomial.
    """
    coeffs(p::Polynomial)::Vector{Int} = [t.coeff for t in p]
    coeffs(p::PolynomialModP)::Vector{Int} = coeffs(p.polynomial)

    """
    The degree of the polynomial.
    """
    degree(p::Polynomial)::Int = leading(p).degree
    degree(p::PolynomialModP)::Int = leading(p.polynomial).degree

    """
    The content of the polynomial is the GCD of its coefficients.
    """
    content(p::Polynomial)::Int = euclid_alg(coeffs(p))
    content(p::PolynomialModP)::Int = content(p.polynomial)

    """
    Evaluate the polynomial at a point `x`.
    """
    evaluate(f::Polynomial, x::T) where T <: Number = sum(evaluate(t,x) for t in f)
    evaluate(f::PolynomialModP, x::T) where T <: Number = evaluate(f.polynomial,x)%f.mod

################################
# Pushing and popping of terms #
################################

    """
    Push a new term into the polynomial.
    """
    #Note that ideally this would throw and error if pushing another term of degree that is already in the polynomial
    function push!(p::Polynomial, t::Term)
        iszero(t) && return #don't push a zero
        Polynomial(push!(p.terms,t))
    end

    function push!(p::PolynomialModP, t::Term)
        iszero(t) && return #don't push a zero
        PolynomialModP(push!(p.polynomial,t),p.mod)
    end

    """
    Pop the leading term out of the polynomial.
    """
    pop!(p::Polynomial)::Term = pop!(p.terms)
    pop!(p::PolynomialModP)::Term = pop!(p.polynomial.terms)

    prepop!(p::Polynomial)::Term = popat!(p.terms,1) #I needed a way to "pop" the first term (When changing to array from heap)
    prepop!(p::PolynomialModP)::Term = popat!(p.polynomial.terms,1)
    """
    Check if the polynomial is zero.
    """
    iszero(p::Polynomial)::Bool = isempty(p.terms)
    iszero(p::PolynomialModP)::Bool = iszero(p.polynomial)

#################################################################
# Transformation of the polynomial to create another polynomial #
#################################################################

    """
    The negative of a polynomial.
    """
    -(p::Polynomial) = Polynomial(map((pt)->-pt, p.terms))
    -(p::PolynomialModP) = PolynomialModP(-p.polynomial,p.mod)

    """
    Create a new polynomial which is the derivative of the polynomial.
    """
    function derivative(p::Polynomial)::Polynomial
        der_p = Polynomial()
        for term in p
            der_term = derivative(term)
            !iszero(der_term) && push!(der_p,der_term)
        end
        return der_p
    end

    derivative(p::PolynomialModP)::PolynomialModP = PolynomialModP(derivative(p.polynomial),p.mod)

    """
    The prim part (multiply a polynomial by the inverse of its content).
    """
    prim_part(p::Polynomial) = p ÷ content(p)
    prim_part(p::PolynomialModP) = p ÷ content(p)

    """
    A square free polynomial.
    """

    square_free(p::PolynomialModP)::PolynomialModP = p ÷ gcd(p,derivative(p))

#################################
# Queries about two polynomials #
#################################

    """
    Check if two polynomials are the same
    """
    ==(p1::Polynomial, p2::Polynomial)::Bool = p1.terms == p2.terms
    ==(p1::PolynomialModP, p2::PolynomialModP)::Bool = p1.polynomial.terms == p2.polynomial.terms

    """
    Check if a polynomial is equal to 0.
    """
    #Note that in principle there is a problem here. E.g The polynomial 3 will return true to equalling the integer 2.
    ==(p::Polynomial, n::T) where T <: Real = iszero(p) == iszero(n)
    ==(p::PolynomialModP, n::T) where T <: Real = ==(p.polynomial,n)

##################################################################
# Operations with two objects where at least one is a polynomial #
##################################################################
    """
    Subtraction of two polynomials.
    """
    -(p1::Polynomial, p2::Polynomial)::Polynomial = p1 + (-p2)
    -(p1::PolynomialModP, p2::PolynomialModP)::PolynomialModP = p1 + (-p2)

    -(p::PolynomialModP, n::Int) = p + PolynomialModP(Term(-n,0),p.mod)
    #Subtract an integer from a Polynomial
    -(p::Polynomial, n::Int)::Polynomial = p + Polynomial(Term(-n,0))

    """
    Multiplication of polynomial and term.
    """
    *(t::Term,p1::Polynomial)::Polynomial = iszero(t) ? Polynomial() : Polynomial(map((pt)->t*pt, p1.terms))
    *(p1::Polynomial, t::Term)::Polynomial = t*p1

    *(t::Term,p1::PolynomialModP) = iszero(t) ? PolynomialModP(p1.mod) : PolynomialModP(map((pt)->mod(t*pt,p1.mod), p1.polynomial.terms),p1.mod)
    *(p1::PolynomialModP, t::Term) = t*p1

    """
    Multiplication of polynomial and an integer.
    """
    *(n::Int,p::Polynomial)::Polynomial = p*Term(n,0)
    *(p::Polynomial,n::Int)::Polynomial = n*p

    *(n::Int,p::PolynomialModP) = p*Term(n,0)
    *(p::PolynomialModP,n::Int) = n*p

    """
    Integer division of a polynomial by an integer.

    Warning this may not make sense if n does not divide all the coefficients of p.
    """
    ÷(p::Polynomial,n::Int) = (prime)->Polynomial(map((pt)->((pt ÷ n)(prime)), p.terms))
    ÷(p::PolynomialModP, n::Int) = PolynomialModP(÷(p.polynomial,n)(p.mod),p.mod)

    """
    Take the s of a polynomial with an integer.
    """
    function mod(f::Polynomial, p::Int)::Polynomial
        p_out = Polynomial()
        for t in f
            push!(p_out, mod(t, p)) #if coeff reduced to zero, push! will handle it
        end
        return p_out
    end

" ---- Polynomial Run ----"


text:
using Distributions
using Distances
include("./Spelliptic.jl")
using Spelliptic
include("./CovFuncs.jl")

#setting the random seed
srand(123)

#Setting problem parameters
d = 2
N = 50000
l = 1.0
nu = 0.5
rho = 4.

#defining covariance function
cov = r -> matern( r, l, nu )

#generating point cloud
dUniform = Uniform(-3,3)
x = rand( dUniform, d, N )

println( "Test without nugget, using SpelMat:" )
#Generating a SpelMat structure
@time K = SpelMat( x , rho )
#Setting the content of the spelliptic struct according to the covariance
#matrix and point cloud:
@time setCovStat!( K, x, cov )

#Measuring elementwise approximation error
trueSol, appSol, error, distance = testAccuracy( K, x, cov, 500 )
println( "mean error = ", mean( error ), " median error = ", median( error ) )

vec = randn( N )
#Matrix vector product:
K * vec;
#Solution of linear system
K \ vec
#computing the logdeterminant
println( "logdet(K) = ", logdet( K ))

println( "Test with large nugget, using SpelMatIter:")
sigma = 10.
@time KIter = SpelMatIter( x , rho )
@time setCovStat!( KIter, x, cov, sigma )

#Matrix vector product:
KIter * vec;
#Solution of linear system
KIter \ vec

max_stars_repo_name: f-t-s/Spelliptic

max_issues_repo_name: f-t-s/Spelliptic

max_forks_repo_name: f-t-s/Spelliptic


text:
using SimpleDifferentialOperators
using Test, LinearAlgebra, DualNumbers

@testset "Accuracy & regression test" begin
    regular_grid = 0:3
    irregular_grid = [0.0; 1.0; 2.0; 4.0; 6.0]
    ## regular grids
    L₁₋_regular = L₁₋(regular_grid)
    L̄₁₊_regular = L₁₊(regular_grid)
    L̄₂_regular = L₂(regular_grid)
    x_regular = interiornodes(regular_grid)
    @test Array(L₁₋_regular) == [-1. 1. 0. 0.; 0. -1. 1. 0.]
    @test Array(L̄₁₊_regular) == [0. -1. 1. 0.; 0. 0. -1. 1.]
    @test Array(L̄₂_regular) == [1. -2. 1. 0.; 0. 1. -2. 1.]
    @test Array(x_regular) == [1; 2]

    @test @inferred(L₁₋(regular_grid)) == L₁₋_regular
    @test @inferred(L₁₊(regular_grid)) == L̄₁₊_regular
    @test @inferred(L₂(regular_grid)) == L̄₂_regular
    @test @inferred(interiornodes(regular_grid)) == x_regular

    ## irregular grids
    L₁₋_irregular = L₁₋(irregular_grid)
    L̄₁₊_irregular = L₁₊(irregular_grid)
    L̄₂_irregular = L₂(irregular_grid)
    x_irregular = interiornodes(irregular_grid)
    @test Array(L₁₋_irregular) == [-1. 1. 0. 0. 0.; 0. -1. 1. 0. 0.; 0. 0. -1/2 1/2 0.]
    @test Array(L̄₁₊_irregular) == [0. -1. 1. 0. 0.; 0. 0. -1/2 1/2 0.; 0. 0. 0. -1/2 1/2]
    @test Array(L̄₂_irregular) == [1. -2. 1. 0. 0.; 0 2/((1+2)*1) -2/(2*1) 2/((1+2)*2) 0; 0. 0. 1/4 -2/(2*2) 1/4]
    @test Array(x_irregular) == [1; 2; 4]

    @test @inferred(L₁₋(irregular_grid)) == L₁₋_irregular
    @test @inferred(L₁₊(irregular_grid)) == L̄₁₊_irregular
    @test @inferred(L₂(irregular_grid)) == L̄₂_irregular
    @test @inferred(interiornodes(irregular_grid)) == x_irregular
end

@testset "Consistency" begin
    regular_grids = [range(0.0, 1.0, length = 5), range(-1.0, 1.0, length = 5)]
    irregular_grids = [[-1.0; 0.0; 2.0; 5.0; 9.0], [0.1; 0.3; 0.9; 1.9; 2.5; 9.0]]
    # Check consistency
    for x̄ in [regular_grids; irregular_grids]
        # setup
        f(x) = x^2
        μ = -0.1 # constant negative drift
        σ = 0.1
        ρ = 0.05

        x = interiornodes(x̄)
        M = length(x)
        bc = (Reflecting(), Reflecting()) # specify BC (reflecting barrier)

        # operators with reflecting boundary conditions
        L = μ*L₁₋bc(x̄, bc) + σ^2 / 2 * L₂bc(x̄, bc)
        ## solve the value function
        v = (I * ρ - L) \ f.(x)

        # operators without boundary conditions, adding extra two rows for boundary conditions
        L̄ = μ*L₁₋(x̄) + σ^2 / 2 * L₂(x̄)
        B = transpose([[-1; 1; zeros(M)] [zeros(M); -1; 1]])
        L = [([zeros(M) Diagonal(ones(M,M)) zeros(M)] * ρ - L̄); B]

        ## solve the value function including the boundary conditions
        v̄ = L \ [f.(x); 0.0; 0.0]

        # currently broken on regular grids
        @test v ≈ v̄[2:end-1]
    end
end


text:
# this is series of tests for all methods 'get_right_endpoint'

for method = methods_list
  str = "Using " * String(method)
  @testset "$str" begin
      @testset "f_2p_1im" begin
          res0 = [get_right_endpoint(
              [3., 1.],
              i,
              (x::Vector{Float64}) -> f_2p_1im(x) - 9.,
              Val(method)
          ) for i in 1:2]
          @test isapprox(res0[1][1], 5., atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[2][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_2p_1im max_iter=5" begin
          res0 = [get_right_endpoint(
              [3., 1.],
              i,
              (x::Vector{Float64}) -> f_2p_1im(x) - 9.,
              Val(method),
              max_iter = 5
          ) for i in 1:2]
          @test res0[1][1] == nothing
          @test res0[1][3] == :MAX_ITER_STOP
          @test res0[2][1] == nothing
          @test res0[2][3] == :MAX_ITER_STOP
      end

      @testset "f_2p" begin
          res0 = [get_right_endpoint(
              [3., 1.1],
              i,
              (x::Vector{Float64}) -> f_2p(x) - 9.,
              Val(method);
              scan_tol = 1e-6
          ) for i in 1:2]
          @test isapprox(res0[1][1], 5., atol=1e-6)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test isapprox(res0[2][1], 6., atol=1e-6)
          @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
      end

      @testset "f_3p_1im" begin
          res0 = [get_right_endpoint(
              [3., 8., 2.1],
              i,
              (x::Vector{Float64}) -> f_3p_1im(x) - 9.,
              Val(method)
          ) for i in 1:3]
          @test isapprox(res0[1][1], 5., atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[2][3] == :SCAN_BOUND_REACHED
          @test res0[3][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_3p_1im restricted" begin
          res0 = [get_right_endpoint(
              [3., 8., 2.1],
              i,
              (x::Vector{Float64}) -> f_3p_1im(x) - 9.,
              Val(method);
              scan_bound = 4.
          ) for i in 1:3]
          @test res0[1][3] == :SCAN_BOUND_REACHED
          @test res0[2][3] == :SCAN_BOUND_REACHED
          @test res0[3][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_3p_1im_dep" begin
          res0 = [get_right_endpoint(
              [3., 2., 2.1],
              i,
              (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
              Val(method)
          ) for i in 1:3]
          @test isapprox(res0[1][1], 5., atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2)
          @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[3][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_3p_1im_dep scan_bound" begin
          res0 = [get_right_endpoint(
              [3., 2., 2.1],
              i,
              (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
              Val(method);
              scan_bound = [4.,10.,10.][i]
          ) for i in 1:3]
          @test res0[1][3] == :SCAN_BOUND_REACHED
          @test isapprox(res0[2][1], 2.0+2.0*sqrt(2.), atol=1e-2)
          @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[3][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_3p_1im_dep scan_tol" begin
          res0 = [get_right_endpoint(
              [3., 2., 2.1],
              1,
              (x::Vector{Float64}) -> f_3p_1im_dep(x) - 9.,
              Val(method);
              scan_tol = [1e-2, 1e-4, 1e-6][i]
          ) for i in 1:3]
          @test isapprox(res0[1][1], 5., atol=1e-1)
          @test isapprox(res0[2][1], 5., atol=1e-3)
          @test isapprox(res0[3][1], 5., atol=1e-5)
      end

      @testset "f_4p_2im" begin
          res0 = [get_right_endpoint(
              [3., 4, 1.1, 10.],
              i,
              (x::Vector{Float64}) -> f_4p_2im(x) - 9.,
              Val(method)
          ) for i in 1:4]
          @test isapprox(res0[1][1], 5., atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test isapprox(res0[2][1], 6., atol=1e-2)
          @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[3][3] == :SCAN_BOUND_REACHED
          @test res0[4][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_4p_3im" begin
          res0 = [get_right_endpoint(
              [3., 4, 1.1, 10.],
              i,
              (x::Vector{Float64}) -> f_4p_3im(x) - 9.,
              Val(method)
          ) for i in 1:4]

          @test isapprox(res0[1][1], 5., atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[2][3] == :SCAN_BOUND_REACHED
          @test res0[3][3] == :SCAN_BOUND_REACHED
          @test res0[4][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_1p_ex" begin
          res0 = [get_right_endpoint(
              [1.5],
              i,
              (x::Vector{Float64}) -> f_1p_ex(x) - 9.,
              Val(method)
          ) for i in 1:1]

          @test isapprox(res0[1][1], 2. + 1e-8, atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
      end

      @testset "f_5p_3im" begin
          res0 = [get_right_endpoint(
              [3., 0.1, 4, 1.1, 8.],
              i,
              (x::Vector{Float64}) -> f_5p_3im(x) - 9.,
              Val(method)
          ) for i in 1:5]

          @test isapprox(res0[1][1], 5., atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test isapprox(res0[2][1], log(3.), atol=1e-2)
          @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[3][3] == :SCAN_BOUND_REACHED
          @test res0[4][3] == :SCAN_BOUND_REACHED
          @test res0[5][3] == :SCAN_BOUND_REACHED
      end

      @testset "f_3p_im" begin
          res0 = [get_right_endpoint(
              [3., 0.1, 8.],
              i,
              (x::Vector{Float64}) -> f_3p_im(x) - 9.,
              Val(method)
          ) for i in 1:3]

          @test isapprox(res0[1][1], 5., atol=1e-2)
          @test res0[1][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test isapprox(res0[2][1], log(3.), atol=1e-2)
          @test res0[2][3] == :BORDER_FOUND_BY_SCAN_TOL
          @test res0[3][3] == :SCAN_BOUND_REACHED
      end
  end
end


text:
#!/usr/bin/env julia

using Keldysh, Test, LinearAlgebra

include("storage.jl")
include("branch.jl")
include("contour.jl")
include("time_grid.jl")
include("dos.jl")
include("time_gf.jl")


text:
import StaticArrays: SVector

@testset "Acoustic constructors" begin
    # 2D Acoustic
    a2 = Acoustic(0.1, 0.1+0.0im, 2)
    @test spatial_dimension(a2) == 2
    @test field_dimension(a2) == 1

    # 3D Acoustic
    a3 = Acoustic(1.0, 1.0+0.0im, 3)
    @test spatial_dimension(a3) == 3
    @test field_dimension(a3) == 1

    # Boundary condition constuctors
    @test sound_hard(2) == hard(a2)
    @test sound_hard(2) == rigid(a2)
    @test sound_hard(2) == zero_neumann(a2)

    @test sound_soft(2) == soft(a2)
    @test sound_soft(2) == pressure_release(a2)
    @test sound_soft(2) == zero_dirichlet(a2)

    @test sound_hard(2)  != sound_soft(2)
end

@testset "Acoustic special functions" begin
    T = Float64
    ms = 0:4
    map(2:3) do Dim
        L1 = basisorder_to_basislength(Acoustic{T,Dim},ms[1])
        L2 = basisorder_to_basislength(Acoustic{T,Dim},ms[end])

        m = basislength_to_basisorder(Acoustic{T,Dim}, L2)

        @test m == ms[end]
        @test L1 == 1
    end

    # Test 3D outgoing translation matrix
    ω = rand() + 0.1
    medium = Acoustic(3; ρ = 1.0, c = 1.0)
    r = rand(3) - [0.5,0.5,0.5];
    d = rand(3) - [0.5,0.5,0.5];
    d = 10 * d * norm(r) / norm(d)

    order = 2
    U = outgoing_translation_matrix(medium, 4*order, ω, d)
    vs = regular_basis_function(medium, ω)(4*order,r)
    us = outgoing_basis_function(medium, ω)(order,r + d)

    L = length(us)
    @test maximum(abs.( (U * vs)[1:L] - us) ./ abs.(us)) < 3e-6 # for linux: 4e-7

    # Test 2D outgoing translation matrix
    ω = rand() + 0.1
    medium = Acoustic(2; ρ = 1.0, c = 1.0)
    r = rand(2) - [0.5,0.5];
    d = rand(2) - [0.5,0.5];
    d = 10 * d * norm(r) / norm(d)

    # Note that to be accurate the order of vs
    order = 4
    larger_order = 3order
    U = outgoing_translation_matrix(medium, larger_order, ω, d)
    vs = regular_basis_function(medium, ω)(larger_order,r)
    us = outgoing_basis_function(medium, ω)(order,r + d)

    L = larger_order+1
    @test maximum(abs.((U * vs)[L-order:L+order] - us) ./ abs.(us)) < 1e-9
end

@testset "Acoustic circle T-matrix" begin

    circle = Sphere((0.0,0.0), 2.0)
    a2 = Acoustic(0.1, 0.1+0.0im, 2)
    a2_host = Acoustic(1.0, 1.0+0.0im, 2)
    ω = 0.5

    N_basis = estimate_outgoing_basisorder(a2_host, Particle(a2,circle), ω)

    # Test return type is satisfied for valid input
    t = t_matrix(Particle(a2,circle), a2_host, ω, N_basis)
    @test typeof(t) <: Diagonal{Complex{Float64}}

    p = Particle(Acoustic(Inf, 0.0im, 2),circle)
    @test_throws DomainError t_matrix(p, Acoustic(1.0, 1.0+0.0im, 2), ω, N_basis)
    p = Particle(Acoustic(1.0, 1.0+0.0im, 2),circle)
    @test_throws DomainError t_matrix(p, Acoustic(0.0, Inf*im, 2), ω, N_basis)
    @test_throws DomainError t_matrix(p, Acoustic(1.0, 0.0im, 2), ω, N_basis)
    p = Particle(Acoustic(0.0, 1.0im, 2),circle)
    @test_throws DomainError t_matrix(p, Acoustic(0.0, 1.0+0.0im, 2), ω, N_basis)
    p = Particle(a2, Sphere(SVector(1.0, 1.0), 0.0))
    @test_throws DomainError t_matrix(p, a2_host, ω, N_basis)

end


text:
using Klara

plogtarget(p::Int, v::Vector) = p*log(v[1])-log(factorial(p))

p = BasicDiscUnvParameter(:p, logtarget=plogtarget, nkeys=2)

model = likelihood_model([Constant(:λ), p], isindexed=false)

psetproposal(i::Int) = (i == 0) ? Binary(0, 1) : Binary(i-1, i+1)

sampler = MH(psetproposal, symmetric=false)

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:λ=>6, :p=>2)

job = BasicMCJob(model, sampler, mcrange, v0)

run(job)

chain = output(job)

mean(chain)

acceptance(chain, diagnostics=false)


text:
using SymPy

function ry_mat(ang)
    return [SymPy.cos(ang/2) -SymPy.sin(ang/2); SymPy.sin(ang/2) SymPy.cos(ang/2)]
end

function rz_mat(ang)
    return [SymPy.exp(-1im*ang/2) 0; 0 SymPy.exp(1im*ang/2)]
end

function gen_mat(a,b,c,d)
    return SymPy.exp(1im*d).*ry_mat(c)*rz_mat(b)*ry_mat(a)
end

function _params_zyz(mat::Matrix{<:Symbol})
    """Return the euler angles and phase for the ZYZ basis."""
    # We rescale the input matrix to be special unitary (det(U) = 1)
    # This ensures that the quaternion representation is real
    coeff = sympy.sqrt(sympy.det(mat))
    phase = -sympy.atan2(real(coeff),imag(coeff))
    su_mat = coeff * mat  # U in SU(2)
    # OpenQASM SU(2) parameterization:
    # U[0, 0] = exp(-i(phi+lambda)/2) * cos(theta/2)
    # U[0, 1] = -exp(-i(phi-lambda)/2) * sin(theta/2)
    # U[1, 0] = exp(i(phi-lambda)/2) * sin(theta/2)
    # U[1, 1] = exp(i(phi+lambda)/2) * cos(theta/2)
    theta = 2 * sympy.atan2(abs(su_mat[2, 1]), abs(su_mat[1, 1]))
    phiplambda = 2 * sympy.atan2(real(su_mat[2, 2]), imag(su_mat[2, 2]))
    phimlambda = 2 * sympy.atan2(real(su_mat[2, 1]), imag(su_mat[2, 1]))
    phi = (phiplambda + phimlambda) / 2
    lam = (phiplambda - phimlambda) / 2
    return theta, phi, lam, phase
end

function _params_zyz(mat::Matrix{<:Number})
    """Return the euler angles and phase for the ZYZ basis."""
    # We rescale the input matrix to be special unitary (det(U) = 1)
    # This ensures that the quaternion representation is real
    coeff = sqrt(det(mat) +0im)
    phase = -atan(real(coeff), imag(coeff))
    su_mat = coeff * mat  # U in SU(2)
    # OpenQASM SU(2) parameterization:
    # U[0, 0] = exp(-i(phi+lambda)/2) * cos(theta/2)
    # U[0, 1] = -exp(-i(phi-lambda)/2) * sin(theta/2)
    # U[1, 0] = exp(i(phi-lambda)/2) * sin(theta/2)
    # U[1, 1] = exp(i(phi+lambda)/2) * cos(theta/2)
    theta = 2 * atan(abs(su_mat[2, 1]), abs(su_mat[1, 1]))
    phiplambda = 2 * atan(real(su_mat[2, 2]), imag(su_mat[2, 2]))
    phimlambda = 2 * atan(real(su_mat[2, 1]), imag(su_mat[2, 1]))
    phi = (phiplambda + phimlambda) / 2
    lam = (phiplambda - phimlambda) / 2
    return theta, phi, lam, phase
end