text:
isnormal(c::SGRColor) = isempty(c.class)

is216color(c::SGRColor) = is216color(c.hex)

function is216color(hex::AbstractString)
    hex == "000" || hex == "fff" || occursin(r"(?:00|5f|87|af|d7|ff){3}$", hex)
end

function codes(c::SGRColor)
    m = match(r"^([345]8)_([25])$", c.class)
    m === nothing && return (parse(Int, c.class),)
    code, sub = m.captures
    codei = parse(Int, code)
    if sub == "5"
        if is216color(c)
            h = parse(UInt32, c.hex, base=16)
            h === 0x00000fff && return (codei, 5, 231)
            r = (h >> 0x10) % UInt8 ÷ 0x30
            g = (h >> 0x08) % UInt8 ÷ 0x30
            b = (h >> 0x00) % UInt8 ÷ 0x30
            return (codei, 5, (r * 0x24 + g * 0x6 + b) + 16)
        else
            h = parse(UInt8, c.hex[1:2], base=16)
            g = (h - 0x8) ÷ 0xa
            return (codei, 5, g + 232)
        end
    else
        if length(c.hex) == 3
            h = parse(UInt16, c.hex, base=16)
            r = (h >> 0x8)
            g = (h >> 0x4) & 0xf
            b = h & 0xf
            return (codei, 2, r * 17, g * 17, b * 17)
        else
            h = parse(UInt32, c.hex, base=16)
            r = (h >> 0x10)
            g = (h >> 0x8) & 0xff
            b = h & 0xff
            return (codei, 2, Int(r), Int(g), Int(b))
        end
    end
end

function short_hex(r::UInt8, g::UInt8, b::UInt8)
    rgb6 = UInt32(r) << 0x10 + UInt32(g) << 0x8 + b
    rgb6 === (rgb6 & 0x0f0f0f) * 0x11 || return string(rgb6, pad=6, base=16)
    string(UInt16(r >> 0x4) << 0x8 + UInt16(g >> 0x4) << 0x4 + b >> 0x4, pad=3, base=16)
end

function set_16colors!(ctx::SGRContext, d::AbstractString)
    if d[1] === '3' || d[1] === '9'
        ctx.fg = SGRColor(d)
    else
        ctx.bg = SGRColor(d)
    end
end

const SCALE_216 = UInt8[0x00, 0x5f, 0x87, 0xaf, 0xd7, 0xff]

function set_256colors!(ctx::SGRContext, d::AbstractString, color::AbstractString)
    fore = d[1] === '3'
    hex = ""
    colorid = isempty(color) ? 0x0 : parse(UInt8, color)
    if colorid < 0x8
        class = d[1] * string(colorid)
    elseif colorid < 0x10
        class = (fore ? "9" : "10") * string(colorid - 0x8)
    else
        if colorid < 0xe8
            c = colorid - 0x10
            r = SCALE_216[c ÷ 0x24 + 1]
            g = SCALE_216[c ÷ 0x6 % 0x6 + 1]
            b = SCALE_216[c % 0x6 + 1]
            hex = short_hex(r, g, b)
        else
            g = (colorid - 0xe8) * 0xa + 0x8
            hex = short_hex(g, g, g)
        end
        class = d * "_5"
    end

    if fore
        ctx.fg = SGRColor(class, hex)
    else
        ctx.bg = SGRColor(class, hex)
    end
end

function set_24bitcolors!(ctx::SGRContext, d::AbstractString,
                          r::AbstractString, g::AbstractString, b::AbstractString)
    r8 = isempty(r) ? 0x0 : parse(UInt8, r)
    g8 = isempty(g) ? 0x0 : parse(UInt8, g)
    b8 = isempty(b) ? 0x0 : parse(UInt8, b)
    hex = short_hex(r8, g8, b8)
    if is216color(hex) || (r8 === g8 === b8 && (r8 - 8) % 10 == 0)
        class = d * "_5"
    else
        class = d * "_2"
    end
    if d == "38"
        ctx.fg = SGRColor(class, hex)
    else
        ctx.bg = SGRColor(class, hex)
    end
end


text:
# Functions for calculating exact overlap between shapes.
#
# Original cython version by Thomas Robitaille released under a
# 3-clause BSD-style license.

# Return area of a circle arc between (x1, y1) and (x2, y2) with radius r
# reference: http://mathworld.wolfram.com/CircularSegment.html
function area_arc(x1, y1, x2, y2, r)
  a = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1))
  theta = 2. * asin(0.5 * a / r)
  return 0.5 * r * r * (theta - sin(theta))
end

# Area of a triangle defined by three verticies
area_triangle(x1, y1, x2, y2, x3, y3) =
    0.5 * abs(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))

# Core of circular overlap routine.
# Assumes that xmax >= xmin >= 0.0, ymax >= ymin >= 0.0.
# (can always modify input to conform to this).
function circoverlap_core(xmin, ymin, xmax, ymax, r)

    xmin2 = xmin * xmin
    ymin2 = ymin * ymin
    r2 = r * r
    (xmin2 + ymin2 > r2) && return 0.0

    xmax2 = xmax * xmax
    ymax2 = ymax * ymax
    (xmax2 + ymax2 < r2) && return (xmax-xmin) * (ymax-ymin)

    a = xmax2 + ymin2  # (corner 1 distance)^2
    b = xmin2 + ymax2  # (corner 2 distance)^2

    if (a < r2 && b < r2)
        x1 = sqrt(r2 - ymax2)
        y1 = ymax
        x2 = xmax
        y2 = sqrt(r2 - xmax2)
        return ((xmax - xmin) * (ymax - ymin) -
                area_triangle(x1, y1, x2, y2, xmax, ymax) +
	            area_arc(x1, y1, x2, y2, r))

    elseif (a < r2)
        x1 = xmin
        y1 = sqrt(r2 - xmin2)
        x2 = xmax
        y2 = sqrt(r2 - xmax2)
        return (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, x1, ymin, xmax, ymin) +
                area_triangle(x1, y1, x2, ymin, x2, y2))

    elseif (b < r2)
        x1 = sqrt(r2 - ymin2)
        y1 = ymin
        x2 = sqrt(r2 - ymax2)
        y2 = ymax
        return (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, xmin, y1, xmin, ymax) +
                area_triangle(x1, y1, xmin, y2, x2, y2))

    else
        x1 = sqrt(r2 - ymin2)
        y1 = ymin
        x2 = xmin
        y2 = sqrt(r2 - xmin2)
        return (area_arc(x1, y1, x2, y2, r) +
                area_triangle(x1, y1, x2, y2, xmin, ymin))
    end
end

# Area of overlap of a rectangle and a circle
function circoverlap(xmin, ymin, xmax, ymax, r)

    # some subroutines demand that r > 0
    (r <= 0.0) && return 0.0

    if (0.0 <= xmin)
        (0.0 <= ymin) && return circoverlap_core(xmin, ymin, xmax, ymax, r)
        (0.0 >= ymax) && return circoverlap_core(-ymax, xmin, -ymin, xmax, r)
        return (circoverlap(xmin, ymin, xmax, 0.0, r) +
                circoverlap(xmin, 0.0, xmax, ymax, r))

    elseif (0.0 >= xmax)
        (0.0 <= ymin) && return circoverlap_core(-xmax, ymin, -xmin, ymax, r)
        (0.0 >= ymax) && return circoverlap_core(-xmax, -ymax, -xmin, -ymin, r)
        return (circoverlap(xmin, ymin, xmax, 0.0, r) +
                circoverlap(xmin, 0.0, xmax, ymax, r))

    else
        (0.0 <= ymin) && return (circoverlap(xmin, ymin, 0.0, ymax, r) +
                                 circoverlap(0.0, ymin, xmax, ymax, r))
        (0.0 >= ymax) && return (circoverlap(xmin, ymin, 0.0, ymax, r) +
                                 circoverlap(0.0, ymin, xmax, ymax, r))
        return (circoverlap(xmin, ymin, 0.0, 0.0, r) +
                circoverlap(0.0, ymin, xmax, 0.0, r) +
                circoverlap(xmin, 0.0, 0.0, ymax, r) +
                circoverlap(0.0, 0.0, xmax, ymax, r))
    end
end

# -----------------------------------------------------------------------------
# ellipse-rectangle overlap

# internal types
immutable Point
    x::Float64
    y::Float64
end

immutable Intersections
    p1::Point
    p2::Point
end

# rotate -> a, b, c = b, c, a

# Check if a point (x,y) is inside a triangle
function in_triangle(x, y, x1, y1, x2, y2, x3, y3)
    c = 0
    c += ((y1 > y) != (y2 > y)) && (x < (x2-x1) * (y-y1) / (y2-y1) + x1)
    c += ((y2 > y) != (y3 > y)) && (x < (x3-x2) * (y-y2) / (y3-y2) + x2)
    c += ((y3 > y) != (y1 > y)) && (x < (x1-x3) * (y-y3) / (y1-y3) + x3)

    return c % 2 == 1
end

# Intersection of a line defined by two points with a unit circle
function circle_line(x1, y1, x2, y2)

    tol = 1.e-10

    dx = x2 - x1
    dy = y2 - y1

    if abs(dx) < tol && abs(dy) < tol
        return Intersections(Point(2.0, 2.0), Point(2.0, 2.0))

    elseif abs(dx) > abs(dy)

        # Find the slope and intercept of the line
        a = dy / dx
        b = y1 - a * x1

        # Find the determinant of the quadratic equation
        delta = 1.0 + a*a - b*b
        if (delta > 0.)  # solutions exist
	    delta = sqrt(delta)
            p1x = (-a*b - delta) / (1. + a*a)
            p2x = (-a*b + delta) / (1. + a*a)
            return Intersections(Point(p1x, a * p1x + b),
                                 Point(p2x, a * p2x + b))
        else             # no solutions exist
            return Intersections(Point(2.0, 2.0), Point(2.0, 2.0))
        end

    else

        # Find the slope and intercept of the line
        a = dx / dy
        b = x1 - a * y1

        # Find the determinant of the quadratic equation
        delta = 1.0 + a*a - b*b
        if (delta > 0.)  # solutions exist
	    delta = sqrt(delta)
	    p1y = (-a*b - delta) / (1. + a*a)
            p2y = (-a*b + delta) / (1.0 + a*a)
            return Intersections(Point(a * p1y + b, p1y),
                                 Point(a * p2y + b, p2y))
        else  # no solution, return values > 1
            return Intersections(Point(2.0, 2.0), Point(2.0, 2.0))
        end
    end
end

# The intersection of a line with the unit circle. The intersection
# closest to (x2, y2) is chosen.
function circle_segment_single2(x1, y1, x2, y2)
    inter = circle_line(x1, y1, x2, y2)

    dx1 = abs(inter.p1.x - x2)
    dy1 = abs(inter.p1.y - y2)
    dx2 = abs(inter.p2.x - x2)
    dy2 = abs(inter.p2.y - y2)

    if (dx1 > dy1)  # compare based on x-axis
        return (dx1 > dx2) ? inter.p2 : inter.p1
    else
        return (dy1 > dy2) ? inter.p2 : inter.p1
    end
end

# Intersection(s) of a segment with the unit circle. Discard any
# solution not on the segment.
function circle_segment(x1, y1, x2, y2)
    inter = circle_line(x1, y1, x2, y2)
    pt1 = inter.p1
    pt2 = inter.p2

    if ((pt1.x > x1 && pt1.x > x2) || (pt1.x < x1 && pt1.x < x2) ||
        (pt1.y > y1 && pt1.y > y2) || (pt1.y < y1 && pt1.y < y2))
        pt1 = Point(2.0, 2.0)
    end

    if ((pt2.x > x1 && pt2.x > x2) || (pt2.x < x1 && pt2.x < x2) ||
        (pt2.y > y1 && pt2.y > y2) || (pt2.y < y1 && pt2.y < y2))
        pt2 = Point(2.0, 2.0)
    end

    if pt1.x > 1. && pt2.x < 2.
        return Intersections(pt1, pt2)
    else
        return Intersections(pt2, pt1)
    end
end

