text:
# polygons, part of Luxor

"""
Draw a polygon.

    poly(pointlist::Array{Point, 1}, action = :none;
        close=false,
        reversepath=false)

A polygon is an Array of Points. By default `poly()` doesn't close or fill the polygon,
to allow for clipping.
"""
function poly(pointlist::Array{Point, 1}, action::Symbol = :none;
        close::Bool=false, reversepath::Bool=false)
    if action != :path
        newpath()
    end
    if reversepath == true
        reverse!(pointlist)
    end
    move(pointlist[1].x, pointlist[1].y)
    @inbounds for p in pointlist[2:end]
        line(p.x, p.y)
    end
    if close==true
        closepath()
    end
    do_action(action)
    return pointlist
end

"""
Find the centroid of simple polygon.

    polycentroid(pointlist)

Returns a point. This only works for simple (non-intersecting) polygons.

You could test the point using `isinside()`.
"""
function polycentroid(pointlist::Array{Point,1})
    # Points are immutable, use separate variables for these calculations
    centroid_x = 0.0
    centroid_y = 0.0
    signedarea = 0.0
    vertexcount = length(pointlist)
    x0 = 0.0 # Current vertex X
    y0 = 0.0 # Current vertex Y
    x1 = 0.0 # Next vertex X
    y1 = 0.0 # Next vertex Y
    a = 0.0  # Partial signed area

    # For all vertices except last
    i = 1
    @inbounds for i in 1:vertexcount-1
        x0 = pointlist[i].x
        y0 = pointlist[i].y
        x1 = pointlist[i+1].x
        y1 = pointlist[i+1].y
        a = x0 * y1 - x1 * y0
        signedarea += a
        centroid_x += (x0 + x1) * a
        centroid_y += (y0 + y1) * a
    end
    # Do last vertex separately to avoid performing an expensive
    # modulus operation in each iteration.
    x0 = pointlist[vertexcount].x
    y0 = pointlist[vertexcount].y
    x1 = pointlist[1].x
    y1 = pointlist[1].y
    a = x0 * y1 - x1 * y0
    signedarea += a
    centroid_x += (x0 + x1) * a
    centroid_y += (y0 + y1) * a
    signedarea *= 0.5
    centroid_x /= (6.0 * signedarea)
    centroid_y /= (6.0 * signedarea)
    return Point(centroid_x, centroid_y)
end

"""
Sort the points of a polygon into order. Points are sorted according to the angle they make
with a specified point.

    polysortbyangle(pointlist::Array, refpoint=minimum(pointlist))

The `refpoint` can be chosen, but the minimum point is usually OK too:

    polysortbyangle(parray, polycentroid(parray))
"""
function polysortbyangle(pointlist::Array{Point, 1}, refpoint=minimum(pointlist))
    angles = Float64[] ; sizehint!(angles, length(pointlist))
    @inbounds for pt in pointlist
        push!(angles, slope(pt, refpoint))
    end
    return pointlist[sortperm(angles)]
end

"""
Sort a polygon by finding the nearest point to the starting point, then
the nearest point to that, and so on.

    polysortbydistance(p, starting::Point)

You can end up with convex (self-intersecting) polygons, unfortunately.
"""
function polysortbydistance(pointlist::Array{Point, 1}, starting::Point)
    route = [starting]; sizehint!(route, length(pointlist))
    # start with the first point in pointlist
    remaining = setdiff(pointlist, route)
    while length(remaining) > 0
        # find the nearest point to the current position on the route
        nearest = first(sort!(remaining, lt = (x, y) -> distance(route[end], x) < distance(route[end], y)))
        # add this to the route and remove from remaining points
        push!(route, nearest)
        popfirst!(remaining)
    end
    return route
end

"""
Use a non-recursive Douglas-Peucker algorithm to simplify a polygon. Used by `simplify()`.

    douglas_peucker(pointlist::Array, start_index, last_index, epsilon)
"""
function douglas_peucker(pointlist::Array{Point, 1}, start_index, last_index, epsilon)
    temp_stack = Tuple{Int, Int}[]
    push!(temp_stack, (start_index, last_index))
    global_start_index = start_index
    keep_list = trues(length(pointlist))
    while length(temp_stack) > 0
        start_index = first(temp_stack[end])
        last_index =  last(temp_stack[end])
        pop!(temp_stack)
        dmax = 0.0
        index = start_index
        for i in index + 1:last_index - 1
            if (keep_list[i - global_start_index])
                d = pointlinedistance(pointlist[i], pointlist[start_index], pointlist[last_index])
                if d > dmax
                    index = i
                    dmax = d
                end
            end
        end
        if dmax > epsilon
            push!(temp_stack, (start_index, index))
            push!(temp_stack, (index, last_index))
        else
            for i in start_index + 2:last_index - 1 # 2 seems to keep the starting point...
                keep_list[i - global_start_index] = false
            end
        end
    end
    return pointlist[keep_list]
end

"""
Simplify a polygon:

    simplify(pointlist::Array, detail=0.1)

`detail` is the maximum approximation error of simplified polygon.
"""
function simplify(pointlist::Array{Point, 1}, detail=0.1)
    douglas_peucker(pointlist, 1, length(pointlist), detail)
end

function det3p(q1::Point, q2::Point, p::Point)
    (q1.x - p.x) * (q2.y - p.y) - (q2.x - p.x) * (q1.y - p.y)
end

"""
    isinside(p, pol; allowonedge=false)

Is a point `p` inside a polygon `pol`? Returns true if it does, or false.

This is an implementation of the Hormann-Agathos (2001) Point in Polygon algorithm.

The classification of points lying on the edges of the target polygon, or coincident with
its vertices is not clearly defined, due to rounding errors or arithmetical
inadequacy. By default these will generate errors, but you can suppress these by setting
`allowonedge` to `true`.
"""
function isinside(p::Point, pointlist::Array{Point, 1};
        allowonedge::Bool=false)
    c = false
    @inbounds for counter in eachindex(pointlist)
        q1 = pointlist[counter]
        # if reached last point, set "next point" to first point
        if counter == length(pointlist)
            q2 = pointlist[1]
        else
            q2 = pointlist[counter + 1]
        end
        if q1 == p
            allowonedge || error("isinside(): VertexException a")
            continue
        end
        if q2.y == p.y
            if q2.x == p.x
                allowonedge || error("isinside(): VertexException b")
                continue
            elseif (q1.y == p.y) && ((q2.x > p.x) == (q1.x < p.x))
                allowonedge || error("isinside(): EdgeException")
                continue
            end
        end
        if (q1.y < p.y) != (q2.y < p.y) # crossing
            if q1.x >= p.x
                if q2.x > p.x
                    c = !c
                elseif ((det3p(q1, q2, p) > 0) == (q2.y > q1.y))
                    c = !c
                end
            elseif q2.x > p.x
                if ((det3p(q1, q2, p) > 0) == (q2.y > q1.y))
                    c = !c
                end
            end
        end
    end
    return c
end

"""
    polysplit(p, p1, p2)

Split a polygon into two where it intersects with a line. It returns two
polygons:

```
(poly1, poly2)
```

This doesn't always work, of course. For example, a polygon the shape of the
letter "E" might end up being divided into more than two parts.
"""
function polysplit(pointlist::Array{Point, 1}, p1::Point, p2::Point)
    # the two-pass version
    # TODO should be one-pass
    newpointlist = Point[]; sizehint!(newpointlist, length(pointlist))
    vertex1 = Point(0, 0)
    vertex2 = Point(0, 0)
    l = length(pointlist)
    @inbounds for i in 1:l
        vertex1 = pointlist[mod1(i, l)]
        vertex2 = pointlist[mod1(i + 1, l)]
        flag, intersectpoint = intersectionlines(vertex1, vertex2, p1, p2, crossingonly=true)
        push!(newpointlist, vertex1)
        if flag
            push!(newpointlist, intersectpoint)
        end
    end
    # close?
    # push!(newpointlist, vertex2)
    # now sort points
    poly1 = Point[]
    poly2 = Point[]
    l = length(newpointlist)
    @inbounds for i in 1:l
        vertex1 = newpointlist[mod1(i, l)]
        d = pointlinedistance(vertex1, p1, p2)
        centerpoint = (p2.x - p1.x) * (vertex1.y - p1.y) > (p2.y - p1.y) * (vertex1.x - p1.x)
        if centerpoint
            push!(poly1, vertex1)
            abs(d) < 0.1 && push!(poly2, vertex1)
        else
            push!(poly2, vertex1)
            abs(d) < 0.1 && push!(poly1, vertex1)
        end
    end
    return(poly1, poly2)
end

"""
    prettypoly(points::Array{Point, 1}, action=:none, vertexfunction = () -> circle(O, 2, :stroke);
        close=false,
        reversepath=false,
        vertexlabels = (n, l) -> ()
        )

Draw the polygon defined by `points`, possibly closing and reversing it, using the current
parameters, and then evaluate the `vertexfunction` function at every vertex of the polygon.

The default vertexfunction draws a 2 pt radius circle.

To mark each vertex of a polygon with a randomly colored filled circle:

    p = star(O, 70, 7, 0.6, 0, vertices=true)
    prettypoly(p, :fill, () ->
        begin
            randomhue()
            circle(O, 10, :fill)
        end,
        close=true)

The optional keyword argument `vertexlabels` lets you supply a function with
two arguments that can access the current vertex number and the total number of vertices
at each vertex. For example, you can label the vertices of a triangle "1 of 3", "2 of 3",
and "3 of 3" using:

    prettypoly(triangle, :stroke,
        vertexlabels = (n, l) -> (text(string(n, " of ", l))))
"""
function prettypoly(pointlist::Array{Point, 1}, action=:none, vertexfunction = () -> circle(O, 2, :stroke);
    close=false,
    reversepath=false,
    vertexlabels = (n, l) -> ()
    )

    if isempty(pointlist)
        return nothing
    end

    if action != :path
        newpath()
    end

    if reversepath
        reverse!(pointlist)
    end

    move(pointlist[1])
    @inbounds for p in pointlist[2:end]
        line(p)
    end

    if close
        closepath()
    end
    do_action(action)
    pointnumber = 1

    for p in pointlist
        gsave()
        translate(p.x, p.y)
        vertexfunction()
        vertexlabels(pointnumber, length(pointlist))
        grestore()
        pointnumber += 1
    end

    if action == :fillpreserve || action == :strokepreserve
        move(pointlist[1])
        @inbounds for p in pointlist[2:end]
            line(p)
        end
    end

    return pointlist
end

function getproportionpoint(point::Point, segment, length, dx, dy)
    scalefactor = segment / length
    return Point((point.x - dx * scalefactor), (point.y - dy * scalefactor))
end

function drawroundedcorner(cornerpoint::Point, p1::Point, p2::Point, radius, path; debug=false)
    dx1 = cornerpoint.x - p1.x     # vector 1
    dy1 = cornerpoint.y - p1.y
    dx2 = cornerpoint.x - p2.x     # vector 2
    dy2 = cornerpoint.y - p2.y

    # Angle between vector 1 and vector 2 divided by 2
    angle2 = (atan(dy1, dx1) - atan(dy2, dx2)) / 2

    #  length of segment between corner point and the
    #  points of intersection with the circle of a given radius
    t = abs(tan(angle2))
    segment = radius / t

    # Check the segment
    length1 = hypot(dx1, dy1)
    length2 = hypot(dx2, dy2)
    seglength = min(length1, length2)
    if segment > seglength
        segment = seglength
        radius = seglength * t
    end

    #  points of intersection are calculated by the proportion between
    #  the coordinates of the vector, length of vector and the length of the segment.
    p1_cross = getproportionpoint(cornerpoint, segment, length1, dx1, dy1)
    p2_cross = getproportionpoint(cornerpoint, segment, length2, dx2, dy2)

    #  calculation of the coordinates of the circle's center by the addition of angular vectors
    dx = cornerpoint.x * 2 - p1_cross.x - p2_cross.x
    dy = cornerpoint.y * 2 - p1_cross.y - p2_cross.y
    L = hypot(dx, dy)
    d = hypot(segment, radius)

    # this prevents impossible constructions; Cairo will crash if L is 0
    if isapprox(L, 0.0)
        L= 0.01
    end

    circlepoint = getproportionpoint(cornerpoint, d, L, dx, dy)

    # if "debugging" or you just like the circles:
    debug && circle(circlepoint, radius, :stroke)

    # start angle and end engle of arc
    startangle = atan(p1_cross.y - circlepoint.y, p1_cross.x - circlepoint.x)
    endangle   = atan(p2_cross.y - circlepoint.y, p2_cross.x - circlepoint.x)

    # add first line segment, up to the start of the arc
    push!(path, (:line, p1_cross)) # draw line to arc start

    # adjust. Cairo also does this when you draw arc()s, btw
    if endangle < 0
        endangle = 2pi + endangle
    end
    if startangle < 0
        startangle = 2pi + startangle
    end
    sweepangle = endangle - startangle

    if abs(sweepangle) > pi
        if startangle < endangle
            push!(path, (:carc, [circlepoint, radius, startangle, endangle]))
        else
            push!(path, (:arc, [circlepoint, radius, startangle, endangle]))
        end
    else
        if startangle < endangle
            push!(path, (:arc, [circlepoint, radius, startangle, endangle]))
        else
            push!(path, (:carc, [circlepoint, radius, startangle, endangle]))
        end
    end
    # line from end of arc to start of next side
    push!(path, (:line,  p2_cross))
