text:
import Base: ==, *, /, +, -
import SparseArrays: sparse

"""
    LazySum([factors,] operators)

Lazy evaluation of sums of operators.

All operators have to be given in respect to the same bases. The field
`factors` accounts for an additional multiplicative factor for each operator
stored in the `operators` field.
"""
mutable struct LazySum{BL,BR,F,T} <: AbstractOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    factors::F
    operators::T
    function LazySum{BL,BR,F,T}(basis_l::BL,basis_r::BR,factors::F,operators::T) where {BL,BR,F,T}
        @assert length(operators)==length(factors)
        new(basis_l, basis_r, factors, operators)
    end
end
LazySum(basis_l::BL,basis_r::BR,factors::F,operators::T) where {BL,BR,F,T} =
    LazySum{BL,BR,F,T}(basis_l,basis_r,factors,operators)

function LazySum(factors, operators)
    for i = 2:length(operators)
        operators[1].basis_l == operators[i].basis_l || throw(IncompatibleBases())
        operators[1].basis_r == operators[i].basis_r || throw(IncompatibleBases())
    end
    Tf = promote_type(eltype(factors), eltype.(operators)...)
    factors_ = Tf.(factors)
    LazySum(operators[1].basis_l, operators[1].basis_r, factors_, operators)
end
LazySum(operators::AbstractOperator...) = LazySum(ones(promote_type(eltype.(operators)...), length(operators)), (operators...,))
#LazySum(factors::Vector{T}, operators::Vector{T2}) where {T<:Number,B1<:Basis,B2<:Basis,T2<:AbstractOperator{B1,B2}} = LazySum(complex(factors), (operators...,))
#LazySum(operators::AbstractOperator...) = LazySum(ones(ComplexF64, length(operators)), (operators...,))
LazySum(factors, operators::Vector) = LazySum(factors, (operators...,))
LazySum() = throw(ArgumentError("LazySum needs at least one operator!"))

Base.copy(x::LazySum) = LazySum(copy(x.factors), ([copy(op) for op in x.operators]...,))
Base.eltype(x::LazySum) = promote_type(eltype(x.factors), eltype.(x.operators)...)

dense(op::LazySum) = sum(op.factors .* dense.(op.operators))
SparseArrays.sparse(op::LazySum) = sum(op.factors .* sparse.(op.operators))

==(x::LazySum, y::LazySum) = (x.basis_l == y.basis_l && x.basis_r == y.basis_r && x.operators==y.operators && x.factors==y.factors)

# Arithmetic operations
+(a::LazySum{B1,B2}, b::LazySum{B1,B2}) where {B1,B2} = LazySum([a.factors; b.factors], (a.operators..., b.operators...))
+(a::LazySum{B1,B2}, b::LazySum{B3,B4}) where {B1,B2,B3,B4} = throw(IncompatibleBases())

-(a::LazySum) = LazySum(-a.factors, a.operators)
-(a::LazySum{B1,B2}, b::LazySum{B1,B2}) where {B1,B2} = LazySum([a.factors; -b.factors], (a.operators..., b.operators...))
-(a::LazySum{B1,B2}, b::LazySum{B3,B4}) where {B1,B2,B3,B4} = throw(IncompatibleBases())

*(a::LazySum, b::Number) = LazySum(b*a.factors, a.operators)
*(a::Number, b::LazySum) = LazySum(a*b.factors, b.operators)

/(a::LazySum, b::Number) = LazySum(a.factors/b, a.operators)

dagger(op::LazySum) = LazySum(conj.(op.factors), dagger.(op.operators))

tr(op::LazySum) = sum(op.factors .* tr.(op.operators))

function ptrace(op::LazySum, indices)
    check_ptrace_arguments(op, indices)
    rank = length(op.basis_l.shape) - length(indices)
    D = ([ptrace(op_i, indices) for op_i in op.operators]...,)
    LazySum(op.factors, D)
end

normalize!(op::LazySum) = (op.factors /= tr(op); op)

permutesystems(op::LazySum, perm) = LazySum(op.factors, ([permutesystems(op_i, perm) for op_i in op.operators]...,))

identityoperator(::Type{<:LazySum}, ::Type{S}, b1::Basis, b2::Basis) where S<:Number = LazySum(identityoperator(S, b1, b2))

# Fast in-place multiplication
function mul!(result::Ket{B1},a::LazySum{B1,B2},b::Ket{B2},alpha,beta) where {B1,B2}
    mul!(result,a.operators[1],b,alpha*a.factors[1],beta)
    for i=2:length(a.operators)
        mul!(result,a.operators[i],b,alpha*a.factors[i],1)
    end
    return result
end

function mul!(result::Bra{B2},a::Bra{B1},b::LazySum{B1,B2},alpha,beta) where {B1,B2}
    mul!(result,a,b.operators[1],alpha*b.factors[1],beta)
    for i=2:length(b.operators)
        mul!(result,a,b.operators[i],alpha*b.factors[i],1)
    end
    return result
end

function mul!(result::Operator{B1,B3},a::LazySum{B1,B2},b::Operator{B2,B3},alpha,beta) where {B1,B2,B3}
    mul!(result,a.operators[1],b,alpha*a.factors[1],beta)
    for i=2:length(a.operators)
        mul!(result,a.operators[i],b,alpha*a.factors[i],1)
    end
    return result
end
function mul!(result::Operator{B1,B3},a::Operator{B1,B2},b::LazySum{B2,B3},alpha,beta) where {B1,B2,B3}
    mul!(result,a,b.operators[1],alpha*b.factors[1],beta)
    for i=2:length(b.operators)
        mul!(result,a,b.operators[i],alpha*b.factors[i],1)
    end
    return result
end

function LinearAlgebra.ishermitian(op::LazySum)
    (all(isreal.(op.factors))&& all(ishermitian.(op.operators))) ||
    all(ishermitian.(op.factors .*op.operators))
end

function LinearAlgebra.issymmetric(op::LazySum)
    (all(isreal.(op.factors))&& all(issymmetric.(op.operators))) ||
    all(issymmetric.(op.factors .*op.operators))
end


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

"""
    calculate_normals(elements::Dict{Int, Vector{Int}},
                      element_types::Dict{Int, Symbol},
                      X::Dict{Int, Vector{Float64})

Given elements, element types and node locations, calculate nodal normals by
first calculating normal directions for each element and then averaging them
in nodes. As a result we get unique normal direction defined to each node.

# Notes
Only linear elements supported.

# Example

```jldoctest
X = Dict(1 => [7.0, 7.0], 2 => [4.0, 3.0], 3 => [0.0, 0.0])
elements = Dict(1 => [1, 2], 2 => [2, 3])
element_types = Dict(1 => :Seg2, 2 => :Seg2)
normals = calculate_normals(elements, element_types, X)

# output

Dict{Int64,Array{Float64,1}} with 3 entries:
  2 => [0.707107, -0.707107]
  3 => [0.6, -0.8]
  1 => [0.8, -0.6]

```

"""
function calculate_normals(elements, element_types, X)
    normals = empty(X)
    for (elid, elcon) in elements
        @assert element_types[elid] == :Seg2
        d = X[elcon[2]] - X[elcon[1]]
        n = [-d[2], d[1]]
        n /= norm(n)
        for nid in elcon
            haskey(normals, nid) || (normals[nid] = zeros(2))
            normals[nid] += n
        end
    end
    for (nid, n) in normals
        normals[nid] /= norm(normals[nid])
    end
    return normals
end

ext: jl


text:
using BuildExecutable
using Base.Test
script = joinpath(dirname(@__FILE__), "test.jl")
@test isfile(script)
targetdir = joinpath(dirname(@__FILE__),"test_dir")
mkdir(targetdir)
exename="test_executable"
try
    @test build_executable(exename, script, targetdir, "native"; force=false) == 0
    mv(targetdir, targetdir*"1")
    targetdir = targetdir*"1"
    out=readall(`$(joinpath(targetdir,exename * @windows? ".exe" : "")) -arg1 arg2 -arg3 arg4`)
    @test contains(out, "-arg1,arg2,-arg3,arg4")
    @test contains(out, "It works")
finally
    rm(targetdir, recursive=true)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "flex"
version = v"2.6.4"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/westes/flex/releases/download/v$(version)/flex-$(version).tar.gz",
                  "e87aae032bf07c26f85ac0ed3250998c37621d95f8bd748b31f15b33c45ee995")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/flex-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-static --enable-shared --disable-bootstrap
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;exclude=Sys.iswindows, experimental=true)

