text:
export rhino

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_api = @remote_functions :CS """
public void RunScript(string script)
public void SetView(Point3d position, Point3d target, double lens, bool perspective, string style)
public void View(Point3d position, Point3d target, double lens)
public void ViewTop()
public Point3d ViewCamera()
public Point3d ViewTarget()
public double ViewLens()
public Guid Point(Point3d p)
public Point3d PointPosition(Guid ent)
public Guid PolyLine(Point3d[] pts)
public Point3d[] LineVertices(RhinoObject id)
public Guid Spline(Point3d[] pts)
public Guid ClosedPolyLine(Point3d[] pts)
public Guid ClosedSpline(Point3d[] pts)
public Guid Circle(Point3d c, Vector3d n, double r)
public Point3d CircleCenter(RhinoObject obj)
public Vector3d CircleNormal(RhinoObject obj)
public double CircleRadius(RhinoObject obj)
public Guid Ellipse(Point3d c, Vector3d n, double radiusX, double radiusY)
public Guid Arc(Point3d c, Vector3d n, double radius, double startAngle, double endAngle)
public Guid JoinCurves(Guid[] objs)
public Guid SurfaceCircle(Point3d c, Vector3d n, double r)
public Guid SurfaceEllipse(Point3d c, Vector3d n, double radiusX, double radiusY)
public Guid SurfaceArc(Point3d c, Vector3d n, double radius, double startAngle, double endAngle)
public Guid SurfaceClosedPolyLine(Point3d[] pts)
public Guid SurfaceFrom(Guid[] objs)
public Guid Sphere(Point3d c, double r)
public Guid Torus(Point3d c, Vector3d vz, double majorRadius, double minorRadius)
public Brep Cylinder(Point3d bottom, double radius, Point3d top)
public Brep Cone(Point3d bottom, double radius, Point3d top)
public Brep ConeFrustum(Point3d bottom, double bottom_radius, Point3d top, double top_radius)
public Brep Box(Point3d corner, Vector3d vx, Vector3d vy, double dx, double dy, double dz)
public Brep XYCenteredBox(Point3d corner, Vector3d vx, Vector3d vy, double dx, double dy, double dz)
public Brep IrregularPyramid(Point3d[] pts, Point3d apex)
public Brep IrregularPyramidFrustum(Point3d[] bpts, Point3d[] tpts)
public Guid SurfaceFromGrid(int nU, int nV, Point3d[] pts, bool closedU, bool closedV, int degreeU, int degreeV)
public Brep[] Thicken(RhinoObject obj, double thickness)
public double[] CurveDomain(RhinoObject obj)
public double CurveLength(RhinoObject obj)
public Plane CurveFrameAt(RhinoObject obj, double t)
public Plane CurveFrameAtLength(RhinoObject obj, double l)
public Brep Extrusion(RhinoObject obj, Vector3d dir)
public Brep[] SweepPathProfile(RhinoObject path, RhinoObject profile, double rotation, double scale)
public Guid[] Intersect(RhinoObject obj0, RhinoObject obj1)
public Guid[] Subtract(RhinoObject obj0, RhinoObject obj1)
public Guid Move(Guid id, Vector3d v)
public Guid Scale(Guid id, Point3d p, double s)
public Guid Rotate(Guid id, Point3d p, Vector3d n, double a)
public Guid Mirror(Guid id, Point3d p, Vector3d n, bool copy)
public void Render(int width, int height, string path)
public void SaveView(int width, int height, string path)
public bool IsPoint(RhinoObject e)
public bool IsCircle(RhinoObject e)
public bool IsPolyLine(RhinoObject e)
public bool IsSpline(RhinoObject e)
public bool IsInterpSpline(RhinoObject e)
public bool IsClosedPolyLine(RhinoObject e)
public bool IsClosedSpline(RhinoObject e)
public bool IsInterpClosedSpline(RhinoObject e)
public bool IsEllipse(RhinoObject e)
public bool IsArc(RhinoObject e)
public bool IsText(RhinoObject e)
public byte ShapeCode(RhinoObject id)
public Guid ClosedPathCurveArray(Point3d[] pts, double[] angles)
public Brep[] PathWall(RhinoObject obj, double lThickness, double rThickness, double height)
public Brep RectangularTable(Point3d c, double angle, double length, double width, double height, double top_thickness, double leg_thickness)
public Brep[] BaseRectangularTable(double length, double width, double height, double top_thickness, double leg_thickness)
public int CreateRectangularTableFamily(double length, double width, double height, double top_thickness, double leg_thickness)
public Guid Table(Point3d c, double angle, int family)
public Brep[] BaseChair(double length, double width, double height, double seat_height, double thickness)
public int CreateChairFamily(double length, double width, double height, double seat_height, double thickness)
public Guid Chair(Point3d c, double angle, int family)
public GeometryBase[] BaseRectangularTableAndChairs(int tableFamily, int chairFamily, double tableLength, double tableWidth, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)
public int CreateRectangularTableAndChairsFamily(int tableFamily, int chairFamily, double tableLength, double tableWidth, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)
public Guid TableAndChairs(Point3d c, double angle, int family)
public String CreateLayer(String name, bool active, byte r, byte g, byte b)
public void SetLayerColor(ObjectId id, byte r, byte g, byte b)
public void SetShapeColor(ObjectId id, byte r, byte g, byte b)
public String CurrentLayer()
public void SetCurrentLayer(String id)
public void SetLayerActive(String name, bool active)
public String ShapeLayer(RhinoObject objId)
public void SetShapeLayer(RhinoObject objId, String layerId)
public Point3d[] CurvePointsAt(RhinoObject obj, double[] ts)
public Vector3d[] CurveTangentsAt(RhinoObject obj, double[] ts)
public Vector3d[] CurveNormalsAt(RhinoObject obj, double[] ts)
public int DeleteAll()
public int DeleteAllInLayer(String name)
public void Delete(Guid id)
public void DeleteMany(Guid[] ids)
public Guid Clone(RhinoObject e)
public Point3d[] GetPosition(string prompt)
public Guid[] GetPoint(string prompt)
public Guid[] GetPoints(string prompt)
public Guid[] GetCurve(string prompt)
public Guid[] GetCurves(string prompt)
public Guid[] GetSurface(string prompt)
public Guid[] GetSurfaces(string prompt)
public Guid[] GetSolid(string prompt)
public Guid[] GetSolids(string prompt)
public Guid[] GetShape(string prompt)
public Guid[] GetShapes(string prompt)
public Guid[] GetAllShapes()
public Guid[] GetAllShapesInLayer(String name)
"""
#=
public byte Sync()
public byte Disconnect()
public Entity InterpSpline(Point3d[] pts, Vector3d tan0, Vector3d tan1)
public Entity InterpClosedSpline(Point3d[] pts)
public Entity Text(string str, Point3d corner, Vector3d vx, Vector3d vy, double height)
public Guid SurfaceFromCurve(Entity curve)
public Guid IrregularPyramidMesh(Point3d[] pts, Point3d apex)
public Entity SolidFromGrid(int m, int n, Point3d[] pts, bool closedM, bool closedN, int level, double thickness)
public Guid NurbSurfaceFrom(ObjectId id)
public double[] SurfaceDomain(Entity ent)
public Frame3d SurfaceFrameAt(Entity ent, double u, double v)
public Guid Sweep(ObjectId pathId, ObjectId profileId, double rotation, double scale)
public Guid Loft(ObjectId[] profilesIds, ObjectId[] guidesIds, bool ruled, bool closed)
public void Unite(ObjectId objId0, ObjectId objId1)
public Guid Revolve(ObjectId profileId, Point3d p, Vector3d n, double startAngle, double amplitude)
public Point3d[] GetPoint(string prompt)
public Point3d[] BoundingBox(ObjectId[] ids)
public void ZoomExtents()
public void SetSystemVariableInt(string name, int value)
public int Command(string cmd)
public void DisableUpdate()
public void EnableUpdate()
public BIMLevel FindOrCreateLevelAtElevation(double elevation)
public BIMLevel UpperLevel(BIMLevel currentLevel, double addedElevation)
public double GetLevelElevation(BIMLevel level)
public Entity MeshFromGrid(int m, int n, Point3d[] pts, bool closedM, bool closedN)

public FloorFamily FloorFamilyInstance(double totalThickness, double coatingThickness)
public Entity LightweightPolyLine(Point2d[] pts, double[] angles, double elevation)
public Entity SurfaceLightweightPolyLine(Point2d[] pts, double[] angles, double elevation)
public ObjectId CreatePathFloor(Point2d[] pts, double[] angles, BIMLevel level, FloorFamily family)
#public TableFamily FindOrCreateTableFamily(double length, double width, double height, double top_thickness, double leg_thickness)
#public ChairFamily FindOrCreateChairFamily(double length, double width, double height, double seat_height, double thickness)
#public TableChairFamily FindOrCreateTableChairFamily(TableFamily tableFamily, ChairFamily chairFamily, int chairsOnTop, int chairsOnBottom, int chairsOnRight, int chairsOnLeft, double spacing)
#public Brep RectangularTable(Point3d c, double angle, TableFamily family)
#public Brep Chair(Point3d c, double angle, ChairFamily family)
#public Guid InstanceChair(Point3d c, double angle, ChairFamily family)
#public Brep[] RowOfChairs(Point3d c, double angle, int n, double spacing, ChairFamily family)
#public Brep[] CenteredRowOfChairs(Point3d c, double angle, int n, double spacing, ChairFamily family)
#public Brep[] RectangularTableAndChairs(Point3d c, double angle, TableChairFamily f)
=#

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), rhino_api)

# 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) =
    @remote(b, Circle(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
        @remote(b, ClosedPolyLine([c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)]))
    end
backend_stroke(b::RH, path::ArcPath) =
  backend_stroke_arc(b, path.center, path.radius, path.start_angle, path.amplitude)

backend_stroke(b::RH, path::OpenPolygonalPath) =
  	@remote(b, PolyLine(path.vertices))
backend_stroke(b::RH, path::ClosedPolygonalPath) =
    @remote(b, ClosedPolyLine(path.vertices))

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

backend_stroke_unite(b::RH, refs) = @remote(b, JoinCurves(refs))

#=backend_fill(b::RH, path::ClosedPolygonalPath) =
    @remote(b, SurfaceClosedPolyLine(path.vertices))
    backend_fill(b::RH, path::RectangularPath) =
        let c = path.corner,
            dx = path.dx,
            dy = path.dy
            @remote(b, SurfaceClosedPolyLine([c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)]))
        end
backend_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) = @remote(b, SurfaceFrom(gs))
backend_fill_curves(b::RH, g::Guid) = @remote(b, SurfaceFrom([g]))

