text:
__precompile__()

"""
The Luxor package provides a set of vector drawing functions for creating graphical documents.
"""
module Luxor

using Colors, Cairo, Compat

include("point.jl")
include("Turtle.jl")
include("polygons.jl")
include("Tiler.jl")
include("arrows.jl")
include("text.jl")

export Drawing, currentdrawing,
    rescale,

    finish, preview,
    origin, axes, background,
    newpath, closepath, newsubpath,
    circle, ellipse, squircle, center3pts,
    rect, box, setantialias, setline, setlinecap, setlinejoin, setdash,
    move, rmove,
    line, rline, curve, arc, carc, arc2r, ngon, sector, pie,
    do_action, stroke, fill, paint, paint_with_alpha, fillstroke,

    poly, simplify, polybbox, polycentroid, polysortbyangle, polysortbydistance, midpoint,
    prettypoly,

    star,

    intersection, polysplit,

    strokepreserve, fillpreserve,
    gsave, grestore,
    scale, rotate, translate,
    clip, clippreserve, clipreset,

    isinside,

    getpath, getpathflat,

    pattern_create_radial, pattern_create_linear,
    pattern_add_color_stop_rgb, pattern_add_color_stop_rgba,
    pattern_set_filter, pattern_set_extend,

    fontface, fontsize, text, textpath,
    textextents, textcurve, textcentred, textcentered, textright,
    setcolor, setopacity, sethue, randomhue, randomcolor, @setcolor_str,
    getmatrix, setmatrix, transform,

    readpng, placeimage,
    Tiler,
    arrow

# as of version 0.4, we still share fill() and scale() with Base.

import Base: fill, scale

type Drawing
    width::Float64
    height::Float64
    filename::String
    surface::CairoSurface
    cr::CairoContext
    surfacetype::String
    redvalue::Float64
    greenvalue::Float64
    bluevalue::Float64
    alpha::Float64
    function Drawing(w=800, h=800, f="luxor-drawing.png")
        global currentdrawing
        (path, ext)         = splitext(f)
        if ext == ".pdf"
            the_surface     =  Cairo.CairoPDFSurface(f, w, h)
            the_surfacetype = "pdf"
            the_cr          =  Cairo.CairoContext(the_surface)
        elseif ext == ".png" || ext == "" # default to PNG
            the_surface     = Cairo.CairoARGBSurface(w,h)
            the_surfacetype = "png"
            the_cr          = Cairo.CairoContext(the_surface)
        elseif ext == ".eps"
            the_surface     = Cairo.CairoEPSSurface(f, w,h)
            the_surfacetype = "eps"
            the_cr          = Cairo.CairoContext(the_surface)
        elseif ext == ".svg"
            the_surface     = Cairo.CairoSVGSurface(f, w,h)
            the_surfacetype = "svg"
            the_cr          = Cairo.CairoContext(the_surface)
        end
        currentdrawing      = new(w, h, f, the_surface, the_cr, the_surfacetype, 0, 0, 0, 1)
        return "drawing '$f' ($w w x $h h) created in $(pwd())"
    end
end

"""
The `paper_sizes` Dictionary holds a few paper sizes, width is first, so default is Portrait:

```
"A0"      => (2384, 3370),
"A1"      => (1684, 2384),
"A2"      => (1191, 1684),
"A3"      => (842, 1191),
"A4"      => (595, 842),
"A5"      => (420, 595),
"A6"      => (298, 420),
"A"       => (612, 792),
"Letter"  => (612, 792),
"Legal"   => (612, 1008),
"Ledger"  => (792, 1224),
"B"       => (612, 1008),
"C"       => (1584, 1224),
"D"       => (2448, 1584),
"E"       => (3168, 2448))
```
"""
paper_sizes = Dict{String, Tuple}(
  "A0" => (2384, 3370),
  "A1" => (1684, 2384),
  "A2" => (1191, 1684),
  "A3" => (842, 1191),
  "A4" => (595, 842),
  "A5" => (420, 595),
  "A6" => (298, 420),
  "A" => (612, 792),
  "Letter" => (612, 792),
  "Legal" => (612, 1008),
  "Ledger" => (792, 1224),
  "B" => (612, 1008),
  "C" => (1584, 1224),
  "D" => (2448, 1584),
  "E" => (3168, 2448))

"""
Create a new drawing, and optionally specify file type (PNG, PDF, SVG, etc) and dimensions.

    Drawing()

creates a drawing, defaulting to PNG format, default filename "luxor-drawing.png",
default size 800 pixels square.

You can specify dimensions, and use the default target filename:

    Drawing(300,300)

creates a drawing 300 by 300 pixels, defaulting to PNG format, default filename
"/tmp/luxor-drawing.png".

    Drawing(300,300, "my-drawing.pdf")

creates a PDF drawing in the file "my-drawing.pdf", 300 by 300 pixels.

    Drawing(800,800, "my-drawing.svg")`

creates an SVG drawing in the file "my-drawing.svg", 800 by 800 pixels.

    Drawing(800,800, "my-drawing.eps")

creates an EPS drawing in the file "my-drawing.eps", 800 by 800 pixels.

    Drawing("A4", "my-drawing.pdf")

creates a drawing in ISO A4 size in the file "my-drawing.pdf". Other sizes available
are:  "A0", "A1", "A2", "A3", "A4", "A5", "A6", "Letter", "Legal", "A", "B", "C", "D", "E".
Append "landscape" to get the landscape version.

    Drawing("A4landscape")

creates the drawing A4 landscape size.

PDF files default to a white background, but PNG defaults to transparent, unless you specify
one using `background()`.
"""
function Drawing(paper_size::String, f="luxor-drawing.png")
  if contains(paper_size, "landscape")
    psize = replace(paper_size, "landscape", "")
    h, w = paper_sizes[psize]
  else
    w, h = paper_sizes[paper_size]
  end
  Drawing(w, h, f)
end

"""
    finish()

Finish the drawing, and close the file. You may be able to open it in an external viewer
application with `preview()`.
"""
function finish()
    if currentdrawing.surfacetype == "png"
        Cairo.write_to_png(currentdrawing.surface, currentdrawing.filename)
        Cairo.finish(currentdrawing.surface)
        Cairo.destroy(currentdrawing.surface)
    else
        Cairo.finish(currentdrawing.surface)
        Cairo.destroy(currentdrawing.surface)
    end
end

"""
    preview()

On macOS, open the file, which probably uses the default, Preview.app.
On Unix, open the file with `xdg-open`.
On Windows, pass the filename to the shell.
"""
function preview()
    if @compat is_apple()
        run(`open $(currentdrawing.filename)`)
    elseif @compat is_windows()
        run(`$(currentdrawing.filename)`)
    elseif @compat is_unix()
        run(`xdg-open $(currentdrawing.filename)`)
    end
end

"""
    origin()

Reset the current matrix, and then set the 0/0 origin to the center of the drawing (otherwise
it will stay at the top left corner, the default).

You can refer to the 0/0 point as `O`. (O = `Point(0, 0)`),
"""
function origin()
    setmatrix([1.0, 0.0, 0.0, 1.0, 0.0, 0.0])
    Cairo.translate(currentdrawing.cr, currentdrawing.width/2., currentdrawing.height/2.)
end

"""
Convert or rescale a value between `oldmin`/`oldmax` to the equivalent value between
`newmin`/`newmax`.

    rescale(value, oldmin, oldmax, newmin, newmax)

For example, to convert 42 lying on a scale from 0 and 100 to the equivalent number
between 1 and 0 (inverting the direction):

    rescale(42, 0, 100, 1, 0)

returns 0.5800000000000001
"""
rescale(value, oldmin, oldmax, newmin=0, newmax=1) =
   ((value - oldmin) / (oldmax - oldmin)) * (newmax - newmin) + newmin

"""
Draw two axes lines starting at `O`, the current 0/0, and continuing out along the current
positive x and y axes.
"""
function axes()
    # draw axes
    gsave()
    setline(1)
    fontsize(20)
    sethue("gray")
    arrow(O, Point(currentdrawing.width/2. - 30, 0))
    text("x", Point(currentdrawing.width/2. - 30, -15))
    text("x", Point(30, -15))
    arrow(O, Point(0, currentdrawing.height/2. - 30))
    text("y", Point(5, currentdrawing.width/2. - 30))
    text("y", Point(5, 30))
    grestore()
end

"""
    background(color)

Fill the canvas with a single color. Returns the (red, green, blue, alpha) values.

Examples:

    background("antiquewhite")
    background("ivory")
    background(Colors.RGB(0, 0, 0))
    background(Colors.Luv(20, -20, 30))

If you don't specify a background color for a PNG drawing, the background will be
transparent.  You can set a partial of completely transparent background for PNG files by
passing a color with an alpha value, such as this transparent black:

    background(RGBA(0, 0, 0, 0))

"""
function background(col::String)
   setcolor(col)
   paint()
   return (currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, currentdrawing.alpha)
end

function background(col::ColorTypes.Colorant)
    temp = convert(RGBA,  col)
    setcolor(temp.r, temp.g, temp.b, temp.alpha)
    paint()
    return (currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, currentdrawing.alpha)
end

function background(r, g, b)
    sethue(r, g, b)
    paint()
    return (currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, currentdrawing.alpha)
end

function background(r, g, b, a)
    setcolor(r, g, b, a)
    paint()
    return (currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, currentdrawing.alpha)
end

# does this do anything in Cairo?
setantialias(n) = Cairo.set_antialias(currentdrawing.cr, n)

"""
    newpath()

Create a new path. This is Cairo's `new_path()` function.
"""
newpath() = Cairo.new_path(currentdrawing.cr)

"""
    newsubpath()

Add a new subpath to the current path. This is Cairo's `new_sub_path()` function. It can
be used for example to make holes in shapes.
"""
newsubpath() = Cairo.new_sub_path(currentdrawing.cr)

"""
    closepath()

Close the current path. This is Cairo's `close_path()` function.
"""
closepath() = Cairo.close_path(currentdrawing.cr)

"""
Stroke the current path with the current line width, line join, line cap, and dash settings.
The current path is then cleared.

    stroke()

"""
stroke() = Cairo.stroke(currentdrawing.cr)

"""
Fill the current path with current settings. The current path is then cleared.

    fill()

"""
fill() = Cairo.fill(currentdrawing.cr)

