text:
#
# TriangulateIO structure and public interface methods
#
# Public interface methods have docstrings, private ones
# have comments.

#########################################################
"""
$(TYPEDEF)

Julia version of Triangle's triangulateio structure.

Used to  pass data into and  out of the triangulate()  procedure.  The
arrays are stored  in column major order with the first index being
the local number in a triangle or segment and the second index being the
index in the element count. This exactly corresponds
to  the layout  used in  Triangle. This  means that  all data  in this
structure are passed from/to triangle without copying.

Arrays are used to store points, triangles, markers, and so forth.

Description of fields:

$(TYPEDFIELDS)

"""
mutable struct TriangulateIO

    """
    An array of point coordinates with `size(pointlist,1)==2`.
    `pointlist' must always point to a list of points. __Mandatory__.
    """
    pointlist :: Array{Cdouble,2}

    """
    An array of point attributes. There can be several attributes per point.
    Optional for input.
    """
    pointattributelist :: Array{Cdouble,2}

    """
    An array of point markers. Optional for input.
    """
    pointmarkerlist :: Array{Cint,1}

    """
    An array of triangle corners. The first three entries
    of each column describe the three nodes of the triangle
    in counterclockwise manner. They are followd by  any other nodes if the triangle
    represents a nonlinear element. This means that without specfiying the 'o2' switch,
    each column contains three node indices, while in the opposite case, each column
    contains six entries: the three triangle vertices and the three edge mindpoint nodes.

    __Mandatory if the 'r' switch is used__. In this case `trianglelist'
    must point to a list of triangles with optional higher order nodes.
    """
    trianglelist :: Array{Cint,2}

    """
    An array of triangle attributes. There can be several attributes per triangle.
    Optional on input.
    """
    triangleattributelist :: Array{Cdouble,2}

    """
    An array of triangle area constraints. Input only.

    __Mandatory if both the 'r' and the 'a'  switch (with no number following) are used.__
    """
    trianglearealist :: Array{Cdouble,1}

    """
    An array of triangle neighbors. `size(neighborlist,1)==3`
    triangle.  Output only.
    """
    neighborlist :: Array{Cint,2}

    """
    An array of segment endpoints. `size(segmentlist,1)==2`

    __Mandatory if the 'p' switch is used__.
    """
    segmentlist :: Array{Cint,2}

    """
    An array of segment markers. Optional on input. If not set then
    segment markers on output default to zero.
    """
    segmentmarkerlist :: Array{Cint,1}

    """
    An array of holes. Holes are marked by some point from within the hole.
    Input only, although the array pointer is copied to the
    output structure for convenience.

    Used if  the 'p' switch is used without the 'r' switch.
    """
    holelist :: Array{Cdouble,2}

    """
    An array of regional attributes and area constraints.

    Used if  the 'p' switch is used without the 'r' switch.

    Each of the columns of this array contains
    the constraint's x and y coordinates are at indices [1] and [2],
    followed by the regional attribute at index [3], followed by the
    maximum area at index [4]. So we have `size(regionlist,1)==4`.
    Note that each regional attribute is
    used only if you select the 'A' switch, and each area constraint is
    used only if you select the 'a' switch (with no number following), but
    omitting one of these switches does not change the memory layout.
    Input only, although the pointer is copied to the output structure for
    convenience.
    """
    regionlist ::  Array{Cdouble,2}

    """
    An array of edge endpoints. `sizeof(edgelist,1)==2`.  Output only.
    """
    edgelist :: Array{Cint,2}

    """
    An array of edge markers; Output   only.
    """
    edgemarkerlist :: Array{Cint,1}

    """
    An array of normal vectors, used for infinite rays in
    Voronoi diagrams. For each finite edge in a Voronoi diagram,
    the normal vector written is the
     zero vector.  `sizeof(normlist,1)==2`. Output only.
    """
    normlist :: Array{Cdouble,2}
end

##########################################################
"""
$(TYPEDSIGNATURES)
Return number of points in triangulatio structure.
"""
numberofpoints(tio::TriangulateIO)=size(tio.pointlist,2)

##########################################################
"""
$(TYPEDSIGNATURES)

Return number of segments in triangulateio structure.
"""
numberofsegments(tio::TriangulateIO)=size(tio.segmentlist,2)

##########################################################
"""
$(TYPEDSIGNATURES)

Return number of triangles in triangulateio structure.
"""
numberoftriangles(tio::TriangulateIO)=size(tio.trianglelist,2)

##########################################################
"""
$(TYPEDSIGNATURES)

Return number of triangles in triangulateio structure.
"""
numberofedges(tio::TriangulateIO)=size(tio.edgelist,2)

##########################################################
"""
$(TYPEDSIGNATURES)

Return number of regions in triangulateio structure.
"""
numberofregions(tio::TriangulateIO)=size(tio.regionlist,2)

##########################################################
"""
$(TYPEDSIGNATURES)

Return number of holes in triangulateio structure.
"""
numberofholes(tio::TriangulateIO)=size(tio.holelist,2)

##########################################################
"""
$(TYPEDSIGNATURES)

Return number of holes in triangulateio structure.
"""
numberoftriangleattributes(tio::TriangulateIO)=size(tio.triangleattributelist,1)

##########################################################
"""
$(TYPEDSIGNATURES)

Create TriangulateIO structure with empty data.
"""
function TriangulateIO(;
                       pointlist=Array{Cdouble,2}(undef,0,0),
                       pointattributelist=Array{Cdouble,2}(undef,0,0),
                       pointmarkerlist=Array{Cint,1}(undef,0),
                       trianglelist=Array{Cint,2}(undef,0,0),
                       triangleattributelist=Array{Cint,2}(undef,0,0),
                       trianglearealist=Array{Cdouble,1}(undef,0),
                       neighborlist=Array{Cint,2}(undef,0,0),
                       segmentlist=Array{Cint,2}(undef,0,0),
                       segmentmarkerlist=Array{Cint,1}(undef,0),
                       holelist=Array{Cdouble,2}(undef,0,0),
                       regionlist=Array{Cdouble,2}(undef,0,0),
                       edgelist=Array{Cint,2}(undef,0,0),
                       edgemarkerlist=Array{Cint,1}(undef,0),
                       normlist=Array{Cdouble,2}(undef,0,0))

    return TriangulateIO(pointlist,
                         pointattributelist,
                         pointmarkerlist,
                         trianglelist,
                         triangleattributelist,
                         trianglearealist,
                         neighborlist,
                         segmentlist,
                         segmentmarkerlist,
                         holelist,
                         regionlist,
                         edgelist,
                         edgemarkerlist,
                         normlist)
