text:
"""
    abstract type AbstractDirectionType

Abstract type for the direction changes in the sampling region.
"""
abstract type AbstractDirectionType end
export AbstractDirectionType

"""
    struct BAT.FlipDirection <: AbstractDirectionType

Flip direction inverts the direction by multiplying the amplitude with -1.
"""
struct FlipDirection <: AbstractDirectionType end
export FlipDirection

"""
    struct BAT.UpdateDirection <: AbstractDirectionType

Updates the direction in the multidimensional space and accepts it with
Metropolis probability.
"""
struct UpdateDirection <: AbstractDirectionType end
export UpdateDirection

"""
    struct ECMCSampler <: AbstractSamplingAlgorithm

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```
"""
@with_kw struct ECMCSampler{TR<:AbstractDensityTransformTarget} <: AbstractSamplingAlgorithm
    trafo::TR = NoDensityTransform()
    nsamples::Int = 10^4
    chain_length::Int = 50
    delta::Float64 = 0.5
    remaining_jumps_before_refresh::Int = 5
    direction_change::AbstractDirectionType = FlipDirection()
end
export ECMCSampler

function bat_sample_impl(
    rng::AbstractRNG,
    target::AnyDensityLike,
    algorithm::ECMCSampler
)
    density_notrafo = convert(AbstractDensity, target)
    shaped_density, trafo = bat_transform(algorithm.trafo, density_notrafo)
    shape = varshape(shaped_density)
    density = unshaped(shaped_density)

    samples = _gen_samples(density, algorithm)

    logvals = map(logdensityof(density), samples)
    weights = ones(length(samples))

    samples_trafo = shape.(DensitySampleVector(samples, logvals, weight = weights))
    samples_notrafo = inv(trafo).(samples_trafo)

    return (result = samples_notrafo, result_trafo = samples_trafo, trafo = trafo, integral = nothing)
end

function _gen_samples(density::AbstractDensity,
    algorithm::ECMCSampler
)
    initial_sample = vec(collect(reduce(vcat, values(rand(BAT.getprior(density).dist)))))
    samples = [initial_sample]

    for i in 1:algorithm.nsamples
        push!(samples, _run_ecmc(density, algorithm, samples[end]))
    end
    return samples
end

function _run_ecmc(density::AbstractDensity,
    algorithm::ECMCSampler,
    C::Vector{Float64};
    moving_dimension::Int = rand(1:totalndof(density)),
    jump_amplitude::Float64 = _jump_amplitude(algorithm)
)
    dim = totalndof(density)
    C = [values(C)[i] for i in 1:length(C)]
    remaining_jumps_before_sample = algorithm.chain_length
    remaining_jumps_before_refresh = algorithm.remaining_jumps_before_refresh

    while remaining_jumps_before_sample >= 0
        remaining_jumps_before_sample -= 1
        remaining_jumps_before_refresh -= 1

        id_matrix = Matrix{Float64}(I, dim, dim)
        new_C = C .+ id_matrix[moving_dimension, :] .* jump_amplitude

        p_acceptance = exp.(-_energy_difference(density, new_C, C))

        rand(Uniform(0.0, 1.0)) <= p_acceptance[1] ? C = new_C : jump_amplitude = _change_direction(algorithm.direction_change, jump_amplitude, C, new_C, algorithm, density)

        if remaining_jumps_before_refresh <= 0
            remaining_jumps_before_refresh = algorithm.remaining_jumps_before_refresh
            moving_dimension = rand(1:dim)
            jump_amplitude = _jump_amplitude(algorithm)
        end
    end
    return (samples = C)
end

_energy_function(density::AbstractDensity, x) = -map(BAT.logvalof(density), [[values(x)[i] for i in 1:length(x)]])

_energy_difference(density::AbstractDensity, new_x::Vector{Float64}, old_x::Vector{Float64}) = _energy_function(density, new_x) .- _energy_function(density, old_x)

_jump_amplitude(algorithm::ECMCSampler) = rand([-1.0, 1.0]) * algorithm.delta * rand(Uniform(0.5, 1.5))

function _change_direction(direction_type::FlipDirection,
    jump_amplitude::Float64,
    C::Vector{Float64},
    new_C::Vector{Float64},
    algorithm::ECMCSampler,
    density::AbstractDensity
)
    return jump_amplitude * (-1)
end

function _change_direction(direction_type::UpdateDirection,
    jump_amplitude::Float64,
    C::Vector{Float64},
    new_C::Vector{Float64},
    algorithm::ECMCSampler,
    density::AbstractDensity
)
    old_vector = jump_amplitude
    new_vector = _jump_amplitude(algorithm)

    prob_current_C = 1.0 .- exp.(-_energy_difference(density, C .+ old_vector .* algorithm.delta, new_C))
    prob_new_C = 1.0 .- exp.(-_energy_difference(density, C .- new_vector .* algorithm.delta, new_C))

    acceptance_probability = minimum([1.0, prob_new_C[1] / prob_current_C[1]])
    rand(Uniform(0, 1)) <= acceptance_probability ? (return new_vector) : (return -old_vector)
end


text:
function write_xsf( filnam, LL_in, atpos; molecule=false )
    #
    f = open(filnam, "w")
    Natoms = size(atpos)[2]
    #
    if molecule
        @printf(f, "MOLECULE\n")
    else
        @printf(f, "CRYSTAL\n")
    end
    LL = LL_in'
    @printf(f, "PRIMVEC\n")
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[1,1], LL[1,2], LL[1,3])
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[2,1], LL[2,2], LL[2,3])
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[3,1], LL[3,2], LL[3,3])
    @printf(f, "PRIMCOORD\n")
    @printf(f, "%8d %8d\n", Natoms, 1)
    for ia = 1:Natoms
        @printf(f, "H  %18.10f %18.10f %18.10f\n", atpos[1,ia], atpos[2,ia], atpos[3,ia])
    end

    close(f)
end

function write_xsf_3d_crystal( filnam, Ns, LL_in, data3d; center=zeros(3) )
    #
    f = open(filnam, "a")  # FIXME: What if filnam is not exist?
    #
    LL = LL_in'
    @printf(f, "BEGIN_BLOCK_DATAGRID_3D\n")
    @printf(f, "made_by_ffr\n")
    @printf(f, "DATAGRID_3D_UNKNOWN\n")
    @printf(f, "%8d %8d %8d\n", Ns[1]+1, Ns[2]+1, Ns[3]+1 )
    @printf(f, "%18.10f %18.10f %18.10f\n", center[1], center[2], center[3])
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[1,1], LL[1,2], LL[1,3])
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[2,1], LL[2,2], LL[2,3])
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[3,1], LL[3,2], LL[3,3])
    #
    rDat3d = reshape( data3d, (Ns[1],Ns[2],Ns[3]) )
    for k = 1:Ns[3]+1
        for j = 1:Ns[2]+1
            for i = 1:Ns[1]+1
                ii = ( i == Ns[1]+1 ? 1 : i )
                jj = ( j == Ns[2]+1 ? 1 : j )
                kk = ( k == Ns[3]+1 ? 1 : k )
                #
                @printf(f, "%18.10f\n", rDat3d[ii,jj,kk])
            end
        end
    end
    @printf(f, "END_DATAGRID_3D\n")
    @printf(f, "END_BLOCK_DATAGRID_3D\n")
    close(f)
end

ext: jl


text:
## Code for optimizers adapted from
## https://github.com/jacobcvt12/GradDescent.jl/
import ChainRulesCore

"""Should we ignore field with this path when updating the struct?"""
ignored_field(st::S, path::Val) where S = false

_dot_path(x::Symbol) = [x]

function _dot_path(ex::Expr)
    subresult = _dot_path(ex.args[1])
    return [subresult; [ex.args[2].value]]
end

"""
Ignore this field when updating the struct. Example:

mutable struct MyModel
    linear::Linear
end

@ignore MyModel.linear.b
"""
macro ignore(ex)
    @assert Meta.isexpr(ex, :(.))
    full_path = _dot_path(ex)
    T, path = full_path[1], (full_path[2:end]...,)
    return esc(:($ignored_field(st::$T, path::Val{$path}) = true))
end

abstract type Optimizer end

function Yota.update!(opt::Optimizer, m, gm; ignore=Set())
    # we use Zero() to designate values that need not to be updated
    gm isa ChainRulesCore.ZeroTangent && return
    for (path, gx) in Yota.path_value_pairs(gm)
        if !in(path, ignore) && !ignored_field(m, Val(path)) && !isa(gx, ChainRulesCore.ZeroTangent)
            x_t0 = Yota.getfield_nested(m, path)
            x_t1 = make_update!(opt, path, x_t0, gx)
            Yota.setfield_nested!(m, path, x_t1)
        end
    end
    opt.t += 1
end

function Yota.update!(opt::Optimizer, x::AbstractArray, gx; ignore=Set())
    x .= make_update!(opt, (), x, gx)
    opt.t += 1
end

################################################################################
#                                    SGD                                       #
################################################################################

mutable struct SGD <: Optimizer
    t::Int
    lr::Float32
    momentum::Float32
    v_t::Dict{Any, Any}   # path => previous velocity
end

function SGD(lr; momentum=0)
    # @assert momentum >= 0.0 "momentum must be >= 0"
    SGD(0, lr, momentum, Dict())
end

Base.show(io::IO, opt::SGD) = print(io, "SGD(lr=$(opt.lr), momentum=$(opt.momentum))")

function make_update!(opt::SGD, path, x, gx)
    v_t = get(opt.v_t, path, zero(gx))
    v_t1 = opt.momentum .* v_t .+ opt.lr .* gx
    opt.v_t[path] = v_t1
    x_t = x
    x_t1 = x_t .- v_t1
    return x_t1
end

################################################################################
#                                RMSprop                                       #
################################################################################

mutable struct RMSprop <: Optimizer
    t::Int32
    eps::Float32
    eta::Float32
    gamma::Float32
    E_g_sq_t::Dict{Any, Any}
end

function RMSprop(; eta::Real=0.001, gamma::Real=0.01, eps::Real=1e-8)
    @assert eta > 0.0 "eta must be greater than 0"
    @assert gamma > 0.0 "gamma must be greater than 0"
    @assert eps > 0.0 "eps must be greater than 0"

    RMSprop(0, eps, eta, gamma, Dict())