# Given a triangle defined by three points (x1, y1), (x2, y2), and
# (x3, y3), find the area of overlap with the unit circle.
function triangle_unitcircle_overlap(x1, y1, x2, y2, x3, y3)

    # Find distance of all vertices to circle center
    d1 = x1*x1 + y1*y1
    d2 = x2*x2 + y2*y2
    d3 = x3*x3 + y3*y3

    # Order vertices by distance from origin
    if (d1 < d2)

        if (d2 < d3)

        elseif (d1 < d3)
	    x2, x3 = x3, x2
            y2, y3 = y3, y2
            d2, d3 = d3, d2

        else
            # rotate -> a, b, c = b, c, a
            x1, x2, x3 = x3, x1, x2
            y1, y2, y3 = y3, y1, y2
            d1, d2, d3 = d3, d1, d2
        end
    else
        if d1 < d3
            x1, x2 = x2, x1
            y1, y2 = y2, y1
            d1, d2 = d2, d1
        elseif d2 < d3
            x1, x2, x3 = x2, x3, x1
            y1, y2, y3 = y2, y3, y1
            d1, d2, d3 = d2, d3, d1
        else
            x1, x3 = x3, x1
            y1, y3 = y3, y1
            d1, d3 = d3, d1
        end
    end

    # Determine number of vertices inside circle
    in1 = d1 < 1.0
    in2 = d2 < 1.0
    in3 = d3 < 1.0

    # Determine which vertices are on the circle
    on1 = abs(d1 - 1.0) < 1e-10
    on2 = abs(d2 - 1.0) < 1e-10
    on3 = abs(d3 - 1.0) < 1e-10

    if (on3 || in3)  # triangle completely within circle
        area = area_triangle(x1, y1, x2, y2, x3, y3)

    elseif (in2 || on2)
        # If vertex 1 or 2 are on the edge of the circle, then we use
        # the dot product to vertex 3 to determine whether an
        # intersection takes place.
        intersect13 = !on1 || (x1*(x3-x1) + y1*(y3-y1) < 0.0)
        intersect23 = !on2 || (x2*(x3-x2) + y2*(y3-y2) < 0.0)

        if (intersect13 && intersect23)
	    pt1 = circle_segment_single2(x1, y1, x3, y3)
	    pt2 = circle_segment_single2(x2, y2, x3, y3)
	    area = (area_triangle(x1, y1, x2, y2, pt1.x, pt1.y) +
		    area_triangle(x2, y2, pt1.x, pt1.y, pt2.x, pt2.y) +
		    area_arc(pt1.x, pt1.y, pt2.x, pt2.y, 1.0))
        elseif (intersect13)
	  pt1 = circle_segment_single2(x1, y1, x3, y3)
	  area = (area_triangle(x1, y1, x2, y2, pt1.x, pt1.y) +
		  area_arc(x2, y2, pt1.x, pt1.y, 1.0))
        elseif (intersect23)
	  pt2 = circle_segment_single2(x2, y2, x3, y3)
	  area = (area_triangle(x1, y1, x2, y2, pt2.x, pt2.y) +
		  area_arc(x1, y1, pt2.x, pt2.y, 1.0))
        else
	    area = area_arc(x1, y1, x2, y2, 1.0)
        end

    elseif in1
        # Check for intersections of far side with circle
        inter = circle_segment(x2, y2, x3, y3)
        pt1 = inter.p1
        pt2 = inter.p2
        pt3 = circle_segment_single2(x1, y1, x2, y2)
        pt4 = circle_segment_single2(x1, y1, x3, y3)

        if pt1.x > 1.0  # indicates no intersection

	    # check if the pixel vertex (x1, y2) and the origin are on
	    # different sides of the circle segment. If they are, the
	    # circle segment spans more than pi radians.
	    # We use the formula (y-y1) * (x2-x1) > (y2-y1) * (x-x1)
	    # to determine if (x, y) is on the left of the directed
	    # line segment from (x1, y1) to (x2, y2)
	    if (((0.-pt3.y) * (pt4.x-pt3.x) > (pt4.y-pt3.y) * (0.-pt3.x)) !=
	        ((y1-pt3.y) * (pt4.x-pt3.x) > (pt4.y-pt3.y) * (x1-pt3.x)))
	        area = (area_triangle(x1, y1, pt3.x, pt3.y, pt4.x, pt4.y) +
		        pi - area_arc(pt3.x, pt3.y, pt4.x, pt4.y, 1.0))
	    else
	        area = (area_triangle(x1, y1, pt3.x, pt3.y, pt4.x, pt4.y) +
		        area_arc(pt3.x, pt3.y, pt4.x, pt4.y, 1.0))
	    end
        else
            # ensure that pt1 is the point closest to (x2, y2)
	    if (((pt2.x-x2)*(pt2.x-x2) + (pt2.y-y2)*(pt2.y-y2)) <
	        ((pt1.x-x2)*(pt1.x-x2) + (pt1.y-y2)*(pt1.y-y2)))
                pt1, pt2 = pt2, pt2
            end
	    area = (area_triangle(x1, y1, pt3.x, pt3.y, pt1.x, pt1.y) +
		    area_triangle(x1, y1, pt1.x, pt1.y, pt2.x, pt2.y) +
		    area_triangle(x1, y1, pt2.x, pt2.y, pt4.x, pt4.y) +
		    area_arc(pt1.x, pt1.y, pt3.x, pt3.y, 1.0) +
		    area_arc(pt2.x, pt2.y, pt4.x, pt4.y, 1.0))
	end

    else
        inter = circle_segment(x1, y1, x2, y2)
        pt1 = inter.p1
        pt2 = inter.p2
        inter = circle_segment(x2, y2, x3, y3)
        pt3 = inter.p1
        pt4 = inter.p2
        inter = circle_segment(x3, y3, x1, y1)
        pt5 = inter.p1
        pt6 = inter.p2

        if (pt1.x <= 1.0)
	    xp = 0.5 * (pt1.x + pt2.x)
	    yp = 0.5 * (pt1.y + pt2.y)
	    area = (triangle_unitcircle_overlap(x1, y1, x3, y3, xp, yp) +
		    triangle_unitcircle_overlap(x2, y2, x3, y3, xp, yp))

        elseif (pt3.x <= 1.)
	    xp = 0.5 * (pt3.x + pt4.x)
	    yp = 0.5 * (pt3.y + pt4.y)
	    area = (triangle_unitcircle_overlap(x3, y3, x1, y1, xp, yp) +
		    triangle_unitcircle_overlap(x2, y2, x1, y1, xp, yp))

        elseif (pt5.x <= 1.)
	    xp = 0.5 * (pt5.x + pt6.x)
	    yp = 0.5 * (pt5.y + pt6.y)
	    area = (triangle_unitcircle_overlap(x1, y1, x2, y2, xp, yp) +
		    triangle_unitcircle_overlap(x3, y3, x2, y2, xp, yp))

        else  # no intersections
	  return in_triangle(0., 0., x1, y1, x2, y2, x3, y3) ? pi : 0.0
        end
    end

    return area
end

# exact overlap between a rectangle defined by (xmin, ymin, xmax,
# ymax) and an ellipse with major and minor axes rx and ry
# respectively and position angle theta.
function ellipoverlap(xmin, ymin, xmax, ymax, a, b, theta)
    cos_m_theta = cos(-theta)
    sin_m_theta = sin(-theta)

    scale = a * b  # areas will be shrunk by this scale

    # Reproject rectangle to a frame in which ellipse is a unit circle
    x1 = (xmin * cos_m_theta - ymin * sin_m_theta) / a
    y1 = (xmin * sin_m_theta + ymin * cos_m_theta) / b
    x2 = (xmax * cos_m_theta - ymin * sin_m_theta) / a
    y2 = (xmax * sin_m_theta + ymin * cos_m_theta) / b
    x3 = (xmax * cos_m_theta - ymax * sin_m_theta) / a
    y3 = (xmax * sin_m_theta + ymax * cos_m_theta) / b
    x4 = (xmin * cos_m_theta - ymax * sin_m_theta) / a
    y4 = (xmin * sin_m_theta + ymax * cos_m_theta) / b

    # Divide resulting quadrilateral into two triangles and find
    # intersection with unit circle
    return scale * (triangle_unitcircle_overlap(x1, y1, x2, y2, x3, y3) +
		    triangle_unitcircle_overlap(x1, y1, x4, y4, x3, y3))
end

ext: jl


text:
module SkinModHelperSkinModSwapTrigger
using ..Ahorn, Maple

@mapdef Trigger "SkinModHelper/SkinModSwapTrigger" SkinModSwapTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	skinId::String="default_skin", revertOnLeave::Bool=false)

const placements = Ahorn.PlacementDict(
	"Skin Mod Swap Trigger (Skin Mod Helper)" => Ahorn.EntityPlacement(
		SkinModSwapTrigger,
		"rectangle"
	)
)

end

max_stars_repo_name: coloursofnoise/SkinModHelper

max_issues_repo_name: coloursofnoise/SkinModHelper

max_forks_repo_name: coloursofnoise/SkinModHelper


text:
"""
meshgrid(x)
meshgrid(x, y)
meshgrid(x, y, z)

Copy of Matlab's meshgrid function.
Returns 2D grid coordinates based on x and y.

# Examples
```jldoctest
julia> R = meshgrid(1:3, 1:3)
([1 2 3; 1 2 3; 1 2 3], [1 1 1; 2 2 2; 3 3 3])
```

In 3D, meshgrid returns vectors rather than tensors.
```jldoctest
julia> meshgrid(1:3)
([1 2 3; 1 2 3; 1 2 3], [1 1 1; 2 2 2; 3 3 3])
```
"""

function meshgrid(x,y)
        X = repeat(transpose(x[:]), size(y[:], 1), 1)
        Y = repeat(y, 1, size(x, 1))
        return X, Y
end

function meshgrid(x)
        return meshgrid(x,x)
end

function meshgrid(x1D,y1D,z1D)
        Np = length(x1D)*length(y1D)*length(z1D)
        x = zeros(Np)
        y = zeros(Np)
        z = zeros(Np)
        sk = 1
        for k = 1:length(z1D)
                for j = 1:length(y1D)
                        for i = 1:length(x1D)
                                x[sk] = x1D[i]
                                y[sk] = y1D[j]
                                z[sk] = z1D[k]
                                sk += 1
                        end
                end
        end
        return x,y,z
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
module LCIOUtils
export getpositions, mappositions, mapparticles, maptruthparticles, numfaces, hcal_outer, hcal_inner, ecal_inner, ecal_outer, hcal_length, ecal_length, radiationLength, interactionLength, inertiatensor, append

using LCIO
using DataStructures

#TODO: procure the correct constants for the ILC
const hcal_outer = 2495
const hcal_inner = 1419
const hcal_length = 6036
const numfaces = 8

const ecal_inner = 1265
const ecal_outer = 1403.5
const ecal_length = 3530
const radiationLength = 25.695807
const interactionLength = 1.041129

"Append contents to the a file with the given filename."
function append(filename, contents)
    file = open(filename, "a")
    write(file, contents)
    close(file)
end

function inertiatensor(positions::Matrix)
    m = zeros(Float64, 3, 3)
    for j=1:size(positions)[2]
        pos = positions[:,j]
        m[1, 1] += (pos[2]^2 + pos[3]^2)*pos[4]
        m[2, 2] += (pos[1]^2 + pos[3]^2)*pos[4]
        m[3, 3] += (pos[1]^2 + pos[2]^2)*pos[4]
        m[1, 2] -= pos[1]*pos[2]*pos[4]
        m[2, 3] -= pos[2]*pos[3]*pos[4]
        m[3, 1] -= pos[1]*pos[3]*pos[4]
    end
    m[2, 1] = m[1, 2]
    m[3, 2] = m[2, 3]
    m[1, 3] = m[3, 1]

    inertia = Symmetric(m)
    return inertia
end

function getParticleHits(p)
    hitList = CalHit[]
    for c in getClusters(p)
        for h in getCalorimeterHits(c)
            push!(hitList, CalHit(h))
        end
    end
    hitList
end

"""
  getPositions(fileName, collectionName, eventNum=1)

Get the position/energy matrix of calorimeter hits in the slcio file, where each output column is [x, y, z, E].
The collection to read from is `collectionName`, and the event index is `eventNum`.
Currently only works for PandoraPFOCollection.
"""
function getpositions(filename, collectionName="PandoraPFOCollection", eventNum=1)
    LCIO.open(filename) do reader
        for (idx, event) in zip(1:eventNum, reader)
            if (idx == eventNum)
                for name in LCIO.getCollectionNames(event)
                    if (name == collectionName)
                        # collection = LCIO.getCollection(event, name)
                        collection = LCIO.getCollection(event, collectionName)
                        positions = Array(Float64, 4, 0)

                        for (hitnum, c) in enumerate(collection)
                            for hit in getParticleHits(c)
                                positions = hcat(positions, [hit.x, hit.y, hit.z, hit.E])
                            end
                        end
                        return positions
                    end
                end
                println("No collection named $(collectionName)")
                throw(KeyError(collectionName))
            end
        end
    end
end

"""
  mapcollections(f, filename, startevent, endevent, collections...)

Map a function f(positionvec4, i) over the specified collections from the specified
start to end event.
"""
function mapcollections(f, filename, startevent, endevent, collections...)
    LCIO.open(filename) do reader
        for (i, event) in zip(1:endevent, reader)
            if i < startevent continue end
            positions = nil(Vector{Float64})
            for name in LCIO.getCollectionNames(event)
                for tname in collections
                    if name == tname
                        collection = LCIO.getCollection(event, tname)
                        for hit in collection
                            positions = cons([hit.x, hit.y, hit.z, hit.E], positions)
                        end
                    end
                end
            end
            f(hcat(positions...), i)
        end
    end
end

"""
  mappositions(f, filename, startevent, endevent, collectionname="PandoraPFOCollection")

Applies the function `f(positionvec4, i, maxenergy, momentum, particletype)` to the matrix of position-energy 4-vectors,
event number `i`, reconstructed energy `maxenergy`, 3-momentum, and particle type for events `startevent:endevent`.
`positionvec4` is a 4xn matrix where n is the number of hits, and the last component is the energy.
Currently only works for PandoraPFOCollection. Only considers one particle per event (with maximum energy)
"""
function mappositions(f, filename, startevent, endevent, collectionname="PandoraPFOCollection")
    LCIO.open(filename) do reader
        for (i, event) in zip(1:endevent, reader)
            if i < startevent continue end
            for name in LCIO.getCollectionNames(event)
                if name == collectionname
                    collection = LCIO.getCollection(event, collectionname)
                    positions = nil(Vector{Float64})
                    maxEnergy = 0.0
                    momentum = Float64[0, 0, 0]
                    particletype = 0
                    for c in collection
                        p4vec = LCIO.getP4(c)
                        if p4vec.t > maxEnergy
                            particletype = getType(c)
                            maxEnergy = p4vec.t
                            momentum = [p4vec.x, p4vec.y, p4vec.z]
                        end
                        for hit in getParticleHits(c)
                            positions = cons([hit.x, hit.y, hit.z, hit.E], positions)
                        end
                    end
                    f(hcat(positions...), i, maxEnergy, momentum, particletype)
                end
            end
        end
    end
end

