text:
using LinearAlgebra
using SparseArrays

# Calculate the objective of logistic regression:
# \frac{1}{n} \sum_{i=1}^n \log( 1 + \exp( -w^T x_i y_i ) ) + \lambda/2 \|w\|^2.
function obj(
    X::SparseMatrixCSC{Float64, Int64},
    y::Vector{Int64},
    W::Matrix{Float64},
    λ::Float64
)
    n, _ = size(X)
    XW = X*W
    objval = 0.0
    for i = 1:n
        prob = softmax(XW[i,:])
        objval += -log(prob[ y[i] ])
    end
    objval /= n
    objval += λ/2*norm(W)^2
    return objval
end

# Calculate the objective of neural network model:
function obj(
    Xt::SparseMatrixCSC{Float64, Int64},
    Y::Flux.OneHotArray,
    W::Flux.Chain,
    λ::Float64
)
    loss(x,y) = Flux.Losses.crossentropy(W(x), y)
    l = 0.0
    num_data = size(Xt, 2)
    for i = 1:num_data
        l += loss(Xt[:,i], Y[:,i])
    end
    sqnorm(w) = sum(abs2, w)
    return l/num_data + (λ/2) * sum(sqnorm, params(W))
end

# Line-search for softmax
function lineSearch(
    X::SparseMatrixCSC{Float64, Int64},
    y::Vector{Int64},
    D::Matrix{Float64},
    W::Matrix{Float64},
    G::Matrix{Float64},
    λ::Float64
)
    maxLineSearchSteps = 100
    GTD = dot(G, D)
    # @show(GTD)
    Wnorm = norm(W)
    XW = X*W
    XD = X*D

    n = size(X, 1)
    objval = 0.0
    @inbounds for i = 1:n
        prob = softmax(XW[i,:])
        objval += -log(prob[ y[i] ])
    end
    objval /= n
    objval += λ/2*Wnorm^2

    Wnew = zeros(Float64, size(W))
    XWnew = zeros(Float64, size(XW))

    η = 1.0
    β = 1e-2
    for t = 1:maxLineSearchSteps
        Wnew .= W - η*D
        XWnew .= XW - η*XD
        objNew = 0.0
        @inbounds for i = 1:n
            prob = softmax(XWnew[i,:])
            objNew += -log(prob[ y[i] ])
        end
        objNew /= n
        objNew += λ/2*norm(Wnew)^2
        if objNew > objval - β*η*GTD
            η *= 0.5
        else
            break
        end
        if t == maxLineSearchSteps
            @warn("Reached maximum linesearch steps.")
        end
    end
    return η
end

# Line-search for FedDCD
function lineSearch2(
    X::SparseMatrixCSC{Float64, Int64},
    Y::Vector{Int64},
    y::Matrix{Float64},  # dual variable
    D::Matrix{Float64},
    W::Matrix{Float64},
    G::Matrix{Float64},
    λ::Float64
)
    maxLineSearchSteps = 100
    GTD = dot(G, D)
    # @show(GTD)
    Wnorm = norm(W)
    XW = X*W
    XD = X*D
    WTy = dot(W, y)
    DTy = dot(D, y)

    n = size(X, 1)
    objval = 0.0
    @inbounds for i = 1:n
        prob = softmax(XW[i,:])
        objval += -log(prob[ Y[i] ])
    end
    objval /= n
    objval -= WTy   # objective have an additional term -W^Ty
    objval += λ/2*Wnorm^2

    Wnew = zeros(Float64, size(W))
    XWnew = zeros(Float64, size(XW))

    η = 1.0
    β = 1e-2
    for t = 1:maxLineSearchSteps
        Wnew .= W - η*D
        XWnew .= XW - η*XD
        objNew = 0.0
        @inbounds for i = 1:n
            prob = softmax(XWnew[i,:])
            objNew += -log(prob[ Y[i] ])
        end
        objNew /= n
        objNew -= ( WTy - η*DTy )
        objNew += λ/2*norm(Wnew)^2
        if objNew > objval - β*η*GTD
            η *= 0.5
        else
            break
        end
        if t == maxLineSearchSteps
            @warn("Reached maximum linesearch steps.")
        end
    end
    return η
end

# Calculate the accuracy
function accuracy(
    X::SparseMatrixCSC{Float64, Int64},
    y::Vector{Int64},
    W::Matrix{Float64},
)
    n = size(X, 1)
    XW = X*W
    ret = 0
    for i = 1:n
        prob = softmax( XW[i,:] )
        if argmax(prob) == y[i]
            ret += 1
        end
    end
    return ret/n
end

