text:
module DirectQhull

# MIT License
#
# Copyright (c) 2021 Juha Tapio Heiskala
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# functions "qh_get_extremes_2d" and "qh_get_simplex_facet_arrays" defined later
# licensed under BSD license from Scipy

import Qhull_jll

import Base.getproperty
import Base.iterate
import Base.getindex
import Base.length

# define Qhull types

QHboolT = Cuint
QHrealT = Cdouble
QHcoordT = QHrealT
QHpointT = QHcoordT
QHintT = Cint
QHprintT = Cint
QHcharT = Cchar
QHuintT = Cuint
QHulongT = Culong
QHcenterT = Cint
QHdoubleT = Cdouble
QHvoidT = Cvoid
QHfileT = Cvoid

export ConvexHull

# holds pointer to qhT structure
mutable struct qhT
end

# defines clockwise or counterclockwise orientation e.g. 2d-convex hull vertex ordering
const qh_ORIENTclock = 0

const qh_lib = Qhull_jll.get_libqhull_r_path()

function qh_alloc_qh(err_file::Ptr{Cvoid}=C_NULL)
    qh_ptr = ccall((:qh_alloc_qh, qh_lib), Ptr{qhT}, (Ptr{Cvoid},), err_file)
    (qh_ptr != C_NULL) ? qh_ptr : throw(ErrorException("qhT initialization failure."))
end

# This should always be true (tested for the QHsetT so that Union{Int, Pointers} as QHsetT element is valid)
@assert(sizeof(Ptr{Cvoid}) == sizeof(Int))

abstract type QHsetelemT end

# this is not an exact representation of QHull's set type
# Qhull set is dynamically allocated so exact representation is not straighforward as Julia type.
# (would require NTuple type for elements with tuple length equal to number elements in the list)
# so the set elements are wrapped by Julia array
# QH set element type is defined as union of void pointer and integer
struct QHsetT{T<:Union{QHintT, Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}}}
    maxsize::QHintT          # /* maximum number of elements (except NULL) */
    e::Array{T}              # /* array of pointers, tail is NULL */
                             # /* last slot (unless NULL) is actual size+1
                             # /* e[maxsize]==NULL or e[e[maxsize]-1]==NULL */
                             # /* this may generate a warning since e[] contains  maxsize elements */
    function QHsetT{T}(ptr::Ptr{QHsetT{T}}) where T<:Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}
        max_size = unsafe_load(Ptr{QHintT}(ptr))
        # with passing C_NULL as qh_ptr below the call will crash, if the setsize is invalid
        # however, qhull would anyway exit in this case with internal error, though in a more gracious manner
        # (this removes the need to pass the qhT pointer to this constructor)
        set_size = qh_setsize(Ptr{qhT}(C_NULL), ptr)
        # assumed here QHsetT field e is offset Ptr-size from the maxsize field
        ptr_array = unsafe_wrap(Array, Ptr{T}(ptr+sizeof(Ptr)), set_size)
        new(max_size, ptr_array)
    end
end

function Base.getindex(set::QHsetT{T}, idx::Int) where T<:Union{QHintT, Ptr{<:QHsetelemT}}
    ptr = set.e[idx]
    unsafe_load(ptr)
end

function Base.length(set::QHsetT{T}) where T<:Union{QHintT, Ptr{<:QHsetelemT}}
    return length(set.e)
end

# QH set as providing raw pointers
struct QHsetPtrT{T<:Union{QHintT, Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}}}
    maxsize::QHintT          # /* maximum number of elements (except NULL) */
    e::Array{T}              # /* array of pointers, tail is NULL */
                             # /* last slot (unless NULL) is actual size+1
                             # /* e[maxsize]==NULL or e[e[maxsize]-1]==NULL */
                             # /* this may generate a warning since e[] contains  maxsize elements */
    function QHsetPtrT{T}(ptr::Ptr{QHsetPtrT{T}}) where T<:Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}
        max_size = unsafe_load(Ptr{QHintT}(ptr))
        # with passing C_NULL as qh_ptr below the call will crash, if the setsize is invalid
        # however, qhull would anyway exit in this case with internal error, though in a more gracious manner
        # (this removes the need to pass the qhT pointer to this constructor)
        set_size = qh_setsize(Ptr{qhT}(C_NULL), ptr)
        # assumed here QHsetT field e is offset Ptr-size from the maxsize field
        ptr_array = unsafe_wrap(Array, Ptr{T}(ptr+sizeof(Ptr)), set_size)
        new(max_size, ptr_array)
    end
end

function Base.getindex(set::QHsetPtrT{T}, idx::Int) where T<:Union{QHintT, Ptr{<:QHsetelemT}}
    ptr = set.e[idx]
    ptr
end

function Base.length(set::QHsetPtrT{T}) where T<:Union{QHintT, Ptr{<:QHsetelemT}}
    return length(set.e)
end

# Qhull vertex type
# Comments are from original qhull code.
mutable struct QHvertexT{HD} <: QHsetelemT
    next::Ptr{QHvertexT{HD}}                # /* next vertex in vertex_list or vertex_tail */
    previous::Ptr{QHvertexT{HD}}            # /* previous vertex in vertex_list or NULL, for C++ interface */
    point::Ptr{NTuple{HD, QHpointT}}        # /* hull_dim coordinates (coordT) */
    neighbors::Ptr{QHsetT{<:QHsetelemT}}    # /* neighboring facets of vertex, qh_vertexneighbors()
                                            # initialized in io_r.c or after first merge
                              # qh_update_vertices for qh_addpoint or qh_triangulate
                              # updated by merges
                              # qh_order_vertexneighbors by 2-d (orientation) 3-d (adjacency), n-d (f.visitid,id) */
    id::QHuintT               # /* unique identifier, 1..qh.vertex_id,  0 for sentinel, printed as 'r%d' */
    visitid::QHuintT          # /* for use with qh.vertex_visit, size must match */
    flags::QHcharT
    # seen:1;      /* used to perform operations only once */
    #flagT    seen2:1;     /* another seen flag */
    #flagT    deleted:1;   /* vertex will be deleted via qh.del_vertices */
    #flagT    delridge:1;  /* vertex belonged to a deleted ridge, cleared by qh_reducevertices */
    #flagT    newfacet:1;  /* true if vertex is in a new facet
    #                       vertex is on qh.newvertex_list and it has a facet on qh.newfacet_list
    #                       or vertex is on qh.newvertex_list due to qh_newvertices while merging
    #                       cleared by qh_resetlists */
    #flagT    partitioned:1; /* true if deleted vertex has been partitioned */
end

# Qhull facet type
# Comments are from original qhull code.
mutable struct QHfacetT{HD} <: QHsetelemT
    furthestdist::QHcoordT  # distance to furthest point of outsideset
    maxoutside::QHcoordT    # max computed distance of point to facet
                            # Before QHULLfinished this is an approximation
                            # since maxdist not always set for qh_mergefacet
                            # Actual outer plane is +DISTround and
                            # computed outer plane is +2*DISTround.
                            # Initial maxoutside is qh.DISTround, otherwise distance tests need to account for DISTround */

    offset::QHcoordT        # exact offset of hyperplane from origin
    normal::Ptr{NTuple{HD, QHcoordT}}   # normal of hyperplane, hull_dim coefficients
                            # if f.tricoplanar, shared with a neighbor

    #union {                # in order of testing */
    area::QHrealT           # area of facet, only in io_r.c if  f.isarea */
    #facetT *replace    # replacement facet for qh.NEWfacets with f.visible
                            # NULL if qh_mergedegen_redundant, interior, or !NEWfacets */
                            # facetT *samecycle;   /* cycle of facets from the same visible/horizon intersection,
                            # if ->newfacet */
    # facetT *newcycle;    /*  in horizon facet, current samecycle of new facets */
    # facetT *trivisible;  /* visible facet for ->tricoplanar facets during qh_triangulate() */
    # facetT *triowner;    /* owner facet for ->tricoplanar, !isarea facets w/ ->keepcentrum */
    # }f;
    center::Ptr{QHcoordT}  # set according to qh.CENTERtype */
                           # qh_ASnone:    no center (not MERGING) */
                           # qh_AScentrum: centrum for testing convexity (qh_getcentrum) */
                           #               assumed qh_AScentrum while merging */
                           # qh_ASvoronoi: Voronoi center (qh_facetcenter) */
                           # after constructing the hull, it may be changed (qh_clearcenter) */
                           # if tricoplanar and !keepcentrum, shared with a neighbor */
    previous::Ptr{QHfacetT{HD}} # previous facet in the facet_list or NULL, for C++ interface */
    next::Ptr{QHfacetT{HD}}     # next facet in the facet_list or facet_tail */
    vertices::Ptr{QHsetT{Ptr{QHvertexT{HD}}}}   # vertices for this facet, inverse sorted by ID
                            # if simplicial, 1st vertex was apex/furthest
                            # qh_reduce_vertices removes extraneous vertices via qh_remove_extravertices
                            # if f.visible, vertices may be on qh.del_vertices */
    ridges::Ptr{QHsetT}     # explicit ridges for nonsimplicial facets or nonsimplicial neighbors.
                            # For simplicial facets, neighbors define the ridges
                            # qh_makeridges() converts simplicial facets by creating ridges prior to merging
                            # If qh.NEWtentative, new facets have horizon ridge, but not vice versa
                            # if f.visible && qh.NEWfacets, ridges is empty */
    neighbors::Ptr{QHsetT{Ptr{QHfacetT{HD}}}}  # neighbors of the facet.  Neighbors may be f.visible
                            # If simplicial, the kth neighbor is opposite the kth vertex and the
                            # first neighbor is the horizon facet for the first vertex.
                            # dupridges marked by qh_DUPLICATEridge (0x01) and qh_MERGEridge (0x02)
                            # if f.visible && qh.NEWfacets, neighbors is empty */
    outsideset::Ptr{QHsetT} # set of points outside this facet
                            # if non-empty, last point is furthest
                            # if NARROWhull, includes coplanars (less than qh.MINoutside) for partitioning*/
    coplanarset::Ptr{QHsetT{Ptr{NTuple{HD, QHpointT}}}} # set of points coplanar with this facet
                             # >= qh.min_vertex and <= facet->max_outside
                             # a point is assigned to the furthest facet
                             # if non-empty, last point is furthest away */
    visitid::QHuintT         # visit_id, for visiting all neighbors, all uses are independent */
    id::QHuintT              # unique identifier from qh.facet_id, 1..qh.facet_id, 0 is sentinel, printed as 'f%d' */
    flags::QHuintT

    # unsigned int nummerge:9; /* number of merges */
    # define qh_MAXnummerge 511 /* 2^9-1 */
    #                        /* 23 flags (at most 23 due to nummerge), printed by "flags:" in io_r.c */
    # flagT    tricoplanar:1; /* True if TRIangulate and simplicial and coplanar with a neighbor */
    #                      /*   all tricoplanars share the same apex */
    #                      /*   all tricoplanars share the same ->center, ->normal, ->offset, ->maxoutside */
    #                      /*     ->keepcentrum is true for the owner.  It has the ->coplanareset */
    #                      /*   if ->degenerate, does not span facet (one logical ridge) */
    #                      /*   during qh_triangulate, f.trivisible points to original facet */
    # flagT    newfacet:1;  /* True if facet on qh.newfacet_list (new/qh.first_newfacet or merged) */
    # flagT    visible:1;   /* True if visible facet (will be deleted) */
    # flagT    toporient:1; /* True if created with top orientation
    #                       after merging, use ridge orientation */
    # flagT    simplicial:1;/* True if simplicial facet, ->ridges may be implicit */
    # flagT    seen:1;      /* used to perform operations only once, like visitid */
    # flagT    seen2:1;     /* used to perform operations only once, like visitid */
    # flagT    flipped:1;   /* True if facet is flipped */
    # flagT    upperdelaunay:1; /* True if facet is upper envelope of Delaunay triangulation */
    # flagT    notfurthest:1; /* True if last point of outsideset is not furthest */
    #
    # /*-------- flags primarily for output ---------*/
    # flagT    good:1;      /* True if a facet marked good for output */
    # flagT    isarea:1;    /* True if facet->f.area is defined */

    # /*-------- flags for merging ------------------*/
    # flagT    dupridge:1;  /* True if facet has one or more dupridge in a new facet (qh_matchneighbor),
    #                         a dupridge has a subridge shared by more than one new facet */
    # flagT    mergeridge:1; /* True if facet or neighbor has a qh_MERGEridge (qh_mark_dupridges)
    #                        ->normal defined for mergeridge and mergeridge2 */
    # flagT    mergeridge2:1; /* True if neighbor has a qh_MERGEridge (qh_mark_dupridges) */
    # flagT    coplanarhorizon:1;  /* True if horizon facet is coplanar at last use */
    # flagT     mergehorizon:1; /* True if will merge into horizon (its first neighbor w/ f.coplanarhorizon). */
    # flagT     cycledone:1;/* True if mergecycle_all already done */
    # flagT    tested:1;    /* True if facet convexity has been tested (false after merge */
    # flagT    keepcentrum:1; /* True if keep old centrum after a merge, or marks owner for ->tricoplanar
    #                          Set by qh_updatetested if more than qh_MAXnewcentrum extra vertices
    #                          Set by qh_mergefacet if |maxdist| > qh.WIDEfacet */
    # flagT    newmerge:1;  /* True if facet is newly merged for reducevertices */
    # flagT    degenerate:1; /* True if facet is degenerate (degen_mergeset or ->tricoplanar) */
    # flagT    redundant:1;  /* True if facet is redundant (degen_mergeset)
    #                     Maybe merge degenerate and redundant to gain another flag */
