text:
export unity, fast_unity,
       unity_material_family

# We need some additional Encoders
encode_GameObject = encode_int
decode_GameObject = decode_int_or_error
encode_GameObject_array = encode_int_array
decode_GameObject_array = decode_int_array
encode_Material = encode_int
decode_Material = decode_int_or_error

# Must convert from local to world coordinates
encode_Vector3(c::IO, v::Union{XYZ, VXYZ}) = begin
  v = in_world(v)
#  encode_float3(c, v.x, v.y, v.z)
  encode_float3(c, v.x, v.z, v.y)
end
decode_Vector3(c::IO) =
  let x = decode_float(c)
      z = decode_float(c)
      y = decode_float(c)
    xyz(x, y, z, world_cs)
  end
encode_Vector3_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_Vector3(c, e) end
end
decode_Vector3_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{XYZ}(undef, len)
  for i in 1:len
    r[i] = decode_Vector3(c)
  end
  r
end

encode_Vector3_array_array(c::IO, v::Vector) = begin
  encode_int(c, length(v))
  for e in v encode_Vector3_array(c, e) end
end
decode_Vector3_array_array(c::IO) = begin
  len = decode_int_or_error(c)
  r = Vector{Vector{XYZ}}(undef, len)
  for i in 1:len
    r[i] = decode_Vector3_array(c)
  end
  r
end

encode_Color(c::IO, v::RGB) = begin
  encode_float3(c, v.r, v.g, v.b)
end
decode_Color(c::IO) =
  let r = decode_float(c)
      g = decode_float(c)
      b = decode_float(c)
    rgb(r, g, b)
  end

#=
encode_Quaternion(c::IO, pv::Union{XYZ, VXYZ}) =
  let tr = pv.cs.transform
    for i in 1:3
      for j in 1:3
        encode_float(c, tr[i,j])
      end
    end
  end
=#

unity_api = @remote_functions :CS """
public void SetApplyMaterials(bool apply)
public void SetApplyColliders(bool apply)
public GameObject SurfacePolygon(Vector3[] ps)
public GameObject SurfacePolygonNamed(String name, Vector3[] ps, Material material)
public GameObject SurfacePolygonWithMaterial(Vector3[] ps, Material material)
public GameObject SurfaceMeshNamed(String name, Vector3[] vertices, int[] triangles, Material material)
public GameObject SurfaceMeshWithMaterial(Vector3[] vertices, int[] triangles, Material material)
public GameObject SurfaceMesh(Vector3[] vertices, int[] triangles)
public GameObject Text(string txt, Vector3 position, Vector3 vx, Vector3 vy, string fontName, int fontSize)
public GameObject Sphere(Vector3 center, float radius)
public GameObject SphereWithMaterial(Vector3 center, float radius, Material material)
public GameObject Pyramid(Vector3[] ps, Vector3 q)
public GameObject PyramidFrustum(Vector3[] ps, Vector3[] qs)
public GameObject PyramidFrustumWithMaterial(Vector3[] ps, Vector3[] qs, Material material)
public GameObject ExtrudeContour(Vector3[] contour, Vector3[][] holes, Vector3 v, Material material)
public GameObject RightCuboid(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz, float angle)
public GameObject RightCuboidWithMaterial(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz, float angle, Material material)
public GameObject Box(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz)
public GameObject Cylinder(Vector3 bottom, float radius, Vector3 top)
public GameObject CylinderWithMaterial(Vector3 bottom, float radius, Vector3 top, Material material)
public GameObject Unite(GameObject s0, GameObject s1)
public GameObject Intersect(GameObject s0, GameObject s1)
public GameObject Subtract(GameObject s0, GameObject s1)
public void SubtractFrom(GameObject s0, GameObject s1)
public GameObject Canonicalize(GameObject s)
public void Move(GameObject s, Vector3 v)
public void Scale(GameObject s, Vector3 p, float scale)
public void Rotate(GameObject s, Vector3 p, Vector3 n, float a)
public GameObject SurfaceFromGrid(int m, int n, Vector3[] pts, bool closedM, bool closedN, int level)
public GameObject LoadResource(String name)
public Material LoadMaterial(String name)
public Material CurrentMaterial()
public void SetCurrentMaterial(Material material)
public GameObject InstantiateResource(GameObject family, Vector3 pos, Vector3 vx, Vector3 vy, float scale)
public GameObject InstantiateBIMElement(GameObject family, Vector3 pos, float angle)
public GameObject Window(Vector3 position, Quaternion rotation, float dx, float dy, float dz)
public GameObject Shelf(Vector3 position, int rowLength, int lineLength, float cellWidth, float cellHeight, float cellDepth)
public GameObject Slab(Vector3[] contour, Vector3[][] holes, float h, Material material)
public GameObject BeamRectSection(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz, float angle, Material material)
public GameObject BeamCircSection(Vector3 bot, float radius, Vector3 top, Material material)
public GameObject Panel(Vector3[] pts, Vector3 n, Material material)
public void SetView(Vector3 position, Vector3 target, float lens)
public Vector3 ViewCamera()
public Vector3 ViewTarget()
public float ViewLens()
public void DeleteAll()
public void DeleteMany(GameObject[] objs)
public GameObject CreateParent(String name, bool active)
public GameObject CurrentParent()
public GameObject SetCurrentParent(GameObject newParent)
public void SetActive(GameObject obj, bool state)
public void DeleteAllInParent(GameObject parent)
public void SwitchToParent(GameObject newParent)
public int SetMaxNonInteractiveRequests(int n)
public void SetNonInteractiveRequests()
public void SetInteractiveRequests()
public GameObject CreateBlockInstance(GameObject block, Vector3 position, Vector3 vx, Vector3 vy, float scale)
public GameObject CreateBlockFromShapes(String name, GameObject[] objs)
public GameObject PointLight(Vector3 position, Color color, float range, float intensity)
public Point3d[] GetPosition(string prompt)
public ObjectId[] GetPoint(string prompt)
public ObjectId[] GetCurve(string prompt)
public ObjectId[] GetSurface(string prompt)
public ObjectId[] GetSolid(string prompt)
public ObjectId[] GetShape(string prompt)
public long GetHandleFromShape(Entity e)
public ObjectId GetShapeFromHandle(long h)
public void RegisterForChanges(ObjectId id)
public void UnregisterForChanges(ObjectId id)
public ObjectId[] ChangedShape()
public void DetectCancel()
public void UndetectCancel()
public bool WasCanceled()
public ObjectId[] GetAllShapes()
public ObjectId[] GetAllShapesInLayer(ObjectId layerId)
public void SetResolution(int width, int height)
public void ScreenShot(String path)
public void SelectGameObjects(GameObject[] objs)
public void StartSelectingGameObject()
public void StartSelectingGameObjects()
public bool EndedSelectingGameObjects()
public int[] SelectedGameObjectsIds(bool existing)
public void SetSun(float altitude, float azimuth)
public string GetRenderResolution()
public float GetCurrentFPS()
public int GetViewTriangleCount()
public int GetViewVertexCount()
public String ShapeType(GameObject s)
public Vector3 SphereCenter(GameObject s)
public float SphereRadius(GameObject s)
"""

abstract type UnityKey end
const UnityId = Int
const UnityIds = Vector{UnityId}
const UnityRef = GenericRef{UnityKey, UnityId}
const UnityRefs = Vector{UnityRef}
const UnityEmptyRef = EmptyRef{UnityKey, UnityId}
const UnityUniversalRef = UniversalRef{UnityKey, UnityId}
const UnityNativeRef = NativeRef{UnityKey, UnityId}
const UnityUnionRef = UnionRef{UnityKey, UnityId}
const UnitySubtractionRef = SubtractionRef{UnityKey, UnityId}
const Unity = SocketBackend{UnityKey, UnityId}

void_ref(b::Unity) = UnityNativeRef(-1)

create_Unity_connection() =
    begin
        #check_plugin()
        create_backend_connection("Unity", 11002)
    end

const unity = Unity(LazyParameter(TCPSocket, create_Unity_connection), unity_api)

backend_name(b::Unity) = "Unity"

(backend::Unity)(; apply_materials=true, apply_colliders=true) =
  begin
    @remote(backend, SetApplyMaterials(apply_materials))
    @remote(backend, SetApplyColliders(apply_colliders))
    backend
  end

realize(b::Unity, s::EmptyShape) =
  UnityEmptyRef()
realize(b::Unity, s::UniversalShape) =
  UnityUniversalRef()

fast_unity() =
  begin
    @remote(unity, SetApplyMaterials(false))
    @remote(unity, SetApplyColliders(false))
  end

slow_unity() =
  begin
    @remote(unity, SetApplyMaterials(true))
    @remote(unity, SetApplyColliders(true))
  end
#=
realize(b::Unity, s::Point) =
  @remote(b, Point(s.position))
realize(b::Unity, s::Line) =
  @remote(b, PolyLine(s.vertices))
realize(b::Unity, s::Spline) =
  if (s.v0 == false) && (s.v1 == false)
    #@remote(b, Spline(s.points))
    @remote(b, InterpSpline(connection(b),
                     s.points,
                     s.points[2]-s.points[1],
                     s.points[end]-s.points[end-1]))
  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::Unity, s::ClosedSpline) =
  @remote(b, InterpClosedSpline(s.points))
realize(b::Unity, s::Circle) =
  @remote(b, Circle(s.center, vz(1, s.center.cs), s.radius))
realize(b::Unity, 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::Unity, s::Ellipse) =
  if s.radius_x > s.radius_y
    @remote(b, Ellipse(s.center, vz(1, s.center.cs), vxyz(s.radius_x, 0, 0, s.center.cs), s.radius_y/s.radius_x))
  else
    @remote(b, Ellipse(s.center, vz(1, s.center.cs), vxyz(0, s.radius_y, 0, s.center.cs), s.radius_x/s.radius_y))
  end
realize(b::Unity, s::EllipticArc) =
  error("Finish this")

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

#realize(b::Unity, s::SurfaceElliptic_Arc) = @remote(b, Circle(connection(b),
#realize(b::Unity, s::SurfaceEllipse) = @remote(b, Circle(connection(b),
=#

backend_surface_polygon(b::Unity, vs::Locs) =
  @remote(b, SurfacePolygon(vs))

backend_fill(b::Unity, path::ClosedPolygonalPath) =
  @remote(b, SurfacePolygon(path.vertices))

#=
realize(b::Unity, s::Surface) =
  let #ids = map(r->@remote(b, NurbSurfaceFrom(connection(b),r), @remote(b, SurfaceFromCurves(collect_ref(s.frontier))))
      ids = @remote(b, SurfaceFromCurves(collect_ref(s.frontier)))
    foreach(mark_deleted, s.frontier)
    ids
  end
backend_surface_boundary(b::Unity, s::Shape2D) =
    map(shape_from_ref, @remote(b, CurvesFromSurface(ref(s).value)))
=#
backend_fill(b::Unity, path::ClosedPathSequence) =
  backend_fill(b, convert(ClosedPolygonalPath, path))

#=
# Iterating over curves and surfaces

Unity"public double[] CurveDomain(Entity ent)"
Unity"public double CurveLength(Entity ent)"
Unity"public Frame3d CurveFrameAt(Entity ent, double t)"
Unity"public Frame3d CurveFrameAtLength(Entity ent, double l)"
=#

backend_map_division(b::Unity, f::Function, s::Shape1D, n::Int) =
  let (t1, t2) = curve_domain(s)
    map_division(t1, t2, n) do t
      f(frame_at(s, t))
    end
  end
#=
Unity"public Vector3d RegionNormal(Entity ent)"
Unity"public Point3d RegionCentroid(Entity ent)"
Unity"public double[] SurfaceDomain(Entity ent)"
Unity"public Frame3d SurfaceFrameAt(Entity ent, double u, double v)"

backend_surface_domain(b::Unity, s::Shape2D) =
    tuple(@remote(b, SurfaceDomain(ref(s).value)...))

backend_map_division(b::Unity, 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

# The previous method cannot be applied to meshes in AutoCAD, which are created by surface_grid

backend_map_division(b::Unity, f::Function, s::SurfaceGrid, nu::Int, nv::Int) =
    let (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::Unity, s::Text) =
  @remote(b, Text(s.str, s.corner, vz(-1, s.corner.cs), vy(1, s.corner.cs), "Fonts/Inconsolata-Regular", s.height))

realize(b::Unity, s::Sphere) =
  @remote(b, Sphere(s.center, s.radius))

#=
realize(b::Unity, s::Torus) =
  @remote(b, Torus(s.center, vz(1, s.center.cs), s.re, s.ri))
=#

realize(b::Unity, s::Cuboid) =
  @remote(b, PyramidFrustum([s.b0, s.b1, s.b2, s.b3], [s.t0, s.t1, s.t2, s.t3]))

realize(b::Unity, s::RegularPyramidFrustum) =
    @remote(b, PyramidFrustum(
                    regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
                    regular_polygon_vertices(s.edges, add_z(s.cb, s.h), s.rt, s.angle, s.inscribed)))

realize(b::Unity, s::RegularPyramid) =
  @remote(b, Pyramid(
               regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
               add_z(s.cb, s.h)))

realize(b::Unity, s::IrregularPyramid) =
  @remote(b, Pyramid(s.bs, s.t))

realize(b::Unity, s::RegularPrism) =
  let bs = regular_polygon_vertices(s.edges, s.cb, s.r, s.angle, s.inscribed)
    @remote(b, PyramidFrustum(bs, map(p -> add_z(p, s.h), bs)))
  end

realize(b::Unity, s::IrregularPyramidFrustum) =
    @remote(b, PyramidFrustum(s.bs, s.ts))

realize(b::Unity, s::IrregularPrism) =
  @remote(b, PyramidFrustum(s.bs, map(p -> (p + s.v), s.bs)))

backend_right_cuboid(b::Unity, cb, width, height, h, angle, material) =
  isnothing(material) ?
    @remote(b, RightCuboid(cb, vz(1, cb.cs), vx(1, cb.cs), height, width, h, angle)) :
    @remote(b, RightCuboidWithMaterial(cb, vz(1, cb.cs), vx(1, cb.cs), height, width, h, angle, material))

#unity"public GameObject Box2(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz)"

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

realize(b::Unity, s::Cone) =
  @remote(b, Pyramid(regular_polygon_vertices(64, s.cb, s.r), add_z(s.cb, s.h)))

realize(b::Unity, s::ConeFrustum) =
  @remote(b, PyramidFrustum(
    regular_polygon_vertices(64, s.cb, s.rb),
    regular_polygon_vertices(64, s.cb + vz(s.h, s.cb.cs), s.rt)))

backend_cylinder(b::Unity, cb::Loc, r::Real, h::Real) =
  @remote(b, Cylinder(cb, r, add_z(cb, h)))

backend_cylinder(b::Unity, cb::Loc, r::Real, h::Real, material) =
  @remote(b, CylinderWithMaterial(cb, r, add_z(cb, h), material))

#=
backend_extrusion(b::Unity, s::Shape, v::Vec) =
    and_mark_deleted(
        map_ref(s) do r
            @remote(b, Extrude(r, v))
        end,
        s)

backend_sweep(b::Unity, path::Shape, profile::Shape, rotation::Real, scale::Real) =
  map_ref(profile) do profile_r
    map_ref(path) do path_r
      @remote(b, Sweep(path_r, profile_r, rotation, scale))
    end
  end

realize(b::Unity, 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_curves(b::Unity, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  and_delete_shapes(UnityLoft(connection(b),
                             collect_ref(profiles),
                             collect_ref(rails),
                             ruled, closed),
                    vcat(profiles, rails))

            MAYBE USE THIS
            ruled_surface(s1, s2) =
                let pts1 = map_division(in_world, s1, 20),
                    pts2 = map_division(in_world, s2, 20)
                  iterate_quads((p0, p1, p2, p3)->(surface_polygon([p0,p1,p3]), surface_polygon([p1,p2,p3])),
                                [pts1, pts2])
                end

            ruled_surface(s1, s2)

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

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

backend_loft_surface_point(b::Unity, profile::Shape, point::Shape) =
    backend_loft_curve_point(b, profile, point)

=#

###
unite_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    ensure_ref(b, @remote(b, Unite(r0.value, r1.value)))

intersect_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    ensure_ref(b, @remote(b, Intersect(r0.value, r1.value)))

subtract_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    let r = @remote(b, Subtract(r0.value, r1.value))
      @remote(b, DeleteMany([r0.value, r1.value]))
      r
    end

#=
subtract_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    begin
      @remote(b, SubtractFrom(r0.value, r1.value))
      r0.value
    end
=#

#=
slice_ref(b::Unity, r::UnityNativeRef, p::Loc, v::Vec) =
    (@remote(b, Slice(r.value, p, v); r))

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

=#
unite_refs(b::Unity, refs::Vector{<:UnityRef}) =
    UnityUnionRef(tuple(refs...))

#
realize(b::Unity, s::UnionShape) =
  let r = foldl((r0,r1)->unite_ref(b,r0,r1), map(ref, s.shapes),
                init=UnityEmptyRef())
    delete_shapes(s.shapes)
    #@remote(b, Canonicalize(r.value))
    r
  end

realize(b::Unity, s::IntersectionShape) =
  let r = foldl((r0,r1)->intersect_ref(b,r0,r1), map(ref, s.shapes),
                init=UnityUniversalRef())
    delete_shapes(s.shapes)
    r
  end

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

realize(b::Unity, 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::Unity, s::Transform) =
  let r = map_ref(s.shape) do r
            @remote(b, Transform(r, s.xform))
            r
          end
    mark_deleted(s.shape)
    r
  end
=#
realize(b::Unity, 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::Unity, 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::Unity, 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::Unity, s::UnionMirror) =
  let r0 = ref(s.shape),
      r1 = map_ref(s.shape) do r
            @remote(b, Mirror(r, s.p, s.n, true))
          end
    UnionRef((r0,r1))
  end
=#

backend_surface_grid(b::Unity, points, closed_u, closed_v, smooth_u, smooth_v) =
  # we create two surfaces to have normals on both sides
  let ptss = points,
      s1 = size(ptss,1),
      s2 = size(ptss,2),
      refs = UnityId[]
    if smooth_u && smooth_v
      push!(refs, @remote(b, SurfaceFromGrid(s2, s1, reshape(ptss,:), closed_u, closed_v, 2)))
    elseif smooth_u
      for i in 1:(closed_v ? s1 : s1-1)
        push!(refs, @remote(b, SurfaceFromGrid(s2, 2, reshape(ptss[[i,i%s1+1],:],:), closed_u, false, 2)))
      end
    elseif smooth_v
      for i in 1:(closed_u ? s2 : s2-1)
        push!(refs, @remote(b, SurfaceFromGrid(2, s1, reshape(ptss[:,[i,i%s1+1]],:), false, closed_v, 2)))
      end
    else
      for i in 1:(closed_v ? s1 : s1-1)
        for j in 1:(closed_u ? s2 : s2-1)
          push!(refs, @remote(b, SurfaceFromGrid(2, 2, reshape(ptss[[i,i%s1+1],[j,j%s2+1]],:), false, false, 2)))
        end
      end
    end
    refs
  end

backend_surface_mesh(b::Unity, vertices, faces) =
  @remote(b, SurfaceMesh(vertices, vcat(faces...)))

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

# backend_frame_at
backend_frame_at(b::Unity, s::Circle, t::Real) = add_pol(s.center, s.radius, t)

backend_frame_at(b::Unity, c::Shape1D, t::Real) = @remote(b, CurveFrameAt(ref(c).value, t))

#backend_frame_at(b::Unity, s::Surface, u::Real, v::Real) =
    #What should we do with v?
#    backend_frame_at(b, s.frontier[1], u)

#backend_frame_at(b::Unity, s::SurfacePolygon, u::Real, v::Real) =

backend_frame_at(b::Unity, s::Shape2D, u::Real, v::Real) = @remote(b, SurfaceFrameAt(ref(s).value, u, v))

=#

# BIM

# Families

abstract type UnityFamily <: Family end

struct UnityMaterialFamily <: UnityFamily
  name::String
end

unity_material_family(name, pairs...) = UnityMaterialFamily(name)
backend_get_family_ref(b::Unity, f::Family, uf::UnityMaterialFamily) = @remote(b, LoadMaterial(uf.name))

struct UnityResourceFamily <: UnityFamily
  name::String
end

unity_resource_family(name, pairs...) = UnityResourceFamily(name)
backend_get_family_ref(b::Unity, f::Family, uf::UnityResourceFamily) = @remote(b, LoadResource(uf.name))

set_backend_family(default_wall_family(), unity, unity_material_family("Default/Materials/Plaster"))
set_backend_family(default_slab_family(), unity, unity_material_family("Default/Materials/Concrete"))
set_backend_family(default_roof_family(), unity, unity_material_family("Default/Materials/Concrete"))
set_backend_family(default_beam_family(), unity, unity_material_family("Default/Materials/Aluminum"))
set_backend_family(default_column_family(), unity, unity_material_family("Default/Materials/Concrete"))
set_backend_family(default_door_family(), unity, unity_material_family("Default/Materials/Wood"))
set_backend_family(default_panel_family(), unity, unity_material_family("Default/Materials/Glass"))
set_backend_family(default_truss_node_family(), unity, unity_material_family("Default/Materials/Steel"))
set_backend_family(default_truss_bar_family(), unity, unity_material_family("Default/Materials/Steel"))

set_backend_family(default_table_family(), unity, unity_resource_family("Default/Prefabs/Table"))
set_backend_family(default_chair_family(), unity, unity_resource_family("Default/Prefabs/Chair"))
set_backend_family(default_table_chair_family(), unity, unity_resource_family("Default/Prefabs/TableChair"))

set_backend_family(default_curtain_wall_family().panel, unity, unity_material_family("Default/Materials/Glass"))
set_backend_family(default_curtain_wall_family().boundary_frame, unity, unity_material_family("Default/Materials/Steel"))
set_backend_family(default_curtain_wall_family().transom_frame, unity, unity_material_family("Default/Materials/Steel"))
set_backend_family(default_curtain_wall_family().mullion_frame, unity, unity_material_family("Default/Materials/Steel"))

backend_rectangular_table(b::Unity, c, angle, family) =
    @remote(b, InstantiateBIMElement(family_ref(b, family), c, -angle))

backend_chair(b::Unity, c, angle, family) =
    @remote(b, InstantiateBIMElement(family_ref(b, family), c, -angle))

backend_rectangular_table_and_chairs(b::Khepri.Unity, c, angle, family) =
    @remote(b, InstantiateBIMElement(family_ref(b, family), c, -angle))

backend_slab(b::Unity, profile, holes, thickness, family) =
  let bot_vs = path_vertices(profile)
    @remote(b, Slab(bot_vs, map(path_vertices, holes), thickness, family_ref(b, family)))
  end

## HACK: Is this still needed? backend_cylinder takes care of this, doesn't it?
realize_beam_profile(b::Unity, s::Union{Beam,FreeColumn,Column}, profile::CircularPath, cb::Loc, length::Real) =
  @remote(b, BeamCircSection(
    cb,
    profile.radius,
    add_z(cb, length*support_z_fighting_factor), #We reduce height just a bit to avoid Z-fighting
    family_ref(b, s.family)))

realize_beam_profile(b::Unity, s::Union{Beam,FreeColumn,Column}, profile::RectangularPath, cb::Loc, length::Real) =
  let profile_u0 = profile.corner,
      c = add_xy(cb, profile_u0.x + profile.dx/2, profile_u0.y + profile.dy/2)
    @remote(b, BeamRectSection(
      c, vz(1, c.cs), vx(1, c.cs), profile.dy, profile.dx,
      length*support_z_fighting_factor,
      -s.angle,
      family_ref(b, s.family)))
  end

#=
sweep_fractions(b, verts, height, l_thickness, r_thickness) =
  let p = add_z(verts[1], height/2),
      q = add_z(verts[2], height/2)
    if distance(p, q) > 1e-14
      let (c, h) = position_and_height(p, q),
          thickness = r_thickness + l_thickness,
          s = UnityNativeRef(@remote(b, RightCuboid(c, vz(1, c.cs), vx(1, c.cs), height, thickness, h, 0)))
        if length(verts) > 2
          (s, sweep_fractions(b, verts[2:end], height, l_thickness, r_thickness)...)
        else
          (s, )
        end
      end
    else
      sweep_fractions(b, verts[2:end], height, l_thickness, r_thickness)
    end
  end

backend_wall(b::Unity, path, height, l_thickness, r_thickness, family) =
  path_length(path) < path_tolerance() ?
    UnityEmptyRef() :
    begin
      @remote(b, SetCurrentMaterial(realize(b, family)))
      backend_wall_path(
          b,
          path,
          height*0.999, #We reduce height just a bit to avoid Z-fighting
          l_thickness, r_thickness)
    end

backend_wall_path(b::Unity, path::OpenPolygonalPath, height, l_thickness, r_thickness) =
    UnityUnionRef(sweep_fractions(b, path.vertices, height, l_thickness, r_thickness))

backend_wall_path(b::Unity, path::Path, height, l_thickness, r_thickness) =
    backend_wall_path(b, convert(OpenPolygonalPath, path), height, l_thickness, r_thickness)
=#

realize(b::Unity, w::Window) = void_ref(b)
realize(b::Unity, w::Door) = void_ref(b)

realize_pyramid_frustum(b::Unity, top_mat, bot_mat, side_mat, bot_vs::Locs, top_vs::Locs) =
  UnityNativeRef(@remote(b, PyramidFrustumWithMaterial(bot_vs, top_vs, top_mat)))

realize_prism(b::Unity, top, bot, side, path::PathSet, h::Real) =
  let v = planar_path_normal(path)*h,
      ptss = path_vertices.(path.paths)
    @remote(b, ExtrudeContour(ptss[1], ptss[2:end], v, top))
  end

#

realize(b::Unity, s::TrussNode) =
    @remote(b, SphereWithMaterial(s.p, s.family.radius, family_ref(b, s.family)))
realize(b::Unity, s::TrussBar) =
    @remote(b, CylinderWithMaterial(s.p0, s.family.radius, s.p1, family_ref(b, s.family)))

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

set_view(camera::Loc, target::Loc, lens::Real, aperture::Real, b::Unity) =
  let c = connection(b)
    @remote(b, SetView(camera, target, lens))
    interrupt_processing(c)
  end

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

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

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

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

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

set_length_unit(unit::String, b::Unity) = nothing # Unused, for now

#=
# Dimensions

const UnityDimensionStyles = Dict(:architectural => "_ARCHTICK", :mechanical => "")

dimension(p0::Loc, p1::Loc, p::Loc, scale::Real, style::Symbol, b::Unity=current_backend()) =
    @remote(b, CreateAlignedDimension(p0, p1, p,
        scale,
        UnityDimensionStyles[style]))

dimension(p0::Loc, p1::Loc, sep::Real, scale::Real, style::Symbol, b::Unity=current_backend()) =
    let v = p1 - p0
        angle = pol_phi(v)
        dimension(p0, p1, add_pol(p0, sep, angle + pi/2), scale, style, b)
    end

=#

# Layers
# Experiment for multiple, simultaneous, alternative layers
# Layers

UnityLayer = Int

current_layer(b::Unity)::UnityLayer =
  @remote(b, CurrentParent())

current_layer(layer::UnityLayer, b::Unity) =
  @remote(b, SetCurrentParent(layer))

backend_create_layer(b::Unity, name::String, active::Bool, color::RGB) =
  let layer = @remote(b, CreateParent(name, active))
    @warn "Ignoring color in create_layer for Unity"
    #@remote(b, SetLayerColor(layer, color.r, color.g, color.b))
    layer
  end

set_layer_active(layer::UnityLayer, status, b::Unity) =
  let c = connection(b)
    @remote(b, SetActive(layer, status))
    interrupt_processing(c)
  end

delete_all_shapes_in_layer(layer::UnityLayer, b::Unity) =
  @remote(b, DeleteAllInParent(layer))

switch_to_layer(layer::UnityLayer, b::Unity) =
  @remote(b, SwitchToParent(layer))

# To preserve interactiveness during background

preserving_interactiveness(f, b::Unity=current_backend()) =
  let prev = @remote(b, SetMaxNonInteractiveRequests(0))
    f()
    @remote(b, SetMaxNonInteractiveRequests(prev))
  end

# Experiment to speed up things

canonicalize_layer(layer::UnityLayer, b::Unity) =
  @remote(b, Canonicalize(layer))

# Materials

UnityMaterial = Int

current_material(b::Unity)::UnityMaterial =
  @remote(b, CurrentMaterial())

current_material(material::UnityMaterial, b::Unity) =
  @remote(b, SetCurrentMaterial(material))

get_material(name::String, b::Unity) =
  @remote(b, LoadMaterial(name))

# Blocks

realize(b::Unity, s::Block) =
  s.shapes == [] ?
    @remote(b, LoadResource(s.name)) :
    @remote(b, Canonicalize(@remote(b, CreateBlockFromShapes(s.name, collect_ref(s.shapes)))))

realize(b::Unity, s::BlockInstance) =
    @remote(b, CreateBlockInstance(
        ref(s.block).value,
        s.loc, vy(1, s.loc.cs), vz(1, s.loc.cs), s.scale))
#=

# Manual process
@time for i in 1:1000 for r in 1:10 circle(x(i*10), r) end end

# Create block...
Khepri.create_block("Foo", [circle(radius=r) for r in 1:10])

# ...and instantiate it
@time for i in 1:1000 Khepri.instantiate_block("Foo", x(i*10), 0) end

=#

# Lights

backend_pointlight(b::Unity, loc::Loc, color::RGB, range::Real, intensity::Real) =
    @remote(b, PointLight(loc, color, range, intensity))
#=
backend_spotlight(b::Unity, loc::Loc, dir::Vec, hotspot::Real, falloff::Real) =
    @remote(b, SpotLight(loc, hotspot, falloff, loc + dir))

backend_ieslight(b::Unity, file::String, loc::Loc, dir::Vec, alpha::Real, beta::Real, gamma::Real) =
    @remote(b, IESLight(file, loc, loc + dir, vxyz(alpha, beta, gamma)))

# User Selection
=#

shape_from_ref(r, b::Unity) =
  let idx = findfirst(s -> r in collect_ref(s), collected_shapes())
    if isnothing(idx)
      let kind = @remote(b, ShapeType(r))
        if kind == "Sphere"
          sphere(@remote(b, SphereCenter(r)), @remote(b, SphereRadius(r)),
                 backend=b, ref=LazyRef(b, UnityNativeRef(r)))
        else
          @warn "No shapes were previously collected (see in_shape_collection)"
          unknown(r, backend=b, ref=LazyRef(b, UnityNativeRef(r), 0, 0))
          #code = @remote(b, ShapeCode(r)),
          #ref = LazyRef(b, UnityNativeRef(r))
          #error("Unknown shape with code $(code)")
        end
      end
    else
      collected_shapes()[idx]
    end
  end
#
#=

select_position(prompt::String, b::Unity) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = @remote(b, GetPosition(prompt))
      length(ans) > 0 && ans[1]
    end
  end

select_with_prompt(prompt::String, b::Backend, f::Function) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = f(connection(b), prompt)
      length(ans) > 0 && shape_from_ref(ans[1], b)
    end
  end

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

select_point(prompt::String, b::Unity) =
  select_with_prompt(prompt, b, UnityGetPoint)

select_curve(prompt::String, b::Unity) =
  select_with_prompt(prompt, b, UnityGetCurve)

select_surface(prompt::String, b::Unity) =
  select_with_prompt(prompt, b, UnityGetSurface)

select_solid(prompt::String, b::Unity) =
  select_with_prompt(prompt, b, UnityGetSolid)

select_shape(prompt::String, b::Unity) =
  select_with_prompt(prompt, b, UnityGetShape)

captured_shape(b::Unity, handle) =
  shape_from_ref(@remote(b, GetShapeFromHandle(handle)),
                 b)

generate_captured_shape(s::Shape, b::Unity) =
    println("captured_shape(autocad, $(@remote(b, GetHandleFromShape(ref(s).value)))"))

# Register for notification

register_for_changes(s::Shape, b::Unity) =
    negin
        @remote(b, RegisterForChanges(ref(s).value))
        @remote(b, DetectCancel())
        s
    end

unregister_for_changes(s::Shape, b::Unity) =
    begin
        @remote(b, UnregisterForChanges(ref(s).value))
        @remote(b, UndetectCancel())
        s
    end

waiting_for_changes(s::Shape, b::Unity) =
    ! @remote(b, WasCanceled())

changed_shape(ss::Shapes, b::Unity) =
    let changed = []
        while length(changed) == 0 && ! @remote(b, WasCanceled())
            changed =  @remote(b, ChangedShape())
            sleep(0.1)
        end
        if length(changed) > 0
            shape_from_ref(changed[1], b)
        else
            nothing
        end
    end

# HACK: This should be filtered on the plugin, not here.
all_shapes(b::Unity) =
  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::Unity) =
  Shape[shape_from_ref(r, b) for r in @remote(b, GetAllShapesInLayer(layer))]

disable_update(b::Unity) =
  @remote(b, DisableUpdate())

enable_update(b::Unity) =
  @remote(b, EnableUpdate())

# Render

=#

#render exposure: [-3, +3] -> [-6, 21]
convert_render_exposure(b::Unity, v::Real) = -4.05*v + 8.8
#render quality: [-1, +1] -> [+1, +50]
convert_render_quality(b::Unity, v::Real) = round(Int, 25.5 + 24.5*v)

render_view(path::String, b::Unity) =
    let c = connection(b)
      @remote(b, SetResolution(render_width(), render_height()))
      interrupt_processing(c)
      @remote(b, ScreenShot(path))
      interrupt_processing(c)
      path
    end

highlight_shape(s::Shape, b::Unity) =
    @remote(b, SelectGameObjects(collect_ref(s)))

highlight_shapes(ss::Shapes, b::Unity) =
    @remote(b, SelectGameObjects(collect_ref(ss)))

select_position(prompt::String, b::Unity) =
  begin
    @info "$(prompt) on the $(b) backend."
    @remote(b, StartSelectingPosition())
    let s = u0() # Means not found
      while s == u0()
        sleep(0.1)
        s = @remote(b, SelectedPosition())
      end
      s
    end
  end

selected_game_objects(b) =
  begin
    while ! @remote(b, EndedSelectingGameObjects())
      sleep(0.1)
    end
    @remote(b, SelectedGameObjectsIds(true))
  end

select_shape(prompt::String, b::Unity) =
  select_one_with_prompt(prompt, b, (c, prompt) ->
    begin
      @remote(b, StartSelectingGameObject())
      selected_game_objects(b)
    end)

select_shapes(prompt::String, b::Unity) =
  select_many_with_prompt(prompt, b, (c, prompt) ->
  begin
    @remote(b, StartSelectingGameObjects())
    selected_game_objects(b)
  end)

set_sun(altitude::Real, azimuth::Real, b::Unity) =
  @remote(b, SetSun(altitude, azimuth))

#=
function SunPos(year, month, day, hour, minute, Lstm, latitude, longitude)
  if abs(longitude-Lstm)>30
     @info("Longitude $(longitude) differs by more than 30 degrees from timezone meridian $(Lstm).")
  end
  # Calculate universal time (UT)
  T = hour+(minute/60);
  UT = T-Lstm/15;
  if 0 > UT
     day = day-1;
     UT = 24+UT;
  end
  if UT > 24
     day = day+1;
     UT = UT-24;
  end
  int(x) = floor(Int, x)
  radians(x) = pi*x/180
  degrees(x) = 180*x/pi
  # Amount of days to, or from, the year 2000
  d = 367*year-int((7*int((year+int((month+9))/12)))/4)+int((275*month)/9)+day-730530+UT/24;
  # Longitude of perihelion (w), eccentricity (e)
  w = 282.9404+4.70935E-5*d;
  e = 0.016709-1.151E-9*d;
  # Mean anomaly (M), sun's mean longitude (L)
  M = 356.0470+0.9856002585*d;
  if 0 < M < 360
     M = M-floor(M/360)*360;
  end
  L = w+M;
  if (0<L<360)
     L = L-floor(L/360)*360;
  end
  # Obliquity of the ecliptic, eccentric anomaly (E)
  oblecl = 23.4393-3.563E-7*d;
  E = M+(180/pi)*e*sin(radians(M))*(1+e*cos(radians(M)));
  # Sun's rectangular coordinates in the plane of ecliptic (A,B)
  A = cos(radians(E))-e;
  B = sin(radians(E))*sqrt(1-e*e);
  # Distance (r), true anomaly (V), longitude of the sun (lon)
  r = sqrt(A*A+B*B);
  V = degrees(atan(radians(B),radians(A)));
  lon = V+w;
  if 0 < lon < 360
     lon = lon-floor(lon/360)*360;
  end
  # Calculate declination and right ascension
  decl = asin(sin(radians(oblecl))*sin(radians(lon)));
  RA = degrees(atan(sin(radians(lon))*cos(radians(oblecl)),cos(radians(lon))))/15;
  # Greenwich meridian siderial time at 00:00 (GMST0),siderial time (SIDTIME), hour angle (HA)
  GMST0 = L/15+12;
  SIDTIME = GMST0+UT+longitude/15;
  HA = (SIDTIME-RA)*15;
  # This is what we're looking for: Altitude & Azimuth
  Al = degrees(asin(sin(radians(latitude))*sin(decl)+cos(radians(latitude))*cos(decl)*cos(radians(HA))));
  Az = degrees(atan(sin(radians(HA)),cos(radians(HA))*sin(radians(latitude))-tan(decl)*cos(radians(latitude))))+180;
  Al, Az
end

SunPos(2000, 6, 21, 4, 0, 0, 51, 0)

for i in 0:23
  set_sun(SunPos(2000, 6, 21, i, 0, 0, 51, 0)..., unity)
  sleep(0.2)
end

=#


text:
using jInv.Mesh
using Test

# setup 2D problem for Nodal DivSigGrad
nc = [5, 7]
x0 = rand(2)
domain = [x0[1], 4, x0[2], 2]
M = getRegularMesh(domain,nc)
print("\ttesting forward operators for $(typeof(M))...")
Lap = getNodalLaplacianMatrix(M)
Lsig = getNodalDivSigGradMatrix(M,ones(prod(M.nc)); saveMat = false, avN2C = avN2C_Nearest);
@test norm(Lsig- Lap,Inf) <1e-12
print("passed!\n")

# setup 3D problem for Nodal DivSigGrad
nc = [5, 7, 8]
x0 = rand(3)
domain = [x0[1], 4, x0[2], 2, x0[3], 6]
h   = (domain[2:2:end]-domain[1:2:end])./nc
h1  = h[1]*ones(nc[1])
h2  = h[2]*ones(nc[2])
h3  = h[3]*ones(nc[3])

Mt = getTensorMesh3D(h1,h2,h3,x0)
Mr = getRegularMesh(domain,nc)
Mt2 = getTensorMesh3D(h1+rand(nc[1]),h2+rand(nc[2]),h3+rand(nc[3]),x0)

Meshes = (Mt,Mr,Mt2)

for k=1:length(Meshes)
	Mk = Meshes[k]
	print("\ttesting forward operators for $(typeof(M))...")
	Lap_k = getNodalLaplacianMatrix(M)
	Lsig_k = getNodalDivSigGradMatrix(M,ones(prod(M.nc)); saveMat = false, avN2C = avN2C_Nearest);
	@test norm(Lsig_k- Lap_k,Inf) <1e-12
	print("passed!\n")
end

# setup 2D problem for Face Elasticity

print("\ttesting 2D elasticity operators for $(typeof(M))...")
nc = [5, 7]
x0 = rand(2)
domain = [x0[1], 4, x0[2], 2]
M = getRegularMesh(domain,nc)
mu = ones(prod(M.n));
lambda = ones(prod(M.n));
L = GetLinearElasticityOperator(M, mu,lambda);
L2 = GetLinearElasticityOperatorFullStrain(M, mu,lambda);
println("Diff between two elasticity ops (not sure if needs to be zero):",norm(L-L2,1));
# @test norm(L- L2,Inf) <1e-8
print("passed!\n")

# setup 3D problem for Face Elasticity
print("\ttesting 3D elasticity operators for $(typeof(M))...")
n = [5,7,8];
domain = [0.0,1.0,0.0,2.0,0.0,1.0];
M = getRegularMesh(domain,n);
mu = ones(prod(M.n));
lambda = ones(prod(M.n));
L = GetLinearElasticityOperator(M, mu,lambda)
L2 = GetLinearElasticityOperatorFullStrain(M, mu,lambda)

println("Diff between two elasticity ops (not sure if needs to be zero): ",norm(L-L2,1));
# @test norm(L- L2,Inf) <1e-8
print("passed!\n")

ext: jl

max_stars_repo_name: JuliaInv/jInv

max_issues_repo_name: JuliaInv/jInv

max_forks_repo_name: JuliaInv/jInv


text:
function StopOnIncreasing()
    return StopOnIncreasing(())
end

function event_occurred(obj::StopOnIncreasing, arg0::SpacecraftState, arg1::EventDetector, arg2::jboolean)
    return jcall(obj, "eventOccurred", Action, (SpacecraftState, EventDetector, jboolean), arg0, arg1, arg2)
end


text:
module DFGPlots

using Colors
using LightGraphs
using ..GraphPlot
using DocStringExtensions
import ..GraphPlot: gplot

using ...DistributedFactorGraphs

export
    dfgplot,
    DFGPlotProps

struct DFGPlotProps
    nodefillc::NamedTuple{(:var,:fac),Tuple{RGB,RGB}}
    nodesize::NamedTuple{(:var,:fac),Tuple{Float64,Float64}}
    shape::NamedTuple{(:var,:fac),Tuple{Symbol,Symbol}} #not upported yet

    layout::Function #spring_layout, spectral_layout
    drawlabels::Bool
end

DFGPlotProps() = DFGPlotProps(  (var=colorant"seagreen", fac=colorant"cyan3"),
                                (var=1.0, fac=0.3),
                                (var=:box, fac=:elipse),
                                spring_layout,
                                true)

"""
    $(SIGNATURES)
Plots the structure of the factor graph. GraphPlot must be imported before DistributedFactorGraphs for these functions to be available.
Returns the plot context.

E.g.
```
using GraphPlot
using DistributedFactorGraphs, DistributedFactorGraphs.DFGPlots
# ... Make graph...
# Using GraphViz plotting
dfgplot(fg)
# Save to PDF
using Compose
draw(PDF("/tmp/graph.pdf", 16cm, 16cm), dfgplot(fg))
```

More information at [GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl)
"""
function dfgplot(dfg::LightDFG, p::DFGPlotProps = DFGPlotProps())

    nodetypes = [haskey(dfg.g.variables, s) for s in dfg.g.labels]

    nodesize = [isVar ? p.nodesize.var : p.nodesize.fac for isVar in nodetypes]
    # nodelabel = [isVar ? string(get_prop(dfg.g,i,:label)) : "" for (i,isVar) in enumerate(nodetypes)]
    if p.drawlabels
        nodelabel = [nodetypes[i] ? string(s) : "" for (i,s) in enumerate(dfg.g.labels)]
    else
        nodelabel = nothing
    end

    nodefillc = [isVar ? p.nodefillc.var : p.nodefillc.fac for isVar in nodetypes]

    gplot(dfg.g, nodelabel=nodelabel, nodesize=nodesize, nodefillc=nodefillc, layout=p.layout)

end

function dfgplot(dfg::AbstractDFG, p::DFGPlotProps = DFGPlotProps())
    # TODO implement convert functions
    ldfg = LightDFG{NoSolverParams}()
    copyGraph!(ldfg, dfg, listVariables(dfg), listFactors(dfg), copyGraphMetadata=false)
    dfgplot(ldfg, p)
end

function gplot(dfg::LightDFG; keyargs...)
    gplot(dfg.g; keyargs...)
end

#TODO decide if we want to overload show for display in juno, It's a bit annoying with development
# function Base.show(io::IO, ::MIME"application/prs.juno.plotpane+html", dfg::AbstractDFG)
function dfgplot(io::IO, ::MIME"application/prs.juno.plotpane+html", dfg::AbstractDFG)

    if length(ls(dfg)) != 0
        size = get(io, :juno_plotsize, [100, 100])

        plot_output = IOBuffer()
        GraphPlot.draw(GraphPlot.SVGJS(plot_output, GraphPlot.Compose.default_graphic_width,
                   GraphPlot.Compose.default_graphic_width, false), dfgplot(dfg))
        plotsvg = String(take!(plot_output))

        print(io,
            """

                $(typeof(dfg))

$(dfg.userId)

$(dfg.robotId)

$(dfg.sessionId)

$(dfg.description)

                $(plotsvg)

            """)
    else
         Base.show(io,  dfg)
    end

end

end


text:
ENV["JULIA_CUDA_SILENT"] = true
# Change the working directory to the file location
cd(@__DIR__)

using DrMZ
using Parameters: @with_kw
using Flux.NNlib
using LaTeXStrings
using ColorSchemes
using Flux
using Plots; pyplot()
using Plots.PlotMeasures
using Random
using Printf

default(fontfamily="serif",frame=:box,grid=:hide,palette=:viridis,markeralpha=0.4,dpi=200,legendfontsize=6);
PyPlot.matplotlib.rc("mathtext",fontset="cm")

@with_kw mutable struct Args
    num_sensors::Int = 128;
    L1::Float64 = 0.0;
    L2::Float64 = 2*pi;
    n_epoch::Int = 50000;
    N::Int = 128;
end

function generate_basis_function_results(pde_function;kws...)

    args = Args(;);

    M = args.num_sensors;

    # Load the trained trunk
    trunk = load_trunk(args.n_epoch,pde_function)

    # Generate quadrature grid
    x, w = gauss_legendre(args.N,args.L1,args.L2);

    # Generate custom basis functions
    basis, S = build_basis(trunk,args.L1,args.L2,M,x,w);

    # Plot singular values
    pltsing = plot((1:1:length(S)),S,label=false,seriestype=:scatter,yaxis = :log10,xlabel = L"k", ylabel = L"$\sigma_k$",xlims=(1,length(S)))
    display(pltsing)
    savefig(pltsing, @sprintf("singular_values_%s_%i_functions.png",pde_function,M))

    # Plot custom basis functions
    x_plot = (args.L1:0.01:args.L2);
    pltbasis = plot(x_plot,basis[1].(x_plot),legend=false,label=L"$\phi_1$",xlabel=L"x",xlims=(x[1],x[end]),linewidth=1.5)
    for i in 2:M
      plot!(pltbasis,x_plot,basis[i].(x_plot),legend=false,label=L"$\phi_{$i}$",linewidth=1.5)
    end
    display(pltbasis)
    savefig(pltbasis, @sprintf("basis_functions_%s_%i_functions.png",pde_function,M))

    # Compute the expansion coefficients
    ak = expansion_coefficients(basis,exp.(sin.(x)),x,w);

    # Plot the expansion coefficients
    pltak = plot((1:1:length(basis)),abs.(ak),label=false,seriestype=:scatter,yaxis = :log10,xlabel = L"k", ylabel = L"|$a_k$|")
    display(pltak)
    savefig(pltak, @sprintf("expansion_coefficients_%s_%i_functions_exp_sinx.png",pde_function,M))

end

# Generate the custom basis functions from the trained trunk functions
generate_basis_function_results("advection_equation")


text:
using ScikitLearn
using ScikitLearn.GridSearch
using PyPlot
using GaussianProcesses: GP, MeanZero, SE

srand(42)
# Training data
n = 10
x = 2π * rand(n, 1)
y = sin(x[:, 1]) + 0.05*randn(n)

# Select mean and covariance function
mZero = MeanZero()                   # Zero mean function
kern = SE(0.0,0.0)                   # Squared exponential kernel with parameters
                                     # log(ℓ) = 0.0, log(σ) = 0.0
gp = fit!(GP(m=mZero,k=kern, logNoise=-1.0), x,y);

gp_cv = fit!(GridSearchCV(GP(m=mZero,k=SE(0.0,0.0)), Dict(:logNoise=>collect(-10:0.3:10), :k_lσ=>collect(0:0.1:5))), x, y);
best_gp = gp_cv.best_estimator_;
@show get_params(best_gp)[:logNoise] get_params(best_gp)[:k_lσ]

xx = -5:0.1:10
plot(xx, predict(gp, reshape(collect(xx), length(xx), 1)), label="hand-specified")
plot(xx, predict(best_gp, reshape(collect(xx), length(xx), 1)), label="gridsearch-optimized")
plot(x, y, "bo")
legend();

get_params(gp)

#Cross-validation

using ScikitLearn, PyPlot
using ScikitLearn.CrossValidation: cross_val_score

srand(2)

n_samples = 30
degrees = [1, 4, 15]

true_fun(X) = cos(1.5 * pi * X)
X = rand(n_samples)
y = true_fun(X) + randn(n_samples) * 0.1

gp = fit!(GP(logNoise=-10.0), X'', y)

X_test = linspace(0, 1, 100)
plot(X_test, predict(gp, X_test''), label="Model")
plot(X_test, true_fun(X_test), label="True function")
scatter(X, y, label="Samples")
xlabel("x")
ylabel("y")
xlim((0, 1))
ylim((-2, +2))
legend(loc="best");

mean(cross_val_score(gp, X'', y, cv=10))


text:
@testset "unit tests" begin

    @testset "tmpdir" begin

        #=
        Leftover files from failed tests or errors can cause testing errors.
        Therefore, we remove any files in the `images` directory not pertinent to testing before executing a test.
        =#

        EXCLUDE_FILES = [".keep"]
        for f in readdir("images")
            !(f in EXCLUDE_FILES) && rm("images/$f")
        end
        @test length(readdir("images")) == length(EXCLUDE_FILES)
    end

    @testset "projection" begin
        x0 = Line(Point(0, 10), O)
        p = Point(10, 20)
        @test projection(p, x0) == Point(0, 20)

        y0 = Line(Point(10, 0), O)
        @test projection(p, y0) == Point(10, 0)
    end

    @testset "translation" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Translation(Point(1, 1), Point(100, 100)))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))
        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == Point(1, 1)
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].by == Point(50, 50)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].by == Point(100, 100)

        # with easing function
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, sineio(), Translation(Point(1, 1), Point(100, 100)))
        anim = Animation([0.0, 1.0], [1.0, 100.0], [sineio()])
        m = 49 / 99
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))
        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == Point(1, 1)
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].by == Point(at(anim, m), at(anim, m))
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].by == Point(100, 100)

        # with animation function
        anim_01 = Animation([0.0, 1.0], [0.0, 1.0], [sineio()])
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, anim_01, Translation(Point(1, 1), Point(100, 100)))
        anim = Animation([0.0, 1.0], [1.0, 100.0], [sineio()])
        m = 49 / 99
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))
        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == Point(1, 1)
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].by == Point(at(anim, m), at(anim, m))
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].by == Point(100, 100)
    end

    @testset "translation subaction" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(
            1:100,
            () -> 1;
            subactions = [SubAction(1:100, Translation(Point(1, 1), Point(100, 100)))],
        )
        # needs internal translation as well
        push!(action.subactions[1].internal_transitions, Javis.InternalTranslation(O))
        Javis.compute_transition!(action.subactions[1], video, 1)
        @test action.subactions[1].internal_transitions[1].by == Point(1, 1)
        Javis.compute_transition!(action.subactions[1], video, 50)
        @test action.subactions[1].internal_transitions[1].by == Point(50, 50)
        Javis.compute_transition!(action.subactions[1], video, 100)
        @test action.subactions[1].internal_transitions[1].by == Point(100, 100)

        ######## testing with easing function ##############
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(
            1:100,
            () -> 1;
            subactions = [SubAction(
                1:100,
                polyout(10),
                Translation(Point(1, 1), Point(100, 100)),
            )],
        )
        anim = Animation([0.0, 1.0], [1.0, 100.0], [polyout(10)])
        m = 49 / 99
        # needs internal translation as well
        push!(action.subactions[1].internal_transitions, Javis.InternalTranslation(O))
        Javis.compute_transition!(action.subactions[1], video, 1)
        @test action.subactions[1].internal_transitions[1].by == Point(1, 1)
        Javis.compute_transition!(action.subactions[1], video, 50)
        @test action.subactions[1].internal_transitions[1].by ==
              Point(at(anim, m), at(anim, m))
        Javis.compute_transition!(action.subactions[1], video, 100)
        @test action.subactions[1].internal_transitions[1].by == Point(100, 100)

        ######## testing with animation ##############
        video = Video(500, 500)
        anim_01 = Animation([0.0, 1.0], [0.0, 1.0], [polyin(10)])
        # dummy action doesn't need a real function
        action = Action(
            1:100,
            () -> 1;
            subactions = [SubAction(
                1:100,
                anim_01,
                Translation(Point(1, 1), Point(100, 100)),
            )],
        )
        anim = Animation([0.0, 1.0], [1.0, 100.0], [polyin(10)])
        m = 49 / 99
        # needs internal translation as well
        push!(action.subactions[1].internal_transitions, Javis.InternalTranslation(O))
        Javis.compute_transition!(action.subactions[1], video, 1)
        @test action.subactions[1].internal_transitions[1].by == Point(1, 1)
        Javis.compute_transition!(action.subactions[1], video, 50)
        @test action.subactions[1].internal_transitions[1].by ==
              Point(at(anim, m), at(anim, m))
        Javis.compute_transition!(action.subactions[1], video, 100)
        @test action.subactions[1].internal_transitions[1].by == Point(100, 100)
    end

    @testset "Relative frames" begin
        video = Video(500, 500)
        action = Action(Rel(10), (args...) -> 1, Translation(Point(1, 1), Point(100, 100)))
        # dummy action doesn't need a real function
        test_file = javis(
            video,
            [
                Action(1:100, (args...) -> 1, Translation(Point(1, 1), Point(100, 100))),
                action,
            ],
        )
        @test Javis.get_frames(action) == 101:110
        rm(test_file)
    end

    @testset "translation from origin" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Translation(Point(99, 99)))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == O
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].by == Point(49, 49)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].by == Point(99, 99)

        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Translation(99, 99))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == O
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].by == Point(49, 49)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].by == Point(99, 99)
    end

    @testset "rotations" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Rotation(2π))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalRotation(0.0, O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].angle == 0.0
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].angle == 2π

        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Rotation(2π, Point(2.0, 5.0)))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalRotation(0.0, O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].angle == 0.0
        @test action.internal_transitions[1].center == Point(2.0, 5.0)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].angle == 2π
        @test action.internal_transitions[1].center == Point(2.0, 5.0)
    end

    @testset "scaling" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(0:100, () -> 1, Scaling(0.0, 1.0))
        # needs internal scaling as well
        push!(action.internal_transitions, Javis.InternalScaling((0, 0)))
        Javis.compute_transition!(action, video, 0)
        @test action.internal_transitions[1].scale == (0.0, 0.0)
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].scale == (0.5, 0.5)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].scale == (1.0, 1.0)
    end

    @testset "Relative frames" begin
        video = Video(500, 500)
        action = Action(Rel(10), (args...) -> 1, Translation(Point(1, 1), Point(100, 100)))
        # dummy action doesn't need a real function
        test_file = javis(
            video,
            [
                Action(1:100, (args...) -> 1, Translation(Point(1, 1), Point(100, 100))),
                action,
            ],
        )
        @test Javis.get_frames(action) == 101:110
        rm(test_file)
    end

    @testset "translation from origin" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Translation(Point(99, 99)))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == O
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].by == Point(49, 49)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].by == Point(99, 99)

        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Translation(99, 99))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == O
        Javis.compute_transition!(action, video, 50)
        @test action.internal_transitions[1].by == Point(49, 49)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].by == Point(99, 99)
    end

    @testset "rotations" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Rotation(2π))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalRotation(0.0, O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].angle == 0.0
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].angle == 2π

        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:100, () -> 1, Rotation(2π, Point(2.0, 5.0)))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalRotation(0.0, O))

        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].angle == 0.0
        @test action.internal_transitions[1].center == Point(2.0, 5.0)
        Javis.compute_transition!(action, video, 100)
        @test action.internal_transitions[1].angle == 2π
        @test action.internal_transitions[1].center == Point(2.0, 5.0)
    end

    @testset "action with a single frame" begin
        video = Video(500, 500)
        # dummy action doesn't need a real function
        action = Action(1:1, () -> 1, Translation(Point(10, 10)))
        # needs internal translation as well
        push!(action.internal_transitions, Javis.InternalTranslation(O))
        Javis.compute_transition!(action, video, 1)
        @test action.internal_transitions[1].by == Point(10, 10)
    end

    @testset "Frames errors" begin
        # throws because a Video object was not previously defined
        empty!(Javis.CURRENT_VIDEO)
        @test_throws ErrorException Action(1:10, (args...) -> star(O, 20, 5, 0.5, 0, :fill))
        # throws because the frames of the first action need to be defined explicitly
        video = Video(500, 500)
        @test_throws ArgumentError javis(
            video,
            [Action((args...) -> 1, Translation(Point(1, 1), Point(100, 100)))],
        )
        # throws because the frames of the first action need to be defined explicitly
        @test_throws ArgumentError javis(
            video,
            [Action(Rel(10), (args...) -> 1, Translation(Point(1, 1), Point(100, 100)))],
        )
        # throws because :some is not supported as Symbol for `frames`
        @test_throws ArgumentError javis(
            video,
            [
                Action(1:100, (args...) -> 1, Translation(Point(1, 1), Point(100, 100))),
                Action(
                    :some,
                    :id,
                    (args...) -> 1,
                    Translation(Point(1, 1), Point(100, 100)),
                ),
            ],
        )
    end

    @testset "Unspecified symbol error" begin
        video = Video(500, 500)
        # throws because the frames of the first action need to be defined explicitly
        @test_throws ErrorException javis(
            video,
            [
                Action(1:100, (args...) -> 1),
                Action(1:100, (args...) -> line(O, pos(:non_existent), :stroke)),
            ],
        )

        video = Video(500, 500)
        # throws because the frames of the first action need to be defined explicitly
        @test_throws ErrorException javis(
            video,
            [
                Action(1:100, (args...) -> 1),
                Action(1:100, (args...) -> line(O, ang(:non_existent), :stroke)),
            ],
        )
    end

    @testset "Frame computation" begin

        actions = [
            BackgroundAction(1:50, (args...) -> 1),
            Action(:atom, (args...) -> 1, subactions = [SubAction(Scaling(1, 2))]),
        ]
        demo = Video(500, 500)
        javis(demo, actions; pathname = "")
        @test Javis.get_frames(actions[1]) == 1:50
        @test Javis.get_frames(actions[2]) == 1:50
        @test Javis.get_frames(actions[2].subactions[1]) == 1:50

        actions = [
            BackgroundAction(1:50, (args...) -> 1),
            Action(:atom, (args...) -> 1, subactions = [SubAction((args...) -> 1)]),
        ]
        demo = Video(500, 500)
        javis(demo, actions; pathname = "")
        @test Javis.get_frames(actions[1]) == 1:50
        @test Javis.get_frames(actions[2]) == 1:50
        @test Javis.get_frames(actions[2].subactions[1]) == 1:50

        actions = [
            BackgroundAction(1:50, (args...) -> 1),
            Action(
                Rel(-19:0),
                :atom,
                (args...) -> 1,
                subactions = [
                    SubAction(1:10, Scaling(1, 2)),
                    SubAction(Rel(10), Scaling(1, 2)),
                ],
            ),
        ]
        demo = Video(500, 500)
        javis(demo, actions; pathname = "")
        @test Javis.get_frames(actions[1]) == 1:50
        @test Javis.get_frames(actions[2]) == 31:50
        @test Javis.get_frames(actions[2].subactions[1]) == 1:10
        @test Javis.get_frames(actions[2].subactions[2]) == 11:20
    end