"""
  mapparticles(f, filename, startevent, endevent, collectionname="PandoraPFOCollection")

Applies the function `f(positionvec4, i, maxenergy, momentum, particletype)` to the matrix of position-energy 4-vectors,
event number `i`, reconstructed energy `maxenergy`, momentum, and particle type for all particles in events `startevent:endevent`.
`positionvec4` is a 4xn matrix where n is the number of hits, and the last component is the energy.
Currently only works for PandoraPFOCollection.
"""
function mapparticles(f, filename, startevent, endevent, collectionname="PandoraPFOCollection")
    LCIO.iterate(filename) do reader
        for (i, event) in zip(1:endevent, reader)
            if i < startevent continue end
            for name in LCIO.getCollectionNames(event)
                if name == collectionname
                    collection = LCIO.getCollection(event, collectionname)
                    for c in collection
                        positions = nil(Vector{Float64})
                        p4vec = LCIO.getP4(c)
                        maxEnergy = p4vec.t
                        momentum = [p4vec.x, p4vec.y, p4vec.z]
                        particletype = getType(c)
                        for hit in getParticleHits(c)
                            positions = cons([hit.x, hit.y, hit.z, hit.E], positions)
                        end
                        f(hcat(positions...), i, maxEnergy, momentum, particletype)
                    end
                end
            end
        end
    end
end

"""
  maptruthparticles(f, filename, startevent, endevent)

Applies the function `f(positionvec4, i, maxenergy, momentum, particletype)` to the matrix of position-energy 4-vectors,
event number `i`, reconstructed energy `maxenergy`, momentum, and particle type for all particles in events `startevent:endevent`.
Hits are assigned to their true particles, with energies scaled to account for shared hits.
`positionvec4` is a 4xn matrix where n is the number of hits, and the last component is the energy.
"""
function maptruthparticles(f, filename, startevent, endevent)
    LCIO.open(filename) do event
        for (i, event) in zip(1:endevent, reader)
            if startevent < i continue end
            particleHits = Dict{Ptr{Void}, Dict}()
            ee = getCollection(i, "EcalEndcapHits")
            eb = getCollection(i, "EcalBarrelHits")
            he = getCollection(i, "HcalEndcapHits")
            hb = getCollection(i, "HcalBarrelHits")
            for relation in getCollection(i, "CalorimeterHitRelations")
            	calibratedHit = getRelationFrom(relation)
            	simCaloHit = getRelationTo(relation)
            	pList = getHitMCParticles(simCaloHit)
            	p4 = getCaloHit(calibratedHit)
            	energyList = getHitEnergyList(simCaloHit)
            	simCaloHitEnergy = sum(energyList)
            	# TODO map the particles to the Hits, but correct the hit energy by the sampling fraction
            	# maybe look at energy spread within the cluster
                if (simCaloHit in eb) || (simCaloHit in ee)
                    collectionName = "ECAL"
                elseif (simCaloHit in hb) || (simCaloHit in he)
                    collectionName = "HCAL"
                else
                    # skip muonhits
                    continue
                end
                for (kdx, particle) in enumerate(pList)
                    h = LCIO.CalHit(p4, energyList[kdx]*p4.E/simCaloHitEnergy)
                    if haskey(particleHits, particle)
                        push!(particleHits[particle][collectionName], h)
                    else
                        particleHits[particle] = Dict("ECAL" => LCIO.CalHit[], "HCAL" => LCIO.CalHit[])
                    end
                end
            end

            for (p, hitDictionary) in particleHits
                p4 = getMCP4(p)
                pdg = getMCPDGid(p)
                len1 = length(hitDictionary["ECAL"])
                len2 = length(hitDictionary["HCAL"])
                positions = Array(Float64, 4, len1+len2)
                ind = 1
                for hit in hitDictionary["ECAL"]
                    positions[:, ind] = [hit.x, hit.y, hit.z, hit.E]
                    ind += 1
                end
                for hit in hitDictionary["HCAL"]
                    positions[:, ind] = [hit.x, hit.y, hit.z, hit.E]
                    ind += 1
                end
                f(positions, idx, p4.t, [p4.x, p4.y, p4.z], pdg)
            end
        end
    end

end
end


text:
## Load the packages.
using LoopVectorization
using Statistics
using BenchmarkTools
using PyPlot
using Interpolations

## Create the refine function.
function refine_field_int!(hi, hi_tmp, lo, n_hi, n_lo)
    if n_hi ÷ n_lo == 2
        @tturbo for i in 0:2:size(hi, 1)-3
            i_lo = i÷2+2
            i_hi = i+2
            hi_tmp[i_hi  ] = lo[i_lo]
            hi_tmp[i_hi+1] = 1//2*lo[i_lo] + 1//2*lo[i_lo+1]
        end

        hi_tmp[1] = hi_tmp[end-1]
        hi_tmp[end] = hi_tmp[2]

        hi[:] = hi_tmp[:]

        # @tturbo for i in 2:size(hi, 1)-1
        #     hi[i] = 1//4*hi_tmp[i-1] + 1//2*hi_tmp[i] + 1//4*hi_tmp[i+1]
        # end
    elseif n_hi ÷ n_lo == 3
        @tturbo for i in 0:3:size(hi, 1)-3
            i_lo = i÷3+2
            i_hi = i+2
            hi_tmp[i_hi  ] = lo[i_lo  ]
            hi_tmp[i_hi+1] = lo[i_lo  ]
            hi_tmp[i_hi+2] = lo[i_lo+1]
        end

        hi_tmp[1] = hi_tmp[end-1]
        hi_tmp[end] = hi_tmp[2]

        @tturbo for i in 2:size(hi, 1)-1
            hi[i] = 1//3*hi_tmp[i-1] + 1//3*hi_tmp[i] + 1//3*hi_tmp[i+1]
        end
    else
        throw(DomainError(n_hi/n_lo, "Refinement should be 2 or 3"))
    end
end

## Set up the grids.
n_hi = 32
n_lo = n_hi ÷ 2

a_lo = rand(n_lo)
a_hi = zeros(n_hi)

a_hi_gc = zeros(n_hi+2)
a_lo_gc = zeros(n_lo+2)

a_hi_gc[2:n_hi+1] = a_hi[:]
a_lo_gc[2:n_lo+1] = a_lo[:]

a_hi_gc[1] = a_hi_gc[end-1]
a_hi_gc[end] = a_hi_gc[2]
a_lo_gc[1] = a_lo_gc[end-1]
a_lo_gc[end] = a_lo_gc[2]

dx_hi = 1/n_hi; dx_lo = 1/n_lo
x_hi = range(0, step=dx_hi, length=n_hi)
x_lo = range(0, step=dx_lo, length=n_lo)

## Compute.
a_hi_tmp = zeros(size(a_hi_gc))
@btime refine_field_int!(a_hi_gc, a_hi_tmp, a_lo_gc, n_hi, n_lo)
a_hi_int = a_hi_gc[2:end-1]

println(mean(a_lo) ≈ mean(a_hi_int))

## Plot the output.
figure()
plot(x_hi, a_hi_int, "C1-o")
plot(x_lo, a_lo, "k:+")
display(gcf())
show()


text:
_subscript(i::Integer) = i<0 ? error("$i is negative") : join('₀'+d for d in reverse(digits(i)))

function _parameters(x::RiskParameters{S}) where {S <: DiseaseStateSequence}
  inds = _indices(x, cumulative=false)
  params1 = ["ϵ", "Ωs", "κ", "Ωt", "Ωl", "Ωr"]
  params2 = String[]
  for k = eachindex(inds)
    for l = 1:inds[k]
      push!(params2, params1[k] * _subscript(l))
    end
  end
  return params2
end

function _parameters(x::N) where {N <: NucleicAcidSubstitutionModel}
  return ["$i" for i in propertynames(x)]
end

"""
Produces a data frame summarizing MCMC results.
"""
function Base.summary(
  mcmc::MCMC{S, M};
  burnin::Int64=0,
  thin::Int64=1,
  bychain::Bool=true,
  credibleinterval::Float64=0.95) where {
  S <: DiseaseStateSequence,
  M <: ILM}

  !(0.0 <= credibleinterval <= 1.0) && error("`credibleinterval` must be between 0 and 1")
  ci = (1-credibleinterval)/2
  params = _parameters(mcmc.markov_chains[1].risk_parameters[1])
  if M == PhyloILM
    append!(params, _parameters(mcmc.markov_chains[1].substitution_model[1]))
  end
  if bychain && length(mcmc.markov_chains) > 1
    df = DataFrame(parameter=String[],
                   chain=Int64[],
                   mean=Float64[],
                   var=Float64[],
                   CI=Tuple{Float64, Float64}[])
    for i = 1:length(mcmc.markov_chains)
      trace_array = convert(Array{Float64, 2}, mcmc.markov_chains[i].risk_parameters[1+burnin:thin:end])
      if M == PhyloILM
        trace_array = hcat(trace_array, convert(Array{Float64, 2}, mcmc.markov_chains[i].substitution_model[1+burnin:thin:end]))
      end
      append!(df,
        DataFrame(parameter = params,
                  chain = fill(i, size(trace_array, 2)),
                  mean = mean(trace_array, dims=1)[:],
                  var = var(trace_array, dims=1)[:],
                  CI = [Tuple(quantile(trace_array[:,j], [ci, 1-ci])) for j = 1:size(trace_array, 2)]))
    end
  else
    trace_array = convert(Array{Float64, 2}, mcmc.markov_chains[1].risk_parameters[1+burnin:thin:end])
    if M == PhyloILM
      trace_array = hcat(trace_array, convert(Array{Float64, 2}, mcmc.markov_chains[1].substitution_model[1+burnin:thin:end]))
    end
    for i = 2:length(mcmc.markov_chains)
      if M == PhyloILM
        vcat(trace_array,
          hcat(
            convert(Array{Float64, 2}, mcmc.markov_chains[i].risk_parameters[1+burnin:thin:end]),
            convert(Array{Float64, 2}, mcmc.markov_chains[i].substitution_model[1+burnin:thin:end])))
      else
        vcat(trace_array,
          convert(Array{Float64, 2}, mcmc.markov_chains[i].risk_parameters[1+burnin:thin:end]))
      end
    end
    df = DataFrame(parameter = params,
                   mean = mean(trace_array, dims=1)[:],
                   var = var(trace_array, dims=1)[:],
                   CI = [Tuple(quantile(trace_array[:,j], [ci, 1-ci])) for j = 1:size(trace_array, 2)])
  end
  return df
end


text:
# Test task copying

using Libtask
using Test

# Test case 1: stack allocated objects are deep copied.
function f_ct()
  t = 0;
  while true
    produce(t)
    t = 1 + t
  end
end

t = CTask(f_ct)

@test consume(t) == 0
@test consume(t) == 1
a = copy(t);
@test consume(a) == 2
@test consume(a) == 3
@test consume(t) == 2
@test consume(t) == 3

# Test case 2: heap allocated objects are shallowly copied.

function f_ct2()
  t = [0 1 2];
  while true
    produce(t[1])
    t[1] = 1 + t[1]
  end
end

t = CTask(f_ct2)

@test consume(t) == 0
@test consume(t) == 1
a = copy(t);
@test consume(a) == 2
@test consume(a) == 3
@test consume(t) == 4
@test consume(t) == 5
@test consume(a) == 6
@test consume(a) == 7

# Breaking test
function g_break()
    t = 0
    while true
        t[3] = 1
        produce(t)
        t = t + 1
    end
end

t = CTask(g_break)
@test_throws Libtask.CTaskException consume(t)


text:
"""
    benchmarkpkg(pkg, [target]::Union{String, BenchmarkConfig}; kwargs...)

Run a benchmark on the package `pkg` using the [`BenchmarkConfig`](@ref) or git identifier `target`.
Examples of git identifiers are commit shas, branch names, or e.g. `"HEAD~1"`.
Return a [`BenchmarkResults`](@ref).

The argument `pkg` can be a name of a package or a path to a directory to a package.

**Keyword arguments**:

* `script` - The script with the benchmarks, if not given, defaults to `benchmark/benchmarks.jl` in the package folder.
* `postprocess` - A function to post-process results. Will be passed the `BenchmarkGroup`, which it can modify, or return a new one.
* `resultfile` - If set, saves the output to `resultfile`
* `retune` - Force a re-tune, saving the new tuning to the tune file.
* `progressoptions` - Options (a `NamedTuple`) to be passed as keyword arguments to
  `ProgressMeter.Progress`.

The result can be used by functions such as [`judge`](@ref). If you choose to, you can save the results manually using
[`writeresults`](@ref) where `results` is the return value of this function. It can be read back with [`readresults`](@ref).

**Example invocations:**

```julia
using PkgBenchmark

import MyPkg
benchmarkpkg(pathof(MyPkg)) # run the benchmarks at the current state of the repository
benchmarkpkg(pathof(MyPkg), "my-feature") # run the benchmarks for a particular branch/commit/tag
benchmarkpkg(pathof(MyPkg), "my-feature"; script="/home/me/mycustombenchmark.jl")
benchmarkpkg(pathof(MyPkg), BenchmarkConfig(id = "my-feature",
                                            env = Dict("JULIA_NUM_THREADS" => 4),
                                            juliacmd = `julia -O3`))
benchmarkpkg(pathof(MyPkg),  # Run the benchmarks and divide the (median of) results by 1000
    postprocess=(results)->(results["g"] = median(results["g"])/1_000)
```
"""
function benchmarkpkg(
        pkg::String,
        target=BenchmarkConfig();
        script=nothing,
        postprocess=nothing,
        resultfile=nothing,
        retune=false,
        progressoptions=NamedTuple(),
        custom_loadpath="" #= used in tests =#
    )
    target = BenchmarkConfig(target)

    pkgid = Base.identify_package(pkg)
    pkgfile_from_pkgname = pkgid === nothing ? nothing : Base.locate_package(pkgid)

    if pkgfile_from_pkgname===nothing
        if isdir(pkg)
            pkgdir = pkg
        else
            error("No package '$pkg' found.")
        end
    else
        pkgdir = normpath(joinpath(dirname(pkgfile_from_pkgname), ".."))
    end

    # Locate script
    if script === nothing
        script = joinpath(pkgdir, "benchmark", "benchmarks.jl")
    elseif !isabspath(script)
        script = joinpath(pkgdir, script)
    end

    if !isfile(script)
        error("benchmark script at $script not found")
    end

    # Locate pacakge
    tunefile = joinpath(pkgdir, "benchmark", "tune.json")

    isgitrepo = ispath(joinpath(pkgdir, ".git"))
    if isgitrepo
        isdirty = LibGit2.with(LibGit2.isdirty, LibGit2.GitRepo(pkgdir))
        original_sha = _shastring(pkgdir, "HEAD")
    end

    # In this function the package is at the commit we want to benchmark
    function do_benchmark()
        shastring = begin
            if isgitrepo
                isdirty ? "dirty" : _shastring(pkgdir, "HEAD")
            else
                "non gitrepo"
            end
        end

        local results
        results_local = _withtemp(tempname()) do f
            _benchinfo("Running benchmarks...")
            _runbenchmark(script, f, target, tunefile;
                          retune = retune,
                          custom_loadpath = custom_loadpath,
                          progressoptions = progressoptions)
        end
        io = IOBuffer(results_local["results"])
        seek(io, 0)
        resgroup = BenchmarkTools.load(io)[1]
        if postprocess != nothing
            retval = postprocess(resgroup)
            if retval != nothing
                resgroup = retval
            end
        end
        juliasha = results_local["juliasha"]
        vinfo = results_local["vinfo"]
        results = BenchmarkResults(pkg, shastring, resgroup, now(), juliasha, vinfo, target)
        return results
    end

    if target.id !== nothing
        if !isgitrepo
            error("$pkgdir is not a git repo, cannot benchmark at $(target.id)")
        elseif isdirty
            error("$pkgdir is dirty. Please commit/stash your ",
                  "changes before benchmarking a specific commit")
        end
        results = _withcommit(do_benchmark, LibGit2.GitRepo(pkgdir), target.id)
    else
        results = do_benchmark()
    end

    if resultfile != nothing
        writeresults(resultfile, results)
        _benchinfo("benchmark results written to $resultfile")
    end
    if isgitrepo
        after_sha = _shastring(pkgdir, "HEAD")
        if original_sha != after_sha
            @warn("Failed to return back to original sha $original_sha, package now at $after_sha")
        end
    end
    return results