end

function Base.show(io::IO, tio::TriangulateIO)
    nonempty(a)=size(a,ndims(a))>0
    println(io,"TriangulateIO(")
    for name in fieldnames(typeof(tio))
        a=getfield(tio,name)
        if nonempty(a)
            print(io,"$(name)=")
            print(io,a)
            println(io,",")
        end
    end
    println(io,")")
end

#
# Construct CTriangulateIO from TriangulateIO
#
function CTriangulateIO(tio::TriangulateIO)
    ctio=CTriangulateIO()
    ctio.numberofpoints=size(tio.pointlist,2)
    @assert ctio.numberofpoints>0
    @assert  size(tio.pointlist,1)==2
    ctio.pointlist=pointer(tio.pointlist)

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

    if size(tio.pointmarkerlist,1)>0
        @assert size(tio.pointmarkerlist,1)==ctio.numberofpoints
        ctio.pointmarkerlist=pointer(tio.pointmarkerlist)
    end

    ctio.numberoftriangles=size(tio.trianglelist,2)
    if ctio.numberoftriangles>0
        ctio.trianglelist=pointer(tio.trianglelist)
        ctio.numberoftriangleattributes=size(tio.triangleattributelist,1)
        ctio.numberofcorners=size(tio.trianglelist,1) # JRS: "Each triangle occupies   `numberofcorners' ints."
        if ctio.numberoftriangleattributes>0
            @assert size(tio.triangleattributelist,2)==ctio.numberoftriangles
            ctio.triangleattributes=pointer(tio.triangleattributes)
        end
        if size(tio.trianglearealist,1)>0
            @assert size(tio.trianglearealist,1)==ctio.numberoftriangles
            ctio.trianglearealist=pointer(tio.trianglearealist)
        end
    end

    ctio.numberofsegments=size(tio.segmentlist,2)
    if ctio.numberofsegments>0
        ctio.segmentlist=pointer(tio.segmentlist)

        if size(tio.segmentmarkerlist, 1) > 0
            @assert size(tio.segmentmarkerlist, 1) == ctio.numberofsegments
            ctio.segmentmarkerlist = pointer(tio.segmentmarkerlist)
        end
    end

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

    ctio.numberofregions=size(tio.regionlist,2)
    if ctio.numberofregions>0
        @assert size(tio.regionlist,1)==4
        ctio.regionlist=pointer(tio.regionlist)
    end
    return ctio
end

#
# Construct TriangulateIO from CTriangulateIO
#
function TriangulateIO(ctio::CTriangulateIO)
    tio=TriangulateIO()

    if ctio.numberofpoints>0 && ctio.pointlist!=C_NULL
        tio.pointlist = convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio.pointlist, (2,Int(ctio.numberofpoints)), own=true))
    end
    if ctio.numberofpointattributes>0  && ctio.pointattributelist!=C_NULL
        tio.pointattributelist=convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio.pointattributelist, (Int(ctio.numberofpointattributes),Int(ctio.numberofpoints)), own=true))
    end

    if ctio.pointmarkerlist!=C_NULL
        tio.pointmarkerlist=convert(Array{Cint,1}, Base.unsafe_wrap(Array, ctio.pointmarkerlist, (Int(ctio.numberofpoints)), own=true))
    end

    if ctio.numberoftriangles>0 && ctio.trianglelist!=C_NULL
        tio.trianglelist=convert(Array{Cint,2}, Base.unsafe_wrap(Array, ctio.trianglelist, (Int(ctio.numberofcorners),Int(ctio.numberoftriangles)), own=true))
    end

    if ctio.numberoftriangleattributes>0  && ctio.triangleattributelist!=C_NULL
        tio.triangleattributelist=convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio.triangleattributelist, (Int(ctio.numberoftriangleattributes),Int(ctio.numberoftriangles)), own=true))
    end

    if ctio.trianglearealist!=C_NULL
        tio.trianglearealist=convert(Array{Cdouble,1}, Base.unsafe_wrap(Array, ctio.trianglearealist, (Int(ctio.numberoftriangles)), own=true))
    end

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

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

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

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

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

##########################################################
"""
$(TYPEDEF)

Exception thrown if triangulate() encounters an error
"""
struct TriangulateError <: Exception
end

##########################################################
"""
$(TYPEDSIGNATURES)

Create triangulation. Returns tuple `(out::TriangulateIO, vor_out::TriangulateIO)`
containing the output triangulation and the optional Voronoi tesselation.

After a call to triangulate(), the valid fields of `out` and `vorout`
will depend, in an obvious way, on the choice of switches used.  Note
that when the 'p' switch is used, the pointers `holelist` and
`regionlist` are copied from `tri_in` to the output, but no new space is
allocated;  On   the other hand, Triangle will never copy the `pointlist` pointer (or any
others); new space is allocated for `out.pointlist`, or if the 'N'
switch is used, `out.pointlist` remains uninitialized.

This is the list of switches used by triangle:

| Switch | Meaning                                                       |
|--------|:--------------------------------------------------------------|
| -p     | Triangulates a Planar Straight Line Graph (.poly file).       |
| -r     | Refines a previously generated mesh.                          |
| -q     | Quality mesh generation.  A minimum angle may be specified.   |
| -a     | Applies a maximum triangle area constraint.                   |
| -u     | Applies a user-defined triangle constraint.                   |
| -A     | Applies attributes to identify triangles in certain regions.  |
| -c     | Encloses the convex hull with segments.                       |
| -D     | Conforming Delaunay:  all triangles are truly Delaunay.       |
| -j     | Jettison unused vertices from output .node file.              |
| -e     | Generates an edge list.                                       |
| -v     | Generates a Voronoi diagram.                                  |
| -n     | Generates a list of triangle neighbors.                       |
| -g     | Generates an .off file for Geomview.                          |
| -B     | Suppresses output of boundary information.                    |
| -P     | Suppresses output of .poly file.                              |
| -N     | Suppresses output of .node file.                              |
| -E     | Suppresses output of .ele file.                               |
| -I     | Suppresses mesh iteration numbers.                            |
| -O     | Ignores holes in .poly file.                                  |
| -X     | Suppresses use of exact arithmetic.                           |
| -z     | Numbers all items starting from zero (rather than one).       |
| -o2    |  Generates second-order subparametric elements.               |
| -Y     | Suppresses boundary segment splitting.                        |
| -S     | Specifies maximum number of added Steiner points.             |
| -i     | Uses incremental method, rather than divide-and-conquer.      |
| -F     | Uses Fortune's sweepline algorithm, rather than d-and-c.      |
| -l     | Uses vertical cuts only, rather than alternating cuts.        |
| -s     | Force segments into mesh by splitting (instead of using CDT). |
| -C     | Check consistency of final mesh.                              |
| -Q     | Quiet:  No terminal output except errors.                     |
| -V     | Verbose:  Detailed information on what I'm doing.             |

"""
function triangulate(switches::String, tri_in::TriangulateIO)
    ctio_in=CTriangulateIO(tri_in)
    ctio_out=CTriangulateIO()
    cvor_out=CTriangulateIO()
    rc=triangulate(switches,ctio_in,ctio_out,cvor_out)
    if rc!=0
        throw(TriangulateError())
    end
    out=TriangulateIO(ctio_out)
    vor_out=TriangulateIO(cvor_out)
    return out,vor_out
