text:
using Interpolations

export Shape,
       Path,
       backend,
       backend_name,
       current_backend,
       has_current_backend,
       switch_to_backend,
       void_ref,
       delete_shape, delete_shapes,
       delete_all_shapes,
       set_length_unit,
       collecting_shapes,
       collected_shapes,
       with_transaction,
       surface_boundary,
       curve_domain,
       surface_domain,
       get_layer,
       create_layer,
       get_or_create_layer,
       current_layer,
       get_material,
       create_material,
       get_or_create_material,
       current_material,
       create_block,
       instantiate_block,
       reset_backend,
       connection,
       Backend,
       @deffamily,
       @defproxy,
       dimension,
       force_creation,
       subpath,
       subpath_starting_at,
       subpath_ending_at,
       bounding_box,
       capture_shape, capture_shapes,
       captured_shape, captured_shapes,
       revolve

#Backends are types parameterized by a key identifying the backend (e.g., AutoCAD) and by the type of reference they use

abstract type Backend{K,R} end
Base.show(io::IO, b::Backend{K,R}) where {K,R} = print(io, backend_name(b))

backend_name(b::Backend{K,R}) where {K,R} = "Backend"

#References can be (single) native references or union or substraction of References
#Unions and subtractions are needed because actual backends frequently fail those operations
abstract type GenericRef{K,T} end

struct EmptyRef{K,T} <: GenericRef{K,T} end
struct UniversalRef{K,T} <: GenericRef{K,T} end

struct NativeRef{K,T} <: GenericRef{K,T}
  value::T
end
struct UnionRef{K,T} <: GenericRef{K,T}
  values::Tuple{Vararg{GenericRef{K,T}}}
end
struct SubtractionRef{K,T} <: GenericRef{K,T}
  value::GenericRef{K,T}
  values::Tuple{Vararg{GenericRef{K,T}}}
end

ensure_ref(b::Backend{K,T}, v::GenericRef{K,T}) where {K,T} = v
ensure_ref(b::Backend{K,T}, v::T) where {K,T} = NativeRef{K,T}(v)
ensure_ref(b::Backend{K,T}, v::Vector{T}) where {K,T} =
  length(v) == 1 ?
    NativeRef{K,T}(v[1]) :
    UnionRef{K,T}(([NativeRef{K,T}(vi) for vi in v]...,))
ensure_ref(b::Backend{K,T}, v::Vector{<:GenericRef{K, T}}) where {K,T} =
  length(v) == 1 ?
    v[1] :
    UnionRef{K,T}((v...,))

# currying
map_ref(b::Backend{K,T}, f::Function) where {K,T} = r -> map_ref(b, f, r)

map_ref(b::Backend{K,T}, f::Function, r::NativeRef{K,T}) where {K,T} = ensure_ref(b, f(r.value))
map_ref(b::Backend{K,T}, f::Function, r::UnionRef{K,T}) where {K,T} = UnionRef{K,T}(map(map_ref(b, f), r.values))
map_ref(b::Backend{K,T}, f::Function, r::SubtractionRef{K,T}) where {K,T} = SubtractionRef{K,T}(map_ref(b, f, r.value), map(map_ref(b, f), r.values))

# currying
collect_ref(b::Backend{K,T}) where {K,T} = r -> collect_ref(b, r)

collect_ref(b::Backend{K,T}, r::EmptyRef{K,T}) where {K,T} = []
collect_ref(b::Backend{K,T}, r::NativeRef{K,T}) where {K,T} = [r.value]
collect_ref(b::Backend{K,T}, r::UnionRef{K,T}) where {K,T} = mapreduce(collect_ref(b), vcat, r.values, init=[])
collect_ref(b::Backend{K,T}, r::SubtractionRef{K,T}) where {K,T} = vcat(collect_ref(b, r.value), mapreduce(collect_ref(b), vcat, r.values, init=[]))

# Boolean algebra laws
unite_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::UniversalRef{K,T}) where {K,T} = r1
unite_ref(b::Backend{K,T}, r0::UniversalRef{K,T}, r1::GenericRef{K,T}) where {K,T} = r0

#To avoid ambiguity
unite_ref(b::Backend{K,T}, r0::UnionRef{K,T}, r1::UnionRef{K,T}) where {K,T} =
  unite_ref(b, unite_refs(b, r0), unite_refs(b, r1))
unite_ref(b::Backend{K,T}, r0::EmptyRef{K,T}, r1::EmptyRef{K,T}) where {K,T} = r0
unite_ref(b::Backend{K,T}, r0::UnionRef{K,T}, r1::EmptyRef{K,T}) where {K,T} = r0
unite_ref(b::Backend{K,T}, r0::EmptyRef{K,T}, r1::UnionRef{K,T}) where {K,T} = r1
unite_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::EmptyRef{K,T}) where {K,T} = r0
unite_ref(b::Backend{K,T}, r0::EmptyRef{K,T}, r1::GenericRef{K,T}) where {K,T} = r1

unite_refs(b::Backend{K,T}, r::UnionRef{K,T}) where {K,T} =
  foldr((r0,r1)->unite_ref(b,r0,r1), r.values, init=EmptyRef{K,T}())
unite_ref(b::Backend{K,T}, r0::UnionRef{K,T}, r1::GenericRef{K,T}) where {K,T} =
  unite_ref(b, unite_refs(b, r0), r1)
unite_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::UnionRef{K,T}) where {K,T} =
  unite_ref(b, r0, unite_refs(b, r1))

intersect_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::UniversalRef{K,T}) where {K,T} = r0
intersect_ref(b::Backend{K,T}, r0::UniversalRef{K,T}, r1::GenericRef{K,T}) where {K,T} = r1
intersect_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::EmptyRef{K,T}) where {K,T} = r1
intersect_ref(b::Backend{K,T}, r0::EmptyRef{K,T}, r1::GenericRef{K,T}) where {K,T} = r0
intersect_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::UnionRef{K,T}) where {K,T} =
  intersect_ref(b, r0, unite_refs(b, r1))
intersect_ref(b::Backend{K,T}, r0::UnionRef{K,T}, r1::GenericRef{K,T}) where {K,T} =
  intersect_ref(b, unite_refs(b, r0), r1)

#To avoid ambiguity
subtract_ref(b::Backend{K,T}, r0::UnionRef{K,T}, r1::UnionRef{K,T}) where {K,T} =
  subtract_ref(b, unite_refs(b, r0), unite_refs(b, r1))
subtract_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::UniversalRef{K,T}) where {K,T} = EmptyRef{K,T}()
subtract_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::EmptyRef{K,T}) where {K,T} = r0
subtract_ref(b::Backend{K,T}, r0::EmptyRef{K,T}, r1::GenericRef{K,T}) where {K,T} = r0
subtract_ref(b::Backend{K,T}, r0::GenericRef{K,T}, r1::UnionRef{K,T}) where {K,T} =
  subtract_ref(b, r0, unite_refs(b, r1))
subtract_ref(b::Backend{K,T}, r0::UnionRef{K,T}, r1::GenericRef{K,T}) where {K,T} =
  subtract_ref(b, unite_refs(b, r0), r1)

# References need to be created, deleted, and recreated, depending on the way the backend works
# For example, each time a shape is consumed, it becomes deleted and might need to be recreated
mutable struct LazyRef{K,R}
  backend::Backend{K,R}
  value::GenericRef{K,R}
  created::Int
  deleted::Int
end

LazyRef(backend::Backend{K,R}) where {K,R} = LazyRef{K,R}(backend, void_ref(backend), 0, 0)
LazyRef(backend::Backend{K,R}, v::GenericRef{K,R}) where {K,R} = LazyRef{K,R}(backend, v, 1, 0)

abstract type Proxy end

backend(s::Proxy) = s.ref.backend
realized(s::Proxy) = s.ref.created == s.ref.deleted + 1
# This is so stupid. We need call-next-method.
really_mark_deleted(s::Proxy) = really_mark_deleted(s.ref)
really_mark_deleted(ref::LazyRef) = ref.deleted = ref.created
really_mark_deleted(s::Any) = nothing
mark_deleted(s::Proxy) = really_mark_deleted(s)
# We also need to propagate this to all dependencies
mark_deleted(ss::Array{<:Proxy}) = foreach(mark_deleted, ss)
mark_deleted(s::Any) = nothing
marked_deleted(s::Proxy) = s.ref.deleted == s.ref.created

ref(s::Proxy) =
  if s.ref.created == s.ref.deleted
    s.ref.value = ensure_ref(s.ref.backend, realize(s.ref.backend, s))
    s.ref.created += 1
    s.ref.value
  elseif s.ref.created == s.ref.deleted + 1
    s.ref.value
  else
    error("Inconsistent creation and deletion")
  end

# We can also use a shape as a surrogate for another shape

ensure_ref(b::Backend{K,T}, v::Proxy) where {K,T} = ref(v)

#This is a dangerous operation. I'm not sure it should exist.
set_ref!(s::Proxy, value) = s.ref.value = value

abstract type Shape <: Proxy end
Base.show(io::IO, s::Shape) =
    print(io, "$(typeof(s))(...)")

Shapes = Vector{<:Shape}

map_ref(f::Function, s::Shape) = map_ref(s.ref.backend, f, ref(s))
collect_ref(s::Shape) = collect_ref(s.ref.backend, ref(s))
collect_ref(ss::Shapes) = mapreduce(collect_ref, vcat, ss, init=[])

#=
Whenever a shape is created, it might be replaced by another shape (e.g., when
it is necessary to avoid collisions), it might be eagerly realized in its backend,
depending on the kind of shape and on the kind of backend (and/or its current state).
Another possibility is for the shape to be saved in some container.
It might also be necessary to record the control flow that caused the shape to be created.
This means that we need to control what happens immediately after a shape is initialized.
The protocol after_init takes care of that.
=#

after_init(a::Any) = a
after_init(s::Shape) = maybe_trace(maybe_collect(maybe_realize(maybe_replace(s))))

#=
Shapes might need to be replaced with other shapes. For example, when elements
join at precise locations, it might be easier, from the programming point of
view, to just duplicate the joining element. However, from the constructive
point of view, just one joining element needs to be produced.
=#

maybe_replace(s::Any) = s

#=
Backends might need to immediately realize a shape while supporting further modifications
e.g., using boolean operations. Others, however, cannot do that and can only realize
shapes by request, presumably, when they have complete information about them.
A middle term might be a backend that supports both modes.
=#

delay_realize(b::Backend, s::Shape) = s
force_realize(b::Backend, s::Shape) = (ref(s); s)

maybe_realize(s::Shape, b::Backend=backend(s)) = maybe_realize(b, s)

#=
Even if a backend is eager, it might be necessary to temporarily delay the
realization of shapes, particularly, when the construction is incremental.
=#

delaying_realize = Parameter(false)
maybe_realize(b::Backend, s::Shape) =
  delaying_realize() ?
    delay_realize(b, s) :
    force_realize(b, s)

abstract type LazyBackend{K,T} <: Backend{K,T} end
maybe_realize(b::LazyBackend, s::Shape) = delay_realize(b, s)
delay_realize(b::LazyBackend, s::Shape) = (push!(b.shapes, s); s)

with_transaction(fn) =
  maybe_realize(with(fn, delaying_realize, true))

#=
Frequently, we need to collect all shapes that are created:
=#

in_shape_collection = Parameter(false)
collected_shapes = Parameter(Shape[])
collect_shape!(s::Shape) = (push!(collected_shapes(), s); s)
collecting_shapes(fn) =
    with(collected_shapes, Shape[]) do
        with(in_shape_collection, true) do
            fn()
        end
        collected_shapes()
    end
maybe_collect(s::Shape) = (in_shape_collection() && collect_shape!(s); s)

######################################################
#Traceability
traceability = Parameter(false)
trace_depth = Parameter(1000)
excluded_modules = Parameter([Base, Base.CoreLogging, Khepri])
# We a dict from shapes to file locations
# and a dict from file locations to shapes
shape_to_file_locations = IdDict()
file_location_to_shapes = Dict()

export traceability, trace_depth, excluded_modules, clear_trace!, shape_source, source_shapes

shape_source(s) = get(shape_to_file_locations, s, [])
source_shapes(file, line) = get(file_location_to_shapes, (file, line), [])

clear_trace!() =
  begin
    empty!(shape_to_file_locations)
    empty!(file_location_to_shapes)
  end
#=
We do not care about frames that are unrelated to the application.
=#
interesting_locations(frames) =
  let locations = [],
      max_depth = min(trace_depth(), length(frames)-0)#14)
    for i in 2:max_depth
      let frame = frames[i],
          linfo = frame.linfo
        if linfo isa Core.CodeInfo ||
           (linfo isa Core.MethodInstance &&
            ! (linfo.def.module in excluded_modules()))
          push!(locations, (frame.file, frame.line))
        end
      end
    end
    locations
  end

trace!(s) =
  let frames = stacktrace(),
      locations = interesting_locations(frames)
    shape_to_file_locations[s] = locations
    for location in locations
      file_location_to_shapes[location] = Shape[get(file_location_to_shapes, location, [])..., s]
    end
    s
  end

maybe_trace(s) = (traceability() && trace!(s); s)

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

replace_in(expr::Expr, replacements) =
    if expr.head == :.
        Expr(expr.head,
             replace_in(expr.args[1], replacements), expr.args[2])
    elseif expr.head == :quote
        expr
    else
        Expr(expr.head,
             map(arg -> replace_in(arg, replacements), expr.args) ...)
    end
replace_in(expr::Symbol, replacements) =
    get(replacements, expr, esc(expr))
replace_in(expr::Any, replacements) =
    expr

showit(s, a) = begin
    print(s)
    print(":")
    println(a)
    a
end

