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

# magic rounding constant: 1.5*2^52 Adding, then subtracting it from a float rounds it to an Int.
# This works because eps(MAGIC_ROUND_CONST(T)) == one(T), so adding it to a smaller number aligns the lsb to the 1s place.
# Values for which this trick doesn't work are going to have outputs of 0 or Inf.
MAGIC_ROUND_CONST(::Type{Float64}) = 6.755399441055744e15
MAGIC_ROUND_CONST(::Type{Float32}) = 1.048576f7

# max, min, and subnormal arguments
# max_exp = T(exponent_bias(T)*log(base, big(2)) + log(base, 2 - big(2.0)^-significand_bits(T)))
MAX_EXP(n::Val{2}, ::Type{Float32}) = 128.0f0
MAX_EXP(n::Val{2}, ::Type{Float64}) = 1024.0
MAX_EXP(n::Val{:ℯ}, ::Type{Float32}) = 88.72284f0
MAX_EXP(n::Val{:ℯ}, ::Type{Float64}) = 709.7827128933841
MAX_EXP(n::Val{10}, ::Type{Float32}) = 38.53184f0
MAX_EXP(n::Val{10}, ::Type{Float64}) = 308.25471555991675

# min_exp = T(-(exponent_bias(T)+significand_bits(T)) * log(base, big(2)))
MIN_EXP(n::Val{2}, ::Type{Float32}) = -150.0f0
MIN_EXP(n::Val{2}, ::Type{Float64}) = -1075.0
MIN_EXP(n::Val{:ℯ}, ::Type{Float32}) = -103.97208f0
MIN_EXP(n::Val{:ℯ}, ::Type{Float64}) = -745.1332191019412
MIN_EXP(n::Val{10}, ::Type{Float32}) = -45.1545f0
MIN_EXP(n::Val{10}, ::Type{Float64}) = -323.60724533877976

# subnorm_exp = abs(log(base, floatmin(T)))
# these vals are positive since it's easier to take abs(x) than -abs(x)
SUBNORM_EXP(n::Val{2}, ::Type{Float32}) = 126.00001f0
SUBNORM_EXP(n::Val{2}, ::Type{Float64}) = 1022.0
SUBNORM_EXP(n::Val{:ℯ}, ::Type{Float32}) = 87.33655f0
SUBNORM_EXP(n::Val{:ℯ}, ::Type{Float64}) = 708.3964185322641
SUBNORM_EXP(n::Val{10}, ::Type{Float32}) = 37.92978f0
SUBNORM_EXP(n::Val{10}, ::Type{Float64}) = 307.6526555685887

# 256/log(base, 2) (For Float64 reductions)
LogBo256INV(::Val{2}, ::Type{Float64}) = 256.0
LogBo256INV(::Val{:ℯ}, ::Type{Float64}) = 369.3299304675746
LogBo256INV(::Val{10}, ::Type{Float64}) = 850.4135922911647

# -log(base, 2)/256 in upper and lower bits
# Upper is truncated to only have 34 bits of significand since N has at most
# ceil(log2(-MIN_EXP(base, Float64)*LogBo256INV(Val(2), Float64))) = 19 bits.
# This ensures no rounding when multiplying LogBo256U*N for FMAless hardware
LogBo256U(::Val{2}, ::Type{Float64}) = -0.00390625
LogBo256U(::Val{:ℯ}, ::Type{Float64}) = -0.002707606173999011
LogBo256U(::Val{10}, ::Type{Float64}) = -0.0011758984204561784
LogBo256L(::Val{2}, ::Type{Float64}) = 0.0
LogBo256L(::Val{:ℯ}, ::Type{Float64}) = -6.327543041662719e-14
LogBo256L(::Val{10}, ::Type{Float64}) = -1.0624811566412999e-13

# 1/log(base, 2) (For Float32 reductions)
LogBINV(::Val{2}, ::Type{Float32}) = 1.0f0
LogBINV(::Val{:ℯ}, ::Type{Float32}) = 1.442695f0
LogBINV(::Val{10}, ::Type{Float32}) = 3.321928f0

# -log(base, 2) in upper and lower bits
# Upper is truncated to only have 16 bits of significand since N has at most
# ceil(log2(-MIN_EXP(n, Float32)*LogBINV(Val(2), Float32))) = 8 bits.
# This ensures no rounding when multiplying LogBU*N for FMAless hardware
LogBU(::Val{2}, ::Type{Float32}) = -1.0f0
LogBU(::Val{:ℯ}, ::Type{Float32}) = -0.69314575f0
LogBU(::Val{10}, ::Type{Float32}) = -0.3010254f0
LogBL(::Val{2}, ::Type{Float32}) = 0.0f0
LogBL(::Val{:ℯ}, ::Type{Float32}) = -1.4286068f-6
LogBL(::Val{10}, ::Type{Float32}) = -4.605039f-6

# -log(base, 2) as a Float32 for Float16 version.
LogB(::Val{2}, ::Type{Float16}) = -1.0f0
LogB(::Val{:ℯ}, ::Type{Float16}) = -0.6931472f0
LogB(::Val{10}, ::Type{Float16}) = -0.30103f0

# Range reduced kernels
@inline function expm1b_kernel(::Val{2}, x::Float64)
    return x * evalpoly(x, (0.6931471805599393, 0.24022650695910058,
                            0.05550411502333161, 0.009618129548366803))
end
@inline function expm1b_kernel(::Val{:ℯ}, x::Float64)
    return x * evalpoly(x, (0.9999999999999912, 0.4999999999999997,
                            0.1666666857598779, 0.04166666857598777))
end

@inline function expm1b_kernel(::Val{10}, x::Float64)
    return x * evalpoly(x, (2.3025850929940255, 2.6509490552391974,
                            2.034678825384765, 1.1712552025835192))
end

@inline function expb_kernel(::Val{2}, x::Float32)
    return evalpoly(x, (1.0f0, 0.6931472f0, 0.2402265f0,
                        0.05550411f0, 0.009618025f0,
                        0.0013333423f0, 0.00015469732f0, 1.5316464f-5))
end
@inline function expb_kernel(::Val{:ℯ}, x::Float32)
    return evalpoly(x, (1.0f0, 1.0f0, 0.5f0, 0.16666667f0,
                        0.041666217f0, 0.008333249f0,
                        0.001394858f0, 0.00019924171f0))
end
@inline function expb_kernel(::Val{10}, x::Float32)
    return evalpoly(x, (1.0f0, 2.3025851f0, 2.650949f0,
                        2.0346787f0, 1.1712426f0, 0.53937745f0,
                        0.20788547f0, 0.06837386f0))
end

# Table stores data with 60 sig figs by using the fact that the first 12 bits of all the
# values would be the same if stored as regular Float64.
# This only gains 8 bits since the least significant 4 bits of the exponent
# of the small part are not the same for all table entries
const JU_MASK = typemax(UInt64)>>12
const JL_MASK = typemax(UInt64)>>8
const JU_CONST = 0x3FF0000000000000
const JL_CONST = 0x3C00000000000000

#function make_table(size)
#    t_array = zeros(UInt64, size);
#    for j in 1:size
#        val = 2.0^(BigFloat(j-1)/size)
#        valU = Float64(val, RoundDown)
#        valL = Float64(val-valU)
#        valU = reinterpret(UInt64, valU) & JU_MASK
#        valL = ((reinterpret(UInt64, valL) & JL_MASK)>>44)<<52
#        t_array[j] = valU | valL
#    end
#    return Tuple(t_array)
#end
#const J_TABLE = make_table(256);
const J_TABLE = (0x0000000000000000, 0xaac00b1afa5abcbe, 0x9b60163da9fb3335, 0xab502168143b0280, 0xadc02c9a3e778060,
                 0x656037d42e11bbcc, 0xa7a04315e86e7f84, 0x84c04e5f72f654b1, 0x8d7059b0d3158574, 0xa510650a0e3c1f88,
                 0xa8d0706b29ddf6dd, 0x83207bd42b72a836, 0x6180874518759bc8, 0xa4b092bdf66607df, 0x91409e3ecac6f383,
                 0x85d0a9c79b1f3919, 0x98a0b5586cf9890f, 0x94f0c0f145e46c85, 0x9010cc922b7247f7, 0xa210d83b23395deb,
                 0x4030e3ec32d3d1a2, 0xa5b0efa55fdfa9c4, 0xae40fb66affed31a, 0x8d41073028d7233e, 0xa4911301d0125b50,
                 0xa1a11edbab5e2ab5, 0xaf712abdc06c31cb, 0xae8136a814f204aa, 0xa661429aaea92ddf, 0xa9114e95934f312d,
                 0x82415a98c8a58e51, 0x58f166a45471c3c2, 0xab9172b83c7d517a, 0x70917ed48695bbc0, 0xa7718af9388c8de9,
                 0x94a1972658375d2f, 0x8e51a35beb6fcb75, 0x97b1af99f8138a1c, 0xa351bbe084045cd3, 0x9001c82f95281c6b,
                 0x9e01d4873168b9aa, 0xa481e0e75eb44026, 0xa711ed5022fcd91c, 0xa201f9c18438ce4c, 0x8dc2063b88628cd6,
                 0x935212be3578a819, 0x82a21f49917ddc96, 0x8d322bdda27912d1, 0x99b2387a6e756238, 0x8ac2451ffb82140a,
                 0x8ac251ce4fb2a63f, 0x93e25e85711ece75, 0x82b26b4565e27cdd, 0x9e02780e341ddf29, 0xa2d284dfe1f56380,
                 0xab4291ba7591bb6f, 0x86129e9df51fdee1, 0xa352ab8a66d10f12, 0xafb2b87fd0dad98f, 0xa572c57e39771b2e,
                 0x9002d285a6e4030b, 0x9d12df961f641589, 0x71c2ecafa93e2f56, 0xaea2f9d24abd886a, 0x86f306fe0a31b715,
                 0x89531432edeeb2fd, 0x8a932170fc4cd831, 0xa1d32eb83ba8ea31, 0x93233c08b26416ff, 0xab23496266e3fa2c,
                 0xa92356c55f929ff0, 0xa8f36431a2de883a, 0xa4e371a7373aa9ca, 0xa3037f26231e7549, 0xa0b38cae6d05d865,
                 0xa3239a401b7140ee, 0xad43a7db34e59ff6, 0x9543b57fbfec6cf4, 0xa083c32dc313a8e4, 0x7fe3d0e544ede173,
                 0x8ad3dea64c123422, 0xa943ec70df1c5174, 0xa413fa4504ac801b, 0x8bd40822c367a024, 0xaf04160a21f72e29,
                 0xa3d423fb27094689, 0xab8431f5d950a896, 0x88843ffa3f84b9d4, 0x48944e086061892d, 0xae745c2042a7d231,
                 0x9c946a41ed1d0057, 0xa1e4786d668b3236, 0x73c486a2b5c13cd0, 0xab1494e1e192aed1, 0x99c4a32af0d7d3de,
                 0xabb4b17dea6db7d6, 0x7d44bfdad5362a27, 0x9054ce41b817c114, 0x98e4dcb299fddd0d, 0xa564eb2d81d8abfe,
                 0xa5a4f9b2769d2ca6, 0x7a2508417f4531ee, 0xa82516daa2cf6641, 0xac65257de83f4eee, 0xabe5342b569d4f81,
                 0x879542e2f4f6ad27, 0xa8a551a4ca5d920e, 0xa7856070dde910d1, 0x99b56f4736b527da, 0xa7a57e27dbe2c4ce,
                 0x82958d12d497c7fd, 0xa4059c0827ff07cb, 0x9635ab07dd485429, 0xa245ba11fba87a02, 0x3c45c9268a5946b7,
                 0xa195d84590998b92, 0x9ba5e76f15ad2148, 0xa985f6a320dceb70, 0xa60605e1b976dc08, 0x9e46152ae6cdf6f4,
                 0xa636247eb03a5584, 0x984633dd1d1929fd, 0xa8e6434634ccc31f, 0xa28652b9febc8fb6, 0xa226623882552224,
                 0xa85671c1c70833f5, 0x60368155d44ca973, 0x880690f4b19e9538, 0xa216a09e667f3bcc, 0x7a36b052fa75173e,
                 0xada6c012750bdabe, 0x9c76cfdcddd47645, 0xae46dfb23c651a2e, 0xa7a6ef9298593ae4, 0xa9f6ff7df9519483,
                 0x59d70f7466f42e87, 0xaba71f75e8ec5f73, 0xa6f72f8286ead089, 0xa7a73f9a48a58173, 0x90474fbd35d7cbfd,
                 0xa7e75feb564267c8, 0x9b777024b1ab6e09, 0x986780694fde5d3f, 0x934790b938ac1cf6, 0xaaf7a11473eb0186,
                 0xa207b17b0976cfda, 0x9f17c1ed0130c132, 0x91b7d26a62ff86f0, 0x7057e2f336cf4e62, 0xabe7f3878491c490,
                 0xa6c80427543e1a11, 0x946814d2add106d9, 0xa1582589994cce12, 0x9998364c1eb941f7, 0xa9c8471a4623c7ac,
                 0xaf2857f4179f5b20, 0xa01868d99b4492ec, 0x85d879cad931a436, 0x99988ac7d98a6699, 0x9d589bd0a478580f,
                 0x96e8ace5422aa0db, 0x9ec8be05bad61778, 0xade8cf3216b5448b, 0xa478e06a5e0866d8, 0x85c8f1ae99157736,
                 0x959902fed0282c8a, 0xa119145b0b91ffc5, 0xab2925c353aa2fe1, 0xae893737b0cdc5e4, 0xa88948b82b5f98e4,
                 0xad395a44cbc8520e, 0xaf296bdd9a7670b2, 0xa1797d829fde4e4f, 0x7ca98f33e47a22a2, 0xa749a0f170ca07b9,
                 0xa119b2bb4d53fe0c, 0x7c79c49182a3f090, 0xa579d674194bb8d4, 0x7829e86319e32323, 0xaad9fa5e8d07f29d,
                 0xa65a0c667b5de564, 0x9c6a1e7aed8eb8bb, 0x963a309bec4a2d33, 0xa2aa42c980460ad7, 0xa16a5503b23e255c,
                 0x650a674a8af46052, 0x9bca799e1330b358, 0xa58a8bfe53c12e58, 0x90fa9e6b5579fdbf, 0x889ab0e521356eba,
                 0xa81ac36bbfd3f379, 0x97ead5ff3a3c2774, 0x97aae89f995ad3ad, 0xa5aafb4ce622f2fe, 0xa21b0e07298db665,
                 0x94db20ce6c9a8952, 0xaedb33a2b84f15fa, 0xac1b468415b749b0, 0xa1cb59728de55939, 0x92ab6c6e29f1c52a,
                 0xad5b7f76f2fb5e46, 0xa24b928cf22749e3, 0xa08ba5b030a10649, 0xafcbb8e0b79a6f1e, 0x823bcc1e904bc1d2,
                 0xafcbdf69c3f3a206, 0xa08bf2c25bd71e08, 0xa89c06286141b33c, 0x811c199bdd85529c, 0xa48c2d1cd9fa652b,
                 0x9b4c40ab5fffd07a, 0x912c544778fafb22, 0x928c67f12e57d14b, 0xa86c7ba88988c932, 0x71ac8f6d9406e7b5,
                 0xaa0ca3405751c4da, 0x750cb720dcef9069, 0xac5ccb0f2e6d1674, 0xa88cdf0b555dc3f9, 0xa2fcf3155b5bab73,
                 0xa1ad072d4a07897b, 0x955d1b532b08c968, 0xa15d2f87080d89f1, 0x93dd43c8eacaa1d6, 0x82ed5818dcfba487,
                 0x5fed6c76e862e6d3, 0xa77d80e316c98397, 0x9a0d955d71ff6075, 0x9c2da9e603db3285, 0xa24dbe7cd63a8314,
                 0x92ddd321f301b460, 0xa1ade7d5641c0657, 0xa72dfc97337b9b5e, 0xadae11676b197d16, 0xa42e264614f5a128,
                 0xa30e3b333b16ee11, 0x839e502ee78b3ff6, 0xaa7e653924676d75, 0x92de7a51fbc74c83, 0xa77e8f7977cdb73f,
                 0xa0bea4afa2a490d9, 0x948eb9f4867cca6e, 0xa1becf482d8e67f0, 0x91cee4aaa2188510, 0x9dcefa1bee615a27,
                 0xa66f0f9c1cb64129, 0x93af252b376bba97, 0xacdf3ac948dd7273, 0x99df50765b6e4540, 0x9faf6632798844f8,
                 0xa12f7bfdad9cbe13, 0xaeef91d802243c88, 0x874fa7c1819e90d8, 0xacdfbdba3692d513, 0x62efd3c22b8f71f1, 0x74afe9d96b2a23d9)

# XXX we want to mark :consistent-cy here so that this function can be concrete-folded,
# because the effect analysis currently can't prove it in the presence of `@inbounds` or
# `:boundscheck`, but still the access to `J_TABLE` is really safe here
Base.@assume_effects :consistent @inline function table_unpack(ind::Int32)
    ind = ind & 255 + 1 # 255 == length(J_TABLE) - 1
    j = @inbounds J_TABLE[ind]
    jU = reinterpret(Float64, JU_CONST | (j&JU_MASK))
    jL = reinterpret(Float64, JL_CONST | (j>>8))
    return jU, jL
end

# Method for Float64
# 1. Argument reduction: Reduce x to an r so that |r| <= log(b, 2)/512. Given x, base b,
#    find r and integers k, j such that
#       x = (k + j/256)*log(b, 2) + r,  0 <= j < 256, |r| <= log(b,2)/512.
#
# 2. Approximate b^r-1 by 3rd-degree minimax polynomial p_b(r) on the interval [-log(b,2)/512, log(b,2)/512].
#    Since the bounds on r are very tight, this is sufficient to be accurate to floating point epsilon.
#
# 3. Scale back: b^x = 2^k * 2^(j/256) * (1 + p_b(r))
#    Since the range of possible j is small, 2^(j/256) is stored for all possible values in slightly extended precision.

# Method for Float32
# 1. Argument reduction: Reduce x to an r so that |r| <= log(b, 2)/2. Given x, base b,
#    find r and integer N such that
#       x = N*log(b, 2) + r,  |r| <= log(b,2)/2.
#
# 2. Approximate b^r by 7th-degree minimax polynomial p_b(r) on the interval [-log(b,2)/2, log(b,2)/2].
# 3. Scale back: b^x = 2^N * p_b(r)
# For both, a little extra care needs to be taken if b^r is subnormal.
# The solution is to do the scaling back in 2 steps as just messing with the exponent wouldn't work.

@inline function exp_impl(x::Float64, base)
    T = Float64
    N_float = muladd(x, LogBo256INV(base, T), MAGIC_ROUND_CONST(T))
    N = reinterpret(UInt64, N_float) % Int32
    N_float -=  MAGIC_ROUND_CONST(T) #N_float now equals round(x*LogBo256INV(base, T))
    r = muladd(N_float, LogBo256U(base, T), x)
    r = muladd(N_float, LogBo256L(base, T), r)
    k = N >> 8
    jU, jL = table_unpack(N)
    small_part =  muladd(jU, expm1b_kernel(base, r), jL) + jU

    if !(abs(x) <= SUBNORM_EXP(base, T))
        x >= MAX_EXP(base, T) && return Inf
        x <= MIN_EXP(base, T) && return 0.0
        if k <= -53
            # The UInt64 forces promotion. (Only matters for 32 bit systems.)
            twopk = (k + UInt64(53)) << 52
            return reinterpret(T, twopk + reinterpret(UInt64, small_part))*(2.0^-53)
        end
        #k == 1024 && return (small_part * 2.0) * 2.0^1023
    end
    twopk = Int64(k) << 52
    return reinterpret(T, twopk + reinterpret(Int64, small_part))
end
# Computes base^(x+xlo). Used for pow.
@inline function exp_impl(x::Float64, xlo::Float64, base)
    T = Float64
    N_float = muladd(x, LogBo256INV(base, T), MAGIC_ROUND_CONST(T))
    N = reinterpret(UInt64, N_float) % Int32
    N_float -=  MAGIC_ROUND_CONST(T) #N_float now equals round(x*LogBo256INV(base, T))
    r = muladd(N_float, LogBo256U(base, T), x)
    r = muladd(N_float, LogBo256L(base, T), r)
    k = N >> 8
    jU, jL = table_unpack(N)
    very_small = muladd(jU, expm1b_kernel(base, r), jL)
    small_part =  muladd(jU,xlo,very_small) + jU
    if !(abs(x) <= SUBNORM_EXP(base, T))
        x >= MAX_EXP(base, T) && return Inf
        x <= MIN_EXP(base, T) && return 0.0
        if k <= -53
            # The UInt64 forces promotion. (Only matters for 32 bit systems.)
            twopk = (k + UInt64(53)) << 52
            return reinterpret(T, twopk + reinterpret(UInt64, small_part))*(2.0^-53)
        end
        #k == 1024 && return (small_part * 2.0) * 2.0^1023
    end
    twopk = Int64(k) << 52
    return reinterpret(T, twopk + reinterpret(Int64, small_part))
end
@inline function exp_impl_fast(x::Float64, base)
    T = Float64
    x >= MAX_EXP(base, T) && return Inf
    x <= -SUBNORM_EXP(base, T) && return 0.0
    N_float = muladd(x, LogBo256INV(base, T), MAGIC_ROUND_CONST(T))
    N = reinterpret(UInt64, N_float) % Int32
    N_float -=  MAGIC_ROUND_CONST(T) #N_float now equals round(x*LogBo256INV(base, T))
    r = muladd(N_float, LogBo256U(base, T), x)
    r = muladd(N_float, LogBo256L(base, T), r)
    k = N >> 8
    jU = reinterpret(Float64, JU_CONST | (@inbounds J_TABLE[N&255 + 1] & JU_MASK))
    small_part =  muladd(jU, expm1b_kernel(base, r), jU)
    twopk = Int64(k) << 52
    return reinterpret(T, twopk + reinterpret(Int64, small_part))
end

@inline function exp_impl(x::Float32, base)
    T = Float32
    N_float = round(x*LogBINV(base, T))
    N = unsafe_trunc(Int32, N_float)
    r = muladd(N_float, LogBU(base, T), x)
    r = muladd(N_float, LogBL(base, T), r)
    small_part = expb_kernel(base, r)
    power = (N+Int32(127))
    x > MAX_EXP(base, T) && return Inf32
    x < MIN_EXP(base, T) && return 0.0f0
    if x <= -SUBNORM_EXP(base, T)
        power += Int32(24)
        small_part *= Float32(0x1p-24)
    end
    if N == 128
        power -= Int32(1)
        small_part *= 2f0
    end
    return small_part * reinterpret(T, power << Int32(23))
end

@inline function exp_impl_fast(x::Float32, base)
    T = Float32
    x >= MAX_EXP(base, T) && return Inf32
    x <= -SUBNORM_EXP(base, T) && return 0f0
    N_float = round(x*LogBINV(base, T))
    N = unsafe_trunc(Int32, N_float)
    r = muladd(N_float, LogBU(base, T), x)
    r = muladd(N_float, LogBL(base, T), r)
    small_part = expb_kernel(base, r)
    twopk = reinterpret(T, (N+Int32(127)) << Int32(23))
    return twopk*small_part
end

@inline function exp_impl(a::Float16, base)
    T = Float32
    x = T(a)
    N_float = round(x*LogBINV(base, T))
    N = unsafe_trunc(Int32, N_float)
    r = muladd(N_float, LogB(base, Float16), x)
    small_part = expb_kernel(base, r)
    if !(abs(x) <= 25)
        x > 16 && return Inf16
        x < 25 && return zero(Float16)
    end
    twopk = reinterpret(T, (N+Int32(127)) << Int32(23))
    return Float16(twopk*small_part)
end

for (func, fast_func, base) in ((:exp2,  :exp2_fast,  Val(2)),
                                (:exp,   :exp_fast,   Val(:ℯ)),
                                (:exp10, :exp10_fast, Val(10)))
    @eval begin
        $func(x::Union{Float16,Float32,Float64}) = exp_impl(x, $base)
        $fast_func(x::Union{Float32,Float64}) = exp_impl_fast(x, $base)
    end
end

@doc """
    exp(x)

Compute the natural base exponential of `x`, in other words ``ℯ^x``.

See also [`exp2`](@ref), [`exp10`](@ref) and [`cis`](@ref).

# Examples
```jldoctest
julia> exp(1.0)
2.718281828459045

julia> exp(im * pi) ≈ cis(pi)
true
```
""" exp(x::Real)

"""
    exp2(x)

Compute the base 2 exponential of `x`, in other words ``2^x``.

See also [`ldexp`](@ref), [`<<`](@ref).

# Examples
```jldoctest
julia> exp2(5)
32.0

julia> 2^5
32

julia> exp2(63) > typemax(Int)
true
```
"""
exp2(x)

"""
    exp10(x)

Compute the base 10 exponential of `x`, in other words ``10^x``.

# Examples
```jldoctest
julia> exp10(2)
100.0

julia> 10^2
100
```
"""
exp10(x)

# functions with special cases for integer arguments
@inline function exp2(x::Base.BitInteger)
    if x > 1023
        Inf64
    elseif x <= -1023
        # if -1073 < x <= -1023 then Result will be a subnormal number
        # Hex literal with padding must be used to work on 32bit machine
        reinterpret(Float64, 0x0000_0000_0000_0001 << ((x + 1074) % UInt))
    else
        # We will cast everything to Int64 to avoid errors in case of Int128
        # If x is a Int128, and is outside the range of Int64, then it is not -1023<x<=1023
        reinterpret(Float64, (exponent_bias(Float64) + (x % Int64)) << (significand_bits(Float64) % UInt))
    end
end

# min and max arguments for expm1 by type
MAX_EXP(::Type{Float64}) =  709.7827128933845   # log 2^1023*(2-2^-52)
MIN_EXP(::Type{Float64}) = -37.42994775023705   # log 2^-54
MAX_EXP(::Type{Float32}) =  88.72284f0          # log 2^127 *(2-2^-23)
MIN_EXP(::Type{Float32}) = -17.32868f0          # log 2^-25
MAX_EXP(::Type{Float16}) =  Float16(11.09)      # log 2^15 *(2-2^-10)
MIN_EXP(::Type{Float16}) = -Float16(8.32)       # log 2^-12

Ln2INV(::Type{Float64}) = 1.4426950408889634
Ln2(::Type{Float64}) = -0.6931471805599453
Ln2INV(::Type{Float32}) = 1.442695f0
Ln2(::Type{Float32}) = -0.6931472f0

# log(.75) <= x <= log(1.25)
@inline function expm1_small(x::Float64)
    p = evalpoly(x, (0.16666666666666632, 0.04166666666666556, 0.008333333333401227,
                     0.001388888889068783, 0.00019841269447671544, 2.480157691845342e-5,
                     2.7558212415361945e-6, 2.758218402815439e-7, 2.4360682937111612e-8))
    p2 = exthorner(x, (1.0, .5, p))
    return fma(x, p2[1], x*p2[2])
end
@inline function expm1_small(x::Float32)
    p = evalpoly(x, (0.16666666f0, 0.041666627f0, 0.008333682f0,
                     0.0013908712f0, 0.0001933096f0))
    p2 = exthorner(x, (1f0, .5f0, p))
    return fma(x, p2[1], x*p2[2])
end

function expm1(x::Float64)
    T = Float64
    if -0.2876820724517809 <= x <= 0.22314355131420976
        return expm1_small(x)
    elseif !(abs(x)<=MIN_EXP(Float64))
        isnan(x) && return x
        x > MAX_EXP(Float64) && return Inf
        x < MIN_EXP(Float64) && return -1.0
    end

    N_float = muladd(x, LogBo256INV(Val(:ℯ), T), MAGIC_ROUND_CONST(T))
    N = reinterpret(UInt64, N_float) % Int32
    N_float -=  MAGIC_ROUND_CONST(T) #N_float now equals round(x*LogBo256INV(Val(:ℯ), T))
    r = muladd(N_float, LogBo256U(Val(:ℯ), T), x)
    r = muladd(N_float, LogBo256L(Val(:ℯ), T), r)
    k = Int64(N >> 8)
    jU, jL = table_unpack(N)
    p = expm1b_kernel(Val(:ℯ), r)
    twopk  = reinterpret(Float64, (1023+k) << 52)
    twopnk = reinterpret(Float64, (1023-k) << 52)
    k>=106 && return reinterpret(Float64, (1022+k) << 52)*(jU + muladd(jU, p, jL))*2
    k>=53 && return twopk*(jU + muladd(jU, p, (jL-twopnk)))
    k<=-2 && return twopk*(jU + muladd(jU, p, jL))-1
    return twopk*((jU-twopnk) + fma(jU, p, jL))
end

function expm1(x::Float32)
    x > MAX_EXP(Float32) && return Inf32
    x < MIN_EXP(Float32) && return -1f0
    if -0.2876821f0 <=x <= 0.22314355f0
        return expm1_small(x)
    end
    x = Float64(x)
    N_float = round(x*Ln2INV(Float64))
    N = unsafe_trunc(UInt64, N_float)
    r = muladd(N_float, Ln2(Float64), x)
    hi = evalpoly(r, (1.0, .5, 0.16666667546642386, 0.041666183019487026,
                      0.008332997481506921, 0.0013966479175977883, 0.0002004037059220124))
    small_part = r*hi
    twopk = reinterpret(Float64, (N+1023) << 52)
    return Float32(muladd(twopk, small_part, twopk-1.0))
end

function expm1(x::Float16)
    x > MAX_EXP(Float16) && return Inf16
    x < MIN_EXP(Float16) && return Float16(-1.0)
    x = Float32(x)
    if -0.2876821f0 <=x <= 0.22314355f0
        return Float16(x*evalpoly(x, (1f0, .5f0, 0.16666628f0, 0.04166785f0, 0.008351848f0, 0.0013675707f0)))
    end
    N_float = round(x*Ln2INV(Float32))
    N = unsafe_trunc(UInt32, N_float)
    r = muladd(N_float, Ln2(Float32), x)
    hi = evalpoly(r, (1f0, .5f0, 0.16666667f0, 0.041665863f0, 0.008333111f0, 0.0013981499f0, 0.00019983904f0))
    small_part = r*hi
    twopk = reinterpret(Float32, (N+Int32(127)) << Int32(23))
    return Float16(muladd(twopk, small_part, twopk-1f0))
end

"""
    expm1(x)

Accurately compute ``e^x-1``. It avoids the loss of precision involved in the direct
evaluation of exp(x)-1 for small values of x.
# Examples
```jldoctest
julia> expm1(1e-16)
1.0e-16

julia> exp(1e-16) - 1
0.0
```
"""
expm1(x)

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: Nexuscompute/julia


text: Dict{String, Any}("title" => Any[Dict{String, Any}("value" => " \", ", "type" => "string")])

ext: jl


text:
import JSON
import MCTS: AbstractTreeVisualizer, node_tag, tooltip_tag, create_json, blink

mutable struct DESPOTVisualizer <: AbstractTreeVisualizer
    root::VNode
end

DESPOTVisualizer(solver::DESPOTSolver) = DESPOTVisualizer(solver.root)

blink(solver::DESPOTSolver) = blink(DESPOTVisualizer(solver))

const NodeDict = Dict{Int, Dict{String, Any}}

function create_json(v::DESPOTVisualizer)
    node_dict = NodeDict()
    dict = recursive_push!(node_dict, v.root, "root (Q values at action nodes are upper bounds)")
    json = JSON.json(node_dict)
    return (json, 1)
end

function recursive_push!(nd::NodeDict, n::VNode, obs, parent_id=-1)
    id = length(nd) + 1
    if parent_id > 0
        push!(nd[parent_id]["children_ids"], id)
    end
    @assert n.n_visits==0
    nd[id] = Dict("id"=>id,
                  "type"=>:V,
                  "children_ids"=>Array{Int}(0),
                  "tag"=>node_tag(obs),
                  "tt_tag"=>tooltip_tag(obs),
                  "N"=>length(n.particles)
                  )
    for (a,c) in n.q_nodes
        recursive_push!(nd, c, id)
    end
    return nd
end

function recursive_push!{S,A,O,B}(nd::NodeDict, n::QNode{S,A,O,B}, parent_id=-1)
    id = length(nd) + 1
    if parent_id > 0
        push!(nd[parent_id]["children_ids"], id)
    end
    nd[id] = Dict("id"=>id,
                  "type"=>:action,
                  "children_ids"=>Array{Int}(0),
                  "tag"=>node_tag(n.action),
                  "tt_tag"=>tooltip_tag(n.action),
                  "N"=>sum(length(last(pair).particles) for pair in n.obs_and_nodes),
                  # "Q"=>"between $(get_lower_bound(n)) and $(get_upper_bound(n))"
                  "Q"=>get_upper_bound(n)
                  )
    for (o,c) in n.obs_and_nodes
        recursive_push!(nd, c, o, id)
    end
    return nd
end


text:
@testset "Result" begin
    @testset "Result+PathResult" begin
        R = solve(equations(heart()), save_all_paths=true, show_progress=false)
        @test R isa Result
        @test nat_infinity(R) ≤ 572
        @test nfinite(R) == 4
        @test length(collect(R)) == 576
        @test finite(R) isa Vector{<:PathResult}

        @test length(finite(R, only_nonsingular=false)) == 4
        @test length(finite(R, only_nonsingular=true)) == 4
        @test length(finite(R, only_singular=true)) == 0
        @test 572 - length(failed(R)) == nat_infinity(R)
        @test length(real(R, tol=1e-6)) == 2
        @test nreal(R, tol=1e-6) == 2
        @test length(real_solutions(R)) == 2
        @test_deprecated realsolutions(R)
        @test length(at_infinity(R)) ≤ 572
        @test length(results(R, only_real=true, real_tol=1e-8)) == 2
        @test length(results(R, only_nonsingular=true, singular_tol=1e9)) == 4
        @test length(finite(results(R, only_real=true))) == 2
        @test nresults(R, only_nonsingular=true, singular_tol=1e9) == 4
        @test length(results(R, onlyfinite=false)) == 576
        @test nresults(R, onlyfinite=false) == 576
        @test nnonsingular(R) == 4
        @test length(nonsingular(R)) == 4
        @test length(singular(R)) == 0
        @test length(singular(real(R))) == 0

        @test_nowarn mapresults(solution, R)
        @test_nowarn mapresults(start_solution, R)
        @test_nowarn mapresults(residual, R)
        @test_nowarn mapresults(is_singular, R)
        @test count(mapresults(is_real, R)) == 2
        # test fallback
        @test count(results(is_real, R)) == 2

        @test length(solutions(R)) == 4
        @test solutions(R) isa Vector{Vector{ComplexF64}}
        @test real_solutions(R, real_tol=1e-8) isa Vector{Vector{Float64}}
        @test length(real_solutions(R, real_tol=1e-8)) == 2

        @test_nowarn string(R)
        @test_nowarn string(R[end])
        @test_nowarn show(IOContext(devnull, :compact=>true), R)
        @test_nowarn show(IOContext(devnull, :compact=>true), R[end])

        test_treeviews(R)

        @polyvar x y z
        R = solve([(x-3)^3,(y-2)], affine_tracking=true)
        @test R isa Result
        test_treeviews(R)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 1
        @test nsingular(R, counting_multiplicities=true) == 3
        @test_nowarn sprint(show, R)

        R = solve([(x-3)^3,(y-2)], affine_tracking=true, system_scaling=nothing)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 1
        @test nsingular(R, counting_multiplicities=true) == 3
        @test_nowarn sprint(show, R)

        @polyvar x y z
        R = solve([(x-3z),(y-2z)])
        @test R isa Result{<:ProjectiveVectors.PVector}
        @test_nowarn string(R)
        test_treeviews(R)
        @test nnonsingular(R) == 1

        @polyvar x y z
        R = solve([(x-3z)^3,(y-2z)])
        @test R isa Result{<:ProjectiveVectors.PVector}
        test_treeviews(R)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 1
        @test nsingular(R, counting_multiplicities=true) == 3

        @polyvar x y z
        V = x^2 + y^2 - z^2
        L = randn(1,3) * [x, y, z]
        R = solve([V; L])
        @test_nowarn sprint(show, R)
    end

    @testset "singular result" begin
        @polyvar x
        f = (x-3)^3*(x-2)
        R = solve([f])
        @test length(solutions(R)) == 2
        @test sort(results(multiplicity, R)) == [1,3]

        @polyvar x y
        f = (x-3y)^3*(x-2y)
        R = solve([f])
        @test nsolutions(R) == 2
        @test sort(results(multiplicity, R)) == [1,3]
    end

    @testset "multiplicities" begin
        @polyvar x y z
        f = (x-1*y)*(x-1.01*y)
        g = x - y - z
        S = solve([f,g], threading=false)
        @test length(multiplicities(S, tol=1e-1)) == 1
        @test length(multiplicities(S, tol=1e-5)) == 0
        @test length(multiplicities(S.pathresults, tol=1e-1)) == 1
        @test length(multiplicities(S.pathresults, tol=1e-5)) == 0

        @polyvar x z
        f = (x-1)*(x-1.01)
        g = x - 1 - z
        S = solve([f,g], threading=false)
        @test length(multiplicities(S, tol=1e-1)) == 1
        @test length(multiplicities(S, tol=1e-5)) == 0
        @test length(multiplicities(S.pathresults, tol=1e-1)) == 1
        @test length(multiplicities(S.pathresults, tol=1e-5)) == 0
        @test norm(solution(S[1]) - solution(S[2])) < 1e-1
    end