end

function _runbenchmark(file::String, output::String, benchmarkconfig::BenchmarkConfig, tunefile::String;
                      retune = false, custom_loadpath = nothing, progressoptions = NamedTuple())
    color = Base.have_color ? "--color=yes" : "--color=no"
    compilecache = "--compiled-modules=" * (Bool(Base.JLOptions().use_compiled_modules) ? "yes" : "no")
    _file, _output, _tunefile, _custom_loadpath = map(escape_string, (file, output, tunefile, custom_loadpath))
    codecov_option = Base.JLOptions().code_coverage
    coverage = if codecov_option == 0
        "none"
    elseif codecov_option == 1
        "user"
    else
        "all"
    end
    exec_str = isempty(_custom_loadpath) ? "" : "push!(LOAD_PATH, \"$(_custom_loadpath)\")\n"
    exec_str *=
        """
        using PkgBenchmark
        PkgBenchmark._runbenchmark_local($(repr(_file)), $(repr(_output)), $(repr(_tunefile)), $(repr(retune)), $(repr(progressoptions)))
        """

    target_env = [k => v for (k, v) in benchmarkconfig.env]
    withenv(target_env...) do
        env_to_use = dirname(Pkg.Types.Context().env.project_file)
        run(`$(benchmarkconfig.juliacmd) --project=$env_to_use --depwarn=no --code-coverage=$coverage $color $compilecache -e $exec_str`)
    end
    return JSON.parsefile(output)
end

function _runbenchmark_local(file, output, tunefile, retune, progressoptions)
    # Loading
    Base.include(Main, file)
    if !isdefined(Main, :SUITE)
        error("`SUITE` variable not found, make sure the BenchmarkGroup is named `SUITE`")
    end
    suite = Main.SUITE

    # Tuning
    if isfile(tunefile) && !retune
        _benchinfo("using benchmark tuning data in $(abspath(tunefile))")
        BenchmarkTools.loadparams!(suite, BenchmarkTools.load(tunefile)[1], :evals, :samples);
    else
        _benchinfo("creating benchmark tuning file $(abspath(tunefile))...")
        mkpath(dirname(tunefile))
        _tune!(suite, progressoptions = progressoptions)
        BenchmarkTools.save(tunefile, params(suite));
    end

    # Running
    results = _run(suite, progressoptions = progressoptions)

    # Output
    vinfo = first(split(sprint((io) -> versioninfo(io; verbose=true)), "Environment"))
    juliasha = Base.GIT_VERSION_INFO.commit

    open(output, "w") do iof
        JSON.print(iof, Dict(
            "results"  => sprint(BenchmarkTools.save, results),
            "vinfo"    => vinfo,
            "juliasha" => juliasha,
        ))
    end
    return nothing
end

function _tune!(group::BenchmarkTools.BenchmarkGroup; verbose::Bool = false, root = true,
                progressoptions = NamedTuple(),
                prog = Progress(length(BenchmarkTools.leaves(group)); desc = "Tuning: ", progressoptions...),
                hierarchy = [], kwargs...)
    BenchmarkTools.gcscrub() # run GC before running group, even if individual benchmarks don't manually GC
    i = 1
    for id in keys(group)
        _tune!(group[id]; verbose = verbose, prog = prog, hierarchy = push!(copy(hierarchy), (repr(id), i, length(keys(group)))), kwargs...)
        i += 1
    end
    return group
end

function _tune!(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters = b.params;
               prog = nothing, verbose::Bool = false, pad = "", hierarchy = [], kwargs...)
    BenchmarkTools.warmup(b, verbose=false)
    estimate = ceil(Int, minimum(BenchmarkTools.lineartrial(b, p; kwargs...)))
    b.params.evals = BenchmarkTools.guessevals(estimate)
    if prog != nothing
        indent = 0
        ProgressMeter.next!(prog; showvalues = [map(id -> ("  "^(indent += 1) * "[$(id[2])/$(id[3])]", id[1]), hierarchy)...])
    end
    return b
end

function _run(group::BenchmarkTools.BenchmarkGroup, args...;
              progressoptions = NamedTuple(),
              prog = Progress(length(BenchmarkTools.leaves(group)); desc = "Benchmarking: ", progressoptions...), hierarchy = [], kwargs...)
    result = similar(group)
    BenchmarkTools.gcscrub() # run GC before running group, even if individual benchmarks don't manually GC
    i = 1
    for id in keys(group)
        result[id] = _run(group[id], args...; prog = prog, hierarchy = push!(copy(hierarchy), (repr(id), i, length(keys(group)))), kwargs...)
        i += 1
    end
    return result
end

function _run(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters = b.params;
                   prog = nothing, verbose::Bool = false, pad = "", hierarchy = [], kwargs...)
    res = BenchmarkTools.run_result(b, p; kwargs...)[1]
    if prog != nothing
        indent = 0
        ProgressMeter.next!(prog; showvalues = [map(id -> ("  "^(indent += 1) * "[$(id[2])/$(id[3])]", id[1]), hierarchy)...])
    end
    return res
end


text:
#__precompile__()
module FTD2XX

@linux? begin const d2xx = "/usr/local/lib/libftd2xx.so" end: nothing
@osx? begin const d2xx = "/usr/local/lib/libftd2xx.dylib" end: nothing
@windows? begin const d2xx = "ftd2xx.dll" end: nothing

include("checkstatus.jl")
include("common_api.jl")
@linux? include("osx_linux_api.jl"):nothing
@osx? include("osx_linux_api.jl"):nothing
@windows? include("windows_api.jl"):nothing
include("IOftuart.jl")  # create a subtype of IO for the UART
end # module


text:
#Written By Johnathan Bizzano
export ReadArrayToDataSet

function ReadArrayToDataSet(inputs::AbstractArray{I}, outputs::AbstractArray{O}) where {I, O}
    data = MemoryDataSet{I, O}(size(inputs, 1))
    inputs = [row for row in eachrow(inputs)]
    outputs = [row for row in eachrow(outputs)]
    for i in eachindex(inputs)
        data[i] = DataEntry{I, O}(Array(inputs[i][1:end]), Array(outputs[i][1:end]))
    end
    return data
end


text:
## SymbolicRandVar
## ===============
"Return a Set of dimension indices of a random variable"
function dims(X::SymbolicRandVar)
  # Do a depth first search and find union of dimensions of all OmegaRandVars
  dimensions = Set{Int}()
  visited = Set{RandVar}()
  tovisit = Set{RandVar}([X])
  while !isempty(tovisit)
    v = pop!(tovisit)
    if has_single_dim(v)
      push!(dimensions, v.dim)
    end
    for arg in args(v)
      arg ∉ visited && push!(tovisit,arg)
    end
  end
  dimensions
end

# "Apply a random variable to some randomness"
# call(X::SymbolicRandVar,ω) = lambda(X)(ω)

# Base case is that symbolic rand vars have multiple dims
has_single_dim(X::SymbolicRandVar) = false

function isequal(X::SymbolicRandVar, Y::SymbolicRandVar)
  # Equivalent Random variables should (at least) have same type and #args
  typeof(X) != typeof(Y) && (return false)
  x_args = fields(X)
  y_args = fields(Y)
  length(x_args) != length(y_args) && (return false)
  for i = 1:length(x_args)
    !isequal(x_args[i],y_args[i]) && (return false)
  end
  true
end

in{T}(X::SymbolicRandVar, bounds::Tuple{Lift{T},Lift{T}}) = (X >= bounds[1]) & (X <=  bounds[2])

## Constant Random Variable
## ========================
"A constant value. A constant function which 'ignores' input, e.g. ω->5"
immutable ConstantRandVar{T} <: SymbolicRandVar{T}
  val::T
end

args(X::ConstantRandVar) = tuple()
dims(X::ConstantRandVar) = Set{Int}()
# QUESTION: Would this ever be called? Do we need to overload functiosn for all constants?
(==)(X::ConstantRandVar, Y::ConstantRandVar) = ConstantRandVar{Bool}(X.val == Y.val)
isequal(X::ConstantRandVar, Y::ConstantRandVar) = isequal(X.val,Y.val)

# FIXME: Hack? Maybe the need for this will go away if RandArrays are of SymbolicRandVar
one{T}(::Type{RandVar{T}}) = ConstantRandVar(one(T))
zero{T}(::Type{RandVar{T}}) = ConstantRandVar(zero(T))
zero{T}(::RandVar{T}) = ConstantRandVar(zero(T))
norm(X::ConstantRandVar) = ConstantRandVar(norm(X.val))

one{T}(::Type{SymbolicRandVar{T}}) = ConstantRandVar(one(T))
zero{T}(::Type{SymbolicRandVar{T}}) = ConstantRandVar(zero(T))

real{T}(X::SymbolicRandVar{T}) = X

## Omega Random Variable
## =====================
"Simplest RandVar: ω->ω[dim] - extracts dim component of omega"
immutable OmegaRandVar{T} <: SymbolicRandVar{T}
  dim::Int
end

args(X::OmegaRandVar) = Set{RandVar}()
dims(X::OmegaRandVar) = Set(X.dim)
omega_component{T<:Real}(i,OmegaType::Type{T}=Float64) = OmegaRandVar{OmegaType}(i)
omega_component{T<:Real}(OmegaType::Type{T}=Float64) = OmegaRandVar{OmegaType}(genint())
isequal(X::OmegaRandVar,Y::OmegaRandVar) = isequal(X.dim,Y.dim)
has_single_dim(X::OmegaRandVar) = true

## Real × Real -> Real
## ===================
real_real_real = ((:PlusRandVar,:+),(:MinusRandVar,:-),(:TimesRandVar,:*),(:DivideRandVar,:/), (:PowRandVar,:(^)))
for (name,op) in real_real_real
  eval(
  quote
  immutable $name{T<:Real,A1<:Real,A2<:Real} <: SymbolicRandVar{T}
    args::Tuple{SymbolicRandVar{A1},SymbolicRandVar{A2}}
  end
  # (^) Fixes ambiguities. Redefined here in each loop iteration but shouldn't matter
  # (^){T1<:Real,T2<:Integer}(X::SymbolicRandVar{T1},c::T2) = PowRandVar{promote_type(T1, T2),T1,T2}((X,ConstantRandVar(c)))
  ($op){T1<:Real, T2<:Real}(X::SymbolicRandVar{T1}, Y::SymbolicRandVar{T2}) = $name{promote_type(T1, T2),T1,T2}((X,Y))
  ($op){T1<:Real, T2<:Real}(X::SymbolicRandVar{T1}, c::T2) = $name{promote_type(T1, T2),T1,T2}((X,ConstantRandVar(c)))
  ($op){T1<:Real, T2<:Real}(c::T1, X::SymbolicRandVar{T2}) = $name{promote_type(T1, T2),T1,T2}((ConstantRandVar(c),X))
  end)
end

# Real -> Real
## ===========
real_real = ((:UnaryPlusRandVar,:+),(:UnaryMinusRandVar,:-),(:AbsRandVar,:abs))
for (name,op) in real_real
  eval(
  quote
  immutable $name{T<:Real, A1<:Real} <: SymbolicRandVar{T}
    args::Tuple{SymbolicRandVar{A1}}
  end
  ($op){T<:Real}(X::SymbolicRandVar{T}) = $name{T,T}((X,))
  end)
end