# The undefined backend
struct UndefinedBackend <: Backend{Int,Int} end
connection(b::UndefinedBackend) = throw(UndefinedBackendException())
void_ref(b::UndefinedBackend) = EmptyRef{Int,Int}()
const undefined_backend = UndefinedBackend()

const current_backend = Parameter{Backend}(undefined_backend)
has_current_backend() = current_backend() != undefined_backend

# Side-effect full operations need to have a backend selected and will generate an exception if there is none

struct UndefinedBackendException <: Exception end
Base.show(io::IO, e::UndefinedBackendException) = print(io, "No current backend.")

realize(::UndefinedBackend, ::Shape) = throw(UndefinedBackendException())

# Many functions default the backend to the current_backend and throw an error if there is none.
# We will simplify their definition with a macro:
# @defop delete_all_shapes()
# that expands into
# delete_all_shapes(backend::Backend=current_backend()) = throw(UndefinedBackendException())
# Note that according to Julia semantics the previous definition actually generates two different ones:
# delete_all_shapes() = delete_all_shapes(current_backend())
# delete_all_shapes(backend::Backend) = throw(UndefinedBackendException())
# Hopefully, backends will specialize the function for each specific backend

macro defop(name_params)
    name, params = name_params.args[1], name_params.args[2:end]
    quote
        export $(esc(name))
        $(esc(name))($(map(esc,params)...), backend::Backend=current_backend()) =
            throw(UndefinedBackendException())
    end
end

macro defshapeop(name_params)
    name, params = name_params.args[1], name_params.args[2:end]
    quote
        export $(esc(name))
        $(esc(name))(s::Shape, $(map(esc,params)...), b::Backend=backend(s)) =
            throw(UndefinedBackendException())
    end
end

backend(backend::Backend) = switch_to_backend(current_backend(), backend)
switch_to_backend(from::Backend, to::Backend) = current_backend(to)

@defop current_backend_name()
@defop delete_all_shapes()
@defop set_length_unit(unit::String)
@defop reset_backend()
@defop save_as(pathname::String, format::String)

struct WrongTypeForParam <: Exception
  param::Symbol
  value::Any
  expected_type::Type
end
Base.showerror(io::IO, e::WrongTypeForParam) =
  print(io, "$(e.param) expected a $(e.expected_type) but got $(e.value) of type $(typeof(e.value))")

macro defproxy(name, parent, fields...)
  name_str = string(name)
  struct_name = esc(Symbol(string(map(uppercasefirst,split(name_str,'_'))...)))
  field_names = map(field -> field.args[1].args[1], fields)
  field_types = map(field -> esc(field.args[1].args[2]), fields)
  field_inits = map(field -> field.args[2], fields)
#  field_renames = map(esc ∘ Symbol ∘ uppercasefirst ∘ string, field_names)
  field_renames = map(Symbol ∘ string, field_names)
  field_replacements = Dict(zip(field_names, field_renames))
  struct_fields = map((name,typ) -> :($(name) :: $(typ)), field_names, field_types)
#  opt_params = map((name,typ,init) -> :($(name) :: $(typ) = $(init)), field_renames, field_types, field_inits)
#  key_params = map((name,typ,rename) -> :($(name) :: $(typ) = $(rename)), field_names, field_types, field_renames)
#  mk_param(name,typ) = Expr(:kw, Expr(:(::), name, typ))
  mk_param(name,typ,init) = Expr(:kw, name, init) #Expr(:kw, Expr(:(::), name, typ), init)
  opt_params = map(mk_param, field_renames, field_types, map(init -> replace_in(init, field_replacements), field_inits))
  key_params = map(mk_param, field_names, field_types, field_renames)
  constructor_name = esc(name)
  predicate_name = esc(Symbol("is_", name_str))
  #mk_convert(name,typ) = :(isa($(esc(name)), $(typ)) ? $(esc(name)) : throw(WrongTypeForParam($(QuoteNode(name)), $(esc(name)), $(typ))))
  mk_convert(name,typ) = :($(esc(name)))
  field_converts = map(mk_convert, field_names, field_types)
  selector_names = map(field_name -> esc(Symbol(name_str, "_", string(field_name))), field_names)
  quote
    export $(constructor_name), $(struct_name), $(predicate_name) #, $(selector_names...)
    struct $struct_name <: $parent
      ref::LazyRef
      $(struct_fields...)
    end
    # we don't need to convert anything because Julia already does that with the default constructor
    # and, by the same idea, we don't need to define parameter types.
    @noinline $(constructor_name)($(opt_params...); $(key_params...), backend::Backend=current_backend(), ref::LazyRef=LazyRef(backend)) =
      after_init($(struct_name)(ref, $(field_converts...)))
    $(predicate_name)(v::$(struct_name)) = true
    $(predicate_name)(v::Any) = false
    $(map((selector_name, field_name) -> :($(selector_name)(v::$(struct_name)) = v.$(field_name)),
          selector_names, field_names)...)
    Khepri.mark_deleted(v::$(struct_name)) =
      if ! marked_deleted(v)
        really_mark_deleted(v)
        $(map(field_name -> :(mark_deleted(v.$(field_name))), field_names)...)
      end
    Khepri.meta_program(v::$(struct_name)) =
        Expr(:call, $(Expr(:quote, name)), $(map(field_name -> :(meta_program(v.$(field_name))), field_names)...))
  end
end

abstract type Shape0D <: Shape end
abstract type Shape1D <: Shape end
abstract type Shape2D <: Shape end
abstract type Shape3D <: Shape end

is_curve(s::Shape) = false
is_surface(s::Shape) = false
is_solid(s::Shape) = false

is_curve(s::Shape1D) = true
is_surface(s::Shape2D) = true
is_solid(s::Shape3D) = true

# HACK: Fix element type
Shapes0D = Vector{<:Any}
Shapes1D = Vector{<:Any}
Shapes2D = Vector{<:Any}

@defproxy(empty_shape, Shape0D)
@defproxy(universal_shape, Shape3D)
@defproxy(point, Shape0D, position::Loc=u0())
@defproxy(line, Shape1D, vertices::Locs=[u0(), ux()])
line(v0::Loc, v1::Loc, vs...) = line([v0, v1, vs...])
@defproxy(closed_line, Shape1D, vertices::Locs=[u0(), ux(), uy()])
closed_line(v0::Loc, v1::Loc, vs...) = closed_line([v0, v1, vs...])
@defproxy(spline, Shape1D, points::Locs=[u0(), ux(), uy()], v0::Union{Bool,Vec}=false, v1::Union{Bool,Vec}=false,
          interpolator::Parameter{Any}=Parameter{Any}(missing))
spline(v0::Loc, v1::Loc, vs...) = spline([v0, v1, vs...])

#=
evaluate(s::Spline, t::Real) =
  let interpolator = s.interpolator
    if ismissing(interpolator())
      interpolator(curve_interpolator(s.points))
    end
    let p = interpolator()(t),
        vt = Interpolations.gradient(interpolator(), t)[1],
        vn = Interpolations.hessian(interpolator(), t)[1]
      loc_from_o_vx_vy(
        xyz(p[1], p[2], p[3], world_cs),
        vxyz(vt[1], vt[2], vt[3], world_cs),
        vxyz(vn[1], vn[2], vn[3], world_cs))
    end
  end

=#
evaluate(s::Spline, t::Real) =
  let interpolator = s.interpolator
    if ismissing(interpolator())
      interpolator(curve_interpolator(s.points))
    end
    let p = interpolator()(t),
        vt = Interpolations.gradient(interpolator(), t)[1],
        vn = Interpolations.hessian(interpolator(), t)[1],
        vy = cross(vt, vn)
      loc_from_o_vx_vy(
        xyz(p[1], p[2], p[3], world_cs),
        vxyz(vn[1], vn[2], vn[3], world_cs),
        vxyz(vy[1], vy[2], vy[3], world_cs))
    end
  end

curve_domain(s::Spline) = (0.0, 1.0)
frame_at(s::Spline, t::Real) = evaluate(s, t)
map_division(f::Function, s::Spline, n::Int, backend::Backend=backend(s)) =
  backend_map_division(backend, f, s, n)
#=HACK, THIS IS NOT READY, YET. COMPARE WITH THE BACKEND VERSION!!!!!!
  let (t1, t2) = curve_domain(s)
    map_division(t1, t2, n) do t
        f(frame_at(s, t))
    end
  end
=#
#(def-base-shape 1D-shape (spline* [pts : (Listof Loc) (list (u0) (ux) (uy))] [v0 : (U Boolean Vec) #f] [v1 : (U Boolean Vec) #f]))

@defproxy(closed_spline, Shape1D, points::Locs=[u0(), ux(), uy()])
closed_spline(v0, v1, vs...) = closed_spline([v0, v1, vs...])
@defproxy(circle, Shape1D, center::Loc=u0(), radius::Real=1)
@defproxy(arc, Shape1D, center::Loc=u0(), radius::Real=1, start_angle::Real=0, amplitude::Real=pi)
@defproxy(elliptic_arc, Shape1D, center::Loc=u0(), radius_x::Real=1, radius_y::Real=1, start_angle::Real=0, amplitude::Real=pi)
@defproxy(ellipse, Shape1D, center::Loc=u0(), radius_x::Real=1, radius_y::Real=1)
@defproxy(polygon, Shape1D, vertices::Locs=[u0(), ux(), uy()])
polygon(v0, v1, vs...) = polygon([v0, v1, vs...])
@defproxy(regular_polygon, Shape1D, edges::Integer=3, center::Loc=u0(), radius::Real=1, angle::Real=0, inscribed::Bool=true)
@defproxy(rectangle, Shape1D, corner::Loc=u0(), dx::Real=1, dy::Real=1)
rectangle(p::Loc, q::Loc) =
  let v = in_cs(q - p, p.cs)
    rectangle(p, v.x, v.y)
  end
@defproxy(surface_circle, Shape2D, center::Loc=u0(), radius::Real=1)
@defproxy(surface_arc, Shape2D, center::Loc=u0(), radius::Real=1, start_angle::Real=0, amplitude::Real=pi)
@defproxy(surface_elliptic_arc, Shape2D, center::Loc=u0(), radius_x::Real=1, radius_y::Real=1, start_angle::Real=0, amplitude::Real=pi)
@defproxy(surface_ellipse, Shape2D, center::Loc=u0(), radius_x::Real=1, radius_y::Real=1)
@defproxy(surface_polygon, Shape2D, vertices::Locs=[u0(), ux(), uy()])
surface_polygon(v0, v1, vs...) = surface_polygon([v0, v1, vs...])
@defproxy(surface_regular_polygon, Shape2D, edges::Integer=3, center::Loc=u0(), radius::Real=1, angle::Real=0, inscribed::Bool=true)
@defproxy(surface_rectangle, Shape2D, corner::Loc=u0(), dx::Real=1, dy::Real=1)
surface_rectangle(p::Loc, q::Loc) =
  let v = in_cs(q - p, p.cs)
    surface_rectangle(p, v.x, v.y)
  end
@defproxy(surface, Shape2D, frontier::Shapes1D=[circle()])
surface(c0::Shape, cs...) = surface([c0, cs...])
#To be removed
surface_from = surface

@defproxy(surface_path, Shape2D, path::ClosedPath=[circular_path()])
realize(b::Backend, s::SurfacePath) = backend_fill(b, s.path)

surface_boundary(s::Shape2D, backend::Backend=current_backend()) =
    backend_surface_boundary(backend, s)

curve_domain(s::Shape1D, backend::Backend=current_backend()) =
    backend_curve_domain(backend, s)
map_division(f::Function, s::Shape1D, n::Int, backend::Backend=current_backend()) =
    backend_map_division(backend, f, s, n)

surface_domain(s::Shape2D, backend::Backend=current_backend()) =
    backend_surface_domain(backend, s)
map_division(f::Function, s::Shape2D, nu::Int, nv::Int, backend::Backend=current_backend()) =
    backend_map_division(backend, f, s, nu, nv)

@defproxy(text, Shape0D, str::String="", corner::Loc=u0(), height::Real=1)

export text_centered
text_centered(str::String="", center::Loc=u0(), height::Real=1) =
  text(str, add_xy(center, -length(str)*height*0.85/2, -height/2), height)

# This is for unknown shapes (they are opaque, the only thing you can do with then
# might be just delete them)
@defproxy(unknown, Shape3D, baseref::Any=required())

@defproxy(sphere, Shape3D, center::Loc=u0(), radius::Real=1)
@defproxy(torus, Shape3D, center::Loc=u0(), re::Real=1, ri::Real=1/2)
@defproxy(cuboid, Shape3D,
  b0::Loc=u0(),        b1::Loc=add_x(b0,1), b2::Loc=add_y(b1,1), b3::Loc=add_x(b2,-1),
  t0::Loc=add_z(b0,1), t1::Loc=add_x(t0,1), t2::Loc=add_y(t1,1), t3::Loc=add_x(t2,-1))

@defproxy(regular_pyramid_frustum, Shape3D, edges::Integer=4, cb::Loc=u0(), rb::Real=1, angle::Real=0, h::Real=1, rt::Real=1, inscribed::Bool=false)
regular_pyramid_frustum(edges::Integer, cb::Loc, rb::Real, angle::Real, ct::Loc, rt::Real=1, inscribed::Bool=false) =
  let (c, h) = position_and_height(cb, ct)
    regular_pyramid_frustum(edges, c, rb, angle, h, rt, inscribed)
  end

@defproxy(regular_pyramid, Shape3D, edges::Integer=3, cb::Loc=u0(), rb::Real=1, angle::Real=0, h::Real=1, inscribed::Bool=false)
regular_pyramid(edges::Integer, cb::Loc, rb::Real, angle::Real, ct::Loc, inscribed::Bool=false) =
  let (c, h) = position_and_height(cb, ct)
    regular_pyramid(edges, c, rb, angle, h, inscribed)
  end

