text:
"""
`Colorant{T,N}` is the abstract super-type of all types in ColorTypes,
and refers to both (opaque) colors and colors-with-transparency (alpha
channel) information.  `T` is the element type (extractable with
`eltype`) and `N` is the number of *meaningful* entries (extractable
with `length`), i.e., the number of arguments you would supply to the
constructor.
"""
abstract type Colorant{T,N} end

# Colors (without transparency)
"""
`Color{T,N}` is the abstract supertype for a color (or
grayscale) with no transparency.
"""
abstract type Color{T, N} <: Colorant{T,N} end

"""
`AbstractGray{T}` is an abstract supertype for gray types which corresponds to
real numbers, where `0` means black and `1` means white.

!!! compat "ColorTypes v0.12"
    Prior to ColorTypes v0.12, `AbstractGray{T}` was the alias for `Color{T,1}`;
    after ColorTypes v0.12, `Color{T,1}` does not necessarily correspond to the
    black-gray-white colors.
"""
abstract type AbstractGray{T} <: Color{T,1} end

"""
`AbstractRGB{T}` is an abstract supertype for red/green/blue color types that
can be constructed as `C(r, g, b)` and for which the elements can be
extracted as `red(c)`, `green(c)`, `blue(c)`. You should *not* make
assumptions about internal storage order, the number of fields, or the
representation. One `AbstractRGB` color-type, `RGB24`, is not
parametric and does not have fields named `r`, `g`, `b`.
"""
abstract type AbstractRGB{T} <: Color{T,3} end

# Types with transparency
"""
`TransparentColor{C,T,N}` is the abstract type for any
color-with-transparency.  The `C` parameter refers to the type of the
pure color (without transparency) and can be extracted with
`color_type`. `T` is the element type of both `C` and the `alpha`
channel, and `N` has the same meaning as in `Colorant` (and is 1 larger
than the corresponding color type).

All transparent types should support two modes of construction:

    P(color, alpha)
    P(component1, component2, component3, alpha) (assuming a 3-component color)

For a `Colorant` `c`, the color component can be extracted with
`color(c)`, and the alpha channel with `alpha(c)`. Note that types
such as `ARGB32` do not have a field named `alpha`.

Most concrete types, like `RGB`, have both `ARGB` and `RGBA`
transparent analogs.  These two indicate different internal storage
order (see `AlphaColor` and `ColorAlpha`, and the `alphacolor` and
`coloralpha` functions).
"""
abstract type TransparentColor{C<:Color,T,N} <: Colorant{T,N} end

"""
`AlphaColor` is an abstract supertype for types like `ARGB`, where the
alpha channel comes first in the internal storage order. **Note** that
the constructor order is still `(color, alpha)`.
"""
abstract type AlphaColor{C<:Color,T,N} <: TransparentColor{C,T,N} end

"""
`ColorAlpha` is an abstract supertype for types like `RGBA`, where the
alpha channel comes last in the internal storage order.
"""
abstract type ColorAlpha{C<:Color,T,N} <: TransparentColor{C,T,N} end

# These are types we'll dispatch on.
Color3{T}                          = Color{T,3}
TransparentGray{C<:AbstractGray,T} = TransparentColor{C,T,2}
Transparent3{C<:Color3,T}          = TransparentColor{C,T,4}
TransparentRGB{C<:AbstractRGB,T}   = TransparentColor{C,T,4}
ColorantNormed{T<:Normed,N}        = Colorant{T,N}
AbstractAGray{C<:AbstractGray,T}   = AlphaColor{C,T,2}
AbstractGrayA{C<:AbstractGray,T}   = ColorAlpha{C,T,2}
AbstractARGB{C<:AbstractRGB,T}     = AlphaColor{C,T,4}
AbstractRGBA{C<:AbstractRGB,T}     = ColorAlpha{C,T,4}

# With reordered type parameters (also useful for dispatch)
ColorantN{N,T}                  = Colorant{T,N}
ColorN{N,T}                     = Color{T,N}
TransparentColorN{N,C<:Color,T} = TransparentColor{C,T,N}
AlphaColorN{N,C<:Color,T}       = AlphaColor{C,T,N}
ColorAlphaN{N,C<:Color,T}       = ColorAlpha{C,T,N}

"""
`RGB` is the standard Red-Green-Blue (sRGB) colorspace.  Values of the
individual color channels range from 0 (black) to 1 (saturated). If
you want "Integer" storage types (e.g., 255 for full color), use `N0f8(1)`
instead (see FixedPointNumbers).
"""
struct RGB{T<:Fractional} <: AbstractRGB{T}
    r::T # Red [0,1]
    g::T # Green [0,1]
    b::T # Blue [0,1]

    RGB{T}(r::T, g::T, b::T) where {T <: Fractional} = new{T}(r, g, b)
end

"""
`BGR` is a variant of `RGB` with the opposite storage order.  Note
that the constructor is still called in the order `BGR(r, g, b)`.
This storage order is noteworthy because on little-endian machines,
`BGRA` (with transparency) can be `reinterpret`ed to the `UInt32`
color format used by libraries such as Cairo and OpenGL.
"""
struct BGR{T<:Fractional} <: AbstractRGB{T}
    b::T
    g::T
    r::T

    BGR{T}(r::T, g::T, b::T) where {T <: Fractional} = new{T}(b, g, r)
end

"""
`XRGB` is a variant of `RGB` which has a padding element inserted at
the beginning. In some applications it may have useful
memory-alignment properties.

Like all other AbstractRGB objects, the constructor is still called
`XRGB(r, g, b)`.
"""
struct XRGB{T<:Fractional} <: AbstractRGB{T}
    alphadummy::T
    r::T
    g::T
    b::T

    XRGB{T}(r::T, g::T, b::T) where {T <: Fractional} = new{T}(oneunit(T), r, g, b)
end

"""
`RGBX` is a variant of `RGB` which has a padding element inserted at
the end. In some applications it may have useful
memory-alignment properties.

Like all other AbstractRGB objects, the constructor is still called
`RGBX(r, g, b)`.
"""
struct RGBX{T<:Fractional} <: AbstractRGB{T}
    r::T
    g::T
    b::T
    alphadummy::T

    RGBX{T}(r::T, g::T, b::T) where {T <: Fractional} = new{T}(r, g, b, oneunit(T))
end

"`HSV` is the Hue-Saturation-Value colorspace."
struct HSV{T<:AbstractFloat} <: Color{T,3}
    h::T # Hue in [0,360]
    s::T # Saturation in [0,1]
    v::T # Value in [0,1]
end

"`HSB` (Hue-Saturation-Brightness) is an alias for `HSV`."
const HSB = HSV

"`HSL` is the Hue-Saturation-Lightness colorspace."
struct HSL{T<:AbstractFloat} <: Color{T,3}
    h::T # Hue in [0,360]
    s::T # Saturation in [0,1]
    l::T # Lightness in [0,1]
end

"`HSI` is the Hue-Saturation-Intensity colorspace."
struct HSI{T<:AbstractFloat} <: Color{T,3}
    h::T # Hue in [0,360]
    s::T # Saturation in [0,1]
    i::T # Intensity in [0,1]
end

"""
`XYZ` is the CIE 1931 XYZ colorspace. It is a linear colorspace,
meaning that mathematical operations such as addition, subtraction,
and scaling make "colorimetric sense" in this colorspace.
"""
struct XYZ{T<:AbstractFloat} <: Color{T,3}
    x::T
    y::T
    z::T
end

"`xyY` is the CIE 1931 xyY (chromaticity + luminance) space"
struct xyY{T<:AbstractFloat} <: Color{T,3}
    x::T
    y::T
    Y::T
end

"`Lab` is the CIELAB colorspace."
struct Lab{T<:AbstractFloat} <: Color{T,3}
    l::T # Lightness in [0,100]
    a::T # Red/Green
    b::T # Blue/Yellow
end

"`LCHab` is the Luminance-Chroma-Hue, Polar-Lab colorspace"
struct LCHab{T<:AbstractFloat} <: Color{T,3}
    l::T # Lightness in [0,100]
    c::T # Chroma
    h::T # Hue in [0,360]
end

"`Luv` is the CIELUV colorspace"
struct Luv{T<:AbstractFloat} <: Color{T,3}
    l::T # Lightness in [0,100]
    u::T # Red/Green
    v::T # Blue/Yellow
end

"`LCHuv` is the Luminance-Chroma-Hue, Polar-Luv colorspace"
struct LCHuv{T<:AbstractFloat} <: Color{T,3}
    l::T # Lightness in [0,100]
    c::T # Chroma
    h::T # Hue in [0,360]
end

"`Oklab` is the Oklab colorspace."
struct Oklab{T<:AbstractFloat} <: Color{T,3}
    l::T # Lightness in [0,100]
    a::T # Red/Green
    b::T # Blue/Yellow
end

"`LCHOklab` is the Luminance-Chroma-Hue, Polar-Oklab colorspace"
struct LCHOklab{T<:AbstractFloat} <: Color{T,3}
    l::T # Lightness in [0,100]
    c::T # Chroma
    h::T # Hue in [0,360]
end

"`DIN99` is the (L99, a99, b99) adaptation of CIELAB"
struct DIN99{T<:AbstractFloat} <: Color{T,3}
    l::T # L99
    a::T # a99
    b::T # b99
end

"`DIN99d` is the (L99d, a99d, b99d) improvement on DIN99"
struct DIN99d{T<:AbstractFloat} <: Color{T,3}
    l::T # L99d
    a::T # a99d
    b::T # b99d
end

"`DIN99o` is the (L99o, a99o, b99o) adaptation of CIELAB"
struct DIN99o{T<:AbstractFloat} <: Color{T,3}
    l::T # L99o
    a::T # a99o
    b::T # b99o
end

"""
`LMS` is the Long-Medium-Short colorspace based on activation of the
three cone photoreceptors.  Like `XYZ`, this is a linear color space.
"""
struct LMS{T<:AbstractFloat} <: Color{T,3}
    l::T # Long
    m::T # Medium
    s::T # Short
end

"`YIQ` is a color encoding, for example used in NTSC transmission."
struct YIQ{T<:AbstractFloat} <: Color{T,3}
    y::T
    i::T
    q::T
end

"`YCbCr` is the Y'CbCr color encoding often used in digital photography or video"
struct YCbCr{T<:AbstractFloat} <: Color{T,3}
    y::T
    cb::T
    cr::T
end

"""
`RGB24` uses a `UInt32` representation of color, 0xAARRGGBB, where
R=red, G=green, B=blue and A is irrelevant. This format is often used
by external libraries such as Cairo.

`RGB24` colors do not have fields named `r`, `g`, `b`, but you can
still extract the individual components with `red(c)`, `green(c)`,
`blue(c)`.  You can construct them directly from a `UInt32`, or as
`RGB(r, g, b)`.
"""
struct RGB24 <: AbstractRGB{N0f8}
    color::UInt32

    RGB24(col::UInt32, ::Type{Val{true}}) = new(col)
end
# The Val{true} constructor solves a consistency problem, one that's
# most obvious for Gray24: Gray24(x) interprets x as a value between 0
# and 1; since a UInt32 can have value 0 or 1 (among other
# possibilities), having Gray24(x::UInt32) mean something entirely
# different (interpreting the UInt32 as a bit pattern rather than a
# value) would be inconsistent, confusing, and error-prone. Because we
# can construct an RGB24 from a grayscale value (just repeating the
# value for all 3 color channels), the same problem extends to all the
# "24/32" color types.

# However, some constructors "build" RGB24s as a UInt32 bit pattern,
# so we need some way of turning a UInt32 directly into an RGB24.
# This is the role of the Val{true} constructor. This is an internal
# implementation detail, and user code should not use it directly:
# use reinterpret(RGB24, x::UInt32) instead. (The Val{true}
# constructor is used to implement reinterpret, see traits.jl.)
function RGB24(r::N0f8, g::N0f8, b::N0f8)
    r32 = UInt32(reinterpret(r)) << 16
    g32 = UInt32(reinterpret(g)) << 8
    b32 = UInt32(reinterpret(b))
    reinterpret(RGB24, r32 | g32 | b32)
end

"""
`ARGB32` uses a `UInt32` representation of color, 0xAARRGGBB, where
R=red, G=green, B=blue and A is the alpha channel. This format is
often used by external libraries such as Cairo.  On a little-endian
machine, this type has the exact same storage format as `BGRA{N0f8}`.

`ARGB32` colors do not have fields named `alpha`, `r`, `g`, `b`, but
you can still extract the individual components with `alpha(c)`,
`red(c)`, `green(c)`, `blue(c)`.  You can construct them directly from
a `UInt32`, or as `ARGB32(r, g, b, alpha)`.
"""
struct ARGB32 <: AbstractARGB{RGB24, N0f8}
    color::UInt32

    ARGB32(col::UInt32, ::Type{Val{true}}) = new(col)
end
function ARGB32(r::N0f8, g::N0f8, b::N0f8, alpha::N0f8 = N0f8(1))
    a32 = UInt32(reinterpret(alpha)) << 24
    r32 = UInt32(reinterpret(r)) << 16
    g32 = UInt32(reinterpret(g)) << 8
    b32 = UInt32(reinterpret(b))
    reinterpret(ARGB32, a32 | r32 | g32 | b32)
