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

#=
`convert(C, c)` might be called as `convert(RGB, c)` or
`convert(RGB{Float32}, c)`.
This is handled in ColorTypes, which calls functions
```
    _convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c)
    _convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c, alpha)
```
Here are the argument types:
- `Cdest` may be any concrete `Colorant{T,N}` type. For parametric Color types
  it _always_ has the desired element type (e.g., `Float32`), so it's safe to
  dispatch on `Cdest <: Colorant{T}`.
- `Odest` and `Osrc` are opaque `Color` subtypes, i.e., things like `RGB` or
  `HSV`. They have no element type.
- `c` is the `Colorant` object you wish to convert.
- `alpha`, if present, is a user-supplied alpha value (to be used in place of
  any default alpha or alpha present in `c`).

The original motivation for this arrangement was that Julia "did not" support
"triangular dispatch", e.g.,
```
    convert(::Type{C{T}}, c) where {C, T}
```
On Julia v0.6 or later, parameter constraints can refer to previous parameters.
Threfore, we can use:
```
    convert(::Type{C}, c) where {T, C <: Colorant{T}}
```
However, the example above does not match `convert(RGB, c)`. Also, we should
catch all the various alpha variants (e.g. `ARGB`/`RGBA` with/without element
type).
The various arguments of `_convert` "peel off" element types (or guarantee them)
so that comparisons may be made via dispatch. Therefore, this arrangement is
still helpful.

Note that ColorTypes handles the cases where `Odest == Osrc`, or they are both
subtypes of `AbstractRGB` or `AbstractGray`. Therefore, here we only have to
deal with conversions between different spaces.
=#

function ColorTypes._convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, p, alpha=alpha(p)) where {Cdest<:TransparentColor,Odest,Osrc}
    # Convert the base color
    c = cnvt(color_type(Cdest), color(p))
    # Append the alpha
    ColorTypes._convert(Cdest, Odest, Odest, c, alpha)
end

ColorTypes._convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, c) where {Cdest<:Color,Odest,Osrc} = cnvt(Cdest, c)

# with the whitepoint `wp` as the third argument
convert(::Type{XYZ}, c, wp::XYZ) = cnvt(XYZ{eltype(wp)}, c, wp)
convert(::Type{Lab}, c, wp::XYZ) = cnvt(Lab{eltype(wp)}, c, wp)
convert(::Type{Luv}, c, wp::XYZ) = cnvt(Luv{eltype(wp)}, c, wp)

convert(::Type{XYZ{T}}, c, wp::XYZ) where {T} = cnvt(XYZ{T}, c, wp)
convert(::Type{Lab{T}}, c, wp::XYZ) where {T} = cnvt(Lab{T}, c, wp)
convert(::Type{Luv{T}}, c, wp::XYZ) where {T} = cnvt(Luv{T}, c, wp)

# FIXME: inference helpers for LCH <--> RGB conversions
convert(::Type{RGB},    c::Union{LCHab{T}, LCHuv{T}}) where {T} = cnvt(RGB{T}, cnvt(XYZ{T}, c))
convert(::Type{RGB{T}}, c::Union{LCHab{T}, LCHuv{T}}) where {T} = cnvt(RGB{T}, cnvt(XYZ{T}, c))
convert(::Type{Lab{T}}, c::RGB{T}) where {T} = cnvt(Lab{T}, cnvt(XYZ{T}, c))
convert(::Type{Luv{T}}, c::RGB{T}) where {T} = cnvt(Luv{T}, cnvt(XYZ{T}, c))

# Fallback to catch undefined operations
cnvt(::Type{C}, c::TransparentColor) where {C<:Color} = cnvt(C, color(c))
cnvt(::Type{C}, c) where {C} = convert(C, convert(RGB{eltype(C)}, c)::RGB{eltype(C)})

# Conversions from grayscale
# --------------------------
function ColorTypes._convert(::Type{Cdest}, ::Type{Odest}, ::Type{Osrc}, g) where {Cdest<:Color,Odest,Osrc<:AbstractGray}
    cnvt(Cdest, convert(RGB{eltype(Cdest)}, g))
end

# Everything to RGB
# -----------------

correct_gamut(c::CV) where {CV<:AbstractRGB} = CV(clamp01(red(c)), clamp01(green(c)), clamp01(blue(c)))
correct_gamut(c::CV) where {T<:Union{N0f8,N0f16,N0f32,N0f64},
                            CV<:Union{AbstractRGB{T},TransparentRGB{T}}} = c
correct_gamut(c::CV) where {CV<:TransparentRGB} =
    CV(clamp01(red(c)), clamp01(green(c)), clamp01(blue(c)), clamp01(alpha(c))) # for `hex`

@inline function srgb_compand(v::T) where T
    F = typeof(0.5f0v) === Float32 ? Float32 : promote_type(Float64, T)
    vf = F(v)
    vc = @fastmath max(vf, F(0.0031308))
    # `pow5_12` is an optimized function to get `v^(1/2.4)`
    vf > F(0.0031308) ? muladd(F(1.055), F(pow5_12(vc)), F(-0.055)) : F(12.92) * vf
end

function _hsx_to_rgb(im::UInt8, v, n, m)
    #=
    if     hue <  60; im = 0b000001 # ---------+
    elseif hue < 120; im = 0b000010 # --------+|
    elseif hue < 180; im = 0b000100 # -------+||
    elseif hue < 240; im = 0b001000 # ------+|||
    elseif hue < 300; im = 0b010000 # -----+||||
    else            ; im = 0b100000 # ----+|||||
    end                             #     ||||||
    (hue < 60 || hue >= 300) === ((im & 0b100001) != 0x0)
    =#
    r = ifelse((im & 0b100001) == 0x0, ifelse((im & 0b010010) == 0x0, m, n), v)
    g = ifelse((im & 0b000110) == 0x0, ifelse((im & 0b001001) == 0x0, m, n), v)
    b = ifelse((im & 0b011000) == 0x0, ifelse((im & 0b100100) == 0x0, m, n), v)
    return (r, g, b)
end
function _hsx_to_rgb(im::UInt8, v::T, n::T, m::T) where T <:Union{Float16, Float32, Float64}
    vu, nu, mu = reinterpret.(Unsigned, (v, n, m)) # prompt the compiler to use conditional moves
    r = ifelse((im & 0b100001) == 0x0, ifelse((im & 0b010010) == 0x0, mu, nu), vu)
    g = ifelse((im & 0b000110) == 0x0, ifelse((im & 0b001001) == 0x0, mu, nu), vu)
    b = ifelse((im & 0b011000) == 0x0, ifelse((im & 0b100100) == 0x0, mu, nu), vu)
    return reinterpret.(T, (r, g, b))
end

function cnvt(::Type{CV}, c::HSV) where {T, CV<:AbstractRGB{T}}
    F = promote_type(T, eltype(c))
    h, s, v = div60(F(c.h)), clamp01(F(c.s)), clamp01(F(c.v))
    hi = unsafe_trunc(Int32, h) # instead of floor
    i = h < 0 ? hi - one(hi) : hi
    f = i & one(i) == zero(i) ? 1 - (h - i) : h - i
    im = 0x1 << (mod6(UInt8, i) & 0x07)
    # use `@fastmath` just to reduce the estimated costs for inlining
    @fastmath m = v * (1 - s)
    @fastmath n = v * (1 - s * f)

    r, g, b = _hsx_to_rgb(im, v, n, m)
    T <: FixedPoint && typemax(T) >= 1 ? CV(r % T, g % T, b % T) : CV(r, g, b)
end

function cnvt(::Type{CV}, c::HSL) where {T, CV<:AbstractRGB{T}}
    F = promote_type(T, eltype(c))
    h, s, l = div60(F(c.h)), clamp01(F(c.s)), clamp01(F(c.l))
    a = @fastmath min(l, 1 - l) * s
    v = l + a
    hi = unsafe_trunc(Int32, h) # instead of floor
    i = h < 0 ? hi - one(hi) : hi
    f = i & one(i) == zero(i) ? 1 - (h - i) : h - i
    im = 0x1 << (mod6(UInt8, i) & 0x07)
    # use `@fastmath` just to reduce the estimated costs for inlining
    @fastmath m = l - a # v - 2 * a
    @fastmath n = v - 2 * a * f

    r, g, b = _hsx_to_rgb(im, v, n, m)
    T <: FixedPoint && typemax(T) >= 1 ? CV(r % T, g % T, b % T) : CV(r, g, b)
end

function cnvt(::Type{CV}, c::HSI) where {T, CV<:AbstractRGB{T}}
    F = promote_type(T, eltype(c))
    h, s, i = deg2rad(normalize_hue(F(c.h))), clamp01(F(c.s)), clamp01(F(c.i))
    is = i * s
    if h < F(2π/3)
        @fastmath cosr = cos(h) / cos(F(π/3)-h)
        r0, g0, b0 = muladd(is, cosr, i), muladd(is, 1-cosr, i), i - is
    elseif h < F(4π/3)
        @fastmath cosr = cos(h-F(2π/3)) / cos(F(π)-h)
        r0, g0, b0 = i - is, muladd(is, cosr, i), muladd(is, 1-cosr, i)
    else
        @fastmath cosr = cos(h-F(4π/3)) / cos(F(5π/3)-h)
        r0, g0, b0 = muladd(is, 1-cosr, i), i - is, muladd(is, cosr, i)
    end
    r, g, b = min(r0, oneunit(F)), min(g0, oneunit(F)), min(b0, oneunit(F))
    T <: FixedPoint && typemax(T) >= 1 ? CV(r % T, g % T, b % T) : CV(r, g, b)
end

# the following matrix is based on the sRGB color primaries in `xy` and D65 whitepoint in `XYZ`
const M_XYZ2RGB = Mat3x3([ 3.2404541621141054   -1.5371385127977166   -0.4985314095560162
                          -0.9692660305051868    1.8760108454466942    0.04155601753034984
                           0.05564343095911469  -0.20402591351675387   1.0572251882231791 ])
function xyz_to_linear_rgb(c::XYZ{T}) where T
    @mul3x3 RGB M_XYZ2RGB c.x c.y c.z
end
function cnvt(::Type{CV}, c::XYZ) where CV<:AbstractRGB
    rgb = xyz_to_linear_rgb(c)
    CV(clamp01(srgb_compand(rgb.r)),
       clamp01(srgb_compand(rgb.g)),
       clamp01(srgb_compand(rgb.b)))
end

const M_YIQ2RGB = Mat3x3([1.0   0.9563   0.621
                          1.0  -0.2721  -0.6474
                          1.0  -1.107    1.7046 ])
function cnvt(::Type{CV}, c::YIQ) where CV<:AbstractRGB
    cc = correct_gamut(c)
    rgb = @mul3x3 RGB M_YIQ2RGB cc.y cc.i cc.q
    CV(clamp01(rgb.r), clamp01(rgb.g), clamp01(rgb.b))
end

# FIXME
function cnvt(::Type{CV}, c::YCbCr) where CV<:AbstractRGB
    cc = correct_gamut(c)
    ny = cc.y - 16
    ncb = cc.cb - 128
    ncr = cc.cr - 128
    CV(clamp01(0.004567ny - 1.39135e-7ncb + 0.0062586ncr),
       clamp01(0.004567ny - 0.00153646ncb - 0.0031884ncr),
       clamp01(0.004567ny + 0.00791058ncb - 2.79201e-7ncr))
end

# To avoid stack overflow, the source types which do not support direct or
# indirect conversion to RGB should be rejected.
cnvt(::Type{CV}, c::Union{LMS, xyY}              ) where {CV<:AbstractRGB} = cnvt(CV, cnvt(XYZ{eltype(c)}, c))
cnvt(::Type{CV}, c::Union{Lab, Luv, LCHab, LCHuv}) where {CV<:AbstractRGB} = cnvt(CV, cnvt(XYZ{eltype(c)}, c))
cnvt(::Type{CV}, c::Union{DIN99d, DIN99o, DIN99} ) where {CV<:AbstractRGB} = cnvt(CV, cnvt(XYZ{eltype(c)}, c))
@noinline function cnvt(::Type{CV}, @nospecialize(c::Color)) where {CV<:AbstractRGB}
    error("No conversion of ", c, " to ", CV, " has been defined")
end

# AbstractGray --> AbstractRGB conversions are implemented in ColorTypes.jl

# Everything to HSV
# -----------------

function cnvt(::Type{HSV{T}}, c::AbstractRGB) where T
    F = promote_type(T, eltype(c))
    r, g, b = F.((red(c), green(c), blue(c)))
    c_min = @fastmath min(min(r, g), b)
    c_max = @fastmath max(max(r, g), b)
    s0 = c_max - c_min
    s0 == zero(F) && return HSV{T}(zero(T), zero(T), T(c_max))
    s = @fastmath s0 / c_max

    # In general, it is dangerous to compare floating point numbers with `===`.
    diff = ifelse(c_max === r,  g - b,         ifelse(c_max === g,  b - r,  r - g))
    ofs  = ifelse(c_max === r, (g < b)*F(360), ifelse(c_max === g, F(120), F(240)))
    h0 = @fastmath diff * F(60) / s0

    HSV{T}(h0 + ofs, s, c_max)
end

cnvt(::Type{HSV{T}}, c::Color) where {T} = cnvt(HSV{T}, convert(RGB{T}, c)::RGB{T})

# Everything to HSL
# -----------------

function cnvt(::Type{HSL{T}}, c::AbstractRGB) where T
    F = promote_type(T, eltype(c))
    r, g, b = F(red(c)), F(green(c)), F(blue(c))
    c_min = @fastmath min(min(r, g), b)
    c_max = @fastmath max(max(r, g), b)
    l0 = c_max + c_min
    s0 = c_max - c_min
    l = l0 * F(0.5)
    s0 == zero(F) && return HSL{T}(zero(T), zero(T), T(l))
    s = @fastmath s0 / min(l0, F(2) - l0)

    # In general, it is dangerous to compare floating point numbers with `===`.
    diff = ifelse(c_max === r,  g - b,         ifelse(c_max === g,  b - r,  r - g))
    ofs  = ifelse(c_max === r, (g < b)*F(360), ifelse(c_max === g, F(120), F(240)))
    h0 = @fastmath diff * F(60) / s0

    HSL{T}(h0 + ofs, s, l)