@defproxy(irregular_pyramid_frustum, Shape3D, bs::Locs=[ux(), uy(), uxy()], ts::Locs=[uxz(), uyz(), uxyz()])
@defproxy(irregular_pyramid, Shape3D, bs::Locs=[ux(), uy(), uxy()], t::Loc=uz())

@defproxy(regular_prism, Shape3D, edges::Integer=3, cb::Loc=u0(), r::Real=1, angle::Real=0, h::Real=1, inscribed::Bool=false)
regular_prism(edges::Integer, cb::Loc, r::Real, angle::Real, ct::Loc, inscribed::Bool=false) =
  let (c, h) = position_and_height(cb, ct)
    regular_prism(edges, c, r, angle, h, inscribed)
  end
@defproxy(irregular_prism, Shape3D, bs::Locs=[ux(), uy(), uxy()], v::Vec=vz(1))
irregular_prism(bs::Locs, h::Real) =
  irregular_prism(bs, vz(h))

@defproxy(right_cuboid, Shape3D, cb::Loc=u0(), width::Real=1, height::Real=1, h::Real=1, angle::Real=0)
right_cuboid(cb::Loc, width::Real, height::Real, ct::Loc, angle::Real=0; backend::Backend=current_backend()) =
  let (c, h) = position_and_height(cb, ct)
    right_cuboid(c, width, height, h, angle, backend=backend)
  end
@defproxy(box, Shape3D, c::Loc=u0(), dx::Real=1, dy::Real=dx, dz::Real=dy)
box(c0::Loc, c1::Loc) =
  let v = in_cs(c1, c0)-c0
    box(c0, v.x, v.y, v.z)
  end
@defproxy(cone, Shape3D, cb::Loc=u0(), r::Real=1, h::Real=1)
cone(cb::Loc, r::Real, ct::Loc) =
  let (c, h) = position_and_height(cb, ct)
    cone(c, r, h)
  end
@defproxy(cone_frustum, Shape3D, cb::Loc=u0(), rb::Real=1, h::Real=1, rt::Real=1)
cone_frustum(cb::Loc, rb::Real, ct::Loc, rt::Real) =
  let (c, h) = position_and_height(cb, ct)
    cone_frustum(c, rb, h, rt)
  end
@defproxy(cylinder, Shape3D, cb::Loc=u0(), r::Real=1, h::Real=1)
cylinder(cb::Loc, r::Real, ct::Loc) =
  let (c, h) = position_and_height(cb, ct)
    cylinder(c, r, h)
  end

@defproxy(extrusion, Shape3D, profile::Shape=point(), v::Vec=vz(1))
extrusion(profile, h::Real) =
  extrusion(profile, vz(h))

realize(b::Backend, s::Extrusion) =
  backend_extrusion(backend(s), s.profile, s.v)

backend_extrusion(b::Backend, p::Point, v::Vec) =
  realize_and_delete_shapes(line([p.position, p.position + v], backend=b), [p])

@defproxy(sweep, Shape3D, path::Shape1D=circle(), profile::Shape=point(), rotation::Real=0, scale::Real=1)

realize(b::Backend, s::Sweep) =
  backend_sweep(backend(s), s.path, s.profile, s.rotation, s.scale)

@defproxy(revolve_point, Shape1D, profile::Shape0D=point(), p::Loc=u0(), n::Vec=vz(1,p.cs), start_angle::Real=0, amplitude::Real=2*pi)
@defproxy(revolve_curve, Shape2D, profile::Shape1D=line(), p::Loc=u0(), n::Vec=vz(1,p.cs), start_angle::Real=0, amplitude::Real=2*pi)
@defproxy(revolve_surface, Shape3D, profile::Shape2D=circle(), p::Loc=u0(), n::Vec=vz(1,p.cs), start_angle::Real=0, amplitude::Real=2*pi)
revolve(profile::Shape=point(x(1)), p::Loc=u0(), n::Vec=vz(1,p.cs), start_angle::Real=0, amplitude::Real=2*pi) =
  if is_point(profile)
    revolve_point(profile, p, n, start_angle, amplitude)
  elseif is_curve(profile)
    revolve_curve(profile, p, n, start_angle, amplitude)
  elseif is_surface(profile)
    revolve_surface(profile, p, n, start_angle, amplitude)
  elseif is_union_shape(profile)
    union(map(s->revolve(s, p, n, start_angle, amplitude), profile.shapes))
  elseif is_empty_shape(profile)
    profile
  else
    error("Profile is neither a point nor a curve nor a surface")
  end

realize(b::Backend, s::RevolvePoint) =
  backend_revolve_point(b, s.profile, s.p, s.n, s.start_angle, s.amplitude)
realize(b::Backend, s::RevolveCurve) =
  backend_revolve_curve(b, s.profile, s.p, s.n, s.start_angle, s.amplitude)
realize(b::Backend, s::RevolveSurface) =
  backend_revolve_surface(b, s.profile, s.p, s.n, s.start_angle, s.amplitude)

@defproxy(loft_points, Shape1D, profiles::Shapes0D=Shape[], rails::Shapes=Shape[], ruled::Bool=false, closed::Bool=false)
@defproxy(loft_curves, Shape2D, profiles::Shapes1D=Shape[], rails::Shapes=Shape[], ruled::Bool=false, closed::Bool=false)
@defproxy(loft_surfaces, Shape3D, profiles::Shapes2D=Shape[], rails::Shapes=Shape[], ruled::Bool=false, closed::Bool=false)
@defproxy(loft_curve_point, Shape2D, profile::Shape1D=circle(), point::Shape0D=point(z(1)))
@defproxy(loft_surface_point, Shape3D, profile::Shape2D=surface_circle(), point::Shapes=point(z(1)))

loft(profiles::Shapes=Shape[], rails::Shapes=Shape[], ruled::Bool=false, closed::Bool=false) =
  if all(is_point, profiles)
    loft_points(profiles, rails, ruled, closed)
  elseif all(is_curve, profiles)
    loft_curves(profiles, rails, ruled, closed)
  elseif all(is_surface, profiles)
    loft_surfaces(profiles, rails, ruled, closed)
  elseif length(profiles) == 2
    let (p, sh) = if is_point(profiles[1])
                    (profiles[1], profiles[2])
                  elseif is_point(profiles[2])
                    (profiles[2], profiles[1])
                  else
                    error("Cross sections are neither points nor curves nor surfaces")
                  end
      if is_curve(sh)
        loft_curve_point(sh, p)
      elseif is_surface(sh)
        loft_surface_point(sh, p)
      else
        error("Can't loft the shapes")
      end
    end
  else
    error("Cross sections are neither points nor curves nor surfaces")
  end

loft_ruled(profiles::Shapes=Shape[]) = loft(profiles, Shape[], true, false)
export loft, loft_ruled

realize(b::Backend, s::LoftPoints) = backend_loft_points(backend(s), s.profiles, s.rails, s.ruled, s.closed)
realize(b::Backend, s::LoftCurves) = backend_loft_curves(backend(s), s.profiles, s.rails, s.ruled, s.closed)
realize(b::Backend, s::LoftSurfaces) = backend_loft_surfaces(backend(s), s.profiles, s.rails, s.ruled, s.closed)
realize(b::Backend, s::LoftCurvePoint) = backend_loft_curve_point(backend(s), s.profile, s.point)
realize(b::Backend, s::LoftSurfacePoint) = backend_loft_surface_point(backend(s), s.profile, s.point)

backend_loft_points(b::Backend, profiles::Shapes, rails::Shapes, ruled::Bool, closed::Bool) =
  let f = (ruled ? (closed ? polygon : line) : (closed ? closed_spline : spline))
    and_delete_shapes(ref(f(map(point_position, profiles), backend=b)),
                      vcat(profiles, rails))
  end

@defproxy(move, Shape3D, shape::Shape=point(), v::Vec=vx())
@defproxy(scale, Shape3D, shape::Shape=point(), s::Real=1, p::Loc=u0())
@defproxy(rotate, Shape3D, shape::Shape=point(), angle::Real=0, p::Loc=u0(), v::Vec=vz(1,p.cs))
@defproxy(transform, Shape3D, shape::Shape=point(), xform::Loc=u0())

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

# We can also translate some shapes
translate(s::Line, v::Vec) = line(map(p -> p+v, s.vertices))
translate(s::Polygon, v::Vec) = polygon(map(p -> p+v, s.vertices))
translate(s::Circle, v::Vec) = circle(s.center+v, s.radius)
translate(s::Text, v::Vec) = text(s.str, s.c+v, s.h)

# We can translate arrays of Shapes
translate(ss::Shapes, v::Vec) = translate.(ss, v)

# We can compute the length of shapes as long as we can convert them
curve_length(s::Shape) = curve_length(convert(Path, s))

# We will also need to compute a bounding rectangle
bounding_rectangle(s::Union{Line, Polygon}) =
    bounding_rectangle(s.vertices)

bounding_rectangle(pts::Locs) =
    let min_p = pts[1]
        max_p = min_p
        for i in 2:length(pts)
            min_p = min_loc(min_p, pts[i])
            max_p = max_loc(max_p, pts[i])
        end
        [min_p, max_p]
    end

bounding_rectangle(ss::Shapes) =
    bounding_rectangle(mapreduce(bounding_rectangle, vcat, ss))

#####################################################################
## We might also be insterested in seeing paths
# This is more for debuging

stroke(path::Path, backend::Backend=current_backend()) = backend_stroke(backend, path)
# We also need a colored stroke (and probably, something that changes line thickness)
stroke(path::Path, color::RGB, backend::Backend=current_backend()) = backend_stroke_color(backend, path, color)
# By default, we ignore the color
backend_stroke_color(backend::Backend, path::Path, color::RGB) = backend_stroke(backend, path)

# The default implementation for stroking segmented path in the backend relies on two
# dedicated functions backend_stroke_arc and backend_stroke_line
stroke(path::PathOps, backend::Backend=current_backend()) = backend_stroke(backend, path)

backend_stroke(b::Backend, path::PathOps) =
    begin
        start, curr, refs = path.start, path.start, []
        for op in path.ops
            start, curr, refs = backend_stroke_op(b, op, start, curr, refs)
        end
        if path.closed
            push!(refs, backend_stroke_line(b, [curr, start]))
        end
        backend_stroke_unite(b, refs)
    end
#=
backend_stroke_op(b::Backend, op::MoveToOp, start::Loc, curr::Loc, refs) =
    (op.loc, op.loc, refs)
backend_stroke_op(b::Backend, op::MoveOp, start::Loc, curr::Loc, refs) =
    (start, curr + op.vec, refs)
backend_stroke_op(b::Backend, op::LineToOp, start::Loc, curr::Loc, refs) =
    (start, op.loc, push!(refs, backend_stroke_line(b, [curr, op.loc])))
=#
backend_stroke_op(b::Backend, op::LineOp, start::Loc, curr::Loc, refs) =
    (start, curr + op.vec, push!(refs, backend_stroke_line(b, [curr, curr + op.vec])))
#backend_stroke_op(b::Backend, op::CloseOp, start::Loc, curr::Loc, refs) =
#    (start, start, push!(refs, backend_stroke_line(b, [curr, start])))
backend_stroke_op(b::Backend, op::ArcOp, start::Loc, curr::Loc, refs) =
    let center = curr - vpol(op.radius, op.start_angle)
        (start,
         center + vpol(op.radius, op.start_angle + op.amplitude),
         push!(refs, backend_stroke_arc(b, center, op.radius, op.start_angle, op.amplitude)))
     end
#=
backend_stroke_op(b::Backend, op::LineXThenYOp, start::Loc, curr::Loc, refs) =
    (start,
     start + op.vec,
     push!(refs, backend_stroke_line(b, [curr, curr + vec_in(op.vec, curr.cs).x, curr + op.vec])))

backend_stroke_op(b::Backend, op::LineYThenXOp, start::Loc, curr::Loc, refs) =
    (start,
     start + op.vec,
     push!(refs, backend_stroke_line(b, [curr, curr + vec_in(op.vec, curr.cs).y, curr + op.vec])))
backend_stroke_op(b::Backend, op::LineToXThenToYOp, start::Loc, curr::Loc, refs) =
    (start, op.loc, push!(refs, backend_stroke_line(b, [curr, xy(curr.x, loc_in(op.loc, curr.cs).x, curr.cs), op.loc])))
backend_stroke_op(b::Backend, op::LineToYThenToXOp, start::Loc, curr::Loc, refs) =
    (start, op.loc, push!(refs, backend_stroke_line(b, [curr, xy(curr.x, loc_in(op.loc, curr.cs).y, curr.cs), op.loc])))
=#

backend_stroke(b::Backend, path::PathSet) =
    for p in path.paths
        backend_stroke(b, p)
    end

# The default implementation for filling segmented path in the backend relies on
# a dedicated function backend_fill_curves

import Base.fill
fill(path, backend=current_backend()) = backend_fill(backend, path)
backend_fill(b, path) = backend_fill_curves(b, backend_stroke(b, path))

backend_stroke(b::Backend, path::Union{OpenPathSequence,ClosedPathSequence}) =
    backend_stroke_unite(b, map(path->backend_stroke(b, path), path.paths))
backend_fill(b::Backend, path::ClosedPathSequence) =
    backend_fill_curves(b, map(path->backend_stroke(b, path), path.paths))

#####################################################################
## Conversions

convert(::Type{ClosedPath}, s::Rectangle) =
  closed_path([MoveToOp(s.corner), RectOp(vxy(s.dx, s.dy))])
convert(::Type{ClosedPath}, s::Circle) =
  closed_path([CircleOp(s.center, s.radius)])
convert(::Type{Path}, s::Line) = convert(OpenPath, s.vertices)

convert(::Type{ClosedPath}, s::Polygon) =
    let vs = polygon_vertices(s)
        closed_path(vcat([MoveToOp(vs[1])], map(LineToOp, vs[2:end]), [CloseOp()]))
    end