backend_stroke_line(b::RH, vs) = @remote(b, PolyLine(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
            @remote(b, Arc(center, vz(1, center.cs), radius, start_angle, end_angle))
        else
            @remote(b, Arc(center, vz(1, center.cs), radius, end_angle, start_angle))
        end
    end

realize(b::RH, s::Point) =
  @remote(b, Point(s.position))

realize(b::RH, s::Line) =
  @remote(b, PolyLine(s.vertices))

realize(b::RH, s::Spline) =
  if (s.v0 == false) && (s.v1 == false)
    @remote(b, Spline(s.points))
  elseif (s.v0 != false) && (s.v1 != false)
    @remote(b, InterpSpline(s.points, s.v0, s.v1))
  else
    @remote(b, InterpSpline(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) =
  @remote(b, InterpClosedSpline(s.points))

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

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

realize(b::RH, s::Ellipse) =
    @remote(b, Ellipse(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) =
  @remote(b, ClosedPolyLine(s.vertices))

realize(b::RH, s::RegularPolygon) =
  @remote(b, ClosedPolyLine(regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed)))

realize(b::RH, s::Rectangle) =
  @remote(b, ClosedPolyLine([s.corner, add_x(s.corner, s.dx), add_xy(s.corner, s.dx, s.dy), add_y(s.corner, s.dy)]))

backend_map_division(b::RH, f::Function, s::Shape1D, n::Int) =
  let conn = connection(b),
      r = ref(s).value,
      (t1, t2) = @remote(b, CurveDomain(r)),
      ti = division(t1, t2, n),
      ps = @remote(b, CurvePointsAt(r, ti)),
      ts = @remote(b, CurveTangentsAt(r, ti)),
      #ns = ACADCurveNormalsAt(conn, r, ti),
      frames = rotation_minimizing_frames(@remote(b, CurveFrameAt(r, t1)), ps, ts)
    map(f, frames)
  end

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

realize(b::RH, s::SurfaceArc) =
  @remote(b, SurfaceArc(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) =
  @remote(b, SurfaceEllipse(s.center, vz(1, s.center.cs), s.radius_x, s.radius_y))

backend_surface_polygon(b::RH, vs::Locs) =
  @remote(b, SurfaceClosedPolyLine(vs))

realize(b::RH, s::Surface) =
  let ids = @remote(b, SurfaceFrom(collect_ref(s.frontier)))
    foreach(mark_deleted, s.frontier)
    ids
  end
#=
realize(b::RH, s::Text) =
  @remote(b, Text(s.str, s.c, vx(1, s.c.cs), vy(1, s.c.cs), s.h))
=#

backend_surface_domain(b::RH, s::Shape2D) =
    tuple(@remote(b, SurfaceDomain(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) = @remote(b, SurfaceDomain(r))
        map_division(u1, u2, nu) do u
            map_division(v1, v2, nv) do v
                f(@remote(b, SurfaceFrameAt(r, u, v)))
            end
        end
    end

realize(b::RH, s::Sphere) =
    @remote(b, Sphere(s.center, s.radius))
realize(b::RH, s::Torus) =
    @remote(b, Torus(s.center, vz(1, s.center.cs), s.re, s.ri))
realize(b::RH, s::Cuboid) =
    @remote(b, IrregularPyramidFrustum([s.b0, s.b1, s.b2, s.b3], [s.t0, s.t1, s.t2, s.t3]))
realize(b::RH, s::RegularPyramidFrustum) =
    @remote(b, IrregularPyramidFrustum(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)))
backend_pyramid(b::RH, bs::Locs, t::Loc) =
  @remote(b, IrregularPyramid(bs, t))
backend_pyramid_frustum(b::RH, bs::Locs, ts::Locs) =
  @remote(b, IrregularPyramidFrustum(bs, ts))
backend_right_cuboid(b::RH, cb, width, height, h, material) =
  @remote(b, XYCenteredBox(cb, vx(1, cb.cs), vy(1, cb.cs), width, height, h))

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

realize(b::RH, s::Cone) =
  @remote(b, Cone(add_z(s.cb, s.h), s.r, s.cb))

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

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

backend_extrusion(b::RH, s::Shape, v::Vec) =
    and_mark_deleted(
        map_ref(s) do r
            @remote(b, Extrusion(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
      @remote(b, SweepPathProfile(path_r, profile_r, rotation, scale))
    end
  end

#=
realize(b::RH, s::Revolve) =
  and_delete_shape(
    map_ref(s.profile) do r
      @remote(b, Revolve(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) =
  (@remote(b, Unite(r0.value, r1.value)); r0)

intersect_ref(b::RH, r0::RHNativeRef, r1::RHNativeRef) =
    let refs = @remote(b, Intersect(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 = try @remote(b, Subtract(r0.value, r1.value)); catch e; [nothing] end,
      n = length(refs)
    if n == 0
        RHEmptyRef()
    elseif n == 1
        if isnothing(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::SubtractionRef, r1::RHNativeRef) =
  subtract_ref(b,
    subtract_ref(b, r0.value, r1),
    length(r0.values) == 1 ? r0.values[1] : RHUnionRef(r0.values))

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) =
  (@remote(b, Slice(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
            @remote(b, Move(r, s.v))
            r
          end
    mark_deleted(s.shape)
    r
  end

realize(b::RH, s::Scale) =
  let r = map_ref(s.shape) do r
            @remote(b, Scale(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
            @remote(b, Rotate(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
                    @remote(b, Mirror(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
            @remote(b, Mirror(r, s.p, s.n, true))
          end
    UnionRef((r0,r1))
  end

backend_surface_grid(b::RH, points, closed_u, closed_v, smooth_u, smooth_v) =
  let (nu, nv) = size(points),
      order(n) = min(max(2*floor(Int,n/30) + 1, 2), 5)
    @remote(b, SurfaceFromGrid(nu, nv,
                               reshape(permutedims(points), :),
                               closed_u, closed_v,
                               smooth_u ? order(nu) : 1,
                               smooth_v ? order(nv) : 1))
  end

realize(b::RH, s::Thicken) =
  and_delete_shape(
    map_ref(s.shape) do r
      @remote(b, Thicken(r, s.thickness))
    end,
    s.shape)

# BIM

# Families

realize(b::RH, f::TableFamily) =
    @remote(b, CreateRectangularTableFamily(f.length, f.width, f.height, f.top_thickness, f.leg_thickness))
realize(b::RH, f::ChairFamily) =
    @remote(b, CreateChairFamily(f.length, f.width, f.height, f.seat_height, f.thickness))
realize(b::RH, f::TableChairFamily) =
    @remote(b, CreateRectangularTableAndChairsFamily(
        realize(b, f.table_family), realize(b, 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) =
    @remote(b, Table(c, angle, realize(b, family)))

backend_chair(b::RH, c, angle, family) =
    @remote(b, Chair(c, angle, realize(b, family)))

backend_rectangular_table_and_chairs(b::RH, c, angle, family) =
    @remote(b, TableAndChairs(c, angle, realize(b, family)))

backend_slab(b::RH, profile, holes, thickness, family) =
    map_ref(b,
            r->@remote(b, Extrusion(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)
    @remote(b, XYCenteredBox(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)
    @remote(b, XYCenteredBox(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)
    @remote(b, XYCenteredBox(o, vx(1, o.cs), vy(1, o.cs), profile.dx, profile.dy, height))
  end
=#

backend_wall(b::RH, path, height, l_thickness, r_thickness, family) =
    @remote(b, PathWall(backend_stroke(b, offset(path, l_thickness - r_thickness)), #offset(path, (l_thickness - r_thickness)/2)),
                        l_thickness, r_thickness,
                        height))

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

backend_bounding_box(b::RH, shapes::Shapes) =
  @remote(b, BoundingBox(collect_ref(shapes)))

import Base.view
set_view(camera::XYZ, target::XYZ, lens::Real, aperture::Real, b::RH) =
  @remote(b, View(camera, target, lens))

get_view(b::RH) =
  @remote(b, ViewCamera()), @remote(b, ViewTarget()), @remote(b, ViewLens())

zoom_extents(b::RH) =
  @remote(b, ZoomExtents())

view_top(b::RH) =
  @remote(b, ViewTop())

backend_delete_shapes(b::RH, shapes::Shapes) =
  @remote(b, DeleteMany(collect_ref(shapes)))

delete_all_shapes(b::RH) =
  @remote(b, DeleteAll())

# Layers
RHLayer = String

current_layer(b::RH)::RHLayer =
  @remote(b, CurrentLayer())

current_layer(layer::RHLayer, b::RH) =
  @remote(b, SetCurrentLayer(layer))

backend_create_layer(b::RH, name::String, active::Bool, color::RGB) =
  let to255(x) = round(UInt8, x*255)
    @remote(b, CreateLayer(name, true, to255(red(color)), to255(green(color)), to255(blue(color))))
  end

delete_all_shapes_in_layer(layer::RHLayer, b::RH) =
  @remote(b, DeleteAllInLayer(layer))

shape_from_ref(r, b::RH) =
  let c = connection(b)
    let code = @remote(b, ShapeCode(r)),
        ref = LazyRef(b, RHNativeRef(r))
      if code == 1
        point(@remote(b, PointPosition(r)),
              backend=b, ref=ref)
      elseif code == 2
        circle(maybe_loc_from_o_vz(@remote(b, CircleCenter(r)), @remote(b, CircleNormal(r))),
               @remote(b, CircleRadius(r)),
               backend=b, ref=ref)
      elseif 3 <= code <= 6
        line(@remote(b, LineVertices(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(@remote(b, LineVertices(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

#

select_position(prompt::String, b::RH) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = @remote(b, GetPosition(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 = @remote(b, GetPosition(prompt)) |> length)) > 0
        push!(ps, p...)
    end
    ps
  end

# 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, @get_remote b GetPoint)

select_points(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, @get_remote b GetPoints)

select_curve(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, @get_remote b GetCurve)

select_curves(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, @get_remote b GetCurves)

select_surface(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, @get_remote b GetSurface)

select_surfaces(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, @get_remote b GetSurfaces)

select_solid(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, @get_remote b GetSolid)

select_solids(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, @get_remote b GetSolids)

select_shape(prompt::String, b::RH) =
  select_one_with_prompt(prompt, b, @get_remote b GetShape)

select_shapes(prompt::String, b::RH) =
  select_many_with_prompt(prompt, b, @get_remote b GetShapes)

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(rhino, [")
    for s in ss
      print(ref(s).value)
      print(", ")
    end
    println("])")
  end

all_shapes(b::RH) =
  Shape[shape_from_ref(r, b)
        for r in filter(r -> @remote(b, ShapeCode(r)) != 0, @remote(b, GetAllShapes()))]

all_shapes_in_layer(layer, b::RH) =
  Shape[shape_from_ref(r, b) for r in @remote(b, GetAllShapesInLayer(layer))]

render_view(path::String, b::RH) =
  @remote(b, Render(render_width(), render_height(), path))

save_view(path::String, b::RH) =
  @remote(b, SaveView(render_width(), render_height(), path))

#=
BIM families for Rhino
=#

abstract type RhinoFamily <: Family end

struct RhinoLayerFamily <: RhinoFamily
  name::String
  color::RGB
  ref::Parameter{Any}
end

rhino_layer_family(name, color::RGB=rgb(1,1,1)) =
  RhinoLayerFamily(name, color, Parameter{Any}(nothing))

backend_get_family_ref(b::RH, f::Family, af::RhinoLayerFamily) =
  backend_create_layer(b, af.name, true, af.color)

set_backend_family(default_wall_family(), rhino, rhino_layer_family("Wall"))
set_backend_family(default_slab_family(), rhino, rhino_layer_family("Slab"))
set_backend_family(default_roof_family(), rhino, rhino_layer_family("Roof"))
set_backend_family(default_beam_family(), rhino, rhino_layer_family("Beam"))
set_backend_family(default_column_family(), rhino, rhino_layer_family("Column"))
set_backend_family(default_door_family(), rhino, rhino_layer_family("Door"))
set_backend_family(default_panel_family(), rhino, rhino_layer_family("Panel"))

set_backend_family(default_table_family(), rhino, rhino_layer_family("Table"))
set_backend_family(default_chair_family(), rhino, rhino_layer_family("Chair"))
set_backend_family(default_table_chair_family(), rhino, rhino_layer_family("TableChairs"))

set_backend_family(default_curtain_wall_family(), rhino, rhino_layer_family("CurtainWall"))
set_backend_family(default_curtain_wall_family().panel, rhino, rhino_layer_family("CurtainWall-Panel"))
set_backend_family(default_curtain_wall_family().boundary_frame, rhino, rhino_layer_family("CurtainWall-Boundary"))
set_backend_family(default_curtain_wall_family().transom_frame, rhino, rhino_layer_family("CurtainWall-Transom"))
set_backend_family(default_curtain_wall_family().mullion_frame, rhino, rhino_layer_family("CurtainWall-Mullion"))
#current_backend(rhino)

set_backend_family(default_truss_node_family(), rhino, rhino_layer_family("TrussNodes"))
set_backend_family(default_truss_bar_family(), rhino, rhino_layer_family("TrussBars"))


text:
using SailRoute, Test

@testset "Distances" begin
    dist, bearing = SailRoute.haversine(-99.436554, 41.507483, -98.315949, 38.504048)
    @test dist ≈ 187.595 atol=0.01
    dist, bearing = SailRoute.euclidean(0.0, 0.0, 10.0, 10.0)
    @test dist ≈ 14.142 atol=0.01
    @test bearing ≈ 45.0
end

ext: jl


text:
@testset "irbuilder" begin

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule"; ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx), [LLVM.Int32Type(ctx), LLVM.Int32Type(ctx),
                                                LLVM.FloatType(ctx), LLVM.FloatType(ctx),
                                                LLVM.PointerType(LLVM.Int32Type(ctx)),
                                                LLVM.PointerType(LLVM.Int32Type(ctx))])
    fn = LLVM.Function(mod, "SomeFunction", ft)

    entrybb = BasicBlock(fn, "entry"; ctx)
    position!(builder, entrybb)
    @assert position(builder) == entrybb

    @test debuglocation(builder) === nothing
    loc = DILocation(ctx, 1, 1)
    debuglocation!(builder, loc)
    @test debuglocation(builder) == loc
    debuglocation!(builder)
    @test debuglocation(builder) == nothing

    retinst1 = ret!(builder)
    @check_ir retinst1 "ret void"
    debuglocation!(builder, retinst1)

    retinst2 = ret!(builder, ConstantInt(LLVM.Int32Type(ctx), 0))
    @check_ir retinst2 "ret i32 0"

    retinst3 = ret!(builder, Value[])
    @check_ir retinst3 "ret void undef"

    thenbb = BasicBlock(fn, "then"; ctx)
    elsebb = BasicBlock(fn, "else"; ctx)

    brinst1 = br!(builder, thenbb)
    @check_ir brinst1 "br label %then"

    cond1 = isnull!(builder, parameters(fn)[1], "cond")
    brinst2 = br!(builder, cond1, thenbb, elsebb)
    @check_ir brinst2 "br i1 %cond, label %then, label %else"

    resumeinst = resume!(builder, UndefValue(LLVM.Int32Type(ctx)))
    @check_ir resumeinst "resume i32 undef"

    unreachableinst = unreachable!(builder)
    @check_ir unreachableinst "unreachable"

    int1 = parameters(fn)[1]
    int2 = parameters(fn)[2]

    float1 = parameters(fn)[3]
    float2 = parameters(fn)[4]

    binopinst = binop!(builder, LLVM.API.LLVMAdd, int1, int2)
    @check_ir binopinst "add i32 %0, %1"

    addinst = add!(builder, int1, int2)
    @check_ir addinst "add i32 %0, %1"

    nswaddinst = nswadd!(builder, int1, int2)
    @check_ir nswaddinst "add nsw i32 %0, %1"

    nuwaddinst = nuwadd!(builder, int1, int2)
    @check_ir nuwaddinst "add nuw i32 %0, %1"

    faddinst = fadd!(builder, float1, float2)
    @check_ir faddinst "fadd float %2, %3"

    subinst = sub!(builder, int1, int2)
    @check_ir subinst "sub i32 %0, %1"

    nswsubinst = nswsub!(builder, int1, int2)
    @check_ir nswsubinst "sub nsw i32 %0, %1"

    nuwsubinst = nuwsub!(builder, int1, int2)
    @check_ir nuwsubinst "sub nuw i32 %0, %1"

    fsubinst = fsub!(builder, float1, float2)
    @check_ir fsubinst "fsub float %2, %3"

    mulinst = mul!(builder, int1, int2)
    @check_ir mulinst "mul i32 %0, %1"

    nswmulinst = nswmul!(builder, int1, int2)
    @check_ir nswmulinst "mul nsw i32 %0, %1"

    nuwmulinst = nuwmul!(builder, int1, int2)
    @check_ir nuwmulinst "mul nuw i32 %0, %1"

    fmulinst = fmul!(builder, float1, float2)
    @check_ir fmulinst "fmul float %2, %3"

    udivinst = udiv!(builder, int1, int2)
    @check_ir udivinst "udiv i32 %0, %1"

    sdivinst = sdiv!(builder, int1, int2)
    @check_ir sdivinst "sdiv i32 %0, %1"

    exactsdivinst = exactsdiv!(builder, int1, int2)
    @check_ir exactsdivinst "sdiv exact i32 %0, %1"

    fdivinst = fdiv!(builder, float1, float2)
    @check_ir fdivinst "fdiv float %2, %3"

    ureminst = urem!(builder, int1, int2)
    @check_ir ureminst "urem i32 %0, %1"

    sreminst = srem!(builder, int1, int2)
    @check_ir sreminst "srem i32 %0, %1"

    freminst = frem!(builder, float1, float2)
    @check_ir freminst "frem float %2, %3"

    shlinst = shl!(builder, int1, int2)
    @check_ir shlinst "shl i32 %0, %1"

    lshrinst = lshr!(builder, int1, int2)
    @check_ir lshrinst "lshr i32 %0, %1"

    ashrinst = ashr!(builder, int1, int2)
    @check_ir ashrinst "ashr i32 %0, %1"

    andinst = and!(builder, int1, int2)
    @check_ir andinst "and i32 %0, %1"

    orinst = or!(builder, int1, int2)
    @check_ir orinst "or i32 %0, %1"

    xorinst = xor!(builder, int1, int2)
    @check_ir xorinst "xor i32 %0, %1"

    allocainst = alloca!(builder, LLVM.Int32Type(ctx))
    @check_ir allocainst "alloca i32"

    array_allocainst = array_alloca!(builder, LLVM.Int32Type(ctx), int1)
    @check_ir array_allocainst "alloca i32, i32 %0"

    # mallocinst = malloc!(builder, LLVM.Int32Type(ctx))
    # @check_ir mallocinst "bitcast i8* %malloccall to i32*"

    ptr1 = parameters(fn)[5]

    freeinst = free!(builder, ptr1)
    @check_ir freeinst "tail call void @free"

    loadinst = load!(builder, ptr1)
    @check_ir loadinst "load i32, i32* %4"
    alignment!(loadinst, 4)
    @test alignment(loadinst) == 4

    storeinst = store!(builder, int1, ptr1)
    @check_ir storeinst "store i32 %0, i32* %4"

    fenceinst = fence!(builder, LLVM.API.LLVMAtomicOrderingNotAtomic)
    @check_ir fenceinst "fence"

    gepinst = gep!(builder, ptr1, [int1])
    @check_ir gepinst "getelementptr i32, i32* %4, i32 %0"

    gepinst1 = inbounds_gep!(builder, ptr1, [int1])
    @check_ir gepinst1 "getelementptr inbounds i32, i32* %4, i32 %0"

    truncinst = trunc!(builder, int1, LLVM.Int16Type(ctx))
    @check_ir truncinst "trunc i32 %0 to i16"

    zextinst = zext!(builder, int1, LLVM.Int64Type(ctx))
    @check_ir zextinst "zext i32 %0 to i64"

    sextinst = sext!(builder, int1, LLVM.Int64Type(ctx))
    @check_ir sextinst "sext i32 %0 to i64"

    fptouiinst = fptoui!(builder, float1, LLVM.Int32Type(ctx))
    @check_ir fptouiinst "fptoui float %2 to i32"

    fptosiinst = fptosi!(builder, float1, LLVM.Int32Type(ctx))
    @check_ir fptosiinst "fptosi float %2 to i32"

    uitofpinst = uitofp!(builder, int1, LLVM.FloatType(ctx))
    @check_ir uitofpinst "uitofp i32 %0 to float"

    sitofpinst = sitofp!(builder, int1, LLVM.FloatType(ctx))
    @check_ir sitofpinst "sitofp i32 %0 to float"

    fptruncinst = fptrunc!(builder, float1, LLVM.HalfType(ctx))
    @check_ir fptruncinst "fptrunc float %2 to half"

    fpextinst = fpext!(builder, float1, LLVM.DoubleType(ctx))
    @check_ir fpextinst "fpext float %2 to double"

    ptrtointinst = ptrtoint!(builder, ptr1, LLVM.Int32Type(ctx))
    @check_ir ptrtointinst "ptrtoint i32* %4 to i32"

    inttoptrinst = inttoptr!(builder, int1, LLVM.PointerType(LLVM.Int32Type(ctx)))
    @check_ir inttoptrinst "inttoptr i32 %0 to i32*"

    bitcastinst = bitcast!(builder, int1, LLVM.FloatType(ctx))
    @check_ir bitcastinst "bitcast i32 %0 to float"

    ptr1typ = llvmtype(ptr1)
    ptr2typ = LLVM.PointerType(ptr1typ, 2)

    addrspacecastinst = addrspacecast!(builder, ptr1, ptr2typ)
    @check_ir addrspacecastinst "addrspacecast i32* %4 to i32* addrspace(2)*"

    zextorbitcastinst = zextorbitcast!(builder, int1, LLVM.FloatType(ctx))
    @check_ir zextorbitcastinst "bitcast i32 %0 to float"

    sextorbitcastinst = sextorbitcast!(builder, int1, LLVM.FloatType(ctx))
    @check_ir sextorbitcastinst "bitcast i32 %0 to float"

    truncorbitcastinst = truncorbitcast!(builder, int1, LLVM.FloatType(ctx))
    @check_ir truncorbitcastinst "bitcast i32 %0 to float"

    castinst = cast!(builder, LLVM.API.LLVMBitCast, int1, LLVM.FloatType(ctx))
    @check_ir castinst "bitcast i32 %0 to float"

    ptr3typ = LLVM.PointerType(LLVM.FloatType(ctx))

    pointercastinst = pointercast!(builder, ptr1, ptr3typ)
    @check_ir pointercastinst "bitcast i32* %4 to float*"

    intcastinst = intcast!(builder, int1, LLVM.Int64Type(ctx))
    @check_ir intcastinst "sext i32 %0 to i64"

    fpcastinst = fpcast!(builder, float1, LLVM.DoubleType(ctx))
    @check_ir fpcastinst "fpext float %2 to double"

    icmpinst = icmp!(builder, LLVM.API.LLVMIntEQ, int1, int2)
    @check_ir icmpinst "icmp eq i32 %0, %1"

    fcmpinst = fcmp!(builder, LLVM.API.LLVMRealOEQ, float1, float2)
    @check_ir fcmpinst "fcmp oeq float %2, %3"

    phiinst = phi!(builder, LLVM.Int32Type(ctx))
    @check_ir phiinst "phi i32 "

    selectinst = LLVM.select!(builder, cond1, int1, int2)
    @check_ir selectinst "select i1 %cond, i32 %0, i32 %1"

    trap = LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    callinst = call!(builder, trap)
    @check_ir callinst "call void @llvm.trap()"
    @test called_value(callinst) == trap

    neginst = neg!(builder, int1)
    @check_ir neginst "sub i32 0, %0"

    nswneginst = nswneg!(builder, int1)
    @check_ir nswneginst "sub nsw i32 0, %0"

    nuwneginst = nuwneg!(builder, int1)
    @check_ir nuwneginst "sub nuw i32 0, %0"

    fneginst = fneg!(builder, float1)
    @check_ir fneginst "fneg float %2"

    notinst = not!(builder, int1)
    @check_ir notinst "xor i32 %0, -1"

    strinst = globalstring!(builder, "foobar")
    @check_ir strinst "private unnamed_addr constant [7 x i8] c\"foobar\\00\""

    strptrinst = globalstring_ptr!(builder, "foobar")
    @check_ir strptrinst "i8* getelementptr inbounds ([7 x i8], [7 x i8]* @1, i32 0, i32 0)"

    isnullinst = isnull!(builder, int1)
    @check_ir isnullinst "icmp eq i32 %0, 0"

    isnotnullinst = isnotnull!(builder, int1)
    @check_ir isnotnullinst "icmp ne i32 %0, 0"

    ptr2 = parameters(fn)[6]

    ptrdiffinst = ptrdiff!(builder, ptr1, ptr2)
    @check_ir ptrdiffinst "sdiv exact i64 %71, ptrtoint (i32* getelementptr (i32, i32* null, i32 1) to i64)"

    position!(builder)
end
end
end

end


text:
using CSV
using JuMP
using DataFrames
using Ipopt

function calibrate(d::Dict, set::Dict; save = true, overwrite = false)

    io_cal = read_build("calibrate"; save = save, overwrite = overwrite)
    !isempty(io_cal) && (return io_cal)

    # Copy the relevant input DataFrames before making any changes.
    set[:cal] = [:a0,:fd0,:fs0,:id0,:m0,:md0,:ms0,:ta0,:tm0,:va0,:x0,:y0,:ys0]
    d = Dict(k => copy(d[k]) for k in set[:cal])

    # Set all values to be at least zero for final demand and tax rates.
    [d[k][d[k][:,:value] .< 0, :value] .= 0 for k in [:fd0, :ta0, :tm0]]

    # Initialize a DataFrame to contain results.
    # io_cal = Dict(k => DataFrame() for k in setdiff(set[:cal], [:ta0,:tm0]))
    io_cal = Dict(k => DataFrame() for k in set[:cal])

    for year in set[:yr]
        io_cal_temp = calibrate(year, d, set)
        [io_cal[k] = [io_cal[k]; io_cal_temp[k]] for k in keys(io_cal)]
    end

    write_build("calibrate", io_cal; save = save)

    return io_cal
end

function calibrate(year::Int, d::Dict, set::Dict)
    println("  Calibrating $year data")

    # Initialize resultant DataFrame.
    io_cal = Dict(k => edit_with(filter_with(copy(d[k]), (yr = year,)),
        Rename(:i,:g)) for k in [:ta0,:tm0])

    cal = Dict(k => convert_type(Dict, edit_with(filter_with(d[k],
        (yr = year,)), Drop.([:yr,:value], ["all", 0.0], "=="))) for k in set[:cal])
    cal[:tm0] = fill_zero((set[:i],), cal[:tm0])
    cal[:ta0] = fill_zero((set[:i],), cal[:tm0])

    # with_optimizer(Ipopt.Optimizer, max_cpu_time=60.0)` becomes `optimizer_with_attributes(Ipopt.Optimizer, "max_cpu_time" => 60.0)
    # calib = Model(with_optimizer(Ipopt.Optimizer, nlp_scaling_method="gradient-based"))
    calib = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_cpu_time" => 60.0))

    @variable(calib,ys0_est[j in set[:j], i in set[:i]] >= 0,start = 0);
    @variable(calib,fs0_est[i in set[:i]] >= 0,start = 0);
    @variable(calib,ms0_est[i in set[:i],m in set[:m]] >= 0,start = 0);
    @variable(calib,y0_est[i in set[:i]] >= 0,start = 0);
    @variable(calib,id0_est[i in set[:i],j in set[:j]] >= 0,start = 0);
    @variable(calib,fd0_est[i in set[:i],fd in set[:fd]] >= 0,start = 0);
    @variable(calib,va0_est[va in set[:va],j in set[:j]] >= 0,start = 0);
    @variable(calib,a0_est[i in set[:i]] >= 0,start = 0);
    @variable(calib,x0_est[i in set[:i]] >= 0,start = 0);
    @variable(calib,m0_est[i in set[:i]] >= 0,start = 0);
    @variable(calib,md0_est[m in set[:m],i in set[:i]] >= 0,start = 0);

    # ```
    # mkt_py(i)..	sum(j, ys0_(j,i)) +  fs0_(i) =E= sum(m, ms0_(i,m)) + y0_(i);
    # mkt_pa(i)..	a0_(i) =E= sum(j, id0_(i,j)) + sum(fd,fd0_(i,fd));
    # mkt_pm(m)..	sum(i,ms0_(i,m)) =E= sum(i, md0_(m,i));
    # prf_y(j)..	sum(i, ys0_(j,i)) =E= sum(i, id0_(i,j)) + sum(va,va0_(va,j));
    # prf_a(i)..	a0_(i)*(1-ta0(i)) + x0_(i) =E= y0_(i) + m0_(i)*(1+tm0(i)) + sum(m, md0_(m,i));
    # ```

    @constraint(calib,mkt_py[i in set[:i]],
        sum(ys0_est[j,i] for j in set[:j]) + fs0_est[i] == sum(ms0_est[i,m] for m in set[:m]) + y0_est[i]
    );

    @constraint(calib,mkt_pa[i in set[:i]],
        a0_est[i] == sum(id0_est[i,j] for j in set[:j]) + sum(fd0_est[i,fd] for fd in set[:fd])
    );

    @constraint(calib,mkt_pm[m in set[:m]],
        sum(ms0_est[i,m] for i in set[:i]) == sum(md0_est[m,i] for i in set[:i])
    );

    @constraint(calib,prf_y[j in set[:j]],
        sum(ys0_est[j,i] for i in set[:i]) == sum(id0_est[i,j] for i in set[:i]) + sum(va0_est[va,j] for va in set[:va])
    );

    @constraint(calib,prf_a[i in set[:i]],
        a0_est[i] * (1 - cal[:ta0][i]) + x0_est[i] == y0_est[i] + m0_est[i] * (1 + cal[:tm0][i]) + sum(md0_est[m,i] for m in set[:m])
    );

    penalty_nokey = 1e4

    @objective(calib, Min,
        + sum(abs(cal[:ys0][j,i]) * (ys0_est[j,i] / cal[:ys0][j,i] - 1)^2 for i in set[:i] for j in set[:j] if haskey(cal[:ys0], (j, i)))
        + sum(abs(cal[:id0][i,j]) * (id0_est[i,j] / cal[:id0][i,j] - 1)^2 for i in set[:i] for j in set[:j] if haskey(cal[:id0], (i, j)))
        + sum(abs(cal[:fs0][i]) * (fs0_est[i] / cal[:fs0][i] - 1)^2 for i in set[:i] if haskey(cal[:fs0], i))
        + sum(abs(cal[:ms0][i,m]) * (ms0_est[i,m] / cal[:ms0][i,m] - 1)^2 for i in set[:i] for m in set[:m] if haskey(cal[:ms0], (i, m)))
        + sum(abs(cal[:y0][i]) * (y0_est[i] / cal[:y0][i] - 1)^2 for i in set[:i] if haskey(cal[:y0], i))
        + sum(abs(cal[:fd0][i,fd]) * (fd0_est[i,fd] / cal[:fd0][i,fd] - 1)^2 for i in set[:i] for fd in set[:fd] if haskey(cal[:fd0], (i, fd)))
        + sum(abs(cal[:va0][va,j]) * (va0_est[va,j] / cal[:va0][va,j] - 1)^2 for va in set[:va] for j in set[:j] if haskey(cal[:va0], (va, j)))
        + sum(abs(cal[:a0][i]) * (a0_est[i] / cal[:a0][i] - 1)^2  for i in set[:i] if haskey(cal[:a0], i))
        + sum(abs(cal[:x0][i]) * (x0_est[i] / cal[:x0][i] - 1)^2  for i in set[:i] if haskey(cal[:x0], i))
        + sum(abs(cal[:m0][i]) * (m0_est[i] / cal[:m0][i] - 1)^2  for i in set[:i] if haskey(cal[:m0], i))
        + sum(abs(cal[:md0][m,i]) * (md0_est[m,i] / cal[:md0][m,i] - 1)^2 for m in set[:m] for i in set[:i] if haskey(cal[:md0], (m, i)))

    + penalty_nokey * (
        + sum(ys0_est[j,i]^2 for i in set[:i] for j in set[:j] if !haskey(cal[:ys0], (j, i)))
        + sum(id0_est[i,j]^2  for i in set[:i] for j in set[:j] if !haskey(cal[:id0], (i, j)))
        + sum(fs0_est[i]^2  for i in set[:i] if !haskey(cal[:fs0], i))
        + sum(ms0_est[i,m]^2  for i in set[:i] for m in set[:m] if !haskey(cal[:ms0], (i, m)))
        + sum(y0_est[i]^2 for i in set[:i] if !haskey(cal[:y0], i))
        + sum(fd0_est[i,fd]^2 for i in set[:i] for fd in set[:fd] if !haskey(cal[:fd0], (i, fd)))
        + sum(va0_est[va,j]^2 for va in set[:va] for j in set[:j] if !haskey(cal[:va0], (va, j)))
        + sum(a0_est[i]^2 for i in set[:i] if !haskey(cal[:a0], i))
        + sum(x0_est[i]^2 for i in set[:i] if !haskey(cal[:x0], i))
        + sum(m0_est[i]^2 for i in set[:i] if !haskey(cal[:m0], i))
        + sum(md0_est[m,i]^2  for m in set[:m] for i in set[:i] if !haskey(cal[:md0], (m, i)))
        )
    );

    [set_start_value(ys0_est[j,i], cal[:ys0][j,i]) for i in set[:i] for j in set[:j] if haskey(cal[:ys0], (j, i)) ] ;
    [set_start_value(id0_est[i,j], cal[:id0][i,j]) for i in set[:i] for j in set[:j] if haskey(cal[:id0], (i, j)) ] ;
    [set_start_value(fs0_est[i],   cal[:fs0][i]) for i in set[:i] if haskey(cal[:fs0], i) ];
    [set_start_value(ms0_est[i,m], cal[:ms0][i,m]) for i in set[:i] for m in set[:m] if haskey(cal[:ms0], (i, m)) ];
    [set_start_value(y0_est[i],    cal[:y0][i]) for i in set[:i] if haskey(cal[:y0], i) ];
    [set_start_value(fd0_est[i,fd],cal[:fd0][i,fd]) for i in set[:i] for fd in set[:fd] if haskey(cal[:fd0], (i, fd)) ];
    [set_start_value(va0_est[va,j],cal[:va0][va,j]) for va in set[:va] for j in set[:j] if haskey(cal[:va0], (va, j)) ];
    [set_start_value(a0_est[i],    cal[:a0][i]) for i in set[:i] if haskey(cal[:a0], (i)) ];
    [set_start_value(x0_est[i],    cal[:x0][i]) for i in set[:i] if haskey(cal[:x0], i) ];
    [set_start_value(m0_est[i],    cal[:m0][i]) for i in set[:i] if haskey(cal[:m0], i) ];
    [set_start_value(md0_est[m,i], cal[:md0][m,i]) for m in set[:m] for i in set[:i] if haskey(cal[:md0], (m, i)) ] ;

    # multipliers for lower and upper bound relative
    # to each respective variables reference parameter
    lb = 0.1
    ub = 5

    [set_lower_bound(ys0_est[j,i], max(0, lb * cal[:ys0][j,i])) for i in set[:i] for j in set[:j] if haskey(cal[:ys0], (j, i)) ] ;
    [set_lower_bound(id0_est[i,j], max(0, lb * cal[:id0][i,j])) for i in set[:i] for j in set[:j] if haskey(cal[:id0], (i, j)) ] ;
    [set_lower_bound(fs0_est[i], max(0, lb * cal[:fs0][i])) for i in set[:i] if haskey(cal[:fs0], i) ];
    [set_lower_bound(ms0_est[i,m], max(0, lb * cal[:ms0][i,m])) for i in set[:i] for m in set[:m] if haskey(cal[:ms0], (i, m)) ];
    [set_lower_bound(y0_est[i], max(0, lb * cal[:y0][i])) for i in set[:i] if haskey(cal[:y0], i) ];
    [set_lower_bound(fd0_est[i,fd], max(0, lb * cal[:fd0][i,fd])) for i in set[:i] for fd in set[:fd] if haskey(cal[:fd0], (i, fd)) ];
    [set_lower_bound(va0_est[va,j], max(0, lb * cal[:va0][va,j])) for va in set[:va] for j in set[:j] if haskey(cal[:va0], (va, j)) ];
    [set_lower_bound(a0_est[i], max(0, lb * cal[:a0][i])) for i in set[:i] if haskey(cal[:a0], (i)) ];
    [set_lower_bound(x0_est[i], max(0, lb * cal[:x0][i])) for i in set[:i] if haskey(cal[:x0], i) ];
    [set_lower_bound(m0_est[i], max(0, lb * cal[:m0][i])) for i in set[:i] if haskey(cal[:m0], i) ];
    [set_lower_bound(md0_est[m,i], max(0, lb * cal[:md0][m,i])) for m in set[:m] for i in set[:i] if haskey(cal[:md0], (m, i)) ] ;

    [set_upper_bound(ys0_est[j,i], abs(ub * cal[:ys0][j,i])) for i in set[:i] for j in set[:j] if haskey(cal[:ys0], (j, i)) ] ;
    [set_upper_bound(id0_est[i,j], abs(ub * cal[:id0][i,j])) for i in set[:i] for j in set[:j] if haskey(cal[:id0], (i, j)) ] ;
    [set_upper_bound(fs0_est[i], abs(ub * cal[:fs0][i])) for i in set[:i] if haskey(cal[:fs0], i) ];
    [set_upper_bound(ms0_est[i,m], abs(ub * cal[:ms0][i,m])) for i in set[:i] for m in set[:m] if haskey(cal[:ms0], (i, m)) ];
    [set_upper_bound(y0_est[i], abs(ub * cal[:y0][i])) for i in set[:i] if haskey(cal[:y0], i) ];
    [set_upper_bound(fd0_est[i,fd], abs(ub * cal[:fd0][i,fd])) for i in set[:i] for fd in set[:fd] if haskey(cal[:fd0], (i, fd)) ];
    [set_upper_bound(va0_est[va,j], abs(ub * cal[:va0][va,j])) for va in set[:va] for j in set[:j] if haskey(cal[:va0], (va, j)) ];
    [set_upper_bound(a0_est[i], abs(ub * cal[:a0][i])) for i in set[:i] if haskey(cal[:a0], (i)) ];
    [set_upper_bound(x0_est[i], abs(ub * cal[:x0][i])) for i in set[:i] if haskey(cal[:x0], i) ];
    [set_upper_bound(m0_est[i], abs(ub * cal[:m0][i])) for i in set[:i] if haskey(cal[:m0], i) ];
    [set_upper_bound(md0_est[m,i], abs(ub * cal[:md0][m,i])) for m in set[:m] for i in set[:i] if haskey(cal[:md0], (m, i)) ] ;

    [fix(fs0_est[i], cal[:fs0][i], force=true) for i in set[:i] if haskey(cal[:fs0], i)];
    [fix(fs0_est[i], 0, force=true) for i in set[:i] if !haskey(cal[:fs0], i)];

    [fix(va0_est[va,j], cal[:va0][va,j], force=true) for va in set[:va] for j in set[:j] if haskey(cal[:va0], (va, j))];
    [fix(va0_est[va,j], 0, force=true) for va in set[:va] for j in set[:j] if !haskey(cal[:va0], (va, j))];

    [fix(m0_est[i], cal[:m0][i], force=true) for i in set[:i] if haskey(cal[:m0], i)];
    [fix(m0_est[i], 0, force=true) for i in set[:i] if !haskey(cal[:m0], i)];

    # fd_temp = filter!(x->x≠"pce",set[:fd])
    # [fix(fd0_est[i,fd],cal[:fd0][i,fd],force=true) for i in set[:i] for fd in fd_temp if haskey(cal[:fd0],(i,fd))]
    # [fix(fd0_est[i,fd],0,force=true) for i in set[:i] for fd in fd_temp if !haskey(cal[:fd0],(i,fd))]

    [fix(ys0_est[j,i], 0, force=true) for j in set[:oth,:use] for i in set[:i]]

    JuMP.optimize!(calib)

    # Populate resultant Dictionary.
    # io_cal = Dict()

    # Using original (i,j) notation...
    # io_cal[:ys0] = convert_type(DataFrame, ys0_est; cols=[:j,:i])
    # io_cal[:fs0] = convert_type(DataFrame, fs0_est; cols=[:i])
    # io_cal[:ms0] = convert_type(DataFrame, ms0_est; cols=[:i,:m])
    # io_cal[:y0]  = convert_type(DataFrame, y0_est;  cols=[:i])
    # io_cal[:id0] = convert_type(DataFrame, id0_est; cols=[:i,:j])
    # io_cal[:fd0] = convert_type(DataFrame, fd0_est; cols=[:i,:fd])
    # io_cal[:va0] = convert_type(DataFrame, va0_est; cols=[:va,:j])
    # io_cal[:a0]  = convert_type(DataFrame, a0_est;  cols=[:i])
    # io_cal[:x0]  = convert_type(DataFrame, x0_est;  cols=[:i])
    # io_cal[:m0]  = convert_type(DataFrame, m0_est;  cols=[:i])
    # io_cal[:md0] = convert_type(DataFrame, md0_est; cols=[:m,:i])

    # Using (g,s) notation...
    io_cal[:ys0] = convert_type(DataFrame, ys0_est; cols=[:s,:g])
    io_cal[:fs0] = convert_type(DataFrame, fs0_est; cols=[:g])
    io_cal[:ms0] = convert_type(DataFrame, ms0_est; cols=[:s,:m])
    io_cal[:y0]  = convert_type(DataFrame, y0_est;  cols=[:s])
    io_cal[:id0] = convert_type(DataFrame, id0_est; cols=[:g,:s])
    io_cal[:fd0] = convert_type(DataFrame, fd0_est; cols=[:s,:fd])
    io_cal[:va0] = convert_type(DataFrame, va0_est; cols=[:va,:s])
    io_cal[:a0]  = convert_type(DataFrame, a0_est;  cols=[:g])
    io_cal[:x0]  = convert_type(DataFrame, x0_est;  cols=[:g])
    io_cal[:m0]  = convert_type(DataFrame, m0_est;  cols=[:g])
    io_cal[:md0] = convert_type(DataFrame, md0_est; cols=[:m,:g])

    [io_cal[k] = edit_with(df, Add(:yr, year))[:,unique([:yr; propertynames(df)])]
        for (k, df) in io_cal]
    return io_cal
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
@testset "Sphere" begin
    # Make some points between (0,0,0) and (1,1,1) and put in a point cloud
    N = 1000
    pos = rand(SVector{3,Float32}, N)
    sphere = Sphere(0.5, 0.5, 0.5, 0.1)

    sphere_search_correct = true
    for i = 1:length(pos)
        v = pos[i]
        dv = v - SVector(0.5,0.5,0.5)
        distance = sum(abs2, dv)
        expected = distance < 0.01
        actual = v ∈ sphere
        sphere_search_correct &= actual == expected
    end
    @test sphere_search_correct

    @test boundingbox(sphere) == BoundingBox(0.4, 0.4, 0.4, 0.6, 0.6, 0.6)

    @test volume(sphere) ≈ 0.04188790204786391

    @test might_intersect(
        BoundingBox(0.0, 0.0, 0.0, 0.1, 0.1, 0.1),
        Sphere(0.5, 0.5, 0.5, 0.5)
    ) == false

    @test might_intersect(
        BoundingBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0),
        Sphere(0.5, 0.5, 0.5, 0.2)
    ) == true

    @test might_intersect(
        BoundingBox(0.0, 0.0, 0.0, 0.4, 0.4, 0.4),
        Sphere(0.5, 0.5, 0.5, 0.2)
    ) == true
end


text:
using Oceananigans
using Oceananigans.Units
using OceanTurbulenceParameterEstimation
using Distributions
using DataDeps
using GLMakie

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

include("best_catke_parameters.jl")

θ = [θ_constant_Pr, θ_variable_Pr, θ_conv_variable_Pr]
Nensemble = length(θ)

function perturbation_prior(θ★, ϵ=0.1)
    L = (1 - ϵ) * θ★
    U = (1 + ϵ) * θ★
    return ScaledLogitNormal(bounds=(L, U))
end

θ★ = θ_constant_Pr
priors = NamedTuple(name => perturbation_prior(θ★[name]) for name in keys(θ★))
free_parameters = FreeParameters(priors)

suites = [
          "two_day_suite_1m",
          "four_day_suite_1m",
          "six_day_suite_1m",
         ]

end_times = [
             2days,
             4days,
             6days
            ]

for (suite, end_time) in zip(suites, end_times)
    case_path(case) = @datadep_str("$suite/$(case)_instantaneous_statistics.jld2")

    cases = [
             "free_convection",
             "weak_wind_strong_cooling",
             "strong_wind_weak_cooling",
             "strong_wind",
             "strong_wind_no_rotation",
            ]

    Δz = 1
    times = [6hours, end_time]
    field_names = (:u, :v, :T, :e)
    regrid_size = (1, 1, Int(256/Δz))
    observations = [SyntheticObservations(case_path(case); field_names, regrid_size, times) for case in cases]
    obs = first(observations)
    α = obs.metadata.parameters.thermal_expansion_coefficient
    equation_of_state = LinearEquationOfState(; α)
    buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=true)

    mixing_length = MixingLength(Cᴬu   = 0.0,
                                 Cᴬc   = 0.0,
                                 Cᴬe   = 0.0,
                                 CᴷRiᶜ = 2.0,
                                 Cᴷuʳ  = 0.0,
                                 Cᴷcʳ  = 0.0,
                                 Cᴷeʳ  = 0.0)

    catke = CATKEVerticalDiffusivity(; mixing_length)

    simulation = ensemble_column_model_simulation(observations; Nensemble, buoyancy,
                                                  architecture = GPU(),
                                                  tracers = (:T, :e),
                                                  closure = catke)

    simulation.Δt = 1.0
    progress(sim) = @info "Iter: $(iteration(sim)), time: $(prettytime(sim))"
    simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

    Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
    Qᵀ = simulation.model.tracers.T.boundary_conditions.top.condition
    dTdz = simulation.model.tracers.T.boundary_conditions.bottom.condition

    for (case, obs) in enumerate(observations)
        @show case cases[case]
        @show obs.metadata.parameters.momentum_flux
        @show obs.metadata.parameters.temperature_flux
        @show obs.metadata.parameters.buoyancy_flux
        @show f = obs.metadata.parameters.coriolis_parameter

        view(Qᵘ, :, case) .= obs.metadata.parameters.momentum_flux
        view(Qᵀ, :, case) .= obs.metadata.parameters.temperature_flux
        view(dTdz, :, case) .= obs.metadata.parameters.dθdz_deep
        view(simulation.model.coriolis, :, case) .= Ref(FPlane(f=f))
    end

    calibration = InverseProblem(observations, simulation, free_parameters)

    model = simulation.model

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, merge(model.velocities, model.tracers, model.diffusivity_fields),
                         prefix = string("catke_simulation_", suite[1:end-2], Δz, "m"),
                         schedule = SpecifiedTimes(times[1]:10minutes:end_time...),
                         force = true)

    @show simulation.output_writers[:fields]

    forward_run!(calibration, θ; suppress=false)
end

max_stars_repo_name: glwagner/LocalOceanClosureCalibration

max_issues_repo_name: glwagner/LocalOceanClosureCalibration

max_forks_repo_name: glwagner/LocalOceanClosureCalibration


text:
using DifferentialEquations
using FastClosures
"""
	comm!(commAB, A, B)
Calculate the commutator of the operators A and B, and save the result in commAB.
The operators are written in the basis { a^2 + a†^2, -i(a^2 - a†^2), 2a†a+1 }.
"""
function comm!(commAB, A, B)
	commAB[1] = 4*(A[3]*B[2] - A[2]*B[3])
	commAB[2] = 4*(A[1]*B[3] - A[3]*B[1])
	commAB[3] = 4*(A[1]*B[2] - A[2]*B[1])
	return nothing
end
"""
	Λ(t, tf)
Definite integral of the driving envelope.
"""
function Λ(t, tf)
	x = t/tf
	return (x - sin(2*pi*x)/(2*pi))
end
"""
	λ(t, tf)
Driving envelope.
"""
function λ(t, tf)
	x = t/tf
	return (1 - cos(2*pi*x))/tf
end
"""
	H(t, tf; cr_terms::Bool=true)
Hamiltonian of a parametrically driven cavity.
"""
function H(t, tf; cr_terms::Bool=true)
	λt = λ(t, tf)
	Ht = zeros(3)
	if cr_terms == true
		Ht[1] = λt*sin(4*t)/2
		Ht[2] = λt*(1 - cos(4*t))/2
		Ht[3] = λt*sin(2*t)
	else
		Ht[2] = λt/2
	end
	return Ht
end
"""
	V_I(t, tf)
Error Hamiltonian V in the interaction picture.
"""
function V_I(t, tf)
	λt = λ(t, tf)
	Λt = Λ(t, tf)

	V_It = zeros(3)
	V_It[1] = λt*(sin(2*t)*sinh(2*Λt) + sin(4*t)*cosh(2*Λt)/2)
	V_It[2] = -λt*cos(4*t)/2
	V_It[3] = λt*(sin(2*t)*cosh(2*Λt) + sin(4*t)*sinh(2*Λt)/2)
	return V_It
end
"""
	dv!(dv, v, args, t)
Heisenberg equations for v = [a, a†].
"""
function dv!(dv, v, args, t)
	H = args[1]
	tf = args[2]
	cr_terms = args[3]
	Ht = H(t)
	a = view(v, :, :, 1)
	adagger = view(v, :, :, 2)
	dv[:, :, 1] .= (-2im*Ht[1] + 2*Ht[2])*adagger .- 2im*Ht[3]*a
	dv[:, :, 2] .= (+2im*Ht[1] + 2*Ht[2])*a .+ 2im*Ht[3]*adagger
end
"""
	d2moments!(dv, v, args, t)
Differential equations for the second moments (v = [a^2 + a†^2, -i(a^2 - a†^2),
aa† + a†a]). Dissipation is taken in account. `H` is a function the returns the Hamiltonian
at time `t`.
"""
function d2moments!(dv, v, args, t)
	H = args[1]
	tf = args[2]
	κ = args[3]
	Ht = args[4]

	Ht .= H(t)
	freal = 2*Ht[2]
	fimag = -2*Ht[1]
	#greal = -κ/2.0
	gimag = -2*Ht[3]
	dv[1] = -κ*v[1] - 2*gimag*v[2] + 2*freal*v[3]
	dv[2] = 2*gimag*v[1] - κ*v[2] + 2*fimag*v[3]
	dv[3] = 2*freal*v[1] + 2*fimag*v[2] - κ*v[3] + κ
end
"""
	solve_2moments(tf, κ, cr_terms, coeffs)
Solve equations of motions for the second moments.
"""
function solve_2moments(tf, κ, cr_terms, coeffs)
	Hmod = @closure t-> H(t, tf; cr_terms=cr_terms) .+ WharmonicR(t, tf, coeffs)
	f0 = [0.0, 0.0, 1.0]
	tspan = (0.0, tf)
	args = (Hmod, tf, κ, zeros(3))
	prob = ODEProblem(d2moments!, f0, tspan, args)
	sol = solve(prob, Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-12)
	return sol(tf)
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
mutable struct AugmentedLagrangianCosts{T,C,CX,CU}
    costs::Objective{T}
    constraint_data::ConstraintsData{T,C,CX,CU}
    constraint_penalty::Vector{Vector{T}}
    constraint_penalty_matrix::Vector{Diagonal{T,Vector{T}}}
    constraint_dual::Vector{Vector{T}}
    active_set::Vector{Vector{Int}}
    constraint_tmp::Vector{Vector{T}}
    constraint_jacobian_state_tmp::Vector{Matrix{T}}
    constraint_jacobian_action_tmp::Vector{Matrix{T}}
end

function augmented_lagrangian(model::Model{T}, costs::Objective{T}, constraints::Constraints{T}) where T
    # horizon
    H = length(model) + 1
    # penalty
    constraint_penalty = [ones(c.num_constraint) for c in constraints]
    constraint_penalty_matrix = [Diagonal(ones(c.num_constraint)) for c in constraints]
    # duals
    constraint_dual = [zeros(c.num_constraint) for c in constraints]
    # active set
    active_set = [ones(Int, c.num_constraint) for c in constraints]
    # pre-allocated memory
    constraint_tmp = [zeros(c.num_constraint) for c in constraints]
    constraint_jacobian_state_tmp = [zeros(c.num_constraint, t < H ? model[t].num_state : model[H-1].num_next_state) for (t, c) in enumerate(constraints)]
    constraint_jacobian_action_tmp = [zeros(c.num_constraint, t < H ? model[t].num_action : 0) for (t, c) in enumerate(constraints)]
    data = constraint_data(model, constraints)
    AugmentedLagrangianCosts(costs,
        data,
        constraint_penalty,
        constraint_penalty_matrix,
        constraint_dual,
        active_set,
        constraint_tmp,
        constraint_jacobian_state_tmp,
        constraint_jacobian_action_tmp)
end

function cost(objective::AugmentedLagrangianCosts, states, actions, parameters)
    # costs
    J = cost(objective.costs, states, actions, parameters)

    # constraints
    c = objective.constraint_data.violations
    ρ = objective.constraint_penalty
    λ = objective.constraint_dual
    a = objective.active_set

    # horizon
    H = length(c)

    constraint!(objective.constraint_data, states, actions, parameters)
    active_set!(a, objective.constraint_data, λ)

    for t = 1:H
        J += λ[t]' * c[t]
        num_constraint = objective.constraint_data.constraints[t].num_constraint
        for i = 1:num_constraint
            if a[t][i] == 1
                J += 0.5 * ρ[t][i] * c[t][i]^2.0
            end
        end
    end

    return J
end

function active_set!(a, data::ConstraintsData, λ)
    # violations
    c = data.violations

    # horizon
    H = length(c)

    for t = 1:H
        # set all constraints active
        fill!(a[t], 1)

        # check inequality constraints
        for i in data.constraints[t].indices_inequality
            # check active-set criteria
            (c[t][i] < 0.0 && λ[t][i] == 0.0) && (a[t][i] = 0)
        end
    end
end

function augmented_lagrangian_update!(objective::AugmentedLagrangianCosts;
        scaling_penalty=10.0,
        max_penalty=1.0e12)

    # constraints
    c = objective.constraint_data.violations
    constraints = objective.constraint_data.constraints
    ρ = objective.constraint_penalty
    λ = objective.constraint_dual

    # horizon
    H = length(c)

    for t = 1:H
        num_constraint = constraints[t].num_constraint
        for i = 1:num_constraint
            λ[t][i] += ρ[t][i] * c[t][i]
            if i in constraints[t].indices_inequality
                λ[t][i] = max(0.0, λ[t][i])
            end
            ρ[t][i] = min(scaling_penalty * ρ[t][i], max_penalty)
        end
    end
end

Base.length(objective::AugmentedLagrangianCosts) = length(objective.costs)


text:
"""
    eigHAD_Distance(𝚽, 𝛌; indexEigs = 1:size(𝚽,2))

compute HAD "distance" (not really a distance) between pairwise graph Laplacian
eigenvectors, i.e., d_HAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁) = √(1 - a_HAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁)²).

# Input Arguments
- `𝚽::Matrix{Float64}`: matrix of graph Laplacian eigenvectors, 𝜙ⱼ₋₁ (j = 1,...,size(𝚽,1)).
- `𝛌::Array{Float64}`: array of eigenvalues. (ascending order)
- `indexEigs::Int`: default is all eigenvectors, indices of eigenvectors considered.

# Output Argument
- `dis::Matrix{Float64}`: the HAD distance matrix, dis[i,j] = d_HAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).
"""
function eigHAD_Distance(𝚽, 𝛌; indexEigs = 1:size(𝚽,2))
    n = length(indexEigs)
    A = eigHAD_Affinity(𝚽, 𝛌; indexEigs = indexEigs)
    dis = sqrt.(ones(n, n) - A.^2)
    dis[diagind(dis)] .= 0
    return dis
end

function eigHAD_Distance_neglog(𝚽, 𝛌; indexEigs = 1:size(𝚽,2))
    A = eigHAD_Affinity(𝚽, 𝛌; indexEigs = indexEigs)
    n = size(A,1)
    dis = zeros(n,n)
    for i = 1:n, j = 1:n
        if A[i,j] == 0
            dis[i,j] = 1e9
        else
            dis[i,j] = -log(A[i,j])
        end
    end
    return dis
end

"""
    eigHAD_Affinity(𝚽, 𝛌; indexEigs = 1:size(𝚽,2))

EIGHAD_AFFINITY compute Hadamard (HAD) affinity between pairwise graph Laplacian eigenvectors.

# Input Arguments
- `𝚽::Matrix{Float64}`: matrix of graph Laplacian eigenvectors, 𝜙ⱼ₋₁ (j = 1,...,size(𝚽,1)).
- `𝛌::Array{Float64}`: array of eigenvalues. (ascending order)
- `indexEigs::Int`: default is all eigenvectors, indices of eigenvectors considered.

# Output Argument
- `A::Matrix{Float64}`: a numEigs x numEigs affinity matrix, A[i,j] = a_HAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).
"""
function eigHAD_Affinity(𝚽, 𝛌; indexEigs = 1:size(𝚽,2))
    N, numEigs = size(𝚽,1), length(indexEigs)
    indNoDC = setdiff(indexEigs, 1) # get rid of DC component
    J = length(indNoDC)
    A = zeros(J, J)
    for a in 1:J, b in a:J
        i, j = indNoDC[a], indNoDC[b]
        hadamardProd = 𝚽[:,i] .* 𝚽[:,j]
        if norm(hadamardProd,2) < 0.01/sqrt(N)
            continue
        end
        λ, μ = 𝛌[i], 𝛌[j]
        x₀ = 1 ./ (max(λ, μ))
        # Find minimizer t
        result = optimize(t -> abs(exp(-t[1]*λ) + exp(-t[1]*μ) - 1), [x₀], BFGS());
        t = Optim.minimizer(result)[1]
        # Compute Hadamard affinity
        heatEvolution = 𝚽 * Diagonal(exp.(-t .* 𝛌)) * 𝚽' * hadamardProd
        A[a,b] = norm(heatEvolution,2) / (norm(hadamardProd,2) + 1e-6)
    end
    A = A + A'; for i in 1:J; A[i,i] /= 2; end

    if 1 in indexEigs
        # Set affinity measure of 𝜙₀ with itself to be the maximum and equals to 1.
        A = matrix_augmentation(A)
        A[1,1] = maximum(A)
    end
    return A ./ maximum(A)
end

function matrix_augmentation(A)
    m, n = size(A)
    B = zeros(m+1, n+1)
    B[2:end, 2:end] = A
    return B
end


text:
@testset "Serialized modifications" begin
    m = test_LP()

    sm = serialize_model(m, tmpfile("recon.serialized"))
    m2 = unwrap_serialized(sm)

    @test typeof(m2) == typeof(m)

    sm = serialize_model(m, tmpfile("recon.serialized"))
    m2 = remove_reaction(sm, reactions(m)[3])

    @test typeof(m2) == typeof(m)
    @test !(reactions(m)[3] in reactions(m2))
end


text:
using Revise
using PyCall
using PyPlot
using Bem2d
ow = pyimport("okada_wrapper") # from okada_wrapper import dc3dwrapper

"""
    thrustfaultfreesurface()

Comparison of surface displacements near a thrust fault dipping
at 45 degrees.  Includes both constant and quadratic elements.
"""
function thrustfaultfreesurface()
    mu = 30e9
    nu = 0.25

    # Free surface
    els = Elements(Int(1e5))
    nfreesurf = 200
    x1, y1, x2, y2 = discretizedline(-50, 0, 50, 0, nfreesurf)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "freesurf"
    end
    standardize_elements!(els)

    # 45 degree dipping fault
    nfault = 1
    x1, y1, x2, y2 = discretizedline(-1, -1, 0, 0, nfault)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "fault"
    end
    standardize_elements!(els)

    # Conviencience structures
    idx = getidxdict(els)
    partialsconst = initpartials(els)
    partialsquad = initpartials(els)

    # Constant slip fault
    partialsconst["disp"]["fault"]["freesurf"], partialsconst["stress"]["fault"]["freesurf"], partialsconst["trac"]["fault"]["freesurf"] = partialsconstdispstress(slip2dispstress, els, idx["fault"], idx["freesurf"], mu, nu)
    partialsconst["disp"]["freesurf"]["freesurf"], partialsconst["stress"]["freesurf"]["freesurf"], partialsconst["trac"]["freesurf"]["freesurf"] = partialsconstdispstress(slip2dispstress, els, idx["freesurf"], idx["freesurf"], mu, nu)
    faultslipconst = sqrt(2) / 2 * [1 ; 1]
    ufullspaceconst = partialsconst["disp"]["fault"]["freesurf"] * faultslipconst
    ufreesurfaceconst = inv(partialsconst["trac"]["freesurf"]["freesurf"]) * (partialsconst["trac"]["fault"]["freesurf"] * faultslipconst)
    xplotconst = els.xcenter[idx["freesurf"]]

    # Quadratic slip fault
    partialsquad["disp"]["fault"]["freesurf"], partialsquad["stress"]["fault"]["freesurf"], partialsquad["trac"]["fault"]["freesurf"] = partialsquaddispstress(slip2dispstress, els, idx["fault"], idx["freesurf"], mu, nu)
    partialsquad["disp"]["freesurf"]["freesurf"], partialsquad["stress"]["freesurf"]["freesurf"], partialsquad["trac"]["freesurf"]["freesurf"] = partialsquaddispstress(slip2dispstress, els, idx["freesurf"], idx["freesurf"], mu, nu)
    xplotquad = sort(els.xnodes[idx["freesurf"], :][:])
    faultslipquad = sqrt(2) ./ 2 * ones(6)
    ufullspacequad = partialsquad["disp"]["fault"]["freesurf"] * faultslipquad
    ufreesurfacequad = inv(partialsquad["trac"]["freesurf"]["freesurf"]) * (partialsquad["trac"]["fault"]["freesurf"] * faultslipquad)

    # Okada solution
    xokada = collect(LinRange(-5, 5, 1000))
    uxokada = zeros(length(xokada))
    uyokada = zeros(length(xokada))

    for i in 1:length(xokada)
        # Fault dipping at 45 degrees
        _, u, _ = ow.dc3dwrapper(
            2.0 / 3.0,
            [0, xokada[i] + 0.5, 0],
            0.5,
            45,  # 135
            [-1000, 1000],
            [-sqrt(2) / 2, sqrt(2) / 2],
            [0.0, 1.0, 0.0],
        )
        uxokada[i] = u[2]
        uyokada[i] = u[3]
    end

    fontsize = 20
    markersize = 12
    linewidth = 2.0
    close("all")
    figure(figsize = (12, 12))

    ax = subplot(2, 1, 1)
    plot(xokada, uxokada, "-k", linewidth=linewidth, label="Okada")
    plot(xplotconst, ufreesurfaceconst[1:2:end], "bx", markeredgewidth=linewidth, markersize=markersize, label = "const halfspace")
    plot(xplotquad, ufreesurfacequad[1:2:end], "r+", markeredgewidth=linewidth, markersize=markersize, label = "quad halfspace")
    gca().set_xlim([-5, 5]); gca().set_ylim([-1.0, 1.0])
    gca().set_xticks([-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5])
    gca().set_yticks([-1.00, -0.75, -0.50, -0.25, 0.0, 0.25, 0.50, 0.75, 1.00])
    legend(fontsize=fontsize)
    ax.tick_params("both", labelsize = fontsize)
    ylabel(L"$u_x$ (m)", fontsize=fontsize)

    ax = subplot(2, 1, 2)
    plot(xokada, uyokada, "-k", linewidth=linewidth, label="Okada")
    plot(xplotconst, ufreesurfaceconst[2:2:end], "bx", markeredgewidth=linewidth, markersize=markersize, label = "const halfspace")
    plot(xplotquad, ufreesurfacequad[2:2:end], "r+", markeredgewidth=linewidth, markersize=markersize, label = "quad halfspace")
    gca().set_xlim([-5, 5]); gca().set_ylim([-1.0, 1.0])
    gca().set_xticks([-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5])
    gca().set_yticks([-1.00, -0.75, -0.50, -0.25, 0.0, 0.25, 0.50, 0.75, 1.00])
    legend(fontsize=fontsize)
    ax.tick_params("both", labelsize = fontsize)
    xlabel(L"$x$ (m)", fontsize=fontsize); ylabel(L"$u_y$ (m)", fontsize=fontsize)
    show()
end
thrustfaultfreesurface()


text:
"""
    OrderedSamplingAlgorithm

An abstract type for ordered sampling algorithms.
"""
abstract type OrderedSamplingAlgorithm end

"""
    StartToEnd

An ordered sampling algorithm indicating that values should be
treated consecutively from start to finish of the dataset.
"""
struct StartToEnd <: OrderedSamplingAlgorithm end

"""
    EndToStart

An ordered sampling algorithm indicating that the values should be
treated consecutively from the end to the start of the dataset.
"""
struct EndToStart <: OrderedSamplingAlgorithm end

"""
RandPtOutwards

An ordered sampling algorithm indicating that the values should be
divided into two groups, separating the values at some midpoint
of the dataset. The two groups of values are then treated separately.
"""
struct RandPtOutwards <: OrderedSamplingAlgorithm
    midpoint_idx::Int
end

"""
    ChuncksForwards

An ordered sampling algorithm indicating that the values should be
divided into multiple (`n_chunks`) groups. The groups of values
are then treated separately, treating values from the start to the end of
each group.
"""
struct ChunksForwards <: OrderedSamplingAlgorithm
    n_chunks::Int
end

"""
    ChuncksBackwards

An ordered sampling algorithm indicating that the values should be
divided into multiple (`n_chunks`) groups. The groups of values
are then treated separately, treating values from the end to the start of
each group.
"""
struct ChunksBackwards <: OrderedSamplingAlgorithm
    n_chunks::Int
end

export
OrderedSamplingAlgorithm,
StartToEnd,
EndToStart,
RandPtOutwards,
ChunksForwards,
ChunksBackwards


text:
# dgs = data generated stefan
# only final landmark positions observed

using BridgeLandmarks
using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD2
using NPZ

workdir = @__DIR__
cd(workdir)

nshapes = 1

testshapes = npzread("match.npy.npz")
xobs0vec =  get(testshapes,"q0",0)
xobsTvec =  get(testshapes,"v",0)
nb = div(length(xobs0vec),2)

subs = 1:1:nb
xobs0 = [PointF(xobs0vec[2i-1],xobs0vec[2i]) for i in subs]
xobsT = [PointF(xobsTvec[2i-1],xobsTvec[2i]) for i in subs]
n = length(subs)

JLD2.@save "data_exp-dgs.jld2" xobs0 xobsT n nshapes


text:
using DelimitedFiles

function calculate_cost(position, start_pos)
    return sum(1:abs(position-start_pos))
end

function test_calculate_cost()
    @assert calculate_cost(5, 16) == 66
    @assert calculate_cost(5, 1) == 10
end

test_calculate_cost()

function calculate_total_cost(position, position_vec)
    diff_vec = map(x->calculate_cost(position,x), position_vec)
    return sum(diff_vec)
end

function test_calculate_total_cost()
    test_list = [16,1,2,0,4,2,7,1,2,14]
    @assert calculate_total_cost(5, test_list) == 168
    @assert calculate_total_cost(2, test_list) == 206
end

test_calculate_total_cost()

# read input file, note single quotes around the delimiter
input = readdlm("input.txt", ',',Int)
# convert to a column vector
input = vec(input)

min_cost = typemax(Int)
min_cost_position = 0
for i in minimum(input):maximum(input)
    cost = calculate_total_cost(i, input)
    if cost < min_cost
        global min_cost = cost
        global min_cost_position = i
    end
end

println("Best position ", min_cost_position, " cost ", min_cost)


text:
using  Arborist
using Arborist: determine_type_params
using InteractiveUtils
using Test

@testset "determine_type_params" begin
    boop(x::Vector{T}) where {T} = string(T)
    meth = @which boop([1,2,3])
    @test [Int] == collect(determine_type_params(meth, Tuple{Vector{Int}}))

    bing(x::T, y::S) where {S,T} = x+y
    meth = @which bing(1, 1f0)
    @test [Float32, Int64] == collect(determine_type_params(meth, Tuple{Int, Float32}))
end


text:
# Invertible Glow network from Kingma et. al (2018)
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August 2020

export NetworkGlow1D, clear_grad!, get_params

"""
    H = NetworkGlow1D(nx, ny, n_in, batchsize, n_hidden, depth; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)
 Create a HINT network for data-driven generative modeling based
 on the change of variables formula.
 *Input*:
 - `nx`, `ny`, `n_in`, `batchsize`: spatial dimensions, number of channels and batchsize of input tensors `X`
 - `n_hidden`: number of hidden units in residual blocks
 - `depth`: number network layers
 - `k1`, `k2`: kernel size for first and third residual layer (`k1`) and second layer (`k2`)
 - `p1`, `p2`: respective padding sizes for residual block layers
 - `s1`, `s2`: respective strides for residual block layers
 *Output*:
 - `H`: HINT network
 *Usage:*
 - Forward mode: `Z, logdet = H.forward(X)`
 - Inverse mode: `X = H.inverse(Z)`
 - Backward mode: `ΔX, X = H.backward(ΔZ, Z)`
 *Trainable parameters:*
 - None in `H` itself
 - Trainable parameters in activation normalizations `H.AN[i]`, and in coupling layers `H.CL[i]`, where `i` ranges from `1` to `depth`.
 See also: [`ActNorm`](@ref), [`CouplingLayerHINT`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct NetworkGlow1D <: InvertibleNetwork
    AN::AbstractArray{ActNorm, 1}
    CL::AbstractArray{CouplingLayerGlow, 1}
    logdet::Bool
    is_reversed::Bool
end

@Flux.functor NetworkGlow1D

# Constructor
function NetworkGlow1D(n_in::Int64, batchsize::Int64, n_hidden::Int64, depth::Int64; logdet=true)

    # AN = Array{ActNorm}(undef, depth)
    AN = Array{ActNorm}(undef, depth+1)
    CL = Array{CouplingLayerGlow}(undef, depth)

    # Create layers
    for j=1:depth
        AN[j] = ActNorm(n_in; logdet=logdet)
        CL[j] = CouplingLayerGlow(1, 1, n_in, n_hidden, batchsize; k1=1, k2=1, p1=0, p2=0, logdet=logdet)
    end
    AN[depth+1] = ActNorm(n_in; logdet=logdet)

    return NetworkGlow1D(AN, CL, logdet, false)
end

# Forward pass and compute logdet
function forward(X, H::NetworkGlow1D; logdet=nothing)
    isnothing(logdet) ? logdet = (H.logdet && ~H.is_reversed) : logdet = logdet

    depth = length(H.CL)
    logdet_ = 0f0
    for j=1:depth
        logdet ? (X, logdet_an) = H.AN[j].forward(X) : X = H.AN[j].forward(X)
        logdet ? (X, logdet_cl) = H.CL[j].forward(X) : X = H.CL[j].forward(X)
        logdet && (logdet_ += (logdet_an + logdet_cl))
    end
    logdet ? (X, logdet_an) = H.AN[depth+1].forward(X) : X = H.AN[depth+1].forward(X)
    logdet && (logdet_ += logdet_an)
    logdet ? (return X, logdet_) : (return X)
end

# Inverse pass and compute gradients
function inverse(Z, H::NetworkGlow1D; logdet=nothing)
    isnothing(logdet) ? logdet = (H.logdet && H.is_reversed) : logdet = logdet

    depth = length(H.CL)
    logdet_ = 0f0
    logdet ? (Z, logdet_an) = H.AN[depth+1].inverse(Z; logdet=true) : Z = H.AN[depth+1].inverse(Z; logdet=false)
    logdet && (logdet_ += logdet_an)
    for j=depth:-1:1
        logdet ? (Z, logdet_an) = H.CL[j].inverse(Z) : Z = H.CL[j].inverse(Z)
        logdet ? (Z, logdet_cl) = H.AN[j].inverse(Z; logdet=true) : Z = H.AN[j].inverse(Z; logdet=false)
        logdet && (logdet_ += (logdet_an + logdet_cl))
    end
    logdet ? (return Z, logdet_) : (return Z)
end

# Backward pass and compute gradients
function backward(ΔZ, Z, H::NetworkGlow1D)
    depth = length(H.CL)
    ΔZ, Z = H.AN[depth+1].backward(ΔZ, Z)
    for j=depth:-1:1
        ΔZ, Z = H.CL[j].backward(ΔZ, Z)
        ΔZ, Z = H.AN[j].backward(ΔZ, Z)
    end
    return ΔZ, Z
end

# Clear gradients
function clear_grad!(H::NetworkGlow1D)
    depth = length(H.CL)
    for j=1:depth
        clear_grad!(H.AN[j])
        clear_grad!(H.CL[j])
    end
    clear_grad!(H.AN[depth+1])
end

# Get parameters
function get_params(H::NetworkGlow1D)
    depth = length(H.CL)
    p = []
    for j=1:depth
        p = cat(p, get_params(H.AN[j]); dims=1)
        p = cat(p, get_params(H.CL[j]); dims=1)
    end
    p = cat(p, get_params(H.AN[depth+1]); dims=1)
    return p
end

# Set is_reversed flag in full network tree
function tag_as_reversed!(H::NetworkGlow1D, tag::Bool)
    depth = length(H.CL)
    H.is_reversed = tag
    for j=1:depth
        tag_as_reversed!(H.AN[j], tag)
        tag_as_reversed!(H.CL[j], tag)
    end
    tag_as_reversed!(H.AN[depth+1], tag)
    return H
end


text:
using OrdinaryDiffEq, DiffEqProblemLibrary
prob = prob_ode_2Dlinear

integrator = init(prob,Euler();dt=1//2^(4))

sol =solve(prob,Euler();dt=1//2^(4),maxiters=1e6)

integrator = init(prob,Tsit5();dt=1//2^(4))

sol =solve(prob,Tsit5();dt=1//2^(4))

# plot(sol,plot_analytic=true)

sol =solve(prob,ExplicitRK();dt=1//2^(4))


text:
## --------------------------------------------------------------------
# Load exchange NCI60 data
const NCI60_exchanges = Dict()
function _load_NCI60_exchanges_data()
    empty!(NCI60_exchanges)
    for model_type in [EC_MODEL_TYPE_KEY, MODEL_TYPE_KEY]
        dict = get!(NCI60_exchanges, model_type, Dict())
        for const_level in CONSTRAINTS_LEVELS
            filename = string(model_type, "s_const_", const_level, "_exchangeFluxesComp.txt")
            path = joinpath(HUMAN1_NCI60_RESULTS_DIR, filename)
            !isfile(path) && continue
            dict[const_level] = CSV.read(path, DataFrame)
        end
    end
end

## --------------------------------------------------------------------
# INTERFACE
get_NCI60_msd_exchanges(model_type, const_level) = String.(NCI60_exchanges[model_type][const_level][:, :exchangeIDs])
get_NCI60_msd_exchanges_mets(model_type, const_level) = String.(NCI60_exchanges[model_type][const_level][:, :exchangeMets])

function get_NCI60_exchange(model_type, cell_line, const_level, data_type)
    df = NCI60_exchanges[model_type][const_level]
    k = data_type == :exp ? Symbol(string("exp_", cell_line)) : Symbol(cell_line)
    Float64.(df[:, k])
end

get_NCI60_exchange(model_type, cell_line, const_level, data_type, idx::Int) =
    get_NCI60_exchange(model_type, cell_line, const_level, data_type)[idx]

function get_NCI60_exchange(model_type, cell_line, const_level, data_type, ider::String)
    exchs = get_NCI60_msd_exchanges(model_type, const_level)
    mets = get_NCI60_msd_exchanges_mets(model_type, const_level)
    idx = something(findfirst(isequal(ider), exchs), findfirst(isequal(ider), mets), -1)
    idx == -1 && error("'$ider' not found")
    get_NCI60_exchange(model_type, cell_line, const_level, data_type, idx)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayConnectionDraining <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    drainTimeoutInSec::Any # spec type: Union{ Nothing, Int32 } # spec name: drainTimeoutInSec

    function ApplicationGatewayConnectionDraining(;enabled=nothing, drainTimeoutInSec=nothing)
        o = new()
        validate_property(ApplicationGatewayConnectionDraining, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(ApplicationGatewayConnectionDraining, Symbol("drainTimeoutInSec"), drainTimeoutInSec)
        setfield!(o, Symbol("drainTimeoutInSec"), drainTimeoutInSec)
        o
    end
end # type ApplicationGatewayConnectionDraining

const _property_map_ApplicationGatewayConnectionDraining = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("drainTimeoutInSec")=>Symbol("drainTimeoutInSec"))
const _property_types_ApplicationGatewayConnectionDraining = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("drainTimeoutInSec")=>"Int32")
Base.propertynames(::Type{ ApplicationGatewayConnectionDraining }) = collect(keys(_property_map_ApplicationGatewayConnectionDraining))
Swagger.property_type(::Type{ ApplicationGatewayConnectionDraining }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayConnectionDraining[name]))}
Swagger.field_name(::Type{ ApplicationGatewayConnectionDraining }, property_name::Symbol) =  _property_map_ApplicationGatewayConnectionDraining[property_name]

function check_required(o::ApplicationGatewayConnectionDraining)
    (getproperty(o, Symbol("enabled")) === nothing) && (return false)
    (getproperty(o, Symbol("drainTimeoutInSec")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ApplicationGatewayConnectionDraining }, name::Symbol, val)
    if name === Symbol("drainTimeoutInSec")
        Swagger.validate_param(name, "ApplicationGatewayConnectionDraining", :maximum, val, 3600, false)
        Swagger.validate_param(name, "ApplicationGatewayConnectionDraining", :minimum, val, 1, false)
    end
end


text:
using Test

using PyPlot
using LaTeXStrings

using Oceananigans.Advection

# Define a few utilities for running tests and unpacking and plotting results
include("ConvergenceTests/ConvergenceTests.jl")

using .ConvergenceTests
using .ConvergenceTests.OneDimensionalGaussianAdvectionDiffusion: run_test
using .ConvergenceTests.OneDimensionalUtils: unpack_errors, defaultcolors, removespines

""" Run advection test for all Nx in resolutions. """
function run_convergence_test(κ, U, resolutions, advection_scheme)

    # Determine safe time-step
           Lx = 2.5
            h = Lx / maximum(resolutions)
           Δt = min(0.01 * h / U, 0.1 * h^2 / κ)
    stop_time = Δt

    # Run the tests
    results = [run_test(Nx=Nx, Δt=Δt, advection=advection_scheme, stop_iteration=1,
                        U=U, κ=κ) for Nx in resolutions]

    return results
end

#####
##### Run test
#####

advection_schemes = (CenteredSecondOrder(), CenteredFourthOrder(), UpwindBiasedThirdOrder(),
                     UpwindBiasedFifthOrder(), WENO5())

U = 1
κ = 1e-8
Nx = [8, 16, 32, 64, 96, 128, 192, 256, 384, 512]

tolerance(::CenteredSecondOrder)    = 0.05
tolerance(::CenteredFourthOrder)    = 0.05
tolerance(::UpwindBiasedThirdOrder) = 0.30
tolerance(::UpwindBiasedFifthOrder) = 0.30
tolerance(::WENO5)                  = 0.40

test_resolution(::CenteredSecondOrder)    = 512
test_resolution(::CenteredFourthOrder)    = 512
test_resolution(::UpwindBiasedThirdOrder) = 128
test_resolution(::UpwindBiasedFifthOrder) = 128
test_resolution(::WENO5)                  = 512

rate_of_convergence(::CenteredSecondOrder) = 2
rate_of_convergence(::CenteredFourthOrder) = 4
rate_of_convergence(::UpwindBiasedThirdOrder) = 3
rate_of_convergence(::UpwindBiasedFifthOrder) = 5
rate_of_convergence(::WENO5) = 5
rate_of_convergence(::WENO{K}) where K = 2K-1

results = Dict()
for scheme in advection_schemes
    t_scheme = typeof(scheme)
    results[t_scheme] = run_convergence_test(κ, U, Nx, scheme)
end

colors = ("xkcd:royal blue", "xkcd:light red")

@testset "tmp" begin
    for scheme in advection_schemes

        fig, ax = subplots()

        t_scheme = typeof(scheme)
        name = string(t_scheme)
        roc = rate_of_convergence(scheme)
        atol = tolerance(scheme)
        atol = tolerance(scheme)
        Ntest = test_resolution(scheme)
        itest = searchsortedfirst(Nx, Ntest)

        (cx_L₁, cy_L₁, cz_L₁, uy_L₁, uz_L₁, vx_L₁, vz_L₁, wx_L₁, wy_L₁,
         cx_L∞, cy_L∞, cz_L∞, uy_L∞, uz_L∞, vx_L∞, vz_L∞, wx_L∞, wy_L∞) = unpack_errors(results[typeof(scheme)])

        test_rate_of_convergence(cx_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cx_L₁")
        test_rate_of_convergence(cy_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cy_L₁")
        test_rate_of_convergence(cz_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cz_L₁")

        test_rate_of_convergence(uy_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uy_L₁")
        test_rate_of_convergence(uz_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uz_L₁")

        test_rate_of_convergence(vx_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vx_L₁")
        test_rate_of_convergence(vz_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vz_L₁")

        test_rate_of_convergence(wx_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wx_L₁")
        test_rate_of_convergence(wy_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wy_L₁")

        #=
        test_rate_of_convergence(cx_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cx_L∞")
        test_rate_of_convergence(cy_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cy_L∞")
        test_rate_of_convergence(cz_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cz_L∞")

        test_rate_of_convergence(uy_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uy_L∞")
        test_rate_of_convergence(uz_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uz_L∞")

        test_rate_of_convergence(vx_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vx_L∞")
        test_rate_of_convergence(vz_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vz_L∞")

        test_rate_of_convergence(wx_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wx_L∞")
        test_rate_of_convergence(wy_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wy_L∞")
        =#

        @test cx_L₁ ≈ cy_L₁
        @test cx_L₁ ≈ cz_L₁
        @test uy_L₁ ≈ uz_L₁
        @test vx_L₁ ≈ vz_L₁
        @test wx_L₁ ≈ wy_L₁

        @test cx_L∞ ≈ cy_L∞
        @test cx_L∞ ≈ cz_L∞
        @test uy_L∞ ≈ uz_L∞
        @test vx_L∞ ≈ vz_L∞
        @test wx_L∞ ≈ wy_L∞

        common_kwargs = (linestyle="None", color=colors[1], mfc="None", alpha=0.8)

        loglog(Nx, cx_L₁; basex=2, marker="*", label="\$L_1\$-norm, \$c(x)\$ $name", common_kwargs...)
        loglog(Nx, cy_L₁; basex=2, marker="+", label="\$L_1\$-norm, \$c(y)\$ $name", common_kwargs...)
        loglog(Nx, cz_L₁; basex=2, marker="_", label="\$L_1\$-norm, \$c(z)\$ $name", common_kwargs...)

        loglog(Nx, uy_L₁; basex=2, marker="1", label="\$L_1\$-norm, \$u(y)\$ $name", common_kwargs...)
        loglog(Nx, uz_L₁; basex=2, marker="^", label="\$L_1\$-norm, \$u(z)\$ $name", common_kwargs...)

        loglog(Nx, vx_L₁; basex=2, marker="s", label="\$L_1\$-norm, \$v(x)\$ $name", common_kwargs...)
        loglog(Nx, vz_L₁; basex=2, marker="v", label="\$L_1\$-norm, \$v(z)\$ $name", common_kwargs...)

        loglog(Nx, wx_L₁; basex=2, marker="X", label="\$L_1\$-norm, \$w(x)\$ $name", common_kwargs...)
        loglog(Nx, wy_L₁; basex=2, marker="D", label="\$L_1\$-norm, \$w(y)\$ $name", common_kwargs...)

        common_kwargs = (linestyle="None", color=colors[2], mfc="None", alpha=0.8)

        loglog(Nx, cx_L∞; basex=2, marker="*", label="\$L_\\infty\$-norm, \$c(x)\$ $name", common_kwargs...)
        loglog(Nx, cy_L∞; basex=2, marker="+", label="\$L_\\infty\$-norm, \$c(y)\$ $name", common_kwargs...)
        loglog(Nx, cz_L∞; basex=2, marker="_", label="\$L_\\infty\$-norm, \$c(z)\$ $name", common_kwargs...)

        loglog(Nx, uy_L∞; basex=2, marker="1", label="\$L_\\infty\$-norm, \$u(y)\$ $name", common_kwargs...)
        loglog(Nx, uz_L∞; basex=2, marker="^", label="\$L_\\infty\$-norm, \$u(z)\$ $name", common_kwargs...)

        loglog(Nx, vx_L∞; basex=2, marker="s", label="\$L_\\infty\$-norm, \$v(x)\$ $name", common_kwargs...)
        loglog(Nx, vz_L∞; basex=2, marker="v", label="\$L_\\infty\$-norm, \$v(z)\$ $name", common_kwargs...)

        loglog(Nx, wx_L∞; basex=2, marker="X", label="\$L_\\infty\$-norm, \$w(x)\$ $name", common_kwargs...)
        loglog(Nx, wy_L∞; basex=2, marker="D", label="\$L_\\infty\$-norm, \$w(y)\$ $name", common_kwargs...)

        label = raw"\sim N_x^{-" * "$roc" * raw"}" |> latexstring

        loglog(Nx[itest-3:itest], uy_L₁[itest] .* (Nx[itest] ./ Nx[itest-3:itest]) .^ roc, color=colors[1], basex=2, alpha=0.8, label=label)

        xlabel(L"N_x")
        ylabel("\$L\$-norms of \$ | c_\\mathrm{sim} - c_\\mathrm{analytical} |\$")
        removespines("top", "right")
        lgd = legend(loc="upper right", bbox_to_anchor=(1.4, 1.0), prop=Dict(:size=>6))

        filepath = joinpath(@__DIR__, "figs", "one_dimensional_convergence_$name.png")
        savefig(filepath, dpi=480, bbox_extra_artists=(lgd,), bbox_inches="tight")
        close(fig)
    end
end


text:
"""
    Atom(T <: ThreeD, ::Int64, ::Union{Real, Integer})
"""
function Atom(geometry::Sphere, N::Int64, kR::Union{Real, Integer}; createFunction=ftn_AtomsOnSphere::Function)
    @debug "start: Shape - Sphere"

    dimensions = 3
    ρ = 3N / (4π*kR^3)
    rₘᵢₙ = get_rₘᵢₙ(ρ)
    if rₘᵢₙ ≥ kR / 10
        rₘᵢₙ = kR / 100
    end

    r = get_atoms(dimensions, N, rₘᵢₙ; createFunction, kR)

    @debug "end  : Shape - Sphere"
    return Atom(Sphere(), r, N, Float64(kR))
end

function Atom(geometry::Sphere, r::Matrix, kR::Union{Real, Integer})
    N = size(r, 2) # remember to use each collum as a atom position
    return Atom(Sphere(), r, N, Float64(kR))
end

function ftn_AtomsOnSphere(; kwargs...)
    kR = kwargs[:kR]
    new_atom = zeros(3)
    new_atom[1] = 2π * rand() # azimuth
    new_atom[2] = asin(2 * rand() - 1) # elevation
    new_atom[3] = kR * (rand()^(1 ./ 3.0)) # radii
    new_atom[:] = sph2cart(new_atom)
    return new_atom
end
"""
    spherical_coordinate=[azimuth, elevation, r]

    azimuth = projection on XY-plane (in radians)
    "elevation" or "Polar" = projection on Z-axis (in radians)
    r = radius

	ref: https://www.mathworks.com/help/matlab/ref/sph2cart.html
"""
function sph2cart(spherical_coordinate)
    azimuth = spherical_coordinate[1]
    elevation = spherical_coordinate[2]
    radius = spherical_coordinate[3]
    x = radius * cos(elevation) * cos(azimuth)
    y = radius * cos(elevation) * sin(azimuth)
    z = radius * sin(elevation)
    return [x, y, z]
end