end

# Iteration for facet list
function iterate(first_fct::QHfacetT{HD}) where HD
    return (first_fct, first_fct.next_ptr)
end

function iterate(first_fct::QHfacetT{HD}, next_fct_ptr::Ptr{QHfacetT{HD}}) where HD
    next_fct = unsafe_load(next_fct_ptr)
    ## ID=0 is dummy facet that indicates end of the list
    if next_fct.id == QHuintT(0)
        return nothing
    else
        return (next_fct, next_fct.next_ptr)
    end
end

# Iteration for vertex list
function iterate(first_vtx::QHvertexT{HD}) where HD
    return (first_vtx, first_vtx.next_ptr)
end

function iterate(first_vtx::QHvertexT{HD}, next_vtx_ptr::Ptr{QHvertexT{HD}}) where HD
    next_vtx = unsafe_load(next_vtx_ptr)
    ## ID=0 is dummy vertext that indicates end of the list
    if next_vtx.id == QHuintT(0)
        return nothing
    else
        return (next_vtx, next_vtx.next_ptr)
    end
end

# Iterate QHsetT with pointer types
function iterate(set::QHsetT{T}) where T<:Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}

    if length(set.e) == 0
        return nothing
    else
        return (unsafe_load(set.e[1]), 2)
    end
end

function iterate(set::QHsetT{T}, idx::Int) where T<:Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}
    if idx > length(set.e)
        return nothing
    else
        return (unsafe_load(set.e[idx]), idx+1)
    end
end

# Iterate QHsetT with pointer types
function iterate(set::QHsetPtrT{T}) where T<:Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}

    if length(set.e) == 0
        return nothing
    else
        return (set.e[1], 2)
    end
end

function iterate(set::QHsetPtrT{T}, idx::Int) where T<:Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}
    if idx > length(set.e)
        return nothing
    else
        return (set.e[idx], idx+1)
    end
end

# Build convex hull from a set points
struct ConvexHull
    qh_ptr::Ptr{qhT}
    points::Matrix{QHcoordT}
    vertices::Vector{QHuintT}
    simplices::Matrix{QHuintT}
    neighbors::Matrix{QHuintT}
    equations::Matrix{QHrealT}
    coplanar::Matrix{QHintT}
    good::Vector{QHboolT}
    area::QHrealT
    volume::QHrealT
    max_bound::Vector{QHrealT}
    min_bound::Vector{QHrealT}

    function ConvexHull(pnts::Matrix{QHcoordT}, qhull_options::Vector{String}=Vector{String}())
        qh_ptr = qh_alloc_qh()

        pushfirst!(qhull_options, "Qt")

        if size(pnts, 2)>=5
            push!(qhull_options, "Qx")
        end

        # make options string
        qh_opts_str = foldl((l,r)->l*" "*r, qhull_options)

        # calculate new qhull
        res = qh_new_qhull(qh_ptr, pnts, qh_opts_str)

        hd = qh_get_hull_dim(qh_ptr)

        if hd == 2
            vertices = qh_get_extremes_2d(qh_ptr)
        else
            vertices = unique((vtx)->qh_pointid(qh_ptr, vtx.point_ptr)+1)
        end

        # collect convex hull points
        # for some reason using hd from above in Val(hd) crashes Julia
        simplices = qh_get_convex_hull_pnts(qh_ptr, Val(size(pnts,1)))

        # facet neighbors, equations, good
        neighbors, equations, coplanar, good = qh_get_simplex_facet_arrays(qh_ptr, Val(size(pnts,1)))

        if ("QG" in qhull_options || "QG4" in qhull_options)
            Bool.(good)
        else
            good = Vector{QHboolT}()
        end

        # calculate total area and volume
        qh_getarea(qh_ptr, Val(size(pnts,1)))
        area = qh_get_totarea(qh_ptr)
        vol = qh_get_totvol(qh_ptr)

        # max and min bounds
        max_bound = maximum(pnts, dims=2)[:]
        min_bound = minimum(pnts, dims=2)[:]

        # the new Qhull value
        new(qh_ptr, pnts, vertices, simplices, neighbors, equations, coplanar,
            good, area, vol, max_bound, min_bound)
    end
end

# calculate new convex hull from the given points and Qhull options
function qh_new_qhull(qh::Ptr{qhT}, pnts::StridedMatrix{Float64}, opts::String)
    ok = ccall((:qh_new_qhull, qh_lib), Cint,
               (Ptr{qhT}, Cint, Cint, Ref{QHcoordT}, QHboolT, Ptr{QHcharT}, Ptr{QHfileT}, Ptr{QHfileT}),
               qh, size(pnts, 1), size(pnts, 2), pnts, false, "qhull " * opts, C_NULL, C_NULL)

    return ok
end

function qh_triangulate(qh::Ptr{qhT})
    ccall((:qh_triangulate, qh_lib), Cvoid, (Ptr{qhT},), qh)
end

function qh_getarea(qh::Ptr{qhT}, ::Val{HD}) where HD
    ccall((:qh_getarea, qh_lib), Cvoid, (Ptr{qhT}, Ptr{QHfacetT}), qh, qh_get_facet_list_ptr(qh, Val(HD)))
end

# retrieve Qhull internal point id
function qh_pointid(qh::Ptr{qhT}, pnt::Ptr{NTuple{N, QHpointT}}) where N
    id = ccall((:qh_pointid, qh_lib), Cuint,
               (Ptr{qhT}, Ptr{QHpointT}), qh, pnt)
    return id
end
function qh_pointid(qh::Ptr{qhT}, pnt::Ptr{QHpointT}) where N
    id = ccall((:qh_pointid, qh_lib), Cuint,
               (Ptr{qhT}, Ptr{QHpointT}), qh, pnt)
    return id
end

# retrieve Qhull internal point id
function qh_nearvertex(qh::Ptr{qhT}, fct::Ptr{QHfacetT{HD}}, pnt::Ptr{QHpointT}, dist::Array{QHrealT,1}) where HD
    return ccall((:qh_nearvertex, qh_lib), Ptr{QHvertexT{HD}},
                 (Ptr{qhT}, Ptr{QHfacetT{HD}}, Ptr{QHpointT}, Ptr{QHrealT}), qh, fct, pnt, dist)
end

# retrieve Qhull set size
function qh_setsize(qh::Ptr{qhT}, set::Ptr{QHsetT{T}}) where T<:Union{QHintT, Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}}
    ccall((:qh_setsize, qh_lib), Cint, (Ptr{qhT}, Ptr{QHsetT}), qh, set)
end
function qh_setsize(qh::Ptr{qhT}, set::Ptr{QHsetPtrT{T}}) where T<:Union{QHintT, Ptr{<:Union{QHsetelemT, NTuple, QHpointT}}}
    ccall((:qh_setsize, qh_lib), Cint, (Ptr{qhT}, Ptr{QHsetPtrT}), qh, set)
end

# GETTER accessors for plain data types
for (T, getter) in ((:QHintT, :qh_get_hull_dim), (:QHintT, :qh_get_num_facets), (:QHintT, :qh_get_num_points),
                    (:QHintT, :qh_get_num_vertices), (:QHintT, :qh_get_visit_id), (:QHintT, :qh_get_vertex_visit),
                    (:QHrealT, :qh_get_totarea), (:QHrealT, :qh_get_totvol), (:QHuintT, :qh_get_facet_id))
    @eval begin
        function ($getter)(qh::Ptr{qhT})
            ccall(($(QuoteNode(getter)), qh_lib), $T, (Ptr{qhT},), qh)
        end
    end
end

# GETTER for facet list pointer
function qh_get_facet_list_ptr(qh::Ptr{qhT}, ::Val{N}) where N
    ccall((:qh_get_facet_list, qh_lib), Ptr{QHfacetT{N}}, (Ptr{qhT},), qh)
end

# GETTER for facet list as Julia QHfacetT type
@noinline function qh_get_facet_list(qh::Ptr{qhT}, ::Val{N}) where N
    local ptr

    # this is workaround for Julia crashes that seem to happen if Val(N) has N defined in runtime
    if N==1
        ptr = qh_get_facet_list_ptr(qh, Val(1))
    elseif N==2
        ptr = qh_get_facet_list_ptr(qh, Val(2))
    elseif N==3
        ptr = qh_get_facet_list_ptr(qh, Val(3))
    elseif N==4
        ptr = qh_get_facet_list_ptr(qh, Val(4))
    elseif N==5
        ptr = qh_get_facet_list_ptr(qh, Val(5))
    else
        error("Too large hull dimension")
    end

    if ptr != C_NULL
        return unsafe_load(ptr)
    else
        return nothing
    end
end

# GETTER for facet list pointer
function qh_get_vertex_list_ptr(qh::Ptr{qhT}, ::Val{N}) where N
    ccall((:qh_get_vertex_list, qh_lib), Ptr{QHvertexT{N}}, (Ptr{qhT},), qh)
end

# GETTER for facet list as Julia QHfacetT type
@noinline function qh_get_vertex_list(qh::Ptr{qhT}, ::Val{N}) where N
    local ptr
    # this is workaround for Julia crashes that seem to happen if Val(N) has N defined in runtime
    if N==1
        ptr = qh_get_vertex_list_ptr(qh, Val(1))
    elseif N==2
        ptr = qh_get_vertex_list_ptr(qh, Val(2))
    elseif N==3
        ptr = qh_get_vertex_list_ptr(qh, Val(3))
    elseif N==4
        ptr = qh_get_vertex_list_ptr(qh, Val(4))
    elseif N==5
        ptr = qh_get_vertex_list_ptr(qh, Val(5))
    else
        error("Too large hull dimension")
    end
    if ptr != C_NULL
        return unsafe_load(ptr)
    else
        return nothing
    end
end

function Base.getproperty(qh::ConvexHull, fld::Symbol)
    if fld === :hull_dim
        return qh_get_hull_dim(qh.qh_ptr)
    elseif fld === :num_facets
        return qh_get_num_facets(qh.qh_ptr)
    elseif fld === :num_points
        return qh_get_num_points(qh.qh_ptr)
    elseif fld === :num_vertices
        return qh_get_num_vertices(qh.qh_ptr)
    elseif fld === :visit_id
        return qh_get_visit_id(qh.qh_ptr)
    elseif fld === :vertex_visit
        return qh_get_vertex_visit(qh.qh_ptr)
    elseif fld === :facet_list
        # Hull dimension given to facet type so that hull dimension array size is known
        return qh_get_facet_list(qh.qh_ptr, Val(qh.hull_dim))
    elseif fld === :vertex_list
        # Hull dimension given to facet type so that hull dimension array size is known
        return qh_get_vertex_list(qh.qh_ptr, Val(qh.hull_dim))
    else
        return getfield(qh, fld)
    end
end

