text:
"""
$(TYPEDEF)

A complex facet as part to the input to TetGen.

$(TYPEDFIELDS)
"""
struct RawFacet{T}
    """
    Polygons given as arrays of indices which point into the
    pointlist array describing the input points.
    """
    polygonlist::Array{Array{Cint,1},1}

    """
    Array of points given by their coordinates
    marking polygons describing holes in the facet.
    """
    holelist::Array{T,2}
end

"""
$(TYPEDEF)

A  structure for  transferring  data  into and  out  of TetGen's
internal representation.

The input of TetGen is either a 3D point set, or a 3D piecewise linear
complex (PLC), or  a tetrahedral mesh.  Depending on  the input object
and the specified  options, the output of TetGen is  either a Delaunay
(or weighted Delaunay) tetrahedralization, or a constrained (Delaunay)
tetrahedralization, or a quality tetrahedral mesh.

A piecewise  linear complex  (PLC) represents  a 3D  polyhedral domain
with  possibly internal  boundaries(subdomains). It  is introduced  in
[Miller  et al,  1996].   Basically  it is  a  set  of "cells",  i.e.,
vertices, edges, polygons, and polyhedra,  and the intersection of any
two of its cells is the union of other cells of it.

The 'RawTetGenIO' structure  is a collection of arrays  of data, i.e.,
points, facets, tetrahedra,  and so forth. All data  are compatible to
the representation in C++ and can be used without copying.

$(TYPEDFIELDS)
"""
mutable struct RawTetGenIO{T}
    """
    'pointlist':  Array of point coordinates with `size(pointlist,1)==3`.
    """
    pointlist::Array{T,2}

    """
    'pointattributelist':  Array of point attributes. The number of
     attributes per point is determined by  `size(pointattributelist,1)`
    """
    pointattributelist::Array{T,2}

    """
    'pointmtrlist': An array of metric tensors at points.
    """
    pointmtrlist::Array{T,2}

    """
    'pointmarkerlist':  An array of point markers; one integer per point.
    """
    pointmarkerlist::Array{Cint,1}

    """ 'tetrahedronlist': An array of tetrahedron corners represented
    by indices of points in `pointlist`. Unless option `-o2` is given,
    one has  `size(tetrahedronlist,1)==4`, i.e. each  column describes
    the    four     corners    of    a     tetrahedron.     Otherwise,
    `size(tetrahedronlist,1)==10` and the 4  corners are followed by 6
    edge midpoints.
    """
    tetrahedronlist::Array{Cint,2}

    """
    'tetrahedronattributelist':  An array of tetrahedron attributes.
    """
    tetrahedronattributelist::Array{T,2}

    """
    'tetrahedronvolumelist':  An array of constraints, i.e. tetrahedron's
     volume;  Input only. This can be used for triggering local refinement.
    """
    tetrahedronvolumelist::Array{T,1}

    """
    'neighborlist':  An array of tetrahedron neighbors; 4 ints per element.
     Output only.
    """
    neighborlist::Array{Cint,2}

    """
    'facetlist':  An array of facets.  Each entry is a structure of facet.
    """
    facetlist::Array{RawFacet{T},1}

    """
    'facetmarkerlist':  An array of facet markers; one int per facet.
    """
    facetmarkerlist::Array{Cint,1}

    """
    'holelist':  An array of holes (in volume).  Each hole is given by a
     point which lies strictly inside it.
    """
    holelist::Array{T,2}

    """
     'regionlist': An array of regions (subdomains).  Each region is given by
     a seed (point) which lies strictly inside it. For each column,
     the point coordinates ade  at indices [1], [2] and [3], followed by the regional
     attribute at index [4], followed by the maximum volume at index [5].

     ote that each regional attribute is used only if you select the 'A'
     switch, and each volume constraint is used only if you select the
     'a' switch (with no number following).
    """
    regionlist::Array{T,2}

    """
    'facetconstraintlist':  An array of facet constraints.  Each constraint
    specifies a maximum area bound on the subfaces of that facet.  Each column
    contains  a facet marker at index [1] and its
    maximum area bound at index [2]. Note: the facet marker is
    actually an integer.
    """
    facetconstraintlist::Array{T,2}

    """
    'segmentconstraintlist': An array of segment constraints. Each constraint
    specifies a maximum length bound on the subsegments of that segment.
    Each columb consists of the  two endpoints of the segment at
    index [1] and [2], and the maximum length bound at index [3].
    Note the segment endpoints are actually integers.
    """
    segmentconstraintlist::Array{T,2}

    """
    'trifacelist':  An array of face (triangle) corners.
    """
    trifacelist::Array{Cint,2}

    """
    'trifacemarkerlist':  An array of face markers; one int per face.
    """
    trifacemarkerlist::Array{Cint,1}

    """
    'edgelist':  An array of edge endpoints.
    """
    edgelist::Array{Cint,2}

    """
    'edgemarkerlist':  An array of edge markers.
    """
    edgemarkerlist::Array{Cint}
end

"""
$(TYPEDSIGNATURES)

Create RawTetGenIO structure with empty data.
"""
function RawTetGenIO{T}(;
                        pointlist=Array{T,2}(undef,0,0),
                        pointattributelist=Array{T,2}(undef,0,0),
                        pointmtrlist=Array{T,2}(undef,0,0),
                        pointmarkerlist=Array{Cint,1}(undef,0),
                        tetrahedronlist=Array{Cint,2}(undef,0,0),
                        tetrahedronattributelist=Array{T,2}(undef,0,0),
                        tetrahedronvolumelist=Array{T,1}(undef,0),
                        neighborlist=Array{Cint,2}(undef,0,0),
                        facetlist=Array{RawFacet{T},1}(undef,0),
                        facetmarkerlist=Array{Cint,1}(undef,0),
                        holelist=Array{T,2}(undef,0,0),
                        regionlist=Array{T,2}(undef,0,0),
                        facetconstraintlist=Array{T,2}(undef,0,0),
                        segmentconstraintlist=Array{T,2}(undef,0,0),
                        trifacelist=Array{Cint,2}(undef,0,0),
                        trifacemarkerlist=Array{Cint,1}(undef,0),
                        edgelist=Array{Cint,2}(undef,0,0),
                        edgemarkerlist=Array{Cint}(undef,0)
                        ) where T
    RawTetGenIO{T}(
        pointlist,
        pointattributelist,
        pointmtrlist,
        pointmarkerlist,
        tetrahedronlist,
        tetrahedronattributelist,
        tetrahedronvolumelist,
        neighborlist,
        facetlist,
        facetmarkerlist,
        holelist,
        regionlist,
        facetconstraintlist,
        segmentconstraintlist,
        trifacelist,
        trifacemarkerlist,
        edgelist,
        edgemarkerlist
    )
end

function Base.show(io::IO, tio::RawTetGenIO)
    nonempty(a)=size(a,ndims(a))>0
    println(io,"RawTetGenIO(")
    print(io,"numberofpoints=")
    print(io,numberofpoints(tio))
    println(io,",")
    print(io,"numberofedges=")
    print(io,numberofedges(tio))
    println(io,",")
    print(io,"numberoftrifaces=")
    print(io,numberoftrifaces(tio))
    println(io,",")
    print(io,"numberoftetrahedra=")
    print(io,numberoftetrahedra(tio))
    println(io,",")
    for name in fieldnames(typeof(tio))
#        a=round.(getfield(tio,name),sigdigits=3)'
        a=getfield(tio,name)
        if nonempty(a)
            print(io,"$(name)'=")
            print(io,a)
            println(io,",")
        end
    end
    println(io,")")
end

"""
    $(TYPEDSIGNATURES)
    Set list of input facets from AbstractMatrix desribing polygons of the same
    size (e.g. triangles)
"""
function facetlist!(tio::RawTetGenIO{T}, facets::AbstractMatrix) where {T}
    numberoffacets=size(facets,2)
    tio.facetlist=Array{RawFacet{T},1}(undef,numberoffacets)
    for ifacet=1:numberoffacets
        tio.facetlist[ifacet]=RawFacet{T}([Vector{Cint}(facets[:,ifacet])],Array{T,2}(undef,0,0))
    end
    tio
end

"""
    $(TYPEDSIGNATURES)
    Set list of input facets from a vector of polygons of different size
"""
function facetlist!(tio::RawTetGenIO{T}, facets::Vector) where {T}
    numberoffacets=size(facets,1)
    tio.facetlist=Array{RawFacet{T},1}(undef,numberoffacets)
    for ifacet=1:numberoffacets
        tio.facetlist[ifacet]=RawFacet{T}([Vector{Cint}(facets[ifacet])],Array{T,2}(undef,0,0))
    end
    tio
end

#
# Create CPPTetGenIO from RawTetGenIO
#
function CPPTetGenIO(tio::RawTetGenIO{T}) where T

    firstnumber = 1
    mesh_dim = 3

    # Set dummy defaults
    pointlist = C_NULL
    pointattributelist = C_NULL
    pointmtrlist = C_NULL

    pointmarkerlist = C_NULL
    numberofpoints = zero(Cint)
    numberofpointattributes = zero(Cint)
    numberofpointmtrs = zero(Cint)

    tetrahedronlist = C_NULL
    tetrahedronattributelist = C_NULL
    tetrahedronvolumelist = C_NULL

    neighborlist = C_NULL
    numberoftetrahedra = zero(Cint)
    numberofcorners = zero(Cint)
    numberoftetrahedronattributes = zero(Cint)

    facetlist =  C_NULL
    facetmarkerlist = C_NULL
    numberoffacets = zero(Cint)

    holelist = C_NULL
    numberofholes = zero(Cint)

    regionlist = C_NULL
    numberofregions = zero(Cint)

    facetconstraintlist = C_NULL
    numberoffacetconstraints = zero(Cint)

    segmentconstraintlist = C_NULL
    numberofsegmentconstraints = zero(Cint)

    trifacelist = C_NULL
    trifacemarkerlist = C_NULL
    numberoftrifaces = zero(Cint)

    edgelist = C_NULL
    edgemarkerlist = C_NULL
    numberofedges = zero(Cint)

    # these need to be returned in order
    # to prevent their content  from being garbage collected
    facetlist_array=Array{CFacet{T},1}(undef,numberoffacets)
    polygonlist_array=Array{Array{CPolygon,1},1}(undef,0)

    # Override defaults
    numberofpoints=size(tio.pointlist,2)
    @assert numberofpoints>zero(Cint)
    @assert size(tio.pointlist,1) == 3
    pointlist=pointer(tio.pointlist)

    numberofpointattributes=size(tio.pointattributelist,1)
    if numberofpointattributes>0
        @assert size(tio.pointattributelist,2)==numberofpoints
        pointattributelist=pointer(tio.pointattributelist)
    end

    numberofpointmtrs=size(tio.pointmtrlist,1)
    if numberofpointmtrs>0
        @assert size(tio.pointmtrlist,2)==numberofpoints
        pointmtrlist=pointer(tio.pointmtrlist)
    end

    numberoftetrahedra=size(tio.tetrahedronlist,2)
    if numberoftetrahedra>0
        tetrahedronlist=pointer(tio.tetrahedronlist)
        numberoftetrahedronattributes=size(tio.tetrahedronattributelist,1)
        numberofcorners=size(tio.tetrahedronlist,1)
        if numberoftetrahedronattributes>0
            @assert size(tio.tetrahedronattributelist,2)==numberoftetrahedra
            tetrahedronattributelist=pointer(tio.tetrahedronattributelist)
        end
        if size(tio.tetrahedronvolumelist,1)>0
            @assert size(tio.tetrahedronvolumelist,1)==numberoftetrahedra
            tetrahedronvolumelist=pointer(tio.tetrahedronvolumelist)
        end
    end

    numberoffacets=size(tio.facetlist,1)
    if numberoffacets>0
        for ifacet=1:numberoffacets
            facet=tio.facetlist[ifacet]
            numberofpolygons=size(facet.polygonlist,1)
            polygonlist=Array{CPolygon,1}(undef, numberofpolygons)
            for ipolygon=1:numberofpolygons
                polygonlist[ipolygon]=CPolygon(pointer(facet.polygonlist[ipolygon]),size(facet.polygonlist[ipolygon],1))
            end
            numberofholes=size(facet.holelist,2)
            push!(polygonlist_array,polygonlist)
            push!(facetlist_array,CFacet(pointer(polygonlist),Cint(numberofpolygons),pointer(facet.holelist),Cint(numberofholes)))
        end
        facetlist=pointer(facetlist_array)
    end

    if size(tio.facetmarkerlist,1)>0
        @assert size(tio.facetmarkerlist,1)==numberoffacets
        facetmarkerlist=pointer(tio.facetmarkerlist)
    end

    numberofholes=size(tio.holelist,2)
    if numberofholes>0
        @assert size(tio.holelist,1)==3
        holelist=pointer(tio.holelist)
    end

    numberofregions=size(tio.regionlist,2)
    if numberofregions>0
        @assert size(tio.regionlist,1)==5 # coord+attribute + volume
        regionlist=pointer(tio.regionlist)
    end

    if size(tio.facetconstraintlist,2)>0
        @assert size(tio.facetconstraintlist,1) == 2
        numberoffacetconstraints=size(tio.facetconstraintlist,2)
        facetconstraintlist=pointer(tio.facetconstraintlist)
    end

    if size(tio.segmentconstraintlist,2)>0
        @assert size(tio.segmentconstraintlist,1) == 3
        numberofsegmentconstraints=size(tio.segmentconstraintlist,2)
        segmentconstraintlist=pointer(tio.segmentconstraintlist)
    end

    if size(tio.trifacelist,2)>0
        @assert size(tio.trifacelist,1) == 3
        numberoftrifaces=size(tio.trifacelist,2)
        trifacelist=pointer(tio.trifacelist)
        if size(tio.trifacemarkerlist,1) > 0
            @assert size(tio.trifacemarkerlist,1) == numberoftrifaces
            trifacemarkerlist=pointer(tio.trifacemarkerlist)
        end
    end

    if size(tio.edgelist,2)>0
        @assert size(tio.edgelist,1) == 2
        numberofedges=size(tio.edgelist,2)
        edgelist=pointer(tio.edgelist)
        if size(tio.edgemarkerlist,1) > 0
            @assert size(tio.edgemarkerlist,1) == numberofedges
            edgemarkerlist=pointer(tio.edgemarkerlist)
        end
    end

    # Create struct
    CPPTetGenIO{T}(firstnumber,
                   mesh_dim,
                   pointlist,
                   pointattributelist,
                   pointmtrlist,
                   pointmarkerlist,
                   numberofpoints,
                   numberofpointattributes,
                   numberofpointmtrs,
                   tetrahedronlist,
                   tetrahedronattributelist,
                   tetrahedronvolumelist,
                   neighborlist,
                   numberoftetrahedra,
                   numberofcorners,
                   numberoftetrahedronattributes,
                   facetlist,
                   facetmarkerlist,
                   numberoffacets,
                   holelist,
                   numberofholes,
                   regionlist,
                   numberofregions,
                   facetconstraintlist,
                   numberoffacetconstraints,
                   segmentconstraintlist,
                   numberofsegmentconstraints,
                   trifacelist,
                   trifacemarkerlist,
                   numberoftrifaces,
                   edgelist,
                   edgemarkerlist,
                   numberofedges),facetlist_array, polygonlist_array
