text:
export notebook

####################################################
abstract type PlotKey end
const PlotId = Any #PlotlyJS.GenericTrace
const PlotRef = GenericRef{PlotKey, PlotId}
const PlotNativeRef = NativeRef{PlotKey, PlotId}
const PlotUnionRef = UnionRef{PlotKey, PlotId}
const PlotSubtractionRef = SubtractionRef{PlotKey, PlotId}

const PlotMaterial = String

struct PlotBackend{K,T} <: Backend{K,T}
  connection::LazyParameter{PlotlyJS.SyncPlot}
end

connection(b::PlotBackend{K,T}) where {K,T} = b.connection()

const PLOT = PlotBackend{PlotKey, PlotId}

has_boolean_ops(::Type{PLOT}) = HasBooleanOps{false}()
void_ref(b::PLOT) = PlotNativeRef(-1)

create_plot_connection() =
  let (width, height) = render_size(),
      layout =
        Layout(autosize=true, width=width, height=height,
               margin=attr(l=0, r=0, b=0, t=0),
               showlegend=false,
               scene=attr(
                 #camera=attr(
                #   center=attr(x=target.x,
                #               y=target.y,
                #               z=target.z),
                #   eye=attr(x=camera.x-target.x,
                #            y=camera.y-target.y,
                #            z=camera.z-target.z)),
                 aspectmode="data",
                 xaxis=attr(showticklabels=false),
                 yaxis=attr(showticklabels=false),
                 zaxis=attr(showticklabels=false))),
      s = PlotlyJS.SyncPlot(Plot(AbstractTrace[], layout))
    display(s)
    s
  end

notebook = PLOT(LazyParameter(PlotlyJS.SyncPlot, create_plot_connection))

reset_backend(b::PLOT) =
  let param = b.connection
    if ! isnothing(param.value)
      param.value = (copy(param.value))
      display(param.value)
    end
  end

new_backend(b::PLOT) =
  reset(b.connection)

#

set_view(camera::Loc, target::Loc, lens::Real, aperture::Real, b::Plot) =
  begin
    b.camera = camera
    b.target = target
    b.lens = lens
  end

get_view(b::Plot) =
  b.camera, b.target, b.lens

###################################
#
delete_all_shapes(b::PLOT) =
  let c = connection(b),
      idxs = collect(1:length(c.plot.data))
    if length(idxs) > 0
      deletetraces!(c, idxs...)
    end
  end

backend_delete_shapes(b::PLOT, shapes::Shapes) =
  let c = connection(b),
      idxs = [findfirst(==(ref(s).value), c.plot.data) for s in shapes]
    if length(idxs) > 0
      deletetraces!(c, idxs...)
    end
  end

#=
backend_stroke(b::ACAD, path::CircularPath) =
    @remote(b, Circle(path.center, vz(1, path.center.cs), path.radius))
backend_stroke(b::ACAD, path::RectangularPath) =
    let c = path.corner,
        dx = path.dx,
        dy = path.dy
        @remote(b, ClosedPolyLine([c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)]))
    end
backend_stroke(b::ACAD, path::ArcPath) =
    backend_stroke_arc(b, path.center, path.radius, path.start_angle, path.amplitude)

backend_stroke(b::ACAD, path::OpenPolygonalPath) =
  	@remote(b, PolyLine(path.vertices))
backend_stroke(b::ACAD, path::ClosedPolygonalPath) =
    @remote(b, ClosedPolyLine(path.vertices))
backend_fill(b::ACAD, path::ClosedPolygonalPath) =
    @remote(b, SurfaceClosedPolyLine(path.vertices))
backend_fill(b::ACAD, path::RectangularPath) =
    let c = path.corner,
        dx = path.dx,
        dy = path.dy
        @remote(b, SurfaceClosedPolyLine([c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)]))
    end
backend_stroke(b::ACAD, path::OpenSplinePath) =
  if (path.v0 == false) && (path.v1 == false)
    #@remote(b, Spline(path.vertices))
    @remote(b, InterpSpline(
                     path.vertices,
                     path.vertices[2]-path.vertices[1],
                     path.vertices[end]-path.vertices[end-1]))
  elseif (path.v0 != false) && (path.v1 != false)
    @remote(b, InterpSpline(path.vertices, path.v0, path.v1))
  else
    @remote(b, InterpSpline(
                     path.vertices,
                     path.v0 == false ? path.vertices[2]-path.vertices[1] : path.v0,
                     path.v1 == false ? path.vertices[end-1]-path.vertices[end] : path.v1))
  end
backend_stroke(b::ACAD, path::ClosedSplinePath) =
    @remote(b, InterpClosedSpline(path.vertices))
backend_fill(b::ACAD, path::ClosedSplinePath) =
    backend_fill_curves(b, @remote(b, InterpClosedSpline(path.vertices)))

backend_fill_curves(b::ACAD, refs::ACADIds) = @remote(b, SurfaceFromCurves(refs))
backend_fill_curves(b::ACAD, ref::ACADId) = @remote(b, SurfaceFromCurves([ref]))

backend_stroke_line(b::ACAD, vs) = @remote(b, PolyLine(vs))

backend_stroke_arc(b::ACAD, center::Loc, radius::Real, start_angle::Real, amplitude::Real) =
  let end_angle = start_angle + amplitude
    @remote(b, Arc(center, vz(1, center.cs), radius, start_angle, end_angle))
  end
backend_stroke_unite(b::ACAD, refs) = @remote(b, JoinCurves(refs))

realize(b::ACAD, s::EmptyShape) =
  ACADEmptyRef()
realize(b::ACAD, s::UniversalShape) =
  ACADUniversalRef()
realize(b::ACAD, s::Point) =
  @remote(b, Point(s.position))

=#

backend_line(b::PLOT, vs::Locs) =
  let pts = map(in_world, vs),
      r = PlotlyJS.scatter3d(
         x=map(cx, pts),
         y=map(cy, pts),
         z=map(cz, pts),
         line_shape="linear",
         marker_size=2,
         autocolorscale=false,
         showscale=false,
         hoverinfo="skip")
    PlotlyJS.addtraces!(connection(b), r)
    PlotNativeRef(r)
  end

realize(b::PLOT, s::Line) =
  backend_line(p, s.vertices)

realize(b::PLOT, s::Spline) = # This should be merged with opensplinepath
  if (s.v0 == false) && (s.v1 == false)
    let pts = map_division(in_world, spline_path(s.points), length(s.points)*4),
        r = PlotlyJS.scatter3d(
           x=map(cx, pts),
           y=map(cy, pts),
           z=map(cz, pts),
           line_shape="spline",
           marker_size=2,
           autocolorscale=false,
           showscale=false,
           hoverinfo="skip")
      PlotlyJS.addtraces!(connection(b), r)
      PlotNativeRef(r)
    end
  elseif (s.v0 != false) && (s.v1 != false)
    @remote(b, InterpSpline(s.points, s.v0, s.v1))
  else
    @remote(b, InterpSpline(
                     s.points,
                     s.v0 == false ? s.points[2]-s.points[1] : s.v0,
                     s.v1 == false ? s.points[end-1]-s.points[end] : s.v1))
  end

#=
realize(b::ACAD, s::ClosedSpline) =
  @remote(b, InterpClosedSpline(s.points))
realize(b::ACAD, s::Circle) =
  @remote(b, Circle(s.center, vz(1, s.center.cs), s.radius))
realize(b::ACAD, s::Arc) =
  if s.radius == 0
    @remote(b, Point(s.center))
  elseif s.amplitude == 0
    @remote(b, Point(s.center + vpol(s.radius, s.start_angle, s.center.cs)))
  elseif abs(s.amplitude) >= 2*pi
    @remote(b, Circle(s.center, vz(1, s.center.cs), s.radius))
  else
    end_angle = s.start_angle + s.amplitude
    if end_angle > s.start_angle
      @remote(b, Arc(s.center, vz(1, s.center.cs), s.radius, s.start_angle, end_angle))
    else
      @remote(b, Arc(s.center, vz(1, s.center.cs), s.radius, end_angle, s.start_angle))
    end
  end

realize(b::ACAD, s::Ellipse) =
  if s.radius_x > s.radius_y
    @remote(b, Ellipse(s.center, vz(1, s.center.cs), vxyz(s.radius_x, 0, 0, s.center.cs), s.radius_y/s.radius_x))
  else
    @remote(b, Ellipse(s.center, vz(1, s.center.cs), vxyz(0, s.radius_y, 0, s.center.cs), s.radius_x/s.radius_y))
  end
realize(b::ACAD, s::EllipticArc) =
  error("Finish this")
=#

realize(b::PLOT, s::Polygon) =
  backend_line(p, [s.vertices..., s.vertices[1]])

#=
realize(b::ACAD, s::RegularPolygon) =
  @remote(b, ClosedPolyLine(regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed)))
realize(b::ACAD, s::Rectangle) =
  @remote(b, ClosedPolyLine(
    [s.corner,
     add_x(s.corner, s.dx),
     add_xy(s.corner, s.dx, s.dy),
     add_y(s.corner, s.dy)]))
realize(b::ACAD, s::SurfaceCircle) =
  @remote(b, SurfaceCircle(s.center, vz(1, s.center.cs), s.radius))
realize(b::ACAD, s::SurfaceArc) =
    #@remote(b, SurfaceArc(s.center, vz(1, s.center.cs), s.radius, s.start_angle, s.start_angle + s.amplitude))
    if s.radius == 0
        @remote(b, Point(s.center))
    elseif s.amplitude == 0
        @remote(b, Point(s.center + vpol(s.radius, s.start_angle, s.center.cs)))
    elseif abs(s.amplitude) >= 2*pi
        @remote(b, SurfaceCircle(s.center, vz(1, s.center.cs), s.radius))
    else
        end_angle = s.start_angle + s.amplitude
        if end_angle > s.start_angle
            @remote(b, SurfaceFromCurves(
                [@remote(b, Arc(s.center, vz(1, s.center.cs), s.radius, s.start_angle, end_angle)),
                 @remote(b, PolyLine([add_pol(s.center, s.radius, end_angle),
                                              add_pol(s.center, s.radius, s.start_angle)]))]))
        else
            @remote(b, SurfaceFromCurves(
                [@remote(b, Arc(s.center, vz(1, s.center.cs), s.radius, end_angle, s.start_angle)),
                 @remote(b, PolyLine([add_pol(s.center, s.radius, s.start_angle),
                                              add_pol(s.center, s.radius, end_angle)]))]))
        end
    end

realize(b::ACAD, s::SurfaceEllipse) =
  if s.radius_x > s.radius_y
    @remote(b, SurfaceEllipse(s.center, vz(1, s.center.cs), vxyz(s.radius_x, 0, 0, s.center.cs), s.radius_y/s.radius_x))
  else
    @remote(b, SurfaceEllipse(s.center, vz(1, s.center.cs), vxyz(0, s.radius_y, 0, s.center.cs), s.radius_x/s.radius_y))
  end

realize(b::ACAD, s::SurfacePolygon) =
  @remote(b, SurfaceClosedPolyLine(s.vertices))
realize(b::ACAD, s::SurfaceRegularPolygon) =
  @remote(b, SurfaceClosedPolyLine(regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed)))
realize(b::ACAD, s::SurfaceRectangle) =
  @remote(b, SurfaceClosedPolyLine(
    [s.corner,
     add_x(s.corner, s.dx),
     add_xy(s.corner, s.dx, s.dy),
     add_y(s.corner, s.dy)]))
realize(b::ACAD, s::Surface) =
  let #ids = map(r->@remote(b, NurbSurfaceFrom(r)), @remote(b, SurfaceFromCurves(collect_ref(s.frontier))))
      ids = @remote(b, SurfaceFromCurves(collect_ref(s.frontier)))
    foreach(mark_deleted, s.frontier)
    ids
  end
backend_surface_boundary(b::ACAD, s::Shape2D) =
    map(c -> shape_from_ref(c, b), @remote(b, CurvesFromSurface(ref(s).value)))

realize(b::Plot, s::Sphere) =
  let mat = get_material(b, s)
    write_Plot_object(buffer(b), "sphere", mat, in_world(s.center), s.radius)
    void_ref(b)
  end

realize(b::Plot, s::Torus) =
  let buf = buffer(b)
    write_Plot_object(buf, "torus", get_material(b, s), s.re, s.ri) do
      let p = in_world(s.center),
          t = s.center.cs.transform
        write_Plot_object(buf, "matrix", nothing,
                            [t[1,1], t[2,1], t[3,1],
                             t[1,2], t[2,2], t[3,2],
                             t[1,3], t[2,3], t[3,3],
                             p.x, p.y, p.z])
      end
    end
    void_ref(b)
  end

#=
realize(b::ACAD, s::Cuboid) =
  ACADIrregularPyramidFrustum(connection(b), [s.b0, s.b1, s.b2, s.b3], [s.t0, s.t1, s.t2, s.t3])
realize(b::ACAD, s::RegularPyramidFrustum) =
    ACADIrregularPyramidFrustum(connection(b),
                                regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                                regular_polygon_vertices(s.edges, add_z(s.cb, s.h), s.rt, s.angle, s.inscribed))
realize(b::ACAD, s::RegularPyramid) =
  ACADIrregularPyramid(connection(b),
                          regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                          add_z(s.cb, s.h))
realize(b::ACAD, s::IrregularPyramid) =
  ACADIrregularPyramid(connection(b), s.bs, s.t)
realize(b::ACAD, s::RegularPrism) =
  let ps = regular_polygon_vertices(s.edges, s.cb, s.r, s.angle, s.inscribed)
    ACADIrregularPyramidFrustum(connection(b),
                                   ps,
                                   map(p -> add_z(p, s.h), ps))
  end
realize(b::ACAD, s::IrregularPyramidFrustum) =
    ACADIrregularPyramidFrustum(connection(b), s.bs, s.ts)

realize(b::ACAD, s::IrregularPrism) =
  write_Plot_object(buffer(b), "prism", get_material(b, s), 0, norm(s.v)) do
    s.bs,
                              map(p -> (p + s.v), s.bs))
## FIXME: deal with the rotation angle
realize(b::ACAD, s::RightCuboid) =
  ACADCenteredBox(connection(b), s.cb, s.width, s.height, s.h)
=#

realize(b::Plot, s::Box) =
  let buf = buffer(b),
      bot = in_world(s.c),
      top = in_world(s.c + vxyz(s.dx, s.dy, s.dz, s.c.cs)),
      mat = get_material(b, s)
    write_Plot_object(buf, "box", mat, bot, top)
    void_ref(b)
  end

realize(b::Plot, s::Cone) =
  let buf = buffer(b),
      bot = in_world(s.cb),
      top = in_world(s.cb + vz(s.h, s.cb.cs)),
      mat = get_material(b, s)
    write_Plot_object(buf, "cone", mat, bot, s.r, top, 0)
    void_ref(b)
  end

realize(b::Plot, s::ConeFrustum) =
  let buf = buffer(b),
      bot = in_world(s.cb),
      top = in_world(s.cb + vz(s.h, s.cb.cs)),
      mat = get_material(b, s)
    write_Plot_object(buf, "cone", mat, bot, s.rb, top, s.rt)
    void_ref(b)
  end

realize(b::Plot, s::Cylinder) =
  let buf = buffer(b),
      bot = in_world(s.cb),
      top = in_world(s.cb + vz(s.h, s.cb.cs)),
      mat = get_material(b, s)
    write_Plot_object(buf, "cylinder", mat, bot, top, s.r)
    void_ref(b)
  end

write_Plot_mesh(buf::IO, mat, points, closed_u, closed_v, smooth_u, smooth_v) =
  let si = size(points, 1),
      sj = size(points, 2),
      pts = in_world.(points),
      vcs = in_world.(add_z.(points, 1) .- points),
      idxs = quad_grid_indexes(si, sj, closed_u, closed_v)
#    let ps = reshape(permutedims(pts), :)
#      delete_all_shapes(autocad)
#      for tr in idxs
#        surface_polygon(ps[map(x->x+1,tr)], backend=autocad)
#      end
#    end
    write_Plot_object(buf, "mesh2", mat) do
      write_Plot_object(buf, "vertex_vectors", nothing, si*sj, reshape(permutedims(pts), :)...)
      # Must understand how to handle smoothness along one direction
      if smooth_u && smooth_v
        write_Plot_object(buf, "normal_vectors", nothing, si*sj, reshape(permutedims(vcs), :)...)
      end
      write_Plot_object(buf, "face_indices", nothing, length(idxs), idxs...)
      # Must understand how to handle smoothness along one direction
      #write_Plot_object(buf, "normal_indices", nothing, length(idxs), idxs...)
    end
  end
=#

realize(b::PLOT, s::SurfaceGrid) =
  let mat = 1, #get_material(b, s)
      pts = map(in_world, s.points),
      pts = [pts[i,:] for i in 1:size(pts,1)],
      r = PlotlyJS.surface(
         x=map(r->map(cx, r), pts),
         y=map(r->map(cy, r), pts),
         z=map(r->map(cz, r), pts),
         autocolorscale=false,
         showscale=false,
         hoverinfo="skip",
         colorscale=[[0, "rgb(40,40,40)"], [1, "rgb(210,210,210)"]])
    PlotlyJS.addtraces!(connection(b), r)
    PlotNativeRef(r)
  end

#=
realize(b::Plot, s::SweepPath) =
  let vertices = in_world.(path_vertices(s.profile)),
      frames = map_division(identity, s.path, 20),
      buf = buffer(b),
      mat = get_material(b, s)
    write_Plot_mesh(
      buf,
      mat,
      [xyz(cx(p), cy(p), cz(p), frame.cs) for p in vertices, frame in frames],
      is_closed_path(s.profile),
      is_closed_path(s.path),
      is_smooth_path(s.profile),
      is_smooth_path(s.path))
    void_ref(b)
  end

=#
# HACK: JUST FOR TESTING
realize(b::PLOT, s::Thicken) =
  realize(b, s.shape)

realize(b::PLOT, s::EmptyShape) = void_ref(b)

realize(b::PLOT, s::UniversalShape) = void_ref(b)

#=
#=
realize(b::Plot, s::Move) =
    let r = map_ref(s.shape) do r
                PlotMove(connection(b), r, s.v)
                r
            end
        mark_deleted(s.shape)
        r
    end

realize(b::Plot, s::Scale) =
    let r = map_ref(s.shape) do r
                PlotScale(connection(b), r, s.p, s.s)
                r
            end
        mark_deleted(s.shape)
        r
    end

realize(b::Plot, s::Rotate) =
    let r = map_ref(s.shape) do r
                PlotRotate(connection(b), r, s.p, s.v, s.angle)
                r
            end
        mark_deleted(s.shape)
        r
    end

=#

realize(b::Plot, s::UnionShape) =
  let shapes = filter(! is_empty_shape, s.shapes)
    length(shapes) == 1 ?
      (ref(shapes[1]); delete_shape(shapes[1])) :
      write_Plot_object(buffer(b), "union", get_material(b, s)) do
        for ss in shapes
          ref(ss)
          delete_shape(ss)
        end
      end
    void_ref(b)
  end

realize(b::Plot, s::IntersectionShape) =
  write_Plot_object(buffer(b), "intersection", get_material(b, s)) do
    for ss in s.shapes
      ref(ss)
      delete_shape(ss)
    end
    void_ref(b)
  end

realize(b::Plot, s::SubtractionShape3D) =
  write_Plot_object(buffer(b), "difference", get_material(b, s)) do
    ref(s.shape)
    delete_shape(s.shape)
    for ss in s.shapes
      ref(ss)
      delete_shape(ss)
    end
    void_ref(b)
  end

# BIM

realize_box(b::Plot, mat, p, dx, dy, dz) =
  let buf = buffer(b),
      bot = in_world(p),
      top = in_world(add_xyz(p, dx, dy, dz))
    write_Plot_object(buf, "box", mat, bot, top)
    void_ref(b)
  end

realize_prism(b::Plot, top, bot, side, path::PathSet, h::Real) =
  # PathSets require a different approach
  let buf = buffer(b),
      bot_vss = map(path_vertices, path.paths),
      top_vss = map(path_vertices, translate(path, vz(h)).paths)
    write_Plot_polygons(buf, bot, map(reverse, bot_vss))
    write_Plot_polygons(buf, top, top_vss)
    for (bot_vs, top_vs) in zip(bot_vss, top_vss)
      for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
        write_Plot_polygon(buf, side, vs)
      end
    end
    void_ref(b)
  end

realize_pyramid_frustum(b::Plot, top, bot, side, bot_vs::Locs, top_vs::Locs, closed=true) =
  let buf = buffer(b)
    if closed
      write_Plot_polygon(buf, bot, reverse(bot_vs))
      write_Plot_polygon(buf, top, top_vs)
    end
    for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
      write_Plot_polygon(buf, side, vs)
    end
    void_ref(b)
  end

backend_surface_polygon(b::Plot, mat, path::PathSet, acw=true) =
  acw ?
    write_Plot_polygons(buffer(b), mat, map(path_vertices, path.paths)) :
    write_Plot_polygons(buffer(b), mat, map(reverse ∘ path_vertices, path.paths))

backend_surface_polygon(b::Plot, mat, vs::Locs, acw=true) =
  let buf = buffer(b)
    write_Plot_polygon(buf, mat, acw ? vs : reverse(vs))
  end

write_Plot_polygon(io::IO, mat, vs) =
  write_Plot_object(io, "polygon", mat, length(vs)+1, vs..., vs[1])

write_Plot_polygons(io::IO, mat, vss) =
  write_Plot_object(io, "polygon", mat,
    mapreduce(length, +, vss) + length(vss),
    mapreduce(vs->[vs..., vs[1]], vcat, vss)...)

# Polygons with holes need a PathSets in Plot

subtract_paths(b::Plot, c_r_w_path::PathSet, c_l_w_path::PathSet, c_r_op_path, c_l_op_path) =
  path_set(c_r_w_path.paths..., c_r_op_path),
  path_set(c_l_w_path.paths..., c_l_op_path)

subtract_paths(b::Plot, c_r_w_path::Path, c_l_w_path::Path, c_r_op_path, c_l_op_path) =
  path_set(c_r_w_path, c_r_op_path),
  path_set(c_l_w_path, c_l_op_path)

#=
Plot families need to know the different kinds of materials
that go on each surface.
In some cases it might be the same material, but in others, such
as slabs, outside walls, etc, we will have different materials.
=#

const PlotMaterialFamily = BackendMaterialFamily{PlotMaterial}
Plot_material_family(mat::PlotMaterial) =
  PlotMaterialFamily(mat)

const PlotSlabFamily = BackendSlabFamily{PlotMaterial}
Plot_slab_family(top::PlotMaterial, bot::PlotMaterial=top, side::PlotMaterial=bot) =
  PlotSlabFamily(top, bot, side)

const PlotRoofFamily = BackendRoofFamily{PlotMaterial}
Plot_roof_family(top::PlotMaterial, bot::PlotMaterial=top, side::PlotMaterial=bot) =
  PlotRoofFamily(top, bot, side)

const PlotWallFamily = BackendWallFamily{PlotMaterial}
Plot_wall_family(right::PlotMaterial, left::PlotMaterial=right) =
  PlotWallFamily(right, left)

export Plot_material_family,
       Plot_slab_family,
       Plot_roof_family,
       Plot_wall_family,
       default_Plot_material

# Layers
current_layer(b::Plot) =
  default_Plot_material()

current_layer(layer, b::Plot) =
  default_Plot_material(layer)

backend_create_layer(b::Plot, name::String, active::Bool, color::RGB) =
  begin
    @assert active
    Plot_material(name, red=red(color), green=green(color), blue=blue(color))
  end

#=
create_ground_plane(shapes, material=default_Plot_ground_material()) =
  if shapes == []
    error("No shapes selected for analysis. Use add-Plot-shape!.")
  else
    let (p0, p1) = bounding_box(union(shapes)),
        (center, ratio) = (quad_center(p0, p1, p2, p3),
                  distance(p0, p4)/distance(p0, p2));
     ratio == 0 ?
      error("Couldn"t compute height. Use add-Plot-shape!.") :
      let pts = map(p -> intermediate_loc(center, p, ratio*10), [p0, p1, p2, p3]);
         create_surface_layer(pts, 0, ground_layer(), material)
        end
       end
  end
        w = max(floor_extra_factor()*distance(p0, p1), floor_extra_width())
        with(current_layer,floor_layer()) do
          box(xyz(min(p0.x, p1.x)-w, min(p0.y, p1.y)-w, p0.z-1-floor_distance()),
              xyz(max(p0.x, p1.x)+w, max(p0.y, p1.y)+w, p0.z-0-floor_distance()))
        end
      end
    end

=#

#=
#FIXME define the family parameters for beams
realize(b::Plot, s::Beam) =
    ref(right_cuboid(s.p0, 0.2, 0.2, s.p1, 0))
=#
realize(b::Plot, s::Union{Door, Window}) =
  void_ref(b)

used_materials(b::Plot) =
  unique(map(f -> realize(s.family, b), b.shapes))

####################################################

export export_to_jupyter
export_to_jupyter(b::Jupyter=current_backend()) =
  let camera = in_world(b.camera),
      target = in_world(b.target),
      (width, height) = render_size(),
      layout =
    PlotlyJS.Layout(autosize=true, width=width, height=height,
                    margin=attr(l=0, r=0, b=0, t=0),
                    camera_center_x=target.x,
                    camera_center_y=target.y,
                    camera_center_z=target.z,
                    camera_eye_x=camera.x,
                    camera_eye_y=camera.y,
                    camera_eye_z=camera.z,
                    scene_aspectmode="data",
                    scene_xaxis_showticklabels=false,
                    scene_yaxis_showticklabels=false,
                    scene_zaxis_showticklabels=false,
                    )
    PlotlyJS.plot([Khepri.ref(s).value for s in b.shapes if !is_empty_shape(s)], layout)
  end

=#


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

function SCIPexprtreeGetVars(tree)
    ccall((:SCIPexprtreeGetVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_EXPRTREE},), tree)
end

function SCIPexprtreeSetVars(tree, nvars, vars)
    ccall((:SCIPexprtreeSetVars, libscip), SCIP_RETCODE, (Ptr{SCIP_EXPRTREE}, Cint, Ptr{Ptr{SCIP_VAR}}), tree, nvars, vars)
end

function SCIPexprtreeAddVars(tree, nvars, vars)
    ccall((:SCIPexprtreeAddVars, libscip), SCIP_RETCODE, (Ptr{SCIP_EXPRTREE}, Cint, Ptr{Ptr{SCIP_VAR}}), tree, nvars, vars)
end

function SCIPexprtreePrintWithNames(tree, messagehdlr, file)
    ccall((:SCIPexprtreePrintWithNames, libscip), SCIP_RETCODE, (Ptr{SCIP_EXPRTREE}, Ptr{SCIP_MESSAGEHDLR}, Ptr{FILE}), tree, messagehdlr, file)
end

function SCIPexprtreeFindVar(tree, var)
    ccall((:SCIPexprtreeFindVar, libscip), Cint, (Ptr{SCIP_EXPRTREE}, Ptr{SCIP_VAR}), tree, var)
end

function SCIPnlrowGetConstant(nlrow)
    ccall((:SCIPnlrowGetConstant, libscip), Cdouble, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetNLinearVars(nlrow)
    ccall((:SCIPnlrowGetNLinearVars, libscip), Cint, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetLinearVars(nlrow)
    ccall((:SCIPnlrowGetLinearVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetLinearCoefs(nlrow)
    ccall((:SCIPnlrowGetLinearCoefs, libscip), Ptr{Cdouble}, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetNQuadVars(nlrow)
    ccall((:SCIPnlrowGetNQuadVars, libscip), Cint, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetQuadVars(nlrow)
    ccall((:SCIPnlrowGetQuadVars, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowSearchQuadVar(nlrow, var)
    ccall((:SCIPnlrowSearchQuadVar, libscip), Cint, (Ptr{SCIP_NLROW}, Ptr{SCIP_VAR}), nlrow, var)
end

function SCIPnlrowGetNQuadElems(nlrow)
    ccall((:SCIPnlrowGetNQuadElems, libscip), Cint, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetQuadElems(nlrow)
    ccall((:SCIPnlrowGetQuadElems, libscip), Ptr{SCIP_QUADELEM}, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetQuadData(nlrow, nquadvars, quadvars, nquadelems, quadelems)
    ccall((:SCIPnlrowGetQuadData, libscip), Cvoid, (Ptr{SCIP_NLROW}, Ptr{Cint}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Cint}, Ptr{Ptr{SCIP_QUADELEM}}), nlrow, nquadvars, quadvars, nquadelems, quadelems)
end

function SCIPnlrowGetExprtree(nlrow)
    ccall((:SCIPnlrowGetExprtree, libscip), Ptr{SCIP_EXPRTREE}, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetLhs(nlrow)
    ccall((:SCIPnlrowGetLhs, libscip), Cdouble, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetRhs(nlrow)
    ccall((:SCIPnlrowGetRhs, libscip), Cdouble, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetCurvature(nlrow)
    ccall((:SCIPnlrowGetCurvature, libscip), SCIP_EXPRCURV, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowSetCurvature(nlrow, curvature)
    ccall((:SCIPnlrowSetCurvature, libscip), Cvoid, (Ptr{SCIP_NLROW}, SCIP_EXPRCURV), nlrow, curvature)
end

function SCIPnlrowGetName(nlrow)
    ccall((:SCIPnlrowGetName, libscip), Cstring, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetNLPPos(nlrow)
    ccall((:SCIPnlrowGetNLPPos, libscip), Cint, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowIsInNLP(nlrow)
    ccall((:SCIPnlrowIsInNLP, libscip), UInt32, (Ptr{SCIP_NLROW},), nlrow)
end

function SCIPnlrowGetDualsol(nlrow)
    ccall((:SCIPnlrowGetDualsol, libscip), Cdouble, (Ptr{SCIP_NLROW},), nlrow)
end

ext: jl


text:
function cancel_notification(params::CancelParams, server::LanguageServerInstance, conn)
end

# TODO Provide type for params
function setTraceNotification_notification(params, server::LanguageServerInstance, conn)
end

function julia_getCurrentBlockRange_request(tdpp::VersionedTextDocumentPositionParams, server::LanguageServerInstance, conn)
    fallback = (Position(0, 0), Position(0, 0), tdpp.position)
    uri = URI2(tdpp.textDocument.uri)

    hasdocument(server, uri) || return nodocuemnt_error(uri)

    doc = getdocument(server, uri)

    if doc._version !== tdpp.version
        return mismatched_version_error(uri, doc, tdpp, "getCurrentBlockRange")
    end

    offset = get_offset(doc, tdpp.position)
    x = getcst(doc)
    loc = 0

    if typof(x) === CSTParser.FileH
        for a in x.args
            if loc <= offset <= loc + a.span
                if CSTParser.defines_module(a) # Within module at the top-level, lets see if we can select on of the block arguments
                    if loc <= offset <= loc + a.args[1].span # Within `module` keyword, so return entire expression
                        return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                    end
                    if loc + a.args[1].fullspan <= offset <= loc + a.args[1].fullspan + a.args[2].span # Within name of the module, so return entire expression
                        return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                    end

                    if loc + a.args[1].fullspan + a.args[2].fullspan <= offset <= loc + a.args[1].fullspan + a.args[2].fullspan + a.args[3].span # Within the body of a module
                        loc += a.args[1].fullspan + a.args[2].fullspan
                        for b in a.args[3].args
                            if loc <= offset <= loc + b.span
                                return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + b.span)...), Position(get_position_at(doc, loc + b.fullspan)...)
                            end
                            loc += b.fullspan
                        end
                    elseif loc + a.args[1].fullspan + a.args[2].fullspan + a.args[3].fullspan < offset < loc + a.args[1].fullspan + a.args[2].fullspan + a.args[3].fullspan + a.args[4].span # Within `end` of the module, so return entire expression
                        return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                    end
                elseif typof(a) === CSTParser.TopLevel
                    for b in a.args
                        if loc <= offset <= loc + b.span
                            return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + b.span)...), Position(get_position_at(doc, loc + b.fullspan)...)
                        end
                        loc += b.fullspan
                    end
                else
                    return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                end
            end
            loc += a.fullspan
        end
    end
    return fallback
end

function julia_activateenvironment_notification(params::String, server::LanguageServerInstance, conn)
    server.env_path = params

    trigger_symbolstore_reload(server)
end

julia_refreshLanguageServer_notification(_, server::LanguageServerInstance, conn) =
    trigger_symbolstore_reload(server)


text:
using ArcadeLearningEnvironment
using CartesianGeneticProgramming
using IICGP
using Test
using Statistics

# Global test parameters
GAME_NAMES = ["freeway", "centipede", "pong"]
N_OUT_ENCO = 2
N_STEPS = 3

function enco_cont_from_reducer(r::AbstractReducer, game_name::String)
    # Temporarily open a game to retrieve parameters
    game = Game(game_name, 0)
    n_in = 3  # RGB images
    n_out = length(getMinimalActionSet(game.ale))  # One output per legal action
    rgb = get_rgb(game)
    img_size = size(rgb[1])
    d_fitness = 1
    close!(game)

    # Encoder
    enco_nodes = [
        Node(1, 1, IICGP.CGPFunctions.f_dilate, [0.5], false)
    ]
    enco_outputs = convert(
        Array{Int16},
        ceil.((n_in + length(enco_nodes)) * rand(N_OUT_ENCO))
    )
    enco_cfg = cfg_from_info(enco_nodes, n_in, enco_outputs, IICGP.CGPFunctions,
                             d_fitness)
    enco = IPCGPInd(enco_nodes, enco_cfg, enco_outputs, img_size)

    # Forward pass to retrieve the number of input of the controller
    enco_out = CartesianGeneticProgramming.process(enco, rgb)
    features = r.reduct(enco_out, r.parameters)
    features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))

    # Controller
    cont_nodes = [
        Node(1, 2, IICGP.CGPFunctions.f_subtract, [0.5], false),
        Node(1, 2, IICGP.CGPFunctions.f_add, [0.5], false),
        Node(3, 3, IICGP.CGPFunctions.f_cos, [0.6], false)
    ]
    cont_n_in = length(features_flatten)
    cont_outputs = convert(
        Array{Int16},
        ceil.((cont_n_in + length(cont_nodes)) * rand(n_out))
    )
    cont_cfg = cfg_from_info(cont_nodes, cont_n_in, cont_outputs,
                             IICGP.CGPFunctions, d_fitness)
    cont = CGPInd(cont_nodes, cont_cfg, cont_outputs)

    enco, cont, img_size
end

@testset "Processing function" begin
    for game_name in GAME_NAMES
        # Pooling reducer
        features_size = 5
        r = PoolingReducer(Statistics.mean, features_size)
        enco, cont, img_size = enco_cont_from_reducer(r, game_name)
        game = Game(game_name, 0)
        n_out = length(getMinimalActionSet(game.ale))
        for step in 1:N_STEPS
            rgb = get_rgb(game)
            features, out = IICGP.process_f(enco, r, cont, rgb)
            @test length(features) == N_OUT_ENCO
            for i in eachindex(features)
                @test typeof(features[i]) == Array{Float64, 2}
                @test size(features[i]) == (features_size, features_size)
                @test all(f -> (0.0 <= f <= 1.0), features[i])
                @test all(f -> (!isnan(f)), features[i])
            end
            @test length(out) == n_out
        end
        close!(game)

        # Centroid reducer
        n_centroids = 20
        r = CentroidReducer(n_centroids, N_OUT_ENCO, img_size)
        enco, cont, img_size = enco_cont_from_reducer(r, game_name)
        game = Game(game_name, 0)
        for step in 1:N_STEPS
            rgb = get_rgb(game)
            features, out = IICGP.process_f(enco, r, cont, rgb)
            action = game.actions[argmax(output)]
            act(game.ale, action)
            @test length(features) == N_OUT_ENCO
            for i in eachindex(features)
                @test typeof(features[i]) == Array{Tuple{Float64,Float64},1}
                @test length(features[i]) == n_centroids
                @test all(f -> ((0.0, 0.0) <= f <= (1.0, 1.0)), features[i])
                @test all(f -> (!isnan(f[1]) && !isnan(f[2])), features[i])
            end
            @test length(out) == n_out
        end
        close!(game)
    end
end

cfg_filename = string(@__DIR__, "/dualcgp_test.yaml")
rom = "assault"

@testset "CGP controller with constant input values" begin
    mcfg, ecfg, ccfg, reducer, _ = IICGP.dualcgp_config(cfg_filename, rom)
    game = Game(rom, 0)
    enco = IICGP.IPCGPInd(ecfg)
    cont = CGPInd(ccfg)
    cstes = collect(0.0:1.0/(ccfg.n_cst_inputs-1):1.0)
    for t in 1:10
        s = get_state(game, true, true)
        output = IICGP.process(enco, reducer, cont, ccfg, s)
        @test cont.buffer[ccfg.n_in-ccfg.n_cst_inputs+1:ccfg.n_in] == cstes
    end
    close!(game)
end


text:
Base.sort(s::SubString{String}) = reduce(*, sort(split(s, "")))

function puzzle04(path::String=joinpath(@__DIR__, "..", "data/04.txt"))
    passphrases = readlines(path)
    partone = allunique.(split.(passphrases, " ")) |> sum
    parttwo = allunique.(map.(sort, split.(passphrases, " "))) |> sum
    return [partone, parttwo]
end


text:
@testset "CompactSizeUInt" begin
    @testset "Read" begin
        tests = [([0x01], 1),
                 ([0xfd, 0xd0, 0x24], 9424),
                 ([0xff, 0x70, 0x9a, 0xeb, 0xb4, 0xbb, 0x7f, 0x00, 0x00], 140444170951280)]
        for t in tests
            n = Bitcoin.read(IOBuffer(t[1]))::CompactSizeUInt
            @test n.value == t[2]
        end
    end
    @testset "Serialize" begin
        want = [0xfd, 0xfe, 0x00]
        @test Bitcoin.serialize(CompactSizeUInt(254)) == want
        want = [0xfe, 0xff, 0xff, 0xff, 0xff]
        @test Bitcoin.serialize(CompactSizeUInt(0x100000000-1)) == want
        want = [0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff]
        @test Bitcoin.serialize(CompactSizeUInt(0x10000000000000000-1)) == want
        @test_throws ErrorException CompactSizeUInt(0x10000000000000000)
    end
end


text:
export IselRobot, IselRobotParams
export initZYX, refZYX, initRefZYX, simRefZYX, prepareRobot
export setZeroPoint, setBrake, setFree, setStartStopFreq, setAcceleration
export iselErrorCodes
export readIOInput, writeIOOutput

"""Errorcodes Isel Robot """
const iselErrorCodes = Dict(
"0" => "HandShake",
"1" => "Error in Number, forbidden Character",
"2" => "Endschalterfehler, NEU Initialisieren, Neu Referenzieren",
"3" => "unzulässige Achsenzahl",
"4" => "keine Achse definiert",
"5" => "Syntax Fehler",
"6" => "Speicherende",
"7" => "unzulässige Parameterzahl",
"8" => "zu speichernder Befehl inkorrekt",
"9" => "Anlagenfehler",
"D" => "unzulässige Geschwindigkeit",
"F" => "Benutzerstop",
"G" => "ungültiges Datenfeld",
"H" => "Haubenbefehl",
"R" => "Referenzfehler",
"A" => "von dieser Steuerung nicht benutzt",
"B" => "von dieser Steuerung nicht benutzt",
"C" => "von dieser Steuerung nicht benutzt",
"E" => "von dieser Steuerung nicht benutzt",
"=" => "von dieser Steuerung nicht benutzt"
)

Base.@kwdef struct IselRobotParams <: DeviceParams
  axisRange::Vector{Vector{typeof(1.0u"mm")}} = [[0,420],[0,420],[0,420]]u"mm"
  defaultVel::Vector{typeof(1.0u"mm/s")} = [10,10,10]u"mm/s"
  defaultRefVel::Vector{typeof(1.0u"mm/s")} = [10,10,10]u"mm/s"
  invertAxes::Vector{Bool} = [false, false, false]

  minMaxVel::Vector{Int64} = [30,10000] # velocity in steps/s
  minMaxAcc::Vector{Int64} = [1,4000] # acceleration in (steps/s)/ms
  minMaxFreq::Vector{Int64} = [20,4000] # initial speed of acceleration ramp in steps/s
  stepsPermm::Float64 = 100

  serial_port::String = "/dev/ttyUSB0"
  pause_ms::Int = 200
  timeout_ms::Int = 40000
  delim_read::String = "\r"
  delim_write::String = "\r"
  baudrate::Int = 9600
  namedPositions::Dict{String,Vector{typeof(1.0u"mm")}} = Dict("origin" => [0,0,0]u"mm")
  referenceOrder::String = "zyx"
  movementOrder::String = "zyx"
  coordinateSystem::ScannerCoordinateSystem = ScannerCoordinateSystem(3)
end

IselRobotParams(dict::Dict) = params_from_dict(IselRobotParams, prepareRobotDict(dict))

Base.@kwdef mutable struct IselRobot <: Robot
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::IselRobotParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  "Current state of the robot"
  state::RobotState = INIT
  "SerialDevice for communication with the robot"
  sd::Union{SerialDevice,Nothing} = nothing
  "State variable indicating reference status"
  isReferenced::Bool = false
  "Version of the Isel Controller 1=C142; 2=newer"
  controllerVersion::Int = 1
end

abstract type IselControllerVersion end
struct IselC142 <: IselControllerVersion end
struct IseliMCS8 <: IselControllerVersion end

controllverVersion(rob::IselRobot) = rob.controllerVersion == 1 ? IselC142() : IseliMCS8()
functionUnavailable(func::AbstractString, version::IselControllerVersion) =   @error "The $func function is not available for Isel Controller version $(string(typeof(version)))"

Base.close(rob::IselRobot) = close(rob.sd.sp)

# TODO: make Isel robots with less axes possible
dof(rob::IselRobot) = 3
defaultVelocity(rob::IselRobot) = rob.params.defaultVel
axisRange(rob::IselRobot) = rob.params.axisRange
movementOrder(rob::IselRobot) = rob.params.movementOrder

function _getPosition(robot::IselRobot)
  ret = queryIsel(robot, "@0P", 19)
  checkIselError(string(ret[1]))
  pos = _parseIselPos(ret[2:19])
  return steps2mm.(pos, robot.params.stepsPermm)
end

function _setup(rob::IselRobot)
  sp = LibSerialPort.open(rob.params.serial_port, rob.params.baudrate)
  rob.sd = SerialDevice(sp, rob.params.pause_ms, rob.params.timeout_ms, rob.params.delim_read, rob.params.delim_write)

  # TODO: verify the way to identify the controller version
  if queryIsel(rob, "@0Id 1600,1600,1600,1600") == "5"
    rob.controllerVersion = 1
  else
    rob.controllerVersion = 2
  end

  initAxes(rob, 3)
  _setup(rob, controllverVersion(rob))
end
function _setup(rob::IselRobot, version::IseliMCS8)
  invertAxes(rob, rob.params.invertAxes) # only with newer version of controller
  setRefVelocity(rob, rob.params.defaultRefVel)
end
function _setup(rob::IselRobot, version::IselC142)
  _setMotorCurrent(rob, false)
end

function _enable(robot::IselRobot)
  writeIOOutput(robot, ones(Bool, 8))
  _enable(robot, controllverVersion(robot))
end
function _enable(robot::IselRobot, version::IseliMCS8)
  # NOP
end
function _enable(robot::IselRobot, version::IselC142)
  _setMotorCurrent(robot, true)
end

function _disable(robot::IselRobot)
  writeIOOutput(robot, zeros(Bool, 8))
  _disable(robot, controllverVersion(robot))
end
function _disable(robot::IselRobot, version::IseliMCS8)
  # NOP
end
function _disable(robot::IselRobot, version::IselC142)
  _setMotorCurrent(robot, false)
end

function prepareReferenceDrive(rob::IselRobot, version::IselC142)
  _moveRel(rob, [0.5u"mm", 0.5u"mm", 0.5u"mm"], nothing)
end
function prepareReferenceDrive(rob::IselRobot, version::IseliMCS8)
  # NOP
end

function _doReferenceDrive(rob::IselRobot)
  # Minor shift for not hitting the limit switch
  prepareReferenceDrive(rob, controllverVersion(rob))

  tempTimeout = rob.sd.timeout_ms
  try
    rob.sd.timeout_ms = 180000
    refAxis(rob, params(rob).referenceOrder[1])
    refAxis(rob, params(rob).referenceOrder[2])
    refAxis(rob, params(rob).referenceOrder[3])
  finally
    rob.sd.timeout_ms = tempTimeout
  end
  rob.isReferenced = true
end

function _isReferenced(robot::IselRobot)
  return robot.isReferenced && (robot.state != MPIMeasurements.ERROR)
  # TODO: find out if there is way to correctly identify the reference status

  #   currPos = getPosition(robot)
  #   currPos[1] += 0.01u"mm"
  # # need to add 0.01mm, otherwise moveAbs returns 0 although it is no longer referenced
  #   moveRes = _moveAbs(robot, currPos, nothing)
  #   if moveRes == "0"
  #     return true
  #   elseif moveRes == "R"
  #     return false
  #   elseif moveRes == "2"
  #     return false
  #   else
  #     error("Not expected \"$(moveRes)\" feedback from robot")
  #   end

  #   return false
end

function _reset(rob::IselRobot)
  close(rob)
  rob.isReferenced = false
end

function _moveAbs(rob::IselRobot, pos::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
  # for z-axis two steps and velocities are needed, compare documentation
  # set second z steps to zero
  steps = mm2steps.(pos, rob.params.stepsPermm)
  if speed === nothing
    speed = defaultVelocity(rob)
  end
  vel = mm2steps.(speed, rob.params.stepsPermm)
  if all(rob.params.minMaxVel[1] .<= vel .<= rob.params.minMaxVel[2])
    cmd = string("@0M", " ", steps[1], ",", vel[1], ",", steps[2], ",", vel[2], ",", steps[3], ",", vel[3], ",", 0, ",", 30)
    ret = queryIsel(rob, cmd)
    checkIselError(ret)
  else
    error("Velocities set not in the range of $(steps2mm.(rob.params.minMaxVel, rob.params.stepsPermm)/u"s"), you are trying to set: $speed")
  end
end

function _moveRel(rob::IselRobot, dist::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
  # for z-axis two steps and velocities are needed, compare documentation
  # set second z steps to zero
  steps = mm2steps.(dist, rob.params.stepsPermm)
  if speed === nothing
    speed = defaultVelocity(rob)
  end
  vel = mm2steps.(speed, rob.params.stepsPermm)

  if all(rob.params.minMaxVel[1] .<= vel .<= rob.params.minMaxVel[2])
    cmd = string("@0A"," ",steps[1],",",vel[1], ",",steps[2],",",vel[2], ",",steps[3],",",vel[3], ",",0,",",30)
    ret = queryIsel(rob, cmd)
    checkIselError(ret)
  else
    error("Velocities set not in the range of $(steps2mm.(rob.params.minMaxVel, rob.params.stepsPermm)/u"s"), you are trying to set: $speed")
  end
end

macro minimumISELversion(version::Int)
  return esc(quote
    if rob.controllerVersion < $version
        @error "The desired function $(var"#self#") is not available for ISEL version $(rob.controllerVersion), the minimum version is $($version)"
        return nothing
    end
  end)
end

function _setMotorCurrent(rob::IselRobot, power::Bool)
  ret = queryIsel(rob, string("@0B65529,", Int(!power)))
  checkIselError(ret)
end

function _parseIselPos(ret::AbstractString)
  # 18 hex values, 6 digits per Axis order XYZ
  xPos = reinterpret(Int32, parse(UInt32, string("0x", ret[1:6])) << 8) >> 8
  yPos = reinterpret(Int32, parse(UInt32, string("0x", ret[7:12])) << 8) >> 8
  zPos = reinterpret(Int32, parse(UInt32, string("0x", ret[13:18])) << 8) >> 8
  return [xPos,yPos,zPos]
end

function mm2steps(len::Unitful.Length, stepsPermm::Real)
  temp = round(ustrip(u"mm", len), digits=1) # round to 100um due to step error after setting powerless
  return round(Int64, temp * stepsPermm)
end

function mm2steps(len::Unitful.Velocity, stepsPermm::Real)
  temp = round(ustrip(u"mm/s", len), digits=1) # round to 100um due to step error after setting powerless
  return round(Int64, temp * stepsPermm)
end

steps2mm(steps::Integer, stepsPermm::Real) = Int64(steps) / stepsPermm * u"mm"

""" Sets the Reference velocities of the axes x,y,z """
function setRefVelocity(rob::IselRobot, vel::Vector{<:Unitful.Velocity})
  vel = mm2steps.(vel, rob.params.stepsPermm)
  minVel = rob.params.minMaxVel[1]
  maxVel = rob.params.minMaxVel[2]

  if all(minVel .<= vel .<= maxVel)
    if rob.controllerVersion < 2
      cmd = string("@0d", vel[1], ",", vel[2], ",", vel[3], ",", vel[3])
    else
      # TODO: check if this distinction is really necessary, both should be able to use @0d
      cmd = string("@0Id", " ", vel[1], ",", vel[2], ",", vel[3], ",", vel[3])
    end
    ret = queryIsel(rob, cmd)
    checkIselError(ret)
  else
      error("Velocities set not in the range of [$minVel,$maxVel], you are trying to set vel: $vel")
  end
end

""" Initializes all axes in order Z,Y,X """
function initAxes(robot::IselRobot, numAxes::Int=3)
  if numAxes == 1
    cmd = "@01"
  elseif numAxes == 2
    cmd = "@03"
  elseif numAxes == 3
    cmd = "@07"
  else
    error("Invalid number of axes, has to be in range 1-3")
  end

  ret = queryIsel(robot, cmd)
  checkIselError(ret)
end

function refAxis(robot::IselRobot, ax::Char)
  if ax == 'x'
    cmd = "@0R1"
  elseif ax == 'y'
    cmd = "@0R2"
  elseif ax == 'z'
    cmd = "@0R4"
  else
    error("Invalid axis, has to be one of x,y,z")
  end
  ret = queryIsel(robot, cmd)
  checkIselError(ret)
end

""" Sets the zero position for absolute moving at current axes position Z,Y,X """
function setZeroPoint(robot::IselRobot)
  ret = queryIsel(robot, "@0n7")
  checkIselError(ret)
  rob.isReferenced = true
end

""" Simulates Reference Z,Y,X """
simRefZYX(robot::IselRobot, version::IselControllerVersion) = functionUnavailable("simRefZYX", version)
function simRefZYX(robot::IselRobot, version::IseliMCS8)
  ret = queryIsel(robot, "@0N7")
  checkIselError(ret)
  rob.isReferenced = true
end
function simRefZYX(robot::IselRobot)
  simRefZYX(robot, controllverVersion(robot))
end

""" Sets Acceleration """
setAcceleration(robot::IselRobot, version::IselControllerVersion, acceleration) = functionUnavailable("setAcceleration", version)
function setAcceleration(robot::IselRobot, version::IseliMCS8, acceleration)
  if robot.params.minMaxAcc[1] <= acceleration <= robot.params.minMaxAcc[2]
    ret = queryIsel(robot, string("@0J", acceleration))
    checkIselError(ret)
  else
    error("Acceleration set not in the range of $(robot.params.minMaxAcc), you are trying to set acc: $acceleration")
  end
end
function setAcceleration(robot::IselRobot, acceleration)
  setAcceleration(robot, controllverVersion(robot), acceleration)
end

""" Sets StartStopFrequency"""
setStartStopFreq(robot::IselRobot, version::IselControllerVersion, freqeuency) = functionUnavailable("setStartStopFreq", version)
function setStartStopFreq(robot::IselRobot, version::IseliMCS8, frequency)
  if robot.params.minMaxFreq[1] <= frequency <= robot.params.minMaxFreq[2]
    ret = queryIsel(robot, string("@0j", frequency))
    checkIselError(ret)
  else
    error("Frequency set not in the range of $(robot.params.minMaxFreq), you are trying to set acc: $frequency")
  end
end
function setStartStopFreq(robot::IselRobot, frequency)
  setStartStopFreq(robot, controllverVersion(robot), frequency)
end

""" Sets brake, brake=false no current on brake , brake=true current on brake """
function setBrake(robot::IselRobot, brake::Bool)
  flag = brake ? 1 : 0
  ret = queryIsel(robot, string("@0g", flag))
  checkIselError(ret)
end

""" Sets free, Freifahren axis, wenn Achse über den Referenzpunkt gefahren ist"""
setFree(robot::IselRobot, version::IselControllerVersion, axis) = functionUnavailable("setFree", version)
function setFree(robot::IselRobot, version::IseliMCS8, axis)
  ret = queryIsel(robot,  string("@0F", axis))
  checkIselError(ret)
end
function setFree(robot::IselRobot, axis)
  setFree(robot, controllverVersion(robot), axis)
end

invertAxes(robot::IselRobot, version::IselControllerVersion, axes) = functionUnavailable("invertAxes", version)
function invertAxes(robot::IselRobot, version::IseliMCS8, axes::Array{Bool})
  num = 0
  for i in 1:dof(robot)
    num += axes[i] * 2^(i - 1)
  end
  ret = queryIsel(robot, string("@0ID", num))
  checkIselError(ret)
end
function invertAxes(robot::IselRobot, axes::Array{Bool})
  invertAxes(robot, controllverVersion(robot), axes)
end

""" Inverts the axes for y,z """
invertAxesYZ(robot::IselRobot) = invertAxes(robot, [false, true, true])

""" Inverts the axis for z """
invertAxisZ(robot::IselRobot) = invertAxes(robot, [false, false, true])

function checkIselError(ret::AbstractString)
  if ret != "0"
    error("Command failed: ", iselErrorCodes[ret])
  end
  return nothing
end

""" queryIsel(sd::SerialDevice,cmd::String) """
function queryIsel(rob::IselRobot, cmd::String, byteLength=1)
  sd = rob.sd
  @debug "queryIsel: " cmd

  flush(sd.sp)
  send(sd, cmd)
  i, c = LibSerialPort.sp_blocking_read(sd.sp.ref, byteLength, sd.timeout_ms)
  if i != byteLength
    error("Isel Robot did not respond!")
  end
  out = String(c)
  flush(sd.sp)
  return out
end

""" `readIOInput(robot::IselRobot)` returns an `Array{Bool,1}` of the 1-8 input slots"""
function readIOInput(robot::IselRobot)
  ret = queryIsel(robot, "@0b0", 3)
  checkIselError(string(ret[1]))
  inputBin = bin(parse(UInt8, string("0x", ret[2:3])))
  inputArray = map(x -> parse(Int64, x), collect(inputBin))
  return convert(Array{Bool,1}, inputArray)
end

""" `readIOInput(robot::IselRobot,input::Int64)` returns an Bool for the `input` slot"""
function readIOInput(robot::IselRobot, input::Int64)
  if 1 <= input && input <= 8
    return readIOInput(robot)[input]
  else
    error("input: $(input) needs to be between 1-8")
  end
end

function _writeIOOutput(robot::IselRobot, output::String)
  cmd = string("@0B0,", output)
  ret = queryIsel(robot, cmd)
  checkIselError(ret)
end

""" `writeIOOutput(robot::IselRobot,output::Array{Bool,1})` output represents 1-8 output slots"""
function writeIOOutput(robot::IselRobot, output::Array{Bool,1})
  outputInt = convert(Array{Int64,1}, output)
  outputStrings = map(x -> string(x), outputInt)
  outputBin = parse(UInt8, string("0b", string(outputStrings...)))
  _writeIOOutput(robot, string(outputBin))
end


text:
module Sundials

if isfile(joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Sundials not properly installed. Please run Pkg.build(\"Sundials\")")
end

##################################################################
# Deprecations
##################################################################

@deprecate ode cvode
@deprecate ida idasol

##################################################################
#
# Read in the wrapping code generated by the Clang.jl package.
#
##################################################################

typealias __builtin_va_list Ptr{:Void}

if isdefined(:libsundials_cvodes)
    libsundials_cvode = libsundials_cvodes
    libsundials_ida = libsundials_idas
end

shlib = libsundials_nvecserial
include("nvector.jl")
shlib = libsundials_cvode
include("libsundials.jl")
include("cvode.jl")
if isdefined(:libsundials_cvodes)
    shlib = libsundials_cvodes
    include("cvodes.jl")
end
shlib = libsundials_ida
include("ida.jl")
if isdefined(:libsundials_cvodes)
    shlib = libsundials_idas
    include("idas.jl")
end
shlib = libsundials_kinsol
include("kinsol.jl")

include("constants.jl")

##################################################################
#
# Methods to convert between Julia Vectors and Sundials N_Vectors.
#
##################################################################

nvlength(x::N_Vector) = unsafe_load(unsafe_load(convert(Ptr{Ptr{Clong}}, x)))
asarray(x::N_Vector) = pointer_to_array(N_VGetArrayPointer_Serial(x), (nvlength(x),))
asarray(x::Vector{realtype}) = x
asarray(x::Ptr{realtype}, dims::Tuple) = pointer_to_array(x, dims)
asarray(x::N_Vector, dims::Tuple) = reinterpret(realtype, asarray(x), dims)
nvector(x::Vector{realtype}) = N_VMake_Serial(length(x), x)
nvector(x::N_Vector) = x

##################################################################
#
# Methods following the C API that allow the direct use of Julia
# Vectors instead of Sundials N_Vectors and Functions in place of
# CFunctions.
#
##################################################################

# KINSOL
KINInit(mem, sysfn::Function, y) =
    KINInit(mem, cfunction(sysfn, Int32, (N_Vector, N_Vector, Ptr{Void})), nvector(y))
KINSetConstraints(mem, constraints::Vector{realtype}) =
    KINSetConstraints(mem, nvector(constraints))
KINSol(mem, u::Vector{realtype}, strategy, u_scale::Vector{realtype}, f_scale::Vector{realtype}) =
    KINSol(mem, nvector(u), strategy, nvector(u_scale), nvector(f_scale))

# IDA
IDAInit(mem, res::Function, t0, yy0, yp0) =
    IDAInit(mem, cfunction(res, Int32, (realtype, N_Vector, N_Vector, N_Vector, Ptr{Void})), t0, nvector(yy0), nvector(yp0))
IDARootInit(mem, nrtfn, g::Function) =
    IDARootInit(mem, nrtfn, cfunction(g, Int32, (realtype, N_Vector, N_Vector, Ptr{realtype}, Ptr{Void})))
IDASVtolerances(mem, reltol, abstol::Vector{realtype}) =
    IDASVtolerances(mem, reltol, nvector(abstol))
IDADlsSetDenseJacFn(mem, jac::Function) =
    IDADlsSetDenseJacFn(mem, cfunction(jac, Int32, (Int32, realtype, realtype, N_Vector, N_Vector, N_Vector, DlsMat, Ptr{Void}, N_Vector, N_Vector, N_Vector)))
IDASetId(mem, id::Vector{realtype}) =
    IDASetId(mem, nvector(id))
IDASetConstraints(mem, constraints::Vector{realtype}) =
    IDASetConstraints(mem, nvector(constraints))
IDASolve(mem, tout, tret, yret::Vector{realtype}, ypret::Vector{realtype}, itask) =
    IDASolve(mem, tout, tret, nvector(yret), nvector(ypret), itask)

# CVODE
CVodeInit(mem, f::Function, t0, y0) =
    CVodeInit(mem, cfunction(f, Int32, (realtype, N_Vector, N_Vector, Ptr{Void})), t0, nvector(y0))
CVodeReInit(mem, t0, y0::Vector{realtype}) =
    CVodeReInit(mem, t0, nvector(y0))
CVodeSVtolerances(mem, reltol, abstol::Vector{realtype}) =
    CVodeSVtolerances(mem, reltol, nvector(abstol))
CVodeGetDky(mem, t, k, dky::Vector{realtype}) =
    CVodeGetDky(mem, t, k, nvector(dky))
CVodeGetErrWeights(mem, eweight::Vector{realtype}) =
    CVodeGetErrWeights(mem, nvector(eweight))
CVodeGetEstLocalErrors(mem, ele::Vector{realtype}) =
    CVodeGetEstLocalErrors(mem, nvector(ele))
CVodeRootInit(mem, nrtfn, g::Function) =
    CVodeRootInit(mem, nrtfn, cfunction(g, Int32, (realtype, N_Vector, Ptr{realtype}, Ptr{Void})))
CVDlsSetDenseJacFn(mem, jac::Function) =
    CVDlsSetDenseJacFn(mem, cfunction(jac, Int32, (Int32, realtype, N_Vector, N_Vector, DlsMat, Ptr{Void}, N_Vector, N_Vector, N_Vector)))
CVode(mem, tout, yout::Vector{realtype}, tret, itask) =
    CVode(mem, tout, nvector(yout), tret, itask)

if isdefined(:libsundials_cvodes)
# CVODES
CVodeQuadInit(mem, fQ::Function, yQ0) =
    CVodeQuadInit(mem, cfunction(fQ, Int32, (realtype, N_Vector, N_Vector, Ptr{Void})), nvector(yQ0))
CVodeQuadReInit(mem, yQ0::Vector{realtype}) =
    CVodeQuadReInit(mem, nvector(yQ0))
CVodeQuadSVtolerances(mem, reltolQ, abstolQ::Vector{realtype}) =
    CVodeQuadSVtolerances(mem, reltolQ, nvector(abstolQ))
CVodeSensInit(mem, Ns, ism, fS::Function, yS0) =
    CVodeSensInit(mem, Ns, ism, cfunction(fS, Int32, (int32, realtype, N_Vector, N_Vector, N_Vector, N_Vector, Ptr{Void}, N_Vector, N_Vector)), nvector(yS0))
CVodeSensInit1(mem, Ns, ism, fS1::Function, yS0) =
    CVodeSensInit1(mem, Ns, ism, cfunction(fS1, Int32, (int32, realtype, N_Vector, N_Vector, int32, N_Vector, N_Vector, Ptr{Void}, N_Vector, N_Vector)), nvector(yS0))
CVodeSensReInit(mem, ism, yS0::Vector{realtype}) =
    CVodeSensReInit(mem, ism, nvector(yS0))
CVodeSensSVtolerances(mem, reltolS, abstolS::Vector{realtype}) =
    CVodeSensSVtolerances(mem, reltolS, nvector(abstolS))
CVodeQuadSensInit(mem, fQS::Function, yQS0) =
    CVodeQuadSensInit(mem, cfunction(fQS, Int32, (int32, realtype, N_Vector, N_Vector, N_Vector, N_Vector, Ptr{Void}, N_Vector, N_Vector)), nvector(yQS0))
CVodeQuadSensReInit(mem, yQS0::Vector{realtype}) =
    CVodeQuadSensReInit(mem, nvector(yQS0))
CVodeQuadSensSVtolerances(mem, reltolQS, abstolQS::Vector{realtype}) =
    CVodeQuadSensSVtolerances(mem, reltolQS, nvector(abstolQS))
CVodeGetQuad(mem, tret, yQout::Vector{realtype}) =
    CVodeGetQuad(mem, tret, nvector(yQout))
CVodeGetQuadDky(mem, t, k, dky::Vector{realtype}) =
    CVodeGetQuadDky(mem, t, k, nvector(dky))
CVodeGetSens(mem, tret, ySout::Vector{realtype}) =
    CVodeGetSens(mem, tret, nvector(ySout))
CVodeGetSens1(mem, tret, is, ySout::Vector{realtype}) =
    CVodeGetSens1(mem, tret, is, nvector(ySout))
CVodeGetSensDky(mem, t, k, dkyA::Vector{realtype}) =
    CVodeGetSensDky(mem, t, k, nvector(dkyA))
CVodeGetSensDky1(mem, t, k, is, dky::Vector{realtype}) =
    CVodeGetSensDky1(mem, t, k, is, nvector(dky))
CVodeGetQuadSens(mem, tret, yQSout::Vector{realtype}) =
    CVodeGetQuadSens(mem, tret, nvector(yQSout))
CVodeGetQuadSens1(mem, tret, is, yQSout::Vector{realtype}) =
    CVodeGetQuadSens1(mem, tret, is, nvector(yQSout))
CVodeGetQuadSensDky(mem, t, k, kdyQS_all::Vector{realtype}) =
    CVodeGetQuadSensDky(mem, t, k, nvector(kdyQS_all))
CVodeGetQuadSensDky1(mem, t, k, is, kdyQS::Vector{realtype}) =
    CVodeGetQuadSensDky1(mem, t, k, is, nvector(kdyQS))
CVodeGetQuadErrWeights(mem, eQweight::Vector{realtype}) =
    CVodeGetQuadErrWeights(mem, nvector(eQweight))
CVodeGetSensErrWeights(mem, eSweight::Vector{realtype}) =
    CVodeGetSensErrWeights(mem, nvector(eSweight))
CVodeGetQuadSensErrWeights(mem, eQSweight::Vector{realtype}) =
    CVodeGetQuadSensErrWeights(mem, nvector(eQSweight))
CVodeInitB(mem, which, fB::Function, tB0, yB0) =
    CVodeInitB(mem, which, cfunction(fB, Int32, (realtype, N_Vector, N_Vector, N_Vector, Ptr{Void})), tB0, nvector(yB0))
CVodeInitBS(mem, which, fBs::Function, tB0, yB0) =
    CVodeInitBS(mem, which, cfunction(fBs, Int32, (realtype, N_Vector, N_Vector, N_Vector, Ptr{Void})), tB0, nvector(yB0))
CVodeReInitB(mem, which, tB0, yB0::Vector{realtype}) =
    CVodeReInitB(mem, which, tB0, nvector(yB0))
CVodeSVtolerancesB(mem, which, reltolB, abstolB::Vector{realtype}) =
    CVodeSVtolerancesB(mem, which, reltolB, nvector(abstolB))
CVodeQuadInitB(mem, which, fQB::Function, yQB0) =
    CVodeQuadInitB(mem, which, cfunction(fQB, Int32, (realtype, N_Vector, N_Vector, N_Vector, Ptr{Void})), nvector(yQB0))
CVodeQuadInitBS(mem, which, fQBs::Function, yQB0) =
    CVodeQuadInitBS(mem, which, cfunction(fQBs, Int32, (realtype, N_Vector, N_Vector, N_Vector, N_Vector, Ptr{Void})), nvector(yQB0))
CVodeQuadReInitB(mem, which, yQB0::Vector{realtype}) =
    CVodeQuadReInitB(mem, which, nvector(yQB0))
CVodeQuadSVtolerancesB(mem, which, reltolQB, abstolQB::Vector{realtype}) =
    CVodeQuadSVtolerancesB(mem, which, reltolQB, nvector(abstolQB))
CVodeF(mem, tout, yout::Vector{realtype}, tret, itask, ncheckPtr) =
    CVodeF(mem, tout, nvector(yout), tret, itask, ncheckPtr)
CVodeGetB(mem, which, tBret, yB::Vector{realtype}) =
    CVodeGetB(mem, which, tBret, nvector(yB))
CVodeGetQuadB(mem, which, tBret, qB::Vector{realtype}) =
    CVodeGetQuadB(mem, which, tBret, nvector(qB))
CVodeGetAdjY(mem, which, t, y::Vector{realtype}) =
    CVodeGetAdjY(mem, which, t, nvector(y))
CVodeGetAdjDataPointHermite(mem, which, t, y::Vector{realtype}, yd::Vector{realtype}) =
    CVodeGetAdjDataPointHermite(mem, which, t, nvector(y), nvector(yd))
CVodeGetAdjDataPointPolynomial(mem, which, t, y::Vector{realtype}) =
    CVodeGetAdjDataPointPolynomial(mem, which, t, nvector(y))
CVodeWFtolerances(mem, efun::Function) =
    CVodeWFtolerances(mem, cfunction(efun, Int32, (N_Vector, N_Vector, Ptr{Void})))
CVodeSetErrHandlerFn(mem, ehfun::Function, eh_data) =
    CVodeSetErrHandlerFn(mem, cfunction(ehfun, Void, (Int32, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Void})), eh_data)

# IDAS  (still incomplete)
IDAReInit(mem, t0, yy0::Vector{realtype}, yp0::Vector{realtype}) =
    IDAReInit(mem, t0, nvector(yy0), nvector(yp0))
IDAQuadInit(mem, rhsQ::Function, yQ0) =
    IDAQuadInit(mem, cfunction(rhsQ, Int32, (realtype, N_Vector, N_Vector, N_Vector, Ptr{Void})), nvector(yQ0))
IDAQuadReInit(mem, yQ0::Vector{realtype}) =
    IDAQuadReInit(mem, nvector(yQ0))

    ## IDAQuadSVtolerances(mem, reltol, abstol::Vector{realtype}) =
    ## IDAQuadSVtolerances(mem, reltol, nvector(abstol))
end

##################################################################
#
# Simplified interfaces.
#
##################################################################

@c Int32 KINSetUserData (Ptr{:None},Any) libsundials_kinsol  ## needed to allow passing a Function through the user data

function kinsolfun(y::N_Vector, fy::N_Vector, userfun::Function)
    y = asarray(y)
    fy = asarray(fy)
    userfun(y, fy)
    return int32(0)
end

function kinsol(f::Function, y0::Vector{Float64})
    # f, Function to be optimized of the form f(y::Vector{Float64}, fy::Vector{Float64})
    #    where `y` is the input vector, and `fy` is the result of the function
    # y0, Vector of initial values
    # return: the solution vector
    neq = length(y0)
    kmem = KINCreate()
    # use the user_data field to pass a function
    #   see: https://github.com/JuliaLang/julia/issues/2554
    flag = KINInit(kmem, cfunction(kinsolfun, Int32, (N_Vector, N_Vector, Function)), nvector(y0))
    flag = KINDense(kmem, neq)
    flag = KINSetUserData(kmem, f)
    ## Solve problem
    scale = ones(neq)
    strategy = 0   # KIN_NONE
    y = copy(y0)
    flag = Sundials.KINSol(kmem,
                           y,
                           strategy,
                           scale,
                           scale)
    if flag != 0
        println("KINSol error found")
    end
    return y
end

@c Int32 CVodeSetUserData (Ptr{:None},Any) libsundials_cvode  ## needed to allow passing a Function through the user data

function cvodefun(t::Float64, y::N_Vector, yp::N_Vector, userfun::Function)
    y = Sundials.asarray(y)
    yp = Sundials.asarray(yp)
    userfun(t, y, yp)
    return int32(0)
end

function cvode(f::Function, y0::Vector{Float64}, t::Vector{Float64}; reltol::Float64=1e-4, abstol::Float64=1e-6)
    # f, Function to be optimized of the form f(y::Vector{Float64}, fy::Vector{Float64}, t::Float64)
    #    where `y` is the input vector, and `fy` is the
    # y0, Vector of initial values
    # t, Vector of time values at which to record integration results
    # reltol, Relative Tolerance to be used (default=1e-4)
    # abstol, Absolute Tolerance to be used (default=1e-6)
    # return: a solution matrix with time steps in `t` along rows and
    #         state variable `y` along columns
    neq = length(y0)
    mem = CVodeCreate(CV_BDF, CV_NEWTON)
    flag = CVodeInit(mem, cfunction(cvodefun, Int32, (realtype, N_Vector, N_Vector, Function)), t[1], nvector(y0))
    flag = CVodeSetUserData(mem, f)
    flag = CVodeSStolerances(mem, reltol, abstol)
    flag = CVDense(mem, neq)
    yres = zeros(length(t), length(y0))
    yres[1,:] = y0
    y = copy(y0)
    tout = [0.0]
    for k in 2:length(t)
        flag = CVode(mem, t[k], y, tout, CV_NORMAL)
        yres[k,:] = y
    end
    return yres
end

@c Int32 IDASetUserData (Ptr{:None},Any) libsundials_ida  ## needed to allow passing a Function through the user data

function idasolfun(t::Float64, y::N_Vector, yp::N_Vector, r::N_Vector, userfun::Function)
    y = Sundials.asarray(y)
    yp = Sundials.asarray(yp)
    r = Sundials.asarray(r)
    userfun(t, y, yp, r)
    return int32(0)   # indicates normal return
end

function idasol(f::Function, y0::Vector{Float64}, yp0::Vector{Float64}, t::Vector{Float64}; reltol::Float64=1e-4, abstol::Float64=1e-6)
    # f, Function to be optimized of the form f(y::Vector{Float64}, fy::Vector{Float64})
    #    where `y` is the input vector, and `fy` is the
    # y0, Vector of initial values
    # yp0, Vector of initial values of the derivatives
    # reltol, Relative Tolerance to be used (default=1e-4)
    # abstol, Absolute Tolerance to be used (default=1e-6)
    # return: (y,yp) two solution matrices representing the states and state derivatives
    #         with time steps in `t` along rows and state variable `y` or `yp` along columns
    neq = length(y0)
    mem = IDACreate()
    flag = IDAInit(mem, cfunction(idasolfun, Int32, (realtype, N_Vector, N_Vector, N_Vector, Function)), t[1], nvector(y0), nvector(yp0))
    flag = IDASetUserData(mem, f)
    flag = IDASStolerances(mem, reltol, abstol)
    flag = IDADense(mem, neq)
    rtest = zeros(neq)
    f(t[1], y0, yp0, rtest)
    if any(abs(rtest) .>= reltol)
        flag = IDACalcIC(mem, Sundials.IDA_YA_YDP_INIT, t[1] + tstep)
    end
    yres = zeros(length(t), length(y0))
    ypres = zeros(length(t), length(y0))
    yres[1,:] = y0
    ypres[1,:] = yp0
    y = copy(y0)
    yp = copy(yp0)
    tout = [0.0]
    for k in 2:length(t)
        retval = Sundials.IDASolve(mem, t[k], tout, y, yp, IDA_NORMAL)
        yres[k,:] = y
        ypres[k,:] = yp
    end
    return yres, ypres
end

end # module


text:
using iLQGames
import iLQGames: dx

# parameters: number of states, number of inputs, sampling time, horizon
nx, nu, ΔT, game_horizon = 4, 2, 0.1, 200

# setup the dynamics
struct Unicycle <: ControlSystem{ΔT,nx,nu} end
# the differential equation of a uncycle with state: (px, py, phi, v)
dx(cs::Unicycle, x, u, t) = SVector(x[4]cos(x[3]), x[4]sin(x[3]), u[1], u[2])
dynamics = Unicycle()

# player-1 wants the unicycle to stay close to the origin,
# player-2 wants to keep close to 1 m/s
costs = (FunctionPlayerCost((g, x, u, t) -> (x[1]^2 + x[2]^2 + u[1]^2)),
         FunctionPlayerCost((g, x, u, t) -> ((x[4] - 1)^2 + u[2]^2)))
# indices of inputs that each player controls
player_inputs = (SVector(1), SVector(2))

g = GeneralGame(game_horizon, player_inputs, dynamics, costs)

# get a solver, choose initial conditions and solve (in about 9 ms with automatic differentiation)
solver = iLQSolver(g)
x0 = SVector(1, 1, 0, 0.5)
converged, trajectory, strategies = solve(g, solver, x0)

# for visualization we need to state which state indices correspond to px and py
position_indices = tuple(SVector(1,2))
# Note: you can use the `plot_traj` call without @animated to get a non-animated plot instead.
@animated(plot_traj(trajectory, g, [:red, :green], player_inputs),
          1:game_horizon, "minimal_example.gif")


text:
@testset "Non-commutative Monomial" begin
    Mod.@ncpolyvar x
    X_ = constantmonomial(typeof(x))
    @test nvariables(X_) == 0
    @test isempty(X_.vars)
    @test isempty(X_.z)
    X0 = constantmonomial(x)
    for Y in [x^0, X0]
        @test nvariables(Y) == 1
        @test variables(Y)[1] == x
        @test collect(exponents(Y)) == [0]
    end
    for X in [X_, X0]
        for Y in [monomial(x), x^1, x * X, X * x]
            @test nvariables(Y) == 1
            @test variables(Y)[1] == x
            @test collect(exponents(Y)) == [1]
        end
    end
    for X in [X_, X0]
        for Y in [x^2, X * x^2, x^2 * X]
            @test nvariables(Y) == 1
            @test variables(Y)[1] == x
            @test collect(exponents(Y)) == [2]
        end
    end
    @testset "Issue #148" begin
        Mod.@ncpolyvar x y z
        m = x * y^4 * x^2
        @test 3 == @inferred degree(m, x)
        @test 4 == @inferred degree(m, y)
        @test 0 == @inferred degree(m, z)
    end
    @testset "Issue #71 of DynamicPolynomials" begin
        Mod.@ncpolyvar x y
        @test x^0 * y == y * x^0
    end
end
@testset "Non-commutative MonomialVector" begin
    Mod.@ncpolyvar x y
    X = emptymonovec(typeof(x))
    @test iszero(nvariables(X))
    @test isempty(variables(X))
    X = emptymonovec(x)
    @test nvariables(X) == 1
    @test variables(X)[1] == x
    X = monovec([y, x])
    @test X[1] == x
    @test X[2] == y
    # `variables` may return `[x, y, x, y]` if the polynomial has, e.g., the monomial `x * y * x * y`.
    X = monomials([y, x, y, x, y], 0:2)
    Y = monomials([x, y], 2)
    for M in [X, Y]
        @test variables(X)[1] == x
        @test variables(X)[2] == y
        @test variables(X)[3] == x
        @test collect(exponents(X[1])) == [2, 0, 0]
        @test collect(exponents(X[2])) == [1, 1, 0]
        @test collect(exponents(X[3])) == [0, 2, 0]
        @test collect(exponents(X[4])) == [0, 1, 1]
    end
    @test collect(exponents(X[5])) == [1, 0, 0]
    @test collect(exponents(X[6])) == [0, 1, 0]
    @test collect(exponents(X[7])) == [0, 0, 0]
    X0 = [x^3, x^2*y, x*y^2, x*y*x, y^3, y^2*x, y*x^2, y*x*y]
    X1 = monomials([x, y], 3)
    X2 = monovec(X0)
    Z = [[3,0,0,0],[2,1,0,0],[1,2,0,0],[1,1,1,0],[0,3,0,0],[0,2,1,0],[0,1,2,0],[0,1,1,1]]
    @test X1 isa AbstractVector{<:AbstractMonomial}
    @test X2 isa AbstractVector{<:AbstractMonomial}
    @test length(X1) == length(X2) == length(Z)
    @test issorted(X1, rev=true)
    @test issorted(X2, rev=true)
    @test X1 == X0
    for i in 1:length(Z)
        @test variables(X1[i])[1] == variables(X2[i])[1] == x
        @test variables(X1[i])[2] == variables(X2[i])[2] == y
        @test variables(X1[i])[3] == variables(X2[i])[3] == x
        @test variables(X1[i])[4] == variables(X2[i])[4] == y
        @test collect(exponents(X1[i])) == collect(exponents(X2[i])) == Z[i]
    end
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.
##########################################################################

## renamed "deployment policy" to "deployment"

const Depol = Deployment # not sure how to deprecate a const

@deprecate(
makeRandDeploymentPolicy(numAmbs::Int, numStations::Int; rng::AbstractRNG = GLOBAL_RNG),
makeRandDeployment(numAmbs, numStations; rng = rng))

@deprecate(
makeRandDeploymentPolicies(numAmbs::Int, numStations::Int, numDeployments::Int; rng::AbstractRNG = GLOBAL_RNG),
makeRandDeployments(numAmbs, numStations, numDeployments; rng = rng))

@deprecate(
applyDeploymentPolicy!(sim::Simulation, deployment::Deployment),
applyDeployment!(sim, deployment))

@deprecate(
simulateDeploymentPolicy!(sim::Simulation, deployment::Deployment),
simulateDeployment!(sim, deployment))

@deprecate(
simulateDeploymentPolicies!(sim::Simulation, deployments::Vector{Deployment}, f::Function; showEta::Bool = false),
simulateDeployments!(sim, deployments, f; showEta = showEta))

# write deployments to file
function writeDeploymentPoliciesFile(filename::String, deployments::Vector{Deployment}, numStations::Int)
	@warn("'writeDeploymentPoliciesFile' is deprecated, should use 'writeDeploymentsFile' instead")
	numAmbs = length(deployments[1])
	@assert(numStations >= maximum([maximum(d) for d in deployments]))
	numDeployments = length(deployments)

	miscTable = Table("miscData", ["numStations", "numDepols"]; rows = [[numStations, numDeployments]])
	deploymentPoliciesTable = Table("deploymentPolicies",
		["ambIndex", ["policy_$i stationIndex" for i = 1:numDeployments]...];
		cols = [collect(1:numAmbs), deployments...])
	writeTablesToFile(filename, [miscTable, deploymentPoliciesTable])
end

function readDeploymentPoliciesFile(filename::String)
	@warn("'readDeploymentPoliciesFile' is deprecated, should save deployment with 'writeDeploymentsFile' and load again with 'readDeploymentsFile'")
	tables = readTablesFromFile(filename)

	# get counts
	table = tables["miscData"]
	numStations = table.columns["numStations"][1]
	numDeployments = table.columns["numDepols"][1]

	# deployments
	table = tables["deploymentPolicies"]
	columns = table.columns # shorthand
	# check number of ambulances
	ambIndexCol = columns["ambIndex"]
	@assert(ambIndexCol == collect(1:length(ambIndexCol)))
	# check that table header has all deployments
	for i = 1:numDeployments
		@assert(in("policy_$i stationIndex", table.header), "Missing deployment $i")
	end
	# create deployments from table
	deployments = Vector{Deployment}(undef, numDeployments)
	for i = 1:numDeployments
		deployments[i] = columns["policy_$i stationIndex"]

		for value in deployments[i]
			@assert(in(value, 1:numStations))
		end
	end

	return deployments, numStations
end

##

@deprecate(
initSimulation(configFilename::String; allowResim::Bool = false, createBackup::Bool = true, allowWriteOutput::Bool = false, doPrint::Bool = false),
initSim(configFilename; allowResim = allowResim, createBackup = createBackup, allowWriteOutput = allowWriteOutput, doPrint = doPrint))

@deprecate(
initResimulation!(sim::Simulation),
initResim!(sim))

@deprecate(
calcBatchMeans(; values::Vector{Float} = [], times::Vector{Float} = [], batchTime::Float = nullTime, startTime::Float = nullTime, endTime::Float = nullTime),
calcBatchMeans(values, times, batchTime; startTime = startTime, endTime = endTime, rmPartialBatch = true, returnBatchSizes = true))

@deprecate(
writeStatsFiles!(sim::Simulation),
writeMiscOutputFiles(sim))

##

# renamed "response time" to "response duration"
@deprecate(
getCallResponseTimes(sim::Simulation),
getCallResponseDurations(sim))

@deprecate(
getAvgCallResponseTime(sim::Simulation; useMinutes::Bool = false),
getAvgCallResponseDuration(sim; useMinutes = useMinutes))

@deprecate(
calcBatchMeanResponseTimes(sim::Simulation; batchTime::Float = nullTime, warmUpTime::Float = nullTime, coolDownTime::Float = nullTime, rmPartialBatch::Bool = true, returnBatchSizes::Bool = true),
calcBatchMeanResponseDurations(sim; batchTime = batchTime, warmUpTime = warmUpTime, coolDownTime = coolDownTime, rmPartialBatch = rmPartialBatch, returnBatchSizes = returnBatchSizes))

@deprecate(
writeBatchMeanResponseTimesFile(filename::String, batchMeanResponseTimes::Array{Float,2}; batchTime = nullTime, startTime = nullTime, endTime = nullTime, responseTimeUnits::String = "minutes"),
writeBatchMeanResponseDurationsFile(filename, batchMeanResponseTimes; batchTime = batchTime, startTime = startTime, endTime = endTime, responseDurationUnits = responseTimeUnits))

##

@deprecate(
isSameLocation(loc1::Location, loc2::Location),
(==)(loc1, loc2))

##

@deprecate(
findNearestNodeInGrid(map::Map, grid::Grid, nodes::Vector{Node}, location::Location),
findNearestNode(map, grid, nodes, location))

##


text:
# Note: some of CodeTracking's functionality can only be tested by Revise

using CodeTracking
using Test, InteractiveUtils
# Note: ColorTypes needs to be installed, but note the intentional absence of `using ColorTypes`

isdefined(Main, :Revise) ? includet("script.jl") : include("script.jl")

@testset "CodeTracking.jl" begin
    m = first(methods(f1))
    file, line = whereis(m)
    scriptpath = normpath(joinpath(@__DIR__, "script.jl"))
    @test file == scriptpath
    @test line == 4
    trace = try
        call_throws()
    catch
        stacktrace(catch_backtrace())
    end
    @test whereis(trace[2]) == (scriptpath, 9)
    @test whereis(trace[3]) === nothing

    src, line = definition(String, m)
    @test src == chomp("""
    function f1(x, y)
        # A comment
        return x + y
    end
    """)
    @test line == 2

    m = first(methods(f2))
    src, line = definition(String, m)
    @test src == "f2(x, y) = x + y"
    @test line == 14

    m = first(methods(throws))
    src, line = definition(String, m)
    @test startswith(src, "@noinline")
    @test line == 7

    m = first(methods(multilinesig))
    src, line = definition(String, m)
    @test startswith(src, "@inline")
    @test line == 16

    info = CodeTracking.PkgFiles(Base.PkgId(CodeTracking))
    @test Base.PkgId(info) === info.id
    @test CodeTracking.basedir(info) == dirname(@__DIR__)

    io = IOBuffer()
    show(io, info)
    str = String(take!(io))
    @test startswith(str, "PkgFiles(CodeTracking [da1fd8a2-8d9e-5ec2-8556-3022fb5608a2]):\n  basedir:")

    @test pkgfiles("ColorTypes") === nothing
    @test_throws ErrorException pkgfiles("NotAPkg")

    # Test a method marked as missing
    m = @which sum(1:5)
    CodeTracking.method_info[m.sig] = missing
    @test whereis(m) == (CodeTracking.maybe_fix_path(String(m.file)), m.line)
    @test definition(m) === nothing

    # Test that definitions at the REPL work with `whereis`
    ex = Base.parse_input_line("replfunc(x) = 1"; filename="REPL[1]")
    eval(ex)
    m = first(methods(replfunc))
    @test whereis(m) == ("REPL[1]", 1)
    # Test with broken lookup
    oldlookup = CodeTracking.method_lookup_callback[]
    CodeTracking.method_lookup_callback[] = m -> error("oops")
    @test whereis(m) == ("REPL[1]", 1)
    # Test with definition(String, m)
    if isdefined(Base, :active_repl)
        hp = Base.active_repl.interface.modes[1].hist
        fstr = "__fREPL__(x::Int16) = 0"
        histidx = length(hp.history) + 1 - hp.start_idx
        ex = Base.parse_input_line(fstr; filename="REPL[$histidx]")
        f = Core.eval(Main, ex)
        push!(hp.history, fstr)
        @test definition(String, first(methods(f))) == (fstr, 1)
        pop!(hp.history)
    end
    CodeTracking.method_lookup_callback[] = oldlookup

    m = first(methods(Test.eval))
    @test occursin(Sys.STDLIB, whereis(m)[1])

    # https://github.com/JuliaDebug/JuliaInterpreter.jl/issues/150
    function f150()
        x = 1 + 1
        @info "hello"
    end
    m = first(methods(f150))
    src = Base.uncompressed_ast(m)
    idx = findfirst(lin -> String(lin.file) != @__FILE__, src.linetable)
    lin = src.linetable[idx]
    file, line = whereis(lin, m)
    @test endswith(file, String(lin.file))

    # Issues raised in #48
    m = @which(sum([1]; dims=1))
    if !isdefined(Main, :Revise)
        def = definition(String, m)
        @test def === nothing || isa(def[1], AbstractString)
        def = definition(Expr, m)
        @test def === nothing || isa(def, Expr)
    else
        def = definition(String, m)
        @test isa(def[1], AbstractString)
        def = definition(Expr, m)
        @test isa(def, Expr)
    end
end

@testset "With Revise" begin
    if isdefined(Main, :Revise)
        m = @which gcd(10, 20)
        sigs = signatures_at(Base.find_source_file(String(m.file)), m.line)
        @test !isempty(sigs)

        m = first(methods(edit))
        sigs = signatures_at(String(m.file), m.line)
        @test !isempty(sigs)
        sigs = signatures_at(Base.find_source_file(String(m.file)), m.line)
        @test !isempty(sigs)

        # issue #23
        @test !isempty(signatures_at("script.jl", 9))
    end
end

(a_34)(x::T, y::T) where {T<:Integer} = no_op_err("&", T)
(b_34)(x::T, y::T) where {T<:Integer} = no_op_err("|", T)
c_34(x::T, y::T) where {T<:Integer} = no_op_err("xor", T)

(d_34)(x::T, y::T) where {T<:Number} = x === y
(e_34)(x::T, y::T) where {T<:Real} = no_op_err("<" , T)
(f_34)(x::T, y::T) where {T<:Real} = no_op_err("<=", T)
l = @__LINE__
@testset "#34 last character" begin
    def, line = definition(String, @which d_34(1, 2))
    @test line == l - 3
    @test def == "(d_34)(x::T, y::T) where {T<:Number} = x === y"
end

function g()
    Base.@_inline_meta
    print("hello")
end
@testset "inline macros" begin
    def, line = CodeTracking.definition(String, @which g())
    @test def == """
    function g()
        Base.@_inline_meta
        print("hello")
    end"""
end


text:
# Julia wrapper for header: /usr/local/include/libavfilter/buffersink.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_buffersink_read,
    av_buffersink_read_samples,
    av_buffersink_get_frame,
    av_buffersink_get_samples

function av_buffersink_read(ctx,buf)
    ccall((:av_buffersink_read,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{Ptr{AVFilterBufferRef}}),ctx,buf)
end

function av_buffersink_read_samples(ctx,buf,nb_samples::Integer)
    ccall((:av_buffersink_read_samples,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{Ptr{AVFilterBufferRef}},Cint),ctx,buf,nb_samples)
end

function av_buffersink_get_frame(ctx,frame)
    ccall((:av_buffersink_get_frame,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{AVFrame}),ctx,frame)
end

function av_buffersink_get_samples(ctx,frame,nb_samples::Integer)
    ccall((:av_buffersink_get_samples,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{AVFrame},Cint),ctx,frame,nb_samples)
end


text:
using AE
using Test

@testset "AE.jl" begin

end


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

# 3D Lattice Boltzmann (BGK) model of a fluid.
# (http://exolete.com/lbm/)
# D3Q19 model. At each timestep, particle densities propagate
# outwards in the directions indicated in the figure. An
# equivalent 'equilibrium' density is found, and the densities
# relax towards that state, in a proportion governed by omega.
#               Iain Haslam, March 2006.
#
# Julia version: devectorized and threaded. The in-place
# circshift implementation is needed because of issue #10317.
# Despite that, GC time is ~22% of runtime. relax! and
# calc_equi! are broken out because of issue #10718. The
# threaded call form is used and the for loop form is is
# commented out in both because of issue #10527.

#using PyPlot
using Base.Threads

include("circshift.jl")

const t1 = 1/3
const t2 = 1/18
const t3 = 1/36

const prop_shifts = tuple([ 0  0  1], [ 0  0 -1], [ 0  1  0],
                    	  [ 0 -1  0], [ 1  0  0], [-1  0  0],
                     	  [ 1  1  0], [ 1 -1  0], [-1  1  0],
                     	  [-1 -1  0], [ 1  0  1], [ 1  0 -1],
                     	  [-1  0  1], [-1  0 -1], [ 0  1  1],
                     	  [ 0  1 -1], [ 0 -1  1], [ 0 -1 -1]);

const fourths = tuple([ 6, 8, 9,12,13], [ 7,10,11,14,15],
                      [ 4, 8,10,16,17], [ 5, 9,11,18,19],
                      [ 2,12,14,16,18], [ 3,13,15,17,19])

function relax!(F, UX, UY, UZ, nx, ny, nz, deltaU, t1D, t2D, t3D, sSQU, chunkid, nchunk)
    outerrange = Base.splitrange(nx, nchunk)
    for i = outerrange[chunkid]
    #@threads all for i = 1:nx
        for j = 1:ny
            for k = 1:nz
                density = 0.0
                for l = 1:size(F,4)
                    density = density + F[i,j,k,l]
                end
                fs = Array{Float64}(6)
                for l = 1:6
                    fs[l] = 0.0
                    for m = 1:5
                        fs[l] = fs[l] + F[i,j,k,fourths[l][m]]
                    end
                end
                UX[i,j,k] = (fs[1] - fs[2]) / density
                UY[i,j,k] = (fs[3] - fs[4]) / density
                UZ[i,j,k] = (fs[5] - fs[6]) / density

                if i == 1
                    UX[i,j,k] = UX[i,j,k] + deltaU # Increase inlet pressure
                end

                t1D[i,j,k] = t1 * density
                t2D[i,j,k] = t2 * density
                t3D[i,j,k] = t3 * density
                sSQU[i,j,k] = 3/2 * (UX[i,j,k]^2 + UY[i,j,k]^2 + UZ[i,j,k]^2)
            end
        end
    end
end

function calc_equi!(F, FEQ, t1D, t2D, t3D, U, UX, UY, UZ, sSQU, nx, ny, nz, omega, chunkid, nchunk)
    outerrange = Base.splitrange(nx, nchunk)
    for i = outerrange[chunkid]
    #@threads all for i = 1:nx
        #tid = threadid()
        for j = 1:ny
            for k = 1:nz

                FEQ[i,j,k,1] = t1D[i,j,k,1] * (1 - sSQU[i,j,k,1])

                # nearest neighbors
                FEQ[i,j,k,2] = t2D[i,j,k,1] * (1 + 3*UZ[i,j,k,1] + 9/2*UZ[i,j,k,1]^2 - sSQU[i,j,k,1])
                FEQ[i,j,k,3] = t2D[i,j,k,1] * (1 - 3*UZ[i,j,k,1] + 9/2*UZ[i,j,k,1]^2 - sSQU[i,j,k,1])
                FEQ[i,j,k,4] = t2D[i,j,k,1] * (1 + 3*UY[i,j,k,1] + 9/2*UY[i,j,k,1]^2 - sSQU[i,j,k,1])
                FEQ[i,j,k,5] = t2D[i,j,k,1] * (1 - 3*UY[i,j,k,1] + 9/2*UY[i,j,k,1]^2 - sSQU[i,j,k,1])
                FEQ[i,j,k,6] = t2D[i,j,k,1] * (1 + 3*UX[i,j,k,1] + 9/2*UX[i,j,k,1]^2 - sSQU[i,j,k,1])
                FEQ[i,j,k,7] = t2D[i,j,k,1] * (1 - 3*UX[i,j,k,1] + 9/2*UX[i,j,k,1]^2 - sSQU[i,j,k,1])

                U[1,chunkid]  = UX[i,j,k,1] + UY[i,j,k,1]
                U[2,chunkid]  = UX[i,j,k,1] - UY[i,j,k,1]
                U[3,chunkid]  = -UX[i,j,k,1] + UY[i,j,k,1]
                U[4,chunkid]  = -U[1,chunkid]
                U[5,chunkid]  = UX[i,j,k,1] + UZ[i,j,k,1]
                U[6,chunkid]  = UX[i,j,k,1] - UZ[i,j,k,1]
                U[7,chunkid]  = -U[6,chunkid]
                U[8,chunkid]  = -U[5,chunkid]
                U[9,chunkid]  = UY[i,j,k,1] + UZ[i,j,k,1]
                U[10,chunkid] = UY[i,j,k,1] - UZ[i,j,k,1]
                U[11,chunkid] = -U[10,chunkid]
                U[12,chunkid] = -U[9,chunkid]

                # next-nearest neighbors
                for l = 1:12
                    FEQ[i,j,k,l+7] = t3D[i,j,k,1] * (1 + 3*U[l,chunkid] + 9/2*(U[l,chunkid]^2) - sSQU[i,j,k,1])
                end

                for l = 1:19
                    F[i,j,k,l] = omega * FEQ[i,j,k,l] + (1 - omega) * F[i,j,k,l]
                end

            end
        end
    end
end

precompile(calc_equi!, (Array{Float64,4}, Array{Float64,4}, Array{Float64,3}, Array{Float64,3}, Array{Float64,3}, Array{Float64,2}, Array{Float64,3}, Array{Float64,3}, Array{Float64,3}, Array{Float64,3}, Int64, Int64, Int64, Float64))

function lbm3d(n)
    nx = n
    ny = nx
    nz = nx
    omega = 1.0
    density = 1.0

    # Implementation note: setting nchunk to nthreads() is a hack
    # to simulate the previous implementation's use of parallel regions.
    nchunk = nthreads()

    tprop = 0
    trelax = 0
    tequi = 0

    F = repeat([density/19], outer=[nx, ny, nz, 19])
    FEQ = F
    matsize = nx*ny*nz

    CI = [0:matsize:matsize*19;]'

    BOUND = Array{Float64}(nx,ny,nz)

    for i=1:nx, j=1:ny, k=1:nz
        BOUND[i,j,k] = ((i-5)^2 + (j-6)^2 + (k-7)^2) < 6
    end

    BOUND[:,:,1] = 1
    BOUND[:,1,:] = 1

    ON = find(BOUND); # matrix offset of each Occupied Node

    TO_REFLECT = [ON+CI[2] ON+CI[3] ON+CI[4] ON+CI[5] ON+CI[6] ON+CI[7] ON+CI[8] ON+CI[9] ON+CI[10] ON+CI[11] ON+CI[12] ON+CI[13] ON+CI[14] ON+CI[15] ON+CI[16] ON+CI[17] ON+CI[18] ON+CI[19]]
    REFLECTED = [ON+CI[3] ON+CI[2] ON+CI[5] ON+CI[4] ON+CI[7] ON+CI[6] ON+CI[11] ON+CI[10] ON+CI[9] ON+CI[8] ON+CI[15] ON+CI[14] ON+CI[13] ON+CI[12] ON+CI[19] ON+CI[18] ON+CI[17] ON+CI[16]]

    UX = Array{Float64}(nx,ny,nz)
    UY = Array{Float64}(nx,ny,nz)
    UZ = Array{Float64}(nx,ny,nz)
    U  = Array{Float64}(12,nchunk)
    t1D = Array{Float64}(nx,ny,nz)
    t2D = Array{Float64}(nx,ny,nz)
    t3D = Array{Float64}(nx,ny,nz)
    sSQU = Array{Float64}(nx,ny,nz)

    avu = 1
    prevavu = 1
    ts = 0
    deltaU = 1e-7
    numactivenodes = sum(1-BOUND)

    @time while (ts < 4000  &&  (1e-10 < abs((prevavu-avu)/avu)))  ||  ts < 100
        tprop += @elapsed begin
            # Propagate -- nearest and next-nearest neighbors
            for i = 2:19
                circshift3d1!(F, i, prop_shifts[i-1])
            end
        end

        # Densities bouncing back at next timestep
        BOUNCEDBACK = F[TO_REFLECT]

        trelax += @elapsed begin
            # Relax; calculate equilibrium state (FEQ) with equivalent speed and density to F
            @threads for chunk=1:nchunk
                relax!(F, UX, UY, UZ, nx, ny, nz, deltaU, t1D, t2D, t3D, sSQU, chunkid, nchunk)
            end
            for o in ON
                UX[o] = UY[o] = UZ[o] = t1D[o] = t2D[o] = t3D[o] = sSQU[o] = 0.0
            end
        end

        tequi += @elapsed begin
            # Calculate equilibrium distribution: stationary
            @threads for chunk=1:nchunk
                calc_equi!(F, FEQ, t1D, t2D, t3D, U, UX, UY, UZ, sSQU, nx, ny, nz, omega)
            end
        end

        F[REFLECTED] = BOUNCEDBACK

        prevavu = avu
        avu = sum(UX) / numactivenodes
        ts = ts + 1
        #println(avu)
    end

    println("ts: $(ts)")
    println("propagation time: $tprop")
    println("relaxation time: $trelax")
    println("equilibrium time: $tequi")

    #zcut=5
    #quiver(UY[:,:,zcut], UX[:,:,zcut])
    #xlabel("y"); ylabel("x")
    #title("Flow field at z = $(zcut)), after $(ts)")
end

@time lbm3d(36)
#ccall(:jl_threading_profile, Void, ())

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
#####
##### Converting samples
#####

# A samples collection is represented on the learning side as a (W, X, A, P, V)
# named-tuple. Each component is a `Float32` tensor whose last dimension corresponds
# to the sample index. Writing `n` the number of samples and `a` the total
# number of actions:
# - W (size 1×n) contains the samples weights
# - X (size …×n) contains the board representations
# - A (size a×n) contains the action masks (values are either 0 or 1)
# - P (size a×n) contains the recorded MCTS policies
# - V (size 1×n) contains the recorded values
# Note that the weight of a sample is computed as an increasing
# function of its `n` field.

function convert_sample(
    gspec::AbstractGameSpec,
    wp::SamplesWeighingPolicy,
    e::TrainingSample)

  if wp == CONSTANT_WEIGHT
    w = Float32[1]
  elseif wp == LOG_WEIGHT
    w = Float32[log2(e.n) + 1]
  else
    @assert wp == LINEAR_WEIGHT
    w = Float32[n]
  end
  x = GI.vectorize_state(gspec, e.s)
  a = GI.actions_mask(GI.init(gspec, e.s))
  p = zeros(size(a))
  p[a] = e.π
  v = [e.z]
  return (; w, x, a, p, v)
end

function convert_samples(
    gspec::AbstractGameSpec,
    wp::SamplesWeighingPolicy,
    es::Vector{<:TrainingSample})

  ces = [convert_sample(gspec, wp, e) for e in es]
  W = Flux.batch((e.w for e in ces))
  X = Flux.batch((e.x for e in ces))
  A = Flux.batch((e.a for e in ces))
  P = Flux.batch((e.p for e in ces))
  V = Flux.batch((e.v for e in ces))
  f32(arr) = convert(AbstractArray{Float32}, arr)
  return map(f32, (; W, X, A, P, V))
end

#####
##### Loss Function
#####

# Surprisingly, Flux does not like the following code (scalar operations):
# mse_wmean(ŷ, y, w) = sum((ŷ .- y).^2 .* w) / sum(w)
mse_wmean(ŷ, y, w) = sum((ŷ .- y) .* (ŷ .- y) .* w) / sum(w)

klloss_wmean(π̂, π, w) = -sum(π .* log.(π̂ .+ eps(eltype(π))) .* w) / sum(w)

entropy_wmean(π, w) = -sum(π .* log.(π .+ eps(eltype(π))) .* w) / sum(w)

wmean(x, w) = sum(x .* w) / sum(w)

function losses(nn, regws, params, Wmean, Hp, (W, X, A, P, V))
  # `regws` must be equal to `Network.regularized_params(nn)`
  creg = params.l2_regularization
  cinv = params.nonvalidity_penalty
  P̂, V̂, p_invalid = Network.forward_normalized(nn, X, A)
  V = V ./ params.rewards_renormalization
  V̂ = V̂ ./ params.rewards_renormalization
  Lp = klloss_wmean(P̂, P, W) - Hp
  Lv = mse_wmean(V̂, V, W)
  Lreg = iszero(creg) ?
    zero(Lv) :
    creg * sum(sum(w .* w) for w in regws)
  Linv = iszero(cinv) ?
    zero(Lv) :
    cinv * wmean(p_invalid, W)
  L = (mean(W) / Wmean) * (Lp + Lv + Lreg + Linv)
  return (L, Lp, Lv, Lreg, Linv)
end

#####
##### Trainer Utility
#####

struct Trainer
  network :: AbstractNetwork
  samples :: Vector{<:TrainingSample}
  params :: LearningParams
  data :: NamedTuple # (W, X, A, P, V) tuple obtained after converting `samples`
  Wmean :: Float32
  Hp :: Float32
  batches_stream # infinite stateful iterator of training batches
  function Trainer(gspec, network, samples, params; test_mode=false)
    if params.use_position_averaging
      samples = merge_by_state(samples)
    end
    data = convert_samples(gspec, params.samples_weighing_policy, samples)
    network = Network.copy(network, on_gpu=params.use_gpu, test_mode=test_mode)
    W, X, A, P, V = data
    Wmean = mean(W)
    Hp = entropy_wmean(P, W)
    # Create a batches stream
    batchsize = min(params.batch_size, length(W))
    batches = Flux.Data.DataLoader(data; batchsize, partial=false, shuffle=true)
    batches_stream = map(batches) do b
      Network.convert_input(network, b)
    end |> Util.cycle_iterator |> Iterators.Stateful
    return new(network, samples, params, data, Wmean, Hp, batches_stream)
  end
end

data_weights(tr::Trainer) = tr.data.W

num_samples(tr::Trainer) = length(data_weights(tr))

num_batches_total(tr::Trainer) = num_samples(tr) ÷ tr.params.batch_size

function get_trained_network(tr::Trainer)
  return Network.copy(tr.network, on_gpu=false, test_mode=true)
end

function batch_updates!(tr::Trainer, n)
  regws = Network.regularized_params(tr.network)
  L(batch...) = losses(tr.network, regws, tr.params, tr.Wmean, tr.Hp, batch)[1]
  data = Iterators.take(tr.batches_stream, n)
  ls = Vector{Float32}()
  Network.train!(tr.network, tr.params.optimiser, L, data, n) do i, l
    push!(ls, l)
  end
  Network.gc(tr.network)
  return ls
end

#####
##### Generating debugging reports
#####

function mean_learning_status(reports, ws)
  L     = wmean([r.loss.L     for r in reports], ws)
  Lp    = wmean([r.loss.Lp    for r in reports], ws)
  Lv    = wmean([r.loss.Lv    for r in reports], ws)
  Lreg  = wmean([r.loss.Lreg  for r in reports], ws)
  Linv  = wmean([r.loss.Linv  for r in reports], ws)
  Hpnet = wmean([r.Hpnet      for r in reports], ws)
  Hp    = wmean([r.Hp         for r in reports], ws)
  return Report.LearningStatus(Report.Loss(L, Lp, Lv, Lreg, Linv), Hp, Hpnet)
end

function learning_status(tr::Trainer, samples)
  # As done now, this is slighly inefficient as we solve the
  # same neural network inference problem twice
  W, X, A, P, V = samples
  regws = Network.regularized_params(tr.network)
  Ls = losses(tr.network, regws, tr.params, tr.Wmean, tr.Hp, samples)
  Ls = Network.convert_output_tuple(tr.network, Ls)
  Pnet, _ = Network.forward_normalized(tr.network, X, A)
  Hpnet = entropy_wmean(Pnet, W)
  Hpnet = Network.convert_output(tr.network, Hpnet)
  return Report.LearningStatus(Report.Loss(Ls...), tr.Hp, Hpnet)
end

function learning_status(tr::Trainer)
  batchsize = min(tr.params.loss_computation_batch_size, num_samples(tr))
  batches = Flux.Data.DataLoader(tr.data; batchsize, partial=true)
  reports = map(batches) do batch
    batch = Network.convert_input(tr.network, batch)
    return learning_status(tr, batch)
  end
  ws = [sum(batch.W) for batch in batches]
  return mean_learning_status(reports, ws)
end

function samples_report(tr::Trainer)
  status = learning_status(tr)
  # Samples in `tr.samples` can be merged by board or not
  num_samples = sum(e.n for e in tr.samples)
  num_boards = length(merge_by_state(tr.samples))
  Wtot = sum(data_weights(tr))
  return Report.Samples(num_samples, num_boards, Wtot, status)
end

function memory_report(
    mem::MemoryBuffer,
    nn::AbstractNetwork,
    learning_params::LearningParams,
    params::MemAnalysisParams
  )
  # It is important to load the neural network in test mode so as to not
  # overwrite the batch norm statistics based on biased data.
  Tr(samples) = Trainer(mem.gspec, nn, samples, learning_params, test_mode=true)
  all_samples = samples_report(Tr(get_experience(mem)))
  latest_batch = isempty(last_batch(mem)) ?
    all_samples :
    samples_report(Tr(last_batch(mem)))
  per_game_stage = begin
    es = get_experience(mem)
    sort!(es, by=(e->e.t))
    csize = ceil(Int, length(es) / params.num_game_stages)
    stages = collect(Iterators.partition(es, csize))
    map(stages) do es
      ts = [e.t for e in es]
      stats = samples_report(Tr(es))
      Report.StageSamples(minimum(ts), maximum(ts), stats)
    end
  end
  return Report.Memory(latest_batch, all_samples, per_game_stage)
end


text:
"""
    gauss(x, y)

Input: A basis {x,y} of a lattice L in R^n.
Output: A minimal basis v1, v2 of the lattice L.

The basis v1, v2 is minimal if
    - v1 is the shortest nonzero vector in L and
    - v2 is the shortest vector in L which is not a multiple of v1.

# Example 2.9 (Bremner)
```
julia> x = [-56, 43]
julia> y = [95, -73]
julia> gauss(x,y)
([1,1], [-1,2])
```
"""
function gauss(x, y)
    if norm(x) <= norm(y)
        v1, v2 = x, y
    else
        v1, v2 = y, x
    end

    finished = false
    while (!finished)
        m = round_ties_down(dot(v2, v1) / dot(v1, v1))
        v2 = v2 - m * v1

        if norm(v1) <= norm(v2)
            finished = true
        else
            v1, v2 = v2, v1
        end
    end

    return v1, v2
end


text:
using Test
using Smoothers
import Smoothers: hmaSymmetricWeights

@testset "hma" begin

    # Types' promotion
    fx = hma(Float32.(1:10),7);
    @test eltype(fx) == Float32
    fx = hma(Float64.(1:10),7);
    @test eltype(fx) == Float64
    fx = hma(BigFloat.(1:10),7);
    @test eltype(fx) == BigFloat
    fx = hma(1:10,7);
    @test eltype(fx) == Float64

    # weights
    w = hmaSymmetricWeights(8,Float64)
    @test ismissing(w[end])

    w = hmaSymmetricWeights(7,Float64)
    @test sum(w) ≈ 1.0

    w = hmaSymmetricWeights(13,Float64)
    @test sum(w) ≈ 1.0
    # (https://www.mathworks.com/help/econ/seasonal-adjustment-using-snxd7m-seasonal-filters.html)
    @test length(w) == 13
    @test round.(w,digits=3) == [-0.019, -0.028, 0.0, 0.065, 0.147, 0.214, 0.24,
                                  0.214, 0.147, 0.065, 0.0, -0.028, -0.019]
    # 13-term moving average
    x = hma(sin.(1:100), 13)
    @test length(x) == 100
    @test first(x)  ≈ 0.5636212576875559
    @test last(x)   ≈ -0.6658500507547161
    @test x[50]     ≈ -0.043655947941143455

end


text:
module LayeredVectors

export LayeredVector,
       push

struct LayeredVector{T} <: AbstractVector{T}
    length::Int
    element::T
    previous::LayeredVector{T}

    function LayeredVector{T}() where {T}
        new{T}(0)
    end

    function LayeredVector{T}(len::Integer, element, previous::LayeredVector) where {T}
        new{T}(len, element, previous)
    end
end

function LayeredVector{T}(x::AbstractVector) where {T}
    v = LayeredVector{T}()
    for element in x
        v = push(v, element)
    end
    v
end

Base.convert(::Type{LayeredVector{T}}, v::LayeredVector{T}) where {T} = v
Base.convert(::Type{LayeredVector{T}}, v::AbstractVector) where {T} = LayeredVector{T}(v)

Base.size(v::LayeredVector) = (v.length,)

function Base.iterate(v::LayeredVector, state=1)
    if state > v.length
        return nothing
    elseif state == v.length
        return (v.element, state + 1)
    else
        return iterate(v.previous, state)
    end
end

function Base.getindex(v::LayeredVector, i::Integer)
    if i > v.length
        throw(BoundsError())
    elseif i == v.length
        return v.element
    else
        return v.previous[i]
    end
end

function push(v::LayeredVector{T}, element::T) where {T}
    LayeredVector{T}(length(v) + 1, element, v)
end

function push(v::AbstractVector{T}, element::T) where {T}
    x = copy(v)
    push!(x, element)
    x
end

end


text:
using Documenter, Rebugger

makedocs(
    modules = [Rebugger],
    clean = false,
    format = :html,
    sitename = "Rebugger.jl",
    authors = "Tim Holy",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "usage.md",
        "config.md",
        "limitations.md",
        "internals.md",
        "reference.md",
    ],
    # # Use clean URLs, unless built as a "local" build
    # html_prettyurls = !("local" in ARGS),
#    html_canonical = "https://juliadocs.github.io/Rebugger.jl/stable/",
)

deploydocs(
    repo = "github.com/timholy/Rebugger.jl.git",
    target = "build",
    julia = "1.0",
    deps = nothing,
    make = nothing,
)