end

function make_update!(opt::RMSprop, path, x, gx)
    # resize accumulated and squared updates
    E_g_sq_t = get(opt.E_g_sq_t, path, zero(gx))
    # accumulate gradient
    E_g_sq_t = opt.gamma * E_g_sq_t + (1f0 - opt.gamma) * (gx .^ 2)
    # compute update
    RMS_g_t = sqrt.(E_g_sq_t .+ opt.eps)
    opt.E_g_sq_t[path] = E_g_sq_t
    return x - opt.eta .* gx ./ RMS_g_t
end

################################################################################
#                                   Adam                                       #
################################################################################

mutable struct Adam <: Optimizer
    t::Int32
    eps::Float32
    lr::Float32
    beta1::Float32
    beta2::Float32
    m_t::Dict{Any, Any}
    v_t::Dict{Any, Any}
end

function Adam(;lr::Real=0.001, beta1::Real=0.9, beta2::Real=0.999, eps::Real=10e-8)
    @assert lr > 0.0 "lr must be greater than 0"
    @assert beta1 > 0.0 "beta1 must be greater than 0"
    @assert beta2 > 0.0 "beta2 must be greater than 0"
    @assert eps > 0.0 "eps must be greater than 0"
    Adam(1, eps, lr, beta1, beta2, Dict(), Dict())
end

function make_update!(opt::Adam, path, x, gx)
    # no update for symbolic Zero
    gx isa ChainRulesCore.ZeroTangent && return x
    # resize biased moment estimates if first iteration
    m_t = get(opt.m_t, path, zero(gx))
    v_t = get(opt.v_t, path, zero(gx))
    # update biased first moment estimate
    m_t = opt.beta1 * m_t + (1f0 - opt.beta1) * gx
    # update biased second raw moment estimate
    v_t = opt.beta2 * v_t + (1f0 - opt.beta2) * (gx .^ 2)
    # save moments for the next iteration
    opt.m_t[path] = m_t
    opt.v_t[path] = v_t
    # compute bias corrected first moment estimate
    m̂_t = m_t ./ (1f0 - opt.beta1^opt.t)
    # compute bias corrected second raw moment estimate
    v̂_t = v_t ./ (1f0 - opt.beta2^opt.t)
    # apply update
    x_t1 = x - opt.lr * m̂_t ./ (sqrt.(v̂_t) .+ opt.eps)
    return x_t1
end


text:
@testset "domainmodel" begin
    @testset "verify_ipfshash" begin
        # Should fail because not long enough
        hash = "Qmaaa"
        @test !verify_ipfshash(hash)

        # Should fail because doesn't start with "Qm"
        hash = "AmauYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk"
        @test !verify_ipfshash(hash)

        # Should pass
        hash = "QmauYgPmss6CEZXtaRvvGW2oiyLqxpoCkWNCTmFPVTFDfk"
        @test verify_ipfshash(hash)
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/julia-beginner-from-python-numba-outperforms-julia-in-rewrite-any-tips-to-improve-performance/66414

# %%
using BenchmarkTools
#using Distributed
#using Plots
#using Profile
#using ProfileVega
#using PyCall
using ScikitLearn
#using Traceur
using TimerOutputs

# %%
@sk_import datasets: (fetch_covtype)

# %%
digits_data = fetch_covtype();

# %%
X_digits = abs.(digits_data["data"]);
X_digits = transpose(X_digits);

# %%
d, n = size(X_digits)

# %%
struct NaiveGreedy
    X::Matrix{Float64}
end

struct Result
    ranking::Vector{Int32}
    gains::Vector{Float64}
end

# %%
function lexsort(a, b, rev=false)
    idxs = sortperm(a, alg=MergeSort, rev=rev)
    return idxs[sortperm(b[idxs], alg=MergeSort, rev=rev)]
end

function get_gains!(X, current_values, idxs, gains)
    @inbounds Threads.@threads for i in eachindex(idxs)
        s = 0.0
        for j in eachindex(current_values)
            s += @fastmath sqrt(current_values[j] + X[j, idxs[i]])
        end
        gains[i] = s
    end
end

function calculate_gains!(X, gains, current_values, idxs, current_concave_value_sum)
    get_gains!(X, current_values, idxs, gains)

    gains .-= current_concave_value_sum
    return gains
end

# %%
function fit(optimizer::NaiveGreedy, k, sample_cost)
    @timeit to "intro" begin
        d, n = size(optimizer.X)

        cost = 0.0

        ranking = Int32[]
        total_gains = Float64[]

        mask = zeros(Int8, n)
        current_values = zeros(Float64, d)
        current_concave_values = sqrt.(current_values)
        current_concave_values_sum = sum(current_concave_values)

        idxs = 1:n
    end

    gains = zeros(Float64, size(idxs)[1])
    @timeit to "while loop" begin
        while cost < k
            gains = @timeit to "calc_gains" calculate_gains!(optimizer.X, gains, current_values, idxs, current_concave_values_sum)

            if sample_cost != nothing
                gains ./= sample_cost[idxs]
                idx_idxs = @timeit to "lexsort" lexsort(gains, 1:size(gains)[1], true)

                @timeit to "select_idx" begin
                    for i in 1:size(idx_idxs)[1]
                        global idx = idx_idxs[i]
                        global best_idx = idxs[idx]
                        if cost + sample_cost[best_idx] <= k
                            break
                        end
                    end
                end
                curr_cost = sample_cost[best_idx]
            else
                global idx = argmax(gains)
                global best_idx = idxs[idx]
                curr_cost = 1.
            end

            if cost + curr_cost > k
                break
            end

            @timeit to "select_idx" begin
                cost += curr_cost
                # Calculate gains
                gain = gains[idx] * curr_cost
            end

            @timeit to "select_next" begin
                # Select next
                current_values += view(optimizer.X, :, best_idx)
                current_concave_values .= sqrt.(current_values)
                current_concave_values_sum = sum(current_concave_values)

                push!(ranking, best_idx)
                push!(total_gains, gain)

                mask[best_idx] = 1
                idxs = findall(==(0), mask)
            end
        end
    end
    return Result(ranking, total_gains)
end

# %%
const to = TimerOutput()
k = 1000
opt1 = NaiveGreedy(X_digits)
res1 = fit(opt1, k, nothing)
to

# %%
function fit(optimizer::NaiveGreedy, k, sample_cost)
    @timeit to "intro" begin
        d, n = size(optimizer.X)

        cost = 0.0

        ranking = Int32[]
        total_gains = Float64[]

        #mask = zeros(Int8, n)
        current_values = zeros(Float64, d)
        current_concave_values = sqrt.(current_values)
        current_concave_values_sum = sum(current_concave_values)

        #idxs = 1:n
        idxs = collect(1:n)
    end

    gains = zeros(Float64, size(idxs)[1])
    @timeit to "while loop" begin
        while cost < k
            gains = @timeit to "calc_gains" calculate_gains!(optimizer.X, gains, current_values, idxs, current_concave_values_sum)

            if sample_cost != nothing
                gains ./= sample_cost[idxs]
                idx_idxs = @timeit to "lexsort" lexsort(gains, 1:size(gains)[1], true)

                @timeit to "select_idx" begin
                    for i in 1:size(idx_idxs)[1]
                        global idx = idx_idxs[i]
                        global best_idx = idxs[idx]
                        if cost + sample_cost[best_idx] <= k
                            break
                        end
                    end
                end
                curr_cost = sample_cost[best_idx]
            else
                global idx = argmax(gains)
                global best_idx = idxs[idx]
                curr_cost = 1.
            end

            if cost + curr_cost > k
                break
            end

            @timeit to "select_idx" begin
                cost += curr_cost
                # Calculate gains
                gain = gains[idx] * curr_cost
            end

            @timeit to "select_next" begin
                # Select next
                #current_values += view(optimizer.X, :, best_idx)
                #current_concave_values .= sqrt.(current_values)
                #current_concave_values_sum = sum(current_concave_values)
                current_values .+= view(optimizer.X, :, best_idx)
                current_concave_values_sum = sum(sqrt, current_values)

                push!(ranking, best_idx)
                push!(total_gains, gain)

                #mask[best_idx] = 1
                #idxs = findall(==(0), mask)
                popat!(idxs, findfirst(==(best_idx), idxs))
            end
        end
    end
    return Result(ranking, total_gains)
end

# %%
const to = TimerOutput()
k = 1000
opt1rev = NaiveGreedy(X_digits)
res1rev = fit(opt1rev, k, nothing)
to

# %%
@show res1rev.ranking == res1.ranking
@show res1rev.gains ≈ res1.gains;

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
"""
    struct SymSparseMatrixCSR{T,Ti<:Integer} <: AbstractSparseMatrix{T,Ti}

Matrix type for storing symmetric sparse matrices in the
Compressed Sparse Row format. The standard way
of constructing SparseMatrixCSR is through the
[`symsparsecsr`](@ref) function.
"""
struct SymSparseMatrixCSR{Bi,T,Ti<:Integer} <: AbstractSparseMatrix{T,Ti}
    uppertrian :: SparseMatrixCSR{Bi,T,Ti}
end

