text:
#math_graphics.jl:
#-------------------------------------------------------------------------------

#==Aliases
===============================================================================#
#Real data type used for display functionnality (low-level):
const DReal = Float64

#==Constants
===============================================================================#
const DInf = convert(DReal, Inf)
const DNaN = convert(DReal, NaN)

const COLOR_TRANSPARENT = ARGB32(0, 0, 0, 0)

const COLOR_BLACK = RGB24(0, 0, 0)
const COLOR_WHITE = RGB24(1, 1, 1)

const COLOR_RED = RGB24(1, 0, 0)
const COLOR_GREEN = RGB24(0, 1, 0)
const COLOR_BLUE = RGB24(0, 0, 1)

#==Data structures.
===============================================================================#
#TODO: Use library version?
abstract type Point end #So we can use as a function
struct Point2D <: Point
	x::DReal
	y::DReal
end

abstract type DirectionalVector end
struct Vector2D <: DirectionalVector
	x::DReal
	y::DReal
end
Vector2D(p::Point2D) = Vector2D(p.x, p.y)
Point2D(p::Vector2D) = Point2D(p.x, p.y)

#Plot extents along one dimension:
struct PExtents1D
	min::DReal
	max::DReal
end
PExtents1D(;min::Real=DNaN, max::Real=DNaN) = PExtents1D(min, max)

#Could use BoundingBox, but cannot control type (and is technically not a BB).
#TODO: PExtents2D from PExtents1D?
struct PExtents2D
	xmin::DReal
	xmax::DReal
	ymin::DReal
	ymax::DReal
end
PExtents2D(;xmin::Real=DNaN, xmax::Real=DNaN, ymin::Real=DNaN, ymax::Real=DNaN) =
	PExtents2D(xmin, xmax, ymin, ymax)
PExtents2D(xext::PExtents1D, yext::PExtents1D) =
	PExtents2D(xext.min, xext.max, yext.min, yext.max)

abstract type LinearTransform end

#Linear 1D transform
#	v' = (v + v0) * vs
#NOTE: Format used instead of v'*vs + v0 for improved numeric stability?
struct LXF1D <: LinearTransform
	vs::DReal #scaling factor
	v0::DReal #offset
end

#2D linear transform (only scale/offset - no rotation).
#	x' = (x + x0) * xs
#	y' = (y + y0) * ys
#where: xu/yu: user coordinates & xd/yd: device coordinates
#NOTE: Format used instead of x'*xs + x0 for improved numeric stability?
#TODO: Does struct cause copy on function calls instead of using pointers??
#TODO: Maybe rename gvien
struct Transform2D <: LinearTransform
	xs::DReal
	x0::DReal
	ys::DReal
	y0::DReal
end

struct LineStyle
	style::Symbol
	width::Float64 #Device units (~pixels)
	color::Colorant
end

#==Constructor-like functions
===============================================================================#
#=Inputs
 -ext: Visible plot extents
 -inputb: Bounding box of display area.
NOTE:
 -Assumes device y-coordinate increases as we descend on the screen.
=#
#TODO: rename Transform2D->transform (result depends on input params)
function Transform2D(ext::PExtents2D, inputb::BoundingBox)
	xs = (inputb.xmax-inputb.xmin) / (ext.xmax-ext.xmin)
	ys = -(inputb.ymax-inputb.ymin) / (ext.ymax-ext.ymin)
	x0 = + inputb.xmin/xs - ext.xmin
	y0 = + inputb.ymax/ys - ext.ymin

	return Transform2D(xs, x0, ys, y0)
end

function LinearTransform(ext::PExtents1D, vmin::DReal, vmax::DReal)
	vs = (vmax-vmin) / (ext.max-ext.min)
	v0 = + vmin/vs - ext.min
	return LXF1D(vs, v0)
end

#==Accessors
===============================================================================#
xvalues(ext::PExtents2D) = PExtents1D(ext.xmin, ext.xmax)
yvalues(ext::PExtents2D) = PExtents1D(ext.ymin, ext.ymax)

#==Basic operations
===============================================================================#
#Compute vector norm:
function vecnorm(v::Vector2D)
	return sqrt(v.x*v.x+v.y*v.y)
end

#Simple scaling - no commutativity:
Base.:*(s::DReal, pt::Point2D) = Point2D(s*pt.x, s*pt.y)
Base.:*(s::Real, pt::Point2D) = DReal(s)*pt
Base.:+(p1::Point2D, p2::Point2D) = Point2D(p1.x+p2.x, p1.y+p2.y)
Base.:-(p1::Point2D, p2::Point2D) = Point2D(p1.x-p2.x, p1.y-p2.y)

Base.:*(s::DReal, v::Vector2D) = Vector2D(s*v.x, s*v.y)
Base.:*(s::Real, v::Vector2D) = DReal(s)*v
Base.:+(v1::Vector2D, v2::Vector2D) = Vector2D(v1.x+v2.x, v1.y+v2.y)

Base.:+(p::Point2D, v::Vector2D) = Point2D(p.x+v.x, p.y+v.y)
Base.:-(p::Point2D, v::Vector2D) = Point2D(p.x-v.x, p.y-v.y)
Base.:+(v::Vector2D, p::Point2D) = p+v

function union(e1::PExtents1D, e2::PExtents1D)
	(_min, _ign) = extents_finite([e1.min, e2.min])
	(_ign, _max) = extents_finite([e1.max, e2.max])
	return PExtents1D(_min, _max)
end
function union(elist::Vector{PExtents1D})
	result = PExtents1D(DNaN, DNaN)
	for ext in elist
		result = union(result, ext)
	end
	return result
end

function union(e1::PExtents2D, e2::PExtents2D)
	(xmin, _ign) = extents_finite([e1.xmin, e2.xmin])
	(_ign, xmax) = extents_finite([e1.xmax, e2.xmax])
	(ymin, _ign) = extents_finite([e1.ymin, e2.ymin])
	(_ign, ymax) = extents_finite([e1.ymax, e2.ymax])

	return PExtents2D(xmin, xmax, ymin, ymax)
end
function union(elist::Vector{PExtents2D})
	result = PExtents2D(DNaN, DNaN, DNaN, DNaN)
	for ext in elist
		result = union(result, ext)
	end
	return result
end

#Overwrite with new extents, if defined
#-------------------------------------------------------------------------------
function Base.merge(base::PExtents1D, new::PExtents1D)
	baseifnan(bv, nv) = isnan(nv) ? bv : nv
	emin = baseifnan(base.min, new.min)
	emax = baseifnan(base.max, new.max)
	return PExtents1D(emin, emax)
end

function Base.merge(base::PExtents2D, new::PExtents2D)
	baseifnan(bv, nv) = isnan(nv) ? bv : nv
	xmin = baseifnan(base.xmin, new.xmin)
	xmax = baseifnan(base.xmax, new.xmax)
	ymin = baseifnan(base.ymin, new.ymin)
	ymax = baseifnan(base.ymax, new.ymax)
	return PExtents2D(xmin, xmax, ymin, ymax)
end

Base.isfinite(ext::PExtents1D) = isfinite(ext.min) && isfinite(ext.max)

Base.isfinite(ext::PExtents2D) =
	isfinite(ext.xmin) && isfinite(ext.xmax) &&
	isfinite(ext.ymin) && isfinite(ext.ymax)

#==Transform computations
===============================================================================#

apply(xf::LXF1D, v::DReal) = (v + xf.v0)*xf.vs

function apply(xf::Transform2D, pt::Point2D)
	x = (pt.x + xf.x0)*xf.xs
	y = (pt.y + xf.y0)*xf.ys
	return Point2D(x, y)
end
function apply(xf::Transform2D, pt::Vector2D)
	x = pt.x*xf.xs
	y = pt.y*xf.ys
	return Vector2D(x, y)
end

#Compute inverse transformation
#-------------------------------------------------------------------------------
function apply_inv(xf::Transform2D, pt::Point2D)
	x = pt.x/xf.xs - xf.x0
	y = pt.y/xf.ys - xf.y0
	return Point2D(x, y)
end
function apply_inv(xf::Transform2D, pt::Vector2D)
	x = pt.x/xf.xs
	y = pt.y/xf.ys
	return Vector2D(x, y)
end

#==Interpolation functions
===============================================================================#

#Interpolate between two points.
#-------------------------------------------------------------------------------
function interpolate(p1::Point2D, p2::Point2D, x::DReal)
	m = (p2.y-p1.y) / (p2.x-p1.x)
	return m*(x-p1.x)+p1.y
end

#Last line


text:
module TestCLBLAS
using Base.Test

using OpenCL: cl
import OpenCL.cl.CLArray
using CLBLAS

CLBLAS.setup()
dev, ctx, q = cl.create_compute_context()
const testtypes = (Float32, Float64)
include("test_l1.jl")
include("test_l2.jl")
include("test_l3.jl")

end

ext: jl


text:
# Module written by Ed Scheinerman, ers@jhu.edu
# distributed under terms of the MIT license

module ShowSet

import Base: string, show, AbstractSet
import Base.Multimedia.display

export set_empty, get_empty

EMPTY_SET = "Ø"  # this is LaTeX \O

"""
The string form of the empty set is initially set to `"Ø"` (LaTeX `\\O`).
To change it, use `set_empty(str)`.
Suggested alternatives are `"{}"` or `"∅"` (LaTeX `\\emptyset`).

See `get_empty`.
"""
function set_empty(str::String)::String
    global EMPTY_SET = str
end

"""
`get_empty()` returns the current string that represents the empty set.

See `set_empty`.
"""
function get_empty()::String
    global EMPTY_SET
    return EMPTY_SET
end

function string(A::AbstractSet)::String
    elements = collect(A)

    if length(elements) == 0
        return get_empty()
    end

    try
        sort!(elements)
    catch
        p = sortperm(hash.(elements))
        elements = elements[p]
    end
    return "{" * join(elements,",") * "}"
end

show(io::IO, A::Set)    = print(io,string(A))
show(io::IO, A::BitSet) = print(io,string(A))
display(A::Set) = print(string(A))
display(A::BitSet) = print(string(A))

end # module ShowsSet


text:
#!/usr/bin/env julia

# Draw lines with various dash styles.

using Compose
using Color

function draw_lines(joins)
    if length(joins) == 0
        canvas()
    else
        compose(compose(canvas(), lines((5mm, 5mm), (10mm, 10mm), (15mm, 5mm)), strokelinejoin(joins[1])),
                compose(canvas(0, 5mm), draw_lines(joins[2:end])))
    end
end

joins = [LineJoinRound(), LineJoinMiter(), LineJoinBevel()]
c = draw_lines(joins)

imgs = [SVG("linejoins.svg", 2cm, 3cm),
        PDF("linejoins.pdf", 2cm, 3cm),
        D3("linejoins.js", 2cm, 3cm)]
for img = imgs
    draw(img, compose(c, stroke(color("black")), linewidth(2mm), fill(Nothing())))
end