#####################################################################

export curve_domain, surface_domain, frame_at
surface_domain(s::SurfaceRectangle) = (0, s.dx, 0, s.dy)
surface_domain(s::SurfaceCircle) = (0, s.radius, 0, 2pi)
surface_domain(s::SurfaceArc) = (0, s.radius, s.start_angle, s.amplitude)

frame_at(c::Shape1D, t::Real) = backend_frame_at(backend(c), c, t)
frame_at(s::Shape2D, u::Real, v::Real) = backend_frame_at(backend(s), s, u, v)

#Some specific cases can be handled in an uniform way without the backend
frame_at(s::SurfaceRectangle, u::Real, v::Real) = add_xy(s.corner, u, v)
frame_at(s::SurfaceCircle, u::Real, v::Real) = add_pol(s.center, u, v)

import Base.union
export union, intersection, subtraction

@defproxy(union_shape, Shape3D, shapes::Shapes=Shape[])
union(shapes::Shapes) =
  let non_empty_shapes = filter(s -> !is_empty_shape(s), shapes)
    non_empty_shapes == [] ?
      empty_shape() :
      union_shape(non_empty_shapes)
  end

union(shape::Shape, shapes...) = union([shape, shapes...])

@defproxy(intersection_shape, Shape3D, shapes::Shapes=Shape[])
intersection(shapes::Shapes) = intersection_shape(shapes)
intersection(shape::Shape, shapes...) = intersection_shape([shape, shapes...])

@defproxy(subtraction_shape2D, Shape2D, shape::Shape=surface_circle(), shapes::Shapes=Shape[])
@defproxy(subtraction_shape3D, Shape3D, shape::Shape=surface_sphere(), shapes::Shapes=Shape[])
subtraction(shape::Shape2D, shapes...) = subtraction_shape2D(shape, [shapes...])
subtraction(shape::Shape3D, shapes...) = subtraction_shape3D(shape, [shapes...])

@defproxy(slice, Shape3D, shape::Shape=sphere(), p::Loc=u0(), n::Vec=vz(1))

@defproxy(mirror, Shape3D, shape::Shape=sphere(), p::Loc=u0(), n::Vec=vz(1))
@defproxy(union_mirror, Shape3D, shape::Shape=sphere(), p::Loc=u0(), n::Vec=vz(1))

@defproxy(surface_grid, Shape2D, points::AbstractMatrix{<:Loc}=zeros(Loc,(2,2)), closed_u::Bool=false, closed_v::Bool=false,
          interpolator::Parameter{Any}=Parameter{Any}(missing))
surface_interpolator(pts::AbstractMatrix{<:Loc}) =
    let pts = map(pts) do p
                let v = in_world(p).raw
                  SVector{3,Float64}(v[1], v[2], v[3])
                end
              end
        Interpolations.scale(
            interpolate(pts, BSpline(Cubic(Natural(OnGrid())))),
            range(0,stop=1,length=size(pts, 1)),
            range(0,stop=1,length=size(pts, 2)))
    end

# For interpolator to work, we need this:

convert(::Type{AbstractMatrix{<:Loc}}, pts::Vector{<:Vector{<:Loc}}) =
  permutedims(hcat(pts...))

evaluate(s::SurfaceGrid, u::Real, v::Real) =
  let interpolator = s.interpolator
    if ismissing(interpolator())
      interpolator(surface_interpolator(s.points))
    end
    let p = interpolator()(u,v)
        v = Interpolations.gradient(interpolator(), u, v)
      loc_from_o_vx_vy(
        xyz(p[1], p[2], p[3], world_cs),
        vxyz(v[1][1], v[1][2], v[1][3], world_cs),
        vxyz(v[2][1], v[2][2], v[2][3], world_cs))
    end
  end

surface_domain(s::SurfaceGrid) = (0.0, 1.0, 0.0, 1.0)
frame_at(s::SurfaceGrid, u::Real, v::Real) = evaluate(s, u, v)
map_division(f::Function, s::SurfaceGrid, nu::Int, nv::Int, backend::Backend=current_backend()) =
  let (u1, u2, v1, v2) = surface_domain(s)
    map_division(u1, u2, nu) do u
      map_division(v1, v2, nv) do v
        f(frame_at(s, u, v))
      end
    end
  end

@defproxy(thicken, Shape3D, shape::Shape=surface_circle(), thickness::Real=1)

# Blocks

@defproxy(block, Shape, name::String="Block", shapes::Shapes = Shape[])
@defproxy(block_instance, Shape, block::Block=required(), loc::Loc=u0(), scale::Real=1.0)

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

#Backends might use different communication mechanisms, e.g., sockets, COM, RMI, etc

#We start with socket-based communication
struct SocketBackend{K,T} <: Backend{K,T}
  connection::LazyParameter{TCPSocket}
  remote::NamedTuple
end

# To simplify remote calls
macro remote(b, call)
  let op = call.args[1],
      args = map(esc, call.args[2:end]),
      b = esc(b)
    :(call_remote(getfield(getfield($(b), :remote), $(QuoteNode(op))), $(b).connection(), $(args...)))
  end
end

macro get_remote(b, op)
  let b = esc(b)
    :(getfield(getfield($(b), :remote), $(QuoteNode(op))))
  end
end

SocketBackend{K,T}(c::LazyParameter{TCPSocket}) where {K,T} =
  SocketBackend{K,T}(c, NamedTuple{}())

connection(b::SocketBackend{K,T}) where {K,T} = b.connection()

reset_backend(b::SocketBackend) =
  begin
    for f in b.remote
      reset_opcode(f)
    end
    close(b.connection())
    reset(b.connection)
  end

#One less dynamic option is to use a file-based backend. To that end, we implement
#the IOBuffer_Backend

struct IOBufferBackend{K,T} <: Backend{K,T}
  out::IOBuffer
end
connection(backend::IOBufferBackend) = backend.out

################################################################################
bounding_box(shape::Shape) =
  bounding_box([shape])

bounding_box(shapes::Shapes=Shape[]) =
  if isempty(shapes)
    [u0(), u0()]
  else
    backend_bounding_box(backend(shapes[1]), shapes)
  end

delete_shape(shape::Shape) =
  delete_shapes([shape])

delete_shapes(shapes::Shapes=Shape[]) =
  if ! isempty(shapes)
    to_delete = filter(realized, shapes)
    backend_delete_shapes(backend(shapes[1]), to_delete)
    foreach(mark_deleted, shapes)
  end

and_delete_shape(r::Any, shape::Shape) =
  begin
    delete_shape(shape)
    r
  end

and_delete_shapes(r::Any, shapes::Shapes) =
  begin
    delete_shapes(shapes)
    r
  end

and_mark_deleted(r::Any, shape) =
  begin
    mark_deleted(shape)
    r
  end

realize_and_delete_shapes(shape::Shape, shapes::Shapes) =
    and_delete_shapes(ref(shape), shapes)

# Common implementations for realize function

realize(b::Backend, s::UnionShape) =
    unite_refs(b, map(ref, s.shapes))

realize(b::Backend, s::SubtractionShape2D) =
    subtract_ref(b, ref(s.shape), unite_refs(b, map(ref, s.shapes)))
realize(b::Backend, s::SubtractionShape3D) =
    subtract_ref(b, ref(s.shape), unite_refs(b, map(ref, s.shapes)))

function startSketchup(port)
  ENV["ROSETTAPORT"] = port
  args = "C:\\Users\\aml\\Dropbox\\AML\\Projects\\rosetta\\sketchup\\rosetta.rb"
  println(args)
  run(`cmd /C Sketchup -RubyStartup $args`)
  #Start listening for Sketchup
  listener = listen(port)
  connection = listener.accept()
  readline(connection) == "connected" ? connection : error("Could not connect!")
end

# CAD
@defop all_shapes()
@defop all_shapes_in_layer(layer)
@defop delete_all_shapes_in_layer(layer)

@defop disable_update()
@defop enable_update()
@defop set_view(camera::Loc, target::Loc, lens::Real)
@defop get_view()
@defop zoom_extents()
@defop view_top()
@defop get_layer(name::String)
@defop create_layer(name::String)
@defop create_layer(name::String, color::RGB)
@defop current_layer()
@defop current_layer(layer)
@defop set_layer_active(layer, status)
@defop switch_to_layer(layer)
@defop get_material(name::String)
@defop create_material(name::String)
@defop current_material()
@defop current_material(material)

angle_of_view(size, focal_length) = 2atan(size/2focal_length)

function dolly_effect(camera, target, lens, new_camera)
  cur_dist = distance(camera, target)
  new_dist = distance(new_camera, target)
  new_lens = lens*new_dist/cur_dist
  view(new_camera, target, new_lens)
end

dolly_effect_pull_back(delta) = begin
  camera, target, lens = get_view()
  d = distance(camera, target)
  new_camera = target + (camera-target)*(d+delta)/d
  dolly_effect(camera, target, lens, new_camera)
end

@defop select_position(prompt::String="Select a position")
@defop select_positions(prompt::String="Select positions")
@defop select_point(prompt::String="Select a point")
@defop select_points(prompt::String="Select points")
@defop select_curve(prompt::String="Select a curve")
@defop select_curves(prompt::String="Select curves")
@defop select_surface(prompt::String="Select a surface")
@defop select_surfaces(prompt::String="Select surfaces")
@defop select_solid(prompt::String="Select a solid")
@defop select_solids(prompt::String="Select solids")
@defop select_shape(prompt::String="Select a shape")
@defop select_shapes(prompt::String="Select shapes")
@defop highlight_shapes(shapes::Shapes)
@defshapeop highlight_shape()
@defshapeop register_for_changes()
@defshapeop unregister_for_changes()
@defshapeop waiting_for_changes()
@defop changed_shape(shapes::Shapes)

capture_shape(s=select_shape("Select shape to be captured")) =
  if s != nothing
    generate_captured_shape(s, backend(s))
  end

capture_shapes(ss=select_shapes("Select shapes to be captured")) =
    generate_captured_shapes(ss, backend(ss[1]))

register_for_changes(shapes::Shapes) =
  map(shapes) do shape
    register_for_changes(shape, backend(shape))
  end

unregister_for_changes(shapes::Shapes) =
  map(shapes) do shape
    unregister_for_changes(shape, backend(shape))
  end

waiting_for_changes(shapes::Shapes) =
  waiting_for_changes(shapes[1], backend(shapes[1]))

export on_change
on_change(f, shape::Shape) = on_change(f, [shape])
on_change(f, shapes) =
  let registered = register_for_changes(shapes)
    try
      while waiting_for_changes(shapes)
        let changed = changed_shape(shapes)
          f()
        end
      end
    finally
      unregister_for_changes(registered)
    end
  end

#
export with_shape_dependency
with_shape_dependency(f, ss) =
    let shapes = collecting_shapes() do
                    f()
                end
        on_change(ss) do
            try
                delete_shapes(shapes)
            catch e
            end
            shapes = collecting_shapes() do
                f()
            end
        end
    end

#
export internalize_shape, internalize_shapes

internalize_shape(s=select_shape("Select shape to be internalized")) =
  if s != nothing
    println(meta_program(s))
  end

internalize_shapes(ss=select_shapes("Select shapes to be internalized")) =
    println(meta_program(ss))

# Later, this will be used to create images.
export to_render
to_render(f, name) =
  begin
    delete_all_shapes()
    f()
  end

# Seletion

select_one_with_prompt(prompt::String, b::Backend, f::Function) =
  let ans = select_many_with_prompt(prompt, b, f)
    length(ans) > 0 ? ans[1] : nothing
  end

select_many_with_prompt(prompt::String, b::Backend, f::Function) =
  begin
    @info "$(prompt) on the $(b) backend."
    map(id -> shape_from_ref(id, b), f(connection(b), prompt))
  end

export render_view
render_view(name::String) =
  let path = prepare_for_saving_file(render_pathname(name))
    render_view(path, current_backend())
    path
  end


text:
"""
    EuclideanTensorSpace

A constant which is defined as `TensorSpace{<:EuclideanSpace}`, where
`TensorSpace{S<:ElementarySpace} = Union{S, ProductSpace{S}}`.
"""
const EuclideanTensorSpace = TensorSpace{<:EuclideanSpace}

"""
    EuclideanTensorMapSpace

A constant which is defined as `TensorMapSpace{<:EuclideanSpace}`, where
`TensorMapSpace{S<:ElementarySpace, N₁, N₂} =
    HomSpace{S, ProductSpace{S, N₁}, ProductSpace{S, N₂}}`.
"""
const EuclideanTensorMapSpace = TensorMapSpace{<:EuclideanSpace}

"""
    AbstractEuclideanTensorMap

A constant which is defined as `AbstractTensorMap{<:EuclideanTensorSpace}`.
"""
const AbstractEuclideanTensorMap = AbstractTensorMap{<:EuclideanTensorSpace}

"""
    EuclideanTensorMap

A constant which is defined as `TensorMap{<:EuclideanTensorSpace}`.
"""
const EuclideanTensorMap = TensorMap{<:EuclideanTensorSpace}

# Wrapping the blocks in a StridedView enables multithreading if JULIA_NUM_THREADS > 1.

"""
    copy!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap) -> tdst

Extend `Base.copy!`. In-place copy of `tsrc` to `tdst`.
"""
function Base.copy!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap)
    space(tdst) == space(tsrc) || throw(SpaceMismatch())
    for c in blocksectors(tdst)
        copy!(StridedView(block(tdst, c)), StridedView(block(tsrc, c)))
    end
    return tdst
end

"""
    copy(t::AbstractTensorMap)

Creat a new tenor map with data copied from the tensor map `t`.
"""
Base.copy(t::AbstractTensorMap) = Base.copy!(similar(t), t)