end


text:
# This file was generated, do not modify it. # hide
#hideall
println(read("tables/reg-voc.md", String))


text:
# Functions for generating gridworld PDDL problems
using PDDL

"Converts ASCII gridworlds to PDDL problem."
function ascii_to_pddl(str::String, name="gridworld-problem")
    rows = split(str, "\n", keepempty=false)
    width, height = maximum(length.(strip.(rows))), length(rows)
    walls = parse_pddl("(= walls (new-bit-matrix false $width $height))")
    init = Term[walls]
    start, goal = Term[], pddl"(true)"
    for (y, row) in enumerate(rows)
        for (x, char) in enumerate(strip(row))
            if char == '.' # Unoccupied
                continue
            elseif char == 'W' # Wall
                wall = parse_pddl("(= walls (set-index walls true $y $x))")
                push!(init, wall)
            elseif char == 's' # Start position
                start = parse_pddl("(= xpos $x)", "(= ypos $y)")
            elseif char == 'g' # Goal position
                goal = parse_pddl("(and (= xpos $x) (= ypos $y))")
            end
        end
    end
    append!(init, start)
    problem = GenericProblem(
        Symbol(name), :gridworld, Const[], Dict{Const,Symbol}(),
        init, goal, nothing)
    return problem
end

function load_ascii_problem(path::AbstractString)
    str = open(f->read(f, String), path)
    return ascii_to_pddl(str)