end

"""
    polysmooth(points, radius, action=:action; debug=false)

Make a closed path from the `points` and round the corners by making them arcs with the
given radius. Execute the action when finished.

The arcs are sometimes different sizes: if the given radius is bigger than the length of the
shortest side, the arc can't be drawn at its full radius and is therefore drawn as large as
possible (as large as the shortest side allows).

The `debug` option also draws the construction circles at each corner.
"""
function polysmooth(points::Array{Point, 1}, radius, action=:action; debug=false)
    temppath = Tuple[]
    l = length(points)
    if l < 3
        # there are less than three points to smooth
        return nothing
    else
        @inbounds for i in 1:l
            p1 = points[mod1(i, l)]
            p2 = points[mod1(i + 1, l)]
            p3 = points[mod1(i + 2, l)]
            drawroundedcorner(p2, p1, p3, radius, temppath, debug=debug)
        end
    end
    # need to close by joining to first point
    push!(temppath, temppath[1])
    # draw the path
    for (c, p) in temppath
        if c == :line
            line(p)    # add line segment
        elseif c == :arc
            arc(p...)  # add clockwise arc segment
        elseif c == :carc
            carc(p...) # add counterclockwise arc segment
        end
    end
    do_action(action)
end

"""
    offsetpoly(plist::Array{Point, 1}, d) where T<:Number

Return a polygon that is offset from a polygon by `d` units.

The incoming set of points `plist` is treated as a polygon, and another set of
points is created, which form a polygon lying `d` units away from the source
poly.

Polygon offsetting is a topic on which people have written PhD theses and
published academic papers, so this short brain-dead routine will give good
results for simple polygons up to a point (!). There are a number of issues to
be aware of:

- very short lines tend to make the algorithm 'flip' and produce larger lines

- small polygons that are counterclockwise and larger offsets may make the new
polygon appear the wrong side of the original

- very sharp vertices will produce even sharper offsets, as the calculated intersection point veers off to infinity

- duplicated adjacent points might cause the routine to scratch its head and wonder how to draw a line parallel to them
"""
function offsetpoly(plist::Array{Point, 1}, d::T) where T<:Number
    # don't try to calculate offset of two identical points
    l = length(plist)
    resultpoly = Array{Point}(undef, l)
    previouspoint = plist[end]
    for i in 1:l
        plist[i] == previouspoint && continue
        p1 = plist[mod1(i, l)]
        p2 = plist[mod1(i + 1, l)]
        p3 = plist[mod1(i + 2, l)]

        p1 == p2 && continue
        p2 == p3 && continue

        L12 = distance(p1, p2)
        L23 = distance(p2, p3)
        # the offset line of p1 - p2
        x1p = p1.x + (d * (p2.y - p1.y))/ L12
        y1p = p1.y + (d * (p1.x - p2.x))/ L12
        x2p = p2.x + (d * (p2.y - p1.y))/ L12
        y2p = p2.y + (d * (p1.x - p2.x))/ L12

        # the offset line of p2 - p3
        x3p = p2.x + (d * (p3.y - p2.y))/ L23
        y3p = p2.y + (d * (p2.x - p3.x))/ L23
        x4p = p3.x + (d * (p3.y - p2.y))/ L23
        y4p = p3.y + (d * (p2.x - p3.x))/ L23

        intersectionpoint = intersectionlines(
            Point(x1p, y1p),
            Point(x2p, y2p),
            Point(x3p, y3p),
            Point(x4p, y4p), crossingonly=false)

        if intersectionpoint[1]
            resultpoly[i] = intersectionpoint[2]
        end
        previouspoint = plist[i]
    end
    return resultpoly
end

"""
    offsetlinesegment(p1, p2, p3, d1, d2)

Given three points, find another 3 points that are offset by
d1 at the start and d2 at the end.

Negative d values put the offset on the left.
"""
function offsetlinesegment(p1, p2, p3, d1, d2)
    if p1 == p2 || p2 == p3 || p1 == p3
        throw(error("offsetlinesegment: the three points must be different"))
    end

    pt1 = perpendicular(p1, p2, -d1)
    pt2 = perpendicular(p3, p2, d2)

    L12 = distance(p1, p2)
    L23 = distance(p2, p3)

    d = (d1 + d2) / 2
    # the offset line of p1 - p2
    x1p = p1.x + (d * (p2.y - p1.y))/ L12
    y1p = p1.y + (d * (p1.x - p2.x))/ L12
    x2p = p2.x + (d * (p2.y - p1.y))/ L12
    y2p = p2.y + (d * (p1.x - p2.x))/ L12

    # the offset line of p2 - p3
    x3p = p2.x + (d * (p3.y - p2.y))/ L23
    y3p = p2.y + (d * (p2.x - p3.x))/ L23
    x4p = p3.x + (d * (p3.y - p2.y))/ L23
    y4p = p3.y + (d * (p2.x - p3.x))/ L23

    intersectionpoint = intersectionlines(
        Point(x1p, y1p),
        Point(x2p, y2p),
        Point(x3p, y3p),
        Point(x4p, y4p), crossingonly=false)

    if first(intersectionpoint)
        pt3 = intersectionpoint[2]
    else
        #  collinear probably
        pt3 = midpoint(pt1, pt2)
    end
    return pt1, pt3, pt2
end

"""
    offsetpoly(plist;
        startoffset = 10,
        endoffset   = 10,
        easingfunction = lineartween)

Return a closed polygon that is offset from and encloses an
open polygon.

The incoming set of points `plist` is treated as an open
polygon, and another set of points is created, which form a
polygon lying `...offset` units away from the source poly.

This method for `offsetpoly()` treats the list of points as
`n` vertices connected with `n - 1` lines. It allows you to
vary the offset from the start of the line to the end.

The other method `offsetpoly(plist, d)` treats the list of
points as `n` vertices connected with `n` lines.

# Extended help

This function accepts a keyword argument that allows you to
control the offset using a function, using the easing
functionality built in to Luxor. By default the function is
`lineartween()`, so the offset changes linearly between the
`startoffset` and the `endoffset`. The function:

```
f(a, b, c, d) = 2sin((a * π))
```

runs from 0 to 2 and back as `a` runs from 0 to 1.
The offsets are scaled by this amount.

"""
function offsetpoly(plist;
        startoffset = 10,
        endoffset   = 10,
        easingfunction = lineartween)

    l = length(plist)

    # TODO: special case a plist with 2 points
    l < 3 && throw(error("variableoffsetline: not enough points, need 3 or more"))
    # can't do 3 points properly, just insert a few extras
    if l == 3
        plist = vcat(plist[1], midpoint(plist[1], plist[2]), plist[2], midpoint(plist[2], plist[3]), plist[3])
    end

    # build the poly in two halves
    leftcurve  = Point[]
    rightcurve = Point[]

    pt1 = perpendicular(plist[1], plist[2], -startoffset)
    pt2 = perpendicular(plist[1], plist[2], startoffset)

    # start the curves off
    push!(leftcurve, pt1)
    push!(rightcurve, pt2)

    for i in 1:l-2
        # allow for the easing function that rescales the offset
        k1 = easingfunction(rescale(i,     0, l ), 0.0, 1.0, 1.0)
        k2 = easingfunction(rescale(i + 1, 0, l ), 0.0, 1.0, 1.0)

        # because easing functions are 0 to 1, use 0 here
        d1 = rescale(k1 * l, 0, l, startoffset, endoffset)
        d2 = rescale(k2 * l, 0, l, startoffset, endoffset)

        p1, mpt, p3 = offsetlinesegment(plist[i], plist[i + 1], plist[i + 2],  d1,  d2)
        push!(leftcurve, mpt)
        p1, mpt, p3 = offsetlinesegment(plist[i], plist[i + 1], plist[i + 2], -d1, -d2)
        push!(rightcurve, mpt)
    end
    # final point
    k = easingfunction(1, 0.0, 1.0, 1.0)
    d = rescale(k * l, 0, l, startoffset, endoffset)
    pt1 = perpendicular(plist[end], plist[end - 1], -d)
    pt2 = perpendicular(plist[end], plist[end - 1],  d)
    push!(leftcurve, pt2)
    push!(rightcurve, pt1)
    return vcat(leftcurve, reverse(rightcurve))
end

# third method

"""
offsetpoly(plist, shape::Function)

Return a closed polygon that is offset from and encloses an
polyline.

The incoming set of points `plist` is treated as an
polyline, and another set of points is created, which form a
closed polygon offset from the source poly.

This method for `offsetpoly()` treats the list of points as
`n` vertices connected with `n - 1` lines. (The other method
`offsetpoly(plist, d)` treats the list of points as `n`
vertices connected with `n` lines.)

The supplied function determines the width of the line.
`f(0, θ)` gives the width at the start (the slope of
the curve at that point is supplied in θ), `f(1, θ)` provides
the width at the end, and `f(n, θ)` is the width of point
`n/l`.

# Examples

This example draws a tilde, with the ends starting at 20
(10 + 10) units wide, swelling to 50 (10 + 10 + 15 + 15) in
the middle, as f(0.5) = 25.

```
f(x, θ) =  10 + 15sin(x * π)
sinecurve = [Point(50x, 50sin(x)) for x in -π:π/24:π]
pgon = offsetpoly(sinecurve, f)
poly(pgon, :fill)
```

This example enhances the vertical part of the curve, and
thins the horizontal parts.

```
g(x, θ) = rescale(abs(sin(θ)), 0, 1, 0.1, 30)
sinecurve = [Point(50x, 50sin(x)) for x in -π:π/24:π]
pgon = offsetpoly(sinecurve, g)
poly(pgon, :fill)
```
"""
function offsetpoly(plist, shape::Function)

    l = length(plist)

    # TODO: special case a plist with 2 points
    l < 3 && throw(error("variableoffsetline: not enough points, need 3 or more"))
    # can't do 3 points properly, just insert a few extras
    if l == 3
        plist = vcat(plist[1], midpoint(plist[1], plist[2]), plist[2], midpoint(plist[2], plist[3]), plist[3])
    end

    # build the poly in two halves
    leftcurve  = Point[]
    rightcurve = Point[]

    θ = slope(plist[1], plist[2])
    d = shape(0.0, θ)
    pt1 = perpendicular(plist[1], plist[2], -d)
    pt2 = perpendicular(plist[1], plist[2], d)

    # start the curves off
    push!(leftcurve, pt1)
    push!(rightcurve, pt2)

    for i in 1:l-2
        # allow for the easing function that rescales the offset
        θ = slope(plist[i], plist[i + 1])
        d = shape(rescale(i, 0, l-2), θ)
        p1, mpt, p3 = Luxor.offsetlinesegment(plist[i], plist[i + 1], plist[i + 2],  d,  d)
        push!(leftcurve, mpt)
        p1, mpt, p3 = Luxor.offsetlinesegment(plist[i], plist[i + 1], plist[i + 2], -d, -d)
        push!(rightcurve, mpt)
    end
    # final point
    θ = mod2pi(2π - slope(plist[end], plist[end - 1]))
    d = shape(1, θ)
    pt1 = perpendicular(plist[end], plist[end - 1], -d)
    pt2 = perpendicular(plist[end], plist[end - 1],  d)
    push!(leftcurve, pt2)
    push!(rightcurve, pt1)
    return vcat(leftcurve, reverse(rightcurve))
