text:
# defining constants
amu = 931.4940954E6   # eV / c^2
ħc = 0.19732697E-6    # eV m

export potential, pot_from_dict, fit_pot!, solve_pot!, find_crossing
export Plotter
# export solve1D_ev_amu
# export sqwell, harmonic, double_well, polyfunc, morse
# export get_bspline, get_spline

"""
Stores a potential in one-dimensional space Q, with discreet points (E0, Q0) and fitting function func.

## Fields

- `name` -- the name of potential.
- `QE_data`   -- the (n X 2) DataFrame of data points (Q vs Energy).
- `E0`, `Q0`  -- the minimum point of the potential [`Q0`, `E0`].
- `func_type`     -- the type of fitting function ("bspline", "spline", "harmonic", "polyfunc", "morse_poly", "morse").
- `params`    -- the list of hyper paramters for the fitting function.
- `Q`, `E`  -- `Q` and `E`=`func(Q, p_opt; params)`.
- `nev`  -- the number of eigenvalues to be evaluated.
- `ϵ` -- the list of eigenvalues

## Constructor
    potential()

"""
mutable struct potential
    name::String
    QE_data::DataFrame
    E0::Float64; Q0::Float64
    func_type::String
    func # can be either Function or Interpolations.ScaledInterpolation
    params::Dict{String, Any}
    Q::Array{Float64,1}; E::Array{Float64,1}
    nev::Int
    # ϵ includes E0
    ϵ::Array{Float64,1}; χ::Array{Float64,2}
    potential() = new("", DataFrame([0 0], [:Q, :E]), Inf, 0,
                      "func_type", x->0, Dict(),
                      [], [],
                      0, [], Array{Float64}(undef, 0, 2))
end

"""
    fit_pot!(pot::potential, Q; params = nothing)

Fit a function `pot.func_type` to `QE_data` on the domain `Q`.

## parameters

- `pot`: `potential`
    - `pot.func_type`: the fitting function; `{"spline" (preferred), "bspline", "harmonic", "polyfunc", "morse_poly", "morse"}`.
- `Q`: the spatial domain (1-d array).
- `params`: the hyperparameters.

### Hyperparameters `params`

- spline (preferred)

    Spline interpolation. See more detail in [Dierckx.jl](https://github.com/kbarbary/Dierckx.jl).
    - `order`: spline order (between 1 and 5; default 2).
    - `smoothness`: the amount of smoothness is determined by the condition that `sum((w[i]*(y[i]-spline(x[i])))**2) <= s`
    - `weight`: the weight applied to each `QE_data` point (length m 1-d array).

- bspline

    Basic spline interpolation. See more detail in [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl).

- harmonic

    Harmonic function whose minimum is at [`pot.Q0`, `pot.E0`].
    - `hw`: the energy quanta of the harmonic oscillator.

- polyfunc

    Polynomial function;
        `y = E₀ + Σ coeffs[i].* (x .- Q₀) .^(i-1)`.
    - `poly_order`: the maximum order of polynomials.
    - `p0`: the initial paramters for the fitting function.

## Example
- Spline fit
```julia
nev = 60
name = "DX-+h"

Q1 = [30.66721918 29.860133 29.05306268 28.24612992 27.43911044 26.63197583 25.82507425 25.01797262 24.21096115 23.40387798 22.59690043 21.78991117 20.98281579 20.17570172 19.36884219 18.56169249 17.75463179 16.94772679 16.14061031 15.33347439 14.52663309 13.71945696 12.91240658 12.10544441 11.29841968 10.4913338 9.684370388 8.877289725 8.070184138]
E1 = [8.0902 7.5970 7.0749 6.5242 5.9461 5.3451 4.7300 4.1147 3.5182 2.9637 2.4769 2.0819 1.7972 1.6315 1.5800 1.6237 1.7301 1.8586 1.9687 2.0283 2.0190 2.0673 1.9910 2.0528 1.9730 2.0857 2.4550 3.1653 4.3448]

pot = potential(); pot.name = name
pot.nev = nev
pot.Q0 = Q1[findmin(E1)[2]]; pot.E0 = 1.69834
pot.QE_data = DataFrame(Q = Q1[:], E = E1[:])
pot.QE_data.E .+= - minimum(pot.QE_data.E) + pot.E0
pot.Q = Q

pot.func_type = "spline"
# spline fitting parameters
params = Dict()
params["order"] = 4
params["smoothness"] = 0.001
params["weight"] = [1 1 1 1 1 0.5 0.4 0.4 0.5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
fit_pot!(pot, Q; params=params)
```

- Harmonic
```julia
nev = 40
name = "SnIII+h"

pot = potential(); pot.name = name
pot.Q0 = 1.90291674728; pot.E0 = 0.585005
pot.nev = nev
pot.func_type = "harmonic"
params = Dict()
params["hw"] = 0.0281812646475
fit_pot!(pot, Q; params = params)
```

"""
function fit_pot!(pot::potential, Q; params = nothing)
    # pot.params["E0"] = pot.E0
    # pot.params["Q0"] = pot.QE_data.Q[findmin(pot.QE_data.E)[2]]

    E_CUT = 2 # defines upper energy limit on data used for fitting
    e_cut_ind = pot.QE_data.E .< E_CUT + pot.E0 # boolean

    params = if params == nothing pot.params else params end
    pot.Q = Q

    println("Potential fitting: $(pot.name)")

    if pot.func_type == "bspline"
        println("=========bspline=========\n")

        bspline = get_bspline(pot.QE_data.Q, pot.QE_data.E)
        pot.E = bspline(Q)
        pot.func = bspline

    elseif pot.func_type == "spline"
        println("=========spline==========\n")

        weight = get(params, "weight", nothing)
        if isa(weight, String)
            weight = parse.(Float64, split(weight))
        end
        if weight ≠ nothing
            weight = vec(weight)
        end
        smoothness = get(params, "smoothness", 0)
        order = get(params, "order", 2)

        spline = get_spline(pot.QE_data.Q, pot.QE_data.E; weight = weight, smoothness = smoothness, order = order)
        pot.E = spline(Q)
        pot.func = spline

    elseif pot.func_type == "harmonic"
        println("========harmonic=========\n")
        # func = x -> harmonic(x, params["hw"]; param = params)
        pot.E = harmonic.(Q, params["hw"]; E₀ = pot.E0, Q₀ = pot.Q0)
        pot.func = x -> harmonic(x, params["hw"]; E₀ = pot.E0, Q₀ = pot.Q0)

    else
        if pot.func_type == "polyfunc"
            println("========polynomial========\n")
            func = (x, p) -> polyfunc(x, p; E₀ = pot.E0, Q₀ = pot.Q0, poly_order = params["poly_order"])
        elseif pot.func_type == "morse_poly"
            println("========morse polynomial========\n")
            func = (x, p) -> morse_poly(x, p; E₀ = pot.E0, Q₀ = pot.Q0, poly_order = params["poly_order"])
        elseif pot.func_type == "morse"
            func = (x, p) -> morse(x, p; E₀ = pot.E0, Q₀ = pot.Q0)
        end

        fit = curve_fit(func, pot.QE_data.Q[e_cut_ind], pot.QE_data.E[e_cut_ind], Float64.(params["p0"])) # curve_fit : model, x data, y data, p0
        pot.E = func.(Q, Ref(fit.param)) # when calling function, need to Ref() so that it can deal with array mismatch
        pot.func = x -> func(x, fit.param)

        println("===========Fit===========")
        println("Function: $(pot.func_type)")
        println("Best fit: $(fit.param)")
        println("=========================\n")
    end
end

"""
    solve_pot!(pot::potential)

Solve 1D Shrödinger equation for `potential`. `pot.ϵ` and `pot.χ` store the eigenvalues and eigenvectors, respectively.
"""
function solve_pot!(pot::potential)
    # pot.ϵ, pot.χ = solve1D_ev_amu(pot.func, pot.Q; nev=pot.nev)
    # Q = pot.Q
    # NQ=length(Q); Qi=minimum(Q); Qf=maximum(Q)
    pot.ϵ, pot.χ = solve1D_ev_amu(pot.func, pot.Q; nev=pot.nev)
end

"""
Solve 1D Shrödinger equation. The Brooglie wrapper in the unit of `eV` and `amu`.
"""
function solve1D_ev_amu(func, Q::AbstractArray{<:Real,1}; nev=30, maxiter=nev*length(Q))
    NQ=length(Q); Qi=minimum(Q); Qf=maximum(Q)
    return solve1D_ev_amu(func, NQ, Qi, Qf; nev=nev, maxiter=maxiter)
end

