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 Color{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{T}.
# - Odest and Osrc are 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 motivation for this arrangement is that Julia doesn't (yet) support
# "triangular dispatch", e.g.,
#     convert{T,C}(::Type{C{T}}, c)
# The various arguments of _convert therefore "peel off" element types
# (or guarantee them) so that comparisons may be made via
# dispatch. The alternative design is
#    for C in parametric_colors
#       @eval convert{T}(::Type{$C{T}}, c) = ...
#       @eval convert(   ::Type{$C},    c) = convert($C{eltype(c)}, c)
#       ...
#    end
# but this requires a fair amount of codegen (especially for all
# the various alpha variants) and can break if not all C support
# the same eltypes.
#
# Note that ColorTypes handles the cases where Odest == Osrc, or they
# are both subtypes of AbstractRGB.  Therefore, here we only have to
# deal with conversions between different spaces.

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

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

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

# Conversions from grayscale
# --------------------------
cnvt{C<:Color3}(::Type{C}, g::AbstractGray)  = cnvt(C, convert(RGB{eltype(C)}, g))

macro mul3x3(T, M, c1, c2, c3)
    esc(quote
        @inbounds ret = $T($M[1,1]*$c1 + $M[1,2]*$c2 + $M[1,3]*$c3,
                           $M[2,1]*$c1 + $M[2,2]*$c2 + $M[2,3]*$c3,
                           $M[3,1]*$c1 + $M[3,2]*$c2 + $M[3,3]*$c3)
        ret
        end)
end

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

correct_gamut{CV<:AbstractRGB}(c::CV) = CV(clamp01(red(c)), clamp01(green(c)), clamp01(blue(c)))
clamp01{T<:Fractional}(v::T) = ifelse(v < zero(T), zero(T), ifelse(v > one(T), one(T), v))

function srgb_compand(v)
    v <= 0.0031308 ? 12.92v : 1.055v^(1/2.4) - 0.055
end

cnvt{CV<:AbstractRGB}(::Type{CV}, c::AbstractRGB) = CV(red(c), green(c), blue(c))

function cnvt{CV<:AbstractRGB}(::Type{CV}, c::HSV)
    h = c.h / 60
    i = floor(Int, h)
    f = h - i
    if i & 1 == 0
        f = 1 - f
    end
    m = c.v * (1 - c.s)
    n = c.v * (1 - c.s * f)
    if i == 6 || i == 0; CV(c.v, n, m)
    elseif i == 1;       CV(n, c.v, m)
    elseif i == 2;       CV(m, c.v, n)
    elseif i == 3;       CV(m, n, c.v)
    elseif i == 4;       CV(n, m, c.v)
    else;                CV(c.v, m, n)
    end
end

function qtrans(u, v, hue)
    if     hue > 360; hue -= 360
    elseif hue < 0;   hue += 360
    end

    if     hue < 60;  u + (v - u) * hue / 60
    elseif hue < 180; v
    elseif hue < 240; u + (v - u) * (240 - hue) / 60
    else;             u
    end
end

function cnvt{CV<:AbstractRGB}(::Type{CV}, c::HSL)
    v = c.l <= 0.5 ? c.l * (1 + c.s) : c.l + c.s - (c.l * c.s)
    u = 2 * c.l - v

    if c.s == 0; CV(c.l, c.l, c.l)
    else;        CV(qtrans(u, v, c.h + 120),
                    qtrans(u, v, c.h),
                    qtrans(u, v, c.h - 120))
    end
end

function cnvt{CV<:AbstractRGB}(::Type{CV}, c::HSI)
    h, s, i = c.h, c.s, c.i
    while(h > 360) h -= 360 end
    while(h < 0) h += 360 end
    is = i*s
    if h < 120
        cosr = cosd(h) / cosd(60-h)
        CV(i+is*cosr, i+is*(1-cosr), i-is)
    elseif h < 240
        cosr = cosd(h-120) / cosd(180-h)
        CV(i-is, i+is*cosr, i+is*(1-cosr))
    else
        cosr = cosd(h-240) / cosd(300-h)
        CV(i+is*(1-cosr), i-is, i+is*cosr)
    end
end

function cnvt{CV<:AbstractRGB}(::Type{CV}, c::XYZ)
    r =  3.2404542*c.x - 1.5371385*c.y - 0.4985314*c.z
    g = -0.9692660*c.x + 1.8760108*c.y + 0.0415560*c.z
    b =  0.0556434*c.x - 0.2040259*c.y + 1.0572252*c.z
    CV(clamp01(srgb_compand(r)),
       clamp01(srgb_compand(g)),
       clamp01(srgb_compand(b)))
end

function cnvt{CV<:AbstractRGB}(::Type{CV}, c::YIQ)
    cc = correct_gamut(c)
    CV(clamp01(cc.y+0.9563*cc.i+0.6210*cc.q),
       clamp01(cc.y-0.2721*cc.i-0.6474*cc.q),
       clamp01(cc.y-1.1070*cc.i+1.7046*cc.q))
end

function cnvt{CV<:AbstractRGB}(::Type{CV}, c::YCbCr)
    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

cnvt{CV<:AbstractRGB}(::Type{CV}, c::LCHab)  = cnvt(CV, convert(Lab{eltype(c)}, c))
cnvt{CV<:AbstractRGB}(::Type{CV}, c::LCHuv)  = cnvt(CV, convert(Luv{eltype(c)}, c))
cnvt{CV<:AbstractRGB}(::Type{CV}, c::Color3)    = cnvt(CV, convert(XYZ{eltype(c)}, c))

cnvt{CV<:AbstractRGB{UFixed8}}(::Type{CV}, c::RGB24) = CV(UFixed8(c.color&0x00ff0000>>>16,0), UFixed8(c.color&0x0000ff00>>>8,0), UFixed8(c.color&0x000000ff,0))
cnvt{CV<:AbstractRGB}(::Type{CV}, c::RGB24) = CV((c.color&0x00ff0000>>>16)/255, ((c.color&0x0000ff00)>>>8)/255, (c.color&0x000000ff)/255)

function cnvt{CV<:AbstractRGB}(::Type{CV}, c::AbstractGray)
    g = convert(eltype(CV), gray(c))
    CV(g, g, g)
end

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

function cnvt{T}(::Type{HSV{T}}, c::AbstractRGB)
    c_min = Float64(min(red(c), green(c), blue(c)))
    c_max = Float64(max(red(c), green(c), blue(c)))
    if c_min == c_max
        return HSV{T}(zero(T), zero(T), c_max)
    end

    if c_min == red(c)
        f = Float64(green(c)) - Float64(blue(c))
        i = 3
    elseif c_min == green(c)
        f = Float64(blue(c)) - Float64(red(c))
        i = 5
    else
        f = Float64(red(c)) - Float64(green(c))
        i = 1
    end

    HSV{T}(60 * (i - f / (c_max - c_min)),
        (c_max - c_min) / c_max,
        c_max)
end

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

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

function cnvt{T}(::Type{HSL{T}}, c::AbstractRGB)
    r, g, b = T(red(c)), T(green(c)), T(blue(c))
    c_min = min(r, g, b)
    c_max = max(r, g, b)
    l = (c_max + c_min) / 2

    if c_max == c_min
        return HSL(zero(T), zero(T), l)
    end

    if l < 0.5; s = (c_max - c_min) / (c_max + c_min)
    else;       s = (c_max - c_min) / (convert(T, 2) - c_max - c_min)
    end

    if c_max == red(c)
        h = (g - b) / (c_max - c_min)
    elseif c_max == green(c)
        h = convert(T, 2) + (b - r) / (c_max - c_min)
    else
        h = convert(T, 4) + (r - g) / (c_max - c_min)
    end

    h *= 60
    if h < 0
        h += 360
    elseif h > 360
        h -= 360
    end

    HSL{T}(h,s,l)
end

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

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

function cnvt{T}(::Type{HSI{T}}, c::AbstractRGB)
    rgb = correct_gamut(c)
    r, g, b = red(rgb), green(rgb), blue(rgb)
    isum = r+g+b
    i = isum/3
    m = min(r, g, b)
    s = i > 0 ? 1-m/i : 1-one(i)/one(i) # the latter is a type-stable 0
    val = (r-(g+b)/2)/sqrt(r^2+g^2+b^2-r*g-r*b-g*b)
    val = clamp(val, -one(val), one(val))
    h = acosd(val)
    if b > g
        h = 360-h
    end
    HSI{T}(h, s, i)
end

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

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

function invert_rgb_compand(v)
    v <= 0.04045 ? v/12.92 : ((v+0.055) /1.055)^2.4
end

function cnvt{T}(::Type{XYZ{T}}, c::AbstractRGB)
    r, g, b = invert_rgb_compand(red(c)), invert_rgb_compand(green(c)), invert_rgb_compand(blue(c))
    XYZ{T}(0.4124564*r + 0.3575761*g + 0.1804375*b,
           0.2126729*r + 0.7151522*g + 0.0721750*b,
           0.0193339*r + 0.1191920*g + 0.9503041*b)
end

cnvt{T}(::Type{XYZ{T}}, c::HSV) = cnvt(XYZ{T}, convert(RGB{T}, c))
cnvt{T}(::Type{XYZ{T}}, c::HSL) = cnvt(XYZ{T}, convert(RGB{T}, c))
cnvt{T}(::Type{XYZ{T}}, c::HSI) = cnvt(XYZ{T}, convert(RGB{T}, c))

function cnvt{T}(::Type{XYZ{T}}, c::xyY)
    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
const xyz_kappa   = 24389 / 27

convert(::Type{XYZ}, c, wp::XYZ) = convert(XYZ{eltype(wp)}, c, wp)
convert{T}(::Type{XYZ{T}}, c, wp::XYZ) = cnvt(XYZ{T}, c, wp)
function cnvt{T}(::Type{XYZ{T}}, c::Lab, wp::XYZ)
    fy = (c.l + 16) / 116
    fx = c.a / 500 + fy
    fz = fy - c.b / 200

    fx3 = fx^3
    fz3 = fz^3

    x = fx3 > xyz_epsilon ? fx3 : (116fx - 16) / xyz_kappa
    y = c.l > xyz_kappa * xyz_epsilon ? ((c. l+ 16) / 116)^3 : c.l / xyz_kappa
    z = fz3 > xyz_epsilon ? fz3 : (116fz - 16) / xyz_kappa

    XYZ{T}(x*wp.x, y*wp.y, z*wp.z)
end

cnvt{T}(::Type{XYZ{T}}, c::Lab)   = convert(XYZ{T}, c, WP_DEFAULT)
cnvt{T}(::Type{XYZ{T}}, c::LCHab) = cnvt(XYZ{T}, convert(Lab{T}, c))
cnvt{T}(::Type{XYZ{T}}, c::DIN99) = cnvt(XYZ{T}, convert(Lab{T}, c))
cnvt{T}(::Type{XYZ{T}}, c::DIN99o) = cnvt(XYZ{T}, convert(Lab{T}, c))

cnvt{T<:UFixed}(::Type{XYZ{T}}, c::LCHab) = cnvt(XYZ{T}, convert(Lab{eltype(c)}, c))
cnvt{T<:UFixed}(::Type{XYZ{T}}, c::DIN99) = cnvt(XYZ{T}, convert(Lab{eltype(c)}, c))
cnvt{T<:UFixed}(::Type{XYZ{T}}, c::DIN99o) = cnvt(XYZ{T}, convert(Lab{eltype(c)}, c))

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

function cnvt{T}(::Type{XYZ{T}}, c::Luv, wp::XYZ = WP_DEFAULT)
    (u_wp, v_wp) = xyz_to_uv(wp)

    a = (52 * (c.l==0 ? zero(T) : c.l / (c.u + 13 * c.l * u_wp)) - 1) / 3
    y = c.l > xyz_kappa * xyz_epsilon ? wp.y * ((c.l + 16) / 116)^3 : wp.y * c.l / xyz_kappa
    b = -5y
    d = y * (39 * (c.l==0 ? zero(T) : c.l / (c.v + 13 * c.l * v_wp)) - 5)
    x = d==b ? zero(T) : (d - b) / (a + 1/3)
    z = a * x + b + zero(T)

    XYZ{T}(x, y, z)
end

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

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

    # Go back to C-h space
    # FIXME: Clean this up (why is there no atan2d?)
    h = rad2deg(atan2(c.b,c.a)) + 50
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end

    C = sqrt(c.a^2 + c.b^2)

    # Intermediate terms
    G = (exp(C/22.5)-1)/0.06
    f = G*sind(h - 50)
    ee = G*cosd(h - 50)

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

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

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

end

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

cnvt{T}(::Type{XYZ{T}}, c::YIQ) = cnvt(XYZ{T}, convert(RGB{T}, c))
cnvt{T}(::Type{XYZ{T}}, c::YCbCr) = cnvt(XYZ{T}, convert(RGB{T}, c))

cnvt{T}(::Type{XYZ{T}}, c::RGB24) = cnvt(XYZ{T}, convert(RGB{T}, c))

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

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

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

    xyY{T}(x, y, convert(typeof(x), c.y))

end

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

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

cnvt{T}(::Type{Lab{T}}, c::AbstractRGB) = convert(Lab{T}, convert(XYZ{T}, c))
cnvt{T}(::Type{Lab{T}}, c::HSV) = cnvt(Lab{T}, convert(RGB{T}, c))
cnvt{T}(::Type{Lab{T}}, c::HSL) = cnvt(Lab{T}, convert(RGB{T}, c))

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

function fxyz2lab(v)
    v > xyz_epsilon ? cbrt(v) : (xyz_kappa * v + 16) / 116
end
function cnvt{T}(::Type{Lab{T}}, c::XYZ, wp::XYZ)
    fx, fy, fz = fxyz2lab(c.x / wp.x), fxyz2lab(c.y / wp.y), fxyz2lab(c.z / wp.z)
    Lab{T}(116fy - 16, 500(fx - fy), 200(fy - fz))
end

cnvt{T}(::Type{Lab{T}}, c::XYZ{T}) = cnvt(Lab{T}, c, WP_DEFAULT)

function cnvt{T}(::Type{Lab{T}}, c::LCHab)
    hr = deg2rad(c.h)
    Lab{T}(c.l, c.c * cos(hr), c.c * sin(hr))
end

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

    # 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 = sqrt(c.a^2 + c.b^2)

    # NOTE: This is calculated in degrees, against the standard, to save
    # computation steps later.
    if (c.a > 0 && c.b >= 0)
        h = atand(c.b/c.a)
    elseif (c.a == 0 && c.b > 0)
        h = 90
    elseif (c.a < 0)
        h = 180+atand(c.b/c.a)
    elseif (c.a == 0 && c.b < 0)
        h = 270
    elseif (c.a > 0 && c.b <= 0)
        h = 360 + atand(c.b/c.a)
    else
        h = 0
    end

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

    # Temporary redness
    ee = g*cosd(h)

    # Temporary yellowness
    f = g*sind(h)

    # 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 = (e^(c.l*ke/105.51)-1)/0.0158

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

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

    # 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 = sqrt(c.a^2 + c.b^2)

    # hue angle h99o
    h = atan2(c.b, c.a)

    # revert rotation by 26°
    ho= rad2deg(h)-26

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

    # Temporary redness
    eo = g*cosd(ho)

    # Temporary yellowness
    fo = g*sind(ho)

    # 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 = (e^(c.l*ke/303.67)-1)/0.0039

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

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

# Everything to Luv
# -----------------

cnvt{T}(::Type{Luv{T}}, c::AbstractRGB) = cnvt(Luv{T}, convert(XYZ{T}, c))
cnvt{T}(::Type{Luv{T}}, c::HSV) = cnvt(Luv{T}, convert(RGB{T}, c))
cnvt{T}(::Type{Luv{T}}, c::HSL) = cnvt(Luv{T}, convert(RGB{T}, c))

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

function cnvt{T}(::Type{Luv{T}}, c::XYZ, wp::XYZ = WP_DEFAULT)
    (u_wp, v_wp) = xyz_to_uv(wp)
    (u_, v_) = xyz_to_uv(c)

    y = c.y / wp.y

    l = y > xyz_epsilon ? 116 * cbrt(y) - 16 : xyz_kappa * y
    u = 13 * l * (u_ - u_wp) + zero(T)
    v = 13 * l * (v_ - v_wp) + zero(T)

    Luv{T}(l, u, v)
end

function cnvt{T}(::Type{Luv{T}}, c::LCHuv)
    hr = deg2rad(c.h)
    Luv{T}(c.l, c.c * cos(hr), c.c * sin(hr))
end

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

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

function cnvt{T}(::Type{LCHuv{T}}, c::Luv)
    h = rad2deg(atan2(c.v, c.u))
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end
    LCHuv{T}(c.l, sqrt(c.u^2 + c.v^2), h)
end

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

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

function cnvt{T}(::Type{LCHab{T}}, c::Lab)
    h = rad2deg(atan2(c.b, c.a))
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end
    LCHab{T}(c.l, sqrt(c.a^2 + c.b^2), h)
end

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

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

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

    # 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)

    # Hue angle
    # Calculated in degrees, against the specification.
    if (ee > 0 && f >= 0)
        h = atand(f/ee)
    elseif (ee == 0 && f > 0)
        h = 90
    elseif (ee < 0)
        h = 180+atand(f/ee)
    elseif (ee == 0 && f < 0)
        h = 270
    elseif (ee > 0 && f <= 0)
        h = 360 + atand(f/ee)
    else
        h = 0
    end

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

    # DIN99 chromaticities
    a99, b99 = cc*cosd(h), cc*sind(h)

    DIN99{T}(l99, a99, b99)

end

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

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

function cnvt{T}(::Type{DIN99d{T}}, c::XYZ{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.22*log(1+0.0036*lab.l)

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

    # Calculate hue/chroma
    C = 22.5*log(1+0.06*G)
    # FIXME: Clean this up (why is there no atan2d?)
    h = rad2deg(atan2(f,ee)) + 50
    while h > 360; h -= 360; end
    while h < 0;   h += 360; end

    DIN99d{T}(adj_L, C*cosd(h), C*sind(h))

end

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

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

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

    # 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)
    ho = atan2(fo,eo)
    # rotation of the color space by 26°
    h  = rad2deg(ho) + 26

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

    # DIN99o chromaticities
    a99, b99 = cc*cosd(h), cc*sind(h)

    DIN99o{T}(l99, a99, b99)

end

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

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

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

const CAT97s_INV = inv(CAT97s)

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

const CAT02_INV = inv(CAT02)

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

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

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

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

function cnvt{T}(::Type{YIQ{T}}, c::AbstractRGB)
    rgb = correct_gamut(c)
    YIQ{T}(0.299*red(rgb)+0.587*green(rgb)+0.114*blue(rgb),
           0.595716*red(rgb)-0.274453*green(rgb)-0.321263*blue(rgb),
           0.211456*red(rgb)-0.522591*green(rgb)+0.311135*blue(rgb))
end

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

# Everything to YCbCr
# -------------------

correct_gamut{T}(c::YCbCr{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)))

function cnvt{T}(::Type{YCbCr{T}}, c::AbstractRGB)
    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{T}(::Type{YCbCr{T}}, c::Color3) = cnvt(YCbCr{T}, convert(RGB{T}, c))

# Everything to RGB24
# -------------------

convert(::Type{RGB24}, c::RGB24) = c
convert(::Type{RGB24}, c::AbstractRGB{UFixed8}) = RGB24(red(c), green(c), blue(c))
convert(::Type{RGB24}, c::AbstractRGB) = RGB24((round(UInt8, 255*red(c))%UInt32)<<16 +
                                               (round(UInt8, 255*green(c))%UInt32)<<8 +
                                                round(UInt8, 255*blue(c))%UInt32)

convert(::Type{RGB24}, c::Color) = convert(RGB24, convert(RGB{UFixed8}, c))

# To ARGB32
# ----------------

convert(::Type{ARGB32}, c::ARGB32) = c
convert{CV<:AbstractRGB{UFixed8}}(::Type{ARGB32}, c::TransparentColor{CV}) =
    ARGB32(red(c), green(c), blue(c), alpha(c))
convert(::Type{ARGB32}, c::TransparentColor) =
    ARGB32(convert(RGB24, c).color | (round(UInt8, 255*alpha(c))%UInt32)<<24)
convert(::Type{ARGB32}, c::Color) = ARGB32(convert(RGB24, c).color | 0xff000000)
convert(::Type{ARGB32}, c::Color, alpha) = ARGB32(convert(RGB24, c).color | round(UInt32, 255*alpha)<<24)

# To Gray
# -------

# Rec 601 luma conversion
if VERSION < v"0.4.0-dev+1827"
    const unsafe_trunc = itrunc
else
    const unsafe_trunc = Base.unsafe_trunc
end

convert{T<:UFixed}(::Type{Gray{T}}, x::AbstractRGB{T}) = (TU = FixedPointNumbers.rawtype(T); Gray{T}(T(round(TU, min(typemax(TU), 0.299f0*reinterpret(x.r) + 0.587f0*reinterpret(x.g) + 0.114f0*reinterpret(x.b))), 0)))
convert{T}(::Type{Gray{T}}, x::AbstractRGB) = convert(Gray{T}, 0.299f0*x.r + 0.587f0*x.g + 0.114f0*x.b)


text:
using Test
using Distributed, Base.Threads
using Base.Iterators: product

exeflags = ("--startup-file=no",
            "--check-bounds=yes",
            "--depwarn=error",
            "--threads=2")

function call_on(f, wid, tid)
  remotecall(wid) do
    t = Task(f)
    ccall(:jl_set_task_tid, Cvoid, (Any, Cint), t, tid-1)
    schedule(t)
    @assert threadid(t) == tid
    t
  end
end

# Run function on process holding the data to only serialize the result of f.
# This becomes useful for things that cannot be serialized (e.g. running tasks)
# or that would be unnecessarily big if serialized.
fetch_from_owner(f, rr) = remotecall_fetch(f∘fetch, rr.where, rr)

isdone(rr) = fetch_from_owner(istaskdone, rr)
isfailed(rr) = fetch_from_owner(istaskfailed, rr)

@testset "RemoteChannel allows put!/take! from thread other than 1" begin
  ws = ts = product(1:2, 1:2)
  timeout = 10.0
  @testset "from worker $w1 to $w2 via 1" for (w1, w2) in ws
    @testset "from thread $w1.$t1 to $w2.$t2" for (t1, t2) in ts
      # We want (the default) lazyness, so that we wait for `Worker.c_state`!
      procs_added = addprocs(2; exeflags, lazy=true)
      @everywhere procs_added using Base.Threads
      p1 = procs_added[w1]
      p2 = procs_added[w2]
      chan_id = first(procs_added)
      chan = RemoteChannel(chan_id)
      send = call_on(p1, t1) do
        put!(chan, nothing)
      end
      recv = call_on(p2, t2) do
        take!(chan)
      end
      timedwait(() -> isdone(send) && isdone(recv), timeout)
      @test isdone(send)
      @test isdone(recv)
      @test !isfailed(send)
      @test !isfailed(recv)
      rmprocs(procs_added)
    end
  end
end

ext: jl

max_stars_repo_name: sudoWin/julia

max_issues_repo_name: sudoWin/julia

max_forks_repo_name: sudoWin/julia


text:
SUITE["arithmetic"] = BenchmarkGroup()

SUITE["arithmetic"]["DatePeriod"] = begin
    zdt = ZonedDateTime(2000, 1, 2, tz"America/Winnipeg")
    @benchmarkable $zdt + $(Day(1))
end
SUITE["arithmetic"]["TimePeriod"] = begin
    zdt = ZonedDateTime(2000, 1, 2, tz"America/Winnipeg")
    @benchmarkable $zdt + $(Hour(1))
end

SUITE["arithmetic"]["broadcast"] = BenchmarkGroup()

SUITE["arithmetic"]["broadcast"]["VariableTimeZone/TimePeriod"] = begin
    @benchmarkable $(gen_range(Hour(100), tz"America/Winnipeg")) .+ Hour(1)
end
SUITE["arithmetic"]["broadcast"]["VariableTimeZone/DatePeriod"] = begin
    @benchmarkable $(gen_range(Day(100), tz"America/Winnipeg")) .+ Day(1)
end
SUITE["arithmetic"]["broadcast"]["FixedTimeZone/TimePeriod"] = begin
    @benchmarkable $(gen_range(Hour(100), tz"UTC")) .+ Hour(1)
end
SUITE["arithmetic"]["broadcast"]["FixedTimeZone/DatePeriod"] = begin
    @benchmarkable $(gen_range(Day(100), tz"UTC")) .+ Day(1)
end


text:
#some useful functions
glorot_uniform(dims...) = (rand(Float64, dims...) .- 0.5) .* sqrt(24.0/sum(dims))

relu(x::Real) = max(zero(x), x)

sigmoid(x::Real) = one(x) / (one(x) + exp(-x))

softplus(x::Real) = ifelse(x > 0, x + log1p(exp(-x)), log1p(exp(x)))


text:
mutable struct SolveState{T,matT<:FieldMatrix, vecT<:FieldVector}
    A::matT
    iD::vecT
    x::vecT
    r::vecT
    ϵ::vecT
    child::Union{SolveState, Nothing}
    P::Union{FieldMatrix, Nothing}
    function SolveState(A::FieldMatrix,x::FieldVector{T},r::FieldVector,invtol=1e-8) where T
        iD = zero(x)
        @loop iD[I] = abs(A.D[I])>invtol ? inv(A.D[I]) : zero(T)
        new{T,typeof(A),typeof(x)}(A,iD,x,r,zero(x),nothing,nothing)
    end
end
Base.show(io::IO, ::MIME"text/plain", st::SolveState{T}) where T = print(io, "SolveState{",T,"}:\n   ", st)
Base.show(io::IO, st::SolveState) = print(io, "residual=",norm(st.r),"\n   ", st.child)

@fastmath resid!(r,A,x) = (@loop r[I] = r[I]-mult(I,A.L,A.D,x);r)
function residual(A,x,b)
    r=zero(x)
    @loop r[I] = b[I]
    resid!(r,A,x)
end
@fastmath function increment!(st;resid=true,kw...)
    @loop st.x[I] = st.x[I]+st.ϵ[I]
    resid && resid!(st.r,st.A,st.ϵ)
    st
end

epsr(st) = eps(real(eltype(st.r)))
function iterate!(st::SolveState,iterator!::Function;
                  abstol::Real = 20*epsr(st),reltol::Real = √epsr(st),
                  mxiter::Int = size(st.A,2), log::Bool = false, kw...)
    res0,i = norm(st.r),0
    log && (hist = Vector{eltype(st.r)}(undef,mxiter+1); hist[1] = res0)
    res = res0
    while res>max(abstol,reltol*res0) && i<mxiter
        iterator!(st;kw...)
        res,i = norm(st.r),i+1
        log && (hist[i+1] = res)
    end
    return log ? resize!(hist,i+1) : i
end

gs!(st::SolveState;kw...) = iterate!(st,GS!;kw...)
function gs(A::AbstractMatrix,b::AbstractVector;inner=8,kw...)
    x = zero(b)
    return x,gs!(SolveState(A,x,copy(b));inner,kw...)
end

@fastmath function GS!(st;inner=2,kw...)
    @loop st.ϵ[I] = st.r[I]*st.iD[I]
    @inbounds for _ ∈ 1:inner, I ∈ st.iD.R
        st.ϵ[I] = st.iD[I]*(st.r[I]-multL(I,st.A.L,st.ϵ)-multU(I,st.A.L,st.ϵ))
    end
    increment!(st;kw...)
end


text:
# Define some example data
import KernelDensity.UnivariateKDE
o1 = UncertainValue(Normal, 0, 0.5)
o2 = UncertainValue(Normal, 2.0, 0.1)
o3 = UncertainValue(Uniform, 0, 4)
o4 = UncertainValue(Uniform, rand(100))
o5 = UncertainValue(Beta, 4, 5)
o6 = UncertainValue(Gamma, 4, 5)
o7 = UncertainValue(Frechet, 1, 2)
o8 = UncertainValue(BetaPrime, 1, 2)
o9 = UncertainValue(BetaBinomial, 10, 3, 2)
o10 = UncertainValue(Binomial, 10, 0.3)
o11 = UncertainValue(rand(100), rand(100))
o12 = UncertainValue(2)
o13 = UncertainValue(2.3)
o14 = UncertainValue([2, 3, 4], [0.3, 0.4, 0.3])
o15 = UncertainValue([rand() for i = 1:100])

uvals1 = [o1, o2, o3, o4, o5, o6, o7, o8, o9, o11, o12, o13, o14, o1, o2, o3, o15, o4, o5,
    o6, o7, o8, o9, o11, o12, o13, o14, o1, o2, o3, o4, o5, o6, o7, o8, o9, o11, o12, o13,
    o14, o1, o2, o3, o4, o5, o6, o7, o8, o9, o11, o12, o13, o14]

uvals2 = [uvals1[i] for i in [rand(1:length(uvals1)) for i = 1:length(uvals1)]]

idxs = UncertainIndexDataset([UncertainValue(Normal, i, 0.8) for i = 1:length(uvals1)])

UVD = UncertainValueDataset(uvals1)
UD = UncertainDataset(uvals2)
UIDX = UncertainIndexDataset(idxs)
UV = UncertainValueDataset(uvals1)

# Resampling interface should work for all subtypes of AbstractUncertainValueDataset,
# both with and without sampling constraints.
# We need to check UncertainDataset, UncertainIndexDataset, and UncertainValueDataset.

n = 3

#resample_elwise(uvd::AbstractUncertainValueDataset, n::Int)
@test length(resample_elwise(UIDX, n)) == length(UIDX)
@test length(resample_elwise(UIDX, n)[1]) == n
@test length(resample_elwise(UD, n)) == length(UD)
@test length(resample_elwise(UD, n)[1]) == n
@test length(resample_elwise(UVD, n)) == length(UVD)
@test length(resample_elwise(UVD, n)[1]) == n

@test length(resample_elwise(UVD)[1]) == 1
@test length(resample_elwise(UVD)) == length(UVD)

n = 5
@test length(resample_elwise(UD, [TruncateQuantiles(0.1, 0.9) for i = 1:length(UD)], n)) == length(UD)
@test length(resample_elwise(UD, [TruncateQuantiles(0.1, 0.9) for i = 1:length(UD)], n)[1]) == n
@test length(resample_elwise(UVD, [TruncateQuantiles(0.1, 0.9) for i = 1:length(UVD)], n)) == length(UD)
@test length(resample_elwise(UVD, [TruncateQuantiles(0.1, 0.9) for i = 1:length(UVD)], n)[1]) == n
@test length(resample_elwise(UIDX, [TruncateQuantiles(0.1, 0.9) for i = 1:length(UIDX)], n)) == length(UD)
@test length(resample_elwise(UIDX, [TruncateQuantiles(0.1, 0.9) for i = 1:length(UIDX)], n)[1]) == n

# resample_elwise(uvd::AbstractUncertainValueDataset, constraint::SamplingConstraint, n::Int)
@test length(resample_elwise(UIDX, TruncateQuantiles(0.1, 0.9), n)) == length(UIDX)
@test length(resample_elwise(UIDX, TruncateQuantiles(0.1, 0.9), n)[1]) == n
@test length(resample_elwise(UD, TruncateQuantiles(0.1, 0.9), n)) == length(UD)
@test length(resample_elwise(UD, TruncateQuantiles(0.1, 0.9), n)[1]) == n
@test length(resample_elwise(UVD, TruncateQuantiles(0.1, 0.9), n)) == length(UVD)
@test length(resample_elwise(UVD, TruncateQuantiles(0.1, 0.9), n)[1]) == n


text:
##==============================================================================
## AbstractDataEntry - compare
##==============================================================================

import Base: ==

@generated function ==(x::T, y::T) where T <: AbstractDataEntry
    mapreduce(n -> :(x.$n == y.$n), (a,b)->:($a && $b), fieldnames(x))
end

##==============================================================================
## AbstractDataEntry - common
##==============================================================================

"""
    $(SIGNATURES)
Function to generate source string - userId|robotId|sessionId|varLabel
"""
buildSourceString(dfg::AbstractDFG, label::Symbol) =
    "$(dfg.userId)|$(dfg.robotId)|$(dfg.sessionId)|$label"

##==============================================================================
## AbstractDataEntry - CRUD
##==============================================================================

"""
    $(SIGNATURES)
Get data entry
"""
function getDataEntry(var::AbstractDFGVariable, key::Symbol)
    !hasDataEntry(var, key) && error("Data entry $(key) does not exist in variable $(getLabel(var))")
    return var.dataDict[key]
end

function getDataEntry(dfg::AbstractDFG, label::Symbol, key::Symbol)
    return getDataEntry(getVariable(dfg, label), key)
end

"""
    $(SIGNATURES)
Add Data Entry to a DFG variable
"""
function addDataEntry!(var::AbstractDFGVariable, bde::AbstractDataEntry)
    haskey(var.dataDict, bde.label) && error("Data entry $(bde.label) already exists in variable $(getLabel(var))")
    var.dataDict[bde.label] = bde
    return bde
end

"""
    $(SIGNATURES)
Add Data Entry to distributed factor graph.
Should be extended if DFG variable is not returned by reference.
"""
function addDataEntry!(dfg::AbstractDFG, label::Symbol, bde::AbstractDataEntry)
    return addDataEntry!(getVariable(dfg, label), bde)
end

"""
    $(SIGNATURES)
Update data entry

DevNote
- DF, unclear if `update` verb is applicable in this case, see #404
"""
function updateDataEntry!(var::AbstractDFGVariable,  bde::AbstractDataEntry)
    !haskey(var.dataDict, bde.label) && (@warn "$(bde.label) does not exist in variable $(getLabel(var)), adding")
    var.dataDict[bde.label] = bde
    return bde
end
function updateDataEntry!(dfg::AbstractDFG, label::Symbol,  bde::AbstractDataEntry)
    # !isVariable(dfg, label) && return nothing
    return updateDataEntry!(getVariable(dfg, label), bde)
end

"""
    $(SIGNATURES)
Delete data entry from the factor graph.
Note this doesn't remove it from any data stores.

Notes:
- users responsibility to delete data in db before deleting entry
"""
function deleteDataEntry!(var::AbstractDFGVariable, key::Symbol)
    return pop!(var.dataDict, key)
end
function deleteDataEntry!(dfg::AbstractDFG, label::Symbol, key::Symbol)
    #users responsibility to delete data in db before deleting entry
    # !isVariable(dfg, label) && return nothing
    return deleteDataEntry!(getVariable(dfg, label), key)
end

function deleteDataEntry!(var::AbstractDFGVariable, entry::AbstractDataEntry)
    #users responsibility to delete data in db before deleting entry
    return deleteDataEntry!(var, entry.label)
end

##==============================================================================
## AbstractDataEntry - Helper functions, Lists, etc
##==============================================================================

"""
    $SIGNATURES

Does a data entry (element) exist at `key`.
"""
hasDataEntry(var::DFGVariable, key::Symbol) = haskey(var.dataDict, key)

"""
    $(SIGNATURES)
Get data entries, Vector{AbstractDataEntry}
"""
function getDataEntries(var::AbstractDFGVariable)
    #or should we return the iterator, Base.ValueIterator{Dict{Symbol,AbstractDataEntry}}?
    collect(values(var.dataDict))
end
function getDataEntries(dfg::AbstractDFG, label::Symbol)
    # !isVariable(dfg, label) && return nothing
    #or should we return the iterator, Base.ValueIterator{Dict{Symbol,AbstractDataEntry}}?
    getDataEntries(getVariable(dfg, label))
end

"""
    $(SIGNATURES)
listDataEntries
"""
function listDataEntries(var::AbstractDFGVariable)
    collect(keys(var.dataDict))
end

function listDataEntries(dfg::AbstractDFG, label::Symbol)
    # !isVariable(dfg, label) && return nothing
    listDataEntries(getVariable(dfg, label))
end


text:
#=
# 125: Terminal flux calculation via test functions
([source code](SOURCE_URL))

For a rather comprehensive explanation
see [225: Terminal flux calculation via test functions, nD](@ref)
=#

module Example125_TestFunctions1D
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

function main(;n=100,Plotter=nothing,verbose=false,unknown_storage=:sparse)
    h=1/n
    grid=VoronoiFVM.Grid(collect(0:h:1))

    eps::Vector{Float64}=[1,1.0e-1]

    physics=VoronoiFVM.Physics(
         num_species=2,
    reaction=function(f,u,node)
        f[1]=10*(u[1]-u[2])
        f[2]=10*(u[2]-u[1])
    end,

    flux=function(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=eps[1]*(u[1,1]-u[1,2])
        f[2]=eps[2]*(u[2,1]-u[2,2])
    end,

    storage=function(f,u,node)
        f[1]=u[1]
        f[2]=u[2]
    end
    )
    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])

    boundary_neumann!(sys,1,1,0.01)
    boundary_dirichlet!(sys,2,2,0.0)

    factory=VoronoiFVM.TestFunctionFactory(sys)
    tf1=testfunction(factory,[2],[1])
    tf2=testfunction(factory,[1],[2])

    U=unknowns(sys)
    inival=unknowns(sys)
    inival[2,:].=0.1
    inival[1,:].=0.1

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    control.damp_initial=0.1
    I1=0
    p=GridVisualizer(Plotter=Plotter,layout=(2,1))
    for xeps in [1.0,0.1,0.01]
        eps=[xeps,xeps]
        solve!(U,inival,sys,control=control)
        I1=integrate(sys,tf1,U)
        coord=coordinates(grid)
        inival.=U
        scalarplot!(p[1,1],grid,U[1,:])
        scalarplot!(p[2,1],grid,U[2,:])
        reveal(p)
        u5=U[5]
    end
    return I1[1]
end

function test()
    main(unknown_storage=:sparse) ≈ 0.01 &&
        main(unknown_storage=:dense) ≈ 0.01
end
end


text:
canvas.is_realized = true
set_gtk_property!(canvas, :hexpand, true)
set_gtk_property!(canvas, :vexpand, true)
set_gtk_property!(canvas, :double_buffered, false)

scrollableWindowMaterialList = nothing
scrollableWindowRoomList = nothing

lastCanvasDraw = time()

function generateMainGrid()
    grid = Grid()

    global scrollableWindowRoomList = ScrolledWindow(vexpand=true, hscrollbar_policy=Gtk.GtkPolicyType.NEVER)
    push!(scrollableWindowRoomList, roomList.tree)

    global scrollableWindowMaterialList = ScrolledWindow(vexpand=true, hscrollbar_policy=Gtk.GtkPolicyType.NEVER)
    push!(scrollableWindowMaterialList, materialList.tree)

    grid[1:6, 1] = menubar
    grid[1, 2:3] = scrollableWindowRoomList
    grid[2, 2:3] = canvas
    grid[3, 2:3] = toolList.tree
    grid[4, 2:3] = layersList.tree
    grid[5, 2:3] = subtoolList.tree
    grid[6, 2] = scrollableWindowMaterialList
    grid[6, 3] = materialFilterEntry

    return grid
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
export corconvcnp_1d

"""
    corconvcnp_1d(;
        receptive_field_μ::Float32,
        receptive_field_Σ::Float32,
        num_layers::Integer,
        num_channels::Integer,
        points_per_unit_μ::Float32,
        points_per_unit_Σ::Float32,
        margin::Float32=receptive_field
    )

Construct a CorrelatedConvCNP for one-dimensional data.

# Keywords
- `receptive_field_μ::Float32`: Width of the receptive field for the mean.
- `receptive_field_Σ::Float32`: Width of the receptive field for the kernel.
- `num_layers::Integer`: Number of layers of the CNN, excluding an initial
    and final pointwise convolutional layer to change the number of channels
    appropriately.
- `num_channels::Integer`: Number of channels of the CNN.
- `points_per_unit_μ::Float32`: Points per unit for the discretisation for the mean. See
     `UniformDiscretisation1D`.
- `points_per_unit_Σ::Float32`: Points per unit for the discretisation for the kernel. See
     `UniformDiscretisation1D`.
- `margin::Float32=receptive_field`: Margin for the discretisation. See
    `UniformDiscretisation1D`.
"""
function corconvcnp_1d(;
    receptive_field_μ::Float32,
    receptive_field_Σ::Float32,
    num_layers::Integer,
    num_channels::Integer,
    points_per_unit_μ::Float32,
    points_per_unit_Σ::Float32,
    margin::Float32=receptive_field
)
    dim_x = 1
    dim_y = 1
    scale_μ = 2 / points_per_unit_μ
    scale_Σ = 2 / points_per_unit_Σ
    return Model(
        Parallel(
            FunctionalCoder(
                UniformDiscretisation1D(points_per_unit_μ, margin),
                Chain(
                    set_conv(dim_y, scale_μ; density=true),
                    DeterministicLikelihood()
                )
            ),
            FunctionalCoder(
                UniformDiscretisation1D(points_per_unit_Σ, margin),
                Chain(
                    set_conv(dim_y, scale_Σ; density=true, pd=true),
                    DeterministicLikelihood()
                )
            ),
        ),
        Chain(
            Parallel(
                Chain(
                    build_conv(
                        receptive_field_μ,
                        num_layers,
                        num_channels,
                        points_per_unit =points_per_unit_μ,
                        dimensionality  =1,
                        num_in_channels =dim_y + 1,  # Account for density channel.
                        num_out_channels=dim_y
                    ),
                    set_conv(dim_y, scale_μ)
                ),
                Chain(
                    build_conv(
                        receptive_field_Σ,
                        num_layers,
                        num_channels,
                        points_per_unit =points_per_unit_Σ,
                        dimensionality  =2,
                        # Account for density and identity channel.
                        num_in_channels =dim_y + 2,
                        num_out_channels=dim_y
                    ),
                    set_conv(dim_y, scale_Σ; pd=true)
                )
            ),
            CorrelatedGaussianLikelihood()
        )
    )
end


text:
using Pkg

Pkg.activate(".") # make the current folder the working environmen

import Counter: Domain, API # import the sub-modules Domain and API

using .Domain, .API # load the sub-modules Domain and API

item = create("ar", 1000) # create an item

println(item)


text:
using Pkg; Pkg.activate(pwd())
using DEMC
using Test

@testset "DEMCz Sampling" begin
    # Write your own tests here.
    println("test with norm pdf in serial computation")
    include("example_normpdf.jl")
    println("test with norm pdf in parallel computation")
    include("example_normpdf_parallel.jl")
end

@testset "DEMCz Annealing" begin
    println("test annealing with quadratic optimization")
    include("test_anneal.jl")
    println("test annealing with quadratic optimization in parallel")
    include("test_anneal_parallel.jl")
end


text:
using VortexDistributions
using JLD2
using GLMakie
using Interpolations
using NearestNeighbors
using LinearAlgebra
using Distances
## Include util functions
include("../util.jl")

## Load psi simulation
# @load "data/box_rings.jld2" psi_ring1 psi_ring2 psi_ring3 X # make sure julia repl is in base directory
@load "data/box_vorts.jld2" psi_chaos1 psi_chaos2 psi_knots1 psi_knots2 psi_knots3 psi_knots4 psi_tubes1 psi_tubes2 psi_ringtube X
# @load "data/slab_tubes.jld2" psi_tubes1 psi_tubes2 psi_tubes3 X

## Set psi
psi = psi_tubes1;

## Plot iso surface
plot_iso(psi, X)

## Find vortices
@time vorts_3d = findvortices3D_itp(psi, X, 2);
num_vorts = length(vcat(vorts_3d...)[:,1])
print("Number of vortices found: ", length(vcat(vorts_3d...)[:,1]))

## Scatter vortices found
scatterVortsOnIso(vorts_3d)


text:
import Splittings: UniformMesh, advection!, BSpline

@testset "BSL advections" begin

  p = 5
  n1, n2 = 128, 128
  x1min, x1max = -10, 10
  x2min, x2max = -10, 10
  mesh1 = UniformMesh(x1min, x1max, n1; endpoint=false)
  mesh2 = UniformMesh(x2min, x2max, n2; endpoint=false)

  f  = zeros(Complex{Float64},(n1,n2))
  f .= exp.(-mesh1.points.^2) * transpose(exp.(-mesh2.points.^2))
  fᵗ = zeros(Complex{Float64},(n2,n1))

  dt =  0.5

  e = ones(Float64, n1)
  v = ones(Float64, n2)

  advection!(f,  mesh1,  v, n2, dt, BSpline(p))
  transpose!(fᵗ, f)
  advection!(fᵗ, mesh2,  e, n1, dt, BSpline(p))
  transpose!(f,  fᵗ)
  advection!(f,  mesh1, -v, n2, dt, BSpline(p))
  transpose!(fᵗ, f)
  advection!(fᵗ, mesh2, -e, n1, dt, BSpline(p))
  transpose!(f,  fᵗ)

  f0 =  exp.(-mesh1.points.^2) * transpose(exp.(-mesh2.points.^2))
  println( maximum( abs.(real(f) -f0)))

  @test real(f) ≈ f0 atol=1e-6

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UniformFolding(k, shuffle=true)

A method for creating `k` random folds from a spatial object
by sampling elements of the object with a uniform distribution.
Optionally `shuffle` the object before creating the folds.
"""
struct UniformFolding <: FoldingMethod
  k::Int
  shuffle::Bool
end

UniformFolding(k::Int) = UniformFolding(k, true)

function folds(object, method::UniformFolding)
  # retrieve parameters
  k, shuffle = method.k, method.shuffle

  # partition object
  p = partition(object, UniformPartition(k, shuffle))
  s = indices(p)
  n = length(p)

  function pair(i)
    # source and target subsets
    source = [1:i-1; i+1:n]
    target = [i]

    # indices within subsets
    sinds = reduce(vcat, s[source])
    tinds = reduce(vcat, s[target])

    sinds, tinds
  end

  (pair(i) for i in 1:n)
end


text:
# m2.1t.jl

using Pkg, DrWatson
@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking

Turing.setprogress!(false)

# Define the data

k = 6; n = 9;

# Define the model

@model function ppl2_1(n, k)
  θ ~ Beta(1, 1) # prior
  k ~ Binomial(n, θ) # model
  return k, θ
end;

# Use Turing mcmc

m2_1t = ppl2_1(n, k)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns2_1t = mapreduce(c -> sample(m2_1t, sampler, nsamples), chainscat, 1:nchains)
part2_1t = Particles(chns2_1t[:θ])

# End of m2.1t.jl


text:
struct StrayCopse end

struct StrayCopseResult{T, S}
    sc::Vector{Node{T}}
    d::Vector{S}
end

function fit(::StrayCopse, X, y, rng = Random.GLOBAL_RNG; max_depth = 6, min_node_records = 1,
        n_features_per_node = Int(floor(sqrt(size(X, 2)))), n_trees = 100)

    if n_features_per_node > size(X, 2)
        @warn "n_features_per_node is $n_features_per_node which is larger than the overall features number $(size(X, 2))."
        n_features_per_node = size(X, 2)
    end

    out_classes = unique(y)
    n_classes = length(out_classes)
    in_classes = collect(1:n_classes)
    d = Dict(zip(out_classes, in_classes))
    target = map(z -> d[z], y)

    T = eltype(X)
    sc = Node{T}[]
    nrow = size(X, 1)
    for i in 1:n_trees
        root = Node{T}()
        push!(sc, root)
        ids = sample(1:nrow, nrow)
        X1 = X[ids, :]
        target1 = target[ids]
        dtc = DecisionTreeContainer{T}(root, n_features_per_node, n_classes, max_depth, min_node_records)
        process_node(dtc, root, X1, target1, rng)
    end

    return StrayCopseResult(sc, out_classes)
end

function predict(scr::StrayCopseResult, X)
    d = scr.d

    n = length(d)
    vals = Vector{Int}(undef, n)
    res = Vector{eltype(d)}(undef, size(X, 1))
    for i in axes(X, 1)
        vals .= 0
        for tree in scr.sc
            vals[predict(tree, X, i)] += 1
        end
        res[i] = d[argmax(vals)]
    end

    return res
end


text:
###########################################################################
# Lazy wrapper around a PyObject to load a module on demand.

mutable struct LazyPyModule
    name::String
    o::PyObject
    LazyPyModule(n) = new(n, PyNULL())
end
PyObject(m::LazyPyModule) = ispynull(getfield(m, :o)) ? copy!(getfield(m, :o), pyimport(getfield(m, :name))) : getfield(m, :o)
pycall(m::LazyPyModule, args...; kws...) = pycall(PyObject(m), args...; kws...)
(m::LazyPyModule)(args...; kws...) = pycall(PyObject(m), PyAny, args...; kws...)
Base.Docs.doc(m::LazyPyModule) = Base.Docs.doc(PyObject(m))

# define each of these separately for Symbol and AbstractString to avoid method ambiguity:
Base.getproperty(m::LazyPyModule, x::Symbol) = getproperty(PyObject(m), x)
Base.getproperty(m::LazyPyModule, x::AbstractString) = getproperty(PyObject(m), x)
Base.setproperty!(m::LazyPyModule, x::Symbol, v) = setproperty!(PyObject(m), x, v)
Base.setproperty!(m::LazyPyModule, x::AbstractString, v) = setproperty!(PyObject(m), x, v)
PyCall.hasproperty(m::LazyPyModule, x::Symbol) = PyCall.hasproperty(PyObject(m), x)
PyCall.hasproperty(m::LazyPyModule, x::AbstractString) = PyCall.hasproperty(PyObject(m), x)

Base.propertynames(m::LazyPyModule) = propertynames(PyObject(m))
keys(m::LazyPyModule) = keys(PyObject(m))

# deprecated methods:
getindex(m::LazyPyModule, x) = getindex(PyObject(m), x)
setindex!(m::LazyPyModule, v, x) = setindex!(PyObject(m), v, x)
haskey(m::LazyPyModule, x) = haskey(PyObject(m), x)

###########################################################################
# Lazily load mplot3d modules.  This (slightly) improves load time of PyPlot,
# and it also allows PyPlot to load on systems where mplot3d is not installed.

const axes3D = LazyPyModule("mpl_toolkits.mplot3d.axes3d")
const art3D = LazyPyModule("mpl_toolkits.mplot3d.art3d")

"""
    using3D()

This function ensures that the `mplot3d` module is loaded for 3d
plotting.   This occurs automatically if you call any of the
3d plotting functions like `plot3D` or `surf`, but it may be
necessary to call this function manually if you are passing
`projection="3d"` explicitly to axes or subplot objects.
"""
using3D() = (PyObject(axes3D); nothing)

###########################################################################
# 3d plotting functions from mplot3d

export art3D, Axes3D, using3D, surf, mesh, bar3D, contour3D, contourf3D, plot3D, plot_surface, plot_trisurf, plot_wireframe, scatter3D, text2D, text3D, zlabel, zlim, zscale, zticks

const mplot3d_funcs = (:bar3d, :contour3D, :contourf3D, :plot3D, :plot_surface,
                       :plot_trisurf, :plot_wireframe, :scatter3D,
                       :text2D, :text3D, :view_init, :voxels)

for f in mplot3d_funcs
    fs = string(f)
    @eval @doc LazyHelp(axes3D,"Axes3D", $fs) function $f(args...; kws...)
        using3D() # make sure mplot3d is loaded
        if PyPlot.version <= v"3.4"
            ax = gca(projection="3d")
        else
            ax = gca()._subplotspec.get_geometry() == (1,1,0,0) ? plt."subplot"(projection="3d") : gca()
        end
        pycall(ax.$fs, PyAny, args...; kws...)
    end
end

@doc LazyHelp(axes3D,"Axes3D") Axes3D(args...; kws...) = pycall(axes3D."Axes3D", PyAny, args...; kws...)

# correct for annoying mplot3d inconsistency
@doc LazyHelp(axes3D,"Axes3D", "bar3d") bar3D(args...; kws...) = bar3d(args...; kws...)

# it's annoying to have xlabel etc. but not zlabel
const zlabel_funcs = (:zlabel, :zlim, :zscale, :zticks)
for f in zlabel_funcs
    fs = string("set_", f)
    @eval @doc LazyHelp(axes3D,"Axes3D", $fs) function $f(args...; kws...)
        using3D() # make sure mplot3d is loaded
        if PyPlot.version <= v"3.4"
            ax = gca(projection="3d")
        else
            ax = gca()._subplotspec.get_geometry() == (1,1,0,0) ? plt."subplot"(projection="3d") : gca()
        end
        pycall(ax.$fs, PyAny, args...; kws...)
    end
end

# export Matlab-like names

function surf(Z::AbstractMatrix; kws...)
    plot_surface([1:size(Z,1);]*ones(1,size(Z,2)),
                 ones(size(Z,1))*[1:size(Z,2);]', Z; kws...)
end

@doc LazyHelp(axes3D,"Axes3D", "plot_surface") function surf(X, Y, Z::AbstractMatrix, args...; kws...)
    plot_surface(X, Y, Z, args...; kws...)
end

function surf(X, Y, Z::AbstractVector, args...; kws...)
    plot_trisurf(X, Y, Z, args...; kws...)
end

@doc LazyHelp(axes3D,"Axes3D", "plot_wireframe") mesh(args...; kws...) = plot_wireframe(args...; kws...)

function mesh(Z::AbstractMatrix; kws...)
    plot_wireframe([1:size(Z,1);]*ones(1,size(Z,2)),
                   ones(size(Z,1))*[1:size(Z,2);]', Z; kws...)
end


text:
# irreducible unitary representations of the groups SO(3) and SE(3) and representations of
# their Lie algebras.
# For SO(3), see
# Varshalovich, D. A., Moskalev, A. N. & Khersonskii, V. K. Quantum Theory of Angular
# Momentum. 1988. World Scientific, Singapore. ISBN: 9971-5-0107-4.
# For SO(3) and SE(3), see §12.9 and 12.12 of
# Chirikjian, G. S. Stochastic Models, Information Theory, and Lie Groups, Volume 2. 2012.
# ISBN: 978-0-8176-4943-2. doi: 10.1007/978-0-8176-4944-9.

# diagonals of infinitesimal rotation elements of Varshalovich §4.18.3.
function representation_diag(::typeof(so3), ::typeof(E1), ℓ, ::Val{i}) where {i}
    ℓ2 = Int(2ℓ)
    absi = abs(i)
    isone(absi) && return @. im * (-sqrt((ℓ2:-1:1) * (1:ℓ2)) / 2)
    return zeros(complex(float(eltype(ℓ))), ℓ2 + 1 - absi)
end
function representation_diag(::typeof(so3), ::typeof(E2), ℓ, ::Val{i}) where {i}
    ℓ2 = Int(2ℓ)
    absi = abs(i)
    isone(absi) && return @. i * sqrt((ℓ2:-1:1) * (1:ℓ2)) / 2
    return zeros(float(eltype(ℓ)), ℓ2 + 1 - absi)
end
function representation_diag(::typeof(so3), ::typeof(E3), ℓ, ::Val{i}) where {i}
    iszero(i) && return Vector(im .* (ℓ:-1:(-ℓ)))
    return zeros(complex(eltype(ℓ)), Int(2ℓ) + 1 - abs(i))
end

"""
    representation_block(::typeof(so3), X, ℓ)

Compute the size `(2ℓ+1, 2ℓ+1)` block ``ℓ`` of the representation of the element ``X`` of
the element of the Lie algebra ``𝔰𝔬(3)``.

The representation is defined as

```math
u^ℓ(X) = \\lim_{t → 0} \\frac{d}{dt} U^ℓ(\\exp(t X)),
```

where ``U^ℓ(R)`` is the representation of the rotation ``R``, and ``\\exp(R)`` is the group
exponential map on ``\\mathrm{SO}(3)``.
"""
representation_block(::typeof(so3), X, ℓ)

function representation_block(::typeof(so3), Ei::BasisVector, ℓ)
    dl = representation_diag(so3, Ei, ℓ, Val(-1))
    d = representation_diag(so3, Ei, ℓ, Val(0))
    du = representation_diag(so3, Ei, ℓ, Val(1))
    return Tridiagonal(dl, d, du)
end
function representation_block(::typeof(so3), ::typeof(E3), ℓ)
    return Diagonal(representation_diag(so3, E3, ℓ, Val(0)))
end
function representation_block(::typeof(so3), Xⁱ::AbstractVector{<:Real}, ℓ)
    # because representation_element returns a result with real eltype Float64, we use `T`
    # to avoid promoting if the eltype of Xⁱ is lower precision
    T = complex(float(eltype(Xⁱ)))
    uE1 = representation_block(so3, E1, ℓ)
    uE2 = representation_block(so3, E2, ℓ)
    uE3 = representation_block(so3, E3, ℓ)
    return Tridiagonal(
        Xⁱ[1] .* T.(uE1.dl) .+ Xⁱ[2] .* T.(uE2.dl),
        Xⁱ[3] .* T.(uE3.diag),
        Xⁱ[1] .* T.(uE1.du) .+ Xⁱ[2] .* T.(uE2.du),
    )
end
function representation_block(::typeof(so3), X::AbstractMatrix, ℓ)
    # represent the vector X as an actual Vector using the basis
    Xⁱ = Manifolds.vee(SO3, Manifolds.Identity(SO3, one(X)), X)
    return representation_block(so3, Xⁱ, ℓ)
end

"""
    representation_block(::typeof(SO3), p, ℓ) -> Matrix{Complex}

Compute the size `(2ℓ+1, 2ℓ+1)` block ``ℓ`` of the representation of the rotation ``p``.

This block is identical to the ``D^ℓ`` block of the Wigner ``D``-matrix using the
conventions of [^Varshalovich1988].

[^Varshalovich1988]: Varshalovich, D. A., Moskalev, A. N. & Khersonskii, V. K. Quantum Theory of Angular
    Momentum. 1988. World Scientific, Singapore. ISBN: 9971-5-0107-4.
# Examples

To create the "little-``d``" block ``d^ℓ(\\frac{π}{4}) = D^ℓ(0, \\frac{π}{4}, 0)``,

```julia-repl
julia> using EnsembleInference, Rotations

julia> EnsembleInference.representation_block(SO3, RotZXZ(0, π/2, 0), 1)
3×3 Matrix{ComplexF64}:
  0.5+0.0im              0.0-0.707107im  -0.5+0.0im
  0.0-0.707107im  1.2175e-16+0.0im        0.0-0.707107im
 -0.5+0.0im              0.0-0.707107im   0.5+0.0im

julia> EnsembleInference.representation_block(SO3, RotZXZ(0, π/2, 0), 2)
5×5 Matrix{ComplexF64}:
      0.25+0.0im   0.0-0.5im          -0.612372+0.0im           0.0+0.5im               0.25+0.0im
       0.0-0.5im  -0.5-0.0im                0.0-3.24998e-16im  -0.5-0.0im                0.0+0.5im
 -0.612372-0.0im   0.0-2.41689e-16im       -0.5-0.0im           0.0-3.62534e-16im  -0.612372-0.0im
       0.0+0.5im  -0.5+0.0im                0.0-2.03875e-16im  -0.5+0.0im                0.0-0.5im
      0.25+0.0im   0.0+0.5im          -0.612372+0.0im           0.0-0.5im               0.25+0.0im
```
"""
function representation_block(::typeof(SO3), p, ℓ)
    X = Manifolds.group_log(SO3, p)
    u = representation_block(so3, X, ℓ)
    return exp(Matrix(u))
end


text:
using FMMLIB2D
using SpecialFunctions
using LinearAlgebra
using Random
using Test

Random.seed!(0)

# TODO: Compute gradient and hessian

function direct_self(source, charge, dipvec, dipstr, zk)
    pot = zero(charge)
    N = size(source, 2)
    for i=1:N
        for j=1:N
            if i == j
                continue
            end
            r1 = source[1,j]-source[1,i]
            r2 = source[2,j]-source[2,i]
            r = sqrt(r1^2 + r2^2)
            rdotq = r1*dipvec[1,j] + r2*dipvec[2,j]
            pot[i] += charge[j]*besselh(0, zk*r)*im/4 -
                dipstr[j]*rdotq*1im/4*besselh(1, zk*r)*zk/r
        end
    end
    return pot
end
function direct_targ(source, charge, dipvec, dipstr, target, zk)
    M = size(target, 2)
    N = size(source, 2)
    pot = zeros(typeof(charge[1]), M)
    for i=1:M
        for j=1:N
            r1 = source[1,j]-target[1,i]
            r2 = source[2,j]-target[2,i]
            r = sqrt(r1^2 + r2^2)
            rdotq = r1*dipvec[1,j] + r2*dipvec[2,j]
            pot[i] += charge[j]*besselh(0, zk*r)*im/4 -
                dipstr[j]*rdotq*1im/4*besselh(1, zk*r)*zk/r
        end
    end
    return pot
end

@testset "Helmholtz FMM 2D (complex)" begin
    iprec = 5
    N = 40
    M = 30
    zk = rand() + 1im*rand()
    source = rand(2, N)
    charge = rand(N) + 1im*rand(N)
    dipstr = rand(N) + 1im*rand(N)
    dipvec = rand(2, N)
    target = rand(2, M)
    ifcharge = 1
    ifdipole = 1
    ifpot = 1
    ifgrad = 1
    ifhess = 1
    ifpottarg = 1
    ifgradtarg = 1
    ifhesstarg = 1

    refpottarg = direct_targ(source, charge, dipvec, dipstr, target, zk)
    refpot = direct_self(source, charge, dipvec, dipstr, zk)

    @testset "Direct interface: hfmm2dparttarg" begin
        pot, grad, hess, pottarg, gradtarg, hesstarg =
            hfmm2dparttarg(iprec, zk, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
                           target, ifpottarg, ifgradtarg, ifhesstarg)
        pot_relerr = norm(refpot-pot, Inf) / norm(refpot, Inf)
        pottarg_relerr = norm(refpottarg-pottarg, Inf) / norm(refpottarg, Inf)
        @test pot_relerr < 1e-15
        @test pottarg_relerr < 1e-15
    end

    @testset "Keyword interface" begin
        U = hfmm2d(zk=zk, source=source, charge=charge, dipstr=dipstr, dipvec=dipvec)
        @test norm(U.pot-refpot, Inf) / norm(refpot, Inf) < 1e-14
        U = hfmm2d(zk=zk, source=source, charge=charge, dipstr=dipstr, dipvec=dipvec, target=target)
        @test norm(U.pottarg-refpottarg, Inf) / norm(refpot, Inf) < 1e-14
    end

    @testset "Tolerance test" begin
        N = 1000
        source = rand(2, N)
        charge = rand(N) + 1im*rand(N)
        dipstr = rand(N) + 1im*rand(N)
        dipvec = rand(2, N)
        refpot = direct_self(source, charge, dipvec, dipstr, zk)
        for tolexp=-14:-1
            TOL = 0.49*10.0^tolexp
            U = hfmm2d(zk=zk, source=source, charge=charge, tol=TOL, dipstr=dipstr, dipvec=dipvec)
            relerr = norm(U.pot-refpot) / norm(refpot)
            @test relerr < TOL
        end
    end
end


text:
erf1arg = [
(:erf, :(exp(-abs2(x))*2/sqrt(pi)), (-Inf,Inf)),     # \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt
(:erfc, :(-exp(-abs2(x))*2/sqrt(pi)), (-Inf,Inf)),   # 1-erf(x)
(:erfcx, :(2y.*x-2/sqrt(pi)), (-Inf,Inf)),           # erfc(x)*exp(x^2)
(:erfi, :(exp(abs2(x))*2/sqrt(pi)), (-Inf,Inf)),     # -i*erf(ix)
(:dawson, :(-2y.*x+1), (-Inf,Inf)),                  # \frac{\sqrt{\pi}}{2} e^{-x^2} erfi(x).
(:erfinv, :(exp(abs2(y))*sqrt(pi)/2), (-1,1)),       # erf(erfinv(x)) = x
(:erfcinv, :(-exp(abs2(y))*sqrt(pi)/2), (0,2)),      # erfc(erfcinv(x)) = x
]

for (f,g,r) in erf1arg
    @eval @primitive $f(x),dy,y  (dy.*($g))
    addtest1(f,r)
end


text:
const study_time = 25*60

const break_times = (15 * 60, 5 * 60)

"function pomodoro_stuff( n_cycles::Integer; cycle_time::Integer = study_time, rest_times::Tuple = break_times)"
function pomodoro_stuff( n_cycles::Integer; cycle_time::Integer = study_time, rest_times::Tuple = break_times)

    times = 0
    long_rest, short_rest = rest_times

    function one_cycle(times)
        @showprogress 1 "Study Time:"  for i in 1:cycle_time
            sleep(1)
        end

        @uhlarm @printf("ok dude; %d cycles \n", times+1)
    end

    function break_time(times)
        if times % 4 ==0
            rest_time = long_rest
            println("longer break \n ")
        else
            println("smaller break \n ")
            rest_time = short_rest
        end
        @showprogress 1 "Break Time:"  for i in 1:rest_time
            sleep(1)
        end

        @uhlarm println("end of break \n")
    end

    function responsefrom_me()

    println("Back to work?")
        response = readline(STDIN);

    if response == "yes"
        println("oh yeah \n")
    else
        error("gotta rest too motherfucker")
    end

    end

    while times < n_cycles
        one_cycle(times)
        times+=1
        break_time(times)
        responsefrom_me()
    end

    if times == n_cyles
        error("too much coffe, arguably")
    end

end

max_stars_repo_name: marcelovmaciel/Pomodoro

max_issues_repo_name: marcelovmaciel/Pomodoro

max_forks_repo_name: marcelovmaciel/Pomodoro


text:
function create_ray(x1::T, y1::T, z1::T, x2::T, y2::T, z2::T) where {T<:AbstractFloat}
    dx = x2 - x1
    dy = y2 - y1
    dz = z2 - z1
    l = sqrt((dx^dx) + (dy*dy) + (dz*dz))
    dx = dx / l
    dy = dy / l
    dz = dz / l
    ray = Ray(x1, y1, z1, dx, dy, dz)
    return ray
end

function create_rays(minx::AbstractFloat, miny::AbstractFloat, maxx::AbstractFloat, maxy::AbstractFloat, height::AbstractFloat, dx::AbstractFloat, dy::AbstractFloat, dz::AbstractFloat, n::Int)
    x = rand(Uniform(minx, maxx), n)
    y = rand(Uniform(miny, maxy), n)
    rays = Ray[]
    for i=1:n
        push!(rays, Ray(x[i], y[i], height, dx, dy, dz))
    end
    return rays
end

function create_ray(walls::Walls, dx::AbstractFloat, dy::AbstractFloat, dz::AbstractFloat)
    minx, maxx = walls.west + SMALL_NUM, walls.east - SMALL_NUM
    miny, maxy = walls.south + SMALL_NUM, walls.north - SMALL_NUM # check this though!
    height = walls.top - SMALL_NUM
    x = rand(Uniform(minx, maxx))
    y = rand(Uniform(miny, maxy))
    ray = Ray(x, y, height, dx, dy, dz)
    return ray
end

function create_rays(walls::Walls, dx::AbstractFloat, dy::AbstractFloat, dz::AbstractFloat, n::Int)
    minx, maxx = walls.west + SMALL_NUM, walls.east - SMALL_NUM
    miny, maxy = walls.south + SMALL_NUM, walls.north - SMALL_NUM # check this though!
    height = walls.top - SMALL_NUM
    rays = create_rays(minx, miny, maxx, maxy, height, dx, dy, dz, n)
    return rays
end

function create_coords(camera::Camera{T,S,A}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}}
    # construct an array of coordinate, just for the shape...
    coords = Array{Coord{S}}(undef, camera.xResolution, camera.yResolution)
    for x::S=1:camera.xResolution
        for y::S=1:camera.yResolution
            coords[x,y] = Coord(x, y)
        end
    end
    return coords
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
using Mimi

@defcomp neteconomy begin
    regions = Index()

    YNET = Variable(index=[time, regions]) # Output net of damages equation (trillions 2005 USD per year)
    Y = Variable(index=[time, regions]) # Gross world product net of abatement and damages (trillions 2005 USD per year)
    I = Variable(index=[time, regions]) # Investment (trillions 2005 USD per year)
    C = Variable(index=[time, regions]) # Consumption (trillions 2005 US dollars per year)
    CPC = Variable(index=[time, regions]) # Per capita consumption (thousands 2005 USD per year)

    YGROSS = Parameter(index=[time, regions]) # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    DAMFRAC = Parameter(index=[time, regions]) # Damages as fraction of gross output
    DAMAGES = Parameter(index=[time, regions]) # Damages (trillions 2005 USD per year)
    ABATECOST = Parameter(index=[time, regions]) # Cost of emissions reductions  (trillions 2005 USD per year)
    S = Parameter(index=[time, regions]) # Gross savings rate as fraction of gross world product
    l = Parameter(index=[time, regions]) # Level of population and labor

    function run_timestep(p, v, d, t)

        #Define function for YNET
        for r in d.regions
            if is_first(t)
                v.YNET[t,r] = p.YGROSS[t,r]/(1+p.DAMFRAC[t,r])
            else
                v.YNET[t,r] = p.YGROSS[t,r] - p.DAMAGES[t,r]
            end
        end

        #Define function for Y
        for r in d.regions
            v.Y[t,r] = v.YNET[t,r] - p.ABATECOST[t,r]
        end

        #Define function for I
        for r in d.regions
            v.I[t,r] = p.S[t,r] * v.Y[t,r]
        end

        #Define function for C
        for r in d.regions
            if t.t != 60
                v.C[t,r] = v.Y[t,r] - v.I[t,r]
            else
                v.C[t,r] = v.C[t-1, r]
            end
        end

        #Define function for CPC
        for r in d.regions
            v.CPC[t,r] = 1000 * v.C[t,r] / p.l[t,r]
        end
    end
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo
using Plots

q1 = nominal_configuration(box)
q1 = [0.0; 0.5; 0.0]
v1 = zeros(box.nq)
u1 = zeros(box.nu)

T = 1
h = 0.1

u = range(0.0, stop=2.0, length=1000)
κ = [10.0^(-i) for i = 2:8]

x = []
∂x∂u = []
κtol = []

for (i, κi) in enumerate(κ)
    xi = []
    ∂xi∂u = []
    κitol = []
    for u in ux
        s = Simulator(box, T, h=h,
            solver_opts=InteriorPointOptions(
                undercut=2.0,
                r_tol=1e-10,
                κ_tol=κi,
                diff_sol=true,
                verbose=false))

        step!(s, q1, v1, [u; 0.0; 0.0], 1)

        push!(xi, s.traj.q[end][1])
        push!(∂xi∂u, s.grad.∂q3∂u1[end][1, 1])
        push!(κitol, s.ip.κ[1])
    end
    push!(x, xi)
    push!(∂x∂u, ∂xi∂u)
    push!(κtol, κitol)
end

plt_x = plot(; xlabel="F", ylabel="Δx", label="");
for (i, xi) in enumerate(x)
    plt_x = plot!(u, xi, label="κ=$(scn(κ[i]))", legend=:topleft)
end
display(plt_x)

plt_∂x∂u = plot(; xlabel="F", ylabel="∂x∂F", label="");
for (i, ∂xi∂u) in enumerate(∂x∂u)
    plt_∂x∂u = plot!(u, ∂xi∂u, label="κ=$(scn(κ[i]))", legend=:topleft)
end
display(plt_∂x∂u)

# ## Visualizer
vis = Visualizer()
render(vis)

# ## Visualize
visualize!(vis, s)


text:
using PkgBenchmark
benchmarkpkg(
    dirname(@__DIR__),
    BenchmarkConfig(
        env = Dict(
            "JULIA_NUM_THREADS" => "4",
        ),
    ),
    resultfile = joinpath(@__DIR__, "result.json")
)


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames, Plots, Random
gr(legend=false)

outputdir = joinpath(tempdir(),"juliaplot" * randstring())
mkdir(outputdir)
filename = joinpath(outputdir, "plot.png")

ui = Builder(filename=joinpath(@__DIR__, "ui.glade"))

if isinteractive()
    set_gtk_property!(ui["preview"], :file, "Julia_prog_language_logo.svg")
else
    set_gtk_property!(ui["preview"], :file, joinpath(dirname(PROGRAM_FILE), "Julia_prog_language_logo.svg"))
end

showall(ui["win"])

mutable struct Variables
    eqn::String
	x::String
	xs::Float64
	xf::Float64
	y::String
    ys::Float64
    yf::Float64
    t::String
    ts::Float64
    tf::Float64
    xlog::Bool
    ylog::Bool
end

variables = Variables(get_gtk_property(ui["eqn"], :text, String),
          get_gtk_property(ui["x"], :text, String),
          eval(Meta.parse(get_gtk_property(ui["xs"], :text, String))),
          eval(Meta.parse(get_gtk_property(ui["xf"], :text, String))),
          get_gtk_property(ui["y"], :text, String),
          eval(Meta.parse(get_gtk_property(ui["ys"], :text, String))),
          eval(Meta.parse(get_gtk_property(ui["yf"], :text, String))),
          get_gtk_property(ui["t"], :text, String),
          eval(Meta.parse(get_gtk_property(ui["ts"], :text, String))),
          eval(Meta.parse(get_gtk_property(ui["tf"], :text, String))),
          get_gtk_property(ui["xlog"], :active, Bool),
          get_gtk_property(ui["ylog"], :active, Bool))

function update_variables()
    variables.eqn = get_gtk_property(ui["eqn"], :text, String)
    variables.x = get_gtk_property(ui["x"], :text, String)
    variables.xs = eval(Meta.parse(get_gtk_property(ui["xs"], :text, String)))
    variables.xf = eval(Meta.parse(get_gtk_property(ui["xf"], :text, String)))
    variables.y = get_gtk_property(ui["y"], :text, String)
    variables.ys = eval(Meta.parse(get_gtk_property(ui["ys"], :text, String)))
    variables.yf = eval(Meta.parse(get_gtk_property(ui["yf"], :text, String)))
    variables.t = get_gtk_property(ui["t"], :text, String)
    variables.ts = eval(Meta.parse(get_gtk_property(ui["ts"], :text, String)))
    variables.tf = eval(Meta.parse(get_gtk_property(ui["tf"], :text, String)))
    variables.xlog = get_gtk_property(ui["xlog"], :active, Bool)
    variables.ylog = get_gtk_property(ui["ylog"], :active, Bool)

    return nothing
end

function plot_option!()
    plot!(title=get_gtk_property(ui["title"], :text, String))
    plot!(xlabel=get_gtk_property(ui["xlabel"], :text, String))
    plot!(ylabel=get_gtk_property(ui["ylabel"], :text, String))
    plot!(zlabel=get_gtk_property(ui["zlabel"], :text, String))
    if variables.xlog
        plot!(xscale=:log10)
    end
    if variables.ylog
        plot!(yscale=:log10)
    end
end

function genfn_single()
    if occursin("=", variables.eqn)
        eqn = split(variables.eqn, "=")[2]
    else
        eqn = variables.eqn
    end
    f = eval(Meta.parse(variables.x * "->" * eqn))

    return f
end

function genfn_double()
    x = variables.xs:variables.xf
    y = variables.ys:variables.yf
    if occursin("=", variables.eqn)
        eqn = split(variables.eqn, "=")[2]
    else
        eqn = variables.eqn
    end
    f = eval(Meta.parse( "(" * variables.x * "," * variables.y * ")" * "->" * eqn))

    return x, y, f
end

function genfn_para()
    if length(split(variables.eqn, ",")) == 2
        de = split(variables.eqn, ",")
        if occursin("=", de[1])
            e1 = split(de[1], "=")[2]
        else
            e1 = de[1]
        end

        if occursin("=", de[2])
            e2 = split(de[2], "=")[2]
        else
            e2 = de[2]
        end

        f = eval(Meta.parse(variables.t * "->" * e1))
        g = eval(Meta.parse(variables.t * "->" * e2))

        return f, g
    else
        de = split(variables.eqn, ",")
        if occursin("=", de[1])
            e1 = split(de[1], "=")[2]
        else
            e1 = de[1]
        end

        if occursin("=", de[2])
            e2 = split(de[2], "=")[2]
        else
            e2 = de[2]
        end

        if occursin("=", de[3])
            e3 = split(de[3], "=")[2]
        else
            e3 = de[3]
        end

        f = eval(Meta.parse(variables.t * "->" * e1))
        g = eval(Meta.parse(variables.t * "->" * e2))
        h = eval(Meta.parse(variables.t * "->" * e3))

        return f, g, h
    end
end

function plot2d_button()
    update_variables()
    f = genfn_single()
    plot(x-> Base.invokelatest(f,x), variables.xs, variables.xf)
    plot_option!()
    savefig(filename)

    set_gtk_property!(ui["preview"], :file, filename)
    println("plot 2d function")

    return nothing
end

function contour_button()
    update_variables()
    x, y, f = genfn_double()
    contour(x, y, (x,y) -> Base.invokelatest(f, x, y), fill=true, colorbar=true)
    plot_option!()
    savefig(filename)

    set_gtk_property!(ui["preview"], :file, filename)
    println("contour")

    return nothing
end

function para2d_button()
    update_variables()
    f, g = genfn_para()
    plot(t1 -> Base.invokelatest(f,t1), t2 -> Base.invokelatest(g,t2), variables.ts, variables.tf)
    plot_option!()
    savefig(filename)
    set_gtk_property!(ui["preview"], :file, filename)
    println("parametric function 2d")

    return nothing
end

function plot3d_button()
    update_variables()
    x, y, f = genfn_double()
    surface(x, y, (x,y) -> Base.invokelatest(f, x, y))
    plot_option!()
    savefig(filename)

    set_gtk_property!(ui["preview"], :file, filename)
    println("plot 3d function")

    return nothing
end

function para3d_button()
    update_variables()
    f, g, h = genfn_para()
    plot3d(t -> Base.invokelatest(f,t),
         t -> Base.invokelatest(g,t),
         t -> Base.invokelatest(h,t),
         variables.ts, variables.tf)
    plot_option!()
    savefig(filename)
    set_gtk_property!(ui["preview"], :file, filename)
    println("parametric function 3d")

    return nothing
end

function savefig_button()
    dstname = save_dialog("Save figure")
    cp(filename, dstname, remove_destination=true)

    return nothing
end

signal_connect(x -> plot2d_button(), ui["plot2d"], "clicked")
signal_connect(x -> contour_button(), ui["contour"], "clicked")
signal_connect(x -> para2d_button(), ui["para2d"], "clicked")
signal_connect(x -> plot3d_button(), ui["plot3d"], "clicked")
signal_connect(x -> para3d_button(), ui["para3d"], "clicked")
signal_connect(x -> savefig_button(), ui["savefig"], "clicked")

if !isinteractive()
    c = Condition()
    signal_connect(ui["win"], :destroy) do widget
        notify(c)
    end
    wait(c)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
using ClangCompiler
using Test

@testset "function call" begin
    src = joinpath(@__DIR__, "code", "main.cpp")
    args = get_compiler_args()
    haskey(ENV, "CI") && push!(args, "-v")
    irgen = IRGenerator(src, args)
    f = lookup_function(get_module(irgen), "main")
    @eval main() = $(call_function(f, Cint))
    @test main() == 42
    dispose(irgen)
end


text:
@testset "ShearX" begin
    @test typeof(@inferred(ShearX(1))) <: ShearX <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test_throws MethodError ShearX()
        @test_throws MethodError ShearX(:a)
        @test_throws MethodError ShearX([:a])
        @test_throws ArgumentError ShearX(3:1)
        @test_throws ArgumentError ShearX(-71:1)
        @test_throws ArgumentError ShearX(1:71)
        @test_throws ArgumentError ShearX([2,71,-4])
        @test_throws ArgumentError ShearX(71)
        @test @inferred(ShearX(70)) === ShearX(70:70)
        @test @inferred(ShearX(-70)) === ShearX(-70:-70)
        @test @inferred(ShearX(0.7)) === ShearX(0.7:1:0.7)
        @test str_show(ShearX(0.7)) == "Augmentor.ShearX(0.7)"
        @test str_showconst(ShearX(0.7)) == "ShearX(0.7)"
        @test str_showcompact(ShearX(0.7)) == "ShearX 0.7 degree"
        @test @inferred(ShearX(10)) === ShearX(10:10)
        @test str_show(ShearX(10)) == "Augmentor.ShearX(10)"
        @test str_showconst(ShearX(10)) == "ShearX(10)"
        @test str_showcompact(ShearX(10)) == "ShearX 10 degree"
        op = @inferred(ShearX(-1:1))
        @test str_show(op) == "Augmentor.ShearX(-1:1)"
        @test str_showconst(op) == "ShearX(-1:1)"
        @test str_showcompact(op) == "ShearX by ϕ ∈ -1:1 degree"
        op = @inferred(ShearX([2,30]))
        @test op.degree == [2,30]
        @test str_show(op) == "Augmentor.ShearX([2, 30])"
        @test str_showconst(op) == "ShearX([2, 30])"
        @test str_showcompact(op) == "ShearX by ϕ ∈ [2, 30] degree"
    end
    @testset "randparam" begin
        @test @inferred(Augmentor.randparam(ShearX(45), rect)) === 45.0
        @test @inferred(Augmentor.randparam(ShearX(20.0), rect)) === 20.0
        @test @inferred(Augmentor.randparam(ShearX(1:10), rect)) in 1:10
        @test_throws MethodError Augmentor.toaffinemap(ShearX(45), rect)
        @test @inferred(Augmentor.toaffinemap(ShearX(45), rect, 45)) ≈ AffineMap([1. 0.; -1. 1.], [0.,1.5])
        @test @inferred(Augmentor.toaffinemap(ShearX(-60:60), rect, 45)) ≈ AffineMap([1. 0.; -1. 1.], [0.,1.5])
        @test @inferred(Augmentor.toaffinemap(ShearX(-45), rect, -45)) ≈ AffineMap([1. 0.; 1. 1.], [0.,-1.5])
        @test @inferred(Augmentor.toaffinemap(ShearX(-60:1:60), rect, -45)) ≈ AffineMap([1. 0.; 1. 1.], [0.,-1.5])
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(ShearX(10), nothing)
        @test Augmentor.supports_eager(ShearX) === false
        # TODO: actual content tests (maybe test_reference)
        img_out1 = @inferred Augmentor.applyeager(ShearX(45), square)
        img_out2 = @inferred Augmentor.applyeager(ShearX(-45), square)
        @test axes(img_out1) == (1:3, 0:4)
        @test axes(img_out1) == axes(img_out2)
        imgs = [
            (square),
            (view(square, :, :)),
            (Augmentor.prepareaffine(square)),
            (OffsetArray(square, -1, -1)),
            (view(square, IdentityRange(1:3), IdentityRange(1:3))),
        ]
        @testset "fixed parameter" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(ShearX(45), img_in))
                @test parent(res) == parent(img_out1)
                @test typeof(res) == typeof(img_out1)
                res = @inferred(Augmentor.applyeager(ShearX(-45), img_in))
                @test parent(res) == parent(img_out2)
                @test typeof(res) == typeof(img_out2)
                # test same with tuple of images
                res1, res2 = @inferred(Augmentor.applyeager(ShearX(45), (img_in, N0f8.(img_in))))
                @test parent(res1) == parent(img_out1)
                @test parent(res2) == parent(img_out1)
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: OffsetArray{N0f8}
                res1, res2 = @inferred(Augmentor.applyeager(ShearX(-45), (img_in, N0f8.(img_in))))
                @test parent(res1) == parent(img_out2)
                @test parent(res2) == parent(img_out2)
                @test typeof(res1) == typeof(img_out2)
                @test typeof(res2) <: OffsetArray{N0f8}
            end
            # check that the affine map is computed for each image
            res1, res2 = @inferred(Augmentor.applyeager(ShearX(45), (square, OffsetArray(square,-5,-5))))
            @test collect(res1) == collect(res2)
            @test axes(res1) != axes(res2)
        end
        @testset "random parameter" begin
            for img_in in imgs
                res1, res2 = @inferred(Augmentor.applyeager(ShearX(1:60), (img_in, N0f8.(img_in))))
                # make sure same angle is used
                @test res1 == res2
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: OffsetArray{N0f8}
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(ShearX) === true
        @test_throws MethodError Augmentor.applyaffine(ShearX(45), nothing)
        wv = @inferred Augmentor.applyaffine(ShearX(45), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test axes(wv) == (1:3,0:4)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        wv2 = @inferred Augmentor.applyaffine(ShearX(-45), wv)
        @test parent(wv).itp.coefs === square
        @test wv2[1:3,1:3] == square
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(ShearX) === true
        @test_throws MethodError Augmentor.applyaffineview(ShearX(45), nothing)
        wv = @inferred Augmentor.applyaffineview(ShearX(45), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test axes(wv) == (1:3,0:4)
        wv2 = @inferred Augmentor.applyaffineview(ShearX(-45), wv)
        @test wv2[1:3,1:3] == square
        @test typeof(wv2) <: SubArray{eltype(square),2}
        @test typeof(parent(wv2)) <: InvWarpedView
        @test parent(parent(wv2)).itp.coefs === square
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(ShearX(45)) === true
        wv = @inferred Augmentor.applylazy(ShearX(45), square)
        @test parent(wv).itp.coefs === square
        @test axes(wv) == (1:3,0:4)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        wv2 = @inferred Augmentor.applylazy(ShearX(-45), wv)
        @test parent(wv).itp.coefs === square
        @test wv2[1:3,1:3] == square
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "view" begin
        @test Augmentor.supports_view(ShearX) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(ShearX) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(ShearX) === false
    end
end

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

@testset "ShearY" begin
    @test typeof(@inferred(ShearY(1))) <: ShearY <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test_throws MethodError ShearY()
        @test_throws MethodError ShearY(:a)
        @test_throws MethodError ShearY([:a])
        @test_throws ArgumentError ShearY(3:1)
        @test_throws ArgumentError ShearY(-71:1)
        @test_throws ArgumentError ShearY(1:71)
        @test_throws ArgumentError ShearY([2,71,-4])
        @test_throws ArgumentError ShearY(71)
        @test @inferred(ShearY(70)) === ShearY(70:70)
        @test @inferred(ShearY(-70)) === ShearY(-70:-70)
        @test @inferred(ShearY(0.7)) === ShearY(0.7:1:0.7)
        @test str_show(ShearY(0.7)) == "Augmentor.ShearY(0.7)"
        @test str_showconst(ShearY(0.7)) == "ShearY(0.7)"
        @test str_showcompact(ShearY(0.7)) == "ShearY 0.7 degree"
        @test @inferred(ShearY(10)) === ShearY(10:10)
        @test str_show(ShearY(10)) == "Augmentor.ShearY(10)"
        @test str_showconst(ShearY(10)) == "ShearY(10)"
        @test str_showcompact(ShearY(10)) == "ShearY 10 degree"
        op = @inferred(ShearY(-1:1))
        @test str_show(op) == "Augmentor.ShearY(-1:1)"
        @test str_showconst(op) == "ShearY(-1:1)"
        @test str_showcompact(op) == "ShearY by ψ ∈ -1:1 degree"
        op = @inferred(ShearY([2,30]))
        @test op.degree == [2,30]
        @test str_show(op) == "Augmentor.ShearY([2, 30])"
        @test str_showconst(op) == "ShearY([2, 30])"
        @test str_showcompact(op) == "ShearY by ψ ∈ [2, 30] degree"
    end
    @testset "randparam" begin
        @test @inferred(Augmentor.randparam(ShearY(45), rect)) === 45.0
        @test @inferred(Augmentor.randparam(ShearY(20.0), rect)) === 20.0
        @test @inferred(Augmentor.randparam(ShearY(1:10), rect)) in 1:10
        @test_throws MethodError Augmentor.toaffinemap(ShearY(45), rect)
        @test @inferred(Augmentor.toaffinemap(ShearY(45), rect, 45)) ≈ AffineMap([1. -1.; 0. 1.], [2.,0.])
        @test @inferred(Augmentor.toaffinemap(ShearY(-60:60), rect, 45)) ≈ AffineMap([1. -1.; 0. 1.], [2.,0.])
        @test @inferred(Augmentor.toaffinemap(ShearY(-45), rect, -45)) ≈ AffineMap([1. 1.; 0. 1.], [-2.,0.])
        @test @inferred(Augmentor.toaffinemap(ShearY(-60:1:60), rect, -45)) ≈ AffineMap([1. 1.; 0. 1.], [-2.,0.])
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(ShearY(10), nothing)
        @test Augmentor.supports_eager(ShearY) === false
        # TODO: actual content tests (maybe test_reference)
        img_out1 = @inferred Augmentor.applyeager(ShearY(45), square)
        img_out2 = @inferred Augmentor.applyeager(ShearY(-45), square)
        @test axes(img_out1) == (0:4, 1:3)
        @test axes(img_out1) == axes(img_out2)
        imgs = [
            (square),
            (view(square, :, :)),
            (Augmentor.prepareaffine(square)),
            (OffsetArray(square, -1, -1)),
            (view(square, IdentityRange(1:3), IdentityRange(1:3))),
        ]
        @testset "fixed parameter" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(ShearY(45), img_in))
                @test parent(res) == parent(img_out1)
                @test typeof(res) == typeof(img_out1)
                res = @inferred(Augmentor.applyeager(ShearY(-45), img_in))
                @test parent(res) == parent(img_out2)
                @test typeof(res) == typeof(img_out2)
                # test same with tuple of images
                res1, res2 = @inferred(Augmentor.applyeager(ShearY(45), (img_in, N0f8.(img_in))))
                @test parent(res1) == parent(img_out1)
                @test parent(res2) == parent(img_out1)
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: OffsetArray{N0f8}
                res1, res2 = @inferred(Augmentor.applyeager(ShearY(-45), (img_in, N0f8.(img_in))))
                @test parent(res1) == parent(img_out2)
                @test parent(res2) == parent(img_out2)
                @test typeof(res1) == typeof(img_out2)
                @test typeof(res2) <: OffsetArray{N0f8}
            end
            # check that the affine map is computed for each image
            res1, res2 = @inferred(Augmentor.applyeager(ShearY(45), (square, OffsetArray(square,-5,-5))))
            @test collect(res1) == collect(res2)
            @test axes(res1) != axes(res2)
        end
        @testset "random parameter" begin
            for img_in in imgs
                res1, res2 = @inferred(Augmentor.applyeager(ShearY(1:60), (img_in, N0f8.(img_in))))
                # make sure same angle is used
                @test res1 == res2
                @test typeof(res1) == typeof(img_out1)
                @test typeof(res2) <: OffsetArray{N0f8}
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(ShearY) === true
        @test_throws MethodError Augmentor.applyaffine(ShearY(45), nothing)
        wv = @inferred Augmentor.applyaffine(ShearY(45), Augmentor.prepareaffine(square))
        @test parent(wv).itp.coefs === square
        @test axes(wv) == (0:4,1:3)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        wv2 = @inferred Augmentor.applyaffine(ShearY(-45), wv)
        @test parent(wv).itp.coefs === square
        @test wv2[1:3,1:3] == square
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(ShearY) === true
        @test_throws MethodError Augmentor.applyaffineview(ShearY(45), nothing)
        wv = @inferred Augmentor.applyaffineview(ShearY(45), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test axes(wv) == (0:4,1:3)
        wv2 = @inferred Augmentor.applyaffineview(ShearY(-45), wv)
        @test typeof(wv2) <: SubArray{eltype(square),2}
        @test typeof(parent(wv2)) <: InvWarpedView
        @test parent(parent(wv2)).itp.coefs === square
        @test wv2[1:3,1:3] == square
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(ShearY(45)) === true
        wv = @inferred Augmentor.applylazy(ShearY(45), square)
        @test parent(wv).itp.coefs === square
        @test axes(wv) == (0:4,1:3)
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
        wv2 = @inferred Augmentor.applylazy(ShearY(-45), wv)
        @test parent(wv).itp.coefs === square
        @test wv2[1:3,1:3] == square
        @test typeof(wv) <: InvWarpedView{eltype(square),2}
    end
    @testset "view" begin
        @test Augmentor.supports_view(ShearY) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(ShearY) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(ShearY) === false
    end
end


text:
function build_simple_lp!(optim)
    MOI.empty!(optim)
    @test MOI.is_empty(optim)

    # add 10 variables - only diagonal is relevant
    X = MOI.add_variables(optim, 2)

    # add sdp constraints - only ensuring positivenesse of the diagonal
    vov = MOI.VectorOfVariables(X)

    c1 = MOI.add_constraint(optim,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X[1]),
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    c2 = MOI.add_constraint(optim,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1]),
            MOI.ScalarAffineTerm(2.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    b1 = MOI.add_constraint(optim,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1])
        ], 0.0), MOI.GreaterThan(0.0))

    b2 = MOI.add_constraint(optim,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.GreaterThan(0.0))

    MOI.set(optim,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
        )
    MOI.set(optim, MOI.ObjectiveSense(), MOI.MIN_SENSE)
end

const optimizer_maxiter = MOI.instantiate(
    ()->ProxSDP.Optimizer(max_iter = 1, log_verbose = false),
    with_bridge_type = Float64)
const optimizer_timelimit = MOI.instantiate(
    ()->ProxSDP.Optimizer(
        tol_gap = 1e-16, tol_feasibility = 1e-16,
        tol_primal = 1e-16, tol_dual = 1e-16,
        tol_feasibility_dual = 1e-16, tol_psd = 1e-16,
        time_limit = 0.0, check_dual_feas = false,
        check_dual_feas_freq = 1),
    with_bridge_type = Float64)

@testset "MOI status" begin
    @testset "MOI.OPTIMAL" begin
        build_simple_lp!(optimizer_bridged)
        MOI.optimize!(optimizer_bridged)
        @test MOI.get(optimizer_bridged, MOI.TerminationStatus()) == MOI.OPTIMAL
        MOI.empty!(optimizer_bridged)
    end

    @testset "MOI.ITERATION_LIMIT" begin
        build_simple_lp!(optimizer_maxiter)
        MOI.optimize!(optimizer_maxiter)
        @test MOI.get(optimizer_maxiter, MOI.TerminationStatus()) == MOI.ITERATION_LIMIT
        MOI.empty!(optimizer_maxiter)
    end

    @testset "MOI.TIME_LIMIT" begin
        build_simple_lp!(optimizer_timelimit)
        MOI.optimize!(optimizer_timelimit)
        @test MOI.get(optimizer_timelimit, MOI.TerminationStatus()) == MOI.TIME_LIMIT
        MOI.empty!(optimizer_timelimit)
    end
end


text:
const res_desc_key = "desc";
const res_desc_comment_key = "comment";;
const res_desc_log_key = "log";

struct ResDesc
    desc_dict::Dict
    function ResDesc()
        desc_dict = Dict();
        desc_dict[res_desc_comment_key] = "";
        desc_dict[res_desc_log_key] = "";
        return new(desc_dict);
    end
    function ResDesc(desc_dict::Dict)
        validate_res_desc_dict(desc_dict);
        return new(desc_dict)
    end
end

get_desc_log(resdesc::ResDesc) = resdesc.desc_dict[res_desc_log_key];
get_desc_comment(resdesc::ResDesc) = resdesc.desc_dict[res_desc_comment_key];

add_desc_log_line!(resdesc::ResDesc, log_line::AbstractString) =
    resdesc.desc_dict[res_desc_log_key] *= "$(now()): $(log_line)\n";

set_desc_comment!(resdesc::ResDesc, new_comment::AbstractString) =
    resdesc.desc_dict[res_desc_comment_key] = new_comment;

load_res_desc(repo::ResRepoTree) =
    FileIO.load(repo.desc_file)[res_desc_key];

save_res_desc!(repo::ResRepoTree, resdesc::ResDesc) =
    FileIO.save(repo.desc_file, res_desc_key, resdesc);


text:
module CommunalHelperSoundAreaTrigger

using ..Ahorn, Maple

@mapdef Trigger "CommunalHelper/SoundAreaTrigger" SoundAreaTrigger(
    x::Integer,
    y::Integer,
    width::Integer = Maple.defaultTriggerWidth,
    height::Integer = Maple.defaultTriggerHeight,
    nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[],
    event::String="",
)

const placements = Ahorn.PlacementDict(
    "Sound Area Trigger (Communal Helper)" => Ahorn.EntityPlacement(
        SoundAreaTrigger,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            x, y = Ahorn.position(entity)
            width = Int(get(entity.data, "width", 8))
            height = Int(get(entity.data, "height", 8))
            entity.data["nodes"] = [(x + width / 2, y + height / 2)]
        end
    ),
)

Ahorn.nodeLimits(entity::SoundAreaTrigger) = 1, 1

end


text:
function createmesh!(bhadata::Dict{Symbol, Any})

  seg = bhadata[:segs]
  bhadata[:materials] = materials
    # Future: Use updatematerialtable() to add or change an entry in the materials table
  bhadata[:media] = media
  # Future: Use updatemediatable() to add or change an entry in the media table
  bhadata[:medium] = :mud
  bhadata[:noofelements] = Int(sum([bha[:segs][i][3] for i in 1:length(bha[:segs])]))
  bhadata[:noofnodes] = bha[:noofelements] + 1

  elements = [BHAPtfem.Element() for i in 1:bhadata[:noofelements]]
  nodes = [BHAPtfem.Node() for i in 1:bhadata[:noofnodes]]
  properties = fill(0.0, length(seg), 3)

  el = 0                                            # Element index
  propertiesindex = 0                      # Row in properties matrix

  sgmedium = bhadata[:media][bhadata[:medium]] / 231.     # in [in^3]

  for i in 1:length(seg)
    mat = bhadata[:materials][seg[i][2]]
    radius = seg[i][5] / 2.0
    if seg[i][1] in [:bit, :stabilizer]
      if (el+1) < bhadata[:noofelements] && radius > nodes[el+1].bharadius
        nodes[el+1].fc = seg[i][6]
        nodes[el+1].bharadius = radius
      end
    else
      propertiesindex += 1
      nels = Int(seg[i][3])       # Divide bha in 1' elements
      elength = 12.0      # [in]

      # Cross sectional area: csa = pi * d^2 / 4

      csa = pi * (seg[i][5]^2 - seg[i][4]^2) / 4.0
      eweight = elength * csa * ( mat.sg - sgmedium)

      # See http://en.wikipedia.org/wiki/Area_moment_of_inertia
      # Area moment of inertia: Pi * r^4 / 4 or Pi * d^4 / 64

      # See http://en.wikipedia.org/wiki/Polar_moment_of_inertia
      # Polar moment of inertia: Pi * r^4 / 2 or Pi * d^4 / 32

      d4 = seg[i][5]^4 - seg[i][4]^4

      # properties = [1:nels ; ea er gj]
      properties[propertiesindex,:] = [csa * mat.ea pi * mat.er * d4 / 64 pi * mat.gj * d4 / 32]

      for j in 1:nels
      	el += 1

        # Set nodes data

		# Use largest bharadius if 2 elements differ. Could also have been influenced by stabilizers above.

		nodes[el].bharadius = max(radius, nodes[el].bharadius)

        # Preset value for next element

		nodes[el+1].fc = seg[i][6]
		nodes[el+1].bharadius = radius

		# Pattern always the same, asign 1/2 to this node and 1/2 to next node,
		# during next loop add 1/2 of next element.

		nodes[el].weight += eweight / 2.0
		nodes[el+1].weight = eweight / 2.0

		# Set elements data

		elements[el].length = elength
		elements[el].id = seg[i][4]
		elements[el].od = seg[i][5]
        elements[el].npind = propertiesindex

      end
    end
  end

  properties = properties[1:propertiesindex, :]
  bha[:properties] = properties
  bha[:elements] = elements
  bha[:nodes] = nodes

  (properties, nodes, elements)
end