"""
    fill!(t::AbstractTensorMap, value::Number)

Extend `Base.fill!`. Fill the all data of the tensor map `t` by the `value`.
"""
function Base.fill!(t::AbstractTensorMap, value::Number)
    for (c, b) in blocks(t)
        fill!(b, value)
    end
    return t
end

"""
    adjoint!(tdst::AbstractEuclideanTensorMap, tsrc::AbstractEuclideanTensorMap)

Extend `LinearAlgebra.adjoint!`. Conjugate transpose the tensor map `tsrc` and store the
result in the preallocated tensor map `tdst`.
"""
function LinearAlgebra.adjoint!(tdst::AbstractEuclideanTensorMap,
                                tsrc::AbstractEuclideanTensorMap)
    space(tdst) == adjoint(space(tsrc)) || throw(SpaceMismatch())
    for c in blocksectors(tdst)
        adjoint!(StridedView(block(tdst, c)), StridedView(block(tsrc, c)))
    end
    return tdst
end

"""
    mul!(t1::AbstractTensorMap, t2::AbstractTensorMap, α::Number)

Extend `LinearAlgebra.mul!`. Multiply number `α` from right to `t2` and store the result in
`t1`.
"""
function LinearAlgebra.mul!(t1::AbstractTensorMap, t2::AbstractTensorMap, α::Number)
    space(t1) == space(t2) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        mul!(StridedView(block(t1, c)), StridedView(block(t2, c)), α)
    end
    return t1
end

"""
    mul!(t1::AbstractTensorMap, α::Number, t2::AbstractTensorMap)

Extend `LinearAlgebra.mul!`. Multiply number `α` from left to `t2` and store the result in
`t1`.
"""
function LinearAlgebra.mul!(t1::AbstractTensorMap, α::Number, t2::AbstractTensorMap)
    space(t1) == space(t2) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        mul!(StridedView(block(t1, c)), α, StridedView(block(t2, c)))
    end
    return t1
end

"""
    mul!(tC::AbstractTensorMap, tA::AbstractTensorMap,
            tB::AbstractTensorMap, α = true, β = false)

Extend `LinearAlgebra.mul!`. Overwrite `tC` by `tA*tB*α + tC*β`. Return `tC`.
"""
function LinearAlgebra.mul!(tC::AbstractTensorMap,
                            tA::AbstractTensorMap,
                            tB::AbstractTensorMap, α = true, β = false)
    if !(codomain(tC) == codomain(tA) && domain(tC) == domain(tB) && domain(tA) == codomain(tB))
        throw(SpaceMismatch())
    end
    for c in blocksectors(tC)
        if hasblock(tA, c)
            @assert hasblock(tB,c) # tB should have such a block
            A = block(tA, c)
            B = block(tB, c)
            C = block(tC, c)
            if isbitstype(eltype(A)) && isbitstype(eltype(B)) && isbitstype(eltype(C))
                @unsafe_strided A B C mul!(C, A, B, α, β)
            else
                mul!(StridedView(C), StridedView(A), StridedView(B), α, β)
            end
        elseif β != one(β)
            rmul!(block(tC, c), β)
        end
    end
    return tC
end

"""
    -(t::AbstractTensorMap)

Extend `Base.:-`. Creat a new tensor map with data corresponds to the minus of the tensor
map `t`.
"""
Base.:-(t::AbstractTensorMap) = mul!(similar(t), t, -one(eltype(t)))

"""
    *(t::AbstractTensorMap, α::Number)

Extend `Base.:*`. Creat a new tensor map that equals to `t*α`.
"""
Base.:*(t::AbstractTensorMap, α::Number) =
    mul!(similar(t, promote_type(eltype(t), typeof(α))), t, α)

"""
    *(α::Number, t::AbstractTensorMap)

Extend `Base.:*`. Creat a new tensor map that equals to `α*t`.
"""
Base.:*(α::Number, t::AbstractTensorMap) =
    mul!(similar(t, promote_type(eltype(t), typeof(α))), α, t)

"""
    *(t1::AbstractTensorMap, t2::AbstractTensorMap)

Extend `Base.:*`. Creat a new tensor map that equals to `t1*t2`.
"""
Base.:*(t1::AbstractTensorMap, t2::AbstractTensorMap) =
    mul!(similar(t1, promote_type(eltype(t1), eltype(t2)), codomain(t1)←domain(t2)), t1, t2)

"""
    rmul!(t::AbstractTensorMap, α::Number)

Extend `LinearAlgebra.rmul!`. Multiply number `α` from right to `t` and overwrite `t`
in-place.
"""
LinearAlgebra.rmul!(t::AbstractTensorMap, α::Number) = mul!(t, t, α)

"""
    lmul!(α::Number, t::AbstractTensorMap)

Extend `LinearAlgebra.lmul!`. Multiply number `α` from left to `t` and overwrite `t`
in-place.
"""
LinearAlgebra.lmul!(α::Number, t::AbstractTensorMap) = mul!(t, α, t)

"""
    axpy!(α::Number, t1::AbstractTensorMap, t2::AbstractTensorMap)

Extend `LinearAlgebra.axpy!`. Overwrite `t2` with `t1*α + t2`.
"""
function LinearAlgebra.axpy!(α::Number, t1::AbstractTensorMap, t2::AbstractTensorMap)
    space(t1) == space(t2) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        axpy!(α, StridedView(block(t1, c)), StridedView(block(t2, c)))
    end
    return t2
end

"""
    +(t1::AbstractTensorMap, t2::AbstractTensorMap)

Extend `Base.:+`. Creat a new tensor map that equals to `t1 + t2`.
"""
function Base.:+(t1::AbstractTensorMap, t2::AbstractTensorMap)
    T = promote_type(eltype(t1), eltype(t2))
    return axpy!(one(T), t2, copy!(similar(t1, T), t1))
end

"""
    -(t1::AbstractTensorMap, t2::AbstractTensorMap)

Extend `Base.:-`: Creat a new tensor map that equals to `t1 - t2`.
"""
function Base.:-(t1::AbstractTensorMap, t2::AbstractTensorMap)
    T = promote_type(eltype(t1), eltype(t2))
    return axpy!(-one(T), t2, copy!(similar(t1, T), t1))
end

"""
    axpby!(α::Number, t1::AbstractTensorMap, β::Number, t2::AbstractTensorMap)

Extend `LinearAlgebra.axpby!`. Overwrite `t2` with `t1*α + t2*β`. Return `t2`.
"""
function LinearAlgebra.axpby!(α::Number, t1::AbstractTensorMap,
                                β::Number, t2::AbstractTensorMap)
    space(t1) == space(t2) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        axpby!(α, StridedView(block(t1, c)), β, StridedView(block(t2, c)))
    end
    return t2
end

"""
    ^(t::AbstractTensorMap, p::Integer)

Extend `Base.^`. Creat a new tensor map that equals to `t*t*t ...` (`p` multiplications).

Note: `Base.power_by_squaring(x,p)` is calculating `^` for any `x` supporting `*`.
"""
Base.:^(t::AbstractTensorMap, p::Integer) =
    p < 0 ? Base.power_by_squaring(inv(t), -p) : Base.power_by_squaring(t, p)

"""
    tr(t::AbstractTensorMap)

Extend `LinearAlgebra.tr`. Return the trace of the true block diagonal matrix that
represent the tensor map.
"""
function LinearAlgebra.tr(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SpaceMismatch("Trace of a tensor only exist when domain == codomain"))
    return sum(dim(c)*tr(b) for (c, b) in blocks(t))
end

"""
    exp!(t::TensorMap)

Overwrite each matrix in the data of the tensor map `t` with its exponential.
"""
function exp!(t::TensorMap)
    domain(t) == codomain(t) ||
        error("Exponentional of a tensor only exist when domain == codomain.")
    for (c, b) in blocks(t)
        copy!(b, LinearAlgebra.exp!(b))
    end
    return t
end

"""
    exp(t::AbstractTensorMap)

Extend `Base.exp`. Return a new tensor map with data matrices which are the
exponentials of that of `t`.
"""
Base.exp(t::AbstractTensorMap) = exp!(copy(t))

"""
    inv(t::AbstractTensorMap)

Extend `Base.inv`. Creat a new tensor map that similar to `t` with data matries correspond
to the inverse of that of `t`.
"""
function Base.inv(t::AbstractTensorMap)
    cod = codomain(t)
    dom = domain(t)
    for c in union(blocksectors(cod), blocksectors(dom))
        blockdim(cod, c) == blockdim(dom, c) ||
            throw(SpaceMismatch("codomain $cod and domain $dom are not isomorphic: no inverse"))
    end
    if sectortype(t) === Trivial
        return TensorMap(inv(block(t, Trivial())), domain(t)←codomain(t))
    else
        data = empty(t.data)
        for (c, b) in blocks(t)
            data[c] = inv(b)
        end
        return TensorMap(data, domain(t)←codomain(t))
    end
end

"""
    pinv(t::AbstractTensorMap; kwargs...)

Extend `LinearAlgebra.pinv`. Return a new tensor map with data which is the Moore-Penrose
pseudoinverse of the corresponding data matrices of tensor map `t`.
"""
function LinearAlgebra.pinv(t::AbstractTensorMap; kwargs...)
    if sectortype(t) === Trivial
        return TensorMap(pinv(block(t, Trivial()); kwargs...), domain(t)←codomain(t))
    else
        data = empty(t.data)
        for (c, b) in blocks(t)
            data[c] = pinv(b; kwargs...)
        end
        return TensorMap(data, domain(t)←codomain(t))
    end
end

"""
    \\(t1::AbstractTensorMap, t2::AbstractTensorMap)

Extend `Base.:(\\)`.

It means matrix division using a polyalgorithm. For input matrices `A` and `B`, the
result `X` is such that `A*X == B` when `A` is square.

Return the tensor map from `domain(t2)` to `domain(t1)` with the data matrices which is the
matrix division between the correpsonding matrices of `t1` and `t2`.
"""
function Base.:(\)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    codomain(t1) == codomain(t2) ||
        throw(SpaceMismatch("non-matching codomains in t1 \\ t2"))
    if sectortype(t1) === Trivial
        data = block(t1, Trivial()) \ block(t2, Trivial())
        return TensorMap(data, domain(t1)←domain(t2))
    else
        cod = codomain(t1)
        data = SectorDict(c=>block(t1, c) \ block(t2, c) for c in blocksectors(codomain(t1)))
        return TensorMap(data, domain(t1)←domain(t2))
    end
end

"""
    /(t1::AbstractTensorMap, t2::AbstractTensorMap)

Extend `Base.:(/)`.

`/(A, B)` means matrix division using a polyalgorithm. For input matrices `A` and `B`, the
result `X` is such that `A == X*B` when `A` is square.

Return the tensor map from `codomain(t2)` to `codomain(t1)` with the data matrices which is
the matrix division `/` between the correpsonding matrices of `t1` and `t2`.
"""
function Base.:(/)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    domain(t1) == domain(t2) ||
        throw(SpaceMismatch("non-matching domains in t1 / t2"))
    if sectortype(t1) === Trivial
        data = block(t1, Trivial()) / block(t2, Trivial())
        return TensorMap(data, codomain(t1)←codomain(t2))
    else
        data = SectorDict(c=>block(t1, c) / block(t2, c) for c in blocksectors(domain(t1)))
        return TensorMap(data, codomain(t1)←codomain(t2))
    end
end

"""
    /(t::AbstractTensorMap, α::Number)

Extend `Base.:/`. Return a new tensor map with data matrices as `1/α` of that of `t`.
"""
Base.:/(t::AbstractTensorMap, α::Number) = *(t, one(eltype(t))/α)

"""
    \\(α::Number, t::AbstractTensorMap)

Return a new tensor map with data matrices as `1/α` of that of `t`.
"""
Base.:\(α::Number, t::AbstractTensorMap) = *(t, one(eltype(t))/α)

# functions that map ℝ to (a subset of) ℝ
for f in (:cos, :sin, :tan, :cot, :cosh, :sinh, :tanh, :coth, :atan, :acot, :asinh)
    sf = string(f)
    @eval function Base.$f(t::AbstractTensorMap)
        domain(t) == codomain(t) ||
            error("$sf of a tensor only exist when domain == codomain.")
        I = sectortype(t)
        T = similarstoragetype(t, float(eltype(t)))
        if sectortype(t) === Trivial
            local data::T
            if eltype(t) <: Real
                data = real($f(block(t, Trivial())))
            else
                data = $f(block(t, Trivial()))
            end
            return TensorMap(data, codomain(t), domain(t))
        else
            if eltype(t) <: Real
                datadict = SectorDict{I, T}(c=>real($f(b)) for (c, b) in blocks(t))
            else
                datadict = SectorDict{I, T}(c=>$f(b) for (c, b) in blocks(t))
            end
            return TensorMap(datadict, codomain(t), domain(t))
        end
    end
end

# functions that don't map ℝ to (a subset of) ℝ
for f in (:sqrt, :log, :asin, :acos, :acosh, :atanh, :acoth)
    sf = string(f)
    @eval function Base.$f(t::AbstractTensorMap)
        domain(t) == codomain(t) ||
            error("$sf of a tensor only exist when domain == codomain.")
        I = sectortype(t)
        T = similarstoragetype(t, complex(float(eltype(t))))
        if sectortype(t) === Trivial
            data::T = $f(block(t, Trivial()))
            return TensorMap(data, codomain(t), domain(t))
        else
            datadict = SectorDict{I, T}(c=>$f(b) for (c, b) in blocks(t))
            return TensorMap(datadict, codomain(t), domain(t))
        end
    end
end