function solve1D_ev_amu(func, NQ::Int, Qi::Real, Qf::Real; nev=30, maxiter=nev*length(Q))
    factor = (1/amu) * (ħc*1E10)^2

    ϵ1, χ1 = Brooglie.solve(x -> func.(x*factor^0.5);
                  N=NQ, a=Qi/factor^0.5, b=Qf/factor^0.5, m=1, nev=nev, maxiter=maxiter)
    return ϵ1, vcat(χ1'...)/factor^0.25
end

"""
    find_crossing(pot_1::potential, pot_2::potential)

Find the crossing point between two potential energy surfaces `pot1` and `pot2`.
`Qx, Ex = find_crossing(pot1, pot2)`.
"""
function find_crossing(pot_1::potential, pot_2::potential)
    # find root of pot_1 - pot_2 = 0
    diff_func = x-> pot_1.func(x) - pot_2.func(x)
    Q = pot_1.Q
    rts = find_zero(diff_func, Q[length(Q)÷2])
    return rts, pot_1.func.(rts)
end

# read potential
"""
Depreciated.
Construct `potential` from `QE_data` and configure dictionalry `cfg`.
"""
function pot_from_dict(QE_data::DataFrame, cfg::Dict)::potential
    pot = potential()
    pot.name = cfg["name"]
    pot.nev = cfg["nev"]
    pot.func_type = cfg["function"]["type"]
    # pot.p0 =  parse.(Float64, split(get(cfg["function"], "p0", "1 1")))
    pot.E0 = get(cfg, "E0", Inf)
    pot.QE_data = QE_data

    pot.params = get(cfg["function"], "params",  Dict("E0" => pot.E0))
    pot.params["p0"] =  parse.(Float64, split(get(cfg["function"], "p0", "1 1")))

    pot.Q0 = pot.QE_data.Q[findmin(pot.QE_data.E)[2]]

    if pot.E0 < Inf
        pot.QE_data.E .+= - minimum(pot.QE_data.E) + pot.E0
    end

    return pot
end

############################################################
# Set of potentials
function sqwell(x, width, depth; x0=0.)
    x = x .- x0
    pot_well = (x .< -width/2.) .| (x .> width/2.)
    return pot_well .* depth
end

function harmonic(x, ħω; E₀, Q₀)
    a = amu/2 * (ħω/ħc/1E10)^2
    return a*(x - Q₀)^2 + E₀
end

function double_well(x, ħω1::Real, ħω2::Real; param)
    a1 = amu / 2. * (ħω1/ħc/1E10)^2
    a2 = amu / 2. * (ħω2/ħc/1E10)^2
    return - a1*x.*x + a2*x.*x.*x.*x
end

# Set of potentials for fitting
function polyfunc(x, coeffs; E₀::Real, Q₀::Real, poly_order::Int)
    # E0 = param["E0"]
    # Q₀ = param["Q0"]
    # poly_order = Int(param["poly_order"])
    y = 0 .* x # array of size x
    y = y .+ E₀

    for i = 2:poly_order + 1
        y = y .+ coeffs[i].* (x .- Q₀) .^(i-1)
    end

    return y
end

function morse(x, coeffs; E₀, Q₀)
    y = 0 .* x # array of size x
    y = y .+ E₀

    y = y .+ coeffs[1] .* ((1 .- exp.(-coeffs[2].*(x.-Q₀))).^2)

    return y
end

function morse_poly(x, coeffs; E₀, Q₀, poly_order)
    # TODO: CHECK # of coeffs and param["poly_order"], and give a WARNING
    # E₀ = param["E0"]
    # Q₀ = param["Q0"]
    # po = param["poly_order"]

    orders = if (typeof(po) == Int64) Array([po]) else parse.(Int, split(po)) end

    A = abs(coeffs[1])
    a = coeffs[2]

    r₀ = coeffs[3]

    morse = x -> A*(exp(-2*a*(x-Q₀-r₀)) - 2*exp(-a*(x-Q₀-r₀)))

    coeff = zeros(maximum(orders)+1)
    for (index, poly_order) in enumerate(orders)
        coeff[poly_order+1] = coeffs[3+index]
    end

    coeff[end] = abs(coeff[end])

    poly = Poly(coeff)
    tot_derv = polyder(poly, 1) -  A*(-2a)*(exp(2a*r₀) - exp(a*r₀))

    r₁ = Polynomials.roots(tot_derv)
    r₁ = real.(r₁[imag.(r₁) .== 0])
    r₁ = r₁[argmin(abs.(r₁))]

    return morse(x) + poly(x-Q₀-r₁) + E₀ - morse(Q₀) - poly(-r₁)
end

function get_spline(Qs, Es; weight = nothing, smoothness = 0, order = 2)
    if Qs[1] > Qs[end]
        Qs = reverse(Qs)
        Es = reverse(Es)
    end

    if weight == nothing
        spl = Spline1D(Qs, Es, k = order, bc = "extrapolate", s = smoothness)
    else
        spl = Spline1D(Qs, Es, w = weight, k = order, bc = "extrapolate", s = smoothness)
    end
    return spl
end

function get_bspline(Qs, Es)
    # BSplines assume your data is uniformly spaced on the grid
    # Qs, Es have to be eqaully-spaced (Range)
    if Qs[1] > Qs[end]
        Qs = reverse(Qs)
        Es = reverse(Es)
    end
    Qs_range = range(Qs[1], stop = Qs[end], length = length(Qs))

    itp = interpolate(Es, BSpline(Quadratic(Line(OnCell()))))
    etpf = extrapolate(itp, Line())
    setpf = scale(etpf, Qs_range)

    return setpf
end


text:
include("linear_regressor.jl")
include("ridge_regressor.jl")
include("bayesian_regressor.jl")
include("empirical_bayesian_regressor.jl")

ext: jl

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
using AlgebraOfGraphics, Test
using AlgebraOfGraphics: data,
                         mapping,
                         Spec,
                         spec,
                         layers,
                         positional,
                         keyword,
                         dims,
                         run_pipeline

using CategoricalArrays: categorical
using NamedDims
using RDatasets: dataset
using Observables: to_value
using AbstractPlotting: default_palettes

@testset "product" begin
    s = dims() * mapping(1:2, ["a", "b"], color = dims(1))
    ps = pairs(s)
    @test ps[1] == Pair((color = 1,), mapping(1, "a"))
    @test ps[2] == Pair((color = 2,), mapping(2, "b"))
end

@testset "contexts" begin
    mpg = dataset("ggplot2", "mpg")
    d = mapping(:Cyl, :Hwy, color = :Year => categorical)
    s = visual(color = :red, font = 10) + mapping(markersize = :Year)
    res = data(mpg) * d * s
    @test res[1].options == (color = :red, font = 10)

    idx1 = mpg.Year .== 1999
    idx2 = mpg.Year .== 2008

    mappings = pairs.(getproperty.(res, :mapping))
    @test last(mappings[1][1]).value == mapping(mpg[idx1, :Cyl], mpg[idx1, :Hwy]).value
    @test last(mappings[1][2]).value == mapping(mpg[idx2, :Cyl], mpg[idx2, :Hwy]).value
    @test last(mappings[2][1]).value == mapping(mpg[idx1, :Cyl], mpg[idx1, :Hwy], markersize = mpg[idx1, :Year]).value
    @test last(mappings[2][2]).value == mapping(mpg[idx2, :Cyl], mpg[idx2, :Hwy], markersize = mpg[idx2, :Year]).value

    @test first(mappings[1][1]).color isa NamedDimsArray
    @test only(first(mappings[1][1]).color) == 1999
    @test dimnames(first(mappings[1][1]).color) == (:Year,)

    @test first(mappings[1][2]).color isa NamedDimsArray
    @test only(first(mappings[1][2]).color) == 2008
    @test dimnames(first(mappings[1][2]).color) == (:Year,)

    @test first(mappings[2][1]).color isa NamedDimsArray
    @test only(first(mappings[2][1]).color) == 1999
    @test dimnames(first(mappings[2][1]).color) == (:Year,)

    @test first(mappings[2][2]).color isa NamedDimsArray
    @test only(first(mappings[2][2]).color) == 2008
    @test dimnames(first(mappings[2][2]).color) == (:Year,)

    @test_throws ArgumentError data(mpg) * mapping(:Cyll)

    x = rand(5, 3, 2)
    y = rand(5, 3)
    s = dims(1) * mapping(x, y, color = dims(2))

    res = pairs(s)
    for (i, r) in enumerate(res)
        group, st = r
        @test group == (; color = mod1(i, 3))
        xsl = x[:, mod1(i, 3), (i > 3) + 1]
        ysl = y[:, mod1(i, 3)]
        @test st.value == mapping(xsl, ysl).value
    end
end

@testset "compute specs" begin
    wong = default_palettes[:color][]
    t = (x = [1, 2], y = [10, 20], z = [3, 4], c = ["a", "b"])
    d = mapping(:x, :y, color = :c)
    s = visual(:log) * visual(font = 10) + mapping(size = :z)
    ds = data(t) * d
    sl = ds * s
    @test length(sl) == 2

    res = run_pipeline(sl)

    @test res[1].options.color == wong[1]
    @test res[2].options.color == wong[2]
    @test dimnames(res[1].pkeys.color) == (:c,)
    @test dimnames(res[2].pkeys.color) == (:c,)
    @test res[1].mapping.value == mapping([1], [10]).value
    @test res[2].mapping.value == mapping([2], [20]).value

    @test res[3].options.color == wong[1]
    @test res[4].options.color == wong[2]
    @test dimnames(res[3].pkeys.color) == (:c,)
    @test dimnames(res[4].pkeys.color) == (:c,)
    @test res[3].mapping.value == mapping([1], [10], size = [3]).value
    @test res[4].mapping.value == mapping([2], [20], size = [4]).value
end


text:
function LutherIntegrator(arg0::jdouble)
    return LutherIntegrator((jdouble,), arg0)
end

function get_a(obj::LutherIntegrator)
    return jcall(obj, "getA", Vector{Vector{jdouble}}, ())
end

function get_b(obj::LutherIntegrator)
    return jcall(obj, "getB", Vector{jdouble}, ())
end

function get_c(obj::LutherIntegrator)
    return jcall(obj, "getC", Vector{jdouble}, ())
end


text:
using Documenter, BufferedFiles

makedocs(sitename = "BufferedFiles Documentation")


text:
function HACCache(X::AbstractMatrix{T}; prewhiten::Bool = false) where {T<:Int}
    HACCache(convert(Matrix{WFLOAT}, X), prewhiten = prewhiten)
end

function HACCache(X::AbstractMatrix{T}; prewhiten::Bool = false) where {T<:Real}
    nr, p = size(X)
    TYPE = prewhiten ? Prewhitened() : Unwhitened()
    n = prewhiten ? nr-1 : nr
    if prewhiten
    return HACCache(TYPE,
                     X,
                     Array{T}(undef, n, p),
                     Array{T}(undef, n, p),
                     Array{T}(undef, n-1, p),
                     Array{T}(undef, n-1, p),
                     Array{T}(undef, 1, p),
                     Array{T}(undef, p, p),
                     Array{T}(undef, p, p),
                     Array{T}(undef, p, p),   ## This sometime host the ldiv! which may have larger type
                     Array{T}(undef, n-1),
                     Array{T}(undef, p),
                     Array{T}(undef, p),
                     Array{T}(undef, n, p))
    else
        return HACCache(TYPE,
                        X,
                        Array{T}(undef, 0, 0),
                        Array{T}(undef, n, p),
                        Array{T}(undef, n-1, p),
                        Array{T}(undef, n-1, p),
                        Array{T}(undef, 1, p),
                        Array{T}(undef, p, p),
                        Array{T}(undef, p, p),
                        Matrix(one(T)I, p, p),
                        Array{T}(undef, n-1),
                        Array{T}(undef, p),
                        Array{T}(undef, p),
                        Array{T}(undef, 0, 0))
    end
end

function HACCache(X::AbstractMatrix, k::HAC; kwargs...)
    ip = isprewhiten(k)
    HACCache(X, prewhiten = ip; kwargs...)
end

check_cache_consistenty(k::HAC, cache::HACCache{T}) where T<:Prewhitened = isprewhiten(k) ? nothing : error("Inconstent cache type")
check_cache_consistenty(k::HAC, cache::HACCache{T}) where T<:Unwhitened = !isprewhiten(k) ? nothing : error("Inconstent cache type")

Optimal() = Optimal{Andrews}()

TruncatedKernel(;prewhiten=false) = TRK(Optimal(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
BartlettKernel(;prewhiten=false) = BTK(Optimal(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
ParzenKernel(;prewhiten=false) = PRK(Optimal(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
TukeyHanningKernel(;prewhiten=false) = THK(Optimal(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
QuadraticSpectralKernel(;prewhiten=false) = QSK(Optimal(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)

BartlettKernel(x::Type{NeweyWest};prewhiten=false) = BTK(Optimal{NeweyWest}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
ParzenKernel(x::Type{NeweyWest};prewhiten=false) = PRK(Optimal{NeweyWest}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
QuadraticSpectralKernel(x::Type{NeweyWest};prewhiten=false) = QSK(Optimal{NeweyWest}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
TukeyHanningKernel(x::Type{NeweyWest};prewhiten=false) = error("Newey-West optimal bandwidth does not support TukeyHanningKernel")
TruncatedKernel(x::Type{NeweyWest};prewhiten=false) = error("Newey-West optimal bandwidth does not support TuncatedKernel")

TruncatedKernel(x::Type{Andrews};prewhiten=false) = TRK(Optimal{Andrews}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
BartlettKernel(x::Type{Andrews};prewhiten=false) = BTK(Optimal{Andrews}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
ParzenKernel(x::Type{Andrews};prewhiten=false) = PRK(Optimal{Andrews}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
TukeyHanningKernel(x::Type{Andrews};prewhiten=false) = THK(Optimal{Andrews}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)
QuadraticSpectralKernel(x::Type{Andrews};prewhiten=false) = QSK(Optimal{Andrews}(), Array{WFLOAT}(undef,1), Array{WFLOAT}(undef,0), prewhiten)

TruncatedKernel(bw::Number;prewhiten=false) = TRK(Fixed(), [float(bw)], Array{WFLOAT}(undef,0), prewhiten)
BartlettKernel(bw::Number;prewhiten=false) = BTK(Fixed(), [float(bw)], Array{WFLOAT}(undef,0), prewhiten)
ParzenKernel(bw::Number;prewhiten=false) = PRK(Fixed(), [float(bw)], Array{WFLOAT}(undef,0), prewhiten)
TukeyHanningKernel(bw::Number;prewhiten=false) = THK(Fixed(), [float(bw)], Array{WFLOAT}(undef,0), prewhiten)
QuadraticSpectralKernel(bw::Number;prewhiten=false) = QSK(Fixed(), [float(bw)], Array{WFLOAT}(undef,0), prewhiten)

bandwidth(k::HAC{G}, X::AbstractMatrix) where {G<:Fixed} = k.bw
bandwidth(k::HAC{Optimal{G}}, X::AbstractMatrix) where {G<:Andrews} = bwAndrews(k, )

function bandwidth(k::QuadraticSpectralKernel, X::AbstractMatrix)
    return k.bw(X, k)
end

isprewhiten(k::HAC) = k.prewhiten

getcovindeces(k::T, n) where T<:QuadraticSpectralKernel = Iterators.filter(x -> x!=0, -n:n)
getcovindeces(k::HAC, n) = Iterators.filter(x -> x!=0, -floor(Int, k.bw[1]):floor(Int, k.bw[1]))

function Γ!(cache, j)
    X = cache.XX
    T, p = size(X)
    Q = fill!(cache.Q, zero(eltype(X)))
    if j >= 0
        for h=1:p, s = 1:h
            for t = j+1:T
                @inbounds Q[s, h] = Q[s, h] + X[t, s]*X[t-j, h]
            end
        end
    elseif j<0
        for h=1:p, s = 1:h
            for t = -j+1:T
                @inbounds Q[s,h] = Q[s ,h] + X[t+j, s]*X[t,h]
            end
        end
    end
    return cache.Q
end

function demean!(cache::HACCache, X, ::Type{Val{true}})
    sum!(cache.μ, X)
    rmul!(cache.μ, 1/size(X,1))
    cache.X_demean .= X .- cache.μ
end

function demean!(cache::HACCache, X, ::Type{Val{false}})
    copyto!(cache.X_demean, X)
end

prewhiten!(cache::HACCache{T}) where T<:Unwhitened = copyto!(cache.XX, cache.X_demean)
prewhiten!(cache::HACCache{T}) where T<:Prewhitened = fit_var!(cache)
swhiten!(cache::HACCache{T}) where T<:Unwhitened = nothing

function swhiten!(cache::HACCache{T}) where T<:Prewhitened
    fill!(cache.Q, zero(eltype(cache.Q)))
    for i = 1:size(cache.Q, 2)
        cache.Q[i,i] = one(eltype(cache.Q))
    end
    v = ldiv!(qr(I-cache.D'), cache.Q)
    cache.V .= v*cache.V*v'
end

# makecholesky!(cache, ::Type{Nothing}) = nothing

# function makecholesky!(cache, ::Type{Cholesky})
#     chol = LinearAlgebra.cholesky(Symmetric(cache.V), check = false)
#     copyto!(cache.chol.UL.data, chol.UL.data)
#     copyto!(cache.chol.U.data, chol.U.data)
#     copyto!(cache.chol.L.data, chol.L.data)
# end

# function makecholesky!(cache, ::Type{PositiveFactorizations.Positive})
#     chol = LinearAlgebra.cholesky(Positive, Symmetric(cache.V))
#     copyto!(cache.chol.UL.data, chol.UL.data)
#     copyto!(cache.chol.U.data, chol.U.data)
#     copyto!(cache.chol.L.data, chol.L.data)
# end

##############################################################################
##
## Kernel methods
##
##############################################################################

kernel(k::HAC, x::Real) = isnan(x) ? (return 1.0) : kernel(k, float(x))
kernel(k::TruncatedKernel, x::Real)    = (abs(x) <= 1.0) ? 1.0 : 0.0
kernel(k::BartlettKernel, x::Real)     = (abs(x) <= 1.0) ? (1.0 - abs(x)) : 0.0
kernel(k::TukeyHanningKernel, x::Real) = (abs(x) <= 1.0) ? 0.5 * (1.0 + cospi(x)) : 0.0

function kernel(k::ParzenKernel, x::Real)
    ax = abs(x)
    if ax > 1.0
        0.0
    elseif ax <= 0.5
        1.0 - 6.0 * ax^2 + 6.0 * ax^3
    else
        2.0 * (1.0 - ax)^3
    end
end

function kernel(k::QuadraticSpectralKernel, x::Real)
    iszero(x) ? 1.0 : (z = 1.2*π*x; 3*(sin(z)/z-cos(z))*(1/z)^2)
end

function setupkernelweights!(k, p, xtype)
    if isempty(k.weights)
        for j in 1:p
            push!(k.weights, one(xtype))
        end
    elseif all(iszero.(k.weights))
        fill!(k.weights, one(xtype))
    end
end

##############################################################################
##
## Fit functions
##
##############################################################################

 function fit_var!(cache::HACCache)
     X, Y, Z, u, D = cache.XX, cache.YY, cache.X_demean, cache.u, cache.D
     n, p = size(Z)
     @inbounds for j in 1:p, i = 1:n-1
         X[i,j] = Z[i,  j]
         Y[i,j] = Z[i+1,j]
     end
     QX = qr(X)
    ldiv!(D, QX, convert(Matrix{eltype(QX)}, Y))
     @inbounds for j in 1:p, i = 1:n-1
         Y[i,j] = Z[i+1,j]
     end
     mul!(u, X, D)
     broadcast!(-, X, Y, u)
 end

 function fit_ar!(cache)
     ## Estimate
     ##
     ## y_{t,j} = ρ y_{t-1,j} + ϵ
     σ⁴ = cache.σ⁴
     ρ = cache.ρ
     U = cache.U
     n, p = size(cache.XX)
     lag!(cache)
     Y = cache.Y_lagged
     X = cache.X_lagged
     for j in 1:p
         y = view(Y, :, j)
         x = view(X, :, j)
         x .= x .- mean(x)
         y .= y .- mean(y)
         ρ[j] = sum(broadcast!(*, cache.U, x, y))/sum(abs2, x)
         copyto!(U, y)
         x .= x.*ρ[j]
         broadcast!(-, U, U, x)
         σ⁴[j]  = (dot(U, U)/(n-1))^2
     end
 end

 function lag!(cache)
     ## This construct two matrices
     ## Z_lagged we store X_demean[1:n-1, :]
     nl, pl = size(cache.Y_lagged)
     n, p  = size(cache.XX)
     for ic in 1:p
         for i = 2:n
             @inbounds cache.Y_lagged[i-1, ic] = cache.XX[i, ic]
             @inbounds cache.X_lagged[i-1, ic] = cache.XX[i-1, ic]
         end
     end
  end

##############################################################################
##
## Optimal bandwidth
##
##############################################################################

optimal_bw!(cache, k::HAC, optype::T) where T<:NeweyWest = bwNeweyWest(cache, k)
optimal_bw!(cache, k::HAC, opttype::T) where T<:Andrews = bwAndrews(cache, k)

function bwAndrews(cache, k::HAC)
    isempty(k.weights) && (fill!(k.weights, 1.0))
    n, p  = size(cache.XX)
    a1, a2 = getalpha!(cache, k.weights)
    k.bw[1] = bw_andrews(k, a1, a2, n)
end

## ---> Andrews Optimal bandwidth <---
d_bw_andrews = Dict(:TruncatedKernel         => :(0.6611*(a2*n)^(0.2)),
                    :BartlettKernel          => :(1.1447*(a1*n)^(1/3)),
                    :ParzenKernel            => :(2.6614*(a2*n)^(0.2)),
                    :TukeyHanningKernel      => :(1.7462*(a2*n)^(0.2)),
                    :QuadraticSpectralKernel => :(1.3221*(a2*n)^(0.2)))

for kerneltype in [:TruncatedKernel, :BartlettKernel, :ParzenKernel, :TukeyHanningKernel, :QuadraticSpectralKernel]
    @eval $:(bw_andrews)(k::($kerneltype), a1, a2, n) = $(d_bw_andrews[kerneltype])
end

function getalpha!(cache, w)
    fit_ar!(cache)
    σ⁴, ρ = cache.σ⁴, cache.ρ
    nm = 4.0.*(ρ.^2).*σ⁴./(((1.0.-ρ).^6).*((1.0.+ρ).^2))
    dn = σ⁴./(1.0.-ρ).^4
    α₁ = sum(w.*nm)/sum(w.*dn)
    nm = 4.0.*(ρ.^2).*σ⁴./((1.0.-ρ).^8)
    α₂ = sum(w.*nm)/sum(w.*dn)
    return α₁, α₂
end

function bwNeweyWest(cache, k::HAC)
    n, p = size(cache.XX)
    l = getrates(cache, k)
    w = k.weights
    xm = cache.XX*w
    a = map(j -> dot(xm[1:n-j], xm[j+1:n])/n, 0:l)::Array{Float64, 1}
    aa = view(a, 2:l+1)
    a0 = a[1] + 2*sum(aa)
    a1 = 2*sum((1:l) .* aa)
    a2 = 2*sum((1:l).^2 .* aa)
    k.bw[1] = bwnw(k, a0, a1, a2)*(n+isprewhiten(k))^growthrate(k)
end

function getrates(cache, k)
    n, p = size(cache.X_demean)
    lrate = lagtruncation(k)
    adj = isprewhiten(k) ? 3 : 4
    floor(Int, adj*(n/100)^lrate)
end

@inline bwnw(k::BartlettKernel, s0, s1, s2) = 1.1447*((s1/s0)^2)^growthrate(k)
@inline bwnw(k::ParzenKernel, s0, s1, s2) = 2.6614*((s2/s0)^2)^growthrate(k)
@inline bwnw(k::QuadraticSpectralKernel, s0, s1, s2) = 1.3221*((s2/s0)^2)^growthrate(k)

## --> Newey-West Optimal bandwidth <---
@inline growthrate(k::HAC) = 1/5
@inline growthrate(k::BartlettKernel) = 1/3

@inline lagtruncation(k::BartlettKernel) = 2/9
@inline lagtruncation(k::ParzenKernel) = 4/25
@inline lagtruncation(k::QuadraticSpectralKernel) = 2/25


text:
################################################################################
# CutPruners
# A package to manage polyhedral convex functions
################################################################################

VERSION < v"0.7.0-beta2.199" && __precompile__()

module CutPruners

using Compat, Compat.LinearAlgebra, Compat.SparseArrays
using MathProgBase

# Redudancy checking
include("redund.jl")

include("abstract.jl")
include("avg.jl")
include("decay.jl")
include("dematos.jl")
include("exact.jl")

end # module


text:
"""
    FEMMDeforNonlinearModule

Module for nonlinear deformation models. Standard isoparametric finite elements.
"""
module FEMMDeforNonlinearModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, gradN!, nodesperelem, manifdim
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
import ..FEMMDeforNonlinearBaseModule: AbstractFEMMDeforNonlinear
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed2DAxisymm, nstressstrain, nthermstrain, Blmat!
import FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, locjac!, add_nnt_ut_only!
import FinEtoolsDeforLinear.MatDeforModule: rotstressvec!
import FinEtools.MatModule: massdensity
import ..MatDeforNonlinearModule: AbstractMatDeforNonlinear, tangentmoduli!, update!, newstate
import FinEtools.SurfaceNormalModule: SurfaceNormal, updatenormal!
import LinearAlgebra: Transpose, mul!
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: norm, dot

"""
    FEMMDeforNonlinear{MR<:AbstractDeforModelRed,  S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear} <: AbstractFEMMDeforNonlinear

Class for nonlinear deformation finite element modeling machine.
"""
mutable struct FEMMDeforNonlinear{MR<:AbstractDeforModelRed,  S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear} <: AbstractFEMMDeforNonlinear
    mr::Type{MR} # model reduction type
    integdomain::IntegDomain{S, F} # integration domain data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    statev::Vector{Vector{FFltVec}} # vector of state vectors, one for each element, and then for each integration point
end

"""
    FEMMDeforNonlinear(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}

Constructor of nonlinear deformation finite element modeling machine.

The material coordinate system is the default.
"""
function FEMMDeforNonlinear(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (integdomain.axisymmetric) || (mr != DeforModelRed2DAxisymm) "Axially symmetric requires axisymmetric to be true"
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(integdomain);
    statev = Vector{FFltVec}[]
    for i in 1:count(integdomain.fes)
        push!(statev, [newstate(material) for j in 1:npts])
    end
    return FEMMDeforNonlinear(mr, integdomain, CSys(manifdim(integdomain.fes)), material, statev)
end

"""
    FEMMDeforNonlinear(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}

Constructor of nonlinear deformation finite element modeling machine.
"""
function FEMMDeforNonlinear(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (integdomain.axisymmetric) || (mr != DeforModelRed2DAxisymm) "Axially symmetric requires axisymmetric to be true"
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(integdomain);
    statev = Vector{FFltVec}[]
    for i in 1:count(integdomain.fes)
        push!(statev, [newstate(material) for j in 1:npts])
    end
    return FEMMDeforNonlinear(mr, integdomain, mcsys, material, statev)
end

end


text:
mutable struct State
  field::Array{Complex{Float64}, 2}
  fieldIterations::Array{Int64, 2}

  iterations::Int64

  res::Tuple{Int64, Int64}

  bounds::Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}
  lengths::Tuple{Float64, Float64}

  delta::Tuple{Float64, Float64}
  deltas::Array{Complex{Float64}, 2}

  op::Function
end

""" Propagate changes in the core state (res, bounds) to dependents"""
function updateState!(state::State)
  state.field .= 0
  state.fieldIterations .= 0

  state.lengths = map(b -> b[2] - b[1], state.bounds)

  state.delta = state.lengths ./ state.res
  for i=1:state.res[1]
    for j=1:state.res[2]
      state.deltas[i, j] = complex(state.bounds[1][1] + i * state.delta[1], state.bounds[2][1] + j * state.delta[2])
    end
  end
end

function State(; iterations=20.0, res=(32, 32), bounds=((-2, 2), (-2, 2)), op=(z,c)->z^2+c)
  field = Array{Complex{Float64}}(undef, res)
  fieldIterations = Array{Int64}(undef, res)

  deltas = Array{Complex{Float64}}(undef, res)

  state = State(field, fieldIterations, iterations, res, bounds, (-1.0, -1.0), (-1.0, -1.0), deltas, op)

  updateState!(state)

  return state
end

max_stars_repo_name: benchislett/Skunkworks

max_issues_repo_name: benchislett/Skunkworks

max_forks_repo_name: benchislett/Skunkworks


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Graphviz"
version = v"2.42.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www2.graphviz.org/Packages/stable/portable_source/graphviz-$(version).tar.gz",
                  "8faf3fc25317b1d15166205bf64c1b4aed55a8a6959dcabaa64dbad197e47add"),

    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/graphviz-*/

if [[ "${target}" == *-mingw* ]]; then
    # We need a `regex.h` header.  MinGW doesn't have one,
    # let's use `pcreposix.h` instead.
    cp ${prefix}/include/pcreposix.h ${prefix}/include/regex.h

    # Apply some fun patches
    atomic_patch -p1 ../patches/0003-sfsetbuf_c_Stat_t_no_st_blksize.patch
    atomic_patch -p1 ../patches/0004-win32_dllexport_dllimport.patch
    atomic_patch -p1 ../patches/0005-missing_libs.patch
    atomic_patch -p1 ../patches/0006-export_neatogen.patch
    atomic_patch -p1 ../patches/0007-remove_missing_def.patch
    atomic_patch -p1 ../patches/0008-export_gvc.patch

    export LDFLAGS="-lpcreposix -lexpat"
    export EXTRA_LDFLAGS="-no-undefined"

    # Remove wrong libtool archives
    rm ${prefix}/lib/libharfbuzz*.la
fi

# Do not build with -ffast-math
atomic_patch -p1 ../patches/1001-no-ffast-math.patch

# Rebuild the configure script
autoreconf -fiv

# Apply patch to build a native `mkdefs` utility that can be run within the
# build environment.
atomic_patch -p1 ../patches/0001-gvpr-build-native-mkdefs.patch

# This patch disable generation of dot's configuration
atomic_patch -p1 ../patches/0002-do-not-build-dot-config.patch

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc} LDFLAGS="${LDFLAGS} ${EXTRA_LDFLAGS}"
make install

if [[ "${target}" == *-mingw* ]]; then
    # Cover up the traces of the hack
    rm ${prefix}/include/regex.h
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(filter!(p -> !isa(p, FreeBSD) & !isa(p, Windows), supported_platforms()))

# The products that we will ensure are always built
products = [
    LibraryProduct("libcdt", :libcdt),
    LibraryProduct("libcgraph", :libcgraph),
    LibraryProduct("libgvc", :libgvc),
    LibraryProduct("libgvpr", :libgvpr),
    LibraryProduct("liblab_gamut", :liblab_gamut),
    LibraryProduct("libpathplan", :libpathplan),
    LibraryProduct("libxdot", :libxdot),
    ExecutableProduct("acyclic", :acyclic),
    ExecutableProduct("bcomps", :bcomps),
    ExecutableProduct("ccomps", :ccomps),
    ExecutableProduct("circo", :circo),
    ExecutableProduct("cluster", :cluster),
    ExecutableProduct("dijkstra", :dijkstra),
    ExecutableProduct("dot", :dot),
    ExecutableProduct("dot2gxl", :dot2gxl),
    ExecutableProduct("dot_builtins", :dot_builtins),
    ExecutableProduct("edgepaint", :edgepaint),
    ExecutableProduct("fdp", :fdp),
    ExecutableProduct("gc", :gc),
    ExecutableProduct("gml2gv", :gml2gv),
    ExecutableProduct("graphml2gv", :graphml2gv),
    ExecutableProduct("gv2gml", :gv2gml),
    ExecutableProduct("gv2gxl", :gv2gxl),
    ExecutableProduct("gvcolor", :gvcolor),
    ExecutableProduct("gvgen", :gvgen),
    ExecutableProduct("gvmap", :gvmap),
    ExecutableProduct("gvmap.sh", :gvmap_sh),
    ExecutableProduct("gvpack", :gvpack),
    ExecutableProduct("gvpr", :gvpr),
    ExecutableProduct("gxl2dot", :gxl2dot),
    ExecutableProduct("gxl2gv", :gxl2gv),
    ExecutableProduct("mm2gv", :mm2gv),
    ExecutableProduct("neato", :neato),
    ExecutableProduct("nop", :nop),
    ExecutableProduct("osage", :osage),
    ExecutableProduct("patchwork", :patchwork),
    ExecutableProduct("prune", :prune),
    ExecutableProduct("sccmap", :sccmap),
    ExecutableProduct("sfdp", :sfdp),
    ExecutableProduct("tred", :tred),
    ExecutableProduct("twopi", :twopi),
    ExecutableProduct("unflatten", :unflatten),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Cairo_jll", uuid="83423d85-b0ee-5818-9007-b63ccbeb887a")),
    Dependency(PackageSpec(name="Expat_jll", uuid="2e619515-83b5-522b-bb60-26c02a35a201")),
    Dependency(PackageSpec(name="Pango_jll", uuid="36c8627f-9965-5494-a995-c6b170f724f3")),
    # PCRE is needed only for Windows.  Maybe it's only a build dependency?
    # Dependency(PackageSpec(name="PCRE_jll",  uuid="2f80f16e-611a-54ab-bc61-aa92de5b98fc")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jpsamaroo/Yggdrasil

max_issues_repo_name: jpsamaroo/Yggdrasil

max_forks_repo_name: jpsamaroo/Yggdrasil


text:
module EchogramImages

using Images
using EchogramColorSchemes
using ColorSchemes

export imagesc, mat2gray, equalizedbins, quantize

"""
    imagesc(A;
             vmin = nothing,
             vmax = nothing,
             cmap= nothing,
             size=(480,640))

A is any numeric array

cmap can be a list of Colors or come from EchogramColorSchemes.jl,
ColorSchemes.jl or PerceptualColourMaps.jl. The default `nothing`
currently means use EK500 colours.

Set size = nothing for full resolution.

vmin and vmax are minimum and maximum values.

"""
function imagesc(A::AbstractArray;
                  vmin = nothing,
                  vmax = nothing,
                  cmap= nothing,
                  size=(480,640))

    if cmap == nothing
        cmap = addwhite(EK500)
    end

    if isa(cmap, ColorScheme)
        cmap = cmap.colors
    end

    if size != nothing
        A = imresize(A, size)
    end

    if vmin == nothing
        vmin = minimum(A[.!isnan.(A)])
    end

    if vmax == nothing
        vmax = maximum(A[.!isnan.(A)])
    end

    vmin = float(vmin)
    vmax = float(vmax)

    n = length(cmap)

    g = x -> (clamp(isnan(x) ? vmin : x, vmin, vmax) - vmin) / (vmax - vmin)
    #f = s->clamp(round(Int, (n-1)*g(s))+1, 1, n)
    f = s->cmap[trunc(Int, (n-1)*g(s)) + 1]

    f.(A)

end

"""
    mat2gray(A; amin=nothing, amax=nothing)

Convert matrix `A` to grayscale with pixel values in the range 0.0
(black) to 1.0 (white).

NB Intentionally, NaNs get replaced by 0.0 to facilitate the use of
the wider Julia Images ecosystem.

"""
function mat2gray(A; amin=nothing, amax=nothing)

    B = A[.!isnan.(A)]

    if amin == nothing
        amin = minimum(B)
    end
    if amax == nothing
        amax = maximum(B)
    end
    f = scaleminmax(amin, amax)
    C = f.(A)
    C[isnan.(C)] .= 0.0
    return C
end

"""
    equalizedbins(A; nbins = 256, amin=nothing, amax=nothing)

Using histogram equalisation, return the bin edges (the cut off
points) that would allow `A` to be quantized into `nbins`.

"""
function equalizedbins(A; nbins = 256, amin=nothing, amax=nothing)
    a = sort(vec(A))
    if amin != nothing
        filter!(x-> x>amin , a)
        # Reserve a bin for everything below minA
        nbins = nbins -1
    end
    if amax != nothing
        filter!(x-> x<amax , a)
        # Reserve a bin for everything above maxA
        nbins = nbins -1
    end

    # Equalise bin spacing
    d = length(a) / nbins
    r = [a[ceil(Int,i)] for i in d:d:d*(nbins-1)]

    if amin != nothing
        # if specified, minA is the first bin dimension
        r = vcat([amin],r)
    end
    if amax != nothing
        # if specified, maxA is the last bin dimension
        r = vcat(r, [amax])
    end
    return r

end

"""

    quantize(A, edges)

Quantise the 2D array`A` into `length(edges)` levels where `edges` is
the list of bin cut-off points.

Quantization, in mathematics and digital signal processing, is the
process of mapping input values from a large set (often a continuous
set) to output values in a (countable) smaller set, often with a
finite number of elements.

"""
function quantize(A, edges)
    m,n = size(A)
    if length(edges) < 256
        B= Matrix{UInt8}(undef,m,n)
    else
        B= Matrix(undef,m,n)
    end
    B .= 0
    for i in 1:length(edges)
        B[A .> edges[i]] .= i
    end
    return B
end

end # module


text:
using DRIPs
using Test
using LinearAlgebra

@testset "MW2009.jl" begin
	## Parameters
	ρ  = 0.95;
	σq = 0.01;
	σz = 11.8*σq;
	κ  = 3;
	ξ  = 1;
	α  = 1 - 0.15;

	## Primitives of drip
	L  = 21; # length of trunction
	A  = [zeros(1,L);[Matrix(I,L-1,L-1);zeros(1,L-1)]']; # MW truncate the state space with linear irfs of length 20
	Qq = zeros(L,1); Qq[1]=σq;
	Qz = zeros(L,1); Qz[1]=σz;
	H  = zeros(L,1); H[1:21] = Array(1:-1/20:0);

	## Function for solving the fixed point of the aggregate shock given ω
	# Note that for a certain range of ω the fixed point is not unique:
	# there is a fixed point with zero capacity and one with positive capacity.
	# In this range, we would like to find the fixed point in which the firms are
	# acqurining a positive amount of information.

	function agg_drip(ω,A,Qq,        #primitives of drip except for H because H is endogenous
	                  α,             #strategic complementarity
	                  H;             #state space rep. of q
	                  β     = 1,     #optional: discount factor, MW's parameterization implies β = 1
	                  H0    = H,     #optional: initial guess for HΔ (H is the true solution when α=0)
	                  maxit = 10000, #optional: max number of iterations for GE code
	                  tol   = 1e-4,  #optional: tolerance for iterations
	                  w     = 1)     #optional: update weight for RI
	    # set primitives
	    errmin= 1;
	    err   = 1;
	    iter  = 0;
	    L     = length(H);
	    while (err > tol) & (iter < maxit)
	            if iter == 0
	                global agg  = Drip(ω,β,A,Qq,H0;w = w);
	            else
	                global agg  = Drip(ω,β,A,Qq,H0;Ω0 = agg.ss.Ω , Σ0 = agg.ss.Σ_1,w = w);
	            end

	            XFUN(jj) = ((I-agg.ss.K*agg.ss.Y')*agg.A)^jj * (agg.ss.K*agg.ss.Y') * (agg.A')^jj
	            X = DRIPs.infinitesum(XFUN; maxit=200, start = 0);  #E[x⃗]=X×x⃗

	            XpFUN(jj) = α^jj * X^(jj)
	            Xp = DRIPs.infinitesum(XpFUN; maxit=200, start = 0);

	            H1 = (1-α)*Xp'*H;

	            err= 0.5*norm(H1-H0,2)/norm(H0)+0.5*err;
	            if DRIPs.capacity(agg) < 1e-2 # perturb the initial guess if solution is the zero capacity one
	                H0 = H0+rand(L).*(H-H0);
	            else # store the solution if it has positive capacity
	                H0 = H1;
	                if err < errmin
	                    global aggmin = agg;
	                    errmin = err;
	                end
	            end
	            iter += 1;
	    end
	    return(aggmin,min(err,errmin))
	end;
	println("Solving Mackowiack and Wiederholt (2009) ...")
	## Solve for κ = 3
	agg, err = agg_drip(2.4σq^2,A,Qq,α,H; H0 = rand(L), maxit = 500, w = 0.95)
	@test agg.ss.err < 1e-4
	agg, err = agg_drip(2.82σq^2,A,Qq,α,H; H0 = rand(L), maxit = 500, w = 0.95)
	@test agg.ss.err < 1e-4
	agg, err = agg_drip(0.01*σq^2,A,Qq,α,H; H0 = rand(L), maxit = 500, w = 0.95)
	@test agg.ss.err < 1e-4
	agg, err = agg_drip(10*σq^2,A,Qq,α,H; H0 = rand(L), maxit = 500, w = 0.95)
	@test agg.ss.err < 1e-4
	idi  = Drip(2.4σq^2,1,A,Qz,H,w = 0.9);
	@test idi.ss.err < 1e-4
end


text:
using Test
using st16SseqJuliaTools
using BioSequences
@testset "Utils" begin
    s1 = dna"CGCAGTCCATGCGATCATACGCAAGAAGGAAGGTCCCCATACACCGACGCACCAGTTTACACGCCGTATGCATAAACGAGCTGCACAAACGAGAGTGCTTGAACTGGACCTCTAGTTCCTCTACAAAGAACAGGTTGACCTGTCGCGAAGTTGCCTTGCCTAGATGCAATGTCGGACGTATTACTTTTGCCTCAACGGCTCCTGCTTTCGCTGAAACCCAAGACAGGCAACAGTAACCGCCTTTTGAAGGCGAGTCCTTCGTCTGTGACTAACTGTGCCAAATCGTCTTCCAAACTCCTA"
    s2 = dna"ACAAAGAACAGGTTGACCTGTCGCGAAGTTGCCTTGCCTAGATGCAATGTCGGACGTATTACTTTTGCCTCAACGGCTCCTGCTTTCGCTGAAACCCAAGACAGGCAACAGTAACCGCCTTTAAAAAAAAAAAAAAAAAAAa"
    @test (0.984, 0.873) == check_pair(s1,s2)
end


text:
@testset "unsqueeze" begin
    x = randn(2, 3, 2)
    @test @inferred(unsqueeze(x; dims=1)) == reshape(x, 1, 2, 3, 2)
    @test @inferred(unsqueeze(x; dims=2)) == reshape(x, 2, 1, 3, 2)
    @test @inferred(unsqueeze(x; dims=3)) == reshape(x, 2, 3, 1, 2)
    @test @inferred(unsqueeze(x; dims=4)) == reshape(x, 2, 3, 2, 1)

    @test unsqueeze(dims=2)(x) == unsqueeze(x, dims=2)
end

@testset "stack and unstack" begin
    x = randn(3,3)
    stacked = stack([x, x], dims=2)
    @test size(stacked) == (3,2,3)

    stacked_array=[ 8 9 3 5; 9 6 6 9; 9 1 7 2; 7 4 10 6 ]
    unstacked_array=[[8, 9, 9, 7], [9, 6, 1, 4], [3, 6, 7, 10], [5, 9, 2, 6]]
    @test unstack(stacked_array, dims=2) == unstacked_array
    @test stack(unstacked_array, dims=2) == stacked_array
    @test stack(unstack(stacked_array, dims=1), dims=1) == stacked_array
end

@testset "batch and unbatch" begin
    stacked_array=[ 8 9 3 5
                    9 6 6 9
                    9 1 7 2
                    7 4 10 6 ]
    unstacked_array=[[8, 9, 9, 7], [9, 6, 1, 4], [3, 6, 7, 10], [5, 9, 2, 6]]
    @test unbatch(stacked_array) == unstacked_array
    @test batch(unstacked_array) == stacked_array

    # no-op for vector of non-arrays
    @test batch([1,2,3]) == [1,2,3]
    @test unbatch([1,2,3]) == [1,2,3]

    # generic iterable
    @test batch(ones(2) for i=1:3) == ones(2, 3)
    @test unbatch(ones(2, 3)) == [ones(2) for i=1:3]

    @testset "tuple" begin
        @test batch([(1,2), (3,4)]) == ([1,3], [2,4])
        @test batch([([1,2], [3,4]), ([5,6],[7,8])]) == ([1 5
                                                          2 6],
                                                         [3 7
                                                          4 8])
    end

    @testset "named tuple" begin
        @test batch([(a=1,b=2), (a=3,b=4)]) == (a=[1,3], b=[2,4])
        @test batch([(a=1,b=2), (b=4,a=3)]) == (a=[1,3], b=[2,4])
        nt = [(a=[1,2], b=[3,4]), (a=[5,6],b=[7,8])]
        @test batch(nt) == (a = [1 5
                                2 6],
                            b = [3 7
                                4 8])
    end

    @testset "dict" begin
        @test batch([Dict(:a=>1,:b=>2), Dict(:a=>3,:b=>4)]) == Dict(:a=>[1,3], :b=>[2,4])
        @test batch([Dict(:a=>1,:b=>2), Dict(:b=>4,:a=>3)]) == Dict(:a=>[1,3], :b=>[2,4])
        d = [Dict(:a=>[1,2], :b=>[3,4]), Dict(:a=>[5,6],:b=>[7,8])]
        @test batch(d) == Dict(:a => [1 5
                                      2 6],
                                :b => [3 7
                                      4 8])
    end
end

@testset "flatten" begin
    x = randn(Float32, 10, 10, 3, 2)
    @test size(flatten(x)) == (300, 2)
end

@testset "normalise" begin
    x = randn(Float32, 3, 2, 2)
    @test normalise(x) == normalise(x; dims=3)
end

@testset "chunk" begin
    cs = chunk(collect(1:10), 3)
    @test length(cs) == 3
    @test cs[1] == [1, 2, 3, 4]
    @test cs[2] == [5, 6, 7, 8]
    @test cs[3] == [9, 10]
end

@testset "group_counts" begin
    d = group_counts(['a','b','b'])
    @test d == Dict('a' => 1, 'b' => 2)
end

@testset "rpad" begin
    @test rpad([1, 2], 4, 0) == [1, 2, 0, 0]
    @test rpad([1, 2, 3], 2, 0) == [1,2,3]
end

@testset "batchseq" begin
    bs = batchseq([[1, 2, 3], [4, 5]], 0)
    @test bs[1] == [1, 4]
    @test bs[2] == [2, 5]
    @test bs[3] == [3, 0]

    bs = batchseq([[1, 2, 3], [4, 5]], -1)
    @test bs[1] == [1, 4]
    @test bs[2] == [2, 5]
    @test bs[3] == [3, -1]
end


text:
module Reqs

using ..Types

# representing lines of REQUIRE files

abstract Line
immutable Comment <: Line
    content::String
end
immutable Requirement <: Line
    content::String
    package::String
    versions::VersionSet
    system::Vector{String}

    function Requirement(content::String)
        fields = split(replace(content, r"#.*$", ""))
        system = String[]
        while !isempty(fields) && fields[1][1] == '@'
            push!(system,shift!(fields)[2:end])
        end
        isempty(fields) && error("invalid requires entry: $content")
        package = shift!(fields)
        all(field->ismatch(Base.VERSION_REGEX, field), fields) ||
            error("invalid requires entry for $package: $content")
        versions = [ convert(VersionNumber, field) for field in fields ]
        issorted(versions) || error("invalid requires entry for $package: $content")
        new(content, package, VersionSet(versions), system)
    end
    function Requirement(package::String, versions::VersionSet, system::Vector{String}=String[])
        content = ""
        for os in system
            content *= "@$os "
        end
        content *= package
        if versions != VersionSet()
            for ival in versions.intervals
                (content *= " $(ival.lower)")
                ival.upper < typemax(VersionNumber) &&
                (content *= " $(ival.upper)")
            end
        end
        new(content, package, versions, system)
    end
end

# TODO: shouldn't be neccessary #4648
==(a::Line, b::Line) = a.content == b.content

# general machinery for parsing REQUIRE files

function read(readable::Union(IO,Base.AbstractCmd))
    lines = Line[]
    for line in eachline(readable)
        line = chomp(line)
        push!(lines, ismatch(r"^\s*(?:#|$)", line) ? Comment(line) : Requirement(line))
    end
    return lines
end
read(file::String) = isfile(file) ? open(read,file) : Line[]

function write(io::IO, lines::Vector{Line})
    for line in lines
        println(io, line.content)
    end
end
function write(io::IO, reqs::Requires)
    for pkg in sort!([keys(reqs)...], by=lowercase)
        println(io, Requirement(pkg, reqs[pkg]).content)
    end
end
write(file::String, r::Union(Vector{Line},Requires)) = open(io->write(io,r), file, "w")

function parse(lines::Vector{Line})
    reqs = Requires()
    for line in lines
        if isa(line,Requirement)
            if !isempty(line.system)
                applies = false
                @windows_only applies |=  ("windows"  in line.system)
                @unix_only    applies |=  ("unix"     in line.system)
                @osx_only     applies |=  ("osx"      in line.system)
                @linux_only   applies |=  ("linux"    in line.system)
                @windows_only applies &= !("!windows" in line.system)
                @unix_only    applies &= !("!unix"    in line.system)
                @osx_only     applies &= !("!osx"     in line.system)
                @linux_only   applies &= !("!linux"   in line.system)
                applies || continue
            end
            reqs[line.package] = haskey(reqs, line.package) ?
                intersect(reqs[line.package], line.versions) : line.versions
        end
    end
    return reqs
end
parse(x) = parse(read(x))

# add & rm – edit the content a requires file

function add(lines::Vector{Line}, pkg::String, versions::VersionSet=VersionSet())
    v = VersionSet[]
    filtered = filter(lines) do line
        if !isa(line,Comment) && line.package == pkg && isempty(line.system)
            push!(v, line.versions)
            return false
        end
        return true
    end
    length(v) == 1 && v[1] == intersect(v[1],versions) && return copy(lines)
    versions = reduce(intersect, versions, v)
    push!(filtered, Requirement(pkg, versions))
end

rm(lines::Vector{Line}, pkg::String) = filter(lines) do line
    isa(line,Comment) || line.package != pkg
end

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demohelix()

Demonstration of a 3D helical curve.

# Examples:
```julia
julia> demohelix()
```
"""
function demohelix()
coefs =[ 6.0  0.0  6.0  1;
        -5.5  0.5  5.5  1;
        -5.0  1.0 -5.0  1;
         4.5  1.5 -4.5  1;
         4.0  2.0  4.0  1;
        -3.5  2.5  3.5  1;
        -3.0  3.0 -3.0  1;
         2.5  3.5 -2.5  1;
         2.0  4.0  2.0  1;
        -1.5  4.5  1.5  1;
        -1.0  5.0 -1.0  1;
         0.5  5.5 -0.5  1;
         0.0  6.0  0.0  1];
coefs = permutedims(coefs,[2,1]);
knots = vec([0 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 1 1]);

crv = nrbmak(coefs,knots);
nrbplot(crv,100,
        title = "3D helical curve.";
        legend = false);
end # demohelix


text:
# Reverse-time migration of the Sigsbee2A model: time-domain vs. frequency domain w/ on-the-fly Fourier transforms
# Author: Philipp Witte, pwitte.slim@gmail.com
# Date: May 2018
#

# TO DO:
# Replace w/ full path the observed data directory
path_to_data = "/path/to/directory/"
data_name = "sigsbee2A_marine"  # common base name of all shots

using JUDI.TimeModeling, PyPlot, JLD, SegyIO

# Load Sigsbee model
M = load("sigsbee2A_model.jld")

# Setup info and model structure
model0 = Model(M["n"], M["d"], M["o"], M["m0"])
dm = vec(M["dm"])

# Set up out-of-core data container
container = segy_scan(path_to_data, data_name, ["GroupX","GroupY","RecGroupElevation","SourceSurfaceElevation","dt"])
d_lin = judiVector(container)

# Set up source
src_geometry = Geometry(container; key="source")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.015)  # 15 Hz peak frequency
q = judiVector(src_geometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(q.geometry,d_lin.geometry, model0)
info = Info(prod(model0.n), d_lin.nsrc, ntComp)

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

opt = Options(isic=true, optimal_checkpointing=true)    # use impedance imaging

# Setup operators
Pr = judiProjection(info, d_lin.geometry)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, q.geometry)
J = judiJacobian(Pr*F0*adjoint(Ps), q)

# Time-domain RTM w/ optimal checkpointing
rtm_time = adjoint(J)*d_lin

# Save time-domain result
save("sigsbee2A_rtm_time_domain", "rtm", reshape(rtm_time, model0.n))

# Generate probability density function from source spectrum
q_dist = generate_distribution(q)

# Select 20 random frequencies per source location
nfreq = 20
J.options.optimal_checkpointing = false
J.options.frequencies = Array{Any}(undef, d_lin.nsrc)
for j=1:d_lin.nsrc
    J.options.frequencies[j] = select_frequencies(q_dist; fmin=0.003, fmax=0.04, nf=nfreq)
end
J.options.dft_subsampling_factor = 8    # perform dft every 8 time steps only

# Frequency-domain RTM w/ on-the-fly DFTs
rtm_freq = adjoint(J)*d_lin

# Save frequency-domain result
save("sigsbee2A_rtm_frequency_domain", "rtm", reshape(rtm_freq, model0.n))


text:
int main()
{
   int a[125];
   int i;
   for (i = 0; i < 125; i++)
   {
      a[i] = i*i;
   }
   for (i = 0; i < 125; i++)
   {
      a[i] = a[a[i] % 125] * 7;
   }
   printInt(a[0]);
   printInt(a[13]);
   printInt(a[124]);
   return 0;
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
# Conversions
# -----------

#=
`convert(C, c)` might be called as `convert(RGB, c)` or
`convert(RGB{Float32}, c)`.
This is handled in ColorTypes, which calls functions
```
    _convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c)
    _convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c, alpha)
```
Here are the argument types:
- `Cdest` may be any concrete `Colorant{T,N}` type. For parametric Color types
  it _always_ has the desired element type (e.g., `Float32`), so it's safe to
  dispatch on `Cdest <: Colorant{T}`.
- `Odest` and `Osrc` are opaque `Color` subtypes, i.e., things like `RGB` or
  `HSV`. They have no element type.
- `c` is the `Colorant` object you wish to convert.
- `alpha`, if present, is a user-supplied alpha value (to be used in place of
  any default alpha or alpha present in `c`).

The original motivation for this arrangement was that Julia "did not" support
"triangular dispatch", e.g.,
```
    convert(::Type{C{T}}, c) where {C, T}
```
On Julia v0.6 or later, parameter constraints can refer to previous parameters.
Threfore, we can use:
```
    convert(::Type{C}, c) where {T, C <: Colorant{T}}
```
However, the example above does not match `convert(RGB, c)`. Also, we should
catch all the various alpha variants (e.g. `ARGB`/`RGBA` with/without element
type).
The various arguments of `_convert` "peel off" element types (or guarantee them)
so that comparisons may be made via dispatch. Therefore, this arrangement is
still helpful.

Note that ColorTypes handles the cases where `Odest == Osrc`, or they are both
subtypes of `AbstractRGB` or `AbstractGray`. Therefore, here we only have to
deal with conversions between different spaces.
=#

function ColorTypes._convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, p, alpha=alpha(p)) where {Cdest<:TransparentColor,Odest,Osrc}
    # Convert the base color
    c = cnvt(color_type(Cdest), color(p))
    # Append the alpha
    ColorTypes._convert(Cdest, Odest, Odest, c, alpha)
end

ColorTypes._convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c) where {Cdest<:Color,Odest,Osrc} = cnvt(Cdest, c)

# with the whitepoint `wp` as the third argument
convert(::Type{XYZ}, c, wp::XYZ) = cnvt(XYZ{eltype(wp)}, c, wp)
convert(::Type{Lab}, c, wp::XYZ) = cnvt(Lab{eltype(wp)}, c, wp)
convert(::Type{Luv}, c, wp::XYZ) = cnvt(Luv{eltype(wp)}, c, wp)

convert(::Type{XYZ{T}}, c, wp::XYZ) where {T} = cnvt(XYZ{T}, c, wp)
convert(::Type{Lab{T}}, c, wp::XYZ) where {T} = cnvt(Lab{T}, c, wp)
convert(::Type{Luv{T}}, c, wp::XYZ) where {T} = cnvt(Luv{T}, c, wp)

# Fallback to catch undefined operations
cnvt(::Type{C}, c::TransparentColor) where {C<:Color} = cnvt(C, color(c))
cnvt(::Type{C}, c) where {C} = convert(C, c)

# Conversions from grayscale
# --------------------------
cnvt(::Type{C}, g::AbstractGray) where {C<:Color3}  = cnvt(C, convert(RGB{eltype(C)}, g))

macro mul3x3(T, M, c1, c2, c3)
    esc(quote
        @inbounds ret = $T($M[1,1]*$c1 + $M[1,2]*$c2 + $M[1,3]*$c3,
                           $M[2,1]*$c1 + $M[2,2]*$c2 + $M[2,3]*$c3,
                           $M[3,1]*$c1 + $M[3,2]*$c2 + $M[3,3]*$c3)
        ret
        end)
end

# Everything to RGB
# -----------------

correct_gamut(c::CV) where {CV<:AbstractRGB} = CV(clamp01(red(c)), clamp01(green(c)), clamp01(blue(c)))
correct_gamut(c::CV) where {T<:Union{N0f8,N0f16,N0f32,N0f64},
                            CV<:Union{AbstractRGB{T},TransparentRGB{T}}} = c
correct_gamut(c::CV) where {CV<:TransparentRGB} =
    CV(clamp01(red(c)), clamp01(green(c)), clamp01(blue(c)), clamp01(alpha(c))) # for `hex`

function srgb_compand(v::Fractional)
    # `pow5_12` is an optimized function to get `v^(1/2.4)`
    v <= 0.0031308 ? 12.92v : 1.055 * pow5_12(v) - 0.055
end

cnvt(::Type{CV}, c::AbstractRGB) where {CV<:AbstractRGB} = CV(red(c), green(c), blue(c))

function _hsx_to_rgb(im::UInt8, v, n, m)
    #=
    if     hue <  60; im = 0b000001 # ---------+
    elseif hue < 120; im = 0b000010 # --------+|
    elseif hue < 180; im = 0b000100 # -------+||
    elseif hue < 240; im = 0b001000 # ------+|||
    elseif hue < 300; im = 0b010000 # -----+||||
    else            ; im = 0b100000 # ----+|||||
    end                             #     ||||||
    (hue < 60 || hue >= 300) === ((im & 0b100001) != 0x0)
    =#
    r = ifelse((im & 0b100001) == 0x0, ifelse((im & 0b010010) == 0x0, m, n), v)
    g = ifelse((im & 0b000110) == 0x0, ifelse((im & 0b001001) == 0x0, m, n), v)
    b = ifelse((im & 0b011000) == 0x0, ifelse((im & 0b100100) == 0x0, m, n), v)
    return (r, g, b)
end
function _hsx_to_rgb(im::UInt8, v::T, n::T, m::T) where T <:Union{Float16, Float32, Float64}
    vu, nu, mu = reinterpret.(Unsigned, (v, n, m)) # prompt the compiler to use conditional moves
    r = ifelse((im & 0b100001) == 0x0, ifelse((im & 0b010010) == 0x0, mu, nu), vu)
    g = ifelse((im & 0b000110) == 0x0, ifelse((im & 0b001001) == 0x0, mu, nu), vu)
    b = ifelse((im & 0b011000) == 0x0, ifelse((im & 0b100100) == 0x0, mu, nu), vu)
    return reinterpret.(T, (r, g, b))
end

function cnvt(::Type{CV}, c::HSV) where {T, CV<:AbstractRGB{T}}
    F = promote_type(T, eltype(c))
    h, s, v = div60(F(c.h)), clamp01(F(c.s)), clamp01(F(c.v))
    hi = unsafe_trunc(Int32, h) # instead of floor
    i = h < 0 ? hi - one(hi) : hi
    f = i & one(i) == zero(i) ? 1 - (h - i) : h - i
    im = 0x1 << (mod6(UInt8, i) & 0x07)
    # use `@fastmath` just to reduce the estimated costs for inlining
    @fastmath m = v * (1 - s)
    @fastmath n = v * (1 - s * f)

    r, g, b = _hsx_to_rgb(im, v, n, m)
    T <: FixedPoint && typemax(T) >= 1 ? CV(r % T, g % T, b % T) : CV(r, g, b)
end

function cnvt(::Type{CV}, c::HSL) where {T, CV<:AbstractRGB{T}}
    F = promote_type(T, eltype(c))
    h, s, l = div60(F(c.h)), clamp01(F(c.s)), clamp01(F(c.l))
    a = @fastmath min(l, 1 - l) * s
    v = l + a
    hi = unsafe_trunc(Int32, h) # instead of floor
    i = h < 0 ? hi - one(hi) : hi
    f = i & one(i) == zero(i) ? 1 - (h - i) : h - i
    im = 0x1 << (mod6(UInt8, i) & 0x07)
    # use `@fastmath` just to reduce the estimated costs for inlining
    @fastmath m = l - a # v - 2 * a
    @fastmath n = v - 2 * a * f

    r, g, b = _hsx_to_rgb(im, v, n, m)
    T <: FixedPoint && typemax(T) >= 1 ? CV(r % T, g % T, b % T) : CV(r, g, b)
end

function cnvt(::Type{CV}, c::HSI) where {T, CV<:AbstractRGB{T}}
    F = promote_type(T, eltype(c))
    h, s, i = deg2rad(normalize_hue(F(c.h))), clamp01(F(c.s)), clamp01(F(c.i))
    is = i * s
    if h < F(2π/3)
        @fastmath cosr = cos(h) / cos(F(π/3)-h)
        r0, g0, b0 = muladd(is, cosr, i), muladd(is, 1-cosr, i), i - is
    elseif h < F(4π/3)
        @fastmath cosr = cos(h-F(2π/3)) / cos(F(π)-h)
        r0, g0, b0 = i - is, muladd(is, cosr, i), muladd(is, 1-cosr, i)
    else
        @fastmath cosr = cos(h-F(4π/3)) / cos(F(5π/3)-h)
        r0, g0, b0 = muladd(is, 1-cosr, i), i - is, muladd(is, cosr, i)
    end
    r, g, b = min(r0, oneunit(F)), min(g0, oneunit(F)), min(b0, oneunit(F))
    T <: FixedPoint && typemax(T) >= 1 ? CV(r % T, g % T, b % T) : CV(r, g, b)
end

function cnvt(::Type{CV}, c::XYZ) where CV<:AbstractRGB
    r =  3.2404542*c.x - 1.5371385*c.y - 0.4985314*c.z
    g = -0.9692660*c.x + 1.8760108*c.y + 0.0415560*c.z
    b =  0.0556434*c.x - 0.2040259*c.y + 1.0572252*c.z
    CV(clamp01(srgb_compand(r)),
       clamp01(srgb_compand(g)),
       clamp01(srgb_compand(b)))
end

function cnvt(::Type{CV}, c::YIQ) where CV<:AbstractRGB
    cc = correct_gamut(c)
    CV(clamp01(cc.y+0.9563*cc.i+0.6210*cc.q),
       clamp01(cc.y-0.2721*cc.i-0.6474*cc.q),
       clamp01(cc.y-1.1070*cc.i+1.7046*cc.q))
end

function cnvt(::Type{CV}, c::YCbCr) where CV<:AbstractRGB
    cc = correct_gamut(c)
    ny = cc.y - 16
    ncb = cc.cb - 128
    ncr = cc.cr - 128
    CV(clamp01(0.004567ny - 1.39135e-7ncb + 0.0062586ncr),
       clamp01(0.004567ny - 0.00153646ncb - 0.0031884ncr),
       clamp01(0.004567ny + 0.00791058ncb - 2.79201e-7ncr))
end

cnvt(::Type{CV}, c::Color3) where {CV<:AbstractRGB} = cnvt(CV, convert(XYZ{eltype(c)}, c))

# AbstractGray --> AbstractRGB conversions are implemented in ColorTypes.jl

# Everything to HSV
# -----------------

function cnvt(::Type{HSV{T}}, c::AbstractRGB) where T
    F = promote_type(T, eltype(c))
    r, g, b = F.((red(c), green(c), blue(c)))
    c_min = @fastmath min(min(r, g), b)
    c_max = @fastmath max(max(r, g), b)
    s0 = c_max - c_min
    s0 == zero(F) && return HSV{T}(zero(T), zero(T), T(c_max))
    s = @fastmath s0 / c_max

    # In general, it is dangerous to compare floating point numbers with `===`.
    diff = ifelse(c_max === r,  g - b,         ifelse(c_max === g,  b - r,  r - g))
    ofs  = ifelse(c_max === r, (g < b)*F(360), ifelse(c_max === g, F(120), F(240)))
    h0 = @fastmath diff * F(60) / s0

    HSV{T}(h0 + ofs, s, c_max)
end

cnvt(::Type{HSV{T}}, c::Color3) where {T} = cnvt(HSV{T}, convert(RGB{T}, c))

# Everything to HSL
# -----------------

function cnvt(::Type{HSL{T}}, c::AbstractRGB) where T
    F = promote_type(T, eltype(c))
    r, g, b = F(red(c)), F(green(c)), F(blue(c))
    c_min = @fastmath min(min(r, g), b)
    c_max = @fastmath max(max(r, g), b)
    l0 = c_max + c_min
    s0 = c_max - c_min
    l = l0 * F(0.5)
    s0 == zero(F) && return HSL{T}(zero(T), zero(T), T(l))
    s = @fastmath s0 / min(l0, F(2) - l0)

    # In general, it is dangerous to compare floating point numbers with `===`.
    diff = ifelse(c_max === r,  g - b,         ifelse(c_max === g,  b - r,  r - g))
    ofs  = ifelse(c_max === r, (g < b)*F(360), ifelse(c_max === g, F(120), F(240)))
    h0 = @fastmath diff * F(60) / s0

    HSL{T}(h0 + ofs, s, l)
end

cnvt(::Type{HSL{T}}, c::Color3) where {T} = cnvt(HSL{T}, convert(RGB{T}, c))

# Everything to HSI
# -----------------

# Since acosd() is slow, the following is "inline-worthy".
@inline function cnvt(::Type{HSI{T}}, c::AbstractRGB) where T
    rgb = correct_gamut(c)
    F = promote_type(T, eltype(c))
    r, g, b = F(red(rgb)), F(green(rgb)), F(blue(rgb))
    dnorm = @fastmath sqrt(((r-g)^2 + (r-b)^2 + (g-b)^2) * F(0.5))
    isum = r + g + b
    i = isum / F(3)
    dnorm == zero(F) && return HSI{T}(T(90), zero(T), T(i))
    val = muladd(g + b, F(-0.5), r) / dnorm
    h = @fastmath acosd(clamp(val, -oneunit(F), oneunit(F)))
    m = @fastmath min(min(r, g), b)
    s = oneunit(F) - m/i
    HSI{T}(b > g ? F(360) - h : h, s, i)
end

cnvt(::Type{HSI{T}}, c::Color3) where {T} = cnvt(HSI{T}, convert(RGB{T}, c))

# Everything to XYZ
# -----------------

function invert_srgb_compand(v::Fractional)
    # `pow12_5` is an optimized function to get `x^2.4`
    v <= 0.04045 ? 1/12.92 * v : pow12_5(1/1.055 * (v + 0.055))
end

# lookup table for `N0f8` (the extra two elements are for `Float32` splines)
const invert_srgb_compand_n0f8 = [invert_srgb_compand(v/255.0) for v = 0:257]

function invert_srgb_compand(v::N0f8)
    @inbounds invert_srgb_compand_n0f8[reinterpret(UInt8, v) + 1]
end

function invert_srgb_compand(v::Float32)
    i = unsafe_trunc(Int32, v * 255)
    (i < 13 || i > 255) && return invert_srgb_compand(Float64(v))
    @inbounds y = view(invert_srgb_compand_n0f8, i:i+3)
    dv = v * 255.0 - i
    dv == 0.0 && @inbounds return y[2]
    if v < 0.38857287f0
        return @fastmath(y[2]+0.5*dv*((-2/3*y[1]- y[2])+(2y[3]-1/3*y[4])+
                                  dv*((     y[1]-2y[2])+  y[3]-
                                  dv*(( 1/3*y[1]- y[2])+( y[3]-1/3*y[4]) ))))
    else
        return @fastmath(y[2]+0.5*dv*((4y[3]-3y[2])-y[4]+dv*((y[4]-y[3])+(y[2]-y[3]))))
    end
end

function cnvt(::Type{XYZ{T}}, c::AbstractRGB) where T
    r, g, b = invert_srgb_compand(red(c)), invert_srgb_compand(green(c)), invert_srgb_compand(blue(c))
    XYZ{T}(0.4124564*r + 0.3575761*g + 0.1804375*b,
           0.2126729*r + 0.7151522*g + 0.0721750*b,
           0.0193339*r + 0.1191920*g + 0.9503041*b)
end

function cnvt(::Type{XYZ{T}}, c::xyY) where T
    X = c.Y*c.x/c.y
    Z = c.Y*(1-c.x-c.y)/c.y
    XYZ{T}(X, c.Y, Z)
end

const xyz_epsilon = 216 / 24389
const xyz_kappa   = 24389 / 27

function cnvt(::Type{XYZ{T}}, c::Lab, wp::XYZ = WP_DEFAULT) where T
    fy = (c.l + 16) / 116
    fx = c.a / 500 + fy
    fz = fy - c.b / 200

    fx3 = fx^3
    fz3 = fz^3

    x = fx3 > xyz_epsilon ? fx3 : (116fx - 16) / xyz_kappa
    y = c.l > xyz_kappa * xyz_epsilon ? ((c. l+ 16) / 116)^3 : c.l / xyz_kappa
    z = fz3 > xyz_epsilon ? fz3 : (116fz - 16) / xyz_kappa

    XYZ{T}(x*wp.x, y*wp.y, z*wp.z)
end

function xyz_to_uv(c::XYZ)
    d = c.x + 15c.y + 3c.z
    d==0 && return (d, d)
    u = 4c.x / d
    v = 9c.y / d
    return (u, v)
end

function cnvt(::Type{XYZ{T}}, c::Luv, wp::XYZ = WP_DEFAULT) where T
    (u_wp, v_wp) = xyz_to_uv(wp)

    a = (52 * (c.l==0 ? zero(T) : c.l / (c.u + 13 * c.l * u_wp)) - 1) / 3
    y = c.l > xyz_kappa * xyz_epsilon ? wp.y * ((c.l + 16) / 116)^3 : wp.y * c.l / xyz_kappa
    b = -5y
    d = y * (39 * (c.l==0 ? zero(T) : c.l / (c.v + 13 * c.l * v_wp)) - 5)
    x = d==b ? zero(T) : (d - b) / (a + 1/3)
    z = a * x + b + zero(T)

    XYZ{T}(x, y, z)
end

function cnvt(::Type{XYZ{T}}, c::DIN99d) where T

    # Go back to C-h space
    C = chroma(c)
    h = atan(c.b, c.a) - 50π/180

    # Intermediate terms
    G = (exp(C/22.5)-1)/0.06
    f, ee = G .* sincos(h)

    l = (exp(c.l/325.221)-1)/0.0036
    # a = ee*cosd(50) - f/1.14*sind(50)
    a = ee*0.6427876096865394 - f/1.14*0.766044443118978
    # b = ee*sind(50) - f/1.14*cosd(50)
    b = ee*0.766044443118978 + f/1.14*0.6427876096865394

    adj = convert(XYZ, Lab(l, a, b))

    XYZ{T}((adj.x + 0.12*adj.z)/1.12, adj.y, adj.z)

end

function cnvt(::Type{XYZ{T}}, c::LMS) where T
    @mul3x3 XYZ{T} CAT02_INV c.l c.m c.s
end

cnvt(::Type{XYZ{T}}, c::Union{LCHab, DIN99, DIN99o}) where {T} = cnvt(XYZ{T}, convert(Lab{T}, c))
cnvt(::Type{XYZ{T}}, c::LCHuv) where {T} = cnvt(XYZ{T}, convert(Luv{T}, c))
cnvt(::Type{XYZ{T}}, c::Color3) where {T} = cnvt(XYZ{T}, convert(RGB{T}, c))

# Everything to xyY
# -----------------

function cnvt(::Type{xyY{T}}, c::XYZ) where T

    x = c.x/(c.x + c.y + c.z)
    y = c.y/(c.x + c.y + c.z)

    xyY{T}(x, y, convert(typeof(x), c.y))

end

cnvt(::Type{xyY{T}}, c::Color3) where {T} = cnvt(xyY{T}, convert(XYZ{T}, c))

# Everything to Lab
# -----------------

function fxyz2lab(v)
    v > xyz_epsilon ? cbrt(v) : (xyz_kappa * v + 16) / 116
end
function cnvt(::Type{Lab{T}}, c::XYZ, wp::XYZ = WP_DEFAULT) where T
    fx, fy, fz = fxyz2lab(c.x / wp.x), fxyz2lab(c.y / wp.y), fxyz2lab(c.z / wp.z)
    Lab{T}(116fy - 16, 500(fx - fy), 200(fy - fz))
end

function cnvt(::Type{Lab{T}}, c::LCHab) where T
    b, a = c.c .* sincos(deg2rad(c.h))
    Lab{T}(c.l, a, b)
end

function cnvt(::Type{Lab{T}}, c::DIN99) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate Chroma (C99) in the DIN99 space
    cc = chroma(c)

    # Temporary variable for chroma
    g = (exp(0.045*cc*kch*ke)-1)/0.045

    # Temporary redness
    ee = cc > 0 ? g * c.a / cc : zero(g)

    # Temporary yellowness
    f = cc > 0 ? g * c.b / cc : zero(g)

    # CIELAB a*b*
    # ciea = ee*cosd(16) - (f/0.7)*sind(16)
    ciea = ee*0.9612616959383189 - (f/0.7)*0.27563735581699916
    # cieb = ee*sind(16) + (f/0.7)*cosd(16)
    cieb = ee*0.27563735581699916 + (f/0.7)*0.9612616959383189

    # CIELAB L*
    ciel = (exp(c.l*ke/105.51)-1)/0.0158

    Lab{T}(ciel, ciea, cieb)
end

function cnvt(::Type{Lab{T}}, c::DIN99o) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate Chroma (C99) in the DIN99o space
    co = chroma(c)

    # hue angle h99o
    h = atan(c.b, c.a) - 26π/180

    # revert logarithmic chroma compression
    g = (exp(co*kch*ke/23.0)-1)/0.075

    # Temporary yellowness and redness
    fo, eo = g .* sincos(h)

    # CIELAB a*b* (revert b* axis compression)
    # ciea = eo*cosd(26) - (fo/0.83)*sind(26)
    ciea = eo*0.898794046299167 - (fo/0.83)*0.4383711467890774
    # cieb = eo*sind(26) + (fo/0.83)*cosd(26)
    cieb = eo*0.4383711467890774 + (fo/0.83)*0.898794046299167

    # CIELAB L* (revert logarithmic lightness compression)
    ciel = (exp(c.l*ke/303.67)-1)/0.0039

    Lab{T}(ciel, ciea, cieb)
end

cnvt(::Type{Lab{T}}, c::Color3) where {T} = cnvt(Lab{T}, convert(XYZ{T}, c))

# Everything to Luv
# -----------------

function cnvt(::Type{Luv{T}}, c::XYZ, wp::XYZ = WP_DEFAULT) where T
    (u_wp, v_wp) = xyz_to_uv(wp)
    (u_, v_) = xyz_to_uv(c)

    y = c.y / wp.y

    l = y > xyz_epsilon ? 116 * cbrt(y) - 16 : xyz_kappa * y
    u = 13 * l * (u_ - u_wp) + zero(T)
    v = 13 * l * (v_ - v_wp) + zero(T)

    Luv{T}(l, u, v)
end

function cnvt(::Type{Luv{T}}, c::LCHuv) where T
    v, u = c.c .* sincos(deg2rad(c.h))
    Luv{T}(c.l, u, v)
end

cnvt(::Type{Luv{T}}, c::Color3) where {T} = cnvt(Luv{T}, convert(XYZ{T}, c))

# Everything to LCHuv
# -------------------

function cnvt(::Type{LCHuv{T}}, c::Luv) where T
    LCHuv{T}(c.l, chroma(c), hue(c))
end

cnvt(::Type{LCHuv{T}}, c::Color3) where {T} = cnvt(LCHuv{T}, convert(Luv{T}, c))

# Everything to LCHab
# -------------------

function cnvt(::Type{LCHab{T}}, c::Lab) where T
    LCHab{T}(c.l, chroma(c), hue(c))
end

cnvt(::Type{LCHab{T}}, c::Color3) where {T} = cnvt(LCHab{T}, convert(Lab{T}, c))

# Everything to DIN99
# -------------------

function cnvt(::Type{DIN99{T}}, c::Lab) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate DIN99 L
    l99 = (1/ke)*105.51*log(1+0.0158*c.l)

    # Temporary value for redness and yellowness
    # ee = c.a*cosd(16) + c.b*sind(16)
    ee = c.a*0.9612616959383189 + c.b*0.27563735581699916
    # f = -0.7*c.a*sind(16) + 0.7*c.b*cosd(16)
    f = -0.7*c.a*0.27563735581699916 + 0.7*c.b*0.9612616959383189

    # Temporary value for chroma
    g = sqrt(ee^2 + f^2)

    # DIN99 chroma
    cc = log(1+0.045*g)/(0.045*kch*ke)

    # DIN99 chromaticities
    a99 = g > 0 ? cc * ee / g : zero(T)
    b99 = g > 0 ? cc * f / g : zero(T)

    DIN99{T}(l99, a99, b99)

end

cnvt(::Type{DIN99{T}}, c::Color3) where {T} = cnvt(DIN99{T}, convert(Lab{T}, c))

# Everything to DIN99d
# --------------------

function cnvt(::Type{DIN99d{T}}, c::XYZ{T}) where T

    # Apply tristimulus-space correction term
    adj_c = XYZ(1.12*c.x - 0.12*c.z, c.y, c.z)

    # Apply L*a*b*-space correction
    lab = convert(Lab, adj_c)
    adj_L = 325.221*log(1+0.0036*lab.l)

    # Calculate intermediate parameters
    # ee = lab.a*cosd(50) + lab.b*sind(50)
    ee = lab.a*0.6427876096865394 + lab.b*0.766044443118978
    # f = 1.14*(lab.b*cosd(50) - lab.a*sind(50))
    f = 1.14*(lab.b*0.6427876096865394 - lab.a*0.766044443118978)
    G = sqrt(ee^2+f^2)

    # Calculate hue/chroma
    C = 22.5*log(1+0.06*G)
    h = atan(f, ee) + 50π/180

    b99, a99 = C .* sincos(h)

    DIN99d{T}(adj_L, a99, b99)

end

cnvt(::Type{DIN99d{T}}, c::Color3) where {T} = cnvt(DIN99d{T}, convert(XYZ{T}, c))

# Everything to DIN99o
# -------------------

function cnvt(::Type{DIN99o{T}}, c::Lab) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate DIN99o L (logarithmic compression)
    l99 = 303.67/ke*log(1+0.0039*c.l)

    # Temporary value for redness and yellowness
    # including rotation by 26°
    # eo = c.a*cosd(26) + c.b*sind(26)
    eo = c.a*0.898794046299167 + c.b*0.4383711467890774
    # compression along the yellowness (blue-yellow) axis
    # fo = 0.83 * (c.b*cosd(26) - c.a*sind(26))
    fo = 0.83 * (c.b*0.898794046299167 - c.a*0.4383711467890774)

    # Temporary value for chroma
    go = sqrt(eo^2 + fo^2)
    h = atan(fo, eo) + 26π/180

    # DIN99o chroma (logarithmic compression)
    cc = 23.0*log(1+0.075*go)/(kch*ke)

    # DIN99o chromaticities
    b99, a99 = cc .* sincos(h)

    DIN99o{T}(l99, a99, b99)

end

cnvt(::Type{DIN99o{T}}, c::Color3) where {T} = cnvt(DIN99o{T}, convert(Lab{T}, c))

# Everything to LMS
# -----------------

# Chromatic adaptation from CIECAM97s
const CAT97s = [ 0.8562  0.3372 -0.1934
                -0.8360  1.8327  0.0033
                 0.0357 -0.0469  1.0112 ]

const CAT97s_INV = inv(CAT97s)

# Chromatic adaptation from CIECAM02
const CAT02 = [ 0.7328 0.4296 -0.1624
               -0.7036 1.6975  0.0061
                0.0030 0.0136  0.9834 ]

const CAT02_INV = inv(CAT02)

function cnvt(::Type{LMS{T}}, c::XYZ) where T
    @mul3x3 LMS{T} CAT02 c.x c.y c.z
end

cnvt(::Type{LMS{T}}, c::Color3) where {T} = cnvt(LMS{T}, convert(XYZ{T}, c))

# Everything to YIQ
# -----------------

correct_gamut(c::YIQ{T}) where {T} = YIQ{T}(clamp(c.y, zero(T), one(T)),
                                     clamp(c.i, convert(T,-0.5957), convert(T,0.5957)),
                                     clamp(c.q, convert(T,-0.5226), convert(T,0.5226)))

function cnvt(::Type{YIQ{T}}, c::AbstractRGB) where T
    rgb = correct_gamut(c)
    YIQ{T}(0.299*red(rgb)+0.587*green(rgb)+0.114*blue(rgb),
           0.595716*red(rgb)-0.274453*green(rgb)-0.321263*blue(rgb),
           0.211456*red(rgb)-0.522591*green(rgb)+0.311135*blue(rgb))
end

cnvt(::Type{YIQ{T}}, c::Color3) where {T} = cnvt(YIQ{T}, convert(RGB{T}, c))

# Everything to YCbCr
# -------------------

correct_gamut(c::YCbCr{T}) where {T} = YCbCr{T}(clamp(c.y, convert(T,16), convert(T,235)),
                                         clamp(c.cb, convert(T,16), convert(T,240)),
                                         clamp(c.cr, convert(T,16), convert(T,240)))

function cnvt(::Type{YCbCr{T}}, c::AbstractRGB) where T
    rgb = correct_gamut(c)
    YCbCr{T}(16+65.481*red(rgb)+128.553*green(rgb)+24.966*blue(rgb),
             128-37.797*red(rgb)-74.203*green(rgb)+112*blue(rgb),
             128+112*red(rgb)-93.786*green(rgb)-18.214*blue(rgb))
end

cnvt(::Type{YCbCr{T}}, c::Color3) where {T} = cnvt(YCbCr{T}, convert(RGB{T}, c))

# To Gray
# -------
# AbstractGray --> AbstractRGB conversions are implemented in ColorTypes.jl, but
# AbstractRGB --> AbstractGray conversions should be implemented here.

# Rec 601 luma conversion

function cnvt(::Type{G}, x::AbstractRGB{T}) where {G<:AbstractGray,T<:Normed}
    TU, Tf = FixedPointNumbers.rawtype(T), floattype(T)
    if sizeof(TU) < sizeof(UInt)
        val = Tf(0.001)*(299*reinterpret(red(x)) + 587*reinterpret(green(x)) + 114*reinterpret(blue(x)))
    else
        val = Tf(0.299)*reinterpret(red(x)) + Tf(0.587)*reinterpret(green(x)) + Tf(0.114)*reinterpret(blue(x))
    end
    return G(reinterpret(T, round(TU, val)))
end
cnvt(::Type{G}, x::AbstractRGB) where {G<:AbstractGray} =
    G(0.299f0*red(x) + 0.587f0*green(x) + 0.114f0*blue(x))

cnvt(::Type{G}, x::Color) where {G<:AbstractGray} = convert(G, convert(RGB, x))


text:
struct Node
    X::Float64
    Y::Float64
end

struct Support
    Node::Int
    XRelease::Int
    YRelease::Int
end

struct Material
    E::Float64
end

struct Section
    Area::Float64
end

struct Element
    NodeB::Int
    NodeE::Int
    Material
    Section
end

struct Load
    Node::Int
    Px::Float64
    Py::Float64
end

function gK(Xb, Yb, Xe, Ye, E, A)
  δx = Xe - Xb
  δy = Ye - Yb
  L = sqrt(δx^2 + δy^2)
  _cos = δx/L
  _sin = δy/L
  _cos² = _cos * _cos
  _sin² = _sin * _sin
  _cos_sin = _cos * _sin
  (A*E/L)*[_cos² _cos_sin -_cos² -_cos_sin
   _cos_sin _sin² -_cos_sin -_sin²
   -_cos² -_cos_sin _cos² _cos_sin
   -_cos_sin -_sin² _cos_sin _sin²]
end

function T(Xb, Yb, Xe, Ye)
    δx = Xe - Xb
    δy = Ye - Yb
    L = sqrt(δx^2 + δy^2)
    _cos = δx/L
    _sin = δy/L

    [_cos _sin 0 0
     -_sin _cos 0 0
     0 0 _cos _sin
     0 0 -_sin _cos]
end

function k(Xb, Yb, Xe, Ye, E, A)
    δx = Xe - Xb
    δy = Ye - Yb
    L = sqrt(δx^2 + δy^2)

    (E*A/L)*[1. 0 -1 0
            0 0 0 0
            -1 0 1 0
            0 0 0 0]
end

COORD = [Node(0, 0)
         Node(288, 0)
         Node(576, 0)
         Node(864, 0)
         Node(288, 216)
         Node(576, 216)]

MSUP = [Support(1, 1, 1)
            Support(3, 0, 1)
            Support(4, 0, 1)]

EM = Dict( 1 => Material(29000),
           2 => Material(10000))

CP = Dict( 1 => Section(8),
           2 => Section(12),
           3 => Section(16))

MPRP = [Element(1, 2, 1, 1)
        Element(2, 3, 1, 1)
        Element(3, 4, 2, 3)
        Element(5, 6, 1, 1)
        Element(2, 5, 1, 1)
        Element(3, 6, 1, 1)
        Element(1, 5, 1, 2)
        Element(2, 6, 1, 2)
        Element(3, 5, 1, 2)
        Element(4, 6, 2, 3)]

loads = Dict(2 => (0.0, -75.),
             5 => (25., 0.),
             6 => (0., -60.))

#=
COORD = [Node(12*12, 16*12)
         Node(0, 0)
         Node(12*12, 0)
         Node(24*12, 0)]

MSUP = [Support(2, 1, 1)
        Support(3, 1, 1)
        Support(4, 1, 1)]

EM = Dict( 1 => Material(29000))

CP = Dict( 1 => Section(8),
           2 => Section(6))

MPRP = [Element(2, 1, 1, 1)
        Element(3, 1, 1, 2)
        Element(4, 1, 1, 1)]

loads = Dict(1 => (150, -300))
=#
#=
COORD = [Node(0,0)
         Node(10,0)
         Node(0,8)
         Node(6,8)]

MSUP = [Support(1, 1, 1)
        Support(2, 1, 1)
        Support(3, 1, 0)]

EM = Dict( 1 => Material(70e6))

CP = Dict( 1 => Section(.004))

MPRP = [Element(1, 3, 1, 1)
        Element(3, 4, 1, 1)
        Element(1, 4, 1, 1)
        Element(2, 3, 1, 1)
        Element(2, 4, 1, 1)]

loads = Dict(3 => (0.0, -400.0),
             4 => (800.0, -400.0))
=#

function nsc(COORD, MSUP, NCJT)
    #Determine number of joints (NJ)
    NJ = size(COORD, 1)

    #Determine Number of Restraints (NR)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))

    numDOF = NJ * NCJT - NR

    NSC = Array{Int}(2 * NJ)                                  # Create the NSC with uninitialized values

    freeCoord = 1
    restraintCoord = numDOF + 1
    for nodeIndex in 1:size(COORD, 1)                           # Loop through all of the COORD
        supportIndexes = [s.Node for s in MSUP]
        if in(nodeIndex, supportIndexes)                        #   if this node is a support then...
            supportIndex = findfirst(supportIndexes, nodeIndex) #     get the support Index for this node
            if MSUP[supportIndex].XRelease == 1             #     if node is restrained in the X direction...
                 NSC[2*nodeIndex - 1] = restraintCoord           #       NSC X value for this node is the next restraint coordinate
                restraintCoord += 1
            else                                                #     otherwise, this node is free in the X direction
                 NSC[2*nodeIndex - 1] = freeCoord         #       NSC X value for this node is the next free coordinate
                freeCoord += 1
            end
            if MSUP[supportIndex].YRelease == 1             #     if node is restrained in the Y direction...
                NSC[2*nodeIndex] = restraintCoord         #       NSC Y value for this node is the next restraint coordinate
                restraintCoord += 1
            else                                                #     otherwise, this node is free in the Y direction
                 NSC[2*nodeIndex] = freeCoord             #       NSC Y value for this node is the next free coordinate
                freeCoord += 1
            end
        else                                                    #   otherwise, this node is free
             NSC[2*nodeIndex - 1] = freeCoord             #     NSC X value for this node is the next free coordinate
             NSC[2*nodeIndex] = freeCoord + 1             #     NSC Y value for this node is the next free coordinate
            freeCoord += 2
        end
    end
    NSC
end

#NCJT = 2
#NSC = nsc(COORD, MSUP, NCJT)

function gs(COORD, MSUP, MPRP, EM, CP, NSC, NCJT)
    #Determine number of joints (NJ)
    NJ = size(COORD, 1)

    #Determine Number of Restraints (NR)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))

    numDOF = NJ * NCJT - NR

    S = zeros(numDOF, numDOF)

    for memberIndex in 1:size(MPRP, 1)
        memberInfo = MPRP[memberIndex]
        nodeBegin, nodeEnd = memberInfo.NodeB, memberInfo.NodeE
        Xb = COORD[nodeBegin].X
        Yb = COORD[nodeBegin].Y
        Xe = COORD[nodeEnd].X
        Ye = COORD[nodeEnd].Y
        E = EM[memberInfo.Material].E
        A = CP[memberInfo.Section].Area
        GK = gK(Xb, Yb, Xe, Ye, A, E)

        NSC_Indecies =
        vcat([NSC[(nodeBegin - 1)*NCJT + x] for x in 1:NCJT],
             [NSC[(nodeEnd - 1)*NCJT + x] for x in 1:NCJT])

        iCount = 0
        for i in NSC_Indecies
            if i == NSC_Indecies[1]
                iCount = 0
            end
            iCount += 1
            if i <= numDOF
                jCount = 0
                for j in NSC_Indecies
                    jCount +=1
                    if j <= numDOF
                         S[j,i] = S[j,i] + GK[jCount, iCount]
                    end
                end
            end
        end
    end
    S
end
#S = gs(COORD, MSUP, MPRP, EM, CP, NSC, NCJT)

function p(COORD, MSUP, NSC, loads, NCJT)
    NJ = size(COORD, 1)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))
    numDOF = NJ * NCJT - NR
    NJL = length(loads)

    P = zeros(numDOF)::Array{Float64,1}

    dofIndecies = find(x -> x <= numDOF, NSC)

    for i in 1:size(dofIndecies,1)
        dofIndex = dofIndecies[i]
        remainder = dofIndex%NCJT
        dofOfNode = (remainder == 0 ? NCJT : remainder)
        node = Int(1 + (dofIndex - dofOfNode)/NCJT)
        if haskey(loads, node)
            P[i] = loads[node][dofOfNode]
        else
            P[i] = 0.0
        end
    end
    P
end

#P = p(COORD, MSUP, NSC, loads, NCJT)

function get_d(S, P)
    inv(S)*P
end

#d = get_d(S, P)

function forces_reactions(COORD, NSC, MSUP, EM, CP, MPRP, d, NCJT)
    NJ = size(COORD, 1)
    NR = sum(map(s -> s.XRelease + s.YRelease, MSUP))
    numDOF = NJ * NCJT - NR

    R = zeros(NR)
    forces = Array{Float64}(size(MPRP,1))

    for member in 1:size(MPRP,1)
        memberInfo = MPRP[member]
        nodeBegin, nodeEnd = memberInfo.NodeB, memberInfo.NodeE
        Xb = COORD[nodeBegin].X
        Yb = COORD[nodeBegin].Y
        Xe = COORD[nodeEnd].X
        Ye = COORD[nodeEnd].Y
        E = EM[memberInfo.Material].E
        A = CP[memberInfo.Section].Area

        v = zeros(NCJT*2)

        NSC_Indecies =
        vcat([NSC[(nodeBegin - 1)*NCJT + x] for x in 1:NCJT],
             [NSC[(nodeEnd - 1)*NCJT + x] for x in 1:NCJT])

        for i in 1:(NCJT*2)
            NSC_Index = NSC_Indecies[i]
            if NSC_Index <= numDOF
                v[i] = d[NSC_Index]
            end
        end

        _T = T(Xb, Yb, Xe, Ye)

        U = _T * v

        BK = k(Xb, Yb, Xe, Ye, E, A)

        Q = BK * U
        forces[member] = Q[1]

        F = _T'Q

        for i in 1:(NCJT*2)
            nscIndex =
                if i<=NCJT
                    (nodeBegin -1)*NCJT + i
                else
                    (nodeEnd - 1)*NCJT + (i - NCJT)
                end
            N = NSC[nscIndex]
            if N > numDOF
                R[N - numDOF] = R[N - numDOF] + F[i]
            end
        end
    end
    forces, R
end

#forces, reactions = forces_reactions(COORD, MSUP, EM, CP, MPRP, d, NCJT)

function run_analysis(COORD, MSUP, EM, CP, MPRP, loads)
    NCJT = 2
    NSC = nsc(COORD, MSUP, NCJT)
    GS = gs(COORD, MSUP, MPRP, EM, CP, NSC, NCJT)
    P = p(COORD, MSUP, NSC, loads, NCJT)
    d = get_d(GS, P)
    forces, reactions = forces_reactions(COORD, NSC, MSUP, EM, CP, MPRP, d, NCJT)
end

forces, reactions = run_analysis(COORD, MSUP, EM, CP, MPRP, loads)
println(forces)
println(reactions)

#=
NCJT = 2
NSC = nsc(COORD, MSUP, 2)
GS = gs(COORD, MSUP, MPRP, EM, CP, NSC, NCJT)
P = p(COORD, MSUP, NSC, loads, NCJT)
d = get_d(GS, P)

@time [nsc(COORD, MSUP, 2) for i in 1:10000]
@time [gs(COORD, MSUP, MPRP, EM, CP, NSC, NCJT) for i in 1:10000]
@time [p(COORD, MSUP, NSC, loads, NCJT) for i in 1:10000]
@time [get_d(GS, P) for i in 1:10000]
@time [forces_reactions(COORD, NSC, MSUP, EM, CP, MPRP, d, NCJT) for i in 1:10000]
@time [run_analysis(COORD, MSUP, EM, CP, MPRP, loads) for i in 1:10000]
=#

max_stars_repo_name: Julia-SE/matrix-analysis

max_issues_repo_name: Julia-SE/matrix-analysis

max_forks_repo_name: Julia-SE/matrix-analysis


text:
filePath = @__DIR__
(@isdefined(SphericalHarmonics_Mod)) ? nothing : include(string(filePath,"\\SphericalHarmonics.jl"))
using .SphericalHarmonics_Mod
using Plots
using Interpolations

savepath = string(filePath,"\\anim.gif")

function main(savepath)

    #pyplot is very slow, use gr for faster (but not as pretty plots)
    pyplot(size=(1080,1080))

    # l-values to be plotted
    lmax = 20
    sz = sum((lmax:lmax).+1) #Number of m-values

    N = 500; #Meshgrid in plot is NxN

    x = Array{Float64}(undef,N,N,sz)
    y = Array{Float64}(undef,N,N,sz)
    z = Array{Float64}(undef,N,N,sz)
    Ylm = Array{Float64}(undef,N,N,sz)

    println("Calculating Ylm")

    idx = 1
    for l in lmax:lmax, m = 0:l
        x[:,:,idx],y[:,:,idx],z[:,:,idx],Ylm[:,:,idx] = Calc_SphericalHarmonics(l,m,N)
        idx += 1
    end

    sz = idx-1
    x = x[:,:,begin:sz]
    y = y[:,:,begin:sz]
    z = z[:,:,begin:sz]
    Ylm = Ylm[:,:,begin:sz]

    interpVals = 49 #Number of values used for interpolation between each m-value. Use 49 to get 1 m-value plotted per 2 seconds in the gif.

    X = Array{Float64}(undef,N,N,sz*interpVals)
    Y = Array{Float64}(undef,N,N,sz*interpVals)
    Z = Array{Float64}(undef,N,N,sz*interpVals)
    YLM = Array{Float64}(undef,N,N,sz*interpVals)

    println("Interpolating")
    for i = 1:N, j = 1:N

        xs = 1:sz
        xintp = range(1,sz,sz*interpVals)

        interp  = LinearInterpolation(xs,x[i,j,:])
        X[i,j,:] .= interp.(xintp)

        interp  = LinearInterpolation(xs,y[i,j,:])
        Y[i,j,:] .= interp.(xintp)

        interp  = LinearInterpolation(xs,z[i,j,:])
        Z[i,j,:] .= interp.(xintp)

        interp  = LinearInterpolation(xs,Ylm[i,j,:])
        YLM[i,j,:] .= interp.(xintp)

    end

    println("Animating")
    anim = @animate for i = 1:sz*interpVals
        println(i,"/",sz*interpVals)

        surface(X[:,:,i],Y[:,:,i],Z[:,:,i], legend = nothing, fill_z = abs.(YLM[:,:,i]), axis=([], false), show = false)

    end

    println("Saving gif")
    gif(anim,savepath,fps=24)

    return

end

main(savepath)

max_stars_repo_name: philipbergander/SphericalHarmonics

max_issues_repo_name: philipbergander/SphericalHarmonics

max_forks_repo_name: philipbergander/SphericalHarmonics