end

##########################################################
"""
$(TYPEDSIGNATURES)

Set triunsuitable callback used by Triangle if the '-u' flag is set.

This is a function called by Triangle with the coordinates of the vertices
of a triangle in order to learn if that triangle needs further
refinement (i.e. 'true' returned) or not ('false' returned).

Other checks (e.g. maximum edge lengths) are possible here as well.

Note, that the handling of this function is currently not thread safe.
````
function unsuitable(x1,y1,x2,y2,x3,y3,area)
    myarea=locally_desired_area(x1,y1,x2,y2,x3,y3)
    if area>myarea
       return true
    else
       return false
    end
end
````

"""
function triunsuitable(unsuitable::Function;check_signature=true)
    global triunsuitable_func
    if check_signature
        unsuitable(1,2,3,4,5,6,7)
    end
    triunsuitable_func=unsuitable
    nothing
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7fb69cf0-fc66-11eb-256f-cb6e9d2db156
using PlutoUI, HypertextLiteral, LightGraphs

# ╔═╡ 20ad770b-34a1-49e7-b948-e12f46bde5ee
using HypertextLiteral: JavaScript

# ╔═╡ 59cef39a-3577-4726-aeff-960ad5b580c5
tprint(input) = with_terminal() do
	println(input)
end

# ╔═╡ c8ba27b1-4d2b-4afc-89a9-806be91534dc
TableOfContents()

# ╔═╡ ef077759-88ad-4655-b04f-83c80270d320
md"""
## Functions and Constants
"""

# ╔═╡ a626d3c1-ba41-4b91-a2a2-6ff69cd13bce
const BTN_STYLE = """
.btn-group, .btn-group-vertical {
    position: relative;
    display: inline-block;
    vertical-align: middle;
}

.btn-group button {
	background-color: #04AA6D;
	border: 1px solid green;
	color: white;
	padding: 10px 10px;
	cursor: pointer;
	float: left;
	border-radius: 10px;
}

.btn-group button:not(:last-child) {
	border-right: none;
}

.btn-group button:hover {
	background-color: #3e8e41;
}

.btn-group > button:not(:first-child):not(:last-child) {
	border-radius: 0;
}

	.btn-group > button:last-child:not(:first-child){
	border-top-left-radius: 0;
	border-bottom-left-radius: 0;
}

.btn-group > button:first-child:not(:last-child) {
	border-top-right-radius: 0;
	border-bottom-right-radius: 0;
}
"""

# ╔═╡ 1b685538-828a-4c58-bfbf-845e4bd7c26e
function graphdata(g::AbstractGraph)
	vs = ""
	for v in vertices(g)
		vs *= "\t\t{ data: {id: '$v', label: '$v'} },\n"
	end

	es = ""
	for e in edges(g)
		es *= "\t\t{ data: {id: '$(e.src)-->$(e.dst)', source: '$(e.src)', target: '$(e.dst)' } },\n"
	end

	output = """
	{
		nodes: [
	$vs
		],
		edges: [
	$es
		]
	}
	"""
end

# ╔═╡ 2eb9f81b-792c-4e45-a4e5-b3364aa3b391
function minify(input::String)
	input = replace(input, "\n"=>"")
	input = replace(input, "\t"=>" ")
	input = replace(input, r"\s{2,}"=>"")
	input = replace(input, '"'=>"")
	return input
end

# ╔═╡ 6c871573-e5b0-49c4-8b38-4800922483dc
md"""
## Interactive GraphViz
"""

# ╔═╡ aa0a5db5-1c71-4756-8c7f-821b85a06ba7
small_graphs = Dict(
	"bull graph" => :bull,
	"Chvátal graph" => :chvatal,
	"Platonic cubical graph" => :cubical,
	"Desarguesgraph" => :desargues,
	"diamond graph" => :diamond,
	"Platonic dodecahedral graph" => :dodecahedral,
	"Frucht graph" => :frucht,
	"Heawood graph" => :heawood,
	"graph mimicing the classic outline of a house" => :house,
	"house graph, with two edges crossing the bottom square" => :housex,
	"Platonic icosahedral graph" => :icosahedral,
	"Krackhardt-Kite social network graph" => :krackhardtkite,
	"Möbius-Kantor graph" => :moebiuskantor,
	"Platonic octahedral graph" => :octahedral,
	"Pappus graph" => :pappus,
	"Petersen graph" => :petersen,
	"Maze graph used in Sedgewick's Algorithms" => :sedgewickmaze,
	"Platonic tetrahedral graph" => :tetrahedral,
	"skeleton of the truncated cube graph" => :truncatedcube,
	"skeleton of the truncated tetrahedron graph" => :truncatedtetrahedron,
	"skeleton of the truncated tetrahedron digraph" => :truncatedtetrahedron_dir,
	"Tutte graph" => :tutte
)

# ╔═╡ d8a11fab-d87e-434f-9833-2192ca0b302f
md"""
## Graph Styles
"""

# ╔═╡ 7741e7a6-931b-4f47-91bb-2ac61e0e1342
style1 = """
node {
	content: data(id);
	text-valign: center;
	text-halign: center;
	opacity: 0.75;
}

edge {
	curve-style: bezier;
	target-arrow-shape: triangle;
	opacity: 0.5;
}
"""

