text:
export unreal
using KhepriBase
#=
=#

#=
julia_khepri = dirname(dirname(abspath(@__FILE__)))
=#

# Unreal is a subtype of CPP
parse_signature(::Val{:UE}, sig::T) where {T} = parse_signature(Val(:CPP), sig)

# parse_signature(::Val{:UE}, sig::AbstractString) =
#   let func_name(name) = replace(name, ":" => "_"),
#       type_name(name) = replace(name, r"[:<>]" => "_"),
#       m = match(r"^ *(public|) *(\w+) *([\[\]]*) +((?:\w|:|<|>)+) *\( *(.*) *\)", sig),
#       ret = type_name(m.captures[2]),
#       array_level = count(c -> c=='[', something(m.captures[3], "")),
#       name = m.captures[4],
#       params = split(m.captures[5], r" *, *", keepempty=false),
#       parse_c_decl(decl) =
#         let m = match(r"^ *((?:\w|:|<|>)+) *([\[\]]*) *(\w+)$", decl)
#           (type_name(m.captures[1]), count(c -> c=='[', something(m.captures[2], "")), Symbol(m.captures[3]))
#         end
#     (func_name(name), name, [parse_c_decl(decl) for decl in params], (ret, array_level))
#   end

encode(::Val{:UE}, t::Val{T}, c::IO, v) where {T} = encode(Val(:CPP), t, c, v)
decode(::Val{:UE}, t::Val{T}, c::IO) where {T} = decode(Val(:CPP), t, c)

# We need some additional Encoders
@encode_decode_as(:UE, Val{:AActor}, Val{:size})
@encode_decode_as(:UE, Val{:UStaticMesh}, Val{:size})
@encode_decode_as(:UE, Val{:UMaterial}, Val{:size})

encode(ns::Val{:UE}, t::Vector{T}, c::IO, v) where {NS,T} = begin
  sub = T()
  encode(ns, Val(:size), c, length(v))
  for e in v encode(ns, sub, c, e) end
end
decode(ns::Val{:UE}, t::Vector{T}, c::IO) where {NS,T} = begin
  sub = T()
  len = decode(ns, Val(:size), c)
  [decode(ns, sub, c) for i in 1:len]
end

# Unreal uses a left-handed CS. We need to swap X and Y.
encode(::Val{:UE}, t::Val{:FVector}, c::IO, p) =
  let (x,y,z) = raw_point(p)
    encode(Val(:CPP), Val(:float3), c, (y*100, x*100, z*100))
  end
decode(::Val{:UE}, t::Val{:FVector}, c::IO) =
  let (y,x,z) = decode(Val(:CPP), Val(:float3), c)
    xyz(x/100, y/100, z/100, world_cs)
  end

#=
encode(ns::Val{:UE}, t::Val{:Frame3d}, c::IO, v) = begin
  encode(ns, Val(:Point3d), c, v)
  t = v.cs.transform
  encode(Val(:CS), Val(:double3), c, (t[1,1], t[2,1], t[3,1]))
  encode(Val(:CS), Val(:double3), c, (t[1,2], t[2,2], t[3,2]))
  encode(Val(:CS), Val(:double3), c, (t[1,3], t[2,3], t[3,3]))
end

decode(ns::Val{:UE}, t::Val{:Frame3d}, c::IO) =
  u0(cs_from_o_vx_vy_vz(
      decode(ns, Val(:Point3d), c),
      decode(ns, Val(:Vector3d), c),
      decode(ns, Val(:Vector3d), c),
      decode(ns, Val(:Vector3d), c)))
=#
#############################################

export unreal, fast_unreal,
       unreal_material_family

# We need some additional Encoders
@encode_decode_as(:UE, Val{:AActor_array}, Vector{Val{:AActor}})
@encode_decode_as(:UE, Val{:std__vector_double_}, Vector{Float64})
@encode_decode_as(:UE, Val{:FLinearColor}, Val{RGB})
@encode_decode_as(:UE, Val{:TArray_FVector_}, Vector{Val{:FVector}})
@encode_decode_as(:UE, Val{:TArray_AActor_}, Vector{Val{:AActor}})
@encode_decode_as(:UE, Val{:TArray_TArray_FVector__}, Vector{Vector{Val{:FVector}}})

unreal_api = @remote_functions :UE """
public AActor Primitive::Sphere(FVector center, float radius)
public AActor Primitive::Box(FVector pos, FVector vx, FVector vy, FVector size)
public AActor Primitive::RightCuboid(FVector pos, FVector vx, FVector vy, float sx, float sy, float sz, float angle)
public AActor Primitive::Cylinder(FVector bottom, float radius, FVector top)
public AActor Primitive::Pyramid(TArray<FVector> ps, FVector q)
public AActor Primitive::PyramidFrustum(TArray<FVector> ps, TArray<FVector> q)
public AActor Primitive::PyramidFrustumWithMaterial(TArray<FVector> ps, TArray<FVector> q,UMaterial material)
public AActor Primitive::Slab(TArray<FVector> contour, TArray<TArray<FVector>> holes, float h, UMaterial material)
public AActor Primitive::CurrentParent()
public AActor Primitive::SetCurrentParent(AActor newParent)
public UMaterial Primitive::LoadMaterial(String name)
public UMaterial Primitive::CurrentMaterial()
public int Primitive::SetCurrentMaterial(UMaterial material)
public int Primitive::DeleteAll()
public UStaticMesh Primitive::LoadResource(String name)
public AActor Primitive::CreateBlockInstance(UStaticMesh mesh, FVector position, FVector vx, FVector vy, float scale)
public AActor Primitive::BeamRectSection(FVector position, FVector vx, FVector vy, float dx, float dy, float dz, float angle, UMaterial material)
public AActor Primitive::BeamCircSection(FVector bot, float radius, FVector top, UMaterial material)
public AActor Primitive::Panel(TArray<FVector> pts, FVector n, UMaterial material)
public AActor Primitive::Chair(FVector pos, float angle)
public AActor Primitive::InstantiateBIMElement(UStaticMesh family, FVector pos, float angle)
public AActor Primitive::Subtract(AActor ac1, AActor ac2)
public AActor Primitive::Unite(AActor ac1, AActor ac2)
public int Primitive::DeleteMany(TArray<AActor> acs)
public AActor Primitive::PointLight(FVector position, FLinearColor color, float range, float intensity)
public AActor Primitive::SetView(FVector position, FVector target, float lens, float aperture)
public FVector Primitive::ViewCamera()
public FVector Primitive::ViewTarget()
public float Primitive::ViewLens()
public int Primitive::RenderView(int width, int height, String name, String path, int frame)
public AActor Primitive::Spotlight(FVector position, FVector dir, FLinearColor color, float range, float intensity, float hotspot, float falloff);
"""

abstract type UEKey end
const UEId = Int
const UEIds = Vector{UEId}
const UERef = GenericRef{UEKey, UEId}
const UERefs = Vector{UERef}
const UEEmptyRef = EmptyRef{UEKey, UEId}
const UEUniversalRef = UniversalRef{UEKey, UEId}
const UENativeRef = NativeRef{UEKey, UEId}
const UEUnionRef = UnionRef{UEKey, UEId}
const UESubtractionRef = SubtractionRef{UEKey, UEId}
const UE = SocketBackend{UEKey, UEId}

KhepriBase.void_ref(b::UE) = UENativeRef(-1)

KhepriBase.before_connecting(b::UE) = nothing #check_plugin()
KhepriBase.after_connecting(b::UE) =
  begin
    # set_material(unity, material_basic, "Default/Materials/White")
    # set_material(unity, material_metal, "Default/Materials/Steel")
    # set_material(unity, material_glass, "Default/Materials/Glass")
    # set_material(unity, material_wood, "Default/Materials/Wood")
    # set_material(unity, material_concrete, "Default/Materials/Concrete")
    # set_material(unity, material_plaster, "Default/Materials/Plaster")
    # set_material(unity, material_grass, "Default/Materials/Grass")
  end

const unreal = UE("Unreal", unreal_port, unreal_api)

realize(b::UE, s::EmptyShape) =
  UEEmptyRef()
realize(b::UE, s::UniversalShape) =
  UEUniversalRef()

#=

realize(b::UE, s::Point) =
  UnrealPoint(connection(b), s.position)
realize(b::UE, s::Line) =
  UnrealPolyLine(connection(b), s.vertices)
realize(b::UE, s::Spline) =
  if (s.v0 == false) && (s.v1 == false)
    #UnrealSpline(connection(b), s.points)
    UnrealInterpSpline(connection(b),
                     s.points,
                     s.points[2]-s.points[1],
                     s.points[end]-s.points[end-1])
  elseif (s.v0 != false) && (s.v1 != false)
    UnrealInterpSpline(connection(b), s.points, s.v0, s.v1)
  else
    UnrealInterpSpline(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::UE, s::ClosedSpline) =
  UnrealInterpClosedSpline(connection(b), s.points)
realize(b::UE, s::Circle) =
  UnrealCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)
realize(b::UE, s::Arc) =
  if s.radius == 0
    UnrealPoint(connection(b), s.center)
  elseif s.amplitude == 0
    UnrealPoint(connection(b), s.center + vpol(s.radius, s.start_angle, s.center.cs))
  elseif abs(s.amplitude) >= 2*pi
    UnrealCircle(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
      UnrealArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, end_angle)
    else
      UnrealArc(connection(b), s.center, vz(1, s.center.cs), s.radius, end_angle, s.start_angle)
    end
  end

realize(b::UE, s::Ellipse) =
  if s.radius_x > s.radius_y
    UnrealEllipse(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
    UnrealEllipse(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::UE, s::EllipticArc) =
  error("Finish this")

realize(b::UE, s::Polygon) =
  UnrealClosedPolyLine(connection(b), s.vertices)
realize(b::UE, s::RegularPolygon) =
  UnrealClosedPolyLine(connection(b), regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed))
realize(b::UE, s::Rectangle) =
  UnrealClosedPolyLine(
    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::UE, s::SurfaceCircle) =
  UnrealSurfaceCircle(connection(b), s.center, vz(1, s.center.cs), s.radius)
realize(b::UE, s::SurfaceArc) =
    #UnrealSurfaceArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, s.start_angle + s.amplitude)
    if s.radius == 0
        UnrealPoint(connection(b), s.center)
    elseif s.amplitude == 0
        UnrealPoint(connection(b), s.center + vpol(s.radius, s.start_angle, s.center.cs))
    elseif abs(s.amplitude) >= 2*pi
        UnrealSurfaceCircle(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
            UnrealSurfaceFromCurves(connection(b),
                [UnrealArc(connection(b), s.center, vz(1, s.center.cs), s.radius, s.start_angle, end_angle),
                 UnrealPolyLine(connection(b), [add_pol(s.center, s.radius, end_angle),
                                              add_pol(s.center, s.radius, s.start_angle)])])
        else
            UnrealSurfaceFromCurves(connection(b),
                [UnrealArc(connection(b), s.center, vz(1, s.center.cs), s.radius, end_angle, s.start_angle),
                 UnrealPolyLine(connection(b), [add_pol(s.center, s.radius, s.start_angle),
                                              add_pol(s.center, s.radius, end_angle)])])
        end
    end

#realize(b::UE, s::SurfaceElliptic_Arc) = UnrealCircle(connection(b),
#realize(b::UE, s::SurfaceEllipse) = UnrealCircle(connection(b),

realize(b::UE, s::SurfacePolygon) =
  UnrealSurfacePolygon(connection(b), reverse(s.vertices))

backend_fill(b::UE, path::ClosedPolygonalPath) =
  UnrealSurfacePolygon(connection(b), path.vertices)

realize(b::UE, s::SurfaceRegularPolygon) =
  UnrealSurfaceClosedPolyLine(connection(b), regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed))
realize(b::UE, s::SurfaceRectangle) =
  UnrealSurfaceClosedPolyLine(
    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::UE, s::Surface) =
  let #ids = map(r->UENurbSurfaceFrom(connection(b),r), UnrealSurfaceFromCurves(connection(b), collect_ref(s.frontier)))
      ids = UnrealSurfaceFromCurves(connection(b), collect_ref(s.frontier))
    foreach(mark_deleted, s.frontier)
    ids
  end
backend_surface_boundary(b::UE, s::Shape2D) =
    map(shape_from_ref, UnrealCurvesFromSurface(connection(b), ref(s).value))

backend_fill(b::UE, path::ClosedPathSequence) =
  backend_fill(b, convert(ClosedPolygonalPath, path))

# Iterating over curves and surfaces

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

Unreal"public Vector3d RegionNormal(Entity ent)"
Unreal"public Point3d RegionCentroid(Entity ent)"
Unreal"public double[] SurfaceDomain(Entity ent)"
Unreal"public Frame3d SurfaceFrameAt(Entity ent, double u, double v)"

backend_surface_domain(b::UE, s::Shape2D) =
    tuple(UnrealSurfaceDomain(connection(b), ref(s).value)...)

backend_map_division(b::UE, f::Function, s::Shape2D, nu::Int, nv::Int) =
    let conn = connection(b)
        r = ref(s).value
        (u1, u2, v1, v2) = UnrealSurfaceDomain(conn, r)
        map_division(u1, u2, nu) do u
            map_division(v1, v2, nv) do v
                f(UnrealSurfaceFrameAt(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::UE, f::Function, s::SurfaceGrid, nu::Int, nv::Int) =
    let (u1, u2, v1, v2) = UnrealSurfaceDomain(conn, r)
        map_division(u1, u2, nu) do u
            map_division(v1, v2, nv) do v
                f(UnrealSurfaceFrameAt(conn, r, u, v))
            end
        end
    end

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

backend_sphere(b::UE, c, r) =
  @remote(b, Primitive__Sphere(c, r))

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

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

KhepriBase.b_pyramid_frustum(b::UE, bs, ts, bmat, tmat, smat) =
  @remote(b, Primitive__PyramidFrustum(bs, ts))

KhepriBase.b_pyramid(b::UE, bs, t, bmat, smat) =
  @remote(b, Primitive__Pyramid(bs, t))

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

realize(b::UE, s::IrregularPyramidFrustum) =
    UnrealPyramidFrustum(connection(b), s.bs, s.ts)

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

unreal"public AActor RightCuboid(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz, float angle)"
=#
#realize(b::UE, s::RightCuboid) =
#  @remote(b, Primitive__RightCuboid(s.cb, vx(1, s.cb.cs), vy(1, s.cb.cs), s.h, s.width, s.height,s.angle))

KhepriBase.b_box(b::UE, c, dx, dy, dz, mat) =
  @remote(b, Primitive__Box(c, vx(1, c.cs), vy(1, c.cs), xyz(dx, dy, dz, world_cs)))

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

realize(b::UE, s::ConeFrustum) =
  UnrealPyramidFrustum(connection(b),
    regular_polygon_vertices(64, s.cb, s.rb),
    regular_polygon_vertic@remote(b, Primitive__Sphere(s.center, s.radius))es(64, s.cb + vz(s.h, s.cb.cs), s.rt))

unreal"public AActor Cylinder(Vector3 bottom, float radius, Vector3 top)"
=#
KhepriBase.b_cylinder(b::UE, c, r, h, bmat, tmat, smat) =
  @remote(b, Primitive__Cylinder(c, r, c + vz(h, c.cs)))
#=
backend_extrusion(b::UE, s::Shape, v::Vec) =
    and_mark_deleted(
        map_ref(s) do r
            UnrealExtrude(connection(b), r, v)
        end,
        s)

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

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

backend_loft_curves(b::UE, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  and_delete_shapes(UnrealLoft(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::UE, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
    backend_loft_curves(b, profiles, rails, ruled, closed)

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

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

=#

unite_ref(b::UE, r0::UENativeRef, r1::UENativeRef) =
    ensure_ref(b, @remote(b, Primitive__Unite(r0.value, r1.value)))
#=
intersect_ref(b::UE, r0::UENativeRef, r1::UENativeRef) =
    ensure_ref(b, UnrealIntersect(connection(b), r0.value, r1.value))
=#
subtract_ref(b::UE, r0::UENativeRef, r1::UENativeRef) =
  let r = @remote(b, Primitive__Subtract(r0.value, r1.value))
    @remote(b, Primitive__DeleteMany([r0.value, r1.value]))
    r
  end

#=
slice_ref(b::UE, r::UENativeRef, p::Loc, v::Vec) =
    (UnrealSlice(connection(b), r.value, p, v); r)

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

unite_refs(b::UE, refs::Vector{<:UnrealRef}) =
    UnrealUnionRef(tuple(refs...))

#
realize(b::UE, s::UnionShape) =
  let r = foldl((r0,r1)->unite_ref(b,r0,r1), map(ref, s.shapes),
                init=UEEmptyRef())
    delete_shapes(s.shapes)
    #UnrealCanonicalize(connection(b), r.value)
    r
  end

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

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

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

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

realize(b::UE, s::Transform) =
  let r = map_ref(s.shape) do r
            UnrealTransform(connection(b), r, s.xform)
            r
          end
    mark_deleted(s.shape)
    r
  end

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

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

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

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

unreal"public AActor SurfaceFromGrid(int m, int n, Vector3[] pts, bool closedM, bool closedN, int level)"

realize(b::UE, s::SurfaceGrid) =
    UnrealSurfaceFromGrid(
        connection(b),
        size(s.points,1),
        size(s.points,2),
        reshape(s.points,:),
        s.closed_u,
        s.closed_v,
        2)

realize(b::UE, s::Thicken) =
  and_delete_shape(
    map_ref(s.shape) do r
      UnrealThicken(connection(b), r, s.thickness)
    end,
    s.shape)

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

backend_frame_at(b::UE, c::Shape1D, t::Real) = UnrealCurveFrameAt(connection(b), ref(c).value, t)

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

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

backend_frame_at(b::UE, s::Shape2D, u::Real, v::Real) = UnrealSurfaceFrameAt(connection(b), ref(s).value, u, v)

=#

# BIM

# Families

abstract type UEFamily <: Family end

struct UEMaterialFamily <: UEFamily
  name::String
  ref::Parameter{Any}
end

unreal_material_family(name, pairs...) =
  UEMaterialFamily(name, Parameter{Any}(nothing))
backend_get_family_ref(b::UE, f::Family, uf::UEMaterialFamily) =
  @remote(b, Primitive__LoadMaterial(uf.name))

struct UEResourceFamily <: UEFamily
  name::String
  parameter_map::Dict{Symbol,String}
  ref::Parameter{Any}
end

unreal_resource_family(name, pairs...) = UEResourceFamily(name, Dict(pairs...), Parameter{Any}(nothing))
backend_get_family_ref(b::UE, f::Family, uf::UEResourceFamily) = @remote(b, Primitive__LoadResource(uf.name))

backend_rectangular_table(b::UE, c, angle, family) =
    @remote(b, Primitive__InstantiateBIMElement(realize(b, family), c, -angle))

backend_chair(b::UE, c, angle, family) =
    @remote(b, Primitive__InstantiateBIMElement(realize(b, family), c, -angle))
#=
backend_rectangular_table_and_chairs(b::Khepri.Unreal, c, angle, family) =
    UnrealInstantiateBIMElement(connection(b), realize(b, family), c, -angle)

unreal"public AActor Slab(Vector3[] contour, Vector3[][] holes, float h, Material material)"
=#
backend_slab(b::UE, profile, holes, thickness, family) =
  let bot_vs = path_vertices(profile)
    @remote(b, Primitive__Slab(bot_vs, map(path_vertices, holes), thickness,realize(b, family)))
  end

#=
unreal"public AActor BeamRectSection(Vector3 position, Vector3 vx, Vector3 vy, float dx, float dy, float dz, float angle, Material material)"
unreal"public AActor BeamCircSection(Vector3 bot, float radius, Vector3 top, Material material)"
=#
realize(b::UE, 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)
    @remote(b, Primitive__BeamRectSection(
      c, vz(1, c.cs), vx(1, c.cs),
      profile.dx, profile.dy, s.h, -s.angle,
      realize(b, s.family)))
  end

realize_beam_profile(b::UE, s::Union{Beam,FreeColumn,Column}, profile::CircularPath, cb::Loc, length::Real) =
    @remote(b, Primitive__BeamCircSection(
      cb,
      profile.radius,
      add_z(cb, length*support_z_fighting_factor), #We reduce height just a bit to avoid Z-fighting
      realize(b, s.family)))

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

#=
unreal"public AActor Panel(Vector3[] pts, Vector3 n, Material material)"
=#
realize(b::UE, 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)
    @remote(b, Primitive__Panel(
      map(p -> p - n, verts),
      n*2,
      realize(b, s.family)))
  end

#
backend_wall(b::UE, w_path, w_height, l_thickness, r_thickness, family) =
  path_length(w_path) < path_tolerance() ?
    UEEmptyRef() :
    let w_paths = subpaths(w_path),
        r_w_paths = subpaths(offset(w_path, -r_thickness)),
        l_w_paths = subpaths(offset(w_path, l_thickness)),
        w_height = w_height*wall_z_fighting_factor,
        prevlength = 0,
        material = realize(b, family),
        refs = UENativeRef[]
      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)
          push!(refs, realize_pyramid_frustum(b, c_l_w_path, c_r_w_path, material))
          prevlength = currlength
        end
      end
      refs
    end
realize_pyramid_frustum(b::UE, bot_path::Path, top_path::Path, material) =
  realize_pyramid_frustum(b, path_vertices(bot_path), path_vertices(top_path), material)
realize_pyramid_frustum(b::UE, bot_vs, top_vs, material) =
    UENativeRef(@remote(b, Primitive__PyramidFrustumWithMaterial(bot_vs, top_vs, material)))

#=
sweep_fractions(b, verts, height, l_thickness, r_thickness) =
  let p = add_z(verts[1], height/2),
      q = add_z(verts[2], height/2),
      (c, h) = position_and_height(p, q),
      thickness = r_thickness + l_thickness, # HACK THIS IS WRONG!
      s = UENativeRef(@remote(b, Primitive__RightCuboid(c, vz(1, c.cs), vx(1, c.cs), height, thickness, h, 0)))
    if length(verts) > 2
      (s, sweep_fractions(b, verts[2:end], height, l_thickness, r_thickness)...)
    else
      (s, )
    end
  end

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

backend_wall_path(b::UE, path::OpenPolygonalPath, height, l_thickness, r_thickness) =
    UEUnionRef(sweep_fractions(b, path.vertices, height, l_thickness, r_thickness))

backend_wall_path(b::UE, path::Path, height, l_thickness, r_thickness) =
    backend_wall_path(b, convert(OpenPolygonalPath, path), height, l_thickness, r_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
=#
=#

backend_curtain_wall(b::UE, 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::UE, shapes::Shapes) =
  UEBoundingBox(connection(b), collect_ref(shapes))
=#
backend_set_view(b::UE, camera::Loc, target::Loc, lens::Real, aperture::Real) =
  let c = connection(b)
    @remote(b, Primitive__SetView(camera, target, lens,aperture))
  end

backend_get_view(b::UE) =
    (@remote(b, Primitive__ViewCamera()), @remote(b, Primitive__ViewTarget()), @remote(b, Primitive__ViewLens()))
#=
  zoom_extents(b::UE) = @remote(b, ZoomExtents())

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

unreal"public void DeleteAll()"
unreal"public void DeleteMany(AActor[] objs)"
=#
backend_delete_all_shapes(b::UE) = @remote(b, Primitive__DeleteAll())

backend_delete_shapes(b::UE, shapes::Shapes) =
    @remote(b, Primitive__DeleteMany(collect_ref(shapes)))
#=
set_length_unit(unit::String, b::UE) = nothing # Unused, for now

#=
# Dimensions

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

dimension(p0::Loc, p1::Loc, p::Loc, scale::Real, style::Symbol, b::UE=current_backend()) =
    UECreateAlignedDimension(connection(b), p0, p1, p,
        scale,
        UEDimensionStyles[style])

dimension(p0::Loc, p1::Loc, sep::Real, scale::Real, style::Symbol, b::UE=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
unreal"public AActor CreateParent(String name)"
unreal"public AActor CurrentParent()"
unreal"public AActor SetCurrentParent(AActor newParent)"
unreal"public void SetActive(AActor obj, bool state)"
unreal"public void DeleteAllInParent(AActor parent)"
unreal"public void SwitchToParent(AActor newParent)"
=#
UELayer = Int

current_layer(b::UE)::UELayer =
  @remote(b, Primitive__CurrentParent())

current_layer(layer::UELayer, b::UE) =
  @remote(b, Primitive__SetCurrentParent(layer))
#=
create_layer(name::String, b::UE) =
  UECreateParent(connection(b), name)

set_layer_active(layer::UELayer, status, b::UE) =
  let c = connection(b)
    UESetActive(c, layer, status)
    interrupt_processing(c)
  end

delete_all_shapes_in_layer(layer::UELayer, b::UE) =
  UEDeleteAllInParent(connection(b), layer)

switch_to_layer(layer::UELayer, b::UE) =
  UESwitchToParent(connection(b), layer)

# Experiment to speed up things

canonicalize_layer(layer::UELayer, b::UE) =
  UECanonicalize(connection(b), layer)

# Materials
=#
UEMaterial = Int
#=
unreal"public Material LoadMaterial(String name)"
unreal"public void SetCurrentMaterial(Material material)"
unreal"public Material CurrentMaterial()"
=#
current_material(b::UE)::UEMaterial =
  @remote(b, Primitive__CurrentMaterial())

current_material(material::UEMaterial, b::UE) =
  @remote(b, Primitive__SetCurrentMaterial(material))

get_material(name::String, b::UE) =
  @remote(b, Primitive__LoadMaterial(name))

#=
# Blocks

unreal"public AActor CreateBlockInstance(AActor block, Vector3 position, Vector3 vx, Vector3 vy, float scale)"
unreal"public AActor CreateBlockFromShapes(String name, AActor[] objs)"
=#
realize(b::UE, s::Block) =
    @remote(b, Primitive__LoadResource(s.name))

realize(b::UE, s::BlockInstance) =
    @remote(b, Primitive__CreateBlockInstance(
        ref(s.block).value,
        s.loc, vx(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
unreal"public AActor PointLight(Vector3 position, Color color, float range, float intensity)"
=#
backend_pointlight(b::UE, loc::Loc, color::RGB, range::Real, intensity::Real) =
    @remote(b, Primitive__PointLight(loc, color, range, intensity))

  backend_spotlightricardo(b::UE, loc::Loc,dir::Vec, color::RGB, range::Real, intensity::Real,hotspot::Real,falloff::Real) =
  @remote(b, Primitive__Spotlight( loc,dir,color,range,intensity, hotspot, falloff))
#=
backend_spotlight(b::UE, loc::Loc, dir::Vec, hotspot::Real, falloff::Real) =
      @remote(b, SpotLight( loc, hotspot, falloff, loc + dir))

backend_ieslight(b::UE, file::String, loc::Loc, dir::Vec, alpha::Real, beta::Real, gamma::Real) =
    UEIESLight(connection(b), file, loc, loc + dir, vxyz(alpha, beta, gamma))

#=
# User Selection
=#

shape_from_ref(r, b::UE) =
  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, UENativeRef(r), 0, 0))
          #code = UEShapeCode(c, r),
          #ref = LazyRef(b, UENativeRef(r))
          #error("Unknown shape with code $(code)")
      end
    else
      collected_shapes()[idx]
    end
  end
#
#=
UE"public Point3d[] GetPosition(string prompt)"

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

UE"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::UE) =
  select_with_prompt(prompt, b, UEGetPoint)

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

select_curve(prompt::String, b::UE) =
  select_with_prompt(prompt, b, UEGetCurve)

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

select_surface(prompt::String, b::UE) =
  select_with_prompt(prompt, b, UEGetSurface)

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

select_solid(prompt::String, b::UE) =
  select_with_prompt(prompt, b, UEGetSolid)

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

select_shape(prompt::String, b::UE) =
  select_with_prompt(prompt, b, UEGetShape)

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

captured_shape(b::UE, handle) =
  shape_from_ref(UEGetShapeFromHandle(connection(b), handle),
                 b)

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

# Register for notification

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

register_for_changes(s::Shape, b::UE) =
    let conn = connection(b)
        UERegisterForChanges(conn, ref(s).value)
        UEDetectCancel(conn)
        s
    end

unregister_for_changes(s::Shape, b::UE) =
    let conn = connection(b)
        UEUnregisterForChanges(conn, ref(s).value)
        UEUndetectCancel(conn)
        s
    end

waiting_for_changes(s::Shape, b::UE) =
    ! UEWasCanceled(connection(b))

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

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

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

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

disable_update(b::UE) =
    UEDisableUpdate(connection(b))

enable_update(b::UE) =
    UEEnableUpdate(connection(b))
# Render

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

#render exposure: [-3, +3] -> [-6, 21]
convert_render_exposure(b::UE, v::Real) = -4.05*v + 8.8
#render quality: [-1, +1] -> [+1, +50]
convert_render_quality(b::UE, v::Real) = round(Int, 25.5 + 24.5*v)
=#
render_view(path::String, b::UE) =
    begin
      @remote(b, Primitive__RenderView(render_width(), render_height(),film_filename(),path,film_frame()))
      path
    end
#=
unreal"public void SelectAActors(AActor[] objs)"

highlight_shape(s::Shape, b::UE) =
    UESelectAActors(connection(b), collect_ref(s))

highlight_shapes(ss::Shapes, b::UE) =
    UESelectAActors(connection(b), collect_ref(ss))

unreal"public void StartSelectingAActor()"
unreal"public int SelectedAActorId(bool existing)"

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


text:
include("../physical_params.jl")
using Plots

##################
#  MMS Function  #
##################

function ϕ(x, y, MMS)
    h = MMS.H
    return (h .* (h .+ x)) ./
        ((h .+ x).^2 .+ y.^2)
end

function ϕ_x(x, y, MMS)
    h = MMS.H
    return (h .* (y.^2 .- (h .+ x).^2)) ./ (((h .+ x).^2 .+ y.^2).^2)
end

function ϕ_xx(x, y, MMS)
    h = MMS.H
    return (2*h .* (h .+ x) .* (h^2 .+ 2*h .* x .+ x.^2 .- 3*y.^2)) ./
        ((h^2 .+ 2*h .* x .+ x.^2 .+ y.^2).^3)
end

function ϕ_y(x, y, MMS)
    h = MMS.H
    return - (2 * h .* y .* (h .+ x)) ./
        ((h .+ x).^2 .+ y.^2).^2
end

function ϕ_yy(x, y, MMS)
    h = MMS.H
    return (2*h .* (h .+ x) .* (3 .* y.^2 .- (h .+ x).^2)) ./
        (((h .+ x).^2 .+ y.^2).^3)
end

function K(t, MMS)

    t̄ = MMS.t̄
    tw = MMS.t_w
    δ = MMS.δ_e
    Vm = MMS.Vmin

    return 1/π * (atan((t - t̄)/tw) + π/2) + Vm/δ * t
end

function K_t(t, MMS)

    t̄ = MMS.t̄
    tw = MMS.t_w
    δ = MMS.δ_e
    Vm = MMS.Vmin

    return tw / (π * ((t - t̄)^2 + tw^2)) + Vm/δ
end

function K_tt(t, MMS)

    t̄ = MMS.t̄
    tw = MMS.t_w
    δ = MMS.δ_e
    Vm = MMS.Vmin

    return -(2 * tw * (t - t̄)) / (π * ((t - t̄)^2 + tw^2)^2)
end

he(x, y, t, MMS) = MMS.δ_e/2 .* K(t, MMS) .* ϕ(x, y, MMS) .+ MMS.Vp/2 .* t .* (1 .- ϕ(x,y,MMS)) .+
    MMS.τ∞/24 .* x

he_t(x, y, t, MMS) = MMS.δ_e/2 .* K_t(t, MMS) .* ϕ(x, y, MMS) .+ MMS.Vp/2 .* (1 .- ϕ(x,y,MMS))

he_tt(x, y, t, MMS) = MMS.δ_e/2 .* K_tt(t, MMS) .* ϕ(x, y, MMS)

he_x(x, y, t, MMS) =  MMS.δ_e/2 .* K(t, MMS) .* ϕ_x(x, y, MMS) .- MMS.Vp/2 .* t .* ϕ_x(x, y, MMS) .+ MMS.τ∞/24

he_xx(x, y, t, MMS) =  MMS.δ_e/2 .* K(t, MMS) .* ϕ_xx(x, y, MMS) .- MMS.Vp/2 .* t .* ϕ_xx(x, y, MMS)

he_y(x, y, t, MMS) =  MMS.δ_e/2 .* K(t, MMS) .* ϕ_y(x, y, MMS) .- MMS.Vp/2 .* t .* ϕ_y(x, y, MMS)

he_yy(x, y, t, MMS) =  MMS.δ_e/2 .* K(t, MMS) .* ϕ_yy(x, y, MMS) .- MMS.Vp/2 .* t .* ϕ_yy(x, y, MMS)

he_xt(x, y, t, MMS) = MMS.δ_e/2 .* K_t(t, MMS) .* ϕ_x(x, y, MMS) .- MMS.Vp/2 .* ϕ_x(x, y, MMS)

function τhe(fx, fy, t, fnum, B_p, MMS)

    if fnum == 1
        τ = -μ(fx, fy, B_p) .* he_x(fx, fy, t, MMS)
    elseif fnum == 2
        τ = μ(fx, fy, B_p) .* he_x(fx, fy, t, MMS)
    elseif fnum == 3
        τ = -μ(fx, fy, B_p) .* he_y(fx, fy, t, MMS)
    elseif fnum == 4
        τ = μ(fx, fy, B_p) .* he_y(fx, fy, t, MMS)
    end
    return τ

end

### Dynamic h solution forcing functions

function Forcing_hd(x, y, t, B_p, MMS)

    Force = ρ(x, y, B_p) .* he_tt(x, y, t, MMS) .-
        (μ_x(x, y, B_p) .* he_x(x, y, t, MMS) .+ μ(x, y, B_p) .* he_xx(x, y, t, MMS) .+
         μ_y(x, y, B_p) .* he_y(x, y, t, MMS) .+ μ(x, y, B_p) .* he_yy(x, y, t, MMS))

    return Force
end

function S_ch(fx, fy, t, fnum, R, B_p, MMS)

    Z = sqrt.(μ(fx, fy, B_p) .* ρ(fx, fy, B_p))
    v = he_t(fx, fy, t, MMS)
    τ = τhe(fx, fy, t, fnum, B_p, MMS)

    return Z .* v .+ τ .- R .* (Z .* v .- τ)
end

function ψe_hd(x, y, t, B_p, RS, MMS)

    τe = τhe(x, y, t, 1, B_p, MMS)
    Ve = 2 .* he_t(x, y, t, MMS)

    return RS.a .* log.((2 * RS.V0 ./ Ve) .* sinh.(-τe ./ (RS.a .* RS.σn)))
end

function ψe_thd(x, y, t, B_p, RS, MMS)

    τe = τhe(x, y, t, 1, B_p, MMS)
    Ve = 2 * he_t(x, y, t, MMS)
    Ve_t = 2 * he_tt(x, y, t, MMS)
    τe_t = -μ(x, y, B_p) .* he_xt(x, y, t, MMS)

    ψ_t = τe_t .* coth.(τe ./ (RS.a * RS.σn)) ./ RS.σn .- RS.a .* Ve_t ./ Ve

    return ψ_t

end

function S_rsdh(fx, fy, b, t, B_p, RS, MMS)

    ψ = ψe_hd(fx, fy, t, B_p, RS, MMS)
    V = 2*he_t(fx, fy, t, MMS)
    G = (b .* RS.V0 ./ RS.Dc) .* (exp.((RS.f0 .- ψ) ./ b) .- abs.(V) / RS.V0)
    ψ_t = ψe_thd(fx, fy, t, B_p, RS, MMS)
    return  ψ_t .- G

end

### Quasi-dynamic h solution forcing functions

Forcing_h(x, y, t, B_p, MMS) = -(μ_x(x, y, B_p) .* he_x(x, y, t, MMS) .+
    μ(x, y, B_p) .* he_xx(x, y, t, MMS) .+
    μ_y(x, y, B_p) .* he_y(x, y, t, MMS) .+
    μ(x, y, B_p) .* he_yy(x, y, t, MMS))

function ψe_h(x, y, t, B_p, RS, MMS)

    τ = τhe(x, y, t, 1, B_p, MMS)
    Ve = 2 .* he_t(x, y, t, MMS)

    return RS.a .* log.((2 * RS.V0 ./ Ve) .* sinh.((-τ #=.- η(y, B_p) .* Ve=#) ./ (RS.a .* RS.σn)))

end

function ψe_th(x, y, t, B_p, RS, MMS)

    τ = τhe(x, y, t, 1, B_p, MMS)
    Ve = 2 * he_t(x, y, t, MMS)
    Ve_t = 2 * he_tt(x, y, t, MMS)
    τe_t = - μ(x, y, B_p) .* he_xt(x, y, t, MMS)

    ψ_t = #=(=#-τe_t #=.- η(y, B_p) .* Ve_t)=# .* coth.((-τ #=.- η(y, B_p) .* Ve=#) ./ (RS.a * RS.σn)) ./ RS.σn .- RS.a .* Ve_t ./ Ve

    return ψ_t[1]

end

function S_rsh(x, y, t, b, B_p, RS, MMS)

    ψ = ψe_h(x, y, t, B_p, RS, MMS)
    Ve = 2 * he_t(x, y, t, MMS)
    G = (b .* RS.V0 ./ RS.Dc) .* (exp.((RS.f0 .- ψ) ./ b) .- abs.(Ve) / RS.V0)

    s_rs = ψe_th(x, y, t, B_p, RS, MMS) .- G

    return s_rs[1]

end

f(a, MMS) = MMS.amp * sin.(π.*(a)/MMS.wl)
fp(a, MMS) = MMS.amp * π/MMS.wl*cos.(π.*(a)/MMS.wl)
fpp(a, MMS) = MMS.amp * (-(π/MMS.wl)^2) .* sin.(π.*(a)/MMS.wl)

ue(x, y, t, MMS) = f(x .+ y .- t, MMS) .+ (MMS.amp * π/MMS.wl + MMS.ϵ)*t .+ (MMS.amp * π/MMS.wl
                                                                             + MMS.ϵ)*x
ue_t(x, y, t, MMS) = -fp(x .+ y .- t, MMS) .+ (MMS.amp * π/MMS.wl + MMS.ϵ)
ue_tt(x, y, t, MMS) = fpp(x .+ y .- t, MMS)

ue_x(x, y, t, MMS) = fp(x .+ y .- t, MMS) .+ (MMS.amp * π/MMS.wl + MMS.ϵ)
ue_y(x, y, t, MMS) = fp(x .+ y .- t, MMS)
ue_xy(x,y, t, MMS) = fpp(x .+ y .- t, MMS)
ue_xx(x, y, t, MMS) = fpp(x .+ y .- t, MMS)
ue_yy(x, y, t, MMS) = fpp(x .+ y .- t, MMS)
ue_xt(x, y, t, MMS) = -fpp(x .+ y .- t, MMS)
ue_yt(x, y, t, MMS) = -fpp(x .+ y .- t, MMS)

function τe(fx, fy, t, fnum, B_p, MMS)

    if fnum == 1
        τ = -μ(fx, fy, B_p) .* ue_x(fx, fy, t, MMS)
    elseif fnum == 2
        τ = μ(fx, fy, B_p) .* ue_x(fx, fy, t, MMS)
    elseif fnum == 3
        τ = -μ(fx, fy, B_p) .* ue_y(fx, fy, t, MMS)
    elseif fnum == 4
        τ = μ(fx, fy, B_p) .* ue_y(fx, fy, t, MMS)
    end
    return τ

end

function τe_t(fx, fy, t, fnum, B_p, MMS)

    if fnum == 1
        τ = -μ(fx, fy, B_p) .* ue_xt(fx, fy, t, MMS)
    elseif fnum == 2
        τ = μ(fx, fy, B_p) .* ue_xt(fx, fy, t, MMS)
    elseif fnum == 3
        τ = -μ(fx, fy, B_p) .* ue_yt(fx, fy, t, MMS)
    elseif fnum == 4
        τ = μ(fx, fy, B_p) .* ue_yt(fx, fy, t, MMS)
    end

    return τ

end

### Dynamic u solution forcing functions

function ψe_d(fx, fy, t, B_p, RS, MMS)

    Ve = 2*ue_t(fx, fy, t, MMS)
    τf = τe(fx, fy, t, 1, B_p, MMS)

    return RS.a .* log.((2 * RS.V0 ./ Ve) .* sinh.(-τf ./ (RS.a .* RS.σn)))

end

function ψe_dt(fx, fy, t, B_p, RS, MMS)

    Ve = 2*ue_t(fx, fy, t, MMS)
    Ve_t = 2*ue_tt(fx, fy, t, MMS)
    τf = τe(fx, fy, t, 1, B_p, MMS)
    τf_t = τe_t(fx, fy, t, 1, B_p, MMS)

    return τf_t ./ RS.σn .* coth.(τf ./ (RS.σn .* RS.a)) - RS.a .* Ve_t ./ Ve

end

function Forcing_ud(x, y, t, B_p, MMS)

    Force = ρ(x, y, B_p) .* ue_tt(x, y, t, MMS) .-
        (μ_x(x, y, B_p) .* ue_x(x, y, t, MMS) .+ μ(x, y, B_p) .* ue_xx(x, y, t, MMS) .+
         μ_y(x, y, B_p) .* ue_y(x, y, t, MMS) .+ μ(x, y, B_p) .* ue_yy(x, y, t, MMS))

    return Force
end

function S_cd(fx, fy, t, fnum, R, B_p, MMS)

    Z = sqrt.(μ(fx, fy, B_p) .* ρ(fx, fy, B_p))
    v = ue_t(fx, fy, t, MMS)
    τ = τe(fx, fy, t, fnum, B_p, MMS)

    return Z .* v .+ τ .- R .* (Z .* v .- τ)
end

function S_rsd(fx, fy, b, t, B_p, RS, MMS)

    ψ = ψe(fx, fy, t, B_p, RS, MMS)
    V = 2*ue_t(fx, fy, t, MMS)
    G = (b .* RS.V0 ./ RS.Dc) .* (exp.((RS.f0 .- ψ) ./ b) .- abs.(V) / RS.V0)
    ψ_t = ψe_t(fx, fy, t, B_p, RS, MMS)
    return  ψ_t .- G

end

### Quasi-dynamic u solution forcing functions

Forcing_u(x, y, t, B_p, MMS) = -(μ_x(x, y, B_p) .* ue_x(x, y, t, MMS) .+
    μ(x, y, B_p) .* ue_xx(x, y, t, MMS) .+
    μ_y(x, y, B_p) .* ue_y(x, y, t, MMS) .+
    μ(x, y, B_p) .* ue_yy(x, y, t, MMS))

function ψe(x, y, t, B_p, RS, MMS)

    τ = τe(x, y, t, 1, B_p, MMS)
    Ve = 2 .* ue_t(x, y, t, MMS)

    return RS.a .* log.((2 * RS.V0 ./ Ve) .* sinh.((-τ .- η(y, B_p) .* Ve) ./ (RS.a .* RS.σn)))

end

function ψe_t(x, y, t, B_p, RS, MMS)

    τ = τe(x, y, t, 1, B_p, MMS)
    Ve = 2 * ue_t(x, y, t, MMS)
    Ve_t = 2 * ue_tt(x, y, t, MMS)
    τe_t = - μ(x, y, B_p) .* ue_xt(x, y, t, MMS)

    ψ_t = (-τe_t .- η(y, B_p) .* Ve_t) .* coth.((-τ .- η(y, B_p) .* Ve) ./ (RS.a * RS.σn)) ./ RS.σn .- RS.a .* Ve_t ./ Ve

    return ψ_t

end

function S_rs(x, y, t, b, B_p, RS, MMS)

    ψ = ψe(x, y, t, B_p, RS, MMS)
    Ve = 2 * ue_t(x, y, t, MMS)
    G = (b .* RS.V0 ./ RS.Dc) .* (exp.((RS.f0 .- ψ) ./ b) .- abs.(Ve) / RS.V0)

    s_rs = ψe_t(x, y, t, B_p, RS, MMS) .- G

    return s_rs

end

function u_face2(x, y, t, MMS, RS, μf2)
    return ue(x, y, t, MMS) .- (MMS.Vp/2 * t .+ (RS.τ_inf * MMS.Lw) ./ μf2)

end

ext: jl


text:
# ---
# title: 717. 1-bit and 2-bit Characters
# id: problem717
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/1-bit-and-2-bit-characters/description/>
# hidden: true
# ---
#
# We have two special characters. The first character can be represented by one
# bit `0`. The second character can be represented by two bits (`10` or `11`).
#
# Now given a string represented by several bits. Return whether the last
# character must be a one-bit character or not. The given string will always end
# with a zero.
#
# **Example 1:**
#
#
#
#     Input:
#     bits = [1, 0, 0]
#     Output: True
#     Explanation:
#     The only way to decode it is two-bit character and one-bit character. So the last character is one-bit character.
#
#
# **Example 2:**
#
#
#
#     Input:
#     bits = [1, 1, 1, 0]
#     Output: False
#     Explanation:
#     The only way to decode it is two-bit character and two-bit character. So the last character is NOT one-bit character.
#
#
# **Note:**
#
# * `1 <= len(bits) <= 1000`.
# * `bits[i]` is always `0` or `1`.
#
#
## @lc code=start
using LeetCode

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


text:
##########################################################
#
#  Proximal Operators
#
##########################################################

abstract type ProximableFunction end

prox!(g::ProximableFunction, hat_x::AbstractArray, x::AbstractArray) = prox!(g, hat_x, x, one(eltype(x)))
prox(g::ProximableFunction,                        x::AbstractArray, γ) = prox!(g, similar(x), x, γ)
prox(g::ProximableFunction,                        x::AbstractArray) = prox!(g, similar(x), x, one(eltype(x)))

##########################################################
######  zero
##########################################################

struct ProxZero <: ProximableFunction end

value(g::ProxZero, x::AbstractVecOrMat{T}) where {T} = zero(T)

function prox!(::ProxZero, out_x::AbstractVecOrMat{T}, x::AbstractVecOrMat{T}, γ::T) where {T}
  size(out_x) == size(x) || throw(ArgumentError("Sizes of the input and ouput need to be the same."))
  copyto!(out_x, x)
end

cdprox!(g::ProxZero, x::Union{SparseIterate{T},SymmetricSparseIterate{T}}, k::Int, γ::T) where {T} = x[k]

##########################################################
######  L1 norm  g(x) = λ0 ⋅ \sum_j λ_j |x_j|
##########################################################

struct ProxL1{T<:AbstractFloat, S} <: ProximableFunction
  λ0::T
  λ::S

  ProxL1{T, S}(λ0::T, λ::Union{Nothing, AbstractArray{T}}) where {T <: AbstractFloat, S} = new(λ0, λ)
end

ProxL1(λ0::T        ) where {T <: AbstractFloat} = ProxL1{T, Nothing}(λ0, nothing)
ProxL1(λ0::T, ::Nothing) where {T <: AbstractFloat} = ProxL1{T, Nothing}(λ0, nothing)
ProxL1(λ0::T, λ::AbstractArray{T}) where {T <: AbstractFloat} = ProxL1{T, typeof(λ)}(λ0, λ)

function value(g::ProxL1{T, S}, x::AbstractArray{T}) where S <: AbstractArray{T} where {T <: AbstractFloat}
  size(g.λ) == size(x) || throw(DimensionMismatch("Sizes of g.λ and x need to be the same"))
  v = zero(T)
  @inbounds @simd for i in eachindex(x)
    v += abs(x[i]) * g.λ[i]
  end
  v * g.λ0
end
prox!(g::ProxL1{T, S}, out_x::AbstractArray{T}, x::AbstractArray{T}, γ::T) where S <: AbstractArray{T} where {T <: AbstractFloat} =
    out_x .= shrink.(x, (g.λ0 * γ) * g.λ)
@inline function cdprox!(g::ProxL1{T, S}, x::Union{SparseIterate{T},SymmetricSparseIterate{T}}, k::Int, γ::T) where S <: AbstractArray{T} where {T <: AbstractFloat}
  length(g.λ) == numCoordinates(x) || throw(DimensionMismatch())
  @boundscheck checkbounds(x, k)
  x[k] = shrink(x[k], g.λ[k] * γ * g.λ0)
end

value(g::ProxL1{T, Nothing}, x::AbstractArray{T}) where {T<:AbstractFloat} = g.λ0 * sum(abs, x)
prox!(g::ProxL1{T, Nothing}, out_x::AbstractArray{T}, x::AbstractArray{T}, γ::T) where {T<:AbstractFloat} =
  out_x .= shrink.(x, γ * g.λ0)
cdprox!(g::ProxL1{T, Nothing}, x::Union{SparseIterate{T},SymmetricSparseIterate{T}}, k::Int, γ::T) where {T} =
  x[k] = shrink(x[k], g.λ0 * γ)

##########################################################
###### L1 + Fused  g(x1,x2) = λ1*(|x1|+|x2|) + λ2*|x1-x2|
##########################################################

struct ProxL1Fused{T<:AbstractFloat} <: ProximableFunction
  λ1::T
  λ2::T
end

value(g::ProxL1Fused{T}, x1::T, x2::T) where {T<:AbstractFloat} = g.λ1*(abs(x1)+abs(x2)) + g.λ2*abs(x1-x2)
function prox!(g::ProxL1Fused{T}, out_x::Tuple{StridedVector{T},StridedVector{T}}, x::Tuple{StridedVector{T},StridedVector{T}}, γ::T) where {T<:AbstractFloat}
  @assert size(out_x[1]) == size(x[1])
  @assert size(out_x[2]) == size(x[2])
  λ1 = g.λ1
  λ2 = g.λ2
  @inbounds @simd for i in eachindex(x[1])
    out_x[1][i], out_x[2][i] = proxL1Fused(x[1][i], x[2][i], λ1, λ2, γ)
  end
  out_x
end

##########################################################
###### L2 norm   g(x) = λ * ||x||_2
##########################################################

struct ProxL2{T<:AbstractFloat, S} <: ProximableFunction
  λ0::T
  λ::S

  ProxL2{T, S}(λ0::T, λ::Union{Nothing, AbstractVector{T}}) where {T <: AbstractFloat, S} = new(λ0, λ)
end

ProxL2(λ0::T        ) where {T <: AbstractFloat} = ProxL2{T, Nothing}(λ0, nothing)
ProxL2(λ0::T, ::Nothing) where {T <: AbstractFloat} = ProxL2{T, Nothing}(λ0, nothing)
ProxL2(λ0::T, λ::AbstractVector{T}) where {T <: AbstractFloat} = ProxL2{T, typeof(λ)}(λ0, λ)

value(g::ProxL2{T, Nothing}, x::AbstractArray{T}) where {T<:AbstractFloat} = g.λ0 * norm(x)
prox!(g::ProxL2{T, Nothing}, out_x::AbstractArray{T}, x::AbstractArray{T}, γ::T) where {T<:AbstractFloat} =
  shrinkL2!(out_x, x, g.λ0*γ)

# λ0⋅sum_k λ_k g(x_k)
function value(g::ProxL2{T, S}, x::AtomIterate{T}) where S <: AbstractVector{T} where {T<:AbstractFloat}
  length(x.atoms) == length(g.λ) || throw(DimensionMismatch())

  v = zero(T)
  @inbounds for i=1:length(x.atoms)
    v += norm(x.atoms[i]) * g.λ[i]
  end
  v * g.λ0
end

function prox!(g::ProxL2{T, S}, out::AtomIterate{T}, x::AtomIterate{T}, γ::T) where S <: AbstractVector{T} where {T<:AbstractFloat}
  length(x.atoms) == length(g.λ) || throw(DimensionMismatch())

  for i=1:length(x.atoms)
    shrinkL2!(out.atoms[i], x.atoms[i], γ*g.λ[i]*g.λ0)
  end
  out
end
cdprox!(g::ProxL2{T, S}, x::AtomIterate{T}, k::Int, γ::T) where S <: AbstractVector{T} where {T<:AbstractFloat} =
  shrinkL2!(x.atoms[k], x.atoms[k], g.λ[k] * γ * g.λ0)

##########################################################
###### L2 norm squared g(x) = λ * ||x||_2^2
##########################################################

struct ProxL2Sq{T<:AbstractFloat} <: ProximableFunction
  λ::T
end
value(g::ProxL2Sq, x::StridedArray) = g.λ * sum(abs2, x)
function prox!(g::ProxL2Sq{T}, out_x::AbstractVecOrMat{T}, x::AbstractVecOrMat{T}, γ::T) where {T<:AbstractFloat}
  size(out_x) == size(x) || throw(ArgumentError("Sizes of the input and ouput need to be the same."))
  c = g.λ * γ
  c = 1. / (1. + 2. * c)
  copyto!(out_x, x)
  rmul!(out_x, c)
end

##########################################################
####### nuclear norm  g(x) = λ * ||x||_*
###########################################################

struct ProxNuclear{T<:AbstractFloat} <: ProximableFunction
  λ::T
end
value(g::ProxNuclear, x::StridedMatrix) = g.λ * sum(svdvals(x))
function prox!(g::ProxNuclear{T}, out_x::StridedMatrix{T}, x::StridedMatrix{T}, γ::T) where {T<:AbstractFloat}
  size(out_x) == size(x) || throw(ArgumentError("Sizes of the input and ouput need to be the same."))
  U, S, V = svd(x)
  c = g.λ * γ
  S .= shrink.(S, c)
  rmul!(U, Diagonal(S))
  mul!(out_x, U, transpose(V))
end

##################################
#
#  sum_k g(x_k)
#
##################################

##########################################################
# Gaussian likelihood prox
# f(X) = tr(SX) - log deg(X)
##########################################################

struct ProxGaussLikelihood{T} <: ProximableFunction
  S::Symmetric{T}
  tmp::Matrix{T}
end
ProxGaussLikelihood(S::Symmetric) = ProxGaussLikelihood{eltype(S)}(S, zeros(eltype(S),size(S)))

value(g::ProxGaussLikelihood{T}, X::StridedMatrix{T}) where {T<:AbstractFloat} = tr(g.S*X) - logdet(X)
function prox!(g::ProxGaussLikelihood{T}, hX::StridedMatrix{T}, X::StridedMatrix{T}, γ::T) where {T<:AbstractFloat}
  S = g.S
  tmp = g.tmp

  ρ = one(T) / γ
  @. tmp = X * ρ - S
  ef = eigen!(Symmetric(tmp))
  d =  ef.values::Vector{T}
  U = ef.vectors::Matrix{T}
  @inbounds @simd for i in eachindex(d)
    t = d[i]
    d[i] = sqrt( (t + sqrt(t^2. + 4. * ρ)) / (2. * ρ) )
  end
  rmul!(U, Diagonal(d))
  mul!(hX, U, transpose(U))
end

###


text:
using NumDisplay

println("Starting...")

d = DisplayBCD([5, 6, 13, 19], (2, 3, 4, 17, 27))
sleep(3)

shutdown_mode_off(d)
sleep(3)

###

println("Low level method write_digit()...")

println("write 3210")
write_digit(d, UInt8(0), 1)
write_digit(d, UInt8(1), 2)
write_digit(d, UInt8(2), 3)
write_digit(d, UInt8(3), 4)

sleep(3)
println("write 7654")
write_digit(d, UInt8(4), 1)
write_digit(d, UInt8(5), 2)
write_digit(d, UInt8(6), 3)
write_digit(d, UInt8(7), 4)
sleep(3)

println("write ??98")
write_digit(d, UInt8(8), 1)
write_digit(d, UInt8(9), 2)
write_digit(d, UInt8(10), 3)
write_digit(d, UInt8(11), 4)
sleep(3)

println("write  ???")
write_digit(d, UInt8(12), 1)
write_digit(d, UInt8(13), 2)
write_digit(d, UInt8(14), 3)
write_digit(d, UInt8(15), 4)
sleep(3)

###

println("High level method write_number()...")

data = [nothing,6,6,6]
write_number(d, data)
println("$data\n")
sleep(3)

println("Maximal intensity 16")
set_intensity(d, 16)
sleep(3)

println("change intensity from 16 to 1")
for i in collect(16:(-1):1)
    set_intensity(d, i)
    println("$i/16")
    sleep(0.5)
end
sleep(3)

println("shutdown_mode_on()")
shutdown_mode_on(d)
sleep(3)

println("shutdown_mode_off()")
shutdown_mode_off(d)
sleep(3)

println("test_mode_on()")
test_mode_on(d)
sleep(3)

println("test_mode_off()")
test_mode_off(d)
sleep(3)

println("Stop.")
shutdown_mode_off(d)


text:
using MultipleViewGeometry, Test, Random
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using MultipleViewGeometry.ModuleConstraints
using MultipleViewGeometry.ModuleConstruct
using LinearAlgebra
using StaticArrays
using GeometryTypes

# Fix random seed.
Random.seed!(1234)
# Construct two camera matrices and parametrise two planar surfaces.
f = 50
image_width = 640 / 10
image_height = 480 / 10
𝐊₁ = @SMatrix [f 0 0 ;
               0 f 0 ;
               0 0 1 ]
𝐑₁ = SMatrix{3,3,Float64,9}(rotxyz(0, 25*(pi/180), 0))
𝐭₁ = [-30.0, 0.0, -5.0]

𝐊₂ = @SMatrix [f 0 0 ;
               0 f 0 ;
               0 0 1 ]

𝐑₂ = SMatrix{3,3,Float64,9}(rotxyz(0, -25*(pi/180), 0))
𝐭₂ = [30.0, 0.0, 5.0]

# Normals and distance from origin
𝐧₁ = [0.0, 0.0, 1.0]
d₁ = 55.0

# Normals and distance from origin
𝐧₁ = [0.0, 0.0, 1.0]
d₁ = 55.0

𝐧₂ = [0.5, -0.2, 2.0]
d₂ = 145.0

𝒳₁ = generate_planar_points(𝐧₁,d₁, 20, 50)
𝒳₂ = generate_planar_points(𝐧₂,d₂, 20, 50)

world_basis = (Vec(1.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(0.0, 0.0, 1.0))
camera_basis = (Point(0.0, 0.0, 0.0), Vec(-1.0, 0.0, 0.0), Vec(0.0, -1.0, 0.0), Vec(0.0, 0.0, 1.0))
picture_basis = (Point(0.0, 0.0), Vec(-1.0, 0.0), Vec(0.0, -1.0))

camera₁ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)
camera₂ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)
relocate!(camera₁, 𝐑₁, 𝐭₁)
relocate!(camera₂, 𝐑₂, 𝐭₂)

𝐑₁′, 𝐭₁′ = ascertain_pose(camera₁, world_basis... )
𝐊₁′ = obtain_intrinsics(camera₁, CartesianSystem())
𝐑₂′, 𝐭₂′ = ascertain_pose(camera₂, world_basis... )
𝐊₂′ = obtain_intrinsics(camera₂, CartesianSystem())

𝐏₁ = construct(ProjectionMatrix(),𝐊₁′,𝐑₁′,𝐭₁′)
𝐏₂ = construct(ProjectionMatrix(),𝐊₂′,𝐑₂′,𝐭₂′)

# Set of corresponding points.
ℳ₁ = project(camera₁,𝐏₁,𝒳₁)
ℳ₁ʹ= project(camera₂,𝐏₂,𝒳₁)
ℳ₂ = project(camera₁,𝐏₁,𝒳₂)
ℳ₂ʹ= project(camera₂,𝐏₂,𝒳₂)

# 3D points corresponding to the first and second planar surface
𝒴₁ = triangulate(DirectLinearTransform(),𝐏₁,𝐏₂,(ℳ₁,ℳ₁ʹ))
𝒴₂ = triangulate(DirectLinearTransform(),𝐏₁,𝐏₂,(ℳ₂,ℳ₂ʹ))

# Triangulating with the same projection matrices that were used to construct
# (ℳ,ℳʹ) should yield 3D points that lie on a plane.
N = length(𝒴₁)
for n = 1:N
    X = 𝒴₁[n]
    @test isapprox(dot(𝐧₁,X) - d₁, 0.0; atol = 1e-11)
end

N = length(𝒴₂)
for n = 1:N
    X = 𝒴₂[n]
    @test isapprox(dot(𝐧₂,X) - d₂, 0.0; atol = 1e-11)
end

# Triangulating with the same projection matrices that were used to construct
# (ℳ,ℳʹ) should yield the same 3D points as the original 𝒳.
N = length(𝒴₁)
for n = 1:N
    @test  isapprox(sum(abs.(𝒳₁[n]-𝒴₁[n])/3), 0.0; atol = 1e-12)
end

N = length(𝒴₂)
for n = 1:N
    @test  isapprox(sum(abs.(𝒳₂[n]-𝒴₂[n])/3), 0.0; atol = 1e-12)
end


text:
export
dense,
flatten,
ae,
ae_num,
ae_init,
fc,
fc_num,
fc_init,
fc_to_code,
ae_to_code,
fcx,
bn,
sparse_softmax_cross_entropy_with_logits,
Resnet1D

#----------------------------------------------------
# activation functions

_string2fn = Dict(
    "relu" => relu,
    "tanh" => tanh,
    "sigmoid" => sigmoid,
    "leakyrelu" => leaky_relu,
    "leaky_relu" => leaky_relu,
    "relu6" => relu6,
    "softmax" => softmax,
    "softplus" => softplus,
    "selu" => selu,
    "elu" => elu,
    "concat_elu"=>concat_elu,
    "concat_relu"=>concat_relu,
    "hard_sigmoid"=>hard_sigmoid,
    "swish"=>swish,
    "hard_swish"=>hard_swish,
    "concat_hard_swish"=>concat_hard_swish,
    "sin"=>sin,
    "fourier"=>fourier
)

function get_activation_function(a::Union{Function, String, Nothing})
    if isnothing(a)
        return x->x
    end
    if isa(a, String)
        if haskey(_string2fn, lowercase(a))
            return _string2fn[lowercase(a)]
        else
            error("Activation function $a is not supported, available activation functions:\n$(collect(keys(_string2fn))).")
        end
    else
        return a
    end
end

@doc raw"""
    fcx(x::Union{Array{Float64,2},PyObject}, output_dims::Array{Int64,1},
    θ::Union{Array{Float64,1}, PyObject};
    activation::String = "tanh")

Creates a fully connected neural network with output dimension `o` and inputs $x\in \mathbb{R}^{m\times n}$.

$$x \rightarrow o_1 \rightarrow o_2 \rightarrow \ldots \rightarrow o_k$$

`θ` is the weights and biases of the neural network, e.g., `θ = ae_init(output_dims)`.

`fcx` outputs two tensors:

- the output of the neural network: $u\in \mathbb{R}^{m\times o_k}$.

- the sensitivity of the neural network per sample: $\frac{\partial u}{\partial x}\in \mathbb{R}^{m \times o_k \times n}$

"""
function fcx(x::Union{Array{Float64,2},PyObject}, output_dims::Array{Int64,1},
    θ::Union{Array{Float64,1}, PyObject};
    activation::String = "tanh")
    pth = install("ExtendedNN")
    extended_nn_ = load_op_and_grad(pth, "extended_nn"; multiple=true)
    config = [size(x,2);output_dims]
    x_,config_,θ_ = convert_to_tensor([x,config,θ], [Float64,Int64,Float64])
    x_ = reshape(x_, (-1,))
    u, du = extended_nn_(x_,config_,θ_,activation)
    reshape(u, (size(x,1), config[end])), reshape(du, (size(x,1), config[end], size(x,2)))
end

"""
    ae(x::PyObject, output_dims::Array{Int64}, scope::String = "default";
        activation::Union{Function,String} = "tanh")

Alias: `fc`

Creates a neural network with intermediate numbers of neurons `output_dims`.
"""
function ae(x::Union{Array{Float64},PyObject}, output_dims::Array{Int64}, scope::String = "default";
        activation::Union{Function,String} = "tanh")
    if isa(x, Array)
        x = constant(x)
    end
    flag = false
    if length(size(x))==1
        x = reshape(x, length(x), 1)
        flag = true
    end
    net = x
    variable_scope(scope, reuse=AUTO_REUSE) do
        for i = 1:length(output_dims)-1
            net = dense(net, output_dims[i], activation=activation)
        end
        net = dense(net, output_dims[end])
    end
    if flag && (size(net,2)==1)
        net = squeeze(net)
    end
    return net
end

"""
    ae(x::Union{Array{Float64}, PyObject}, output_dims::Array{Int64}, θ::Union{Array{Float64}, PyObject};
    activation::Union{Function,String, Nothing} = "tanh")

Alias: `fc`

Creates a neural network with intermediate numbers of neurons `output_dims`. The weights are given by `θ`

# Example 1: Explicitly construct weights and biases
```julia
x = constant(rand(10,2))
n = ae_num([2,20,20,20,2])
θ = Variable(randn(n)*0.001)
y = ae(x, [20,20,20,2], θ)
```

# Example 2: Implicitly construct weights and biases
```julia
θ = ae_init([10,20,20,20,2])
x = constant(rand(10,10))
y = ae(x, [20,20,20,2], θ)
```

See also [`ae_num`](@ref), [`ae_init`](@ref).
"""
function ae(x::Union{Array{Float64}, PyObject}, output_dims::Array{Int64}, θ::Union{Array{Float64}, PyObject};
    activation::Union{Function,String, Nothing} = "tanh")
    activation = get_activation_function(activation)
    x = convert_to_tensor(x)
    θ = convert_to_tensor(θ)
    flag = false
    if length(size(x))==1
        x = reshape(x, length(x), 1)
        flag = true
    end
    offset = 0
    net = x

    output_dims = [size(x,2); output_dims]
    for i = 1:length(output_dims)-2
        m = output_dims[i]
        n = output_dims[i+1]
        net = net * reshape(θ[offset+1:offset+m*n], (m, n))  + θ[offset+m*n+1:offset+m*n+n]
        net = activation(net)
        offset += m*n+n
    end
    m = output_dims[end-1]
    n = output_dims[end]
    net = net * reshape(θ[offset+1:offset+m*n], (m, n))  + θ[offset+m*n+1:offset+m*n+n]

    offset += m*n+n
    if offset!=length(θ)
        error("ADCME: the weights and configuration does not match. Required $offset but given $(length(θ)).")
    end
    if flag && (size(net,2)==1)
        net = squeeze(net)
    end
    return net
end

"""
    ae(x::Union{Array{Float64}, PyObject},
        output_dims::Array{Int64},
        θ::Union{Array{Array{Float64}}, Array{PyObject}};
        activation::Union{Function,String} = "tanh")

Alias: `fc`

Constructs a neural network with given weights and biases `θ`

# Example
```julia
x = constant(rand(10,30))
θ = ae_init([30, 20, 20, 5])
y = ae(x, [20, 20, 5], θ) # 10×5
```
"""
function ae(x::Union{Array{Float64}, PyObject},
    output_dims::Array{Int64},
    θ::Union{Array{Array{Float64}}, Array{PyObject}};
    activation::Union{Function,String} = "tanh")
    if isa(θ, Array{Array{Float64}})
        val = []
        for t in θ
            push!(val, θ'[:])
        end
        val = vcat(val...)
    else
        val = []
        for t in θ
            push!(val, reshape(θ, (-1,)))
        end
        vcat(val...)
    end
    ae(x, output_dims, θ, activation=activation)
end

@doc raw"""
    ae_init(output_dims::Array{Int64}; T::Type=Float64, method::String="xavier")
    fc_init(output_dims::Array{Int64})

Return the initial weights and bias values by TensorFlow as a vector. The neural network architecture is

$$o_1 (Input layer) \rightarrow o_2 \rightarrow \cdots \rightarrow o_n (Output layer)$$

Three types of
random initializers are provided

- `xavier` (default). It is useful for `tanh` fully connected neural network.
```math
W^l_i \sim \sqrt{\frac{1}{n_{l-1}}}
```
- `xavier_avg`. A variant of `xavier`
```math
W^l_i \sim \sqrt{\frac{2}{n_l + n_{l-1}}}
```
- `he`. This is the activation aware initialization of weights and helps mitigate the problem
of vanishing/exploding gradients.
```math
W^l_i \sim \sqrt{\frac{2}{n_{l-1}}}
```

# Example
```julia
x = constant(rand(10,30))
θ = ae_init([30, 20, 20, 5])
y = ae(x, [20, 20, 5], θ) # 10×5
```
"""
function ae_init(output_dims::Array{Int64}; T::Type=Float64, method::String="xavier")
    N = ae_num(output_dims)
    W = zeros(T, N)
    offset = 0
    for i = 1:length(output_dims)-1
        m = output_dims[i]
        n = output_dims[i+1]
        if method=="xavier"
            W[offset+1:offset+m*n] = randn(T, m*n) * T(sqrt(1/m))
        elseif method=="xavier_normal"
            W[offset+1:offset+m*n] = randn(T, m*n) * T(sqrt(2/(n+m)))
        elseif method=="xavier_uniform"
            W[offset+1:offset+m*n] = rand(T, m*n) * T(sqrt(6/(n+m)))
        elseif method=="he"
            W[offset+1:offset+m*n] = randn(T, m*n) * T(sqrt(2/(m)))
        else
            error("Method $method not understood")
        end
        offset += m*n+n
    end
    W
end

"""
    ae_num(output_dims::Array{Int64})
    fc_num(output_dims::Array{Int64})

Estimates the number of weights and biases for the neural network. Note the first dimension
should be the feature dimension (this is different from [`ae`](@ref) since in `ae` the feature
dimension can be inferred), and the last dimension should be the output dimension.

# Example
```julia
x = constant(rand(10,30))
θ = ae_init([30, 20, 20, 5])
@assert ae_num([30, 20, 20, 5])==length(θ)
y = ae(x, [20, 20, 5], θ)
```
"""
function ae_num(output_dims::Array{Int64})
    offset = 0
    for i = 1:length(output_dims)-2
        m = output_dims[i]
        n = output_dims[i+1]
        offset += m*n+n
    end
    m = output_dims[end-1]
    n = output_dims[end]
    offset += m*n+n
    return offset
end

function _ae_to_code(d::Dict, scope::String; activation::String)
    i = 0
    nn_code = ""
    while true
        si = i==0 ? "" : "_$i"
        Wkey = "$(scope)backslashfully_connected$(si)backslashweightscolon0"
        bkey = "$(scope)backslashfully_connected$(si)backslashbiasescolon0"
        if haskey(d, Wkey)
            if i!=0
                nn_code *= "    isa(net, Array) ? (net = $activation.(net)) : (net=$activation(net))\n"
                nn_code *= "    #-------------------------------------------------------------------\n"
            end
            nn_code *= "    W$i = aedict$scope[\"$Wkey\"]\n    b$i = aedict$scope[\"$bkey\"];\n"
            nn_code *= "    isa(net, Array) ? (net = net * W$i .+ b$i') : (net = net *W$i + b$i)\n"
            i += 1
        else
            break
        end
    end
    nn_code = """  global nn$scope\n  function nn$scope(net)
$(nn_code)    return net\n  end """
    nn_code
end

"""
    ae_to_code(file::String, scope::String; activation::String = "tanh")

Return the code string from the feed-forward neural network data in `file`. Usually we can immediately evaluate
the code string into Julia session by
```julia
eval(Meta.parse(s))
```
If `activation` is not specified, `tanh` is the default.
"""
function ae_to_code(file::String, scope::String = "default"; activation::String = "tanh")
    d = matread(file)
    s = "let aedict$scope = matread(\"$file\")\n"*_ae_to_code(d, scope; activation = activation)*"\nend\n"
    return s
end

fc_to_code = ae_to_code

# tensorflow layers from contrib
for (op, tfop) in [(:avg_pool2d, :avg_pool2d), (:avg_pool3d, :avg_pool3d),
        (:flatten, :flatten), (:max_pool2d, :max_pool2d), (:max_pool3d, :max_pool3d)]
    @eval begin
        export $op
        $op(args...; kwargs...) = tf.contrib.layers.$tfop(args...; kwargs...)
    end
end

for (op, tfop) in [(:conv1d, :conv1d), (:conv2d, :conv2d), (:conv2d_in_plane, :conv2d_in_plane),
    (:conv2d_transpose, :conv2d_transpose), (:conv3d, :conv3d), (:conv3d_transpose, :conv3d_transpose)]
    @eval begin
        export $op
        function $op(args...;activation = nothing, bias=false,  kwargs...)
            activation = get_activation_function(activation)
            if bias
                biases_initializer = tf.zeros_initializer()
            else
                biases_initializer = nothing
            end
            tf.contrib.layers.$tfop(args...; activation_fn = activation, biases_initializer=biases_initializer, kwargs...)
        end
    end
end

"""
    dense(inputs::Union{PyObject, Array{<:Real}}, units::Int64, args...;
        activation::Union{String, Function} = nothing, kwargs...)

Creates a fully connected layer with the activation function specified by `activation`
"""
function dense(inputs::Union{PyObject, Array{<:Real}}, units::Int64, args...; activation::Union{String, Function, Nothing} = nothing, kwargs...)
    inputs = constant(inputs)
    activation = get_activation_function(activation)
    tf.contrib.layers.fully_connected(inputs, units,  args...; activation_fn=activation, kwargs...)
end

"""
    bn(args...;center = true, scale=true, kwargs...)

`bn` accepts a keyword parameter `is_training`.
# Example
```julia
bn(inputs, name="batch_norm", is_training=true)
```

!!! note
    `bn` should be used with `control_dependency`
    ```julia
    update_ops = get_collection(UPDATE_OPS)
    control_dependencies(update_ops) do
        global train_step = AdamOptimizer().minimize(loss)
    end
    ```
"""
function bn(args...;center = true, scale=true, kwargs...)
    @warn  """
`bn` should be used with `control_dependency`
```julia
update_ops = get_collection(UPDATE_OPS)
control_dependencies(update_ops) do
    global train_step = AdamOptimizer().minimize(loss)
end
```
""" maxlog=1
    kwargs = Dict{Any, Any}(kwargs)
    if :is_training in keys(kwargs)
        kwargs[:training] = kwargs[:is_training]
        delete!(kwargs, :is_training)
    end
    if :scope in keys(kwargs)
        kwargs[:name] = kwargs[:scope]
        delete!(kwargs, :scope)
    end
    tf.layers.batch_normalization(args...;center = center, scale=scale, kwargs...)
end
sparse_softmax_cross_entropy_with_logits(args...;kwargs...) = tf.nn.sparse_softmax_cross_entropy_with_logits(args...;kwargs...)

export group_conv2d
function group_conv2d(inputs::PyObject, filters::Int64, args...;  groups = 1, scope=scope, kwargs...)
    if groups==1
        return conv2d(inputs, filters, args...;scope=scope, kwargs...)
    else
        dims = size(inputs)
        if mod(dims[end], groups)!=0 || mod(filters, groups)!=0
            error("channels and outputs must be the multiples of `groups`")
        end
        n = div(filters, groups)
        in_ = Array{PyObject}(undef, groups)
        out_ = Array{PyObject}(undef, groups)
        for i = 1:groups
py"""
temp = $inputs[:,:,:,$((i-1)*n):$(i*n)]
"""
            in_[i] = py"temp"
            out_[i] = conv2d(in_[i], n, args...;scope=scope*"_group$i", kwargs...)
        end
        out = concat(out_, dims=4)
        return out
    end
end

export separable_conv2d
function separable_conv2d(inputs, num_outputs, args...; activation=nothing, bias=false, kwargs...)
    activation = get_activation_function(activation)
    if bias
        biases_initializer = tf.zeros_initializer()
    else
        biases_initializer = nothing
    end
    tf.contrib.layers.separable_conv2d(inputs, num_outputs, args...;activation_fn = activation,
                        biases_initializer=biases_initializer, kwargs...)
end

export depthwise_conv2d
function depthwise_conv2d(input, num_outputs, args...;
    kernel_size = 3, channel_multiplier = 1, stride = 1, padding="SAME", bias=false, scope="default", reuse=AUTO_REUSE, kwargs...)
    local res, strides
    if isa(kernel_size, Int64)
        kernel_size = (3,3)
    end
    if isa(stride, Int64)
        strides = (1, stride, stride, 1)
    end
    variable_scope(scope, reuse=reuse) do
        filter = get_variable("dconv2d",
                shape=[kernel_size[1], kernel_size[2], size(input,4), channel_multiplier],
                initializer=tf.contrib.layers.xavier_initializer())
        res = tf.nn.depthwise_conv2d(input, filter, strides, padding, args...; kwargs...)
        if bias
            res = tf.contrib.layers.bias_add(res, scope=scope, reuse=AUTO_REUSE)
        end
    end
    return res
end

"""
$(@doc ae)
"""
fc = ae
"""
$(@doc ae_num)
"""
fc_num = ae_num
"""
$(@doc ae_init)
"""
fc_init = ae_init

#------------------------------------------------------------------------------------------
export dropout

"""
    dropout(x::Union{PyObject, Real, Array{<:Real}},
    rate::Union{Real, PyObject}, training::Union{PyObject,Bool} = true; kwargs...)

Randomly drops out entries in `x` with a rate of `rate`.
"""
function dropout(x::Union{PyObject, Real, Array{<:Real}},
    rate::Union{Real, PyObject}, training::Union{PyObject,Bool, Nothing} = nothing ; kwargs...)
    x = constant(x)
    if isnothing(training)
        training = options.training.training
    else
        training = constant(training)
    end
    tf.keras.layers.Dropout(rate, kwargs...)(x, training)
end

export BatchNormalization
mutable struct BatchNormalization
    dims::Int64
    o::PyObject
end

"""
    BatchNormalization(dims::Int64=2; kwargs...)

Creates a batch normalization layer.
# Example
```julia
b = BatchNormalization(2)
x = rand(10,2)
training = placeholder(true)
y = b(x, training)
run(sess, y)
```
"""
function BatchNormalization(dims::Int64=-1; kwargs...)
    local o
    if dims>=1
        o = tf.keras.layers.BatchNormalization(dims-1, kwargs...)
    else
        o = tf.keras.layers.BatchNormalization(dims, kwargs...)
    end
    BatchNormalization(dims, o)
end

function Base.:show(io::IO, b::BatchNormalization)
    print("<BatchNormalization normalization_dim=$(b.dims)>")
end

function (o::BatchNormalization)(x, training=ADCME.options.training.training)
    flag = false
    if get_dtype(x)==Float64
        x = cast(x, Float32)
        flag = true
    end
    out = o.o(x, training)
    if flag
        out = cast(Float64, out)
    end
    return out
end

export Dense, Conv1D, Conv2D, Conv3D, Conv2DTranspose

mutable struct Dense
    hidden_dim::Int64
    activation::Union{String, Function}
    o::PyObject
end

"""
    Dense(units::Int64, activation::Union{String, Function, Nothing} = nothing,
        args...;kwargs...)

Creates a callable dense neural network.
"""
function Dense(units::Int64, activation::Union{String, Function, Nothing} = nothing,
    args...;kwargs...)
    activation = get_activation_function(activation)
    o = tf.keras.layers.Dense(units, activation, args...;kwargs...)
    Dense(units, activation, o)
end

function Base.:show(io::IO, o::Dense)
    print("<Fully connected neural network with $(o.hidden_dim) hidden units and activation function \"$(o.activation)\">")
end

(o::Dense)(x) = o.o(x)

mutable struct Conv1D
    filters
    kernel_size
    strides
    activation
    o::PyObject
end
"""
    Conv1D(filters, kernel_size, strides, activation, args...;kwargs...)

```julia
c = Conv1D(32, 3, 1, "relu")
x = rand(100, 6, 128) # 128-length vectors with 6 timesteps ("channels")
y = c(x) # shape=(100, 4, 32)
```
"""
function Conv1D(filters, kernel_size, strides=1, activation=nothing, args...;kwargs...)
    activation = get_activation_function(activation)
    o = tf.keras.layers.Conv1D(filters, kernel_size, strides,args...; activation = activation, kwargs...)
    Conv1D(filters, kernel_size, strides, activation, o)
end
function Base.:show(io::IO, o::Conv1D)
    print("<Conv1D filters=$(o.filters) kernel_size=$(o.kernel_size) strides=$(o.strides) activation=$(o.activation)>")
end
function (o::Conv1D)(x::Union{PyObject, Array{<:Real,3}})
    x = constant(x)
    @assert length(size(x))==3
    o.o(x)
end

mutable struct Conv2D
    filters
    kernel_size
    strides
    activation
    o::PyObject
end
"""
    Conv2D(filters, kernel_size, strides, activation, args...;kwargs...)

The arrangement is (samples, rows, cols, channels) (data_format='channels_last')
```julia
Conv2D(32, 3, 1, "relu")
```
"""
function Conv2D(filters, kernel_size, strides=1, activation=nothing, args...;kwargs...)
    activation = get_activation_function(activation)
    o = tf.keras.layers.Conv2D(filters, kernel_size, strides,args...; activation = activation, kwargs...)
    Conv2D(filters, kernel_size, strides, activation, o)
end
function Base.:show(io::IO, o::Conv2D)
    print("<Conv2D filters=$(o.filters) kernel_size=$(o.kernel_size) strides=$(o.strides) activation=$(o.activation)>")
end
function (o::Conv2D)(x::Union{PyObject, Array{<:Real,4}})
    x = constant(x)
    @assert length(size(x))==4
    o.o(x)
end

mutable struct Conv2DTranspose
    filters
    kernel_size
    strides
    activation
    o::PyObject
end

function Conv2DTranspose(filters, kernel_size, strides=1, activation=nothing, args...;kwargs...)
    activation = get_activation_function(activation)
    o = tf.keras.layers.Conv2DTranspose(
        filters, kernel_size, strides=strides; activation = activation, kwargs...
    )
    Conv2DTranspose(filters, kernel_size, strides, activation, o)
end

function (o::Conv2DTranspose)(x::Union{PyObject, Array{<:Real,4}})
    x = constant(x)
    @assert length(size(x))==4
    o.o(x)
end

mutable struct Conv3D
    filters
    kernel_size
    strides
    activation
    o::PyObject
end
"""
    Conv3D(filters, kernel_size, strides, activation, args...;kwargs...)

The arrangement is (samples, rows, cols, channels) (data_format='channels_last')
```julia
c = Conv3D(32, 3, 1, "relu")
x = constant(rand(100, 10, 10, 10, 16))
y = c(x)
# shape=(100, 8, 8, 8, 32)
```
"""
function Conv3D(filters, kernel_size, strides=1, activation=nothing, args...;kwargs...)
    activation = get_activation_function(activation)
    o = tf.keras.layers.Conv3D(filters, kernel_size, strides,args...; activation = activation, kwargs...)
    Conv3D(filters, kernel_size, strides, activation, o)
end
function Base.:show(io::IO, o::Conv3D)
    print("<Conv3D filters=$(o.filters) kernel_size=$(o.kernel_size) strides=$(o.strides) activation=$(o.activation)>")
end
function (o::Conv3D)(x::Union{PyObject, Array{<:Real,5}})
    x = constant(x)
    @assert length(size(x))==5
    o.o(x)
end

#------------------------------------------------------------------------------------------
# resnet, adapted from https://github.com/bayesiains/nsf/blob/master/nn/resnet.py
mutable struct ResnetBlock
    use_batch_norm::Bool
    activation::Union{String,Function}
    linear_layers::Array{Dense}
    bn_layers::Array{BatchNormalization}
    dropout_probability::Float64
end

function ResnetBlock(features::Int64; dropout_probability::Float64=0., use_batch_norm::Bool,
     activation::String="relu")
    activation = get_activation_function(activation)
    bn_layers = []
    if use_batch_norm
        bn_layers = [
            BatchNormalization(),
            BatchNormalization()
        ]
    end
    linear_layers = [
        Dense(features, nothing)
        Dense(features, nothing)
    ]
    ResnetBlock(use_batch_norm, activation, linear_layers, bn_layers, dropout_probability)
end

function (res::ResnetBlock)(input)
    x = input
    if res.use_batch_norm
        x = res.bn_layers[1](x)
    end
    x = res.activation(x)
    x = res.linear_layers[1](x)
    if res.use_batch_norm
        x = res.bn_layers[2](x)
    end
    x = res.activation(x)
    x = dropout(x, res.dropout_probability, options.training.training)
    x = res.linear_layers[2](x)
    return x + input
end

mutable struct Resnet1D
    initial_layer::Dense
    blocks::Array{ResnetBlock}
    final_layer::Dense
end

"""
    Resnet1D(out_features::Int64, hidden_features::Int64;
        num_blocks::Int64=2, activation::Union{String, Function, Nothing} = "relu",
        dropout_probability::Float64 = 0.0, use_batch_norm::Bool = false, name::Union{String, Missing} = missing)

Creates a 1D residual network. If `name` is not missing, `Resnet1D` does not create a new entity.
# Example
```julia
resnet = Resnet1D(20)
x = rand(1000,10)
y = resnet(x)
```

# Example: Digit recognition
```
using MLDatasets
using ADCME

# load data
train_x, train_y = MNIST.traindata()
train_x = reshape(Float64.(train_x), :, size(train_x,3))'|>Array
test_x, test_y = MNIST.testdata()
test_x = reshape(Float64.(test_x), :, size(test_x,3))'|>Array

# construct loss function
ADCME.options.training.training = placeholder(true)
x = placeholder(rand(64, 784))
l = placeholder(rand(Int64, 64))
resnet = Resnet1D(10, num_blocks=10)
y = resnet(x)
loss = mean(sparse_softmax_cross_entropy_with_logits(labels=l, logits=y))

# train the neural network
opt = AdamOptimizer().minimize(loss)
sess = Session(); init(sess)
for i = 1:10000
    idx = rand(1:60000, 64)
    _, loss_ = run(sess, [opt, loss], feed_dict=Dict(l=>train_y[idx], x=>train_x[idx,:]))
    @info i, loss_
end

# test
for i = 1:10
    idx = rand(1:10000,100)
    y0 = resnet(test_x[idx,:])
    y0 = run(sess, y0, ADCME.options.training.training=>false)
    pred = [x[2]-1 for x in argmax(y0, dims=2)]
    @info "Accuracy = ", sum(pred .== test_y[idx])/100
end
```
![](https://github.com/ADCMEMarket/ADCMEImages/tree/master/ADCME/assets/resnet.png?raw=true)
"""
function Resnet1D(out_features::Int64, hidden_features::Int64 = 20;
     num_blocks::Int64=2, activation::Union{String, Function, Nothing} = "relu",
     dropout_probability::Float64 = 0.0, use_batch_norm::Bool = false,
     name::Union{String, Missing} = missing)
     if haskey(ADCME.STORAGE, name)
        @info "Reusing $name..."
        return ADCME.STORAGE[name]
     end
     initial_layer = Dense(hidden_features)
     blocks = ResnetBlock[]
     for i = 1:num_blocks
        push!(blocks, ResnetBlock(
            hidden_features,
            dropout_probability = dropout_probability,
            use_batch_norm = use_batch_norm
        ))
     end
     final_layer = Dense(out_features)
     if ismissing(name)
        name = "Resnet1D_"*randstring(10)
     end
     res = Resnet1D(initial_layer, blocks, final_layer)
     ADCME.STORAGE[name] = res
     return res
end

function (res::Resnet1D)(x)
    x = res.initial_layer(x)
    for b in res.blocks
        x = b(x)
    end
    x = res.final_layer(x)
    x
end

function Base.:show(io::IO, res::Resnet1D)
    println("( Input )")
    println("\t↓")
    show(io, res.initial_layer)
    println("\n")
    for i = 1:length(res.blocks)
        show(io, res.blocks[i])
        println("\n")
        println("\t↓")
    end
    show(io, res.final_layer)
    println("\n( Output )")
end


text:
using JuliaWeBWorK

qs = JuliaWeBWorK.QUESTIONS()
## --------------------------------------------------

meta=(Project="Some sample statistics questions",
     Question="3"
      )

intro = jmt"""

Here are two sample   statistics questions.

"""

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

# need raw, dollar sign
numericq(raw"""

### Problem 1

From the  following `R`  output, what is a  95% confidence interval?

```
t.test(Data$Steps, conf.level=0.95)

One Sample t-test

95 percent confidence interval:
 7171.667 8212.949

mean of x
 7692.308
```
""",
         () -> begin
         List([Interval(7171.667, 8212.949)])
         end,
         ()
         ) |> qs

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

r = randomizer(20:25) |> qs

numericq(raw"""

### Problem 2

Output from an  `R` command is

```
    One Sample t-test
data:  my_data$weight
t = -9.0783, df = 9, p-value = 7.953e-06
alternative hypothesis: true mean is not equal to {{:a1}}
95 percent confidence interval:
 17.8172 20.6828
sample estimates:
mean of x
    19.25
```

(a) What is the observed value of the statistic?
""",
         (μ) ->  -9.0783,
         r) |> qs

numericq(jmt"""

(b) What was the sample  size?
""",
         (μ) -> 9 + 1,
         r) |> qs

numericq(jmt"""

(c)  Using  the fact that

```math
t = \frac{\bar{x} - \mu}{SE},
```

Find \( SE \).
""",
         (μ) -> begin
         t, x̄ = -9.0783, 19.25
         (x̄ -μ)/t
         end,
         r) |> qs

## --------------------------------------------------
p  =  Page(intro, qs; context="Interval",  meta...)


text:
mutable struct SmallLLLRelationsCtx
  A::NfOrdIdl
  b::Array{nf_elem, 1}
  bd::Int
  cnt::Int
  elt::nf_elem
  function SmallLLLRelationsCtx()
    n = new()
    n.bd = 1
    n.cnt = 0
    return n
  end
end

function class_group_small_lll_elements_relation_start{T}(clg::ClassGrpCtx{T},
                O::NfOrd; prec::Int = 200, val::Int = 0,
                limit::Int = 0)
  return class_group_small_lll_elements_relation_start(clg, hecke.ideal(O, parent(basis_mat(O).num)(1)), prec = prec)
end

function class_group_small_lll_elements_relation_start{T}(clg::ClassGrpCtx{T},
                A::NfOrdIdl; prec::Int = 200, val::Int = 0,
                limit::Int = 0)
  global _start
  K = nf(order(A))
  local rt::UInt
  @v_do :ClassGroup_time 2 rt = time_ns()

  local I, S, bd

  while true
    try
      L::FakeFmpqMat, Tr::fmpz_mat = lll(A, prec = prec)
      @v_do :ClassGroup_time 2 _start += time_ns()-rt
      I = SmallLLLRelationsCtx()
      S::FakeFmpqMat = FakeFmpqMat(Tr)*basis_mat(A, Val{false})*basis_mat(order(A), Val{false})
      bd::fmpz = abs(discriminant(order(A)))*norm(A)^2
      bd = root(bd, degree(K))::fmpz
      bd *= L.den
      f = find(i-> cmpindex(L.num, i, i, bd) < 0, 1:degree(K))
      m = div(degree(K), 4)
      if m < 2
        m = degree(K)
      end
      while length(f) < m
        f = find(i-> cmpindex(L.num, i, i, bd) < 0, 1:degree(K))
        bd *= 2
      end
      I.b = nf_elem[]
      for i=f
        push!(I.b, elem_from_mat_row(K, S.num, i, S.den))
      end
      #println([Float64(numerator(L)[i,i]//denominator(L)*1.0) for i=1:degree(K)])
      #now select a subset that can yield "small" relations, where
      #small means of effective norm <= sqrt(disc)
      I.A = A
      I.elt = K()
      return I
    catch e
      if isa(e, LowPrecisionLLL)
        print_with_color(:red, "prec too low in LLL,")
        prec = Int(ceil(1.2*prec))
#        println(" increasing to ", prec)
        if prec > 1000
          error("2:too much prec")
        end
      else
        rethrow(e)
      end
    end
  end
end

function class_group_small_lll_elements_relation_next(I::SmallLLLRelationsCtx)
  #the 1st n relations are the basis elements
  if I.cnt < length(I.b)
    I.cnt += 1
    I.elt = I.b[I.cnt]
    return deepcopy(I.b[I.cnt])
  end
  #the next 2*n*(n-1)/2 relations are the ones of weight 2
  #1st b[i] + b[j], all combinations, then b[i]-b[j]
  #it may (or may not) be helpful
  if I.cnt <= length(I.b)^2
    n = length(I.b)
    c = I.cnt - n # we already did n relations in the 1st if
    if c > n*(n-1)/2
      c -= div(n*(n-1), 2)
      s = -1
    else
      s = 1
    end
    i = 1
    while c+i-1 >= n
      c -= (n-i)
      i += 1
    end
    I.cnt += 1
    I.elt = I.b[i] + s*I.b[c+i]
    return I.elt
  end

  if I.cnt > (2*I.bd+1)^div(length(I.b), 2)
    I.bd += 1
  end

  I.cnt += 1
  while true
    rand!(I.elt, I.b, -I.bd:I.bd, min(length(I.b), 5))
    if !iszero(I.elt)
      return deepcopy(I.elt)
    end
  end
end


text:
module RulesBetaMarginalsTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_marginalrules

@testset "marginalrules:Beta" begin

     @testset "out_a_b: (m_out::Beta, m_a::PointMass, m_b::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] Beta(:out_a_b) [
            (input = (m_out = Beta(1.0, 2.0), m_a = PointMass(1.0), m_b = PointMass(2.0)), output = (out = Beta(1.0, 3.0), a = PointMass(1.0), b = PointMass(2.0))),
            (input = (m_out = Beta(2.0, 2.0), m_a = PointMass(2.0), m_b = PointMass(3.0)), output = (out = Beta(3.0, 4.0), a = PointMass(2.0), b = PointMass(3.0))),
            (input = (m_out = Beta(2.0, 3.0), m_a = PointMass(1.0), m_b = PointMass(3.0)), output = (out = Beta(2.0, 5.0), a = PointMass(1.0), b = PointMass(3.0))),
        ]

    end
end
end


text:
using SafeTestsets

@safetestset "Common Integrator Functionality                                                 " begin include("integrators_common_tests.jl") end
@safetestset "Extrapolation Methods                                                           " begin include("extrapolation_tests.jl") end
@safetestset "Initial Guesses                                                                 " begin include("initial_guess_tests.jl") end
@safetestset "Runge-Kutta Integrators                                                         " begin include("rk_integrators_tests.jl") end
@safetestset "VPRK Integrators                                                                " begin include("vprk_integrators_tests.jl") end
@safetestset "SPARK Integrators                                                               " begin include("spark_integrators_tests.jl") end
@safetestset "Degenerate Variational Integrators                                              " begin include("dvi_integrators_tests.jl") end
@safetestset "Special Integrators                                                             " begin include("special_integrators_tests.jl") end
@safetestset "Splitting Integrators                                                           " begin include("splitting_integrators_tests.jl") end
@safetestset "Galerkin Integrators                                                            " begin include("galerkin_integrators_tests.jl") end


text:
module TestNaturalSelection

using Base.Test

using BioSequences, NaturalSelection

include("dNdS.jl")
include("tajima.jl")
include("mkt.jl")

end


text:
#=
  Structs for the approxChol solver
=#

"""
  LLp elements are all in the same column.
  row tells us the row, and val is the entry.
  val is set to zero for some edges that we should remove.
  next gives the next in the column.  It points to itself to terminate.
  reverse is the index into lles of the other copy of this edge,
  since every edge is stored twice as we do not know the order of elimination in advance.
"""
mutable struct LLp{Tind,Tval}
    row::Tind
    val::Tval
    next::LLp{Tind,Tval}
    reverse::LLp{Tind,Tval}

    LLp{Tind,Tval}() where {Tind,Tval} = (x = new(zero(Tind), zero(Tval)); x.next = x; x.reverse = x)
    LLp{Tind,Tval}(row, val, next, rev) where {Tind,Tval} = new(row, val, next, rev)
    LLp{Tind,Tval}(row, val) where {Tind,Tval} = (x = new(row, val); x.next = x; x.reverse = x)
    LLp{Tind,Tval}(row, val, next) where {Tind,Tval} = (x = new(row, val, next); x.reverse = x)
end

"""
  LLmatp is the data structure used to maintain the matrix during elimination.
  It stores the elements in each column in a singly linked list (only next ptrs)
  Each element is an LLp (linked list pointer).
  The head of each column is pointed to by cols.

  We probably can get rid of degs - as it is only used to store initial degrees.
"""
mutable struct LLmatp{Tind,Tval}
    n::Int64
    degs::Array{Tind,1}
    cols::Array{LLp{Tind,Tval},1}
    lles::Array{LLp{Tind,Tval},1}
end

# these are the types we use with a fixed ordering
struct LLord{Tind,Tval}
    row::Tind
    next::Tind
    val::Tval
end

mutable struct LLMatOrd{Tind,Tval}
    n::Int64
    cols::Array{Tind,1}
    lles::Array{LLord{Tind,Tval},1}
end

struct LLcol{Tind,Tval}
      row::Tind
      ptr::Tind
      val::Tval
  end

#=============================================================

LDLinv

=============================================================#

"""
  LDLinv contains the information needed to solve the Laplacian systems.
  It does it by applying Linv, then Dinv, then Linv (transpose).
  But, it is specially constructed for this particular solver.
  It does not explicitly make the matrix triangular.
  Rather, col[i] is the name of the ith col to be eliminated
"""
mutable struct LDLinv{Tind,Tval}
    col::Array{Tind,1}
    colptr::Array{Tind,1}
    rowval::Array{Tind,1}
    fval::Array{Tval,1}
    d::Array{Tval,1}
end

#=============================================================

ApproxCholPQ
the data strcture we use to keep track of degrees

=============================================================#

struct ApproxCholPQElem{Tind}
    prev::Tind
    next::Tind
    key::Tind
end

"""
  An approximate priority queue.
  Items are bundled together into doubly-linked lists with all approximately the same key.
  minlist is the min list we know to be non-empty.
  It should always be a lower bound.
  keyMap maps keys to lists
"""
mutable struct ApproxCholPQ{Tind}
    elems::Array{ApproxCholPQElem{Tind},1} # indexed by node name
    lists::Array{Tind,1}
    minlist::Int
    nitems::Int
    n::Int
end


text:
# Metagraphs files are simply JLD2 files.

struct MGFormat <: AbstractGraphFormat end
struct DOTFormat <: AbstractGraphFormat end

function loadmg(fn::AbstractString)
    @load fn g
    return g
end

function savemg(fn::AbstractString, g::AbstractMetaGraph)
    @save fn g
    return 1
end

loadgraph(fn::AbstractString, gname::String, ::MGFormat) = loadmg(fn)
savegraph(fn::AbstractString, g::AbstractMetaGraph) =  savemg(fn, g)

function savedot(io::IO, g::MetaDiGraph)
    write(io, "digraph G {\n")
    for p in props(g)
        write(io, "$(p[1])=$(p[2]);\n")
    end

    for v in vertices(g)
        write(io, "$v")
        if length(props(g, v)) > 0
            write(io, " [ ")
        end
        for p in props(g, v)
            key = p[1]
            write(io, "$key=\"$(p[2])\",")
        end
        if length(props(g, v)) > 0
            write(io, "];")
        end
        write(io, "\n")
    end

    for e in edges(g)
        write(io, "$(src(e)) -> $(dst(e)) [ ")
        for p in props(g,e)
            write(io, "$(p[1])=$(p[2]), ")
        end
        write(io, "]\n")
    end
    write(io, "}\n")
end

function savegraph(fn::AbstractString, g::AbstractMetaGraph, ::DOTFormat)
    open(fn, "w") do fp
        savedot(fp, g)
    end
end