# Real -> _<:Floating
## ==================
real_floating = (
  (:SqrtRandVar, :sqrt), (:ExpRandVar,:exp), (:LogRandVar,:log), (:SinRandVar,:sin),
  (:CosRandVar,:cos), (:TanRandVar,:tan), (:AsinRandVar,:asin),
  (:AcosRandVar,:acos), (:AtanRandVar,:atan), (:SinhRandVar,:sinh),
  (:CoshRandVar,:cosh), (:TanhRandVar,:tanh), (:Atan2RandVar,:atan2))

for (name,op) in real_floating
  eval(
  quote
  immutable $name{T<:Real,A1<:Real} <: SymbolicRandVar{T}
    args::Tuple{SymbolicRandVar{A1}}
  end
  ($op){T<:Real}(X::SymbolicRandVar{T}, returntype::DataType = Float64) = $name{returntype,T}((X,))
  end)
end

# Real × Real -> Bool
## ===================
real_real_bool = ((:GTRandVar, :>), (:GTERandVar,:>=), (:LTERandVar,:<=), (:LTRandVar,:<),
                  (:EqRandVar, :(==)), (:NeqRandVar, :!=))

for (name,op) in real_real_bool
  eval(
  quote
  immutable $name{T<:Real,A1<:Real,A2<:Real} <: SymbolicRandVar{T}
    args::Tuple{SymbolicRandVar{A1},SymbolicRandVar{A2}}
  end
  ($op){T1<:Real, T2<:Real}(X::SymbolicRandVar{T1}, Y::SymbolicRandVar{T2}) = $name{Bool,T1,T2}((X,Y))
  ($op){T1<:Real, T2<:Real}(X::SymbolicRandVar{T1}, c::T2) = $name{Bool,T1,T2}((X,ConstantRandVar(c)))
  ($op){T1<:Real, T2<:Real}(c::T1, X::SymbolicRandVar{T2}) = $name{Bool,T1,T2}((ConstantRandVar(c),X))
  end)
end

## Real × Real -> Bool
## ===================
bool_bool_bool = ((:OrRandVar, :|), (:AndRandVar,:&), (:BicondRandVar, :(==)))
for (name,op) in bool_bool_bool
  eval(
  quote
  immutable $name{T,A1,A2} <: SymbolicRandVar{Bool}
    args::Tuple{SymbolicRandVar{A1},SymbolicRandVar{A2}}
  end
  ($op)(X::SymbolicRandVar{Bool}, Y::SymbolicRandVar{Bool}) = $name{Bool,Bool,Bool}((X,Y))
  ($op)(X::SymbolicRandVar{Bool}, c::Bool) = $name{Bool,Bool,Bool}((X,ConstantRandVar(c)))
  ($op)(c::Bool, X::SymbolicRandVar{Bool}) = $name{Bool,Bool,Bool}((ConstantRandVar(c),X))
  end)
end

## Bool -> Bool
## ============
immutable NotRandVar{T,A1} <: SymbolicRandVar{Bool}
  args::Tuple{SymbolicRandVar{A1}}
end
!(X::SymbolicRandVar{Bool})= NotRandVar{Bool,Bool}((X,))

immutable IfElseRandVar{T,A1,A2,A3} <: SymbolicRandVar{T}
  args::Tuple{SymbolicRandVar{A1},SymbolicRandVar{A2},SymbolicRandVar{A3}}
end

## ifelse
## ======
ifelse{T}(A::SymbolicRandVar{Bool}, B::SymbolicRandVar{T}, C::SymbolicRandVar{T}) = IfElseRandVar{T,Bool,T,T}((A,B,C))
ifelse{T<:Real}(A::SymbolicRandVar{Bool}, B::T, C::T) = IfElseRandVar{T,Bool,T,T}((A,ConstantRandVar(B),ConstantRandVar(C)))
ifelse{T<:Real}(A::SymbolicRandVar{Bool}, B::SymbolicRandVar{T}, C::T) = IfElseRandVar{T,Bool,T,T}((A,B,ConstantRandVar(C)))
ifelse{T<:Real}(A::SymbolicRandVar{Bool}, B::T, C::SymbolicRandVar{T}) = IfElseRandVar{T,Bool,T,T}((A,ConstantRandVar(B),C))

# Unions
## =====
BinaryRealExpr = Union{PlusRandVar, MinusRandVar, TimesRandVar,DivideRandVar,PowRandVar}
UnaryRealExpr = Union{UnaryPlusRandVar,UnaryMinusRandVar,AbsRandVar}
TrigExpr = Union{ExpRandVar,LogRandVar,SinRandVar,CosRandVar,TanRandVar,AsinRandVar,
                 AcosRandVar,AtanRandVar,SinhRandVar,CoshRandVar,TanhRandVar,Atan2RandVar}
IneqExpr = Union{GTRandVar,GTERandVar,LTERandVar, LTRandVar,EqRandVar,NeqRandVar}
LogicalExpr = Union{OrRandVar,AndRandVar, BicondRandVar, NotRandVar}

# All Functional expressions
CompositeRandVar = Union{BinaryRealExpr, UnaryRealExpr, TrigExpr, IneqExpr, LogicalExpr, SqrtRandVar, IfElseRandVar}

args{T<:CompositeRandVar}(X::T) = X.args

## Printing
## ========
function to_dimacs(Y::SymbolicRandVar{Bool})
  cmap, cnf, ω, aux_vars = analyze(Y)
  indep_vars = join(values(cmap), " ")
  last_var = var(cnf[length(cnf)][1])
  println("cnf $last_var $(length(cnf))")
  println(string(cnf))
  println("c ind $indep_vars")
end


text:
export perturb

# Markov derivatives
@compat struct ExpandingBranchDerivative{B<:ExpandingBranch}
  b::B
end
(bd::ExpandingBranchDerivative)(x::Number) = mapD(bd.b,x)
Base.adjoint(b::ExpandingBranch) = ExpandingBranchDerivative(b)
# Base.transpose(b::ExpandingBranch) = ExpandingBranchDerivative(b)

@compat struct ExpandingBranchInverse{B<:ExpandingBranch}
  b::B
end
(bi::ExpandingBranchInverse)(x::Number) = mapinv(bi.b,x)
inv(b::ExpandingBranch) = ExpandingBranchInverse(b)
## TODO one day: convert routine to ExpandingBranch?

@compat struct ExpandingBranchDerivativeInverse{B<:ExpandingBranch}
  b::B
end
(bi::ExpandingBranchDerivativeInverse)(x::Number) = mapinvD(bd.b,x)
Base.adjoint(bi::ExpandingBranchInverse) = ExpandingBranchDerivativeInverse(bi.b)
# Base.transpose(b::ExpandingBranch) = ExpandingBranchDerivative(b)

@compat struct MarkovMapDerivative{M<:AbstractMarkovMap}
  m::M
end
(md::MarkovMapDerivative)(x::Number) = mapD(md.m,x)
Base.adjoint(b::AbstractMarkovMap) = MarkovMapDerivative(b)

# Linear response perturbations
"""
    perturb(d, X, ϵ)

Construct a self-map on domain d: x ↦ x + ϵ X(x)
"""
perturb(d,X,ϵ) = perturb(convert(Domain,d),X,ϵ)
perturb(d::AbstractInterval,X,ϵ) = MarkovMap([x->x+ϵ*X(x)],[d],d)
perturb(d::PeriodicDomain,X,ϵ) = FwdCircleMap([x->x+ϵ*X(x)],d)

"""
    perturb(m::AbstractMarkovMap, X, ϵ)

Output perturbation of m: x ↦ m(x) + ϵ X(m(x))
"""
perturb(m::AbstractMarkovMap,X,ϵ) = perturb(rangedomain(m),X,ϵ)∘m

# Eigvals overloads
"""
    eigvals(m::AbstractMarkovMap, n)

Output eigenvalues of Transfer(m) using n×n Galerkin discretisation.

Calls directly to ApproxFun: you can also call eigvals(Transfer(m), n)
"""
@compat LinearAlgebra.eigvals(m::AbstractMarkovMap,n::Int64) = LinearAlgebra.eigvals(Transfer(m),n)
ApproxFun.eigs(m::AbstractMarkovMap,n::Int64) = ApproxFun.eigs(Transfer(m),n)

"""
    eigvecs(m::AbstractMarkovMap, n)

Output eigenfunctions of Transfer(m) using n×n Galerkin discretisation.

Calls directly to ApproxFun: you can also call eigvecs(Transfer(m), n)
"""
@compat LinearAlgebra.eigvecs(m::AbstractMarkovMap,n::Int64) = ApproxFun.eigvecs(Transfer(m),n)

# # plotting
# function plot(m::MarkovMap)
#   pts = cfstype(m)[]
#   vals = cfstype(m)[]
#   sp = sortperm([minimum(∂(domain(b))) for b in branches(m)])
#   for b in branches(m)[sp]
#     append!(pts,points(domain(b),100))
#     append!(vals,broadcast(b,points(domain(b),100)))
#   end
#   p = PyPlot.plot(pts,vals)
#   xlim(∂(domain(m)))
#   ylim(∂(rangedomain(m)))
#   xlabel("\$x\$")
#   ylabel("\$f(x)\$")
#   p
# end

# function plot(m::AbstractMarkovMap)
#   pts = points(domain(m),500)
#   vals = broadcast(m,pts)
#   p = PyPlot.plot(pts,vals)
#   xlim(∂(domain(m)))
#   ylim(∂(rangedomain(m)))
#   xlabel("\$x\$")
#   ylabel("\$f(x)\$")
#   p
# end


text:
#
# nctypes.jl -
#
# Definitions of types for the interface to the Nüvü Camēras SDK.
#

import Base.Libc: TmStruct

"""

All C functions to the Nüvü Camēras SDK yield a value of type `Status`.  Which
can be `NC_SUCCESS` or some other value which is an error code.

"""
struct Status; code::Cint; end

"""

`NuvuCameraError` is the type of exception thrown by calls to the C functions
of the Nüvü Camēras SDK in case of error.

"""
struct NuvuCameraError <: Exception
    func::Symbol
    status::Status
end

"""

`Name` is the union of valid types for the name of a parameter or a file.  In
practice, these are the types which can be automatically converted to a
`Cstring` by `ccall`.

"""
const Name = Union{AbstractString,Symbol}

"""

Abstract type `Handle` is the super type of all types used to store the
address of an opaque sructure.  All such types have a field `ptr` with the
pointer.

"""
abstract type Handle end

for T in (:Cam, :Grab, :Proc, :CtrlList, :CropModeSolutions,
          :ImageSaved, :StatsCtx)
    @eval begin

        struct $T <: Handle; ptr::Ptr{Void}; end

        @doc @doc(Handle) $T

    end
end

struct ImageParams{T<:Union{Cam,Grab}} <: Handle
    ptr::Ptr{Void}
end

# Pixel type for images (`unsigned short` in C code).
const Image = UInt16

const ImageBuffer{T} = Union{Ptr{T},DenseVector{T},DenseMatrix{T}}


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

mutable struct MD
    content::Vector{Any}
    meta::Dict{Symbol, Any}

    MD(content::AbstractVector, meta::Dict = Dict()) =
        new(content, meta)
end

MD(xs...) = MD(vcat(xs...))

function MD(cfg::Config, xs...)
    md = MD(xs...)
    md.meta[:config] = cfg
    return md
end

config(md::MD) = md.meta[:config]::Config

# Forward some array methods

Base.push!(md::MD, x) = push!(md.content, x)
Base.getindex(md::MD, args...) = md.content[args...]
Base.setindex!(md::MD, args...) = setindex!(md.content, args...)
Base.lastindex(md::MD) = lastindex(md.content)
Base.firstindex(md::MD) = firstindex(md.content)
Base.length(md::MD) = length(md.content)
Base.isempty(md::MD) = isempty(md.content)
Base.copy(md::MD) = MD(copy(md.content), copy(md.meta))

==(a::MD, b::MD) = (html(a) == html(b))

# Parser functions:
#   md – should be modified appropriately
#   return – basically, true if parse was successful
#     false uses the next parser in the queue, true
#     goes back to the beginning
#
# Inner parsers:
#   return – element to use or nothing

# Inner parsing

function parseinline(stream::IO, md::MD, parsers::Vector{Function})
    for parser in parsers
        inner = parser(stream, md)
        inner ≡ nothing || return inner
    end
end

function parseinline(stream::IO, md::MD, config::Config)
    content = []
    buffer = IOBuffer()
    while !eof(stream)
        char = peek(stream, Char)
        if haskey(config.inner, char) &&
                (inner = parseinline(stream, md, config.inner[char])) !== nothing
            c = String(take!(buffer))
            !isempty(c) && push!(content, c)
            buffer = IOBuffer()
            push!(content, inner)
        else
            write(buffer, read(stream, Char))
        end
    end
    c = String(take!(buffer))
    !isempty(c) && push!(content, c)
    return content
end

parseinline(s::AbstractString, md::MD, c::Config) =
    parseinline(IOBuffer(s), md, c)

parseinline(s, md::MD) = parseinline(s, md, config(md))

# Block parsing

function parse(stream::IO, block::MD, config::Config; breaking = false)
    skipblank(stream)
    eof(stream) && return false
    for parser in (breaking ? config.breaking : [config.breaking; config.regular])
        parser(stream, block) && return true
    end
    return false
end

parse(stream::IO, block::MD; breaking = false) =
  parse(stream, block, config(block), breaking = breaking)

function parse(stream::IO; flavor = julia)
    isa(flavor, Symbol) && (flavor = flavors[flavor])
    markdown = MD(flavor)
    while parse(stream, markdown, flavor) end
    return markdown
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
using ExprOptimization, ExprRules
using Test, Random

