text:
export unsafe_delete!,
       personality, personality!,
       callconv, callconv!,
       gc, gc!, intrinsic_id,
       entry

# forward declaration of Function in src/core/basicblock.jl
identify(::Type{Value}, ::Val{API.LLVMFunctionValueKind}) = Function

Function(mod::Module, name::String, ft::FunctionType) =
    Function(API.LLVMAddFunction(ref(mod), name, ref(ft)))

unsafe_delete!(::Module, f::Function) = API.LLVMDeleteFunction(ref(f))

function personality(f::Function)
    has_personality = convert(Core.Bool, API.LLVMHasPersonalityFn(ref(f)))
    return has_personality ? Function(API.LLVMGetPersonalityFn(ref(f))) : nothing
end
personality!(f::Function, persfn::Union{Nothing,Function}) =
    API.LLVMSetPersonalityFn(ref(f), persfn===nothing ? C_NULL : ref(persfn))

intrinsic_id(f::Function) = API.LLVMGetIntrinsicID(ref(f))

callconv(f::Function) = API.LLVMGetFunctionCallConv(ref(f))
callconv!(f::Function, cc) = API.LLVMSetFunctionCallConv(ref(f), cc)

function gc(f::Function)
  ptr = API.LLVMGetGC(ref(f))
  return ptr==C_NULL ? "" :  unsafe_string(ptr)
end
gc!(f::Function, name::String) = API.LLVMSetGC(ref(f), name)

entry(f::Function) = BasicBlock(API.LLVMGetEntryBasicBlock(ref(f)))

# attributes

export function_attributes, parameter_attributes, return_attributes

struct FunctionAttrSet
    f::Function
    idx::API.LLVMAttributeIndex
end

function_attributes(f::Function) = FunctionAttrSet(f, API.LLVMAttributeFunctionIndex)
parameter_attributes(f::Function, idx::Integer) = FunctionAttrSet(f, API.LLVMAttributeIndex(idx))
return_attributes(f::Function) = FunctionAttrSet(f, API.LLVMAttributeReturnIndex)

Base.eltype(::FunctionAttrSet) = Attribute

function Base.collect(iter::FunctionAttrSet)
    elems = Vector{API.LLVMAttributeRef}(undef, length(iter))
    if length(iter) > 0
      # FIXME: this prevents a nullptr ref in LLVM similar to D26392
      API.LLVMGetAttributesAtIndex(ref(iter.f), iter.idx, elems)
    end
    return Attribute[Attribute(elem) for elem in elems]
end

Base.push!(iter::FunctionAttrSet, attr::Attribute) =
    API.LLVMAddAttributeAtIndex(ref(iter.f), iter.idx, ref(attr))

Base.delete!(iter::FunctionAttrSet, attr::EnumAttribute) =
    API.LLVMRemoveEnumAttributeAtIndex(ref(iter.f), iter.idx, kind(attr))

function Base.delete!(iter::FunctionAttrSet, attr::StringAttribute)
    k = kind(attr)
    API.LLVMRemoveStringAttributeAtIndex(ref(iter.f), iter.idx, k, length(k))
end

function Base.length(iter::FunctionAttrSet)
    API.LLVMGetAttributeCountAtIndex(ref(iter.f), iter.idx)
end

# parameter iteration

export Argument, parameters

@checked struct Argument <: Value
    ref::reftype(Value)
end
identify(::Type{Value}, ::Val{API.LLVMArgumentValueKind}) = Argument

struct FunctionParameterSet
    f::Function
end

parameters(f::Function) = FunctionParameterSet(f)

Base.eltype(::FunctionParameterSet) = Argument

function Base.getindex(iter::FunctionParameterSet, i)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return Argument(API.LLVMGetParam(ref(iter.f), i-1))
end

function Base.iterate(iter::FunctionParameterSet, state=API.LLVMGetFirstParam(ref(iter.f)))
    state == C_NULL ? nothing : (Argument(state), API.LLVMGetNextParam(state))
end

Base.last(iter::FunctionParameterSet) =
    Argument(API.LLVMGetLastParam(ref(iter.f)))

Base.isempty(iter::FunctionParameterSet) =
    API.LLVMGetLastParam(ref(iter.f)) == C_NULL

Base.length(iter::FunctionParameterSet) = API.LLVMCountParams(ref(iter.f))

# NOTE: optimized `collect`
function Base.collect(iter::FunctionParameterSet)
    elems = Vector{API.LLVMValueRef}(undef, length(iter))
    API.LLVMGetParams(ref(iter.f), elems)
    return map(el->Argument(el), elems)
end

# basic block iteration

export blocks

struct FunctionBlockSet
    f::Function
end

blocks(f::Function) = FunctionBlockSet(f)

Base.eltype(::FunctionBlockSet) = BasicBlock

function Base.iterate(iter::FunctionBlockSet, state=API.LLVMGetFirstBasicBlock(ref(iter.f)))
    state == C_NULL ? nothing : (BasicBlock(state), API.LLVMGetNextBasicBlock(state))
end

Base.last(iter::FunctionBlockSet) =
    BasicBlock(API.LLVMGetLastBasicBlock(ref(iter.f)))

Base.isempty(iter::FunctionBlockSet) =
    API.LLVMGetLastBasicBlock(ref(iter.f)) == C_NULL

Base.length(iter::FunctionBlockSet) = API.LLVMCountBasicBlocks(ref(iter.f))

# NOTE: optimized `collect`
function Base.collect(iter::FunctionBlockSet)
    elems = Vector{API.LLVMValueRef}(undef, length(iter))
    API.LLVMGetBasicBlocks(ref(iter.f), elems)
    return BasicBlock[BasicBlock(elem) for elem in elems]
end


text:
using LasIO
using LazIO
using FileIO
using DataFrames
using ProgressMeter
using StaticArrays
using BenchmarkTools
using PointCloudRasterizers
using NearestNeighbors #use it for spatial search

#get the pulse
function classify(points)
    prev_point = points[1]
    pulse =  Dict{Int,String}() #dictionary with indexes and tags
    n = 0
    for (i,p) in enumerate(points)
        point = p
        inten = intensity(p)
        rn = return_number(p)
        nr = number_of_returns(p)
        dt = gps_time(point) - gps_time(prev_point)
        #rn = return number, nr = number_of_returns
        if (rn == n + 1) && (nr > 1)
            if rn == 1
                pulse[i] = "firstpoint"
            elseif rn < nr && dt < 1e-7
                pulse[i] = "nextpoint"
            elseif rn == nr && dt < 1e-7
                pulse[i] = "lastpoint"
            end
            if nr == rn
                n = 0
            end
        elseif nr == 1 && rn == 1
            pulse[i] = "firstpoint"
            n = 0
        else
            pulse[i] ="unclassified"
            n = 0
        end
        n+=1
        prev_point = p
    end
    pulse
end

#write point cloud
function write_output(writefiles,ds, pulse)
    LazIO.write(writefiles, ds.header) do io
        for (i,p) in enumerate(ds) #I need to interate only to dataset, since it is Lazpoints
            allkeys = keys(pulse)
            if i in allkeys
                if pulse[i] == "unclassified"
                    p.classification = UInt(31)
                elseif pulse[i] == "firstpoint"
                    p.classification = UInt(3)
                elseif pulse[i] == "nextpoint"
                    p.classification = UInt(4)
                elseif pulse[i] == "lastpoint"
                    p.classification = UInt(5)
                end
                LazIO.writepoint(io,p)
            end
        end
    end
end

##MAIN##
#working with 4 different datasets
workdir = dirname(@__FILE__)
# dataname_1 = "NL3_clipped_water_cropheight_pointformat3_classification_sorted"
# dataname_2 = "3_splitted_sorted_pointformat3"
# dataname_3 = "12_splitted_sorted_pointformat3_clipZ"
dataname_4 = "14_splitted_sorted_pointformat3_clipZ"

const filenn_NL1 = joinpath(workdir,dataname_4 * ".laz")
const filenn_NL1_las = File{format"LAZ_"}(filenn_NL1)

#read specific LAS file
header, points = LazIO.load(filenn_NL1_las) #point cloud as LAS format
n = length(points) #length

#output files
outputname = "_classified"
writefiles = joinpath(workdir, dataname_4 * outputname * ".laz")

#read dataset
ds = LazIO.open(filenn_NL1)#dataset

#pulse
classified = classify(points)

#write output
write_output(writefiles,ds,classified)

ext: jl


text:
export meshcat,
       meshcat_material,
       meshcat_glass_material,
       meshcat_metal_material,
       meshcat_line_material

#=
ATTENTION!!!

Use Chrome debugger to visualize threeJS examples, break the execution and
evaluate in the console:

JSON.stringify(scene.toJSON(), undefined, 2)

=#

## Primitives
using UUIDs, MsgPack

send_meshcat(vis::Visualizer, obj) =
  write(vis.core, pack(obj))

send_setobject(vis, path, obj) =
  let msg = (type="set_object", path=path, object=obj),
      data = pack(msg)
    #vis.core.tree[path].object = data
    write(vis.core, data)
    path
  end

send_settransform(vis, path, transform) =
  let msg = (type="set_transform", path=path, matrix=convert(Vector{Float32}, transform)),
      data = pack(msg)
    #vis.core.tree[path].transform = data
    write(vis.core, data)
    nothing
  end

send_setproperty(vis, path, property, value) =
  let msg = (type="set_property", path=path, property=property, value=value),
      data = pack(msg)
    #vis.core.tree[path].properties[property] = data
    write(vis.core, data)
    nothing
  end

send_setvisible(vis, path, value) =
  send_meshcat(vis, (type="set_property", path=path, property="visible", value=value))

#send_setproperty(connection(meshcat), "/Background/hide_background", "value", true)
#send_meshcat(connection(meshcat), (type="hide_background",))

send_delobject(vis, path) =
  send_meshcat(vis, (type="delete", path=path))

meshcat_point(p::Loc) =
  let p = in_world(p)
    [cx(p),cz(p),cy(p)]
  end

meshcat_transform(p::Loc) =
  # MeshCat swaps Y with Z
  let m = translated_cs(p.cs, p.x, p.y, p.z).transform*[1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]
    m[:]
  end

meshcat_metadata(type="Object") =
  (type=type, version=4.5)

meshcat_material(color) =
  (uuid=string(uuid1()),
   type="MeshLambertMaterial",
    #"vertexColors" => 0,
    #"transparent" => false,
    #"opacity" => 1.0,
    #"depthTest" => true,
    #"linewidth" => 1.0,
    #"depthFunc" => 3,
    side=2,
    color="0x$(hex(color))",
    #color="0xAAAAAA",
    #"reflectivity" => 0.5,
    #depthWrite=true
    )
meshcat_glass_material(opacity=0.5, color=RGB(0.9,0.9,1.0)) =
  (uuid=string(uuid1()),
   type="MeshPhysicalMaterial",
    #"vertexColors" => 0,
   transparent=true,
   opacity=opacity,
    #"depthTest" => true,
    #"linewidth" => 1.0,
    #"depthFunc" => 3,
   side=2,
   color="0x$(hex(color))",
   reflectivity=0.1,
    #depthWrite=true
  )
meshcat_metal_material(roughness=0.5, color=RGB(1.0,1.0,1.0)) =
  (uuid=string(uuid1()),
   type="MeshStandardMaterial",
   metalness=1,
   roughness=roughness,
   side=2,
   color="0x$(hex(color))",
  )

meshcat_line_material(color) =
  (uuid=string(uuid1()),
   type="LineBasicMaterial",
   color="0x$(hex(color))",
   #linewidth=2, Due to limitations of the OpenGL Core Profile with the WebGL renderer on most platforms linewidth will always be 1 regardless of the set value.
   #depthFunc=3,
   #depthTest=true,
   depthWrite=true,
   #stencilWrite=false,
   #stencilWriteMask=255,
   #stencilFunc=519,
   #stencilRef=0,
   #stencilFuncMask=255,
   #stencilFail=7680,
   #stencilZFail=7680,
   #stencilZPass=7680
   )

meshcat_object(type, geom, material, p=u0(world_cs)) =
  (metadata=meshcat_metadata(),
   geometries=[geom],
   materials=[material],
   object=(uuid=string(uuid1()),
           type=type,
           geometry=geom.uuid,
           material=material.uuid,
           matrix=meshcat_transform(p)))

meshcat_object_2D(type, geom, shapes, material) =
 (metadata=meshcat_metadata(),
  shapes=shapes,
  geometries=[geom],
  materials=[material],
  object=(uuid=string(uuid1()),
          type=type,
          geometry=geom.uuid,
          material=material.uuid,
          matrix=(1, 0, 0, 0,
                  0, 1, 0, 0,
                  0, 0,-1, 0,
                  0, 0, 0, 1)))

meshcat_object_shapes(type, geom, shapes, material, p=u0(world_cs)) =
 (metadata=meshcat_metadata(),
  shapes=shapes,
  geometries=[geom],
  materials=[material],
  object=(uuid=string(uuid1()),
          type=type,
          geometry=geom.uuid,
          material=material.uuid,
          matrix=(translated_cs(p.cs, p.x, p.y, p.z).transform*[1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 1])[:])) # meshcat_transform(p))) #(1, 0, 0, 0, 0, 1, 0, 0, 0, 0,-1, 0, 0, 0, 0, 1)))

meshcat_buffer_geometry_attributes_position(vertices) =
  (itemSize=3,
   type="Float32Array",
   array=convert(Vector{Float32}, reduce(vcat, meshcat_point.(vertices))))

meshcat_line(vertices, material) =
  let geom = (uuid=string(uuid1()),
              type="BufferGeometry",
              data=(
                attributes=(
                  position=meshcat_buffer_geometry_attributes_position(vertices),),))
    meshcat_object("Line", geom, material)
  end

#=
Three.js uses 2D locations and 3D locations
=#

abstract type Meshcat2D end

convert(::Type{Meshcat2D}, p::Loc) = (cx(p),cy(p))
meshcat_2d(p::Loc) = let z =  @assert(abs(cz(p)) < 1e-10); (cx(p),cy(p)) end
meshcat_3d(p::Loc) = (cx(p),cy(p),cz(p))
meshcat_line_curve_2d(v1::Loc, v2::Loc) =
  (type="LineCurve", v1=meshcat_2d(v1), v2=meshcat_2d(v2))
meshcat_line_curve_3d(v1::Loc, v2::Loc) =
  (type="LineCurve3", v1=meshcat_3d(v1), v2=meshcat_3d(v2))

#=
Three.js provides a hierarchy of curves.
Curve - Abstract
2D curves:
  ArcCurve
  CubicBezierCurve
  EllipseCurve
  LineCurve
  QuadraticBezierCurve
  SplineCurve