function show(io::IO, ::MIME"text/plain", S::SymSparseMatrixCSR)
    xnnz = nnz(S)
    print(io, S.uppertrian.m, "×", S.uppertrian.n, " ",
              typeof(S), " with ", xnnz, " stored ",
              xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        print(io, ":")
        show(IOContext(io, :typeinfo => eltype(S)), S)
    end
end
show(io::IO, S::SymSparseMatrixCSR) = show(io, S.uppertrian)

size(A::SymSparseMatrixCSR) = size(A.uppertrian)
getindex(A::SymSparseMatrixCSR, x::Integer, y::Integer) = getindex(A.uppertrian,min(x,y),max(x,y))

"""
    nnz(S::SymSparseMatrixCSR)

Returns the number of stored (filled) elements in a sparse array.
"""
nnz(S::SymSparseMatrixCSR) = nnz(S.uppertrian)

"""
    count(pred, S::SymSparseMatrixCSR) -> Integer

Count the number of elements in nonzeros(S) for which predicate pred returns true.
"""
count(pred, S::SymSparseMatrixCSR) = count(pred, S.uppertrian)

"""
    nonzeros(S::SymSparseMatrixCSR)

Return a vector of the structural nonzero values in sparse array S.
This includes zeros that are explicitly stored in the sparse array.
The returned vector points directly to the internal nonzero storage of A,
and any modifications to the returned vector will mutate A as well.
"""
nonzeros(S::SymSparseMatrixCSR) = nonzeros(S.uppertrian)

"""
    nzrange(S::SymSparseMatrixCSR, row::Integer)

Return the range of indices to the structural nonzero values of a
sparse matrix row.
"""
nzrange(S::SymSparseMatrixCSR, row::Integer) = nzrange(S.uppertrian, row)

"""
    findnz(S::SymSparseMatrixCSR)

Return a tuple (I, J, V) where I and J are the row and column indices
of the stored ("structurally non-zero") values in sparse matrix A,
and V is a vector of the values.
"""
findnz(S::SymSparseMatrixCSR) = findnz(S.uppertrian)

"""
    rowvals(S::SymSparseMatrixCSR)

Return an error.
CSR sparse matrices does not contain raw row values.
It contains row pointers instead that can be accessed
by using [`nzrange`](@ref).
"""

rowvals(S::SymSparseMatrixCSR) = rowvals(S.uppertrian)

"""
    colvals(S::SparseMatrixCSR)

Return a vector of the col indices of S.
Any modifications to the returned vector will mutate S as well.
Providing access to how the co,l indices are stored internally
can be useful in conjunction with iterating over structural
nonzero values. See also [`nonzeros`](@ref) and [`nzrange`](@ref).
"""

colvals(S::SymSparseMatrixCSR) = colvals(S.uppertrian)

"""
    symsparsecsr(I, J, V, [m, n, combine])

Create a Symmetric sparse matrix S of dimensions m x n
 such that S[I[k], J[k]] = V[k], and m=n.
The combine function is used to combine duplicates.
If m and n are not specified, they are set to
maximum(I) and maximum(J) respectively.
If the combine function is not supplied, combine defaults to +.
All elements of I must satisfy 1 <= I[k] <= m,
and all elements of J must satisfy 1 <= J[k] <= n.
Numerical zeros in (I, J, V) are retained as structural nonzeros.
"""
function symsparsecsr(T::Type{SymSparseMatrixCSR{Bi,Tv,Ti}},I::Vector{Ti},J::Vector{Ti},V::Vector{Tv},args...)  where {Bi,Tv,Ti}
    m = length(args)>0 ? args[1] : isempty(I) ? 0 : Int(maximum(I))
    n = length(args)>1 ? args[2] : m
    c = length(args)>2 ? args[3] : +
    symsparsecsr(SymSparseMatrixCSR{Bi},I,J,V,m,n,c)
end

function symsparsecsr(T::Type{SymSparseMatrixCSR{Bi}},I,J,V,args...)  where {Bi}
    m = length(args)>0 ? args[1] : isempty(I) ? 0 : Int(maximum(I))
    n = length(args)>1 ? args[2] : m
    c = length(args)>2 ? args[3] : +
    symsparsecsr(T,I,J,V,m,n,c)
end

function symsparsecsr(::Type{SymSparseMatrixCSR{Bi}},I,J,V,m,n,combine) where {Bi}
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are ($m, $n)"))
    SymSparseMatrixCSR(sparsecsr(SparseMatrixCSR{Bi},I,J,V,m,n,combine))
end

symsparsecsr(::Type{SymSparseMatrixCSR}, I,J,V,args...) =
    symsparsecsr(SymSparseMatrixCSR{1},I,J,V,args...)

symsparsecsr(I,J,V,args...) =
    symsparsecsr(SymSparseMatrixCSR,I,J,V,args...)

"""
    function push_coo!(::Type{SymSparseMatrixCSR},I,J,V,ik,jk,vk)

Inserts entries in COO vectors for further building a SymSparseMatrixCSR.
It stores only the upper triangle, ignoring entries with (ik>jk) coordinates.
"""
function push_coo!(::Type{SymSparseMatrixCSR},
        I::Vector,J::Vector,V::Vector,ik::Integer,jk::Integer,vk::Number)
    (ik>jk) && return
    (push!(I, ik), push!(J, jk), push!(V, vk))
end

push_coo!(::Type{SymSparseMatrixCSR{Bi}}, I, J, V, ik, jk, vk) where {Bi} =
    push_coo!(SymSparseMatrixCSR, I, J, V, ik, jk, vk)

push_coo!(::Type{SymSparseMatrixCSR{Bi,Tv,Ti}}, I, J, V, ik, jk, vk) where {Bi,Tv,Ti} =
    push_coo!(SymSparseMatrixCSR, I, J, V, ik, jk, vk)

"""
    function finalize_coo!(::Type{SymSparseMatrixCSR},I,J,V,m,n)

Finalize COO arrays for building a SymSparseMatrixCSR.
Check and insert diagonal entries in COO vectors if needed.
"""
function finalize_coo!(T::Type{SymSparseMatrixCSR},
        I::Vector,J::Vector,V::Vector, m::Integer, n::Integer)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are ($m, $n)"))
    touched = zeros(Bool,m)
    for k in 1:length(I)
        Ik = I[k]
        Jk = J[k]
        if Ik == Jk
            touched[Ik] = true
        end
    end
    for k in 1:m
        if ! touched[k]
            push_coo!(T,I,J,V,k,k,zero(eltype(V)))
        end
    end
end

finalize_coo!(T::Type{SymSparseMatrixCSR{Bi}}, I , J, V, m, n)  where {Bi} =
    finalize_coo!(SymSparseMatrixCSR, I, J, V, m, n)

finalize_coo!(T::Type{SymSparseMatrixCSR{Bi,Tv,Ti}}, I , J, V, m, n)  where {Bi,Tv,Ti} =
    finalize_coo!(SymSparseMatrixCSR, I, J, V, m, n)

"""
    function mul!(y::AbstractVector,A::SymSparseMatrixCSR,v::AbstractVector{T}) where {T}

Calculates the matrix-vector product ``Av`` and stores the result in `y`,
overwriting the existing value of `y`.
"""
function mul!(y::AbstractVector,A::SymSparseMatrixCSR,v::AbstractVector{T}) where {T}
    A.uppertrian.n == size(v, 1) || throw(DimensionMismatch())
    A.uppertrian.m == size(y, 1) || throw(DimensionMismatch())

    y .= zero(T)
    for row = 1:size(y, 1)
        @inbounds for nz in nzrange(A,row)
            col = A.uppertrian.colval[nz]-A.uppertrian.offset
            y[row] += A.uppertrian.nzval[nz]*v[col]
            row != col && (y[col] += A.uppertrian.nzval[nz]*v[row])
        end
    end
    return y
end

*(A::SymSparseMatrixCSR, v::Vector) = (y = similar(v,A.uppertrian.n);mul!(y,A,v))

"""
    function hasrowmajororder(::Type{SymSparseMatrixCSR})

Check if values are stored in row-major order.
Return true.
"""
hasrowmajororder(::Type{SymSparseMatrixCSR}) = true
hasrowmajororder(a::SymSparseMatrixCSR) = hasrowmajororder(SymSparseMatrixCSR)

"""
    function hascolmajororder(::Type{SymSparseMatrixCSR})

Check if values are stored in col-major order.
Return false.
"""
hascolmajororder(::Type{SymSparseMatrixCSR}) = false
hascolmajororder(a::SymSparseMatrixCSR) = hascolmajororder(SymSparseMatrixCSR)

"""
    function getptr(S::SymSparseMatrixCSR)

Return rows pointer.
"""
getptr(S::SymSparseMatrixCSR) = getptr(S.uppertrian)

"""
    function getindices(S::SymSparseMatrixCSR)

Return column indices.
"""
getindices(S::SymSparseMatrixCSR) = colvals(S)

"""
    function convert(::Type{SymSparseMatrixCSR}, x::SymSparseMatrixCSR)

Convert x to a value of type SymSparseMatrixCSR.
"""
convert(::Type{SymSparseMatrixCSR}, x::SymSparseMatrixCSR)  = convert(SymSparseMatrixCSR{1}, x)

function convert(::Type{SymSparseMatrixCSR{Bi}}, x::SymSparseMatrixCSR{xBi}) where {Bi,xBi}
    if Bi == xBi
        return x
    else
        return SymSparseMatrixCSR(convert(SparseMatrixCSR{Bi}, x.uppertrian))
    end
end

function convert(::Type{SymSparseMatrixCSR{Bi,Tv,Ti}}, x::SymSparseMatrixCSR{xBi,xTv,xTi}) where {Bi,Tv,Ti,xBi,xTv,xTi}
    if (Bi,Tv,Ti) == (xBi,xTv,xTi)
        return x
    else
        return SymSparseMatrixCSR(convert(SparseMatrixCSR{Bi,Tv,Ti}, x.uppertrian))
    end
end


text:
module BilinearCheck

using Random

import Random: rand

export randComplex,inner,isLeftLinear,isRightLinear,isHomogenousBi,
       isHomogenousSus,nonDegenFirst,nonDegenSecond,nonDegen,
       isSymmetric,isHermitian,isBilinear,isSusquilinear

function randComplex(max)
    r = rand(Int64,3,1) .% max
    i = rand(Int64,3,1) .% max
    return r .+ i.*im
end

inner(v,w) = (w'*v)[1]

function isLeftLinear(β,a,b,c)
    β(a+b,c) == β(a,c) + β(b,c)
end

function isRightLinear(β,a,b,c)
    β(c,a+b) == β(c,a) + β(c,b)
end

function isHomogenousBi(β,a,b,λ)
    (β(λ*a,b) == λ*β(a,b)) && (β(a,λ*b) == λ*β(a,b))
end

function isHomogenousSus(β,a,b,λ)
    (β(λ*a,b) == λ*β(a,b)) && (β(a,λ*b) == conj(λ)*β(a,b))
end

function nonDegenFirst(β,a,b)
    a = zeros(size(a))
    norm(β(a,b)) == 0
end

function nonDegenSecond(β,a,b)
    b = zeros(size(b))
    norm(β(a,b)) == 0
end

function nonDegen(β,a,b)
    return nonDegenFirst(β,a,b) && nonDegenSecond(β,a,b)
end

function isSymmetric(β,a,b)
    β(a,b) == β(b,a)
end

function isHermitian(β,a,b)
    β(a,b) == conj(β(b,a))
end

function isBilinear(β,a,b,c,λ)
    return isLeftLinear(β,a,b,c) && isRightLinear(β,a,b,c) && isHomogenousBi(β,a,b,λ)
end

function isSusquilinear(β,a,b,c,λ)
    return isLeftLinear(β,a,b,c) && isRightLinear(β,a,b,c) && isHomogenousSus(β,a,b,λ)
end

end #module BilinearCheck

max_stars_repo_name: sphanna/BilinearCheck

max_issues_repo_name: sphanna/BilinearCheck

max_forks_repo_name: sphanna/BilinearCheck


text:
# An alternative resampling algorithm described in Fearnhead (2004):
#
# for each particle, generate putatives (z_1:n^i, j)
# sample one successor proportional to p((z_1:n^i, j) | x_1:n+1)
# assign weight based on sum_1^k_i+1(p((z_1:n^i, j) | x_1:n+1)/p(z_1:n^i | x_1:n))

mutable struct ChenLiuParticles{P} <: ParticleFilter
    particles::Vector{P}
    N::Int
    rejuvination_threshold::Float64
end

function Base.show(io::IO, ps::ChenLiuParticles)
    n = length(ps.particles)
    if n < ps.N
        println(io, "Particle filter with $n (up to $(ps.N)) particles:")
    else
        println(io, "Particle filter with $n particles:")
    end
    show(IOContext(io, :compact=>true), ps.particles)
end

# initialize with full population of empty particles, because this method
# doens't benefit from redundancy in the same way as the Fearnhead method does
ChenLiuParticles(n::Int, prior::Union{Tuple,<:Distribution}, stateprior::T; rejuv::Float64=50.) where T<:StatePrior =
    ChenLiuParticles([InfiniteParticle(prior, stateprior) for _ in 1:n], n, rejuv)

function propogate_chenliu(p::P, y) where P<:AbstractParticle
    ps = collect(putatives(p, y))
    # sample next based on updated weights (which are proportional to the
    # posterior p( (z_1:n, j) | x_1:n+1 ) because they've been updated based on
    # the same ancestor, multiplying by
    # p( (z_1:n,j) | x_1:n+1) / p(z_1:n | x_1:n) ∝ p( (z_1:n,j) | x_1:n+1 )
    weights = weight.(ps)
    next = instantiate(wsample(ps, weights))
    # now to update the weight.  we need it to be
    # w_n+1 = w_n × sum( p(z_1:n,j | x_1:n+1) / p(z_1:n | x_1:n) )
    # call the putative weight of particle j v_j.
    # v_j = w_n × p((z_1:n,j) | x_1:n+1) / p(z_1:n | x_1:n), so we just need to
    # add those up
    return weight(next, sum(weights))
end

"""
    coefvar(x)

Calculat the coefficient of variation for x: σ / μ * 100
"""
function coefvar(x)
    μ, σ = StatsBase.mean_and_std(x)
    σ / μ * 100
end

function fit!(ps::ChenLiuParticles, y)
    map!(p->propogate_chenliu(p, y), ps.particles, ps.particles)
    # rejuvinate if necessary
    ws = weight.(ps.particles)
    cv = coefvar(ws)
    if cv > ps.rejuvination_threshold
        @debug "  Rejuvinating ($cv > $(ps.rejuvination_threshold))"
        ps.particles = weight.(wsample(ps.particles, weight.(ps.particles), ps.N, replace=true), 1/ps.N)
    else
        # normalize weights to prevent underflow
        # TODO: use log-weight instead to avoid this altogether
        total_w = sum(ws)
        ps.particles = weight.(ps.particles, ws ./ total_w)
    end
    ps
end


text:
"""微分運動学

"""
module DifferentialKinematics

export J_0
export J_1
export J_2
export J_3
export J_4

const N = 5

function J_0!(q::Vector{Float64}, ξ::Float64, out::Matrix{Float64})
    ccall(
        (:J_0, "o/soft_robot/derivation_of_kinematics/derived/so/J_0.so"),
        Cvoid,
        (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}),
        q, ξ, out
    )
end

"""ヤコビ行列"""
function J_0(q::Vector{Float64}, ξ::Float64, n::Int64)
    Z = Matrix{Float64}(undef, 3*N, 3)
    J_0!(q, ξ, Z)
    return Z'[:, 1:3*n]
end

function J_1!(q::Vector{Float64}, ξ::Float64, out::Matrix{Float64})
    ccall(
        (:J_1, "o/soft_robot/derivation_of_kinematics/derived/so/J_1.so"),
        Cvoid,
        (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}),
        q, ξ, out
    )