let
    grammar = @grammar begin
        R = R + R
        R = |(1:3)
    end

    function loss(node::RuleNode, grammar::Grammar)
        Core.eval(node, grammar)
    end

    Random.seed!(0)
    p = PIPE(PPT(0.8),20,2,0.2,0.1,0.05,1,0.2,0.6,0.999,10)
    res = optimize(p, grammar, :R, loss)
    @test res.expr == 1
    @test Core.eval(res.tree, grammar) == 1
    @test res.loss == 1

    iter = ExpressionIterator(grammar, 2, :R)
    pop = collect(iter)

    losses = Vector{Float64}(undef,length(pop))
    (best_tree, best_loss) = PIPEs.evaluate!(loss, grammar, pop, losses, pop[1], Inf)
    @test Core.eval(best_tree, grammar) == 1
    @test best_loss == 1

    ppt = PIPEs.PPTNode(p.ppt_params, grammar)
    PIPEs.update!(p, ppt, grammar, pop[1], losses[1], best_loss)
    PIPEs.p_target(p.ppt_params, ppt, grammar, pop[1], losses[1], best_loss, p.α, p.ϵ)
    PIPEs.mutate!(ppt, grammar, pop[1], p.p_mutation, p.β)
end


text:
let # preventing global scope
    # ===================================================================
    # Discretization of a continuous-time system without input (ZeroSet)
    # ===================================================================
    A = sparse([1, 1, 2, 3, 4], [1, 2, 2, 4, 3], [1., 2., 3., 4., 5.], 4, 4)
    X0 = BallInf(zeros(4), 0.1)
    cont_sys_homog = IVP(CLCS(A, Universe(4)), X0)
    δ = 0.01

    # no bloating, use default options
    discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="nobloating")
    @test inputdim(discr_sys_homog) == 0 # there is no input set!

    # no bloating, use Pade approximation
    discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="nobloating", exp_method="pade")
    @test inputdim(discr_sys_homog) == 0

    # bloating, default options
    discr_sys_homog = discretize(cont_sys_homog, δ)
    @test inputdim(discr_sys_homog) == 0

    # bloating, first order
    discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="firstorder")
    @test inputdim(discr_sys_homog) == 0

    # bloating, forward interpolation
    discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="forward")
    @test inputdim(discr_sys_homog) == 0

    # bloating, backward interpolation
    discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="backward")
    @test inputdim(discr_sys_homog) == 0

    # bloating, forward interpolation using lazy operations (default)
    discr_sys_homog = discretize(cont_sys_homog, δ, algorithm="forward", set_operations="lazy")
    @test discr_sys_homog.x0 isa ConvexHull
    @test inputdim(discr_sys_homog) == 0

    # using concrete zonotopes
    discr_sys_homog = discretize(cont_sys_homog, δ, set_operations="zonotope")
    @test discr_sys_homog.x0 isa Zonotope
    @test inputdim(discr_sys_homog) == 0

    # ===============================================================
    # Discretization of a continuous-time system with constant input
    # ===============================================================
    U = ConstantInput(Ball2(ones(4), 0.5))
    B = Matrix{Float64}(1.0I, 4, 4)
    cont_sys = IVP(CLCCS(A, B, Universe(4), U), X0)

    # no bloating
    discr_sys = discretize(cont_sys, δ, algorithm="nobloating")
    @test inputdim(discr_sys) == 4

    inputs = next_set(inputset(discr_sys))
    @test dim(inputs) == 4
    @test isa(inputs, LinearMap)
    @test isa(inputs.X, Ball2) && inputs.X.center == ones(4) && inputs.X.radius == 0.5

    # no bloating, use Pade approximation
    discr_sys = discretize(cont_sys, δ, algorithm="nobloating", exp_method="pade")
    @test inputdim(discr_sys) == 4

    # bloating, use scaling and squaring method
    discr_sys = discretize(cont_sys, δ, exp_method="base")
    @test inputdim(discr_sys) == 4

    inputs = next_set(inputset(discr_sys))
    @test dim(inputs) == 4
    @test isa(inputs, MinkowskiSum)

    # bloating, use Pade approximation
    discr_sys = discretize(cont_sys, δ, exp_method="pade")
    @test inputdim(discr_sys) == 4

    discr_sys = discretize(cont_sys, δ, algorithm="firstorder")
    @test inputdim(discr_sys) == 4

    # using concrete zonotopes
    # since a Ball2 cannot be converted to a Zonotope, we use a box instead
    U = ConstantInput(BallInf(ones(4), 0.5))
    cont_sys = IVP(CLCCS(A, B, Universe(4), U), X0)
    discr_sys = discretize(cont_sys, δ, set_operations="zonotope")
    @test discr_sys.x0 isa Zonotope
    @test inputdim(discr_sys) == 4

    # ===================================================================
    # Discretization of a continuous-time system with time-varying input
    # ===================================================================
    Ui = [Ball2(0.01*i*ones(4), i*0.2) for i in 1:3]
    cont_sys = IVP(CLCCS(A, I(4), Universe(4), VaryingInput(Ui)), X0)

    # no bloating
    discr_sys = discretize(cont_sys, δ, algorithm="nobloating")
    Ui_d = inputset(discr_sys)
    @test length(Ui_d) == 3
    for (i, inputs) in enumerate(discr_sys.s.U)
        if i == 1
            @test dim(inputs) == 4
            @test isa(inputs, LinearMap)
            @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*ones(4) && inputs.X.radius == 0.2
        elseif i == 2
            @test dim(inputs) == 4
            @test isa(inputs, LinearMap)
            @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*2*ones(4) && inputs.X.radius == 0.2*2
        else
            @test dim(inputs) == 4
            @test isa(inputs, LinearMap)
            @test isa(inputs.X, Ball2) && inputs.X.center == 0.01*3*ones(4) && inputs.X.radius == 0.2*3
        end
    end

    # no bloating, use Pade approximation
    discr_sys = discretize(cont_sys, δ, algorithm="nobloating", exp_method="pade")

    # bloating
    discr_sys = discretize(cont_sys, δ, exp_method="base")
    Ui_d = inputset(discr_sys)
    @test length(Ui_d) == 3
    for (i, inputs) in enumerate(discr_sys.s.U)
        if i == 1
            @test dim(inputs) == 4
            @test isa(inputs, MinkowskiSum)
            X = inputs.X.X
            @test isa(X, Ball2) && X.center == 0.01*ones(4) && X.radius == 0.2
        elseif i == 2
            @test dim(inputs) == 4
            @test isa(inputs, MinkowskiSum)
            X = inputs.X.X
            @test isa(X, Ball2) && X.center == 0.01*2*ones(4) && X.radius == 0.2*2
        else
            @test dim(inputs) == 4
            @test isa(inputs, MinkowskiSum)
            X = inputs.X.X
            @test isa(X, Ball2) && X.center == 0.01*3*ones(4) && X.radius == 0.2*3
        end
    end

    # bloating, use Pade approximation
    discr_sys = discretize(cont_sys, δ, exp_method="pade")

    # lazy symmetric interval hull
    discr_sys = discretize(cont_sys, δ, sih_method="lazy")

    # concrete symmetric interval hull
    discr_sys = discretize(cont_sys, δ, sih_method="concrete")

    # ================================================================
    # Discretization of a continuous-time system with interval matrix
    # ================================================================
    A_itv = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    X0 = BallInf(zeros(2), 0.1)
    sys_homog = IVP(CLCS(A_itv, Universe(2)), X0)
    B = IntervalMatrix([1.0..1.0 0.0..0.0; 0.0..0.0 1.0..1.0])
    U = ConstantInput(BallInf(ones(2), 0.5))
    sys_heterog = IVP(CLCCS(A_itv, B, Universe(2), U), X0)
    algorithm = "interval_matrix"

    for cont_sys in [sys_homog, sys_heterog]
        # default order
        discr_sys = discretize(cont_sys, δ; algorithm=algorithm)

        # specific order
        discr_sys = discretize(cont_sys, δ; algorithm=algorithm, order=4)
    end
end


text:
function place_invest(invest, client, invest_path)
    #
    ## Get Symbol Info
    buy_precision, sell_precision, min_qty, stepSize, tickSize = Pinguu.get_precision(client, invest["symbol"])

    #
    ## Get price/qty
    if haskey(invest, "price")
        current = invest["price"]
        ## ADD FORMAT CHECK
    else
        sleep(0.5)
        current = client.get_ticker(symbol=invest["symbol"])["lastPrice"];
    end
    current_val = parse(Float64, current)
    #
    ## Apply mode
    if invest["mode"] == "falling"
        current_val = current_val * (1 - parse(Float64, invest["fall_by"])/100)
        current = Pinguu.price2string(current_val, sell_precision, round_up=false)
    end
    qty = floor(parse(Float64, invest["amount"]) / current_val, digits=buy_precision)

    if qty * current_val <= min_qty
        min_val = min_qty + ceil(min_qty - qty * current_val, digits=buy_precision)
        invest["status"] = "Wrong QTY"
        invest["message"] = "Minimum von $min_qty nicht erreicht. Klicke auf Stornieren und versuche es mit einem höheren Kaufwert erneut.<br> Minimum: $min_val €";
        open(invest_path, "w") do f
            JSON.print(f, invest, 4)
        end
        return invest
    end
    @info("Price: " * current)
    @info("QTY: $qty")
    sleep(0.5)
    order = client.order_limit_buy(
        symbol=invest["symbol"],
        quantity=qty,
        price=current)

    if order !== nothing
        invest["order"] = order
        invest["status"] = "Order placed"
    else
        invest["status"] = "Placing order failed"
    end

    return invest
end
export place_invest


text:
function ismembern(v1,v2)
    flagvec=isa.(indexin(v1,v2), Number)
end


text:
# Driver script for all unit tests

using Distributed
@everywhere using Pkg

if abspath(PROGRAM_FILE) == @__FILE__
    # when run as a script, activate the test environment of this package so that test dependencies are available
    @everywhere Pkg.activate(dirname(@__FILE__))
    @everywhere Pkg.instantiate()
end

using Jive

runtests(@__DIR__, skip=["runtests.jl", "helper"])


text:
using TermInterface

function match(p::PatVar, x, mem)
    if isassigned(mem, p.idx)
        return x == mem[p.idx]
    end
    mem[p.idx] = x
    true
end

match(p::PatLiteral{T}, x::T, mem) where {T} = (p.val == x)
match(p::PatLiteral, x, mem) = false

function match(p::PatTypeAssertion, x::T, mem) where {T}
    if T <: p.type
        return match(p.var, x, mem)
    end
    false
end

# TODO PatSplatVar

match(p::PatEquiv, x, mem) = error("PatEquiv can only be used in EGraphs rewriting")

function match(p::PatTerm, x, mem)
    !istree(typeof(x)) && (return false)
    if exprhead(p) == exprhead(x) && operation(p) == operation(x) && arity(p) == arity(x)
        p_args, x_args = arguments(p), arguments(x)
        for i in 1:arity(p)
            !match(p_args[i], x_args[i], mem) && (return false)
        end
        return true
    end
    false
end

function (r::RewriteRule)(x)
    mem = Vector(undef, length(r.patvars))
    if match(r.left, x, mem)
        return instantiate(x, r.right, mem)
    end
    return nothing
end

function (r::EqualityRule)(x)
    mem = Vector(undef, length(r.patvars))
    if match(r.left, x, mem)
        return instantiate(x, r.right, mem)
    end
    return nothing
end

function (r::DynamicRule)(x)
    # print("matching ")
    # display(r)
    # println(" against $x")
    mem = Vector(undef, length(r.patvars))
    if match(r.left, x, mem)
        # println("matched")
        return r.rhs_fun(x, mem, nothing, collect(mem)...)
    end
    # println("failed")
    return nothing
end

# TODO revise
function instantiate(left, pat::PatTerm, mem)
    ar = arguments(pat)
    similarterm(typeof(left), operation(pat),
        [instantiate(left, ar[i], mem) for i in 1:length(ar)]; exprhead=exprhead(pat))
end

instantiate(left, pat::PatLiteral, mem) = pat.val

function instantiate(left, pat::PatVar, mem)
    # println(left)
    # println(pat)
    # println(mem)
    mem[pat.idx]
end


text:
const samplerate = 44100
const windowlength = 256
const windowstep = 10
const lowharmonic = 10
const highharmonic = 64
const envelope_bandwidth = 10
const lowCutOff = 3 # When making sonograms, discard anything that is this many decibels below the loudest harmonic.


text:
include("../base.jl")
include("algorithms.jl")
include("fitnessFunctions.jl")

function f(x, y)
    (y/5.5)+x^2
end

function g(x, y)
    2*x + y
end

x = collect(-5:0.5:5)
y = collect(0:0.5:10)
nvalues = length(x)
objs1 = Array{Number}(undef, nvalues)
for i=1:nvalues
    objs1[i] = f(x[i], y[i])
end

nvalues = length(x)
objs2 = Array{Number}(undef, nvalues)
for i=1:nvalues
    objs2[i] = g(x[i], y[i])
end

mirng = Random.MersenneTwister(5)

@time experimentos=runGenJ("ECJ/tests/example.json", verbose=false)

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using Compat

function gdbm_open(name::String, flag::String="r")
  handle = ccall((:gdbm_open, libgdbm), Ptr{Cvoid}, (Cstring, Int32, Int32, Int32, Ptr{Cvoid}), name, 0, OPENFLAGS[flag], 420, C_NULL)
  handle == C_NULL && error("File could not be opened!")
  handle
end

function gdbm_close(handle::Ptr{Cvoid})
  ccall((:gdbm_close, libgdbm), Cvoid, (Ptr{Cvoid},), handle)
end

function gdbm_store(handle::Ptr{Cvoid}, key::Datum, value::Datum, flag::String="r")
  ret = ccall((:gdbm_store, libgdbm), Int32, (Ptr{Cvoid}, Datum, Datum, Int32), handle, key, value, STOREFLAGS[flag])
  ret == -1 && error("Database is not writable or key/value is not a valid string.")
  ret == 1 && error("Key is already in database.")
  nothing
end