3D curves:
  CatmullRomCurve3
  CubicBezierCurve3
  LineCurve3
  QuadraticBezierCurve3
Sequences:
  CurvePath - Abstract
    Path
      Shape
=#
abstract type MeshcatCurve end
abstract type MeshcatCurve2D <: MeshcatCurve end
abstract type MeshcatCurve3D <: MeshcatCurve end
abstract type MeshcatCurvePath <: MeshcatCurve end
abstract type MeshcatPath <: MeshcatCurvePath end
abstract type MeshcatShape <: MeshcatPath end

abstract type MeshcatCurves end

meshcat_curve(path) = convert(MeshcatCurve, path)
convert(::Type{MeshcatCurve}, p::CircularPath) =
  (type="EllipseCurve",
   aX=cx(p.center), aY=cy(p.center),
   xRadius=p.radius, yRadius=p.radius,
   aStartAngle=0, aEndAngle=2π,
   aClockwise=false,
   aRotation=0)
convert(::Type{MeshcatCurve}, p::OpenPolygonalPath) =
  let ps = path_vertices(p)
    length(ps) == 2 ?
      meshcat_line_curve_2d(ps[1], ps[2]) :
      convert(MeshcatPath, p)
  end

meshcat_path(path) = convert(MeshcatPath, path)
convert(::Type{MeshcatPath}, path) =
  (type="Path",
   curves=meshcat_curves(path),
   autoclose=false,
   currentPoint=(0,0))

meshcat_curves(path) = convert(MeshcatCurves, path)
convert(::Type{MeshcatCurves}, path) =
  [meshcat_curve(path)]
convert(::Type{MeshcatCurves}, vs::Locs) =
  [meshcat_line_curve_2d(v1, v2)
   for (v1,v2) in zip(vs, circshift(vs, -1))]
convert(::Type{MeshcatCurves}, p::Union{RectangularPath, ClosedPolygonalPath}) =
  convert(MeshcatCurves, path_vertices(p))

meshcat_shape(path) = convert(MeshcatShape, path)
convert(::Type{MeshcatShape}, p::PathSet) =
  (uuid=string(uuid1()),
   type="Shape",
   curves=meshcat_curves(p.paths[1]),
   autoclose=false,
   currentPoint=(0,0),
   holes=meshcat_path.(p.paths[2:end]))
convert(::Type{MeshcatShape}, p::Path) =
  (uuid=string(uuid1()),
   type="Shape",
   curves=meshcat_curves(p),
   autoclose=false,
   currentPoint=(0,0),
   holes=[])

meshcat_surface_2d(path, material, p=u0(world_cs)) =
  let shape = meshcat_shape(path),
      geom = (uuid=string(uuid1()),
              type="ShapeBufferGeometry",
              shapes=[shape.uuid],
              curveSegments=64)
    meshcat_object_shapes("Mesh", geom, [shape], material, p)
  end

#=
backend(meshcat)
delete_all_shapes()
add_object(
  meshcat,
  meshcat_surface_polygon_2d(
    path_vertices(closed_polygonal_path(regular_polygon_vertices(5, u0(), 3))),
    [path_vertices(closed_polygonal_path(regular_polygon_vertices(4, u0(), 1)))],
  material(meshcat)))

dump(meshcat_surface_polygon_2d(
  path_vertices(closed_polygonal_path(regular_polygon_vertices(5, u0(), 3))),
  [path_vertices(closed_polygonal_path(regular_polygon_vertices(4, u0(), 1)))],
  material(meshcat)))
dump(meshcat_surface_2d(
    path_set(closed_polygonal_path(regular_polygon_vertices(5, u0(), 3)),
             closed_polygonal_path(regular_polygon_vertices(4, u0(), 1))),
    material(meshcat)))

add_object(meshcat, meshcat_surface_polygon(regular_polygon_vertices(5), material(meshcat)))
add_object(
  meshcat,
  meshcat_surface_2d(
    path_set(
      circular_path(u0(), 2),
      circular_path(u0(), 1)),
    material(meshcat)))

add_object(
 meshcat,
  meshcat_surface_2d(circular_path(u0(), 3), material(meshcat)))

add_object(
  meshcat,
  meshcat_surface_2d(
    closed_polygonal_path(regular_polygon_vertices(5, u0(), 5)),
    material(meshcat)))
add_object(
  meshcat,
  meshcat_surface_2d(
    closed_polygonal_path(regular_polygon_vertices(4, u0(), 1)),
    material(meshcat)))

add_object(
  meshcat,
  meshcat_surface_2d(
    path_set(
      closed_polygonal_path(regular_polygon_vertices(5, u0(), 3)),
      closed_polygonal_path(regular_polygon_vertices(4, u0(), 1))),
    material(meshcat)))

dump(meshcat_surface_polygon(regular_polygon_vertices(5), material(meshcat)))

dump(meshcat_surface_2d(
  path_set(
    circular_path(u0(), 5),
    circular_path(u0(), 1)),
  material(meshcat)))
=#

meshcat_surface_polygon(vertices, material) =
  let n = length(vertices)
    n <= 0 ?
      meshcat_mesh(vertices, n < 4 ? [(0,1,2)] : [(0,1,2),(2,3,0)], material) :
      let ps = in_world.(vertices),
          n = vertices_normal(ps),
          cs = cs_from_o_vz(ps[1], n),
          vs = [in_cs(v, cs) for v in vertices]
        meshcat_surface_2d(closed_polygonal_path(vs), material, u0(cs))
      end
  end

#=
reset_backend()
delete_all_shapes()
add_object(meshcat, meshcat_surface_polygon([x(1), x(3), xz(3,4), xz(1,4)], material(meshcat)))
=#

meshcat_circle_mesh(center, radius, start_angle, amplitude, material) =
  let geom = (uuid=string(uuid1()),
              type="CircleBufferGeometry",
              radius=radius,
              segments=64,
              thetaStart=start_angle,
              thetaLength=amplitude),
      cs = cs_from_o_vz(center, vy(1, center.cs))
    meshcat_object("Mesh", geom, material, u0(cs))
  end

meshcat_centered_box(p, dx, dy, dz, material) =
  let geom = (uuid=string(uuid1()),
              type="BoxBufferGeometry",
              width=dx,
              depth=dy,
              height=dz)
    meshcat_object("Mesh", geom, material, p)
  end

meshcat_box(p, dx, dy, dz, material) =
  meshcat_centered_box(p+vxyz(dx/2, dy/2, dz/2, p.cs), dx, dy, dz, material)

#=
send_delobject(v, "/meshcat")
send_setobject(v, "/meshcat/box1", meshcat_box(xyz(1,1,1), 1, 2, 3))
send_setobject(v, "/meshcat/box2", meshcat_box(loc_from_o_vx_vy(u0(), vxy(1,1), vxy(-1,1)), 1, 2, 3))
send_delobject(v, "/meshcat/box2")
=#

meshcat_sphere(p, r, material) =
  let geom = (uuid=string(uuid1()),
              type="SphereBufferGeometry",
              radius=r,
              widthSegments=64,
              heightSegments=64)
    meshcat_object("Mesh", geom, material, p)
  end

meshcat_torus(p, re, ri, material) =
  let geom = (uuid=string(uuid1()),
              type="TorusBufferGeometry",
              radius=re,
              tube=ri,
              radialSegments=64,
              tubularSegments=16)
    meshcat_object_shapes("Mesh", geom, [], material, p)
  end

meshcat_centered_cone(p, r, h, material) =
  let geom = (uuid=string(uuid1()),
              type="ConeBufferGeometry",
              radius=r,
              height=h,
              radialSegments=64)
    meshcat_object("Mesh", geom, material, p)
  end

meshcat_cone(p, r, h, material) =
  meshcat_centered_cone(p+vz(h/2, p.cs), r, h, material)

meshcat_centered_cone_frustum(p, rb, rt, h, material) =
  let geom = (uuid=string(uuid1()),
              type="CylinderBufferGeometry",
              radiusTop=rt,
              radiusBottom=rb,
              height=h,
              radialSegments=64)
    meshcat_object("Mesh", geom, material, p)
  end

meshcat_cone_frustum(p, rb, rt, h, material) =
  meshcat_centered_cone_frustum(p+vz(h/2, p.cs), rb, rt, h, material)

meshcat_cylinder(p, r, h, material) =
  meshcat_centered_cone_frustum(p+vz(h/2, p.cs), r, r, h, material)

meshcat_extrusion_z(profile, h, material, p=u0(world_cs)) =
  let shape = meshcat_shape(profile),
      geom = (uuid=string(uuid1()),
              type="ExtrudeBufferGeometry",
              shapes=[shape.uuid],
              options=(
                #steps=2,
                depth=-h,
                bevelEnabled=false,
                #bevelThickness=1,
                #bevelSize=1,
                #bevelOffset=0,
                #bevelSegments=1,
                #extrudePath=meshcat_path(path),
                curveSegments=64
                ))
    meshcat_object_shapes("Mesh", geom, [shape], material, p)
  end
#

#=
delete_all_shapes()
add_object(
  meshcat,
  meshcat_extrusion_z(
    path_set(
      circular_path(u0(), 2),
      circular_path(u0(), 1)),
    10,
    material(meshcat),
    y(10)))

delete_all_shapes()
add_object(
  meshcat,
  meshcat_extrusion(open_polygonal_path([x(0), z(10)]),
                    closed_polygonal_path(regular_polygon_vertices(5, u0(), 3)),
                    material(meshcat)))
dump(meshcat_extrusion(open_polygonal_path([x(0), z(10)]),
                  closed_polygonal_path(regular_polygon_vertices(5, u0(), 3)),
                  material(meshcat)))
=#

#=
send_setobject(v, "/meshcat/sphere1", meshcat_sphere(xyz(-2,-3,0), 1))
send_setobject(v, "/meshcat/cylinder1", meshcat_cylinder(loc_from_o_vx_vy(x(-3), vxy(1,1), vxz(-1,1)), 1, 5))
=#
meshcat_mesh(vertices, faces, material) =
  let geom = (uuid=string(uuid1()),
              type="BufferGeometry",
              data=(
                attributes=(
                  position=meshcat_buffer_geometry_attributes_position(vertices),
                  #uv=?
                  ),
                index=(
                  itemSize=3,
                  type="Uint32Array",
                  array=convert(Vector{UInt32}, collect(Iterators.flatten(faces))))))
    meshcat_object("Mesh", geom, material, u0(world_cs))
  end
#=
meshcat_extruded_surface_polygon(outer_vertices, holes_vertices, material) =
  let ps = in_world.(outer_vertices),
      n = vertices_normal(ps),
      cs = cs_from_o_vz(ps[1], n),
      vs = [in_cs(v, cs) for v in outer_vertices],
      hsvs = [[in_cs(v, cs) for v in hvs] for hvs in holes_vertices],
      curves = [(#metadata=meshcat_metadata("Curve"),
                 type="LineCurve",
                 v1=(cx(v1),cy(v1)), v2=(cx(v2),cy(v2)))
                for (v1,v2) in zip(vs, circshift(vs, -1))],
      shape = (uuid=string(uuid1()),
               type="Shape",
               curves=curves,
               autoclose=false,
               currentPoint=(cx(vs[1]), cy(vs[1])),
               holes=[],
               ),
      geom = (uuid=string(uuid1()),
              type="ShapeBufferGeometry",
              shapes=[shape.uuid],
              #curveSegments=12,
              )
    #meshcat_object_2D("Mesh", geom, [shape], material, )
    meshcat_object_shapes("Mesh", geom, [shape], material, u0(cs))
  end
=#

