text:
export empty_path,
       open_path,
       closed_path,
       open_path_ops,
       closed_path_ops,
       MoveOp,
       MoveToOp,
       LineOp,
       ArcOp,
       CloseOp,
       ArcPath,
       arc_path,
       CircularPath,
       circular_path,
       RectangularPath,
       rectangular_path,
       centered_rectangular_path,
       OpenPolygonalPath,
       open_polygonal_path,
       ClosedPolygonalPath,
       closed_polygonal_path,
       PolygonalPath,
       polygonal_path,
       OpenSplinePath,
       open_spline_path,
       ClosedSplinePath,
       closed_spline_path,
       SplinePath,
       spline_path,
       PathSet,
       path_set,
       OpenPathSequence,
       open_path_sequence,
       ClosedPathSequence,
       closed_path_sequence,
       PathSequence,
       path_sequence,
       translate,
       stroke,
       fill,
       path_length,
       location_at,
       location_at_length,
       path_start,
       path_end,
       path_domain,
       planar_path_normal,
       subpath,
       subpaths,
       join_paths,
       subtract_paths,
       path_vertices,
       path_frames,
       mirrored_path,
       mirrored_on_x,
       mirrored_on_y,
       mirrored_on_z,
       ClosedPath,
       OpenPath,
       Region,
       region,
       outer_path,
       inner_paths,
       path_tolerance

path_tolerance = Parameter(1e-10)
coincident_path_location(p1::Loc, p2::Loc) = distance(p1, p2) < path_tolerance()

abstract type Path end

struct EmptyPath <: Path
end

empty_path() = EmptyPath()
is_empty_path(path::Path) = false
is_empty_path(path::EmptyPath) = true

struct PointPath <: Path
    location::Loc
end

getindex(p::Path, i::Real) = location_at_length(p, i)
firstindex(p::Path) = 0
lastindex(p::Path) = path_length(p)
#getindex(p::Path, i::ClosedInterval) = subpath(p, i.left, i.right)
path_start(p::Path) = location_at_length(p, 0)
path_end(p::Path) = location_at_length(p, path_length(p))

abstract type OpenPath <: Path end
abstract type ClosedPath <: Path end

is_closed_path(path::Path) = false
is_closed_path(path::ClosedPath) = true
is_open_path(path::Path) = ! is_closed_path(path)
# To avoid duplicate convertions, we also deal with Locs
is_closed_path(pts::Locs) = coincident_path_location(pts[1], pts[end])

# Default path domain, probably only makes sense for Splines.
path_domain(path::Path) = (0, 1)

# Generic map_division using path default domain divided into equal parts
map_division(f::Function, path::Path, n::Integer) =
  map_division(ϕ->f(location_at(path, ϕ)), path_domain(path)..., n, is_open_path(path))
# Even more generic map_division, with default domain and divided "naturaly"
map_division(f::Function, path::Path) =
  f.(path_interpolated_frames(path))

path_interpolated_frames(path::Path, (t0, t1)=path_domain(path), epsilon=collinearity_tolerance(), min_recursion=1) =
  let p0 = location_at(path, t0),
      p1 = location_at(path, t1),
      tm = (t0 + t1)/2.0,
      pm = location_at(path, tm)
    min_recursion < 0 &&
    collinear_points(p0, pm, p1, epsilon) &&
    let tr = (t0 + t1)/(1.99+rand()*0.02),
        pr = location_at(path, tr) # To avoid coincidences
      collinear_points(p0, pr, p1, epsilon)
    end ?
      [p0, pm, p1] :
      [path_interpolated_frames(path, (t0, tm), epsilon, min_recursion - 1)...,
       path_interpolated_frames(path, (tm, t1), epsilon, min_recursion - 1)[2:end]...]
  end

## Arc path
struct ArcPath <: OpenPath
    center::Loc
    radius::Real
    start_angle::Real
    amplitude::Real
end
arc_path(center::Loc=u0(), radius::Real=1, start_angle::Real=0, amplitude::Real=pi*1.0) =
  # pi*1.0 just to circunvent a Julia bug https://github.com/JuliaLang/julia/issues/31949.
  false ? #amplitude < 0 ?
    ArcPath(center, radius, start_angle + amplitude, - amplitude) :
    ArcPath(center, radius, start_angle, amplitude)
path_domain(path::ArcPath) = (0, path.amplitude)
location_at(path::ArcPath, ϕ::Real) =
  let s = sign(path.amplitude*1.0), # just to circunvent a Julia bug https://github.com/JuliaLang/julia/issues/31949.
      ϕ = ϕ*s
    loc_from_o_vx_vy(add_pol(path.center, path.radius, path.start_angle + ϕ),
                     vpol(1, path.start_angle + ϕ + π, path.center.cs),
                     vz(1, path.center.cs))
  end
planar_path_normal(path::ArcPath) = uvz(path.center.cs)

## Circular path
struct CircularPath <: ClosedPath
    center::Loc
    radius::Real
end
circular_path(Center::Loc=u0(), Radius::Real=1; center::Loc=Center, radius::Real=Radius) =
  CircularPath(center, radius)
path_domain(path::CircularPath) = (0, 2π)
location_at(path::CircularPath, ϕ::Real) =
  loc_from_o_vx_vy(add_pol(path.center, path.radius, ϕ),
                   vpol(1, ϕ + π, path.center.cs),
                   vz(1, path.center.cs))
planar_path_normal(path::CircularPath) = uvz(path.center.cs)

## Rectangular path
struct RectangularPath <: ClosedPath
    corner::Loc
    dx::Real
    dy::Real
end
rectangular_path(corner::Loc=u0(), dx::Real=1, dy::Real=1) =
  RectangularPath(corner, dx, dy)
centered_rectangular_path(p, dx, dy) =
  rectangular_path(p-vxy(dx/2, dy/2), dx, dy)
path_domain(path::RectangularPath) = (0, path_length(path))
location_at(path::RectangularPath, d::Real) = location_at_length(path, d)
planar_path_normal(path::RectangularPath) = uvz(path.corner.cs)

struct OpenPolygonalPath <: OpenPath
    vertices::Locs
end
open_polygonal_path(vertices=[u0(), x(), xy(), y()]) =
  OpenPolygonalPath(vertices)

struct ClosedPolygonalPath <: ClosedPath
    vertices::Locs
end
closed_polygonal_path(vertices=[u0(), x(), xy(), y()]) =
  ClosedPolygonalPath(ensure_no_repeated_locations(vertices))

PolygonalPath = Union{OpenPolygonalPath, ClosedPolygonalPath}

polygonal_path(vertices::Locs=[u0(), x(), xy(), y(), u0()]) =
  coincident_path_location(vertices[1], vertices[end]) ?
    closed_polygonal_path(vertices[1:end-1]) :
    open_polygonal_path(vertices)
polygonal_path(v::Loc, vs...) = polygonal_path([v, vs...])

ensure_no_repeated_locations(locs) =
    begin
        @assert ! coincident_path_location(locs[1], locs[end])
        locs
    end

path_start(path::PolygonalPath) = path.vertices[1]
path_end(path::OpenPolygonalPath) = path.vertices[end]
path_end(path::ClosedPath) = path_start(path)
path_domain(path::PolygonalPath) = (0, path_length(path))
location_at(path::PolygonalPath, ϕ::Real) = location_at_length(path, ϕ)
planar_path_normal(path::PolygonalPath) = vertices_normal(path_vertices(path))

join_paths(p1::OpenPolygonalPath, p2::OpenPolygonalPath) =
    coincident_path_location(p1.vertices[end], p2.vertices[1]) ?
      collinear_points(p1.vertices[end-1], p1.vertices[end], p2.vertices[2]) ?
        polygonal_path([p1.vertices[1:end-1]..., p2.vertices[2:end]...]) :
        polygonal_path([p1.vertices[1:end-1]..., p2.vertices...]) :
      polygonal_path([p1.vertices..., p2.vertices...])

# Splines

struct OpenSplinePath <: OpenPath
    vertices::Locs
    v0::Union{Bool,Vec}
    v1::Union{Bool,Vec}
    interpolator
end
open_spline_path(vertices=[u0(), x(), xy(), y()], v0=false, v1=false) =
  OpenSplinePath(vertices, v0, v1, curve_interpolator(vertices, false))

struct ClosedSplinePath <: ClosedPath
    vertices::Locs
    interpolator
end
closed_spline_path(vertices=[u0(), x(), xy(), y()]) =
  let vertices = ensure_no_repeated_locations(vertices)
    ClosedSplinePath(vertices, curve_interpolator(vertices, true))
  end

spline_path(vertices::Locs) =
  coincident_path_location(vertices[1], vertices[end]) ?
    closed_spline_path(vertices[1:end-1]) :
    open_spline_path(vertices)
spline_path(v::Loc, vs...) = spline_path([v, vs...])

location_at(path::Union{OpenSplinePath,ClosedSplinePath}, ϕ::Real) =
  location_at(path.interpolator, ϕ)
location_at(interpolator::ParametricSpline, ϕ) =
  let f = interpolator(ϕ),
      d1 = derivative(interpolator, ϕ), #Interpolations.gradient(interpolator, ϕ)[1],
      d2 = derivative(interpolator, ϕ, nu=2),#u0 = path.vertices[1],
      p = xyz(f[1], f[2], f[3], world_cs),
      t = vxyz(d1[1], d1[2], d1[3], world_cs),
      n = vxyz(d2[1], d2[2], d2[3], world_cs)#,
    norm(n) <= 1e-9 ? loc_from_o_vz(p, t) : loc_from_o_vx_vy(p, n, cross(t, n))
  end

map_division(func::Function, path::OpenSplinePath, n::Integer) =
  let interpol = path.interpolator,
      fs = map_division(interpol, 0.0, 1.0, n, true),
      f = fs[1],
      d1s = map_division(t->derivative(interpol, t), 0.0, 1.0, n, true),
      d1 = d1s[1],
      d2 = derivative(interpol, 0.0, nu=2),
      p = xyz(f[1], f[2], f[3], world_cs),
      t = vxyz(d1[1], d1[2], d1[3], world_cs),
      n = vxyz(d2[1], d2[2], d2[3], world_cs)
    func.(rotation_minimizing_frames(
            norm(n) <= 1e-9 ? loc_from_o_vz(p, t) : loc_from_o_vx_vy(p, n, cross(t, n)),
            map(p->xyz(p[1], p[2], p[3], world_cs), fs),
            map(t->vxyz(t[1], t[2], t[3], world_cs), d1s)))
  end

map_division(f::Function, path::ClosedSplinePath, n::Integer) =
  let interpolator = path.interpolator,
      ps = map_division(interpolator, 0.0, 1.0, n, false),
      vts = map_division(t->derivative(interpolator, t), 0.0, 1.0, n, false)
      #vns = map_division(t->Interpolations.hessian(interpolator, t)[1], 0.0, 1.0, n)
    f.(rotation_minimizing_frames(
         path.vertices[1],
         map(p->xyz(p[1], p[2], p[3], world_cs), ps),
         map(vy->vxyz(vy[1], vy[2], vy[3], world_cs), vts)))
  end

SplinePath = Union{OpenSplinePath, ClosedSplinePath}

path_start(path::SplinePath) = path.vertices[1]
path_end(path::OpenSplinePath) = path.vertices[end]
path_end(path::ClosedSplinePath) = path.vertices[1]

convert(::Type{ClosedSplinePath}, path::Path) =
  closed_spline_path(path_vertices(path))