end

"""ヤコビ行列"""
function J_1(q::Vector{Float64}, ξ::Float64, n::Int64)
    Z = Matrix{Float64}(undef, 3*N, 3)
    J_1!(q, ξ, Z)
    return Z'[:, 1:3*n]
end

function J_2!(q::Vector{Float64}, ξ::Float64, out::Matrix{Float64})
    ccall(
        (:J_2, "o/soft_robot/derivation_of_kinematics/derived/so/J_2.so"),
        Cvoid,
        (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}),
        q, ξ, out
    )
end

"""ヤコビ行列"""
function J_2(q::Vector{Float64}, ξ::Float64, n::Int64)
    Z = Matrix{Float64}(undef, 3*N, 3)
    J_2!(q, ξ, Z)
    return Z'[:, 1:3*n]
end

function J_3!(q::Vector{Float64}, ξ::Float64, out::Matrix{Float64})
    ccall(
        (:J_3, "o/soft_robot/derivation_of_kinematics/derived/so/J_3.so"),
        Cvoid,
        (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}),
        q, ξ, out
    )
end

"""ヤコビ行列"""
function J_3(q::Vector{Float64}, ξ::Float64, n::Int64)
    Z = Matrix{Float64}(undef, 3*N, 3)
    J_3!(q, ξ, Z)
    return Z'[:, 1:3*n]
end

function J_4!(q::Vector{Float64}, ξ::Float64, out::Matrix{Float64})
    ccall(
        (:J_4, "o/soft_robot/derivation_of_kinematics/derived/so/J_4.so"),
        Cvoid,
        (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}),
        q, ξ, out
    )
end

"""ヤコビ行列"""
function J_4(q::Vector{Float64}, ξ::Float64, n::Int64)
    Z = Matrix{Float64}(undef, 3*N, 3)
    J_4!(q, ξ, Z)
    return Z'[:, 1:3*n]
end

end

# using .DifferentialKinematics

# using LinearAlgebra

# q = [
#     0, 0, 0,
#     0, 0, 0,
#     0.1, 0.2, 0.1
# ]

# q0 = [0.01, 0.05, 0.1]

# J = J_2(q, 1.0, 3)
# J0 = J_0(q0, 1.0, 1)


text:
@testset "lookup" begin
    a = V"a"
    b = V"b"
    c = V"c"
    bindings = Bindings(a, 1,
                        Bindings(b, 2,
                                 EmptyBindings()))
    @test lookupfirst(bindings, a) == (1, true)
    @test lookupfirst(bindings, b) == (2, true)
    @test lookupfirst(bindings, c) == (nothing, false)
end

@testset "ubind" begin
    a = V"a"
    b = V"b"
    c = V"c"
    ubind(a, 1) do bindings
        ubind(b, 2, bindings) do bindings
            ubind(c, 3, bindings) do bindings
                @test lookupfirst(bindings, a) == (1, true)
                @test lookupfirst(bindings, b) == (2, true)
                @test lookupfirst(bindings, c) == (3, true)
            end
        end
    end
    ubind([a=>b]) do bindings
        # Binding two variables together is symetric:
        @test lookupfirst(bindings, a) == (b, true)
        @test lookupfirst(bindings, b) == (a, true)
    end
end

@testset "bindings iteration" begin
    @test Base.IteratorSize(EmptyBindings) == Base.SizeUnknown()
    @test Base.IteratorSize(Bindings) == Base.SizeUnknown()
    @test Base.IteratorEltype(EmptyBindings) == Tuple{AbstractVar, Any}
    @test Base.IteratorEltype(Bindings) == Tuple{AbstractVar, Any}
    ubind([V"a"=>1, V"b"=>2, V"c"=>3]) do bindings
        @test map(p -> p[1], bindings) == [V"a", V"b", V"c"]
        @test map(p -> p[2], bindings) == [1, 2, 3]
    end
end

@testset "lookupall" begin
    a = V"a"
    b = V"b"
    ubind([a=>1, b=>2, a=>3]) do bindings
        @test lookupall(bindings, b) == Set([2])
        @test lookupall(bindings, a) == Set([3, 1])
    end
    ubind([a=>b]) do bindings
        # Binding two variables together is symetric:
        @test lookupall(bindings, a) == Set([b])
        @test lookupall(bindings, b) == Set([a])
    end
end

