text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## RandomDevice

# SamplerUnion(X, Y, ...}) == Union{SamplerType{X}, SamplerType{Y}, ...}
SamplerUnion(U...) = Union{Any[SamplerType{T} for T in U]...}
const SamplerBoolBitInteger = SamplerUnion(Bool, BitInteger_types...)

if Sys.iswindows()
    struct RandomDevice <: AbstractRNG
        buffer::Vector{UInt128}

        RandomDevice() = new(Vector{UInt128}(undef, 1))
    end

    function rand(rd::RandomDevice, sp::SamplerBoolBitInteger)
        rand!(rd, rd.buffer)
        @inbounds return rd.buffer[1] % sp[]
    end
else # !windows
    struct RandomDevice <: AbstractRNG
        file::IOStream
        unlimited::Bool

        RandomDevice(; unlimited::Bool=true) =
            new(open(unlimited ? "/dev/urandom" : "/dev/random"), unlimited)
    end

    rand(rd::RandomDevice, sp::SamplerBoolBitInteger) = read( rd.file, sp[])

    function serialize(s::AbstractSerializer, rd::RandomDevice)
        Serialization.serialize_type(s, typeof(rd))
        serialize(s, rd.unlimited)
    end
    function deserialize(s::AbstractSerializer, t::Type{RandomDevice})
        unlimited = deserialize(s)
        return RandomDevice(unlimited=unlimited)
    end

end # os-test

# NOTE: this can't be put within the if-else block above
for T in (Bool, BitInteger_types...)
    if Sys.iswindows()
        @eval function rand!(rd::RandomDevice, A::Array{$T}, ::SamplerType{$T})
            ccall((:SystemFunction036, :Advapi32), stdcall, UInt8, (Ptr{Cvoid}, UInt32),
                  A, sizeof(A))
            A
        end
    else
        @eval rand!(rd::RandomDevice, A::Array{$T}, ::SamplerType{$T}) = read!(rd.file, A)
    end
end

# RandomDevice produces natively UInt64
rng_native_52(::RandomDevice) = UInt64

"""
    RandomDevice()

Create a `RandomDevice` RNG object.
Two such objects will always generate different streams of random numbers.
The entropy is obtained from the operating system.
"""
RandomDevice

RandomDevice(::Nothing) = RandomDevice()
seed!(rng::RandomDevice) = rng

## MersenneTwister

const MT_CACHE_F = 501 << 1 # number of Float64 in the cache
const MT_CACHE_I = 501 << 4 # number of bytes in the UInt128 cache

@assert dsfmt_get_min_array_size() <= MT_CACHE_F

mutable struct MersenneTwister <: AbstractRNG
    seed::Vector{UInt32}
    state::DSFMT_state
    vals::Vector{Float64}
    ints::Vector{UInt128}
    idxF::Int
    idxI::Int

    function MersenneTwister(seed, state, vals, ints, idxF, idxI)
        length(vals) == MT_CACHE_F && 0 <= idxF <= MT_CACHE_F ||
            throw(DomainError((length(vals), idxF),
                      "`length(vals)` and `idxF` must be consistent with $MT_CACHE_F"))
        length(ints) == MT_CACHE_I >> 4 && 0 <= idxI <= MT_CACHE_I ||
            throw(DomainError((length(ints), idxI),
                      "`length(ints)` and `idxI` must be consistent with $MT_CACHE_I"))
        new(seed, state, vals, ints, idxF, idxI)
    end
end

MersenneTwister(seed::Vector{UInt32}, state::DSFMT_state) =
    MersenneTwister(seed, state,
                    Vector{Float64}(undef, MT_CACHE_F),
                    Vector{UInt128}(undef, MT_CACHE_I >> 4),
                    MT_CACHE_F, 0)

"""
    MersenneTwister(seed)
    MersenneTwister()

Create a `MersenneTwister` RNG object. Different RNG objects can have
their own seeds, which may be useful for generating different streams
of random numbers.
The `seed` may be a non-negative integer or a vector of
`UInt32` integers. If no seed is provided, a randomly generated one
is created (using entropy from the system).
See the [`seed!`](@ref) function for reseeding an already existing
`MersenneTwister` object.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> x1 = rand(rng, 2)
2-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592

julia> rng = MersenneTwister(1234);

julia> x2 = rand(rng, 2)
2-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592

julia> x1 == x2
true
```
"""
MersenneTwister(seed=nothing) =
    seed!(MersenneTwister(Vector{UInt32}(), DSFMT_state()), seed)

function copy!(dst::MersenneTwister, src::MersenneTwister)
    copyto!(resize!(dst.seed, length(src.seed)), src.seed)
    copy!(dst.state, src.state)
    copyto!(dst.vals, src.vals)
    copyto!(dst.ints, src.ints)
    dst.idxF = src.idxF
    dst.idxI = src.idxI
    dst
end

copy(src::MersenneTwister) =
    MersenneTwister(copy(src.seed), copy(src.state), copy(src.vals), copy(src.ints),
                    src.idxF, src.idxI)

==(r1::MersenneTwister, r2::MersenneTwister) =
    r1.seed == r2.seed && r1.state == r2.state &&
    isequal(r1.vals, r2.vals) &&
    isequal(r1.ints, r2.ints) &&
    r1.idxF == r2.idxF && r1.idxI == r2.idxI

hash(r::MersenneTwister, h::UInt) =
    foldr(hash, (r.seed, r.state, r.vals, r.ints, r.idxF, r.idxI); init=h)

function fillcache_zeros!(r::MersenneTwister)
    # the use of this function is not strictly necessary, but it makes
    # comparing two MersenneTwister RNGs easier
    fill!(r.vals, 0.0)
    fill!(r.ints, zero(UInt128))
    r
end

### low level API

#### floats

mt_avail(r::MersenneTwister) = MT_CACHE_F - r.idxF
mt_empty(r::MersenneTwister) = r.idxF == MT_CACHE_F
mt_setfull!(r::MersenneTwister) = r.idxF = 0
mt_setempty!(r::MersenneTwister) = r.idxF = MT_CACHE_F
mt_pop!(r::MersenneTwister) = @inbounds return r.vals[r.idxF+=1]

function gen_rand(r::MersenneTwister)
    GC.@preserve r dsfmt_fill_array_close1_open2!(r.state, pointer(r.vals), length(r.vals))
    mt_setfull!(r)
end

reserve_1(r::MersenneTwister) = (mt_empty(r) && gen_rand(r); nothing)
# `reserve` allows one to call `rand_inbounds` n times
# precondition: n <= MT_CACHE_F
reserve(r::MersenneTwister, n::Int) = (mt_avail(r) < n && gen_rand(r); nothing)

#### ints

logsizeof(::Type{<:Union{Bool,Int8,UInt8}}) = 0
logsizeof(::Type{<:Union{Int16,UInt16}}) = 1
logsizeof(::Type{<:Union{Int32,UInt32}}) = 2
logsizeof(::Type{<:Union{Int64,UInt64}}) = 3
logsizeof(::Type{<:Union{Int128,UInt128}}) = 4

idxmask(::Type{<:Union{Bool,Int8,UInt8}}) = 15
idxmask(::Type{<:Union{Int16,UInt16}}) = 7
idxmask(::Type{<:Union{Int32,UInt32}}) = 3
idxmask(::Type{<:Union{Int64,UInt64}}) = 1
idxmask(::Type{<:Union{Int128,UInt128}}) = 0

mt_avail(r::MersenneTwister, ::Type{T}) where {T<:BitInteger} =
    r.idxI >> logsizeof(T)

function mt_setfull!(r::MersenneTwister, ::Type{<:BitInteger})
    rand!(r, r.ints)
    r.idxI = MT_CACHE_I
end

mt_setempty!(r::MersenneTwister, ::Type{<:BitInteger}) = r.idxI = 0

function reserve1(r::MersenneTwister, ::Type{T}) where T<:BitInteger
    r.idxI < sizeof(T) && mt_setfull!(r, T)
    nothing
end

function mt_pop!(r::MersenneTwister, ::Type{T}) where T<:BitInteger
    reserve1(r, T)
    r.idxI -= sizeof(T)
    i = r.idxI
    @inbounds x128 = r.ints[1 + i >> 4]
    i128 = (i >> logsizeof(T)) & idxmask(T) # 0-based "indice" in x128
    (x128 >> (i128 * (sizeof(T) << 3))) % T
end

# not necessary, but very slightly more efficient
function mt_pop!(r::MersenneTwister, ::Type{T}) where {T<:Union{Int128,UInt128}}
    reserve1(r, T)
    @inbounds res = r.ints[r.idxI >> 4]
    r.idxI -= 16
    res % T
end

### seeding

#### make_seed()

# make_seed produces values of type Vector{UInt32}, suitable for MersenneTwister seeding
function make_seed()
    try
        return rand(RandomDevice(), UInt32, 4)
    catch
        println(stderr,
                "Entropy pool not available to seed RNG; using ad-hoc entropy sources.")
        seed = reinterpret(UInt64, time())
        seed = hash(seed, UInt64(getpid()))
        try
            seed = hash(seed, parse(UInt64,
                                    read(pipeline(`ifconfig`, `sha1sum`), String)[1:40],
                                    base = 16))
        catch
        end
        return make_seed(seed)
    end
end

function make_seed(n::Integer)
    n < 0 && throw(DomainError(n, "`n` must be non-negative."))
    seed = UInt32[]
    while true
        push!(seed, n & 0xffffffff)
        n >>= 32
        if n == 0
            return seed
        end
    end
end

#### seed!()

function seed!(r::MersenneTwister, seed::Vector{UInt32})
    copyto!(resize!(r.seed, length(seed)), seed)
    dsfmt_init_by_array(r.state, r.seed)
    mt_setempty!(r)
    mt_setempty!(r, UInt128)
    fillcache_zeros!(r)
    return r
end

seed!(r::MersenneTwister=default_rng()) = seed!(r, make_seed())
seed!(r::MersenneTwister, n::Integer) = seed!(r, make_seed(n))
seed!(seed::Union{Integer,Vector{UInt32}}) = seed!(default_rng(), seed)

### Global RNG

const THREAD_RNGs = MersenneTwister[]
@inline default_rng() = default_rng(Threads.threadid())
@noinline function default_rng(tid::Int)
    @assert 0 < tid <= length(THREAD_RNGs)
    if @inbounds isassigned(THREAD_RNGs, tid)
        @inbounds MT = THREAD_RNGs[tid]
    else
        MT = MersenneTwister()
        @inbounds THREAD_RNGs[tid] = MT
    end
    return MT
end
function __init__()
    resize!(empty!(THREAD_RNGs), Threads.nthreads()) # ensures that we didn't save a bad object
end

struct _GLOBAL_RNG <: AbstractRNG
    global const GLOBAL_RNG = _GLOBAL_RNG.instance
end

copy!(dst::MersenneTwister, ::_GLOBAL_RNG) = copy!(dst, default_rng())
copy!(::_GLOBAL_RNG, src::MersenneTwister) = copy!(default_rng(), src)
copy(::_GLOBAL_RNG) = copy(default_rng())

seed!(::_GLOBAL_RNG, seed::Vector{UInt32}) = seed!(default_rng(), seed)
seed!(::_GLOBAL_RNG, n::Integer) = seed!(default_rng(), n)
seed!(::_GLOBAL_RNG, ::Nothing) = seed!(default_rng(), nothing)
seed!(::_GLOBAL_RNG) = seed!(default_rng(), nothing)

rng_native_52(::_GLOBAL_RNG) = rng_native_52(default_rng())
rand(::_GLOBAL_RNG, sp::SamplerBoolBitInteger) = rand(default_rng(), sp)
for T in (:(SamplerTrivial{UInt52Raw{UInt64}}),
          :(SamplerTrivial{UInt2x52Raw{UInt128}}),
          :(SamplerTrivial{UInt104Raw{UInt128}}),
          :(SamplerTrivial{CloseOpen12_64}),
          :(SamplerUnion(Int64, UInt64, Int128, UInt128)),
          :(SamplerUnion(Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32)),
         )
    @eval rand(::_GLOBAL_RNG, x::$T) = rand(default_rng(), x)