function Base.getproperty(fct::QHfacetT{HD}, fld::Symbol) where HD
    if fld === :next
        ptr = getfield(fct, :next)
        if ptr == C_NULL
            return nothing
        else
            return unsafe_load(ptr)
        end
    elseif fld === :next_ptr
        return getfield(fct, :next)
    elseif fld === :self_ptr
        # ID=0 is dummy facet, with next=C_NULL
        # (so the facet cannot be updated in this case)
        if (fct.id != 0)
            return fct.next.previous
        else
            throw(ErrorException("Cannot get self pointer for last facet in facet list."))
        end
    elseif fld === :vertices
        return QHsetT{Ptr{QHvertexT{HD}}}(fct.vertices_ptr)
    elseif fld === :vertices_ptr
        return getfield(fct, :vertices)
    elseif fld === :neighbors
        return QHsetT{Ptr{QHfacetT{HD}}}(fct.neighbors_ptr)
    elseif fld === :neighbors_ptr
        return getfield(fct, :neighbors)
    elseif fld === :coplanarset
        ptr = fct.coplanarset_ptr
        if ptr != C_NULL
            return QHsetT{Ptr{NTuple{HD, QHpointT}}}(ptr)
        else
            return nothing
        end
    elseif fld === :coplanarset_ptr
        return getfield(fct, :coplanarset)
    elseif fld === :coplanar_ptr_set
        ptr = fct.coplanarset_ptr
        if ptr != C_NULL
            return QHsetPtrT{Ptr{QHpointT}}(Ptr{QHsetPtrT{Ptr{QHpointT}}}(ptr))
        else
            return nothing
        end
    elseif fld === :normal
        return unsafe_load(fct.normal_ptr)
    elseif fld === :normal_ptr
        return getfield(fct, :normal)
    elseif fld == :toporient
        return QHboolT( (getfield(fct, :flags)>>12)&1)  # toporient is 13th bit in the flags field
    elseif fld == :simplicial
        return QHboolT( (getfield(fct, :flags)>>13)&1)  # simplicial is 14th bit in the flags field
    elseif fld == :good
        return QHboolT( (getfield(fct, :flags)>>19)&1)  # good is 20th bit in the flags field
    else
        return getfield(fct, fld)
    end
end

function Base.setproperty!(fct::QHfacetT{HD}, fld::Symbol, value) where HD
    setfield!(fct, fld, value)
    # store back to qhull, pointer to self is obtained (somewhat dangerously) utilizing the facet linked list
    unsafe_store!(fct.self_ptr, fct)
    return value
end

function Base.getproperty(vtx::QHvertexT, fld::Symbol)
    if fld === :next
        ptr = getfield(vtx, :next)
        if ptr == C_NULL
            return nothing
        else
            return unsafe_load(ptr)
        end
    elseif fld === :next_ptr
        return getfield(vtx, :next)
    elseif fld === :point
        ptr = getfield(vtx, :point)
        if ptr == C_NULL
            return nothing
        else
            return unsafe_load(ptr)
        end
    elseif fld === :point_ptr
        return getfield(vtx, :point)
    else
        return getfield(vtx, fld)
    end
end

function Base.setproperty!(vtx::QHvertexT{HD}, fld::Symbol, value) where HD
    setfield!(vtx, fld, value)
    # store back to qhull, pointer to self is obtained (somewhat dangerously) utilizing the vertex linked list
    if (vtx.id != 0) # ID=0 is dummy vertex, with next=C_NULL (so the vertex cannot be updated in this case)
        unsafe_store!(vtx.next.previous, vtx)
    end
end

# get calculated convex hull points as Julia Int Array
function qh_get_convex_hull_pnts(qh_ptr::Ptr{qhT}, ::Val{HD}) where HD
    n_facets = qh_get_num_facets(qh_ptr)

    facet_list = qh_get_facet_list(qh_ptr, Val(HD))

    # convex hull point ids (i.e. indexes to input points to qhull)
    pnts = Matrix{QHuintT}(undef, HD, n_facets)

    facet_ix = 1
    for facet in facet_list
        vtxSet = facet.vertices

        vtx_ix = 1
        for vtx in vtxSet
            pnt_id = qh_pointid(qh_ptr, vtx.point_ptr)

            # +1 to change to 1-based index for Julia from C 0-based index
            pnts[vtx_ix, facet_ix] = pnt_id + 1
            vtx_ix+=1
        end
        facet_ix+=1
    end

    return pnts
end

# get calculated convex hull vertices as Julia Int Array
function qh_get_convex_hull_vertices(qh_ptr::Ptr{qhT}, ::Val{HD}) where HD
    n_vertices = qh_get_num_vertices(qh_ptr)

    vertex_list = qh_get_vertex_list(qh_ptr, Val(HD))

    # convex hull point ids (i.e. indexes to input points to qhull)
    vertices = Vector{QHuintT}(undef, n_vertices)

    vertex_ix = 1
    for vtx in vertex_list
        pnt_id = qh_pointid(qh_ptr, vtx.point_ptr)
        vertices[vertex_ix] = pnt_id + 1
        vertex_ix+=1
    end

    return vertices
end

# Below functions "qh_get_extremes_2d" and "qh_get_simplex_facet_arrays" are
# adapted from Qhull/io.c and Scipy/_qhull.pyx/get_extremes_2d with the
# below BSD license from _qhull.pyx/Scipy
#
# Copyright (C)  Pauli Virtanen, 2010.
#
# Distributed under the same BSD license as Scipy.
#
#
# Copyright (c) 2001-2002 Enthought, Inc.  2003-2019, SciPy Developers.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above
#   copyright notice, this list of conditions and the following
#   disclaimer in the documentation and/or other materials provided
#   with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its
#   contributors may be used to endorse or promote products derived
#   from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

function qh_get_extremes_2d(qh_ptr::Ptr{qhT})

    # qhull io.c has the below call, Scipy does not
    # This is called from ConvexHull constructor, hence assumed visit_id is valid
    # qh_countfacets(facetlist, facets, printall, &numfacets, &numsimplicial,
    #    &totneighbors, &numridges, &numcoplanars, &numtricoplanars); /* marks qh visit_id */
    # vertices= qh_facetvertices(facetlist, facets, printall);
    # qh_fprintf(fp, 9088, "%d\n", qh_setsize(vertices));
    # qh_settempfree(&vertices);

    if (qh_get_num_facets(qh_ptr) == 0)
        return Vector{QHintT}()
    end

    # qhull/Scipy update the internal qhT structure fields for the visit ids
    # as the values are only read in the rest of the code, local variables used here
    # (assumed the qh internal value update is not necessary)
    qh_vertex_visit = qh_get_vertex_visit(qh_ptr) + 1
    qh_visit_id = qh_get_visit_id(qh_ptr) + 1

    # Init result array
    extremes = zeros(QHintT, 100)
    n_extremes = 0

    # get first facet in facet list
    facet = qh_get_facet_list(qh_ptr, Val(qh_get_hull_dim(qh_ptr)))

    # use facet id instead of pointer comparision for ending the while loop
    start_facet_id = facet.id

    while !isnothing(facet)
        if facet.visitid == qh_visit_id
            throw(ErrorException("Internal Qhull error, loop in facet list"))
        end

        if xor(facet.toporient, qh_ORIENTclock) != 0
            vertexA = facet.vertices[1]
            vertexB = facet.vertices[2]
            nextfacet = facet.neighbors[1]
        else
            vertexA = facet.vertices[2]
            vertexB = facet.vertices[1]
            nextfacet = facet.neighbors[2]
        end

        # check if result array needs resizing
        if n_extremes + 2 > length(extremes)
            resize!(extremes, 2*length(extremes)+1)
        end
        if facet.visitid != 0  # qhull has this check, Scipy does not
            if (vertexA.visitid != qh_vertex_visit)
                # this updates also internal qhull vertex state
                vertexA.visitid = QHuintT(qh_vertex_visit)
                n_extremes += 1
                extremes[n_extremes] = qh_pointid(qh_ptr, vertexA.point_ptr) + 1 # +1 to change to 1-based index
            end
            if vertexB.visitid != qh_vertex_visit
                vertexB.visitid = QHuintT(qh_vertex_visit)
                n_extremes += 1
                extremes[n_extremes] = qh_pointid(qh_ptr, vertexB.point_ptr) + 1 # +1 to change to 1-based index
            end
        end
        # this updates also internal qhull facet state
        facet.visitid = QHuintT(qh_visit_id);
        facet = nextfacet

        if facet.id == start_facet_id
            break
        end
    end

    resize!(extremes, n_extremes)
    return extremes
end

# get calculated convex hull points as Julia Int Array
function qh_get_simplex_facet_arrays(qh_ptr::Ptr{qhT}, ::Val{HD}) where HD
    n_facets = qh_get_num_facets(qh_ptr)

    facet_list = qh_get_facet_list(qh_ptr, Val(HD))

    id_map = fill!(Vector{QHintT}(undef, qh_get_facet_id(qh_ptr)), QHintT(-1))

    j = 1 # 1 based index in Julia
    for facet in facet_list
        if (facet.simplicial != 0) &&
            (length(facet.vertices) != HD) || (length(facet.neighbors) != HD)
            throw(ErrorException("Non-simplical facet encountered."))
        end
        id_map[facet.id] = j
        j += 1
    end

    # facet neighors
    neighbors = Matrix{QHuintT}(undef, HD, n_facets)
    facets = Matrix{QHintT}(undef, HD, n_facets)
    good = Vector{QHintT}(undef, n_facets)
    equations = Matrix{QHrealT}(undef, HD+1, n_facets)
    coplanar = zeros(QHintT, 10, 3)
    ncoplanar = 1

    facet_ix = 1
    for facet in facet_list
        neighborSet = facet.neighbors

        neighbor_ix = 1
        for neighbor in neighborSet
            neighbors[neighbor_ix, facet_ix] = id_map[neighbor.id]
            neighbor_ix+=1
        end

        for ix in 1:HD
            equations[ix, facet_ix] = facet.normal[ix]
        end
        equations[HD+1, facet_ix] = facet.offset

        dist = Array{QHrealT}(undef, 1)

        # Save coplanar info
        if !isnothing(facet.coplanar_ptr_set)
            for point_ptr in facet.coplanar_ptr_set
                vertex = unsafe_load(qh_nearvertex(qh_ptr, facet.self_ptr, point_ptr, dist))
                if ncoplanar >= size(coplanar, 1)
                    # The array is always safe to resize
                    coplanar = cat(coplanar, zeros(QHintT, ncoplanar + 1, 3), dims=1)
                end
                coplanar[ncoplanar, 1] = qh_pointid(qh_ptr, point_ptr)
                coplanar[ncoplanar, 2] = id_map[facet.id]
                coplanar[ncoplanar, 3] = qh_pointid(qh_ptr, vertex.point_ptr)
                ncoplanar += 1
            end
        end

        # save good info
        good[facet_ix] = facet.good
        facet_ix+=1
    end

    # resize
    coplanar = coplanar[1:ncoplanar-1, :]

    return (neighbors, equations, coplanar, good)
end

end


text:
using ITensors
using Printf
using Random

# Use DMRG to solve the spin 1, 1D Heisenberg model on 100 sites
# For the Heisenberg model in one dimension
# H = J ∑ᵢ(SᶻᵢSᶻᵢ₊₁ + SˣᵢSˣᵢ₊₁ + SʸᵢSʸᵢ₊₁ )
#   = J ∑ᵢ[SᶻᵢSᶻᵢ₊₁ + ½(S⁺ᵢS⁻ᵢ₊₁ + S⁻ᵢS⁺ᵢ₊₁)]
# We'll work in units where J=1

let
  N = 100

  # Create N spin-one degrees of freedom
  sites = siteinds("S=1",N)
  # Alternatively can make spin-half sites instead
  #sites = siteinds("S=1/2",N)

  # Input operator terms which define a Hamiltonian
  ampo = AutoMPO()
  for j=1:N-1
    ampo += "Sz",j,"Sz",j+1
    ampo += 0.5,"S+",j,"S-",j+1
    ampo += 0.5,"S-",j,"S+",j+1
  end
  # Convert these terms to an MPO tensor network
  H = MPO(ampo,sites)

  # Create an initial random matrix product state
  psi0 = randomMPS(sites,10)

  # Plan to do 5 DMRG sweeps:
  sweeps = Sweeps(5)
  # Set maximum MPS bond dimensions for each sweep
  maxdim!(sweeps, 10,20,100,100,200)
  # Set maximum truncation error allowed when adapting bond dimensions
  cutoff!(sweeps, 1E-11)
  @show sweeps

  # Run the DMRG algorithm, returning energy and optimized MPS
  energy, psi = dmrg(H,psi0, sweeps)
  @printf("Final energy = %.12f\n",energy)