end

"""
    polyfit(plist::Array, npoints=30)

Build a polygon that constructs a B-spine approximation to it. The resulting list of points
makes a smooth path that runs between the first and last points.
"""
function polyfit(plist::Array{Point, 1}, npoints=30)
    l = length(plist)
    resultpoly = Array{Point}(undef, 0) ; sizehint!(resultpoly, npoints)
    # start at first point
    push!(resultpoly, plist[1])
    # skip the first point
    @inbounds for i in 2:l-1
        p1 = plist[mod1(i - 1, l)]
        p2 = plist[mod1(i, l)]
        p3 = plist[mod1(i + 1, l)]
        p4 = plist[mod1(i + 2, l)]
        a3 = (-p1.x + 3 * (p2.x - p3.x) + p4.x) / 6.0
        b3 = (-p1.y + 3 * (p2.y - p3.y) + p4.y) / 6.0
        a2 = (p1.x - 2p2.x + p3.x) / 2.0
        b2 = (p1.y - 2p2.y + p3.y) / 2.0
        a1 = (p3.x - p1.x) / 2.0
        b1 = (p3.y - p1.y) / 2.0
        a0 = (p1.x + 4p2.x + p3.x) / 6.0
        b0 = (p1.y + 4p2.y + p3.y) / 6.0
        for i in 1:l-1
            t = i/npoints
            x = ((((a3 * t + a2) * t) + a1) * t) + a0
            y = ((((b3 * t + b2) * t) + b1) * t) + b0
            push!(resultpoly, Point(x, y))
        end
    end
    # finish at last point
    push!(resultpoly, plist[end])
    return resultpoly
end

"""
    pathtopoly()

Convert the current path to an array of polygons.

Returns an array of polygons, corresponding to the paths and subpaths of the original path.
"""
function pathtopoly()
    originalpath = getpathflat()
    polygonlist = Array{Point, 1}[] ; sizehint!(polygonlist, length(originalpath))
    if length(originalpath) > 0
        pointslist = Point[]
        for e in originalpath
            if e.element_type == Cairo.CAIRO_PATH_MOVE_TO                # 0
                push!(pointslist, Point(first(e.points), last(e.points)))
            elseif e.element_type == Cairo.CAIRO_PATH_LINE_TO            # 1
                push!(pointslist, Point(first(e.points), last(e.points)))
            elseif e.element_type == Cairo.CAIRO_PATH_CLOSE_PATH         # 3
                if last(pointslist) == first(pointslist)
                    # don’t repeat first point, we can close it ourselves
                    if length(pointslist) > 2
                        pop!(pointslist)
                    end
                end
                push!(polygonlist, pointslist)
                pointslist = Point[]
            else
                error("pathtopoly(): unknown CairoPathEntry " * repr(e.element_type))
                error("pathtopoly(): unknown CairoPathEntry " * repr(e.points))
            end
        end
        # the path was never closed, so flush
        if length(pointslist) > 0
            push!(polygonlist, pointslist)
        end
    end
    return polygonlist
end

"""
    polydistances(p::Array{Point, 1}; closed=true)

Return an array of the cumulative lengths of a polygon.
"""
function polydistances(p::Array{Point, 1}; closed=true)
    l = length(p)
    r = Float64[0.0] ; sizehint!(r, l)
    t = 0.0
    @inbounds for i in 1:l - 1
        t += distance(p[i], p[i + 1])
        push!(r, t)
    end
    if closed
        t += distance(p[end], p[1])
        push!(r, t)
    end
    return r
end

"""
    polyperimeter(p::Array{Point, 1}; closed=true)

Find the total length of the sides of polygon `p`.
"""
function polyperimeter(p::Array{Point, 1}; closed=true)
    return polydistances(p, closed=closed)[end]
end

"""
    nearestindex(polydistancearray, value)

Return a tuple of the index of the largest value in `polydistancearray` less
than `value`, and the difference value. Array is assumed to be sorted.

(Designed for use with `polydistances()`).
"""
function nearestindex(a::Array{T, 1} where T <: Real, val)
    ind = findlast(v -> (v < val), a)
    surplus = 0.0
    if ind > 0.0
        surplus = val - a[ind]
    else
        ind = 1
    end
    return (ind, surplus)
end

"""
    polyportion(p::Array{Point, 1}, portion=0.5; closed=true, pdist=[])

Return a portion of a polygon, starting at a value between 0.0 (the beginning)
and 1.0 (the end). 0.5 returns the first half of the polygon, 0.25 the first
quarter, 0.75 the first three quarters, and so on.

Use `closed=false` to exclude the line joining the
final point to the first point from the calculations.

If you already have a list of the distances between each point in the polygon
(the "polydistances"), you can pass them in `pdist`, otherwise they'll be
calculated afresh, using `polydistances(p, closed=closed)`.

Use the complementary `polyremainder()` function to return the other part.
"""
function polyportion(p::Array{Point, 1}, portion=0.5; closed=true, pdist=[])
    # portion is 0 to 1
    if isempty(pdist)
        pdist = polydistances(p, closed=closed)
    end
    portion = clamp(portion, 0.0, 1.0)
    # don't bother to do 0.0
    isapprox(portion, 0.0, atol=0.00001) && return p[1:1]
    # don't bother to do 1.0
    if closed == false && isapprox(portion, 1.0, atol=0.00001)
        return p
    elseif isapprox(portion, 1.0, atol=0.00001)
        return p
    end
    ind, surplus = nearestindex(pdist, portion * pdist[end])
    if surplus > 0.0
        nextind = mod1(ind + 1, length(p))
        overshootpoint = between(p[ind], p[nextind], surplus/distance(p[ind], p[nextind]))
        return vcat(p[1:ind], overshootpoint)
    else
        return p[1:end]
    end
end

"""
    polyremainder(p::Array{Point, 1}, portion=0.5; closed=true, pdist=[])

Return the rest of a polygon, starting at a value between 0.0 (the beginning)
and 1.0 (the end). 0.5 returns the last half of the polygon, 0.25 the last three
quarters, 0.75 the last quarter, and so on.

Use `closed=false` to exclude the line joining the
final point to the first point from the calculations.

If you already have a list of the distances between each point in the polygon
(the "polydistances"), you can pass them in `pdist`, otherwise they'll be
calculated afresh, using `polydistances(p, closed=closed)`.

Use the complementary `polyportion()` function to return the other part.
"""
function polyremainder(p::Array{Point, 1}, portion=0.5; closed=true, pdist=[])
    # portion is 0 to 1
    if isempty(pdist)
        pdist = polydistances(p, closed=closed)
    end
    portion = clamp(portion, 0.0, 1.0)

    # don't bother to do 0.0
    isapprox(portion, 0.0, atol=0.00001) && return p

    # don't bother to do 1.0
    if isapprox(portion, 1.0, atol=0.00001)
        return p[1:1]
    end
    ind, surplus = nearestindex(pdist, portion * pdist[end])
    if surplus > 0.0
        nextind = mod1(ind + 1, length(p))
        overshootpoint = between(p[ind], p[nextind], surplus/distance(p[ind], p[nextind]))
        return vcat(overshootpoint, p[nextind:end])
    else
        return p[1:end]
    end
end

"""
    polysample(p::Array{Point, 1}, npoints::T where T <: Integer;
            closed=true)

Sample the polygon `p`, returning a polygon with `npoints` to represent it. The
first sampled point is:

     1/`npoints` * `perimeter of p`

away from the original first point of `p`.

If `npoints` is the same as `length(p)` the returned polygon is the same as the
original, but the first point finishes up at the end (so `new=circshift(old,
1)`).

If `closed` is true, the entire polygon (including the edge joining the last
point to the first point) is sampled.
"""
function polysample(p::Array{Point, 1}, npoints::T where T <: Integer;
        closed=true)
    l = length(p)
    l < 2 && error("polysample(): not enough points in polygon to take samples")
    npoints < 2  && return p[[1, end]]
    distances = polydistances(p, closed=closed)
    result = Point[] ; sizehint!(result, l)
    for i in 1:npoints
        ind, surplus = nearestindex(distances, (i/npoints) * distances[end])
        if surplus > 0.0
            nextind = mod1(ind + 1, l)
            overshootpoint = between(p[ind], p[nextind], surplus/distance(p[ind], p[nextind]))
            push!(result, overshootpoint)
        else
            push!(result, p[i])
        end
    end
    return result
end

"""
    polyarea(p::Array)

Find the area of a simple polygon. It works only for polygons that don't
self-intersect. See also `polyorientation()`.
"""
function polyarea(plist::Array{Point, 1})
    l = length(plist)
    area = 0.0
    @inbounds for i in eachindex(plist)
        j = mod1(i + 1, l)
        area += plist[i].x * plist[j].y
        area -= plist[j].x * plist[i].y
    end
    area = abs(area) / 2.0
    return area
end

"""
    intersectlinepoly(pt1::Point, pt2::Point, C)

Return an array of the points where a line between pt1 and pt2 crosses polygon C.
"""
function intersectlinepoly(pt1::Point, pt2::Point, C::Array{Point, 1})
    intersectingpoints = Point[]
    @inbounds for j in eachindex(C)
        Cpointpair = (C[j], C[mod1(j + 1, length(C))])
        flag, pt = intersectionlines(pt1, pt2, Cpointpair..., crossingonly=true)
        if flag
            push!(intersectingpoints, pt)
        end
    end
    # sort by distance from pt1
    sort!(intersectingpoints, lt = (p1, p2) -> distance(p1, pt1) < distance(p2, pt1))
    return intersectingpoints
end

"""
    polyintersections(S::Array{Point, 1}, C::Array{Point, 1})

Return an array of the points in polygon S plus the points where polygon S crosses
polygon C. Calls `intersectlinepoly()`.

TODO This code is experimental...
"""
function polyintersections(S::Array{Point, 1}, C::Array{Point, 1})
    Splusintersectionpoints = Point[] ; sizehint!(Splusintersectionpoints, length(S) + length(C))
    @inbounds for i in eachindex(S)
        Spointpair = (S[i], S[mod1(i+1, length(S))])
        push!(Splusintersectionpoints, S[i])
        for pt in intersectlinepoly(Spointpair..., C)
            push!(Splusintersectionpoints, pt)
        end
    end
    return Splusintersectionpoints
end

# TODO these experimental functions don't work all the time
# use with caution... :)

"""
    polyorientation(pgon)

Returns a number which is positive if the polygon is clockwise in Luxor...

TODO This code is still experimental...
"""
function polyorientation(pgon::Array{Point, 1})
    # in Luxor polys are usually clockwise
    # perhaps this is because the Y axis goes down...
    sum = 0.0
    @inbounds for i in 1:length(pgon)
        sum += crossproduct(pgon[i], pgon[mod1(i + 1, end)])
    end
    return sum
end

polyorientation(pt1, pt2, pt3) = polyorientation(Point[pt1, pt2, pt3])

"""
    ispolyclockwise(pgon)

Returns true if polygon is clockwise. WHEN VIEWED IN A LUXOR DRAWING...?

TODO This code is still experimental...
"""
function ispolyclockwise(pgon::Array{Point, 1})
    polyorientation(pgon) > 0.0
end

"""
    ispointinsidetriangle(p, p1, p2, p3)
    ispointinsidetriangle(p, triangle::Array{Point, 1})

Returns false if `p` is not inside triangle p1 p2 p3.
"""