#=
      fixed_normal(vn, vt) = norm(vn) < path_tolerance() ? SVector{3}(vpol(1, sph_phi(xyz(vt[1],vt[2],vt[3], world_cs))+pi/2).raw[1:3]) : vn
    map_division(
      t-> let p = interpolator(t),
              vt = Interpolations.gradient(interpolator, t)[1],
              vn = Interpolations.hessian(interpolator, t)[1],
              vy = cross(vt, vn)
            f(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,
     0.0, 1.0, n)
 end
=#

# There is a set of operations over Paths:
# 1. translate a path a given vector
# 2. stroke a path
# 3. fill a (presumably closed) path
# 4. compute a path location given a length from the path beginning
# 5. compute a sub path from a path, a length from the path begining and a length increment
# 6. Produce the meta representation of a path

translate(path::PointPath, v::Vec) = PointPath(path.location + v)
translate(path::CircularPath, v::Vec) = CircularPath(path.center + v, path.radius)
translate(path::ArcPath, v::Vec) = ArcPath(path.center + v, path.radius, path.start_angle, path.amplitude)
translate(path::RectangularPath, v::Vec) = RectangularPath(path.corner + v, path.dx, path.dy)
translate(path::OpenPolygonalPath, v::Vec) = OpenPolygonalPath(translate(path.vertices, v))
translate(path::ClosedPolygonalPath, v::Vec) = ClosedPolygonalPath(translate(path.vertices, v))
translate(path::OpenSplinePath, v::Vec) = OpenSplinePath(translate(path.vertices, v), path.v0, path.v1, path.interpolator)
translate(path::ClosedSplinePath, v::Vec) = ClosedSplinePath(translate(path.vertices, v), path.interpolator)
translate(ps::Locs, v::Vec) = map(p->p+v, ps)

in_cs(path::PointPath, cs::CS) = PointPath(in_cs(path.location, cs))
in_cs(path::CircularPath, cs::CS) = CircularPath(in_cs(path.center, cs), path.radius)
in_cs(path::ArcPath, cs::CS) = ArcPath(in_cs(path.center, cs), path.radius, path.start_angle, path.amplitude)
in_cs(path::RectangularPath, cs::CS) = RectangularPath(in_cs(path.corner, cs), path.dx, path.dy)
in_cs(path::OpenPolygonalPath, cs::CS) = OpenPolygonalPath(in_cs(path.vertices, cs))
in_cs(path::ClosedPolygonalPath, cs::CS) = ClosedPolygonalPath(in_cs(path.vertices, cs))
in_cs(path::OpenSplinePath, cs::CS) = OpenSplinePath(in_cs(path.vertices, cs), path.v0, path.v1)
in_cs(path::ClosedSplinePath, cs::CS) = ClosedSplinePath(in_cs(path.vertices, cs))
in_cs(ps::Locs, cs::CS) = map(p->in_cs(p, cs), ps)

scale(path::CircularPath, s::Real, p::Loc=u0()) =
  let v = path.center - p
    circular_path(p + v*s, path.radius*s)
  end
scale(path::RectangularPath, s::Real, p::Loc=u0()) =
  let v = path.corner - p
    rectangular_path(p + v*s, path.dx*s, path.dy*s)
  end
scale(path::OpenPolygonalPath, s::Real, p::Loc=u0()) =
  open_polygonal_path([q + (q-p)*s for q in path.vertices])
scale(path::ClosedPolygonalPath, s::Real, p::Loc=u0()) =
  closed_polygonal_path([q + (q-p)*s for q in path.vertices])

path_length(path::CircularPath) = 2*pi*path.radius
path_length(path::ArcPath) = path.radius*abs(path.amplitude)
path_length(path::RectangularPath) = 2*(path.dx + path.dy)
path_length(path::OpenPolygonalPath) = path_length(path.vertices)
path_length(path::ClosedPolygonalPath) = path_length(path.vertices) + distance(path.vertices[end], path.vertices[1])
path_length(ps::Locs) =
  let p = ps[1]
      l = 0.0
    for i in 2:length(ps)
        pp = ps[i]
        l += distance(p, pp)
        p = pp
    end
    l
  end

location_at_length(path::CircularPath, d::Real) =
  location_at(path, d/path.radius)
location_at_length(path::ArcPath, d::Real) =
  let Δα = d/path.radius,
      s = sign(path.amplitude)
    Δα <= abs(path.amplitude) + path_tolerance() ?
      location_at(path, Δα*s) :
      error("Exceeded path length by ", Δα - path.amplitude)
  end
location_at_length(path::RectangularPath, d::Real) =
    let d = d % (2*(path.dx + path.dy)) # remove multiple periods
        p = path.corner
        for (delta, phi) in zip([path.dx, path.dy, path.dx, path.dy], [0, pi/2, pi, 3pi/2])
            if d - delta < path_tolerance()
                return loc_from_o_phi(add_pol(p, d, phi), phi)
            else
                p = add_pol(p, delta, phi)
                d -= delta
            end
        end
    end
location_at_length(path::OpenPolygonalPath, d::Real) =
  let p = path.vertices[1]
    for i in 2:length(path.vertices)
      pp = path.vertices[i]
      delta = distance(p, pp)
      if d - delta < path_tolerance()
        v = unitized(pp - p)
        return loc_from_o_vz(p+v*d, v)
      else
        p = pp
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end
location_at_length(path::ClosedPolygonalPath, d::Real) =
  let p = path.vertices[1]
    for i in countfrom(1)
      pp = path.vertices[i%length(path.vertices)+1]
      delta = distance(p, pp)
      if d - delta < path_tolerance()
        v = unitized(pp - p)
        return loc_from_o_vz(p+v*d, v)
      else
        p = pp
        d -= delta
      end
    end
  end

subpath_starting_at(path::Path, d::Real) = subpath(path, d, path_length(path))
subpath_ending_at(path::Path, d::Real) = subpath(path, 0, d)

subpath(path::CircularPath, a::Real, b::Real) =
  arc_path(path.center, path.radius, a/path.radius, (b-a)/path.radius)
subpath(path::ArcPath, a::Real, b::Real) =
  b <= path_length(path) + path_tolerance() ?
      arc_path(path.center,
               path.radius,
               path.start_angle + a/path.radius*sign(path.amplitude),
               (b - a)/path.radius*sign(path.amplitude)) :
    error("Exceeded path length by ", path.amplitude - b/path.radius)
subpath(path::RectangularPath, a::Real, b::Real) =
  subpath(convert(ClosedPolygonalPath, path), a, b)
subpath(path::ClosedPolygonalPath, a::Real, b::Real) =
  subpath(convert(OpenPolygonalPath, path), a, b)
subpath(path::OpenPolygonalPath, a::Real, b::Real) =
  subpath_starting_at(subpath_ending_at(path, b), a)

subpath_starting_at(path::OpenPolygonalPath, d::Real) =
    let pts = path.vertices
        p1 = pts[1]
        for i in 2:length(pts)
            p2 = pts[i]
            delta = distance(p1, p2)
            diff = d - delta
            if diff < 0
                mp = p1 + (p2 - p1)*d/delta
                return open_polygonal_path([mp, pts[i:end]...])
            elseif abs(diff) < 1e-9
                return open_polygonal_path(i == length(pts) ? [pts[i], pts[i]] : pts[i:end])
            else
                p1 = p2
                d -= delta
            end
        end
        abs(d) < path_tolerance() ?
          path :
          error("Exceeded path length by ", d)
    end

subpath_ending_at(path::OpenPolygonalPath, d::Real) =
    let pts = path.vertices
        p1 = pts[1]
        for i in 2:length(pts)
            p2 = pts[i]
            delta = distance(p1, p2)
            diff = d - delta
            if diff < 0
                mp = p1 + (p2 - p1)*d/delta
                return open_polygonal_path([pts[1:i-1]..., mp])
            elseif abs(diff) < path_tolerance()
                return open_polygonal_path(pts[1:i])
            else
              p1 = p2
              d -= delta
            end
        end
        abs(d) < path_tolerance() ?
          path :
          error("Exceeded path length by ", d)
    end

#=
collect_vertices_length(p::Loc, vs::Locs, d::Real) =
    let pp = vs[1]
        dist = distance(p, pp)
        d <= dist ?
        push!(Vector{Loc}(), p + (pp - p)/dist*d) :
        unshift!(collect_vertices_length(pp, vs[2:end], d - dist), pp)
    end

subpath(path::OpenPolygonalPath, a::Real, b::Real) =
    let p = path.vertices[1]
        pts = []
        for i in 2:length(path.vertices)
            pp = path.vertices[i]
            delta = distance(p, pp)
            if a <= delta
                phi = pol_phi(pp - p)
                p0 = add_pol(p, d, phi)
                return open_polygonal_path(
                        unshift!(collect_vertices_length(p0,
                                                         vcat(path.vertices[i:end], path.vertices),
                                                         delta_d),
                                 p0))
            else
                p = pp
                d -= delta
            end
        end
        error("Exceeded path length")
    end

=#

meta_program(p::OpenPolygonalPath) =
    Expr(:call, :open_polygonal_path, meta_program(p.vertices))

# Path can be made of subparts
abstract type PathOp end
#struct MoveToOp <: PathOp loc::Loc end
#struct MoveOp <: PathOp vec::Vec end
struct LineOp <: PathOp
  vec::Vec
end
#struct CloseOp <: PathOp end
struct ArcOp <: PathOp
    radius::Real
    start_angle::Real
    amplitude::Real
end
struct LineXThenYOp <: PathOp
  vec::Vec
end
struct LineYThenXOp <: PathOp
  vec::Vec
end

struct PathOps <: Path
  start::Loc
  ops::Vector{<:PathOp}
  closed::Bool
end
open_path_ops(start, ops...) = PathOps(start, [ops...], false)
closed_path_ops(start, ops...) = PathOps(start, [ops...], true)

path_length(path::PathOps) =
  let len = mapreduce(path_length, +, path.ops, init=0)
    path.closed ?
      len + distance(path.start, location_at_length(path, len)) :
      len
  end
path_length(op::LineOp) = norm(op.vec)
path_length(op::ArcOp) = op.radius*abs(op.amplitude)

translate(path::PathOps, v::Vec) =
  PathOps(path.start + v, path.ops, path.closed)

#translate_op(op::MoveToOp, v) = MoveToOp(op.loc + v)
#translate_op(op::LineToOp, v) = LineToOp(op.loc + v)
#translate_op(op::LineToXThenToYOp, v) = LineToXThenToYOpp(op.loc + v)
#translate_op(op::LineToYThenToXOp, v) = LineToYThenToXOp(op.loc + v)
#translate_op(op::PathOp, v) = op

location_at_length(path::PathOps, d::Real) =
  let ops = path.ops,
      start = path.start
    for op in ops
      delta = path_length(op)
      if d < delta + path_tolerance()
        return location_at_length(op, start, d)
      else
        start = location_at_length(op, start, delta)
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end

location_at_length(op::LineOp, start::Loc, d::Real) =
  start + op.vec*d/path_length(op)
location_at_length(op::ArcOp, start::Loc, d::Real) =
  let center = start - vpol(op.radius, op.start_angle),
      a = d/op.radius*sign(op.amplitude)
    center + vpol(op.radius, op.start_angle + a)
  end

subpath_starting_at(path::PathOps, d::Real) =
  let ops = path.ops,
      start = location_at_length(path, d)
    for i in 1:length(ops)
      op = ops[i]
      delta = path_length(op)
      if d == delta
        return PathOps(start, ops[i+1:end], false)
      elseif d < delta + path_tolerance()
        op = subpath_starting_at(op, d)
        return PathOps(start, [op, ops[i+1:end]...], false)
      else
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end
subpath_ending_at(path::PathOps, d::Real) =
  let ops = path.ops
    for i in 1:length(ops)
      op = ops[i]
      delta = path_length(op)
      if d == delta
        return PathOps(path.start, ops[1:i], false)
      elseif d < delta + path_tolerance()
        return PathOps(path.start, [ops[1:i-1]..., subpath_ending_at(op, d)], false)
      else
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end

subpath_starting_at(pathOp::LineOp, d::Real) =
  let len = path_length(pathOp.vec)
    LineOp(pathOp.vec*(len-d)/len)
  end
subpath_starting_at(pathOp::ArcOp, d::Real) =
  let a = d/pathOp.radius
    ArcOp(pathOp.radius, pathOp.start_angle + a, pathOp.amplitude - a)
  end

subpath_ending_at(pathOp::LineOp, d::Real) =
  let len = path_length(pathOp.vec)
    LineOp(pathOp.vec*d/len)
  end
subpath_ending_at(pathOp::ArcOp, d::Real) =
  let a = d/pathOp.radius
    ArcOp(pathOp.radius, pathOp.start_angle, a)
  end

subpath(path::Path, a::Real, b::Real) =
  subpath_starting_at(subpath_ending_at(path, b), a)

# A path sequence is a sequence of paths where the next element of the sequence
# starts at the same place where the previous element ends.

struct OpenPathSequence <: OpenPath
  paths::Vector{<:Path}
end
open_path_sequence(paths...) =
  OpenPathSequence(ensure_connected_paths([paths...]))
struct ClosedPathSequence <: ClosedPath
  paths::Vector{<:Path}
end
closed_path_sequence(paths...) =
  ClosedPathSequence(ensure_connected_paths([paths..., polygonal_path(path_end(paths[end]), path_start(paths[1]))]))

PathSequence = Union{OpenPathSequence, ClosedPathSequence}

# HACK : Include treatment of empty paths.

path_sequence(paths...) =
  let paths = ensure_connected_paths([paths...])
    coincident_path_location(path_start(paths[1]), path_end(paths[end])) ?
      ClosedPathSequence(paths) :
      OpenPathSequence(paths)
  end

ensure_connected_paths(paths) = filter(!is_empty_path, paths)

path_length(path::PathSequence) =
  sum(map(path_length, path.paths))
location_at_length(path::PathSequence, d::Real) =
  let paths = path.paths
    for path in paths
      delta = path_length(path)
      if d <= delta + path_tolerance()
          return location_at_length(path, d)
      else
          d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end
subpath_starting_at(path::PathSequence, d::Real) =
  let subpaths = path.paths
    for i in 1:length(subpaths)
      subpath = subpaths[i]
      delta = path_length(subpath)
      if d == delta
        return OpenPathSequence(ops[i+1:end])
      elseif d < delta + path_tolerance()
        subpath = subpath_starting_at(subpath, d)
        return OpenPathSequence([subpath, subpaths[i+1:end]...])
      else
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end
subpath_ending_at(path::PathSequence, d::Real) =
  let subpaths = path.paths
    for i in 1:length(subpaths)
      subpath = subpaths[i]
      delta = path_length(subpath)
      if d == delta
        return OpenPathSequence(subpaths[1:i])
      elseif d < delta + path_tolerance()
        return OpenPathSequence([subpaths[1:i-1]..., subpath_ending_at(subpath, d)])
      else
        d -= delta
      end
    end
    error("Exceeded path length by ", d)
  end

convert(::Type{PathOps}, path::PathSequence) =
  let start = in_world(location_at_length(path, 0)),
      ops = mapfoldl(convert_to_path_ops, vcat, path.paths, init=PathOp[])
    PathOps(start, ops, coincident_path_location(start, path_end(path)))
  end
convert_to_path_ops(path::ArcPath) =
  [ArcOp(path.radius, path.start_angle, path.amplitude)]
convert_to_path_ops(path::PolygonalPath) =
  [LineOp(in_world(p)-in_world(q)) for (p,q) in zip(path.vertices[[2:end;1]], path.vertices)]
convert_to_path_ops(path::RectangularPath) =
  let p0 = in_world(path.corner),
      p1 = in_world(add_x(path.corner, path.dx)),
      p2 = in_world(add_xy(path.corner, path.dx, path.dy)),
      p3 = in_world(add_y(path.corner, path.dy))
    [LineOp(p1-p0), LineOp(p2-p1), LineOp(p3-p2), LineOp(p0-p3)]
  end

# A path set is a set of independent paths.
struct PathSet <: Path
  paths::Vector{<:Path}
end

# Should we just use tuples instead of arrays?
path_set(paths...) =
  PathSet([paths...])

# Regions are areas delimited by paths. They are assumed to be planar and might contain holes.
struct Region
  paths::Vector{<:ClosedPath}
end

region(outer, inners...) =
  Region([convert(ClosedPath, outer), convert.(ClosedPath, inners)...])
outer_path(region::Region) = region.paths[1]
inner_paths(region::Region) = region.paths[2:end]

region(r::Region, inners...) =
  region(outer_path(r), inner_paths(r)..., inners...)

# Convertions from/to paths
convert(::Type{Region}, vs::Locs) =
  region(closed_polygonal_path(vs))
convert(::Type{Region}, path::ClosedPath) =
  region(path)

# Operations on path containers
translate(path::T, v::Vec) where T<:Union{PathSequence,PathSet,Region} =
  T(translate.(path.paths, v))
in_cs(path::T, cs::CS) where T<:Union{PathSequence,PathSet,Region} =
  T(map(p->in_cs(p, cs), path.paths))
planar_path_normal(path::T) where T<:Union{PathSequence,PathSet,Region} =
  planar_path_normal(path.paths[1])

# Convertions from/to paths
convert(::Type{OpenPath}, vs::Locs) =
  open_polygonal_path(vs)
convert(::Type{ClosedPath}, vs::Locs) =
  closed_polygonal_path(vs)
convert(::Type{Vector{<:ClosedPath}}, ps::Vector{<:Any}) =
  ClosedPath[convert(ClosedPath, p) for p in ps]

convert(::Type{Path}, vs::Locs) =
  coincident_path_location(vs[1], vs[end]) ?
    closed_polygonal_path(vs[1:end-1]) :
    open_polygonal_path(vs)
convert(::Type{ClosedPath}, p::OpenPolygonalPath) =
  closed_polygonal_path(coincident_path_location(path_start(p), path_end(p)) ?
    path_vertices(p)[1:end-1] :
    path_vertices(p))
convert(::Type{ClosedPath}, p::OpenPath) =
  if isa(p.ops[end], CloseOp) || coincident_path_location(path_start(p), path_end(p))
    closed_path(p.ops)
  else
    error("Can't convert to a Closed Path: $p")
  end
convert(::Type{ClosedPath}, ops::Vector{<:PathOp}) =
  if isa(ops[end], CloseOp)
    closed_path(ops)
  else
    error("Can't convert to a Closed Path: $ops")
  end
convert(::Type{OpenPath}, ops::Vector{<:PathOp}) = open_path(ops)
convert(::Type{Path}, ops::Vector{<:PathOp}) =
  if isa(ops[end], CloseOp)
    closed_path(ops)
  else
    open_path(ops)
  end
convert(::Type{ClosedPolygonalPath}, path::RectangularPath) =
  let p = path.corner,
      dx = path.dx,
      dy = path.dy
    closed_polygonal_path([p, add_x(p, dx), add_xy(p, dx, dy), add_y(p, dy)])
  end

path_interpolated_lengths(path, t0=0, t1=path_length(path), epsilon=collinearity_tolerance(), min_recursion=1) =
  let p0 = location_at_length(path, t0),
      p1 = location_at_length(path, t1),
      tm = (t0 + t1)/2.0,
      pm = location_at_length(path, tm)
    min_recursion < 0 && collinear_points(p0, pm, p1, epsilon) ?
      [t0, tm, t1] :
      [path_interpolated_lengths(path, t0, tm, epsilon, min_recursion - 1)...,
       path_interpolated_lengths(path, tm, t1, epsilon, min_recursion - 1)[2:end]...]
  end

convert(::Type{ClosedPolygonalPath}, path::CircularPath) =
  closed_polygonal_path(path_interpolated_frames(path)[1:end-1])
convert(::Type{OpenPolygonalPath}, path::CircularPath) =
  open_polygonal_path(path_interpolated_frames(path))
convert(::Type{OpenPolygonalPath}, path::ArcPath) =
  open_polygonal_path(path_interpolated_frames(path))
convert(::Type{OpenPolygonalPath}, path::ClosedPolygonalPath) =
  open_polygonal_path(vcat(path.vertices, [path.vertices[1]]))
convert(::Type{OpenPolygonalPath}, path::RectangularPath) =
  convert(OpenPolygonalPath, convert(ClosedPolygonalPath, path))
convert(::Type{OpenPolygonalPath}, path::OpenSplinePath) =
  open_polygonal_path(map_division(identity, path))

curve_interpolator(pts::Locs, closed::Bool) =
  ParametricSpline(
    [pt.raw[i] for i in 1:3, pt in in_world.(closed ? [pts..., pts[1]] : pts)],
    k=length(pts)<=3 ? 2 : 3,
    periodic=closed)
curve_control_points(interpolator) =
  let cps = get_coeffs(interpolator)
    [xyz(cps[1,j], cps[2,j], cps[3,j], world_cs)
     for j in 1:size(cps, 2)]
  end
curve_knots(interpolator) =
  get_knots(interpolator)

convert(::Type{ClosedPolygonalPath}, path::ClosedPathSequence) =
  let paths = path.paths,
      vertices = []
    for path in paths
      append!(vertices, convert(OpenPolygonalPath, path).vertices[1:end-1])
    end
    closed_polygonal_path(vertices)
  end

convert(::Type{OpenPolygonalPath}, path::OpenPathSequence) =
  let paths = path.paths,
      vertices = []
    for path in paths[1:end-1]
      append!(vertices, convert(OpenPolygonalPath, path).vertices[1:end-1])
    end
    append!(vertices, convert(OpenPolygonalPath, paths[end]).vertices)
    open_polygonal_path(vertices)
  end

# It is possible to convert a PathSet to a singleton path
# by considering the first path as the outer path and all the
# others as inner paths
convert(::Type{ClosedPath}, pset::PathSet) =
  foldl(subtract_paths, pset.paths)

#### Utilities
path_vertices(path::Union{PolygonalPath,SplinePath}) =
  path.vertices
path_vertices(path::Path) =
  path_vertices(convert(is_closed_path(path) ? ClosedPolygonalPath : OpenPolygonalPath,
                        path))

iterate(path::PolygonalPath, i=1) =
  i > length(path.vertices) ?
    nothing :
    (path.vertices[i], i+1)

loc_from_o_vx_vz(o, vx, vz) =
  norm(vx) < min_norm ?
    loc_from_o_vz(o, vz) :
    let vy = -cross(vx, vz)
      loc_from_o_vx_vy(o, cross(vy, vz), vy)
    end

loc_from_o_px_pz(o, px, pz) =
  loc_from_o_vx_vz(o, px - o, pz - o)

path_frames(path::Path) =
  path_interpolated_frames(path)

path_frames(path::RectangularPath) =
  path_frames(convert(ClosedPolygonalPath, path))

average_frame(f1, f2) =
  let vx1 = in_world(vx(1, f1.cs)),
      vy1 = in_world(vy(1, f1.cs)),
      vx2 = in_world(vx(1, f2.cs)),
      vy2 = in_world(vy(1, f2.cs))
    loc_from_o_vx_vy(intermediate_loc(f1, f2), (vx1 + vx2)/2, (vy1 + vy2)/2)
  end

path_frames(path::OpenPolygonalPath) =
  let pts = path.vertices,
      pt1 = pts[1],
      pt2 = pts[2],
      ptn1 = pts[end-1],
      ptn = pts[end]
    length(pts) < 3 ?
      loc_from_o_vz(pt1, pt1-pt2) :
      let normal(pprev, p, pnext) = (unitized(pnext-p) - unitized(p-pprev))
        [loc_from_o_vx_vz(pt1, normal(pt1, pt2, pts[3]), pt2-pt1),
         [loc_from_o_vx_vz(p, normal(pprev, p, pnext), ((p-pprev) + (pnext-p))/2)
          for (pprev, p, pnext) in zip(pts, drop(pts, 1), drop(pts, 2))]...,
         loc_from_o_vx_vz(ptn, normal(pts[end-2], ptn1, ptn), ptn - ptn1)]
      end
  end

path_frames(path::ClosedPolygonalPath) =
  let pts1 = path.vertices,
      pts2 = drop(cycle(pts1), 1),
      pts3 = drop(pts2, 1)
    [loc_from_o_py_pz(p0, py, pz) for (p0, py, pz) in zip(pts2, pts1, pts3)]
  end

subpaths(path::OpenPolygonalPath) =
  let ps = path.vertices
    map((p0, p1)->open_polygonal_path([p0, p1]), ps[1:end-1], ps[2:end])
  end
subpaths(path::ClosedPolygonalPath) =
  let ps = path.vertices
      map((p0, p1)->open_polygonal_path([p0, p1]), ps, [ps[2:end]..., ps[1]])
  end
subpaths(path::Path) = subpaths(convert(OpenPolygonalPath, path))

subtract_paths(path1::Path, path2::Path) =
  subtract_paths(convert(ClosedPolygonalPath, path1), convert(ClosedPolygonalPath, path2))
subtract_paths(path1::ClosedPolygonalPath, path2::ClosedPolygonalPath) =
  closed_polygonal_path(
    subtract_polygon_vertices(
      path_vertices(path1),
      path_vertices(path2)))

## Smoothnesss
# A smooth curve means that the curve is differentiable up to an intended order
# In practice, it means that it does not have "corners"
is_smooth_path(path::Path) = false
is_smooth_path(pts::Locs) = false
is_smooth_path(path::Union{ArcPath,CircularPath,SplinePath}) = true

##Profiles are just predefined paths used for sections, usually centered at the origin
export rectangular_profile,
       circular_profile,
       top_aligned_rectangular_profile,
       bottom_aligned_rectangular_profile,
       i_profile, plus_profile

rectangular_profile(Width::Real=0.1, Height::Real=0.1; width::Real=Width, height::Real=Height) =
  centered_rectangular_path(u0(), width, height)

circular_profile(Radius::Real=0.1; radius::Real=Radius) =
  circular_path(u0(), radius)

top_aligned_rectangular_profile(Width::Real=0.1, Height::Real=0.1; width::Real=Width, height::Real=Height) =
  rectangular_path(xy(-width/2,-height), width, height)

bottom_aligned_rectangular_profile(Width::Real=0.1, Height::Real=0.1; width::Real=Width, height::Real=Height) =
  rectangular_path(xy(-width/2, 0), width, height)

i_profile(Width::Real=0.1, Height::Real=Width; width::Real=Width, height::Real=Height, web_thickness::Real=width/5, flange_thickness::Real=width/10) =
    let w2 = width/2,
        h2 = height/2,
        wt2 = web_thickness/2,
        ft = flange_thickness
      mirrored_on_x(
        mirrored_on_y(
          open_polygonal_path([
            x(wt2), xy(wt2, h2-ft), xy(w2, h2-ft), xy(w2, h2), y(h2)])))
    end

plus_profile(Width::Real=0.1, Height::Real=Width; width::Real=Width, height::Real=Height, thickness::Real=width/5) =
    let w2 = width/2,
        h2 = height/2,
        t2 = thickness/2
      mirrored_on_x(
        mirrored_on_y(
          open_polygonal_path([
            x(w2), xy(w2, t2), xy(t2, t2), xy(t2, h2), y(h2)])))
    end

export path_vertices_on
path_vertices_on(path::Path, p) =
  path_vertices_on(path_vertices(path), p)
path_vertices_on(vs::Locs, p) =
  [add_xyz(p, raw_point(v)...) for v in vs]

## Utility operations
Base.reverse(path::CircularPath) =
  circular_path(loc_from_o_vz(path.center, vz(-1)), path.radius)
Base.reverse(path::RectangularPath) =
  rectangular_path(loc_from_o_vz(path.corner, vz(-1)), path.dy, path.dx)
Base.reverse(path::OpenPolygonalPath) =
  open_polygonal_path(reverse(path_vertices(path)))
Base.reverse(path::ClosedPolygonalPath) =
  closed_polygonal_path(reverse(path_vertices(path)))
Base.reverse(path::Region) =
  region(reverse(outer_path(path)), reverse.(inner_paths(path))...)

mirrored_path(path::Path, p::Loc, v::Vec) =
  error("Mush be finished")

mirrored_on_x(path::OpenPolygonalPath) =
  join_paths(path, open_polygonal_path(reverse(map(p->xyz(p.x, -p.y, -p.z, p.cs), path_vertices(path)))))
mirrored_on_y(path::OpenPolygonalPath) =
  join_paths(path, open_polygonal_path(reverse(map(p->xyz(-p.x, p.y, -p.z, p.cs), path_vertices(path)))))
mirrored_on_z(path::OpenPolygonalPath) =
  join_paths(path, open_polygonal_path(reverse(map(p->xyz(-p.x, -p.y, p.z, p.cs), path_vertices(path)))))

export mesh
struct Mesh <: Path
  vertices::Locs
  faces::Vector{Vector{Int}}
end
mesh(vertices::Locs=[u0(), ux(), uy()], faces::Vector{Vector{Int}}=[[0,1,2]]) =
  Mesh(vertices, faces)

union(m::Mesh, ms...) =
  length(ms) == 0 ?
    m :
    let n = length(m.vertices)
      union(
        mesh(vcat(m.vertices, ms[1].vertices),
            vcat(m.faces, [face.+n for face in ms[1].faces])),
        ms[2:end]...)
    end

# The inverse of location_at_length
export length_at_location
length_at_location(path::Path, p::Loc, t0=0, t1=path_length(path)) =
  t1 - t0 < epsilon() ?
    (t0+t1)/2 :
    let ts = division(t0, t1, 10),
        pts = map(t->location_at_length(path, t), ts),
        idx = argmin(map(p1 -> distance(p, p1), pts))
      length_at_location(path, p, ts[max(1, idx-1)], ts[min(length(ts), idx+1)])
    end

export Grid, grid, grid_interpolator, location_at

grid_interpolator(ptss) =
  let (nu, nv) = size(ptss).-1,
      us = division(0, 1, nu),
      vs = division(0, 1, nv),
      rawptss = map(p->in_world(p).raw, ptss),
      order(n) = min(n, 5),
      ou = order(nu),
      ov = order(nv),
      xinterpolator = Spline2D(us, vs, map(p->p[1], rawptss), kx=ou, ky=ov),
      yinterpolator = Spline2D(us, vs, map(p->p[2], rawptss), kx=ou, ky=ov),
      zinterpolator = Spline2D(us, vs, map(p->p[3], rawptss), kx=ou, ky=ov)
    (xinterpolator, yinterpolator, zinterpolator)
  end
smooth_grid(ptss, level_u=2, level_v=2) =
  let (nu, nv) = size(ptss)
    if level_u > 1 && level_v > 1
	    let interpolator = grid_interpolator(ptss)
        [location_at(interpolator, u, v)
	       for u in division(0, 1, 4*nu-7), v in division(0, 1, 4*nv-7)] # 2->1, 3->5, 4->9, 5->13
      end
    elseif level_u > 1
      let interpolators = map(col->curve_interpolator(collect(col), false), eachcol(ptss))
        [location_at(interpolator, u) for u in division(0, 1, 4*nu-7), interpolator in interpolators]
      end
    elseif level_v > 1
      let interpolators = map(row->curve_interpolator(collect(row), false), eachrow(ptss))
        [location_at(interpolator, v) for interpolator in interpolators, v in division(0, 1, 4*nv-7)]
      end
    else
      ptss
    end
  end

struct Grid
  vertices::Array{Loc,2}
  interpolator::LazyParameter{Any}
end
grid(ptss) = Grid(ptss, LazyParameter(Any, ()->grid_interpolator(ptss)))

location_at(interpolator::Tuple{Spline2D,Spline2D,Spline2D}, u, v) =
  let o = xyz(evaluate.(interpolator, u, v)..., world_cs),
      vx = unitized(vxyz(evaluate_derivative.(interpolator, u, v, 1, 0)..., world_cs)),
      vy = unitized(vxyz(evaluate_derivative.(interpolator, u, v, 0, 1)..., world_cs))
    loc_from_o_vx_vy(o, vx, vy)
  end
location_at(grid::Grid, u, v) =
  location_at(grid.interpolator(), u, v)

# We need to support convertion to OBJ format
export write_obj
write_obj(io::IO, ptss::AbstractMatrix{<:Loc}, closed_u, closed_v, smooth_u, smooth_v, interpolator) =
  let (nu, nv) = size(ptss)
    if smooth_u || smooth_v
      ptss = [location_at(interpolator, u, v)
              for u in division(0, 1, nu*4), v in division(0, 1, nv*4)]
      (nu, nv) = size(ptss)
    end
    vertex(p) =
      println(io,"v  ",p.x," ",p.y," ",p.z)
    normal(n) =
      println(io,"vn ",n.x," ",n.y," ",n.z)
    face(p0, p1, p2, p3) = # OBJ format uses 1-based indexes
      let p0 = p0+1, p1=p1+1, p2=p2+1, p3=p3+1
        smooth_u || smooth_v ?
          println(io,"f ",p0,"//",p0," ",p1,"//",p1," ",p2,"//",p2," ",p3,"//",p3) :
          println(io,"f ",p0," ",p1," ",p2," ",p3)
      end
    quad_strip_faces(s, n) =
      foreach(face, zip(s:s+n-1, s+1:s+n, s+n+1:s+2*n-1, s+n:s+2*n-2))
    quad_strip_closed_faces(s, n) =
      begin
        quad_strip_faces(s, n)
        face(s+n-1, s, s+n, s+2*n-1)
      end
    # Vertices
    for j in 1:nv
      for i in 1:nu
        vertex(in_world(ptss[i,j]))
      end
    end
    # Normals
    if smooth_u || smooth_v
      for j in 1:nv
        for i in 1:nu
          normal(in_world(uz(ptss[i,j].cs)))
        end
      end
    end
    # Faces
    if closed_u
      for i in 0:nv-2
        quad_strip_closed_faces(i*nu, nu)
      end
      if closed_v
        foreach((p, q)->face(p, p+1, q+1, q), zip((nv-1)*nu:nv*nu-2, 0:nu-2))
        face(nv*nu-1, (nv-1)*nu, 0, nu-1)
      end
    else
      for i in 0:nv-2
        quad_strip_faces(i*nu, nu)
        if closed_v
          foreach((p, q)->face(p, p+1, q+1, q), zip((nv-1)*nu:nv*nu-1, 0:nu-1))
        end
      end
    end
  end


text:
# nested archimedean copulas

# Algorithms from:
# M. Hofert, `Efficiently sampling nested Archimedean copulas` Computational Statistics and Data Analysis 55 (2011) 57–70
# M. Hofert, 'Sampling  Archimedean copulas', Computational Statistics & Data Analysis, Volume 52, 2008
# McNeil, A.J., 2008. 'Sampling nested Archimedean copulas'. Journal of Statistical Computation and Simulation 78, 567–581.

#Basically we use Alg. 5 of McNeil, A.J., 2008. 'Sampling nested Archimedean copulas'.

"""
    Nested_Clayton_cop

Fields:
- children::Vector{Clayton_cop}  vector of children copulas
- m::Int ≧ 0 - number of additional marginals modeled by the parent copula only
- θ::Real - parameter of parent copula, domain θ > 0.

Nested Clayton copula: C_θ(C_ϕ₁(u₁₁, ..., u₁,ₙ₁), ..., C_ϕₖ(uₖ₁, ..., uₖ,ₙₖ), u₁ , ... uₘ).
If m > 0, the last m variables will be modeled by the parent copula only.

Constructor

    Nested_Clayton_cop(children::Vector{Clayton_cop}, m::Int, θ::Real)

Let ϕ be the vector of parameter of children copula, sufficient nesting condition requires
θ <= minimum(ϕ)

Constructor

    Nested_Clayton_cop(children::Vector{Clayton_cop}, m::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest
julia> a = Clayton_cop(2, 2.)
Clayton_cop(2, 2.0)

julia> Nested_Clayton_cop([a], 2, 0.5)
Nested_Clayton_cop(Clayton_cop[Clayton_cop(2, 2.0)], 2, 0.5)

julia> Nested_Clayton_cop([a, a], 2, 0.5)
Nested_Clayton_cop(Clayton_cop[Clayton_cop(2, 2.0), Clayton_cop(2, 2.0)], 2, 0.5)

```
"""
struct Nested_Clayton_cop{T}
  children::Vector{Clayton_cop{T}}
  m::Int
  θ::T
  function(::Type{Nested_Clayton_cop})(children::Vector{Clayton_cop{T}}, m::Int, θ::T) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      testθ(θ, "clayton")
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      maximum(ϕ) < θ+2*θ^2+750*θ^5 || @warn("θ << ϕ, marginals may not be uniform")
      new{T}(children, m, θ)
  end
  function(::Type{Nested_Clayton_cop})(children::Vector{Clayton_cop{T}}, m::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      θ = getθ4arch(ρ, "clayton", cor)
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      maximum(ϕ) < θ+2*θ^2+750*θ^5 || @warn("θ << ϕ, marginals may not be uniform")
      new{T}(children, m, θ)
  end
end

"""
    Nested_AMH_cop

Nested Ali-Mikhail-Haq copula, fields:
- children::Vector{AMH _cop}  vector of children copulas
- m::Int ≧ 0 - number of additional marginals modeled by the parent copula only
- θ::Real - parameter of parent copula, domain θ ∈ (0,1).

Nested Ali-Mikhail-Haq copula: C _θ(C _ϕ₁(u₁₁, ..., u₁,ₙ₁), ..., C _ϕₖ(uₖ₁, ..., uₖ,ₙₖ), u₁ , ... uₘ).
If m > 0, the last m variables will be modeled by the parent copula only.

Constructor

    Nested_AMH_cop(children::Vector{AMH_cop}, m::Int, θ::Real)

Let ϕ be the vector of parameter of children copula, sufficient nesting condition requires
θ <= minimum(ϕ)

Constructor

    Nested_AMH_cop(children::Vector{AMH_cop}, m::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest

julia> a = AMH_cop(2, .2)
AMH_cop(2, 0.2)

julia> Nested_AMH_cop([a, a], 2, 0.1)
Nested_AMH_cop(AMH_cop[AMH_cop(2, 0.2), AMH_cop(2, 0.2)], 2, 0.1)

```
"""
struct Nested_AMH_cop{T}
  children::Vector{AMH_cop{T}}
  m::Int
  θ::T
  function(::Type{Nested_AMH_cop})(children::Vector{AMH_cop{T}}, m::Int, θ::T) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      testθ(θ, "amh")
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, m, θ)
  end
  function(::Type{Nested_AMH_cop})(children::Vector{AMH_cop{T}}, m::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      θ = getθ4arch(ρ, "amh", cor)
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, m, θ)
  end
end

"""
    Nested_Frank_cop

Fields:
- children::Vector{Frank_cop}  vector of children copulas
- m::Int ≧ 0 - number of additional marginals modeled by the parent copula only
- θ::Real - parameter of parent copula, domain θ ∈ (0,∞).

Nested Frank copula: C _θ(C _ϕ₁(u₁₁, ..., u₁,ₙ₁), ..., C _ϕₖ(uₖ₁, ..., uₖ,ₙₖ), u₁ , ... uₘ).
If m > 0, the last m variables will be modeled by the parent copula only.

Constructor

    Nested_Frank_cop(children::Vector{Frank_cop}, m::Int, θ::Real)

Let ϕ be the vector of parameter of children copula, sufficient nesting condition requires
θ <= minimum(ϕ)

Constructor

    Nested_Frank_cop(children::Vector{Frank_ cop}, m::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctests

julia> a = Frank_cop(2, 5.)
Frank_cop(2, 5.0)

julia> Nested_Frank_cop([a, a], 2, 0.1)
Nested_Frank_cop(Frank_cop[Frank_cop(2, 5.0), Frank_cop(2, 5.0)], 2, 0.1)
```
"""
struct Nested_Frank_cop{T}
  children::Vector{Frank_cop{T}}
  m::Int
  θ::T
  function(::Type{Nested_Frank_cop})(children::Vector{Frank_cop{T}}, m::Int, θ::T) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      testθ(θ, "frank")
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, m, θ)
  end
  function(::Type{Nested_Frank_cop})(children::Vector{Frank_cop{T}}, m::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      θ = getθ4arch(ρ, "frank", cor)
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, m, θ)
  end