end

#
# Create  RawTetGenIO from CPPTetGenIO
#
function RawTetGenIO(ctio::CPPTetGenIO{T}) where T

    tio=RawTetGenIO{T}()
    if ctio.numberofpoints>0 && ctio.pointlist!=C_NULL
        tio.pointlist = convert(Array{T,2}, Base.unsafe_wrap(Array, ctio.pointlist, (3,Int(ctio.numberofpoints)), own=true))
    end
    if ctio.numberofpointattributes>0  && ctio.pointattributelist!=C_NULL
        tio.pointattributelist=convert(Array{T,2}, Base.unsafe_wrap(Array, ctio.pointattributelist, (Int(ctio.numberofpointattributes),Int(ctio.numberofpoints)), own=true))
    end
    if ctio.numberofpointmtrs>0  && ctio.pointmtrlist!=C_NULL
        tio.pointmtrlist=convert(Array{T,2}, Base.unsafe_wrap(Array, ctio.pointmtrlist, (Int(ctio.numberofpointmtrs),Int(ctio.numberofpoints)), own=true))
    end
    if ctio.pointmarkerlist!=C_NULL
        tio.pointmarkerlist=convert(Array{Cint,1}, Base.unsafe_wrap(Array, ctio.pointmarkerlist, (Int(ctio.numberofpoints)), own=true))
    end
    @assert ctio.numberofcorners==4
    if ctio.numberoftetrahedra>0 && ctio.tetrahedronlist!=C_NULL
        tio.tetrahedronlist=convert(Array{Cint,2}, Base.unsafe_wrap(Array, ctio.tetrahedronlist, (4,Int(ctio.numberoftetrahedra)), own=true))
    end
    if ctio.numberoftetrahedronattributes>0  && ctio.tetrahedronattributelist!=C_NULL
        tio.tetrahedronattributelist=convert(Array{T,2}, Base.unsafe_wrap(Array, ctio.tetrahedronattributelist, (Int(ctio.numberoftetrahedronattributes),Int(ctio.numberoftetrahedra)), own=true))
    end
    if ctio.tetrahedronvolumelist!=C_NULL
        tio.tetrahedronvolumelist=convert(Array{T,1}, Base.unsafe_wrap(Array, ctio.tetrahedronvolumelist, (Int(ctio.numberoftetrahedra)), own=true))
    end
    if ctio.numberoftetrahedra>0 && ctio.neighborlist!=C_NULL
        tio.neighborlist=convert(Array{Cint,2}, Base.unsafe_wrap(Array, ctio.neighborlist, (4,Int(ctio.numberoftetrahedra)), own=true))
    end

    # Essentially, facetlist appears to be used only during input, so we can skip the
    # conversion here (and increase test coverage...).
    # The created facets are in the trifacelist further down this code.
    # When uncommenting, test...
    # if ctio.numberoffacets>0
    #     facetlist=convert(Array{CFacet,1}, Base.unsafe_wrap(Array, ctio.facetlist, (Int(ctio.numberoffacets)), own=true))
    #     tio.facetlist=Array{RawFacet{T},1}(undef,ctio.numberoffacets)
    #     for ifacet=1:numberoffacets
    #         # Possibly this is copied from the input, so we
    #         # might not want to own the pointers here.
    #         facet=facetlist[i]
    #         polygonlist=Array{Array{Cint,1},1}(undef,facet.numberofpolygons)
    #         for ipolygon=1:facet.numberofpolygons
    #             polygonlist[ipolygon]=convert(Array{Cint,1},Base.unsafe_wrap(Array, facet.polygonlist[i].vertexlist , (Int(facet.polygonlist[i].numberofvertices)), own=true))
    #         end
    #         holelist=convert(Array{T,1},Base.unsafe_wrap(Array, facet.holelist , (Int(facet.numberofholes)), own=true))
    #         tio.facetlist=RawFacet{T}(polygonlist,holelist)
    #     end
    #     if ctio.facetmarkerlist != C_NULL
    #         tio.facetmarkerlist=convert(Array{Cint,1}, Base.unsafe_wrap(Array, ctio.facetmarkerlist, (Int(ctio.numberoffacets)), own=true))
    #     end
    # end

    # Usually, this ctio comes from the output of tetrahedralize(). In this case, the holelist pointer is copied from the input
    # so we would get a double free corruption if we own it here.
    if ctio.numberofholes>0 && ctio.holelist!=C_NULL
        tio.holelist=convert(Array{T,2}, Base.unsafe_wrap(Array, ctio.holelist, (3,Int(ctio.numberofholes)), own=false))
    end

    # Usually, this ctio comes from the output of tetrahedralize(). In this case, the regionlist pointer is copied from the input
    # so we would get a double free corruption if we own it here.
    if ctio.numberofregions>0  && ctio.regionlist!=C_NULL
        tio.regionlist=convert(Array{T,2}, Base.unsafe_wrap(Array, ctio.regionlist, (5,Int(ctio.numberofregions)), own=false))
    end

    # own ? May be comes from input
    if ctio.numberoffacetconstraints>0
        tio.facetconstraintlist=convert(Array{T,1}, Base.unsafe_wrap(Array, ctio.facetconstraintlist, (2,ctio.numberoffacetconstraints), own=false))
    end

    # own ? May be comes from input
    if ctio.numberofsegmentconstraints>0
        tio.segmentconstraintlist=convert(Array{T,1}, Base.unsafe_wrap(Array, ctio.segmentconstraintlist, (3,ctio.numberofsegmentconstraints), own=false))
    end

    if ctio.numberoftrifaces>0  && ctio.trifacelist!=C_NULL && ctio.trifacemarkerlist!=C_NULL
        tio.trifacelist=convert(Array{Cint,2}, Base.unsafe_wrap(Array, ctio.trifacelist, (3,Int(ctio.numberoftrifaces)), own=true))
        tio.trifacemarkerlist=convert(Array{Cint,1}, Base.unsafe_wrap(Array, ctio.trifacemarkerlist, (Int(ctio.numberoftrifaces)), own=true))
    end

    if ctio.numberofedges>0  && ctio.edgelist!=C_NULL && ctio.edgemarkerlist!=C_NULL
        tio.edgelist=convert(Array{Cint,2}, Base.unsafe_wrap(Array, ctio.edgelist, (2,Int(ctio.numberofedges)), own=true))
        tio.edgemarkerlist=convert(Array{Cint,1}, Base.unsafe_wrap(Array, ctio.edgemarkerlist, (Int(ctio.numberofedges)), own=true))
    end

    return tio
end

"""
$(TYPEDSIGNATURES)

Number of points in tetrahedralization
"""
numberofpoints(tio::RawTetGenIO{T}) where T= size(tio.pointlist,2)

"""
$(TYPEDSIGNATURES)

Number of tetrahedra in tetrahedralization
"""
numberoftetrahedra(tio::RawTetGenIO{T}) where T= size(tio.tetrahedronlist,2)

"""
$(TYPEDSIGNATURES)

Number of triangle faces in tetrahedralization
"""
numberoftrifaces(tio::RawTetGenIO{T}) where T= size(tio.trifacelist,2)

"""
$(TYPEDSIGNATURES)

Number of edges in tetrahedralization
"""
numberofedges(tio::RawTetGenIO{T}) where T= size(tio.edgelist,2)

"""
$(TYPEDSIGNATURES)

Tetrahedralize input.

````
  flags: -pYrq_Aa_miO_S_T_XMwcdzfenvgkJBNEFICQVh
    -p  Tetrahedralizes a piecewise linear complex (PLC).
    -Y  Preserves the input surface mesh (does not modify it).
    -r  Reconstructs a previously generated mesh.
    -q  Refines mesh (to improve mesh quality).
    -R  Mesh coarsening (to reduce the mesh elements).
    -A  Assigns attributes to tetrahedra in different regions.
    -a  Applies a maximum tetrahedron volume constraint.
    -m  Applies a mesh sizing function.
    -i  Inserts a list of additional points.
    -O  Specifies the level of mesh optimization.
    -S  Specifies maximum number of added points.
    -T  Sets a tolerance for coplanar test (default 1e-8).
    -X  Suppresses use of exact arithmetic.
    -M  No merge of coplanar facets or very close vertices.
    -w  Generates weighted Delaunay (regular) triangulation.
    -c  Retains the convex hull of the PLC.
    -d  Detects self-intersections of facets of the PLC.
    -z  Numbers all output items starting from zero.
    -f  Outputs all faces to .face file.
    -e  Outputs all edges to .edge file.
    -n  Outputs tetrahedra neighbors to .neigh file.
    -v  Outputs Voronoi diagram to files.
    -g  Outputs mesh to .mesh file for viewing by Medit.
    -k  Outputs mesh to .vtk file for viewing by Paraview.
    -J  No jettison of unused vertices from output .node file.
    -B  Suppresses output of boundary information.
    -N  Suppresses output of .node file.
    -E  Suppresses output of .ele file.
    -F  Suppresses output of .face and .edge file.
    -I  Suppresses mesh iteration numbers.
    -C  Checks the consistency of the final mesh.
    -Q  Quiet:  No terminal output except errors.
    -V  Verbose:  Detailed information, more terminal output.
    -h  Help:  A brief instruction for using TetGen.
````
"""
function tetrahedralize(input::RawTetGenIO{Float64}, flags::String)
    cinput,flist,plist=CPPTetGenIO(input)
    rc=Cint[0]
    coutput = ccall((:tetrahedralize2_f64, libtet), CPPTetGenIO{Float64}, (CPPTetGenIO{Float64}, Cstring, Ptr{Cint}),
                    cinput, flags, rc)
    if rc[1]!=0
        throw(TetGenError(rc[1]))
    end
    RawTetGenIO(coutput)
end

"""
$(TYPEDSIGNATURES)

Create GeometryBasics.Mesh from the triface list
(for quick visualization purposes using Makie's wireframe).
"""
function surfacemesh(tgio::RawTetGenIO)
    points=[ Point3f0(tgio.pointlist[:,i]...) for i=1:size(tgio.pointlist,2)]
    faces= [TriangleFace(tgio.trifacelist[:,i]...) for i=1:size(tgio.trifacelist,2)]
    mesh=GeometryBasics.Mesh(points,faces)
end

"""
$(TYPEDSIGNATURES)

Create GeometryBasics.Mesh of all tetrahedron faces
(for quick visualization purposes using Makie's wireframe).
"""
function volumemesh(tgio::RawTetGenIO)
    points=[ Point3f0(tgio.pointlist[:,i]...) for i=1:size(tgio.pointlist,2)]
    faces=Array{NgonFace{3,Int32},1}(undef,0)
    tetlist=tgio.tetrahedronlist
    for itet=1:size(tetlist,2)
        tet=view(tetlist,:, itet)
        push!(faces,TriangleFace(tet[1],tet[2],tet[3]))
        push!(faces,TriangleFace(tet[1],tet[2],tet[4]))
        push!(faces,TriangleFace(tet[2],tet[3],tet[4]))
        push!(faces,TriangleFace(tet[3],tet[1],tet[4]))
    end
    mesh=GeometryBasics.Mesh(points,faces)