end

function convert_ascii_problem(path::String)
    str = open(f->read(f, String), path)
    str = ascii_to_pddl(str)
    new_path = join(split(path, ".")[1:end-1]) * ".pddl"
    write(new_path, write_problem(str))
    return new_path
end


text:
using InfiniteLinearAlgebra, BandedMatrices, LazyArrays, Test
import InfiniteLinearAlgebra: UpperHessenbergQ, LowerHessenbergQ, tail_de, _BandedMatrix, QL, InfToeplitz

@testset "HessenbergQ" begin
    @testset "finite UpperHessenbergQ" begin
        c,s = cos(0.1), sin(0.1); q = [c s; s -c];
        Q = UpperHessenbergQ(Fill(q,1))
        @test size(Q) == (size(Q,1),size(Q,2)) == (2,2)
        @test Q ≈ q
        Q = UpperHessenbergQ(Fill(q,2))
        @test bandwidths(Q) == (1,2)
        @test Q ≈ [q zeros(2); zeros(1,2) 1] * [1 zeros(1,2); zeros(2) q]
        @test Q' isa LowerHessenbergQ
        @test Q' ≈ [1 zeros(1,2); zeros(2) q'] * [q' zeros(2); zeros(1,2) 1]

        q = qr(randn(2,2) + im*randn(2,2)).Q
        Q = UpperHessenbergQ(Fill(q,1))
        @test Q ≈ q
        Q = UpperHessenbergQ(Fill(q,2))
        @test bandwidths(Q) == (1,2)
        @test Q ≈ [q zeros(2); zeros(1,2) 1] * [1 zeros(1,2); zeros(2) q]
        @test Q' isa LowerHessenbergQ
        @test Q' ≈ [1 zeros(1,2); zeros(2) q'] * [q' zeros(2); zeros(1,2) 1]
    end

    @testset "finite LowerHessenbergQ" begin
        c,s = cos(0.1), sin(0.1); q = [c s; s -c];
        Q = LowerHessenbergQ(Fill(q,1))
        @test size(Q) == (size(Q,1),size(Q,2)) == (2,2)
        @test Q ≈ q
        Q = LowerHessenbergQ(Fill(q,2))
        @test bandwidths(Q) == (2,1)
        @test Q ≈ [1 zeros(1,2); zeros(2) q] * [q zeros(2); zeros(1,2) 1]
        @test Q' isa UpperHessenbergQ
        @test Q' ≈ [q' zeros(2); zeros(1,2) 1]  * [1 zeros(1,2); zeros(2) q']

        q = qr(randn(2,2) + im*randn(2,2)).Q
        Q = LowerHessenbergQ(Fill(q,1))
        @test Q ≈ q
        Q = LowerHessenbergQ(Fill(q,2))
        @test bandwidths(Q) == (2,1)
        @test Q ≈ [1 zeros(1,2); zeros(2) q] * [q zeros(2); zeros(1,2) 1]
        @test Q' isa UpperHessenbergQ
        @test Q' ≈ [q' zeros(2); zeros(1,2) 1]  * [1 zeros(1,2); zeros(2) q']
    end

    @testset "infinite-Q" begin
        c,s = cos(0.1), sin(0.1); q = [c s; s -c];
        Q = LowerHessenbergQ(Fill(q,∞))
        @test Q[1,1] ≈ 0.9950041652780258
        Q = UpperHessenbergQ(Fill(q,∞))
        @test Q[1,1] ≈ 0.9950041652780258
    end

    @testset "QRPackedQ,QLPackedQ -> LowerHessenbergQ,UpperHessenbergQ" begin
        for T in (Float64, ComplexF64)
            A = brand(T,10,10,1,1)
            Q,R = qr(A)
            @test UpperHessenbergQ(Q) ≈ Q
            Q,L = ql(A)
            @test LowerHessenbergQ(Q) ≈ Q
        end
    end
end


text:
# ---
# title: 886. Possible Bipartition
# id: problem886
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Graph
# link: <https://leetcode.com/problems/possible-bipartition/description/>
# hidden: true
# ---
#
# Given a set of `N` people (numbered `1, 2, ..., N`), we would like to split
# everyone into two groups of **any** size.
#
# Each person may dislike some other people, and they should not go into the
# same group.
#
# Formally, if `dislikes[i] = [a, b]`, it means it is not allowed to put the
# people numbered `a` and `b` into the same group.
#
# Return `true` if and only if it is possible to split everyone into two groups
# in this way.
#
#
#
# **Example 1:**
#
#
#
#     Input: N = 4, dislikes = [[1,2],[1,3],[2,4]]
#     Output: true
#     **Explanation** : group1 [1,4], group2 [2,3]
#
#
# **Example 2:**
#
#
#
#     Input: N = 3, dislikes = [[1,2],[1,3],[2,3]]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: N = 5, dislikes = [[1,2],[2,3],[3,4],[4,5],[1,5]]
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= N <= 2000`
#   * `0 <= dislikes.length <= 10000`
#   * `dislikes[i].length == 2`
#   * `1 <= dislikes[i][j] <= N`
#   * `dislikes[i][0] < dislikes[i][1]`
#   * There does not exist `i != j` for which `dislikes[i] == dislikes[j]`.
#
#
## @lc code=start
using LeetCode

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


text:
import Base:
    +, show, *

doc"""N-dimensional Volume with lower and upper bounds"""
immutable Vol{N,T}
    bounds::Interval{T}
end

Vol_name(i::Integer) = i == 1 ? "length" : i == 2 ? "area" : i == 3 ? "Volume" : "hyper-Volume"

show{N,T}(io::IO, v::Vol{N,T}) = print(io, "$N-dimensional $(Vol_name(N)): $(v.bounds)")

Vol{N,T}(X::IntervalBox{N,T}) = Vol{N,T}(prod([Interval(x.hi) - Interval(x.lo) for x in X]))

Vol{T}(x::Interval{T}) = Vol(IntervalBox(x))

+{N,T}(v::Vol{N,T}, w::Vol{N,T}) = Vol{N,T}(v.bounds + w.bounds)

*{N1,N2,T}(v::Vol{N1,T}, w::Vol{N2,T}) = Vol{N1+N2,T}(v.bounds * w.bounds)

Vol{N,T}(XX::SubPaving{N,T}) = sum([Vol(X) for X in XX])

function Vol{N,T}(P::Paving{N,T})
    Vin = Vol(P.inner)
    Vout = Vin + Vol(P.boundary)

    return Vol{N,T}(hull(Vin.bounds, Vout.bounds))
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_nw_expansion(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	L = inputs["L"]     # Number of transmission lines

	# Transmission network reinforcements
	transcap = zeros(L)
	for i in 1:L
		if i in inputs["EXPANSION_LINES"]
			transcap[i] = value.(EP[:vNEW_TRANS_CAP][i])
		end
	end
	dfTransCap = DataFrame(
	Line = 1:L, New_Trans_Capacity = convert(Array{Union{Missing,Float32}}, transcap),
	Cost_Trans_Capacity = convert(Array{Union{Missing,Float32}}, transcap.*inputs["pC_Line_Reinforcement"]),
	)
	CSV.write(joinpath(path, "network_expansion.csv"), dfTransCap)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
struct AntiHermitianStorage{T,scalar} <: AbstractStorage{T, scalar}
  data::Array{T,3}
  norb::Int
  N::Int
  M::Int

  function AntiHermitianStorage{T, scalar}(data::AbstractArray{T,3}) where {T, scalar}
    norb = size(data,1)
    m = size(data,3)
    @assert norb == size(data,1) == size(data,2)
    @assert !(scalar && norb > 1)
    N = div(-1 + isqrt(1 + 8m), 2)
    @assert div(N*(N+1),2) == m
    new(data,norb,N,N)
  end
end

function AntiHermitianStorage(data::AbstractArray{T,3}, scalar=false) where T
  AntiHermitianStorage{T,scalar}(data)
end

function AntiHermitianStorage(data::AbstractArray{T,4}, scalar=false) where T
  norb = size(data,1)
  N = size(data,3)
  @assert size(data) == (norb, norb, N, N)
  M = div(N*(N+1), 2)

  lo = zeros(T, norb, norb, M)
  k = 1
  for i in 1:N
    for j in 1:i
      lo[:,:,k] .= data[:,:,i,j]
      k += 1
    end
  end

  AntiHermitianStorage(lo, scalar)
end

function AntiHermitianStorage(::Type{T}, N::Integer, norb=1, scalar=false) where T <: Number
  data = zeros(T, norb, norb, div(N*(N+1), 2))
  return AntiHermitianStorage(data, scalar)
end

AntiHermitianStorage(N::Integer, norb=1, scalar=false) = AntiHermitianStorage(ComplexF64, N, norb, scalar)

function Base.getindex(X::AntiHermitianStorage{T,scalar}, k, l, i::Int, j::Int) where {T, scalar}
  i < j && return copy(-adjoint(X[l,k,j,i])) #FIXME k,l must BOTH be integer or colons for this to make sense
  n = div((i-1)*i,2)
  n += j
  return X.data[k,l,n]
end

function Base.setindex!(X::AntiHermitianStorage{T,scalar}, v, k, l, i::Int, j::Int) where {T, scalar}
  i < j && return X[l,k,j,i] = -adjoint(v)
  n = div((i-1)*i,2)
  n += j
  X.data[k,l,n] = v
end


text:
#####
##### Canvas
#####

@testset "Canvas -- push!" begin
    c = Canvas()
    @test Variable(1) == push!(c, Expr(:call, rand))
    @test Variable(2) == push!(c, Expr(:call, rand))
    @test Variable(3) == push!(c, Expr(:call, rand))
    @test Variable(4) == push!(c, Expr(:call, rand))
    @test Variable(5) == push!(c, Expr(:call, rand))
    @test length(c) == 5
    for i in 1 : 5
        @test c.defs[i] == (i, i)
        @test c.codelocs[i] == Int32(1)
        @test getindex(c, i) == Expr(:call, rand)
    end
end

@testset "Canvas -- insert!" begin
    c = Canvas()
    @test Variable(1) == push!(c, Expr(:call, rand))
    @test Variable(2) == insert!(c, 1, Expr(:call, rand))
    @test Variable(3) == insert!(c, 1, Expr(:call, rand))
    @test Variable(4) == insert!(c, 1, Expr(:call, rand))
    @test Variable(5) == insert!(c, 1, Expr(:call, rand))
    @test length(c) == 5
    v = [5, 1, 2, 3, 4]
    for i in 1 : 5
        @test c.defs[i][2] == v[i]
    end
    c = Canvas()
    @test Variable(1) == push!(c, Expr(:call, rand))
    @test Variable(2) == pushfirst!(c, Expr(:call, rand))
    @test Variable(3) == pushfirst!(c, Expr(:call, rand))
    @test Variable(4) == pushfirst!(c, Expr(:call, rand))
    @test Variable(5) == pushfirst!(c, Expr(:call, rand))
    @test length(c) == 5
    for i in 1 : 5
        @test c.defs[i][1] == i
    end
end

@testset "Canvas -- delete!" begin
    c = Canvas()
    push!(c, Expr(:call, +, 5, 5))
    push!(c, Expr(:call, +, 10, 10))
    push!(c, Expr(:call, +, 15, 15))
    push!(c, Expr(:call, +, 20, 20))
    push!(c, Expr(:call, +, 25, 25))
    delete!(c, 4)
    @test c.defs == [(1,1), (2,2), (3,3), (4, -1), (5, 5)]
    @test c[1] == Expr(:call, +, 5, 5)
    @test c[2] == Expr(:call, +, 10, 10)
    @test c[3] == Expr(:call, +, 15, 15)
    @test c[5] == Expr(:call, +, 25, 25)
    push!(c, Expr(:call, +, 10, 10))
    delete!(c, 4)
    @test c.defs == [(1,1), (2,2), (3,3), (4, -1), (5, 5), (6, 6)]
    @test c[1] == Expr(:call, +, 5, 5)
    @test c[2] == Expr(:call, +, 10, 10)
    @test c[3] == Expr(:call, +, 15, 15)
    @test c[6] == Expr(:call, +, 10, 10)
end

@testset "Canvas -- setindex!" begin
    c = Canvas()
    push!(c, Expr(:call, +, 5, 5))
    pushfirst!(c, Expr(:call, +, 10, 10))
    pushfirst!(c, Expr(:call, +, 15, 15))
    pushfirst!(c, Expr(:call, +, 20, 20))
    pushfirst!(c, Expr(:call, +, 25, 25))

    @test getindex(c, 1) == Expr(:call, +, 5, 5)
    setindex!(c, Expr(:call, *, 10, 10), 1)
    @test getindex(c, 1) == Expr(:call, *, 10, 10)

    @test getindex(c, 2) == Expr(:call, +, 10, 10)
    setindex!(c, Expr(:call, *, 5, 5), 2)
    @test getindex(c, 2) == Expr(:call, *, 5, 5)

    @test getindex(c, 3) == Expr(:call, +, 15, 15)
    setindex!(c, Expr(:call, *, 10, 10), 3)
    @test getindex(c, 3) == Expr(:call, *, 10, 10)
end

@testset "Canvas -- renumber" begin
    c = Canvas()
    push!(c, Expr(:call, +, 5, 5))
    pushfirst!(c, Expr(:call, +, 10, 10))
    pushfirst!(c, Expr(:call, +, 15, 15))
    pushfirst!(c, Expr(:call, +, 20, 20))
    pushfirst!(c, Expr(:call, +, 25, 25))
    delete!(c, 4)
    display(c)
    println()
    c = renumber(c)
    display(c)
    println()
    @test c.defs == [(1, 1), (2, 2), (3, 3), (4, 4)]
    @test c[1] == Expr(:call, +, 10, 10)
    @test c[2] == Expr(:call, +, 15, 15)
    @test c[3] == Expr(:call, +, 25, 25)
    @test c[4] == Expr(:call, +, 5, 5)
end

@testset "Canvas -- misc." begin
    c = Canvas()
    push!(c, Expr(:call, +, 5, 5))
    pushfirst!(c, Expr(:call, +, 10, 10))
    pushfirst!(c, Expr(:call, +, 15, 15))
    pushfirst!(c, Expr(:call, +, 20, 20))
    pushfirst!(c, Expr(:call, +, 25, 25))
    insert!(c, var(1), Expr(:call, rand))
end


text:
include(joinpath(@__DIR__, ".."))

using FileIO
using ForwardDiff
using IfElse
using JLD2
using LinearAlgebra
using Rotations
using RoboDojo
import RoboDojo: LinearSolver, LUSolver, Model, ResidualMethods, Space, Disturbances, IndicesZ, InteriorPoint, EmptySolver, Policy, Trajectory, GradientTrajectory, InteriorPointOptions, IndicesOptimization, interior_point, interior_point_solve!, bilinear_violation, residual_violation, general_correction_term!, r!, rz!, rθ!, linear_solve!, lu_solver, empty_policy, empty_disturbances, friction_coefficients, SimulatorStatistics, SimulatorOptions, indices_θ, num_data, initialize_z!, initialize_θ!, indices_z, indices_θ, simulate!, policy, process!, Simulator

using StaticArrays
using Symbolics

# Utilities
include("utils.jl")

# Solver

# Environment
include("simulator/environment.jl")

# Dynamics
include("dynamics/model.jl")

# Simulator
include("simulation/index.jl")

# Simulator
include("simulation/contact_methods.jl")
include("simulation/simulation.jl")

include("dynamics/code_gen_dynamics.jl")
include("dynamics/fast_methods_dynamics.jl")

# # Models
include("dynamics/particle_2D/model.jl")
include("dynamics/particle/model.jl")
include("dynamics/hopper_2D/model.jl")
include("dynamics/hopper_3D/model.jl")
include("dynamics/quadruped/model.jl")
include("dynamics/flamingo/model.jl")
include("dynamics/pushbot/model.jl")
include("dynamics/walledcartpole/model.jl")
include("dynamics/centroidal_quadruped/model.jl")
include("dynamics/centroidal_quadruped_wall/model.jl")
include("dynamics/centroidal_quadruped_box/model.jl")

# Simulation
include("simulation/environments/flat.jl")
include("simulation/environments/piecewise.jl")
include("simulation/environments/quadratic.jl")
include("simulation/environments/slope.jl")
include("simulation/environments/sinusoidal.jl")
include("simulation/environments/stairs.jl")

include("simulation/residual_approx.jl")
include("simulation/code_gen_simulation.jl")


text:
using Mocking

create_environment_ec2_1_patch = @patch function create_environment_ec2(a...; b...)
    return Dict{String, Any}("environmentId" => "8d9967e2f0624182b74e7690ad69ebEX")
end

create_environment_membership_1_patch = @patch function create_environment_membership(a...; b...)
    return Dict{String, Any}("membership" => Dict{String, Any}("userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser", "userId" => "AIDAJ3BA6O2FMJWCWXHEX", "permissions" => "read-write", "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX"))
end

delete_environment_1_patch = @patch function delete_environment(a...; b...)
    return Dict{String, Any}()
end

delete_environment_membership_1_patch = @patch function delete_environment_membership(a...; b...)
    return Dict{String, Any}()
end

describe_environment_memberships_1_patch = @patch function describe_environment_memberships(a...; b...)
    return Dict{String, Any}("memberships" => Any[Dict{String, Any}("userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser", "userId" => "AIDAJ3BA6O2FMJWCWXHEX", "permissions" => "read-write", "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX"), Dict{String, Any}("userArn" => "arn:aws:iam::123456789012:user/MyDemoUser", "userId" => "AIDAJNUEDQAQWFELJDLEX", "permissions" => "owner", "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX")])
end

describe_environment_memberships_2_patch = @patch function describe_environment_memberships(a...; b...)
    return Dict{String, Any}("memberships" => Any[Dict{String, Any}("userArn" => "arn:aws:iam::123456789012:user/MyDemoUser", "userId" => "AIDAJNUEDQAQWFELJDLEX", "permissions" => "owner", "environmentId" => "8d9967e2f0624182b74e7690ad69ebEX")])
end

describe_environment_memberships_3_patch = @patch function describe_environment_memberships(a...; b...)
    return Dict{String, Any}("memberships" => Any[Dict{String, Any}("userArn" => "arn:aws:iam::123456789012:user/MyDemoUser", "userId" => "AIDAJNUEDQAQWFELJDLEX", "lastAccess" => "2018-01-19T11:06:13Z", "permissions" => "owner", "environmentId" => "10a75714bd494714929e7f5ec4125aEX"), Dict{String, Any}("userArn" => "arn:aws:iam::123456789012:user/MyDemoUser", "userId" => "AIDAJNUEDQAQWFELJDLEX", "lastAccess" => "2018-01-19T11:39:19Z", "permissions" => "owner", "environmentId" => "12bfc3cd537f41cb9776f8af5525c9EX")])
end

describe_environment_status_1_patch = @patch function describe_environment_status(a...; b...)
    return Dict{String, Any}("status" => "ready", "message" => "Environment is ready to use")
end

describe_environments_1_patch = @patch function describe_environments(a...; b...)
    return Dict{String, Any}("environments" => Any[Dict{String, Any}("lifecycle" => Dict{String, Any}("status" => "CREATED"), "name" => "my-demo-environment", "ownerArn" => "arn:aws:iam::123456789012:user/MyDemoUser", "id" => "8d9967e2f0624182b74e7690ad69ebEX", "type" => "ec2", "arn" => "arn:aws:cloud9:us-east-2:123456789012:environment:8d9967e2f0624182b74e7690ad69ebEX", "description" => "This is my demonstration environment."), Dict{String, Any}("lifecycle" => Dict{String, Any}("status" => "CREATED"), "name" => "another-demo-environment", "ownerArn" => "arn:aws:sts::123456789012:assumed-role/AnotherDemoUser/AnotherDemoUser", "id" => "349c86d4579e4e7298d500ff57a6b2EX", "type" => "ssh", "arn" => "arn:aws:cloud9:us-east-2:123456789012:environment:349c86d4579e4e7298d500ff57a6b2EX", "description" => "")])
end

list_environments_1_patch = @patch function list_environments(a...; b...)
    return Dict{String, Any}("environmentIds" => Any["349c86d4579e4e7298d500ff57a6b2EX", "45a3da47af0840f2b0c0824f5ee232EX"])
end

update_environment_1_patch = @patch function update_environment(a...; b...)
    return Dict{String, Any}()
end

update_environment_membership_1_patch = @patch function update_environment_membership(a...; b...)
    return Dict{String, Any}("membership" => Dict{String, Any}("userArn" => "arn:aws:iam::123456789012:user/AnotherDemoUser", "userId" => "AIDAJ3BA6O2FMJWCWXHEX", "permissions" => "read-only", "environmentId" => "8d9967e2f0624182b74e7690ad69eb31"))
end


text:
using CollegeStratBase, CollegeStratData, Test

csd = CollegeStratData;

function regr_files_test()
    @testset "Regression files" begin
        @test csd.validate_regressor_name(:inc_quartile2)
        @test csd.validate_regressor_name(:quartile3)
        @test csd.validate_regressor_name(:quality4)
        @test !csd.validate_regressor_name(:qq4)
    end
end

function wage_regr_test(wr :: csd.WageRegressions)
    @testset "Wage regressions" begin
        @test csd.max_exper(wr) > 10
        @test 2 ≤ csd.max_exper_exponent(wr) ≤ 4
        @test length(ed_symbols(wr)) == n_school(wr)
        @test isa(ed_symbols(wr), Vector{Symbol})
	end
end

@testset "Wage Regressions" begin
    wage_regr_test(csd.default_wage_regressions());
    regr_files_test()
end

# ------------


text:
# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory

using MUMPS
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end
include("getDivGrad.jl");

A  = getDivGrad(24,23,25);
A2 = getDivGrad(34,32,36);

n  = size(A,1);
n2 = size(A2,1);
A  = A + im*spdiagm(rand(n),0);

nrhs = 10;
rhs = randn(n,nrhs) + im*randn(n,nrhs);
rhs2 = randn(n2,nrhs);

println("Factorize complex matrix")
tic();
F1 = factorMUMPS(A,1);
toc();

println("Factorize real matrix")
tic();
F2 = factorMUMPS(A2,1);
toc();

println("Solve complex system")
tic();
x = applyMUMPS(F1,rhs);
toc();
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A*x[:,i]-rhs[:,i]) / norm(rhs[:,i]);
end
@test maximum(err) < 1e-14

println("Solve real system")
tic();
x2 = applyMUMPS(F2,rhs2);
toc();
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A2*x2[:,i]-rhs2[:,i]) / norm(rhs2[:,i]);
end
@test maximum(err) < 1e-14

println("Free memory")
destroyMUMPS(F1);
destroyMUMPS(F2);

println("DONE!")


text:
module GLFWBackend

using Libdl
using GLFW
using ..LibCImGui
using CImGui: FLT_MAX, AddInputCharacter
using CImGui: GetIO, Set_KeyMap, Set_KeysDown, Get_KeysDown, Set_MouseDown, GetMouseCursor

include("callbacks.jl")
export SetCustomMouseButtonCallback
export SetCustomScrollCallback
export SetCustomKeyCallback
export SetCustomCharCallback

include("impl.jl")
export GlfwClientApi, GlfwClientApi_Unknown, GlfwClientApi_OpenGL, GlfwClientApi_Vulkan
export ImGui_ImplGlfw_Init
export ImGui_ImplGlfw_InitForOpenGL, ImGui_ImplGlfw_InitForVulkan
export ImGui_ImplGlfw_Shutdown
export ImGui_ImplGlfw_NewFrame

function __init__()
    global g_Window = C_NULL
    global g_ClientApi = GlfwClientApi_Unknown
    global g_Time = 0.0
    global g_MouseJustPressed = [false, false, false, false, false]
    global g_MouseCursors = [GLFW.Cursor(C_NULL) for i = 1:ImGuiMouseCursor_COUNT]
    global g_ImplGlfw_GetClipboardText = dlsym(dlopen(GLFW.libglfw), :glfwGetClipboardString)
    global g_ImplGlfw_SetClipboardText = dlsym(dlopen(GLFW.libglfw), :glfwSetClipboardString)
    global g_CustomCallbackMousebutton = C_NULL
    global g_CustomCallbackScroll = C_NULL
    global g_CustomCallbackKey = C_NULL
    global g_CustomCallbackChar = C_NULL
end

end # module


text:
abstract type AbstractInstruction end
abstract type AbstractBFInstruction <: AbstractInstruction end

const OP_CODES = Symbol.(['+', '-', '>', '<', '.', ',', '[', ']'])
const INSTRUCTIONS = (:Increment, :Decrement, :RightShift, :LeftShift, :Write, :Read, :LoopStart, :LoopEnd)

for (op_code, instruction) in zip(OP_CODES, INSTRUCTIONS)
    @eval begin
        struct $instruction <: AbstractBFInstruction end
        lex(::Val{Symbol($(string(op_code)))}) = $instruction()
    end
end

struct UnknownInstruction <: AbstractInstruction end
lex(::Val{x}) where {x} = UnknownInstruction()

lex(program::Vector{Char}) = @. lex(Val(Symbol(program)))


text:
"""
```
mutation(m, data, p, d, blocks_free, blocks_all, ϕ_n, ϕ_n1; c = 1., α = 1., old_data)
```

Execute one proposed move of the Metropolis-Hastings algorithm for a given parameter

### Arguments:
- `m::AbstractDSGEModel`: Model of type AbstractDSGEModel being estimated.
- `data::Matrix{Float64}`: Matrix of data
- `p::Vector{Float64}`: Initial particle value
- `d::Distribution`: A distribution with μ = the weighted mean, and Σ = the weighted variance/covariance matrix
- `blocks_free::Vector{Vector{Int64}}`: A vector of index blocks, where the indices in each block corresponds to the ordering of free parameters only (e.g. all the indices will be ∈ 1:n_free_parameters)
- `blocks_all::Vector{Vector{Int64}}`: A vector of index blocks, where the indices in each block corresponds to the ordering of all parameters (e.g. all the indices will be in ∈ 1:n_para, free and fixed)
- `ϕ_n::Float64`: The current tempering factor
- `ϕ_n1::Float64`: The previous tempering factor

### Keyword Arguments:
- `c::Float64`: The scaling parameter for the proposed covariance matrix
- `α::Float64`: The mixing proportion
- `old_data::Matrix{Float64}`: The matrix of old data to be used in calculating the old_loglh, old_logpost in time tempering

### Outputs:

- `p::Vector{Float64}`: An updated particle containing updated parameter values, log-likelihood, prior, and acceptance indicator.

"""
function mutation(m::AbstractDSGEModel, data::Matrix{S}, p::Vector{S},
                  d_μ::Vector{S}, d_Σ::Matrix{S},#d::Distribution,
                  blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}},
                  ϕ_n::S, ϕ_n1::S; c::S = 1., α::S = 1.,
                  old_data::T = T(undef, size(data, 1), 0),
                  use_chand_recursion::Bool = false,
                  verbose::Symbol = :low) where {S<:Float64, T<:AbstractMatrix}

    n_steps     = get_setting(m, :n_mh_steps_smc)
    n_free_para = length([!θ.fixed for θ in m.parameters])
    step_prob   = rand() # Draw initial step probability

    N         = length(p)
    para      = p[1:ind_para_end(N)]
    like      = p[ind_loglh(N)]
    logprior  = p[ind_logprior(N)]
    like_prev = p[ind_old_loglh(N)] # Likelihood evaluated at the old data (for time tempering)
    accept    = 0.0

    d = MvNormal(d_μ, d_Σ)

    for step in 1:n_steps
        for (block_f, block_a) in zip(blocks_free, blocks_all)
            # Index out parameters corresponding to given random block, create distribution
            # centered at weighted mean, with Σ corresponding to the same random block
            para_subset = para[block_a]
            d_subset    = MvNormal(d.μ[block_f], d.Σ.mat[block_f, block_f])

            para_draw, para_new_density, para_old_density = mvnormal_mixture_draw(para_subset,
                                                                       d_subset; c = c, α = α)
            para_new = deepcopy(para)
            para_new[block_a] = para_draw

            like_init  = like
            prior_init = logprior

            q0 = α * exp(logpdf(MvNormal(para_draw,   c^2 * d_subset.Σ.mat), para_subset))
            q1 = α * exp(logpdf(MvNormal(para_subset, c^2 * d_subset.Σ.mat), para_draw))

            ind_pdf = 1.0

            for i=1:length(block_a)
                Σ_ii    = sqrt(d_subset.Σ.mat[i,i])
                zstat   = (para_subset[i] - para_draw[i]) / Σ_ii
                ind_pdf = ind_pdf / (Σ_ii * sqrt(2.0 * π)) * exp(-0.5 * zstat^2)
            end

            q0 += (1.0 - α) / 2.0 * ind_pdf
            q1 += (1.0 - α) / 2.0 * ind_pdf

            q0 += (1.0 - α) / 2.0 * exp(logpdf(MvNormal(d_subset.μ, c^2 * d_subset.Σ.mat),
                                               para_subset))
            q1 += (1.0 - α) / 2.0 * exp(logpdf(MvNormal(d_subset.μ, c^2 * d_subset.Σ.mat),
                                               para_draw))
            q0 = log(q0)
            q1 = log(q1)

            prior_new = like_new = like_old_data = -Inf

            try
                update!(m, para_new)
                para_new  = [θ.value for θ in m.parameters]
                prior_new = prior(m)
                like_new  = likelihood(m, data; sampler = true,
                                       use_chand_recursion = use_chand_recursion,
                                       verbose = verbose)
                if like_new == -Inf
                    prior_new = like_old_data = -Inf
                end
                like_old_data = isempty(old_data) ? 0. : likelihood(m, old_data; sampler = true,
                                                      use_chand_recursion = use_chand_recursion,
                                                      verbose = verbose)
            catch err
                if isa(err, ParamBoundsError) || isa(err, LinearAlgebra.LAPACKException)
                    prior_new = like_new = like_old_data = -Inf
                elseif isa(err, PosDefException) || isa(err, SingularException) || isa(err, DSGE.SteadyStateConvergenceError)
                    prior_new = like_new = like_old_data = -Inf
                else
                    throw(err)
                end
            end

            if (q0 == Inf && q1 == Inf)
                q0 = 0.0
            end

            η = exp(ϕ_n * (like_new - like_init) + (1 - ϕ_n) * (like_old_data - like_prev) + (prior_new - prior_init) + (q0 - q1))

            if step_prob < η
                para      = para_new
                like      = like_new
                logprior  = prior_new
                like_prev = like_old_data
                accept   += length(block_a)
            end
            step_prob = rand() # Draw again for next step
        end
    end
    update_mutation!(p, para, like, logprior, like_prev, accept / n_free_para)
    return p
end


text:
import Oceananigans.Architectures: architecture

const R_Earth = 6371.0e3    # Mean radius of the Earth [m] https://en.wikipedia.org/wiki/Earth

struct RegularLatitudeLongitudeGrid{FT, TX, TY, TZ, A} <: AbstractHorizontallyCurvilinearGrid{FT, TX, TY, TZ}
        Nx :: Int
        Ny :: Int
        Nz :: Int
        Hx :: Int
        Hy :: Int
        Hz :: Int
        Lx :: FT
        Ly :: FT
        Lz :: FT
        Δλ :: FT
        Δφ :: FT
        Δz :: FT
      λᶠᵃᵃ :: A
      λᶜᵃᵃ :: A
      φᵃᶠᵃ :: A
      φᵃᶜᵃ :: A
      zᵃᵃᶠ :: A
      zᵃᵃᶜ :: A
    radius :: FT
end

function RegularLatitudeLongitudeGrid(FT=Float64; size, latitude, longitude, z, radius=R_Earth, halo=(1, 1, 1))
    @assert length(latitude) == 2
    @assert length(longitude) == 2
    @assert length(z) == 2

    λ₁, λ₂ = longitude
    @assert λ₁ < λ₂ && λ₂ - λ₁ ≤ 360

    φ₁, φ₂ = latitude
    @assert -90 <= φ₁ < φ₂ <= 90

    (φ₁ == -90 || φ₂ == 90) &&
        @warn "Are you sure you want to use a latitude-longitude grid with a grid point at the pole?"

    z₁, z₂ = z
    @assert z₁ < z₂

    Lλ = λ₂ - λ₁
    Lφ = φ₂ - φ₁
    Lz = z₂ - z₁

    TX = Lλ == 360 ? Periodic : Bounded
    TY = Bounded
    TZ = Bounded
    topo = (TX, TY, TZ)

    Nλ, Nφ, Nz = N = validate_size(TX, TY, TZ, size)
    Hλ, Hφ, Hz = H = validate_halo(TX, TY, TZ, halo)

            Λ₁ = (λ₁, φ₁, z₁)
            L  = (Lλ, Lφ, Lz)
    Δλ, Δφ, Δz = Δ = @. L / N

    # Calculate end points for cell faces and centers
    λF₋, φF₋, zF₋ = ΛF₋ = @. Λ₁ - H * Δ
    λF₊, φF₊, zF₊ = ΛF₊ = @. ΛF₋ + total_extent(topo, H, Δ, L)

    λC₋, φC₋, zC₋ = ΛC₋ = @. ΛF₋ + Δ / 2
    λC₊, φC₊, zC₊ = ΛC₊ = @. ΛC₋ + L + Δ * (2H - 1)

    TFλ, TFφ, TFz = total_length.(Face, topo, N, H)
    TCλ, TCφ, TCz = total_length.(Center, topo, N, H)

    λᶠᵃᵃ = range(λF₋, λF₊, length = TFλ)
    φᵃᶠᵃ = range(φF₋, φF₊, length = TFφ)
    zᵃᵃᶠ = range(zF₋, zF₊, length = TFz)

    λᶜᵃᵃ = range(λC₋, λC₊, length = TCλ)
    φᵃᶜᵃ = range(φC₋, φC₊, length = TCφ)
    zᵃᵃᶜ = range(zC₋, zC₊, length = TCz)

    λᶠᵃᵃ = OffsetArray(λᶠᵃᵃ, -Hλ)
    φᵃᶠᵃ = OffsetArray(φᵃᶠᵃ, -Hφ)
    zᵃᵃᶠ = OffsetArray(zᵃᵃᶠ, -Hz)

    λᶜᵃᵃ = OffsetArray(λᶜᵃᵃ, -Hλ)
    φᵃᶜᵃ = OffsetArray(φᵃᶜᵃ, -Hφ)
    zᵃᵃᶜ = OffsetArray(zᵃᵃᶜ, -Hz)

    return RegularLatitudeLongitudeGrid{FT, TX, TY, TZ, typeof(λᶠᵃᵃ)}(Nλ, Nφ, Nz, Hλ, Hφ, Hz, Lλ, Lφ, Lz, Δλ, Δφ, Δz, λᶠᵃᵃ, λᶜᵃᵃ, φᵃᶠᵃ, φᵃᶜᵃ, zᵃᵃᶠ, zᵃᵃᶜ, radius)
end

function domain_string(grid::RegularLatitudeLongitudeGrid)
    λ₁, λ₂ = domain(topology(grid, 1), grid.Nx, grid.λᶠᵃᵃ)
    φ₁, φ₂ = domain(topology(grid, 2), grid.Ny, grid.φᵃᶠᵃ)
    z₁, z₂ = domain(topology(grid, 3), grid.Nz, grid.zᵃᵃᶠ)
    return "longitude λ ∈ [$λ₁, $λ₂], latitude ∈ [$φ₁, $φ₂], z ∈ [$z₁, $z₂]"
end

function show(io::IO, g::RegularLatitudeLongitudeGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ}
    print(io, "RegularLatitudeLongitudeGrid{$FT, $TX, $TY, $TZ}\n",
              "                   domain: $(domain_string(g))\n",
              "                 topology: ", (TX, TY, TZ), '\n',
              "  resolution (Nx, Ny, Nz): ", (g.Nx, g.Ny, g.Nz), '\n',
              "   halo size (Hx, Hy, Hz): ", (g.Hx, g.Hy, g.Hz), '\n',
              "grid spacing (Δλ, Δφ, Δz): ", (g.Δλ, g.Δφ, g.Δz))
end

#####
##### << Spherical nodes >>
#####

# Node by node
@inline xnode(::Center, i, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.λᶜᵃᵃ[i]
@inline xnode(::Face,   i, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.λᶠᵃᵃ[i]

@inline ynode(::Center, j, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.φᵃᶜᵃ[j]
@inline ynode(::Face,   j, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.φᵃᶠᵃ[j]

@inline znode(::Center, k, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.zᵃᵃᶠ[k]
@inline znode(::Face,   k, grid::RegularLatitudeLongitudeGrid) = @inbounds grid.zᵃᵃᶜ[k]

all_x_nodes(::Type{Center}, grid::RegularLatitudeLongitudeGrid) = grid.λᶜᵃᵃ
all_x_nodes(::Type{Face},   grid::RegularLatitudeLongitudeGrid) = grid.λᶠᵃᵃ
all_y_nodes(::Type{Center}, grid::RegularLatitudeLongitudeGrid) = grid.φᵃᶜᵃ
all_y_nodes(::Type{Face},   grid::RegularLatitudeLongitudeGrid) = grid.φᵃᶠᵃ
all_z_nodes(::Type{Center}, grid::RegularLatitudeLongitudeGrid) = grid.zᵃᵃᶜ
all_z_nodes(::Type{Face},   grid::RegularLatitudeLongitudeGrid) = grid.zᵃᵃᶠ

architecture(::RegularLatitudeLongitudeGrid) = nothing


text:
export my_spy

"""
`my_spy(A::Matrix)` creates a black-and-white image of `A` in which
zeros are white squares and nonzeros are black.
"""
function my_spy(A::Matrix{T}, borders::Bool = false) where {T<:Number}
    newdraw()
    edge_color = :white
    if borders
        edge_color = :black
    end
    r, c = size(A)
    for i = 1:r
        for j = 1:c
            x = j
            y = -i
            if A[i, j] == 0
                draw_rectangle(
                    x - 1,
                    y + 1,
                    x,
                    y,
                    line = borders,
                    color = edge_color,
                    fill = true,
                    fillcolor = :white,
                )
            else
                draw_rectangle(
                    x - 1,
                    y + 1,
                    x,
                    y,
                    line = borders,
                    color = edge_color,
                    fill = true,
                    fillcolor = :black,
                )
            end
        end
    end
    finish()
end


text:
function runcase(t::Tuple)

  N = 20
  Np1 = N + 1
  Nhp1 = Int(N/2) + 1

  data = Dict(
    # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
    :struc_el => Frame(N, Np1, 3, 1, 1, Line(2, 3)),
    :properties => [1.0e6 1.0e6 1.0e6 3.0e5;],
    :x_coords => range(0, stop=4, length=Np1),
    :y_coords => zeros(Np1),
    :z_coords => zeros(Np1),
    :g_num => [
      collect(1:N)';
      collect(2:Np1)'],
    :support => [
      (1, [1 0 0 0 0 0]),
      (Np1, [0 0 0 0 0 0]),
      ],
    :loaded_nodes => [
      (Nhp1, [0.0 -10000.0 0.0 0.0 0.0 0.0])]
  )

  println("Case wob = $(t[2]), incl = $(t[3]) done!")

  t[4](data)

 end