text:
abstract type ApproxCDF end
abstract type InterpolateCDF <: ApproxCDF end
abstract type SmoothCDF <: ApproxCDF end

struct weights_values <: AbstractArray{Float64,1}
  weights::Vector{Float64}
  values::Vector{Float64}
end
struct Grid <: InterpolateCDF
  weights::Array{Float64,1}
  values::Array{Float64,1}
end
struct NestedPolyGLM{D <: ContinuousUnivariateDistribution} <: SmoothCDF
  β::Vector{Float64}
  θ::Vector{Float64}
  d::D
end

function simultaneous_sort!(wv::weights_values)
  si = sortperm(wv.values)
  wv.values .= wv.values[si]
  wv.weights .= wv.weights[si]
  wv
end

function interpolate_weight_values(wv::weights_values)
  simultaneous_sort!(wv)
  interpolate((wv.values,), cumsum(wv.weights), Gridded(Linear()))
end

OneParamReal = Union{TDist, VonMises}
OneParamRealScale = Union{Cauchy, Gumbel, Laplace, Logistic, Normal}
#ThreeParamReal = Union{NormalInverseGaussian}#Skip
#TwoParamProb = Union{Beta}
OneParamPositive = Union{Chi, Chisq, Exponential, Rayleigh}
OneParamPositiveScale = Union{InverseGaussian, Levy, LogNormal}
TwoParamPositive = Union{BetaPrime, FDist, Frechet, Gamma, InverseGamma, Pareto, Weibull}
#Do GeneralizedExtremeValue seperately, because it needs 3 params?
#Do Arcsine, SymTriangularDist, TriangularDist, Uniform seperately

RealDist = Union{OneParamReal, OneParamRealScale}
PositiveDist = Union{OneParamPositive, OneParamPositiveScale, TwoParamPositive}

OneParam = Union{OneParamReal, OneParamRealScale, OneParamPositive, OneParamPositiveScale}
TwoParam = Union{TwoParamPositive, Beta}

# The problem has already been centered and scaled.
function unconstrained_cdf(d::Distributions.Normal, x::T) where T <: Number
  (1 + erf(x/√2)) / 2
end
function unconstrained_pdf(d::Distributions.Normal, x::T) where T <: Number
  exp(-x^2/2)*(2π)^-0.5
end

function calculate_common!(x, m, d)
    if x != m.ϕ
        copy!(m.ϕ, x)
        update_αβ!(m)
        At_mul_B!(m.Vβ, m.V, m.β)
        m.δ .= unconstrained_cdf.(d, m.Vβ) .- m.w
        ToeplitzSymTriQuadForm!(m)
    end
end

function mul_tstd_x!(z, ρ, x)

    lag, current, next = 0.0, x[1], x[2]
    z[1] = current - ρ * next
    @inbounds for i ∈ 3:length(x)
        lag, current, next = current, next, x[i]
        z[i-1] = current - ρ * (lag + next)
    end
    z[end] = next - ρ * current

end

function ntscore!(stor, x, m, d)
    calculate_common!(x, m, d)

    ρ = m.Λ[2]

    mul_tstd_x!(m.nbuff1, ρ, m.δ)
    m.Vfβ .= unconstrained_pdf.(d, m.Vβ) .* m.nbuff1
    A_mul_B!(m.∇β, m.V, m.Vfβ)
    @inbounds for i ∈ 1:10
      m.∇β[i] = 2m.∇β[i] / m.Λ[1]#+ 10m.β[i]
    end
  #  m.∇β[10] -= 10
    A_mul_B!(m.∇α, m.α, m.∇β)

    for (i, grad) ∈ enumerate(m.∇α)
      stor[i] = grad
    end
    stor[1] += 2m.ϕ[1]
    stor[8] = m.Q[2] + m.n
    epsilon = eps()
    stor[9] = (m.Q[3] - imag(log_det_tstd(ρ+im*epsilon, m.n)) / epsilon) * logit_lj(m.ϕ[9])/4
    nlj_grad!(stor, m)
    stor ./= m.n
end

function ntl_likelihood!(x, m, d)
    calculate_common!(x, m, d)
    ( m.Q[1] - log_det_tstd(m.Λ[2], m.n) - 2log_jacobian(m) + abs2(m.ϕ[1])  ) / m.n + m.ϕ[8] #5dot(m.β,m.β) - 10m.β[10]
end
#+ dot(m.β, m.β) - 2m.β[10]

function ToeplitzSymTriQuadForm(δ::Vector{T}, ρ::Number, σ²::Number) where T
    δ²_sum = zero(T)
    δcross_sum = zero(T)
    δ_lag = zero(T)
    @simd for δᵢ ∈ δ
        δ²_sum += δᵢ^2
        δcross_sum += δᵢ * δ_lag
        δ_lag = δᵢ
    end
    (δ²_sum - 2ρ*δcross_sum) / 2σ²
end
function ToeplitzSymTriQuadForm!(m::MarginalBuffer)
    δ²_sum = 0.0
    δcross_sum = 0.0
    δ_lag = 0.0
    @simd for δᵢ ∈ m.δ
        δ²_sum += δᵢ^2
        δcross_sum += δᵢ * δ_lag
        δ_lag = δᵢ
    end
    σ² = m.Λ[1]
    out = (δ²_sum - 2m.Λ[2]*δcross_sum) / σ²
    m.Q[1] = out
    m.Q[2] = - out #derivative with respect to m.ϕ[8] = log(σ²)
    m.Q[3] = - 2δcross_sum / σ² #derivative with respect to ρ
    out
end

function log_det_tstd(ρ::T, n::Int) where T
    lag1 = out = one(T)
    lag2 = zero(T)
    ρ² = ρ^2
    for i ∈ 1:n
        out -= ρ² * lag2
        lag2, lag1 = lag1, out
    end
    log(out)
end

function update_βθ!(mb::MarginalBuffer)
    a = exp(mb.ϕ[1])
    c = exp(mb.ϕ[2])
    eϕ₃ = exp(mb.ϕ[3])
    b = √(3a*c)*(eϕ₃ - 1)/(eϕ₃ + 1)
    m = exp(mb.ϕ[5])
    eϕ₆ = exp(mb.ϕ[6])
    l = √(3m)*(eϕ₆ - 1)/(eϕ₆ + 1)
    mb.Λ[1] = exp(mb.ϕ[8])
    eϕ₉ = exp(mb.ϕ[9])
    mb.Λ[2] = 0.25eϕ₉/(1+eϕ₉)
    d = mb.ϕ[4]
    n = mb.ϕ[7]
    l² = l^2
    l³ = l*l²
    m² = m^2
    m³ = m*m²
    n² = n^2
    n³ = n*n²
    bmn = b*m*n; aln² = a*l*n²; amn = a*m*n; bln = b*l*n; alm = a*l*m; blm = b*l*m
    al²m = alm*l
    alm² = a*l*m²; aln = a*l*n
    mb.β[1] = a*n³ + b*n² + c*n + d
    mb.β[2] = 3amn*n + 2b*m*n + c*m
    mb.β[3] = 3aln² + 3amn*m + 2bln + b*m² + c*l
    mb.β[4] = 6alm*n + a*m³ + 3a*n² + 2blm + 2b*n + c
    mb.β[5] = 3aln*l + 3alm*m + 6amn + b*l² + 2b*m
    mb.β[6] = 3alm*l + 6aln + 3a*m² + 2b*l
    mb.β[7] = a*l³ + 6alm + 3a*n + b
    mb.β[8] = 3a*l² + 3a*m
    mb.β[9] = 3a*l
    mb.θ[1] = mb.β[10] = a
    mb.θ[2] = c
    mb.θ[3] = b
    mb.θ[4] = d
    mb.θ[5] = m
    mb.θ[6] = l
    mb.θ[7] = n
end

function update_αβ!(mb::MarginalBuffer)
    a = exp(mb.ϕ[1])
    c = exp(mb.ϕ[2])
    eϕ₃ = exp(mb.ϕ[3])
    b = √(3a*c)*(eϕ₃ - 1)/(eϕ₃ + 1)
    m = exp(mb.ϕ[5])
    eϕ₆ = exp(mb.ϕ[6])
    l = √(3m)*(eϕ₆ - 1)/(eϕ₆ + 1)
    mb.Λ[1] = exp(mb.ϕ[8])
    mb.Λ[2] = 1/( 4 * (1+exp(-mb.ϕ[9])) )
    d = mb.ϕ[4]
    n = mb.ϕ[7]
    l² = l^2
    l³ = l*l²
    m² = m^2
    m³ = m*m²
    n² = n^2
    n³ = n*n²
    bmn = b*m*n; aln² = a*l*n²; amn = a*m*n; bln = b*l*n; alm = a*l*m; blm = b*l*m
    al²m = alm*l
    alm² = a*l*m²; aln = a*l*n
    mb.β[1] = a*n³ + b*n² + c*n + d
    mb.β[2] = 3amn*n + 2b*m*n + c*m
    mb.β[3] = 3aln² + 3amn*m + 2bln + b*m² + c*l
    mb.β[4] = 6alm*n + a*m³ + 3a*n² + 2blm + 2b*n + c
    mb.β[5] = 3aln*l + 3alm*m + 6amn + b*l² + 2b*m
    mb.β[6] = 3alm*l + 6aln + 3a*m² + 2b*l
    mb.β[7] = a*l³ + 6alm + 3a*n + b
    mb.β[8] = 3a*l² + 3a*m
    mb.β[9] = 3a*l
    mb.β[10] = a
    bt = √(12a*c)*eϕ₃/(eϕ₃+1)^2
    lt = √(12m)*eϕ₆/(eϕ₆+1)^2
    mb.α[1] = a*n³ + b*n²/2
    mb.α[2] = b*n²/2 + n*c
    mb.α[3] = n²*bt
    mb.α[4] = 1
    mb.α[7] = 3n²*a + 2n*b + c

    mb.α[8] = 3amn*n + bmn
    mb.α[9] = bmn + c*m
    mb.α[10] = 2m*n*bt
    mb.α[12] = 3amn*n + 2bmn + c*m
    mb.α[14] = 6amn +2b*m

    mb.α[15] = 3aln² + 3amn*m + bln + b*m²/2
    mb.α[16] = bln + b*m²/2 + c*l
    mb.α[17] = (2l*n + m²)*bt
    mb.α[19] = 3aln²/2 + bln + c*l/2 + m*(6amn + 2b*m)
    mb.α[20] = lt*(3a*n² + 2b*n + c)
    mb.α[21] = 6a*l*n + 3a*m² + 2b*l

    mb.α[22] = 6l*amn + a*m³ + 3a*n² + blm + b*n
    mb.α[23] = blm + b*n + c
    mb.α[24] = (2l*m + 2n)*bt
    mb.α[26] = 9l*amn + 3a*m³ + 3blm
    mb.α[27] = lt*(6a*m*n + 2b*m)
    mb.α[28] = 6a*l*m + 6a*n + 2b

    mb.α[29] = 3aln*l + 3alm² + 6amn + b*l²/2 + b*m
    mb.α[30] = b*l²/2 + b*m
    mb.α[31] = (l² + 2m)*bt
    mb.α[33] = 3aln*l + 3alm²/2 + b*l² + 6alm² + 6amn + 2b*m
    mb.α[34] = lt*(6aln + 3a*m² + 2b*l)
    mb.α[35] = 3a*l² + 6a*m
    mb.α[36] = 3alm*l + 6aln + 3a*m² + b*l
    mb.α[37] = b*l
    mb.α[38] = 2l*bt
    mb.α[40] = 6alm*l + 3aln + b*l + 6a*m²
    mb.α[41] = (6alm + 6a*n + 2b)*lt
    mb.α[42] = 6a*l
    mb.α[43] = a*l³ + 6alm + 3a*n + b/2
    mb.α[44] = b/2
    mb.α[45] = bt
    mb.α[47] = l/2*(3a*l² + 6a*m) + m*6a*l
    mb.α[48] = lt*(3a*l² + 6a*m)
    mb.α[49] = 3a
    mb.α[50] = a*(3l² + 3m)
    mb.α[54] = 3a*l² + 3a*m
    mb.α[55] = 6a*l*lt
    mb.α[57] = 3a*l
    mb.α[61] = 3a*l/2
    mb.α[62] = lt*3a
    mb.α[64] = a
end

function log_jacobian(m::MarginalBuffer)
  3(m.ϕ[1]+m.ϕ[2]+m.ϕ[5])/2 - nlogit_lj(m.ϕ[3]) - nlogit_lj(m.ϕ[6]) + m.ϕ[8]-nlogit_lj(m.ϕ[9])
end
function nlj_grad!(score::AbstractVector, m::MarginalBuffer)
  score[1] -= 3.0
  score[2] -= 3.0
  score[3] -= 2dnlogit_lj(m.ϕ[3])
  score[5] -= 3.0
  score[6] -= 2dnlogit_lj(m.ϕ[6])
  score[8] -= 2.0
  score[9] -= 2dnlogit_lj(m.ϕ[9])
  score
end
#See for cubic version.
function log_jacobian_C(m::MarginalBuffer)
  3(m.ϕ[1]+m.ϕ[2])/2 - nlogit_lj(m.ϕ[3]) + m.ϕ[5]-nlogit_lj(m.ϕ[6])
end
function nlj_grad_C!(score::AbstractVector, m::MarginalBuffer)
  score[1] -= 3.0
  score[2] -= 3.0
  score[3] -= 2dnlogit_lj(m.ϕ[3])
  score[5] -= 2.0
  score[6] -= 2dnlogit_lj(m.ϕ[9])
end
function nlogit_lj(l, u, x)
  eˣ = exp(x)
  log(2 + eˣ + 1/eˣ) - log(u - l)
end
function nlogit_lj(b, x)
  eˣ = exp(x)
  log(2 + eˣ + 1/eˣ) - log(b)
end
function nlogit_lj(x)
  eˣ = exp(x)
  log(2 + eˣ + 1/eˣ)
end
function logit_lj(x)
  eˣ = exp(x)
  1/(2 + eˣ + 1/eˣ)
end
function dnlogit_lj(x)
  eˣ = exp(x)
  e²ˣ = eˣ^2
  (1 - e²ˣ) / (2eˣ + e²ˣ + 1)
end

function polyexpreval(β::AbstractVector{S}, x::T) where T where S
  x_i = convert(promote_type(Float64, T), x)
  out = β[1] + x_i*β[2]
  @inbounds for i ∈ 3:length(β)
    x_i *= x
    out += x_i * β[i]
  end
  out
end

function ∂polyexpreval(β::AbstractVector{S}, x::T) where T where S
  fₓ = convert(promote_type(S, T), β[1])
  ∂fₓ = zero(promote_type(S, T))
  x_i = one(promote_type(S, T))
  @inbounds for i ∈ 1:length(β)-1
    βᵢ = β[i+1]
    ∂fₓ += x_i * i * βᵢ
    x_i *= x
    fₓ += x_i * βᵢ
  end
  fₓ, ∂fₓ
end

function Distributions.cdf(itp::NestedPolyGLM{Distributions.Normal{Float64}}, x::Real)
  unconstrained_cdf( itp.d, polyexpreval( itp.β, (x-itp.d.μ)/itp.d.σ ) )
end
function Base.quantile(itp::NestedPolyGLM{Distributions.Normal{Float64}}, x::Real)
  nested_root( itp.θ, √(2)*erfinv(2x-1) )*itp.d.σ + itp.d.μ
end
function Distributions.pdf(itp::NestedPolyGLM{Distributions.Normal{Float64}}, x::Real)
  fₓ, ∂fₓ = ∂polyexpreval(itp.β, (x-itp.d.μ)/itp.d.σ )
  unconstrained_pdf(itp.d, fₓ) * ∂fₓ / itp.d.σ
end

function NestedPolyGLM(m::LogDensities.MarginalBuffer, d::ContinuousUnivariateDistribution)
  f = x -> ntl_likelihood!(x, m, d)
  g! = (stor, x) -> ntscore!(stor, x, m, d)
#  copy!(m.ϕ, Optim.minimizer(Optim.optimize(f, g!, m.init, NewtonTrustRegion())))
  copy!(m.init, Optim.minimizer(Optim.optimize(f, g!, m.init, BFGS(; linesearch = BackTracking()))))
  update_βθ!(m)
  NestedPolyGLM(m.β, m.θ, d)
end

function one_cubic_root(a::Real, c::Real, b::Real, d::Real)
    Δₒ = b^2 - 3a*c
    Δ₁ = 2b^3 - 9a*b*c + 27a^2*d
#    Δ₂ = Δ₁^2 - 4Δ₀^3
    C = cbrt( (Δ₁ + sqrt(Δ₁^2 - 4Δₒ^3)) / 2 )
    -(b + C + Δₒ / C) / 3a
end
function one_cubic_root(c::Real, b::Real, d::Real)
    Δₒ = b^2 - 3c
    Δ₁ = 2b^3 - 9b*c + 27d
#    Δ₂ = Δ₁^2 - 4Δ₀^3
    C = cbrt( (Δ₁ + sqrt(Δ₁^2 - 4Δₒ^3)) / 2 )
    -(b + C + Δₒ / C) / 3
end
function nested_root(β::Vector{Float64}, δ::Float64)
    r1 = one_cubic_root(β[1], β[2], β[3], β[4] - δ)
    one_cubic_root(β[5], β[6], β[7] - r1)
end
#Leaving this, because a brute force version could be useful.
#Although, I don't realistically see the case where I'd allow
#a cubic polynomial to have more than 1 root.
function find_root(β::Vector{Float64}, δ::Float64)
  a = β[4]
  b = β[3]
  c = β[2]
  d = β[1] - δ
  Δ0 = b^2 - 3a*c
  Δ1 = 2b^3 - 9a*b*c + 27a^2*d
  Δreal = Δ1^2 - 4Δ0^3
  if Δreal > 0
    Δ1prΔreal = (Δ1 + √Δreal)/2
    if Δ1prΔreal > 0
      Creal = cbrt(Δ1prΔreal)
      return - (b + Creal + Δ0/Creal)/3a
    else
      C = complex(Δ1prΔreal)^(1/3)
    end
  else
    rΔ = sqrt(complex(Δreal))
    C = complex((Δ1 + rΔ)/2)^(1/3)
  end

  ξ = -0.5 + √3/2im
  x = Vector{Complex{Float64}}(3)
  x[1] = - (b + C + Δ0/C)/3a
  for i ∈ 2:3
    C *= ξ
    x[i] = - (b + C + Δ0/C)/3a
  end
  real_roots = imag(x) .≈ 0
  if sum(real_roots) == 1
    return real(x[real_roots][1])
  elseif sum(real_roots) > 1
    rrs = real.(x[real_roots])
    mrrs = mean(rrs)
    if std(rrs) / mrrs > 1e-10
      warn("Multiple real roots in quantile estimation; returning average.")
    end
    return mrrs
  else
    throw("No real roots found.")
  end
end

function Grid(wv::weights_values)
  itp = interpolate_weight_values(wv)
  value_nodes = [linspace(itp.knots[1][1], itp.knots[1][end], 100)...]
  weight_nodes = zeros(100)
  weight_nodes[end] = 1
  for i ∈ 2:99
    weight_nodes[i] = itp[value_nodes[i]]
  end
  Grid(weight_nodes, value_nodes)
end
function Distributions.cdf(itp::Grid, x::Real)
  if x < itp.values[1]
    return 0.0
  elseif x > itp.values[end]
    return 1.0
  else
    return grid_interp(itp.values, itp.weights, searchsortedfirst(itp.values, x), x)
  end
end
function Base.quantile(itp::Grid, x::Real)
  if x <= 0
    return - Inf
  elseif x >= 1
    return Inf
  elseif x < 0.5
    i = searchsortedfirst(itp.weights, x)
  else
    i = searchsortedlast(itp.weights, x) + 1
  end
  grid_interp(itp.weights, itp.values, i, x)
end
function grid_interp(x::Array{Float64,1}, y::Array{Float64,1}, i::Int, z::Real)
  y[i-1] + (z - x[i-1]) * (y[i] - y[i-1])/(x[i]-x[i-1])
end


text:
"""
`module JAC.AngularMomentum`
    ... a submodel of JAC that contains various methods to calculate (standard) symbols, coefficients and functions from
        the theory of angular momentum and Racah's algebra.
"""
module AngularMomentum

    using SpecialFunctions, ..Basics
    using GSL: sf_coupling_3j, sf_coupling_6j, sf_coupling_9j

    """
    `AngularMomentum.allowedDoubleKappaCouplingSequence(syma::LevelSymmetry, symb::LevelSymmetry, maxKappa::Int64)`
        ... to determine all allowed coupling sequences that fulfill
            syma + (kappa1, symx, kappa2) --> symb  ==  syma + kappa1 --> symx + kappa2 --> symb,
            and where kappa1|, |kappa2| <= maxKappa. A list::Array{Tuple{Int64,LevelSymmetry,Int64},1}) of (kappa1, symx, kappa2)
            is returned.
    """
    function  allowedDoubleKappaCouplingSequence(syma::LevelSymmetry, symb::LevelSymmetry, maxKappa::Int64)
        couplings = Tuple{Int64,LevelSymmetry,Int64}[];     kx = abs(maxKappa)
        for  kappa1 = -kx:kx
            if  kappa1 == 0          continue    end
            for  kappa2 = -kx:kx
                if  kappa2 == 0      continue    end
                symx = allowedDoubleKappaSymmetries(syma, kappa1, kappa2, symb)
                for  sx in  symx    push!( couplings, (kappa1, sx, kappa2))      end
            end
        end

        return( couplings )
    end

    """
    `AngularMomentum.allowedDoubleKappaSymmetries(syma::LevelSymmetry, kappa1::Int64, kappa2::Int64, symb::LevelSymmetry)`
        ... to determine all allowed level symmetries symx that can be coupled to the sequence
            syma + kappa1 --> {symx} + kappa2 --> symb. A list::Array{LevelSymmetry,1} of symx is returned.
    """
    function  allowedDoubleKappaSymmetries(syma::LevelSymmetry, kappa1::Int64, kappa2::Int64, symb::LevelSymmetry)
        symx1 = AngularMomentum.allowedTotalSymmetries(syma, kappa1)
        symx2 = AngularMomentum.allowedTotalSymmetries(symb, kappa2)
        symx  = intersect(symx1, symx2)
        return( symx )
    end

    """
    `AngularMomentum.allowedKappaSymmetries(syma::LevelSymmetry, symb::LevelSymmetry)`
        ... to determine all allowed single-electron symmetries/partial waves kappa (l,j) that can be coupled to the given
            level symmetries. A list::Array{Int64,1} of kappa-values is returned.
    """
    function  allowedKappaSymmetries(syma::LevelSymmetry, symb::LevelSymmetry)
        kappaList = Int64[];    JList = Basics.oplus(syma.J, symb.J);    kappaMax = syma.J.num + symb.J.num + 2
        for  kappa = -kappaMax:kappaMax
             kappa == 0  &&  continue
             if  syma.parity == Basics.plus   la = 0               else   la = 1       end
             if  symb.parity == Basics.plus   lb = 0               else   lb = 1       end
             if  kappa < 0                 l  = abs(kappa) -1   else   l  = kappa   end
             isodd( la + lb + l )          && continue
             for  j in JList
                 if  j == AngularJ64(abs(kappa) - 1//2)   push!( kappaList, kappa)  end
             end
        end
        return( kappaList )
    end

   """
    `AngularMomentum.allowedMultipoleSymmetries(syma::LevelSymmetry, multipole::EmMultipole)`
        ... to determine all allowed level symmetries for which the given multipole can give rise to a non-zero (transition)
            amplitude; a symList::Array{LevelSymmetry,1} is returned.
    """
    function  allowedMultipoleSymmetries(syma::LevelSymmetry, multipole::EmMultipole)
        symList = LevelSymmetry[];    JList = Basics.oplus(syma.J, AngularJ64(multipole.L) )
        for  J in JList
            if      parityEmMultipolePi(syma.parity, multipole, Basics.plus)     push!( symList, LevelSymmetry(J, Basics.plus) )
            elseif  parityEmMultipolePi(syma.parity, multipole, Basics.minus)    push!( symList, LevelSymmetry(J, Basics.minus) )
            else    error("stop a")
            end
        end
        return( symList )
    end

    """
    `AngularMomentum.allowedTotalSymmetries(syma::LevelSymmetry, kappa::Int64)`
        ... to determine all allowed total symmetries J^P that can be constructed by coupling a partial wave kappa (l,j)
            to the given level symmetry syma. A list::Array{LevelSymmetry,1} of total symmetries is returned.
    """
    function allowedTotalSymmetries(syma::LevelSymmetry, kappa::Int64)
        symtList = LevelSymmetry[]
        if  kappa < 0    l  = abs(kappa) -1   else   l  = kappa   end
        j = AngularJ64( abs(kappa) - 1//2 )
        #
        JList = Basics.oplus(syma.J, j)
        for J in JList
            if    iseven(l)   push!( symtList, LevelSymmetry(J, syma.parity) )
            else              push!( symtList, LevelSymmetry(J, Basics.invertParity(syma.parity)) )
            end
        end
        return( symtList )
    end

    """
    `AngularMomentum.allowedTotalSymmetries(symf::LevelSymmetry, mp2::EmMultipole, mp1::EmMultipole, symi::LevelSymmetry)`
        ... to determine all allowed total symmetries J^P that can be constructed by coupling a multipole wave mp1 to the
            initial symmetry symi, and which can be further coupled with mp2 to the final symmetry symf.
            A list::Array{LevelSymmetry,1} of total symmetries is returned.
    """
    function allowedTotalSymmetries(symf::LevelSymmetry, mp2::EmMultipole, mp1::EmMultipole, symi::LevelSymmetry)
        symtList = LevelSymmetry[]
        waList   = oplus(symi.J, mp1.L);    wbList   = oplus(symf.J, mp2.L)
        for  wa in waList
            if  wa in wbList
                if      AngularMomentum.parityEmMultipolePi(symi.parity, mp1, Basics.plus)    &&
                        AngularMomentum.parityEmMultipolePi(symf.parity, mp2, Basics.plus)    push!( symtList, LevelSymmetry(wa, Basics.plus) )
                elseif  AngularMomentum.parityEmMultipolePi(symi.parity, mp1, Basics.minus)   &&
                        AngularMomentum.parityEmMultipolePi(symf.parity, mp2, Basics.minus)   push!( symtList, LevelSymmetry(wa, Basics.minus) )
                end
            end
        end

        return( symtList )
    end

    """
    `AngularMomentum.bracket(jList::Array{AngularJ64,1})`
        ... to compute the bracket [a, b, c, ... ] = (2a+1) * (2b+1) * (2b+1) * ... of the given angular momenta.
            A value::Int64 is returned.
    """
    function bracket(jList::Array{AngularJ64,1})
        value = 1;    for  j in jList    value = value * (Basics.twice(j) + 1)    end
        return( value )
    end

    """
    `AngularMomentum.ChengI

    + (kapa::Int64, ma::AngularM64, kapb::Int64, mb::AngularM64, L::AngularJ64, M::AngularM64)`
        ... evaluates the angular I (kappa m, kappa' m', LM) integral as defined by Cheng, NATO summerschool (198x),
            Eq. (A4.5), and including the full magnetic (orientational) dependence. A value::Float64 is returned.
    """
    function ChengI(kapa::Int64, ma::AngularM64, kapb::Int64, mb::AngularM64, L::AngularJ64, M::AngularM64)
        ja = Basics.subshell_j( Subshell(9, kapa) );    jb = Basics.subshell_j( Subshell(9, kapb) )   # Use principal QN n=9 arbitrarely here
        la = Basics.subshell_l( Subshell(9, kapa) );    lb = Basics.subshell_l( Subshell(9, kapb) )
        #
        # Test for parity
        if  L.den != 1   error("stop a")                end
        if  isodd( la + lb + L.num )    return( 0. )    end
        #
        wa = AngularMomentum.phaseFactor([jb, +1, L, -1, ja]) * sqrt( (Basics.twice(jb)+1)*(Basics.twice(JL)+1) / (4*pi*(Basics.twiceJ(ja)+1) ) ) *
             AngularMomentum.ClebschGordan(jb, AngularM64(1//2), L, AngularM64(0), ja, AngularM64(1//2)) *
             AngularMomentum.ClebschGordan(jb, mb, L, M, ja, ma)
        return( wa )
    end

    """
    + (kapa::Int64, kapb::Int64, L::AngularJ64)`
        ... evaluates the same angular I (kappa m, kappa' m', LM) integral but without the magnetic (orientational) dependence.
            A value::Float64 is returned.
    """
    function ChengI(kapa::Int64, kapb::Int64, L::AngularJ64)
        ja = Basics.subshell_j( Subshell(9, kapa) );    jb = Basics.subshell_j( Subshell(9, kapb) )   # Use principal QN n=9 arbitrarely here
        la = Basics.subshell_l( Subshell(9, kapa) );    lb = Basics.subshell_l( Subshell(9, kapb) )
        #
        # Test for parity
        if  L.den != 1   error("stop a")                end
        if  isodd( la + lb + L.num )    return( 0. )    end
        #
        # There occurs here two changes with regard to Cheng formulas:
        # i)  We now divide by sqrt(twoJ(jb)+1) instead of sqrt(twoJ(ja)+1) ... which is likely related to change emission - absorption
        # ii) The phase (-1)^(ja + L - jb) is replaced by (-1)^L  to get a proper phase between 1/2 --> 3/2 ME (compared to 3/2 --> 3/2) ...
        #     but which already comes from the Clebsch-Gordan
        ## wa = AngularMomentum.phaseFactor([ja, +1, L, -1, jb]) *
        wa = (-1)^L.num  * sqrt( (Basics.twice(jb)+1)*(Basics.twice(L)+1) / (4*pi) ) / sqrt(Basics.twice(jb)+1) *
             AngularMomentum.ClebschGordan(jb, AngularM64(1//2), L, AngularM64(0), ja, AngularM64(1//2))
        return( wa )
    end

    """
    `AngularMomentum.ClebschGordan(ja, ma, jb, mb, Jab, Mab)`
        ... calculates the Clebsch-Gordan coefficient  <ja, ma, jb, mb; Jab, Mab> for given quantum numbers by
            a proper call to a Wigner 3-j symbol. A value::Float64 is returned.
    """
    function ClebschGordan(ja, ma, jb, mb, Jab, Mab)
        mab = - Basics.twice(Mab) / 2
        pp  = (Basics.twice(ja) - Basics.twice(jb) + Basics.twice(Jab))/2
        cg  = (-1)^pp * sqrt(Basics.twice(Jab) + 1) * sf_coupling_3j(Basics.twice(ja), Basics.twice(jb), Basics.twice(Jab),
                                                                     Basics.twice(ma), Basics.twice(mb), Basics.twice(mab))
        return( cg )
    end

    """
    `AngularMomentum.ClebschGordan_old(ja::AngularJ64, ma::AngularM64, jb::AngularJ64, mb::AngularM64, Jab::AngularJ64, Mab::AngularM64)`
        ... calculates the Clebsch-Gordan coefficient  <ja, ma, jb, mb; Jab, Mab> for given quantum numbers by
            a proper call to a Wigner 3-j symbol; a value::Float64 is returned.
    """
    function ClebschGordan_old(ja::AngularJ64, ma::AngularM64, jb::AngularJ64, mb::AngularM64, Jab::AngularJ64, Mab::AngularM64)
        cg = 0.   ## ja-jb+Mab must be integer
        mab = AngularM64(-Mab.num, Mab.den)
        cg = AngularMomentum.phaseFactor([ja, -1, jb, +1, Mab]) * sqrt( (Basics.twice(Jab)+1) ) *
             AngularMomentum.Wigner_3j(ja, jb, Jab, ma, mb, mab)
    end

    """
    `AngularMomentum.CL_reduced_me(suba::Subshell, L::Int64, subb::Subshell)`
        ... calculates the reduced matrix element of the C^L spherical tensor <suba || C^(L) || subb>; a value::Float64
            is returned.
    """
    function  CL_reduced_me(suba::Subshell, L::Int64, subb::Subshell)
        la = Basics.subshell_l(suba);    ja2 = Basics.subshell_2j(suba);
        lb = Basics.subshell_l(subb);    jb2 = Basics.subshell_2j(subb)
        rem(ja2+1, 2) != 0    &&    error("stop a")

        redme = ((-1)^((ja2+1)/2)) * sqrt( (ja2+1)*(jb2+1) ) *
                Wigner_3j(AngularJ64(ja2//2), AngularJ64(L), AngularJ64(jb2//2),  AngularM64(1//2), AngularM64(0), AngularM64(-1//2) )

        return( redme )
    end

    """
    `AngularMomentum.CL_reduced_me_rb(suba::Subshell, L::Int64, subb::Subshell)`
        ... calculates the reduced matrix element of the C^L spherical tensor <suba || C^(L) || subb>; a value::Float64 is returned.
    """
    function  CL_reduced_me_rb(suba::Subshell, L::Int64, subb::Subshell)
        la = Basics.subshell_l(suba);    ja2 = Basics.subshell_2j(suba);
        lb = Basics.subshell_l(subb);    jb2 = Basics.subshell_2j(subb)
        if rem(la + lb + L, 2) != 0   return 0.     end

        redme = ((-1)^((jb2-2L-1)/2)) * sqrt( (jb2+1) ) *
                Wigner_3j(AngularJ64(ja2//2), AngularJ64(jb2//2), AngularJ64(L), AngularM64(1//2), AngularM64(-1//2), AngularM64(0) )

        return( redme )
    end

    """
    `AngularMomentum.isAllowedMultipole(syma::LevelSymmetry, multipole::EmMultipole, symb::LevelSymmetry)`
        ... evaluates to true if the given multipole may connect the two level symmetries, and false otherwise.
    """
    function  isAllowedMultipole(syma::LevelSymmetry, multipole::EmMultipole, symb::LevelSymmetry)
        if  isTriangle(syma.J, AngularJ64(multipole.L), symb.J)  &&
            parityEmMultipolePi(syma.parity, multipole, symb.parity)    return( true )
        else                                                            return( false )
        end
    end

    """
    `AngularMomentum.isTriangle(ja::AngularJ64, jb::AngularJ64, jc::AngularJ64)`
        ... evaluates to true if Delta(ja,jb,jc) = 1, ie. if the angular momenta ja, jb and jc can couple to each other,
            and false otherwise.
    """
    function isTriangle(ja::AngularJ64, jb::AngularJ64, jc::AngularJ64)
        if  ja.den == 1   ja2 = 2ja.num   else   ja2 = ja.num   end
        if  jb.den == 1   jb2 = 2jb.num   else   jb2 = jb.num   end
        if  jc.den == 1   jc2 = 2jc.num   else   jc2 = jc.num   end
        isodd(ja2 + jb2 + jc2)   &&    error("Angular momenta do no fullfill proper coupling rules; 2ja = $ja2, 2jb = $jb2, 2jc = $jc2")
        if  ja2 + jb2 >= jc2     &&    jc2 + ja2 >= jb2   &&   jb2 + jc2 >= ja2    return( true )
        else                                                                       return( false )
        end
    end

    """
    `AngularMomentum.isTriangle(ja::Int64, jb::Int64, jc::Int64)`
        ... evaluates to true if Delta(ja,jb,jc) = 1, ie. if the three integer (length) ja, jb and jc can form a triangle,
            and false otherwise.
    """
    function isTriangle(ja::Int64, jb::Int64, jc::Int64)
        if  ja + jb >= jc     &&    jc + ja >= jb   &&   jb + jc >= ja    return( true )
        else                                                              return( false )
        end
    end

    """
    `AngularMomentum.kappa_j(kappa::Int64)`  ... calculates the j::AngularJ64 value of a given kappa.
    """
    function  kappa_j(kappa::Int64)
        j = AngularJ64( abs(kappa) - 1//2 )
        return( j )
    end

    """
    `AngularMomentum.kappa_l(kappa::Int64)`  ... calculates the l::AngularJ64 value of a given kappa.
    """
    function  kappa_l(kappa::Int64)
        if  kappa < 0    l  = abs(kappa) -1   else   l  = kappa   end
        return( AngularJ64(l) )
    end

    """
    `AngularMomentum.m_values(j::AngularJ64)`  ... returns a list of m-values for given j::AngularJ64.
    """
    function  m_values(j::AngularJ64)
        mList = AngularM64[]
        if       j.den == 1   for  jx = -j.num:j.num      push!( mList, AngularM64(jx))       end
        elseif   j.den == 2   for  jx = -j.num:2:j.num    push!( mList, AngularM64(jx//2))    end
        else     error("stop a")
        end

        return( mList )
    end

    """
    `AngularMomentum.oneJ(ja::AngularJ64)`  ... calculates ja; a (positive) value::Float64 is returned.
    """
    function  oneJ(ja::AngularJ64)
        if  ja.den  == 1    ja1 = 1.0 * ja.num   else   ja1 = ja.num / 2.   end
        return( ja1 )
    end

    """
    `AngularMomentum.parityEmMultipolePi(pa::Parity, multipole::EmMultipole, pb::Parity)`
        ... evaluates to true if the given multipole fullfills the parity selection rule pi(a, multipole, b) = 1,
            and false otherwise. This includes a proper test for both, electric and magnetic multipoles, based on
            multipole.electric.
    """
    function  parityEmMultipolePi(pa::Parity, multipole::EmMultipole, pb::Parity)
        # The the multipolarity into account to 'interprete' electric multipoles
        if       isodd(multipole.L)    Le = multipole.electric    else   Le = !(multipole.electric)   end

        if        pa == Basics.plus    &&   Le   &&   pb == Basics.minus        return( true )
        elseif    pa == Basics.minus   &&   Le   &&   pb == Basics.plus         return( true )
        elseif    pa == pb          &&   !(Le)                                  return( true )
        else                                                                    return( false )
        end
    end

    """
    `AngularMomentum.phaseFactor(list::Array{Any,1})`
        ... checks and calculates the phase factor (-1)^(ja + mb -jc ...) that occur frequently in angular momentum theory;
            a value +1. or -1. is returned. Use phaseFactor([ja::Union{AngularJ64,AngularM64), -1, mb::Union{AngularJ64,AngularM64),
            ..., 1, jc::Union{AngularJ64,AngularM64)]) to specify the phase.
    """
    function phaseFactor(list::Array{Any,1})
        if   iseven( length(list) )                                                   error("Wrong number of arguments.")   end
        for  i = 1:2:length(list)
             if  !(typeof(list[i]) == AngularJ64  || typeof(list[i]) == AngularM64)   error("Wrong type of argument $i")    end
             if  i == length(list)                                                    break                                 end
             if  list[i+1] != 1   &&  list[i+1] != -1                error("Wrong type of argument $(i+1); must be +-1")    end
        end
        #
        ja = list[1];        if  ja.den  == 1    jm2 = 2ja.num                  else   jm2 = ja.num                         end
        for  i = 2:2:length(list)
           ja = list[i+1];   if  ja.den  == 1    jm2 = jm2 + list[i]*2ja.num    else   jm2 = jm2 + list[i]*ja.num           end
        end
        #
        if rem(jm2,2) != 0    error("Improper combination of angular momenta")  end

        return( (-1.)^(jm2/2) )
    end

    """
    `AngularMomentum.sigma_reduced_me(suba::Subshell, subb::Subshell)`
        ... calculates the reduced matrix element of the sigma^(1) spherical tensor <suba || sigma^(1) || subb>;
            a value::Float64 is returned.
    """
    function  sigma_reduced_me(suba::Subshell, subb::Subshell)
        la = Basics.subshell_l(suba);    ja2 = Basics.subshell_2j(suba);
        lb = Basics.subshell_l(subb);    jb2 = Basics.subshell_2j(subb)
        if rem(la + lb, 2) != 0   return 0.     end

        ## redme = ((-1)^((jb2-2L-1)/2)) * sqrt( (jb2+1) ) *
        ##         Wigner_3j(AngularJ64(ja2//2), AngularJ64(jb2//2), AngularJ64(L), AngularM64(1//2), AngularM64(-1//2), AngularM64(0) )
        redme = 0.

        return( redme )
    end

    """
    `AngularMomentum.sphericalYlm(l::Int64, m::Int64, theta::Float64, phi::Float64)`
        ... calculates the spherical harmonics for low l-values explicitly. A value::Complex{Float64} is returned.
            Ylm = sqrt( (2*l+1) / (two*two*pi) ) * spherical_Clm(l,m,theta,phi).
    """
    function sphericalYlm(l::Int64, m::Int64, theta::Float64, phi::Float64)
        one = 1.0 + 0.0im;    iphi = 0. + phi*im
        #
        if      l < abs(m)              ylm = 0 * one
        elseif  l ==  0  &&   m == 0    ylm = one / (2*sqrt(pi))
        elseif  l ==  1  &&   m ==  1   ylm = - sqrt(3 /(2*pi)) / 2 * sin(theta) * exp(iphi)
        elseif  l ==  1  &&   m ==  0   ylm = sqrt(3/pi) / 2 * cos(theta) * one
        elseif  l ==  1  &&   m == -1   ylm = sqrt(3/(2*pi)) / 2 * sin(theta) * exp(-iphi)
        elseif  l ==  2  &&   m ==  2   ylm = sqrt(3*5/(2*pi)) / (2*2) * sin(theta) * sin(theta) * exp(2*iphi)
        elseif  l ==  2  &&   m ==  1   ylm =  - sqrt(3*5/(2*pi)) / 2 * cos(theta) * sin(theta) * exp(iphi)
        elseif  l ==  2  &&   m ==  0   ylm = sqrt(5/pi) / 4 * (3 * cos(theta) * cos(theta) - 1) * one
        elseif  l ==  2  &&   m == -1   ylm = sqrt(3*5/(2*pi)) / 2 * cos(theta) * sin(theta) * exp(-iphi)
        elseif  l ==  2  &&   m == -2   ylm = sqrt(3*5/(2*pi)) / 4 * sin(theta) * sin(theta) * exp(-2*iphi)
        elseif  l ==  3  &&   m ==  3   ylm = - sqrt(5*7/pi) / 8 * sin(theta) * sin(theta) * sin(theta) * exp(3*iphi)
        elseif  l ==  3  &&   m ==  2   ylm = sqrt(3*5*7/(2*pi)) / 4 * cos(theta) * sin(theta) * sin(theta) * exp(2*iphi)
        elseif  l ==  3  &&   m ==  1   ylm = - sqrt(3*7/pi) / 8 * (5*cos(theta) * cos(theta) - 1) * sin(theta) * exp(iphi)
        elseif  l ==  3  &&   m ==  0   ylm = sqrt(7/pi) / 4 * (5*cos(theta) * cos(theta) - 3) * cos(theta) * one
        elseif  l ==  3  &&   m == -1   ylm = sqrt(3*7/pi) / 8 * (5*cos(theta) * cos(theta) - 1) * sin(theta) * exp(-iphi)
        elseif  l ==  3  &&   m == -2   ylm = sqrt(3*5*7/(2*pi)) / 4 * cos(theta) * sin(theta) * sin(theta) * exp(-2*iphi)
        elseif  l ==  3  &&   m == -3   ylm = sqrt(5*7/pi) / 8 * sin(theta) * sin(theta) * sin(theta) * exp(-3*iphi)
        elseif  l ==  4  &&   m ==  4   ylm =  3 * sqrt(5*7/(2*pi)) / 16 * (sin(theta)^4) * exp(4*iphi)
        elseif  l ==  4  &&   m ==  3   ylm = -3 * sqrt(5*7/pi) / 8 * (sin(theta)^3) * cos(theta) * exp(3*iphi)
        elseif  l ==  4  &&   m ==  2   ylm =  3 * sqrt(5/(2*pi)) / 8 * sin(theta) * sin(theta) * (7*cos(theta)*cos(theta) - 1) * exp(2*iphi)
        elseif  l ==  4  &&   m ==  1   ylm = -3 * sqrt(5/pi) / 8 * sin(theta) * (7*(cos(theta)^3) - 3*cos(theta)) * exp(iphi)
        elseif  l ==  4  &&   m ==  0   ylm =  3 * sqrt(1/pi) / 16 * (35*(cos(theta)^4) - 30*cos(theta)*cos(theta) + 3) * one
        elseif  l ==  4  &&   m == -1   ylm =  3 * sqrt(5/pi) / 8 * sin(theta) * (7*(cos(theta)^3) - 3*cos(theta))* exp(-iphi)
        elseif  l ==  4  &&   m == -2   ylm =  3 * sqrt(5/(2*pi)) / 8 * sin(theta) * sin(theta) * (7*cos(theta)*cos(theta) - 1) * exp(-2*iphi)
        elseif  l ==  4  &&   m == -3   ylm =  3 * sqrt(5*7/pi) / 8 * (sin(theta)^3) * cos(theta) * exp(-3*iphi)
        elseif  l ==  4  &&   m == -4   ylm =  3 * sqrt(5*7/(2*pi)) / 16 * (sin(theta)^4) * exp(-4*iphi)
        else    error("stop a")
            # Further values can be find in /home/fritzsch/fri/rabs-jena/rabs_dirac_orbital.f90
        end

        return( ylm )
    end

    """
    `AngularMomentum.triangularDelta(ia2::Int64, ib2::Int64, ic2::Int64)`
        ... calculates the tringular Delta(ja,jb,jc). The arguments in this integer function are i2a = 2*ja+1, ...
            The result is 0 if the triangular condition failes and 1 otherwise.
    """
    function triangularDelta(ia2::Int64, ib2::Int64, ic2::Int64)
        i = ib2 - ic2
        if  ia2 >= abs(i) + 1   &&   ia2 <= ib2 + ic2 - 1    return( 1 )   else    return( 0 )    end
    end

    """
    `AngularMomentum.triangularDelta(ja::AngularJ64, jb::AngularJ64, jc::AngularJ64)`
        ... calculates the tringular Delta(ja,jb,jc). The result is 0 if the triangular condition failes and 1 otherwise.
    """
    function triangularDelta(ja::AngularJ64, jb::AngularJ64, jc::AngularJ64)
        if  abs(ja.num//ja.den - jb.num//jb.den) <= jc.num//jc.den <= ja.num//ja.den + jb.num//jb.den  return( 1 )   else    return( 0 )    end
    end

    #==
    """
    `AngularMomentum.twoJ(ja::AngularJ64)`  ... calculates 2*ja; an (positive) value::Int64 is returned.
    """
    function  twoJ(ja::AngularJ64)
        if  ja.den  == 1    ja2 = 2ja.num   else   ja2 = ja.num   end
        return( ja2 )
    end

    """
    `AngularMomentum.twoM(ma::AngularM64)`  ... calculates 2*ma; a value::Int64 is returned.
    """
    function  twoM(ma::AngularM64)
        if  ma.den  == 1    ma2 = 2ma.num   else   ma2 = ma.num   end
        return( ma2 )
    end  ==#

    """
    `AngularMomentum.Wigner_DFunction(j, p, q, alpha::Float64, beta::Float64, gamma::Float64)`
        ... calculates the value of a Wigner D^j_pq (alpha, beta, gamma) for given quantum numbers and (Euler) angles
            (alpha, beta, gamma). It makes use of the small Wigner d(beta) matrix as the key part.
            A value::Float64 is returned.
    """
    function Wigner_DFunction(j, p, q, alpha::Float64, beta::Float64, gamma::Float64)
        wa = exp(-im*p*alpha -im*q*gamma) * AngularMomentum.Wigner_dmatrix(j, p, q, beta)
    end

    """
    `AngularMomentum.Wigner_dmatrix(j, mp, m, beta::Float64)`
        ... calculates the value of the small Wigner d^j_m',m (beta) for given quantum numbers and the angle beta.
            Wigner's formula is applied to evaluate the small Wigner matrix; a value::Float64 is returned.
    """
    function Wigner_dmatrix(jj, mmp, mm, beta::Float64)
        j = Float64(jj);      mp = Float64(mmp);    m = Float64(mm)
        if     mod(j+j + mp+mp + m+m, 2.0) != 0.       error("Inappropriate quantum numbers j=$j, mp=$mp, m=$m")
        elseif abs(mp) > j    ||   abs(m) > j          return( 0. )
        end
        #
        factor = sqrt( factorial(j+mp) * factorial(j-mp) * factorial(j+m) * factorial(j-m) )
        wa = 0.
        for  s = 0:100
            if  j+m-s < 0  ||  j-mp-s < 0   break       end
            if  mp - m + s < 0              continue    end
            wa = wa + (-1)^(mp-m+s) * cos(beta/2.)^(2j+m-mp-2s) * sin(beta/2.)^(mp-m+2s) /
                      ( factorial(j+m-s) * factorial(s) * factorial(mp-m+s) * factorial(j-mp-s) )
        end

        return( factor*wa )
    end

    """
    `AngularMomentum.Wigner_3j(a, b, c, m_a, m_b, m_c)`
        ... calculates the value of a Wigner 3-j symbol for given quantum numbers as displayed in many texts on the theory of
            angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981, p. 142);
            it calls the corresponding function from the GNU Scientific Library. A value::Float64 is returned.
    """
    function Wigner_3j(a, b, c, m_a, m_b, m_c)
        sf_coupling_3j(Basics.twice(a),   Basics.twice(b),   Basics.twice(c),
                       Basics.twice(m_a), Basics.twice(m_b), Basics.twice(m_c))
    end

    """
    `AngularMomentum.Wigner_3j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, m_a::AngularM64, m_b::AngularM64, m_c::AngularM64)`
         ... calculates the value of a Wigner 3-j symbol  for given quantum numbers by its algebraic formulae as displayed in many texts on the
             theory of angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981,
             p. 142); a value::Float64 is returned.
    """
    function Wigner_3j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, m_a::AngularM64, m_b::AngularM64, m_c::AngularM64)
        # Use twice the given arguments and integer arithmetik throughout.
        if  a.den == 1      ja = a.num + a.num    else    ja = a.num      end
        if  b.den == 1      jb = b.num + b.num    else    jb = b.num      end
        if  c.den == 1      jc = c.num + c.num    else    jc = c.num      end
        if  m_a.den == 1    ma = 2m_a.num         else    ma = m_a.num    end
        if  m_b.den == 1    mb = 2m_b.num         else    mb = m_b.num    end
        if  m_c.den == 1    mc = 2m_c.num         else    mc = m_c.num    end

        # Test the triangular condition and that for magnetic quantum numbers
        if      ma+mb+mc != 0  ||   triangularDelta(ja+1,jb+1,jc+1) == 0      return( 0. )
        elseif  abs(ma) > ja   ||   abs(mb) > jb   ||  abs(mc) > jc           return( 0. )
        elseif  rem(ma+ja+ja,2) != rem(ja,2)   ||   rem(mb+jb+jb,2) != rem(jb,2)   ||   rem(mc+jc+jc,2) != rem(jc,2)    error("stop a")
        end

        ik = [ ja + jb - jc,    ja - jb + jc,    -ja + jb + jc,    ja + jb + jc + 2,    ja - ma,         ja + ma,         jb - mb,
               jb + mb,         jc - mc,          jc + mc,         jb - jc - ma,        ja - jc + mb,    jc - jb + ma,    ja - jb - mc ]
        for  i  in  1:14
            rem(ik[i],2) == 1   &&   error("stop b")
            ik[i] = ik[i] / 2
        end
        # Calculate the 3-j delta factor
        delta = (factorial( ik[1] ) * factorial( ik[2] ) * factorial( ik[3]) * factorial( ik[5]) * factorial( ik[6] ) * factorial( ik[7] ) *
                 factorial( ik[8] ) * factorial( ik[9] ) * factorial(ik[10] )  /  factorial( ik[4] ) )
        # Find out the intervall of summation  k  and sum up
        kmin = max(0, ik[11], ik[12])
        kmax = min(ik[1], ik[5], ik[8] )
        qsum = 0.
        for  k in kmin:kmax
            sumk = (factorial(k) * factorial( ik[1]-k )  * factorial( ik[5]-k ) * factorial( ik[8]-k ) *
                    factorial( ik[13]+k ) * factorial( k-ik[12]) )
            sumk = 1.0 / sumk
            if  rem(k,2) == 0   qsum = qsum + sumk    else   qsum = qsum - sumk    end
        end

        if  rem(ik[14],2) != 0    w3j = - sqrt(delta) * qsum    else    w3j = sqrt(delta) * qsum    end
        return( w3j )
    end

    """
    `AngularMomentum.Wigner_6j(a, b, c, d, e, f)`
        ... calculates the value of a Wigner 6-j symbol for given quantum numbers as displayed in many texts on the theory of
            angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981, p. 142);
            it calls the corresponding function from the GNU Scientific Library. A value::Float64 is returned.
    """
    function Wigner_6j(a, b, c, d, e, f)
        sf_coupling_6j(Basics.twice(a), Basics.twice(b), Basics.twice(c),
                       Basics.twice(d), Basics.twice(e), Basics.twice(f))
    end

    """
    `AngularMomentum.Wigner_6j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, e::AngularJ64, f::AngularJ64)`
         ... calculates the value of a Wigner 6-j symbol  for given quantum numbers by its algebraic formulae as displayed in many texts on the
             theory of angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981,
             p. 142); a value::Float64 is returned.
    """
    function Wigner_6j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, ee::AngularJ64, f::AngularJ64)
        # Use twice the given arguments and integer arithmetik throughout.
        if  a.den == 1       ja = a.num + a.num    else    ja = a.num      end
        if  b.den == 1       jb = b.num + b.num    else    jb = b.num      end
        if  c.den == 1       jc = c.num + c.num    else    jc = c.num      end
        if  d.den == 1       jd = d.num + d.num    else    jd = d.num      end
        if  ee.den == 1      je = ee.num + ee.num  else    je = ee.num     end
        if  f.den == 1       jf = f.num + f.num    else    jf = f.num      end
        #
        wa = AngularMomentum.triangularDelta(ja+1, jb+1, jc+1) * AngularMomentum.triangularDelta(ja+1, je+1, jf+1) *
             AngularMomentum.triangularDelta(jd+1, jb+1, jf+1) * AngularMomentum.triangularDelta(jd+1, je+1, jc+1)
        if  wa == 0    return( 0. )    end
        #
        # Evaluate the algebraic formulae
        n1 = Int64( (ja + jb + jc)//2 );    n2 = Int64( (je + jd + jc)//2 );       n3 = Int64( (ja + je + jf)//2 )
        n4 = Int64( (jb + jd + jf)//2 );    n5 = Int64( (ja + jb + je + jd)//2 );  n6 = Int64( (ja + jd + jc + jf)//2 )
        n7 = Int64( (jb + je + jc + jf)//2 )
        kmin = max(n1,n2,n3,n4) + 1;   kmax = min(n5,n6,n7) + 1
        #
        w6j = 1.0;   icount = 0
        if  kmin != kmax
            for k = kmin+1:kmax
                ki = kmax - icount
	            w6j = 1.0 - (w6j * ki*(n5 - ki + 2.) * (n6 - ki + 2.) * (n7 - ki + 2.)) /
                            ((ki - 1. - n1) * (ki - 1. - n2) * (ki - 1. - n3) * (ki - 1. - n4))
                icount = icount + 1
            end
        end
        #
        w6j = w6j * exp( (lgamma(kmin+1) - lgamma(kmin-n1) - lgamma(kmin-n2) - lgamma(kmin-n3) - lgamma(kmin-n4) - lgamma(n5+2-kmin) -
	                      lgamma(n6+2-kmin) - lgamma(n7+2-kmin))  +
	                   ( (lgamma(n1+1-ja) + lgamma(n1+1-jb) + lgamma(n1+1-jc) - lgamma(n1+2) + lgamma(n2+1-je) + lgamma(n2+1-jd) +
	                      lgamma(n2+1-jc) - lgamma(n2+2) + lgamma(n3+1-ja) + lgamma(n3+1-je) + lgamma(n3+1-jf) - lgamma(n3+2) +
                          lgamma(n4+1-jb) + lgamma(n4+1-jd) + lgamma(n4+1-jf) - lgamma(n4+2)) / 2. ) )
        #
        if   rem(n5+kmin,2) == 0       w6j = -w6j   end
        if   rem(ja+jb+je+jd,4) != 0   w6j = -w6j   end

        return( w6j )
    end

    """
    `AngularMomentum.Wigner_9j(a, b, c, d, e, f, g, h, i)`
        ... calculates the value of a Wigner 3-j symbol for given quantum numbers as displayed in many texts on the theory of
            angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981, p. 142);
            it calls the corresponding function from the GNU Scientific Library. A value::Float64 is returned.
    """
    function Wigner_9j(a, b, c, d, e, f, g, h, i)
        sf_coupling_9j(Basics.twice(a), Basics.twice(b), Basics.twice(c),
                       Basics.twice(d), Basics.twice(e), Basics.twice(f),
                       Basics.twice(g), Basics.twice(h), Basics.twice(i))
    end

    """
    `AngularMomentum.Wigner_9j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, e::AngularJ64, f::AngularJ64,
                                   g::AngularJ64, h::AngularJ64, i::AngularJ64)`
         ... calculates the value of a Wigner 6-j symbol  for given quantum numbers by its algebraic formulae as displayed in many texts on the
             theory of angular momentum (see R. D. Cowan, The Theory of Atomic Structure and Spectra; University of California Press, 1981,
             p. 142); a value::Float64 is returned.
    """
    function Wigner_9j_old(a::AngularJ64, b::AngularJ64, c::AngularJ64, d::AngularJ64, ee::AngularJ64, f::AngularJ64,
                       g::AngularJ64, h::AngularJ64, i::AngularJ64)
        Defaults.warn(AddWarning, "AngularMomentum.Wigner_9j(): Here, we wish to call an external Julia package/interface; now set to 2.0")
        w9j = 2.0
        return( w9j )
    end

end # module

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_IdealClutch

Simulate models/IdealClutch.jl model.

Demonstrates/tests varying number of constraint equations fc leading to Dirac impulses at events
"""
module Simulate_IdealClutch

import ..ModiaMath
include(joinpath("models", "IdealClutch.jl"))
import .IdealClutch

model  = IdealClutch.Model()
result = ModiaMath.simulate!(model, stopTime=500.0; log=true)
ModiaMath.plot(result, [("inertia1.w", "inertia2.w", "capacitor.v", "engaged")], heading="Simulate_IdealClutch.jl")

        w1_end = result["inertia1.w"][end]
		w2_end = result["inertia2.w"][end]

        println("w1_end = ", w1_end, ", w2_end = ", w2_end)
end


text:
#  Well trying to use Varpro for Gaussian Mixture models is problematic.  We can
#  get away with that here because we use only use 3 single variate Gaussians
#  that are well separated and have comparable magnitudes.  Even then, I had to
#  run the random intial value version about a dozen times to find an initial starting
#  point that converges to the correct value

using Varpro
using Base.Iterators

phi(x, μ, σ) = exp(-((x - μ)/σ)^2/2)
g(x, μ, σ)   = phi(x, μ, σ) / (σ * sqrt(2π))
dμ(x, μ, σ) = (x - μ)   * phi(x, μ, σ) / (σ^3 * sqrt(2π))
dσ(x, μ, σ) = (x - μ)^2 * phi(x, μ, σ) / (σ^4 * sqrt(2π))

# We will use these params to form the synthetic gaussian mixture data.
# We have a mixture of 3 single variate gaussians and these are the
# parameters that we expect to recover using Varpro.
const mu   = [2.0, 3.0, 0.0]
const sig  = [0.5, 0.2, 1.0]
const cc   = [0.25, 0.25, 0.50]

# for creating synthetic data
function mixn(x)
    sum = 0.0
    for i in 1:length(cc)
        sum +=  cc[i] * g(x, mu[i], sig[i])
    end
    return sum
end

# The phi matrix.  Each row is for a different data point.  Each column
# is the contribution from a single gaussian
function f_gaussian(alpha, ctx)
    for j = 1:ctx.n
        for i in 1:ctx.m
           ctx.phi[i, j] = g(ctx.t[i], alpha[2j-1], alpha[2j])
        end
    end
    ctx.phi
end

function g_gaussian(alpha, ctx)
    for j = 1:ctx.n
        for i in 1:ctx.m
            ctx.dphi[i, 2j-1] = dμ(ctx.t[i], alpha[2j-1], alpha[2j])
            ctx.dphi[i, 2j]   = dσ(ctx.t[i], alpha[2j-1], alpha[2j])
        end
    end
    ctx.dphi
end

function testgauss(;verbose=true)
    n = 3  # length(c)  # the number of gaussian components

    # Our synthetic data
    t = [-10.0:0.01:10.0;]
    y = mixn.(t)
    m = length(t)

    # we might want to use different weights for different problems here
    w = ones(m)

    # this _sometimes_ converges to correct result
    #x_init = 2.0 * rand(2n)
    #@show x_init

    # carefully chosen starting point (actually one of the random ones that worked)
    x_init = [1.680, 0.744, 0.030, 0.312, 0.259, 1.427]

    # tell varpro how the jacobian is structured.
    ind = [collect(flatten((i,i) for i in 1:n))'; [1:2n;]']

    ctx = FitContext(y, t, w, x_init, n, ind, f_gaussian, g_gaussian)
    ctx.skip_stats = true
    ctx.verbose = verbose
    alpha, cfit, wresid, resid_norm, y_est, reg = varpro(ctx)
    return all(isapprox.(cc, sort(cfit))) && all(isapprox.(sort([mu;sig]), sort(alpha), atol=1e-13))
end


text:
using LinearAlgebra, SparseArrays

# State ∈ [0, up, down, updown] = [00, 10, 01, 11]
# 2x2 lattice -> 4 (exponent) -> 4^4 states = 256

function state_from_integer(int, sites=4, substates_per_site=2)
    int > (2substates_per_site)^sites-1 && return
    out = BitArray(undef, (substates_per_site, sites))
    out.chunks[1] = int
    out
end

function create(state, site, substate)
    # No state case
    state == 0 && return 0, 0.0
    # create(|1⟩) -> no state
    state[substate, site] && return 0, 0.0

    # fermionic sign
    lindex = 2(site-1) + substate
    _sign = iseven(sum(state[1:lindex-1])) ? +1.0 : -1.0

    # create(|0⟩) -> |1⟩
    s = copy(state)
    s[substate, site] = true
    return _sign, s
end

function annihilate(state, site, substate)
    # No state
    state == 0 && return 0, 0.0
    # annihilate(|0⟩) -> no state
    !state[substate, site] && return 0, 0.0

    # fermionic sign
    lindex = 2(site-1) + substate
    _sign = iseven(sum(state[1:lindex-1])) ? +1.0 : -1.0

    # annihilate(|1⟩) -> |0⟩
    s = copy(state)
    s[substate, site] = false
    return _sign, s
end

function HamiltonMatrix(model::HubbardModelAttractive)
    lattice = model.l
    t = model.t
    U = model.U
    mu = model.mu

    H = zeros(Float64, 4^lattice.sites, 4^lattice.sites)

    # -t ∑_ijσ c_iσ^† c_jσ
    # +U ∑_i (n_i↑ - 1/2)(n_i↓ - 1/2)
    # -μ ∑_i n_i
    for i in 1:4^lattice.sites
        lstate = state_from_integer(i-1, lattice.sites)
        for j in 1:4^lattice.sites
            rstate = state_from_integer(j-1, lattice.sites)

            E = 0
            # hopping (hermitian conjugate implied/included by lattice generation)
            for substate in [1, 2]
                for source in 1:lattice.sites
                    for target in lattice.neighs[:, source]
                        _sign1, state = annihilate(rstate, source, substate)
                        _sign2, state = create(state, target, substate)
                        if state != 0 && lstate == state
                            E -= _sign1 * _sign2 * t
                        end
                    end
                end
            end

            # # U, μ terms
            for p in 1:4
                up_occ = rstate[1, p]
                down_occ = rstate[2, p]
                if lstate == rstate
                    E += U * (up_occ - 0.5) * (down_occ - 0.5)
                    E -= mu * (up_occ + down_occ)
                end
            end

            H[i, j] = E
        end
    end
    H
end

# Greens function

function Greens(site1, site2, substate1, substate2)
    s -> begin
        _sign1, state = create(s, site1, substate1)
        state == 0 && return typeof(s)[], Float64[]
        _sign2, state = annihilate(state, site2, substate2)
        state == 0 && return typeof(s)[], Float64[]
        return [state], [_sign1 * _sign2]
    end
end

# According to p188 QMCM this is equivalent
# can be used to check if Greens/ED is correct
function Greens_permuted(site1, site2, substate1, substate2)
    s -> begin
        _sign1, state = annihilate(s, site2, substate2)
        _sign2, state = create(state, site1, substate1)
        delta = ((site1 == site2) && (substate1 == substate2)) ? 1.0 : 0.0

        if state == 0 && delta == 0.0
            # off-diagonal
            return typeof(s)[], Float64[]
        elseif state == 0
            # only delta function triggers (|0⟩ state)
            return [s], [delta]
        else
            # both trigger (|1⟩ state)
            return [state], [delta - _sign1 * _sign2]
        end
    end
end

function expectation_value(
        observable::Function,
        H;
        T=1.0, beta = 1.0 / T,
        N_sites = 4,
        N_substates = 2
    )

    vals, vecs = eigen(H)
    Z = 0.0
    O = 0.0
    for i in eachindex(vals)
        # exp(βEᵢ)
        temp = exp(-beta * vals[i])
        Z += temp

        # ⟨ψᵢ|Ô|ψᵢ⟩
        right_coefficients = zeros(eltype(vecs), size(vecs, 1))
        for j in 1:size(vecs, 1)
            state = state_from_integer(j-1, N_sites, N_substates)
            states, values = observable(state)
            for (s, v) in zip(states, values)
                # Assuming no (s, v) pair if state destroyed
                k = s.chunks[1]+1
                right_coefficients[k] += v * vecs[j, i]
            end
        end
        O += temp * dot(vecs[:, i], right_coefficients)
    end
    O / Z
end

function calculate_Greens_matrix(H, lattice; beta=1.0, N_substates=2)
    G = Matrix{Float64}(
        undef,
        lattice.sites*N_substates,
        lattice.sites*N_substates
    )
    for substate1 in 1:N_substates, substate2 in 1:N_substates
        for site1 in 1:lattice.sites, site2 in 1:lattice.sites
            G[
                lattice.sites * (substate1-1) + site1,
                lattice.sites * (substate2-1) + site2
            ] = expectation_value(
                Greens(site1, site2, substate1, substate2),
                H,
                beta = beta,
                N_sites = lattice.sites,
                N_substates=N_substates
            )
        end
    end
    G
end

# utility

function state2string(state)
    sub, L = size(state)
    spin = ["↑", "↓"]
    chunks = String[]
    for p in 1:L
        str = ""
        for s in 1:2
            str *= state[s, p] ? spin[s] : "⋅"
        end
        push!(chunks, str)
    end
    join(chunks, " ")
end


text:
using Base.Test
using Attributes

g = 5;
threshold = 0.55
v = 3
b = UnorderedAttributes(g, threshold, v)
nodes = 3;
attr = get_attributes(b, nodes)

@test min(attr...) == 1 #test rozmiaru
@test max(attr...) == v #test wartosci powinny byc +-1

attr = ones(Int, nodes, g)
attr[1, 1:4] = 3
attr[2, 1:3] = 3
attr[2, 4:5] = 2
#attr is:
# [3  3  3  3  1
#  3  3  3  2  2
#  1  1  1  1  1]
#This should give following weights:
# [0.0  0.1  -0.7
#  0.0  0.0  -1.0
#  0.0  0.0   0.0]
exp_weights = zeros(3, 3)
exp_weights[1, 2] = 0.1;
exp_weights[1, 3] = -0.7;
exp_weights[2, 3] = -1.1;

weights = get_attribute_layer_weights(b, attr)

@test size(weights) == (nodes, nodes)
@test weights ≈ exp_weights

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
using Intervals
using TimeZones
using Compat
using Compat.Test
using Compat.Dates

# The method based on keyword arguments is not in Compat, so to avoid
# deprecation warnings on 0.7 we need this little definition.
if VERSION < v"0.7.0-DEV.4524"
    function sprint(f::Function, args...; context=nothing)
        if context !== nothing
            Base.sprint((io, args...) -> f(IOContext(io, context), args...), args...)
        else
            Base.sprint(f, args...)
        end
    end
end

# The name of the module is no longer printed for types reachable from Main
if VERSION < v"0.7.0-DEV.2657"
    mod_prefix = "Intervals."
    tz_prefix = "TimeZones."
else
    mod_prefix = tz_prefix = ""
end

@testset "Intervals" begin
    include("inclusivity.jl")
    include("endpoint.jl")
    include("interval.jl")
    include("anchoredinterval.jl")
    include("comparisons.jl")
end


text:
############################################################################################
"""
$(TYPEDEF)

Holds information about updating tempering parameter, parameter for basic logistic function.

# Fields
$(TYPEDFIELDS)
"""
struct TemperingParameter{T<:AbstractFloat}
    L::T
    k::T
    x₀::T
    function TemperingParameter(L::T, k::T, x₀::T) where {T<:AbstractFloat}
        ArgCheck.@argcheck 0 < L
        ArgCheck.@argcheck 0 < k
        return new{T}(L, k, x₀)
    end
end
function TemperingParameter(::Type{T}, growthrate::F, iterations::Integer; normalizer::Float64 = 100.0
) where {F<:AbstractFloat, T<:AbstractFloat}
    return TemperingParameter(T(1.0), T(growthrate*normalizer/iterations), T(growthrate*iterations/normalizer))
end

############################################################################################
"""
$(SIGNATURES)
Initialize temperature.

# Examples
```julia
```

"""
function init(::Type{T}, adaption::UpdateFalse, val::F, parameter::TemperingParameter
) where {T<:AbstractFloat, F<:AbstractFloat}
    return T(val)
end
function init(::Type{T}, adaption::UpdateTrue, val::F, parameter::TemperingParameter
) where {T<:AbstractFloat, F<:AbstractFloat}
    return update(parameter, 1)
end

############################################################################################
"""
$(SIGNATURES)
Calculate new temperature based on current iteration.

# Examples
```julia
```

"""
function update(param::TemperingParameter{T}, index::Integer) where {T<:AbstractFloat}
    @unpack L, k, x₀ = param
    return L / (1 + exp(-k*(T(index) - x₀)))
end

############################################################################################
"""
$(SIGNATURES)
Helper function to determine temperature shape given param variables.

# Examples
```julia
```

"""
function checktemperature(param::TemperingParameter, iterations::Integer)
    return map(iter -> update(param, iter), Base.OneTo(iterations))
end

############################################################################################
"""
$(TYPEDEF)

Tuning container for iteration uptdates of temperature.

# Fields
$(TYPEDFIELDS)
"""
struct IterationTempering{A<:UpdateBool, T<:AbstractFloat} <: TemperingMethod
    "Checks if temperature will be updated after new iterations"
    adaption::A
    "Initial temperature for tempering."
    initial::ValueHolder{T}
    "Tuning parameter for temperature adjustment."
    parameter::TemperingParameter{T}
    function IterationTempering(adaption::A, val::ValueHolder{T}, parameter::TemperingParameter{T}
    ) where {A<:UpdateBool,T<:AbstractFloat}
        @argcheck 0.0 < val.current <= 1.0
        return new{A,T}(adaption, val, parameter)
    end
end

function IterationTempering(::Type{T}, adaption::B, val::F, idx::Integer
) where {T<:AbstractFloat, B<:UpdateBool, F<:AbstractFloat}
    # Assign temperature adaption
    parameter = TemperingParameter(T, val, idx)
    # Assign initial temperature
    val₀ = init(T, adaption, val, parameter)
    # Return tuning struct
    return IterationTempering(adaption, ValueHolder(val₀), parameter)
end

############################################################################################
function update!(tempering::IterationTempering, adaption::UpdateTrue, index::Integer)
    # Compute new temperature
    return update(tempering.parameter, index)
end
function update!(tempering::IterationTempering, adaption::UpdateFalse, index::Integer)
    return tempering.initial.current
end
function update!(tempering::IterationTempering, index::Integer)
    return update!(tempering, tempering.adaption, index)
end

############################################################################################
"""
$(SIGNATURES)
Return initial temperature of tuning container.

# Examples
```julia
```

"""
function initial(tempertune::IterationTempering)
    return tempertune.initial.current
end

############################################################################################
# Export
export
    IterationTempering,
    TemperingParameter,
    init,
    initial,
    update,
    update!


text:
export Map_Quadratic
export eval_map

# ****************************************************************

@doc raw"""
```julia
Map_Quadratic(X_hat_A => X_A, X_hat_B => X_B, X_hat_C => X_C)

Map_Quadratic(X_hat_A, X_hat_B, X_hat_C)
```

Same semantic as [`Map_Affine`](@ref)
"""
struct Map_Quadratic{T} <: Abstract_Map
    _hat_to_A::Pair{T,T}
    _hat_to_B::Pair{T,T}
    _hat_to_C::Pair{T,T}
end

function Map_Quadratic(x_hat_A::T,x_hat_B::T,x_hat_C::T) where {T}
    Map_Quadratic(x_hat_A=>one(T),
                  x_hat_B=>one(T),
                  x_hat_C=>one(T))
end

parameter_size(m::Map_Quadratic) = 3

function eval_map(m::Map_Quadratic{T},
                  X_hat::AbstractVector{T_HAT},
                  θ::AbstractVector{T_θ}) where {T_HAT,T,T_θ}

    @assert length(θ) == parameter_size(m)

    X_T = promote_type(T_HAT,T,T_θ)

    θ_A = θ[1]
    θ_B = θ[2]
    θ_C = θ[3]

    X_hat_A, X_A = m._hat_to_A
    X_hat_B, X_B = m._hat_to_B
    X_hat_C, X_C = m._hat_to_C

    basis_A(x) = (X_hat_B - x)*(X_hat_C - x)/((X_hat_B - X_hat_A)*(X_hat_C - X_hat_A))
    basis_B(x) = (x - X_hat_A)*(X_hat_C - x)/((X_hat_B - X_hat_A)*(X_hat_C - X_hat_B))
    basis_C(x) = (x - X_hat_A)*(x - X_hat_B)/((X_hat_C - X_hat_A)*(X_hat_C - X_hat_B))

    transform(x) = basis_A(x)*X_A*θ_A + basis_B(x)*X_B*θ_B + basis_C(x)*X_C*θ_C

    X=similar(X_hat,X_T)
    @. X = transform(X_hat)

    X
end


text:
import NonUniformRandomVariateGeneration.samplePoisson
import Random.GLOBAL_RNG

## linear is essentially an implementation of the Huber (2016) algorithm
"""
    linear(f::F, C::Float64, ϵ::Float64, rng::RNG=GLOBAL_RNG) where {F<:Function, RNG<:AbstractRNG}

Input:
  - ```f()``` simulates Bernoulli(p) random variates (of type Bool)
  - ```C >= 0```
  - ```ϵ ∈ (0,1)``` satisfying C*p < 1-ϵ.
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```X, flips```
  - ```X``` is a Bernoulli(C*p) variate
  - ```flips``` is the number of calls of f()
"""
function linear(f::F, C::Float64, ϵ::Float64, rng::RNG=GLOBAL_RNG) where
  {F<:Function, RNG<:AbstractRNG}
  ϵ <= 0.0 || ϵ >= 1.0 && throw(DomainError(ϵ, "ϵ must be in (0,1)"))
  if C > 1.0
    return _huber2016(f, C, ϵ, rng)
  end
  if rand(rng) < C
    return f(), 1
  else
    return false, 0
  end
end

## inverse is an implementation of Algorithm 6 of
## Lee, A., Doucet, A. and Łatuszyński, K., 2014. Perfect simulation using
## atomic regeneration with application to Sequential Monte Carlo.
## arXiv:1407.5770
"""
    inverse(f::F, C::Float64, ϵ::Float64, rng::RNG=GLOBAL_RNG) where {F<:Function, RNG<:AbstractRNG}

Input:
  - ```f()``` simulates Bernoulli(p) random variates (of type Bool)
  - ```C ∈ [0,p)```
  - ```ϵ ∈ (0,1)``` satisfying p > C + ϵ.
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```X, flips```
  - ```X``` is a Bernoulli(C/p) variate
  - ```flips``` is the number of calls of f()
"""
function inverse(f::F, C::Float64, ϵ::Float64, rng::RNG=GLOBAL_RNG) where
  {F<:Function, RNG<:AbstractRNG}

  (C < 0.0 || C >= 1.0) && throw(DomainError(C, "C must be in [0,1)"))
  ϵ <= 0.0 && throw(DomainError(ϵ, "ϵ must be positive"))
  C+ϵ > 1.0 && throw(DomainError(C+ϵ, "C+ϵ cannot exceed 1.0"))

  fq = !f
  Cprime = 1.0/(1.0-C)
  ϵprime = ϵ/(1.0-C)
  flips::Int64 = 0
  k::Int64 = _sampleGeometric(C, rng)
  i::Int64 = 1
  while i < k
    X, newFlips = linear(fq, Cprime, ϵprime, rng)
    flips += newFlips
    if !X
      return false, flips
    end
    i += 1
  end
  return true, flips
end

## use Mendo's power algorithm
"""
    power(f::F, a::Float64, rng::RNG=GLOBAL_RNG) where {F<:Function, RNG<:AbstractRNG}

Input:
  - ```f()``` simulates Bernoulli(p) random variates (of type Bool)
  - ```a >= 0```
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```X, flips```
  - ```X``` is a Bernoulli(p^a) variate
  - ```flips``` is the number of calls of f()
"""
function power(f::F, a::Float64, rng::RNG=GLOBAL_RNG) where
  {F<:Function, RNG<:AbstractRNG}
  a < 0.0 && throw(DomainError(a, "a cannot be negative"))
  a == 0.0 && return true, 0
  a < 1.0 && return _mendoPower(f, a, rng)
  a == 1.0 && return f(), 1
  k = ceil(Int64, a)
  flips::Int64 = 0
  for i in 1:k
    X, newFlips = _mendoPower(f, a/k, rng)
    flips += newFlips
    !X && return false, flips
  end
  return true, flips
end

## use Mendo's power algorithm
"""
    sqrt(f::F, rng::RNG=GLOBAL_RNG) where {F<:Function, RNG<:AbstractRNG}

Input:
  - ```f()``` simulates Bernoulli(p) random variates (of type Bool)
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```X, flips```
  - ```X``` is a Bernoulli(sqrt(p)) variate
  - ```flips``` is the number of calls of f()
"""
function sqrt(f::F, rng::RNG=GLOBAL_RNG) where {F<:Function, RNG<:AbstractRNG}
  return _mendoPower(f, 0.5, rng)
end

## expMinus is based on series expansion of exp(-λp) = exp(λ(1-p))/exp(λ),
## which is a natural extension of Wastlund's observation for λ = 1
"""
    expMinus(f::F, λ::Float64, rng::RNG=GLOBAL_RNG) where {F<:Function, RNG<:AbstractRNG}

Input:
  - ```f()``` simulates Bernoulli(p) random variates (of type Bool)
  - ```λ >= 0```
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```X, flips```
  - ```X``` is a Bernoulli(exp(-λ*p)) variate
  - ```flips``` is the number of calls of f()
"""
function expMinus(f::F, λ::Float64, rng::RNG=GLOBAL_RNG) where
  {F<:Function, RNG<:AbstractRNG}
  λ < 0.0 && throw(DomainError(λ, "λ cannot be negative"))
  flips::Int64 = 0
  K::Int64 = samplePoisson(λ, rng)
  for i in 1:K
    if f()
      return false, flips + 1
    end
    flips += 1
  end
  return true, flips
end

## use the algorithm described in Huber (2017)
"""
    logistic(f::F, C::Float64, rng::RNG=GLOBAL_RNG) where {F<:Function, RNG<:AbstractRNG}

Input:
  - ```f()``` simulates Bernoulli(p) random variates (of type Bool)
  - ```C >= 0```
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```X, flips```
  - ```X``` is a Bernoulli(C*p/(1.0+C*p)) variate
  - ```flips``` is the number of calls of f()
"""
function logistic(f::F, C::Float64, rng::RNG=GLOBAL_RNG) where
  {F<:Function, RNG<:AbstractRNG}
  return _logistic(f, C, rng)
end

## described in
## Gonçalves, F.B., Łatuszyński, K.G. and Roberts, G.O., 2017. Exact Monte Carlo
## likelihood-based inference for jump-diffusion processes. arXiv:1707.00332
"""
    twocoin(f1::F1, f2::F2, C1::Float64, C2::Float64, rng::RNG=GLOBAL_RNG) where {F1<:Function, F2<:Function, RNG<:AbstractRNG}

Input:
  - ```f1()``` simulates Bernoulli(p1) random variates (of type Bool)
  - ```f2()``` simulates Bernoulli(p2) random variates (of type Bool)
  - ```C1, C2 >= 0```
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```X, flips```
  - ```X``` is a Bernoulli(C1*p1/(C1*p1+C2*p2)) variate
  - ```flips``` is the number of calls of f()
"""
function twocoin(f1::F1, f2::F2, C1::Float64, C2::Float64,
  rng::RNG=GLOBAL_RNG) where {F1<:Function, F2<:Function, RNG<:AbstractRNG}
  Cprob::Float64 = C1/(C1+C2)
  flips1::Int64 = 0
  flips2::Int64 = 0
  while true
    if rand(rng) < Cprob
      flips1 += 1
      if f1()
        return true, flips1, flips2
      end
    else
      flips2 += 1
      if f2()
        return false, flips1, flips2
      end
    end
  end
end

"""
signedEstimate(μ::F1, φ::F2, c::Float64, δ::Float64, n::Int64, rng::RNG=GLOBAL_RNG) where {F1<:Function, F2<:Function, RNG<:AbstractRNG} where {F1<:Function, F2<:Function, RNG<:AbstractRNG}

Input:
  - ```μ()``` simulates a random variate X
  - ```φ(X)``` is real-valued
  - ```c``` satisfies sup_x |φ(x)| <= c
  - ```δ > 0``` satisfies μ(φ) >= δ
  - ```n``` specifies a number of X variables to average
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```Y, flips, calls```
  - ```Y``` is almost surely valued in [0.0,c] with E[Y] = μ(φ)
  - ```flips``` is the number of coin tosses by the Bernoulli Factory algorithm
  - ```calls``` is the number of calls of μ() in total
"""
function signedEstimate(μ::F1, φ::F2, c::Float64, δ::Float64, n::Int64,
  rng::RNG=GLOBAL_RNG) where {F1<:Function, F2<:Function, RNG<:AbstractRNG}
  δ >= c && throw(DomainError("requires δ < c"))
  return _signedEstimate(μ, φ, c, δ, n, rng)
end

"""
signedEstimate(μ::F1, φ::F2, a::Float64, b::Float64, δ::Float64, c::Float64, n::Int64, rng::RNG=GLOBAL_RNG) where {F1<:Function, F2<:Function, RNG<:AbstractRNG}

Input:
  - ```μ()``` simulates a random variate X
  - ```φ(X)``` is real-valued
  - ```a,b,δ,c``` satisfy a <= inf_x φ(x) < b < δ <= μ(φ) < sup_x φ(x) <= c
  - ```n``` specifies a number of X variables to average
  - ```[rng]``` the RNG to used by the algorithm (does not affect f)
Output: ```Y, flips, calls```
  - ```Y``` is almost surely valued in [b, max{2b-a,c}] with E[Y] = μ(φ)
  - ```flips``` is the number of coin tosses by the Bernoulli Factory algorithm
  - ```calls``` is the number of calls of μ() in total
"""
function signedEstimate(μ::F1, φ::F2, a::Float64, b::Float64, δ::Float64,
  c::Float64, n::Int64, rng::RNG=GLOBAL_RNG) where {F1<:Function, F2<:Function,
  RNG<:AbstractRNG}
  !(a < b < δ < c) && throw(DomainError("requires a < b < δ < c"))
  return _signedEstimate(μ, φ, a, b, δ, c, n, rng)
end


text:
mulexprcost(::Number) = 0
mulexprcost(::Symbol) = 1
function mulexprcost(ex::Expr)
    base = ex.head === :call ? 10 : 1
    base + length(ex.args)
end
function mul_fast_expr(args)
    b = Expr(:call, :mul_fast)
    for i ∈ 2:length(args)
        push!(b.args, args[i])
    end
    b
end
function mulexpr(mulexargs)
    a = (mulexargs[1])::Union{Symbol,Expr,Number}
    if length(mulexargs) == 2
        return (a, mulexargs[2]::Union{Symbol,Expr,Number})
    elseif length(mulexargs) == 3
        # We'll calc the product between the guesstimated cheaper two args first, for better out of order execution
        b = (mulexargs[2])::Union{Symbol,Expr,Number}
        c = (mulexargs[3])::Union{Symbol,Expr,Number}
        ac = mulexprcost(a)
        bc = mulexprcost(b)
        cc = mulexprcost(c)
        maxc = max(ac, bc, cc)
        if ac == maxc
            return (a, Expr(:call, :mul_fast, b, c))
        elseif bc == maxc
            return (b, Expr(:call, :mul_fast, a, c))
        else
            return (c, Expr(:call, :mul_fast, a, b))
        end
    else
        return (a, mul_fast_expr(mulexargs))
    end
    a = (mulexargs[1])::Union{Symbol,Expr,Number}
    b = if length(mulexargs) == 2 # two arg mul
        (mulexargs[2])::Union{Symbol,Expr,Number}
    else
        mul_fast_expr(mulexargs)
    end
    a, b
end
function append_args_skip!(call, args, i, mod)
    for j ∈ eachindex(args)
        j == i && continue
        push!(call.args, args[j])
    end
    call
end

fastfunc(f) = get(VectorizationBase.FASTDICT, f, f)
function muladd_arguments!(argv, mod, f = first(argv))
    if f === :*
        argv[1] = :mul_fast
    else
        argv[1] = fastfunc(f)
    end
    for i ∈ 2:length(argv)
        a = argv[i]
        a isa Expr || continue
        argv[i] = capture_muladd(a::Expr, mod)
    end
end

function recursive_muladd_search!(call, argv, mod, cnmul::Bool = false, csub::Bool = false)
    if length(argv) < 3
        muladd_arguments!(argv, mod)
        return length(call.args) == 4, cnmul, csub
    end
    fun = first(argv)
    isadd = fun === :+ || fun === :add_fast || fun === :vadd || (fun == :(Base.FastMath.add_fast))::Bool
    issub = fun === :- || fun === :sub_fast || fun === :vsub || (fun == :(Base.FastMath.sub_fast))::Bool
    if !(isadd | issub)
        muladd_arguments!(argv, mod, fun)
        return length(call.args) == 4, cnmul, csub
    end
    exargs = @view(argv[2:end])
    issub && @assert length(exargs) == 2
    for (i,ex) ∈ enumerate(exargs)
        if ex isa Expr && ex.head === :call
            exa = ex.args
            f = first(exa)
            exav = @view(exa[2:end])
            if f === :* || f === :mul_fast || f === :vmul || (f == :(Base.FastMath.mul_fast))::Bool
                a, b = mulexpr(exav)
                call.args[2] = a
                call.args[3] = b
                if length(exargs) == 2
                    push!(call.args, exargs[3 -  i])
                else
                    push!(call.args, append_args_skip!(Expr(:call, :add_fast), exargs, i, mod))
                end
                if issub
                    csub = i == 1
                    cnmul = !csub
                end
                return true, cnmul, csub
            elseif isadd
                found, cnmul, csub = recursive_muladd_search!(call, exa, mod)
                if found
                    if csub
                        call.args[4] = if length(exargs) == 2
                            Expr(:call, :sub_fast, exargs[3 - i], call.args[4])
                        else
                            Expr(:call, :sub_fast, append_args_skip!(Expr(:call, :add_fast), exargs, i, mod), call.args[4])
                        end
                    else
                        call.args[4] = append_args_skip!(Expr(:call, :add_fast, call.args[4]), exargs, i, mod)
                    end
                    return true, cnmul, false
                end
            elseif issub
                found, cnmul, csub = recursive_muladd_search!(call, exa, mod)
                if found
                    if i == 1
                        if csub
                            call.args[4] = Expr(:call, :add_fast, call.args[4], exargs[3 - i])
                        else
                            call.args[4] = Expr(:call, :sub_fast, call.args[4], exargs[3 - i])
                        end
                    else
                        cnmul = !cnmul
                        if csub
                            call.args[4] = Expr(:call, :add_fast, exargs[3 - i], call.args[4])
                        else
                            call.args[4] = Expr(:call, :sub_fast, exargs[3 - i], call.args[4])
                        end
                        csub = false
                    end
                    return true, cnmul, csub
                end
            end
        end
    end
    length(call.args) == 4, cnmul, csub
end

function capture_a_muladd(ex::Expr, mod)
    call = Expr(:call, Symbol(""), Symbol(""), Symbol(""))
    found, nmul, sub = recursive_muladd_search!(call, ex.args, mod)
    found || return false, ex
    # found || return ex
    # a, b, c = call.args[2], call.args[3], call.args[4]
    # call.args[2], call.args[3], call.args[4] = c, a, b
    f = if nmul && sub
        :vfnmsub_fast
    elseif nmul
        :vfnmadd_fast
    elseif sub
        :vfmsub_fast
    else
        :vfmadd_fast
    end
    if mod === nothing
        call.args[1] = f
    else
        call.args[1] = Expr(:(.), mod, QuoteNote(f))#_fast))
    end
    true, call
end
function capture_muladd(ex::Expr, mod)
    while true
        ex.head === :ref && return ex
        found, ex = capture_a_muladd(ex, mod)
        found || return ex
    end
end

function append_update_args!(call::Expr, ex::Expr)
    for i ∈ 2:length(ex.args)
        push!(call.args, ex.args[i])
    end
    push!(call.args, ex.args[1])
    nothing
end
contract_pass!(::Any, ::Any) = nothing
function contract!(expr::Expr, ex::Expr, i::Int, mod)
    # if ex.head === :call
        # expr.args[i] = capture_muladd(ex, mod)
    if ex.head === :(+=)
        call = Expr(:call, :add_fast)
        append_update_args!(call, ex)
        expr.args[i] = ex = Expr(:(=), first(ex.args), call)
    elseif ex.head === :(-=)
        call = Expr(:call, :sub_fast)
        append!(call.args, ex.args)
        expr.args[i] = ex = Expr(:(=), first(ex.args), call)
    elseif ex.head === :(*=)
        call = Expr(:call, :mul_fast)
        append_update_args!(call, ex)
        expr.args[i] = ex = Expr(:(=), first(ex.args), call)
    elseif ex.head === :(/=)
        call = Expr(:call, :div_fast)
        append!(call.args, ex.args)
        expr.args[i] = ex = Expr(:(=), first(ex.args), call)
    end
    if ex.head === :(=)
        RHS = ex.args[2]
        # @show ex
        if RHS isa Expr && RHS.head === :call
            ex.args[2] = capture_muladd(RHS, mod)
        end
    end
    contract_pass!(expr.args[i], mod)
end
# contract_pass(x) = x # x will probably be a symbol
function contract_pass!(expr::Expr, mod = nothing)
    for (i,ex) ∈ enumerate(expr.args)
        ex isa Expr || continue
        contract!(expr, ex, i, mod)
    end
end


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), May 23, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames # load package

# ------------------------------------------------------------------------------------------
# ## Getting basic information about a data frame
#
# Let's start by creating a `DataFrame` object, `x`, so that we can learn how to get
# information on that data frame.
# ------------------------------------------------------------------------------------------

x = DataFrame(A = [1, 2], B = [1.0, missing], C = ["a", "b"])

# ------------------------------------------------------------------------------------------
# The standard `size` function works to get dimensions of the `DataFrame`,
# ------------------------------------------------------------------------------------------

size(x), size(x, 1), size(x, 2)

# ------------------------------------------------------------------------------------------
# as well as `nrow` and `ncol` from R; `length` gives number of columns.
# ------------------------------------------------------------------------------------------

nrow(x), ncol(x), length(x)

# ------------------------------------------------------------------------------------------
# `describe` gives basic summary statistics of data in your `DataFrame`.
# ------------------------------------------------------------------------------------------

describe(x)

# ------------------------------------------------------------------------------------------
# Use `showcols` to get informaton about columns stored in a DataFrame.
# ------------------------------------------------------------------------------------------

showcols(x)

# ------------------------------------------------------------------------------------------
# `names` will return the names of all columns,
# ------------------------------------------------------------------------------------------

names(x)

# ------------------------------------------------------------------------------------------
# and `eltypes` returns their types.
# ------------------------------------------------------------------------------------------

eltypes(x)

# ------------------------------------------------------------------------------------------
# Here we create some large DataFrame
# ------------------------------------------------------------------------------------------

y = DataFrame(rand(1:10, 1000, 10));

# ------------------------------------------------------------------------------------------
# and then we can use `head` to peek into its top rows
# ------------------------------------------------------------------------------------------

head(y)

# ------------------------------------------------------------------------------------------
# and `tail` to see its bottom rows.
# ------------------------------------------------------------------------------------------

tail(y, 3)

# ------------------------------------------------------------------------------------------
# ### Most elementary get and set operations
#
# Given the `DataFrame`, `x`, here are three ways to grab one of its columns as a `Vector`:
# ------------------------------------------------------------------------------------------

x[1], x[:A], x[:, 1]

# ------------------------------------------------------------------------------------------
# To grab one row as a DataFrame, we can index as follows.
# ------------------------------------------------------------------------------------------

x[1, :]

# ------------------------------------------------------------------------------------------
# We can grab a single cell or element with the same syntax to grab an element of an array.
# ------------------------------------------------------------------------------------------

x[1, 1]

# ------------------------------------------------------------------------------------------
# Assignment can be done in ranges to a scalar,
# ------------------------------------------------------------------------------------------

x[1:2, 1:2] = 1
x

# ------------------------------------------------------------------------------------------
# to a vector of length equal to the number of assigned rows,
# ------------------------------------------------------------------------------------------

x[1:2, 1:2] = [1,2]
x

# ------------------------------------------------------------------------------------------
# or to another data frame of matching size.
# ------------------------------------------------------------------------------------------

x[1:2, 1:2] = DataFrame([5 6; 7 8])
x

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using BioSequences
function orf()
    jay = 0
    genome = randdnaseq(500)
    orfs = Dict()
    for i = 1:498
        codon = translate(genome[i:i+2])
        if codon == aa"M"
            for j = i:3:498
                codon = translate(genome[j:j+2])
                if codon == aa"*"
                    jay = j
                    break
                end
            end
            length_ = jay - i
            if haskey(orfs, length_)
                orfs[length_]+= 1
            else
                orfs[length_] =1
            end
        end
    end
    println(orfs)
end
orf()
"""
Variables in loops aren't global --> created "jay" to track "j"
Want to look for other BioSequences functions to improve this program
"""

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
using Test, YaoBlocks, YaoArrayRegister

@testset "Yao/#186" begin
    @test getiparams(phase(0.1)) == 0.1
    @test getiparams(2 * phase(0.1)) == ()
end

@testset "block to matrix conversion" begin
    for each in [X, Y, Z, H]
        Matrix{ComplexF64}(each) == Matrix{ComplexF64}(mat(each))
    end
end

@testset "apply lambda" begin
    r = rand_state(3)
    @test apply!(copy(r), put(1=>X)) ≈ apply!(copy(r), put(3, 1=>X))
    f(x::Float32) = x
    @test_throws ErrorException apply!(r, f)
end

# copy return itself by default
@test copy(X) === X

# block type can be used as traits
@test nqubits(typeof(X)) == 1


text:
using MAT

@testset "exampleHBF" begin
    # MATLAB output
    test = matread("resources/exampleHBF.mat")

    # HBFex1
    F1, F2 = exampleHBF(1)
    for i = 1:length(F1)
        @test test["F1_1"]["val"][i] == F1[1][i]
        @test test["F1_1"]["csd"][i] == F1[2][i]
    end
    for i = 1:length(F2)
        @test test["F2_1"]["val"][i] == F2[1][i]
        @test test["F2_1"]["csd"][i] == F2[2][i]
    end

    # HBFex2
    F1, F2 = exampleHBF(2)
    for i = 1:length(F1)
        @test test["F1_2"]["val"][i] == F1[1][i]
        @test test["F1_2"]["csd"][i] == F1[2][i]
    end
    for i = 1:length(F2)
        @test test["F2_2"]["val"][i] == F2[1][i]
        @test test["F2_2"]["csd"][i] == F2[2][i]
    end

    # Saramaki88
    F1, F2 = exampleHBF(3)
    for i = 1:length(F1)
        @test test["F1_3"]["val"][i] == F1[1][i]
        @test test["F1_3"]["csd"][i] == F1[2][i]
    end
    for i = 1:length(F2)
        @test test["F2_3"]["val"][i] == F2[1][i]
        @test test["F2_3"]["csd"][i] == F2[2][i]
    end

    # Saramaki90
    F1, F2 = exampleHBF(4)
    for i = 1:length(F1)
        @test test["F1_4"]["val"][i] == F1[1][i]
        @test test["F1_4"]["csd"][i] == F1[2][i]
    end
    for i = 1:length(F2)
        @test test["F2_4"]["val"][i] == F2[1][i]
        @test test["F2_4"]["csd"][i] == F2[2][i]
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


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

export find_sources

"""
`find_sources([basedir]; recursive = true, exclude_self = true)`

Finds all Julia source (".jl") files in the specified directory (`pwd()`, by
default).

Subdirectories are listed before files, directory and file names are
sorted alphabetically.
"""

function find_sources(
    basedir::AbstractString = pwd();
    recursive::Bool = true, exclude_self::Bool = true
)
    const current_src = current_src_file()

    function impl(subdir::AbstractString)
        const fullpath = joinpath(basedir, subdir)

        const result = Array(AbstractString, 0)

        const files = Array(AbstractString, 0)
        const dirs = Array(AbstractString, 0)
        for f in sort(readdir(fullpath))
            const fname = joinpath(fullpath, f)
            if isdir(fname)
                push!(dirs, f)
            elseif isfile(fname) && endswith(fname, ".jl")
                if (
                    !exclude_self || (
                        # Make double sure to exclude calling file:
                        ( fname != current_src ) &&
                        ( subdir != "" || f != basename(current_src) )
                    )
                )
                    push!(files, joinpath(subdir, f))
                end
            end
        end

        if (recursive)
            for d in dirs
                append!(result, impl(joinpath(subdir, d)))
            end
        end

        append!(result, files)

        result
    end

    impl("")
end


text:
using Base: splat

test = """
    5483143223
    2745854711
    5264556173
    6141336146
    6357385478
    4167524645
    2176841721
    6882881134
    4846848554
    5283751526
    """

function parse_input(input::AbstractString)::Matrix
    return [[parse(Int, c) for c in row]
            for row in split(input, r"\s+"; keepempty = false)] |>
                splat(hcat) |> permutedims
end

function neighbors(M::AbstractMatrix, ci::CartesianIndex)::Vector
    ci1 = CartesianIndex(1, 1)
    ns0 = (ci - ci1):(ci + ci1)
    return setdiff(vec(ns0) ∩ keys(M), Ref(ci))
end

function step!(M::AbstractMatrix)::Int
    M .+= 1
    flash = findall(>(9), M)
    for f in flash
        M[neighbors(M, f)] .+= 1
        union!(flash, findall(>(9), M))
    end
    M[flash] .= 0
    return length(flash)
end

function q1(octopuses::AbstractMatrix)::Int
    oct = copy(octopuses)
    return sum(_ -> step!(oct), 1:100)
end

function q2(octopuses::AbstractMatrix)::Int
    oct = copy(octopuses)
    steps = 0
    while !all(iszero, oct)
        step!(oct)
        steps += 1
    end
    return steps
end

let octopuses = parse_input(test)
    @assert q1(octopuses) == 1656
    @assert q2(octopuses) == 195
end

octopuses = parse_input(read("day11.in", String))
println("Q1: ", q1(octopuses))
println("Q2: ", q2(octopuses))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using Serialization
using Random

Random.seed!(1234)

include("models.jl")
include("integrators.jl")
include("enoc.jl")
using .Models
using .Integrators
using .enoc

M = 60
D = 3
k = 40
k_r = 30

window = 20
da = false
inflation = 1.0

osc_vars = 1:D
modes = 3:4
model = Models.chua_true
model_err = Models.chua_err
integrator = Integrators.rk4
outfreq = 1
Δt = 0.1
m = 20
cycles = 1000
record_length = 25000
ens_err_pct = 0.2
obs_err_pct = 0.1
transient = 3000
mp = 5
y0 = [-1.06095, 0.160678, 0.267729]
varimax = true
check_bounds = true
test_time = 10.0
y_fcst = true
α = 0.5

if !da
    info, ssa_info = enoc.run(model=model, model_err=model_err, M=M, D=D, k=k,
                              k_r=k_r, modes=modes, osc_vars=osc_vars,
                              integrator=integrator, outfreq=outfreq, Δt=Δt,
                              m=m, cycles=1000, window=window,
                              record_length=record_length,
                              ens_err_pct=ens_err_pct, obs_err_pct=obs_err_pct,
                              transient=transient, y0=y0, mp=mp,
                              varimax=varimax, da=false, inflation=inflation,
                              check_bounds=check_bounds, y_fcst=y_fcst, α=α,
                              preload="chua")

    mp = argmin(enoc.optimal_ens(info)[1])[1]
end

Random.seed!(1234)

info, ssa_info = enoc.run(model=model, model_err=model_err, M=M, D=D, k=k,
                          k_r=k_r, modes=modes, osc_vars=osc_vars,
                          integrator=integrator, outfreq=outfreq, Δt=Δt, m=m,
                          cycles=cycles, window=window,
                          record_length=record_length, ens_err_pct=ens_err_pct,
                          obs_err_pct=obs_err_pct, transient=transient, y0=y0,
                          mp=mp, varimax=varimax, da=da, inflation=inflation,
                          check_bounds=check_bounds, y_fcst=y_fcst, α=α,
                          preload="chua")

serialize(open(string("out_chua_", ARGS[1]), "w"), info)


text:
grouped_bar_defaults() = Dict([
    :strokewidth => 0,
]);

"""
	$(SIGNATURES)

Grouped bar graph. Rows of `dataM` are groups.

# Example

```julia
grouped_bar_graph(["1", "2"], rand(2, 3); pos = (1,2), yerror = rand(2, 3))
```
"""
function grouped_bar_graph(
    xStrV, dataM :: AbstractMatrix{F};
    fig = Figure(), pos = (1,1),
    yerror = nothing,
    kwargs...) where F

    ng = length(xStrV);
    @assert size(xStrV, 1) == size(dataM, 1)  "Size mismatch: $(size(xStrV)), $(size(dataM))";
    args = merge(grouped_bar_defaults(), kwargs);
    ax = make_axis(fig, pos;
        xticks = (1 : ng, xStrV), args...);

    grouped_bar_graph!(ax, dataM; yerror = yerror, kwargs...);
    return fig, ax
end

"""
	$(SIGNATURES)

Grouped bar graph for a given axis.
"""
function grouped_bar_graph!(ax :: Axis, dataM :: AbstractMatrix{F};
    yerror = nothing, kwargs...) where F

    if !isnothing(yerror)
        @assert size(yerror) == size(dataM)  "Size mismatch yerror";
    end
    args = merge(grouped_bar_defaults(), kwargs);
    nr, nc = size(dataM);
    xV = repeat(1 : nr; outer = (nc,));
    grps = groups_from_rows(nr, nc);
    bars = bars_from_columns(nr, nc);
    barplot!(ax, xV, vec(dataM);
        dodge = grps, color = get_colors(bars, nc),
        # color = bars, colorrange = (1, nc),
        args...);

    add_grouped_errorbars!(xV, dataM, yerror, grps);
end

# should be improved and be more robust +++
function add_grouped_errorbars!(xV, dataM, yerror, grps)
    crossbar!(xV, vec(dataM), vec(dataM) .- vec(yerror), vec(dataM) .+ vec(yerror);
        dodge = grps, color = (:gray50, 0.35));
end

add_grouped_errorbars!(xV, dataM, yerror :: Nothing, grps) = nothing;

groups_from_rows(nr, nc) = repeat(1 : nc; inner = (nr,));
bars_from_columns(nr, nc) = repeat(1 : nc; inner = (nr,));

"""
	$(SIGNATURES)

Places a legend into Axis `ax`. Uses the same colors as for grouped bar graph.
`title` cannot be "". That throws an error. Useful not just for grouped bar graphs, but also for any kind of subplot.

This needs to be plotted into a blank axis. Otherwise the legend position is odd.

For horizontal legends: `tellheight = true` prevents vertical spacing issues.

Should be improved +++.

# Example

```julia
fig, ax = grouped_bar_graph(["1", "2"], rand(2, 3));
grouped_bar_legend(fig[2, 1], colLabelV; orientation = :horizontal, tellheight = true);
```
"""
function grouped_bar_legend(ax, labelV; title = " ", kwargs...)
    n = length(labelV);
    elements = [PolyElement(polycolor = get_colors(i, n)) for i in 1:n];
    Legend(ax, elements, labelV, title; kwargs...);
end

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


text:
using ImageView
using Images, OffsetArrays, Reactive
using Test

@testset "Simple grayscale" begin
    imshow_now(rand(Gray{N0f8}, 10, 10))
    # clamping checks
    imshow_now(rand(UInt8, 10, 10))  # most/all pixels should be white
    A = randn(10, 10)
    imshow_now(A)
    A[2,2] = NaN
    A[3,3] = -Inf
    A[4,4] = Inf
    imshow_now(A)

    # default contrast setting with a homogenous image
    imgdict = imshow_now(zeros(3, 3))
    @test value(imgdict["clim"]) == ImageView.CLim(0.0,1.0)
end

@testset "Simple RGB" begin
    imshow_now(rand(RGB{N0f8}, 10, 10))
    # clamping checks
    A = randn(3, 10, 10)
    imshow_now(colorview(RGB, A))
    A[1,2,2] = NaN
    A[1,3,3] = -Inf
    A[1,4,4] = Inf
    imshow_now(colorview(RGB, A))

end

@testset "Non-1 indices" begin
    A = OffsetArray(rand(11, 10), -5:5, 0:9)
    ret = imshow_now(A)
end


text:
using MonteCarloTreeSearch
using Documenter

DocMeta.setdocmeta!(MonteCarloTreeSearch, :DocTestSetup, :(using MonteCarloTreeSearch); recursive=true)

makedocs(;
    modules=[MonteCarloTreeSearch],
    authors="Kyu-Young Kim <kykim@cs.stanford.edu>, Zachary Sunberg <zachary.sunberg@colorado.edu>, and contributors",
    repo="https://github.com/JuliaPOMDP/MonteCarloTreeSearch.jl/blob/{commit}{path}#{line}",
    sitename="MonteCarloTreeSearch.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaPOMDP.github.io/MonteCarloTreeSearch.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaPOMDP/MonteCarloTreeSearch.jl",
    devbranch="main",
)