end

cnvt(::Type{HSL{T}}, c::Color) where {T} = cnvt(HSL{T}, convert(RGB{T}, c)::RGB{T})

# Everything to HSI
# -----------------

# Since acosd() is slow, the following is "inline-worthy".
@inline function cnvt(::Type{HSI{T}}, c::AbstractRGB) where T
    rgb = correct_gamut(c)
    F = promote_type(T, eltype(c))
    r, g, b = F(red(rgb)), F(green(rgb)), F(blue(rgb))
    dnorm = @fastmath sqrt(((r-g)^2 + (r-b)^2 + (g-b)^2) * F(0.5))
    isum = r + g + b
    i = isum / F(3)
    dnorm == zero(F) && return HSI{T}(T(90), zero(T), T(i))
    val = muladd(g + b, F(-0.5), r) / dnorm
    h = @fastmath acosd(clamp(val, -oneunit(F), oneunit(F)))
    m = @fastmath min(min(r, g), b)
    s = oneunit(F) - m/i
    HSI{T}(b > g ? F(360) - h : h, s, i)
end

cnvt(::Type{HSI{T}}, c::Color) where {T} = cnvt(HSI{T}, convert(RGB{T}, c)::RGB{T})

# Everything to XYZ
# -----------------

@inline function invert_srgb_compand(v::T) where T
    F = typeof(0.5f0v) === Float32 ? Float32 : promote_type(Float64, T)
    vf = F(v)
    # `pow12_5` is an optimized function to get `x^2.4`
    vf > F(0.04045) ? pow12_5(muladd(F(1000/1055), vf, F(55/1055))) : F(100/1292) * vf
end

# lookup table for `N0f8` (the extra two elements are for `Float32` splines)
const invert_srgb_compand_n0f8 = Float32[invert_srgb_compand(v/255.0) for v = 0:257]

function invert_srgb_compand(v::N0f8)
    @inbounds invert_srgb_compand_n0f8[reinterpret(v) + 1]
end

function invert_srgb_compand(v::Float32)
    i = unsafe_trunc(Int32, v * 255)
    (i < 13 || i > 255) && return Float32(invert_srgb_compand(Float64(v)))
    @inbounds y = view(invert_srgb_compand_n0f8, i:i+3)
    dv = v * 255.0f0 - i
    dv == 0.0f0 && @inbounds return y[2]
    if v < 0.38857287f0
        return @fastmath(y[2]+0.5f0*dv*((-2/3f0*y[1]- y[2])+(2y[3]-1/3f0*y[4])+
                                    dv*((       y[1]-2y[2])+  y[3]-
                                    dv*(( 1/3f0*y[1]- y[2])+( y[3]-1/3f0*y[4]) ))))
    else
        return @fastmath(y[2]+0.5f0*dv*((4y[3]-3y[2])-y[4]+dv*((y[4]-y[3])+(y[2]-y[3]))))
    end
end

# the following matrix is based on the sRGB color primaries in `xy` and D65 whitepoint in `XYZ`
const M_RGB2XYZ = Mat3x3([0.4124564390896921    0.357576077643909  0.18043748326639894
                          0.21267285140562248   0.715152155287818  0.07217499330655958
                          0.019333895582329317  0.119192025881303  0.9503040785363677 ])

function linear_rgb_to_xyz(c::AbstractRGB{T}) where T
    @mul3x3 XYZ M_RGB2XYZ red(c) green(c) blue(c)
end
function cnvt(::Type{XYZ{T}}, c::AbstractRGB) where T
    r = invert_srgb_compand(red(c))
    g = invert_srgb_compand(green(c))
    b = invert_srgb_compand(blue(c))
    @mul3x3 XYZ{T} M_RGB2XYZ r g b
end

function cnvt(::Type{XYZ{T}}, c::xyY) where T
    X = c.Y*c.x/c.y
    Z = c.Y*(1-c.x-c.y)/c.y
    XYZ{T}(X, c.Y, Z)
end

const xyz_epsilon   = 216 / 24389 # (6/29)^3
const xyz_kappa     = 24389 / 27  # (29/6)^3*8
const xyz_kappa_inv = 27 / 24389

function lab2xyz(::Type{XYZ{T}}, c) where T
    F = promote_type(T, eltype(c))

    fy1 = c.l * F(0x1p-7)
    fy2 = muladd(c.l, F(3 / 3712), F(16 / 116))
    fy = fy1 + fy2 # (c.l + 16) / 116
    fx = fy1 + muladd(c.a, F( 0x1p-9), muladd(c.a, F(3 / 64000), fy2)) # fy + c.a / 500
    fz = fy1 + muladd(c.b, F(-0x1p-8), muladd(c.b, F(-7 / 6400), fy2)) # fy - c.b / 200

    fx3 = fx^3
    fy3 = fy^3
    fz3 = fz^3

    epsilon = oftype(fx3, xyz_epsilon)
    kappa_inv = oftype(fx3, xyz_kappa_inv)

    x = fx3 > epsilon ? fx3 : muladd(116, fx, -16) * kappa_inv
    y = fy3 > epsilon ? fy3 : c.l * kappa_inv
    z = fz3 > epsilon ? fz3 : muladd(116, fz, -16) * kappa_inv
    return XYZ{T}(x, y, z)
end

function cnvt(::Type{XYZ{T}}, c::Lab) where T
    xyz = lab2xyz(XYZ{T}, c)
    XYZ{T}(xyz.x * T(WP_DEFAULT.x), xyz.y * T(WP_DEFAULT.y), xyz.z * T(WP_DEFAULT.z))
end
function cnvt(::Type{XYZ{T}}, c::Lab, wp::XYZ) where T
    xyz = lab2xyz(XYZ{T}, c)
    XYZ{T}(xyz.x * T(wp.x), xyz.y * T(wp.y), xyz.z * T(wp.z))
end

function xyz_to_uv(c::XYZ)
    F = typeof(4c.x / oneunit(c.x))
    d = c.x + 15c.y + 3c.z
    d == zero(d) && return (zero(F), zero(F))
    return (4c.x, 9c.y) ./ d # (u, v)
end

@inline function luv2xyz(::Type{XYZ{T}}, c, y_wp::F, uv_wp::Tuple{F, F}) where {T, F}
    c.l == zero(c.l) && return XYZ{T}(zero(T), zero(T), zero(T))
    ls = c.l * F(xyz_kappa_inv)
    clc = @fastmath max(F(c.l), F(8))
    lc = muladd(clc, F(0x1p-7), muladd(clc, F(3 / 3712), F(16 / 116)))^3 # (c.l + 16) / 116)^3
    y = y_wp * ifelse(c.l > 8, lc, ls)
    u, v = (c.u, c.v) ./ (13c.l) .+ uv_wp
    ys = F(0.25) * y / v
    x = ys * 9u
    z = ys * (12 - 3u - 20v)
    XYZ{T}(x, y, z)
end
function cnvt(::Type{XYZ{T}}, c::Luv) where T
    F = promote_type(T, eltype(c))
    luv2xyz(XYZ{T}, c, F(WP_DEFAULT.y), F.(xyz_to_uv(WP_DEFAULT)))
end
function cnvt(::Type{XYZ{T}}, c::Luv, wp::XYZ) where T
    F = promote_type(T, eltype(c))
    luv2xyz(XYZ{T}, c, F(wp.y), F.(xyz_to_uv(wp)))
end

function cnvt(::Type{XYZ{T}}, c::DIN99d) where T

    # Go back to C-h space
    C = chroma(c)
    h = atan(c.b, c.a) - 50π/180

    # Intermediate terms
    G = (exp(C/22.5)-1)/0.06
    f, ee = G .* sincos(h)

    l = (exp(c.l/325.221)-1)/0.0036
    # a = ee*cosd(50) - f/1.14*sind(50)
    a = ee*0.6427876096865394 - f/1.14*0.766044443118978
    # b = ee*sind(50) - f/1.14*cosd(50)
    b = ee*0.766044443118978 + f/1.14*0.6427876096865394

    adj = convert(XYZ, Lab(l, a, b))

    XYZ{T}((adj.x + 0.12*adj.z)/1.12, adj.y, adj.z)
end

function cnvt(::Type{XYZ{T}}, c::LMS) where T
    @mul3x3 XYZ{T} CAT02_INV c.l c.m c.s
end

cnvt(::Type{XYZ{T}}, c::Union{LCHab, DIN99, DIN99o}) where {T} = cnvt(XYZ{T}, cnvt(Lab{T}, c))
cnvt(::Type{XYZ{T}}, c::LCHuv) where {T} = cnvt(XYZ{T}, cnvt(Luv{T}, c))
cnvt(::Type{XYZ{T}}, c::Color) where {T} = cnvt(XYZ{T}, convert(RGB{T}, c)::RGB{T})

# Everything to xyY
# -----------------

function cnvt(::Type{xyY{T}}, c::XYZ) where T

    x = c.x/(c.x + c.y + c.z)
    y = c.y/(c.x + c.y + c.z)

    xyY{T}(x, y, oftype(x, c.y))
end

cnvt(::Type{xyY{T}}, c::Color) where {T} = cnvt(xyY{T}, convert(XYZ{T}, c)::XYZ{T})

# Everything to Lab
# -----------------

@inline function fxyz2lab(v)
    ka = oftype(v, 841 / 108) # (29/6)^2 / 3 = xyz_kappa / 116
    kb = oftype(v, 16 / 116) # 4/29
    vc = @fastmath max(v, oftype(v, xyz_epsilon))
    @fastmath min(cbrt01(vc), muladd(ka, v, kb))
end
@inline function xyz2lab(::Type{Lab{T}}, c::XYZ) where T
    f = XYZ(fxyz2lab(c.x), fxyz2lab(c.y), fxyz2lab(c.z)) # mapc(fxyz2lab, c)
    Lab{T}(116f.y - 16, 500(f.x - f.y), 200(f.y - f.z))
end
function cnvt(::Type{Lab{T}}, c::XYZ) where T
    wp = WP_DEFAULT
    F = promote_type(T, eltype(c))
    xyz2lab(Lab{T}, XYZ(c.x * F(1 / wp.x), c.y * F(1 / wp.y), c.z * F(1 / wp.z)))
end
function cnvt(::Type{Lab{T}}, c::XYZ, wp::XYZ) where T
    F = promote_type(T, eltype(c))
    xyz2lab(Lab{T}, XYZ(c.x / F(wp.x), c.y / F(wp.y), c.z / F(wp.z)))
end

function cnvt(::Type{Lab{T}}, c::LCHab) where T
    Lab{T}(c.l, polar_to_cartesian(c.c, c.h)...)
end

function cnvt(::Type{Lab{T}}, c::DIN99) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate Chroma (C99) in the DIN99 space
    cc = chroma(c)

    # Temporary variable for chroma
    g = (exp(0.045*cc*kch*ke)-1)/0.045

    # Temporary redness
    ee = cc > 0 ? g * c.a / cc : zero(g)

    # Temporary yellowness
    f = cc > 0 ? g * c.b / cc : zero(g)

    # CIELAB a*b*
    # ciea = ee*cosd(16) - (f/0.7)*sind(16)
    ciea = ee*0.9612616959383189 - (f/0.7)*0.27563735581699916
    # cieb = ee*sind(16) + (f/0.7)*cosd(16)
    cieb = ee*0.27563735581699916 + (f/0.7)*0.9612616959383189

    # CIELAB L*
    ciel = (exp(c.l*ke/105.51)-1)/0.0158

    Lab{T}(ciel, ciea, cieb)
end

function cnvt(::Type{Lab{T}}, c::DIN99o) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate Chroma (C99) in the DIN99o space
    co = chroma(c)

    # hue angle h99o
    h = atan(c.b, c.a) - 26π/180

    # revert logarithmic chroma compression
    g = (exp(co*kch*ke/23.0)-1)/0.075

    # Temporary yellowness and redness
    fo, eo = g .* sincos(h)

    # CIELAB a*b* (revert b* axis compression)
    # ciea = eo*cosd(26) - (fo/0.83)*sind(26)
    ciea = eo*0.898794046299167 - (fo/0.83)*0.4383711467890774
    # cieb = eo*sind(26) + (fo/0.83)*cosd(26)
    cieb = eo*0.4383711467890774 + (fo/0.83)*0.898794046299167

    # CIELAB L* (revert logarithmic lightness compression)
    ciel = (exp(c.l*ke/303.67)-1)/0.0039

    Lab{T}(ciel, ciea, cieb)
end

cnvt(::Type{Lab{T}}, c::Color) where {T} = cnvt(Lab{T}, convert(XYZ{T}, c)::XYZ{T})

# Everything to Luv
# -----------------
@inline function xyz2luv(::Type{Luv{T}}, c::XYZ, y::F, uv_wp::Tuple{F, F}) where {T, F}
    yc = @fastmath max(y, F(xyz_epsilon))
    l = @fastmath min(muladd(116, cbrt01(yc), -16), F(xyz_kappa) * y)
    uv = (13 * l) .* (xyz_to_uv(c) .- uv_wp)
    Luv{T}(l, uv...)
end

function cnvt(::Type{Luv{T}}, c::XYZ) where T
    F = promote_type(T, eltype(c))
    y = c.y * F(1 / WP_DEFAULT.y)
    xyz2luv(Luv{T}, c, F(y), F.(xyz_to_uv(WP_DEFAULT)))
end
function cnvt(::Type{Luv{T}}, c::XYZ, wp::XYZ) where T
    F = promote_type(T, eltype(c))
    y = c.y / F(wp.y)
    xyz2luv(Luv{T}, c, F(y), F.(xyz_to_uv(wp)))
end

function cnvt(::Type{Luv{T}}, c::LCHuv) where T
    Luv{T}(c.l, polar_to_cartesian(c.c, c.h)...)
end

cnvt(::Type{Luv{T}}, c::Color) where {T} = cnvt(Luv{T}, convert(XYZ{T}, c)::XYZ{T})

# Everything to LCHuv
# -------------------

function cnvt(::Type{LCHuv{T}}, c::Luv) where T
    LCHuv{T}(c.l, chroma(c), hue(c))
end

cnvt(::Type{LCHuv{T}}, c::Color) where {T} = cnvt(LCHuv{T}, convert(Luv{T}, c)::Luv{T})

