text:
#### Math with images ####

(+)(img::AbstractImageDirect{Bool}, n::Bool) = img .+ n
(+)(img::AbstractImageDirect, n::Number) = img .+ n
(+)(n::Bool, img::AbstractImageDirect{Bool}) = n .+ img
(+)(n::Number, img::AbstractImageDirect) = n .+ img
(.+)(img::AbstractImageDirect, n::Number) = limadj(copy(img, data(img).+n), limplus(limits(img), n))
(.+)(n::Number, img::AbstractImageDirect) = limadj(copy(img, data(img).+n), limplus(limits(img), n))
(+)(img::AbstractImageDirect, A::BitArray) = limadj(copy(img, data(img)+A), limplus(limits(img), Bool))
if isdefined(:UniformScaling)
    (+){Timg,TA<:Number}(img::AbstractImageDirect{Timg,2}, A::UniformScaling{TA}) = limadj(copy(img, data(img)+A), limplus(limits(img), A.λ))
    (-){Timg,TA<:Number}(img::AbstractImageDirect{Timg,2}, A::UniformScaling{TA}) = limadj(copy(img, data(img)-A), limminus(limits(img), A.λ))
end
(+)(img::AbstractImageDirect, A::AbstractArray) = limadj(copy(img, data(img)+data(A)), limplus(limits(img), limits(A)))
(.+)(img::AbstractImageDirect, A::BitArray) = limadj(copy(img, data(img).+A), limplus(limits(img), Bool))
(.+)(img::AbstractImageDirect, A::AbstractArray) = limadj(copy(img, data(img).+data(A)), limplus(limits(img), limits(A)))
(-)(img::AbstractImageDirect{Bool}, n::Bool) = img .- n
(-)(img::AbstractImageDirect, n::Number) = img .- n
(.-)(img::AbstractImageDirect, n::Number) = limadj(copy(img, data(img).-n), limminus(limits(img), n))
(-)(n::Bool, img::AbstractImageDirect{Bool}) = n .- img
(-)(n::Number, img::AbstractImageDirect) = n .- img
(.-)(n::Number, img::AbstractImageDirect) = limadj(copy(img, n.-data(img)), limminus(n, limits(img)))
(-)(img::AbstractImageDirect, A::BitArray) = limadj(copy(img, data(img)-A), limminus(limits(img), Bool))
(-){T}(img::AbstractImageDirect{T,2}, A::Diagonal) = limadj(copy(img, data(img)-A), limminus(limits(img), limits(A))) # fixes an ambiguity warning
(-)(img::AbstractImageDirect, A::AbstractArray) = limadj(copy(img, data(img)-data(A)), limminus(limits(img), limits(A)))
# (-)(A::AbstractArray, img::AbstractImageDirect) = limadj(copy(img, data(A) - data(img)), limminus(limits(A), limits(img)))
(-)(img::AbstractImageDirect) = limadj(copy(img, -data(img)), limtimes(limits(img), -1))
(.-)(img::AbstractImageDirect, A::BitArray) = limadj(copy(img, data(img).-A), limminus(limits(img), Bool))
(.-)(img::AbstractImageDirect, A::AbstractArray) = limadj(copy(img, data(img).-data(A)), limminus(limits(img), limits(A)))
(*)(img::AbstractImageDirect, n::Number) = (.*)(img, n)
(*)(n::Number, img::AbstractImageDirect) = (.*)(n, img)
(.*)(img::AbstractImageDirect, n::Number) = limadj(copy(img, data(img).*n), limtimes(limits(img), n))
(.*)(n::Number, img::AbstractImageDirect) = limadj(copy(img, data(img).*n), limtimes(limits(img), n))
(/)(img::AbstractImageDirect, n::Number) = limadj(copy(img, data(img)/n), limdivide(limits(img), n))
(.*)(img1::AbstractImageDirect, img2::AbstractImageDirect) = limadj(copy(img1, data(img1).*data(img2)), limtimes(limits(img1), limits(img2)))
(.*)(img::AbstractImageDirect, A::BitArray) = limadj(copy(img, data(img).*A), limtimes(limits(img), Bool))
(.*)(A::BitArray, img::AbstractImageDirect) = limadj(copy(img, data(img).*A), limtimes(limits(img), Bool))
(.*)(img::AbstractImageDirect{Bool}, A::BitArray) = limadj(copy(img, data(img).*A), (false,true))
(.*)(A::BitArray, img::AbstractImageDirect{Bool}) = limadj(copy(img, data(img).*A), (false,true))
(.*)(img::AbstractImageDirect, A::AbstractArray) = limadj(copy(img, data(img).*A), limtimes(limits(img), limits(A)))
(.*)(A::AbstractArray, img::AbstractImageDirect) = limadj(copy(img, data(img).*A), limtimes(limits(img), limits(A)))
(./)(img::AbstractImageDirect, A::BitArray) = limadj(copy(img, data(img)./A), limdivide(limits(img), Bool))  # needed to avoid ambiguity warning
(./)(img1::AbstractImageDirect, img2::AbstractImageDirect) = limadj(copy(img1, data(img1)./data(img2)), limdivide(limits(img1), limits(img2)))
(./)(img::AbstractImageDirect, A::AbstractArray) = limadj(copy(img, data(img)./A), limdivide(limits(img), limits(A)))
# (./)(A::AbstractArray, img::AbstractImageDirect) = limadj(copy(img, A./data(img))
(.^)(img::AbstractImageDirect, p::Number) = limadj(copy(img, data(img).^p), limpower(limits(img), p))
sqrt(img::AbstractImageDirect) = limadj(copy(img, sqrt(data(img))), limpower(limits(img), 0.5))
atan2(img1::AbstractImageDirect, img2::AbstractImageDirect) = (img = copy(img1, atan2(data(img1), data(img2))); img["limits"] = (-float(pi),float(pi)); img)
hypot(img1::AbstractImageDirect, img2::AbstractImageDirect) = limadj(copy(img1, hypot(data(img1), data(img2))), limpower(limplus(limpower(limits(img1),2),limpower(limits(img2),2)),0.5))

function limadj(img::AbstractImageDirect, newlim)
    if haskey(img, "limits")
        T = eltype(img)
        if method_exists(typemin, (T,)) && method_exists(typemax, (T,))
            img["limits"] = (minclamp(typemin(T), newlim[1]), maxclamp(typemax(T), newlim[2]))
        else
            img["limits"] = (convert(T, newlim[1]), convert(T, newlim[2]))
        end
    end
    img
end

limplus(a::Tuple, n::Number) = a[1]+n, a[2]+n
limplus(a::Tuple, ::Type{Bool}) = a[1], a[2]+1
limplus(a::Tuple, b::Tuple) = a[1]+b[1], a[2]+b[2]
limminus(a::Tuple, n::Number) = a[1]-n, a[2]-n
limminus(n::Number, a::Tuple) = n-a[2], n-a[1]
limminus(a::Tuple, ::Type{Bool}) = a[1]-1, a[2]
limminus(a::Tuple, b::Tuple) = a[1]-b[2], a[2]-b[1]
limtimes(a::Tuple, n::Number) = n > 0 ? (a[1]*n, a[2]*n) : (a[2]*n, a[1]*n)
limtimes(a::Tuple, ::Type{Bool}) = min(a[1], false*a[1]), max(a[2], false*a[2])
limtimes(a::Tuple, b::Tuple) = min(a[1]*b[1],a[1]*b[2],a[2]*b[1],a[2]*b[2]), max(a[1]*b[1],a[1]*b[2],a[2]*b[1],a[2]*b[2])
limdivide(a::Tuple, n::Number) = n > 0 ? (a[1]/n, a[2]/n) : (a[2]/n, a[1]/n)
limdivide(a::Tuple, ::Type{Bool}) = min(a[1],a[1]/0), max(a[2],a[2]/0)
limdivide(a::Tuple, b::Tuple) = min(a[1]/b[1],a[1]/b[2],a[2]/b[1],a[2]/b[2]), max(a[1]/b[1],a[1]/b[2],a[2]/b[1],a[2]/b[2])
limpower(a::Tuple, p::Number) = (l = a[1]^p; u = a[2]^p; (min(l,u), max(l,u)))

function sum(img::AbstractImageDirect, region::Union(AbstractVector,Tuple,Integer))
    f = prod(size(img)[[region...]])
    imgs = sum(data(img), region)
    l = map(x->convert(eltype(imgs),x), limits(img))
    out = limadj(copy(img, imgs), (f*l[1], f*l[2]))
    if in(colordim(img), region)
        out["colorspace"] = "Unknown"
    end
    out
end

# Logical operations
(.<)(img::AbstractImageDirect, n::Number) = data(img) .< n
(.>)(img::AbstractImageDirect, n::Number) = data(img) .> n
(.<)(img::AbstractImageDirect{Bool}, A::AbstractArray{Bool}) = data(img) .< A
(.<)(img::AbstractImageDirect, A::AbstractArray) = data(img) .< A
(.>)(img::AbstractImageDirect, A::AbstractArray) = data(img) .> A
(.==)(img::AbstractImageDirect, n::Number) = data(img) .== n
(.==)(img::AbstractImageDirect{Bool}, A::AbstractArray{Bool}) = data(img) .== A
(.==)(img::AbstractImageDirect, A::AbstractArray) = data(img) .== A

#### Overlay AbstractArray implementation ####
length(o::Overlay) = isempty(o.channels) ? 0 : length(o.channels[1])
size(o::Overlay) = isempty(o.channels) ? (0,) : size(o.channels[1])
size(o::Overlay, d::Integer) = isempty(o.channels) ? 0 : size(o.channels[1],d)
eltype(o::Overlay) = RGB
nchannels(o::Overlay) = length(o.channels)

similar(o::Overlay) = Array(RGB, size(o))
similar(o::Overlay, ::NTuple{0}) = Array(RGB, size(o))
similar{T}(o::Overlay, ::Type{T}) = Array(T, size(o))
similar{T}(o::Overlay, ::Type{T}, sz::Int64) = Array(T, sz)
similar{T}(o::Overlay, ::Type{T}, sz::Int64...) = Array(T, sz)
similar{T}(o::Overlay, ::Type{T}, sz) = Array(T, sz)

function getindex(o::Overlay, indexes::Integer...)
    rgb = RGB(0,0,0)
    for i = 1:nchannels(o)
        if o.visible[i]
            rgb += scale(o.scalei[i], o.channels[i][indexes...])*o.colors[i]
        end
    end
    clip(rgb)
end

# Fix ambiguity warning
getindex(o::Overlay, i::Real) = getindex_overlay(o, i)
getindex(o::Overlay, indexes::Union(Real,AbstractVector)...) = getindex_overlay(o, indexes...)

function getindex_overlay(o::Overlay, indexes::Union(Real,AbstractVector)...)
    len = [length(i) for i in indexes]
    n = length(len)
    while len[n] == 1
        pop!(len)
        n -= 1
    end
    rgb = fill(RGB(0,0,0), len...)
    for i = 1:nchannels(o)
        if o.visible[i]
            tmp = scale(o.scalei[i], o.channels[i][indexes...])
            accumulate!(rgb, tmp, o.colors[i])
        end
    end
    clip!(rgb)
end

getindex(o::Overlay, indexes::(AbstractVector...)) = getindex(o, indexes...)

setindex!(o::Overlay, x, index::Real) = error("Overlays are read-only")
setindex!(o::Overlay, x, indexes...) = error("Overlays are read-only")

# Identical to getindex except it saves a call to each array's getindex
function convert(::Type{Array{RGB}}, o::Overlay)
    rgb = fill(RGB(0,0,0), size(o))
    for i = 1:length(o.channels)
        if o.visible[i]
            tmp = scale(o.scalei[i], o.channels[i])
            accumulate!(rgb, tmp, o.colors[i])
        end
    end
    clip!(rgb)
end

function accumulate!(rgb::Array{RGB}, A::Array{Float64}, color::RGB)
    for j = 1:length(rgb)
        rgb[j] += A[j]*color
    end
end

for N = 1:4
    @eval begin
        function accumulate!{T}(rgb::Array{RGB,$N}, A::AbstractArray{T,$N}, color::RGB, scalei::ScaleInfo)
            k = 0
            @inbounds @nloops $N i A begin
                rgb[k+=1] += scale(scalei, (@nref $N A i))*color
            end
        end
    end
end

(*)(f::FloatingPoint, c::RGB) = RGB(f*c.r, f*c.g, f*c.b)
(*)(f::Uint8, c::RGB) = (f/255)*c
(/)(c::RGB, f::Real) = (1.0/f)*c
(.*)(f::AbstractArray, c::RGB) = [x*c for x in f]
(+)(a::RGB, b::RGB) = RGB(a.r+b.r, a.g+b.g, a.b+b.b)
convert(::Type{Uint32}, c::ColorValue) = convert(RGB24, c).color

#### Converting images to uint32 color ####
# This is the crucial operation in image display

function uint32color(img, args...)
    sz = size(img)
    ssz = sz[coords_spatial(img)]
    buf = Array(Uint32, ssz)
    uint32color!(buf, img, args...)
    buf
end

function uint32color!(buf::Array{Uint32}, img, args...)
    cdim = colordim(img)
    if cdim != 0
        sz = size(img)
        if size(buf) != tuple(sz[1:cdim-1]..., sz[cdim+1:end]...)
            error("Size mismatch")
        end
        _uint32color!(buf, img, colorspace(img), cdim, args...)
    else
        if size(buf) != size(img)
            error("Size mismatch")
        end
        if colorspace(img) == "RGB24"
            _uint32color_rgb24!(buf, img, args...)
        else
            _uint32color_gray!(buf, img, args...)
        end
    end
    buf
end

# Indexed images (colormaps)
function uint32color!(buf::Array{Uint32}, img::AbstractImageIndexed)
    dat = data(img)
    cmap = img.cmap
    for i = 1:length(buf)
        buf[i] = convert(Uint32, cmap[dat[i]])
    end
    buf