@testset "lookupequiv" begin
    a = V"a"
    b = V"b"
    c = V"c"
    ubind([a=>b, b=>a, c=>3, a=>c]) do bindings
        found, vars = lookupequiv(bindings, a)
        @test found == Set{Any}(3)
        @test vars == Set{AbstractVar}([a, b, c])
    end
    ubind([a=>b, a=>c]) do bindings
        found, vars = lookupequiv(bindings, c)
        @test found == Set{Any}()
        @test vars == Set{AbstractVar}([a, b, c])
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1GlusterfsPersistentVolumeSource <: SwaggerModel
    endpoints::Any # spec type: Union{ Nothing, String } # spec name: endpoints
    endpointsNamespace::Any # spec type: Union{ Nothing, String } # spec name: endpointsNamespace
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly

    function IoK8sApiCoreV1GlusterfsPersistentVolumeSource(;endpoints=nothing, endpointsNamespace=nothing, path=nothing, readOnly=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1GlusterfsPersistentVolumeSource, Symbol("endpoints"), endpoints)
        setfield!(o, Symbol("endpoints"), endpoints)
        validate_property(IoK8sApiCoreV1GlusterfsPersistentVolumeSource, Symbol("endpointsNamespace"), endpointsNamespace)
        setfield!(o, Symbol("endpointsNamespace"), endpointsNamespace)
        validate_property(IoK8sApiCoreV1GlusterfsPersistentVolumeSource, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(IoK8sApiCoreV1GlusterfsPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        o
    end
end # type IoK8sApiCoreV1GlusterfsPersistentVolumeSource

const _property_map_IoK8sApiCoreV1GlusterfsPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("endpoints")=>Symbol("endpoints"), Symbol("endpointsNamespace")=>Symbol("endpointsNamespace"), Symbol("path")=>Symbol("path"), Symbol("readOnly")=>Symbol("readOnly"))
const _property_types_IoK8sApiCoreV1GlusterfsPersistentVolumeSource = Dict{Symbol,String}(Symbol("endpoints")=>"String", Symbol("endpointsNamespace")=>"String", Symbol("path")=>"String", Symbol("readOnly")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1GlusterfsPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1GlusterfsPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1GlusterfsPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1GlusterfsPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1GlusterfsPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1GlusterfsPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1GlusterfsPersistentVolumeSource)
    (getproperty(o, Symbol("endpoints")) === nothing) && (return false)
    (getproperty(o, Symbol("path")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1GlusterfsPersistentVolumeSource }, name::Symbol, val)
end


text:
# This code tests two functions: keplerij! and kep_elliptic!
#using PyPlot
#include("../src/kepler_solver_derivative.jl")
#include("../src/ttv.jl")

@testset "keplerij" begin

# Next, try computing two-body Keplerian Jacobian:

n = 3
t0 = 7257.93115525
h  = 0.05
hbig  = big(h)
tmax = 600.0
#dlnq = 1e-8
dlnq = 1e-20
H = [3,1,1]

elements = readdlm("elements.txt",',')
#elements[2,1] = 0.75
elements[2,1] = 1.0
elements[3,1] = 1.0

m =zeros(n)
x0=zeros(NDIM,n)
v0=zeros(NDIM,n)

for k=1:n
  m[k] = elements[k,1]
end

for iter = 1:2

init = ElementsIC(elements,H,t0)
x0,v0,_ = init_nbody(init)
 if iter == 2
   # Reduce masses to trigger hyperbolic routine:
    m[1:n] *= 1e-1
    h = 0.05
    hbig = big(h)
 end
# Tilt the orbits a bit:
x0[2,1] = 5e-1*sqrt(x0[1,1]^2+x0[3,1]^2)
x0[2,2] = -5e-1*sqrt(x0[1,2]^2+x0[3,2]^2)
v0[2,1] = 5e-1*sqrt(v0[1,1]^2+v0[3,1]^2)
v0[2,2] = -5e-1*sqrt(v0[1,2]^2+v0[3,2]^2)

jac_ij = zeros(Float64,14,14)
dqdt_ij = zeros(Float64,14)
dqdt_num = zeros(BigFloat,14)
i=1 ; j=2
x = copy(x0) ; v=copy(v0)
# Predict values of s:
#println("Before first step: ",x," ",v)
xerror = zeros(NDIM,n); verror = zeros(NDIM,n)
d = Derivatives(Float64)
keplerij!(m,x,v,xerror,verror,i,j,h,jac_ij,dqdt_ij,d)
#println("After first step: ",x," ",v)
x0 = copy(x) ; v0 = copy(v)
xbig = big.(x) ; vbig=big.(v); mbig = big.(m)
xerr_big = zeros(BigFloat,NDIM,n) ; verr_big= zeros(BigFloat,NDIM,n)
xerror = zeros(NDIM,n); verror = zeros(NDIM,n)
keplerij!(m,x,v,xerror,verror,i,j,h,jac_ij,dqdt_ij,d)
# Now compute Jacobian with BigFloat precision:
jac_ij_big = zeros(BigFloat,14,14)
dqdt_ij_big = zeros(BigFloat,14)
KEPLER_TOL = sqrt(eps(big(1.0)))
dBig = Derivatives(BigFloat)
keplerij!(mbig,xbig,vbig,xerr_big,verr_big,i,j,hbig,jac_ij_big,dqdt_ij_big,dBig)
#println("jac_ij: ",convert(Array{Float64,2},jac_ij_big))
#println("jac_ij - jac_ij_big: ",convert(Array{Float64,2},jac_ij_big)-jac_ij)
println("max(jac_ij - jac_ij_big): ",maxabs(convert(Array{Float64,2},jac_ij_big)-jac_ij))

# xtest = copy(x0) ; vtest=copy(v0)
# keplerij!(m,xtest,vtest,i,j,h,jac_ij)
# println("Test of jacobian vs. none: ",maximum(abs(x-xtest)),maximum(abs(v-vtest)))

# Now, compute the derivatives numerically:
jac_ij_num = zeros(BigFloat,14,14)
xsave = big.(x)
vsave = big.(v)
msave = big.(m)

# Compute the time derivatives:
# Initial positions, velocities & masses:
xm = big.(x0)
vm = big.(v0)
mm = big.(msave)
dq = dlnq * hbig
hbig -= dq
keplerij!(mm,xm,vm,i,j,hbig)
xp = big.(x0)
vp = big.(v0)
hbig += 2dq
keplerij!(mm,xp,vp,i,j,hbig)
# Now x & v are final positions & velocities after time step
for k=1:3
  dqdt_num[   k] = .5*(xp[k,i]-xm[k,i])/dq
  dqdt_num[ 3+k] = .5*(vp[k,i]-vm[k,i])/dq
  dqdt_num[ 7+k] = .5*(xp[k,j]-xm[k,j])/dq
  dqdt_num[10+k] = .5*(vp[k,j]-vm[k,j])/dq
end
hbig = big(h)

for jj=1:3
  # Initial positions, velocities & masses:
  xm = big.(x0)
  vm = big.(v0)
  mm = big.(msave)
  dq = dlnq * xm[jj,i]
  if xm[jj,i] != 0.0
    xm[jj,i] -=  dq
  else
    dq = dlnq
    xm[jj,i] = -dq
  end
  keplerij!(mm,xm,vm,i,j,hbig)
  xp = big.(x0)
  vp = big.(v0)
  if xm[jj,i] != 0.0
    xp[jj,i] +=  dq
  else
    dq = dlnq
    xp[jj,i] = dq
  end
  keplerij!(mm,xp,vp,i,j,hbig)
  # now x & v are final positions & velocities after time step
  for k=1:3
    jac_ij_num[   k,  jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,  jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,  jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,  jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
  xm = big.(x0)
  vm = big.(v0)
  mm  = big.(msave)
  dq = dlnq * vm[jj,i]
  if vm[jj,i] != 0.0
    vm[jj,i] -=  dq
  else
    dq = dlnq
    vm[jj,i] = -dq
  end
  keplerij!(mm,xm,vm,i,j,hbig)
  xp = big.(x0)
  vp = big.(v0)
  mm  = big.(msave)
  if vp[jj,i] != 0.0
    vp[jj,i] +=  dq
  else
    dq = dlnq
    vp[jj,i] = dq
  end
  keplerij!(mm,xp,vp,i,j,hbig)
  for k=1:3
    jac_ij_num[   k,3+jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,3+jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,3+jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,3+jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
end

# Now vary mass of inner planet:
xm= big.(x0)
vm= big.(v0)
mm= big.(msave)
dq = mm[i]*dlnq
mm[i] -= dq
keplerij!(mm,xm,vm,i,j,hbig)
xp= big.(x0)
vp= big.(v0)
mp= big.(msave)
dq = mp[i]*dlnq
mp[i] += dq
keplerij!(mp,xp,vp,i,j,hbig)
for k=1:3
  jac_ij_num[   k,7] = .5*(xp[k,i]-xm[k,i])/dq
  jac_ij_num[ 3+k,7] = .5*(vp[k,i]-vm[k,i])/dq
  jac_ij_num[ 7+k,7] = .5*(xp[k,j]-xm[k,j])/dq
  jac_ij_num[10+k,7] = .5*(vp[k,j]-vm[k,j])/dq
end
# The mass doesn't change:
jac_ij_num[7,7] =  1.0

for jj=1:3
  # Now vary parameters of outer planet:
  xm = big.(x0)
  vm = big.(v0)
  mm = big.(msave)
  dq = dlnq * xm[jj,j]
  if xm[jj,j] != 0.0
    xm[jj,j] -=  dq
  else
    dq = dlnq
    xm[jj,j] = -dq
  end
  keplerij!(mm,xm,vm,i,j,hbig)
  xp = big.(x0)
  vp = big.(v0)
  if xp[jj,j] != 0.0
    xp[jj,j] +=  dq
  else
    dq = dlnq
    xp[jj,j] = dq
  end
  keplerij!(mm,xp,vp,i,j,hbig)
  for k=1:3
    jac_ij_num[   k,7+jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,7+jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,7+jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,7+jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
  xm= big.(x0)
  vm= big.(v0)
  mm = big.(msave)
  dq = dlnq * vm[jj,j]
  if vm[jj,j] != 0.0
    vm[jj,j] -=  dq
  else
    dq = dlnq
    vm[jj,j] = -dq
  end
  keplerij!(mm,xm,vm,i,j,hbig)
  xp= big.(x0)
  vp= big.(v0)
  if vp[jj,j] != 0.0
    vp[jj,j] +=  dq
  else
    dq = dlnq
    vp[jj,j] = dq
  end
  keplerij!(mm,xp,vp,i,j,hbig)
  for k=1:3
    jac_ij_num[   k,10+jj] = .5*(xp[k,i]-xm[k,i])/dq
    jac_ij_num[ 3+k,10+jj] = .5*(vp[k,i]-vm[k,i])/dq
    jac_ij_num[ 7+k,10+jj] = .5*(xp[k,j]-xm[k,j])/dq
    jac_ij_num[10+k,10+jj] = .5*(vp[k,j]-vm[k,j])/dq
  end
end

# Now vary mass of outer planet:
xm = big.(x0)
vm = big.(v0)
mm = big.(msave)
dq = mm[j]*dlnq
mm[j] -= dq
keplerij!(mm,xm,vm,i,j,hbig)
xp = big.(x0)
vp = big.(v0)
mp = big.(msave)
dq = mp[j]*dlnq
mp[j] += dq
keplerij!(mp,xp,vp,i,j,hbig)
for k=1:3
  jac_ij_num[   k,14] = .5*(xp[k,i]-xm[k,i])/dq
  jac_ij_num[ 3+k,14] = .5*(vp[k,i]-vm[k,i])/dq
  jac_ij_num[ 7+k,14] = .5*(xp[k,j]-xm[k,j])/dq
  jac_ij_num[10+k,14] = .5*(vp[k,j]-vm[k,j])/dq
end
# The mass doesn't change:
jac_ij_num[14,14] =  1.0

#println(jac_ij)
#println(jac_ij_num)
#println(jac_ij./jac_ij_num)
emax = 0.0; imax = 0; jmax = 0
for i=1:14, j=1:14
  if jac_ij[i,j] != 0.0
    diff = abs(convert(Float64,jac_ij_num[i,j])/jac_ij[i,j]-1.0)
    if  diff > emax
      emax = diff; imax = i; jmax = j
    end
  end
end
println("Maximum fractional error: ",emax," ",imax," ",jmax)
#println(jac_ij)
#println(convert(Array{Float64,2},jac_ij_num))
println("Maximum jac_ij error:   ",maxabs(convert(Array{Float64,2},jac_ij_num)-jac_ij))
println("Maximum jac_ij_big-jac_ij_num:   ",maxabs(convert(Array{Float64,2},jac_ij_num-jac_ij_big)))
println("Max dqdt error: ",maxabs(dqdt_ij-convert(Array{Float64,1},dqdt_num)))

@test isapprox(jac_ij_num,jac_ij;norm=maxabs)
@test isapprox(dqdt_ij,convert(Array{Float64,1},dqdt_num);norm=maxabs)
#println("dqdt: ",dqdt_ij," ",dqdt_num," Max error: ",maximum(dqdt_ij-dqdt_num))
end
end


text:
import JuMP
function simplextest(lib::Polyhedra.Library)
    hsim = HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0) ∩ HyperPlane([1, 1], 1)
    vsim = convexhull([0, 1], [1, 0])
    poly1 = polyhedron(hsim, lib)
    @test !isempty(poly1)
    center, radius = chebyshevcenter(poly1)
    @test center ≈ [1/2, 1/2]
    @test radius ≈ 1/2
    @test dim(poly1) == 1 # FIXME doing dim earlier makes chebyshevcenter fail
    inequality_fulltest(poly1, hsim)
    generator_fulltest(poly1, vsim)

    # Test incidence
    hpidx = first(eachindex(hyperplanes(poly1)))
    for ps in (incidentpoints(poly1, hpidx), get.(poly1, incidentpointindices(poly1, hpidx)))
        @test (ps == [[0, 1], [1, 0]] || ps == [[1, 0], [0, 1]])
    end
    for hsidx in eachindex(halfspaces(poly1))
        h = get(poly1, hsidx)
        if dot(h.a, [1, 0]) ≈ h.β
            expps = [[1, 0]]
        else
            @assert dot(h.a, [0, 1]) ≈ h.β
            expps = [[0, 1]]
        end
        for ps in (incidentpoints(poly1, hsidx), get.(poly1, incidentpointindices(poly1, hsidx)))
            @test ps == expps
        end
    end
    for pidx in eachindex(points(poly1))
        for hps in (incidenthyperplanes(poly1, pidx), get.(poly1, incidenthyperplaneindices(poly1, pidx)))
            @test hps == [HyperPlane([1, 1], 1)]
        end
        for hss in (incidenthalfspaces(poly1, pidx), get.(poly1, incidenthalfspaceindices(poly1, pidx)))
            h = hss[1]
            @test dot(h.a, get(poly1, pidx)) ≈ h.β
        end
    end

    @testset "Optimize with objective : max 2x_1" begin
        model, T = Polyhedra.layered_optimizer(Polyhedra.linear_objective_solver(poly1))
        x = MOI.add_variables(model, fulldim(poly1))
        MOI.add_constraint(model, MOI.VectorOfVariables(x),
                           Polyhedra.PolyhedraOptSet(poly1))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        @test MOI.MAX_SENSE == MOI.get(model, MOI.ObjectiveSense())
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm{T}(2, x[1])], zero(T)))
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.ObjectiveValue()) == 2
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.VariablePrimal(), x[1]) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x[2]) == 0
    end
    @testset "Optimize with objective : max x_1 + 3x_2" begin
        model, T = Polyhedra.layered_optimizer(Polyhedra.linear_objective_solver(poly1))
        x = MOI.add_variables(model, fulldim(poly1))
        MOI.add_constraint(model, MOI.VectorOfVariables(x),
                           Polyhedra.PolyhedraOptSet(poly1))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        @test MOI.MAX_SENSE == MOI.get(model, MOI.ObjectiveSense())
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm{T}(1, x[1]),
                     MOI.ScalarAffineTerm{T}(3, x[2])], zero(T)))
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.ObjectiveValue()) == 3
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.VariablePrimal(), x[1]) == 0
        @test MOI.get(model, MOI.VariablePrimal(), x[2]) == 1
    end

    poly2 = polyhedron(vsim, lib)
    @test dim(poly2) == 1
    @test !isempty(poly2)
    inequality_fulltest(poly2, hsim)
    generator_fulltest(poly2, vsim)

    @testset "x_1 cannot be 2" begin
        hempty = hsim ∩ HyperPlane([1, 0], 2)
        poly = polyhedron(hempty, lib)
        @test isempty(poly)
    end

    # We now add the vertex (0, 0)
    ext0 = convexhull([0, 0])
    @test collect(points(translate(ext0, [1, 0]))) == [[1, 0]]

    htri = HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0) ∩ HalfSpace([1, 1], 1)
    vtri = convexhull(vsim, ext0)
    convexhull!(poly1, ext0)
    inequality_fulltest(poly1, htri)
    generator_fulltest(poly1, vtri)
    convexhull!(poly2, ext0)
    inequality_fulltest(poly2, htri)
    generator_fulltest(poly2, vtri)

    # nonnegative orthant cut by x_1 + x_2 = 1
    vray = conichull([1, 0], [0, 1])
    poly3 = polyhedron(vray, lib)
    @test_throws ErrorException chebyshevcenter(poly3)
    @test dim(poly3) == 2

    @testset "Optimize with objective : min x_1 + x_2" begin
        model, T = Polyhedra.layered_optimizer(Polyhedra.linear_objective_solver(poly3))
        x = MOI.add_variables(model, fulldim(poly3))
        MOI.add_constraint(model, MOI.VectorOfVariables(x),
                           Polyhedra.PolyhedraOptSet(poly3))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        @test MOI.MIN_SENSE == MOI.get(model, MOI.ObjectiveSense())
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm{T}(1, x[1]),
                     MOI.ScalarAffineTerm{T}(1, x[2])], zero(T)))
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.ObjectiveValue()) == 0
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.VariablePrimal(), x[1]) == 0
        @test MOI.get(model, MOI.VariablePrimal(), x[2]) == 0
    end

    @testset "Optimize with objective : max x_2" begin
        model, T = Polyhedra.layered_optimizer(Polyhedra.linear_objective_solver(poly3))
        x = MOI.add_variables(model, fulldim(poly3))
        MOI.add_constraint(model, MOI.VectorOfVariables(x),
                           Polyhedra.PolyhedraOptSet(poly3))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        @test MOI.MAX_SENSE == MOI.get(model, MOI.ObjectiveSense())
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm{T}(1, x[2])], zero(T)))
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
        @test MOI.get(model, MOI.ObjectiveValue()) == 1
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        @test MOI.get(model, MOI.VariablePrimal(), x[1]) == 0
        @test MOI.get(model, MOI.VariablePrimal(), x[2]) == 1
    end

    hcutel = HyperPlane([1, 1], 1)
    hcut = intersect(hcutel)
    vcut = convexhull([1, 0]) + conichull(Line([1, -1]))
    @test !ininterior([1/2, 1/2], hcut)
    @test inrelativeinterior([1/2, 1/2], hcut)

    poly4 = copy(poly3)

    polycut3 = poly3 ∩ hcutel
    @test dim(polycut3) == 1
    inequality_fulltest(polycut3, hsim)
    generator_fulltest(polycut3, vsim)
    intersect!(poly3, hcutel)
    @test dim(poly3) == 1
    inequality_fulltest(poly3, hsim)
    generator_fulltest(poly3, vsim)

    # It should not have been cut as it is a copy of poly3
    @test dim(poly4) == 2

    polycut4 = poly4 ∩ hcut
    @test dim(polycut4) == 1
    inequality_fulltest(polycut4, hsim)
    generator_fulltest(polycut4, vsim)
    intersect!(poly4, hcut)
    @test dim(poly4) == 1
    inequality_fulltest(poly4, hsim)
    generator_fulltest(poly4, vsim)

    # FIXME needs float currently but should be updated
    # poly4 = project(poly1, [1; 0])
    # inequality_fulltest(poly4, [-1; 1], [0, 1], BitSet())
    # generator_fulltest(poly4, [0; 1], [])

    #\
    # \
    # |\
    # |_\
    #    \
    hlin = HalfSpace([1, 1], 1) ∩ HalfSpace([-1, -1], -1)
    plin = polyhedron(hlin, lib)
    @test dim(plin) == 1
    inequality_fulltest(plin, hcut)
    generator_fulltest(plin, vcut)
    #ineout = hrep(plin)
    #@test linset(ineout) == BitSet(1)
    vlin = convexhull([1, 0]) + conichull([1, -1], [-1, 1])
    plin = polyhedron(vlin, lib)
    inequality_fulltest(plin, hcut)
    generator_fulltest(plin, vcut)
    #extout = vrep(plin)
    #@test linset(extout) == BitSet(1)