# inner product and norm only valid for spaces with Euclidean inner product:
"""
    dot(t1::AbstractEuclideanTensorMap, t2::AbstractEuclideanTensorMap)

Extend `LinearAlgebra.dot`. Return the elementwise dot product of the data between two
tensor maps. The elements of `t1` are conjugated before multiplication. And the dimension
of blocksector `c` is multiplied in the after we make dot operation between the two
matrices corresponding to block `c` from `t1` and `t2`. Thus, the dot operation is between
the true block diagonalized tensor map data.
"""
function LinearAlgebra.dot(t1::AbstractEuclideanTensorMap, t2::AbstractEuclideanTensorMap)
    space(t1) == space(t2) || throw(SpaceMismatch())
    T = promote_type(eltype(t1), eltype(t2))
    s = zero(T)
    for c in blocksectors(t1)
        s += convert(T, dim(c)) * dot(block(t1, c), block(t2, c))
    end
    return s
end

"""
    norm(t::AbstractEuclideanTensorMap, p::Real = 2)

Extend `LinearAlgebra.norm`. Return the norm of the tensor map `t` as the norm of the true
block diagonal matrix which representing the tensor map.

For any iterable container `A` (including arrays of any dimension) of numbers, the `p`-norm
is defined as
```math
\\|A\\|_p = \\left(\\sum_{i=1}^n |a_i|^p \\right)^{1/p}
```
with ``a_i`` the entries of ``A``, ``|a_i|`` the `norm` of ``a_i``, and
``n`` the length of ``A``.

The `norm(A, Inf)` returns the largest value in `abs.(A)`. If `A` is a matrix and `p=2`,
then this is equivalent to the Frobenius norm.
"""
LinearAlgebra.norm(t::AbstractEuclideanTensorMap, p::Real = 2) =
    _norm(blocks(t), p, float(zero(real(eltype(t)))))

function _norm(blockiter, p::Real, init::Real)
    if p == Inf
        return mapreduce(max, blockiter; init = init) do (c, b)
            isempty(b) ? init : oftype(init, LinearAlgebra.normInf(b))
        end
    elseif p == 2
        return sqrt(mapreduce(+, blockiter; init = init) do (c, b)
            isempty(b) ? init : oftype(init, dim(c)*LinearAlgebra.norm2(b)^2)
        end)
    elseif p == 1
        return mapreduce(+, blockiter; init = init) do (c, b)
            isempty(b) ? init : oftype(init, dim(c)*sum(abs, b))
        end
    elseif p > 0
        s = mapreduce(+, blockiter; init = init) do (c, b)
            isempty(b) ? init : oftype(init, dim(c)*LinearAlgebra.normp(b, p)^p)
        end
        return s^inv(oftype(s, p))
    else
        msg = "Norm with non-positive p is not defined for `AbstractTensorMap`"
        throw(ArgumentError(msg))
    end
end

"""
    normalize!(t::AbstractTensorMap, p::Real = 2)

Extend `LinearAlgebra.normalize!`. Replace the tensor map `t` with the normalize one which
has `norm(t,p)==1`.
"""
LinearAlgebra.normalize!(t::AbstractTensorMap, p::Real = 2) = rmul!(t, inv(norm(t, p)))

"""
    normalize(t::AbstractTensorMap, p::Real = 2)

Extend `LinearAlgebra.normalize`. Creat a new tensor map that is similar to `t` and has the
same data with  `normalize!(t, p)`.
"""
LinearAlgebra.normalize(t::AbstractTensorMap, p::Real = 2) =
    mul!(similar(t), t, inv(norm(t, p)))

"""
    sylvester(A::AbstractTensorMap, B::AbstractTensorMap, C::AbstractTensorMap)

Extend `LinearAlgebra.sylvester`.

For matrices `A`, `B`, and `C`, it computes the solution `X` to the Sylvester equation
`AX + XB + C = 0`, where `A`, `B` and `C` have compatible dimensions and `A` and `-B`
have no eigenvalues with equal real part.

Return a new tensor map with each block matrix corresponds to the solution `X` obtained
from the corresponding matrices of tensor maps `A`, `B` and `C`.
"""
function LinearAlgebra.sylvester(A::AbstractTensorMap, B::AbstractTensorMap,
                                    C::AbstractTensorMap)
    (codomain(A) == domain(A) == codomain(C) && codomain(B) == domain(B) == domain(C)) ||
        throw(SpaceMismatch())
    cod = domain(A)
    dom = codomain(B)
    sylABC(c) = sylvester(block(A, c), block(B, c), block(C, c))
    data = SectorDict(c=>sylABC(c) for c in blocksectors(cod ← dom))
    return TensorMap(data, cod ← dom)
end