letter_glyph = Dict{Char, NamedTuple}(
  ' '=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[]),
  '!'=>(bb=[xy(0, 0), xy(0, 1)], vss=[[xy(0, 0), xy(0, 1/6)], [xy(0, 1/3), xy(0, 1)]]),
  '"'=>(bb=[xy(0, 2/3), xy(1/3, 1)], vss=[[xy(0, 2/3), xy(1/6, 1)], [xy(1/3, 1), xy(1/6, 2/3)]]),
  '#'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1/3), xy(2/3, 1/3)], [xy(2/3, 2/3), xy(0, 2/3)], [xy(1/6, 1), xy(1/6, 0)], [xy(1/2, 0), xy(1/2, 1)]]),
  '$'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1/6), xy(1/2, 1/6), xy(2/3, 1/3), xy(1/2, 1/2), xy(1/6, 1/2), xy(0, 2/3), xy(1/6, 5/6), xy(2/3, 5/6)], [xy(1/3, 1), xy(1/3, 0)]]),
  '%'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(0, 5/6), xy(1/6, 5/6), xy(1/6, 1), xy(0, 1)], [xy(2/3, 1), xy(0, 0)], [xy(2/3, 0), xy(1/2, 0), xy(1/2, 1/6), xy(2/3, 1/6), xy(2/3, 0)]]),
  '&'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(2/3, 1/3), xy(1/3, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 1/3), xy(1/3, 2/3), xy(1/3, 5/6), xy(1/6, 1), xy(0, 5/6), xy(0, 2/3), xy(2/3, 0)]]),
  '\''=>(bb=[xy(0, 2/3), xy(1/6, 1)], vss=[[xy(0, 2/3), xy(1/6, 1)]]),
  '('=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(1/3, 1), xy(0, 2/3), xy(0, 1/3), xy(1/3, 0)]]),
  ')'=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(0, 1), xy(1/3, 2/3), xy(1/3, 1/3), xy(0, 0)]]),
  '*'=>(bb=[xy(0, 1/6), xy(2/3, 5/6)], vss=[[xy(1/3, 1/6), xy(1/3, 5/6)], [xy(2/3, 1/2), xy(0, 1/2)], [xy(2/3, 5/6), xy(0, 1/6)], [xy(0, 5/6), xy(2/3, 1/6)]]),
  '+'=>(bb=[xy(0, 1/6), xy(2/3, 5/6)], vss=[[xy(1/3, 1/6), xy(1/3, 5/6)], [xy(2/3, 1/2), xy(0, 1/2)]]),
  ','=>(bb=[xy(0, -1/6), xy(1/6, 1/6)], vss=[[xy(1/6, 1/6), xy(1/6, 0), xy(0, -1/6)]]),
  '-'=>(bb=[xy(0, 1/2), xy(2/3, 1/2)], vss=[[xy(0, 1/2), xy(2/3, 1/2)]]),
  '.'=>(bb=[xy(0, 0), xy(0, 1/6)], vss=[[xy(0, 0), xy(0, 1/6)]]),
  '/'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(2/3, 1)]]),
  '0'=>(bb=[xy(0, 0), xy(1/2, 1)], vss=[[xy(1/6, 0), xy(0, 1/6), xy(0, 5/6), xy(1/6, 1), xy(1/3, 1), xy(1/2, 5/6), xy(1/2, 1/6), xy(1/3, 0), xy(1/6, 0)]]),
  '1'=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(0, 5/6), xy(1/6, 1), xy(1/6, 0)], [xy(0, 0), xy(1/3, 0)]]),
  '2'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 5/6), xy(1/6, 1), xy(1/2, 1), xy(2/3, 5/6), xy(2/3, 2/3), xy(1/2, 1/2), xy(1/6, 1/2), xy(0, 1/3), xy(0, 0), xy(2/3, 0)]]),
  '3'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 5/6), xy(1/6, 1), xy(1/2, 1), xy(2/3, 5/6), xy(2/3, 2/3), xy(1/2, 1/2), xy(1/3, 1/2)], [xy(1/2, 1/2), xy(2/3, 1/3), xy(2/3, 1/6), xy(1/2, 0), xy(1/6, 0), xy(0, 1/6)]]),
  '4'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(2/3, 1/3), xy(0, 1/3), xy(1/2, 1), xy(1/2, 0)]]),
  '5'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1/6), xy(1/6, 0), xy(1/2, 0), xy(2/3, 1/6), xy(2/3, 1/2), xy(1/2, 2/3), xy(0, 2/3), xy(0, 1), xy(2/3, 1)]]),
  '6'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1/2), xy(1/2, 1/2), xy(2/3, 1/3), xy(2/3, 1/6), xy(1/2, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 2/3), xy(1/3, 1), xy(1/2, 1)]]),
  '7'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(2/3, 1), xy(1/6, 0)]]),
  '8'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(1/6, 0), xy(0, 1/6), xy(0, 1/3), xy(1/6, 1/2), xy(1/2, 1/2), xy(2/3, 2/3), xy(2/3, 5/6), xy(1/2, 1), xy(1/6, 1), xy(0, 5/6), xy(0, 2/3), xy(1/6, 1/2)], [xy(1/2, 1/2), xy(2/3, 1/3), xy(2/3, 1/6), xy(1/2, 0), xy(1/6, 0)]]),
  '9'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(1/6, 0), xy(1/3, 0), xy(2/3, 1/3), xy(2/3, 5/6), xy(1/2, 1), xy(1/6, 1), xy(0, 5/6), xy(0, 2/3), xy(1/6, 1/2), xy(2/3, 1/2)]]),
  ':'=>(bb=[xy(0, 1/6), xy(0, 2/3)], vss=[[xy(0, 2/3), xy(0, 1/2)], [xy(0, 1/3), xy(0, 1/6)]]),
  ';'=>(bb=[xy(0, -1/6), xy(1/6, 2/3)], vss=[[xy(1/6, 2/3), xy(1/6, 1/2)], [xy(1/6, 1/3), xy(1/6, 0), xy(0, -1/6)]]),
  '<'=>(bb=[xy(0, 0), xy(1/2, 1)], vss=[[xy(1/2, 1), xy(0, 1/2), xy(1/2, 0)]]),
  '='=>(bb=[xy(0, 1/3), xy(2/3, 2/3)], vss=[[xy(0, 2/3), xy(2/3, 2/3)], [xy(2/3, 1/3), xy(0, 1/3)]]),
  '>'=>(bb=[xy(0, 0), xy(1/2, 1)], vss=[[xy(0, 1), xy(1/2, 1/2), xy(0, 0)]]),
  '?'=>(bb=[xy(0, 0), xy(1/2, 1)], vss=[[xy(0, 5/6), xy(1/6, 1), xy(1/3, 1), xy(1/2, 5/6), xy(1/2, 2/3), xy(1/3, 1/2), xy(1/3, 1/3)], [xy(1/3, 1/6), xy(1/3, 0)]]),
  '@'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(1/2, 1/2), xy(1/3, 1/3), xy(1/6, 1/3), xy(1/6, 1/2), xy(1/3, 2/3), xy(1/2, 2/3), xy(1/2, 1/3), xy(2/3, 1/2), xy(2/3, 5/6), xy(1/2, 1), xy(1/6, 1), xy(0, 5/6), xy(0, 1/6), xy(1/6, 0), xy(2/3, 0)]]),
  'A'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1/3), xy(1/3, 1), xy(2/3, 1/3), xy(2/3, 0)], [xy(0, 1/3), xy(2/3, 1/3)]]),
  'B'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(1/2, 0), xy(2/3, 1/6), xy(2/3, 1/3), xy(1/2, 1/2), xy(1/6, 1/2)], [xy(1/2, 1/2), xy(2/3, 2/3), xy(2/3, 5/6), xy(1/2, 1), xy(0, 1)], [xy(1/6, 1), xy(1/6, 0)]]),
  'C'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(2/3, 1/6), xy(1/2, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 5/6), xy(1/6, 1), xy(1/2, 1), xy(2/3, 5/6)]]),
  'D'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(1/2, 0), xy(2/3, 1/6), xy(2/3, 5/6), xy(1/2, 1), xy(0, 1)], [xy(1/6, 1), xy(1/6, 0)]]),
  'E'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(2/3, 1)], [xy(0, 1/2), xy(1/3, 1/2)], [xy(0, 0), xy(2/3, 0)]]),
  'F'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(2/3, 1)], [xy(0, 1/2), xy(1/3, 1/2)]]),
  'G'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(1/2, 1/2), xy(2/3, 1/2), xy(2/3, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 5/6), xy(1/6, 1), xy(2/3, 1)]]),
  'H'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1)], [xy(0, 1/2), xy(2/3, 1/2)], [xy(2/3, 1), xy(2/3, 0)]]),
  'I'=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(0, 1), xy(1/3, 1)], [xy(1/6, 1), xy(1/6, 0)], [xy(0, 0), xy(1/3, 0)]]),
  'J'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1/6), xy(1/6, 0), xy(1/2, 0), xy(2/3, 1/6), xy(2/3, 1)]]),
  'K'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1)], [xy(2/3, 1), xy(1/6, 1/2), xy(0, 1/2)], [xy(1/6, 1/2), xy(2/3, 0)]]),
  'L'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(0, 0), xy(2/3, 0)]]),
  'M'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(1/3, 1/3), xy(2/3, 1), xy(2/3, 0)]]),
  'N'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(2/3, 0), xy(2/3, 1)]]),
  'O'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(2/3, 1), xy(2/3, 0), xy(0, 0)]]),
  'P'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(1/2, 1), xy(2/3, 5/6), xy(2/3, 2/3), xy(1/2, 1/2), xy(0, 1/2)]]),
  'Q'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(1/3, 1/3), xy(1/2, 1/6), xy(1/3, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 5/6), xy(1/6, 1), xy(1/2, 1), xy(2/3, 5/6), xy(2/3, 1/3), xy(1/2, 1/6), xy(2/3, 0)]]),
  'R'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(1/2, 1), xy(2/3, 5/6), xy(2/3, 2/3), xy(1/2, 1/2), xy(0, 1/2)], [xy(1/6, 1/2), xy(2/3, 0)]]),
  'S'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1/6), xy(1/6, 0), xy(1/2, 0), xy(2/3, 1/6), xy(0, 5/6), xy(1/6, 1), xy(1/2, 1), xy(2/3, 5/6)]]),
  'T'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(2/3, 1)], [xy(1/3, 1), xy(1/3, 0)]]),
  'U'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(0, 1/6), xy(1/6, 0), xy(1/2, 0), xy(2/3, 1/6), xy(2/3, 1)]]),
  'V'=>(bb=[xy(0, 0), xy(1, 1)], vss=[[xy(0, 1), xy(1/2, 0), xy(1, 1)]]),
  'W'=>(bb=[xy(0, 0), xy(1, 1)], vss=[[xy(0, 1), xy(1/3, 0), xy(1/2, 1/2), xy(2/3, 0), xy(1, 1)]]),
  'X'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(2/3, 1)], [xy(0, 1), xy(2/3, 0)]]),
  'Y'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(1/3, 1/2), xy(1/3, 0)], [xy(1/3, 1/2), xy(2/3, 1)]]),
  'Z'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(2/3, 1), xy(0, 0), xy(2/3, 0)]]),
  '['=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(0, 0), xy(0, 1), xy(1/3, 1)], [xy(1/3, 0), xy(0, 0)]]),
  '\\'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1), xy(2/3, 0)]]),
  ']'=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(0, 1), xy(1/3, 1), xy(1/3, 0), xy(0, 0)]]),
  '^'=>(bb=[xy(0, 2/3), xy(2/3, 1)], vss=[[xy(0, 2/3), xy(1/3, 1), xy(2/3, 2/3)]]),
  '_'=>(bb=[xy(0, -1/6), xy(2/3, -1/6)], vss=[[xy(0, -1/6), xy(2/3, -1/6)]]),
  '`'=>(bb=[xy(0, 2/3), xy(1/6, 1)], vss=[[xy(0, 1), xy(1/6, 2/3)]]),
  'a'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(1/3, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 1/2), xy(1/6, 2/3), xy(1/3, 2/3), xy(1/2, 1/2), xy(1/2, 1/6), xy(1/3, 0)], [xy(1/2, 1/6), xy(2/3, 0)]]),
  'b'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1)], [xy(0, 1/3), xy(1/3, 2/3), xy(1/2, 2/3), xy(2/3, 1/2), xy(2/3, 1/6), xy(1/2, 0), xy(1/3, 0), xy(0, 1/3)]]),
  'c'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(2/3, 2/3), xy(1/6, 2/3), xy(0, 1/2), xy(0, 1/6), xy(1/6, 0), xy(2/3, 0)]]),
  'd'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(2/3, 1/3), xy(1/3, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 1/2), xy(1/6, 2/3), xy(1/3, 2/3), xy(2/3, 1/3)], [xy(2/3, 1), xy(2/3, 0)]]),
  'e'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 1/3), xy(1/2, 1/3), xy(2/3, 1/2), xy(1/2, 2/3), xy(1/6, 2/3), xy(0, 1/2), xy(0, 1/6), xy(1/6, 0), xy(1/2, 0)]]),
  'f'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 1/2), xy(1/2, 1/2)], [xy(2/3, 5/6), xy(1/2, 1), xy(1/3, 1), xy(1/6, 5/6), xy(1/6, 0)]]),
  'g'=>(bb=[xy(0, -1/3), xy(2/3, 2/3)], vss=[[xy(0, -1/6), xy(1/6, -1/3), xy(1/2, -1/3), xy(2/3, -1/6), xy(2/3, 1/2), xy(1/2, 2/3), xy(1/6, 2/3), xy(0, 1/2), xy(0, 1/6), xy(1/6, 0), xy(2/3, 0)]]),
  'h'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1)], [xy(0, 1/3), xy(1/3, 2/3), xy(1/2, 2/3), xy(2/3, 1/2), xy(2/3, 0)]]),
  'i'=>(bb=[xy(0, 0), xy(0, 1)], vss=[[xy(0, 0), xy(0, 2/3)], [xy(0, 5/6), xy(0, 1)]]),
  'j'=>(bb=[xy(0, -1/3), xy(1/2, 1)], vss=[[xy(0, -1/6), xy(1/6, -1/3), xy(1/3, -1/3), xy(1/2, -1/6), xy(1/2, 2/3)], [xy(1/2, 5/6), xy(1/2, 1)]]),
  'k'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 0), xy(0, 1)], [xy(0, 1/3), xy(1/3, 1/3), xy(2/3, 2/3)], [xy(1/3, 1/3), xy(2/3, 0)]]),
  'l'=>(bb=[xy(0, 0), xy(1/6, 1)], vss=[[xy(0, 1), xy(0, 1/6), xy(1/6, 0)]]),
  'm'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 0), xy(0, 2/3)], [xy(0, 1/2), xy(1/6, 2/3), xy(1/3, 1/2), xy(1/3, 1/3)], [xy(1/3, 1/2), xy(1/2, 2/3), xy(2/3, 1/2), xy(2/3, 0)]]),
  'n'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 0), xy(0, 2/3)], [xy(0, 1/3), xy(1/3, 2/3), xy(1/2, 2/3), xy(2/3, 1/2), xy(2/3, 0)]]),
  'o'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(1/2, 0), xy(1/6, 0), xy(0, 1/6), xy(0, 1/2), xy(1/6, 2/3), xy(1/2, 2/3), xy(2/3, 1/2), xy(2/3, 1/6), xy(1/2, 0)]]),
  'p'=>(bb=[xy(0, -1/3), xy(2/3, 2/3)], vss=[[xy(0, -1/3), xy(0, 2/3)], [xy(0, 1/2), xy(1/6, 2/3), xy(1/2, 2/3), xy(2/3, 1/2), xy(2/3, 1/6), xy(1/2, 0), xy(0, 0)]]),
  'q'=>(bb=[xy(0, -1/3), xy(2/3, 2/3)], vss=[[xy(2/3, -1/3), xy(2/3, 2/3)], [xy(2/3, 1/2), xy(1/2, 2/3), xy(1/6, 2/3), xy(0, 1/2), xy(0, 1/6), xy(1/6, 0), xy(2/3, 0)]]),
  'r'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 0), xy(0, 2/3)], [xy(0, 1/3), xy(1/3, 2/3), xy(1/2, 2/3), xy(2/3, 1/2)]]),
  's'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 0), xy(1/2, 0), xy(2/3, 1/6), xy(1/2, 1/3), xy(1/6, 1/3), xy(0, 1/2), xy(1/6, 2/3), xy(2/3, 2/3)]]),
  't'=>(bb=[xy(0, 0), xy(2/3, 1)], vss=[[xy(0, 2/3), xy(2/3, 2/3)], [xy(1/3, 1), xy(1/3, 1/6), xy(1/2, 0), xy(2/3, 1/6)]]),
  'u'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 2/3), xy(0, 1/6), xy(1/6, 0), xy(1/3, 0), xy(2/3, 1/3)], [xy(2/3, 2/3), xy(2/3, 0)]]),
  'v'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 2/3), xy(1/3, 0), xy(2/3, 2/3)]]),
  'w'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 2/3), xy(1/6, 0), xy(1/3, 2/3), xy(1/2, 0), xy(2/3, 2/3)]]),
  'x'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 0), xy(2/3, 2/3)], [xy(0, 2/3), xy(2/3, 0)]]),
  'y'=>(bb=[xy(0, -1/3), xy(2/3, 2/3)], vss=[[xy(0, 2/3), xy(1/3, 0)], [xy(2/3, 2/3), xy(1/6, -1/3), xy(0, -1/3)]]),
  'z'=>(bb=[xy(0, 0), xy(2/3, 2/3)], vss=[[xy(0, 2/3), xy(2/3, 2/3), xy(0, 0), xy(2/3, 0)]]),
  '{'=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(1/3, 1), xy(1/6, 5/6), xy(1/6, 2/3), xy(0, 1/2), xy(1/6, 1/3), xy(1/6, 1/6), xy(1/3, 0)]]),
  '|'=>(bb=[xy(0, 0), xy(0, 1)], vss=[[xy(0, 0), xy(0, 1)]]),
  '}'=>(bb=[xy(0, 0), xy(1/3, 1)], vss=[[xy(0, 0), xy(1/6, 1/6), xy(1/6, 1/3), xy(1/3, 1/2), xy(1/6, 2/3), xy(1/6, 5/6), xy(0, 1)]]),
  '~'=>(bb=[xy(0, 1/2), xy(2/3, 2/3)], vss=[[xy(0, 1/2), xy(1/6, 2/3), xy(1/2, 1/2), xy(2/3, 2/3)]])
)

inter_letter_spacing_factor = 1/3

meshcat_text(str::String, p::Loc, size::Real) =
  for c in str
    let glyph = letter_glyph[c]
      for vs in glyph.vss
        line([add_xy(p, q.x*size, q.y*size) for q in vs])
      end
      p = add_x(p, (glyph.bb[2].x-glyph.bb[1].x + inter_letter_spacing_factor)*size)
    end
  end

send_setview(v, camera::Loc, target::Loc, lens::Real, aperture::Real) =
  let (x1,y1,z1) = raw_point(camera),
      (x2,y2,z2) = raw_point(target)
    send_settransform(v, "/Cameras/default", [
      1, 0, 0, 0,
      0, 1, 0, 0,
      0, 0, 1, 0,
      x2, y2, z2, 1])
    send_setproperty(v, "/Cameras/default/rotated/<object>", "zoom", 1)
    send_setproperty(v, "/Cameras/default/rotated/<object>", "fov", view_angles(lens)[2])
    send_setproperty(v, "/Cameras/default/rotated/<object>", "position", [x1-x2,z1-z2,y2-y1])
  end

#send_setproperty(connection(meshcat), "/Orbit/<object>", "target", [2,0,0])
####################################################
abstract type MCATKey end
const MCATId = String #MCATlyJS.GenericTrace
const MCATRef = GenericRef{MCATKey, MCATId}
const MCATNativeRef = NativeRef{MCATKey, MCATId}
const MCATUnionRef = UnionRef{MCATKey, MCATId}
const MCATSubtractionRef = SubtractionRef{MCATKey, MCATId}

struct MCATLayer
  name
  material
  line_material
end

mcat_layer(name, color) =
  MCATLayer(name, meshcat_material(color), meshcat_line_material(color))

mutable struct MCATBackend{K,T} <: Backend{K,T}
  connection::LazyParameter{Visualizer}
  count::Int64
  layer::MCATLayer
  camera::Loc
  target::Loc
  lens::Real
  sun_altitude::Real
  sun_azimuth::Real

end

const MCAT = MCATBackend{MCATKey, MCATId}

connection(b::MCAT) = b.connection()
material(b::MCAT) = b.layer.material
line_material(b::MCAT) = b.layer.line_material

create_MCAT_connection() =
  let (width, height) = render_size(),
      vis = Visualizer()
    display(MCATViewer(vis))
    send_setproperty(vis, "/Cameras/default/rotated/<object>", "far", 50000)
    vis
  end

meshcat = MCAT(LazyParameter(Visualizer, create_MCAT_connection),
               0,
               mcat_layer("default", RGB(1,1,1)),
               xyz(10,10,10),
               xyz(0,0,0),
               35,
               90,
               0)

#=
To visualize, we piggyback on Julia's display mechanisms
=#

display_meshcat(io, vis, (w, h) = render_size()) =
  let frame = vis.core
    print(io, """
    <div style="height: $(h)px; width: $(w)px; overflow-x: auto; overflow-y: hidden; resize: both">
    <iframe src="$(MeshCat.url(frame))" style="width: 100%; height: 100%; border: none"></iframe>
    </div>
    """)
    MeshCat.wait_for_server(frame)
  end

struct MCATViewer
  visualizer
end

Base.show(
  io::IO,
  ::Union{
    MIME"text/html",
    MIME"juliavscode/html",
    MIME"application/prs.juno.plotpane+html"},
  v::MCATViewer) = display_meshcat(io, v.visualizer)

export display_view
display_view(b::Backend=current_backend()) = error("Needs specialization")
display_view(b::MCAT) = MCATViewer(connection(b))

const meshcat_root_path = "/Khepri"

next_id(b::MCATBackend{K,T}) where {K,T} =
  begin
      b.count += 1
      string(meshcat_root_path, "/", b.layer.name, "/", b.count - 1)
  end

add_object(b::MCAT, obj) =
  send_setobject(connection(b), next_id(b), obj)

has_boolean_ops(::Type{MCAT}) = HasBooleanOps{false}()
void_ref(b::MCAT) = MCATNativeRef("")

reset_backend(b::MCAT) =
  begin
    display(render(connection(b)))
    set_view(get_view()...)
  end

new_backend(b::MCAT) =
  begin
    reset(b.connection)
    display_view(b)
  end

#

#=
setprop!(
  connection(meshcat)["/Cameras/default/rotated/<object>"], #"Cameras","default","rotated","<object>"],
  "zoom",
  1.0)
=#

#=
function MeshCat.lower(s::Sphere)
  error("BUM")
    Dict{String, Any}(
        "uuid" => string(uuid1()),
        "type" => "SphereGeometry",
        "radius" => s.radius,
        "widthSegments" => 20,
        "heightSegments" => 20,
    )
end

add_object!(meshcat, s)
=#
#=
MeshCat.setcontrol!(
  connection(meshcat)["Background"], #"Cameras","default","rotated","<object>"],
  "hide_background")

{
    type: "set_property",
    path: "/Cameras/default/rotated/<object>",
    property: "zoom",
    value: 2.0
}
=#

set_view(camera::Loc, target::Loc, lens::Real, aperture::Real, b::MCAT) =
  let v = connection(b)
    send_setview(v, camera, target, lens, aperture)
    b.camera = camera
    b.target = target
    b.lens = lens
  end

get_view(b::MCAT) =
  b.camera, b.target, b.lens

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

render_view(path::String, b::MCAT) = false
  #=
  let v = connection(b)
      @remote(b, Render(
               render_width(), render_height(),
               path,
               convert_render_quality(b, render_quality()),
               convert_render_exposure(b, render_exposure())))
=#
###################################
#
delete_all_shapes(b::MCAT) =
  let c = connection(b)
    send_delobject(c, meshcat_root_path)
  end

backend_delete_shapes(b::MCAT, shapes::Shapes) =
  let c = connection(b)
    for r in collect_ref(shapes)
      send_delobject(c, r)
    end
  end

realize(b::MCAT, s::Text) =
  (meshcat_text(s.str, s.corner, s.height); MCATNativeRef(""))

backend_stroke(b::MCAT, path::CircularPath) =
  add_object(b, meshcat_line(path_frames(path), line_material(b)))

backend_stroke(b::MCAT, path::RectangularPath) =
  let c = path.corner,
      dx = path.dx,
      dy = path.dy
    add_object(
      b,
      meshcat_line([c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy), c],
                   line_material(b)))
  end
backend_polygon(b::MCAT, vs::Locs) =
  backend_stroke_line(b, [vs..., vs[1]])
backend_surface_polygon(b::MCAT, vs::Locs) =
  add_object(b, meshcat_surface_polygon(vs, material(b)))
backend_stroke(b::MCAT, path::ArcPath) =
  add_object(
    b,
    meshcat_line(
      path_frames(arc_path(path.center, path.radius, path.start_angle, path.amplitude)),
      line_material(b)))

backend_fill(b::MCAT, path::ClosedPolygonalPath) =
  backend_surface_polygon(b, path.vertices)

backend_fill(b::MCAT, path::RectangularPath) =
  let c = path.corner,
      dx = path.dx,
      dy = path.dy
    add_object(b, meshcat_surface_polygon([c, add_x(c, dx), add_xy(c, dx, dy), add_y(c, dy)], material(b)))
  end

backend_stroke(b::MCAT, path::OpenSplinePath) =
  if (path.v0 == false) && (path.v1 == false)
    add_object(b, meshcat_line(path_frames(path), line_material(b)))
  elseif (path.v0 != false) && (path.v1 != false)
    @remote(b, InterpSpline(path.vertices, path.v0, path.v1))
  else
    @remote(b, InterpSpline(
                     path.vertices,
                     path.v0 == false ? path.vertices[2]-path.vertices[1] : path.v0,
                     path.v1 == false ? path.vertices[end-1]-path.vertices[end] : path.v1))
  end
backend_stroke(b::MCAT, path::ClosedSplinePath) =
    add_object(b, meshcat_line(path_frames(path), line_material(b)))
backend_fill(b::MCAT, path::ClosedSplinePath) =
    add_object(b, meshcat_surface_polygon(path_frames(path), material(b)))

#=
backend_fill_curves(b::MCAT, refs::MCATIds) = @remote(b, SurfaceFromCurves(refs))
backend_fill_curves(b::MCAT, ref::MCATId) = @remote(b, SurfaceFromCurves([ref]))
=#
backend_stroke_line(b::MCAT, vs) =
	add_object(b, meshcat_line(vs, line_material(b)))
#=
realize(b::MCAT, s::Point) =
  @remote(b, Point(s.position))
  =#
realize(b::MCAT, s::Line) =
  backend_stroke_line(b, s.vertices)

realize(b::MCAT, s::Spline) =
 # This should be merged with opensplinepath
  if (s.v0 == false) && (s.v1 == false)
    backend_stroke(b, open_spline_path(s.points))
  elseif (s.v0 != false) && (s.v1 != false)
    error("Finish this")
  else
    error("Finish this")
  end

realize(b::MCAT, s::ClosedSpline) =
  backend_stroke(b, closed_spline_path(s.points))

realize(b::MCAT, s::Circle) =
  backend_stroke(b, circular_path(s.center, s.radius))
realize(b::MCAT, s::Arc) =
  backend_stroke(b, arc_path(s.center, s.radius, s.start_angle, s.amplitude))
#=
realize(b::MCAT, s::Ellipse) =
  if s.radius_x > s.radius_y
    @remote(b, Ellipse(s.center, vz(1, s.center.cs), vxyz(s.radius_x, 0, 0, s.center.cs), s.radius_y/s.radius_x))
  else
    @remote(b, Ellipse(s.center, vz(1, s.center.cs), vxyz(0, s.radius_y, 0, s.center.cs), s.radius_x/s.radius_y))
  end
realize(b::MCAT, s::EllipticArc) =
  error("Finish this")
=#
realize(b::MCAT, s::SurfaceCircle) =
  add_object(b, meshcat_circle_mesh(s.center, s.radius, 0, 2pi, material(b)))

realize(b::MCAT, s::Surface) =
  let #ids = map(r->@remote(b, NurbSurfaceFrom(r)), @remote(b, SurfaceFromCurves(collect_ref(s.frontier))))
      ids = @remote(b, SurfaceFromCurves(collect_ref(s.frontier)))
    foreach(mark_deleted, s.frontier)
    ids
  end
backend_surface_boundary(b::MCAT, s::Shape2D) =
    map(c -> shape_from_ref(c, b), @remote(b, CurvesFromSurface(ref(s).value)))

realize(b::MCAT, s::Sphere) =
  add_object(b, meshcat_sphere(s.center, s.radius, material(b)))

realize(b::MCAT, s::Torus) =
  add_object(b, meshcat_torus(s.center, s.re, s.ri, material(b)))

backend_right_cuboid(b::MCAT, cb, width, height, h, material) =
  add_object(b, meshcat_box(add_xy(cb, -width/2, -height/2), width, height, h, material))

realize(b::MCAT, s::Box) =
  add_object(b, meshcat_box(s.c, s.dx, s.dy, s.dz, material(b)))

realize(b::MCAT, s::Cone) =
  add_object(b, meshcat_cone(s.cb, s.r, s.h, material(b)))

realize(b::MCAT, s::ConeFrustum) =
  add_object(b, meshcat_cone_frustum(s.cb, s.rb, s.rt, s.h, material(b)))

backend_cylinder(b::MCAT, cb::Loc, r::Real, h::Real, material=material(b)) =
  send_setobject(connection(b), next_id(b), meshcat_cylinder(cb, r, h, material))

meshcat_faces(si, sj, closed_u, closed_v) =
  let idx(i,j) = (i-1)*sj+(j-1),
      idxs = [],
      quad(a,b,c,d) = (push!(idxs, (a, b, d)); push!(idxs, (d, b, c)))
    for i in 1:si-1
      for j in 1:sj-1
        quad(idx(i,j), idx(i+1,j), idx(i+1,j+1), idx(i,j+1))
      end
      if closed_v
        quad(idx(i,sj), idx(i+1,sj), idx(i+1,1), idx(i,1))
      end
    end
    if closed_u
      for j in 1:sj-1
        quad(idx(si,j), idx(1,j), idx(1,j+1), idx(si,j+1))
      end
      if closed_v
        quad(idx(si,sj), idx(1,sj), idx(1,1), idx(si,1))
      end
    end
    idxs
  end

realize(b::MCAT, s::SurfaceGrid) =
  backend_surface_grid(b, s.points, s.closed_u, s.closed_v, s.smooth_u, s.smooth_v)

backend_surface_grid(b::MCAT, pts, closed_u, closed_v, smooth_u, smooth_v) =
  let si = size(pts, 1),
      sj = size(pts, 2),
      idxs = meshcat_faces(si, sj, closed_u, closed_v)
    add_object(b, meshcat_mesh(reshape(permutedims(pts),:), idxs, material(b)))
  end

smooth_pts(pts) = in_world.(path_frames(open_spline_path(pts)))

backend_surface_mesh(b::MCAT, vertices, faces) =
  add_object(b, meshcat_mesh(vertices, faces, material(b)))

#=
realize(b::MCAT, s::SweepPath) =
  let vertices = in_world.(path_vertices(s.profile)),
      frames = map_division(identity, s.path, 20),
      buf = buffer(b),
      mat = get_material(b, s)
    write_MCAT_mesh(
      buf,
      mat,
      [xyz(cx(p), cy(p), cz(p), frame.cs) for p in vertices, frame in frames],
      is_closed_path(s.profile),
      is_closed_path(s.path),
      is_smooth_path(s.profile),
      is_smooth_path(s.path))
    void_ref(b)
  end

=#
# HACK: JUST FOR TESTING
realize(b::MCAT, s::Thicken) =
  realize(b, s.shape)

realize(b::MCAT, s::EmptyShape) = void_ref(b)

realize(b::MCAT, s::UniversalShape) = void_ref(b)

#=
#=
realize(b::MCAT, s::Move) =
    let r = map_ref(s.shape) do r
                MCATMove(connection(b), r, s.v)
                r
            end
        mark_deleted(s.shape)
        r
    end

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

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

