text:
export unity, fast_unity,
       unity_material_family

macro unity_str(str)
    rpc("Unity", str)
end

# 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
=#

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

backend_name(b::Unity) = "Unity"

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

unity"public void SetApplyMaterials(bool apply)"
unity"public void SetApplyColliders(bool apply)"

fast_unity() =
  let c = connection(unity)
    UnitySetApplyMaterials(c, false)
    UnitySetApplyColliders(c, false)
    c
  end

#=
realize(b::Unity, s::Point) =
  UnityPoint(connection(b), s.position)
realize(b::Unity, s::Line) =
  UnityPolyLine(connection(b), s.vertices)
realize(b::Unity, s::Spline) =
  if (s.v0 == false) && (s.v1 == false)
    #UnitySpline(connection(b), s.points)
    UnityInterpSpline(connection(b),
                     s.points,
                     s.points[2]-s.points[1],
                     s.points[end]-s.points[end-1])
  elseif (s.v0 != false) && (s.v1 != false)
    UnityInterpSpline(connection(b), s.points, s.v0, s.v1)
  else
    UnityInterpSpline(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) =
  UnityInterpClosedSpline(connection(b), s.points)
realize(b::Unity, s::Circle) =
  UnityCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)
realize(b::Unity, s::Arc) =
  if s.radius == 0
    UnityPoint(connection(b), s.center)
  elseif s.amplitude == 0
    UnityPoint(connection(b), s.center + vpol(s.radius, s.start_angle, s.center.cs))
  elseif abs(s.amplitude) >= 2*pi
    UnityCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)
  else
    end_angle = s.start_angle + s.amplitude
    if end_angle > s.start_angle
      UnityArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, end_angle)
    else
      UnityArc(connection(b), 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
    UnityEllipse(connection(b), s.center, vz(1, s.center.cs), vxyz(s.radius_x, 0, 0, s.center.cs), s.radius_y/s.radius_x)
  else
    UnityEllipse(connection(b), 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) =
  UnityClosedPolyLine(connection(b), s.vertices)
realize(b::Unity, s::RegularPolygon) =
  UnityClosedPolyLine(connection(b), regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed))
realize(b::Unity, s::Rectangle) =
  UnityClosedPolyLine(
    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) =
  UnitySurfaceCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)