end


text:
import Pkg

# parse command-line parameters
if length(ARGS) < 5
    println("Usage: octree.jl <input dir> <output dir> <levels> <channel> <voxelsize: x,y,z>")
    exit(1)
end

frompath = ARGS[1]
topath = ARGS[2]
nlevels = parse(Int,ARGS[3])
channel = parse(Int,ARGS[4])
voxelsize_um = map(x->parse(Float64,x), split(ARGS[5], ","))

import OctreeBuilder

OctreeBuilder.build_octree(frompath, topath, nlevels, channel)
OctreeBuilder.write_result(topath, voxelsize_um, nlevels)

max_stars_repo_name: JaneliaSciComp/jacs-tools-docker

max_issues_repo_name: JaneliaSciComp/jacs-tools-docker

max_forks_repo_name: JaneliaSciComp/jacs-tools-docker


text:
import GPUArrays

struct CuArrayBackend <: GPUArrays.GPUBackend end
GPUArrays.backend(::Type{<:CuArray}) = CuArrayBackend()

#Abstract GPU interface
struct CuKernelState end

@inline function GPUArrays.LocalMemory(::CuKernelState, ::Type{T}, ::Val{N}, ::Val{id}) where {T, N, id}
    ptr = CUDAnative._shmem(Val(id), T, Val(N))
    CuDeviceArray(N, DevicePtr{T, CUDAnative.AS.Shared}(ptr))