end

# ColorValue arrays
for N = 1:4
    @eval begin
        function uint32color!{C<:ColorValue}(buf::Array{Uint32}, img::AbstractArray{C,$N})
            if size(buf) != size(img)
                error("Size mismatch")
            end
            dat = data(img)
            k = 0
            @inbounds @nloops $N i dat begin
                val = @nref $N dat i
                buf[k+=1] = convert(RGB24, val)
            end
            buf
        end
    end
end

# Grayscale arrays
for N = 1:4
    @eval begin
        function _uint32color_gray!{T}(buf::Array{Uint32}, A::AbstractArray{T,$N}, scalei::ScaleInfo = scaleinfo(Uint8, A))
            scalei_t = take(scalei, A)
            Adat = data(A)
            k = 0
            @inbounds @nloops $N i A begin
                val = @nref $N Adat i
                gr = scale(scalei_t, val)
                buf[k+=1] = rgb24(gr, gr, gr)
            end
            buf
        end
    end
end

# RGB24 arrays (just a copy operation)
for N = 1:4
    @eval begin
        function _uint32color_rgb24!{T}(buf::Array{Uint32}, A::AbstractArray{T,$N})
            k = 0
            Adat = data(A)
            @inbounds @nloops $N i Adat begin
                buf[k+=1] = @nref $N Adat i
            end
            buf
        end
    end
end

# Arrays where one dimension encodes color or transparency
for N = 1:5
    @eval begin
        function _uint32color!{T}(buf::Array{Uint32}, A::AbstractArray{T,$N}, cs::String, cdim::Int, scalei::ScaleInfo = scaleinfo(Uint8, A))
            scalei_t = take(scalei, A)
            k = 0
            Adat = data(A)
            # Because cdim is not known at compile-time, we have to prepare size and step variables
            cstep = zeros(Int, $N)
            cstep[cdim] = 1
            szA = [size(Adat, d) for d = 1:$N]
            szA[cdim] = 1  # don't loop over color dimension
            @nextract $N szA szA
            @nextract $N cstep cstep
            if cs == "GrayAlpha"
                @inbounds @nloops $N i d->(1:szA_d) begin
                    gr = @nref $N Adat i
                    alpha = @nref $N Adat d->(i_d+cstep_d)
                    buf[k+=1] = argb32(scalei_t, alpha, gr, gr, gr)
                end
            elseif cs == "RGB"
                @inbounds @nloops $N i d->(1:szA_d) begin
                    r = @nref $N Adat i
                    g = @nref $N Adat d->(i_d+cstep_d)
                    b = @nref $N Adat d->(i_d+2cstep_d)
                    buf[k+=1] = rgb24(scalei_t, r, g, b)
                end
            elseif cs == "ARGB"
                @inbounds @nloops $N i d->(1:szA_d) begin
                    a = @nref $N Adat i
                    r = @nref $N Adat d->(i_d+cstep_d)
                    g = @nref $N Adat d->(i_d+2cstep_d)
                    b = @nref $N Adat d->(i_d+3cstep_d)
                    buf[k+=1] = argb32(scalei_t, a, r, g, b)
                end
            elseif cs == "RGBA"
                @inbounds @nloops $N i d->(1:szA_d) begin
                    r = @nref $N Adat i
                    g = @nref $N Adat d->(i_d+cstep_d)
                    b = @nref $N Adat d->(i_d+2cstep_d)
                    a = @nref $N Adat d->(i_d+3cstep_d)
                    buf[k+=1] = argb32(scalei_t, a, r, g, b)
                end
            else
                error("colorspace ", cs, " not yet supported")
            end
        end
    end
end

# Overlays
_wrap(A::Array, B::AbstractArray) = B
_wrap(S::SubArray, B::AbstractArray) = sub(B, S.indexes...)
function uint32color!{O<:Overlay,N,IT<:(RangeIndex...,)}(buf::Array{Uint32}, img::Union(Overlay, Image{RGB,N,O}, SubArray{RGB,N,O,IT}, Image{RGB,N,SubArray{RGB,N,O,IT}}))
    if size(buf) != size(img)
        error("Size mismatch")
    end
    rgb = fill(RGB(0,0,0), size(buf))
    dat = data(img)
    p = parent(dat)
    for i = 1:nchannels(p)
        if p.visible[i]
            accumulate!(rgb, _wrap(dat, data(p.channels[i])), p.colors[i], p.scalei[i])
        end
    end
    clip!(rgb)
    for i = 1:length(buf)
        buf[i] = convert(RGB24, rgb[i])
    end
    buf
end

# Signed grayscale arrays colorized magenta (positive) and green (negative
for N = 1:4
    @eval begin
        function _uint32color_gray!{T}(buf::Array{Uint32}, A::AbstractArray{T,$N}, scalei::ScaleSigned)
            k = 0
            Adat = data(A)
            @inbounds @nloops $N i Adat begin
                val = @nref $N Adat i
                gr = scale(scalei, val)
                if isfinite(gr)
                    if gr >= 0
                        gr8 = iround(Uint8, 255.0*gr)
                        buf[k+=1] = rgb24(gr8, 0x00, gr8)
                    else
                        gr8 = iround(Uint8, -255.0*gr)
                        buf[k+=1] = rgb24(0x00, gr8, 0x00)
                    end
                else
                    buf[k+=1] = zero(Uint32)
                end
            end
            buf
        end
    end
end

function rgb24(r::Uint8, g::Uint8, b::Uint8)
    ret::Uint32
    ret = convert(Uint32,r)<<16 | convert(Uint32,g)<<8 | convert(Uint32,b)
end

function argb32(a::Uint8, r::Uint8, g::Uint8, b::Uint8)
    ret::Uint32
    ret = convert(Uint32,a)<<24 | convert(Uint32,r)<<16 | convert(Uint32,g)<<8 | convert(Uint32,b)
end

function rgb24{T}(scalei::ScaleInfo{Uint8}, r::T, g::T, b::T)
    ret::Uint32
    ret = convert(Uint32,scale(scalei,r))<<16 | convert(Uint32,scale(scalei,g))<<8 | convert(Uint32,scale(scalei,b))
end

function argb32{T}(scalei::ScaleInfo{Uint8}, a::T, r::T, g::T, b::T)
    ret::Uint32
    ret = convert(Uint32,scale(scalei,a))<<24 | convert(Uint32,scale(scalei,r))<<16 | convert(Uint32,scale(scalei,g))<<8 | convert(Uint32,scale(scalei,b))
end

#### Color palettes ####

function lut(pal::Vector, a)
    out = similar(a, eltype(pal))
    n = length(pal)
    for i=1:length(a)
        out[i] = pal[clamp(a[i], 1, n)]
    end
    out
end

function indexedcolor(data, pal)
    mn = minimum(data); mx = maximum(data)
    indexedcolor(data, pal, mx-mn, (mx+mn)/2)
end

function indexedcolor(data, pal, w, l)
    n = length(pal)-1
    if n == 0
        return fill(pal[1], size(data))
    end
    w_min = l - w/2
    scale = w==0 ? 1 : w/n
    lut(pal, iround((data - w_min)./scale) + 1)
end

const palette_gray32  = [0xff000000,0xff080808,0xff101010,0xff181818,0xff202020,0xff292929,0xff313131,0xff393939,
                         0xff414141,0xff4a4a4a,0xff525252,0xff5a5a5a,0xff626262,0xff6a6a6a,0xff737373,0xff7b7b7b,
                         0xff838383,0xff8b8b8b,0xff949494,0xff9c9c9c,0xffa4a4a4,0xffacacac,0xffb4b4b4,0xffbdbdbd,
                         0xffc5c5c5,0xffcdcdcd,0xffd5d5d5,0xffdedede,0xffe6e6e6,0xffeeeeee,0xfff6f6f6,0xffffffff]

const palette_gray64  = [0xff000000,0xff040404,0xff080808,0xff0c0c0c,0xff101010,0xff141414,0xff181818,0xff1c1c1c,
                         0xff202020,0xff242424,0xff282828,0xff2c2c2c,0xff303030,0xff343434,0xff383838,0xff3c3c3c,
                         0xff404040,0xff444444,0xff484848,0xff4c4c4c,0xff505050,0xff555555,0xff595959,0xff5d5d5d,
                         0xff616161,0xff656565,0xff696969,0xff6d6d6d,0xff717171,0xff757575,0xff797979,0xff7d7d7d,
                         0xff818181,0xff858585,0xff898989,0xff8d8d8d,0xff919191,0xff959595,0xff999999,0xff9d9d9d,
                         0xffa1a1a1,0xffa5a5a5,0xffaaaaaa,0xffaeaeae,0xffb2b2b2,0xffb6b6b6,0xffbababa,0xffbebebe,
                         0xffc2c2c2,0xffc6c6c6,0xffcacaca,0xffcecece,0xffd2d2d2,0xffd6d6d6,0xffdadada,0xffdedede,
                         0xffe2e2e2,0xffe6e6e6,0xffeaeaea,0xffeeeeee,0xfff2f2f2,0xfff6f6f6,0xfffafafa,0xffffffff]

const palette_fire    = [0xff5a5a5a,0xff636058,0xff6c6757,0xff756e56,0xff7e7455,0xff877b54,0xff908253,0xff998851,
                         0xffa28f50,0xffab964f,0xffb49c4e,0xffbda34d,0xffc6aa4c,0xffcfb04a,0xffd8b749,0xffe1be48,
                         0xffeac447,0xfff3cb46,0xfffdd245,0xfffccc42,0xfffcc640,0xfffcc03d,0xfffbba3b,0xfffbb438,
                         0xfffbae36,0xfffaa833,0xfffaa231,0xfffa9c2e,0xfffa962c,0xfff99029,0xfff98a27,0xfff98424,
                         0xfff87e22,0xfff8781f,0xfff8721d,0xfff76c1a,0xfff76618,0xfff76015,0xfff75a13,0xfff65410,
                         0xfff64e0e,0xfff6480b,0xfff54209,0xfff53c06,0xfff53604,0xfff53102,0xffe72e01,0xffd92b01,
                         0xffcc2801,0xffbe2601,0xffb02301,0xffa32001,0xff951d01,0xff881b01,0xff7a1801,0xff6c1500,
                         0xff5f1300,0xff511000,0xff440d00,0xff360a00,0xff280800,0xff1b0500,0xff0d0200,0xff000000]

const palette_rainbow = [0xff0e46e9,0xff0d58ea,0xff0c6bec,0xff0c7eee,0xff0b91f0,0xff0ba4f1,0xff0ab7f3,0xff0acaf5,
                         0xff09ddf7,0xff09f0f9,0xff06efbd,0xff04ef81,0xff02ee45,0xff00ee0a,0xff1cee08,0xff38ee07,
                         0xff54ee06,0xff70ee05,0xff8dee04,0xffa9ee03,0xffc5ee02,0xffe1ee01,0xfffeee00,0xfffcd401,
                         0xfffbba02,0xfffaa104,0xfff98705,0xfff86e07,0xfff75408,0xfff63b0a,0xfff5210b,0xfff4080d]

redval(p)   = (p>>>16)&0xff
greenval(p) = (p>>>8)&0xff
blueval(p)  = p&0xff
alphaval(p)   = (p>>>24)&0xff

function imadjustintensity{T}(img::AbstractArray{T}, range)
    assert_scalar_color(img)
    if length(range) == 0
        range = [minimum(img) maximum(img)]
    elseif length(range) == 1
        error("incorrect range")
    end
    tmp = (img .- range[1])/(range[2] - range[1])
    tmp[tmp .> 1] = 1
    tmp[tmp .< 0] = 0
    out = tmp
end

function red{T<:ColorValue}(img::AbstractArray{T})
    out = Array(Float64, size(img))
    for i = 1:length(img)
        out[i] = convert(RGB, img[i]).r
    end
    out
end

function red(img::AbstractArray)
    if colorspace(img) == "RGB" || colorspace(img) == "RGBA"
        ret = sliceim(img, "color", 1)
    else
        error("Not yet implemented")
    end
end

function green{T<:ColorValue}(img::AbstractArray{T})
    out = Array(Float64, size(img))
    for i = 1:length(img)
        out[i] = convert(RGB, img[i]).g
    end
    out
end

function green(img::AbstractArray)
    if colorspace(img) == "RGB" || colorspace(img) == "RGBA"
        ret = sliceim(img, "color", 2)
    else
        error("Not yet implemented")
    end
end

function blue{T<:ColorValue}(img::AbstractArray{T})
    out = Array(Float64, size(img))
    for i = 1:length(img)
        out[i] = convert(RGB, img[i]).b
    end
    out
end

function blue(img::AbstractArray)
    if colorspace(img) == "RGB" || colorspace(img) == "RGBA"
        ret = sliceim(img, "color", 3)
    else
        error("Not yet implemented")
    end
end

for N = 1:4
    N1 = N-1
    @eval begin
        function rgb2gray{T}(img::AbstractArray{T,$N})
            cs = colorspace(img)
            if cs == "Gray"
                return img
            end
            if cs != "RGB"
                error("Color space of image is $cs, not RGB")
            end
            strds = strides(img)
            cd = colordim(img)
            cstrd = strds[cd]
            sz = [size(img,d) for d=1:$N]
            sz[cd] = 1
            szgray = sz[setdiff(1:$N,cd)]
            out = Array(T, szgray...)::Array{T,$N1}
            wr, wg, wb = 0.299, 0.587, 0.114    # Rec 601 luma conversion
            dat = data(img)
            indx = 0
            @nexprs $N d->(strd_d = strds[d])
            @nexprs $N d->(sz_d = sz[d])
            @nexprs $N d->(k_d = 1)
            @nloops $N i d->1:sz_d d->(k_{d-1}=k_d+strd_d*(i_d-1)) begin
                out[indx+=1] = truncround(T,wr*dat[k_0] + wg*dat[k_0+cstrd] + wb*dat[k_0+2cstrd])
            end
            p = copy(properties(img))
            p["colorspace"] = "Gray"
            p["colordim"] = 0
            Image(out, p)
        end
    end
end

# average filter
function imaverage(filter_size=[3,3])
    if length(filter_size) != 2
        error("wrong filter size")
    end
    m, n = filter_size
    if mod(m, 2) != 1 || mod(n, 2) != 1
        error("filter dimensions must be odd")
    end
    f = ones(Float64, m, n)/(m*n)
end

# laplacian filter kernel
function imlaplacian(diagonals::String="nodiagonals")
    if diagonals == "diagonals"
        return [1.0 1.0 1.0; 1.0 -8.0 1.0; 1.0 1.0 1.0]
    elseif diagonals == "nodiagonals"
        return [0.0 1.0 0.0; 1.0 -4.0 1.0; 0.0 1.0 0.0]
    else
        error("Expected \"diagnoals\" or \"nodiagonals\" or Number, got: \"$diagonals\"")
    end
end

# more general version
function imlaplacian(alpha::Number)
    lc = alpha/(1 + alpha)
    lb = (1 - alpha)/(1 + alpha)
    lm = -4/(1 + alpha)
    return [lc lb lc; lb lm lb; lc lb lc]
end

# 2D gaussian filter kernel
function gaussian2d(sigma::Number=0.5, filter_size=[])
    if length(filter_size) == 0
        # choose 'good' size
        m = 4*ceil(sigma)+1
        n = m
    elseif length(filter_size) != 2
        error("wrong filter size")
    else
        m, n = filter_size[1], filter_size[2]
    end
    if mod(m, 2) != 1 || mod(n, 2) != 1
        error("filter dimensions must be odd")
    end
    g = [exp(-(X.^2+Y.^2)/(2*sigma.^2)) for X=-floor(m/2):floor(m/2), Y=-floor(n/2):floor(n/2)]
    return g/sum(g)
end

# difference of gaussian
function imdog(sigma::Number=0.5)
    m = 4*ceil(sqrt(2)*sigma)+1
    return gaussian2d(sqrt(2)*sigma, [m m]) - gaussian2d(sigma, [m m])
end

# laplacian of gaussian
function imlog(sigma::Number=0.5)
    m = ceil(8.5sigma)
    m = m % 2 == 0 ? m + 1 : m
    return [(1/(2pi*sigma^4))*(2 - (x^2 + y^2)/sigma^2)*exp(-(x^2 + y^2)/(2sigma^2))
            for x=-floor(m/2):floor(m/2), y=-floor(m/2):floor(m/2)]
end

# Sum of squared differences
ssd{T}(A::AbstractArray{T}, B::AbstractArray{T}) = sum((data(A)-data(B)).^2)

# normalized by Array size
ssdn{T}(A::AbstractArray{T}, B::AbstractArray{T}) = ssd(A, B)/length(A)

# sum of absolute differences
sad{T}(A::AbstractArray{T}, B::AbstractArray{T}) = sum(abs(data(A)-data(B)))

# normalized by Array size
sadn{T}(A::AbstractArray{T}, B::AbstractArray{T}) = sad(A, B)/length(A)

# normalized cross correlation
function ncc{T}(A::AbstractArray{T}, B::AbstractArray{T})
    Am = (data(A).-mean(data(A)))[:]
    Bm = (data(B).-mean(data(B)))[:]
    return dot(Am,Bm)/(norm(Am)*norm(Bm))
end

# Array padding
function padarray{T,n}(img::AbstractArray{T,n}, prepad::Union(Vector{Int},Dims), postpad::Union(Vector{Int},Dims), border::String)
    I = Array(Vector{Int}, n)
    for d = 1:n
        M = size(img, d)
        I[d] = [(1 - prepad[d]):(M + postpad[d])]
        if border == "replicate"
            I[d] = min(max(I[d], 1), M)
        elseif border == "circular"
            I[d] = 1 .+ mod(I[d] .- 1, M)
        elseif border == "symmetric"
            I[d] = [1:M, M:-1:1][1 .+ mod(I[d] .- 1, 2 * M)]
        elseif border == "reflect"
            I[d] = [1:M, M-1:-1:2][1 .+ mod(I[d] .- 1, 2 * M - 2)]
        else
            error("unknown border condition")
        end
    end
    img[I...]::Array{T,n}
end

function padarray{T,n}(img::AbstractArray{T,n}, prepad::Union(Vector{Int},Dims), postpad::Union(Vector{Int},Dims), border::String, value)
    if border != "value"
        return padarray(img, prepad, postpad, border)
    end
    A = Array(T, ntuple(n, d->size(img,d)+prepad[d]+postpad[d]))
    fill!(A, value)
    I = Vector{Int}[1+prepad[d]:size(A,d)-postpad[d] for d = 1:n]
    A[I...] = img
    A::typeof(img)
end

padarray{T,n}(img::AbstractArray{T,n}, padding::Union(Vector{Int},Dims), border::String = "replicate") = padarray(img, padding, padding, border)
# Restrict the following to Number to avoid trouble when img is an Array{String}
padarray{T<:Number,n}(img::AbstractArray{T,n}, padding::Union(Vector{Int},Dims), value::T) = padarray(img, padding, padding, "value", value)

function padarray{T,n}(img::AbstractArray{T,n}, padding::Union(Vector{Int},Dims), border::String, direction::String)
    if direction == "both"
        return padarray(img, padding, padding, border)
    elseif direction == "pre"
        return padarray(img, padding, zeros(Int, n), border)
    elseif direction == "post"
        return padarray(img, zeros(Int, n), padding, border)
    end
end

function padarray{T<:Number,n}(img::AbstractArray{T,n}, padding::Vector{Int}, value::T, direction::String)
    if direction == "both"
        return padarray(img, padding, padding, "value", value)
    elseif direction == "pre"
        return padarray(img, padding, zeros(Int, n), "value", value)
    elseif direction == "post"
        return padarray(img, zeros(Int, n), padding, "value", value)
    end
end

function prep_kernel(img::AbstractArray, kern::AbstractArray)
    sc = coords_spatial(img)
    if ndims(kern) > length(sc)
        error("""The kernel has $(ndims(kern)) dimensions, more than the $(sdims(img)) spatial dimensions of img.
                 You'll need to set the dimensions and type of the kernel to be the same as the image.""")
    end
    sz = ones(Int, ndims(img))
    for i = 1:ndims(kern)
        sz[sc[i]] = size(kern,i)
    end
    reshape(kern, sz...)
end

###
### imfilter
###
imfilter(img, kern, border, value) = imfilter_inseparable(img, kern, border, value)
# Do not combine these with the previous using a default value (see the 2d specialization below)
imfilter(img, filter) = imfilter(img, filter, "replicate", 0)
imfilter(img, filter, border) = imfilter(img, filter, border, 0)

imfilter_inseparable{T,K,N,M}(img::AbstractArray{T,N}, kern::AbstractArray{K,M}, border::String, value) =
    imfilter_inseparable(img, prep_kernel(img, kern), border, value)

function imfilter_inseparable{T,K,N}(img::AbstractArray{T,N}, kern::AbstractArray{K,N}, border::String, value)
    prepad  = [div(size(kern,i)-1, 2) for i = 1:N]
    postpad = [div(size(kern,i),   2) for i = 1:N]
    A = padarray(img, prepad, postpad, border, convert(T, value))
    result = _imfilter!(Array(typeof(one(T)*one(K)), size(img)), A, data(kern))
    copy(img, result)
end

# Special case for 2d kernels: check for separability
function imfilter{T}(img::AbstractArray{T}, kern::AbstractMatrix, border::String, value)
    sc = coords_spatial(img)
    if length(sc) < 2
        imfilter_inseparable(img, kern, border, value)
    end
    SVD = svdfact(kern)
    U, S, Vt = SVD[:U], SVD[:S], SVD[:Vt]
    separable = true
    EPS = sqrt(eps(eltype(S)))
    for i = 2:length(S)
        separable &= (abs(S[i]) < EPS)
    end
    if !separable
        return imfilter_inseparable(img, kern, border, value)
    end
    s = S[1]
    u,v = U[:,1],Vt[1,:]
    ss = sqrt(s)
    sz1 = ones(Int, ndims(img)); sz1[sc[1]] = size(kern, 1)
    sz2 = ones(Int, ndims(img)); sz2[sc[2]] = size(kern, 2)
    tmp = imfilter_inseparable(data(img), reshape(u*ss, sz1...), border, value)
    copy(img, imfilter_inseparable(tmp, reshape(v*ss, sz2...), border, value))
end

for N = 1:5
    @eval begin
        function _imfilter!{T,K}(B, A::AbstractArray{T,$N}, kern::AbstractArray{K,$N})
            for i = 1:$N
                if size(B,i)+size(kern,i) > size(A,i)+1
                    throw(DimensionMismatch("Output dimensions $(size(B)) and kernel dimensions $(size(kern)) do not agree with size of padded input, $(size(A))"))
                end
            end
            @nloops $N i B begin
                tmp = zero(T)
                @inbounds begin
                    @nloops $N j kern begin
                        tmp += (@nref $N A d->(i_d+j_d-1))*(@nref $N kern j)
                    end
                    (@nref $N B i) = tmp
                end
            end
            B
        end
    end
end

###
### imfilter_fft
###
imfilter_fft(img, kern, border, value) = copy(img, imfilter_fft_inseparable(img, kern, border, value))
imfilter_fft(img, filter) = imfilter_fft(img, filter, "replicate", 0)
imfilter_fft(img, filter, border) = imfilter_fft(img, filter, border, 0)

imfilter_fft_inseparable{T,K,N,M}(img::AbstractArray{T,N}, kern::AbstractArray{K,M}, border::String, value) =
    imfilter_fft_inseparable(img, prep_kernel(img, kern), border, value)

function imfilter_fft_inseparable{T,K,N}(img::AbstractArray{T,N}, kern::AbstractArray{K,N}, border::String, value)
    prepad  = [div(size(kern,i)-1, 2) for i = 1:N]
    postpad = [div(size(kern,i),   2) for i = 1:N]
    fullpad = [nextprod([2,3], size(img,i) + prepad[i] + postpad[i]) - size(img, i) - prepad[i] for i = 1:N]
    A = padarray(img, prepad, fullpad, border, convert(T, value))
    krn = zeros(typeof(one(T)*one(K)), size(A))
    indexesK = ntuple(N, d->[size(krn,d)-prepad[d]+1:size(krn,d),1:size(kern,d)-prepad[d]])
    krn[indexesK...] = rot180(kern)
    AF = ifft(fft(A).*fft(krn))
    out = Array(T, size(img))
    indexesA = ntuple(N, d->postpad[d]+1:size(img,d)+postpad[d])
    copyreal!(out, AF, indexesA)
    out
end

for N = 1:5
    @eval begin
        function copyreal!{T<:Real}(dst::Array{T,$N}, src, I::(Range1{Int}...))
            @nexprs $N d->(I_d = I[d])
            @nloops $N i dst d->(j_d = first(I_d)+i_d-1) begin
                (@nref $N dst i) = real(@nref $N src j)
            end
            dst
        end
        function copyreal!{T<:Complex}(dst::Array{T,$N}, src, I::(Range1{Int}...))
            @nexprs $N d->I_d = I[d]
            @nloops $N i dst d->(j_d = first(I_d)+i_d-1) begin
                (@nref $N dst i) = @nref $N src j
            end
            dst
        end
    end
end

function imfilter_fft_old{T}(img::StridedMatrix{T}, filter::Matrix{T}, border::String, value)
    si, sf = size(img), size(filter)
    fw = iceil(([sf...] - 1) / 2)
    A = padarray(img, fw, fw, border, convert(T, value))
    # correlation instead of convolution
    filter = rot180(filter)
    # check if separable
    SVD = svdfact(filter)
    U, S, Vt = SVD[:U], SVD[:S], SVD[:Vt]
    separable = true
    for i = 2:length(S)
        separable &= (abs(S[i]) < sqrt(eps(T)))
    end
    if separable
        # conv2 isn't suitable for this (kernel center should be the actual center of the kernel)
        y = U[:,1]*sqrt(S[1])
        x = vec(Vt[1,:])*sqrt(S[1])
        sa = size(A)
        m = length(y)+sa[1]
        n = length(x)+sa[2]
        B = zeros(T, m, n)
        B[int(length(y)/2)+1:sa[1]+int(length(y)/2),int(length(x)/2)+1:sa[2]+int(length(x)/2)] = A
        yp = zeros(T, m)
        halfy = int((m-length(y)-1)/2)
        yp[halfy+1:halfy+length(y)] = y
        y = fft(yp)
        xp = zeros(T, n)
        halfx = int((n-length(x)-1)/2)
        xp[halfx+1:halfx+length(x)] = x
        x = fft(xp)
        C = fftshift(ifft(fft(B) .* (y * x.')))
        if T <: Real
            C = real(C)
        end
    else
        #C = conv2(A, filter)
        sa, sb = size(A), size(filter)
        At = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1)
        Bt = zeros(T, sa[1]+sb[1]-1, sa[2]+sb[2]-1)
        halfa1 = ifloor((size(At,1)-sa[1])/2)
        halfa2 = ifloor((size(At,2)-sa[2])/2)
        halfb1 = ifloor((size(Bt,1)-sb[1])/2)
        halfb2 = ifloor((size(Bt,2)-sb[2])/2)
        At[halfa1+1:halfa1+sa[1], halfa2+1:halfa2+sa[2]] = A
        Bt[halfb1+1:halfb1+sb[1], halfb2+1:halfb2+sb[2]] = filter
        C = fftshift(ifft(fft(At).*fft(Bt)))
        if T <: Real
            C = real(C)
        end
    end
    sc = size(C)
    out = C[int(sc[1]/2-si[1]/2):int(sc[1]/2+si[1]/2)-1, int(sc[2]/2-si[2]/2):int(sc[2]/2+si[2]/2)-1]
end

# IIR filtering with Gaussians
# See
#  Young, van Vliet, and van Ginkel, "Recursive Gabor Filtering",
#    IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 50: 2798-2805.
# and
#  Triggs and Sdika, "Boundary Conditions for Young - van Vliet
#    Recursive Filtering, IEEE TRANSACTIONS ON SIGNAL PROCESSING,
# Here we're using NA boundary conditions, so we set i- and i+
# (in Triggs & Sdika notation) to zero.
# Note these two papers use different sign conventions for the coefficients.

# Note: astype is ignored for FloatingPoint input
function imfilter_gaussian{T<:FloatingPoint}(img::AbstractArray{T}, sigma::Vector; emit_warning = true, astype::Type=Float64)
    if all(sigma .== 0)
        return img
    end
    A = copy(data(img))
    nanflag = isnan(A)
    hasnans = any(nanflag)
    if hasnans
        A[nanflag] = zero(T)
        validpixels = convert(Array{T}, !nanflag)
        imfilter_gaussian!(A, validpixels, sigma; emit_warning=emit_warning)
        A[nanflag] = nan(T)
    else
        imfilter_gaussian_no_nans!(A, sigma; emit_warning=emit_warning)
    end
    share(img, A)
end

function imfilter_gaussian{T<:Integer,TF<:FloatingPoint}(img::AbstractArray{T}, sigma::Vector; emit_warning = true, astype::Type{TF}=Float64)
    A = convert(Array{astype}, data(img))
    if all(sigma .== 0)
        return share(img, A)
    end
    imfilter_gaussian_no_nans!(A, sigma; emit_warning=emit_warning)
    share(img, A)
end

# This version is in-place, and destructive
# Any NaNs have to already be removed from data (and marked in validpixels)
function imfilter_gaussian!{T<:FloatingPoint}(data::Array{T}, validpixels::Array{T}, sigma::Vector; emit_warning = true)
    nd = ndims(data)
    if length(sigma) != nd
        error("Dimensionality mismatch")
    end
    _imfilter_gaussian!(data, sigma, emit_warning=emit_warning)
    _imfilter_gaussian!(validpixels, sigma, emit_warning=false)
    for i = 1:length(data)
        data[i] /= validpixels[i]
    end
    return data
end

# When there are no NaNs, the normalization is separable and hence can be computed far more efficiently
# This speeds the algorithm by approximately twofold
function imfilter_gaussian_no_nans!{T<:FloatingPoint}(data::Array{T}, sigma::Vector; emit_warning = true)
    nd = ndims(data)
    if length(sigma) != nd
        error("Dimensionality mismatch")
    end
    _imfilter_gaussian!(data, sigma, emit_warning=emit_warning)
    denom = Array(Vector{T}, nd)
    for i = 1:nd
        denom[i] = ones(T, size(data, i))
        if sigma[i] > 0
            _imfilter_gaussian!(denom[i], sigma[i:i], emit_warning=false)
        end
    end
    imfgnormalize!(data, denom)
    return data
end

for N = 1:5
    @eval begin
        function imfgnormalize!{T}(data::Array{T,$N}, denom)
            @nextract $N denom denom
            @nloops $N i data begin
                den = one(T)
                @nexprs $N d->(den *= denom_d[i_d])
                (@nref $N data i) /= den
            end
        end
    end
end

function iir_gaussian_coefficients(T::Type, sigma::Number; emit_warning::Bool = true)
    if sigma < 1 && emit_warning
        warn("sigma is too small for accuracy")
    end
    m0 = convert(T,1.16680)
    m1 = convert(T,1.10783)
    m2 = convert(T,1.40586)
    q = convert(T,1.31564*(sqrt(1+0.490811*sigma*sigma) - 1))
    scale = (m0+q)*(m1*m1 + m2*m2  + 2m1*q + q*q)
    B = m0*(m1*m1 + m2*m2)/scale
    B *= B
    # This is what Young et al call b, but in filt() notation would be called a
    a1 = q*(2*m0*m1 + m1*m1 + m2*m2 + (2*m0+4*m1)*q + 3*q*q)/scale
    a2 = -q*q*(m0 + 2m1 + 3q)/scale
    a3 = q*q*q/scale
    a = [-a1,-a2,-a3]
    Mdenom = (1+a1-a2+a3)*(1-a1-a2-a3)*(1+a2+(a1-a3)*a3)
    M = [-a3*a1+1-a3^2-a2      (a3+a1)*(a2+a3*a1)  a3*(a1+a3*a2);
          a1+a3*a2            -(a2-1)*(a2+a3*a1)  -(a3*a1+a3^2+a2-1)*a3;
          a3*a1+a2+a1^2-a2^2   a1*a2+a3*a2^2-a1*a3^2-a3^3-a3*a2+a3  a3*(a1+a3*a2)]/Mdenom;
    return a, B, M
end

function _imfilter_gaussian!{T<:FloatingPoint}(A::Array{T}, sigma::Vector; emit_warning::Bool = true)
    nd = ndims(A)
    szA = [size(A,i) for i = 1:nd]
    strdsA = [stride(A,i) for i = 1:nd]
    for d = 1:nd
        if sigma[d] == 0
            continue
        end
        if size(A, d) < 3
            error("All filtered dimensions must be of size 3 or larger")
        end
        a, B, M = iir_gaussian_coefficients(T, sigma[d], emit_warning=emit_warning)
        a1 = a[1]
        a2 = a[2]
        a3 = a[3]
        n1 = size(A,1)
        keepdims = [false,trues(nd-1)]
        if d == 1
            x = zeros(T, 3)
            vstart = zeros(T, 3)
            szhat = szA[keepdims]
            strdshat = strdsA[keepdims]
            if isempty(szhat)
                szhat = [1]
                strdshat = [1]
            end
            @inbounds @forcartesian c szhat begin
                coloffset = offset(c, strdshat)
                A[2+coloffset] -= a1*A[1+coloffset]
                A[3+coloffset] -= a1*A[2+coloffset] + a2*A[1+coloffset]
                for i = 4+coloffset:n1+coloffset
                    A[i] -= a1*A[i-1] + a2*A[i-2] + a3*A[i-3]
                end
                copytail!(x, A, coloffset, 1, n1)
                A_mul_B!(vstart, M, x)
                A[n1+coloffset] = vstart[1]
                A[n1-1+coloffset] -= a1*vstart[1]   + a2*vstart[2] + a3*vstart[3]
                A[n1-2+coloffset] -= a1*A[n1-1+coloffset] + a2*vstart[1] + a3*vstart[2]
                for i = n1-3+coloffset:-1:1+coloffset
                    A[i] -= a1*A[i+1] + a2*A[i+2] + a3*A[i+3]
                end
            end
        else
            x = Array(T, 3, n1)
            vstart = similar(x)
            keepdims[d] = false
            szhat = szA[keepdims]
            szd = szA[d]
            strdshat = strdsA[keepdims]
            strdd = strdsA[d]
            if isempty(szhat)
                szhat = [1]
                strdshat = [1]
            end
            @inbounds @forcartesian c szhat begin
                coloffset = offset(c, strdshat)  # offset for the remaining dimensions
                for i = 1:n1 A[i+strdd+coloffset] -= a1*A[i+coloffset] end
                for i = 1:n1 A[i+2strdd+coloffset] -= a1*A[i+strdd+coloffset] + a2*A[i+coloffset] end
                for j = 3:szd-1
                    jj = j*strdd+coloffset
                    for i = jj+1:jj+n1 A[i] -= a1*A[i-strdd] + a2*A[i-2strdd] + a3*A[i-3strdd] end
                end
                copytail!(x, A, coloffset, strdd, szd)
                A_mul_B!(vstart, M, x)
                for i = 1:n1 A[i+(szd-1)*strdd+coloffset] = vstart[1,i] end
                for i = 1:n1 A[i+(szd-2)*strdd+coloffset] -= a1*vstart[1,i]   + a2*vstart[2,i] + a3*vstart[3,i] end
                for i = 1:n1 A[i+(szd-3)*strdd+coloffset] -= a1*A[i+(szd-2)*strdd+coloffset] + a2*vstart[1,i] + a3*vstart[2,i] end
                for j = szd-4:-1:0
                    jj = j*strdd+coloffset
                    for i = jj+1:jj+n1 A[i] -= a1*A[i+strdd] + a2*A[i+2strdd] + a3*A[i+3strdd] end
                end
            end
        end
        @inbounds for i = 1:length(A)
            A[i] *= B
        end
    end
    A
end

function offset(c::Vector{Int}, strds::Vector{Int})
    o = 0
    for i = 1:length(c)
        o += (c[i]-1)*strds[i]
    end
    o
end

function copytail!(dest, A, coloffset, strd, len)
    for j = 1:3
        for i = 1:size(dest, 2)
            tmp = A[i + coloffset + (len-j)*strd]
            dest[j,i] = tmp
        end
    end
    dest
end

# Laplacian of Gaussian filter
# Separable implementation from Huertas and Medioni,
# IEEE Trans. Pat. Anal. Mach. Int., PAMI-8, 651, (1986)

function imfilter_LoG{T}(img::AbstractArray{T,2}, σ::Vector, border="replicate")
    # Limited to 2D for now.
    # See Sage D, Neumann F, Hediger F, Gasser S, Unser M.
    # Image Processing, IEEE Transactions on. (2005) 14(9):1372-1383.
    # for 3D.

    # Set up 1D kernels
    @assert length(σ) == 2
    σx, σy = σ[1], σ[2]
    h1(ξ, σ) = sqrt(1/(2π*σ^4))*(1 - ξ^2/σ^2)exp(-ξ^2/(2σ^2))
    h2(ξ, σ) = sqrt(1/(2π*σ^4))*exp(-ξ^2/(2σ^2))

    w = 8.5σx
    kh1x = Float64[h1(i, σx) for i = -floor(w/2):floor(w/2)]
    kh2x = Float64[h2(i, σx) for i = -floor(w/2):floor(w/2)]

    w = 8.5σy
    kh1y = Float64[h1(i, σy) for i = -floor(w/2):floor(w/2)]
    kh2y = Float64[h2(i, σy) for i = -floor(w/2):floor(w/2)]

    # Set up padding index lists
    kernlenx = length(kh1x)
    prepad  = div(kernlenx - 1, 2)
    postpad = div(kernlenx, 2)
    Ix = padindexes(img, 2, prepad, postpad, border)

    kernleny = length(kh1y)
    prepad  = div(kernleny - 1, 2)
    postpad = div(kernleny, 2)
    Iy = padindexes(img, 1, prepad, postpad, border)

    sz = size(img)
    # Store intermediate result in a transposed array
    # Allows column-major second stage filtering
    img1 = Array(Float64, (sz[2], sz[1]))
    img2 = Array(Float64, (sz[2], sz[1]))
    for j in 1:sz[2]
        for i in 1:sz[1]
            tmp1, tmp2 = 0.0, 0.0
            for k in 1:kernlenx
                @inbounds tmp1 += kh1x[k] * img[Ix[i + k - 1], j]
                @inbounds tmp2 += kh2x[k] * img[Ix[i + k - 1], j]
            end
            @inbounds img1[j, i] = tmp1 # Note the transpose
            @inbounds img2[j, i] = tmp2
        end
    end

    img12 = Array(Float64, sz) # Original image dims here
    img21 = Array(Float64, sz)
    for j in 1:sz[1]
        for i in 1:sz[2]
            tmp12, tmp21 = 0.0, 0.0
            for k in 1:kernleny
                @inbounds tmp12 += kh2y[k] * img1[Iy[i + k - 1], j]
                @inbounds tmp21 += kh1y[k] * img2[Iy[i + k - 1], j]
            end
            @inbounds img12[j, i] = tmp12 # Transpose back to original dims
            @inbounds img21[j, i] = tmp21
        end
    end
    img12 + img21
end

imfilter_LoG{T}(img::AbstractArray{T,2}, σ::Real, border="replicate") =
imfilter_LoG(img::AbstractArray{T,2}, [σ, σ], border)

function padindexes{T,n}(img::AbstractArray{T,n}, dim, prepad, postpad, border::String)
    M = size(img, dim)
    I = Array(Int, M + prepad + postpad)
    I = [(1 - prepad):(M + postpad)]
    if border == "replicate"
        I = min(max(I, 1), M)
    elseif border == "circular"
        I = 1 .+ mod(I .- 1, M)
    elseif border == "symmetric"
        I = [1:M, M:-1:1][1 .+ mod(I .- 1, 2 * M)]
    elseif border == "reflect"
        I = [1:M, M-1:-1:2][1 .+ mod(I .- 1, 2 * M - 2)]
    else
        error("unknown border condition")
    end
    I
end

### restrict, for reducing the image size by 2-fold
# This properly anti-aliases. The only "oddity" is that edges tend towards zero under
# repeated iteration.
# This implementation is considerably faster than the one in Grid,
# because it traverses the array in storage order.
if isdefined(:restrict)
    import Grid.restrict
end

restrict(img::AbstractImageDirect, ::()) = img

function restrict(img::AbstractImageDirect, region::Union(Dims, Vector{Int})=coords_spatial(img))
    A = data(img)
    for dim in region
        A = _restrict(A, dim)
    end
    props = copy(properties(img))
    ps = pixelspacing(img)
    ind = findin(coords_spatial(img), region)
    ps[ind] *= 2
    props["pixelspacing"] = ps
    Image(A, props)
end
function restrict(A::AbstractArray, region::Union(Dims, Vector{Int})=coords_spatial(img))
    for dim in region
        A = _restrict(A, dim)
    end
    A
end

function restrict{S<:ByteString}(img::AbstractImageDirect, region::Union((ByteString...), Vector{S}))
    so = spatialorder(img)
    regioni = Int[]
    for i = 1:length(region)
        push!(regioni, require_dimindex(img, region[i], so))
    end
    restrict(img, regioni)
end

function _restrict(A, dim)
    if size(A, dim) <= 2
        return A
    end
    out = Array(typeof(A[1]/4+A[2]/2), ntuple(ndims(A), i->i==dim?restrict_size(size(A,dim)):size(A,i)))
    restrict!(out, A, dim)
    out
end

# out should have efficient linear indexing
for N = 1:5
    @eval begin
        function restrict!{T}(out::AbstractArray{T,$N}, A::AbstractArray, dim)
            if isodd(size(A, dim))
                half = convert(T, 0.5)
                quarter = convert(T, 0.25)
                indx = 0
                if dim == 1
                    @nloops $N i d->(d==1 ? (1:1) : (1:size(A,d))) d->(j_d = d==1 ? i_d+1 : i_d) begin
                        nxt = @nref $N A j
                        out[indx+=1] = half*(@nref $N A i) + quarter*nxt
                        for k = 3:2:size(A,1)-2
                            prv = nxt
                            i_1 = k
                            j_1 = k+1
                            nxt = @nref $N A j
                            out[indx+=1] = quarter*(prv+nxt) + half*(@nref $N A i)
                        end
                        i_1 = size(A,1)
                        out[indx+=1] = quarter*nxt + half*(@nref $N A i)
                    end
                else
                    strd = stride(out, dim)
                    # Must initialize the i_dim==1 entries with zero
                    @nexprs $N d->sz_d=d==dim?1:size(out,d)
                    @nloops $N i d->(1:sz_d) begin
                        (@nref $N out i) = 0
                    end
                    stride_1 = 1
                    @nexprs $N d->(stride_{d+1} = stride_d*size(out,d))
                    @nexprs $N d->offset_d = 0
                    ispeak = true
                    @nloops $N i d->(d==1?(1:1):(1:size(A,d))) d->(if d==dim; ispeak=isodd(i_d); offset_{d-1} = offset_d+(div(i_d+1,2)-1)*stride_d; else; offset_{d-1} = offset_d+(i_d-1)*stride_d; end) begin
                        indx = offset_0
                        if ispeak
                            for k = 1:size(A,1)
                                i_1 = k
                                out[indx+=1] += half*(@nref $N A i)
                            end
                        else
                            for k = 1:size(A,1)
                                i_1 = k
                                tmp = quarter*(@nref $N A i)
                                out[indx+=1] += tmp
                                out[indx+strd] = tmp
                            end
                        end
                    end
                end
            else
                threeeighths = convert(T, 0.375)
                oneeighth = convert(T, 0.125)
                indx = 0
                if dim == 1
                    @nloops $N i d->(d==1 ? (1:1) : (1:size(A,d))) d->(j_d = i_d) begin
                        c = d = zero(T)
                        for k = 1:size(out,1)-1
                            a = c
                            b = d
                            j_1 = 2*k
                            i_1 = j_1-1
                            c = @nref $N A i
                            d = @nref $N A j
                            out[indx+=1] = oneeighth*(a+d) + threeeighths*(b+c)
                        end
                        out[indx+=1] = oneeighth*c+threeeighths*d
                    end
                else
                    fill!(out, 0)
                    strd = stride(out, dim)
                    stride_1 = 1
                    @nexprs $N d->(stride_{d+1} = stride_d*size(out,d))
                    @nexprs $N d->offset_d = 0
                    peakfirst = true
                    @nloops $N i d->(d==1?(1:1):(1:size(A,d))) d->(if d==dim; peakfirst=isodd(i_d); offset_{d-1} = offset_d+(div(i_d+1,2)-1)*stride_d; else; offset_{d-1} = offset_d+(i_d-1)*stride_d; end) begin
                        indx = offset_0
                        if peakfirst
                            for k = 1:size(A,1)
                                i_1 = k
                                tmp = @nref $N A i
                                out[indx+=1] += threeeighths*tmp
                                out[indx+strd] += oneeighth*tmp
                            end
                        else
                            for k = 1:size(A,1)
                                i_1 = k
                                tmp = @nref $N A i
                                out[indx+=1] += oneeighth*tmp
                                out[indx+strd] += threeeighths*tmp
                            end
                        end
                    end
                end
            end
        end
    end
end

restrict_size(len::Integer) = isodd(len) ? (len+1)>>1 : (len>>1)+1

function imlineardiffusion{T}(img::Array{T,2}, dt::FloatingPoint, iterations::Integer)
    u = img
    f = imlaplacian()
    for i = dt:dt:dt*iterations
        u = u + dt*imfilter(u, f, "replicate")
    end
    u
end

function imgaussiannoise{T}(img::AbstractArray{T}, variance::Number, mean::Number)
    return img + sqrt(variance)*randn(size(img)) + mean
end

imgaussiannoise{T}(img::AbstractArray{T}, variance::Number) = imgaussiannoise(img, variance, 0)
imgaussiannoise{T}(img::AbstractArray{T}) = imgaussiannoise(img, 0.01, 0)

function rgb2ntsc{T}(img::Array{T})
    trans = [0.299 0.587 0.114; 0.596 -0.274 -0.322; 0.211 -0.523 0.312]
    out = zeros(T, size(img))
    for i = 1:size(img,1), j = 1:size(img,2)
        out[i,j,:] = trans * vec(img[i,j,:])
    end
    return out
end

function ntsc2rgb{T}(img::Array{T})
    trans = [1 0.956 0.621; 1 -0.272 -0.647; 1 -1.106 1.703]
    out = zeros(T, size(img))
    for i = 1:size(img,1), j = 1:size(img,2)
        out[i,j,:] = trans * vec(img[i,j,:])
    end
    return out
end

function rgb2ycbcr{T}(img::Array{T})
    trans = [65.481 128.533 24.966; -37.797 -74.203 112; 112 -93.786 -18.214]
    offset = [16.0; 128.0; 128.0]
    out = zeros(T, size(img))
    for i = 1:size(img,1), j = 1:size(img,2)
        out[i,j,:] = offset + trans * vec(img[i,j,:])
    end
    return out
end

function ycbcr2rgb{T}(img::Array{T})
    trans = inv([65.481 128.533 24.966; -37.797 -74.203 112; 112 -93.786 -18.214])
    offset = [16.0; 128.0; 128.0]
    out = zeros(T, size(img))
    for i = 1:size(img,1), j = 1:size(img,2)
        out[i,j,:] = trans * (vec(img[i,j,:]) - offset)
    end
    return out
end

function imcomplement{T}(img::AbstractArray{T})
    l = limits(img)
    if l[2] != 1
        error("imcomplement not defined unless upper limit is 1")
    end
    return 1 - img
end

function rgb2hsi{T}(img::Array{T})
    R = img[:,:,1]
    G = img[:,:,2]
    B = img[:,:,3]
    H = acos((1/2*(2*R - G - B)) ./ (((R - G).^2 + (R - B).*(G - B)).^(1/2)+eps(T)))
    H[B .> G] = 2*pi - H[B .> G]
    H /= 2*pi
    rgb_sum = R + G + B
    rgb_sum[rgb_sum .== 0] = eps(T)
    S = 1 - 3./(rgb_sum).*min(R, G, B)
    H[S .== 0] = 0
    I = 1/3*(R + G + B)
    return cat(3, H, S, I)
end

function hsi2rgb{T}(img::Array{T})
    H = img[:,:,1]*(2pi)
    S = img[:,:,2]
    I = img[:,:,3]
    R = zeros(T, size(img,1), size(img,2))
    G = zeros(T, size(img,1), size(img,2))
    B = zeros(T, size(img,1), size(img,2))
    RG = 0 .<= H .< 2*pi/3
    GB = 2*pi/3 .<= H .< 4*pi/3
    BR = 4*pi/3 .<= H .< 2*pi
    # RG sector
    B[RG] = I[RG].*(1 - S[RG])
    R[RG] = I[RG].*(1 + (S[RG].*cos(H[RG]))./cos(pi/3 - H[RG]))
    G[RG] = 3*I[RG] - R[RG] - B[RG]
    # GB sector
    R[GB] = I[GB].*(1 - S[GB])
    G[GB] = I[GB].*(1 + (S[GB].*cos(H[GB] - pi/3))./cos(H[GB]))
    B[GB] = 3*I[GB] - R[GB] - G[GB]
    # BR sector
    G[BR] = I[BR].*(1 - S[BR])
    B[BR] = I[BR].*(1 + (S[BR].*cos(H[BR] - 2*pi/3))./cos(-pi/3 - H[BR]))
    R[BR] = 3*I[BR] - G[BR] - B[BR]
    return cat(3, R, G, B)
end

function imstretch{T}(img::AbstractArray{T}, m::Number, slope::Number)
    assert_scalar_color(img)
    if limits(img) != (0,1)
        warn("Image limits ", limits(img), " are not (0,1)")
    end
    share(img, 1./(1 + (m./(data(img) + eps(T))).^slope))
end

# image gradients

# forward and backward differences
# can be very helpful for discretized continuous models
forwarddiffy{T}(u::Array{T,2}) = [u[2:end,:]; u[end,:]] - u
forwarddiffx{T}(u::Array{T,2}) = [u[:,2:end] u[:,end]] - u
backdiffy{T}(u::Array{T,2}) = u - [u[1,:]; u[1:end-1,:]]
backdiffx{T}(u::Array{T,2}) = u - [u[:,1] u[:,1:end-1]]

function imROF{T}(img::Array{T,2}, lambda::Number, iterations::Integer)
    # Total Variation regularized image denoising using the primal dual algorithm
    # Also called Rudin Osher Fatemi (ROF) model
    # lambda: regularization parameter
    s1, s2 = size(img)
    p = zeros(T, s1, s2, 2)
    u = zeros(T, s1, s2)
    grad_u = zeros(T, s1, s2, 2)
    div_p = zeros(T, s1, s2)
    dt = lambda/4
    for i = 1:iterations
        div_p = backdiffx(p[:,:,1]) + backdiffy(p[:,:,2])
        u = img + div_p/lambda
        grad_u = cat(3, forwarddiffx(u), forwarddiffy(u))
        grad_u_mag = sqrt(grad_u[:,:,1].^2 + grad_u[:,:,2].^2)
        tmp = 1 + grad_u_mag*dt
        p = (dt*grad_u + p)./cat(3, tmp, tmp)
    end
    return u
end

# ROF Model for color images
function imROF(img::AbstractArray, lambda::Number, iterations::Integer)
    cd = colordim(img)
    local out
    if cd != 0
        out = similar(img)
        for i = size(img, cd)
            imsl = img["color", i]
            outsl = slice(out, "color", i)
            copy!(outsl, imROF(imsl, lambda, iterations))
        end
    else
        out = share(img, imROF(data(img), lambda, iterations))
    end
    out
end

### Morphological operations

# Erode and dilate support 3x3 regions only (and higher-dimensional generalizations).
dilate(img::AbstractArray, region=coords_spatial(img)) = dilate!(copy(img), region)
erode(img::AbstractArray, region=coords_spatial(img)) = erode!(copy(img), region)

dilate!(maxfilt, region=coords_spatial(maxfilt)) = extremefilt!(data(maxfilt), Base.Order.Forward, region)
erode!(minfilt, region=coords_spatial(minfilt)) = extremefilt!(data(minfilt), Base.Order.Reverse, region)
function extremefilt!(extrfilt::Array, order::Ordering, region=coords_spatial(extrfilt))
    for d = 1:ndims(extrfilt)
        if size(extrfilt, d) == 1 || !in(d, region)
            continue
        end
        sz = [size(extrfilt,i) for i = 1:ndims(extrfilt)]
        s = stride(extrfilt, d)
        sz[d] = 1
        @forcartesian i sz begin
            k = cartesian_linear(extrfilt, i)
            a2 = extrfilt[k]
            a3 = extrfilt[k+s]
            extrfilt[k] = extr(order, a2, a3)
            for l = 2:size(extrfilt,d)-1
                k += s
                a1 = a2
                a2 = a3
                a3 = extrfilt[k+s]
                extrfilt[k] = extr(order, a1, a2, a3)
            end
            extrfilt[k+s] = extr(order, a2, a3)
        end
    end
    extrfilt
end

opening(img::AbstractArray, region=coords_spatial(img)) = opening!(copy(img), region)
opening!(img::AbstractArray, region=coords_spatial(img)) = dilate!(erode!(img, region),region)
closing(img::AbstractArray, region=coords_spatial(img)) = closing!(copy(img), region)
closing!(img::AbstractArray, region=coords_spatial(img)) = erode!(dilate!(img, region),region)

extr(order::ForwardOrdering, x::Real, y::Real) = max(x,y)
extr(order::ForwardOrdering, x::Real, y::Real, z::Real) = max(x,y,z)
extr(order::ReverseOrdering, x::Real, y::Real) = min(x,y)
extr(order::ReverseOrdering, x::Real, y::Real, z::Real) = min(x,y,z)

extr(order::Ordering, x::RGB, y::RGB) = RGB(extr(order, x.r, y.r), extr(order, x.g, y.g), extr(order, x.b, y.b))
extr(order::Ordering, x::RGB, y::RGB, z::RGB) = RGB(extr(order, x.r, y.r, z.r), extr(order, x.g, y.g, z.g), extr(order, x.b, y.b, z.b))

extr(order::Ordering, x::ColorValue, y::ColorValue) = extr(order, convert(RGB, x), convert(RGB, y))
extr(order::Ordering, x::ColorValue, y::ColorValue, z::ColorValue) = extr(order, convert(RGB, x), convert(RGB, y), convert(RGB, z))

# phantom images

function shepp_logan(M,N; highContrast=true)
  # Initially proposed in Shepp, Larry; B. F. Logan (1974).
  # "The Fourier Reconstruction of a Head Section". IEEE Transactions on Nuclear Science. NS-21.

  P = zeros(M,N)

  x = linspace(-1,1,M)'
  y = linspace(1,-1,N)

  centerX = [0, 0, 0.22, -0.22, 0, 0, 0, -0.08, 0, 0.06]
  centerY = [0, -0.0184, 0, 0, 0.35, 0.1, -0.1, -0.605, -0.605, -0.605]
  majorAxis = [0.69, 0.6624, 0.11, 0.16, 0.21, 0.046, 0.046, 0.046, 0.023, 0.023]
  minorAxis = [0.92, 0.874, 0.31, 0.41, 0.25, 0.046, 0.046, 0.023, 0.023, 0.046]
  theta = [0, 0, -18.0, 18.0, 0, 0, 0, 0, 0, 0]

  # original (CT) version of the phantom
  grayLevel = [2, -0.98, -0.02, -0.02, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]

  if(highContrast)
    # high contrast (MRI) version of the phantom
    grayLevel = [1, -0.8, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  end

  for l=1:length(theta)
    P += grayLevel[l] * (
           ( (cos(theta[l] / 360*2*pi) * (x .- centerX[l]) .+
              sin(theta[l] / 360*2*pi) * (y .- centerY[l])) / majorAxis[l] ).^2 .+
           ( (sin(theta[l] / 360*2*pi) * (x .- centerX[l]) .-
              cos(theta[l] / 360*2*pi) * (y .- centerY[l])) / minorAxis[l] ).^2 .< 1 )
  end

  return P
end

shepp_logan(N;highContrast=true) = shepp_logan(N,N;highContrast=highContrast)


text:
module Grids

export
    Cell, Face,
    AbstractTopology, Periodic, Bounded, Flat, topology,
    AbstractGrid, RegularCartesianGrid, VerticallyStretchedCartesianGrid,
    xnode, ynode, znode, xnodes, ynodes, znodes, nodes,
    xC, xF, yC, yF, zC, zF

import Base: size, length, eltype, show

import Oceananigans: short_show

using Oceananigans
using Oceananigans.Architectures

using OffsetArrays

#####
##### Abstract types
#####

"""
    Cell

A type describing the location at the center of a grid cell.
"""
struct Cell end

"""
	Face

A type describing the location at the face of a grid cell.
"""
struct Face end

"""
    AbstractTopology

Abstract supertype for grid topologies.
"""
abstract type AbstractTopology end

"""
    Periodic

Grid topology for periodic dimensions.
"""
struct Periodic <: AbstractTopology end

"""
    Bounded

Grid topology for bounded dimensions. These could be wall-bounded dimensions
or dimensions
"""
struct Bounded <: AbstractTopology end

"""
    Flat

Grid topology for flat dimensions, generally with one grid point, along which the solution
is uniform and does not vary.
"""
struct Flat <: AbstractTopology end

"""
    AbstractGrid{FT, TX, TY, TZ}

Abstract supertype for grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractGrid{FT, TX, TY, TZ} end

Base.eltype(::AbstractGrid{FT}) where FT = FT
Base.size(grid::AbstractGrid) = (grid.Nx, grid.Ny, grid.Nz)
Base.length(grid::AbstractGrid) = (grid.Lx, grid.Ly, grid.Lz)

halo_size(grid) = (grid.Hx, grid.Hy, grid.Hz)

topology(::AbstractGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = (TX, TY, TZ)
topology(grid, dim) = topology(grid)[dim]

include("grid_utils.jl")
include("input_validation.jl")
include("regular_cartesian_grid.jl")
include("vertically_stretched_cartesian_grid.jl")

end

ext: jl


text:
function sqlite3_table_column_metadata(
    dbptr::Ptr{Void},
    table::String,
    column::String,
    datatype::Vector{Ptr{Uint8}},
    collseq::Vector{Ptr{Uint8}},
    notnull::Vector{Cint},
    primarykey::Vector{Cint},
    autoinc::Vector{Cint}
)
    return ccall(
        (:sqlite3_table_column_metadata, sqlite3_lib),
        Cint,
        (
            Ptr{Void},
            Ptr{Uint8},
            Ptr{Uint8},
            Ptr{Uint8},
            Ptr{Ptr{Uint8}},
            Ptr{Ptr{Uint8}},
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Cint}
        ),
        dbptr,
        convert(Ptr{Uint8}, C_NULL),
        table,
        column,
        datatype,
        collseq,
        notnull,
        primarykey,
        autoinc
    )
end


text:
module AtomicFileWrite

export atomic_write

function atomic_write(f, path; backup=nothing, overwrite_backup=false)
    mktemp() do temp_path, temp_io
        f(temp_io)
        backup == nothing || cp(path, backup; force=overwrite_backup)
        close(temp_io)
        _replace(temp_path, path)
    end
end

# Unlike Base.mv, no fallback to cp && rm on error, to ensure atomicity
function _replace(src, dst)
    err = ccall(:jl_fs_rename, Int32, (Cstring, Cstring), src, dst)
    Base.uv_error("Failed to replace $dst with $src", err)
end

end


text:
import Knet.Ops20: elu, relu, selu, sigm, eluback, reluback, seluback, sigmback
import Base.Broadcast: broadcasted
import Knet
using Knet.KnetArrays: KnetArray, DevArray, Bcasted
using Knet.LibKnet8: @knet8
using CUDA: CuArray, CuPtr
using AutoGrad: AutoGrad, @primitive

# Specialize tanh gradient for DevArrays here. The others have been declared in Ops20 as generic gradients.
tanhback(dyi::T,yi::T) where {T<:Number} = dyi*(T(1)-yi*yi)
@primitive tanh(x::DevArray),dy,y tanhback.(dy,y)
@primitive tanhback(dy,y),ddx  ddx.*(1 .- y.*y)  ddx.*(-2 .* dy.*y)

for (R,P) in ((KnetArray,Ptr), (CuArray,CuPtr)), T in (Float32,Float64); S = sizeof(T) * 8
    for f in ("elu","relu","selu","sigm")
        J, F = Symbol(f), "$(f)_$S"; M = which(@__MODULE__,J)
        @eval begin
            function broadcasted(::typeof($J),x::$R{$T})
                y = similar(x)
                @knet8($F,(Cint,$P{$T},$P{$T}),length(y),x,y)
                return y
            end
            # Bcasted methods -- only needed for KnetArray
            ($M).$J(x::Bcasted{<:$R{$T}}) = broadcasted($J, x.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:$R{$T}}) = broadcasted($J, x.value) |> Bcasted
        end
    end
    for f in ("eluback","reluback","seluback","sigmback","tanhback")
        J, F = Symbol(f), "$(f)_$(S)_11"; M = which(@__MODULE__,J)
        @eval begin
            function broadcasted(::typeof($J),x::$R{$T},y::$R{$T})
                z = similar(x)
                @knet8($F,(Cint,$P{$T},$P{$T},$P{$T}),length(z),x,y,z)
                return z
            end
            # Bcasted methods -- only needed for KnetArray
            ($M).$J(x::Bcasted{<:$R{$T}}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            ($M).$J(x::$R{$T}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            ($M).$J(x::Bcasted{<:$R{$T}}, y::$R{$T}) = broadcasted($J, x.value, y) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:$R{$T}}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x.value, y.value) |> Bcasted
            broadcasted(::typeof($J),x::$R{$T}, y::Bcasted{<:$R{$T}}) = broadcasted($J, x, y.value) |> Bcasted
            broadcasted(::typeof($J),x::Bcasted{<:$R{$T}}, y::$R{$T}) = broadcasted($J, x.value, y) |> Bcasted
        end
    end
end


text:
@testset "λ T" begin
    for n_rule = 1:7
        # table = Radau.RadauTable(k)
        # @test table.T * diagm(0=>table.λ) * table.T⁻¹ ≈ inv(table.A)

        A, b, c, bi, b_hat = radau_butcher_table_plus(2 * n_rule - 1)
        inv_A = inv(A)
        λ, T = big_eigen(inv_A)

        @test T * diagm(0=>λ) * inv(T) ≈ inv_A

        if isodd(n_rule)
            @test λ[1] == real.(λ[1])
        end
    end
end


text:
using HypothesisTests, Base.Test

# Basic exact test
@test abs(pvalue(ExactSignedRankTest([1:10;], [2:2:20;])) - 0.0020) <= 1e-4
@test abs(pvalue(ExactSignedRankTest([2:2:20;], [1:10;])) - 0.0020) <= 1e-4
@test abs(pvalue(ExactSignedRankTest([1:10;], [2:2:16; -1; 1])) - 0.4316) <= 1e-4
@test abs(pvalue(ExactSignedRankTest([2:2:16; -1; 1], [1:10;])) - 0.4316) <= 1e-4
show(IOBuffer(), ExactSignedRankTest([1:10;], [2:2:20;]))

# Exact with ties
@test abs(pvalue(ExactSignedRankTest([1:10;], [1:10;])) - 1) <= 1e-4
@test abs(pvalue(ExactSignedRankTest([1:10;], [2:11;])) - 0.0020) <= 1e-4
@test abs(pvalue(ExactSignedRankTest([2:11;], [1:10;])) - 0.0020) <= 1e-4
@test abs(pvalue(ExactSignedRankTest(1:10, [1:5; ones(5)])) - 0.0625) <= 1e-4
@test abs(pvalue(ExactSignedRankTest([1:5; ones(5)], [1:10;])) - 0.0625) <= 1e-4
show(IOBuffer(), ExactSignedRankTest([1:10;], [1:10;]))

# Approximate test
@test abs(pvalue(ApproximateSignedRankTest([1:10;], [2:2:20;])) - 0.005922) <= 1e-6
@test abs(pvalue(ApproximateSignedRankTest([2:2:20;], [1:10;])) - 0.005922) <= 1e-6
@test abs(pvalue(ApproximateSignedRankTest([1:10;], [2:2:16; -1; 1])) - 0.4148) <= 1e-4
@test abs(pvalue(ApproximateSignedRankTest([2:2:16; -1; 1], [1:10;])) - 0.4148) <= 1e-4
show(IOBuffer(), ApproximateSignedRankTest([1:10;], [2:2:20;]))

# Approximate with ties
@test abs(pvalue(ApproximateSignedRankTest([1:10;], [1:10;])) - 1) <= 1e-4
@test abs(pvalue(ApproximateSignedRankTest([1:10;], [2:11;])) - 0.001904) <= 1e-6
@test abs(pvalue(ApproximateSignedRankTest([2:11;], [1:10;])) - 0.001904) <= 1e-6
@test abs(pvalue(ApproximateSignedRankTest([1:10;], [1:5; ones(5)])) - 0.05906) <= 1e-5
@test abs(pvalue(ApproximateSignedRankTest([1:5; ones(5)], 1:10)) - 0.05906) <= 1e-5
show(IOBuffer(), ApproximateSignedRankTest([1:10;], [1:10;]))

# # Tests for automatic selection
@test abs(pvalue(SignedRankTest([1:10;], [2:2:20;])) - 0.0020) <= 1e-4
@test abs(pvalue(SignedRankTest([1:10;], [2:11;])) - 0.0020) <= 1e-4
show(IOBuffer(), SignedRankTest([1:10;], [2:2:20;]))

# One Sample tests
# P-value computed using R wilcox.test
@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] - 10.1)) - 0.09460449) <= 1e-4
# P-value computed using R wilcox.test
@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16] - 10.1)) - 0.1928101) <= 1e-4

# One Sample tests with ties
# P-value computed using R package exactRankTests wilcox.exact
@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,10,10,10,10,10,13,14,15] - 10.1)) - 0.04052734) <= 1e-4
# P-value computed using R wilcox.test
@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,10,10,10,10,10,13,14,15,16] - 10.1)) - 0.1021964) <= 1e-4

# Test confidence interval
x = [-7.8, -6.9, -4.7, 3.7, 6.5, 8.7, 9.1, 10.1, 10.8, 13.6, 14.4, 16.6, 20.2, 22.4, 23.5]
@test_approx_eq_eps ci(ExactSignedRankTest(x))[1] 3.3 1e-4
@test_approx_eq_eps ci(ExactSignedRankTest(x))[2] 15.5 1e-4
@test_approx_eq_eps ci(ApproximateSignedRankTest(x))[1] 3.3 1e-4
@test_approx_eq_eps ci(ApproximateSignedRankTest(x))[2] 15.5 1e-4
@test_approx_eq_eps ci(SignedRankTest(x); tail=:left)[1] 4.45 1e-4
@test_approx_eq_eps ci(SignedRankTest(x); tail=:right)[2] 14.45 1e-4


text:
@testset "data" begin
	fs = joinpath(PdfFraud.proto_path, "protos/ff953f62f5739a910f96465c34367e31f6f4fdc112b8a27f66a7d88dd519ec91.pb")

	pdf = PdfFraud.read_pb(fs);
	@test typeof(pdf) == PdfSampleRaw
 	vs = map(i -> PdfFraud.visual_representation(pdf,i), 1:length(pdf.pages))
	@test length(vs) == 8
	xs = map(x -> reshape(x, 64, 48), vs)
end


text:
using PyPlot, .ScientificColourMaps

cmaps = [acton,bamako,batlow,berlin,bilbao]

gradient = collect(0:256)
G = cat(dims=2,gradient,gradient)'

fix,axs = subplots(5,1)

for i in 1:5
    axs[i].imshow(G,aspect="auto",cmap=cmaps[i])
    axs[i].set_xticks([])
    axs[i].set_yticks([])
    axs[i].set_ylabel(cmaps[i].name)
end


text:
#!/usr/bin/env julia

# Copyright 2022 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
module project3

using Gridap
using GridapDistributed
using PartitionedArrays
using MPI

function fe_solver(parts,
                   nx::Integer=10,
                   ny::Integer=10,
                   left_bc::Real=15.0,
                   right_bc::Real=5.0,
                   write_plot_files::Bool=false)

    #############################
    ####### ADD CODE HERE #######
    #############################

    # ph = assign the outpt of solve

    # Uncomment the code below once everything is implemented

    # if write_plot_files
        # writevtk(Ω, "results", cellfields=["ph"=>ph])
    # end
    # ph
end

if abspath(PROGRAM_FILE) == @__FILE__
    MPI.Init()
    size_ = MPI.Comm_size(MPI.COMM_WORLD)
    partition = (size_ ÷ 2, size_ ÷ 2)
    ph = prun(x -> fe_solver(x, 100, 100, 10.0, 10.0, true), mpi, partition)
end

export fe_solver

end


text:
export lorentz_factors, LorentzModel, setup!, apply!

# Deprecated, will be moved to Osnaps.jl in the future.
lorentz_factors(::Type{T}; ϵc, ωp, fs, Ωs, γs) where T = lorentz_factors(T, ϵc, ωp, fs, Ωs, γs)

function lorentz_factors(::Type{Vector}, ϵc::Real, ωp::Real, fs::NTuple{N}, Ωs::NTuple{N}, γs::NTuple{N}) where N
    if @generated
        a = Vector{Any}(undef, 3*N+2)
        @inbounds a[1] = :ϵc
        @inbounds a[2] = :ωp
        @inbounds for i in 1:N
            a[i+2]     = :(fs[$i])
            a[i+2+  N] = :(Ωs[$i])
            a[i+2+2*N] = :(γs[$i])
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:vect, a...))
        end
    else
        a = Vector{Float64}(undef, 3*N+2)
        @inbounds a[1] = ϵc
        @inbounds a[2] = ωp
        @inbounds for i in 1:N
            a[i+2], a[i+2+N], a[i+2+2*N] = fs[i], Ωs[i], γs[i]
        end
        return a
    end
end

function lorentz_factors(::Type{Tuple}, ϵc::Real, ωp::Real, fs::NTuple{N}, Ωs::NTuple{N}, γs::NTuple{N}) where N
    if @generated
        a = Vector{Any}(undef, 3*N+2)
        @inbounds a[1] = :ϵc
        @inbounds a[2] = :ωp
        @inbounds for i in 1:N
            a[i+2]     = :(fs[$i])
            a[i+2+  N] = :(Ωs[$i])
            a[i+2+2*N] = :(γs[$i])
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:tuple, a...))
        end
    else
        return ntuple(i -> lorentz_factors(i, ϵc, ωp, fs, Ωs, γs), 3*N+2)
    end
end

function lorentz_factors(ix::Int, ϵc::Real, ωp::Real, fs::NTuple{N}, Ωs::NTuple{N}, γs::NTuple{N}) where N
    ix <= 0  && error("BoundsError: lorentz_factors(ix = $ix, ...)")
    ix == 1  && return ϵc
    ix == 2  && return ωp
    lx  = 2; rx = 2 + N
    ix <= rx && return fs[ix - lx]
    lx += N; rx += N
    ix <= rx && return Ωs[ix - lx]
    lx += N; rx += N
    ix <= rx && return γs[ix - lx]
    error("BoundsError: lorentz_factors(ix = $ix, ...)")
end

struct LorentzModel
    re_kernel::SubArray{Float64, 2, Array{Float64,3}, Tuple{Base.Slice{Base.OneTo{Int}}, Base.Slice{Base.OneTo{Int}}, Int}, true}
    im_kernel::SubArray{Float64, 2, Array{Float64,3}, Tuple{Base.Slice{Base.OneTo{Int}}, Base.Slice{Base.OneTo{Int}}, Int}, true}
    re_params::SubArray{Float64, 1, Array{Float64,2}, Tuple{Base.Slice{Base.OneTo{Int}}, Int}, true}
    im_params::SubArray{Float64, 1, Array{Float64,2}, Tuple{Base.Slice{Base.OneTo{Int}}, Int}, true}
    kernel::Array{Float64,3}
    params::Array{Float64,2}
    ratios::Array{Float64,1}
    nExp::Int
    nLor::Int

    function LorentzModel(nExp::Int, nLor::Int)
        kernel = Array{Float64,3}(undef, nExp, nLor, 2)
        params = Array{Float64,2}(undef, nExp, 2)
        ratios = Array{Float64,1}(undef, nLor)
        return new(view(kernel, :, :, 1), view(kernel, :, :, 2), view(params, :, 1), view(params, :, 2), kernel, params, ratios, nExp, nLor)
    end
end

function setup!(o::LorentzModel, fs::VecI)
    re_params = o.re_params
    im_params = o.im_params
    @simd for i in eachindex(1:o.nExp)
        @inbounds re_params[i], im_params[i] = abs2(fs[i]), fs[i]
    end
    return nothing
end

function setup!(o::LorentzModel, fs::VecI, Ωs::VecI, γs::VecI)
    one2nExp  = eachindex(1:o.nExp)
    ratios    = o.ratios
    re_params = o.re_params
    im_params = o.im_params
    re_kernel = o.re_kernel
    im_kernel = o.im_kernel

    for j in 1:o.nLor
        @inbounds Ωj2 = abs2(Ωs[j])
        @inbounds γj  = γs[j]
        @simd for i in one2nExp
            @inbounds re_kernel[i,j], im_kernel[i,j] = Ωj2 - re_params[i], γj * im_params[i]
        end
        @inbounds for i in one2nExp
            denom = sqr2(re_kernel[i,j], im_kernel[i,j])
            re_kernel[i,j] /= denom
            im_kernel[i,j] /= denom
        end
        @inbounds ratios[j] = fs[j]
    end
    return nothing
end

function apply!(ϵ1::VecIO, ϵ2::VecIO, o::LorentzModel, ϵc::Real, ωp::Real)
    @simd for i in eachindex(1:o.nExp)
        @inbounds ϵ1[i], ϵ2[i] = ϵc, 0.0
    end

    ωp2 = ωp * ωp
    gemv!('N', ωp2, o.re_kernel, o.ratios, true, ϵ1)
    gemv!('N', ωp2, o.im_kernel, o.ratios, true, ϵ2)
    return nothing
end


text:
"""
    val, d = parse_from_keys(d::Dict, symbs)

This function is based on find_in_dict() in GMT.
Check if d (Dict) contains either of symbs (list of symbol).
If true, return corresponding value
"""
function parse_from_keys(d::Dict, symbs)
	for symb in symbs
        if haskey(d, symb)
            val = d[symb]
		    delete!(d,symb)
		    return val, d
		end
	end
	# return if nothing
	return nothing, d
end

# Plots
# -----------------------------
"""
key: seriestype
To parse a value of the key in kwargs
"""
parse_seriestype(d::Dict) =
parse_from_keys(d, [:seriestype, :st, :t, :typ, :linetype, :lt])
# -----------------------------
parse_seriescolor(d::Dict) =
parse_from_keys(d, [:seriescolor, :c, :color, :colour])
# -----------------------------
parse_clims(d::Dict) =
parse_from_keys(d, [:clims, :clim, :cbarlims, :cbar_lims, :climits, :color_limits])
# -----------------------------
parse_bgcolor_inside(d::Dict) =
parse_from_keys(d, [:background_color_inside, :bg_inside, :bginside,
                    :bgcolor_inside, :bg_color_inside, :background_inside,
                    :background_colour_inside, :bgcolour_inside, :bg_colour_inside])
# -----------------------------
parse_xlims(d::Dict) =
parse_from_keys(d, [:xlims, :xlim, :xlimit, :xlimits])
# -----------------------------
parse_ylims(d::Dict) =
parse_from_keys(d, [:ylims, :ylim, :ylimit, :ylimits])
# -----------------------------
parse_colorbar(d::Dict) =
parse_from_keys(d, [:colorbar,:cb, :cbar, :colorkey])
# -----------------------------
parse_colorbar_title(d::Dict) =
parse_from_keys(d, [:colorbar_title])
# -----------------------------
parse_linewidth(d::Dict) =
parse_from_keys(d, [:linewidth, :w, :width, :lw])
# -----------------------------
parse_linestyle(d::Dict) =
parse_from_keys(d, [:linestyle, :style, :s, :ls])
# -----------------------------
parse_linecolor(d::Dict) =
parse_from_keys(d, [:linecolor, :lc, :lcolor, :lcolour, :linecolour])
# -----------------------------
parse_label(d::Dict) =
parse_from_keys(d, [:lab, :labels, :label])
# -----------------------------

# -----------------------------
parse_B(d::Dict) =
parse_from_keys(d, [:B, :frame, :axis])
# -----------------------------

function kwarg_default(d::Dict, func::Function, default_value)
	val, d = func(d)
	if val==nothing; val=default_value; end
	return val, d
end


text:
    function mime_type()::String
        executor(data["faker"]["file"]["mime_type"])
    end
    #param category: application|audio|image|message|model|multipart|text|video

    function file_extension()::String
        executor(data["faker"]["file"]["extension"])
    end

    function file_name(extencion::String="None")::String

        if extencion == "None"
            extencion = file_extension()
        end

        return string(lexify("????"),".",extencion)

    end


text:
using Random
using MinimalRLCore

"""
    Simple reward structure for goals. Can be duck typed.
"""
struct BasicRewFunc
    value::Float64
end

is_terminal(rew::BasicRewFunc) = true
get_reward(rew::BasicRewFunc) = rew.value

module ContGridWorldParams
UP = 1
RIGHT = 2
DOWN = 3
LEFT = 4

ACTIONS = [UP, RIGHT, DOWN, LEFT]

const BODY_RADIUS = 0.2
const STEP = 0.5

const AGENT_BOUNDRIES = [
    [BODY_RADIUS, 0.0],
    [-BODY_RADIUS, 0.0],
    [0.0, BODY_RADIUS],
    [0.0, -BODY_RADIUS]
]
end

"""
    ContGridWorld
A Continuous grid world domain with generic dynamics to allow for different types of mazes.
    - state: [y, x]
# Args
"""
mutable struct ContGridWorld{RF, SF} <: AbstractEnvironment
    # State is stored as (y,x)
    state::Array{Float64, 1}
    walls::Array{Bool, 2}

    # goals stored as (y,x)
    goals::Array{Int, 2}
    # reward_funcs::Dict{Int, F}
    cumulant_schedule::RF
    start_func::SF

    max_action_noise::Float64
    drift_noise::Float64
    normalized::Bool
    collision::Bool

    edge_locs::Array{Array{Int64, 1}, 1}
    collision_check::Array{Bool, 1}
    new_state::Array{Float64, 1}

    per_step_reward::Float64

    ContGridWorld(walls, goals,
                  cumulant_schedule::RF,# rew_funcs::Dict{Int, F},
                  start_func::SF, max_action_noise,
                  drift_noise, normalized, per_step_rew=0.0) where {RF, SF} =
                      new{RF, SF}([0.0, 0.0], walls, goals,
                                  cumulant_schedule,
                                  start_func,
                                  max_action_noise, drift_noise,
                                  normalized, false,
                                  [[0,0] for i in 1:4],
                                  fill(false, 4),
                                  zeros(2),
                                  per_step_rew)
end

ContGridWorld(walls, goals,
              # rew_funcs,
              max_action_noise,
              drift_noise, normalized, per_step_rew = 0.0) =
                  ContGridWorld(walls, goals,
                                nothing,
                                max_action_noise,
                                drift_noise,
                                normalized,
                                per_step_rew = 0.0)

MinimalRLCore.is_terminal(env::ContGridWorld) = begin
    goal_id = which_goal(env)
    if env.cumulant_schedule isa Dict
        if goal_id ∈ keys(env.cumulant_schedule)
            is_terminal(env.cumulant_schedule[goal_id])
        else
            false
        end
    elseif env.cumulant_schedule isa CumulantSchedule
        which_goal(env) > 0
    elseif env.cumulant_schedule isa nothing
        false
    end
end

MinimalRLCore.get_reward(env::ContGridWorld) = begin
    goal_id = which_goal(env)
    if env.cumulant_schedule isa Dict
        if goal_id ∈ keys(env.cumulant_schedule)
            get_reward(env.cumulant_schedule[goal_id])
        else
            env.per_step_reward
        end
    elseif env.cumulant_schedule isa CumulantSchedule
        env.per_step_reward
    elseif env.cumulant_schedule isa nothing
        env.per_step_reward
    end
end

function get_cumulants(env::ContGridWorld, cs::CumulantSchedule)
    goal_id = which_goal(env)
    cumulants = zeros(4)
    if goal_id == 1 #check_goal(env, 1)
        cumulants[1] = get_cumulant(cs, "G1")
    elseif goal_id == 2 #check_goal(env, 2)
        cumulants[2] = get_cumulant(cs, "G2")
    elseif goal_id == 3 #check_goal(env, 3)
        cumulants[3] = get_cumulant(cs, "G3")
    elseif goal_id == 4 #check_goal(env, 4)
        cumulants[4] = get_cumulant(cs, "G4")
    end
    cumulants
end

MinimalRLCore.get_state(env::ContGridWorld) = begin
    s = env.normalized ? env.state./size(env) : env.state
    if env.cumulant_schedule isa CumulantSchedule
        s = vcat(s, get_cumulants(env, env.cumulant_schedule))
    end
    s
end

random_state(env::ContGridWorld, rng) = [rand(rng)*size(env.walls)[1], rand(rng)*size(env.walls)[2]]
function random_start_state(env::ContGridWorld, rng)
    state = random_state(env, rng)
    while is_wall(env, state) || which_goal(env, state) > 0
        state = random_state(env, rng)
    end
    return state
end
Base.size(env::ContGridWorld) = size(env.walls)
num_actions(env::ContGridWorld) = 4
get_states(env::ContGridWorld) = findall(x->x==false, env.walls)
MinimalRLCore.get_actions(env::ContGridWorld) = ContGridWorldParams.ACTIONS

valid_projected_state(env::ContGridWorld, prj) = begin
    gw_size = size(env.walls)
    prj[1] >= 1 && prj[2] >= 1 && prj[1] <= gw_size[1]  && prj[2] <= gw_size[2]
end

project(env::ContGridWorld, state) = [Int64(floor(state[1])) + 1, Int64(floor(state[2])) + 1]
project(env::ContGridWorld, state, loc) = begin; loc[1] = Int64(floor(state[1]) + 1); loc[2] = Int64(floor(state[2]) + 1); end;

function is_wall(env::ContGridWorld, state::Array{Float64, 1})
    prj = project(env, state)
    return is_wall(env, prj)
end

function is_wall(env::ContGridWorld, prj::Array{Int64, 1})
    gw_size = size(env.walls)
    if !valid_projected_state(env, prj)
        return false
    end
    env.walls[prj[1], prj[2]]
end

which_goal(env::ContGridWorld) = which_goal(env, env.state)

function which_goal(env::ContGridWorld, state::Array{Float64, 1})
    prj = project(env, state)
    if !valid_projected_state(env, prj)
        return -1
    end
    env.goals[prj[1], prj[2]]
end

function handle_collision(env::ContGridWorld, state, action)

    #Approximate by the square...
    CRP = ContGridWorldParams
    new_state = copy(state)

    new_state[1] = clamp(new_state[1], CRP.BODY_RADIUS, size(env.walls)[1] - CRP.BODY_RADIUS)
    new_state[2] = clamp(new_state[2], CRP.BODY_RADIUS, size(env.walls)[2] - CRP.BODY_RADIUS)

    # Really basic collision detection for 2-d plane worlds.
    collided = new_state[1] != state[1] || new_state[2] != state[2]

    for i in 1:4
        project(env::ContGridWorld, new_state .+ CRP.AGENT_BOUNDRIES[i], env.edge_locs[i])
        env.collision_check[i] = is_wall(env, env.edge_locs[i])
    end

    collided = collided || any(env.collision_check)

    if env.collision_check[1] && env.collision_check[2]
        wall_piece = project(env, new_state)
        if action == CRP.DOWN
            new_state[1] = (wall_piece[1]) + CRP.BODY_RADIUS
        elseif action == CRP.UP
            new_state[1] = (wall_piece[1] - 1) - CRP.BODY_RADIUS
        end
    elseif env.collision_check[1]
        new_state[1] = (env.edge_locs[1][1] - 1) - CRP.BODY_RADIUS
    elseif env.collision_check[2]
        new_state[1] = (env.edge_locs[2][1]) + CRP.BODY_RADIUS
    end

    if env.collision_check[3] && env.collision_check[4]
        wall_piece = project(env, new_state)
        if action == CRP.LEFT
            new_state[2] = (wall_piece[2]) + CRP.BODY_RADIUS
        elseif action == CRP.RIGHT
            new_state[2] = (wall_piece[2] - 1) - CRP.BODY_RADIUS
        end
    elseif env.collision_check[3]
        new_state[2] = (env.edge_locs[3][2] - 1) - CRP.BODY_RADIUS
    elseif env.collision_check[4]
        new_state[2] = (env.edge_locs[4][2]) + CRP.BODY_RADIUS
    end

    new_state[1] = clamp(new_state[1], CRP.BODY_RADIUS, size(env.walls)[1] - CRP.BODY_RADIUS)
    new_state[2] = clamp(new_state[2], CRP.BODY_RADIUS, size(env.walls)[2] - CRP.BODY_RADIUS)

    return new_state, collided
end

function handle_collision!(env::ContGridWorld, action)
    env.state, collided = handle_collision(env, env.state, action)
    return collided
end

function MinimalRLCore.reset!(env::ContGridWorld, rng::Random.AbstractRNG=Random.GLOBAL_RNG)
    if env.start_func isa Nothing
        env.state = random_start_state(env, rng)
        env.state
    else
        env.state = env.start_func(env, rng)
        env.state
    end
end

function MinimalRLCore.reset!(env::ContGridWorld, state::AbstractArray)
    s = env.normalized ? state.*size(env) : state

    if is_wall(env, s)
        throw("Cannot reset environment to invalid state: $(state)")
    end

    env.state .= s
    # @show env.state
end

function mini_step!(env::ContGridWorld, step, action)
    env.state .+= step
    collision = handle_collision!(env, action)
    return collision
end

function MinimalRLCore.environment_step!(env::ContGridWorld, action, rng=Random.GLOBAL_RNG)

    CRP = ContGridWorldParams
    next_step = zeros(2)
    @assert action ∈ CRP.ACTIONS

    if action == CRP.UP
        next_step[1] = (CRP.STEP - rand(rng)*env.max_action_noise)
        next_step[2] = randn(rng)*env.drift_noise
    elseif action == CRP.DOWN
        next_step[1] = -(CRP.STEP - rand(rng)*env.max_action_noise)
        next_step[2] = randn(rng)*env.drift_noise
    elseif action == CRP.RIGHT
        next_step[1] = randn(rng)*env.drift_noise
        next_step[2] = CRP.STEP - rand(rng)*env.max_action_noise
    elseif action == CRP.LEFT
        next_step[1] = randn(rng)*env.drift_noise
        next_step[2] = -(CRP.STEP - rand(rng)*env.max_action_noise)
    end

    # mini_physics simulation for 1 second (== 10 steps of 0.1 seconds)
    Δt=1.0
    τ=10
    next_step .*= Δt/τ
    collision = false
    for t in 1:τ
        collision = mini_step!(env, next_step, action)
        if collision
            break
        end
    end
    env.collision = collision

    update!(env.cumulant_schedule, env.state)
end

module ContGridWorldPlotParams

using Colors
import ColorSchemes

const SIZE = 10
const BG = Colors.RGB(1.0, 1.0, 1.0)
const WALL = Colors.RGB(0.3, 0.3, 0.3)
const AC = Colors.RGB(0.69921875, 0.10546875, 0.10546875)
const GOAL = Colors.RGB(0.796875, 0.984375, 0.76953125)
const GOAL_PALETTE = ColorSchemes.tol_muted
const AGENT = [AC AC AC AC;
               AC AC AC AC;
               AC AC AC AC;
               AC AC AC AC;]

end

@recipe function f(env::ContGridWorld)
    ticks := nothing
    foreground_color_border := nothing
    grid := false
    legend := false
    aspect_ratio := 1
    xaxis := false
    yaxis := false

    PP = ContGridWorldPlotParams

    s = size(env)
    screen = fill(PP.BG, (s[2] + 2)*PP.SIZE, (s[1] + 2)*PP.SIZE)
    screen[:, 1:PP.SIZE] .= PP.WALL
    screen[1:PP.SIZE, :] .= PP.WALL
    screen[end-(PP.SIZE-1):end, :] .= PP.WALL
    screen[:, end-(PP.SIZE-1):end] .= PP.WALL

    for i ∈ 1:s[1]
        for j ∈ 1:s[2]
            sqr_i = ((i)*PP.SIZE + 1):((i+1)*PP.SIZE)
            sqr_j = ((j)*PP.SIZE + 1):((j+1)*PP.SIZE)
            if env.walls[s[1] - i + 1, j]
                screen[sqr_i, sqr_j] .= PP.WALL
            end
            if env.goals[s[1] - i + 1, j] > 0
                screen[sqr_i, sqr_j] .=
                    PP.GOAL_PALETTE[env.goals[s[1] - i + 1, j]]
            end
        end
    end
    state = ((env.state[1] - 1, env.state[2] + 1).*PP.SIZE)

    screen[collect(-1:2) .+ Int(floor(PP.SIZE*s[1] - state[1])), collect(-1:2) .+ Int(floor(state[2]))] .= PP.AGENT

    screen
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function forward(backend::GPUBackend, state::ConcatLayerState, inputs::Vector{Blob})
  output = state.blobs[1]
  shift = 0
  for i = 1:length(inputs)
    copy_to_shifted!(backend, output, inputs[i], (state.layer.dim, shift))
    shift += size(inputs[i], state.layer.dim)
  end
end

function backward(backend::GPUBackend, state::ConcatLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  grad = state.blobs_diff[1]
  shift = 0
  for i = 1:length(diffs)
    if !isa(diffs[i], NullBlob)
      copy_from_shifted!(backend, diffs[i], grad, (state.layer.dim, shift))
    end
    shift += size(inputs[i], state.layer.dim)
  end
end


text:
using Test
using LinearAlgebra
using DataFrames
# Preprocess Tests:
include("preprocesstests.jl")
# Stats Tests:

# Models Tets:


text:
"""
Author: Eric Ford
Created: September 2020
"""

#module FitTemplateToCCF

#import ..RVFromCCF: AbstractMeasureRvFromCCF, MeasureRvFromCCFQuadratic
# packages
using LinearAlgebra

# TODO: export once passes test that it's working
#export MeasureRvFromCCFTemplate

"""  Functor to estimate RV based on fitting a template and its derivative near minimum of the CCF. """
struct MeasureRvFromCCFTemplate <: AbstractMeasureRvFromCCF
    v_grid::AbstractVector{Float64}
    template::AbstractVector{Float64}
    deriv::AbstractVector{Float64}
	mean_var::Real
    v_idx_to_fit::UnitRange{Int64}
end

# Todo: Replace with version in RvSpectMLBase or a (brute force?/temporal?) GP?
function numerical_deriv( x::AA1, y::AA2 )   where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}  }
	@assert length(x) == length(y)
	dydx = zeros(size(y))
	dydx[1] = (y[2]-y[1])/(x[2]-x[1])
	dydx[2:end-1] .= (y[3:end].-y[1:end-2])./(x[3:end].-x[1:end-2])
	dydx[end] = (y[end]-y[end-1])/(x[end]-x[end-1])
	return dydx
end

"""
Construct functor to estimate RV based on the CCF.
TODO: Implement correctly.  Not yet working/tested.
Optional Arguments:
- `v_grid`:  list of velocities where template is evaluated
- `tempalte`: template ccf evaluated at v_grid
- `frac_of_width_to_fit`: (0.5)
- `measure_width_at_frac_depth`: (0.5)
"""
function MeasureRvFromCCFTemplate(; v_grid::AbstractVector{T1},
                                    template::AbstractVector{T2},
                                    frac_of_width_to_fit::Real = default_frac_of_width_to_fit,
                                    measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth,
									mean_var = 0
                                     ) where { T1<:Real, T2<:Real }
    @assert length(v_grid) == length(template)
    @assert length(v_grid) >= 3
    @assert 0.25 <= measure_width_at_frac_depth <= 0.75
    @assert 0.1 <= frac_of_width_to_fit <= 5.0  # TODO: FIgure out appropriate range
	@assert 0 <= mean_var < Inf
	# Estimate derivative
    deriv = numerical_deriv(v_grid, template)
	# find the min and fit only the part near the minimum of the CCF
    v_min, v_idx_to_fit = find_idx_at_and_around_minimum(v_grid, template, frac_of_width_to_fit=frac_of_width_to_fit, measure_width_at_frac_depth=measure_width_at_frac_depth)
    MeasureRvFromCCFTemplate(v_grid, template, deriv, mean_var, v_idx_to_fit)
end

function set_mean_var(in::MeasureRvFromCCFTemplate, mean_var::Real)
	MeasureRvFromCCFTemplate(v_grid=in.v_grid,template=in.template, frac_of_width_to_fit=in.frac_of_width_to_fit,
		measure_width_at_frac_depth=in.measure_width_at_frac_depth, mean_var=mean_var)
end

function (mrv::MeasureRvFromCCFTemplate)(vels::A1, ccf::A2, ccf_var::A3 = ones(length(ccf)) ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1}  }
        @assert all(vels .== mrv.v_grid)  # Could relax this, but keep it simple for now
        # fit only the part near the minimum of the CCF
		deriv = view(mrv.deriv,mrv.v_idx_to_fit)
		var = view(ccf_var,mrv.v_idx_to_fit)
		# Minus sign due to weird convention of defining CCF to look like an absorption line (e.g., 1 at Inf minimum at best fit velocity)
		rv = -sum( (view(ccf,mrv.v_idx_to_fit) .- view(mrv.template,mrv.v_idx_to_fit) ) .* deriv ./ var)
        # Warning: Uncertaintiy ignores correlations between pixels, particularly problematic when oversample pixels
		denom = sum( abs2.(deriv) ./ var )
		rv *= (1/denom)
		σ_rv =  sqrt(1/denom)
		return (rv=rv, σ_rv=σ_rv)
end

function (mrv::MeasureRvFromCCFTemplate)(vels::A1, ccf::A2, ccf_var::A3, ccf_covar::A4  ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1}, T4<:Real, A4<:AbstractArray{T4,2}  }
        @assert all(vels .== mrv.v_grid)  # Could relax this, but keep it simple for now
		# fit only the part near the minimum of the CCF
		deriv = view(mrv.deriv,mrv.v_idx_to_fit)
		diag_scale_factor = mean(ccf_var)-mrv.mean_var
		covar = view(ccf_covar,mrv.v_idx_to_fit,mrv.v_idx_to_fit) + diagm(diag_scale_factor*ones(length(mrv.v_idx_to_fit)))
        # find the min and fit only the part near the minimum of the CCF
		ccf_inv_covar_times_deriv = covar \ deriv
		denom = deriv' * ccf_inv_covar_times_deriv
		# Minus sign due to weird convention of defining CCF to look like an absorption line (e.g., 1 at Inf minimum at best fit velocity)
		rv = -(view(ccf,mrv.v_idx_to_fit) .- view(mrv.template,mrv.v_idx_to_fit) )' * ccf_inv_covar_times_deriv
		rv *= (1/denom)
		σ_rv = sqrt(1/denom)
		return (rv=rv, σ_rv=σ_rv)
end