text:
###
# Optimization problems for VoPP Upper Bounds
###

####Separation oracles
#hbar(t) = M * abs(t-1)/2 - D
#notice no dependence on mu.
#maximum occurs at an end point or 1
function sep_MAD(a, b, lam1, lam2, M, D)
	lhs = lam1 * a + lam2/2 * abs(a-1) * M - lam2 * D
	rhs = lam1 * b + lam2/2 * abs(b-1) * M - lam2 * D
	mid = a <= 1. <= b ? lam1 - lam2 * D : -Inf
	val = max(lhs, rhs, mid)

	if lhs >= val
		return val, a
	elseif rhs >= val
		return val, b
	elseif mid >= val
		return val, 1.
	end
	throw("Unreachable")
	return -1., -1.
end

#hbar(t) = M^2 (t - 1)^2 - C^2
#Again, no dependence on mu
#maximum occurs at an end point or at unique minimizer
function sep_CV(a, b, lam1, lam2, M, C)
	lhs = lam1 * a + lam2 * M^2 * (a - 1)^2 - lam2 * C^2
	rhs = lam1 * b + lam2 * M^2 * (b - 1)^2 - lam2 * C^2
	vstar = 1 - lam1 / 2 / lam2 / M^2
	mid = a <= vstar <= b ? lam1 * vstar + lam2 * M^2 * (vstar - 1)^2 - lam2 * C^2 : -Inf
	val = max(lhs, rhs, mid)

	if lhs >= val
		return val, a
	elseif rhs >= val
		return val, b
	elseif mid >= val
		return val, vstar
	end
	throw("Unreachable")
	return -1., -1.
end

#h(t) = -log(t/mu) + log (B/mu)   equivalent to Exp[log[V]]= log[B]
#hbar(t)  = -log(M(t -1 ) + 1) + log (B/mu)
function sep_GM(a, b, lam1, lam2, M, mu, B)
	lhs = lam1 * a - lam2 * log(M * (a - 1) + 1 ) + lam2 * log(B/mu)
	rhs = lam1 * b - lam2 * log(M * (b - 1) + 1 ) + lam2 * log(B/mu)

	vstar = lam2 <= 0 ? 1 - 1/M + lam2 / lam1 : a - 1e-6
	mid = a <= vstar <= b ? lam1 * vstar - lam2 * log(M * (vstar - 1) + 1) + lam2 * log(B/mu) : -Inf
	val = max(lhs, rhs, mid)

	if lhs >= val
		return val, a
	elseif rhs >= val
		return val, b
	elseif mid >= val
		return val, vstar
	end
	throw("Unreachable")
	return -1., -1.
end

#### Functions based on a generic separation argument
#h 1d function that represents moment constraint (includes moment):  E[h(V)] = 0
#sep_fun(a, b, lam1, lam2, M, mu) maximizes of lam1 * v  + lam2 h(mu * (Mv + 1-M) ).  returns optval, vstar
function _vopp_moment_dual(S, M, mu, h, sep_fun; N, numCuts=100, TOL =1e-6, print_trace=false)
	#translate everything to the Vc space
	hbar(t) = h( mu * (M * t + (1 - M)) )
	Sc = vopp.comp_Sc(S, M)

	@assert isfinite(Sc) "Currenty only supports finite Sc"

	p_grid = range(0, stop=Sc, length=N+1) #indexes from 1, doc indexes from 0
	m = Model(solver=GurobiSolver(OutputFlag = false))
	@variable(m, theta)
	@variable(m, lam[1:2])
	@variable(m, Qs[1:N+1] >= 0)  #indexes form 1, doc indexes from 0

	@constraint(m, sum(Qs) == 1)
	@constraint(m, theta + lam[1] * Sc + lam[2] * hbar(Sc) <= sum([p_grid[i] * Qs[i] for i = 1:N+1]))
	@objective(m, Max, theta + lam[1])

	for iter = 1:numCuts
		iter == numCuts && throw("Max Iteration Count Reached")
		status = solve(m)  #check status?
		feasible = true
		#check if any cut is violated
		for k = 2:N+1
			val, vstar = sep_fun(p_grid[k-1], p_grid[k], getvalue(lam[1]), getvalue(lam[2]), M, mu)
			lhs = val + getvalue(theta)
			rhs = sum([p_grid[j] * getvalue(Qs[j]) for j = 1:k-1])

			if lhs >= rhs + TOL
				@constraint(m, theta + vstar * lam[1] + hbar(vstar) * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )
				feasible = false
				print_trace && println("iter $(iter): \t $lhs \t $rhs")
				continue
			end
		end
		if feasible
			@assert status == :Optimal "Cutting Procedure failed: $(status)"
			return 1/(getobjectivevalue(m) - TOL)
		end
	end
	##unreachable
	return -1.
end

######################
####
# Specialized function for MAD
# Not using separation.
# Algorithm relies on fact that minimizers are always at endpt or 1.
####
function vopp_ub_MAD_opt(S, M, D; N = 500)
	#Solve the two cases for lambda2
	val_neg = _vopp_ub_MAD(S, M, D, false, N)
	val_pos = _vopp_ub_MAD(S, M, D, true, N)

	#take minimum bc we already took reciprocals
	min(val_neg, val_pos), val_neg, val_pos
end

function _vopp_ub_MAD(S, M, D, is_lam2_pos, N)
	Sc = comp_Sc(S, M)
	@assert isfinite(Sc) "Currenty only supports finite Sc"

	p_grid = range(0, stop=Sc, length=N+1) #indexes from 1, doc indexes from 0
	m = Model(solver=GurobiSolver(OutputFlag=false))
	@variable(m, theta)
	@variable(m, lam[1:2])
	@variable(m, Qs[1:N+1] >= 0)  #indexes form 1, doc indexes from 0

	if is_lam2_pos
		@constraint(m, lam[2] >= 0)
	else
		@constraint(m, lam[2] <= 0)
	end

	hbar(t) = M * abs(t-1)/2 - D

	@constraint(m, sum(Qs) == 1)
	@constraint(m, theta + lam[1] * Sc + lam[2] * hbar(Sc) <= sum([p_grid[i] * Qs[i] for i = 1:N+1]))
	@objective(m, Max, theta + lam[1])

	for k = 2:N+1
		@constraint(m, theta + p_grid[k-1] * lam[1] + hbar(p_grid[k-1]) * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )
		@constraint(m, theta + p_grid[k]   * lam[1] + hbar(p_grid[k])   * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )

		if !is_lam2_pos && p_grid[k-1] <= 1 <= p_grid[k]
			@constraint(m, theta + lam[1] + hbar(1.) * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )
		end
	end
	status = solve(m)
	@assert status == :Optimal "Problem Failed $status"
	return 1/(getobjectivevalue(m))
end

### Specialized function for CV
# useSep determines use separator or else Reformulation
# Notice that the minimizer is not always at 1 because of linear offset
#Weirdly slow for some reason.
function vopp_ub_CV(S, M, mu, C; method=:ConstraintGeneration, N=500, numCuts=100, TOL=1e-6, print_trace=false)
	if method == :Slemma
		return _vopp_CV_SLemma(S, M, mu, C, N)
	elseif method == :ConstraintGeneration
		h(v) = (v - mu)^2 / mu^2 - C^2
		sep_fun(a, b, lam1, lam2, M, mu) = sep_CV(a, b, lam1, lam2, M, C)
		return _vopp_moment_dual(S, M, mu, h, sep_fun, N=N, numCuts=numCuts, TOL=TOL, print_trace=print_trace)
	elseif method == :Reformulation
		return min( _vopp_ub_CV(S, M, mu, C, true, N), _vopp_ub_CV(S, M, mu, C, false, N) )
	else
		throw("Method must be one of :Slemma , :ConstraintGeneration , :Reformulation")
	end
end

function _vopp_ub_CV(S, M, mu, C, is_lam2_pos, N=500)
	#translate everything to the Vc space
	hbar(t) = M^2 * (t - 1)^2 - C^2  ##This functional form implicitly coded into dual constraints
	Sc = vopp.comp_Sc(S, M)
	@assert isfinite(Sc) "Currenty only supports finite Sc"

	p_grid = range(0, stop=Sc, length=N+1) #indexes from 1, doc indexes from 0
	m = Model(solver=GurobiSolver(OutputFlag = false))

	@variable(m, theta)
	@variable(m, lam[1:2])

	if is_lam2_pos
		@constraint(m, lam[2] >= 0)  #means solutions only occur at endpoints
	else
		@constraint(m, lam[2] <= 0)
	end

	@variable(m, Qs[1:N+1] >= 0)  #indexes form 1, doc indexes from 0
	@constraint(m, sum(Qs) == 1)

	@constraint(m, theta + lam[1] * Sc + lam[2] * hbar(Sc) <= sum([p_grid[i] * Qs[i] for i = 1:N+1]))

	#Add the auxiliary variables for each constraint in lam_2_neg
	if !is_lam2_pos
		@variable(m, alphas[2:N+1])
		@variable(m, betas[2:N+1])
		@variable(m, ts[2:N+1] >= 0)
	end

	for k = 2:N+1
		#always need the end points
		@constraint(m, theta + p_grid[k-1] * lam[1] + hbar(p_grid[k-1]) * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )
		@constraint(m, theta + p_grid[k] * lam[1] + hbar(p_grid[k]) * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )

		if !is_lam2_pos #concave case
			#use auxiliary variables to specify dual constraints
			@constraint(m, alphas[k] - lam[1] == betas[k])
			@constraint(m, betas[k]^2 <= -4 * M^2 * lam[2] * ts[k])
			@constraint(m, ts[k] - alphas[k] + lam[1] + alphas[k] * p_grid[k-1] <=
								C^2 * lam[2] - theta + sum([p_grid[i] * Qs[i] for i = 1:k-1]))

			@constraint(m, ts[k] - alphas[k] + lam[1] + alphas[k] * p_grid[k] <=
								C^2 * lam[2] - theta + sum([p_grid[i] * Qs[i] for i = 1:k-1]))
		end
	end
	@objective(m, Max, theta + lam[1])

	##Solve and extract solution
	status = solve(m)
	@assert status == :Optimal "Reformulation failed: $(status)"

	return 1/getobjectivevalue(m)
end

function _vopp_CV_SLemma(S, M, mu, C, N)
	#translate everything to the Vc space
	hbar(t) = M^2 * (t - 1)^2 - C^2  ##This functional form implicitly coded into dual constraints
	Sc = vopp.comp_Sc(S, M)
	@assert isfinite(Sc) "Currenty only supports finite Sc"

	ps = range(0, stop=Sc, length=N+1) #indexes from 1, doc indexes from 0
	m = Model(solver=GurobiSolver(OutputFlag = false))

	@variable(m, theta)
	@variable(m, lam[1:2])
	@variable(m, Qs[1:N+1] >= 0)  #doc indexes from 0 to N
	@constraint(m, sum(Qs) == 1)

	@constraint(m, theta + lam[1] * Sc + lam[2] * hbar(Sc) <= sum([ps[i] * Qs[i] for i = 1:N+1]))

	for k = 2:N+1
		#add auxiliary variables and specify modified S-Lemma
		x = @variable(m)
		y = @variable(m)
		z = @variable(m)

		@constraint(m, x >= 0)
		@constraint(m, y >= 0)
		@constraint(m, y >= -lam[2] * M^2)

		@constraint(m, ps[k-1] * ps[k] * (y + lam[2]*M^2) - x
						>=
						theta - sum(ps[j] * Qs[j] for j = 1:k) +
						lam[2] * M^2 - lam[2] * C^2
					)

		@constraint(m, 4 * y * x >= z^2)
		@constraint(m, z == 2 * lam[2] * M^2 - lam[1] -
							(ps[k-1] + ps[k]) * (y + lam[2] * M^2)
					)
	end
	@objective(m, Max, theta + lam[1])

	##Solve and extract solution
	status = solve(m)
	@assert status == :Optimal "Reformulation failed: $(status)"

	return 1/getobjectivevalue(m)
end

### Specialized function for BM
# useSep determines use separator or else Reformulation
function vopp_ub_GM(S, M, mu, B; useSep=false, N=500, numCuts=100, TOL=1e-6, print_trace=false)
	if useSep
		h(v) = -log(v/mu) + log(B/mu)
		sep_fun(a, b, lam1, lam2, M, mu) = sep_GM(a, b, lam1, lam2, M, mu, B)
		return _vopp_moment_dual(S, M, mu, h, sep_fun, N=N, numCuts=numCuts, TOL=TOL, print_trace=print_trace)
	end
	return min( _vopp_ub_GM(S, M, mu, B, true, N), _vopp_ub_GM(S, M, mu, B, false, N) )
end

### Specialized function for GM using convex conjugates
#uses a combination of IPOPT and prayer.
#h(t) = -log(t/mu) + log(B/mu)
function _vopp_ub_GM(S, M, mu, B, is_lam2_pos, N)
	#translate everything to the Vc space
	hbar(t) = -log(M*t + 1- M) + log(B/mu)  ##This functional form implicitly coded into dual constraints

	Sc = vopp.comp_Sc(S, M)
	@assert isfinite(Sc) "Currenty only supports finite Sc"

	p_grid = range(0, stop=Sc, length=N+1) #indexes from 1, doc indexes from 0
	m = Model(solver=IpoptSolver())

	@variable(m, theta)
	@variable(m, lam[1:2])

	if is_lam2_pos
		@constraint(m, lam[2] >= 0)  #means solutions only occur at endpoints
	else
		@constraint(m, lam[2] <= 0)
	end

	@variable(m, Qs[1:N+1] >= 0)  #indexes form 1, doc indexes from 0
	@constraint(m, sum(Qs) == 1)

	@constraint(m, theta + lam[1] * Sc + lam[2] * hbar(Sc) <= sum([p_grid[i] * Qs[i] for i = 1:N+1]))

	#Add the auxiliary variables for each constraint in lam_2_neg
	if !is_lam2_pos
		@variable(m, rhos[2:N+1])
		@variable(m, ts[2:N+1])
	end

	for k = 2:N+1
		#always need the end points
		@constraint(m, theta + p_grid[k-1] * lam[1] + hbar(p_grid[k-1]) * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )
		@constraint(m, theta + p_grid[k]   * lam[1] + hbar(p_grid[k])   * lam[2] <= sum( [p_grid[j] * Qs[j] for j = 1:k-1]) )

		if !is_lam2_pos #concave case
			#use auxiliary variables to specify dual constraints
			@constraint(m, rhos[k] >= lam[1])
			@NLconstraint(m, -lam[2] * log(-M * lam[2] / (rhos[k] - lam[1])) <= ts[k])
			@constraint(m, (1-M)/M * (rhos[k] - lam[1]) + lam[2] + lam[2]*log(B/mu) + ts[k] +
						rhos[k] * p_grid[k-1] <= sum([p_grid[i] * Qs[i] for i = 1:k-1]) - theta)
			@constraint(m, (1-M)/M * (rhos[k] - lam[1]) + lam[2] + lam[2]*log(B/mu) + ts[k] +
						rhos[k] * p_grid[k] <= sum([p_grid[i] * Qs[i] for i = 1:k-1]) - theta)

			# @constraint(m, -lam[2] * log( -rhos[k] / M /lam[2] ) <= ts[k])
			# @constraint(m, lam[1]* p_grid[k-1] - rhos[k] * p_grid[k-1] - lams[2] * log(B/mu) - lams[2] +
			# 				(1 - 1/M) * rhos[k] + ts[k] <= sum([p_grid[i] * Qs[i] for i = 1:k-1]) - theta )
			# @constraint(m, lam[1]* p_grid[k] - rhos[k] * p_grid[k] - lams[2] * log(B/mu) - lams[2] +
			# 				(1 - 1/M) * rhos[k] + ts[k] <= sum([p_grid[i] * Qs[i] for i = 1:k-1]) - theta )
		end
	end
	@objective(m, Max, theta + lam[1])

	##Solve and extract solution
	status = solve(m)
	@assert status == :Optimal "Reformulation failed: $(status)"

	return 1/getobjectivevalue(m)