end

"""
    Nested_Gumbel_cop

Fields:
- children::Vector{Gumbel_cop}  vector of children copulas
- m::Int ≧ 0 - number of additional marginals modeled by the parent copula only
- θ::Real - parameter of parent copula, domain θ ∈ [1,∞).

Nested Gumbel copula: C _θ(C _ϕ₁(u₁₁, ..., u₁,ₙ₁), ..., C _ϕₖ(uₖ₁, ..., uₖ,ₙₖ), u₁ , ... uₘ).
If m > 0, the last m variables will be modeled by the parent copula only.

Constructor

    Nested_Gumbel_cop(children::Vector{Gumbel_cop}, m::Int, θ::Real)

Let ϕ be the vector of parameter of children copula, sufficient nesting condition requires
θ <= minimum(ϕ)

Constructor

    Nested_Gumbel_cop(children::Vector{Gumbel_cop}, m::Int, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest

julia> a = Gumbel_cop(2, 5.)
Gumbel_cop(2, 5.0)

julia> Nested_Gumbel_cop([a, a], 2, 2.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 5.0)], 2, 2.1)
```
"""
struct Nested_Gumbel_cop{T}
  children::Vector{Gumbel_cop{T}}
  m::Int
  θ::T
  function(::Type{Nested_Gumbel_cop})(children::Vector{Gumbel_cop{T}}, m::Int, θ::T) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      testθ(θ, "gumbel")
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, m, θ)
  end
  function(::Type{Nested_Gumbel_cop})(children::Vector{Gumbel_cop{T}}, m::Int, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      m >= 0 || throw(DomainError("not supported for m  < 0 "))
      θ = getθ4arch(ρ, "gumbel", cor)
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, m, θ)
  end