# The products that we will ensure are always built
products = [
    ExecutableProduct("flex", :flex),
    ExecutableProduct("flex++", :flexpp),
    LibraryProduct("libfl", :libfl)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# runtests.jl : run all MultilevelEstimators tests
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

using Distributed, MultilevelEstimators, Printf, Test

## unit tests ##

include("index.jl")

include("index_set.jl")

include("distribution.jl")

include("lattice.jl")

include("estimator.jl")


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function checkconfig(x, ximages, B2C, C2B)

nd = size(B2C,1);
n = size(x,2);

xx = B2C*x;
pp = B2C*ximages;

if nd ==2
    for i = 1:n
        xt = xx[:,i];
        if (0<=xt[1]) && (xt[1]<=1) && (0<=xt[2]) && (xt[2]<=1)
        else
            xp = xt + pp;
            for j = 1:9
                xq = xp[:,j];
                if (0<=xq[1]) && (xq[1]<=1) && (0<=xq[2]) && (xq[2]<=1)
                    xx[:,i] = xq;
                    break;
                end
            end
        end
    end
else
    for i = 1:n
        xt = xx[:,i];
        if (0<=xt[1]) && (xt[1]<=1) && (0<=xt[2]) && (xt[2]<=1) && (0<=xt[3]) && (xt[3]<=1)
        else
            xp = xt .+ pp;  # periodic images
            for j = 1:27
                xq = xp[:,j];
                if (0<=xq[1]) && (xq[1]<=1) && (0<=xq[2]) && (xq[2]<=1) && (0<=xq[3]) && (xq[3]<=1)
                    xx[:,i] = xq;
                    break;
                end
            end
        end
    end
end

x = C2B*xx;

return x

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
export iauPmat76
"""
Precession matrix from J2000.0 to a specified date, IAU 1976 model.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2 double       ending date, TT (Note 1)

Returned:
   rmatp       double[3][3] precession matrix, J2000.0 -> date1+date2

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix operates in the sense V(date) = RMATP * V(J2000),
   where the p-vector V(J2000) is with respect to the mean
   equatorial triad of epoch J2000.0 and the p-vector V(date)
   is with respect to the mean equatorial triad of the given
   date.

   3. Though the matrix method itself is rigorous, the precession
   angles are expressed through canonical polynomials which are
   valid only for a limited time span.  In addition, the IAU 1976
   precession rate is known to be imperfect.  The absolute accuracy
   of the present formulation is better than 0.1 arcsec from
   1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD,
   and remains below 3 arcsec for the whole of the period
   500BC to 3000AD.  The errors exceed 10 arcsec outside the
   range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to
   5600AD and exceed 1000 arcsec outside 6800BC to 8200AD.

Called:
   iauPrec76    accumulated precession angles, IAU 1976
   iauIr        initialize r-matrix to identity
   iauRz        rotate around Z-axis
   iauRy        rotate around Y-axis
   iauCr        copy r-matrix

References:

   Lieske, J.H., 1979, Astron.Astrophys. 73, 282.
   equations (6) & (7), p283.

   Kaplan,G.H., 1981. USNO circular no. 163, pA2.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPmat76(double date1, double date2, double rmatp[3][3])
function iauPmat76(date1::Real, date2::Real)
   rmatp = zeros(Float64, 3, 3)

   ccall((:iauPmat76, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Ptr{Cdouble}),
       convert(Float64, date1),
       convert(Float64, date2),
       rmatp)

   return SMatrix{3,3}(rmatp')
end


text:
# syntax: proto2
using ProtoBuf

export Tile_GeomType, Tile_Value, Tile_Feature, Tile_Layer, Tile

struct __enum_Tile_GeomType <: ProtoEnum
    UNKNOWN::Int32
    POINT::Int32
    LINESTRING::Int32
    POLYGON::Int32
    __enum_Tile_GeomType() = new(0,1,2,3)
end #struct __enum_Tile_GeomType
const Tile_GeomType = __enum_Tile_GeomType()

struct Tile_Value
    string_value::AbstractString
    float_value::Float32
    double_value::Float64
    int_value::Int64
    uint_value::UInt64
    sint_value::Int64
    bool_value::Bool
    Tile_Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #struct Tile_Value
const __wtype_Tile_Value = Dict(:sint_value => :sint64)
ProtoBuf.meta(t::Type{Tile_Value}) = ProtoBuf.meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, __wtype_Tile_Value, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
Base.hash(v::Tile_Value) = ProtoBuf.protohash(v)
Base.isequal(v1::Tile_Value, v2::Tile_Value) = ProtoBuf.protoisequal(v1, v2)
Base.==(v1::Tile_Value, v2::Tile_Value) = ProtoBuf.protoeq(v1, v2)

struct Tile_Feature
    id::UInt64
    tags::Array{UInt32,1}
    _type::Int32
    geometry::Array{UInt32,1}
    Tile_Feature(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #struct Tile_Feature
const __val_Tile_Feature = Dict(:id => 0, :_type => Tile_GeomType.UNKNOWN)
const __pack_Tile_Feature = Symbol[:tags,:geometry]
ProtoBuf.meta(t::Type{Tile_Feature}) = ProtoBuf.meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, __val_Tile_Feature, true, __pack_Tile_Feature, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
Base.hash(v::Tile_Feature) = ProtoBuf.protohash(v)
Base.isequal(v1::Tile_Feature, v2::Tile_Feature) = ProtoBuf.protoisequal(v1, v2)
Base.==(v1::Tile_Feature, v2::Tile_Feature) = ProtoBuf.protoeq(v1, v2)

struct Tile_Layer
    version::UInt32
    name::AbstractString
    features::Array{Tile_Feature,1}
    keys::Array{AbstractString,1}
    values::Array{Tile_Value,1}
    extent::UInt32
    Tile_Layer(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #struct Tile_Layer
const __req_Tile_Layer = Symbol[:version,:name]
const __val_Tile_Layer = Dict(:version => 1, :extent => 4096)
const __fnum_Tile_Layer = Int[15,1,2,3,4,5]
ProtoBuf.meta(t::Type{Tile_Layer}) = ProtoBuf.meta(t, __req_Tile_Layer, __fnum_Tile_Layer, __val_Tile_Layer, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
Base.hash(v::Tile_Layer) = ProtoBuf.protohash(v)
Base.isequal(v1::Tile_Layer, v2::Tile_Layer) = ProtoBuf.protoisequal(v1, v2)
Base.==(v1::Tile_Layer, v2::Tile_Layer) = ProtoBuf.protoeq(v1, v2)

struct Tile
    layers::Array{Tile_Layer,1}
    Tile(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #struct Tile
const __fnum_Tile = Int[3]
ProtoBuf.meta(t::Type{Tile}) = ProtoBuf.meta(t, ProtoBuf.DEF_REQ, __fnum_Tile, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
Base.hash(v::Tile) = ProtoBuf.protohash(v)
Base.isequal(v1::Tile, v2::Tile) = ProtoBuf.protoisequal(v1, v2)
Base.==(v1::Tile, v2::Tile) = ProtoBuf.protoeq(v1, v2)


text:
### ============== ### ============== ### ============== ###
### SHARED ARRAYS VERSION
### COUZIN COLLECTIVE MOTION MODEL
### ============== ### ============== ### ============== ###

### ============ OUTPUT FOLDER STRUCTURE ============ ###

function set_output_data_structure(path, N, o, a)

    parent_folder_path = joinpath(homedir(),"art_DATA",path)
    folder_path        = joinpath(parent_folder_path,"DATA","data_N_$(N)")
    reps_path          = joinpath(folder_path,"data_N_$(N)_o_$(o)_a_$(a)")

    try
        mkdir(joinpath(homedir(),"art_DATA"))
        println(joinpath(homedir(),"art_DATA"))
    catch error
        println("Main data folder already exists")
    end

    try
        mkdir(parent_folder_path)
        println(parent_folder_path)
    catch error
        println("Parent folder already exists")
    end

    try
        mkdir(joinpath(parent_folder_path,"DATA"))
        println(joinpath(parent_folder_path,"DATA"))
    catch error
        println("Parent folder already exists")
    end

    try
        mkdir(folder_path)
        println(folder_path)
    catch error
        println("Folder already exists")
    end

    try
        mkdir(reps_path)
        println(reps_path)
    catch error
        println("Parameter folder already exists")
    end

    return reps_path
end

### ============ COMPUTE RELATIVE (SQUARED) DISTANCES ============ ###

function calc_Rij(Rij::SharedArray, pos::SharedArray)

    @parallel for i in 1:3:length(pos)

        ri = div(i,3) + 1

        for j in (i+3):3:length(pos)

            rj = div(j,3) + 1

            Rij[rj,ri] = sqrt((pos[i]-pos[j])^2 + (pos[i+1]-pos[j+1])^2 + (pos[i+2]-pos[j+2])^2)
        end
    end

end

### ============ COMPUTE INTERACTIONS ============ ###

function compute_interactions(v_int::SharedArray, v_r::SharedArray, v_o::SharedArray, v_a::SharedArray, pos::SharedArray, vel::SharedArray, Rij::SharedArray, N::Int64, zor::Float64, zoo::Float64, zoa::Float64)

    F_Rij = Symmetric(Rij, :L)

    for id in first(localindexes(pos)):3:last(localindexes(pos))

        i = div(id,3)
        i_st = (i*N) # posicion inicial de la columna de distancias de la part i
        # println(i)

        k_int = zeros(Int64, 3) # number of neighbors

        v_int[3i+1] = 0.0
        v_int[3i+2] = 0.0
        v_int[3i+3] = 0.0

        v_r[3i+1] = 0.0
        v_r[3i+2] = 0.0
        v_r[3i+3] = 0.0

        v_o[3i+1] = 0.0
        v_o[3i+2] = 0.0
        v_o[3i+3] = 0.0

        v_a[3i+1] = 0.0
        v_a[3i+2] = 0.0
        v_a[3i+3] = 0.0

        for j in 1:N

            rij = F_Rij[i_st + j]

            if rij > 0.0 && rij <= zor

                v_r[3i+1] -= (pos[3(j-1)+1] - pos[3i+1]) / rij
                v_r[3i+2] -= (pos[3(j-1)+2] - pos[3i+2]) / rij
                v_r[3i+3] -= (pos[3(j-1)+3] - pos[3i+3]) / rij
                k_int[1] += 1

            else
                if rij > zor && rij <= zoo

                    v_o[3i+1] += vel[3(j-1)+1]
                    v_o[3i+2] += vel[3(j-1)+2]
                    v_o[3i+3] += vel[3(j-1)+3]
                    k_int[2] += 1

                elseif rij > zoo && rij <= zoa

                    v_a[3i+1] += (pos[3(j-1)+1] - pos[3i+1]) / rij
                    v_a[3i+2] += (pos[3(j-1)+2] - pos[3i+2]) / rij
                    v_a[3i+3] += (pos[3(j-1)+3] - pos[3i+3]) / rij
                    k_int[3] += 1

                end
            end

        end

        if k_int[1] > 0

            v_int[3i+1] = v_r[3i+1]
            v_int[3i+2] = v_r[3i+2]
            v_int[3i+3] = v_r[3i+3]

            v_norm = sqrt(v_int[3i+1]^2 + v_int[3i+2]^2 + v_int[3i+3]^2)

            v_int[3i+1] /= v_norm
            v_int[3i+2] /= v_norm
            v_int[3i+3] /= v_norm

        else

            if k_int[2] > 0 && k_int[3] > 0

                v_int[3i+1] = 0.5 * (v_o[3i+1] + v_a[3i+1])
                v_int[3i+2] = 0.5 * (v_o[3i+2] + v_a[3i+2])
                v_int[3i+3] = 0.5 * (v_o[3i+3] + v_a[3i+3])

                v_norm = sqrt(v_int[3i+1]^2 + v_int[3i+2]^2 + v_int[3i+3]^2)

                v_int[3i+1] /= v_norm
                v_int[3i+2] /= v_norm
                v_int[3i+3] /= v_norm

            elseif k_int[2] > 0 && k_int[3] == 0

                v_int[3i+1] = v_o[3i+1]
                v_int[3i+2] = v_o[3i+2]
                v_int[3i+3] = v_o[3i+3]

                v_norm = sqrt(v_int[3i+1]^2 + v_int[3i+2]^2 + v_int[3i+3]^2)

                v_int[3i+1] /= v_norm
                v_int[3i+2] /= v_norm
                v_int[3i+3] /= v_norm

            elseif k_int[3] > 0 && k_int[2] == 0

                v_int[3i+1] = v_a[3i+1]
                v_int[3i+2] = v_a[3i+2]
                v_int[3i+3] = v_a[3i+3]

                v_norm = sqrt(v_int[3i+1]^2 + v_int[3i+2]^2 + v_int[3i+3]^2)

                v_int[3i+1] /= v_norm
                v_int[3i+2] /= v_norm
                v_int[3i+3] /= v_norm

            end

        end

    end

end

### ============ UPDATE PARTICLE'S POSITIONS AND VELOCITIES ============ ###

function update_particles(v_r::SharedArray, pos::SharedArray, vel::SharedArray, η::Float64, θ::Float64, v0::Float64)

    q_r = Quaternion(zeros(Float64, 3))

    for id in first(localindexes(pos)):3:last(localindexes(pos))

        i = div(id, 3)
        # println(i)

        signal = [v_r[3i+1],v_r[3i+2],v_r[3i+3]]
        p_vel  = [vel[3i+1],vel[3i+2],vel[3i+3]]

        if norm(signal) > zero(Float64)

            # signal_angle = dot(p_vel, signal) / (norm(signal)*norm(p_vel))
            signal_angle = dot(p_vel, signal)

            signal_angle = ifelse( signal_angle < -1, -1, signal_angle)
            signal_angle = ifelse( signal_angle > 1, 1, signal_angle)

            q_r = qrotation(cross(p_vel, signal),  acos(signal_angle) + η * (2.0 * rand() * pi - pi)) * Quaternion(p_vel)

            # rot_ang = acos(signal_angle) + η * (2.0 * rand() * pi - pi)

            # if rot_ang < θ
            #     q_r = qrotation(cross(p_vel, signal),  rot_ang) * Quaternion(p_vel)
            # else
            #     # q_r = qrotation(cross(vel[i], v_r[i]), θ + η * (2.0 * rand() * pi - pi)) * Quaternion(vel[i])
            #     q_r = qrotation(cross(p_vel, signal), θ) * Quaternion(vel[i])
            # end

        else
            noise = randn(3)
            q_r = qrotation(cross(p_vel, noise),  η * (2.0 * rand() * pi - pi)) * Quaternion(p_vel)
        end

        u_vel = normalize([q_r.v1, q_r.v2, q_r.v3])

        vel[3i+1] = u_vel[1]
        vel[3i+2] = u_vel[2]
        vel[3i+3] = u_vel[3]

        pos[3i+1] += v0 * u_vel[1]
        pos[3i+2] += v0 * u_vel[2]
        pos[3i+3] += v0 * u_vel[3]

    end
end

### ============ SYSTEM UPDATE ============ ###

function evolve_system(pos::SharedArray, vel::SharedArray, v_int::SharedArray, v_r::SharedArray, v_o::SharedArray, v_a::SharedArray, R_ij::SharedArray, N::Int64, zor::Float64, zoo::Float64, zoa::Float64, η::Float64, θ::Float64, v0::Float64)

    calc_Rij(R_ij, pos)

    @sync begin
        for p in workers()
            @async remotecall_wait(compute_interactions, p, v_int, v_r, v_o, v_a, pos, vel, R_ij, N, zor, zoo, zoa)
        end
    end

    @sync begin
        for p in workers()
            @async remotecall_wait(update_particles, p, v_int, pos, vel, η, θ, v0)
        end
    end

end


text:
export Attr, Literal, Field
abstract Attr{T}
eltype{T}(::Attr{T}) = T

immutable Literal{T} <: Attr{T}
    value::T
end
name(l::Literal) = symbol(string(l.value))

convert{T<:Number}(::Type{Attr}, x::T) = Literal(x)

immutable Field{T} <: Attr{T}
    name::Symbol
end
Field(T, name) = Field{T}(name)
name(f::Field) = f.name

immutable As{T} <: Attr{T}
    name::Symbol
    field::Attr{T}
end
name(x::As) = x.name

immutable BinOp{op, L, R, T} <: Attr{T}
    l::Attr{L}
    r::Attr{R}
end
name{op}(f::BinOp{op}) = symbol(string(name(f.l), op, name(f.r)))

for op in [:(==), :(>=), :(<=), :(!=), :($), :(&), :(|)]
    quote
        function Base.$op{L, R}(l::Attr{L}, r::Attr{R})
            BinOp{$(Meta.quot(op)),L,R,Bool}(l,r)
        end
        function Base.$op{L, R}(l::Attr{L}, r::R)
            BinOp{$(Meta.quot(op)),L,R,Bool}(l,Literal(r))
        end
        function Base.$op{L, R}(l::L, r::Attr{R})
            BinOp{$(Meta.quot(op)),L,R,Bool}(Literal(l),r)
        end
    end |> eval
end

for op in [:+, :-, :*]
    @eval begin
        function Base.$op{L, R}(l::Attr{L}, r::Attr{R})
            BinOp{$(Meta.quot(op)),L,R,promote_type(L,R)}(l,r)
        end
        function Base.$op{L, R}(l::Attr{L}, r::R)
            BinOp{$(Meta.quot(op)),L,R,promote_type(L,R)}(l,Literal(r))
        end
        function Base.$op{L, R}(l::L, r::Attr{R})
            BinOp{$(Meta.quot(op)),L,R,promote_type(L,R)}(Literal(l),r)
        end
    end
end

abstract Aggregator <: Attr

immutable Count <: Aggregator end
name(::Count) = symbol("count()")

immutable Avg <: Aggregator
    attr::Attr
end
name(a::Avg) = symbol("avg($(name(a.attr)))")

immutable Sum <: Aggregator
    attr::Attr
end
name(s::Sum) = symbol("sum($(name(s.attr)))")

immutable ReduceAgg <: Aggregator
    f::Function
    attr::Attr
end


text:
using Documenter, DiffusionDefinition

makedocs(;
    modules=[DiffusionDefinition],
    format=Documenter.HTML(
        mathengine = Documenter.MathJax(
            Dict(
                :TeX => Dict(
                    :equationNumbers => Dict(
                        :autoNumber => "AMS"
                    ),
                    :Macros => Dict(
                        :dd => "{\\textrm d}",
                        :RR => "\\mathbb{R}",
                        :wt => ["\\widetilde{#1}", 1]
                    ),
                )
            )
        ),
        collapselevel = 1,
    ),
    pages=[
        "Home" => "index.md",
        "Get started" => joinpath("get_started", "package_overview.md"),
        "User manual" => Any[
            "Defining diffusions" => joinpath("manual", "defining_diffusions.md"),
            "Drift and volatility" => joinpath("manual", "drift_and_volatility.md"),
            "Sampling trajectories" => joinpath("manual", "sampling.md"),
            "Computing path-functionals" => joinpath("manual", "functionals_of_paths.md"),
            "Loading diffusions" => joinpath("manual", "load_diff.md"),
            "Utility functions" => joinpath("manual", "convenience_functions.md"),
            "Buffers" => joinpath("manual", "buffers.md"),
            "State space restrictions" => joinpath("manual", "state_space.md"),
            "Reparameterizations" => joinpath("manual", "reparameterization.md"),
            "Conjugate updates" =>joinpath("manual", "conjugate_updates.md"),
        ],
        "How to..." => Any[
            "... sample Brownian bridges?" => joinpath("how_to_guides", "sample_brownian_bridges.md"),
            "... sample Diffusion bridges?" => joinpath("how_to_guides", "sample_diffusion_bridges.md"),
            "... sample complex-valued diffusions?" => joinpath("how_to_guides", "sample_complex_valued.md"),
            "... customize plots of diffusion trajectories?" => joinpath("how_to_guides", "customize_my_diffusion_plots.md"),
            "... compute expectation of gradients of functionals?" => joinpath("how_to_guides", "expect_of_gradient.md"),
            "... use higher order SDE solvers for simulating trajectories?" => joinpath("how_to_guides", "use_higher_order_solvers.md"),
        ],
        "Tutorials" => Any[
            "Defining a diffusion process for a swinging pendulum" => joinpath("tutorials", "swinging_pendulum.md"),
            "In-place vs out-of-place methods" => joinpath("tutorials", "in_place_vs_out_of_place.md"),
            "Sampling of diffusion bridges" => joinpath("tutorials", "diffusion_bridges.md"),
        ],
        "Predefined Diffusions" => Any[
            "Lotka-Volterra system" => joinpath("predefined_processes", "lotka_volterra.md"),
            "Sine diffusion" => joinpath("predefined_processes", "sine.md"),
            "SIR model" => joinpath("predefined_processes", "sir.md"),
            "Lorenz63 system" => joinpath("predefined_processes", "lorenz63.md"),
            "Lorenz96 system" => joinpath("predefined_processes", "lorenz96.md"),
            "Favetto-Samson model" => joinpath("predefined_processes", "favetto_samson.md"),
            "Prokaryotic autoregulatory gene network" => joinpath("predefined_processes", "prokaryote.md"),
            "FitzHugh-Nagumo model" => joinpath("predefined_processes", "fitzhugh_nagumo.md"),
            "Jansen-Rit model" => joinpath("predefined_processes", "jansen_rit.md"),
        ],
        "Index" => "post_index.md",
    ],
    repo="https://github.com/JuliaDiffusionBayes/DiffusionDefinition.jl/blob/{commit}{path}#L{line}",
    sitename="DiffusionDefinition.jl",
    authors="Sebastiano Grazzi, Frank van der Meulen, Marcin Mider, Moritz Schauer",
    #assets=String[],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/DiffusionDefinition.jl",
)


text:
type SAGcbsd <: StepData
	d
	y
	numClasses::Int64
	c::Vector{Int64}
	batchSize::Int64
	getStep::Function
	stepString::AbstractString
	m
	function SAGcbsd(numVars::Int64,y, stepString::AbstractString, batchSize::Int64,c,m)
		d= sum(y)
		new(d, y, size(y)[1], c,batchSize, SAGcbComputation, stepString,m)
	end
end

function SAGcb(numVars::Int64, y, batchSize::Int64,c::Vector{Int64},m)
	SAGcbsd(numVars,y, "SAGcb.$batchSize",  batchSize,c,m)
end

function SAGcbComputation(x, k, gnum, sd::SAGcbsd, problem::Problem; outputLevel =0)

	outputLevel >0 && println("  Starting SAG CB step computation")
	indices = sample(1:sd.m,sd.batchSize;replace=false)
	outputLevel >1 && println("   indices = $indices")

	for i in indices

		outputLevel >1 && println("    Working on index $i")

		# step 1
		sd.d -= sd.y[sd.c[i]]

		(f,sampleG) = (problem.getSampleFunctionAt(i))(x)
		sd.y[sd.c[i]]=sampleG
		gnum += 1

		#step 3
		sd.d += sd.y[sd.c[i]]

	end

	outputLevel >0 && println("  Finished SAG CB step computation")

	(sd.d/sd.numClasses, gnum)
end


text:
# class constructors
/(::Type{T}, m::Integer) where T<:Integer = new_class(ZZmod{sintern(m),T}, T(m))
/(::Type{ZZ{T}}, m::Integer) where T<:Integer = T / T(m)
/(::Type{ZZ}, m::Integer) = mintype_for(m, 1, false) / m

# construction
isprimemod(Z::Type{<:ZZmod}) = isprime(modulus(Z))
basetype(::Type{<:ZZmod{m,T}}) where {m,T} = ZZ{wsigned(T)}
depth(::Type{<:ZZmod}) = depth(ZZ) + 1
_lcunit(a::ZZmod) = one(a)
issimpler(a::T, b::T) where T<:ZZmod = deg(a) < deg(b)

wsigned(::Type{T}) where T<:Signed = T
wsigned(::Type{T}) where T<:Unsigned = widen(signed(T))
wsigned(a::T) where T<:Integer = convert(wsigned(T), a)

function ZZmod{m,T}(a::Integer) where {m,T}
    mo = modulus(ZZmod{m,T})
    mo > 0 || throw(DomainError(m, "modulus > 0 required."))
    ZZmod{m,T}(T(mod(a, T(mo))), NOCHECK)
end
#ZZmod{m}(a::Integer) where {m} = ZZmod{m,typeof(m)}(oftype(m, a))
ZZmod(a::T, m::S) where {T,S} = (S/m)(S(a))
ZZmod{m,T}(a::ZZmod{m,T}) where {m,T} = a
ZZmod{m,T}(a::ZZmod{m,S}) where {m,T,S} = ZZmod{m,T}(a.val)
ZZmod{m,T}(a::ZZ) where {m,T} = ZZmod{m,T}(a.val)

copy(p::ZZmod) = typeof(p)(p.val)

#promotion and conversion
function _promote_rule(ZT::Type{ZZmod{m,S}}, ZS::Type{ZZmod{n,T}}) where {n,m,T,S}
    if modulus(ZT) == modulus(ZS)
        R = promote_type(T, S)
        if m == n
            ZZmod{promote(n,m)[1],R}
        elseif n isa Symbol
            ZZmod{n,R}
        else
            ZZmod{m,R}
        end
     else
         Union{}
    end
end
promote_rule(::Type{ZZmod{m,S}}, ::Type{T}) where {m,S,T<:Integer} = ZZmod{m,S}
promote_rule(::Type{ZZmod{m,S}}, ::Type{T}) where {m,S,T<:ZZ} = ZZmod{m,S}

function (::Type{ZT})(a::ZS) where {n,m,T,S,ZT<:ZZmod{n,T},ZS<:ZZmod{m,S}}
    mzt = modulus(ZT)
    mzs = modulus(ZS)
    if mzs % mzt == 0
        ZT(a.val % mzt)
    else
        throw(DomainError((ZT, a), "cannot convert $ZS to $ZT"))
    end
end

(::Type{T})(a::ZZmod) where T<:Integer = T(value(a))

"""
    value(a::ZZmod{m,T})

Return unique signed(T) integer in the range `[-m/2, m/2)` representing `a`, where `m` is the modulus.
"""
function value(a::ZZmod{X,T}) where {X,T}
    S = signed(T)
    v = a.val
    m = modulus(a)
    m1 = m >> 1
    m2 = m - m1
    v < m2 ? S(v) : S(v - m2) - S(m1)
end

# get type variable
modulus(t::Type{<:ZZmod{m,T}}) where {m,T} = m isa Integer ? T(m) : gettypevar(t).modulus
modulus(::T) where T<:ZZmod = modulus(T)

Base.isless(p::T, q::T) where T <: ZZmod = isless(p.val, q.val)

# arithmetic

function +(a::T, b::T) where T<:ZZmod
    mb = modulus(T) - b.val
    s = a.val < mb ? a.val + b.val : a.val - mb
    T(s, NOCHECK)
end
function -(a::T, b::T) where T<:ZZmod
    mb = modulus(T) - b.val
    s = a.val < b.val ? a.val + mb : a.val - b.val
    T(s, NOCHECK)
end
-(a::T) where T<:ZZmod = iszero(a.val) ? a : T(modulus(T) - a.val, NOCHECK)
*(a::T, b::T) where T<:ZZmod = T(mult_ZZmod(a.val, b.val, modulus(T)), NOCHECK)
*(a::T, b::Integer) where T<:ZZmod = T(mult_ZZmod(a.val, b, modulus(T)), NOCHECK)
*(a::Integer, b::T) where T<:ZZmod = T(mult_ZZmod(a, b.val, modulus(T)), NOCHECK)

inv(a::T) where T<:ZZmod = T(invmod2(a.val, modulus(T)), NOCHECK)
(a::T, n::Integer) where T<:ZZmod = T(powermod(a.val, n, modulus(T)), NOCHECK)

isunit(x::T) where T<:ZZmod = gcd(modulus(T), x.val) == 1
iszero(x::ZZmod) = iszero(x.val)
isone(x::ZZmod) = isone(x.val)
zero(::Type{<:ZZmod{m,S}}) where {m,S} = ZZmod{m,S}(zero(S), NOCHECK)
zero(::Type{ZZmod{m}}) where {m} = ZZmod{m,typeof(m)}(zero(m), NOCHECK)
one(::Type{<:ZZmod{m,S}}) where {m,S} = ZZmod{m,S}(one(S))
one(::Type{ZZmod{m}}) where {m} = ZZmod{m}(one(m))
==(a::ZZmod{m1},b::ZZmod{m2}) where {m1,m2} = modulus(a) == modulus(b) && a.val == b.val
hash(a::ZZmod, h::UInt) = hash(a.val, hash(modulus(a), h))
dimension(::Type{<:ZZmod}) = 1

# induced homomorphism
function (h::Hom{F,R,S})(p::Z) where {X,F,R,S,Z<:ZZmod{X,<:R}}
    m = F(modulus(Z))
    ZF = Z / m
    ZF(F(a.val))
end

# implementation of checked multiplication
function mult_ZZmod(a::T, b::T, m::T) where T<:Integer
    p, ov = Base.mul_with_overflow(a, b)
    ov || return mod(p, m)
    T(mod(widemul(a, b), m))
end
function mult_ZZmod(a::T1, b::T2, c::T) where {T1<:Integer,T2<:Integer,T<:Integer}
    a = T(mod(a,c))
    b = T(mod(b,c))
    mult_ZZmod(a, b, c)
end

# improved version of invmod
function invmod2(n::T, m::T) where T<:Integer
    iszero(m) && throw(DomainError(m, "`m` must not be 0."))
    p = _unsigned(abs(m))
    q = n >= 0 ? rem(_unsigned(n), p) : p - rem(_unsigned(-n), p)
    !iszero(q) && p != q || throw(DomainError((n, m), "Greatest common divisor is $m"))
    r = invmod_unsigned(q, p)
    m < 0 ? T(r) + m : T(r)
    # The postcondition is: mod(widemul(r, n), m) == mod(T(1), m) && div(r, m) == 0
end
# assumptions is, that  0 <= n < m
function invmod_unsigned(n::T, m::T) where T<:Union{Unsigned,BigInt}
    n = n < m ? n : rem(n, m)
    g, r, y = gcdx(n, m)
    g != 1 && throw(DomainError((n, m), "Greatest common divisor is $g."))
    # For unsigned T, x might be close to typemax; add m to force a wrap-around.
    if _underflow(r)
        r += m
    end
    r
    # The postcondition is: mod(widemul(r, n), m) == mod(T(1), m) && div(r, m) == 0
end

_underflow(r::Signed) = r < 0
_underflow(r::T) where T<:Unsigned = Base.hastypemax(T) && typemax(T)>>1 < r

_unsigned(::Type{BigInt}) = BigInt
_unsigned(T::Type) = unsigned(T)
_unsigned(x::BigInt) = x
_unsigned(x::Integer) = unsigned(x)

function Base.show(io::IO, a::ZZmod)
    v = a.val
    m = modulus(a)
    if m >= 100 && v > m÷2
        print(io, '-', signed(m-v), '°')
    else
        print(io, signed(v), '°')
    end
end

@inline function minlength_for(v::Integer, ::Val{1}, off::Bool)
    ndigits(v - off, base=2, pad=0)
end
@inline function minlength_for(v::Integer, ::Val{ex}, off::Bool) where ex
    l2n = log2(v) * ex
    Int(off ? ceil(l2n) : floor(l2n) + 1)
end

function mintype_for(v::Integer, ex::Integer, off::Bool)
    v isa BigInt && v > typemax(UInt128) && return BigInt
    len = minlength_for(v, Val(ex), off)
    for T in (UInt8, UInt16, UInt32, UInt64, UInt128)
        n = sizeof(T) * 8
        len <= n && return len < n ? signed(T) : T
    end
    BigInt
end


text:
using DataFrames
using CSV
using Statistics
using Base.Threads: @threads
using StatsBase: sample
using StaticArrays
using Revise
using BenchmarkTools
using EvoTrees
using EvoTrees: get_gain, get_edges, binarize, get_max_gain, update_grads!, grow_tree, grow_gbtree, SplitInfo, Tree, TrainNode, TreeNode, EvoTreeRegressor, predict, predict!, sigmoid
using EvoTrees: find_bags, update_bags!
using EvoTrees: find_split_static!, pred_leaf, update_hist!, find_split!, find_split_narrow!

# prepare a dataset
features = rand(100_000, 100)
# features = rand(1_000 10)
# x = cat(ones(20), ones(80)*2, dims=1)
# features =  hcat(x, features)

X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X,1))
𝑗 = collect(1:size(X,2))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]
𝑖 = collect(1:size(X_train,1))

# set parameters
params1 = EvoTreeRegressor(
    loss=:linear, metric=:mae,
    nrounds=1, nbins=32,
    λ = 0.0, γ=0.0, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=1.0, colsample=1.0)

# initial info
@time δ, δ² = zeros(SVector{params1.K, Float64}, size(X_train, 1)), zeros(SVector{params1.K, Float64}, size(X_train, 1))
𝑤 = zeros(SVector{1, Float64}, size(X_train, 1)) .+ 1
pred = zeros(SVector{params1.K,Float64}, size(X_train,1))
@time update_grads!(params1.loss, params1.α, pred, Y_train, δ, δ², 𝑤)
∑δ, ∑δ², ∑𝑤 = sum(δ[𝑖]), sum(δ²[𝑖]), sum(𝑤[𝑖])
@time gain = get_gain(params1.loss, ∑δ, ∑δ², ∑𝑤, params1.λ)
# @btime gain = get_gain($params1.loss, $∑δ, $∑δ², $∑𝑤, $params1.λ)

# initialize train_nodes
train_nodes = Vector{TrainNode{params1.K, Float64, Int}}(undef, 2^params1.max_depth-1)
for node in 1:2^params1.max_depth-1
    train_nodes[node] = TrainNode(0, SVector{params1.K, Float64}(fill(-Inf, params1.K)), SVector{params1.K, Float64}(fill(-Inf, params1.K)), SVector{1, Float64}(fill(-Inf, 1)), -Inf, [0], [0])
    # train_nodes[feat] = TrainNode(0, fill(-Inf, params1.K), fill(-Inf, params1.K), -Inf, -Inf, BitSet([0]), [0])
end

# initializde node splits info and tracks - colsample size (𝑗)
splits = Vector{SplitInfo{params1.K, Float64, Int}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    splits[feat] = SplitInfo{params1.K, Float64, Int}(gain, SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), -Inf, -Inf, 0, feat, 0.0)
end

# binarize data and create bags
@time edges = get_edges(X_train, params1.nbins)
@time X_bin = binarize(X_train, edges)
@time bags = Vector{Vector{BitSet}}(undef, size(𝑗, 1))
function prep(X_bin, bags)
    @threads for feat in 1:size(𝑗, 1)
         bags[feat] = find_bags(X_bin[:,feat])
    end
    return bags
end
@time prep(X_bin, bags);

# initialize histograms
feat=1
hist_δ = Vector{Vector{SVector{params1.K, Float64}}}(undef, size(𝑗, 1))
hist_δ² = Vector{Vector{SVector{params1.K, Float64}}}(undef, size(𝑗, 1))
hist_𝑤 = Vector{Vector{SVector{1, Float64}}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    hist_δ[feat] = zeros(SVector{params1.K, Float64}, length(bags[feat]))
    hist_δ²[feat] = zeros(SVector{params1.K, Float64}, length(bags[feat]))
    hist_𝑤[feat] = zeros(SVector{1, Float64}, length(bags[feat]))
end

# grow single tree
#  0.135954 seconds (717.54 k allocations: 15.219 MiB)
@time train_nodes[1] = TrainNode(1, SVector(∑δ), SVector(∑δ²), SVector(∑𝑤), gain, 𝑖, 𝑗)
@time tree = grow_tree(bags, δ, δ², 𝑤, hist_δ, hist_δ², hist_𝑤, params1, train_nodes, splits, edges, X_bin)
# @btime tree = grow_tree($bags, $δ, $δ², $𝑤, $hist_δ, $hist_δ², $hist_𝑤, $params1, $train_nodes, $splits, $tracks, $edges, $X_bin)
@time pred_train = predict(tree, X_train, params1.K)
# 705.901 μs (18 allocations: 626.08 KiB)
@btime pred_train = predict($tree, $X_train, $params1.K)
@time pred_leaf_ = pred_leaf(params1.loss, train_nodes[1], params1, δ²)
# 1.899 ns (0 allocations: 0 bytes)
@btime pred_leaf_ = pred_leaf($params1.loss, $train_nodes[1], $params1, $δ²)
# @btime pred_train = predict($tree, $X_train, params1.K)

@time model = grow_gbtree(X_train, Y_train, params1, print_every_n = 1)
# @btime model = grow_gbtree($X_train, $Y_train, $params1, print_every_n = 1)
@time pred_train = predict(model, X_train)
@btime pred_train = predict($model, $X_train)

params1 = Params(:linear, 10, λ, γ, 0.1, 5, min_weight, rowsample, colsample, nbins)
@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 1, metric=:mae)
@btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval, print_every_n = 1, metric=:mae)

@time pred_train = predict(model, X_train)
sqrt(mean((pred_train .- Y_train) .^ 2))

#############################################
# Quantiles with turbo
#############################################

𝑖_set = BitSet(𝑖);
@time bags = prep(X_bin, bags);
# target: find_split_turbo! in 0.001 sec for 100_000 observations
feat = 1
typeof(bags[feat][1])
# initialise node, info and tracks
train_nodes[1] = TrainNode(1, ∑δ, ∑δ², ∑𝑤, gain, 𝑖, 𝑗)
splits[feat] = SplitInfo{params1.K, Float64, Int}(gain, SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), ∑δ, ∑δ², ∑𝑤, -Inf, -Inf, 0, feat, 0.0)

# 492.199 μs (343 allocations: 6.83 KiB)
# 1.038 ms (343 allocations: 6.83 KiB) for 200_000
splits[feat]
@time find_split_static!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], train_nodes[1].𝑖)
# @code_warntype find_split_static!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], train_nodes[1].𝑖)
@btime find_split_static!($hist_δ[feat], $hist_δ²[feat], $hist_𝑤[feat], $bags[feat], $view(X_bin,:,feat), $δ, $δ², $𝑤, $∑δ, $∑δ², $∑𝑤, $params1, $splits[feat], $edges[feat], $train_nodes[1].𝑖)

feat = 2
typeof(bags[feat][1])
# initialise node, info and tracks
train_nodes[1] = TrainNode(1, ∑δ, ∑δ², ∑𝑤, gain, BitSet(𝑖), 𝑗)
splits[feat] = SplitInfo{Float64, Int}(gain, SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), ∑δ, ∑δ², ∑𝑤, -Inf, -Inf, 0, feat, 0.0)
@time find_split_static!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], train_nodes[1].𝑖)
@btime find_split_static!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], train_nodes[1].𝑖)

function find_split_static()
    @threads for feat in 𝑗
        find_split_static!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], train_nodes[1].𝑖)
    end
end
@btime find_split_tot()

# find split wide
set = Int64.(train_nodes[1].𝑖)
# @time update_hist!(hist_δ, hist_δ², hist_𝑤, X_bin, δ, δ², 𝑤, set)
@time update_hist!(hist_δ, hist_δ², hist_𝑤, X_bin, δ, δ², 𝑤, set, feat)
# 9.933 ms (1 allocation: 96 bytes)
@btime update_hist!($hist_δ, $hist_δ², $hist_𝑤, $X_bin, $δ, $δ², $𝑤, $set, $feat)
function update_hist_wide()
    @threads for j in 𝑗
        update_hist!(hist_δ, hist_δ², hist_𝑤, X_bin, δ, δ², 𝑤, train_nodes[1].𝑖, j)
    end
end
# 100K: 6.502 ms (0 allocations: 0 bytes)
@btime update_hist_wide()

@time find_split!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], feat)
@btime find_split!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], feat)
function find_split_wide()
    @threads for feat in 𝑗
        find_split!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], feat)
    end
end
@btime find_split_wide()

function find_split_tot()
    @threads for feat in 𝑗
        update_hist!(hist_δ, hist_δ², hist_𝑤, X_bin, δ, δ², 𝑤, train_nodes[1].𝑖, feat)
        find_split!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], feat)
    end
end
@btime find_split_tot()

# find split narrow
@time find_split_narrow!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], set)
@btime find_split_narrow!($hist_δ[feat], $hist_δ²[feat], $hist_𝑤[feat], $bags[feat], $view(X_bin,:,feat), $δ, $δ², $𝑤, $∑δ, $∑δ², $∑𝑤, $params1, $splits[feat], $edges[feat], $set)
function test_narrow()
    # set = Int64.(train_nodes[1].𝑖)
    @threads for feat in 𝑗
        find_split_narrow!(hist_δ[feat], hist_δ²[feat], hist_𝑤[feat], bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], edges[feat], set)
    end
end
# 100K: 18.765 ms (300 allocations: 9.38 KiB)
@btime test_narrow()

using Base.Threads: @threads
function test()
    # x = collect(1:8)
    x = BitSet(1:8)
    # x = Int64.(x)
    res = zero(Int)
    @threads for i in x
        res += i
    end
    return res
end

test()

length(union(train_nodes[1].bags[1][1:13]...))
length(union(train_nodes[1].bags[1][1:13]...))
length(new_bags[2])
length(new_bags[2][1])
length(bags[2][32])
typeof(bags)
@btime update_bags_intersect($new_bags, $bags, $union(train_nodes[1].bags[1][1:13]...))
length(new_bags[2])
length(new_bags[2][2])
length(bags[2][1])

# extract the best feat from bags, and join all the underlying bins up to split point
best_bag = bags[1]
bins_L = union(best_bag[1:4]...)

function set_1(x, y)
    intersect!(x, y)
    return x
end

x = rand(UInt32, 100_000)
y = rand(x, 1000)

x_set = BitSet(x);
y_set = BitSet(y);

@btime set_1(x, y)
@btime set_1(x_set, y)

x = rand([1,2,3,4,5,6,7,8,9,10, 11,12], 1000)
x = rand(1000)
x_edges = quantile(x, (0:10)/10)
x_edges = unique(x_edges)
x_edges = x_edges[2:(end-1)]

length(x_edges)

x_bin = searchsortedlast.(Ref(x_edges), x) .+ 1
using StatsBase
x_map = countmap(x_bin)

x = reshape(x, (1000, 1))
x_edges = get_edges(x)
unique(quantile(view(X, :,i), (0:nbins)/nbins))[2:(end-1)]
x_bin = searchsortedlast.(Ref(x_edges[1]), x[:,1]) .+ 1
x_map = countmap(x_bin)

edges = get_edges(X, 32)
X_bin = zeros(UInt8, size(X))
@btime binindices(X[:,1], edges[1])
@btime X_bin = binarize(X, edges)

using StatsBase
x_map = countmap(x_bin)

x_edges[1]


text:
# Autogenerated wrapper script for GAP_pkg_digraphs_jll for x86_64-linux-gnu
export digraphs

using GAP_jll
using GAP_lib_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_digraphs")
JLLWrappers.@declare_file_product(digraphs)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll)
    JLLWrappers.@init_file_product(
        digraphs,
        "lib/gap/digraphs.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module MathOptInterfaceMosek

import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
using Mosek
#using Mosek.Ext

using Compat # for findall

include("LinkedInts.jl")

const DEBUG = false

mosek_block_type_unallocated = 0
mosek_block_type_zero   = 1
mosek_block_type_nonneg = 2
mosek_block_type_nonpos = 3
mosek_block_type_range  = 4
mosek_block_type_qcone  = 5
mosek_block_type_rqcone = 6
mosek_block_type_psd    = 7
mosek_block_type_integer = 8

problemtype_linear    = 0
problemtype_conic     = 1
problemtype_quadratic = 2

import MathOptInterface

boundflag_lower = 0x1
boundflag_upper = 0x2
boundflag_cone  = 0x4
boundflag_int   = 0x8
boundflag_all   = 0x0f

# Mapping of all constraint types to its index
struct ConstraintMap
    x_lessthan           :: Dict{Int64,Int} # (SingleVariable,LessThan) -> constraint number
    x_greaterthan        :: Dict{Int64,Int} # (SingleVariable,GreaterThan) -> constraint number
    x_equalto            :: Dict{Int64,Int} # (SingleVariable,EqualTo) -> constraint number
    x_interval           :: Dict{Int64,Int} # (SingleVariable,Interval) -> constraint number
    x_nonpositives       :: Dict{Int64,Int} # (SingleVariable,Nonpositives) -> constraint number
    x_nonnegatives       :: Dict{Int64,Int} # (SingleVariable,Nonnegatives) -> constraint number
    x_binary             :: Dict{Int64,Int} # (SingleVariable,ZeroOne) -> constraint number
    x_integer            :: Dict{Int64,Int} # (SingleVariable,Integer) -> constraint number

    xs_nonpositives      :: Dict{Int64,Int} # (VectorOfVariables,Nonpositives) -> constraint number
    xs_nonnegatives      :: Dict{Int64,Int} # (VectorOfVariables,Nonnegatives) -> constraint number
    xs_zeros             :: Dict{Int64,Int} # (VectorOfVariables,Zeros) -> constraint number
    xs_reals             :: Dict{Int64,Int} # (VectorOfVariables,Reals) -> constraint number
    xs_qcone             :: Dict{Int64,Int} # (VectorOfVariables,SecondOrderCone) -> constraint number
    xs_rqcone            :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    xs_psdconetriangle   :: Dict{Int64,Int} # (VectorOfVariables,PositiveSemidefiniteConeTriangle) -> constraint number
    xs_ppowcone          :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    xs_dpowcone          :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    xs_pexpcone          :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    xs_dexpcone          :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number

    axb_lessthan         :: Dict{Int64,Int} # (ScalarAffineFunction,LessThan) -> constraint number
    axb_greaterthan      :: Dict{Int64,Int} # (ScalarAffineFunction,GreaterThan) -> constraint number
    axb_equalto          :: Dict{Int64,Int} # (ScalarAffineFunction,EqualTo) -> constraint number
    axb_interval         :: Dict{Int64,Int} # (ScalarAffineFunction,Interval) -> constraint number
    axb_nonpositives     :: Dict{Int64,Int} # (ScalarAffineFunction,Nonpositives) -> constraint number
    axb_nonnegatives     :: Dict{Int64,Int} # (ScalarAffineFunction,Nonnegatives) -> constraint number
    axb_binary           :: Dict{Int64,Int} # (ScalarAffineFunction,ZeroOne) -> constraint number
    axb_integer          :: Dict{Int64,Int} # (ScalarAffineFunction,Integer) -> constraint number

    axbs_nonpositives    :: Dict{Int64,Int} # (VectorAffineFunction,Nonpositives) -> constraint number
    axbs_nonnegatives    :: Dict{Int64,Int} # (VectorAffineFunction,Nonnegatives) -> constraint number
    axbs_zeros           :: Dict{Int64,Int} # (VectorAffineFunction,Zeros) -> constraint number
    axbs_reals           :: Dict{Int64,Int} # (VectorAffineFunction,Reals) -> constraint number
    axbs_qcone           :: Dict{Int64,Int} # (VectorAffineFunction,SecondOrderCone) -> constraint number
    axbs_rqcone          :: Dict{Int64,Int} # (VectorAffineFunction,RotatedSecondOrderCone) -> constraint number
    axbs_psdconetriangle :: Dict{Int64,Int} # (VectorAffineFunction,PositiveSemidefiniteConeTriangle) -> constraint number
    axbs_ppowcone        :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    axbs_dpowcone        :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    axbs_pexpcone        :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
    axbs_dexpcone        :: Dict{Int64,Int} # (VectorOfVariables,RotatedSecondOrderCone) -> constraint number
end

ConstraintMap() = ConstraintMap([Dict{Int64,Int}() for i in 1:38]...)
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.LessThan{Float64}}) =                cm.x_lessthan
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.GreaterThan{Float64}}) =             cm.x_greaterthan
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.EqualTo{Float64}}) =                 cm.x_equalto
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.Interval{Float64}}) =                cm.x_interval
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.Nonpositives}) =                     cm.x_nonpositives
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.Nonnegatives}) =                     cm.x_nonnegatives
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.ZeroOne}) =                          cm.x_binary
select(cm::ConstraintMap,::Type{MOI.SingleVariable},               ::Type{MOI.Integer}) =                          cm.x_integer
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.Nonpositives}) =                     cm.xs_nonpositives
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.Nonnegatives}) =                     cm.xs_nonnegatives
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.Zeros}) =                            cm.xs_zeros
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.Reals}) =                            cm.xs_reals
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.SecondOrderCone}) =                  cm.xs_qcone
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.RotatedSecondOrderCone}) =           cm.xs_rqcone
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.PositiveSemidefiniteConeTriangle}) = cm.xs_psdconetriangle
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.PowerCone{Float64}}) =               cm.xs_ppowcone
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.DualPowerCone{Float64}}) =           cm.xs_dpowcone
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.ExponentialCone}) =                  cm.xs_pexpcone
select(cm::ConstraintMap,::Type{MOI.VectorOfVariables},            ::Type{MOI.DualExponentialCone}) =              cm.xs_dexpcone
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.LessThan{Float64}}) =                cm.axb_lessthan
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.GreaterThan{Float64}}) =             cm.axb_greaterthan
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.EqualTo{Float64}}) =                 cm.axb_equalto
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.Interval{Float64}}) =                cm.axb_interval
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.ZeroOne}) =                          cm.axb_binary
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.Integer}) =                          cm.axb_integer
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.Nonpositives}) =                     cm.axb_nonpositives
select(cm::ConstraintMap,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.Nonnegatives}) =                     cm.axb_nonnegatives
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.Nonpositives}) =                     cm.axbs_nonpositives
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.Nonnegatives}) =                     cm.axbs_nonnegatives
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.Zeros}) =                            cm.axbs_zeros
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.Reals}) =                            cm.axbs_reals
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.SecondOrderCone}) =                  cm.axbs_qcone
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.RotatedSecondOrderCone}) =           cm.axbs_rqcone
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.PositiveSemidefiniteConeTriangle}) = cm.axbs_psdconetriangle
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.PowerCone{Float64}}) =               cm.axbs_ppowcone
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.DualPowerCone{Float64}}) =           cm.axbs_dpowcone
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.ExponentialCone}) =                  cm.axbs_pexpcone
select(cm::ConstraintMap,::Type{MOI.VectorAffineFunction{Float64}},::Type{MOI.DualExponentialCone}) =              cm.axbs_dexpcone

Base.getindex(cm::ConstraintMap,r :: MOI.ConstraintIndex{F,D}) where {F,D} = select(cm,F,D)[r.value]

struct MosekSolution
    whichsol :: Soltype
    solsta   :: Solsta
    prosta   :: Prosta

    xxstatus :: Vector{Stakey}
    xx       :: Vector{Float64}
    slx      :: Vector{Float64}
    sux      :: Vector{Float64}
    snx      :: Vector{Float64}

    cstatus  :: Vector{Stakey}
    xc       :: Vector{Float64}
    slc      :: Vector{Float64}
    suc      :: Vector{Float64}
    y        :: Vector{Float64}
end

"""
    MosekModel <: MathOptInterface.AbstractModel

Linear variables and constraint can be deleted. For some reason MOSEK
does not support deleting PSD variables.

Note also that adding variables and constraints will permanently add
some (currently between 1 and 3) Int64s that a `delete!` will not
remove. This ensures that Indices (Variable and constraint) that
are deleted are thereafter invalid.
"""
mutable struct MosekModel  <: MOI.AbstractOptimizer
    task :: Mosek.MSKtask
    ## Options passed in `MosekOptimizer` that are used to create a new task
    ## in `MOI.empty!`:
    # Should Mosek output be ignored or printed ?
    be_quiet :: Bool
    # Integer parameters, i.e. parameters starting with `MSK_IPAR_`
    ipars :: Dict{String, Integer}
    # Floating point parameters, i.e. parameters starting with `MSK_DPAR_`
    dpars :: Dict{String, AbstractFloat}
    # String parameters, i.e. parameters starting with `MSK_SPAR_`
    spars :: Dict{String, AbstractString}

    """
    Number of variables explicitly created by user
    """
    publicnumvar :: Int

    """
    """
    constrmap :: ConstraintMap

    """
    """
    constrnames :: Dict{String,Vector{MOI.ConstraintIndex}}

    """
        The total length of `x_block` matches the number of variables in
    the underlying task, and the number of blocks corresponds to the
    number variables allocated in the Model.
    """
    x_block      :: LinkedInts

    """
    One entry per scalar variable in the task indicating which bound
    types are defined
    """
    x_boundflags :: Vector{Int}

    """
    One entry per scalar variable in the task defining the number of
    variable constraints imposed on that variable. It is not allowed
    to delete a variable without deleting the constraints on it first
    (to get around the problem of deleting roots in conic constraints).
    """
    x_numxc :: Vector{Int}

    """
    One entry per variable-constraint
    """
    xc_block     :: LinkedInts

    """
    One entry per variable-constraint block indicating which bound
    types it defines. The values are binary ORed `boundflag_...` values.
    """
    xc_bounds    :: Vector{UInt8}
    """
    One entry per variable-constraint block indicating which cone it belongs to, 0 if none.
    """
    xc_coneid    :: Vector{Int}

    """
    One entry per scalar variable-constraint, defining which native
    variables the bound block covers.
    """
    xc_idxs      :: Vector{Int}

    ###########################
    """
    One scalar entry per constraint in the underlying task. One block
    per constraint allocated in the Model.
    """
    c_block :: LinkedInts

    """
    One entry per allocated scalar constraint. Defines the fixed term
    on the left-hand for each scalar constraint.
    """
    c_constant :: Vector{Float64}

    """
    One entry per allocated scalar constraint.
    Each element is either
    - 0: meaning: no slack, when domain is defined directly as a bound,
    - positive: a `x_block` reference, e.g. for qcones, or
    - negative: Negated index of a PSD variable in the underlying task
    """
    c_block_slack   :: Vector{Int}
    """
    One entry per variable-constraint block indicating which cone it belongs to, 0 if none.
    """
    c_coneid   :: Vector{Int}

    ###########################
    conecounter :: Int

    ###########################
    trm :: Union{Nothing, Rescode}
    solutions :: Vector{MosekSolution}

    ###########################
    """
    Indicating whether the objective sense is MOI.FEASIBILITY_SENSE. It is
    encoded as a MOI.MIN_SENSE with a zero objective internally but this allows
    MOI.get(::MosekModel, ::ObjectiveSense) to still return the right value
    """
    feasibility :: Bool

    fallback :: Union{String,Nothing}
    service  :: Union{String,Nothing}
end

function parametrized_task(be_quiet::Bool,
                           ipars::Dict{String, Integer},
                           dpars::Dict{String, AbstractFloat},
                           spars::Dict{String, AbstractString})
    task = maketask()
    try
        for (name, value) in ipars
            Mosek.putnaintparam(task, name, value)
        end
        for (name, value) in dpars
            Mosek.putnadouparam(task, name, value)
        end
        for (name, value) in spars
            Mosek.putnastrparam(task, name, value)
        end
        if !be_quiet
            Mosek.putstreamfunc(task, Mosek.MSK_STREAM_LOG, m -> print(m))
        end
    catch
        Mosek.deletetask(task)
        rethrow()
    end
    return task
end

function parse_parameters(kws)
    be_quiet = false
    fallback = nothing
    service   = nothing
    ipars = Dict{String, Integer}()
    dpars = Dict{String, AbstractFloat}()
    spars = Dict{String, AbstractString}()
    for (option, val) in kws
        parname = string(option)
        if parname == "QUIET"
            be_quiet = be_quiet || convert(Bool,val)
        elseif parname == "fallback"
            fallback = val
        elseif parname == "service"
            service = val
        elseif parname == "verbose"
            # pass
        elseif startswith(parname, "MSK_IPAR_")
            ipars[parname] = val
        elseif startswith(parname, "MSK_DPAR_")
            dpars[parname] = val
        elseif startswith(parname, "MSK_SPAR_")
            spars[parname] = val
        elseif isa(val, Integer)
            ipars["MSK_IPAR_$(uppercase(parname))"] = val
        elseif isa(val, AbstractFloat)
            dpars["MSK_DPAR_$(uppercase(parname))"] = val
        elseif isa(val, AbstractString)
            spars["MSK_SPAR_$(uppercase(parname))"] = val
        else
            error("Value $val for parameter $option has unrecognized type")
        end
    end
    return be_quiet,fallback,service,ipars, dpars, spars
end

export Mosek
function Mosek.Optimizer(; kws...)
    be_quiet, fallback, service,ipars, dpars, spars = parse_parameters(kws)
    return MosekModel(parametrized_task(be_quiet, ipars, dpars, spars), # task
                      be_quiet,
                      ipars,
                      dpars,
                      spars,
                      0, # public numvar
                      ConstraintMap(), # public constraints
                      Dict{String,MOI.ConstraintIndex}(),
                      LinkedInts(),# c_block
                      Int[], # x_boundflags
                      Int[], # x_numxc
                      LinkedInts(), # xc_block
                      UInt8[], # xc_bounds
                      Int[], # xc_coneid
                      Int[], # xc_idxs
                      LinkedInts(), # c_block
                      Float64[], # c_constant
                      Int[], # c_block_slack
                      Int[], # c_coneid
                      0, # cone counter
                      nothing,# trm
                      MosekSolution[],
                      true,# feasibility_sense
                      fallback,
                      service)
end

function MOI.optimize!(m::MosekModel)
    m.trm =
        if m.service === nothing
            if m.fallback === nothing
                optimize(m.task)
            else
                optimize(m.task,m.fallback)
            end
        else
            r = match(r"mosek://([a-zA-Z-]+(\.[a-zA-Z-]+)*)(:([0-9]+))?$",m.service)
            if r == nothing
                error("Invalid server specification $service")
            else
                server = r.captures[1]
                port   = if r.captures[4] == nothing "30080" else r.captures[4] end
                optimizermt(m.task,server,port)
            end
        end

    m.solutions = MosekSolution[]
    if solutiondef(m.task,MSK_SOL_ITG)
        push!(m.solutions,
              MosekSolution(MSK_SOL_ITG,
                            getsolsta(m.task,MSK_SOL_ITG),
                            getprosta(m.task,MSK_SOL_ITG),
                            getskx(m.task,MSK_SOL_ITG),
                            getxx(m.task,MSK_SOL_ITG),
                            Float64[],
                            Float64[],
                            Float64[],
                            getskc(m.task,MSK_SOL_ITG),
                            getxc(m.task,MSK_SOL_ITG),
                            Float64[],
                            Float64[],
                            Float64[]))
    end
    if solutiondef(m.task,MSK_SOL_BAS)
        push!(m.solutions,
              MosekSolution(MSK_SOL_BAS,
                            getsolsta(m.task,MSK_SOL_BAS),
                            getprosta(m.task,MSK_SOL_BAS),
                            getskx(m.task,MSK_SOL_BAS),
                            getxx(m.task,MSK_SOL_BAS),
                            getslx(m.task,MSK_SOL_BAS),
                            getsux(m.task,MSK_SOL_BAS),
                            Float64[],
                            getskc(m.task,MSK_SOL_BAS),
                            getxc(m.task,MSK_SOL_BAS),
                            getslc(m.task,MSK_SOL_BAS),
                            getsuc(m.task,MSK_SOL_BAS),
                            gety(m.task,MSK_SOL_BAS)))
    end
    if solutiondef(m.task,MSK_SOL_ITR)
        push!(m.solutions,
              MosekSolution(MSK_SOL_ITR,
                            getsolsta(m.task,MSK_SOL_ITR),
                            getprosta(m.task,MSK_SOL_ITR),
                            getskx(m.task,MSK_SOL_ITR),
                            getxx(m.task,MSK_SOL_ITR),
                            getslx(m.task,MSK_SOL_ITR),
                            getsux(m.task,MSK_SOL_ITR),
                            getsnx(m.task,MSK_SOL_ITR),
                            getskc(m.task,MSK_SOL_ITR),
                            getxc(m.task,MSK_SOL_ITR),
                            getslc(m.task,MSK_SOL_ITR),
                            getsuc(m.task,MSK_SOL_ITR),
                            gety(m.task,MSK_SOL_ITR)))
    end
end

function MOI.is_empty(m::MosekModel)
    getnumvar(m.task) == 0 && getnumcon(m.task) == 0 && getnumcone(m.task) == 0 && getnumbarvar(m.task) == 0
end

function MOI.empty!(model::MosekModel)
    model.task          = parametrized_task(model.be_quiet, model.ipars,
                                            model.dpars,    model.spars)
    model.publicnumvar  = 0
    model.constrmap     = ConstraintMap()
    model.x_block       = LinkedInts()
    model.x_boundflags  = Int[]
    model.x_numxc       = Int[]
    model.xc_block      = LinkedInts()
    model.xc_bounds     = UInt8[]
    model.xc_coneid     = Int[]
    model.xc_idxs       = Int[]
    model.c_block       = LinkedInts()
    model.c_constant    = Float64[]
    model.c_block_slack = Int[]
    model.c_coneid      = Int[]
    model.conecounter   = 0
    model.trm           = nothing
    model.solutions     = MosekSolution[]
    model.feasibility   = true
end

MOI.get(::MosekModel, ::MOI.SolverName) = "Mosek"

MOIU.supports_default_copy_to(::MosekModel, copy_names::Bool) = !copy_names
function MOI.copy_to(dest::MosekModel, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

function MOI.write_to_file(m::MosekModel, filename :: String)
    putintparam(m.task,MSK_IPAR_OPF_WRITE_SOLUTIONS, MSK_ON)
    writedata(m.task,filename)
end

# For linear objectives we accept:
# EITER affine left-hand side and ranged, unbounded, half-open, fixed (equality), PSD or SOC domains
# OR affine and quadratic left-hand side, and ranged, unbounded, half-open, fixed (equality) domains (quadratic constraints must be unbounded or half-open)
#
# For non-quadratic problems we allow binary and integer variables (but not constraints)
#function supportsconstraints(m::MosekSolver, constraint_types) :: Bool
#    for (fun,dom) in constraint_types
#        if  fun in [MOI.ScalarAffineFunction{Float64},
#                    MOI.SingleVariable,
#                    MOI.VectorAffineFunction{Float64},
#                    MOI.VectorOfVariables] &&
#            dom in [MOI.Zeros,
#                    MOI.Reals,
#                    MOI.Nonnegatives,
#                    MOI.Nonpositives,
#                    MOI.GreaterThan{Float64},
#                    MOI.LessThan{Float64},
#                    MOI.EqualTo{Float64},
#                    MOI.Interval{Float64},
#                    MOI.SecondOrderCone,
#                    MOI.RotatedSecondOrderCone,
#                    MOI.PositiveSemidefiniteConeTriangle,
#                    MOI.PositiveSemidefiniteConeScaled ]
#            # ok
#        elseif dom in [MOI.ZeroOne,
#                       MOI.Integer] &&
#                           fun in [MOI.SingleVariable,
#                                   MOI.VectorOfVariables]
#            # ok
#        else
#            return false
#        end
#    end
#    true
#end

#MOI.supportsproblem(m::MosekSolver, ::Type{MOI.SingleVariable},                constraint_types) :: Bool = supportsconstraints(m,constraint_types)
#MOI.supportsproblem(m::MosekSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, constraint_types) :: Bool = supportsconstraints(m,constraint_types)
#MOI.supportsproblem{F}(m::MosekSolver, ::Type{F}, constraint_types) :: Bool = false

ref2id(ref :: MOI.VariableIndex) :: Int = Int(ref.value)

ref2id(ref :: MOI.ConstraintIndex) :: Int =
    if ref.value & 1 == 0
        Int(ref.value >> 1)
    else
        - Int(ref.value >> 1)
    end

function id2vref(id :: Int) :: MOI.VariableIndex
    @assert(id > 0)
    MOI.VariableIndex(id)
end

include("objective.jl")
include("variable.jl")
include("constraint.jl")
include("attributes.jl")

export MosekModel

end # module


text:
export SwiftService

"""
Creates a swift service object, which describes the properties of your connection.
Default values for are read from enviroment variables, using the names from your `openstackrc.sh`.
Just like for the swift commandline client.

	- Note that `segment_size` is in bytes -- so `1024^3` is 1Gb
"""
function SwiftService(;auth_url = ENV["OS_AUTH_URL"],
                        username::AbstractString = ENV["OS_USERNAME"],
                        password::AbstractString = ENV["OS_PASSWORD"],
                        tenant_name::AbstractString = ENV["OS_TENANT_NAME"],
                        auth_version::AbstractString = "2",
                        segment_size::Integer = 1024^3, #1GB
                        segment_threads::Integer=24,
                        object_threads::Integer=24)
    service[:SwiftService](Dict(
        "os_auth_url"     => auth_url,
        "os_username"     => username,
        "os_password"     => password,
        "os_tenant_name"  => tenant_name,
        "auth_version"    => auth_version,
        "segment_size"    => segment_size,
        "segment_threads" => segment_threads,
        "object_threads"  => object_threads
    ))
end


text:
# Comparison of two (equivalent) ways of constructing a mechanism with an
# x, y, θ base: an explicit Planar joint vs. two Prismatic joints and one Rotational
# joint in series.

using LCPSim
using Base.Test
using RigidBodyDynamics
using RigidBodyDynamics: Bounds
using StaticArrays: SVector
using Cbc: CbcSolver
using Gurobi: GurobiSolver
using Rotations: RotY
using MeshCat
using MeshCatMechanisms

urdf = joinpath(@__DIR__, "urdf", "box.urdf")

function box_with_planar_base()
    mechanism = parse_urdf(Float64, urdf)
    core = findbody(mechanism, "core")
    fixed_joint = joint_to_parent(core, mechanism)
    floating_base = Joint(fixed_joint.name, frame_before(fixed_joint), frame_after(fixed_joint),
                          Planar([1., 0, 0], [0., 0, 1.]),
                          position_bounds=[Bounds(-5., 5), Bounds(0., 3), Bounds(-2π, 2π)],
                          velocity_bounds=[Bounds(-100., 10), Bounds(-100., 10), Bounds(-100., 100)],
                          effort_bounds=[Bounds(0., 0) for i in 1:3])
    replace_joint!(mechanism, fixed_joint, floating_base)

    x0 = MechanismState{Float64}(mechanism)
    set_velocity!(x0, zeros(num_velocities(x0)))
    set_configuration!(x0, floating_base, [-1, 0.3, -0.3])
    configuration_derivative_to_velocity!(velocity(x0, floating_base), floating_base, configuration(x0, floating_base), [3., 0, 0])
    mechanism, x0
end

function box_with_dummy_links()
    urdf_mech = parse_urdf(Float64, urdf)
    mechanism, base = planar_revolute_base()
    attach!(mechanism, base, urdf_mech)

    x0 = MechanismState{Float64}(mechanism)
    set_velocity!(x0, zeros(num_velocities(x0)))
    set_configuration!(x0, findjoint(mechanism, "base_x"), [-1])
    set_configuration!(x0, findjoint(mechanism, "base_z"), [0.3])
    set_configuration!(x0, findjoint(mechanism, "base_rotation"), [0.3])
    set_velocity!(x0, findjoint(mechanism, "base_x"), [3])
    mechanism, x0
end

function environment_with_floor(mechanism)
    world = root_body(mechanism)
    floor = planar_obstacle(default_frame(world), [0, 0, 1.], [0, 0, 0.], 0.5)

    core = findbody(mechanism, "core")
    env = Environment([
        (core, pt, floor) for pt in [
                    Point3D(default_frame(core), SVector(0.1, 0, 0.2)),
                    Point3D(default_frame(core), SVector(-0.1, 0, 0.2)),
                    Point3D(default_frame(core), SVector(0.1, 0, -0.2)),
                    Point3D(default_frame(core), SVector(-0.1, 0, -0.2)),
                     ]
            ])
    env
end

function withenv(f::Function)
    env = Gurobi.Env()
    try
        f(env)
    finally
        # https://github.com/JuliaOpt/Gurobi.jl/issues/110
        # Gurobi.free_env(env)
    end
end

@testset "planar vs dummy link brick" begin
    mech1, x1 = box_with_planar_base()
    env1 = environment_with_floor(mech1)

    mech2, x2 = box_with_dummy_links()
    env2 = environment_with_floor(mech2)

    controller = LCPSim.passive_controller()
    Δt = 0.01
    N = 600

    results1 = withenv() do env
        LCPSim.simulate(x1, controller, env1, Δt, N, GurobiSolver(env, OutputFlag=0))
    end
    results2 = withenv() do env
        LCPSim.simulate(x2, controller, env2, Δt, N, GurobiSolver(env, OutputFlag=0))
    end

    vis = Visualizer()
    if !haskey(ENV, "CI")
        open(vis)
        wait(vis)
    end

    mv1 = MechanismVisualizer(mech1, URDFVisuals(urdf), vis[:planar])
    mv2 = MechanismVisualizer(mech2, URDFVisuals(urdf), vis[:dummy])

    for i in 1:length(results1)
        set_configuration!(mv1, configuration(results1[i].state))
        set_configuration!(mv2, configuration(results2[i].state))
        sleep(Δt)
    end

    @test length(results1) == N
    @test length(results2) == N

    for i in 1:N
        set_configuration!(x1, configuration(results1[i].state))
        set_velocity!(x1, velocity(results1[i].state))
        set_configuration!(x2, configuration(results2[i].state))
        set_velocity!(x2, velocity(results2[i].state))
        T1 = transform_to_root(x1, findbody(mech1, "core"))
        T2 = transform_to_root(x2, findbody(mech2, "core"))
        atol = N <= 30 ? 1e-2 : 5e-2
        @test rotation(T1) ≈ rotation(T2) atol=atol
        @test translation(T1) ≈ translation(T2) atol=atol
    end
end


text:
using Revise
using VCFTools
using MendelImpute
using GeneticVariation
using Random
using StatsBase
using CodecZlib
using ProgressMeter
using JLD2, FileIO, JLSO
using BenchmarkTools
using GroupSlices
using LinearAlgebra
using Profile

chr = 20
maf = 0.01
width = 64

BLAS.set_num_threads(1)
Threads.nthreads() != 10 && error("not 10 threads!")
Profile.init(n = 10^7, delay = 1.0)

Random.seed!(2020)
tgtfile = "target.chr$chr.typedOnly.maf$maf.masked.vcf.gz"
reffile = "./single_precision/ref.chr$chr.w$width.maf$maf.excludeTarget.jlso"
outfile = "mendel.imputed.dp$width.maf$maf.vcf.gz"

phase(tgtfile, reffile, outfile=outfile, impute=true, width=width,
    dynamic_programming = false)
Profile.clear_malloc_data()
phase(tgtfile, reffile, outfile=outfile, impute=true, width=width,
    dynamic_programming = false)

rm(outfile, force=true)

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using Documenter, PowerSystems, DocStringExtensions, PowerSimulationsDynamics

makedocs(
    modules = [PowerSimulationsDynamics],
    format = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename = "PowerSimulationsDynamics.jl",
    pages = Any[
        "Welcome Page" => "index.md",
        "Quick Start Guide" => "quick_start_guide.md",
        "Tutorials" => "tutorials_page.md",
        "Models" => "models.md",
        "Initialization" => "initialization.md",
        "Small Signal" => "small.md",
        "Reference Frames" => "reference_frames.md",
        "Generator Component Library" => Any[
            "Machine" => "component_models/machines.md",
            "Shaft" => "component_models/shafts.md",
            "AVR" => "component_models/avr.md",
            "PSS" => "component_models/pss.md",
            "Turbine and Governor" => "component_models/turbine_gov.md",
        ],
        "CIG Component Library" => Any[
            "Converter" => "component_models/converter.md",
            "DC Sources" => "component_models/dc_source.md",
            "Filter" => "component_models/filters.md",
            "Frequency Estimators" => "component_models/freq_esti.md",
            "Inner Control" => "component_models/inner_control.md",
            "Outer Control" => "component_models/outer_control.md",
        ],
        "Branch Models" => Any["Network" => "component_models/network.md",],
        "Code Base Developer Guide" => Any["Developer Guide" => "code_base_developer_guide/developer.md",],
        "Public API Reference" => "api/public.md",
        "Internal API Reference" => "api/internal.md",
    ],
)

deploydocs(
    repo = "github.com/NREL-SIIP/PowerSimulationsDynamics.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    push_preview = true,
    versions = ["stable" => "v^", "v#.#"],
)


text:
#######################################################################
#
# An example of creating an in-memory XLS file with XlsxWriter.

using XlsxWriter

function test()

	buff = IOBuffer()
	wb = Workbook(buff)
	ws = add_worksheet!(wb)
	write!(ws, 0, 0, "In memory")
	close(wb)

	# dump it to a file for testing but really for throwing across the wire via HTTP or whatever

	fid = open("in_memory.xlsx", "w+")
	write(fid, take!(buff))
	close(fid)

	isfile("in_memory.xlsx")
end

test()


text:
function GetPath()
    path=readlines(STDIN)
    return path
end

type TMResults
    Dr::Array{Float64,2}
    Dt::Array{Float64,2}
    Nr::Array{Float64,2}
    Rec::Array{Float64,2}
    Prec::Array{Float64,2}
    RecFinal::Array{Float64,1}
    PrecFinal::Array{Float64,1}
    TMResults()=new()
end

type NewDataCorpus
    D_C::Corpus
    Q_C::Corpus
    NewDataCorpus()=new()
end

type NewDataMatrix
    A::Array{Float64,2}
    NumQueries::Int
    NewDataMatrix()=new()
end

#Initializes the Results type with zeros
function InitResults!(R::TMResults,nq::Int,lt::Int)
    R.Dr=zeros(nq,lt)
    R.Dt=zeros(nq,lt)
    R.Nr=zeros(nq,lt)
    R.Rec=zeros(nq,lt)
    R.Prec=zeros(nq,lt)
    R.RecFinal=zeros(lt)
    R.PrecFinal=zeros(lt)
    return R
end

function Normalize(A::Array{Float64,2})
    At=sqrt(diag(A'*A))
    A_norm=A/(diagm(At))
    return A_norm
end

function CosTheta(q::Array{Float64,1},d::Array{Float64,1})
    return q'*d/(norm(q,2)*norm(d,2))
end

function FindResults(k::Int,z::Int,Costheta::Array{Float64,2},tol::Float64)
    function testf1(x)
        x>tol
    end
    ii=find(testf1,Costheta)
    function testf2(x)
        x==z
    end
    b=find(testf2,T[:,1])
    dr=length(intersect(ii,T[b,2]))
    dt=length(ii)
    nr=length(b)
    Reca = (dr/nr)*100
    Preca = (dr/dt)*100
    return dr,dt,nr,Reca,Preca
end


text:
module ChemicalKinetics
  using Distributions
  using Sundials

  include("ode_model.jl")
  include("ode_solver.jl")
  include("parser.jl")
  include("symbolic.jl")

  export
    NSE,
    SE,
    OdeModel,
    CkCvode,
    parse_model,
    replace
end


text:
score_api = "http://api.juliatetris.com"

rounds_endpoint = "$(score_api)/rounds"

scores_endpoint = "$(score_api)/scores"

json_headers = Dict(
  "Content-Type" => "application/json",
  "Accept" => "application/json"
)

gzip_headers = Dict(
  "Content-Type" => "application/gzip",
  "Accept" => "application/json"
)

api_errors = [
  HTTP.ExceptionRequest.StatusError,
  Base.UVError
]

function api_fetch(cur_response::HTTP.Response, cur_key::Union{AbstractString, Int, Void}=nothing)
  cur_fetch = JSON.parse(
    String(cur_response.body)
  )

  ( cur_key == nothing ) ||
    ( cur_fetch = cur_fetch[cur_key] )

  cur_fetch
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using ACE, JuLIP, BenchmarkTools, LinearAlgebra

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N), 0
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

trans = PolyTransform(2, 1.0)
fcut = PolyCutoff2s(2, 0.5, 3.0)
BB = [ SHIPBasis(SparseSHIP(2, 20, wL=1.0), trans, fcut),
       SHIPBasis(SparseSHIP(3, 16, wL=1.5), trans, fcut),
       SHIPBasis(SparseSHIP(4, 12, wL=1.5), trans, fcut),
       SHIPBasis(SparseSHIP(5, 10, wL=1.5), trans, fcut) ]

Nat = 30
Rs, Zs, z0 = randR(Nat)
btmp = ACE.alloc_temp(BB[1])
@info("profile precomputation of A")
@btime ACE.precompute_A!($btmp, $(BB[1]), $Rs, $Zs, 1)

@info("profile ship-basis and fast-ship site energies")
for n = 2:5
   @info("  body-order $(n+1):")
   Rs, Zs, z0 = randR(Nat)
   B = BB[n-1]
   coeffs = randcoeffs(B)
   🚢 = SHIP(B, coeffs)
   b = ACE.alloc_B(B)
   btmp = ACE.alloc_temp(B, length(Rs))
   tmp = ACE.alloc_temp(🚢, length(Rs))
   @info("     evaluate a site energy:")
   print("         SHIPBasis: "); @btime evaluate!($b, $btmp, $B, $Rs, $Zs, $z0)
   print("         SHIP     : "); @btime evaluate!($tmp, $🚢, $Rs, $Zs, $z0)

   dEs = zeros(JVecF, length(Rs))
   db = ACE.alloc_dB(B, length(Rs))
   dbtmp = ACE.alloc_temp_d(B, length(Rs))
   tmp = ACE.alloc_temp_d(🚢, length(Rs))

   @info("     site energy gradient:")
   print("         SHIPBasis: "); @btime evaluate_d!($db, $dbtmp, $B, $Rs, $Zs, $z0)
   print("         SHIP     : "); @btime evaluate_d!($dEs, $tmp, $🚢, $Rs, $Zs, $z0)
end