end

ext: jl


text:
function main()
    country_labels = Dict{String,Int64}()
    city_labels = Dict{String,Int64}()
    open("locations-TrivagoClickout.txt") do f
        open("city-labels.txt", "w") do g1
            open("country-labels.txt", "w") do g2
                for line in eachline(f)
                    city = strip(line)
                    if !haskey(city_labels, city)
                        city_labels[city] = length(city_labels) + 1
                    end
                    country = strip(split(line, ",")[2])
                    if !haskey(country_labels, country)
                        country_labels[country] = length(country_labels) + 1
                    end
                    write(g1, "$(city_labels[city])\n")
                    write(g2, "$(country_labels[country])\n")
                end
            end
        end
    end

    for (labels, name) in [(city_labels, "city-names"),
                           (country_labels, "country-names")]
        info = [(v, k) for (k, v) in labels]
        sort!(info)
        open("$name.txt", "w") do f
            for (v, k) in info
                write(f, "$v $k\n")
            end
        end
    end
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
using LinearAlgebra
using StaticArrays
using Distributions
using PyPlot
const Vec2f = SVector{2, Float64}
const Vec4f = SVector{4, Float64}
const VecFrame = SVector{18, Float64}

include("src/geometry/geometry.jl")
include("src/world.jl")
include("src/double_integrator.jl")

@inline function distance(p::Vec2f, q::Vec2f)
    return sqrt((p[1]-q[1])^2+(p[2]-q[2])^2)
end

function read_scenario(file_name)
    f_r = open(file_name, "r")
    n = countlines(f_r)
    seekstart( f_r )

    scenario = VecFrame[]
    for i = 1:n
        line = split( readline( f_r ), "," )
        frame = VecFrame(parse(Float64, line[1] ),
                    parse(Float64, line[2] ),
                    parse(Float64, line[3] ),
                    parse(Float64, line[4] ),
                    parse(Float64, line[5] ),
                    parse(Float64, line[6] ),
                    parse(Float64, line[7] ),
                    parse(Float64, line[8] ),
                    parse(Float64, line[9] ),
                    parse(Float64, line[10] ),
                    parse(Float64, line[11] ),
                    parse(Float64, line[12] ),
                    parse(Float64, line[13] ),
                    parse(Float64, line[14] ),
                    parse(Float64, line[15] ),
                    parse(Float64, line[16] ),
                    parse(Float64, line[17] ),
                    parse(Float64, line[18] ))
        push!(scenario, frame)
    end
    start_frame=scenario[1][2]
    print(start_frame)
    return scenario, start_frame
end

function read_sample_points(file_name)
    f_r = open(file_name, "r")
    n = countlines(f_r)
    seekstart( f_r )

    sample_points = Vec4f[]
    for i = 1:n
        line = split( readline( f_r ), " " )
        p = Vec4f(parse(Float64, line[1] ),
                    parse(Float64, line[2] ),
                    parse(Float64, line[3] ),
                    parse(Float64, line[4] ))
        push!(sample_points, p)
    end
    return sample_points
end

function get_frame(scenario::Vector{VecFrame}, frame_num::Int64, start_frame)
    # start_frame = 1530
    choose_frame = start_frame + frame_num
    frame = VecFrame[]
    n = length(scenario)
    for i = 1:n
        if convert(Int64, scenario[i][2]) == choose_frame
            push!(frame, scenario[i])
        end
    end
    return frame
end

function find_neighbors(frame::Vector{VecFrame}, center::Vec2f, lane_id::Int64, r::Float64)
    n = length(frame)
    neighbors = VecFrame[]
    front_neighbors = VecFrame[]
    for i = 1:n
        v = Vec2f(frame[i][5], frame[i][6])
        if distance(v, center) < r
            push!(neighbors, frame[i])
            if v[2] >= center[2]
                push!(front_neighbors, frame[i])
            end
        end
    end
    return neighbors, front_neighbors
end

function set_obstacles(neighbors::Vector{VecFrame})
    n = length(neighbors)
    obstacle_set = Rectangle[]
    for i = 1:n
        obs = Rectangle([neighbors[i][5], neighbors[i][6]], neighbors[i][10], neighbors[i][9])
        push!(obstacle_set, obs)
    end
    return obstacle_set
end

function get_boundary(neighbors::Vector{VecFrame}, ego_vehicle::Vec2f, lane_id::Int64, s_goal::Vec4f)
    lane1_c = 7.231027419583393
    lane2_c = 18.84918099161813
    lane3_c = 29.990037979406054
    lane4_c = 40.99582664476257
    lane5_c = 53.099915981198585
    lane12 = (lane1_c + lane2_c)/2
    lane23 = (lane2_c + lane3_c)/2
    lane34 = (lane3_c + lane4_c)/2
    lane45 = (lane4_c + lane5_c)/2
    lane5edge = lane5_c + (lane5_c - lane45)
    center = ego_vehicle
    n = length(neighbors)
    y_max = neighbors[1][6]
    y_min = neighbors[1][6]
    for i = 1:n
        if neighbors[i][6] > y_max
            y_max = neighbors[i][6]
        end
        if neighbors[i][6] < y_min
            y_min = neighbors[i][6]
        end
    end
    if y_min > center[2]
        y_min = center[2]
    end
    if y_max < s_goal[2]
        y_max = s_goal[2] + 10
    end
    xmin = [0.0, 0.0]
    xmax = [0.0, 0.0]
    # println("ego vehicle:", center)
    if lane_id == 1
        global xmin = [0, y_min-10]
        global xmax = [lane23, y_max]
    elseif lane_id == 2
        global xmin = [0, y_min-10]
        global xmax = [lane34, y_max]
    elseif lane_id == 3
        global xmin = [lane12, y_min-10]
        global xmax = [lane45, y_max]
    elseif lane_id == 4
        global xmin = [lane23, y_min-10]
        global xmax = [lane5edge, y_max]
    elseif lane_id == 5
        global xmin = [lane34, y_min-10]
        global xmax = [lane5edge, y_max]
    end
    vmin = [-20.0, 0.0]
    vmax = [20.0, 100.0]
    return xmin, xmax, vmin, vmax
end

function nearest_front_vehicle(neighbors::Vector{VecFrame}, host::Vec2f, lane_id::Int64)
    min_d = 1500
    nearest_v = neighbors[1]
    for v in neighbors
        if v[14] == lane_id && v[6] > host[2]
            if distance(Vec2f(v[5], v[6]), host) < min_d
                min_d = distance(Vec2f(v[5], v[6]), host)
                nearest_v = v
            end
        end
    end
    return min_d, nearest_v[12]
end

function nearest_front_left_vehicle(neighbors::Vector{VecFrame}, host::Vec2f, lane_id::Int64)
    min_d = 1500
    nearest_v = neighbors[1]
    for v in neighbors
        if v[14] == lane_id-1 && v[6] > host[2]
            if 100 < distance(Vec2f(v[5], v[6]), host) < min_d
                min_d = distance(Vec2f(v[5], v[6]), host)
                nearest_v = v
                y_d = v[6] - host[2]
            end
        end
    end
    return y_d, nearest_v[12]
end

function nearest_front_right_vehicle(neighbors::Vector{VecFrame}, host::Vec2f, lane_id::Int64)
    min_d = 1500
    nearest_v = neighbors[1]
    y_d = 0.0
    for v in neighbors
        if v[14] == lane_id+1 && v[6] > host[2]
            if 100 < distance(Vec2f(v[5], v[6]), host) < min_d
                min_d = distance(Vec2f(v[5], v[6]), host)
                nearest_v = v
                y_d = v[6] - host[2]
            end
        end
    end
    return y_d, nearest_v[12]
end

function nearest_back_vehicle(neighbors::Vector{VecFrame}, host::Vec2f, lane_id::Int64)
    min_d = 1500
    nearest_v = neighbors[1]
    y_d = 0.0
    for v in neighbors
        if v[14] == lane_id && v[6] < host[2]
            if distance(Vec2f(v[5], v[6]), host) < min_d
                min_d = distance(Vec2f(v[5], v[6]), host)
                nearest_v = v
            end
        end
    end
    return min_d, nearest_v[12]
end

function get_init_goal(file_name)
    f_host = open(file_name, "r")
    n = countlines(f_host)
    seekstart( f_host )

    host_scenario = VecFrame[]
    for i = 1:n
        line = split( readline( f_host ), "," )
        frame = VecFrame(parse(Float64, line[1] ),
                    parse(Float64, line[2] ),
                    parse(Float64, line[3] ),
                    parse(Float64, line[4] ),
                    parse(Float64, line[5] ),
                    parse(Float64, line[6] ),
                    parse(Float64, line[7] ),
                    parse(Float64, line[8] ),
                    parse(Float64, line[9] ),
                    parse(Float64, line[10] ),
                    parse(Float64, line[11] ),
                    parse(Float64, line[12] ),
                    parse(Float64, line[13] ),
                    parse(Float64, line[14] ),
                    parse(Float64, line[15] ),
                    parse(Float64, line[16] ),
                    parse(Float64, line[17] ),
                    parse(Float64, line[18] ))
        push!(host_scenario, frame)
    end
    init_goal = Vec4f[]
    for i = 1:(n-1)
        cur = host_scenario[i]
        next = host_scenario[i+1]
        theta = atan((next[6]-cur[6])/(next[5]-cur[5]))
        if theta < 0
            theta = pi + theta
        end
        v_x = cur[12] * cos(theta)
        v_y = cur[12] * sin(theta)
        p = Vec4f(cur[5], cur[6], v_x, v_y)
        if i%10 == 1
            push!(init_goal, p)
        end
    end
    cur = host_scenario[n-1]
    next = host_scenario[n]
    theta = atan((next[6]-cur[6])/(next[5]-cur[5]))
    if theta < 0
        theta = pi + theta
    end
    v_x = next[12] * cos(theta)
    v_y = next[12] * sin(theta)
    p = Vec4f(next[5], next[6], v_x, v_y)
    push!(init_goal, p)

    return init_goal
end


text:
"""
    theta(z, R, char=[], derivs=[], derivs_t=[])

Compute the Riemann theta function of the matrix in R at z, with optional characteristics or derivatives.

The input derivatives must be unit vectors, and the input z must be of the form ``z = a + τb, b∈[0,1)^g, a∈ \\mathbb{R}^g``.

# Arguments
- `z::Array{<:Number}`: Array of size g, each entry is a complex number.
- `R::RiemannMatrix`: RiemannMatrix type containing matrix.
- `char::Array{}=[]`: 2-element array containing 2 arrays of size g, where entries are either 0 or 1.
- `derivs::Array{}=[]`: N-element array specifying the derivative to the N-th order with respect to z, where each element is an array representing a unit vector in the direction of the derivative.
- `derivs_t::Array{}=[]`: M-element array specifying the derivative to the M-th order with respect to τ, where each element is an array of size 2 whose entries are the index in τ.

# Examples
```julia
julia> theta([0,0], RiemannMatrix([1+im -1; -1 1+im]), char=[[1,0],[1,1]], derivs=[[1,0]], derivs_t=[[2,1]])
```
"""
function theta(z::Array{<:Number}, R::RiemannMatrix; char::Array{}=[], derivs::Array{}=[], derivs_t::Array{}=[])
    # compute theta function with characteristics
    shift_x = zeros(R.g); # shift z by characteristic
    shift_n = zeros(R.g); # shift points in ellipsoid by characteristic
    if length(char) > 0
        shift_x = 0.5*char[2];
        shift_n = 0.5*char[1];
    end
    # compute derivatives with respect to τ by converting to derivatives with respect to z using the heat equation
    deriv_t_factor = 1; # multiply result by this factor
    for d in derivs_t
        if d[1] == d[2]
            deriv_t_factor *= 1/(4*π*im);
        else
            deriv_t_factor *= 1/(2*π*im);
        end
        deriv = [zeros(R.g), zeros(R.g)];
        deriv[1][d[1]] = 1;
        deriv[2][d[2]] = 1;
        derivs = vcat(derivs, deriv);
    end
    # compute theta function
    x = real(z) + shift_x; # shift z by characteristic
    y = convert(Array{Float64}, imag(z));
    y0 = (R.g > 1 ?  inv(R.Y)*y : inv.(R.Y).*y); # Y⁻¹y
    exp_part = exp((π*transpose(y)*y0)[1]); # exponential part of theta function
    osc_part = oscillatory_part(R, x, y0, shift_n, derivs); # oscillatory part of theta function
    return deriv_t_factor*exp_part*osc_part;