end

rand!(::_GLOBAL_RNG, A::AbstractArray{Float64}, I::SamplerTrivial{<:FloatInterval_64}) = rand!(default_rng(), A, I)
rand!(::_GLOBAL_RNG, A::Array{Float64}, I::SamplerTrivial{<:FloatInterval_64}) = rand!(default_rng(), A, I)
for T in (Float16, Float32)
    @eval rand!(::_GLOBAL_RNG, A::Array{$T}, I::SamplerTrivial{CloseOpen12{$T}}) = rand!(default_rng(), A, I)
    @eval rand!(::_GLOBAL_RNG, A::Array{$T}, I::SamplerTrivial{CloseOpen01{$T}}) = rand!(default_rng(), A, I)
end
for T in BitInteger_types
    @eval rand!(::_GLOBAL_RNG, A::Array{$T}, I::SamplerType{$T}) = rand!(default_rng(), A, I)
end

### generation

# MersenneTwister produces natively Float64
rng_native_52(::MersenneTwister) = Float64

#### helper functions

# precondition: !mt_empty(r)
rand_inbounds(r::MersenneTwister, ::CloseOpen12_64) = mt_pop!(r)
rand_inbounds(r::MersenneTwister, ::CloseOpen01_64=CloseOpen01()) =
    rand_inbounds(r, CloseOpen12()) - 1.0

rand_inbounds(r::MersenneTwister, ::UInt52Raw{T}) where {T<:BitInteger} =
    reinterpret(UInt64, rand_inbounds(r, CloseOpen12())) % T

function rand(r::MersenneTwister, x::SamplerTrivial{UInt52Raw{UInt64}})
    reserve_1(r)
    rand_inbounds(r, x[])
end

function rand(r::MersenneTwister, ::SamplerTrivial{UInt2x52Raw{UInt128}})
    reserve(r, 2)
    rand_inbounds(r, UInt52Raw(UInt128)) << 64 | rand_inbounds(r, UInt52Raw(UInt128))
end

function rand(r::MersenneTwister, ::SamplerTrivial{UInt104Raw{UInt128}})
    reserve(r, 2)
    rand_inbounds(r, UInt52Raw(UInt128)) << 52 ⊻ rand_inbounds(r, UInt52Raw(UInt128))
end

#### floats

rand(r::MersenneTwister, sp::SamplerTrivial{CloseOpen12_64}) =
    (reserve_1(r); rand_inbounds(r, sp[]))

#### integers

rand(r::MersenneTwister, T::SamplerUnion(Int64, UInt64, Int128, UInt128)) =
    mt_pop!(r, T[])

rand(r::MersenneTwister, T::SamplerUnion(Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32)) =
    rand(r, UInt52Raw()) % T[]

#### arrays of floats

##### AbstractArray

function rand!(r::MersenneTwister, A::AbstractArray{Float64},
               I::SamplerTrivial{<:FloatInterval_64})
    region = LinearIndices(A)
    # what follows is equivalent to this simple loop but more efficient:
    # for i=region
    #     @inbounds A[i] = rand(r, I[])
    # end
    m = Base.checked_sub(first(region), 1)
    n = last(region)
    while m < n
        s = mt_avail(r)
        if s == 0
            gen_rand(r)
            s = mt_avail(r)
        end
        m2 = min(n, m+s)
        for i=m+1:m2
            @inbounds A[i] = rand_inbounds(r, I[])
        end
        m = m2
    end
    A
end

##### Array : internal functions

# internal array-like type to circumevent the lack of flexibility with reinterpret
struct UnsafeView{T} <: DenseArray{T,1}
    ptr::Ptr{T}
    len::Int
end

Base.length(a::UnsafeView) = a.len
Base.getindex(a::UnsafeView, i::Int) = unsafe_load(a.ptr, i)
Base.setindex!(a::UnsafeView, x, i::Int) = unsafe_store!(a.ptr, x, i)
Base.pointer(a::UnsafeView) = a.ptr
Base.size(a::UnsafeView) = (a.len,)
Base.elsize(::UnsafeView{T}) where {T} = sizeof(T)

# this is essentially equivalent to rand!(r, ::AbstractArray{Float64}, I) above, but due to
# optimizations which can't be done currently when working with pointers, we have to re-order
# manually the computation flow to get the performance
# (see https://discourse.julialang.org/t/unsafe-store-sometimes-slower-than-arrays-setindex)
function _rand_max383!(r::MersenneTwister, A::UnsafeView{Float64}, I::FloatInterval_64)
    n = length(A)
    @assert n <= dsfmt_get_min_array_size()+1 # == 383
    mt_avail(r) == 0 && gen_rand(r)
    # from now on, at most one call to gen_rand(r) will be necessary
    m = min(n, mt_avail(r))
    GC.@preserve r unsafe_copyto!(A.ptr, pointer(r.vals, r.idxF+1), m)
    if m == n
        r.idxF += m
    else # m < n
        gen_rand(r)
        GC.@preserve r unsafe_copyto!(A.ptr+m*sizeof(Float64), pointer(r.vals), n-m)
        r.idxF = n-m
    end
    if I isa CloseOpen01
        for i=1:n
            A[i] -= 1.0
        end
    end
    A
end

fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::CloseOpen01_64) =
    dsfmt_fill_array_close_open!(s, A, n)

fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::CloseOpen12_64) =
    dsfmt_fill_array_close1_open2!(s, A, n)

function rand!(r::MersenneTwister, A::UnsafeView{Float64},
               I::SamplerTrivial{<:FloatInterval_64})
    # depending on the alignment of A, the data written by fill_array! may have
    # to be left-shifted by up to 15 bytes (cf. unsafe_copyto! below) for
    # reproducibility purposes;
    # so, even for well aligned arrays, fill_array! is used to generate only
    # the n-2 first values (or n-3 if n is odd), and the remaining values are
    # generated by the scalar version of rand
    n = length(A)
    n2 = (n-2) ÷ 2 * 2
    n2 < dsfmt_get_min_array_size() && return _rand_max383!(r, A, I[])

    pA = A.ptr
    align = Csize_t(pA) % 16
    if align > 0
        pA2 = pA + 16 - align
        fill_array!(r.state, pA2, n2, I[]) # generate the data in-place, but shifted
        unsafe_copyto!(pA, pA2, n2) # move the data to the beginning of the array
    else
        fill_array!(r.state, pA, n2, I[])
    end
    for i=n2+1:n
        A[i] = rand(r, I[])
    end
    A
end

# fills up A reinterpreted as an array of Float64 with n64 values
function _rand!(r::MersenneTwister, A::Array{T}, n64::Int, I::FloatInterval_64) where T
    # n64 is the length in terms of `Float64` of the target
    @assert sizeof(Float64)*n64 <= sizeof(T)*length(A) && isbitstype(T)
    GC.@preserve A rand!(r, UnsafeView{Float64}(pointer(A), n64), SamplerTrivial(I))
    A
end

##### Array: Float64, Float16, Float32

rand!(r::MersenneTwister, A::Array{Float64}, I::SamplerTrivial{<:FloatInterval_64}) =
    _rand!(r, A, length(A), I[])

mask128(u::UInt128, ::Type{Float16}) =
    (u & 0x03ff03ff03ff03ff03ff03ff03ff03ff) | 0x3c003c003c003c003c003c003c003c00

mask128(u::UInt128, ::Type{Float32}) =
    (u & 0x007fffff007fffff007fffff007fffff) | 0x3f8000003f8000003f8000003f800000

for T in (Float16, Float32)
    @eval function rand!(r::MersenneTwister, A::Array{$T}, ::SamplerTrivial{CloseOpen12{$T}})
        n = length(A)
        n128 = n * sizeof($T) ÷ 16
        _rand!(r, A, 2*n128, CloseOpen12())
        GC.@preserve A begin
            A128 = UnsafeView{UInt128}(pointer(A), n128)
            for i in 1:n128
                u = A128[i]
                u ⊻= u << 26
                # at this point, the 64 low bits of u, "k" being the k-th bit of A128[i] and "+"
                # the bit xor, are:
                # [..., 58+32,..., 53+27, 52+26, ..., 33+7, 32+6, ..., 27+1, 26, ..., 1]
                # the bits needing to be random are
                # [1:10, 17:26, 33:42, 49:58] (for Float16)
                # [1:23, 33:55] (for Float32)
                # this is obviously satisfied on the 32 low bits side, and on the high side,
                # the entropy comes from bits 33:52 of A128[i] and then from bits 27:32
                # (which are discarded on the low side)
                # this is similar for the 64 high bits of u
                A128[i] = mask128(u, $T)
            end
        end
        for i in 16*n128÷sizeof($T)+1:n
            @inbounds A[i] = rand(r, $T) + one($T)
        end
        A
    end

    @eval function rand!(r::MersenneTwister, A::Array{$T}, ::SamplerTrivial{CloseOpen01{$T}})
        rand!(r, A, CloseOpen12($T))
        I32 = one(Float32)
        for i in eachindex(A)
            @inbounds A[i] = Float32(A[i])-I32 # faster than "A[i] -= one(T)" for T==Float16
        end
        A
    end
end

#### arrays of integers

function rand!(r::MersenneTwister, A::UnsafeView{UInt128}, ::SamplerType{UInt128})
    n::Int=length(A)
    i = n
    while true
        rand!(r, UnsafeView{Float64}(A.ptr, 2i), CloseOpen12())
        n < 5 && break
        i = 0
        while n-i >= 5
            u = A[i+=1]
            A[n]    ⊻= u << 48
            A[n-=1] ⊻= u << 36
            A[n-=1] ⊻= u << 24
            A[n-=1] ⊻= u << 12
            n-=1
        end
    end
    if n > 0
        u = rand(r, UInt2x52Raw())
        for i = 1:n
            A[i] ⊻= u << (12*i)
        end
    end
    A
end

for T in BitInteger_types
    @eval rand!(r::MersenneTwister, A::Array{$T}, sp::SamplerType{$T}) =
        (GC.@preserve A rand!(r, UnsafeView(pointer(A), length(A)), sp); A)

    T == UInt128 && continue

    @eval function rand!(r::MersenneTwister, A::UnsafeView{$T}, ::SamplerType{$T})
        n = length(A)
        n128 = n * sizeof($T) ÷ 16
        rand!(r, UnsafeView{UInt128}(pointer(A), n128))
        for i = 16*n128÷sizeof($T)+1:n
            @inbounds A[i] = rand(r, $T)
        end
        A
    end
end

#### from a range

for T in BitInteger_types, R=(1, Inf) # eval because of ambiguity otherwise
    @eval Sampler(::Type{MersenneTwister}, r::UnitRange{$T}, ::Val{$R}) =
        SamplerRangeFast(r)
end

### randjump

# Old randjump methods are deprecated, the scalar version is in the Future module.

_randjump(r::MersenneTwister, jumppoly::DSFMT.GF2X) =
    fillcache_zeros!(MersenneTwister(copy(r.seed), DSFMT.dsfmt_jump(r.state, jumppoly)))

max_stars_repo_name: mschauer/julia


text:
using BenchmarkTools, FinancialMonteCarlo, CUDA, ReverseDiff

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 100000;
Nstep = 30;
sigma = 0.2
mc = MonteCarloConfiguration(Nsim, Nstep);

mc_2 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.CudaMode());
toll = 1e-3;

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

f(x) = pricer(BlackScholesProcess(x[1], Underlying(x[2], x[4])), ZeroRate(x[3]), mc, EuropeanOption(x[5], K), FinancialMonteCarlo.CudaMode());
x = Float64[sigma, S0, r, d, T]
g = x -> ReverseDiff.gradient(f, x);
y0 = g(x);
@btime f(x);
@btime g(x);

f_(x::Vector) = pricer(BlackScholesProcess(x[1], Underlying(x[2], x[4])), ZeroRate(x[3]), mc, AmericanOption(x[5], K), FinancialMonteCarlo.CudaMode());
g_ = x -> ReverseDiff.gradient(f_, x);
y0 = g(x);
@btime f_(x);
@btime g_(x);

ext: jl


text:
################################################################################
#
#   Struct for holding all parameters associated with an md simulation
#
################################################################################
abstract MDTrial

################################################################################
#
#   AT SOME POINT, IT WOULD BE GOOD TO HAVE A MACRO FOR AUTOMATICALLY
#   POPULATING MolecularDynamicsTrial TYPE FIELDS
#
################################################################################

################################################################################
#
#   Struct for holding all parameters associated with an md verlet simulation
#
################################################################################
type MolecularDynamicsTrial <: MDTrial
    # Parameters
    numBodies::Int64                # Number of particles
    ρ::Float64                      # Density of particles
    Td::Float64                     # Desired temperature
    dx::Float64                     # Time step; defaults to Verlet's value
    V::Float64                      # Volume of box
    L::Float64                      # Side length of box
    steps::Int64                    # How long to run the simulation
    σ::Float64                      # Length scale, from Verlet, in m
    ϵ::Float64                      # Energy scale, in Kelvin
    m::Float64                      # Argon atom mass, in kg
    τ::Float64                      # Time scale, in seconds

    # Dynamic Quantities
    y::Array{Float64, 3}            # Positions
    v::Array{Float64, 3}            # Velocities
    f::Array{Float64, 3}            # Accelerations
    pe::Array{Float64, 2}           # Potential Energies
    ke::Array{Float64, 2}           # Kinetic Energies
    pet::Array{Float64, 1}          # Total Potential Energy
    ket::Array{Float64, 1}          # Total Kinetic Energy
    T::Array{Float64, 1}            # Temperature at each step
    e::Array{Float64, 1}            # Total Energy
    vir::Array{Float64, 1}          # Virial, ∑∑r(∂V/∂r)
    currentStep::Int64              # Where this simulation left off
end

################################################################################
#
#   Struct for holding all parameters associated with an md metropolis
#   simulation
#
################################################################################
type MDMetropolisTrial <: MDTrial
    # Parameters
    numBodies::Int64                # Number of particles
    ρ::Float64                      # Density of particles
    Td::Float64                     # Desired temperature
    # dx::Float64                     # Time step; defaults to Verlet's value
    V::Float64                      # Volume of box
    L::Float64                      # Side length of box
    steps::Int64                    # How long to run the simulation
    σ::Float64                      # Length scale, from Verlet, in m
    ϵ::Float64                      # Energy scale, in Kelvin
    m::Float64                      # Argon atom mass, in kg
    τ::Float64                      # Time scale, in seconds

    # Dynamic Quantities
    y::Array{Float64, 3}            # Positions
    # v::Array{Float64, 3}            # Velocities
    # f::Array{Float64, 3}            # Accelerations
    pe::Array{Float64, 2}           # Potential Energies
    ke::Array{Float64, 2}           # Kinetic Energies
    pet::Array{Float64, 1}          # Total Potential Energy
    ket::Array{Float64, 1}          # Total Kinetic Energy
    # T::Array{Float64, 1}            # Temperature at each step
    e::Array{Float64, 1}            # Total Energy
    vir::Array{Float64, 1}          # Virial, ∑∑r(∂V/∂r)
    currentStep::Int64              # Where this simulation left off
end

# Outer constructor just for input parameters
MolecularDynamicsTrial( numBodies, ρ, Td, dx, steps, σ, ϵ, m ) =
    MolecularDynamicsTrial(
        numBodies,
        ρ,
        Td,
        dx,
        ( numBodies / ρ ),                      # Volume
        ( numBodies / ρ )^( 1 / 3 ),            # Side length
        steps,
        σ,                                      # Length scale, in m
        ϵ,                                      # Energy scale, in Kelvin
        m,                                      # Argon atom mass, in kg
        sqrt( m * σ^2 / ( 48ϵ * 1.3806e-23 ) ), # Time scale, in seconds
        zeros( 3, numBodies, steps ),           # Positions
        zeros( 3, numBodies, steps ),           # Velocities
        zeros( 3, numBodies, steps ),           # Accelerations
        zeros( numBodies, steps ),              # Potential Energies
        zeros( numBodies, steps ),              # Kinetic Energies
        zeros( steps ),                         # Total Potential Energy
        zeros( steps ),                         # Total Kinetic Energy
        zeros( steps ),                         # Temperature at each step
        zeros( steps ),                         # Total Energy
        zeros( steps ),                         # Virial
        0 );                                    # Completed steps (start at 0)

# Outer constructor for empty object
MolecularDynamicsTrial() = MolecularDynamicsTrial(1,1,0,0,0,0,1,0);

# Comparator for MolecularDynamicsTrial type
function ==( r::MolecularDynamicsTrial, p::MolecularDynamicsTrial )
    for field in names( MolecularDynamicsTrial )
        if r.( field ) != p.( field )
            return false
        end
    end
    return true
end

# Check if trial is finished
function isfinished( r::MolecularDynamicsTrial )
    return r.currentStep == r.steps
end

# Duplicate a trial
function duplicate( r::MolecularDynamicsTrial )
    p = MolecularDynamicsTrial()
    for field in names( MolecularDynamicsTrial )
        p.(field) = r.(field)
    end
    assert(p == r)
    return p
end

# Extract a sub-trial TODO : reimplement using slicedim(A,d,i)
function subtrial( r::MolecularDynamicsTrial, n::Int, m::Int )
    if r.steps < m
        error( "Out of bounds: final step must be less than r.steps" )
    elseif m < n
        error( "first step must be less than final step" )
    elseif n < 1
        error( "Out of bounds: first step must be greater than or equal to 1")
    elseif n > r.currentStep
        error( "first step must be less than current step" )
    else
        steps = m - n + 1
        numBodies = r.numBodies
        p = MolecularDynamicsTrial(numBodies,1,0,0,steps,0,1,0);
        for field in names( MolecularDynamicsTrial )
            if :steps == field
                continue                            # clobber not total steps
            elseif isa( r.(field), Array )
                fsize = size( p.(field) )           # dimensions of the field
                perstep = prod(fsize[1:(end-1)])    # elements per step
                first = (n-1)*perstep + 1           # first element to copy
                last = m*perstep                    # last element to copy
                assert((last-first+1)==perstep*steps)   # correct num elements
                i = 1
                for j in first:last                 # linear indexing to set p
                    p.(field)[i] = r.(field)[j]
                    i+=1
                end
            else
                p.(field) = r.(field)
            end
        end
        p.currentStep = r.currentStep + 1 - n
        return p
    end
end

# Lop off unfinished steps
function mdtruncate( r::MolecularDynamicsTrial )
    if isfinished( r )
        return r
    else
        return subtrial( r, 1, r.currentStep )      # go till last finished step
    end
end

# Add steps, if you need to run for longer
function addsteps( r::MolecularDynamicsTrial, n::Int )
    if n<0
        error( "Must specify non-negative number of steps to add" )
    else
        steps = r.steps + n
        numBodies = r.numBodies
        p = MolecularDynamicsTrial(numBodies,1,0,0,steps,0,1,0);
        for field in names( MolecularDynamicsTrial )
            if :steps == field
                continue                            # clobber not total steps
            elseif isa( r.(field), Array )
                i = 1
                for val in r.(field)[1:end]         # linear indexing to set p
                    p.(field)[i] = val
                    i+=1
                end
            else
                p.(field) = r.(field)               # copy scalar fields
            end
        end
        return p
    end
end


text:
struct ModalInstance <: AbstractVector{SubDataFrame}
    # worlds::Set{<:AbstractOntology} # TODO add set of worlds, for the semantics (DT)
    rows::Vector{DataFrameRow}
end

ModalInstance() = ModalInstance(DataFrameRow[])
DataFrames.nrow(::Array{DataFrameRow,1}) = 1

function show(io::IO, ::MIME"text/plain", mi::ModalInstance)
    println(io, "Modal instance with $(length(mi.rows)) frames")
    for i in 1:length(mi.rows)
        println(io, "Frame #$i")
        println(io, mi.rows[i])
    end
end

size(mi::ModalInstance) = (nrow(mi.rows),)
getindex(mi::ModalInstance, i::Int) = mi.rows[i]
length(mi::ModalInstance, fid::Int) = max([length(mi[fid][i]) for i in 1:length(mi[fid])]...)

mutable struct ModalFrame <: AbstractVector{SubDataFrame}
    # name::String # TODO add name?
    dimension::Int  # (0=static, 1=timeseries, 2=images, ecc.)
    data::DataFrame
    views::Vector{DataFrameRow} # intraframe instances (i.e., first, second, .., series)
end

dimension(mf::ModalFrame) = mf.dimension

@inline function ModalFrame(data::DataFrame)
    colnames = names(data)

    @assert length(colnames) > 0 && nrow(data) > 0 "There must be at least one column and one instance in data to compute its dimension."

    dimension = ndims(data[1, Symbol(colnames[1])])

    for colname in colnames
        curr_dimension = ndims(data[1, Symbol(colname)])
        @assert dimension == curr_dimension "Each attribute must have the same dimension."
    end

    # if dimension > 0
    #     typ = eltype(data[1, Symbol(colnams[1])])
    #     for colname in colnames
    #         curr_typ = eltype(data[1, Symbol(colname)])
    #         @assert typ == cur_typ "Each modal frame with higher-than-zero dimension must have the same element types."
    #     end
    # end

    views = Vector{DataFrameRow}(undef, nrow(data))
    for i in 1:nrow(data)
        views[i] = view(data, i, :)
    end
    ModalFrame(dimension, data, views)
end

size(mf::ModalFrame) = (nrow(mf.data),)
getindex(mf::ModalFrame, i::Int) = mf.views[i]

struct ClassificationDataset
    ldim::Int
    hdim::Int
    frames::Vector{ModalFrame}
    instances::Vector{ModalInstance}
    domains::Dict{Int, Dict{Symbol, Vector{Number}}}
    classes::CategoricalArray
    unique_classes::Vector
end

function ClassificationDataset(frames::Vector{ModalFrame}, classes::CategoricalArray)
    instances = ModalInstance[]
    domains = Dict{Int, Dict{Symbol, Vector{Number}}}()

    n_instances = size(frames[1])

    for i in 1:length(frames)
        curr_n_instances = size(frames[i])
        @assert n_instances == curr_n_instances "Each frame must have the same number of (intraframe) instances!"
    end

    @assert n_instances[1] == length(classes) "Each instance must have one and only one class!"

    # Computing the instances.
    # for each instance
    for i in 1:n_instances[1]
        instance = ModalInstance()
        # for each frame
        for f in 1:length(frames)
            # push the instance
            push!(instance.rows, frames[f][i])
        end
        push!(instances, instance)
    end

    # Computing the domains.
    # for each frame
    for f in 1:length(frames)
        attrs = names(frames[f][1])
        # for each attribute
        attr_domain = Dict{Symbol, Vector{Number}}()
        for attr in attrs
            # convert to symbol
            A = Symbol(attr)
            # for each instance
            # for i in 1:n_instances[1]

            # end
            attr_domain[A] = collect(Set([v for i in 1:n_instances[1] for v in frames[f][i][A]]))
            domains[f] = attr_domain
        end
    end

    ldim, hdim = extrema([dimension(i) for i in frames])

    ClassificationDataset(ldim, hdim, frames, instances, domains, classes, unique(classes))
end

instance(cs::ClassificationDataset, i::Int) = cs.instances[i] # TODO cs -> ds (either sup or unsupervised)
instance(cs::ClassificationDataset, i::Int, f::Int) = cs.instances[i][f]
function attributes(ds::ClassificationDataset)
    d = Dict{Int,Array{String,1}}()
    for (fid,frame) in enumerate(ds.frames)
        d[fid] = names(frame.data)
    end
    return d
end
attributes(ds::ClassificationDataset, fid::Int) = names(ds.frames[fid].data)

length(ds::ClassificationDataset) = length(ds.instances)

function show(io::IO, ::MIME"text/plain", ds::ClassificationDataset)
    println(io, "Classification dataset with $(length(ds.instances)) instances")

    for (i, inst) in enumerate(ds.instances)
        println(io, "Instance #$(i):")
        println(io, instance(ds, i))
        println(io)
    end
end

# TODO
# possible test: all([auslan.instances[i].rows[1][attr] === auslan.frames[1].data[i,attr] for i in 1:length(auslan.instances) for attr in attributes(auslan, 1)])
function transform!(ds::ClassificationDataset, f::Function, fid::Int; kwargs...)
    for attr in attributes(ds, fid)
        for i in 1:nrow(ds.frames[fid].data)
            ds.frames[fid].data[i,attr] = f(ds.frames[fid].data[i,attr]; kwargs...)
        end
    end
    return ds
end

# function Base.minimum(ds::ClassificationDataset)
#     d = Dict{Int,Array{Float64,1}}()

#     for (fid,frame) in enumerate(ds.frames)
#         frame = ds.frames[fid]
#         if frame.dimension == 0
#             d[fid] = missing
#         else
#             d[fid] = [Base.minimum(vcat(frame.data[!,attr]...)) for attr in attributes(ds, fid)]
#         end
#     end
#     return d
# end

# function Base.maximum(ds::ClassificationDataset)
#     d = Dict{Int,Array{Float64,1}}()

#     for (fid,frame) in enumerate(ds.frames)
#         frame = ds.frames[fid]
#         if frame.dimension == 0
#             d[fid] = missing
#         else
#             d[fid] = [Base.maximum(vcat(frame.data[!,attr]...)) for attr in attributes(ds, fid)]
#         end
#     end
#     return d
# end

# TODO fix problem with select_attributes!
# TODO make it with bang!
function select_attributes(ds::ClassificationDataset, frame_attributes::Dict{Int, Array{Int, 1}})
    frames = collect(keys(frame_attributes))

    # Frame level.
    for i in 1:length(frames)
        ds.frames[frames[i]].data = ds.frames[frames[i]].data[:, frame_attributes[frames[i]]]
    end

    modal_frames = ModalFrame[]
    for i in 1:length(frames)
        push!(modal_frames, ModalFrame(ds.frames[frames[i]].data))
    end

    return ClassificationDataset(modal_frames, ds.classes)
end

# function transform!(ds::ClassificationDataset, f::Function; kwargs...)

# using MultivariateTimeSeries
# X, y = read_data_labeled(joinpath(dirname(pathof(GBDTs)), "..", "data", "auslan_youtube8"));
#
# # adesso voglio passare da MTS a ClassificationDataset
#
# my_frame = DataFrame()
#
# for attr in names(X)
#     insertcols!(my_frame, attr => Array{Float64,1}[])
# end
#
# # for each instance
# for i in 1:length(X)
#     array_data = Array{Float64, 2}(X[i])
#     push!(my_frame, [array_data[:, j] for j in 1:length(names(X))])
# end
#
# auslan = ClassificationDataset([ModalFrame(my_frame)], CategoricalArray(y))

function sample(ds::ClassificationDataset,
        range::Union{UnitRange{Int},Nothing};
        invert::Bool=false)

    n_instances = length(ds.instances)

    @assert range ⊆ 1:n_instances "The range must be a subset of 1:$(n_instances)."

    if ismissing(range)
        return ds
    else

        output_classes = Int[]
        timeseries_df = DataFrame()

        attrs = attributes(ds, 1) # TODO up to now, only timeseries

        for attr in attrs
            insertcols!(timeseries_df, attr => Array{Float64, 1}[]) # TODO also here Array{Float64, 1}[], i.e., timeseries
        end

        if invert == false
            for idx in range
                push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                push!(output_classes, ds.classes[idx])
            end
        else
            for idx in setdiff(1:n_instances,range)
                push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                push!(output_classes, ds.classes[idx])
            end
        end
    end

    return ClassificationDataset([ModalFrame(timeseries_df)], CategoricalArray(output_classes))
end

function sample(ds::ClassificationDataset,
    percentage::Float64;
    seed::Int=1,
    bias::Float64=0.0, # bias to uniform
    replace::Bool=true,
    invert::Bool=false)

    @assert 0 ≤ percentage && percentage ≤ 1 "Percentage p must be 0 ≤ p ≤ 1."

    n_instances = length(ds.instances)
    n_instances_per_class = countmap(ds.classes)
    classes = keys(n_instances_per_class)
    n_classes = length(n_instances_per_class)

    instances_per_class = Dict{CategoricalArrays.CategoricalValue{String, UInt32}, Array{Int, 1}}()

    for class in classes
        instances_per_class[class] = ModalInstance[]
    end

    instances_to_sample = Dict{CategoricalArrays.CategoricalValue{String, UInt32}, Int64}()

    for class in classes
        sample_size = Int(round(percentage * ((1-bias)*n_instances_per_class[class] + bias*n_instances / n_classes)))
        instances_to_sample[class] = sample_size
    end

    for i in 1:n_instances
        push!(instances_per_class[ds.classes[i]], i)
    end

    r = Random.seed!(seed)

    output_classes = String[]
    timeseries_df = DataFrame()

    attrs = attributes(ds, 1) # TODO up to now, only timeseries

    for attr in attrs
        insertcols!(timeseries_df, attr => Array{Float64, 1}[]) # TODO also here Array{Float64, 1}[], i.e., timeseries
    end

    if replace
        for class in classes
            for i in 1:instances_to_sample[class]
                idx = instances_per_class[class][rand(r, 1:n_instances_per_class[class])]
                push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                push!(output_classes, ds.classes[idx])
                # @show idx
                # @show ds.classes[idx]
            end
        end
    else
        for class in classes
            n_eligible = length(instances_per_class[class])
            selected = instances_per_class[class]

            # the selected elements are those at the end of then selected array
            for i in 1:instances_to_sample[class]
                idx = rand(r, 1:n_eligible)
                selected[idx], selected[n_eligible] = selected[n_eligible], selected[idx]
                n_eligible -= 1
            end

            # take the last instances_to_sample[class] instances
            if invert == false
                for i in n_instances_per_class[class]-instances_to_sample[class]+1:n_instances_per_class[class]
                    idx = selected[i]
                    push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                    push!(output_classes, ds.classes[idx])
                end
            # take the first n_instances_per_class[class]-instances_to_sample[class] instances
            else
                for i in 1:n_instances_per_class[class]-instances_to_sample[class]
                    idx = selected[i]
                    push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                    push!(output_classes, ds.classes[idx])
                end
            end
        end
    end

    return ClassificationDataset([ModalFrame(timeseries_df)], CategoricalArray(output_classes))
end

function paa(x::AbstractArray{T} where T <: Real;
        n_chunks::Union{Missing,Int}=missing,
        f::Function=mean,
        decdigits::Int=4,
        kwargs...)
    N = length(x)
    if ismissing(n_chunks) || n_chunks ≥ N return x
    else
        @assert 0 ≤ n_chunks && n_chunks ≤ N "The number of chunks must be in [0,$(N)]"

        y = Array{Float64}(undef, n_chunks) # TODO Float64?
        for i in 1:n_chunks
            l = Int(ceil((N*(i-1)/n_chunks) + 1))
            h = Int(ceil(N*i/n_chunks))
            y[i] = round(f(x[l:h]; kwargs...), digits=decdigits)
        end
        return y
    end
end

######################################################
const _ARFF_SPACE       = UInt8(' ')
const _ARFF_COMMENT     = UInt8('%')
const _ARFF_AT          = UInt8('@')
const _ARFF_SEP         = UInt8(',')
const _ARFF_NEWLINE     = UInt8('\n')
const _ARFF_NOMSTART    = UInt8('{')
const _ARFF_NOMEND      = UInt8('}')
const _ARFF_ESC         = UInt8('\\')
const _ARFF_MISSING     = UInt8('?')
const _ARFF_RELMARK     = UInt8('\'')

function readARFF(p::String)
    open(p, "r") do io
        df = DataFrame()
        classes = String[]
        lines = readlines(io)
        for i in 1:length(lines)
            line = lines[i]
            # If not empty line or comment
            if !isempty(line)
                if UInt8(line[1]) != _ARFF_COMMENT
                    sline = split(line, " ")
                    # println(sline[1][1])
                    # If the first symbol is @
                    if UInt8(sline[1][1]) == _ARFF_AT
                        # If @relation
                        if sline[1][2:end] == "relation"
                            # println("Relation: " * sline[2])
                        end

                        # if sline[1][2:end] == "attribute" && sline[2] == "class"
                        #     classes = sline[3][2:end-1]
                        #     println(classes)
                        # end
                    # data, first char is '
                    elseif UInt8(sline[1][1]) == _ARFF_RELMARK
                        sline[1] = sline[1][2:end]
                        data_and_class = split(sline[1],"\'")
                        string_data = split(data_and_class[1], "\\n")
                        class = data_and_class[2][2:end]

                        if isempty(names(df))
                            for i in 1:length(string_data)
                                insertcols!(df, Symbol("A" * string(i)) => Array{Float64, 1}[]) # add the attributes as 1,2,3,ecc.
                            end
                        end

                        float_data = Dict{Int,Vector{Float64}}()

                        for i in 1:length(string_data)
                            float_data[i] = map(x->parse(Float64,x), split(string_data[i], ","))
                        end

                        # @show float_data

                        push!(df, [float_data[i] for i in 1:length(string_data)])
                        push!(classes, class)
                        # @show data
                        # @show class
                    end
                end
            end
            #error(1)
        end
        return ClassificationDataset([ModalFrame(df)], CategoricalArray(classes))
    end
end


text:
function get_cigar_funcs(dimensions)
    fs = Dict()

    function cigar()
        f(x) = x[1]^2 + 10^6 * sum(x[2:end].^2)
        return f
    end

    for d in dimensions
        fs[d] = cigar()
    end

    return fs
end

function get_linear_funcs(dimensions)
    fs = Dict()

    function linear()
        f(x) = sum(x)
        return f
    end

    for d in dimensions
        fs[d] = linear()
    end

    return fs
end

function get_box_funcs(dimensions, alpha=1/2)
    fs = Dict()

    function box()
        f(x) = sum(abs.(x)) <= alpha
        return f
    end

    for d in dimensions
        fs[d] = box()
    end

    return fs
end

function get_step_funcs(dimensions)
    fs = Dict()

    function step()
        f(x) = begin
            v = sum(x)
            if v < 0
                return 0
            elseif v < 2
                return 2
            elseif v < 4
                return 4
            else
                return 5
            end
        end

        return f
    end

    for d in dimensions
        fs[d] = step()
    end

    return fs
end

function get_sinc_funcs(dimensions)
    fs = Dict()

    function msinc()
        f(x) = *(sinc.(x)...)
        return f
    end

    for d in dimensions
        fs[d] = msinc()
    end

    return fs
end

function get_sphere_funcs(dimensions)
    fs = Dict()

    function sphere()
        f(x) = sum(x.^2)
        return f
    end

    for d in dimensions
        fs[d] = sphere()
    end

    return fs
end

function get_rastrigin_funcs(dimensions, A=10.0)
    fs = Dict()

    function rastrigrin(n)
        f(x) = A*n + sum(x.^2 - A * cos.(2π * x))
        return f
    end

    for d in dimensions
        fs[d] = rastrigrin(d)
    end

    return fs
end

function get_rosenbrock_funcs(dimensions)
    fs = Dict()

    function rosenbrock(n)
        f(x) = sum([100 * (x[i+1] - x[i]^2) + (1 - x[i])^2 for i in 1:(n-1)])
        return f
    end

    for d in dimensions
        fs[d] = rosenbrock(d)
    end

    return fs
end

function get_styblinskitang_funcs(dimensions)
    fs = Dict()

    function styblinskitang()
        f(x) = sum(x.^4 - 16x.^2 + 5x) / 2
        return f
    end

    for d in dimensions
        fs[d] = styblinskitang()
    end

    return fs
end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
@testset "Update" begin
    r1 = range(String, :r1; values=["a", "b", "c"])
    r2 = range(Int, :r2; lower=1, upper=3, scale=exp10)
    r3 = range(Float64, :r3; lower=1, upper=Inf, origin=4, unit=1)
    r4 = range(Float32, :r4; lower=-Inf, upper=Inf, origin=0, unit=1)

    @testset "Move Swarm" begin
        X = [0.14280010160187237 0.49409071076694583 0.3631091876311819  2.7429797605672808 3.9372495283243105 -0.8067647083847199;
             0.295010672512568   0.4534584876713112  0.25153083981612073 2.3976392099947    3.6569395920512977  0.420991611378423 ;
             0.12881930462550284 0.2617407494915029  0.6094399458829942  2.5742724788985445 3.6354556967115146  0.673601904658013 ]
        gbest_X = [-0.8571998983981276   0.49409071076694583  1.3631091876311818  1.7429797605672808  7.93724952832431   -0.8067647083847199;
                    1.295010672512568   -0.5465415123286887   0.25153083981612073 5.397639209994701   3.6569395920512977 -99.57900838862157 ;
                    0.12881930462550284  1.261740749491503   -0.3905600541170058  2.5742724788985445 -0.3645443032884854  100.67360190465801]
        dummy_state = PSO.ParticleSwarmState(
            (r1, r2, r3, r4),
            nothing,
            (1:3, 4, 5, 6),
            X,
            zeros(3, 6),
            X,
            gbest_X,
            Float64[],
            Float64[]
        )
        PSO.move!(StableRNG(8888), dummy_state, 1.0, 2.0, 2.0)

        @test all(0 .<= X[:, 1:3] .<= 1)
        @test all(sum(X[:, 1:3]; dims=2) .== 1)
        @test all(r2.lower .<= X[:, 4] .<= r2.upper)
        @test all(r3.lower .<= X[:, 5] .<= r3.upper)
        @test all(r4.lower .<= X[:, 6] .<= r4.upper)
        @test X ≈ [1.4861521012405697e-16 0.3306965950637089  0.6693034049362909    2.426443736254841  9.889381738850881  -0.8067647083847199;
                   0.6561834218566733     0.17876621097030176 0.16505036717302488   3.0                3.6569395920512977 -40.736143561117814;
                   0.11411862296972314    0.8858813770302767  1.967051803731305e-16 2.5742724788985445 1.0                 73.175744055348   ]
    end

    rng = StableRNG(1234)
    state = PSO.initialize(rng, [r1, r2, r3, r4], 3)
    PSO.retrieve!(rng, state)
    measurements = [0.5, 0.1, 0.3]
    PSO.pbest!(state, measurements, 0.25)
    PSO.gbest!(state)

    @testset "Update personal best" begin
        @test state.pbest == measurements
        @test state.pbest_X ≈ [0.35710007620140427 0.37056803307520936 0.2723318907233864  2.7429797605672808 3.9372495283243105 -0.8067647083847199;
                               0.471258004384426   0.3400938657534834  0.18864812986209056 2.3976392099947    3.6569395920512977  0.420991611378423 ;
                               0.09661447846912713 0.19630556211862715 0.7070799594122457  2.5742724788985445 3.6354556967115146  0.6736019046580138]
    end

    @testset "Update global best" begin
        @test all(state.gbest .== 0.1)
        @test all(state.gbest_X .== state.pbest_X[2, :]')
    end
end


text:
# Julia wrapper for header: sundials_band.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_config.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_dense.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_direct.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_fconfig.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_fnvector.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_futils.h
# Automatically generated using Clang.jl

function SUNDIALSFileOpen(filename, modes)
    ccall((:SUNDIALSFileOpen, libsundials_sundials), Ptr{FILE}, (Cstring, Cstring), filename, modes)
end

function SUNDIALSFileClose(fp)
    ccall((:SUNDIALSFileClose, libsundials_sundials), Cvoid, (Ptr{FILE},), fp)
end
# Julia wrapper for header: sundials_iterative.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_lapack.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_linearsolver.h
# Automatically generated using Clang.jl

function SUNLinSolNewEmpty()
    ccall((:SUNLinSolNewEmpty, libsundials_sundials), SUNLinearSolver, ())
end

function SUNLinSolFreeEmpty(S::SUNLinearSolver)
    ccall((:SUNLinSolFreeEmpty, libsundials_sundials), Cvoid, (SUNLinearSolver,), S)
end

function SUNLinSolGetType(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType, libsundials_sundials), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolGetID(S::SUNLinearSolver)
    ccall((:SUNLinSolGetID, libsundials_sundials), SUNLinearSolver_ID, (SUNLinearSolver,), S)
end

function SUNLinSolSetATimes(S, A_data, ATimes::ATimesFn)
    ccall((:SUNLinSolSetATimes, libsundials_sundials), Cint, (SUNLinearSolver, Ptr{Cvoid}, ATimesFn), S, A_data, ATimes)
end

function SUNLinSolSetPreconditioner(S, P_data, Pset::PSetupFn, Psol::PSolveFn)
    ccall((:SUNLinSolSetPreconditioner, libsundials_sundials), Cint, (SUNLinearSolver, Ptr{Cvoid}, PSetupFn, PSolveFn), S, P_data, Pset, Psol)
end

function SUNLinSolSetScalingVectors(S::SUNLinearSolver, s1::N_Vector, s2::N_Vector)
    ccall((:SUNLinSolSetScalingVectors, libsundials_sundials), Cint, (SUNLinearSolver, N_Vector, N_Vector), S, s1, s2)
end

function SUNLinSolSetScalingVectors(S, s1, s2)
    __s1 = convert(NVector, s1)
    __s2 = convert(NVector, s2)
    SUNLinSolSetScalingVectors(S, convert(N_Vector, __s1), convert(N_Vector, __s2))
end

function SUNLinSolInitialize(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize, libsundials_sundials), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetup(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup, libsundials_sundials), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function SUNLinSolSolve(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve, libsundials_sundials), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    SUNLinSolSolve(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolNumIters(S::SUNLinearSolver)
    ccall((:SUNLinSolNumIters, libsundials_sundials), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolResNorm(S::SUNLinearSolver)
    ccall((:SUNLinSolResNorm, libsundials_sundials), realtype, (SUNLinearSolver,), S)
end

function SUNLinSolResid(S::SUNLinearSolver)
    ccall((:SUNLinSolResid, libsundials_sundials), N_Vector, (SUNLinearSolver,), S)
end

function SUNLinSolLastFlag(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag, libsundials_sundials), sunindextype, (SUNLinearSolver,), S)
end

function SUNLinSolSpace(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace, libsundials_sundials), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree(S::SUNLinearSolver)
    ccall((:SUNLinSolFree, libsundials_sundials), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: sundials_math.h
# Automatically generated using Clang.jl

function SUNRpowerI(base::realtype, exponent::Cint)
    ccall((:SUNRpowerI, libsundials_sundials), realtype, (realtype, Cint), base, exponent)
end

function SUNRpowerI(base::realtype, exponent)
    SUNRpowerI(base, convert(Cint, exponent))
end

function SUNRpowerR(base::realtype, exponent::Cint)
    ccall((:SUNRpowerR, libsundials_sundials), realtype, (realtype, realtype), base, exponent)
end

function SUNRpowerR(base::realtype, exponent)
    SUNRpowerR(base, convert(Cint, exponent))
end
# Julia wrapper for header: sundials_matrix.h
# Automatically generated using Clang.jl

function SUNMatNewEmpty()
    ccall((:SUNMatNewEmpty, libsundials_sundials), SUNMatrix, ())
end

function SUNMatFreeEmpty(A::SUNMatrix)
    ccall((:SUNMatFreeEmpty, libsundials_sundials), Cvoid, (SUNMatrix,), A)
end

function SUNMatCopyOps(A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatCopyOps, libsundials_sundials), Cint, (SUNMatrix, SUNMatrix), A, B)
end

function SUNMatGetID(A::SUNMatrix)
    ccall((:SUNMatGetID, libsundials_sundials), SUNMatrix_ID, (SUNMatrix,), A)
end

function SUNMatClone(A::SUNMatrix)
    ccall((:SUNMatClone, libsundials_sundials), SUNMatrix, (SUNMatrix,), A)
end

function SUNMatDestroy(A::SUNMatrix)
    ccall((:SUNMatDestroy, libsundials_sundials), Cvoid, (SUNMatrix,), A)
end

function SUNMatZero(A::SUNMatrix)
    ccall((:SUNMatZero, libsundials_sundials), Cint, (SUNMatrix,), A)
end

function SUNMatCopy(A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatCopy, libsundials_sundials), Cint, (SUNMatrix, SUNMatrix), A, B)
end

function SUNMatScaleAdd(c::realtype, A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatScaleAdd, libsundials_sundials), Cint, (realtype, SUNMatrix, SUNMatrix), c, A, B)
end

function SUNMatScaleAddI(c::realtype, A::SUNMatrix)
    ccall((:SUNMatScaleAddI, libsundials_sundials), Cint, (realtype, SUNMatrix), c, A)
end

function SUNMatMatvecSetup(A::SUNMatrix)
    ccall((:SUNMatMatvecSetup, libsundials_sundials), Cint, (SUNMatrix,), A)
end

function SUNMatMatvec(A::SUNMatrix, x::N_Vector, y::N_Vector)
    ccall((:SUNMatMatvec, libsundials_sundials), Cint, (SUNMatrix, N_Vector, N_Vector), A, x, y)
end