"""
Paint the current clip region with the current settings.

    paint()

"""
paint() = Cairo.paint(currentdrawing.cr)

"""
Stroke the current path with current line width, line join, line cap, and dash
settings, but then keep the path current.

    strokepreserve()

"""
strokepreserve()    = Cairo.stroke_preserve(currentdrawing.cr)

"""
Fill the current path with current settings, but then keep the path current.

    fillpreserve()

"""
fillpreserve()      = Cairo.fill_preserve(currentdrawing.cr)

"""
Fill and stroke the current path.
"""
function fillstroke()
    fillpreserve()
    stroke()
end

"""
    do_action(action)

This is usually called by other graphics functions. Actions for graphics commands include
`:fill`, `:stroke`, `:clip`, `:fillstroke`, `:fillpreserve`, `:strokepreserve` and `:path`.
"""
function do_action(action)
    if action == :fill
        fill()
    elseif action == :stroke
        stroke()
    elseif action == :clip
        clip()
    elseif action == :fillstroke
        fillstroke()
    elseif action == :fillpreserve
        fillpreserve()
    elseif action == :strokepreserve
        strokepreserve()
    elseif action == :none
    end
    return true
end

"""
Establish a new clipping region by intersecting the current clipping region with the
current path and then clearing the current path.

    clip()
"""
clip() = Cairo.clip(currentdrawing.cr)

"""
Establish a new clipping region by intersecting the current clipping region with the current
path, but keep the current path.

    clippreserve()
"""
clippreserve() = Cairo.clip_preserve(currentdrawing.cr)

"""
Reset the clipping region to the current drawing's extent.

    clipreset()
"""
clipreset() = Cairo.reset_clip(currentdrawing.cr)

"""
Make a circle of radius `r` centred at `x`/`y`.

    circle(x, y, r, action=:nothing)

`action` is one of the actions applied by `do_action`, defaulting to `:nothing`. You can
also use `ellipse()` to draw circles and place them by their centerpoint.
"""
function circle(x, y, r, action=:nothing)
    if action != :path
      newpath()
    end
    Cairo.circle(currentdrawing.cr, x, y, r)
    do_action(action)
end

"""
Make a circle centred at `pt`.

    circle(pt, r, action)

"""
circle(centerpoint::Point, r, action=:nothing) =
  circle(centerpoint.x, centerpoint.y, r, action)

"""
Make a circle that passes through two points that define the diameter:

    circle(pt1::Point, pt2::Point, action=:nothing)
"""
function circle(pt1::Point, pt2::Point, action=:nothing)
  center = midpoint(pt1, pt2)
  radius = norm(pt1, pt2)/2
  circle(center, radius, action)
end

"""
Find the radius and center point for three points lying on a circle.

    center3pts(a::Point, b::Point, c::Point)

returns (centerpoint, radius) of a circle. Then you can use `circle()` to place a
circle, or `arc()` to draw an arc passing through those points.

If there's no such circle, then you'll see an error message in the console and the function
returns `(Point(0,0), 0)`.
"""
function center3pts(a::Point, b::Point, c::Point)
# Find perpendicular bisectors of the segments connecting the first two and last two
# points. If they bisectors intersect, that's the center of the circle. If they don't,
# the points are colinear so they do not define a circle.

  # Get the perpendicular bisector of (x1, y1) and (x2, y2).
  x1 = (b.x + a.x) / 2
  y1 = (b.y + a.y) / 2
  dy1 = b.x - a.x
  dx1 = -(b.y - a.y)

  # Get the perpendicular bisector of (x2, y2) and (x3, y3).
  x2 = (c.x + b.x) / 2
  y2 = (c.y + b.y) / 2
  dy2 = c.x - b.x
  dx2 = -(c.y - b.y)

  # See where the lines intersect.

  dxy = (dx1 * dy2 - dy1 * dx2)
  if isapprox(dxy, 0)
    info("no circle passes through all three points")
    return Point(0,0), 0
  end

  ox = (y1 * dx1 * dx2 + x2 * dx1 * dy2 - x1 * dy1 * dx2 - y2 * dx1 * dx2) / dxy
  oy = (ox - x1) * dy1 / dx1 + y1
  dx = ox - a.x
  dy = oy - a.y
  radius = sqrt(dx * dx + dy * dy)
  return Point(ox, oy), radius
end

"""
Make an ellipse, centered at `xc/yc`, fitting in a box of width `w` and height `h`.

    ellipse(xc, yc, w, h, action=:none)
"""
function ellipse(xc, yc, w, h, action=:none)
    x  = xc - w/2
    y  = yc - h/2
    kappa = .5522848 # ??? http://www.whizkidtech.redprince.net/bezier/circle/kappa/
    ox = (w / 2) * kappa  # control point offset horizontal
    oy = (h / 2) * kappa  # control point offset vertical
    xe = x + w            # x-end
    ye = y + h            # y-end
    xm = x + w / 2        # x-middle
    ym = y + h / 2        # y-middle
    move(x, ym)
    curve(x, ym - oy, xm - ox, y, xm, y)
    curve(xm + ox, y, xe, ym - oy, xe, ym)
    curve(xe, ym + oy, xm + ox, ye, xm, ye)
    curve(xm - ox, ye, x, ym + oy, x, ym)
    do_action(action)
end

"""
Make an ellipse, centered at point `c`, with width `w`, and height `h`.

    ellipse(cpt, w, h, action=:none)
"""
ellipse(c::Point, w, h, action=:none) = ellipse(c.x, c.y, w, h, action)

"""
Make a squircle (basically a rectangle with rounded corners). Specify the center position, horizontal radius (distance from center to a side), and vertical radius (distance from center to top or bottom):

    squircle(center::Point, hradius, vradius, action=:none; rt = 0.5, vertices=false)

The `rt` option defaults to 0.5, and gives an intermediate shape. Values less than 0.5 make the shape more square. Values above make the shape more round.
"""

function squircle(center::Point, hradius, vradius, action=:none; rt = 0.5, vertices=false, reversepath=false)
  points = Point[]
  for theta in 0:pi/40:2pi
      xpos = center.x + ^(abs(cos(theta)), rt) * hradius * sign(cos(theta))
      ypos = center.y + ^(abs(sin(theta)), rt) * vradius * sign(sin(theta))
      push!(points, Point(xpos, ypos))
  end
  if vertices
    return points
  end
  poly(points, action, close=true, reversepath=reversepath)
end

"""
Add an arc to the current path from `angle1` to `angle2` going clockwise.

    arc(xc, yc, radius, angle1, angle2, action=:nothing)

Angles are defined relative to the x-axis, positive clockwise.

"""
function arc(xc, yc, radius, angle1, angle2, action=:nothing)
  Cairo.arc(currentdrawing.cr, xc, yc, radius, angle1, angle2)
  do_action(action)
end

"""
Add an arc to the current path from `angle1` to `angle2` going clockwise.

    arc(centerpoint::Point, radius, angle1, angle2, action=:nothing)
"""

arc(centerpoint::Point, radius, angle1, angle2, action=:nothing) =
  arc(centerpoint.x, centerpoint.y, radius, angle1, angle2, action)

"""
Add an arc to the current path from `angle1` to `angle2` going counterclockwise.

    carc(xc, yc, radius, angle1, angle2, action=:nothing)

Angles are defined relative to the x-axis, positive clockwise.

"""
function carc(xc, yc, radius, angle1, angle2, action=:nothing)
  Cairo.arc_negative(currentdrawing.cr, xc, yc, radius, angle1, angle2)
  do_action(action)
end

carc(centerpoint::Point, radius, angle1, angle2, action=:nothing) =
  carc(centerpoint.x, centerpoint.y, radius, angle1, angle2, action)

"""
      arc2r(c1, p2, p3, action=:nothing)

Make a circular arc centered at `c1` that starts at `p2` and ends at `p3`, going clockwise.

`c1`-`p2` really determines the radius. If `p3` doesn't lie on the circular path, it will be used only as an indication of the arc's length, rather than its position.
"""

function arc2r(c1, p2, p3, action=:nothing)
    r = norm(c1, p2)
    startangle = atan2(p2.y - c1.y, p2.x - c1.x)
    endangle   = atan2(p3.y - c1.y, p3.x - c1.x)
    if startangle < endangle
      startangle = mod2pi(startangle + 2pi)
    end
    arc(c1, r, startangle, endangle, action)
end

"""
Create a rectangle with one corner at (`xmin`/`ymin`) with width `w` and height `h` and do
an action.

    rect(xmin, ymin, w, h, action)

See `box()` for more ways to do similar things, such as supplying two opposite corners,
placing by centerpoint and dimensions.
"""

function rect(xmin, ymin, w, h, action=:nothing)
    if action != :path
        newpath()
    end
    Cairo.rectangle(currentdrawing.cr, xmin, ymin, w, h)
    do_action(action)
end

"""
Create a rectangle with one corner at `cornerpoint` with width `w` and height `h` and do an
action.

    rect(cornerpoint, w, h, action)
"""
function rect(cornerpoint::Point, w, h, action)
    rect(cornerpoint.x, cornerpoint.y, w, h, action)
end

"""
Create a rectangle between two points and do an action.

    box(cornerpoint1, cornerpoint2, action=:nothing)

"""
function box(corner1::Point, corner2::Point, action=:nothing)
    rect(corner1.x, corner1.y, corner2.x - corner1.x, corner2.y - corner1.y, action)
end

"""
Create a box/rectangle using the first two points of an array of Points to defined
opposite corners.

    box(points::Array, action=:nothing)
"""
function box(bbox::Array, action=:nothing)
    box(bbox[1], bbox[2], action)
end

"""
Create a box/rectangle centered at point `pt` with width and height.

    box(pt::Point, width, height, action=:nothing)
"""
function box(pt::Point, width, height, action=:nothing)
    rect(pt.x - width/2, pt.y - height/2, width, height, action)
end

"""
Create a box/rectangle centered at point `x/y` with width and height.

    box(x, y, width, height, action=:nothing)
"""
function box(x, y, width, height, action=:nothing)
    rect(x - width/2, y - height/2, width, height, action)
end