end

"""
    simulate_copula(t::Int, copula::Nested_Clayton_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations of data generated using Nested Clayton copula

```jldoctest

julia> Random.seed!(43);

julia> c1 = Clayton_cop(2, 2.)
Clayton_cop(2, 2.0)

julia> c2 = Clayton_cop(2, 3.)
Clayton_cop(2, 3.0)

julia> cp = Nested_Clayton_cop([c1, c2], 1, 1.1)
Nested_Clayton_cop(Clayton_cop[Clayton_cop(2, 2.0), Clayton_cop(2, 3.0)], 1, 1.1)

julia> simulate_copula(4, cp)
4×5 Array{Float64,2}:
0.514118  0.84089   0.870106  0.906233  0.739349
0.588245  0.85816   0.935308  0.944444  0.709009
0.59625   0.665947  0.483649  0.603074  0.153501
0.200051  0.304099  0.242572  0.177836  0.0851603
```
"""
function simulate_copula(t::Int, copula::Nested_Clayton_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
     n = [ch.n for ch in copula.children]
     n2 = sum(n)+copula.m
     U = zeros(T, t, n2)
     simulate_copula!(U, copula; rng = rng)
     return U
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Nested_Clayton_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the nested Clayton copula - Nested_Clayton_cop
N.o. marginals is size(U,2), these must be euqal to n.o. marginals of the copula

```jldoctest

julia> c1 = Clayton_cop(2, 2.)
Clayton_cop(2, 2.0)

julia> c2 = Clayton_cop(2, 3.)
Clayton_cop(2, 3.0)

julia> cp = Nested_Clayton_cop([c1, c2], 1, 1.1)
Nested_Clayton_cop(Clayton_cop[Clayton_cop(2, 2.0), Clayton_cop(2, 3.0)], 1, 1.1)

julia> U = zeros(5,5)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(U, cp)

julia> U
5×5 Array{Float64,2}:
 0.514118  0.84089    0.870106   0.906233  0.739349
 0.588245  0.85816    0.935308   0.944444  0.709009
 0.59625   0.665947   0.483649   0.603074  0.153501
 0.200051  0.304099   0.242572   0.177836  0.0851603
 0.120914  0.0683055  0.0586907  0.126257  0.519241
```
"""
function simulate_copula!(U::Matrix{T}, copula::Nested_Clayton_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
      m = copula.m
      θ = copula.θ
      children = copula.children
      ϕ = [ch.θ for ch in children]
      n = [ch.n for ch in children]
      n1 = vcat([collect(1:n[1])], [collect(cumsum(n)[i]+1:cumsum(n)[i+1]) for i in 1:length(n)-1])
      n2 = sum(n)+m
      size(U, 2) == n2 || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
      for j in 1:size(U,1)
         rand_vec = rand(rng, T, n2+1)
         U[j,:] = nested_clayton_gen(n1, ϕ, θ, rand_vec; rng=rng)
     end
end

"""
     simulate_copula(t::Int, copula::Nested_AMH_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations of data generated using Nested AMH copula

```jldoctest

julia> c1 = AMH_cop(2, .7)
AMH_cop(2, 0.7)

julia> c2 = AMH_cop(2, .8)
AMH_cop(2, 0.8)

julia> cp = Nested_AMH_cop([c1, c2], 1, 0.2)
Nested_AMH_cop(AMH_cop[AMH_cop(2, 0.7), AMH_cop(2, 0.8)], 1, 0.2)

julia> Random.seed!(43);

julia> simulate_copula(4, cp)
4×5 Array{Float64,2}:
 0.557393  0.902767  0.909853  0.938522  0.586068
 0.184204  0.866664  0.699134  0.226744  0.102932
 0.268634  0.383355  0.179023  0.533749  0.995958
 0.578143  0.840169  0.743728  0.963226  0.576695
```
"""
function simulate_copula(t::Int, copula::Nested_AMH_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
     n = [ch.n for ch in copula.children]
     n2 = sum(n)+copula.m
     U = zeros(T, t, n2)
     simulate_copula!(U, copula; rng = rng)
     return U
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Nested_AMH_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the nested AMH copula - Nested_AMH_cop
N.o. marginals is size(U,2), these must be euqal to n.o. marginals of the copula

```jldoctest

julia> c1 = AMH_cop(2, .7)
AMH_cop(2, 0.7)

julia> c2 = AMH_cop(2, .8)
AMH_cop(2, 0.8)

julia> cp = Nested_AMH_cop([c1, c2], 1, 0.2)
Nested_AMH_cop(AMH_cop[AMH_cop(2, 0.7), AMH_cop(2, 0.8)], 1, 0.2)

julia> Random.seed!(43);

julia> U = zeros(4,5)
4×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> simulate_copula!(U, cp)

julia> U
4×5 Array{Float64,2}:
 0.557393  0.902767  0.909853  0.938522  0.586068
 0.184204  0.866664  0.699134  0.226744  0.102932
 0.268634  0.383355  0.179023  0.533749  0.995958
 0.578143  0.840169  0.743728  0.963226  0.576695
```
"""
function simulate_copula!(U::Matrix{T}, copula::Nested_AMH_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    m = copula.m
    θ = copula.θ
    children = copula.children
    ϕ = [ch.θ for ch in children]
    n = [ch.n for ch in children]
    n1 = vcat([collect(1:n[1])], [collect(cumsum(n)[i]+1:cumsum(n)[i+1]) for i in 1:length(n)-1])
    n2 = sum(n)+m
    size(U, 2) == n2 || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
    for j in 1:size(U,1)
       rand_vec = rand(rng, T, n2+1)
       U[j,:] = nested_amh_gen(n1, ϕ, θ, rand_vec; rng=rng)
   end
end

"""
    simulate_copula(t::Int, copula::Nested_Frank_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations of data generated using Nested Frank copula
```jldoctest

julia> c1 = Frank_cop(2, 4.)
Frank_cop(2, 4.0)

julia> c2 = Frank_cop(2, 5.)
Frank_cop(2, 5.0)

julia> c = Nested_Frank_cop([c1, c2],1, 2.0)
Nested_Frank_cop(Frank_cop[Frank_cop(2, 4.0), Frank_cop(2, 5.0)], 1, 2.0)

julia> Random.seed!(43);

julia> simulate_copula(1, c)
1×5 Array{Float64,2}:
 0.642765  0.901183  0.969422  0.9792  0.74155
```
"""
function simulate_copula(t::Int, copula::Nested_Frank_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
     n = [ch.n for ch in copula.children]
     n2 = sum(n)+copula.m
     U = zeros(T, t, n2)
     simulate_copula!(U, copula; rng = rng)
     return U
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Nested_Frank_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the nested Frank copula a - Nested_Frank_cop
N.o. marginals is size(U,2), these must be euqal to n.o. marginals of the copula

```jldoctest

julia> c1 = Frank_cop(2, 4.)
Frank_cop(2, 4.0)

julia> c2 = Frank_cop(2, 5.)
Frank_cop(2, 5.0)

julia> c = Nested_Frank_cop([c1, c2],1, 2.0)
Nested_Frank_cop(Frank_cop[Frank_cop(2, 4.0), Frank_cop(2, 5.0)], 1, 2.0)

julia> U = zeros(1,5)
1×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(U, c)

julia> U
1×5 Array{Float64,2}:
 0.642765  0.901183  0.969422  0.9792  0.74155

```
"""
function simulate_copula!(U::Matrix{T}, copula::Nested_Frank_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    m = copula.m
    θ = copula.θ
    children = copula.children
    ϕ = [ch.θ for ch in children]
    n = [ch.n for ch in children]
    n2 = sum(n)+m
    size(U, 2) == n2 || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))

    ws = [logseriescdf(1-exp(theta)) for theta in ϕ]
    n1 = vcat([collect(1:n[1])], [collect(cumsum(n)[i]+1:cumsum(n)[i+1]) for i in 1:length(n)-1])
    w = logseriescdf(1-exp(-θ))
    for j in 1:size(U,1)
       rand_vec = rand(rng, T, n2+1)
       U[j,:] = nested_frank_gen(n1, ϕ, θ, rand_vec, w, ws; rng=rng)
   end
end

"""
    simulate_copula(t::Int, copula::Nested_Gumbel_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations of data generated using  Nested Gumbel copula

```jldoctest

julia> Random.seed!(43);

julia> c1 = Gumbel_cop(2, 2.)
Gumbel_cop(2, 2.0)

julia> c2 = Gumbel_cop(2, 3.)
Gumbel_cop(2, 3.0)

julia> cp = Nested_Gumbel_cop([c1, c2], 1, 1.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 2.0), Gumbel_cop(2, 3.0)], 1, 1.1)

julia> simulate_copula(4, cp)
4×5 Array{Float64,2}:
 0.387085   0.693399   0.94718   0.953776  0.583379
 0.0646972  0.0865914  0.990691  0.991127  0.718803
 0.966896   0.709233   0.788019  0.855622  0.755476
 0.272487   0.106996   0.756052  0.834068  0.661432
```
"""
function simulate_copula(t::Int, copula::Nested_Gumbel_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
     n = [ch.n for ch in copula.children]
     n2 = sum(n)+copula.m
     U = zeros(T, t, n2)
     simulate_copula!(U, copula; rng = rng)
     return U
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Nested_Gumbel_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the nested Gumbel copula - Nested_Gumbel_cop
N.o. marginals is size(U,2), these must be euqal to n.o. marginals of the copula

```jldoctest
julia> c1 = Gumbel_cop(2, 2.)
Gumbel_cop(2, 2.0)

julia> c2 = Gumbel_cop(2, 3.)
Gumbel_cop(2, 3.0)

julia> cp = Nested_Gumbel_cop([c1, c2], 1, 1.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 2.0), Gumbel_cop(2, 3.0)], 1, 1.1)

julia> u = zeros(4,5)
4×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(u, cp)

julia> u
4×5 Array{Float64,2}:
 0.387085   0.693399   0.94718   0.953776  0.583379
 0.0646972  0.0865914  0.990691  0.991127  0.718803
 0.966896   0.709233   0.788019  0.855622  0.755476
 0.272487   0.106996   0.756052  0.834068  0.661432
```
"""
function simulate_copula!(U::Matrix{T}, copula::Nested_Gumbel_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    m = copula.m
    θ = copula.θ
    children = copula.children
    ϕ = [ch.θ for ch in children]
    n = [ch.n for ch in children]
    n1 = vcat([collect(1:n[1])], [collect(cumsum(n)[i]+1:cumsum(n)[i+1]) for i in 1:length(n)-1])
    n2 = sum(n)+m
    size(U, 2) == n2 || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))

    for j in 1:size(U,1)
       rand_vec = rand(rng, T, n2)
       U[j,:] = nested_gumbel_gen(n1, ϕ, θ, rand_vec; rng=rng)
   end
end

"""
    Double_Nested_Gumbel_cop

Fields:
- children::Vector{Nested _Gumbel _cop}  vector of children copulas
- θ::Real - parameter of parent copula, domain θ ∈ [1,∞).

Constructor

    Double_Nested_Gumbel _cop(children::Vector{Nested_Gumbel_cop}, θ::Real)
requires sufficient nesting condition for θ and child copulas.

Constructor

    Doulbe_Nested_Gumbel_cop(children::Vector{Nested_Gumbel_cop}, θ::Real, cor::Type{<:CorrelationType})

For computing copula parameter from expected correlation use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlation in the place of θ  in the constructor.
The copula parameter will be computed then. The correlation must be greater than zero.

```jldoctest

julia> a = Gumbel_cop(2, 5.)
Gumbel_cop(2, 5.0)

julia> b = Gumbel_cop(2, 6.)
Gumbel_cop(2, 6.0)

julia> c = Gumbel_cop(2, 5.5)
Gumbel_cop(2, 5.5)

julia> p1 = Nested_Gumbel_cop([a,b], 1, 2.)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0)

julia> p2 = Nested_Gumbel_cop([c], 2, 2.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)

julia> Double_Nested_Gumbel_cop([p1, p2], 1.5)
Double_Nested_Gumbel_cop(Nested_Gumbel_cop[Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0), Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)], 1.5)
```
"""
struct Double_Nested_Gumbel_cop{T}
  children::Vector{Nested_Gumbel_cop{T}}
  θ::T
  function(::Type{Double_Nested_Gumbel_cop})(children::Vector{Nested_Gumbel_cop{T}}, θ::T) where T <: Real
      testθ(θ, "gumbel")
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, θ)
  end
  function(::Type{Double_Nested_Gumbel_cop})(children::Vector{Nested_Gumbel_cop{T}}, ρ::T, cor::Type{<:CorrelationType}) where T <: Real
      θ = getθ4arch(ρ, "gumbel", cor)
      ϕ = [ch.θ for ch in children]
      θ <= minimum(ϕ) || throw(DomainError("violated sufficient nesting condition"))
      new{T}(children, θ)
  end
end

"""
    simulate_copula(t::Int, copula::Double_Nested_Gumbel_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Simulate t realization of the Double Nested Gumbel copula i.e.

    Double_Nested_Gumbel_cop(vec_of_children, θ)

```jldoctest
julia> a = Gumbel_cop(2, 5.)
Gumbel_cop(2, 5.0)

julia> b = Gumbel_cop(2, 6.)
Gumbel_cop(2, 6.0)

julia> c = Gumbel_cop(2, 5.5)
Gumbel_cop(2, 5.5)

julia> p1 = Nested_Gumbel_cop([a,b], 1, 2.)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0)

julia> p2 = Nested_Gumbel_cop([c], 2, 2.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)

julia> copula = Double_Nested_Gumbel_cop([p1, p2], 1.5)
Double_Nested_Gumbel_cop(Nested_Gumbel_cop[Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0), Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)], 1.5)

julia> Random.seed!(43);

julia> simulate_copula(5, copula)
5×9 Array{Float64,2}:
 0.598555   0.671584   0.8403     0.846844  0.634609  0.686927  0.693906  0.651968    0.670812
 0.0518892  0.191236   0.0803859  0.104325  0.410727  0.529354  0.557387  0.370518    0.592302
 0.367914   0.276196   0.382616   0.470171  0.264135  0.144503  0.13097   0.00687015  0.01417
 0.632727   0.596879   0.244176   0.338809  0.58771   0.147539  0.219099  0.287937    0.0569943
 0.310365   0.0483216  0.119312   0.107155  0.336619  0.279602  0.262756  0.438432    0.403061
```
"""
function simulate_copula(t::Int, copula::Double_Nested_Gumbel_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    v = copula.children
    ns = [[ch.n for ch in vs.children] for vs in v]
    dims = sum([sum(ns[i])+v[i].m for i in 1:length(v)])
     U = zeros(T, t, dims)
     simulate_copula!(U, copula; rng = rng)
     return U
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Double_Nested_Gumbel_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the double nested Gumbel copula - Double_Nested_Gumbel_cop
N.o. marginals is size(U,2), these must be euqal to n.o. marginals of the copula

```jldoctest
julia> a = Gumbel_cop(2, 5.)
Gumbel_cop(2, 5.0)

julia> b = Gumbel_cop(2, 6.)
Gumbel_cop(2, 6.0)

julia> c = Gumbel_cop(2, 5.5)
Gumbel_cop(2, 5.5)

julia> p1 = Nested_Gumbel_cop([a,b], 1, 2.)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0)

julia> p2 = Nested_Gumbel_cop([c], 2, 2.1)
Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)

julia> copula = Double_Nested_Gumbel_cop([p1, p2], 1.5)
Double_Nested_Gumbel_cop(Nested_Gumbel_cop[Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.0), Gumbel_cop(2, 6.0)], 1, 2.0), Nested_Gumbel_cop(Gumbel_cop[Gumbel_cop(2, 5.5)], 2, 2.1)], 1.5)

julia> u = zeros(3,9)
3×9 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(u, copula)

julia> u
3×9 Array{Float64,2}:
 0.598555   0.671584  0.8403     0.846844  0.634609  0.686927  0.693906  0.651968    0.670812
 0.0518892  0.191236  0.0803859  0.104325  0.410727  0.529354  0.557387  0.370518    0.592302
 0.367914   0.276196  0.382616   0.470171  0.264135  0.144503  0.13097   0.00687015  0.01417
```
"""

function simulate_copula!(U::Matrix{T}, copula::Double_Nested_Gumbel_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    θ = copula.θ
    v = copula.children
    ns = [[ch.n for ch in vs.children] for vs in v]
    Ψs = [[ch.θ for ch in vs.children] for vs in v]
    dims = sum([sum(ns[i])+v[i].m for i in 1:length(v)])
    size(U, 2) == dims || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))

    for j in 1:size(U,1)
        X = T[]
        for k in 1:length(v)
            n = ns[k]
            n1 = vcat([collect(1:n[1])], [collect(cumsum(n)[i]+1:cumsum(n)[i+1]) for i in 1:length(n)-1])
            n2 = sum(n)+v[k].m
            rand_vec = rand(rng, T, n2)
            X = vcat(X, nested_gumbel_gen(n1, Ψs[k], v[k].θ./θ, rand_vec; rng = rng))
        end
        X = -log.(X)./levyel(θ, rand(rng), rand(rng))
        U[j,:] = exp.(-X.^(1/θ))
    end
end

"""
    Hierarchical_Gumbel_cop

Fields:
- n::Int - number of marginals
- θ::Vector{Real} - vector of parameters, must be decreasing  and θ[end] ≧ 1, for the
sufficient nesting condition to be fulfilled.

The hierarchically nested Gumbel copula C_θₙ₋₁(C_θₙ₋₂( ... C_θ₂(C_θ₁(u₁, u₂), u₃)...uₙ₋₁) uₙ)

Constructor

    Hierarchical_Gumbel_cop(θ::Vector{Real})

Constructor

    Hierarchical_Gumbel_cop(ρ::Vector{Real}, cor::Type{<:CorrelationType})

For computing copula parameters from expected correlations use empty type cor::Type{<:CorrelationType} where
SpearmanCorrelation <:CorrelationType and KendallCorrelation<:CorrelationType. If used cor put expected correlations in the place of θ  in the constructor.
The copula parameters will be computed then. The correlation must be greater than zero.

```jldoctest

julia> c = Hierarchical_Gumbel_cop([5., 4., 3.])
Hierarchical_Gumbel_cop(4, [5.0, 4.0, 3.0])

julia> c = Hierarchical_Gumbel_cop([0.95, 0.5, 0.05], KendallCorrelation)
Hierarchical_Gumbel_cop(4, [19.999999999999982, 2.0, 1.0526315789473684])
```
"""
struct Hierarchical_Gumbel_cop{T}
  n::Int
  θ::Vector{T}
  function(::Type{Hierarchical_Gumbel_cop})(θ::Vector{T}) where T <: Real
      testθ(θ[end], "gumbel")
      issorted(θ; rev=true) || throw(DomainError("violated sufficient nesting condition, parameters must be descending"))
      new{T}(length(θ)+1, θ)
  end
  function(::Type{Hierarchical_Gumbel_cop})(ρ::Vector{T}, cor::Type{<:CorrelationType}) where T <: Real
      θ = map(i -> getθ4arch(ρ[i], "gumbel", cor), 1:length(ρ))
      issorted(θ; rev=true) || throw(DomainError("violated sufficient nesting condition, parameters must be descending"))
      new{T}(length(θ)+1, θ)
  end
end

"""
    simulate_copula(t::Int, copula::Hierarchical_Gumbel_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Returns t realizations of multivariate data from hierarchically nested Gumbel copula, i.e.

    Hierarchical_Gumbel_cop(θ)

```jldoctest
julia> using Random

julia> Random.seed!(43);

julia> c = Hierarchical_Gumbel_cop([5., 4., 3.])
Hierarchical_Gumbel_cop(4, [5.0, 4.0, 3.0])

julia> simulate_copula(3, c)
3×4 Array{Float64,2}:
 0.100353  0.207903  0.0988337  0.0431565
 0.347417  0.217052  0.223734   0.042903
 0.73617   0.347349  0.168348   0.410963
```
"""
function simulate_copula(t::Int, copula::Hierarchical_Gumbel_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
    U = zeros(T, t, copula.n)
    simulate_copula!(U, copula; rng = rng)
    return U
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Hierarchical_Gumbel_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the hierachically nested Gumbel copula - Hierarchical_Gumbel_cop
N.o. marginals is size(U,2), these must be euqal to n.o. marginals of the copula i.e. copula.n

```jldoctest

julia> c = Hierarchical_Gumbel_cop([5., 4., 3.])
Hierarchical_Gumbel_cop(4, [5.0, 4.0, 3.0])

julia> Random.seed!(43);

julia> u = zeros(3,4)
3×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> simulate_copula!(u, c)

julia> u
3×4 Array{Float64,2}:
 0.100353  0.207903  0.0988337  0.0431565
 0.347417  0.217052  0.223734   0.042903
 0.73617   0.347349  0.168348   0.410963
```
"""
function simulate_copula!(U::Matrix{T}, copula::Hierarchical_Gumbel_cop{T}; rng::AbstractRNG = Random.GLOBAL_RNG) where T <: Real
  θ = copula.θ
  θ = vcat(θ, [1.])
  size(U, 2) == copula.n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))

  for j in 1:size(U,1)
      X = rand(rng, T)
      for i in 1:(copula.n-1)
          X = gumbel_step(vcat(X, rand(rng)), θ[i], θ[i+1]; rng = rng)
      end
     U[j,:] = X
    end
end

"""
    nested_gumbel_gen(n::Vector{Vector{Int}}, ϕ::Vector{Real},
                         θ::Real, rand_vec::Vector{Real}; rng::AbstractRNG)

Convert a vector of random independnet elements to such sampled from the
Nested Gumbel copula
"""
function nested_gumbel_gen(n::Vector{Vector{Int}}, ϕ::Vector{T},
                         θ::T, rand_vec::Vector{T}; rng::AbstractRNG) where T <: Real
    V0 = levyel(θ, rand(rng, T), rand(rng, T))
    u = copy(rand_vec)
    for i in 1:length(n)
      u[n[i]] = gumbel_step(rand_vec[n[i]], ϕ[i], θ; rng = rng)
    end
    u = -log.(u)./V0
    return exp.(-u.^(1/θ))
end

"""
    nested_amh_gen(n::Vector{Vector{Int}}, ϕ::Vector{Real},
                         θ::Real, rand_vec::Vector{Real}; rng::AbstractRNG)
"""
function nested_amh_gen(n::Vector{Vector{Int}}, ϕ::Vector{T},
                         θ::T, rand_vec::Vector{T}; rng::AbstractRNG) where T <: Real
    V0 = 1 .+ quantile.(Geometric(1-θ), rand_vec[end])
    u = copy(rand_vec[1:end-1])
    for i in 1:length(n)
      u[n[i]] = amh_step(rand_vec[n[i]], V0, ϕ[i], θ; rng = rng)
    end
    u = -log.(u)./V0
    return (1-θ) ./(exp.(u) .-θ)
end

"""
    nested_frank_gen(n::Vector{Vector{Int}}, ϕ::Vector{Real}, θ::Real, rand_vec::Vector{Real}, logseries::Vector{Real},
                         logseries_children::Vector{Vector{Real}};
                         rng::AbstractRNG)
"""
function nested_frank_gen(n::Vector{Vector{Int}}, ϕ::Vector{T},
                         θ::T, rand_vec::Vector{T}, logseries::Vector{T},
                         logseries_children::Vector{Vector{T}};
                         rng::AbstractRNG) where T <: Real
    V0 = findlast(logseries .< rand_vec[end])
    u = copy(rand_vec[1:end-1])
    for i in 1:length(n)
      u[n[i]] = frank_step(rand_vec[n[i]], V0, ϕ[i], θ, logseries_children[i]; rng = rng)
    end
    u = -log.(u)./V0
    return -log.(1 .+exp.(-u) .*(exp(-θ)-1)) ./θ
end

"""
    nested_clayton_gen(n::Vector{Vector{Int}}, ϕ::Vector{Real}, θ::Real, rand_vec::Vector{Real}; rng::AbstractRNG = Random.GLOBAL_RNG)
"""
function nested_clayton_gen(n::Vector{Vector{Int}}, ϕ::Vector{T},
                         θ::T, rand_vec::Vector{T}; rng::AbstractRNG) where T <: Real
    V0 = gamma_inc_inv(1/θ, T(1.), T(1.)-rand_vec[end])
    #quantile.(Gamma(1/θ, 1), rand_vec[end])
    u = copy(rand_vec[1:end-1])
    for i in 1:length(n)
      u[n[i]] = clayton_step(rand_vec[n[i]], V0, ϕ[i], θ; rng = rng)
    end
    u = -log.(u)./V0
    return (1 .+ u).^(-1/θ)
end

"""
    gumbel_step(u::Vector{Real}, ϕ::Real, θ::Real; rng::AbstractRNG)
"""
function gumbel_step(u::Vector{T}, ϕ::T, θ::T; rng::AbstractRNG) where T <: Real
    u = -log.(u)./levyel(ϕ/θ, rand(rng, T), rand(rng, T))
    return exp.(-u.^(θ/ϕ))
end

"""
    clayton_step(u::Vector{Real}, V0::Real, ϕ::Real, θ::Real; rng::AbstractRNG)
"""
function clayton_step(u::Vector{T}, V0::T, ϕ::T, θ::T; rng::AbstractRNG) where T <: Real
    u = -log.(u)./tiltedlevygen(V0, ϕ/θ; rng = rng)
    return exp.(V0.-V0.*(1 .+u).^(θ/ϕ))
end

"""
    frank_step(u::Vector{Real}, V0::Int, ϕ::Real, θ::Real, logseries_child::Vector{Real}; rng::AbstractRNG)
"""
function frank_step(u::Vector{T}, V0::Int, ϕ::Real, θ::T, logseries_child::Vector{T}; rng::AbstractRNG) where T <: Real
    u = -log.(u)./nestedfrankgen(ϕ, θ, V0, logseries_child; rng = rng)
    X = (1 .-(1 .-exp.(-u)*(1-exp(-ϕ))).^(θ/ϕ))./(1-exp(-θ))
    return X.^V0
end
"""
    amh_step(u::Vector{Real}, V0::Real, ϕ::Real, θ::Real; rng::AbstractRNG)
"""
function amh_step(u::Vector{T}, V0::T, ϕ::T, θ::T; rng::AbstractRNG) where T <: Real
    # TODO this need to be changed for BigFloat
    w = quantile(NegativeBinomial(V0, (1-ϕ)/(1-θ)), rand(rng, T))
    u = -log.(u)./(V0 + w)
    X = ((exp.(u) .-ϕ) .*(1-θ) .+θ*(1-ϕ)) ./(1-ϕ)
    return X.^(-V0)
end

"""
nestedcopulag(copula::String, ns::Vector{Vector{Int}}, ϕ::Vector{Real}, θ::Real, r::Matrix{Real})

Given [0,1]ᵗˣˡ ∋ r, returns t realizations of l-1 variate data from nested archimedean copula

```jldoctest
julia> Random.seed!(43)

julia> nestedcopulag("clayton", [[1,2],[3,4]], [2., 3.], 1.1, [0.1 0.2 0.3 0.4 0.5; 0.2 0.3 0.4 0.5 0.6])
julia> nestedcopulag("clayton", [[1,2],[3,4]], [2., 3.], 1.1, [0.1 0.2 0.3 0.4 0.5; 0.2 0.3 0.4 0.5 0.6])
2×4 Array{Float64,2}:
 0.153282  0.182421  0.374228  0.407663
 0.69035   0.740927  0.254842  0.279192
```
"""
function nestedcopulag(copula::String, ns::Vector{Vector{Int}}, ϕ::Vector{T}, θ::T,
                                                        r::Matrix{T}; rng::AbstractRNG) where T <: Real
    t = size(r,1)
    n = size(r,2)-1
    u = zeros(T, t, n)
    if copula == "clayton"
        for j in 1:t
            u[j,:] = nested_clayton_gen(ns, ϕ, θ, r[j,:]; rng = rng)
        end
    elseif copula == "amh"
        for j in 1:t
            u[j,:] = nested_amh_gen(ns, ϕ, θ, r[j,:]; rng = rng)
        end
    elseif copula == "frank"
        ws = [logseriescdf(1-exp(theta)) for theta in ϕ]
        w = logseriescdf(1-exp(-θ))
        for j in 1:t
            u[j,:] = nested_frank_gen(ns, ϕ, θ, r[j,:], w, ws; rng = rng)
        end
    elseif copula == "gumbel"
        v = r[:,end]
        p = invperm(sortperm(v))
        V0 = [levyel(θ, rand(rng), rand(rng)) for i in 1:t]
        V0 = sort(V0)[p]
        for j in 1:t
            rand_vec = r[j,1:end-1]
            x = copy(rand_vec)
            for i in 1:length(ϕ)
              x[ns[i]] = gumbel_step(rand_vec[ns[i]], ϕ[i], θ; rng = rng)
            end
            x = -log.(x)./V0[j]
            u[j,:] = exp.(-x.^(1/θ))
        end
    end
    return u
end

ext: jl


text:
module ToggleableAsserts

export @toggled_assert, toggle

assert_toggle() = true

macro toggled_assert(cond, text=nothing)
    if text==nothing
        assert_stmt = esc(:(@assert $cond))
    else
        assert_stmt = esc(:(@assert $cond $text))
    end
    :(assert_toggle() ? $assert_stmt  : nothing)
end

const toggle_lock = ReentrantLock()

function toggle(enable::Bool)
    lock(toggle_lock) do
        @eval ToggleableAsserts assert_toggle() = $enable
        on_or_off = enable ? "on." : "off."
        @info "Toggleable asserts turned "*on_or_off
    end
end

end # module


text:
"""
    loadnmr(template::String)

Main function for loading nmrPipe or bruker format NMR data.

Returns an NMRData array, or throws an NMRToolsException is there is a problem.

# Examples

nmrPipe import:

```julia
loadnmr("exampledata/1D_1H/test.ft1");
loadnmr("exampledata/1D_19F/test.ft1");
loadnmr("exampledata/2D_HN/test.ft2");
loadnmr("exampledata/pseudo2D_XSTE/test.ft1");
loadnmr("exampledata/pseudo3D_HN_R2/ft/test%03d.ft2");
```

bruker import:

```julia
loadnmr("exampledata/1D_19F/pdata/1/1r");
```
"""
function loadnmr(filename::String, acqusfilename::Union{String,Nothing}=nothing)
    # load the data
    ispipe = occursin(r"[a-z0-9]+(\%03d)?\.ft[123]?$", filename)
    if ispipe
        # match files ending in .ft, .ft1, .ft2 or .ft3 => NMRPipe
        spectrum = loadnmrpipe(filename)
        metadata(spectrum)[:format] = :NMRPipe
    else
        # TODO bruker import
        throw(NMRToolsException("unknown file format for loadnmr\ntemplate = " * template))
    end

    metadata(spectrum)[:filename] = filename
    metadata(spectrum)[:NMRTools] = true

    # locate the acqus file
    if isnothing(acqusfilename)
        base = basename(filename)
        dir = dirname(filename)
        if ispipe
            # assume that acqus is in same directory as test.ft2 files, and one up from ft/test%03d.ft2 files
            if occursin("%", filename)
                dir = dirname(dir) # this moves up one directory
            end
            acqusfilename = joinpath(dir, "acqus")
        end
    end

    # parse the acqus file
    if ispath(acqusfilename)
        acqusmetadata = parseacqus(acqusfilename)
        metadata(spectrum)[:acqus] = acqusmetadata
        metadata(spectrum)[:ns] = acqusmetadata["NS"]
        metadata(spectrum)[:rg] = acqusmetadata["RG"]
        metadata(spectrum)[:pulseprogram] = acqusmetadata["PULPROG"]
        metadata(spectrum)[:acqusfilename] = acqusfilename
    else
        @warn "cannot locate acqus file for $(filename) - some metadata will be missing"
    end

    # load the title file
    titlefilename = joinpath(dirname(acqusfilename), "pdata", "1", "title")
    if ispath(titlefilename)
        title = read(titlefilename, String)
        metadata(spectrum)[:title] = strip(title)
        titleline1 = split(title, "\n")[1]
        label!(spectrum, titleline1)
    else
        @warn "cannot locate title file for $(filename) - some metadata will be missing"
    end

    # populate the spectrum noise level
    estimatenoise!(spectrum)

    return spectrum
end

expandpipetemplate(template, n) = replace(template, "%03d"=>lpad(n,3,"0"))

function loadnmrpipe(filename::String)
    if occursin("%", filename)
        # something like test%03d.ft2 - pseudo3D
        filename1 = expandpipetemplate(filename, 1)
    else
        filename1 = filename
    end
    # parse the header
    md, mdax = loadnmrpipeheader(filename1)

    ndim = md[:ndim]
    if ndim == 1
        loadnmrpipe1d(filename, md, mdax)
    elseif ndim == 2
        loadnmrpipe2d(filename, md, mdax)
    elseif ndim == 3
        loadnmrpipe3d(filename, md, mdax)
    else
        throw(NMRToolsException("can't load data $(filename), unsupported number of dimensions."))
    end
end

function loadnmrpipeheader(filename::String)
    # check file exists
    isfile(filename) || throw(NMRToolsException("cannot load $(filename), not a recognised file."))

    # preallocate header
    header = zeros(Float32, 512)

    # read the file
    open(filename) do f
        read!(f, header)
    end

    # parse the header, returning md and mdax
    parsenmrpipeheader(header)
end

"""
    parsenmrpipeheader(header)

Pass a 512 x 4 byte array containing an nmrPipe header file, and returns dictionaries of metadata.
The nmrPipe header format is defined in `fdatap.h`.

# Return values
- `md`: dictionary of spectrum metadata
- `mdax`: array of dictionaries containing axis metadata

# Examples
```julia
md, mdax = parsenmrpipeheader(header)
```
"""
function parsenmrpipeheader(header::Vector{Float32})
    # declare constants for clarity in accessing nmrPipe header

    # general parameters (independent of number of dimensions)
    genpars = Dict(
        :FDMAGIC => 1,
        :FDFLTFORMAT => 2,
        :FDFLTORDER => 3,
        :FDPIPEFLAG => (Int, 58),  #Dimension code of data stream
        :FDCUBEFLAG => (Int, 448),
        :FDSIZE => (Int, 100),
        :FDSPECNUM => (Int, 220),
        :FDQUADFLAG => (Int, 107),
        :FD2DPHASE => (Int, 257),
        :FDTRANSPOSED => (Int, 222),
        :FDDIMCOUNT => (Int, 10),
        :FDDIMORDER => (Int, 25:28),
        :FDFILECOUNT => (Int, 443),
        :FDSIZE => (Int, (100, 220, 16, 33)),
    )
    # dimension-specific parameters
    dimlabels = (17:18, 19:20, 21:22, 23:24) # UInt8
    dimpars = Dict(
        :FDAPOD => (Int, (96, 429, 51, 54, )),
        :FDSW => (101, 230, 12, 30, ),
        :FDOBS => (120, 219, 11, 29, ),
        :FDOBSMID => (379, 380, 381, 382, ),
        :FDORIG => (102, 250, 13, 31, ),
        :FDUNITS => (Int, (153, 235, 59, 60, )),
        :FDQUADFLAG => (Int, (57, 56, 52, 55, )),
        :FDFTFLAG => (Int, (221, 223, 14, 32, )),
        :FDCAR => (67, 68, 69, 70, ),
        :FDCENTER => (80, 81, 82, 83, ),
        :FDOFFPPM => (481, 482, 483, 484, ),
        :FDP0 => (110, 246, 61, 63, ),
        :FDP1 => (111, 247, 62, 64, ),
        :FDAPODCODE => (Int, (414, 415, 401, 406, )),
        :FDAPODQ1 => (416, 421, 402, 407, ),
        :FDAPODQ2 => (417, 422, 403, 408, ),
        :FDAPODQ3 => (418, 423, 404, 409, ),
        :FDLB => (112, 244, 373, 374, ),
        :FDGB => (375, 276, 377, 378, ),
        :FDGOFF => (383, 384, 385, 386, ),
        :FDC1 => (419, 424, 405, 410, ),
        :FDZF => (Int, (109, 438, 439, 440, )),
        :FDX1 => (Int, (258, 260, 262, 264, )),
        :FDXN => (Int, (259, 261, 263, 265, )),
        :FDFTSIZE => (Int, (97, 99, 201, 202, )),
        :FDTDSIZE => (Int, (387, 388, 389, 390, ))
    )

    # check correct endian-ness
    if header[genpars[:FDFLTORDER]] ≉ 2.345
        header = bswap.(header)
    end

    # populate spectrum metadata
    md = Dict{Symbol,Any}()
    for k in keys(genpars)
        x = genpars[k]
        # cast to new type if needed
        if x[1] isa Type
            if length(x[2]) == 1
                md[k] = x[1].(header[x[2]])
            else
                # array parameter
                md[k] = [x[1].(header[i]) for i in x[2]]
            end
        else
            md[k] = Float64(header[x])
        end
    end

    # add some metadata in nice format
    ndim = md[:FDDIMCOUNT]
    md[:ndim] = ndim

    # populate metadata for each dimension
    axesmd = []
    for i=1:ndim
        d = Dict{Symbol,Any}()
        for k in keys(dimpars)
            x = dimpars[k]
            # cast to new type if needed
            if x[1] isa Type
                d[k] = x[1].(header[x[2][i]])
            else
                d[k] = Float64(header[x[i]])
            end
        end
        # add label, removing null characters
        tmp = reinterpret(UInt8, header[dimlabels[i]])
        d[:label] = String(tmp[tmp .≠ 0])

        # add some data in a nice format
        if d[:FDFTFLAG] == 0
            d[:pseudodim] = true
            d[:npoints] = d[:FDTDSIZE]
            d[:val] = 1:d[:npoints] # coordinates for this dimension are just 1 to N
        else # frequency domain
            d[:pseudodim] = false
            d[:td] = d[:FDTDSIZE]
            d[:tdzf] = d[:FDFTSIZE]
            d[:bf] = d[:FDOBS]
            d[:swhz] = d[:FDSW]
            d[:swppm] = d[:swhz] / d[:bf]
            d[:offsetppm] = d[:FDCAR]
            d[:offsethz] = (d[:offsetppm]*1e-6 + 1) * d[:bf]
            d[:sf] = d[:offsethz]*1e-6 + d[:bf]
            if d[:FDX1] == 0 && d[:FDXN] == 0
                d[:region] = missing
                d[:npoints] = d[:tdzf]
            else
                d[:region] = d[:FDX1] : d[:FDXN]
                d[:npoints] = length(d[:region])
            end

            edge_frq = d[:FDORIG]
            # calculate chemical shift values
            cs_at_edge = edge_frq / d[:bf]
            cs_at_other_edge = cs_at_edge + d[:swppm]
            x = range(cs_at_other_edge, cs_at_edge, length=d[:npoints]+1);
            d[:val] = x[2:end];
            # # alternative calculation (without extraction) - this agrees OK
            # x = range(d[:offsetppm] + 0.5*d[:swppm], d[:offsetppm] - 0.5*d[:swppm], length=d[:npoints]+1)
            # d[:val2] = x[1:end-1]

            # create a representation of the window function
            # calculate acquisition time = td / sw
            taq = d[:td] / d[:swhz]
            #:FDAPODCODE => "Window function used (0=none, 1=SP, 2=EM, 3=GM, 4=TM, 5=ZE, 6=TRI)",
            w = d[:FDAPODCODE]
            q1 = d[:FDAPODQ1]
            q2 = d[:FDAPODQ2]
            q3 = d[:FDAPODQ3]
            if w == 0
                window = NullWindow(taq)
            elseif w == 1
                window = SineWindow(q1, q2, q3, taq)
            elseif w == 2
                window = ExponentialWindow(q1, taq)
            elseif w == 3
                window = GaussWindow(q1, q2, q3, taq)
            else
                window = UnknownWindow(taq)
            end
            d[:window] = window
        end

        push!(axesmd, d)
    end

    return md, axesmd
end

"""
    loadnmrpipe1d(filename, md, mdax)

Return an NMRData array containing spectrum and associated metadata.
"""
function loadnmrpipe1d(filename::String, md, mdax)
    npoints = mdax[1][:npoints]
    # preallocate data (and dummy header)
    header = zeros(Float32, 512)
    y = zeros(Float32, npoints)

    # read the file
    open(filename) do f
        read!(f, header)
        read!(f, y)
    end
    y = Float64.(y)

    valx = mdax[1][:val]
    delete!(mdax[1],:val) # remove values from metadata to prevent confusion when slicing up
    xaxis = X(valx, metadata=mdax[1])

    NMRData(y, (xaxis, ), metadata=md)
end

"""
    loadnmrpipe2d(filename, md, mdax)

Return NMRData containing spectrum and associated metadata.
"""
function loadnmrpipe2d(filename::String, md, mdax)
    npoints1 = mdax[1][:npoints]
    npoints2 = mdax[2][:npoints]
    transposeflag = md[:FDTRANSPOSED] > 0
    dimorder = md[:FDDIMORDER]

    # preallocate data (and dummy header)
    header = zeros(Float32, 512)
    if transposeflag
        y = zeros(Float32, npoints2, npoints1)
    else
        y = zeros(Float32, npoints1, npoints2)
    end

    # read the file
    open(filename) do f
        read!(f, header)
        read!(f, y)
    end
    if transposeflag
        y = transpose(y)
    end
    y = Float64.(y)

    valx = mdax[1][:val]
    delete!(mdax[1],:val) # remove values from metadata to prevent confusion when slicing up
    valy = mdax[2][:val]
    delete!(mdax[2],:val) # remove values from metadata to prevent confusion when slicing up
    ax2 = mdax[2][:pseudodim] ? Ti : Y
    xaxis = X(valx, metadata=mdax[1])
    yaxis = ax2(valy, metadata=mdax[2])

    NMRData(y, (xaxis, yaxis), metadata=md)
end

"""
    loadnmrpipe3d(filename, md, mdax)

Return NMRData containing spectrum and associated metadata.
"""
function loadnmrpipe3d(filename::String, md, mdax)
    datasize = md[:FDSIZE][1:3]

    # load data
    header = zeros(Float32, 512)
    y = zeros(Float32, datasize...)
    if md[:FDPIPEFLAG] == 0
        # series of 2D files
        y2d = zeros(Float32, datasize[1], datasize[2])
        for i = 1:datasize[3]
            filename1 = expandpipetemplate(filename, i)
            open(filename1) do f
                read!(f, header)
                read!(f, y2d)
            end
            y[:,:,i] = y2d
        end
    else
        # single stream
        open(filename) do f
            read!(f, header)
            read!(f, y)
        end
    end
    y = Float64.(y)

    # get axis values
    val1 = mdax[1][:val]
    val2 = mdax[2][:val]
    val3 = mdax[3][:val]
    delete!(mdax[1],:val) # remove values from metadata to prevent confusion when slicing up
    delete!(mdax[2],:val)
    delete!(mdax[3],:val)

    # y is currently in order of FDSIZE - we need to rearrange
    dimorder = md[:FDDIMORDER][1:3]
    tr = md[:FDTRANSPOSED]

    # figure out current ordering of data matrix
    # e.g. order = [2 1 3] indicates the data matrix is axis2 x axis1 x axis3
    if tr == 0
        order = [findfirst(x->x.==i, dimorder) for i=[2,1,3]]
    else # transpose flag
        # swap first two entries of dimorder round
        dimorder2 = [dimorder[2], dimorder[1], dimorder[3]]
        order = [findfirst(x->x.==i, dimorder2) for i=[1,2,3]]
    end
    # calculate the permutation required to bring the data matrix into the order axis1 x axis2 x axis3
    unorder = [findfirst(x->x.==i, order) for i=[1,2,3]]
    y = permutedims(y, unorder)

    # finally rearrange data into a useful order - always place pseudo-dimension last - and generate axes
    # 1 is always direct, and should be placed first (i.e. 1 x x)
    # is there is a pseudodimension, put that last (i.e. 1 y p)
    # otherwise, return data in the order 1 2 3 = order of fid.com, inner -> outer loop
    pdim = [mdax[i][:pseudodim] for i in 1:3]
    if pdim[2]
        # dimensions are x p y => we want ordering 1 3 2
        xaxis = X(val1, metadata=mdax[1])
        yaxis = Y(val3, metadata=mdax[3])
        zaxis = Ti(val2, metadata=mdax[2])
        y = permutedims(y, [1, 3, 2])
    elseif pdim[3]
        # dimensions are x y p => we want ordering 1 2 3
        xaxis = X(val1, metadata=mdax[1])
        yaxis = Y(val2, metadata=mdax[2])
        zaxis = Ti(val3, metadata=mdax[3])
    else
        # no pseudodimension, use Z axis not Ti
        # dimensions are x y z => we want ordering 1 2 3
        xaxis = X(val1, metadata=mdax[1])
        yaxis = Y(val2, metadata=mdax[2])
        zaxis = Z(val3, metadata=mdax[3])
    end

    NMRData(y, (xaxis, yaxis, zaxis), metadata=md)
end

"""
    metadatahelp(entry::Symbol)

Get a brief description of what's represented by metadata entry.

# Examples

```julia
metadatahelp(:FDQUADFLAG);
```
"""
function metadatahelp(entry::Symbol)
    refdict = Dict(
        # spectrum parameters
        :ndim => "Number of dimensions",
        :format => "Original file format, :NMRPipe or :bruker",
        :acqus => "Dictionary containing parsed acqus file",
        :ns => "Number of scans (NS) from acqus file",
        :rg => "Receiver gain (RG) from acqus file",
        :pulseprogram => "Pulse program (PULPROG) from acqus file",
        :filename => "Original filename or template",
        :acqusfilename => "Path to imported acqus file",
        :noise => "rms noise level",

        # dimension parameters
        :pseudodim => "Flag indicating non-frequency domain data",
        :npoints => "Number of (real) data points in dimension",
        :td => "Number of complex points acquired",
        :tdzf => "Number of complex points when FT executed, including linear prediction and zero filling",
        :bf => "Base frequency, in MHz",
        :swhz => "Spectrum width, in Hz",
        :swppm => "Spectrum width, in ppm",
        :offsetppm => "Carrier offset from bf, in ppm",
        :offsethz => "Carrier offset from bf, in Hz",
        :sf => "Carrier frequency, in MHz",
        :region => "Extracted region, expressed as a range in points, otherwise missing",
        :window => "WindowFunction structure indicating applied apodization",

        # NMRPipe spectrum parameters
        :FDMAGIC => "Should be zero in valid NMRPipe data",
        :FDFLTFORMAT => "Constant defining floating point format",
        :FDFLTORDER => "Constant defining byte order",
        :FDSIZE => "Number of points in current dim R|I",
        :FDSPECNUM => "Number of complex 1D slices in file",
        :FDQUADFLAG => "0=quad/complex, 1=real, 2=pseudoquad, 3=se, 4=grad",
        :FD2DPHASE => "0=magnitude, 1=tppi, 2=states,3=image, 4=array",
        :FDTRANSPOSED => "1=Transposed, 0=Not Transposed",
        :FDDIMCOUNT => "Number of dimensions in complete data",
        :FDDIMORDER => "Dimension stored in X/Y/Z/A axes",
        :FDFILECOUNT => "Number of files in complete data",
        :FDPIPEFLAG => "Dimension code of data stream: 0=single or multi-file; 2 1 3 or 4 = CUR_XDIM of a datastream",
        :FDCUBEFLAG => "1=3D cube multi-file series (4D only); 0 otherwise",

        # axis parameters
        :FDSIZE => "Actual size of data matrix (dim 1 = FDSIZE, dim 2 = FDSPECNUM, dim 3/4 = FDFxSIZE)"
        :FDAPOD => "Current valid time-domain size (complex points, before ZF, after extraction)",
        :FDSW => "Sweep Width Hz (after extraction)",
        :FDOBS => "Obs Freq MHz",
        :FDOBSMID => "Original Obs Freq before 0.0ppm adjust",
        :FDORIG => "Axis Origin (Last Point, After Extraction), Hz",
        :FDUNITS => "Axis Units (1=sec, 2=Hz, 3=ppm, 4=pts)",
        #:FDQUADFLAG => defined above
        :FDFTFLAG => "1=Freq Domain 0=Time Domain",
        :FDCAR => "Carrier Position, PPM",
        :FDCENTER => "Point Location of Zero Freq",
        :FDOFFPPM => "Additional PPM offset (for alignment)",
        :FDP0 => "Zero Order Phase, Degrees",
        :FDP1 => "First  Order Phase, Degrees",
        :FDAPODCODE => "Window function used (0=none, 1=SP, 2=EM, 3=GM, 4=TM, 5=ZE, 6=TRI)",
        :FDAPODQ1 => "Window parameter 1",
        :FDAPODQ2 => "Window parameter 2",
        :FDAPODQ3 => "Window parameter 3",
        :FDLB => "Extra Exponential Broadening, Hz",
        :FDGB => "Extra Gaussian Broadening, Hz",
        :FDGOFF => "Offset for Gaussian Broadening, 0 to 1",
        :FDC1 => "Add 1.0 to get First Point Scale",
        :FDZF => "Negative of Zero Fill Size",
        :FDX1 => "Extract region origin, pts (if not zero)",
        :FDXN => "Extract region endpoint, pts (if not zero)",
        :FDFTSIZE => "Size of data when FT performed",
        :FDTDSIZE => "Original valid time-domain size"
    )
    entry in keys(refdict) || throw(NMRToolsException("symbol :$(entry) not found in reference dictionary"))
    return refdict[entry]
end

function parseacqus(acqusfilename::String)
    ispath(acqusfilename) || throw(NMRToolsException("loading acqus data: $(acqusfilename) is not a valid file"))

    dat = open(acqusfilename) do f
        read(f, String)
    end

    fields = [strip(strip(x),'$') for x in split(dat, "##")]
    dic = Dict{String,Any}()
    for field in fields
        field == "" && continue # skip empty lines
        x = split(field, "= ")
        if length(x) > 1
            # store keys as uppercase - there's only a couple of annoying exceptions like NusAMOUNT
            # and forcing uppercase makes it easier to access, e.g. as :vclist not :VCLIST
            dic[uppercase(x[1])] = parseacqusentry(x[2])
        end
    end

    # lastly, check for referenced files like vclist, fq1list, and load these in place of filename
    parseacqusauxfiles!(dic, dirname(acqusfilename))

    return dic
end

function parseacqusentry(dat)
    if dat[1] == '('
        # array data - split into fields
        fields = split(split(dat, ")\n")[2])
        parsed = map(x -> parseacqusfield(x), fields)
        if parsed isa Vector{Real}
            parsed = float.(parsed)
        end
        # convert to a zero-based array to match expected bruker notation
        parsed = OffsetArray(parsed, 0:length(parsed)-1)
    else
        parsed = parseacqusfield(dat)
    end
    return parsed
end

function parseacqusfield(dat)
    if dat[1] == '<' # <string>
        dat = dat[2:end-1]
    else
        x = tryparse(Int64, dat)
        if isnothing(x)
            dat = tryparse(Float64, dat)
        else
            dat = x
        end
    end
    return dat
end

function parseacqusauxfiles!(dic, basedir)
    for k in ("VCLIST",) # integer lists
        dic[k] == "" && continue
        # note that filenames aren't actually used - vclist is always stored as vclist
        filename = joinpath(basedir, lowercase(k))
        ispath(filename) || continue
        x = readlines(filename)
        xi = tryparse.(Int,x)
        if any(xi.==nothing)
            @warn "Unable to parse format of list $(k) when opening $(filename)"
            dic[k] = x
        else
            dic[k] = xi
        end
    end

    for k in ("FQ1LIST", "FQ2LIST", "FQ3LIST", "FQ4LIST", "FQ5LIST", "FQ6LIST", "FQ7LIST", "FQ8LIST",
                "VALIST", "VDLIST", "VPLIST", "VTLIST") # float lists
        dic[k] == "" && continue
        # note that filenames aren't actually used - vclist is always stored as vclist
        filename = joinpath(basedir, lowercase(k))
        ispath(filename) || continue
        x = readlines(filename)
        xp = [replace(replace(xs, "u"=>"e-6"),"m"=>"e-3") for xs in x]
        xf = tryparse.(Float64,xp)
        if any(xf.==nothing)
            @warn "Unable to parse format of list $(k) when opening $(filename)"
            dic[k] = x
        else
            dic[k] = xf
        end
    end
end


text:
function part1(input)
    reindeers = parse_line.(eachline(input))
    return foldl(max, flying_distance.(reindeers, 2503))
end

function part2(input)
    reindeers = parse_line.(eachline(input))
    return maximum(scores(reindeers, 2503))
end

parse_line(line) = parse.(Int, split(line)[[4, 7, 14]])

function flying_distance(reindeer, t)
    v, t1, t2 = reindeer
    return v * (t1 * (t ÷ (t1 + t2)) + min(t1, t % (t1 + t2)))
end

function scores(reindeers, t)
    n = length(reindeers)
    positions = fill(0, n)
    scores = fill(0, n)
    for t = 1:t
        for i = 1:n
            v, t1, t2 = reindeers[i]
            if mod1(t, t1 + t2) <= t1
                positions[i] += v
            end
        end
        scores += positions .== maximum(positions)
    end
    return scores
end


text:
# TODO: documment all these methods

degree(g::Graph{N,E,D}, n::Integer) where {N,E,D} = ptr_length(g.nodes[n], dir_xor(D, Both))
degree(g::Graph) = [degree(g, i) for i = nodes(g)]

indegree(g::Graph) = [indegree(g, i) for i = nodes(g)]
indegree(g::Graph{N,E,D}, n::Integer) where {N,E,D} = ptr_length(g.nodes[n], dir_xor(D, Reverse))

outdegree(g::Graph{N,E,D}, n::Integer) where {N,E,D} = ptr_length(g.nodes[n], dir_xor(D, Forward))
outdegree(g::Graph) = [outdegree(g, i) for i = nodes(g)]

"""
    is_isolate(g::Graph, n::Integer)

Is node `n` an isolate, i.e. not connected to any other node?
"""
is_isolate(g::Graph, n::Integer) = degree(g, n) == 0
isolates(g::Graph) = [i for i = nodes(g) if is_isolate(g, i)]
isolatecount(g::Graph) = count(i -> is_isolate(g, i), nodes(g))

"""
    is_leaf(g::Graph, n::Integer)

Is node `n` a leaf node (i.e. a source or sink)?
See the [reference](https://en.wikipedia.org/wiki/Vertex_(graph_theory)#Types_of_vertices)
for more information.
"""
is_leaf(g::Graph, n::Integer) = degree(g, n) == 1
leaves(g::Graph) = [i for i = nodes(g) if is_leaf(g, i)]
leafcount(g::Graph) = count(i -> is_leaf(g, i), nodes(g))

"""
    is_source(g::Graph, n::Integer)

Is node `n` a source? Sources can only have outgoing edges.
In undirected graphs, there can't be any sources.
"""
is_source(g::DirectedGraph, n::Integer)  = indegree(g, n) == 0 && outdegree(g, n) > 0
is_source(g::UndirectedGraph, ::Integer) = false
sources(g::Graph) = [i for i = nodes(g) if is_source(g, i)]
sourcecount(g::Graph) = count(i -> is_source(g, i), nodes(g))

"""
    is_sink(g::Graph, n::Integer)

Is node `n` a sink? Sinks can only have incoming edges.
In undirected graphs, there can't be any sinks.
"""
is_sink(g::DirectedGraph, n::Integer)  = outdegree(g, n) == 0 && indegree(g, n) > 0
is_sink(g::UndirectedGraph, ::Integer) = false
sinks(g::Graph) = [i for i = nodes(g) if is_sink(g, i)]
sinkcount(g::Graph) = count(i -> is_sink(g, i), nodes(g))


text:
cond1 = make_condition_dict(num_nets=5,mcts_d=[8],seed=5,role=:trusted)
cond2 = make_condition_dict(num_nets=250,mcts_d=[3],seed=4,role=:candidate)
cond3 = make_condition_dict(num_nets=250,mcts_d=[1],seed=3,role=:candidate)

exp_dict = Dict(:name=>"x4_test", :conditions=>Dict(1=>cond1, 2=>cond2, 3=> cond3))

fldr = "logs"
for x in keys(exp_dict[:conditions])
    fname = make_useful_fname(exp_dict,x)
    exp_dict[:conditions][x][:fname] = fname
    exp_dict[:conditions][x][:fldr] = fldr
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using Test, LinearAlgebra, STMO

include("helpers.jl")
include("quadratic.jl")
include("unconstrained.jl")
include("testfuns.jl")
include("optimaltransport.jl")
include("mst.jl")

max_stars_repo_name: mferrodo/STMO

max_issues_repo_name: mferrodo/STMO

max_forks_repo_name: mferrodo/STMO