end


text:
@noinline nonatomic_field_error(name::Symbol) =
    ConcurrencyViolationError("property $name cannot be accessed atomically")

@inline invalid_ordering(name::Symbol) =
    ConcurrencyViolationError("invalid ordering: $name")

let body = foldr(_JULIA_ORDERINGS, init = :(invalid_ordering(order))) do name, ex
        quote
            if order === $(QuoteNode(name))
                f(UnsafeAtomics.$name)
            else
                $ex
            end
        end
    end

    @eval @inline function with_llvm_order_from_julia_order(f, order::Symbol)
        $body
    end
end

@inline function Base.getproperty(ref::AtomicRef, name::Symbol, order::Symbol)
    if name === :_
        @inline _get(order::Ordering) = UnsafeAtomics.load(ref, order)
        return with_llvm_order_from_julia_order(_get, order)
    end
    nonatomic_field_error(name)
end

@inline function Base.setproperty!(ref::AtomicRef, name::Symbol, x, order::Symbol)
    if name === :_
        @inline _set!(order::Ordering) = UnsafeAtomics.store!(ref, x, order)
        return with_llvm_order_from_julia_order(_set!, order)
    end
    nonatomic_field_error(name)
end

let body = foldr(_JULIA_ORDERINGS, init = :(invalid_ordering(order))) do succ, ex1
        ex2 = foldr(_JULIA_ORDERINGS, init = :(invalid_ordering(order))) do fail, ex
            quote
                if fail_order === $(QuoteNode(fail))
                    return UnsafeAtomics.cas!(
                        ref,
                        expected,
                        desired,
                        UnsafeAtomics.$succ,
                        UnsafeAtomics.$fail,
                    )
                else
                    $ex
                end
            end
        end
        quote
            if success_order === $(QuoteNode(succ))
                $ex2
            else
                $ex1
            end
        end
    end

    @eval @inline function Base.replaceproperty!(
        ref::AtomicRef,
        name::Symbol,
        expected,
        desired,
        success_order::Symbol,
        fail_order::Symbol,
    )
        if name === :_
            $body
        end
        nonatomic_field_error(name)
    end
end

@inline function Base.modifyproperty!(
    ref::AtomicRef,
    name::Symbol,
    op::OP,
    x,
    order::Symbol,
) where {OP}
    if name === :_
        @inline _modify!(order::Ordering) = UnsafeAtomics.modify!(ref, op, x, order)
        return with_llvm_order_from_julia_order(_modify!, order)
    end
    nonatomic_field_error(name)
end

@inline Base.swapproperty!(ref::AtomicRef, name::Symbol, x, order::Symbol) =
    first(modifyproperty!(ref, name, right, x, order))

ext: jl


text:
module PrecompileFuncs

using TuringDash
using TuringDash.Turing
using TuringDash.Random
using TuringDash.CSV, TuringDash.DataFrames

function precompile_turing()
    chn = Base.Channel{Vector{Float64}}(Inf)
    df1 = nothing
    global df1
    df1 = CSV.read(download("https://raw.githubusercontent.com/efmanu/TuringDashApp.jl/master/TuringDash/datasets/data1.csv"), DataFrame)
    model_str1 = """
    @model regression_model(x, y) = begin
        a ~ Normal()
        b ~ Normal()
        for i in 1:10
            y[i] ~ Normal(a + b * x[i], 1.0)
        end
    end
    """
    model_str2 = "regression_model(x, y)"
    for colname in Symbol.(names(df1))
        @eval $colname = df1.$colname
    end

    eval(Meta.parse(model_str1))
    model = eval(Meta.parse(model_str2))

    alg = MH(
        [0.25 0.05;
         0.05 0.50]
    )
    rng = Random.GLOBAL_RNG

    nperiteration = 50
    nsamples = 20

    r = sample(rng, model, alg, nperiteration; chain_type=MCMCChains.Chains, save_state=true, progress=false)
    put!(chn, Array(r)[end,1:2])
    for i in (nperiteration + 1):nperiteration:nsamples
        r = Turing.Inference.resume(r, nperiteration, save_state=true, progress=false)
        put!(chn, Array(r)[end,1:2])
    end
    return "success"
end

end


text:
module JpegTurbo

include("../lib/LibJpeg.jl")
using .LibJpeg
include("libjpeg_utils.jl")

end


text:
module Dagger

using Distributed, SharedArrays

using MemPool

import Base: collect, adjoint, reduce
import Distributed: procs

using LinearAlgebra
import LinearAlgebra: transpose

using UUIDs

import ContextVariablesX

using Requires

const PLUGINS = Dict{Symbol,Any}()
const PLUGIN_CONFIGS = Dict{Symbol,String}(
    :scheduler => "Dagger.Sch"
)

include("lib/util.jl")
include("lib/logging.jl")

# Distributed data
include("options.jl")
include("processor.jl")
include("scopes.jl")
include("thunk.jl")
include("chunks.jl")

# Task scheduling
include("compute.jl")
include("utils/clock.jl")
include("utils/system_uuid.jl")
include("sch/Sch.jl"); using .Sch

# Array computations
include("array/darray.jl")
include("array/alloc.jl")
include("array/map-reduce.jl")

# File IO
include("file-io.jl")

include("array/operators.jl")
include("array/getindex.jl")
include("array/setindex.jl")
include("array/matrix.jl")
include("array/sparse_partition.jl")
include("array/sort.jl")

# Other
include("ui/graph.jl")
include("ui/gantt-common.jl")
include("ui/gantt-text.jl")

include("table/dtable.jl")
include("table/gdtable.jl")
include("table/tables.jl")
include("table/operations.jl")
include("table/groupby.jl")
include("table/join_interface.jl")
include("table/join.jl")

include("lib/logging-extras.jl")

function __init__()
    # Initialize system UUID
    system_uuid()

    @require Luxor="ae8d54c2-7ccd-5906-9d76-62fc9837b5bc" begin
        # Gantt chart renderer
        include("ui/gantt-luxor.jl")
    end
    @require Mux="a975b10e-0019-58db-a62f-e48ff68538c9" begin
        # Gantt chart HTTP server
        include("ui/gantt-mux.jl")
    end
    @require ProfileSVG="132c30aa-f267-4189-9183-c8a63c7e05e6" begin
        # Profile renderer
        include("ui/profile-profilesvg.jl")
    end
    @require FFMPEG="c87230d0-a227-11e9-1b43-d7ebe4e7570a" begin
        @require FileIO="5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" begin
            # Video generator
            include("ui/video.jl")
        end
    end
    for tid in 1:Threads.nthreads()
        add_processor_callback!("__cpu_thread_$(tid)__") do
            ThreadProc(myid(), tid)
        end
    end
end

end # module


text:
AbstractAE = Union{AbstractVLAE,} # add more here
Base.Dict(nt::NamedTuple) = Dict(pairs(nt))
# for show functions
function trim(s::String)
    if length(s) > 70
        return "$(s[1:67])..."
    else
        return s
    end
end

### represents the outputs of training ###
mutable struct TrainingOutputs{AE<:Dict, D<:Dict, C<:Dict}
    autoencoder::AE
    detectors::D
    classifier::C
end
function Base.show(io::IO, to::TrainingOutputs)
    msg = """autoencoder =
        """
    ae_outs = to.autoencoder
    for (k, v) in pairs(ae_outs)
        v = trim(repr(v))
        msg = msg*"""    $k = $v
        """
    end

    msg = msg*"""detectors =
        """
    d_outs = to.detectors
    for (k, v) in pairs(d_outs)
        v = trim(repr(v))
        msg = msg*"""    $k = $v
        """
    end

    msg = msg*"""classifier =
        """
    c_outs = to.classifier
    for (k, v) in pairs(c_outs)
        v = trim(repr(v))
        msg = msg*"""    $k = $v
        """
    end

    print(io, msg)
end

### represents the input parameters ###
mutable struct TrainingParameters{AE<:Dict, D<:Dict, C<:Dict}
    autoencoder::AE
    detectors::D
    classifier::C
end
function Base.show(io::IO, tp::TrainingParameters)
    msg = """autoencoder =
        """
    ae_outs = tp.autoencoder
    for (k, v) in pairs(ae_outs)
       msg = msg*"""    $k = $v
        """
    end

    msg = msg*"""detectors =
        """
    d_outs = tp.detectors
    for (k, v) in pairs(d_outs)
       msg = msg*"""    $k = $v
        """
    end

    msg = msg*"""classifier =
        """
    c_outs = tp.classifier
    for (k, v) in pairs(c_outs)
       msg = msg*"""    $k = $v
        """
    end

    print(io, msg)
end

TrainingOutputs() = TrainingOutputs(Dict(), Dict(), Dict())
TrainingParameters() = TrainingParameters(Dict(), Dict(), Dict())

### Hierarchical Anomaly Detector ###

mutable struct HAD{AE<:AbstractAE,DS<:Vector,CL,PS<:TrainingParameters,TO<:TrainingOutputs}
    autoencoder::AE
    detectors::DS
    classifier::CL
    parameters::PS
    training_output::TO
end
function Base.show(io::IO, m::HAD)
    nl = length(m.detectors)
    ae = repr(m.autoencoder)
    cl = repr(m.classifier)
    d = repr(m.detectors[1])
    outs = repr(m.training_output)

    msg = """Hierarchical Anomaly Detector:

     autoencoder        =
        $(ae)
     detectors          =
        $(nl) x $d
     classifier          =
        $(cl)

     training_output    =
    $outs
    """

    # now parse the output
    print(io, msg)
end

"""
    HAD(
        n_latents::Int,
        autoencoder_params,
        autoencoder_constructor,
        detector_constructor,
        detector_params,
        classifier_params,
        [classifier]
        )

Params must a a Dict/NamedTuple. Default classifier is just a single Dense layer.
"""
HAD(nl::Int, ap, ac, dp, dc, cp, c) = HAD(
        ac(;ap...),
        [dc(;dp...) for _ in 1:nl],
        c,
        TrainingParameters(Dict(ap), Dict(dp), Dict(cp)),
        TrainingOutputs()
        )
function HAD(nl::Int, ap, ac, dp, dc, cp)
    # decide whether to use reconstruction error in the classifier or not
    userec = get(cp, :use_reconstruction, true)
    if userec
        cp = merge(cp, (use_reconstruction=true,))
        return HAD(nl, ap, ac, dp, dc, cp, Dense(nl+1, 2))
    else
        return HAD(nl, ap, ac, dp, dc, cp, Dense(nl, 2))
    end
end

# dense fvlae parameters
"""
zdim,
hdims,
batchsize,
discriminator_hdim,
nepochs,
λ=0.0f0, γ=1.0f0, epochsize = size(tr_x,2), layer_depth=1, lr=0.001f0,
activation="relu", discriminator_nlayers=3,
initial_convergence_threshold=0, initial_convergence_epochs=10,
max_retrain_tries=10
"""
fvlae_constructor(;zdim::Int=1, hdims=(1,), discriminator_hdim=1, datasize=(1,1), kwargs...) =
    FVLAE(zdim, hdims, discriminator_hdim, datasize; kwargs...)

# conv fvlae parameters
conv_fvlae_constructor(;zdim::Int=1, hdim=1, ks=[(1,1)], ncs=[1], strd=1, datasize=(1,1,1,1), kwargs...) =
    gpu(FVLAE(zdim, hdim, ks, ncs, strd, datasize; kwargs...))

# knn parameters
"""
k::Int=1 - number of neighbors
v::Symbol=:kappa - distance computation method, one of [:kappa, :gamma, :delta]
t=:BruteTree - tree that encodes training data, one of [:BruteTree, :KDTree, :BallTree]
"""

### FIT METHODS ###
"""
	fit_autoencoder!(model::HAD, tr_x::AbstractArray, val_x::AbstractArray)
"""
function fit_autoencoder!(model::HAD, tr_x::AbstractArray, val_x::AbstractArray)
    @info "############################# \nStarting autoencoder training"

    # train fvlae
    local args
    if haskey(model.parameters.autoencoder, :hdim) # conv fvlae
        args = (
            model.parameters.autoencoder[:zdim],
            model.parameters.autoencoder[:hdim],
            model.parameters.autoencoder[:batchsize],
            model.parameters.autoencoder[:ks],
            model.parameters.autoencoder[:ncs],
            model.parameters.autoencoder[:strd],
            model.parameters.autoencoder[:nepochs],
            tr_x,
            val_x
            )
    else
        args = (
            model.parameters.autoencoder[:zdim],
            model.parameters.autoencoder[:hdims],
            model.parameters.autoencoder[:batchsize],
            model.parameters.autoencoder[:discriminator_hdim],
            model.parameters.autoencoder[:nepochs],
            tr_x,
            val_x
            )
    end
    ae, hist, rdata, zs, aeopt, copt = train_fvlae(
        args...;
        model.parameters.autoencoder...
    )
    println("Training values at the end of training:")
	for key in keys(hist)
	    v = get(hist, key)[2][end]
	    println("    $key = $v")
	end

    # now assign the right fields
    @info "Autoencoder training finished. \n#############################\n"
    model.autoencoder = ae
    model.training_output.autoencoder = Dict((history = hist, reconstructions = rdata, latents = zs)) # opts?
end

"""
	fit_detectors!(model::HAD, tr_x::AbstractArray)
"""
function fit_detectors!(model::HAD, tr_x::AbstractArray)
    @info "############################# \nStarting detector training"

    # encode the data
    tr_e = encode_all(model.autoencoder, tr_x; batchsize=model.parameters.autoencoder[:batchsize], mean=true)

    # fit the detectors
    map(x->fit!(x[1], x[2]), zip(model.detectors, tr_e))

    # now assign the right fields
    @info "Detector training finished. \n#############################\n"
end
function fit_detectors!(model::HAD, tr_x::AbstractArray{T,4}) where T
    @info "############################# \nStarting detector training"

    # encode the data
    tr_e = cpu(encode_all(model.autoencoder, tr_x; batchsize=model.parameters.autoencoder[:batchsize], mean=true))

    # fit the detectors
    map(x->fit!(x[1], x[2]), zip(model.detectors, tr_e))

    # now assign the right fields
    @info "Detector training finished. \n#############################\n"
end

"""
	train_classifier!(classifier, tr_x, tr_y, val_x, val_y;
        verb = true, λ=0f0, batchsize=128, lr=0.001f0, nepochs=500,
        patience = 10, val_ratio=0.2)
"""
function train_classifier!(classifier, tr_x, tr_y, val_x, val_y;
        verb = true, λ=0f0, batchsize=128, lr=0.001f0, nepochs=500,
        patience = 10, val_ratio=0.2, kwargs...)
    opt = ADAM(lr)
    ps = params(classifier)
    loss(x,y) = Flux.Losses.logitcrossentropy(classifier(x), y) + λ*sum(sqnorm, ps)
    loss(x) = loss(x...)

       batchsize2 = floor(Int, batchsize/2)
    local data_itr_n, data_itr_a
    @suppress begin
        data_itr_n = Flux.Data.DataLoader(tr_x[:,tr_y[1,:]], batchsize=batchsize2)
        data_itr_a = Flux.Data.DataLoader(tr_x[:,tr_y[2,:]], batchsize=batchsize2)
    end

    _val_auc = 0
    _patience = 1
    history = MVHistory()
    for epoch in 1:nepochs
        for (xn, xa) in zip(data_itr_n, data_itr_a)
            y = Flux.onehotbatch(Bool.(vcat(zeros(Int, size(xn,2)), ones(Int, size(xa,2)))), 0:1)
            HierarchicalAD.param_update!(loss, ps, (hcat(xn,xa), y), opt)
         end
        tr_l = loss(tr_x, tr_y)
        val_l = loss(val_x, val_y)
        tr_auc = auc_val(Flux.onecold(tr_y) .-1, classifier_score(classifier, tr_x))
        val_auc = auc_val(Flux.onecold(val_y) .-1, classifier_score(classifier, val_x))
        if verb
            println("Epoch $epoch: loss = $(tr_l), tr AUC = $(tr_auc), val AUC = $(val_auc)")
        end
        push!(history, :tr_loss, epoch, tr_l)
        push!(history, :val_loss, epoch, val_l)
        push!(history, :tr_auc, epoch, tr_auc)
        push!(history, :val_auc, epoch, val_auc)

        if (val_auc < _val_auc)
            if _patience >= patience
                if verb
                    @info "Stopping early after $epoch epochs."
                end
                return history
            else
                _patience += 1
            end
        else
            _val_auc = val_auc
            #_patience = 0
        end
    end

    return history
end

"""
	fit_classifier!(model::HAD, x::AbstractArray, y::Vector)
"""
function fit_classifier!(model::HAD, x::AbstractArray, y::Vector)
    # split data
    vr = model.parameters.classifier[:val_ratio]
    _, (tr_x, tr_y), (val_x, val_y) = _split_data(x, y, vr)
    fit_classifier!(model, tr_x, tr_y, val_x, val_y)
end
"""
    fit_classifier!(model::HAD, tr_x::AbstractArray, tr_y::Vector, val_x::AbstractArray, val_y::Vector)
"""
function fit_classifier!(model::HAD, tr_x::AbstractArray, tr_y::Vector, val_x::AbstractArray, val_y::Vector)
    @info "############################# \nStarting classifier training"

    # create and scale inputs
    cl_tr_x, cl_tr_y = classifier_inputs(model, tr_x, tr_y)
    cl_val_x, cl_val_y = classifier_inputs(model, val_x, val_y)
    norm_coeffs = normalization_coefficients(hcat(cl_tr_x, cl_val_x))
    cl_tr_x = normalize(cl_tr_x, norm_coeffs...)
    cl_val_x = normalize(cl_val_x, norm_coeffs...)
    model.training_output.classifier[:normalization_coefficients] = norm_coeffs

    # train it
    # try out more classifiers
    n_candidates = model.parameters.classifier[:n_candidates]
    candidates = [deepcopy(restart_weights!(model.classifier)) for _ in 1:n_candidates]
    histories = []
    for (i,candidate) in enumerate(candidates)
        history = train_classifier!(candidate, cl_tr_x, cl_tr_y, cl_val_x, cl_val_y;
            verb = false, model.parameters.classifier...)
        push!(histories, history)
        val_auc = get(history[:val_auc])[2][end]
        @info "Trained classifier $i/$(n_candidates), val AUC = $(val_auc)"
    end

    # now assign the right fields
    # select the best one
    best_candidate_ind = argmax([get(history[:val_auc])[2][end] for history in histories])
    model.classifier = candidates[best_candidate_ind]
    model.training_output.classifier[:history] = histories[best_candidate_ind]

    @info "Selected classifier $(best_candidate_ind)"

    @info "Classifier training finished. \n#############################\n"
end
_split_data(x::AbstractArray{T,2}, y, vr) where T =
    train_val_test_split(x[:,y.==0], x[:,y.==1],(0, 1-vr, vr))
_split_data(x::AbstractArray{T,4}, y, vr) where T =
    train_val_test_split(x[:,:,:,y.==0], x[:,:,:,y.==1],(0, 1-vr, vr))

"""
	StatsBase.fit!(model::HAD, tr_x::AbstractArray, val_x::AbstractArray)
"""
function StatsBase.fit!(model::HAD, tr_x::AbstractArray, val_x::AbstractArray, val_y::Vector)
    fit_autoencoder!(model, tr_x, val_x)
    fit_detectors!(model, tr_x)
    fit_classifier!(model, val_x, val_y)
end
"""
	StatsBase.predict(model::HAD, x::AbstractArray; L=10)
"""
function StatsBase.predict(model::HAD, x::AbstractArray; L=10)
    cl_x = classifier_inputs(model, x; L=L)
    cl_x = normalize(cl_x, model.training_output.classifier[:normalization_coefficients]...)
    classifier_score(model.classifier, cl_x)
end

### UTILS ###
"""
	ensemble_scores(detectors, encodings)

Scores of the ensemble of detectors.
"""
function ensemble_scores(detectors, x)
    scores = map(_x->StatsBase.predict(_x[1], _x[2]), zip(detectors, x))
end

"""
	all_scores(model::HAD, x; L::Int=10)

Detector + AE rec. probability, L is number of samples in the computation of rec. probability. Note that
rec. probability is only used if 'model.parameters.classifier.use_reconstruction==true' (default).
"""
function all_scores(model::HAD, x; L::Int=10)
    detector_scores = ensemble_scores(
        model.detectors,
        cpu(encode_all(model.autoencoder, x; batchsize=model.parameters.autoencoder[:batchsize], mean=true))
        )
    userec = get(model.parameters.classifier, :use_reconstruction, true)
    if !userec
        return detector_scores
    else
        rec_score = reconstruction_probability(model.autoencoder, x, L, model.parameters.autoencoder[:batchsize])
        return cat(detector_scores, [rec_score], dims=1)
    end
end

"""
	classifier_inputs(model::HAD, x[, y]; L=10)

Transforms x (and y) into inputs for the classifier - combination of all scores.
"""
classifier_inputs(model::HAD, x; L=10) = Array(transpose(cat(all_scores(model, x; L=L)..., dims = 2)))
function classifier_inputs(model::HAD, x, y; L=10)
    cl_x = classifier_inputs(model::HAD, x; L=L)
    cl_y = Flux.onehotbatch(Bool.(y), 0:1);
    cl_x, cl_y
end

"""
	classifier_score(classifier, x)

Extracts the probability that a sample is an anomaly.
"""
classifier_score(c, x) = softmax(c(x))[2,:]
classifier_score(c, x::AbstractArray{T,4}) where T = cpu(softmax(c(gpu(x)))[2,:])
classifier_score(c, x, batchsize::Int) =
    vcat(map(b->cpu(classifier_score(c, b)), Flux.Data.DataLoader(x, batchsize=batchsize))...)

### extension of some methods
reconstruct(model::HAD, x) = reconstruct(model.autoencoder, x)

### some other stuff
function autoencoder_params_outputs_from_fvlae_data(autoencoder_data, datasize)
    autoencoder_parameters = (
        zdim=autoencoder_data[:experiment_args].latent_dim,
        hdim=autoencoder_data[:experiment_args].hdim,
        ks=autoencoder_data[:experiment_args].kernelsizes,
        ncs=autoencoder_data[:experiment_args].channels,
        strd=autoencoder_data[:experiment_args].stride,
        datasize=datasize,
        nepochs=autoencoder_data[:experiment_args].nepochs,
        λ=autoencoder_data[:experiment_args].lambda,
        γ=autoencoder_data[:experiment_args].gamma,
        epochsize=autoencoder_data[:experiment_args].epochsize,
        layer_depth=autoencoder_data[:experiment_args].layer_depth,
        lr=autoencoder_data[:experiment_args].lr,
        pad=autoencoder_data[:experiment_args].pad,
        var= (autoencoder_data[:experiment_args].last_conv) ? :conv : :dense,
        activation=autoencoder_data[:experiment_args].activation,
        xdist=autoencoder_data[:experiment_args].xdist,
        discriminator_nlayers=autoencoder_data[:experiment_args].discriminator_nlayers,
        batchsize=autoencoder_data[:experiment_args].batchsize
        )
    autoencoder_outputs = Dict(
        :latents => autoencoder_data[:latent_representations],
        :history => autoencoder_data[:training_history],
        :reconstructions => autoencoder_data[:reconstructions]
        )

    autoencoder_parameters, autoencoder_outputs
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
function computeCI_TStudent(r::River, level::Float64=.95; ci_period::Period=Week(1))
  TS = getScenarios(r)
  period = getPeriod(r)
  time_steps = length(TS[1])
  ci_low = zeros(Float64, time_steps)
  ci_high = zeros(Float64, time_steps)
  means = zeros(Float64, time_steps)
  matrix = vcat([values(ts)' for ts in TS]...) # Indexes: scenario, time step.

  # Compute the confidence interval for these time steps, one at a time.
  for t in 1:size(matrix, 2)
    data = vec(matrix[:, t])
    means[t] = mean(data)
    period_ci = confint(OneSampleTTest(data), 1 - level)
    ci_low[t] = (period_ci[1] < 0) ? 0 : period_ci[1]
    ci_high[t] = (period_ci[2] < 0) ? 0 : period_ci[2]
  end

  # Copy things as required.
  if period != ci_period # Only possible case: get days in, compute on weeks, must output days.
    # Computations on one week, average of seven days.
    means /= 7
    ci_low /= 7
    ci_high /= 7

    means[1:364] = vec(repmat(means[1:52]', 7, 1))
    means[365] = means[364] # Scenarios are one-year long. Ensure a year is 365 days, and not 364 = 7 * 52 days.

    ci_low[1:364] = vec(repmat(ci_low[1:52]', 7))
    ci_high[1:364] = vec(repmat(ci_high[1:52]', 7))
    ci_low[365, :] = ci_low[364, :]
    ci_high[365, :] = ci_high[364, :]
  end

  # Finally, make new time series out of it.
  ts_means   = TimeArray(timestamp(TS[1]), means,   colnames(TS[1]), meta(TS[1]))
  ts_ci_low  = TimeArray(timestamp(TS[1]), ci_low,  colnames(TS[1]), meta(TS[1]))
  ts_ci_high = TimeArray(timestamp(TS[1]), ci_high, colnames(TS[1]), meta(TS[1]))

  return ts_means, ts_ci_low, ts_ci_high
end


text:
using CSV
using Plots
using DataFrames

csv_file_jud = try
    ARGS[1]
catch
    "benchmark_output_jud.csv"
end

csv_file_pul = try
    ARGS[2]
catch
    "benchmark_output_pul.csv"
end

p_threads = plot(
     xlabel="threads",
     ylabel="time (s)",
     yscale=:log,
     legendtitle="resol",
     legend = :outerright,
)

df_jud = DataFrame(CSV.File(csv_file_jud, types=[Int64, Int64, Float64, Float64]))
df_pul = DataFrame(CSV.File(csv_file_pul, types=[Int64, Int64, Float64, Float64]))

for (resol, color) in zip(unique(df_jud.resol), [:blue, :red, :black])

    df_jud_filtered = df_jud[df_jud.resol .== resol, :]
    df_pul_filtered = df_pul[df_pul.resol .== resol, :]

    sort!(df_jud_filtered, :threads)

    df_jud_filtered[!, :theory] = df_jud_filtered[1, :time] .* df_jud_filtered[1, :threads] ./ df_jud_filtered[!, :threads]

    plot!(
        p_threads,
        df_jud_filtered.threads,
        df_jud_filtered.time,
        seriestype = :scatter,
        color=color,
        label = "UltraDark.jl $resol^3",
    )

    plot!(
        p_threads,
        df_pul_filtered.threads,
        df_pul_filtered.time,
        seriestype = :scatter,
        color=color,
        marker=:x,
        label = "PyUltraLight $resol^3",
    )

    plot!(
        p_threads,
        df_jud_filtered.threads,
        df_jud_filtered.theory,
        color=color,
        label = "theory $resol^3",
    )

end

savefig("threads.pdf")

p_resol = plot(
     xlabel="resol",
     ylabel="time (s)",
     yscale=:log,
     legendtitle="threads",
     legend = :outerright,
)

for threads in [1, 8]

    df_jud_filtered = df_jud[df_jud.threads .== threads, :]

    sort!(df_jud_filtered, :resol)

    df_jud_filtered[!, :theory] = df_jud_filtered[1, :time] ./ df_jud_filtered[1, :resol]^3 .* df_jud_filtered[!, :resol].^3

    plot!(
        p_resol,
        df_jud_filtered.resol,
        df_jud_filtered.time,
        seriestype = :scatter,
        label = "$threads",
    )

    plot!(
        p_resol,
        df_jud_filtered.resol,
        df_jud_filtered.theory,
        label = "theory, $threads",
    )

end

savefig("resol.pdf")

p = plot(p_threads, p_resol)
savefig("bench.pdf")


text:
"""
    Script to generate lineplot graphs associated with:
    - Experiment 3 (propV_randE.jl)
    - Experiment 4 (propV_propE05.jl)
"""

using Pkg
Pkg.activate(".")
using LTMSim
using SimpleHypergraphs
using PyPlot
using Statistics
using Serialization
using LaTeXStrings

project_path = dirname(pathof(LTMSim))

fname = "propV_randE_rev.data" #"propV_propE05.data" #"propV_randE.data"
data_path = joinpath(project_path, "..", "res", "journal", fname)
res_path = joinpath(project_path, "..", "res", "journal", "plot", "propV_randE") #propV_randE propV_propE05

hg_files = readdir(joinpath(project_path, "..", "data", "hgs"))
hg_names = [split(file, ".")[1] for file in hg_files]

#
# DATA
#
data = deserialize(data_path)
hgs = [prune_hypergraph!(hg_load(joinpath(project_path, "..", "data", "hgs", hg_file))) for hg_file in hg_files]

#
# PLOTTING INFO
#
algorithms = [
    "BinarySearch(H)", "Greedy([H]₂)", "Greedy(H)", "SubTSS(H)",
    "BinarySearch(H)-noOpt", "Greedy([H]₂)-noOpt", "Greedy(H)-noOpt", "SubTSS(H)-noOpt"
]

labels_dict = Dict{String, String}(
    "BinarySearch(H)" => "StaticGreedy",
    "BinarySearch(H)-noOpt" => "StaticGreedy-noOpt",
    "Greedy(H)" => "DynamicGreedy",
    "Greedy(H)-noOpt" => "DynamicGreedy-noOpt",
    "Greedy([H]₂)" => L"DynamicGreedy_{[H]_2}",
    "Greedy([H]₂)-noOpt" => L"DynamicGreedy_{[H]_2}-noOpt",
    "SubTSS(H)" => "SubTSS",
    "SubTSS(H)-noOpt" => "SubTSS-noOpt"
)

#
# LINEPLOT
#
nvalues = range(0.2, stop=0.8, step=0.1)

colorz=["#2C7BB6", "#D7191C", "#FF8900", "#33CC33", "#cc0099", "#4d4dff", "#008080", "#2C7BB6"]

for (index, hg_name) in enumerate(hg_names)
    println(hg_name)

    clf()
    plt.figure(figsize=(7,5))

    val = -0.4
    c = 1

    labels = Array{String, 1}()

    for algo in algorithms
        y = Array{Float64, 1}()

        for v=1:length(collect((nvalues)))
            yval = mean(data[algo][hg_name][v] / nhv(hgs[index]))
            push!(y, yval)
        end

        push!(
            labels,
            labels_dict[algo]
        )

        b = plt.plot(collect(nvalues), y, color=colorz[c])

        c+=1
    end

    plt.legend(labels, ncol=2)#, fontsize="x-large")

    plt.xticks(nvalues, fontsize="x-large") #range(0, length(ticks) * 2, step=2),
    plt.yticks(fontsize="x-large")

    #plt.xlim(-2, length(ticks)*2)
    plt.ylim(0, 1.05)

    ylabel("Seed set size / n", fontstyle = "italic", fontsize="xx-large", labelpad=10) #fontweight="semibold",
    xlabel("Thresholds", fontstyle = "italic", fontsize="xx-large", labelpad=10) #, fontweight="semibold"

    title("propV_randE_$hg_name", fontstyle = "italic", fontsize="xx-large")

    plt.tight_layout()
    gcf()

    PyPlot.savefig(joinpath(res_path, "propV_randE_$hg_name.png"))

    plt.close()
end

h = hgs[5]

w3c_v_dist = [length(gethyperedges(h, v)) for v in 1:nhv(h)]

sum(w3c_v_dist[w3c_v_dist .== 1]) * 100 / length(w3c_v_dist)


text:
using DeconvOptim, TestImages, Images, FFTW, Noise, ImageView
using Plots

function norm_fft(x)
    x = fftshift(abs.(fft(x)))
    n = div(size(x)[1], 2)
    x ./= maximum(x)
end

function ideal_freq()
    img = zeros((512, 512))
    for i = 1:1
        img[rand((1:512)), rand((1:512))] = 1
    end
    img = fftshift(img)
    N_phot = 376
    img ./= maximum(img)
    img = convert(Array{Float32}, img .* N_phot)
    dist = [sqrt((-1 + i - size(img)[1] / 2)^2 + (-1 + j - size(img)[2] / 2)^2)
                for i = 1:size(img)[1],  j = 1:size(img)[2]]
    psf = ifftshift(exp.(-dist .^2 ./ 5.0 .^2))
    psf ./= sum(psf)
    psf = convert(Array{Float32}, psf)

    #img_b = center_extract(conv_psf(center_set!(copy(z1), img), ifftshift(center_set!(z, fftshift(psf))), [1, 2]), size(img))
    img_b = conv_psf(img, psf, [1, 2])
    img_n = poisson(img_b, N_phot)

    reg = DeconvOptim.TV(num_dims=2, sum_dims=[1, 2])
    @time res, o = deconvolution(img_n, psf, iterations=10, λ=0.001f0,
            loss=Poisson(), regularizer=reg, padding=0.00, plan_fft=true)
    img_ft = norm_fft(img)[:, 257]
    img_n_ft = norm_fft(img_n)[:, 257]
    res_ft = norm_fft(res)[:, 257]
    freq = fftshift(fftfreq(512, 1))
    mtf = norm_fft(psf)[:, 257]
    plot(freq, mtf, xlabel="Frequency in 1/pixel", ylabel = "Normalized intensity in Frequency space", label="OTF", dpi=300)
    plot!(freq, img_ft,label="Image with Constant Frequency content")
    plot!(freq, img_n_ft, label="Blurry image with noise", linestyle = :dot)
    plot!(freq, res_ft, label="Deconvolved image")

    savefig("src/assets/ideal_frequencies.png")
end

ideal_freq()


text:
mutable struct NothingXMLElement <: MyXMLElement
    el::Nothing
    NothingXMLElement() = new(nothing)
end

function make_xml(::NothingXMLElement)
    # do nothing
end


text:
# implementation of the GPUCompiler interfaces for generating WASM code

## target

export WASMCompilerTarget

Base.@kwdef struct WASMCompilerTarget <: AbstractCompilerTarget
end

llvm_triple(::WASMCompilerTarget) = "wasm32-unknown-unknown"

function llvm_machine(target::WASMCompilerTarget)
  triple = llvm_triple(target)

  t = Target(triple=triple)

  cpu = ""
  feat = ""
  tm = TargetMachine(t, triple, cpu, feat)
  asm_verbosity!(tm, true)

  return tm
end

function process_entry!(job::CompilerJob{WASMCompilerTarget}, mod::LLVM.Module, entry::LLVM.Function)
  push!(function_attributes(entry), StringAttribute("wasm-export-name", String(chop(LLVM.name(entry), tail=4)), context(mod)))
  invoke(process_entry!, Tuple{CompilerJob, LLVM.Module, LLVM.Function}, job, mod, entry)
end

## job

runtime_slug(job::CompilerJob{WASMCompilerTarget}) = "wasm"


text:
# symmetric, skew-symmetric and symmetric checks
"""
    symmetric(::SecondOrderTensor)
    symmetric(::FourthOrderTensor)

Computes the (minor) symmetric part of a second or fourth order tensor.
Return a `SymmetricTensor`.

# Examples
```jldoctest
julia> A = rand(Tensor{2,2})
2×2 Tensor{2, 2, Float64, 4}:
 0.590845  0.566237
 0.766797  0.460085

julia> symmetric(A)
2×2 SymmetricTensor{2, 2, Float64, 3}:
 0.590845  0.666517
 0.666517  0.460085
```
"""
@inline symmetric(S1::SymmetricTensors) = S1

@inline function symmetric(S::Tensor{2, dim}) where {dim}
    SymmetricTensor{2, dim}(@inline function(i, j) @inbounds i == j ? S[i,j] : (S[i,j] + S[j,i]) / 2 end)
end

"""
    minorsymmetric(::FourthOrderTensor)

Compute the minor symmetric part of a fourth order tensor, return a `SymmetricTensor{4}`.
"""
@inline function minorsymmetric(S::Tensor{4, dim}) where {dim}
    SymmetricTensor{4, dim}(
        @inline function(i, j, k, l)
            @inbounds if i == j && k == l
                return S[i,j,k,l]
            else
                return (S[i,j,k,l] + S[j,i,k,l] + S[i,j,k,l] + S[i,j,l,k]) / 4
            end
        end
    )
end

@inline minorsymmetric(S::SymmetricTensors) = S

@inline symmetric(S::Tensor{4}) = minorsymmetric(S)

"""
    majorsymmetric(::FourthOrderTensor)

Compute the major symmetric part of a fourth order tensor.
"""
@inline function majorsymmetric(S::T) where {T <: FourthOrderTensor}
    get_base(T)(
        @inline function(i, j, k, l)
            @inbounds if i == j == k == l || i == k && j == l
                return S[i,j,k,l]
            else
                return (S[i,j,k,l] + S[k,l,i,j]) / 2
            end
        end
    )
end

# Unsafe versions that just truncates
@inline unsafe_symmetric(S::SymmetricTensors) = S
@inline function unsafe_symmetric(S::Tensor{2, dim}) where {dim}
    SymmetricTensor{2, dim}(@inline function(i, j) @inbounds S[i,j] end)
end
@inline function unsafe_symmetric(S::Tensor{4, dim}) where {dim}
    SymmetricTensor{4, dim}(@inline function(i, j, k, l) @inbounds return S[i,j,k,l] end)
end

"""
    skew(::SecondOrderTensor)

Computes the skew-symmetric (anti-symmetric) part of a second order tensor, returns a `Tensor{2}`.
"""
@inline skew(S1::Tensor{2}) = (S1 - S1') / 2
@inline skew(S1::SymmetricTensor{2,dim,T}) where {dim, T} = zero(Tensor{2,dim,T})

# Symmetry checks
@inline LinearAlgebra.issymmetric(t::Tensor{2, 1}) = true
@inline LinearAlgebra.issymmetric(t::Tensor{2, 2}) = @inbounds t[1,2] == t[2,1]

@inline function LinearAlgebra.issymmetric(t::Tensor{2, 3})
    return @inbounds t[1,2] == t[2,1] && t[1,3] == t[3,1] && t[2,3] == t[3,2]
end

function isminorsymmetric(t::Tensor{4, dim}) where {dim}
    @inbounds for l in 1:dim, k in l:dim, j in 1:dim, i in j:dim
        if t[i,j,k,l] != t[j,i,k,l] || t[i,j,k,l] != t[i,j,l,k]
            return false
        end
    end
    return true
end

isminorsymmetric(::SymmetricTensor{4}) = true

function ismajorsymmetric(t::FourthOrderTensor{dim}) where {dim}
    @inbounds for l in 1:dim, k in l:dim, j in 1:dim, i in j:dim
        if t[i,j,k,l] != t[k,l,i,j]
            return false
        end
    end
    return true
end

LinearAlgebra.issymmetric(t::Tensor{4}) = isminorsymmetric(t)

LinearAlgebra.issymmetric(::SymmetricTensors) = true


text:
r=1.0;s=1.0;t=1.0;

Relaxation_success = true

# psi(x::Any,r::Float64,s::Float64,t::Float64)
@test norm(MPCCSolver.psi(s,r,s,t)-s) <= eps(Float64)*10
# x = linspace(-0.5,10,1000); y = MPCCSolver.psi(x,r,s,t)
# PyPlot.plot(x, y, color="red", linewidth=2.0, linestyle="-")
# x = linspace(-0.5,10,1000); y = MPCCSolver.psi(x,r,s,t)
# PyPlot.plot(y, x, color="blue", linewidth=2.0, linestyle="-")

#invpsi(x::Any,r::Float64,s::Float64,t::Float64)
@test norm(MPCCSolver.invpsi(s,r,s,t)-s) <= eps(Float64)*10

# x = linspace(-0.5,10,1000); y = MPCCSolver.psi(x,r,s,t)
# PyPlot.plot(x, y, color="red", linewidth=2.0, linestyle="-")
# x = linspace(-0.5,10,1000); y = MPCCSolver.psi(x,r,s,t)
# PyPlot.plot(y, x, color="blue", linewidth=2.0, linestyle="-")

#dpsi(x::Any,r::Float64,s::Float64,t::Float64)
#phi(yg::Any,yh::Any,r::Float64,s::Float64,t::Float64)
#phi(x::Any,nb_comp::Int64,r::Float64,s::Float64,t::Float64)
#dphi(yg::Any,yh::Any,r::Float64,s::Float64,t::Float64)
#dphi(x::Any,nb_comp::Int64,r::Float64,s::Float64,t::Float64)

#alpha_max(x::Float64,dx::Float64,y::Float64,dy::Float64,r::Float64,s::Float64,t::Float64)
#alpha_max : résoud l'équation en alpha où psi(x,r,s,t)=s+t*theta(x,r)
r=1.0;s=1.0;t=0.0 #résoud l'éq.

x=1.0;y=2.0;dx=1.0;dy=-1.0;
#doit rendre (0,1)
@test isequal(MPCCSolver.alpha_max(x,dx,y,dy,r,s,t),[0.0,1.0])

x=0.0;y=2.0;dx=1.0;dy=-2.0;
#solution : (1,0.5)
@test isequal(MPCCSolver.alpha_max(x,dx,y,dy,r,s,t),[1.0,0.5])

r=1.0;s=1.0;t=1.0
#on tape la contrainte :
x=1.0;y=2.0;dx=1.0;dy=0.0;

@test minimum(MPCCSolver.alpha_max(x,dx,y,dy,r,s,t))>=0.0
@test maximum(MPCCSolver.alpha_max(x,dx,y,dy,r,s,t))==Inf

x=2.0;y=1.0;dx=-1.0;dy=2.0;

@test minimum(MPCCSolver.alpha_max(x,dx,y,dy,r,s,t)) >= 0.0
@test maximum(MPCCSolver.alpha_max(x,dx,y,dy,r,s,t)) <  Inf

#on est sous (s,s) et donc on ne tape pas la contrainte virtuelle
x=0.0;y=0.0;dx=0.0;dy=-2.0;

@test isequal(maximum(MPCCSolver.alpha_max(x,dx,y,dy,r,s,t)),Inf)

#Conclusion:
@test Relaxation_success


text:
module DobotMagician

using LibSerialPort: LibSerialPort
using LibSerialPort.Lib:
    SP_TRANSPORT_USB,
    SP_MODE_READ_WRITE,
    SP_PARITY_NONE,
    SP_FLOWCONTROL_NONE,
    sp_list_ports,
    sp_get_port_transport,
    sp_get_port_usb_vid_pid,
    sp_get_port_name,
    sp_free_port_list,
    sp_get_port_by_name,
    sp_open,
    sp_set_baudrate,
    sp_set_bits,
    sp_set_stopbits,
    sp_set_parity,
    sp_set_flowcontrol,
    sp_close,
    sp_free_port,
    sp_blocking_write,
    sp_blocking_read,
    sp_drain
using DocStringExtensions: SIGNATURES

const Port = LibSerialPort.Lib.Port

mutable struct Magician
    port::Port
    timeout::Int32
    isopen::Bool
    function Magician(port::Port, timeout::Integer)
        dobot = new(port, timeout, false)
        finalizer(destroy!, dobot)
        connect(dobot)
        return dobot
    end
end

function Magician(port::AbstractString; timeout=1000)
    return Magician(sp_get_port_by_name(port), timeout)
end

function Magician(; kwargs...)
    port = find_magician()
    if port === nothing
        throw(ErrorException("Dobot Magician not found"))
    else
        return Magician(port; kwargs...)
    end
end

function destroy!(dobot::Magician)
    disconnect(dobot)
    return sp_free_port(dobot.port)
end

"""
    $(SIGNATURES)

Return the serial port used by the attached Dobot Magician. Returns nothing if no
appropriate interface is found.

# Keywords
- `nports_guess=64`: maximum number of serial ports to inspect to find the Dobot interface.

# Returns
- `Union{Port, Nothing}`
"""
function find_magician(; nports_guess=64)
    dobot_port = nothing
    ports = sp_list_ports()
    for port in unsafe_wrap(Array, ports, nports_guess; own=false)
        if port == C_NULL
            break
        elseif sp_get_port_transport(port) == SP_TRANSPORT_USB
            if sp_get_port_usb_vid_pid(port) == (0x10c4, 0xea60)
                dobot_port = sp_get_port_name(port)
                break
            end
        end
    end
    sp_free_port_list(ports)
    return dobot_port
end

"""
    $(SIGNATURES)

Open the serial port to connect to the Dobot Magician.

# Returns
- `Magician`: a connection to the Dobot Magician.
"""
function connect(dobot::Magician)
    if dobot.isopen
        disconnect(dobot)
    end
    sp_open(dobot.port, SP_MODE_READ_WRITE)
    dobot.isopen = true
    # Values taken from the Dobot Magician Communication Protocol document
    sp_set_baudrate(dobot.port, 115200)
    sp_set_bits(dobot.port, 8)
    sp_set_stopbits(dobot.port, 1)
    sp_set_parity(dobot.port, SP_PARITY_NONE)
    sp_set_flowcontrol(dobot.port, SP_FLOWCONTROL_NONE)
    return dobot
end

"""
    $(SIGNATURES)

Close the serial port to the Dobot Magician.
"""
function disconnect(dobot::Magician)
    if dobot.isopen
        sp_close(dobot.port)
        dobot.isopen = false
    end
    return dobot
end

"""
    Command{P}
"""
struct Command{P}
    id::UInt8
    rw::Bool
    allow_queue::Bool
    description::String
    payload::Vector{Union{Tuple{Symbol,DataType},Tuple{Symbol,UnionAll}}}
end

function Command(id, rw, allow_queue, payload, description)
    return Command{Tuple{values(payload)...}}(
        id, rw, allow_queue, description, [(name, type) for (name, type) in pairs(payload)]
    )
end

function (cmd::Command)(magician::Magician, args...; kwargs...)
    return execute_command(magician, cmd, args...; kwargs...)
end

function Base.show(io::IO, ::MIME"text/plain", cmd::Command)
    return print(
        io,
        """
            Dobot Magician command
            ID:          $(cmd.id)
            R/W:         $(cmd.rw)
            Allow queue: $(cmd.allow_queue)
            Description: $(cmd.description)
            Payload:
                $(join([string(name)*"::"*string(type) for (name, type) in cmd.payload], "\n    "))""",
    )
end

const TIMEOUT_MSG = "Serial port errored or timed out waiting for response"

function execute_command(magician::Magician, cmd::Command, payload=(); queue=false)
    # Construct and send the command package
    package = construct_command(cmd, payload, queue)
    sp_blocking_write(magician.port, package, magician.timeout)
    sp_drain(magician.port)  # ensure that the command package has been sent
    # Read the header + ID + Ctrl
    nbytes, header = sp_blocking_read(magician.port, 5, magician.timeout)
    if nbytes < 5
        throw(ErrorException(TIMEOUT_MSG))
    end
    if (header[1] != 0xaa) || (header[2] != 0xaa)
        throw(ErrorException("Invalid header"))
    end
    sz = header[3]
    id = header[4]
    ctrl = header[5]
    # Check the ID
    if id != cmd.id
        throw(ErrorException("Incorrect ID returned"))
    end
    # Check the Ctrl
    if ((ctrl & 0x1) > 0) != cmd.rw
        throw(ErrorException("Incorrect rw returned"))
    end
    if ((ctrl & 0x2) > 0) != queue
        throw(ErrorException("Incorrect isQueued returned"))
    end
    # Read the rest of the payload
    nbytes, payload = sp_blocking_read(magician.port, sz - 2, magician.timeout)
    if nbytes < (sz - 2)
        throw(ErrorException(TIMEOUT_MSG))
    end
    # Read the checksum
    nbytes, checksum = sp_blocking_read(magician.port, 1, magician.timeout)
    if nbytes < 1
        throw(ErrorException(TIMEOUT_MSG))
    end
    # Check the checksum
    checksum_actual = -id - ctrl
    for x in payload
        checksum_actual -= x
    end
    if checksum_actual != checksum[1]
        throw(ErrorException("Invalid checksum returned"))
    end
    # Unpack the result
    if queue
        return read(IOBuffer(payload), UInt64)
    elseif !cmd.rw
        return unpack_payload(cmd, payload)
    else
        return nothing
    end
end

write_var_as_type(io::IO, T::Type, var) = (write(io, convert(T, var)); nothing)

function write_var_as_type(io::IO, T::Type{String}, var)
    str = string(var)
    if !isascii(str)
        throw(ArgumentError("Strings must only use ASCII characters"))
    end
    if !isempty(str) && str[end] != '\0'
        str *= '\0'
    end
    write(io, str)
    return nothing
end

_write_var_as_type(T::Type, ::Any, varsym) = :(write_var_as_type(io, $T, $varsym))

function _write_var_as_type(T::Type{<:Tuple}, var, varsym)
    if T isa UnionAll
        # This occurs where there can be an unspecified number of repeated units
        # Note the body of the UnionAll of an NTuple is Tuple{Vararg{innerT, N}}
        body = T.body
        if !(body isa UnionAll) && (body.parameters[1] <: Vararg)
            innerT = body.parameters[1].parameters[1]
            N = var <: Tuple ? length(var.parameters) : length(var)  # allows for StaticArrays
            return _write_var_as_type(NTuple{N,innerT}, var, varsym)  # call again with bound N
        else
            throw(ArgumentError("Got a UnionAll that is not an NTuple with unbound N (only)"))
        end
    else
        body = quote end
        for (i, TT) in enumerate(T.parameters)
            push!(
                body.args,
                _write_var_as_type(
                    TT, var <: Tuple ? var.parameters[i] : eltype(var), :($(varsym)[$i])
                ),
            )
        end
        push!(body.args, :(nothing))
        return body
    end
end

@generated function write_var_as_type(io::IO, T::Type{<:Tuple}, var)
    # Use a generated function to unroll the loop over all the components of the Tuple
    return _write_var_as_type(T.parameters[1], var, :var)
end

function construct_command(cmd::Command{P}, payload, queue::Bool) where {P}
    # Uses dispatch to enforce correct payload types
    (!cmd.allow_queue && queue) && throw(ArgumentError("Command cannot be queued"))
    # Construct the message
    io = IOBuffer()
    write(io, 0xaaaa)
    write(io, 0x0) # write len later
    write(io, cmd.id)
    write(io, UInt8(cmd.rw) | (UInt8(queue) << 1))
    if cmd.rw
        write_var_as_type(io, P, payload)
    end
    sz = UInt8(position(io) - 3)
    seek(io, 2)
    write(io, sz)
    checksum = 0x0
    for i in Base.OneTo(sz)
        checksum -= read(io, UInt8)
    end
    write(io, checksum)
    return take!(io)
end

read_var_as_type(io::IO, T::Type) = read(io, T)
read_var_as_type(io::IO, T::Type{Vector{UInt8}}) = read(io)
read_var_as_type(io::IO, T::Type{String}) = rstrip(read(io, String), ('\0',))

_read_var_as_type(T::Type) = :(read_var_as_type(io, $T))

function _read_var_as_type(T::Type{<:Tuple})
    if T isa UnionAll
        throw(ArgumentError("Cannot read varargs"))
    else
        body = :(())
        for TT in T.parameters
            push!(body.args, _read_var_as_type(TT))
        end
        return body
    end
end

@generated function read_var_as_type(io::IO, T::Type{<:Tuple})
    return _read_var_as_type(T.parameters[1])
end

function unpack_payload(cmd::Command{P}, payload) where {P}
    io = IOBuffer(payload)
    return read_var_as_type(io, P)
end

#! format: off

# Commands - Device information
const set_device_sn = Command(0, true, false, (device_sn=String,), "Set device serial number")
const get_device_sn = Command(0, false, false, (device_sn=String,), "Get device serial number")
const set_device_name = Command(1, true, false, (device_name=String,), "Set device name")
const get_device_name = Command(1, false, false, (device_name=String,), "Get device name")
const get_device_version = Command(2, false, false, (major_version=UInt8, minor_version=UInt8, revision=UInt8), "Get device version")
const set_device_with_l = Command(3, true, false, (is_with_l=UInt8, version=UInt8), "Set sliding rail enable status")
const get_device_with_l = Command(3, false, false, (is_with_l=UInt8,), "Get sliding rail enable status")
const get_device_time = Command(4, false, false, (g_systick=UInt32,), "Get device time")
const get_device_info = Command(6, false, false, (runtime=UInt64, power_on=UInt32, power_off=UInt32), "Get device info")

# Commands - Real-time pose
const get_pose = Command(10, false, false, (xyzr=NTuple{4,Float32}, J=NTuple{4,Float32}), "Get the real-time pose")
const reset_pose = Command(11, true, false, (manual=UInt8, rear_arm_angle=Float32, front_arm_angle=Float32), "Reset the real-time pose")
const get_kinematics = Command(12, false, false, (velocity=Float32, acceleration=Float32), "Get the kinematics")
const get_pose_l = Command(13, false, false, (pose_l=Float32,), "Get the real-time pose of the sliding rail")

# Commands - Alarm
const get_alarms_state = Command(20, false, false, (alarms_state=Vector{UInt8},), "Get the alarms state")
const clear_all_alarms_state = Command(20, true, false, (), "Clear the alarms state")

# Commands - Homing
const set_home_params = Command(30, true, true, (x=Float32, y=Float32, z=Float32, r=Float32), "Set the homing position")
const get_home_params = Command(30, false, false, (x=Float32, y=Float32, z=Float32, r=Float32), "Get the homing position")
const set_home_cmd = Command(31, true, true, (reserved=UInt32,), "Execute the homing function")
const set_auto_leveling = Command(32, true, true, (is_autoleveling=UInt8, accuracy=Float32), "Set automatic leveling")
const get_auto_leveling = Command(32, false, false, (result=Float32,), "Get automatic leveling result")

# Commands - Handhold teaching
const set_hht_trig_mode = Command(40, true, false, (hht_trig_mode=UInt8,), "Set hand-hold teaching mode")
const get_hht_trig_mode = Command(40, false, false, (hht_trig_mode=UInt8,), "Get hand-hold teaching mode")
const set_hht_trig_output_enabled = Command(41, true, false, (is_enabled=UInt8,), "Set status of hand-hold teaching mode")
const get_hht_trig_output_enabled = Command(41, false, false, (is_enabled=UInt8,), "Get status of hand-hold teaching mode")
const get_hht_trig_output = Command(42, false, false, (is_triggered=UInt8,), "Get the status of the hand-hold trigger")

# Commands - Arm orientation
const set_arm_orientation = Command(50, true, true, (orientation=UInt8,), "Set the handedness of the arm (0=left, 1=right)")
const get_arm_orientation = Command(50, false, false, (orientation=UInt8,), "Get the handedness of the arm (0=left, 1=right)")

# Commands - End effector
const set_end_effector_params = Command(60, true, true, (x_bias=Float32, y_bias=Float32, z_bias=Float32), "Set the offset of the end effector")
const get_end_effector_params = Command(60, false, false, (x_bias=Float32, y_bias=Float32, z_bias=Float32), "Get the offset of the end effector")
const set_end_effector_laser = Command(61, true, true, (enable_ctrl=UInt8, on=UInt8), "Set the status of the laser")
const get_end_effector_laser = Command(61, false, false, (is_ctrl_enabled=UInt8, is_on=UInt8), "Get the status of the laser")
const set_end_effector_suction_cup = Command(62, true, true, (is_ctrl_enabled=UInt8, is_sucked=UInt8), "Set the status of the suction cup")
const get_end_effector_suction_cup = Command(62, false, false, (is_ctrl_enabled=UInt8, is_sucked=UInt8), "Get the status of the suction cup")
const set_end_effector_gripper = Command(63, true, true, (is_ctrl_enabled=UInt8, is_gripped=UInt8), "Set the status of the gripper")
const get_end_effector_gripper = Command(63, false, false, (is_ctrl_enabled=UInt8, is_gripped=UInt8), "Set the status of the gripper")

# Commands - Jog
const set_jog_joint_params = Command(70, true, true, (velocity=NTuple{4,Float32}, acceleration=NTuple{4,Float32}), "Set the velocity and acceleration of joints in jogging mode")
const get_jog_joint_params = Command(70, false, false, (velocity=NTuple{4,Float32}, acceleration=NTuple{4,Float32}), "Get the velocity and acceleration of joints in jogging mode")
const set_jog_coordinate_params = Command(71, true, true, (velocity=NTuple{4,Float32}, acceleration=NTuple{4,Float32}), "Set the velocity and acceleration in Cartesian coordinates in jogging mode")
const get_jog_coordinate_params = Command(71, false, false, (velocity=NTuple{4,Float32}, acceleration=NTuple{4,Float32}), "Get the velocity and acceleration in Cartesian coordinates in jogging mode")
const set_jog_common_params = Command(72, true, true, (velocity_ratio=Float32, acceleration_ratio=Float32), "Set the velocity and acceleration ratios of the sliding rail")
const get_jog_common_params = Command(72, false, false, (velocity_ratio=Float32, acceleration_ratio=Float32), "Get the velocity and acceleration ratios of the sliding rail")
const set_jog_cmd = Command(73, true, true, (is_joint=UInt8, cmd=UInt8), "Execute the jog command in Cartesian coordinates or joints (cmd=0 stop, cmd=1 or 2 for positive or negative jogging in x/J1, cmd=3 or 4 for y/J2, etc)")
const set_jog_l_params = Command(74, true, true, (velocity=Float32, acceleration=Float32), "Set the velocity and acceleration of the sliding rail in jog mode") # NOTE: there is a discrepancy in the manual regarding the data types and overall length of the message

# Commands - PTP
const set_ptp_joint_params = Command(80, true, true, (velocity=NTuple{4,Float32}, acceleration=NTuple{4,Float32}), "Set the velocity and acceleration of the joint coordinate axes in PTP mode")
const get_ptp_joint_params = Command(80, false, false, (velocity=NTuple{4,Float32}, acceleration=NTuple{4,Float32}), "Get the velocity and acceleration of the joint coordinate axes in PTP mode")
const set_ptp_coordinate_params = Command(81, true, true, (xyz_velocity=Float32, r_velocity=Float32, xyz_acceleration=Float32, r_acceleration=Float32), "Set the velocity and acceleration of the Cartesian coordinate axes in PTP mode")
const get_ptp_coordinate_params = Command(81, false, false, (xyz_velocity=Float32, r_velocity=Float32, xyz_acceleration=Float32, r_acceleration=Float32), "Get the velocity and acceleration of the Cartesian coordinate axes in PTP mode")
const set_ptp_jump_params = Command(82, true, true, (jump_height=Float32, z_limit=Float32), "Set the lifting height and maximum lifting height in JUMP mode")
const get_ptp_jump_params = Command(82, false, false, (jump_height=Float32, z_limit=Float32), "Get the lifting height and maximum lifting height in JUMP mode")
const set_ptp_common_params = Command(83, true, true, (velocity_ratio=Float32, acceleration_ratio=Float32), "Set the velocity and acceleration ratio in PTP mode")
const get_ptp_common_params = Command(83, false, false, (velocity_ratio=Float32, acceleration_ratio=Float32), "Set the velocity and acceleration ratio in PTP mode")
const set_ptp_cmd = Command(84, true, true, (ptp_mode=UInt8, position=NTuple{4,Float32}), "Execute a PTP command (ptp_mode: 0=JUMP_XYZ, 1=MOVJ_XYZ, 2=MOVL_XYZ, 3=JUMP_ANGLE, 4=MOVJ_ANGLE, 5=MOVL_ANGLE, 6=MOVJ_INC, 7=MOVL_INC, 8=MOVJ_XYZ_INC, 9=JUMP_MOVL_XYZ)")
const set_ptp_l_params = Command(85, true, true, (velocity=Float32, acceleration=Float32), "Set the velocity and acceleration of the sliding rail in PTP mode")
const set_ptp_with_l_cmd = Command(86, true, true, (ptp_mode=UInt8, position=NTuple{5,Float32}), "Execute a PTP command with the sliding rail (ptp_mode: 0=JUMP_XYZ, 1=MOVJ_XYZ, 2=MOVL_XYZ, 3=JUMP_ANGLE, 4=MOVJ_ANGLE, 5=MOVL_ANGLE, 6=MOVJ_INC, 7=MOVL_INC, 8=MOVJ_XYZ_INC, 9=JUMP_MOVL_XYZ)")
const set_ptp_jump2_params = Command(87, true, true, (start_jump_height=Float32, end_jump_height=Float32, z_limit=Float32), "Set the extended parameters in JUMP mode")
const get_ptp_jump2_params = Command(87, false, false, (start_jump_height=Float32, end_jump_height=Float32, z_limit=Float32), "Get the extended parameters in JUMP mode")
const set_ptp_po_cmd = Command(88, true, true, (ptp_mode=UInt8, position=NTuple{4,Float32}, po=NTuple{N,Tuple{UInt8,UInt16,UInt8}} where N), "Execute a PTP command while manipulating the I/O at a certain point (po=(ratio, address, level))")
const set_ptp_po_with_l_cmd = Command(89, true, true, (ptp_mode=UInt8, position=NTuple{5,Float32}, po=NTuple{N,Tuple{UInt8,UInt16,UInt8}} where N), "Execute a PTP command with the sliding rail while manipulating the I/O at a certain point (po=(ratio, address, level))")

# Commands - CP
const set_cp_params = Command(90, true, true, (plan_acc=Float32, junction_vel=Float32, acc_period=Float32, real_time_track=UInt8), "Set the parameters of continuous trajectory; acc_period is the maximum actual acceleration in non-real time mode and the interpolation period in real-time mode")
const get_cp_params = Command(90, false, false, (plan_acc=Float32, junction_vel=Float32, acc_period=Float32, real_time_track=UInt8), "Get the parameters of continuous trajectory; acc_period is the maximum actual acceleration in non-real time mode and the interpolation period in real-time mode")
const set_cp_cmd = Command(91, true, true, (cp_mode=UInt8, xyz=NTuple{3,Float32}, velocity=Float32), "Execute the CP command (cp_mode=0 for relative motion or cp_mode=1 for absolute motion; velocity is reserved/ignored?)")
const set_cp_le_cmd = Command(92, true, true, (cp_mode=UInt8, xyz=NTuple{3,Float32}, power=Float32), "Execute the CP command with laser engraving (cp_mode=0 for relative motion or cp_mode=1 for absolute motion)")

# Commands - Wait
const set_wait_cmd = Command(110, true, true, (timeout=UInt32,), "Wait the specified number of milliseconds")

# Commands - Queued execution control commands
const set_queued_cmd_start_exec = Command(240, true, false, (), "Start the command queue")
const set_queued_cmd_stop_exec = Command(241, true, false, (), "Stop the command queue")
const set_queued_cmd_force_stop_exec = Command(242, true, false, (), "Forcibly stop the command queue (interrupt executing commands)")
const set_queued_cmd_start_download = Command(243, true, false, (total_loop=UInt32, line_per_loop=UInt32), "Start downloading")
const set_queued_cmd_stop_download = Command(244, true, false, (), "Stop downloading")
const set_queued_cmd_clear = Command(245, true, false, (), "Clear the command queue")
const get_queued_cmd_current_index = Command(246, false, false, (queued_cmd_current_index=UInt64,), "Get the command index")
const get_queued_cmd_motion_finish = Command(248, false, false, (isfinish=UInt8,), "Get the current motion status")

#! format: on

# Alarms

"""
    $(SIGNATURES)

Return the index of the active alarm. Returns `nothing` if no alarm is active.
"""
function active_alarm(alarms)
    for (i, alarm) in enumerate(alarms)
        for j in 1:8
            if (alarm & 0x1) == 0x1
                return (i - 1) * 8 + j - 1
            end
            alarm = alarm >> 1
        end
    end
    return nothing
end

"""
    $(SIGNATURES)

Return the description of the specified alarm index.
"""
alarm_description(alarm::Integer) = ALARMS[alarm]
alarm_description(alarm::Vector{UInt8}) = alarm_description(active_alarm(alarm))
alarm_description(dobot::Magician) = alarm_description(get_alarms_state(dobot)[1])
alarm_description(::Nothing) = nothing

const ALARMS = Dict{Int,String}(
    0x00 => "The alarm will be triggered after resetting system",
    0x01 => "Receive undefined instruction",
    0x02 => "File system error",
    0x03 =>
        "There is a failed communication between MCU and FPGA when system is initializing",
    0x04 => "Get an error value of angle sensor",
    # Plan error
    0x10 => "The target point is in abnormal position in Cartesian coordinate system",
    0x11 => "The target point is out of the workspace, that causes inverse resolve alarm",
    0x12 => "The inverse resolve of target is out of the limitation",
    0x13 => "There are some repetitive points in ARC or JUMP_MOVL mode",
    0x15 => "Set a wrong parameter in JUMP mode",
    # Move error
    0x20 =>
        "Motion trajectory is in singularity position that causes inverse resolve alarm in Cartesian coordinate system",
    0x21 =>
        "Motion is out of workspace that causes inverse resolve alarm when Dobot Magician moves",
    0x22 => "The motion is out of the limitation when Dobot Magican moves",
    # Over speed error
    0x30 => "Joint 1 is overspeed",
    0x31 => "Joint 2 is overspeed",
    0x32 => "Joint 3 is overspeed",
    0x33 => "Joint 4 is overspeed",
    # Limit error
    0x40 => "Joint 1 moves to the positive limitation area",
    0x41 => "Joint 1 moves to the negative limitation area",
    0x42 => "Joint 2 moves to the positive limitation area",
    0x43 => "Joint 2 moves to the negative limitation area",
    0x44 => "Joint 3 moves to the positive limitation area",
    0x45 => "Joint 3 moves to the negative limitation area",
    0x46 => "Joint 4 moves to the positive limitation area",
    0x47 => "Joint 4 moves to the negative limitation area",
    0x48 => "Parallelogram is stretched to the positive limitation area",
    0x49 => "Parallelogram is stretched to the negative limitation area",
    # Lose Step error
    0x50 => "Joint 1 loses step",
    0x51 => "Joint 2 loses step",
    0x52 => "Joint 3 loses step",
    0x53 => "Joint 4 loses step",
)

include("Direct.jl")
include("Simple.jl")

end


text:
using TemporalGPs: are_harmonised

function test_harmonise(a, b; recurse=true)
    h = harmonise(a, b)
    @test h isa Tuple
    @test length(h) == 2
    @test are_harmonised(h[1], h[2])

    recurse && test_harmonise(b, a; recurse=false)
    h′ = harmonise(b, a)
    @test h isa Tuple
    @test length(h) == 2
    @test are_harmonised(h′[1], h′[2])
    @test are_harmonised(h[1], h′[1])
    @test are_harmonised(h[1], h′[2])
end

@testset "harmonise" begin
    test_harmonise(5.0, 4.0)

    @testset "AbstractZero" begin
        test_harmonise(5.0, ZeroTangent())
        test_harmonise(ZeroTangent(), randn(10))
        test_harmonise(ZeroTangent(), ZeroTangent())
    end

    @testset "Array" begin
        test_harmonise(randn(5), randn(5))
        test_harmonise(
            [(randn(), randn()) for _ in 1:10],
            [Tangent{Any}(randn(), rand()) for _ in 1:10],
        )
    end

    @testset "Tuple / Tangent{Tuple}" begin
        test_harmonise((5, 4), (5, 4))
        test_harmonise(Tangent{Tuple}(5, 4), (5, 4))
        test_harmonise(Tangent{Tuple}(5, 4), Tangent{Tuple}(5, 4))

        test_harmonise((5, Tangent{Tuple}(randn(5))), (5, (randn(5), )))
        test_harmonise(
            Tangent{Any}(Tangent{Any}(randn(5))),
            (Tangent{Any}(randn(5)), ),
        )
    end

    @testset "NamedTuple / Tangent{NamedTuple}" begin
        test_harmonise(Tangent{Any}(; m=4, P=5), Tangent{Gaussian}(; m=5, P=4))
        test_harmonise(Tangent{Any}(; m=4, P=5), Tangent{Any}(; m=4))
        test_harmonise(Tangent{Any}(; m=5), Tangent{Any}(; P=4))

        test_harmonise(Tangent{Any}(; m=(5, 4)), Tangent{Any}(; P=4))

        test_harmonise(Tangent{Any}(; m=5, P=4), Gaussian(5, 4))
        test_harmonise(Tangent{Any}(; P=4), Gaussian(4, 5))
    end
end


text:
module TestExamples2DHypDiff

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Hyperbolic diffusion" begin
  @testset "elixir_hypdiff_lax_friedrichs.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_lax_friedrichs.jl"),
      l2   = [0.00015687751817403066, 0.001025986772216324, 0.0010259867722164071],
      linf = [0.001198695637957381, 0.006423873515531753, 0.006423873515533529])
  end

  @testset "elixir_hypdiff_harmonic_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_harmonic_nonperiodic.jl"),
      l2   = [8.618132355121019e-8, 5.619399844384306e-7, 5.619399844844044e-7],
      linf = [1.1248618588430072e-6, 8.622436487026874e-6, 8.622436487915053e-6])
  end

  @testset "elixir_hypdiff_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic.jl"),
      l2   = [8.523077653954864e-6, 2.8779323653020624e-5, 5.454942769125663e-5],
      linf = [5.522740952468297e-5, 0.00014544895978971679, 0.00032396328684924924])
  end

  @testset "elixir_hypdiff_godunov.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_godunov.jl"),
      l2   = [5.868147556427088e-6, 3.80517927324465e-5, 3.805179273249344e-5],
      linf = [3.701965498725812e-5, 0.0002122422943138247, 0.00021224229431116015])
  end
end

end # module


text:
"""juliaでやる

単一セクション
"""

using Plots
using LinearAlgebra

using StaticArrays
using ArraysOfArrays

function split_vec_of_arrays(u)
    vec.(u) |>
    x ->
    VectorOfSimilarVectors(x).data |>
    transpose |>
    VectorOfSimilarVectors
end

function solve_RungeKutta(dx, x₀::Vector{T}, t_span, Δt::T) where T
    """ルンゲクッタ法（4次）"""

    t = range(t_span..., step = Δt)  # 時間軸
    x = Vector{typeof(x₀)}(undef, length(t))  # 解を格納する1次元配列

    x[1] = x₀  # 初期値
    for i in 1:length(x)-1
        k₁ = dx(x[i])
        k₂ = dx(x[i]+k₁*Δt/2)
        k₃ = dx(x[i]+k₂*Δt/2)
        k₄ = dx(x[i]+k₃*Δt)
        x[i+1] = x[i] .+ (k₁ .+ 2k₂ .+ 2k₃ .+k₄) .* Δt/6
    end

    t, x
end

name = "ikko_dake"

include("derived/" * name * "/eqs/julia_style/C0_0.jl")
include("derived/" * name * "/eqs/julia_style/C0_1.jl")
include("derived/" * name * "/eqs/julia_style/C0_2.jl")
include("derived/" * name * "/eqs/julia_style/C1_0.jl")
include("derived/" * name * "/eqs/julia_style/C1_1.jl")
include("derived/" * name * "/eqs/julia_style/C1_2.jl")
include("derived/" * name * "/eqs/julia_style/C2_0.jl")
include("derived/" * name * "/eqs/julia_style/C2_1.jl")
include("derived/" * name * "/eqs/julia_style/C2_2.jl")

include("derived/" * name * "/eqs/julia_style/M0_0.jl")
include("derived/" * name * "/eqs/julia_style/M0_1.jl")
include("derived/" * name * "/eqs/julia_style/M0_2.jl")
include("derived/" * name * "/eqs/julia_style/M1_0.jl")
include("derived/" * name * "/eqs/julia_style/M1_1.jl")
include("derived/" * name * "/eqs/julia_style/M1_2.jl")
include("derived/" * name * "/eqs/julia_style/M2_0.jl")
include("derived/" * name * "/eqs/julia_style/M2_1.jl")
include("derived/" * name * "/eqs/julia_style/M2_2.jl")

include("derived/" * name * "/eqs/julia_style/G0.jl")
include("derived/" * name * "/eqs/julia_style/G1.jl")
include("derived/" * name * "/eqs/julia_style/G2.jl")

include("derived/" * name * "/eqs/julia_style/Phi0.jl")
include("derived/" * name * "/eqs/julia_style/Theta0.jl")

using .C0_0
using .C0_1
using .C0_2
using .C1_0
using .C1_1
using .C1_2
using .C2_0
using .C2_1
using .C2_2

using .M0_0
using .M0_1
using .M0_2
using .M1_0
using .M1_1
using .M1_2
using .M2_0
using .M2_1
using .M2_2

using .G0
using .G1
using .G2
using .Phi0
using .Theta0

const N = 1  # セクションの数

"""慣性行列"""
function M(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    Z = Matrix{T}(undef, 3*N, 3*N)

    Z[1,1] = M0_0.f(q, q_dot, xi)
    Z[1,2] = M0_1.f(q, q_dot, xi)
    Z[1,3] = M0_2.f(q, q_dot, xi)
    Z[2,1] = M1_0.f(q, q_dot, xi)
    Z[2,2] = M1_1.f(q, q_dot, xi)
    Z[2,3] = M1_2.f(q, q_dot, xi)
    Z[3,1] = M2_0.f(q, q_dot, xi)
    Z[3,2] = M2_1.f(q, q_dot, xi)
    Z[3,3] = M2_2.f(q, q_dot, xi)

    return Z
end

"""コリオリ＋遠心力行列"""
function C(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    Z = Matrix{T}(undef, 3*N, 3*N)

    Z[1,1] = C0_0.f(q, q_dot, xi)
    Z[1,2] = C0_1.f(q, q_dot, xi)
    Z[1,3] = C0_2.f(q, q_dot, xi)
    Z[2,1] = C1_0.f(q, q_dot, xi)
    Z[2,2] = C1_1.f(q, q_dot, xi)
    Z[2,3] = C1_2.f(q, q_dot, xi)
    Z[3,1] = C2_0.f(q, q_dot, xi)
    Z[3,2] = C2_1.f(q, q_dot, xi)
    Z[3,3] = C2_2.f(q, q_dot, xi)

    return Z
end

"""重力行列"""
function G(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    Z = Vector{T}(undef, 3*N)

    Z[1] = G0.f(q, q_dot, xi)
    Z[2] = G1.f(q, q_dot, xi)
    Z[3] = G2.f(q, q_dot, xi)

    return Z
end

"""状態方程式"""
function X_dot(X::Vector{T}) where T
    q = X[1:3]
    q_dot = X[4:6]
    xi = 1.0

    τ = [0.0, 0.0, 0.0]

    inv_M = inv(M(q, q_dot, xi))

    x1_dot = q_dot
    x2_dot = -inv_M * C(q, q_dot, xi) * q_dot .+
    inv_M * (τ .- G(q, q_dot, xi))

    #println(norm([x1_dot; x2_dot]))
    return [x1_dot; x2_dot]

end

"""合ってるかテスト"""
function test()
    println("計算中...")
    TIME_SPAN = 6.0
    TIME_INTERVAL = 0.01

    q = [0.01, 0.02, 0.0]
    xi = 1.0
    q_dot = [0.0, 0.0, 0.0]

    t, X = solve_RungeKutta(
        X_dot,
        [q; q_dot],
        (0.0, TIME_SPAN),
        TIME_INTERVAL
    )

    l1, l2, l3, l1_dot, l2_dot, l3_dot = split_vec_of_arrays(X)
    fig = plot(xlim=(0, TIME_SPAN))
    plot!(fig, t, l1, label="l1")
    plot!(fig, t, l2, label="l2")
    plot!(fig, t, l3, label="l3")

    fig2 = plot(xlim=(0, TIME_SPAN))
    plot!(fig2, t, l1_dot, label="l1_dot")
    plot!(fig2, t, l2_dot, label="l2_dot")
    plot!(fig2, t, l3_dot, label="l3_dot")

    fig_I = plot(fig, fig2, layout=(1, 2), size=(800, 400))
    savefig(fig_I, "julia.png")

    # アニメ制作

end

@time test()


text:
using DataFrames, CSV
using Profile
include("../src/sdistance.jl") # 必ずダブルクオーテーション
include("../test/draw.jl")
include("../src/SignedDistanceFunction.jl")
include("../src/environments.jl")
include("../test/APT.jl")
import .Draw: parformance_graphs
using .SignedDistanceFunction
using .APT

# This script run by test.sh

# ====== Debug ======
# p = signedDistance2D("./test/mock_csv_data/interface.csv", parse(Int, ARGS[1]))
# DataFrame(p, :auto) |> CSV.write("./test/result/interface_result_n"*ARGS[1]*".csv", header=false)
# p = signedDistance2D("./test/mock_csv_data/interface.csv", parse(Int, ARGS[1]), "multi")
# DataFrame(p, :auto) |> CSV.write("./test/result/interface_floodfill_result_n"*ARGS[1]*".csv", header=false)
# p = signedDistance2D("./test/mock_csv_data/multiple_curves.csv", parse(Int, ARGS[1]), "multi")
# DataFrame(p, :auto) |> CSV.write("./test/result/multiple_curves_result_n"*ARGS[1]*".csv", header=false)

# ====== profiling =======

# @profile signedDistance2D("./test/mock_csv_data/interface.csv",parse(Int, ARGS[1]))
# @profile signedDistance2D("./test/mock_csv_data/interface.csv",parse(Int, ARGS[1]), "multi")
# Profile.print()
# open("prof.txt", "w") do s
#     Profile.print(IOContext(s, :displaysize => (24, 500)))
# end

# ====== memory size =======
# p = @allocated signedDistance2D("./test/mock_csv_data/interface.csv", parse(Int, ARGS[1]), "multi")
# p = @allocated signedDistance2D("./test/mock_csv_data/interface.csv", parse(Int, ARGS[1]), "multi")
p = @allocated signedDistance2D("./test/mock_csv_data/multiple_curves.csv", parse(Int, ARGS[1]), "multi")
println("\nUsed memory size: ",p/(1024*1024), " MB")

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function