end

#h(t) = I( V >= phat mu) - q
#hbar(t) = I( t >= v0) - q where v0 = (phat -1)/M + 1
#sep_fun(a, b, lam1, lam2, M, mu) maximizes of lam1 * v  + lam2 hbar(mu * (Mv + 1-M) ).  returns optval, vstar
function sep_IC(a, b, lam1, lam2, M, mu, phat, q; tol=1e-10)
	v0 = (phat - 1)/M + 1
	lhs(v) =  v >= v0 ? lam1 * v + lam2 * (1 - q) : lam1 * v - lam2 * q

	#check end points
	tstar = a
	val = lhs(tstar)

	if lhs(b) > val
		tstar = b
		val = lhs(tstar)
	end

	if a <= v0 < b
		if lam2 >= 0
			#check v0
			if lhs(v0) > val
				tstar = v0
				val = lhs(tstar)
			end
		else #lam2 < 0
			#check v0-tol
			if lhs(v0 -tol) > val
				tstar = v0-tol
				val = lhs(tstar)
			end
		end
	end
	val, tstar
end

### Specialized function for Incumbent Price
function vopp_ub_IC(S, M, mu, phat, q; N=500,numCuts=100, TOL=1e-6, print_trace=false)
	v0 = (phat - 1)/M + 1
	h(v) = v >= v0 ?  1 - q : -q
	sep_fun(a, b, lam1, lam2, M, mu) = sep_IC(a, b, lam1, lam2, M, mu, phat, q)
	return _vopp_moment_dual(S, M, mu, h, sep_fun, N=N, numCuts=numCuts, TOL=TOL, print_trace=print_trace)
end


text:
module CUBLAS

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcublas, unsafe_free!, @retry_reclaim, isdebug

using GPUArrays

using LinearAlgebra

using BFloat16s

using CEnum

# core library
include("libcublas_common.jl")
include("error.jl")
include("libcublas.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

# thread cache for task-local library handles
const thread_handles = Vector{Union{Nothing,cublasHandle_t}}()
const thread_xt_handles = Vector{Union{Nothing,cublasXtHandle_t}}()

function handle()
    tid = Threads.threadid()
    if @inbounds thread_handles[tid] === nothing
        ctx = context()
        thread_handles[tid] = get!(task_local_storage(), (:CUBLAS, ctx)) do
            handle = cublasCreate()
            cublasSetStream_v2(handle, CuStreamPerThread())
            finalizer(current_task()) do task
                CUDA.isvalid(ctx) || return
                context!(ctx) do
                    cublasDestroy_v2(handle)
                end
            end

            # enable tensor math mode if our device supports it, and fast math is enabled
            dev = device()
            if Base.JLOptions().fast_math == 1 && capability(dev) >= v"7.0" && version() >= v"9"
                cublasSetMathMode(handle, CUBLAS_TENSOR_OP_MATH)
            end

            handle
        end
    end
    something(@inbounds thread_handles[tid])
end

function xt_handle()
    tid = Threads.threadid()
    if @inbounds thread_xt_handles[tid] === nothing
        ctx = context()
        thread_xt_handles[tid] = get!(task_local_storage(), (:CUBLASxt, ctx)) do
            handle = cublasXtCreate()
            finalizer(current_task()) do task
                CUDA.isvalid(ctx) || return
                context!(ctx) do
                    cublasXtDestroy(handle)
                end
            end

            # select the devices
            # TODO: this is weird, since we typically use a single device per thread/context
            devs = convert.(Cint, CUDA.devices())
            cublasXtDeviceSelect(handle, length(devs), devs)

            handle
        end
    end
    something(@inbounds thread_xt_handles[tid])
end

function __init__()
    resize!(thread_handles, Threads.nthreads())
    fill!(thread_handles, nothing)

    resize!(thread_xt_handles, Threads.nthreads())
    fill!(thread_xt_handles, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
        thread_xt_handles[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
        thread_xt_handles[tid] = nothing
    end
end

function log_message(cstr)
    str = unsafe_string(cstr)
    print(str)
    # NOTE: we can't `@debug` these messages, because the logging function is called for
    #       every line... also not sure what the i!/I! prefixes mean (info?)
    return
end

function __runtime_init__()
    # enable library logging when launched with JULIA_DEBUG=CUBLAS
    if isdebug(:init, CUBLAS)
        callback = @cfunction(log_message, Nothing, (Cstring,))
        cublasSetLoggerCallback(callback)
    end
end

end

ext: jl


text:
## --- Test `StaticString`s

    # Test StaticString constructors
    str = c"Hello, world! 🌍"
    @test isa(str, StaticString{19})
    @test sizeof(str) == 19
    @test StaticTools.strlen(str) == length(str) == 18
    @test pointer(str) == Base.unsafe_convert(Ptr{UInt8}, str)
    str1 = StaticString(str[1:end])
    @test isa(str1, StaticString{19})

    # Test basic string operations
    @test str == c"Hello, world! 🌍"
    @test str*str == str^2
    @test codeunit(str) === UInt8
    @test codeunit(str, 5) == UInt8('o')
    @test ncodeunits(str) == length(str)+1
    @test codeunits(c"Hello") == codeunits(c"Hello")

    # Test mutability
    str[8] = 'W'
    @test str[8] == 0x57 # W
    str[:] = c"Hello, world! 🌍"
    @test str[8] == 0x77 # w

    # Test indexing
    @test isa(str[1:end], StringView)
    @test str == str[1:end]
    @test str == str[:]
    @test str[1:2] == str[1:2]
    @test str[1:2] != str[1:3]
    @test str == copy(str)

    # Test ascii escaping
    many_escapes = c"\"\0\a\b\f\n\r\t\v\'\"\\"
    @test isa(many_escapes, StaticString{13})
    @test length(many_escapes) == 12
    @test all(codeunits(many_escapes) .== codeunits("\"\0\a\b\f\n\r\t\v'\"\\\0"))

    # Test consistency with other strings
    abc = c"abc"
    @test abc == "abc"
    @test "abc" == abc
    @test abc == c"abc"
    @test abc == abc[1:3]
    @test abc[1:3] == "abc"


text:
using JuMP
using Gurobi
model = Model(Gurobi.Optimizer)
# set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ON)
# set_optimizer_attribute(model, "maxNodes", 5)
month_num = 8
product_num = 2
month = 1:month_num
product = 1:product_num
# paramerters
q = [5000, 3000] # product fix cost
p = [100, 60] # cost per product
d = [400 400 800 800 1200 1200 1200 1200;
     200 200 200 200 200  200  500  500] # demand of product by month
s_init = [200,0] # initial stock
h = [5,3] # stock cost per product month
capacity = 1500

@variable(model, x[i in product, j in month] >= 0)
@variable(model, s[i in product, j in month] >= 0)
@variable(model, y[i in product, j in month], Bin)

# 库存平衡
@constraint(model ,dem_sat1[i in product], s_init[i] + x[i, 1] == d[i, 1] + s[i, 1])
@constraint(model, dem_sat[i in product, j in month[2:end]], s[i, j - 1] + x[i, j] == d[i, j] + s[i, j] )
# 最大产量约束
@constraint(model, vub[i in product, j in month], x[i, j] <= sum(d[i, k]  for k in j:month_num) * y[i, j])
@constraint(model, max_cap[j in month], sum(x[i, j] for i in product) <= capacity)
# 目标函数
@expression(model ,cost, sum(p[i] * x[i, j] + q[i] * y[i, j] for i in product, j in month) + sum(h[i] * s[i, j] for i in product, j in month[1:end - 1]) + sum(h[i] * s[i, month[end]] / 2 for i in product))
@objective(model, Min, cost)
optimize!(model)
# 显示结果
cost_val = value.(cost)
x_val = value.(x)
s_val = value.(s)
y_val = value.(y)
@show cost_val;
@show x_val;
@show y_val;
@show s_val;

max_stars_repo_name: Behemoth-s/ProductPlanningByMIP

max_issues_repo_name: Behemoth-s/ProductPlanningByMIP

max_forks_repo_name: Behemoth-s/ProductPlanningByMIP


text:
"""
Placeholder for a short summary about {PKGNAME}.
"""
module {PKGNAME}

end # module


text:
# converters from selected RSEXPREC to Hash
# They are used to translate SEXPREC attributes into Hash

using Dates

function Base.convert(::Type{Hash}, pl::RPairList)
    res = Hash()
    for i in eachindex(pl.items)
        @inbounds setindex!(res, pl.items[i], pl.tags[i])
    end
    res
end

##############################################################################
##
## Conversion of intermediate R objects into Vector{T} and DataFrame objects
##
##############################################################################

isna(x::Int32) = x == R_NA_INT32
isna(x::Float64) = isna_float64(reinterpret(UInt64, x))
# if re or im is NA, the whole complex number is NA
isna(x::ComplexF64) = isna(real(x)) || isna(imag(x))

# convert R vector into Vector holding elements of type T
# if force_missing is true, the result is always Vector{Union{T,Missing}},
# otherwise it's Vector{T} if `rv` doesn't contain NAs
function jlvec(::Type{T}, rv::RVEC, force_missing::Bool=true) where T <: Number
    anyna = any(isna, rv.data)
    if force_missing || anyna
        res = convert(Vector{Union{T,Missing}}, rv.data)
        if anyna
            @inbounds for (i,x) in enumerate(rv.data)
                isna(x) && (res[i] = missing)
            end
        end
        return res
    else
        return convert(Vector{T}, rv.data)
    end
end

# convert R nullable vector (has an explicit NA mask) into Vector{T[?]}
function jlvec(::Type{T}, rv::RNullableVector{R}, force_missing::Bool=true) where {T, R}
    anyna = any(rv.na)
    if force_missing || anyna
        res = convert(Vector{Union{T,Missing}}, rv.data)
        anyna && @inbounds res[rv.na] .= missing
        return res
    else
        return convert(Vector{T}, rv.data)
    end
end

# convert R vector into Vector of appropriate type
function jlvec(rv::RVEC, force_missing::Bool=true)
    if inherits(rv, R_Date_Class)
        return jlvec(Dates.Date, rv, force_missing)
    elseif inherits(rv, R_POSIXct_Class)
        return jlvec(ZonedDateTime, rv, force_missing)
    else
        return jlvec(isconcretetype(eltype(rv.data)) ? eltype(rv.data) : Any, rv, force_missing)
    end
end

# convert R logical vector (uses Int32 to store values) into Vector{Bool[?]}
function jlvec(rl::RLogicalVector, force_missing::Bool=true)
    anyna = any(isna, rl.data)
    if force_missing || anyna
        return Union{Bool,Missing}[ifelse(isna(x), missing, x != 0) for x in rl.data]
    else
        return Bool[x != 0 for x in rl.data]
    end
end

# kernel method that converts Vector{Int32} into Vector{R} replacing R_NA_INT32 with 0
# it's assumed that v fits into R
na2zero(::Type{R}, v::Vector{Int32}) where R =
    [ifelse(!isna(x), x % R, zero(R)) for x in v]

# convert to CategoricalVector{String[?]} if `ri` is a factor,
# to Dates.Date is it's a vector of dates and
# to Vector{Int32[?]} otherwise
function jlvec(ri::RIntegerVector, force_missing::Bool=true)
    if isfactor(ri)
        return jlvec(CategoricalArray, ri, force_missing)
    elseif inherits(ri, R_Date_Class)
        return jlvec(Dates.Date, ri, force_missing)
    else
        return jlvec(eltype(ri.data), ri, force_missing)
    end
end

# convert R factor to CategoricalArray
function jlvec(::Type{CategoricalArray}, ri::RVEC, force_missing::Bool=true)
    @assert isfactor(ri)

    rlevels0 = getattr(ri, "levels")
    sz0 = length(rlevels0)
    # CategoricalArrays#v0.8 does not allow duplicate levels
    rlevels = unique(rlevels0)
    sz = length(rlevels)
    hasduplicates = sz0 != sz

    REFTYPE = sz0 <= typemax(UInt8)  ? UInt8 :
              sz0 <= typemax(UInt16) ? UInt16 :
              sz0 <= typemax(UInt32) ? UInt32 :
                                      UInt64
    refs = na2zero(REFTYPE, ri.data)

    if hasduplicates
        # map refs with dups to unique refs
        ref_map = REFTYPE.(indexin(rlevels0, rlevels))
        @inbounds for i in eachindex(refs)
            ref = refs[i]
            refs[i] = ref == 0 ? 0 : ref_map[ref]
        end
        @warn "Dropped duplicate factor levels"
    end

    anyna = any(iszero, refs)
    pool = CategoricalPool{String, REFTYPE}(rlevels, inherits(ri, "ordered"))
    if force_missing || anyna
        return CategoricalArray{Union{String, Missing}, 1}(refs, pool)
    else
        return CategoricalArray{String, 1}(refs, pool)
    end
end

# convert R Date to Dates.Date
function jlvec(::Type{Dates.Date}, rv::RVEC, force_missing::Bool=true)
    @assert inherits(rv, R_Date_Class)
    nas = isnan.(rv.data)
    if force_missing || any(nas)
        dates = Union{Dates.Date, Missing}[isna ? missing : rdays2date(dtfloat)
                 for (isna, dtfloat) in zip(nas, rv.data)]
    else
        dates = rdays2date.(rv.data)
    end
    return dates
end

# convert R POSIXct to ZonedDateTime
function jlvec(::Type{ZonedDateTime}, rv::RVEC, force_missing::Bool=true)
    @assert class(rv) == R_POSIXct_Class
    tz, validtz = getjuliatz(rv)
    nas = isnan.(rv.data)
    if force_missing || any(nas)
        datetimes = Union{ZonedDateTime, Missing}[isna ? missing : _unix2zdt(dtfloat, tz=tz)
                     for (isna, dtfloat) in zip(nas, rv.data)]
    else
        datetimes = _unix2zdt.(rv.data, tz=tz)
    end
    return datetimes
end

function simplify_eltype(v::AbstractVector)
    isconcretetype(eltype(v)) && return eltype(v)

    eltypes = Set{DataType}()
    eltyp = Union{}
    try
        for x in v
            xtyp = typeof(x)
            if !(xtyp in eltypes)
                # promote only if arrays have the same eltype, otherwise return Any
                if (eltyp === Union{}) ||
                   (nonmissingtype(xtyp) === nonmissingtype(eltyp)) ||
                   (xtyp <: AbstractArray && eltyp <: AbstractArray &&
                    nonmissingtype(eltype(xtyp)) === nonmissingtype(eltype(eltyp)))
                    eltyp = promote_type(eltyp, xtyp)
                else
                    return Any
                end
                if !isconcretetype(eltyp) || (eltyp === Any)
                    return Any
                end
                push!(eltypes, xtyp)
            end
        end
    catch e # missing promotion rule
        return Any
    end
    return eltyp
end

# generic vector conversion
function jlvec(::Type{T}, rv::RVEC, force_missing::Bool=true) where T
    res = sexp2julia.(rv.data)
    if !isconcretetype(eltype(res))
        return convert(Vector{simplify_eltype(res)}, res)
    else
        return res
    end
end

function sexp2julia(rex::RSEXPREC)
    @warn "Conversion of $(typeof(rex)) to Julia is not implemented" maxlog=1
    return nothing
end

function sexp2julia(rv::RVEC)
    # TODO dimnames?
    # FIXME add force_missing option to control whether always convert to Union{T, Missing}
    jv = jlvec(rv, false)
    if hasnames(rv)
        # if data has no NA, convert to simple Vector
        return DictoVec(jv, names(rv))
    else
        hasdims = hasdim(rv)
        if !hasdims && length(rv.data)==1
            # scalar
            return jv[1]
        elseif !hasdims
            # vectors
            return jv
        else
            # matrices and so on
            dims = tuple(convert(Vector{Int}, getattr(rv, "dim"))...)
            return reshape(jv, dims)
        end
    end
end

function sexp2julia(rl::RList)
    if isdataframe(rl)
        # FIXME add force_missing option to control whether always convert to Union{T, Missing}
        DataFrame(Any[isa(col, RAltRep) ? sexp2julia(col) : jlvec(col, false) for col in rl.data],
                  identifier.(names(rl)), makeunique=true)
    elseif hasnames(rl)
        DictoVec(jlvec(Any, rl), names(rl))
    else
        # FIXME return DictoVec if forceDictoVec is on
        jlvec(Any, rl)
    end
end

function sexp2julia(ar::RAltRep)
    if iswrapped(ar)
        return sexp2julia(unwrap(ar))
    elseif iscompactseq(ar)
        return jlrange(ar)
    else
        @warn unsupported_altrep_message(ar)
        return nothing
    end
end

function rdays2date(days::Real)
    epoch_conv = 719528 # Dates.date2epochdays(Date("1970-01-01"))
    Dates.epochdays2date(days + epoch_conv)
end

# gets R timezone from the data attribute and converts it to TimeZones.TimeZone
# see r2juliatz()
function getjuliatz(rv::RVEC, deftz=tz"UTC")
    tzattr = getattr(rv, "tzone", [""])[1]
    if tzattr == ""
        return deftz, true # R will store a blank for tzone
    else
        return r2juliatz(tzattr, deftz)
    end
end

# converts R timezone code to TimeZones.TimeZone
# returns a tuple:
#  - timezone (or `deftz` if `rtz` is not recognized as a valid time zone)
#  - boolean flag: true if `rtz` is not recognized, false otherwise
function r2juliatz(rtz::AbstractString, deftz=tz"UTC")
    valid = istimezone(rtz)
    if !valid
        @warn "Could not determine the timezone of '$(rtz)', treating as '$deftz'" maxlog=1
        return deftz, false
    else
        return TimeZone(rtz), true
    end
end

# version with user-specified tz (unix2zdt(seconds) is fixed to tz"UTZ")
_unix2zdt(seconds::Real; tz::TimeZone=tz"UTC") =
    ZonedDateTime(Dates.unix2datetime(seconds), tz, from_utc=true)


text:
# From user `Lazaro` on Slack:
set_theme!(
    rowgap = 5,
    colgap = 5,
    fontsize = 20,
    backgroundcolor = :black, # "#212946",
    cycle = :color,
    palette = (color = COLORS,),
    Axis = (
        xgridcolor = :white,
        ygridcolor = :white,
        backgroundcolor = :black, # "#212946",
        xgridwidth=0.35,
        ygridwidth=0.35,
        xtickalign=0,
        ytickalign=0,
        xtickcolor = :white,
        ytickcolor = :white,
        bottomspinecolor = :white,
        topspinecolor =:white,
        leftspinecolor =:white,
        rightspinecolor =:white,
        xlabelcolor = :white,
        ylabelcolor = :white,
        xlabelsize = 24,
        ylabelsize = 24,
        xticklabelcolor = :white,
        yticklabelcolor = :white,
        xgridvisible = true,
        ygridvisible = true,
    ),
)

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
"""
    getNullCursor() -> CXCursor
Return the "NULL" CXCursor.
"""
getNullCursor() = clang_getNullCursor()

"""
    isNull(c::Union{CXCursor,CLCursor}) -> Bool
Return true if cursor is null.
"""
isNull(c::Union{CXCursor,CLCursor}) = clang_Cursor_isNull(c) != 0

# equality
Base.:(==)(c1::CXCursor, c2::CXCursor)::Bool = clang_equalCursors(c1, c2)
Base.hash(c::CXCursor)::UInt = clang_hashCursor(c)
Base.:(==)(c1::CLCursor, c2::CLCursor) = c1.cursor == c2.cursor
Base.hash(x::CLCursor) = hash(x.cursor)
Base.:(==)(c1::CLCursor, c2::CXCursor) = c1.cursor == c2
Base.:(==)(c1::CXCursor, c2::CLCursor) = c1 == c2.cursor

"Free and convert a CXString to String. Note this function will free the given CXString so ensure it is not called twice."
function _cxstring_to_string(cxstr::CXString)
    ptr = clang_getCString(cxstr)
    ptr == C_NULL && return ""
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    kind(c::CXCursor) -> CXCursorKind
Return the kind of the given cursor.
"""
kind(c::CXCursor) = clang_getCursorKind(c)

"""
    kind(c::CLCursor) -> CXCursorKind
Return the kind of the given cursor.
Note this method directly reads CXCursor's `kind` field, which won't invoke additional
`clang_getCursorKind` function calls.
"""
kind(c::CLCursor) = c.cursor.kind

"""
    isDeclaration(k::CXcursorKind) -> Bool
    isDeclaration(c::CLCursor) -> Bool
Return true if the given cursor kind represents a declaration.
"""
isDeclaration(k::CXCursorKind)::Bool = clang_isDeclaration(k)
isDeclaration(c::CLCursor) = isDeclaration(kind(c))

"""
    isInvalidDeclaration(x::CXcursorKind) -> Bool
Return true if the given declaration is invalid.
A declaration is invalid if it could not be parsed successfully.
"""
isInvalidDeclaration(x::Union{CXCursor,CLCursor})::Bool = clang_isInvalidDeclaration(x)

"""
    isReference(k::CXcursorKind) -> Bool
    isReference(c::CLCursor) -> Bool
Return true if the given cursor kind represents a simple reference. Note that other kinds of
cursors (such as expressions) can also refer to other cursors. Use [`getCursorReferenced`](@ref) to determine
whether a particular cursor refers to another entity.
"""
isReference(k::CXCursorKind)::Bool = clang_isReference(k)
isReference(c::CLCursor) = isReference(kind(c))

"""
    isExpression(k::CXcursorKind) -> Bool
    isExpression(c::CLCursor) -> Bool
Return true if the given cursor kind represents an expression.
"""
isExpression(k::CXCursorKind)::Bool = clang_isExpression(k)
isExpression(c::CLCursor) = isExpression(kind(c))

"""
    isStatement(k::CXcursorKind) -> Bool
    isStatement(c::CLCursor) -> Bool
Return true if the given cursor kind represents a statement.
"""
isStatement(k::CXCursorKind)::Bool = clang_isStatement(k)
isStatement(c::CLCursor) = isstat(kind(c))

"""
    isAttribute(k::CXcursorKind) -> Bool
    isAttribute(c::CLCursor) -> Bool
Return true if the given cursor kind represents an attribute.
"""
isAttribute(k::CXCursorKind)::Bool = clang_isAttribute(k)
isAttribute(c::CLCursor) = isAttribute(kind(c))

"""
    hasAttrs(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether the given cursor has any attributes.
"""
hasAttrs(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_hasAttrs(c)

"""
    isInvalid(k::CXcursorKind) -> Bool
    isInvalid(c::CLCursor) -> Bool
Return true if the given cursor kind represents an valid cursor.
"""
isInvalid(k::CXCursorKind) = !(Bool(clang_isInvalid(k)))
isInvalid(c::CLCursor) = isInvalid(kind(c))

"""
    isTranslationUnit(k::CXcursorKind) -> Bool
    isTranslationUnit(c::CLCursor) -> Bool
Return true if the given cursor kind represents a translation unit.
"""
isTranslationUnit(k::CXCursorKind)::Bool = clang_isTranslationUnit(k)
isTranslationUnit(c::CLCursor) = isTranslationUnit(kind(c))

"""
    isPreprocessing(k::CXcursorKind) -> Bool
    isPreprocessing(c::CLCursor) -> Bool
Return true if the given cursor kind represents a preprocessing element, such as a
preprocessor directive or macro instantiation.
"""
isPreprocessing(k::CXCursorKind)::Bool = clang_isPreprocessing(k)
isPreprocessing(c::CLCursor) = isPreprocessing(kind(c))

"""
    isUnexposed(k::CXcursorKind) -> Bool
Return true if the given cursor kind represents a currently unexposed piece of the AST
(e.g., CXCursor_UnexposedStmt).
"""
isUnexposed(k::CXCursorKind)::Bool = clang_isUnexposed(k)
isUnexposed(c::CLCursor) = isUnexposed(kind(c))

"""
    getCursorLinkage(c::Union{CXCursor,CLCursor}) -> CXLinkageKind
Return the linkage of the entity referred to by a given cursor.
"""
getCursorLinkage(c::Union{CXCursor,CLCursor})::CXLinkageKind = clang_getCursorLinkage(c)

## TODO:
# clang_getCursorVisibility
# clang_getCursorAvailability
# clang_getCursorPlatformAvailability
# clang_disposeCXPlatformAvailability

"""
    getCursorLanguage(c::Union{CXCursor,CLCursor}) -> CXLanguageKind
Return the language of the entity referred to by a given cursor.
Return `CXLanguage_Invalid` if the input cursor is not a decl.
Note that, this function has limitations, for example, it cannot distinguish C++ structs
from C structs, in both cases, it returns `CXLanguage_C`.
"""
getCursorLanguage(c::Union{CXCursor,CLCursor})::CXLanguageKind = clang_getCursorLanguage(c)

# clang_getCursorTLSKind

"""
    getTranslationUnit(c::Union{CXCursor,CLCursor}) -> CXTranslationUnit
Returns the translation unit that a cursor originated from.
"""
getTranslationUnit(c::Union{CXCursor,CLCursor}) = clang_Cursor_getTranslationUnit(c)

## TODO: this is mainly for internal debugging?
# clang_createCXCursorSet
# clang_disposeCXCursorSet
# clang_CXCursorSet_contains
# clang_CXCursorSet_insert

"""
    getCursorSemanticParent(c::CXCursor) -> CXCursor
    getCursorSemanticParent(c::CLCursor) -> CLCursor
Return the semantic parent of the given cursor. Please checkout libclang's doc to know more.
"""
getCursorSemanticParent(c::CXCursor) = clang_getCursorSemanticParent(c)
getCursorSemanticParent(c::CLCursor)::CLCursor = clang_getCursorSemanticParent(c)

"""
    getCursorLexicalParent(c::CXCursor) -> CXCursor
    getCursorLexicalParent(c::CLCursor) -> CLCursor
Return the lexical parent of the given cursor. Please checkout libclang's doc to know more.
"""
getCursorLexicalParent(c::CXCursor) = clang_getCursorLexicalParent(c)
getCursorLexicalParent(c::CLCursor)::CLCursor = clang_getCursorLexicalParent(c)

## TODO:
# clang_getOverriddenCursors
# clang_disposeOverriddenCursors

"""
    getIncludedFile(c::Union{CXCursor,CLCursor}) -> CXFile
Return the file that is included by the given inclusion directive cursor.
"""
getIncludedFile(c::Union{CXCursor,CLCursor}) = clang_getIncludedFile(c)

"""
    getCursorLocation(c::Union{CXCursor,CLCursor}) -> CXSourceLocation
Return the physical location of the source constructor referenced by the given cursor.
"""
getCursorLocation(c::Union{CXCursor,CLCursor}) = clang_getCursorLocation(c)

"""
    getCursorExtent(c::Union{CXCursor,CLCursor}) -> CXSourceRange
Return the physical extent of the source construct referenced by the given cursor.

The extent of a cursor starts with the file/line/column pointing at the first character
within the source construct that the cursor refers to and ends with the last character
within that source construct. For a declaration, the extent covers the declaration itself.
For a reference, the extent covers the location of the reference (e.g., where the referenced
entity was actually used).
"""
getCursorExtent(c::Union{CXCursor,CLCursor}) = clang_getCursorExtent(c)

"""
    getCursorType(c::CXCursor) -> CXType
    getCursorType(c::CLCursor) -> CLType
Return the type of a CXCursor (if any). To get the cursor from a type, see [`getTypeDeclaration`](@ref).
"""
getCursorType(c::CXCursor) = clang_getCursorType(c)
getCursorType(c::CLCursor)::CLType = clang_getCursorType(c)

"""
    getTypedefDeclUnderlyingType(c::CLTypedefDecl) -> CLType
Return the underlying type of a typedef declaration.
"""
getTypedefDeclUnderlyingType(c::CXCursor) = clang_getTypedefDeclUnderlyingType(c)
getTypedefDeclUnderlyingType(c::CLTypedefDecl)::CLType = clang_getTypedefDeclUnderlyingType(c)

"""
    getEnumDeclIntegerType(c::CLEnumDecl) -> CLType
Retrieve the integer type of an enum declaration.
"""
getEnumDeclIntegerType(c::CXCursor) = clang_getEnumDeclIntegerType(c)
getEnumDeclIntegerType(c::CLEnumDecl)::CLType = clang_getEnumDeclIntegerType(c)

"""
    value(c::CLCursor) -> Integer
Return the integer value of an enum constant declaration.
"""
function value(c::CLEnumConstantDecl)::Integer
    typeKind = kind(getCursorType(c))
    if typeKind == CXType_Int || typeKind == CXType_Long || typeKind == CXType_LongLong
        return clang_getEnumConstantDeclValue(c)
    elseif typeKind == CXType_UInt ||
           typeKind == CXType_ULong ||
           typeKind == CXType_ULongLong
        return clang_getEnumConstantDeclUnsignedValue(c)
    end
    return error("Unknown EnumConstantDecl type: ", typeKind, " cursor: ", kind(c))
end

"""
    getFieldDeclBitWidth(c::CLFieldDecl) -> Int
Return the bit width of a bit field declaration as an integer.
"""
getFieldDeclBitWidth(c::CXCursor)::Int = clang_getFieldDeclBitWidth(c)
getFieldDeclBitWidth(c::CLFieldDecl) = clang_getFieldDeclBitWidth(c)

"""
    getNumArguments(c::Union{CXCursor,CLFunctionDecl,CLCXXMethod}) -> Int
Return the number of non-variadic arguments associated with a given cursor.
"""
getNumArguments(c::Union{CXCursor,CLFunctionDecl,CLCXXMethod})::Int = clang_Cursor_getNumArguments(c)

"""
    getArgument(c::CXCursor, i::Integer) -> CXCursor
    getArgument(c::Union{CLFunctionDecl,CLCXXMethod}, i::Integer) -> CLCursor
Return the argument cursor of a function or method.
"""
getArgument(c::CXCursor, i::Integer) = clang_Cursor_getArgument(c, Unsigned(i))
getArgument(c::Union{CLFunctionDecl,CLCXXMethod}, i::Integer)::CLCursor =
    clang_Cursor_getArgument(c, Unsigned(i))

## TODO:
# clang_Cursor_getNumTemplateArguments
# clang_Cursor_getTemplateArgumentKind
# clang_Cursor_getTemplateArgumentType
# clang_Cursor_getTemplateArgumentValue
# clang_Cursor_getTemplateArgumentUnsignedValue

"""
    isMacroFunctionLike(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether a CXCursor that is a macro, is function like.
"""
isMacroFunctionLike(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isMacroFunctionLike(c)

"""
    isMacroBuiltin(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether a  CXCursor that is a macro, is a builtin one.
"""
isMacroBuiltin(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isMacroBuiltin(c)

"""
    isFunctionInlined(c::Union{CXCursor,CLCursor}) -> Bool
Determine whether a CXCursor that is a function declaration, is an inline declaration.
"""
isFunctionInlined(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isFunctionInlined(c)

"""
    getCursorResultType(c::CXCursor) -> CXType
    getCursorResultType(c::Union{CLFunctionDecl,CLCXXMethod}) -> CLType
Return the return type associated with a given cursor. This only returns a valid type if
the cursor refers to a function or method.
"""
getCursorResultType(c::CXCursor) = clang_getCursorResultType(c)
getCursorResultType(c::Union{CLFunctionDecl,CLCXXMethod})::CLType = clang_getCursorResultType(c)

## TODO:
# clang_getCursorExceptionSpecificationType
# clang_Cursor_getOffsetOfField
#

"""
    isAnonymous(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the given cursor represents an anonymous record declaration(C++).
"""
isAnonymous(c::Union{CXCursor,CLCursor})::Bool = clang_Cursor_isAnonymous(c)

"""
    isBitField(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the cursor specifies a Record member that is a bitfield.
"""
isBitField(c::Union{CXCursor,CLCursor}) = clang_Cursor_isBitField(c) != 0

"""
    isVirtualBase(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the base class specified by the cursor with kind CX_CXXBaseSpecifier is virtual.
"""
isVirtualBase(c::Union{CXCursor,CLCursor})::Bool = clang_isVirtualBase(c)

## TODO:
# clang_getCXXAccessSpecifier
# clang_Cursor_getStorageClass
# clang_getNumOverloadedDecls
# clang_getOverloadedDecl
# clang_getIBOutletCollectionType
#

"""
    spelling(c::Union{CXCursor,CLCursor}) -> String
Return a name for the entity referenced by this cursor.
"""
function spelling(c::Union{CXCursor,CLCursor})
    cxstr = clang_getCursorSpelling(c)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

# clang_Cursor_getSpellingNameRange

"""
    name(c::Union{CXCursor,CLCursor}) -> String
Return the display name for the entity referenced by this cursor.
"""
function name(c::Union{CXCursor,CLCursor})
    cxstr = clang_getCursorDisplayName(c)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    getCursorReferenced(c::CXCursor) -> CXCursor
    getCursorReferenced(c::CLCursor) -> CLCursor
For a cursor that is a reference, retrieve a cursor representing the entity that it references.
"""
getCursorReferenced(c::CXCursor) = clang_getCursorReferenced(c)
getCursorReferenced(c::CLCursor)::CLCursor = clang_getCursorReferenced(c)

"""
    getCursorDefinition(c::CXCursor) -> CXCursor
    getCursorDefinition(c::CLCursor) -> CLCursor
For a cursor that is either a reference to or a declaration of some entity, retrieve a cursor
that describes the definition of that entity.
"""
getCursorDefinition(c::CXCursor) = clang_getCursorDefinition(c)
getCursorDefinition(c::CLCursor)::CLCursor = clang_getCursorDefinition(c)

"""
    isCursorDefinition(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the declaration pointed to by this cursor is also a definition of that entity.
"""
isCursorDefinition(c::Union{CXCursor,CLCursor})::Bool = clang_isCursorDefinition(c)

"""
    getCanonicalType(c::CXCursor) -> CXCursor
    getCanonicalType(c::CLCursor) -> CLCursor
Return the getCanonicalType cursor corresponding to the given cursor.
"""
getCanonicalType(c::CXCursor) = clang_getCanonicalCursor(c)
getCanonicalType(c::CLCursor)::CLCursor = clang_getCanonicalCursor(c)

## TODO:
# clang_Cursor_getObjCSelectorIndex
# clang_Cursor_isDynamicCall
# clang_Cursor_getReceiverType

"""
    isVariadic(c::Union{CXCursor,CLCursor}) -> Bool
Return true if the given cursor is a variadic function or method.
"""
isVariadic(c::Union{CXCursor,CLCursor}) = clang_Cursor_isVariadic(c) != 0

# clang_Cursor_isExternalSymbol
# clang_Cursor_getCommentRange
getRawCommentText(c::Union{CXCursor, CLCursor}) = _cxstring_to_string(clang_Cursor_getRawCommentText(c))
getBriefCommentText(c::Union{CXCursor, CLCursor}) = _cxstring_to_string(clang_Cursor_getBriefCommentText(c))

function spelling(k::CXCursorKind)
    cxstr = clang_getCursorKindSpelling(k)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    get_filename(x::CXFile) -> String
Return the complete file and path name of the given file
"""
function get_filename(x::CXFile)
    cxstr = clang_getFileName(x)
    ptr = clang_getCString(cxstr)
    ptr == C_NULL && return ""
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    get_filename(c::Union{CXCursor,CLCursor}) -> String
Return the complete file and path name of the given file referenced by the input cursor.
"""
function get_filename(c::Union{CXCursor,CLCursor})
    file = Ref{CXFile}(C_NULL)
    location = clang_getCursorLocation(c)
    clang_getExpansionLocation(location, file, Ref{Cuint}(0), Ref{Cuint}(0), Ref{Cuint}(0))
    if file[] != C_NULL
        str = GC.@preserve file get_filename(file[])
        return str
    else
        return ""
    end
end

"""
    get_file_line_column(c::Union{CXCursor,CLCursor}) -> (String, Int, Int)
Return file name, line and column number.
"""
function get_file_line_column(c::Union{CXCursor,CLCursor})
    file = Ref{CXFile}(C_NULL)
    line = Ref{Cuint}(0)
    column = Ref{Cuint}(0)
    offset = Ref{Cuint}(0)
    location = clang_getCursorLocation(c)
    clang_getExpansionLocation(location, file, line, column, offset)
    if file[] != C_NULL
        cxstr = clang_getFileName(file[])
        ptr = clang_getCString(cxstr)
        s = unsafe_string(ptr)
        clang_disposeString(cxstr)
        return s, Int(line[]), Int(column[])
    else
        return "", 0, 0
    end
end

"""
    search(cursors::Vector{CLCursor}, ismatch::Function) -> Vector{CLCursor}
Return vector of CLCursors that match predicate. `ismatch` is a function that accepts a CLCursor argument.
"""
function search(cursors::Vector{T}, ismatch::Function) where {T<:CLCursor}
    matched = CLCursor[]
    for cursor in cursors
        ismatch(cursor) && push!(matched, cursor)
    end
    return matched
end
search(c::CLCursor, ismatch::Function) = search(children(c), ismatch)
search(c::CLCursor, s::String) = search(c, x -> spelling(x) == s)
search(c::CLCursor, k::CXCursorKind) = search(c, x -> kind(x) == k)

# visitor
function cu_children_visitor(cursor::CXCursor, parent::CXCursor, list)::Cuint
    push!(list, cursor)
    return CXChildVisit_Continue
end

"""
    children(cursor::CXCursor) -> Vector{CXCursor}
    children(cursor::CLCursor) -> Vector{CLCursor}
Return immediate cursors of the given cursor.
"""
function children(cursor::CXCursor)
    list = CXCursor[]
    cu_visitor_cb = @cfunction(
        cu_children_visitor, Cuint, (CXCursor, CXCursor, Ref{Vector{CXCursor}})
    )
    GC.@preserve list ccall(
        (:clang_visitChildren, LibClang.libclang),
        UInt32,
        (CXCursor, CXCursorVisitor, Any),
        cursor,
        cu_visitor_cb,
        list,
    )
    return list
end
children(c::CLCursor)::Vector{CLCursor} = children(c.cursor)

"""
    get_function_args(cursor::CLCursor) -> Vector{CLCursor}
Return function arguments for a given cursor.
"""
get_function_args(cursor::CLCursor) = search(cursor, CXCursor_ParmDecl)

"""
    is_typedef_anon(current::CLCursor, next::CLCursor) -> Bool
Return true if the current cursor is a typedef anonymous struct/enum.
"""
function is_typedef_anon(current::CLCursor, next::CLCursor)
    !isempty(name(current)) && return false
    refback = children(next)
    if kind(next) == CXCursor_TypedefDecl && !isempty(refback) && isempty(name(refback[1]))
        return true
    else
        return false
    end
end

"""
    is_forward_declaration(x::CLCursor) -> Bool
Return true if the cursor is a forward declaration.
Reference: https://joshpeterson.github.io/identifying-a-forward-declaration-with-libclang
"""
function is_forward_declaration(x::CLCursor)
    def = getCursorDefinition(x)
    def == getNullCursor() && return true
    return !(x == def)
end

"""
    is_inclusion_directive(x::CLCursor) -> Bool
Return true if the cursor is an inclusion directive.
"""
function is_inclusion_directive(x::CLCursor)
    k = kind(x)
    return k == CXCursor_InclusionDirective && k == CXCursor_LastPreprocessing
end

## Comment utilities

"""
    getParsedComment(c::CXCursor) -> CXComment
    getParsedComment(c::CLCursor) -> CLComment
Return the parsed AST of doxygen comment associated with the given node.
"""
getParsedComment(c::CXCursor) = clang_Cursor_getParsedComment(c)
getParsedComment(c::CLCursor)::CLComment = getParsedComment(c.cursor)

"""
    getAsHTML(c::Union{CXComment, FullComment})
Convert a comment into an HTML fragment. Check [libclang's documentation](https://clang.llvm.org/doxygen/group__CINDEX__COMMENT.html) for more details.
"""
getAsHTML(c::Union{CXComment, FullComment}) = _cxstring_to_string(clang_FullComment_getAsHTML(c))

"""
    getAsXML(c::Union{CXComment, FullComment})
Convert a comment into an XML document.
"""
getAsXML(c::Union{CXComment, FullComment}) = _cxstring_to_string(clang_FullComment_getAsXML(c))

"""
    kind(c::Union{CXComment, CXComment}) -> CXCommentKind
Return the kind of a comment.
"""
kind(c::Union{CXComment, CXComment}) = clang_Comment_getKind(c)

"""
    isWhiteSpace(c::Union{CXComment, CLComment}) -> Bool
Return whether a CXParagraph or CXText is contains only white space. Return false for other kinds of comment.
"""
isWhiteSpace(c::Union{CXComment, CLComment}) = !iszero(clang_Comment_isWhitespace(c))
hasTrailingNewline(c::Union{CXComment, CLComment}) = !iszero(clang_InlineContentComment_hasTrailingNewline(c))

"""
    getText(c::Text) -> String
    getText(c::VerbatimBlockLine) -> String
    getText(c::VerbatimLine) -> String
Get the text content assiciated with the node.
"""
getText(c::Text) = _cxstring_to_string(clang_TextComment_getText(c))
getText(c::VerbatimBlockLine) = _cxstring_to_string(clang_VerbatimBlockLineComment_getText(c))
getText(c::VerbatimLine) = _cxstring_to_string(clang_VerbatimLineComment_getText(c))

"""
    getCommandName(c::InlineCommand) -> String
    getCommandName(c::BlockCommand) -> String
Return the command name of a command, e.g. "brief" for `\\brief blah`.
"""
getCommandName(c::InlineCommand) = _cxstring_to_string(clang_InlineCommandComment_getCommandName(c))
getCommandName(c::BlockCommand) = _cxstring_to_string(clang_BlockCommandComment_getCommandName(c))

"""
    getParamName(c::ParamCommand) -> String
Return the name of a parameter.
"""
getParamName(c::ParamCommand) = _cxstring_to_string(clang_ParamCommandComment_getParamName(c))
getParamName(c::TParamCommand) = _cxstring_to_string(clang_TParamCommandComment_getParamName(c))

# See getArguments
getNumArgs(c::InlineCommand) = clang_InlineCommandComment_getNumArgs(c)
getNumArgs(c::BlockCommand) = clang_BlockCommandComment_getNumArgs(c)

# See getArguments
getArgText(c::InlineCommand, i) = _cxstring_to_string(clang_InlineCommandComment_getArgText(c, i))
getArgText(c::BlockCommand, i) = _cxstring_to_string(clang_BlockCommandComment_getArgText(c, i))

"""
    getParagraph(c::CXComment)
    getParagraph(c::BlockCommand)
Return the paragraph argument of a given command.
"""
getParagraph(c::CXComment) = clang_BlockCommandComment_getParagraph(c)
getParagraph(c::BlockCommand)::CLComment = getParagraph(c.comment)

"""
    getRenderKind(c::Union{CXComment, InlineCommand}) -> CXCommentInlineCommandRenderKind
Return the most appropriate rendering mode, chosen on command semantics in Doxygen.
"""
getRenderKind(c::Union{CXComment, InlineCommand}) = clang_InlineCommandComment_getRenderKind(c)

"""
    getTagName(c::Union{CXComment, HTMLStartTag, HTMLEndTag})
Return the tag name of a HTML tag.
"""
getTagName(c::Union{CXComment, HTMLStartTag, HTMLEndTag}) = _cxstring_to_string(clang_HTMLTagComment_getTagName(c))
"""
    isSelfClosing(c::Union{CXComment, HTMLStartTag})
Return whether a tag is self-closing (for example, `<br />`).
"""
isSelfClosing(c::Union{CXComment, HTMLStartTag}) = !iszero(clang_HTMLStartTagComment_isSelfClosing(c))

# See getAttributes
getNumAttrs(c::Union{CXComment, HTMLStartTag}) = clang_HTMLStartTag_getNumAttrs(c)

# See getAttributes
getAttrName(c::Union{CXComment, HTMLStartTag}, i) = _cxstring_to_string(clang_HTMLStartTag_getAttrName(c, i))

# See getAttributes
getAttrValue(c::Union{CXComment, HTMLStartTag}, i) = _cxstring_to_string(clang_HTMLStartTag_getAttrName(c, i))

"""
    getAsString(c::Union{CXComment, HTMLStartTag, HTMLEndTag})
Convert an HTML tag to its string representation.
"""
getAsString(c::Union{CXComment, HTMLStartTag, HTMLEndTag}) = _cxstring_to_string(clang_HTMLTagComment_getAsString(c))

"""
    getIndex(c::Union{CXComment, ParamCommand})
Return zero-based parameter index in function prototype or `missing` if the parameter is invalid.
"""
function getIndex(c::Union{CXComment, ParamCommand})
    if iszero(clang_ParamCommandComment_isParamIndexValid(c))
        return missing
    else
        return clang_ParamCommandComment_getParamIndex(c)
    end
end

"""
    getDirection(c::Union{CXComment, ParamCommand})
Return the parameter passing direction (`[in]`, `[out]`, `[in,out]`) or `missing` if the direction is not specified.
"""
function getDirection(c::Union{CXComment, ParamCommand})
    if iszero(clang_ParamCommandComment_isDirectionExplicit(c))
        return missing
    else
        return clang_ParamCommandComment_getDirection(c)
    end
end

getNumChildren(c::Union{CXComment, CLComment}) = clang_Comment_getNumChildren(c)

getChild(c::CXComment, i) = clang_Comment_getChild(c, i)
getChild(c::CLComment, i)::CLComment = getChild(c.comment, i)

function children(c::CXComment)
    num = getNumChildren(c)
    result = Vector{CXComment}(undef, num)
    for i in 1:num
        result[i] = getChild(c, i-1)
    end
    result
end

function children(c::CLComment)
    num = getNumChildren(c)
    result = Vector{CLComment}(undef, num)
    for i in 1:num
        result[i] = getChild(c, i-1)
    end
    result
end

"""
    getArguments(c::Union{BlockCommand, InlineCommand}) -> Vector{String}
Return the arguments of a command.
"""
function getArguments(c::Union{BlockCommand, InlineCommand})
    num = getNumArgs(c)
    result = Vector{String}(undef, num)
    for i in 1:num
        result[i] = getArgText(c, i-1)
    end
    return result
end

"""
    getAttributes(c::Union{CXComment, HTMLStartTag}) -> Vector{Pair{String, String}}
Return the attributes of an HTML tag.
"""
function getAttributes(c::Union{CXComment, HTMLStartTag})
    num = getNumAttrs(c)
    result = Vector{Pair{String, String}}(undef, num)
    for i in 1:num
        result[i] = (getAttrName(c, i-1) => getAttrValue(c, i-1))
    end
    return result
end

# clang_TParamCommandComment_isParamPositionValid
# clang_TParamCommandComment_getDepth
# clang_TParamCommandComment_getIndex

"""
    getSourceCode(cursor::Union{CXCursor, CLCursor}) -> String
Get source code for the cursor.
"""
function getSourceCode(cursor::Union{CXCursor, CLCursor})
    ext = getCursorExtent(cursor)
    range_begin = clang_getRangeStart(ext)
    range_end = clang_getRangeEnd(ext)
    file_begin = Ref{CXFile}(C_NULL)
    file_end = Ref{CXFile}(C_NULL)
    offset_begin = Ref{Cuint}(0)
    offset_end = Ref{Cuint}(0)
    clang_getExpansionLocation(range_begin, file_begin, C_NULL, C_NULL, offset_begin)
    clang_getExpansionLocation(range_end, file_end, C_NULL, C_NULL, offset_end)
    @assert !iszero(clang_File_isEqual(file_begin[], file_end[]))
    tu = getTranslationUnit(cursor)
    size = Ref{Csize_t}(0)
    buf = clang_getFileContents(tu, file_begin[], size)
    @assert offset_begin[] <= offset_end[] <= size[]
    unsafe_string(pointer(buf) + offset_begin[], offset_end[] - offset_begin[])
end


text:
# Genome explicit Metacommunity Model - GeMM
#
# Ludwig Leidinger <l.leidinger@gmx.net>,
# Daniel Vedder <daniel.vedder@stud-mail.uni-wuerzburg.de>
# (c) 2017-2020, licensed under the terms of the MIT license
#
# Contains all custom types and functions necessary for an individual based,
# genome and niche explicit, eco-evo-env metacommunity model.

module GeMM

using
    ArgParse,
    Dates,
    Distributions,
    Random,
    Statistics

export
    rungemm

include("defaults.jl")

include("entities.jl")

include("input.jl")

include("output.jl")

include("constants.jl")

include("initialisation.jl")

include("genetics.jl")

include("reproduction.jl")

include("dispersal.jl")

include("survival.jl")

include("habitatchange.jl")

include("invasion.jl")

include("zosterops.jl")

include("scheduling.jl")

include("run_simulation.jl")

end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
#versioninfo()

#import Pkg; Pkg.add("Arpack")
#import Pkg; Pkg.add("ArgParse")

using LinearAlgebra
using Arpack
using SparseArrays
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--N"
            help = "set Nsize (should be >=4)"
            arg_type = Int
            default = 8
    end
    return parse_args(s)
end

function make_spin()
    S0 = sparse([1,2],[1,2],[1,1])
    Sx = sparse([1,2],[2,1],[1,1])
    Sy = sparse([1,2],[2,1],[-1im,1im])
    Sz = sparse([1,2],[1,2],[1,-1])
    return S0, Sx, Sy, Sz
end

function make_interaction_list(N,Nbond)
    list_site1 = [i for i in 1:N]
    list_site2 = [mod(i,N)+1 for i in 1:N]
    list_Jxx = ones(N)
    list_Jyy = ones(N)
    list_Jzz = ones(N)
    return list_site1, list_site2, list_Jxx, list_Jyy, list_Jzz
end

function make_hamiltonian(S0,Sx,Sy,Sz,N,Nbond,
    list_site1,list_site2,list_Jxx,list_Jyy,list_Jzz)
    Ham = spzeros(2^N,2^N)
    for bond in 1:Nbond
        i1 = list_site1[bond]
        i2 = list_site2[bond]
        Jxx = list_Jxx[bond]
        Jyy = list_Jyy[bond]
        Jzz = list_Jzz[bond]
        SxSx = 1
        SySy = 1
        SzSz = 1
        for site in 1:N
            if (site==i1 || site==i2)
                SxSx = kron(SxSx,Sx)
                SySy = kron(SySy,Sy)
                SzSz = kron(SzSz,Sz)
            else
                SxSx = kron(SxSx,S0)
                SySy = kron(SySy,S0)
                SzSz = kron(SzSz,S0)
            end
#        println(bond,site,SxSx)
#        println(bond,site,SySy)
#        println(bond,site,SzSz)
        end
#        Ham += Jxx * SxSx + Jyy * SySy + Jzz * SzSz
        Ham += real(Jxx * SxSx + Jyy * SySy + Jzz * SzSz)
    end
#    println(Ham)
    return Ham
end

function main()
    parsed_args = parse_commandline()
#    println("Parsed args:")
#    for (arg,val) in parsed_args
#        println("  $arg  =>  $val")
#    end

    N = parsed_args["N"]
    Nbond = N
    println("N=",N)

    S0, Sx, Sy, Sz = make_spin()
    println("S0=",Matrix(S0))
    println("Sx=",Matrix(Sx))
    println("Sy=",Matrix(Sy))
    println("Sz=",Matrix(Sz))

    list_site1, list_site2, list_Jxx, list_Jyy, list_Jzz = make_interaction_list(N,Nbond)
    println("list_site1=",list_site1)
    println("list_site2=",list_site2)
    println("list_Jxx=",list_Jxx)
    println("list_Jyy=",list_Jyy)
    println("list_Jzz=",list_Jzz)

    @time Ham = make_hamiltonian(S0,Sx,Sy,Sz,N,Nbond,list_site1,list_site2,list_Jxx,list_Jyy,list_Jzz)

    @time ene,vec = eigs(Ham, nev=5, which=:LM)
    println("# energy:",ene[1]," ",ene[2]," ",ene[3]," ",ene[4]," ",ene[5])
end

main()


text:
module TestOptimTensorBoard

using Optim
using OptimTensorBoard
using Test

logdir = joinpath(@__DIR__, "tmp", "logdir")
rm(joinpath(@__DIR__, "tmp"); force=true, recursive=true)

f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
function g!(G, x)
    G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    G[2] = 200.0 * (x[2] - x[1]^2)
end
x0 = [0.0, 0.0]
calledwith = []
callback(tr) = (push!(calledwith, tr); false)

@testset "savearrays = true" begin
    empty!(calledwith)
    opt = Optim.Options(
        extended_trace = true,
        callback = optimcallback(
            logdir,
            callback = callback,
            savearrays = true,
        ))
    optimize(f, g!, x0, LBFGS(), opt)
    @test !isempty(calledwith)
end

@testset "savearrays = :image" begin
    empty!(calledwith)
    opt = Optim.Options(
        extended_trace = true,
        callback = optimcallback(
            logdir,
            callback = callback,
            savearrays = :image,
        ))
    optimize(f, g!, x0, LBFGS(), opt)
    @test !isempty(calledwith)
end

end  # module


text:
struct LibSixelEncoder <: AbstractSixelEncoder
    # (Experimental) internal fields

    # We need `allocator` to constructor libsixel objects, however, users of this
    # Julia package is not expected to use this field as it really should just live
    # in the C world.
    allocator::SixelAllocator

    function LibSixelEncoder(allocator=SixelAllocator())
        new(allocator)
    end
end

function sixel_write_callback_function(buffer_ptr::Ptr{Cchar}, sz::Cint, priv::Ref{T})::Cint where {T<:IO}
    unsafe_write(priv[], buffer_ptr, sz)
end

# libsixel backend requires contiguous memeory layout, to avoid unexpected bugs
# we limit ourself to `Matrix` type and let the frontend API `sixel_encode` transforms
# other fancy array types into `Matrix`.
function (enc::LibSixelEncoder)(io::T, bytes::Matrix) where {T<:IO}
    # colorbits = default_colorbits(bytes)
    pixelformat = default_pixelformat(bytes)
    quality_mode = default_quality_mode(bytes)

    height, width = size(bytes)
    depth = 3 # unused

    dither = SixelDither(bytes, height, width, pixelformat, quality_mode; allocator=enc.allocator)

    fn_write_cb = @cfunction(sixel_write_callback_function, Cint, (Ptr{Cchar}, Cint, Ref{T}))
    output = SixelOutput(fn_write_cb, Ref{T}(io); allocator=enc.allocator)

    status = C.sixel_encode(bytes, height, width, depth, dither.ptr, output.ptr)
    check_status(status)

    return nothing
end

"""
    default_pixelformat(img)
    default_pixelformat(CT)

Infer the default libsixel `pixelformat` for image `img` or colorant type `CT`.
"""
default_pixelformat(::AbstractArray{CT}) where CT<:Union{Colorant, Real} = default_pixelformat(CT)
default_pixelformat(::Type{RGB{N0f8}})   = C.SIXEL_PIXELFORMAT_RGB888
default_pixelformat(::Type{BGR{N0f8}})   = C.SIXEL_PIXELFORMAT_BGR888
default_pixelformat(::Type{ARGB{N0f8}})  = C.SIXEL_PIXELFORMAT_ARGB8888
default_pixelformat(::Type{RGBA{N0f8}})  = C.SIXEL_PIXELFORMAT_RGBA8888
default_pixelformat(::Type{ABGR{N0f8}})  = C.SIXEL_PIXELFORMAT_ABGR8888
default_pixelformat(::Type{BGRA{N0f8}})  = C.SIXEL_PIXELFORMAT_BGRA8888
default_pixelformat(::Type{Gray{N0f8}})  = C.SIXEL_PIXELFORMAT_G8
default_pixelformat(::Type{T}) where T<:Union{N0f8, UInt8}  = C.SIXEL_PIXELFORMAT_G8
# RGB555, RGB565, BGR555, BGR565 are not provided by ColorTypes
# https://github.com/JuliaGraphics/ColorTypes.jl/issues/114
# Gray2, Gray4 are not used in JuliaImages
# TODO: Gray1, i.e., Gray{Bool} is not supported yet.

"""
    default_quality_mode(img)
    default_quality_mode(CT)

Infer the default quality mode that used to encode pixel.
"""
default_quality_mode(::AbstractArray{CT}) where CT<:Union{Colorant, Real} = default_quality_mode(CT)
default_quality_mode(::Type{CT}) where CT<:AbstractRGB = C.SIXEL_QUALITY_AUTO
# CHECK: highcolor is needed for iTerm2 on macOS, check if this works on other terminal
default_quality_mode(::Type{CT}) where CT<:AbstractGray = C.SIXEL_QUALITY_HIGHCOLOR
default_quality_mode(::Type{T}) where T<:Union{N0f8, UInt8} = C.SIXEL_QUALITY_HIGHCOLOR

canonical_sixel_eltype(::LibSixelEncoder, ::Type{CT}) where CT<:Colorant = n0f8(CT)
canonical_sixel_eltype(::LibSixelEncoder, ::Type{CT}) where CT<:Color3 = RGB{N0f8}
canonical_sixel_eltype(::LibSixelEncoder, ::Type{T}) where T<:Real = N0f8
canonical_sixel_eltype(::LibSixelEncoder, ::Type{T}) where T<:Integer = UInt8
# strip the alpha channel before sending into libsixel
canonical_sixel_eltype(lib::LibSixelEncoder, ::Type{CT}) where CT<:Union{ColorAlpha, AlphaColor} =
    canonical_sixel_eltype(lib, base_color_type(CT))
# TODO: these special types might have native libsixel support, but I haven't
#       figured it out yet.
# canonical_sixel_eltype(::LibSixelEncoder, ::Type{Bool}) = Gray{N0f8}
# canonical_sixel_eltype(::LibSixelEncoder, ::Type{Gray{Bool}}) = Gray{N0f8}
# canonical_sixel_eltype(::LibSixelEncoder, ::Gray24) = Gray{N0f8}
# canonical_sixel_eltype(::LibSixelEncoder, ::RGB24) = RGB{N0f8}
# canonical_sixel_eltype(::LibSixelEncoder, ::ARGB32) = ARGB{N0f8}
# TODO: For unknown reasons, AGray and GrayA encoded by libsixel is not correctly displayed
#       in iTerm. Thus here we convert it to `ARGB` types.
# canonical_sixel_eltype(::LibSixelEncoder, ::Union{AGray, GrayA}) = ARGB{N0f8}


text:
# workaround for opnorm(::AbstractVector, Inf) not being implemented
_opnormInf(B) = opnorm(B, Inf)
_opnormInf(B::AbstractVector) = norm(B, Inf)

# we only use _opnormInf for control flow, so avoid differentiating through it
ChainRulesCore.@non_differentiable _opnormInf(B)

# TODO: replace with estimate
# see https://github.com/JuliaLang/julia/pull/39058
opnormest1(A) = opnorm(A, 1)


text:
using NumericExtensions
function main()
	N = big(100)
	n = sum([1:N])^2 - sumsq([1:N])
	println(n)
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
import Base: size, length, ndims, eltype, first, last, ==, parent
import Base: getindex, setindex!, iterate, eachindex, IndexStyle, CartesianIndices, tail, copyto!, fill!
using LinearAlgebra
using TupleTools
using Adapt

"size of a single symmetric group with Nsym dimensions and size Nt per dimension"
symgrp_size(Nt,Nsym) = binomial_simple(Nt-1+Nsym, Nsym)
symgrp_size(Nt,::Val{Nsym}) where Nsym = binomial_unrolled(Nt+(Nsym-1),Val(Nsym))

# calculates the length of a SymArray
symarrlength(Nts,Nsyms) = prod(symgrp_size.(Nts,Nsyms))

@generated function _getNts(::Val{Nsyms},size::NTuple{N,Int}) where {Nsyms,N}
    @assert sum(Nsyms)==N
    symdims = cumsum(collect((1,Nsyms...)))
    Nts = [ :( size[$ii] ) for ii in symdims[1:end-1]]
    code = quote
        Nts = ($(Nts...),)
    end
    err = :( error("SymArray: sizes $size not compatible with symmetry numbers $Nsyms") )
    for ii = 1:length(Nsyms)
        dd = [ :( size[$jj] ) for jj=symdims[ii]:symdims[ii+1]-1 ]
        push!(code.args,:( all(($(dd...),) .== Nts[$ii]) || $err ))
    end
    push!(code.args, :( Nts ))
    code
end

struct SymArray{Nsyms,T,N,M,datType<:AbstractArray} <: AbstractArray{T,N}
    data::datType
    size::NTuple{N,Int}
    Nts::NTuple{M,Int}
    function SymArray{Nsyms,T}(::Type{arrType},size::Vararg{Int,N}) where {Nsyms,T,N,arrType}
        Nts = _getNts(Val(Nsyms),size)
        M = length(Nsyms)
        data = arrType{T,M}(undef,symgrp_size.(Nts,Nsyms)...)
        new{Nsyms,T,N,M,typeof(data)}(data,size,Nts)
    end
    function SymArray{Nsyms,T}(size::Vararg{Int,N}) where {Nsyms,T,N}
        SymArray{Nsyms,T}(Array,size...)
    end
    # this creates a SymArray that serves as a view on an existing array
    function SymArray{Nsyms}(data::datType,size::Vararg{Int,N}) where {Nsyms,N,datType<:AbstractArray{T}} where T
        Nts = _getNts(Val(Nsyms),size)
        @assert Base.size(data) == symgrp_size.(Nts,Val.(Nsyms))
        new{Nsyms,T,N,length(Nsyms),datType}(data,size,Nts)
    end
end

parent(A::SymArray) = A.data
size(A::SymArray) = A.size
length(A::SymArray) = length(A.data)

# this is necessary for CUDAnative kernels, but also generally useful
# e.g., adapt(CuArray,S) will return a copy of the SymArray with storage in a CuArray
Adapt.adapt_structure(to, x::SymArray{Nsyms}) where Nsyms = SymArray{Nsyms}(adapt(to,x.data),x.size...)

symgrp_size(S::SymArray{Nsyms}) where Nsyms = symgrp_size.(S.Nts,Nsyms)
symgrp_size(S::SymArray{Nsyms},d::Integer) where Nsyms = symgrp_size(S.Nts[d],Nsyms[d])
symgrps(S) = symgrps(typeof(S))
symgrps(::Type{<:SymArray{Nsyms}}) where Nsyms = Nsyms
nsymgrps(S) = nsymgrps(typeof(S))
nsymgrps(::Type{<:SymArray{Nsyms,T,N,M}}) where {Nsyms,T,N,M} = M

"""
    storage_type(A)

Return the type of the underlying storage array for array wrappers.
"""
storage_type(A) = storage_type(typeof(A))
storage_type(T::Type) = (P = parent_type(T); P===T ? T : storage_type(P))

parent_type(T::Type{<:AbstractArray}) = T
parent_type(::Type{<:PermutedDimsArray{T,N,perm,iperm,AA}}) where {T,N,perm,iperm,AA} = AA
parent_type(::Type{<:LinearAlgebra.Transpose{T,S}}) where {T,S} = S
parent_type(::Type{<:LinearAlgebra.Adjoint{T,S}}) where {T,S} = S
parent_type(::Type{<:SubArray{T,N,P}}) where {T,N,P} = P
parent_type(::Type{<:SymArray{Nsyms,T,N,M,datType}}) where {Nsyms,T,N,M,datType} = datType

copyto!(S::SymArray,A::AbstractArray) = begin
    Ainds, Sinds = LinearIndices(A), LinearIndices(S)
    isempty(Ainds) || (checkbounds(Bool, Sinds, first(Ainds)) && checkbounds(Bool, Sinds, last(Ainds))) || throw(BoundsError(S, Ainds))
    @inbounds for (i,I) in zip(Ainds,eachindex(S))
        S[i] = A[I]
    end
    S
end
copyto!(A::AbstractArray,S::SymArray) = begin
    Ainds, Sinds = LinearIndices(A), LinearIndices(S)
    isempty(Sinds) || (checkbounds(Bool, Ainds, first(Sinds)) && checkbounds(Bool, Ainds, last(Sinds))) || throw(BoundsError(A, Sinds))
    @inbounds for (i,I) in zip(Ainds,CartesianIndices(A))
        A[i] = S[I]
    end
    A
end
copyto!(S::SymArray,Ssrc::SymArray) = begin
    @assert symgrps(S) == symgrps(Ssrc)
    @assert size(S) == size(Ssrc)
    copyto!(S.data,Ssrc.data)
    S
end

Base.similar(src::SymArray{Nsyms}) where Nsyms = SymArray{Nsyms}(similar(parent(src)),size(src)...)
Base.copy(src::SymArray{Nsyms}) where Nsyms = SymArray{Nsyms}(copy(parent(src)),size(src)...)

fill!(S::SymArray,v) = fill!(S.data,v)

==(S1::SymArray,S2::SymArray) = (symgrps(S1),S1.data) == (symgrps(S2),S2.data)

SymArray{Nsyms}(A::AbstractArray{T}) where {Nsyms,T} = (S = SymArray{Nsyms,T}(size(A)...); copyto!(S,A))
# to avoid ambiguity with Vararg "view" constructor above
SymArray{(1,)}(A::AbstractVector{T}) where {T} = (S = SymArray{(1,),T}(size(A)...); copyto!(S,A))

"`SymArr_ifsym(A,Nsyms)` make a SymArray if there is some symmetry (i.e., any of the Nsyms are not 1)"
SymArr_ifsym(A,Nsyms) = all(Nsyms.==1) ? A : SymArray{(Nsyms...,)}(A)

"calculates the contribution of index idim in (i1,...,idim,...,iN) to the corresponding linear index for the group"
symind2ind(i,::Val{dim}) where dim = binomial_unrolled(i+(dim-2),Val(dim))

"calculates the linear index corresponding to the symmetric index group (i1,...,iNsym)"
@inline @generated function symgrp_sortedsub2ind(I::Vararg{T,Nsym})::T where {Nsym,T<:Integer}
    terms2toN = [ :( symind2ind(I[$dim],Val($dim)) ) for dim=2:Nsym ]
    :( +(I[1],$(terms2toN...)) )
end

@generated _sub2grp(A::SymArray{Nsyms,T,N}, I::Vararg{Int,N}) where {Nsyms,T,N} = begin
    result = []
    ii::Int = 0
    for Nsym in Nsyms
        Ilocs = ( :( I[$(ii+=1)] ) for _=1:Nsym)
        push!(result, :( symgrp_sortedsub2ind(TupleTools.sort(($(Ilocs...),))...) ) )
    end
    code = :( ($(result...),) )
    code
end

IndexStyle(::Type{<:SymArray}) = IndexCartesian()
getindex(A::SymArray, i::Int) = A.data[i]
getindex(A::SymArray{Nsyms,T,N}, I::Vararg{Int,N}) where {Nsyms,T,N} = (@boundscheck checkbounds(A,I...); @inbounds A.data[_sub2grp(A,I...)...])
setindex!(A::SymArray, v, i::Int) = A.data[i] = v
setindex!(A::SymArray{Nsyms,T,N}, v, I::Vararg{Int,N}) where {Nsyms,T,N} = (@boundscheck checkbounds(A,I...); @inbounds A.data[_sub2grp(A,I...)...] = v);

eachindex(S::SymArray) = CartesianIndices(S)
CartesianIndices(S::SymArray) = SymArrayIter(S)

@generated function lessnexts(::SymArray{Nsyms}) where Nsyms
    lessnext = ones(Bool,sum(Nsyms))
    istart = 1
    for Nsym in Nsyms
        istart += Nsym
        lessnext[istart-1] = false
    end
    Tuple(lessnext)
end

struct SymArrayIter{N}
    lessnext::NTuple{N,Bool}
    sizes::NTuple{N,Int}
    "create an iterator that gives i1<=i2<=i3 etc for each index group"
    SymArrayIter(A::SymArray{Nsyms,T,N,M}) where {Nsyms,T,N,M} = new{N}(lessnexts(A),A.size)
end
ndims(::SymArrayIter{N}) where N = N
eltype(::Type{SymArrayIter{N}}) where N = NTuple{N,Int}
first(iter::SymArrayIter) = CartesianIndex(map(one, iter.sizes))
last(iter::SymArrayIter) = CartesianIndex(iter.sizes...)

@inline function iterate(iter::SymArrayIter)
    iterfirst = first(iter)
    iterfirst, iterfirst
end
@inline function iterate(iter::SymArrayIter, state)
    valid, I = __inc(state.I, iter.sizes, iter.lessnext)
    valid || return nothing
    return CartesianIndex(I...), CartesianIndex(I...)
end
# increment post check to avoid integer overflow
@inline __inc(::Tuple{}, ::Tuple{}, ::Tuple{}) = false, ()
@inline function __inc(state::Tuple{Int}, size::Tuple{Int}, lessnext::Tuple{Int})
    valid = state[1] < size[1]
    return valid, (state[1]+1,)
end

@inline function __inc(state, sizes, lessnext)
    smax = lessnext[1] ? state[2] : sizes[1]
    if state[1] < smax
        return true, (state[1]+1, tail(state)...)
    end
    valid, I = __inc(tail(state), tail(sizes), tail(lessnext))
    return valid, (1, I...)
end

struct SymIndexIter{Nsym}
    size::Int
    "create an iterator that gives i1<=i2<=i3 etc for one index group"
    SymIndexIter(Nsym,size) = new{Nsym}(size)
end
ndims(::SymIndexIter{Nsym}) where Nsym = Nsym
eltype(::Type{SymIndexIter{Nsym}}) where Nsym = NTuple{Nsym,Int}
length(iter::SymIndexIter{Nsym}) where Nsym = symgrp_size(iter.size,Nsym)
first(iter::SymIndexIter{Nsym}) where Nsym = ntuple(one,Val(Nsym))
last(iter::SymIndexIter{Nsym}) where Nsym = ntuple(i->iter.size,Val(Nsym))

@inline function iterate(iter::SymIndexIter)
    I = first(iter)
    I, I
end
@inline function iterate(iter::SymIndexIter, state)
    valid, I = __inc(state, iter.size)
    ifelse(valid, (I, I), nothing)
end
# increment post check to avoid integer overflow
@inline __inc(::Tuple{}, ::Tuple{}) = false, ()
@inline function __inc(state::Tuple{Int}, size::Int)
    valid = state[1] < size
    return valid, (state[1]+1,)
end
@inline function __inc(state::NTuple{N,Int}, size::Int) where N
    if state[1] < state[2]
        return true, (state[1]+1, tail(state)...)
    end
    valid, I = __inc(tail(state), size)
    return valid, (1, I...)
end

function _find_symind(ind::T, ::Val{dim}, high::T) where {dim,T<:Integer}
    dim==1 ? ind+one(T) : searchlast_func(ind, x->symind2ind(x,Val(dim)),one(T),high)
end

"""convert a linear index for a symmetric index group into a group of subindices"""
@generated function ind2sub_symgrp(SI::SymIndexIter{N},ind::T) where {N,T<:Integer}
    code = quote
        ind -= 1
    end
    kis = Symbol.(:k,1:N)
    for dim=N:-1:2
        push!(code.args,:( $(kis[dim]) = _find_symind(ind,Val($dim),T(SI.size)) ))
        push!(code.args,:( ind -= symind2ind($(kis[dim]),Val($dim)) ))
    end
    push!(code.args, :( k1 = ind + 1 ))
    push!(code.args, :( return ($(kis...),)) )
    #display(code)
    code
end

@generated function _grp2sub(A::SymArray{Nsyms,T,N,M}, I::Vararg{Int,M}) where {Nsyms,T,N,M}
    exs = [:( ind2sub_symgrp(SymIndexIter($Nsym,A.Nts[$ii]),I[$ii]) ) for (ii,Nsym) in enumerate(Nsyms)]
    code = :( TupleTools.flatten($(exs...)) )
    #display(code)
    code
end

@inline ind2sub(A, ii) = Tuple(CartesianIndices(A)[ii])
@inline ind2sub(A::SymArray,ii) = _grp2sub(A,ind2sub(A.data,ii)...)


text:
function SeisSemblance(in,param=Dict())

	dt = get(param,"dt",0.002)
	dx = get(param,"dx",10)
	nt,nx = size(in)
	offset = get(param,"offset",collect(0:1:nx-1)*dx)
	vmin = get(param,"vmin",1500)
	vmax = get(param,"vmax",3500)
	nv = get(param,"nv",nx)
	v = linspace(vmin,vmax,nv)
	tau = collect(0:4:nt-1)*dt
	ntau = length(tau)
	L = get(param,"L",5)

	S = zeros(ntau,nv)
	for it = 1:ntau
		for iv = 1:nv
			t = sqrt( tau[it].^2 + (offset[iv]/v[iv]).^2 )
			s = zeros(2*L+1,nx)
			for ig = -L:L
				ts = t + (ig-1)*dt
				for ix = 1:nx
					is = ts/dt+1
					i1 = convert(Int64,floor(is))
					i2 = i1 + 1
					if (i1>=1 && i2<=nt)
						a = is-i1
						s[ig+L+1,ix] = (1.-a)*in[i1,ix] + a*in[i2,ix]   #Grab sample with Linear interpolation
					end
				end
			end
			#s = s.*H
			s1  = sum( (sum(s,2)).^2)
			s2  = sum( sum(s.^2))
			S[it,iv] = abs(s1-s2)
		end
	end
	S = S/maximum(S[:])

	return S,tau,v
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using Oceananigans.Diagnostics

for closure in closures
    @eval begin
        using Oceananigans.TurbulenceClosures: $closure
    end
end

datatuple(args, names) = NamedTuple{names}(a.data for a in args)

function test_closure_instantiation(FT, closurename)
    closure = getproperty(TurbulenceClosures, closurename)(FT)
    return eltype(closure) == FT
end

function test_constant_isotropic_diffusivity_basic(T=Float64; ν=T(0.3), κ=T(0.7))
    closure = ConstantIsotropicDiffusivity(T; κ=(T=κ, S=κ), ν=ν)
    return closure.ν == ν && closure.κ.T == κ
end

function test_constant_isotropic_diffusivity_fluxdiv(FT=Float64; ν=FT(0.3), κ=FT(0.7))
          arch = CPU()
       closure = ConstantIsotropicDiffusivity(FT, κ=(T=κ, S=κ), ν=ν)
          grid = RegularCartesianGrid(FT, size=(3, 1, 4), extent=(3, 1, 4))
    velocities = VelocityFields(arch, grid)
       tracers = TracerFields(arch, grid, (:T, :S))

    u, v, w = velocities
       T, S = tracers

    for k in 1:4
        interior(u)[:, 1, k] .= [0, -1, 0]
        interior(v)[:, 1, k] .= [0, -2, 0]
        interior(w)[:, 1, k] .= [0, -3, 0]
        interior(T)[:, 1, k] .= [0, -1, 0]
    end

    state = (velocities=datatuple(velocities), tracers=datatuple(tracers), diffusivities=nothing)
    fill_halo_regions!(merge(velocities, tracers), arch, nothing, state)

    U, C = datatuples(velocities, tracers)

    return (   ∇_κ_∇c(2, 1, 3, grid, closure, C.T, Val(1)) == 2κ &&
            ∂ⱼ_2ν_Σ₁ⱼ(2, 1, 3, grid, closure, U) == 2ν &&
            ∂ⱼ_2ν_Σ₂ⱼ(2, 1, 3, grid, closure, U) == 4ν &&
            ∂ⱼ_2ν_Σ₃ⱼ(2, 1, 3, grid, closure, U) == 6ν )
end

function test_anisotropic_diffusivity_fluxdiv(FT=Float64; νh=FT(0.3), κh=FT(0.7), νv=FT(0.1), κv=FT(0.5))
          arch = CPU()
       closure = ConstantAnisotropicDiffusivity(FT, νh=νh, νv=νv, κh=(T=κh, S=κh), κv=(T=κv, S=κv))
          grid = RegularCartesianGrid(FT, size=(3, 1, 4), extent=(3, 1, 4))
           eos = LinearEquationOfState(FT)
      buoyancy = SeawaterBuoyancy(FT, gravitational_acceleration=1, equation_of_state=eos)
    velocities = VelocityFields(arch, grid)
       tracers = TracerFields(arch, grid, (:T, :S))

    u, v, w, T, S = merge(velocities, tracers)

    interior(u)[:, 1, 2] .= [0,  1, 0]
    interior(u)[:, 1, 3] .= [0, -1, 0]
    interior(u)[:, 1, 4] .= [0,  1, 0]

    interior(v)[:, 1, 2] .= [0,  1, 0]
    interior(v)[:, 1, 3] .= [0, -2, 0]
    interior(v)[:, 1, 4] .= [0,  1, 0]

    interior(w)[:, 1, 2] .= [0,  1, 0]
    interior(w)[:, 1, 3] .= [0, -3, 0]
    interior(w)[:, 1, 4] .= [0,  1, 0]

    interior(T)[:, 1, 2] .= [0,  1, 0]
    interior(T)[:, 1, 3] .= [0, -4, 0]
    interior(T)[:, 1, 4] .= [0,  1, 0]

    state = (velocities=datatuple(velocities), tracers=datatuple(tracers), diffusivities=nothing)
    fill_halo_regions!(merge(velocities, tracers), arch, nothing, state)

    U, C = datatuples(velocities, tracers)

    return (   ∇_κ_∇c(2, 1, 3, grid, closure, C.T, Val(1)) == 8κh + 10κv &&
            ∂ⱼ_2ν_Σ₁ⱼ(2, 1, 3, grid, closure, U) == 2νh + 4νv &&
            ∂ⱼ_2ν_Σ₂ⱼ(2, 1, 3, grid, closure, U) == 4νh + 6νv &&
            ∂ⱼ_2ν_Σ₃ⱼ(2, 1, 3, grid, closure, U) == 6νh + 8νv)
end

function test_calculate_diffusivities(arch, closurename, FT=Float64; kwargs...)
      tracernames = (:b,)
          closure = getproperty(TurbulenceClosures, closurename)(FT, kwargs...)
          closure = with_tracers(tracernames, closure)
             grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(3, 3, 3))
    diffusivities = DiffusivityFields(arch, grid, tracernames, closure)
         buoyancy = BuoyancyTracer()
       velocities = VelocityFields(arch, grid)
          tracers = TracerFields(arch, grid, tracernames)

    U, C, K = datatuples(velocities, tracers, diffusivities)
    calculate_diffusivities!(K, arch, grid, closure, buoyancy, U, C)

    return true
end

function time_step_with_tupled_closure(FT, arch)
    closure_tuple = (AnisotropicMinimumDissipation(FT), ConstantAnisotropicDiffusivity(FT))

    model = IncompressibleModel(
        architecture=arch, float_type=FT, closure=closure_tuple,
        grid=RegularCartesianGrid(FT, size=(16, 16, 16), extent=(1, 2, 3))
    )

    time_step!(model, 1, euler=true)
    return true
end

function compute_closure_specific_diffusive_cfl(closurename)
    grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 2, 3))
    closure = getproperty(TurbulenceClosures, closurename)()

    model = IncompressibleModel(grid=grid, closure=closure)
    dcfl = DiffusiveCFL(0.1)
    @test dcfl(model) isa Number

    tracerless_model = IncompressibleModel(grid=grid, closure=closure,
                                           buoyancy=nothing, tracers=nothing)
    dcfl = DiffusiveCFL(0.2)
    @test dcfl(tracerless_model) isa Number

    return nothing
end

@testset "Turbulence closures" begin
    @info "Testing turbulence closures..."

    @testset "Closure instantiation" begin
        @info "  Testing closure instantiation..."
        for T in float_types
            for closure in closures
                @test test_closure_instantiation(T, closure)
            end
        end
    end

    @testset "Constant isotropic diffusivity" begin
        @info "  Testing constant isotropic diffusivity..."
        for T in float_types
            @test test_constant_isotropic_diffusivity_basic(T)
            @test test_constant_isotropic_diffusivity_fluxdiv(T)
        end
    end

    @testset "Constant anisotropic diffusivity" begin
        @info "  Testing constant anisotropic diffusivity..."
        for T in float_types
            @test test_anisotropic_diffusivity_fluxdiv(T, νv=zero(T), νh=zero(T))
            @test test_anisotropic_diffusivity_fluxdiv(T)
        end
    end

    @testset "Calculation of nonlinear diffusivities" begin
        @info "  Testing calculation of nonlinear diffusivities..."
        for FT in [Float64]
            for arch in archs
                for closure in closures
                    @info "    Calculating diffusivities for $closure [$FT, $(typeof(arch))]"
                    @test test_calculate_diffusivities(arch, closure, FT)
                end
            end
        end
    end

    @testset "Closure tuples" begin
        @info "  Testing time-stepping with a tuple of closures..."
        for arch in archs
            for FT in float_types
                @test time_step_with_tupled_closure(FT, arch)
            end
        end
    end

    @testset "Diagnostics" begin
        @info "  Testing turbulence closure diagnostics..."
        for closure in closures
            compute_closure_specific_diffusive_cfl(closure)
        end
    end
end


text:
function adjust_step(tbar :: Float64, t :: Float64, dt :: Float64)

    lstat  = 0

    if (abs(tbar - (t + dt)) / tbar < 1.e-15)

        lstat = 1

    elseif (t + dt > tbar)

        dt    = tbar - t
        lstat = 1

    end

    dt, lstat

end


text:
"""
    blocks(a::AbstractArray{T,N}) :: AbstractArray{<:AbstractArray{T,N},N}

Return the array-of-arrays view to `a` such that

```
blocks(a)[i₁, i₂, ..., iₙ] == a[Block(i₁), Block(i₂), ..., Block(iₙ)]
```

This function does not copy the blocks and give a mutable viwe to the original
array.  This is an "inverse" of [`mortar`](@ref).

# Examples
```jldoctest; setup = quote using BlockArrays end
julia> bs1 = permutedims(reshape([
               1ones(1, 3), 2ones(1, 2),
               3ones(2, 3), 4ones(2, 2),
           ], (2, 2)))
2×2 Matrix{Matrix{Float64}}:
 [1.0 1.0 1.0]               [2.0 2.0]
 [3.0 3.0 3.0; 3.0 3.0 3.0]  [4.0 4.0; 4.0 4.0]

julia> a = mortar(bs1)
2×2-blocked 3×5 BlockMatrix{Float64}:
 1.0  1.0  1.0  │  2.0  2.0
 ───────────────┼──────────
 3.0  3.0  3.0  │  4.0  4.0
 3.0  3.0  3.0  │  4.0  4.0

julia> bs2 = blocks(a)
2×2 Matrix{Matrix{Float64}}:
 [1.0 1.0 1.0]               [2.0 2.0]
 [3.0 3.0 3.0; 3.0 3.0 3.0]  [4.0 4.0; 4.0 4.0]

julia> bs1 == bs2
true

julia> bs2[1, 1] .*= 100;

julia> a  # in-place mutation is reflected to the block array
2×2-blocked 3×5 BlockMatrix{Float64}:
 100.0  100.0  100.0  │  2.0  2.0
 ─────────────────────┼──────────
   3.0    3.0    3.0  │  4.0  4.0
   3.0    3.0    3.0  │  4.0  4.0
```
"""
blocks(a::AbstractArray) = BlocksView(a)
blocks(a::BlockArray) = a.blocks
blocks(A::Adjoint) = adjoint(blocks(parent(A)))
blocks(A::Transpose) = transpose(blocks(parent(A)))

# convert a tuple of BlockRange to a tuple of `AbstractUnitRange{Int}`
_block2int(B::Block{1}) = Int(B):Int(B)
_block2int(B::BlockRange{1}) = Int.(B)
_blockrange2int() = ()
_blockrange2int(A, B...) = tuple(_block2int(A.block), _blockrange2int(B...)...)

blocks(A::SubArray{<:Any,N,<:Any,<:NTuple{N,BlockSlice}}) where N =
    view(blocks(parent(A)), _blockrange2int(parentindices(A)...)...)

struct BlocksView{
    S,                            # eltype(eltype(BlocksView(...)))
    N,                            # ndims
    T<:AbstractArray{S,N},        # eltype(BlocksView(...)), i.e., block type
    B<:AbstractArray{S,N},   # array to be wrapped
} <: AbstractArray{T,N}
    array::B
end

BlocksView(a::AbstractArray{S,N}) where {S,N} =
    BlocksView{S,N,AbstractArray{eltype(a),N},typeof(a)}(a)
# Note: deciding concrete eltype of `BlocksView` requires some extra
# interface for `AbstractBlockArray`.

Base.IteratorEltype(::Type{<:BlocksView}) = Base.EltypeUnknown()

Base.size(a::BlocksView) = blocksize(a.array)
Base.axes(a::BlocksView) = map(br -> only(br.indices), blockaxes(a.array))

#=
This is broken for now. See: https://github.com/JuliaArrays/BlockArrays.jl/issues/120
# IndexLinear implementations
@propagate_inbounds Base.getindex(a::BlocksView, i::Int) = view(a.array, Block(i))
@propagate_inbounds Base.setindex!(a::BlocksView, b, i::Int) = copyto!(a[i], b)
=#

# IndexCartesian implementations
@propagate_inbounds Base.getindex(a::BlocksView{T,N}, i::Vararg{Int,N}) where {T,N} =
    view(a.array, Block.(i)...)
@propagate_inbounds Base.setindex!(a::BlocksView{T,N}, b, i::Vararg{Int,N}) where {T,N} =
    copyto!(a[i...], b)

function Base.showarg(io::IO, a::BlocksView, toplevel::Bool)
    if toplevel
        print(io, "blocks of ")
        Base.showarg(io, a.array, true)
    else
        print(io, "::BlocksView{…,")
        Base.showarg(io, a.array, false)
        print(io, '}')
    end
end


text:
function get_solution_acdc(pm::AbstractPowerModel, sol::Dict{String,Any})
    PowerModels.add_setpoint_bus_voltage!(sol, pm)
    PowerModels.add_setpoint_generator_power!(sol, pm)
    PowerModels.add_setpoint_branch_flow!(sol, pm)
    PowerModels.add_setpoint_dcline_flow!(sol, pm)
    if haskey(pm.data, "convdc")
        add_dc_bus_voltage_setpoint(sol, pm)
        add_dcconverter_setpoint(sol, pm)
        add_dcgrid_flow_setpoint(sol,pm)
        add_dcbranch_losses(sol, pm)
        add_dcconverter_losses(sol, pm)
        add_dcconverter_voltage_setpoint(sol, pm)
        add_dcconverter_firing_angle(sol, pm)
    end
    return sol
end

function add_dcconverter_setpoint(sol, pm::AbstractPowerModel)
    PowerModels.add_setpoint!(sol, pm, "convdc", "pgrid", :pconv_tf_fr)
    PowerModels.add_setpoint!(sol, pm, "convdc", "qgrid", :qconv_tf_fr)
    PowerModels.add_setpoint!(sol, pm, "convdc", "pconv", :pconv_ac)
    PowerModels.add_setpoint!(sol, pm, "convdc", "qconv", :qconv_ac)
    PowerModels.add_setpoint!(sol, pm, "convdc", "pdc", :pconv_dc)
    PowerModels.add_setpoint!(sol, pm, "convdc", "iconv", :iconv_ac)
    PowerModels.add_setpoint!(sol, pm, "convdc", "ptf_fr", :pconv_tf_fr)
    PowerModels.add_setpoint!(sol, pm, "convdc", "ptf_to", :pconv_tf_to)
    add_dcconverter_voltage_setpoint(sol, pm)
end

function add_dcgrid_flow_setpoint(sol, pm::AbstractPowerModel)
    # check the branch flows were requested
    PowerModels.add_setpoint!(sol, pm, "branchdc", "pf", :p_dcgrid,  status_name="status",  var_key = (idx,item) -> (idx, item["fbusdc"], item["tbusdc"]))
    PowerModels.add_setpoint!(sol, pm, "branchdc", "pt", :p_dcgrid;  status_name="status",  var_key = (idx,item) -> (idx, item["tbusdc"], item["fbusdc"]))
end

function add_dcbranch_losses(sol, pm::AbstractPowerModel)
    for (i, branchdc) in sol["branchdc"]
        pf = branchdc["pf"]
        pt = branchdc["pt"]
        sol["branchdc"]["$i"]["ploss"] = pf + pt
    end
end

function add_dcconverter_losses(sol, pm::AbstractPowerModel)
    for (i, convdc) in sol["convdc"]
        pf = convdc["pdc"]
        pt = convdc["pgrid"]
        ptf_fr = convdc["ptf_fr"]
        ptf_to = convdc["ptf_to"]
        sol["convdc"]["$i"]["ploss_tot"] = pf + pt
        sol["convdc"]["$i"]["ploss_tf"] = ptf_fr + ptf_to
        sol["convdc"]["$i"]["ploss_conv"] = sol["convdc"]["$i"]["ploss_tot"] - sol["convdc"]["$i"]["ploss_tf"] # TODO update later with filter etc..
    end
end

function add_dc_bus_voltage_setpoint(sol, pm::AbstractPowerModel)
    PowerModels.add_setpoint!(sol, pm, "busdc", "vm", :vdcm, status_name="Vdc", inactive_status_value = 4)
end

function add_dcconverter_voltage_setpoint(sol, pm::AbstractPowerModel)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vmconv", :vmc, status_name="islcc", inactive_status_value = 0)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vaconv", :vac, status_name="islcc", inactive_status_value = 0)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vmfilt", :vmf, status_name="islcc", inactive_status_value = 0)
    PowerModels.add_setpoint!(sol, pm, "convdc", "vafilt", :vaf, status_name="islcc", inactive_status_value = 0)
end

function add_dcconverter_firing_angle(sol, pm::AbstractPowerModel)
    PowerModels.add_setpoint!(sol, pm, "convdc", "phi", :phiconv, status_name="islcc", inactive_status_value = 0)
end


text:
const BlockSparseMatrix{ElT,StoreT,IndsT} = BlockSparseTensor{ElT,2,StoreT,IndsT}
const DiagBlockSparseMatrix{ElT,StoreT,IndsT} = DiagBlockSparseTensor{ElT,2,StoreT,IndsT}
const DiagMatrix{ElT,StoreT,IndsT} = DiagTensor{ElT,2,StoreT,IndsT}

function _truncated_blockdim(
  S::DiagMatrix, docut::Float64; singular_values=false, truncate=true
)
  !truncate && return diaglength(S)
  newdim = 0
  val = singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1))
  while newdim + 1 ≤ diaglength(S) && val > docut
    newdim += 1
    if newdim + 1 ≤ diaglength(S)
      val =
        singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1))
    end
  end
  return newdim
end

"""
    svd(T::BlockSparseTensor{<:Number,2}; kwargs...)

svd of an order-2 BlockSparseTensor.

This function assumes that there is one block
per row/column, otherwise it fails.
This assumption makes it so the result can be
computed from the dense svds of seperate blocks.
"""
function LinearAlgebra.svd(T::BlockSparseMatrix{ElT}; kwargs...) where {ElT}
  alg::String = get(kwargs, :alg, "divide_and_conquer")

  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)

  #@timeit_debug timer "block sparse svd" begin
  Us = Vector{DenseTensor{ElT,2}}(undef, nnzblocks(T))
  Ss = Vector{DiagTensor{real(ElT),2}}(undef, nnzblocks(T))
  Vs = Vector{DenseTensor{ElT,2}}(undef, nnzblocks(T))

  # Sorted eigenvalues
  d = Vector{real(ElT)}()

  for (n, b) in enumerate(eachnzblock(T))
    blockT = blockview(T, b)
    USVb = svd(blockT; alg=alg)
    if isnothing(USVb)
      return nothing
    end
    Ub, Sb, Vb = USVb
    Us[n] = Ub
    Ss[n] = Sb
    Vs[n] = Vb
    # Previously this was:
    # vector(diag(Sb))
    # But it broke, did `diag(::Tensor)` change types?
    # TODO: call this a function `diagonal`, i.e.:
    # https://github.com/JuliaLang/julia/issues/30250
    # or make `diag(::Tensor)` return a view by default.
    append!(d, data(Sb))
  end

  # Square the singular values to get
  # the eigenvalues
  d .= d .^ 2
  sort!(d; rev=true)

  # Get the list of blocks of T
  # that are not dropped
  nzblocksT = nzblocks(T)

  dropblocks = Int[]
  if truncate
    truncerr, docut = truncate!(d; kwargs...)
    for n in 1:nnzblocks(T)
      blockdim = _truncated_blockdim(Ss[n], docut; singular_values=true, truncate=truncate)
      if blockdim == 0
        push!(dropblocks, n)
      else
        Strunc = tensor(Diag(storage(Ss[n])[1:blockdim]), (blockdim, blockdim))
        Us[n] = Us[n][1:dim(Us[n], 1), 1:blockdim]
        Ss[n] = Strunc
        Vs[n] = Vs[n][1:dim(Vs[n], 1), 1:blockdim]
      end
    end
    deleteat!(Us, dropblocks)
    deleteat!(Ss, dropblocks)
    deleteat!(Vs, dropblocks)
    deleteat!(nzblocksT, dropblocks)
  else
    truncerr, docut = 0.0, 0.0
  end

  # The number of non-zero blocks of T remaining
  nnzblocksT = length(nzblocksT)

  #
  # Make indices of U and V
  # that connect to S
  #
  i1 = ind(T, 1)
  i2 = ind(T, 2)
  uind = dag(sim(i1))
  vind = dag(sim(i2))
  resize!(uind, nnzblocksT)
  resize!(vind, nnzblocksT)
  for (n, blockT) in enumerate(nzblocksT)
    Udim = size(Us[n], 2)
    b1 = block(i1, blockT[1])
    setblock!(uind, resize(b1, Udim), n)
    Vdim = size(Vs[n], 2)
    b2 = block(i2, blockT[2])
    setblock!(vind, resize(b2, Vdim), n)
  end

  #
  # Put the blocks into U,S,V
  #

  nzblocksU = Vector{Block{2}}(undef, nnzblocksT)
  nzblocksS = Vector{Block{2}}(undef, nnzblocksT)
  nzblocksV = Vector{Block{2}}(undef, nnzblocksT)

  for (n, blockT) in enumerate(nzblocksT)
    blockU = (blockT[1], UInt(n))
    nzblocksU[n] = blockU

    blockS = (n, n)
    nzblocksS[n] = blockS

    blockV = (blockT[2], UInt(n))
    nzblocksV[n] = blockV
  end

  indsU = setindex(inds(T), uind, 2)

  indsV = setindex(inds(T), vind, 1)
  indsV = permute(indsV, (2, 1))

  indsS = setindex(inds(T), dag(uind), 1)
  indsS = setindex(indsS, dag(vind), 2)

  U = BlockSparseTensor(ElT, undef, nzblocksU, indsU)
  S = DiagBlockSparseTensor(real(ElT), undef, nzblocksS, indsS)
  V = BlockSparseTensor(ElT, undef, nzblocksV, indsV)

  for n in 1:nnzblocksT
    Ub, Sb, Vb = Us[n], Ss[n], Vs[n]

    blockU = nzblocksU[n]
    blockS = nzblocksS[n]
    blockV = nzblocksV[n]

    if VERSION < v"1.5"
      # In v1.3 and v1.4 of Julia, Ub has
      # a very complicated view wrapper that
      # can't be handled efficiently
      Ub = copy(Ub)
      Vb = copy(Vb)
    end

    blockview(U, blockU) .= Ub
    blockviewS = blockview(S, blockS)
    for i in 1:diaglength(Sb)
      setdiagindex!(blockviewS, getdiagindex(Sb, i), i)
    end

    blockview(V, blockV) .= Vb
  end

  return U, S, V, Spectrum(d, truncerr)
  #end # @timeit_debug
end

_eigen_eltypes(T::Hermitian{ElT,<:BlockSparseMatrix{ElT}}) where {ElT} = real(ElT), ElT

_eigen_eltypes(T::BlockSparseMatrix{ElT}) where {ElT} = complex(ElT), complex(ElT)

function LinearAlgebra.eigen(
  T::Union{Hermitian{ElT,<:BlockSparseMatrix{ElT}},BlockSparseMatrix{ElT}}; kwargs...
) where {ElT<:Union{Real,Complex}}
  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)

  ElD, ElV = _eigen_eltypes(T)

  # Sorted eigenvalues
  d = Vector{real(ElT)}()

  for b in eachnzblock(T)
    all(==(b[1]), b) || error("Eigen currently only supports block diagonal matrices.")
  end

  b = first(eachnzblock(T))
  blockT = blockview(T, b)
  Db, Vb = eigen(blockT)
  Ds = [Db]
  Vs = [Vb]
  append!(d, abs.(data(Db)))
  for (n, b) in enumerate(eachnzblock(T))
    n == 1 && continue
    blockT = blockview(T, b)
    Db, Vb = eigen(blockT)
    push!(Ds, Db)
    push!(Vs, Vb)
    append!(d, abs.(data(Db)))
  end

  dropblocks = Int[]
  sort!(d; rev=true, by=abs)

  if truncate
    truncerr, docut = truncate!(d; kwargs...)
    for n in 1:nnzblocks(T)
      blockdim = _truncated_blockdim(Ds[n], docut)
      if blockdim == 0
        push!(dropblocks, n)
      else
        Dtrunc = tensor(Diag(storage(Ds[n])[1:blockdim]), (blockdim, blockdim))
        Ds[n] = Dtrunc
        Vs[n] = copy(Vs[n][1:dim(Vs[n], 1), 1:blockdim])
      end
    end
    deleteat!(Ds, dropblocks)
    deleteat!(Vs, dropblocks)
  else
    truncerr = 0.0
  end

  # Get the list of blocks of T
  # that are not dropped
  nzblocksT = nzblocks(T)
  deleteat!(nzblocksT, dropblocks)

  # The number of blocks of T remaining
  nnzblocksT = nnzblocks(T) - length(dropblocks)

  #
  # Put the blocks into D, V
  #

  i1, i2 = inds(T)
  l = sim(i1)

  lkeepblocks = Int[bT[1] for bT in nzblocksT]
  ldropblocks = setdiff(1:nblocks(l), lkeepblocks)
  deleteat!(l, ldropblocks)

  # l may have too many blocks
  (nblocks(l) > nnzblocksT) && error("New index l in eigen has too many blocks")

  # Truncation may have changed
  # some block sizes
  for n in 1:nnzblocksT
    setblockdim!(l, minimum(dims(Ds[n])), n)
  end

  r = dag(sim(l))

  indsD = (l, r)
  indsV = (dag(i2), r)

  nzblocksD = Vector{Block{2}}(undef, nnzblocksT)
  nzblocksV = Vector{Block{2}}(undef, nnzblocksT)
  for n in 1:nnzblocksT
    blockT = nzblocksT[n]

    blockD = (n, n)
    nzblocksD[n] = blockD

    blockV = (blockT[1], n)
    nzblocksV[n] = blockV
  end

  D = DiagBlockSparseTensor(ElD, undef, nzblocksD, indsD)
  V = BlockSparseTensor(ElV, undef, nzblocksV, indsV)

  for n in 1:nnzblocksT
    Db, Vb = Ds[n], Vs[n]

    blockD = nzblocksD[n]
    blockviewD = blockview(D, blockD)
    for i in 1:diaglength(Db)
      setdiagindex!(blockviewD, getdiagindex(Db, i), i)
    end

    blockV = nzblocksV[n]
    blockview(V, blockV) .= Vb
  end

  return D, V, Spectrum(d, truncerr)
end

function exp(
  T::Union{BlockSparseMatrix{ElT},Hermitian{ElT,<:BlockSparseMatrix{ElT}}}
) where {ElT<:Union{Real,Complex}}
  expT = BlockSparseTensor(ElT, undef, nzblocks(T), inds(T))
  for b in eachnzblock(T)
    all(==(b[1]), b) || error("exp currently supports only block-diagonal matrices")
  end
  for b in eachdiagblock(T)
    blockT = blockview(T, b)
    if isnothing(blockT)
      # Block was not found in the list, treat as 0
      id_block = Matrix{ElT}(I, blockdims(T, b))
      insertblock!(expT, b)
      blockview(expT, b) .= id_block
    else
      blockview(expT, b) .= exp(blockT)
    end
  end
  return expT
end


text:
using CSDP
const CSDP_PRIMAL_OPT = CSDP.Optimizer(printlevel = 0)
const CSDP_DUAL_OPT = DualOptimizer(CSDP.Optimizer(printlevel = 0))
const CSDP_PRIMAL_FACTORY = with_optimizer(CSDP.Optimizer, printlevel = 0)
const CSDP_DUAL_FACTORY = with_optimizer(DualOptimizer, CSDP.Optimizer(printlevel = 0))

push!(primal_conic_factory, CSDP_PRIMAL_FACTORY)
push!(dual_conic_factory, CSDP_DUAL_FACTORY)
push!(dual_conic_optimizer, CSDP_DUAL_OPT)
push!(primal_conic_optimizer, CSDP_PRIMAL_OPT)

@testset "CSDP SDP triangle Problems" begin
    list_of_sdp_triang_problems = [
        # sdpt1_test, # CSDP is returning SLOW_PROGRESS
        # sdpt2_test, # CSDP is returning SLOW_PROGRESS
        sdpt3_test,
        sdpt4_test
    ]
    test_strong_duality(list_of_sdp_triang_problems, CSDP_PRIMAL_FACTORY)
end


text:
import ADIOS2

function hello_adios2(data::AbstractString)
    adios = ADIOS2.adios_init_mpi(MPI.COMM_WORLD)
    io = ADIOS2.declare_io(adios, "IO")
    open(io, "hello.bp", ADIOS2.mode_write) do engine
        v = ADIOS2.define_variable(io, "hello", data)
        ADIOS2.put!(engine, v, data)
    end
end