=#

realize(b::MCAT, s::UnionShape) =
  let shapes = filter(! is_empty_shape, s.shapes)
    length(shapes) == 1 ?
      (ref(shapes[1]); delete_shape(shapes[1])) :
      write_MCAT_object(buffer(b), "union", get_material(b, s)) do
        for ss in shapes
          ref(ss)
          delete_shape(ss)
        end
      end
    void_ref(b)
  end

realize(b::MCAT, s::IntersectionShape) =
  write_MCAT_object(buffer(b), "intersection", get_material(b, s)) do
    for ss in s.shapes
      ref(ss)
      delete_shape(ss)
    end
    void_ref(b)
  end

realize(b::MCAT, s::SubtractionShape3D) =
  write_MCAT_object(buffer(b), "difference", get_material(b, s)) do
    ref(s.shape)
    delete_shape(s.shape)
    for ss in s.shapes
      ref(ss)
      delete_shape(ss)
    end
    void_ref(b)
  end
=#
# BIM

realize_box(b::MCAT, mat, p, dx, dy, dz) =
  add_object(b, meshcat_box(p, dx, dy, dz, mat))

realize_prism(b::MCAT, top, bot, side, path::PathSet, h::Real) =
  let p = path_start(path.paths[1]),
      n = planar_path_normal(path.paths[1]),
      cs = cs_from_o_vz(p, n)
    add_object(b, meshcat_extrusion_z(in_cs(path, cs), h, top, u0(cs)))
  end
#=
  # PathSets require a different approach
  let buf = buffer(b),
      bot_vss = map(path_vertices, path.paths),
      top_vss = map(path_vertices, translate(path, vz(h)).paths)
    write_MCAT_polygons(buf, bot, map(reverse, bot_vss))
    write_MCAT_polygons(buf, top, top_vss)
    for (bot_vs, top_vs) in zip(bot_vss, top_vss)
      for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
        write_MCAT_polygon(buf, side, vs)
      end
    end
    void_ref(b)
  end
=#

realize_pyramid_frustum(b::MCAT, top, bot, side, bot_vs::Locs, top_vs::Locs, closed=true) =
  let
    if closed
      add_object(b, meshcat_surface_polygon(reverse(bot_vs), bot))
      add_object(b, meshcat_surface_polygon(top_vs, top))
    end
    add_object(b,
      meshcat_mesh([top_vs...,bot_vs...],
                   meshcat_faces(2, length(bot_vs), false, false),
                   material(b)))
#    for vs in zip(bot_vs, circshift(bot_vs, 1), circshift(top_vs, 1), top_vs)
#      add_object(b, meshcat_surface_polygon(vs, side))
#    end
    void_ref(b)
  end

#meshcat_faces(2, 3, false, false)

#regular_pyramid_frustum(5)

backend_surface_polygon(b::MCAT, mat, path::PathSet, acw=true) =
  backend_surface_polygon(b, mat, path_vertices(path.paths[1])) # WARNING: Ignore the remaining paths that should be subtracted

backend_surface_polygon(b::MCAT, mat, vs::Locs, acw=true) =
  add_object(b, meshcat_surface_polygon(vs, mat))
#=
# Polygons with holes need a PathSets in MCAT

subtract_paths(b::MCAT, c_r_w_path::PathSet, c_l_w_path::PathSet, c_r_op_path, c_l_op_path) =
  path_set(c_r_w_path.paths..., c_r_op_path),
  path_set(c_l_w_path.paths..., c_l_op_path)

subtract_paths(b::MCAT, c_r_w_path::Path, c_l_w_path::Path, c_r_op_path, c_l_op_path) =
  path_set(c_r_w_path, c_r_op_path),
  path_set(c_l_w_path, c_l_op_path)
=#

# Layers
current_layer(b::MCAT) =
  b.layer

current_layer(layer, b::MCAT) =
  b.layer = layer

backend_create_layer(b::MCAT, name::String, active::Bool, color::RGB) =
  begin
    @assert active
    mcat_layer(name, color)
  end