function ispointinsidetriangle(p::Point, p1::Point, p2::Point, p3::Point)
    if polyorientation([p1, p2, p3]) < 0.0
        p1, p3 = p3, p1
    end
    # barycentric method?
    s = p1.y * p3.x - p1.x * p3.y + (p3.y - p1.y) * p.x + (p1.x - p3.x) * p.y
    t = p1.x * p2.y - p1.y * p2.x + (p1.y - p2.y) * p.x + (p2.x - p1.x) * p.y

    if s < 0 != t < 0
        return false
    end
    A = -p2.y * p3.x + p1.y * (p3.x - p2.x) + p1.x * (p2.y - p3.y) + p2.x * p3.y
    return A < 0 ?
        (s <= 0 && s + t >= A) :
        (s >= 0 && s + t <= A)
end

ispointinsidetriangle(p::Point, triangle::Array{Point, 1}) =
    ispointinsidetriangle(p, triangle[1], triangle[2], triangle[3])

"""
    polyremovecollinearpoints(pgon::Array{Point, 1})

Return copy of polygon with no collinear points.

Caution: may return an empty polygon... !

TODO This code is still experimental...
"""
function polyremovecollinearpoints(pgon::Array{Point, 1})
    markfordeletion = []
    @inbounds for n in 1:length(pgon)
        p1 = pgon[n]

        pfirst = pgon[mod1(n - 1, length(pgon))]
        plast = pgon[mod1(n + 1, length(pgon))]

        if ispointonline(p1, pfirst, plast, extended=true, atol=0.1)
            push!(markfordeletion, n)
        end
    end
    return pgon[setdiff(1:length(pgon), markfordeletion)]
end

"""
    polymove!(pgon, frompoint::Point, topoint::Point)

Move (permanently) a polygon from `frompoint` to `topoints`.
"""
function polymove!(pgon, frompoint::Point, topoint::Point)
    d = topoint - frompoint
    @inbounds for i in eachindex(pgon)
        pgon[i] = Point(pgon[i].x + d.x, pgon[i].y + d.y)
    end
    return pgon
end

"""
    polyscale!(pgon, s;
       center=O)

Scale (permanently) a polygon by `s`, relative to `center`.
"""
function polyscale!(pgon, s;
        center=O)
    @inbounds for i in eachindex(pgon)
        pgon[i] = between(center, pgon[i], s)
    end
    return pgon
end

"""
    polyscale!(pgon, sh, sv;
        center=O)

Scale (permanently) a polygon by `sh` horizontally and `sv` vertically,
relative to `center`.
"""
function polyscale!(pgon, sh, sv;
        center=O)
    @inbounds for i in eachindex(pgon)
        pgon[i] = (pgon[i] - center) * (sh, sv)
    end
    return pgon
end

"""
    polyrotate!(pgon, θ;
        center=O)

Rotate (permanently) a polygon around `center` by `θ` radians.
"""
function polyrotate!(pgon, θ;
        center=O)
    costheta = cos(θ)
    sintheta = sin(θ)
    @inbounds for i in eachindex(pgon)
        pgon[i] = Point(
            (pgon[i].x - center.x) * costheta  -  ((pgon[i].y - center.y) * sintheta) + center.x,
            (pgon[i].x - center.x) * sintheta  +  ((pgon[i].y - center.y) * costheta) + center.y)
    end
    return pgon
end

"""
    polyreflect!(pgon, pt1 = O, pt2 = O + (0, 100)

Reflect (permanently) a polygon in a line (default to the y-axis)
joining two points.

"""
function polyreflect!(pgon, pt1=O, pt2 = O + (0, 100))
    @inbounds for i in eachindex(pgon)
        gnp = getnearestpointonline(pt1, pt2, pgon[i])
        pgon[i] = between(pgon[i], gnp, 2.0)
    end
    return pgon
end

"""
    insertvertices!(pgon;
        ratio=0.5)

Insert a new vertex into each edge of a polygon `pgon`. The default `ratio` of
0.5 divides the original edge of the polygon into half.
"""
function insertvertices!(pgon;
        ratio=0.5)
    counter = 1
    vertexnumber = 2
    originallength = length(pgon)
    while true
        insert!(pgon, vertexnumber,
        between(pgon[mod1(vertexnumber - 1, length(pgon))],
            pgon[mod1(vertexnumber, length(pgon))],
            ratio))
        vertexnumber += 2
        counter += 1
        if counter > originallength
            break
        end
    end
    return pgon
end

"""
    polyintersect(p1::AbstractArray{Point, 1}, p2::AbstractArray{Point, 1};
        closed=true)

TODO: Fix/test/improve this experimental polygon intersection routine.

Return the points where polygon p1 and polygon p2 cross.

If `closed` is false, the intersection points must lie on the first `n - 1` lines of each polygon.

"""
function polyintersect(p1::AbstractArray{Point, 1}, p2::AbstractArray{Point, 1};
    closed=true)
    length(p1) < 3 || length(p2) < 3 && error("polyintersect(): not enough points")
    temp = Point[]
    @inbounds for i in eachindex(p1)
        Spointpair = (p1[i], p1[mod1(i + 1, length(p1))])
        for pt in intersectlinepoly(Spointpair..., p2)
            push!(temp, pt)
        end
    end
    # if not closed polygons, remove ipts that are on close side
    if closed == false
        ripts = Point[]
        for ipt in temp
            if !(ispointonline(ipt, p1[end], p1[1]) || ispointonline(ipt, p2[end], p2[1]) )
                push!(ripts, ipt)
            end
        end
        return ripts
    else
        return temp
    end
end

# triangulation functions

function _smallesttriangle(bb::BoundingBox)
    # find the smallest triangle that completely encloses bounding box
    diag = boxdiagonal(bb)

    # circle that encloses this box
    circ = (center = boxmiddlecenter(bb), radius = diag / 2)

    # three equally spaced points on circle
    side1pt = circ.center + polar(circ.radius, 0)
    side2pt = circ.center + polar(circ.radius, 2π / 3)
    side3pt = circ.center + polar(circ.radius, 4π / 3)

    # make them clockwise
    pgon = [side1pt, side2pt, side3pt]
    if !ispolyclockwise(pgon)
        reverse!(pgon)
    end
    return offsetpoly(pgon, circ.radius / 2)
end

mutable struct TriEdge
    spt::Point
    ept::Point
    valid::Bool
end

function _edgesequal(edge1::TriEdge, edge2::TriEdge)
    !edge1.valid && return false
    !edge2.valid && return false
    if ((edge1.spt ≈ edge2.spt) && (edge1.ept ≈ edge2.ept)) ||
       ((edge1.spt ≈ edge2.ept) && (edge1.ept ≈ edge2.spt))
        return true
    else
        return false
    end
end

"""
    polytriangulate(plist::Array{Point,1}; epsilon = -0.01)

Triangulate the polygon in `plist`.

This uses the Bowyer–Watson/Delaunay algorithm to make triangles. It returns an array of triangular polygons.

TODO: This experimental polygon function is not very efficient, because it first copies the list of points (to avoid modifying the original), and sorts it, before making triangles.
"""
function polytriangulate(plist::Array{Point,1}; epsilon = -0.001)
    trianglelist = Vector{Point}[]
    pointlist = deepcopy(plist)
    vertexcount = length(pointlist)

    # make enclosing supertriangle
    bb = BoundingBox(pointlist)
    supertriangle = Luxor._smallesttriangle(bb)

    # add the supertriangle to the trianglelist
    push!(trianglelist, supertriangle)

    # add the supertriangle’s vertices to the point list
    push!(pointlist, supertriangle[1])
    push!(pointlist, supertriangle[2])
    push!(pointlist, supertriangle[3])

    # sorting the list of points by x coordinate
    sort!(pointlist)

    @inbounds for point in pointlist
        edgebuffer = Vector{Luxor.TriEdge}() # to store edges
        tobedeleted = Int64[]

        for ntriangle in eachindex(trianglelist)
            # for each triangle currently in the trianglelist
            # calculate the triangle circumcircle center and radius
            triangle = trianglelist[ntriangle]
            cp, r = center3pts(triangle[1], triangle[2], triangle[3])

            #as soon as the x component of the distance from the current point to the
            #circumcircle center is greater than the circumcircle radius, that triangle need
            #never be considered for later points, as further points will never again be on
            #the interior of that triangles circumcircle.
            abs(point.x - cp.x) > r && continue

            # if this fails, r is zero, so no problem with adding edges/deleting triangles
            if distance(point, cp) < (r + epsilon)
                # if the point lies in the triangle’s circumcircle,
                # add the triangle’s edges to the edge buffer
                push!(edgebuffer, Luxor.TriEdge(triangle[1], triangle[2], true))
                push!(edgebuffer, Luxor.TriEdge(triangle[2], triangle[3], true))
                push!(edgebuffer, Luxor.TriEdge(triangle[3], triangle[1], true))
                # and mark the triangle to be deleted
                push!(tobedeleted, ntriangle)
            end
        end
        deleteat!(trianglelist, tobedeleted)

        # find all shared edges in the edge buffer
        # when they’re removed, the edges of the enclosing polygon are left
        j = 1
        while j <= length(edgebuffer)
            k = j + 1
            while k <= length(edgebuffer)
                if Luxor._edgesequal(edgebuffer[j], edgebuffer[k]) == true
                    edgebuffer[j].valid = false
                    edgebuffer[k].valid = false
                end
                k += 1
            end
            j += 1
        end

        # make new triangles around the current point
        # skip over any tagged edges
        j = 1
        while j <= length(edgebuffer)
            if edgebuffer[j].valid
                # the order is important
                pgon = sort([edgebuffer[j].ept, point, edgebuffer[j].spt])
                # make sure every polygon is clockwise
                if !ispolyclockwise(pgon)
                    reverse!(pgon)
                end
                push!(trianglelist, pgon)
            end
            j += 1
        end
    end # do next point

    # # tidy up
    # # remove any triangles that use the supertriangle’s vertices
    indexes = Int64[]
    for n in eachindex(trianglelist)
        triangle = trianglelist[n]
        for tpt in triangle
            for stpt in supertriangle
                if tpt ≈ stpt
                    push!(indexes, n) # mark for deletion
                end
            end
        end
    end
    deleteat!(trianglelist, unique(indexes))
    return trianglelist
end


text:
include("runtests.jl")
include("test_extensions.jl") # NB will currently clean out the 'output' folder! - tests number of model extension output files for PAGE2020 update

ext: jl


text:
using PackageCompiler

create_app(
   @__DIR__,
   join((@__DIR__, "Compiled")),
   precompile_execution_file="precompile_app.jl",
   audit=true, # Warn about eventual relocatability problems with the app
   force=false,
)

max_stars_repo_name: arthursoprana/FlowPattern

max_issues_repo_name: arthursoprana/FlowPattern

max_forks_repo_name: arthursoprana/FlowPattern


text:
using Test

@testset "Autotuning.jl" begin
end


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

# get next call from list
function getNextCall!(queuedCallList::Vector{Call})
	return length(queuedCallList) > 0 ? pop!(queuedCallList) : nothing
end

# reset calls
function reset!(calls::Vector{Call})
	# do not reset fixed call parameters
	fnamesFixed = (:index, :priority, :transport, :hospitalIndex, :location,
		:arrivalTime, :dispatchDelay, :onSceneDuration, :handoverDuration,
		:nearestNodeIndex, :nearestNodeDist)
	fnames = setdiff(fieldnames(Call), fnamesFixed)

	# reset calls
	nullCall = Call()
	recentCallIndex = something(findlast(call -> call.status != callNullStatus, calls), 0) # index of last call that needs reset
	for fname in fnames
		val = getfield(nullCall, fname)
		if isprimitivetype(fieldtype(Call, fname))
			for i = 1:recentCallIndex
				setfield!(calls[i], fname, val)
			end
		elseif fieldtype(Call, fname) == Location
			for i = 1:recentCallIndex
				copy!(getfield(calls[i], fname), val) # faster than deepcopy for Location
			end
		else
			for i = 1:recentCallIndex
				setfield!(calls[i], fname, deepcopy(val))
			end
		end
	end
end
resetCalls!(sim::Simulation) = reset!(sim.calls) # compat