# Calculate the accuracy
function accuracy(
    X::SparseMatrixCSC{Float64, Int64},
    Y::Flux.OneHotArray,
    W::Flux.Chain,
)
    Xt = copy(X')
    acc(x,y) = 1.0* ( Flux.onecold(W(x)) == Flux.onecold(y) )
    num_data = size(Xt, 2)
    out = 0.0
    for i = 1:num_data
        out += acc(Xt[:,i], Y[:,i])
    end
    return out/num_data
end

# Calculate the stochastic gradient
function getStochasticGrad(
    Xt::SparseMatrixCSC{Float64, Int64},
    y::Vector{Int64},
    W::Matrix{Float64},
    idx::Int64
)
    x = Xt[:,idx]
    Wtx = W'*x
    s = softmax( Wtx )
    s[ y[idx] ] -= 1
    g = x * s'
    return g
end

# Calculate gradient
function getGradient(
    X::SparseMatrixCSC{Float64, Int64},
    Xt::SparseMatrixCSC{Float64, Int64},
    y::Vector{Int64},
    W::Matrix{Float64},
    λ::Float64
)
    n, d = size(X)
    _, K = size(W)
    g = zeros(Float64, d, K)
    XW = X*W
    @inbounds for i = 1:n
        s = softmax( XW[i,:] )
        s[ y[i] ] -= 1
        x = Xt[:,i]
        I, V = findnz(x)
        @inbounds for col = 1:K
            @inbounds for j = 1:length(I)
                row = I[j]
                val = V[j]
                g[row, col] += val*s[col]
            end
        end
    end
    g ./= n
    g += λ.*W
    return g
end

# Add a dense matrix with a sparse matrix in place, A = A + x * b', x sparse, b dense.
function FastHvMatrixUpdate!(
    A::Matrix{Float64},
    x::SparseVector{Float64,Int64},
    b::Vector{Float64}
)
    I, V = findnz(x)
    K = length(b)
    @inbounds for col = 1:K
        @inbounds for idx = 1:length(I)
            row = I[idx]
            v = V[idx]
            A[row, col] += v*b[col]
        end
    end
    return nothing
end

# Fast Hessian-vector product for softmax classification
function FastHv(
    X::SparseMatrixCSC{Float64, Int64},
    Xt::SparseMatrixCSC{Float64, Int64},
    W::Matrix{Float64},
    λ::Float64,
    V::Matrix{Float64}
)
    # Time complexity is O( nnz(X)*K )
    n, d = size(X)
    _, K = size(W)
    XW = X*W
    XV = X*V
    P = zeros(Float64, n, K)
    ret = zeros(Float64, d, K)
    for i = 1:n
        P[i,:] = softmax( XW[i,:] )
    end
    @inbounds for i = 1:n
        # Complexity of the loop: nnz(xi)*K
        # Calculate (x_i x_i^T V ( Λ - p_i p_i^T ) )
        xi = Xt[:, i]
        ppi = P[i, :]
        # First calculate x_i^T V
        xtV = XV[i, :]
        # # Calculate x_i^T V ( Λ - p_i p_i^T )
        xtVppt = xtV .* ppi
        xtVppt -= dot(xtV, ppi) .* ppi
        # # Calculate x_i x_i^T V ( Λ - p_i p_i^T )
        FastHvMatrixUpdate!( ret, xi, xtVppt)
        # ret += xi*xtVppt'
    end
    ret ./= n
    ret += λ.*V
    return ret
end

# Compute Newton direction. Use cg to solve the linear system H D = g.
function ComputeNewtonDirection(
    X::SparseMatrixCSC{Float64, Int64},
    Xt::SparseMatrixCSC{Float64, Int64},
    W::Matrix{Float64},
    λ::Float64,
    g::Matrix{Float64}
)
    n, d = size(X)
    _, K = size(W)

    H = LinearMap(v->vec(FastHv(X, Xt, W, λ, reshape(v,d,K))), d*K, issymmetric=true, isposdef=true)
    D = cg(H, vec(g), abstol=1e-8, reltol=1e-4, maxiter=1000)

    return reshape(D, d, K)
end

# Using Newton's method to solve the softmax classification problem.
function SoftmaxNewtonMethod(
    X::SparseMatrixCSC{Float64, Int64},
    Xt::SparseMatrixCSC{Float64, Int64},
    y::Vector{Int64},
    W::Matrix{Float64},
    λ::Float64,
    maxIter::Int64=20,
    tol::Float64=1e-4
)
    startTime = time()
    for iter = 1:maxIter
        objval = obj(X, y, W, λ)
        g = getGradient(X, Xt, y, W, λ)
        gnorm = norm(g)
        @printf("Iter %3d, obj: %4.5e, gnorm: %4.5e, time: %4.2f\n", iter, objval, gnorm, time()-startTime)
        if gnorm < tol
            break
        end
        # Compute Newton direction.
        D = ComputeNewtonDirection( X, Xt, W, λ, g)
        η = lineSearch(X, Y, D, W, g, λ)
        W .-= η*D
    end
    return W
end


text:
@enum ENTITY_TYPES begin
    is_function = 1
    is_type = 2
    is_abstract_type = 3
    is_untyped = 4
    not_found = 5
end

struct Detective
    moduleinst::Module
    modulename::Symbol
    allnames::Vector{Symbol}
    functions::Vector{Symbol}
    types::Vector{Symbol}
    abstracttypes::Vector{Symbol}
    undefined_exports::Vector{Symbol}
    graph::SimpleDiGraph
    nv::Int
    tag::Dict{Int,Symbol}
    lookup::Dict{Symbol,Int}
end

"""
    functions(d::Detective)::Vector{Symbol}

Returns a list of functions of a module in a `Detective` object.

"""
functions(d::Detective) = d.functions

"""
    types(d::Detective)::Vector{Symbol}

Returns a list of types of a module in a `Detective` object.

"""
types(d::Detective) = d.types

"""
    abstracttypes(d::Detective)::Vector{Symbol}

Returns a list of abstract types of a module in a `Detective` object.

"""
abstracttypes(d::Detective) = d.abstracttypes

"""
    undefined(d::Detective)::Vector{Symbol}

Returns a list of undefined objects in a module in a `Detective` object.
This could be exported items that don't exist, or objects without explicit types,
until compile time.

"""
undefined(d::Detective) = d.undefined_exports

"""
    inquire( d::Detective, s::Union{Symbol,String} )

This function will search a loaded module for an instance or method of symbol `s`.
If it finds it, it will, return an enumeration of its category. This is useful
for quickly categorizing an item in a struct for display purposes.

"""
function inquire(d::Detective, s::Symbol)
    typeis = nothing
    if s in d.types
        typeis = is_type
    elseif s in d.functions
        typeis = is_function
    elseif s in d.abstracttypes
        typeis = is_abstract_type
    elseif s in d.undefined_exports
        typeis = is_untyped
    else
        typeis = not_found
    end
    return typeis
end
inquire(d::Detective, s::String) = inquire(d, Symbol(s))

function safeisfield(m::Module, s::Symbol, t::Type)
    try
        return isa(getfield(m, s), t)
    catch
        return false
    end
end

function safeisnotabstract(m::Module, s::Symbol, t::Type)
    return safeisfield(m, s, t) && !isabstracttype(getfield(m, s))
end

function safeisabstract(m::Module, s::Symbol, t::Type)
    return safeisfield(m, s, t) && isabstracttype(getfield(m, s))
end

"""
    ismacro(s::Symbol)

checks the first character of an input symbol to see if it contains a `@` or not.

"""
ismacro(s::Symbol) = first(string(s)) == '@'

"""
    Detective( moduleinst::Module )

Instantiate Detective object via a module.

"""
function Detective(moduleinst::Module)
    modname = Symbol(moduleinst)
    allnames = [n for n in names(moduleinst)]
    nv = length(allnames)
    graph = SimpleDiGraph(nv)
    tags = Dict(1:nv .=> allnames)
    lookup = Dict(allnames .=> 1:nv)
    #remove package name from list
    allnames = allnames[allnames .!= modname]
    fns = [safeisfield(moduleinst, curname, Function) for curname in allnames]
    types = [
        safeisnotabstract(moduleinst, curname, Type) && !ismacro(curname) for
        curname in allnames
    ]
    #remove constructors from functions?
    fns[fns .& types] .= false
    abstypes = [
        safeisabstract(moduleinst, curname, Type) && !ismacro(curname) for
        curname in allnames
    ]
    others = (fns .+ types .+ abstypes) .== 0
    return Detective(
        moduleinst,
        modname,
        allnames,
        allnames[fns],
        allnames[types],
        allnames[abstypes],
        allnames[others],
        graph,
        nv,
        tags,
        lookup,
    )
end

ext: jl


text:
#! /usr/bin/julia

# Rosetta Code, Ranking methods

function ties{T<:Real}(a::Array{T,1})
    unique(a[2:end][a[2:end] .== a[1:end-1]])
end

function rankstandard{T<:Real}(a::Array{T,1})
    r = collect(1:length(a))
    1 < r[end] || return r
    for i in ties(a)
        r[a.==i] = r[a.==i][1]
    end
    return r
end

function rankmodified{T<:Real}(a::Array{T,1})
    indexin(a, a)
end

function rankdense{T<:Real}(a::Array{T,1})
    indexin(a, unique(a))
end

function rankordinal{T<:Real}(a::Array{T,1})
    collect(1:length(a))
end

function rankfractional{T<:Real}(a::Array{T,1})
    r = float64(collect(1:length(a)))
    1.0 < r[end] || return r
    for i in ties(a)
        r[a.==i] = mean(r[a.==i])
    end
    return r
end

scores = [44, 42, 42, 41, 41, 41, 39]
names = ["Solomon", "Jason", "Errol", "Garry",
         "Bernard", "Barry", "Stephen"]

srank = rankstandard(scores)
mrank = rankmodified(scores)
drank = rankdense(scores)
orank = rankordinal(scores)
frank = rankfractional(scores)

println("    Name    Score  Std  Mod  Den  Ord  Frac")
for i in 1:length(scores)
    print(@sprintf("   %-7s", names[i]))
    print(@sprintf("%5d ", scores[i]))
    print(@sprintf("%5d", srank[i]))
    print(@sprintf("%5d", mrank[i]))
    print(@sprintf("%5d", drank[i]))
    print(@sprintf("%5d", orank[i]))
    print(@sprintf("%7.2f", frank[i]))
    println()
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
"""
    AshPriors

Empirical Bayes priors that are used in the simulations of:
> Stephens, M., 2017. False discovery rates: a new deal. Biostatistics, 18(2), pp.275-294.
"""
const AshPriors = Dict(
    :Spiky => MixtureModel(
        [
            Normal(0, 0.25)
            Normal(0, 0.5)
            Normal(0, 1)
            Normal(0, 2)
        ],
        [0.4; 0.2; 0.2; 0.2],
    ),
    :NearNormal => MixtureModel([Normal(0, 1); Normal(0, 2)], [1 / 3; 2 / 3]),
    :FlatTop => MixtureModel([Normal(x, 0.5) for x = -1.5:0.5:1.5]),
    :Skew => MixtureModel(
        [Normal(-2, 2), Normal(-1, 1.5), Normal(0, 1), Normal(1, 1)],
        [1 / 4; 1 / 4; 1 / 3; 1 / 6],
    ),
    :BigNormal => Normal(0, 4),
    :Bimodal => MixtureModel([Normal(-2, 1); Normal(2, 1)], [0.5; 0.5]),
)

"""
    MarronWandGaussianMixtures

Flexible Gaussian Mixture distributions described in

> Marron, J. Steve, and Matt P. Wand. Exact mean integrated squared error.
The Annals of Statistics (1992): 712-736.
"""
const MarronWandGaussianMixtures = Dict(
    :KurtoticUnimodal => MixtureModel([Normal(0,1), Normal(0, 0.1)],
                                          [2/3, 1/3]),
    :SkewedUnimodal => MixtureModel([Normal(0,1), Normal(0.5, 2/3), Normal(13/12, 5/9)],
                                        [1/5, 1/5, 3/5]),
    :Bimodal => MixtureModel([Normal(-1.0, 2/3), Normal(1.0, 2/3)]),
    :SeparatedBimodal => MixtureModel([Normal(-1.5, 0.5), Normal(1.5, 0.5)]),
    :SkewedBimodal =>  MixtureModel([Normal(0,1), Normal(3/2, 1/3)],
                                    [0.75, 0.25]),
    :Outlier => MixtureModel([Normal(0,1), Normal(0, 1/10)], [1/10, 9/10])
)

"""
    IWPriors

Empirical Bayes priors that are used in the simulations of:
> Ignatiadis, N. and Wager, S., 2019. Bias-aware confidence intervals for empirical Bayes analysis. arXiv preprint arXiv:1902.02774.
"""
const IWPriors = Dict(
    :Unimodal => MixtureModel([ Normal(-0.2,.2), Normal(0,0.9)],[0.7, 0.3]),
    :Bimodal => MixtureModel([Normal(-1.5,.2), Normal(1.5, .2)]),
    :NegSpiky => MixtureModel([Normal(-0.25, 0.25), Normal(0, 1.0)], [0.8, 0.2]))

"""
    EfronPriors

Empirical Bayes priors that are used in the simulations of:
> Efron, B., 2016. Empirical Bayes deconvolution estimates. Biometrika, 103(1), pp.1-20.
"""
const EfronPriors = Dict(
    :UniformNormal => MixtureModel([Uniform(-3.0,3.0), Normal(0.0, 0.5)], [1/8, 7/8]))


text:
function momsos!(blasso::BLASSO{T, FT, DT},
                 ϵ::Real,
                 maxdeg::Union{Int, String}) where {T,
                                                    PT<:AbstractPolynomialLike,
                                                    FT<:Array{PT, 1},
                                                    DT<:AbstractSemialgebraicSet}

    model = SOSModel(with_optimizer(Mosek.Optimizer, QUIET=true))

    y = blasso.y
    φ = blasso.φ
    m = length(y)

    @variable(model, p[1:m])  # dual variable

    if maxdeg == "defaut"
        maxdeg = max([maxdegree(monomials(φ[i])) for i=1:m]...)
    end

    @constraint(model, cminus,
                1 + sum(φ[i] * p[i] for i=1:m) in SOSCone(),
                domain=blasso.domain, maxdegree=maxdeg)

    @constraint(model, cplus,
                1 - sum(φ[i] * p[i] for i=1:m) in SOSCone(),
                domain=blasso.domain, maxdegree=maxdeg)

    if blasso.λ == 0  # noiseless formulation
        @objective(model, Max, sum(y[i] * p[i] for i=1:m))
    else  # regularized formulation
        # workaround to avoid MosekTools error with quadratic objective
        @variable(model, t)
        @constraint(model, sum(p[i]^2 for i=1:m) <= t)

        @objective(model, Max, sum(y[i] * p[i] for i=1:m) - blasso.λ / 2 * t)
    end

    optimize!(model)
    blasso.objvalue = objective_value(model)
    blasso.dualsol = [value(p[i]) for i=1:m]

    # finds a primal solution given the dual solution
    νplus = moment_matrix(cplus)
    μplus = extractatoms(νplus, ϵ)

    νminus = moment_matrix(cminus)
    μminus = extractatoms(νminus, ϵ)

    # this will raise an error if extraction fails
    # (i.e. μplus or μminus is of type Nothing)
    atoms = μplus.atoms
    append!(atoms, [WeightedDiracMeasure(atom.center, -atom.weight)
                    for atom in μminus.atoms])

    blasso.μ = AtomicMeasure(μplus.variables, atoms)
end


text:
using ISA #References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
using Plots

sig = cos.(0.0:π/100:10pi) +  cos.( 5*(0.0:π/100:10pi))+  cos.( 20*(0.0:π/100:10pi))
φ₁ = SIFT(sig)
plot(sig)
plot!(φ₁)
φ₂ = SIFT(sig-φ₁)
plot!(φ₂)
plot!(sig-φ₁-φ₂)

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

sig = cos.(0.0:π/100:10pi) +  cos.( 5*(0.0:π/100:10pi))+  cos.( 20*(0.0:π/100:10pi))
IMF = EMD(sig)
plot(sig)
plot!(IMF)

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

sig = exp.( 1im*(0.0:π/100:10pi)) +  exp.( 5im*(0.0:π/100:10pi)) +  exp.( 20im*(0.0:π/100:10pi))
φ₁ = ℂSIFT(sig)
plot(real(sig))
plot!(real(φ₁))
φ₂ = ℂSIFT(sig-φ₁)
plot!(real(φ₂))
plot!(real(sig-φ₁-φ₂))

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

sig = exp.( 1im*(0.0:π/100:10pi)) +  exp.( 5im*(0.0:π/100:10pi)) +  exp.( 20im*(0.0:π/100:10pi))
IMF = ℂEMD(sig)
plot(real(sig))
plot!(real.(IMF))

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
using DelimitedFiles
using Plots; gr()

function plot_contour(
                        func_list::Array{Function,1},
                        a::Vector{Float64},
                        b::Vector{Float64},
                        figname::String;
                        p_step::Float64=1.0e-2,
                        p_levels::Int64=30,
                        p_fill::Bool=true
                        )

    # Creates 'x' and 'y'
    x = a[1]:p_step:b[1]
    y = a[2]:p_step:b[2]

    # fmin function
    function fmin(x, y)

        return min( func_list[1](x, y), func_list[2](x, y) )

    end

    plot_problem = contour(x, y, fmin, levels = p_levels, fill = p_fill)
    plot( plot_problem, aspec_ratio = :equal )
    png( figname )

end

# -----------------------------------------------
# Problem definition
# -----------------------------------------------

function f(x, y)

    return - 0.5 * x + 10 * y

end

function g(x, y)

    return x ^ 2.0 + y ^ 2.0

end

f_list = [f, g]
a = [ 0.0, 0.0 ]
b = [ 5.0, 5.0 ]
path_to_file = "./images/test_01/fig_contour"

# -----------------------------------------------
# Function call
# -----------------------------------------------

plot_contour( f_list, a, b, path_to_file )


text:
int main() {
   int a^9;
   return 1;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# Julia wrapper for header: /usr/include/scip/pub_branch.h
# Automatically generated using Clang.jl wrap_c

function SCIPbranchruleComp(elem1, elem2)
    ccall((:SCIPbranchruleComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPbranchruleCompName(elem1, elem2)
    ccall((:SCIPbranchruleCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPbranchruleGetData(branchrule)
    ccall((:SCIPbranchruleGetData, libscip), Ptr{SCIP_BRANCHRULEDATA}, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleSetData(branchrule, branchruledata)
    ccall((:SCIPbranchruleSetData, libscip), Cvoid, (Ptr{SCIP_BRANCHRULE}, Ptr{SCIP_BRANCHRULEDATA}), branchrule, branchruledata)
end

function SCIPbranchruleGetName(branchrule)
    ccall((:SCIPbranchruleGetName, libscip), Cstring, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetDesc(branchrule)
    ccall((:SCIPbranchruleGetDesc, libscip), Cstring, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetPriority(branchrule)
    ccall((:SCIPbranchruleGetPriority, libscip), Cint, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetMaxdepth(branchrule)
    ccall((:SCIPbranchruleGetMaxdepth, libscip), Cint, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetMaxbounddist(branchrule)
    ccall((:SCIPbranchruleGetMaxbounddist, libscip), Cdouble, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetSetupTime(branchrule)
    ccall((:SCIPbranchruleGetSetupTime, libscip), Cdouble, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetTime(branchrule)
    ccall((:SCIPbranchruleGetTime, libscip), Cdouble, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNLPCalls(branchrule)
    ccall((:SCIPbranchruleGetNLPCalls, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNExternCalls(branchrule)
    ccall((:SCIPbranchruleGetNExternCalls, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNPseudoCalls(branchrule)
    ccall((:SCIPbranchruleGetNPseudoCalls, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNCutoffs(branchrule)
    ccall((:SCIPbranchruleGetNCutoffs, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNCutsFound(branchrule)
    ccall((:SCIPbranchruleGetNCutsFound, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNConssFound(branchrule)
    ccall((:SCIPbranchruleGetNConssFound, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNDomredsFound(branchrule)
    ccall((:SCIPbranchruleGetNDomredsFound, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleGetNChildren(branchrule)
    ccall((:SCIPbranchruleGetNChildren, libscip), Clonglong, (Ptr{SCIP_BRANCHRULE},), branchrule)
end

function SCIPbranchruleIsInitialized(branchrule)
    ccall((:SCIPbranchruleIsInitialized, libscip), UInt32, (Ptr{SCIP_BRANCHRULE},), branchrule)
end


text:
using ChaosTools, DynamicalSystemsBase
using Test
using Statistics
using ChaosTools.DynamicalSystemsBase
using ChaosTools.DelayEmbeddings

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

println("\nTesting correlation sum...")
@testset "Correlation dim" begin
    @testset "Henon Map" begin
        ds = Systems.henon()
        ts = trajectory(ds, 5000)
        R = 1.5maximum(maxima(ts) - minima(ts))
        # check for right normalisation with various combinations of `q` and `w`
        @test correlationsum(ts, R; q = 2, w = 0) ≈ 1
        @test correlationsum(ts, R; q = 2, w = 1) ≈ 1
        @test correlationsum(ts, R; q = 2.2, w = 0) ≈ 1
        @test correlationsum(ts, R; q = 2.2, w = 1) ≈ 1
        es = 10 .^ range(-3, stop = 0, length = 7)
        cs1 = correlationsum(ts, es)
        dim1 = linear_region(log.(es), log.(cs1))[2]
        test_value(dim1, 1.1, 1.3)
        cs2 = correlationsum(ts, es; q = 2.0001)
        dim2 = linear_region(log.(es), log.(cs2))[2]
        test_value(dim2, 1.1, 1.3)
        test_value(grassberger_dim(ts), 1.1, 1.3)
    end
    @testset "Lorenz System" begin
        ds = Systems.lorenz()
        ts = trajectory(ds, 1000; Δt = 0.1)
        es = 10 .^ range(-3, stop = 1, length = 8)
        cs1 = correlationsum(ts, es)
        dim1 = linear_region(log.(es), log.(cs1))[2]
        test_value(dim1, 1.85, 2.2)
        cs2 = correlationsum(ts, es; q = 2.001, w = 5)
        dim2 = linear_region(log.(es), log.(cs2))[2]
        test_value(dim2, 1.85, 2.2)
        test_value(grassberger_dim(ts), 1.85, 2.2)
    end
end

println("\nTesting box-assisted correlation sum...")
@testset "Box-assisted correlation sum" begin
    @testset "Henon Map" begin
        ds = Systems.henon()
        ts = standardize(trajectory(ds, 10000))
        r0 = estimate_r0_buenoorovio(ts)[1]
        es = r0 .* 10 .^ range(-2, stop = 0, length = 10)
        C = correlationsum(ts, es)
        @test boxed_correlationsum(ts, es, r0) ≈ C
        @test boxed_correlationsum(ts, es) ≈ C
        @test boxed_correlationsum(ts, es, r0; q = 2.3) ≈ correlationsum(ts, es, q = 2.3)
        @test boxed_correlationsum(ts, es, r0; w = 10) ≈ correlationsum(ts, es; w = 10)
        @test boxed_correlationsum(ts, es, r0; q = 2.3, w = 10) ≈ correlationsum(ts, es; q = 2.3, w = 10)
    end
    @testset "Lorenz System" begin
        ds = Systems.lorenz()
        ts = standardize(trajectory(ds, 1000; Δt = 0.1))
        r0 = estimate_r0_buenoorovio(ts)[1]
        es = r0 .* 10 .^ range(-2, stop = 0, length = 10)
        C = correlationsum(ts, es)
        @test boxed_correlationsum(ts, es, r0) ≈ C
        @test boxed_correlationsum(ts, es) ≈ C
        @test boxed_correlationsum(ts, es, r0; q = 2.3) ≈ correlationsum(ts, es; q = 2.3)
        @test boxed_correlationsum(ts, es, r0; w = 10) ≈ correlationsum(ts, es; w = 10)
        @test boxed_correlationsum(ts, es, r0; q = 2.3, w = 10) ≈ correlationsum(ts, es; q = 2.3, w = 10)
    end
end

println("\nTesting the fixed mass correlation sum...")
@testset "Fixed mass correlation sum" begin
    @testset "Henon Map" begin
        ds = Systems.henon()
        ts = trajectory(ds, 10000)
        rs, ys = correlationsum_fixedmass(ts, 30)
        test_value(linear_region(rs, ys)[2], 1.0, 1.3)
        rs, ys = correlationsum_fixedmass(ts, 30; M = 2000)
        test_value(linear_region(rs, ys)[2], 1.0, 1.3)
        rs, ys = correlationsum_fixedmass(ts, 30; w = 5)
        test_value(linear_region(rs, ys)[2], 1.0, 1.3)
    end
    @testset "Lorenz System" begin
        ds = Systems.lorenz()
        ts = trajectory(ds, 1000; Δt = 0.1)
        rs, ys = correlationsum_fixedmass(ts, 30)
        test_value(linear_region(rs, ys)[2], 1.7, 2.2)
        rs, ys = correlationsum_fixedmass(ts, 30; M = 2000)
        test_value(linear_region(rs, ys)[2], 1.7, 2.2)
        rs, ys = correlationsum_fixedmass(ts, 30; w = 5)
        test_value(linear_region(rs, ys)[2], 1.7, 2.2)
    end
end

println("\nTesting Takens' best estimate")
@testset "Takens best" begin
    @testset "Henon Map" begin
        ds = Systems.henon()
        ts = trajectory(ds, 5000)
        x = ts[:, 1]
        X = embed(x, 2, 1)
        D_C, D_C_95u, D_C_95l = takens_best_estimate(X, std(x)/4)
        @test 1.15 < D_C < 1.25
        @test D_C_95u < 1.05*D_C
        @test D_C_95l > 0.95*D_C
    end
    @testset "Lorenz System" begin
        ds = Systems.lorenz()
        tr = trajectory(ds, 2000; Δt = 0.1)
        x = tr[:, 1]
        τ = estimate_delay(x, "mi_min", 1:20)
        X = embed(x, 4, τ)
        D_C, D_C_95u, D_C_95l = takens_best_estimate(X, std(x)/4)
        @test 1.85 < D_C < 2.1
        @test D_C_95u < 1.05*D_C
        @test D_C_95l > 0.95*D_C
    end
end


text:
using MyPackage2
using Test

@testset "MyPackage2.jl" begin

    x,y=3,4

    #测试正常情况
    @test add_func(x,y) == 7

    #测试异常情况
    x = "blah"
    @test_throws MethodError add_func(x, y)

    #测试错误情况
    #@test subtract_func(x,y) == 0

end


text:
using LinearAlgebra

export BinTree, isleaf, box, crd, ind, sub, leaf, maxdepth

struct BinTree{T<:Real}
    m_isleaf::Bool
    m_box::BBox{T}
    m_crd::Matrix{T}
    m_ind::Tuple{Vector{UInt32},Vector{UInt32}}
    m_sub::Vector{BinTree{T}}

end

# non-trivial constructor
function BinTree(X::Matrix{T},n) where {T<:Real}
    return builderBinTree(X,n)
end

# recursive build of the BinTree
function builderBinTree(X::Matrix{T},n) where {T<:Real}
    if size(X,2) != 3
        throw(DomainError("size(X,2) = $(size(X,2)) != 3. The input array should only have three columns."))
    end
    #
    N = size(X,1)
    isleaf = N < n
    bb     = BBox(X)
    # recursion
    if isleaf
        # if current is a leave, the set of child indices is empty, so is the children set
        return BinTree{T}(isleaf,bb,X,(UInt32[],UInt32[]),BinTree{T}[])
    else isleaf
        # if not a leave, get max edge length for splitting
        dim  = argmax(get_edg(bb))
        iloc = sortperm(view(X,1:N,dim)) # dispatch with same number in each new cell
        ind  = (iloc[1:(N÷2)],iloc[(N÷2+1):N])
        return BinTree{T}(isleaf,bb,X,ind,[
            builderBinTree(X[ind[1],:],n),
            builderBinTree(X[ind[2],:],n)
        ]) # create the children
    end
end

# accessors
isleaf(bt::BinTree) = bt.m_isleaf
box(bt::BinTree)    = bt.m_box
crd(bt::BinTree)    = bt.m_crd
ind(bt::BinTree)    = bt.m_ind
sub(bt::BinTree)    = bt.m_sub

function leaf(bt::BinTree{T}) where {T<:Real}
    v = zeros(T,size(crd(bt),1))
    if length(sub(bt)) == 2
        v[ind(bt)[1]] .= leaf(sub(bt)[1])
        v[ind(bt)[2]] .= leaf(sub(bt)[2])
    else
        v .= rand(T,1)
    end
    return v
end

# get max depth of the BinTree
function maxdepth(bt::BinTree)
    if isleaf(bt)
        return 1
    else
        return max(maxdepth(sub(bt)[1]),maxdepth(sub(bt)[2]))+1
    end
end

max_stars_repo_name: marcbakry/HierarchicalMatrices

max_issues_repo_name: marcbakry/HierarchicalMatrices

max_forks_repo_name: marcbakry/HierarchicalMatrices


text:
export nanfill!

# replace NaNs with the mean
function nanfill!(x::Array{T,1}) where T<: Real
  J = findall(isnan.(x))
  if !isempty(J)
    if length(J) == length(x)
      fill!(x, zero(T))
    else
      x[J] .= T(mean(findall(isnan.(x).==false)))
    end
  end
  return length(J)
end

"""
  nanfill!(S::SeisData)
  nanfill!(C::SeisChannel)

Replace NaNs in `:x` with mean of non-NaN values.
"""
function nanfill!(S::GphysData)
  for i = 1:S.n
    if !isempty(S.x[i])
      nn = nanfill!(S.x[i])
      if nn > 0
        proc_note!(S, i, "nanfill!(S)", "replaced NaNs with the mean of all non-NaN values")
      end
    end
  end
  return nothing
end

function nanfill!(C::GphysChannel)
  nn = nanfill!(C.x)
  if nn > 0
    proc_note!(C, "nanfill!(C)", "replaced NaNs with the mean of all non-NaN values")
  end
  return nothing
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SARSOP_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SARSOP")
JLLWrappers.@generate_main_file("SARSOP", UUID("e478675d-f996-546e-baa1-da59d9de2062"))
end  # module SARSOP_jll


text:
using Test, InMemoryDatasets, Random, CategoricalArrays

@testset "duplicates, unique, unique! with extra argument" begin
    ds1 = Dataset(a = Union{String, Missing}["a", "b", "a", "b", "a", "b"],
                    b = Vector{Union{Int, Missing}}(1:6),
                    c = Union{Int, Missing}[1:3;1:3])
    ds = vcat(ds1, ds1)
    @test findall(duplicates(ds)) == collect(7:12)
    @test findall(duplicates(ds, leave = :last)) == collect(1:6)
    @test findall(duplicates(ds, :)) == collect(7:12)
    @test findall(duplicates(ds, Colon())) == collect(7:12)
    @test findall(duplicates(ds, :a)) == collect(3:12)
    @test findall(duplicates(ds, "a")) == collect(3:12)
    @test findall(duplicates(ds, [:a, :c])) == collect(7:12)
    @test findall(duplicates(ds, ["a", "c"])) == collect(7:12)
    @test findall(duplicates(ds, r"[ac]")) == collect(7:12)
    @test findall(duplicates(ds, Not(2))) == collect(7:12)
    @test findall(duplicates(ds, Not([2]))) == collect(7:12)
    @test findall(duplicates(ds, Not(:b))) == collect(7:12)
    @test findall(duplicates(ds, Not([:b]))) == collect(7:12)
    @test findall(duplicates(ds, Not([false, true, false]))) == collect(7:12)
    @test findall(duplicates(ds, [1, 3])) == collect(7:12)
    @test findall(duplicates(ds, 1)) == collect(3:12)

    @test findall(duplicates(ds, threads = false)) == collect(7:12)
    @test findall(duplicates(ds, leave = :last, threads = false)) == collect(1:6)
    @test findall(duplicates(ds, :, threads = false)) == collect(7:12)
    @test findall(duplicates(ds, Colon(), threads = false)) == collect(7:12)
    @test findall(duplicates(ds, :a, threads = false)) == collect(3:12)
    @test findall(duplicates(ds, "a", threads = false)) == collect(3:12)
    @test findall(duplicates(ds, [:a, :c], threads = false)) == collect(7:12)
    @test findall(duplicates(ds, ["a", "c"], threads = false)) == collect(7:12)
    @test findall(duplicates(ds, r"[ac]", threads = false)) == collect(7:12)
    @test findall(duplicates(ds, Not(2), threads = false)) == collect(7:12)
    @test findall(duplicates(ds, Not([2]), threads = false)) == collect(7:12)
    @test findall(duplicates(ds, Not(:b), threads = false)) == collect(7:12)
    @test findall(duplicates(ds, Not([:b]), threads = false)) == collect(7:12)
    @test findall(duplicates(ds, Not([false, true, false]), threads = false)) == collect(7:12)
    @test findall(duplicates(ds, [1, 3], threads = false)) == collect(7:12)
    @test findall(duplicates(ds, 1, threads = false)) == collect(3:12)

    fmt(x)=1
    setformat!(ds, :a=>fmt)
    @test findall(duplicates(ds, :a, mapformats = true)) == 2:12
    @test findall(duplicates(ds, :a, mapformats = true, threads = false)) == 2:12

    @test unique(ds) == ds1
    @test unique(ds, keep = :last) == ds1
    @test unique(ds, :) == ds1
    @test unique(ds, Colon()) == ds1
    @test unique(ds, 2:3) == ds1
    @test unique(ds, 3) == ds1[1:3, :]
    @test unique(ds, [1, 3]) == ds1
    @test unique(ds, [:a, :c]) == ds1
    @test unique(ds, ["a", "c"]) == ds1
    @test unique(ds, r"[ac]") == ds1
    @test unique(ds, Not(2)) == ds1
    @test unique(ds, Not([2])) == ds1
    @test unique(ds, Not(:b)) == ds1
    @test unique(ds, Not([:b])) == ds1
    @test unique(ds, Not([false, true, false])) == ds1
    @test unique(ds, :a) == ds1[1:2, :]
    @test unique(ds, "a") == ds1[1:2, :]
    @test unique(ds, :a, mapformats = true) == ds[1:1, :]

    @test unique(ds, threads = false) == ds1
    @test unique(ds, keep = :last, threads = false) == ds1
    @test unique(ds, :, threads = false) == ds1
    @test unique(ds, Colon(), threads = false) == ds1
    @test unique(ds, 2:3, threads = false) == ds1
    @test unique(ds, 3, threads = false) == ds1[1:3, :]
    @test unique(ds, [1, 3], threads = false) == ds1
    @test unique(ds, [:a, :c], threads = false) == ds1
    @test unique(ds, ["a", "c"], threads = false) == ds1
    @test unique(ds, r"[ac]", threads = false) == ds1
    @test unique(ds, Not(2), threads = false) == ds1
    @test unique(ds, Not([2]), threads = false) == ds1
    @test unique(ds, Not(:b), threads = false) == ds1
    @test unique(ds, Not([:b]), threads = false) == ds1
    @test unique(ds, Not([false, true, false]), threads = false) == ds1
    @test unique(ds, :a, threads = false) == ds1[1:2, :]
    @test unique(ds, "a", threads = false) == ds1[1:2, :]
    @test unique(ds, :a, mapformats = true, threads = false) == ds[1:1, :]

    @test_throws ArgumentError unique(Dataset())
    @test_throws ArgumentError duplicates(Dataset())

    @test_throws ArgumentError unique(Dataset(), threads = false)
    @test_throws ArgumentError duplicates(Dataset(), threads = false)

    @test unique(copy(ds1), "a") == unique(copy(ds1), :a) == unique(copy(ds1), 1) ==
          ds1[1:2, :]
  @test unique(copy(ds1), "a", threads = false) == unique(copy(ds1), :a, threads = false) == unique(copy(ds1), 1, threads = false) ==
        ds1[1:2, :]
    unique!(ds, [1, 3])
    @test ds == ds1
    for cols in (r"[ac]", Not(:b), Not(2), Not([:b]), Not([2]), Not([false, true, false]))
        ds = vcat(ds1, ds1)
        unique!(ds, cols)
        @test ds == ds1
    end

    for cols in (r"[ac]", Not(:b), Not(2), Not([:b]), Not([2]), Not([false, true, false]))
        ds = vcat(ds1, ds1)
        unique!(ds, cols, threads = false)
        @test ds == ds1
    end

    ds = Dataset([Union{Missing, Int64}[3, 3, 1, 2, 2, 3, 3, 1, 1, 3], Union{Missing, Int64}[1, 1, 1, 1, 2, 1, 2, 2, 1, 3], Union{Missing, Int64}[1, 3, 1, 1, 2, 3, 2, 2, 1, 1]], :auto)
    ds1 = unique(ds)
    ds2 = unique(ds, keep = :last)
    ds3 = unique(ds, keep = :none)
    ds4 = unique(ds, 1, keep = :last)
    ds5 = unique(ds, 1, keep = :none)

    @test ds1 == unique(ds, threads = false)
    @test ds2 == unique(ds, keep = :last, threads = false)
    @test ds3 == unique(ds, keep = :none, threads = false)
    @test ds4 == unique(ds, 1, keep = :last, threads = false)
    @test ds5 == unique(ds, 1, keep = :none, threads = false)

    ds1_t = Dataset([Union{Missing, Int64}[3, 3, 1, 2, 2, 3, 1, 3], Union{Missing, Int64}[1, 1, 1, 1, 2, 2, 2, 3], Union{Missing, Int64}[1, 3, 1, 1, 2, 2, 2, 1]], :auto)
    ds2_t = Dataset([Union{Missing, Int64}[3, 2, 2, 3, 3, 1, 1, 3], Union{Missing, Int64}[1, 1, 2, 1, 2, 2, 1, 3], Union{Missing, Int64}[1, 1, 2, 3, 2, 2, 1, 1]], :auto)
    ds3_t = Dataset([Union{Missing, Int64}[3, 2, 2, 3, 1, 3], Union{Missing, Int64}[1, 1, 2, 2, 2, 3], Union{Missing, Int64}[1, 1, 2, 2, 2, 1]], :auto)
    ds4_t = Dataset([Union{Missing, Int64}[2, 1, 3], Union{Missing, Int64}[2, 1, 3], Union{Missing, Int64}[2, 1, 1]], :auto)
    ds5_t = Dataset([Union{Missing, Int64}[], Union{Missing, Int64}[], Union{Missing, Int64}[]], :auto)
    ds6 = unique!(copy(ds))
    ds7 = unique!(copy(ds), keep = :last)
    ds8 = unique!(copy(ds), keep = :none)
    ds9 = unique!(copy(ds), 1, keep = :last)
    ds10 = unique!(copy(ds), 1, keep = :none)

    @test ds6 == unique!(copy(ds), threads = false)
    @test ds7 == unique!(copy(ds), keep = :last, threads = false)
    @test ds8 == unique!(copy(ds), keep = :none, threads = false)
    @test ds9 == unique!(copy(ds), 1, keep = :last, threads = false)
    @test ds10 == unique!(copy(ds), 1, keep = :none, threads = false)

    @test ds1 == ds1_t
    @test ds2 == ds2_t
    @test ds3 == ds3_t
    @test ds4 == ds4_t
    @test ds5 == ds5_t
    @test ds6 == ds1_t
    @test ds7 == ds2_t
    @test ds8 == ds3_t
    @test ds9 == ds4_t
    @test ds10 == ds5_t
    @test unique(ds, 1, keep = :only) == ds
    @test unique(ds, 1, keep = :only, threads = false) == ds
    mft(x) = x == 1 ? missing : x==2 ? 10 : 4
    setformat!(ds, 1:3=>mft)
    ds1 = unique(ds, 2:3, mapformats = true)
    @test ds1 == unique(ds, 2:3, mapformats = true, threads = false)
    ds1_t = ds[[1,2,5,10], :]
    ds2 = unique(ds, 1, mapformats = true)
    @test ds2 == unique(ds, 1, mapformats = true, threads = false)
    ds2_t = ds[[1,3,4], :]
    ds3 = unique(ds, 1, mapformats = true, keep = :last)
    @test ds3 == unique(ds, 1, mapformats = true, keep = :last, threads = false)
    ds3_t = ds[[5,9,10], :]
    ds4 = unique(ds,1,mapformats = true, keep = :none)
    @test ds4 == unique(ds,1,mapformats = true, keep = :none, threads = false)
    ds4_t = ds[[], :]
    ds5 = unique(ds,1:2, mapformats = true, keep = :none)
    @test ds5 == unique(ds,1:2, mapformats = true, keep = :none, threads = false)
    ds5_t = ds[[4,5,7,8,10], :]
    @test ds1 == ds1_t
    @test ds2 == ds2_t
    @test ds3 == ds3_t
    @test ds4 == ds4_t
    @test ds5 == ds5_t
    @test unique(ds, 1:2, keep = :only, mapformats = true) == ds[[1,2,3,6,9], :]
    @test unique(ds, 1, keep = :only, mapformats = true) == ds
    @test unique(ds, 1:3, keep = :only, mapformats = true) == ds[[2,3,6,9], :]

    @test unique(ds, 1:2, keep = :only, mapformats = true, threads = false) == ds[[1,2,3,6,9], :]
    @test unique(ds, 1, keep = :only, mapformats = true, threads = false) == ds
    @test unique(ds, 1:3, keep = :only, mapformats = true, threads = false) == ds[[2,3,6,9], :]

    @test byrow(compare(ds1, ds1_t, mapformats =true), all)|>all
    @test byrow(compare(ds2, ds2_t, mapformats =true), all)|>all
    @test byrow(compare(ds3, ds3_t, mapformats =true), all)|>all
    @test byrow(compare(ds4, ds4_t, mapformats =true), all)|>all
    @test byrow(compare(ds5, ds5_t, mapformats =true), all)|>all

    for cols in [1, 1:2, 1:3], keepval in (:first, :last, :none, :only), mfmt in [true, false]
        @test unique(ds, cols, keep = keepval, mapformats = mfmt) == unique(view(ds, :, :), cols, keep = keepval, mapformats = mfmt)
        @test unique(ds, cols, keep = keepval, mapformats = mfmt) == unique(view(ds, :, :), cols, keep = keepval, mapformats = mfmt, threads = false)
        ds2 = ds[nrow(ds):-1:1, ncol(ds):-1:1]
        @test unique(ds2, cols, keep = keepval, mapformats = mfmt) == unique(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), cols, keep = keepval, mapformats = mfmt)
        @test unique(ds2, cols, keep = keepval, mapformats = mfmt) == unique(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), cols, keep = keepval, mapformats = mfmt, threads = false)
    end
    ds1 = Dataset(a = Union{String, Missing}["a", "b", "a", "b", "a", "b"],
                   b = Vector{Union{Int, Missing}}(1:6),
                   c = Union{Int, Missing}[1:3;1:3])
   ds = vcat(ds1, ds1)
   unique!(ds, 1, keep = :random)
   @test ds == unique(ds, 1)
end

@testset "completecases and dropmissing" begin
    ds1 = Dataset([Vector{Union{Int, Missing}}(1:4), Vector{Union{Int, Missing}}(1:4)],
                    :auto)
    ds2 = Dataset([Union{Int, Missing}[1, 2, 3, 4], ["one", "two", missing, "four"]],
                    :auto)
    ds3 = Dataset(x = Int[1, 2, 3, 4], y = Union{Int, Missing}[1, missing, 2, 3],
                    z = Missing[missing, missing, missing, missing])

    @test completecases(ds2) == .!ismissing.(ds2.x2)
    @test completecases(ds3, :x) == trues(nrow(ds3))
    @test completecases(ds3, :y) == .!ismissing.(ds3.y)
    @test completecases(ds3, :z) == completecases(ds3, [:z, :x]) ==
          completecases(ds3, [:x, :z]) == completecases(ds3, [:y, :x, :z]) ==
          falses(nrow(ds3))
    @test completecases(ds3, [:y, :x]) ==
          completecases(ds3, [:x, :y]) == .!ismissing.(ds3.y)
    @test dropmissing(ds2) == ds2[[1, 2, 4], :]
    returned = dropmissing(ds1)
    @test ds1 == returned && ds1 !== returned
    ds2b = copy(ds2)
    @test dropmissing!(ds2b) === ds2b
    @test ds2b == ds2[[1, 2, 4], :]
    ds1b = copy(ds1)
    @test dropmissing!(ds1b) === ds1b
    @test ds1b == ds1

    @test completecases(Dataset()) == Bool[]
    @test_throws MethodError completecases(Dataset(x=1), true)
    @test_throws ArgumentError completecases(ds3, :a)

    for cols in (:x2, "x2", [:x2], ["x2"], [:x1, :x2], ["x1", "x2"], 2, [2], 1:2,
                 [true, true], [false, true], :,
                 r"x2", r"x", Not(1), Not([1]), Not(Int[]), Not([]), Not(Symbol[]),
                 Not(1:0), Not([true, false]), Not(:x1), Not([:x1]))
        @test ds2[completecases(ds2, cols), :] == ds2[[1, 2, 4], :]
        @test dropmissing(ds2, cols) == ds2[[1, 2, 4], :]
        returned = dropmissing(ds1, cols)
        @test ds1 == returned && ds1 !== returned
        ds2b = copy(ds2)
        @test dropmissing!(ds2b, cols) === ds2b
        @test ds2b == ds2[[1, 2, 4], :]
        @test dropmissing(ds2, cols) == ds2b
        @test ds2 != ds2b
        ds1b = copy(ds1)
        @test dropmissing!(ds1b, cols) === ds1b
        @test ds1b == ds1
    end

    ds = Dataset(a=[1, missing, 3])
    sds = view(ds, :, :)
    @test dropmissing(sds) == Dataset(a=[1, 3])
    @test eltype(dropmissing(ds).a) == Union{Int, Missing}
    @test eltype(dropmissing(sds).a) == Union{Int, Missing}
    @test ds == Dataset(a=[1, missing, 3]) # make sure we did not mutate ds

    @test_throws MethodError dropmissing!(sds)

    ds2 = copy(ds)
    @test dropmissing!(ds) === ds
    @test dropmissing!(ds2) === ds2
    @test eltype(ds.a) == Union{Int, Missing}
    @test eltype(ds2.a) == Union{Int, Missing}
    @test ds.a.val == ds2.a.val == [1, 3]

    a = [1, 2]
    ds = Dataset(a=a)
    @test dropmissing!(ds) === ds
    @test a == ds.a.val
    dsx = dropmissing(ds)
    @test ds == ds
    @test dsx !== ds
    @test dsx.a !== ds.a.val
    @test a == ds.a.val # we did not touch ds

    b = Union{Int, Missing}[1, 2]
    ds = Dataset(b=b)
    @test eltype(dropmissing(ds).b) == Union{Int, Missing}
    @test eltype(dropmissing!(ds).b) == Union{Int, Missing}

    ds = Dataset(x1 = [1,2,3,4,missing, 1], x2 = [4,4,4,missing, missing, missing])
    fmt(x) = isequal(x, 1) ? missing : x
    setformat!(ds, 1=>fmt)
    @test completecases(ds, 1) == [true, true, true, true, false, true]
    @test completecases(ds, 1, mapformats = true) == [false, true, true, true, false, false]
    for i in 1:100
        @test completecases(ds, 1, mapformats = true, threads = true) == [false, true, true, true, false, false]
    end
    @test completecases(ds, 1:2) == [true, true, true, false, false, false]
    @test completecases(ds, 1:2, mapformats = true) == [false, true, true, false, false, false]
    @test completecases(ds, [2,1]) == [true, true, true, false, false, false]
    @test completecases(ds, [2,1], mapformats = true) == [false, true, true, false, false, false]
    @test dropmissing(ds, :) == ds[[1,2,3], :]
    @test dropmissing(ds, :, mapformats = true) == ds[[2,3], :]
    @test dropmissing(ds, :, mapformats = true, view = true) == view(ds, 2:3, :)

end

@testset "bigint tests" begin
    ds = Dataset(x = big.([1,2,1,5]), y = [1,2,1,6])
    @test unique(ds, 1) == ds[[1,2,4], :]
    @test unique(ds, 1, keep = :last) == ds[[2,3,4], :]
    @test unique(ds, 1, keep = :none) == ds[[2,4], :]
    ds[1,1]=missing
    @test completecases(ds, 1) == [false, true, true, true]
    @test completecases(ds, 1:2) == [false, true, true, true]
    @test completecases(ds, 1:2; mapformats = true) == [false, true, true, true]

    ds = Dataset(x = [89789834638463864324, 89789834638463864324, 89789834638463864325, 89789834638463864326])
    @test completecases(ds, 1) == [true, true, true, true]
    @test unique(ds, 1) == ds[[1,3,4], :]
    setformat!(ds, 1=>isodd)
    @test unique(ds, 1, mapformats = true) == ds[[1,3], :]
    @test unique(ds, 1, mapformats = true, keep = :last) == ds[[3, 4], :]
    ds[2,1]=missing
    @test completecases(ds, 1) == [true, false, true, true]
    @test dropmissing(ds, 1) == ds[[1,3,4], :]

    ds = Dataset(x = big.([1,1,5,5,1,3,2,2,2]), y=rand(9))
    @test unique(ds, 1) == ds[[1,3,6,7], :]
    @test unique(ds, 1, keep = :last) == ds[[4,5,6,9], :]

    ds = Dataset(x = [0xfffffffffffffff3, 0xfffffffffffffff1, 0xfffffffffffffff4, 0xfffffffffffffff1], y = [1,1,2,2])
    @test unique(ds, 1) == ds[[1,2,3], :]
    @test unique(ds, 1, keep = :last) == ds[[1,3,4], :]
    setformat!(ds, 1=>isodd)
    @test unique(ds, 1:2, mapformats = true) == ds[[1,3,4], :]
    @test unique(ds, 1:2, mapformats = true, keep = :none) == ds[[3,4], :]

end

@testset "filtering" begin
    ds = Dataset(x=[3, 1, 2, 1], y=["b", "c", "a", "b"])
    @test filter(ds, 1, by = >(1)) == Dataset(x=[3, 2], y=["b", "a"])
    @test filter!(ds, 1, by = >(1)) === ds == Dataset(x=[3, 2], y=["b", "a"])

    ds = Dataset(x=[3, 1, 2, 1], y=["b", "c", "a", "b"])
    @test filter(ds, :x, by = >(1)) == Dataset(x=[3, 2], y=["b", "a"])
    @test filter!(ds, :x, by = >(1)) === ds == Dataset(x=[3, 2], y=["b", "a"])

    ds = Dataset(x = [1,2,missing,1], y = ["a", "d", "c", "f"])
    @test byrow(ds, all, :, by = [isequal(1), >("a")]) == [false, false, false, true]
    setformat!(ds, 1=>isodd)
    @test byrow(mask(ds, [isequal(1), >("a")], :), all) == [false, false, false, true]
    @test byrow(mask(ds, [isequal(1), >("a")], :, mapformats = true), all) == [false, false, false, true]
    @test byrow(mask(ds, [isequal(1), ==("a")], :), all) == [true, false, false, false]
    @test byrow(mask(ds, [isequal(1), ==("a")], :, mapformats = true), all) == [true, false, false, false]
    @test byrow(mask(ds, isequal(1), 1, mapformats = true), all) == [true, false, false, true]
    @test byrow(mask(ds, ==(1), 1, mapformats = true, missings = true), all) == [true, false, true, true]
    @test isequal(mask(ds, ==(1), 1, mapformats = true, missings = missing).x , [true, false, missing, true])
    @test byrow(mask(ds, isequal(1), 1, mapformats = false), all) == [true, false, false, true]
    @test byrow(mask(ds, ==(1), 1, mapformats = false, missings = true), all) == [true, false, true, true]
    @test isequal(mask(ds, ==(1), 1, mapformats = false, missings = missing).x , [true, false, missing, true])
    setformat!(ds, 1=>iseven)
    @test byrow(mask(ds, isequal(1), 1, mapformats = true), all) == [false, true, false, false]
    @test byrow(mask(ds, ==(1), 1, mapformats = true, missings = true), all) == [false, true, true, false]
    @test isequal(mask(ds, ==(1), 1, mapformats = true, missings = missing).x , [false, true, missing, false])

    ds = Dataset(x = repeat(1:10, inner = 100), y = 10)
    @test byrow(ds, all, :, by = [>(5), ==(10)]) == [falses(500);trues(500)]
    @test byrow(mask(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), [>(5), ==(10)], [2,1]), all) == [trues(500);falses(500)]
    @test byrow(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), all, [2,1], by = [>(5), ==(10)]) == [trues(500);falses(500)]
    @test byrow(ds, all, :, by = [>(5), ==(10)], threads = false) == [falses(500);trues(500)]
    @test byrow(mask(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), [>(5), ==(10)], [2,1], threads = false), all, threads = false) == [trues(500);falses(500)]
    @test byrow(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), all, [2,1], by = [>(5), ==(10)], threads = false) == [trues(500);falses(500)]
end


text:
export Size, Rect, width, height, isinside, isintersecting,
       origin, center

import Base: size, +, *, /

struct Size{T <: Number}
  width::T
  height::T
end

struct Rect{T <: Number}
  x::T
  y::T
  width::T
  height::T
end

size(r::Rect) = Size(r.width, r.height)
width(r::Union{Rect, Size}) = r.width
height(r::Union{Rect, Size}) = r.height

+(p::Point, sz::Size) = Point(p.x + sz.width, p.y + sz.height)
+(sz::Size, p::Point) = p + sz
*(sz::Size, k::Number) = Size(sz.width * k, sz.height * k)
*(k::Number, sz::Size) = sz * k
/(sz::Size, k::Number) = Size(sz.width/k, sz.height/k)

origin(r::Rect) = Point(r.x, r.y)
center(r::Rect) = origin(r) + 0.5*size(r)

function isinside(r::Rect, p::Point)
    origin(r) ≤ p ≤ origin(r) + size(r)
end

function isintersecting(r::Rect, s::Rect)
	d = abs(center(s) - center(r))
	h1 = size(r)/2
	h2 = size(s)/2

	d.dx <= h1.width + h2.width &&
	d.dy <= h1.height + h2.height
end


text: function f 1 end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
import Base: size

"""
$(TYPEDEF)

A system type that utilizes a grid of `NX` x `NY` dual cells and `N` Lagrange forcing
points to solve the discrete Navier-Stokes equations in vorticity form. The
parameter `isstatic` specifies whether the forcing points remain static in the
grid.

# Fields
- `Re`: Reynolds number
- `U∞`: Tuple of components of free-stream velocity
- `Δx`: Size of each side of a grid cell
- `I0`: Tuple of indices of the primal node corresponding to physical origin
- `Δt`: Time step
- `rk`: Runge-Kutta coefficients
- `L`: Pre-planned discrete Laplacian operator and inverse
- `X̃`: Lagrange point coordinate data (if present), expressed in inertial coordinates
        (if static) or in body-fixed coordinates (if moving)
- `Hmat`: Pre-computed regularization matrix (if present)
- `Emat`: Pre-computed interpolation matrix (if present)
- `Vb`: Buffer space for vector data on Lagrange points
- `Fq`: Buffer space for primal cell edge data
- `Ww`: Buffer space for dual cell edge data
- `Qq`: More buffer space for dual cell edge data
- `_isstore`: flag to specify whether to store regularization/interpolation matrices

# Constructors:

`NavierStokes(Re,Δx,xlimits,ylimits,Δt
              [,U∞ = (0.0, 0.0)][,X̃ = VectorData{0}()]
              [,isstore=false][,isstatic=true][,isfilter=false]
              [,rk=TimeMarching.RK31]
              [,ddftype=Fields.Yang3])` specifies the Reynolds number `Re`, the grid
              spacing `Δx`, the dimensions of the domain in the tuples `xlimits`
              and `ylimits` (excluding the ghost cells), and the time step size `Δt`.
              The other arguments are optional. Note that `isstore` set to `true`
              would store matrix versions of the operators. This makes the method
              faster, at the cost of storage. If `isfilter` is set to true, then
              the regularization relies on a filtered version.

"""
mutable struct NavierStokes{NX, NY, N, isstatic}  #<: System{Unconstrained}
    # Physical Parameters
    "Reynolds number"
    Re::Float64
    "Free stream velocities"
    U∞::Tuple{Float64, Float64}

    # Discretization
    "Grid metadata"
    grid::Fields.PhysicalGrid{2}
    #"Grid spacing"
    #Δx::Float64
    #"Indices of the primal node corresponding to the physical origin"
    #I0::Tuple{Int,Int}
    "Time step"
    Δt::Float64
    "Runge-Kutta method"
    rk::TimeMarching.RKParams

    # Operators
    "Laplacian operator"
    L::Fields.Laplacian{NX,NY}

    # Body coordinate data, if present
    # if a static problem, these coordinates are in inertial coordinates
    # if a non-static problem, in their own coordinate systems
    X̃::VectorData{N,Float64}

    # Pre-stored regularization and interpolation matrices (if present)
    Hmat::Union{RegularizationMatrix,Nothing}
    Emat::Union{InterpolationMatrix,Nothing}
    Hmat_grad::Union{RegularizationMatrix,Nothing}
    Emat_grad::Union{InterpolationMatrix,Nothing}

    # Conditioner matrices
    Cmat::Union{AbstractMatrix,Nothing}
    Cmat_grad::Union{AbstractMatrix,Nothing}

    # Scratch space

    ## Pre-allocated space for intermediate values
    Vb::VectorData{N,Float64}
    Fq::Edges{Primal, NX, NY, Float64}
    Ww::Edges{Dual, NX, NY,Float64}
    Qq::Edges{Dual, NX, NY,Float64}

    # Flags
    _isstore :: Bool

end

function NavierStokes(Re, Δx, xlimits::Tuple{Real,Real},ylimits::Tuple{Real,Real}, Δt;
                       U∞ = (0.0, 0.0), X̃ = VectorData{0,Float64}(),
                       isstore = false,
                       isstatic = true,
                       isasymptotic = false,
                       isfilter = false,
                       rk::TimeMarching.RKParams=TimeMarching.RK31,
                       ddftype=Fields.Yang3)

    g = PhysicalGrid(xlimits,ylimits,Δx)
    NX, NY = size(g)

    α = Δt/(Re*Δx^2)

    L = plan_laplacian((NX,NY),with_inverse=true)

    Vb = VectorData(X̃)
    Fq = Edges{Primal,NX,NY,Float64}()
    Ww = Edges{Dual, NX, NY,Float64}()
    Qq = Edges{Dual, NX, NY,Float64}()
    N = length(X̃)÷2

    Hmat = nothing
    Emat = nothing
    Cmat = nothing

    Hmat_grad = nothing
    Emat_grad = nothing
    Cmat_grad = nothing

    if length(N) > 0 && isstore && isstatic
      # in this case, X̃ is assumed to be in inertial coordinates

      regop = Regularize(X̃,Δx;I0=Fields.origin(g),issymmetric=true,ddftype=ddftype)
      Hmat, Emat = RegularizationMatrix(regop,Vb,Fq)
      if isfilter
        regopfilt = Regularize(X̃,Δx;I0=Fields.origin(g),filter=true,weights=Δx^2,ddftype=ddftype)
        Ẽmat = InterpolationMatrix(regopfilt,Fq,Vb)
        Cmat = sparse(Ẽmat*Hmat)
      end
      if isasymptotic
        Hmat_grad, Emat_grad = RegularizationMatrix(regop,TensorData{N}(),grad(Fq))
        if isfilter
          Ẽmat_grad = InterpolationMatrix(regopfilt,grad(Fq),TensorData{N}())
          Cmat_grad = sparse(Ẽmat_grad*Hmat_grad)
        end
      end
    end

    # should be able to set up time marching operator here...

    #NavierStokes{NX, NY, N, isstatic}(Re, U∞, Δx, I0, Δt, rk, L, X̃, Hmat, Emat, Vb, Fq, Ww, Qq, isstore)
    NavierStokes{NX, NY, N, isstatic}(Re, U∞, g, Δt, rk, L, X̃, Hmat, Emat,
                                        Hmat_grad, Emat_grad,
                                        Cmat,Cmat_grad,
                                        Vb, Fq, Ww, Qq, isstore)
end

function Base.show(io::IO, sys::NavierStokes{NX,NY,N,isstatic}) where {NX,NY,N,isstatic}
    print(io, "Navier-Stokes system on a grid of size $NX x $NY")
end

# some convenience functions
"""
    size(sys::NavierStokes,d::Int) -> Int

Return the number of indices of the grid used by `sys` along dimension `d`.
"""
size(sys::NavierStokes{NX,NY},d::Int) where {NX,NY} = d == 1 ? NX : NY

"""
    size(sys::NavierStokes) -> Tuple{Int,Int}

Return a tuple of the number of indices of the grid used by `sys`
"""
size(sys::NavierStokes{NX,NY}) where {NX,NY} = (size(sys,1),size(sys,2))

"""
    cellsize(sys::NavierStokes) -> Float64

Return the grid cell size of system `sys`
"""
Fields.cellsize(sys::NavierStokes) = cellsize(sys.grid)

"""
    origin(sys::NavierStokes) -> Tuple{Int,Int}

Return a tuple of the indices of the primal node that corresponds to the
physical origin of the coordinate system used by `sys`. Note that these
indices need not lie inside the range of indices occupied by the grid.
For example, if the range of physical coordinates occupied by the grid
is (1.0,3.0) x (2.0,4.0), then the origin is not inside the grid.
"""
Fields.origin(sys::Systems.NavierStokes) = origin(sys.grid)

# Basic operators for any Navier-Stokes system

# Integrating factor -- rescale the time-step size
Fields.plan_intfact(Δt,w,sys::NavierStokes{NX,NY}) where {NX,NY} =
        Fields.plan_intfact(Δt/(sys.Re*cellsize(sys)^2),w)

# RHS of Navier-Stokes (non-linear convective term)
function TimeMarching.r₁(w::Nodes{Dual,NX,NY,T},t,sys::NavierStokes{NX,NY}) where {NX,NY,T}

  Ww = sys.Ww
  Qq = sys.Qq
  L = sys.L
  Δx⁻¹ = 1/cellsize(sys)

  grid_interpolate!(Qq,curl(L\w)) # -velocity, on dual edges
  Qq.u .-= sys.U∞[1]
  Qq.v .-= sys.U∞[2]

  return rmul!(divergence(Qq∘grid_interpolate!(Ww,w)),Δx⁻¹) # -∇⋅(wu)

end

# RHS of Navier-Stokes (non-linear convective term)
function TimeMarching.r₁(w::Nodes{Dual,NX,NY,T},t,sys::NavierStokes{NX,NY},U∞::RigidBodyMotions.RigidBodyMotion) where {NX,NY,T}

  Ww = sys.Ww
  Qq = sys.Qq
  L = sys.L
  Δx⁻¹ = 1/cellsize(sys)

  grid_interpolate!(Qq,curl(L\w)) # -velocity, on dual edges
  _,ċ,_,_,_,_ = U∞(t)
  Qq.u .-= real(ċ)
  Qq.v .-= imag(ċ)

  return rmul!(divergence(Qq∘grid_interpolate!(Ww,w)),Δx⁻¹) # -∇⋅(wu)

end

# Operators for a system with a body

# RHS of a stationary body with no surface velocity
function TimeMarching.r₂(w::Nodes{Dual,NX,NY,T},t,sys::NavierStokes{NX,NY,N,true}) where {NX,NY,N,T}
    ΔV = VectorData(sys.X̃)
    ΔV.u .-= sys.U∞[1]
    ΔV.v .-= sys.U∞[2]
    return ΔV
end

function TimeMarching.r₂(w::Nodes{Dual,NX,NY,T},t,sys::NavierStokes{NX,NY,N,true},U∞::RigidBodyMotions.RigidBodyMotion) where {NX,NY,N,T}
    ΔV = VectorData(sys.X̃)
    _,ċ,_,_,_,_ = U∞(t)
    ΔV.u .-= real(ċ)
    ΔV.v .-= imag(ċ)
    return ΔV
end

# Constraint operators, using stored regularization and interpolation operators
# B₁ᵀ = CᵀEᵀ, B₂ = -ECL⁻¹
TimeMarching.B₁ᵀ(f::VectorData{N},sys::NavierStokes{NX,NY,N,C}) where {NX,NY,N,C} = Curl()*(sys.Hmat*f)
TimeMarching.B₂(w::Nodes{Dual,NX,NY,T},sys::NavierStokes{NX,NY,N,C}) where {NX,NY,T,N,C} = -(sys.Emat*(Curl()*(sys.L\w)))

# Constraint operators, using non-stored regularization and interpolation operators
TimeMarching.B₁ᵀ(f::VectorData{N},regop::Regularize,sys::NavierStokes{NX,NY,N,false}) where {NX,NY,N} = Curl()*regop(sys.Fq,f)
TimeMarching.B₂(w::Nodes{Dual,NX,NY,T},regop::Regularize,sys::NavierStokes{NX,NY,N,false}) where {NX,NY,T,N} = -(regop(sys.Vb,Curl()*(sys.L\w)))

# Constraint operator constructors
# Constructor using stored operators
TimeMarching.plan_constraints(w::Nodes{Dual,NX,NY,T},t,sys::NavierStokes{NX,NY,N,true}) where {NX,NY,T,N} =
                    (f -> TimeMarching.B₁ᵀ(f,sys),w -> TimeMarching.B₂(w,sys))

# Constructor using non-stored operators
function TimeMarching.plan_constraints(w::Nodes{Dual,NX,NY,T},t,sys::NavierStokes{NX,NY,N,false}) where {NX,NY,T,N}
  regop = Regularize(sys.X̃,Fields.cellsize(sys);I0=Fields.origin(sys),issymmetric=true)

  return f -> TimeMarching.B₁ᵀ(f,regop,sys),w -> TimeMarching.B₂(w,regop,sys)
end

# compute physical values of fields
vorticity(w::Nodes{Dual},sys) = w/cellsize(sys)
velocity(w::Nodes{Dual},sys) = -curl(sys.L\w)
streamfunction(w::Nodes{Dual},sys) = -cellsize(sys)*(sys.L\w)
nl(w::Nodes{Dual},sys) = TimeMarching.r₁(w,0.0,sys)/cellsize(sys)

include("navierstokes/systemutils.jl")

include("navierstokes/movingbody.jl")


text:
ALL_TASK_OFFSETS[("darwin-x86_64", v"1.6.1")] = Dict(
    :END => 264,
    :_isexception => 58,
    :_state => 56,
    :bufsz => 240,
    :copy_stack => 248,
    :ctx => 80,
    :donenotify => 24,
    :eh => 72,
    :excstack => 64,
    :gcstack => 256,
    :logstate => 40,
    :next => 0,
    :prio => 62,
    :queue => 8,
    :result => 32,
    :sizeof_ctx => 152,
    :start => 48,
    :sticky => 57,
    :stkbuf => 232,
    :tid => 60,
    :tls => 16,
    :tls_base_context => 6680,
)


text:
using BenchmarkTools, CvxCompress, LinearAlgebra, Pkg

wavefields = Pkg.artifact"wavefields"
p3d = read!(joinpath(wavefields, "field_3D_for_compression_Float32_201x401x101.bin"), Array{Float32,3}(undef,201,401,101))
p2d = read!(joinpath(wavefields, "field_2D_for_compression_Float32_352x688.bin"), Array{Float32,2}(undef,352,688))

p3d_compressed = Vector{UInt32}(undef, length(p3d))
p2d_compressed = Vector{UInt32}(undef, length(p2d))

p3d_decompressed = similar(p3d)
p2d_decompressed = similar(p2d)

compressor_3d = CvxCompressor((32,32,32), 1e-2)
compressor_2d = CvxCompressor((32,32), 1e-2)

const SUITE = BenchmarkGroup()

n_2d = compress!(p2d_compressed, compressor_2d, p2d)
decompress!(p2d_decompressed, compressor_2d, p2d_compressed, n_2d)
ratio_2d = length(p2d) / n_2d
snr_2d = 20*log10(norm(p2d) / norm(p2d - p2d_decompressed))

SUITE["compression, 2D"] = BenchmarkGroup([Dict("compression ratio"=>ratio_2d, "signal-to-noise ratio"=>snr_2d, "ncells"=>length(p2d))])
SUITE["compression, 2D"]["compress, F32"] = @benchmarkable compress!($p2d_compressed, $compressor_2d, $p2d)
SUITE["compression, 2D"]["decompress, F32"] = @benchmarkable decompress!($p2d_decompressed, $compressor_2d, $p2d_compressed, $n_2d)

n_3d = compress!(p3d_compressed, compressor_3d, p3d)
decompress!(p3d_decompressed, compressor_3d, p3d_compressed, n_3d)
ratio_3d = length(p3d) / n_3d
snr_3d = 20*log10(norm(p3d) / norm(p3d - p3d_decompressed))

SUITE["compression, 3D"] = BenchmarkGroup([Dict("compression ratio"=>ratio_3d, "signal-to-noise ratio"=>snr_3d, "ncells"=>length(p3d))])
SUITE["compression, 3D"]["compress, F32"] = @benchmarkable compress!($p3d_compressed, $compressor_3d, $p3d)
SUITE["compression, 3D"]["decompress, F32"] = @benchmarkable decompress!($p3d_decompressed, $compressor_3d, $p3d_compressed, $n_3d)

include(joinpath(pkgdir(CvxCompress), "benchmark", "mcells_per_second.jl"))

SUITE


text:
"""
    ComputeBias(z::Float64, ::PiecewiseBias,
    ConvolvedDensity::AbstractConvolvedDensity)

This function evaluate the piecewise bias, given by ``\\sqrt{1+z}``, for a given
redshift.

"""
function ComputeBias(z::Float64, Piecewise::PiecewiseBias,
    ConvolvedDensity::AbstractConvolvedDensity)
    idx = BinSearch(z, ConvolvedDensity.ZBinArray)
    bias = sqrt(1+(ConvolvedDensity.ZBinArray[idx]+
    ConvolvedDensity.ZBinArray[idx+1])/2) * Piecewise.BiasMultiplier[idx]
    return bias
end

function ComputeBias(z::Float64, Bias::EuclidBias,
    ConvolvedDensity::AbstractConvolvedDensity)
    bias = Bias.A+Bias.B/(1+exp((Bias.D-z)*Bias.C))
    return bias
end

"""
    ComputeBiasGrid!(CosmologicalGrid::CosmologicalGrid,
    gcWeightFunction::GCWeightFunction, Bias::AbstractBias,
    ConvolvedDensity::AbstractConvolvedDensity)

This function evaluate the bias, over the cosmological redshift grid.
"""
function ComputeBiasGrid!(cosmologicalGrid::CosmologicalGrid,
    gcWeightFunction::GCWeightFunction,
    ConvolvedDensity::AbstractConvolvedDensity)
    for (zidx, zvalue) in enumerate(cosmologicalGrid.ZArray)
        gcWeightFunction.BiasArray[:, zidx] .=
        ComputeBias(zvalue, gcWeightFunction.BiasKind, ConvolvedDensity)
    end
end


text:
@inline function _degenerate(h::T, sinφ1::T, sinφ2::T) where T
    CUDA.max(zero(h), h) < _etol(h) ||
    1 - _abs(sinφ1) < _etol(h) ||
    1 - _abs(sinφ2) < _etol(h) ||
    _abs(sinφ1 - sinφ2) < _etol(h)
end

@inline function _distance(ξ::CuPosition{T}, n::CuPosition{T}, dist::T) where T
    _dot(ξ, n) - dist
end

@inline function _logterm(χ2::T, sinφ::T) where T
    t1 = _sqrt(one(χ2) - χ2 * sinφ * sinφ)
    t2 = _sqrt(one(χ2) - χ2) * sinφ
    (t1 + t2) / (t1 - t2)
end

@inline function _project(ξ::CuPosition{T}, n::CuPosition{T}, dist::T) where T
    _sub(ξ, _smul(n, dist))
end

@inline function _laplacepot(
    ξ   ::CuPosition{T},
    elem::CuTriangle{T},
    dist    ::T,
    pot     ::F
) where {T, F <: Function}
    _laplacepot(ξ, elem.v1, elem.v2, elem.normal, dist, pot) +
    _laplacepot(ξ, elem.v2, elem.v3, elem.normal, dist, pot) +
    _laplacepot(ξ, elem.v3, elem.v1, elem.normal, dist, pot)
end

function _laplacepot(
    ξ   ::CuPosition{T},
    x1  ::CuPosition{T},
    x2  ::CuPosition{T},
    n   ::CuPosition{T},
    dist::T,
    pot ::F
) where {T, F <: Function}
    u1  = _sub(x1, ξ)
    u1n = _norm(u1)

    u2  = _sub(x2, ξ)
    u2n = _norm(u2)

    v  = _sub(u2, u1)
    vn = _norm(v)

    sinφ1 = _clamp(_cos(u1, u1n, v, vn), -one(vn), one(vn))
    sinφ2 = _clamp(_cos(u2, u2n, v, vn), -one(vn), one(vn))

    h = _cathethus(u1n, sinφ1)

    _degenerate(h, sinφ1, sinφ2) ?
        zero(h) :
        _sign(_dot(_cross(u1, u2), n)) * pot(sinφ1, sinφ2, h, dist)
end

@inline function _laplacepot_single_plane(sinφ1::T, sinφ2::T, h::T, ::T) where T
    2 \ h * _log((1+sinφ2) * (1-sinφ1) / ((1-sinφ2) * (1+sinφ1)))
end

@inline function _laplacepot_single_space(sinφ1::T, sinφ2::T, h::T, d::T) where T
    χ2 = d * d / (d * d + h * h)
    χ  = _sqrt(χ2)
    2 \ h * _log(_logterm(χ2, sinφ2) / _logterm(χ2, sinφ1)) +
        _abs(d) * (_asin(χ * sinφ2) - _asin(sinφ2) - _asin(χ * sinφ1) + _asin(sinφ1))
end

@inline function _laplacepot_double_space(sinφ1::T, sinφ2::T, h::T, d::T) where T
    χ = _abs(d) / _sqrt(d * d + h * h)
    _sign(d) * (_asin(χ * sinφ1) - _asin(sinφ1) - _asin(χ * sinφ2) + _asin(sinφ2))
end

@inline function laplacepot_single(
    ξ   ::CuPosition{T},
    elem::CuTriangle{T}
) where T
    dist = _distance(ξ, elem.normal, elem.distorig)
    _abs(dist) < _etol(dist) ?
        _laplacepot(ξ, elem, dist, _laplacepot_single_plane) :
        _laplacepot(_project(ξ, elem.normal, dist), elem, dist, _laplacepot_single_space)
end

@inline function laplacepot_double(
    ξ   ::CuPosition{T},
    elem::CuTriangle{T}
) where T
    dist = _distance(ξ, elem.normal, elem.distorig)
    _abs(dist) < _etol(dist) ?
        zero(dist) :
        _laplacepot(_project(ξ, elem.normal, dist), elem, dist, _laplacepot_double_space)
end


text:
using Plots

struct Line
    pos
    dir
    len
end

@recipe function f(line::Line)
    # layout := (2, 2)

    @series begin
        # subplot := 1
        [line.pos[1], line.pos[1] + line.dir[1] * line.len],
        [line.pos[2], line.pos[2] + line.dir[2] * line.len]
    end
end

line = Line([0,0,0], [1,10,1], 10)

plot(line)

plot( plot([1, 2], [3, 4]), plot([2, 4], [1, 4]) )


text:
function test_dual_problems()
    n  = rand(50:100)
    m  = rand(50:100)
    X1 = rand(Normal(-1.5,1), n, 2)
    X2 = rand(Normal(1.5,1), m, 2)
    X  = vcat(X1, X2)
    y  = vcat(ones(n), zeros(m)) .== 1

    τ    = 0.2*rand() + 0.3
    K    = rand(1:10)
    C    = 0.5*rand() + 1
    ϑ1   = 0.4*rand() + 0.5
    ϑ2   = 0.5*rand() + 1

    @testset "PatMat with $surrogate loss" for surrogate in [Hinge, Quadratic]
        l1    = surrogate(ϑ1);
        l2    = surrogate(ϑ2);
        model = PatMat(τ, C, l1, l2);
        data  = Dual(model, X, y);

        test_dual(model, data, 30000, 200000, Descent(0.001))
    end

    @testset "PatMatNP with $surrogate loss" for surrogate in [Hinge, Quadratic]
        l1    = surrogate(ϑ1);
        l2    = surrogate(ϑ2);
        model = PatMatNP(τ, C, l1, l2);
        data  = Dual(model, X, y);

        test_dual(model, data, 30000, 200000, Descent(0.001))
    end

    @testset "TopPushK with $surrogate loss" for surrogate in [Hinge, Quadratic]
        l     = surrogate(ϑ1);
        model = TopPushK(K, C, l);
        data  = Dual(model, X, y);

        test_dual(model, data, 30000, 200000, Descent(0.001))
    end

    @testset "TopPush with $surrogate loss" for surrogate in [Hinge, Quadratic]
        l     = surrogate(ϑ1);
        model = TopPush(C, l);
        data  = Dual(model, X, y);

        test_dual(model, data, 30000, 200000, Descent(0.001))
    end
end

function test_dual(model::AbstractModel, data::Dual{<:DTrain}, maxiter::Integer, maxiter2::Integer, optimizer::Any; atol::Real = 1e-2)

    sol1 = solve(Gradient(maxiter = maxiter, optimizer = optimizer, verbose = false), model, data)
    L1   = AccuracyAtTopKernels.objective(model, data, sol1)

    sol2 = solve(Coordinate(maxiter = maxiter2, verbose = false), model, data)
    L2   = AccuracyAtTopKernels.objective(model, data, sol2)

    sol3 = solve(General(), model, data)
    L3   = AccuracyAtTopKernels.objective(model, data, sol3)

    @testset "gradient solver" begin
        @testset "feasibility" begin
            isfeasible(model, data, sol1)
        end
        @testset "optimality" begin
            @test L1 ≈ L3 atol = atol
        end
    end
    @testset "coordinate solver" begin
        @testset "feasibility" begin
            isfeasible(model, data, sol2)
        end
        @testset "optimality" begin
            @test L2 ≈ L3 atol = atol
        end
    end
    @testset "general solver feasibility" begin
        isfeasible(model, data, sol3)
    end
end

function isfeasible(model::AbstractPatMat{<:Hinge}, data::Dual{<:DTrain}, solution::NamedTuple; atol::Real = 1e-5)
    α, β, δ = solution.α, solution.β, solution.δ
    @test sum(α) ≈ sum(β) atol = atol
    @test maximum(α) <= model.l1.ϑ*model.C + atol
    @test maximum(β) <= model.l2.ϑ*δ + atol
    @test minimum(α) >= - atol
    @test minimum(β) >= - atol
end

function isfeasible(model::AbstractPatMat{<:Quadratic}, data::Dual{<:DTrain}, solution::NamedTuple; atol::Real = 1e-5)
    α, β, δ = solution.α, solution.β, solution.δ
    @test sum(α) ≈ sum(β) atol = atol
    @test minimum(α) >= - atol
    @test minimum(β) >= - atol
    @test δ >= - atol
end

function isfeasible(model::M, data::Dual{<:DTrain}, solution::NamedTuple; atol::Real = 1e-5) where {M<:AbstractTopPushK{<:Hinge}}
    α, β = solution.α, solution.β
    @test sum(α) ≈ sum(β) atol = atol
    @test maximum(α) <= model.l.ϑ*model.C + atol
    @test minimum(α) >= 0 - atol
    @test minimum(β) >= 0 - atol
    M <: TopPushK &&  @test maximum(β) <= sum(α)/model.K + atol
end

function isfeasible(model::M, data::Dual{<:DTrain}, solution::NamedTuple; atol::Real = 1e-5) where {M<:AbstractTopPushK{<:Quadratic}}
    α, β = solution.α, solution.β
    @test sum(α) ≈ sum(β) atol = atol
    @test minimum(α) >= 0 - atol
    @test minimum(β) >= 0 - atol
    M <: TopPushK && @test maximum(β) <= sum(α)/model.K + atol
end


text:
using myJuliaUtils
using Test
using LinearAlgebra
using Random
using Distributions
# using BenchmarkTools
using DelimitedFiles

# comparing the procedure to matlab

path = dirname(@__FILE__)
# path = "/Users/zymantas/TVPVARPkg/test/TestWish"

T = readdlm("$path/T.txt", ',');
T = Int(T[1])
if T > 80
    error("matlab uses a different procedure for high degrees of freedom")
end
draw = readdlm("$path/draw.txt", ',');
aa = readdlm("$path/WW.txt", ',');
MatlabInvWish = readdlm("$path/matlabinvwish.txt", ',');
MatlabWish = readdlm("$path/matlabwish.txt", ',');

A = draw * cholesky(aa).U
JuliaWish = A'A

B = copy(draw)
A = cholesky(aa).L / qr!(B).R
JuliaInvWish = A * A'

# C2 = inv(Symmetric(aa));
# C3= draw * cholesky(C2).U
# C3=C3'C3
# D2 = inv(Symmetric(C3))
# @test D2 ≈ MatlabInvWish

@testset "Compared to Matlab" begin
    @test JuliaWish ≈ MatlabWish
    @test JuliaInvWish ≈ MatlabInvWish
end

# testing wishart and inverse wishart for univariate compared to multivariate
WW = randn(1, 1)
WW = WW'WW
T = rand(1:1000)

SomeSeed = rand(1:10000)

Random.seed!(SomeSeed);
C = wish(WW, T)

Random.seed!(SomeSeed);
C1 = wish(vec(WW), T)

Random.seed!(SomeSeed);
D = iwish(WW, T)

Random.seed!(SomeSeed);
D1 = iwish(vec(WW), T)

@testset "matching univariate to multivariate" begin
    @test C ≈ [C1]
    @test D ≈ [D1]
end

# testing wishart and inverse wishart for univariate
WW = randn(1, 1)
WW = WW'WW
T = rand(1:1000)

Random.seed!(SomeSeed);
A = rand(Wishart(T, WW))
Random.seed!(SomeSeed);
C = wish(WW, T)

Random.seed!(SomeSeed);
A = rand(InverseWishart(T, WW))

Random.seed!(SomeSeed);
B = iwish(WW, T)

Random.seed!(SomeSeed);
C2 = inv(WW);
C3 = wish(C2, T);
C = inv(C3)

# for univariate case where WW is a vector
WW = [1 + rand() * 10000]
T = rand(1:1000)

SomeSeed = rand(1:10000)

Random.seed!(SomeSeed);
A_v = rand(InverseWishart(T, reshape(WW, 1, 1)))

Random.seed!(SomeSeed);
B_v = iwish(WW, T)

Random.seed!(SomeSeed);
C2 = inv.(WW);
C3 = wish(C2, T);
C_v = inv.(C3)

@testset "univariate" begin
    # @test A ≈ B
    @test B ≈ C

    # @test A_v ≈ B_v
    @test B_v ≈ C_v
end

# # testing for multivariate wishart

# K = round(Int, 1 + rand() * 6)

# WW = randn(K, K)
# WW = WW'WW
# T = round(Int, 1 + rand() * 1000)

# # testing the equality between wish and iwish
# Random.seed!(SomeSeed);
# A = rand(InverseWishart(T, WW))
# Random.seed!(SomeSeed);
# B = iwish(WW, T)
# Random.seed!(SomeSeed);
# C2 = inv(Symmetric(WW));
# C3 = wish(C2, T);
# C = inv(C3)

# @testset "multivariate" begin
#     @test_broken A ≈ B
#     @test_broken B ≈ C
# end

# using Law of large numbers

reps = 2000000
K = rand(1:5)

WW = randn(K, K)
WW = WW'WW
T = rand(K+1:300)

SomeSeed = rand(1:10000)
Random.seed!(SomeSeed)
a1 = wish(WW, T)

Random.seed!(SomeSeed);
a2 = rand(Wishart(T, WW))

cont1 = zeros(K, K);
cont2 = copy(cont1);
cont3 = copy(cont1);
cont4 = copy(cont1);

wD=Wishart(T, WW)
iwD=InverseWishart(T, WW)

for i = 1:reps
    global cont1 .+= wish(WW, T);
    global cont2 .+= rand(wD);

    global cont3 .+= iwish(WW, T);
    global cont4 .+= rand(iwD);
end
cont1 ./= reps;
cont2 ./= reps;
cont3 ./= reps;
cont4 ./= reps;

@testset "based on LLN" begin

    # @test maximum(abs.(round.(cont1 - cont2; digits=0)))==0.0
    @test maximum(abs.(round.(cont1 - mean(Wishart(T, WW)); digits=0)))==0.0
    @test maximum(abs.(round.(mean(Wishart(T, WW)) - cont2; digits=0)))==0.0

    # @test maximum(abs.(round.(cont3 - cont4; digits=4)))==0.0
    @test maximum(abs.(round.(cont3 - mean(InverseWishart(T, WW)); digits=3)))==0.0
    @test maximum(abs.(round.(mean(InverseWishart(T, WW)) - cont4; digits=3)))==0.0

end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
using DataFrames, CSV
using Profile
include("./sdistance.jl") # 必ずダブルクオーテーション
include("./draw.jl")
include("./LevelSet.jl")
include("./environments.jl")
include("../test/APT.jl")
import .Sdistance: computing_bench
import .Draw: parformance_graphs
using .LevelSet
using .APT

# This script run by test.sh

# === profiling ===

# @profile signedDistance2D("./test/mock_csv_data/interface.csv",parse(Int, ARGS[1]))
# @profile signedDistance2D("./test/mock_csv_data/interface.csv",parse(Int, ARGS[1]), "multi")
# Profile.print()
# open("prof.txt", "w") do s
#     Profile.print(IOContext(s, :displaysize => (24, 500)))
# end

# === memory size ===
# p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/interface.csv")
# p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/interface.csv", "multi")
p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/multiple_curves.csv", "multi")

# p = @allocated signedDistance2D("./test/mock_csv_data/interface.csv", parse(Int, ARGS[1]), "multi")
# p =  signedDistance2D("./test/mock_csv_data/interface.csv", parse(Int, ARGS[1]))
# p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/interface.csv", "multi")
# p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/mini_interface.csv")
# p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/mini_interface.csv", "multi")
# p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/infinity_shaped.csv", "multi")
# p = @allocated computing_bench(parse(Int, ARGS[1]), "./test/mock_csv_data/double_circle.csv", "multi")
# q = signedDistance2D( "./test/mock_csv_data/multiple_curves.csv", parse(Int, ARGS[1]), "multi")
# plots_contours([i for i = 50:50:300], "./test/mock_csv_data/interface.csv", "multi")
# plots_contours([i for i = 50:50:300], "./test/mock_csv_data/interface.csv")
# plots_wireframe([i for i = 50:50:300], "./test/mock_csv_data/interface.csv", "multi")

# DataFrame(p, :auto) |> CSV.write("./test/result/interface_result_n500.csv", header=false)
# println("\nmemory size: ",p/(1024*1024), " MB")

max_stars_repo_name: jabelic/senior-thesis-julia

max_issues_repo_name: jabelic/senior-thesis-julia

max_forks_repo_name: jabelic/senior-thesis-julia


text:
@testset "XDAReader" begin
    mesh = readXDAMesh("data/square.xda")

    @test mesh.nodes[3].coords[1] == 0.5
    @test mesh.nodes[3].coords[2] == 0.5
end


text:
import Base.show
function name(h::Hid)
    iob = IOBuffer()
    print(iob, h.vendor_name, " ")
    print(iob, h.product_name)
    print(iob, "(vID=0x", string(h.vID, base=16), ", ")
    print(iob, "pID=0x", string(h.pID, base=16), ")")
    String(take!(iob))
end

# Long form, as in display(h) or REPL: h enter
function Base.show(io::IO, ::MIME"text/plain", h::Hid)
    print(io, "Hid(")
    print(io, "\tvendor_name:\t ", h.vendor_name, "\n")
    print(io, "\tproduct_name: ", h.product_name, "\n")
    print(io, "\tvID:\t 0x", string(h.vID, base=16), "\n")
    print(io, "\tpID:\t 0x", string(h.pID, base=16), "\n")
    print(io, "\tdevice_path:\t ", h.device_path, ")")
end

# Short form, as in print(stdout, h)
function Base.show(io::IO, h::Hid)
    ioc = IOContext(io)
    if get(ioc, :color, false)
        printstyled(ioc, h.vendor_name, color=:green, bold = false)
        print(ioc, " ")
        printstyled(ioc, h.product_name, color=:green, bold = true)
        print(ioc, "(vID=0x", string(h.vID, base=16), ", ")
        print(ioc, "pID=0x", string(h.pID, base=16), ")")
    else
        print(ioc, name(h))
    end
end

# The default Juno / Atom display works nicely with standard output
Base.show(io::IO, ::MIME"application/prs.juno.inline", hid::Hid) = Base.show(io, hid)

# Vector
function Base.show(io::IO, ::MIME"text/plain", hids::Array{Hid,1})
    print(io, length(hids), "-element ")
    println(io, typeof(hids), ":")
    for (i, h) in enumerate(hids)
        println(io, " [", i, "]\t", h)
    end
end


text:
module problem124

description = """
Ordered radicals
The radical of n, rad(n), is the product of the distinct prime factors of n. For example, 504 = 23 × 32 × 7, so rad(504) = 2 × 3 × 7 = 42.
Let E(k) be the kth element in the sorted n column; for example, E(4) = 8 and E(6) = 9.

If rad(n) is sorted for 1 ≤ n ≤ 100000, find E(10000).
""

# Calculate rad for all values between 1 and n using algorithm similar to sieve of eratosthenes
function radsieve(n)
  rad = ones(Int, n)
  for i = 2:n
    if rad[i] == 1
      for j = i:i:n
        rad[j] *= i
      end
    end
  end
  return rad
end

function solve(n=100000, k=10000)
  radsorted = sort(collect(enumerate(radsieve(n))), by = pr -> pr[2], alg=MergeSort)
  return radsorted[k][1]
end

using Base.Test

@test radsieve(10) == [1, 2, 3, 2, 5, 6, 7, 2, 3, 10]
@test solve(10, 4) == 8
@test solve(10, 6) == 9

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


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

"""
    nrbkntins(nurbs::NURBS{I,F}, iknots::Union{Vector{F}, Vector{Vector{F}}}
             )::NURBS{I,F} where {I<:Integer,F<:AbstractFloat}

Insert a single or multiple knots into a NURBS curve or surface.

Inserts knots into the NURBS data structure, these can be knots at new
positions or at the location of existing knots to increase the multiplicity.
Note that the knot multiplicity cannot be increased beyond the order of the
spline. Knots along the V direction can only inserted into NURBS surfaces,
not curves that are always defined along the U direction. This function uses
the B-Spline function [`bspkntins`](@ref).

# Arguments:
- `nurbs`: NURBS structure (curve or surface see [`nrbmak`](@ref))
- `iknots`: vektor of knots or pair of vector of knots to insert into
  respective U,V directions of the structure

# Output:
- `enurbs`: new NURBS structure with degrees elevated

# Examples:
```julia
julia> icrv = nrbkntins(crv,iknots);
julia> isrf = nrbkntins(srf,[iuknots,ivknots]);
```

!!! info "Note:"
    No knot multiplicity will be increased beyond the order of the spline.

Insert two knots into a curve, one at 0.3 and another twice at 0.4.
```julia
julia> icrv = nrbkntins(crv, vec([0.3 0.4 0.4]))
```

Insert into a surface two knots as (1) into the U knot sequence and one knot
into the V knot sequence at 0.5.
```julia
julia> isrf = nrbkntins(srf, [vec([0.3 0.4 0.4]),[0.5]])
```

# Reference
- see also: [`bspkntins`](@ref)
"""
function nrbkntins(nurbs::NURBS{I,F},iknots::Union{Vector{F},Vector{Vector{F}}}
                  )::NURBS{I,F} where {I<:Integer,F<:AbstractFloat}
# TODO: this ported version does not do all the tests of the matlab version,
#       instead the type features of julia are used -> check if the matlab
#       errors get caught

degree = nurbs.order .- 1;

if typeof(nurbs)<:NURBS2D
  if typeof(iknots)<:Vector{F} || length(iknots) != 2
    throw(ArgumentError("iknots must be a vector of two Vector{F}!"));
  end # if
  # NURBS represents a surface
  num1 = nurbs.number[1];
  num2 = nurbs.number[2];
  knots = Vector{Vector{F}}(undef,2);

  # Degree elevate along the v direction
  if isempty(iknots[2])
    # TODO: i think julia might not even allow empty vectors
    coefs = nurbs.coefs;
    knots[2] = nurbs.knots[2];
  else
    coefs = reshape(nurbs.coefs,(4*num1,num2));
    coefs, knots[2] = bspkntins(degree[2],coefs,nurbs.knots[2],iknots[2]);
    num2 = size(coefs,2);
    coefs = reshape(coefs,(4, num1, num2));
  end # if

  # Degree elevate along the u direction
  if isempty(iknots[1])
    # TODO: i think julia might not even allow empty vectors
    knots[1] = nurbs.knots[1];
  else
    coefs = permutedims(coefs,[1,3,2]);
    coefs = reshape(coefs,(4*num2,num1));
    coefs, knots[1] = bspkntins(degree[1],coefs,nurbs.knots[1],iknots[1]);
    coefs = reshape(coefs,(4,num2,size(coefs,2)));
    coefs = permutedims(coefs,[1,3,2]);
  end # if

  # build the new nurbs surface
  return nrbmak(coefs,knots);
else
  if !(typeof(iknots)<:Vector{F})
    throw(ArgumentError("Valid input of iknots for curves is Vector{F}!"));
  end # if

  # NURBS represents a curve
  if isempty(iknots)
    # TODO: i think julia might not even allow empty vectors
    coefs = nurbs.coefs;
    knots = nurbs.knots;
  else
    coefs, knots = bspkntins(degree,nurbs.coefs,nurbs.knots,iknots);
  end # if

  # build the new nurbs curve
  return nrbmak(coefs,knots);
end # if
end # nrbkntins


text:
function binomial_coefficient(::Type{CombinatoricsUtils}, arg0::jint, arg1::jint)
    return jcall(CombinatoricsUtils, "binomialCoefficient", jlong, (jint, jint), arg0, arg1)
end

function binomial_coefficient_double(::Type{CombinatoricsUtils}, arg0::jint, arg1::jint)
    return jcall(CombinatoricsUtils, "binomialCoefficientDouble", jdouble, (jint, jint), arg0, arg1)
end

function binomial_coefficient_log(::Type{CombinatoricsUtils}, arg0::jint, arg1::jint)
    return jcall(CombinatoricsUtils, "binomialCoefficientLog", jdouble, (jint, jint), arg0, arg1)
end

function check_binomial(::Type{CombinatoricsUtils}, arg0::jint, arg1::jint)
    return jcall(CombinatoricsUtils, "checkBinomial", void, (jint, jint), arg0, arg1)
end

function combinations_iterator(::Type{CombinatoricsUtils}, arg0::jint, arg1::jint)
    return jcall(CombinatoricsUtils, "combinationsIterator", JIterator, (jint, jint), arg0, arg1)
end

function factorial(::Type{CombinatoricsUtils}, arg0::jint)
    return jcall(CombinatoricsUtils, "factorial", jlong, (jint,), arg0)
end

function factorial_double(::Type{CombinatoricsUtils}, arg0::jint)
    return jcall(CombinatoricsUtils, "factorialDouble", jdouble, (jint,), arg0)
end

function factorial_log(::Type{CombinatoricsUtils}, arg0::jint)
    return jcall(CombinatoricsUtils, "factorialLog", jdouble, (jint,), arg0)
end

function stirling_s2(::Type{CombinatoricsUtils}, arg0::jint, arg1::jint)
    return jcall(CombinatoricsUtils, "stirlingS2", jlong, (jint, jint), arg0, arg1)
end


text:
if !@isdefined(Jpycachesim)
    include("Jpycachesim.jl")
    using .Jpycachesim
end

#pcs_get_cachesim_from_file("cachesim/cachedef");

# This is my laptop
# l3 = Cache("L3", 12288, 16, 64, "LRU");
# l2 = Cache("L2", 256, 4, 64, "LRU");
# l1 = Cache("L1", 32, 8, 64, "LRU");

al3 = Cache("L3", 32, 1, 64, "LRU");
al2 = Cache("L2", 16, 1, 64, "LRU");
al1 = Cache("L1", 8, 1, 64, "LRU");

bl3 = Cache("L3", 32, 1, 64, "LRU");
bl2 = Cache("L2", 16, 1, 64, "LRU");
bl1 = Cache("L1", 8, 1, 64, "LRU");

c1 = pcs_build_cache([al1,al2,al3]);
c2 = pcs_build_cache([bl1,bl2,bl3]);

pcs_load(0,1,c1);
(m1,h1) = pcs_get_stats(c1);

pcs_load(0,130,c2);
(m2,h2) = pcs_get_stats(c2);

pcs_load(0,1,c1);
(m3,h3) = pcs_get_stats(c1);

println("test stats: should be miss=1, ?, 0, hit = 0, ?, 1");
println("miss= "*string(m1[1])*", "*string(m2[1])*", "*string(m3[1]));
println("hit= "*string(h1[1])*", "*string(h2[1])*", "*string(h3[1]));

pcs_load(0,8,c1);
pcs_store(0,8,c1);
pcs_load(28, 32,c1);

pcs_load(0,128,c2);
pcs_store(0,16,c2);
pcs_load(28, 1024,c2);

pcs_print_stats(c1);
pcs_print_stats(c2);

cobj = unsafe_load(c1);
lev2 = unsafe_load(cobj.load_from);
lev3 = unsafe_load(lev2.load_from);

println("miss: "*string(Int(cobj.MISS.count)));
println("hit: "*string(Int(cobj.HIT.count)));
println("load: "*string(Int(cobj.LOAD.count)));
println("store: "*string(Int(cobj.STORE.count)));


text:
# Types and Structs
"""
# Description

An abstract type which, at the moment, has only the abstract `AbstractParametersBinary` as subtype.
Thus, it works only with binary (dichotomous) responses. The available mutables for the latter are `Parameters1PL`, `Parameters2PL`, `Parameters3PL`.
They contain the details about the item parameters under the 1-parameter logistic (1PL), 2-parameters logistic (2PL) and 3-parameter logistic (3PL) IRT models, respectively.

# Example

The object `Parameters2PL` has the fields difficulty `b::Float64` and discrimination `a::Float64`. It is possible to define the Bayesian priors and posterior by assigning multivariate distributions from the package `Distributions` to the fields `prior::Distributions.MultivariateDistribution` and `posterior::Distributions.MultivariateDistribution`.
"""
abstract type AbstractParameters end

include("binary/binary.jl")

# Assign prior

"""
```julia
add_prior!(
    parameters::AbstractParameters,
    prior::Distributions.MultivariateDistribution,
)
```

# Description

It assigns the prior `prior` to a `AbstractParameters` instance.

# Arguments

- **`parameters::AbstractParameters`** : Required. Any type of parameters object.
- **`prior::Distributions.Distribution`** : Required. A <n>-variate probability distribution where <n> > 1 and is the numebr of item parameters in `parameters`.

# Examples

    parameters2PL = Parameters2PL()
    bivariate_normal = Distributions.MultivariateNormal([0,0], LinearAlgebra.I(2))
    _add_prior!(parameters2PL, bivariate_normal)
"""
function add_prior!(
    parameters::AbstractParameters,
    prior::Distributions.MultivariateDistribution,
)
    parameters.prior .= prior
end

"""
```julia
add_prior!(
    parameters::AbstractParameters,
    priors::Vector{Distributions.UnivariateDistribution},
)
```

# Description

It transforms the vector `priors` of univariate distributions to their products and assign it to `AbstractParameters` instance.

# Arguments

- **`parameters::AbstractParameters`** : Required. Any type of parameters object.
- **`priors::Vector{Distributions.UnivariateDistribution}`** : Required. A vector of probability distributions. The size of the vector must be the same as the number of item parameters.

# Examples

```@example
parameters2PL = Parameters2PL()
a_dist = Distributions.Normal(0,1)
b_dist = Distributions.Normal(0,1)
_add_prior!(parameters2PL, [a_dist, b_dist])
```
"""
function add_prior!(
    parameters::AbstractParameters,
    priors::Vector{Distributions.UnivariateDistribution},
)
    parameters.prior .= Distributions.Product(priors)
end

# Assign posterior distribution

"""
```julia
_add_posterior!(
    parameters::AbstractParameters,
    posterior::Distributions.Distribution,
)
```

# Description

It assigns the <n>-variate `posterior` distribution to a `AbstractParameters` instance with <n> parameters.

# Arguments

- **`parameters::AbstractParameters`** : Required. Any type of parameters object.
- **`posterior::Distributions.Distribution`** : Required. A <n>-variate probability distribution where <n> > 1 and is the numebr of item parameters in `parameters`.

# Examples

```@example
parameters2PL = Parameters2PL()
bivariate_normal = Distributions.MultivariateNormal([0,0], LinearAlgebra.I(2))
_add_posterior!(parameters2PL, bivariate_normal)
```
"""
function _add_posterior!(
    parameters::AbstractParameters,
    posterior::Distributions.Distribution,
)
    parameters.posterior .= posterior
end

"""
```julia
_add_posterior!(
    parameters::AbstractParameters,
    posteriors::Vector{Distributions.UnivariateDistribution},
)
```

# Description

It transforms the vector `posteriors` of univariate distributions to their products and assign it to `AbstractParameters` instance.

# Arguments

- **`parameters::AbstractParameters`** : Required. Any type of parameters object.
- **`posteriors::Vector{Distributions.UnivariateDistribution}`** : Required. A vector of probability distributions. The size of the vector must be the same as the number of `parameters`.

# Examples

```@example
parameters2PL = Parameters2PL()
a_dist = Distributions.Normal(0,1)
b_dist = Distributions.Normal(0,1)
_add_posterior!(parameters2PL, [a_dist, b_dist])
```
"""
function _add_posterior!(
    parameters::AbstractParameters,
    posteriors::Vector{Distributions.UnivariateDistribution},
)
    parameters.posterior .= Distributions.Product(posteriors)
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Include utilities
include("interpolation.jl")
include("l2projection.jl")
include("basis_lobatto_legendre.jl")

"""
    DGSEM(; RealT=Float64, polydeg::Integer,
            surface_flux=flux_central,
            surface_integral=SurfaceIntegralWeakForm(surface_flux),
            volume_integral=VolumeIntegralWeakForm(),
            mortar=MortarL2(basis))

Create a discontinuous Galerkin spectral element method (DGSEM) using a
[`LobattoLegendreBasis`](@ref) with polynomials of degree `polydeg`.
"""
const DGSEM = DG{Basis} where {Basis<:LobattoLegendreBasis}

# TODO: Deprecated in v0.3 (no longer documented)
function DGSEM(basis::LobattoLegendreBasis,
               surface_flux=flux_central,
               volume_integral=VolumeIntegralWeakForm(),
               mortar=MortarL2(basis))

  surface_integral = SurfaceIntegralWeakForm(surface_flux)
  return DG{typeof(basis), typeof(mortar), typeof(surface_integral), typeof(volume_integral)}(
    basis, mortar, surface_integral, volume_integral)
end

# TODO: Deprecated in v0.3 (no longer documented)
function DGSEM(basis::LobattoLegendreBasis,
               surface_integral::AbstractSurfaceIntegral,
               volume_integral=VolumeIntegralWeakForm(),
               mortar=MortarL2(basis))

  return DG{typeof(basis), typeof(mortar), typeof(surface_integral), typeof(volume_integral)}(
    basis, mortar, surface_integral, volume_integral)
end

# TODO: Deprecated in v0.3 (no longer documented)
function DGSEM(RealT, polydeg::Integer,
               surface_flux=flux_central,
               volume_integral=VolumeIntegralWeakForm(),
               mortar=MortarL2(LobattoLegendreBasis(RealT, polydeg)))
  basis = LobattoLegendreBasis(RealT, polydeg)

  return DGSEM(basis, surface_flux, volume_integral, mortar)
end

DGSEM(polydeg, surface_flux=flux_central, volume_integral=VolumeIntegralWeakForm()) = DGSEM(Float64, polydeg, surface_flux, volume_integral)

# The constructor using only keyword arguments is convenient for elixirs since
# it allows to modify the polynomial degree and other parameters via
# `trixi_include`.
function DGSEM(; RealT=Float64,
                 polydeg::Integer,
                 surface_flux=flux_central,
                 surface_integral=SurfaceIntegralWeakForm(surface_flux),
                 volume_integral=VolumeIntegralWeakForm())
  basis = LobattoLegendreBasis(RealT, polydeg)
  return DGSEM(basis, surface_integral, volume_integral)
end

@inline polydeg(dg::DGSEM) = polydeg(dg.basis)

Base.summary(io::IO, dg::DGSEM) = print(io, "DGSEM(polydeg=$(polydeg(dg)))")

end # @muladd


text:
# Autogenerated wrapper script for Objconv_jll for i686-linux-musl
export objconv

JLLWrappers.@generate_wrapper_header("Objconv")
JLLWrappers.@declare_executable_product(objconv)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        objconv,
        "bin/objconv",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@info "Generating Input Configuration Descriptor Table"
function create_md()
    descriptor = PowerSystems._read_config_file(joinpath(
        dirname(pathof(PowerSystems)),
        "descriptors",
        "power_system_inputs.json",
    ))

    columns = [
        "name",
        "description",
        "unit",
        "unit_system",
        "base_reference",
        "default_value",
        "value_options",
        "value_range",
    ]
    header = "| " * join(columns, " | ") * " |\n" * repeat("|----", length(columns)) * "|\n"

    s = "## [`PowerSystemTableData` Accepted CSV Columns](@id tabledata_input_config) \n\n"
    s = string(
        s,
        "The following tables describe default CSV column definitions accepted by the ",
    )
    s = string(
        s,
        "`PowerSystemeTableData` parser defined by `src/descriptors/power_system_inputs.json`:\n\n"
    )
    for (cat, items) in descriptor
        csv = ""
        for name in PowerSystems.INPUT_CATEGORY_NAMES
            if name[2] == cat
                csv = name[1]
                break
            end
        end

        csv == "" && continue

        s = string(s, "### $csv.csv:\n\n")
        s = string(s, header)
        for item in items
            extra_cols = setdiff(keys(item), columns)
            if !isempty(extra_cols)
                # make sure that there arent unexpected entries
                throw(@error "config file fields not included in header" extra_cols)
            end
            row = []
            for col in columns
                val = string(get(item, col, " "))
                if col == "default_value" && val == " "
                    val = "*REQUIRED*"
                end
                push!(row, val)
            end
            s = string(s, "|" * join(row, "|") * "|\n")
        end
        s = string(s, "\n")
    end
    s = replace(s, r"[_$]" => s"\\\g<0>")
    return s
end

txt = create_md()

open("/Users/cbarrows/Documents/repos/PowerSystems.jl/docs/src/modeler_guide/markdown.txt", "w") do f
      write(f, txt)
end