"""
    catdomain(t1::AbstractTensorMap{S, N₁, 1},
                t2::AbstractTensorMap{S, N₁, 1}) where {S, N₁}

Return the horizontally concatenation of two tensor maps with length `1` domain and the same
length of codomain. A new tensor map is created, whose domain is a product space with one
vector space which is the direct sum of two vector spaces in the domain of `t1` and `t2`,
and the codomain is still the same with `t1` (or `t2`).
"""
function catdomain(t1::AbstractTensorMap{S, N₁, 1},
                    t2::AbstractTensorMap{S, N₁, 1}) where {S, N₁}

    codomain(t1) == codomain(t2) || throw(SpaceMismatch())

    V1, = domain(t1)
    V2, = domain(t2)
    isdual(V1) == isdual(V2) ||
        throw(SpaceMismatch("cannot horizontally concatenate tensors whose domain has
                                non-matching duality"))

    V = V1 ⊕ V2
    t = TensorMap(undef, promote_type(eltype(t1), eltype(t2)), codomain(t1), V)
    for c in sectors(V)
        block(t, c)[:, 1:dim(V1, c)] .= block(t1, c)
        block(t, c)[:, dim(V1, c) .+ (1:dim(V2, c))] .= block(t2, c)
    end
    return t
end

"""
    catcodomain(t1::AbstractTensorMap{S, 1, N₂},
                t2::AbstractTensorMap{S, 1, N₂}) where {S, N₂}

Return the vertically concatenation of two tensor maps with length `1` codomain and the same
length of domain. A new tensor map is created, whose codomain is a product space with one
vector space which is the direct sum of two vector spaces in the codomain of `t1` and `t2`,
and the domain is still the same with `t1` (or `t2`).
"""
function catcodomain(t1::AbstractTensorMap{S, 1, N₂},
                        t2::AbstractTensorMap{S, 1, N₂}) where {S, N₂}
    domain(t1) == domain(t2) || throw(SpaceMismatch())

    V1, = codomain(t1)
    V2, = codomain(t2)
    isdual(V1) == isdual(V2) ||
        throw(SpaceMismatch("cannot vertically concatenate tensors whose codomain has
                                non-matching duality"))

    V = V1 ⊕ V2
    t = TensorMap(undef, promote_type(eltype(t1), eltype(t2)), V, domain(t1))
    for c in sectors(V)
        block(t, c)[1:dim(V1, c), :] .= block(t1, c)
        block(t, c)[dim(V1, c) .+ (1:dim(V2, c)), :] .= block(t2, c)
    end
    return t
end

"""
    ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}, ...) -> TensorMap{S}

Compute the tensor product between two `AbstractTensorMap` instances, which results in a
new `TensorMap` instance whose codomain is `codomain(t1) ⊗ codomain(t2)` and whose domain
is `domain(t1) ⊗ domain(t2)`.
"""
function ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}) where S
    cod1, cod2 = codomain(t1), codomain(t2)
    dom1, dom2 = domain(t1), domain(t2)
    cod = cod1 ⊗ cod2
    dom = dom1 ⊗ dom2
    t = TensorMap(zeros, promote_type(eltype(t1), eltype(t2)), cod, dom)
    if sectortype(S) === Trivial
        d1 = dim(cod1)
        d2 = dim(cod2)
        d3 = dim(dom1)
        d4 = dim(dom2)
        m1 = reshape(t1[], (d1, 1, d3, 1))
        m2 = reshape(t2[], (1, d2, 1, d4))
        m = reshape(t[], (d1, d2, d3, d4))
        m .= m1 .* m2
    else
        for (f1l, f1r) in fusiontrees(t1)
            for (f2l, f2r) in fusiontrees(t2)
                c1 = f1l.coupled # = f1r.coupled
                c2 = f2l.coupled # = f2r.coupled
                for c in c1 ⊗ c2
                    degeneracyiter = FusionStyle(c) isa GenericFusion ?
                                        (1:Nsymbol(c1, c2, c)) : (nothing,)
                    for μ in degeneracyiter
                        for (fl, coeff1) in merge(f1l, f2l, c, μ)
                            for (fr, coeff2) in merge(f1r, f2r, c, μ)
                                d1 = dim(cod1, f1l.uncoupled)
                                d2 = dim(cod2, f2l.uncoupled)
                                d3 = dim(dom1, f1r.uncoupled)
                                d4 = dim(dom2, f2r.uncoupled)
                                m1 = reshape(t1[f1l, f1r], (d1, 1, d3, 1))
                                m2 = reshape(t2[f2l, f2r], (1, d2, 1, d4))
                                m = reshape(t[fl, fr], (d1, d2, d3, d4))
                                m .+= coeff1 .* coeff2 .* m1 .* m2
                            end
                        end
                    end
                end
            end
        end
    end
    return t
end

"""
    ⊠(t1::AbstractTensorMap{<:EuclideanSpace{ℂ}},
        t2::AbstractTensorMap{<:EuclideanSpace{ℂ}})

Return the deligne product of tensors.
"""
function ⊠(t1::AbstractTensorMap{<:EuclideanSpace{ℂ}},
            t2::AbstractTensorMap{<:EuclideanSpace{ℂ}})
    S1 = spacetype(t1)
    I1 = sectortype(S1)
    S2 = spacetype(t2)
    I2 = sectortype(S2)
    codom1 = codomain(t1) ⊠ one(S2)
    dom1 = domain(t1) ⊠ one(S2)
    data1 = SectorDict{I1 ⊠ I2, storagetype(t1)}(c ⊠ one(I2) => b for (c,b) in blocks(t1))
    t1′ = TensorMap(data1, codom1, dom1)
    codom2 = one(S1) ⊠ codomain(t2)
    dom2 = one(S1) ⊠ domain(t2)
    data2 = SectorDict{I1 ⊠ I2, storagetype(t2)}(one(I1) ⊠ c => b for (c,b) in blocks(t2))
    t2′ = TensorMap(data2, codom2, dom2)
    return t1′ ⊗ t2′
end

# Special purpose constructors:
"""
    zero(t::AbstractTensorMap)

Extend `Base.zero`. Creat a tensor that is similar to the tensor map `t` with all `0` in the
data.
"""
Base.zero(t::AbstractTensorMap) = fill!(similar(t), 0)

"""
    one!(t::AbstractTensorMap)

Overwrite the tensor map `t` by a tensor map in which every matrix in data is an identity
matrix.
"""
function one!(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SectorMismatch("no identity if domain and codomain are different"))
    for (c, b) in blocks(t)
        _one!(b)
    end
    return t
end

"""
    one(t::AbstractTensorMap)

Extend `Base.one`. Creat a tensor map that similar to tensor map `t` and with identity
matrices in data.
"""
function Base.one(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SectorMismatch("no identity if domain and codomain are different"))
    return one!(similar(t))
end

"""
    id([A::Type{<:DenseMatrix} = Matrix{Float64},] space::VectorSpace) -> TensorMap

Construct the identity endomorphism on space `space`, i.e. return a `t::TensorMap` with
`domain(t) == codomain(t) == V`, where `storagetype(t) = A` can be specified.
"""
id(::Type{A}, P::ProductSpace) where {A<:DenseMatrix} =one!(TensorMap(s->A(undef, s), P, P))
id(V::VectorSpace) = id(Matrix{Float64}, V)
id(A, V::ElementarySpace) = id(A, ProductSpace(V))

"""
    isomorphism([A::Type{<:DenseMatrix} = Matrix{Float64},]
                    cod::VectorSpace, dom::VectorSpace) -> TensorMap

Return a `t::TensorMap` that implements a specific isomorphism between the codomain `cod`
and the domain `dom`, and for which `storagetype(t)` can optionally be chosen to be of type
`A`. If the two spaces do not allow for such an isomorphism, and are thus not isomorphic,
and error will be thrown. When they are isomorphic, there is no canonical choice for a
specific isomorphism, but the current choice is such that
`isomorphism(cod, dom) == inv(isomorphism(dom, cod))`.

See also [`unitary`](@ref) when `spacetype(cod) isa EuclideanSpace`.
"""
function isomorphism(::Type{A}, cod::ProductSpace, dom::ProductSpace) where {A<:DenseMatrix}
    cod ≅ dom || throw(SpaceMismatch("codomain $cod and domain $dom are not isomorphic"))
    t = TensorMap(s->A(undef, s), cod, dom)
    for (c, b) in blocks(t)
        _one!(b)
    end
    return t
end
isomorphism(cod::TensorSpace, dom::TensorSpace) = isomorphism(Matrix{Float64}, cod, dom)
isomorphism(P::TensorMapSpace) = isomorphism(codomain(P), domain(P))
isomorphism(A::Type{<:DenseMatrix}, P::TensorMapSpace) =
    isomorphism(A, codomain(P), domain(P))
isomorphism(A::Type{<:DenseMatrix}, cod::TensorSpace, dom::TensorSpace) =
    isomorphism(A, convert(ProductSpace, cod), convert(ProductSpace, dom))

"""
    unitary([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace)
    -> TensorMap

Return a `t::TensorMap` that implements a specific unitary isomorphism between the codomain
`cod` and the domain `dom`, for which `spacetype(dom)` (`== spacetype(cod)`) must be a
subtype of `EuclideanSpace`. Furthermore, `storagetype(t)` can optionally be chosen to be
of type `A`. If the two spaces do not allow for such an isomorphism, and are thus not
isomorphic, and error will be thrown. When they are isomorphic, there is no canonical choice
for a specific isomorphism, but the current choice is such that
`unitary(cod, dom) == inv(unitary(dom, cod)) = adjoint(unitary(dom, cod))`.
"""
unitary(cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) = isomorphism(cod, dom)
unitary(P::EuclideanTensorMapSpace) = isomorphism(P)
unitary(A::Type{<:DenseMatrix}, P::EuclideanTensorMapSpace) = isomorphism(A, P)
unitary(A::Type{<:DenseMatrix}, cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) =
    isomorphism(A, cod, dom)

"""
    isometry([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace)
    -> TensorMap

Return a `t::TensorMap` that implements a specific isometry that embeds the domain `dom`
into the codomain `cod`, and which requires that `spacetype(dom)` (`== spacetype(cod)`) is
a subtype of `EuclideanSpace`. An isometry `t` is such that its adjoint `t'` is the left
inverse of `t`, i.e. `t'*t = id(dom)`, while `t*t'` is some idempotent endomorphism of
`cod`, i.e. it squares to itself. When `dom` and `cod` do not allow for such an isometric
inclusion, an error will be thrown.
"""
function isometry(::Type{A}, cod::ProductSpace{S},
                    dom::ProductSpace{S}) where {A<:DenseMatrix, S<:EuclideanSpace}
    dom ≾ cod || throw(SpaceMismatch("codomain $cod and domain $dom do not allow for an
                                        isometric mapping"))
    t = TensorMap(s->A(undef, s), cod, dom)
    for (c, b) in blocks(t)
        _one!(b)
    end
    return t
end
isometry(cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) =
    isometry(Matrix{Float64}, cod, dom)
isometry(P::EuclideanTensorMapSpace) = isometry(codomain(P), domain(P))
isometry(A::Type{<:DenseMatrix}, P::EuclideanTensorMapSpace) =
    isometry(A, codomain(P), domain(P))
isometry(A::Type{<:DenseMatrix}, cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) =
    isometry(A, convert(ProductSpace, cod), convert(ProductSpace, dom))

ext: jl


text:
function templateTestTmap(arr::AbstractArray)::Bool
	return all(map(x->2x, arr) .== tmap(x->2x, arr));
end

function testTmapStable()::Bool
	return templateTestTmap(randArray());
end

function testTmapUnstable()::Bool
	arr = randStrArray();
	str = lipsum();
	return all(map(x->str[1:x], arr) .== tmap(x->str[1:x], arr));
end

function testTmapOffset()::Bool
	return templateTestTmap(offsetArray());
end

function testTmapEmpty()::Bool
	return templateTestTmap([]);
end

#####

function templateTestTmapMultiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	return all(map((x, y)->2x + 3y, arr1, arr2) .== tmap((x, y)->2x + 3y, arr1, arr2));
end

function testTmapStableMultiple()::Bool
	return templateTestTmapMultiple(rand2Arrays()...);
end

function testTmapUnstableMultiple()::Bool
	arr1, arr2 = rand2StrArrays();
	str = lipsum();
	return all(map((x, y)->str[1:x + y], arr1, arr2) .== tmap((x, y)->str[1:x + y], arr1, arr2));
end

function testTmapMultipleOffset()::Bool
	return templateTestTmapMultiple(offsetArray(), offsetArray());
end

function testTmapMultipleEmpty()::Bool
	return templateTestTmapMultiple([], []);
end

#####

function testTmap!()::Bool
	n = rand(1:2);
	srcDim = rand(1:3, n);
	destDim = srcDim[randperm(n)];
	srcArr = rand(1:1000, srcDim...);
	destArr = rand(1:1000, destDim...);
	tSrcArr = deepcopy(srcArr);
	tDestArr = deepcopy(destArr);
	map!(x->2x, destArr, srcArr);
	tmap!(x->2x, tDestArr, tSrcArr);
	return all(destArr .== tDestArr);
end

function testTmap!Empty()::Bool
	destDim = rand(1:3, rand(1:2));
	destArr = rand(1:1000, destDim...);
	tDestArr = deepcopy(destArr);
	map!(x->2x, destArr, []);
	tmap!(x->2x, tDestArr, []);
	return all(destArr .== tDestArr);
end

#####

function testTmap!Multiple()::Bool
	n = rand(1:2);
	srcDim = rand(1:3, n);
	destDim = srcDim[randperm(n)];
	srcArr1 = rand(1:1000, srcDim...);
	srcArr2 = rand(1:1000, srcDim...);
	destArr = rand(1:1000, destDim...);
	tSrcArr1 = deepcopy(srcArr1);
	tSrcArr2 = deepcopy(srcArr2);
	tDestArr = deepcopy(destArr);
	map!((x, y)->2x + 3y, destArr, srcArr1, srcArr2);
	tmap!((x, y)->2x + 3y, tDestArr, tSrcArr1, tSrcArr2);
	return all(destArr .== tDestArr);
end

function testTmap!MultipleEmpty()::Bool
	destDim = rand(1:3, rand(1:2));
	destArr = rand(1:1000, destDim...);
	tDestArr = deepcopy(destArr);
	map!((x, y)->2x + 3y, destArr, [], []);
	tmap!((x, y)->2x + 3y, tDestArr, [], []);
	return all(destArr .== tDestArr);
end

function testTmapMultipleDimensionMismatch()
	tmap((x, y)->2x + 3y, [1, 2], [3, 4, 5]);
end


text:
# these functions are indicators, an event occurs between their value change
nbSigns(x, ::AbstractContinuousEvent) = mapreduce(x -> x > 0, +, x)
nbSigns(x, ::AbstractDiscreteEvent) = x

function nbSigns(x, event::PairOfEvents)
	xc = x[1:event.eventC.nb]
	xd = x[event.eventC.nb+1:end]
	res = (nbSigns(xc, event.eventC)..., nbSigns(xd, event.eventD)...)
	return res
end

function nbSigns(x, event::SetOfEvents)
	nb = 0
	nC = length(event.eventC)
	nD = length(event.eventD)
	nCb = length(x)
	@inbounds for i in 1:nC
		nb += nbSigns(x[i], event.eventC[i])
	end
	@inbounds for i in nC+1:nCb
		nb += nbSigns(x[i], event.eventC[i-nC])
	end
	return nb
end
####################################################################################################
# Function to locate precisely an Event using a bisection algorithm. We make sure that, at the end of the algorithm, the state is just after the event (in the s coordinate).
# I put the event in first argument even if it is in `iter` in order to allow for easier dispatch
function locateEvent!(event::AbstractEvent, iter, _state, verbose::Bool = true)
	@assert isnothing(_state.eventValue) == false "Empty event value, this should not be happening. Please open an issue."

	# type of scalars in iter
	_T = eltype(iter)
	if abs(_state.ds) < iter.contParams.dsmin; return :none, (_T(0), _T(0)); end

	# get continuation parameters
	contParams = iter.contParams

	n2, n1 = nbSigns(_state.eventValue[1], event), nbSigns(_state.eventValue[2], event)
	verbose && println("----> Entering [Event], indicator of 2 last events = ", (n2, n1))
	verbose && println("----> [Bisection] initial ds = ", _state.ds)

	# we create a new state copy for stepping through the continuation routine
	after = copy(_state)	# after the bifurcation point
	state = copy(_state)	# current state of the bisection
	before = copy(_state)	# before the bifurcation point

	# we reverse some indicators for `before`. It is OK, it will never be used other than for getp(before)
	before.n_unstable = (before.n_unstable[2], before.n_unstable[1])
	before.n_imag = (before.n_imag[2], before.n_imag[1])
	before.eventValue = (before.eventValue[2], before.eventValue[1])
	before.z_pred.p, before.z.p = before.z.p, before.z_pred.p

	# the bifurcation point is before the current state
	# so we want to first iterate backward
	# we turn off stepsizecontrol because it would not be a
	# bisection otherwise
	state.ds *= -1
	state.step = 0
	state.stepsizecontrol = false

	next = (state, state)

	# record sequence of event indicators
	nsigns = [n2]

	# interval which contains the event
	interval = getinterval(getp(state), state.z_pred.p)
	# index of active index in the bisection interval, allows to track interval
	indinterval = interval[1] == getp(state) ? 1 : 2

	verbose && println("----> [Bisection] state.ds = ", state.ds)

	# we put this to be able to reference it at the end of this function
	# we don't know its type yet
	eiginfo = nothing

	# we compute the number of changes in event indicator
	n_inversion = 0
	status = :guess

	eventlocated = false

	# for a polynomial tangent predictor, we disable the update of the predictor parameters
	# TODO Find better way to do this
	if iter.tangentAlgo isa PolynomialPred; iter.tangentAlgo.update = false; end

	verbose && printstyled(color=:green, "--> eve (initial) ",
		_state.eventValue[2], " --> ",  _state.eventValue[1], "\n")
	if verbose && ~isnothing(_state.eigvals)
		printstyled(color=:green, "\n--> eigvals = \n")
		displayEV(_state.eigvals, :green)
		# calcul des VP et determinant
	end

	# emulate a do-while
	while true
		if ~state.isconverged
			@error "Newton failed to fully locate bifurcation point using bisection parameters!"
			break
		 end

		# if PALC stops, break the bisection
		if isnothing(next)
			break
		end

		# perform one continuation step
		(_, state) = next

		# the eigenelements have been computed/stored in state during the call iterate(iter, state)
		updateEvent!(iter, state)
		push!(nsigns, nbSigns(state.eventValue[1], event))
		verbose && printstyled(color=:green, "\n----> eve (current) ",
			state.eventValue[2], " --> ", state.eventValue[1], "\n")

		if verbose && ~isnothing(state.eigvals)
			printstyled(color=:blue, "----> eigvals = \n")
			displayEV(state.eigvals, :blue)
		end

		if nsigns[end] == nsigns[end-1]
			# event still after current state, keep going
			state.ds /= 2
		else
			# we passed the event, reverse continuation
			state.ds /= -2
			n_inversion += 1
			indinterval = (indinterval == 2) ? 1 : 2
		end

		if iseven(n_inversion)
			copyto!(after, state)
		else
			copyto!(before, state)
		end

		# update the interval containing the event
		state.step > 0 && (@set! interval[indinterval] = getp(state))

		# we call the finalizer
		state.stopcontinuation = ~iter.finaliseSolution(state.z, state.τ, state.step, nothing; bisection = true, state = state)

		if verbose
			printstyled(color=:blue, bold = true, "----> ", state.step,
				" - [Bisection] (n1, n_current, n2) = ", (n1, nsigns[end], n2),
				"\n\t\t\tds = ", state.ds, ", p = ", getp(state), ", #reverse = ", n_inversion,
				"\n----> event ∈ ", getinterval(interval...),
				", precision = ", @sprintf("%.3E", interval[2] - interval[1]), "\n")
		end

		eventlocated = (isEventCrossed(event, iter, state) &&
				abs(interval[2] - interval[1]) < contParams.tolParamBisectionEvent)

		# condition for breaking the while loop
		if (isnothing(next) == false &&
				abs(state.ds) >= contParams.dsminBisection &&
				state.step < contParams.maxBisectionSteps &&
				n_inversion < contParams.nInversion &&
				eventlocated == false) == false
			break
		end

		next = iterate(iter, state; _verbosity = 0)
	end

	if verbose
		printstyled(color=:red, "----> Found at p = ", getp(state), " ∈ $interval, \n\t\t\t  δn = ", abs.(2 .* nsigns[end] .- n1 .- n2), ", from p = ", getp(_state), "\n")
		printstyled(color=:blue, "-"^40*"\n----> Stopping reason:\n------> isnothing(next)           = ", isnothing(next),
				"\n------> |ds| < dsminBisection     = ", abs(state.ds) < contParams.dsminBisection,
				"\n------> step >= maxBisectionSteps = ", state.step >= contParams.maxBisectionSteps,
				"\n------> n_inversion >= nInversion = ", n_inversion >= contParams.nInversion,
				"\n------> eventlocated              = ", eventlocated == true, "\n")

	end

	if iter.tangentAlgo isa PolynomialPred
		iter.tangentAlgo.update = true
	end

	######## update current state ########
	# So far we have (possibly) performed an even number of event crossings.
	# We started at the right of the event point. The current state is thus at the
	# right of the event point if iseven(n_inversion) == true. Otherwise, the event
	# point is deemed undetected.
	# In the case n_inversion = 0, we are still on the right of the event point
	if iseven(n_inversion)
		status = n_inversion >= contParams.nInversion ? :converged : :guess
		copyto!(_state.z_pred, state.z_pred)
		copyto!(_state.z,  state.z)
		copyto!(_state.τ, state.τ)
		# if there is no inversion, the eventValue will possibly be constant like (0, 0). Hence

		if computeEigenElements(iter.event)
			# save eigen-elements
			_state.eigvals = state.eigvals
			if saveEigenvectors(contParams)
				_state.eigvecs = state.eigvecs
			end
			# to prevent spurious event detection, update the following numbers carefully
			_state.n_unstable = (state.n_unstable[1], before.n_unstable[1])
			_state.n_imag = (state.n_imag[1], before.n_imag[1])
		end

		_state.eventValue = (state.eventValue[1], before.eventValue[1])
		interval = (getp(state), getp(before))
	else
		status = :guessL
		copyto!(_state.z_pred, after.z_pred)
		copyto!(_state.z,  after.z)
		copyto!(_state.τ, after.τ)

		if computeEigenElements(iter.event)
			# save eigen-elements
			_state.eigvals = after.eigvals
			if contParams.saveEigenvectors
				_state.eigvecs = after.eigvecs
			end
			# to prevent spurious event detection, update the following numbers carefully
			_state.n_unstable = (after.n_unstable[1], state.n_unstable[1])
			_state.n_imag = (after.n_imag[1], state.n_imag[1])
		end

		_state.eventValue = (after.eventValue[1], state.eventValue[1])
		interval = (getp(state), getp(after))
	end
	verbose && println("----> Leaving [Loc-Bif]")
	return status, getinterval(interval...)
end
####################################################################################################
####################################################################################################
# because of the way the results are recorded, with state corresponding to the (continuation) step = 0 saved in br.branch[1], it means that br.eig[k] corresponds to state.step = k-1. Thus, the eigen-elements (and other information) corresponding to the current event point are saved in br.eig[step+1]
EventSpecialPoint(state::ContState, Utype::Symbol, status::Symbol, printsolution, normC, interval) = SpecialPoint(state, Utype, status, printsolution, normC, interval; idx = state.step + 1)

# I put the callback in first argument even if it is in iter in order to allow for dispatch
# function to tell the event type based  on the coordinates of the zero
function getEventType(event::AbstractEvent, iter::AbstractContinuationIterable, state, verbosity, status::Symbol, interval::Tuple{T, T}, ind = :) where T
	# record information about the event point
	userpoint = EventSpecialPoint(state, :user, status, iter.recordFromSolution, iter.normC, interval)
	(verbosity > 0) && printstyled(color=:red, "!! User point at p ≈ $(getp(state)) \n")
	return true, userpoint
end
####################################################################################################
function getEventType(event::AbstractContinuousEvent, iter::AbstractContinuationIterable, state, verbosity, status::Symbol, interval::Tuple{T, T}, ind = :; typeE = "userC") where T
	event_index_C = Int32[]
	if state.eventValue[1] isa Real
		if (state.eventValue[1] * state.eventValue[2] < 0)
			push!(event_index_C,1)
		end
	elseif state.eventValue[1][ind] isa Real
		if state.eventValue[1][ind] * state.eventValue[2][ind] < 0
			push!(event_index_C,1)
		end
	else
		for ii in eachindex(state.eventValue[1][ind])
			if state.eventValue[1][ind][ii] * state.eventValue[2][ind][ii] < 0
				push!(event_index_C,ii)
				typeE = typeE * "-$ii"
			end
		end
	end
	if isempty(event_index_C) == true
		@error "Error, no event was characterized whereas one was detected. Please open an issue at https://github.com/rveltz/BifurcationKit.jl/issues. \n The events are eventValue = $(state.eventValue)"
		# we halt continuation as it will mess up the detection of events
		state.stopcontinuation = true
		return false, EventSpecialPoint(state, Symbol(typeE), status, iter.recordFromSolution, iter.normC, interval)
	end

	if hasCustomLabels(event)
		typeE = labels(event, event_index_C)
	end
	# record information about the event point
	userpoint = EventSpecialPoint(state, Symbol(typeE), status, iter.recordFromSolution, iter.normC, interval)
	(verbosity > 0) && printstyled(color=:red, "!! Continuous user point at p ≈ $(getp(state)) \n")
	return true, userpoint
end
####################################################################################################
function getEventType(event::AbstractDiscreteEvent, iter::AbstractContinuationIterable, state, verbosity, status::Symbol, interval::Tuple{T, T}, ind = :; typeE = "userD") where T
	event_index_D = Int32[]
	if state.eventValue[1] isa Real && (abs(state.eventValue[1] - state.eventValue[2]) > 0)
		push!(event_index_D,1)
	elseif state.eventValue[1][ind] isa Real && (abs(state.eventValue[1][ind] - state.eventValue[2][ind]) > 0)
		push!(event_index_D,1)
	else
		for ii in eachindex(state.eventValue[1][ind])
			if abs(state.eventValue[1][ind][ii] - state.eventValue[2][ind][ii]) > 0
				push!(event_index_D,ii)
				typeE = typeE * "-$ii"
			end
		end
	end
	if isempty(event_index_D) == true
		@error "Error, no event was characterized whereas one was detected. Please open an issue at https://github.com/rveltz/BifurcationKit.jl/issues. \n The events are eventValue = $(state.eventValue)"
		# we halt continuation as it will mess up the detection of events
		state.stopcontinuation = true
		return false, EventSpecialPoint(state, Symbol(typeE), status, iter.recordFromSolution, iter.normC, interval)
	end
	if hasCustomLabels(event)
		typeE = labels(event, event_index_D)
	end
	# record information about the ev point
	userpoint = EventSpecialPoint(state, Symbol(typeE), status, iter.recordFromSolution, iter.normC, interval)
	(verbosity > 0) && printstyled(color=:red, "!! Discrete user point at p ≈ $(getp(state)) \n")
	return true, userpoint
end
####################################################################################################
function getEventType(event::PairOfEvents, iter::AbstractContinuationIterable, state, verbosity, status::Symbol, interval::Tuple{T, T}) where T
	nC = length(event.eventC)
	n = length(event)

	(isEventCrossed(event.eventC, iter, state, 1:nC) && isEventCrossed(event.eventD, iter, state, nC+1:n)) && @warn "More than one Event was detected. We call the continuous event to save data in the branch."
	if isEventCrossed(event.eventC, iter, state, 1:nC)
		return getEventType(event.eventC, iter, state, verbosity, status, interval, 1:nC; typeE = "userC")
	elseif isEventCrossed(event.eventD, iter, state, nC+1:n)
		return getEventType(event.eventD, iter, state, verbosity, status, interval, nC+1:n; typeE = "userD")
	else
		@error "Error, no event was characterized whereas one was detected. Please open an issue at https://github.com/rveltz/BifurcationKit.jl/issues. \n The events are eventValue = $(state.eventValue)"
		# we halt continuation as it will mess up the detection of events
		state.stopcontinuation = true
		return false, EventSpecialPoint(state, :PairOfEvents, status, iter.recordFromSolution, iter.normC, interval)
	end
end

####################################################################################################
function getEventType(event::SetOfEvents, iter::AbstractContinuationIterable, state, verbosity, status::Symbol, interval::Tuple{T, T}) where T
	# find the active events
	event_index_C = Int32[]
	event_index_D = Int32[]
	for (ind, eve) in enumerate(event.eventC)
		if isEventCrossed(eve, iter, state, ind)
			push!(event_index_C, ind)
		end
	end

	nC = length(event.eventC)
	for (ind, eve) in enumerate(event.eventD)
		if isEventCrossed(eve, iter, state, nC + ind)
			push!(event_index_D, ind)
		end
	end

	length(event_index_C) + length(event_index_D) > 1 && @warn "More than one event in `SetOfEvents` was detected. We take the first in the list to save data in the branch."

	if isempty(event_index_C) == false
		indC = event_index_C[1]
		return getEventType(event.eventC[indC], iter, state, verbosity, status, interval, indC; typeE = "userC$indC")
	elseif isempty(event_index_D) == false
		indD = event_index_D[1]
		return getEventType(event.eventD[indD], iter, state, verbosity, status, interval, indD+nC; typeE = "userD$indD")
	else
		@error "Error, no event was characterized whereas one was detected. Please open an issue at https://github.com/rveltz/BifurcationKit.jl/issues. \n The events are eventValue = $(state.eventValue)"
		# we halt continuation as it will mess up the detection of events
		state.stopcontinuation = true
		return false, EventSpecialPoint(state, :SetOfEvents, status, iter.recordFromSolution, iter.normC, interval)
	end
end


text:
# Adapted from:
# SOSDEMO9 --- Matrix SOS decomposition
# Section 3.9 of SOSTOOLS User's Manual

@testset "SOSDEMO9 with $(factory.constructor)" for factory in sdp_factories
    @polyvar x1 x2 x3

    P = [x1^4+x1^2*x2^2+x1^2*x3^2 x1*x2*x3^2-x1^3*x2-x1*x2*(x2^2+2*x3^2);
         x1*x2*x3^2-x1^3*x2-x1*x2*(x2^2+2*x3^2) x1^2*x2^2+x2^2*x3^2+(x2^2+2*x3^2)^2]

    # Test if P(x1,x2,x3) is an SOS matrix
    m = SOSModel(factory)
    # TODO return H so that P = H.'*H
    @SDconstraint m P ⪰ 0

    JuMP.optimize!(m)
    @test JuMP.primal_status(m) == MOI.FeasiblePoint
end


text:
# https://github.com/dpsanders/hands_on_julia/blob/master/during_sessions/Fractale%20de%20Newton.ipynb
newton(f, f′, x) = x - f(x) / f′(x)
function solve(g, g′, x0, max_iterations = 100)
    x = x0

    tolerance = 2 * eps(x0)
    iteration = 0

    while (abs(g(x) - 0) > tolerance && iteration < max_iterations)
        x = newton(g, g′, x)
        iteration += 1
    end

    return x
end

abstract type InterpolationKind end

struct CubicSpline <: InterpolationKind end
struct LinearSpline <: InterpolationKind end

"""
    bootstrap(rates, maturities, settlement_frequency, interpolation::CubicSpline())

Bootstrap the rates with the given maturities, treating the rates according to the periodic frequencies in settlement_frequency.

`interpolator` is any function that will take two vectors of inputs and output points and return a function that will estimate an output given a scalar input. That is
`interpolator` should be: `interpolator(xs, ys) -> f(x)` where `f(x)` is the interpolated value of `y` at `x`.

Built in `interpolator`s in Yields are:
- `CubicSpline()`: Cubic spline interpolation.
- `LinearSpline()`: Linear spline interpolation.

The default is `CubicSpline()`.
"""
function bootstrap(rates, maturities, settlement_frequency, interpolation::InterpolationKind=CubicSpline())
    return _bootstrap_choose_interp(rates, maturities, settlement_frequency, interpolation)
end

# the fall-back if user provides own interpolation function
function bootstrap(rates, maturities, settlement_frequency, interpolation)
    return _bootstrap_inner(rates, maturities, settlement_frequency, interpolation)
end

# dispatch on the user-exposed InterpolationKind to the right
# internally named interpolation function
function _bootstrap_choose_interp(rates, maturities, settlement_frequency, i::CubicSpline)
    return _bootstrap_inner(rates, maturities, settlement_frequency, cubic_interp)
end

function _bootstrap_choose_interp(rates, maturities, settlement_frequency, i::LinearSpline)
    return _bootstrap_inner(rates, maturities, settlement_frequency, linear_interp)
end

function _bootstrap_inner(rates, maturities, settlement_frequency, interpolation_function)
    discount_vec = zeros(length(rates)) # construct a placeholder discount vector matching maturities
    # we have to take the first rate as the starting point
    discount_vec[1] = discount(Constant(rates[1]), maturities[1])

    for t = 2:length(maturities)
        if isnothing(settlement_frequency[t])
            # no settlement before maturity
            discount_vec[t] = discount(Constant(rates[t]), maturities[t])
        else
            # need to account for the interim cashflows settled
            times = settlement_frequency[t]:settlement_frequency[t]:maturities[t]
            cfs = [rate(rates[t]) * settlement_frequency[t] for s in times]
            cfs[end] += 1

            function pv(v_guess)
                v = interpolation_function([[0.0]; maturities[1:t]], vcat(1.0, discount_vec[1:t-1], v_guess...))
                return sum(v.(times) .* cfs)
            end
            target_pv = sum(map(t2 -> discount(Constant(rates[t]), t2), times) .* cfs)
            root_func(v_guess) = pv(v_guess) - target_pv
            root_func′(v_guess) = ForwardDiff.derivative(root_func, v_guess)
            discount_vec[t] = solve(root_func, root_func′, rate(rates[t]))
        end

    end
    zero_vec = -log.(discount_vec) ./ maturities
    return linear_interp([0.0; maturities], [first(zero_vec); zero_vec])
end

# the ad-hoc approach to extrapoliatons is based on suggestion by author of
# BSplineKit at https://github.com/jipolanco/BSplineKit.jl/issues/19
# this should not be exposed directly to user
struct _Extrap{I,L,R}
	int::I # the BSplineKit interpolation
	left::L # a tuple of (boundary, extrapolation function)
	right::R # a tuple of (boundary, extrapolation function)
end

function _wrap_spline(itp)

	S = BSplineKit.spline(itp)  # spline passing through data points
	B = BSplineKit.basis(S)     # B-spline basis

	a, b = BSplineKit.boundaries(B)  # left and right boundaries

	# For now, we construct the full spline S′(x).
	# There are faster ways of doing this that should be implemented...
	S′ = diff(S, BSplineKit.Derivative(1))

	return _Extrap(itp,
		(boundary = a, func = x->S(a) + S′(a)*(x-a)),
		(boundary = b, func = x->S(b) + S′(b)*(x-b)),

	)
end

function _interp(e::_Extrap,x)
	if x <= e.left.boundary
		return e.left.func(x)
	elseif x >= e.right.boundary
		return e.right.func(x)
	else
		return e.int(x)
	end
end

function linear_interp(xs, ys)
    int = BSplineKit.interpolate(xs, ys, BSplineKit.BSplineOrder(2))
    e = _wrap_spline(int)
    return x -> _interp(e, x)
end

function cubic_interp(xs, ys)
    int = BSplineKit.interpolate(xs, ys, BSplineKit.BSplineOrder(2))
    e = _wrap_spline(int)
    return x -> _interp(e, x)
end

# used to display simple type name in show method
# https://stackoverflow.com/questions/70043313/get-simple-name-of-type-in-julia?noredirect=1#comment123823820_70043313
name(::Type{T}) where {T} = (isempty(T.parameters) ? T : T.name.wrapper)

function Base.show(io::IO, curve::T) where {T<:AbstractYield}
    println() # blank line for padding
    r = zero(curve, 1)
    ylabel = isa(r.compounding, Continuous) ? "Continuous" : "Periodic($(r.compounding.frequency))"
    kind = name(typeof(curve))
    l = lineplot(
        t -> rate(zero(curve, t)),
        0.0, #from
        30.0,  # to
        xlabel = "time",
        ylabel = ylabel,
        compact = true,
        name = "Zero rates",
        width = 60,
        title = "Yield Curve ($kind)"
    )
    show(io, l)
end