end

GPUArrays.AbstractDeviceArray(A::CUDAnative.CuDeviceArray, shape) = CUDAnative.CuDeviceArray(shape, pointer(A))

@inline GPUArrays.synchronize_threads(::CuKernelState) = CUDAnative.sync_threads()

GPUArrays.blas_module(::CuArray) = CuArrays.CUBLAS
GPUArrays.blasbuffer(x::CuArray) = x

"""
Blocks until all operations are finished on `A`
"""
GPUArrays.synchronize(A::CuArray) =
    CUDAdrv.synchronize()

for (i, sym) in enumerate((:x, :y, :z))
    for (f, fcu) in (
            (:blockidx, :blockIdx),
            (:blockdim, :blockDim),
            (:threadidx, :threadIdx),
            (:griddim, :gridDim)
        )
        fname = Symbol(string(f, '_', sym))
        cufun = Symbol(string(fcu, '_', sym))
        @eval GPUArrays.$fname(::CuKernelState) = CUDAnative.$cufun()
    end
end

# devices() = CUDAdrv.devices()
GPUArrays.device(A::CuArray) = CUDAdrv.device(CUDAdrv.CuCurrentContext())
GPUArrays.is_gpu(dev::CUDAdrv.CuDevice) = true
GPUArrays.name(dev::CUDAdrv.CuDevice) = string("CU ", CUDAdrv.name(dev))
GPUArrays.threads(dev::CUDAdrv.CuDevice) = CUDAdrv.attribute(dev, CUDAdrv.MAX_THREADS_PER_BLOCK)

GPUArrays.blocks(dev::CUDAdrv.CuDevice) =
    (CUDAdrv.attribute(dev, CUDAdrv.MAX_BLOCK_DIM_X),
     CUDAdrv.attribute(dev, CUDAdrv.MAX_BLOCK_DIM_Y),
     CUDAdrv.attribute(dev, CUDAdrv.MAX_BLOCK_DIM_Z))

GPUArrays.free_global_memory(dev::CUDAdrv.CuDevice) = CUDAdrv.Mem.info()[1]
GPUArrays.global_memory(dev::CUDAdrv.CuDevice) = CUDAdrv.totalmem(dev)
GPUArrays.local_memory(dev::CUDAdrv.CuDevice) = CUDAdrv.attribute(dev, CUDAdrv.TOTAL_CONSTANT_MEMORY)

function GPUArrays._gpu_call(::CuArrayBackend, f, A, args::Tuple,
                             blocks_threads::Tuple{T, T}) where {N, T <: NTuple{N, Integer}}
    blk, thr = blocks_threads
    @cuda blocks=blk threads=thr f(CuKernelState(), args...)
end

# Save reinterpret and reshape implementation use this in GPUArrays
GPUArrays.unsafe_reinterpret(::Type{T}, A::CuArray, size::NTuple{N, Integer}) where {T, N} =
    CuArray{T, N}(A.buf, size)


text:
function initialize!(p::Printer, problem, gradient, ξ, t)
    @info "Solving for" typeof(problem) gradient
end

function initialize!(writer::VTKWriter, problem, gradient, ξ, t)
    (; dir, filename) = writer
    ispath(dir) || mkdir(dir)
    writer.pvd = paraview_collection(joinpath(dir, filename))
    update!(writer, problem, gradient, ξ, t)
end

function initialize!(p::Plotter{T}, problem, gradient, ξ, t) where {T}
    femesh = problem.model.mesh
    M = problem.matrices.M
    ncompartment, nboundary = size(femesh.facets)
    ξ_cmpts = split_field(femesh, ξ)
    ξmax = maximum(abs, ξ)
    p.S₀ = sum(abs, M * ξ)
    TE = echotime(gradient)

    p.n = 1
    p.t[] = [t]
    p.ξ[] = copy(ξ)
    p.attenuation = Node(T[1])
    p.magnitude = [Node(T[]) for _ = 1:ncompartment, _ = 1:nboundary]
    p.phase = [Node(T[]) for _ = 1:ncompartment, _ = 1:nboundary]

    p.fig = Figure()

    if gradient isa ScalarGradient
        p.f[] = [gradient.profile(t)]
        ax = Axis(p.fig[1, 1]; xlabel = "t [μs]", title = "Time profile")
        xlims!(ax, 0, TE)
        ylims!(ax, -1.1, 1.1)
        lines!(ax, p.t, p.f)
    else
        grads = mapreduce(gradient, hcat, LinRange(0, TE, 200))
        pmin = min.(minimum(grads; dims = 2), zero(T))
        pmax = max.(maximum(grads; dims = 2), zero(T))
        inds = pmin .≈ pmax
        gmax = maximum(norm, grads)
        pmin[inds] .-= gmax / 2
        pmax[inds] .+= gmax / 2
        grad = Vec3f(gradient(t))
        p.g⃗ = Node([grad])
        p.g⃗_hist = Node([grad])
        ax = Axis3(p.fig[1, 1]; title = "Gradient [T/m]")
        ax.aspect = :data
        xlims!(ax, pmin[1], pmax[1])
        ylims!(ax, pmin[2], pmax[2])
        zlims!(ax, pmin[3], pmax[3])
        lines!(ax, p.g⃗_hist)
        arrows!(ax, [Point3f(0, 0, 0)], p.g⃗)
    end

    ax = Axis(p.fig[2, 1]; xlabel = "t [μs]", title = "Signal attenuation")
    xlims!(ax, 0, TE)
    # ylims!(ax, 1e-8, 1)
    ylims!(ax, 0, 1.1)
    lines!(ax, p.t, p.attenuation)

    gm = p.fig[1, 2] = GridLayout()
    ax = Axis3(gm[1, 1], title = "Magnetization (magnitude)")
    ax.aspect = :data
    m = nothing
    first = true
    colorrange = (0, ξmax)
    for icmpt = 1:ncompartment, iboundary = 1:nboundary
        facets = femesh.facets[icmpt, iboundary]
        points = femesh.points[icmpt]
        color = p.magnitude[icmpt, iboundary]
        color[] = abs.(ξ_cmpts[icmpt])
        if first
            m = mesh!(ax, points', facets'; color, shading = false, colorrange)
            first = false
        else
            mesh!(ax, points', facets'; color, shading = false, colorrange)
        end
    end
    Colorbar(gm[1, 2]; limits = colorrange)

    gm = p.fig[2, 2] = GridLayout()
    ax = Axis3(gm[1, 1], title = "Magnetization (phase-shift)")
    ax.aspect = :data
    m = nothing
    first = true
    # colormap = :cyclic_wrwbw_40_90_c42_n256
    colormap = :cyclic_mrybm_35_75_c68_n256
    colorrange = (-π, π)
    for icmpt = 1:ncompartment, iboundary = 1:nboundary
        facets = femesh.facets[icmpt, iboundary]
        points = femesh.points[icmpt]
        color = p.phase[icmpt, iboundary]
        color[] = angle.(ξ_cmpts[icmpt])
        if first
            m = mesh!(ax, points', facets'; color, shading = false, colorrange, colormap)
            first = false
        else
            mesh!(ax, points', facets'; color, shading = false, colorrange, colormap)
        end
    end
    Colorbar(gm[1, 2]; limits = colorrange, colormap)

    display(p.fig)

    p
end


text:
include("library.jl")
include("struct.jl")
include("Normalization.jl")
include("miscellaneous.jl")

max_stars_repo_name: Cortalak/Pissaro

max_issues_repo_name: Cortalak/Pissaro

max_forks_repo_name: Cortalak/Pissaro


text:
module Sub2Index

export @_

const SUB2NDIC = Dict(zip("₍₎₊₋₀₁₂₃₄₅₆₇₈₉₌ₐᵦᵪₑᵧₕᵢⱼₖₗₘₙₒₚᵩᵣᵨₛₔₜᵤᵥₓ", "()+-0123456789=aβχeγhijklmnopφrρsətuvx"))

function _sub2ind(src::AbstractString)
    idx = 1
    substrings = Any[]
    for m in eachmatch(r"[ᵢ-ᵪ₀-₎ₐ-ₜⱼ]+", src)
        offset = m.offset
        offset > idx && push!(substrings, src[idx:prevind(src, offset)])
        push!(substrings, '[')
        prechar = '\0'
        for c in m.match
            nextchar = SUB2NDIC[c]
            prechar == ')' && nextchar == '(' && push!(substrings, ',')
            push!(substrings, nextchar)
            prechar = nextchar
        end
        push!(substrings, ']')
        idx = offset + ncodeunits(m.match)
    end
    idx <= length(src) && push!(substrings, src[idx:end])
    join(substrings)
end

function sub2ind(src::Symbol)
    srcstr = string(src)
    convertedstr = _sub2ind(srcstr)
    esc(Meta.parse(convertedstr))
end

function sub2ind(src::Expr)
    Expr(src.head, (sub2ind(arg) for arg in src.args)...)
end

sub2ind(src) = src

"""
    @_(expr)

Replace all subscript charcters to by normal characters and encloses them in square brackets (i.e. indexing expression).
For example, `@_ Fₙ = Fₙ₋₁ + Fₙ₋₂` is converted to `F[n] = F[n-1] + F[n-2]`.

# Example

```julia-repl
julia> F = ones(Int, 10);
julia> @_ for n = 3:10 Fₙ = Fₙ₋₁ + Fₙ₋₂ end;
       # equivalent to `for n = 3:10 F[n] = F[n-1] + F[n-2] end`
julia> F
10-element Array{Int64,1}:
  1
  1
  2
  3
  5
  8
 13
 21
 34
 55

```

If each subscript is put in parentheses, it will be converted to the indexing of each dimension.
For example, `@_ Π₍ᵢ₎₍ⱼ₎` is converted to `Π[i, j]`.

# Example

```julia-repl

julia> Π = zeros(Int, (10, 10));
julia> @_ for i = 1:10 Π₍ᵢ₎₍₁₎ = Π₍ᵢ₎₍ᵢ₎ = 1 end;
julia> @_ for i = 3:10, j = 2:i-1 Π₍ᵢ₎₍ⱼ₎ = Π₍ᵢ₋₁₎₍ⱼ₋₁₎ + Π₍ᵢ₋₁₎₍ⱼ₎ end;
julia> Π
10×10 Array{Int64,2}:
 1  0   0   0    0    0   0   0  0  0
 1  1   0   0    0    0   0   0  0  0
 1  2   1   0    0    0   0   0  0  0
 1  3   3   1    0    0   0   0  0  0
 1  4   6   4    1    0   0   0  0  0
 1  5  10  10    5    1   0   0  0  0
 1  6  15  20   15    6   1   0  0  0
 1  7  21  35   35   21   7   1  0  0
 1  8  28  56   70   56  28   8  1  0
 1  9  36  84  126  126  84  36  9  1

```
"""
macro _(ex)
    sub2ind(ex)
end

end # module


text:
# Autogenerated wrapper script for embree_jll for x86_64-linux-musl
export libembree3

using oneTBB_jll
using GLFW_jll
using JpegTurbo_jll
using libpng_jll
JLLWrappers.@generate_wrapper_header("embree")
JLLWrappers.@declare_library_product(libembree3, "libembree3.so.3")
function __init__()
    JLLWrappers.@generate_init_header(oneTBB_jll, GLFW_jll, JpegTurbo_jll, libpng_jll)
    JLLWrappers.@init_library_product(
        libembree3,
        "lib/libembree3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    gravy(sequence)

Calculate GRAVY index of `sequence`.
Source: https://www.ncbi.nlm.nih.gov/pubmed/7108955
"""
function gravy(sequence)
    aacomposition = BioSequences.composition(sequence)
    gravy = 0.0
    for (aa, n) in aacomposition
        gravy += n * get(hydropathy, aa, 0.0)
    end
    return gravy / length(sequence)
end


text:
using DashBootstrapComponents

collapse = html_div([
    dbc_button(
        "Open collapse",
        id = "collapse-button",
        className = "mb-3",
        color = "primary",
        n_clicks = 0,
    ),
    dbc_collapse(
        dbc_card(dbc_cardbody("This content is hidden in the collapse")),
        id = "collapse",
        is_open = false,
    ),
]);

callback!(
    app,
    Output("collapse", "is_open"),
    Input("collapse-button", "n_clicks"),
    State("collapse", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using Printf
using Glob
using JLD2
using PyCall

plt = pyimport("matplotlib.pyplot")
gridspec = pyimport("matplotlib.gridspec")
ccrs = pyimport("cartopy.crs")
cmocean = pyimport("cmocean")

## Extract data

file = jldopen("tracer_advection_over_the_poles.jld2")

iterations = parse.(Int, keys(file["timeseries/t"]))

Nf = file["grid/faces"] |> keys |> length
Nx = file["grid/faces/1/Nx"]
Ny = file["grid/faces/1/Ny"]
Nz = file["grid/faces/1/Nz"]

## Plot!

for (n, i) in enumerate(iterations)
    @info "Plotting iteration $i/$(iterations[end]) (frame $n/$(length(iterations)))..."

    # fig = plt.figure(constrained_layout=true)
    fig = plt.figure()

    ax1 = plt.subplot2grid((3, 4), (2, 0))
    ax2 = plt.subplot2grid((3, 4), (2, 1))
    ax3 = plt.subplot2grid((3, 4), (1, 1))
    ax4 = plt.subplot2grid((3, 4), (1, 2))
    ax5 = plt.subplot2grid((3, 4), (0, 2))
    ax6 = plt.subplot2grid((3, 4), (0, 3))

    ax1.set_title("face 1")
    ax1.pcolormesh(file["timeseries/h/$i"][1][:, :, 1]', cmap=cmocean.cm.thermal, vmin=0, vmax=1000)
    ax1.get_xaxis().set_visible(false)
    ax1.get_yaxis().set_visible(false)

    ax2.set_title("face 2")
    ax2.pcolormesh(file["timeseries/h/$i"][2][:, :, 1]', cmap=cmocean.cm.thermal, vmin=0, vmax=1000)
    ax2.get_xaxis().set_visible(false)
    ax2.get_yaxis().set_visible(false)

    ax3.set_title("face 3")
    ax3.pcolormesh(file["timeseries/h/$i"][3][:, :, 1]', cmap=cmocean.cm.thermal, vmin=0, vmax=1000)
    ax3.get_xaxis().set_visible(false)
    ax3.get_yaxis().set_visible(false)

    ax4.set_title("face 4")
    ax4.pcolormesh(file["timeseries/h/$i"][4][:, :, 1]', cmap=cmocean.cm.thermal, vmin=0, vmax=1000)
    ax4.get_xaxis().set_visible(false)
    ax4.get_yaxis().set_visible(false)

    ax5.set_title("face 5")
    ax5.pcolormesh(file["timeseries/h/$i"][5][:, :, 1]', cmap=cmocean.cm.thermal, vmin=0, vmax=1000)
    ax5.get_xaxis().set_visible(false)
    ax5.get_yaxis().set_visible(false)

    ax6.set_title("face 6")
    ax6.pcolormesh(file["timeseries/h/$i"][6][:, :, 1]', cmap=cmocean.cm.thermal, vmin=0, vmax=1000)
    ax6.get_xaxis().set_visible(false)
    ax6.get_yaxis().set_visible(false)

    filename = @sprintf("tracer_avection_over_the_poles_%04d.png", n)
    plt.savefig(filename, dpi=200)
    plt.close(fig)
end

close(file)

run(`ffmpeg -y -i tracer_avection_over_the_poles_%04d.png -c:v libx264 -vf fps=10 -pix_fmt yuv420p out.mp4`)

# [rm(f) for f in glob("*.png")];


text:
module TextStylometry

include("corpus.jl")
include("features.jl")
include("measures.jl")
include("bootstrap.jl")

export Corpus, GenericDocument, SimpleDocument, ComplexityDocument,
    tokens, tokenfrequencies, equal_windows, bootstrap_window,
    bootstrap_ci, ttr, sttr, mattr, V, P, MWF, guiraud_r, herdan_c,
    dugast_k, a², dugast_U, LN, brunet_W

end # module


text:
# How to print SimplePolynomials? We take advantage of the
# excellent printing available in the Polynomials module!

import Base: string, show
export integerize

# The all_ones_denom function checks if all the coefficients
# are (Gaussian) integers, but perhaps of Rational type.

function all_ones_denom(list::Vector{T})::Bool where {T<:IntegerX}
    return true
end
#
# function all_ones_denom(list::Vector{Complex{T}})::Bool where T <: Integer
#     return true
# end

function all_ones_denom(list::Vector{T})::Bool where {T<:Real}
    return all(denominator.(list) .== 1)
end

function all_ones_denom(list::Vector{T})::Bool where {T<:Complex}
    return all_ones_denom(real.(list)) && all_ones_denom(imag.(list))
end

function all_ones_denom(p::SimplePolynomial)::Bool
    return all_ones_denom(p.data)
end

function _top(x::T) where {T<:Union{Rational,Integer}}
    return numerator(x)
end

function _top(x::T) where {T<:Complex}
    return numerator(real(x)) + numerator(imag(x)) * im
end

_top(x::AbstractMod) = x

"""
`integerize(p::SimplePolynomial)` is useful when `p`
has `Rational` coefficients all of whose denominators are `1`.
In that case, `integerize(p)` returns a polynomial equal to the
original, but with coefficients that are integers (or Gaussian integers).
However, if any denominators are not `1`, the original polynomial
is returned.
"""
function integerize(p::SimplePolynomial)
    if all_ones_denom(p)
        return SimplePolynomial(_top.(p.data))
    end
    return p
end

"""
`string(p::SimplePolynomial)` returns a nice `String` representation
of `p`.
"""
function string(p::SimplePolynomial)
    return string(Polynomial(integerize(p)))
end

function show(io::IO, p::SimplePolynomial)
    print(io, string(p))
end

export string3

function pstring(p::SimplePolynomial)
    nterms = count((x) -> (x != 0), p.data)
    if nterms <= 1
        return string(p)
    end
    return "(" * string(p) * ")"
end

function string(f::SimpleRationalFunction)
    return pstring(f.num) * " / " * pstring(f.den)
end

function show(io::IO, f::SimpleRationalFunction)
    print(io, string(f))
end

function string3(f::SimpleRationalFunction)
    top = string(numerator(f))
    bot = string(denominator(f))

    n1 = length(top)
    n2 = length(bot)
    n = max(n1, n2)

    pad1 = 0
    pad2 = 0

    if n1 < n
        pad1 = div(n - n1, 2)
    else
        pad2 = div(n - n2, 2)
    end

    line1 = " "^pad1 * top
    line2 = "-"^n
    line3 = " "^pad2 * bot
    NL = "\n"

    return line1 * NL * line2 * NL * line3
end


text:
"""
    ispangram(input)

Return `true` if `input` contains every alphabetic character (case insensitive).

"""
function ispangram(input)
    input_low = lowercase(input)
    lower_alpha = 'a':'z'

    return all(
        [(c in input_low) for c in lower_alpha])

end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions