text:
# Rank-based correlations
#
# - Spearman's correlation
# - Kendall's correlation
#

#######################################
#
#   Spearman correlation
#
#######################################

corspearman(x::RealVector, y::RealVector) = cor(tiedrank(x), tiedrank(y))

corspearman(X::RealMatrix, Y::RealMatrix) = cor(mapslices(tiedrank, X, 1), mapslices(tiedrank, Y, 1))
corspearman(X::RealMatrix, y::RealVector) = cor(mapslices(tiedrank, X, 1), tiedrank(y))
corspearman(x::RealVector, Y::RealMatrix) = cor(tiedrank(x), mapslices(tiedrank, Y, 1))

corspearman(X::RealMatrix) = (Z = mapslices(tiedrank, X, 1); cor(Z, Z))

#######################################
#
#   Kendall correlation
#
#######################################

# Knigh JASA (1966)

function corkendall!(x::RealVector, y::RealVector)
    if any(isnan(x)) || any(isnan(y)) return NaN end
    n = length(x)
    if n != length(y) error("Vectors must have same length") end

    # Initial sorting
    pm = sortperm(y)
    x[:] = x[pm]
    y[:] = y[pm]
    pm[:] = sortperm(x)
    x[:] = x[pm]

    # Counting ties in x and y
    iT = 1
    nT = 0
    iU = 1
    nU = 0
    for i = 2:n
        if x[i] == x[i-1]
            iT += 1
        else
            nT += iT*(iT - 1)
            iT = 1
        end
        if y[i] == y[i-1]
            iU += 1
        else
            nU += iU*(iU - 1)
            iU = 1
        end
    end
    if iT > 1 nT += iT*(iT - 1) end
    nT = div(nT,2)
    if iU > 1 nU += iU*(iU - 1) end
    nU = div(nU,2)

    # Sort y after x
    y[:] = y[pm]

    # Calculate double ties
    iV = 1
    nV = 0
    jV = 1
    for i = 2:n
        if x[i] == x[i-1] && y[i] == y[i-1]
            iV += 1
        else
            nV += iV*(iV - 1)
            iV = 1
        end
    end
    if iV > 1 nV += iV*(iV - 1) end
    nV = div(nV,2)

    nD = div(n*(n - 1),2)
    return (nD - nT - nU + nV - 2swaps!(y)) / (sqrt(nD - nT) * sqrt(nD - nU))
end

corkendall(x::RealVector, y::RealVector) = corkendall!(float(copy(x)), float(copy(y)))

corkendall(X::RealMatrix, y::RealVector) = Float64[corkendall!(float(X[:,i]), float(copy(y))) for i in 1:size(X, 2)]

corkendall(x::RealVector, Y::RealMatrix) = (n = size(Y,2); reshape(Float64[corkendall!(float(copy(x)), float(Y[:,i])) for i in 1:n], 1, n))

corkendall(X::RealMatrix, Y::RealMatrix) = Float64[corkendall!(float(X[:,i]), float(Y[:,j])) for i in 1:size(X, 2), j in 1:size(Y, 2)]

function corkendall(X::RealMatrix)
    n = size(X, 2)
    C = eye(n)
    for j = 2:n
        for i = 1:j-1
            C[i,j] = corkendall!(X[:,i],X[:,j])
            C[j,i] = C[i,j]
        end
    end
    return C
end

# Auxilliary functions for Kendall's rank correlation

function swaps!(x::RealVector)
    n = length(x)
    if n == 1 return 0 end
    n2 = div(n, 2)
    xl = sub(x, 1:n2)
    xr = sub(x, n2+1:n)
    nsl = swaps!(xl)
    nsr = swaps!(xr)
    sort!(xl)
    sort!(xr)
    return nsl + nsr + mswaps(xl,xr)
end

function mswaps(x::RealVector, y::RealVector)
    i = 1
    j = 1
    nSwaps = 0
    n = length(x)
    while i <= n && j <= length(y)
        if y[j] < x[i]
            nSwaps += n - i + 1
            j += 1
        else
            i += 1
        end
    end
    return nSwaps
end


text:
# ---
# title: 626. Exchange Seats
# id: problem626
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/exchange-seats/description/>
# hidden: true
# ---
#
# Mary is a teacher in a middle school and she has a table `seat` storing
# students' names and their corresponding seat ids.
#
# The column **id** is continuous increment.
#
#
#
# Mary wants to change seats for the adjacent students.
#
#
#
# Can you write a SQL query to output the result for Mary?
#
#
#
#
#
#     +---------+---------+
#     |    id   | student |
#     +---------+---------+
#     |    1    | Abbot   |
#     |    2    | Doris   |
#     |    3    | Emerson |
#     |    4    | Green   |
#     |    5    | Jeames  |
#     +---------+---------+
#
#
# For the sample input, the output is:
#
#
#
#
#
#     +---------+---------+
#     |    id   | student |
#     +---------+---------+
#     |    1    | Doris   |
#     |    2    | Abbot   |
#     |    3    | Green   |
#     |    4    | Emerson |
#     |    5    | Jeames  |
#     +---------+---------+
#
#
# **Note:**
# If the number of students is odd, there is no need to change the last one's
# seat.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
n_values  = [1, 10, 20, 50, 100, 150, 200]
k_values  = [10, 20, 50, 100, 150, 200]
m_values  = [10, 20, 50, 100, 150, 200]

include("_matmul.jl")
typ = get(ENV, "JULIA_TEST_ELTYPE", "ALL")
types = if typ == "Float64"
  DataType[Float64]
elseif typ == "Float32"
  DataType[Float32]
elseif typ == "Int64"
  DataType[Int64]
elseif typ == "Int32"
  DataType[Int32]
else
  DataType[Float64, Float32, Int64, Int32]
end
testset_name_suffix = "(coverage)"
for T ∈ types
  @time test_complex(T, m_values, k_values, n_values, testset_name_suffix)
  @time test_real(T, m_values, k_values, n_values, testset_name_suffix)
end


text:
using InterpretableModels: slim_scoring_loss01_norm0, column_name_and_index
using DataFrames

@testset "slim_scoring_loss01_norm0" begin
    y = Int[1, -1]
    X = Int[1 0 -1; 0 -1 1]
    λ = [0, 1, 2]
    # sum(y .* (X * λ) .<= 0) but we write it out per element
    @test slim_scoring_loss01_norm0(λ, X, y, 10) ==
        ((y[1] * (X[1, 1] * λ[1] + X[1, 2] * λ[2] + X[1, 3] * λ[3]) <= 0) +
         (y[2] * (X[2, 1] * λ[1] + X[2, 2] * λ[2] + X[2, 3] * λ[3]) <= 0))/2 +
         (10 * sum(λ .!= 0))

    @test slim_scoring_loss01_norm0(λ, X, y, 0.34) ==
         ((y[1] * (X[1, 1] * λ[1] + X[1, 2] * λ[2] + X[1, 3] * λ[3]) <= 0) +
          (y[2] * (X[2, 1] * λ[1] + X[2, 2] * λ[2] + X[2, 3] * λ[3]) <= 0))/2 +
          (0.34 * sum(λ .!= 0))

    @test slim_scoring_loss01_norm0(λ, X, float.(y), 0.34) ==
          slim_scoring_loss01_norm0(λ, X,        y,  0.34)

    @test slim_scoring_loss01_norm0(λ, float.(X), float.(y), 1e-3) ==
          slim_scoring_loss01_norm0(λ,        X,         y,  1e-3)

    @test slim_scoring_loss01_norm0(float.(λ), X, y, 5e-5) ==
          slim_scoring_loss01_norm0(       λ,  X, y, 5e-5)

end

@testset "column_name_and_index" begin
    df = DataFrame(Cl = Int[1,2], F1 = Float64[1.3, 4.5], F2 = Int[1, 3])
    @test column_name_and_index(df, 1) == ("Cl", 1)
    @test column_name_and_index(df, 2) == ("F1", 2)
    @test column_name_and_index(df, 3) == ("F2", 3)

    @test column_name_and_index(df, :Cl) == ("Cl", 1)
    @test column_name_and_index(df, :F2) == ("F2", 3)
    @test_throws AssertionError column_name_and_index(df, :MyClass)

    @test column_name_and_index(df, "Cl") == ("Cl", 1)
    @test column_name_and_index(df, "F2") == ("F2", 3)
    @test_throws AssertionError column_name_and_index(df, "MyF1")

    @test column_name_and_index(df, "Class"[1:2]) == ("Cl", 1)
    @test_throws AssertionError column_name_and_index(df, "Class"[1:3])
end

@testset "SLIMScoringSystem" begin
end