# ╔═╡ 2a119909-c45c-4154-b85d-ef55ba8fbe8a
style2 = """
node {
  background-color: yellow;
  border-color: black;
  width: 20;
  height: 20;
  min-zoomed-font-size: 12;
  color: #fff;
  font-size: 16;
  z-index: 1;
}

node:selected,
node.start,
node.end {
  height: 30;
  width: 30;
  min-zoomed-font-size: 0;
  font-size: 48;
  border-color: #000;
  border-width: 4px;
  text-outline-color: #000;
  text-outline-width: 10px;
  z-index: 9999;
}

edge {
  min-zoomed-font-size: 12;
  font-size: 8;
  color: #fff;
  line-color: green;
  width: 20;
  curve-style: haystack;
  haystack-radius: 0;
  opacity: 0.5;
}

"""

# ╔═╡ 94ada1f4-9484-4f5d-90b2-a99165dd9f8c
style_green = """
node {
	content: data(id);
	text-valign: center;
	text-halign: center;
	background-color: #126814;
}

edge {
	line-color: #126814;
	opacity: 0.5;
}
"""

# ╔═╡ 401148b3-b3ef-46f9-9869-f3da8d54f6ea
style_dict = Dict(
	"Simple Gray (Directed)" => style1,
	"Green and Yellow" => style2,
	"Simple Green (Undirected)" => style_green,
)

# ╔═╡ cc203fb4-4755-45cc-8b9b-d06048623ec9
md"""
### Interactive Demo

**Graph Appearance:** $(@bind style Select(collect(keys(style_dict))))

**Layout Method:** $(@bind layout Select(["cose", "random", "circle", "concentric", "breadthfirst", "grid"])) **Animate:** $(@bind animate CheckBox(false))

**Sample Graph:** $(@bind sg_name Select(collect(keys(small_graphs))))

"""

# ╔═╡ 824c7e26-5bf8-4272-a4f0-e366271c8221
sg = smallgraph(small_graphs[sg_name])

# ╔═╡ 78405683-feee-4ff1-ac93-675c92d01a12
function graphviz(g::AbstractGraph, name="cy", style=style1; layout="cose", animate=false)
	gd = graphdata(sg)

	layout_block = """
	{
		name: '$layout',
		padding: 5,
		animate: $animate
	}
	"""

	# Testing loading styles from file
	# style_text = read(joinpath(pwd(), "tokyo-railways.cycss"), String)
	# style_text = minify(style_text)
	graph_style = minify(style)

	cyjs = JavaScript("""
		var prnt = currentScript.parentElement;
		var fit_button = prnt.querySelector("button#fit");
		var center_button = prnt.querySelector("button#center");
		var layout_button = prnt.querySelector("button#layout");
		var download_button = prnt.querySelector("button#download");

		var graph_style = '$graph_style';

		var cy = window.cy = cytoscape({
			container: document.getElementById('$name'),
			boxSelectionEnabled: false,
			style: graph_style,
			elements: $gd,
			layout: $layout_block
		});
		fit_button.addEventListener("click", (e) => {cy.fit()});
		center_button.addEventListener("click", (e) => { cy.center() });

		var layout = cy.layout($layout_block);
		layout_button.addEventListener("click", (e) => { layout.run() });

		download_button.addEventListener("click", (e) => {
			var png64 = cy.png()
			var aDownloadLink = document.createElement('a');
			aDownloadLink.href = png64;
			aDownloadLink.target = '_blank';
			aDownloadLink.click();
		});
	""")

	output = htl"""
	<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.4/css/all.min.css" integrity="sha512-1ycn6IcaQQ40/MKBW2W4Rhis/DbILU74C1vSrLJxCq57o941Ym01SwNsOMqvEBFlcgUa6xLiPY/NS5R+E6ztJQ==" crossorigin="anonymous" referrerpolicy="no-referrer" />

			<button id="fit"></button>
			<button id="center"></button>
			<button id="layout"></button>
			<button id="download"></button>

	"""
end

# ╔═╡ 546baf39-09e0-46a1-8df3-b39ba04aa11f
graphviz(sg, "demo", style_dict[style], layout=layout, animate=animate)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
HypertextLiteral = "~0.9.0"
LightGraphs = "~1.3.5"
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.1.0"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "79b9563ef3f2cc5fc6d3046a5ee1a57c9de52495"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.33.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "477bf42b4d1496b454c10cce46645bb5b8a0cf2c"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.12"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═7fb69cf0-fc66-11eb-256f-cb6e9d2db156
# ╠═20ad770b-34a1-49e7-b948-e12f46bde5ee
# ╟─59cef39a-3577-4726-aeff-960ad5b580c5
# ╠═c8ba27b1-4d2b-4afc-89a9-806be91534dc
# ╟─ef077759-88ad-4655-b04f-83c80270d320
# ╟─a626d3c1-ba41-4b91-a2a2-6ff69cd13bce
# ╟─1b685538-828a-4c58-bfbf-845e4bd7c26e
# ╟─2eb9f81b-792c-4e45-a4e5-b3364aa3b391
# ╟─78405683-feee-4ff1-ac93-675c92d01a12
# ╟─6c871573-e5b0-49c4-8b38-4800922483dc
# ╟─aa0a5db5-1c71-4756-8c7f-821b85a06ba7
# ╟─824c7e26-5bf8-4272-a4f0-e366271c8221
# ╟─cc203fb4-4755-45cc-8b9b-d06048623ec9
# ╠═546baf39-09e0-46a1-8df3-b39ba04aa11f
# ╟─d8a11fab-d87e-434f-9833-2192ca0b302f
# ╟─401148b3-b3ef-46f9-9869-f3da8d54f6ea
# ╟─7741e7a6-931b-4f47-91bb-2ac61e0e1342
# ╟─2a119909-c45c-4154-b85d-ef55ba8fbe8a
# ╟─94ada1f4-9484-4f5d-90b2-a99165dd9f8c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_name: bradeneliason/PlutoGraphViz

max_issues_repo_name: bradeneliason/PlutoGraphViz

max_forks_repo_name: bradeneliason/PlutoGraphViz


text:
# tested for Julia v1.2
##### this version works with arbitrary precision arithmetic (BigFloat)
##### !!!! and for data with a (slightly) non-uniform time step (for the discrete argument)

# © Siberian Federal University, 2017
# Authors: Sergey P. Tsarev, Alexander S. Pustoshilov
# Paper reference: S.P.Tsarev, A.A.Kytmanov, Discrete orthogonal polynomials as a tool for detection of small anomalies of time series: a case study of GPS final orbits: https://arxiv.org/abs/2004.00414

# This script robustly approximates a (probably vector-valued) time series with polynomials of large degrees
# (up to a few hundred) in the sense of the best RMS (best least squares) approximation
# any one-dimensional set of scalar or vector-valued data (which formally can be considered a function of a discrete argument)
# with a very large number of points (thousands)
# and find the approximation residue.
# You can approximate the entire dataset by polynomials on the entire data (or a part thereof) available in the file,
# or alternatively using a sliding window (possibly with overlapping).

# Requires auxiliary script service_functions-en.jl with some service functions

############################################################
# Script launch parameters:
############################################################

# julia <script-name.jl> <file> <read_columns> <output_degrees> <windowsparams> <prec>
# where:

# <script-name.jl> - the name of this file

# <file> - name of the input file with the data being processed
# The file contains the data for approximation in columns - a plain text file with numeric data in any standard format: integer, floating point, exponential format.
# The file may contain empty and commented lines (the # character is used for commenting), such lines will be skipped in processing.

# <read_columns> - numbers of columns read from the file;
# numbers are indicated in square brackets,
# the first is the number of the column containing the discrete argument ("time"),
# if the line number should be used as a discrete argument,
# you have to specify 0 as the argument column (see examples);
# second and subsequent column numbers are columns that contain (vector)-function of the disctere argument.
# Examples of specifying column numbers:
# [1,2] - the first column of the file will be read as a discrete argument, the second - as a function of this argument;
# [0,1] - the line number will act as a discrete argument, the first column of the file as its function;
# [2,5:9] - the second column will act as a discrete argument, the  vector-function columns will be 5,6,7,8,9

# <output_degrees> - indicates degrees of approximating polynomials that will be used to output the results.
# Best least squares approximation is performed from degree 0 to the maximum degree specified in this parameter, but only the specified degrees will produce the output files.
# Examples of indicating the degrees:
# [1,2,3,4,5] - approximation by polynomials up to the 5th degree is performed, and the result is output after approximation by polynomials of degrees (1,2,3,4,5);
# [5] - approximation is carried out again by polynomials up to the 5th degree, but the result is output only for approximation by the best approximation polynomial of the 5th degree;
# [1:4, 5:5:20] - approximation by polynomials up to the 20th degree is performed, and the result is output for approximations by polynomials of degrees (1,2,3,4,5,10,15,20).
# The notations 1:4 and 5:5:20 have the same meaning as for loops in Julia.

# <windowsparams> - data processing windows parameters
# The script has the ability to process data in three ways:
# M:N - one window is given;
# M:W:N - windows with one point overlapping;
# M:W:S:N - with windows with a given overlap.
# Here:
# M - starting point of the first window,
# N - last point of the last window (as the last point
# you can specify the end of the file, the "end" keyword is used for this),
# W - window size
# S - windows overlap.
# !!! IMPORTANT NOTE - 1 !!!
# When forming windows, the last window is formed so that
# the last point of this window coincided with the last given point N,
# at the same time, the overlap with the penultimate window may be greater than specified in the parameter S.
# !!! IMPORTANT NOTE - 2 !!!
# If the specified window size is larger than the amount of data available, then the script exits logging the error in the log file.
# *************************************************** *********
# Examples of specifying window parameters:
# 1:end - one window is formed, the entire input file will be processed;
# 1:1000 - one window with a size of 1000 points is formed, the starting point;
# matches the first data point, the window end matches 1000th data point;
# 1:100:1000 - windows with the size of 100 points are formed at the following points (1, 100), (100, 199), (199, 298), ..., (901, 1000);
# 1:100:10:1000 - windows with a size of 100 points are formed with 10 points overlapping: (1, 100), (91, 190), (181, 280), ..., (901, 1000).

# <prec> - precision for BigFloat calculations (in bits)

####################################################### ###########
# SCRIPT RESULTS
####################################################### ###########

# As a result of the script, APPROX directory is formed in the starting directory,
# in which one subdirectory is created with the name <out_dir_name>,
# then in this subdirectory individually for each approximated column from <read_columns>,
# and each degree from <output_degrees> and each processing window,
# a file containing 4 columns is created. Its columns are:
# <discrete arguments>,
# <initial approximated values>,
# <values of the approximating polynomial>,
# <approximation residue = difference between the second and third columns>.
# The name of each file containd the number of the approximated column, the approximation window (line number of the file,
# not including comments) and the degree of approximating polynomial.
# All results in the files are printed with @sprintf ("% 1.15e", ...) - you can change the number of printed digits if necessary.
# For convenience of plotting the results using gnuplot we set the file extension: *.gplot_4
# **********************************************************
# Additionally, a file is created in the APPROX directory (appending new lines if the file exists) result_max_resid.txt,
# which contains the maximum(abs(residues)) at the maximal degree approximation for each window for each of the input columns,
# and the following data is written as one new line:
# - <out_dir_name> name of the output approximation directory;
# - <current_window> window being processed;
# - maximum(abs(residues)) for each approximated column in the current window with the maximal degree of approximation.

# To keep the errors, the errorlog_approx.log file is created in the script launch directory,
# which displays the following types of errors with the name of the experiment and the start time:
# "Incorrect data, 0 lines read" - occurs if the script could not
# read the source file or invalid column numbers were specified;
# "Invalid data, window length greater than data size" -
# informs that when setting windows parameters the window size was incorrectly selected
# and there is not enough data to cover at least 1 window.

####################################################### ###########
# Launch example:
############################################################
# julia <script.jl> F-testdata.txt [1:4] [0,1,3,5:5:20] 1:end 256
# 5:5:20 in <output_degrees> means a list of degrees with start 5, step 5 and end 20
# 1:end means processing the entire "time" interval, from the beginning of the file to the end as a single window.
# 256 means that all computations will be performed with precision 256 bits
############################################################

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

if length(ARGS) != 5
    println(length(ARGS))
    println("warning: wrong number of script parameters")
    exit()
end

out_dir_name = "nonequ_$(ARGS[1])_col$(ARGS[2])_degs$(ARGS[3])_w[$(ARGS[4])]_prec$(ARGS[5])"
out_dir_name = replace(out_dir_name, ":"=>"-")

errorlog = "errorlog_approx.log"
max_resid_file = "APPROX/result_max_resid.txt"

using DelimitedFiles
using LinearAlgebra
using Printf
using Dates

# Loading the module with service functions

include("service_functions-en.jl")

# Beginning of the main approximation script

filename = ARGS[1]
read_columns = make_numberlist(ARGS[2])
prec = parse(Int, ARGS[5])
setprecision(prec)

DATA = readbigfile(filename, read_columns, BigFloat)
length_DATA = length(DATA[:,1])

if length_DATA == 0
    add_log_messange(errorlog, "$(out_dir_name) Incorrect data, 0 lines read")
    exit()
end

output_degrees = make_numberlist(ARGS[3])
windows_list = make_windows(ARGS[4], length_DATA)

if length(windows_list) == 0
    add_log_messange(errorlog, "$(out_dir_name) Incorrect data, window size larger then the data size")
    exit()
end

println("output dir: ", out_dir_name)
flush(stdout)

mkpath("APPROX/$out_dir_name")

out_approx_max_resid = open(max_resid_file, "a")

cd("APPROX/$out_dir_name")

# maxdeg = output_degrees[end]
maxdeg = maximum(output_degrees)
sym_num_max_deg = length("$maxdeg")

# orthogonality tolerance:
orttol = big(2.0^(-prec * 90 / 100))
println("orttol=", orttol)
## 1 and 2 into BigFloat:
ONE = big(1.0)
TWO = big(2.0)

# Processing the sliding windows

for current_window in windows_list
    println("Processing window: $current_window")
    flush(stdout)
    # Preparing the table of discrete orthogonal polynomial values for the CURRENT window
    # pnts array contains the values of the discrete argument (the first specified column) of the one-dimensional discrete (vector) function (the remaining specified columns)
    pnts = DATA[current_window[1]:current_window[2],1]
    NN = length(pnts)
    p1 = pnts[1]
    pNN = pnts[end]
    # normalize the discrete argument to the interval [-1,1]:
    pnts11 = zeros(BigFloat, NN)
    for k=1:NN
        pnts11[k] = TWO*(pnts[k] - p1) / (pNN - p1) - ONE
    end
    pols = zeros(BigFloat, NN, maxdeg + 1)
    # pols[j, k] contains the value of orthogonal polynomial of degree (k-1) at the jth point of the discrete argument (specified by the first column)
    # At first we fill the array pols with the values of Legendre polynomials (they are NOT orthogonal on the discrete lattice!), then we use a special variant of Gram-Schmidt orthogonalization (and the we normalize the polynomials)
    # first polynomial is constant:
    pols[:,1] = ones(BigFloat, NN)
    # the second is a linear function of the discrete argument (specified by the first column)
    pols[:,2] = pnts11
    for d = 3:(maxdeg + 1)
        n = big(d - 2.0)
        # pols[:,d] =   (((2 * n + 1) * (pnts11 .* pols[:,d - 1])) - (n * pols[:,d - 2])) / (n + 1);
        pols[:,d] =   (((TWO * n + ONE) * (pnts11 .* pols[:,d - 1])) - (n * pols[:,d - 2])) / (n + ONE);
        println("degree = ", d - 1)
        flush(stdout)
    end
    println("Legendre polynomials done")
    norm0 = sqrt(sum(abs2, pols[:,1]))
    pols[:,1] = pols[:,1] / norm0
    println("Orthogonalization started")
    flush(stdout)
    # The obtained Legendre polynomials now need to be orthogonalized and normalized.
    # pols[:,2] SHOULD be orthogonalized for non-equidistant points  !!!!
    for n = 2:(maxdeg + 1)
        u = pols[:,n]
        repN = 0
        nev = big(1000000.0)
        while ( nev > 2 * NN * orttol )
            repN = repN + 1
            for i = 1:(n - 1)
                cf = sum(u .* pols[:,i])
                u = u - cf * pols[:,i]  # Wikipedia version
                # - this version of Gram-Schmidt algorithm turned out to be robust even for Float64 computations!!!!
            end
            nev = big(0.0)
            norm0 = sqrt(sum(abs2, u))
            for i = 1:(n - 1)
                nev = max(nev, abs(sum(u .* pols[:,i])))
                # nev = max(nev, abs(sum(u .* pols[:,i]))/norm0)
            end
        end
        norm0 = sqrt(sum(abs2, u))
        pols[:,n] = u / norm0
        println("degree = ", n - 1, "   max nev after orthogonalization = ", @sprintf("%1.15e",nev),
        "   repN = ", repN, "   norm was: ", @sprintf("%1.15e",norm0))
        flush(stdout)
    end
    # End of orthogonalization - in the array pols[:,:] we have the set of values of the normalized discrete orthogonal Hahn polynomials for the CURRENT window
    println("Orthogonalization done")
    flush(stdout)
    max_resid_in_col = zeros(BigFloat, length(read_columns)-1)
    for ind_cn = 2:length(read_columns)
        data_col = DATA[current_window[1]:current_window[2],ind_cn]
        data_arg = DATA[current_window[1]:current_window[2], 1]
        data_cn = copy(data_col)
        pol_approx = zeros(BigFloat, NN)
        # approximation of the input data (individually for each column) by polynomials - since they are orthonormal, a simple scalar product is enough!
        for dn in 0:maxdeg
            p = pols[:,dn + 1]
            cf = sum(data_cn .* p)
            data_cn = data_cn .- (cf * p)
            pol_approx = pol_approx .+ (cf * p)
            print("col: ", read_columns[ind_cn], "   win: ", current_window, "   deg: ", dn)
            flush(stdout)
            if (dn in output_degrees)
                output_filename = "col$(read_columns[ind_cn])_win$(current_window)_deg$(print_with_zero(dn, sym_num_max_deg)).gplot_4"
                output_filename = replace(output_filename, " "=>"_")
                output_filename = replace(output_filename, ","=>"")
                out_file = open(output_filename, "w")
                for k = 1:length(pnts)
                    @printf(out_file, "%1.15e\t%1.15e\t%1.15e\t%1.15e\n", data_arg[k], data_col[k], pol_approx[k], data_cn[k])
                end
                println(" written!")
                close(out_file)
            else
                println()
                flush(stdout)
            end
        end # dn in 0:maxdeg
        max_resid_in_col[ind_cn-1] = maximum(abs.(data_cn))
    end # ind_cn
    println("$out_dir_name $(current_window) $(join(map(x->@sprintf("%1.15e",x),max_resid_in_col), "\t"))")
    println(out_approx_max_resid, "$out_dir_name $(current_window) $(join(map(x->@sprintf("%1.15e",x),max_resid_in_col), "\t"))")
    flush(out_approx_max_resid)
end # current_window

close(out_approx_max_resid)
exit()

max_stars_repo_name: sptsarev/high-deg-polynomial-fitting

max_issues_repo_name: sptsarev/high-deg-polynomial-fitting

max_forks_repo_name: sptsarev/high-deg-polynomial-fitting


text:
immutable Conv2dParameters
    x_max :: Float64
    y_max :: Float64
    filter :: Function
    filter_dx :: Function
    filter_dy :: Function
    pix_size_x :: Float64
    pix_size_y :: Float64
    guess_error_x :: Float64
    guess_error_y :: Float64
end
immutable Conv2d <: SuperRes
    x_max :: Float64
    y_max :: Float64
    filter :: Function
    filter_dx :: Function
    filter_dy :: Function
    eval_grid_x :: Array{Float64}
    eval_grid_y :: Array{Float64}
    guess_grid_x :: Array{Float64}
    guess_grid_y :: Array{Float64}
    n_x :: Int64
    n_y :: Int64
    function Conv2d(params :: Conv2dParameters)
        eval_grid_x_tmp = 0:params.pix_size_x:params.x_max
        eval_grid_y_tmp = 0:params.pix_size_y:params.y_max
        eval_grid_x = [x for x in eval_grid_x_tmp, y in eval_grid_y_tmp][:]
        eval_grid_y = [y for x in eval_grid_x_tmp, y in eval_grid_y_tmp][:]
        guess_grid_x_tmp = 0:params.guess_error_x:params.x_max
        guess_grid_y_tmp = 0:params.guess_error_y:params.y_max
        guess_grid_x = [x for x in guess_grid_x_tmp, y in guess_grid_y_tmp][:]
        guess_grid_y = [y for x in guess_grid_x_tmp, y in guess_grid_y_tmp][:]
        return new(params.x_max, params.y_max,
                   params.filter, params.filter_dx, params.filter_dy,
                   eval_grid_x, eval_grid_y,
                   guess_grid_x, guess_grid_y,
                   length(eval_grid_x_tmp), length(eval_grid_y_tmp))
    end
end
function psi(model :: Conv2d, theta :: Vector{Float64})
    # This function computes the direct problem for a single point theta
    return Float64[model.filter(model.eval_grid_x[i] - theta[1], model.eval_grid_y[i] - theta[2]) for i in 1:length(model.eval_grid_x)]
end

function dpsi(model :: Conv2d, theta :: Vector{Float64})
    # This function computes the gradient of psi and adds it to eval_array_grad field
    dx = Float64[- model.filter_dx(model.eval_grid_x[i] - theta[1], model.eval_grid_y[i] - theta[2]) for i in 1:length(model.eval_grid_x)]
    dy = Float64[- model.filter_dy(model.eval_grid_x[i] - theta[1], model.eval_grid_y[i] - theta[2]) for i in 1:length(model.eval_grid_x)]
    return [dx dy]
end

function phi(model :: Conv2d, parameters :: Matrix{Float64}, weights :: Vector{Float64})
    # Evaluates the direct problem
    return sum([weights[i] * psi(model, vec(parameters[:, i])) for i in 1:length(weights)])
end

function getStartingPoint(model :: Conv2d, v :: Vector{Float64})
    # This function gets an initial value for the location
    guess_array = map((x,y) -> dot(v, psi(model, [x;y])), model.guess_grid_x, model.guess_grid_y)
    ind = indmin(guess_array);
    res = [model.guess_grid_x[ind]; model.guess_grid_y[ind]]
    return res
end

parameterBounds(model :: Conv2d) = [0.0; 0.0], [model.x_max; model.y_max] # Bounds for the parameters
dim(model :: Conv2d) = 2

type DynamicConv2dParameters
    K :: Int64
    tau :: Float64
    v_max :: Float64
    num_v :: Int64
end

type DynamicConv2d<: SuperRes
    static :: Conv2d
    v_max :: Float64
    eval_grid_x :: Array{Float64}
    eval_grid_y :: Array{Float64}
    eval_grid_t :: Array{Float64}
    guess_grid_x :: Array{Float64}
    guess_grid_y :: Array{Float64}
    guess_grid_vx :: Array{Float64}
    guess_grid_vy :: Array{Float64}
    eval_array :: Vector{Float64}
    eval_array_grad :: Array{Float64}
    function DynamicConv2d(static :: Conv2d, params::DynamicConv2dParameters)
        times = linspace(-params.K*params.tau, params.K*params.tau, 2*params.K + 1)
        eval_grid_x = [x for x in static.eval_grid_x, t in times][:]
        eval_grid_y = [y for y in static.eval_grid_y, t in times][:]
        eval_grid_t = [t for x in static.eval_grid_y, t in times][:]
        velocities = linspace(-params.v_max, params.v_max, params.num_v)
        guess_grid_x = [x for x in static.guess_grid_x, vx in velocities, vy in velocities][:]
        guess_grid_y = [y for y in static.guess_grid_y, vx in velocities, vy in velocities][:]
        guess_grid_vx = [vx for y in static.guess_grid_y, vx in velocities, vy in velocities][:]
        guess_grid_vy = [vy for y in static.guess_grid_y, vx in velocities, vy in velocities][:]
        return new(static, params.v_max,
                   eval_grid_x, eval_grid_y, eval_grid_t,
                   guess_grid_x, guess_grid_y, guess_grid_vx, guess_grid_vy)
    end
end
function psi(model :: DynamicConv2d, theta :: Vector{Float64}, w = 1.0, acc=false)
    # This function computes the direct problem for a single point theta
    return Float64[model.static.filter(model.eval_grid_x[i] - theta[1] - theta[3]*model.eval_grid_t[i],
                                       model.eval_grid_y[i] - theta[2] - theta[4]*model.eval_grid_t[i]) for i in 1:length(model.eval_grid_x)]

end

function dpsi(model :: DynamicConv2d, theta :: Vector{Float64})
    # This function computes the gradient of psi and adds it to eval_array_grad field
    dx =  Float64[- model.static.filter_dx(model.eval_grid_x[i] - theta[1] - theta[3]*model.eval_grid_t[i],
                                           model.eval_grid_y[i] - theta[2] - theta[4]*model.eval_grid_t[i]) for i in 1:length(model.eval_grid_x)]
    dy =  Float64[- model.static.filter_dy(model.eval_grid_x[i] - theta[1] - theta[3]*model.eval_grid_t[i],
                                           model.eval_grid_y[i] - theta[2] - theta[4]*model.eval_grid_t[i]) for i in 1:length(model.eval_grid_x)]
    dvx =  Float64[-model.eval_grid_t[i] * model.static.filter_dx(model.eval_grid_x[i] - theta[1] - theta[3]*model.eval_grid_t[i],
                                                                  model.eval_grid_y[i] - theta[2] - theta[4]*model.eval_grid_t[i]) for i in 1:length(model.eval_grid_x)]
    dvy =  Float64[-model.eval_grid_t[i] * model.static.filter_dy(model.eval_grid_x[i] - theta[1] - theta[3]*model.eval_grid_t[i],
                                                                  model.eval_grid_y[i] - theta[2] - theta[4]*model.eval_grid_t[i]) for i in 1:length(model.eval_grid_x)]
    return [dx dy dvx dvy]
end

function phi(model :: DynamicConv2d, parameters :: Matrix{Float64}, weights :: Vector{Float64})
    # Evaluates the direct problem
    return sum([weights[i] * psi(model, vec(parameters[:, i])) for i in 1:length(weights)])
end

function getStartingPoint(model :: DynamicConv2d, v :: Vector{Float64})
    # This function gets an initial value for the location
    guess_array = map((x,y, vx, vy) -> dot(v, psi(model, [x;y;vx;vy])), model.guess_grid_x, model.guess_grid_y, model.guess_grid_vx, model.guess_grid_vy)
    ind = indmin(guess_array);
    return vec([model.guess_grid_x[ind]; model.guess_grid_y[ind]; model.guess_grid_vx[ind]; model.guess_grid_vy[ind]])
end

parameterBounds(model :: DynamicConv2d) = [0.0; 0.0; -model.v_max; -model.v_max], [model.static.x_max; model.static.y_max; model.v_max; model.v_max] # Bounds for the parameters


text:
module Utils

using Lynx
using Luxor
using MacroTools

export triangle
export Mouse
export Map
export @init
export @setup
export @var

mutable struct Mouse
    x::Real
    y::Real
    dx::Real
    dy::Real
    zoom::Real
    Mouse() = new()
end

function Mouse(canvas::Canvas)
    this = Mouse()
    this.zoom = 1
    onmousemove(canvas) do event
        if isdefined(this, :x)
            this.dx = event.x - this.x
            this.dy = event.y - this.y
        end
        this.x = event.x
        this.y = event.y
    end
    onscroll(canvas) do event
        dir = event.direction
        this.zoom += 0.1 * dir
    end
    return this
end

macro var(expr)
    if @capture(expr, name_(args__) = block_)
        var = Symbol("@$name")
        return esc( :( $(var)(source::LineNumberNode, ::Module, $(args...)) = $(block) ) )
    end
end

macro init(expr)
    if @capture(expr, function Name_(args__) body__ end)
        return quote
            function $(Name)($(args...))
                this = new()
                $(body...)
                return this
            end
        end |> esc
    elseif Meta.isexpr(expr, :block)
        return quote
            this = new()
            $(expr.args...)
            return this
        end |> esc
    end
end

macro setup(block)
    globals = Symbol[:width, :height]
    res = Expr[
        :(width = Lynx.@width),
        :(height = Lynx.@height)
    ]
    if @capture(block, function() body__ end | begin body__ end)
        foreach(body) do expr
            if @capture(expr, var_ = x_)
                push!(globals, var)
            end
            push!(res, expr)
        end
    end
    return esc(:(
        $(map(x -> :( $x = nothing ), globals)...);
        function setup()
            global $(globals...);
            $(res...)
        end
    ))
end

triangle(x::Real, y::Real, d::Real) = Point[
    Point(x, y - d),
    Point(x - d, y + d),
    Point(x + d, y + d),
]

triangle(center::Point, d::Real) = triangle(center.x, center.y, d)

function triangle(x::Real, y::Real, d::Real, action::Symbol)
    vertices = triangle(x, y, d)
    poly(vertices, action, close = true)
end

end # Utils

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
module Types
    using Configurations
    import Pluto: Pluto, Token, Notebook
    export NotebookSession, RunningNotebookSession, QueuedNotebookSession, FinishedNotebookSession, SliderServerSettings, ExportSettings, PlutoDeploySettings
    ###
    # SESSION DEFINITION

    abstract type NotebookSession end

    Base.@kwdef struct RunningNotebookSession <: NotebookSession
        path::String
        current_hash::String
        notebook::Pluto.Notebook
        original_state
        token::Token=Token()
        bond_connections::Dict{Symbol,Vector{Symbol}}
    end

    Base.@kwdef struct QueuedNotebookSession <: NotebookSession
        path::String
        current_hash::String
    end

    Base.@kwdef struct FinishedNotebookSession <: NotebookSession
        path::String
        current_hash::String
        original_state
    end

    ###
    # CONFIGURATION

    @option struct SliderServerSettings
        exclude::Vector=String[]
        port::Integer=2345
        host="127.0.0.1"
        simulated_lag::Real=0
        serve_static_export_folder::Bool=true
    end

    @option struct ExportSettings
        output_dir::Union{Nothing,String}=nothing
        exclude::Vector=String[]
        ignore_cache::Vector=String[]
        pluto_cdn_root::Union{Nothing,String}=nothing
        baked_state::Bool=true
        baked_notebookfile::Bool=true
        offer_binder::Bool=true
        disable_ui::Bool=true
        cache_dir::Union{Nothing,String}=nothing
        slider_server_url::Union{Nothing,String}=nothing
        binder_url::Union{Nothing,String}=nothing
        create_index::Bool=true
    end

    @option struct PlutoDeploySettings
        SliderServer::SliderServerSettings=SliderServerSettings()
        Export::ExportSettings=ExportSettings()
        Pluto::Pluto.Configuration.Options=Pluto.Configuration.Options()
    end

end


text:
## New Drive for the CLT parallel
tic()
n_grid = [2^i for i = 8:18]
numRuns = 50
N = parse(Int, ARGS[1])

a = @spawn test_CLTExp("", numRuns, n_grid, 8675309000, N, 2.)
b = @spawn test_CLTExp("", numRuns, n_grid, 5164174290, N, 2.)
c = @spawn test_CLTExp("", numRuns, n_grid, 5167462266, N, 2.)

######
file_a = fetch(a)
file_b = fetch(b)
file_c = fetch(c)

println( file_a )
println( file_b )
println( file_c )

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)
data_t[:, 1] += numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] += 2numRuns
data = vcat(data, data_t)

#strip the name of file_a to make the numbers better
indx = rsearch(file_a, "_")[1]
f = open("$(file_a[1:indx])parallel_results.csv", "w")
writedlm(f,  header)
writedlm(f,  data)
close(f)

println("Num Paths: \t $(numRuns) \t Time:", toc() )