end
function ARGB32(c::RGB24, alpha::Real = N0f8(1))
    checkval(ARGB32, alpha)
    a32 = UInt32(reinterpret(_rem(alpha, N0f8))) << 24
    reinterpret(ARGB32, a32 | (c.color & 0xFFFFFF))
end

"""
`Gray` is a grayscale object. You can extract its value with `gray(c)`.
"""
struct Gray{T<:Union{Fractional,Bool}} <: AbstractGray{T}
    val::T
end

"""
`Gray24` uses a `UInt32` representation of color, 0xAAIIIIII, where
I=intensity (grayscale value) and A is irrelevant. Each II pair is
assumed to be the same.  This format is often used by external
libraries such as Cairo.

You can extract the single gray value with `gray(c)`.  You can
construct them directly from a `UInt32`, or as `Gray24(i)`. Note that
`i` is interpreted on a scale from 0 (black) to 1 (white).
"""
struct Gray24 <: AbstractGray{N0f8}
    color::UInt32

    Gray24(c::UInt32, ::Type{Val{true}}) = new(c)
end
Gray24(val::N0f8) = reinterpret(Gray24, reinterpret(val) * 0x010101)

"""
`AGray32` uses a `UInt32` representation of color, 0xAAIIIIII, where
I=intensity (grayscale value) and A=alpha. Each II pair is
assumed to be the same.  This format is often used by external
libraries such as Cairo.

You can extract the single gray value with `gray(c)` and the alpha as
`alpha(c)`.  You can construct them directly from a `UInt32`, or as
`AGray32(i,alpha)`. Note that `i` and `alpha` are interpreted on a
scale from 0 (black) to 1 (white).
"""
struct AGray32 <: AbstractAGray{Gray24, N0f8}
    color::UInt32

    AGray32(c::UInt32, ::Type{Val{true}}) = new(c)
end
function AGray32(val::N0f8, alpha::N0f8 = N0f8(1))
    a32 = UInt32(reinterpret(alpha)) << 24
    reinterpret(AGray32, a32 | reinterpret(val) * 0x010101)
end
function AGray32(g::Gray24, alpha::Real= N0f8(1))
    checkval(AGray32, alpha)
    a32 = UInt32(reinterpret(_rem(alpha, N0f8))) << 24
    reinterpret(AGray32, a32 | (g.color & 0xffffff))
end

# Generated code:
#   - TransparentColor definitions (e.g., ARGB and RGBA) with inner constructors
#   - `export`s
#   - `coloralpha(::Color)` and `alphacolor(::Color)` traits for corresponding types

# Note: with the exceptions of `alphacolor` and `coloralpha`, all
# traits in the rest of this file are intended just for internal use

# Provide the field names in the order expected by the constructor
colorfields(::Type{C}) where {C<:Color} = (fieldnames(C)...,)
colorfields(::Type{C}) where {C<:XRGB} = (:r, :g, :b)
colorfields(::Type{C}) where {C<:RGBX} = (:r, :g, :b)
colorfields(::Type{C}) where {C<:BGR } = (:r, :g, :b)
colorfields(::Type{P}) where {P<:TransparentColor} = tuple(colorfields(color_type(P))..., :alpha)
colorfields(c::Colorant) = colorfields(typeof(c))

eltype_default(::Type) = Float32
eltype_default(::Type{C}) where {C<:AbstractRGB } = N0f8
eltype_default(::Type{C}) where {C<:AbstractGray} = N0f8
eltype_default(::Type{C}) where {C<:TransparentColor} = eltype_default(color_type(C))

# TODO: Generalize the promotion rules so that they are reasonable even for custom types.
@inline function promote_args_type(::Type{C}, args...) where C <: Color
    T = promote_type(map(typeof, args)...)
    _promote_args_type(eltype_default(C), T)
end
@inline function promote_args_type(::Type{C}, args...) where {C <: TransparentColor}
    Ta = typeof(last(args)) # alpha
    T = _promote_wol(map(typeof, args)...)
    if T <: Union{Integer, FixedPoint} && Ta <: Integer
        _promote_args_type(eltype_default(C), T)
    else
        _promote_args_type(eltype_default(C), promote_type(T, Ta))
    end
end
# a variant of `promote_type` that ignores the last (i.e. alpha) element
@inline _promote_wol(t, tail...) = length(tail) == 1 ? t : promote_type(t, _promote_wol(tail...))

_promote_args_type(::Type{Tdef}, ::Type{T}) where {Tdef<:AbstractFloat, T<:AbstractFloat} = T
_promote_args_type(::Type{Tdef}, ::Type{T}) where {Tdef<:AbstractFloat, T<:Real} = Tdef
_promote_args_type(::Type{Tdef}, ::Type{T}) where {Tdef<:FixedPoint, T<:Fractional} = T
_promote_args_type(::Type{Tdef}, ::Type{T}) where {Tdef<:FixedPoint, T<:Real} = Tdef
_promote_args_type(::Type{Tdef}, ::Type{T}) where {Tdef, T} = promote_type(Tdef, T)

# Generate transparent versions
macro make_alpha(C, acol, cola, fields, constrfields, ub, elty)
    # ub = upper-bound on T in C{T}
    # elty = default element type when supplied with Integer arguments
    fields = fields.args
    constrfields = constrfields.args
    N = length(fields) + 1
    Tfields = Expr[:($f::T) for f in fields]
    Targs   = Expr[:($f::T) for f in constrfields]
    Tconstr = Expr(:<:, :T, ub)
    esc(quote
        struct $acol{$Tconstr} <: AlphaColor{$C{T}, T, $N}
            alpha::T
            $(Tfields...)

            $acol{T}($(Targs...), alpha::T=oneunit(T)) where {$Tconstr} = new{T}(alpha, $(fields...))
        end
        struct $cola{$Tconstr} <: ColorAlpha{$C{T}, T, $N}
            $(Tfields...)
            alpha::T

            $cola{T}($(Targs...), alpha::T=oneunit(T)) where {$Tconstr} = new{T}($(fields...), alpha)
        end
        export $acol, $cola
        alphacolor(::Type{C}) where {C<:$C} = $acol
        coloralpha(::Type{C}) where {C<:$C} = $cola
    end)
end

const color3types = map(s->getfield(ColorTypes,s),
                        filter(names(ColorTypes, all=false)) do s
                            isdefined(ColorTypes, s) || return false
                            t = getfield(ColorTypes, s)
                            isa(t, Type{<:Color3}) && !isabstracttype(t)
                        end |> unique
                        )

# The above should have filtered out every non-DataType that's not also a
# wrapped UnionAll-wrapped DataType.
const parametric3 = filter(x->!isa(x, DataType) || !isempty(x.parameters), color3types)

# For searchability, explicitly enumerate the type names. (issue #16)
for (C, acol, cola) in [(DIN99d, :ADIN99d, :DIN99dA),
                        (DIN99o, :ADIN99o, :DIN99oA),
                        (DIN99, :ADIN99, :DIN99A),
                        (HSI, :AHSI, :HSIA),
                        (HSL, :AHSL, :HSLA),
                        (HSV, :AHSV, :HSVA),
                        (LCHab, :ALCHab, :LCHabA),
                        (LCHuv, :ALCHuv, :LCHuvA),
                        (LCHOklab, :ALCHOklab, :LCHOklabA),
                        (LMS, :ALMS, :LMSA),
                        (Lab, :ALab, :LabA),
                        (Luv, :ALuv, :LuvA),
                        (Oklab, :AOklab, :OklabA),
                        (XYZ, :AXYZ, :XYZA),
                        (YCbCr, :AYCbCr, :YCbCrA),
                        (YIQ, :AYIQ, :YIQA),
                        (xyY, :AxyY, :xyYA),
                        (BGR, :ABGR, :BGRA),
                        (RGB, :ARGB, :RGBA),
                        (Gray, :AGray, :GrayA)]
    fn  = Expr(:tuple, fieldnames(C)...)
    cfn = Expr(:tuple, colorfields(C)...)
    elty = eltype_default(C)
    ub   = elty <: FixedPoint ? Fractional : AbstractFloat
    Csym = nameof(C)
    @eval @make_alpha $Csym $acol $cola $fn $cfn $ub $elty
end

const GrayLike = Union{Number, AbstractGray}

_real(x) = real(x)
@noinline function _real(@nospecialize(x::Colorant))
    throw(ArgumentError("""
        Color objects other than `AbstractGray` cannot be used as if they were `Real` arguments."""))
end
_real(x::AbstractGray) = real(x)

function (::Type{C})() where {N, C <: ColorantN{N}}
    d0 = zero(eltype_default(C))
    dx = C <: TransparentColor ? oneunit(eltype_default(C)) : d0
    C(ntuple(_ -> d0, Val(N - 1))..., dx)
end

(::Type{C})(x         ) where {C <: Color    } = _new_colorant(C, x)
(::Type{C})(x, y      ) where {C <: ColorN{2}} = _new_colorant(C, x, y)
(::Type{C})(x, y, z   ) where {C <: ColorN{3}} = _new_colorant(C, x, y, z)
(::Type{C})(x, y, z, w) where {C <: ColorN{4}} = _new_colorant(C, x, y, z, w)

(::Type{C})(x                  ) where {C <: TransparentColor    } = _new_colorant(C, x)
(::Type{C})(x,          alpha  ) where {C <: TransparentColor    } = _new_colorant(C, x, alpha)
(::Type{C})(x, y,       alpha=1) where {C <: TransparentColorN{3}} = _new_colorant(C, x, y, alpha)
(::Type{C})(x, y, z,    alpha=1) where {C <: TransparentColorN{4}} = _new_colorant(C, x, y, z, alpha)
(::Type{C})(x, y, z, w, alpha=1) where {C <: TransparentColorN{5}} = _new_colorant(C, x, y, z, w, alpha)

(::Type{C})(g::GrayLike) where {C <: TransparentColorN{2}} = _new_colorant(C, g, 1)
(::Type{C})(g::GrayLike) where {C <: AbstractRGB} = _new_colorant(C, g, g, g)
(::Type{C})(g::GrayLike) where {C <: TransparentRGB} = _new_colorant(C, g, g, g, oneunit(g))
(::Type{C})(g::GrayLike, alpha) where {C <: TransparentRGB} = _new_colorant(C, g, g, g, alpha)

function _new_colorant(::Type{C}, args::Vararg{Any,N}) where {N, C <: ColorantN{N}}
    rargs = _real.(args)
    base_colorant_type(C){promote_args_type(C, rargs...)}(rargs...)
end

function _new_colorant(::Type{C}, c::Colorant) where {N, C <: ColorantN{N}}
    convert(C, c)
end

function _new_colorant(::Type{C}, c::Colorant, alpha) where {N, C <: TransparentColorN{N}}
    # TODO: Make the following a single `convert` call
    cc = convert(base_color_type(C), color(c))
    convert(C, cc, alpha)
end

function (::Type{C})(x::UInt32) where C <: Union{RGB24, Gray24, ARGB32, AGray32}
    x <= UInt32(1) || throw_bit_pattern_error(C, x)
    reinterpret(C, C <: Color ? (-x) & 0xffffff : (-x) | 0xff000000)
end
function (::Type{C})(x::GrayLike) where C <: Union{RGB24, Gray24, ARGB32, AGray32}
    checkval(C, real(x))
    v = _rem(real(x), N0f8)
    return C <: Union{RGB24, ARGB32} ? C(v, v, v) : C(v)
end

# T might be a Normed, and so some inputs will result in an error.
# Try to make it a nice error.
function _new_colorant(::Type{C}, args::Vararg{GrayLike,N}) where {N, T, C <: Color{T,N}}
    isconcretetype(C) || throw(MethodError(C, args))
    r = real.(args)
    checkval(C, r...)
    C(_rem.(r, T)...)
end
function _new_colorant(::Type{TC}, args::Vararg{GrayLike,N}) where {N, T, C, TC <: TransparentColor{C,T,N}}
    r = real.(args)
    checkval(TC, r...)
    TC(_rem.(r, T)...)
end
function _new_colorant(::Type{TC}, args::Vararg{GrayLike,2}) where {T, C, TC <: TransparentColor{C,T,2}}
    r = real.(args)
    checkval(TC, r...)
    TC(_rem.(r, T)...)
end

"""
    alphacolor(::Type{<:Colorant})
    alphacolor(::Colorant)

Return the corresponding transparent color type/instance with storage order
(alpha, color).

# Examples
```jldocest; setup = :(using ColorTypes)
julia> alphacolor(RGB)
ARGB

julia> alphacolor(RGBA{Float32})
ARGB

julia> alphacolor(Gray(0.8)) === AGray(0.8, 1.0)
true
```
"""
alphacolor(::Type{C}) where {C<:AlphaColor} = base_colorant_type(C)
alphacolor(::Type{C}) where {C<:ColorAlpha} = alphacolor(base_color_type(C))

"""
    coloralpha(::Type{<:Colorant})
    coloralpha(::Colorant)

Return the corresponding transparent color type/instance with storage order
(color, alpha).

# Examples
```jldocest; setup = :(using ColorTypes)
julia> coloralpha(RGB)
RGBA

julia> coloralpha(ARGB{Float32})
RGBA

julia> coloralpha(Gray(0.8)) === GrayA(0.8, 1.0)
true
```
"""
coloralpha(::Type{C}) where {C<:ColorAlpha} = base_colorant_type(C)
coloralpha(::Type{C}) where {C<:AlphaColor} = coloralpha(base_color_type(C))

alphacolor(::Type{C}) where {C<:Gray24} = AGray32
alphacolor(::Type{C}) where {C<:RGB24} = ARGB32
alphacolor(::Type{C}) where {C<:XRGB} = ARGB
alphacolor(::Type{C}) where {C<:RGBX} = ARGB
coloralpha(::Type{C}) where {C<:XRGB} = RGBA
coloralpha(::Type{C}) where {C<:RGBX} = RGBA

### Validating the inputs for Normed constructors

# Throw helpful errors in case of trouble. The inlining here is
# carefully designed to reduce the impact on runtime performance, and
# we also avoid splatting.

@inline function isok(::Type{T}, x) where T<:Normed
    Δ = eps(T)/2 # as long as the number rounds to a valid number, that's OK
    (-Δ <= x) & (x < typemax(T)+Δ)
end
@inline function isok(::Type{T}, x, y) where T<:Normed
    Δ = eps(T)/2
    @fastmath n, m = minmax(x, y)
    (-Δ <= n) & (m < typemax(T)+Δ)
end

@inline function checkval(::Type{C}, a) where {T<:Normed, C<:Colorant{T}}
    isok(T, a) || throw_colorerror(C, (a,))
end
@inline function checkval(::Type{C}, a, b) where {T<:Normed, C<:Colorant{T}}
    isok(T, a, b) || throw_colorerror(C, (a, b))
end
@inline function checkval(::Type{C}, a, b, c) where {T<:Normed, C<:Colorant{T}}
    isok(T, a, b) & isok(T, b, c) || throw_colorerror(C, (a, b, c))
end
@inline function checkval(::Type{C}, a, b, c, d) where {T<:Normed, C<:Colorant{T}}
    Δ = eps(T)/2
    @fastmath min_ab, max_ab = minmax(a, b)
    @fastmath min_cd, max_cd = minmax(c, d)
    @fastmath n = min(min_ab, min_cd)
    @fastmath m = max(max_ab, max_cd)
    (-Δ <= n) & (m < typemax(T)+Δ) || throw_colorerror(C, (a, b, c, d))
end
@inline function checkval(::Type{C}, a, b, c, d, e) where {T<:Normed, C<:Colorant{T}}
    isok(T, a, b) & isok(T, c, d) & isok(T, e) || throw_colorerror(C, (a, b, c, d, e))
end

checkval(::Type{C}, args::Vararg{T}) where {T<:Normed, C<:Colorant{T}} = nothing
checkval(::Type{C}, args...) where {C<:Colorant} = nothing

function throw_colorerror_(::Type{T}, values) where T<:Normed
    Tmin = repr(typemin(T), context=:compact=>true)
    Tmax = repr(typemax(T), context=:compact=>true)
    bitstring = sizeof(T) == 1 ? "an 8-bit" : "a $(8*sizeof(T))-bit"
    throw(ArgumentError("""
component type $T is $bitstring type representing $(2^(8*sizeof(T))) values from $Tmin to $Tmax,
  but the values $values do not lie within this range.
  See the READMEs for FixedPointNumbers and ColorTypes for more information."""))
end

function throw_colorerror(::Type{C}, values::Tuple{Vararg{Integer}}) where C<:Colorant{N0f8}
    # Let's try to read the user's mind
    all(x->0<=x<=255, values) || throw_colorerror_(N0f8, values)
    if length(values) == 1
        vstr = "$(values[1]) is an integer"
    else
        vstr = "$values are integers"
    end
    Cstr = colorant_string_with_eltype(C)
    args = join(map(v -> "$v/255", values), ", ")
    throw(ArgumentError("""
The components of $(nameof(C)) are normalized to the range 0-1,
  but $vstr in the range 0-255.
  Consider dividing your input values by 255, for example: $Cstr($args)
  The component type N0f8 is an 8-bit type representing 256 values from 0 to 1.
  You can also use `reinterpret(N0f8, x % UInt8)` to encode the input into N0f8.
  See the READMEs for FixedPointNumbers and ColorTypes for more information."""))
end
function throw_colorerror(::Type{C}, values::Tuple) where {T<:Normed, C<:Colorant{T}}
    throw_colorerror_(T, values)
end

function throw_bit_pattern_error(@nospecialize(C), value)
    hex = string(value, base=16, pad=8)
    throw(ArgumentError("""
$C cannot be constructed or converted directly from a UInt32 input as a bit pattern.
  Use `reinterpret($C, 0x$hex)` instead."""))
end

_rem(x, ::Type{T}) where {T<:Normed} = x % T
_rem(x, ::Type{T}) where {T}         = convert(T, x)

struct ColorTypeResolutionError <: Exception
    func::Symbol
    msg::String
    C1
    C2
end

function Base.showerror(io::IO, ex::ColorTypeResolutionError)
    print(io, "in ", ex.func, ", ", ex.msg, ' ', ex.C1, " and ", ex.C2)
end


text:
using DependencyParser.Core

# Sentence: He wrote her a letter.
function build_sentence()
  words = ["He", "wrote", "her", "a", "letter", "."]
  tags = ["PRP", "VBN", "PRP", "DT", "NN", "JJ"]

  Sentence(zip(words, tags) |> collect)
end

ext: jl

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualMachineScaleSetExtensionsApi <: SwaggerApi
    client::Swagger.Client
end

"""

The operation to create or update an extension.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: vmssExtensionName::String (required)
Param: extensionParameters::VirtualMachineScaleSetExtension (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineScaleSetExtension
"""
function virtualMachineScaleSetExtensionsCreateOrUpdate(_api::VirtualMachineScaleSetExtensionsApi, resourceGroupName::String, vmScaleSetName::String, vmssExtensionName::String, extensionParameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", VirtualMachineScaleSetExtension, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/extensions/{vmssExtensionName}", ["azure_auth"], extensionParameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "vmssExtensionName", vmssExtensionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

The operation to delete the extension.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: vmssExtensionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: OperationStatusResponse
"""
function virtualMachineScaleSetExtensionsDelete(_api::VirtualMachineScaleSetExtensionsApi, resourceGroupName::String, vmScaleSetName::String, vmssExtensionName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", OperationStatusResponse, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/extensions/{vmssExtensionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "vmssExtensionName", vmssExtensionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

The operation to get the extension.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: vmssExtensionName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: VirtualMachineScaleSetExtension
"""
function virtualMachineScaleSetExtensionsGet(_api::VirtualMachineScaleSetExtensionsApi, resourceGroupName::String, vmScaleSetName::String, vmssExtensionName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineScaleSetExtension, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/extensions/{vmssExtensionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "vmssExtensionName", vmssExtensionName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a list of all extensions in a VM scale set.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineScaleSetExtensionListResult
"""
function virtualMachineScaleSetExtensionsList(_api::VirtualMachineScaleSetExtensionsApi, resourceGroupName::String, vmScaleSetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineScaleSetExtensionListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/extensions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualMachineScaleSetExtensionsCreateOrUpdate, virtualMachineScaleSetExtensionsDelete, virtualMachineScaleSetExtensionsGet, virtualMachineScaleSetExtensionsList


text:
#function f(x, y) println(x, " ", y) end
#f((println(1); 1), (println(2); 2))
#println(typeof(42))
#println(typeof(nothing))
#println(typeof(true))
#println(typeof(false))
#println(typeof(print()))
x = true
if x
	for i=3:5
		for j=(-9):3
			println(i, " ", j)
		end
	end
else
	for j=(-9):3
		println(j)
	end
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using OperatorFlux
using Test
using LinearAlgebra

@testset "ChebyshevTransform" begin
    # auxiliary data
    M = 32

    # test constructors
    trafo = ChebyshevTransform(modes = (3, 7))
    @test trafo.modes == (3, 7)

    # test forward & inverse
    for mode in [(4,), (4, 6), (4, 7, 9), (5,), (5, 6), (5, 7, 9)]
        z = rand(3, mode..., 5)
        trafo = ChebyshevTransform(modes = mode)
        @test norm(
            OperatorFlux.inverse(trafo, OperatorFlux.forward(trafo, z)) - z,
        ) ≤ 20 * eps(norm(z))
    end

    # test truncate_modes & pad_modes
    M = 32
    trafo = ChebyshevTransform(modes = (3, 7, 6))
    c = rand(7, M, M, M, 3)
    inds = [[1, 2, 3], [1, 2, 3, 4, 5, 6, 7], [1, 2, 3, 4, 5, 6]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = ChebyshevTransform(modes = (3, 6))
    c = rand(7, M, M, 3)
    inds = [[1, 2, 3], [1, 2, 3, 4, 5, 6]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = ChebyshevTransform(modes = (3,))
    c = rand(7, M, 3)
    inds = [[1, 2, 3]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = ChebyshevTransform(modes = (12, 5))
    c = rand(27, 16, 14, 1)
    inds = [collect(1:m) for m in size(c)[2:(end - 1)]]
    @test all(OperatorFlux.pad_modes(trafo, c, (32, 14))[:, inds..., :] .== c)

    # rountrip integration test
    for M in [32, 33]
        a, b = -π, π
        x = @. (cos(pi * (0:M) / M) + 1) / 2
        x = (b - a) .* x .+ a
        y = copy(x)
        z = x * y'
        z = z[:, :, :, :] # emulate channels and batches
        z = permutedims(z, (3, 1, 2, 4))
        trafo = ChebyshevTransform(modes = (12, 17))
        a = OperatorFlux.forward(trafo, z)
        b = OperatorFlux.truncate_modes(trafo, a)
        c = OperatorFlux.pad_modes(trafo, b, size(a)[2:(end - 1)])
        d = OperatorFlux.inverse(trafo, c)
        @test norm(d - z) ≤ eps(10 * norm(z))
    end

    # test Base extensions
    trafo = ChebyshevTransform(modes = (5, 11))
    @test ndims(trafo) == 2
    @test eltype(trafo) == Float32
    @test size(trafo) == [5, 11]
end


text:
#####################################################################################################
# Read BAM
#####################################################################################################
using CpelNano
using FASTX.FASTA

# IO
dataDir = "$(dirname(dirname(pathof(CpelNano))))/examples/wgbs/"
fasta = "$(dataDir)/reference/ref.fa"
bam = "$(dataDir)/bam/g1_s1.bam"
chr = "chrTest"

# Get FASTA data
fa_rec = CpelNano.get_chr_fa_rec(chr, fasta)

# Get size of chromosome
dna_seq = FASTA.sequence(fa_rec)
chr_size = length(dna_seq)

# Get position CpG sites in analysis region
diff = 151
chr_st = 1
chr_end = 1000 - diff
dna_seq = FASTA.sequence(String, fa_rec, chr_st:chr_end)
cpg_pos = map(x -> getfield(x, :offset), eachmatch(r"[Cc][Gg]", dna_seq)) .+ chr_st .- 1

# Read bam (get calls)
calls = CpelNano.get_calls_bam(bam, chr, chr_st, chr_end, cpg_pos, false, (0, 0, 0, 0))

#####################################################################################################
# Split BAM file
#####################################################################################################
using CpelNano

# IO
dataDir = "$(dirname(dirname(pathof(CpelNano))))/examples/wgbs/"
fasta = "$(dataDir)/reference/ref.fa"
bam = "$(dataDir)/bam/meth_sample.bam"
win_size = 500

# Split BAM
CpelNano.split_bam(bam,fasta,win_size)

#####################################################################################################
# ESTIMATION
#####################################################################################################

## Process chromosome
using CpelNano

# IO
dataDir = "$(dirname(dirname(pathof(CpelNano))))/examples/wgbs/"
fasta = "$(dataDir)/reference/ref.fa"
bam = "$(dataDir)/bam/meth_sample.bam"
outdir = "$(dataDir)/cpelnano/"
outdir = "/Users/jordiabante/Desktop/"

# Configuration
min_cov = 2.5; max_em_init = 5; max_em_iters = 100;
config = CpelNano.CpelNanoConfig();
config.verbose = true;
config.out_dir = outdir;
config.out_prefix = "test";

# Analyze each region in chr
CpelNano.analyze_bam(bam,fasta,config)

# Output files
mml_path = "$(outdir)/test_mml.bedGraph"
nme_path = "$(outdir)/test_nme.bedGraph"
theta_path = "$(outdir)/test_theta.bedGraph"

# Print files
# println(String(read(mml_path)))
# println(String(read(nme_path)))
println(String(read(theta_path)))


text:
# m6.6s.jl

using Pkg, DrWatson

begin
    using MonteCarloMeasurements
    using StanSample
    using StatisticalRethinking
end

begin
    N = 100
    df = DataFrame(
    :h0 => rand(Normal(10,2 ), N),
    :treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2)))
    );
    df[!, :fungus] = [rand(Binomial(1, 0.5 - 0.4 * df[i, :treatment]), 1)[1] for i in 1:N]
    df[!, :h1] = [df[i, :h0] + rand(Normal(5 - 3 * df[i, :fungus]), 1)[1] for i in 1:N]
    scale!(df, [:h0, :treatment, :fungus, :h1])
end

sim_p = DataFrame(:sim_p => rand(LogNormal(0, 0.25), 10000))

stan6_6 = "
data {
    int <lower=1> N;
    vector[N] h0;
    vector[N] h1;
}
parameters{
    real<lower=0> p;
    real<lower=0> sigma;
}
model {
    vector[N] mu;
    p ~ lognormal(0, 0.25);
    sigma ~ exponential(1);
    mu = h0 * p;
    h1 ~ normal(mu, sigma);
}
";

m6_6_data = Dict(
    :N => nrow(df),
    :h0 => df[:, :h0],
    :h1 => df[:, :h1]
)
m6_6s = SampleModel("m6.6s", stan6_6)
rc6_6s = stan_sample(m6_6s; data=m6_6_data)

if success(rc6_6s)
    part6_6s = read_samples(m6_6s, :particles);
    part6_6s |> display
end

# End of m6.6s.jl


text:
"""
$(TYPEDEF)

Solver parameters for the DistributedFactoGraph.

Dev Notes
- FIXME change to using kwargs from Parameters.jl
- TODO remove NothingUnion
- TODO Upgrade to common @kwargs struct approach
"""
Base.@kwdef mutable struct SolverParams <: DFG.AbstractParams
  dimID::Int = 0
  reference::NothingUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}} = nothing
  stateless::Bool = false
  qfl::Int = (2^(Sys.WORD_SIZE-1)-1)# Quasi fixed length
  isfixedlag::Bool = false          # true when adhering to qfl window size for solves
  limitfixeddown::Bool = false      # if true, then fixed lag will not update marginalized during down pass on tree
  incremental::Bool = true          # use incremental tree updates, TODO consolidate with recycling
  useMsgLikelihoods::Bool = false   # Experimental, insert differential factors from upward joints
  upsolve::Bool = true              # do tree upsolve
  downsolve::Bool = true            # do tree downsolve
  drawtree::Bool = false            # draw tree during solve
  drawCSMIters::Bool = true         # show CSM iteration count on tree visualization
  showtree::Bool = false
  drawtreerate::Float64 = 0.5       # how fast should the tree vis file be redrawn
  dbg::Bool = false                 # Experimental, enable additional tier debug features
  async::Bool = false               # do not block on CSM tasks
  limititers::Int = 500             # limit number of steps CSMs can take
  N::Int = 100                      # default number of particles
  multiproc::Bool = 1 < nprocs()    # should Distributed.jl tree solve compute features be used
  logpath::String = "/tmp/caesar/$(now())" # unique temporary file storage location for a solve
  graphinit::Bool = true            # default to graph-based initialization of variables
  treeinit::Bool =false             # init variables on the tree
  limittreeinit_iters::Int = 10
  algorithms::Vector{Symbol} = [:default] # list of algorithms to run [:default] is mmisam
  spreadNH::Float64 = 3.0           # entropy spread adjustment used for both null hypo cases.
  inflation::Float64 = 5.0          # how much to disperse particles before convolution solves, #1051
  nullSurplusAdd::Float64 = 0.3     # minimum nullhypo for relative factors sibling to multihypo factors onto a specific variable.
  inflateCycles::Int = 3            # repeat convolutions for inflation to occur
  gibbsIters::Int = 3               # number of Gibbs cycles to take per clique iteration variables
  maxincidence::Int = 500           # maximum incidence to a variable in an effort to enhance sparsity
  alwaysFreshMeasurements::Bool = true # Development feature on whether new samples should be sampled at each Gibbs cycle convolution
  attemptGradients::Bool = false    # should factor gradients be calculated or attempted (UNDER DEVELOPMENT, 21Q3)
  devParams::Dict{Symbol,String} = Dict{Symbol,String}() # empty container for new features, allowing workaround for breaking changes and legacy
  #
end

#


text:
function build_marketclearance!(m, jm)
    # Add market clearance constraints

    commodities = Set{CommodityRef}()
    for pf in m._productions
        for input in pf.inputs
            push!(commodities, input.commodity)
        end

        for output in pf.outputs
            push!(commodities, output.commodity)
        end
    end

    # Loop over commodities
    for commodity in commodities
        # Find all endowments for current commodity
        endows = []
        for demand_function in m._demands
            for en in demand_function.endowments
                if en.commodity==commodity
                    push!(endows, :($(en.quantity)))
                end
            end
        end

        # Find all intermediate supplies for current commodity
        comp_supplies = []
        for production_function in m._productions
            for output in production_function.outputs
                if output.commodity==commodity
                    push!(comp_supplies, :($(get_jump_variable_for_sector(jm, production_function.sector)) * $(get_jump_variable_for_intermediate_supply(jm, output))))
                end
            end
        end

        # Find all final demand for current commodity
        final_demand = []
        for demand_function in m._demands
            for demand in demand_function.demands
                if demand.commodity == commodity
                    push!(final_demand, :($(get_jump_variable_for_final_demand(jm, demand))))
                end
            end
        end

        # Find all the intermediate demands for current commodity
        comp_demands = []
        for production_function in m._productions
            for input in production_function.inputs
                if input.commodity==commodity
                    push!(comp_demands, :($(get_jump_variable_for_sector(jm, production_function.sector)) * $(get_jump_variable_for_intermediate_demand(jm, input))))
                end
            end
        end

        exa = :(
            JuMP.@NLexpression(
                $jm,
                +(0., $(endows...), $(comp_supplies...)) - +(0., $(final_demand...), $(comp_demands...))
            )
        )
        exb = eval(swap_our_param_with_jump_param(jm, exa))

        Complementarity.add_complementarity(jm, get_jump_variable_for_commodity(jm, commodity), exb, string("F_", get_name(commodity, true)))
        push!(m._nlexpressions, exb)
    end
end


text:
function maketoc(path::String; depth=2, title = "# Table of Contents")
    s = "$title\n\n"
    for item in Markdown.parse(read(path, String)).content
        for level in 1:depth
            if isa(item, Markdown.Header{level})
                header = strip(replace(Markdown.plain(item), "#" => ""))
                link = make_toc_link(header)
                s *= repeat(' ', (level-1)*4) * "- [$header](#$link)\n"
            end
        end
    end
    s * "\n\n---\n\n"
end

function make_toc_link(header)
    s = lowercase(replace(header, " " => "-"))
    s = lowercase(replace(s, "?" => ""))
end


text:
Pkg.add("Plots")
Pkg.add("DataFrames")
Pkg.add("ScikitLearn")
Pkg.add("CSV")
Pkg.add("PyCall")
using PyCall
using Plots
using ScikitLearn
using CSV
using DataFrames
np=pyimport("numpy")
data=CSV.read("D:\\Datasets-master\\archive\\iris.csv",DataFrame)
select!(data,Not("Id"))
names(data)
describe(data)
plot(data.SepalLengthCm)
plot(data.SepalWidthCm)
plot(data.PetalLengthCm)
plot(data.PetalWidthCm)
scatter(data.SepalLengthCm,data.SepalWidthCm)
scatter(data.SepalLengthCm,data.PetalLengthCm)
scatter(data.SepalLengthCm,data.PetalWidthCm)
scatter(data.SepalWidthCm,data.Species)
pie(data.SepalLengthCm)
iris_x=convert(Array,data[:,[1,2,3,4]])
iris_y=convert(Array,data[:,5])
@sk_import model_selection: train_test_split
x_train, x_test, y_train, y_test=train_test_split(iris_x,iris_y,random_state=1)
@sk_import linear_model:LogisticRegression
log_reg=LogisticRegression()
fit!(log_reg,x_train,y_train)
prediction=predict(log_reg,x_test)
@sk_import metrics: accuracy_score
accuracy=accuracy_score(prediction,y_test)
print(accuracy)

max_stars_repo_name: abhishek-parashar/Right-From-Scratch

max_issues_repo_name: abhishek-parashar/Right-From-Scratch

max_forks_repo_name: abhishek-parashar/Right-From-Scratch


text:
using LiterateOrg
tangle_package(joinpath(dirname(@__FILE__), "..", "src", "PPT.org"), "PPT")


text:
using Turing
using Test

xs = [1.5 2.0]
# xx = 1

@model fggibbstest(xs) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  # xx ~ Normal(m, sqrt(s)) # this is illegal

  for i = 1:length(xs)
    xs[i] ~ Normal(m, sqrt(s))
  # for xx in xs
    # xx ~ Normal(m, sqrt(s))
  end
  s, m
end

gibbs = Gibbs(2, PG(10, 2, :s), HMC(1, 0.4, 8, :m))
chain = sample(fggibbstest(xs), gibbs);

#
#
# using Turing
# using Test
#
# @model ttt(a) = begin
#   a ~ Normal(0, 1)
#   bb ~ Normal(0, 1)
#   bb
# end
#
# a = 1
#
# sample(ttt(a), SMC(10));


text:
using Revise, FourierGPE
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

# Units
# this example works in oscillator units

function jrphase(x,y)
    return 4x/(1.5 + (2y^2)+x^2)
end
## Initialize simulation

L = (30.0,30.0)
N = (128,128)
sim = Sim(L,N)
@unpack_Sim sim

## simulation parameters
μ = 60.0
ξ = 1/sqrt(μ)
offset = 2.0
# potential
import FourierGPE.V
V1 = 5*μ
sig = 3*ξ
V(x,y,t) = 0.5*(x^2 + y^2) + V1*exp(-((x-offset)^2 + y^2)/sig^2)

# TF state
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

# make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)

@pack_Sim! sim

## evolve and plot final state
sol = runsim(sim)
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)

## imprint and evolve JR for short time
# shift = sig
# ψjr = @. ψg*exp(-im*jrphase((x - offset-shift)/ξ,y'/ξ))
# showpsi(x,y,ψjr)
# ϕi = kspace(ψjr,sim)
# γ = 0.0
# tf = .3
# t = LinRange(ti,tf,Nt)
#
# simjr = Sim(L,N)
# @pack! simjr = tf,t,γ,ϕi

## try jump of obstacle
shift = -2*sig
# ψjr = @. ψg*exp(-im*jrphase((x - offset-shift)/ξ,y'/ξ))

V(x,y,t) = 0.5*(x^2 + y^2) + V1*exp(-((x-offset-shift)^2 + y^2)/sig^2)

ϕi = kspace(ψg,sim)
γ = 0.0
tf = .3
t = LinRange(ti,tf,Nt)

simjr = Sim(L,N)
@pack! simjr = tf,t,γ,ϕi

## evolve
soljr = runsim(simjr)

ϕf = soljr[end]
ψf = xspace(ϕf,simjr)
showpsi(x,y,ψf)

anim = @animate for i=1:Nt-6
    ψ = xspace(soljr[i],simjr)
    showpsi(x,y,ψ)
end

gif(anim,"./examples/jr.gif"),fps=30)


text:
function createClassLabels(labels; outputLevel = 0)

	numDatapoints = length(labels)

	outputLevel > 0  && println("Creating class labels")

	classesDict=Dict{Any,Tuple{Int64,Int64}}()

	currentClass =0
	classLabels = zeros(Float64,numDatapoints )
	for i in 1:numDatapoints
		if haskey(classesDict, labels[i])
			classesDict[labels[i]]=(classesDict[labels[i]][1],classesDict[labels[i]][2]+ 1)
			classLabels[i] = classesDict[labels[i]][1]
		else
			currentClass +=1
			classesDict[labels[i]] = (currentClass,1.0)
			classLabels[i] = currentClass
		end
	end
	outputLevel > 1  && println(classesDict)
	outputLevel > 0  && println("Number of classes is $(length(classesDict))")

	(classLabels, length(classesDict))
end

# Normalize to [-1,1]. This affects sparsity!
function normalizeFeatures(features, numFeatures)
	for j in 1:numFeatures
		thismax = maximum(features[:,j])
		thismin = minimum(features[:,j])
		features[:,j] =(2* features[:,j]-thismax-thismin)/(thismax-thismin)
	end
end

function trainTestRandomSeparate(features,labels::MinusPlusOneVector, labels2::Vector{Float64}; outputLevel = 0)
	srand(1)
	shuffledIndices = shuffle(collect(1 : length(labels)))
	numTP = div(length(labels)*3,4)
	outputLevel > 0  && println("numTP = $numTP")
	outputLevel > 0  && println("numTestingPoints = $(length(labels)-numTP)")
	(copy(features[shuffledIndices[1:numTP      ],:]),
	 MinusPlusOneVector(labels[shuffledIndices[1:numTP      ]  ]),
     copy(labels2[shuffledIndices[1:numTP      ]  ]),
	 numTP,
	 copy(features[shuffledIndices[(numTP+1):end],:]),
	 MinusPlusOneVector(labels[  shuffledIndices[(numTP+1):end]  ]),
	copy(labels2[  shuffledIndices[(numTP+1):end]  ]))
end

# function StrainTestRandomSeparate(features,labels)
# 	srand(1)
# 	shuffledIndices = shuffle([1 : size(labels)[1];])
# 	numTP = div(size(labels)[1]*3,4)
# 	# println(numTP)
# 	# println(shuffledIndices[1:numTP      ])
# 	# println(features)
# 	(copy(features[shuffledIndices[1:numTP      ]]),
# 	copy(labels  [shuffledIndices[1:numTP      ]  ]),
# 	numTP,
# 	copy(features[shuffledIndices[(numTP+1):end]]),
# 	copy(labels[  shuffledIndices[(numTP+1):end]  ]))
# end


text:
using BenchmarkTools

if length(ARGS) != 2
    println("Usage: $PROGRAM_FILE <filename> <count>")
	exit()
end

versioninfo()
println()

tic()
using SASLib
@printf "Loaded library in %.3f seconds\n" toq()

b = @benchmark readsas($ARGS[1], verbose_level=0, string_array_fn=Dict(:_all_ => REGULAR_STR_ARRAY)) samples=parse(Int, ARGS[2])
display(b)

println()
println()


text:
################################################################################
# Types of Potentials
################################################################################

include("lj.jl")
include("bm.jl")
include("coulomb.jl")
include("zbl.jl")
include("morse.jl")

export LennardJones, BornMayer, Coulomb, ZBL
export get_parameters, get_hyperparameters, serialize_parameters, serialize_hyperparameters

################################################################################
# InteratomicPotentials API implmentations for emperical potentials
################################################################################

function energy_and_force(A::AbstractSystem, p::EmpiricalPotential)
    nnlist = neighborlist(A, p.rcutoff)

    e = 0.0
    f = fill(SVector{3}(zeros(3)), length(A))
    for ii in 1:length(A)
        for (jj, r, R) in zip(nnlist.j[ii], nnlist.r[ii], nnlist.R[ii])
            species = unique([atomic_symbol(A, ii), atomic_symbol(A, jj)])
            if (intersect(species, p.species) == species)
                e += potential_energy(R, p)
                fo = force(R, r, p)
                f[ii] = f[ii] + fo
                f[jj] = f[jj] - fo
            end
        end
    end
    (; e, f)
end

force(r::SVector{3,<:AbstractFloat}, p::EmpiricalPotential) = force(norm(r), r, p)

function virial_stress(A::AbstractSystem, p::EmpiricalPotential)
    nnlist = neighborlist(A, p.rcutoff)

    v = SVector{6}(zeros(6))
    for ii in 1:length(A)
        for (r, R) in zip(nnlist.r[ii], nnlist.R[ii])
            fo = force(R, r, p)
            vi = r * fo'
            v = v + [vi[1, 1], vi[2, 2], vi[3, 3], vi[3, 2], vi[3, 1], vi[2, 1]]
        end
    end
    v
end


text:
using Knet, ArgParse, Base.Test

# Uncomment these if you want lots of messages:
# import Base.Test: default_handler, Success, Failure, Error
# default_handler(r::Success) = info("$(r.expr)")
# default_handler(r::Failure) = warn("FAIL: $(r.expr)")
# default_handler(r::Error)   = warn("$(r.err): $(r.expr)")

load_only = true

s = ArgParseSettings()
@add_arg_table s begin
    ("--all"; action=:store_true)
    ("--twice"; action=:store_true)
    ("--gcheck"; arg_type=Int; default=0)
    ("--linreg"; action=:store_true)
    ("--mnist2d"; action=:store_true)
    ("--mnist2dy"; action=:store_true)
    ("--mnist2dx"; action=:store_true)
    ("--mnist2dxy"; action=:store_true)
    ("--mnist4d"; action=:store_true)
    ("--addirnn"; action=:store_true)
    ("--addlstm"; action=:store_true)
    ("--mnistpixels"; action=:store_true)
    ("--rnnlm"; action=:store_true)
    ("--copyseq"; action=:store_true)
    ("--ncelm"; action=:store_true)
    ("--ner"; action=:store_true)
end
opts = parse_args(isempty(ARGS) ? ["--all"] : ARGS, s)
gcheck = opts["gcheck"]
twice = opts["twice"]

if opts["all"] || opts["linreg"]
    include("linreg.jl")
    @time @show test1 = LinReg.main("--gcheck $gcheck")
    #@test test1 == (0.0005497372347062405,32.77256166946498,0.11244349406523031)
    #@test test1 == (0.0005497372347062409,32.77256166946497,0.11244349406522969) # Mon Oct 26 11:10:17 PDT 2015: update uses axpy to scale with gclip&lr
    #@test test1 == (0.0005497846637255734,32.77257400591496,0.11265426200775067) # Wed Nov 18 21:39:18 PST 2015: xavier init
    @test test1  == (0.0005497846637255735,32.77257400591496,0.11265426200775067) # Mon Dec  7 22:40:27 PST 2015: cpu quadloss
    twice && (gc(); @time @show test1 = LinReg.main(""))
    # 0.739858 seconds (394.09 k allocations: 71.335 MB, 1.23% gc time) Tue Oct 20 18:29:41 PDT 2015
    # 0.731515 seconds (391.62 k allocations: 71.451 MB, 1.21% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 0.837992 seconds (690.89 k allocations: 86.433 MB, 0.80% gc time) Mon Dec  7 22:43:42 PST 2015: knet7
    # 0.661330 seconds (408.42 k allocations: 77.432 MB, 0.89% gc time) 466c40e not using stack for fnn
    # 0.762190 seconds (391.60 k allocations: 71.451 MB, 1.50% gc time) 5f60a84 latest master
    # 0.674939 seconds (418.42 k allocations: 77.355 MB, 1.24% gc time) ad5d969 latest Knet7
    # 0.577099 seconds (356.36 k allocations: 73.864 MB, 0.94% gc time) 53bfcd1 broadcast.jl
end

if opts["all"] || opts["mnist2d"]
    include("mnist2d.jl")
    @time @show test2 = MNIST2D.main("--gcheck $gcheck")
    # @test test2 == (0.10628127f0,24.865438f0,3.5134742f0)
    # @test test2 == (0.10626979f0,24.866688f0,3.5134728f0) # softloss with mask
    # @test test2 == (0.10610757f0,24.87226f0,3.3128357f0)  # 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    # @test test2 == (0.108679175f0,43.344624f0,3.0013776f0) # Wed Nov 18 21:39:18 PST 2015: xavier init
    # @test test2 == (0.10911515f0,43.322735f0,3.2024312f0)  # improved softmax (dx=q-p)
    # @test test2 == (0.10843158f0,36.997395f0,2.9633076f0)    # switched back to gaussian, xavier slows down xsparse!?
    @test test2 == (0.10843158f0,36.997395f0,2.963311f0)     # softmaxtemp introduced
    twice && (gc(); @time @show test2 = MNIST2D.main("--gcheck $gcheck"))
    # 6.941715 seconds (3.35 M allocations: 151.876 MB, 1.33% gc time) Tue Oct 20 19:15:59 PDT 2015
    # 6.741272 seconds (3.35 M allocations: 151.858 MB, 1.41% gc time) Mon Oct 26 11:10:17 PDT 2015: update uses axpy to scale with gclip&lr
    # 7.031675 seconds (3.47 M allocations: 158.983 MB, 1.31% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 6.730379 seconds (3.61 M allocations: 161.998 MB, 1.42% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
    # 4.832591 seconds (4.67 M allocations: 407.652 MB, 4.33% gc time) Tue Dec  8 09:36:22 PST 2015: knet7
    # 4.096904 seconds (3.28 M allocations: 362.141 MB, 4.84% gc time) Wed Dec  9 16:10:19 PST 2015: 466c40e not using stack for fnn
    # 7.656945 seconds (3.54 M allocations: 159.374 MB, 1.28% gc time) 5f60a84 latest master (xavier)
    # 4.219116 seconds (3.35 M allocations: 364.417 MB, 5.91% gc time) ad5d969 latest Knet7
    # 3.611838 seconds (2.95 M allocations: 346.479 MB, 5.23% gc time) 53bfcd1 broadcast.jl
end

if opts["all"] || opts["mnist2dy"]
    isdefined(:mnist2d) || include("mnist2d.jl")
    @time @show test3 = MNIST2D.main("--ysparse --gcheck $gcheck")
    #@test test3 == (0.1062698f0,24.866688f0,3.513474f0)
    #@test test3 == (0.10610757f0,24.87226f0,3.3128357f0)  # 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    #@test test3 == (0.10906915f0,43.341377f0,3.1931002f0)  # Wed Nov 18 21:39:18 PST 2015: xavier init
    #@test test3 == (0.10911515f0,43.322735f0,3.2024312f0) # improved softmax (dx=q-p)
    #@test test3 == (0.10843158f0,36.997395f0,2.9633076f0)    # switched back to gaussian, xavier slows down xsparse!?
    @test test3 == (0.10843156f0,36.997395f0,2.963311f0)     # softmaxtemp introduced
    twice && (gc(); @time @show test3 = MNIST2D.main("--ysparse --gcheck $gcheck"))
    # 8.478264 seconds (3.59 M allocations: 173.689 MB, 2.06% gc time) Tue Oct 20 19:14:45 PDT 2015
    # 8.205758 seconds (3.59 M allocations: 173.636 MB, 2.14% gc time) Mon Oct 26 11:10:17 PDT 2015: update uses axpy to scale with gclip&lr
    # 8.542426 seconds (3.73 M allocations: 181.290 MB, 2.06% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 8.073397 seconds (3.86 M allocations: 184.237 MB, 2.15% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
    # 5.237889 seconds (4.86 M allocations: 432.598 MB, 5.14% gc time) Tue Dec  8 08:06:30 PST 2015: knet7
    # 4.370380 seconds (3.48 M allocations: 387.233 MB, 5.65% gc time) 466c40e not using stack for fnn
    # 9.067960 seconds (3.71 M allocations: 184.288 MB, 1.74% gc time) 5f60a84 latest master (xavier)
    # 4.605455 seconds (3.56 M allocations: 389.701 MB, 6.76% gc time) ad5d969 latest Knet7
    # 3.959572 seconds (3.31 M allocations: 371.077 MB, 6.99% gc time) 53bfcd1 broadcast.jl
end

if opts["all"] || opts["mnist2dx"]
    isdefined(:mnist2d) || include("mnist2d.jl")
    @time @show test4 = MNIST2D.main("--xsparse --gcheck $gcheck")

    # @test isapprox(test4[1], 0.10628127f0; rtol=0.005)
    # @test isapprox(test4[2], 24.865437f0; rtol=0.002)
    # @test isapprox(test4[3], 3.5134742f0; rtol=0.02) # cannot compute csru vecnorm

    # 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    # @test isapprox(test4[1], 0.10610757f0; rtol=0.01)
    # @test isapprox(test4[2], 24.87226f0; rtol=0.001)
    # @test isapprox(test4[3], 3.3128357f0; rtol=0.1)

    # (0.109181836f0,43.36849f0,3.2359037f0) # Wed Nov 18 21:39:18 PST 2015: xavier init
    # @test isapprox(test4[1], 0.109181836f0; rtol=0.01)
    # @test isapprox(test4[2], 43.36849f0; rtol=0.001)
    # @test isapprox(test4[3], 3.2359037f0; rtol=0.1)

    # (0.109118156f0,43.322895f0,3.2024305f0) # improved softmax (dx=q-p)
    # @test isapprox(test4[1], 0.109118156f0; rtol=0.01)
    # @test isapprox(test4[2], 43.322895f0; rtol=0.001)
    # @test isapprox(test4[3], 3.2024305f0; rtol=0.1)

    # (0.10815071f0,37.009987f0,2.8580344f0) TODO: switched back to gaussian, xavier in final layer slows down xsparse!?
    @test isapprox(test4[1], 0.10815071f0; rtol=0.01)
    @test isapprox(test4[2], 37.009987f0; rtol=0.001)
    @test isapprox(test4[3], 2.8580344f0; rtol=0.1)

    twice && (gc(); @time @show test4 = MNIST2D.main("--xsparse --gcheck $gcheck"))
    # 12.362125 seconds (3.81 M allocations: 753.744 MB, 1.87% gc time) Tue Oct 20 19:13:25 PDT 2015
    # 11.751002 seconds (3.84 M allocations: 753.959 MB, 1.95% gc time) Mon Oct 26 11:10:17 PDT 2015: update uses axpy to scale with gclip&lr
    # 12.005169 seconds (3.95 M allocations: 761.003 MB, 1.90% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 12.647064 seconds (4.10 M allocations: 764.156 MB, 1.71% gc time) 2a9c883 2015-11-18 fixed lstm infersize issue
    # 19.643433 seconds (4.08 M allocations: 750.880 MB, 1.12% gc time) 405b7d7 2015-11-18 Changed default weight init to Xavier
    # 16.538108 seconds (5.16 M allocations: 909.466 MB, 1.65% gc time) Tue Dec  8 08:07:17 PST 2015: knet7 xavier
    # 9.608017 seconds (5.16 M allocations: 921.080 MB, 3.53% gc time)  Tue Dec  8 09:39:05 PST 2015: knet7 gaussian
    # 7.946223 seconds (3.76 M allocations: 875.557 MB, 3.11% gc time)  466c40e not using stack for fnn
    # 19.994696 seconds (4.01 M allocations: 748.440 MB, 1.39% gc time) 5f60a84 latest master (xavier)
    # 8.290566 seconds (3.83 M allocations: 877.837 MB, 3.95% gc time)  ad5d969 latest Knet7
    # 7.527577 seconds (3.50 M allocations: 861.132 MB, 3.36% gc time)  53bfcd1 broadcast.jl
end

if opts["all"] || opts["mnist2dxy"]
    isdefined(:mnist2d) || include("mnist2d.jl")
    @time @show test5 = MNIST2D.main("--xsparse --ysparse --gcheck $gcheck")
    # @test isapprox(test5[1], 0.10628127f0; rtol=0.005)
    # @test isapprox(test5[2], 24.865437f0; rtol=0.002)
    # @test isapprox(test5[3], 3.5134742f0; rtol=0.02) # cannot compute csru vecnorm

    # 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    # @test isapprox(test5[1], 0.10610757f0; rtol=0.01)
    # @test isapprox(test5[2], 24.87226f0; rtol=0.001)
    # @test isapprox(test5[3], 3.3128357f0; rtol=0.1)

    # # (0.1091015f0,43.368706f0,3.2359257f0); Wed Nov 18 21:39:18 PST 2015: xavier init
    # @test isapprox(test5[1], 0.1091015f0; rtol=0.01)
    # @test isapprox(test5[2], 43.368706f0; rtol=0.001)
    # @test isapprox(test5[3], 3.2359257f0; rtol=0.1)

    # (0.109118156f0,43.322895f0,3.2024305f0) # improved softmax (dx=q-p)
    # @test isapprox(test5[1], 0.109118156f0; rtol=0.01)
    # @test isapprox(test5[2], 43.322895f0; rtol=0.001)
    # @test isapprox(test5[3], 3.2024305f0; rtol=0.1)

    # (0.10815071f0,37.009987f0,2.8580344f0) switched back to gaussian, xavier in final layer slows down xsparse!?
    @test isapprox(test5[1], 0.10815071f0; rtol=0.01)
    @test isapprox(test5[2], 37.009987f0; rtol=0.001)
    @test isapprox(test5[3], 2.8580344f0; rtol=0.1)

    twice && (gc(); @time @show test5 = MNIST2D.main("--xsparse --ysparse --gcheck $gcheck"))
    # 14.077099 seconds (4.09 M allocations: 776.263 MB, 2.22% gc time) Tue Oct 20 19:11:52 PDT 2015
    # 13.320959 seconds (4.11 M allocations: 776.397 MB, 2.29% gc time) Mon Oct 26 11:10:17 PDT 2015: update uses axpy to scale with gclip&lr
    # 13.339761 seconds (4.23 M allocations: 783.602 MB, 2.27% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 13.421199 seconds (4.37 M allocations: 786.728 MB, 2.27% gc time) Wed Nov 18 21:28:22 PST 2015: lcn gaussian
    # 16.959799 seconds (5.43 M allocations: 935.729 MB, 1.98% gc time) Tue Dec  8 08:07:43 PST 2015: knet7 xavier
    # 10.177581 seconds (5.42 M allocations: 947.233 MB, 4.17% gc time) Tue Dec  8 09:40:31 PST 2015: knet7 gaussian
    # 8.360521 seconds (4.01 M allocations: 901.532 MB, 3.63% gc time)  466c40e not using stack for fnn
    # 21.294241 seconds (4.22 M allocations: 774.002 MB, 1.58% gc time) 5f60a84 latest master (xavier)
    # 8.688249 seconds (4.10 M allocations: 904.181 MB, 4.55% gc time)  ad5d969 latest Knet7
    # 8.237475 seconds (3.87 M allocations: 885.827 MB, 4.08% gc time)  53bfcd1 broadcast.jl
end

if opts["all"] || opts["mnist4d"]
    include("mnist4d.jl")
    @time @show test6 = MNIST4D.main("--gcheck $gcheck")

    # @test isapprox(test6[1], 0.050180204f0; rtol=.01)
    # @test isapprox(test6[2], 25.783848f0;   rtol=.01)
    # @test isapprox(test6[3], 9.420588f0;    rtol=.1)

    # 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    # @test isapprox(test6[1], .050181; rtol=0.01)
    # @test isapprox(test6[2], 25.7783; rtol=0.001)
    # @test isapprox(test6[3], 9.59026; rtol=0.1)

    # # (0.02938571214979068,65.9176025390625,7.652309894561768); Wed Nov 18 21:39:18 PST 2015: xavier init
    # @test isapprox(test6[1], .029385; rtol=0.01)
    # @test isapprox(test6[2], 65.9176; rtol=0.001)
    # @test isapprox(test6[3], 7.65231; rtol=0.1)

    # (0.03091877909648853,66.5221939086914,4.423511505126953); improved softmax (dx=q-p)
    @test isapprox(test6[1], .030919; rtol=0.02)
    @test isapprox(test6[2], 66.5222; rtol=0.001)
    @test isapprox(test6[3], 5.0; rtol=0.2)

    twice && (gc(); @time @show test6 = MNIST4D.main("--gcheck $gcheck"))
    # 17.093371 seconds (10.15 M allocations: 479.611 MB, 1.11% gc time) Tue Oct 20 19:09:19 PDT 2015
    # 17.135514 seconds (10.38 M allocations: 494.816 MB, 1.11% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 17.002958 seconds (10.58 M allocations: 499.822 MB, 1.11% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
    # 14.898975 seconds (10.11 M allocations: 697.407 MB, 1.97% gc time) 466c40e not using stack for fnn
    # 16.956578 seconds (10.41 M allocations: 495.519 MB, 1.28% gc time) 5f60a84 latest master (xavier)
    # 14.990415 seconds (10.23 M allocations: 701.862 MB, 2.59% gc time) ad5d969 latest Knet7
    # 14.402732 seconds (9.66 M allocations: 674.339 MB, 1.97% gc time)  53bfcd1 broadcast.jl
end

if opts["all"] || opts["addirnn"]
    include("adding.jl")
    @time @show test8b = Adding.main("") # gcheck fails with relu

    # @test test8b == (0.04885713f0, 5.6036315f0,3.805253f0)  	# --epochs 20 --nettype irnn
    # @test test8b == (0.04885713f0, 5.6057444f0, 3.805253f0) 	# measuring wnorm after update now
    # @test test8b == (0.05627571f0,5.484082f0,4.1594324f0)    	# new generator
    # @test test8b == (0.056275677f0,5.484083f0,4.159457f0)	# 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    @test test8b ==   (0.05627568f0,5.484083f0,4.159457f0) 	# 0c1148e switched to cpu quadloss

    twice && (gc(); @time @show test8b = Adding.main(""))
    # 9.114330 seconds (16.23 M allocations: 704.629 MB, 1.80% gc time) # --epochs 20 --nettype irnn
    # 10.703243 seconds (20.59 M allocations: 863.693 MB, 2.41% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 10.528267 seconds (21.14 M allocations: 876.542 MB, 2.01% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
    # 12.071895 seconds (21.19 M allocations: 922.258 MB, 1.80% gc time) 0c1148e knet7 expensive initforw
    # 11.802089 seconds (21.81 M allocations: 901.486 MB, 3.06% gc time) 5f60a84 latest master
    # 12.800159 seconds (21.51 M allocations: 935.733 MB, 2.57% gc time) ad5d969 latest Knet7
    # 11.881547 seconds (21.80 M allocations: 935.974 MB, 1.90% gc time) 53bfcd1 broadcast.jl
end

if opts["all"] || opts["addlstm"]
    include("adding.jl")
    @time @show test8 = Adding.main("--gcheck $gcheck --epochs 1 --nettype lstm")

    @test test8 == (0.24768005f0,3.601481f0,1.2290705f0)      # switched to --epochs 1 --nettype lstm, gradcheck does not work with irnn/relu

    twice && (gc(); @time @show test8 = Adding.main("--gcheck $gcheck --epochs 1 --nettype lstm"))
    # 2.028728 seconds (3.82 M allocations: 164.958 MB, 1.95% gc time)  Tue Oct 20 19:03:01 PDT 2015
    # 2.246450 seconds (3.90 M allocations: 169.582 MB, 2.18% gc time)  Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 2.557893 seconds (4.20 M allocations: 189.816 MB, 1.68% gc time)  08758b7 Knet7
    # 2.482413 seconds (4.18 M allocations: 178.685 MB, 2.97% gc time)  5f60a84 latest master (xavier)
    # 2.676177 seconds (4.27 M allocations: 192.726 MB, 2.35% gc time)  ad5d969 latest Knet7
    # 2.403824 seconds (4.06 M allocations: 181.446 MB, 1.85% gc time)  53bfcd1 broadcast.jl
end

if opts["all"] || opts["mnistpixels"]
    include("mnistpixels.jl")

    # @time @show test7 = MNISTPixels.main("--gcheck $gcheck")
    # 9.909841 seconds (45.76 M allocations: 1.208 GB, 3.52% gc time)
    # 8.877034 seconds (43.27 M allocations: 1.099 GB, 4.33% gc time)
    # @test test7  == (0.1216,2.3023171f0,10.4108f0,30.598776f0)
    # @test test7 == (0.12159999999999982,2.3023171f0,10.4108f0,30.598776f0) # switched to itembased
    # @test test7 == (0.12159999999999982,2.3023171f0,10.412794f0,30.598776f0) # measuring wnorm after update now

    # switch to lstm so we can gradcheck, too slow for gcheck>1
    @time @show test7 = MNISTPixels.main("--gcheck $gcheck --nettype lstm --testfreq 2 --epochs 1 --batchsize 64 --epochsize 128")
    @test test7 == (0,2.3025737f0,14.70776f0,0.12069904f0) # switched to --gcheck 1 --nettype lstm --testfreq 2 --epochs 1 --batchsize 64 --epochsize 128

    twice && (gc(); @time @show test7 = MNISTPixels.main("--gcheck $gcheck --nettype lstm --testfreq 2 --epochs 1 --batchsize 64 --epochsize 128"))
    # 2.599979 seconds (5.19 M allocations: 212.248 MB, 2.77% gc time)  Tue Oct 20 19:07:11 PDT 2015
    # 2.713217 seconds (5.28 M allocations: 217.967 MB, 2.56% gc time)  Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 2.889495 seconds (5.09 M allocations: 208.266 MB, 2.43% gc time)  f826bce: Knet7 small diff in result
    # 2.856504 seconds (5.09 M allocations: 208.169 MB, 2.33% gc time)  08758b7: Knet7 diff fixed
    # 3.021688 seconds (5.64 M allocations: 229.431 MB, 3.30% gc time)  5f60a84 latest master (xavier)
    # 3.117822 seconds (4.95 M allocations: 212.537 MB, 4.01% gc time)  ad5d969 latest Knet7
    # 2.726979 seconds (4.58 M allocations: 196.119 MB, 3.06% gc time)  53bfcd1 broadcast.jl
end

const knetdata = Pkg.dir("Knet/data")

function download_mikolov_examples()
    url = "http://www.fit.vutbr.cz/~imikolov/rnnlm/simple-examples.tgz"
    if !isfile("$knetdata/ptb.valid.txt")
        info("Downloading $url to $knetdata")
        curr = pwd()
        cd(knetdata)
	run(pipeline(`wget -q -O- http://www.fit.vutbr.cz/~imikolov/rnnlm/simple-examples.tgz`,
                     `tar --strip-components 3 -xvzf - ./simple-examples/data/ptb.valid.txt ./simple-examples/data/ptb.test.txt`))
        cd(curr)
    end
end

if opts["all"] || opts["rnnlm"]
    include("rnnlm.jl")
    download_mikolov_examples()
    @time @show test9 = RNNLM.main("$knetdata/ptb.valid.txt $knetdata/ptb.test.txt --gcheck $gcheck")

    # This is for: Float64
    # @test isapprox(test9[1], 814.9780887272417;  rtol=.0001)
    # @test isapprox(test9[2], 541.2457922913605;  rtol=.0001)
    # @test isapprox(test9[3], 267.626257438979;   rtol=.005)
    # @test isapprox(test9[4], 120.16170771885587; rtol=.0001)

    # Changing to: Float32
    # @test isapprox(test9[1], 825.336, rtol=0.05)
    # @test isapprox(test9[2], 531.640, rtol=0.05)
    # @test isapprox(test9[3], 267.337, rtol=.005)
    # @test isapprox(test9[4], 136.923, rtol=0.0001)

    # # 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    # @test isapprox(test9[1], 822.177, rtol=0.01)
    # @test isapprox(test9[2], 535.746, rtol=0.1)
    # @test isapprox(test9[3], 267.272, rtol=0.001)
    # @test isapprox(test9[4], 136.923, rtol=0.0001)

    # 5b1eca0 Sat Dec 19 09:30:14 PST 2015 Knet7 passes gcheck
    @test isapprox(test9[1], 824.207, rtol=0.01)
    @test isapprox(test9[2], 532.005, rtol=0.1)
    @test isapprox(test9[3], 267.190, rtol=0.001)
    @test isapprox(test9[4], 136.923, rtol=0.0001)

    twice && (gc(); @time @show test9 = RNNLM.main("$knetdata/ptb.valid.txt $knetdata/ptb.test.txt --gcheck $gcheck"))
    # 32.368835 seconds (22.35 M allocations: 2.210 GB, 1.56% gc time)   for Float64
    # 22.892147 seconds (22.46 M allocations: 945.257 MB, 2.17% gc time) after switching to Float32
    # 21.982870 seconds (20.64 M allocations: 866.929 MB, 3.08% gc time) Tue Oct 20 19:00:29 PDT 2015
    # 22.972406 seconds (21.15 M allocations: 893.519 MB, 3.08% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 23.269379 seconds (21.74 M allocations: 902.787 MB, 3.01% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
    # 23.352192 seconds (19.00 M allocations: 2.850 GB, 3.80% gc time)   9312287: Knet7 with small diff
    # 23.227981 seconds (19.04 M allocations: 2.851 GB, 3.81% gc time)   08758b7: Knet7 fixed bug
    # 23.439377 seconds (21.80 M allocations: 2.970 GB, 4.00% gc time)   5f60a84 latest master (xavier)
    # 26.134600 seconds (18.88 M allocations: 2.859 GB, 4.20% gc time)   ad5d969 latest Knet7
    # 23.591260 seconds (18.73 M allocations: 2.852 GB, 4.09% gc time)   5b1eca0 Sat Dec 19 09:30:14 PST 2015
    # 26.092232 seconds (18.77 M allocations: 2.853 GB, 4.11% gc time)   when run with runtests.jl --all ??
    # 23.319572 seconds (19.01 M allocations: 795.725 MB, 3.02% gc time) fixed sparse softloss, ran by itself
    # 25.047980 seconds (20.10 M allocations: 844.618 MB, 2.97% gc time) 53bfcd1 broadcast.jl
end

if opts["all"] || opts["copyseq"]
    download_mikolov_examples()
    include("copyseq.jl")
    @time @show test10 = CopySeq.main("--epochs 1 --gcheck $gcheck $knetdata/ptb.valid.txt $knetdata/ptb.test.txt")

    # @test isapprox(test10[1], 3143.22; rtol=.001)
    # @test isapprox(test10[2], 1261.19; rtol=.0001)
    # @test isapprox(test10[3], 106.760; rtol=.0001)
    # @test isapprox(test10[4], 206.272; rtol=.0001)

    # 51a1bc1 v0.6.8 improved nan-proofing of softmax Fri Nov  6 12:53:16 PST 2015
    # @test isapprox(test10[1], 3130.86 ; rtol=.0001)
    # @test isapprox(test10[2], 1363.82 ; rtol=.0001)
    # @test isapprox(test10[3],  105.026; rtol=.0001)
    # @test isapprox(test10[4],  184.931; rtol=.0001)

    # (4248.32913889032,959.0360211411156,102.61302185058594,145.56500244140625); Wed Nov 18 21:39:18 PST 2015: xavier init
    # How did this happen when we have a winit argument?
    # @test isapprox(test10[1], 4248.329; rtol=.0001)
    # @test isapprox(test10[2],  959.036; rtol=.0001)
    # @test isapprox(test10[3],  102.613; rtol=.0001)
    # @test isapprox(test10[4],  145.565; rtol=.0001)

    # (2624.353349018743,1325.262392101933,103.99244689941406,224.86605834960938); ad5d969 latest Knet7
    @test isapprox(test10[1], 2624.353; rtol=.001)
    @test isapprox(test10[2], 1325.262; rtol=.001)
    @test isapprox(test10[3],  103.992; rtol=.001)
    @test isapprox(test10[4],  224.866; rtol=.001)

    twice && (gc(); @time @show test10 = CopySeq.main("--epochs 1 --gcheck $gcheck $knetdata/ptb.valid.txt $knetdata/ptb.test.txt"))
    # 5.984980 seconds (8.33 M allocations: 353.611 MB, 4.15% gc time) Tue Oct 20 18:58:25 PDT 2015
    # 11.230476 seconds (16.29 M allocations: 701.612 MB, 4.05% gc time) Wed Oct 21 23:19:24 PDT 2015 (unsorted input)
    # 11.658034 seconds (17.49 M allocations: 752.336 MB, 4.47% gc time) Fri Nov  6 12:53:16 PST 2015: new add kernels
    # 11.743344 seconds (17.82 M allocations: 749.691 MB, 4.43% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
    # 11.432132 seconds (15.06 M allocations: 4.608 GB, 5.12% gc time)   e7e299b Knet7 results differ.
    # 15.257276 seconds (18.37 M allocations: 4.756 GB, 6.38% gc time)   5f60a84 latest master (sparse softloss dx=q-p copies ygold to dense)
    # 15.604067 seconds (16.58 M allocations: 4.699 GB, 7.37% gc time)   ad5d969 latest Knet7 (probably same problem)
    # 12.093022 seconds (16.60 M allocations: 713.095 MB, 4.80% gc time) fixed sparse softloss
    # 11.591728 seconds (15.00 M allocations: 651.078 MB, 3.43% gc time) 53bfcd1 broadcast.jl
end

# if opts["all"] || opts["ncelm"] # TODO: port ncelm to Knet7
#     include("ncelm.jl")
#     if !isfile("ptb.valid.txt")
#         info("Downloading ptb...")
# 	run(pipeline(`wget -q -O- http://www.fit.vutbr.cz/~imikolov/ncelm/simple-examples.tgz`,
#                      `tar --strip-components 3 -xvzf - ./simple-examples/data/ptb.valid.txt ./simple-examples/data/ptb.test.txt`))
#     end
#     @time @show test11 = ncelm("ptb.valid.txt ptb.test.txt --gcheck $gcheck")
#     @test isapprox(test11[1], 1.04277, rtol=0.0001)
#     @test isapprox(test11[2], 1411.14, rtol=0.0001)
#     @test isapprox(test11[3], 968.846, rtol=0.0001)
#     @test isapprox(test11[4], 31.8226, rtol=0.0001)
#     twice && (gc(); @time @show test11 = ncelm("ptb.valid.txt ptb.test.txt --gcheck $gcheck"))
#     # 6.069526 seconds (5.36 M allocations: 204.252 MB, 2.06% gc time)
#     # 6.079928 seconds (5.45 M allocations: 204.452 MB, 2.21% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
#     # 6.146882 seconds (5.60 M allocations: 209.355 MB, 2.60% gc time) 5f60a84 latest master (xavier)
# end

# if (opts["all"] || opts["ner"]) && isfile("ner.jld") # TODO: gcheck for ner, port to Knet7
#     include("ner.jl")
#     @time @show test12 = ner("--devfortrn --epochs 1 --batchsize 128")
#     @test test12 == (1,5.391641813553446,5.146268547771243,0.8020976309565352)
#     twice && (gc(); @time @show test12 = ner("--devfortrn --epochs 1 --batchsize 128"))
#     # 20.866555 seconds (37.80 M allocations: 1.829 GB, 6.12% gc time)
#     # 20.983972 seconds (38.88 M allocations: 1.854 GB, 6.21% gc time) Wed Nov 18 21:28:22 PST 2015: lcn
#     # 20.359684 seconds (40.20 M allocations: 1.908 GB, 5.04% gc time) 5f60a84 latest master (xavier)
# end

# 1    # 0.577099 seconds (356.36 k allocations: 73.864 MB, 0.94% gc time) 53bfcd1 broadcast.jl
# 2    # 3.611838 seconds (2.95 M allocations: 346.479 MB, 5.23% gc time) 53bfcd1 broadcast.jl
# 3    # 3.959572 seconds (3.31 M allocations: 371.077 MB, 6.99% gc time) 53bfcd1 broadcast.jl
# 4    # 7.527577 seconds (3.50 M allocations: 861.132 MB, 3.36% gc time)  53bfcd1 broadcast.jl
# 5    # 8.237475 seconds (3.87 M allocations: 885.827 MB, 4.08% gc time)  53bfcd1 broadcast.jl
# 6    # 14.402732 seconds (9.66 M allocations: 674.339 MB, 1.97% gc time)  53bfcd1 broadcast.jl
# 8b   # 11.881547 seconds (21.80 M allocations: 935.974 MB, 1.90% gc time) 53bfcd1 broadcast.jl
# 8    # 2.403824 seconds (4.06 M allocations: 181.446 MB, 1.85% gc time)  53bfcd1 broadcast.jl
# 7    # 2.726979 seconds (4.58 M allocations: 196.119 MB, 3.06% gc time)  53bfcd1 broadcast.jl
# 9    # 25.047980 seconds (20.10 M allocations: 844.618 MB, 2.97% gc time) 53bfcd1 broadcast.jl
# 10   # 11.591728 seconds (15.00 M allocations: 651.078 MB, 3.43% gc time) 53bfcd1 broadcast.jl


text:
@testset "684.redundant-connection.jl" begin
    @test findRedundantConnection([[3,4],[1,2],[2,4],[3,5],[2,5]]) == [2, 5]
    @test findRedundantConnection([[1,2], [1,3], [2,3]]) == [2, 3]
    @test findRedundantConnection([[1,2], [2,3], [3,4], [1,4], [1,5]]) == [1, 4]
end


text:
# constructs a semantic theory about a commutative monoid
# A monoid whose operation is commutative is called a
# commutative monoid (or, less commonly, an abelian monoid).

include("docstrings.jl")

module Library

using Metatheory.Patterns
using Metatheory.Rules

macro associativity(op)
	quote
		[
			(@left_associative $op),
			(@right_associative $op)
		]
	end
end

macro monoid(op, id)
	quote
		[
			(@left_associative(op)),
			(@right_associative(op)),
			(@identity_left(op, id)),
			(@identity_right(op, id))
		]
	end
end

macro commutative_monoid(op, id)
	quote
		[
			(@commutativity $op),
			(@left_associative $op),
			(@right_associative $op),
			(@identity_left $op $id)
		]
	end
end

# constructs a semantic theory about a an abelian group
# The definition of a group does not require that a ⋅ b = b ⋅ a
# for all elements a and b in G. If this additional condition holds,
# then the operation is said to be commutative, and the group is called an abelian group.
macro commutative_group(op, id, invop)
	# @assert Base.isbinaryoperator(op)
	# @assert Base.isunaryoperator(invop)
	quote
		(@commutative_monoid $op $id) ∪ [@inverse_right $op $id $invop]
	end
end

macro distrib(outop, inop)
	quote
		[
			(@distrib_left $outop $inop),
			(@distrib_right $outop $inop),
		]
	end
end

macro commutativity(op)
	RewriteRule(
		PatTerm(:call, op, [PatVar(:a), PatVar(:b)], __module__),
		PatTerm(:call, op, [PatVar(:b), PatVar(:a)], __module__))
end

macro right_associative(op)
	RewriteRule(
		PatTerm(:call, op, [PatVar(:a),
			PatTerm(:call, op, [PatVar(:b), PatVar(:c)], __module__)], __module__),
		PatTerm(:call, op, [
			PatTerm(:call, op, [PatVar(:a), PatVar(:b)], __module__),
			PatVar(:c),
			], __module__))
end
macro left_associative(op)
	RewriteRule(
		PatTerm(:call, op, [
			PatTerm(:call, op, [PatVar(:a), PatVar(:b)], __module__),
			PatVar(:c),
			], __module__),
		PatTerm(:call, op, [PatVar(:a),
			PatTerm(:call, op, [PatVar(:b), PatVar(:c)], __module__)], __module__))
end

macro identity_left(op, id)
	RewriteRule(PatTerm(:call, op, [id, PatVar(:a)], __module__), PatVar(:a))
end

macro identity_right(op, id)
	RewriteRule(PatTerm(:call, op, [PatVar(:a), id], __module__), PatVar(:a))
end

macro inverse_left(op, id, invop)
	RewriteRule(PatTerm(:call, op, [
		PatTerm(:call, invop, [PatVar(:a)], __module__), PatVar(:a)], __module__), id)
end
macro inverse_right(op, id, invop)
	RewriteRule(PatTerm(:call, op, [
		PatVar(:a),
		PatTerm(:call, invop, [PatVar(:a)], __module__)], __module__), id)
end

# distributivity of two operations
# example: `@distrib (⋅) (⊕)`
macro distrib_left(outop, inop)
	EqualityRule(
		# left
		PatTerm(:call, outop, [
			PatVar(:a),
			PatTerm(:call, inop, [PatVar(:b), PatVar(:c)], __module__)
		], __module__),
		# right
		PatTerm(:call, inop, [
			PatTerm(:call, outop, [PatVar(:a), PatVar(:b)], __module__),
			PatTerm(:call, outop, [PatVar(:a), PatVar(:c)], __module__),
		], __module__))

end

macro distrib_right(outop, inop)
	EqualityRule(
		# left
		PatTerm(:call, outop, [
			PatTerm(:call, inop, [PatVar(:a), PatVar(:b)], __module__),
			PatVar(:c)
		], __module__),
		# right
		PatTerm(:call, inop, [
			PatTerm(:call, outop, [PatVar(:a), PatVar(:c)], __module__),
			PatTerm(:call, outop, [PatVar(:b), PatVar(:c)], __module__),
		], __module__))
end

# theory generation macros
export @commutativity
export @associativity
export @identity_left
export @identity_right
export @distrib_left
export @distrib_right
export @distrib
export @monoid
export @commutative_monoid
export @commutative_group

end


text:
using Sass

function compile(name)
    src = joinpath(@__DIR__, "scss", "$(name).scss")
    dst = joinpath(@__DIR__, "$(name).css")
    isfile(src) || error("$name not at $src")
    Sass.compile_file(src, dst)
end

compile("documenter")
compile("darkly")


text:
# Testing:
#
#  - computation of sufficient statistics
#  - distribution fitting (i.e. estimation)
#

using Distributions
using Base.Test

n0 = 100
N = 10^5

w = rand(n0)

# DiscreteUniform

x = rand(DiscreteUniform(10, 15), n0)
d = fit(DiscreteUniform, x)
@test isa(d, DiscreteUniform)
@test minimum(d) == minimum(x)
@test maximum(d) == maximum(x)

d = fit(DiscreteUniform, rand(DiscreteUniform(10, 15), N))
@test minimum(d) == 10
@test maximum(d) == 15

# Bernoulli

x = rand(Bernoulli(0.7), n0)

ss = suffstats(Bernoulli, x)
@test isa(ss, Distributions.BernoulliStats)
@test ss.cnt0 == n0 - countnz(x)
@test ss.cnt1 == countnz(x)

ss = suffstats(Bernoulli, x, w)
@test isa(ss, Distributions.BernoulliStats)
@test_approx_eq ss.cnt0 sum(w[x .== 0])
@test_approx_eq ss.cnt1 sum(w[x .== 1])

d = fit(Bernoulli, x)
p = countnz(x) / n0
@test isa(d, Bernoulli)
@test_approx_eq mean(d) p

d = fit(Bernoulli, x, w)
p = sum(w[x .== 1]) / sum(w)
@test isa(d, Bernoulli)
@test_approx_eq mean(d) p

d = fit(Bernoulli, rand(Bernoulli(0.7), N))
@test isa(d, Bernoulli)
@test_approx_eq_eps mean(d) 0.7 0.01

# Beta

d = fit(Beta, rand(Beta(1.3, 3.7), N))
@test isa(d, Beta)
@test_approx_eq_eps d.α 1.3 0.1
@test_approx_eq_eps d.β 3.7 0.1

# Binomial

x = rand(Binomial(100, 0.3), n0)

ss = suffstats(Binomial, (100, x))
@test isa(ss, Distributions.BinomialStats)
@test_approx_eq ss.ns sum(x)
@test ss.ne == n0
@test ss.n == 100

ss = suffstats(Binomial, (100, x), w)
@test isa(ss, Distributions.BinomialStats)
@test_approx_eq ss.ns dot(Float64[xx for xx in x], w)
@test_approx_eq ss.ne sum(w)
@test ss.n == 100

d = fit(Binomial, (100, x))
@test isa(d, Binomial)
@test ntrials(d) == 100
@test_approx_eq succprob(d) sum(x) / (n0 * 100)

d = fit(Binomial, (100, x), w)
@test isa(d, Binomial)
@test ntrials(d) == 100
@test_approx_eq succprob(d) dot(x, w) / (sum(w) * 100)

d = fit(Binomial, 100, rand(Binomial(100, 0.3), N))
@test isa(d, Binomial)
@test ntrials(d) == 100
@test_approx_eq_eps succprob(d) 0.3 0.01

# Categorical

p = [0.2, 0.5, 0.3]
x = rand(Categorical(p), n0)

ss = suffstats(Categorical, (3, x))
h = Float64[countnz(x .== i) for i = 1 : 3]
@test isa(ss, Distributions.CategoricalStats)
@test_approx_eq ss.h h

d = fit(Categorical, (3, x))
@test isa(d, Categorical)
@test d.K == 3
@test_approx_eq probs(d) h / sum(h)

d2 = fit(Categorical, x)
@test isa(d2, Categorical)
@test probs(d2) == probs(d)

ss = suffstats(Categorical, (3, x), w)
h = Float64[sum(w[x .== i]) for i = 1 : 3]
@test isa(ss, Distributions.CategoricalStats)
@test_approx_eq ss.h h

d = fit(Categorical, (3, x), w)
@test isa(d, Categorical)
@test_approx_eq probs(d) h / sum(h)

d = fit(Categorical, suffstats(Categorical, 3, x, w))
@test isa(d, Categorical)
@test_approx_eq probs(d) (h / sum(h))

d = fit(Categorical, rand(Categorical(p), N))
@test isa(d, Categorical)
@test_approx_eq_eps probs(d) p 0.01

# Exponential

x = rand(Exponential(0.5), n0)

ss = suffstats(Exponential, x)
@test isa(ss, Distributions.ExponentialStats)
@test_approx_eq ss.sx sum(x)
@test ss.sw == n0

ss = suffstats(Exponential, x, w)
@test isa(ss, Distributions.ExponentialStats)
@test_approx_eq ss.sx dot(x, w)
@test ss.sw == sum(w)

d = fit(Exponential, x)
@test isa(d, Exponential)
@test_approx_eq scale(d) mean(x)

d = fit(Exponential, x, w)
@test isa(d, Exponential)
@test_approx_eq scale(d) dot(x, w) / sum(w)

d = fit(Exponential, rand(Exponential(0.5), N))
@test isa(d, Exponential)
@test_approx_eq_eps scale(d) 0.5 0.01

# Normal

μ = 11.3
σ = 3.2

x = rand(Normal(μ, σ), n0)

ss = suffstats(Normal, x)
@test isa(ss, Distributions.NormalStats)
@test_approx_eq ss.s sum(x)
@test_approx_eq ss.m mean(x)
@test_approx_eq ss.s2 sum((x .- ss.m).^2)
@test_approx_eq ss.tw n0

ss = suffstats(Normal, x, w)
@test isa(ss, Distributions.NormalStats)
@test_approx_eq ss.s dot(x, w)
@test_approx_eq ss.m dot(x, w) / sum(w)
@test_approx_eq ss.s2 dot((x .- ss.m).^2, w)
@test_approx_eq ss.tw sum(w)

d = fit(Normal, x)
@test isa(d, Normal)
@test_approx_eq d.μ mean(x)
@test_approx_eq d.σ sqrt(mean((x .- d.μ).^2))

d = fit(Normal, x, w)
@test isa(d, Normal)
@test_approx_eq d.μ dot(x, w) / sum(w)
@test_approx_eq d.σ sqrt(dot((x .- d.μ).^2, w) / sum(w))

d = fit(Normal, rand(Normal(μ, σ), N))
@test isa(d, Normal)
@test_approx_eq_eps d.μ μ 0.1
@test_approx_eq_eps d.σ σ 0.1

import Distributions.NormalKnownMu, Distributions.NormalKnownSigma

ss = suffstats(NormalKnownMu(μ), x)
@test isa(ss, Distributions.NormalKnownMuStats)
@test ss.μ == μ
@test_approx_eq ss.s2 sum((x .- μ).^2)
@test_approx_eq ss.tw n0

ss = suffstats(NormalKnownMu(μ), x, w)
@test isa(ss, Distributions.NormalKnownMuStats)
@test ss.μ == μ
@test_approx_eq ss.s2 dot((x .- μ).^2, w)
@test_approx_eq ss.tw sum(w)

d = fit_mle(Normal, x; mu=μ)
@test isa(d, Normal)
@test d.μ == μ
@test_approx_eq d.σ sqrt(mean((x .- d.μ).^2))

d = fit_mle(Normal, x, w; mu=μ)
@test isa(d, Normal)
@test d.μ == μ
@test_approx_eq d.σ sqrt(dot((x .- d.μ).^2, w) / sum(w))

ss = suffstats(NormalKnownSigma(σ), x)
@test isa(ss, Distributions.NormalKnownSigmaStats)
@test ss.σ == σ
@test_approx_eq ss.sx sum(x)
@test_approx_eq ss.tw n0

ss = suffstats(NormalKnownSigma(σ), x, w)
@test isa(ss, Distributions.NormalKnownSigmaStats)
@test ss.σ == σ
@test_approx_eq ss.sx dot(x, w)
@test_approx_eq ss.tw sum(w)

d = fit_mle(Normal, x; sigma=σ)
@test isa(d, Normal)
@test d.σ == σ
@test_approx_eq d.μ mean(x)

d = fit_mle(Normal, x, w; sigma=σ)
@test isa(d, Normal)
@test d.σ == σ
@test_approx_eq d.μ dot(x, w) / sum(w)

# Uniform

x = rand(Uniform(1.2, 5.8), n0)
d = fit(Uniform, x)
@test isa(d, Uniform)
@test 1.2 <= minimum(d) <= maximum(d) <= 5.8
@test minimum(d) == minimum(x)
@test maximum(d) == maximum(x)

d = fit(Uniform, rand(Uniform(1.2, 5.8), N))
@test 1.2 <= minimum(d) <= maximum(d) <= 5.8
@test_approx_eq_eps minimum(d) 1.2 0.02
@test_approx_eq_eps maximum(d) 5.8 0.02

# Gamma

x = rand(Gamma(3.9, 2.1), n0)

ss = suffstats(Gamma, x)
@test isa(ss, Distributions.GammaStats)
@test_approx_eq ss.sx sum(x)
@test_approx_eq ss.slogx sum(log(x))
@test_approx_eq ss.tw n0

ss = suffstats(Gamma, x, w)
@test isa(ss, Distributions.GammaStats)
@test_approx_eq ss.sx dot(x, w)
@test_approx_eq ss.slogx dot(log(x), w)
@test_approx_eq ss.tw sum(w)

d = fit(Gamma, rand(Gamma(3.9, 2.1), N))
@test isa(d, Gamma)
@test_approx_eq_eps shape(d) 3.9 0.1
@test_approx_eq_eps scale(d) 2.1 0.2

# Geometric

x = rand(Geometric(0.3), n0)

ss = suffstats(Geometric, x)
@test isa(ss, Distributions.GeometricStats)
@test_approx_eq ss.sx sum(x)
@test_approx_eq ss.tw n0

ss = suffstats(Geometric, x, w)
@test isa(ss, Distributions.GeometricStats)
@test_approx_eq ss.sx dot(x, w)
@test_approx_eq ss.tw sum(w)

d = fit(Geometric, x)
@test isa(d, Geometric)
@test_approx_eq succprob(d) inv(1. + mean(x))

d = fit(Geometric, x, w)
@test isa(d, Geometric)
@test_approx_eq succprob(d) inv(1. + dot(x, w) / sum(w))

d = fit(Geometric, rand(Geometric(0.3), N))
@test isa(d, Geometric)
@test_approx_eq_eps succprob(d) 0.3 0.01

# Laplace

d = fit(Laplace, rand(Laplace(5.0, 3.0), N))
@test isa(d, Laplace)
@test_approx_eq_eps location(d) 5.0 0.1
@test_approx_eq_eps scale(d) 3.0 0.2

# Poisson

x = rand(Poisson(8.2), n0)

ss = suffstats(Poisson, x)
@test isa(ss, Distributions.PoissonStats)
@test_approx_eq ss.sx sum(x)
@test_approx_eq ss.tw n0

ss = suffstats(Poisson, x, w)
@test isa(ss, Distributions.PoissonStats)
@test_approx_eq ss.sx dot(x, w)
@test_approx_eq ss.tw sum(w)

d = fit(Poisson, x)
@test isa(d, Poisson)
@test_approx_eq mean(d) mean(x)

d = fit(Poisson, x, w)
@test isa(d, Poisson)
@test_approx_eq mean(d) dot(Float64[xx for xx in x], w) / sum(w)

d = fit(Poisson, rand(Poisson(8.2), N))
@test isa(d, Poisson)
@test_approx_eq_eps mean(d) 8.2 0.2


text:
# somehow get all AT detections
# you get some kind of list/dict of AT detections per pose
for p in poses
  # p is symbol such as :x45
  for detection in dets[p]
    # detection has (id, homography)
    lmsyms = ls(fg)[2]
    atsym = Symbol("l$(detection.id)")
    if !(atsym in lmsyms)
      addVariable!(fg, atsym, Point2(labels=["LANDMARK"; "APRILTAG"]))
    end

    # reverse homography to range, bearing
    pp2 = Point2Point2BearingRange(...)
    addFactor!(fg, [p; atsym], pp2)
  end
end