text:
function simula01()
	#
	# simulated model to test the functions
	#
	# simulated mechanical system
	#
	# |    k1    ------     k2    ------
	# |---###---|  m1  |---###---|  m2  |
	# |    c1    ------     c2    ------
	# |           @  @             @  @
	# |========================================
	#
    # javier.cara@upm.es, 2018-04-16

	# time step
	dt = 0.02

	# matrices for the mechanical system M-K-G
	# ----------------------------------------
	# mass matrix
	m1 = 35
	m2 = 17.5
	M = [m1 0;0 m2]
	Minv = [1/m1 0;0 1/m2] # inv(M)

	# stiffness matrix
	k1 = 8750
	k2 = 3500
	K = [k1+k2 -k2;-k2 k2]
	# theoretical natural frecuencies: 1.71 Hz, 3.31 Hz

	# eigenvalues and eigenvectors
	D,V = eigen(K,M)
	W = sqrt.(D) # modal frequencies rad/s

	# modal matrices
	Mm = V'*M*V
	Km = V'*K*V

	# damping ratios
	zm = [0.02, 0.02]

	# modal damping matrix
	Gm = zeros(2,2)
	for i in 1:2
		Gm[i,i] = 2*sqrt(Mm[i,i]*Km[i,i])*zm[i]
	end

	# damping matrix
	G = (V')\Gm/V

	# matrices for the discrete state space model
	# --------------------------------------------------------------
	# Matrix A
	Ac = [zeros(2,2) Matrix{Float64}(I,2,2);-Minv*K -Minv*G] # Continuous Matrix A
	A = exp(Ac*dt) # discrete Matrix A

	# Matrix C (continuous = discrete)
	C = [-Minv*K -Minv*G]

	# Matrix B
	# Js is the input location matrix
	J = Matrix{Float64}(I,2,2)
	Bc = [zeros(2,2);Minv*J] # continuous matrix B
	B = (A - Matrix{Float64}(I,4,4))*inv(Ac)*Bc # discrete matrix B

	# matrix D (continuous = discrete)
	D = Minv*J

	# eigenvectors with maximum component = 1
	V1 = zeros(2,2)
	for j in 1:2
		maxval = V[1,j]
		if abs(V[2,j]) > abs(maxval)
			maxval = V[2,j]
		end
		V1[:,j] = V[:,j]/maxval
	end
	# modal masses corresponding to V1
	Mm1 = V1'*M*V1
	mm1 = diag(Mm1)

	m = Dict("dt"=>dt,"M"=>M,"K"=>K,"G"=>G,
	"wm"=>W,"zm"=>zm,"Vm"=>V1,"mm"=>mm1,
	"A"=>A,"B"=>B,"C"=>C,"D"=>D
	)

	return m
end


text:
cint = BasicLU.cint
cdbl = BasicLU.cdbl
spvector = BasicLU.spvector
spmatrix = BasicLU.spmatrix
BASICLU_OK = BasicLU.BASICLU_OK
BASICLU_UPDATE_COST = BasicLU.BASICLU_UPDATE_COST
BASICLU_NFACTORIZE = BasicLU.BASICLU_NFACTORIZE
BASICLU_NFORREST_TOTAL = BasicLU.BASICLU_NFORREST_TOTAL
BASICLU_NUPDATE_TOTAL = BasicLU.BASICLU_NUPDATE_TOTAL

# =============================================================================
# readmat: read variables from MATLAB file
# =============================================================================

function readmat(file, names...)
    n = length(names)
    objects = Array{Any}(undef,n)
    f1 = matopen(file)
    for k = 1:n
        objects[k] = read(f1, names[k])
    end
    close(f1)
    return tuple(objects...)
end

# =============================================================================
# test_factorize
# =============================================================================
"""
    test_factorize(testdir, trans=false)

For all *.mat files in @testdir read matrix B and factorize it. Monitor residual
of factorization and forward/backward solves.

@trans specifies if B or its transposed are factorized.
"""
function test_factorize(testdir::String, trans=false)
    files = readdir(testdir)
    for f in files
        if length(f) < 4 || f[end-3:end] != ".mat"
            continue
        end
        @printf(" %-24s", f)
        B, = readmat(string(testdir, f), "B")
        Bt = sparse(transpose(B))
        if trans
            (B,Bt) = (Bt,B)
        end
        m = size(B,1)
        blu = initialize(m)
        err = factorize(blu, B)
        if err != BASICLU_OK
            @printf(" failed (%d)\n", err)
            continue
        end
        (L,U,p,q) = get_factors(blu)
        res = norm(L*U-B[p,q], Inf)
        rhs = ones(m)
        lhs = solve(blu, rhs, 'N')
        res = max(res, norm(B*lhs-rhs, Inf))
        lhs = solve(blu, rhs, 'T')
        res = max(res, norm(Bt*lhs-rhs, Inf))
        rhs = spvector(m, [1;m], [1.0;1.0])
        lhs = solve(blu, rhs, 'N')
        res = max(res, norm(B*lhs-rhs, Inf))
        lhs = solve(blu, rhs, 'T')
        res = max(res, norm(Bt*lhs-rhs, Inf))
        lhs = solve4update(blu, rhs, true)
        res = max(res, norm(B*lhs-rhs, Inf))
        lhs = solve4update(blu, m-1, true)
        rhs = zeros(m); rhs[m-1] = 1
        res = max(res, norm(Bt*lhs-rhs, Inf))
        @printf("%.2e\n", res)
    end
    nothing
end

# =============================================================================
# test_update
# =============================================================================
"""
    test_update(testdir)

For all *.mat files in @testdir read matrix A and the pivot sequence. Starting
from the slack basis apply pivot operations. Monitor residuals to
forward/backward solves after each 100 updates.
"""
function test_update(testdir::String)
    files = readdir(testdir)
    for f in files
        if length(f) < 4 || f[end-3:end] != ".mat"
            continue
        end
        @printf(" %-24s", f)
        A,invar,outvar = readmat(string(testdir, f), "A", "invar", "outvar")
        invar = invar[:]
        outvar = outvar[:]
        m,n = size(A)
        A1 = [A sparse(1.0I, m, m)]
        basis = Array{cint}(undef,m)
        basis[:] = collect(1:m) .+ n # slack basis
        res,nfactor,nforrest,nperm = test_update(A1, basis, invar, outvar)
        @printf("%.2e      %5d factor, %6d forrest, %6d perm\n",
                res, nfactor, nforrest, nperm)
    end
    nothing
end

function test_update(A::spmatrix, basis::Array{cint,1}, invar::Array{cint,1},
                     outvar::Array{cint,1})
    m,n = size(A)
    niter = length(invar)
    map2basis = zeros(cint, n)
    map2basis[basis] = 1:m
    blu = initialize(m)
    # blu.xstore[BASICLU_REL_PIVOT_TOLERANCE] = 0.5
    err = factorize(blu, A[:,basis])
    if err != BASICLU_OK
        return (NaN, 0, 0, 0)
    end
    res = 0.0
    for k = 1:niter
        if invar[k] == outvar[k] # bound flip
            continue
        end
        j = invar[k]
        p = map2basis[outvar[k]]
        @assert p > 0
        lhs = solve4update(blu, A[:,j], true)
        xtbl = lhs[p]
        solve4update(blu, p)
        piverr = update(blu, xtbl)
        basis[p] = j
        map2basis[j] = p
        map2basis[outvar[k]] = 0
        cost = blu.xstore[BASICLU_UPDATE_COST]
        if piverr > 1e-10 || cost > 1.0
            err = factorize(blu, A[:,basis])
            if err != BASICLU_OK
                return (NaN, 0, 0, 0)
            end
        end
        if k%100 == 0
            B = A[:,basis]
            rhs = spvector(m, [1;m], [1.0;1.0])
            lhs = solve(blu, rhs, 'N')
            res = max(res, norm(B*lhs-rhs, Inf))
            lhs = solve(blu, rhs, 'T')
            res = max(res, norm(B'*lhs-rhs, Inf))
        end
    end
    nfactor = convert(cint, blu.xstore[BASICLU_NFACTORIZE])
    nforrest = convert(cint, blu.xstore[BASICLU_NFORREST_TOTAL])
    nperm = convert(cint, blu.xstore[BASICLU_NUPDATE_TOTAL]) - nforrest
    return (res, nfactor, nforrest, nperm)
end


text:
using Test
using Printf: @printf
using QuantumLattices.Mathematics.AlgebraOverFields
using QuantumLattices.Interfaces: dimension, rank, add!, sub!, mul!, div!, sequence
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Mathematics.Combinatorics: Combinations
using QuantumLattices.Mathematics.VectorSpaces: SimpleVectorSpace
import QuantumLattices.Interfaces: ⊗, ⋅, permute

struct AOFID{O<:Real, S<:Real} <: SimpleID
    orbital::O
    nambu::S
end

@testset "ID" begin
    @test ID{SimpleID, 1}|>rank == 1
    @test ID|>rank == Any
    @test promote_type(ID{AOFID{Int, Int}, 1}, ID{AOFID{Int, Int}}) == ID{AOFID{Int, Int}}
    @test promote_type(ID{AOFID{Int, Int}}, ID{AOFID{Int, Int}, 1}) == ID{AOFID{Int, Int}}
    @test promote_type(ID{AOFID{Int, Int}}, ID{AOFID{Int, Float}}) == ID{AOFID{Int, <:Real}}
    @test promote_type(Tuple{}, ID{AOFID{Int, Int}}) == ID{AOFID{Int, Int}}
    @test promote_type(Tuple{}, ID{AOFID{Int, Int}, 2}) == ID{AOFID{Int, Int}}

    cid = AOFID(2, 1) * AOFID(1, Inf)
    @test cid == ID(AOFID(2, 1), AOFID(1, Inf))
    @test cid|>string == "ID(AOFID(2, 1), AOFID(1, Inf))"
    @test cid|>eltype == AOFID{Int, <:Real}
    @test cid|>rank == 2
    @test cid|>typeof|>rank == 2

    sid = AOFID(1, 1)
    @test sid * cid == ID(AOFID(1, 1), AOFID(2, 1), AOFID(1, Inf))
    @test cid * sid == ID(AOFID(2, 1), AOFID(1, Inf), AOFID(1, 1))
    @test cid * cid == ID(AOFID(2, 1), AOFID(1, Inf), AOFID(2, 1), AOFID(1, Inf))

    @test deepcopy(ID(AOFID(1, 1))) == ID(AOFID(1, 1))
    @test idisless(ID(AOFID(1, 2)), ID(AOFID(1, Inf)))
    @test idisless(ID(AOFID(2, 1)), ID(AOFID(1, 2), AOFID(1, Inf)))

    cid = AOFID(2, 1) * AOFID(3, 4.0)
    @test cid == ID(AOFID, (2, 3), (1, 4.0))
    @test cid|>typeof|>idpropertynames == (:orbitals, :nambus)
    @test cid.orbitals == (2, 3)
    @test cid.nambus == (1, 4.0)
end

struct AOFOperator{V<:Number, I<:ID{SimpleID}} <: Element{V, I}
    value::V
    id::I
end
⊗(m1::AOFOperator, m2::AOFOperator) = m1 * m2
⋅(m1::AOFOperator, m2::AOFOperator) = m1 * m2
Base.show(io::IO, opt::AOFOperator) = @printf io "AOFOperator(value=%s, id=%s)" opt.value opt.id
Base.repr(opt::AOFOperator) = "AOFOperator($(opt.value), $(opt.id))"

function permute(::Type{<:AOFOperator}, id1::AOFID, id2::AOFID, ::Any = nothing)
    @assert id1 ≠ id2 "permute error: permuted ids should not be equal to each other."
    if (id1.nambu == 3-id2.nambu) && (id1.orbital == id2.orbital)
        if id1.nambu == 2
            return (AOFOperator(1, ID()), AOFOperator(1, ID(id2, id1)))
        else
            return (AOFOperator(-1, ID()), AOFOperator(1, ID(id2, id1)))
        end
    else
        return (AOFOperator(1, ID(id2, id1)),)
    end
end

@testset "Elements" begin
    @test valtype(Element) == Any
    @test valtype(Element{Int}) == Int
    @test idtype(Element{<:Number}) == ID{SimpleID}
    @test idtype(Element{<:Number, ID{AOFID}}) == ID{AOFID}
    @test scalartype(AOFOperator{Number}) == AOFOperator{Number, Tuple{}}
    @test promote_type(AOFOperator{Int}, AOFOperator) == AOFOperator
    @test promote_type(AOFOperator, AOFOperator{Int}) == AOFOperator
    @test promote_type(AOFOperator{Int}, AOFOperator{Float}) == AOFOperator{Float}
    @test promote_type(AOFOperator{Int, ID{AOFID{Int, Int}, 2}}, AOFOperator{Float, ID{AOFID{Int, Int}, 2}}) == AOFOperator{Float, ID{AOFID{Int, Int}, 2}}
    @test promote_type(AOFOperator{Int, ID{AOFID}}, AOFOperator{Float, ID{AOFID}}) == AOFOperator{Float, <:ID{AOFID}}

    opt = AOFOperator(1.0, ID(AOFID(1, 1)))
    @test opt|>deepcopy == opt
    @test isequal(opt|>deepcopy, opt)
    @test isapprox(opt, replace(opt, value=opt.value+10^-6); atol=10^-5)
    @test opt|>valtype == opt|>typeof|>valtype == Float
    @test opt|>idtype == opt|>typeof|>idtype == ID{AOFID{Int, Int}, 1}
    @test opt|>rank == opt|>typeof|>rank == 1
    @test opt|>typeof|>one == AOFOperator(1.0, ID())
    @test +opt == opt
    @test opt+zero(Elements) == zero(Elements)+opt == opt-zero(Elements)
    @test -opt == AOFOperator(-1.0, ID(AOFID(1, 1)))
    @test zero(Elements)-opt == -opt
    @test opt*2 == 2*opt == opt*AOFOperator(2, ID()) == AOFOperator(2, ID())*opt == AOFOperator(2.0, ID(AOFID(1, 1)))
    @test opt*zero(Elements) == zero(Elements)*opt == zero(Elements) == nothing
    @test opt/2 == opt/AOFOperator(2, ID()) == AOFOperator(0.5, ID(AOFID(1, 1)))
    @test opt^2 == opt*opt

    @test AOFOperator(1, ID())+1 == 1+AOFOperator(1, ID()) == AOFOperator(1, ID())+AOFOperator(1, ID()) == AOFOperator(2, ID())
    @test opt+1 == 1+opt == opt+AOFOperator(1, ID()) == AOFOperator(1, ID())+opt
    @test AOFOperator(1, ID())-2 == 1-AOFOperator(2, ID()) == AOFOperator(1, ID())-AOFOperator(2, ID()) == AOFOperator(-1, ID())
    @test AOFOperator(2, ID())*AOFOperator(3, ID()) == AOFOperator(6, ID())

    opt1 = AOFOperator(1.0, ID(AOFID(1, 1)))
    opt2 = AOFOperator(2.0, ID(AOFID(1, 2)))
    opts = Elements(opt1.id=>opt1, opt2.id=>opt2)
    @test string(opts) == "Elements with 2 entries:\n  AOFOperator(value=2.0, id=ID(AOFID(1, 2)))\n  AOFOperator(value=1.0, id=ID(AOFID(1, 1)))\n"
    @test repr(opts) == "Elements with 2 entries:\n  AOFOperator(2.0, ID(AOFID(1, 2)))\n  AOFOperator(1.0, ID(AOFID(1, 1)))"
    @test opts|>deepcopy == opts
    @test opts|>zero == opts|>typeof|>zero == nothing
    @test add!(deepcopy(opts)) == opts
    @test add!(deepcopy(opts), zero(Elements)) == opts
    @test sub!(deepcopy(opts)) == opts
    @test sub!(deepcopy(opts), zero(Elements)) == opts
    @test mul!(deepcopy(opts), 2.0) == mul!(deepcopy(opts), AOFOperator(2, ID())) == opts*2 == opts*AOFOperator(2, ID())
    @test div!(deepcopy(opts), 2.0) == div!(deepcopy(opts), AOFOperator(2, ID())) == opts/2 == opts/AOFOperator(2, ID())
    @test +opts == opts+zero(Elements) == zero(Elements)+opts == opts == opts-zero(Elements)
    @test opt1+opt2 == opts
    @test opt1+opts == opts+opt1 == Elements(opt1*2, opt2)
    @test opts+opts == Elements(opt1*2, opt2*2)
    @test -opts == Elements(-opt1, -opt2) == zero(Elements)-opts
    @test 1.0-opts == AOFOperator(1.0, ID())-opts == -opts+1.0
    @test opts-1.0 == opts-AOFOperator(1.0, ID())
    @test opt-1.0 == opt-AOFOperator(1.0, ID()) == opt+AOFOperator(-1.0, ID())
    @test 1.0-opt == AOFOperator(1.0, ID())-opt == 1.0+(-opt)
    @test zero(Elements) == opt1-opt1 == opts-opts
    @test opts-opts == zero(Elements)
    @test isequal(zero(Elements), opts-opts)
    @test isequal(opts-opts, zero(Elements))
    @test opts-opt1 == Elements(opt2)
    @test opt1-opts == Elements(-opt2)
    @test 2*opts == opts*2 == Elements(opt1*2, opt2*2) == AOFOperator(2, ID())*opts
    @test opts/2 == Elements(opt1/2, opt2/2)
    @test opts^2 == opts*opts
    @test opts*zero(Elements) == zero(Elements)*opts == nothing

    temp = Elements{ID{AOFID{Int, Int}}, AOFOperator{Float}}(opts)
    @test add!(deepcopy(temp), 1) == add!(deepcopy(temp), AOFOperator(1, ID())) == opts+1 == 1+opts == opts+AOFOperator(1, ID()) == AOFOperator(1, ID())+opts
    @test sub!(deepcopy(temp), 1) == sub!(deepcopy(temp), AOFOperator(1, ID())) == opts-1 == opts-AOFOperator(1, ID())

    opt3 = AOFOperator(3, ID(AOFID(1, 3)))
    @test opt1*opt2 == opt1⊗opt2 == opt1⋅opt2 == AOFOperator(2.0, ID(AOFID(1, 1), AOFID(1, 2)))
    @test opts*opt3 == opts⊗opt3 == opts⋅opt3 == Elements(AOFOperator(3.0, ID(AOFID(1, 1), AOFID(1, 3))), AOFOperator(6.0, ID(AOFID(1, 2), AOFID(1, 3))))
    @test opt3*opts == opt3⊗opts == opt3⋅opts == Elements(AOFOperator(3.0, ID(AOFID(1, 3), AOFID(1, 1))), AOFOperator(6.0, ID(AOFID(1, 3), AOFID(1, 2))))
    @test opts*Elements(opt3) == opts⊗Elements(opt3) == opts⋅Elements(opt3) == Elements(AOFOperator(3.0, ID(AOFID(1, 1), AOFID(1, 3))), AOFOperator(6.0, ID(AOFID(1, 2), AOFID(1, 3))))

    table = Dict(AOFID(1, 1)=>1, AOFID(1, 2)=>2)
    opt = AOFOperator(2.0, ID(AOFID(1, 1), AOFID(1, 2)))
    @test sequence(opt, table) == (1, 2)
    @test split(opt) == (2.0, AOFOperator(1.0, ID(AOFID(1, 1))), AOFOperator(1.0, ID(AOFID(1, 2))))

    opt1 = AOFOperator(1, ID(AOFID(1, 1)))
    opt2 = AOFOperator(2, ID(AOFID(1, 2)))
    opts = Elements(opt1.id=>opt1, opt2.id=>opt2)
    opt = AOFOperator(3, ID())
    @test opt1//3 == opt1//opt == AOFOperator(1//3, ID(AOFID(1, 1)))
    @test opt2//3 == opt2//opt == AOFOperator(2//3, ID(AOFID(1, 2)))
    @test opts//3 == opts//opt == Elements(opt1.id=>AOFOperator(1//3, ID(AOFID(1, 1))), opt2.id=>AOFOperator(2//3, ID(AOFID(1, 2))))
end

@testset "replace" begin
    id1, id2 = AOFID(1, 1), AOFID(1, 2)
    opt = AOFOperator(1.5, ID(id1, id2))
    id3, id4 = AOFID(2, 1), AOFID(2, 2)
    opt1 = AOFOperator(2.0, ID(id3))+AOFOperator(3.0, ID(id4))
    opt2 = AOFOperator(2.0, ID(id3))-AOFOperator(3.0, ID(id4))
    @test replace(opt, id1=>opt1, id2=>opt2) == 1.5*opt1*opt2
    @test replace(Elements(opt), id2=>opt2, id1=>opt1) == 1.5*opt1*opt2
end

@testset "permute" begin
    id1, id2 = AOFID(1, 1), AOFID(1, 2)
    opt1 = AOFOperator(1.5, ID(id1, id2))
    opt2 = AOFOperator(1.5, ID(id2, id1))
    opt0 = AOFOperator(1.5, ID())
    @test permute(opt1, Dict(id1=>1, id2=>2)) == opt2-opt0
    @test permute(opt1+opt2, Dict(id1=>1, id2=>2)) == 2*opt2-opt0
    @test permute(opt2, Dict(id1=>2, id2=>1)) == opt1+opt0
    @test permute(opt1+opt2, Dict(id1=>2, id2=>1)) == 2*opt1+opt0
end


text:
A = rand(3, 3)
B = ones(4, 4)
C = randn(5, 6)
K = A ⊗ B
X = collect(K)
v = rand(12)

K3 = A ⊗ B ⊗ C

@testset "vectrick" begin
    @testset "Vec trick" begin

        @test K * v ≈ X * v

        V = rand(4, 3)
        @test K * vec(V) ≈ X * vec(V)
        @test_throws DimensionMismatch K * V
        @test_throws DimensionMismatch K * reshape(V, 2, 6)

        v3 = randn(size(K3, 2))
        @test K3 * v3 ≈ collect(K3) * v3
        u = similar(v)
        @test mul!(u, K, v) ≈ X * v
    end

    @testset "Reshaped vec trick" begin

        @test K * v ≈ X * v

        V = sprand(4, 3, 1.0)
        @test K * vec(V) ≈ X * vec(V)
        @test_throws DimensionMismatch K * V
        @test_throws DimensionMismatch K * reshape(V, 2, 6)
        K3 = A ⊗ B ⊗ C
        v3 = sprand(size(K3, 2),1.0)
        @test K3 * vec(v3) ≈ collect(K3) * v3
        u = similar(v)
        @test mul!(u, K, v) ≈ X * v

    end

    @testset "AbstractKroneckerProduct * AbstractMatrix" begin
        rng = MersenneTwister(123456)
        a, b, x = randn(rng, 30, 20), randn(rng, 40, 50), randn(rng, 1000, 1100)

        @test kron(a, Diagonal(ones(50))) * x ≈ Kronecker.kron_a_id(a, x)
        @test kron(Diagonal(ones(20)), b) * x ≈ Kronecker.kron_id_a(b, x)
        @test kron(a, b) * x ≈ Kronecker.kron_a_b(a, b, x)

        @test kron(a, b) * x ≈ (a ⊗ b) * x
        @test kron(a, Eye(50)) * x ≈ (a ⊗ Eye(50)) * x
        @test kron(Eye(20), b) * x ≈ (Eye(20) ⊗ b) * x
    end

    @testset "sum" begin
        @test sum(K) ≈ sum(X)
        @test sum(K3) ≈ sum(collect(K3))

        @test sum(K, dims=1) ≈ sum(X, dims=1)
        @test sum(K3, dims=2) ≈ sum(collect(K3), dims=2)
        @test sum(K3, dims=2) isa AbstractKroneckerProduct

        @test sum(kronecker(A, 3)) ≈ sum(kron(A, A, A))
    end
end


text:
"""
    transcription_model(model::InfiniteOpt.InfiniteModel)::JuMP.Model

Return the transcription model stored in `model` if that is what is stored in
`model.optimizer_model`.
"""
function transcription_model(model::InfiniteOpt.InfiniteModel)::JuMP.Model
    if !is_transcription_model(InfiniteOpt.optimizer_model(model))
        error("The model does not contain a transcription model.")
    end
    return InfiniteOpt.optimizer_model(model)
end

"""
    InfiniteOpt.build_optimizer_model!(model::InfiniteOpt.InfiniteModel,
                                       key::Val{:TransData})

Transcribe `model` and store it as a `TranscriptionModel` in the
`model.optimizer_model` field which can be accessed with `transcription_model`.
"""
function InfiniteOpt.build_optimizer_model!(model::InfiniteOpt.InfiniteModel,
                                            key::Val{:TransData})
    # clear the optimzier model contents
    trans_model = InfiniteOpt.clear_optimizer_model_build!(model)
    # build the transcription model based on model
    _build_transcription_model!(trans_model, model)
    # update the optimizer model status
    InfiniteOpt.set_optimizer_model_ready(model, true)
    return
end


text:
using Test

    @testset "MyExample.jl" begin
        # Write your tests here.
        @test my_f(2,1) == 5
        @test my_f(2,4) == 8
        @test my_f(3,1) == 7

end


text:
using Distributions

abstract(::AType{typeof(rand)}) = Float64
abstract(::AType{typeof(randn)}) = Float64
abstract(::AType{typeof(rand)}, ::AType{<:Type{Bool}}) where T = Bool

abstract(::AType{typeof(rand)}, T::AType{<:Distribution}) =
  Core.Compiler.return_type(rand, Tuple{widen(T)})


text:
using XGPaint
using Healpix

## Load halos from HDF5 files, establish a CIB model and cosmology
halo_pos, halo_mass = read_halo_catalog_hdf5(
    ENV["SCRATCH"] * "/websky_halos-light.hdf5")
# @load "/home/zequnl/testsamp.jld2" halo_mass halo_pos

cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0)
model = CIB_Planck2013{Float32}()

## Allocate some arrays and file them up for centrals and satellites
@time sources = generate_sources(model, cosmo, halo_pos, halo_mass);

## Deposit the sources into maps
fluxes_cen = Array{Float32, 1}(undef, sources.N_cen)
fluxes_sat = Array{Float32, 1}(undef, sources.N_sat)
m = Map{Float64,RingOrder}(model.nside)

for freq in [
    "18.7", "21.6",
    "24.5", "27.3", "30.0", "35.9",
    "41.7", "44.0", "47.4", "63.9", "67.8",
    "70.0", "73.7", "79.6", "90.2", "100", "111",
    "129", "143", "153", "164", "189", "210", "217",
    "232", "256", "275", "294", "306", "314", "340",
    "353", "375", "409", "467", "525", "545", "584",
    "643", "729", "817", "857", "906", "994", "1080"]

    @time begin
        XGPaint.paint!(m, parse(Float32, freq) * 1.0f9, model, sources,
            fluxes_cen, fluxes_sat)
        Healpix.saveToFITS(m, "!/global/cscratch1/sd/xzackli/cib/cib$(freq).fits")
    end
end


text:
@testset "Ranges" begin
    # @test range(half(1), half(2), length=3)[2]

    @testset "UnitRange" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test UnitRange{$T}(1//2, 5//1) isa UnitRange{$T}
            @eval @test UnitRange($T(1/2), $T(5)) isa UnitRange{$T}
            @eval @test $T(1/2):$T(5) isa UnitRange{$T}
            @eval @test $T(1/2):big(5) isa UnitRange{BigHalfInt}
            @eval @test first($T(1/2):$T(5)) == $T(1/2)
            @eval @test last($T(1/2):$T(5)) == $T(9/2)
            @eval @test length($T(1/2):$T(5)) == 5
            @eval @test length($T(1/2):$T(0)) == 0
            @eval @test length($T(1/2):$T(1/2)) == 1
            @eval @test length($T(1/2):$T(1)) == 1
            @eval @test $T(5/2) ∉ $T(1):$T(5)
            @eval @test $T(5/2) ∈ $T(1/2):$T(5)
            @eval @test $T(2) ∈ $T(1):$T(5)
            @eval @test $T(2) ∉ $T(1/2):$T(5)
            @eval @test 2.5 ∉ $T(1):$T(5)
            @eval @test 2.5 ∈ $T(1/2):$T(5)
            @eval @test 5//2 ∉ $T(1):$T(5)
            @eval @test 5//2 ∈ $T(1/2):$T(5)
            @eval @test 2 ∈ $T(1):$T(5)
            @eval @test 2 ∉ $T(1/2):$T(5)

            @eval @test any(isinteger, $T(1):$T(5))
            @eval @test !any(isinteger, $T(1):$T(0))
            @eval @test !any(isinteger, $T(1/2):$T(5))
            @eval @test !any(isinteger, $T(3/2):$T(0))
            @eval @test !any(!isinteger, $T(1):$T(5))
            @eval @test !any(!isinteger, $T(1):$T(0))
            @eval @test any(!isinteger, $T(1/2):$T(5))
            @eval @test !any(!isinteger, $T(3/2):$T(0))

            @eval @test all(isinteger, $T(1):$T(5))
            @eval @test all(isinteger, $T(1):$T(0))
            @eval @test !all(isinteger, $T(1/2):$T(5))
            @eval @test all(isinteger, $T(3/2):$T(0))
            @eval @test !all(!isinteger, $T(1):$T(5))
            @eval @test all(!isinteger, $T(1):$T(0))
            @eval @test all(!isinteger, $T(1/2):$T(5))
            @eval @test all(!isinteger, $T(3/2):$T(0))
        end
        for T in (halfinttypes..., :BigHalfInt)
            @eval @test first($T(-1/2):$T(1)) == $T(-1/2)
            @eval @test last($T(-1/2):$T(1)) == $T(1/2)
            @eval @test length($T(-1/2):$T(1)) == 2
        end
    end

    @testset "StepRange" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test StepRange{$T,$T}($T(2), $T(1/2), $T(5)) isa StepRange{$T,$T}
            @eval @test StepRange{$T,Int8}($T(2), Int8(2), $T(9/2)) isa StepRange{$T,Int8}
            @eval @test StepRange($T(2), $T(1/2), $T(5)) isa StepRange{$T,$T}
            @eval @test StepRange($T(2), Int8(2), $T(9/2)) isa StepRange{$T,Int8}
            @eval @test StepRange($T(2), big(2), $T(9/2)) isa StepRange{$T,BigInt}

            @eval @test $T(2):$T(1/2):$T(5) isa StepRange{$T,$T}
            @eval @test $T(2):$T(2):$T(9/2) isa StepRange{$T,$T}
            @eval @test $T(2):big(2):$T(9/2) isa StepRange{BigHalfInt,BigInt}
            @eval @test big(2):2:$T(9/2) isa StepRange{BigHalfInt,Int}

            @eval @test $T(2):$T(2):$T(9/2) == $T(2):$T(2):$T(4)
            @eval @test first($T(2):$T(2):$T(9/2)) == 2
            @eval @test step($T(2):$T(2):$T(9/2)) == 2
            @eval @test last($T(2):$T(2):$T(9/2)) == 4
            @eval @test length($T(2):$T(2):$T(9/2)) == 2

            @eval @test first($T(2):$T(1/2):$T(5)) == 2
            @eval @test step($T(2):$T(1/2):$T(5)) == $T(1/2)
            @eval @test last($T(2):$T(1/2):$T(5)) == 5
            @eval @test length($T(2):$T(1/2):$T(5)) == 7

            @eval @test $T(5/2) ∈ $T(1/2):$T(2):$T(7)
            @eval @test $T(5/2) ∉ $T(1/2):$T(5/2):$T(7)
            @eval @test $T(2) ∈ $T(0):$T(1):$T(5/2)
            @eval @test $T(2) ∉ $T(1/2):$T(1):$T(5/2)
            @eval @test 2.5 ∈ $T(1/2):$T(2):$T(7)
            @eval @test 2.5 ∉ $T(1/2):$T(5/2):$T(7)
            @eval @test 5//2 ∈ $T(1/2):$T(2):$T(7)
            @eval @test 5//2 ∈ $T(1/2):2:$T(7)
            @eval @test 5//2 ∈ $T(1/2):2:7
            @eval @test 5//2 ∉ $T(1/2):$T(5/2):$T(7)
            @eval @test 2 ∈ $T(0):$T(1):$T(5/2)
            @eval @test 2 ∉ $T(1/2):$T(1):$T(5/2)
        end
        for T in (halfinttypes..., :BigHalfInt)
            @eval @test $T(2):$T(-1):$T(1/2) == $T(2):$T(-1):$T(1)
            @eval @test first($T(2):$T(-1):$T(1/2)) == 2
            @eval @test step($T(2):$T(-1):$T(1/2)) == -1
            @eval @test last($T(2):$T(-1):$T(1/2)) == 1
            @eval @test length($T(2):$T(-1):$T(1/2)) == 2

            @eval @test $T(5/2) ∉ $T(1/2):$T(-1):$T(-1/2)
            @eval @test 2.5 ∉ $T(1/2):$T(-1):$T(-1/2)
            @eval @test 5//2 ∉ $T(1/2):$T(-1):$T(-1/2)
            @eval @test 5//2 ∉ $T(1/2):-1:$T(-1/2)
        end
    end

    @testset "StepRangeLen" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test StepRangeLen{$T,$T,$T}($T(2), $T(1/2), 4) isa StepRangeLen{$T,$T,$T}
            @eval @test StepRangeLen{$T,$T,$T}($T(2), $T(1/2), 4, 3) isa StepRangeLen{$T,$T,$T}
            @eval @test StepRangeLen{$T,$T,Int8}($T(9/2), Int8(2), 4) isa StepRangeLen{$T,$T,Int8}
            @eval @test StepRangeLen{$T,$T,Int8}($T(9/2), Int8(2), 4, 3) isa StepRangeLen{$T,$T,Int8}
            @eval @test StepRangeLen{$T,Int,Int}(8, 2, 4) isa StepRangeLen{$T,Int,Int}
            @eval @test StepRangeLen{$T,Int,Int}(8, 2, 4, 3) isa StepRangeLen{$T,Int,Int}
            @eval @test StepRangeLen{$T}($T(2), $T(1/2), 4) isa StepRangeLen{$T,$T,$T}
            @eval @test StepRangeLen{$T}($T(2), $T(1/2), 4, 3) isa StepRangeLen{$T,$T,$T}
            @eval @test StepRangeLen{$T}($T(9/2), Int8(2), 4) isa StepRangeLen{$T,$T,Int8}
            @eval @test StepRangeLen{$T}($T(9/2), Int8(2), 4, 3) isa StepRangeLen{$T,$T,Int8}
            @eval @test StepRangeLen{$T}(8, 2, 4) isa StepRangeLen{$T,Int,Int}
            @eval @test StepRangeLen{$T}(8, 2, 4, 3) isa StepRangeLen{$T,Int,Int}
            @eval @test_skip StepRangeLen($T(2), $T(1/2), 4) isa StepRangeLen{$T,$T,$T}
            @eval @test_skip StepRangeLen($T(2), $T(1/2), 4, 3) isa StepRangeLen{$T,$T,$T}
            @eval @test_skip StepRangeLen($T(9/2), Int8(2), 4) isa StepRangeLen{$T,$T,Int8}
            @eval @test_skip StepRangeLen($T(9/2), Int8(2), 4, 3) isa StepRangeLen{$T,$T,Int8}
            @eval @test StepRangeLen($T(9/2), big(2), 4) isa StepRangeLen{BigHalfInt,$T,BigInt}
            @eval @test StepRangeLen($T(9/2), big(2), 4, 3) isa StepRangeLen{BigHalfInt,$T,BigInt}

            @eval @test StepRangeLen{$T}($T(2), $T(1/2), 5) == $T[2, 5/2, 3, 7/2, 4]
            @eval @test first(StepRangeLen{$T}($T(2), $T(1/2), 5)) == $T(2)
            @eval @test step(StepRangeLen{$T}($T(2), $T(1/2), 5)) == $T(1/2)
            @eval @test last(StepRangeLen{$T}($T(2), $T(1/2), 5)) == $T(4)
            @eval @test length(StepRangeLen{$T}($T(2), $T(1/2), 5)) == 5

            @eval @test StepRangeLen{$T}($T(2), $T(1/2), 5, 3) == $T[1, 3/2, 2, 5/2, 3]
            @eval @test first(StepRangeLen{$T}($T(2), $T(1/2), 5, 3)) == $T(1)
            @eval @test step(StepRangeLen{$T}($T(2), $T(1/2), 5, 3)) == $T(1/2)
            @eval @test last(StepRangeLen{$T}($T(2), $T(1/2), 5, 3)) == $T(3)
            @eval @test length(StepRangeLen{$T}($T(2), $T(1/2), 5, 3)) == 5
        end
    end

    @testset "LinRange" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test LinRange{$T}(1/2, 3/2, 3) isa LinRange{$T}
            @eval @test LinRange($T(1/2), $T(3/2), 3) isa LinRange{float($T)}

            @eval @test LinRange{$T}(1/2, 3/2, 3) == $T[1/2, 1, 3/2]
            @eval @test first(LinRange{$T}(1/2, 3/2, 3)) == $T(1/2)
            @eval @test step(LinRange{$T}(1/2, 3/2, 3)) == 0.5
            @eval @test last(LinRange{$T}(1/2, 3/2, 3)) == $T(3/2)
            @eval @test length(LinRange{$T}(1/2, 3/2, 3)) == 3
        end
    end

    @testset "range" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test @inferred(range($T(1/2), stop=$T(5))) isa UnitRange{$T}
            @eval @test @inferred(range($T(1/2), stop=Int8(5))) isa UnitRange{$T}
            @eval @test @inferred(range(Int8(1), stop=$T(7/2))) isa UnitRange{$T}
            @eval @test @inferred(range($T(1/2), stop=big(5))) isa UnitRange{BigHalfInt}
            @eval @test @inferred(range(big(1), stop=$T(7/2))) isa UnitRange{BigHalfInt}
            @eval @test range($T(1/2), stop=$T(5)) == $T[1/2, 3/2, 5/2, 7/2, 9/2]
            @eval @test range($T(1/2), stop=Int8(5)) == $T[1/2, 3/2, 5/2, 7/2, 9/2]
            @eval @test range(Int8(1), stop=$T(7/2)) == $T[1, 2, 3]
            @eval @test range($T(1/2), stop=big(5)) == BigHalfInt[1/2, 3/2, 5/2, 7/2, 9/2]
            @eval @test range(big(1), stop=$T(7/2)) == BigHalfInt[1, 2, 3]
            @static if VERSION ≥ v"1.7.0-DEV.263"
                @eval @test @inferred(range(start=$T(1/2), stop=$T(5))) isa UnitRange{$T}
                @eval @test @inferred(range(start=$T(1/2), stop=Int8(5))) isa UnitRange{$T}
                @eval @test @inferred(range(start=Int8(1), stop=$T(7/2))) isa UnitRange{$T}
                @eval @test @inferred(range(start=$T(1/2), stop=big(5))) isa UnitRange{BigHalfInt}
                @eval @test @inferred(range(start=big(1), stop=$T(7/2))) isa UnitRange{BigHalfInt}
                @eval @test range(start=$T(1/2), stop=$T(5)) == $T[1/2, 3/2, 5/2, 7/2, 9/2]
                @eval @test range(start=$T(1/2), stop=Int8(5)) == $T[1/2, 3/2, 5/2, 7/2, 9/2]
                @eval @test range(start=Int8(1), stop=$T(7/2)) == $T[1, 2, 3]
                @eval @test range(start=$T(1/2), stop=big(5)) == BigHalfInt[1/2, 3/2, 5/2, 7/2, 9/2]
                @eval @test range(start=big(1), stop=$T(7/2)) == BigHalfInt[1, 2, 3]
            end

            @eval @test @inferred(range($T(1/2), length=5)) isa UnitRange{$T}
            @eval @test range($T(1/2), length=5) == $T[1/2, 3/2, 5/2, 7/2, 9/2]
            @static if VERSION ≥ v"1.7.0-DEV.263"
                @eval @test @inferred(range(start=$T(1/2), length=5)) isa UnitRange{$T}
                @eval @test range(start=$T(1/2), length=5) == $T[1/2, 3/2, 5/2, 7/2, 9/2]

                @static if VERSION ≥ v"1.7.0-DEV.349"
                    @eval @test @inferred(range(stop=$T(11/2), length=5)) isa UnitRange{$T}
                end
                @eval @test range(stop=$T(11/2), length=5) == $T[3/2, 5/2, 7/2, 9/2, 11/2]
            end

            @eval @test @inferred(range($T(2), step=$T(1/2), stop=$T(5))) isa StepRange{$T,$T}
            @eval @test @inferred(range($T(2), step=Int8(2), stop=$T(9/2))) isa StepRange{$T,Int8}
            @eval @test @inferred(range($T(2), step=big(2), stop=$T(9/2))) isa StepRange{BigHalfInt,BigInt}
            @eval @test @inferred(range(Int8(2), step=$T(1/2), stop=Int8(5))) isa StepRange{$T,$T}
            @eval @test range($T(2), step=$T(1/2), stop=$T(5)) == $T[2, 5/2, 3, 7/2, 4, 9/2, 5]
            @eval @test range($T(2), step=Int8(2), stop=$T(9/2)) == $T[2, 4]
            @eval @test range($T(2), step=big(2), stop=$T(9/2)) == BigHalfInt[2, 4]
            @eval @test range(Int8(2), step=$T(1/2), stop=Int8(5)) == $T[2, 5/2, 3, 7/2, 4, 9/2, 5]
            @static if VERSION ≥ v"1.7.0-DEV.263"
                @eval @test @inferred(range(start=$T(2), step=$T(1/2), stop=$T(5))) isa StepRange{$T,$T}
                @eval @test @inferred(range(start=$T(2), step=Int8(2), stop=$T(9/2))) isa StepRange{$T,Int8}
                @eval @test @inferred(range(start=$T(2), step=big(2), stop=$T(9/2))) isa StepRange{BigHalfInt,BigInt}
                @eval @test @inferred(range(start=Int8(2), step=$T(1/2), stop=Int8(5))) isa StepRange{$T,$T}
                @eval @test range(start=$T(2), step=$T(1/2), stop=$T(5)) == $T[2, 5/2, 3, 7/2, 4, 9/2, 5]
                @eval @test range(start=$T(2), step=Int8(2), stop=$T(9/2)) == $T[2, 4]
                @eval @test range(start=$T(2), step=big(2), stop=$T(9/2)) == BigHalfInt[2, 4]
                @eval @test range(start=Int8(2), step=$T(1/2), stop=Int8(5)) == $T[2, 5/2, 3, 7/2, 4, 9/2, 5]
            end

            if T === :BigHalfInt
                @eval @test @inferred(range($T(1/2), stop=$T(9/2), length=5)) isa LinRange{BigFloat}
                @eval @test @inferred(range($T(2), stop=Int8(10), length=5)) isa LinRange{BigFloat}
                @eval @test range($T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                @eval @test range($T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                @static if VERSION ≥ v"1.7.0-DEV.263"
                    @eval @test @inferred(range(start=$T(1/2), stop=$T(9/2), length=5)) isa LinRange{BigFloat}
                    @eval @test @inferred(range(start=$T(2), stop=Int8(10), length=5)) isa LinRange{BigFloat}
                    @eval @test range(start=$T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                    @eval @test range(start=$T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                end
            elseif T === :HalfUInt64
                @eval @test_broken @inferred(range($T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                @eval @test_broken range($T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                @static if VERSION ≥ v"1.7.0-DEV.263"
                    @eval @test_broken @inferred(range(start=$T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test_broken range(start=$T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                end
                if Int === Int32
                    @eval @test_broken @inferred(range($T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test_broken range($T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                    @static if VERSION ≥ v"1.7.0-DEV.263"
                        @eval @test_broken @inferred(range(start=$T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                        @eval @test_broken range(start=$T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                    end
                else
                    @eval @test @inferred(range($T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test range($T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                    @static if VERSION ≥ v"1.7.0-DEV.263"
                        @eval @test @inferred(range(start=$T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                        @eval @test range(start=$T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                    end
                end
            elseif T === :HalfUInt128
                @eval @test_broken @inferred(range($T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                @eval @test_broken @inferred(range($T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                @eval @test_broken range($T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                @eval @test_broken range($T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                @static if VERSION ≥ v"1.7.0-DEV.263"
                    @eval @test_broken @inferred(range(start=$T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test_broken @inferred(range(start=$T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test_broken range(start=$T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                    @eval @test_broken range(start=$T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                end
            else
                if Int === Int32 && T ∈ (:HalfUInt32, :HalfUInt64, :HalfUInt128)
                    @eval @test_broken @inferred(range($T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test_broken range($T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                    @static if VERSION ≥ v"1.7.0-DEV.263"
                        @eval @test_broken @inferred(range(start=$T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                        @eval @test_broken range(start=$T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                    end
                else
                    @eval @test @inferred(range($T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test range($T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                    @static if VERSION ≥ v"1.7.0-DEV.263"
                        @eval @test @inferred(range(start=$T(1/2), stop=$T(9/2), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                        @eval @test range(start=$T(1/2), stop=$T(9/2), length=5) == [0.5, 1.5, 2.5, 3.5, 4.5]
                    end
                end
                @eval @test @inferred(range($T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                @eval @test range($T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                @static if VERSION ≥ v"1.7.0-DEV.263"
                    @eval @test @inferred(range(start=$T(2), stop=Int8(10), length=5)) isa StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
                    @eval @test range(start=$T(2), stop=Int8(10), length=5) == [2, 4, 6, 8, 10]
                end
            end
            @eval @test @inferred(range(big(10), stop=$T(2), length=5)) isa LinRange{BigFloat}
            @eval @test range(big(10), stop=$T(2), length=5) == [10, 8, 6, 4, 2]
            @static if VERSION ≥ v"1.7.0-DEV.263"
                @eval @test @inferred(range(start=big(10), stop=$T(2), length=5)) isa LinRange{BigFloat}
                @eval @test range(start=big(10), stop=$T(2), length=5) == [10, 8, 6, 4, 2]
            end

            @eval @test @inferred(range($T(2), step=$T(1/2), length=5)) isa StepRange{$T,$T}
            @eval @test @inferred(range($T(1/2), step=Int8(2), length=5)) isa StepRange{$T,Int8}
            @eval @test range($T(2), step=$T(1/2), length=5) == $T[2, 5/2, 3, 7/2, 4]
            @eval @test range($T(1/2), step=Int8(2), length=5) == $T[1/2, 5/2, 9/2, 13/2, 17/2]
            @eval @test range($T(1/2), step=big(2), length=5) == BigHalfInt[1/2, 5/2, 9/2, 13/2, 17/2]
            @static if VERSION ≥ v"1.7.0-DEV.16"
                @eval @test @inferred(range($T(1/2), step=big(2), length=5)) isa StepRange{BigHalfInt,BigInt}
                @eval @test @inferred(range(Int8(2), step=$T(1/2), length=5)) isa StepRange{$T,$T}
                @eval @test range(Int8(2), step=$T(1/2), length=5) == $T[2, 5/2, 3, 7/2, 4]
            else
                if T === :BigHalfInt
                    @eval @test @inferred(range($T(1/2), step=big(2), length=5)) isa StepRange{BigHalfInt,BigInt}
                else
                    @eval @test_broken @inferred(range($T(1/2), step=big(2), length=5)) isa StepRange{BigHalfInt,BigInt}
                end
                @eval @test_broken @inferred(range(Int8(2), step=$T(1/2), length=5)) isa StepRange{$T,$T}
                @eval @test_broken range(Int8(2), step=$T(1/2), length=5) == $T[2, 5/2, 3, 7/2, 4]
            end
            @static if VERSION ≥ v"1.7.0-DEV.263"
                @eval @test @inferred(range(start=$T(2), step=$T(1/2), length=5)) isa StepRange{$T,$T}
                @eval @test @inferred(range(start=$T(1/2), step=Int8(2), length=5)) isa StepRange{$T,Int8}
                @eval @test @inferred(range(start=$T(1/2), step=big(2), length=5)) isa StepRange{BigHalfInt,BigInt}
                @eval @test @inferred(range(start=Int8(2), step=$T(1/2), length=5)) isa StepRange{$T,$T}
                @eval @test range(start=$T(2), step=$T(1/2), length=5) == $T[2, 5/2, 3, 7/2, 4]
                @eval @test range(start=$T(1/2), step=Int8(2), length=5) == $T[1/2, 5/2, 9/2, 13/2, 17/2]
                @eval @test range(start=$T(1/2), step=big(2), length=5) == BigHalfInt[1/2, 5/2, 9/2, 13/2, 17/2]
                @eval @test range(start=Int8(2), step=$T(1/2), length=5) == $T[2, 5/2, 3, 7/2, 4]
            end
        end

        @static if VERSION ≥ v"1.7.0-DEV.263"
            for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
                if T in (:HalfUInt64, :HalfUInt128) || (Int === Int32 && T === :HalfUInt32)
                    @eval @test_broken @inferred(range(stop=$T(23/2), step=Int8(2), length=5)) isa StepRange{$T,Int8}
                    @eval @test_broken range(stop=$T(23/2), step=Int8(2), length=5) == $T[7/2, 11/2, 15/2, 19/2, 23/2]
                else
                    @eval @test @inferred(range(stop=$T(23/2), step=Int8(2), length=5)) isa StepRange{$T,Int8}
                    @eval @test range(stop=$T(23/2), step=Int8(2), length=5) == $T[7/2, 11/2, 15/2, 19/2, 23/2]
                end
                @eval @test @inferred(range(stop=$T(23/2), step=big(2), length=5)) isa StepRange{BigHalfInt,BigInt}
                @eval @test range(stop=$T(23/2), step=big(2), length=5) == BigHalfInt[7/2, 11/2, 15/2, 19/2, 23/2]
            end
            for T in (halfinttypes..., :BigHalfInt)
                @eval @test @inferred(range(stop=$T(2), step=$T(1/2), length=6)) isa StepRange{$T,$T}
                @eval @test @inferred(range(stop=Int8(2), step=$T(1/2), length=6)) isa StepRange{$T,$T}
                @eval @test range(stop=$T(2), step=$T(1/2), length=6) == $T[-1/2, 0, 1/2, 1, 3/2, 2]
                @eval @test range(stop=Int8(2), step=$T(1/2), length=6) == $T[-1/2, 0, 1/2, 1, 3/2, 2]
            end
        end
    end

    @testset "reverse" begin
        for T in (halfinttypes..., :BigHalfInt)
            # There is currently a discussion on how to treat Unsigned types:
            # * https://github.com/JuliaLang/julia/issues/29576
            # * https://github.com/JuliaLang/julia/pull/29842
            @eval @test @inferred(reverse($T(1/2):$T(9/2))) isa StepRange{$T}
            @eval @test @inferred(reverse($T(9/2):$T(-1):$T(1/2))) isa StepRange{$T}
            @eval @test reverse($T(1/2):$T(9/2)) == $T(9/2):$T(-1):$T(1/2)
            @eval @test reverse($T(9/2):$T(-1):$T(1/2)) == $T(1/2):$T(1):$T(9/2)
            @eval @test reverse($T(1/2):$T(1/2):$T(3)) == $T(3):$T(-1/2):$T(1/2)
            @eval @test reverse($T(-1/2):$T(2):$T(3/2)) == $T(3/2):$T(-2):$T(-1/2)
            @eval @test isempty(reverse($T(3/2):$T(1/2)))
            @eval @test isempty(reverse($T(3/2):$T(1):$T(1/2)))
            @eval @test isempty(reverse($T(3/2):$T(1/2):$T(1)))
            @eval @test isempty(reverse($T(7/2):$T(-2):$T(11/2)))
        end
    end

    @testset "intersect" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test @inferred(($T(1):$T(3)) ∩ ($T(1/2):$T(3))) isa UnitRange{$T}
            @eval @test @inferred((big(1):big(3)) ∩ ($T(1/2):$T(3))) isa UnitRange{BigHalfInt}
            @eval @test @inferred(($T(1/2):$T(3)) ∩ (big(1):big(3))) isa UnitRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(1):$T(5)) ∩ ($T(1):$T(3))) isa StepRange{$T}
            @eval @test @inferred(($T(3/2):$T(1):$T(5)) ∩ (big(1):big(3))) isa StepRange{BigHalfInt}
            @eval @test @inferred((big(1):big(3)) ∩ ($T(3/2):$T(1):$T(5))) isa StepRange{BigHalfInt}
            @eval @test @inferred((big(1):big(3):big(7)) ∩ ($T(3/2):$T(5))) isa StepRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(5)) ∩ (big(1):big(3):big(7))) isa StepRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(1/2):$T(5)) ∩ ($T(1):$T(3):$T(5))) isa StepRange{$T}
            @eval @test @inferred((big(-1):big(3):big(5)) ∩ ($T(3/2):$T(1/2):$T(5))) isa StepRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(1/2):$T(5)) ∩ (big(-1):big(3):big(5))) isa StepRange{BigHalfInt}
            @eval @test isempty(($T(1):$T(3)) ∩ ($T(1/2):$T(3)))
            @eval @test ($T(1/2):$T(2)) ∩ ($T(3/2):$T(3)) == $T(3/2):$T(3/2)
            @eval @test isempty(($T(1/2):$T(2)) ∩ (1:2))
            @eval @test isempty((1:2) ∩ ($T(1/2):$T(2)))
            @eval @test ($T(2):$T(5)) ∩ (1:3) == 2:3
            @eval @test (1:3) ∩ ($T(2):$T(5)) == 2:3
            @eval @test ($T(2):$T(1):$T(5)) ∩ (1:3) == 2:1:3
            @eval @test (1:3) ∩ ($T(2):$T(1):$T(5)) == 2:1:3
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ (1:3))
            @eval @test isempty((1:3) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test ($T(2):$T(1/2):$T(5)) ∩ (1:3) == 2:1:3
            @eval @test (1:3) ∩ ($T(2):$T(1/2):$T(5)) == 2:1:3
            @eval @test ($T(3/2):$T(1/2):$T(5)) ∩ (1:3) == 2:1:3
            @eval @test (1:3) ∩ ($T(3/2):$T(1/2):$T(5)) == 2:1:3
            @eval @test ($T(2):$T(1):$T(5)) ∩ ($T(1):$T(3)) == $T(2):$T(1):$T(3)
            @eval @test ($T(1):$T(3)) ∩ ($T(2):$T(1):$T(5)) == $T(2):$T(1):$T(3)
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ ($T(1):$T(3)))
            @eval @test isempty(($T(1):$T(3)) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test ($T(2):$T(1/2):$T(5)) ∩ ($T(1):$T(3)) == $T(2):$T(1):$T(3)
            @eval @test ($T(1):$T(3)) ∩ ($T(2):$T(1/2):$T(5)) == $T(2):$T(1):$T(3)
            @eval @test ($T(3/2):$T(1/2):$T(5)) ∩ ($T(1):$T(3)) == $T(2):$T(1):$T(3)
            @eval @test ($T(1):$T(3)) ∩ ($T(3/2):$T(1/2):$T(5)) == $T(2):$T(1):$T(3)
            @eval @test ($T(3/2):$T(1):$T(6)) ∩ ($T(1):$T(3/2):$T(6)) == $T(5/2):$T(3):$T(11/2)
        end
        for T in (halfinttypes..., :BigHalfInt)
            @eval @test ($T(3/2):$T(1/2):$T(5)) ∩ (-1:3:5) == 2:3:5
            @eval @test (-1:3:5) ∩ ($T(3/2):$T(1/2):$T(5)) == 2:3:5
            @eval @test ($T(2):$T(1):$T(5)) ∩ (-1:3:5) == 2:3:5
            @eval @test (-1:3:5) ∩ ($T(2):$T(1):$T(5)) == 2:3:5
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ (-1:3:5))
            @eval @test isempty((-1:3:5) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ ($T(-1):$T(3):$T(5)))
            @eval @test isempty(($T(-1):$T(3):$T(5)) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test ($T(3/2):$T(1):$T(5)) ∩ ($T(-1):$T(3/2):$T(5)) == $T(7/2):$T(3):$T(7/2)
            # These tests should in theory also work for Unsigned types, but due to some
            # bugs with ranges of Unsigned integers, they currently don’t (see below)
            @eval @test ($T(5/2):$T(3/2):$T(7)) ∩ ($T(1):$T(8)) == $T(4):$T(3):$T(7)
            @eval @test ($T(1):$T(8)) ∩ ($T(5/2):$T(3/2):$T(7)) == $T(4):$T(3):$T(7)
            @eval @test isempty((1:3:7) ∩ ($T(3/2):$T(5)))
            @eval @test isempty(($T(3/2):$T(5)) ∩ (1:3:7))
            @eval @test (1:3:7) ∩ ($T(3):$T(5)) == 4:3:4
            @eval @test ($T(3):$T(5)) ∩ (1:3:7) == 4:3:4
        end
        # Ranges of UInts are currently broken, cf. e.g.
        # * https://github.com/JuliaLang/julia/issues/29576
        # * https://github.com/JuliaLang/julia/issues/29801
        # * https://github.com/JuliaLang/julia/issues/29810
        # These (or other) bugs related to ranges of Unsigned integers prevent some of the
        # following tests from passing
        for T in halfuinttypes
            @eval @test_broken ($T(5/2):$T(3/2):$T(7)) ∩ ($T(1):$T(8)) == $T(4):$T(3):$T(7)
            @eval @test_broken ($T(1):$T(8)) ∩ ($T(5/2):$T(3/2):$T(7)) == $T(4):$T(3):$T(7)
        end
        for T in (:HalfUInt8, :HalfUInt16)
            @eval @test isempty((1:3:7) ∩ ($T(3/2):$T(5)))
            @eval @test isempty(($T(3/2):$T(5)) ∩ (1:3:7))
            @eval @test (1:3:7) ∩ ($T(3):$T(5)) == 4:3:4
            @eval @test ($T(3):$T(5)) ∩ (1:3:7) == 4:3:4
        end
        @static if Sys.WORD_SIZE == 64
            @test isempty((1:3:7) ∩ (HalfUInt32(3/2):HalfUInt32(5)))
            @test isempty((HalfUInt32(3/2):HalfUInt32(5)) ∩ (1:3:7))
            @test (1:3:7) ∩ (HalfUInt32(3):HalfUInt32(5)) == 4:3:4
            @test (HalfUInt32(3):HalfUInt32(5)) ∩ (1:3:7) == 4:3:4
        else
            @test_broken isempty((1:3:7) ∩ (HalfUInt32(3/2):HalfUInt32(5)))
            @test_broken isempty((HalfUInt32(3/2):HalfUInt32(5)) ∩ (1:3:7))
            @test_broken (1:3:7) ∩ (HalfUInt32(3):HalfUInt32(5)) == 4:3:4
            @test_broken (HalfUInt32(3):HalfUInt32(5)) ∩ (1:3:7) == 4:3:4
        end
        for T in (:HalfUInt64, :HalfUInt128)
            @eval @test_broken isempty((1:3:7) ∩ ($T(3/2):$T(5)))
            @eval @test_broken isempty(($T(3/2):$T(5)) ∩ (1:3:7))
            @eval @test_broken (1:3:7) ∩ ($T(3):$T(5)) == 4:3:4
            @eval @test_broken ($T(3):$T(5)) ∩ (1:3:7) == 4:3:4
        end
    end
end


text:
module TestExamplesPrimes
include("../examples/primes.jl")
end  # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -24.652335523355834 <= q <= 21.32472335033785)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.6982139593401175 + -0.19504968807074974*1/(1 + exp(-(0.6793369023557885 + 0.6523875537646115*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.5142038874540145*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.0456439541419269*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6096854328317232*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + -0.41803298572193937*1/(1 + exp(-(-0.7048442383777926 + 0.47781322839334983*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.713369260368808*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.44964139924772173*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.9876898123617992*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + 0.9542284451941447*1/(1 + exp(-(0.9969387928417852 + -0.011395717836928387*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + 0.8103499615759371*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.12667586987008894*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6252342222679741*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + 0.7916680056214345*1/(1 + exp(-(-0.5156629799588863 + -0.22745979325095078*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.8222624364306839*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.4937419854366709*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.1046356011041838*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))))))))) + 1/(1 + exp(-(-0.4539405071565561 + 0.5320780591355936*1/(1 + exp(-(0.6793369023557885 + 0.6523875537646115*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.5142038874540145*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.0456439541419269*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6096854328317232*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + -0.8489788623944037*1/(1 + exp(-(-0.7048442383777926 + 0.47781322839334983*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.713369260368808*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.44964139924772173*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.9876898123617992*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + -0.25620536200506905*1/(1 + exp(-(0.9969387928417852 + -0.011395717836928387*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + 0.8103499615759371*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.12667586987008894*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6252342222679741*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + 0.9264345634068367*1/(1 + exp(-(-0.5156629799588863 + -0.22745979325095078*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.8222624364306839*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.4937419854366709*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.1046356011041838*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))))))))) + 1/(1 + exp(-(-0.4017849565754932 + 0.7092320156791541*1/(1 + exp(-(0.6793369023557885 + 0.6523875537646115*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.5142038874540145*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.0456439541419269*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6096854328317232*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + -0.20657569919094643*1/(1 + exp(-(-0.7048442383777926 + 0.47781322839334983*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.713369260368808*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.44964139924772173*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.9876898123617992*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + 0.6260976974386958*1/(1 + exp(-(0.9969387928417852 + -0.011395717836928387*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + 0.8103499615759371*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.12667586987008894*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6252342222679741*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + 0.026773003462286393*1/(1 + exp(-(-0.5156629799588863 + -0.22745979325095078*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.8222624364306839*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.4937419854366709*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.1046356011041838*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))))))))) + 1/(1 + exp(-(0.8541923459688499 + -0.7979821016654509*1/(1 + exp(-(0.6793369023557885 + 0.6523875537646115*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.5142038874540145*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.0456439541419269*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6096854328317232*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + 0.2172669330681649*1/(1 + exp(-(-0.7048442383777926 + 0.47781322839334983*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.713369260368808*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.44964139924772173*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.9876898123617992*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + -0.6626207135044284*1/(1 + exp(-(0.9969387928417852 + -0.011395717836928387*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + 0.8103499615759371*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.12667586987008894*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + 0.6252342222679741*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + -0.5735396155482761*1/(1 + exp(-(-0.5156629799588863 + -0.22745979325095078*1/(1 + exp(-(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3])))) + -0.8222624364306839*1/(1 + exp(-(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3])))) + -0.4937419854366709*1/(1 + exp(-(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3])))) + -0.1046356011041838*1/(1 + exp(-(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "Test function `ENG.pluralize`" begin
    PostgresORMUtil.Pluralize.ENG.pluralize("analysis")
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Abaqus
using JuliaFEM.Testing

@testset "parse abaqus inp file to AbaqusModel" begin
    fn = Pkg.dir("JuliaFEM") * "/test/testdata/cube_tet4.inp"
    model = abaqus_read_model(fn)

    @test length(model.properties) == 1
    section = first(model.properties)
    @test section.element_set == :CUBE
    @test section.material_name == :MAT

    @test haskey(model.materials, :MAT)
    material = model.materials[:MAT]
    @test isapprox(first(material.properties).E, 208.0e3)

    @test length(model.steps) == 1
    step = first(model.steps)
    @test length(step.boundary_conditions) == 2

    bc = step.boundary_conditions[1]
    @test bc.data[1] == [:SYM12, 3]
    @test bc.data[2] == [:SYM23, 1]
    @test bc.data[3] == [:SYM13, 2]

    load = step.boundary_conditions[2]
    @test load.data[1] == [:LOAD, :P, 1.00000]
end

#=
@testset "given abaqus model solve field" begin
    fn = Pkg.dir("JuliaFEM") * "/test/testdata/cube_tet4.inp"
    model = abaqus_read_model(fn)
    problems = model()
    body = first(problems)
    info(body("displacement", 0.0))
    result = XDMF()
    xdmf_new_result!(result, body, 0.0)
    xdmf_save_field!(result, body, 0.0, "displacement"; field_type="Vector")
    xdmf_save!(result, "/tmp/cube_tet4.xmf")
end
=#


text:
InputFolder = './Images/CellsSpots/';
OutputFolder = './Results/Images/CellsSpots/';
Lbl = 1;
@iA = '*.tif';

@fxg_sLoGLocMax [iA] > [M];
params.GRad = 2;
params.LocalMaxBox = [5 5];
params.ThrLocalMax = 2;
params.MinLocThr = 0;
/endf

@fxs_lMeanShift3D [M] > [S L];
params.Bandwidth = 40;
params.MinPts = 3;
/endf

/show iA > L;
/keep L > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
"""
    chroma(signal, sr, window=ones, windowsize=2048, overlap=512, γ=0.1, ϵ=0.1)

Compute the chroma features of a signal.

The return will be a matrix of size 12xN where N is the number of frames of the
STFT and each column is a vector of normalized chroma features.

**Parameters:**

`signal`: The signal on which the coefficients will be computed.

`sr`: The sampling rate of the signal.

`window`: The window function to be applied on the signal by the STFT.

`windowsize`: The size of the STFT window.

`overlap`: The overlap of two adjacent windows of the STFT in samples.

`γ`: The compression factor to be used in the chroma features extraction.

`ϵ`: The smallest norm that a chroma vector needs to have to not be replaced by
a stub one.
"""
function chroma(signal, sr, window=ones, windowsize=2048, overlap=512, γ=0.1, ϵ=0.1)

	sigstft = stft(signal, windowsize, overlap, window=window)

	# If the signal is real then the absolute value of it's spectrum is symmetric
	eltype(signal)<:Real && (windowsize = windowsize ÷ 2 + 1)

	# The number of chroma vectors is defined by the frames in the STFT
	nofframes = size(sigstft)[2]
	pitchspectrum = Matrix{Real}(undef, 128, nofframes)

	# Relation of a pitch with a frequency defined by the MIDI standard
	pitch2freq(p) = round(Cint, (windowsize/sr)*2^((p - 69)/12)*440)

	# The bound frequencies of each one of the pitch classes
	bounds = [1 + pitch2freq(p - 0.5) for p in 0:128]

	for n in 1:nofframes
		# The pitch spectrum is a log spaced spectrum where a component is
		# proportional to the sum of the power of all stft components inside the
		# bounds of a pitch class
		pitchspectrum[:, n] =
			[sum(abs2.(sigstft[bounds[i]:bounds[i + 1],n])) for i in 1:128]
	end

	chroma = Matrix{Real}(undef, 12, nofframes)

	for n in 1:nofframes
		# A chroma vector is given by the sum of the pitch classes of the same
		# note in all octaves
		chroma[:, n] =
			[sum(pitchspectrum[(class % 12 + 1):12:end, n]) for class in 0:11]
	end

	# Compresses the chroma vectors
	chroma = log10.(1 .+ γ .* chroma)

	stubFeature = ones(12)./√12

	# Evaluates the norm of each chroma vector...
	for i in 1:nofframes
		# ...and if it's smaller than the ϵ parameter normalizes the vector
		if ((norm = √(sum(chroma[:, i].^2))) > ϵ)
			chroma[:, i] ./= norm
		# Otherwise replaces the vector with a stub one
		else
			chroma[:, i] = stubFeature
		end
	end

	return chroma

end

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
using KiteUtils
using Test

cd("..")

@testset "KiteUtils.jl: Settings       " begin
    @test se().project == "settings.yaml"
    @test se().log_file == joinpath("data", "log_8700W_8ms")
    @test se().time_lapse == 1.0
    @test se().sim_time == 100.0
    @test length(se().alpha_cl) == 12
    set_data_path(tempdir())
    @test KiteUtils.DATA_PATH[1] == tempdir()
    set_data_path("data")
    set2 = load_settings(joinpath("data", "settings.yaml"))
    @test set2.project == "settings.yaml"
end
@testset "KiteUtils.jl: Log files      " begin
    state = KiteUtils.demo_state(7)
    @test typeof(state) == SysState{7}
    @test state.X[end] == 10.0
    @test all(state.pos[end] .≈ [10, 0, 6.0])
    @test repr(state) == "time      [s]:       0.0\norient    [w,x,y,z]: Float32[0.5, 0.5, -0.5, -0.5]\nelevation [rad]:     0.5404195\nazimuth   [rad]:     0.0\nl_tether  [m]:       0.0\nv_reelout [m/s]:     0.0\nforce     [N]:       0.0\ndepower   [-]:       0.0\nv_app     [m/s]:     0.0\nX         [m]:       Float32[0.0, 1.6666666, 3.3333333, 5.0, 6.6666665, 8.333333, 10.0]\nY         [m]:       Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\nZ         [m]:       Float32[0.0, 0.15380114, 0.6194867, 1.4100224, 2.5474184, 4.063342, 6.0000005]\n"
    state = KiteUtils.demo_state_4p(7)
    @test typeof(state) == SysState{11}
    @test state.X[end] ≈ 12.45
    @test state.Y[end] ≈ -1.115
    @test state.Y[end-1] ≈ 1.115
    set_data_path(tempdir())
    log = KiteUtils.test(true)
    @test typeof(log) == SysLog{7}
    @test log.syslog.Z[end][7] ≈ 6 # height of the last particle which represents the kite
    @test log.z[end] ≈ 6.0
    @test log.y[end] ≈ 0.0
    @test log.x[end] ≈ 10.0
    @test export_log(log) == joinpath(tempdir(), "Test_flight.csv")
end
@testset "KiteUtils.jl: Transformations" begin
    ax, ay, az = [1, 0, 0], [0, 1, 0],  [0, 0, 1]
    bx, by, bz = [0, 1, 0], [-1, 0, 0], [0, 0, 1]
    res = rot3d(ax, ay, az, bx, by, bz)
    vec = [1, 0, 0]
    @test res * vec ≈ [0, 1, 0]
    pos_kite = [1.0, 1, 10]
    pos_before =  [1.0, 1, 9]
    v_app = [10, 0, 0.0]
    m = rot(pos_kite, pos_before, v_app)
    @test m == [0.0 0.0 -1.0; -1 0.0 0.0; 0 1 0]
    @test ground_dist(pos_kite) ≈ 1.4142135623730951
    @test calc_elevation(pos_kite) ≈ 1.4303066250413763
    @test azimuth_east(pos_kite) ≈ -0.7853981633974483
    @test acos2(1.0001) == 0.0
    @test acos2(-1.0001) ≈ π
end
nothing


text:
# using TuringModels
using CSV, DataFrames, Query, Turing, Optim, StatsBase

# define precis from StatisticalRethinking module
using NamedArrays
function precis(df::DataFrame; digits=3, depth=Inf, alpha=0.11)
  m = zeros(length(names(df)), 5)
  for (indx, col) in enumerate(names(df))
    m[indx, 1] = mean(df[:, col])
    m[indx, 2] = std(df[:, col])
    q = quantile(df[:, col], [alpha/2, 0.5, 1-alpha/2])
    m[indx, 3] = q[1]
    m[indx, 4] = q[2]
    m[indx, 5] = q[3]
  end
  NamedArray(
    round.(m, digits=digits),
    (names(df), ["mean", "sd", "5.5%", "50%", "94.5%"]),
    ("Rows", "Cols"))
end

# ### snippet 5.3

df = CSV.read(joinpath(@__DIR__, "data", "Howell1.csv"))

# create categorical sex variable
df = df |> @mutate(sexcat = _.male == 1 ? "male" : "female") |> DataFrame
categorical!(df, :sexcat)

```
get levels code to create index variable:
this could have been done like this:
    df = df |> @mutate(sexcat = _.male == 1 ? 2 : 1) |> DataFrame
for the same result, but this code will do the job if you already
have a categorial variable.
```
df = df |> @mutate(sex = levelcode(_.sexcat)) |> DataFrame

@model model5_8(y, group) = begin
    n_groups = length(unique(group))
    σ ~ Uniform(0, 50)
    α ~ filldist(Normal(178, 20), n_groups)
    μ = α[group]
    @. y ~ Normal(μ, σ)
end

y = df.height
group = df.sex

model = model5_8(y, group)
m5_8 = optimize(model, MAP())

m5_8 = sample(model, NUTS(), MCMCThreads(), 2000, 4)
describe(m5_8)

post = DataFrame(m5_8)

# post |> @rename("α[1]" => "male")

# retrieve category names
dict = zip(["α[1]", "α[2]"], levels(df.sexcat)) |> Dict
rename!(post, dict)

post = post |> @mutate(diff = _.female - _.male) |> DataFrame

precis(post)


text:
#find how many different letters are there
#counting same letters just once
function find_yes(group_ans)
   #remove spaces
   group_ans_no_spaces = replace(group_ans, " " => "")
   #get size of array with unique occurrences in the string
   number_of_y = length(unique(group_ans_no_spaces))
   return number_of_y
end

f = open("input_day6","r")

temp_s = ""
sum_y = 0

#note that I modified my inputs to have an empty line at the end
for lines in readlines(f)
   #if not empty line concatenate de string with space
   if lines != ""
      global temp_s = temp_s*" "*lines
   else
      #println(temp_s)
      n = find_yes(temp_s)
      global sum_y += n
      global temp_s = ""

   end
end

println(sum_y)


text:
# implementation of 5 stage, 4th order Low Storage Explicit Runge Kutta
# from Carpenter and Kennedy's Fourth-Order 2N-Storage Runge Kutta Schemes
# NASA Technical Memorandum 109112

# this code borrows some infrastructure (pde_pre_func, pde_post_func) from
# the classical rk4 file
# Also borrows the RK$CheckpointData

"""
  This function implements the 5 stage, 4th order Low Storage Explicit Runge Kutta scheme of
  Carpenter and Kennedy

  Arguments:
    f: a function that evalutes du/dt = f(q, t)
    delta_t: the time step
    t_max: the maximum time value
    q_vec: vector (of length numDof) containing initial solution.  Will contain final solution
           at exit
    res_vec: vector to store the residual in during evaluation of f.  The contents of this vector
             at exit is undefined
    pre_func: function to call after new values are written into q_vec but before f is called
    post_func: function to call after f is called but before res_vec is accessed
    ctx: tuple arguments of f (ie. f = f(ctx...))
    opts: options dictionary
    timing: a Timings object

  Keyword Arguments:
    majorIterationCallback: function to call after first function evaluation of each time step, ie.
                            when q_vec and res_vec have been updated.  Useful for logging.  Defaults
                            to no-op
    res_tol: stopping tolerance for residual (useful for pseudo-timestepping), default -1.0
    real_time: whether or not to advance time (ie. pseudo timestepping or not) default faulse

  See the documentation for rk4.
"""

function lserk54(f::Function, delta_t::AbstractFloat, t_max::AbstractFloat,
             q_vec::AbstractVector, res_vec::AbstractVector, pre_func,
             post_func, ctx, opts, timing::Timings=Timings();
             majorIterationCallback=((a...) -> (a...)),
             res_tol = -1.0, real_time=false)

  # LSERK coefficients
  const a_coeffs = [0;
                    -567301805773.0/1357537059087.0;
                    -2404267990393.0/2016746695238.0;
                    -3550918686646.0/2091501179385.0;
                    -1275806237668.0/842570457699.0]

  const b_coeffs = [1432997174477.0/9575080441755.0;
                    5161836677717.0/13612068292357.0;
                    1720146321549.0/2090206949498.0;
                    3134564353537.0/4481467310338.0;
                    2277821191437.0/14882151754819.0]

  const c_coeffs = [0;
                    1432997174477/9575080441755;
                    2526269341429/6820363962896;
                    2006345519317/3224310063776;
                    2802321613138/2924317926251]

  myrank = MPI.Comm_rank(MPI.COMM_WORLD)
#  MPI.Barrier(MPI.COMM_WORLD)
  if myrank == 0
    println(BSTDOUT, "\nEntered lserk54")
    println(BSTDOUT, "res_tol = ", res_tol)
  end
#  flush(BSTDOUT)
#  MPI.Barrier(MPI.COMM_WORLD)
# res_tol is alternative stopping criteria

  # unpack options
  output_freq = opts["output_freq"]::Int
  use_itermax = opts["use_itermax"]::Bool
  if use_itermax
    itermax = opts["itermax"]
  end

  use_checkpointing = opts["use_checkpointing"]::Bool
  chkpoint_freq = opts["checkpoint_freq"]::Int
  ncheckpoints = opts["ncheckpoints"]::Int

  t = 0.0  # timestepper time
  treal = 0.0  # real time (as opposed to pseudo-time)
  t_steps = round(Int, t_max/delta_t)
  println(BSTDOUT, "t_steps: ",t_steps)
  println(BSTDOUT, "delta_t = ", delta_t)

  (m,) = size(q_vec)

  # allocate storage
  # this is actually a 3N scheme because the function f cannot overwrite its
  # input
  dq_vec = zeros(q_vec)

  if myrank == 0
    _f1 = open("convergence.dat", "a")
    f1 = BufferedIO(_f1)
  end

  # setup all the checkpointing related data
  chkpointer, chkpointdata, skip_checkpoint = explicit_checkpoint_setup(opts, myrank)
  istart = chkpointdata.i

  flush(BSTDOUT)
  #-----------------------------------------------------
  # Main timestepping loop
  timing.t_timemarch += @elapsed for i=istart:(t_steps + 1)

    t = (i-2)*delta_t

    @mpi_master if i % output_freq == 0
       println(BSTDOUT, "\ntimestep ",i)
       if i % output_freq == 0
         flush(BSTDOUT)
       end
    end

    if use_checkpointing && i % chkpoint_freq == 0
      if skip_checkpoint    # skip only the first checkpoint
        skip_checkpoint = false
      else
        @mpi_master println(BSTDOUT, "Saving checkpoint at timestep ", i)
        skip_checkpoint = false
        # save all needed variables to the chkpointdata
        chkpointdata.i = i

        if countFreeCheckpoints(chkpointer) == 0
          freeOldestCheckpoint(chkpointer)  # make room for a new checkpoint
        end

        # save the checkpoint
        saveNextFreeCheckpoint(chkpointer, ctx..., opts, chkpointdata)
      end   # end of if skip_checkpoint check
    end   # end of if use_checkpointing check

    #--------------------------------------------------------------------------
    # stage 1
#    f(params, u, F_vals, t_i)

    pre_func(ctx..., opts)
    if real_time treal = t end
    timing.t_func += @elapsed f(ctx..., opts, treal)
    sol_norm = post_func(ctx..., opts)

    #--------------------------------------------------------------------------
    # callback and logging
    timing.t_callback += @elapsed majorIterationCallback(i, ctx..., opts, BSTDOUT)

    # logging
    @mpi_master if i % 1 == 0
      println(f1, i, " ", sol_norm)
    end

    @mpi_master if i % output_freq == 0
      println(BSTDOUT, "flushing convergence.dat to disk")
      flush(f1)
    end

    # check stopping conditions
    if (sol_norm < res_tol) && !real_time
      if myrank == 0
        println(BSTDOUT, "breaking due to res_tol, res norm = $sol_norm")
        close(f1)
        flush(BSTDOUT)
      end
      break
    end

    if use_itermax && i > itermax
      if myrank == 0
        println(BSTDOUT, "breaking due to itermax")
        close(f1)
        flush(BSTDOUT)
      end
      break
    end

    #--------------------------------------------------------------------------
    # remaining stages

    # stage 1 update

    fac = b_coeffs[1]
    for j=1:length(q_vec)
      dq_vec[j] = delta_t*res_vec[j]
      q_vec[j] += fac*dq_vec[j]
    end

    # loop over remaining stages
    for stage=2:5
      pre_func(ctx..., opts)
      if real_time
        treal = t + c_coeffs[stage]*delta_t
      end
      timing.t_func += @elapsed f( ctx..., opts, treal)
      post_func(ctx..., opts, calc_norm=false)

      # update
      fac = a_coeffs[stage]
      fac2 = b_coeffs[stage]
      for j=1:length(q_vec)
        dq_vec[j] = fac*dq_vec[j] + delta_t*res_vec[j]
        q_vec[j] += fac2*dq_vec[j]
      end
    end  # end loop over stages

  end  # end loop over timesteps

  # final update
  t += delta_t

  @mpi_master println("---------------------------------------------")
  @mpi_master println("   LSERK: final time step reached. t = $t")
  @mpi_master println("---------------------------------------------")

  if myrank == 0
    close(f1)
  end

  flush(BSTDOUT)

  return t
end  # end lserk54

"""
  See rk4 method with same signature
"""
function lserk54(f::Function, h::AbstractFloat, t_max::AbstractFloat,
             q_vec::AbstractVector, res_vec::AbstractVector, ctx, opts,
             timing::Timings=Timings();
             majorIterationCallback=((a...) -> (a...)), res_tol=-1.0,
             real_time=false)

  t = lserk54(f::Function, h::AbstractFloat, t_max::AbstractFloat,
              q_vec::AbstractVector, res_vec::AbstractVector,
              pde_pre_func, pde_post_func, ctx, opts, timing;
              majorIterationCallback=majorIterationCallback, res_tol=res_tol,
              real_time=real_time)

        return t
end


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    optimize(params::AdmmParams[, settings::Settings[, vars::Vars]])

Carries out the ADMM algorithm by calling `admm_step` until termination conditions are met.
"""
function optimize(params::AdmmParams, settings::Settings, vars::Vars)
    assert_valid(params, settings, vars)
    kkt = factorize_kkt(params, settings)
    pc = ProxCache(params, settings)
    stats = Stats(settings)
    term_cache = get_term_cond_cache(params, settings)
    iter = 0
    while iter < settings.max_iters
        iter += 1
        admm_step!(vars, params, kkt, pc, settings, iter)
        if settings.compute_stats
            compute_stats(vars, params, stats, iter)
        end
        if check_term_conds!(term_cache, vars, params, settings, iter)
            break
        end
    end
    return vars, stats
end
optimize(params::AdmmParams, settings::Settings) = optimize(params, settings, Vars(params))
optimize(params::AdmmParams) = optimize(params, Settings(params))

@doc raw"""
    factorize_kkt(params::AdmmParams, settings::Settings)

Factorize the KKT matrix representing the optimzality equations of the provided problem.

# Extended help
Factorize the coefficient matrix of that KKT system representing the optimality equations
that will be used (and reused) for the constrained least squares problem:

```math
\begin{array}{ll}
\text{minimize} & \displaystyle
\frac12 x^T P x  + q^T x + \sum_{i=1}^n g_i(x_i) \\
\text{subject to} & Ax = b.
\end{array}
```

We can formulate this equivalently as

```math
\begin{array}{ll}
\text{minimize} & \displaystyle
\frac12 \tilde x^T P \tilde x  + q^T \tilde x + \mathcal{I}_\mathcal{A}(\tilde x, \tilde z)
        + g(x) + \mathcal{I}_{\mathcal B}(z) \\
\text{subject to} & \tilde x = x \\
                  & \tilde z = z.
\end{array}
```

where ``\mathcal I_{\mathcal A}`` is the indicator function over the affine
constraints of the original form, *i.e.*,
```math
\mathcal{I_{\mathcal A}}(x, z) = \begin{cases}
0 & Ax = z \\
\infty & {\rm otherwise}.
\end{cases}
```

and ``\mathcal I_{\mathcal B}`` is the indicator function over the singleton
set of vectors equal to ``b``, *i.e.*,

```math
\mathcal{I_{\mathcal B}}(x) = \begin{cases}
0 & x = b \\
\infty & {\rm otherwise}.
\end{cases}
```

We will denote the objective in this forumulation as ``f(x, \tilde x, z, \tilde z)``.

The augmented Lagrangian for this problem is

```math
L_{S,R} (x, \tilde x, z, \tilde z, w, y) =
f(x, \tilde x, z, \tilde z)
+ \frac12 \| \tilde x - x + S^{-1} w \|_S^2
+ \frac12 \| \tilde z - z + R^{-1} y \|_R^2,
```

where ``f`` is the previously mentioned objective function.
The diagonal, positive definite matrices ``S\in \bf S{++}^l`` and ``R\in \bf S_{++}^m``
are algorithm parameters.
We define ``\|\cdot \|_S`` is the ``S``-norm, *i.e.*,
``\|x\|_S = \sqrt{x^T S x}``.

Solving for the optimality conditions, we can find optimal ``\tilde x`` and ``\tilde z`` by solving the linear system

```math
\begin{bmatrix} P + S & A^T \\ A & -R^{-1} \end{bmatrix}
\begin{bmatrix} \tilde x \\ \nu \end{bmatrix}
=
\begin{bmatrix} Sx - q - w \\ z - R^{-1} y \end{bmatrix}.
```

where ``\tilde z`` can be recovered as ``\tilde z = z + R^{-1}(\nu - y)``

In this function, we factor the coefficient matrix on the left hand side so that we repeatedly
solve the linear system efficiently.
"""
function factorize_kkt(params::AdmmParams, settings::Settings)
    AT = sparse(params.A')
    Dσ = sparse(Diagonal(settings.σ))
    Dρ = sparse(Diagonal(1 ./ settings.ρ))

    # construct the block matrix K
    K = [params.P+Dσ AT;
         params.A -Dρ]
    return ldl(K)
end

"""
    solve_kkt!(vars::Vars, kkt::KKT, params::AdmmParams, settings::Settings)

Solve the KKT system set up in the description of `factorize_kkt` for a particular right hand side.
"""
function solve_kkt!(vars::Vars, kkt::KKT, params::AdmmParams, settings::Settings)
    m, n = size(params.A)
    rhs = Vector{Float64}(undef, n + m)
    for i in 1:n
        rhs[i] = settings.σ[i] * vars.x[i] - vars.w[i] - params.q[i]
    end
    for i in 1:m
        rhs[i + n] = vars.z[i] - (vars.y[i] / settings.ρ[i])
    end
    lhs = kkt \ rhs
    xt = lhs[1:n]
    ν = lhs[(n + 1):end]
    zt = vars.z + (ν - vars.y) ./ settings.ρ
    return xt, zt
end

"""
    admm_step!(vars::Vars, params::AdmmParams, kkt::KKT, pc::ProxCache, settings::Settings, iter::Int64; polish::Bool=false)

Carry out a single ADMM iteration.

# Steps:
1. Solve the KKT system to carry out the `xt` and `zt` updates.
2. Use a mixture of the result of (1) and the previous `x` iterate to solve `n` univariate
   subproblems by minimizing the Lagrangian in the description of `factorize_kkt` w.r.t.
   `x`. (Requires evaluation of a proximal operator.)
3. Update `w` and `y` using the residuals associated with the constraints `x = xt` and `z = zt`.

# Note:
* Note that for the entire run of the algorithm, the `z` update is just `z = b`, so we don't need
to do anything here.
"""
function admm_step!(vars::Vars, params::AdmmParams, kkt::KKT, pc::ProxCache,
                    settings::Settings, iter::Int64; polish::Bool=false)
    # update xt and zt by solving the KKT system
    xt, zt = solve_kkt!(vars, kkt, params, settings)

    # get the next x value by using the proximal operators of the separable function pieces
    prox_target = settings.α * xt + (1 - settings.α) * vars.x + vars.w ./ settings.σ
    x = prox_step(pc, params, settings, prox_target)
    @assert all(isfinite.(x))

    # update the dual variables
    z = params.b
    w = vars.w + settings.σ .* (settings.α * xt + (1 - settings.α) * vars.x - x)
    y = vars.y + settings.ρ .* (settings.α * zt + (1 - settings.α) * vars.z - z)

    # save in vars struct
    vars.xt = xt
    vars.zt = zt
    vars.x = x
    vars.w = w
    vars.z = z
    return vars.y = y
end

"""
    compute_stats(vars::Vars, params::AdmmParams, stats::Stats, iter::Int64)

Update the objective, residual, and iteration fields of the Stats object.
"""
function compute_stats(vars::Vars, params::AdmmParams, stats::Stats, iter::Int64)
    x = vars.x
    P = params.P
    q = params.q
    A = params.A
    b = params.b
    g = params.g

    fx = 0.5 * x' * P * x + q' * x
    gx = 0.0
    for i in 1:get_num_vars(params)
        if !isfinite(g[i, 1](x[i]))
        else
            gx += g[i, 1](x[i])
        end
    end
    stats.obj[iter] = fx + gx
    stats.iters = iter

    res = norm(A * x - b)
    return stats.res[iter] = res
end


text:
using Sockets

server = Sockets.listen(8080)
while true
    conn = accept(server) @async begin
        try
            while true
                line = readline(conn)
                println(line)
                write(conn,line)
            end
        catch ex
            println("connection ended with error $ex")
        end
    end # end coroutine block
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function solution()
  sum, curr, next = 0, 1, 2

  while curr < 4000000
    if curr % 2 == 0
      sum += curr
    end

    temp = next
    next += curr
    curr = temp
  end

  sum
end

if length(ARGS) == 1 && ARGS[1] == "-a"
  println(solution())
  exit(0)
end

for line in eachline(STDIN)
  iters = int(line)

  start = time_ns()
  for i in 1:iters
    solution()
  end
  end_ = time_ns()

  println(end_ - start)
end


text:
"""
    LBMCollisionCallback()

Apply the Lattice-Boltzmann method (LBM) collision operator before each time step.
See [`LatticeBoltzmannEquations2D`](@ref) for further details.
"""
function LBMCollisionCallback()
  DiscreteCallback(lbm_collision_callback, lbm_collision_callback,
                   save_positions=(false,false),
                   initialize=initialize!)
end

# Always execute collision step after a time step, but not after the last step
lbm_collision_callback(u, t, integrator) = !isfinished(integrator)

function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:typeof(lbm_collision_callback)})
  @nospecialize cb # reduce precompilation time

  print(io, "LBMCollisionCallback()")
end

function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:typeof(lbm_collision_callback)})
  @nospecialize cb # reduce precompilation time

  if get(io, :compact, false)
    show(io, cb)
  else
    summary_box(io, "LBMCollisionCallback")
  end
end

# Execute collision step once in the very beginning
function initialize!(cb::DiscreteCallback{Condition,Affect!}, u, t, integrator) where {Condition, Affect!<:typeof(lbm_collision_callback)}
  cb.affect!(integrator)
end

# This method is called as callback after the StepsizeCallback during the time integration.
@inline function lbm_collision_callback(integrator)

  dt = get_proposed_dt(integrator)
  semi = integrator.p
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  @unpack collision_op = equations

  u_ode = integrator.u
  u = wrap_array(u_ode, mesh, equations, solver, cache)

  @timed timer() "LBM collision" apply_collision!(u, dt, collision_op, mesh, equations, solver, cache)

  return nothing
end

include("lbm_collision_dg2d.jl")
include("lbm_collision_dg3d.jl")


text:
using Base.Test
using TimeIt
import NestedSampling: logaddexp, nball_vol_factor, sample

# Test logaddexp()
x = log(1.e-100)
y = log(2.e-100)
@test_approx_eq exp(logaddexp(x, y)) 3.e-100

# Test volume prefactor for an N-ball
@test_approx_eq nball_vol_factor(4) pi^2 / 2
@test_approx_eq nball_vol_factor(5) 8pi^2/15
@test_approx_eq nball_vol_factor(6) pi^3 / 6
@test_approx_eq nball_vol_factor(7) 16pi^3/105

# -----------------------------------------------------------------------------
# Simple likelihood function

    # gaussians centered at (1, 1) and (-1, -1) with a width of 0.1.
const mu1 = [1., 1.]
const mu2 = [-1., -1.]
const sigma = 0.1
const invvar = eye(2) / sigma^2

function logl(x)
    dx1 = x .- mu1
    dx2 = x .- mu2
    return logaddexp((dx1' * invvar * dx1)[1] / 2.0,
                     (dx2' * invvar * dx2)[1] / 2.0)
end

# Use a flat prior, over [-5, 5] in both dimensions
prior(x) = 10.0 .* x .- 5.0

srand(0)
res = sample(logl, prior, 2; npoints=100)
@printf "evidence = %6.3f +/- %6.3f\n" res["logz"] res["logzerr"]

#(Approximate) analytic evidence for two identical Gaussian blobs,
# over a uniform prior [-5:5][-5:5] with density 1/100 in this domain:
analytic_logz = log(2.0 * 2.0*pi*sigma*sigma / 100.)
@printf "analytic = %6.3f\n" analytic_logz

    # calculate evidence on fine grid.
    #dx = 0.1
    #xv = np.arange(-5.0 + dx/2., 5., dx)
    #yv = np.arange(-5.0 + dx/2., 5., dx)
    #grid_logz = -1.e300
    #for x in xv:
    #    for y in yv:
    #        grid_logz = np.logaddexp(grid_logz, logl(np.array([x, y])))
    #grid_logz += np.log(dx * dx / 100.)  # adjust for point density
    #print "grid_logz =", grid_logz

    #assert abs(res.logz - analytic_logz) < 2.0 * res.logzerr
    #assert abs(res.logz - grid_logz) < 2.0 * res.logzerr