#=
create_ground_plane(shapes, material=default_MCAT_ground_material()) =
  if shapes == []
    error("No shapes selected for analysis. Use add-MCAT-shape!.")
  else
    let (p0, p1) = bounding_box(union(shapes)),
        (center, ratio) = (quad_center(p0, p1, p2, p3),
                  distance(p0, p4)/distance(p0, p2));
     ratio == 0 ?
      error("Couldn"t compute height. Use add-MCAT-shape!.") :
      let pts = map(p -> intermediate_loc(center, p, ratio*10), [p0, p1, p2, p3]);
         create_surface_layer(pts, 0, ground_layer(), material)
        end
       end
  end
        w = max(floor_extra_factor()*distance(p0, p1), floor_extra_width())
        with(current_layer,floor_layer()) do
          box(xyz(min(p0.x, p1.x)-w, min(p0.y, p1.y)-w, p0.z-1-floor_distance()),
              xyz(max(p0.x, p1.x)+w, max(p0.y, p1.y)+w, p0.z-0-floor_distance()))
        end
      end
    end

=#

#=
#FIXME define the family parameters for beams
realize(b::MCAT, s::Beam) =
    ref(right_cuboid(s.p0, 0.2, 0.2, s.p1, 0))
=#
realize(b::MCAT, s::Union{Door, Window}) =
  void_ref(b)

####################################################
MeshcatMaterial = NamedTuple

meshcat_material_family(mat::MeshcatMaterial) =
  BackendMaterialFamily(mat)

meshcat_slab_family(top::MeshcatMaterial, bot::MeshcatMaterial=top, side::MeshcatMaterial=bot) =
  BackendSlabFamily(top, bot, side)

meshcat_roof_family(top::MeshcatMaterial, bot::MeshcatMaterial=top, side::MeshcatMaterial=bot) =
  BackendRoofFamily(top, bot, side)

meshcat_wall_family(right::MeshcatMaterial, left::MeshcatMaterial=right) =
  BackendWallFamily(right, left)

export meshcat_material_family,
       meshcat_slab_family,
       meshcat_roof_family,
       meshcat_wall_family,
       default_meshcat_material

meshcat_generic_floor_20 = meshcat_material(RGB(0.2,0.2,0.2))
meshcat_generic_ceiling_80 = meshcat_material(RGB(0.8,0.8,0.8))
meshcat_generic_interior_wall_70 = meshcat_material(RGB(0.7,0.7,0.7))
meshcat_outside_facade_30 = meshcat_material(RGB(0.3,0.3,0.3))
meshcat_generic_metal = meshcat_metal_material()
meshcat_generic_furniture_50 = meshcat_material(RGB(0.5,0.5,0.5))
meshcat_glass = meshcat_glass_material()

set_backend_family(default_wall_family(), meshcat, meshcat_wall_family(meshcat_generic_interior_wall_70))
#set_backend_family(default_slab_family(), meshcat, meshcat_slab_family(meshcat_generic_floor_20, meshcat_generic_ceiling_80))
set_backend_family(default_slab_family(), meshcat, meshcat_material_family(meshcat_generic_ceiling_80))
set_backend_family(default_roof_family(), meshcat, meshcat_roof_family(meshcat_generic_floor_20, meshcat_outside_facade_30))
set_backend_family(default_beam_family(), meshcat, meshcat_material_family(meshcat_generic_metal))
set_backend_family(default_column_family(), meshcat, meshcat_material_family(meshcat_generic_metal))
set_backend_family(default_door_family(), meshcat, meshcat_material_family(meshcat_generic_furniture_50))
set_backend_family(default_panel_family(), meshcat, meshcat_material_family(meshcat_glass))

set_backend_family(default_truss_node_family(), meshcat, meshcat_material_family(meshcat_generic_metal))
set_backend_family(default_truss_bar_family(), meshcat, meshcat_material_family(meshcat_generic_metal))

#set_backend_family(default_table_family(), meshcat, meshcat_resource_family("Default/Prefabs/Table"))
#set_backend_family(default_chair_family(), meshcat, meshcat_resource_family("Default/Prefabs/Chair"))
#set_backend_family(default_table_chair_family(), meshcat, meshcat_resource_family("Default/Prefabs/TableChair"))

set_backend_family(default_curtain_wall_family().panel, meshcat, meshcat_material_family(meshcat_glass))
set_backend_family(default_curtain_wall_family().boundary_frame, meshcat, meshcat_material_family(meshcat_generic_metal))
set_backend_family(default_curtain_wall_family().transom_frame, meshcat, meshcat_material_family(meshcat_generic_metal))
set_backend_family(default_curtain_wall_family().mullion_frame, meshcat, meshcat_material_family(meshcat_generic_metal))

realize_slab(b::MCAT, contour::ClosedPath, holes::Vector{<:ClosedPath}, level::Level, family::Family) =
  let base = vz(level.height + slab_family_elevation(b, family)),
      thickness = slab_family_thickness(b, family),
      (mattop, matbot, matside) = slab_materials(b, family_ref(b, family))
    add_object(b, meshcat_extrusion_z(path_set(contour, holes...), thickness, mattop, base))
  end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.08"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
data = CSV.read(datadir("purchaseData.csv"), copycols = true)

# ╔═╡ 9504fbf4-2441-11eb-2269-59ab4e930570
begin
	newCol = DataFrame(Validated=ones(Int, size(data,1)))
	newRow = DataFrame([["JOHN", "JACK"] [123456, 909595]], [:Name, :PhoneNo])
	newData = DataFrame(Name=["JOHN", "ASHELY","MARYANNA"],
						Job=["Lawyer", "Doctor","Lawyer"])
end

# ╔═╡ 95053800-2441-11eb-1ad6-e7f9cf8f1440
begin
	data1 = hcat(data, newCol)
	first(data1, 3)
end

# ╔═╡ 9505e3fc-2441-11eb-29de-b73f5d680690
begin
	data2 = vcat(data1, newRow, cols=:union)
	last(data2, 3)
end

# ╔═╡ 9513dd5e-2441-11eb-1aa1-2351da214449
begin
	data3 = join(data2, newData, on=:Name)
	data3
end

# ╔═╡ 95146f26-2441-11eb-0e37-337a3ffc09c2
begin
	select!(data3,[:Name,:Job])
	data3
end

# ╔═╡ 951fe570-2441-11eb-3ebd-539c016899b4
begin
	unique!(data3,:Job)
	data3
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.08"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═9504fbf4-2441-11eb-2269-59ab4e930570
# ╠═95053800-2441-11eb-1ad6-e7f9cf8f1440
# ╠═9505e3fc-2441-11eb-29de-b73f5d680690
# ╠═9513dd5e-2441-11eb-1aa1-2351da214449
# ╠═95146f26-2441-11eb-0e37-337a3ffc09c2
# ╠═951fe570-2441-11eb-3ebd-539c016899b4
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
@testset "Diagram" begin
    # Diagram = DiagTreeNew.Diagram
    # DiagramId = DiagTreeNew.DiagramId
    # add_subdiagram! = DiagTreeNew.add_subdiagram!

    struct ID <: DiagramId
        uid::Int
    end
    Base.show(io::IO, d::ID) = print(io, d.uid)
    # Base.isequal(a::ID, b::ID) = (a.index == b.index)
    # Base.Dict(d::ID) = Dict(:id => d.index)

    DiagTree.uidreset()

    W = Int
    ll = Diagram{W}(ID(3))
    l = Diagram{W}(ID(1), Sum(), [ll,])
    r = Diagram{W}(ID(2))
    root = Diagram{W}(ID(0), Sum(), [l, r])
    print_tree(root)
    """
    4 : 0=0=⨁ (2, 3)
    ├─ 2 : 1=0=⨁ (1)
    │  └─ 1 : 3=0
    └─ 3 : 2=0
    """

    collect(PostOrderDFS(root))
    @test [node.id.uid for node in PostOrderDFS(root)] == [3, 1, 2, 0]
    @test [node.id.uid for node in PreOrderDFS(root)] == [0, 1, 3, 2]
    @test [node.id.uid for node in Leaves(root)] == [3, 2]

    eval(d::ID) = d.uid
    @test evalDiagTree!(root, eval) == sum(node.id.uid for node in Leaves(root))

    print_tree(root)
    # DiagTreeNew.plot_tree(root)

    println(toDataFrame([root,]))
end

function getdiagram(spin = 2.0, D = 3, Nk = 4, Nt = 2)
    """
        k1-k3                     k2+k3
        |                         |
    t1.L ↑     t1.L       t2.L     ↑ t2.L
        |-------------->----------|
        |       |  k3+k4   |      |
        |   v   |          |  v   |
        |       |    k4    |      |
        |--------------<----------|
    t1.L ↑    t1.L        t2.L     ↑ t2.L
        |                         |
        k1                        k2
    """

    DiagTree.uidreset()
    # We only consider the direct part of the above diagram

    paraG = GenericPara(diagType = GreenDiag,
        innerLoopNum = 0, totalLoopNum = Nk, loopDim = D,
        hasTau = true, totalTauNum = Nt)
    paraV = paraG

    # #construct the propagator table
    gK = [[0.0, 0.0, 1.0, 1.0], [0.0, 0.0, 0.0, 1.0]]
    gT = [(1, 2), (2, 1)]
    g = [Diagram(BareGreenId(paraG, k = gK[i], t = gT[i]), name = :G) for i in 1:2]

    vdK = [[0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0]]
    # vdT = [[1, 1], [2, 2]]
    vd = [Diagram(BareInteractionId(paraV, ChargeCharge, k = vdK[i], permu = Di), name = :Vd) for i in 1:2]

    veK = [[1, 0, -1, -1], [0, 1, 0, -1]]
    # veT = [[1, 1], [2, 2]]
    ve = [Diagram(BareInteractionId(paraV, ChargeCharge, k = veK[i], permu = Ex), name = :Ve) for i in 1:2]

    Id = GenericId(paraV)
    # contruct the tree
    ggn = Diagram(Id, Prod(), [g[1], g[2]])
    vdd = Diagram(Id, Prod(), [vd[1], vd[2]], factor = spin)
    vde = Diagram(Id, Prod(), [vd[1], ve[2]], factor = -1.0)
    ved = Diagram(Id, Prod(), [ve[1], vd[2]], factor = -1.0)
    vsum = Diagram(Id, Sum(), [vdd, vde, ved])
    root = Diagram(Id, Prod(), [vsum, ggn], factor = 1 / (2π)^D, name = :root)

    return root, gK, gT, vdK, veK
end

@testset "Generic Diagrams" begin

    DiagTree.uidreset()
    # We only consider the direct part of the above diagram
    spin = 2.0
    D = 3
    kF, β, mass2 = 1.919, 0.5, 1.0
    Nk, Nt = 4, 2

    root, gK, gT, vdK, veK = getdiagram(spin, D, Nk, Nt)

    #optimize the diagram
    root = DiagTree.optimize(root)

    # autodiff
    droot_dg = DiagTree.derivative(root, BareGreenId)
    droot_dv = DiagTree.derivative(root, BareInteractionId)
    # plot_tree(droot_dg)

    evalDiagTree!(root, x -> 1.0)
    @test root.weight ≈ -2 + spin

    evalDiagTree!(droot_dg, x -> 1.0)
    @test root.weight ≈ (-2 + spin) * 2

    evalDiagTree!(droot_dv, x -> 1.0)
    @test root.weight ≈ (-2 + spin) * 2

    # #more sophisticated test of the weight evaluation
    varK = rand(D, Nk)
    varT = [rand() * β for t in 1:Nt]

    function evalG(K, τBasis, varT, order = 0)
        ϵ = dot(K, K) / 2 - kF^2
        τ = varT[τBasis[2]] - varT[τBasis[1]]
        if order == 0
            return Spectral.kernelFermiT(τ, ϵ, β)
        elseif order == 1
            return Spectral.kernelFermiT(τ, ϵ, β) * 3.1415
        else
            error("not implemented!")
        end
    end

    function evalV(K, order = 0)
        if order == 0
            return 8π / (dot(K, K) + mass2)
        elseif order == 1
            return 8π / (dot(K, K) + mass2) * 3.1415
        else
            error("not implemented!")
        end
    end

    # # getK(basis, varK) = sum([basis[i] * K for (i, K) in enumerate(varK)])
    getK(basis, varK) = varK * basis

    eval(id::BareGreenId, varK, varT) = evalG(getK(id.extK, varK), id.extT, varT, id.order[1])
    eval(id::BareInteractionId, varK, varT) = evalV(getK(id.extK, varK), id.order[2])

    gw = [evalG(getK(gK[i], varK), gT[i], varT) for i = 1:2]
    vdw = [evalV(getK(vdK[i], varK)) for i = 1:2]
    vew = [evalV(getK(veK[i], varK)) for i = 1:2]

    dgw = [evalG(getK(gK[i], varK), gT[i], varT, 1) for i = 1:2]
    dvdw = [evalV(getK(vdK[i], varK), 1) for i = 1:2]
    dvew = [evalV(getK(veK[i], varK), 1) for i = 1:2]

    Vweight = spin * vdw[1] * vdw[2] - vdw[1] * vew[2] - vew[1] * vdw[2]
    Gweight = gw[1] * gw[2]
    Weight = Gweight * Vweight / (2π)^D

    dVweight = spin * (dvdw[1] * vdw[2] + vdw[1] * dvdw[2]) -
               (dvdw[1] * vew[2] + vdw[1] * dvew[2]) -
               (dvew[1] * vdw[2] + vew[1] * dvdw[2])

    dGweight = dgw[1] * gw[2] + gw[1] * dgw[2]
    dWeight_dg = dGweight * Vweight / (2π)^D
    dWeight_dv = Gweight * dVweight / (2π)^D

    # print_tree(root)
    evalDiagTree!(root, eval, varK, varT)
    @test root.weight ≈ Weight

    evalDiagTree!(droot_dg, eval, varK, varT)
    @test droot_dg.weight ≈ dWeight_dg

    evalDiagTree!(droot_dv, eval, varK, varT)
    @test droot_dv.weight ≈ dWeight_dv

    ############### test diagram optimization #################
    uniqueG, uniqueInt = DiagTree.removeDuplicatedLeaves!([root,], verbose = 1)
    @test length(uniqueG) == 2
    @test length(uniqueInt) == 3
    evalDiagTree!(root, eval, varK, varT)
    @test root.weight ≈ Weight
end

@testset "optimize" begin
    DiagTree.uidreset()

    W = Int
    lll = Diagram{W}(ID(5))
    ll = Diagram{W}(ID(3), Prod(), [lll,])
    l = Diagram{W}(ID(1), Sum(), [ll,])
    r = Diagram{W}(ID(2))
    root = Diagram{W}(ID(0), Sum(), [l, r])
    # print_tree(root)
    """
    5 : 0=0=⨁ (3, 4)
    ├─ 3 : 1=0=⨁ (2)
    │  └─ 2 : 3=0=Ⓧ (1)
    │     └─ 1 : 5=0
    └─ 4 : 2=0
    """

    #remove the 2, which only has one child
    DiagTree.removeOneChildParent!([root,])
    """
    4 : 0=0=⨁ (1, 3)
    ├─ 1 : 3=0
    └─ 3 : 2=0
    """

    # print_tree(root)

    @test root.subdiagram[1].hash == 1
    # print_tree(root)

    # root, gK, gT, vdK, veK = getdiagram()
    # uniqueG, uniqueInt = DiagTree.removeDuplicatedLeaves!(root, verbose = 1)
    # @test length(uniqueG) == 2
    # @test length(uniqueInt) == 3
end


text:
export clenshaw, @clenshaw

##
# ClenshawPlan can be used for multiple evaluations with different functions
# that have the same length
##

mutable struct ClenshawPlan{S,T}
    sp::S
    bk::Vector{T}
    bk1::Vector{T}
    bk2::Vector{T}
    A::Vector{T}
    B::Vector{T}
    C::Vector{T}
end

function ClenshawPlan(::Type{T},sp,N::Int,n::Int) where T
    A = T[recA(T,sp,k) for k=0:N-1]
    B = T[recB(T,sp,k) for k=0:N-1]
    C = T[recC(T,sp,k) for k=1:N]
    ClenshawPlan(sp,Array{T}(undef,n),Array{T}(undef,n),Array{T}(undef,n),A,B,C)
end

macro clenshaw(x, c...)
    a, b = :(zero(t)), :(zero(t))
    as = []
    N = length(c)
    for k = N:-1:2
        ak = Symbol("a",k)
        push!(as, :($ak = $a))
        a = :(muladd(t,$a,$(esc(c[k]))-$b))
        b = :($ak)
    end
    ex = Expr(:block,as...,:(muladd(t/2,$a,$(esc(c[1]))-$b)))
    Expr(:block, :(t = $(esc(2))*$(esc(x))), ex)
end

clenshaw(x,c) = clenshaw_halved(2*x, c)

@generated function clenshaw_halved(x, c::StaticArrays.SVector{N,R}) where {N, R}
    a, b = :(zero(R)), :(zero(R))
    as = []
    for k = N:-1:2
        ak = Symbol("a",k)
        push!(as, :($ak = $a))
        a = :(muladd(x,$a,c[$k]-$b))
        b = :($ak)
    end
    Expr(:block,
    as...,
    :(muladd(x/2,$a,c[1]-$b)))
end

for TYP in (:AbstractVector,:AbstractMatrix)
    @eval clenshaw(c::$TYP,x::AbstractArray,plan::ClenshawPlan) = reshape(clenshaw(c,vec(x),plan),size(x))
end

function clenshaw(c::AbstractVector,x::AbstractVector,plan::ClenshawPlan{S,V}) where {S,V}
    N,n = length(c),length(x)
    if isempty(c)
        return zeros(V,length(x))
    end

    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2
    A=plan.A
    B=plan.B
    C=plan.C

    @inbounds for i = 1:n
        bk1[i] = zero(V)
        bk2[i] = zero(V)
    end

    @inbounds for k = N:-1:2
        ck,Ak,Bk,Ck = c[k],A[k],B[k],C[k]
        for i = 1:n
            bk[i] = muladd(muladd(Ak,x[i],Bk),bk1[i],muladd(-Ck,bk2[i],ck))
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    ck,Ak,Bk,Ck = c[1],A[1],B[1],C[1]
    @inbounds for i = 1:n
        bk[i] = muladd(muladd(Ak,x[i],Bk),bk1[i],muladd(-Ck,bk2[i],ck))
    end

    bk
end

#Clenshaw routine for many Funs, x is a vector of same number of funs
#each fun is a column

function clenshaw(c::AbstractMatrix,x,plan::ClenshawPlan{S,V}) where {S,V}
    N,n = size(c)
    if isempty(c)
        return zeros(V,n)
    end

    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2
    A=plan.A
    B=plan.B
    C=plan.C

    @inbounds for i = 1:n
        bk1[i] = zero(V)
        bk2[i] = zero(V)
    end

    @inbounds for k = N:-1:2
        Ak,Bk,Ck = A[k],B[k],C[k]
        for i = 1:n
            cki = c[k,i]
            bk[i] = muladd(muladd(Ak,x,Bk),bk1[i],muladd(-Ck,bk2[i],cki))
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    Ak,Bk,Ck = A[1],B[1],C[1]
    @inbounds for i = 1:n
        ci = c[1,i]
        bk[i] = muladd(muladd(Ak,x,Bk),bk1[i],muladd(-Ck,bk2[i],ci))
    end

    bk
end

function clenshaw(c::AbstractMatrix,x::AbstractVector,plan::ClenshawPlan{S,V}) where {S,V}
    N,n = size(c)
    @assert n == length(x)
    if isempty(c)
        return zeros(V,length(x))
    end

    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2
    A=plan.A
    B=plan.B
    C=plan.C

    @inbounds for i = 1:n
        bk1[i] = zero(V)
        bk2[i] = zero(V)
    end

    @inbounds for k = N:-1:2
        Ak,Bk,Ck = A[k],B[k],C[k]
        for i = 1:n
            cki = c[k,i]
            bk[i] = muladd(muladd(Ak,x[i],Bk),bk1[i],muladd(-Ck,bk2[i],cki))
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    Ak,Bk,Ck = A[1],B[1],C[1]
    @inbounds for i = 1:n
        ci = c[1,i]
        bk[i] = muladd(muladd(Ak,x[i],Bk),bk1[i],muladd(-Ck,bk2[i],ci))
    end

    bk
end

# overwrite x

function clenshaw!(c::AbstractVector,x::AbstractVector,plan::ClenshawPlan{S,V}) where {S,V}
    N,n = length(c),length(x)
    if isempty(c)
        for k=1:n
            x[k]=zero(V)
        end
    end

    bk=plan.bk
    bk1=plan.bk1
    bk2=plan.bk2
    A=plan.A
    B=plan.B
    C=plan.C

    @inbounds for i = 1:n
        bk1[i] = zero(V)
        bk2[i] = zero(V)
    end

    @inbounds for k = N:-1:2
        ck,Ak,Bk,Ck = c[k],A[k],B[k],C[k]
        for i = 1:n
            bk[i] = muladd(muladd(Ak,x[i],Bk),bk1[i],muladd(-Ck,bk2[i],ck))
        end
        bk2, bk1, bk = bk1, bk, bk2
    end

    ck,Ak,Bk,Ck = c[1],A[1],B[1],C[1]
    @inbounds for i = 1:n
        x[i] = muladd(muladd(Ak,x[i],Bk),bk1[i],muladd(-Ck,bk2[i],ck))
    end

    x
end

function sineshaw(c::AbstractVector,θ)
    if isempty(c)
        return zero(θ)
    end

    x = 2cos(θ)
    bk1,bk2 = zero(x),zero(x)
    @inbounds for k = length(c):-1:1
        bk2, bk1 = bk1, muladd(x,bk1,c[k]-bk2)
    end

    sin(θ)*bk1
end
sineshaw(c::AbstractVector,θ::AbstractVector) = promote_type(eltype(c),eltype(θ))[sineshaw(c,θ[k]) for k=1:length(θ)]
sineshaw(c::AbstractVector,θ::AbstractMatrix) = promote_type(eltype(c),eltype(θ))[sineshaw(c,θ[k,j]) for k=1:size(θ,1),j=1:size(θ,2)]


text:
#=
choose the frequency of experiments to minimize a given convex spectral function
of the information matrix and satisfy an experiment budget constraint

adapted from Boyd and Vandenberghe, "Convex Optimization", section 7.5

min f(V × Diagonal(x) × V')
subject to:
x ∈ ℤᵏ
0 ≤ x ≤ xmax
∑x = xsum

where:
- there are k = 2d experiments
- the menu of experiments is V ∈ ℝ^(d × k)
- variable x is the frequency of each experiment
- f is a convex spectral function
- xsum = k is the total number of experiments to run
- xmax = 2 is the maximum frequency of each experiment
=#

struct ExperimentDesign <: ExampleInstance
    d::Int
    f::MatSpecExt # formulation specifier
end

function build(inst::ExperimentDesign)
    d = inst.d
    @assert d >= 1
    @assert is_domain_pos(inst.f)
    k = 2 * d

    # generate data
    V = randn(d, k)
    V .*= sqrt(d / sum(eigvals(Symmetric(V * V'))))

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, x[1:k] >= 0)
    JuMP.@constraint(model, x .<= 2)
    JuMP.set_integer.(x)
    JuMP.@constraint(model, sum(x) == k)

    # vectorized information matrix
    Q_vec = [
        JuMP.@expression(
            model,
            (i == j ? 1.0 : rt2) * sum(V[i, k] * x[k] * V[j, k] for k in 1:k)
        ) for i in 1:d for j in 1:i
    ]

    # convex objective
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)
    add_homog_spectral(inst.f, d, vcat(1.0 * epi, Q_vec), model)

    # save for use in tests
    model.ext[:V] = V
    model.ext[:x] = x

    return model
end

function test_extra(inst::ExperimentDesign, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    tol = eps()^0.2
    x = JuMP.value.(model.ext[:x])
    @test x ≈ round.(Int, x) atol = tol rtol = tol
    @test all(0 .<= x .<= 2)
    # check objective
    V = model.ext[:V]
    λ = eigvals(Symmetric(V * Diagonal(x) * V', :U))
    @test minimum(λ) >= -tol
    obj_result = get_val(pos_only(λ), inst.f)
    @test JuMP.objective_value(model) ≈ obj_result atol = tol rtol = tol
    return
end


text:
#PCC
function ccQuick(k, l=1)
	k=map(deg2rad, k)
	n=length(k) #get num segemnts
	p=zeros(n+1, 2) #[x, y]
	t=zeros(n+1) #[angle]
	#calculate profile
	for i=1:n
		t[i+1]=t[i]+k[i]*l
		if k[i]==0 #handles 0 curvature
			p[i+1, 1]=p[i, 1]+cos(t[i+1])*l #update x
	        p[i+1, 2]=p[i, 2]+sin(t[i+1])*l #update y
		else
			p[i+1, 1]=p[i, 1]+(sin(t[i+1])-sin(t[i]))/k[i] #update x
	        p[i+1, 2]=p[i, 2]+(-cos(t[i+1])+cos(t[i]))/k[i] #update y
	    end
	end
	return p
end

function ccEnd(k, l=1)
	return ccQuick(k, l)[end,:]
end

function ccArm2(k, res=10, l=1)
	l=l/res
	k=map(deg2rad, k)
	n=length(k) #get num segemnts
	k=vec(repmat(k', res )) #set curv for each subsegment
	p=zeros(n*res+1, 2) #[x, y]
	t=zeros(n*res+1) #[angle]
	#calculate profile
	for i=1:(n*res)
		t[i+1]=t[i]+k[i]*l
		if k[i]==0 #handles 0 curvature
			p[i+1, 1]=p[i, 1]+cos(t[i+1])*l #update x
	        p[i+1, 2]=p[i, 2]+sin(t[i+1])*l #update y
		else
			p[i+1, 1]=p[i, 1]+(sin(t[i+1])-sin(t[i]))/k[i] #update x
	        p[i+1, 2]=p[i, 2]+(-cos(t[i+1])+cos(t[i]))/k[i] #update y
	    end
	end
	e=zeros(n+1,2) #get edge values
	for i=1:n+1
		e[i,:]=p[(i-1)*res+1,:]
	end
	#println("Generated Arm wth $(n) segments")
	return (p, e)
end

function drawArm(arm, edge, res=10)
	gcol=["c","b","r", "g", "m"]
	#plot end eff
	plot([arm[end,1]],[arm[end, 2]], "k*")

	#plot profile and edges
	for i=1:nSeg
	    inx=(i-1)*res+1
	    plot(arm[inx:inx+res, 1],arm[inx:inx+res, 2], gcol[i%5+1])
	    plot(edge[i, 1], edge[i, 2], gcol[i%5+1]*"o")
	end
end
function drawArm(state)
	(p,e)=ccArm2(state)
	drawArm(p, e)
end

function selfCollide(p, thre)
	n=size(p,1)
	for i=1:n-1
		for j=i+1:n
			#@show norm(p[i,:]- p[j,:])
			if norm(p[i,:]- p[j,:])<thre
				println("self collsions at: $(i)<->$(j)")
				col=(p[j,:]+p[i,:])/2
				plot(col[1], col[2], "k*")
				return col
			end
		end
	end
	println("no self collsions founds")
	return 0
end

type Zone
	cent
	radius
end

function drawZone(zone::Zone)
	t=[linspace(0,2*pi,10)]
	p=zeros(length(t),2)
	for i=1:length(t)
		p[i,:]=zone.cent+zone.radius*[cos(t[i]) sin(t[i])]
	end
	plot(p[:,1], p[:,2], "ko")
	return p
end

function obsCollide(arm, zone::Zone, thre=1.1, draw=0)
	for i=1:size(arm, 1)
		if norm(arm[i, :]-zone.cent)<(zone.radius*thre)
			if draw==1
				println("obstacle collision at $(arm[i,:])")
				plot(arm[i, 1], arm[i, 2], "k*")
			end
			return arm[i,:]
		end
	end
	println("no obstacle collsions founds")
	return 0
end

function zoneDist(s, zone::Zone)
    return norm(zone.cent- ccEnd(s))
end

###########broken or unused

function ccSeg(k, res, l)
	t=linspace(0, k*l, res+1) #calcuate angle change
	p=zeros(res+1, 2)
	y=zeros(res+1)
	for i=1:res
		p[i+1, 1]=p[i, 1]+(sin(t[i+1])-sin(t[i]))/k #update x
        p[i+1, 2]=p[i, 2]+(-cos(t[i+1])+cos(t[i]))/k #update y
	end
	return (p[:, 1], p[:, 2], t[end])
end

function ccArm(k, res, l=1)
	n=length(k) #number of segments
    #init segments cords
    x=zeros(n, res+1)
    y=zeros(n, res+1)
    t=zeros(n)
    #conected arm

    #calculate segment profile
    for i=1:n
        (x[i,:], y[i,:], t[i])=ccSeg(k[i], res, l)
    end
	arm=posArm(x, y, t)

	return (arm, t)
end

##postion arm segments
function posArm(xSeg, ySeg, t)
	arm=zeros(n*res+1, 2) #[x, y]
	#connect segments
    arm[1:res+1,:]=hcat(xSeg[1,:]', ySeg[1,:]') #set 1st segment
    for i=2:n
        idx=(i-1)*res+1
        @show sum(t[1:(i)])
        st=rotate(sum(t[1:(i-1)]))*vcat(xSeg[1,:], ySeg[1,:]) #rotate segment
        arm[idx:idx+res,:]=st'.+arm[idx,:] #translate segment
    end

    return arm
end
function rotate(t)
    return [cos(t) -sin(t); sin(t) cos(t)]
end

#cant
function curv(x, q, l=1, EI=1)
    return atan(q*(3*l^2*x-3*l*x^2+x^3)/(6*EI))
end

function deflec(x, q, l=1, EI=1)
    return q*x^2*(6*l^2-4*l*x+x^2)/(24*EI)
end

function cantSeg(q, res, l=1, EI=1)
    x=linspace(0, l, res+1)
    y=zeros(length(x))
    for i=1:length(x)
        y[i]=deflec(x[i], q)
    end
    t=curv(l, q)
    return (x, y, t)
end

function smash(e)
	n=size(e,1)
	e=e'
	#checks for colisions
	for i=1:(n-2)
		for j=(i+1):(n-1)
			t=inv(hcat(-(e[:,i+1]-e[:,i]), e[:,j+1]-e[:,j]))*e[:,i]-e[:,j]
			col=vcat(t.<0,t.>1 )
			if norm(col)==0
				println("crash at: ", i, j)
				return e[:,i]+t[1]*(e[:,i+1]-e[:,i])

			end
		end
	end
    return 0
end

println("Import successful: sim.jl")


text:
#=
 # LinkedList of not yet treated vertices
 =#

mutable struct _ListEle{T <: Integer}
    vertex::WordVertex{T}
    next::Union{_ListEle{T}, Nothing}

    # _ListEle constructors
    _ListEle{T}(vertex, next) where {T <: Integer}  = new(vertex, next)
end # struct

mutable struct _UninitializedList{T <: Integer}
    n::Int32
    head::Union{_ListEle{T}, Nothing}

    # _UninitializedList constructors
    _UninitializedList{T}() where {T <: Integer}    = new(0, nothing)
end

function _push!(l::_UninitializedList{T}, v::WordVertex{T}) where {T <: Integer}
    l.n += 1
    l.head = _ListEle{T}(v, l.head) # Create list element
end

function _pop!(l::_UninitializedList{T}) where {T <: Integer}
    if (isnothing(l.head))
        error("pop!: Can't pop from an empty list")
        return
    end

    l.n     -= 1
    v       = l.head.vertex
    l.head  = l.head.next
    return v
end


text:
using ucbature

include("integrands.jl")

const dim = 2
const nb_runs = 10
const nb_evals = 1e4
const xmin = fill(-1.0, dim)
const xmax = fill( 1.0, dim)
const nb_zones_per_dim = 5
const k_exploration = 0.01
const nb_init_evals_per_zone = 10

my_mc(integrand) = ucbature.mc(integrand, xmin, xmax, nb_evals)
my_str(integrand) = ucbature.str(integrand, xmin, xmax, nb_evals, nb_zones_per_dim)
my_ucb(integrand) = ucbature.ucb(integrand, xmin, xmax, nb_evals, nb_zones_per_dim,
                                 k_exploration, nb_init_evals_per_zone)

my_run(integrator, integrand) = mapreduce((x)->integrator(integrand), +, 1:nb_runs)/nb_runs

println("sphere, mc: ", my_run(my_mc, integrands.sphere))
println("sphere, str: ", my_run(my_str, integrands.sphere))
println("sphere, ucb: ", my_run(my_ucb, integrands.sphere))

println("linear, mc: ", my_run(my_mc, integrands.linear))
println("linear, str: ", my_run(my_str, integrands.linear))
println("linear, ucb: ", my_run(my_ucb, integrands.linear))

max_stars_repo_name: florianLepretre/ucbature

max_issues_repo_name: florianLepretre/ucbature

max_forks_repo_name: florianLepretre/ucbature


text:
#module HARPSNData
#export read_data

import ..EchelleObservation, ..EchelleObservationSet, ..OrdersType
import ..AbstractInstrument
#, ..HARPSNType
#import ..harpsn_flux_hdu,  ..harpsn_all_orders
using FITSIO

"""
read_data(filename; orders)
Example:
   data = read_data(filename_list[1])
"""
function read_data(filename::String; orders::OrdersType = harpsn_all_orders )
   local fits
   try
      fits = FITS(filename)
   catch
      @error("# Can't read " * filename)
   end
   if typeof(orders) <: Integer
      orders = orders:orders
   end
   hdr = read_header(fits[1])
   # Get header info for estimating flux uncertainties
   gain = hdr["TNG DRS CCD CONAD"]
   ron = hdr["TNG DRS CCD SIGDET"]
   # Get barycentric Earth RV correction (should be positive in morning when observer is approaching sun)
   berv = hdr["TNG DRS BERV"]
   berv /= speed_of_light_mps
   flux = read(fits[harpsn_flux_hdu])[:,orders]
   flux *= gain # in ADU
   #flux .*= gain
   var = copy(flux)
   flux[flux .< zero(eltype(flux))] .= zero(eltype(flux))
   var .+= ron^2
   npix,nord = size(flux)
   # Extract pixel->wavelegnth polynomials (thank ACC)
   coeff = zeros(4,nord)
   for ord in 0:(nord-1)
      for poly in 0:3
         n = ord*4+poly
         coeff[1+poly,1+ord] = hdr["TNG DRS CAL TH COEFF LL" * string(n)]
         #push!(coeff,hdr["TNG DRS CAL TH COEFF LL" * string(n)])
         #n += 1
      end
   end
   xpix = 0:(npix-1)
   xarray =  vcat(ones(length(xpix))', xpix', (xpix.^2)', (xpix.^3)' )
   lambda = zeros(npix,length(orders))
   for ord in orders
      lambda[:,ord] .= xarray' * coeff[:,ord]
   end
   metadata = Dict{Symbol,Any}(:Filename=>filename,:Instrument=>:HARPSN,:z_bc=>berv)
#   return EchelleObservation{eltype(lambda),eltype(flux),eltype(var),typeof(berv),HARPSNType}(flux=flux,var=var,lambda=lambda,z_bc=berv,metadata=metadata)
   return EchelleObservation(flux=flux,var=var,lambda=lambda,z_bc=berv,instrument=HARPSNType(), metadata=metadata)
end

"""
   read_data!(data_out, obsid, filename; orders)
Writes to data_out[:,:,i]'s flux, var and lambda.
Creates metadata_list including :Filename and :z_bc (redshift of observer relative to solar system barcenter)

Example:
   data = read_data!(data_out, obsid, filename, orders=10:15)
"""
function read_data!(data::EchelleObservationSet{T1,T2,T3,#= T4, =# HARPSNType}, i::Int, filename::String;
      orders::OrdersType = harpsn_all_orders
      ) where {T1<:Number, T2<:Number, T3<:Number, #=T4<:Number,  HARPSNType <:AbstractInstrument =#}
   @assert 1 <= i <= size(data.flux,3)
   local fits
   try
      fits = FITS(filename)
   catch
      @error("# Can't read " * filename)
   end
   if typeof(orders) <: Integer
      orders = orders:orders
   end
   obs = read_data(filename, orders=orders)
   data.flux[:,:,i] .= obs.flux
   data.var[:,:,i] .= obs.var
   data.lambda[:,:,i] .= obs.lambda
   data.metadata_list[i] = Dict{Symbol,Any}(:Filename=>filename,:z_bc=>obs.metadata[:z_bc])
   return data
end

"""
   read_data(filenames; orders, z_bsrv)
Returns struct with flux, var, lambda, and metadata_list including :Filename and :z_bc
If pass z_bsrv, then adds that to the barrycentric correction stored in metadata_list.

Example:
   data = read_data(filename_list, orders=10:15)
"""
function read_data(filename_list::AbstractArray{String,1};
      orders::OrdersType = harpsn_all_orders,
      z_bsrv::AbstractArray = zeros(length(filename_list)) )
   @assert length(filename_list) >= 1
   data_first = read_data(filename_list[1], orders=orders)
   data = EchelleObservationSet{Float64,Float64,Float64,HARPSNType}(size(data_first.flux), length(filename_list), metadata=Dict(:Instrument=>:HARPSN) )
   data.lambda[:,:,1] = data_first.lambda
   data.flux[:,:,1]   = data_first.flux
   data.var[:,:,1]    = data_first.var
   for (i,filename) in enumerate(filename_list)
          read_data!(data, i, filename, orders=orders)
          data.metadata_list[i][:z_bc] += z_bsrv[i]
   end
   return data
end

default_harpsn_blaze_file = joinpath(homedir(),"Downloads","HARPN.2015-07-31T17-02-35.072_blaze_A.fits")
function read_blaze_list(filename::String = default_harpsn_blaze_file,
      pixel_buffer_lo::Integer = default_pixel_buffer_lo,
      pixel_buffer_hi::Integer = default_pixel_buffer_hi,
      log_lambda_lo::AbstractArray{Float64,1} = zeros(length(harpsn_all_orders)),
      log_lambda_hi::AbstractArray{Float64,1} = zeros(length(harpsn_all_orders)) )
   fits = FITS(filename)
   blaze_data = read(fits[1])
   pixels_per_order, num_orders = size(blaze_data)
   pixel_range = (1+default_pixel_buffer_lo):(pixels_per_order-default_pixel_buffer_lo)
   blaze_list = Array{BlazeMap1D{Int64,Float64,Float32},1}(undef,num_orders)
   @assert length(log_lambda_lo) == length(log_lambda_hi) == num_orders
   for i in 1:num_orders
      blaze_list[i] = BlazeMap1D(i,log_lambda_lo[i],log_lambda_hi[i],pixel_range, view(blaze_data,:,i) )
   end

   return BlazeMap2D(blaze_list)
end
export read_blaze_list

function read_blaze_list(data::EchelleObservationSet, filename::String = default_harpsn_blaze_file,
      pixel_buffer_lo::Integer = default_pixel_buffer_lo,
      pixel_buffer_hi::Integer = default_pixel_buffer_hi
      )
   pixels_per_order, num_orders, num_obs = size(data.lambda)
   pixel_range = (1+default_pixel_buffer_lo):(pixels_per_order-default_pixel_buffer_lo)
   log_lambda_lo = map(i->log(minimum(data.lambda[1,i,:])),1:num_orders)
   log_lambda_hi = map(i->log(maximum(data.lambda[end,i,:])),1:num_orders)
   read_blaze_list(filename,pixel_buffer_lo,pixel_buffer_hi, log_lambda_lo, log_lambda_hi )
 end

#end # module


text:
using Flight
using GLFW

function test_xbox()

    #detect joysticks already connected
    init_joysticks()

    #enable detection of further joystick connections or disconnections
    GLFW.SetJoystickCallback(joystick_callback)
    GLFW.PollEvents() #check for newly connected joysticks

    run_time = 10.0
    dt = 0.5

    t_wall_0 = time()
    t_wall = t_wall_0
    while t_wall - t_wall_0 < run_time

        #compute the wall time corresponding to the newly updated simulation
        t_wall_next = t_wall + dt

        #busy wait until the wall time catches up
        while (time() < t_wall_next) end

        t_wall = t_wall_next

        #try to catch potential joystick disconnections
        GLFW.PollEvents()
        for joystick in values(connected_joysticks)
            update_joystick(joystick)
            println(get_axis_data(joystick))
            println()
            println(get_button_state(joystick, :button_A))
            println(get_button_change(joystick, :button_A))
        end

    end

end


text:
#####
# StandardLine vs. StandardLine
#####

function is_colliding(a::StandardLine{T}, b::StandardLine{T}, pos_ba::Vector2D{T}) where {T}
    y = pos_ba[2]
    if iszero(y)
        x = pos_ba[1]
        half_length_a = get_half_length(a)
        half_length_b = get_half_length(b)
        return !((x + half_length_b <= -half_length_a) || (x - half_length_b >= half_length_a))
    else
        return false
    end
end

function is_colliding(l1::StandardLine{T}, l2::StandardLine{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T}
    half_length = get_half_length(l1)

    tail_l2, head_l2 = get_vertices(l2, pos, dir)

    x1 = tail_l2[1]
    y1 = tail_l2[2]

    x2 = head_l2[1]
    y2 = head_l2[2]

    t = y1 / (y1 - y2)

    if isfinite(t)
        x_intersection = x1 + t * (x2 - x1)
        return (zero(t) < t < one(t)) && (-half_length < x_intersection < half_length)
    else
        return (iszero(y1)) && (iszero(y2)) && !((half_length <= x1) || (x2 <= -half_length))
    end
end

#####
# StandardCircle vs. StandardPoint
#####

function is_inside(circle::StandardCircle{T}, pos::Vector2D{T}) where {T}
    radius = get_radius(circle)
    return LA.dot(pos, pos) < radius * radius
end

is_colliding(circle::StandardCircle{T}, point::StandardPoint{T}, pos::Vector2D{T}) where {T} = is_inside(circle, pos)
is_colliding(point::StandardPoint{T}, circle::StandardCircle{T}, pos::Vector2D{T}) where {T} = is_inside(circle, pos) # no need to reverse pos because of symmetry

is_colliding(circle::StandardCircle{T}, point::StandardPoint{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_inside(circle, pos)
is_colliding(point::StandardPoint{T}, circle::StandardCircle{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_inside(circle, pos) # no need to reverse pos because of symmetry

#####
# StandardCircle vs. StandardLine
#####

function get_projection(line::StandardLine{T}, pos::Vector2D{T}) where {T}
    half_length = get_half_length(line)
    x = pos[1]
    if x < -half_length
        return get_tail(line)
    elseif x > half_length
        return get_head(line)
    else
        return typeof(pos)(x, zero(x))
    end
end

function is_colliding(line::StandardLine{T}, circle::StandardCircle{T}, pos::Vector2D{T}) where {T}
    projection = get_projection(line, pos)
    vec = pos - projection
    radius = get_radius(circle)
    return LA.dot(vec, vec) < radius * radius
end

is_colliding(circle::StandardCircle{T}, line::StandardLine{T}, pos::Vector2D{T}) where {T} = is_colliding(line, circle, pos) # no need to reverse pos because of symmetry

is_colliding(line::StandardLine{T}, circle::StandardCircle{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_colliding(line, circle, pos)
is_colliding(circle::StandardCircle{T}, line::StandardLine{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_colliding(line, circle, invert(pos, dir)...)

#####
# StandardCircle vs. StandardCircle
#####

function is_colliding(c1::StandardCircle{T}, c2::StandardCircle{T}, pos::Vector2D{T}) where {T}
    r1 = get_radius(c1)
    r2 = get_radius(c2)
    r = r1 + r2
    return LA.dot(pos, pos) < r * r
end

is_colliding(c1::StandardCircle{T}, c2::StandardCircle{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_colliding(c1, c2, pos)

#####
# StandardRect vs. StandardPoint
#####

function is_inside(rect::StandardRect{T}, pos::Vector2D{T}) where {T}
    half_width = get_half_width(rect)
    half_height = get_half_height(rect)

    x = pos[1]
    y = pos[2]

    return (-half_width < x < half_width) && (-half_height < y < half_height)
end

is_colliding(rect::StandardRect{T}, point::StandardPoint{T}, pos::Vector2D{T}) where {T} = is_inside(rect, pos)
is_colliding(point::StandardPoint{T}, rect::StandardRect{T}, pos::Vector2D{T}) where {T} = is_inside(rect, pos) # no need to reverse pos because of symmetry

is_colliding(rect::StandardRect{T}, point::StandardPoint{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_inside(rect, pos)
is_colliding(point::StandardPoint{T}, rect::StandardRect{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_colliding(rect, point, invert(pos, dir)...)

#####
# StandardRect vs. StandardLine
#####

function separating_axis_exists(rect::StandardRect{T}, line::StandardLine{T}, pos::Vector2D{T}) where {T}
    half_height = get_half_height(rect)
    half_width = get_half_width(rect)

    x = pos[1]
    y = pos[2]
    half_length = get_half_length(line)

    return (y <= -half_height) || (y >= half_height) || (x + half_length <= -half_width) || (x - half_length >= half_width)
end

is_colliding(rect::StandardRect{T}, line::StandardLine{T}, pos::Vector2D{T}) where {T} = !separating_axis_exists(rect, line, pos)
is_colliding(line::StandardLine{T}, rect::StandardRect{T}, pos::Vector2D{T}) where {T} = is_colliding(rect, line, pos) # no need to reverse pos because of symmetry

function separating_axis_exists(rect::StandardRect{T}, line::StandardLine{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T}
    half_width = get_half_width(rect)
    half_height = get_half_height(rect)

    tail, head = get_vertices(line, pos, dir)

    min_x, max_x = minmax(tail[1], head[1])
    min_y, max_y = minmax(tail[2], head[2])

    return (max_x <= -half_width) || (min_x >= half_width) || (max_y <= -half_height) || (min_y >= half_height)
end

function separating_axis_exists(line::StandardLine{T}, rect::StandardRect{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T}
    half_length = get_half_length(line)

    bottom_left, bottom_right, top_right, top_left = get_vertices(rect, pos, dir)

    min_x, max_x = extrema((bottom_left[1], bottom_right[1], top_right[1], top_left[1]))
    min_y, max_y = extrema((bottom_left[2], bottom_right[2], top_right[2], top_left[2]))

    return (max_x <= -half_length) || (min_x >= half_length) || (max_y <= zero(T)) || (min_y >= zero(T))
end

is_colliding(rect::StandardRect{T}, line::StandardLine{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = !(separating_axis_exists(rect, line, pos, dir) || separating_axis_exists(line, rect, invert(pos, dir)...))
is_colliding(line::StandardLine{T}, rect::StandardRect{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = !(separating_axis_exists(line, rect, pos, dir) || separating_axis_exists(rect, line, invert(pos, dir)...))

#####
# StandardRect vs. StandardCircle
#####

get_projection(rect::StandardRect{T}, pos::Vector2D{T}) where {T} = clamp.(pos, get_bottom_left(rect), get_top_right(rect))

function is_colliding(rect::StandardRect{T}, circle::StandardCircle{T}, pos::Vector2D{T}) where {T}
    projection = get_projection(rect, pos)
    vec = pos - projection
    radius = get_radius(circle)
    return LA.dot(vec, vec) < radius * radius
end

is_colliding(circle::StandardCircle{T}, rect::StandardRect{T}, pos::Vector2D{T}) where {T} = is_colliding(rect, circle, pos) # no need to invert pos because of symmetry

is_colliding(rect::StandardRect{T}, circle::StandardCircle{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_colliding(rect, circle, pos)
is_colliding(circle::StandardCircle{T}, rect::StandardRect{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = is_colliding(rect, circle, invert(pos, dir)...)

#####
# StandardRect vs. StandardRect
#####

function separating_axis_exists(r1::StandardRect{T}, r2::StandardRect{T}, pos::Vector2D{T}) where {T}
    half_width_r1 = get_half_width(r1)
    half_height_r1 = get_half_height(r1)

    half_width_r2 = get_half_width(r2)
    half_height_r2 = get_half_height(r2)

    x = pos[1]
    y = pos[2]

    return (x + half_width_r2 <= -half_width_r1) || (x - half_width_r2 >= half_width_r1) || (y + half_height_r2 <= -half_height_r1) || (y - half_height_r2 >= half_height_r1)
end

is_colliding(a::StandardRect{T}, b::StandardRect{T}, pos_ba::Vector2D{T}) where {T} = !separating_axis_exists(a, b, pos_ba)

function separating_axis_exists(r1::StandardRect{T}, r2::StandardRect{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T}
    half_width_r1 = get_half_width(r1)
    half_height_r1 = get_half_height(r1)

    bottom_left_r2, bottom_right_r2, top_right_r2, top_left_r2 = get_vertices(r2, pos, dir)

    min_x_r2, max_x_r2 = extrema((bottom_left_r2[1], bottom_right_r2[1], top_right_r2[1], top_left_r2[1]))
    min_y_r2, max_y_r2 = extrema((bottom_left_r2[2], bottom_right_r2[2], top_right_r2[2], top_left_r2[2]))

    return ((half_width_r1 <= min_x_r2) || (max_x_r2 <= -half_width_r1) || (half_height_r1 <= min_y_r2) || (max_y_r2 <= -half_height_r1))
end

is_colliding(r1::StandardRect{T}, r2::StandardRect{T}, pos::Vector2D{T}, dir::Vector2D{T}) where {T} = !(separating_axis_exists(r1, r2, pos, dir) || separating_axis_exists(r2, r1, invert(pos, dir)...))


text:
# Define a linear model, cross validate it and add some synthetic data.
#

abstract type ExampleModel end


text:
"""
```
sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) + sum(pconvac[c] for c in bus_convs) - pd - gs*1^2
sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) + sum(qconvac[c] for c in bus_convs) - qd + bs*1^2
```
"""
function constraint_kcl_shunt(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs)
    p = PowerModels.var(pm, n, cnd, :p)
    pg = PowerModels.var(pm, n, cnd, :pg)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac)
    pconv_grid_ac = PowerModels.var(pm, n, cnd, :pconv_tf_fr)
    v = 1

    PowerModels.con(pm, n, cnd, :kcl_p)[i] = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac)  == sum(pg[g] for g in bus_gens) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end

"""
Creates Ohms constraints for DC branches

```
p[f_idx] + p[t_idx] == 0)
```
"""
function constraint_ohms_dc_branch(pm::AbstractDCPModel, n::Int, cnd::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = PowerModels.var(pm, n, cnd, :p_dcgrid, f_idx)
    p_dc_to = PowerModels.var(pm, n, cnd, :p_dcgrid, t_idx)

    @constraint(pm.model, p_dc_fr + p_dc_to == 0)
end
"`vdc[i] == vdcm`"
function constraint_dc_voltage_magnitude_setpoint(pm::AbstractDCPModel, n::Int, cnd::Int, i, vdcm)
    # not used
end

function variable_dcgrid_voltage_magnitude(pm::AbstractDCPModel; kwrags...)
    # not used nw::Int=pm.cnw
end

function constraint_dc_branch_current(pm::AbstractDCPModel, n::Int, cnd::Int, f_bus, f_idx, ccm_max, p)
# do nothing
end


text:
@testset "Script" begin
    @testset "parse" begin
        script_pubkey = IOBuffer(hex2bytes("6a47304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a7160121035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937"))
        script = Bitcoin.scriptparse(script_pubkey)
        want = hex2bytes("304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a71601")
        @test script.instructions[1] == want
        want = hex2bytes("035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937")
        @test script.instructions[2] == want
        script_pubkey = IOBuffer(hex2bytes("fdfe0000483045022100c679944ff8f20373685e1122b581f64752c1d22c67f6f3ae26333aa9c3f43d730220793233401f87f640f9c39207349ffef42d0e27046755263c0a69c436ab07febc01483045022100eadc1c6e72f241c3e076a7109b8053db53987f3fcc99e3f88fc4e52dbfd5f3a202201f02cbff194c41e6f8da762e024a7ab85c1b1616b74720f13283043e9e99dab8014c69522102b0c7be446b92624112f3c7d4ffc214921c74c1cb891bf945c49fbe5981ee026b21039021c9391e328e0cb3b61ba05dcc5e122ab234e55d1502e59b10d8f588aea4632102f3bd8f64363066f35968bd82ed9c6e8afecbd6136311bb51e91204f614144e9b53aeffffffff05a08601000000000017a914081fbb6ec9d83104367eb1a6a59e2a92417d79298700350c00000000001976a914677345c7376dfda2c52ad9b6a153b643b6409a3788acc7f341160000000017a914234c15756b9599314c9299340eaabab7f1810d8287c02709000000000017a91469be3ca6195efcab5194e1530164ec47637d44308740420f00000000001976a91487fadba66b9e48c0c8082f33107fdb01970eb80388ac00000000"))
        want = hex2bytes("522102b0c7be446b92624112f3c7d4ffc214921c74c1cb891bf945c49fbe5981ee026b21039021c9391e328e0cb3b61ba05dcc5e122ab234e55d1502e59b10d8f588aea4632102f3bd8f64363066f35968bd82ed9c6e8afecbd6136311bb51e91204f614144e9b53ae")
        script = Bitcoin.scriptparse(script_pubkey)
        @test script.instructions[4] == want
    end
    @testset "Serialize" begin
        want = "6a47304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a7160121035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937"
        script_pubkey = IOBuffer(hex2bytes(want))
        script = Bitcoin.scriptparse(script_pubkey)
        @test bytes2hex(Bitcoin.scriptserialize(script)) == want
    end
    @testset "Evaluate" begin
        modified_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000001976a914a802fc56c704ce87c42d7c92eb75e7896bdc41ae88acfeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac1943060001000000")
        h256 = Bitcoin.hash256(modified_tx)
        z = bytes2int(h256)
        stream = IOBuffer(hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600"))
        tx = txparse(stream)
        tx_in = tx.tx_ins[1]
        combined_script = Bitcoin.Script(nothing)
        append!(combined_script.instructions, copy(tx_in.script_sig.instructions))
        append!(combined_script.instructions, copy(Bitcoin.txin_scriptpubkey(tx_in).instructions))
        @test Bitcoin.scriptevaluate(combined_script, z) == true
    end
end