realize(b::Unity, s::SurfaceArc) =
    #UnitySurfaceArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, s.start_angle + s.amplitude)
    if s.radius == 0
        UnityPoint(connection(b), s.center)
    elseif s.amplitude == 0
        UnityPoint(connection(b), s.center + vpol(s.radius, s.start_angle, s.center.cs))
    elseif abs(s.amplitude) >= 2*pi
        UnitySurfaceCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)
    else
        end_angle = s.start_angle + s.amplitude
        if end_angle > s.start_angle
            UnitySurfaceFromCurves(connection(b),
                [UnityArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, end_angle),
                 UnityPolyLine(connection(b), [add_pol(s.center, s.radius, end_angle),
                                              add_pol(s.center, s.radius, s.start_angle)])])
        else
            UnitySurfaceFromCurves(connection(b),
                [UnityArc(connection(b), s.center, vz(1, s.center.cs), s.radius, end_angle, s.start_angle),
                 UnityPolyLine(connection(b), [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) = UnityCircle(connection(b),
#realize(b::Unity, s::SurfaceEllipse) = UnityCircle(connection(b),
=#

unity"public GameObject SurfacePolygon(Vector3[] ps)"

realize(b::Unity, s::SurfacePolygon) =
  UnitySurfacePolygon(connection(b), reverse(s.vertices))

backend_fill(b::Unity, path::ClosedPolygonalPath) =
  UnitySurfacePolygon(connection(b), path.vertices)

#=
realize(b::Unity, s::SurfaceRegularPolygon) =
  UnitySurfaceClosedPolyLine(connection(b), regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed))
realize(b::Unity, s::SurfaceRectangle) =
  UnitySurfaceClosedPolyLine(
    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::Surface) =
  let #ids = map(r->UnityNurbSurfaceFrom(connection(b),r), UnitySurfaceFromCurves(connection(b), collect_ref(s.frontier)))
      ids = UnitySurfaceFromCurves(connection(b), collect_ref(s.frontier))
    foreach(mark_deleted, s.frontier)
    ids
  end
backend_surface_boundary(b::Unity, s::Shape2D) =
    map(shape_from_ref, UnityCurvesFromSurface(connection(b), 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(UnitySurfaceDomain(connection(b), 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) = UnitySurfaceDomain(conn, r)
        map_division(u1, u2, nu) do u
            map_division(v1, v2, nv) do v
                f(UnitySurfaceFrameAt(conn, 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) = UnitySurfaceDomain(conn, r)
        map_division(u1, u2, nu) do u
            map_division(v1, v2, nv) do v
                f(UnitySurfaceFrameAt(conn, r, u, v))
            end
        end
    end
=#

unity"public GameObject Text(string txt, Vector3 position, Vector3 vx, Vector3 vy, string fontName, int fontSize)"

realize(b::Unity, s::Text) =
  UnityText(
    connection(b),
    s.str, s.corner, vz(-1, s.corner.cs), vy(1, s.corner.cs), "Fonts/Inconsolata-Regular", s.height)

unity"public GameObject Sphere(Vector3 center, float radius)"

realize(b::Unity, s::Sphere) =
  UnitySphere(connection(b), s.center, s.radius)

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

unity"public GameObject Pyramid(Vector3[] ps, Vector3 q)"
unity"public GameObject PyramidFrustum(Vector3[] ps, Vector3[] qs)"

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

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

realize(b::Unity, s::RegularPyramid) =
  UnityPyramid(connection(b),
               regular_polygon_vertices(s.edges, s.cb, s.rb, s.angle, s.inscribed),
               add_z(s.cb, s.h))

realize(b::Unity, s::IrregularPyramid) =
  UnityPyramid(connection(b), s.bs, s.t)

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

realize(b::Unity, s::IrregularPyramidFrustum) =
    UnityPyramidFrustum(connection(b), s.bs, s.ts)

realize(b::Unity, s::IrregularPrism) =
  UnityPyramidFrustum(connection(b),
                      s.bs,
                      map(p -> (p + s.v), s.bs))

unity"public GameObject RightCuboid(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz, float angle)"

realize(b::Unity, s::RightCuboid) =
  UnityRightCuboid(connection(b), s.cb, vz(1, s.cb.cs), vx(1, s.cb.cs), s.height, s.width, s.h, s.angle)

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

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

realize(b::Unity, s::Cone) =
  UnityPyramid(connection(b), regular_polygon_vertices(64, s.cb, s.r), add_z(s.cb, s.h))

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

unity"public GameObject Cylinder(Vector3 bottom, float radius, Vector3 top)"

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

#=
backend_extrusion(b::Unity, s::Shape, v::Vec) =
    and_mark_deleted(
        map_ref(s) do r
            UnityExtrude(connection(b), 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
      UnitySweep(connection(b), path_r, profile_r, rotation, scale)
    end
  end

realize(b::Unity, s::Revolve) =
  and_delete_shape(
    map_ref(s.profile) do r
      UnityRevolve(connection(b), 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)

=#
unity"public GameObject Unite(GameObject s0, GameObject s1)"
unity"public GameObject Intersect(GameObject s0, GameObject s1)"
unity"public GameObject Subtract(GameObject s0, GameObject s1)"
unity"public void SubtractFrom(GameObject s0, GameObject s1)"

unity"public GameObject Canonicalize(GameObject s)"
###
unite_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    ensure_ref(b, UnityUnite(connection(b), r0.value, r1.value))

intersect_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    ensure_ref(b, UnityIntersect(connection(b), r0.value, r1.value))

subtract_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    let r = UnitySubtract(connection(b), r0.value, r1.value)
      UnityDeleteMany(connection(b), [r0.value, r1.value])
      r
    end

#=
subtract_ref(b::Unity, r0::UnityNativeRef, r1::UnityNativeRef) =
    begin
      UnitySubtractFrom(connection(b), r0.value, r1.value)
      r0.value
    end
=#

#=
slice_ref(b::Unity, r::UnityNativeRef, p::Loc, v::Vec) =
    (UnitySlice(connection(b), 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)
    #UnityCanonicalize(connection(b), 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)

unity"public void Move(GameObject s, Vector3 v)"
unity"public void Scale(GameObject s, Vector3 p, float scale)"
unity"public void Rotate(GameObject s, Vector3 p, Vector3 n, float a)"

realize(b::Unity, s::Move) =
  let r = map_ref(s.shape) do r
            UnityMove(connection(b), r, s.v)
            r
          end
    mark_deleted(s.shape)
    r
  end
#=
realize(b::Unity, s::Transform) =
  let r = map_ref(s.shape) do r
            UnityTransform(connection(b), r, s.xform)
            r
          end
    mark_deleted(s.shape)
    r
  end
=#
realize(b::Unity, s::Scale) =
  let r = map_ref(s.shape) do r
            UnityScale(connection(b), 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
            UnityRotate(connection(b), 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
                    UnityMirror(connection(b), 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
            UnityMirror(connection(b), r, s.p, s.n, true)
          end
    UnionRef((r0,r1))
  end
=#
unity"public GameObject SurfaceFromGrid(int m, int n, Vector3[] pts, bool closedM, bool closedN, int level)"

realize(b::Unity, s::SurfaceGrid) =
    UnitySurfaceFromGrid(
        connection(b),
        size(s.points,1),
        size(s.points,2),
        reshape(s.points,:),
        s.closed_u,
        s.closed_v,
        2)
#=
realize(b::Unity, s::Thicken) =
  and_delete_shape(
    map_ref(s.shape) do r
      UnityThicken(connection(b), 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) = UnityCurveFrameAt(connection(b), 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) = UnitySurfaceFrameAt(connection(b), ref(s).value, u, v)

=#

# BIM

unity"public GameObject LoadResource(String name)"
unity"public Material LoadMaterial(String name)"
unity"public void SetCurrentMaterial(Material material)"

unity"public GameObject InstantiateResource(GameObject family, Vector3 pos, Vector3 vx, Vector3 vy, float scale)"
unity"public GameObject InstantiateBIMElement(GameObject family, Vector3 pos, float angle)"

# Families

abstract type UnityFamily <: Family end

struct UnityMaterialFamily <: UnityFamily
  name::String
  parameter_map::Dict{Symbol,String}
  ref::Parameter{Any}
end

unity_material_family(name, pairs...) = UnityMaterialFamily(name, Dict(pairs...), Parameter{Any}(nothing))
backend_get_family_ref(b::Unity, f::Family, uf::UnityMaterialFamily) = UnityLoadMaterial(connection(b), uf.name)

struct UnityResourceFamily <: UnityFamily
  name::String
  parameter_map::Dict{Symbol,String}
  ref::Parameter{Any}
end

unity_resource_family(name, pairs...) = UnityResourceFamily(name, Dict(pairs...), Parameter{Any}(nothing))
backend_get_family_ref(b::Unity, f::Family, uf::UnityResourceFamily) = UnityLoadResource(connection(b), uf.name)

set_backend_family(default_wall_family(), unity, unity_material_family("Materials/Plaster/Plaster1"))
set_backend_family(default_slab_family(), unity, unity_material_family("Materials/Concrete/Concrete2"))
set_backend_family(default_roof_family(), unity, unity_material_family("Materials/Concrete/Concrete2"))
set_backend_family(default_beam_family(), unity, unity_material_family("Materials/Metal/Aluminum"))
set_backend_family(default_column_family(), unity, unity_material_family("Materials/Concrete/Concrete2"))
set_backend_family(default_door_family(), unity, unity_material_family("Materials/Wood/InteriorWood2"))
set_backend_family(default_panel_family(), unity, unity_material_family("Materials/Glass/Glass"))

unity"public GameObject Window(Vector3 position, Quaternion rotation, float dx, float dy, float dz)"
unity"public GameObject Shelf(Vector3 position, int rowLength, int lineLength, float cellWidth, float cellHeight, float cellDepth)"

set_backend_family(default_table_family(), unity, unity_resource_family("Prefabs/Tables/ModernTable/ModernTable"))
set_backend_family(default_chair_family(), unity, unity_resource_family("Prefabs/Chairs/ModernChair/ModernChair"))
set_backend_family(default_table_chair_family(), unity, unity_resource_family("Prefabs/TablesChairs/ModernTableChair/ModernTableChair"))

backend_rectangular_table(b::Unity, c, angle, family) =
    UnityInstantiateBIMElement(connection(b), realize(b, family), c, -angle)

backend_chair(b::Unity, c, angle, family) =
    UnityInstantiateBIMElement(connection(b), realize(b, family), c, -angle)

backend_rectangular_table_and_chairs(b::Khepri.Unity, c, angle, family) =
    UnityInstantiateBIMElement(connection(b), realize(b, family), c, -angle)

unity"public GameObject Slab(Vector3[] contour, Vector3[][] holes, float h, Material material)"

backend_slab(b::Unity, profile, holes, thickness, family) =
  let bot_vs = path_vertices(profile)
      c = connection(b)
    UnitySlab(c, bot_vs, map(path_vertices, holes), thickness, realize(b, family))
  end

unity"public GameObject BeamRectSection(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz, float angle, Material material)"
unity"public GameObject BeamCircSection(Vector3 bot, float radius, Vector3 top, Material material)"

realize(b::Unity, s::Beam) =
  let profile = s.family.profile
      profile_u0 = profile.corner
      c = add_xy(s.cb, profile_u0.x + profile.dx/2, profile_u0.y + profile.dy/2)
    UnityBeamRectSection(
      connection(b),
      c, vz(1, c.cs), vx(1, c.cs),
      profile.dy, profile.dx, s.h, -s.angle,
      realize(b, s.family))
  end

#Columns are aligned along the center axis.
realize(b::Unity, s::FreeColumn) =
  let profile = s.family.profile
      profile_u0 = profile.corner
      c = add_xy(s.cb, profile_u0.x + profile.dx/2, profile_u0.y + profile.dy/2)
    UnityBeamRectSection(
      connection(b),
      c, vz(1, c.cs), vx(1, c.cs),
      profile.dy, profile.dx, s.h, -s.angle,
      realize(b, s.family))
  end

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

unity"public GameObject Panel(Vector3[] pts, Vector3 n, Material material)"

realize(b::Unity, s::Panel) =
  let #p1 = s.vertices[1],
      #p2 = s.vertices[2],
      #p3 = s.vertices[3],
      #n = vz(s.family.thickness, cs_from_o_vx_vy(p1, p2-p1, p3-p1))
      verts = in_world.(s.vertices)
      n = vertices_normal(verts)*(s.family.thickness/2)
    UnityPanel(
      connection(b),
      map(p -> p - n, verts),
      n*2,
      realize(b, s.family))
  end

sweep_fractions(b, verts, height, thickness) =
  let p = add_z(verts[1], height/2),
      q = add_z(verts[2], height/2),
      (c, h) = position_and_height(p, q),
      s = UnityNativeRef(UnityRightCuboid(connection(b), 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, thickness)...)
    else
      (s, )
    end
  end

backend_wall(b::Unity, path, height, thickness, family) =
  path_length(path) < path_tolerance() ?
    UnityEmptyRef() :
    let c = connection(b)
      UnitySetCurrentMaterial(c, realize(b, family))
      backend_wall_path(
          b,
          path,
          height*0.999, #We reduce height just a bit to avoid Z-fighting
          thickness)
    end

backend_wall_path(b::Unity, path::PolygonalPath, height, thickness) =
    UnityUnionRef(sweep_fractions(b, path.vertices, height, thickness))

backend_wall_path(b::Unity, path::RectangularPath, height, thickness) =
    backend_wall_path(b, convert(OpenPolygonalPath, path), height, thickness)

#=
realize(b::Radiance, w::Wall) =
  let w_base_height = w.bottom_level.height,
      w_height = w.top_level.height - w_base_height,
      r_thickness = r_thickness(w),
      l_thickness = l_thickness(w),
      w_path = translate(w.path, vz(w_base_height)),
      w_paths = subpaths(w_path),
      r_w_paths = subpaths(offset(w_path, r_thickness)),
      l_w_paths = subpaths(offset(w_path, l_thickness)),
      openings = [w.doors..., w.windows...],
      prevlength = 0
    for (w_seg_path, r_w_path, l_w_path) in zip(w_paths, r_w_paths, l_w_paths)
      let currlength = prevlength + path_length(w_seg_path),
          c_r_w_path = closed_path_for_height(r_w_path, w_height),
          c_l_w_path = closed_path_for_height(l_w_path, w_height)
        realize_pyramid_fustrum(b, w, "wall", c_l_w_path, c_r_w_path, false)
        openings = filter(openings) do op
          if prevlength <= op.loc.x < currlength ||
             prevlength <= op.loc.x + op.family.width <= currlength # contained (at least, partially)
            let op_height = op.family.height,
                op_at_start = op.loc.x <= prevlength,
                op_at_end = op.loc.x + op.family.width >= currlength,
                op_path = subpath(w_path,
                                  max(prevlength, op.loc.x),
                                  min(currlength, op.loc.x + op.family.width)),
                r_op_path = offset(op_path, r_thickness),
                l_op_path = offset(op_path, l_thickness),
                fixed_r_op_path =
                  open_polygonal_path([path_start(op_at_start ? r_w_path : r_op_path),
                                       path_end(op_at_end ? r_w_path : r_op_path)]),
                fixed_l_op_path =
                  open_polygonal_path([path_start(op_at_start ? l_w_path : l_op_path),
                                       path_end(op_at_end ? l_w_path : l_op_path)]),
                c_r_op_path = closed_path_for_height(translate(fixed_r_op_path, vz(op.loc.y)), op_height),
                c_l_op_path = closed_path_for_height(translate(fixed_l_op_path, vz(op.loc.y)), op_height),
                idxs = closest_vertices_indexes(path_vertices(c_r_w_path), path_vertices(c_r_op_path))
              realize_pyramid_fustrum(b, w, "wall", c_r_op_path, c_l_op_path, false)
              c_r_w_path =
                closed_polygonal_path(
                  inject_polygon_vertices_at_indexes(path_vertices(c_r_w_path), path_vertices(c_r_op_path), idxs))
              c_l_w_path =
                closed_polygonal_path(
                  inject_polygon_vertices_at_indexes(path_vertices(c_l_w_path), path_vertices(c_l_op_path), idxs))
              # preserve if not totally contained
              ! (op.loc.x >= prevlength && op.loc.x + op.family.width <= currlength)
            end
          else
            true
          end
        end
        prevlength = currlength
        realize_polygon(b, w, "wall", c_l_w_path, false)
        realize_polygon(b, w, "wall", c_r_w_path, true)
      end
    end
    void_ref(b)
  end
=#

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

backend_curtain_wall(b::Unity, s, path::Path, bottom::Real, height::Real, thickness::Real, kind::Symbol) =
  backend_wall(b, translate(path, vz(bottom)), height, thickness, getproperty(s.family, kind))

############################################
#=
backend_bounding_box(b::Unity, shapes::Shapes) =
  UnityBoundingBox(connection(b), collect_ref(shapes))
=#

unity"public void SetView(Vector3 position, Vector3 target, float lens)"
unity"public Vector3 ViewCamera()"
unity"public Vector3 ViewTarget()"
unity"public float ViewLens()"

set_view(camera::Loc, target::Loc, lens::Real, b::Unity) =
  let c = connection(b)
    UnitySetView(c, camera, target, lens)
    interrupt_processing(c)
  end

get_view(b::Unity) =
  let c = connection(b)
    UnityViewCamera(c), UnityViewTarget(c), UnityViewLens(c)
  end

zoom_extents(b::Unity) = UnityZoomExtents(connection(b))

view_top(b::Unity) = UnityViewTop(connection(b))

unity"public void DeleteAll()"
unity"public void DeleteMany(GameObject[] objs)"
delete_all_shapes(b::Unity) = UnityDeleteAll(connection(b))

backend_delete_shapes(b::Unity, shapes::Shapes) =
  UnityDeleteMany(connection(b), 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()) =
    UnityCreateAlignedDimension(connection(b), 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
unity"public GameObject CreateParent(String name)"
unity"public GameObject CurrentParent()"
unity"public GameObject SetCurrentParent(GameObject newParent)"
unity"public void SetActive(GameObject obj, bool state)"
unity"public void DeleteAllInParent(GameObject parent)"
unity"public void SwitchToParent(GameObject newParent)"

UnityLayer = Int

current_layer(b::Unity)::UnityLayer =
  UnityCurrentParent(connection(b))

current_layer(layer::UnityLayer, b::Unity) =
  UnitySetCurrentParent(connection(b), layer)

create_layer(name::String, b::Unity) =
  UnityCreateParent(connection(b), name)

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

delete_all_shapes_in_layer(layer::UnityLayer, b::Unity) =
  UnityDeleteAllInParent(connection(b), layer)

switch_to_layer(layer::UnityLayer, b::Unity) =
  UnitySwitchToParent(connection(b), layer)

# Experiment to speed up things

canonicalize_layer(layer::UnityLayer, b::Unity) =
  UnityCanonicalize(connection(b), layer)

# Materials

UnityMaterial = Int

unity"public Material LoadMaterial(String name)"
unity"public void SetCurrentMaterial(Material material)"
unity"public Material CurrentMaterial()"

current_material(b::Unity)::UnityMaterial =
  UnityCurrentMaterial(connection(b))

current_material(material::UnityMaterial, b::Unity) =
  UnitySetCurrentMaterial(connection(b), material)

get_material(name::String, b::Unity) =
  UnityLoadMaterial(connection(b), name)

# Blocks

unity"public GameObject CreateBlockInstance(GameObject block, Vector3 position, Vector3 vx, Vector3 vy, float scale)"
unity"public GameObject CreateBlockFromShapes(String name, GameObject[] objs)"

realize(b::Unity, s::Block) =
  s.shapes == [] ?
    UnityLoadResource(connection(b), s.name) :
    UnityCanonicalize(connection(b), UnityCreateBlockFromShapes(connection(b), s.name, collect_ref(s.shapes)))

realize(b::Unity, s::BlockInstance) =
    UnityCreateBlockInstance(
        connection(b),
        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
unity"public GameObject PointLight(Vector3 position, Color color, float range, float intensity)"

backend_pointlight(b::Unity, loc::Loc, color::RGB, range::Real, intensity::Real) =
    UnityPointLight(connection(b), loc, color, range, intensity)
#=
backend_spotlight(b::Unity, loc::Loc, dir::Vec, hotspot::Real, falloff::Real) =
    UnitySpotLight(connection(b), loc, hotspot, falloff, loc + dir)

backend_ieslight(b::Unity, file::String, loc::Loc, dir::Vec, alpha::Real, beta::Real, gamma::Real) =
    UnityIESLight(connection(b), 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 c = connection(b)
          unknown(r, backend=b, ref=LazyRef(b, UnityNativeRef(r), 0, 0))
          #code = UnityShapeCode(c, r),
          #ref = LazyRef(b, UnityNativeRef(r))
          #error("Unknown shape with code $(code)")
      end
    else
      collected_shapes()[idx]
    end
  end
#
#=
Unity"public Point3d[] GetPosition(string prompt)"

select_position(prompt::String, b::Unity) =
  begin
    @info "$(prompt) on the $(b) backend."
    let ans = UnityGetPosition(connection(b), 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

Unity"public ObjectId[] GetPoint(string prompt)"

# 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)

Unity"public ObjectId[] GetCurve(string prompt)"

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

Unity"public ObjectId[] GetSurface(string prompt)"

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

Unity"public ObjectId[] GetSolid(string prompt)"

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

Unity"public ObjectId[] GetShape(string prompt)"

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

Unity"public long GetHandleFromShape(Entity e)"
Unity"public ObjectId GetShapeFromHandle(long h)"

captured_shape(b::Unity, handle) =
  shape_from_ref(UnityGetShapeFromHandle(connection(b), handle),
                 b)

generate_captured_shape(s::Shape, b::Unity) =
    println("captured_shape(autocad, $(UnityGetHandleFromShape(connection(b), ref(s).value)))")

# Register for notification

Unity"public void RegisterForChanges(ObjectId id)"
Unity"public void UnregisterForChanges(ObjectId id)"
Unity"public ObjectId[] ChangedShape()"
Unity"public void DetectCancel()"
Unity"public void UndetectCancel()"
Unity"public bool WasCanceled()"

register_for_changes(s::Shape, b::Unity) =
    let conn = connection(b)
        UnityRegisterForChanges(conn, ref(s).value)
        UnityDetectCancel(conn)
        s
    end

unregister_for_changes(s::Shape, b::Unity) =
    let conn = connection(b)
        UnityUnregisterForChanges(conn, ref(s).value)
        UnityUndetectCancel(conn)
        s
    end

waiting_for_changes(s::Shape, b::Unity) =
    ! UnityWasCanceled(connection(b))

changed_shape(ss::Shapes, b::Unity) =
    let conn = connection(b)
        changed = []
        while length(changed) == 0 && ! UnityWasCanceled(conn)
            changed =  UnityChangedShape(conn)
            sleep(0.1)
        end
        if length(changed) > 0
            shape_from_ref(changed[1], b)
        else
            nothing
        end
    end

Unity"public ObjectId[] GetAllShapes()"
Unity"public ObjectId[] GetAllShapesInLayer(ObjectId layerId)"

# HACK: This should be filtered on the plugin, not here.
all_shapes(b::Unity) =
    let c = connection(b)
        Shape[shape_from_ref(r, b)
              for r in filter(r -> UnityShapeCode(c, r) != 0, UnityGetAllShapes(c))]
    end

all_shapes_in_layer(layer, b::Unity) =
    let c = connection(b)
        Shape[shape_from_ref(r, b) for r in UnityGetAllShapesInLayer(c, layer)]
    end

disable_update(b::Unity) =
    UnityDisableUpdate(connection(b))

enable_update(b::Unity) =
    UnityEnableUpdate(connection(b))
# Render

=#
unity"public void SetResolution(int width, int height)"
unity"public void ScreenShot(String path)"

#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)
      UnitySetResolution(c, render_width(), render_height())
      interrupt_processing(c)
      UnityScreenShot(c, path)
    end

unity"public void SelectGameObjects(GameObject[] objs)"

highlight_shape(s::Shape, b::Unity) =
    UnitySelectGameObjects(connection(b), collect_ref(s))

highlight_shapes(ss::Shapes, b::Unity) =
    UnitySelectGameObjects(connection(b), collect_ref(ss))

unity"public void StartSelectingGameObject()"
unity"public int SelectedGameObjectId(bool existing)"

select_shape(prompt::String, b::Unity) =
  select_one_with_prompt(prompt, b, (c, prompt) ->
    let s = -2 # Means not found
      UnityStartSelectingGameObject(c)
      while s == -2
        sleep(0.1)
        s = UnitySelectedGameObjectId(c, true)
      end
      [s]
    end)


text:
requiredTypes = [ "AbstractLevel" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include(  reqType * ".jl" )
    end  # if !isdefined( Symbol( ...
end

export Assignment

type Assignment <: AbstractLevel
    # Name of the Assignment ex: VR, CP
    Name::String

    # Maximum and minimum stay in this Assignment
    MinStay::Int
    MaxStay::Int

    # A list of allowed Affiliations
    NextLevel::Array{String}

    Assignment() = new("", 0, 0, Array{String}())
    Assignment(nm::String;
                min::Int = 0,
                max::Int = 0,
                next::Array{String} = Array{String}()) = new(nm, min, max, next)
end

ext: jl

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
header_tmpl = "
<!DOCTYPE html>

<title>{{title}}</title>
<link href=\"./bootstrap/css/bootstrap.min.css\" rel=\"stylesheet\" media=\"screen\">

"

footer_tmpl = "

"

index_tmpl = "$header_tmpl

{{title}}

Name

Summary

{{#d}}

{{title}}

{{summary}}

{{/d}}
$footer_tmpl"

post_tmpl = "$header_tmpl

{{title}}

{{date}}

{{summary}}

{{content}}
$footer_tmpl"


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ d2a71096-feb7-11ea-0738-c170ec8f0f5a
using Pkg, DrWatson

# ╔═╡ d2a74dae-feb7-11ea-1c23-e7050ee80803
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 34c9c198-feb7-11ea-3b09-e5fa493d9fae
md"## Clip-06-06-09s.jl"

# ╔═╡ d2a7e3ea-feb7-11ea-1722-456d20e56b38
begin
	N = 100
	df = DataFrame(
	  :h0 => rand(Normal(10,2 ), N),
	  :treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2)))
	);
	df[!, :fungus] = [rand(Binomial(1, 0.5 - 0.4 * df[i, :treatment]), 1)[1] for i in 1:N]
	df[!, :h1] = [df[i, :h0] + rand(Normal(5 - 3 * df[i, :fungus]), 1)[1] for i in 1:N]
end

# ╔═╡ d2b6d230-feb7-11ea-13de-4fc9a25244fa
stan6_7 = "
data {
  int <lower=1> N;
  vector[N] h0;
  vector[N] h1;
  vector[N] treatment;
  vector[N] fungus;
}
parameters{
  real a;
  real bt;
  real bf;
  real<lower=0> sigma;
}
model {
  vector[N] mu;
  vector[N] p;
  a ~ lognormal(0, 0.2);
  bt ~ normal(0, 0.5);
  bf ~ normal(0, 0.5);
  sigma ~ exponential(1);
  for ( i in 1:N ) {
    p[i] = a + bt*treatment[i] + bf*fungus[i];
    mu[i] = h0[i] * p[i];
  }
  h1 ~ normal(mu, sigma);
}
";

# ╔═╡ d2b75bc2-feb7-11ea-1a1c-19917ffb69d4
begin
	m6_7s = SampleModel("m6.7s", stan6_7)
	m6_7_data = Dict(
	  :N => nrow(df),
	  :h0 => df[:, :h0],
	  :h1 => df[:, :h1],
	  :fungus => df[:, :fungus],
	  :treatment => df[:, :treatment]
	)
	rc6_7s = stan_sample(m6_7s; data=m6_7_data)
	success(rc6_7s) && (post6_7s_df = read_samples(m6_7s, :dataframe))
end;

# ╔═╡ d2c2b71c-feb7-11ea-124d-6114c352b17b
success(rc6_7s) && (part6_7s = Particles(post6_7s_df))

# ╔═╡ d2c34824-feb7-11ea-3c22-f91f4df2b0f5
success(rc6_7s) && (Text(precis(post6_7s_df; io=String)))

# ╔═╡ d2cebcf4-feb7-11ea-2d6a-798f2e408ddf
if success(rc6_7s)
	(s1, p1) = plot_model_coef([m6_7s], [:a, :bt, :bf])
	p1
end

# ╔═╡ 36092dde-00ec-11eb-2149-ad44ba127c89
s1

# ╔═╡ d2cf9264-feb7-11ea-08c7-b1bf1379fa13
md"## End of clip-06-06-09s.jl"

# ╔═╡ Cell order:
# ╠═34c9c198-feb7-11ea-3b09-e5fa493d9fae
# ╠═d2a71096-feb7-11ea-0738-c170ec8f0f5a
# ╠═d2a74dae-feb7-11ea-1c23-e7050ee80803
# ╠═d2a7e3ea-feb7-11ea-1722-456d20e56b38
# ╠═d2b6d230-feb7-11ea-13de-4fc9a25244fa
# ╠═d2b75bc2-feb7-11ea-1a1c-19917ffb69d4
# ╠═d2c2b71c-feb7-11ea-124d-6114c352b17b
# ╠═d2c34824-feb7-11ea-3c22-f91f4df2b0f5
# ╠═d2cebcf4-feb7-11ea-2d6a-798f2e408ddf
# ╠═36092dde-00ec-11eb-2149-ad44ba127c89
# ╟─d2cf9264-feb7-11ea-08c7-b1bf1379fa13


text:
using StanOptimize, Test

if haskey(ENV, "JULIA_CMDSTAN_HOME")

    @testset "Bernoulli optimize" begin
      include(joinpath(@__DIR__, "../examples/Bernoulli/bernoulli.jl"))
      @test optim1["theta"][end] ≈ 0.3 atol=0.1
      @test optim2["theta"][end] ≈ 0.3 atol=0.1
    end

else
  println("\nJULIA_CMDSTAN_HOME not set. Skipping tests")
end


text:
function featurize(df)
  X = zeros(size(df, 1), 0)
  function onehot(column::DataArray)
    categories = unique(column)
    p = length(categories)
    n = size(column, 1)
    result = zeros(n, p)
    for (i, x) in enumerate(column)
      result[i,:] = eye(p)[:,findfirst(categories, x)]
    end
    return result
  end

  for feature in names(df)
    column = df[feature]
    if eltype(column) <: Number
      X = [X (column - mean(column))/std(column)]
    else
      X = [X onehot(column)]
    end
  end
  return X
end

function class(column)
  result = zeros(length(column))
  categories = unique(column)
  @assert length(categories) == 2
  for (i, x) in enumerate(column)
    if findfirst(categories, x) == 1
      result[i] = -1
    else
      @assert findfirst(categories, x) == 2
      result[i] = 1
    end
  end
  return result
end

max_stars_repo_name: pcmoritz/slbfgs

max_issues_repo_name: pcmoritz/slbfgs

max_forks_repo_name: pcmoritz/slbfgs


text:
#=
    Used from https://github.com/lanl-ansi/Juniper.jl
=#

"""
MOI functions, sets and, other type definitions
"""
# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# sets
const BOUNDS = Union{
    MOI.EqualTo{Float64},
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
}

const VAR_TYPES = Union{
    MOI.ZeroOne,
    MOI.Integer
}

# other MOI types
const AFF_TERM = MOI.ScalarAffineTerm{Float64}
const QUAD_TERM = MOI.ScalarQuadraticTerm{Float64}

"""
Variable information struct definition
"""
mutable struct VariableInfo
    lower_bound::Float64  # May be -Inf even if has_lower_bound == true
    has_lower_bound::Bool # false implies lower_bound == -Inf
    upper_bound::Float64  # May be Inf even if has_upper_bound == true
    has_upper_bound::Bool # false implies upper_bound == Inf
    is_fixed::Bool        # Implies lower_bound == upper_bound and !has_lower_bound and !has_upper_bound
    is_binary::Bool       # Implies lower_bound == 0, upper_bound == 1 and is MOI.ZeroOne
    is_integer::Bool      # Implies variable is MOI.Integer
    start::Real           # Primal start
    name::String
end
VariableInfo() = VariableInfo(-Inf, false, Inf, false, false, false, false, 0.0, "")

"""
Optimizer struct
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{OriginalProblem, Nothing}
    variable_info::Vector{VariableInfo}
    nlp_data::MOI.NLPBlockData
    sense::MOI.OptimizationSense
    objective::Union{SVF, SAF, SQF, Nothing}
    linear_le_constraints::Vector{Tuple{SAF, MOI.LessThan{Float64}}}
    linear_ge_constraints::Vector{Tuple{SAF, MOI.GreaterThan{Float64}}}
    linear_eq_constraints::Vector{Tuple{SAF, MOI.EqualTo{Float64}}}
    quadratic_le_constraints::Vector{Tuple{SQF, MOI.LessThan{Float64}}}
    quadratic_ge_constraints::Vector{Tuple{SQF, MOI.GreaterThan{Float64}}}
    quadratic_eq_constraints::Vector{Tuple{SQF, MOI.EqualTo{Float64}}}
    #one field needs to be added to keep track of nlp constraint types
    options::SolverOptions
end

MOI.get(::Optimizer, ::MOI.SolverName) = "DuranSolver"

MOI.supports(::Optimizer, ::MOI.Silent) = true
MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = true
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::Optimizer, ::MOI.RawParameter) = true

function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    if value
        model.options.log_levels = []
    end
    model.options.silent = value
    return
end

function MOI.set(model::Optimizer, ::MOI.NumberOfThreads, value::Union{Nothing,Int})
    if value === nothing
        model.options.processors = 1
    else
        model.options.processors = value
    end
    return
end

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Nothing,Float64})
    if value === nothing
        model.options.time_limit = Inf
    else
        model.options.time_limit = value
    end
    return
end
#this function needs to be updated
function MOI.set(model::Optimizer, p::MOI.RawParameter, value)
    p_symbol = Symbol(p.name)
    if in(p_symbol, fieldnames(SolverOptions))
        type_of_param = fieldtype(SolverOptions, p_symbol)
        if hasmethod(convert, (Type{type_of_param}, typeof(value)))
            passed_checks = true
            if p_symbol == :traverse_strategy
                if !(value in [:BFS, :DFS, :DBFS])
                    passed_checks = false
                    @error "Traverse strategy $(value) is not supported. Use one of `[:BFS, :DFS, :DBFS]`."
                end
            end
            if p_symbol == :branch_strategy
                if !(value in [:StrongPseudoCost, :PseudoCost, :Reliability, :MostInfeasible])
                    passed_checks = false
                    @error "Branch strategy $(value) is not supported. Use one of `[:StrongPseudoCost, :PseudoCost, :Reliability, :MostInfeasible]`."
                end
            end
            passed_checks && setfield!(model.options, p_symbol, convert(type_of_param, value))
        else
            @error "The option $(p.name) has a different type ($(type_of_param))"
        end
    else
        @error "The option $(p.name) doesn't exist."
    end
    return
end

# Returns the number of solutions that can be retrieved
MOI.get(model::Optimizer, ::MOI.ResultCount) = length(model.inner.solutions)

MOI.get(model::Optimizer, ::MOI.NumberOfThreads) = model.options.processors

MOI.get(model::Optimizer, ::MOI.TimeLimitSec) = model.options.time_limit

MOI.get(model::Optimizer, ::MOI.Silent) = model.options.silent

function MOI.get(model::Optimizer, p::MOI.RawParameter)
    if in(p.name, fieldnames(SolverOptions))
        return getfield(model.options, p.name)
    end
    @error "The option $(p.name) doesn't exist."
end

"""
EmptyNLPEvaluator struct and associated functions
"""
struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    @assert length(g) == 0
    return
end
function MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x)
    fill!(g, 0.0)
    return
end
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    @assert length(J) == 0
    return
end
function MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, σ, μ)
    @assert length(H) == 0
    return
end
empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

function register(s::Symbol, dimension::Integer, f::Function; autodiff::Bool=false)
    return RegisteredFunction(s, dimension, f, nothing, nothing, autodiff)
end

function register(s::Symbol, dimension::Integer, f::Function, gradf::Function; autodiff::Bool=false)
    return RegisteredFunction(s, dimension, f, gradf, nothing, autodiff)
end

function register(s::Symbol, dimension::Integer, f::Function, gradf::Function, grad2f::Function; autodiff::Bool=false)
    return RegisteredFunction(s, dimension, f, gradf, grad2f, autodiff)
end

"""
Optimizer struct constructor
"""
function Optimizer(;options...)

    solver_options = combine_options(options)

    return Optimizer(
    nothing,
    [],
    empty_nlp_data(),
    MOI.FEASIBILITY_SENSE,
    nothing,
    [], [], [], # linear constraints
    [], [], [], # quadratic constraints
    solver_options)
end

function Optimizer(options::Vector{Pair{String,Any}})
    symbol_options = Dict{Symbol, Any}()
    for option in options
        symbol_options[Symbol(option.first)] = option.second
    end
    return Optimizer(symbol_options)
end

function Optimizer(options::Dict{Symbol,Any})

    solver_options = combine_options(options)

    return Optimizer(
    nothing,
    [],
    empty_nlp_data(),
    MOI.FEASIBILITY_SENSE,
    nothing,
    [], [], [], # linear constraints
    [], [], [], # quadratic constraints
    solver_options)
end

"""
Printing the optimizer
"""
function Base.show(io::IO, model::Optimizer)
    println("A HybridMINLPSolver MathOptInterface model with backend")
    return
end

"""
Copy constructor for the optimizer
"""
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = true
function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; kws...)
    return MOI.Utilities.automatic_copy_to(model, src; kws...)
end

"""
``MOI.is_empty(model::Optimizer)`` overload for Alpine.Optimizer
"""
function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) &&
        model.nlp_data.evaluator isa EmptyNLPEvaluator &&
        model.sense == MOI.FEASIBILITY_SENSE &&
        isempty(model.linear_le_constraints) &&
        isempty(model.linear_ge_constraints) &&
        isempty(model.linear_eq_constraints) &&
        isempty(model.quadratic_le_constraints) &&
        isempty(model.quadratic_ge_constraints) &&
        isempty(model.quadratic_eq_constraints)
end

"""
``MOI.empty!(model::Optimizer)`` overload for Alpine.Optimizer
"""
function MOI.empty!(model::Optimizer)
    model.inner = nothing
    empty!(model.variable_info)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_ge_constraints)
    empty!(model.linear_eq_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_ge_constraints)
    empty!(model.quadratic_eq_constraints)
end

"""
ordering of constraints provided to HybridMINLPSolver.jl
"""
linear_le_offset(model::Optimizer) = 0
linear_ge_offset(model::Optimizer) = length(model.linear_le_constraints)
linear_eq_offset(model::Optimizer) = linear_ge_offset(model) + length(model.linear_ge_constraints)
quadratic_le_offset(model::Optimizer) = linear_eq_offset(model) + length(model.linear_eq_constraints)
quadratic_ge_offset(model::Optimizer) = quadratic_le_offset(model) + length(model.quadratic_le_constraints)
quadratic_eq_offset(model::Optimizer) = quadratic_ge_offset(model) + length(model.quadratic_ge_constraints)
nlp_constraint_offset(model::Optimizer) = quadratic_eq_offset(model) + length(model.quadratic_eq_constraints)

function info_array_of_variables(variable_info::Vector{VariableInfo}, attr::Symbol)
    len_var_info = length(variable_info)
    type_dict = get_type_dict(variable_info[1])
    result = Array{type_dict[attr], 1}(undef, len_var_info)
    for i = 1:len_var_info
        result[i] = getfield(variable_info[i], attr)
        # if type is binary then set bounds correctly
        if result[i] < 0 && attr == :lower_bound && getfield(variable_info[i], :is_binary)
            result[i] = 0
        end
        if result[i] > 1 && attr == :upper_bound && getfield(variable_info[i], :is_binary)
            result[i] = 1
        end
    end
    return result
end

"""
``MOI.optimize!()`` for Duran Solver
"""
function MOI.optimize!(model::Optimizer)
    MOI.initialize(model.nlp_data.evaluator, [:Grad, :Jac, :ExprGraph])
    if ~isa(model.nlp_data.evaluator, EmptyNLPEvaluator)
    else
        @info "no explicit NLP constraints or objective provided using @NLconstraint or @NLobjective macros"
    end
    model.inner = OriginalProblem()
    @views op = model.inner
    #initiate to replicate original problem
    init_original_problem!(op, model)
    #complete original problem
    create_original_model!(model, op)
    print("This is the original problem \n", op.model)
    #perform outer approximation algorithm
    oa_algorithm(model, jp)
end

getnsolutions(m::OriginalProblem) = m.nsolutions
getsolutions(m::OriginalProblem) = m.solutions

include("moi_wrapper.jl")

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

using UUIDs

export zeropad
function zeropad(x::Vector{Matrix{T}}) where {T}
  Nt_max = maximum(size.(x,1))
  Nf = size(x[1],2)
  Nb = length(x)
  x_batched = zeros(T,Nt_max,Nf,Nb)
  for (i,xi) in enumerate(x)
    for f = 1:Nf, t = 1:size(xi,1)
      x_batched[t,f,i] = xi[t,f]
    end
  end
  return x_batched
end

export posterior2phones
function posterior2phones(ippsym, gamma)
  #z = [replace(ippsym[argmax(gamma[t,:])], r"[0-9\-]"=>"") for t in 1:size(gamma,1)]
  p = [ippsym[argmax(gamma[t,:])] for t in 1:size(gamma,1)] # all phones symbols
  k = [p[1]] # phones without repetitions
  for t in 2:length(p)
    if p[t-1] != p[t]
      push!(k,p[t])
    end
  end
  k2 = String[]
  for ki in k
    if ki[end] == '1' # emitting symbol
      push!(k2,replace(ki, r"[0-9\-]"=>""))
    end
  end
  return k2
end

export levenshtein
function levenshtein(s,t)
  n,m = length(s),length(t)
  D = zeros(Int,n+1,m+1)

  D[:,1] = 0:n # this is the cost we would have for insertion only
  D[1,:] = 0:m # this is the cost we would have for deletion only
  for i = 2:n+1, j = 2:m+1
    # check substition is needed
    cost = s[i-1] == t[j-1] ? 0 : 1

    D[i,j] = min(
                 D[i-1,j] + 1,        # del
                 D[i,j-1] + 1,        # ins
                 D[i-1,j-1] + cost,   # subs / ok
                )
  end
  return D[n+1,m+1]
end

export text2phones
function text2phones(lexicon,text; add_sil=true)
  if add_sil
    phones = [[lexicon[t]...,"<SIL>"] for t in split(text)]
  else
    phones = [lexicon[t] for t in split(text)]
  end
  phones = vcat(phones...)
  if add_sil
    phones = ["<SIL>",phones...]
  end
  return phones
end

export get_error_rate
function get_error_rate(uttID2seq::Dict,
    uttID2seq_dec::Dict; kwargs...)
  seqs, seq_decs = [], []
  for uttID in keys(uttID2seq)
    push!(seqs,uttID2seq[uttID])
    push!(seq_decs,uttID2seq_dec[uttID])
  end
  get_error_rate(seqs,seq_decs; kwargs...)
end

function get_error_rate(seqs::Vector,seq_decs::Vector; is_split=false)
  Nw = 0
  err = 0
  for i in eachindex(seqs)
    seq, seq_dec = seqs[i], seq_decs[i]
    if is_split == false
      seq     = split(seq;keepempty=false)
      seq_dec = split(seq_dec; keepempty=false)
    end
    Nw += length(seq)
    err += levenshtein(seq,seq_dec)
  end
  er = err/Nw
end

export min_dist_word
function min_dist_word(prons,min_dist,word_phones)
  d = [levenshtein(word_phones,pr) for pr in prons]
  idxs = findall(d .<= min_dist)
  if isempty(idxs)
    return "<?>"
  else
    return ilexicon[prons[idxs[argmin(d[idxs])]]]
  end
end

function add_word!(dec,prons,word_phones; min_dist=2)
  t = try
    ilexicon[word_phones]
  catch
    if min_dist == 0
      "<?>"
    else
      min_dist_word(prons,min_dist,word_phones)
    end
  end
  push!(dec,t)
end

function phones2words_greedy(ilexicon,phones; min_dist=2)
  prons = [keys(ilexicon)...]
  word_phones = String[]
  dec = String[]
  for (i,p) in enumerate(phones)
    if i == 1
      word_phones = String[]
      if p != "<SIL>"
        push!(word_phones,p)
      end
    elseif (p == "<SIL>") && i > 1
      if !isempty(word_phones)
        add_word!(dec,prons,word_phones; min_dist=min_dist)
        word_phones = String[]
      end
    else
      if p != "<SIL>"
        push!(word_phones,p)
      end
    end
  end
  if !isempty(word_phones)
    add_word!(dec,prons,word_phones; min_dist=min_dist)
  end
  return dec
end

export check_env
function check_env()
  if !("TIDIGITS_PATH" in keys(ENV))
    @warn "ENV[\"TIDIGITS_PATH\"] not exisitng: `export TIDIGITS_PATH=path/to/dataset` to your env."
  end
  if !("CPU_CMD" in keys(ENV))
    @warn "ENV[\"CPU_CMD\"] not exisitng: `export CPU_CMD='...'` to your env first. Only needed for SGE."
  end
end

export get_feat_dir
function get_feat_dir(setup; root="data")
  uuid_folder = UUID("04a07b93-95e4-4b85-94b9-d3516eb06ea2")
  conf = read("conf/$(setup)/feat_conf.jl",String)
  return joinpath("data", string(uuid5(uuid_folder,conf)))
end


text:
# Autogenerated wrapper script for IntelOpenMP_jll for i686-w64-mingw32
export libiomp

JLLWrappers.@generate_wrapper_header("IntelOpenMP")
JLLWrappers.@declare_library_product(libiomp, "libiomp5md.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libiomp,
        "bin\\libiomp5md.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## Dependencies

#To install Laplacians you need to put in Julia's Prompt "add Laplacians#master"
#For some reason Pkg.add() doesn't work with this.
#Pkg.add("MATLAB")
#using Laplacians
#using Pkg
#include(string(Pkg.dir("Laplacians") , "/src/matlabSolvers.jl"))
getlagged(x) = [NaN; x[1:(end - 1)]]

#Defining types and structures
abstract type AbstractLLSAlgorithm end
abstract type AbstractGraphLLSAlgorithm <: AbstractLLSAlgorithm end  # i.e. requires graph laplacian
struct CMGPreconditionedLLS <: AbstractGraphLLSAlgorithm  end
struct AMGPreconditionedLLS <: AbstractGraphLLSAlgorithm  end
struct DirectLLS <: AbstractLLSAlgorithm  end  # no graph required

abstract type AbstractLeverageAlgorithm end
struct ExactAlgorithm <: AbstractLeverageAlgorithm  end
@with_kw struct JLAAlgorithm <: AbstractLeverageAlgorithm
    num_simulations::Int64 = 0
end

@with_kw struct Settings{LeverageAlgorithm}
    cg_maxiter::Int64 = 300
    leverage_algorithm::LeverageAlgorithm = ExactAlgorithm()
    #clustering_level::String = "obs"
    person_effects::Bool = false
    cov_effects::Bool = false
end

## Methods and Types Definitions

#1) Finds AKM largest connected set
function find_connected_set(y, idvar, firmidvar; verbose=false)

    #Observation identifier to join later the FE
    obs_id = collect(1:size(y,1))

    firms = unique(sort(firmidvar))
    firmid = indexin(firmidvar, firms)

    workers = unique(sort(idvar))
    id = indexin(idvar, workers)

    #Save ids just in case
    firmid_old=firmidvar;
    id_old=idvar;

    n_firms=length(firms)
    n_workers=length(workers)
    firmid=firmid.+n_workers

    graph_size=length(firms)+length(workers)
    G=Graph(graph_size)
    for i in 1:size(y,1)
        add_edge!(G, firmid[i], id[i])
    end

    cs=connected_components(G)
    if verbose == true
        println("Largest connected set has been found")
    end
    pos = indexin(  [  maximum(size.(cs,1))] , size.(cs,1) )[1]
    lcs=cs[pos]

    connected_firms=lcs[lcs.>n_workers]

    sel=findall(in(connected_firms),firmid)

    obs_id = [obs_id[x] for x in sel ]
    yvec = [y[x] for x in sel ]
    firmid = [firmidvar[x] for x in sel ]
    id = [idvar[x] for x in sel ]

    #Relabel firms
    firms = unique(sort(firmid))
    firmid = indexin(firmid, firms)

    #Relabel workers
    ids = unique(sort(id))
    id = indexin(id, ids)

    return (obs_id = obs_id , y = yvec , id = id, firmid = firmid )

end

#2) Pruning and finding Leave-Out Largest connected set
function prunning_connected_set(yvec, idvar, firmidvar, obs_id; verbose=false)

    firms = unique(sort(firmidvar))
    firmid = indexin(firmidvar, firms)
    workers = unique(sort(idvar))
    id = indexin(idvar, workers)

    nbadworkers=1
    while nbadworkers>0

        n_firms=length(firms)
        n_workers=length(workers)
        firmid=firmid.+n_workers

        graph_size=n_workers + n_firms
        lcs_graph=Graph(graph_size)
            for i in 1:size(yvec,1)
                add_edge!(lcs_graph, firmid[i], id[i])
            end

        #get articulation vertex
        artic_vertex = articulation(lcs_graph)

        sel=findall(x->x==nothing, indexin(id,artic_vertex))

        bad_workers = artic_vertex[artic_vertex.<=n_workers]
        nbadworkers = size(bad_workers,1)

        if verbose==true
            println("Number of workers that are articulation points: ", nbadworkers)
        end

        #Restrict the sample
        yvec = [yvec[x] for x in sel ]
        firmid = [firmid[x] for x in sel ]
        id = [id[x] for x in sel ]
        obs_id = [obs_id[x] for x in sel ]

        #Relabel firms
        firms = unique(sort(firmid))
        firmid = indexin(firmid, firms)

        #Relabel workers
        ids = unique(sort(id))
        id = indexin(id, ids)

        n_workers=maximum(id)
        n_firms=maximum(firmid);

        firmid = firmid .+ n_workers

        #Constructing new Graph
        G=Graph(n_workers+n_firms)
        Nprimeprime = size(sel,1)
        for i in 1:Nprimeprime
            add_edge!(G, firmid[i], id[i])
        end

        #Find Largest Connected Set (LCS)
        cs=connected_components(G)

        pos = indexin(  [  maximum(size.(cs,1))] , size.(cs,1) )[1]
        lcs=cs[pos]

        connected_firms=lcs[lcs.>n_workers]

        sel=findall(in(connected_firms),firmid)

        obs_id = [obs_id[x] for x in sel ]
        yvec = [yvec[x] for x in sel ]
        firmid = [firmid[x] for x in sel ]
        id = [id[x] for x in sel ]

        #Relabel firms
        firms = unique(sort(firmid))
        firmid = indexin(firmid, firms)

        #Relabel workers
        ids = unique(sort(id))
        id = indexin(id, ids)

    end

    return (obs_id = obs_id , y = yvec , id = id, firmid = firmid )
end

#3) Drops Single Observations

function drop_single_obs(yvec, idvar, firmidvar,obs_id)

    firms = unique(sort(firmidvar))
    firmid = indexin(firmidvar, firms)
    workers = unique(sort(idvar))
    id = indexin(idvar, workers)

    T = Matlab.accumarray(id,1)
    T = [T[x] for x in id]
    sel = T.>1
    sel = findall(x->x==true,sel )

    obs_id = [obs_id[x] for x in sel ]
    yvec = [yvec[x] for x in sel ]
    firmid = [firmid[x] for x in sel ]
    id = [id[x] for x in sel ]

    #Relabel firms
    firms = unique(sort(firmid))
    firmid = indexin(firmid, firms)

    #Relabel workers
    ids = unique(sort(id))
    id = indexin(id, ids)

    return (obs_id = obs_id , y = yvec , id = id, firmid = firmid )
end

#4) Finds the observations connected for every cluster

function index_constr(clustering_var, id, match_id )
    NT = length(clustering_var)
    counter = ones(size(clustering_var,1));

    #Indexing obs number per worked/id
    gcs = Int.(@transform(groupby(DataFrame(counter = counter, id = id), :id), gcs = cumsum(:counter)).gcs);
    maxD = maximum(gcs);

    index = collect(1:NT);

    #This will be the output, and we will append observations to it in the loop
    list_final=DataFrame(row = Int64[],col = Int64[], match_id = Int64[], id_cluster = Int64[]);

    for t=1:maxD
    rowsel =  findall(x->x==true,gcs.==t);
    list_base = DataFrame( id_cluster= [clustering_var[x] for x in rowsel], row =
    [index[x] for x in rowsel] , match_id = [match_id[x] for x in rowsel]  );

        for tt=t:maxD
            colsel =  findall(x->x==true,gcs.==tt);
            list_sel =  DataFrame( id_cluster= [clustering_var[x] for x in colsel], col =
            [index[x] for x in colsel] );

            merge = outerjoin(list_base, list_sel, on = :id_cluster)
            merge = dropmissing(merge)
            merge = merge[:,[:row, :col, :match_id, :id_cluster]]

            append!(list_final, merge)

        end
    end

    sort!(list_final, (:row));

    return Matrix(list_final)

end

#5) Compute Movers
function compute_movers(id,firmid)

    gcs = [NaN; id[1:end-1]]
    gcs = id.!=gcs

    lagfirmid=[NaN; firmid[1:end-1]]
    for x in findall(x->x ==true , gcs)
        lagfirmid[x] = NaN
    end

    stayer=(firmid.==lagfirmid)
    for x in findall(x->x ==true , gcs)
        stayer[x] = true
    end

    stayer=Int.(stayer)
    stayer=Matlab.accumarray(id,stayer)
    T=Matlab.accumarray(id,1)
    stayer=T.==stayer
    movers=stayer.==false

    movers = [movers[x] for x in id]
    T = [T[x] for x in id]

    return (movers = movers, T = T)
end

#5) Finds observations connected for every cluster

function check_clustering(clustering_var)
    NT = length(clustering_var)
    counter = ones(size(clustering_var,1));

    #Indexing obs number per worked/id
    gcs = Int.(@transform(groupby(DataFrame(counter = counter, clustering_var = clustering_var), :clustering_var), gcs = cumsum(:counter)).gcs)
    maxD = maximum(gcs)

    index = collect(1:NT)

    #This will be the output, and we will append observations to it in the loop
    list_final=DataFrame(row = Int64[],col = Int64[], id_cluster = Int64[])

    for t=1:maxD
        rowsel =  findall(x->x==true,gcs.==t)
        list_base = DataFrame( id_cluster= [clustering_var[x] for x in rowsel], row =
        [index[x] for x in rowsel]  )

            for tt=t:maxD
                colsel =  findall(x->x==true,gcs.==tt)
                list_sel =  DataFrame( id_cluster= [clustering_var[x] for x in colsel], col =
                [index[x] for x in colsel] )

                merge = outerjoin(list_base, list_sel, on = :id_cluster)
                merge = dropmissing(merge)
                merge = merge[:,[:row, :col, :id_cluster]]

                append!(list_final, merge)

            end
        end

        sort!(list_final, (:row))
        list_final = Matrix(list_final)

        return (list_final = list_final , nnz_2 = size(list_final,1) )

end

# Assumes that obs_leaveoutset is an array of zeros
function whole_prunning!(obs_leaveoutset,y,id,firmid;verbose=false)
    obs,  y  , id , firmid  = find_connected_set(y,id,firmid,verbose=verbose)
    obs,  y  , id , firmid  = prunning_connected_set(y,id,firmid,obs,verbose=verbose)
    obs,  y  , id , firmid  = drop_single_obs(y,id,firmid,obs)
    for x in obs
        obs_leaveoutset[x] = 1
    end
    return nothing
end

#6) Eff res : Compute Effective Resistance - Lambda Matrices

function eff_res(::ExactAlgorithm, X,id,firmid,match_id, K, settings)

    #Indexing Observations
    elist = index_constr(collect(1:length(id)), id, match_id )

    #Dimensions
    NT = size(X,1)
    M = size(elist,1)
    J = maximum(firmid)
    N = maximum(id)

    #Define solver
    S_xx = X'*X

    # Create the solvers
    ldli, la = computeLDLinv(S_xx)
    buffs = zeros(size(la)[1],Threads.nthreads())
    compute_sol = []
    for i in 1:Threads.nthreads()
        P = approxcholOperator(ldli,buffs[:,i])
        push!(compute_sol,approxcholSolver(P,la))
    end

    #Initialize output
    Pii = zeros(M)
    Bii_fe = zeros(M)
    Bii_cov= settings.cov_effects ==true ? zeros(M) : nothing
    Bii_pe= settings.person_effects == true ? zeros(M) : nothing

    #No controls case: We compute Pii,Bii for stayers manually
    if K == 0

        #Compute Auxiliaries
        movers , T = compute_movers(id, firmid)

        Nmatches = maximum(match_id)
        match_id_movers = [match_id[x] for x in findall(x->x==true, movers)]
        firmid_movers = [firmid[x] for x in findall(x->x==true, movers)]
        id_movers = [id[x] for x in findall(x->x==true, movers)]

        sel = unique(z -> match_id_movers[z], 1:length(match_id_movers))
        match_id_movers = match_id_movers[sel]
        firmid_movers = firmid_movers[sel]
        id_movers = id_movers[sel]

        maxT = maximum([T[x] for x in findall(x->x == false, movers)])

        counter = ones(Int,NT)
        gcs = Int.(@transform(groupby(DataFrame(counter = counter, id = id), :id), gcs = cumsum(:counter)).gcs)
        sel_stayers = (gcs.==1).*(movers.==false)
        stayers_matches_sel = [match_id[z] for z in findall(x->x == true , sel_stayers)]
        Tinv = 1 ./T
        elist_JLL = [id_movers N.+firmid_movers id_movers N.+firmid_movers]

        M = size(elist_JLL,1)
        Pii_movers = zeros(M)
        Bii_fe_movers = zeros(M)
        Bii_cov_movers= settings.cov_effects ==true ? zeros(M) : nothing
        Bii_pe_movers= settings.person_effects == true ? zeros(M) : nothing

        #Initializing dependent variables for solver
        Xright = sparse(collect(1:M),elist_JLL[:,1],1.0,M,N+J)
        Xright = Xright .+ sparse(collect(1:M),elist_JLL[:,2],-1.0,M,N+J)
        # N+J x N+J-1 restriction matrix
        S= sparse(1.0I, J-1, J-1)
        S=vcat(S,sparse(-zeros(1,J-1)))

        Xright = hcat(Xright[:,1:N], Xright[:,N+1:end]*S)

        Threads.@threads for i=1:M

            #Only one inversion needed for exact alg
            zexact = compute_sol[Threads.threadid()]( [Xright[i,:]...] ; verbose=false)

            #Compute Pii
            Pii_movers[i] = Xright[i,:]'*zexact

            #Compute Bii for firms
            aux_right = zexact[N+1:N+J-1]
            aux_left = zexact[N+1:N+J-1]

            COV = cov(X[:,N+1:N+J-1]*aux_left,X[:,N+1:N+J-1]*aux_right)
            Bii_fe_movers[i] = COV[1]*(NT-1)

            if Bii_pe != nothing
                aux_right = zexact[1:N]
                aux_left = zexact[1:N]
                COV = cov(X[:,1:N]*aux_left,X[:,1:N]*aux_right)
                Bii_pe_movers[i] = COV[1]*(NT-1)
            end

            if Bii_cov != nothing
                aux_right = zexact[N+1:N+J-1]
                aux_left = zexact[1:N]
                COV = cov(X[:,1:N]*aux_left,X[:,N+1:N+J-1]*aux_right)
                Bii_cov_movers[i] = COV[1]*(NT-1)
            end

        end

        #Assign Step
        Pii_movers = sparse(match_id_movers,ones(Int,length(match_id_movers)),Pii_movers[:,1],Nmatches,1)
        Pii_stayers = sparse(stayers_matches_sel,ones(Int,length(stayers_matches_sel)),[Tinv[x] for x in findall(x->x==true,sel_stayers)],Nmatches,1)
        Pii = Pii_movers.+Pii_stayers

        Bii_fe = sparse(match_id_movers,ones(Int,length(match_id_movers)),Bii_fe_movers[:,1],Nmatches,1)

        if settings.cov_effects == true
            Bii_cov = sparse(match_id_movers,ones(Int,length(match_id_movers)),Bii_cov_movers[:,1],Nmatches,1)
        end

        if settings.person_effects == true
            Bii_pe = sparse(match_id_movers,ones(Int,length(match_id_movers)),Bii_pe_movers[:,1],Nmatches,1)
            stayers = .!movers

            Threads.@threads for t=2:maxT #T=1 have Pii=1 so need to be dropped.

                sel = (gcs.==true).*stayers.*(T.==t)
                N_sel = sum(sel)

                if N_sel > 0
                    index_sel = findall(x->x==true,sel)
                    match_sel_aux = Int.([match_id[z] for z in index_sel])
                    first = index_sel[1]
                    Xuse = X[first,:]

                    ztilde = compute_sol[Threads.threadid()]([Xuse...] ;verbose=false)

                    aux_right = ztilde[1:N]
                    aux_left = ztilde[1:N]

                    COV = cov(X[:,1:N]*aux_left,X[:,1:N]*aux_right)
                    Bii_pe_stayers = COV[1]*(NT-1)

                    Bii_pe_stayers = sparse(match_sel_aux,ones(Int,length(match_sel_aux)),Bii_pe_stayers,Nmatches,1)
                    Bii_pe = Bii_pe.+Bii_pe_stayers
                end
            end
        end

    #Controls case
    elseif K> 0
        #AUX = initialize_auxiliary_variables(settings.lls_algorithm, X, elist, M,NT, N, J, K, settings)
        nparameters = N + J + K

        D=sparse(collect(1:NT),id,1)
        F=sparse(collect(1:NT),firmid,1)
        S= sparse(1.0I, J-1, J-1)
        S=vcat(S,sparse(-zeros(1,J-1)))

        Dvar = hcat(  D, spzeros(NT,nparameters-N) )
        Fvar = hcat(spzeros(NT,N), -F*S, spzeros(NT,nparameters-N-J) )
        #Wvar = hcat(spzeros(NT,N+J), controls )
        Xleft = X[elist[:,1],:]
        Xright = X[elist[:,2],:]

        Threads.@threads for i=1:M

                #Again, one inversion needed
                zexact = compute_sol[Threads.threadid()]([Xright[i,:]...];verbose=false)

                #Compute Pii
                Pii[i] = Xleft[i,:]'*zexact

                #Compute Bii for firms
                aux_right = zexact[N+1:N+J-1,:]
                aux_left = zexact[N+1:N+J-1,:]

                COV = cov(X[:,N+1:N+J-1]*aux_left,X[:,N+1:N+J-1]*aux_right)
                Bii_fe[i] = COV[1]*(NT-1)

                if Bii_pe != nothing
                    aux_right = zexact[1:N]
                    aux_left = zexact[1:N]
                    COV = cov(X[:,1:N]*aux_left,X[:,1:N]*aux_right)
                    Bii_pe[i] = COV[1]*(NT-1)
                end

                if Bii_cov != nothing
                    aux_right = zexact[N+1:N+J-1]
                    aux_left = zexact[1:N]
                    COV = cov(X[:,1:N]*aux_left,X[:,N+1:N+J-1]*aux_right)
                    Bii_cov[i] = COV[1]*(NT-1)
                end

        end

    end

    #Create matrices
    rows = elist[:,1]
    cols = elist[:,2]
    index_cluster = match_id

    #Censor
    Pii[ findall(Pii.>=0.99)] .= 0.99

    if K==0
        Pii = [Pii[x] for x in index_cluster]
        Bii_fe = [Bii_fe[x] for x in index_cluster]

        if settings.cov_effects == true
            Bii_cov = [Bii_cov[x] for x in index_cluster]
        end

        if settings.person_effects == true
            Bii_pe = [Bii_pe[x] for x in index_cluster]
        end

    end

    #Lambda P
    Lambda_P=sparse(rows,cols,Pii,NT,NT)
    Lambda_P=Lambda_P+triu(Lambda_P,1)'
    println("Lambda P Computed!\n")

    #Lambda B var(fe)
    Lambda_B_fe=sparse(rows,cols,Bii_fe,NT,NT)
    Lambda_B_fe=Lambda_B_fe+triu(Lambda_B_fe,1)'
    println("Lambda B FE Computed!\n")

    #Lambda B cov(fe,pe)
    if settings.cov_effects == true
        Lambda_B_cov=sparse(rows,cols,Bii_cov,NT,NT)
        Lambda_B_cov=Lambda_B_cov+triu(Lambda_B_cov,1)'
    end

    #Lambda B, var(pe)
    if settings.person_effects == true
        Lambda_B_pe=sparse(rows,cols,Bii_pe,NT,NT)
        Lambda_B_pe=Lambda_B_pe+triu(Lambda_B_pe,1)'
    end

    #TODO: maybe we can make the function to be inplace with those Lambdas
    if settings.person_effects == false & settings.cov_effects == false
        return (Lambda_P = Lambda_P, Lambda_B_fe=Lambda_B_fe)
    elseif settings.person_effects == true & settings.cov_effects == false
        return (Lambda_P = Lambda_P, Lambda_B_fe=Lambda_B_fe, Lambda_B_pe=Lambda_B_pe)
    elseif settings.person_effects == true  & settings.cov_effects == true
        return (Lambda_P = Lambda_P, Lambda_B_fe=Lambda_B_fe, Lambda_B_pe=Lambda_B_pe, Lambda_B_cov=Lambda_B_cov)
    end

end

function eff_res(lev::JLAAlgorithm, X,id,firmid,match_id, K, settings)

    #Indexing Observations
    elist = index_constr(collect(1:length(id)), id, match_id )

    #Dimensions
    NT=size(X,1)
    M=size(elist,1)
    J = maximum(firmid)
    N = maximum(id)
    p = lev.num_simulations == 0 ? ceil(log2(NT)/0.005) : lev.num_simulations

    #Define solver
    S_xx = X'*X
    ldli, la = computeLDLinv(S_xx)
    buffs = zeros(size(la)[1],Threads.nthreads())
    compute_sol = []
    for i in 1:Threads.nthreads()
        P = approxcholOperator(ldli,buffs[:,i])
        push!(compute_sol,approxcholSolver(P,la))
    end

    #Initialize output
    Pii=zeros(M)
    Bii_fe=zeros(M)
    Bii_cov= settings.cov_effects ==true ? zeros(M) : nothing
    Bii_pe= settings.person_effects == true ? zeros(M) : nothing

    #No controls case
    if K == 0

        #Compute Auxiliaries
        movers , T = compute_movers(id, firmid)

        Nmatches = maximum(match_id)
        match_id_movers = [match_id[x] for x in findall(x->x==true, movers)]
        firmid_movers = [firmid[x] for x in findall(x->x==true, movers)]
        id_movers = [id[x] for x in findall(x->x==true, movers)]

        sel = unique(z -> match_id_movers[z], 1:length(match_id_movers))
        match_id_movers=match_id_movers[sel]
        firmid_movers=firmid_movers[sel]
        id_movers=id_movers[sel]

        maxT = maximum([T[x] for x in findall(x->x == false, movers)])

        counter = ones(Int,NT)
        gcs = Int.(@transform(groupby(DataFrame(counter = counter, id = id), :id), gcs = cumsum(:counter)).gcs)
        sel_stayers=(gcs.==1).*(movers.==false)
        stayers_matches_sel=[match_id[z] for z in findall(x->x == true , sel_stayers)]
        Tinv=1 ./T
        elist_JLL=[id_movers N.+firmid_movers id_movers N.+firmid_movers]

        M=size(elist_JLL,1)
        Pii_movers=zeros(M)
        Bii_fe_movers=zeros(M)
        Bii_cov_movers= settings.cov_effects ==true ? zeros(M) : nothing
        Bii_pe_movers= settings.person_effects == true ? zeros(M) : nothing

        #Initializing dependent variables for LSS
        Fvar= hcat(spzeros(NT,N), X[:,N+1:N+J-1])
        Dvar=hcat(X[:,1:N], spzeros(NT,J-1))

        Threads.@threads for i=1:p

            #Draw Rademacher entry
            rademach = rand(1,NT) .> 0.5
            rademach = rademach - (rademach .== 0)
            rademach = rademach ./sqrt(p)

            Z  = compute_sol[Threads.threadid()]( [rademach*X...] ; verbose=false)

            rademach = rademach .- mean(rademach)
            ZB = compute_sol[Threads.threadid()]( [rademach*Fvar...] ; verbose=false)

            if settings.person_effects == true | settings.cov_effects == true
                ZB_pe = compute_sol[Threads.threadid()]( [rademach*Dvar...] ; verbose=false)
            end

            Z = [Z;0.0]
            ZB = [ZB;0.0]
            ZB_pe = [ZB_pe;0.0]

            #Computing
            Pii_movers = Pii_movers .+ ( [Z[j]  for j in elist_JLL[:,1] ]  .- [Z[j]  for j in elist_JLL[:,2] ] ) .* ( [Z[j]  for j in elist_JLL[:,3] ]  .- [Z[j]  for j in elist_JLL[:,4] ] )
            Bii_fe_movers = Bii_fe_movers .+ ( [ZB[j]  for j in elist_JLL[:,1] ]  .- [ZB[j]  for j in elist_JLL[:,2] ] ) .* ( [ZB[j]  for j in elist_JLL[:,3] ]  .- [ZB[j]  for j in elist_JLL[:,4] ] )

            if settings.person_effects == true
                Bii_pe_movers = Bii_pe_movers .+  ( [ZB_pe[j]  for j in elist_JLL[:,1] ]  .- [ZB_pe[j]  for j in elist_JLL[:,2] ] ) .* ( [ZB_pe[j]  for j in elist_JLL[:,3] ]  .- [ZB_pe[j]  for j in elist_JLL[:,4] ] )
            end

            if settings.cov_effects == true
                ZB[N+1:end] .= ZB[N+1:end]
                ZB_pe[N+1:end] .= ZB_pe[N+1:end]
                Bii_cov_movers = Bii_cov_movers .+ ( [ZB[j]  for j in elist_JLL[:,1] ]  .- [ZB[j]  for j in elist_JLL[:,2] ] ) .* ( [ZB_pe[j]  for j in elist_JLL[:,3] ]  .- [ZB_pe[j]  for j in elist_JLL[:,4] ] )
                # Bii_cov_movers = Bii_cov_movers .+ ( [ZB[j]  for j in elist_JLL[:,1] ]  .- [ZB[j]  for j in elist_JLL[:,2] ] ) .* ( [ZB_pe[j]  for j in elist_JLL[:,3] ]  .- [ZB_pe[j]  for j in elist_JLL[:,4] ] )
            end

        end

        println("Computed Pii , Bii for movers")
        #Assign Step
        Pii_movers=sparse(match_id_movers,ones(Int,length(match_id_movers)),Pii_movers[:,1],Nmatches,1)
        Pii_stayers=sparse(stayers_matches_sel,ones(Int,length(stayers_matches_sel)),[Tinv[x] for x in findall(x->x==true,sel_stayers)],Nmatches,1)
        Pii=Pii_movers.+Pii_stayers

        Bii_fe=sparse(match_id_movers,ones(Int,length(match_id_movers)),Bii_fe_movers[:,1],Nmatches,1)

        if settings.cov_effects == true
            Bii_cov=sparse(match_id_movers,ones(Int,length(match_id_movers)),Bii_cov_movers[:,1],Nmatches,1)
        end

        if settings.person_effects == true
            Bii_pe=sparse(match_id_movers,ones(Int,length(match_id_movers)),Bii_pe_movers[:,1],Nmatches,1)
            stayers = .!movers

            Threads.@threads for t=2:maxT #T=1 have Pii=1 so need to be dropped.
                sel=(gcs.==true).*stayers.*(T.==t)
                N_sel=sum(sel)
                if N_sel > 0
                    index_sel=findall(x->x==true,sel)
                    match_sel_aux=Int.([match_id[z] for z in index_sel])
                    first=index_sel[1]
                    Xuse=X[first,:]

                    ztilde = compute_sol[Threads.threadid()]([X[first,:]...] ; verbose=false)

                    aux_right=ztilde[1:N]
                    aux_left=ztilde[1:N]

                    COV=cov(X[:,1:N]*aux_left,X[:,1:N]*aux_right)
                    Bii_pe_stayers=COV[1]*(NT-1)

                    Bii_pe_stayers=sparse(match_sel_aux,ones(Int,length(match_sel_aux)),Bii_pe_stayers,Nmatches,1)
                    Bii_pe=Bii_pe.+Bii_pe_stayers
                end
            end
        end

    elseif K> 0
        #AUX = initialize_auxiliary_variables(settings.lls_algorithm, X, elist, M,NT, N, J, K, settings)
        nparameters = N + J + K

        D=sparse(collect(1:NT),id,1)
        F=sparse(collect(1:NT),firmid,1)
        Dvar = hcat(  D, spzeros(NT,nparameters-N) )
        Fvar = hcat(spzeros(NT,N), F, spzeros(NT,nparameters-N-J) )
        #Wvar = hcat(spzeros(NT,N+J), controls )
        Xleft = X[elist[:,1],:]
        Xright = X[elist[:,2],:]

        Threads.@threads for i=1:p

            #Rademacher Entries
            rademach = rand(1,NT) .> 0.5
            rademach = rademach - (rademach .== 0)
            rademach = rademach ./sqrt(p)

            Zleft = compute_sol[Threads.threadid()]( [rademach*Xleft...] ; verbose=false)
            #Zright = lss(settings.lls_algorithm, X, rademach*Xright, settings)

            Pii = Pii .+ (X*Zleft).^2

            rademach = rademach .- mean(rademach)

            aux = compute_sol[Threads.threadid()]( [rademach*Fvar...] ;verbose=false)
            ZF = X*aux

            Bii_fe = Bii_fe .+ ZF.^2 ./NT

            if settings.person_effects == true |    settings.cov_effects == true
                aux = compute_sol[Threads.threadid()]( [rademach*Dvar...] ;verbose=false)
                ZD = X*aux
            end

            if settings.person_effects==true
                Bii_pe = Bii_pe .+ (ZD).^2 ./ NT
            end

            if settings.cov_effects==true
                Bii_cov = Bii_cov .+ (ZD .* ZF ) ./ NT
            end

        end

    end

    #Create matrices
    rows = elist[:,1]
    cols = elist[:,2]
    index_cluster = elist[:,3]

    #Censor
    Pii[ findall(Pii.>=0.99)] .= 0.99

    if K==0
        Pii = [Pii[x] for x in index_cluster]
        Bii_fe = [Bii_fe[x] for x in index_cluster]

        if settings.cov_effects == true
            Bii_cov = [Bii_cov[x] for x in index_cluster]
        end

        if settings.person_effects == true
            Bii_pe = [Bii_pe[x] for x in index_cluster]
        end

    end

    #Lambda P
    Lambda_P=sparse(rows,cols,Pii,NT,NT)
    Lambda_P=Lambda_P+triu(Lambda_P,1)'
    println("Lambda P Computed!\n")

    #Lambda B var(fe)
    Lambda_B_fe=sparse(rows,cols,Bii_fe,NT,NT)
    Lambda_B_fe=Lambda_B_fe+triu(Lambda_B_fe,1)'
    println("Lambda B FE Computed!\n")

    #Lambda B cov(fe,pe)
    if settings.cov_effects == true
        Lambda_B_cov=sparse(rows,cols,Bii_cov,NT,NT)
        Lambda_B_cov=Lambda_B_cov+triu(Lambda_B_cov,1)'
    end

    #Lambda B, var(pe)
    if settings.person_effects == true
        Lambda_B_pe=sparse(rows,cols,Bii_pe,NT,NT)
        Lambda_B_pe=Lambda_B_pe+triu(Lambda_B_pe,1)'
    end

    if settings.person_effects == false & settings.cov_effects == false
        return (Lambda_P = Lambda_P, Lambda_B_fe=Lambda_B_fe)
    elseif settings.person_effects == true & settings.cov_effects == false
        return (Lambda_P = Lambda_P, Lambda_B_fe=Lambda_B_fe, Lambda_B_pe=Lambda_B_pe)
    elseif settings.person_effects == true  & settings.cov_effects == true
        return (Lambda_P = Lambda_P, Lambda_B_fe=Lambda_B_fe, Lambda_B_pe=Lambda_B_pe, Lambda_B_cov=Lambda_B_cov)
    end

end

#7) Computes exact Pii to do inference (may be deprecated in the future)
function do_Pii(X, clustering_var)
    n=size(X,1)

    #If no clustering, Lambda_P is just diagonal matrix.
    if clustering_var == nothing
        clustering_var = collect(1:n)
    end

    #Set matrices for parallel environment.
    xx=X'*X
    #P = aspreconditioner(ruge_stuben(xx))
    ldli, la = computeLDLinv(xx)
    buffs = zeros(size(la)[1],Threads.nthreads())
    compute_sol = []
    for i in 1:Threads.nthreads()
        P = approxcholOperator(ldli,buffs[:,i])
        push!(compute_sol,approxcholSolver(P,la))
    end

    #Return the structure of the indexes associated with the clustering variable
    elist, nnz_2 = check_clustering(clustering_var)
    M=size(elist,1)

    #Set elist
    elist_1 = elist[:,1]
    elist_2 = elist[:,2]
    Pii=zeros(M)

    Threads.@threads for i=1:M
        #zexact = zeros(size(X,2))
        col = elist_2[i]
        row = elist_1[i]
        #cg!(zexact, xx, X[col,:], Pl = P , log=true, maxiter=300)
        zexact = compute_sol[Threads.threadid()]([X[col,:]...],verbose=false)

        Pii[i]= (SparseMatrixCSC{Float64,Int64}(X[col,:])'*zexact)[1]
    end

    Lambda_P = sparse(elist[:,1],elist[:,2],Pii,n,n)

    return Lambda_P
end

#8) Performs Statistical Inference on Results
function lincom_KSS(y,X,Z,Transform,clustering_var,Lambda_P; joint_test =false, labels=nothing, restrict=nothing, nsim = 10000)
    #SET DIMENSIONS
    n=size(X,1)
    K=size(X,2)
    #Add Constant
    Z=hcat(ones(size(Z,1)), Z)

    # PART 1: ESTIMATE HIGH DIMENSIONAL MODEL
    xx=X'*X
    xy=X'*y
    compute_sol = approxcholSolver(xx;verbose=true)
    beta = compute_sol([xy...];verbose=false)
    eta=y-X*beta

    # PART 1B: VERIFY LEAVE OUT COMPUTATION
    if Lambda_P == nothing
        Lambda_P=do_Pii(X,clustering_var)
    end

    if Lambda_P != nothing && clustering_var !=nothing
        nnz_1=nnz(Lambda_P)
        nnz_2=check_clustering(clustering_var).nnz_2

        if nnz_1 == nnz_2
            println("The structure of the specified Lambda_P is consistent with the level of clustering required by the user.")
        elseif nnz_1 != nnz_2
            error("The user wants cluster robust inference but the Lambda_P provided by the user is not consistent with the level of clustering asked by the user. Try to omit input Lambda_P when running lincom_KSS")
        end
    end
    I_Lambda_P = I-Lambda_P
    eta_h = I_Lambda_P\eta

    #PART 2: SET UP MATRIX FOR SANDWICH FORMULA
    rows,columns, V = findnz(Lambda_P)

    aux= 0.5*(y[rows].*eta_h[columns] + y[columns].*eta_h[rows])
    sigma_i=sparse(rows,columns,aux,n,n)

    aux= 0.5*(eta[rows].*eta[columns] + eta[columns].*eta[rows])
    sigma_i_res=sparse(rows,columns,aux,n,n)

    r=size(Z,2);
    wy=Transform*beta
    zz=Z'*Z

    numerator=Z\wy
    chet=wy-Z*numerator
    aux= 0.5*(chet[rows].*chet[columns] + chet[columns].*chet[rows])
    sigma_i_chet=sparse(rows,columns,aux,n,n)

    #PART 3: COMPUTE
    denominator=zeros(r,1)
    denominator_RES=zeros(r,1)

    for q=1:r
        v=sparse([q],[1.0],[1.0],r,1)
        v=zz\[v...]
        v=Z*v
        v=Transform'*v

        right = compute_sol(v;verbose=false)

        left=right'

        denominator[q]=left*(X'*sigma_i*X)*right
        denominator_RES[q]=left*(X'*sigma_i_res*X)*right
    end

    test_statistic=numerator./(sqrt.(denominator))
    #zz_inv=zz^(-1)
    SE_linear_combination_NAI=zz\(Z'*sigma_i_chet*Z)/zz

    #PART 4: REPORT
    println("Inference on Linear Combinations:")
    if labels == nothing
        for q=2:r
            if q <= r
                println("Linear Combination - Column Number ", q-1," of Z: ", numerator[q] )
                println("Standard Error of the Linear Combination - Column Number ", q-1," of Z: ", sqrt(denominator[q]) )
                println("T-statistic - Column Number ", q-1, " of Z: ", test_statistic[q])
            end
        end
    else
        for q=2:r
            tell_me = labels[q-1]
            println("Linear Combination associated with ", tell_me,": ", numerator[q] )
            println("Standard Error  associated with ", tell_me,": ", sqrt(denominator[q]) )
            println("T-statistic  associated with ", tell_me,": ", test_statistic[q])
        end
    end

    # PART 5: Joint-test. Quadratic form beta'*A*beta
    if joint_test == true

        if restrict  == nothing
            restrict=sparse(collect(1:r-1),collect(2:r),1.0,r-1,r)
        end

        v=restrict*(zz\(Z'*Transform))
        v=v'
        #v=sparse(v) #ldiv doesn't work for sparse RHS
        r=size(v,2)

        #Auxiliary
        aux=xx\v[:,:]
        opt_weight=v'*aux
        opt_weight=opt_weight^(-1)
        opt_weight=(1/r)*(opt_weight+opt_weight')/2

        #Eigenvalues, eigenvectors, and relevant components
        lambda , Qtilde = eigs( v*opt_weight*v', xx; nev=r,ritzvec=true)
        lambda = Real.(lambda)
        Qtilde = Real.(Qtilde)
        #lambdaS, QtildeS = eigs(v*opt_weight*v', xx; nev=1,which=:SM,ritzvec=true)
        #lambdaS = [lambdaL; lambdaS]
        #Qtilde = hcat(QtildeL,QtildeS)

        W=X*Qtilde
        V_b=W'*sigma_i*W
        V_b = (1/2)*(V_b + V_b')

        #Now focus on obtaining matrix Lambda_B with the A test associated with a joint hypothesis testing.
        Bii=opt_weight^(0.5)*aux';
        Bii=Bii*X'
        Bii=Bii'
        Bii = 0.5*(Bii[rows,:].*Bii[columns,:] + Bii[columns,:].*Bii[rows,:])
        Bii = sum(Bii,dims=2)[:]
        Lambda_B=sparse(rows,columns,Bii,n,n)

        #Leave Out Joint-Statistic
        stat=(v'*beta)'*opt_weight*(v'*beta)-y'*Lambda_B*eta_h

        #Now simulate critical values under the null.
        mu=zeros(r)
        sigma = V_b
        b_sim = MvNormal(mu,sigma)
        b_sim = rand(b_sim, nsim)

        #theta_star_sim=sum(lambda'.*(b_sim.^2 - diag(V_b)'),2)
        theta_star_sim = sum(lambda.*b_sim.^2 .- lambda.* diag(V_b),dims=1)
        pvalue=mean(theta_star_sim.>stat)

        #Report
        println("Joint-Test Statistic: ", stat)
        println("p-value: ", pvalue)

    end

    test_statistic=test_statistic[2:end]
    linear_combination=numerator[2:end]
    SE_linear_combination_KSS=sqrt.(denominator[2:end])
    SE_linear_combination_RES=sqrt.(denominator_RES[2:end])
    SE_linear_combination_NAI=diag(SE_linear_combination_NAI)
    SE_linear_combination_NAI=sqrt.(SE_linear_combination_NAI[2:end])

    return nothing
end

#9) Creates match id using firmid id
function compute_matchid(firmid,id)
    match_id2 = string.(id).*string.("+",firmid)
    match_id2 = indexin(match_id2, unique(match_id2))

    return match_id2
end

#10) Leave Out Function
function leave_out_estimation(y,id,firmid,controls,settings)

    #Create matrices for computations
    NT = size(y,1)
    J = maximum(firmid)
    N = maximum(id)
    K = controls ==nothing ? 0 : size(controls,2)
    nparameters = N + J + K

    match = compute_matchid(id, firmid)

    #Worker Dummies
    D = sparse(collect(1:NT),id,1)

    #Firm Dummies
    F = sparse(collect(1:NT),firmid,1)

    # N+J x N+J-1 restriction matrix
    S= sparse(1.0I, J-1, J-1)
    S=vcat(S,sparse(-zeros(1,J-1)))

    X = hcat(D, -F*S)

    #SET DIMENSIONS
    n=size(X,1)

    # PART 1: ESTIMATE HIGH DIMENSIONAL MODEL
    xx=X'*X
    xy=X'*y
    compute_sol = approxcholSolver(xx;verbose=true)
    beta = compute_sol([xy...];verbose=false)
    eta=y-X*beta

    pe=D * beta[1:N]
    fe=F*S * beta[N+1:N+J-1]

    println("\n","Plug-in Variance Components:")

    σ2_ψ_AKM = var(fe)
    println("Plug-in Variance of Firm Effects: ", σ2_ψ_AKM )
    σ2_α_AKM = var(pe)
    println("Plug-in Variance of Person Effects: ", σ2_α_AKM )
    σ2_ψα_AKM = cov(pe,-fe)
    println("Plug-in Covariance of Firm-Person Effects: ", σ2_ψα_AKM, "\n")

    #Part 2: Compute Pii, Bii
    LambdaP, LambdaBfe,LambdaBpe, LambdaBcov = eff_res(settings.leverage_algorithm, X,id,firmid,match, K, settings)

    #Compute Leaveo-out residual
    I_Lambda_P = I-LambdaP
    eta_h = I_Lambda_P\eta

    #Compute bias corrected variance comp of Firm Effects
    θFE = σ2_ψ_AKM -(1/NT)*y'*LambdaBfe*eta_h

    θPE = settings.person_effects==true  ? σ2_α_AKM -(1/NT)*y'*LambdaBpe*eta_h : nothing

    θCOV = settings.cov_effects==true ? σ2_ψα_AKM -(1/NT)*y'*LambdaBcov*eta_h : nothing

    return (θFE = θFE, θPE = θPE, θCOV = θCOV, β = beta, Dalpha = pe, Fpsi = fe, Pii = diag(LambdaP), Bii_pe = diag(LambdaBpe),
            Bii_fe = diag(LambdaBfe), Bii_cov = diag(LambdaBcov))

end

# Do everything naively with no inplace operations, just to get the desired result
function compute_whole(y,id,firmid,controls,settings;verbose=false)

    # compute y, id firmid, controls, settings
    verbose && println("Finding the leave-one-out connected set")
    obs,  y  , id , firmid  = find_connected_set(y,id,firmid;verbose=verbose)
    obs,  y  , id , firmid  = prunning_connected_set(y,id,firmid, obs;verbose=verbose)
    obs,  y  , id , firmid  = drop_single_obs(y,id,firmid, obs)
    controls == nothing ? nothing : controls = controls[obs,:]

    if verbose == true
        # compute the number of movers
        num_movers = length(unique(compute_movers(id,firmid).movers .* id))

        println("\n","Summary statistics of the leave-one-out connected set:")
        println("Number of observations: ", length(obs))
        println("Number of workers: ", length(unique(id)))
        println("Number of movers: ", num_movers)
        println("Mean wage: ", mean(y))
        println("Variance of wages: ", var(y))
    end

    θFE, θPE, θCOV, β, Dalpha, Fpsi, Pii, Bii_pe, Bii_fe, Bii_cov = leave_out_estimation(y,id,firmid,controls,settings)
    return (θFE = θFE, θPE = θPE, θCOV = θCOV, obs = obs, β = β, Dalpha = Dalpha, Fpsi = Fpsi, Pii = Pii, Bii_pe = Bii_pe,
            Bii_fe = Bii_fe, Bii_cov = Bii_cov)
end


text:
mutable struct AffineMean{T<:Real,V<:AbstractVector{<:Real}} <: PriorMean{T}
    w::V
    b::T
    nDim::Int
    opt::Optimizer
end

"""
**AffineMean**
```julia
    AffineMean(A::V,b::V;opt::Optimizer=Adam(α=0.01))
    AffineMean(dims::Int,features::Int;opt::Optimizer=Adam(α=0.01))
```
Construct an affine operation on `X` : `μ₀(X) = X*w + b` where `w` is a vector and `b` a scalar
Optionally give an optimizer `opt` (`Adam(α=0.01)` by default)
"""
function AffineMean(w::V,b::T,;opt::Optimizer=Adam(α=0.01)) where {V<:AbstractVector{<:Real},T<:Real}
    AffineMean{eltype(w),V}(w,b,length(w),opt)
end

function AffineMean(dims::Int;opt::Optimizer=Adam(α=0.01))
    AffineMean{Float64,Vector{Float64}}(randn(dims),0.0,dims,opt)
end

function update!(μ::AffineMean{T},grad::AbstractVector{T},X::AbstractMatrix) where {T<:Real}
    Δ = vcat(X'*grad,sum(grad))
    Δ = update(μ.opt,Δ)
    μ.w .+= Δ[1:μ.nDim]
    μ.b += Δ[end]
end

function (μ::AffineMean{T})(x::AbstractMatrix) where {T<:Real}
    @assert μ.nDim == size(x,2) "Number of dimensions do not match"
    return x*μ.w .+ μ.b
end