"""
    sector(innerradius, outerradius, startangle, endangle, action=:none)

Make an annular sector centered at the current `0/0` point.
"""
function sector(innerradius, outerradius, startangle, endangle, action=:none)
    newpath()
    move(innerradius * cos(startangle), innerradius * sin(startangle))
    line(outerradius * cos(startangle), outerradius * sin(startangle))
    arc(0, 0, outerradius, startangle, endangle,:none)
    line(innerradius * cos(endangle), innerradius * sin(endangle))
    carc(0, 0, innerradius, endangle, startangle, :none)
    closepath()
    do_action(action)
end

"""
    pie(x, y, radius, startangle, endangle, action=:none)

Make a pie shape centered at `x`/`y`. Angles start at the positive x-axis and are measured
clockwise.
"""

function pie(x, y, radius, startangle, endangle, action=:none)
    gsave()
    translate(x, y)
    newpath()
    move(0, 0)
    line(radius * cos(startangle), radius * sin(startangle))
    arc(0, 0, radius, startangle, endangle, :none)
    closepath()
    grestore()
    do_action(action)
end

"""
    pie(centerpoint, radius, startangle, endangle, action=:none)

Make a pie shape centered at `centerpoint`.

Angles start at the positive x-axis and are measured clockwise.
"""

pie(centerpoint::Point, radius, startangle, endangle, action) =
 pie(centerpoint.x, centerpoint.y, radius, startangle, endangle, action)

"""
Set the line width.

    setline(n)

"""
setline(n) = Cairo.set_line_width(currentdrawing.cr, n)

"""
Set the line ends. `s` can be "butt" (the default), "square", or "round".

    setlinecap(s)

    setlinecap("round")
"""
function setlinecap(str="butt")
    if str == "round"
        Cairo.set_line_cap(currentdrawing.cr, Cairo.CAIRO_LINE_CAP_ROUND)
    elseif str == "square"
        Cairo.set_line_cap(currentdrawing.cr, Cairo.CAIRO_LINE_CAP_SQUARE)
    else
        Cairo.set_line_cap(currentdrawing.cr, Cairo.CAIRO_LINE_CAP_BUTT)
    end
end

"""
Set the line join, or how to render the junction of two lines when stroking.

    setlinejoin("round")
    setlinejoin("miter")
    setlinejoin("bevel")
"""
function setlinejoin(str="miter")
    if str == "round"
        Cairo.set_line_join(currentdrawing.cr, Cairo.CAIRO_LINE_JOIN_ROUND)
    elseif str == "bevel"
        Cairo.set_line_join(currentdrawing.cr, Cairo.CAIRO_LINE_JOIN_BEVEL)
    else
        Cairo.set_line_join(currentdrawing.cr, Cairo.CAIRO_LINE_JOIN_MITER)
    end
end

"""
Set the dash pattern to one of: "solid", "dotted", "dot", "dotdashed", "longdashed",
"shortdashed", "dash", "dashed", "dotdotdashed", "dotdotdotdashed"

    setlinedash("dot")
"""
function setdash(dashing)
    Cairo.set_line_type(currentdrawing.cr, dashing)
end

"""
Move to a point.

    move(x, y)
    move(pt)
"""
move(x, y)      = Cairo.move_to(currentdrawing.cr,x, y)
move(pt)        = move(pt.x, pt.y)

"""
Move by an amount from the current point. Move relative to current position by `x` and `y`:

    rmove(x, y)

Move relative to current position by the `pt`'s x and y:

    rmove(pt)
"""
rmove(x, y)     = Cairo.rel_move(currentdrawing.cr,x, y)
rmove(pt)       = rmove(pt.x, pt.y)

"""
Create a line from the current position to the `x/y` position and optionally apply an action:

    line(x, y)
    line(x, y, :action)
    line(pt)

"""
line(x, y)      = Cairo.line_to(currentdrawing.cr,x, y)
line(pt)        = line(pt.x, pt.y)

"""
Make a line between two points, `pt1` and `pt2`.

    line(pt1::Point, pt2::Point, action=:nothing)

"""
function line(pt1::Point, pt2::Point, action=:nothing)
    move(pt1)
    line(pt2)
    do_action(action)
end

"""
Create a line relative to the current position to the `x/y` position and optionally apply an action:

    rline(x, y)
    rline(x, y, :action)
    rline(pt)
"""
rline(x, y)     = Cairo.rel_line_to(currentdrawing.cr,x, y)
rline(pt)       = rline(pt.x, pt.y)

"""
Create a cubic Bézier spline curve.

    curve(x1, y1, x2, y2, x3, y3)
    curve(p1, p2, p3)

The spline starts at the current position, finishing at `x3/y3` (`p3`), following two
control points `x1/y1` (`p1`) and `x2/y2` (`p2`)

"""

curve(x1, y1, x2, y2, x3, y3) = Cairo.curve_to(currentdrawing.cr, x1, y1, x2, y2, x3, y3)
curve(pt1, pt2, pt3)          = curve(pt1.x, pt1.y, pt2.x, pt2.y, pt3.x, pt3.y)

saved_colors = Tuple{Float64,Float64,Float64,Float64}[]

# I originally used simple Cairo save() but the colors/opacity
# thing I've got going didn't save/restore properly, hence the stack
"""
Save the current graphics state on the stack.
"""
function gsave()
    Cairo.save(currentdrawing.cr)
    push!(saved_colors,
        (currentdrawing.redvalue,
         currentdrawing.greenvalue,
         currentdrawing.bluevalue,
         currentdrawing.alpha))
end
"""
Replace the current graphics state with the one on top of the stack.
"""
function grestore()
    Cairo.restore(currentdrawing.cr)
    try
    (currentdrawing.redvalue,
     currentdrawing.greenvalue,
     currentdrawing.bluevalue,
     currentdrawing.alpha) = pop!(saved_colors)
     catch err
         println("$err Not enough colors on the stack to restore.")
    end
end

"""
Scale subsequent drawing in x and y.

Example:

    scale(0.2, 0.3)

"""

scale(sx, sy) = Cairo.scale(currentdrawing.cr, sx, sy)

"""
Rotate subsequent drawing by `a` radians clockwise.

    rotate(a)

"""

rotate(a) = Cairo.rotate(currentdrawing.cr, a)

"""
Translate to new location.

    translate(x, y)

or

    translate(point)
"""

translate(tx, ty)        = Cairo.translate(currentdrawing.cr, tx, ty)
translate(pt::Point)     = translate(pt.x, pt.y)

"""
Get the current path (thanks Andreas Lobinger!)

Returns a CairoPath which is an array of .element_type and .points. With the results you
could typically step through and examine each entry:

```
o = getpath()
for e in o
      if e.element_type == Cairo.CAIRO_PATH_MOVE_TO
          (x, y) = e.points
          move(x, y)
      elseif e.element_type == Cairo.CAIRO_PATH_LINE_TO
          (x, y) = e.points
          # straight lines
          line(x, y)
          stroke()
          circle(x, y, 1, :stroke)
      elseif e.element_type == Cairo.CAIRO_PATH_CURVE_TO
          (x1, y1, x2, y2, x3, y3) = e.points
          # Bezier control lines
          circle(x1, y1, 1, :stroke)
          circle(x2, y2, 1, :stroke)
          circle(x3, y3, 1, :stroke)
          move(x, y)
          curve(x1, y1, x2, y2, x3, y3)
          stroke()
          (x, y) = (x3, y3) # update current point
      elseif e.element_type == Cairo.CAIRO_PATH_CLOSE_PATH
          closepath()
      else
          error("unknown CairoPathEntry " * repr(e.element_type))
          error("unknown CairoPathEntry " * repr(e.points))
      end
  end
```
"""
getpath()      = Cairo.convert_cairo_path_data(Cairo.copy_path(currentdrawing.cr))

"""
Get the current path, like `getpath()` but flattened so that there are no Bezier curves.

Returns a CairoPath which is an array of .element_type and .points.
"""
getpathflat()  = Cairo.convert_cairo_path_data(Cairo.copy_path_flat(currentdrawing.cr))

# patterns not yet looked at these

#=
Cairo.pattern_create_radial(cx0::Real, cy0::Real, radius0::Real, cx1::Real, cy1::Real, radius1::Real)
Cairo.pattern_create_linear(x0::Real, y0::Real, x1::Real, y1::Real)
Cairo.pattern_add_color_stop_rgb(pat::CairoPattern, offset::Real, red::Real, green::Real, blue::Real)
Cairo.pattern_add_color_stop_rgba(pat::CairoPattern, offset::Real, red::Real, green::Real, blue::Real, alpha::Real)
Cairo.pattern_set_filter(p::CairoPattern, f)
Cairo.pattern_set_extend(p::CairoPattern, val)
Cairo.set_source(dest::CairoContext, src::CairoPattern)
=#

"""
    setcolor(col::String)

Set the current color to a named color. This relies on Colors.jl to convert a string to RGBA
eg setcolor("gold") # or "green", "darkturquoise", "lavender" or what have you. The list is at `Colors.color_names`.

    setcolor("gold")
    setcolor("darkturquoise")

Use `sethue()` for changing colors without changing current opacity level.
"""

function setcolor(col::String)
    temp = parse(RGBA, col)
    currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, currentdrawing.alpha = temp.r, temp.g, temp.b, temp.alpha
    Cairo.set_source_rgba(currentdrawing.cr, currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, temp.alpha)
    return (temp.r, temp.g, temp.b, temp.alpha)
end

"""
Set the current color.

    setcolor(r, g, b)
    setcolor(r, g, b, alpha)
    setcolor(color)
    setcolor(col::ColorTypes.Colorant)

Examples:

    setcolor(convert(Colors.HSV, Colors.RGB(0.5, 1, 1)))
    setcolor(.2, .3, .4, .5)
    setcolor(convert(Colors.HSV, Colors.RGB(0.5, 1, 1)))

    for i in 1:15:360
       setcolor(convert(Colors.RGB, Colors.HSV(i, 1, 1)))
       ...
    end
"""
function setcolor(col::ColorTypes.Colorant)
  temp = convert(RGBA, col)
  setcolor(temp.r, temp.g, temp.b)
  return (temp.r, temp.g, temp.b)
end

function setcolor(r, g, b, a=1)
    currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue,
      currentdrawing.alpha = r, g, b, a
    Cairo.set_source_rgba(currentdrawing.cr, r, g, b, a)
    return (r, g, b, a)
end

"""
Set the current color to a string.

For example:

    @setcolor"red"
"""
macro setcolor_str(ex)
    isa(ex, String) || error("colorant requires literal strings")
    col = parse(RGBA, ex)
    quote
    currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue,
      currentdrawing.alpha = $col.r, $col.g, $col.b, $col.alpha
    Cairo.set_source_rgba(currentdrawing.cr, currentdrawing.redvalue,
      currentdrawing.greenvalue, currentdrawing.bluevalue, $col.alpha)
    end
end

"""
Set the color. `sethue()` is like `setcolor()`, but (like Mathematica) we sometimes want to
change the current 'color' without changing alpha/opacity. Using `sethue()` rather than
`setcolor()` doesn't change the current alpha opacity.

    sethue("black")
    sethue(0.3,0.7,0.9)

"""
function sethue(col::String)
    temp = parse(RGBA,  col)
    currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue = temp.r, temp.g, temp.b
    Cairo.set_source_rgba(currentdrawing.cr, currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, currentdrawing.alpha) # use current alpha, not incoming one
    return (temp.r, temp.g, temp.b)
end

"""
Set the color to a named color:

    sethue("red")

"""
function sethue(col::ColorTypes.Colorant)
    temp = convert(RGBA,  col)
    currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue = temp.r, temp.g, temp.b
    # use current alpha
    Cairo.set_source_rgba(currentdrawing.cr, temp.r, temp.g, temp.b, currentdrawing.alpha)
    return (temp.r, temp.g, temp.b, currentdrawing.alpha)
end

"""
Set the color's `r`, `g`, `b` values:

    sethue(0.3, 0.7, 0.9)

Use `setcolor(r,g,b,a)` to set transparent colors.
"""

function sethue(r, g, b)
    currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue = r, g, b
    # use current alpha
    Cairo.set_source_rgba(currentdrawing.cr, r, g, b, currentdrawing.alpha)
    return (r, g, b, currentdrawing.alpha)
end

"""
Set the current opacity to a value between 0 and 1.

    setopacity(alpha)

This modifies the alpha value of the current color.
"""
function setopacity(a)
    # use current RGB values
    currentdrawing.alpha = a
    Cairo.set_source_rgba(currentdrawing.cr, currentdrawing.redvalue, currentdrawing.greenvalue, currentdrawing.bluevalue, currentdrawing.alpha)
    return a
end

"""
Set a random hue.

    randomhue()

Choose a random color without changing the current alpha opacity.
"""

function randomhue()
  rrand, grand, brand = rand(3)
  sethue(rrand, grand, brand)
  return (rrand, grand, brand)
end

"""
Set a random color.

    randomcolor()

This probably changes the current alpha opacity too.

"""
function randomcolor()
  rrand, grand, brand, arand = rand(4)
  setcolor(rrand, grand, brand, arand)
  return (rrand, grand, brand, arand)
end

"""
Get the current matrix.

    getmatrix()

Returns the current Cairo matrix as an array. In Cairo/Luxor, a matrix is an array of six float64 numbers:

- xx component of the affine transformation
- yx component of the affine transformation
- xy component of the affine transformation
- yy component of the affine transformation
- x0 translation component of the affine transformation
- y0 translation component of the affine transformation

Some basic matrix transforms:

- translate
  `transform([1,  0, 0,  1, dx, dy])`
  => shift by `dx`, `dy`

- scale
  `transform([fx, 0, 0, fy,  0, 0])`
  => scale by `fx`, `fy`

- rotate
  `transform([cos(a), sin(a), -cos(a), cos(a), 0, 0])`
  => rotate to `a` radians

- x-skew
  `transform([1,  0, tan(a), 1, 0, 0])`
  => xskew by `a`

- y-skew
  `transform([1, tan(a), 0, 1, 0, 0])`
  => yskew by `a`

- flip
  `transform([fx, 0, 0, fy, centerx * (1 - fx), centery * (fy-1)])`
  => flip with center at `centerx`/`centery`

When a drawing is first created, the matrix looks like this:

    getmatrix() = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0]

When the origin is moved to 400/400, it looks like this:

    getmatrix() = [1.0, 0.0, 0.0, 1.0, 400.0, 400.0]

To reset the matrix to the original:

    setmatrix([1.0, 0.0, 0.0, 1.0, 0.0, 0.0])

"""
function getmatrix()
    gm = Cairo.get_matrix(currentdrawing.cr)
    return([gm.xx, gm.yx, gm.xy, gm.yy, gm.x0, gm.y0])
end

"""
Change the current Cairo matrix to matrix `m`.

    setmatrix(m::Array)

Use `getmatrix()` to get the current matrix.
"""
function setmatrix(m::Array)
    if eltype(m) != Float64
        m = map(Float64,m)
    end
    # some matrices make Cairo freak out and need reset. Not sure what the rules are yet…
    if length(m) < 6
        throw("didn't like that matrix $m: not enough values")
    elseif countnz(m) == 0
        throw("didn't like that matrix $m: too many zeroes")
    else
        cm = Cairo.CairoMatrix(m[1], m[2], m[3], m[4], m[5], m[6])
        Cairo.set_matrix(currentdrawing.cr, cm)
    end
end

"""
Modify the current matrix by multiplying it by matrix `a`.

    transform(a::Array)

For example, to skew the current state by 45 degrees in x and move by 20 in y direction:

    transform([1, 0, tand(45), 1, 0, 20])

Use `getmatrix()` to get the current matrix.
"""
function transform(a::Array)
    b = Cairo.get_matrix(currentdrawing.cr)
    setmatrix([
        (a[1] * b.xx)  + a[2]  * b.xy,             # xx
        (a[1] * b.yx)  + a[2]  * b.yy,             # yx
        (a[3] * b.xx)  + a[4]  * b.xy,             # xy
        (a[3] * b.yx)  + a[4]  * b.yy,             # yy
        (a[5] * b.xx)  + (a[6] * b.xy) + b.x0,     # x0
        (a[5] * b.yx)  + (a[6] * b.yy) + b.y0      # y0
    ])
end

"""
    readpng(pathname)

Read a PNG file into Cairo.

This returns a image object suitable for placing on the current drawing with `placeimage()`.
You can access its width and height properties:

    image = readpng("/tmp/test-image.png")
    w = image.width
    h = image.height
"""
function readpng(pathname)
    return Cairo.read_from_png(pathname)
end

function paint_with_alpha(ctx::Cairo.CairoContext, a = 0.5)
    Cairo.paint_with_alpha(currentdrawing.cr, a)
end

"""
Place a PNG image on the drawing.

    placeimage(img, xpos, ypos)

The image `img` has been previously loaded using `readpng()`.
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos)
    Cairo.set_source_surface(currentdrawing.cr, img, xpos, ypos)
    # no alpha
    Cairo.paint(currentdrawing.cr)
end

"""
Place a PNG image on the drawing.

    placeimage(img, pos)
"""
placeimage(img::Cairo.CairoSurface, pt::Point) =
placeimage(img::Cairo.CairoSurface, pt.x, pt.y)

"""
Place a PNG image on the drawing using alpha transparency.

    placeimage(img, xpos, ypos, a)
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos, alpha)
    Cairo.set_source_surface(currentdrawing.cr, img, xpos, ypos)
    paint_with_alpha(currentdrawing.cr, alpha)
end

"""
Place a PNG image on the drawing using alpha transparency.

    placeimage(img, pos, a)
"""
placeimage(img::Cairo.CairoSurface, pt::Point, alpha) =
  placeimage(img::Cairo.CairoSurface, pt.x, pt.y, alpha)

end
# module


text:
import LinearAlgebra

# TODO: Orientation: should relatedness be calculated for row or column sets
# TODO: Test if dimensions match the input's: for colwise relatedness, output mat
# should be col x col size.
# TODO: Test result

function similarity(O::Array{Int64, 2}; method::String, orientation="colwise")

	# Prepare inputs --------------------------

	# If the sets that should be related are in rows, flip O
	@assert orientation ∈ ["colwise", "rowwise"] "`orientiation` needs to be either `colwise` or `rowwise`"

	O = orientation == "rowwise" ? O' : O

	# Get co-occurrence matrix + dims
	C = O'O
	n = size(O, 2)

	# No sets are allowed to be completely empty
	@assert minimum(LinearAlgebra.diag(C)) != 0 "Each relevant set must have at least one non-zero value"

	# Define possible similarity methods:
	method_dict = Dict("association" => (Cij, Si, Sj) -> Cij / (Si * Sj),
			   "cosine" => (Cij, Si, Sj) -> Cij / sqrt(Si * Sj), # Cosine index
			   "inclusion" => (Cij, Si, Sj) -> Cij / min(Si, Sj), # Inclusion index
			   "jaccard" => (Cij, Si, Sj) -> Cij / (Si + Sj - Cij), # Jaccard index
			   "proximity" => (Cij, Si, Sj) -> Cij / max(Si, Sj)) # Proximity

	# Make sure that the chosen method is available:
	@assert haskey(method_dict, method) "The chosen method, `$method`, is not specified correctly. Options are: `association` for association strength, `cosine` for cosine similarity, `inclusion` for inclusion index, `jaccard` for jaccard index, and `proximity` for proximity similarity."

	# Get similarity matrix: --------------------------

	similarity_fun = method_dict[method] # Select method
	SIM = zeros(Float64, n, n) # Initialize collection matrix

	# Calculate similarity:
	for i in 1:n
		for j in 1:n
			SIM[i, j] = similarity_fun(C[i, j], C[i, i], C[j, j])
		end
	end

	# Finish up: ---------------------------------------
	return SIM
end

function density(;M::Array{Int64, 2}, SIM::Array{Float64, 2})

	# Check inputs ------------------------------------
	# Make sure that SIM is symmetrical
	@assert size(SIM, 1) == size(SIM, 2) "`SIM` needs to be a square matrix"
	# Make sure that the number of products in M matches SIM
	@assert size(M, 2) == sizE(SIM, 2) "`M` and `SIM` needs to have the same number of columns"

	# Get local density -------------------------------
	# Set SIM diagonal to 0: a products similarity to itself is not important
	SIM[LinearAlgebra.diagind(SIM)] .= 0

	# Get country density around each product they export. D[c, p] is the
	# total "internal" similarity between product p to all the other products
	# exported by country c
	DENSITY = M * SIM

	# Get global density ------------------------------
	# Get vector of global similarities
	global_sim = SIM * ones(size(SIM, 2), 1)

	# Normalize each country-density value by global density
	for p in 1:size(global_sim, 1)
		DENSITY[:, p] = DENSITY[:, p] ./ global_sim[p]
	end

	# Finish up ---------------------------------------
	return DENSITY

end

ext: jl

max_stars_repo_name: sbpost/compost

max_issues_repo_name: sbpost/compost

max_forks_repo_name: sbpost/compost


text:
using WAVI, Plots
# Grid and boundary conditions
nx = 320
ny = 40
nσ = 4
x0 = 0.0
y0 = -40000.0
dx = 2000.0
dy = 2000.0
h_mask = trues(nx, ny)
u_iszero = falses(nx + 1, ny); u_iszero[1,:] .= true
v_iszero = falses(nx, ny + 1); v_iszero[:,1] .= true; v_iszero[:,end] .= true
grid = Grid(nx=nx,
            ny=ny,
            nσ=nσ,
            x0=x0,
            y0=y0,
            dx=dx,
            dy=dy,
            h_mask=h_mask,
            u_iszero=u_iszero,
            v_iszero=v_iszero)

# Bed
bed = WAVI.mismip_plus_bed # function definition

# Inputing thickness profile, reading from binary file
#fname = "examples\\Favier2019_melt_params\\data\\MISMIP_ice0_2km_SteadyThickness.bin";
fname = joinpath(dirname(@__FILE__), "data",  "MISMIP_ice0_2km_SteadyThickness.bin")
h = Array{Float64,2}(undef, nx, ny)
read!(fname, h)
h = ntoh.(h)

# make the model
if ~(@isdefined model)
initial_conditions = InitialConditions(initial_thickness=h) # set thickness
model = Model(grid=grid,
            bed_elevation=bed,
            initial_conditions=initial_conditions,
            solver_params=SolverParams(maxiter_picard=1) # we'll update state to get basal melt, but don't need other info accurate
            )

# update the configuration
end
update_state!(model)

ice_front_mask = zeros(nx,ny);
ice_front_mask[end,:] .= 1;
pico = PICO(ice_front_mask = ice_front_mask,
            T0 =1.2,
            S0 = 34.6,
            use_box_mean_depth = true,
            γT = 3.e-5);

#return the gl mask
ice_front_mask = zeros(nx,ny);
ice_front_mask[end,:] .= 1;
zb = model.fields.gh.b .* (model.fields.gh.grounded_fraction .== 1) + - pico.ρi / pico.ρw .* model.fields.gh.h .* (model.fields.gh.grounded_fraction .< 1)
gl_mask = WAVI.get_gl_mask(model.fields.gh.grounded_fraction)
boxes = WAVI.get_pico_boxes(gl_mask,ice_front_mask, model.fields.gh.grounded_fraction, pico.nbox )

WAVI.update_melt_rate!(pico, model.fields,model.grid)

#contour plot the boxes
x = model.grid.xxh[:,1]; y = model.grid.yyh[1,:];
plt1 = contour(x ./ 1e3,y / 1e3, boxes',
                fill=true,
                linewidth=0,
                colorbar=true,
                colorbar_title="pico boxes",
                framestyle=:box)
xlims!((420, 640))
xlabel!("x (km)")
ylabel!("y (km)")
plt1

#contour plot melt rate
m = deepcopy(model.fields.gh.basal_melt);
m[model.fields.gh.grounded_fraction .== 1.] .= NaN;
msat = deepcopy(m)
msat[msat .> 50] .= 50
#msat[zb .> -300] .= NaN;
plt2 = contour(x ./ 1e3,y / 1e3, msat',
                fill=true,
                linewidth=0,
                colorbar=true,
                colorbar_title="melt rate (m/yr)",
                framestyle=:box,
                colorscale = "Jet")
xlims!((420, 640))
xlabel!("x (km)")
ylabel!("y (km)")

m = deepcopy(model.fields.gh.basal_melt);
idx = (zb .< -300) .& (model.fields.gh.basal_melt .> 0)
mean_melt = sum(m[idx])./length(m[idx])
println("mean melt rate for shelf points below -300m depth is $mean_melt")


text:
using LinearAlgebra
using LinearAlgebra.BLAS
import LinearAlgebra.BLAS: gemm!, @blasfunc, libblas, BlasInt

export A_mul_B!, At_mul_B!, A_mul_Bt!, At_mul_B!

function gemm!(ta::Char, tb::Char, alpha::Float64, a::Union{Ref{Float64}, AbstractVecOrMat{Float64}},
               ma::Int64, na::Int64, b::Union{Ref{Float64}, AbstractVecOrMat{Float64}}, nb::Int64,
               beta::Float64, c::Union{Ref{Float64}, AbstractVecOrMat{Float64}})
    ccall((@blasfunc(dgemm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{BlasInt}, Ref{Float64}, Ptr{Float64}, Ref{BlasInt},
           Ptr{Float64}, Ref{BlasInt}, Ref{Float64}, Ptr{Float64},
           Ref{BlasInt}),
          ta, tb, ma, nb,
          na, alpha, a, max(ma, 1),
          b, max(na, 1), beta, c,
          max(ma, 1))
end

function A_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  offset_a::Int64, ma::Int64, na::Int64, b::AbstractVecOrMat{Float64},
                  offset_b::Int64, nb::Int64)

    gemm!('N', 'N', 1.0, Ref(a, offset_a), ma, na, Ref(b, offset_b),
          nb, 0.0, Ref(c, offset_c))
end

function A_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  offset_a::Int64, ma::Int64, na::Int64, b::AbstractVecOrMat{Float64})
    nb = (ndims(b) > 1) ? size(b, 2) : 1
    A_mul_B!(c, offset_c, a, offset_a, ma, na, b, 1, nb)
end

function A_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  b::AbstractVecOrMat{Float64}, offset_b::Int64, nb::Int64)
    ma = size(a, 1)
    na = (ndims(a) > 1) ? size(a, 2) : 1
    A_mul_B!(c, offset_c, a, 1, ma, na, b, offset_b, nb)
end

function At_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  offset_a::Int64, ma::Int64, na::Int64, b::AbstractVecOrMat{Float64},
                   offset_b::Int64, nb::Int64)
    ccall((@blasfunc(dgemm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{BlasInt}, Ref{Float64}, Ptr{Float64}, Ref{BlasInt},
           Ptr{Float64}, Ref{BlasInt}, Ref{Float64}, Ptr{Float64},
           Ref{BlasInt}),
          'T', 'N', ma, nb,
          na, 1.0, Ref(a, offset_a), max(na, 1),
          Ref(b, offset_b), max(na, 1), 0.0, Ref(c, offset_c),
          max(ma, 1))
end

function At_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  b::AbstractVecOrMat{Float64}, offset_b::Int64, nb::Int64)
    na = size(a, 1)
    ma = (ndims(a) > 1) ? size(a, 2) : 1
    At_mul_B!(c, offset_c, a, 1, ma, na, b, offset_b, nb)
end

function A_mul_Bt!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  offset_a::Int64, ma::Int64, na::Int64, b::AbstractVecOrMat{Float64},
                  offset_b::Int64, nb::Int64)
    ccall((@blasfunc(dgemm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{BlasInt}, Ref{Float64}, Ptr{Float64}, Ref{BlasInt},
           Ptr{Float64}, Ref{BlasInt}, Ref{Float64}, Ptr{Float64},
           Ref{BlasInt}),
          'N', 'T', ma, nb,
          na, 1.0, Ref(a, offset_a), max(ma, 1),
          Ref(b, offset_b), max(nb, 1), 0.0, Ref(c, offset_c),
          max(ma, 1))
end

function A_mul_Bt!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  offset_a::Int64, ma::Int64, na::Int64, b::AbstractVecOrMat{Float64})
    nb = size(b, 1)
    A_mul_Bt!(c, offset_c, a, offset_a, ma, na, b, 1, nb)
end

function At_mul_Bt!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64},
                  offset_a::Int64, ma::Int64, na::Int64, b::AbstractVecOrMat{Float64},
                  offset_b::Int64, nb::Int64)
    ccall((@blasfunc(dgemm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{BlasInt}, Ref{Float64}, Ptr{Float64}, Ref{BlasInt},
           Ptr{Float64}, Ref{BlasInt}, Ref{Float64}, Ptr{Float64},
           Ref{BlasInt}),
          'T', 'T', ma, nb,
          na, 1.0, Ref(a, offset_a), max(na, 1),
          Ref(b, offset_b), max(ma, 1), 0.0, Ref(c, offset_c),
          max(ma, 1))
end


text:
using Documenter, Alistair

makedocs()


text:
# Distance metrics between chromosomes
export genetic_distance,
    positional_distance,
    constant_functional_distance,
    random_functional_distance,
    active_distance,
    distance

function genetic_distance(c1::Chromosome, c2::Chromosome)
    distance = 0
    if length(c1.genes) == length(c2.genes)
        distance = sum(abs.(c1.genes .- c2.genes)) / length(c1.genes)
    elseif length(c1.genes) > length(c2.genes)
        distance = (sum(abs.(c1.genes .-
                             [c2.genes; zeros(length(c1.genes)-length(c2.genes))]))/
                    length(c1.genes))
    else
        distance = (sum(abs.([c1.genes; zeros(length(c2.genes)-length(c1.genes))]
                             .- c2.genes))/
                    length(c2.genes))
    end
    distance
end

function positional_distance(c1::Chromosome, c2::Chromosome)
    pc1 = get_positions(c1)
    pc2 = get_positions(c2)
    distance = 0
    if length(pc1) == length(pc2)
        distance = sum(abs.(pc1 .- pc2)) / length(pc1)
    elseif length(pc1) > length(pc2)
        distance = (sum(abs.(pc1 .- [pc2; ones(length(pc1)-length(pc2))]))/length(pc1))
    else
        distance = (sum(abs.([pc1; ones(length(pc2)-length(pc1))] .- pc2))/length(pc2))
    end
    distance
end

function constant_functional_distance(c1::Chromosome, c2::Chromosome)
    distance = 0
    reset!(c1)
    reset!(c2)
    range_inps = range(-1, stop=1, length=10)
    for i=1:10
        inps = repeat([range_inps[i]], c1.nin)
        distance += sum(abs.(process(c1, inps) .- process(c2, inps)))/c1.nout
    end
    distance/10
end

function random_functional_distance(c1::Chromosome, c2::Chromosome)
    distance = 0
    reset!(c1)
    reset!(c2)
    range_inps = range(-1, stop=1, length=10)
    for i=1:10
        inps = 2*rand(c1.nin) .- 1.0
        distance += sum(abs.(process(c1, inps) .- process(c2, inps)))/c1.nout
    end
    distance/10
end

function active_distance(c1::Chromosome, c2::Chromosome)
    diff = 0.0
    smaller = min(length(c1.nodes), length(c2.nodes))
    larger = max(length(c1.nodes), length(c2.nodes))
    for i in 1:smaller
        if c1.nodes[i].active && c2.nodes[i].active
            diff += sum(abs.(get_genes(c1, i) .- get_genes(c2, i)))
        end
    end
    diff /= smaller * node_genes(c1)
end

function distance(c1::Chromosome, c2::Chromosome)
    eval(Meta.parse(string(Config.distance_method)))(c1, c2)
end


text:
using ArbTools
using Nemo
using Test

include("enclosemaximum.jl")
include("isolateroots.jl")
include("integrate-test.jl")


text:
@generated function Base.copy!{T}(dest::AbstractCuArray{T,2}, src::AbstractCuArray{T,2})
    CT = cstring(T)
    f = compile("""
    copy(Array<$CT,2> dest, Array<$CT,2> src) {
        int idx0 = threadIdx.x + blockIdx.x * blockDim.x;
        int idx1 = threadIdx.y + blockIdx.y * blockDim.y;
        if (idx0 >= src.dims[0] || idx1 >= src.dims[1]) return;
        dest(idx0,idx1) = src(idx0,idx1);
    }""")
    quote
        @assert size(dest) == size(src)
        launch($f, size(src), (dest,src))
        dest
    end
end

@generated function Base.copy!{T}(dest::AbstractCuArray{T,3}, src::AbstractCuArray{T,3})
    CT = cstring(T)
    f = compile("""
    copy(Array<$CT,$N> dest, Array<$CT,$N> src) {
        int idx0 = threadIdx.x + blockIdx.x * blockDim.x;
        int idx1 = threadIdx.y + blockIdx.y * blockDim.y;
        int idx2 = threadIdx.z + blockIdx.z * blockDim.z;
        if (idx0 >= src.dims[0] || idx1 >= src.dims[1] || idx2 >= src.dims[2]) return;
        dest(idx0,idx1,idx2) = src(idx0,idx1,idx2);
    }""")
    quote
        @assert size(dest) == size(src)
        launch($f, size(src), (dest,src))
        dest
    end
end

#Base.copy{T}(src::CuSubArray{T}) = copy!(similar(src), src)


text:
######################################################################
# gif.jl
#
# makes gifs of ergodic trajectory generation
######################################################################

using Reel
export gif

"""
`gif(em::ErgodicManager, tm::TrajectoryManager, trajectory_file::String="temp.csv"; show_score=true, fps::Int=20)`

Below is an example of how you might use this. You `must` set `logging` to true during the trajectory generation.

`xd, ud = pto_trajectory(em, tm, logging=true)`

`gif(em, tm, fps=15)`

The gif will be saved at `temp.gif`.
"""
function gif(em::ErgodicManager, tm::TrajectoryManager, trajectory_file::String="temp.csv"; show_score=true, fps::Int=17)
	# First, let's read the trajectories in...
	trajectories = readcsv(trajectory_file)

	frames = Frames(MIME("image/png"), fps=fps)

	# ok, loop through all xd
	num_rows, n = size(trajectories)
	#N = round(Int, num_rows / num_trajectories)
	N = tm.N + 1
	num_trajectories = round(Int, num_rows / N)

	# start with the first trajectory...
	xd = mat2traj(trajectories[1:N, :])
	plot(em, xd, show_score=show_score)
	push!(frames, gcf())
	close()

	# then do the rest
	for traj_idx = 2:num_trajectories
		start_idx = (traj_idx-1)*N + 1
		end_idx = start_idx + N - 1
		xd = mat2traj(trajectories[start_idx:end_idx,:])
		plot(em, xd, show_score=true)
		push!(frames, gcf())
		close()
	end

	write("temp.gif", frames)
end

# watch vehicle progress
function gif(em::ErgodicManager, xd::VVF; show_score=true, fps::Int=5)

	frames = Frames(MIME("image/png"), fps=fps)

	N = length(xd) - 1

	# plot ever-increasing pieces of the trajectory
	for n = 1:(N+1)
		plot(em, xd, n, show_score=show_score)
		push!(frames, gcf())
		close()
	end

	write("temp.gif", frames)
end

# watch vehicle progress
function gif(em::ErgodicManager, xd::VVF, vtm::VTM; show_score=true, fps::Int=5)

	frames = Frames(MIME("image/png"), fps=fps)

	N = length(xd) - 1

	# plot ever-increasing pieces of the trajectory
	for n = 1:(N+1)
		plot(em, xd, n, vtm, show_score=show_score)
		push!(frames, gcf())
		close()
	end

	write("temp.gif", frames)
end


text:
using Chess
function testMove(b::Board, m::Move)
        negamax(b, 3, true; evalFunc = LichessBot.reallySimpleEval) == m
end
function testMove(fen::AbstractString, m::Move)
        testMove(fromfen(fen), m)
end
function testMove(b, m::AbstractString)
        testMove(b, movefromstring(m))
end
boards = [
("5k2/pR4R1/3p3p/6p1/4N3/8/PPr2PPP/7K b - - 2 26", "c2c1"),
("rnbk4/pppp2b1/4pr2/P4p1Q/1q3P2/2N5/3R2P1/2B1KB1R b - - 2 28", "b4c3"),
("rnbk1b1Q/pppp4/4pr2/P4p2/5P2/2q5/3R2P1/2B1KB1R w - - 2 30", "c1b2"),
("rnb1kbr1/pppp1ppp/4p3/8/1qP4Q/1P6/P2PNPPP/RNB1KB1R w KQq - 3 8", "h4h7"),
("1r1kr3/1pp4p/7b/p2b1pp1/3B1PP1/PPnP2Q1/R2N1K1P/5BNR b - - 7 26", "c3d1")
]
@testset "Evaluation" begin
        for (bs, ms) in boards
                b = fromfen(bs)
                @test LichessBot.negamax(b, 3, true) == movefromstring(ms)
        end
        #=begin
                local b = startboard()
                @test !isnothing(LichessBot.oldChooseMove(b))
        end
        begin
                local b = fromfen(
                        "5k2/pR4R1/3p3p/6p1/4N3/8/PPr2PPP/7K b - - 2 26",
                )
                @test LichessBot.negamax(b, 3, true) == Move(SQ_C2, SQ_C1)
        end
        begin
                local b = fromfen(
                        "rnbk4/pppp2b1/4pr2/P4p1Q/1q3P2/2N5/3R2P1/2B1KB1R b - - 2 28",
                )
                @test negamax(b, 3, true) == Move(SQ_B4, SQ_C3)
        end
        begin
                local b = fromfen(
                        "rnbk1b1Q/pppp4/4pr2/P4p2/5P2/2q5/3R2P1/2B1KB1R w - - 2 30",
                )
                @test negamax(
                        b,
                        3,
                        true;
                        evalFunc = LichessBot.reallySimpleEval,
                ) == movefromstring("c1b2")
        end
        begin
                local b = fromfen(
                        "rnb1kbr1/pppp1ppp/4p3/8/1qP4Q/1P6/P2PNPPP/RNB1KB1R w KQq - 3 8",
                )
                @test negamax(
                        b,
                        3,
                        true;
                        evalFunc = LichessBot.reallySimpleEval,
                ) == movefromstring("h4h7")
        end
        begin
                local b = fromfen("1r1kr3/1pp4p/7b/p2b1pp1/3B1PP1/PPnP2Q1/R2N1K1P/5BNR b - - 7 26")
                @test negamax(
                b, 1, true
                ) == movefromstring("c3d1")
        end=#
        for (bs, ms) in boards
                b = fromfen(bs)
                @test negamax(b, 3, true) == negamax(b, 3, true; AB=false)
                @test negamax(b, 3) == negamax(b, 3; AB=false)
        end
end


text:
@testset "Construction" begin

    # basic constructor
    xs = [1, 2, 3]
    f(x) = x^2
    q = fit(Newton, xs, f)
    @test all(f.(xs) .≈ q.(xs))

    # test tableau
    @test all(q.tableau .≈ [1 3 1; 0 4 5; 0 0 9])

    # test one, zero, variable
    x = variable()
    @test one(q)(x) == one(x)
    @test iszero(zero(q))
    @test variable(q)(x) == x
end

@testset "Conversion" begin
    xs = [1, 2, 3]
    f(x) = x^2
    p = fit(Newton, xs, f)
    q = convert(Polynomial, p)
    @test coeffs(q) == [0, 0, 1]
end

@testset "Evaluation" begin
    xs = sort(rand(5))
    ys = sin.(xs)
    p = fit(Newton, xs, sin)
    @test p.(xs) ≈ sin.(xs)
    us = rand(5)
    @test maximum(abs.([p(x) - sin(x) for x in us])) <= 1e-1
end

@testset "Arithmetic" begin
    p = fit(Newton, [1, 2, 3], x -> x^2)
    q = fit(Newton, [2, 3, 4, 5], x -> x^3)

    us = rand(5)

    @test all(p.(us) + q.(us) .≈ (p + q).(us))
    @test p.(us) .* q.(us) ≈ (p * q).(us)
end

@testset "Elementwise operations" begin
    p = fit(Newton, [1, 2, 3, 4], x -> x^2)
    p + 1
    -p
    2p
end

@testset "roots" begin
    q = fit(Newton, [2, 3, 4, 5], x -> (x - 1) * (x - 1 / 2) * (x + 1 / 2))
    @test all(roots(q) .≈ [-1 / 2, 1 / 2, 1])
end


text:
include("utils.jl")

filename = []
title_names = []
# 1
push!(filename, "data/avg_earth_hs_he_12_hp_4_ve_12_vp_4_roefanov.jld2")
push!(title_names, "Earth")

# 2
push!(filename, "data/averages_small_earth_long_hs_he_12_hp_4_ve_12_vp_4_roefanov.jld2")
push!(title_names, "Small Earth")

s_strings = ["wT", "ww", "TT"]
sL_strings = [L"\langle w'T' \rangle", L"\langle w'w' \rangle", L"\langle T'T' \rangle"]

colorrange = (0, 360) # modify to extrema of "truth"
contour_levels = [0, 40, 80, 160, 240, 320]
fig = Figure(resolution = (1700+600, 1000+400))
add_label = false

for i in 1:6
    ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
    jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong

    jl_file = jldopen(filename[jj], "r+")

    s_string = s_strings[ii]
    if length(s_string) == 2
        slice_zonal, s_string = eddy_variance(s_string, jl_file)
        s_string = sL_strings[ii]
    else
        slice_zonal = grab_state(s_string, jl_file)
        s_string = sL_strings[ii]
    end
    colorrange, contour_levels, s_string = plot_helper(s_string, slice_zonal)
    println(colorrange)

    s_string = title_names[jj]
    λ, ϕ, r, p_coord = grab_grid(jl_file)
    ax = fig[jj,ii] = Axis(fig, title = s_string, titlesize = 40)
    if ii !==1
    contour_heatmap!(ax, ϕ, p_coord, slice_zonal, contour_levels, colorrange,
                     add_labels = add_label, colormap = :thermometer)
    else
        contour_heatmap!(ax, ϕ, p_coord, slice_zonal, contour_levels, colorrange,
                     add_labels = add_label, colormap = :thermometer,
                     random_seed = 10)
    end
    if ii !== 1
        hideydecorations!(ax, grid = false)
    end
end

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
# HttpServer module
#
# Serve HTTP requests in Julia.
#
isdefined(Main, :__precompile__) && __precompile__()

module SimpleHttpCommon

const VERSION=v"0.1.0-dev"
const MODULE = current_module()
const NAME = module_name(MODULE)

using SimpleHttpIO

const N = SimpleHttpIO.N
const N_TYPE = SimpleHttpIO.N_TYPE
const STR_TYPE = SimpleHttpIO.STR_TYPE
import URIParser: URI, unescape

include("Chunk.jl")

export readlinebare,
       readlinestr,
       readlinebarestr,
       AbstractIOSocket,
       writeln,
       CRLF,
       AbstractIOSocket,
       IOSocket,
       readbyte,
       BYTE_TYPE

import DataStructures: OrderedDict
import Base: readbytes,
             readline,
             readuntil
using HttpCommon: STATUS_CODES,
       escapeHTML,
       mimetypes

import .Chunk

export IOChunked,
    AbstractIOSocket

include("common.jl")
include("protocol_api.jl")
include("request.jl")
include("response.jl")
include("protocols.jl")

end


text:
module LossFunction

using LinearAlgebra
using Statistics

LOSSES = [
    :mse, :cee, :mae, :huber, :logcosh_loss, :poisson, :hinge, :smooth_hinge
]

for f in LOSSES
    @eval export $(f)
end

safe_log(x, miniv = 1e-8) = map(log, clamp.(x, miniv, Inf))

@doc raw"""
    mse(y, t; reduction="mean")
Mean Square Error. This is the expression:
```math
MSE(y, t) = \frac{\sum_{i=1}^{n} (t_{i}-y_{i})^{2}}{n}
```
"""
function mse(y::AbstractVector, t::AbstractVector; reduction::String="mean")
    loss = (y.-t).^2
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end
mse(y::Number, t::Number) = (y-t)^2

@doc raw"""
    cee(y, t; reduction="mean")
Cross Entropy Error. This is the expression:
```math
CEE(y, t) = \frac{\sum_{i=1}^{n} t\ln y}{n}
```
"""
function cee(y::AbstractVector, t::AbstractVector; reduction::String="mean")
    loss = -t.*safe_log(y)
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end

cee(y::Number, t::Number) = -t*safe_log(y)

@doc raw"""
    mae(y, t)
Mean Absolute Error. This is the expression:
```math
MAE(y, t) = \frac{\sum_{i=1}^{n} |t_{i}-y_{i}|}{n}
```
"""
function mae(y::AbstractVector, t::AbstractVector; reduction::String="mean")
    loss = abs.(t-y)
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end

mae(y::Number, t::Number) = abs(t-y)

@doc raw"""
    huber(y, t; δ=1, reduction="mean")
Huber-Loss. If `δ` is large, it will be a function like [`mse`](@ref), and if it is small, it will be a function like [`mae`](@ref). This is the expression:
```math
a = |t_{i}-y_{i}| \\
Huber(y, t) = \frac{1}{n} \sum_{i=1}^{n} \left\{
\begin{array}{ll}
\frac{1}{2}a^{2} & (a \leq \delta) \\
\delta(a-\frac{1}{2}\delta) & (a \gt \delta)
\end{array}
\right.
```
"""
function huber(y::AbstractVector, t::AbstractVector; reduction::String="mean", δ=1)
    a = abs.(t-y)
    loss = @. ifelse(a<=δ, a^2/2, (a-δ/2)δ)
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end

huber(y::Number, t::Number; δ=1) = ifelse(abs(t-y)<=δ, abs(t-y)^2/2, (abs(t-y)-δ/2)δ)

@doc raw"""
    logcosh_loss(y, t; reduction="mean")
Log Cosh. Basically, it's [`mae`](@ref), but if the loss is small, it will be close to [`mse`](@ref). This is the expression:
```math
Logcosh(y, t) = \frac{\sum_{i=1}^{n} \log(\cosh(t_{i}-y_{i}))}{n}
```
"""
function logcosh_loss(y::AbstractVector, t::AbstractVector; reduction::String="mean")
    loss = @. log(cosh(t-y))
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end

logcosh_loss(y::Number, t::Number) = log(cosh(t-y))

@doc raw"""
    Poisson(y, t; reduction="mean")
Poisson Loss, Distribution of predicted value and loss of Poisson distribution. This is the expression:
```math
Poisson(y, t) = \frac{\sum_{i=1}^{n} y_{i}-t_{i} \ln y}{n}
```
"""
function poisson(y::AbstractVector, t::AbstractVector; reduction::String="mean")
    loss = @. y - t*log(y)
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end

poisson(y::Number, t::Number) = y - t*log(y)

@doc raw"""
    hinge(y, t; reduction="mean")
Hinge Loss, for SVM. This is the expression:
```math
Hinge(y, t) = \frac{\sum_{i=1}^{n} \max(1-y_{i}t_{i}, 0)}{n}
```
"""
function hinge(y::AbstractVector, t::AbstractVector; reduction::String="mean")
    loss = @. max(1-y*t, 0)
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end

hinge(y::Number, t::Number) = max(1-y*t, 0)

@doc raw"""
    smooth_hinge(y, t; reduction="mean")
Smoothing Hinge Loss. This is the expression:
```math
smoothHinge(y, t) = \frac{1}{n} \sum_{i=1}^{n} \left\{
\begin{array}{ll}
0 & (t_{i}y_{i} \geq 1) \\
\frac{1}{2}(1-t_{i}y_{i})^{2} & (0 \lt t_{i}y_{i} \lt 1) \\
\frac{1}{2} - t_{i}y_{i} & (t_{i}y_{i} \leq 0)
\end{array}
\right.
```
"""
function smooth_hinge(y::AbstractVector, t::AbstractVector; reduction::String="mean")
    z = t.*y
    loss = similar(z)
    for i in 1 : length(y)
        if z[i] >= 0
            loss[i] = 0
        elseif 0 < z[i] < 1
            loss[i] = (1-z[i])^2/2
        else
            loss[i] = 1/2-z[i]
        end
    end
    if reduction=="none"
        return loss
    elseif reduction=="sum"
        return sum(loss)
    elseif reduction=="mean"
        return mean(loss)
    else
        throw(ArgumentError("`reduction` must be either `none`, `sum` or `mean`"))
    end
end

function smooth_hinge(y::Number, t::Number)
    if t*y >= 0
        return 0
    elseif 0 < t*y < 1
        return (1-t*y)^2/2
    else
        return 1/2-t*y
    end
end

#TODO:Consider whether to implement Calback liverer information volume

for lossfunc in LOSSES
    #For Vector and Number
    @eval begin
        function $(lossfunc)(y::AbstractVector{T}, t::Number; reduction::String="mean") where {T}
            $(lossfunc)(y, fill(t, length(y)), reduction=reduction)
        end
    end
    #For Matrix and Number
    @eval begin
        function $(lossfunc)(y::AbstractMatrix{T}, t::Number; reduction::String="mean") where {T}
            if length(y) == 1
                return $(lossfunc)(y..., t)
            elseif size(y, 1)==1 || size(y, 2)==1
                return $(lossfunc)(vec(y), t, reduction=reduction)
            else
                throw(DimensionMismatch("LossFunctions don't support for Matrix!"))
            end
        end
    end
end

end


text:
"""
`Dict` containing pK values from different sources. Values are ordered [NTerm, CTerm, C, D, E, H, K, R, Y].
Available options are "EMBOSS", "Solomon", "Rodwell", "IPC protein", and "IPC peptide".
Values were collected from http://isoelectric.org/theory.html
"""
const pKtable = Dict{String, Vector{Float64}}(
    "EMBOSS" => [8.6, 3.6, 8.5, 3.9, 4.1, 6.5, 10.8, 12.5, 10.1],
    "Solomon" => [9.6, 2.4, 8.3, 3.9, 4.3, 6.0, 10.5, 12.5, 10.1],
    "Rodwell" => [8.0, 3.1, 8.33, 3.68, 4.25, 6.0, 11.5, 11.5, 10.07],
    "IPC protein" => [9.094, 2.869, 7.555, 3.872, 4.412, 5.637, 9.052, 11.84, 10.85],
    "IPC peptide" => [9.564, 2.383, 8.297, 3.887, 4.317, 6.018, 10.517, 12.503, 10.071]
)

"""
Amino acid weights (source: http://education.expasy.org/student_projects/isotopident/htdocs/aa-list.html)
"""
const aminoacidweights = Dict(
    AA_A => 71.0788,
    AA_R => 156.1875,
    AA_N => 114.1038,
    AA_D => 115.0886,
    AA_C => 103.1388,
    AA_Q => 128.1307,
    AA_E => 129.1155,
    AA_G => 57.0519,
    AA_H => 137.1411,
    AA_I => 113.1594,
    AA_L => 113.1594,
    AA_K => 128.1741,
    AA_M => 131.1926,
    AA_F => 147.1766,
    AA_P => 97.1167,
    AA_S => 87.0782,
    AA_T => 101.1051,
    AA_W => 186.2132,
    AA_Y => 163.1760,
    AA_V => 99.1326,
    AA_Z => 146.6385,
    AA_J => 131.175,
    AA_B => 132.6115,
    AA_U => 168.064,
    AA_O => 255.313,
    AA_X => 118.88603)

"""
Hydropathy indices (source: https://www.ncbi.nlm.nih.gov/pubmed/7108955)
"""
const hydropathy = Dict(
    AA_A => 1.8,
    AA_R => -4.5,
    AA_N => -3.5,
    AA_D => -3.5,
    AA_C => 2.5,
    AA_Q => -3.5,
    AA_E => -3.5,
    AA_G => -0.4,
    AA_H => -3.2,
    AA_I => 4.5,
    AA_L => 3.8,
    AA_K => -3.9,
    AA_M => 1.9,
    AA_F => 2.8,
    AA_P => -1.6,
    AA_S => -0.8,
    AA_T => -0.7,
    AA_W => -0.9,
    AA_Y => -1.3,
    AA_V => 4.2)


text:
using Xfoil, Test

# read airfoil data
f = open("naca2412.dat", "r")

x = Float64[]
y = Float64[]
for line in eachline(f)
    push!(x, Meta.parse(split(chomp(line))[1]))
    push!(y, Meta.parse(split(chomp(line))[2]))
end

close(f)

# begin tests

set_coordinates(x, y)

@testset "set_coordinates" begin
    for i = 1:length(x)
        @test Xfoil.xfoilglobals.xb[i] == x[i]
        @test Xfoil.xfoilglobals.yb[i] == y[i]
    end
end

pane(npan=140)

@testset "pane" begin
    @test Xfoil.xfoilglobals.n[1]==140
end

cl, cdaf, cdp, cm, conv = solve_alpha(5.0, 1e5, iter=100)

@testset "solve_alpha" begin
    @test isapprox(cl, 0.8029209409078045)
    @test isapprox(cdaf, 0.017368353583814947)
    @test isapprox(cdp, 0.0060056723724223405)
    @test isapprox(cm, -0.053541244481716696)
    @test conv == true
end

sbl, xbl, ybl, uebl, dstarbl, thetabl, cfbl = bldump()

@testset "bldump" begin
    @test isapprox(sum(sbl), 186.17483029646593)
    @test isapprox(sum(xbl), 81.11167420760269)
    @test isapprox(sum(ybl), 1.5695474130166154)
    @test isapprox(sum(uebl), 60.87432109391924)
    @test isapprox(sum(dstarbl), 0.9183533699318299)
    @test isapprox(sum(thetabl), 0.44113905597135494)
    @test isapprox(sum(cfbl), 1.3955896003688344)
end

aoas = -15.0:1.0:15.0

cls, cds, cdps, cms, convs = alpha_sweep(x, y, aoas, 1e5; iter=100,
    npan=140, percussive_maintenance=true, printdata=false, zeroinit=true,
    clmaxstop=true, clminstop=true)

@testset "alpha_sweep" begin
    # aoa=-5
    @test isapprox(cls[11], -0.48488210527938785)
    @test isapprox(cds[11], 0.023137034265852948)
    @test isapprox(cdps[11], 0.009548882396644153)
    @test isapprox(cms[11], -0.026256734143567897)

    # aoa=5
    @test isapprox(cls[21], 0.802920940907813)
    @test isapprox(cds[21], 0.01736835358381533)
    @test isapprox(cdps[21], 0.006005672372422114)
    @test isapprox(cms[21], -0.053541244481718306)
end

Xfoil.set_coordinates_cs(x, y)

@testset "set_coordinates_cs" begin
    for i = 1:length(x)
        @test isapprox(Xfoil.xfoilglobals.xb[i], real(Xfoil.xfoilglobals_cs.xb[i]))
        @test isapprox(Xfoil.xfoilglobals.yb[i], real(Xfoil.xfoilglobals_cs.yb[i]))
    end
end

pane_cs(npan=140)

@testset "pane_cs" begin
    for i = 1:140
        @test isapprox(Xfoil.xfoilglobals.x[i], real(Xfoil.xfoilglobals_cs.x[i]))
        @test isapprox(Xfoil.xfoilglobals.y[i], real(Xfoil.xfoilglobals_cs.y[i]))
    end
end

cl_cs, cd_cs, cdp_cs, cm_cs, conv_cs = solve_alpha_cs(5.0, 1e5, iter=100)

@testset "solve_alpha_cs" begin
    @test isapprox(cl, real(cl_cs))
    @test isapprox(cdaf, real(cd_cs))
    @test isapprox(cdp, real(cdp_cs))
    @test isapprox(cm, real(cm_cs))
    @test conv == conv_cs
end

sbl_cs, xbl_cs, ybl_cs, uebl_cs, dstarbl_cs, thetabl_cs, cfbl_cs = bldump_cs()

@testset "bldump_cs" begin
    @test isapprox(sum(sbl), sum(real(sbl_cs)))
    @test isapprox(sum(xbl), sum(real(xbl_cs)))
    @test isapprox(sum(ybl), sum(real(ybl_cs)))
    @test isapprox(sum(uebl), sum(real(uebl_cs)))
    @test isapprox(sum(dstarbl), sum(real(dstarbl_cs)))
    @test isapprox(sum(thetabl), sum(real(thetabl_cs)))
    @test isapprox(sum(cfbl), sum(real(cfbl_cs)))
end

cls_cs, cds_cs, cdps_cs, cms_cs, convs_cs = alpha_sweep_cs(x, y, aoas, 1e5;
    iter=100, npan=140, percussive_maintenance=false, printdata=false,
    zeroinit=true, clmaxstop=true, clminstop=true)

@testset "alpha_sweep_cs" begin
    # aoa=-5
    @test isapprox(cls[11], real(cls_cs[11]))
    @test isapprox(cds[11], real(cds_cs[11]))
    @test isapprox(cdps[11], real(cdps_cs[11]))
    @test isapprox(cms[11], real(cms_cs[11]))

    # aoa=5
    @test isapprox(cls[21], real(cls_cs[21]))
    @test isapprox(cds[21], real(cds_cs[21]))
    @test isapprox(cdps[21], real(cdps_cs[21]))
    @test isapprox(cms[21], real(cms_cs[21]))
end


text:
# ---
# title: 69. Sqrt(x)
# id: problem69
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/sqrtx/description/>
# hidden: true
# ---
#
# Given a non-negative integer `x`, compute and return _the square root of_ `x`.
#
# Since the return type is an integer, the decimal digits are **truncated** ,
# and only **the integer part** of the result is returned.
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 4
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: x = 8
#     Output: 2
#     Explanation: The square root of 8 is 2.82842..., and since the decimal part is truncated, 2 is returned.
#
#
#
# **Constraints:**
#
#   * `0 <= x <= 231 - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
fs()

@testset "ignore/fs2" begin
    gotd()
    s = """
        @def ignore = ["foo.md", "path/foo.md", "dir/", "path/dir/", r"index2.*"]
        """
    write(joinpath(td, "config.md"), s);
    F.process_config()
    @test F.globvar("ignore") == ["foo.md", "path/foo.md", "dir/", "path/dir/", r"index2.*"]

    write(joinpath(td, "foo.md"), "anything")
    mkpath(joinpath(td, "path"))
    write(joinpath(td, "path", "foo.md"), "anything")
    mkpath(joinpath(td, "dir"))
    write(joinpath(td, "dir", "foo1.md"), "anything")
    mkpath(joinpath(td, "path", "dir"))
    write(joinpath(td, "path", "dir", "foo2.md"), "anything")
    write(joinpath(td, "index.md"), "standard things")
    write(joinpath(td, "index2.md"), "anything")
    watched = F.fd_setup()
    @test length(watched.md) == 1
    @test first(watched.md).first.second == "index.md"
end

@testset "coverage" begin
    # form custom output path
    F.def_LOCAL_VARS!()
    path = F.form_custom_output_path("aa/bb")
    @test endswith(F.locvar(:fd_url), "aa/bb/index.html")
    @test isdir(dirname(path))

    # helper functions around regexes
    r1 = r"foo/bar"
    r2 = "foo/bar"
    @test F._access(r1) == r1.pattern
    @test !F._isempty(r1)
    @test !F._isempty(r2)
    @test F._isempty(r"")
    @test F._isempty("")
    @test !F._endswith(r1)
    @test F._endswith("foo/bar/")
    @test F._endswith(r"foo/bar/")
end


text:
struct Iteration
    i::Int
end

MatrixOf(::Type{Vector{T}}) where T = Matrix{T}
MatrixOf(::Type{CuVector{T}}) where T = CuMatrix{T}
zerosof(::Type{TM}, n...) where TM = (TM(undef, n...) .= 0)
onesof(::Type{TM}, n...) where TM = (TM(undef, n...) .= 1)
infsof(::Type{TM}, n...) where TM = (TM(undef, n...) .= Inf)
ninfsof(::Type{TM}, n...) where TM = (TM(undef, n...) .= -Inf)
nansof(::Type{TM}, n...) where TM = (TM(undef, n...) .= NaN)

@inline minus_plus(a, b) = a - b, a + b

@inline or(a,b) = a || b

macro matdot(v, A, j)
    r = gensym()
    T = gensym()
    esc(quote
        $T = promote_type(eltype($v), eltype($A))
        $r = zero($T)
        for i in 1:length($v)
            $r += $v[i] * $A[$j, i]
        end
        $r
    end)
end

function check_error(m, x0)
    if length(x0) != dim(m)
        throw(ArgumentError("initial variable must have same length as number of design variables"))
    end

    Threads.@threads for j in 1:length(x0)
        # x is not in box
        if !(min(m, j) <= x0[j] <= max(m,j))
            throw(ArgumentError("initial variable at index $j outside box constraint"))
        end
    end
end