end

"""
    oscillatory_part(R, x, y0, shift_n, derivs=[])

Compute the oscillatory part of the Riemann theta function. Helper function to [`theta`](@ref).
"""
function oscillatory_part(R::RiemannMatrix, x::Array{<:Real}, y0::Array{<:Real}, shift_n::Array{<:Real}, derivs::Array{}=[])
    nderivs = size(derivs)[1];
    E = [n + shift_n for n in R.ellipsoid[nderivs + 1]]; # shift points in ellipsoid for computing theta function with characteristics
    w = round.(y0);
    w0 = y0 - w;
    s = (2*π*im)^(nderivs) * sum([
(nderivs > 0 ? prod(transpose(d)*(p-w) for d in derivs) : 1) * exp(2*π*im*((0.5*(transpose(p-w)*R.X*(p-w))[1] + transpose(p-w)*x)[1]) - π*(transpose(p+w0)*R.Y*(p+w0))[1]) for p in E]); # compute sum over points in ellipsoid
    return s;
end


text:
push!(LOAD_PATH, ".")

module GAMutation

using Debug
using Individual
using BinaryIndividual
using RealIndividual
using Parallel

abstract type BitFlipMutation end
abstract type UniformMutation end

function mutate!(ind::BinaryIndividual.AbstractBinaryIndividual, ::Type{BitFlipMutation},  mr::Float64)

  Debug.ga_debug && println("----- Bit Flip -----\n")
  Debug.ga_debug && println("Before mutation: ", ind[:])

  for i in eachindex(ind)
    pr = Parallel.threadRand()
    if pr < mr
      ind[i] = xor(ind[i], 1)
    end
  end

  Debug.ga_debug && println("After mutation: ", ind[:], "\n")
  Debug.ga_debug && println("----- Bit Flip End -----\n")

  return ind
end

function mutate!(ind::RealIndividual.AbstractRealIndividual, ::Type{UniformMutation}, mr::Float64, lb::Float64, ub::Float64)

  Debug.ga_debug && println("----- Uniform -----\n")
  Debug.ga_debug && println("Before mutation: ", ind[:])

  for i in eachindex(ind)
    pr = Parallel.threadRand()
    if pr < mr
      Δ = lb + Parallel.threadRand() * (ub - lb)
      ind[i] += Parallel.threadRand([1, -1]) * Δ
      ind[i] = max(lb, min(ind[i], ub))
    end
  end

  Debug.ga_debug && println("After mutation: ", ind[:], "\n")
  Debug.ga_debug && println("----- Uniform End -----\n")

  return ind
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
struct Word{T}
    alphabet::Vector{T}
    letters::Vector{Int}

    function Word(a::AbstractVector{T}, l::AbstractVector{<:Integer}) where {T}
        all(i -> 1 <= i <= length(a), l) ||
            throw(ArgumentError("Invalid word over alphabet $a: $w"))
        return new{T}(a, l)
    end
end
Base.length(w::Word) = length(w.letters)
function Base.show(io::IO, w::Word)
    if isone(w)
        print(io, "(id)")
    else
        join(io, w.alphabet[w.letters], "·")
    end
end

Base.:(==)(w::Word, v::Word) = w.alphabet == v.alphabet && w.letters == v.letters
Base.hash(w::Word, h::UInt) = hash(w.alphabet, hash(w.letters, hash(Word, h)))

Base.one(w::Word) = Word(w.alphabet, Int[])
Base.isone(w::Word) = length(w) == 0

function Base.:*(w::Word, z::Word)
    @assert w.alphabet == z.alphabet
    return Word(w.alphabet, [w.letters; z.letters])
end

function StarAlgebras.star(w::Word)
    newletters = similar(w.letters)
    A = w.alphabet

    # star(:a) = :b
    # star(:b) = :a
    # star(:c) = :c

    for (i, l) in enumerate(Iterators.reverse(w.letters))
        k = if A[l] === :a
            2
        elseif A[l] === :b
            1
        else
            l
        end
        newletters[i] = k
    end
    return Word(w.alphabet, newletters)
end

function words(alphabet; radius)

    words = [Word(alphabet, Int[])] # word identity

    for r in 1:radius
        append!(
            words,
            [
                Word(alphabet, collect(w)) for
                w in Iterators.product(fill(1:length(alphabet), r)...)
            ],
        )
    end

    return words
end


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

# Load necessary components of Espresso
using Espresso: subs, simplify, ExGraph, ExH, to_expr, sanitize, genname,
    find_vars, findex, find_indices, isindexed

# Define some constants to create the new (parsed) functions
# The (irrelevant) `nothing` below is there to have a :block Expr; deleted later
const _HEAD_PARSEDFN_SCALAR = sanitize(:(
function jetcoeffs!(::Val{__fn}, __tT::Taylor1{_T}, __x::Taylor1{_S}, __params) where
        {_T<:Real, _S<:Number}

    order = __tT.order
    nothing
end)
);

const _HEAD_PARSEDFN_VECTOR = sanitize(:(
function jetcoeffs!( ::Val{__fn}, __tT::Taylor1{_T}, __x::AbstractVector{Taylor1{_S}},
        __dx::AbstractVector{Taylor1{_S}}, __params) where {_T<:Real, _S<:Number}

    order = __tT.order
    nothing
end)
);

# Constants for the initial declaration and initialization of arrays
const _DECL_ARRAY = sanitize( Expr(:block,
    :(__var1 = Array{Taylor1{_S}}(undef, __var2)),
    :(__var1 .= Taylor1( zero(_S), order ))
    ) )

"""
`_make_parsed_jetcoeffs( ex, debug=false )`

This function constructs a new function, equivalent to the
differential equations, which exploits the mutating functions
of TaylorSeries.jl.

"""
function _make_parsed_jetcoeffs(ex::Expr, debug=false)

    # Extract the name, args and body of the function
    fn, fnargs, fnbody = _extract_parts(ex)
    debug && (println("****** _extract_parts ******");
        @show(fn, fnargs, fnbody); println())

    # Set up new function
    newfunction = _newhead(fn, fnargs)

    # for-loop block (recursion); it will contain the parsed body
    forloopblock = Expr(:for, :(ord = 1:order-1),
        Expr(:block, :(ordnext = ord + 1)) )

    #= Transform graph representation of the body of the function
    - `defspreamble` includes the definitions used for zeroth order (preamble)
    - `fnbody` is the transformed function body, using mutating functions
        from TaylorSeries, and is used within the for-loop block
    - `retvar` is the (guessed) return variable, which defines the LHS
        of the ODEs
    =#
    defspreamble, fnbody, retvar = _preamble_body(fnbody, fnargs, debug)
    debug && (println("****** _preamble_body ******");
        @show(defspreamble); println(); @show(fnbody); println();)

    #= Create body of recursion loop; temporary assignements may be needed.
    - rec_preamb
    - rec_fnbody
    =#
    debug && println("****** _recursionloop ******")
    rec_preamb, rec_fnbody = _recursionloop(fnargs, retvar)

    # Add preamble to newfunction
    push!(newfunction.args[2].args, defspreamble..., rec_preamb)

    # Add parsed fnbody to forloopblock
    push!(forloopblock.args[2].args, fnbody.args[1].args..., rec_fnbody)

    # Push preamble and forloopblock to newfunction
    push!(newfunction.args[2].args, forloopblock, Meta.parse("return nothing"))

    # Rename variables of the body of the new function
    if length(fnargs) == 3
        newfunction = subs(newfunction, Dict(:__tT => fnargs[3],
            :__params => fnargs[2], :(__x) => fnargs[1], :(__fn) => fn ))
    elseif length(fnargs) == 4
        newfunction = subs(newfunction, Dict(:__tT => fnargs[4],
            :__params => fnargs[3], :(__x) => fnargs[2], :(__dx) => fnargs[1],
            :(__fn) => fn ))
    end

    return newfunction
end

"""
`_extract_parts(ex::Expr)`

Returns the function name, the function arguments, and the body of
a function passed as an `Expr`. The function may be provided as
a one-line function, or in the long form (anonymous functions
do not work).

"""
function _extract_parts(ex::Expr)

    # Capture name, args and body
    # @capture( shortdef(ex), ffn_(ffnargs__) = fnbody_ ) ||
    #     throw(ArgumentError("It must be a function call:\n $ex"))
    dd = Dict{Symbol, Any}()
    _capture_fn_args_body!(ex, dd)
    ffn, ffnargs, fnbody = dd[:fname], dd[:fnargs], dd[:fnbody]

    # Clean-up `ffn`
    fn = isa(ffn, Symbol) ? ffn : ffn.args[1]

    # Clean-up (from declarations) `ffnargs`
    fnargs = []
    for ff in ffnargs
        if isa(ff, Symbol)
            push!(fnargs, ff)
        elseif isa(ff, Expr)
            push!(fnargs, ff.args[1])
        else
            throw(ArgumentError("$ff is neither a `Symbol` nor an `Expr`"))
        end
    end

    # Standarize fnbody, same structure for one-line or long-form functions
    if length(fnbody.args) > 1
        fnbody.args[1] = Expr(:block, copy(fnbody.args)...)
        deleteat!(fnbody.args, 2:length(fnbody.args))
    end

    # Special case: the last arg of the function is a simple
    # assignement (symbol), a numerical value or simply `nothing`
    if isa(fnbody.args[1].args[end], Symbol)
        # This accouns the case of having `nothing` at the end of the function
        if fnbody.args[1].args[end] == :nothing
            pop!(fnbody.args[1].args)
        else
            fnbody.args[1].args[end] = :( identity($(fnbody.args[1].args[end])) )
        end
    elseif isa(fnbody.args[1].args[end], Number)
        fnbody.args[1].args[end] = :( $(fnbody.args[1].args[end])+zero($(fnargs[1])) )
    end

    return fn, fnargs, fnbody
end

"""
`_capture_fn_args_body!(ex, vout::Dict{Symbol, Any})`

Captures the name of a function, arguments, body and other properties,
returning them as the values of the dictionary `dd`, which is updated
in place.

"""
function _capture_fn_args_body!(ex::Expr, vout::Dict{Symbol, Any} = Dict())
    exh = ExH(ex)
    if exh.head == :(=) || exh.head == :function
        _capture_fn_args_body!.(ex.args, Ref(vout))
    elseif exh.head == :call
        push!(vout, Pair(:fname, exh.args[1]))
        push!(vout, Pair(:fnargs, exh.args[2:end]))
    elseif exh.head == :block
        push!(vout, Pair(:fnbody, ex))
    elseif exh.head == :where
        push!(vout, Pair(:where, ex.args[2:end]))
        _capture_fn_args_body!(ex.args[1], vout)
    end
    vout
end

"""
`_newhead(fn, fnargs)`

Creates the head of the new method of `jetcoeffs!`. Here,
`fn` is the name of the passed function and `fnargs` is a
vector with its arguments (which are two or three).

"""
function _newhead(fn, fnargs)

    # Construct common elements of the new expression
    if length(fnargs) == 3
        newfunction = copy(_HEAD_PARSEDFN_SCALAR)
    elseif length(fnargs) == 4
        newfunction = copy(_HEAD_PARSEDFN_VECTOR)
    else
        throw(ArgumentError(
        "Wrong number of arguments in the definition of the function $fn"))
    end

    # Add `TaylorIntegration` (module name) to create a new method of `jetcoeffs!`
    newfunction.args[1].args[1].args[1] =
        Expr(:., :TaylorIntegration, :(:jetcoeffs!))

    # Delete the irrelevant `nothing`
    pop!(newfunction.args[2].args)

    return newfunction
end

"""
`_preamble_body(fnbody, fnargs, debug=false)`

Returns the preamble, the body and a guessed return
variable, which will be used to build the parsed
function. `fnbody` is the expression with the body of
the original function, `fnargs` is a vector of symbols
of the original diferential equations function.

"""
function _preamble_body(fnbody, fnargs, debug=false)

    # Bookkeeping:
    #   v_vars:   List of symbols with created variables
    #   v_assign: Numeric assignments (to be deleted)
    v_vars = Union{Symbol,Expr}[]
    v_assign = Dict{Union{Symbol,Expr}, Number}()

    #= Rename vars to have the body in non-indexed form
    `d_indx` is a dictionary mapping new variables (symbols) to old
    (perhaps indexed) symbols
    =#
    fnbody, d_indx = _rename_indexedvars(fnbody)
    debug && (println("------ _rename_indexedvars ------");
        @show(fnbody); println(); @show(d_indx); println())

    #= Create newfnbody
    - v_newindx: symbols of auxiliary indexed vars
    - v_arraydecl: symbols which are explicitly declared as Array or Vector
    - `newfnbody` corresponds to `fnbody`, cleaned (without irrelevant comments)
        and with all new variables in place
    =#
    newfnbody, v_newindx, v_arraydecl = _newfnbody(fnbody, fnargs, d_indx)
    debug && (println("------ _newfnbody ------");
        @show(v_newindx); println(); @show(v_arraydecl); println();
        @show(newfnbody); println())

    # Parse `newfnbody` and create `prepreamble`, updating the bookkeeping vectors.
    # The returned `newfnbody` and `prepreamble` use the mutating functions
    # of TaylorSeries.
    prepreamble = Expr(:block,)
    _parse_newfnbody!(newfnbody, prepreamble, v_vars, v_assign, d_indx,
        v_newindx, v_arraydecl)
    preamble = prepreamble.args[1]

    # Substitute numeric assignments in new function's body
    newfnbody = subs(newfnbody, Dict(v_assign))
    preamble = subs(preamble, Dict(v_assign))
    debug && (println("------ _parse_newfnbody! ------");
        @show(newfnbody); println(); @show(preamble); println();
        @show(v_vars); println(); @show(d_indx); println();
        @show(v_newindx); println())

    # Include the assignement of indexed auxiliary variables
    defspreamble = _defs_preamble!(preamble, fnargs,
        d_indx, v_newindx, v_arraydecl, Union{Symbol,Expr}[], empty(d_indx))
    # preamble = subs(preamble, d_indx)

    # Bring back substitutions
    newfnbody = subs(newfnbody, d_indx)

    # Define retvar; for scalar eqs is the last entry included in v_vars
    retvar = length(fnargs) == 3 ? subs(v_vars[end], d_indx) : fnargs[1]

    debug && (println("------ _defs_preamble! ------");
        @show(d_indx); println(); @show(v_newindx); println();
        @show(newfnbody); println(); @show(retvar); println())

    return defspreamble, newfnbody, retvar
end

"""
`_rename_indexedvars(fnbody)`

Renames the indexed variables (using `Espresso.genname()`) that
exists in `fnbody`. Returns `fnbody` with the renamed variables
and a dictionary that links the new variables to the old
indexed ones.

"""
function _rename_indexedvars(fnbody)
    # Thanks to Andrei Zhabinski (@dfdx, see #31) for this
    # implementation, now slightly modified
    indexed_vars = findex(:(_X[_i...]), fnbody)
    st = Dict(ivar => genname() for ivar in indexed_vars)
    new_fnbody = subs(fnbody, st)
    return new_fnbody, Dict(v => k for (k, v) in st)
end

"""
`_newfnbody(fnbody, fnargs, d_indx)`

Returns a new (modified) body of the function, a priori unfolding
the expression graph (AST) as unary and binary calls, and a vector
(`v_newindx`) with the name of auxiliary indexed variables.

"""
function _newfnbody(fnbody, fnargs, d_indx)
    # `fnbody` is assumed to be a `:block` `Expr`
    newfnbody = Expr(:block,)
    v_newindx = Symbol[]
    v_arraydecl = Symbol[]

    # The magic happens HERE!!
    # Each line of fnbody (fnbody.args) is parsed separately
    for (i, ex) in enumerate(fnbody.args)
        if isa(ex, Expr)
            ex_head = ex.head

            # Ignore the following cases
            (ex_head == :return) && continue

            # Treat `for` loops and `if` blocks separately
            if ex_head == :block
                newblock, tmp_newindx, tmp_arraydecl = _newfnbody(ex, fnargs, d_indx)
                push!(newfnbody.args, newblock )
                append!(v_newindx, tmp_newindx)
                append!(v_arraydecl, tmp_arraydecl)
            elseif ex_head == :for
                push!(newfnbody.args, Expr(:for, ex.args[1]))
                loopbody, tmp_newindx, tmp_arraydecl = _newfnbody( ex.args[2], fnargs, d_indx )
                push!(newfnbody.args[end].args, loopbody)
                append!(v_newindx, tmp_newindx)
                append!(v_arraydecl, tmp_arraydecl)
            elseif ex_head == :if
                push!(newfnbody.args, Expr(:if, ex.args[1]))
                for exx in ex.args[2:end]
                    ifbody, tmp_newindx, tmp_arraydecl = _newfnbody( exx, fnargs, d_indx)
                    push!(newfnbody.args[end].args, ifbody)
                    append!(v_newindx, tmp_newindx)
                    append!(v_arraydecl, tmp_arraydecl)
                end
            elseif ex_head == :(=) || ex_head == :call  # assignements or function calls

                ex_lhs = ex.args[1]
                ex_rhs = ex.args[2]

                # Case of explicit declaration of Array or Vector
                # TO-DO: Include cases where the definition uses `similar`
                if !isempty(findex(:(_Array{_TT...}), ex)) ||
                        !isempty(findex(:(Vector{_TT...}), ex))
                    push!(v_arraydecl, ex_lhs)
                    push!(newfnbody.args, ex)
                    continue
                end

                # Case of multiple assignments, e.g., ex = :((σ, ρ, β) = par)
                # Distinguishes if the RHS is the parameter or the dependent vars
                if isa(ex_lhs, Expr)
                    if ex_lhs.head == :tuple && length(fnargs) == 4 && isa(ex_rhs, Symbol)
                        # Construct new expressions
                        if ex_rhs == fnargs[3]  # params
                            for i in eachindex(ex_lhs.args)
                                ex_indx = Expr(:local, Expr(:(=), ex_lhs.args[i], Expr(:ref, ex_rhs, i)))
                                push!(newfnbody.args, ex_indx)
                            end
                        elseif ex_rhs == fnargs[2] # `x`
                            for i in eachindex(ex_lhs.args)
                                ex_indx = Expr(:(=), ex_lhs.args[i], Expr(:ref, ex_rhs, i))
                                push!(newfnbody.args, ex_indx)
                            end
                        else
                            throw(ArgumentError("Assignment no implemented in $(ex)"))
                        end
                        #
                    end
                    continue
                end

                # Unfold AST graph
                nex = deepcopy(ex)
                try
                    nex = to_expr(ExGraph(simplify(ex)))
                catch
                    # copy `ex` as it is, if it is not "recognized"
                    push!(newfnbody.args, ex)
                    continue
                end
                push!(newfnbody.args, nex.args[2:end]...)

                # Bookkeeping of indexed vars, to define assignements
                isindx_lhs = haskey(d_indx, ex_lhs)
                for nexargs in nex.args[2:end]
                    # (nexargs.head == :line ||
                    #     haskey(d_indx, nexargs.args[1])) && continue
                    haskey(d_indx, nexargs.args[1]) && continue
                    vars_nex = find_vars(nexargs)

                    # any(haskey.(d_indx, vars_nex[:])) &&
                    #     !in(vars_nex[1], v_newindx) &&
                    (isindx_lhs || vars_nex[1] != ex_lhs) && push!(v_newindx, vars_nex[1])
                end
            elseif ex_head == :local
                # If declared as `local`, copy `ex` as it is. In some cases this
                # helps performance. Very useful for including (numeric) constants
                push!(newfnbody.args, ex)
                #
            else
                # If not implemented, throw an `ArgumentError`
                throw(ArgumentError("$(ex_head) is not yet implemented; $(typeof(ex))"))
            end
            #
        elseif isa(ex, LineNumberNode)
            continue  # Ignore `LineNumberNode`s
        else
            # In any other case, throw an `ArgumentError`
            throw(ArgumentError("$ex is not an `Expr`; $(typeof(ex))"))
            #
        end
    end

    # Needed, if `newfnbody` consists of a single assignment (unary call)
    if newfnbody.head == :(=)
        newfnbody = Expr(:block, newfnbody)
    end

    return newfnbody, v_newindx, v_arraydecl
end

"""
`_parse_newfnbody!(ex, preex, v_vars, v_assign, d_indx, v_newindx, v_arraydecl,
    [inloop=false])`

Parses `ex` (the new body of the function) replacing the expressions
to use the mutating functions of TaylorSeries, and building the preamble
`preex`. This is done by traversing recursively the args of `ex`, updating
the bookkeeping vectors `v_vars` and `v_assign`. `d_indx` is
used to substitute back the explicit indexed variables.

"""
function _parse_newfnbody!(ex::Expr, preex::Expr,
        v_vars, v_assign, d_indx, v_newindx, v_arraydecl, inloop=false)

    # Numeric assignements to be deleted
    indx_rm = Int[]

    # Magic happens HERE
    # Each line of ex (ex.args) is parsed separately
    for (i, aa) in enumerate(ex.args)

        # Treat for loops, blocks and if blocks separately
        if (aa.head == :for)
            push!(preex.args, Expr(:for, aa.args[1]))
            _parse_newfnbody!(aa, preex.args[end],
                v_vars, v_assign, d_indx, v_newindx, v_arraydecl, true)
            #
        elseif (aa.head == :block)

            push!(preex.args, Expr(:block))
            _parse_newfnbody!(aa, preex.args[end],
                v_vars, v_assign, d_indx, v_newindx, v_arraydecl)
            #
        elseif (aa.head == :if)
            push!(preex.args, Expr(:if, aa.args[1]))
            for exx in aa.args[2:end]
                push!(preex.args[end].args, Expr(:block))
                _parse_newfnbody!(exx, preex.args[end].args[end],
                    v_vars, v_assign, d_indx, v_newindx, v_arraydecl)
            end
            #
        elseif (aa.head == :(=))

            i == 1 && inloop && continue

            # Main case
            aa_lhs = aa.args[1]
            aa_rhs = aa.args[2]

            # Replace expressions when needed, and bookkeeping
            if isa(aa_rhs, Expr)
                (aa_rhs.args[1] == :eachindex || aa_rhs.head == :(:)) && continue
                # Replace expression
                _replace_expr!(ex, preex, i, aa_lhs, aa_rhs, v_vars, d_indx, v_newindx)

                # Remove new aa_lhs declaration from `ex`; it is already declared
                in(aa_lhs, v_arraydecl) && push!(indx_rm, i)
                #
            elseif isa(aa_rhs, Symbol) # occurs when there is a simple assignment

                bb = subs(aa, Dict(aa_rhs => :(identity($aa_rhs))))
                bb_lhs = bb.args[1]
                bb_rhs = bb.args[2]
                # Replace expression
                _replace_expr!(ex, preex, i, bb_lhs, bb_rhs, v_vars, d_indx, v_newindx)
                #
            elseif isa(aa_rhs, Number)  # case of numeric values
                push!(v_assign, aa_lhs => aa_rhs)
                push!(indx_rm, i)
                #
            else # needed?
                error("Either $aa or $typeof(aa_rhs[2]) are different from `Expr`, `Symbol` or `Number`")
                #
            end
            #
        elseif aa.head == :local
            # If declared as `local`, copy `ex` as it is, and delete it from
            # the recursion body
            push!(preex.args, aa)
            push!(indx_rm, i)   # delete associated expr in body function
            #
        else

            # Unrecognized head; pass the expression as it is
            push!(preex.args, aa)
        end
    end

    # Delete trivial numeric assignement statements
    isempty(indx_rm) || deleteat!(ex.args, indx_rm)

    return nothing
end

"""
`_replace_expr!(ex, preex, i, aalhs, aarhs, v_vars, d_indx, v_newindx)`

Replaces the calls in `ex.args[i]`, and updates `preex` with the definitions
of the expressions, based on the the LHS (`aalhs`) and RHS (`aarhs`) of the
base assignment. The bookkeeping vectors (`v_vars`, `v_newindx`)
are updated. `d_indx` is used to bring back the indexed variables.

"""
function _replace_expr!(ex::Expr, preex::Expr, i::Int,
        aalhs, aarhs, v_vars, d_indx, v_newindx)

    # Replace calls
    fnexpr, def_fnexpr, auxfnexpr = _replacecalls!(aarhs, aalhs, v_vars)

    # Bring back indexed variables
    fnexpr = subs(fnexpr, d_indx)
    def_fnexpr = subs(def_fnexpr, d_indx)
    auxfnexpr = subs(auxfnexpr, d_indx)

    # Update `ex` and `preex`
    push!(preex.args, def_fnexpr)
    ex.args[i] = fnexpr

    # Same for the auxiliary expressions
    if auxfnexpr.head != :nothing
        push!(preex.args, auxfnexpr)
        # in(aalhs, v_newindx) && push!(v_newindx, auxfnexpr.args[1])
        push!(v_newindx, auxfnexpr.args[1])
    end

    !in(aalhs, v_vars) && push!(v_vars, subs(aalhs, d_indx))
    return nothing
end

"""
`_replacecalls!(fnold, newvar, v_vars)`

Replaces the symbols of unary and binary calls
of the expression `fnold`, which defines `newvar`,
by the mutating functions in TaylorSeries.jl.
The vector `v_vars` is updated if new auxiliary
variables are introduced.

"""
function _replacecalls!(fnold::Expr, newvar::Symbol, v_vars)

    ll = length(fnold.args)
    dcall = fnold.args[1]
    newarg1 = fnold.args[2]

    # If call is not in mutating functions dictionaries, copy original one
    if !( in(dcall, keys(TaylorSeries._dict_unary_calls)) ||
            in(dcall, keys(TaylorSeries._dict_binary_calls)) )

        fnexpr = :($newvar = $fnold)
        def_fnexpr = fnexpr
        aux_fnexpr = Expr(:nothing)
        return fnexpr, def_fnexpr, aux_fnexpr
    end

    if ll == 2

        # Unary call
        # Replacements
        fnexpr, def_fnexpr, aux_fnexpr = TaylorSeries._dict_unary_calls[dcall]
        fnexpr = subs(fnexpr, Dict(:_res => newvar, :_arg1 => newarg1, :_k => :ord))
        def_fnexpr = :( _res = Taylor1($(def_fnexpr.args[2]), order) )
        def_fnexpr = subs(def_fnexpr,
            Dict(:_res => newvar, :_arg1 => :(constant_term($(newarg1))),
                :_k => :ord))

        # Auxiliary expression
        if aux_fnexpr.head != :nothing
            newaux = genname()
            aux_fnexpr = :( _res = Taylor1($(aux_fnexpr.args[2]), order) )
            aux_fnexpr = subs(aux_fnexpr,
                Dict(:_res => newaux, :_arg1 => :(constant_term($(newarg1))),
                    :_aux => newaux))
            fnexpr = subs(fnexpr, Dict(:_aux => newaux))
            !in(newaux, v_vars) && push!(v_vars, newaux)
        end
    elseif ll == 3

        # Binary call; no auxiliary expressions needed
        newarg2 = fnold.args[3]

        # Replacements
        fnexpr, def_fnexpr, aux_fnexpr = TaylorSeries._dict_binary_calls[dcall]
        fnexpr = subs(fnexpr,
            Dict(:_res => newvar, :_arg1 => newarg1, :_arg2 => newarg2,
                :_k => :ord))

        def_fnexpr = :(_res = Taylor1($(def_fnexpr.args[2]), order) )
        def_fnexpr = subs(def_fnexpr,
            Dict(:_res => newvar,
                :_arg1 => :(constant_term($(newarg1))),
                :_arg2 => :(constant_term($(newarg2))),
                :_k => :ord))
    else
        # Recognized call, but not as a unary or binary call
        fnexpr = :($newvar = $fnold)
        def_fnexpr = fnexpr
        aux_fnexpr = Expr(:nothing)
    end

    return fnexpr, def_fnexpr, aux_fnexpr
end

"""
`_defs_preamble!(preamble, fnargs, d_indx, v_newindx, v_arraydecl, v_preamb,
    d_decl, [inloop=false])`

Returns a vector with expressions defining the auxiliary variables
in the preamble; it may modify `d_indx` if new variables are introduced.
`v_preamb` is for bookkeeping the introduced variables.

"""
function _defs_preamble!(preamble::Expr, fnargs,
        d_indx, v_newindx, v_arraydecl, v_preamb, d_decl, inloop::Bool=false,
        ex_aux::Expr = Expr(:block,))

    # Initializations
    defspreamble = Expr[]

    for (i, ex) in enumerate(preamble.args)

        if isa(ex, Expr)

            # Treat :block and :for separately
            if (ex.head == :block)
                newdefspr = _defs_preamble!(ex, fnargs, d_indx,
                    v_newindx, v_arraydecl, v_preamb, d_decl, inloop, ex_aux)
                append!(defspreamble, newdefspr)
            elseif (ex.head == :for)
                push!(ex_aux.args, ex.args[1])
                newdefspr = _defs_preamble!(ex.args[2], fnargs, d_indx,
                    v_newindx, v_arraydecl, v_preamb, d_decl, true, ex_aux)
                append!(defspreamble, newdefspr)
                pop!(ex_aux.args)
            elseif (ex.head == :if)
                for exx in ex.args[2:end]
                    newdefspr = _defs_preamble!(exx, fnargs, d_indx,
                        v_newindx, v_arraydecl, v_preamb, d_decl, inloop, ex_aux)
                    append!(defspreamble, newdefspr)
                end
                continue
            elseif (ex.head == :(=))

                # `ex.head` is a :(=) of some kind
                alhs = ex.args[1]
                arhs = subs(ex.args[2], d_indx) # substitute updated vars in rhs

                # Outside of a loop
                if !inloop
                    in(alhs, v_preamb) && continue
                    push!(defspreamble, subs(ex, d_decl))
                    push!(v_preamb, alhs)
                    continue
                end

                # Inside a loop
                if isindexed(alhs)

                    # `var1` may be a vector or a matrix; declaring it is subtle
                    var1 = alhs.args[1]
                    (in(var1, fnargs) || in(var1, v_arraydecl) ||
                        in(alhs, v_preamb)) && continue

                    # indices of var1
                    d_subs = Dict(veexx.args[1] => veexx.args[2]
                        for veexx in ex_aux.args)
                    indx1 = alhs.args[2:end]
                    ex_tuple = :( [$(indx1...)] )
                    ex_tuple = subs(ex_tuple, d_subs)
                    ex_tuple = :( size( $(Expr(:tuple, ex_tuple.args...)) ))
                    exx = subs(_DECL_ARRAY, Dict(:__var1 => :($var1),
                        :__var2 => :($ex_tuple)) )
                    push!(defspreamble, exx.args...)
                    push!(v_preamb, var1)
                    continue
                    #
                elseif in(alhs, v_newindx) || isindexed(arhs)
                    # `alhs` is an aux indexed var, so something in `arhs`
                    # is indexed.

                    in(alhs, v_preamb) && continue

                    vars_indexed = findex(:(_X[_i...]), arhs)

                    # NOTE: Uses the size of the var with more indices
                    # to define the declaration of the new array.
                    iimax, ii_indx = findmax(
                        [length(find_indices(aa)[1]) for aa in vars_indexed] )
                    var1 = vars_indexed[ii_indx].args[1]
                    indx1 = vars_indexed[ii_indx].args[2:end]

                    exx_indx = ones(Int, length(indx1))
                    push!(d_indx, alhs => :($alhs[$(indx1...)]) )
                    push!(d_decl, alhs => :($alhs[$exx_indx...]))
                    exx = subs(_DECL_ARRAY,
                        Dict(:__var1 => :($alhs), :__var2 => :(size($var1))) )
                    push!(defspreamble, exx.args...)
                    push!(v_preamb, alhs)
                    continue
                    #
                else

                    # `alhs` is not indexed nor an aux indexed var
                    in(alhs, v_preamb) && continue
                    vars_indexed = findex(:(_X[_i...]), arhs)
                    if !isempty(vars_indexed)
                        ex = subs(ex, Dict(vv => :(one(S)) for vv in vars_indexed))
                    end
                    ex = subs(ex, d_decl)
                    push!(defspreamble, ex)
                    push!(v_preamb, alhs)
                    continue
                    #
                end
                #
            end
        else
            throw(ArgumentError("$ex is not an `Expr`; $(typeof(ex))"))
            #
        end

        exx = subs(ex, d_indx)
        !inloop && push!(defspreamble, exx)

    end
    return defspreamble
end

"""
`_recursionloop(fnargs, retvar)`

Build the expression for the recursion-loop.

"""
function _recursionloop(fnargs, retvar)

    ll = length(fnargs)
    if ll == 3

        rec_preamb = sanitize( :( $(fnargs[2])[1] = $(retvar)[0] ) )
        rec_fnbody = sanitize( :( $(fnargs[2])[ordnext] = $(retvar)[ord]/ordnext ) )
        #
    elseif ll == 4

        retvar = fnargs[1]
        rec_preamb = sanitize(:(
            for __idx in eachindex($(fnargs[2]))
                $(fnargs[2])[__idx].coeffs[2] = $(retvar)[__idx].coeffs[1]
            end
        ))
        rec_fnbody = sanitize(:(
            for __idx in eachindex($(fnargs[2]))
                $(fnargs[2])[__idx].coeffs[ordnext+1] =
                    $(retvar)[__idx].coeffs[ordnext]/ordnext
            end
        ))
        #
    else

        throw(ArgumentError(
        "Wrong number of arguments in the definition of the function $fn"))
    end
    return rec_preamb, rec_fnbody
end

"""
`@taylorize expr`

This macro `eval`s the function given by `expr` and defines a new
method of [`jetcoeffs!`](@ref) which is specialized on that
function. Integrating via [`taylorinteg`](@ref) of
[`lyap_taylorinteg`](@ref) after using the macro yields better performance.

See the [documentation](@ref taylorize) for more details and limitations.

!!! warning
    This macro is on an experimental stage; check the integration
    results carefully.

"""
macro taylorize( ex )
    nex = _make_parsed_jetcoeffs(ex)
    quote
        eval( $(esc(ex)) )  # evals to calling scope the passed function
        eval( $(esc(nex)) ) # evals the new method of `jetcoeffs!`
    end
end


text:
using Plots
using Colors
using Statistics
using MAT
import Random

include("lambs_krylov.jl")

δt = step_size_lamb(load_lambs(mesh=2), 1.0);

# ---------------------------------------------
# Generar figuras condicion inicial singleton
# ---------------------------------------------
generar_figura_malla2_singleton_position()
generar_figura_malla2_singleton_velocity()

# ---------------------------------------------
# Generar datos envolvente
# ---------------------------------------------

vecsamples = [1, 10, 100, 1000, 10_000, 100_000]

# corridas de warm-up
out_flowpipe_pos, out_extrema_malla2_pos, fig_pos = generar_figura_malla2_envolventes_position(vecsamples=[1, 10]);
out_flowpipe_vel, out_extrema_malla2_vel, fig_vel = generar_figura_malla2_envolventes_velocity(vecsamples=[1, 10]);

# corridas largas
Random.seed!(1234)
out_flowpipe_pos, out_extrema_malla2_pos, fig_pos = generar_figura_malla2_envolventes_position(vecsamples=vecsamples);
out_flowpipe_vel, out_extrema_malla2_vel, fig_vel = generar_figura_malla2_envolventes_velocity(vecsamples=vecsamples);

# ---------------------------------------------
# Calculo de normas del flowpipe
# ---------------------------------------------

# dt siempre el mismo para cada tramo, y la ventana de tiempo va hasta 1 segundo,
# entonces el dt se saca de factor comun

# max(max_t |v_max(t)|, max_t |v_min(t)|)
out_flowpipe_vel_norma_inf_envolvente = max(abs(ρ([1.0], out_flowpipe_vel)), abs(ρ([-1.0], out_flowpipe_vel)))
# salida: 0.0012224968252669921

# δt * sum_{i=1}^N max(|v_max(t)|, |v_min(t)|)
δt = step_size_lamb(load_lambs(mesh=2), 1.0)
@assert tspan(out_flowpipe_vel[1]).hi ≈ δt
out_flowpipe_vel_norma_1_envolvente = sum(max(abs(ρ([1.0], set(X))), abs(ρ([-1.0], set(X)))) for X in out_flowpipe_vel) * δt
# salida: 0.0008133221497265329

# ---------------------------------------------
# Calculo de normas de numericas
# ---------------------------------------------

out_extrema_malla2_vel_norma_inf_envolvente = _norma_inf_envolvente_numericas(out_extrema_malla2_vel, vecsamples)
#=
6-element Vector{Float64}:
 0.00056983336205989
 0.0005752808846025316
 0.0005758510647681003
 0.0005821999057203168
 0.0006118261295338235
 0.0006220764784910799
 =#

out_extrema_malla2_vel_norma_1_envolvente = _norma_1_envolvente_numericas(out_extrema_malla2_vel, vecsamples)
#=
julia> out_extrema_malla2_vel_norma_1_envolvente = _norma_1_envolvente_numericas(out_extrema_malla2_vel, vecsamples)
6-element Vector{Float64}:
 9.268765778974471e-5
 0.00013515950042666908
 0.0001661339910196396
 0.00018520059379205858
 0.0001997982127229552
 0.00021423496225308065
 =#

# ------------------------------
# Guardar datos para las tablas
# ------------------------------

# convert the keys of a dictionary into strings; the values might themselves
# be dictionaries, in which case the operation is applied recursively
# moreover, all keys are appended a capital letter N, e.g. "1000", it is converted to "N1000"
function _to_str_keys(d::Dict)
    strkeys = "N" .* string.(keys(d))

    if first(d)[2] isa Dict
        return Dict(strkeys .=> _to_str_keys.(values(d)))
    else
        return Dict(strkeys .=> values(d))
    end
end

using MAT

results  = Dict{String, Any}("vecsamples" => vecsamples,
                             "dt" => δt,
                             "tiempo_flowpipe" => "8.504446 seconds (9.82 k allocations: 292.719 MiB, 1.10% gc time)",

                             "out_flowpipe_pos_lo" => [set(X).dat.lo for X in out_flowpipe_pos],
                             "out_flowpipe_pos_hi" => [set(X).dat.hi for X in out_flowpipe_pos],
                             "out_flowpipe_pos_tstart" => [tstart(X) for X in out_flowpipe_pos],
                             "out_flowpipe_pos_tend" => [tend(X) for X in out_flowpipe_pos],

                             "out_flowpipe_vel_lo" => [set(X).dat.lo for X in out_flowpipe_vel],
                             "out_flowpipe_vel_hi" => [set(X).dat.hi for X in out_flowpipe_vel],
                             "out_flowpipe_vel_tstart" => [tstart(X) for X in out_flowpipe_vel],
                             "out_flowpipe_vel_tend" => [tend(X) for X in out_flowpipe_vel],

                             "out_flowpipe_vel_norma_inf_envolvente" => out_flowpipe_vel_norma_inf_envolvente,
                             "out_flowpipe_vel_norma_1_envolvente" => out_flowpipe_vel_norma_1_envolvente,
                             "out_extrema_malla2_vel_norma_inf_envolvente" => out_extrema_malla2_vel_norma_inf_envolvente,
                             "out_extrema_malla2_vel_norma_1_envolvente" => out_extrema_malla2_vel_norma_1_envolvente,

                             "out_extrema_malla2_pos_times" => out_extrema_malla2_pos.idx,
                             "out_extrema_malla2_pos_U_min" => out_extrema_malla2_pos.U_min,
                             "out_extrema_malla2_pos_U_max" => out_extrema_malla2_pos.U_max,
                             "out_extrema_malla2_pos_Udot_min" => out_extrema_malla2_pos.U′_min,
                             "out_extrema_malla2_pos_Udot_max" => out_extrema_malla2_pos.U′_max,
                             "out_extrema_malla2_pos_hist" => _to_str_keys(out_extrema_malla2_pos.hist),
                             "out_extrema_malla2_pos_t" => collect(out_extrema_malla2_pos.t),
                             "out_extrema_malla2_pos_times" => _to_str_keys(out_extrema_malla2_pos.times),

                             "out_extrema_malla2_vel_times" => out_extrema_malla2_vel.idx,
                             "out_extrema_malla2_vel_U_min" => out_extrema_malla2_vel.U_min,
                             "out_extrema_malla2_vel_U_max" => out_extrema_malla2_vel.U_max,
                             "out_extrema_malla2_vel_Udot_min" => out_extrema_malla2_vel.U′_min,
                             "out_extrema_malla2_vel_Udot_max" => out_extrema_malla2_vel.U′_max,
                             "out_extrema_malla2_vel_hist" => _to_str_keys(out_extrema_malla2_vel.hist),
                             "out_extrema_malla2_vel_t" => collect(out_extrema_malla2_vel.t),
                             "out_extrema_malla2_vel_times" => _to_str_keys(out_extrema_malla2_vel.times)
                             );

matwrite("results.mat", results)

# ---------------
# Generar figura
# ---------------

generar_figura_malla2_envolvente_velocidad()

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
using Reproduce

function reproduce_config_experiment(config_file::AbstractString; tldr="", save_path="")
    experiment = Experiment(config_file, save_path)
    create_experiment_dir(experiment; tldr=tldr)
    add_experiment(experiment)
    ret = job(experiment)
    post_experiment(experiment, ret)
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module AsciiMath

include("Tokens.jl")
include("Lexer.jl")

end

max_stars_repo_name: MhmCats/AsciiMath

max_issues_repo_name: MhmCats/AsciiMath

max_forks_repo_name: MhmCats/AsciiMath


text:
# ---
# title: 301. Remove Invalid Parentheses
# id: problem301
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/remove-invalid-parentheses/description/>
# hidden: true
# ---
#
# Remove the minimum number of invalid parentheses in order to make the input
# string valid. Return all possible results.
#
# **Note:**  The input string may contain letters other than the parentheses `(`
# and `)`.
#
# **Example 1:**
#
#
#
#     Input: "()())()"
#     Output: ["()()()", "(())()"]
#
#
# **Example 2:**
#
#
#
#     Input: "(a)())()"
#     Output: ["(a)()()", "(a())()"]
#
#
# **Example 3:**
#
#
#
#     Input: ")("
#     Output: [""]
#
#
#
## @lc code=start
using LeetCode

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


text:
function cum_loans_by_qual_year()

end

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


text:
# include("../activate_env.jl")

using JLD2
using model1d1d

#= ###########################################################################
Example for Simulation Setup
########################################################################### =#
casename = "my_example_case"

sdm_enabled = false
front_rad_enabled = true
cross_rad_enabled = true

logpath="./"
logfile = joinpath(logpath, casename*"_log.txt")
new_log(logfile)

desc = casename * ", logfile: ", logfile

submodel = model1d1d.SimSettingsSubModels(
        scrap_discretization_model_enabled = false,
        cross_radiation_model_enabled = true,
        furnace_side_rad_model =  model1d1d.SimpleRadiationModel(900+273.15)
        )

gas_scrap_htc = model1d1d.ConstantHeatTransfer(50.0) # in W/(m² K)

#= ###########################################################################
Process Settings
########################################################################### =#
shaft_length  = 7.0 # in m
shaft_width   = 3.0 # in m
shaft_height  = 3.0 # in m

shaft_cross_sect_area = shaft_width * shaft_height

duration = 50.0 * 60.0 # in seconds (48 min) - old 40 / 2

# save results every t_save seconds
t_save = 30.0

# setting general simulation settings
sim_phases = [model1d1d.SimSettings(duration, submodel, t_save)]

#= ###########################################################################
Hot-Gas Settings
########################################################################### =#
gas_volume_flow = 19500 # NCm_per_h

gas_inlet_temperature =  1200 + 273.15 # K

# Material data
gas_density = model1d1d.ReciprocalMaterialPropertyOfX(360.47962061763553, 0.0,
        name = "density", unit = "kg/m^3")

gas_cp_true = model1d1d.PolynomialMaterialPropertyOfX([908.692796686028,
    0.336368996611928, -6.950621118017843e-5], name = "true heat capacity",
    unit = "J/(kg K)", bounds = (1000.0,1500.0))

gas_lambda = model1d1d.PolynomialMaterialPropertyOfX([0.007187645849107943,
        5.784085171400682e-5, 5.428742980927024e-9, -4.200871592280629e-12],
        name = "heat conductivity", unit = "W/(m K)", bounds = (0.01,0.15))

gas_mat = model1d1d.IncompressibleFluidThermalMaterial("off-gas", gas_density,
        gas_cp_true, gas_lambda)

# mass flow for norm temperature
gas_mass_flow = gas_volume_flow/3600 * model1d1d.pX(gas_mat.rho, 25 + 273.15)

offgas = model1d1d.OffGas(gas_mass_flow, gas_inlet_temperature, gas_mat)

#= ###########################################################################
Bulk Settings
########################################################################### =#
# complete mass of charged scrap
scrap_mass_charged = 111.1 * 1000.0 # in kg

scrap_bulk_mass_flow = scrap_mass_charged / duration

scrap_charged_temperature = 25.0 + 273.15 # in K

scrap_bulk_density = 801.0	# kg / m³ at RT
scrap_avg_thickness_of_piece = 0.01 # in m
scrap_avg_length_of_piece = 0.5 # in m
scrap_avg_width_of_piece = 0.4 # in m

# Material data
scrap_density = model1d1d.ConstantMaterialPropertyOfX(7850.0, name = "density",
        unit ="kg/m^3")

scrap_cp_true = model1d1d.ConstantMaterialPropertyOfX(500.0,  name = "true heat capacity",
unit = "J/(kg K)")

scrap_lambda = model1d1d.ConstantMaterialPropertyOfX(30.0, name = "heat conductivity"
, unit = "W/(m K)")

scrap_epsilon = model1d1d.ConstantMaterialPropertyOfX(0.8, name = "emissivity",
        unit ="kg/m^3")

scrap_mat = model1d1d.SolidThermalMaterial("scrap", scrap_density,
        scrap_cp_true, scrap_lambda, scrap_epsilon)

print_log(logfile, "\nBulk properties:")
print_log(logfile, ["\t Scrap bulk density: ", scrap_bulk_density])
print_log(logfile, ["\t Scrap avgerage thickness of piece: ", scrap_avg_thickness_of_piece])
print_log(logfile, ["\t Scrap avgerage length of piece: ", scrap_avg_length_of_piece])
print_log(logfile, ["\t Scrap avgerage width of piece: ", scrap_avg_width_of_piece])

scrap_bulk_porosity = 1.0 - (scrap_bulk_density/model1d1d.pX(scrap_mat.rho, 25 + 273.15))

print_log(logfile, ["\t Scrap bulk porosity: ", scrap_bulk_porosity])

scrap_piece = model1d1d.CharScrapPiece(scrap_avg_thickness_of_piece,
        scrap_avg_length_of_piece, scrap_avg_width_of_piece, scrap_mat)

bulk = model1d1d.Bulk(scrap_bulk_mass_flow,
        scrap_charged_temperature, shaft_length, shaft_width*shaft_height,
        scrap_bulk_density, scrap_bulk_porosity, 0.95, scrap_piece)

#= ###########################################################################
Numerical Settings
########################################################################### =#

no_cells_over_shaft = 120
no_cells_over_scrap_piece = nothing

mesh_settings = model1d1d.M1d1dMeshSettings(no_cells_over_shaft, no_cells_over_scrap_piece)

#= ###########################################################################
Final Model
########################################################################### =#

model = model1d1d.M1d1d(sim_phases, [gas_scrap_htc], [bulk],
        [offgas], [mesh_settings]);

#= ###########################################################################
Solution
########################################################################### =#
results = model1d1d.solve_1d1d_diffeq(model);

#= ###########################################################################
Save Results
########################################################################### =#

@save joinpath("./results", casename*".jld2") model desc results


text:
module CanonicalizeRequest

using ..Messages
using ..Strings: tocameldash

export canonicalizelayer

"""
    canonicalizelayer(req) -> HTTP.Response

Rewrite request and response headers in Canonical-Camel-Dash-Format.
"""
function canonicalizelayer(handler)
    return function(req; canonicalize_headers::Bool=false, kw...)
        if canonicalize_headers
            req.headers = canonicalizeheaders(req.headers)
        end
        res = handler(req; kw...)
        if canonicalize_headers
            res.headers = canonicalizeheaders(res.headers)
        end
        return res
    end
end

canonicalizeheaders(h::T) where {T} = T([tocameldash(k) => v for (k,v) in h])

end # module CanonicalizeRequest