# Everything to LCHab
# -------------------

function cnvt(::Type{LCHab{T}}, c::Lab) where T
    LCHab{T}(c.l, chroma(c), hue(c))
end

cnvt(::Type{LCHab{T}}, c::Color) where {T} = cnvt(LCHab{T}, convert(Lab{T}, c)::Lab{T})

# Everything to DIN99
# -------------------

function cnvt(::Type{DIN99{T}}, c::Lab) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate DIN99 L
    l99 = (1/ke)*105.51*log(1+0.0158*c.l)

    # Temporary value for redness and yellowness
    # ee = c.a*cosd(16) + c.b*sind(16)
    ee = c.a*0.9612616959383189 + c.b*0.27563735581699916
    # f = -0.7*c.a*sind(16) + 0.7*c.b*cosd(16)
    f = -0.7*c.a*0.27563735581699916 + 0.7*c.b*0.9612616959383189

    # Temporary value for chroma
    g = sqrt(ee^2 + f^2)

    # DIN99 chroma
    cc = log(1+0.045*g)/(0.045*kch*ke)

    # DIN99 chromaticities
    a99 = g > 0 ? convert(T, cc / g * ee) : zero(T)
    b99 = g > 0 ? convert(T, cc / g * f ) : zero(T)

    DIN99{T}(l99, a99, b99)

end

cnvt(::Type{DIN99{T}}, c::Color) where {T} = cnvt(DIN99{T}, convert(Lab{T}, c)::Lab{T})

# Everything to DIN99d
# --------------------

function cnvt(::Type{DIN99d{T}}, c::XYZ{T}) where T

    # Apply tristimulus-space correction term
    adj_c = XYZ(1.12*c.x - 0.12*c.z, c.y, c.z)

    # Apply L*a*b*-space correction
    lab = convert(Lab, adj_c)
    adj_L = 325.221*log(1+0.0036*lab.l)

    # Calculate intermediate parameters
    # ee = lab.a*cosd(50) + lab.b*sind(50)
    ee = lab.a*0.6427876096865394 + lab.b*0.766044443118978
    # f = 1.14*(lab.b*cosd(50) - lab.a*sind(50))
    f = 1.14*(lab.b*0.6427876096865394 - lab.a*0.766044443118978)
    G = sqrt(ee^2+f^2)

    # Calculate hue/chroma
    C = 22.5*log(1+0.06*G)
    h = atan(f, ee) + 50π/180

    b99, a99 = C .* sincos(h)

    DIN99d{T}(adj_L, a99, b99)

end

cnvt(::Type{DIN99d{T}}, c::Color) where {T} = cnvt(DIN99d{T}, convert(XYZ{T}, c)::XYZ{T})

# Everything to DIN99o
# -------------------

function cnvt(::Type{DIN99o{T}}, c::Lab) where T

    # We assume the adjustment parameters are always 1; the standard recommends
    # that they not be changed from these values.
    kch = 1
    ke = 1

    # Calculate DIN99o L (logarithmic compression)
    l99 = 303.67/ke*log(1+0.0039*c.l)

    # Temporary value for redness and yellowness
    # including rotation by 26°
    # eo = c.a*cosd(26) + c.b*sind(26)
    eo = c.a*0.898794046299167 + c.b*0.4383711467890774
    # compression along the yellowness (blue-yellow) axis
    # fo = 0.83 * (c.b*cosd(26) - c.a*sind(26))
    fo = 0.83 * (c.b*0.898794046299167 - c.a*0.4383711467890774)

    # Temporary value for chroma
    go = sqrt(eo^2 + fo^2)
    h = atan(fo, eo) + 26π/180

    # DIN99o chroma (logarithmic compression)
    cc = 23.0*log(1+0.075*go)/(kch*ke)

    # DIN99o chromaticities
    b99, a99 = cc .* sincos(h)

    DIN99o{T}(l99, a99, b99)

end

cnvt(::Type{DIN99o{T}}, c::Color) where {T} = cnvt(DIN99o{T}, convert(Lab{T}, c)::Lab{T})

# Everything to LMS
# -----------------

# Chromatic adaptation from CIECAM97s
const CAT97s = Mat3x3([ 0.8562  0.3372 -0.1934
                       -0.8360  1.8327  0.0033
                        0.0357 -0.0469  1.0112 ])

const CAT97s_INV = Mat3x3(inv(Float64.(CAT97s)))

# Chromatic adaptation from CIECAM02
const CAT02 = Mat3x3([ 0.7328 0.4296 -0.1624
                      -0.7036 1.6975  0.0061
                       0.0030 0.0136  0.9834 ])

const CAT02_INV = Mat3x3(inv(Float64.(CAT02)))

function cnvt(::Type{LMS{T}}, c::XYZ) where T
    @mul3x3 LMS{T} CAT02 c.x c.y c.z
end

cnvt(::Type{LMS{T}}, c::Color) where {T} = cnvt(LMS{T}, convert(XYZ{T}, c)::XYZ{T})

# Everything to YIQ
# -----------------

correct_gamut(c::YIQ{T}) where {T} = YIQ{T}(clamp(c.y, zero(T), one(T)),
                                     clamp(c.i, convert(T,-0.595716), convert(T,0.595716)),
                                     clamp(c.q, convert(T,-0.522591), convert(T,0.522591)))

const M_RGB2YIQ = Mat3x3([0.299      0.587      0.114
                          0.595716  -0.274453  -0.321263
                          0.211456  -0.522591   0.311135 ])
function cnvt(::Type{YIQ{T}}, c::AbstractRGB) where T
    rgb = correct_gamut(c)
    @mul3x3 YIQ{T} M_RGB2YIQ red(rgb) green(rgb) blue(rgb)
end

cnvt(::Type{YIQ{T}}, c::Color) where {T} = cnvt(YIQ{T}, convert(RGB{T}, c)::RGB{T})

# Everything to YCbCr
# -------------------
# FIXME
correct_gamut(c::YCbCr{T}) where {T} = YCbCr{T}(clamp(c.y, convert(T,16), convert(T,235)),
                                         clamp(c.cb, convert(T,16), convert(T,240)),
                                         clamp(c.cr, convert(T,16), convert(T,240)))
# FIXME
function cnvt(::Type{YCbCr{T}}, c::AbstractRGB) where T
    rgb = correct_gamut(c)
    YCbCr{T}(16+65.481*red(rgb)+128.553*green(rgb)+24.966*blue(rgb),
             128-37.797*red(rgb)-74.203*green(rgb)+112*blue(rgb),
             128+112*red(rgb)-93.786*green(rgb)-18.214*blue(rgb))
end

cnvt(::Type{YCbCr{T}}, c::Color) where {T} = cnvt(YCbCr{T}, convert(RGB{T}, c)::RGB{T})

# To Gray
# -------
# AbstractGray --> AbstractRGB conversions are implemented in ColorTypes.jl, but
# AbstractRGB --> AbstractGray conversions should be implemented here.

# Rec 601 luma conversion

function cnvt(::Type{G}, x::AbstractRGB{T}) where {G<:AbstractGray,T<:Normed}
    TU, Tf = FixedPointNumbers.rawtype(T), floattype(T)
    if sizeof(TU) < sizeof(UInt)
        val = Tf(0.001)*(299*reinterpret(red(x)) + 587*reinterpret(green(x)) + 114*reinterpret(blue(x)))
    else
        val = Tf(0.299)*reinterpret(red(x)) + Tf(0.587)*reinterpret(green(x)) + Tf(0.114)*reinterpret(blue(x))
    end
    return G(reinterpret(T, round(TU, val)))
end
cnvt(::Type{G}, x::AbstractRGB) where {G<:AbstractGray} =
    G(0.299f0*red(x) + 0.587f0*green(x) + 0.114f0*blue(x))

cnvt(::Type{G}, x::Color) where {G<:AbstractGray} = convert(G, convert(RGB, x))


text:
function solve(laplace::AnalyticalLaplace)
    (; D, T₂, dim, eiglim, eigstep) = laplace

    # TODO: Add support for dim = 1 (PlateSetup)
    L = length(D)

    # Compute radial and angular eigenvalues (λ and ν), with α^2 = λ and n^2 = ν for
    # cylinders and n(n+1) = ν for spheres
    dα = √eigstep * 1e3
    α_min = dα / 100
    α_max = √eiglim * 1e3
    α, n = find_α(laplace, α_min, α_max, dα)

    N = length(α)
    λ = α .^ 2
    Λ = diagm(λ)

    β = zeros(N)
    J = zeros(N)
    for m = 1:N
        β[m] = compute_β(laplace, α[m], n[m])
        if n[m] == 0
            J[m] = compute_int_J(laplace, α[m])[1]
        end
    end

    if dim == 2
        coeff = 2
    else
        coeff = √6
    end
    U = coeff * J .* β

    B = zeros(N, N)
    for m1 = 1:N
        for m2 = m1:N
            if abs(n[m1] - n[m2]) == 1
                K = compute_int_K(laplace, α[m1], n[m1], α[m2], n[m2])
                if dim == 2
                    ε = √(1 + (n[m1] == 0) + (n[m2] == 0))
                else
                    ε = (n[m1] + n[m2] + 1) / √((2 * n[m1] + 1) * (2 * n[m2] + 1))
                end
                B[m1, m2] = ε * β[m1] * β[m2] * K
                B[m2, m1] = B[m1, m2]
            end
        end
    end

    Bri = zeros(N, N, L)
    for m1 = 1:N
        for m2 = m1:N
            if n[m1] == n[m2]
                int_I = compute_int_I(laplace, α[m1], α[m2], n[m1])
                Bri[m1, m2, :] = 2 * β[m1] * β[m2] * int_I
                Bri[m2, m1, :] = Bri[m1, m2, :]
            end
        end
    end

    Br = zeros(N, N)
    for i = 1:L
        Br .+= Bri[:, :, i] / T₂[i]
    end

    (; Λ, β, U, B, Bri, Br)
end

ext: jl


text:
using Test
using TestSetExtensions

include("../src/Groebner.jl")

using .Groebner
using .Groebner.AbstractAlgebra

# Taken from JuMP/test/solvers.jl
function try_import(name::Symbol)
    try
        @eval import $name
        return true
    catch e
        return false
    end
end

⊂(xs, ys) = all(in(ys), xs)
# set equality
≂(xs, ys) = ⊂(xs, ys) && ⊂(ys, xs)

@testset "All tests" begin

    @includetests ["commons", "io", "io_consistency"]

    @includetests ["coeffs", "sorting",
                    "mod_reconstruction", "crt_reconstruction"]

    @includetests ["core_f4_reduce", "core_f4_stress",
                    "core_f4", "rational_f4"]

    @includetests ["core_isgroebner", "core_isgroebner_stress"]

    @includetests ["core_normalform", "core_normalform_stress"]

    @includetests ["hard_problems_f4", "large_problems_f4"]

    if try_import(:DynamicPolynomials)
        @includetests ["io_dynamic", "core_dynamic"]
    end
end


text:
"""
    track(A::Matrix{Float64}, p, h = 10) -> p2, err

Track blurry lightsource by applying a window with half-width `h` at
an approximate location `p = (i,j)` and find the
average weighted location of points with *high* light intensity.

Gives an standard error estimate.
"""
function track(img::Matrix, p, h = 10)
    i, j = round.(Int, p)
    CR = CartesianRange(CartesianIndex(i - h, j - h), CartesianIndex(i + h, j + h))
    μ = mean(img[ci] for ci in CR)
    C = sum(max(img[ci] - μ, 0) for ci in CR)

    xhat = sum(max(img[ci] - μ, 0)*ci[1] for ci in CR)/C
    yhat = sum(max(img[ci] - μ, 0)*ci[2] for ci in CR)/C

    xerr =  sum(max(img[ci] - μ, 0)*(ci[1] - xhat)^2 for ci in CR)/C
    yerr =  sum(max(img[ci] - μ, 0)*(ci[2] - yhat)^2 for ci in CR)/C

    err = sqrt(xerr + yerr)
    (xhat, yhat), err
end


text:
module Utils

using GraphUtils

export
    interpolate,
    linear_interp,
    prune_sorted_array

interpolate(a,b,t) = (1 - t)*a + t*b

function linear_interp(v,t)
    idx = find_index_in_sorted_array(v,t)-1
    idx = max(1, min(length(v)-1, idx))
    if isapprox(v[idx+1] - v[idx],0.0)
        return idx,0.0
    end
    Δt = (t - v[idx]) / (v[idx+1] - v[idx])
    return idx, Δt
end

"""
    `prune_array(arr,v,t1,t2)`

    prunes `arr` on either side based on where t1 and t2 fall in v.
"""
function prune_sorted_array(v,idx1,idx2,Δt1,Δt2)
    @assert [idx1,Δt1] <= [idx2,Δt2] "[idx1,Δt1] > [idx2,Δt2]"
    vp = [interpolate(v[idx1],v[idx1+1],Δt1)]
    for idx in min(idx1+1,length(v)):min(idx2,length(v))
        push!(vp,v[idx])
    end
    push!(vp,interpolate(v[idx2],v[idx2+1],Δt2))
    return vp
end

end


text:
#=
nonlinearize:
- Julia version: 1.3.1
- Author: Berk
- Date: 2020-11-24
=#

vars = flat(expr_vars)
aux_vars = [Symbol("aux$(i)") for i=1:length(vars)]
varmap = get_varmap(expr_vars, vars)
lhs = :(($(aux_vars)...))
vars_copy = deepcopy(expr.args[1].args)
expr_copy = deepcopy(expr.args[2])

for i = 1:length(expr_vars)
    tuple_idxs = findall(x -> x[1] == i, varmap)
    if length(tuple_idxs) == 1 # Scalar inputs to the function
        expr_copy = substitute(expr_copy, vars_copy[i]=> aux_vars[tuple_idxs[1]])
    else
        expr_copy = substitute(expr_copy, vars_copy[i] => aux_vars[tuple_idxs])
    end
end
fn = :($(lhs) -> $(expr_copy))


text:
using Test
using Liblet

test_item_wrong_lhs() = Item(1, ["a"])

test_item_neg_pos() = Item("A", ["B"], -1)

test_item_wrong_pos() = Item("A", ["B"], 3)

test_item_inset() = begin
    P = Item("a",["b","c"])
    Q = Item("a",["b", "c"])
    s = Set()
    push!(s, P)
    push!(s, Q)
    length(s) == 1
end

test_item_totalorder() = Item("A",["B","C"],2) > Item("a",["b","c"],1)

test_item_advance() = Item("A", ("x", "B"), 2) == advance(Item("A", ("x", "B")),"x")

test_item_notadvance() = advance(Item("A",["B","C"]),"x") === nothing

test_item_symbol_after_dot() = "B" == (Item("A",["x","B"],2) |> afterdotsymbol)

test_item_nosymbol_after_dot() = (Item("A", ("B", "C"), 3) |> afterdotsymbol) === nothing

function runitemtests()
    @testset "item_tests" begin
        # constructors
        @test_throws ArgumentError test_item_wrong_lhs()
        @test_throws ArgumentError test_item_neg_pos()
        @test_throws ArgumentError test_item_wrong_pos()
        @test test_item_inset()
        @test test_item_totalorder()

        # advance
        @test test_item_advance()
        @test test_item_notadvance()

        # afterdotsymbol
        @test test_item_symbol_after_dot()
        @test test_item_nosymbol_after_dot()
    end
end


text:
# ---
# title: 504. Base 7
# id: problem504
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/base-7/description/>
# hidden: true
# ---
#
# Given an integer, return its base 7 string representation.
#
# **Example 1:**
#
#
#
#     Input: 100
#     Output: "202"
#
#
# **Example 2:**
#
#
#
#     Input: -7
#     Output: "-10"
#
#
# **Note:** The input will be in range of [-1e7, 1e7].
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
include("bench.jl")
include("visualize.jl")
include("place.jl")
import Profile
import ProfileView

function test_profile()
    Profile.clear()
    Profile.init(n = 10^7, delay = 0.01)

    open("/tmp/prof.txt", "w") do s
        # I want to ignore all low-cost operations (e.g. <10)
        Profile.print(IOContext(s, :displaysize => (24, 200)), mincount=10)
    end

    # https://discourse.julialang.org/t/profiling-code-reveals-that-most-time-is-used-calling-3-functions-in-task-jl/30751
    #
    # profiling is showing mainly task.jl staff. This does not work:
    #
    # ENV["JULIA_NUM_THREADS"] = 1

    # other visualization methods
    #
    # Profile.print()
    # import ProfileView
    ProfileView.view()
end

function test()
    str = open("/tmp/rackematic/out/gh60.json") do io
        read(io, String)
    end;
    jobj = JSON.parse(str)

    xs, ys, ws, hs, Es, mask, diearea, params = decode_place_spec(jobj);

    # place(xs, ys, ws, hs, Es, mask, diearea, vis=true)

    solxs, solys = place(xs, ys, ws, hs, Es, mask, diearea, vis=true, nsteps=50)
    simulated_annealing_legalization(solxs, solys, ws, hs, mask, diearea)

    @time solxs, solys = place(xs, ys, ws, hs, Es, mask, diearea, vis=false, nsteps=100)
    Profile.@profile solxs, solys = place(xs, ys, ws, hs, Es, mask, diearea, nsteps=20)
    # save to json file
    res_payload = Dict("xs"=>solxs, "ys"=>solys) |> JSON.json
    open("/tmp/rackematic/out/gh60-sol.json", "w") do io
        write(io, res_payload)
    end
end

max_stars_repo_name: lihebi/bhdl

max_issues_repo_name: lihebi/bhdl

max_forks_repo_name: lihebi/bhdl


text:
using BERT

config = BertConfig(128, 30022, 256, 512, 4, 2, 8, 2, 3, Array{Float32}, 0.1, 0.1, "relu")

model = BertPreTraining(config)

x = [213 234 7789; 712 9182 8912; 7812 12 432; 12389 1823 8483] # 4x3
segment_ids = [1 1 1;1 2 1;1 2 1;1 1 1]
mlm_labels = [-1 234 -1; -1 -1 8912; -1 -1 -1; 12389 -1 -1]
nsp_labels = [1, 2, 1]

loss = model(x, segment_ids, mlm_labels, nsp_labels)
println(loss)

max_issues_repo_name: OsmanMutlu/BERT-in-Knet

max_forks_repo_name: OsmanMutlu/BERT-in-Knet


text:
##########################################################################
#       Some useful models and observables to get started with           #
##########################################################################

using LinearAlgebra

"""
    getPauliMatrices()

Provide the Pauli matrices as dense matrices.
"""
function getPauliMatrices()
    Id = [1.0 0;0.0 1.0]
    X = [0.0 1.0;1.0 0.0]
    Y = [0.0 -1.0im; 1.0im 0.0]
    Z = [1.0 0.0;0.0 -1.0]
    return Id, X, Y, Z
end

"""
    getHeisenbergMPO(N::Int64,J::Float64,lambda::Float64)

MPO implementation of the Heisenberg Hamiltonian
H = J * sum_{i) (X^{i} X^{i+1} +Y^{i} Y^{i+1} + Z^{i} Z^{i+1}) + λ * sum_i X^i
where J, λ have to be real numbers .
"""
function getHeisenbergMPO(N::Int, J::Real, lambda::Real)::MPO{ComplexF64}
    # Provide the Pauli matrices
    Id, X, Y, Z = getPauliMatrices()
    # Initialize single tensor as complex array and the MPO as an array holding those
    first_tensor = zeros(ComplexF64, 1, 5, 2, 2)
    tensor = zeros(ComplexF64, 5, 5, 2, 2)
    last_tensor = zeros(ComplexF64, 5, 1, 2, 2)
    HeisenbergMPO = MPO{ComplexF64}(undef, N)

    # First MPO tensor
    first_tensor[1,1,:,:] = -lambda * Z
    first_tensor[1,2,:,:] = -J * X
    first_tensor[1,3,:,:] = -J * Y
    first_tensor[1,4,:,:] = -J * Z
    first_tensor[1,5,:,:] = Id

    # Tensors in the center of the MPO
    tensor[1,1,:,:] = Id
    tensor[2,1,:,:] = X
    tensor[3,1,:,:] = Y
    tensor[4,1,:,:] = Z
    tensor[5,1,:,:] = -lambda * Z
    tensor[5,2,:,:] = -J * X
    tensor[5,3,:,:] = -J * Y
    tensor[5,4,:,:] = -J * Z
    tensor[5,5,:,:] = Id

    # Last MPO tensor
    last_tensor[1,1,:,:] = Id
    last_tensor[2,1,:,:] = X
    last_tensor[3,1,:,:] = Y
    last_tensor[4,1,:,:] = Z
    last_tensor[5,1,:,:] = -lambda * Z

    # Fill the MPO
    HeisenbergMPO[1] = first_tensor
    for i = 2:N - 1
        HeisenbergMPO[i] = tensor
    end
    HeisenbergMPO[N] = last_tensor

    return HeisenbergMPO
end

"""
    getIsingMPO(N::Int,J::Real,lambda::Real)

Function to generate the MPO for the Ising Hamilton of the form -J * sum_{i} X^i X^i+1 - lambda * sum_i Z^i.
"""
function getIsingMPO(N::Int, J::Real, lambda::Real)::MPO{Float64}
    # Provide the Pauli matrices
    Id, X, _, Z = getPauliMatrices()
    # Initialize single tensor as complex array and the MPO as an array holding those
    first_tensor = zeros(Float64, 1, 5, 2, 2)
    tensor = zeros(Float64, 5, 5, 2, 2)
    last_tensor = zeros(Float64, 5, 1, 2, 2)
    IsingMPO = MPO{Float64}(undef, N)

    # First MPO tensor
    first_tensor[1,1,:,:] = Id
    first_tensor[1,2,:,:] = -J * X
    first_tensor[1,3,:,:] = -lambda * Z

    # Tensors in the center of the MPO
    tensor[1,1,:,:] = Id
    tensor[1,2,:,:] = -J * X
    tensor[1,3,:,:] = -lambda * Z
    tensor[2,3,:,:] = X
    tensor[3,3,:,:] = Id

    # Last MPO tensor
    last_tensor[1,1,:,:] = -lambda * Z
    last_tensor[2,1,:,:] = X
    last_tensor[3,1,:,:] = Id

    # Fill the MPO
    IsingMPO[1] = first_tensor
    for i = 2:N - 1
        IsingMPO[i] = tensor
    end
    IsingMPO[N] = last_tensor

    return IsingMPO
end

"""
    getTotalSpinMPO(N::Int)

Function to generate the MPO for sum_{i)Z^i.
"""
function getTotalSpinMPO(N::Int)::MPO{Float64}
    # Provide the Pauli matrices
    Id, X, Y, Z = getPauliMatrices()
    # Initialize the result
    TotalSpinMPO = MPO{Float64}(undef, N)

    # Initialize single tensor as complex array and the MPO as an array holding those
    first_tensor = zeros(Float64, 1, 2, 2, 2)
    tensor = zeros(Float64, 2, 2, 2, 2)
    last_tensor = zeros(Float64, 2, 1, 2, 2)

    # First MPO tensor
    first_tensor[1,1,:,:] = Id
    first_tensor[1,2,:,:] = Z

    # Tensors in the center of the MPO
    tensor[1,1,:,:] = Id
    tensor[1,2,:,:] = Z
    tensor[2,2,:,:] = Id

    # Last MPO tensor
    last_tensor[1,1,:,:] = Z
    last_tensor[2,1,:,:] = Id

    # Fill the MPO
    TotalSpinMPO[1] = first_tensor
    for i = 2:N - 1
        TotalSpinMPO[i] = tensor
    end
    TotalSpinMPO[N] = last_tensor

    return TotalSpinMPO
end

"""
    getIdentityMPO(N::Int,d::Int,tensortype::Type{T}=Float64) where T

Function to generate the MPO for the Identity which is useful to check contractions and other functions.
"""
function getIdentityMPO(N::Int, d::Int, tensortype::Type{T}=Float64)::MPO{T} where T
    # Initialize single tensor as complex array and the MPO as an array holding those
    tensor = reshape(Matrix{T}(I, d, d), (1, 1, d, d))
    IdentityMPO = MPO{T}(undef, N)

    # Fill the MPO
    for i = 1:N
        IdentityMPO[i] = tensor
    end

    return IdentityMPO
end


text:
struct HalfEdge
    src::Int
    dst::Int
end
src(he::HalfEdge) = he.src
dst(he::HalfEdge) = he.dst

mutable struct PlanarMultigraph
    v2he::Dict{Int, Int}    # v_id -> he_id
    half_edges::Dict{Int, HalfEdge} # he_id -> he

    f2he::Dict{Int, Int}    # f_id -> he_id
    he2f::Dict{Int, Int}    # he_id -> f_id

    next::Dict{Int, Int}    # he_id -> he_id
    twin::Dict{Int, Int}    # he_id -> he_id

    v_max::Int
    he_max::Int
    f_max::Int
end

vertices(g::PlanarMultigraph) = sort!(collect(keys(g.v2he)))
faces(g::PlanarMultigraph) = sort!(collect(keys(g.f2he)))
half_edges(g::PlanarMultigraph) = sort!(collect(keys(g.half_edges)))

src(g::PlanarMultigraph, he_id::Integer) = src(g.half_edges[he_id])
dst(g::PlanarMultigraph, he_id::Integer) = dst(g.half_edges[he_id])
half_edge(g::PlanarMultigraph, he_id::Integer) = g.half_edges[he_id]
face(g::PlanarMultigraph, he_id::Integer) = g.he2f[he_id]

α(g::PlanarMultigraph, he::Integer) = g.next[he]
next(g::PlanarMultigraph, he::Integer) = α(g, he)
function prev(g::PlanarMultigraph, he::Integer)
    current_he = next(g, he)
    next_he = next(g, current_he)
    while next_he != he
        current_he = next_he
        next_he = next(g, next_he)
    end
    return current_he
end

ϕ(g::PlanarMultigraph, he::Integer) = g.twin[he]
twin(g::PlanarMultigraph, he::Integer) = ϕ(g, he)

σ(g::PlanarMultigraph, he::Integer) = twin(g, prev(g, he))
σ_inv(g::PlanarMultigraph, he::Integer) = next(g, twin(g, he))

nv(g::PlanarMultigraph) = length(g.v2he)
nf(g::PlanarMultigraph) = length(g.f2he)
nhe(g::PlanarMultigraph) = length(g.half_edges)
ne(g::PlanarMultigraph) = nhe(g) ÷ 2

out_half_edge(g::PlanarMultigraph, v::Integer) = g.v2he[v]
surrounding_half_edge(g::PlanarMultigraph, f::Integer) = g.f2he[f]

function trace_orbit(f::Function, a::T; rev::Bool = false) where T
    next = f(a)
    perm = T[a]
    while next != a
        if rev
            pushfirst!(perm, next)
        else
            push!(perm, next)
        end
        next = f(next)
    end
    return perm
end
trace_face(g::PlanarMultigraph, f::Integer) = trace_orbit(h -> g.next[h], surrounding_half_edge(g, f))
trace_vertex(g::PlanarMultigraph, v::Integer) = trace_orbit(h -> σ_inv(g, h), out_half_edge(g, v); rev = true)
neighbors(g::PlanarMultigraph, v::Integer) = [dst(g, he) for he in trace_vertex(g, v)]

is_boundary(g::PlanarMultigraph, he_id::Integer) = (face(g, he_id) == 0)

function rem_vertex!(g::PlanarMultigraph, v::Integer; update::Bool = true)
    for he in trace_vertex(g, v)
        rem_edge!(g, he; update = update)
    end
    delete!(g.v2he, v)
    return g
end

function rem_edge!(g::PlanarMultigraph, he_id::Integer; update::Bool = true)
    # make sure the face of he_id is an inner face
    if is_boundary(g, he_id)
        he_id = twin(g, he_id)
    end

    # handle self-loop
    if next(g, twin(g, he_id)) == twin(g, he_id)
        he_id = twin(g, he_id)
    end

    twin_id = twin(g, he_id)

    if update
        he_next = next(g, he_id)
        he_prev = prev(g, he_id)
        twin_next = next(g, twin_id)
        twin_prev = prev(g, twin_id)

        face_he = face(g, he_id)
        face_twin = face(g, twin_id)

        # remove a inner face
        if face_he != face_twin
            hes_f_he = trace_face(g, face_he)
            rem_face!(g, face_he)
            for he in hes_f_he
                g.he2f[he] = face_twin
            end
        end

        # update f2he
        (surrounding_half_edge(g, face_twin) == twin_id) && (g.f2he[face_twin] = twin_next)

        # update v2he
        (out_half_edge(g, src(g, he_id)) == he_id) && (g.v2he[src(g, he_id)] = twin(g, he_prev))
        (out_half_edge(g, src(g, twin_id)) == twin_id) && (g.v2he[src(g, twin_id)] = twin(g, twin_prev))

        g.next[he_prev] = twin_next
        g.next[twin_prev] = he_next
    end
    delete!(g.next, he_id)
    delete!(g.next, twin_id)
    delete!(g.half_edges, he_id)
    delete!(g.half_edges, twin_id)
    delete!(g.twin, he_id)
    delete!(g.twin, twin_id)
    delete!(g.he2f, he_id)
    delete!(g.he2f, twin_id)

    return g
end

function rem_face!(g::PlanarMultigraph, f::Integer)
    f == 0 && error("Face 0 is for the boundary. It can not be removed.")
    half_edges_f = trace_face(g, f)
    for he in half_edges_f
        delete!(g.he2f, he)
    end
    delete!(g.f2he, f)
    return g
end

function merge_graph!(A::PlanarMultigraph, B::PlanarMultigraph)
    for (v, he_id) in B.v2he
        A.v2he[v + A.v_max] = he_id + A.he_max
    end
    for (he_id, he) in B.half_edges
        A.half_edges[he_id + A.he_max] = HalfEdge(src(he) + A.v_max, dst(he) + A.v_max)
    end
    for (f, he_id) in B.f2he
        if f != 0
            A.f2he[f + A.f_max] = he_id + A.he_max
        end
    end
    for (he_id, f) in B.he2f
        if f != 0
            A.he2f[he_id + A.he_max] = f + A.f_max
        else
            A.he2f[he_id + A.he_max] = 0
        end
    end
    for (curr, next) in B.next
        A.next[curr + A.he_max] = next + A.he_max
    end
    for (curr, twin) in B.twin
        A.twin[curr + A.he_max] = twin + A.he_max
    end

    A.v_max += B.v_max
    A.he_max += B.he_max
    A.f_max += B.f_max

    return A
end

function check_faces(g::PlanarMultigraph)
    for f in faces(g)
        hes_f = trace_face(g, f)
        for he in hes_f
            face(g, he) == f || return false
        end
    end
    return true
end
function check_vertices(g::PlanarMultigraph)
    for v in vertices(g)
        hes_v = trace_vertex(g, v)
        for he in hes_v
            src(g, he) == v || return false
        end
    end
    return true
end
function check_combinatorial_maps(g::PlanarMultigraph)
    for he in half_edges(g)
        (he == α(g, ϕ(g, σ(g, he)))) || return false
    end
    return true
end


text:
using CUDAdrv, CUDAnative

if Sys.iswindows()
	hello_world() =
   	@cuprintf("Greetings from block %lld, thread %lld!\n", Int64(blockIdx().x), Int64(threadIdx().x))
else
	hello_world() =
    	@cuprintf("Greetings from block %ld, thread %ld!\n", Int64(blockIdx().x), Int64(threadIdx().x))
end
@cuda blocks=2 threads=2 hello_world()
synchronize()


text:
using Logging
using SeqLoggers
using Test

@testset "AdvancedLogger" begin

file_name_pattern = raw"\a\c\c\e\s\s-YYYY-mm-dd-HH-MM.\l\o\g"
dir_path = mktempdir()
logger =  AdvancedFileLogger(dir_path, file_name_pattern; min_level=Logging.Debug)
with_logger(logger) do
    @debug "Log to file"
    @info "Log to file"
    @warn "Log to file"
    @error "Log to file"
    sleep(61)
    @info "Log to next file" A = 2
end
log_files = joinpath.(dir_path, readdir(dir_path)) #kwarg join=true no in Julia 1.3
@test length(log_files) == 2

content = """[Debug] Log to file
[Info] Log to file
[Warning] Log to file
[Error] Log to file
"""
@test log_files[1] |> read |> String == content

@test_throws ArgumentError run_with_logger(logger) do
    @debug "Debug Log Message with 1 kwargs" keywarg="A"
    @info "Info Log Message with 2 kwargs" keywarg="A" second=rand(200, 10) third=DataFrame(rand(20,10))
    @warn "Warn Log message without kwargs"
    @error "Error Log Before Exception" A="A" B="AB" C="ABC" D="ABCD"
    throw(ArgumentError("Throw an exception"))
    @info "After Exception" # Will not be excecuted
end

end

@testset "Load AdvancedFileLogger" begin

logging_config = Dict(
    "min_level" => "INFO",
    "dir_path" => "C:\\temp",
    "file_name_pattern" => "\\t\\e\\s\\t_YYYY-mm-dd.\\l\\o\\g",
)

logger = SeqLoggers.load_advanced_filelogger(logging_config)
@test logger.min_level == Logging.Info

logging_config = Dict(
    "min_level" => "INFO",
    "dir_path" => "C:\\temp",
    "file_name_pattern" => "\\t\\e\\s\\t_YYYY-mm-dd.\\l\\o\\g",
    "transformation" => "add_timestamp",
)

tran_logger = SeqLoggers.load_advanced_filelogger(logging_config)
@test tran_logger.logger.min_level == Logging.Info

end

@testset "Load logger from config" begin

config = Dict(
    "logging" => [
        "SeqLogger" => Dict(
            "server_url" => "http://subdn215:5341/",
            "min_level" => "INFO",
        ),
        "ConsoleLogger" => Dict(
            "min_level" => "INFO",
        ),
        "FileLogger" => Dict(
            "min_level" => "INFO",
            "file_path" => raw"C:\Temp\test.txt",
            "append" => true,
        ),
        "AdvancedFileLogger" => Dict(
            "min_level" => "INFO",
            "dir_path" => "C:\\temp",
            "file_name_pattern" => "\\t\\e\\s\\t_YYYY-mm-dd.\\l\\o\\g",
        )
    ]
)

@test SeqLoggers.get_logger(config["logging"][4]...) isa MinLevelLogger # of FileLogger

tee_logger = SeqLoggers.load_logger_from_config(config)
@test tee_logger.loggers[4] isa MinLevelLogger # of FileLogger

file_path = joinpath(@__DIR__, "data", "config_with_advanced_logger.json")
tee_logger = SeqLoggers.load_logger_from_config(file_path)
@test tee_logger.loggers[4] isa MinLevelLogger # of FileLogger

end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(145)
    x, y, yerr = 1:2:20, 5*rand(10), 0.4*abs.(randn(10))
    fig = Figure(resolution = (700, 450), font = "sans")
    ax = Axis(fig, xlabel = "variables", ylabel = "values")
    barplot!(ax, x,y,strokewidth = 1,color = x,colormap = (:Spectral_10, 0.85),
            strokecolor = :black)
    errorbars!(ax, x, y, yerr, whiskerwidth = 12)
    fig[1,1] = ax
    fig
    save(joinpath(@__DIR__, "output", "BoxErrorBarsCmap.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using CSV, DataFrames, Test
using FixedEffectModels
using FixedEffectModels: parse_fixedeffect, _parse_fixedeffect, _multiply
using FixedEffects
import Base: ==

==(x::FixedEffect{R,I}, y::FixedEffect{R,I}) where {R,I} =
    x.refs == y.refs && x.interaction == y.interaction && x.n == y.n

@testset "parse_fixedeffect" begin
    csvfile = CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv"))
    df = DataFrame(csvfile)
    # Any table type supporting the Tables.jl interface should work
    for data in [df, csvfile]
    	@test _parse_fixedeffect(data, term(:Price)) === nothing
        @test _parse_fixedeffect(data, ConstantTerm(1)) === nothing
        @test _parse_fixedeffect(data, fe(:State)) == (FixedEffect(data.State), :fe_State)

    	@test _parse_fixedeffect(data, fe(:State)&term(:Year)) ==
            (FixedEffect(data.State, interaction=_multiply(data, [:Year])), Symbol("fe_State&Year"))
        @test _parse_fixedeffect(data, fe(:State)&fe(:Year)) ==
            (FixedEffect(data.State, data.Year), Symbol("fe_State&fe_Year"))

        @test parse_fixedeffect(data, ()) == (FixedEffect[], Symbol[], ())

        f = @formula(y ~ 1 + Price)
        ts1 = f.rhs
        ts2 = term(1) + term(:Price)
        @test parse_fixedeffect(data, f) == (FixedEffect[], Symbol[], f)
        @test parse_fixedeffect(data, ts1) == (FixedEffect[], Symbol[], ts1)
        @test parse_fixedeffect(data, ts2) == parse_fixedeffect(data, ts1)

        fparsed = term(:y) ~ InterceptTerm{false}() + term(:Price)
        tsparsed = (InterceptTerm{false}(), term(:Price))

        f = @formula(y ~ 1 + Price + fe(State))
        ts1 = f.rhs
        ts2 = term(1) + term(:Price) + fe(:State)
        @test parse_fixedeffect(data, f) == ([FixedEffect(data.State)], [:fe_State], fparsed)
        @test parse_fixedeffect(data, ts1) == ([FixedEffect(data.State)], [:fe_State], tsparsed)
        @test parse_fixedeffect(data, ts2) == parse_fixedeffect(data, ts1)

        f = @formula(y ~ Price + fe(State) + fe(Year))
        ts1 = f.rhs
        ts2 = term(:Price) + fe(:State) + fe(:Year)
        @test parse_fixedeffect(data, f) == ([FixedEffect(data.State), FixedEffect(data.Year)], [:fe_State, :fe_Year], fparsed)
        @test parse_fixedeffect(data, ts1) == ([FixedEffect(data.State), FixedEffect(data.Year)], [:fe_State, :fe_Year], tsparsed)
        @test parse_fixedeffect(data, ts2) == parse_fixedeffect(data, ts1)

        f = @formula(y ~ Price + fe(State)&Year)
        ts1 = f.rhs
        ts2 = term(:Price) + fe(:State)&term(:Year)
        @test parse_fixedeffect(data, f) == ([FixedEffect(data.State, interaction=_multiply(data, [:Year]))], [Symbol("fe_State&Year")], term(:y) ~ (term(:Price),))
        @test parse_fixedeffect(data, ts1) == ([FixedEffect(data.State, interaction=_multiply(data, [:Year]))], [Symbol("fe_State&Year")], (term(:Price),))
        @test parse_fixedeffect(data, ts2) == parse_fixedeffect(data, ts1)

        f = @formula(y ~ Price + fe(State)*fe(Year))
        ts1 = f.rhs
        ts2 = term(:Price) + fe(:State) + fe(:Year) + fe(:State)&fe(:Year)
        @test parse_fixedeffect(data, f) == ([FixedEffect(data.State), FixedEffect(data.Year), FixedEffect(data.State, data.Year)], [:fe_State, :fe_Year, Symbol("fe_State&fe_Year")], fparsed)
        @test parse_fixedeffect(data, ts1) == ([FixedEffect(data.State), FixedEffect(data.Year), FixedEffect(data.State, data.Year)], [:fe_State, :fe_Year, Symbol("fe_State&fe_Year")], tsparsed)
        @test parse_fixedeffect(data, ts2) == parse_fixedeffect(data, ts1)
    end
end


text:
using NESSie.Format
using NESSie.Format: readmcsf_nodes, readmcsf_elements

testfiles = ((mktemp()..., (0,0)), # empty file
             (mktemp()..., (4,2))) # dummy file
write(testfiles[2][2], """
mcsf_begin=1;

      dim = 3;
    dimii = 3;
 vertices = 4;
simplices = 2;
vert=[
%------------------------------------------------------------------------------
%  Node-ID  Chrt        X-Coordinate        Y-coordinate        Z-coordinate
%---------  ----     ----------------     ----------------     ----------------
        0    2      0.0000000000e+00     0.0000000000e+00     0.0000000000e+00
        1    2      1.0000000000e+00     0.0000000000e+00     0.0000000000e+00
        2    2      0.0000000000e+00     1.0000000000e+00     0.0000000000e+00
        3    2      0.0000000000e+00     0.0000000000e+00     1.0000000000e+00
];
simp=[
%------------------------------------------------------------------------------------------
%  Simp-ID Grp    Mat          Face-Types                      Vertex-Numbers
%--------- ---    ---    ---------------------  -------------------------------------------
        0  0      0      0     0     0     0           0          1          2          3
        1  0      1      0     0     0     0           3          2          1          0
];
mcsf_end=1;""")

@testset "readmcsf" begin
    try
        for (fname, fh, len) in testfiles, fn in ("single", "bulk"), T in testtypes
            seekstart(fh)
            nodes = (); elements = ()
            if fn == "single"
                # check single methods
                nodes = readmcsf_nodes(fh, T)
                elements = readmcsf_elements(fh, nodes)
            elseif fn == "bulk"
                # check bulk method
                model = readmcsf(fh, T)
                nodes, elements = (model.nodes, model.elements)
            end
            # check return types
            @test typeof(nodes) == Vector{Vector{T}}
            @test typeof(elements) == Vector{Tetrahedron{T}}
            # check lengths
            @test length(nodes) == len[1]
            @test length(elements) == len[2]

            # check content
            if fname == testfiles[2][1]
                @test nodes[1] == T[0, 0, 0]
                @test nodes[2] == T[1, 0, 0]
                @test nodes[3] == T[0, 1, 0]
                @test nodes[4] == T[0, 0, 1]
                @test elements[1].v1 === nodes[1]
                @test elements[1].v2 === nodes[2]
                @test elements[1].v3 === nodes[3]
                @test elements[1].v4 === nodes[4]
                @test elements[2].v1 === nodes[4]
                @test elements[2].v2 === nodes[3]
                @test elements[2].v3 === nodes[2]
                @test elements[2].v4 === nodes[1]
            end
        end
    finally
        [begin close(fh); rm(fn); end for (fn, fh, _) in testfiles]
    end
end


text:
@testset "$TEST $G" begin
# TODO

end # testset


text:
module Swagger

using HTTP
using JSON
using MbedTLS
using Dates

import Base: convert, show, summary, getindex, keys, length
import JSON: lower

@static if VERSION < v"0.7.0-"
    import Base: start, done, next
else
    import Base: iterate
end

abstract type SwaggerModel end
abstract type SwaggerApi end

include("client.jl")
include("json.jl")
include("val.jl")

end # module Swagger


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

using Markdown
using InteractiveUtils

# ╔═╡ 8c1ba773-e9ca-4269-9519-d893b4e062e8
begin
	let
		# Set up temporary environment
		env = mktempdir()
		import Pkg
		Pkg.activate(env)
		Pkg.Registry.update()

		# Download packages
		Pkg.add("PlutoUI")
		Pkg.add("DICOM")
		Pkg.add(Pkg.PackageSpec(url="https://github.com/Dale-Black/DICOMUtils.jl"))
	end

	# Import packages
	using PlutoUI
	using DICOM
	using DICOMUtils
end

# ╔═╡ bae6661b-1e81-4006-b52a-76938bfafa2a
md"""
## Set up
"""

# ╔═╡ dcf61796-1f8b-414c-9854-5b62ae653ddb
TableOfContents()

# ╔═╡ 2d9983f2-b9ae-424d-bb0e-90ca7b89aee7
filepath_new = raw"Y:\Canon Images for Dynamic Heart Phantom\Dynamic Phantom\clean_data";

# ╔═╡ 6ba2428e-9507-4b74-a815-0a59a4d5feb2
isdir(filepath_new) || mkdir(filepath_new)

# ╔═╡ 0d171ffa-e5f0-4148-b553-e4fa240a91e5
md"""
## 1. Organize by PatientName
Not needed, already completed
"""

# ╔═╡ 5765c970-67b6-4f58-9b38-66dca0681d03
filepath = raw"Y:\Canon Images for Dynamic Heart Phantom\Dynamic Phantom\DD2710";

# ╔═╡ 8f2f872e-2916-4bc5-bd72-1ece24e5c8a7
const PatientName = (0x0010,0x0010)

# ╔═╡ edaf8e6f-edcd-43c6-8c12-de5839a95bf9
# DICOMUtils.sortbytag_move(filepath, filepath_new, PatientName)

# ╔═╡ 2f3f06fe-7db4-42c3-b84e-d7a81da85983
md"""
## 2. Organize by SeriesNumber
"""

# ╔═╡ d13b882a-e91b-451e-bec0-e0e7b1aff6de
const SeriesNumber = (0x0020,0x0011)

# ╔═╡ 7503aac9-f9c2-4c0d-ab1c-42088e82211e
config1_path = raw"Y:\Canon Images for Dynamic Heart Phantom\Dynamic Phantom\clean_data\CONFIG 4^275";

# ╔═╡ 79b490af-264c-4bb8-832f-ca38795ba002
# DICOMUtils.sortbytag_move(config1_path, config1_path, SeriesNumber)

# ╔═╡ 8cb32b19-5d78-4327-af3e-b00d279b55c3
md"""
## 3a. Organize by kVp
This should be done for the static (0 bpm) helical scan folders
"""

# ╔═╡ 0c1379e4-9576-459a-94cb-ad8f190e122e
const kVp = (0x0018,0x0060)

# ╔═╡ a1a24093-ae31-42ec-a8ca-adf197c3232a
kv_path = raw"Y:\Canon Images for Dynamic Heart Phantom\Dynamic Phantom\clean_data\CONFIG 4^275\52"

# ╔═╡ a1fab547-777d-4511-94ea-5a9966c88cea
# DICOMUtils.sortbytag_move(kv_path, kv_path, kVp)

# ╔═╡ 5500af3a-06b3-47c1-8594-6ac5aeff6bc4
md"""
## 3b. Organize by ImageComments
This should be done for the dynamic (60, 80, 100 bpm) helical scan folders
"""

# ╔═╡ ae75f4be-82d9-4a0f-9430-2b72e8c58525
const ImageComments =  (0x0020,0x4000)

# ╔═╡ 71fb31d3-a204-49ac-9683-4265058cd014
hel_path = raw"Y:\Canon Images for Dynamic Heart Phantom\Dynamic Phantom\clean_data\CONFIG 4^350\51"

# ╔═╡ 4563ad8d-e2f3-4cc0-a95f-491d40be2408
 DICOMUtils.sortbytag_move(hel_path, hel_path, ImageComments)

# ╔═╡ 238ee3f5-3f1b-432d-acd2-29035c0a0791
md"""
## (Extra) Sort through every directory at once
"""

# ╔═╡ 3eb6343a-1f78-47ac-96f8-e2bfb82d7b2b
head_dir = raw"Y:\Canon Images for Dynamic Heart Phantom\Dynamic Phantom"

# ╔═╡ 85f3fcd7-d644-4fe0-9233-3374aa2b9bdb
dirs = readdir(head_dir)

# ╔═╡ abc8e616-20ce-487d-aa56-028fe4a5f9e3
head_dir * "/" * dirs[1]

# ╔═╡ 2be55e24-1478-4e94-a02d-71b999c52117
# for i = 1:(length(dirs) - 1)
# 	DICOMUtils.sortbytag_move(head_dir * "/" * dirs[i], filepath_new, PatientName)
# end

# ╔═╡ Cell order:
# ╟─bae6661b-1e81-4006-b52a-76938bfafa2a
# ╠═8c1ba773-e9ca-4269-9519-d893b4e062e8
# ╟─dcf61796-1f8b-414c-9854-5b62ae653ddb
# ╠═2d9983f2-b9ae-424d-bb0e-90ca7b89aee7
# ╠═6ba2428e-9507-4b74-a815-0a59a4d5feb2
# ╟─0d171ffa-e5f0-4148-b553-e4fa240a91e5
# ╠═5765c970-67b6-4f58-9b38-66dca0681d03
# ╠═8f2f872e-2916-4bc5-bd72-1ece24e5c8a7
# ╠═edaf8e6f-edcd-43c6-8c12-de5839a95bf9
# ╟─2f3f06fe-7db4-42c3-b84e-d7a81da85983
# ╠═d13b882a-e91b-451e-bec0-e0e7b1aff6de
# ╠═7503aac9-f9c2-4c0d-ab1c-42088e82211e
# ╠═79b490af-264c-4bb8-832f-ca38795ba002
# ╟─8cb32b19-5d78-4327-af3e-b00d279b55c3
# ╠═0c1379e4-9576-459a-94cb-ad8f190e122e
# ╠═a1a24093-ae31-42ec-a8ca-adf197c3232a
# ╠═a1fab547-777d-4511-94ea-5a9966c88cea
# ╟─5500af3a-06b3-47c1-8594-6ac5aeff6bc4
# ╠═ae75f4be-82d9-4a0f-9430-2b72e8c58525
# ╠═71fb31d3-a204-49ac-9683-4265058cd014
# ╠═4563ad8d-e2f3-4cc0-a95f-491d40be2408
# ╟─238ee3f5-3f1b-432d-acd2-29035c0a0791
# ╠═3eb6343a-1f78-47ac-96f8-e2bfb82d7b2b
# ╠═85f3fcd7-d644-4fe0-9233-3374aa2b9bdb
# ╠═abc8e616-20ce-487d-aa56-028fe4a5f9e3
# ╠═2be55e24-1478-4e94-a02d-71b999c52117


text:
using LinearAlgebra
using Images
using Printf

# load data
xs = load("./examples/ChessBoard/inpainting_sparse.png")
xr = load("./examples/ChessBoard/inpainting_low_rank.png")

# size
m, n = size(xs)

# preprocessing on xs
xs = Gray.(xs); xs = convert(Array{Float64}, xs)
γ = xs[2,2]; xs = xs .- γ # normalize
idx = findall(x->abs(x)<=5e-2, xs); xs[idx] .= 0
idx = findall(x->x!=0, xs); xs[idx] .= maximum(xs)
ks = count(x->x!=0, xs)
xs = 3*vec(xs)

# preprocessing on xr
xr = Gray.(xr); xr = convert(Array{Float64}, xr)
F = svd(xr); kr = 2
xr = F.U[:, 1:kr] * Diagonal(F.S[1:kr]) * F.Vt[1:kr, :]
xr = 0.5*vec(xr)

# observation
b = xs + xr

# write to file
f = open("/home/zhenan/Github/AtomicOpt.jl/examples/ChessBoard/ChessBoardData.txt", "w")
l = length(b)
@printf(f, "xs = ")
for i in 1:l
    @printf(f, "%.3f ", xs[i])
end
println(f, "\n")
@printf(f, "xr = ")
for i in 1:l
    @printf(f, "%.3f ", xr[i])
end
println(f, "\n")
@printf(f, "b = ")
for i in 1:l
    @printf(f, "%.3f ", b[i])
end
println(f, "\n")
@printf(f, "ks = %d\n", ks)
@printf(f, "kr = %d\n", kr)
@printf(f, "m = %d\n", m)
@printf(f, "n = %d\n", n)
close(f)


text:
mutable struct Leadfield{T<:AbstractFloat,S<:AbstractString}
    L::Array{T,3}
    x::Vector{T}
    y::Vector{T}
    z::Vector{T}
    sensors::Vector{S}
end

function Base.show(io::IO, l::Leadfield)
    @printf "Leadfield\n"
    @printf "  Number of sources: %d\n" size(l.L, 1)
    @printf "  Number of dimensions + nulls: %d\n" size(l.L, 2)
    @printf "  Number of sensors: %d\n" size(l.L, 3)
end

function match_leadfield(l::Leadfield, s::SSR)

    @info("Matching leadfield to SSR")

    idx = [something(findfirst(isequal(name), l.sensors), 0) for name in channelnames(s)]

    if length(unique(idx)) < length(idx)
        error(
            "Not all SSR channels mapped to sensor #SSR=$(length(channelnames(s))), #L=$(length(l.sensors))",
        )
    end

    l.L = l.L[:, :, idx]
    l.sensors = l.sensors[idx]

    @debug("matched $(length(idx)) sensors")

    return l
end


text:
module Testing
    export run_test

    include("Utils.jl")
    include("Spa.jl")

    struct SynthesizedSignal
        f::Real
        Fs::Int
        amplitudes::Array{Real}
    end

    struct Test
        f::Real
        signal::Array{Real}
    end

    function training_signal(f, Fs; amplitudes = .01 .* rand(50:100, num_harmonics))
        return SynthesizedSignal(f, Fs, amplitudes)
    end

    function get_signal_step(signal::SynthesizedSignal, n::Int)
        return sum([sin(2*π*n*i * signal.f / signal.Fs) * signal.amplitudes[i] / i for i in 1:length(signal.amplitudes)])
    end

    function generate_test(f, length::Int, Fs; num_harmonics = 20, amplitudes = 0.01 .* rand(50:100, num_harmonics))
        test = zeros(length)
        signal = training_signal(f, Fs, amplitudes = amplitudes)
        for n = 0:length-1
            test[n+1] = get_signal_step(signal, n)
        end
        return Test(f, test)
    end

    function get_error(guess, answer)
        return abs(answer - guess) / answer
    end

    function run_test(detect_pitch::Function; N = 10000, Fs = 48000, Wl = 2048, low_note = 80, high_note = Wl / 2)
        note_range = Utils.getRange(low_note, high_note)

        err = zeros(N)
        Threads.@threads for n in 1:N
            test = generate_test(note_range[rand(1:length(note_range))], Wl, Fs)
            guess, _ = detect_pitch(test.signal, Fs)

            err[n] = get_error(guess, test.f)
        end

        return sum(err) / N
    end
end


text:
addprocs(3)

@everywhere using JuLIP, JuLIP.ASE
@everywhere using TightBinding
# using ProfileView

# test parameters
beta = 20.0        # temperature / smearing paramter
                  # 10 to 50 for room temperature
n0 = 1            # site index where we compute the site energy
NQUAD = 8         # number of contour points
DIM = (7,7,7)

TB=TightBinding

# define the model
tbm = TB.NRLTB.NRLTBModel(elem = TB.NRLTB.Si_sp, nkpoints = (0,0,0))
calc = TB.Contour.ContourCalculator(tbm, NQUAD)
# now the real system to test on
at = DIM * bulk("Si", pbc=true, cubic=true)
JuLIP.rattle!(at, 0.02)
TB.Contour.calibrate2!(calc, at, beta, nkpoints=(6,6,6))
@show length(at)
Nat = length(at)
Isub = 1:(length(at) ÷ 3)

@time Es, dEs = TB.Contour.site_energy(calc, at, n0, true)
@time Es, dEs = TB.Contour.site_energy(calc, at, n0, true)
@time Es, dEs = TB.Contour.partial_energy(calc, at, Isub, true)
@time Es, dEs = TB.Contour.partial_energy(calc, at, Isub, true);

# MEMORY:
#   216 >>>   103,514,112   ~ 100MB
#  1000 >>>   479,232,000   ~ 500MB
#  2744 >>> 1,315,012,608   ~ 1.3GBP
#  5000 >>> 2,500,000,000   ~ 2.5 GB  (guessed)
# 10000 >>>  ~ 5GB   (ca 21^3 cell)
#    i.e. a 41^3 ~ 100k atoms cell will lead to ~ 17 GB, which is still OK!
#    on Apple Laopto: 16 GB  i.e. up to 10k atoms
#    on Kaiju we have 32,897,748k = 32GB   i.e. up to 100k atoms
#
# storage for neighbourlist:
#     4_096_000_000 atoms  >>> not a chance

# @code_warntype TB.Contour.site_energy(calc, at, n0, true);
# Es, dEs = TB.Contour.site_energy(calc, at, n0, true)
# ;

# TESTING HOW MUCH STORAGE WOULD BE REQUIRED TO KEEP THE RESIDUALS!
# PROBABLY TOO MUCH!
# for d in (3, 5, 7)
#    at = (d,d,d) * bulk("Si", pbc=(false,false,false), cubic=true)
#    Isub = 1:(length(at) ÷ 3)
#    @show length(at)
#    Es, dEs = TB.Contour.partial_energy(calc, at, Isub, true)
# end
# quit()

# cross-over between full and sparse seems to be
# around (1,10,10) in 2D and  smaller in 3D

# for n in (4,8,16)
#    @show n
#
#    # now go to the real system
#    at = (1,n,n) * bulk("Si", pbc=(true,true,true), cubic=false)
#    JuLIP.rattle!(at, 0.02)
#
#    H, M = TB.hamiltonian(tbm, at)
#    Hf = full(H); Mf = full(M)
#
#    rhs = rand(size(H,1), size(H,1) ÷ 3)
#
#    println("sparse")
#    @time LU=lufact(H - 0.1im * M)
#    @time LU=lufact(H - 0.1im * M)
#    @time LU \ rhs
#    @time LU \ rhs
#    println("full")
#    @time LU=lufact(Hf - 0.1im * Mf)
#    @time LU=lufact(Hf - 0.1im * Mf)
#    @time LU \ rhs
#    @time LU \ rhs
# end

# for n in (4,6,8)
#    @show n
#
#    # now go to the real system
#    at = (n,n,n) * bulk("Si", pbc=(true,true,true), cubic=false)
#    JuLIP.rattle!(at, 0.02)
#
#    H, M = TB.hamiltonian(tbm, at)
#    Hf = full(H); Mf = full(M)
#
#    rhs = rand(size(H,1))
#
#    println("sparse")
#    @time LU=lufact(H - 0.1im * M)
#    @time LU=lufact(H - 0.1im * M)
#    @time LU \ rhs
#    @time LU \ rhs
#    println("full")
#    @time LU=lufact(Hf - 0.1im * Mf)
#    @time LU=lufact(Hf - 0.1im * Mf)
#    @time LU \ rhs
#    @time LU \ rhs
# end

# # timing test
# println("timing with nquad = ", calc.nquad, "  (ca. 6 digits)")
# @time TB.Contour.site_energy(calc, at, n0)
# @time TB.Contour.site_energy(calc, at, n0)
# @show length(at) * tbm.norbitals


text:
using Revise
using NLS_Fit
using DelimitedFiles
using LinearAlgebra

rootDir = joinpath(dirname(pathof(NLS_Fit)), "..")
dataDir = joinpath(rootDir,"data")
dataFile = joinpath(dataDir,"mapped_parameters.txt")

# Getting started example ****************
#

# Data ================
#
XY=readdlm(dataFile);
X=XY[:,1];
Y=XY[:,2];

# # Model
# #

# # σ(μ)
# map_μ_to_σ = Map_Affine(1.0=>1.0,30.0=>5.0)

# # map (μ,σ(μ))
# struct Varying_σ_Map <: NLS_Fit.Abstract_Map
#     _n_peak::Int
#     _map_μ_to_σ::NLS_Fit.Abstract_Map
# end

# NLS_Fit.parameter_size(varying_σ_map::Varying_σ_Map) =
#     2*varying_σ_map._n_peak + NLS_Fit.parameter_size(varying_σ_map._map_μ_to_σ)

# function NLS_Fit.eval_map(varying_σ_map::Varying_σ_Map,hat_θ::AbstractVector)
#     @assert parameter_size(varying_σ_map) == length(hat_θ)

#     n_peak = varying_σ_map._n_peak
#     map_μ_to_σ = varying_σ_map._map_μ_to_σ

#     input_θ =  @view hat_θ[1:(2*varying_σ_map._n_peak)]
#     map_μ_to_σ_θ = @view hat_θ[(2*varying_σ_map._n_peak + 1):end]

#     θ = Vector{eltype(hat_θ)}(undef,3*varying_σ_map._n_peak)

#     # create input views
#     view_h = @view input_θ[1:2:2*n_peak]
#     view_μ = @view input_θ[2:2:2*n_peak]

#     # copy height, μ and σ(μ)
#     θ[1:3:3*n_peak] .= view_h
#     θ[2:3:3*n_peak] .= view_μ
#     θ[3:3:3*n_peak] .= eval_map(map_μ_to_σ,view_μ,map_μ_to_σ_θ)

#     θ
# end

# t = Varying_σ_Map(4,map_μ_to_σ)
# parameter_size(t)
# eval_map(t,ones(parameter_size(t)))

# model = Gaussian_Peak() + Gaussian_Peak() + Gaussian_Peak()
# model = Model2Fit_Transformed_Parameters(model,Varying_σ_Map(3,map_μ_to_σ))

# parameter_size(model)
# θ_init = ones(8)
# θ_init[2]=4
# θ_init[4]=8
# θ_init[6]=25

# nls = NLS_ForwardDiff_From_Model2Fit(model,X,Y)
# conf = NLS_Solver.LevenbergMarquardt_Conf()
# result = NLS_Solver.solve(nls,θ_init,conf)

# @assert NLS_Solver.converged(result)
# θ_fit = NLS_Solver.solution(result)
# get_model(model)
# get_model_θ(model,θ_fit)

# Y_init = eval_y(model,X,θ_init)
# Y_fit = eval_y(model,X,θ_fit)

# writedlm("test.txt",hcat(X,Y,Y_init,Y_fit))

#
source=collect(1:3:12)
dest=collect(2:3:12)

function src_indices_after_dest_removed(source_dest::Pair{SOURCE,DEST}) where {SOURCE<:AbstractVector{Int},DEST<:AbstractVector{Int}}
    (src,dest) = source_dest

    @assert length(src)==length(dest)

    n = length(src)

    count =  zeros(Int,n)
    i_count = 1

    # note this works even if src, dest are not sorted but this has
    # n^2 complexity
    for src_i in src
        for dest_i in dest
            if dest_i < src_i
                count[i_count] += 1
            end
        end
        i_count += 1
    end
    src - count
end
src_indices_after_dest_removed(dest=>dest)

src_indices_after_dest_removed(source=>dest)

# An helper that  maps
#
# [ θ̂_1, θ̂_2, ... , θ̂_n, θ̂f_1, θ̂f_2, ..., θ̂f_nf ]
#
# to
#
# [ θ_1, θ_2, ... , θ_{n+m} ] (m is dim of src = dim of dest)
#
# where θ̂_i for i ∈ src have be taken, and f_{θf}(θ̂_i) inserted at pos dest
#
struct Transform_Src_Insert_Dest_Map <: NLS_Fit.Abstract_Map
    _map::NLS_Fit.Abstract_Map
    _src_dest::Matrix{Int} # src=col(1), dest=col(2)
end

function Transform_Src_Insert_Dest_Map(map::NLS_Fit.Abstract_Map, source_dest::Pair{SOURCE,DEST}) where {SOURCE<:AbstractVector{Int},DEST<:AbstractVector{Int}}
    @assert length(first(source_dest))==length(last(source_dest))

    src = src_indices_after_dest_removed(source_dest)
    dest = last(source_dest)

    Transform_Src_Insert_Dest_Map(map,hcat(src,dest))
end

NLS_Fit.parameter_size(tm::Transform_Src_Insert_Dest_Map) = NLS_Fit.parameter_size(tm._map)

function NLS_Fit.eval_map(tm::Transform_Src_Insert_Dest_Map,hat_model_θ::AbstractVector,hat_map_θ::AbstractVector)
    @assert parameter_size(tm) == length(hat_map_θ)

    src  = @view tm._src_dest[:,1]
    dest = @view tm._src_dest[:,2]

    θ_mapped = eval_map(tm._map,(@view hat_model_θ[src]),hat_map_θ)

    θ = insert_some_elements(hat_model_θ, dest, θ_mapped)

    θ
end

# Demo ****************
#
model = Gaussian_Peak() + Gaussian_Peak() + Gaussian_Peak()
n = NLS_Fit.parameter_size(model)

map_μ_to_σ = Map_Affine(1.0=>1.0,30.0=>5.0)
source=collect(2:3:n)
dest=collect(3:3:n)

map = Transform_Src_Insert_Dest_Map(map_μ_to_σ,source=>dest)
model = Model2Fit_Transformed_Parameters(model,n-length(dest),map)

NLS_Fit.parameter_size(model)
NLS_Fit.get_model_θ(model,Float64[1:8;])

θ_init = ones(8)
θ_init[2]=4
θ_init[4]=8
θ_init[6]=25

nls = NLS_ForwardDiff_From_Model2Fit(model,X,Y)
conf = NLS_Solver.LevenbergMarquardt_Conf()
result = NLS_Solver.solve(nls,θ_init,conf)

@assert NLS_Solver.converged(result)
θ_fit = NLS_Solver.solution(result)
get_model(model)
get_model_θ(model,θ_fit)

Y_init = eval_y(model,X,θ_init)
Y_fit = eval_y(model,X,θ_fit)

writedlm("test.txt",hcat(X,Y,Y_init,Y_fit))

# parameter_size(t)
# hat_model_θ = Float64[1:12;]
# hat_model_θ = deleteat!(hat_model_θ,dest)
# hat_map_θ = ones(2)
# NLS_Fit.eval_map(t,hat_model_θ,hat_map_θ)


text:
module Mozi

include("./FESparse.jl")
include("./assembly/Enums.jl")
include("./assembly/CoordinateSystem.jl")
include("./structure/FEStructure.jl")
include("./load/LoadCase.jl")
include("./assembly/Assembly.jl")

include("./solver/Solver.jl")
include("./result/Result.jl")

using .FEStructure
using .LoadCase
using .FEAssembly
using .Solver
using .Result

using Serialization

export save,load

function save(assembly)
    open("./test_assembly.mz","w+") do f
        serialize(f,assembly)
    end
end

function load()
    assembly=nothing
    open("./test_assembly.mz") do f
        assembly=deserialize(f)
    end
    return assembly
end

export Structure,LoadCaseSet,Assembly

export add_uniaxial_metal!

export add_general_section!,add_beam_section!

export add_node!,set_nodal_mass!,set_nodal_restraint!,set_nodal_spring!

export add_beam!,set_beam_release!,set_beam_orient!

export add_quad!,set_quad_rotation!

export add_tria!,set_tria_rotation!

export add_static_case!, add_modal_case!, add_buckling_case!, add_time_history_case!

export add_nodal_force!

export add_beam_strain!, add_beam_distributed!

export set_modal_params!, set_buckling_params!,set_time_history_params!

export assemble!,clear_result!

export solve

export result_nodal_displacement,result_nodal_reaction,result_nodal_time_history

export result_beam_force,result_beam_displacement

export result_quad_force,result_quad_displacement

export result_modal_period,result_eigen_value, result_eigen_vector

end


text:
using PersistenceDiagrams
using Aqua

Aqua.test_all(PersistenceDiagrams; ambiguities=false)


text:
include("../utils.jl")

@testset "NotImplementedError" begin
    e = G.NotImplementedError("foo")
    @test e.msg // "[foo] hasn't been implemented."
end

@testset "UnknownOptionError" begin
    e = G.UnknownOptionError("foo", 5)
    @test e.msg // "[foo] is not recognised as a valid option name for Int64."
end

@testset "OptionValueError" begin
    e = G.OptionValueError("foo", 5, "Bool")
    @test e.msg // "[foo] value given (5) didn't meet the expected format Bool."
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

@userplot ReturnPlot

@recipe function f(rp::ReturnPlot)
  # get user input
  obj = rp.args[1]

  if obj isa AbstractVector
    seriestype --> :scatter
    δt, ms = returnlevels(obj)
  elseif obj isa GeneralizedExtremeValue
    seriestype --> :path
    mmin, mmax = rp.args[2:3]
    δt, ms = returnlevels(obj, mmin, mmax)
  end

  xscale --> :log10
  xlabel --> "return period"
  ylabel --> "return level"
  label  --> "return plot"

  δt, ms
end


text:
module MCMC

using Distributions
using LinearAlgebra
include("PropCov.jl")

function _metropolis(x::Float64, logprob::Function, propsd::Float64)
  cand = curr + randn() * propsd
  met_ratio = logprob(cand) - logprob(curr)
  accept = met_ratio > log(rand())
  draw = accept ? cand : curr
  return (draw, accept)
end

function metropolis(x::Float64, logprob::Function, propsd::Float64)
  return _metropolis(x, logprob, propsd)[1]
end

function metropolis(curr::Vector{Float64}, logprob::Function,
                    propcov::Matrix{Float64})
  cand = rand(MvNormal(curr, propcov))
  met_ratio = logprob(cand) - logprob(curr)
  out = met_ratio > log(rand()) ? cand : curr
  return out
end

function mcmc(init::Vector{T}, ll::Function, lp::Function;
              nmcmc::Int, nburn::Int, batchsize::Int=50,
              thin::Int=1,
              propcov=nothing,
              propcov_init=eye(length(init)),
              temper_power=1,
              momentum=0.1, min_temper=1.0, max_temper=1.0,
              verbose=0) where T <: AbstractFloat
  @assert min_temper <= max_temper

  # Initialize MCMC
  state = deepcopy(init)

  # Number of parameters
  nparam = length(state)

  # Preallocate memory
  out = zeros(nmcmc, nparam)

  # Log density of target distribution, for a given temperature
  logprob(s::Vector{Float64}, temper::Float64) = ll(s) / temper + lp(s)

  # Total number of iterations
  total_iters = nmcmc + nburn

  # proposal covariance
  if propcov == nothing
    propcov = PropCov(propcov_init, momentum=momentum, batchsize=batchsize)
  end

  # Update function
  function update(s, iter, propcov)
    # Compute temperature
    p = ((nburn - iter)^temper_power + 1) / (nburn^temper_power - 1)
    temper = iter < nburn ? max_temper^p : min_temper

    # Log density of target distribution
    _logprob(s::Vector{Float64}) = logprob(s, temper)

    # Update state
    s .= metropolis(s, _logprob, propcov)

    # Build a message
    msg = "temper: $(round(temper, digits=3))"

    return msg
  end

  # Run MCMC
  msg = ""
  for i in 1:total_iters
    print("\r$(i)/$(total_iters)")
    # Update state
    for t in 1:thin
      msg = update(state, i, propcov.value)
    end

    # Update proposal covariance
    update!(propcov, state)

    # Sanity check
    if verbose > 0 && mod(i, batchsize) == 0
      println(" | proposal covariance: $(round.(propcov.value[1], digits=3)) | $(msg)")
    end

    # Keep samples
    if i > nburn
      out[i - nburn, :] .= state
    end
  end

  println()
  return out, propcov
end

end # module MCMC


text:
###### This file is part of the MendelImpute.jl package.
###### These are wrapper functions for importing VCF/PLINK/BGEN target
###### and reference files.

"""
    convert_gt(b::Bgen, T=Float32)

Imports dosage information and chr/sampleID/pos/snpID/ref/alt into numeric arrays.
Assumes every variant is biallelic (ie only 1 alt allele).

# Input
- `b`: a `Bgen` object
- `T`: Type for genotype array

# Output
- `G`: a `p × n` matrix of type `T`. Each column is a genotype
- `Gchr`: Vector of `String`s holding chromosome number for each variant
- `Gpos`: Vector of `Int` holding each variant's position
- `GsnpID`: Vector of `String`s holding variant ID for each variant
- `Gref`: Vector of `String`s holding reference allele for each variant
- `Galt`: Vector of `String`s holding alterante allele for each variant
"""
function convert_gt(t::Type{T}, b::Bgen) where T <: Real
    n = n_samples(b)
    p = n_variants(b)

    # return arrays
    G = Matrix{Union{Missing, t}}(undef, p, n)
    Gchr = Vector{String}(undef, p)
    Gpos = Vector{Int}(undef, p)
    GsnpID = [String[] for _ in 1:p] # each variant can have >1 rsid, although we don't presently allow this
    Gref = Vector{String}(undef, p)
    Galt = [String[] for _ in 1:p] # each variant can have >1 alt allele, although we don't presently allow this

    # loop over each variant
    i = 1
    for v in iterator(b; from_bgen_start=true)
        dose = ref_allele_dosage!(b, v; T=t) # this reads REF allele as 1
        BGEN.alt_dosage!(dose, v.genotypes.preamble) # switch 2 and 0 (ie treat ALT as 1)
        copyto!(@view(G[i, :]), dose)
        # store chr/pos/snpID/ref/alt info
        Gchr[i], Gpos[i] = chrom(v), pos(v)
        push!(GsnpID[i], rsid(v))
        ref_alt_alleles = alleles(v)
        length(ref_alt_alleles) > 2 && error("Marker $i of BGEN is not biallelic!")
        Gref[i] = ref_alt_alleles[1]
        push!(Galt[i], ref_alt_alleles[2])
        i += 1
        clear!(v)
    end

    # convert NaN to missing
    replace!(G, NaN => missing)

    return G, b.samples, Gchr, Gpos, GsnpID, Gref, Galt
end

"""
    convert_ht(b::Bgen)

Import phased haplotypes as a `BitMatrix`, and store chr/sampleID/pos/snpID/ref/alt.
Assumes every variant is phased and biallelic (ie only 1 alt allele).

# Input
- `b`: a `Bgen` object. Each variant must be phased and samples must be diploid

# Output
- `H`: a `p × 2n` matrix of type `T`. Each column is a haplotype.
- `Hchr`: Vector of `String`s holding chromosome number for each variant
- `Hpos`: Vector of `Int` holding each variant's position
- `HsnpID`: Vector of `String`s holding variant ID for each variant
- `Href`: Vector of `String`s holding reference allele for each variant
- `Halt`: Vector of `String`s holding alterante allele for each variant
"""
function convert_ht(b::Bgen)
    n = 2n_samples(b)
    p = n_variants(b)

    # return arrays
    H = BitMatrix(undef, p, n)
    Hchr = Vector{String}(undef, p)
    Hpos = Vector{Int}(undef, p)
    HsnpID = [String[] for _ in 1:p] # each variant can have >1 rsid, although we don't presently allow this
    Href = Vector{String}(undef, p)
    Halt = [String[] for _ in 1:p] # each variant can have >1 alt allele, although we don't presently allow this

    # loop over each variant
    i = 1
    for v in iterator(b; from_bgen_start=true)
        dose = probabilities!(b, v)
        phased(v) == true || error("variant $(rsid(v)) at position $(pos(v)) not phased!")
        for j in 1:n_samples(b)
            Hi = @view(dose[:, j])
            H[i, 2j - 1] = read_haplotype1(Hi)
            H[i, 2j] = read_haplotype2(Hi)
        end
        # store chr/pos/snpID/ref/alt info
        Hchr[i], Hpos[i] = chrom(v), pos(v)
        push!(HsnpID[i], rsid(v))
        ref_alt_alleles = alleles(v)
        length(ref_alt_alleles) > 2 && error("Marker $i of BGEN is not biallelic!")
        Href[i] = ref_alt_alleles[1]
        push!(Halt[i], ref_alt_alleles[2])
        i += 1
        clear!(v)
    end
    return H, b.samples, Hchr, Hpos, HsnpID, Href, Halt
end

read_haplotype1(Hi::AbstractVector) = Hi[2] ≥ 0.5 ? true : false
read_haplotype2(Hi::AbstractVector) = Hi[4] ≥ 0.5 ? true : false

isplink(tgtfile::AbstractString) = isfile(tgtfile * ".bed") &&
                                   isfile(tgtfile * ".fam") &&
                                   isfile(tgtfile * ".bim")

function import_target(tgtfile::AbstractString, dosage=false)
    if (endswith(tgtfile, ".vcf") || endswith(tgtfile, ".vcf.gz")) && !dosage
        X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt =
            VCFTools.convert_gt(UInt8, tgtfile, trans=true,
            save_snp_info=true, msg = "Importing genotype file...")
    elseif (endswith(tgtfile, ".vcf") || endswith(tgtfile, ".vcf.gz")) && dosage
        X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt =
            VCFTools.convert_ds(Float32, tgtfile, trans=true,
            save_snp_info=true, msg = "Importing genotype file as dosages...")
    elseif endswith(tgtfile, ".bgen")
        samplefile = isfile(tgtfile[1:end-5] * ".sample") ?
            tgtfile[1:end-5] * ".sample" : nothing
        indexfile = isfile(tgtfile * ".bgi") ? tgtfile * ".bgi" : nothing
        bgen = Bgen(tgtfile; sample_path=samplefile, idx_path=indexfile)
        X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt = convert_gt(Float32, bgen)
    elseif isplink(tgtfile)
        dosage && error("PLINK files detected but dosage = true!")
        # convert SnpArray data to matrix.
        X_snpdata = SnpArrays.SnpData(tgtfile)
        X = convert(Matrix{Union{UInt8, Missing}}, Transpose(X_snpdata.snparray))
        _convert_plink_genotypes!(X)
        # get other relevant information
        X_sampleID = X_snpdata.person_info[!, :iid]
        X_chr = X_snpdata.snp_info[!, :chromosome]
        X_pos = X_snpdata.snp_info[!, :position]
        X_ids = X_snpdata.snp_info[!, :snpid]
        X_ref = X_snpdata.snp_info[!, :allele1]
        X_alt = X_snpdata.snp_info[!, :allele2]
    else
        error("Unrecognized target file format: target file can only be VCF" *
            " files (ends in .vcf or .vcf.gz), BGEN (ends in .bgen) or PLINK" *
            " (do not include.bim/bed/fam) and all trio must exist in 1 directory)")
    end
    return X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt
end

"""
    _convert_plink_genotypes!(X::Matrix{Union{UInt8, Missing})

When PLINK binary files are imported into `Matrix{AbstractInt}`, conversion
conventions are not obeyed. So we need to manually set
0x00 => 0
0x01 => missing
0x02 => 1
0x03 => 2
"""
function _convert_plink_genotypes!(X::Matrix{Union{UInt8, Missing}})
    p, n = size(X)
    @inbounds for j in 1:n, i in 1:p
        if X[i, j] == 0x00
            continue
        elseif X[i, j] == 0x01
            X[i, j] = missing
        elseif X[i, j] == 0x02
            X[i, j] = 0x01
        else
            X[i, j] = 0x02
        end
    end
end

function import_compressed_reference(reffile::AbstractString)
    println("Importing reference haplotype data...")
    flush(stdout)

    if endswith(reffile, ".jlso")
        compressed_Hunique = read_jlso(reffile)
    elseif endswith(reffile, ".vcf") || endswith(reffile, ".vcf.gz")
        error("reference panel is VCF format: please first compress reference" *
            " file to .jlso format using the compress_haplotypes() function.")
    elseif endswith(reffile, ".bgen")
        error("reference panel is BGEN format: please first compress reference" *
            " file to .jlso format using the compress_haplotypes() function.")
    elseif endswith(reffile, ".jld2")
        @load reffile compressed_Hunique
    else
        error("Unrecognized reference file format: only VCF (ends in .vcf" *
            " or .vcf.gz), BGEN (ends in .bgen) `.jlso` files are accepted.")
    end
    return compressed_Hunique
end

function import_reference(reffile::AbstractString)
    if endswith(reffile, ".vcf") || endswith(reffile, ".vcf.gz")
        H, H_sampleID, H_chr, H_pos, H_ids, H_ref, H_alt =
            VCFTools.convert_ht(Bool, reffile, trans=true, save_snp_info=true,
            msg="importing reference data...")
    elseif endswith(reffile, ".bgen")
        samplefile = isfile(reffile[1:end-5] * ".sample") ? reffile[1:end-5] * ".sample" : nothing
        indexfile = isfile(reffile * ".bgi") ? reffile * ".bgi" : nothing
        bgen = Bgen(reffile; sample_path=samplefile, idx_path=indexfile)
        H, H_sampleID, H_chr, H_pos, H_ids, H_ref, H_alt = convert_ht(bgen)
    else
        error("Unrecognized reference file format: only VCF (ends in .vcf" *
            " or .vcf.gz) or BGEN (ends in .bgen) files are accepted.")
    end
    return H, H_sampleID, H_chr, H_pos, H_ids, H_ref, H_alt
end

"""
    read_jlso(file::AbstractString)

Imports a `.jlso`-compressed reference haplotype panel.
"""
function read_jlso(reffile::AbstractString)
    loaded = JLSO.load(reffile)
    return loaded[:compressed_Hunique]
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
module ADI

using LinearAlgebra
using ProgressLogging
using Reexport
using Statistics
using StructArrays

@reexport using HCIToolbox

export reconstruct,
       subtract,
       process,
       Classic,
       PCA,
       GreeDS,
       TPCA,
       NMF,
       Framewise,
       SingleSDI,
       DoubleSDI,
       SliceSDI

# core API and default implementations
include("common.jl")
include("design.jl")

# further techniques
include("framewise.jl")
include("sdi.jl")

# algorithms
include("classic.jl")
include("pca.jl")
include("nmf.jl")
include("greeds.jl")

# metrics
include("metrics/Metrics.jl")
@reexport using .Metrics

end


text:
using TuringModels

Turing.setadbackend(:reversediff)

# Turing model
@model m8_3(y) = begin
    α ~ Normal(1, 10)
    σ ~ truncated(Cauchy(0, 1), 0, Inf)

    y .~ Normal(α, σ)
end

y = [-1,1]

# Sample

chns = sample(m8_3(y), NUTS(0.95), 1000);

# Results rethinking

m83rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
alpha 0.09 1.63 -2.13  2.39   959    1
sigma 2.04 2.05  0.68  4.83  1090    1
";

# Describe the posterior samples

describe(chns)

# End of `08/m8.3t.jl`


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

@recipe function f(solution::EstimationSolution{<:RegularGrid}; variables=nothing)
  # grid dimension and size
  sdomain = domain(solution)
  dim = ndims(sdomain)
  sz  = size(sdomain)
  or  = origin(sdomain)
  sp  = spacing(sdomain)

  # valid variables
  validvars = sort(collect(keys(solution.mean)))

  # plot all variables by default
  variables == nothing && (variables = validvars)
  @assert variables ⊆ validvars "invalid variable name"

  # plot layout: mean and variance for each variable
  layout := (length(variables), 2)

  for (i,var) in enumerate(variables)
    # results in grid format
    M = reshape(solution.mean[var], sz)
    V = reshape(solution.variance[var], sz)

    if dim == 1 # plot a line
      x = or[1]:sp[1]:or[1]+(sz[1]-1)*sp[1]
      @series begin
        subplot := 2i - 1
        seriestype --> :path
        legend --> false
        title --> string(var, " mean")
        x, M
      end
      @series begin
        subplot := 2i
        seriestype --> :path
        legend --> false
        title --> string(var, " variance")
        x, V
      end
    elseif dim == 2 # plot a heat map
      @series begin
        subplot := 2i - 1
        seriestype --> :contourf
        seriescolor --> :bluesreds
        title --> string(var, " mean")
        reverse(rotr90(M), dims=2)
      end
      @series begin
        subplot := 2i
        seriestype --> :contourf
        seriescolor --> :bluesreds
        title --> string(var, " variance")
        reverse(rotr90(V), dims=2)
      end
    elseif dim == 3 # plot a volume
      @series begin
        subplot := 2i - 1
        seriestype --> :volume
        seriescolor --> :bluesreds
        title --> string(var, " mean")
        M
      end
      @series begin
        subplot := 2i
        seriestype --> :volume
        seriescolor --> :bluesreds
        title --> string(var, " variance")
        V
      end
    else
      error("cannot plot solution in more than 3 dimensions")
    end
  end
end


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

# make sure that typeinf is executed before turning on typeinf_ext
# this ensures that typeinf_ext doesn't recurse before it can add the item to the workq
# especially try to make sure any recursive and leaf functions have concrete signatures,
# since we won't be able to specialize & infer them at runtime

let fs = Any[typeinf_ext, typeinf, typeinf_edge, pure_eval_call, optimize, run_passes],
    world = ccall(:jl_get_world_counter, UInt, ())
    for x in T_FFUNC_VAL
        push!(fs, x[3])
    end
    for i = 1:length(T_IFUNC)
        if isassigned(T_IFUNC, i)
            x = T_IFUNC[i]
            push!(fs, x[3])
        else
            println(stderr, "WARNING: tfunc missing for ", reinterpret(IntrinsicFunction, Int32(i)))
        end
    end
    for f in fs
        for m in _methods_by_ftype(Tuple{typeof(f), Vararg{Any}}, 10, typemax(UInt))
            # remove any TypeVars from the intersection
            typ = Any[m[1].parameters...]
            for i = 1:length(typ)
                if isa(typ[i], TypeVar)
                    typ[i] = typ[i].ub
                end
            end
            typeinf_type(m[3], Tuple{typ...}, m[2], true, Params(world))
        end
    end
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork