text:
import Base.rand
import Base.Filesystem
import LinearAlgebra.det
import LinearAlgebra.dot
import LinearAlgebra.norm
import Plots
import JSON
using Logging
using LoggingExtras
import Glob.glob

#------------------------------------------------------------------------------
# Datastructures and functions related to delaunay triangulation
# - Edge: structure used to describe an edge
# - Triangule: structure used to describe a triangle
# - create_nodes: create a set of random nodes
# - in_circumcircle: find if a node is in the circumcircle of a triangle
# - is_clockwise: find if a set of 3 nodes is ordered clockwise
# - create_supertriangle: create a triangle containing all the nodes
# - require_flip: tell if two neighbor triangles need to be flipped because
#   they don't follow the criterion on opposite angles
# - compute_angle: compute the angle (in radians) between two edges
# - get_tri_node_ids: get the ids of the nodes of a triangle
# - get_triangulation_nodes: get the ids of the nodes used in a triangulation
# - flip_triangulation!: flip the triangles of the triangulation when needed
# - save_step: save the details of a step of the triangulation to a file
# - plot_step: create a visualization for a step of the triangulation
# - inside_triangle: find if a node is located inside a triangle
# - triangulate: create the triangulation of a set of nodes
#------------------------------------------------------------------------------

struct Edge
    # Note: possible improvement using static arrays to nodes and in_triangles
    #       in_triangles::StaticArrays.SVector{2,Int}

    # indices of the nodes
    nodes::Array{Int, 1}
    # indices of the triangles containing this edge
    in_triangles::Array{Int, 1}
end

struct Triangle
    # Note: possible improvement using static arrays to nodes and in_triangles
    #       edges::StaticArrays.SVector{3,Int}

    # indices of the edges of the triangle
    edges::Array{Int, 1}
end

function create_nodes(number::Int, x_bounds, y_bounds)
    # number:   number of nodes created
    # x_bounds: X bounds of nodes - list of 2 reals
    # y_bounds: Y bounds of nodes - list of 2 reals
    nodes = [[rand(Float64), rand(Float64)] for i in 1:number]
    # apply X and Y bounds
    for i=1:number
        nodes[i][1] = nodes[i][1]*(x_bounds[2] - x_bounds[1]) + x_bounds[1]
        nodes[i][2] = nodes[i][2]*(y_bounds[2] - y_bounds[1]) + y_bounds[1]
    end
    return nodes
end

function in_circumcircle(triangle_ids, id, nodes)
    # Tell if node with if "id" is in the circumcircle of triangle
    # made of the nodes of "triangle_ids"
    # - triangle_ids : id of the nodes of the triangle - list of 3 int
    # - id : id of the node investigated - int
    # - nodes : coordinates of the nodes - list of [x, y] arrays of floats
    if (length(triangle_ids) != 3)
        error("Argument triangle_ids must be a list of 3 int")
    end
    # vertices of the triangle
    A = nodes[triangle_ids[1]]
    B = nodes[triangle_ids[2]]
    C = nodes[triangle_ids[3]]
    P = nodes[id]
    # switch vertices if the triangle is clockwise
    if (is_clockwise(A, B, C))
        B, C = C, B
    end
    # create the matrix (see https://fr.wikipedia.org/wiki/Triangulation_de_Delaunay)
    mat = zeros(Float64, 3, 3)
    mat[1, 1] = A[1] - P[1]
    mat[2, 1] = B[1] - P[1]
    mat[3, 1] = C[1] - P[1]
    mat[1, 2] = A[2] - P[2]
    mat[2, 2] = B[2] - P[2]
    mat[3, 2] = C[2] - P[2]
    mat[1, 3] = A[1]^2 - P[1]^2 + A[2]^2 - P[2]^2
    mat[2, 3] = B[1]^2 - P[1]^2 + B[2]^2 - P[2]^2
    mat[3, 3] = C[1]^2 - P[1]^2 + C[2]^2 - P[2]^2
    return (det(mat) > 0)
end

function is_clockwise(A, B, C)
    # A, B & C: list of 2 reals representing 2D points
    AB = B - A
    AC = C - A
    determinant = AB[1]*AC[2] - AB[2]*AC[1]
    return (determinant < 0)
end

function create_supertriangle(nodes)
    # Create a triangle containing all the nodes
    # - nodes : coordinates of the nodes - list of [x, y] arrays of floats

    # bounds of the nodes
    nodes_x = [a_node[1] for a_node in nodes]
    nodes_y = [a_node[2] for a_node in nodes]
    x_bounds = [minimum(nodes_x), maximum(nodes_x)]
    y_bounds = [minimum(nodes_y), maximum(nodes_y)]

    # take a margin (relative to the X & Y bounds)
    m = 0.05
    min_length = 1e-4
    lx = x_bounds[2] - x_bounds[1] + min_length
    x_bounds = [x_bounds[1]-m*lx, x_bounds[2]+m*lx]
    ly = y_bounds[2] - y_bounds[1] + min_length
    y_bounds = [y_bounds[1]-m*ly, y_bounds[2]+m*ly]

    # get the circle containing the nodes
    right_top = [x_bounds[2], y_bounds[2]]
    center = [0.5*(x_bounds[1]+x_bounds[2]), 0.5*(y_bounds[1]+y_bounds[2])]
    radius = sqrt( (right_top[1]-center[1])^2 + (right_top[2]-center[2])^2 )

    # build an equilateral triangle with edges adjacent to the circle
    # base's center = lowest point of the circle on the Y axis
    base_mid = [center[1], center[2]-radius]
    # radius projection on an edge other than the base
    edge_mid = [center[1]+radius*cos(pi/6), center[2]+radius*sin(pi/6)]
    edge_l = 4*abs(edge_mid[1] - base_mid[1])
    tri_lower_right = [base_mid[1]+edge_l/2, base_mid[2]]
    tri_lower_left = [base_mid[1]-edge_l/2, base_mid[2]]
    tri_top = [base_mid[1], base_mid[2]+edge_l*sin(pi/3)]

    return [tri_lower_right, tri_top, tri_lower_left]
end

function require_flip(tri1_id, tri2_id, triangles, edges, nodes)
    # Evaluate if the triangles must be flipped, that is to say if the opposite angles don't
    # match the delaunay criterion (sum <= 180°)
    # - tri1_id : id of the first triangle
    # - tri2_id : id of the second triangle
    # - triangles : array(n) containing all the triangles
    # - edges : array(n) containing all the edges
    # - nodes : coordinates of the all nodes - list of [x, y] arrays of floats

    triangles_pair = [triangles[tri1_id], triangles[tri2_id]]

    # common_edge
    common_edge = 0
    for triangle1_edge in triangles_pair[1].edges
        if triangle1_edge in triangles_pair[2].edges
            common_edge = triangle1_edge
            break
        end
    end
    if common_edge == 0
        error("Common edge not found")
    end

    # angle opposite to the common edge
    angles = [NaN, NaN]
    for i in [1, 2]
        angle_edges = [edges[e] for e in triangles_pair[i].edges if e != common_edge]
        angle_node = intersect(angle_edges[1].nodes, angle_edges[2].nodes)[1]

        # ensure that the first point of both edges is their common point
        if (angle_edges[1].nodes[1] != angle_node)
            angle_edges[1] = Edge(reverse(angle_edges[1].nodes), [])
        end
        if (angle_edges[2].nodes[1] != angle_node)
            angle_edges[2] = Edge(reverse(angle_edges[2].nodes), [])
        end

        # compute angle
        angles[i] = compute_angle(angle_edges[1], angle_edges[2], nodes)
    end
    return (abs(angles[1]) + abs(angles[2]) > pi)
end

function compute_angle(edge1, edge2, nodes)
    # - edge1 : first edge (type Edge)
    # - edge2 : second edge (type Edge)
    # - nodes : coordinates of the all nodes - list of [x, y] arrays of floats
    p1 = nodes[edge1.nodes[1]]
    p2 = nodes[edge1.nodes[2]]
    v1 = [p2[1]-p1[1], p2[2]-p1[2]]
    v1 = v1/norm(v1)

    p1 = nodes[edge2.nodes[1]]
    p2 = nodes[edge2.nodes[2]]
    v2 = [p2[1]-p1[1], p2[2]-p1[2]]
    v2 = v2/norm(v2)

    return acos(dot(v1, v2))
end

function get_tri_node_ids(triangle_id, triangles, edges)
    # - triangle_id : id of the triangle
    # - triangles : array(n) containing all the triangles
    # - edges : array(n) containing all the edges
    triangle_edges = [edges[edge_id] for edge_id in triangles[triangle_id].edges]
    all_nodes = [e.nodes for e in triangle_edges]
    return unique(Iterators.flatten(all_nodes))
end

function get_triangulation_nodes(edges)
    # - edges : array(n) containing all the edges
    node_ids = unique(Iterators.flatten([e.nodes for e in values(edges)]))
    return node_ids
end

function flip_triangulation!(triangles, edges, nodes, do_save_steps, step_file_pattern; init_flip_step=0)
    # Flip the triangles of the triangulation that don't passe delaunay's criterion
    # This function modifies the triangles and edges input arguments
    # - triangles : array(n) containing all the triangles
    # - edges : array(n) containing all the edges
    # - nodes : coordinates of the all nodes - list of [x, y] arrays of floats

    flip_step = init_flip_step
    nb_edges_flipped = 0
    for edge_id in keys(edges)
        common_edge = edges[edge_id]
        common_nodes = common_edge.nodes
        # find if neighbor triangles separated by this edge must be flipped
        neighbors_id = common_edge.in_triangles
        if 0 in neighbors_id
            # edge is on the border of the triangulation: no flip
            continue
        elseif require_flip(neighbors_id[1], neighbors_id[2], triangles, edges, nodes)
            flip_step = flip_step + 1
            nb_edges_flipped = nb_edges_flipped + 1
            # do flip
            # 1. analyse the neighbor triangles: opposite nodes and id of other edges
            tri1_nodes = get_tri_node_ids(neighbors_id[1], triangles, edges)
            tri2_nodes = get_tri_node_ids(neighbors_id[2], triangles, edges)
            @debug "triangle 1 : $(neighbor_tri_ids[1]) has nodes $(tri1_nodes)"
            @debug "triangle 2 : $(neighbor_tri_ids[2]) has nodes $(tri2_nodes)"
            opposite_nodes = [setdiff(tri1_nodes, common_nodes)[1], setdiff(tri2_nodes, common_nodes)[1]]
            #top_nodes = setdiff(union(triangles[neighbors_id[1]].nodes, triangles[neighbors_id[1]].nodes)
            #       common_edge.nodes)
            @debug "opposite_nodes : $(opposite_nodes)"
            if length(unique(opposite_nodes)) != 2
                @error "Flip: same opposite nodes"
                error("Flip: same opposite nodes")
            end
            other_edges_ids = setdiff(triangles[neighbors_id[1]].edges, edge_id)
            append!(other_edges_ids, setdiff(triangles[neighbors_id[2]].edges, edge_id))
            # TODO: rewrite with union
            # other_edges_ids = setdiff(union(triangles[neighbors_id[1]].edges, triangles[neighbors_id[2]].edges), common_edge)

            # 2. reconstruct the common edge and the neighbor triangles
            #    -> the common edge is now built with the opposite nodes
            #    -> the nodes of the former common edge are now the new opposite nodes

            # 2.1 common edge
            common_edge = Edge(opposite_nodes, neighbors_id)
            @debug "Common edge $(common_edge_id) updated: nodes = $(opposite_nodes)"
            edges[edge_id] = common_edge

            # 2.2 triangles
            for i_tri in 1:2
                @debug "Rework triangle $(i_tri)"
                tri_id = neighbors_id[i_tri]
                # use a (former) common node and the new common edge
                top_node = common_nodes[i_tri]
                @debug "- Top node = $(top_node)"
                side_edges_ids = [e_id for e_id in other_edges_ids if top_node in edges[e_id].nodes]
                @debug "- Side edge 1 = $(side_edges_ids[1]): $(edges[side_edges_ids[1]])"
                @debug "- Side edge 2 = $(side_edges_ids[2]): $(edges[side_edges_ids[2]])"
                # - create new triangle (replace existing one with same id)
                triangles[tri_id] = Triangle([edge_id, side_edges_ids[1], side_edges_ids[2]])
                @debug "- Triangle $(tri_id) updated: $(triangles[tri_id].edges)"
                # - update the neighbor information (in_triangle) for the side edges
                for i_edge in side_edges_ids
                    # before the flip: side edge is in
                    # - one of the triangles that will be flipped: an id in in neighbor_tri_ids
                    # - a triangle out of the parent triangle unchanged : outer_tri
                    #   -> outer_tri_id is 0 is this outer triangle does not exist
                    # after the flip: side edge is in
                    # - the triangle that has been flipped: id = tri_id
                    # - the outer triangle with id outer_tri_id
                    # => replace id of triangle-to-be-flipped by id of new triangle
                    in_triangles = copy(edges[i_edge].in_triangles)
                    if (in_triangles[1] in neighbors_id)
                        in_triangles[1] = tri_id
                    elseif (in_triangles[2] in neighbors_id)
                        in_triangles[2] = tri_id
                    else
                        @error "- Edge $(i_edge) is not in flipped triangle $(tri_id)"
                        error("- Edge ", i_edge, " is not in flipped triangle ", tri_id)
                    end
                    edge_nodes = edges[i_edge].nodes
                    #edges[i_edge] = Edge(edge_nodes, [tri_id, outer_tri_id])
                    edges[i_edge] = Edge(edge_nodes, in_triangles)
                    @debug "- Edge $(i_edge) now in triangles $(in_triangles)"
                end
            end

            # save result
            if do_save_steps
                new_edges_ids = union(other_edges_ids, edge_id)
                save_step(step_file_pattern, "flip_$(flip_step)", nodes, 0, edges,
                new_edges_ids, [], [], triangles)
            end
        end
    end
    if nb_edges_flipped > 0
        @info "$(nb_edges_flipped) on $(length(edges)) edges have been flipped"
    end

    # check if there remain edges to flip after this first step
    nb_edges_to_flip = 0
    for edge_id in keys(edges)
        common_edge = edges[edge_id]
        # find if neighbor triangles separated by this edge must be flipped
        neighbors_id = common_edge.in_triangles
        if 0 in neighbors_id
            # edge is on the border of the triangulation: no flip
            continue
        elseif require_flip(neighbors_id[1], neighbors_id[2], triangles, edges, nodes)
            nb_edges_to_flip = nb_edges_to_flip + 1
        end
    end
    if (nb_edges_to_flip > 0)
        @info "Recursive call to flip_triangulation to fix $(nb_edges_to_flip) edges"
        flip_triangulation!(triangles, edges, nodes, do_save_steps, step_file_pattern, init_flip_step=flip_step)
    end

end

function save_step(file_pattern, file_suffix, nodes, last_node, edges, last_edges,
                    supertriangle_node_ids, supertriangle_edge_ids, triangles)
    # Create the name of the step from the pattern by inserting a suffix
    # between the pattern and the file extension
    (base, ext) = splitext(file_pattern)
    file = base*file_suffix*ext
    output = Dict("nodes" => nodes,
                  "last_node" => last_node,
                  "edges" => edges,
                  "last_edges" => last_edges,
                  "supertriangle_node_ids" => supertriangle_node_ids,
                  "supertriangle_edge_ids" => supertriangle_edge_ids,
                  "triangles" => triangles,
                  )
    open(file, "w") do io
        write(io, JSON.json(output))
    end
end

function plot_step(step_file, show_non_delaunay_triangles)
    # plot the nodes, edges and triangles for a step
    # warning: an existing file with the same name will be overwritten
    io = open(step_file, "r")
    string_content = read(io, String)
    close(io)
    json_content = JSON.parse(string_content)
    nodes = json_content["nodes"]
    last_node = json_content["last_node"]
    edges = Dict{Int, Edge}()
    for item in collect(json_content["edges"])
        key = parse(Int, item[1])
        edges[key] = Edge(item[2]["nodes"], item[2]["in_triangles"])
    end
    last_edges = json_content["last_edges"]
    supertriangle_edge_ids = json_content["supertriangle_edge_ids"]
    supertriangle_node_ids = json_content["supertriangle_node_ids"]
    triangles = Dict{Int, Triangle}()
    for item in collect(json_content["triangles"])
        key = parse(Int, item[1])
        triangles[key] = Triangle(item[2]["edges"])
    end

    # create a node to init figure - it will be reploted after
    Plots.scatter([nodes[1][1]], [nodes[1][2]], aspect_ratio=:equal, label="", color=:red)

    # non-delaunay triangles
    if show_non_delaunay_triangles
        for (tri_id, tri) in triangles
            tri_node_ids = get_tri_node_ids(tri_id, triangles, edges)
            # find if any node used in the triangulation is in this triangle's circumcircle
            for lookup_node_id in get_triangulation_nodes(edges)
                if !(lookup_node_id in tri_node_ids)
                    if in_circumcircle(tri_node_ids, lookup_node_id, nodes)
                        # the triangle is non-delaunay -> plot it
                        tri_points = [Tuple(nodes[i]) for i in tri_node_ids]
                        #append!(shapes, Plots.Shape(tri_points))
                        Plots.plot!([Plots.Shape(tri_points)], color=:red, label="")
                        # end nodes lookup
                        break
                    end
                end
            end
        end
    end

    # edges
    for e in edges
        if e[1] in last_edges
            # this edge will be plotted with the last edges - skip it
            continue
        elseif e[1] in supertriangle_edge_ids
            # this edge will be plotted with the supertriangle edges - skip it
            continue
        end
        edge_x = [nodes[n][1] for n in e[2].nodes]
        edge_y = [nodes[n][2] for n in e[2].nodes]
        Plots.plot!(edge_x, edge_y, label="", color=:black)
    end

    # supertriangle's edges (if any)
    for edge_id in supertriangle_edge_ids
        if edge_id in last_edges
            # this edge will be plotted with the last edges - skip it
            continue
        end
        edge_nodes = edges[edge_id].nodes
        edge_x = [nodes[n][1] for n in edge_nodes]
        edge_y = [nodes[n][2] for n in edge_nodes]
        Plots.plot!(edge_x, edge_y, label="", color=:chocolate3)
    end

    # last edges (if any)
    for edge_id in last_edges
        edge_nodes = edges[edge_id].nodes
        edge_x = [nodes[n][1] for n in edge_nodes]
        edge_y = [nodes[n][2] for n in edge_nodes]
        Plots.plot!(edge_x, edge_y, label="", color=:green4)
    end

    # nodes
    nodes_x = [n[1] for n in nodes]
    nodes_y = [n[2] for n in nodes]
    Plots.scatter!(nodes_x, nodes_y, label="", color=:deepskyblue)

    # supertriangle's nodes
    nodes_x = [nodes[i][1] for i in supertriangle_node_ids]
    nodes_y = [nodes[i][2] for i in supertriangle_node_ids]
    Plots.scatter!(nodes_x, nodes_y, label="", color=:burlywood1)

    # last node inserted
    if last_node > 0
        Plots.scatter!([nodes[last_node][1]], [nodes[last_node][2]], label="", color=:red)
    end

    # title
    figure_title = basename(splitext(step_file)[1])
    Plots.title!(figure_title)

    # save fig
    figure_file = splitext(step_file)[1]*".png"
    #rm(figure_file, force=true)
    Plots.savefig(figure_file)
end

function inside_triangle(triangle_ids, node_id, nodes)
    # Evaluate if a node is inside a triangle using barycentric coordinates method
    # https://en.wikipedia.org/wiki/Barycentric_coordinate_system
    #   See #Edge_approach for the calculation of the barycentric coordinates
    #   and #Determining_location_with_respect_to_a_triangle for the inside/outside criterion
    # - triangle_ids : id of the nodes of the triangle - list of 3 int
    # - node_id : id of the node investigated - int
    # - nodes : coordinates of the nodes - list of nodes coordinates [x, y]
    (x, y) = nodes[node_id]
    (x1, y1) = nodes[triangle_ids[1]]
    (x2, y2) = nodes[triangle_ids[2]]
    (x3, y3) = nodes[triangle_ids[3]]
    denominator = (y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3)
    lambda1 = ((y2 - y3)*(x - x3) + (x3 - x2)*(y - y3))/denominator
    lambda2 = ((y3 - y1)*(x - x3) + (x1 - x3)*(y - y3))/denominator
    lambda3 = 1 - lambda1 - lambda2
    inside = (0 < lambda1 < 1) && (0 < lambda2 < 1) && (0 < lambda3 < 1)
    return inside
end

function triangulate(input_nodes, must_save_steps=false, step_file_pattern="./steps/step_.txt")
    # - input_nodes : coordinates of the all nodes - list of [x, y] arrays of floats
    # - must_save_steps : save each step to a file (named from pattern)
    # - step_file_pattern : pattern of the step files - ex: "./steps/step_.txt"

    # nodes: work on a shallow copy of input nodes (we will add then remove the 3 supernodes)
    nodes = copy(input_nodes)
    number_of_nodes = length(nodes)
    if number_of_nodes < 3
        # TODO: return empty set instead of printing an error
        @error "The number of nodes must be 3 or more"
        error("The number of nodes must be 3 or more")
    end

    # remove existing step files
    if (must_save_steps)
        @info "Removing existing step files"
        step_file_pattern_basename = splitext(basename(step_file_pattern))[1]
        for step_file in readdir(dirname(step_file_pattern), join=true)
            if startswith(basename(step_file), step_file_pattern_basename)
                rm(step_file)
                @debug "Removed: "*step_file
            end
        end
    end

    # create dicts to store the triangles (1 triangle is described by the indices of its 3 nodes)
    # notice: dicts are used instead of arrays to avoid renumbering elements each time an element
    #         is removed
    triangles = Dict{Int, Triangle}()
    current_triangle_index = 0
    edges = Dict{Int, Edge}()
    current_edge_index = 0
    # create the super triangle
    # -> its id is 1
    # -> its nodes have ids n+1, n+2 and n+3
    # -> its edges are numbered 1, 2 and 3
    supertriangle_nodes = create_supertriangle(nodes)
    supertriangle_node_ids = [number_of_nodes+1, number_of_nodes+2, number_of_nodes+3]
    append!(nodes, supertriangle_nodes)

    edge1 = Edge([supertriangle_node_ids[1], supertriangle_node_ids[2]], [1, 0])
    edges[1] = edge1
    edge2 = Edge([supertriangle_node_ids[2], supertriangle_node_ids[3]], [1, 0])
    edges[2] = edge2
    edge3 = Edge([supertriangle_node_ids[3], supertriangle_node_ids[1]], [1, 0])
    edges[3] = edge3
    current_edge_index = 3

    supertriangle_edge_ids = [1, 2, 3]
    supertriangle = Triangle(supertriangle_edge_ids)
    triangles[1] = supertriangle
    current_triangle_index = 1

    # iterate over the nodes to insert
    for node_id=1:number_of_nodes
        @info "Insert node $(node_id)"

        # find the parent triangle : triangle which circumcircle contains this node
        # --------------------------
        parent_tri_id = 0
        for tri_id in keys(triangles)
            # get nodes of the triangle
            tri_nodes = get_tri_node_ids(tri_id, triangles, edges)
            if inside_triangle(tri_nodes, node_id, nodes)
                parent_tri_id = tri_id
                break
            end
            # Warning: what to do if the node is ON an edge of the triangle ?
        end
        if (parent_tri_id == 0)
            @error "parent triangle not found for node $(node_id)"
            error("parent triangle not found for node ", node_id)
        end
        # DEBUG: remove
        @debug "Id of the parent the triangle: $(parent_tri_id)"

        # update the triangulation: create the new triangles from the parent triangle and the inserted node
        # -------------------------
        # - edges of the parent triangle
        parent_edges = triangles[parent_tri_id].edges
        parent_nodes = get_tri_node_ids(parent_tri_id, triangles, edges)

        # - create the new edges (nodes of the parent tri + new node)
        new_edges = [Edge([parent_nodes[i], node_id], [0, 0]) for i in 1:3]
        new_edges_ids = [0, 0, 0]
        for i_new_edge in 1:3
            if length(unique(new_edges[i_new_edge].nodes)) != 2
                @error "triangulation: error edge with same node twice $(new_edges[i_new_edge])"
                error("triangulation: error edge with same node twice ", new_edges[i_new_edge])
            end
            current_edge_index = current_edge_index + 1
            new_edges_ids[i_new_edge] = current_edge_index
            edges[current_edge_index] = new_edges[i_new_edge]
        end

        # - create the new triangles with a bound edge and the two matching new edges
        new_triangles_indices = []
        for parent_edge_id in parent_edges
            @debug "Create new triangle from parent edge $(parent_edge_id) : $(edges[parent_edge_id])"
            #
            ref_edge = edges[parent_edge_id]
            ref_nodes = ref_edge.nodes
            # find the edges of the new triangle
            tri_edges_ids = [parent_edge_id]
            for i_edge=1:3
                if (new_edges[i_edge].nodes[1] in ref_nodes) || (new_edges[i_edge].nodes[2] in ref_nodes)
                    append!(tri_edges_ids, new_edges_ids[i_edge])
                end
            end
            if length(tri_edges_ids) != 3
                @error "Newly created triangle has not 3 edges"
                error("Newly created triangle has not 3 edges")
            end

            # create the triangle and add it to the table of triangles
            current_triangle_index = current_triangle_index + 1
            new_triangle_id = current_triangle_index
            append!(new_triangles_indices, new_triangle_id)
            triangles[current_triangle_index] = Triangle(tri_edges_ids)
            @debug "Create triangle #$(new_triangle_id) with edges $(tri_edges_ids)"

            # - update neighbor info of the edge from the parent triangle
            ref_edge_triangles = ref_edge.in_triangles
            replace!(ref_edge_triangles, parent_tri_id => new_triangle_id)
            edges[parent_edge_id] = Edge(ref_nodes, ref_edge_triangles)
        end

        # - update the new edges: neighbor triangles
        for i_new_edge in new_edges_ids

            # find the new triangle containing this edge
            tri_neighbors = []
            for tri_id in new_triangles_indices
                if i_new_edge in triangles[tri_id].edges
                    append!(tri_neighbors, tri_id)
                end
            end
            # update the edge
            @debug "Update edge #$(i_new_edge) : in_triangles = $(tri_neighbors)"
            edge_nodes = edges[i_new_edge].nodes
            edges[i_new_edge] = Edge(edge_nodes, tri_neighbors)
        end

        # - remove the parent triangle
        @debug "Delete parent triangle with id $(parent_tri_id)"
        pop!(triangles, parent_tri_id)

        # save step
        if must_save_steps
            save_step(step_file_pattern, string(node_id), nodes, node_id, edges,
                new_edges_ids, supertriangle_node_ids, supertriangle_edge_ids, triangles)
        end

        # flip the new triangles, if needed
        # ---------------------------------
        # we iterate over the edges that made the former parent triangle
        # -> for each, we need to see if the triangles on each side respect the angle criterion

        flip_step_file_pattern = splitext(step_file_pattern)[1]*string(node_id)*"_"*splitext(step_file_pattern)[2]
        flip_triangulation!(triangles, edges, nodes, must_save_steps, flip_step_file_pattern)
    end

    #--------------------------------------------------------------------------
    # Remove the data related to the (now deleted) supertriangle: nodes, edges
    #--------------------------------------------------------------------------

    @info "Cleanup"

    # find the edges containg a supernode
    edges_to_delete_id = []
    triangles_to_delete_id = []
    for (edge_id, edge) in edges
        if intersect(edge.nodes, supertriangle_node_ids) != []
            @debug "Edge $(edge_id) contains a node from the supertriangle -> to delete"
            # mark the edge for deleting
            append!(edges_to_delete_id, edge_id)
            # mark the triangle(s) made with this edge for deletion
            for tri_id in edge.in_triangles
                if tri_id != 0
                    append!(triangles_to_delete_id, tri_id)
                end
            end
        end
    end
    triangles_to_delete_id = unique(triangles_to_delete_id)

    # delete the edges
    for edge_id in edges_to_delete_id
        pop!(edges, edge_id)
    end

    # update information for the edges in triangles about to be deleted
    for (edge_id, edge) in edges
        in_deleted_tri_ids = intersect(edge.in_triangles, triangles_to_delete_id)
        # 3 cases for the current edge:
        # - in 2 triangles marked for deletion -> it has *already* been deleted -> nothing to do
        # - in 1 triangle marked for deletion -> replace its id by 0
        # - in 0 triangles marked for deleting -> nothing to do
        if length(in_deleted_tri_ids) == 1
            in_triangles_updated = replace(edge.in_triangles, in_deleted_tri_ids[1] => 0)
            edges[edge_id] = Edge(edge.nodes, in_triangles_updated)
        end
    end

    # delete the triangles
    for tri_id in triangles_to_delete_id
        @debug "Delete triangle $(tri_id)"
        pop!(triangles, tri_id)
    end

    # delete the nodes of the supertriangle: 3 last elements of "nodes"
    nodes = nodes[1:number_of_nodes]

    # write the final state
    if must_save_steps
        save_step(step_file_pattern, "final", nodes, 0, edges, [], [], [], triangles)
    end

    return edges, triangles
end

#------------------------------------------------------------------------------
# Functions related to the CLI application
# - parse_args: parses the options from a user's input string
# - print_help: displays the help of the CLI
# - get_arg_value: get the value of a couple of options (long/short) from
#   parse_args's output
# - get_main_parameters: get the CLI parameters from the user's input string
# - main: perform the operation requested by the user
#------------------------------------------------------------------------------

function parse_args(args_string::String)
    # parse the input arguments
    output = Dict()
    output["non_options"] = String[]
    option_regex = r"^-[a-zA-Z]"
    long_option_regex = r"^--[a-zA-Z]"
    # note: elements of the result of split are of type SubString{String}
    #       -> convert to String to ease future handling
    args = [string(s) for s in split(args_string)]
    index = 1
    while index <= length(args)
        currentArg = string(args[index])

        if currentArg == "-"
            error("Invalid argument '-'")
        end

        if occursin(option_regex, currentArg)
            # options found -> extract each one are store it in the dict
            number_of_options = length(currentArg) - 1
            if number_of_options == 1
                # single option
                # -> the next argument can hold a value, or be another option
                option = string(currentArg[2])
                if index < length(args)
                    nextArg = args[index+1]
                    if occursin(option_regex, nextArg)
                        # option: skip it
                        value = nothing
                    else
                        # value: store it
                        value = nextArg
                        # increment index because arg[index+1] as been dealt with
                        index = index + 1
                    end
                else
                    # no following argument: empty value
                    value = nothing
                end
                # store option
                output[option] = value
            else
                # several options stacked -> they cannot have values
                for option in currentArg[2:end]
                    output[string(option)] = nothing
                end
            end
        elseif occursin(long_option_regex, currentArg)
            # long option
            if occursin("=", currentArg)
                # the long option is passed with a value
                option, value = split(currentArg[3:end], "=")
                value = string(value)   # convert from SubString{String} to String to ease handling
                if value == ""
                    error("No value for option $(option)")
                end
            else
                option = currentArg[3:end]
                value = nothing
            end
            output[option] = value
        else
            # value not related to an option
            push!(output["non_options"], currentArg)
        end

        # go to the next index in the arguments array
        index = index + 1
    end

    return output
end

function print_help()
    println("Help: Delaunay triangulation of a set of nodes")
    println("-----")
    println("")
    println("Syntax for options: -s short_value --long-option=long_value")
    println("")
    println("Options:")
    println("-c               Create a random set of nodes. The value is the number of nodes")
    println("                 (default: 100")
    println("--create-nodes   Same as -c")
    println("-f               Name of the text file containing the nodes (1 node per line, ")
    println("                 x et y coordinates separated by a comma)")
    println("                 [!] NOT YET IMPLEMENTED")
    println("--from-file      Same as -f")
    println("-h               Display the help")
    println("--help           Same as -h")
    println("-l               When used with -p (or --plot), do not generate a triangulation ")
    println("                 and plot from a result file. The value is the name of the ")
    println("                 result file. Give a regex to plot several files.")
    println("--load-result    Same as -l")
    println("-o               Name of the output file.")
    println("--output         Same as -o")
    println("-p               Plot the triangulation. The possible values are:")
    println("                   'final'   (default) only plot the final triangulation")
    println("                   'all'     plot the final triangulation and each step")
    println("--plot           Same as -p")
    println("-s               Save the steps. The file names are built from the output: ")
    println("                 <output>_step<stepId>_flip<flipId>.<output_ext>")
    println("--save-steps     Same as -s")

end

function get_arg_value(options, args_dict)
    # Get the value of an option which can have several names
    # - options: list of names corresponding to the same option.
    #            eg: ["c", "create"] to get short option "-c" and long option "--create"
    # - args_dict: dictionary containing the options and their values (nothing if no value)
    # Return
    #  - false if none of the options is in the arguments dict
    #  - the value (string or nothing) of the option (or of all options is they match)
    #  - error if several options are in the dict with different values
    options_in_args = [n for n in options if n in keys(args_dict)]
    if options_in_args == []
        # None of the options is in the arguments dict -> use the default value
        return false
    end

    if length(options_in_args) == 1
        return args_dict[options_in_args[1]]
    else
        # Redundant options: check that they have the same values
        if length(unique([args_dict[o] for o in options_in_args])) == 1
            return args_dict[options_in_args[1]]
        else
            # error
            error("Options passed with different value: $(options_in_args)")
        end
    end
end

function get_main_parameters(args_string::String)

    # initialize the parameters' dict with default values
    params = Dict(
        "create_nodes" => false,
        "create_nodes_number" => 100,
        "load_node_file" => false,
        "load_result_files" => false,
        "output" =>  "./triangulation.json",
        "plot_results" => false,
        "save_steps" => false,
        "help" => false
    )

    # parse arguments
    args = parse_args(args_string)

    # handle invalid options: signal, set help and return
    all_options = ["c", "create-nodes", "f", "from-file", "h", "help", "l", "load-result",
                     "o", "output", "p", "plot", "s", "save-steps", "non_options"]
    has_invalid_option = false
    for key in keys(args)
        if ! (key in all_options)
            pop!(args, key)
            if length(key) == 1
                option_name = "-$(key)"
            else
                option_name = "--$(key)"
            end
            println("Invalid option $(option_name)")
            has_invalid_option = true
        end
    end
    if has_invalid_option === true
        params["help"] = true
        return params
    end

    # if there is no option, set help and return
    if length(keys(args)) == 1 && "non_options" in keys(args)
        params["help"] = true
        return params
    end

    # get the parameters from input arguments & default value

    # * create nodes: -c and --create_nodes
    arg_create_nodes = get_arg_value(["c", "create-nodes"], args)
    if arg_create_nodes != false
        # either "c" or "create-nodes" is passed
        params["create_nodes"] = true
        if arg_create_nodes !== nothing
            # passed with a value: convert it to int
            params["create_nodes_number"] = parse(Int, arg_create_nodes)
        end
    end

    # * get nodes from a file: -f and --from-file
    arg_nodes_file = get_arg_value(["f", "from-file"], args)
    if arg_nodes_file === nothing
        error("Options -f and --from-file must be used with a value")
    elseif arg_nodes_file != false
        if Filesystem.isfile(arg_nodes_file)
            params["load_node_file"] = arg_nodes_file
        else
            error("Node file passed to option -f or --from-file does not exist")
        end
    end

    # * display help
    arg_help = get_arg_value(["h", "help"], args)
    if arg_help != false
        params["help"] = true
    end

    # * get the results to plot from a file: -l and --load-result
    arg_load_result = get_arg_value(["l", "load-result"], args)
    if arg_load_result === nothing
        error("Options -l and --load-result must be used with a value")
    elseif arg_load_result != false
        # check the existence of files matching the pattern
        if length(glob(arg_load_result)) > 0
            params["load_result_files"] = arg_load_result
        else
            error("No file match the pattern given for option -l or --load-result")
        end
    end

    # * output file: -o and --output
    arg_output = get_arg_value(["o", "output"], args)
    if typeof(arg_output) === String
        # use the value given
        params["output"] = arg_output
    end

    # * plot results: -p and --plot
    arg_plot_results = get_arg_value(["p", "plot"], args)
    if arg_plot_results === nothing
        params["plot_results"] = "final"
    elseif typeof(arg_plot_results) == String
        if arg_plot_results in ["final", "all"]
            # use the value given
            params["plot_results"] = arg_plot_results
        else
            error("Valid values for -p and --plot are \"final\" and \"all\"")
        end
    end

    # * save the triangulation steps: -s and --save
    arg_save_steps = get_arg_value(["s", "save-steps"], args)
    if arg_save_steps !== false
        params["save_steps"] = true
    end

    return params
end

function main(args)
    # get parameters from the user's arguments
    args_string = join(args, " ")
    params = Dict()     # need to create the variable at this level, before the try/catch
    try
        params = get_main_parameters(args_string)
    catch ErrorException
        # something is invalid with the input parameters
        # -> print help and exit
        print_help()
        return
    end

    # apply a scenario, depending on the inputs
    if params["help"] === true
        # scenario 1: display help
        print_help()
    elseif params["load_result_files"] !== false
        # scenario 2: plot the results of a previous triangulation from result files
        step_file_pattern = params["load_result_files"]
        for filename in glob(step_file_pattern)
            plot_step(filename, true)
        end
    else
        # scenario 3: calculate a triangulation
        # - nodes
        nodes = []
        if params["load_node_file"] !== false
            # load the nodes from a file
            error("Load the nodes from a file: not implemented")
        elseif params["create_nodes"] === true
            # create random nodes
            nodes = create_nodes(params["create_nodes_number"], [-1., 1.], [-1., 1.])
        else
            # no policy for the nodes: notify error and exit
            println("The nodes must be either created or loaded from a file")
            println("  - created: -c or --create-nodes | optional value: number of nodes")
            println("  - loaded: -f or --from-file | value: name of the file")
            return
        end
        # - triangulation
        output = params["output"]
        if params["save_steps"] === true
            do_save_steps = true
            # create the pattern for the step files
            step_file_pattern = splitext(output)[1]*"_step_"*splitext(output)[2]
        else
            do_save_steps = false
            step_file_pattern = ""
        end
        (edges, triangles) = triangulate(nodes, do_save_steps, step_file_pattern)
        save_step(output, "", nodes, 0, edges, [], [], [], triangles)
        # - plot
        if params["plot_results"] !== false
            if params["plot_results"] == "final"
                plot_step(params["output"], true)
            elseif params["plot_results"] == "all"
                output_glob_pattern = splitext(output)[1]*"*"*splitext(output)[2]
                for filename in glob(output_glob_pattern)
                    plot_step(filename, true)
                end
            end
        end
    end
end
;

# CLI call from a terminal: call the main function with CLI args
if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end
;


text:
"""
`module  JAC.Defaults`
    ... a submodel of JAC that contains all constants and global parameters of the program.
"""
module Defaults

    using Dates,  JAC, ..Basics,  ..Radial,  ..Math
    # 2014 CODATA recommended values, obtained from http://physics.nist.gov/cuu/Constants/

    export  convertUnits, getDefaults,  setDefaults

    # Dimensionless constants
    const FINE_STRUCTURE_CONSTANT         = 7.297_352_566_4e-3
    const INVERSE_FINE_STRUCTURE_CONSTANT = 137.035_999_139

    # Constants in SI units
    const BOHR_RADIUS_SI                  = 0.529_177_210_67e-10
    const BOLTZMANN_CONSTANT_SI           = 1.380_648_52e-23
    const ELECTRON_MASS_SI                = 9.109_383_56e-31
    const ELEMENTARY_CHARGE_SI            = 1.602_176_620_8e-19
    const HARTREE_ENERGY_SI               = 4.359_744_650e-18
    const PLANCK_CONSTANT_SI              = 6.626_070_040e-34
    const PLANCK_CONSTANT_OVER_2_PI_SI    = 1.054_571_800e-34
    const SPEED_OF_LIGHT_IN_VACUUM_SI     = 299_792_458

    # Constants with eV instead of J
    const BOLTZMANN_CONSTANT_EV           = 8.617_330_3e-5
    const HARTREE_ENERGY_EV               = 27.211_386_02
    const PLANCK_CONSTANT_EV              = 4.135_667_662e-15
    const PLANCK_CONSTANT_OVER_2_PI_EV    = 6.582_119_514e-16

    # Constants with other atomic units
    const RYDBERG_IN_KAYSERS              = 109.73731534e3
    const ELECTRON_MASS_IN_G              = 9.1093897e-28
    const HBAR_IN_ERGS                    = 1.05457266e-27
    const ELECTRON_CHARGE_IN_ESU          = 4.80320680e-10

    # Constants in unified atomic mass units (u)
    const ELECTRON_MASS_U                 = 5.485_799_090_70e-4
    const NEUTRON_MASS_U                  = 1.008_664_915_88
    const PROTON_MASS_U                   = 1.007_276_466_879

    # Relationships between energy equivalents
    const ELECTRON_VOLT_ATOMIC_MASS_UNIT_RELATIONSHIP = 1.073_544_110_5e-9

    # Predefined conversion factors
    const CONVERT_ENERGY_AU_TO_EV           = HARTREE_ENERGY_EV
    const CONVERT_ENERGY_AU_TO_KAYSERS      = 2.0 * RYDBERG_IN_KAYSERS
    const CONVERT_ENERGY_AU_TO_PER_SEC      = 6.57968974479e15
    const CONVERT_INTENSITY_AU_TO_W_CM2     = 3.50944758e16
    const CONVERT_TIME_AU_TO_SEC            = 2.418_884_254e-17
    const CONVERT_CROSS_SECTION_AU_TO_BARN  = BOHR_RADIUS_SI^2 * 1.0e28
    const CONVERT_RATE_AU_TO_PER_SEC        = (ELECTRON_MASS_IN_G/HBAR_IN_ERGS) * ((ELECTRON_CHARGE_IN_ESU^2/HBAR_IN_ERGS)^2)
    const CONVERT_STRENGTH_AU_TO_BARN_EV    = CONVERT_CROSS_SECTION_AU_TO_BARN * HARTREE_ENERGY_EV
    const CONVERT_STRENGTH_AU_TO_CM2_EV     = BOHR_RADIUS_SI^2 * 1.0e4 * HARTREE_ENERGY_EV
    const CONVERT_LENGTH_AU_TO_FEMTOMETER   = BOHR_RADIUS_SI * 1.0e15

    # Predefined coefficients for numerical integration
    const FINITE_DIFFERENCE_NPOINTS         = 6
    const fivePointCoefficients             = Array{Float64}([2 * 7, 32, 12, 32, 7]) * 2 / 45  # * grid.h
    const newtonCotesCoefficients           = fivePointCoefficients

    n = FINITE_DIFFERENCE_NPOINTS
     const weights = Array{Float64}(undef, 2*n + 1, 2*n + 1)
    for i = 1:2*n + 1
        weights[i,:] = Math.finiteDifferenceWeights(-n + i - 1, 2*n + 1, order=1)[2,:]
    end

    #
    # Global settings that can be (re-) defined by the user.
    GBL_FRAMEWORK               = "relativistic"
    GBL_CONT_SOLUTION           = BsplineGalerkin()         ###  ContBessel(), ContSine(), AsymptoticCoulomb(), NonrelativisticCoulomb(), BsplineGalerkin()
    GBL_CONT_NORMALIZATION      = PureSineNorm()            ###  PureSineNorm(), CoulombSineNorm(), OngRussekNorm()
    GBL_QED_HYDROGENIC_LAMBDAC  = [1.0,  1.0,  1.0,  1.0,  1.0]
    GBL_QED_NUCLEAR_CHARGE      = 0.1
    GBL_WARNINGS                = String[]

    GBL_ENERGY_UNIT             = "eV"
    GBL_CROSS_SECTION_UNIT      = "barn"
    GBL_RATE_UNIT               = "1/s"
    GBL_STRENGTH_UNIT           = "cm^2 eV"
    GBL_TIME_UNIT               = "sec"

    GBL_PRINT_SUMMARY           = false
    GBL_PRINT_TEST              = false
    GBL_PRINT_DEBUG             = false

    GBL_SUMMARY_IOSTREAM        = nothing

    GBL_STANDARD_GRID           = Radial.Grid(true, printout=false)

    #
    # Global variables that are used for internal storage
    GBL_Storage_XL_Coulomb      = Dict{String, Float64}()
    GBL_Storage_XL_Breit        = Dict{String, Float64}()

    """
    `Defaults.convertUnits()`
        ... converts some data from one format/unit into another one; cf. the supported keystrings and return values.

    + `("cross section: from atomic to predefined unit", value::Float64)`  or  `("cross section: from atomic", value::Float64)`
        ... to convert an energy value from atomic to the predefined cross section unit; a Float64 is returned.

    + `("cross section: from atomic to barn", value::Float64)`  or  `("cross section: from atomic to Mbarn", value::Float64)`  or
      `("cross section: from atomic to Hz", value::Float64)`  or  `("energy: from atomic to Angstrom", value::Float64)`
        ... to convert an energy value from atomic to the speficied cross section unit; a Float64 is returned.

    + `("cross section: from predefined to atomic unit", value::Float64)`  or  `("cross section: to atomic", value::Float64)`
        ... to convert a cross section value from the predefined to the atomic cross section unit; a Float64 is returned.

    + `("einstein B: from atomic", value::Float64)`
        ... to convert a Einstein B coefficient from atomic to the speficied energy units; a Float64 is returned.or

    + `("energy: from atomic to eV", value::Float64)`  or  `("energy: from atomic to Kayser", value::Float64)`  or
      `("energy: from atomic to Hz", value::Float64)`  or  `("energy: from atomic to Angstrom", value::Float64)`
        ... to convert an energy value from atomic to the speficied energy unit; a Float64 is returned.

    + `("energy: from predefined to atomic unit", value::Float64)`  or  `("energy: to atomic", value::Float64)`... to convert an energy value
                                                    from the predefined to the atomic energy unit; a Float64 is returned.
    + `("energy: from eV to atomic unit", value::Float64)` ... to convert an energy value from eV to the atomic energy unit; a Float64 is returned.

    + `("energy: from wavelength [nm] to atomic", value::Float64)` ... to convert a wavelength [nm] to the atomic energy unit; a Float64 is returned.

    + `("intensity: from W/cm^2 to atomic", value::Float64)` ... to convert the intensity [in W/cm^2] to the atomic intensity unit; a Float64 is returned.

    + `("kinetic energy to wave number: atomic units", value::Float64)`  ... to convert a kinetic energy value (in a.u.) into a wave number
                                                    k (a.u.); a Float64 is returned.

    + `("kinetic energy to wavelength: atomic units", value::Float64)`  ... to convert a kinetic energy value (in a.u.) into a wavelength (a.u.);
                                                    a Float64 is returned.

    + `("length: from fm to atomic", value::Float64)`  ... to convert a length value (in fm) into a.u.;  a Float64 is returned.
    + `("length: from atomic to fm", value::Float64)`  ... to convert a length value (in Bohr's a.u.) into fm;  a Float64 is returned.

    + `("rate: from atomic to predefined unit", value::Float64)`  or  ("rate: from atomic", value::Float64)  ... to convert a rate value
                                                from atomic to the predefined rate unit; a Float64 is returned.

    + `("rate: from atomic to 1/s", value::Float64)`  ... to convert an rate value from atomic to the speficied rate unit; a Float64 is returned.

    + `("rate: from predefined to atomic unit", value::Float64)`  or  `("rate: to atomic", value::Float64)'... to convert a
                                                    rate value from the predefined to the atomic rate unit; a Float64 is returned.

    + `("strength: from atomic to predefined unit", value::Float64)`  or  ("strength: from atomic", value::Float64)  ... to convert a (resonance)
                                                    strength value from atomic to the predefined rate unit; a Float64 is returned.

    + `("time: from atomic to predefined unit", value::Float64)`  or  ("time: from atomic", value::Float64)  ... to convert an time value
                                                from atomic to the predefined time unit; a Float64 is returned.

    + `("time: from atomic to sec", value::Float64)`  or  `("time: from atomic to fs", value::Float64)'  or
        `("time: from atomic to as", value::Float64)`  ... to convert a time value from atomic to the speficied time unit; a Float64 is returned.

    + `("time: from predefined to atomic unit", value::Float64)`  or  `("time: to atomic", value::Float64)'  ... to convert a
                                                time value from the predefined to the atomic time unit; a Float64 is returned.

    + `("temperature: from Kelvin to (Hartree) units", value::Float64)`  ... to convert a temperature in Kelvin into atomic (Hartree) units;
                                                   a Float64 is returned.

    + `("temperature: from atomic to Kelvin", value::Float64)`  ... to convert an atomic (energy) unit into Kelvin; 1 Hartree = 315774.64 K;
                                                   a Float64 is returned.

    + `("wave number to total electron energy: atomic units", value::Float64)`  ... to convert a wavenumber (a.u.) into the total electron
                                                energy, including the rest energy; a Float64 is returned.

    + `("wave number to kinetic energy: atomic units", value::Float64)`  ... to convert a wavenumber (a.u.) into the kinetic energy;
                                                a Float64 is returned.
    """
    function convertUnits(sa::String, wa::Float64)
        global  CONVERT_ENERGY_AU_TO_EV, CONVERT_ENERGY_AU_TO_KAYSERS, CONVERT_ENERGY_AU_TO_PER_SEC, CONVERT_TIME_AU_TO_SEC,
                CONVERT_CROSS_SECTION_AU_TO_BARN,  CONVERT_RATE_AU_TO_PER_SEC,  CONVERT_LENGTH_AU_TO_FEMTOMETER

        if       sa in ["cross section: from atomic to predefined unit", "cross section: from atomic"]
            if      Defaults.getDefaults("unit: cross section") == "a.u."   return( wa )
            elseif  Defaults.getDefaults("unit: cross section") == "barn"   return( wa * CONVERT_CROSS_SECTION_AU_TO_BARN )
            elseif  Defaults.getDefaults("unit: cross section") == "Mbarn"  return( wa * CONVERT_CROSS_SECTION_AU_TO_BARN * 1.0e-6 )
            else    error("stop a")
            end

        elseif   sa in ["cross section: from atomic to barn"]               return( wa * CONVERT_CROSS_SECTION_AU_TO_BARN )
        elseif   sa in ["cross section: from atomic to Mbarn"]              return( wa * CONVERT_CROSS_SECTION_AU_TO_BARN * 1.0e-6 )

        elseif   sa in ["cross section: from predefined to atomic unit", "cross section: to atomic"]
            if      Defaults.getDefaults("unit: cross section") == "a.u."   return( wa )
            elseif  Defaults.getDefaults("unit: cross section") == "barn"   return( wa / CONVERT_CROSS_SECTION_AU_TO_BARN )
            elseif  Defaults.getDefaults("unit: cross section") == "Mbarn"  return( wa / CONVERT_CROSS_SECTION_AU_TO_BARN * 1.0e6 )
            else    error("stop b")
            end

        elseif   sa in ["Einstein B: from atomic"]                          return( wa )  # Einstein B not yet properly converted.

        elseif   sa in ["energy: from atomic to predefined unit", "energy: from atomic"]
            if      Defaults.getDefaults("unit: energy") == "eV"            return( wa * CONVERT_ENERGY_AU_TO_EV )
            elseif  Defaults.getDefaults("unit: energy") == "Kayser"        return( wa * CONVERT_ENERGY_AU_TO_KAYSERS )
            elseif  Defaults.getDefaults("unit: energy") == "Hartree"       return( wa )
            elseif  Defaults.getDefaults("unit: energy") == "Hz"            return( wa * CONVERT_ENERGY_AU_TO_PER_SEC )
            elseif  Defaults.getDefaults("unit: energy") == "A"             return( 1.0e8 / ( CONVERT_ENERGY_AU_TO_KAYSERS * wa ) )
            else    error("stop c")
            end

        elseif   sa in ["energy: from atomic to eV"]                        return( wa * CONVERT_ENERGY_AU_TO_EV )
        elseif   sa in ["energy: from atomic to Kayser"]                    return( wa * CONVERT_ENERGY_AU_TO_KAYSERS )
        elseif   sa in ["energy: from atomic to Hz"]                        return( wa * CONVERT_ENERGY_AU_TO_PER_SEC )
        elseif   sa in ["energy: from atomic to Angstrom"]                  return( 1.0e8 / ( CONVERT_ENERGY_AU_TO_KAYSERS * wa ) )

        elseif   sa in ["energy: from predefined to atomic unit", "energy: to atomic"]
            if      Defaults.getDefaults("unit: energy") == "eV"            return( wa / CONVERT_ENERGY_AU_TO_EV )
            elseif  Defaults.getDefaults("unit: energy") == "Kayser"        return( wa / CONVERT_ENERGY_AU_TO_KAYSERS )
            elseif  Defaults.getDefaults("unit: energy") == "Hartree"       return( wa )
            elseif  Defaults.getDefaults("unit: energy") == "Hz"            return( wa / CONVERT_ENERGY_AU_TO_PER_SEC )
            elseif  Defaults.getDefaults("unit: energy") == "A"             return( CONVERT_ENERGY_AU_TO_KAYSERS / (1.0e8 * wa) )
            else    error("stop d")
            end

        elseif   sa in ["energy: from eV to atomic"]                        return( wa / CONVERT_ENERGY_AU_TO_EV )
        elseif   sa in ["energy: from wavelength [nm] to atomic"]           return( 1.0e7 / (CONVERT_ENERGY_AU_TO_KAYSERS * wa) )

        elseif   sa in ["intensity: from W/cm^2 to atomic"]                 return( wa / CONVERT_INTENSITY_AU_TO_W_CM2 )
        elseif   sa in ["intensity: from atomic to W/cm^2"]                 return( wa * CONVERT_INTENSITY_AU_TO_W_CM2 )

        elseif    sa in ["rate: from atomic to predefined unit", "rate: from atomic"]
            if       Defaults.getDefaults("unit: rate") == "1/s"            return( wa * CONVERT_RATE_AU_TO_PER_SEC )
            elseif   Defaults.getDefaults("unit: rate") == "a.u."           return( wa )
            else     error("stop e")
            end

        elseif   sa in ["rate: from atomic to 1/s"]                         return( wa * CONVERT_RATE_AU_TO_PER_SEC )

        elseif   sa in ["rate: from predefined to atomic unit", "rate: to atomic"]
            if      Defaults.getDefaults("rate: time") == "1/s"             return( wa / CONVERT_RATE_AU_TO_PER_SEC )
            elseif  Defaults.getDefaults("rate: time") == "a.u."            return( wa  )
            else    error("stop f")
            end

        elseif    sa in ["strength: from atomic to predefined unit", "strength: from atomic"]
            if       Defaults.getDefaults("unit: strength") == "barn eV"    return( wa * CONVERT_STRENGTH_AU_TO_BARN_EV )
            elseif   Defaults.getDefaults("unit: strength") == "cm^2 eV"    return( wa * CONVERT_STRENGTH_AU_TO_CM2_EV )
            elseif   Defaults.getDefaults("unit: strength") == "a.u."       return( wa )
            else     error("stop fx")
            end

        elseif  sa in ["time: from atomic to predefined unit", "time: from atomic"]
            if      Defaults.getDefaults("unit: time") == "sec"             return( wa * CONVERT_TIME_AU_TO_SEC )
            elseif  Defaults.getDefaults("unit: time") == "fs"              return( wa * CONVERT_TIME_AU_TO_SEC * 10.0e15 )
            elseif  Defaults.getDefaults("unit: time") == "as"              return( wa * CONVERT_TIME_AU_TO_SEC * 10.0e18 )
            elseif  Defaults.getDefaults("unit: time") == "a.u."            return( wa )
            else    error("stop g")
            end

        elseif   sa in ["time: from atomic to sec"]                         return( wa * CONVERT_TIME_AU_TO_SEC )
        elseif   sa in ["time: from atomic to fs"]                          return( wa * CONVERT_TIME_AU_TO_SEC * 10.0e15 )
        elseif   sa in ["time: from atomic to as"]                          return( wa * CONVERT_TIME_AU_TO_SEC * 10.0e18 )

        elseif  sa in ["time: from predefined to atomic unit", "time: to atomic"]
            if      Defaults.getDefaults("unit: time") == "sec"             return( wa / CONVERT_TIME_AU_TO_SEC )
            elseif  Defaults.getDefaults("unit: time") == "fs"              return( wa / (CONVERT_TIME_AU_TO_SEC * 10.0e15) )
            elseif  Defaults.getDefaults("unit: time") == "as"              return( wa / (CONVERT_TIME_AU_TO_SEC * 10.0e18) )
            elseif  Defaults.getDefaults("unit: time") == "a.u."            return( wa )
            else    error("stop h")
            end

        elseif   sa in ["temperature: from Kelvin to (Hartree) units"]      return( wa / 315774.64 )
        elseif   sa in ["temperature: from atomic to Kelvin"]               return( wa * 315774.64 )

        elseif  sa in ["length: from fm to atomic"]                         return (wa / CONVERT_LENGTH_AU_TO_FEMTOMETER )
        elseif  sa in ["length: from atomic to fm"]                         return (wa * CONVERT_LENGTH_AU_TO_FEMTOMETER )

        elseif  sa in ["kinetic energy to wave number: atomic units"]
            c = INVERSE_FINE_STRUCTURE_CONSTANT
            wb = sqrt( wa*wa/(c*c) + 2wa );                                 return( wb )

        elseif  sa in ["kinetic energy to wavelength: atomic units"]
            c = INVERSE_FINE_STRUCTURE_CONSTANT
            wb = sqrt( wa*wa/(c*c) + 2wa );                                 return( 2pi / wb )

        elseif  sa in ["wave number to total electron energy: atomic units"]
            c = INVERSE_FINE_STRUCTURE_CONSTANT
            wb = sqrt( wa*wa/(c*c) + 2wa );                                 return( sqrt(wb*wb*c*c + c^4) )

        elseif  sa in ["wave number to kinetic energy: atomic units"]
            c = INVERSE_FINE_STRUCTURE_CONSTANT
            wb = sqrt( wa*wa/(c*c) + 2wa );                                 return( c + sqrt(wb*wb*c*c + c^2) )

        else
            error("Unsupported keystring = $sa")
        end
    end

    """
    `Defaults.setDefaults()`
        ... (re-) defines some 'standard' settings which are common to all the computations with the JAC module, and which can
            be 'overwritten' by the user. --- An improper setting of some variable may lead to an error message, if recognized
            immediately. The following defaults apply if not specified otherwise by the user: the framework is 'relativistic',
            energies are given in eV and cross sections in barn. Note that, internally, atomic units are used throughout for
            all the computations within the program. nothing is returned if not indicated otherwise.

    ## + `("framework: relativistic")`  or  `("framework: non-relativistic")`
    ##     ... to define a relativistic or non-relativistic framework for all subsequent computations.

    + `("method: continuum, spherical Bessel")`  or  `("method: continuum, pure sine")`  or
      `("method: continuum, asymptotic Coulomb")`  or  `("method: continuum, nonrelativistic Coulomb")`  or
      `("method: continuum, Galerkin")`
        ... to define a a method for the generation of the continuum orbitals as (pure) spherical Bessel, pure sine,
            asymptotic Coulomb, nonrelativistic Coulomb orbital or by means of the B-spline-Galerkin method.

    + `("method: normalization, pure sine")`  or  `("method: normalization, pure Coulomb")`  or  `("method: normalization, Ong-Russek")`
        ... to define a method for the normalization of the continuum orbitals as asymptotically (pure) sine or Coulomb
            functions, or following the procedure by Ong & Russek (1978).

    + `("QED model: Petersburg")`  or  `("QED model: Sydney")`
        ... to define a model for the computation of the QED corrections following the work by Shabaev et al. (2011; Petersburg)
            or Flambaum and Ginges (2004; Syney).

    + `("unit: energy", "eV")`  or  `("unit: energy", "Kayser")`  or  `("unit: energy", "Hartree")`  or
      `("unit: energy", "Hz")`  or  `("unit: energy", "Hz")`
        ... to (pre-) define the energy units for all further printouts and communications with the JAC module.

    + `("unit: cross section", "a.u.")`  or  `("unit: cross section", "barn")`  or  `("unit: cross section", "Mbarn")`
        ... to (pre-) define the unit for the printout of cross sections.

    + `("unit: rate", "a.u.")`  or  `("unit: rate", "1/s")`  ... to (pre-) define the unit for the printout of rates.

    + `("unit: resonance strength", "a.u.")`  or  `("unit: resonance strength", "barn eV")`  or
      `("unit: resonance strength", "cm^2 eV")`  ... to (pre-) define the unit for the printout of resonance strengths.

    + `("unit: time", "a.u.")`  or  `("unit: time", "sec")`  or  `("unit: time", "fs")`  or  `("unit: time", "as")`
        ... to (pre-) define the unit for the printout and communications of times with the JAC module.
    """
    function setDefaults(sa::String)
        global GBL_FRAMEWORK, GBL_CONT_SOLUTION, GBL_CONT_NORMALIZATION

        if        sa == "framework: relativistic"                            GBL_FRAMEWORK           = "relativistic"
        elseif    sa == "framework: non-relativistic"                        GBL_FRAMEWORK           = "non-relativistic"
        elseif    sa == "method: continuum, spherical Bessel"                GBL_CONT_SOLUTION       = ContBessel()
        elseif    sa == "method: continuum, pure sine"                       GBL_CONT_SOLUTION       = ContSine()
        elseif    sa == "method: continuum, asymptotic Coulomb"              GBL_CONT_SOLUTION       = AsymptoticCoulomb()
        elseif    sa == "method: continuum, nonrelativistic Coulomb"         GBL_CONT_SOLUTION       = NonrelativisticCoulomb()
        elseif    sa == "method: continuum, Galerkin"                        GBL_CONT_SOLUTION       = BsplineGalerkin()
        elseif    sa == "method: normalization, pure sine"                   GBL_CONT_NORMALIZATION  = PureSineNorm()
        elseif    sa == "method: normalization, pure Coulomb"                GBL_CONT_NORMALIZATION  = CoulombSineNorm()
        elseif    sa == "method: normalization, Ong-Russek"                  GBL_CONT_NORMALIZATION  = OngRussekNorm()
        else      error("Unsupported keystring:: $sa")
        end
        nothing
    end

    function setDefaults(sa::String, sb::String)
        ## global GBL_ENERGY_UNIT, GBL_CROSS_SECTION_UNIT, GBL_RATE_UNIT, GBL_TIME_UNIT, GBL_SUMMARY_IOSTREAM, GBL_PRINT_SUMMARY,
        ##    GBL_TEST_IOSTREAM, GBL_PRINT_TEST

        if        sa == "unit: energy"
            units = ["eV", "Kayser", "Hartree", "Hz", "A"]
            !(sb in units)    &&    error("Currently supported energy units: $(units)")
            global GBL_ENERGY_UNIT = sb
        elseif    sa == "unit: cross section"
            units = ["a.u.", "barn", "Mbarn"]
            !(sb in units)    &&    error("Currently supported cross section units: $(units)")
            global GBL_CROSS_SECTION_UNIT = sb
        elseif    sa == "unit: rate"
            units = ["a.u.", "1/s"]
            !(sb in units)    &&    error("Currently supported rate units: $(units)")
            global GBL_RATE_UNIT = sb
        elseif    sa == "unit: strength"
            units = ["a.u.", "barn eV", "cm^2 eV"]
            !(sb in units)    &&    error("Currently supported resonance strength units: $(units)")
            global GBL_STRENGTH_UNIT = sb
        elseif    sa == "unit: time"
            units = ["a.u.", "sec", "fs", "as"]
            !(sb in units)    &&    error("Currently supported time units: $(units)")
            global GBL_TIME_UNIT = sb
        elseif    sa == "print summary: open"
            global GBL_PRINT_SUMMARY    = true
            global GBL_SUMMARY_IOSTREAM = open(sb, "w")
            println(GBL_SUMMARY_IOSTREAM, "Summary file opened at $( string(now())[1:16] ): \n" *
                                        "========================================   \n")
        elseif    sa == "print summary: append"
            global GBL_PRINT_SUMMARY    = true
            global GBL_SUMMARY_IOSTREAM = open(sb, "a")
            println(GBL_SUMMARY_IOSTREAM, "Summary file re-opened (to append) at $( string(now())[1:16]): \n" *
                                        "======================================================= \n")
        elseif    sa == "print summary: close"
            global GBL_PRINT_SUMMARY    = false
            close(GBL_SUMMARY_IOSTREAM)
        elseif    sa == "print test: open"
            global GBL_PRINT_TEST    = true
            JAC.JAC_TEST_IOSTREAM = open(sb, "w")
            println(JAC.JAC_TEST_IOSTREAM, "Test report file opened at $( string(now())[1:16] ): \n" *
                                    "============================================  \n")
        elseif    sa == "print test: close"
            global GBL_PRINT_TEST    = false
            close(GBL_TEST_IOSTREAM)
        else      error("Unsupported keystring:: $sa")
        end

        nothing
    end

    """
    + `("relativistic subshell list", subshells::Array{Subshell,1}; printout::Bool=true)`
        ... to (pre-) define internally the standard relativistic subshell list on which the standard order of orbitals is based.
    """
    function setDefaults(sa::String, subshells::Array{Subshell,1}; printout::Bool=true)
        if        sa == "relativistic subshell list"
            if printout    println("(Re-) Define a new standard subshell list.")    end
            global GBL_STANDARD_SUBSHELL_LIST = deepcopy(subshells)
        else
            error("Unsupported keystring:: $sa")
        end

        nothing
    end

    """
    + `("standard grid", grid::Radial.Grid; printout::Bool=true)`
        ... to (pre-) define internally the standard radial grid which is used to represent most orbitals.
    """
    function setDefaults(sa::String, grid::Radial.Grid; printout::Bool=true)
        global GBL_STANDARD_GRID

        if        sa == "standard grid"
            if  printout    println("(Re-) Define the standard grid with $(grid.NoPoints) grid points.")    end
            GBL_STANDARD_GRID = grid
        else
        error("Unsupported keystring:: $sa")
        end

        nothing
    end

    """
    + `("QED: damped-hydrogenic", Znuc::Float64, wa::Array{Float64,1})`
        ... to (re-) define the lambda-C damped overlap integrals of the lowest kappa-orbitals
            [ wa_1s_1/2, wa_2p_1/2, wa_2p_3/2, wa_3d_3/2, wa_3d_5/2 ] for the (new) nuclear charge Znuc;
            nothing is returned.
    """
    function setDefaults(sa::String, Znuc::Float64, wa::Array{Float64,1})
        global GBL_QED_HYDROGENIC_LAMBDAC,  GBL_QED_NUCLEAR_CHARGE

        if        sa == "QED: damped-hydrogenic"
            if  length(wa) != 5  error("stop a ")   end
            println("Re-define the damped overlap integrals < a | e^{-r/lambda_C} | a > of the lowest kappa-orbitals for nuclear charge Z = " *
                    "$(Znuc).")
            GBL_QED_NUCLEAR_CHARGE     = Znuc
            GBL_QED_HYDROGENIC_LAMBDAC = wa
        else
        error("Unsupported keystring:: $sa")
        end

        nothing
    end

    """
    `Defaults.getDefaults()`
        ... gives/supplies different information about the (present) framework of the computation or about some
            given data; cf. Defaults.setDefaults().

    + `("alpha")`  or  `("fine-structure constant alpha")`
        ... to get the (current) value::Float64 of the fine-structure constant alpha.

    + `("electron mass: kg")`  or  `("electron mass: amu")`
        ... to get the (current) value::Float64 of the electron mass in the specified unit.

    + `("framework")`  ... to give the (current) setting::String  of the overall framework.

    + `("electron rest energy")`  or  `("mc^2")`  ... to get the electron rest energy.

    + `("electron g-factor")`  ... to give the electron g-factor g_s = 2.00232.

    + `("unit: energy")`  or  `("unit: cross section")`  or  `("unit: rate")`  or  `("unit: strength")`  or  `("unit: time")`
        ... to get the corresponding (user-defined) unit::String for the current computations.

    + `("standard grid")`
        ... to get the (current standard) grid::Array{Float64,1} to which all radial orbital functions usually refer.

    + `("speed of light: c")`  ... to get the speed of light in atomic units.

    + `("summary flag/stream")`
        ... to get the logical flag and stream for printing a summary file; a tupel (flag, iostream) is returned.
    """
    function getDefaults(sa::String)
        global GBL_FRAMEWORK, GBL_CONTINUUM, GBL_ENERGY_UNIT, GBL_CROSS_SECTION_UNIT, GBL_RATE_UNIT, GBL_STRENGTH_UNIT, GBL_TIME_UNIT
        global ELECTRON_MASS_SI, ELECTRON_MASS_U, FINE_STRUCTURE_CONSTANT, GBL_STANDARD_GRID, GBL_SUMMARY_IOSTREAM, GBL_PRINT_SUMMARY,
               GBL_TEST_IOSTREAM, GBL_TEST_SUMMARY

        if        sa in ["alpha", "fine-structure constant alpha"]          return (FINE_STRUCTURE_CONSTANT)
        elseif    sa == "electron mass: kg"                                 return (ELECTRON_MASS_SI)
        elseif    sa == "electron mass: amu"                                return (ELECTRON_MASS_U)
        elseif    sa == "electron g-factor"                                 return (2.00232)
        elseif    sa == "framework"                                         return (GBL_FRAMEWORK)
        elseif    sa == "method: continuum"                                 return (GBL_CONTINUUM)
        elseif    sa in ["electron rest energy", "mc^2"]                    return (1/(FINE_STRUCTURE_CONSTANT^2))
        elseif    sa == "unit: energy"                                      return (GBL_ENERGY_UNIT)
        elseif    sa == "unit: cross section"                               return (GBL_CROSS_SECTION_UNIT)
        elseif    sa == "unit: rate"                                        return (GBL_RATE_UNIT)
        elseif    sa == "unit: strength"                                    return (GBL_STRENGTH_UNIT)
        elseif    sa == "unit: time"                                        return (GBL_TIME_UNIT)
        elseif    sa == "standard grid"                                     return (GBL_STANDARD_GRID)
        elseif    sa == "speed of light: c"                                 return (1.0/FINE_STRUCTURE_CONSTANT)
        elseif    sa == "summary flag/stream"                               return ( (GBL_PRINT_SUMMARY, GBL_SUMMARY_IOSTREAM) )
        elseif    sa == "test flag/stream"                                  return ( (GBL_PRINT_TEST,    JAC.JAC_TEST_IOSTREAM) )
        else      error("Unsupported keystring:: $sa")
        end
    end

    """
    + `("ordered shell list: non-relativistic", n_max::Int64)`
        ... to give an ordered list of non-relativistic shells::Array{Shell,1} up to the (maximum) principal number n_max.

    + `("ordered subshell list: relativistic", n_max::Int64)`
        ... to give an ordered list of relativistic subshells::Array{Subshell,1} up to the (maximum) principal number n_max.
    """
    function getDefaults(sa::String, n_max::Int64)

        !(1 <= n_max < 10)    &&    error("Unsupported value of n_max = $n_max")

        if        sa == "ordered shell list: non-relativistic"
            wa = Shell[]
            for  n = 1:n_max
            for l = 0:n_max-1    push!(wa, Shell(n,l) )    end
            end
        elseif    sa == "ordered subshell list: relativistic"
            wa = Subshell[]
            for  n = 1:n_max
                for l = 0:n_max-1
                    if     l == 0   push!(wa, Subshell(n, -1) )
                    else   j2 = 2l - 1;    kappa = Int64(  (j2+1)/2 );   push!(wa, Subshell(n, kappa) )
                        j2 = 2l + 1;    kappa = Int64( -(j2+1)/2 );   push!(wa, Subshell(n, kappa) )
                    end
                end
            end
        else    error("Unsupported keystring:: $sa")
        end
        return( wa )
    end

    """
    `Defaults.warn()`
        ... deals with warnings that occur during a run and session; it handles the global array GBL_WARNINGS.

    + `(AddWarning, warning::String)`  ... to add warning to the global array GBL_WARNINGS.

    + `(PrintWarnings)`  ... to print all warnings that are currently kept in the global array GBL_WARNINGS.

    + `(ResetWarnings)`  ... to reset the global array GBL_WARNINGS.
    """
    function warn(wa::AbstractWarning)
        global GBL_WARNINGS

        if        wa == PrintWarnings()
            iostream = open("jac-warn.report", "w")
            println(iostream, " ")
            println(iostream, "\n\n",
                            "Warnings from the present sessions or run, printed at $( string(now())[1:16] ): \n",
                            "======================================================================= \n")
            for sa in  GBL_WARNINGS
                println(iostream, "++ " * sa)
            end
            close(iostream)
            #
        elseif    wa == ResetWarnings()
            ## @warn("Reset global array GBL_WARNINGS.")
            printstyled("Constants.warn():  Reset global array GBL_WARNINGS.", color=:light_magenta)
            GBL_WARNINGS = String[]
        else      error("Unsupported Warnings:: $wa")
        end

        nothing
    end

    function warn(wa::AbstractWarning, sa::String)
        global GBL_WARNINGS

        if        wa == AddWarning()
            push!(GBL_WARNINGS, sa)
        else      error("Unsupported Warnings:: $wa")
        end

        nothing
    end

    """
    `Defaults.saRatip()`
        ... returns true or false if the spin-angular coefficients should be calculated by Fortran codes from RATIP/Grasp.
    """
    saRatip() = false
    saGG()    = true

end # module

ext: jl


text:
using TidalFluxConfigurations
using Base.Test


text:
function is_valid_res_repo_name(res_name::AbstractString)
    try validate_res_name(res_name)
        return true;
    catch return false; end
end

function is_valid_res_dir(res_dir::AbstractString)
    try validate_res_dir(res_dir)
        return true;
    catch return false; end
end

function is_valid_res_repo_tree(rr::ResRepoTree)
    try validate_res_repo_tree(rr)
        return true;
    catch return false; end
end
is_valid_res_repo_tree(res_dir::AbstractString) =
    is_valid_res_repo_tree(ResRepoTree(res_dir))

function is_valid_time_tag(time_tag)
    try validate_time_tag(time_tag)
        return true;
    catch return false; end
end

repo_exist(rr::ResRepoTree) =
    is_valid_res_repo_tree(rr);

has_versions(rr::ResRepoTree) =
    length(find_versions(rr.parent_dir, rr.res_name)) > 0;

has_only_one_version(rr::ResRepoTree) =
    length(find_versions(rr.parent_dir, rr.res_name)) == 1;

# titi delicioso te quiero un poquito mas si me dices que me rio bonito mayor igual que Daniela.
# Esos serian los paraolimpicos de la risa, asi q mira a ver que vas a decir.....


text:
export differentiate

#################################################################
#
# differentiate()
#   based on John's differentiate and this code, I think by Miles Lubin:
#     https://github.com/IainNZ/NLTester/blob/master/julia/nlp.jl#L74
#
#################################################################

differentiate(ex::SymbolicVariable, wrt::SymbolicVariable) = (ex == wrt) ? 1 : 0

differentiate(ex::Number, wrt::SymbolicVariable) = 0

function differentiate(ex::Expr,wrt)
    if ex.head != :call
        error("Unrecognized expression $ex")
    end
    simplify(differentiate(SymbolParameter(ex.args[1]), ex.args[2:end], wrt))
end

differentiate{T}(x::SymbolParameter{T}, args, wrt) = error("Derivative of function " * string(T) * " not supported")

# The Power Rule:
function differentiate(::SymbolParameter{:^}, args, wrt)
    x = args[1]
    y = args[2]
    xp = differentiate(x, wrt)
    yp = differentiate(y, wrt)
    if xp == 0 && yp == 0
        return 0
    elseif yp == 0
        return :( $y * $xp * ($x ^ ($y - 1)) )
    else
        return :( $x ^ $y * ($xp * $y / $x + $yp * log($x)) )
    end
end

function differentiate(::SymbolParameter{:+}, args, wrt)
    termdiffs = Any[:+]
    for y in args
        x = differentiate(y, wrt)
        if x != 0
            push!(termdiffs, x)
        end
    end
    if (length(termdiffs) == 1)
        return 0
    elseif (length(termdiffs) == 2)
        return termdiffs[2]
    else
        return Expr(:call, termdiffs...)
    end
end

function differentiate(::SymbolParameter{:-}, args, wrt)
    termdiffs = Any[:-]
    # first term is special, can't be dropped
    term1 = differentiate(args[1], wrt)
    push!(termdiffs, term1)
    for y in args[2:end]
        x = differentiate(y, wrt)
        if x != 0
            push!(termdiffs, x)
        end
    end
    if term1 != 0 && length(termdiffs) == 2 && length(args) >= 2
        # if all of the terms but the first disappeared, we just return the first
        return term1
    elseif (term1 == 0 && length(termdiffs) == 2)
        return 0
    else
        return Expr(:call, termdiffs...)
    end
end

# The Product Rule
# d/dx (f * g) = (d/dx f) * g + f * (d/dx g)
# d/dx (f * g * h) = (d/dx f) * g * h + f * (d/dx g) * h + ...
function differentiate(::SymbolParameter{:*}, args, wrt)
    n = length(args)
    res_args = Array(Any, n)
    for i in 1:n
       new_args = Array(Any, n)
       for j in 1:n
           if j == i
               new_args[j] = differentiate(args[j], wrt)
           else
               new_args[j] = args[j]
           end
       end
       res_args[i] = Expr(:call, :*, new_args...)
    end
    return Expr(:call, :+, res_args...)
end

# The Quotient Rule
# d/dx (f / g) = ((d/dx f) * g - f * (d/dx g)) / g^2
function differentiate(::SymbolParameter{:/}, args, wrt)
    x = args[1]
    y = args[2]
    xp = differentiate(x, wrt)
    yp = differentiate(y, wrt)
    if xp == 0 && yp == 0
        return 0
    elseif xp == 0
        return :( -$yp * $x / $y^2 )
    elseif yp == 0
        return :( $xp / $y )
    else
        return :( ($xp * $y - $x * $yp) / $y^2 )
    end
end

symbolic_derivative_1arg_list = [
    ( :sqrt,        :(  1 / 2 / sqrt(x)                         ))
    ( :cbrt,        :(  1 / 3 / cbrt(x)^2                       ))
    ( :abs2,        :(  1 * 2 * x                               ))
    ( :inv,         :( -1 * abs2(inv(x))                        ))
    ( :log,         :(  1 / x                                   ))
    ( :log10,       :(  1 / x / log(10)                         ))
    ( :log2,        :(  1 / x / log(2)                          ))
    ( :log1p,       :(  1 / (x + 1)                             ))
    ( :exp,         :(  exp(x)                                  ))
    ( :exp2,        :(  log(2) * exp2(x)                        ))
    ( :expm1,       :(  exp(x)                                  ))
    ( :sin,         :(  cos(x)                                  ))
    ( :cos,         :( -sin(x)                                  ))
    ( :tan,         :(  (1 + tan(x)^2)                          ))
    ( :sec,         :(  sec(x) * tan(x)                         ))
    ( :csc,         :( -csc(x) * cot(x)                         ))
    ( :cot,         :( -(1 + cot(x)^2)                          ))
    ( :sind,        :(  pi / 180 * cosd(x)                      ))
    ( :cosd,        :( -pi / 180 * sind(x)                      ))
    ( :tand,        :(  pi / 180 * (1 + tand(x)^2)              ))
    ( :secd,        :(  pi / 180 * secd(x) * tand(x)            ))
    ( :cscd,        :( -pi / 180 * cscd(x) * cotd(x)            ))
    ( :cotd,        :( -pi / 180 * (1 + cotd(x)^2)              ))
    ( :asin,        :(  1 / sqrt(1 - x^2)                       ))
    ( :acos,        :( -1 / sqrt(1 - x^2)                       ))
    ( :atan,        :(  1 / (1 + x^2)                           ))
    ( :asec,        :(  1 / abs(x) / sqrt(x^2 - 1)              ))
    ( :acsc,        :( -1 / abs(x) / sqrt(x^2 - 1)              ))
    ( :acot,        :( -1 / (1 + x^2)                           ))
    ( :asind,       :(  180 / pi / sqrt(1 - x^2)                ))
    ( :acosd,       :( -180 / pi / sqrt(1 - x^2)                ))
    ( :atand,       :(  180 / pi / (1 + x^2)                    ))
    ( :asecd,       :(  180 / pi / abs(x) / sqrt(x^2 - 1)       ))
    ( :acscd,       :( -180 / pi / abs(x) / sqrt(x^2 - 1)       ))
    ( :acotd,       :( -180 / pi / (1 + x^2)                    ))
    ( :sinh,        :(  cosh(x)                                 ))
    ( :cosh,        :(  sinh(x)                                 ))
    ( :tanh,        :(  sech(x)^2                               ))
    ( :sech,        :( -tanh(x) * sech(x)                       ))
    ( :csch,        :( -coth(x) * csch(x)                       ))
    ( :coth,        :( -csch(x)^2                               ))
    ( :asinh,       :(  1 / sqrt(x^2 + 1)                       ))
    ( :acosh,       :(  1 / sqrt(x^2 - 1)                       ))
    ( :atanh,       :(  1 / (1 - x^2)                           ))
    ( :asech,       :( -1 / x / sqrt(1 - x^2)                   ))
    ( :acsch,       :( -1 / abs(x) / sqrt(1 + x^2)              ))
    ( :acoth,       :(  1 / (1 - x^2)                           ))
    ( :deg2rad,     :(  pi / 180                                ))
    ( :rad2deg,     :(  180 / pi                                ))
    ( :erf,         :(  2 * exp(-x*x) / sqrt(pi)                ))
    ( :erfinv,      :(  0.5 * sqrt(pi) * exp(erfinv(x) * erfinv(x))  ))
    ( :erfc,        :( -2 * exp(-x*x) / sqrt(pi)                ))
    ( :erfcinv,     :( -0.5 * sqrt(pi) * exp(erfcinv(x) * erfcinv(x))  ))
    ( :erfi,        :(  2 * exp(x*x) / sqrt(pi)                 ))
    ( :gamma,       :(  digamma(x) * gamma(x)                   ))
    ( :lgamma,      :(  digamma(x)                              ))
    ( :digamma,     :(  trigamma(x)                             ))
    ( :invdigamma,  :(  inv(trigamma(invdigamma(x)))            ))
    ( :trigamma,    :(  polygamma(2, x)                         ))
    ( :airy,        :(  airyprime(x)                            ))  # note: only covers the 1-arg version
    ( :airyprime,   :(  x * airy(x)                             ))
    ( :airyai,      :(  airyaiprime(x)                          ))
    ( :airybi,      :(  airybiprime(x)                          ))
    ( :airyaiprime, :(  x * airyai(x)                           ))
    ( :airybiprime, :(  x * airybi(x)                           ))
    ( :besselj0,    :( -besselj1(x)                             ))
    ( :besselj1,    :(  (besselj0(x) - besselj(2, x)) / 2       ))
    ( :bessely0,    :( -bessely1(x)                             ))
    ( :bessely1,    :(  (bessely0(x) - bessely(2, x)) / 2       ))
    ( :erfcx,       :(  (2 * x * erfcx(x) - 2 / sqrt(pi))       ))
    ( :dawson,      :(  (1 - 2x * dawson(x))                    ))

]

# This is the public interface for accessing the list of symbolic
# derivatives. The format is a list of (Symbol,Expr) tuples
# (:f, deriv_expr), where deriv_expr is a symbolic
# expression for the first derivative of the function f.
# The symbol :x is used within deriv_expr for the point at
# which the derivative should be evaluated.
symbolic_derivatives_1arg() = symbolic_derivative_1arg_list
export symbolic_derivatives_1arg

# deprecated: for backward compatibility with packages that used
# this unexported interface.
derivative_rules = Array(@Compat.compat(Tuple{Symbol,Expr}),0)
for (s,ex) in symbolic_derivative_1arg_list
    push!(derivative_rules, (s, :(xp*$ex)))
end

for (funsym, exp) in symbolic_derivative_1arg_list
    @eval function differentiate(::SymbolParameter{$(Meta.quot(funsym))}, args, wrt)
        x = args[1]
        xp = differentiate(x, wrt)
        if xp != 0
            return @sexpr(xp*$exp)
        else
            return 0
        end
    end
end

derivative_rules_bessel = [
    ( :besselj,    :(    (besselj(nu - 1, x) - besselj(nu + 1, x)) / 2   ))
    ( :besseli,    :(    (besseli(nu - 1, x) + besseli(nu + 1, x)) / 2   ))
    ( :bessely,    :(    (bessely(nu - 1, x) - bessely(nu + 1, x)) / 2   ))
    ( :besselk,    :( -1 * (besselk(nu - 1, x) + besselk(nu + 1, x)) / 2   ))
    ( :hankelh1,   :(    (hankelh1(nu - 1, x) - hankelh1(nu + 1, x)) / 2 ))
    ( :hankelh2,   :(    (hankelh2(nu - 1, x) - hankelh2(nu + 1, x)) / 2 ))
]

# This is the public interface for accessing the list of symbolic
# derivatives. The format is a list of (Symbol,Expr) tuples
# (:f, deriv_expr), where deriv_expr is a symbolic
# expression for the first derivative of the function f with respect to x.
# The symbol :nu and :x are used within deriv_expr
# :nu specifies the first parameter of the bessel
# function (usually written n or alpha)
# :x gives the point at which the derivative should be evaluated.
symbolic_derivative_bessel_list() = derivative_rules_bessel
export symbolic_derivative_bessel_list

# 2-argument bessel functions
for (funsym, exp) in derivative_rules_bessel
    @eval function differentiate(::SymbolParameter{$(Meta.quot(funsym))}, args, wrt)
        nu = args[1]
        x = args[2]
        xp = differentiate(x, wrt)
        if xp != 0
            return @sexpr(xp*$exp)
        else
            return 0
        end
    end
end

### Other functions from julia/base/math.jl we might want to define
### derivatives for. Some have two arguments.

## atan2
## hypot
## beta, lbeta, eta, zeta, digamma

## Differentiate for piecewise functions defined using ifelse
function differentiate(::SymbolParameter{:ifelse}, args, wrt)
    :(ifelse($(args[1]), $(differentiate(args[2],wrt)),$(differentiate(args[3],wrt))))
end

function differentiate(ex::Expr, targets::Vector{Symbol})
    n = length(targets)
    exprs = Array(Any, n)
    for i in 1:n
        exprs[i] = differentiate(ex, targets[i])
    end
    return exprs
end

differentiate(ex::Expr) = differentiate(ex, :x)
differentiate(s::Compat.AbstractString, target...) = differentiate(parse(s), target...)
differentiate(s::Compat.AbstractString, target::Compat.AbstractString) = differentiate(parse(s), symbol(target))
differentiate{T <: Compat.AbstractString}(s::Compat.AbstractString, targets::Vector{T}) = differentiate(parse(s), map(symbol, targets))


text:
@testset "Conversion from and to MATLAB model" begin
    filename = download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.mat",
        joinpath("data", "iJO1366.mat"),
        "b5cfe21b6369a00e45d600b783f89521f5cc953e25ee52c5f1d0a3f83743be30",
    )

    mm = load_mat_model(filename)
    sm = convert(StandardModel, mm)
    mm2 = convert(MATModel, sm)

    @test Set(reactions(mm)) == Set(reactions(sm))
    @test Set(reactions(mm)) == Set(reactions(mm2))
end


text:
using Dagger
import Dagger:DomainBlocks

import Dagger: DArray, chunks

function onchunks(X::Indexing)
    # TODO: handle IndexConsts
    idx = X.idx
    Indexing(map(delayed(x -> Indexing(x, idx)), chunks(X.array)), idx)
end

function onchunks(X::Map)
    f = X.f
    Map(delayed((x...) -> Map(f, x)), map(onchunks, X.arrays))
end

function arrayop!{D<:DArray}(::Type{D}, t::Assign)
    lhs = onchunks(t.lhs)
    rhs = onchunks(t.rhs)

    idx = t.lhs.idx
    combine_chunks = delayed() do a, b
        Map(t.reducefn, (a, b))
    end
    ct = Assign(lhs, rhs, combine_chunks, delayed(()->nothing))
    cs = arrayop_treereduce(ct)
    f = delayed() do l, r
        arrayop!(Assign(l, r, t.reducefn, t.empty))
    end
    t.lhs.array.chunks = map(f, onchunks(t.lhs).array, cs)
    t.lhs.array
end

@generated function arrayop_treereduce{L,R,F,E}(op::Assign{L,R,F,E})

    rspaces = index_spaces(:(op.rhs), R)
    lspaces = index_spaces(:(op.lhs), L)
    expr = kernel_expr(:(op.rhs), R)

    reduceddims = setdiff(keys(rspaces), keys(lspaces))
    for sym in reduceddims
        spaces = rspaces[sym]
        T,dim,nm = first(spaces)
        valx = gensym("valx")
        valy = gensym("valy")
        expr = quote
            Dagger.treereduce(indices($nm, $dim)) do x, y
                x, y
                $sym = x

                $valx = $expr
                $sym = y
                $valy = $expr
                op.reducefn($valx, $valy)
            end
        end
    end

    lhs_expr = kernel_expr(:(op.lhs), L) # :() will be ignored
    expr = :($lhs_expr = $expr)

    checks = :()
    for (sym, spaces) in lspaces
        if length(spaces) > 1
            # check dimensions for equality
            equal_dims = [:(size($(d[3]), $(d[2]))) for d in spaces]
            checks = :($checks; @assert allequal($(equal_dims...),))
        end
        T,dim,nm = first(spaces)
        expr = quote
            for $sym = indices($nm, $dim)
                $expr
            end
        end
    end

    :($checks; $expr; op.lhs.array)
end

function Base.indices(x::DArray)
    Dagger.domainchunks(x)
end

function Base.indices(x::DArray, i)
    idxs = indices(x)
    Dagger.DomainBlocks((idxs.start[i],), (idxs.cumlength[i],))
end

function allocarray{T,N}(::Type{DArray{T,N}}, default, idxs...)
    dmnchunks = DomainBlocks(map(i->1, idxs),
                             map(i->isa(i, DomainBlocks) ?
                                 i.cumlength[1] : (length(i),), idxs))

    chnks = map(delayed(subd -> allocarray(Array{T,1}, default, size(subd))), dmnchunks)
    sz = map((x,y)->x-y+1, map(last, dmnchunks.cumlength), dmnchunks.start)
    dmn = ArrayDomain(map(x->1:x, sz))
    DArray{T, length(idxs)}(dmn, dmnchunks, chnks)
end


text:
using WAV
using ChirpSignal

fs = 16000.0; # sampling frequency
fl = 500.0;   # lower bound frequency
fh = 8000.0;  # upper bound frequency
T  = 5.0 ;    # total time of chirp signal

wav = chirp(T, fs, fl, fh; method="linear");
wavwrite(wav, "./doc/chirpLinear.wav", Fs=fs, nbits=32)

wav = chirp(T, fs, fl, fh; method="quadratic");
wavwrite(wav, "./doc/chirpQuadratic.wav", Fs=fs, nbits=32)

wav = chirp(T, fs, fl, fh; method="logarithmic");
wavwrite(wav, "./doc/chirpLogarithmic.wav", Fs=fs, nbits=32)

wav = chirp(T, fs, fl, fh; method="exponential");
wavwrite(wav, "./doc/chirpExponential.wav", Fs=fs, nbits=32)

fun(t) = 2000*t
wav = chirp(2.0, 16000.0, fun)
wavwrite(wav, "./doc/0_4000hz_linear.wav", Fs=fs, nbits=32)

wav = chirp(2.0, 16000.0, t -> 2000*t^2)
wavwrite(wav, "./doc/0_8000hz_quadratic.wav", Fs=fs, nbits=32)

max_stars_repo_name: sonosole/ChirpSignal

max_issues_repo_name: sonosole/ChirpSignal

max_forks_repo_name: sonosole/ChirpSignal


text:
module Parquet
using StaticArrays, PyCall
using AbstractTrees
using Parameters, Combinatorics
using DataFrames
using ..DiagTree

const DI, EX, BOTH = 1, 2, 3
const INL, OUTL, INR, OUTR = 1, 2, 3, 4
# orginal diagrams T, U, S; particle-hole counterterm Ts, Us; and their counterterm Tc, Uc, Sc, Tsc, Usc
# symmetry factor for Alli, PHr, PHEr, PPr, PHrc, PHErc
const SymFactor = [1.0, -1.0, 1.0, -0.5, +1.0, -1.0]

include("common.jl")
export ParquetBlocks

include("filter.jl")
include("vertex4.jl")

include("sigma.jl")
include("green.jl")
include("vertex3.jl")
include("polarization.jl")

include("benchmark/benchmark.jl")
end