function setCallStatus!(call::Call, status::CallStatus, time::Float)
	@assert(call.statusSetTime <= time)

	# stats for previous status
	statusDuration = time - call.statusSetTime
	prevStatus = call.status # shorthand
	if prevStatus == callQueued
		call.queuedDuration += statusDuration
	elseif prevStatus == callWaitingForAmb
		call.waitingForAmbDuration += statusDuration
		call.ambGoingToCallDuration = statusDuration # overwrite previous value if have attempted dispatch to call multiple times
	elseif prevStatus == callGoingToHospital
		@assert(call.transportDuration == 0.0) # should not have been calculated yet
		call.transportDuration = statusDuration
	end

	# stats for new status
	if status == callQueued
		call.wasQueued = true
	elseif status == callWaitingForAmb
		if prevStatus == callScreening
			@assert(isapprox(time, call.arrivalTime + call.dispatchDelay)) # if dispatching after call screening, only delay should be dispatch delay
		end
		call.dispatchTime = time
	elseif status == callOnSceneTreatment
		@assert(call.ambArrivalTime == nullTime) # value should not have been set yet
		call.ambArrivalTime = time
		call.responseDuration = time - call.arrivalTime
		call.bumpedDuration = call.waitingForAmbDuration - call.ambGoingToCallDuration
	elseif status == callAtHospital
		call.hospitalArrivalTime = time
	elseif status == callProcessed
		@assert(call.serviceDuration == 0.0) # value should not have been set yet
		call.serviceDuration = time - call.arrivalTime
	end

	call.status = status
	call.statusSetTime = time
end

function initCalls!(sim::Simulation, calls::Vector{Call})
	for (i,call) in enumerate(calls)
		call.index = i

		# find nearest nodes
		if call.nearestNodeIndex == nullIndex
			(call.nearestNodeIndex, call.nearestNodeDist) = findNearestNode(sim.map, sim.grid, sim.net.fGraph.nodes, call.location)
		end
	end
end

# Set the calls of the simulation.
# Note that the old sim.calls will not be kept.
# I tried to write this function to allow sim.startTime to be changed,
# but that would also require sim stats capturing to change, so startTime cannot be changed.
# Mutates: sim, calls
function setSimCalls!(sim::Simulation, calls::Vector{Call})
	@assert(length(calls) >= 1)
	@assert(sim.startTime <= calls[1].arrivalTime)

	reset!(sim) # for safety, as changing calls part-way through a simulation could be tricky

	sim.numCalls = length(calls)

	initCalls!(sim, calls)

	# set first event to match first call arrival time; hacky
	firstCallEvents = filter(e -> e.form == callArrives, sim.eventList)
	@assert(length(firstCallEvents) == 1)
	firstCallEvents[1].time = calls[1].arrivalTime

	if isdefined(sim, :backup)
		sim.backup.numCalls = sim.numCalls
		sim.backup.eventList = deepcopy(sim.eventList)
	end

	reset!(calls)

	sim.calls = calls
end


text: const SSubArray{S,T,N,P,I,L} = SizedArray{S,T,N,N,SubArray{T,N,P,I,L}}


text:
# Autogenerated wrapper script for VideoCaptureWrap_jll for x86_64-w64-mingw32-cxx11
export libvideocapture

using libcxxwrap_julia_jll
using OpenCVQt_jll
JLLWrappers.@generate_wrapper_header("VideoCaptureWrap")
JLLWrappers.@declare_library_product(libvideocapture, "libvideocapture.dll")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll, OpenCVQt_jll)
    JLLWrappers.@init_library_product(
        libvideocapture,
        "bin\\libvideocapture.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
Fexp(x)=(exp(x),nothing)

function Fexp_inplace(value,auxvalue,x)
    copy!(value,exp(x))
end

function Dexp(derivativeIDX,f_c,faux_c,grad_c,grad_n,x)
    axpy!(1.0,grad_c.*f_c,grad_n)
end

if PROC=="GPU"

    function Fexp(x::CudaArray)
        tmp=CudaArray(Float64,size(x))
        exp!(x,tmp)
        return (tmp,nothing) # memory leak here
    end

    function Fexp_inplace(value,auxvalue,x::CudaArray)
        exp!(x,value)
    end

    function Dexp(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray)
        vmultupdate!(1.0,grad_c,f_c,grad_n)
    end

end

Derivative[Fexp]=Dexp
Inplace[Fexp]=Fexp_inplace

import Base.exp

exp(n::ADnode)=ADnode(Fexp,n)

export exp

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using SeisConvert
using Test

@testset "SeisConvert.jl" begin
    # Write your own tests here.
end


text:
using Test
#@startexample Global BPS (Truncated Gaussian)
#=
In this example we use the global Bouncy Particle Sampler on 2D Gaussian truncated to the positive orthan to show how to declare a BPS model.

![](../assets/truncatedgaussian.png)

Start by loading the library:
=#
using PDSampler, Random
#=
you then need to define two elements:
1. a geometry (boundaries)
2. an energy (gradient of the log-likelihood of the target)

The positive orthan corresponds to a simple Polygonal domain where the boundaries are the axes. The normal to these boundaries (`ns`) are therefore unit vectors and the intercepts (`a`) are zero.
A polygonal domain is then declared with the constructor `Polygonal`.
=#
p = 2
# normal to faces and intercepts
ns, a = diagm(0=>ones(p)), zeros(p)
geom  = Polygonal(ns, a)
#=
The function `nextboundary` returns a function that can compute the next boundary on the current ray `[x,x+tv]` with `t>0` as well as the time of the hit.
=#
nextbd(x, v) = nextboundary(geom, x, v)
#=
The model then needs to be specified: you need to define a function of the form
`gradll(x)` which can return the gradient of the log-likelihood at some point `x`.
Here, let us consider a 2D gaussian.
=#
# here we build a valid precision matrix. The cholesky decomposition of
# the covariance matrix will be useful later to build a sensible
# starting point for the algorithm.
Random.seed!(12)
P1  = randn(p,p)
P1 *= P1'
P1 += norm(P1)/100*diagm(0=>ones(p))
C1  = inv(P1); C1 += C1'; C1/=2;
L1  = cholesky(C1).L
mu  = zeros(p) .+ 1.
mvg = MvGaussianCanon(mu, P1)
#=
Here, we have defined the gaussian through the *Canonical* representation i.e.: by specifying a mean and a precision matrix.

Every model must implement a `gradloglik` function returning the gradient of the log-likelihood at a point `x`.
=#
gradll(x) = gradloglik(mvg, x)
#=
Next, you need to define the function which can return the first arrival time of
the corresponding Inhomogenous Poisson Process.

Note that you could be using `nextevent_zz` here as well if you wanted to use
the Zig-Zag sampler (and you could implement other kernels as well).
=#
nextev(x, v) = nextevent_bps(mvg, x, v)
#=
For a Gaussian (and some other simple distributions), this is analytical through
an inversion-like method.

Finally, you need to specify the parameters of the simulation such as the
starting point and velocity, the length of the path generated, the rate of
refreshment and the maximum number of gradient evaluations.
=#
T    = 1000.0   # length of path generated
lref = 2.0      # rate of refreshment
x0   = mu+L1*randn(p) # sensible starting point
v0   = randn(p) # starting velocity
v0  /= norm(v0) # put it on the sphere (not necessary)
# Define a simulation
sim = Simulation( x0, v0, T, nextev, gradll,
                  nextbd, lref ; maxgradeval = 10000)
#=
And finally, generate the path and recover some details about the simulation.
=#
(path, details) = simulate(sim)
#=
The `path` object belongs to the type `Path` and can be sampled using
`samplepath`.

A crude sanity check is for example to check that the estimated mean obtained through quadrature along the path yields a similar result as a basic Monte Carlo estimator.
=#
# Building a basic MC estimator
# (taking samples from 2D MVG that are in positive orthan)
sN = 1000
s  = broadcast(+, mu, L1*randn(p,sN))
mt = zeros(2)
np = 0
# Sum for all samples in the positive orthan
ss = [s; ones(sN)']
mt = sum(ss[:,i] for i in 1:sN if !any(e->e<0, ss[1:p,i]))
mt = mt[1:p]/mt[end]
#=
You can now compare the norm of `mt` (a crude MC estimator) to `pathmean(path)` (computing the integrals along the segments of the path) and you will see that
the relative error is below 5%.
=#
#@endexample

# all that's below won't show in doc

@test norm(pathmean(path)-mt)/norm(mt) < 0.05
ess,ns = esspath(path; ns=100)
@test minimum(ess) > 150

# Additional testing: this is for the quadrature test

nextev2(x, v, t) = nextevent_bps_q(gradll, x, v, t; n=20)
sim2 = Simulation( x0, v0, T, nextev2, gradll,
                  nextbd, lref ; maxgradeval = 10000)

(path2, details2) = simulate(sim2)

@test norm(pathmean(path2)-mt)/norm(mt) < 0.05


text:
"""
Plot the performance statistics of greedy-multi-path with respect to edge percolation
rate (The probability that a given edge is removed)
"""

using QuNet
using LaTeXStrings
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using JLD
using Parameters

datafile = "data/percolations"

# Params
num_pairs = 1::Int64
grid_size = 10::Int64
perc_range = (0.0, 0.01, 0.7)::Tuple{Float64, Float64, Float64}
# 5000
num_trials = 5000::Int64

generate_new_data = false
if generate_new_data == true

    net = GridNetwork(grid_size, grid_size)

    perf_data = []
    perf_err = []
    path_data = []
    path_err = []
    ave_path_distance = []
    ave_path_distance_err = []
    postman = []
    postman_err = []

    for p in perc_range[1]:perc_range[2]:perc_range[3]
        println("Collecting for percolation rate: $p")

        # Percolate the network
        # perc_net = QuNet.percolate_edges(net, p)
        # refresh_graph!(perc_net)

        # Collect performance data with error, percolating the network edges
        p, p_e, pat, pat_e, dummy, dummy, apd, apd_err, post, post_err = net_performance(net, num_trials, num_pairs, max_paths=4,
        edge_perc_rate = p, get_apd=true, get_postman=true)
        push!(perf_data, p)
        push!(perf_err, p_e)
        push!(path_data, pat)
        push!(path_err, pat_e)
        push!(ave_path_distance, apd)
        push!(ave_path_distance_err, apd_err)
        push!(postman, post)
        push!(postman_err, post_err)
    end

    # Save data
    d = Dict{Symbol, Any}()
    @pack! d = num_pairs, grid_size, perc_range, num_trials, perf_data, perf_err, path_data, path_err, ave_path_distance,
    ave_path_distance_err, postman, postman_err
    save("$datafile.jld", "data", d)

else
    # Load data
    if !isfile("$datafile.jld")
        error("$datafile.jld not found in working directory")
    end
    d = load("$datafile.jld")["data"]
    @unpack num_pairs, grid_size, perc_range, num_trials, perf_data, perf_err, path_data, path_err, ave_path_distance, ave_path_distance_err,
    postman, postman_err = d
end

# Get values for x axis
x = collect(perc_range[1]:perc_range[2]:perc_range[3])

# Extract performance data
loss = collect(map(x->x["loss"], perf_data))
z = collect(map(x->x["Z"], perf_data))
loss_err = collect(map(x->x["loss"], perf_err))
z_err = collect(map(x->x["Z"], perf_err))

# Extract data from path: PX is the rate of using X paths
P0 = [path_data[i][1]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P1 = [path_data[i][2]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P2 = [path_data[i][3]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P3 = [path_data[i][4]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P4 = [path_data[i][5]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]

# Extract errors from path:
P0e = [path_err[i][1]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P1e = [path_err[i][2]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P2e = [path_err[i][3]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P3e = [path_err[i][4]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]
P4e = [path_err[i][5]/num_pairs for i in 1:length(perc_range[1]:perc_range[2]:perc_range[3])]

pathratio =  postman ./ ave_path_distance

##### Plot cost_percolation #####
msize = 5
# color_palette = palette(:lightrainbow, 4),
plot(x, loss, ylims=(0,1), seriestype = :scatter, yerror = loss_err, label=L"\eta",
legend=:top, guidefontsize=14, tickfontsize=12, legendfontsize=8, fontfamily="computer modern",
markersize=msize, size=(600,400), color = :Dodgerblue, markershape=:utriangle)
ylabel!("Costs")
plot!(x, z, seriestype = :scatter, yerror = z_err, markersize=msize, label=L"F",
color = :Crimson)
# Adding new ave_shortest_path costs
f = dB_to_Z.(ave_path_distance)
e = dB_to_P.(ave_path_distance)
plot!(x, f, linewidth=4, markersize=msize, label=L"$\textrm{Conditional Manhattan } F$",
color=:Tomato)
plot!(x, e, linewidth=4, markersize=msize, label=L"$\textrm{Conditional Manhattan } \eta$",
color=:LightSkyBlue)
# TODO: Try cost scaled by path_ratio
# xaxis!(L"$\textrm{Probability of Edge Removal}$")
savefig("plots/cost_percolation.pdf")
savefig("plots/cost_percolation.png")

##### Plot path_percolation #####
plot(x, P0, ylims=(0,1), seriestype=:scatter, yerr = P0e, label=L"$P_0$", legend= :right,
guidefontsize=14, tickfontsize=12, legendfontsize=8, fontfamily="computer modern",
color_palette = palette(:plasma, 5), markersize=5)
yaxis!("Path probability")
plot!(x, P1, seriestype=:scatter, yerr = P1e, label=L"$P_1$", markersize=5)
plot!(x, P2, seriestype=:scatter, yerr = P2e, label=L"$P_2$", markersize=5)
plot!(x, P3, seriestype=:scatter, yerr = P3e, label=L"$P_3$", markersize=5)
plot!(x, P4, seriestype=:scatter, yerr = P3e, label=L"$P_4$", markersize=5)
xaxis!(L"$\textrm{Probability of Edge Removal}$")
savefig("plots/path_percolation.pdf")
savefig("plots/path_percolation.png")

# Plot ave_path_distance
# println(ave_path_distance)
# println(ave_path_distance_err)
# f = dB_to_Z.(ave_path_distance)
# e = dB_to_P.(ave_path_distance)
# plot(x, f, seriestype = :scatter, label="f", legend= true)
# plot!(x, e, seriestype = :scatter, label="e", legend= true)
# plot!(x, loss, ylims=(0,1), seriestype = :scatter, yerror = loss_err, label=L"$\eta$")
# plot!(x, z, seriestype = :scatter, yerror = z_err, label=L"$F$")

##### Plot path_distance_percolation #####
plot(x, ave_path_distance, yerr = ave_path_distance_err, seriestype = :scatter, label="Average Manhattan distance", legend = :bottomleft,
ylims = (0,10), right_margin = 15Plots.mm, guidefontsize=14, tickfontsize=12,
legendfontsize=8, fontfamily="computer modern", markersize=5, color=:orange)
ylabel!("Path distance")
plot!(x, postman, yerr=postman_err, seriestype = :scatter, label="Average shortest path distance",
markersize=5, color=:magenta)

# Plot ratio of postman with man
pathratio =  postman ./ ave_path_distance
pathratio_err = ((ave_path_distance_err ./ ave_path_distance) + (postman_err ./ postman)) .* pathratio
plot!(twinx(), x, pathratio, yerr = pathratio_err, seriestype = :scatter, color= :pink,
ylabel="Ratio", label="Path ratio", legend = :bottomright, guidefontsize=14, xticks=false, tickfontsize=12,
ylims = (0,1.6), markersize=4, markershape=:square)
# ylabel!("Ratio")

# Plot horizontal line for ave manhattan distance for 10x10 grid:
mandist = ones(length(x)) .* 6.67
plot!(x, mandist, linewidth=2, label="Analytic Manhattan distance")
savefig("plots/path_distance_percolation.pdf")
savefig("plots/path_distance_percolation.png")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
module TestConfigs

using Test
using GarishPrint
using Configurations

struct OptionA
    a::Int
    b::Int
end

struct OptionB
    a::Vector{Int}
    b::OptionA
end

Base.show(io::IO, x::OptionA) = GarishPrint.pprint_struct(io, x)
Base.show(io::IO, x::OptionB) = GarishPrint.pprint_struct(io, x)

opt = OptionB([1, 2, 3], OptionA(2, 2))
pprint_struct(opt)

@option "geometry" struct Geometry
    L::Int
    graph::Int
    radius::Float64 = 1.5
end

@option "emulation" struct Emulation
    pulse::String
    algo::String # = "Vern8"
    precision::String # = "float32"
    reltol::Maybe{Float64} = nothing
    abstol::Maybe{Float64} = nothing
    dt::Maybe{Float64} = nothing
    total_time::Vector{Int} = collect(1:7)
    postprocess::Bool = false
    gen_subspace::Bool = false
    geometry::Geometry
end

@option "postprocess" struct PostProcess
    pulse::String
    algo::String
    precision::String
    total_time::Vector{Int} = collect(1:7)
    gen_mappings::Bool=false
    skip_done::Bool=false
    geometry::Geometry
end

@option "cache" struct Cache
    subspace::Maybe{Geometry} = nothing
    postprocess::Maybe{Geometry} = nothing
end

@option struct JobInstance
    project::String
    device::Maybe{Int} = nothing
    emulation::Vector{Emulation} = Emulation[]
    postprocess::Vector{PostProcess} = PostProcess[]
    cache::Vector{Cache} = Cache[]
end

Base.show(io::IO, x::JobInstance) = pprint_struct(io, x; include_defaults=get(io, :include_defaults, false))
Base.show(io::IO, x::Cache) = pprint_struct(io, x; include_defaults=false)
Base.show(io::IO, x::PostProcess) = pprint_struct(io, x; include_defaults=false)
Base.show(io::IO, x::Emulation) = pprint_struct(io, x; include_defaults=false)
Base.show(io::IO, x::Geometry) = pprint_struct(io, x; include_defaults=false)

opt = JobInstance(
    project="test",
    emulation=[
        Emulation(;pulse="linear", algo="Vern8", precision="float32", geometry=Geometry(8, 188, 1.5)),
        Emulation(;pulse="linear", algo="Vern8", precision="float32", geometry=Geometry(8, 188, 1.5)),
    ],
    postprocess=[
        PostProcess(pulse="linear", algo="Vern8", precision="float32", geometry=Geometry(8, 188, 1.5)),
        PostProcess(pulse="linear", algo="Vern8", precision="float32", geometry=Geometry(8, 188, 1.5)),
    ]
)

show(stdout, MIME"text/plain"(), opt)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -24.652335523355834 <= q <= 21.32472335033785)

                     add_NL_constraint(m, :(log(1 + exp(-0.6982139593401175 + -0.19504968807074974*log(1 + exp(0.6793369023557885 + 0.6523875537646115*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.5142038874540145*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.0456439541419269*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6096854328317232*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + -0.41803298572193937*log(1 + exp(-0.7048442383777926 + 0.47781322839334983*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.713369260368808*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.44964139924772173*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.9876898123617992*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + 0.9542284451941447*log(1 + exp(0.9969387928417852 + -0.011395717836928387*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + 0.8103499615759371*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.12667586987008894*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6252342222679741*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + 0.7916680056214345*log(1 + exp(-0.5156629799588863 + -0.22745979325095078*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.8222624364306839*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.4937419854366709*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.1046356011041838*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + log(1 + exp(-0.4539405071565561 + 0.5320780591355936*log(1 + exp(0.6793369023557885 + 0.6523875537646115*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.5142038874540145*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.0456439541419269*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6096854328317232*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + -0.8489788623944037*log(1 + exp(-0.7048442383777926 + 0.47781322839334983*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.713369260368808*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.44964139924772173*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.9876898123617992*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + -0.25620536200506905*log(1 + exp(0.9969387928417852 + -0.011395717836928387*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + 0.8103499615759371*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.12667586987008894*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6252342222679741*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + 0.9264345634068367*log(1 + exp(-0.5156629799588863 + -0.22745979325095078*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.8222624364306839*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.4937419854366709*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.1046356011041838*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + log(1 + exp(-0.4017849565754932 + 0.7092320156791541*log(1 + exp(0.6793369023557885 + 0.6523875537646115*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.5142038874540145*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.0456439541419269*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6096854328317232*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + -0.20657569919094643*log(1 + exp(-0.7048442383777926 + 0.47781322839334983*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.713369260368808*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.44964139924772173*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.9876898123617992*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + 0.6260976974386958*log(1 + exp(0.9969387928417852 + -0.011395717836928387*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + 0.8103499615759371*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.12667586987008894*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6252342222679741*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + 0.026773003462286393*log(1 + exp(-0.5156629799588863 + -0.22745979325095078*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.8222624364306839*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.4937419854366709*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.1046356011041838*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) + log(1 + exp(0.8541923459688499 + -0.7979821016654509*log(1 + exp(0.6793369023557885 + 0.6523875537646115*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.5142038874540145*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.0456439541419269*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6096854328317232*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + 0.2172669330681649*log(1 + exp(-0.7048442383777926 + 0.47781322839334983*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.713369260368808*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.44964139924772173*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.9876898123617992*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + -0.6626207135044284*log(1 + exp(0.9969387928417852 + -0.011395717836928387*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + 0.8103499615759371*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.12667586987008894*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + 0.6252342222679741*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))) + -0.5735396155482761*log(1 + exp(-0.5156629799588863 + -0.22745979325095078*log(1 + exp(0.1880507757579264 + -0.8703531991881865*$(x[1]) + 0.8423062039785143*$(x[2]) + 0.27626604928751597*$(x[3]))) + -0.8222624364306839*log(1 + exp(0.8004597365811503 + 0.4093664270553261*$(x[1]) + -0.7983161406747903*$(x[2]) + -0.755209029087681*$(x[3]))) + -0.4937419854366709*log(1 + exp(-0.8428376848805481 + -0.4971126206674761*$(x[1]) + -0.0245959153233013*$(x[2]) + 0.4386350201089879*$(x[3]))) + -0.1046356011041838*log(1 + exp(-0.26516744777719925 + 0.12012704290515952*$(x[1]) + -0.4833932117044406*$(x[2]) + -0.25804110713494666*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#= Given a number n, print its prime factorization. We could do it
normal way which will be having time complexity of O (sqrt(N))
but we will be using Sieve to lessen the time complexity
We will make a sieve to get all the prime numbers and then get the
prime factorization of the number using that.=#

## Function

function PrimeFactor(n)
    a = zeros(Int64, n)
    for i = 2:n
        a[i] = i
    end
    for i = 2:n
        if (a[i] == i)
            for j = (i*i):i:n
                if (a[j] == j)
                    a[j] = i
                end
            end
        end
    end
    while (n != 1)
        println(a[n])
        n = div(n, a[n])
    end
end

## Input

n = readline()
n = parse(Int64, n)

#Calling the function

PrimeFactor(n)

#=
Sample Test case
Input:
    n = 495
Output:
    3 3 5 11

Time Complexity: O( log N )
=#


text:
using Test

@info "CANalyze.Frames tests..."
@testset "equal" begin
    using CANalyze.Frames

    @testset "equal_1" begin
        frame1 = CANFrame(0x14, Integer[]; is_extended=true)
        frame2 = CANFrame(0x14; is_extended=true)
        @test frame1 == frame2
    end

    @testset "equal_2" begin
        frame1 = CANFrame(0x14, Integer[]; is_extended=true)
        frame2 = CANFrame(0x15, Integer[]; is_extended=true)
        @test frame1 != frame2
    end

    @testset "equal_3" begin
        frame1 = CANFrame(0x14, Integer[]; is_extended=false)
        frame2 = CANFrame(0x14, Integer[]; is_extended=true)
        @test frame1 != frame2
    end

    @testset "equal_4" begin
        frame1 = CANFrame(0x14, Integer[]; is_extended=true)
        frame2 = CANFrame(0x15, Integer[]; is_extended=true)
        @test frame1 != frame2
    end

    @testset "equal_5" begin
        frame1 = CANFrame(0x14, Integer[1,2,3,4]; is_extended=true)
        frame2 = CANFrame(0x14, 1, 2, 3, 4; is_extended=true)
        @test frame1 == frame2
    end
end

@testset "frame_id" begin
    using CANalyze.Frames

    @testset "frame_id_1" begin
        frame = CANFrame(0x0AFF, Integer[1,2,3,4]; is_extended=true)
        @test frame_id(frame) == (0x0AFF & 0x01_FF_FF_FF)
    end

    @testset "frame_id_1" begin
        frame = CANFrame(0x0AFF, Integer[1,2,3,4]; is_extended=false)
        @test frame_id(frame) == (0x0AFF & 0x7FF)
    end
end

@testset "data" begin
    using CANalyze.Frames

    for i=0:8
        frame = CANFrame(0x0AFF, Integer[j for j=1:i]; is_extended=true)
        @test data(frame) == UInt8[j for j=1:i]
    end
end

@testset "dlc" begin
    using CANalyze.Frames

    for i=0:8
        frame = CANFrame(0x0AFF, Integer[j for j=1:i]; is_extended=true)
        @test dlc(frame) == i
    end
end

@testset "is_extended" begin
    using CANalyze.Frames

    @testset "is_extended_1" begin
        frame = CANFrame(0x0AFF; is_extended=true)
        @test is_extended(frame) == true
        @test is_standard(frame) == false
    end

    @testset "is_extended_2" begin
        frame = CANFrame(0x0AFF; is_extended=false)
        @test is_extended(frame) == false
        @test is_standard(frame) == true
    end
end

@testset "max_size" begin
    using CANalyze.Frames

    @testset "max_size_1" begin
        frame = CANFrame(0x0AFF; is_extended=true)
        @test max_size(typeof(frame)) == 8
    end

    @testset "max_size_2" begin
        frame = CANFdFrame(0x0AFF; is_extended=true)
        @test max_size(typeof(frame)) == 64
    end
end

@testset "too_much_data" begin
    using CANalyze.Frames

    @testset "too_much_data_1" begin
        @test_throws DomainError CANFrame(0x0AFF, [i for i=1:9]; is_extended=true)
    end

    @testset "too_much_data_2" begin
        @test_throws DomainError CANFdFrame(0x0AFF, [i for i=1:65]; is_extended=true)
    end
end


text:
using Base.Filesystem

function getDiskSpace( diskName  )
    cmdstr =  `cmd /c dir $diskName`
    # println( cmds )
    aa=read( cmdstr,String)
    bb= split( aa,'\n')
    for b1 in ( bb[end-1],)
        c= split(strip(b1), " ");
        for c1 in (c[3],)
            a11= c1
            # println( a11 )

            c0 = split(a11,",")
            w=""
            for c11 in c0
                w = string(w,c11)
            end
            # println( w ,"  ,    ",a11)
            space=  parse(UInt, w)
            return space
        end
    end
end
function getSizes( root,RAW1MaxSpace=3 , RAW2MaxSpace=3, RAW3MaxSpace=1, RAW4MaxSpace=1,Unit= 1024*1024*1024*1024)
    row1 = 0x0000
    row2= 0x0000
    row3 = 0x0000
    row4 = 0x0000
    row1fs=[]
    row2fs=[]
    row3fs=[]
    row4fs=[]
    a = Filesystem.readdir(root)
    for b in a
        m,e = splitext( b  )
        fsize = Base.Filesystem.filesize( joinpath(root,b) )
        # fsize =
        # println(m ," === ",e ,"  ", fsize)
        if occursin( "1",e)
            row1 += fsize
            append!(row1fs ,[b])
        elseif occursin( "2",e)
            row2 += fsize
            append!(row2fs ,[b])
        elseif occursin( "3",e)
            append!(row3fs ,[b])
            row3 += fsize
        elseif occursin( "4",e)
            append!(row4fs ,[b])
            row4 += fsize
        end
    end
    println("---------------------------")

    for (row,fs,max ) in ( (row1,row1fs,RAW1MaxSpace),(row2,row2fs,RAW2MaxSpace),(row3,row3fs,RAW3MaxSpace) ,(row4,row4fs,RAW4MaxSpace)   )
        sort!( fs,rev=true )
        # row *= Unit #in kbytes
        max *= Unit # in kbytes
        while row > max && length( fs ) >2
            if length( fs ) >1
                fs1 = pop!(fs)
                fullnamefs = joinpath(root,fs1)
                fsize = Base.Filesystem.filesize( fullnamefs )
                # fsize =round(UInt, s  )
                Filesystem.rm(fullnamefs )
                row -= fsize

                println( "$fullnamefs, Cursize(kbytes):$(row/Unit) ,Max(kbytes):$(max/Unit)")
            end
        end
        println()
    end
    # return ( row1,row2,row3,row4 )
end

function testMain(root)
    RAW1MaxSpace=0
    RAW2MaxSpace=0
    RAW3MaxSpace=0
    RAW4MaxSpace=1
    Unit= 1024*1024*1024
    while true
        rowsizes  = getSizes( root ,RAW1MaxSpace, RAW2MaxSpace, RAW3MaxSpace, RAW4MaxSpace,Unit)
        exit(1)
        sleep(10)
        println( "start clear expired files ")
    end
end

# testMain()

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
## Note:
## This is the `ChebyshevT.jl` test suite from the `Polynomials.jl` package by Miles Lucas

@testset "Construction" for coeff in [
    Int64[1, 1, 1, 1],
    Float32[1, -4, 2],
    ComplexF64[1 - 1im, 2 + 3im],
    [3 // 4, -2 // 1, 1 // 1]
]
    p = Chebyshev(coeff)
    @test p.coeffs == coeff
    @test coeffs(p) == coeff
    @test degree(p) == length(coeff) - 1
    @test Polynomials.indeterminate(p) == :x
    @test length(p) == length(coeff)
    @test size(p) == size(coeff)
    @test size(p, 1) == size(coeff, 1)
    @test typeof(p).parameters[1] == eltype(coeff)
    @test eltype(p) == eltype(coeff)
end

@testset "Other Construction" begin
    # Leading 0s
    p = Chebyshev([1, 2, 0, 0])
    @test p.coeffs == [1, 2]
    @test length(p) == 2

    # Different type
    p = Chebyshev{Float64}(ones(Int32, 4))
    @test p.coeffs == ones(Float64, 4)

    p = Chebyshev(30)
    @test p.coeffs == [30]

    p = zero(Chebyshev{Int})
    @test p.coeffs == Int[]

    p = one(Chebyshev{Int})
    @test p.coeffs == [1]

    pNULL = Chebyshev(Int[])
    @test iszero(pNULL)
    @test degree(pNULL) == -1

    p0 = Chebyshev([0])
    @test iszero(p0)
    @test degree(p0) == -1
end

@testset "Roots $i" for i in 1:5
    roots = cos.(range(-π, 0, length = 2i + 1)[2:2:end])
    target = Chebyshev(vcat(zeros(i), 1))
    res = fromroots(ChebyshevT, roots) .* 2^(i - 1)
    @test res == target
end

@testset "Roots" begin
    r = [-1, 1, 0]
    c = fromroots(ChebyshevT, r)
    @test c == Chebyshev([0, -0.25, 0, 0.25])
    @test roots(c) ≈ r

    r = [1im, -1im]
    c = fromroots(ChebyshevT, r)
    @test c ≈ Chebyshev([1.5 + 0im, 0 + 0im, 0.5 + 0im])
    @test roots(c) ≈ r

end

@testset "Values" begin
    c1 = Chebyshev([2.5, 2.0, 1.5]) # 1 + 2x + 3x^2
    x = -1:0.5:1
    expected = [2.0, 0.75, 1.0, 2.75, 6.0]
    @test c1.(x) == expected
end

@testset "Conversions" begin
    c1 = Chebyshev([2.5, 2.0, 1.5])
    @test c1 ≈ Polynomial([1, 2, 3])
    p = convert(Polynomial{Float64}, c1)
    @test p == Polynomial{Float64}([1, 2, 3])
    @test convert(ChebyshevT, p) == c1

end

@testset "Companion" begin
    c_null = Chebyshev(Int[])
    c_1 = Chebyshev([1])
    @test_throws ArgumentError companion(c_null)
    @test_throws ArgumentError companion(c_1)
    for i in 1:5
        coef = vcat(zeros(i), 1)
        c = Chebyshev(coef)
        @test size(companion(c)) == (i, i)
    end
    c = Chebyshev([1, 2])
    @test companion(c)[1, 1] == -0.5
end

@testset "Vander" begin
    x = 0:0.1:1
    v = vander(ChebyshevT, x, 5)
    @test size(v) == (length(x), 6)
    @inbounds for i in 1:6
        coef = vcat(zeros(i - 1), 1)
        c = Chebyshev(coef)
        @test v[:, i] ≈ c.(x)
    end
end

@testset "Arithmetic $i, $j" for i in 0:5, j in 0:5
    # multiplication
    target = zeros(i + j + 1)
    target[end] += 0.5
    target[abs(i - j) + 1] += 0.5
    c1 = Chebyshev(vcat(zeros(i), 1))
    c2 = Chebyshev(vcat(zeros(j), 1))
    @test c1 * c2 ≈ Chebyshev(target)

    # divrem
    target = c1 + c2
    quo, rem = divrem(target, c1)
    res = quo * c1 + rem
    @test res ≈ target
end

@testset "Mapdomain" begin
    x = -30:20
    mx = mapdomain(ChebyshevT, x)
    @test extrema(mx) == (-1, 1)

    x = 0.5:0.01:0.6
    mx = mapdomain(ChebyshevT, x)
    @test extrema(mx) == (-1, 1)
end

@testset "Arithmetic" begin
    # multiplication
    c1 = Chebyshev([1, 2, 3])
    c2 = Chebyshev([3, 2, 1])
    @test c1 * c2 == Chebyshev([6.5, 12, 12, 4, 1.5])

    c1 = Chebyshev([1, 2, 3])
    c2 = Chebyshev([3, 2, 1])
    d, r = divrem(c1, c2)
    @test d.coeffs ≈ [3]
    @test r.coeffs ≈ [-8, -4]

    c2 = Chebyshev([0, 1, 2, 3])
    d, r = divrem(c2, c1)

    @test d.coeffs ≈ [0, 2]
    @test r.coeffs ≈ [-2, -4]

    # GCD
    c1 = Chebyshev([1, 2, 3])
    c2 = Chebyshev([3, 2, 1])
    @test gcd(c1, c2) ≈ Chebyshev(6)
end

@testset "integrals derivatives" begin
    c1 = one(Chebyshev{Int})
    @test integrate(c1) == Chebyshev([0, 1])
    for k in [-3, 0, 2]
        @test integrate(c1, k) == Chebyshev([k, 1])
    end

    for i in 0:4
        scl = i + 1
        p = Polynomial(vcat(zeros(i), 1))
        target = Polynomial(vcat(i, zeros(i), 1 / scl))
        cheb = convert(ChebyshevT, p)
        cint = integrate(cheb, i)
        res = convert(Polynomial, cint)
        @test res ≈ target
        @test derivative(cint) == cheb
    end
end

@testset "z-series" for i in 0:5
    # c to z
    input = vcat(2, ones(i))
    target = vcat(fill(0.5, i), 2, fill(0.5, i))
    zs = Polynomials._c_to_z(input)
    @test zs == target
    c = Polynomials._z_to_c(zs)
    @test c == input

    # div
    z1 = [0.5, 0.0, 0.5]
    z2 = [0.5, 0.0, 0.5]
    q, r = Polynomials._z_division(z1, z2)
    @test q == [1]
    @test r == [0]
end

@testset "Fit" begin

    f(x) = exp(x)
    p = fit(Chebyshev{Float64}, f, 16)
    @test maximum(abs.(p(x) - f(x) for  x in range(0, 1, length=100))) <= 100*eps()

    f(x) = sin(6x) - sin(60*exp(x)) # chebfun example
    p = fit(Chebyshev{Float64}, f, 150)
    @test maximum(abs.(p(x) - f(x) for  x in range(0, 1, length=100))) <= 1000*eps()

end


text:
include("utils.jl") # contains most dependencies and the saving function
using PDMats
using FillArrays

s = ArgParseSettings()
@add_arg_table! s begin
	"dataset"
		default = "MNIST"
		arg_type = String
		help = "dataset"
end
parsed_args = parse_args(ARGS, s)
@unpack dataset = parsed_args

#######################################################################################
################ THIS PART IS TO BE PROVIDED FOR EACH MODEL SEPARATELY ################

modelname = "fAnoGAN"

"""
This returns encodings, parameters and scoring functions in order to reconstruct the experiment.
This is a slightly updated version of the original run script.
"""
function evaluate(model_data, data, parameters)
	# load the model file, extract params and model
	model = model_data["model"] |> gpu

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = get(model_data, "fit_t", nothing),
		history = get(model_data, "history", nothing),
		npars = get(model_data, "npars", nothing),
		model = model |> cpu
	)

	# now return the different scoring functions
	training_info, [(x -> GenerativeAD.Models.anomaly_score_gpu(model, x), parameters)]
end

##################
# this is common #
##################
include("run_loop.jl")


text:
module KernelAbstractions

export @kernel
export @Const, @localmem, @private, @uniform, @synchronize, @index, groupsize, @print
export Device, GPU, CPU, CUDA, Event, MultiEvent, NoneEvent
export async_copy!

using MacroTools
using StaticArrays
using Cassette
using Requires
using Adapt

"""
   @kernel function f(args) end

Takes a function definition and generates a Kernel constructor from it.
The enclosed function is allowed to contain kernel language constructs.
In order to call it the kernel has first to be specialized on the backend
and then invoked on the arguments.

# Kernel language

- [`@Const`](@ref)
- [`@index`](@ref)
- [`@localmem`](@ref)
- [`@private`](@ref)
- [`@uniform`](@ref)
- [`@synchronize`](@ref)
- [`@print`](@ref)

# Example:

@kernel function vecadd(A, @Const(B))
    I = @index(Global)
    @inbounds A[I] += B[I]
end

A = ones(1024)
B = rand(1024)
event = vecadd(CPU(), 64)(A, B, ndrange=size(A))
wait(event)
"""
macro kernel(expr)
    __kernel(expr)
end

"""
   @Const(A)

`@Const` is an argument annotiation that asserts that the memory reference
by `A` is both not written to as part of the kernel and that it does not alias
any other memory in the kernel.

!!! danger
    Violating those constraints will lead to arbitrary behaviour.

as an example given a kernel signature `kernel(A, @Const(B))`, you are not
allowed to call the kernel with `kernel(A, A)` or `kernel(A, view(A, :))`.
"""
macro Const end

abstract type Event end
import Base.wait

struct NoneEvent <: Event end

struct MultiEvent{T} <: Event
    events::T
    MultiEvent() = new{Tuple{}}(())
    function MultiEvent(events::Tuple{Vararg{<:Event}})
        evs = tuplejoin(map(flatten, events)...)
        new{typeof(evs)}(evs)
    end
    function MultiEvent(event::E) where {E<:Event}
        new{Tuple{E}}((event,))
    end
end
MultiEvent(::Nothing) = MultiEvent()
MultiEvent(ev::MultiEvent) = ev

@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

flatten(ev::MultiEvent) = tuplejoin(map(flatten, ev.events)...)
flatten(ev::NoneEvent) = ()
flatten(ev::Event) = (ev,)

"""
    async_copy!(::Device, dest::AbstractArray, src::AbstractArray; dependencies = nothing)

Perform an asynchronous copy on the device. Returns an event that can be waited upon.
"""
function async_copy! end

###
# Kernel language
# - @localmem
# - @private
# - @uniform
# - @synchronize
# - @index
# - groupsize
###

"""
    groupsize()

Query the workgroupsize on the device. This function returns
a tuple corresponding to kernel configuration. In order to get
the total size you can use `prod(groupsize())`.
"""
function groupsize end

"""
    @localmem T dims
"""
macro localmem(T, dims)
    # Stay in sync with CUDAnative
    id = gensym("static_shmem")

    quote
        $SharedMemory($(esc(T)), Val($(esc(dims))), Val($(QuoteNode(id))))
    end
end

"""
    @private T dims
"""
macro private(T, dims)
    quote
        $Scratchpad($(esc(T)), Val($(esc(dims))))
    end
end

"""
    @uniform value
"""
macro uniform(value)
    esc(value)
end

"""
   @synchronize()

After a `@synchronize` statement all read and writes to global and local memory
from each thread in the workgroup are visible in from all other threads in the
workgroup.
"""
macro synchronize()
    quote
        $__synchronize()
    end
end

"""
   @synchronize(cond)

After a `@synchronize` statement all read and writes to global and local memory
from each thread in the workgroup are visible in from all other threads in the
workgroup. `cond` is not allowed to have any visible sideffects.

# Platform differences
  - `GPU`: This synchronization will only occur if the `cond` evaluates.
  - `CPU`: This synchronization will always occur.
"""
macro synchronize(cond)
    quote
        $(esc(cond)) && $__synchronize()
    end
end

"""
   @print(items...)

This is a unified print statement.

# Platform differences
  - `GPU`: This will reorganize the items to print via @cuprintf
  - `CPU`: This will call `print(items...)`
"""
macro print(items...)

    args = Union{Val,Expr,Symbol}[]

    items = [items...]
    while true
        isempty(items) && break

        item = popfirst!(items)

        # handle string interpolation
        if isa(item, Expr) && item.head == :string
            items = vcat(item.args, items)
            continue
        end

        # expose literals to the generator by using Val types
        if isbits(item) # literal numbers, etc
            push!(args, Val(item))
        elseif isa(item, QuoteNode) # literal symbols
            push!(args, Val(item.value))
        elseif isa(item, String) # literal strings need to be interned
            push!(args, Val(Symbol(item)))
        else # actual values that will be passed to printf
            push!(args, item)
        end
    end

    quote
        $__print($(map(esc,args)...))
    end
end

"""
   @index

The `@index` macro can be used to give you the index of a workitem within a kernel
function. It supports both the production of a linear index or a cartesian index.
A cartesian index is a general N-dimensional index that is derived from the iteration space.

# Index granularity

  - `Global`: Used to access global memory.
  - `Group`: The index of the `workgroup`.
  - `Local`: The within `workgroup` index.

# Index kind

  - `Linear`: Produces an `Int64` that can be used to linearly index into memory.
  - `Cartesian`: Produces a `CartesianIndex{N}` that can be used to index into memory.
  - `NTuple`: Produces a `NTuple{N}` that can be used to index into memory.

If the index kind is not provided it defaults to `Linear`, this is suspect to change.

# Examples

```julia
@index(Global, Linear)
@index(Global, Cartesian)
@index(Local, Cartesian)
@index(Group, Linear)
@index(Local, NTuple)
@index(Global)
```
"""
macro index(locale, args...)
    if !(locale === :Global || locale === :Local || locale === :Group)
        error("@index requires as first argument either :Global, :Local or :Group")
    end

    if length(args) >= 1
        if args[1] === :Cartesian ||
           args[1] === :Linear ||
           args[1] === :NTuple
            indexkind = args[1]
            args = args[2:end]
        else
            indexkind = :Linear
        end
    else
        indexkind = :Linear
    end

    index_function = Symbol(:__index_, locale, :_, indexkind)
    Expr(:call, GlobalRef(KernelAbstractions, index_function), map(esc, args)...)
end

###
# Internal kernel functions
###

function __index_Local_Linear end
function __index_Group_Linear end
function __index_Global_Linear end

function __index_Local_Cartesian end
function __index_Group_Cartesian end
function __index_Global_Cartesian end

__index_Local_NTuple(I...) = Tuple(__index_Local_Cartesian(I...))
__index_Group_NTuple(I...) = Tuple(__index_Group_Cartesian(I...))
__index_Global_NTuple(I...) = Tuple(__index_Global_Cartesian(I...))

struct ConstAdaptor end

Adapt.adapt_storage(to::ConstAdaptor, a::Array) = Base.Experimental.Const(a)

constify(arg) = adapt(ConstAdaptor(), arg)

###
# Backend hierarchy
###

abstract type Device end
abstract type GPU <: Device end

struct CPU <: Device end
struct CUDA <: GPU end
# struct AMD <: GPU end
# struct Intel <: GPU end

include("nditeration.jl")
using .NDIteration
import .NDIteration: get

###
# Kernel closure struct
###

"""
    Kernel{Device, WorkgroupSize, NDRange, Func}

Kernel closure struct that is used to represent the device
kernel on the host. `WorkgroupSize` is the number of workitems
in a workgroup.
"""
struct Kernel{Device, WorkgroupSize<:_Size, NDRange<:_Size, Fun}
    f::Fun
end

workgroupsize(::Kernel{D, WorkgroupSize}) where {D, WorkgroupSize} = WorkgroupSize
ndrange(::Kernel{D, WorkgroupSize, NDRange}) where {D, WorkgroupSize,NDRange} = NDRange

function partition(kernel, ndrange, workgroupsize)
    static_ndrange = KernelAbstractions.ndrange(kernel)
    static_workgroupsize = KernelAbstractions.workgroupsize(kernel)

    if ndrange === nothing && static_ndrange <: DynamicSize ||
       workgroupsize === nothing && static_workgroupsize <: DynamicSize
        errmsg = """
            Can not partition kernel!

            You created a dynamically sized kernel, but forgot to provide runtime
            parameters for the kernel. Either provide them statically if known
            or dynamically.
            NDRange(Static):  $(static_ndrange)
            NDRange(Dynamic): $(ndrange)
            Workgroupsize(Static):  $(static_workgroupsize)
            Workgroupsize(Dynamic): $(workgroupsize)
        """
        error(errmsg)
    end

    if static_ndrange <: StaticSize
        if ndrange !== nothing && ndrange != get(static_ndrange)
            error("Static NDRange ($static_ndrange) and launch NDRange ($ndrange) differ")
        end
        ndrange = get(static_ndrange)
    end

    if static_workgroupsize <: StaticSize
        if workgroupsize !== nothing && workgroupsize != get(static_workgroupsize)
            error("Static WorkgroupSize ($static_workgroupsize) and launch WorkgroupSize $(workgroupsize) differ")
        end
        workgroupsize = get(static_workgroupsize)
    end

    @assert workgroupsize !== nothing
    @assert ndrange !== nothing
    blocks, workgroupsize, dynamic = NDIteration.partition(ndrange, workgroupsize)

    if static_ndrange <: StaticSize
        static_blocks = StaticSize{blocks}
        blocks = nothing
    else
        static_blocks = DynamicSize
        blocks = CartesianIndices(blocks)
    end

    if static_workgroupsize <: StaticSize
        static_workgroupsize = StaticSize{workgroupsize} # we might have padded workgroupsize
        workgroupsize = nothing
    else
        workgroupsize = CartesianIndices(workgroupsize)
    end

    iterspace = NDRange{length(ndrange), static_blocks, static_workgroupsize}(blocks, workgroupsize)
    return iterspace, dynamic
end

###
# Compiler/Cassette
###

include("compiler.jl")

###
# Compiler/Frontend
###

function __workitems_iterspace end
function __validindex end

include("macros.jl")

###
# Backends/Interface
###

function Scratchpad(::Type{T}, ::Val{Dims}) where {T, Dims}
    throw(MethodError(Scratchpad, (T, Val(Dims))))
end

function SharedMemory(::Type{T}, ::Val{Dims}, ::Val{Id}) where {T, Dims, Id}
    throw(MethodError(SharedMemory, (T, Val(Dims), Val(Id))))
end

function __synchronize()
    error("@synchronize used outside kernel or not captured")
end

@generated function __print(items...)
    str = ""
    args = []

    for i in 1:length(items)
        item = :(items[$i])
        T = items[i]
        if T <: Val
            item = QuoteNode(T.parameters[1])
        end
        push!(args, item)
    end

    quote
        print($(args...))
    end
end

###
# Backends/Implementation
###

# Utils
__size(args::Tuple) = Tuple{args...}
__size(i::Int) = Tuple{i}

include("backends/cpu.jl")
@init @require CUDAnative="be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    include("backends/cuda.jl")
end

###
# Extras
# - LoopInfo
###

include("extras/extras.jl")
end #module