end


text:
##

using ACE
using StaticArrays, Random, Printf, Test, LinearAlgebra, ACE.Testing
using ACE: evaluate, evaluate_d, SymmetricBasis, PIBasis, O3
using ACE.Random: rand_rot, rand_refl
using ACEbase.Testing: fdtest, println_slim
using ACE.Testing: __TestSVec

# using Profile, ProfileView

# Extra using Wigner for computing Wigner Matrix
using ACE.Wigner
using ACE.Wigner: get_orbsym

# construct the 1p-basis
maxdeg = 6
ord = 3
Bsel = SimpleSparseBasis(ord, maxdeg)

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)

# generate a configuration
nX = 10
Xs = rand(PositionState{Float64}, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

##

@info("SymmetricBasis construction and evaluation: Invariant Scalar")

φ = ACE.Invariant()
basis = SymmetricBasis(φ, B1p, Bsel)
@time SymmetricBasis(φ, B1p, Bsel);

BB = evaluate(basis, cfg)

# a stupid but necessary test
AA = evaluate(basis.pibasis, cfg)
BB1 = basis.A2Bmap * AA
println_slim(@test isapprox(BB, BB1, rtol=1e-10))

@info("evaluate with vector vs config")
println_slim(@test BB ≈ evaluate(basis, Xs))
println_slim(@test evaluate_d(basis, cfg) ≈ evaluate_d(basis, Xs))

# check there are no superfluous columns
Iz = findall(iszero, sum(norm, basis.A2Bmap, dims=1)[:])
if !isempty(Iz)
   @warn("The A2B map for Invariants has $(length(Iz))/$(length(basis.pibasis)) zero-columns!!!!")
end

for ntest = 1:30
      Xs1 = shuffle(rand_refl(rand_rot(Xs)))
      local BB1 = evaluate(basis, ACEConfig(Xs1))
      print_tf(@test isapprox(BB, BB1, rtol=1e-10))
end
println()

##

@info("Test Invariant Basis without Constant")
basis_noc = SymmetricBasis(φ, B1p, Bsel; filterfun = ACE.NoConstant())
println_slim(@test length(basis_noc) == length(basis) - 1)
BB_noc = evaluate(basis_noc, cfg)
println_slim(@test (BB_noc == BB[2:end]))

##

@info("Test what happens with an empty configuration")

Xs_empty = Vector{eltype(Xs)}(undef, 0)
cfg_empty = ACEConfig(Xs_empty)
B_empty = evaluate(basis, cfg_empty)
println(@test( all(iszero, B_empty[2:end]) ))
println(@test( B_empty[1] == ACE.Invariant(1.0) ) )

##
import ACEbase
@info("Test FIO")
let basis1 = basis
   println_slim(@test(all(ACEbase.Testing.test_fio(basis1; warntype=true))))
end

##

@info("Test linear independence of the basis")
# generate some random configurations; ord^2 + 1 sounds good :)
cfgs = [ ACEConfig(rand(PositionState{Float64}, B1p.bases[1], nX))
         for _ = 1:(3*length(basis)) ]
A = zeros(length(cfgs), length(basis))
for (i, cfg) in enumerate(cfgs)
   A[i, :] = getproperty.(evaluate(basis, cfg), :val)
end
println_slim(@test rank(A) == length(basis))

# ## Keep for futher profiling
# φ = ACE.Invariant()
# pibasis = PIBasis(B1p, ord, maxdeg; property = φ)
# basis = SymmetricBasis(pibasis, φ)
# @time SymmetricBasis(pibasis, φ);
#
# Profile.clear()
# @profile SymmetricBasis(pibasis, φ);
# ProfileView.view()
#

## Testing derivatives

for ntest = 1:30
   _rrval(x::ACE.XState) = x.rr
   Us = randn(SVector{3, Float64}, length(Xs))
   c = randn(length(basis))
   F = t -> sum(c .* ACE.evaluate(basis, ACEConfig(Xs + t[1] * Us))).val
   dF = t -> [ Us' * _rrval.(sum(c .* ACE.evaluate_d(basis, ACEConfig(Xs + t[1] * Us)), dims=1)[:]) ]
   print_tf(@test fdtest(F, dF, [0.0], verbose=false))
end
println()

##
@info("SymmetricBasis construction and evaluation: Spherical Vector")

for L = 0:3
   @info "Tests for L = $L ⇿ $(get_orbsym(0))-$(get_orbsym(L)) block"
   local φ, pibasis, basis, BB, Iz
   φ = ACE.SphericalVector(L; T = ComplexF64)
   basis = SymmetricBasis(φ, B1p, Bsel)
   @time SymmetricBasis(φ, B1p, Bsel)
   BB = evaluate(basis, cfg)

   @info("check that no-const gives the same (sub-)basis")
   basis0 = SymmetricBasis(φ, B1p, Bsel; filterfun = ACE.NoConstant())
   BB0 = evaluate(basis0, cfg)
   if L == 0
      println_slim(@test BB[2:end] ≈ BB0)
   else
      println_slim(@test BB ≈ BB0)
   end

   Iz = findall(iszero, sum(norm, basis.A2Bmap, dims = 1))
   if !isempty(Iz)
      @warn("The A2B map for SphericalVector has $(length(Iz))/$(length(basis.pibasis)) zero-columns!!!!")
   end

   @info("check symmetry ")
   for ntest = 1:30
      local Q, D, BB1
      Q, D = ACE.Wigner.rand_QD(L)
      cfg1 = ACEConfig( shuffle(Ref(Q) .* Xs) )
      BB1 = evaluate(basis, cfg1)
      DtxBB1 = Ref(D') .* BB1
      print_tf(@test isapprox(DtxBB1, BB, rtol=1e-10))
   end
   println()

   @info(" .... derivatives")
   for ntest = 1:30
      _rrval(x::ACE.XState) = x.rr
      Us = __TestSVec.(randn(SVector{3, Float64}, length(Xs)))
      C = randn(typeof(φ.val), length(basis))
      F = t -> sum( sum(c .* b.val)
                    for (c, b) in zip(C, ACE.evaluate(basis, ACEConfig(Xs + t[1] * Us))) )
      dF = t -> [ sum( sum(c .* db)
                  for (c, db) in zip(C, _rrval.(ACE.evaluate_d(basis, ACEConfig(Xs + t[1] * Us))) * Us) ) ]
      print_tf(@test fdtest(F, dF, [0.0], verbose=false))
   end
   println()
end

# ## Keep for futher profiling
# L = 1
# φ = ACE.SphericalVector(L; T = ComplexF64)
# pibasis = PIBasis(B1p, 4, 8; property = φ, isreal = false)
# basis = SymmetricBasis(pibasis, φ)
# @time SymmetricBasis(pibasis, φ);
#
# Profile.clear(); # Profile.init(; delay = 0.0001)
# @profile SymmetricBasis(pibasis, φ);
# ProfileView.view()

##

@info("SymmetricBasis construction and evaluation: Spherical Matrix")

for L1 = 0:2, L2 = 0:2
   @info "Tests for L₁ = $L1, L₂ = $L2 ⇿ $(get_orbsym(L1))-$(get_orbsym(L2)) block"
   local φ, pibasis, basis, BB, Iz
   φ = ACE.SphericalMatrix(L1, L2; T = ComplexF64)
   basis = SymmetricBasis(φ, B1p, Bsel)
   @time basis = SymmetricBasis(φ, B1p, Bsel)
   BB = evaluate(basis, cfg)

   @info("check that no-const gives the same (sub-)basis")
   basis0 = SymmetricBasis(φ, B1p, Bsel; filterfun = ACE.NoConstant())
   BB0 = evaluate(basis0, cfg)
   if L1 == L2
      println_slim(@test BB[2:end] ≈ BB0)
   else
      println_slim(@test BB ≈ BB0)
   end

   @info("check symmetry ")
   for ntest = 1:30
      local Q, D1, D2, BB1
      Q, D1, D2 = ACE.Wigner.rand_QD(L1, L2)
      cfg1 = ACEConfig( shuffle(Ref(Q) .* Xs) )
      BB1 = evaluate(basis, cfg1)
      D1txBB1xD2 = Ref(D1') .* BB1 .* Ref(D2)
      print_tf(@test isapprox(D1txBB1xD2, BB, rtol=1e-10))
   end
   println()

   @info(" .... derivatives")
   for ntest = 1:30
      _rrval(x::ACE.XState) = x.rr
      Us = __TestSVec.(randn(SVector{3, Float64}, length(Xs)))
      C = randn(typeof(φ.val), length(basis))
      F = t -> sum( sum(c .* b.val)
                    for (c, b) in zip(C, ACE.evaluate(basis, ACEConfig(Xs + t[1] * Us))) )
      dF = t -> [ sum( sum(c .* db)
                       for (c, db) in zip(C, _rrval.(ACE.evaluate_d(basis, ACEConfig(Xs + t[1] * Us))) * Us) ) ]
      print_tf(@test fdtest(F, dF, [0.0], verbose=false))
   end
   println()
end

##
@info("Consistency between SphericalVector & SphericalMatrix")

for L = 0:3
   @info "L = $L"
   local Xs, cfg
   φ1 = ACE.SphericalVector(L; T = ComplexF64)
   basis1 = SymmetricBasis(φ1, B1p, Bsel)

   φ2 = ACE.SphericalMatrix(L, 0; T = ComplexF64)
   basis2 = SymmetricBasis(φ2, B1p, Bsel)

   for ntest = 1:30
      Xs = rand(PositionState{Float64}, B1p.bases[1], nX)
      cfg = ACEConfig(Xs)
      BBvec = evaluate(basis1, cfg)
      value1 = [reshape(BBvec[i].val, 2L+1, 1) for i in 1:length(BBvec)]
      BBmat = evaluate(basis2, cfg)
      value2 = [BBmat[i].val for i in 1:length(BBvec)]
      print_tf(@test isapprox(value1, value2, rtol=1e-10))
   end
   println()
end

##
@info("Consistency between Invariant Scalar & SphericalMatrix")

φ = ACE.Invariant()
basis = SymmetricBasis(φ, B1p, Bsel)
φ2 = ACE.SphericalMatrix(0, 0; T = ComplexF64)
basis2 = SymmetricBasis(φ2, B1p, Bsel)

for ntest = 1:30
   local Xs, cfg, BB
   Xs = rand(PositionState{Float64}, B1p.bases[1], nX)
   cfg = ACEConfig(Xs)
   BB = evaluate(basis, cfg)
   BBsca = [BB[i].val for i in 1:length(BB)]
   BB2 =  evaluate(basis2, cfg)
   BBCFlo = [ComplexF64(BB2[i].val...) for i in 1:length(BB2)]
   print_tf(@test isapprox(BBsca, BBCFlo, rtol=1e-10))
end
println()

# ## Keep for futher profiling
#
# L1 = 1; L2 = 1
# φ = ACE.SphericalMatrix(L1, L2; T = ComplexF64)
# pibasis = PIBasis(B1p, ord, maxdeg; property = φ, isreal = false)
# basis = SymmetricBasis(pibasis, φ)
# @time SymmetricBasis(pibasis, φ);
#
# Profile.clear()
# @profile SymmetricBasis(pibasis, φ);
# ProfileView.view()

##


text:
# name = "Bonmin_jll"
# Bonmin_jll = "29cba6d7-6840-5cf2-a2fa-9bdfccfe29ea"
# name = "AmplNLWriter"
# AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"

using AmplNLWriter, Bonmin_jll

BONMIN_OPT = AmplNLWriter.Optimizer(Bonmin_jll.amplexe)
# MOI.set(BONMIN_OPT, MOI.Silent(), true)
# MOI.set(BONMIN_OPT, MOI.TimeLimitSec(), 5)
# BONMIN_OPT = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(BONMIN_OPT)

push!(solvers_nlp, (opt = BONMIN_OPT, mode = BilevelJuMP.ProductMode(1e-5)))


text:
abstract type ResourceKey end

abstract type AbstractResource end

function show(io::IO, res::AbstractResource)
  print(io, "$(typeof(res))")
end

abstract type ResourceEvent <: AbstractEvent end

struct Put <: ResourceEvent
  bev :: BaseEvent
  function Put(env::Environment)
    new(BaseEvent(env))
  end
end

struct Get <: ResourceEvent
  bev :: BaseEvent
  function Get(env::Environment)
    new(BaseEvent(env))
  end
end

function isless(a::ResourceKey, b::ResourceKey)
  (a.priority < b.priority) || (a.priority == b.priority && a.id < b.id)
end

function trigger_put(put_ev::ResourceEvent, res::AbstractResource)
  queue = DataStructures.PriorityQueue(res.put_queue)
  while length(queue) > 0
    (put_ev, key) = DataStructures.peek(queue)
    proceed = do_put(res, put_ev, key)
    state(put_ev) == scheduled && DataStructures.dequeue!(res.put_queue, put_ev)
    proceed ? DataStructures.dequeue!(queue) : break
  end
end

function trigger_get(get_ev::ResourceEvent, res::AbstractResource)
  queue = DataStructures.PriorityQueue(res.get_queue)
  while length(queue) > 0
    (get_ev, key) = DataStructures.peek(queue)
    proceed = do_get(res, get_ev, key)
    state(get_ev) == scheduled && DataStructures.dequeue!(res.get_queue, get_ev)
    proceed ? DataStructures.dequeue!(queue) : break
  end
end

function cancel(res::AbstractResource, put_ev::Put)
  DataStructures.dequeue!(res.put_queue, put_ev)
end

function cancel(res::AbstractResource, get_ev::Get)
  DataStructures.dequeue!(res.get_queue, get_ev)
end


text:
using DataFrames, GLM

df = simulate_data_correlation([[0.2, 0.5], [0.4]])

@test nrow(datagrid(df, n=3)) == 54
@test nrow(datagrid(df, cols=["Var1", "Var2"], n=3)) == 9
@test nrow(datagrid(df, cols=[:Var1, :Var2], n=3)) == 9
@test nrow(datagrid(df, cols=[1, 2], n=3)) == 9
@test_throws ArgumentError datagrid(df, cols=[1.5, 2.2], n=3)
@test nrow(datagrid(df, cols=[1], n=3, fix_num=minimum)) == 3
@test nrow(datagrid(df, cols=[1], n=3, fix_num=3)) == 3
@test nrow(datagrid(df, cols=[1], n=3, fix_num=3, fix_fac=2)) == 3
@test length(datagrid([1, 3, missing])) == 10
@test length(datagrid(["dupa", "a", missing])) == 2

model = GLM.lm(GLM.@formula(y ~ Var1), df)

@test nrow(datagrid(model)) == 10


text:
# This file is a part of RadiationDetectorDSP.jl, licensed under the MIT License (MIT).

__precompile__(true)

module RadiationDetectorDSP

using ArgCheck
using ArraysOfArrays
using Distributions
using DSP
using ElasticArrays
using LinearAlgebra
using ParallelProcessingTools
using Parameters
using RadiationDetectorSignals
using Random
using RecipesBase
using StaticArrays
using Statistics
using StatsBase
using TypedTables

# package code goes here

end # module


text:
using Random
using WindingNumbers
using WindingNumbers: angle, Point
Random.seed!(10)

using Test

@testset "Winding Numbers tests" begin

@testset "Basic angle tests" begin
  a = Point(1.0, 0.0)
  b = Point(0.0, 1.0)
  @test π/2 ≈ angle(b) - angle(a)

  a = Point(-1.0, 0.0)
  b = Point(1.0, 0.0)
  @test -π ≈ angle(b) - angle(a)

  a = Point(1.0, 0.0)
  b = Point(-1.0, 0.0)
  @test π ≈ angle(b) - angle(a)

  a = Point(1.0, 0.0)
  b = Point(1.0, 1.0)
  @test π/4 ≈ angle(b) - angle(a)
 end

function mock(x::Vector{Float64}, n, r=zeros(2))
  θ = atan(x[2] - r[2], x[1] - r[1])
  return [cos(n*θ), sin(n*θ)]
end

turnings(n::Int) = -n:n

@testset "Basic understanding - angles" begin
  a = Point(1.0, 1.0)
  b = Point(-1.0, 1.0)
  c = Point(-1.0, -1.0)
  d = Point(1.0, -1.0)
  @test angle(a, b) ≈ π/2
  @test angle(b, c) ≈ π/2
  @test angle(c, d) ≈ π/2
  @test angle(d, a) ≈ π/2
  @test angle(b, a) ≈ -π/2
  @test angle(c, b) ≈ -π/2
  @test angle(d, c) ≈ -π/2
  @test angle(a, d) ≈ -π/2
  ϕ = angle(Point(-1.0, 2*eps()), Point(-1.0, -2*eps()))
  @test ϕ ≈ 4*eps()
  @test ϕ > 0.0
  ϕ = angle(Point(-1.0, -2*eps()), Point(-1.0, 2*eps()))
  @test ϕ ≈ -4*eps()
  @test ϕ < 0.0
end

@testset "Basic understanding - square with pole" begin
  a = Point(1.0, 1.0)
  b = Point(-1.0, 1.0)
  c = Point(-1.0, -1.0)
  d = Point(1.0, -1.0)
  ϕ = angle(a, b) + angle(b, c) + angle(c, d) + angle(d, a)
  @test ϕ ≈ 2π
  ϕ = angle(b, a) + angle(c, b) + angle(d, c) + angle(a, d)
  @test ϕ ≈ -2π
end

@testset "Basic understanding - square without pole" begin
  a = Point(1.0, 1.0)
  b = Point(2.0, 1.0)
  c = Point(2.0, 2.0)
  d = Point(1.0, 2.0)
  ϕ = angle(a, b) + angle(b, c) + angle(c, d) + angle(d, a)
  @test ϕ ≈ 0.0
  ϕ = angle(b, a) + angle(c, b) + angle(d, c) + angle(a, d)
  @test ϕ ≈ 0.0
end

@testset "Pole at origin in grid" begin
  for i in turnings(5)
    n = WindingNumbers.windingnumber(x->mock(x, i), -ones(2), ones(2))
    @test Int(round(n)) == i
  end
end

@testset "Random location of pole inside grid" begin
  for i in turnings(2)
    r = rand(2) # solution
    a, b = r - rand(2), r + rand(2)
    @assert all(a .< r .<b)
    n = WindingNumbers.windingnumber(x->mock(x, i, r), a, b)
    @test Int(round(n)) == i
  end
end

@testset "Unphysical pole" begin
  for i in turnings(5)
    s = rand(2) # square
    @assert s[1] > 0.0
    @assert s[2] > 0.0
    n = WindingNumbers.windingnumber(x->[0.0, 0.0], s, 2*s)
    @test Int(round(n)) == 0
  end
end

@testset "Physical pole not in grid" begin
  for i in turnings(5)
    s = rand(2) # square
    @assert s[1] > 0.0
    @assert s[2] > 0.0
    n = WindingNumbers.windingnumber(x->mock(x, i), s, 2*s)
    @test Int(round(n)) == 0
  end
end

@testset "solve" begin
  for i ∈ 1:100
    rtol = 10^Float64(-rand(1:14))
    root = 2 .* rand(2) .- 1
    lower = root .- rand(2)
    upper = root .+ rand(2)
    allpointsinside = true
    function solvemock(x)
      allpointsinside &= all(lower .<= x .<= upper)
      y = mock(x, 1, root)
      return y[1] + im * y[2]
    end
    solutions = WindingNumbers.boundingboxes(solvemock, lower, upper,
                                             xtol_rel=rtol, stopval=eps())
    @test length(solutions) == 1
    @test allpointsinside
    length(solutions) == 1 || continue
    A, B = solutions[1]
    @test all(A .<= root .<= B)
    @test all((B .- A) ./ (upper .- lower) .< rtol)
  end
end

end

max_stars_repo_name: jwscook/WindingNumbers

max_issues_repo_name: jwscook/WindingNumbers

max_forks_repo_name: jwscook/WindingNumbers


text:
"""
    TypeN(var_1::TypeM)

Lorem ipsum dolor sit amet.
"""
struct TypeN

  var_1::TypeM

end


text:
include("../tcm_code/tcm_estimateForest_L1.jl")

function tcm_estimateForest_L1_RTS(N, assortments, choice_probs, depth, K)

	M = size(assortments,1)

	forest_left = Array{Int64,1}[]
    forest_right = Array{Int64,1}[]
    forest_product = Array{Int64,1}[]
    forest_isLeaf = Array{Bool,1}[]

	v_by_class = Array{Int64,2}[];

    start_time = time()

    for k in 1:K
        tree_left, tree_right, tree_product, tree_isLeaf = tcm_randomCompleteTree( N, depth-1 )

        v_one_class = tcm_treeToA( tree_left, tree_right, tree_product, tree_isLeaf, assortments)
        push!(v_by_class, v_one_class)

        push!(forest_left, tree_left)
        push!(forest_right, tree_right)
        push!(forest_product, tree_product)
        push!(forest_isLeaf, tree_isLeaf)
    end

	m_primal = Model(solver = GurobiSolver())
	@variable(m_primal, lambda_var[1:K] >= 0)
	@variable(m_primal, eps_plus[1:M,1:N] >= 0)
	@variable(m_primal, eps_minus[1:M,1:N] >= 0)

	for m in 1:M
		for i in 1:N
			@constraint(m_primal, sum( v_by_class[k][m,i] * lambda_var[k] for k in 1:K) + eps_plus[m,i] - eps_minus[m,i] == choice_probs[m,i])
		end
	end
	@constraint(m_primal, sum(lambda_var[k] for k in 1:K) == 1)

	@objective(m_primal, Min, sum( eps_minus[m,i] for m in 1:M, i in 1:N) + sum(eps_plus[m,i] for m in 1:M, i in 1:N))

	solve(m_primal)

    end_time = time()
    elapsed_time = end_time - start_time

	lambda_val = getvalue(lambda_var)

	primal_objval = getobjectivevalue(m_primal)

	return lambda_val, forest_left, forest_right, forest_product, forest_isLeaf, primal_objval, elapsed_time

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
#Implementation of the EM solver.
struct EMproto <: Solver end #prototyping code to check for correctness
struct EM <: Solver end

function integrate(a!::Function, b!::Function, x0::AbstractVector{T}, t::AbstractVector{<:Real},
    noise::Noise, em::EMproto, callback::Function) where {T<:Number}
    W = W0(noise)
    x = copy(x0)

    n = length(x0)
    f = Vector{T}(n)
    g = Matrix{T}(n, noise.m)
    ΔW = copy(W)

    #save the first data point
    callback(t[1], x, W)

    for i in 1:length(t)-1
        Δt = t[i+1]-t[i]
        randn!(noise, ΔW)
        ΔW = ΔW*sqrt(Δt)

        a!(t[i], x, f)
        b!(t[i], x, g)

        x .= x + f*Δt + g*ΔW
        W .= W + ΔW
        callback(t, x, W)
    end
    return nothing
end

"""
explicitEM function that is optimized
ab! : bias & diffusion function
x0 : Initial state
t : times to use for time evolution
noise : Noise object
callback : a function that runs on t, x, W
"""
function integrate(a!::Function, b!::Function, x0::AbstractVector{T}, t::AbstractVector{<:Real},
    noise::Noise, em::EM, callback::Function) where {T<:Number}
    W = W0(noise)
    x = copy(x0)

    n = length(x0)
    fΔt = Vector{T}(n)
    g = Matrix{T}(n, noise.m)
    ΔW = copy(W)
    gΔW = Vector{T}(n)

    #save the first data point
    callback(t[1], x, W)

    for i in 1:length(t)-1
        Δt = t[i+1]-t[i]
        randn!(noise, ΔW)
        ΔW .*= sqrt(Δt)

        a!(t[i], x, fΔt)
        b!(t[i], x, g)

        fΔt .*= Δt
        A_mul_B!(gΔW, g, ΔW) #This line performs better than the more readable version!

        x .+= fΔt
        x .+= gΔW
        W .+= ΔW
        callback(t, x, W)
    end
    return nothing
end

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
#------------------------------------------------------------------------------

function Base.convert(::Type{Int}, x::gfp_elem)
    return Int(x.data)
end

function Base.convert(::Type{BigInt}, x::gfp_fmpz_elem)
    return BigInt(x.data)
end

function Base.convert(::Type{FlintIntegerRing}, x::gfp_fmpz_elem)
    return ZZ(x.data)
end

function Base.convert(::Type{FlintIntegerRing}, x::gfp_elem)
    return ZZ(x.data)
end

# Bad
# Base.zero(::Type{gfp_elem}) = GF(2)(2)

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

# extends random generator to Rational field
Base.rand(::FlintRationalField) = Nemo.QQ(rand(-2^16:2^16))
Base.rand(::FlintRationalField, n::Int) = [rand(Nemo.QQ) for _ in 1:n]

Base.rand(::FlintIntegerRing) = Nemo.ZZ(rand(-2^16:2^16))
Base.rand(::FlintIntegerRing, n::Int) = [rand(Nemo.ZZ) for _ in 1:n]

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

function normtime(x)
     x/1e9
end

# estimates the elapsed time for `ex` and stores the result into the `storage`
macro savetime(ex, storage)
    return quote
        local t0 = time_ns()
        local val = $(esc(ex))
        local t1 = time_ns()
        push!($storage, normtime(t1-t0))
        val
    end
end

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

# conveniece wrapper for `from_dense` function,
# which returns a sparse representation of an array
# and generally defaults to
#   1D vector --> Sparsik
#   2D vector --> DOK_Sparsik
# over the Rational field of Nemo.QQ instance
"""
    @sparse(ex)

Converts an `Array` in `ex` to sparse representation
"""
macro sparse(ex)
    return quote
        from_dense($(esc(ex)), Nemo.QQ)
    end
end

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

# make the given path absolute w.r.t the project and normalize it
function normalizepath(path)
    replace("$(normpath(joinpath(@__FILE__, "..", "..")))$path", "\\" => "/")
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: schemas
using AWS.Compat
using AWS.UUIDs
"""
    CreateDiscoverer()

Creates a discoverer.

# Required Parameters
- `SourceArn`: The ARN of the event bus.

# Optional Parameters
- `Description`: A description for the discoverer.
- `tags`: Tags associated with the resource.
"""

create_discoverer(SourceArn; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discoverers", Dict{String, Any}("SourceArn"=>SourceArn); aws_config=aws_config)
create_discoverer(SourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discoverers", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("SourceArn"=>SourceArn), args)); aws_config=aws_config)

"""
    CreateRegistry()

Creates a registry.

# Required Parameters
- `registryName`: The name of the registry.

# Optional Parameters
- `Description`: A description of the registry to be created.
- `tags`: Tags to associate with the registry.
"""

create_registry(registryName; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/registries/name/$(registryName)"; aws_config=aws_config)
create_registry(registryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/registries/name/$(registryName)", args; aws_config=aws_config)

"""
    CreateSchema()

Creates a schema definition. Inactive schemas will be deleted after two years.

# Required Parameters
- `Content`: The source of the schema definition.
- `Type`: The type of schema.
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

# Optional Parameters
- `Description`: A description of the schema.
- `tags`: Tags associated with the schema.
"""

create_schema(Content, Type, registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)", Dict{String, Any}("Content"=>Content, "Type"=>Type); aws_config=aws_config)
create_schema(Content, Type, registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Content"=>Content, "Type"=>Type), args)); aws_config=aws_config)

"""
    DeleteDiscoverer()

Deletes a discoverer.

# Required Parameters
- `discovererId`: The ID of the discoverer.

"""

delete_discoverer(discovererId; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/discoverers/id/$(discovererId)"; aws_config=aws_config)
delete_discoverer(discovererId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/discoverers/id/$(discovererId)", args; aws_config=aws_config)

"""
    DeleteRegistry()

Deletes a Registry.

# Required Parameters
- `registryName`: The name of the registry.

"""

delete_registry(registryName; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/registries/name/$(registryName)"; aws_config=aws_config)
delete_registry(registryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/registries/name/$(registryName)", args; aws_config=aws_config)

"""
    DeleteResourcePolicy()

Delete the resource-based policy attached to the specified registry.

# Optional Parameters
- `registryName`: The name of the registry.
"""

delete_resource_policy(; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/policy"; aws_config=aws_config)
delete_resource_policy(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/policy", args; aws_config=aws_config)

"""
    DeleteSchema()

Delete a schema definition.

# Required Parameters
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

"""

delete_schema(registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)"; aws_config=aws_config)
delete_schema(registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)", args; aws_config=aws_config)

"""
    DeleteSchemaVersion()

Delete the schema version definition

# Required Parameters
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.
- `schemaVersion`: The version number of the schema

"""

delete_schema_version(registryName, schemaName, schemaVersion; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/version/$(schemaVersion)"; aws_config=aws_config)
delete_schema_version(registryName, schemaName, schemaVersion, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/version/$(schemaVersion)", args; aws_config=aws_config)

"""
    DescribeCodeBinding()

Describe the code binding URI.

# Required Parameters
- `language`: The language of the code binding.
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

# Optional Parameters
- `schemaVersion`: Specifying this limits the results to only this schema version.
"""

describe_code_binding(language, registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/language/$(language)"; aws_config=aws_config)
describe_code_binding(language, registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/language/$(language)", args; aws_config=aws_config)

"""
    DescribeDiscoverer()

Describes the discoverer.

# Required Parameters
- `discovererId`: The ID of the discoverer.

"""

describe_discoverer(discovererId; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/discoverers/id/$(discovererId)"; aws_config=aws_config)
describe_discoverer(discovererId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/discoverers/id/$(discovererId)", args; aws_config=aws_config)

"""
    DescribeRegistry()

Describes the registry.

# Required Parameters
- `registryName`: The name of the registry.

"""

describe_registry(registryName; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)"; aws_config=aws_config)
describe_registry(registryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)", args; aws_config=aws_config)

"""
    DescribeSchema()

Retrieve the schema definition.

# Required Parameters
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

# Optional Parameters
- `schemaVersion`: Specifying this limits the results to only this schema version.
"""

describe_schema(registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)"; aws_config=aws_config)
describe_schema(registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)", args; aws_config=aws_config)

"""
    ExportSchema()

# Required Parameters
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.
- `type`:

# Optional Parameters
- `schemaVersion`: Specifying this limits the results to only this schema version.
"""

export_schema(registryName, schemaName, type; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/export", Dict{String, Any}("type"=>type); aws_config=aws_config)
export_schema(registryName, schemaName, type, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/export", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("type"=>type), args)); aws_config=aws_config)

"""
    GetCodeBindingSource()

Get the code binding source URI.

# Required Parameters
- `language`: The language of the code binding.
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

# Optional Parameters
- `schemaVersion`: Specifying this limits the results to only this schema version.
"""

get_code_binding_source(language, registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/language/$(language)/source"; aws_config=aws_config)
get_code_binding_source(language, registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/language/$(language)/source", args; aws_config=aws_config)

"""
    GetDiscoveredSchema()

Get the discovered schema that was generated based on sampled events.

# Required Parameters
- `Events`: An array of strings where each string is a JSON event. These are the events that were used to generate the schema. The array includes a single type of event and has a maximum size of 10 events.
- `Type`: The type of event.

"""

get_discovered_schema(Events, Type; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discover", Dict{String, Any}("Events"=>Events, "Type"=>Type); aws_config=aws_config)
get_discovered_schema(Events, Type, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discover", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Events"=>Events, "Type"=>Type), args)); aws_config=aws_config)

"""
    GetResourcePolicy()

Retrieves the resource-based policy attached to a given registry.

# Optional Parameters
- `registryName`: The name of the registry.
"""

get_resource_policy(; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/policy"; aws_config=aws_config)
get_resource_policy(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/policy", args; aws_config=aws_config)

"""
    ListDiscoverers()

List the discoverers.

# Optional Parameters
- `discovererIdPrefix`: Specifying this limits the results to only those discoverer IDs that start with the specified prefix.
- `limit`:
- `nextToken`: The token that specifies the next page of results to return. To request the first page, leave NextToken empty. The token will expire in 24 hours, and cannot be shared with other accounts.
- `sourceArnPrefix`: Specifying this limits the results to only those ARNs that start with the specified prefix.
"""

list_discoverers(; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/discoverers"; aws_config=aws_config)
list_discoverers(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/discoverers", args; aws_config=aws_config)

"""
    ListRegistries()

List the registries.

# Optional Parameters
- `limit`:
- `nextToken`: The token that specifies the next page of results to return. To request the first page, leave NextToken empty. The token will expire in 24 hours, and cannot be shared with other accounts.
- `registryNamePrefix`: Specifying this limits the results to only those registry names that start with the specified prefix.
- `scope`: Can be set to Local or AWS to limit responses to your custom registries, or the ones provided by AWS.
"""

list_registries(; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries"; aws_config=aws_config)
list_registries(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries", args; aws_config=aws_config)

"""
    ListSchemaVersions()

Provides a list of the schema versions and related information.

# Required Parameters
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

# Optional Parameters
- `limit`:
- `nextToken`: The token that specifies the next page of results to return. To request the first page, leave NextToken empty. The token will expire in 24 hours, and cannot be shared with other accounts.
"""

list_schema_versions(registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/versions"; aws_config=aws_config)
list_schema_versions(registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/versions", args; aws_config=aws_config)

"""
    ListSchemas()

List the schemas.

# Required Parameters
- `registryName`: The name of the registry.

# Optional Parameters
- `limit`:
- `nextToken`: The token that specifies the next page of results to return. To request the first page, leave NextToken empty. The token will expire in 24 hours, and cannot be shared with other accounts.
- `schemaNamePrefix`: Specifying this limits the results to only those schema names that start with the specified prefix.
"""

list_schemas(registryName; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas"; aws_config=aws_config)
list_schemas(registryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas", args; aws_config=aws_config)

"""
    ListTagsForResource()

Get tags for resource.

# Required Parameters
- `resource-arn`: The ARN of the resource.

"""

list_tags_for_resource(resource_arn; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/tags/$(resource-arn)", args; aws_config=aws_config)

"""
    PutCodeBinding()

Put code binding URI

# Required Parameters
- `language`: The language of the code binding.
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

# Optional Parameters
- `schemaVersion`: Specifying this limits the results to only this schema version.
"""

put_code_binding(language, registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/language/$(language)"; aws_config=aws_config)
put_code_binding(language, registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)/language/$(language)", args; aws_config=aws_config)

"""
    PutResourcePolicy()

The name of the policy.

# Required Parameters
- `Policy`: The resource-based policy.

# Optional Parameters
- `RevisionId`: The revision ID of the policy.
- `registryName`: The name of the registry.
"""

put_resource_policy(Policy; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/policy", Dict{String, Any}("Policy"=>Policy); aws_config=aws_config)
put_resource_policy(Policy, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/policy", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Policy"=>Policy), args)); aws_config=aws_config)

"""
    SearchSchemas()

Search the schemas

# Required Parameters
- `keywords`: Specifying this limits the results to only schemas that include the provided keywords.
- `registryName`: The name of the registry.

# Optional Parameters
- `limit`:
- `nextToken`: The token that specifies the next page of results to return. To request the first page, leave NextToken empty. The token will expire in 24 hours, and cannot be shared with other accounts.
"""

search_schemas(keywords, registryName; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/search", Dict{String, Any}("keywords"=>keywords); aws_config=aws_config)
search_schemas(keywords, registryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("GET", "/v1/registries/name/$(registryName)/schemas/search", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("keywords"=>keywords), args)); aws_config=aws_config)

"""
    StartDiscoverer()

Starts the discoverer

# Required Parameters
- `discovererId`: The ID of the discoverer.

"""

start_discoverer(discovererId; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discoverers/id/$(discovererId)/start"; aws_config=aws_config)
start_discoverer(discovererId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discoverers/id/$(discovererId)/start", args; aws_config=aws_config)

"""
    StopDiscoverer()

Stops the discoverer

# Required Parameters
- `discovererId`: The ID of the discoverer.

"""

stop_discoverer(discovererId; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discoverers/id/$(discovererId)/stop"; aws_config=aws_config)
stop_discoverer(discovererId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/v1/discoverers/id/$(discovererId)/stop", args; aws_config=aws_config)

"""
    TagResource()

Add tags to a resource.

# Required Parameters
- `resource-arn`: The ARN of the resource.
- `tags`: Tags associated with the resource.

"""

tag_resource(resource_arn, tags; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes tags from a resource.

# Required Parameters
- `resource-arn`: The ARN of the resource.
- `tagKeys`: Keys of key-value pairs.

"""

untag_resource(resource_arn, tagKeys; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateDiscoverer()

Updates the discoverer

# Required Parameters
- `discovererId`: The ID of the discoverer.

# Optional Parameters
- `Description`: The description of the discoverer to update.
"""

update_discoverer(discovererId; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/discoverers/id/$(discovererId)"; aws_config=aws_config)
update_discoverer(discovererId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/discoverers/id/$(discovererId)", args; aws_config=aws_config)

"""
    UpdateRegistry()

Updates a registry.

# Required Parameters
- `registryName`: The name of the registry.

# Optional Parameters
- `Description`: The description of the registry to update.
"""

update_registry(registryName; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/registries/name/$(registryName)"; aws_config=aws_config)
update_registry(registryName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/registries/name/$(registryName)", args; aws_config=aws_config)

"""
    UpdateSchema()

Updates the schema definition Inactive schemas will be deleted after two years.

# Required Parameters
- `registryName`: The name of the registry.
- `schemaName`: The name of the schema.

# Optional Parameters
- `ClientTokenId`: The ID of the client token.
- `Content`: The source of the schema definition.
- `Description`: The description of the schema.
- `Type`: The schema type for the events schema.
"""

update_schema(registryName, schemaName; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)", Dict{String, Any}("ClientTokenId"=>string(uuid4())); aws_config=aws_config)
update_schema(registryName, schemaName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = schemas("PUT", "/v1/registries/name/$(registryName)/schemas/name/$(schemaName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientTokenId"=>string(uuid4())), args)); aws_config=aws_config)


text:
using HorizonSideRobots

function movements!(r::Robot, side::HorizonSide, n::Int)
    while (n > 0 && !ismarker(r))
        move!(r, side)
        n -= 1
    end
end

function find_marker(r::Robot)
    i = 3
    keep_const = 0
    counter = 0
    while !ismarker(r)
        movements!(r, HorizonSide(i % 4), counter + keep_const)
        i += 1
        keep_const = 1 - 1 * keep_const
        counter += keep_const
    end
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
# return the subset of expressed markers
function filterdf(df,colnames)
    filtnames = []
    for cn in colnames
        ff = findfirst(n->n==cn,names(df))
        if !isnothing(ff)
            push!(filtnames,cn)
        end
    end
    unique(filtnames)
end

# get the average over a subset of markers
function filteravg(df,colnames)
    cn = filterdf(df,colnames)
    mean(eachcol(df[!,Symbol.(cn)]))
end

# get the sum over a subset of markers
function filtersum(df,colnames)
    cn = filterdf(df,colnames)
    sum(eachcol(df[!,Symbol.(cn)]))
end

# nms is feature columns for umap, all others interpreted as labels during stack
function umapexpr(df,nms,params)
    inds = [all(map(n->n!=nm,nms)) for nm in names(df)]
    um = umap(Matrix(df[!,inds])',
                n_neighbors=params[:nneigh],
                min_dist=params[:mindist],
                n_epochs=params[:nepochs])

    # scatter(res_jl[1,:], res_jl[2,:], zcolor=mnist_y,
    #         title="MNIST: Julia UMAP", marker=(2, 2, :auto, stroke(0)))
    #         inds = [all(map(n->n!=nm,nms)) for nm in names(df)]

end

# feature plot with dimensionality reduction
# df is a dataframe with shape samples x features
# ids are identifiers, such that names(df) \ ids  = the features we care about
# trans is a transformation of the samples x features passed to the dimred call
# dimred is a function that outputs a 2d feature vector
# the feature vector is added to df as feat_x and feat_y
function fplot(df,ids,dimred,trans,args)
    # NOT IMPLEMENTED
    dat = df[!,[all(map(n->n!=nm,ids)) for nm in names(df)]]
    reduced = eval(Expr(:call,dimred,trans(dat),args))
end

# given a list of at least 2 vectors, iteratively replace nonmissing values with right associativity
function sublabels(labels...)
    @assert length(labels) >= 2 "less than 2 label vectors provided"
    lab = labels[1]
    for l in 2:length(labels)
        for ll in 1:length(lab)
            if !ismissing(labels[l][ll])
                local lab[ll] = labels[l][ll]
            end
        end
    end
    lab
end


text:
dat = readdlm("mds_G7-1950.csv",';')
cnames = dat[1,:]
dat    = dat[2:end,:]

nrow,ncol = size(dat)
pdat      = hcat([log.(dat[:,i]) - log.(dat[:,j]) for j=1:7, i=1:7 if i<j]...)
pnames    = [cnames[i] * " - " * cnames[j] for j=1:7, i=1:7 if i<j]

w = pdat[:,1]
b = [.4 ,.8,   1.1,    .3,    .3,.3,.3,.3,.3,.4,.1,.2,.3]
function lnviD3(b,w)

  # b      = [.4 ,.8,   1.1,    .3,    .3,.3,.3,.3,.3,.4,.1,.2,.3]
  # b = inits
  # w = randn(1000)
  datsize = length(w)

  d1hat, d2hat, d3hat  = b[1:3];
  p11hat, p22hat, p33hat = b[4:6];
  p12hat, p21hat, p31hat = b[7:9];
  sigma1hat, sigma2hat, sigma3hat =  b[10], b[10], b[10];
  mu1hat, mu2hat, mu3hat = b[11:13];

  rho1hat=0; theta1hat=0;
  # tranmax = Array(Float64,3,3);
  # [tranmax[i,i] = b[i+3]  for i = 1:3]
  tranmax = diagm(b[4:6])

  tranmax[[4 2 3]] = [p12hat,p21hat,p31hat]
  tranmax[[7 8 6]] = 1-sum(tranmax, 2)

  mu     = [mu1hat,mu2hat,mu3hat]
  sigma  = [sigma1hat,sigma2hat,sigma3hat]
  coefd1 = cumprod((d1hat+(0:(datsize-2)))./(1:(datsize-1)))
  coefd2 = cumprod((d2hat+(0:(datsize-2)))./(1:(datsize-1)))
  coefd3 = cumprod((d3hat+(0:(datsize-2)))./(1:(datsize-1)))
  coeff  = [coefd1 coefd2 coefd3]

  # muhat is to store the selected path of state: from date 2 */
  eta1  = zeros(9,datsize); # row 1:3 for path1, 4:6 for path2 */
  eta   = zeros(9,datsize); # same as above */
  xihat = zeros(9,datsize); # same as above */
  ehat  = zeros(3,datsize); # row 1, for path1, row 2 for path2 */

  # to keep the long memory part, not iid */
  zsigmahat = zeros(3,datsize) # row 1, for path1, row 2 for path2 */
  zhat      = zeros(3,datsize) # row 1, for path1, row 2 for path2 */
  lnlk      = zeros(3,datsize)

  path1 = zeros(Int8,3,datsize)
  path2 = zeros(Int8,3,datsize)
  path3 = zeros(Int8,3,datsize)
  path1[:,1] = [1,0,0]
  path2[:,1] = [0,1,0]
  path3[:,1] = [0,0,1]

  zsigmahat[:,1] = w[1] - mu
  # zsigmahat is de-drifted series for each given path; or it is z*sigmahat

  zhat[:,1] = zsigmahat[:,1] ./ sigma;

  ehat[:,1] = zhat[:,1];

  limprobs  = (tranmax^7)[1,:]
  lnlk[:,1] = log.(exp.(-ehat[:,1].^2./2 + 1e-50 ) ./ (sqrt(2*pi)).*limprobs);

  j = 2;

  ####################!!!!!!!!!!!!!!!!!!!!!
  while j <= datsize
    # tempsMUZ and temps are needed for likelihood possibility

    tempsMUZ = vec(transpose([(zsigmahat[1,1:(j-1)])[end:-1:1] (zsigmahat[2,1:(j-1)])[end:-1:1] (zsigmahat[3,1:(j-1)])[end:-1:1]]) * coeff[1:(j-1),:])

    temps    = repeat(mu; inner=[3]) + rho1hat * repeat(sigma;inner=[3]) .* repeat(zhat[:,j-1];outer=[3]) +
      theta1hat*repeat(sigma,inner=[3]) .* repeat(ehat[:,j-1], outer=[3])

    # likelihood possibility
    eta1[:,j]  = exp.(-(w[j] - temps - tempsMUZ).^2 ./ (2*repeat(sigma;inner=[3]).^2+1e-50)) ./
      (sqrt.(2*pi*repeat(sigma,inner=[3]).^2))

    eta[:,j] = eta1[:,j] .* vec(tranmax)

    # Case is 1 for 1->1; 2 for 2->1; 3 for 1->2; 4 for 2->2
    case  = reshape([log(eta[x,j]+1e-50) + lnlk[(x-1) % 3 + 1,j-1] for x in 1:9],3,3)
    case  = findmax(case,1)[2]
    case2 = (case-1) % 3 + 1

    # holding previous paths in temp memory
    oldlnlk = lnlk
    immutable oldPath2
      path
    end

    oldzsigmahat = zsigmahat
    oldzhat = zhat
    oldehat = ehat

    # recalculating values
    # path ends up with state 1:
    lnlk[1,j]    = log(eta[case[1],j]+1e-50) + oldlnlk[case2[1],j-1]
    path1[:,1:j] = [oldPath[((case2[1]-1)*3+1):(case2[1]*3),1:(j-1)] [1,0,0]]
    zsigmahat[1,1:j] = append!(oldzsigmahat[case2[1],1:(j-1)], w[j] - mu[1] - tempsMUZ[case[1]])
    zhat[1,1:j]  = append!(oldzhat[case2[1],1:(j-1)],zsigmahat[1,j]/(sigma1hat));
    ehat[1,1:j]  = append!(oldehat[case2[1],1:(j-1)],(w[j] - tempsMUZ[case[1]] -temps[case[1]])/sigma1hat);

    # path ends up with state 2:
    lnlk[2,j]    = log(eta[case[2],j]+1e-50) + oldlnlk[case2[2],j-1]
    path2[:,1:j] = [oldPath[((case2[2]-1)*3+1):(case2[2]*3),1:(j-1)] [0,1,0]];
    zsigmahat[2,1:j] = append!(oldzsigmahat[case2[2],1:(j-1)], w[j] - mu[2] - tempsMUZ[case[2]])
    zhat[2,1:j]  = append!(oldzhat[case2[2],1:(j-1)],zsigmahat[2,j]/(sigma2hat));
    ehat[2,1:j]  = append!(oldehat[case2[2],1:(j-1)],(w[j] - tempsMUZ[case[2]] -temps[case[2]])/sigma2hat)

    # path ends up with state 3:
    lnlk[3,j]    = log(eta[case[3],j]+1e-50) + oldlnlk[case2[3],j-1]
    path3[:,1:j] = [oldPath[((case2[3]-1)*3+1):(case2[3]*3),1:(j-1)] [0,0,1]]
    zsigmahat[3,1:j] = append!(oldzsigmahat[case2[3],1:(j-1)], w[j] - mu[3] - tempsMUZ[case[3]])
    zhat[3,1:j]  = append!(oldzhat[case2[3],1:(j-1)],zsigmahat[3,j]/(sigma3hat));
    ehat[3,1:j]  = append!(oldehat[case2[3],1:(j-1)],(w[j] - tempsMUZ[case[3]] -temps[case[3]])/sigma3hat)

    # for(x in 0:2){
    #   lnlk[x+1,j]   = log(eta[case[x+1],j]+1e-50) + oldlnlk[case2[x+1],j-1]
    #   paths[(3*x+1):(3*x+3),1:j] = cbind(oldPath[(3*x+1):(3*x+3),1:(j-1)],(1:3 == (x+1)) * 1)
    #   zsigmahat[x+1,1:j] = c(oldzsigmahat[case2[x+1],1:(j-1)],  w[j] - mu[case2[x+1]] - tempsMUZ[case[x+1]])
    #   zhat[x+1,1:j] = c(oldzhat[case2[x+1],1:(j-1)], zsigmahat[x+1,j]/(sigma[x+1]))
    #   ehat[x+1,1:j] = c(oldehat[case2[x+1],1:(j-1)], (w[j] - tempsMUZ[case[x+1]] - temps[case[x+1]])/sigma[x+1]) # ONEMLI!!!
    # }

    j=j+1

  end
  #   lnlipath1 = t(lnlk1);
  #   lnlipath2 = t(lnlk2);
  lnlkRes = findmax(lnlk[:,datsize])[1][1]
  # if(lnlk[1,datsize] > lnlk[2,datsize]){
  #   lnlkRes = lnlk[1,datsize]} else {lnlkRes = lnlk[2,datsize]}
  #
  return lnlkRes
end

dat = readdlm("mds_G7-1950.csv",';')
cnames = dat[1,:]
dat    = dat[2:end,:]

nrow,ncol = size(dat)
pdat      = hcat([log.(dat[:,i]) - log.(dat[:,j]) for j=1:7, i=1:7 if i<j]...)
pnames    = [cnames[i] * " - " * cnames[j] for j=1:7, i=1:7 if i<j]

source("~/Documents/threestate/arfimaSimH.R")

genseries2 = function(n){
  n1 = ceil(n/3)
  n2 = ceil(n/3)
  n3 = n - n1 - n2
  ser1 = arfima.simH(n1, model = list(dfrac = .9))
  ser2 = arfima.simH(n2, model = list(dfrac = .6)) + tail(ser1,1)
  ser3 = arfima.simH(n3, model = list(dfrac = .2)) + tail(ser2,1)
  c(ser1,ser2,ser3)
}

# inits      = c(.4 ,.8,1.1,.8,.8,.8,.1,.1,.1,.01,.1)
# inits      = c(runif(3,-2,2),runif(3,.3,.6))
# inits      = c(inits,.8-inits[4:6],.01,.1)
# lnviD3(inits, ww)
#
# ww = log(read.csv("../Downloads/g7_1950.csv",sep=";")[:,-1])
# ww = ww[:,1] - ww[:,2]
# lowerV = c(-2,-2,-2,.7,.7,.7,0,0,0,0.001,-5,-5,-5)
# upperV = c(2,2,2,.999,.999,.999,.3,.3,.3,5,5,5,5)
# inits  = c(.4 ,   .8,   1.1,    .3,    .3,.3,.3,.3,.3,.01,.1)
# optim(inits, function(p) -lnviD3(p,ww),method="L-BFGS-B",control = list(trace=3))
#
# aa =  c(2,-2,2,0.201,0.201,0.201,0.999,0.999,0.201,0.001,-5 )
# lnviD3(aa, ww)

ser = genseries2(200)
plot(ser)

lowerV = c(-2,-2,-2, .6 , .6 , .6 , 0, 0, 0,0.001,-5,-5,-5)
inits  = c( 1,.5, 0, .61, .61, .61,.1,.1,.1,   .4,.1,.1,.1)
upperV = c( 2, 2, 2,.999,.999,.999,.3,.3,.3,    5, 5, 5, 5)

const_mat = matrix(0,13,13)
diag(const_mat) = 1
const_mat = rbind(const_mat,-const_mat)
dmmvec = -c(0,0,0,1,0,0,1,0,0,0,0,0,0)
dmmvec = c(dmmvec,0,dmmvec,0,dmmvec)[1:39]
const_mat = rbind(const_mat,matrix(dmmvec,3,,T))
const_mat = cbind(const_mat,c(lowerV,-upperV,rep(-1,3)))

res = constrOptim(inits, function(p) -lnviD3(p,ser), NULL, ui = const_mat[:,-14], const_mat[:,14])
thepath = return_path(res$par,ser)
plot(apply(res$par[1:3]*thepath,2,sum),type="l")

max_stars_repo_name: fuadcan/MSReplication

max_issues_repo_name: fuadcan/MSReplication

max_forks_repo_name: fuadcan/MSReplication


text:
export Message
"""
A [`Message`](@ref) activity.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-activity-structure).
"""
struct MessageActivity
    type::Int
    party_id::Optional{String}
end
@boilerplate MessageActivity :constructors :docs :lower :merge :mock

"""
A Rich Presence [`Message`](@ref)'s application information.
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object-message-application-structure).
"""
struct MessageApplication <: DiscordObject
    id::Snowflake
    cover_image::Optional{String}
    description::String
    icon::String
    name::String
end
@boilerplate MessageApplication :constructors :docs :lower :merge :mock
struct MessageReference <: DiscordObject
    message_id::Snowflake
    channel_id::Optional{Snowflake}
    guild_id::Optional{Snowflake}
    fail_if_not_exists::Optional{Bool}
end
@boilerplate MessageReference :constructors :lower :merge

"""
A message sent to a [`DiscordChannel`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/channel#message-object).
"""
struct Message <: DiscordObject
    id::Snowflake
    channel_id::Snowflake
    # MessageUpdate only requires the ID and channel ID.
    guild_id::Optional{Snowflake}
    author::Optional{User}
    member::Optional{Member}
    content::Optional{String}
    timestamp::Optional{DateTime}
    edited_timestamp::OptionalNullable{DateTime}
    tts::Optional{Bool}
    mention_everyone::Optional{Bool}
    mentions::Optional{Vector{User}}
    mention_roles::Optional{Vector{Snowflake}}
    attachments::Optional{Vector{Attachment}}
    message_reference::Optional{MessageReference}
    embeds::Optional{Vector{Embed}}
    reactions::Optional{Vector{Reaction}}
    nonce::OptionalNullable{Snowflake}
    pinned::Optional{Bool}
    webhook_id::Optional{Snowflake}
    type::Optional{Int}
    activity::Optional{MessageActivity}
    application::Optional{MessageApplication}
end
@boilerplate Message :constructors :docs :lower :merge :mock


text:
function ec_candle_init()
    opt = ec_init_cartesian()
    opt["legend"]["data"] = [""]
    opt["tooltip"]["trigger"] = "axis"
    opt["toolbox"]["feature"]["dataZoom"] = true
    opt["yAxis"] = [
        Dict("splitArea" => Dict("show" => true),
         "scale" => true,
         "boundaryGap" => [0.03,0.03],
         "splitNumber" => 9,
         "precision" => 2)
        ]
    opt["xAxis"] = [
        Dict("boundaryGap" => true)
        ]
    opt["dataZoom"] = Dict("show" => true,
                       "realtime" => true,
                       "start" => 50,
                       "end" => 100)

    return opt
end

# data=[open close low high]
# colors = [up_bar_body; down_bar_body;
#            up_bar_shadow; down_bar_shadow;
#            up_bar_selected; down_bar_selected]

_color_candle_US = ["green", "white", "green", "red", "green", "white"]
_color_candle_CN = ["red", "green", "red", "green", "white", "black"]

function ec_candle{S<:String, C<:String}(data::Array, labels::Vector{S};
colors::Vector{C}=_color_candle_US)
    opt = ec_candle_init()

    merge!(opt["xAxis"][1], Dict("data" => labels))

    cl = length(colors)
    if cl == 2
        colors = [colors; colors]
    elseif cl == 4
        colors = [colors; colors[1:2]]
    elseif cl != 6
        error("colors should be defined correctly")
    end

    data_dict =
    Dict(
        "name" => "",
        "type" => "k",

        "itemStyle" =>
        Dict(
            "normal" =>
            Dict(
                "color" => colors[1],
                "color0" => colors[2],
                "lineStyle" =>
                Dict(
                    "width" => 2,
                    "color" => colors[3],
                    "color0" => colors[4]
                )
            ),
            "emphasis" =>
            Dict(
                "color" => colors[5],
                "color0" => colors[6]
            )
        ),

        "data" => data.'
    )

    push!(opt["series"], data_dict)

    return opt
end


text:
__precompile__()
module ml_for_pt

# export package name as MLP
export MLP
const MLP = ml_for_pt

# load packages
using Flux
using StaticArrays
using UnPack
using StatsBase
using DelimitedFiles
using Random
using Base.Threads

# include additional files
include("supervised_learning.jl")
include("learning_by_confusion.jl")
include("prediction_based_method.jl")
include("distributions.jl")
include("utils.jl")

end


text:
module PlotTools

using PyPlot

function pointexclusion1D(ys, range=0.04, xspacing=0.04)
    incluster = BitArray(length(ys))
    clusters = BitArray[]
    for j in 1:length(ys)
        if !incluster[j]
            incluster[j] = true
            c = BitArray(length(ys))
            c[j] = true
            for k in 1:length(ys)
                if !incluster[k] && any(abs(ys[c] - ys[k]) .< range)
                    incluster[k] = true
                    c[k] = true
                end
            end

            push!(clusters, c)
        end
    end
    xs = zeros(ys)
    for c in clusters
        clen = countnz(c)
        if clen > 1
            offs = Float64[j * xspacing for j in 0:clen-1]
            offs = offs - mean(offs)
            xs[c] = offs
        end
    end
    xs
end

function pointexclusion1D{T}(xs::Vector{T}, ys::Vector{T}, yrange=0.04, xspacing=0.04, x_round_digits=6)
    xs_u = unique(round(xs, x_round_digits))

    xs_new = zeros(xs)

    for x in xs_u
        flt = abs(xs - x) .< 10.0^(-x_round_digits)
        ys_x = ys[flt]

        x_shifts = pointexclusion1D(ys_x, yrange, xspacing)
        xs_new[flt] = xs[flt] + x_shifts
    end

    xs_new, ys
end

function savepdf_cropped(fn, args...; kwargs...)
    savefig(fn, args...; kwargs...)
    run(`pdfcrop $fn $fn`)
end

function savepdf_cropped_compfonts(fn, args...; kwargs...)
    savefig(fn, args...; kwargs...)
    fntmp = fn * ".tmp"
    run(`pdfcrop $fn $fntmp`)
    run(`gs -sDEVICE=pdfwrite -dNOPAUSE -dQUIET -dBATCH -dCompressFonts=true -sOutputFile=$fn $fntmp`)
    run(`rm $fntmp`)
end

function savepdf_cropped_viaeps(fn, args...; kwargs...)
    fneps = fn * ".eps"
    savefig(fneps, args...; kwargs...)
    fnpdf = fn * ".pdf"
    run(`epstopdf $fneps $fnpdf`)
    run(`pdfcrop $fnpdf $fnpdf`)
end

end

max_stars_repo_name: amilsted/spinchain-virasoro-demo

max_issues_repo_name: amilsted/spinchain-virasoro-demo

max_forks_repo_name: amilsted/spinchain-virasoro-demo


text:
@testset "zp" begin
    @test iszero(ERFA._zp(randn(3)))
    @test iszero(fill!(randn(3), 0.0))
end

@testset "zpv" begin
    pv1 = [randn(3), randn(3)]
    pv2 = [randn(3), randn(3)]
    @test iszero(ERFA._zpv(pv1))
    @test iszero(fill!.(pv2, 0.0))
end

@testset "zr" begin
    r1 = randn(3, 3)
    r2 = randn(3, 3)
    @test iszero(ERFA._zr(r1))
    @test iszero(fill!(r2, 0.0))
end