text:
@testset "canny" begin
    @info "Test: Canny"

    @testset "API" begin
        img_gray = Gray{N0f8}.(load("algorithms/References/circle.png"))
        img = copy(img_gray)

        # Detect edges
        f = Canny()
        edges_img_1 = detect_edges(img, f)
        @test img == img_gray # img unchanged
        @test eltype(edges_img_1) == Gray{N0f8}

        edges_img_2 = detect_edges(Gray{Bool}, img, f)
        @test img == img_gray # img unchanged
        @test eltype(edges_img_2) == Gray{Bool}

        edges_img_3 = similar(img, Bool)
        detect_edges!(edges_img_3, img, f)
        @test img == img_gray # img unchanged
        @test eltype(edges_img_3) == Bool

        edges_img_4 = copy(img_gray)
        detect_edges!(edges_img_4, f)
        @test eltype(edges_img_4) == Gray{N0f8}

        @test edges_img_1 == edges_img_2
        @test edges_img_1 == edges_img_3
        @test edges_img_1 == edges_img_4

        for T in generate_test_types([Float32, N0f8, Bool], [Gray])
            @test eltype(detect_edges(T, img, f)) == T
        end

        # Detect subpixel edges
        img_gray = Gray{N0f8}.(load("algorithms/References/circle.png"))
        img = copy(img_gray)
        f = Canny(thinning_algorithm = SubpixelNonmaximaSuppression())
        edges_img_1, offsets_1 = detect_subpixel_edges(img, f)
        @test img == img_gray # img unchanged
        @test eltype(edges_img_1) == Gray{N0f8}
        @test eltype(offsets_1) == SVector{2, Float64}

        edges_img_2, offsets_2 = detect_subpixel_edges(Gray{Bool}, img, f)
        @test img == img_gray # img unchanged
        @test eltype(edges_img_2) == Gray{Bool}
        @test eltype(offsets_2) == SVector{2, Float64}

        edges_img_3, offsets_3 = detect_subpixel_edges(Gray{Bool}, Float32, img, f)
        @test img == img_gray # img unchanged
        @test eltype(edges_img_3) == Gray{Bool}
        @test eltype(offsets_3) == SVector{2, Float32}

        edges_img_4 = similar(img, Bool)
        offsets_4 = similar(offsets_3)
        detect_subpixel_edges!(edges_img_4, offsets_4, img, f)
        @test img == img_gray # img unchanged
        @test eltype(edges_img_4) == Bool
        @test offsets_4  == offsets_3
        @test eltype(offsets_4) == SVector{2, Float32}

        edges_img_5 = copy(img_gray)
        offsets_5 = similar(offsets_1)
        detect_subpixel_edges!(edges_img_5, offsets_5, img, f)
        @test eltype(edges_img_5) == Gray{N0f8}
        @test eltype(offsets_5) == SVector{2, Float64}

        @test edges_img_1 == edges_img_2
        @test edges_img_1 == edges_img_3
        @test edges_img_1 == edges_img_4
        @test edges_img_1 == edges_img_5

        @test offsets_1  == offsets_2
        @test offsets_1  ≈ offsets_3
        @test offsets_1  ≈ offsets_4
        @test offsets_1  == offsets_5
    end

    @testset "Offset Arrays" begin
        img_gray = Gray{N0f8}.(load("algorithms/References/circle.png"))
        img_gray_offset = OffsetArray(img_gray, -25:25, -25:25)

        f = Canny()
        edges_img_1 = detect_edges(img_gray_offset, f)
        @test eltype(edges_img_1) == Gray{N0f8}
        @test axes(edges_img_1) == axes(img_gray_offset)

        g = Canny(thinning_algorithm = SubpixelNonmaximaSuppression())
        edges_img_2, offsets = detect_subpixel_edges(img_gray_offset, g)

        @test edges_img_1 == edges_img_2
        @test axes(edges_img_2) == axes(img_gray_offset)
        @test axes(offsets) == axes(img_gray_offset)
    end

    @testset "Keywords" begin
        img_gray = Gray{N0f8}.(load("algorithms/References/circle.png"))
        img = copy(img_gray)
        low = [0.053374808162753785, Percentile(80)]
        high = [0.16280777841581243, Percentile(90)]
        spatial_scale = [1.4, 1.4]
        # Detect edges
        for i = 1:2
            f = Canny(spatial_scale = spatial_scale[i], low = low[i], high = high[i])
            @test_reference "References/circle_edge.png" Gray.(detect_edges(img, f)) by=edge_detection_equality()
            if i == 2
                @test_reference "References/circle_edge.png" Gray.(detect_edges(img * 0.5, f)) by=edge_detection_equality()
            end
        end

        # Detect subpixel edges
        for i = 1:2
            g = Canny(spatial_scale = spatial_scale[i], low = low[i], high = high[i], thinning_algorithm = SubpixelNonmaximaSuppression())
            out1, offsets1 = detect_subpixel_edges(img, g)
            out2, offsets2 = detect_subpixel_edges(img * 0.5, g)
            @test_reference "References/circle_edge.png" Gray.(out1) by=edge_detection_equality()
            if i == 2
                @test_reference "References/circle_edge.png" Gray.(out2) by=edge_detection_equality()
            end
        end
    end

    @testset "Types" begin
        # Gray
        img_gray = Gray{N0f8}.(load("algorithms/References/circle.png"))
        f = Canny(low = Percentile(80), high = Percentile(90),
                  spatial_scale = 1.4)
        g = Canny(low = Percentile(80), high = Percentile(90),
                  spatial_scale = 1.4,
                  thinning_algorithm = SubpixelNonmaximaSuppression())

        type_list = generate_test_types([Float32, N0f8], [Gray])
        for T in type_list
            img = T.(img_gray)
            out1 = Gray.(detect_edges(img, f))
            out2, offsets2 = detect_subpixel_edges(img, g)
            @test_reference "References/circle_edge.png" out1 by=edge_detection_equality()
            @test_reference "References/circle_edge.png" Gray.(out2) by=edge_detection_equality()
        end

        # Color3
        img_color = RGB{Float64}.(load("algorithms/References/circle.png"))
        f = Canny(low = Percentile(80), high = Percentile(90),
                  spatial_scale = 1.4)
        g = Canny(low = Percentile(80), high = Percentile(90),
                  spatial_scale = 1.4,
                  thinning_algorithm = SubpixelNonmaximaSuppression())

        type_list = generate_test_types([Float32, N0f8], [RGB, Lab])
        for T in type_list
            img = T.(img_color)
            out1 = Gray.(detect_edges(img, f))
            out2, offsets2 = detect_subpixel_edges(img, g)
            @test_reference "References/circle_edge.png" out1 by=edge_detection_equality()
            @test_reference "References/circle_edge.png" Gray.(out2) by=edge_detection_equality()
        end
    end

    @testset "Default Values" begin
        img = Gray{N0f8}.(testimage("cameraman"))
        out, offsets = detect_subpixel_edges(img)
        @test_reference "References/cameraman_edge.png" Gray.(detect_edges(img)) by=edge_detection_equality()
        @test_reference "References/cameraman_edge.png" Gray.(out) by=edge_detection_equality()
    end

    @testset "Numerical" begin
        # Check that the image only has ones or zeros.
        img = Gray{N0f8}.(load("algorithms/References/circle.png"))
        f = Canny()
        img₀₁ = detect_edges(img, f)
        non_zeros = findall(x -> x != 0.0 && x != 1.0, img₀₁)
        @test length(non_zeros) == 0
    end

    @testset "Subpixel Accuracy on Circle Image" begin
        # Equation of circle (x-a)^2 + (y - b)^2 = r^2 and corresponding image.
        a = 26
        b = 26
        r = 15
        img = Gray{N0f8}.(load("algorithms/References/circle.png"))
        algo = Canny(spatial_scale = 1.4,
                     thinning_algorithm = SubpixelNonmaximaSuppression(),
                     low = Percentile(75), high = Percentile(80))
        nms, offsets = detect_subpixel_edges(img, algo)

        # Verify that the subpixel coordinates more accurately satisfy the
        # circle equation.
        total₁ = 0.0
        total₂ = 0.0
        N = sum(nms)
        for i in CartesianIndices(nms)
            y,x = i.I
            δy, δx = offsets[i]
            y′ = y + δy
            x′ = x + δx
            if nms[i] == 1
                # Degree to which pixel coordinates fail to satisfy circle
                # equation.
                val₁ = abs((x-a)^2 + (y - b)^2 - r^2)
                total₁ = total₁ + val₁

                # Degree to which subpixel coordinates fail to satisfy circle
                # equation.
                val₂ = abs((x′-a)^2 + (y′ - b)^2 - r^2)
                total₂ = total₂ + val₂
            end
        end
        # The subpixel coordinates yield a better fit.
        @test total₂ <  total₁
        @test total₂ / N < 5.54
    end

    @testset "Subpixel Accuracy on Synthetic Image" begin
        # Test image for which we want to detect edges.
        img  = [0.5  0.5   0.7  0.6  0.6
                0.5  0.5   0.7  0.6  0.6
                0.5  0.5   0.7  0.6  0.6
                0.5  0.5   0.7  0.6  0.6
                0.5  0.5   0.7  0.6  0.6]

        thinning_algorithm = SubpixelNonmaximaSuppression()

        # Calculate the gradient vector at each position of the filtered image.
        # The derivatives are taken with respect to the first and second
        # dimension (rows first, then columns).
        g₁, g₂ = imgradients(img, KernelFactors.scharr)
        # Gradient magnitude
        mag = hypot.(g₁, g₂)
        # Isolate local maxima of gradient magnitude by “non-maximum
        # suppression” along the local gradient direction and store location of
        # edge to subpixel precision.
        nms, offsets = thin_subpixel_edges(mag, g₁, g₂, thinning_algorithm)

        # This is the expected result for the non-maxima suppression.
        target_nms = [0.0  0.1   0.0  0.05  0.0
                      0.0  0.1   0.0  0.05  0.0
                      0.0  0.1   0.0  0.05  0.0
                      0.0  0.1   0.0  0.05  0.0
                      0.0  0.1   0.0  0.05  0.0]

        for i in CartesianIndices(target_nms)
            @test nms[i] ≈ target_nms[i]
        end

        # Verify the correctness of the subpixel refinement.
        for r = 1:5
            @test all(offsets[r, 1] .≈ SVector(0.0, 0.0))
            @test all(offsets[r, 2] .≈ SVector(0.0, 0.16666666666666666))
            @test all(offsets[r, 3] .≈ SVector(0.0, 0.0))
            @test all(offsets[r, 4] .≈ SVector(0.0, -0.5))
            @test all(offsets[r, 5] .≈ SVector(0.0, 0.0))
        end

        # Verify correctness on a 90 degree rotated variant of the synthetic
        # image
        img = img'
        # Calculate the gradient vector at each position of the filtered image.
        # The derivatives are taken with respect to the first and second
        # dimension (rows first, then columns).
        g₁, g₂ = imgradients(img, KernelFactors.scharr)
        # Gradient magnitude
        mag = hypot.(g₁, g₂)
        # Isolate local maxima of gradient magnitude by “non-maximum
        # suppression” along the local gradient direction and store location of
        # edge to subpixel precision.
        nms, offsets = thin_subpixel_edges(mag, g₁, g₂, thinning_algorithm)
        target_nms = target_nms'
        for i in CartesianIndices(target_nms)
            @test nms[i] ≈ target_nms[i]
        end
        for c = 1:5
            @test all(offsets[1, c] .≈ SVector(0.0, 0.0))
            @test all(offsets[2, c] .≈ SVector(0.16666666666666666, 0.0))
            @test all(offsets[3, c] .≈ SVector(0.0, 0.0))
            @test all(offsets[4, c] .≈ SVector(-0.5, 0.0))
            @test all(offsets[5, c] .≈ SVector(0.0, 0.0))
        end
    end

    @testset "NaNs" begin
        img = Gray.(ones(15, 15)) .* NaN
        img[4:12, 4:12] .= 0
        img[6:10, 6:10] .= 1
        img[7:9, 7:9] .= NaN

        f = Canny(spatial_scale = 1,
                  high = ImageEdgeDetection.Percentile(80),
                  low = ImageEdgeDetection.Percentile(20))
        out = detect_edges(img, f)
        @test_reference "References/edges_from_image_with_nan.png" Gray.(out) by=edge_detection_equality()

        f = Canny(spatial_scale = 1,
                  high = ImageEdgeDetection.Percentile(80),
                  low = ImageEdgeDetection.Percentile(20),
                   thinning_algorithm = SubpixelNonmaximaSuppression())
        out, offsets  = detect_subpixel_edges(img, f)
        @test_reference "References/edges_from_image_with_nan.png" Gray.(out) by=edge_detection_equality()
    end
end


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

nelms(obj) = nelms(domain(obj))

ncoords(obj) = ncoords(domain(obj))

coordtype(obj) = coordtype(domain(obj))

coordinates!(buff, obj, ind) =
  coordinates!(buff, domain(obj), ind)

"""
    coordinates!(buff, object, inds)

Non-allocating version of [`coordinates`](@ref)
"""
function coordinates!(buff, obj, inds::AbstractVector{Int})
  for j in 1:length(inds)
    coordinates!(view(buff,:,j), obj, inds[j])
  end
end

"""
    coordinates(object, ind)

Return the coordinates of the `ind` in the `object`.
"""
function coordinates(obj, ind::Int)
  N = ncoords(obj)
  T = coordtype(obj)
  x = MVector{N,T}(undef)
  coordinates!(x, obj, ind)
  x
end

"""
    coordinates(object, inds)

Return the coordinates of `inds` in the `object`.
"""
function coordinates(obj, inds::AbstractVector{Int})
  N = ncoords(obj)
  T = coordtype(obj)
  X = Matrix{T}(undef, N, length(inds))
  coordinates!(X, obj, inds)
  X
end

"""
    coordinates(object)

Return the coordinates of all indices in `object`.
"""
coordinates(obj) = coordinates(obj, 1:nelms(obj))


text:
"""
	geometric_spreading(r_ps::T, m::S, geo::GeometricSpreadingParameters, sat::NearSourceSaturationParameters) where {S<:Real, T<:Real}

Geometric spreading function, switches between different approaches on `path.geo_model`.
"""
function geometric_spreading(r_ps::T, m::S, geo::GeometricSpreadingParameters, sat::NearSourceSaturationParameters) where {S<:Real,T<:Real}
    if geo.model == :Piecewise
        return geometric_spreading_piecewise(r_ps, geo)
    elseif geo.model == :CY14
        return geometric_spreading_cy14(r_ps, geo)
    elseif geo.model == :CY14mod
        return geometric_spreading_cy14mod(r_ps, m, geo, sat)
    else
        return NaN * oneunit(get_parametric_type(geo))
    end
end

geometric_spreading(r_ps::T, m::S, path::PathParameters) where {S<:Real,T<:Real} = geometric_spreading(r_ps, m, path.geometric, path.saturation)
geometric_spreading(r_ps::T, m::S, fas::FourierParameters) where {S<:Real,T<:Real} = geometric_spreading(r_ps, m, fas.path)
# special variants that don't require a magnitude input
geometric_spreading(r_ps::T, path::PathParameters) where {T<:Real} = geometric_spreading(r_ps, NaN, path.geometric, path.saturation)
geometric_spreading(r_ps::T, fas::FourierParameters) where {T<:Real} = geometric_spreading(r_ps, fas.path)

"""
	geometric_spreading_piecewise(r_ps::V, geo::GeometricSpreadingParameters{S,T,U}) where {S<:Float64, T<:Real, U<:AbstractVector{Bool}, V<:Real}

Piecewise linear (in log-log space) geometric spreading function.
Makes use of the reference distances `Rrefi` and spreading rates `γi` in `path`.
"""
function geometric_spreading_piecewise(r_ps::V, geo::GeometricSpreadingParameters{S,T,U}) where {S<:Float64,T<:Real,U<:AbstractVector{Bool},V<:Real}
    z_r = oneunit(T)
    j = 1
    k = 1
    for i = 1:length(geo.γfree)
        @inbounds Rr0 = geo.Rrefi[i]
        @inbounds Rr1 = geo.Rrefi[i+1]
        γ_r = oneunit(T)

        if geo.γfree[i] == 0
            @inbounds γ_r *= geo.γconi[j]
            j += 1
        else
            @inbounds γ_r *= geo.γvari[k]
            k += 1
        end

        if r_ps < Rr1
            z_r *= (Rr0 / r_ps)^γ_r
            return z_r
        else
            z_r *= (Rr0 / Rr1)^γ_r
        end
    end
    return z_r
end

geometric_spreading_piecewise(r_ps, path::PathParameters) = geometric_spreading_piecewise(r_ps, path.geometric)
geometric_spreading_piecewise(r_ps, fas::FourierParameters) = geometric_spreading_piecewise(r_ps, fas.path)

"""
	geometric_spreading_cy14(r_ps::V, geo::GeometricSpreadingParameters{S,T,U}) where {S<:Float64, T<:Real, U<:AbstractVector{Bool}, V<:Real}

Geometric spreading function from Chiou & Youngs (2014).
Defines a smooth transition from one rate `γi[1]` to another `γi[2]`, with a spreading bandwidth of `Rrefi[2]` km.
"""
function geometric_spreading_cy14(r_ps::V, geo::GeometricSpreadingParameters{S,T,U}) where {S<:Float64,T<:Real,U<:AbstractVector{Bool},V<:Real}
    unit = oneunit(T)
    j = 1
    k = 1
    if geo.γfree[1] == 0
        γ1 = geo.γconi[j] * unit
        j += 1
    else
        γ1 = geo.γvari[k]
        k += 1
    end
    if geo.γfree[2] == 0
        γ2 = geo.γconi[j] * unit
    else
        γ2 = geo.γvari[k]
    end
    R0sq = (geo.Rrefi[1])^2
    Rrsq = (geo.Rrefi[2])^2
    ln_z_r = -γ1 * log(r_ps) + (-γ2 + γ1) / 2 * log((r_ps^2 + Rrsq) / (R0sq + Rrsq))
    z_r = exp(ln_z_r)
    return z_r
end

geometric_spreading_cy14(r_ps, path::PathParameters) = geometric_spreading_cy14(r_ps, path.geometric)
geometric_spreading_cy14(r_ps, fas::FourierParameters) = geometric_spreading_cy14(r_ps, fas.path)

"""
	geometric_spreading_cy14mod(r_ps::V, m::S, geo::GeometricSpreadingParameters{S,T,U}, sat::NearSourceSaturationParameters) where {S<:Float64, T<:Real, U<:AbstractVector{Bool}, V<:Real, W<:Real}

Geometric spreading function from Chiou & Youngs (2014).
Modified to make use of both `r_ps` and `r_rup` so that only the first saturation term contaminates the source amplitudes.
Defines a smooth transition from one rate `γi[1]` to another `γi[2]`, with a spreading bandwidth of `Rrefi[2]` km.
"""
function geometric_spreading_cy14mod(r_ps::V, m::W, geo::GeometricSpreadingParameters{S,T,U}, sat::NearSourceSaturationParameters) where {S<:Float64,T<:Real,U<:AbstractVector{Bool},V<:Real,W<:Real}
    unit = oneunit(T)
    j = 1
    k = 1
    if geo.γfree[1] == 0
        γ1 = geo.γconi[j] * unit
        j += 1
    else
        γ1 = geo.γvari[k]
        k += 1
    end
    if geo.γfree[2] == 0
        γ2 = geo.γconi[j] * unit
    else
        γ2 = geo.γvari[k]
    end
    R0sq = (geo.Rrefi[1])^2
    Rrsq = (geo.Rrefi[2])^2
    r_rup = rupture_distance_from_equivalent_point_source_distance(r_ps, m, sat)
    ln_z_r = -γ1 * log(r_ps) + (-γ2 + γ1) / 2 * log((r_rup^2 + Rrsq) / (R0sq + Rrsq))
    z_r = exp(ln_z_r)
    return z_r
end

geometric_spreading_cy14mod(r_ps, m, path::PathParameters) = geometric_spreading_cy14mod(r_ps, m, path.geometric, path.saturation)
geometric_spreading_cy14mod(r_ps, m, fas::FourierParameters) = geometric_spreading_cy14mod(r_ps, m, fas.path)

"""
	near_source_saturation(m, sat::NearSourceSaturationParameters)

Near-source saturation term. Used to create equivalent point-source distance.
Switches methods based upon `sat.model`.

# Arguments
Options for `sat.model` are:
- `:BT15` for Boore & Thompson (2015) finite fault factor
- `:YA15` for Yenier & Atkinson (2014) finite fault factor
- `:CY14` for a model fitted to the Chiou & Youngs (2014) saturation lengths (over all periods)
- `:None` for zero saturation length
- `:ConstantConstrained` for a constant value, `sat.hconi[1]`, not subject to AD operations
- `:ConstantVariable` for a constant value, `sat.hvari[1]`, that is subject to AD operations

Any other symbol passed will return `NaN`.

See also: [`near_source_saturation`](@ref)
"""
function near_source_saturation(m, sat::NearSourceSaturationParameters)
    unit = oneunit(get_parametric_type(sat))
    if sat.model == :BT15
        # use the Boore & Thompson (2015) finite fault factor
        return finite_fault_factor_bt15(m) * unit
    elseif sat.model == :YA15
        # use the Yenier & Atkinson (2015) finite fault factor
        return finite_fault_factor_ya15(m) * unit
    elseif sat.model == :CY14
        # use the fitted model averaging over Chiou & Youngs (2014)
        return finite_fault_factor_cy14(m) * unit
    elseif sat.model == :None
        return zero(get_parametric_type(sat))
    elseif sat.model == :ConstantConstrained
        return sat.hconi[1] * unit
    elseif sat.model == :ConstantVariable
        return sat.hvari[1] * unit
    else
        return NaN * unit
    end
end

"""
	near_source_saturation(m, path::PathParameters)

Near-source saturation term taking a `PathParameters` struct.

See also: [`near_source_saturation`](@ref)
"""
near_source_saturation(m, path::PathParameters) = near_source_saturation(m, path.saturation)

"""
	near_source_saturation(m, fas::FourierParameters)

Near-source saturation term taking a `FourierParameters` struct.

See also: [`near_source_saturation`](@ref)
"""
near_source_saturation(m, fas::FourierParameters) = near_source_saturation(m, fas.path)

"""
	finite_fault_factor_bt15(m::T) where T<:Real

Finite fault factor from Boore & Thompson (2015) used to create equivalent point-source distance.

See also: [`near_source_saturation`](@ref), [`finite_fault_factor_ya15`](@ref), [`finite_fault_factor_cy14`](@ref)
"""
function finite_fault_factor_bt15(m::T) where {T<:Real}
    Mt1 = 5.744
    Mt2 = 7.744
    if m <= Mt1
        a1 = 0.7497
        b1 = 0.4300
        h = a1 + b1 * (m - Mt1)
    elseif m >= Mt2
        a2 = 1.4147
        b2 = 0.2350
        h = a2 + b2 * (m - Mt2)
    else
        c0 = 0.7497
        c1 = 0.4300
        c2 = -0.04875
        c3 = 0.0
        h = c0 + c1 * (m - Mt1) + c2 * (m - Mt1)^2 + c3 * (m - Mt1)^3
    end
    return 10.0^h
end

"""
	finite_fault_factor_ya15(m::T) where T<:Real

Finite fault factor from Yenier & Atkinson (2015) used to create equivalent point-source distance.

See also: [`near_source_saturation`](@ref), [`finite_fault_factor_bt15`](@ref), [`finite_fault_factor_cy14`](@ref)
"""
function finite_fault_factor_ya15(m::T) where {T<:Real}
    return 10.0^(-0.405 + 0.235 * m)
end

"""
	finite_fault_factor_cy14(m::T) where T<:Real

Finite fault factor for Chiou & Youngs (2014) used to create equivalent point-source distance.
This is a model developed to match the average saturation length over the full period range.

See also: [`near_source_saturation`](@ref), [`finite_fault_factor_bt15`](@ref), [`finite_fault_factor_ya15`](@ref)
"""
function finite_fault_factor_cy14(m::T) where {T<:Real}
    hα = 7.308
    hβ = 0.4792
    hγ = 3.556
    return hα * cosh(hβ * max(m - hγ, 0.0))
end

"""
	equivalent_point_source_distance(r, m, sat::NearSourceSaturationParameters)

Compute equivalent point source distance
- `r` is ``r_{hyp}`` or ``r_{rup}`` (depending upon the size of the event -- but is nominally ``r_{rup}``)
- `m` is magnitude
- `sat` contains the `NearSourceSaturationParameters`

See also: [`near_source_saturation`](@ref)
"""
function equivalent_point_source_distance(r, m, sat::NearSourceSaturationParameters)
    h = near_source_saturation(m, sat)
    n = sat.exponent
    return (r^n + h^n)^(1 / n)
end

equivalent_point_source_distance(r, m, path::PathParameters) = equivalent_point_source_distance(r, m, path.saturation)
equivalent_point_source_distance(r, m, fas::FourierParameters) = equivalent_point_source_distance(r, m, fas.path)

"""
	rupture_distance_from_equivalent_point_source_distance(r_ps, m, sat::NearSourceSaturationParameters)

Compute rupture distance from equivalent point source distance
- `r_ps` is the equivalent point source distance
- `m` is magnitude
- `sat` contains the `NearSourceSaturationParameters`
"""
function rupture_distance_from_equivalent_point_source_distance(r_ps, m, sat::NearSourceSaturationParameters)
    h = near_source_saturation(m, sat)
    n = sat.exponent
    return (r_ps^n - h^n)^(1 / n)
end

rupture_distance_from_equivalent_point_source_distance(r_ps, m, path::PathParameters) = rupture_distance_from_equivalent_point_source_distance(r_ps, m, path.saturation)
rupture_distance_from_equivalent_point_source_distance(r_ps, m, fas::FourierParameters) = rupture_distance_from_equivalent_point_source_distance(r_ps, m, fas.path)

"""
	anelastic_attenuation(f::S, r::T, anelastic::AnelasticAttenuationParameters) where {S<:Float64,T<:Real}

Anelastic attenuation filter, computed using equivalent point source distance metric or a standard rupture distance.
"""
function anelastic_attenuation(f::S, r::T, anelastic::AnelasticAttenuationParameters) where {S<:Float64,T<:Real}
    PT = get_parametric_type(anelastic)
    Qfilt = oneunit(PT)
    jq = jη = 1
    kq = kη = 1
    for i = 1:length(anelastic.Qfree)
        @inbounds Rr0 = anelastic.Rrefi[i]
        @inbounds Rr1 = anelastic.Rrefi[i+1]

        Q0_r = oneunit(PT)
        η_r = oneunit(PT)
        cQ_r = anelastic.cQ[i]

        # get the relevant constrained or free quality factor for this path segment
        if anelastic.Qfree[i] == 0
            @inbounds Q0_r *= anelastic.Q0coni[jq]
            jq += 1
        else
            @inbounds Q0_r *= anelastic.Q0vari[kq]
            kq += 1
        end
        # get the relevant constrained of free quality exponent for this path segment
        if anelastic.ηfree[i] == 0
            @inbounds η_r *= anelastic.ηconi[jη]
            jη += 1
        else
            @inbounds η_r *= anelastic.ηvari[kη]
            kη += 1
        end

        if r < Rr1
            Qfilt *= exp(-π * f^(1.0 - η_r) * (r - Rr0) / (Q0_r * cQ_r))
            return Qfilt
        else
            Qfilt *= exp(-π * f^(1.0 - η_r) * (Rr1 - Rr0) / (Q0_r * cQ_r))
        end
    end
end

anelastic_attenuation(f, r, path::PathParameters) = anelastic_attenuation(f, r, path.anelastic)
anelastic_attenuation(f, r, fas::FourierParameters) = anelastic_attenuation(f, r, fas.path)

"""
	effective_quality_parameters(r::T, anelastic::AnelasticAttenuationParameters) where {T<:Real}

Distance weighted quality parameters for segmented model.
Returns a tuple of effective Q0, η & cQ values from a weighted sum, weighted by relative distance in each path segment.
"""
function effective_quality_parameters(r::T, anelastic::AnelasticAttenuationParameters) where {T<:Real}
    PT = get_parametric_type(anelastic)
    Q_eff = zero(PT)
    η_eff = zero(PT)
    cQ_eff = zero(PT)
    jq = jη = 1
    kq = kη = 1
    for i = 1:length(anelastic.Qfree)
        @inbounds Rr0 = anelastic.Rrefi[i]
        @inbounds Rr1 = anelastic.Rrefi[i+1]

        Q0_r = oneunit(PT)
        η_r = oneunit(PT)
        cQ_r = anelastic.cQ[i]

        # get the relevant constrained or free quality factor for this path segment
        if anelastic.Qfree[i] == 0
            @inbounds Q0_r *= anelastic.Q0coni[jq]
            jq += 1
        else
            @inbounds Q0_r *= anelastic.Q0vari[kq]
            kq += 1
        end
        # get the relevant constrained of free quality exponent for this path segment
        if anelastic.ηfree[i] == 0
            @inbounds η_r *= anelastic.ηconi[jη]
            jη += 1
        else
            @inbounds η_r *= anelastic.ηvari[kη]
            kη += 1
        end

        if r < Rr1
            rfrac = (r - Rr0) / r
            Q_eff += rfrac * Q0_r
            η_eff += rfrac * η_r
            cQ_eff += rfrac * cQ_r
            return (Q_eff, η_eff, cQ_eff)
        else
            rfrac = (Rr1 - Rr0) / r
            Q_eff += rfrac * Q0_r
            η_eff += rfrac * η_r
            cQ_eff += rfrac * cQ_r
        end
    end
end


text:
include("../batch/figure-whitenoise.jl")

include("../batch/figure-mle-accuracy.jl")

include("../batch/figure-cmp-estimators.jl")

include("../batch/figure-fit-kld.jl")

choice = copy_icmla_figures(figure_dir)

include("figure-rainfall.jl")

copy_icmla_figures(results_dir,choice)
nothing


text:
module CorpusLoaders
using DataDeps
using WordTokenizers
using MultiResolutionIterators
using InternedStrings
using Glob
using StringEncodings

export Document, TaggedWord, SenseAnnotatedWord, PosTaggedWord
export title, sensekey, word
export load

export WikiCorpus, SemCor, Senseval3

function __init__()

    include("WikiCorpus_DataDeps.jl")
    include("SemCor_DataDeps.jl")
    include("SemEval2007Task7_DataDeps.jl")
    include("Senseval3_DataDeps.jl")
end

include("types.jl")
include("apply_subparsers.jl")
include("WikiCorpus.jl")
include("SemCor.jl")
include("SemEval2007Task7.jl")
include("Senseval3.jl")

end


text: Dict{String, Any}("h1" => Dict{String, Any}("value" => "0", "type" => "integer"), "a2" => Dict{String, Any}("value" => "0", "type" => "integer"), "b2" => Dict{String, Any}("value" => "0", "type" => "integer"), "h3" => Dict{String, Any}("value" => "0", "type" => "integer"), "o1" => Dict{String, Any}("value" => "0", "type" => "integer"), "d2" => Dict{String, Any}("value" => "0", "type" => "integer"), "a3" => Dict{String, Any}("value" => "0", "type" => "integer"), "d1" => Dict{String, Any}("value" => "0", "type" => "integer"), "h2" => Dict{String, Any}("value" => "0", "type" => "integer"), "b3" => Dict{String, Any}("value" => "0", "type" => "integer"), "d3" => Dict{String, Any}("value" => "0", "type" => "integer"), "b1" => Dict{String, Any}("value" => "0", "type" => "integer"))


text:
# STO-nG basis definition

abstract type BasisFunction end
abstract type GTO <: BasisFunction end
abstract type GaussianFunction <: GTO end

"Gaussian s wave function"
struct Gaussian_s <: GaussianFunction
    α::Float64       # gussian function exponent
    R::Vec3{Float64} # center (or actually nuclues' position)
end

"Gaussian p wave function"
struct Gaussian_p <:GaussianFunction
    α::Float64       # gussian function exponent
    R::Vec3{Float64} # center
    n::Vec3{Float64} # orbital direction
end

"STO-nG function"
struct STONG{N} <: GTO
    d::NTuple{N, Float64}           # coefficients
    g::NTuple{N, GaussianFunction}  # gaussian functions
    STONG{N}(d::NTuple{N, Float64}, g::NTuple{N, GaussianFunction}) where N = new(d, g)
end

STONG(d::NTuple{N, Float64}, g::NTuple{N, GaussianFunction}) where N = STONG{N}(d, g)

struct STONG_coeff{N}
    coefficients::NTuple{N, Float64}
    exponents::NTuple{N,Float64}
    STONG_coeff{N}(coeff::NTuple{N, Float64}, _exp::NTuple{N, Float64}) where N = new(coeff, _exp)
end

"construct STONG from STONG_coeff"
STONG(coeff::STONG_coeff{N}, center::Vec3{Float64}; type::Symbol=:s) where N = begin
    if type === :s
        return STONG{N}(
            coeff.coefficients,
            ntuple(i->Gaussian_s(coeff.exponents[i], center), N)
        )
    elseif type === :px
        return STONG{N}(
            coeff.coefficients,
            ntuple(i->Gaussian_p(coeff.exponents[i], center, vec3(1., 0., 0.)), N)
        )
    elseif type === :py
        return STONG{N}(
            coeff.coefficients,
            ntuple(i->Gaussian_p(coeff.exponents[i], center, vec3(0., 1., 0.)), N)
        )
    elseif type === :pz
        return STONG{N}(
            coeff.coefficients,
            ntuple(i->Gaussian_p(coeff.exponents[i], center, vec3(0., 0., 1.)), N)
        )
    end
    error("orbital type $type is not support yet!")
end

const STO3G = STONG{3}
const STO3G_coeff = STONG_coeff{3}

# the data comes from basis_set_exchange
# https://github.com/MolSSI-BSE/basis_set_exchange
const coefficient_data = Dict(
    "H" => Dict(
        "1s" => STO3G_coeff(
            (0.1543289673E+00, 0.5353281423E+00, 0.4446345422E+00), # coefficients
            (0.3425250914E+01, 0.6239137298E+00, 0.1688554040E+00)  # exponents
        )
    ),
    "He" => Dict(
        "1s" => STO3G_coeff(
            (0.1543289673E+00, 0.5353281423E+00, 0.4446345422E+00), # coefficients
            (0.6362421394E+01, 0.1158922999E+01, 0.3136497915E+00)  # exponents
        )
    ),
    "C" => Dict(
        "1s" => STO3G_coeff(
            (0.1543289673E+00, 0.5353281423E+00, 0.4446345422E+00), # coefficients
            (0.7161683735E+02, 0.1304509632E+02, 0.3530512160E+01)  # exponents
        ),
        "2s" => STO3G_coeff(
            (-0.9996722919E-01, 0.3995128261E+00, 0.7001154689E+00),# coefficients
            (0.2941249355E+01, 0.6834830964E+00, 0.2222899159E+00)  # exponents
        ),
        "2p" => STO3G_coeff(
            (0.1559162750E+00, 0.6076837186E+00, 0.3919573931E+00), # coefficients
            (0.2941249355E+01, 0.6834830964E+00, 0.2222899159E+00)  # exponents
        )
    ),
    "N" => Dict(
        "1s" => STO3G_coeff(
            (0.1543289673E+00, 0.5353281423E+00, 0.4446345422E+00), # coefficients
            (0.9910616896E+02, 0.1805231239E+02, 0.4885660238E+01)  # exponents
        ),
        "2s" => STO3G_coeff(
            (-0.9996722919E-01, 0.3995128261E+00, 0.7001154689E+00),# coefficients
            (0.3780455879E+01, 0.8784966449E+00, 0.2857143744E+00)  # exponents
        ),
        "2p" => STO3G_coeff(
            (0.1559162750E+00, 0.6076837186E+00, 0.3919573931E+00), # coefficients
            (0.3780455879E+01, 0.8784966449E+00, 0.2857143744E+00)  # exponents
        )
    ),
    "O" => Dict(
        "1s" => STO3G_coeff(
            (0.1543289673E+00, 0.5353281423E+00, 0.4446345422E+00), # coefficients
            (0.1307093214E+03, 0.2380886605E+02, 0.6443608313E+01)  # exponents
        ),
        "2s" => STO3G_coeff(
            (-0.9996722919E-01, 0.3995128261E+00, 0.7001154689E+00),# coefficients
            (0.5033151319E+01, 0.1169596125E+01, 0.3803889600E+00)  # exponents
        ),
        "2p" => STO3G_coeff(
            (0.1559162750E+00, 0.6076837186E+00, 0.3919573931E+00), # coefficients
            (0.5033151319E+01, 0.1169596125E+01, 0.3803889600E+00)  # exponents
        )
    ),
    "F" => Dict(
        "1s" => STO3G_coeff(
            (0.1543289673E+00, 0.5353281423E+00, 0.4446345422E+00), # coefficients
            (0.1666791340E+03, 0.3036081233E+02, 0.8216820672E+01)  # exponents
        ),
        "2s" => STO3G_coeff(
            (-0.9996722919E-01, 0.3995128261E+00, 0.7001154689E+00),# coefficients
            (0.6464803249E+01, 0.1502281245E+01, 0.4885884864E+00)  # exponents
        ),
        "2p" => STO3G_coeff(
            (0.1559162750E+00, 0.6076837186E+00, 0.3919573931E+00), # coefficients
            (0.6464803249E+01, 0.1502281245E+01, 0.4885884864E+00)  # exponents
        )
    )
)

Base.length(::STONG{N}) where N = N

Base.show(io::IO, gs::Gaussian_s) = begin
    print(io, "gs(α=$(gs.α), R=$(gs.R))")
end

Base.show(io::IO, gp::Gaussian_p) = begin
    print(io, "gp(α=$(gp.α), R=$(gp.R), n=$(gp.n))")
end

Base.show(io::IO, stong::STONG{N}) where N = begin
    println(io, "STO-$(N)G")
    for i = 1:N
        println(io, "  $(stong.d[i]) × $(stong.g[i])")
    end
end


text:
export least_squares

function least_squares(Dpred,Dtrue)
    r = Dpred-Dtrue
    f = 0.5*vecnorm(r)^2
    g = r
    h = ones(eltype(Dpred),size(Dpred))
    return (f,g,h)
end


text:
# Julia wrapper for header: /usr/include/scip/scip_prob.h
# Automatically generated using Clang.jl wrap_c

function SCIPcreateProb(scip, name, probdelorig, probtrans, probdeltrans, probinitsol, probexitsol, probcopy, probdata)
    ccall((:SCIPcreateProb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_PROBDATA}), scip, name, probdelorig, probtrans, probdeltrans, probinitsol, probexitsol, probcopy, probdata)
end

function SCIPcreateProbBasic(scip, name)
    ccall((:SCIPcreateProbBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPsetProbDelorig(scip, probdelorig)
    ccall((:SCIPsetProbDelorig, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}), scip, probdelorig)
end

function SCIPsetProbTrans(scip, probtrans)
    ccall((:SCIPsetProbTrans, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}), scip, probtrans)
end

function SCIPsetProbDeltrans(scip, probdeltrans)
    ccall((:SCIPsetProbDeltrans, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}), scip, probdeltrans)
end

function SCIPsetProbInitsol(scip, probinitsol)
    ccall((:SCIPsetProbInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}), scip, probinitsol)
end

function SCIPsetProbExitsol(scip, probexitsol)
    ccall((:SCIPsetProbExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}), scip, probexitsol)
end

function SCIPsetProbCopy(scip, probcopy)
    ccall((:SCIPsetProbCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}), scip, probcopy)
end

function SCIPreadProb(scip, filename, extension)
    ccall((:SCIPreadProb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring), scip, filename, extension)
end

function SCIPwriteOrigProblem(scip, filename, extension, genericnames)
    ccall((:SCIPwriteOrigProblem, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, UInt32), scip, filename, extension, genericnames)
end

function SCIPwriteTransProblem(scip, filename, extension, genericnames)
    ccall((:SCIPwriteTransProblem, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, UInt32), scip, filename, extension, genericnames)
end

function SCIPfreeProb(scip)
    ccall((:SCIPfreeProb, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPpermuteProb(scip, randseed, permuteconss, permutebinvars, permuteintvars, permuteimplvars, permutecontvars)
    ccall((:SCIPpermuteProb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, randseed, permuteconss, permutebinvars, permuteintvars, permuteimplvars, permutecontvars)
end

function SCIPgetProbData(scip)
    ccall((:SCIPgetProbData, libscip), Ptr{SCIP_PROBDATA}, (Ptr{SCIP_},), scip)
end

function SCIPsetProbData(scip, probdata)
    ccall((:SCIPsetProbData, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PROBDATA}), scip, probdata)
end

function SCIPgetProbName(scip)
    ccall((:SCIPgetProbName, libscip), Cstring, (Ptr{SCIP_},), scip)
end

function SCIPsetProbName(scip, name)
    ccall((:SCIPsetProbName, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPchgReoptObjective(scip, objsense, vars, coefs, nvars)
    ccall((:SCIPchgReoptObjective, libscip), SCIP_RETCODE, (Ptr{SCIP_}, SCIP_OBJSENSE, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint), scip, objsense, vars, coefs, nvars)
end

function SCIPgetObjsense(scip)
    ccall((:SCIPgetObjsense, libscip), SCIP_OBJSENSE, (Ptr{SCIP_},), scip)
end

function SCIPsetObjsense(scip, objsense)
    ccall((:SCIPsetObjsense, libscip), SCIP_RETCODE, (Ptr{SCIP_}, SCIP_OBJSENSE), scip, objsense)
end

function SCIPaddObjoffset(scip, addval)
    ccall((:SCIPaddObjoffset, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble), scip, addval)
end

function SCIPaddOrigObjoffset(scip, addval)
    ccall((:SCIPaddOrigObjoffset, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble), scip, addval)
end

function SCIPgetOrigObjoffset(scip)
    ccall((:SCIPgetOrigObjoffset, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetOrigObjscale(scip)
    ccall((:SCIPgetOrigObjscale, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetTransObjoffset(scip)
    ccall((:SCIPgetTransObjoffset, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetTransObjscale(scip)
    ccall((:SCIPgetTransObjscale, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPsetObjlimit(scip, objlimit)
    ccall((:SCIPsetObjlimit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble), scip, objlimit)
end

function SCIPgetObjlimit(scip)
    ccall((:SCIPgetObjlimit, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPsetObjIntegral(scip)
    ccall((:SCIPsetObjIntegral, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPisObjIntegral(scip)
    ccall((:SCIPisObjIntegral, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPgetObjNorm(scip)
    ccall((:SCIPgetObjNorm, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPaddVar(scip, var)
    ccall((:SCIPaddVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPaddPricedVar(scip, var, score)
    ccall((:SCIPaddPricedVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, score)
end

function SCIPdelVar(scip, var, deleted)
    ccall((:SCIPdelVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{UInt32}), scip, var, deleted)
end

function SCIPgetVarsData(scip, vars, nvars, nbinvars, nintvars, nimplvars, ncontvars)
    ccall((:SCIPgetVarsData, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), scip, vars, nvars, nbinvars, nintvars, nimplvars, ncontvars)
end

function SCIPgetVars(scip)
    ccall((:SCIPgetVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNVars(scip)
    ccall((:SCIPgetNVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNBinVars(scip)
    ccall((:SCIPgetNBinVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNIntVars(scip)
    ccall((:SCIPgetNIntVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNImplVars(scip)
    ccall((:SCIPgetNImplVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNContVars(scip)
    ccall((:SCIPgetNContVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNObjVars(scip)
    ccall((:SCIPgetNObjVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetFixedVars(scip)
    ccall((:SCIPgetFixedVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNFixedVars(scip)
    ccall((:SCIPgetNFixedVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetOrigVarsData(scip, vars, nvars, nbinvars, nintvars, nimplvars, ncontvars)
    ccall((:SCIPgetOrigVarsData, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), scip, vars, nvars, nbinvars, nintvars, nimplvars, ncontvars)
end

function SCIPgetOrigVars(scip)
    ccall((:SCIPgetOrigVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNOrigVars(scip)
    ccall((:SCIPgetNOrigVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNOrigBinVars(scip)
    ccall((:SCIPgetNOrigBinVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNOrigIntVars(scip)
    ccall((:SCIPgetNOrigIntVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNOrigImplVars(scip)
    ccall((:SCIPgetNOrigImplVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNOrigContVars(scip)
    ccall((:SCIPgetNOrigContVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNTotalVars(scip)
    ccall((:SCIPgetNTotalVars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetSolVarsData(scip, sol, vars, nvars, nbinvars, nintvars, nimplvars, ncontvars)
    ccall((:SCIPgetSolVarsData, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SOL}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), scip, sol, vars, nvars, nbinvars, nintvars, nimplvars, ncontvars)
end

function SCIPfindVar(scip, name)
    ccall((:SCIPfindVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPallVarsInProb(scip)
    ccall((:SCIPallVarsInProb, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPaddCons(scip, cons)
    ccall((:SCIPaddCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPdelCons(scip, cons)
    ccall((:SCIPdelCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPfindOrigCons(scip, name)
    ccall((:SCIPfindOrigCons, libscip), Ptr{SCIP_CONS}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPfindCons(scip, name)
    ccall((:SCIPfindCons, libscip), Ptr{SCIP_CONS}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetNUpgrConss(scip)
    ccall((:SCIPgetNUpgrConss, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNConss(scip)
    ccall((:SCIPgetNConss, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetConss(scip)
    ccall((:SCIPgetConss, libscip), Ptr{Ptr{SCIP_CONS}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNOrigConss(scip)
    ccall((:SCIPgetNOrigConss, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetOrigConss(scip)
    ccall((:SCIPgetOrigConss, libscip), Ptr{Ptr{SCIP_CONS}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNCheckConss(scip)
    ccall((:SCIPgetNCheckConss, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPaddConflict(scip, node, cons, validnode, conftype, iscutoffinvolved)
    ccall((:SCIPaddConflict, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Ptr{SCIP_CONS}, Ptr{SCIP_NODE}, SCIP_CONFTYPE, UInt32), scip, node, cons, validnode, conftype, iscutoffinvolved)
end

function SCIPclearConflictStore(scip, event)
    ccall((:SCIPclearConflictStore, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_EVENT}), scip, event)
end

function SCIPaddConsNode(scip, node, cons, validnode)
    ccall((:SCIPaddConsNode, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Ptr{SCIP_CONS}, Ptr{SCIP_NODE}), scip, node, cons, validnode)
end

function SCIPaddConsLocal(scip, cons, validnode)
    ccall((:SCIPaddConsLocal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_NODE}), scip, cons, validnode)
end

function SCIPdelConsNode(scip, node, cons)
    ccall((:SCIPdelConsNode, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Ptr{SCIP_CONS}), scip, node, cons)
end

function SCIPdelConsLocal(scip, cons)
    ccall((:SCIPdelConsLocal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLocalOrigEstimate(scip)
    ccall((:SCIPgetLocalOrigEstimate, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLocalTransEstimate(scip)
    ccall((:SCIPgetLocalTransEstimate, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLocalDualbound(scip)
    ccall((:SCIPgetLocalDualbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLocalLowerbound(scip)
    ccall((:SCIPgetLocalLowerbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetNodeDualbound(scip, node)
    ccall((:SCIPgetNodeDualbound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_NODE}), scip, node)
end

function SCIPgetNodeLowerbound(scip, node)
    ccall((:SCIPgetNodeLowerbound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_NODE}), scip, node)
end

function SCIPupdateLocalDualbound(scip, newbound)
    ccall((:SCIPupdateLocalDualbound, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble), scip, newbound)
end

function SCIPupdateLocalLowerbound(scip, newbound)
    ccall((:SCIPupdateLocalLowerbound, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble), scip, newbound)
end

function SCIPupdateNodeDualbound(scip, node, newbound)
    ccall((:SCIPupdateNodeDualbound, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Cdouble), scip, node, newbound)
end

function SCIPupdateNodeLowerbound(scip, node, newbound)
    ccall((:SCIPupdateNodeLowerbound, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Cdouble), scip, node, newbound)
end

function SCIPchgChildPrio(scip, child, priority)
    ccall((:SCIPchgChildPrio, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Cdouble), scip, child, priority)
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, Threadicare

makedocs(
     modules = [Threadicare],
       clean = true,
     doctest = true,
   checkdocs = :all,
      format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
     authors = "Andre N. Souza",
    sitename = "Threadicare.jl",

       pages = Any[
                "Home" => "index.md"
                      ]
)

deploydocs(repo = "github.com/sandreza/Threadicare.jl.git")

max_stars_repo_name: sandreza/Threadicare


text:
" Allometry. Scaling rules to relate size to mass. "
abstract type AbstractAllometry end

@columns @flattenable struct SqrtAllometry{B0,B1} <: AbstractAllometry
    # Field       | Flatn | Default | Unit | Pror            | Limits         | Log  | Description
    β0::B0        | false | 1e-4*24 | g    | Gamma(2.0, 0.2) | [1e-6, 10.00]  | true | "Intercept"
    β1::B1        | true  | 0.1     | m    | Gamma(2.0, 0.2) | [1e-2, 1.0]    | true | "Scalar"
end

@columns @flattenable struct Allometry{B0,B1,A} <: AbstractAllometry
    β0::B0        | false | 1e-4*24 | g    | Gamma(2.0, 0.2) | [1e-6, 10.00]  | true | "Intercept"
    β1::B1        | true  | 0.1     | m    | Gamma(2.0, 0.2) | [1e-5, 10.00]  | true | "Scalar"
    α::A          | true  | 0.1     | _    | Gamma(2.0, 0.2) | [1e-3, 100.00] | true | "Exponent"
end

@columns struct FixedAllometry{M} <: AbstractAllometry
    height::M             | 1.0     | m    | Gamma(2.0, 0.2) | [1e-3, 100.00] | true | "Fixed height or depth"
end

update_height!(o, u) = update_height!(allometry_pars(o), o, u)
update_height!(f::Nothing, o, u) = nothing
update_height!(f, o, u) = set_height!(o, allometric_height(f, mass(o, u)))

allometric_height(f::FixedAllometry, mass) = f.height
allometric_height(f::SqrtAllometry, mass) = sqrt((mass - f.β0) / unit(f.β0)) * f.β1
allometric_height(f::Allometry, mass) = begin
    f.β1 * ((max(zero(mass), mass - f.β0))/unit(f.β0))^f.α
end


text:
using Test
using CcallMacros: @ccall, parsecall, CcallError

@testset "test basic parsecall functionality" begin
    callexpr = :(
        libc.printf("%d"::Cstring, value::Cuint)::Cvoid
    )
    @test parsecall(callexpr) == (
        :(:printf, libc),     # function
        :Cvoid,               # return type
        :((Cstring, Cuint)),  # argument types
        ["%d", :value]        # argument symbols
    )
end

@testset "ensure the base-case of @ccall works, including library name" begin
    call = @macroexpand @ccall libstring.func(
        str::Cstring,
        num1::Cint,
        num2::Cint
    )::Cstring

    @test call == :(
        ccall((:func, libstring), Cstring,
            (Cstring, Cint, Cint), str, num1, num2)
    )
end

@testset "ensure @ccall handles varargs correctly" begin
    call = @macroexpand @ccall printf(
        "%d, %d, %d\n"::Cstring ; 1::Cint, 2::Cint, 3::Cint
    )::Cint
    @test call == :(
        ccall(:printf, Cint, (Cstring, Cint...), "%d, %d, %d\n", 1, 2, 3)
    )
end

@testset "ensure parsecall throws errors appropriately" begin
    # missing return type
    @test_throws CcallError parsecall(:( foo(4.0::Cdouble )))
    # not a function call
    @test_throws CcallError parsecall(:( foo::Type ))
    # mismatched types on varargs
    @test_throws CcallError parsecall(:( foo(x::Cint; y::Cstring, z::Cint)::Cvoid ))
    # missing type annotations on arguments
    @test_throws CcallError parsecall(:( foo(x)::Cint ))
end

# call some c functions
@testset "run @ccall with C standard library functions" begin
    @test @ccall(sqrt(4.0::Cdouble)::Cdouble) == 2.0

    STRING = "hello"
    BUFFER = Ptr{UInt8}(Libc.malloc((length(STRING) + 1) * sizeof(Cchar)))
    @ccall strcpy(BUFFER::Cstring, STRING::Cstring)::Cstring
    @test unsafe_string(BUFFER) == STRING

    # let's write C in Julia. Uppercasing the hard way.
    let buffer = BUFFER
        while (byte = unsafe_load(buffer)) != 0x00
            bigger = @ccall toupper(byte::Cint)::Cint
            unsafe_store!(buffer, bigger)
            buffer += 1
        end
    end

    @test unsafe_string(BUFFER) == uppercase(STRING)
    Libc.free(BUFFER)
end


text:
function test_permutation_permutations()
    @test collect(permutations(Int[], Val{:unique})) == Array{Int,1}[[]]
    @test length(permutations(Int[], Val{:unique})) == 1
    @test eltype(permutations(Int[], Val{:unique})) == Array{Int,1}

    @test collect(permutations(["a", "b", "c"], Val{:unique})) == Array[
        ["a", "b", "c"],
        ["a", "c", "b"],
        ["b", "a", "c"],
        ["b", "c", "a"],
        ["c", "a", "b"],
        ["c", "b", "a"],
    ]
    @test length(permutations(["a", "b", "c"], Val{:unique})) == 6
    @test eltype(permutations(["a", "b", "c"], Val{:unique})) == Array{String,1}
end

function test_permutation_permutations_with_repetition()
    @test collect(permutations(Int[], Val{:repeated})) == Array{Int,1}[[]]
    @test length(permutations(Int[], Val{:repeated})) == 1
    @test eltype(permutations(Int[], Val{:repeated})) == Array{Int,1}

    @test collect(permutations(["a", "b", "c"], Val{:repeated})) == Array[
        ["a", "a", "a"],
        ["a", "a", "b"],
        ["a", "a", "c"],
        ["a", "b", "a"],
        ["a", "b", "b"],
        ["a", "b", "c"],
        ["a", "c", "a"],
        ["a", "c", "b"],
        ["a", "c", "c"],
        ["b", "a", "a"],
        ["b", "a", "b"],
        ["b", "a", "c"],
        ["b", "b", "a"],
        ["b", "b", "b"],
        ["b", "b", "c"],
        ["b", "c", "a"],
        ["b", "c", "b"],
        ["b", "c", "c"],
        ["c", "a", "a"],
        ["c", "a", "b"],
        ["c", "a", "c"],
        ["c", "b", "a"],
        ["c", "b", "b"],
        ["c", "b", "c"],
        ["c", "c", "a"],
        ["c", "c", "b"],
        ["c", "c", "c"],
    ]
    @test length(permutations(["a", "b", "c"], Val{:repeated})) == 27
    @test eltype(permutations(["a", "b", "c"], Val{:repeated})) == Array{String,1}
end

function test_permutation_permutations_unknown_mode()
    @test_throws MethodError permutations(["a", "b", "c"], Val{:unknown})
end

function test_permutation_kpermutations()
    @test collect(permutations(Int[], 0, Val{:unique})) == Array{Int,1}[[]]
    @test length(permutations(Int[], 0, Val{:unique})) == 1
    @test eltype(permutations(Int[], 0, Val{:unique})) == Array{Int,1}

    @test collect(permutations(Int[], 1, Val{:unique})) == Array{Int,1}[]
    @test length(permutations(Int[], 1, Val{:unique})) == 0
    @test eltype(permutations(Int[], 1, Val{:unique})) == Array{Int,1}

    @test collect(permutations(["a", "b", "c", "d"], 2, Val{:unique})) == Array[
        ["a", "b"],
        ["a", "c"],
        ["a", "d"],
        ["b", "a"],
        ["b", "c"],
        ["b", "d"],
        ["c", "a"],
        ["c", "b"],
        ["c", "d"],
        ["d", "a"],
        ["d", "b"],
        ["d", "c"],
    ]
    @test length(permutations(["a", "b", "c", "d"], 2, Val{:unique})) == 12
    @test eltype(permutations(["a", "b", "c", "d"], 2, Val{:unique})) == Array{String,1}

    @test collect(permutations(["a", "b", "c", "d"], 3, Val{:unique})) == Array[
        ["a", "b", "c"],
        ["a", "b", "d"],
        ["a", "c", "b"],
        ["a", "c", "d"],
        ["a", "d", "b"],
        ["a", "d", "c"],
        ["b", "a", "c"],
        ["b", "a", "d"],
        ["b", "c", "a"],
        ["b", "c", "d"],
        ["b", "d", "a"],
        ["b", "d", "c"],
        ["c", "a", "b"],
        ["c", "a", "d"],
        ["c", "b", "a"],
        ["c", "b", "d"],
        ["c", "d", "a"],
        ["c", "d", "b"],
        ["d", "a", "b"],
        ["d", "a", "c"],
        ["d", "b", "a"],
        ["d", "b", "c"],
        ["d", "c", "a"],
        ["d", "c", "b"],
    ]
    @test length(permutations(["a", "b", "c", "d"], 3, Val{:unique})) == 24
    @test eltype(permutations(["a", "b", "c", "d"], 3, Val{:unique})) == Array{String,1}

    @test collect(permutations(["a", "b"], 3, Val{:unique})) == Array{String,1}[]
    @test length(permutations(["a", "b"], 3, Val{:unique})) == 0
end

function test_permutation_kpermutations_with_repetition()
    @test collect(permutations(Int[], 0, Val{:repeated})) == Array{Int,1}[[]]
    @test length(permutations(Int[], 0, Val{:repeated})) == 1
    @test eltype(permutations(Int[], 0, Val{:repeated})) == Array{Int,1}

    @test collect(permutations(Int[], 1, Val{:repeated})) == Array{Int,1}[]
    @test length(permutations(Int[], 1, Val{:repeated})) == 0
    @test eltype(permutations(Int[], 1, Val{:repeated})) == Array{Int,1}

    @test collect(permutations(["a", "b", "c", "d"], 2, Val{:repeated})) == Array[
        ["a", "a"],
        ["a", "b"],
        ["a", "c"],
        ["a", "d"],
        ["b", "a"],
        ["b", "b"],
        ["b", "c"],
        ["b", "d"],
        ["c", "a"],
        ["c", "b"],
        ["c", "c"],
        ["c", "d"],
        ["d", "a"],
        ["d", "b"],
        ["d", "c"],
        ["d", "d"],
    ]
    @test length(permutations(["a", "b", "c", "d"], 2, Val{:repeated})) == 4^2
    @test eltype(permutations(["a", "b", "c", "d"], 2, Val{:repeated})) == Array{String,1}

    @test collect(permutations(["a", "b"], 3, Val{:repeated})) == Array[
        ["a", "a", "a"],
        ["a", "a", "b"],
        ["a", "b", "a"],
        ["a", "b", "b"],
        ["b", "a", "a"],
        ["b", "a", "b"],
        ["b", "b", "a"],
        ["b", "b", "b"],
    ]
    @test length(permutations(["a", "b"], 3, Val{:repeated})) == 2^3
end

function test_permutation_kpermutations_unknown_mode()
    @test_throws MethodError permutations(["a", "b", "c", "d"], 2, Val{:unknown})
end

function test_permutation_all()
    print(rpad("Math.Combinatorics.Permutation...", 50, ' '))

    test_permutation_permutations()
    test_permutation_permutations_with_repetition()
    test_permutation_permutations_unknown_mode()
    test_permutation_kpermutations()
    test_permutation_kpermutations_with_repetition()
    test_permutation_kpermutations_unknown_mode()

    println("PASS")
end


text:
function raw_terminal(f, terminal=terminal, hide_cursor=true)
  cur_enabled = _enable_raw_terminal(terminal)

  hide_cursor &= cur_enabled

  hide_cursor && _cursor_hide(terminal.out_stream)

  try
    f()
  finally
    cur_enabled && _disable_raw_terminal(terminal)
    hide_cursor && _cursor_show(terminal.out_stream)
  end
end

function _enable_raw_terminal(terminal)
  try
    REPL.Terminals.raw!(terminal, true)
    return true
  catch err
    warn("TerminalMenus: Unable to enable raw mode: $err")
  end

  false
end

function _disable_raw_terminal(terminal)
  try
    REPL.Terminals.raw!(terminal, false)
    return true
  catch err
    warn("TerminalMenus: Unable to disable raw mode: $err")
  end

  false
end

_cursor_hide(cur_buffer::IO=STDOUT) = print(cur_buffer, "\x1b[?25l")
_cursor_show(cur_buffer::IO=STDOUT) = print(cur_buffer, "\x1b[?25h")


text:
# ordered sets

# This was largely copied and modified from Base

# TODO: Most of these functions should be removed once AbstractSet is introduced there
# (see https://github.com/JuliaLang/julia/issues/5533)

immutable OrderedSet{T}
    dict::OrderedDict{T,Void}

    OrderedSet() = new(OrderedDict{T,Void}())
    OrderedSet(xs) = union!(new(OrderedDict{T,Void}()), xs)
end
OrderedSet() = OrderedSet{Any}()
OrderedSet(xs) = OrderedSet{eltype(xs)}(xs)

show(io::IO, s::OrderedSet) = (show(io, typeof(s)); print(io, "("); !isempty(s) && Base.show_comma_array(io, s,'[',']'); print(io, ")"))

@delegate OrderedSet.dict [isempty, length]

sizehint!(s::OrderedSet, sz::Integer) = (sizehint!(s.dict, sz); s)
eltype{T}(s::OrderedSet{T}) = T

in(x, s::OrderedSet) = haskey(s.dict, x)

push!(s::OrderedSet, x) = (s.dict[x] = nothing; s)
pop!(s::OrderedSet, x) = (pop!(s.dict, x); x)
pop!(s::OrderedSet, x, deflt) = pop!(s.dict, x, deflt) == deflt ? deflt : x
delete!(s::OrderedSet, x) = (delete!(s.dict, x); s)

getindex(x::OrderedSet,i::Int) = x.dict.keys[i]
endof(x::OrderedSet) = endof(x.dict.keys)

union!(s::OrderedSet, xs) = (for x in xs; push!(s,x); end; s)
setdiff!(s::OrderedSet, xs) = (for x in xs; delete!(s,x); end; s)
setdiff!(s::Set, xs::OrderedSet) = (for x in xs; delete!(s,x); end; s)

similar{T}(s::OrderedSet{T}) = OrderedSet{T}()
copy(s::OrderedSet) = union!(similar(s), s)

empty!{T}(s::OrderedSet{T}) = (empty!(s.dict); s)

start(s::OrderedSet)       = start(s.dict)
done(s::OrderedSet, state) = done(s.dict, state)
# NOTE: manually optimized to take advantage of OrderedDict representation
next(s::OrderedSet, i)     = (s.dict.keys[i], i+1)

pop!(s::OrderedSet) = pop!(s.dict)[1]

union(s::OrderedSet) = copy(s)
function union(s::OrderedSet, sets...)
    u = OrderedSet{Base.join_eltype(s, sets...)}()
    union!(u,s)
    for t in sets
        union!(u,t)
    end
    return u
end

intersect(s::OrderedSet) = copy(s)
function intersect(s::OrderedSet, sets...)
    i = copy(s)
    for x in s
        for t in sets
            if !in(x,t)
                delete!(i,x)
                break
            end
        end
    end
    return i
end

function setdiff(a::OrderedSet, b)
    d = similar(a)
    for x in a
        if !(x in b)
            push!(d, x)
        end
    end
    d
end

==(l::OrderedSet, r::OrderedSet) = (length(l) == length(r)) && (l <= r)
<(l::OrderedSet, r::OrderedSet) = (length(l) < length(r)) && (l <= r)
<=(l::OrderedSet, r::OrderedSet) = issubset(l, r)

function filter!(f::Function, s::OrderedSet)
    for x in s
        if !f(x)
            delete!(s, x)
        end
    end
    return s
end
filter(f::Function, s::OrderedSet) = filter!(f, copy(s))

const orderedset_seed = UInt === UInt64 ? 0x2114638a942a91a5 : 0xd86bdbf1
function hash(s::OrderedSet, h::UInt)
    h = hash(orderedset_seed, h)
    s.dict.ndel > 0 && rehash!(s.dict)
    hash(s.dict.keys, h)
end


text:
A = UInt8[200 150; 50 1]
img = grayim(A)
img2 = grayim(A+10)
B = UInt8[200 150 100; 80 50 10; 5 3 0]
img3 = grayim(B)

@testset "ImageTransformation" begin
    @test_throws ArgumentError multiplier(FaultyOp())
end

@testset "NoOp" begin
    @test NoOp <: ImageTransformation
    op = NoOp()
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 1
    @test typeof(op) <: NoOp
    @test transform(NoOp(), img) == img
    @test transform(NoOp(), img2) == img2
    @test transform(NoOp(), img3) == img3
end

@testset "FlipX" begin
    @test FlipX <: ImageTransformation
    op = FlipX()
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 1
    @test typeof(op) <: FlipX
    op = FlipX(0.7)
    @test multiplier(op) == 2
    @test typeof(op) <: Either
    @test_approx_eq op.chance [0.7, .3]
    @test typeof(last(transform(FlipX(), img)["transformations"])) <: FlipX
    @test transform(FlipX(0), img) == img
    @test transform(FlipX(1), img) == flipdim(img, "x")
    @imagetest "FlipX" transform(FlipX(), testimg)
end

@testset "FlipY" begin
    @test FlipY <: ImageTransformation
    op = FlipY()
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 1
    @test typeof(op) <: FlipY
    op = FlipY(0.7)
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 2
    @test typeof(op) <: Either
    @test_approx_eq op.chance [0.7, .3]
    @test typeof(last(transform(FlipY(), img)["transformations"])) <: FlipY
    @test transform(FlipY(0), img) == img
    @test transform(FlipY(1), img) == flipdim(img, "y")
    @imagetest "FlipY" transform(FlipY(), testimg)
end

@testset "Rotate90" begin
    @test Rotate90 <: ImageTransformation
    op = Rotate90()
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 1
    @test typeof(op) <: Rotate90
    op = Rotate90(0.7)
    @test multiplier(op) == 2
    @test typeof(op) <: Either
    @test_approx_eq op.chance [0.7, .3]
    @test typeof(last(transform(Rotate90(), img)["transformations"])) <: Rotate90
    @test transform(Rotate90(0), img) == img
    @test transform(Rotate90(1), img) == rotate_expand(img, deg2rad(90))
    @imagetest "Rotate90" transform(Rotate90(), testimg)
end

@testset "Rotate180" begin
    @test Rotate180 <: ImageTransformation
    op = Rotate180()
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 1
    @test typeof(op) <: Rotate180
    op = Rotate180(0.7)
    @test multiplier(op) == 2
    @test typeof(op) <: Either
    @test_approx_eq op.chance [0.7, .3]
    @test typeof(last(transform(Rotate180(), img)["transformations"])) <: Rotate180
    @test transform(Rotate180(0), img) == img
    @test transform(Rotate180(1), img) == rotate_expand(img, deg2rad(180))
    @imagetest "Rotate180" transform(Rotate180(), testimg)
end

@testset "Rotate270" begin
    @test Rotate270 <: ImageTransformation
    op = Rotate270()
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 1
    @test typeof(op) <: Rotate270
    op = Rotate270(0.7)
    @test multiplier(op) == 2
    @test typeof(op) <: Either
    @test_approx_eq op.chance [0.7, .3]
    @test typeof(last(transform(Rotate270(), img)["transformations"])) <: Rotate270
    @test transform(Rotate270(0), img) == img
    @test transform(Rotate270(1), img) == rotate_expand(img, deg2rad(270))
    @imagetest "Rotate270" transform(Rotate270(), testimg)
end

@testset "Resize" begin
    @test Resize <: ImageTransformation
    @test multiplier(Resize()) == 1
    op = Resize()
    show(op); println()
    showcompact(op); println()
    @test op.width == 64
    @test op.height == 64
    op = Resize(width = 23, height = 12)
    @test op.width == 23
    @test op.height == 12
    @test typeof(last(transform(Resize(), img)["transformations"])) <: Resize
    @test size(transform(op, img)) == (23, 12)
    @imagetest "Resize" transform(Resize(160, 80), testimg)
    @imagetest "Resize_y" permutedims(transform(Resize(160, 80), permutedims(testimg, [2, 1])), [2, 1])
end

@testset "CropRatio" begin
    @test CropRatio <: ImageTransformation
    @test multiplier(CropRatio()) == 1
    op = CropRatio(1//2)
    show(op); println()
    showcompact(op); println()
    @test op.ratio == .5
    op = CropRatio(1/2)
    @test op.ratio == .5
    op = CropRatio(ratio = 1)
    @test op.ratio == 1.
    op = CropRatio(1)
    @test op.ratio == 1.
    @test typeof(last(transform(CropRatio(), img)["transformations"])) <: CropRatio

    @test_throws ArgumentError transform(CropRatio(0.), img)
    @test_throws ArgumentError transform(CropRatio(-1.), img)

    @test transform(CropRatio(1), img)  == img
    @test transform(CropRatio(2), img)  == img[:, 1:1]
    @test transform(CropRatio(10), img) == img[:, 1:1]
    @test transform(CropRatio(.5), img) == img[1:1, :]
    @test transform(CropRatio(.1), img) == img[1:1, :]

    @test transform(CropRatio(1), img3)  == img3
    @test transform(CropRatio(2), img3)  == img3[:, 2:2]
    @test transform(CropRatio(10), img3) == img3[:, 2:2]
    @test transform(CropRatio(.5), img3) == img3[2:2, :]
    @test transform(CropRatio(.1), img3) == img3[2:2, :]

    img_1 = grayim(rand(UInt8, 50, 20))
    @test size(transform(op, img_1)) == (20, 20)
    img_2 = grayim(rand(UInt8, 10, 20))
    @test size(transform(op, img_2)) == (10, 10)
    @imagetest "CropRatio2to1" transform(CropRatio(2.), testimg)
    @imagetest "CropRatio1to2" transform(CropRatio(.5), testimg)
end

@testset "CropSize" begin
    @test CropSize <: ImageTransformation
    @test multiplier(CropSize()) == 1
    op = CropSize(2, 3)
    show(op); println()
    showcompact(op); println()
    @test op.width == 2
    @test op.height == 3
    op = CropSize(width = 4, height = 5)
    @test op.width == 4
    @test op.height == 5
    op = CropSize()
    @test op.width == 64
    @test op.height == 64
    @test typeof(last(transform(CropSize(1,1), img)["transformations"])) <: CropSize

    @test_throws ArgumentError transform(CropSize(3,2), img)
    @test_throws ArgumentError transform(CropSize(2,3), img)
    @test_throws ArgumentError transform(CropSize(0,2), img)
    @test_throws ArgumentError transform(CropSize(2,0), img)
    @test_throws ArgumentError transform(CropSize(-1,2), img)
    @test_throws ArgumentError transform(CropSize(2,-1), img)

    @test transform(CropSize(2,2), img)  == img
    @test transform(CropSize(1,1), img)  == img[2:2, 2:2]
    @test transform(CropSize(3,3), img3) == img3
    @test transform(CropSize(1,3), img3) == img3[2:2, :]
    @test transform(CropSize(3,1), img3) == img3[:, 2:2]
    @test transform(CropSize(2,2), img3) == img3[1:2, 1:2]
    @test transform(CropSize(1,1), img3) == img3[2:2, 2:2]

    op = CropSize(20, 10)
    img_1 = grayim(rand(UInt8, 50, 20))
    @test size(transform(op, img_1)) == (20, 10)
    img_2 = grayim(rand(UInt8, 30, 40))
    @test size(transform(op, img_2)) == (20, 10)
    @imagetest "CropSize" transform(CropSize(64, 32), testimg)
end

@testset "RCropSize" begin
    @test RCropSize <: ImageTransformation
    @test multiplier(RCropSize()) == 1
    op = RCropSize(2, 3)
    show(op); println()
    showcompact(op); println()
    @test op.width == 2
    @test op.height == 3
    op = RCropSize(width = 4, height = 5)
    @test op.width == 4
    @test op.height == 5
    op = RCropSize()
    @test op.width == 64
    @test op.height == 64
    @test typeof(last(transform(RCropSize(1,1), img)["transformations"])) <: RCropSize

    @test_throws ArgumentError transform(RCropSize(3,2), img)
    @test_throws ArgumentError transform(RCropSize(2,3), img)
    @test_throws ArgumentError transform(RCropSize(0,2), img)
    @test_throws ArgumentError transform(RCropSize(2,0), img)
    @test_throws ArgumentError transform(RCropSize(-1,2), img)
    @test_throws ArgumentError transform(RCropSize(2,-1), img)

    @test size(transform(RCropSize(1,1), img)) == (1, 1)
    @test size(transform(RCropSize(1,2), img)) == (1, 2)
    @test size(transform(RCropSize(2,1), img)) == (2, 1)
    @test size(transform(RCropSize(2,2), img)) == (2, 2)

    @test size(transform(RCropSize(1,1), img3)) == (1, 1)
    @test size(transform(RCropSize(1,2), img3)) == (1, 2)
    @test size(transform(RCropSize(2,1), img3)) == (2, 1)
    @test size(transform(RCropSize(2,2), img3)) == (2, 2)
    @test size(transform(RCropSize(1,3), img3)) == (1, 3)
    @test size(transform(RCropSize(3,1), img3)) == (3, 1)
    @test size(transform(RCropSize(2,3), img3)) == (2, 3)
    @test size(transform(RCropSize(3,2), img3)) == (3, 2)
    @test size(transform(RCropSize(3,3), img3)) == (3, 3)

    op = RCropSize(20, 10)
    img_1 = grayim(rand(UInt8, 50, 20))
    @test size(transform(op, img_1)) == (20, 10)
    img_2 = grayim(rand(UInt8, 30, 40))
    @test size(transform(op, img_2)) == (20, 10)
end

@testset "RCropRatio" begin
    @test RCropRatio <: ImageTransformation
    @test multiplier(RCropRatio()) == 1
    op = RCropRatio(1//2)
    show(op); println()
    showcompact(op); println()
    @test op.ratio == .5
    op = RCropRatio(1/2)
    @test op.ratio == .5
    op = RCropRatio(ratio = 1)
    @test op.ratio == 1.
    op = RCropRatio(1)
    @test op.ratio == 1.
    @test typeof(last(transform(RCropRatio(), img)["transformations"])) <: RCropRatio

    @test_throws ArgumentError transform(RCropRatio(0.), img)
    @test_throws ArgumentError transform(RCropRatio(-1.), img)

    @test size(transform(RCropRatio(1), img)) == (2, 2)
    @test size(transform(RCropRatio(2), img)) == (2, 1)
    @test size(transform(RCropRatio(10), img)) == (2, 1)
    @test size(transform(RCropRatio(.5), img)) == (1, 2)
    @test size(transform(RCropRatio(.1), img)) == (1, 2)

    @test size(transform(RCropRatio(1), img3)) == (3, 3)
    @test size(transform(RCropRatio(1.5), img3)) == (3, 2)
    @test size(transform(RCropRatio(2), img3)) == (3, 1)
    @test size(transform(RCropRatio(10), img3)) == (3, 1)
    @test size(transform(RCropRatio(.7), img3)) == (2, 3)
    @test size(transform(RCropRatio(.5), img3)) == (1, 3)
    @test size(transform(RCropRatio(.1), img3)) == (1, 3)

    img_1 = grayim(rand(UInt8, 50, 20))
    @test size(transform(op, img_1)) == (20, 20)
    img_2 = grayim(rand(UInt8, 10, 20))
    @test size(transform(op, img_2)) == (10, 10)
end

@testset "Crop" begin
    @test Crop <: ImageTransformation
    @test multiplier(Crop()) == 1
    op = Crop(1, 2, 3, 4)
    show(op); println()
    showcompact(op); println()
    @test op.x == 1
    @test op.y == 2
    @test op.width == 3
    @test op.height == 4
    op = Crop(x = 2, y = 3, width = 4, height = 5)
    @test op.x == 2
    @test op.y == 3
    @test op.width == 4
    @test op.height == 5
    op = Crop()
    @test op.x == 1
    @test op.y == 1
    @test op.width == 64
    @test op.height == 64
    @test typeof(last(transform(Crop(1,1,1,1), img)["transformations"])) <: Crop

    @test_throws ArgumentError transform(Crop(0,1,1,1), img)
    @test_throws ArgumentError transform(Crop(1,0,1,1), img)
    @test_throws ArgumentError transform(Crop(1,1,0,1), img)
    @test_throws ArgumentError transform(Crop(1,1,1,0), img)
    @test_throws ArgumentError transform(Crop(3,1,1,1), img)
    @test_throws ArgumentError transform(Crop(1,3,1,1), img)
    @test_throws ArgumentError transform(Crop(1,1,3,1), img)
    @test_throws ArgumentError transform(Crop(1,1,1,3), img)

    @test transform(Crop(1,1,1,1), img) == img[1:1, 1:1]
    @test transform(Crop(2,2,1,1), img) == img[2:2, 2:2]
    @test transform(Crop(1,1,2,2), img) == img
    @test transform(Crop(1,1,2,1), img) == img[1:2, 1:1]
    @test transform(Crop(1,1,1,2), img) == img[1:1, 1:2]
    @test transform(Crop(2,2,1,1), img3) == img3[2:2, 2:2]

    op = Crop(10, 5, 20, 10)
    img_1 = grayim(rand(UInt8, 50, 20))
    @test transform(op, img_1) == img_1[10:29, 5:14]
    @test size(transform(op, img_1)) == (20, 10)
    img_2 = grayim(rand(UInt8, 30, 40))
    @test transform(op, img_2) == img_2[10:29, 5:14]
    @test size(transform(op, img_2)) == (20, 10)
    @imagetest "Crop" transform(Crop(45, 10, 64, 32), testimg)
end

@testset "Zoom" begin
    @test Zoom <: ImageTransformation
    @test multiplier(Zoom()) == 1
    op = Zoom(1//2)
    show(op); println()
    showcompact(op); println()
    @test op.factor == .5
    op = Zoom(1/2)
    @test op.factor == .5
    op = Zoom(factor = 1)
    @test op.factor == 1.
    op = Zoom(1)
    @test op.factor == 1.
    @test typeof(last(transform(Zoom(), img)["transformations"])) <: Zoom

    @imagetest "Zoom08" transform(Zoom(0.8), testimg)
    @imagetest "Zoom12" transform(Zoom(1.2), testimg)
end

@testset "Scale" begin
    @test Scale <: ImageTransformation
    @test multiplier(Scale()) == 1
    op = Scale(2, .5)
    show(op); println()
    showcompact(op); println()
    @test op.width == 2
    @test op.height == .5
    op = Scale()
    @test op.width == 1
    @test op.height == 1
    op = Scale(width = 2)
    @test op.width == 2
    @test op.height == 1
    op = Scale(height = 2)
    @test op.width == 1
    @test op.height == 2
    op = Scale(width = 1.5, height = 0.5)
    @test op.width == 1.5
    @test op.height == 0.5
    @test typeof(last(transform(Scale(), img)["transformations"])) <: Scale

    @imagetest "Scale_x" transform(Scale(0.8, 1.2), testimg)
    @imagetest "Scale_y" permutedims(transform(Scale(0.8, 1.2), permutedims(testimg, [2, 1])), [2, 1])
end

@testset "Either" begin
    @test Either <: ImageTransformation
    @test_throws ArgumentError Either()
    op = Either(NoOp())
    show(op); println()
    showcompact(op); println()
    @test multiplier(op) == 1
    @test op.transformations == [NoOp()]
    @test op.chance == [1.]
    @test op.cum_chance == [1.]
    @test transform(op, img) == img
    op = Either(FlipX(), FlipY())
    @test multiplier(op) == 2
    @test op.transformations == [FlipX(), FlipY()]
    @test op.chance == [.5, .5]
    @test op.cum_chance == [.5, 1.]
    op = Either([FlipX(), FlipY()])
    @test op.transformations == [FlipX(), FlipY()]
    @test op.chance == [.5, .5]
    @test op.cum_chance == [.5, 1.]
    @test_throws AssertionError Either(FlipX(), FlipY(), chance = [0., 0.])
    @test_throws ArgumentError Either(FlipX(), FlipY(), Rotate90(), chance = [.5, .5])
    op = Either(FlipX(), FlipY(), Rotate90())
    @test multiplier(op) == 3
    @test op.transformations == [FlipX(), FlipY(), Rotate90()]
    @test op.chance == [1/3, 1/3, 1/3]
    @test op.cum_chance == [1/3, 2/3, 1.]
    op = Either([FlipX(), FlipY(), Rotate90()])
    @test op.transformations == [FlipX(), FlipY(), Rotate90()]
    @test op.chance == [1/3, 1/3, 1/3]
    @test op.cum_chance == [1/3, 2/3, 1.]
    op = Either(FlipX(), FlipY(), Rotate90(), chance = [1., 3, 1])
    @test op.transformations == [FlipX(), FlipY(), Rotate90()]
    @test op.chance == [.2, .6, .2]
    @test op.cum_chance == [.2, .8, 1.]
    op = Either(FlipX(), FlipY(), chance = [1., 0.])
    @test multiplier(op) == 1
    @test op.transformations == [FlipX(), FlipY()]
    @test op.chance == [1., 0]
    @test op.cum_chance == [1., 1.]
    @test transform(op, deepcopy(img)) == transform(FlipX(), deepcopy(img))
    op = Either(FlipX(), FlipY(), chance = [0., 1.])
    @test multiplier(op) == 1
    @test op.transformations == [FlipX(), FlipY()]
    @test op.chance == [0., 1]
    @test op.cum_chance == [0., 1.]
    @test transform(op, deepcopy(img)) == transform(FlipY(), deepcopy(img))
end

@testset "RandomDisplacement" begin
    @test RandomDisplacement <: ImageTransformation
    @test multiplier(RandomDisplacement(10,10)) == 1
    op = RandomDisplacement(4,5)
    @test op == RandomDisplacement(4., 5.)
    show(op); println()
    showcompact(op); println()
    @test op.gridwidth == 4
    @test op.gridheight == 5
    @test op.scale == .2
    @test op.static_border == true
    @test op.normalize == true

    @test typeof(last(transform(op, testimg)["transformations"])) <: DisplacementField

    op = RandomDisplacement(10,12, scale = .4, normalize = false)
    @test op.gridwidth == 10
    @test op.gridheight == 12
    @test op.scale == .4
    @test op.static_border == true
    @test op.normalize == false
    op = RandomDisplacement(10,12, scale = .4, static_border = false)
    @test op.gridwidth == 10
    @test op.gridheight == 12
    @test op.scale == .4
    @test op.static_border == false
    @test op.normalize == true

    @test_throws ArgumentError RandomDisplacement(2,5)
    @test_throws ArgumentError RandomDisplacement(4,2)

    srand(123)
    @imagetest "RandomDisplacement" transform(RandomDisplacement(4,5), testimg)
end

@testset "SmoothedDisplacement" begin
    @test SmoothedDisplacement <: ImageTransformation
    @test multiplier(SmoothedDisplacement(10,10)) == 1
    op = SmoothedDisplacement(4,5)
    @test op == SmoothedDisplacement(4., 5.)
    show(op); println()
    showcompact(op); println()
    @test op.gridwidth == 4
    @test op.gridheight == 5
    @test op.scale == .2
    @test op.sigma == 2.
    @test op.iterations == 1
    @test op.static_border == true
    @test op.normalize == true

    @test typeof(last(transform(op, testimg)["transformations"])) <: DisplacementField

    op = SmoothedDisplacement(10,12, scale = .4, sigma = 4, iterations = 9, normalize = false)
    @test op.gridwidth == 10
    @test op.gridheight == 12
    @test op.scale == .4
    @test op.sigma == 4.
    @test op.iterations == 9
    @test op.static_border == true
    @test op.normalize == false
    op = SmoothedDisplacement(10,12, scale = 2, sigma = 3, iterations = 2, static_border = false)
    @test op.gridwidth == 10
    @test op.gridheight == 12
    @test op.scale == 2.
    @test op.sigma == 3.
    @test op.iterations == 2
    @test op.static_border == false
    @test op.normalize == true

    @test_throws ArgumentError SmoothedDisplacement(2,5)
    @test_throws ArgumentError SmoothedDisplacement(4,2)
    @test_throws ArgumentError SmoothedDisplacement(4,5, sigma = 0.)
    @test_throws ArgumentError SmoothedDisplacement(4,5, sigma = -1.)
    @test_throws ArgumentError SmoothedDisplacement(4,5, iterations = 0)

    srand(123)
    @imagetest "SmoothedRandomDisplacement" transform(SmoothedDisplacement(4,5), testimg)
end

@testset "Tuple of Image" begin
    imgs = (img, img2)
    op = FlipX()
    out1, out2 = transform(op, imgs)
    @test typeof(last(out1["transformations"])) <: FlipX
    @test typeof(last(out2["transformations"])) <: FlipX
    @test out1 == flipx(img)
    @test out2 == flipx(img2)
    op = FlipX(1)
    out1, out2 = transform(op, imgs)
    @test out1 == flipx(img)
    @test out2 == flipx(img2)
end


text:
module DataStructs

import Base: <, <=, ==, length, isempty, iterate, show,
             last, first, eltype, map, filter, size,
             maximum, minimum

include("list.jl")
# Write your package code here.
export LinkedList, Nil, Cons, nil, cons, head, tail, list

end


text:
using ParameterSchedulers
using Flux
using LinearAlgebra
using Random
using Zygote
using BSON: @save
using BSON: @load
using CUDA
using Printf

ndigit = 2

struct CausalSelfAttention
  key::Dense
  query::Dense
  value::Dense
  attn_drop::Dropout
  resid_drop::Dropout
  proj::Dense
  n_head::Int64
#   mask::Any
end
function CausalSelfAttention(config::Dict)
  n = config["n_embed"]
  key = Dense(n, n)
  query = Dense(n, n)
  value = Dense(n, n)
  # regularization
  attn_drop = Dropout(config["attn_pdrop"])
  resid_drop = Dropout(config["resid_pdrop"])
  # output projection
  proj = Dense(n, n)

CausalSelfAttention(key, query, value, attn_drop, resid_drop, proj, config["n_head"])
end

function (m::CausalSelfAttention)(x)
  C, T, B = size(x)

  perhead = Int64(C ÷ m.n_head)

  kj = m.key(x)
  kj = reshape(kj, m.n_head, perhead, T, B)
  kj = permutedims(kj, (1, 3, 2, 4))

  qj = m.query(x)
  qj = reshape(qj, m.n_head, perhead, T, B)
  qj = permutedims(qj, (1, 3, 2, 4))

  vj = m.value(x)
  vj = reshape(vj, m.n_head, perhead, T, B)
  vj = permutedims(vj, (1, 3, 2, 4))

  att1j = reshape(
    batched_mul(batched_transpose(reshape(kj, m.n_head, T, :)), reshape(qj, m.n_head, T, :)) ./
    Float32(sqrt(size(kj, 1))),
    (T, T, perhead, B),
  )

  mask =
    (1 .- reshape(triu(ones(Float32, T, T)), Val(4))) .* (-1f9) |>
    gpu
  att1j = att1j .+ mask
  att1j = Flux.softmax(att1j)
  att1j = m.attn_drop(att1j)

  y1j = reshape(
    batched_mul(reshape(vj, m.n_head, T, :), reshape(att1j, T, T, :)),
    (m.n_head, T, perhead, B),
  )
  y1j = reshape(permutedims(y1j, (1, 3, 2, 4)), :, T, B)
  y1j = m.proj(y1j)
  y1j = m.resid_drop(y1j)
end

Flux.@functor CausalSelfAttention

Flux.trainable(a::CausalSelfAttention) = (
  a.key.weight,
  a.key.bias,
  a.query.weight,
  a.query.bias,
  a.value.weight,
  a.value.bias,
  a.proj.weight,
  a.proj.bias,
)
decayed_trainable(a::CausalSelfAttention) =
  (a.key.weight, a.query.weight, a.value.weight, a.proj.weight)

struct Block
  ln1::LayerNorm
  ln2::LayerNorm
  attn::CausalSelfAttention
  d1::Dense
  d2::Dense
  drop::Dropout
  mlp::Chain
end

function Block(config::Dict)
  n = config["n_embed"]
  ln1 = LayerNorm(n)
  ln2 = LayerNorm(n)

  attn = CausalSelfAttention(config)

  d1 = Dense(n, 4 * n, gelu)
  d2 = Dense(4 * n, n)
  drop = Dropout(config["resid_pdrop"])
  mlp = Chain(d1, d2, drop)
  Block(ln1, ln2, attn, d1, d2, drop, mlp)
end

function (m::Block)(x)
  out = m.ln1(x)
  out = m.attn(out)
  xplus = x .+ out

  out = m.ln2(xplus)
  out = m.mlp(out)
  return xplus .+ out
end

Flux.@functor Block
Flux.trainable(a::Block) = (
  a.ln1.diag.α,
  a.ln1.diag.β,
  a.ln2.diag.α,
  a.ln2.diag.β,
  Flux.trainable(a.attn)...,
  a.d1.weight,
  a.d1.bias,
  a.d2.weight,
  a.d2.bias,
)
decayed_trainable(a::Block) = (decayed_trainable(a.attn)..., a.d1.weight, a.d2.weight)

struct Gpt
  tok_emb::Flux.Embedding
  pos_emb::Array{Float32,3}
  dropout::Dropout
  blocks::Chain
  ln::LayerNorm
  head::Dense
end

function Gpt(config::Dict)
  tok_emb = Flux.Embedding(config["vocab_size"], config["n_embed"])
  pos_emb = zeros(Float32, (config["n_embed"], config["block_size"], 1))
  dropout = Dropout(config["embd_pdrop"])
  blocks = Chain([Block(config) for i = 1:config["n_layer"]]...)
  ln = LayerNorm(config["n_embed"])
  head = Dense(config["n_embed"], config["vocab_size"], bias = false)
  Gpt(tok_emb, pos_emb, dropout, blocks, ln, head)
end

function (m::Gpt)(idx)
  t, b = size(idx)
  token_embeddings = m.tok_emb(idx)
  position_embeddings = m.pos_emb[:, 1:t, :] |> gpu
  x = m.dropout(token_embeddings .+ position_embeddings)
  x = m.blocks(x)
  x = m.ln(x)
  return m.head(x)
end

Flux.@functor Gpt
function Flux.trainable(a::Gpt)
  p = []
  for b in a.blocks
    for e in Flux.trainable(b)
      push!(p, e)
    end
  end
  push!(p, a.tok_emb.weight, a.pos_emb, a.ln.diag.α, a.ln.diag.β, a.head.weight, a.head.bias)
    # push!(p, a.tok_emb.weight, a.ln.diag.α, a.ln.diag.β, a.head.weight, a.head.bias)
  return Tuple(p)
end

function decayed_trainable(a::Gpt)
  p = []
  for b in a.blocks
    for e in decayed_trainable(b)
      push!(p, e)
    end
  end
  push!(p, a.head.weight)
  return Tuple(p)
end

function dataItem(ix, ndigit = 2)
  nd = 10^ndigit
  a = ix ÷ nd
  b = ix % nd
  c = a + b

  render = lpad(a, 2, "0") * lpad(b, 2, "0") * lpad(c, 2 + 1, "0")

  dix = [parse(Int, i) for i in render]
  x = dix[1:end-1]
  y = dix[2:end]
  y[1:(ndigit*2-1)] .= -100
  x, y
end

function makeData(ndigit = 2, test_pct = 0.2)
  num_all = (10^ndigit)^2
  all_data_idx = shuffle(0:num_all-1)
  test_end_idx = (Int)(num_all * test_pct)
  test_idx = all_data_idx[1:test_end_idx]
  train_idx = all_data_idx[test_end_idx+1:end]

  test_data_x = Array{Int32,2}(undef, ndigit * 3, test_end_idx)
  test_data_y = Array{Int32,2}(undef, ndigit * 3, test_end_idx)
  for (i, ix) in enumerate(test_idx)
    x, y = dataItem(ix, ndigit)
    test_data_x[:, i] = x
    test_data_y[:, i] = y
  end

  train_data_x = Array{Int32,2}(undef, ndigit * 3, num_all - test_end_idx)
  train_data_y = Array{Int32,2}(undef, ndigit * 3, num_all - test_end_idx)
  for (i, ix) in enumerate(train_idx)
    x, y = dataItem(ix, ndigit)
    train_data_x[:, i] = x
    train_data_y[:, i] = y
  end
  return (train_data_x, train_data_y), (test_data_x, test_data_y)
end

function addOneForJulia(data)
  data .= data .+ 1
end

mutable struct MYADAM
  eta::Float64
  beta::Tuple{Float64,Float64}
  epsilon::Float64
  wd::Float64
  state::IdDict
  sched_mult::Real
end

MYADAM(η::Real, β::Tuple, wd::Real, state::IdDict) = MYADAM(η, β, Flux.Optimise.EPS, wd, state, 1f0)

function MYADAM(
  allParams,
  decayedParams,
  η::Real = 0.001f0,
  β::Tuple = (0.9f0, 0.999f0),
  wd::Real = 0.1f0,
)
  myiddict = IdDict()
  for p in allParams
    myiddict[p] = (zero(p), zero(p), Float64[β[1], β[2]], p in decayedParams ? true : false)
  end
  MYADAM(η, β, wd, myiddict)
end

function Flux.Optimise.apply!(o::MYADAM, x, Δ)
    η, β = o.eta, o.beta
    wd = o.wd

    mt, vt, βp, decay = get(o.state, x, (zero(x), zero(x), Float64[β[1], β[2]], false))

    ## Clip norm
    Δnrm = norm(Δ)
    if Δnrm > 1.0f0
        rmul!(Δ, 1.0f0 / Δnrm)
    end

    @. mt = β[1] * mt + (1 - β[1]) * Δ
    @. vt = β[2] * vt + (1 - β[2]) * Δ * conj(Δ)
    @. Δ = o.sched_mult * (mt / (1 - βp[1]) / (√(vt / (1 - βp[2])) + o.epsilon))
    βp .= βp .* β

    if decay
        @. Δ += o.sched_mult * wd * x
    end

    ## Descent
    Δ .*= o.eta

  return Δ
end

function myloss(gpt, x, allRowsButLast, tot)
  logits = gpt(x)
  loss = -sum(sum(allRowsButLast .* Flux.Losses.logsoftmax(logits, dims = 1), dims = 1)) / tot
  return loss
end

function mytraining(trnx, trny, tstx, tsty, config)
  gpt = Gpt(config) |> gpu
  allParams = Flux.params(gpt)
  decayedParams = Flux.Params()
  push!(decayedParams, decayed_trainable(gpt)...)

  lr = config["learning_rate"]
  opt = MYADAM(allParams, decayedParams, lr, config["betas"])
    # opt = ADAM(config[learning_rate], config[betas])
    # opt = Flux.Optimiser(adam, WeightDecay(0.1))

  train_loader =
    Flux.DataLoader((data = trnx, label = trny), batchsize = config["batch_size"], shuffle = true)
  test_loader =
    Flux.DataLoader((data = tstx, label = tsty), batchsize = config["batch_size"], shuffle = true)

  labels = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, -99]
  tokens = 0
  best_loss = Inf32
  for epoch = 1:config["max_epochs"]

    for (it, (x, y)) in enumerate(train_loader)
      x = x |> gpu

      allRowsButLast =
        Flux.onehotbatch(y, labels)[1:end-1, :, :] |> gpu

      total_num_of_entries = reduce(*, size(y))

      tot = total_num_of_entries - count(i -> i == -99, y) #|> gpu

      train_loss, back = Zygote.pullback(() -> myloss(gpt, x, allRowsButLast, tot), allParams)

      gs = back(one(train_loss))

      Flux.Optimise.update!(opt, allParams, gs)

      if config["lr_decay"]
        tokens += count(i -> i >= 0, y)
        if tokens < config["warmup_tokens"]
          lr_mult = Float32(tokens / max(1, config["warmup_tokens"]))
        else
          progress = Float32(
            (tokens - config["warmup_tokens"]) / max(1, config["final_tokens"] - config["warmup_tokens"]),
          )
          lr_mult = max(0.1f0, 0.5f0 * (1.0 + cos(pi * progress)))
        end
        # lr = config[learning_rate] * lr_mult
        opt.sched_mult = lr_mult

      end

      if it % 10 == 1
        @printf "Epoch: %d Iter: %d Train Loss: %.2f lr_mult: %.2f tokens: %d\n" epoch it train_loss lr_mult tokens
        # println(
        #   "Epoch: $(epoch) Iter: $(it) Train Loss: $(train_loss) lr: $(lr) lr_mult: $(lr_mult) tokens: $(tokens)",
        # )
      end
    end

    test_loss = epochOnTest(test_loader, gpt, labels)
    if test_loss < best_loss
      best_loss = test_loss
      @save "mymodel.bson" gpt
    end
  end
  return gpt
end

function epochOnTest(test_loader, model, labels)
  losses = Vector{Float32}(undef, 0)
  for (it, (x, y)) in enumerate(test_loader)
    x = x |> gpu
    allRowsButLast = Flux.onehotbatch(y, labels)[1:end-1, :, :] |> gpu
    total_num_of_entries = reduce(*, size(y))
    tot = Float32(total_num_of_entries - count(i -> i == -99, y))
    loss = myloss(model, x, allRowsButLast, tot)
    push!(losses, loss)
  end
  test_loss = Flux.mean(losses)
  println("Test Loss: $(test_loss)")
  return test_loss
end

function sample(config, model, x, steps, temperature = 1.0, sample = false, top_k = nothing)
  bs = config["block_size"]

  for k = 1:steps
    S, B = size(x)
    x_cond = nothing
    if S <= bs
      x_cond = x
    else
      x_cond = x[end-bs:end, :]
    end
    x_cond = x_cond |> gpu
    logits = model(x_cond)
    logits = logits[:, end, :] / temperature

    probs = Flux.softmax(logits)

    ix = getindex.(argmax(probs, dims = 1), 1)
    x = vcat(x, ix)
  end
  return x
end

function give_exam(gpt, x, y, config)
  data_loader =
    Flux.DataLoader((data = x, label = y), batchsize = config["batch_size"], shuffle = true)

  tot_correct = 0
  tot = 0
  for (x, y) in data_loader

    d1d2 = x[1:ndigit*2, :]
    d1d2d3 = sample(config, gpt, d1d2, ndigit + 1)
    d3 = d1d2d3[end-(ndigit):end, :]

    factors = [10^i for i = ndigit:-1:0]
    d1i = sum((d1d2.-1)[1:ndigit, :] .* factors[2:end], dims = 1)
    d2i = sum((d1d2.-1)[ndigit+1:ndigit*2, :] .* factors[2:end], dims = 1)
    d3i_pred = sum((d3 .- 1) .* factors, dims = 1)
    d3_gt = d1i .+ d2i
    correct = sum(d3_gt .== d3i_pred)
    tot += size(x, 2)
    tot_correct += correct
  end
  println("tot: $(tot) tot_correct: $(tot_correct)")
end

max_stars_repo_name: cancandan/mingpt-julia

max_issues_repo_name: cancandan/mingpt-julia

max_forks_repo_name: cancandan/mingpt-julia


text:
using CairoMakie, Distributions, Random
Random.seed!(124)
colors = ["#FF410D", "#6EE2FF", "#F7C530", "#95CC5E", "#D0DFE6", "#F79D1E"]
μσpairs = [[2,0.5], [-1,2], [0.25,1], [1,0.1], [1, 0.05], [1.2,0.1]]

with_theme(theme_dark()) do
    fig = Figure(resolution = (600,400))
    ax = Axis(fig[1,1]; palette = (; patchcolor = tuple.(colors, 0.45)))
    for (idx, μσ) in enumerate(μσpairs)
        density!(rand(Normal(μσ...), 1000), strokewidth = 1.5,
            strokecolor = :white, direction = idx > 3 ? :x : :y,
            linestyle = idx > 3 ? :dash : nothing,
            label = "$(μσ[1]),  $(μσ[2])")
    end
    axislegend(L"\mu,\quad\sigma"; position= :cb, titlesize= 22)
    hidedecorations!(ax; grid = false)
    save("density.png", fig, px_per_unit = 2)
    display(fig)
end