function gdbm_fetch(handle::Ptr{Cvoid}, key::Datum)
  datum = ccall((:gdbm_fetch, libgdbm), Datum, (Ptr{Cvoid}, Datum), handle, key)
  datum.dptr == C_NULL && throw(KeyError(key))
  datum
end

function gdbm_exists(handle::Ptr{Cvoid}, key::Datum)
  ret = ccall((:gdbm_exists, libgdbm), Int32, (Ptr{Cvoid}, Datum), handle, key)
  ret == 0 && return false
  true
end

function gdbm_count(handle::Ptr{Cvoid})
  count = Ref(UInt(0))
  ret = ccall((:gdbm_count, libgdbm), Int32, (Ptr{Cvoid}, Ref{UInt}), handle, count)
  ret == -1 && error("Error reading database.")
  Int(count[])
end

function gdbm_delete(handle::Ptr{Cvoid}, key::Datum)
  ret = ccall((:gdbm_delete, libgdbm), Int32, (Ptr{Cvoid}, Datum), handle, key)
  ret ≠ 0 && error("Database is not writable or key not found.")
  nothing
end

function gdbm_firstkey(handle::Ptr{Cvoid})
  ccall((:gdbm_firstkey, libgdbm), Datum, (Ptr{Cvoid}, ), handle)
end

function gdbm_nextkey(handle::Ptr{Cvoid}, prev::Datum)
  ccall((:gdbm_nextkey, libgdbm), Datum, (Ptr{Cvoid}, Datum), handle, prev)
end


text:
# demo
# ====

# fibonacci
# ---------

fib(n) = n ≤ 2 ? n : fib(n-1) + fib(n-2)

fib(1000)   # never terminates in ordinal execution
fib(m)      # undef var
fib("1000") # obvious type error

# language features
# -----------------

# user-defined types
struct Ty{T}
    fld::T
end

function foo(a)
    v = Ty(a)
    return bar(v)
end

# macros will be expanded
@inline bar(n::T)     where {T<:Number} = n < 0 ? zero(T) : one(T)
@inline bar(v::Ty{T}) where {T<:Number} = bar(v.fdl) # typo "fdl"
@inline bar(v::Ty)                      = bar(convert(Number, v.fld))

foo(1.2)
foo("1") # `String` can't be converted to `Number`


text:
# Single-network containing only fixed data (i.e. data that does not depend on year), method without candidates
function nw(source::AbstractDict, lookup::AbstractDict, y::Int)
    target = Dict{String,Any}(
        "branch"       => Dict{String,Any}(),
        "branchdc"     => Dict{String,Any}(),
        "bus"          => Dict{String,Any}(),
        "busdc"        => Dict{String,Any}(),
        "convdc"       => Dict{String,Any}(),
        "dcline"       => Dict{String,Any}(),
        "gen"          => Dict{String,Any}(),
        "load"         => Dict{String,Any}(),
        "shunt"        => Dict{String,Any}(),
        "storage"      => Dict{String,Any}(),
        "switch"       => Dict{String,Any}(),
        "dcpol"        => 2, # Assumption: DC grid has 2 poles.
        "per_unit"     => true,
        "time_elapsed" => 1.0, # Assumption: each period lasts 1 hour.
    )
    if haskey(source["genericParameters"], "basePower")
        target["baseMVA"] = source["genericParameters"]["basePower"]
    else
        Memento.error(_LOGGER, "\"genericParameters.basePower\" is a required parameter.")
    end

    # AC branches are split between `acBranches` and `transformers` dicts in JSON files
    branch_path = ["gridModelInputFile", "acBranches"]
    for comp in walkpath(source, branch_path)
        index     = lookup["acBranches"][comp["id"]]
        source_id = push!(copy(branch_path), comp["id"])
        f_bus     = lookup["acBuses"][comp["acBusOrigin"]]
        t_bus     = lookup["acBuses"][comp["acBusExtremity"]]
        target["branch"]["$index"] = make_branch(comp, index, source_id, f_bus, t_bus, y; transformer=false)
    end
    branch_path = ["gridModelInputFile", "transformers"]
    for comp in walkpath(source, branch_path)
        index     = lookup["transformers"][comp["id"]]
        source_id = push!(copy(branch_path), comp["id"])
        f_bus     = lookup["acBuses"][comp["acBusOrigin"]]
        t_bus     = lookup["acBuses"][comp["acBusExtremity"]]
        target["branch"]["$index"] = make_branch(comp, index, source_id, f_bus, t_bus, y; transformer=true)
    end

    branchdc_path = ["gridModelInputFile", "dcBranches"]
    for comp in walkpath(source, branchdc_path)
        index     = lookup["dcBranches"][comp["id"]]
        source_id = push!(copy(branchdc_path), comp["id"])
        fbusdc    = lookup["dcBuses"][comp["dcBusOrigin"]]
        tbusdc    = lookup["dcBuses"][comp["dcBusExtremity"]]
        target["branchdc"]["$index"] = make_branchdc(comp, index, source_id, fbusdc, tbusdc, y)
    end

    bus_path = ["gridModelInputFile", "acBuses"]
    for comp in walkpath(source, bus_path)
        index     = lookup["acBuses"][comp["id"]]
        source_id = push!(copy(bus_path), comp["id"])
        target["bus"]["$index"] = make_bus(comp, index, source_id)
    end

    busdc_path = ["gridModelInputFile", "dcBuses"]
    for comp in walkpath(source, busdc_path)
        index     = lookup["dcBuses"][comp["id"]]
        source_id = push!(copy(busdc_path), comp["id"])
        target["busdc"]["$index"] = make_busdc(comp, index, source_id)
    end

    convdc_path = ["gridModelInputFile", "converters"]
    for comp in walkpath(source, convdc_path)
        index     = lookup["converters"][comp["id"]]
        source_id = push!(copy(convdc_path), comp["id"])
        busac     = lookup["acBuses"][comp["acBusConnected"]]
        busdc     = lookup["dcBuses"][comp["dcBusConnected"]]
        typeac    = walkpath(source, bus_path)[busac]["busType"]
        target["convdc"]["$index"] = make_convdc(comp, index, source_id, busac, busdc, typeac, y)
    end

    gen_path = ["gridModelInputFile", "generators"]
    for comp in walkpath(source, gen_path)
        index     = lookup["generators"][comp["id"]]
        source_id = push!(copy(gen_path), comp["id"])
        gen_bus   = lookup["acBuses"][comp["acBusConnected"]]
        target["gen"]["$index"] = make_gen(comp, index, source_id, gen_bus, y)
    end

    load_path = ["gridModelInputFile", "loads"]
    for comp in walkpath(source, load_path)
        index     = lookup["loads"][comp["id"]]
        source_id = push!(copy(load_path), comp["id"])
        load_bus  = lookup["acBuses"][comp["acBusConnected"]]
        target["load"]["$index"] = make_load(comp, index, source_id, load_bus, y)
    end

    storage_path = ["gridModelInputFile", "storage"]
    for comp in walkpath(source, storage_path)
        index       = lookup["storage"][comp["id"]]
        source_id   = push!(copy(storage_path), comp["id"])
        storage_bus = lookup["acBuses"][comp["acBusConnected"]]
        target["storage"]["$index"] = make_storage(comp, index, source_id, storage_bus, y)
    end

    return target
end

# Single-network containing only fixed data (i.e. data that does not depend on year), method with candidates
function nw(source::AbstractDict, lookup::AbstractDict, cand_availability::AbstractDict, y::Int)

    target = nw(source, lookup, y)
    target["branchdc_ne"] = Dict{String,Any}()
    target["busdc_ne"]    = Dict{String,Any}()
    target["convdc_ne"]   = Dict{String,Any}()
    target["ne_branch"]   = Dict{String,Any}()
    target["ne_storage"]  = Dict{String,Any}()

    bus_path = ["gridModelInputFile", "acBuses"] # Needed by convdc_ne

    ne_branch_path = ["candidatesInputFile", "acBranches"]
    for cand in walkpath(source, ne_branch_path)
        comp = cand["acBranch"]
        if cand_availability["acBranches"][comp["id"]][y]
            index     = lookup["cand_acBranches"][comp["id"]]
            source_id = push!(copy(ne_branch_path), comp["id"])
            f_bus     = lookup["acBuses"][comp["acBusOrigin"]]
            t_bus     = lookup["acBuses"][comp["acBusExtremity"]]
            t = make_branch(comp, index, source_id, f_bus, t_bus, y; transformer=false)
            t["construction_cost"] = cand["invCost"][y]
            t["lifetime"]          = cand["lifetime"]
            t["replace"]           = replace(cand, comp["id"], lookup["acBranches"]) # Assumption: specified id is that of the branch that connects the same buses.
            target["ne_branch"]["$index"] = t
        end
    end
    ne_branch_path = ["candidatesInputFile", "transformers"]
    for cand in walkpath(source, ne_branch_path)
        comp = cand["acBranch"]
        if cand_availability["transformers"][comp["id"]][y]
            index     = lookup["cand_transformers"][comp["id"]]
            source_id = push!(copy(ne_branch_path), comp["id"])
            f_bus     = lookup["acBuses"][comp["acBusOrigin"]]
            t_bus     = lookup["acBuses"][comp["acBusExtremity"]]
            t = make_branch(comp, index, source_id, f_bus, t_bus, y; transformer=true)
            t["construction_cost"] = cand["invCost"][y]
            t["lifetime"]          = cand["lifetime"]
            t["replace"]           = replace(cand, comp["id"], lookup["transformers"]) # Assumption: specified id is that of the branch that connects the same buses.
            target["ne_branch"]["$index"] = t
        end
    end

    branchdc_ne_path = ["candidatesInputFile", "dcBranches"]
    for cand in walkpath(source, branchdc_ne_path)
        comp = cand["dcBranch"]
        if cand_availability["dcBranches"][comp["id"]][y]
            index     = lookup["cand_dcBranches"][comp["id"]]
            source_id = push!(copy(branchdc_ne_path), comp["id"])
            fbusdc    = lookup["dcBuses"][comp["dcBusOrigin"]]
            tbusdc    = lookup["dcBuses"][comp["dcBusExtremity"]]
            t = make_branchdc(comp, index, source_id, fbusdc, tbusdc, y)
            t["cost"]     = cand["invCost"][y]
            t["lifetime"] = cand["lifetime"]
            target["branchdc_ne"]["$index"] = t
        end
    end

    convdc_ne_path = ["candidatesInputFile", "converters"]
    for cand in walkpath(source, convdc_ne_path)
        comp = cand["converter"]
        if cand_availability["converters"][comp["id"]][y]
            index     = lookup["cand_converters"][comp["id"]]
            source_id = push!(copy(convdc_ne_path), comp["id"])
            busac     = lookup["acBuses"][comp["acBusConnected"]]
            busdc     = lookup["dcBuses"][comp["dcBusConnected"]]
            typeac    = walkpath(source, bus_path)[busac]["busType"]
            t = make_convdc(comp, index, source_id, busac, busdc, typeac, y)
            t["cost"]     = cand["invCost"][y]
            t["lifetime"] = cand["lifetime"]
            target["convdc_ne"]["$index"] = t
        end
    end

    load_path = ["candidatesInputFile", "flexibleLoads"]
    for cand in walkpath(source, load_path)
        comp = cand["load"]
        if cand_availability["loads"][comp["id"]][y]
            index     = lookup["loads"][comp["id"]] # Candidate loads have same ids as existing loads
            source_id = push!(copy(load_path), comp["id"])
            load_bus  = lookup["acBuses"][comp["acBusConnected"]]
            t = make_load(comp, index, source_id, load_bus, y)
            t["cost_inv"] = cand["invCost"][y]
            t["lifetime"] = cand["lifetime"]
            target["load"]["$index"] = t # The candidate load replaces the existing load that has the same id and is assumed to have the same parameters as that existing load.
        end
    end

    ne_storage_path = ["candidatesInputFile", "storage"]
    for cand in walkpath(source, ne_storage_path)
        comp = cand["storage"]
        if cand_availability["storage"][comp["id"]][y]
            index     = lookup["cand_storage"][comp["id"]]
            source_id = push!(copy(ne_storage_path), comp["id"])
            storage_bus = lookup["acBuses"][comp["acBusConnected"]]
            t = make_storage(comp, index, source_id, storage_bus, y)
            t["eq_cost"]   = cand["invCost"][y]
            t["inst_cost"] = 0.0
            t["lifetime"]  = cand["lifetime"]
            target["ne_storage"]["$index"] = t
        end
    end

    return target
end

function walkpath(node::AbstractDict, path::Vector{String})
    for key in path
        node = node[key]
    end
    return node
end

function optional_value(target::AbstractDict, target_key::String, source::AbstractDict, source_key::String)
    if haskey(source, source_key)
        target[target_key] = source[source_key]
    end
end

function optional_value(target::AbstractDict, target_key::String, source::AbstractDict, source_key::String, y::Int)
    if haskey(source, source_key) && !isempty(source[source_key])
        target[target_key] = source[source_key][y]
    end
end

function replace(cand::AbstractDict, id::String, comp_lookup::AbstractDict)
    if haskey(cand, "replace")
        if haskey(comp_lookup, cand["replace"])
            # FlexPlan.jl only supports replacement of the only branch that connects the same buses as the candidate.
            # That branch will be replaced, regardless of which branch is specified by `cand["replace"]`.
            return true
        else
            Memento.warn(_LOGGER, "Cannot set \"$id\" to replace \"" * cand["replace"] * "\" because \"" * cand["replace"] * "\" does not exist.")
            return false
        end
    else
        return false
    end
end

function make_branch(source::AbstractDict, index::Int, source_id::Vector{String}, f_bus::Int, t_bus::Int, y::Int; transformer::Bool)
    target = Dict{String,Any}(
        "index"       => index,
        "source_id"   => source_id,
        "f_bus"       => f_bus,
        "t_bus"       => t_bus,
        "br_status"   => 1, # Assumption: all branches defined in JSON file are in service.
        "transformer" => transformer,
        "b_fr"        => 0.0, # Assumption: all branches defined in JSON file have zero shunt susceptance.
        "b_to"        => 0.0, # Assumption: all branches defined in JSON file have zero shunt susceptance.
        "g_fr"        => 0.0, # Assumption: all branches defined in JSON file have zero shunt conductance.
        "g_to"        => 0.0, # Assumption: all branches defined in JSON file have zero shunt conductance.
        "rate_a"      => source["ratedApparentPower"][y],
        "rate_c"      => source["emergencyRating"],
        "tap"         => source["voltageTapRatio"],
        "shift"       => 0.0, # Assumption: all branches defined in JSON file have zero shift.
        "angmin"      => source["minAngleDifference"],
        "angmax"      => source["maxAngleDifference"],
    )
    optional_value(target, "br_r", source, "resistance")
    if source["isTransmission"]
        target["br_r"] = 0.0
        target["br_x"] = 1/source["susceptance"]
    else
        target["br_r"] = source["resistance"]
        target["br_x"] = source["reactance"]
    end
    return target
end

function make_branchdc(source::AbstractDict, index::Int, source_id::Vector{String}, fbusdc::Int, tbusdc::Int, y::Int)
    target = Dict{String,Any}(
        "index"     => index,
        "source_id" => source_id,
        "fbusdc"    => fbusdc,
        "tbusdc"    => tbusdc,
        "status"    => 1, # Assumption: all branches defined in JSON file are in service.
        "rateA"     => source["ratedActivePower"][y],
        "rateC"     => source["emergencyRating"],
        "r"         => 0.0, # Assumption: zero resistance (the parameter is required by PowerModelsACDC but unused in lossless models).
    )
    return target
end

function make_bus(source::AbstractDict, index::Int, source_id::Vector{String})
    target = Dict{String,Any}(
        "index"     => index,
        "bus_i"     => index,
        "source_id" => source_id,
        "vm"        => source["nominalVoltageMagnitude"],
    )
    optional_value(target, "bus_type", source, "busType")
    if haskey(target, "bus_type") && target["bus_type"] == 3
        target["va"] = 0.0 # Set voltage angle of reference bus to 0.0
    end
    optional_value(target, "base_kv",  source, "baseVoltage")
    optional_value(target, "vmax",     source, "maxVoltageMagnitude")
    optional_value(target, "vmin",     source, "minVoltageMagnitude")
    if haskey(source, "location")
        target["lat"] = source["location"][1]
        target["lon"] = source["location"][2]
    end
    return target
end

function make_busdc(source::AbstractDict, index::Int, source_id::Vector{String})
    target = Dict{String,Any}(
        "index"     => index,
        "busdc_i"   => index,
        "source_id" => source_id,
        "Vdc"       => source["nominalVoltageMagnitude"],
        "Vdcmin"    => 0.9, # Assumption: minimum DC voltage is 0.9 p.u. for every DC bus
        "Vdcmax"    => 1.1, # Assumption: maximum DC voltage is 1.1 p.u. for every DC bus
        "Pdc"       => 0.0, # Assumption: power withdrawn from DC bus is 0.0 p.u.
    )
    optional_value(target, "basekVdc", source, "baseVoltage")
    return target
end

function make_convdc(source::AbstractDict, index::Int, source_id::Vector{String}, busac::Int, busdc::Int, typeac::Int, y::Int)
    target = Dict{String,Any}(
        "index"       => index,
        "source_id"   => source_id,
        "status"      => 1, # Assumption: all converters defined in JSON file are in service.
        "busac_i"     => busac,
        "busdc_i"     => busdc,
        "type_ac"     => typeac,
        "type_dc"     => 3, # Assumption: all converters defined in JSON file have DC droop.
        "Vmmin"       => 0.9, # Required by PowerModelsACDC, but not relevant, since we use an approximation where voltage magnitude is 1.0 p.u. at each AC transmission network bus
        "Vmmax"       => 1.1, # Required by PowerModelsACDC, but not relevant, since we use an approximation where voltage magnitude is 1.0 p.u. at each AC transmission network bus
        "Pacrated"    => source["ratedActivePowerAC"][y],
        "Pacmin"      => -source["ratedActivePowerAC"][y],
        "Pacmax"      => source["ratedActivePowerAC"][y],
        "Qacrated"    => 0.0, # Required by PowerModelsACDC, but unused in active power only models.
        "LossA"       => source["auxiliaryLosses"][y],
        "LossB"       => source["linearLosses"][y],
        "LossCinv"    => 0.0,
        "Imax"        => 0.0, # Required by PowerModelsACDC, but unused in lossless models.
        "transformer" => false, # Assumption: the converter is not a transformer.
        "tm"          => 0.0, # Required by PowerModelsACDC, but unused, provided that the converter is not a transformer.
        "rtf"         => 0.0, # Required by PowerModelsACDC, but unused, provided that the converter is not a transformer.
        "xtf"         => 0.0, # Required by PowerModelsACDC, but unused, provided that the converter is not a transformer.
        "reactor"     => false, # Assumption: the converter is not a reactor.
        "rc"          => 0.0, # Required by PowerModelsACDC, but unused, provided that the converter is not a reactor.
        "xc"          => 0.0, # Required by PowerModelsACDC, but unused, provided that the converter is not a reactor.
        "filter"      => false, # Required by PowerModelsACDC, but unused, provided that the model is lossless.
        "bf"          => 0.0, # Required by PowerModelsACDC, but unused, provided that the model is lossless.
        "islcc"       => 0.0, # Required by PowerModelsACDC, but unused, provided that the model is DC.
    )
    return target
end

function make_gen(source::AbstractDict, index::Int, source_id::Vector{String}, gen_bus::Int, y::Int)
    target = Dict{String,Any}(
        "index"        => index,
        "source_id"    => source_id,
        "gen_status"   => 1, # Assumption: all generators defined in JSON file are in service.
        "gen_bus"      => gen_bus,
        "qmin"         => source["minReactivePower"][y],
        "qmax"         => source["maxReactivePower"][y],
        "vg"           => 1.0,
        "model"        => 2, # Polynomial cost model
        "ncost"        => 2, # 2 cost coefficients: c1 and c0
        "cost"         => [source["generationCosts"][y], 0.0], # [c1, c0]
    )
    pmin = source["minActivePower"][y]
    pmax = source["maxActivePower"][y]
    target["pmax"] = pmax
    if pmin == pmax # Non-dispatchable generators are characterized in JSON file by having coincident power bounds
        target["dispatchable"] = false
        target["pmin"] = 0.0 # Must be zero to allow for curtailment
        target["cost_curt"] = source["curtailmentCosts"][y]
    else
        target["dispatchable"] = true
        target["pmin"] = pmin
    end
    return target
end

function make_load(source::AbstractDict, index::Int, source_id::Vector{String}, load_bus::Int, y::Int)
    target = Dict{String,Any}(
        "index"     => index,
        "source_id" => source_id,
        "status"    => 1, # Assumption: all loads defined in JSON file are in service.
        "load_bus"  => load_bus,
        "flex"      => get(source, "isFlexible", false),
    )
    if haskey(source, "powerFactor")
        target["pf_angle"] = acos(source["powerFactor"])
    end
    optional_value(target, "pred_rel_max",        source, "superiorBoundNCP",        y)
    optional_value(target, "ered_rel_max",        source, "maxEnergyNotConsumed",    y)
    optional_value(target, "pshift_up_rel_max",   source, "superiorBoundUDS",        y)
    optional_value(target, "pshift_down_rel_max", source, "superiorBoundDDS",        y)
    optional_value(target, "tshift_up",           source, "gracePeriodUDS",          y)
    optional_value(target, "tshift_down",         source, "gracePeriodDDS",          y)
    optional_value(target, "eshift_rel_max",      source, "maxEnergyShifted",        y)
    optional_value(target, "cost_curt",           source, "valueOfLossLoad",         y)
    optional_value(target, "cost_red",            source, "compensationConsumeLess", y)
    optional_value(target, "cost_shift",          source, "compensationDemandShift", y)
    return target
end

function make_storage(source::AbstractDict, index::Int, source_id::Vector{String}, storage_bus::Int, y::Int)
    target = Dict{String,Any}(
        "index"                 => index,
        "source_id"             => source_id,
        "status"                => 1, # Assumption: all generators defined in JSON file are in service.
        "storage_bus"           => storage_bus,
        "energy_rating"         => source["maxEnergy"][y],
        "charge_rating"         => source["maxAbsActivePower"][y],
        "discharge_rating"      => source["maxInjActivePower"][y],
        "charge_efficiency"     => source["absEfficiency"][y],
        "discharge_efficiency"  => source["injEfficiency"][y],
        "qmin"                  => source["minReactivePowerExchange"][y],
        "qmax"                  => source["maxReactivePowerExchange"][y],
        "self_discharge_rate"   => source["selfDischargeRate"][y],
        "r"                     => 0.0, # JSON API does not support `r`. Neither Flexplan.jl does (in lossless models), however a value is required by the constraint templates `_PM.constraint_storage_losses` and `_FP.constraint_storage_losses_ne`.
        "x"                     => 0.0, # JSON API does not support `x`. Neither Flexplan.jl does (in lossless models), however a value is required by the constraint templates `_PM.constraint_storage_losses` and `_FP.constraint_storage_losses_ne`.
        "p_loss"                => 0.0, # JSON API does not support `p_loss`. Neither Flexplan.jl does, however a value is required by the constraint templates `_PM.constraint_storage_losses` and `_FP.constraint_storage_losses_ne`.
        "q_loss"                => 0.0, # JSON API does not support `q_loss`. Neither Flexplan.jl does, however a value is required by the constraint templates `_PM.constraint_storage_losses` and `_FP.constraint_storage_losses_ne`.
    )
    # JSON API does not support storage thermal rating, but a value is required by
    # `FlexPlan.constraint_storage_thermal_limit`. The following expression prevents it from
    # limiting active or reactive power, even in the case of octagonal approximation of
    # apparent power.
    target["thermal_rating"] = 2 * max(target["charge_rating"], target["discharge_rating"], target["qmax"], -target["qmin"])
    optional_value(target, "max_energy_absorption", source, "maxEnergyYear", y)
    return target
end


text:
using Documenter

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

# pages="Phasing.md"  # for testing

makedocs(
    sitename = "Harmonie wiki",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = ["Home" => "index.md"; pages],
#   repo="https://hirlam.org/trac/wiki/HarmonieSystemDocumentation"
)

deploydocs(
    repo = "github.com/Hirlam/Trac2Markdown.jl.git",
)


text:
# Saner way to get modifier keys in tools
modifierControl() = keyHeld(Gtk.GdkKeySyms.Control_L) || keyHeld(Gtk.GdkKeySyms.Control_R)
modifierShift() = keyHeld(Gtk.GdkKeySyms.Shift_L) || keyHeld(Gtk.GdkKeySyms.Shift_R)
modifierMeta() = keyHeld(Gtk.GdkKeySyms.Meta_L) || keyHeld(Gtk.GdkKeySyms.Meta_R)
modifierAlt() = keyHeld(Gtk.GdkKeySyms.Alt_L) || keyHeld(Gtk.GdkKeySyms.Alt_R)
modifierSuper() = keyHeld(Gtk.GdkKeySyms.Super_L) || keyHeld(Gtk.GdkKeySyms.Super_R)
modifierHyper() = keyHeld(Gtk.GdkKeySyms.Hyper_L) || keyHeld(Gtk.GdkKeySyms.Hyper_R)

const modifierNames = Dict{String, Function}(
    "ctrl" => modifierControl,
    "shift" => modifierShift,
    "alt" => modifierAlt,
    "super" => modifierSuper,
    "win" => modifierSuper,
    "hyper" => modifierHyper,
    "meta" => modifierMeta
)

const specialNames = Dict{String, Int}(
    "plus" => Int('+'),
    "minus" => Int('-'),
)

struct Hotkey
    key::Integer
    callback::Function
    conditions::Array{Function, 1}

    Hotkey(key::Integer, callback::Function, modifiers::Array{Function, 1}) = new(key, callback, modifiers)
    Hotkey(key::Char, callback::Function, modifiers::Array{Function, 1}) = new(Int(key), callback, modifiers)
end

# Allow "keyword" charcters like Delete, Enter, etc
function Hotkey(s::String, callback::Function,  modifiers::Array{Function, 1}=Function[])
    key = -1

    words = strip.(split(s, "+"))

    for (i, word) in enumerate(words)
        if i < length(words)
            if haskey(modifierNames, word)
                push!(modifiers, modifierNames[word])

            else
                println(Base.stderr, "Invalid hotkey sequence")

                return nothing
            end
        end

        if i == length(words)
            if length(word) == 1
                key = Int(word[1])

            elseif haskey(specialNames, word)
                key = specialNames[word]

            # This is case sensitive!
            # Use values from Gtk.GdkKeySyms
            elseif isdefined(Gtk.GdkKeySyms, Symbol(word))
                key = getfield(Gtk.GdkKeySyms, Symbol(word))

            else
                println(Base.stderr, "Invalid hotkey sequence")

                return nothing
            end
        end
    end

    return Hotkey(key, callback, modifiers)
end

callback(h::Hotkey, args...) = h.callback(args...)

function active(h::Hotkey, event::Gtk.GdkEventKey)
    keyMeet = lowercase(Char(h.key)) == lowercase(Char(event.keyval))
    allModifiers = all((f -> f()).(h.conditions))

    return keyMeet && allModifiers
end

function callbackFirstActive(hotkeys::Array{Hotkey, 1}, event::Gtk.GdkEventKey)
    for hotkey in hotkeys
        if active(hotkey, event)
            return callback(hotkey)
        end
    end

    return false
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn