text:
export rhino,
       reset_rhino_connection

macro rhino_str(str)
    rpc("RH", str)
end

# We need some additional Encoders
encode_Brep = encode_Guid
decode_Brep = decode_Guid
decode_Brep_array = decode_Guid_array
encode_RhinoObject = encode_Guid
encode_RhinoObject_array = encode_Guid_array
decode_RhinoObject = decode_Guid

encode_Plane(c::IO, v::XYZ) =
  let o = in_world(v)
      vx = in_world(vx(1, c.cs))
      vy = in_world(vy(1, c.cs))
    encode_Point3d(c, o)
    encode_Vector3d(c, vx)
    encode_Vector3d(c, vy)
  end
decode_Plane(c::IO) =
    loc_from_o_vx_vy(decode_Point3d(c), decode_Vector3d(c), decode_Vector3d(c))

encode_TableFamily = encode_int
decode_TableFamily = decode_int
encode_ChairFamily = encode_int
decode_ChairFamily = decode_int
encode_TableChairFamily = encode_int
decode_TableChairFamily = decode_int

rhino"public void SetView(Point3d position, Point3d target, double lens, bool perspective, string style)"
rhino"public void View(Point3d position, Point3d target, double lens)"
rhino"public void ViewTop()"
rhino"public Point3d ViewCamera()"
rhino"public Point3d ViewTarget()"
rhino"public double ViewLens()"
#=
rhino"public byte Sync()"
rhino"public byte Disconnect()"
=#
rhino"public Guid Point(Point3d p)"
rhino"public Point3d PointPosition(Guid ent)"

rhino"public Guid PolyLine(Point3d[] pts)"
rhino"public Point3d[] LineVertices(RhinoObject id)"
rhino"public Guid Spline(Point3d[] pts)"
#=
rhino"public Entity InterpSpline(Point3d[] pts, Vector3d tan0, Vector3d tan1)"
=#
rhino"public Guid ClosedPolyLine(Point3d[] pts)"
rhino"public Guid ClosedSpline(Point3d[] pts)"
#=
rhino"public Entity InterpClosedSpline(Point3d[] pts)"
=#
rhino"public Guid Circle(Point3d c, Vector3d n, double r)"
rhino"public Point3d CircleCenter(RhinoObject obj)"
rhino"public Vector3d CircleNormal(RhinoObject obj)"
rhino"public double CircleRadius(RhinoObject obj)"
rhino"public Guid Ellipse(Point3d c, Vector3d n, double radiusX, double radiusY)"
rhino"public Guid Arc(Point3d c, Vector3d n, double radius, double startAngle, double endAngle)"
rhino"public Guid JoinCurves(Guid[] objs)"
#=

rhino"public Entity Text(string str, Point3d corner, Vector3d vx, Vector3d vy, double height)"

rhino"public Guid SurfaceFromCurve(Entity curve)"
=#
rhino"public Guid SurfaceCircle(Point3d c, Vector3d n, double r)"
rhino"public Guid SurfaceEllipse(Point3d c, Vector3d n, double radiusX, double radiusY)"
rhino"public Guid SurfaceArc(Point3d c, Vector3d n, double radius, double startAngle, double endAngle)"
rhino"public Guid SurfaceClosedPolyLine(Point3d[] pts)"
rhino"public Guid SurfaceFrom(Guid[] objs)"
rhino"public Guid Sphere(Point3d c, double r)"
rhino"public Guid Torus(Point3d c, Vector3d vz, double majorRadius, double minorRadius)"
rhino"public Brep Cylinder(Point3d bottom, double radius, Point3d top)"
rhino"public Brep Cone(Point3d bottom, double radius, Point3d top)"
rhino"public Brep ConeFrustum(Point3d bottom, double bottom_radius, Point3d top, double top_radius)"
rhino"public Brep Box(Point3d corner, Vector3d vx, Vector3d vy, double dx, double dy, double dz)"
rhino"public Brep XYCenteredBox(Point3d corner, Vector3d vx, Vector3d vy, double dx, double dy, double dz)"
#=
rhino"public Guid IrregularPyramidMesh(Point3d[] pts, Point3d apex)"
=#
rhino"public Brep IrregularPyramid(Point3d[] pts, Point3d apex)"
rhino"public Brep IrregularPyramidFrustum(Point3d[] bpts, Point3d[] tpts)"
#rhino"public Entity MeshFromGrid(int m, int n, Point3d[] pts, bool closedM, bool closedN)"
rhino"public Guid SurfaceFromGrid(int nU, int nV, Point3d[] pts, bool closedU, bool closedV, int degreeU, int degreeV)"

#=
rhino"public Entity SolidFromGrid(int m, int n, Point3d[] pts, bool closedM, bool closedN, int level, double thickness)"
=#
rhino"public Brep[] Thicken(RhinoObject obj, double thickness)"
rhino"public double[] CurveDomain(RhinoObject obj)"
rhino"public double CurveLength(RhinoObject obj)"
rhino"public Plane CurveFrameAt(RhinoObject obj, double t)"
rhino"public Plane CurveFrameAtLength(RhinoObject obj, double l)"
#=
rhino"public Guid NurbSurfaceFrom(ObjectId id)"
rhino"public double[] SurfaceDomain(Entity ent)"
rhino"public Frame3d SurfaceFrameAt(Entity ent, double u, double v)"
=#
rhino"public Brep Extrusion(RhinoObject obj, Vector3d dir)"
rhino"public Brep[] SweepPathProfile(RhinoObject path, RhinoObject profile, double rotation, double scale)"
#=
rhino"public Guid Sweep(ObjectId pathId, ObjectId profileId, double rotation, double scale)"
rhino"public Guid Loft(ObjectId[] profilesIds, ObjectId[] guidesIds, bool ruled, bool closed)"
rhino"public void Unite(ObjectId objId0, ObjectId objId1)"
=#
rhino"public Guid[] Intersect(RhinoObject obj0, RhinoObject obj1)"
rhino"public Guid[] Subtract(RhinoObject obj0, RhinoObject obj1)"
#=
rhino"public Guid Revolve(ObjectId profileId, Point3d p, Vector3d n, double startAngle, double amplitude)"
=#
rhino"public Guid Move(Guid id, Vector3d v)"
rhino"public Guid Scale(Guid id, Point3d p, double s)"
rhino"public Guid Rotate(Guid id, Point3d p, Vector3d n, double a)"
rhino"public Guid Mirror(Guid id, Point3d p, Vector3d n, bool copy)"
#=
rhino"public Point3d[] GetPoint(string prompt)"
rhino"public Point3d[] BoundingBox(ObjectId[] ids)"
rhino"public void ZoomExtents()"
rhino"public void SetSystemVariableInt(string name, int value)"
rhino"public int Render(int width, int height, string path)"
rhino"public int Command(string cmd)"
rhino"public void DisableUpdate()"
rhino"public void EnableUpdate()"

=#
rhino"public bool IsPoint(RhinoObject e)"
rhino"public bool IsCircle(RhinoObject e)"
rhino"public bool IsPolyLine(RhinoObject e)"
rhino"public bool IsSpline(RhinoObject e)"
rhino"public bool IsInterpSpline(RhinoObject e)"
rhino"public bool IsClosedPolyLine(RhinoObject e)"
rhino"public bool IsClosedSpline(RhinoObject e)"
rhino"public bool IsInterpClosedSpline(RhinoObject e)"
rhino"public bool IsEllipse(RhinoObject e)"
rhino"public bool IsArc(RhinoObject e)"
rhino"public bool IsText(RhinoObject e)"
rhino"public byte ShapeCode(RhinoObject id)"
#=
rhino"public BIMLevel FindOrCreateLevelAtElevation(double elevation)"
rhino"public BIMLevel UpperLevel(BIMLevel currentLevel, double addedElevation)"
rhino"public double GetLevelElevation(BIMLevel level)"

rhino"public FloorFamily FloorFamilyInstance(double totalThickness, double coatingThickness)"
rhino"public Entity LightweightPolyLine(Point2d[] pts, double[] angles, double elevation)"
rhino"public Entity SurfaceLightweightPolyLine(Point2d[] pts, double[] angles, double elevation)"
rhino"public ObjectId CreatePathFloor(Point2d[] pts, double[] angles, BIMLevel level, FloorFamily family)"
=#
rhino"public Guid ClosedPathCurveArray(Point3d[] pts, double[] angles)"
rhino"public Brep[] PathWall(RhinoObject obj, double thickness, double height)"
rhino"public Brep RectangularTable(Point3d c, double angle, double length, double width, double height, double top_thickness, double leg_thickness)"

#rhino"public TableFamily FindOrCreateTableFamily(double length, double width, double height, double top_thickness, double leg_thickness)"
#rhino"public ChairFamily FindOrCreateChairFamily(double length, double width, double height, double seat_height, double thickness)"
#rhino"public TableChairFamily FindOrCreateTableChairFamily(TableFamily tableFamily, ChairFamily chairFamily, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)"
#rhino"public Brep RectangularTable(Point3d c, double angle, TableFamily family)"
#rhino"public Brep Chair(Point3d c, double angle, ChairFamily family)"
#rhino"public Guid InstanceChair(Point3d c, double angle, ChairFamily family)"
#rhino"public Brep[] RowOfChairs(Point3d c, double angle, int n, double spacing, ChairFamily family)"
#rhino"public Brep[] CenteredRowOfChairs(Point3d c, double angle, int n, double spacing, ChairFamily family)"
#rhino"public Brep[] RectangularTableAndChairs(Point3d c, double angle, TableChairFamily f)"

rhino"public Brep[] BaseRectangularTable(double length, double width, double height, double top_thickness, double leg_thickness)"
rhino"public int CreateRectangularTableFamily(double length, double width, double height, double top_thickness, double leg_thickness)"
rhino"public Guid Table(Point3d c, double angle, int family)"
rhino"public Brep[] BaseChair(double length, double width, double height, double seat_height, double thickness)"
rhino"public int CreateChairFamily(double length, double width, double height, double seat_height, double thickness)"
rhino"public Guid Chair(Point3d c, double angle, int family)"
rhino"public GeometryBase[] BaseRectangularTableAndChairs(int tableFamily, int chairFamily, double tableLength, double tableWidth, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)"
rhino"public int CreateRectangularTableAndChairsFamily(int tableFamily, int chairFamily, double tableLength, double tableWidth, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)"
rhino"public Guid TableAndChairs(Point3d c, double angle, int family)"

abstract type RHKey end
const RHId = Guid
const RHRef = GenericRef{RHKey, RHId}
const RHRefs = Vector{RHRef}
const RHEmptyRef = EmptyRef{RHKey, RHId}
const RHNativeRef = NativeRef{RHKey, RHId}
const RHUnionRef = UnionRef{RHKey, RHId}
const RHSubtractionRef = SubtractionRef{RHKey, RHId}
const RH = SocketBackend{RHKey, RHId}

void_ref(b::RH) = RHNativeRef(zeros(UInt8, 16))

create_RH_connection() = create_backend_connection("Rhinoceros", 12000)

const rhino = RH(LazyParameter(TCPSocket, create_RH_connection))
reset_rhino_connection() = reset(rhino.connection)

# This should not be done automatically
backend_name(b::RH) = "Rhino"

realize(b::RH, s::EmptyShape) =
  RHEmptyRef()
realize(b::RH, s::UniversalShape) =
  RHUniversalRef()

backend_stroke(b::RH, path::CircularPath) =
    RHCircle(connection(b), path.center, vz(1, path.center.cs), path.radius)
backend_stroke(b::RH, path::RectangularPath) =
    let c = path.corner,
        dx = path.dx,
        dy = path.dy
        RHClosedPolyLine(connection(b), [c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)])
    end
backend_stroke(b::RH, path::OpenPolygonalPath) =
  	RHPolyLine(connection(b), path.vertices)
backend_stroke(b::RH, path::ClosedPolygonalPath) =
    RHClosedPolyLine(connection(b), path.vertices)

backend_stroke(b::RH, path::OpenSplinePath) =
  if (path.v0 == false) && (path.v1 == false)
    #ACADSpline(connection(b), path.vertices)
    RHSpline(connection(b), path.vertices)
  elseif (path.v0 != false) && (path.v1 != false)
    RHInterpSpline(connection(b), path.vertices, path.v0, path.v1)
  else
    RHInterpSpline(connection(b),
                     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::RH, path::ClosedSplinePath) =
    RHInterpClosedSpline(connection(b), path.vertices)
backend_fill(b::RH, path::ClosedSplinePath) =
    backend_fill_curves(b, RHInterpClosedSpline(connection(b), path.vertices))

backend_stroke_unite(b::RH, refs) = RHJoinCurves(connection(b), refs)

#=backend_fill(b::RH, path::ClosedPolygonalPath) =
    RHSurfaceClosedPolyLine(connection(b), path.vertices)
    backend_fill(b::RH, path::RectangularPath) =
        let c = path.corner,
            dx = path.dx,
            dy = path.dy
            RHSurfaceClosedPolyLine(connection(b), [c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)])
        end
backend_fill(b::RH, path::RectangularPath) =
    let c = path.corner,
        dx = path.dx,
        dy = path.dy
        SurfaceClosedPolyLine(connection(b), [c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)])
    end
=#

backend_fill_curves(b::RH, gs::Guids) = RHSurfaceFrom(connection(b), gs)
backend_fill_curves(b::RH, g::Guid) = RHSurfaceFrom(connection(b), [g])

backend_stroke_line(b::RH, vs) = RHPolyLine(connection(b), vs)

backend_stroke_arc(b::RH, center::Loc, radius::Real, start_angle::Real, amplitude::Real) =
    let end_angle = start_angle + amplitude
        if end_angle > start_angle
            RHArc(connection(b), center, vz(1, center.cs), radius, start_angle, end_angle)
        else
            RHArc(connection(b), center, vz(1, center.cs), radius, end_angle, start_angle)
        end
    end

realize(b::RH, s::Point) =
  RHPoint(connection(b), s.position)

realize(b::RH, s::Line) =
  RHPolyLine(connection(b), s.vertices)

realize(b::RH, s::Spline) =
  if (s.v0 == false) && (s.v1 == false)
    RHSpline(connection(b), s.points)
  elseif (s.v0 != false) && (s.v1 != false)
    RHInterpSpline(connection(b), s.points, s.v0, s.v1)
  else
    RHInterpSpline(connection(b),
                     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::RH, s::ClosedSpline) =
  RHInterpClosedSpline(connection(b), s.points)

realize(b::RH, s::Circle) =
  RHCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)

realize(b::RH, s::Arc) =
  if s.radius == 0
    RHPoint(connection(b), s.center)
  elseif s.amplitude == 0
    RHPoint(connection(b), s.center + vpol(s.radius, s.start_angle, s.center.cs))
  elseif abs(s.amplitude) >= 2*pi
    RHCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)
  else
    end_angle = s.start_angle + s.amplitude
    if end_angle > s.start_angle
      RHArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, end_angle)
    else
      RHArc(connection(b), s.center, vz(1, s.center.cs), s.radius, end_angle, s.start_angle)
    end
  end

realize(b::RH, s::Ellipse) =
    RHEllipse(connection(b), s.center, vz(1, s.center.cs), s.radius_x, s.radius_y)

realize(b::RH, s::EllipticArc) =
  error("Finish this")

realize(b::RH, s::Polygon) =
  RHClosedPolyLine(connection(b), s.vertices)

realize(b::RH, s::RegularPolygon) =
  RHClosedPolyLine(connection(b), regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed))

realize(b::RH, s::Rectangle) =
  RHClosedPolyLine(connection(b), [s.corner, add_x(s.corner, s.dx), add_xy(s.corner, s.dx, s.dy), add_y(s.corner, s.dy)])
#
rhino"public Point3d[] CurvePointsAt(RhinoObject obj, double[] ts)"
rhino"public Vector3d[] CurveTangentsAt(RhinoObject obj, double[] ts)"
rhino"public Vector3d[] CurveNormalsAt(RhinoObject obj, double[] ts)"

backend_map_division(b::RH, f::Function, s::Shape1D, n::Int) =
  let conn = connection(b),
      r = ref(s).value,
      (t1, t2) = RHCurveDomain(conn, r),
      ti = division(t1, t2, n),
      ps = RHCurvePointsAt(conn, r, ti),
      ts = RHCurveTangentsAt(conn, r, ti),
      #ns = ACADCurveNormalsAt(conn, r, ti),
      frames = rotation_minimizing_frames(RHCurveFrameAt(conn, r, t1), ps, ts)
    map(f, frames)
  end

rhino"public Guid SurfaceCircle(Point3d c, Vector3d n, double r)"

realize(b::RH, s::SurfaceCircle) =
  RHSurfaceCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)

realize(b::RH, s::SurfaceArc) =
  RHSurfaceArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, s.start_angle + s.amplitude)

realize(b::RH, s::SurfaceEllipticArc) =
  error("Finish this")

realize(b::RH, s::SurfaceEllipse) =
    RHSurfaceEllipse(connection(b), s.center, vz(1, s.center.cs), s.radius_x, s.radius_y)

rhino"public Guid SurfaceClosedPolyLine(Point3d[] pts)"

realize(b::RH, s::SurfacePolygon) =
  RHSurfaceClosedPolyLine(connection(b), s.vertices)
realize(b::RH, s::SurfaceRegularPolygon) =
  RHSurfaceClosedPolyLine(connection(b), regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed))
realize(b::RH, s::SurfaceRectangle) =
  let c = s.corner
    RHSurfaceClosedPolyLine(connection(b), [c, add_x(c, s.dx), add_xy(c, s.dx, s.dy), add_y(c, s.dy)])
  end
#=
realize(b::RH, s::Surface) =
  let ids = RHSurfaceFromCurves(connection(b), collect_ref(s.frontier))
    foreach(mark_deleted, s.frontier)
    ids
  end
realize(b::RH, s::Text) =
  RHText(connection(b), s.str, s.c, vx(1, s.c.cs), vy(1, s.c.cs), s.h)
=#

rhino"public double[] SurfaceDomain(RhinoObject obj)"
rhino"public Plane SurfaceFrameAt(RhinoObject obj, double u, double v)"

backend_surface_domain(b::RH, s::Shape2D) =
    tuple(RHSurfaceDomain(connection(b), ref(s).value)...)

backend_map_division(b::RH, f::Function, s::Shape2D, nu::Int, nv::Int) =
    let conn = connection(b)
        r = ref(s).value
        (u1, u2, v1, v2) = RHSurfaceDomain(conn, r)
        map_division(u1, u2, nu) do u
            map_division(v1, v2, nv) do v
                f(RHSurfaceFrameAt(conn, r, u, v))
            end
        end
    end

realize(b::RH, s::Sphere) =
    RHSphere(connection(b), s.center, s.radius)
realize(b::RH, s::Torus) =
    RHTorus(connection(b), s.center, vz(1, s.center.cs), s.re, s.ri)
realize(b::RH, s::Cuboid) =
    RHIrregularPyramidFrustum(connection(b), [s.b0, s.b1, s.b2, s.b3], [s.t0, s.t1, s.t2, s.t3])
realize(b::RH, s::RegularPyramidFrustum) =
    RHIrregularPyramidFrustum(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::RH, s::RegularPyramid) =
    RHIrregularPyramid(connection(b),
                       regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                       add_z(s.cb, s.h))
realize(b::RH, s::IrregularPyramid) =
    RHIrregularPyramid(connection(b), s.bs, s.t)
realize(b::RH, s::RegularPrism) =
    let cbs = regular_polygon_vertices(s.edges, s.cb, s.r, s.angle, s.inscribed)
        RHIrregularPyramidFrustum(connection(b),
                                  cbs,
                                  map(p -> add_z(p, s.h), cbs))
    end
realize(b::RH, s::IrregularPrism) =
    RHIrregularPyramidFrustum(connection(b),
                              s.bs,
                              map(p -> (p + s.v), s.bs))

realize(b::RH, s::RightCuboid) =
    RHXYCenteredBox(connection(b), s.cb, vx(1, s.cb.cs), vy(1, s.cb.cs), s.width, s.height, s.h)

realize(b::RH, s::Box) =
    RHBox(connection(b), s.c, vx(1,s.c.cs), vy(1,s.c.cs), s.dx, s.dy, s.dz)

realize(b::RH, s::Cone) =
  RHCone(connection(b), add_z(s.cb, s.h), s.r, s.cb)

realize(b::RH, s::ConeFrustum) =
  RHConeFrustum(connection(b), s.cb, s.rb, s.cb + vz(s.h, s.cb.cs), s.rt)

realize(b::RH, s::Cylinder) =
  RHCylinder(connection(b), s.cb, s.r, s.cb + vz(s.h, s.cb.cs))

#realize(b::RH, s::Circle) = RHCircle(connection(b),

backend_extrusion(b::RH, s::Shape, v::Vec) =
    and_mark_deleted(
        map_ref(s) do r
            RHExtrusion(connection(b), r, v)
        end,
        s)

backend_sweep(b::RH, path::Shape, profile::Shape, rotation::Real, scale::Real) =
  map_ref(profile) do profile_r
    map_ref(path) do path_r
      RHSweepPathProfile(connection(b), path_r, profile_r, rotation, scale)
    end
  end

#=
realize(b::RH, s::Revolve) =
  and_delete_shape(
    map_ref(s.profile) do r
      RHRevolve(connection(b), r, s.p, s.n, s.start_angle, s.amplitude)
    end,
    s.profile)

backend_loft_points(b::Backend, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  let f = (ruled ? (closed ? polygon : line) : (closed ? closed_spline : spline))
    and_delete_shapes(ref(f(map(point_position, profiles), backend=b)),
                      vcat(profiles, rails))
  end

backend_loft_curves(b::RH, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  and_delete_shapes(RHLoft(connection(b),
                             collect_ref(profiles),
                             collect_ref(rails),
                             ruled, closed),
                    vcat(profiles, rails))

backend_loft_surfaces(b::RH, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  backend_loft_curves(b, profiles, rails, ruled, closed)

backend_loft_curve_point(b::RH, profile::Shape, point::Shape) =
  and_delete_shapes(RHLoft(connection(b),
                             vcat(collect_ref(profile), collect_ref(point)),
                             [],
                             true, false),
                    [profile, point])

backend_loft_surface_point(b::RH, profile::Shape, point::Shape) =
  backend_loft_curve_point(b, profile, point)
=#
unite_refs(b::RH, refs::Vector{<:RHRef}) =
    RHUnionRef(tuple(refs...))

unite_refs(b::RH, r::RHUnionRef) =
  r

unite_ref(b::RH, r0::RHNativeRef, r1::RHNativeRef) =
  (RHUnite(connection(b), r0.value, r1.value); r0)

intersect_ref(b::RH, r0::RHNativeRef, r1::RHNativeRef) =
    let refs = RHIntersect(connection(b), r0.value, r1.value)
        n = length(refs)
        if n == 0
            RHEmptyRef()
        elseif n == 1
            RHNativeRef(refs[1])
        else
            RHUnionRef(map(RHNativeRef, tuple(refs...)))
        end
    end

subtract_ref(b::RH, r0::RHNativeRef, r1::RHNativeRef) =
    let refs = RHSubtract(connection(b), r0.value, r1.value)
        n = length(refs)
        if n == 0
            RHEmptyRef()
        elseif n == 1
            if is_empty_guid(refs[1]) # failed
                RHSubtractionRef(r0, tuple(r1))
            else
                RHNativeRef(refs[1])
            end
        else
            RHUnionRef(map(RHNativeRef, tuple(refs...)))
        end
    end

subtract_ref(b::RH, r0::RHRef, r1::RHUnionRef) =
  foldl((r0,r1)->subtract_ref(b,r0,r1), r1.values, init=r0)

#=
slice_ref(b::RH, r::RHNativeRef, p::Loc, v::Vec) =
  (RHSlice(connection(b), r.value, p, v); r)

slice_ref(b::RH, r::RHUnionRef, p::Loc, v::Vec) =
  map(r->slice_ref(b, r, p, v), r.values)

=#

#=
realize(b::RH, s::IntersectionShape) =
  foldl((r0,r1)->intersect_ref(b,r0,r1), UniversalRef{RHId}(), map(ref, s.shapes))
  =#

realize(b::RH, s::Slice) =
  slice_ref(b, ref(s.shape), s.p, s.n)

realize(b::RH, s::Move) =
  let r = map_ref(s.shape) do r
            RHMove(connection(b), r, s.v)
            r
          end
    mark_deleted(s.shape)
    r
  end

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

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

realize(b::RH, s::Mirror) =
  and_delete_shape(map_ref(s.shape) do r
                    RHMirror(connection(b), r, s.p, s.n, false)
                   end,
                   s.shape)

realize(b::RH, s::UnionMirror) =
  let r0 = ref(s.shape),
      r1 = map_ref(r0) do r
            RHMirror(connection(b), r, s.p, s.n, true)
          end
    UnionRef((r0,r1))
  end

realize(b::RH, s::SurfaceGrid) =
  let (nu, nv) = size(s.points)
    RHSurfaceFromGrid(connection(b), nu, nv, reshape(s.points,:), s.closed_u, s.closed_v,
      max(2*floor(Int,nu/30)+1, 2), max(2*floor(Int, nv/30)+1, 2))
  end
#=
realize(b::RH, s::Thicken) =
  and_delete_shape(
    map_ref(s.shape) do r
      RHThicken(connection(b), r, s.thickness)
    end,
    s.shape)

=#
# BIM

# Families

backend_get_family(b::RH, f::TableFamily) =
    RHCreateRectangularTableFamily(connection(b), f.length, f.width, f.height, f.top_thickness, f.leg_thickness)
backend_get_family(b::RH, f::ChairFamily) =
    RHCreateChairFamily(connection(b), f.length, f.width, f.height, f.seat_height, f.thickness)
backend_get_family(b::RH, f::TableChairFamily) =
    RHCreateRectangularTableAndChairsFamily(connection(b),
        ref(f.table_family), ref(f.chair_family),
        f.table_family.length, f.table_family.width,
        f.chairs_top, f.chairs_bottom, f.chairs_right, f.chairs_left,
        f.spacing)

backend_rectangular_table(b::RH, c, angle, family) =
    RHTable(connection(b), c, angle, ref(family))

backend_chair(b::RH, c, angle, family) =
    RHChair(connection(b), c, angle, ref(family))

backend_rectangular_table_and_chairs(b::RH, c, angle, family) =
    RHTableAndChairs(connection(b), c, angle, ref(family))

backend_slab(b::RH, profile, holes, thickness, family) =
    map_ref(b,
            r->RHExtrusion(connection(b), r, vz(thickness)),
            ensure_ref(b, backend_fill(b, profile)))

#Beams are aligned along the top axis.
realize(b::RH, s::Beam) =
  let profile = s.family.profile
      profile_u0 = profile.corner
      c = add_xy(s.cb, profile_u0.x + profile.dx/2, profile_u0.y + profile.dy/2)
      # need to test whether it is rotation on center or on axis
      o = add_y(loc_from_o_phi(s.cb, s.angle), -profile.dy/2)
    RHXYCenteredBox(connection(b), o, vx(1, o.cs), vy(1, o.cs), profile.dx, profile.dy, s.h)
  end

#Columns are aligned along the center axis.
realize(b::RH, s::FreeColumn) =
  let profile = s.family.profile
      profile_u0 = profile.corner
      c = add_xy(s.cb, profile_u0.x + profile.dx/2, profile_u0.y + profile.dy/2)
      # need to test whether it is rotation on center or on axis
      o = loc_from_o_phi(c, s.angle)
    RHXYCenteredBox(connection(b), o, vx(1, o.cs), vy(1, o.cs), profile.dx, profile.dy, s.h)
  end

realize(b::RH, s::Column) =
  let profile = s.family.profile,
      profile_u0 = profile.corner,
      c = add_xy(s.cb, profile_u0.x + profile.dx/2, profile_u0.y + profile.dy/2),
      base_height = s.bottom_level.height,
      height = s.top_level.height - base_height,
      # need to test whether it is rotation on center or on axis
      o = loc_from_o_phi(s.cb + vz(base_height), s.angle)
    RHXYCenteredBox(connection(b), o, vx(1, o.cs), vy(1, o.cs), profile.dx, profile.dy, height)
  end

backend_wall(b::RH, path, height, thickness, family) =
    RHPathWall(connection(b), backend_stroke(b, path), thickness, height)

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

backend_bounding_box(b::RH, shapes::Shapes) =
  RHBoundingBox(connection(b), collect_ref(shapes))

import Base.view
set_view(camera::XYZ, target::XYZ, lens::Real, b::RH) =
  RHView(connection(b), camera, target, lens)

get_view(b::RH) =
  let c = connection(b)
    RHViewCamera(c), RHViewTarget(c), RHViewLens(c)
  end

zoom_extents(b::RH) =
  RHZoomExtents(connection(b))

view_top(b::RH) =
    RHViewTop(connection(b))

rhino"public int DeleteAll()"
rhino"public int DeleteAllInLayer(String name)"
rhino"public void Delete(Guid id)"
rhino"public void DeleteMany(Guid[] ids)"

backend_delete_shapes(b::RH, shapes::Shapes) =
  RHDeleteMany(connection(b), collect_ref(shapes))

delete_all_shapes(b::RH) =
  RHDeleteAll(connection(b))

# Layers
rhino"public String CreateLayer(string name)"
rhino"public void SetLayerColor(ObjectId id, byte r, byte g, byte b)"
rhino"public void SetShapeColor(ObjectId id, byte r, byte g, byte b)"
rhino"public String CurrentLayer()"
rhino"public void SetCurrentLayer(String id)"
rhino"public String ShapeLayer(RhinoObject objId)"
rhino"public void SetShapeLayer(RhinoObject objId, String layerId)"

RHLayer = String

current_layer(b::RH)::RHLayer =
  RHCurrentLayer(connection(b))

current_layer(layer::RHLayer, b::RH) =
  RHSetCurrentLayer(connection(b), layer)

create_layer(name::String, b::RH) =
  RHCreateLayer(connection(b), name)

create_layer(name::String, color::RGB, b::RH) =
  let layer = RHCreateLayer(connection(b), name)
    RHSetLayerColor(connection(b), layer, color.r, color.g, color.b)
    layer
  end

delete_all_shapes_in_layer(layer::RHLayer, b::RH) =
  RHDeleteAllInLayer(connection(b), layer)

rhino"public Guid Clone(RhinoObject e)"

shape_from_ref(r, b::RH) =
  let c = connection(b)
    let code = RHShapeCode(c, r),
        ref = LazyRef(b, RHNativeRef(r))
              #ref = LazyRef(b, RHNativeRef(RHClone(connection(b), r))) # HACK CLONING!!!!!
      if code == 1
        point(RHPointPosition(c, r),
              backend=b, ref=ref)
      elseif code == 2
        circle(maybe_loc_from_o_vz(RHCircleCenter(c, r), RHCircleNormal(c, r)),
               RHCircleRadius(c, r),
               backend=b, ref=ref)
      elseif 3 <= code <= 6
        line(RHLineVertices(c, r),
             backend=b, ref=ref)
      elseif code == 7
        spline([xy(0,0)], false, false, #HACK obtain interpolation points
               backend=b, ref=ref)
      elseif 103 <= code <= 106
        polygon(RHLineVertices(c, r),
                backend=b, ref=ref)
      elseif code == 40
        # FIXME: frontier is missing
        surface(frontier=[], backend=b, ref=ref)
      elseif code == 41
        surface(frontier=[], backend=b, ref=ref)
      elseif code == 81
        sphere(backend=b, ref=ref)
      elseif code == 82
        cylinder(backend=b, ref=ref)
      elseif code == 83
        cone(backend=b, ref=ref)
      elseif code == 84
        torus(backend=b, ref=ref)
      else
        error("Unknown shape")
      end
    end
  end

"""            obj = _geometry_from_id(r)
            ref = native_ref(r)
            if isinstance(obj, geo.Point):
                return point.new_ref(ref, fromPt(obj.Location))
            elif isinstance(obj, geo.Curve):
                if rh.IsLine(r) or rh.IsPolyline(r):
                    if rh.IsCurveClosed(r):
                        return polygon.new_ref(ref, [fromPt(p) for p in rh.CurvePoints(r)[:-1]])
                    else:
                        return line.new_ref(ref, [fromPt(p) for p in rh.CurvePoints(r)])
                elif obj.IsCircle(Rhino.RhinoMath.ZeroTolerance):
                    return circle.new_ref(ref, fromPt(rh.CircleCenterPoint(r)), rh.CircleRadius(r))
                elif rh.IsCurveClosed(r):
                    return closed_spline.new_ref(ref, [fromPt(p) for p in rh.CurvePoints(r)])
                else:
                    return spline.new_ref(ref, [fromPt(p) for p in rh.CurvePoints(r)])
            elif rh.IsObject(r) and rh.IsObjectSolid(r):
                return solid(native_ref(r))
            elif rh.IsSurface(r) or rh.IsPolysurface(r):
                return surface(native_ref(r))
            else:
                raise RuntimeError("{0}: Unknown Rhino object {1}".format('shape_from_ref', r))
                """

#
rhino"public Point3d[] GetPosition(string prompt)"

select_position(prompt::String, b::RH) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = RHGetPosition(connection(b), prompt)
      length(ans) > 0 ? ans[1] : nothing
    end
  end

select_positions(prompt::String, b::RH) =
  let ps = Loc[]
      p = nothing
    @info "$(prompt) on the $(b) backend."
    while ((p = RHGetPosition(connection(b), prompt)) |> length) > 0
        push!(ps, p...)
    end
    ps
  end

rhino"public Guid[] GetPoint(string prompt)"
rhino"public Guid[] GetPoints(string prompt)"

# HACK: The next operations should receive a set of shapes to avoid re-creating already existing shapes

select_point(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, RHGetPoint)

select_points(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, RHGetPoints)

rhino"public Guid[] GetCurve(string prompt)"
rhino"public Guid[] GetCurves(string prompt)"

select_curve(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, RHGetCurve)

select_curves(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, RHGetCurves)

rhino"public Guid[] GetSurface(string prompt)"
rhino"public Guid[] GetSurfaces(string prompt)"

select_surface(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, RHGetSurface)

select_surfaces(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, RHGetSurfaces)

rhino"public Guid[] GetSolid(string prompt)"
rhino"public Guid[] GetSolids(string prompt)"

select_solid(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, RHGetSolid)

select_solids(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, RHGetSolids)

rhino"public Guid[] GetShape(string prompt)"
rhino"public Guid[] GetShapes(string prompt)"

select_shape(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, RHGetShape)

select_shapes(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, RHGetShapes)

captured_shape(b::RH, handle) =
  shape_from_ref(handle, b)
#
captured_shapes(b::RH, handles) =
  map(handles) do handle
      shape_from_ref(handle, b)
  end

generate_captured_shape(s::Shape, b::RH) =
    println("captured_shape(rhino, $(ref(s).value))")

generate_captured_shapes(ss::Shapes, b::RH) =
  begin
    print("captured_shapes(autocad, [")
    for s in ss
      print(ref(s).value)
      print(", ")
    end
    println("])")
  end

rhino"public Guid[] GetAllShapes()"
rhino"public Guid[] GetAllShapesInLayer(String name)"

all_shapes(b::RH) =
  [shape_from_ref(r) for r in RHGetAllShapes(connection(b))]

all_shapes_in_layer(layer, b::RH) =
  [shape_from_ref(r) for r in RHGetAllShapesInLayer(connection(b), layer)]


text:
empty!(InternedStrings.pool)

"This function makes use of `xs` in a way no optimizer can possibly remove"
function use(xs...)
    mktemp() do fn, fh
        print(fh, xs)
    end
end

@testset "Basic String Functionality" begin let
    empty!(InternedStrings.pool)

    s = intern("Hello My Friends1")

    @test length(s) == length("Hello My Friends1")

    @test startswith(s, "Hello")

    @test s == s
    @test s == "Hello My Friends1"

    @test addr_eq(intern(s), s)
end end

@testset "Interning" begin
    let a = "Hello My Friends2",
        b = join(["Hello", "My", "Friends2"], " ")

        empty!(InternedStrings.pool)

        # sanity check that strings are not already Interning
        @test !addr_eq(a, b)

        ai = intern(a)
        bi = intern(b)
        @test addr_eq(ai, bi)

        @test intern("a $(2*54) c") == "a 108 c"
    end
end

@testset "ID check" begin
    let a = "Gold", b = String(b"Gold")

        empty!(InternedStrings.pool)

        target_addr = pointer(a)

        a = intern(a)
        @test pointer(a) == target_addr

        @test pointer(b) != target_addr

        b = intern(b)
        @test pointer(b) == target_addr

        @test pointer(intern(SubString("Gold", 1))) == target_addr

        use(a,b)
    end
end

using Base.GC
@testset "Garbage Collection 1" begin let
    empty!(InternedStrings.pool)
    @test length(InternedStrings.pool)==0
    ai =  intern("Hello My Friends3")
    ai = [44] #remove the reference
    GC.gc();
    @test 0<=length(InternedStrings.pool)<=1 #May or may not have been collected yet
end end

@testset "Garbage Collection 2" begin let
    empty!(InternedStrings.pool)
    @test length(InternedStrings.pool)==0
    ai = intern("Hello My Friends4")
    bi = intern(join(["Hello", "My", "Friends4"], " "))
    @test addr_eq(ai, bi)
    @test length(InternedStrings.pool)==1
    use(ai,bi)
    ai = [44]
    GC.gc()
    @test length(InternedStrings.pool)==1 #don't collect when only one reference is gone
    use(bi)
    bi=[32]
    GC.gc()
    @test 0<=length(InternedStrings.pool)<=1
end end

using Random
Random.seed!(1)
@testset "Garbage Collection stress test" begin let
    empty!(InternedStrings.pool)
    oldpoolsize = length(InternedStrings.pool)
    function checkpool(op)
        GC.gc()
        @test op(length(InternedStrings.pool), oldpoolsize)
        oldpoolsize = length(InternedStrings.pool)
    end

    originals = [Random.randstring(rand(1:1024)) for _ in 1:10^5]
    n_orginals = length(originals)

    interns = intern.(originals);
    checkpool(>)

    for ii in 1:10^5
        push!(interns, intern(rand(originals)))
    end
    checkpool(==)
    originals = nothing
    checkpool(==)

    for ii in 1:30
        Random.shuffle!(interns)
        for jj in 1:1000
            pop!(interns)
        end
        checkpool(<=)
    end

    # This one matters:
    @test length(InternedStrings.pool) < n_orginals
end end

ext: jl


text:
@expectation function radon_hierarchical(
    log_radon,
    county_idx,
    floor,
    num_counties;
    log_sigma_alpha=-2.1,
    log_sigma_beta=-1.5,
    cost_fn=radon_target_f
)
    mu_alpha ~ Normal(0, 10)
    sigma_alpha = exp(log_sigma_alpha)

    mu_beta ~ Normal(0, 10)
    sigma_beta = exp(log_sigma_beta)

    alpha ~ filldist(Normal(mu_alpha, sigma_alpha), num_counties)
    beta ~ filldist(Normal(mu_beta, sigma_beta), num_counties)

    log_eps ~ transformed(truncated(Cauchy(0, 5), 0, Inf))
    eps = exp(log_eps)

    radon_est = alpha[county_idx] .+ beta[county_idx] .* floor
    log_radon ~ MvNormal(radon_est, eps*LinearAlgebra.I)

    # NOTE: Here the county idx and floor are supposed to be fixed.
    return cost_fn(
        alpha, beta, POST_PRED_COUNTY_IDX, POST_PRED_FLOOR)
end

function radon_target_f(alpha, beta, county_idx, floor)
    radon_level(a, b) = exp(a + b * floor)
    # Note we do not have a negative sign to "flip" the logistic function
    step_fun(a, b) = 1 / (1 + exp(5 * (radon_level(a, b) - 4)))
    vals = step_fun.(alpha, beta)
    return max(prod(vals), 1e-200)
end

const RADON_COST_FNS = (
    radon_target_f = radon_target_f,
)

function preprocess_radon_df(radon_df, data_file)
    num_counties = length(unique(radon_df[:,:county]))
    if data_file != "radon.csv"
        # Have N<84 counties but county_code is still between 0 and 84.
        # Need to map all the old county codes to new ones between 1 and N
        old_county_idx = unique(radon_df[:,:county_code])
        oldix2newix = Dict(old_county_idx[i] => i for i in 1:num_counties)
        county_idx = map(ix -> oldix2newix[ix], radon_df[:,:county_code])
    else
        county_idx = radon_df[:,:county_code] .+ 1
    end
    log_radon = radon_df[:,:log_radon]
    floor = radon_df[:,:floor]
    return log_radon, county_idx, floor, num_counties
end


text:
# The complete functionality of LibSixel consists of three parts:
#
# 1. [libsixel_jll], the pre-built [libsixel] artifacts using [BinaryBuilder]
# 2. `CModule.jl`, the C wrapper generated by [Clang.jl] using codes in `$PROJECT_ROOT/gen/`
# 3. The julia high-level interfaces that hides the C details.
#
# # Links
#
# [libsixel_jll]: https://github.com/JuliaBinaryWrappers/libsixel_jll.jl
# [libsixel]: https://github.com/saitoha/libsixel
# [BinaryBuilder]: https://github.com/JuliaPackaging/BinaryBuilder.jl
# [Clang.jl]: https://github.com/JuliaInterop/Clang.jl

# The high-level libsixel APIs do not support custom IO types. Hence here we provide our Julia
# flavor high-level API by wrapping the low-level libsixel API `sixel_encode` and
# `sixel_decode`/`sixel_decode_raw`

module LibSixel

export LibSixelEncoder, LibSixelDecoder

using ..SixelInterface
import ..SixelInterface: canonical_sixel_eltype

# This file is auto-generated by Clang.jl and is not expected to be modified directly.
# However, there're still some methods are modified:
# - `sixel_output_new`: ccall argtype for `priv` is relaxed from `Ptr{Cvoid}` to `Any`
include("CModule.jl")
using .C
include("types.jl")

using ImageCore, OffsetArrays
include("encoder.jl")  # high-level encoder API
include("decoder.jl") # high-level decoder API

end # module


text:
using TrainView

using CSV, DataFrames
using Statistics, LinearAlgebra

file = "../data/rail_detection_results_20210127.csv";
df = CSV.read(file, DataFrame);

Lu1 = df." Left_Rail_Bottom_X";
Lv1 = df." Left_Rail_Bottom_Y";
Lu2 = df." Left_Rail_Top_X";
Lv2 = df." Left_Rail_Top_Y";

Ru1 = df." Right_Rail_Bottom_X";
Rv1 = df." Right_Rail_Bottom_Y";
Ru2 = df." Right_Rail_Top_X";
Rv2 = df." Right_Rail_Top_Y";

# max_v = max(maximum(Lv1),maximum(Rv1))
max_v = 720
vs = [Lv1, Lv2, Rv1, Rv2];
for v in vs
    v[1:end] = v .- (max_v / 2.0)
end

# max_u = max(maximum(Lu2),maximum(Ru2))
max_u = 1280
us = [Lu1, Lu2, Ru1, Ru2];
for v in us
    v[1:end] = v .- (max_u / 2.0)
end

Lu = LinRange(mean(Lu1),mean(Lu2),20);
Lv = LinRange(mean(Lv1),mean(Lv2),20);
Ru = LinRange(mean(Ru1),mean(Ru2),20);
Rv = LinRange(mean(Rv1),mean(Rv2),20);

left_uvs =  [[Lu[i],Lv[i]] for i in eachindex(Lu)];
right_uvs = [[Ru[i],Rv[i]] for i in eachindex(Ru)];

# Initial guess at reference camera

# Richard Shenton email:
# Camera height: 2.065m
#
# Camera offset from centre line: 0.6m
# Pitch: 19.16 degrees down
# Yaw: 7.28 degrees to the left
# Roll: Assumed to be zero

trackprop = TrackProperties(track_gauge = 1.435 + 0.065)

cameraposition_reference = [0.0,0.604,-2.165];
cameraposition_reference = [0.0,0.6,-2.165];
ψθφ_ref = [0.0,19.16,-7.28] .* (pi/180.0)
ψθφ_ref = [0.0,-19.1,-7.28] .* (pi/180.0)

# cameraposition_reference = [0.0,-0.604,2.065];
# ψθφ_ref = [0.0,0.0,0.0] .* (pi/180.0)
# ψθφ_ref = [0.0,11.0,-13.0] .* (pi/180.0)
# ψθφ_ref = [0.0,10.5,-16.0] .* (pi/180.0)

camera_reference = VideoCamera(cameraposition_reference;
    focal_length = 5.8e-3,
    pixelspermeter = 1 / 5.5e-6,
    ψθφ = ψθφ_ref
)

uLs = left_track_image_u(camera_reference, trackprop, Lv);
uRs = right_track_image_u(camera_reference, trackprop, Rv);

using Plots
gr(yflip = true)
plot(Lu,Lv, lab = "")
plot!(Ru,Rv, lab = "")

plot!(uLs,Lv,linestyle=:dash, lab = "")
plot!(uRs,Rv,linestyle=:dash, lab = "")

using Optim

function calibration_objective(v::Vector)
    cameraposition = [0.0,v[1],v[2]];
    ψθφ = [v[3],v[4],v[5]] .* (pi/180.0)

    camera = VideoCamera(cameraposition;
        focal_length = 5.8e-3,
        pixelspermeter = 1 / 5.5e-6,
        ψθφ = ψθφ
    )

    uLs = left_track_image_u(camera, trackprop, Lv);
    uRs = right_track_image_u(camera, trackprop, Rv);

    scale = 10.0;
    # Y Z ψ θ φ
    mins = [-1.0,1.0,-3.0,-1.0,-15.0]
    maxs = [1.0 ,2.6, 3.0,50.0,15.0]
    constraints = map(1:length(v)) do j
        (v[j] > maxs[j]) * (exp(scale*(v[j]-maxs[j])^2) - 1.0) + (v[j] < mins[j]) * (exp(scale*(v[j] - mins[j])^2) - 1.0)
    end

    # return sum(constraints) + sum(abs2.(uLs - Lu)) / sum(abs2.(Lu) + abs2.(Lv)) + sum(abs2.(uRs - Ru)) / sum(abs2.(Ru) + abs2.(Rv))
    return sum(abs2.(uLs - Lu)) / sum(abs2.(Lu) + abs2.(Lv)) + sum(abs2.(uRs - Ru)) / sum(abs2.(Ru) + abs2.(Rv))
end

v = [cameraposition_reference[2:3]; ψθφ_ref[1:3] .* (180.0 / pi)];
calibration_objective(v)

tol = 1e-5
res = optimize(calibration_objective, v, Optim.Options(x_tol=tol));
res.minimum
calibration_objective(v)

camerapos = [0.0,v[1],v[2]]
ψθφ = [v[3],v[4],v[5]] .* (pi/180.0)

camera_reference = VideoCamera(camerapos;
    focal_length = 5.8e-3,
    pixelspermeter = 1 / 5.5e-6,
    ψθφ = ψθφ
)

uLs = left_track_image_u(camera_reference, trackprop, Lv);
uRs = right_track_image_u(camera_reference, trackprop, Rv);

plot(Lu,Lv, lab = "")
plot!(Ru,Rv, lab = "")

plot!(uLs,Lv,linestyle=:dash, lab = "")
plot!(uRs,Rv,linestyle=:dash, lab = "")

camera_reference = VideoCamera(cameraposition_reference;
    focal_length = 5.8e-3,
    pixelspermeter = 1 / 5.5e-6,
    ψθφ = ψθφ_ref
)

distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop;
    # choose_distortions = [:Y,:f],
    choose_distortions = [:Y,:Z,:f,:θ,:φ],
    # choose_distortions = [:Y,:Z,:θ,:φ,:ψ,:α,:β],
    iterations = 3)

camera = VideoCamera(camera_reference,distortion)

uLs2 = left_track_image_u(camera, trackprop, Lv);
uRs2 = right_track_image_u(camera, trackprop, Rv);

plot(Lu,Lv, lab = "")
plot!(Ru,Rv, lab = "")

plot!(uLs2,Lv,linestyle=:dash, lab = "")
plot!(uRs2,Rv,linestyle=:dash, lab = "")


text:
using Distributions, Random, DelimitedFiles, LinearAlgebra

include("KPNormal.jl")
using .KP
########################################################
###
# Robustness to Non_Normality
##########
mutable struct CLTExp
	cs
	thetas
	vs   #always represents precisions of muhat
	dist

	data  #contains the raw runs, always size (n, S)
	vs_orig  #original precisions pre-scaling by S
	useConstantPrecision #ensures muhat has same precision for all S
end

function CLTExp(cs, thetas, vs_, dist, S, useConstantPrecision)
	vs = copy(vs_)
	if !useConstantPrecision
		vs[:] .*= S
	end
	CLTExp(cs, thetas, vs, dist, zeros(length(cs), S),
		copy(vs), useConstantPrecision)
end

#Simulates by drawing S obs according to dist, and scaling by sqrt S
#Assumes object pre-configured!
function sim!(o::CLTExp, muhat)
	S = size(o.data, 2)
	n = size(o.data, 1)
	shifts = mean(o.dist)
	scales = 1/ std(o.dist) ./ sqrt.(o.vs_orig)
	if o.useConstantPrecision
		scales[:] .*= sqrt(S)
	end

	for i = 1:S
		#first generate variables with centered, scaled
		o.data[:, i] .= (rand(o.dist, n) .- mean(o.dist)) .* scales .+ o.thetas
	end
	muhat[:] = vec(mean(o.data, dims=2))
end

#sets object up to be used for size S
function resize!(o::CLTExp, S)
	o.data = zeros(length(o.cs), S)
	o.vs = copy(o.vs_orig)
	if !o.useConstantPrecision
		o.vs[:] .*= S
	end
end

######################
####
#For increasing N, fixed n, many simulations compare
#		#SAA
#		#FullInfo
#		#EB MM
#		#EB MLE
#		#SURE MSE
#		#OR MSE
#		#Diract Stein
#		#Box Stein
#		#tauOR
#Compare on metrics
#		#Value wrt theta
#		#Time to compute
#		#optimal value of tau0
function test_CLTharness(f, numRuns, o, S_grid, Gamma_min, Gamma_max, Gamma_step; onlyOracle = false)
	n = length(o.cs)
	muhat = zeros(Float64, n)
	xs = zeros(Float64, n)
	lam_t = 0.

	#write a header
	writedlm(f,  ["Run" "N" "Method" "thetaVal" "time" "tau0"], ',')

	for iRun = 1:numRuns
		for S in S_grid
			#reset the object and simulate
			resize!(o, S)
			#o.data = zeros(n, S)
			sim!(o, muhat)

			#fullInfo val
			t =
			  @elapsed xs[:] = x(o.cs, o.thetas)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "FullInfo" thetaval t 0.], ',')

			#Oracle
			t =
			  @elapsed xs[:], vals, objs = best_x_tau(o.cs, muhat, o.vs, o.thetas)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "OR" thetaval t vals[argmax(objs)]], ',')

			#Box with the optimized rate, i.e. h_n = n^-1/6
			h = n^-.16666
			t =
			  @elapsed xs[:], vals, objs = x_stein_box(o.cs, muhat, o.vs, h, tau_step = .05)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "BoxStein" thetaval t vals[argmax(objs)]], ',')

			#Oracle
			t =
			  @elapsed xs[:], Gamma_grid, objs = KP.x_l2reg_CV(o.cs, muhat, o.vs, o.thetas,
															Gamma_min=Gamma_min, Gamma_max=Gamma_max, Gamma_step=Gamma_step)

			Gammahat = Gamma_grid[argmax(objs)]
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "OracleReg" thetaval t Gammahat], ',')

			#Our Stein Approach to Regularization
			t =
			  @elapsed xs[:], Gamma_grid, objs = KP.x_stein_reg(o.cs, muhat, o.vs,
											Gamma_min=Gamma_min, Gamma_max=Gamma_max, Gamma_step=Gamma_step)

			Gammahat = Gamma_grid[argmax(objs)]
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "SteinReg" thetaval t Gammahat], ',')

			if onlyOracle
				continue
			end

			#SAA
			t =
			  @elapsed xs[:] = KP.x(o.cs, muhat)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "SAA" thetaval t 0.], ',')

			#Tau MLE
			t =
			  @elapsed tauMLE, xs[:] = x_MLE(o.cs, muhat, o.vs)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "EB_MLE" thetaval t tauMLE], ',')

			#Tau MM
			t =
			  @elapsed tauMM, xs[:] = x_MM(o.cs, muhat, o.vs)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "EB_MM" thetaval t tauMM], ',')

			#Oracle MSE
			t =
			  @elapsed xs[:], tau_CV = x_OR_MSE(o.cs, muhat, o.thetas, o.vs)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "OR_MSE" thetaval t tau_CV], ',')

			#Sure MSE
			t =
			  @elapsed xs[:], tau_CV = x_sure_MSE(o.cs, muhat, o.vs)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "SURE_MSE" thetaval t tau_CV], ',')

			#Hold-Out validation
			if S >= 2
				t =
				  @elapsed xs[:], tau_grid, objs = KP.x_kFoldCV(o.cs, o.vs/S, o.data, 2; tau_grid=collect(0.:.05:5))

				thetaval = dot(o.thetas, xs)/n
				writedlm(f,  [iRun S "HO_EB" thetaval t tau_grid[argmax(objs)]], ',')

			end

			#5-Fold validation
			if S >= 5
				t =
				  @elapsed xs[:], tau_grid, objs = KP.x_kFoldCV(o.cs, o.vs/S, o.data, 5; tau_grid=collect(0.:.05:5))

				thetaval = dot(o.thetas, xs)/n
				writedlm(f,  [iRun S "K5_EB" thetaval t tau_grid[argmax(objs)]], ',')
			end

			#LOO validation
			if S >= 2
				t =
				  @elapsed xs[:], tau_grid, objs = KP.x_kFoldCV(o.cs, o.vs/S, o.data, S; tau_grid=collect(0.:.05:5))

				thetaval = dot(o.thetas, xs)/n
				writedlm(f,  [iRun S "LOO_EB" thetaval t tau_grid[argmax(objs)]], ',')
			end

			####
			# Regularization
			####
			#Robust values
			thresh = sqrt(2*log(1/.1))
			t =
			  @elapsed xs[:] = KP.x_robFW(o.cs, muhat, o.vs, thresh, TOL=1e-4)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "FWRO_Eps_.1" thetaval t thresh], ',')

			thresh = sqrt(2*log(1/.05))
			t =
			  @elapsed xs[:] = KP.x_robFW(o.cs, muhat, o.vs, thresh, TOL=1e-4)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "FWRO_Eps_.05" thetaval t thresh], ',')

			thresh = sqrt(2*log(1/.01))
			t =
			  @elapsed xs[:] = KP.x_robFW(o.cs, muhat, o.vs, thresh, TOL=1e-4)
			thetaval = dot(o.thetas, xs)/n
			writedlm(f,  [iRun S "FWRO_Eps_.01" thetaval t thresh], ',')

			#####
			## Add K Fold values IF S is adequately large
			#####
			#Hold Out Validation choice
			if S >= 2
				t =
				  @elapsed xs[:], Gamma_grid, objs = KP.x_l2reg_kFoldCV(o.cs, o.vs/S, o.data, 2;
	                    Gamma_step=Gamma_step, Gamma_min=Gamma_min, Gamma_max=Gamma_max)

				Gammahat = Gamma_grid[argmax(objs)]
				thetaval = dot(o.thetas, xs)/n
				writedlm(f,  [iRun S "HO_Reg" thetaval t Gammahat], ',')
			end

			# #5-Fold Validation
			if S >= 5
				t =
				  @elapsed xs[:], Gamma_grid, objs = KP.x_l2reg_kFoldCV(o.cs, o.vs/S, o.data, 5;
	                    Gamma_step=Gamma_step, Gamma_min=Gamma_min, Gamma_max=Gamma_max)

				Gammahat = Gamma_grid[argmax(objs)]
				thetaval = dot(o.thetas, xs)/n
				writedlm(f,  [iRun S "K5_Reg" thetaval t Gammahat], ',')
			end

			#LOO Validation
			if S >= 2
				t =
				  @elapsed xs[:], Gamma_grid, objs = KP.x_l2reg_kFoldCV(o.cs, o.vs/S, o.data, S;
		                Gamma_step=Gamma_step, Gamma_min=Gamma_min, Gamma_max=Gamma_max)
				Gammahat = Gamma_grid[argmax(objs)]
				thetaval = dot(o.thetas, xs)/n
				writedlm(f,  [iRun S "LOO_Reg" thetaval t Gammahat], ',')
			end
		end
		flush(f)
	end
end

#helper
function get_dist(dist_type)
	if dist_type == "uniform"
		dist = Uniform()
	elseif dist_type == "bernoulli"  #discrete, no tV
		dist = Bernoulli(.5)
	elseif dist_type == "beta"
		dist = Beta(2, 2)
	elseif dist_type == "beta5"  #skewed
		dist = Beta(2, 5)
	elseif dist_type == "beta2"  #skewed
		dist = Beta(5, 2)
	elseif dist_type == "arcsine"
		dist = Arcsine()
	elseif dist_type == "triangle"
		dist = TriangularDist(-1, 1)
	elseif dist_type == "exponential" #not subgaussian
		dist = Exponential()
	elseif dist_type == "t"  #no higher moments
		dist = TDist(3)
	elseif dist_type == "pareto" #no higher moments
		dist = Pareto(3)
	else
		throw("Distribution Type Not Recognized: $(dist_type)")
	end
	return dist
end

function test_POAPCLT(file_out, param_path, numRuns, n, S_grid, seed, dist_type, Gammamin, Gammamax, Gamma_step; onlyOracle=false)
	Random.seed!(seed)
	dat, header = readdlm(param_path, ',', header=true)

	@assert n <= size(dat, 1) "Param file too short for n"
	cs = vec(dat[1:n, 3])
	thetas = vec(dat[1:n, 1])
	vs = vec(dat[1:n, 2])

	dist = get_dist(dist_type)
	o = CLTExp(cs, thetas, vs, dist, S_grid[1], true) #true forces constant precision

	file_name = "$(file_out)_$(dist_type)_$(seed).csv"
	f = open(file_name, "w")
	test_CLTharness(f, numRuns, o, S_grid, Gammamin, Gammamax, Gamma_step, onlyOracle=onlyOracle)
	close(f)
	return file_name
end

function test_POAPLargeSample(file_out, param_path, numRuns, n, S_grid, seed, dist_type, Gammamin, Gammamax, Gamma_step)
	Random.seed!(seed)
	dat, header = readdlm(param_path, ',', header=true)

	@assert n <= size(dat, 1) "Param file too short for n"
	cs = vec(dat[1:n, 3])
	thetas = vec(dat[1:n, 1])
	vs = vec(dat[1:n, 2])

	dist = get_dist(dist_type)
	o = CLTExp(cs, thetas, vs, dist, S_grid[1], false) #false forces constant precision

	file_name = "$(file_out)_$(dist_type)_$(seed).csv"
	f = open(file_name, "w")
	test_CLTharness(f, numRuns, o, S_grid, Gammamin, Gammamax, Gamma_step)
	close(f)
	return file_name
end

########################################################
#########
#Small run for pre0compilation
#test_POAPCLT("Results/temp_POAPCLT", "Results/param_portExp_mtn2.csv", 2, 100, [2 3], 8675309, "exponential")

max_stars_repo_path: Normal Model/testCLTHarness.jl

max_issues_repo_path: Normal Model/testCLTHarness.jl

max_forks_repo_path: Normal Model/testCLTHarness.jl


text:
module Determinism
export randBool, resetSeed

"""
seed used for the deterministic stochastic arithmetic
should be changed from one sample to another
"""
DSAseed = rand(Int64)

"""
changes the seed in order to run the computation and get a new result
"""
function resetSeed()
  global DSAseed = rand(Int64)
end

"""
Function that returns a boolean pseudo-randomly
garantees *determinism*: the same inputs will always produce the same output
and *diffusion*: all bits in the output have a 50% probability of being flipped if any bit of an input change
"""
function randBool(args...)
  global DSAseed # put global variable in scope
  pseudoRandomNumber = hash((DSAseed, args...)) # hash used as a random number
  return isodd(pseudoRandomNumber) # turn hash into a bool
end

end # module

# NOTES:
# diffusion is an important property in cryptographic cyphers (see paper "random as easy as 1, 2, 3" for concept and fast PRNG that have this property)
# means that a change of 1 bit in the input guarantees that all bits in the output have a 50% proba of changing
# splittable RNG paper (hash functions used as RNG): http://www.thesalmons.org/john/random123/papers/random123sc11.pdf


text:
#=
Equal number of elements in each bin.
=#
function bin_equalheight(X::Array;binsize,mergelast::Bool=false)
	inds = sortperm(X)

	inds_bybin = Array{Int}[]

	a = 1
	b = binsize
	while b < length(inds)
		push!(inds_bybin,inds[a:b])
		a += binsize
		b += binsize
	end

	if mergelast == false
		push!(inds_bybin,inds[a:end])
	else
		inds_bybin[end] = vcat(inds_bybin[end],inds[a:end])
	end

	return inds_bybin
end


text:
using PyPlot
PyPlot.ioff()
push!(LOAD_PATH, "/home/berceanu/Development/topo-photon/anc/scripts")
using PyCall
@pyimport mpl_toolkits.axes_grid1.inset_locator as axloc
@pyimport matplotlib.gridspec as gspec
import BP

# system parameters
const N = 45
const q = 11
const κ = 0.02
const γ = 0.001
const ν = collect(linspace(-3.45,-2.47,981))

const r = 11 # points in MBZ

# Nk should be an odd multiple of q
const Nk = r*q # zero-padded system size = no of k points
const l = div(Nk-1,2)
const x = -l:l
const δk = 2π/Nk # resolution in mom space
const k = x * δk

# full plot range, both in x and y
const xm = collect(-div(N-1,2):div(N-1,2))
# zoom in
const edge = 10
const st = findin(xm, -edge)[1]
const en = findin(xm,  edge)[1]

# exact spectrum, first 29 eigenvalues
exstates = BP.ExactStates(29, :landau, N, 1/q, κ)
# half the level spacing
const hf = (exstates.νs[2] - exstates.νs[1])/2

#for plotting filter markers
βlan = [0,2,4,6]
βsym = [0,1,9,20]
βreal = [3,6,15,26]
#we filter state η
ηlan = βlan + 1
ηsym = βsym + 1
ηreal = βreal + 1
#at energy
sω0lan = [exstates.νs[state]::Float64 for state in ηlan]
sω0sym = [exstates.νs[state]::Float64 for state in ηsym]
sω0real= [exstates.νs[state]::Float64 for state in ηreal]
# energy boundaries for state with β=4
ω₁= exstates.νs[4] - 0.005
ω₂= exstates.νs[6] + 0.005

δpmp(n₀::Int,m₀::Int) = BP.δpmp(N; n0=n₀, m0=m₀)
gausspmp(n₀::Int,m₀::Int) = BP.gausspmp(N; σ=1., n0=n₀, m0=m₀)
homopmp() = BP.homopmp(N)
randpmp(s::Int) = BP.randpmp(N; seed=s) #1234

prm = (1/q,γ,κ);

spδl = BP.Spectrum(ν,δpmp(5,5), :landau, prm...)
spgaussl = BP.Spectrum(ν,gausspmp(5,5), :landau, prm...)
sphoml = BP.Spectrum(ν,homopmp(), :landau, prm...)

# calculating all real space wfs
ψ = Array(Float64, (N, N, length(βreal)))
for (i,ω) in enumerate(sω0real)
    ψ[:,:,i] = abs2(BP.getstate(spδl, ω))
end

spgausss = BP.Spectrum(ν,gausspmp(5,5), :symmetric, prm...)
sphoms = BP.Spectrum(ν,homopmp(), :symmetric, prm...)

# calculating all mom space wfs
ψL = Array(Float64, (Nk, Nk, length(βlan)))
ψS = Array(Float64, (Nk, Nk, length(βsym)))
for (i,ω) in enumerate(sω0lan)
    ψL[:,:,i] = abs2(BP.myfft2(BP.getstate(sphoml, ω), k, k))
end
for (i,ω) in enumerate(sω0sym)
    ψS[:,:,i] = abs2(BP.myfft2(BP.getstate(spgausss, ω), k, k))
end

## averaging over 100 random phase distributions ##
intvec = zeros(Float64, length(ν));
A = spzeros(Complex{Float64}, N^2,N^2);
for j=1:100
    P=randpmp(j)
    for (i,ω) in enumerate(ν)
        BP.buildham_landau!(A, N,1/q,κ,γ,ω)
        intvec[i] += sum(abs2(A\P))
    end
end
sprandl = intvec./100;

# matplotlib parameters
matplotlib["rcParams"][:update](Dict("axes.labelsize" => 22,
                                     "axes.titlesize" => 20,
                                     "font.size" => 18,
                                     "legend.fontsize" => 14,
                                     "axes.linewidth" => 1.5,
                                     "font.family" => "serif",
                                     "font.serif" => "Computer Modern Roman",
                                     "xtick.labelsize" => 20,
                                     "xtick.major.size" => 5.5,
                                     "xtick.major.width" => 1.5,
                                     "ytick.labelsize" => 20,
                                     "ytick.major.size" => 5.5,
                                     "ytick.major.width" => 1.5,
                                     "text.usetex" => true,
                                     "figure.autolayout" => true))

ics=3.0
el = 0.86/(3/ics + 5)
b1 = 0.97-el
t2 = b1 - el/ics
b2 = t2 - 2el
t3 = b2 - el/ics

fig = plt[:figure]()
gs1 = gspec.GridSpec(1, 1)
gs1[:update](top=0.97, bottom=b1, right=0.98, left=0.16)
ax1 = plt[:subplot](get(gs1, (0, 0)))
ax1[:plot](spδl.νs,spδl.intensity,"k")
for (i,ω) in enumerate(sω0real)
    ax1[:axvline](x = ω, color="orange", ls="dotted", lw=2.0)
    ax1[:text](ω + hf/4, 8e3, string(βreal[i]))
end
ax1[:set_ylim](0, 1e4)
ax1[:yaxis][:set_ticks]([0, 1e4])
ax1[:yaxis][:set_ticklabels]([L"$0$", L"$10^4$"])
ax1[:text](ν[1] + hf, 5e3, "(a)")
gs2 = gspec.GridSpec(2, 1)
gs2[:update](top=t2, bottom=b2, hspace=0.0, right=0.98, left=0.16)
ax2 = plt[:subplot](get(gs2, (0, 0)))
ax3 = plt[:subplot](get(gs2, (1, 0)))
ax2[:plot](spgaussl.νs,spgaussl.intensity,"k")
ax2[:set_ylim](0, 3.5e3)
ax2[:yaxis][:set_ticks]([0, 3.5e3])
ax2[:yaxis][:set_ticklabels]([L"$0$", L"$3.5\!\! \times\!\! 10^3$"])
ax2[:text](ν[1] + hf, 1.75e3, "(b)")
ax3[:plot](spgausss.νs,spgausss.intensity, color="green", ls="dashed",
   linewidth=1.5)
for (i,ω) in enumerate(sω0sym)
    ax3[:axvline](x = ω, color="orange", ls="dotted", lw=2.0)
    ax3[:text](ω + hf/4, 1.6e4, string(βsym[i]))
end
ax3[:set_ylim](0, 2e4)
ax3[:yaxis][:set_ticks]([0, 1e4])
ax3[:yaxis][:set_ticklabels]([L"$0$", L"$10^4$"])
gs3 = gspec.GridSpec(2, 1)
gs3[:update](top=t3, bottom=0.11, hspace=0.4, right=0.98, left=0.16)
ax4 = plt[:subplot](get(gs3, (0, 0)))
ax4[:plot](sphoml.νs,sphoml.intensity,"k")
ax4[:plot](sphoms.νs,sphoms.intensity, color="green", ls="dashed",
   linewidth=1.5)
# insert with zoom of peak β=4
axins = axloc.inset_axes(ax4,
                        width="30%", # width = 30% of parent_bbox
                        height="50%",
                        loc=9) # located at upper middle part
# plot same thing as in parent box
axins[:plot](sphoms.νs,sphoms.intensity, color="green", ls="dashed",
   linewidth=1.5)
# but set much narrower limits
axins[:set_xlim]([ω₁, ω₂])
axins[:xaxis][:set_ticks]([ω₁, ω₂])
axins[:xaxis][:set_ticklabels]([string(round(ω₁,2)),
   string(round(ω₂,2))], fontsize=8)
axins[:set_ylim]([1e3, 1e4])
axins[:yaxis][:set_ticks]([1e3, 1e4])
axins[:yaxis][:set_ticklabels]([L"$10^3$", L"$10^4$"], fontsize=8)
# draw vertical lines at position of every exact eigenstate
axins[:axvline](x = exstates.νs[5], color="orange", ls="dotted", lw=2.0)
axins[:text](exstates.νs[5] + hf/8, 8e3, string(4), fontsize=8)
# draw a bbox of the region of the inset axes in the parent axes and
# connecting lines between the bbox and the inset axes area
axloc.mark_inset(ax4, axins, loc1=2, loc2=4, ec="0.", fc="none")
for (i,ω) in enumerate(sω0lan)
    ax4[:axvline](x = ω, color="orange", ls="dotted", lw=2.0)
    ax4[:text](ω + hf/4, 2e7, string(βlan[i]))
end
ax4[:set_ylim](0, 2.5e7)
ax4[:yaxis][:set_ticks]([0, 2.5e7])
ax4[:yaxis][:set_ticklabels]([L"$0$", L"$2.5\!\! \times\!\! 10^7$"])
ax4[:text](ν[1] + hf, 1.25e7, "(c)")
ax5 = plt[:subplot](get(gs3, (1, 0)))
ax5[:plot](ν,sprandl,"k")
ax5[:set_xlabel](L"$\omega_0/J$")
ax5[:set_ylim](0, 1.2e6)
ax5[:yaxis][:set_ticks]([0, 1.2e6])
ax5[:yaxis][:set_ticklabels]([L"$0$", L"$1.2\!\! \times\!\! 10^6$"])
ax5[:text](ν[1] + hf, 6e5, "(d)")
for (i, ax) in enumerate([ax1,ax2,ax3,ax4,ax5])
    ax[:set_xlim](ν[1], ν[end])
    i != 5 && ax[:set_xticklabels]([])
end
# set common y label to all subplots
fig[:text](0.022, 0.5, L"$\sum_{m,n} |a_{m,n}|^2$ (arb. units)",
   ha="center", va="center", rotation="vertical")
fig[:savefig]("../../figures/selection.pdf", pad_inches=0.0, transparent=true, bbox_inches="tight")
plt[:close](fig)

# plot w.fs. in real space
fig, axes = plt[:subplots](1,length(βreal), figsize=(10, 5))
for (i,ax) in enumerate(axes)

    ax[:imshow](ψ[st:en,st:en,i], origin="upper",
                     ColorMap("gist_heat_r"), interpolation="none",
                     extent=[-edge, edge, -edge, edge])
    ax[:set_xlabel](L"$m$")
    if i == 1 # leftmost panel
        ax[:set_ylabel](L"$n$")
    else
        ax[:set_yticklabels]([])
    end
end
fig[:savefig]("../../figures/real.pdf", transparent=true, pad_inches=0.0, bbox_inches="tight")
plt[:close](fig)

# plot w.fs. in  mom space
fig, axes = plt[:subplots](2,length(βlan), figsize=(10, 5))
for i = 1:length(βlan) # loop over columns
    # top row
    ax = axes[1,i]
    ax[:imshow](ψL[:,:,i], origin="upper", ColorMap("gist_heat_r"),
                     interpolation="none",
                     extent=[-π, π, -π, π])
    ax[:set_xticklabels]([])
    ax[:set_xticks]([-π,0,π])
    ax[:set_yticks]([-π,0,π])
    if i == 1 # leftmost panel
        ax[:set_ylabel](L"$k_y$")
        ax[:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    else
        ax[:set_yticklabels]([])
    end
    # bottom row
    ax = axes[2,i]
    if i == 3 # third pannel
        ax[:imshow](ψS[:,:,i], origin="upper",
                         ColorMap("gist_heat_r"), interpolation="none",
                         extent=[-π, π, -π, π],
                         vmin=0, vmax=270000)
    else
        ax[:imshow](ψS[:,:,i], origin="upper",
                         ColorMap("gist_heat_r"), interpolation="none",
                         extent=[-π, π, -π, π])
    end
    ax[:set_xlabel](L"$k_x$")
    ax[:set_xticks]([-π,0,π])
    ax[:set_yticks]([-π,0,π])
    ax[:set_xticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    if i == 1 # leftmost panel
        ax[:set_ylabel](L"$k_y$")
        ax[:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    else
        ax[:set_yticklabels]([])
    end
end
fig[:savefig]("../../figures/momentum.pdf", transparent=true, pad_inches=0.0, bbox_inches="tight")
plt[:close](fig)

max_stars_repo_name: berceanu/topo-photon

max_issues_repo_name: berceanu/topo-photon

max_forks_repo_name: berceanu/topo-photon


text:
# Solution found.  To write solution and print statistic to STDOUT.
function processSolution!(gs::GameState)
    writeSolution!(gs.solutions, gs.cells)
    PrintArray(gs)
    PrintStats(gs)
end

function PrintArray(gs::GameState)
    for i in 1:9
        (i % 3 == 1) && println()
        for j in 1:9
            (j % 3 == 1) && print(' ')
            number = gs.cells[Cell(i, j)]
            ch = number ∈ [1:9] ? '0' + number : '-'
            print(ch)
        end
        println()
    end
end

function PrintStats(gs::GameState)
    @printf "\nLevel Counts:\n\n"

    S = 1
    while gs.levelCount[S] == 0
        S += 1
    end

    i = 1

    while S ≤ 81
        cell = gs.sequence[S]
        @printf "(%d, %d):%4d " cell.row cell.col gs.levelCount[S]
        if (i > 5)
            @printf "\n"
            i = 1
        else
            i += 1
        end
        S += 1
    end

    @printf "\n\nCOUNT = %d\n" getTotalCount(gs.levelCount)
end

max_stars_repo_path: src/process solution.jl

max_stars_repo_name: Teo-ShaoWei/SudokuSolverSG

max_issues_repo_path: src/process solution.jl

max_issues_repo_name: Teo-ShaoWei/SudokuSolverSG

max_forks_repo_path: src/process solution.jl

max_forks_repo_name: Teo-ShaoWei/SudokuSolverSG


text:
module Interfaces

include("ephemerides.jl")
include("coords.jl")

end


text:
# VectorOfDecisions #
# ========================== #
struct VectorDecisionConstraintBridge{T, S <: MOI.AbstractVectorSet} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.VectorOfVariables, S}
    fixing_constraints::Vector{FixingConstraint{T}}
    decisions::VectorOfDecisions
end

function MOIB.Constraint.bridge_constraint(::Type{VectorDecisionConstraintBridge{T,S}},
                                           model,
                                           f::VectorOfDecisions,
                                           set::S) where {T, S <: MOI.AbstractVectorSet}
    # Perform the bridge mapping manually
    mapped_variables = map(f.decisions) do decision
        bridged = MOIB.bridged_function(
            model,
            AffineDecisionFunction{T}(SingleDecision(decision)))
        bridged.decision_part.terms[1].variable_index
    end
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    MOI.VectorOfVariables(mapped_variables),
                                    set)
    # Save the constraint indices and the decision to allow modifications
    return VectorDecisionConstraintBridge{T,S}(constraint, fill(FixingConstraint{T}(0), MOI.dimension(set)), f)
end

function MOI.supports_constraint(::Type{<:VectorDecisionConstraintBridge{T}},
                                 ::Type{VectorOfDecisions},
                                 ::Type{<:MOI.AbstractVectorSet}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:VectorDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:VectorDecisionConstraintBridge{T,S}}) where {T,S}
    return [(MOI.VectorOfVariables, S), (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T})]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:VectorDecisionConstraintBridge{T}},
                              ::Type{VectorOfDecisions},
                              S::Type{<:MOI.AbstractVectorSet}) where T
    return VectorDecisionConstraintBridge{T,S}
end

MOI.get(b::VectorDecisionConstraintBridge{T,S}, ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where {T,S} = 1
MOI.get(b::VectorDecisionConstraintBridge{T}, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}) where T =
    count(fc -> fc.value != 0, b.fixing_constraints)
MOI.get(b::VectorDecisionConstraintBridge{T}, ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where {T,S} = [b.constraint]
MOI.get(b::VectorDecisionConstraintBridge{T}, ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}) where {T} =
    filter(fc -> fc.value != 0, b.fixing_constraints)

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::VectorDecisionConstraintBridge)
    return bridge.decisions
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.delete(model::MOI.ModelLike, bridge::VectorDecisionConstraintBridge)
    MOI.delete(model, bridge.constraint)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::VectorDecisionConstraintBridge{T,S}, change::S) where {T,S}
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint, change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorDecisionConstraintBridge{T,S}, change::DecisionStateChange) where {T,S}
    i = findfirst(d -> d.decision == change.decision, bridge.decisions.decisions)
    # Switch on state transition
    if change.new_state == NotTaken
        if bridge.fixing_constraints[i].value != 0
            # Remove the fixing constraint
            MOI.delete(model, bridge.fixing_constraints[i])
            bridge.fixing_constraints[i] = FixingConstraint{T}(0)
        end
    end
    if change.new_state == Taken
        if bridge.fixing_constraints[i].value != 0
            # Remove any existing fixing constraint
            MOI.delete(model, bridge.fixing_constraints[i])
        end
        # Perform the bridge mapping manually
        aff = MOIB.bridged_function(model, AffineDecisionFunction{T}(bridge.decisions.decisions[i]))
        f = MOI.ScalarAffineFunction{T}(MOI.SingleVariable(aff.decision_part.terms[1].variable_index))
        # Get the decision value
        set = MOI.EqualTo(aff.decision_part.constant)
        # Add a fixing constraint to ensure that fixed decision is feasible.
        bridge.fixing_constraints[i] = MOI.add_constraint(model, f, set)
    end
    return nothing
end


text:
"""
    Duration{PeriodInSeconds}(value)

Construct a `Duration` object which stores an amount of time in units
`PeriodInSeconds` seconds. For instance, `Duration{60}(3)` represents 3 minutes
(three intervals of 60 seconds long).

The common units `weeks`, `days`, `hours`, `minutes`, `seconds`, `milliseconds`,
`microseconds` and `nanonseconds` are exported (as well as their non-plural
forms, `week`, etc) to make it convenient to construct any length of time,
such as `10seconds` or `10 * second`.
"""
struct Duration{PeriodInSeconds, T}
    value::T
end
Duration{P}(v::T) where {P,T} = Duration{P,T}(v)
Duration{P}(d::Duration) where {P} = convert(Duration{P}, d)
Duration{P,T}(d::Duration) where {P,T} = convert(Duration{P,T}, d)

@pure period(::Type{Duration{P}}) where {P} = P

@inline Base.get(d::Duration) = d.value

#@inline seconds(d::Duration) = get(d) * period(typeof(d))

# Conversion and promotion designed to maximize chance of staying "rational" if
# the user requests that (though Float64 is "default").
Base.convert(::Type{Duration{P}}, d::Duration{P,T}) where {P,T} = Duration{P,T}(get(d))
Base.convert(::Type{Duration{P1}}, d::Duration{P2}) where {P1,P2} = Duration{P1}((get(d) * P2) / P1)

Base.convert(::Type{Duration{P,T}}, d::Duration{P,T}) where {P,T} = Duration{P,T}(get(d))
Base.convert(::Type{Duration{P,T}}, d::Duration{P,T2}) where {P,T,T2} = Duration{P,T}(convert(T, get(d)))
Base.convert(::Type{Duration{P1,T}}, d::Duration{P2}) where {P1,P2,T} = Duration{P1,T}(convert(T, (get(d) * P2) / P1))

@pure function Base.promote_type(::Type{Duration{P1,T1}}, ::Type{Duration{P2,T2}}) where {P1,P2,T1,T2}
    # Note: this relation is assymetric if P1 == P2 but !(P1 === P2) (e.g. 1 vs 1.0)
    # Using promote_rule() results in stack overflows in this case
    if P1 === P2
        Duration{P1, promote_type(T1, T2)}
    elseif P1 < P2
        Duration{P1, Base.promote_op(/, promote_type(promote_type(T1, T2), typeof(P2)), typeof(P1))}
    else
        Duration{P2, Base.promote_op(/, promote_type(promote_type(T1, T2), typeof(P1)), typeof(P2))}
    end
end

@pure function Base.promote_op(f::F, ::Type{Duration{P1,T1}}, ::Type{Duration{P2,T2}}) where {F,P1,P2,T1,T2}
    if P1 === P2
        Duration{P1, Base.promote_op(f, T1, T2)}
    elseif P1 < P2
        Duration{P1, Base.promote_op(/, Base.promote_op(f, promote_type(T1, T2), typeof(P2)), typeof(P1))}
    else
        Duration{P2, Base.promote_op(/, Base.promote_op(f, promote_type(T1, T2), typeof(P1)), typeof(P2))}
    end
end

#Base.promote_op{P,T1,T2,F}(::F, ::Type{Duration{P,T1}}, ::Type{Duration{P,T2}}) = Duration{P, Base.promote_op(F,T1,T2)}
#Base.promote_op{P1,P2,T1,T2,F}(::F, ::Type{Duration{P1,T1}}, ::Type{Duration{P2,T2}}) = Duration{P1, promote_type(Base.promote_op(F,T1,T2),Base.promote_op(/, typeof(P1), typeof(P2)))}

function Base.show(io::IO, d::Duration{P}) where P
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " period of ") : print(io, " periods of ")
    print(io, P)
    (isa(P, Integer) && P == 1) ? print(io, " seconds") : print(io, " seconds")
end

# Mathematical operations
@inline Base.:-(d::Duration{P,T}) where {P,T} = Duration{P}(-get(d))

@inline Base.:+(d1::Duration, d2::Duration) = ((d1p, d2p) = promote(d1, d2); d1p + d2p)
@inline Base.:+(d1::Duration{P,T}, d2::Duration{P,T}) where {P,T} = Duration{P,T}(get(d1) + get(d2))

@inline Base.:-(d1::Duration, d2::Duration) = ((d1p, d2p) = promote(d1, d2); d1p - d2p)
@inline Base.:-(d1::Duration{P,T}, d2::Duration{P,T}) where {P,T} = Duration{P,T}(get(d1) - get(d2))

@inline Base.:/(d1::Duration{P1}, d2::Duration{P2}) where {P1,P2} = (get(d1) * P1) / (get(d2) * P2)
@inline Base.://(d1::Duration{P1}, d2::Duration{P2}) where {P1,P2} = (get(d1) * P1) // (get(d2) * P2)
@inline Base.div(d1::Duration{P1}, d2::Duration{P2}) where {P1,P2} = div((get(d1) * P1), (get(d2) * P2)) # Euclidean division, pairs with rem()
@inline Base.rem(d1::Duration{P1}, d2::Duration{P2}) where {P1,P2} = Duration{P1}(rem((get(d1) * P1), (get(d2) * P2)) / P1) # Negative iff d1.value is negative
@inline Base.fld(d1::Duration{P1}, d2::Duration{P2}) where {P1,P2} = fld((get(d1) * P1), (get(d2) * P2)) # Flooring division, pairs with mod()
@inline Base.mod(d1::Duration{P1}, d2::Duration{P2}) where {P1,P2} = Duration{P1}(mod((get(d1) * P1), (get(d2) * P2)) / P1) # Always positive

# Scaling
@inline Base.:*(d::Duration{P}, s::Number) where {P} = ((dp, sp) = promote(get(d), s); Duration{P}(dp * s))
@inline Base.:*(s::Number, d::Duration{P}) where {P} = ((dp, sp) = promote(get(d), s); Duration{P}(s * dp))

@inline Base.:/(d::Duration{P}, s::Number) where {P} = ((dp, sp) = promote(get(d), s); Duration{P}(dp / s))
@inline Base.://(d::Duration{P}, s::Number) where {P} = ((dp, sp) = promote(get(d), s); Duration{P}(dp // s))

# Interface with Base.Dates. Autoconversion will result in Floats, not Rationals
Base.convert(::Type{Duration}, d::Base.Dates.Week) = Duration{604800//1,Int64}(d.value)
Base.convert(::Type{Duration}, d::Base.Dates.Day) = Duration{86400//1,Int64}(d.value)
Base.convert(::Type{Duration}, d::Base.Dates.Hour) = Duration{3600//1,Int64}(d.value)
Base.convert(::Type{Duration}, d::Base.Dates.Minute) = Duration{60//1,Int64}(d.value)
Base.convert(::Type{Duration}, d::Base.Dates.Second) = Duration{1//1,Int64}(d.value)
Base.convert(::Type{Duration}, d::Base.Dates.Millisecond) = Duration{1//1000,Int64}(d.value)

Base.convert(::Type{Duration{P}}, d::Base.Dates.Week) where {P} = Duration{P}(d.value * P * 604800)
Base.convert(::Type{Duration{P}}, d::Base.Dates.Day) where {P} = Duration{P}(d.value * P * 86400)
Base.convert(::Type{Duration{P}}, d::Base.Dates.Hour) where {P} = Duration{P}(d.value * P * 3600)
Base.convert(::Type{Duration{P}}, d::Base.Dates.Minute) where {P} = Duration{P}(d.value * P * 60)
Base.convert(::Type{Duration{P}}, d::Base.Dates.Second) where {P} = Duration{P}(d.value * P)
Base.convert(::Type{Duration{P}}, d::Base.Dates.Millisecond) where {P} = Duration{P}(d.value * P / 1000)

Base.convert(::Type{Duration{P,T}}, d::Base.Dates.Week) where {P,T} = Duration{P,T}(d.value * P * 604800)
Base.convert(::Type{Duration{P,T}}, d::Base.Dates.Day) where {P,T} = Duration{P,T}(d.value * P * 86400)
Base.convert(::Type{Duration{P,T}}, d::Base.Dates.Hour) where {P,T} = Duration{P,T}(d.value * P * 3600)
Base.convert(::Type{Duration{P,T}}, d::Base.Dates.Minute) where {P,T} = Duration{P,T}(d.value * P * 60)
Base.convert(::Type{Duration{P,T}}, d::Base.Dates.Second) where {P,T} = Duration{P,T}(d.value * P)
Base.convert(::Type{Duration{P,T}}, d::Base.Dates.Millisecond) where {P,T} = Duration{P,T}(d.value * P / 1000)

Base.convert(::Type{Base.Dates.Week}, d::Duration{P}) where {P} = Base.Dates.Week(div(Int64(d.value * P), 604800))
Base.convert(::Type{Base.Dates.Day}, d::Duration{P}) where {P} = Base.Dates.Day(div(Int64(d.value * P), 86400))
Base.convert(::Type{Base.Dates.Hour}, d::Duration{P}) where {P} = Base.Dates.Hour(div(Int64(d.value * P), 3600))
Base.convert(::Type{Base.Dates.Minute}, d::Duration{P}) where {P} = Base.Dates.Minute(div(Int64(d.value * P), 60))
Base.convert(::Type{Base.Dates.Second}, d::Duration{P}) where {P} = Base.Dates.Second(Int(d.value * P))
Base.convert(::Type{Base.Dates.Millisecond}, d::Duration{P}) where {P} = Base.Dates.Millisecond(Int(d.value * P) * 1000)

Base.promote_rule(::Type{Duration{P,T}}, ::Type{DatesPeriod}) where {P,T, DatesPeriod <: Base.Dates.Period} = Duration{P, T}
Base.promote_rule(::Type{Duration{P,T}}, ::Type{DatesPeriod}) where {P,T, DatesPeriod <: Base.Dates.Millisecond} = Duration{P, promote_op(/, T, Int)}

# Export some convenience timescales as efficient singletons, e.g. week or weeks
const week = Duration{604800//1, Unit}(unit)
const day = Duration{86400//1, Unit}(unit)
const hour = Duration{3600//1, Unit}(unit)
const minute = Duration{60//1, Unit}(unit)
const second = Duration{1//1, Unit}(unit)
const millisecond = Duration{1//1000, Unit}(unit)
const microsecond = Duration{1//1_000_000, Unit}(unit)
const nanosecond = Duration{1//1_000_000_000, Unit}(unit)

const weeks = Duration{604800//1, Unit}(unit)
const days = Duration{86400//1, Unit}(unit)
const hours = Duration{3600//1, Unit}(unit)
const minutes = Duration{60//1, Unit}(unit)
const seconds = Duration{1//1, Unit}(unit)
const milliseconds = Duration{1//1000, Unit}(unit)
const microseconds = Duration{1//1_000_000, Unit}(unit)
const nanoseconds = Duration{1//1_000_000_000, Unit}(unit)

# Some nice outputs

function Base.show(io::IO, d::Duration{604800//1})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " week") : print(io, " weeks")
end

function Base.show(io::IO, d::Duration{86400//1})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " day") : print(io, " days")
end

function Base.show(io::IO, d::Duration{3600//1})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " hour") : print(io, " hours")
end

function Base.show(io::IO, d::Duration{60//1})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " minute") : print(io, " minutes")
end

function Base.show(io::IO, d::Duration{1//1})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " second") : print(io, " seconds")
end

function Base.show(io::IO, d::Duration{1/1000})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " millisecond") : print(io, " milliseconds")
end

function Base.show(io::IO, d::Duration{1//1_000_000})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " microsecond") : print(io, " microseconds")
end

function Base.show(io::IO, d::Duration{1//1_000_000_000})
    print(io, d.value)
    (isa(d.value, Integer) && d.value == 1) ? print(io, " nanosecond") : print(io, " nanoseconds")
end


text:
## Node-based analysis of species distributions

# This example demonstrates how to do a node-based comparison of species
# distributions, as described in Borregaard et al., MEE 2014.

# We will work with the distributions and phylogenetic relationships for all
# species of the family `Tyrannidae` in the Americas. The species occurrences
# are defined on a regular grid with a cellsize of 1 lat/long degree. This is
# one of the datasets used in the Borregaard _et al._ (2014) paper.

### Load data and create objects

# Load relevant packages
using CSV, DataFrames, SpatialEcology, Plots, Phylo

# Load the datasets
phylocom = CSV.read("tyrann_phylocom.tsv", DataFrame)
coord = CSV.read("tyrann_coords.tsv", DataFrame)

# Exercise
# Use the various inspection functions on the DataFrames
# What is the format?
# How many unique species are there?

# Exercise
# The site columns are used to match the two DataFrames together, but for that they
# need to be strings. What type are they now?.
# We can transform a single number to a String with e.g. `string(3)`.
# How can you use this knowledge to change all elements of the column with site names
# in each of the dataframes to be Strings, so they can be matched?

# Create the assemblage
tyrants = Assemblage(phylocom, coord)

# Exercise
# Create a nice species richness plot with a good color scale

# Next, we'll read in the phylogenetic tree
tree = open(parsenewick, "tyrannid_tree.tre") # most file opening methods in Julia start with
                                              # the `open` function - the first argument is a
                                              # function to apply to the file contents
                                              # `parsenewick` is an inbuilt function to parse
                                              # newick files

# Exercise
# plot the tree, experimenting with the keywords `showtips` (true or false),
# `treetype` (:fan or :dendrogram) or `tipfont` (try `(5,)`, or add a name of a
# font or a color to the tuple), or try a relevant keyword from here
# https://docs.juliaplots.org/latest/generated/attributes_series/ (many will not work)

# Exercise
# What happens if you call `sort!` on the tree, and then plot it?

# Phylo has functions to create iterators over nodes (`nodes` in Phylo include both
# tips/leaves and internal nodes). Using `nodenamefilter` you can get a lazy
# iterator over the internal nodes on the tree
# `isleaf` is an anonymous function used to filter
nodes = nodenamefilter(!isleaf, tree)

# Exercise
# We are interested in getting the 131'st node. How do you make the `nodes` iterator
# into a vector, and extract the 131'st element?

# Exercise
# Here's a more complex exercise
# we want to create a function that gives us a list of the species that
# descend from a given node
# use the functions `filter`, `isleaf`, and `getdescendants`
# Note that `isleaf` normally takes two arguments (try `?isleaf`), so
# should be called as `isleaf(tree, node)`

# Exercise
# Use your function to get the name of all species desciending from our 131st node

# Exercise
# Now use that species list to subset an `Assemblage` object that
# only contains the species descended from our node.

# Exercise
# Plot the species richness for that clade

# Exercise
# Take the steps you did above and merge them into a function that can be used
# to create the sub-Assemblage for any node on the tree

### Comparing the richness of sister clades
# The question we are interested in addressing here is: At a given node where the
# lineage splits into two sister clades - are the two descendant clades distributed
# differently? This could be an indication that an evolutionary or biogeographic
# event happened at that time, of consequence for the current distribution of
# the species.

# Exercise
# Use the `getchildren` function to get the two daughter nodes from our focal node
# Plot the species richness of the parent clade, and those of it's two descendents
# next to each other.

# Comment on the comparison - did something interesting happen here? What?

max_stars_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_issues_repo_name: EcoJulia/Julia-PhD-course-Copenhagen

max_forks_repo_name: EcoJulia/Julia-PhD-course-Copenhagen


text:
@with_kw mutable struct AdadeltaParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    ρ::T = 0.95
    ϵ::T = 1e-6
end

struct Adadelta <: AbstractSmoothing
    params::AdadeltaParameters{Float64}
    grms::Vector{Float64}
    xrms::Vector{Float64}
    xprev::Vector{Float64}

    function (::Type{Adadelta})(; kwargs...)
        return new(AdadeltaParameters(; kwargs...),Vector{Float64}(),Vector{Float64}(),Vector{Float64}())
    end
end

function initialize!(adadelta::Adadelta,x₀::Vector{Float64})
    resize!(adadelta.grms,length(x₀))
    resize!(adadelta.xrms,length(x₀))
    resize!(adadelta.xprev,length(x₀))
    adadelta.grms .= zero(length(x₀))
    adadelta.xrms .= zero(length(x₀))
    adadelta.xprev .= x₀
end

function smooth!(adadelta::Adadelta,klocal::Integer,kglobal::Integer,x::AbstractVector,gprev::AbstractVector,gcurr::AbstractVector)
    @unpack ρ,ϵ = adadelta.params
    Δx = x - adadelta.xprev
    adadelta.grms .= ρ*adadelta.grms + (1-ρ)*(gprev .* gprev)
    adadelta.xrms .= ρ*adadelta.xrms + (1-ρ)*(Δx .* Δx)
    gcurr .= gprev .* ((.√adadelta.xrms .+ ϵ) ./ (.√adadelta.grms .+ ϵ))
    adadelta.xprev .= x
end


text:
module PandorasBoxIntroCar

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/introCar" IntroCar(x::Integer, y::Integer, acceleration::Number=256.0, deceleration::Number=384.0, maxSpeed::Number=384.0, brokenDoor::Bool=false, keepCarSpeedOnExit::Bool=true, nitroAcceleration::Number=448.0, nitroMaxDuration::Number=3.0, nitroRegenMultiplier::Number=0.2, facing::Integer=1)

const placements = Ahorn.PlacementDict(
    "Intro Car (Pandora's Box)" => Ahorn.EntityPlacement(
        IntroCar
    )
)

Ahorn.editingOptions(entity::IntroCar) = Dict{String, Any}(
    "facing" => Dict{String, Any}(
        "Left" => -1,
        "Right" => 1
    )
)

const bodySprite = "scenery/car/body"
const wheelsSprite = "scenery/car/wheels"

function Ahorn.selection(entity::IntroCar)
    x, y = Ahorn.position(entity)

    return Ahorn.coverRectangles([
        Ahorn.getSpriteRectangle(bodySprite, x, y, jx=0.5, jy=1.0),
        Ahorn.getSpriteRectangle(wheelsSprite, x, y, jx=0.5, jy=1.0),
    ])
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::IntroCar, room::Maple.Room)
    x, y = Ahorn.position(entity)

    scaleX = get(entity.data, "facing", 1)

    Ahorn.drawSprite(ctx, wheelsSprite, x, y, jx=0.5, jy=1.0, sx=scaleX, sy=1)
    Ahorn.drawSprite(ctx, bodySprite, x, y, jx=0.5, jy=1.0, sx=scaleX, sy=1)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
# indexnotation/indexedobject.jl
#
# Defines a parameteric type to store indices at compile-time, as well as a
# hierarchy of types to assign these indices to objects and implement tensor
# operations with them.

import Base: +, -, *

immutable Indices{I}
end

abstract AbstractIndexedObject
indices(a::AbstractIndexedObject) = indices(typeof(a))

immutable IndexedObject{I,C,A,T} <: AbstractIndexedObject
    object::A
    α::T
    function IndexedObject(object::A, α::T)
        checkindices(object, I)
        new(object, α)
    end
end
Base.call{I,C, A, T}(::Type{IndexedObject{I,C}}, object::A, α::T=1) = IndexedObject{I,C, A, T}(object, α)

Base.conj{I}(a::IndexedObject{I,:N}) = IndexedObject{I,:C}(a.object, conj(a.α))
Base.conj{I}(a::IndexedObject{I,:C}) = IndexedObject{I,:N}(a.object, conj(a.α))

Base.eltype(A::IndexedObject) = promote_type(eltype(A.object), typeof(A.α))
Base.eltype{I,C, A, T}(::Type{IndexedObject{I,C, A, T}}) = promote_type(eltype(A), T)

*{I,C}(a::IndexedObject{I,C}, β::Number) = IndexedObject{I,C}(a.object, a.α*β)
*(β::Number, a::IndexedObject) = *(a, β)
-(a::IndexedObject) = *(a, -1)

@generated function indices{I,C,A,T}(::Type{IndexedObject{I,C,A,T}})
    J = tuple(unique2(I)...)
    meta = Expr(:meta, :inline)
    Expr(:block, meta, :($J))
end

indexify{I}(object, ::Indices{I}) = IndexedObject{I,:N}(object)

deindexify{I}(A::IndexedObject{I,:N}, ::Indices{I}) = A.α == 1 ? A.object : A.α*A.object
deindexify{I}(A::IndexedObject{I,:C}, ::Indices{I}) = A.α == 1 ? conj(A.object) : A.α*conj(A.object)

@generated function deindexify{I,C,J}(A::IndexedObject{I,C}, ::Indices{J}, T::Type = eltype(A))
    meta = Expr(:meta, :inline)
    indCinA, = trace_indices(I,J)
    conj = Val{C}
    quote
        $meta
        deindexify!(similar_from_indices(T, $indCinA, A.object, $conj), A, Indices{$J}())
    end
end

@generated function deindexify!{Idst, Isrc, C}(dst, src::IndexedObject{Isrc, C}, ::Indices{Idst}, β=0)
    meta = Expr(:meta, :inline)
    Jdst = unique2(Idst)
    length(Jdst) == length(Idst) || throw(IndexError("left-hand side cannot have partial trace: $Idst"))
    if length(Isrc) == length(Jdst)
        indCinA = add_indices(Isrc, Idst)
        :($meta;add!(src.α, src.object, Val{C}, β, dst, $indCinA))
    else
        indCinA, cindA1, cindA2 = trace_indices(Isrc, Idst)
        return :($meta;trace!(src.α, src.object, Val{C}, β, dst, $indCinA, $cindA1, $cindA2))
    end
end


text:
"""
    AbstractElements

Abstract base type for two-body element sets.

## Two Body Element Interface

- [`semi_major_axis`](@ref)
- [`eccentricity`](@ref)
- [`right_ascension`](@ref)
- [`argument_of_periapsis`](@ref)
- [`true_anomaly`](@ref)
"""
abstract type AbstractElements end

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# Interface definition
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #

"""
    semi_major_axis(::AbstractElements)

Retrieve the semi major axis orbital element.
"""
function semi_major_axis(::AbstractElements) end

"""
    eccentricity(::AbstractElements)

Retrieve the eccentricity orbital element.
"""
function eccentricity(::AbstractElements) end

"""
    right_ascension(::AbstractElements)

Retrieve the right ascension orbital element.
"""
function right_ascension(::AbstractElements) end

"""
    argument_of_periapsis(::AbstractElements)

Retrieve the argument of periapsis orbital element.
"""
function argument_of_periapsis(::AbstractElements) end

"""
    true_anomaly(::AbstractElements)

Retrieve the true anomaly orbital element.
"""
function true_anomaly(::AbstractElements) end

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# Default derived values
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# Note, when implementing a new orbital element set, it's a good idea to see if there are
# faster ways to calculate any of these. For example, in the modified equinoctial
# elements, the parameter (semi-latus rectum) is an orbital element. If a new `parameter`
# method is not defined, it will be calculated as
#
#        |-- semi_major_axis(m) -> a -|
#   m ---|                            |-- parameter(a, e)
#        |-- eccentricity(m) -> e ----|
#
# which is relatively expensive to a single field retrieval.

# Un-exported helper functions so that we don't have to write them everywhere
const _sma = semi_major_axis
const _ecc = eccentricity
const _raan = right_ascension
const _aop = argument_of_periapsis
const _ta = true_anomaly

# Semi-Latus Rectum
parameter(k::AbstractElements) = parameter(_sma(k), _ecc(k))

# Radii
periapsis_radius(k::AbstractElements) = _sma(k) * (1 - _ecc(k))
apoapsis_radius(k::AbstractElements) = _sma(k) * (1 + _ecc(k))
radius(k::AbstractElements) = radius(_sma(k), _ecc(k), _ta(k))

# True Longitude: See https://en.wikipedia.org/wiki/True_longitude
true_longitude(k::AbstractElements) = true_longitude(_raan(k), _aop(k), _ta(k))

# Flight Path Angle: angle between the local horizon and the velocity vector
flight_path_angle(k::AbstractElements) = flight_path_angle(_sma(k), _ecc(k), _ta(k))

# Closure Properties of Orbit
# Note: We don't need to define `isopen` as that is defined in terms of `isclosed`
isclosed(k::AbstractElements) = _ecc(k) |> isclosed

# Orbit Type Querying
iscircular(k::AbstractElements) = _ecc(k) |> iscircular
iselliptical(k::AbstractElements) = _ecc(k) |> iselliptical
isparabolic(k::AbstractElements) = _ecc(k) |> isparabolic
ishyperbolic(k::AbstractElements) = _ecc(k) |> ishyperbolic

# Mean Motion and Period
mean_motion(cb::CelestialBody, e::AbstractElements) = mean_motion(cb, semi_major_axis(e))
period(cb::CelestialBody, e::AbstractElements) = period(cb, _sma(e))

# Angular Momentum
angular_momentum(gm, e::AbstractElements) = angular_momentum(gm, _sma(e), _ecc(e))
angular_momentum(cb::CelestialBody, e::AbstractElements) = angular_momentum(_GM(cb), e)

# Dynamic Quantitie
circular_velocity(gm, e::AbstractElements) = circular_velocity(gm, _sma(e))
mean_motion(gm, e::AbstractElements) = mean_motion(gm, _sma(e))