function SUNMatMatvec(A, x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    SUNMatMatvec(A, convert(N_Vector, __x), convert(N_Vector, __y))
end

function SUNMatSpace(A::SUNMatrix, lenrw, leniw)
    ccall((:SUNMatSpace, libsundials_sundials), Cint, (SUNMatrix, Ptr{Clong}, Ptr{Clong}), A, lenrw, leniw)
end
# Julia wrapper for header: sundials_mpi_types.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_nonlinearsolver.h
# Automatically generated using Clang.jl

function SUNNonlinSolNewEmpty()
    ccall((:SUNNonlinSolNewEmpty, libsundials_sundials), SUNNonlinearSolver, ())
end

function SUNNonlinSolFreeEmpty(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolFreeEmpty, libsundials_sundials), Cvoid, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolGetType(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolGetType, libsundials_sundials), SUNNonlinearSolver_Type, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolInitialize(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolInitialize, libsundials_sundials), Cint, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolSetup(NLS::SUNNonlinearSolver, y::N_Vector, mem)
    ccall((:SUNNonlinSolSetup, libsundials_sundials), Cint, (SUNNonlinearSolver, N_Vector, Ptr{Cvoid}), NLS, y, mem)
end

function SUNNonlinSolSetup(NLS, y, mem)
    __y = convert(NVector, y)
    SUNNonlinSolSetup(NLS, convert(N_Vector, __y), mem)
end

function SUNNonlinSolSolve(NLS::SUNNonlinearSolver, y0::N_Vector, y::N_Vector, w::N_Vector, tol::realtype, callLSetup::Cint, mem)
    ccall((:SUNNonlinSolSolve, libsundials_sundials), Cint, (SUNNonlinearSolver, N_Vector, N_Vector, N_Vector, realtype, Cint, Ptr{Cvoid}), NLS, y0, y, w, tol, callLSetup, mem)
end

function SUNNonlinSolSolve(NLS, y0, y, w, tol, callLSetup, mem)
    __y0 = convert(NVector, y0)
    __y = convert(NVector, y)
    __w = convert(NVector, w)
    SUNNonlinSolSolve(NLS, convert(N_Vector, __y0), convert(N_Vector, __y), convert(N_Vector, __w), tol, convert(Cint, callLSetup), mem)
end

function SUNNonlinSolFree(NLS::SUNNonlinearSolver)
    ccall((:SUNNonlinSolFree, libsundials_sundials), Cint, (SUNNonlinearSolver,), NLS)
end

function SUNNonlinSolSetSysFn(NLS::SUNNonlinearSolver, SysFn::SUNNonlinSolSysFn)
    ccall((:SUNNonlinSolSetSysFn, libsundials_sundials), Cint, (SUNNonlinearSolver, SUNNonlinSolSysFn), NLS, SysFn)
end

function SUNNonlinSolSetLSetupFn(NLS::SUNNonlinearSolver, SetupFn::SUNNonlinSolLSetupFn)
    ccall((:SUNNonlinSolSetLSetupFn, libsundials_sundials), Cint, (SUNNonlinearSolver, SUNNonlinSolLSetupFn), NLS, SetupFn)
end

function SUNNonlinSolSetLSolveFn(NLS::SUNNonlinearSolver, SolveFn::SUNNonlinSolLSolveFn)
    ccall((:SUNNonlinSolSetLSolveFn, libsundials_sundials), Cint, (SUNNonlinearSolver, SUNNonlinSolLSolveFn), NLS, SolveFn)
end

function SUNNonlinSolSetConvTestFn(NLS::SUNNonlinearSolver, CTestFn::SUNNonlinSolConvTestFn, ctest_data)
    ccall((:SUNNonlinSolSetConvTestFn, libsundials_sundials), Cint, (SUNNonlinearSolver, SUNNonlinSolConvTestFn, Ptr{Cvoid}), NLS, CTestFn, ctest_data)
end

function SUNNonlinSolSetMaxIters(NLS::SUNNonlinearSolver, maxiters::Cint)
    ccall((:SUNNonlinSolSetMaxIters, libsundials_sundials), Cint, (SUNNonlinearSolver, Cint), NLS, maxiters)
end

function SUNNonlinSolSetMaxIters(NLS, maxiters)
    SUNNonlinSolSetMaxIters(NLS, convert(Cint, maxiters))
end

function SUNNonlinSolGetNumIters(NLS::SUNNonlinearSolver, niters)
    ccall((:SUNNonlinSolGetNumIters, libsundials_sundials), Cint, (SUNNonlinearSolver, Ptr{Clong}), NLS, niters)
end

function SUNNonlinSolGetCurIter(NLS::SUNNonlinearSolver, iter)
    ccall((:SUNNonlinSolGetCurIter, libsundials_sundials), Cint, (SUNNonlinearSolver, Ptr{Cint}), NLS, iter)
end

function SUNNonlinSolGetNumConvFails(NLS::SUNNonlinearSolver, nconvfails)
    ccall((:SUNNonlinSolGetNumConvFails, libsundials_sundials), Cint, (SUNNonlinearSolver, Ptr{Clong}), NLS, nconvfails)
end
# Julia wrapper for header: sundials_nvector.h
# Automatically generated using Clang.jl

function N_VNewEmpty()
    ccall((:N_VNewEmpty, libsundials_sundials), N_Vector, ())
end

function N_VFreeEmpty(v::N_Vector)
    ccall((:N_VFreeEmpty, libsundials_sundials), Cvoid, (N_Vector,), v)
end

function N_VFreeEmpty(v)
    __v = convert(NVector, v)
    N_VFreeEmpty(convert(N_Vector, __v))
end

function N_VCopyOps(w::N_Vector, v::N_Vector)
    ccall((:N_VCopyOps, libsundials_sundials), Cint, (N_Vector, N_Vector), w, v)
end

function N_VCopyOps(w, v)
    __w = convert(NVector, w)
    __v = convert(NVector, v)
    N_VCopyOps(convert(N_Vector, __w), convert(N_Vector, __v))
end

function N_VGetVectorID(w::N_Vector)
    ccall((:N_VGetVectorID, libsundials_sundials), N_Vector_ID, (N_Vector,), w)
end

function N_VGetVectorID(w)
    __w = convert(NVector, w)
    N_VGetVectorID(convert(N_Vector, __w))
end

function N_VClone(w::N_Vector)
    ccall((:N_VClone, libsundials_sundials), N_Vector, (N_Vector,), w)
end

function N_VClone(w)
    __w = convert(NVector, w)
    N_VClone(convert(N_Vector, __w))
end

function N_VCloneEmpty(w::N_Vector)
    ccall((:N_VCloneEmpty, libsundials_sundials), N_Vector, (N_Vector,), w)
end

function N_VCloneEmpty(w)
    __w = convert(NVector, w)
    N_VCloneEmpty(convert(N_Vector, __w))
end

function N_VDestroy(v::N_Vector)
    ccall((:N_VDestroy, libsundials_sundials), Cvoid, (N_Vector,), v)
end

function N_VDestroy(v)
    __v = convert(NVector, v)
    N_VDestroy(convert(N_Vector, __v))
end

function N_VSpace(v::N_Vector, lrw, liw)
    ccall((:N_VSpace, libsundials_sundials), Cvoid, (N_Vector, Ptr{sunindextype}, Ptr{sunindextype}), v, lrw, liw)
end

function N_VSpace(v, lrw, liw)
    __v = convert(NVector, v)
    N_VSpace(convert(N_Vector, __v), lrw, liw)
end

function N_VGetArrayPointer(v::N_Vector)
    ccall((:N_VGetArrayPointer, libsundials_sundials), Ptr{realtype}, (N_Vector,), v)
end

function N_VGetArrayPointer(v)
    __v = convert(NVector, v)
    N_VGetArrayPointer(convert(N_Vector, __v))
end

function N_VSetArrayPointer(v_data, v::N_Vector)
    ccall((:N_VSetArrayPointer, libsundials_sundials), Cvoid, (Ptr{realtype}, N_Vector), v_data, v)
end

function N_VSetArrayPointer(v_data, v)
    __v = convert(NVector, v)
    N_VSetArrayPointer(v_data, convert(N_Vector, __v))
end

function N_VGetCommunicator(v::N_Vector)
    ccall((:N_VGetCommunicator, libsundials_sundials), Ptr{Cvoid}, (N_Vector,), v)
end

function N_VGetCommunicator(v)
    __v = convert(NVector, v)
    N_VGetCommunicator(convert(N_Vector, __v))
end

function N_VGetLength(v::N_Vector)
    ccall((:N_VGetLength, libsundials_sundials), sunindextype, (N_Vector,), v)
end

function N_VGetLength(v)
    __v = convert(NVector, v)
    N_VGetLength(convert(N_Vector, __v))
end

function N_VLinearSum(a::realtype, x::N_Vector, b::realtype, y::N_Vector, z::N_Vector)
    ccall((:N_VLinearSum, libsundials_sundials), Cvoid, (realtype, N_Vector, realtype, N_Vector, N_Vector), a, x, b, y, z)
end

function N_VLinearSum(a, x, b, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    N_VLinearSum(a, convert(N_Vector, __x), b, convert(N_Vector, __y), convert(N_Vector, __z))
end

function N_VConst(c::realtype, z::N_Vector)
    ccall((:N_VConst, libsundials_sundials), Cvoid, (realtype, N_Vector), c, z)
end

function N_VConst(c, z)
    __z = convert(NVector, z)
    N_VConst(c, convert(N_Vector, __z))
end

function N_VProd(x::N_Vector, y::N_Vector, z::N_Vector)
    ccall((:N_VProd, libsundials_sundials), Cvoid, (N_Vector, N_Vector, N_Vector), x, y, z)
end

function N_VProd(x, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    N_VProd(convert(N_Vector, __x), convert(N_Vector, __y), convert(N_Vector, __z))
end

function N_VDiv(x::N_Vector, y::N_Vector, z::N_Vector)
    ccall((:N_VDiv, libsundials_sundials), Cvoid, (N_Vector, N_Vector, N_Vector), x, y, z)
end

function N_VDiv(x, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    N_VDiv(convert(N_Vector, __x), convert(N_Vector, __y), convert(N_Vector, __z))
end

function N_VScale(c::realtype, x::N_Vector, z::N_Vector)
    ccall((:N_VScale, libsundials_sundials), Cvoid, (realtype, N_Vector, N_Vector), c, x, z)
end

function N_VScale(c, x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    N_VScale(c, convert(N_Vector, __x), convert(N_Vector, __z))
end

function N_VAbs(x::N_Vector, z::N_Vector)
    ccall((:N_VAbs, libsundials_sundials), Cvoid, (N_Vector, N_Vector), x, z)
end

function N_VAbs(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    N_VAbs(convert(N_Vector, __x), convert(N_Vector, __z))
end

function N_VInv(x::N_Vector, z::N_Vector)
    ccall((:N_VInv, libsundials_sundials), Cvoid, (N_Vector, N_Vector), x, z)
end

function N_VInv(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    N_VInv(convert(N_Vector, __x), convert(N_Vector, __z))
end

function N_VAddConst(x::N_Vector, b::realtype, z::N_Vector)
    ccall((:N_VAddConst, libsundials_sundials), Cvoid, (N_Vector, realtype, N_Vector), x, b, z)
end

function N_VAddConst(x, b, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    N_VAddConst(convert(N_Vector, x), b, convert(N_Vector, z))
end

function N_VDotProd(x::N_Vector, y::N_Vector)
    ccall((:N_VDotProd, libsundials_sundials), realtype, (N_Vector, N_Vector), x, y)
end

function N_VDotProd(x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    N_VDotProd(convert(N_Vector, __x), convert(N_Vector, __y))
end

function N_VMaxNorm(x::N_Vector)
    ccall((:N_VMaxNorm, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VMaxNorm(x)
    __x = convert(NVector, x)
    N_VMaxNorm(convert(N_Vector, __x))
end

function N_VWrmsNorm(x::N_Vector, w::N_Vector)
    ccall((:N_VWrmsNorm, libsundials_sundials), realtype, (N_Vector, N_Vector), x, w)
end

function N_VWrmsNorm(x, w)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    N_VWrmsNorm(convert(N_Vector, __x), convert(N_Vector, __w))
end

function N_VWrmsNormMask(x::N_Vector, w::N_Vector, id::N_Vector)
    ccall((:N_VWrmsNormMask, libsundials_sundials), realtype, (N_Vector, N_Vector, N_Vector), x, w, id)
end

function N_VWrmsNormMask(x, w, id)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __id = convert(NVector, id)
    N_VWrmsNormMask(convert(N_Vector, __x), convert(N_Vector, __w), convert(N_Vector, __id))
end

function N_VMin(x::N_Vector)
    ccall((:N_VMin, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VMin(x)
    __x = convert(NVector, x)
    N_VMin(convert(N_Vector, __x))
end

function N_VWL2Norm(x::N_Vector, w::N_Vector)
    ccall((:N_VWL2Norm, libsundials_sundials), realtype, (N_Vector, N_Vector), x, w)
end

function N_VWL2Norm(x, w)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    N_VWL2Norm(convert(N_Vector, __x), convert(N_Vector, __w))
end

function N_VL1Norm(x::N_Vector)
    ccall((:N_VL1Norm, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VL1Norm(x)
    __x = convert(NVector, x)
    N_VL1Norm(convert(N_Vector, __x))
end

function N_VCompare(c::realtype, x::N_Vector, z::N_Vector)
    ccall((:N_VCompare, libsundials_sundials), Cvoid, (realtype, N_Vector, N_Vector), c, x, z)
end

function N_VCompare(c, x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    N_VCompare(c, convert(N_Vector, __x), convert(N_Vector, __z))
end

function N_VInvTest(x::N_Vector, z::N_Vector)
    ccall((:N_VInvTest, libsundials_sundials), Cint, (N_Vector, N_Vector), x, z)
end

function N_VInvTest(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    N_VInvTest(convert(N_Vector, __x), convert(N_Vector, __z))
end

function N_VConstrMask(c::N_Vector, x::N_Vector, m::N_Vector)
    ccall((:N_VConstrMask, libsundials_sundials), Cint, (N_Vector, N_Vector, N_Vector), c, x, m)
end

function N_VConstrMask(c, x, m)
    __c = convert(NVector, c)
    __x = convert(NVector, x)
    __m = convert(NVector, m)
    N_VConstrMask(convert(N_Vector, __c), convert(N_Vector, __x), convert(N_Vector, __m))
end

function N_VMinQuotient(num::N_Vector, denom::N_Vector)
    ccall((:N_VMinQuotient, libsundials_sundials), realtype, (N_Vector, N_Vector), num, denom)
end

function N_VMinQuotient(num, denom)
    __num = convert(NVector, num)
    __denom = convert(NVector, denom)
    N_VMinQuotient(convert(N_Vector, __num), convert(N_Vector, __denom))
end

function N_VLinearCombination(nvec::Cint, c, X, z)
    ccall((:N_VLinearCombination, libsundials_sundials), Cint, (Cint, Ptr{realtype}, Ptr{N_Vector}, N_Vector), nvec, c, X, z)
end

function N_VLinearCombination(nvec, c, X, z)
    N_VLinearCombination(convert(Cint, nvec), c, X, z)
end

function N_VScaleAddMulti(nvec::Cint, a, x::N_Vector, Y, Z)
    ccall((:N_VScaleAddMulti, libsundials_sundials), Cint, (Cint, Ptr{realtype}, N_Vector, Ptr{N_Vector}, Ptr{N_Vector}), nvec, a, x, Y, Z)
end

function N_VScaleAddMulti(nvec, a, x, Y, Z)
    __x = convert(NVector, x)
    N_VScaleAddMulti(convert(Cint, nvec), a, convert(N_Vector, __x), Y, Z)
end

function N_VDotProdMulti(nvec::Cint, x::N_Vector, Y, dotprods)
    ccall((:N_VDotProdMulti, libsundials_sundials), Cint, (Cint, N_Vector, Ptr{N_Vector}, Ptr{realtype}), nvec, x, Y, dotprods)
end

function N_VDotProdMulti(nvec, x, Y, dotprods)
    __x = convert(NVector, x)
    N_VDotProdMulti(convert(Cint, nvec), convert(N_Vector, __x), Y, dotprods)
end

function N_VLinearSumVectorArray(nvec::Cint, a::realtype, X, b::realtype, Y, Z)
    ccall((:N_VLinearSumVectorArray, libsundials_sundials), Cint, (Cint, realtype, Ptr{N_Vector}, realtype, Ptr{N_Vector}, Ptr{N_Vector}), nvec, a, X, b, Y, Z)
end

function N_VLinearSumVectorArray(nvec, a, X, b, Y, Z)
    N_VLinearSumVectorArray(convert(Cint, nvec), a, X, b, Y, Z)
end

function N_VScaleVectorArray(nvec::Cint, c, X, Z)
    ccall((:N_VScaleVectorArray, libsundials_sundials), Cint, (Cint, Ptr{realtype}, Ptr{N_Vector}, Ptr{N_Vector}), nvec, c, X, Z)
end

function N_VScaleVectorArray(nvec, c, X, Z)
    N_VScaleVectorArray(convert(Cint, nvec), c, X, Z)
end

function N_VConstVectorArray(nvec::Cint, c::realtype, Z)
    ccall((:N_VConstVectorArray, libsundials_sundials), Cint, (Cint, realtype, Ptr{N_Vector}), nvec, c, Z)
end

function N_VConstVectorArray(nvec, c, Z)
    N_VConstVectorArray(convert(Cint, nvec), c, Z)
end

function N_VWrmsNormVectorArray(nvec::Cint, X, W, nrm)
    ccall((:N_VWrmsNormVectorArray, libsundials_sundials), Cint, (Cint, Ptr{N_Vector}, Ptr{N_Vector}, Ptr{realtype}), nvec, X, W, nrm)
end

function N_VWrmsNormVectorArray(nvec, X, W, nrm)
    N_VWrmsNormVectorArray(convert(Cint, nvec), X, W, nrm)
end

function N_VWrmsNormMaskVectorArray(nvec::Cint, X, W, id::N_Vector, nrm)
    ccall((:N_VWrmsNormMaskVectorArray, libsundials_sundials), Cint, (Cint, Ptr{N_Vector}, Ptr{N_Vector}, N_Vector, Ptr{realtype}), nvec, X, W, id, nrm)
end

function N_VWrmsNormMaskVectorArray(nvec, X, W, id, nrm)
    __id = convert(NVector, id)
    N_VWrmsNormMaskVectorArray(convert(Cint, nvec), X, W, convert(N_Vector, __id), nrm)
end

function N_VScaleAddMultiVectorArray(nvec::Cint, nsum::Cint, a, X, Y, Z)
    ccall((:N_VScaleAddMultiVectorArray, libsundials_sundials), Cint, (Cint, Cint, Ptr{realtype}, Ptr{N_Vector}, Ptr{Ptr{N_Vector}}, Ptr{Ptr{N_Vector}}), nvec, nsum, a, X, Y, Z)
end

function N_VScaleAddMultiVectorArray(nvec, nsum, a, X, Y, Z)
    N_VScaleAddMultiVectorArray(convert(Cint, nvec), convert(Cint, nsum), a, X, Y, Z)
end

function N_VLinearCombinationVectorArray(nvec::Cint, nsum::Cint, c, X, Z)
    ccall((:N_VLinearCombinationVectorArray, libsundials_sundials), Cint, (Cint, Cint, Ptr{realtype}, Ptr{Ptr{N_Vector}}, Ptr{N_Vector}), nvec, nsum, c, X, Z)
end

function N_VLinearCombinationVectorArray(nvec, nsum, c, X, Z)
    N_VLinearCombinationVectorArray(convert(Cint, nvec), convert(Cint, nsum), c, X, Z)
end

function N_VDotProdLocal(x::N_Vector, y::N_Vector)
    ccall((:N_VDotProdLocal, libsundials_sundials), realtype, (N_Vector, N_Vector), x, y)
end

function N_VDotProdLocal(x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    N_VDotProdLocal(convert(N_Vector, __x), convert(N_Vector, __y))
end

function N_VMaxNormLocal(x::N_Vector)
    ccall((:N_VMaxNormLocal, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VMaxNormLocal(x)
    __x = convert(NVector, x)
    N_VMaxNormLocal(convert(N_Vector, __x))
end

function N_VMinLocal(x::N_Vector)
    ccall((:N_VMinLocal, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VMinLocal(x)
    __x = convert(NVector, x)
    N_VMinLocal(convert(N_Vector, __x))
end

function N_VL1NormLocal(x::N_Vector)
    ccall((:N_VL1NormLocal, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VL1NormLocal(x)
    __x = convert(NVector, x)
    N_VL1NormLocal(convert(N_Vector, __x))
end

function N_VWSqrSumLocal(x::N_Vector, w::N_Vector)
    ccall((:N_VWSqrSumLocal, libsundials_sundials), realtype, (N_Vector, N_Vector), x, w)
end

function N_VWSqrSumLocal(x, w)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    N_VWSqrSumLocal(convert(N_Vector, __x), convert(N_Vector, __w))
end

function N_VWSqrSumMaskLocal(x::N_Vector, w::N_Vector, id::N_Vector)
    ccall((:N_VWSqrSumMaskLocal, libsundials_sundials), realtype, (N_Vector, N_Vector, N_Vector), x, w, id)
end

function N_VWSqrSumMaskLocal(x, w, id)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __id = convert(NVector, id)
    N_VWSqrSumMaskLocal(convert(N_Vector, __x), convert(N_Vector, __w), convert(N_Vector, __id))
end

function N_VInvTestLocal(x::N_Vector, z::N_Vector)
    ccall((:N_VInvTestLocal, libsundials_sundials), Cint, (N_Vector, N_Vector), x, z)
end

function N_VInvTestLocal(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    N_VInvTestLocal(convert(N_Vector, __x), convert(N_Vector, __z))
end

function N_VConstrMaskLocal(c::N_Vector, x::N_Vector, m::N_Vector)
    ccall((:N_VConstrMaskLocal, libsundials_sundials), Cint, (N_Vector, N_Vector, N_Vector), c, x, m)
end

function N_VConstrMaskLocal(c, x, m)
    __c = convert(NVector, c)
    __x = convert(NVector, x)
    __m = convert(NVector, m)
    N_VConstrMaskLocal(convert(N_Vector, __c), convert(N_Vector, __x), convert(N_Vector, __m))
end

function N_VMinQuotientLocal(num::N_Vector, denom::N_Vector)
    ccall((:N_VMinQuotientLocal, libsundials_sundials), realtype, (N_Vector, N_Vector), num, denom)
end

function N_VMinQuotientLocal(num, denom)
    __num = convert(NVector, num)
    __denom = convert(NVector, denom)
    N_VMinQuotientLocal(convert(N_Vector, __num), convert(N_Vector, __denom))
end

function N_VNewVectorArray(count::Cint)
    ccall((:N_VNewVectorArray, libsundials_sundials), Ptr{N_Vector}, (Cint,), count)
end

function N_VNewVectorArray(count)
    N_VNewVectorArray(convert(Cint, count))
end

function N_VCloneEmptyVectorArray(count::Cint, w::N_Vector)
    ccall((:N_VCloneEmptyVectorArray, libsundials_sundials), Ptr{N_Vector}, (Cint, N_Vector), count, w)
end

function N_VCloneEmptyVectorArray(count, w)
    __w = convert(NVector, w)
    N_VCloneEmptyVectorArray(convert(Cint, count), convert(N_Vector, __w))
end

function N_VCloneVectorArray(count::Cint, w::N_Vector)
    ccall((:N_VCloneVectorArray, libsundials_sundials), Ptr{N_Vector}, (Cint, N_Vector), count, w)
end

function N_VCloneVectorArray(count, w)
    __w = convert(NVector, w)
    N_VCloneVectorArray(convert(Cint, count), convert(N_Vector, __w))
end

function N_VDestroyVectorArray(vs, count::Cint)
    ccall((:N_VDestroyVectorArray, libsundials_sundials), Cvoid, (Ptr{N_Vector}, Cint), vs, count)
end

function N_VDestroyVectorArray(vs, count)
    N_VDestroyVectorArray(vs, convert(Cint, count))
end

function N_VGetVecAtIndexVectorArray(vs, index::Cint)
    ccall((:N_VGetVecAtIndexVectorArray, libsundials_sundials), N_Vector, (Ptr{N_Vector}, Cint), vs, index)
end

function N_VGetVecAtIndexVectorArray(vs, index)
    N_VGetVecAtIndexVectorArray(vs, convert(Cint, index))
end

function N_VSetVecAtIndexVectorArray(vs, index::Cint, w::N_Vector)
    ccall((:N_VSetVecAtIndexVectorArray, libsundials_sundials), Cvoid, (Ptr{N_Vector}, Cint, N_Vector), vs, index, w)
end

function N_VSetVecAtIndexVectorArray(vs, index, w)
    __w = convert(NVector, w)
    N_VSetVecAtIndexVectorArray(vs, convert(Cint, index, convert(N_Vector, __w)))
end
# Julia wrapper for header: sundials_types.h
# Automatically generated using Clang.jl

# Julia wrapper for header: sundials_version.h
# Automatically generated using Clang.jl

function SUNDIALSGetVersion(version, len::Cint)
    ccall((:SUNDIALSGetVersion, libsundials_sundials), Cint, (Cstring, Cint), version, len)
end

function SUNDIALSGetVersion(version, len)
    SUNDIALSGetVersion(version, convert(Cint, len))
end

function SUNDIALSGetVersionNumber(major, minor, patch, label, len::Cint)
    ccall((:SUNDIALSGetVersionNumber, libsundials_sundials), Cint, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cstring, Cint), major, minor, patch, label, len)
end

function SUNDIALSGetVersionNumber(major, minor, patch, label, len)
    SUNDIALSGetVersionNumber(major, minor, patch, label, convert(Cint, len))
end


text:
using LinearAlgebra
using Random
using Distributions
using Roots
using Cubature

include("../computation/post_prob_comp.jl")
include("../root_finding/my_bisection.jl")
include("../root_finding/zero_finding.jl")

# single point prediction
# given x_pred, compute the pdf of Z at each index
# also compute median and 95% CI and each index
function distribution_comp(x_pred, trainBasicInfo, sideSetInfo, param_gridInfo, alpha, Gamma)
    m = size(x_pred, 1)
    @assert m == 1
    dim = size(x_pred, 2)

    # compute the pdf
    pdf = z0 -> prob = pdf_z0(x_pred, z0, trainBasicInfo, sideSetInfo, param_gridInfo, alpha, Gamma)
    # find the upper bound for the pdf integral interval
    b0 = 30.0 # initial guess of the interval [0, b0] for the pdf integral
    kmax = 20
    tol = 1e-9
    b, iter = zero_finding(pdf, kmax, tol, b0)

    # compute mean, median, standard deviation, etc.
    # mean_z, median_z, stdv_z = median_comp(pdf, b)
    mean_z, stdv_z = median_comp(pdf, b)

    # compute confidence interval
    # p = 0.95 # by default we compute the 95% confidence interval
    # q = [0.8, 1.2]
    # eps = 1e-6
    # z2_initl = median_z + stdv_z
    # l_initl = pdf(z2_initl)
    # z1_initl = my_bisection(x -> pdf(x) - l_initl, 1e-3, median_z)
    # # z1_initl = find_zero(x -> pdf_temp(x) - l_initl, median_z - stdv_z)
    # initl = [l_initl, z1_initl, z2_initl]
    # CI = CI_comp_BTG(initl, p, q, pdf, eps, zu_new, median_z)

    # currently, median, mode and CI are not included
    dist_z = (mean = mean_z, stdv =  stdv_z, pdf = pdf, upperbd = b)
    return dist_z
end

function median_comp(pdf_z, b)
    # compute the mean
    mean_z = hquadrature(x -> x * pdf_z(x), 1e-5, b; reltol=1e-8, abstol=0, maxevals=0)[1]

    # compute the median
    # medianfun = x ->  hquadrature(pdf_z, 1e-5, x; reltol=1e-8, abstol=0, maxevals=0)[1] - 0.5
    # median_z = find_zero(medianfun, mean_z)

    # standard deviation
    variancefun = x -> (x - mean_z)^2 * pdf_z(x)
    variance_z = hquadrature(variancefun, 1e-5, b; reltol=1e-8, abstol=0, maxevals=0)[1]
    stdv_z = sqrt(variance_z)

    # return mean_z, median_z, stdv_z
    return mean_z, stdv_z

end

# compute confidence interval
function CI_comp_BTG(initl, p, q, pdf_z, eps, zu_new, median_z)
    l = initl[1]
    y1 = initl[2]
    y2 = initl[3]
    int = hquadrature(x -> pdf_z(x), y1, y2; reltol=1e-8, abstol=0, maxevals=0)
    Int = int[1]
    println("initial Int = $Int")
    flush(stdout)
    if int[1] < p
       while int[1] < p
            l_new = q[1] * l
            # println("1, l_new = $l_new")
            # flush(stdout1)
            # y1_new = find_zero(x -> pdf_z(x) - l_new, y1)
            # y2_new = find_zero(x -> pdf_z(x) - l_new, y2)
            y1_new = my_bisection(x -> pdf_z(x) - l_new, 1e-3, y1)
            y2_new = my_bisection(x -> pdf_z(x) - l_new, y2, zu_new)
            int = hquadrature(x -> pdf_z(x), y1_new, y2_new;
                                reltol=1e-8, abstol=0, maxevals=0)
            l = l_new
            y1 = y1_new
            y2 = y2_new
            CurentInt = int[1]
            println("Current Int = $CurentInt")
            flush(stdout)
        end
        l_low = l
        y1_low = y1
        y2_low = y2
        l_up = initl[1]
        y1_up = initl[2]
        y2_up = initl[3]
    elseif int[1] > p
        while  int[1] > p
            l_new = q[2] * l
            # println("2, l_new = $l_new")
            # flush(stdout1)
            # y1_new = find_zero(x -> pdf_z(x) - l_new, y1)
            # y2_new = find_zero(x -> pdf_z(x) - l_new, y2)
            y1_new = my_bisection(x -> pdf_z(x) - l_new, y1, median_z)
            y2_new = my_bisection(x -> pdf_z(x) - l_new, median_z, y2)
            int = hquadrature(x -> pdf_z(x), y1_new, y2_new;
                                reltol=1e-8, abstol=0, maxevals=0)
            l = l_new
            y1 = y1_new
            y2 = y2_new
        end
        l_low = initl[1]
        y1_low = initl[2]
        y2_low = initl[3]
        l_up = l
        y1_up = y1
        y2_up = y2
    end
    println("l_low = $l_low, y1_low = $y1_low, y2_low = $y2_low")
    flush(stdout)
    println("l_l_uplow = $l_up, y1_up = $y1_up, y2_up = $y2_up")
    flush(stdout)

    #  initialize middle points
    l_mid = 0
    y1_mid = 0
    y2_mid = 0
    N = 0
    # now with lower and upper bound, do sth like bisection...
    while (abs(int[1] - p) > eps) && N < 100
        l_mid = (l_up + l_low)/2
        # y1_mid = find_zero(x -> pdf_z(x) - l_mid, (y1_up + y1_low)/2)
        # y2_mid = find_zero(x -> pdf_z(x) - l_mid, (y2_up + y2_low)/2)
        y1_mid = my_bisection(x -> pdf_z(x) - l_mid, y1_low, y1_up)
        y2_mid = my_bisection(x -> pdf_z(x) - l_mid, y2_up, y2_low)
        int = hquadrature(x -> pdf_z(x), y1_mid, y2_mid; reltol=1e-8, abstol=0, maxevals=0)
        if int[1] > p
            l_low = l_mid
            y1_low = y1_mid
            y2_low = y2_mid
        else
            l_up = l_mid
            y1_up = y1_mid
            y2_up = y2_mid
        end
        N += 1
    end
    CI = [l_mid, y1_mid, y2_mid]
    return CI
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@testset "modules" begin
    @testset "find module definition" begin
        using Atom: moduledefinition

        let (path, line) = moduledefinition(Atom)
            @test path == atomjlfile
            @test line == 4
        end
        let (path, line) = moduledefinition(Junk)
            @test path == junkpath
            @test line == 1
        end
        let (path, line) = moduledefinition(Junk.SubJunk)
            @test path == subjunkspath
            @test line == 4
        end
    end

    @testset "find module files" begin
        using Atom: modulefiles, use_compiled_modules

        ## Revise-like module file detection
        # NOTE: only test when using compiled modules
        if use_compiled_modules()
            # # XXX modules.jl is unmaintained now and kinda broken for v1.6
            # # works for precompiled packages
            # let (parentfile, included_files) = modulefiles(Atom)
            #     for f in included_files
            #         @test f in atommodfiles
            #     end
            #     # NOTE: `+ 1` stands for Atom.jl (`parentfile`)
            #     @test length(included_files) + 1 === length(atommodfiles)
            #
            #     # TODO:
            #     # Revise-like module traverse fails to detect lazily loaded files even if
            #     # they are actually loaded.
            #     @test_broken webiofile in included_files
            #     @test_broken traceurfile in included_files
            # end

            # fails for non-precompiled packages
            @test_broken junkpath == modulefiles(Junk)[1]
        else
            @warn "skipped Revise-like module file detection test (precompiled modules are currently not used)"
        end

        ## CSTPraser-based module file detection
        # basic
        let included_files = normpath.(modulefiles("Atom", atomjlfile))
            for f in included_files
                @test f in atommodfiles
            end
            @test length(included_files) === length(atommodfiles)

            # TODO:
            # CSTParser-based module traverse is currently not able to find non-toplevel
            # `include` calls
            @test_broken webiofile in included_files
            @test_broken traceurfile in included_files
        end

        # safety checks for recursive `include`s
        @test (modulefiles("Main", joinpath′(fixturedir, "self_recur.jl"); inmod = true); true)
        @test (modulefiles("Main", joinpath′(fixturedir, "mutual_recur.jl"); inmod = true); true)
    end
end


text:
@set basis cc-pvdz
@test (@get basis) == "cc-pvdz"
@reset
@test (@get basis) == "sto-3g"


text:
function process_params()

    mi = zeros(Float64, nspec)
    mi_g = zeros(Float64, nspec)
    qi = zeros(Float64, nspec)

    dr = L / nz
    if geom=="spherical"
        dr = -dr #i.e., dr is negative
    end

    for ispec = 1:nspec
        mi_g[ispec] = Ai[ispec] * mp_g #mass of ion in grams
        mi[ispec] = 1.e-3 * mi_g[ispec]
        qi[ispec] = Zi[ispec] * qe #charge of species in statcoulomb
    end

    nquiet::Int = floor(1.18e-2 / abs(dr))
    cs = 7985.8074573355016 #sqrt(2 * g * maximum(temp0) / mean(mi)) #it's an "average" speed of sound
    println("---WARNING!!!! cs needs to be replaced!!!!!!-----")
    CFL = abs(dr) / cs
    dtm = coeff_CFL * CFL #may need to be varied to avoid NaN when tm > tm_coll
    lm = dtm / dr

    nsmooth::Int32 = floor(abs(dr_smooth/dr)) #only used if smoothing==.true.
    vsmooth = nsmooth  #smoothing for velocity
    tsmooth = nsmooth  #smoothing for temperature

    return dr, mi_g, mi, qi, nquiet, cs, CFL, dtm, lm,
            nsmooth, vsmooth, tsmooth
end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
include("splinesSerial.jl")

#controldata = [[0,1],[0,0],[1,1],[1,0]]
controldata = [[0,1],[0,0],[1,1],[1,0],[0.1,0.1],[0.2,0.2],[0.3,0.3],[0.41,0.4],[0.89,0.73],[0.5,0.5],[0.6,0.6],[0.7,0.7],[0.8,0.8],[0.9,0.9]]
point = [0.5, 0.7, 1.2]
u = S2
larBezier(u, controldata)(point)


text:
module ThermProps

    #using Plots
    #using PyPlot
    #using LaTeXStrings
    #using JLD

    export  LiquidDensity
    export  LiquidViscosity
    export  GasViscosity

    include("Data.jl")
    include("LiquidDensity.jl")
    include("LiquidViscosity.jl")
    include("GasViscosity.jl")
end


text:
"""
    funarray(fun, x, dim; axis)

Helper function to execute a function on each column (row if axis==1) and
return the results.
"""
function funarray( fun::Function, x::AbstractArray{Float},
                   dim::Int, axis::Int=2 )::Union{Float,Vector{Float}}
    res = nothing
    if length(size(x))==1
        res = fun(x)
    else
        @assert axis in [1,2] "unknown axis index"
        @assert size(x)[axis==2?1:2]==dim "dimensions don't match"
        if axis==1
            # read by rows
            res = [fun(x[i,:]) for i in 1:size(x)[1]]
        else
            # read by columns
            res = [fun(x[:,i]) for i in 1:size(x)[2]]
        end
    end
    res
end