text:
# math functionality

using Base: FastMath

using SpecialFunctions

## helpers

within(lower, upper) = (val) -> lower <= val <= upper

## trigonometric

@device_override Base.cos(x::Float64) = ccall("extern __nv_cos", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.cos(x::Float32) = ccall("extern __nv_cosf", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.cos_fast(x::Float32) = ccall("extern __nv_fast_cosf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.cospi(x::Float64) = ccall("extern __nv_cospi", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.cospi(x::Float32) = ccall("extern __nv_cospif", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.sin(x::Float64) = ccall("extern __nv_sin", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.sin(x::Float32) = ccall("extern __nv_sinf", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.sin_fast(x::Float32) = ccall("extern __nv_fast_sinf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.sinpi(x::Float64) = ccall("extern __nv_sinpi", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.sinpi(x::Float32) = ccall("extern __nv_sinpif", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.tan(x::Float64) = ccall("extern __nv_tan", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.tan(x::Float32) = ccall("extern __nv_tanf", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.tan_fast(x::Float32) = ccall("extern __nv_fast_tanf", llvmcall, Cfloat, (Cfloat,), x)

## inverse trigonometric

@device_override Base.acos(x::Float64) = ccall("extern __nv_acos", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.acos(x::Float32) = ccall("extern __nv_acosf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.asin(x::Float64) = ccall("extern __nv_asin", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.asin(x::Float32) = ccall("extern __nv_asinf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.atan(x::Float64) = ccall("extern __nv_atan", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.atan(x::Float32) = ccall("extern __nv_atanf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.atan(x::Float64, y::Float64) = ccall("extern __nv_atan2", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_override Base.atan(x::Float32, y::Float32) = ccall("extern __nv_atan2f", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

## hyperbolic

@device_override Base.cosh(x::Float64) = ccall("extern __nv_cosh", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.cosh(x::Float32) = ccall("extern __nv_coshf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.sinh(x::Float64) = ccall("extern __nv_sinh", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.sinh(x::Float32) = ccall("extern __nv_sinhf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.tanh(x::Float64) = ccall("extern __nv_tanh", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.tanh(x::Float32) = ccall("extern __nv_tanhf", llvmcall, Cfloat, (Cfloat,), x)

# TODO: enable once PTX > 7.0 is supported
# @device_override Base.tanh(x::Float16) = @asmcall("tanh.approx.f16 \$0, \$1", "=h,h", Float16, Tuple{Float16}, x)

## inverse hyperbolic

@device_override Base.acosh(x::Float64) = ccall("extern __nv_acosh", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.acosh(x::Float32) = ccall("extern __nv_acoshf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.asinh(x::Float64) = ccall("extern __nv_asinh", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.asinh(x::Float32) = ccall("extern __nv_asinhf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.atanh(x::Float64) = ccall("extern __nv_atanh", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.atanh(x::Float32) = ccall("extern __nv_atanhf", llvmcall, Cfloat, (Cfloat,), x)

## logarithmic

@device_override Base.log(x::Float64) = ccall("extern __nv_log", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.log(x::Float32) = ccall("extern __nv_logf", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.log_fast(x::Float32) = ccall("extern __nv_fast_logf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.log10(x::Float64) = ccall("extern __nv_log10", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.log10(x::Float32) = ccall("extern __nv_log10f", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.log10_fast(x::Float32) = ccall("extern __nv_fast_log10f", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.log1p(x::Float64) = ccall("extern __nv_log1p", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.log1p(x::Float32) = ccall("extern __nv_log1pf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.log2(x::Float64) = ccall("extern __nv_log2", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.log2(x::Float32) = ccall("extern __nv_log2f", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.log2_fast(x::Float32) = ccall("extern __nv_fast_log2f", llvmcall, Cfloat, (Cfloat,), x)

@device_function logb(x::Float64) = ccall("extern __nv_logb", llvmcall, Cdouble, (Cdouble,), x)
@device_function logb(x::Float32) = ccall("extern __nv_logbf", llvmcall, Cfloat, (Cfloat,), x)

@device_function ilogb(x::Float64) = ccall("extern __nv_ilogb", llvmcall, Int32, (Cdouble,), x)
@device_function ilogb(x::Float32) = ccall("extern __nv_ilogbf", llvmcall, Int32, (Cfloat,), x)

## exponential

@device_override Base.exp(x::Float64) = ccall("extern __nv_exp", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.exp(x::Float32) = ccall("extern __nv_expf", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.exp_fast(x::Float32) = ccall("extern __nv_fast_expf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.exp2(x::Float64) = ccall("extern __nv_exp2", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.exp2(x::Float32) = ccall("extern __nv_exp2f", llvmcall, Cfloat, (Cfloat,), x)
# TODO: enable once PTX > 7.0 is supported
# @device_override Base.exp2(x::Float16) = @asmcall("ex2.approx.f16 \$0, \$1", "=h,h", Float16, Tuple{Float16}, x)

@device_override Base.exp10(x::Float64) = ccall("extern __nv_exp10", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.exp10(x::Float32) = ccall("extern __nv_exp10f", llvmcall, Cfloat, (Cfloat,), x)
@device_override FastMath.exp10_fast(x::Float32) = ccall("extern __nv_fast_exp10f", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.expm1(x::Float64) = ccall("extern __nv_expm1", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.expm1(x::Float32) = ccall("extern __nv_expm1f", llvmcall, Cfloat, (Cfloat,), x)
@device_override Base.expm1(x::Float16) = Float16(CUDA.expm1(Float32(x)))

@device_override Base.ldexp(x::Float64, y::Int32) = ccall("extern __nv_ldexp", llvmcall, Cdouble, (Cdouble, Int32), x, y)
@device_override Base.ldexp(x::Float32, y::Int32) = ccall("extern __nv_ldexpf", llvmcall, Cfloat, (Cfloat, Int32), x, y)

## error

@device_override SpecialFunctions.erf(x::Float64) = ccall("extern __nv_erf", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erf(x::Float32) = ccall("extern __nv_erff", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfinv(x::Float64) = ccall("extern __nv_erfinv", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfinv(x::Float32) = ccall("extern __nv_erfinvf", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfc(x::Float64) = ccall("extern __nv_erfc", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfc(x::Float32) = ccall("extern __nv_erfcf", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfcinv(x::Float64) = ccall("extern __nv_erfcinv", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfcinv(x::Float32) = ccall("extern __nv_erfcinvf", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfcx(x::Float64) = ccall("extern __nv_erfcx", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfcx(x::Float32) = ccall("extern __nv_erfcxf", llvmcall, Cfloat, (Cfloat,), x)

## integer handling (bit twiddling)

@device_function brev(x::Union{Int32, UInt32}) =   ccall("extern __nv_brev", llvmcall, UInt32, (UInt32,), x)
@device_function brev(x::Union{Int64, UInt64}) =   ccall("extern __nv_brevll", llvmcall, UInt64, (UInt64,), x)

@device_function clz(x::Union{Int32, UInt32}) =
    assume(within(UInt32(0), UInt32(32)),
           ccall("extern __nv_clz", llvmcall, Int32, (UInt32,), x))
@device_function clz(x::Union{Int64, UInt64}) =
    assume(within(UInt64(0), UInt64(64)),
           ccall("extern __nv_clzll", llvmcall, Int32, (UInt64,), x))

@device_function ffs(x::Union{Int32, UInt32}) =
    assume(within(UInt32(0), UInt32(32)),
           ccall("extern __nv_ffs", llvmcall, Int32, (UInt32,), x))
@device_function ffs(x::Union{Int64, UInt64}) =
    assume(within(UInt64(0), UInt64(64)),
           ccall("extern __nv_ffsll", llvmcall, Int32, (UInt64,), x))

@device_function popc(x::Union{Int32, UInt32}) =
    assume(within(UInt32(0), UInt32(32)),
           ccall("extern __nv_popc", llvmcall, Int32, (UInt32,), x))
@device_function popc(x::Union{Int64, UInt64}) =
    assume(within(UInt64(0), UInt64(64)),
           ccall("extern __nv_popcll", llvmcall, Int32, (UInt64,), x))

@device_function byte_perm(x::Union{Int32, UInt32}, y::Union{Int32, UInt32}, z::Union{Int32, UInt32}) =
    ccall("extern __nv_byte_perm", llvmcall, Int32, (UInt32, UInt32, UInt32), x, y, z)

## floating-point handling

@device_override Base.isfinite(x::Float32) = (ccall("extern __nv_finitef", llvmcall, Int32, (Cfloat,), x)) != 0
@device_override Base.isfinite(x::Float64) = (ccall("extern __nv_isfinited", llvmcall, Int32, (Cdouble,), x)) != 0

@device_override Base.isinf(x::Float64) = (ccall("extern __nv_isinfd", llvmcall, Int32, (Cdouble,), x)) != 0
@device_override Base.isinf(x::Float32) = (ccall("extern __nv_isinff", llvmcall, Int32, (Cfloat,), x)) != 0

@device_override Base.isnan(x::Float64) = (ccall("extern __nv_isnand", llvmcall, Int32, (Cdouble,), x)) != 0
@device_override Base.isnan(x::Float32) = (ccall("extern __nv_isnanf", llvmcall, Int32, (Cfloat,), x)) != 0

@device_function nearbyint(x::Float64) = ccall("extern __nv_nearbyint", llvmcall, Cdouble, (Cdouble,), x)
@device_function nearbyint(x::Float32) = ccall("extern __nv_nearbyintf", llvmcall, Cfloat, (Cfloat,), x)

@device_function nextafter(x::Float64, y::Float64) = ccall("extern __nv_nextafter", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_function nextafter(x::Float32, y::Float32) = ccall("extern __nv_nextafterf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

## sign handling

@device_override Base.signbit(x::Float64) = (ccall("extern __nv_signbitd", llvmcall, Int32, (Cdouble,), x)) != 0
@device_override Base.signbit(x::Float32) = (ccall("extern __nv_signbitf", llvmcall, Int32, (Cfloat,), x)) != 0

@device_override Base.copysign(x::Float64, y::Float64) = ccall("extern __nv_copysign", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_override Base.copysign(x::Float32, y::Float32) = ccall("extern __nv_copysignf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@device_override Base.abs(x::Int32) =   ccall("extern __nv_abs", llvmcall, Int32, (Int32,), x)
@device_override Base.abs(f::Float64) = ccall("extern __nv_fabs", llvmcall, Cdouble, (Cdouble,), f)
@device_override Base.abs(f::Float32) = ccall("extern __nv_fabsf", llvmcall, Cfloat, (Cfloat,), f)
# TODO: enable once PTX > 7.0 is supported
# @device_override Base.abs(x::Float16) = @asmcall("abs.f16 \$0, \$1", "=h,h", Float16, Tuple{Float16}, x)
@device_override Base.abs(x::Int64) =   ccall("extern __nv_llabs", llvmcall, Int64, (Int64,), x)

## roots and powers

@device_override Base.sqrt(x::Float64) = ccall("extern __nv_sqrt", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.sqrt(x::Float32) = ccall("extern __nv_sqrtf", llvmcall, Cfloat, (Cfloat,), x)

@device_function rsqrt(x::Float64) = ccall("extern __nv_rsqrt", llvmcall, Cdouble, (Cdouble,), x)
@device_function rsqrt(x::Float32) = ccall("extern __nv_rsqrtf", llvmcall, Cfloat, (Cfloat,), x)
@device_function rsqrt(x::Float16) = Float16(rsqrt(Float32(x)))

@device_override Base.cbrt(x::Float64) = ccall("extern __nv_cbrt", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.cbrt(x::Float32) = ccall("extern __nv_cbrtf", llvmcall, Cfloat, (Cfloat,), x)

@device_function rcbrt(x::Float64) = ccall("extern __nv_rcbrt", llvmcall, Cdouble, (Cdouble,), x)
@device_function rcbrt(x::Float32) = ccall("extern __nv_rcbrtf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.:(^)(x::Float64, y::Float64) = ccall("extern __nv_pow", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_override Base.:(^)(x::Float32, y::Float32) = ccall("extern __nv_powf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)
@device_override FastMath.pow_fast(x::Float32, y::Float32) = ccall("extern __nv_fast_powf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)
@device_override Base.:(^)(x::Float64, y::Int32) = ccall("extern __nv_powi", llvmcall, Cdouble, (Cdouble, Int32), x, y)
@device_override Base.:(^)(x::Float32, y::Int32) = ccall("extern __nv_powif", llvmcall, Cfloat, (Cfloat, Int32), x, y)
@device_override Base.:(^)(x::Float64, y::Int64) = x ^ Float64(y)
@device_override Base.:(^)(x::Float32, y::Int64) = x ^ Float32(y)

## rounding and selection

# TODO: differentiate in return type, map correctly
#@device_override Base.round(x::Float64) = ccall("extern __nv_llround", llvmcall, Int64, (Cdouble,), x)
#@device_override Base.round(x::Float32) = ccall("extern __nv_llroundf", llvmcall, Int64, (Cfloat,), x)
#@device_override Base.round(x::Float64) = ccall("extern __nv_round", llvmcall, Cdouble, (Cdouble,), x)
#@device_override Base.round(x::Float32) = ccall("extern __nv_roundf", llvmcall, Cfloat, (Cfloat,), x)

# TODO: differentiate in return type, map correctly
#@device_override Base.rint(x::Float64) = ccall("extern __nv_llrint", llvmcall, Int64, (Cdouble,), x)
#@device_override Base.rint(x::Float32) = ccall("extern __nv_llrintf", llvmcall, Int64, (Cfloat,), x)
#@device_override Base.rint(x::Float64) = ccall("extern __nv_rint", llvmcall, Cdouble, (Cdouble,), x)
#@device_override Base.rint(x::Float32) = ccall("extern __nv_rintf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.trunc(x::Float64) = ccall("extern __nv_trunc", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.trunc(x::Float32) = ccall("extern __nv_truncf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.ceil(x::Float64) = ccall("extern __nv_ceil", llvmcall, Cdouble, (Cdouble,), x)
@device_override Base.ceil(x::Float32) = ccall("extern __nv_ceilf", llvmcall, Cfloat, (Cfloat,), x)

@device_override Base.floor(f::Float64) = ccall("extern __nv_floor", llvmcall, Cdouble, (Cdouble,), f)
@device_override Base.floor(f::Float32) = ccall("extern __nv_floorf", llvmcall, Cfloat, (Cfloat,), f)

#@device_override Base.min(x::Int32, y::Int32) = ccall("extern __nv_min", llvmcall, Int32, (Int32, Int32), x, y)
#@device_override Base.min(x::Int64, y::Int64) = ccall("extern __nv_llmin", llvmcall, Int64, (Int64, Int64), x, y)
#@device_override Base.min(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umin", llvmcall, Int32, (Int32, Int32), x, y))
#@device_override Base.min(x::UInt64, y::UInt64) = convert(UInt64, ccall("extern __nv_ullmin", llvmcall, Int64, (Int64, Int64), x, y))
@device_override Base.min(x::Float64, y::Float64) = ccall("extern __nv_fmin", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_override Base.min(x::Float32, y::Float32) = ccall("extern __nv_fminf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

#@device_override Base.max(x::Int32, y::Int32) = ccall("extern __nv_max", llvmcall, Int32, (Int32, Int32), x, y)
#@device_override Base.max(x::Int64, y::Int64) = ccall("extern __nv_llmax", llvmcall, Int64, (Int64, Int64), x, y)
#@device_override Base.max(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umax", llvmcall, Int32, (Int32, Int32), x, y))
#@device_override Base.max(x::UInt64, y::UInt64) = convert(UInt64, ccall("extern __nv_ullmax", llvmcall, Int64, (Int64, Int64), x, y))
@device_override Base.max(x::Float64, y::Float64) = ccall("extern __nv_fmax", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_override Base.max(x::Float32, y::Float32) = ccall("extern __nv_fmaxf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@device_function saturate(x::Float32) = ccall("extern __nv_saturatef", llvmcall, Cfloat, (Cfloat,), x)

## division and remainder

# NOTE: CUDA follows fmod, which behaves differently than Base.mod for negative numbers
#@device_override Base.mod(x::Float64, y::Float64) = ccall("extern __nv_fmod", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
#@device_override Base.mod(x::Float32, y::Float32) = ccall("extern __nv_fmodf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@device_override Base.rem(x::Float64, y::Float64) = ccall("extern __nv_remainder", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_override Base.rem(x::Float32, y::Float32) = ccall("extern __nv_remainderf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)
@device_override Base.rem(x::Float16, y::Float16) = Float16(rem(Float32(x), Float32(y)))

@device_override FastMath.div_fast(x::Float32, y::Float32) = ccall("extern __nv_fast_fdividef", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

## gamma function

@device_override SpecialFunctions.lgamma(x::Float64) = ccall("extern __nv_lgamma", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.lgamma(x::Float32) = ccall("extern __nv_lgammaf", llvmcall, Cfloat, (Cfloat,), x)

@device_function tgamma(x::Float64) = ccall("extern __nv_tgamma", llvmcall, Cdouble, (Cdouble,), x)
@device_function tgamma(x::Float32) = ccall("extern __nv_tgammaf", llvmcall, Cfloat, (Cfloat,), x)

## Bessel

@device_override SpecialFunctions.besselj0(x::Float64) = ccall("extern __nv_j0", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.besselj0(x::Float32) = ccall("extern __nv_j0f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.besselj1(x::Float64) = ccall("extern __nv_j1", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.besselj1(x::Float32) = ccall("extern __nv_j1f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.besselj(n::Int32, x::Float64) = ccall("extern __nv_jn", llvmcall, Cdouble, (Int32, Cdouble), n, x)
@device_override SpecialFunctions.besselj(n::Int32, x::Float32) = ccall("extern __nv_jnf", llvmcall, Cfloat, (Int32, Cfloat), n, x)

@device_override SpecialFunctions.bessely0(x::Float64) = ccall("extern __nv_y0", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.bessely0(x::Float32) = ccall("extern __nv_y0f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.bessely1(x::Float64) = ccall("extern __nv_y1", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.bessely1(x::Float32) = ccall("extern __nv_y1f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.bessely(n::Int32, x::Float64) = ccall("extern __nv_yn", llvmcall, Cdouble, (Int32, Cdouble), n, x)
@device_override SpecialFunctions.bessely(n::Int32, x::Float32) = ccall("extern __nv_ynf", llvmcall, Cfloat, (Int32, Cfloat), n, x)

## distributions

# TODO: override StatsFun.jl?

@device_function normcdf(x::Float64) = ccall("extern __nv_normcdf", llvmcall, Cdouble, (Cdouble,), x)
@device_function normcdf(x::Float32) = ccall("extern __nv_normcdff", llvmcall, Cfloat, (Cfloat,), x)

@device_function normcdfinv(x::Float64) = ccall("extern __nv_normcdfinv", llvmcall, Cdouble, (Cdouble,), x)
@device_function normcdfinv(x::Float32) = ccall("extern __nv_normcdfinvf", llvmcall, Cfloat, (Cfloat,), x)

#
# Unsorted
#

@device_override Base.hypot(x::Float64, y::Float64) = ccall("extern __nv_hypot", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_override Base.hypot(x::Float32, y::Float32) = ccall("extern __nv_hypotf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@device_override Base.fma(x::Float64, y::Float64, z::Float64) = ccall("extern __nv_fma", llvmcall, Cdouble, (Cdouble, Cdouble, Cdouble), x, y, z)
@device_override Base.fma(x::Float32, y::Float32, z::Float32) = ccall("extern __nv_fmaf", llvmcall, Cfloat, (Cfloat, Cfloat, Cfloat), x, y, z)
@device_override Base.fma(x::Float16, y::Float16, z::Float16) = ccall("llvm.fma.f16", llvmcall, Float16, (Float16, Float16, Float16), x, y, z)

@device_function sad(x::Int32, y::Int32, z::Int32) = ccall("extern __nv_sad", llvmcall, Int32, (Int32, Int32, Int32), x, y, z)
@device_function sad(x::UInt32, y::UInt32, z::UInt32) = convert(UInt32, ccall("extern __nv_usad", llvmcall, Int32, (Int32, Int32, Int32), x, y, z))

@device_function dim(x::Float64, y::Float64) = ccall("extern __nv_fdim", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@device_function dim(x::Float32, y::Float32) = ccall("extern __nv_fdimf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@device_function mul24(x::Int32, y::Int32) = ccall("extern __nv_mul24", llvmcall, Int32, (Int32, Int32), x, y)
@device_function mul24(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umul24", llvmcall, Int32, (Int32, Int32), x, y))

@device_function mul64hi(x::Int64, y::Int64) = ccall("extern __nv_mul64hi", llvmcall, Int64, (Int64, Int64), x, y)
@device_function mul64hi(x::UInt64, y::UInt64) = convert(UInt64, ccall("extern __nv_umul64hi", llvmcall, Int64, (Int64, Int64), x, y))
@device_function mulhi(x::Int32, y::Int32) = ccall("extern __nv_mulhi", llvmcall, Int32, (Int32, Int32), x, y)
@device_function mulhi(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umulhi", llvmcall, Int32, (Int32, Int32), x, y))

@device_function hadd(x::Int32, y::Int32) = ccall("extern __nv_hadd", llvmcall, Int32, (Int32, Int32), x, y)
@device_function hadd(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_uhadd", llvmcall, Int32, (Int32, Int32), x, y))

@device_function rhadd(x::Int32, y::Int32) = ccall("extern __nv_rhadd", llvmcall, Int32, (Int32, Int32), x, y)
@device_function rhadd(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_urhadd", llvmcall, Int32, (Int32, Int32), x, y))

@device_function scalbn(x::Float64, y::Int32) = ccall("extern __nv_scalbn", llvmcall, Cdouble, (Cdouble, Int32), x, y)
@device_function scalbn(x::Float32, y::Int32) = ccall("extern __nv_scalbnf", llvmcall, Cfloat, (Cfloat, Int32), x, y)


text:
# Activating Julia environment
using Pkg
Pkg.activate(".")

# Packages
using JuMP, Gurobi, Random, Distributions, LinearAlgebra, CSV, DataFrames

# CRUNCHING THE DATA
# NODES (I = International supplier; R = Regional supplier; L = Local market (both supply and deliver); D = delivery point)
N = []
N_I = []     # set of international suppliers
N_R = []     # set of regional suppliers
N_L = []     # set of local markets
N_D = []     # delivery points
dem = Dict() # set of demands
file = CSV.File("syria_nodes.csv")
for row in file
    push!(N, row.Name)
    if !ismissing(row.Demand)
        dem[row.Name] = row.Demand
    end
    if row.Type == "I"
        push!(N_I, row.Name)
    elseif row.Type == "R"
        push!(N_R, row.Name)
    elseif row.Type == "L"
        push!(N_L, row.Name) # Note: local markets supply and deliver goods at the given cost.
    elseif row.Type == "D"
        push!(N_D, row.Name)
    else
        throw(ErrorException("rowType $(row.Type) not supported."))
    end
end

# EDGES
hc = DataFrame(CSV.File("syria_edges.csv"))

# FOOD NUTRITION AND INTERNATIONAL COSTS
fooddata = DataFrame(CSV.File("syria_foodnutrition.csv"))
intfoodcosts = select(fooddata, [:Food, :InternationalPrice])
commodities = sort(Array(intfoodcosts.Food)) # Commodities
select!(fooddata, Not([14,15]))
fooddata = Dict(fooddata.Food .=> eachrow(fooddata))
# Note: fooddata contains the nutrients provided by 100g of a commodity!

# FOOD COST ($/metric ton for regional suppliers)
pc = DataFrame(CSV.File("syria_foodcost.csv"))
pc.Food = convert.(String63, pc.Food)
for int_supply_node in N_I # adding international prices to pc for easier processing
    for row in eachrow(intfoodcosts)
        append!(pc, DataFrame(:A => N_I, :Food => row.Food, :Price => row.InternationalPrice))
    end
end
pc = unique(pc)
international_items = DataFrame([r for r in eachrow(pc) if r.A in N_I])
regional_items =  DataFrame([r for r in eachrow(pc) if r.A in N_R])

# FOOD REQUIREMENTS (avg. per person per day)
foodreqs = DataFrame(CSV.File("syria_foodreq.csv"))
select!(foodreqs, Not(:Type))
nutrients = String.(propertynames(foodreqs))
foodreqs = Dict(string(pptname) => foodreqs[1, pptname] for pptname in propertynames(foodreqs))

# FINALLY CREATING THE MODEL
m = Model(Gurobi.Optimizer)

# Procurement and delivery
procurement_links = unique([row.A => row.Food for row in eachrow(pc)])       # all places where we can procure food
@variable(m, procurement[A = N, Food = commodities; (A => Food) in procurement_links] >= 0) # procurement in tons
@variable(m, delivery[N_D, commodities] >= 0)                                               # delivery in tons

# Total procurement cost
@variable(m, procurement_cost >= 0)
@constraint(m, procurement_cost == sum(r[:Price] * procurement[r.A, r.Food] for r in eachrow(pc)))

# Transportation
transportation_links = unique([row.A => row.B for row in eachrow(hc)]) # all possible edges
@variable(m, transportation[A = N, B = N; (A => B) in transportation_links] >= 0)     # transportation in tons...
@variable(m, F[A = N, B = N, W = commodities; (A => B) in transportation_links] >= 0) # linked directly to F, also in tons.
for r in eachrow(hc) # Linking transportation cost to total food transported on an edge
    @constraint(m, transportation[r.A, r.B] == sum(F[r.A, r.B, commodity] for commodity in commodities))
end
# Total transportation cost
@variable(m, transportation_cost >= 0)
@constraint(m, transportation_cost == sum(r.tCost * transportation[r.A, r.B] for r in eachrow(hc)))

# Flow constraints
for node in N
    valid_sources = [link.first for link in transportation_links if link.second == node]
    valid_sinks = [link.second for link in transportation_links if link.first == node]
    for commodity in commodities
        if (node =>  commodity) in procurement_links
            @constraint(m, procurement[node, commodity] + sum(F[source, node, commodity] for source in valid_sources) ==
                                sum(F[node, sink, commodity] for sink in valid_sinks))
        elseif node in N_D
            @assert length(valid_sinks) == 0
            @constraint(m, delivery[node, commodity] == sum(F[source, node, commodity] for source in valid_sources))
        else
            @constraint(m, sum(F[source, node, commodity] for source in valid_sources) ==
                                sum(F[node, sink, commodity] for sink in valid_sinks))
        end
    end
end

# Serving demand
@variable(m, ration_pp[commodities] >= 0) # Rations (kg/person) of commodities
@variable(m, nutrients_pp[nutrients] >= 0) # Total nutrients per person

# Making sure the rations are good nutritionally
for nutrient in nutrients # Note the factor of 10 for conversion of 100g to kg (since rations are in kg/pp)
    @constraint(m, nutrients_pp[nutrient] == 10 * sum(ration_pp[commodity] * fooddata[commodity][nutrient] for commodity in commodities))
    @constraint(m, nutrients_pp[nutrient] >= foodreqs[nutrient])
end
for node in N_D
    for commodity in commodities
        @constraint(m, 1000*delivery[node, commodity] >= dem[node] * ration_pp[commodity])
    end
end

# Diet constraints
# Achieving a greater than 4:1 ratio by mass of carbohydrates to protein, and a greater than 4:1 ratio by mass of carbohydrates to fats
# for the same total energy intake. The energy stored in a gram of carbohydrate, protein and fat are 4kcal/g, 4kcal/g and 9kcal/g respectively.
@variable(m, carbs_pp >= 0)
@constraint(m, 4*carbs_pp == nutrients_pp["Energy(kcal)"] - 4*nutrients_pp["Protein(g)"] - 9*nutrients_pp["Fat(g)"])
@constraint(m, carbs_pp >= 4*nutrients_pp["Protein(g)"])
@constraint(m, carbs_pp >= 4*nutrients_pp["Fat(g)"])

# Setting objectives
@objective(m, Min, procurement_cost + transportation_cost)

# Solving
optimize!(m)

ext: jl


text:
using ActifMPCCmod

f(x)=x[1]-x[2]
r=1.0;s=1.0;t=1.0;
ActifMPCCmod_success=true
println("On teste le module ActifMPCCmod")

nlp=ADNLPModel(f,ones(2), lvar=-Inf*ones(2), uvar=Inf*ones(2), c=c,y0=ones(1),lcon=lcon,ucon=Inf*ones(1))
mpcc = MPCCmod.MPCC(nlp,G,H,1)

##############################################################
#Test constructeur MPCC_actif :
nlp_ma=ADNLPModel(x->mpcc.mp.f(x), [mpcc.mp.meta.x0;-r;1],lvar=[mpcc.mp.meta.lvar;-r*ones(2*mpcc.nb_comp)], uvar=[mpcc.mp.meta.uvar;Inf*ones(2*mpcc.nb_comp)])
ma=ActifMPCCmod.MPCC_actif(nlp_ma,r,s,t,1)
##############################################################
#Test setw et updatew

w_save=sparse(zeros(2,2))
w_save[1,1]=0.0
w_save[1,2]=1.0
ActifMPCCmod.setw(ma,w_save)
ActifMPCCmod.updatew(ma)

w_save2=sparse(zeros(2,2))
w_save2[2,1]=1.0
w_save2[1,2]=0.0
ma.w=w_save2
ActifMPCCmod.updatew(ma)

##############################################################
# Exemple sans contrainte active :
#println("Exemple 1")
xk=[0.8;0.6;0.5;0.45]
dxk=[0.0;0.0;-1.0;0.0]
ActifMPCCmod.setw(ma,sparse(zeros(2,2)))

 #[alpha,w_save]=ActifMPCCmod.PasMax(ma,xk,dxk)
if !(ma.bar_w==[1;2] && ActifMPCCmod.evalx(ma,xk)==xk && norm(ActifMPCCmod.obj(ma,xk)-0.2)<precision && ActifMPCCmod.grad(ma,xk)==[1;-1;0;0] && norm(ActifMPCCmod.PasMax(ma,xk,dxk)[1]-0.3525,Inf)<precision)
 println("Exemple 1")
 println(ma.bar_w) #[1,2]
 println(ActifMPCCmod.evalx(ma,xk)) #doit rentre la même chose que xk
 println(ActifMPCCmod.obj(ma,xk)) #doit donner xk[1]-xk[2]
 println(ActifMPCCmod.grad(ma,xk)) #doit donner [1,-1,0,0]
 println(ActifMPCCmod.PasMax(ma,xk,dxk))

 ActifMPCCmod_success=false
end

#println("Exemple 2")
nlp_ma2=ADNLPModel(x->mpcc.mp.f(x)+x[3]^2+x[4]^2, [1;0;-r;1],lvar=[mpcc.mp.meta.lvar;-r*ones(2*mpcc.nb_comp)], uvar=[mpcc.mp.meta.uvar;Inf*ones(2*mpcc.nb_comp)])
ma2=ActifMPCCmod.MPCC_actif(nlp_ma2,r,s,t,1)

if !(ma2.w13c==[] && ma2.w24c==[1] && ActifMPCCmod.evalx(ma2,[1;0;1])==[1;0;-r;1] && ActifMPCCmod.obj(ma2,[1;0;1])==3 && ActifMPCCmod.grad(ma2,[1;0;1])==[1;-1;2] && norm(ActifMPCCmod.PasMax(ma2,[1;0;1],[0.0;0.0;-1.0])[1]-1.0,Inf)<precision)
 println("Exemple 2")
 println(ma2.w13c) # []
 println(ma2.w24c) # [1]
 println(ActifMPCCmod.evalx(ma2,[1;0;1])) #doit rentre la même chose que xkp
 println(ActifMPCCmod.obj(ma2,[1;0;1])) #doit donner xkp[1]-xkp[2]+xkp[3]^2+xkp[4]^2=3
 #gradient : [1;-1;2x3;2x4]=[1;-1;-2;2]
 println(ActifMPCCmod.grad(ma2,[1;0;1])) #doit donner [1,-1,2]
 println(ActifMPCCmod.PasMax(ma2,[1;0;1],[0.0;0.0;-1.0]))

 ActifMPCCmod_success=false
end

#println("Exemple 3")
t=0.0;s=2.0;r=2.0;
nlp_ma3=ADNLPModel(x->mpcc.mp.f(x)+x[3]^2+x[4]^2, [1;0;s;3],lvar=[-Inf;-Inf;-r;-r],uvar=ones(4)*Inf)
ma3=ActifMPCCmod.MPCC_actif(nlp_ma3,r,s,t,1)

if !(ma3.w[2,1]==1 && ActifMPCCmod.evalx(ma3,[1;0;3])==[1;0;s;3] && ActifMPCCmod.obj(ma3,[1;0;3])==14 && ActifMPCCmod.grad(ma3,[1;0;3])==[1;-1;6])
 println("Exemple 3")
 println(ma3.w) # (2,1)=1
 println(ActifMPCCmod.evalx(ma3,[1;0;3])) #doit rentre la même chose que xkf
 println(ActifMPCCmod.obj(ma3,[1;0;3])) #doit donner xkf[1]-xkf[2]+xkf[3]^2+xkf[4]^2=14
 #gradient : [1;-1;2x3;2x4]=[1;-1;4;6]
 println(ActifMPCCmod.grad(ma3,[1;0;3])) #doit donner [1,-1,6]

 ActifMPCCmod_success=false
end

#println("Exemple 4")
t=0.0;s=2.0;r=2.0;
nlp_ma4=ADNLPModel(x->mpcc.mp.f(x)+x[3]^2+x[4]^2, [1;0;0;s],lvar=[-Inf;-Inf;-r;-r],uvar=ones(4)*Inf)
ma4=ActifMPCCmod.MPCC_actif(nlp_ma4,r,s,t,1)
xkf=[1;0;0;s]
if !(ma4.w[2,2]==1 && ActifMPCCmod.evalx(ma4,[1;0;0])==xkf && norm(ActifMPCCmod.obj(ma4,[1;0;0])-5)<precision && ActifMPCCmod.grad(ma4,[1;0;0])==[1;-1;0])
 println("Exemple 4")
 println(ma4.w) # (2,2)=1
 println(ActifMPCCmod.evalx(ma4,[1;0;0])) #doit rentre la même chose que xkf
 println(ActifMPCCmod.obj(ma4,[1;0;0])) #doit donner xkf[1]-xkf[2]+xkf[3]^2+xkf[4]^2=5
 #gradient : [1;-1;2x3;2x4]=[1;-1;0;4]
 println(ActifMPCCmod.grad(ma4,[1;0;0])) #doit donner [1,-1,0]

 ActifMPCCmod_success=false
end

#Bilan : donner une sortie bilan
if ActifMPCCmod_success==true
 println("ActifMPCCmod.jl passes the test !")
else
 println("ActifMPCCmod.jl contains some error")
end


text:
discretization(ss::SamplingStyle, ap::ApproximationProblem; options...) =
    discretization(ss, ap, haskey(options,:S) ? options[:S] : samplingoperator(ss, ap; options...); options...)
default_discretization(::SamplingStyle, dict::Dictionary, S::AbstractOperator; options...) =
    apply(S, dict; options...)
discretization(ss::ProductSamplingStyle, ap::ApproximationProblem, S::AbstractOperator; options...) =
    tensorproduct( elements(discretization, ss, ap,
        BasisFunctions.productelements(S); options...)... )

# The discretization routine can also take f as an argument
discretization(f, dict::Dictionary, args...; options...) =
    discretization(f, approximationproblem(dict, args...); options...)
discretization(f, platform::Platform, args...; options...) =
    discretization(f, approximationproblem(platform, args...); options...)
discretization(f, ap::ApproximationProblem; samplingstyle=SamplingStyle(ap), options...) =
    discretization(f, samplingstyle, ap; options...)
discretization(f, ss::SamplingStyle, ap::ApproximationProblem; options...) =
    discretization(f, ss, ap, haskey(options,:S) ? options[:S] : samplingoperator(ss, ap; options...); options...)
function discretization(f, ss::SamplingStyle, ap::ApproximationProblem, S::AbstractOperator; options...)
    A = discretization(ss, ap, S; options...)
    B = apply(S, f; options...)
    A, B
end
function discretization(f::AbstractArray, ss::SamplingStyle, ap::ApproximationProblem, S::AbstractOperator; verbose = false, options...)
    A = discretization(ss, ap, S; options...)
    if size(S,1) == length(f)
        verbose && println("Discretization: sampled data is given, omitting sampling operation")
        A, f
    else
        error("Sample data is given, but the array size is not compatible with the sampling operator.")
    end
end

function dualdiscretization(ss::SamplingStyle, ap::ApproximationProblem; options...)
    dualdict = haskey(options, :dualdict) ? options[:dualdict] : azdual_dict(ss, ap; options...)
    S = haskey(options, :S) ? options[:S] : dualsamplingoperator(ss, ap; dualdict=dualdict, options...)
    dualdiscretization(ss, ap, S, dualdict; options...)
end

function dualdiscretization(ss::SamplingStyle, ap::ApproximationProblem, S::AbstractOperator; options...)
    dualdict = haskey(options, :dualdict) ? options[:dualdict] : azdual_dict(ss, ap; options...)
    dualdiscretization(ss, ap, S, dualdict; options...)
end

default_dualdiscretization(ss::SamplingStyle, dict::Dictionary, Stilde::AbstractOperator, dualdict::Dictionary; options...) =
    apply(Stilde, dualdict; options...)
dualdiscretization(ss::ProductSamplingStyle, ap::ApproximationProblem, S::AbstractOperator, dualdict::TensorProductDict; options...) =
    tensorproduct( elements(dualdiscretization, ss, ap,
        BasisFunctions.productelements(S), elements(dualdict); options...)... )


text:
using CircGeometry, Test

# error for too large vf
vf = 0.9
@test_throws ErrorException CircGeometry.check_vf(0.85)

# error for max iterations reached
vf = 0.79
n_bodies = 400
material = MaterialParameters(vf,n_bodies)
radius = 1.5
center = Point(-0.5,1.0)
outline = OutlineCircle(radius,center)
between_buffer = 100
@test_throws MethodError generate_porous_structure(outline,material,between_buffer;log=false)

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
# # Symmetry reducting a strong relaxation of the quadratic assigment problem
# Here, we are going to show how to load a QAP from QABLib, formulate
# a strong semidefinite relaxation of it, symmetry reduce it, and finally solve it.

# ## Quadratic assigment problems
# QAPs are given by two quadratic matrices $A$ and $B$. The objective is to permute
# the rows and columns of $B$, such that the inner product between the matrices is
# minimized.
#
# ``\mathrm{QAP}(A,B) = \min_{\phi\in S_n} \sum_{i,j=1}^n a_{ij}b_{\phi(i)\phi(j)}``
#
# QAPs are notoriously hard to solve exactly, but there exist strong polynomial time
# relaxations, such as the following semidefinite programming relaxation:
#
# ```math
# \begin{aligned}
# \min\enspace & \langle B\otimes A ,Y\rangle\\
# \mathrm{s.t.}\enspace & \langle I_n\otimes E_{jj},Y\rangle=1 \text{ for }j\in [n],\\
# & \langle E_{jj}\otimes I_n,Y\rangle=1 \text{ for }j\in [n],\\
# & \langle I_n\otimes (J_n-I_n)+(J_n-I_n)\otimes I_n,Y\rangle =0, \\
# & \langle J_{n^2},Y\rangle = n^2,\\
# & Y\in D^{n^2},
# \end{aligned}
# ```
#
# But in practice this relaxation is often too big to be solved directly.

# ## Loading the data of a QAP
using SparseArrays, LinearAlgebra
using Test #src

file = joinpath(@__DIR__, "esc16j.dat")
data = open(file) do file
    read(file, String)
end
data = split(data, [' ', '\n', '\r'], keepempty = false)

n = parse(Int64, data[1])
A = zeros(Int64, n, n)
B = zeros(Int64, n, n)

pos = 2
for x = 1:n
    for y = 1:n
        A[x, y] = parse(Int64, data[pos])
        global pos += 1
    end
end
for x = 1:n
    for y = 1:n
        B[x, y] = parse(Int64, data[pos])
        global pos += 1
    end
end

# ## Building the SDP (in vectorized standard form)
n = size(A, 1)

## Objective
CPrg = sparse(kron(B, A))

In = sparse(I, n, n)
Jn = ones(n, n)

## Vectorised constraint matrices as rows of large matrix APrg
APrg = spzeros(2n + 1, n^4)
bPrg = zeros(2n + 1)
currentRow = 1

for j = 1:n
    Ejj = spzeros(n, n)
    Ejj[j, j] = 1.0
    APrg[currentRow, :] = vec(kron(In, Ejj))
    bPrg[currentRow] = 1
    global currentRow += 1
    ## Last constraint is linearly dependent on others
    if (j < n)
        APrg[currentRow, :] = vec(kron(Ejj, In))
        bPrg[currentRow] = 1
        global currentRow += 1
    end
end

APrg[currentRow, :] = vec(kron(In, Jn - In) + kron(Jn - In, In))
bPrg[currentRow] = 0
currentRow += 1
APrg[currentRow, :] = vec(ones(n^2, n^2))
bPrg[currentRow] = n^2

CPrg = sparse(vec(0.5 * (CPrg + CPrg')));

# ## Symmetry reducing the SDP

# We first determine an optimal admissible partition subspace
using SDPSymmetryReduction
P = admPartSubspace(CPrg, APrg, bPrg, true)
P.n

@test P.n == 150 #src
# And then we block-diagonalize it
blkD = blockDiagonalize(P, true);
@test sort(blkD.blkSizes) == sort([7, 7, 7, 7, 7, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) #src

# ## Determining the coefficients of the reduced SDP
PMat = hcat([sparse(vec(P.P .== i)) for i = 1:P.n]...)
newA = APrg * PMat
newB = bPrg
newC = CPrg' * PMat;

## Removing linearly dependent constraints #src
# using RowEchelon  #src
# T = rref!(Matrix(hcat(newA,newB))) #src
# r = rank(T) #src
# newA = T[1:r,1:end-1] #src
# newB = T[1:r,end] #src
# length(newB) #src

# ## Solving the reduced SDP with JuMP and CSDP

using JuMP, CSDP
m = Model(CSDP.Optimizer)

## Initialize variables corresponding parts of the partition P
## >= 0 because the original SDP-matrices are entry-wise nonnegative
x = @variable(m, x[1:P.n] >= 0)

@constraint(m, newA * x .== newB)
@objective(m, Max, newC * x)

psdBlocks = sum(blkD.blks[i] .* x[i] for i = 1:P.n)
for blk in psdBlocks
    if size(blk, 1) > 1
        @constraint(m, blk in PSDCone())
    else
        @constraint(m, blk .>= 0)
    end
end

optimize!(m)
#
termination_status(m)
#
objective_value(m)

@test objective_value(m) ≈ 48.3042 atol = 5#src


text:
using AstroTime
using AstroTime.Epochs
using Base.Test
using ERFA

AstroTime.update()

function fractionofday(dt)
    Dates.hour(dt)/24 + Dates.minute(dt)/(24*60) + Dates.second(dt)/86400 + Dates.millisecond(dt)/8.64e7
end

@testset "AstroTime" begin
    @testset "Periods" begin
        s = 1.0seconds
        m = 1.0minutes
        h = 1.0hours
        d = 1.0days
        y = 1.0years
        c = 1.0centuries
        @test s == Period{Second}(1.0)
        @test m == Period{Minute}(1.0)
        @test h == Period{Hour}(1.0)
        @test d == Period{Day}(1.0)
        @test y == Period{Year}(1.0)
        @test c == Period{Century}(1.0)

        @test seconds(s) == 1.0seconds
        @test seconds(m) == 60.0seconds
        @test seconds(h) == 3600.0seconds
        @test seconds(d) == 86400.0seconds
        @test seconds(y) == 3.15576e7seconds
        @test seconds(c) == 3.15576e9seconds

        @test minutes(s) == (1.0 / 60.0)minutes
        @test minutes(m) == 1.0minutes
        @test minutes(h) == 60.0minutes
        @test minutes(d) == 1440.0minutes
        @test minutes(y) == 525960.0minutes
        @test minutes(c) == 5.2596e7minutes

        @test hours(s) == (1.0 / 3600.0)hours
        @test hours(m) == (1.0 / 60.0)hours
        @test hours(h) == 1.0hours
        @test hours(d) == 24.0hours
        @test hours(y) == 8766.0hours
        @test hours(c) == 876600.0hours

        @test days(s) == (1.0 / 86400.0)days
        @test days(m) == (1.0 / 1440.0)days
        @test days(h) == (1.0 / 24.0)days
        @test days(d) == 1.0days
        @test days(y) == 365.25days
        @test days(c) == 36525.0days

        @test years(s) == (1.0 / 3.15576e7)years
        @test years(m) == (1.0 / 525960.0)years
        @test years(h) == (1.0 / 8766.0)years
        @test years(d) == (1.0 / 365.25)years
        @test years(y) == 1.0years
        @test years(c) == 100.0years

        @test centuries(s) == (1.0 / 3.15576e9)centuries
        @test centuries(m) == (1.0 / 5.2596e7)centuries
        @test centuries(h) == (1.0 / 876600.0)centuries
        @test centuries(d) == (1.0 / 36525.0)centuries
        @test centuries(y) == (1.0 / 100.0)centuries
        @test centuries(c) == 1.0centuries
    end
    @testset "Epoch Type" begin
        ep0 = TTEpoch(0.0)

        ep1 = ep0 - 1days
        @test ep1.jd1 == -1
        ep1 = ep0 - 2minutes
        @test ep1.jd2 == -2 / 1440

        ep1 = ep0 + 1days
        @test ep1.jd1 == 1
        ep1 = ep0 + 2minutes
        @test ep1.jd2 == 2 / 1440

        ep0 = TTEpoch(2000, 1, 1)
        @test ep0.jd1 == J2000 - 0.5

        @test TTEpoch(2000,1,1) < TTEpoch(2001,1,1)
        @test TTEpoch(2000,1,1) <= TTEpoch(2001,1,1)
        @test TTEpoch(2001,1,1) > TTEpoch(2000,1,1)
        @test TTEpoch(2001,1,1) >= TTEpoch(2000,1,1)

        ep0 = TTEpoch(2000,1,1,12,0,0,123)
        @test string(ep0) == "2000-01-01T12:00:00.123 TT"
    end
    @testset "Constructors" begin
        dt = DateTime(2000, 1, 1, 12, 0, 0.0)
        tt = TTEpoch(2000, 1, 1, 12, 0, 0.0)
        @test string(tt) == "2000-01-01T12:00:00.000 TT"

        tdb = TDBEpoch(tt)
        tcb = TCBEpoch(tt)
        tcg = TCGEpoch(tt)
        tai = TAIEpoch(tt)
        utc = UTCEpoch(tt)
        ut1 = UT1Epoch(tt)

        @test TTEpoch(J2000) ≈ tt
        @test jd2000(tt) == 0
        @test jd1950(TTEpoch(1950, 1, 1, 12)) == 0
        @test mjd(TTEpoch(1858, 11, 17)) == 0

        @test seconds(tt + 1seconds, J2000) ≈ 1.0seconds
        @test minutes(tt + 1minutes, J2000) ≈ 1.0minutes
        @test hours(tt + 1hours, J2000) ≈ 1.0hours
        @test days(tt + 1days, J2000) ≈ 1.0days
        @test years(tt + 1years, J2000) ≈ 1.0years
        @test centuries(tt + 1centuries, J2000) ≈ 1.0centuries

        @test tai ≈ TAIEpoch(utc)
        @test utc ≈ UTCEpoch(tai)
        @test utc ≈ UTCEpoch(ut1)
        @test ut1 ≈ UT1Epoch(utc)
        @test tai ≈ TAIEpoch(ut1)
        @test ut1 ≈ UT1Epoch(tai)
        @test tt ≈ TTEpoch(ut1)
        @test ut1 ≈ UT1Epoch(tt)
        @test tt ≈ TTEpoch(tai)
        @test tai ≈ TAIEpoch(tt)
        @test tt ≈ TTEpoch(tcg)
        @test tcg ≈ TCGEpoch(tt)
        @test tt ≈ TTEpoch(tdb)
        @test tdb ≈ TDBEpoch(tt)
        @test tdb ≈ TDBEpoch(tcb)
        @test tcb ≈ TCBEpoch(tdb)

        @test tt ≈ TTEpoch(tcb)
        @test tcb ≈ TCBEpoch(tt)
        @test tt == TTEpoch(tt)

        # Reference values from Orekit
        ref = TDBEpoch(2013, 3, 18, 12)
        @test UT1Epoch(ref) == UT1Epoch("2013-03-18T11:58:52.994")
        @test UTCEpoch(ref) == UTCEpoch("2013-03-18T11:58:52.814")
        @test TAIEpoch(ref) == TAIEpoch("2013-03-18T11:59:27.814")
        @test TTEpoch(ref) == TTEpoch("2013-03-18T11:59:59.998")
        @test TCBEpoch(ref) == TCBEpoch("2013-03-18T12:00:17.718")
        @test TCGEpoch(ref) == TCGEpoch("2013-03-18T12:00:00.795")
        @test ref == TDBEpoch(UT1Epoch("2013-03-18T11:58:52.994"))
        @test ref == TDBEpoch(UTCEpoch("2013-03-18T11:58:52.814"))
        @test ref == TDBEpoch(TAIEpoch("2013-03-18T11:59:27.814"))
        @test ref == TDBEpoch(TTEpoch("2013-03-18T11:59:59.998"))
        @test ref == TDBEpoch(TCBEpoch("2013-03-18T12:00:17.718"))
        @test ref == TDBEpoch(TCGEpoch("2013-03-18T12:00:00.795"))
    end
    @testset "Conversions" begin
        tai = TAIEpoch(2000, 1, 1, 12, 0, 0.0)
        utc = UTCEpoch(2000, 1, 1, 12, 0, 0.0)
        ut1 = UT1Epoch(2000, 1, 1, 12, 0, 0.0)
        tcg = TCGEpoch(2000, 1, 1, 12, 0, 0.0)
        tt = TTEpoch(2000, 1, 1, 12, 0, 0.0)
        tdb = TDBEpoch(2000, 1, 1, 12, 0, 0.0)
        tcb = TCBEpoch(2000, 1, 1, 12, 0, 0.0)
        Δtr(ep) = Epochs.diff_tdb_tt(julian1(ep), julian2(ep))
        dut1(ep) = Epochs.dut1(ep)
        dat(ep) = dut1(ep) - Epochs.leapseconds(julian(ep))

        @test Epochs.tttai(julian1(tt), julian2(tt)) == ERFA.tttai(julian1(tt), julian2(tt))
        @test Epochs.tttai(julian2(tt), julian1(tt)) == ERFA.tttai(julian2(tt), julian1(tt))
        @test Epochs.taitt(julian1(tai), julian2(tai)) == ERFA.taitt(julian1(tai), julian2(tai))
        @test Epochs.taitt(julian2(tai), julian1(tai)) == ERFA.taitt(julian2(tai), julian1(tai))

        @test Epochs.ut1tai(julian1(ut1), julian2(ut1), dat(ut1)) == ERFA.ut1tai(julian1(ut1), julian2(ut1), dat(ut1))
        @test Epochs.ut1tai(julian2(ut1), julian1(ut1), dat(ut1)) == ERFA.ut1tai(julian2(ut1), julian1(ut1), dat(ut1))
        @test Epochs.taiut1(julian1(tai), julian2(tai), dat(tai)) == ERFA.taiut1(julian1(tai), julian2(tai), dat(tai))
        @test Epochs.taiut1(julian2(tai), julian1(tai), dat(tai)) == ERFA.taiut1(julian2(tai), julian1(tai), dat(tai))

        @test Epochs.tcgtt(julian1(tcg), julian2(tcg)) == ERFA.tcgtt(julian1(tcg), julian2(tcg))
        @test Epochs.tcgtt(julian2(tcg), julian1(tcg)) == ERFA.tcgtt(julian2(tcg), julian1(tcg))
        @test Epochs.tttcg(julian1(tt), julian2(tt)) == ERFA.tttcg(julian1(tt), julian2(tt))
        @test Epochs.tttcg(julian2(tt), julian1(tt)) == ERFA.tttcg(julian2(tt), julian1(tt))

        @test Epochs.taiutc(julian1(tai), julian2(tai)) == ERFA.taiutc(julian1(tai), julian2(tai))
        @test Epochs.taiutc(julian2(tai), julian1(tai)) == ERFA.taiutc(julian2(tai), julian1(tai))
        @test Epochs.utctai(julian1(utc), julian2(utc)) == ERFA.utctai(julian1(utc), julian2(utc))
        @test Epochs.utctai(julian2(utc), julian1(utc)) == ERFA.utctai(julian2(utc), julian1(utc))

        leap = UTCEpoch(2016, 12, 31, 23, 59, 60)
        tai1, tai2 = ERFA.utctai(julian1(leap), julian2(leap))
        @test Epochs.utctai(julian1(leap), julian2(leap)) == ERFA.utctai(julian1(leap), julian2(leap))
        let (jd2, jd1) = Epochs.taiutc(tai1, tai2)
            erfa_jd2, erfa_jd1 =  ERFA.taiutc(tai1, tai2)
            @test jd2 ≈ erfa_jd2
            @test jd1 ≈ erfa_jd1
        end
        @test Epochs.diff_tdb_tt(julian1(tdb), julian2(tdb), 1.0, 2.0, 3.0, 4.0) == ERFA.dtdb(julian1(tdb), julian2(tdb), 1.0, 2.0, 3.0, 4.0)
        @test Epochs.tdbtt(julian1(tdb), julian2(tdb), Δtr(tdb)) == ERFA.tdbtt(julian1(tdb), julian2(tdb), Δtr(tdb))
        @test Epochs.tdbtt(julian2(tdb), julian1(tdb), Δtr(tdb)) == ERFA.tdbtt(julian2(tdb), julian1(tdb), Δtr(tdb))
        @test Epochs.tttdb(julian1(tt), julian2(tt), Δtr(tdb)) == ERFA.tttdb(julian1(tt), julian2(tt), Δtr(tdb))
        @test Epochs.tttdb(julian2(tt), julian1(tt), Δtr(tdb)) == ERFA.tttdb(julian2(tt), julian1(tt), Δtr(tdb))

        dt(ep) = Epochs.diff_ut1_tt(ep)
        @test Epochs.ttut1(julian1(tt), julian2(tt), dt(tt)) == ERFA.ttut1(julian1(tt), julian2(tt), dt(tt))
        @test Epochs.ttut1(julian2(tt), julian1(tt), dt(tt)) == ERFA.ttut1(julian2(tt), julian1(tt), dt(tt))
        @test Epochs.ut1tt(julian1(ut1), julian2(ut1), dt(ut1)) == ERFA.tttdb(julian1(ut1), julian2(ut1), dt(ut1))
        @test Epochs.ut1tt(julian2(ut1), julian1(ut1), dt(ut1)) == ERFA.tttdb(julian2(ut1), julian1(ut1), dt(ut1))

        # Doing approximate checking due to small machine epsilon. (fails on windows 32-bit)
        let (jd1, jd2) = Epochs.tdbtcb(julian1(tdb), julian2(tdb))
            erfa_jd1, erfa_jd2 = ERFA.tdbtcb(julian1(tdb), julian2(tdb))
            @test jd1 ≈ erfa_jd1
            @test jd2 == erfa_jd2
        end

        let (jd2, jd1) = Epochs.tdbtcb(julian2(tdb), julian1(tdb))
            erfa_jd2, erfa_jd1 = ERFA.tdbtcb(julian2(tdb), julian1(tdb))
            @test jd2 ≈ erfa_jd2
            @test jd1 == erfa_jd1
        end

        @test Epochs.tcbtdb(julian1(tcb), julian2(tcb)) == ERFA.tcbtdb(julian1(tcb), julian2(tcb))
        @test Epochs.tcbtdb(julian2(tcb), julian1(tcb)) == ERFA.tcbtdb(julian2(tcb), julian1(tcb))

        for jd in 2414105.0:10.0:2488985.0
            @test Epochs.diff_tdb_tt(jd, 0.5) ≈ Epochs.diff_tdb_tt(jd,0.5,0.0,0.0,0.0,0.0) atol=40e-6
        end

        @test Epochs.cal2jd(2000, 1, 1) == ERFA.cal2jd(2000, 1, 1)
        @test Epochs.cal2jd(2016, 2, 29) == ERFA.cal2jd(2016, 2, 29)
        @test TTEpoch(Epochs.cal2jd(2000, 1, 1)...) == TTEpoch(2000, 1, 1)
        @test_throws ArgumentError Epochs.cal2jd(-4800, 1, 1)
        @test_throws ArgumentError Epochs.cal2jd(2000, 15, 1)
        @test_throws ArgumentError Epochs.cal2jd(2000, 1, 40)

        @test Epochs.jd2cal(julian1(tt), julian2(tt)) == ERFA.jd2cal(julian1(tt), julian2(tt))
        @test Epochs.jd2cal(julian2(tt), julian1(tt)) == ERFA.jd2cal(julian2(tt), julian1(tt))

        @test Epochs.d2tf(1, -1.7) == ERFA.d2tf(1, -1.7)
        @test Epochs.d2tf(-1, 1.7) == ERFA.d2tf(-1, 1.7)

        @test Epochs.utcut1(julian1(utc), julian2(utc), dut1(utc),leapseconds(julian(utc))) == ERFA.utcut1(julian1(utc), julian2(utc), dut1(utc))
        @test Epochs.utcut1(julian2(utc), julian1(utc),dut1(utc),leapseconds(julian(utc))) == ERFA.utcut1(julian2(utc), julian1(utc), dut1(utc))
        ut1_nearleap = UT1Epoch(2016, 12, 31, 23, 59, 59)
        @test Epochs.ut1utc(julian1(ut1), julian2(ut1), dut1(ut1)) == ERFA.ut1utc(julian1(ut1), julian2(ut1), dut1(ut1))
        @test Epochs.ut1utc(julian2(ut1), julian1(ut1), dut1(ut1)) == ERFA.ut1utc(julian2(ut1), julian1(ut1), dut1(ut1))
        let (jd1, jd2) = Epochs.ut1utc(julian1(ut1_nearleap), julian2(ut1_nearleap), dut1(ut1_nearleap))
            erfa_jd1, erfa_jd2 = ERFA.ut1utc(julian1(ut1_nearleap), julian2(ut1_nearleap), dut1(ut1_nearleap))
            @test jd1 == erfa_jd1
            @test jd2 ≈ erfa_jd2
        end
        let (jd2, jd1) = Epochs.ut1utc(julian2(ut1_nearleap), julian1(ut1_nearleap), dut1(ut1_nearleap))
            erfa_jd2, erfa_jd1 = ERFA.ut1utc(julian2(ut1_nearleap), julian1(ut1_nearleap), dut1(ut1_nearleap))
            @test jd2 ≈ erfa_jd2
            @test jd1 == erfa_jd1
        end

        @test Epochs.datetime2julian(UTC, 2016, 12, 31, 23, 59, 60) == ERFA.dtf2d("UTC", 2016, 12, 31, 23, 59, 60)
        @test Epochs.datetime2julian(TT, 2016, 12, 31, 23, 59, 59) == ERFA.dtf2d("TT", 2016, 12, 31, 23, 59, 59)
        @test_throws ArgumentError Epochs.datetime2julian(TT, 2016, 12, 31, 23, 59, 60)
        @test_throws ArgumentError Epochs.datetime2julian(UTC, 2016, 12, 31, 25, 59, 60)
        @test_throws ArgumentError Epochs.datetime2julian(UTC, 2016, 12, 31, 23, 61, 60)
        @test_throws ArgumentError Epochs.datetime2julian(UTC, 2016, 12, 31, 23, 59, 61)

        utc2 = UTCEpoch(2001, 03, 01, 23, 30, 0)
        @test Epochs.julian2datetime(timescale(leap), 3, leap.jd1, leap.jd2) == ERFA.d2dtf("UTC", 3, leap.jd1, leap.jd2)
        @test Epochs.julian2datetime(timescale(tt), 3, tt.jd1, tt.jd2) == ERFA.d2dtf("TT", 3, tt.jd1, tt.jd2)
        @test Epochs.julian2datetime(timescale(utc2), 3, utc2.jd1, utc2.jd2) == ERFA.d2dtf("UTC", 3, utc2.jd1, utc2.jd2)
    end
    @testset "Leap Seconds" begin
        @test leapseconds(TTEpoch(1959,1,1)) == 0
        @test fractionofday(DateTime(1959,1,1)) == 0.0
        # Doing approximate checking due to small machine epsilon. (fails on windows 32-bit)
        for year = 1960:Dates.year(now())
            @test leapseconds(TTEpoch(year, 4, 1)) ≈ ERFA.dat(year, 4, 1, 0.0)
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.19.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 1f0cb625-d1ec-4775-bc75-59d4675d181e
# ╠═╡ show_logs = false
begin
	using Logging
	global_logger(NullLogger())
	display("")
end

# ╔═╡ 25448e4a-345b-4bd6-9cdb-a6a280cfd22c
# ╠═╡ show_logs = false
#Set-up packages
begin

	using PlutoUI, DataFrames, HTTP, CSV, Dates, Printf, LaTeXStrings, HypertextLiteral, XLSX

	# gr();
	# Plots.GRBackend()

	#Define html elements
	nbsp = html"&nbsp" #non-breaking space
	vspace = html"""

"""
	br = html"
"

	#Sets the width of cells, caps the cell width by 90% of screen width
	#(setting overwritten by cell below)
	# @bind screenWidth @htl("""
	#

	#
	#

	# """)

	# cellWidth= min(1000, screenWidth*0.9)
	# @htl("""
	#
	# """)

	#Sets the width of the cells
	#begin
	#	html"""


text:
using Base.Test
using DataFrames
using Celeste.SDSSIO: PlainFITSStrategy
import WCS

const rcf = RunCamcolField(3900, 6, 269)

"""
Get the PSF located at a particular world location in an image.

Args:
 - world_loc: A location in world coordinates.
 - img: An Image

Returns:
 - An array of PsfComponent objects that represents the PSF as a mixture
     of Gaussians.
"""
function get_source_psf(world_loc::Vector{Float64}, img::Image, psf_K::Int)
    # Some stamps or simulated data have no raw psf information. In that case,
    # just use the psf from the image.
    if size(img.raw_psf_comp.rrows) == (0, 0)
        # Also return a vector of empty psf params
        return img.psf, fill(fill(NaN, length(PsfParams)), psf_K)
    else
        pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
        psfstamp = Model.eval_psf(img.raw_psf_comp, pixel_loc[1], pixel_loc[2])
        return PSF.fit_raw_psf_for_celeste(psfstamp, psf_K)
    end
end

function test_images()
    # A lot of tests are in a single function to avoid having to reload
    # the full image multiple times.
    strategy = PlainFITSStrategy(datadir)
    images = SDSSIO.load_field_images(strategy, [rcf])

    cat_entries = convert(Vector{CatalogEntry}, SDSSIO.read_photoobj(strategy, rcf))

    ea = make_elbo_args(images, cat_entries, patch_radius_pix=1e-6)

    # Just check some basic facts about the catalog.
    @test length(cat_entries) > 500
    @test 0 < sum([ce.is_star for ce in cat_entries]) < 805

    # Find an object near the middle of the image.
    ctr = WCS.pix_to_world(images[3].wcs, [images[3].H / 2, images[3].W / 2])
    dist = [sqrt((ce.pos[1] - ctr[1])^2 + (ce.pos[2] - ctr[2])^2)
                    for ce in cat_entries]
    obj_index = findmin(dist)[2]    # index of closest object
    obj_loc = cat_entries[obj_index].pos    # location of closest object

    # Test get_source_psf at point while we have the images loaded.
    test_b = 3
    img = ea.images[test_b]
    ea_obj = make_elbo_args(images, cat_entries[obj_index:obj_index])
    vs_obj = DeterministicVI.catalog_init_source(cat_entries[obj_index])
    pixel_loc = WCS.world_to_pix(img.wcs, obj_loc)
    original_psf_val = Model.eval_psf(img.raw_psf_comp, pixel_loc[1], pixel_loc[2])

    original_psf_celeste =
        PSF.fit_raw_psf_for_celeste(original_psf_val, ea.psf_K)[1]
    fit_original_psf_val = PSF.get_psf_at_point(original_psf_celeste)

    obj_psf = get_source_psf(vs_obj[ids.pos], img, ea.psf_K)[1]
    obj_psf_val = PSF.get_psf_at_point(obj_psf)

    # The fits should match exactly.
    @test isapprox(obj_psf_val, fit_original_psf_val, atol=1e-6)

    # The raw psf will not be as good.
    @test isapprox(obj_psf_val, original_psf_val, atol=1e-2)

    cat_several = [cat_entries[1]; cat_entries[obj_index]]
    ea_several = make_elbo_args(ea_obj.images, cat_several)

    # The second set of vp is the object of interest
    point_patch_psf = PSF.get_psf_at_point(ea_several.patches[2, test_b].psf)
    # The threshold for the test below was formerly 1e-6.
    # Is it a problem I needed to increase it?
    @test isapprox(obj_psf_val, point_patch_psf, atol=5e-4)
end

test_images()


text:
"""
USE THIS AS A TEMPLATE WHEN GENERATING NEW PACKAGE!!
"""
using PkgTemplates

t = Template(;
     user = "Shoram444",
     license = "MIT",
     authors = ["Maros Petro"],
     plugins = [
         TravisCI(),
         Codecov(),
         Coveralls(),
         AppVeyor(),
     ],
)

generate("PackageName", t)


text:
function JuliaBenchSIMD( operationMode )

    allFunctions = [MatrixAddition, MatrixMultiplication, ElementWiseOperations]; # only SIMD functions to run

    allFunctionsString = [ "Matrix Addition", "Matrix Multiplication", "Element Wise Operations"];

	if (operationMode == 1) # partial benchmark
		vMatrixSize =  dropdims(readdlm(joinpath("Inputs","vMatrixSizePartial.csv"), ',',Int64), dims=1);
		numIterations = dropdims(readdlm(joinpath("Inputs","numIterationsPartial.csv"), ',',Int64), dims=1);

	elseif (operationMode == 2) # full benchmark
		vMatrixSize = dropdims(readdlm(joinpath("Inputs","vMatrixSizeFull.csv"), ',',Int64), dims=1);
		numIterations =  dropdims(readdlm(joinpath("Inputs","numIterationsFull.csv"), ',',Int64), dims=1);

	elseif (operationMode == 0) # Test benchmark
		vMatrixSize = 2;
		numIterations =  1;

	end

    numIterations = numIterations[1]; # It is 1x1 Array -> Scalar

    mRunTime = zeros(length(vMatrixSize), length(allFunctions), numIterations);
    tRunTime= Array{Any}(undef,length(allFunctions)+1,length(vMatrixSize)+1)# a table containing all the information
	tRunTime[1,1]="FunctionName\\MatrixSize";

    for ii = 1:length(vMatrixSize)
        matrixSize = vMatrixSize[ii];
        mX = randn(matrixSize, matrixSize);
        mY = randn(matrixSize, matrixSize);
        println("Matrix Size - $matrixSize");

        jj=1;
        for fun in allFunctions
            println("Processing $(allFunctionsString[jj]) - MatrixSize= $matrixSize");

            for kk = 1:numIterations;

                benchIJK =@benchmark $fun($matrixSize, $mX, $mY)
                # t =@benchmarkable $fun($matrixSize, $mX, $mY);
                # tune!(t)
                # run(t)

                mRunTime[ii, jj, kk]=median(benchIJK).time/1e3;
                # println("$(mRunTime[ii, jj, kk])")

            end
			tRunTime[jj+1,1]="$(allFunctionsString[jj])";
			tRunTime[1,ii+1]="$matrixSize";
			tRunTime[jj+1,ii+1]=mean(mRunTime[ii, jj,:]);
            jj+=1;
        end

    end

    return tRunTime, mRunTime;
end

#=
function MatrixGeneration( matrixSize, mX, mY )

  mA = randn(matrixSize, matrixSize);
  mB = rand(matrixSize, matrixSize);

  return mA;
end
=#
function MatrixAddition( matrixSize, mX, mY )

  scalarA = rand();
  scalarB = rand();

  # mA = (scalarA .* mX) .+ (scalarB .* mY);
  mA = Array{Float64}(undef, matrixSize, matrixSize);
  @simd for ii = 1:(matrixSize * matrixSize)
    @inbounds mA[ii] = (scalarA * mX[ii]) + (scalarB * mY[ii]);
  end

  return mA;
end

function MatrixMultiplication( matrixSize, mX, mY )

  scalarA = rand();
  scalarB = rand();

  # mA = (scalarA .+ mX) * (scalarB .+ mY);
  mA = Array{Float64}(undef, matrixSize, matrixSize);
  @simd for ii = 1:(matrixSize * matrixSize)
    @inbounds mA[ii] = (scalarA + mX[ii]) * (scalarB + mY[ii]);
  end

  return mA;
end
#=
function MatrixQuadraticForm( matrixSize, mX, mY )

  vX = randn(matrixSize);
  vB = randn(matrixSize);
  sacalrC = rand();

  mA = (transpose(mX * vX) * (mX * vX)) .+ (transpose(vB) * vX) .+ sacalrC;

  return mA;
end

function MatrixReductions( matrixSize, mX, mY )

  mA = sum(mX, dims=1) .+ minimum(mY, dims=2); #Broadcasting

  return mA;
end
=#
function ElementWiseOperations( matrixSize, mX, mY )

  mA = rand(matrixSize, matrixSize);
  mB = 3 .+ rand(matrixSize, matrixSize);
  mC = rand(matrixSize, matrixSize);

  # mD = abs.(mA) .+ sin.(mA);
  mD = Array{Float64}(undef, matrixSize, matrixSize);
  @simd for ii = 1:(matrixSize * matrixSize)
    @inbounds mD[ii] = abs(mA[ii]) + sin(mA[ii]);
  end

  # mE = exp.(-(mA .^ 2));
  mE = Array{Float64}(undef, matrixSize, matrixSize);
  @simd for ii = 1:(matrixSize * matrixSize)
    @inbounds mE[ii] = exp(- (mA[ii] * mA[ii]));
  end

  # mF = (-mB .+ sqrt.((mB .^ 2) .- (4 .* mA .* mC))) ./ (2 .* mA);
  mF = Array{Float64}(undef, matrixSize, matrixSize);
  @simd for ii = 1:(matrixSize * matrixSize)
    @inbounds mF[ii] = (-mB[ii] + sqrt( (mB[ii] * mB[ii]) - (4 * mA[ii] * mC[ii]) )) ./ (2 * mA[ii]);
  end

  mA = mD .+ mE .+ mF;

  return mA;
end

#=
function MatrixExp( matrixSize, mX, mY )

    mA = exp(mX);

    return mA;
end

function MatrixSqrt( matrixSize, mX, mY )

    mY = transpose(mX) * mX;

    mA = sqrt(mY);

  return mA;
end

function Svd( matrixSize, mX, mY )

    F = svd(mX, full = false); # F is SVD object
    mU, mS, mV = F;

    return mA=0;
end

function Eig( matrixSize, mX, mY )

    F  = eigen(mX); # F is eigen object
    mD, mV = F;

    return mA=0;
end

function CholDec( matrixSize, mX, mY )

  mY = transpose(mX) * mX;

  mA = cholesky(mY);

  return mA;
end

function MatInv( matrixSize, mX, mY )

  mY = transpose(mX) * mX;

  mA = inv(mY);
  mB = pinv(mX);

  mA = mA .+ mB;

  return mA;
end

function LinearSystem( matrixSize, mX, mY )

  mB = randn(matrixSize, matrixSize);
  vB = randn(matrixSize);

  vA = mX \ vB;
  mA = mX \ mB;

  mA = mA .+ vA;

  return mA;
end

function LeastSquares( matrixSize, mX, mY )

  mB = randn(matrixSize, matrixSize);
  vB = randn(matrixSize);

  vA = (transpose(mX) * mX) \ (transpose(mX) * vB);
  mA = (transpose(mX) * mX) \ (transpose(mX) * mB);

  mA = mA .+ vA;

  return mA;
end

function CalcDistanceMatrix( matrixSize, mX, mY )

  mY = randn(matrixSize, matrixSize);

  mA = transpose(sum(mX .^ 2, dims=1)) .- (2 .* transpose(mX) * mY) .+ sum(mY .^ 2, dims=1);

  return mA;
end

function KMeans( matrixSize, mX, mY )

    # Assuming Samples are slong Columns (Rows are features)
    numClusters  = Int64( max( round(matrixSize / 100), 1 ) ); # % max between 1 and round(...)
    numIterations = 10;

    # http://stackoverflow.com/questions/36047516/julia-generating-unique-random-integer-array
    mA = mX[:, randperm(matrixSize)[1:numClusters]]; #<! Cluster Centroids

    for ii = 1:numIterations
        vMinDist, mClusterId = findmin( transpose(sum(mA .^ 2, dims=1)) .- (2 .* transpose(mA)* mX), dims=1); #<! Is there a `~` equivalent in Julia?
        vClusterId = LinearIndices( dropdims(mClusterId, dims=1) ); # to be able to access it later

        for jj = 1:numClusters
            mA[:, jj] = mean( mX[:, vClusterId .== jj ], dims=2 );
        end
    end

    mA = mA[:, 1] .+ transpose(mA[:, end]);

  return mA;
end
=#

max_stars_repo_name: aminya/Julia-Matlab-Benchmark

max_issues_repo_name: aminya/Julia-Matlab-Benchmark

max_forks_repo_name: aminya/Julia-Matlab-Benchmark


text:
# these are used internally
const double_π = DoubleFloat64(pi)
const double_2pi = convert(FastDouble, 2 * BigFloat(Base.pi))
const double_pi = convert(FastDouble, Base.pi)
const double_pi2 = convert(FastDouble, BigFloat(Base.pi) * 0.5)
const double_pi4 = convert(FastDouble, BigFloat(Base.pi) * 0.25)
const double_pi16 = convert(FastDouble, BigFloat(Base.pi) * (1/16))
const double_3pi4 = convert(FastDouble, BigFloat(Base.pi) * 0.75)
const double_nan = DoubleFloat64{ComputeFast}(NaN, NaN)
const double_inf = DoubleFloat64{ComputeFast}(Inf)
const double_e = convert(FastDouble, Base.e)

const double_log2 = convert(FastDouble, log(BigFloat(2.0)))
const double_log10 = convert(FastDouble, log(BigFloat(10.0)))

const double_eps = 4.93038065763132e-32 # 2^-104


text:
# LI7200Diagnostic.jl
#
#	Parse LI-7200 diagnostic values
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
#
# Julia 1.6.2
# 17.12.2020
# Last Edit: 13.08.2021

"""# LI7200Diagnostic

Parse the LI-7200 diagnostic values

---

### Example

```jldoctest
julia> x = LI7200Diagnostic(8190)
LI-7200 Diagnistic Values (8190)
		AGC (dirtiness): 87.5% - Optics need cleaning
		Sync OK
		PLL OK - Optical wheel rotating correctly
		Detector Temperature OK
		Chopper OK - Wheel temperature OK
		Pressure OK - Differential pressure sensor OK
		Voltages OK - Internal reference voltages OK
		Tinlet OK - Inlet temperature OK
		Toutlet OK - Outlet temperature OK
		Head detected
```

```jldoctest
julia> x.value
8190
```

```jldoctest
julia> x.TinletOK
true
```

---

`diag = LI7200Diagnostic(x);`\n
* **x**::Int or Float64 = Unparsed diagnostic value
* **diag**::LI7200Diagnostic = Parsed diagnostic value

---

## Fields
* **value**::Int
* **AGC**::Int
* **SyncError**::Bool
* **PLL_OK**::Bool
* **DetectorOK**::Bool
* **ChopperWheelOK**::Bool
* **PressureOK**::Bool
* **VoltageOK**::Bool
* **TinletOK**::Bool
* **ToutletOK**::Bool
* **HeadDetected**::Bool
* **AGC**::Int
"""
struct LI7200Diagnostic
	"Original diagnostic value"
	value::Int
	"AGC - Automatic Gain Control (dirtiness)"
	AGC::Int
	"Sync Error"
	SyncError::Bool
	"Optical wheel rotation"
	PLL_OK::Bool
	"Detector Temperature"
	DetectorOK::Bool
	"Chopper wheel"
	ChopperWheelOK::Bool
	"Differential pressure sensor"
	PressureOK::Bool
	"Internal reference voltage"
	VoltageOK::Bool
	"Inlet Temperature"
	TinletOK::Bool
	"Outlet Temperature"
	ToutletOK::Bool
	"Sensor Head Detection"
	HeadDetected::Bool

	function LI7200Diagnostic(val::Float64)
		LI7200Diagnostic(Int(val))
	end

	function LI7200Diagnostic(val::Int)
		value = val

		AGC = Int(UInt16(val & 2^(0) > 0) << 0 | UInt16(val & 2^(1) > 0) << 1 | UInt16(val & 2^(2) > 0) << 2 | UInt16(val & 2^(3) > 0) << 3)
		SyncError = val & 2^(4) > 0
		PLL_OK = val & 2^(5) > 0
		DetectorOK = val & 2^(6) > 0
		ChopperWheelOK = val & 2^(7) > 0
		PressureOK = val & 2^(8) > 0
		VoltageOK = val & 2^(9) > 0
		TinletOK = val & 2^(10) > 0
		ToutletOK = val & 2^(11) > 0
		HeadDetected = val & 2^(12) > 0

		new(value ,AGC, SyncError, PLL_OK, DetectorOK, ChopperWheelOK, PressureOK, VoltageOK, TinletOK, ToutletOK, HeadDetected)
	end # End of constructor
end # End of type

function Base.show(io::IO,::MIME"text/plain", diag::LI7200Diagnostic)
	# AGC and Errors only
	print(string(diag.value),
		" - AGC(" * string(diag.AGC) * "): ",diag.AGC < 14 ? @sprintf("%2.0f",diag.AGC * 6.25) * "%" : RED_FG(@sprintf("%2.0f",diag.AGC * 6.25) * "% - Optics need cleaning"))
	diag.SyncError ? nothing : print(", Sync: ",RED_FG("Error"))
	diag.PLL_OK ? nothing : print(", PLL: ",RED_FG("Error"))
	diag.DetectorOK ? nothing : print(", Detector: ",RED_FG("Error"))
	diag.ChopperWheelOK ? nothing : print(", Chopper: ",RED_FG("Error"))
	diag.PressureOK ? nothing : print(", Pressure: ",RED_FG("Error"))
	diag.VoltageOK ? nothing : print(", Voltage: ",RED_FG("Error"))
	diag.TinletOK ? nothing : print(", Tinlet: ",RED_FG("Error"))
	diag.ToutletOK ? nothing : print(", Toutlet: ",RED_FG("Error"))
	diag.HeadDetected ? nothing : print(", Head: ",RED_FG("Absent"))

	# Full Report Version
	#=println(string(diag.value),
		" - AGC: ",diag.AGC < 14 ? @sprintf("%2.0f",diag.AGC * 6.25) * "%" : RED_FG(@sprintf("%2.0f",diag.AGC * 6.25) * "% - Optics need cleaning"),
		", Sync: ",diag.SyncError ? GREEN_FG("+") : RED_FG("Error"),
		", PLL: ",diag.PLL_OK ? GREEN_FG("+") : RED_FG("Error"),
		", Detector: ",diag.DetectorOK ? GREEN_FG("+") : RED_FG("Error"),
		", Chopper: ",diag.ChopperWheelOK ? GREEN_FG("+") : RED_FG("Error"),
		", Pressure: ",diag.PressureOK ? GREEN_FG("+") : RED_FG("Error"),
		", Volt: ",diag.VoltageOK ? GREEN_FG("+") : RED_FG("Error"),
		", Tin:",diag.TinletOK ? GREEN_FG("+") : RED_FG("Error"),
		", Tout: ",diag.ToutletOK ? GREEN_FG("+") : RED_FG("Error"),
		", Head: ",diag.HeadDetected ? GREEN_FG("Present") : RED_FG("Absent"))=#
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
using HTTP
using JSON

function weather_past_week(api_key)
    req_url = string("https://api.nasa.gov/insight_weather/?api_key=", api_key, "?&feedtype=json&ver=1.0")
    req = HTTP.request("GET", req_url)
    json = JSON.parse(String(req.body))
    return json
end


text:
#
# This script documents the cost of mapping from a structure
# to an environment, which doubles the amount of memory used
# (but this can of course be allocated), so primarily it is
# the cost of copying the memory.
#
using ACE, StaticArrays, BenchmarkTools
using ACE: State

randstate() =
           ( rr = randn(SVector{3, Float64}),
             μ = rand([1,2,3,4,5]),
             q = randn(),
             p = randn(SVector{3, Float64}) )

function struct2env(Xs, i, J)
   xi = Xs[i]
   return [ ( xj = Xs[j];
                   (rr = xj.rr - xi.rr,
                    μi = xi.μ, μj = xj.μ,
                    qi = xi.q, qj = xj.q,
                    pi = xi.p, pj = xj.p) )  for j in J ]
end

function struct2env!(Xs_i, Xs, i, J)
   xi = Xs[i]
   for (ij, j) in enumerate(J)
      xj = Xs[j]
      Xs_i[ij] =      (rr = xj.rr - xi.rr,
                       μi = xi.μ, μj = xj.μ,
                       qi = xi.q, qj = xj.q,
                       pi = xi.p, pj = xj.p )
   end
   return Xs_i
end

Xs = [ randstate() for _ = 1:1000 ]
i = rand(1:1000)
J = setdiff(unique(rand(1:1000, 35)), [i])

Xs_i = struct2env(Xs, i, J)
struct2env!(Xs_i, Xs, i, J)

@btime struct2env($Xs, $i, $J)
@btime struct2env!($Xs_i, $Xs, $i, $J)

##

struct Wrapper{SYMS, TT}
   x::NamedTuple{SYMS, TT}
end
randnt() = Wrapper((a = randn(), b = rand(), c = rand(1:5),
                   d = randn(SVector{3, Float64})))
X = randnt()
@allocated randnt()
isbits(X)
isbits(X.x)

struct Wrapper1{NT}
   x::NT
end
randnt1() = Wrapper1((a = randn(), b = rand(), c = rand(1:5), d = randn(SVector{3, Float64})))
X1 = randnt1()
@allocated randnt1()
isbits(X1)
isbits(X1.x)

##


text:
"""
    RexHopper <: Environment

    hopping robot designed and build by the Robotic Exploration Lab (Benjamin Bokser)
"""
struct RexHopper end

function rexhopper(;
    representation=:minimal,
    model=:rexhopper,
    timestep=0.05,
    seed=1,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=1.0,
    spring=0.0,
    damper=1.0,
    contact_type=:nonlinear,
    contact_foot=true,
    contact_body=true,
    limits=true,
    info=nothing,
    infeasible_control=false,
    vis=Visualizer(),
    name=:robot,
    opts_step=SolverOptions(rtol=1.0e-4, btol=1.0e-4, undercut=10.0),
    opts_grad=SolverOptions(rtol=1.0e-4, btol=1.0e-4, undercut=10.0),
    T=Float64)

    mechanism = get_rexhopper(
        model=model,
        limits=limits,
        contact_type=contact_type,
        timestep=timestep,
        gravity=gravity,
        friction_coefficient=friction_coefficient,
        spring=spring,
        damper=damper,
        contact_foot=contact_foot,
        contact_body=contact_body)

    initialize_rexhopper!(mechanism)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end
    nu_inf = input_dimension(mechanism)
    nu = infeasible_control ? nu_inf : 5
    no = nx

    # values taken from Mujoco's model, combining the control range -1, 1 and the motor gears.
    aspace = BoxSpace(nu,
        low=(-Inf * ones(nu)),
        high=(Inf * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)

    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z

    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_map = infeasible_control ? 1.0 * I(nu) : cat(zeros(3, 3), 1.0 * I(3), 1.0, 0.0, 1.0, zeros(5, 5), dims=(1,2))

    build_robot(mechanism, vis=vis, name=name)

    TYPES = [RexHopper, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    env = Environment{TYPES...}(mechanism,
        representation,
        aspace,
        ospace,
        x, fx, fu,
        u_prev, control_map,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)

    return env
end

function Base.reset(env::Environment{RexHopper};
    x=nothing)

    if x != nothing
        env.state .= x
    else
        # initialize above the ground to make sure that with random initialization we do not violate the ground constraint.
        initialize!(env.mechanism, :rexhopper)
        x = get_minimal_state(env.mechanism)
        nx = minimal_dimension(env.mechanism)
        z = minimal_to_maximal(env.mechanism, x)
        set_maximal_state!(env.mechanism, z)
        if env.representation == :minimal
            env.state .= get_minimal_state(env.mechanism)
        elseif env.representation == :maximal
            env.state .= get_maximal_state(env.mechanism)
        end
        env.input_previous .= 0.0
    end
    return get_observation(env)
end


text:
module SALTBase

using LinearAlgebra
using CatViews
using ForwardDiff: derivative  # used in gain.jl
using AbbreviatedTypes
using Printf: @printf

export AbstractMaxwellSolver

# Each concrete SALT solver package (e.g., MaxwellSALT) should define a concrete subtype of
# AbstractMaxwellSolver to store the information (e.g., the Maxwell operator A) needed by the
# specific linear Maxwell solver to use in the SALT solver package.
abstract type AbstractMaxwellSolver end

# Each concrete SALT solver package (e.g., MaxwellSALT) should extend the following Base
# functions:
#
# - similar(mxsolver:AbstractMaxwellSolver)
# - size(mxsolver::AbstractMaxwellSolver)

# Each concrete SALT solver package (e.g., MaxwellSALT) should implement a concrete subtype of
# AbstractMaxwellSolver and extend the following functions to be used in lasing.jl and nonlasing.jl:
#
# init_mxsolver!(mxsolver::AbstractMaxwellSolver, ω::Number, ε::AbsVecNumber)
# mxsolve!(x::AbsVecComplexF, mxsolver::AbstractMaxwellSolver, b::AbsVecNumber)
# mxsolve_transpose!(x::AbsVecComplexF, mxsolver::AbstractMaxwellSolver, b::AbsVecNumber)
# mxapply!(b::AbsVecComplexF, mxsolver::AbstractMaxwellSolver, x::AbsVecNumber)
# abs2_interp!(abs2ψ::AbsVecReal, mxsolver::AbstractMaxwellSolver, ψ::AbsVecNumber)
function init_mxsolver! end  # initialize mxsolver for current ω and ε
function mxsolve! end  # x = A⁻¹ b
function mxsolve_transpose! end  # x = A⁻ᵀ b
function mxapply! end  # b = A x
function abs2_interp! end  # abs2ψ = abs2.(interpolate(ψ))

include("util.jl")
include("gain.jl")
include("lasing.jl")
include("nonlasing.jl")
include("anderson.jl")
include("switching.jl")
include("simulation.jl")

end # module


text:
module Inference

using Gen
using ..MangroveBase: MangroveData, Rating, Sub, mean, subs, kids
import Gen.get_traces

function get_traces(
    model::DynamicDSLFunction,
    ratings::MangroveData,
    amount_of_computation::Int
  )::Vector{Gen.Trace}
  println("Running inference...")
  observations = Gen.choicemap()
  subjects = subs(ratings)
  reviewers = kids(ratings)
  for sub in subjects
    for kid in reviewers
      observations[:selected => sub => kid] = false
    end
  end
  for (k, v) in ratings
    observations[k => :rating] = v
    observations[:selected => k.sub => k.kid] = true
  end

  trace, = Gen.importance_resampling(
    model,
    (subjects, reviewers),
    observations,
    100
  )

  traces = Gen.Trace[]

  for i in 1:amount_of_computation
    if i%(amount_of_computation / 20) == 0
      println(
        "MH iteration ", i, " out of ", amount_of_computation,
        " (", round(Int, i/amount_of_computation*100), "%)"
      )
    end
    trace, = mh(trace, select(:quality))
    push!(traces, trace)
  end

  return traces
end

function qualities(t::Gen.Trace)::Dict{Sub, Rating}
  Dict(sub => t[:quality => sub] for sub in get_args(t)[1])
end

function mean_l2(data::MangroveData, traces::Vector{Gen.Trace})::Real
    subjects = subs(data)
  data_mean = mean(data)
  trace_qualities = (qualities(trace) for trace in traces)
  sum((data_mean[sub] - mean(qualities[sub] for qualities in trace_qualities))^2 for sub in subjects) / length(subjects)
end

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
# capsule contact point origin sphere-capsule
function contact_point_segment(ps, xc, qc, ka, kb)
    # capsule origin a
    ca = contact_point_origin(xc, qc, ka)

    # capsule origin b
    cb = contact_point_origin(xc, qc, kb)

    # point differences
    dab = cb - ca
    das = ps - ca

    t = dot(das, dab) / dot(dab, dab)
    # t_clamp = min(max(t, 0.0), 1.0)

    # closest point
    if t <= 0.0
        p = ca
    elseif t >= 1.0
        p = cb
    else
        p = ca + t * dab
    end

    return p
end

# capsule contact point origin sphere-capsule
function ∂contact_point_segment∂p(ps, xc, qc, ka, kb)
    # capsule origin a
    ca = contact_point_origin(xc, qc, ka)

    # capsule origin b
    cb = contact_point_origin(xc, qc, kb)

    # point differences
    dab = cb - ca
    das = ps - ca

    t = dot(das, dab) ./ dot(dab, dab)
    # t_clamp = min(max(t, 0.0), 1.0)

    # closest point
    if t <= 0.0
        # p = ca
        return szeros(3, 3)
    elseif t >= 1.0
        # p = cb
        return szeros(3, 3)
    else
        # p = ca + t * dab
        return dab * dab' ./ dot(dab, dab)
    end

    return p
end

# capsule contact point origin sphere-capsule
function ∂contact_point_segment∂x(ps, xc, qc, ka, kb)
    # capsule origin a
    ca = contact_point_origin(xc, qc, ka)

    # capsule origin b
    cb = contact_point_origin(xc, qc, kb)

    # point differences
    dab = cb - ca
    das = ps - ca

    t = dot(das, dab) ./ dot(dab, dab)
    # t_clamp = min(max(t, 0.0), 1.0)

    # closest point
    if t <= 0.0
        # p = ca
        return ∂contact_point_origin∂x(xc, qc, ka)
    elseif t >= 1.0
        # p = cb
        return ∂contact_point_origin∂x(xc, qc, kb)
    else
        # p = ca + t * dab
        X = ∂contact_point_origin∂x(xc, qc, ka)
        X += t * (∂contact_point_origin∂x(xc, qc, kb) - ∂contact_point_origin∂x(xc, qc, ka))

        t1 = 1.0 / (dab' * dab) * dab' * (-∂contact_point_origin∂x(xc, qc, ka))
        t2 = 1.0 / (dab' * dab) * das' * (∂contact_point_origin∂x(xc, qc, kb) - ∂contact_point_origin∂x(xc, qc, ka))
        t3 = -2.0 * das' * dab / (dab' * dab)^2 * dab' * (∂contact_point_origin∂x(xc, qc, kb) - ∂contact_point_origin∂x(xc, qc, ka))

        X += dab * (t1 + t2 + t3)
        return X
    end
end

# capsule contact point origin sphere-capsule
function ∂contact_point_segment∂q(ps, xc, qc, ka, kb)
    # capsule origin a
    ca = contact_point_origin(xc, qc, ka)

    # capsule origin b
    cb = contact_point_origin(xc, qc, kb)

    # point differences
    dab = cb - ca
    das = ps - ca

    t = dot(das, dab) ./ dot(dab, dab)
    # t_clamp = min(max(t, 0.0), 1.0)

    # closest point
    if t <= 0.0
        # p = ca
        return ∂contact_point_origin∂q(xc, qc, ka)
    elseif t >= 1.0
        # p = cb
        return ∂contact_point_origin∂q(xc, qc, kb)
    else
        # p = ca + t * dab
        Q = ∂contact_point_origin∂q(xc, qc, ka)
        Q += t * (∂contact_point_origin∂q(xc, qc, kb) - ∂contact_point_origin∂q(xc, qc, ka))

        t1 = 1.0 / (dab' * dab) * dab' * (-∂contact_point_origin∂q(xc, qc, ka))
        t2 = 1.0 / (dab' * dab) * das' * (∂contact_point_origin∂q(xc, qc, kb) - ∂contact_point_origin∂q(xc, qc, ka))
        t3 = -2.0 * das' * dab ./ (dab' * dab)^2 * dab' * (∂contact_point_origin∂q(xc, qc, kb) - ∂contact_point_origin∂q(xc, qc, ka))
        Q += dab * (t1 + t2 + t3)

        return Q
    end
end

# # case 1
# ps = [-2.0 - 0.0e-8; 0.0; 1.0]
# xc = [0.0; 0.0; 0.0]
# qc = Quaternion([1.0; 0.0; 0.0; 0.0])
# ka = [1.0; 0.0; 0.0]
# kb = [-1.0; 0.0; 0.0]

# contact_point_segment(ps, xc, qc, ka, kb)

# ∂contact_point_segment∂p(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(p -> contact_point_segment(p, xc, qc, ka, kb), ps)

# ∂contact_point_segment∂x(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(x -> contact_point_segment(ps, x, qc, ka, kb), xc)

# ∂contact_point_segment∂q(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(q -> contact_point_segment(ps, xc, Quaternion(q..., false), ka, kb), vector(qc))

# # case 2
# ps = [1.0 + 100.0e-8; 0.0; 1.0]
# xc = [0.0; 0.0; 0.0]
# qc = Quaternion([1.0; 0.0; 0.0; 0.0])
# ka = [1.0; 0.0; 0.0]
# kb = [-1.0; 0.0; 0.0]

# contact_point_segment(ps, xc, qc, ka, kb)

# ∂contact_point_segment∂p(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(p -> contact_point_segment(p, xc, qc, ka, kb), ps)

# ∂contact_point_segment∂x(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(x -> contact_point_segment(ps, x, qc, ka, kb), xc)

# ∂contact_point_segment∂q(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(q -> contact_point_segment(ps, xc, Quaternion(q..., false), ka, kb), vector(qc))

# # case 3
# ps = [0.0; 0.0; 1.0]
# xc = [0.0; 0.0; 0.0]
# qc = Quaternion([1.0; 0.0; 0.0; 0.0])
# ka = [1.0; 0.0; 0.0]
# kb = [-1.0; 0.0; 0.0]

# contact_point_segment(ps, xc, qc, ka, kb)

# ∂contact_point_segment∂p(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(p -> contact_point_segment(p, xc, qc, ka, kb), ps)

# ∂contact_point_segment∂x(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(x -> contact_point_segment(ps, x, qc, ka, kb), xc)

# ∂contact_point_segment∂q(ps, xc, qc, ka, kb)
# FiniteDiff.finite_difference_jacobian(q -> contact_point_segment(ps, xc, Quaternion(q..., false), ka, kb), vector(qc))

# function get_t(a, b)
#     return a' * b ./ (b' * b)
# end

# function dget_tda(a, b)
#     return b' ./ (b' * b)
# end

# function dget_tdb(a, b)
#     D = a' ./ (b' * b)
#     D += -1.0 * a' * b ./ (b' * b)^2 * 2.0 * b'
#     return D
# end

# a = rand(3)
# b = rand(3)

# get_t(a, b)

# dget_tda(a, b)
# FiniteDiff.finite_difference_jacobian(x -> get_t(x, b), a)

# dget_tdb(a, b)
# FiniteDiff.finite_difference_jacobian(x -> get_t(a, x), b)


text:
@testset "NMU" begin
    nmu = NMU(2,3)
    x = rand(2,10)
    z = nmu(x)

    @test size(z) == (3,10)
    ps = params(nmu)
    @test length(ps) == 1

    loss() = sum(abs2, nmu(x) .- reshape(x[1,:],1,:))
    opt = ADAM()
    train_data = Iterators.repeated((), 10000)
    Flux.train!(loss, ps, train_data, opt)

    @test all(isapprox.(nmu.W, [1.0 0.0; 1.0 0.0; 1.0 0.0], atol=1e-4))
end