text:
"""
    struct IntSetVar <: AbstractIntSetVar

A set variable, that can be a subset of an integer interval.
"""
struct IntSetVar <: AbstractVar
    onDomainChange      ::Array{Constraint}
    domain              ::IntSetDomain
    id                  ::String
end

"""
    function IntSetVar(min::Int, max::Int, id::String, trailer::Trailer)

Create an `IntSetVar` with a domain being the subsets of the integer range [`min`, `max`] with the `id` string identifier
and that will be backtracked by `trailer`.
"""
function IntSetVar(min::Int, max::Int, id::String, trailer::Trailer)
    offset = min - 1

    dom = IntSetDomain(trailer, min, max)

    return IntSetVar(Constraint[], dom, id)
end

"""
    SetModification(; excluded::Array{Int}, required::Array{Int})

Stock the two types of modification of a IntSetVar.
"""
struct SetModification
    excluded::Array{Int}
    required::Array{Int}

    function SetModification(;excluded::Union{Nothing,Array{Int}}=nothing, required::Union{Nothing,Array{Int}}=nothing)
        if isnothing(excluded)
            excluded = Int[]
        end
        if isnothing(required)
            required = Int[]
        end
        return new(excluded, required)
    end
end

"""
    Base.length(set::SetModification)

a generic function length is needed for all modifications.
"""
function Base.length(set::SetModification)
    return length(set.excluded) + length(set.required)
end

"""
    mergeSetModifications!(modif1, modif2)

Add the modifications in `modif2` at the end of `modif1`
"""
function mergeSetModifications!(modif1::SetModification, modif2::SetModification)
    append!(modif1.required, modif2.required)
    append!(modif1.excluded, modif2.excluded)
    return
end

"""
    isbound(x::IntSetVar)

Check whether x has an assigned value (meaning its domain only contains one subset)
"""
isbound(x::IntSetVar) = x.domain.requiring_index.value - 1 == x.domain.excluding_index.value

"""
    assignedValue(x::IntSetVar)

Return the assigned value of `x`, i.e. the only subset it contains. Throw an error if `x` is not bound.
"""
function assignedValue(x::IntSetVar)
    @assert isbound(x)

    return required_values(x.domain)
end

function Base.show(io::IO, var::IntSetVar)
    print(io, var.id, ": ", var.domain)
end

function Base.show(io::IO, ::MIME"text/plain", var::IntSetVar)
    println(io, typeof(var), ": ", var.id)
    print(io, "   ", var.domain)
end


text:
# 2D LS-RTM gradient test with 1 source
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020
#
# Ziyi Yin, ziyi.yin@gatech.edu
# Updated July 2021

parsed_args = parse_commandline()

nlayer = parsed_args["nlayer"]
tti = parsed_args["tti"]
fs =  parsed_args["fs"]

### Model
model, model0, dm = setup_model(tti, 4)
q, srcGeometry, recGeometry, info = setup_geom(model)
dt = srcGeometry.dt[1]

opt = Options(sum_padding=true, free_surface=fs)
F = judiModeling(info, model, srcGeometry, recGeometry; options=opt)
F0 = judiModeling(info, model0, srcGeometry, recGeometry; options=opt)
J = judiJacobian(F0, q)

# Observed data
dobs = F*q
dobs0 = F0*q

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

@testset "LSRTM gradient test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
	# Gradient test
	ftol = (tti && fs) ? 1f-1 : 5f-2
	h = 5f-2
	maxiter = 5
	err1 = zeros(maxiter, 2)
	err2 = zeros(maxiter, 2)
	h_all = zeros(maxiter)

	# LS-RTM gradient and function value for m0
	Jm0, grad = lsrtm_objective(model0, q, dobs, dm; options=opt)
	Jm01, grad1 = lsrtm_objective(model0, q, dobs, dm; options=opt, nlind=true)

	# Perturbation
	dmp = 2f0*circshift(dm, 10)
	dJ = dot(grad, dmp)
	dJ1 = dot(grad1, dmp)

	for j=1:maxiter
		dmloc = dm + h*dmp
		# LS-RTM gradient and function falue for m0 + h*dm
		Jm, _ = lsrtm_objective(model0, q, dobs, dmloc; options=opt)
		Jm1, _ = lsrtm_objective(model0, q, dobs, dmloc; options=opt, nlind=true)
		@printf("h = %2.2e, J0 = %2.2e, Jm = %2.2e \n", h, Jm0, Jm)
		@printf("h = %2.2e, J01 = %2.2e, Jm1 = %2.2e \n", h, Jm01, Jm1)
		# Check convergence
		err1[j, 1] = abs(Jm - Jm0)
		err1[j, 2] = abs(Jm1 - Jm01)
		err2[j, 1] = abs(Jm - Jm0 - h*dJ)
		err2[j, 2] = abs(Jm1 - Jm01 - h*dJ1)

		j == 1 ? prev = 1 : prev = j - 1
		for i=1:2
			@printf("h = %2.2e, e1 = %2.2e, rate = %2.2e", h, err1[j, i], err1[prev, i]/err1[j, i])
			@printf(", e2  = %2.2e, rate = %2.2e \n", err2[j, i], err2[prev, i]/err2[j, i])
		end
		h_all[j] = h
		h = h * .5f0
	end

	for i=1:2
		# CHeck convergence rates
		rate_1 = sum(err1[1:end-1, i]./err1[2:end, i])/(maxiter - 1)
		rate_2 = sum(err2[1:end-1, i]./err2[2:end, i])/(maxiter - 1)

		# This is a linearized problem, so the whole expansion is O(dm) and
		# "second order error" should be first order
		@test isapprox(rate_1, 2f0; rtol=ftol)
		@test isapprox(rate_2, 4f0; rtol=ftol)
	end

	# test that with zero dm we get the same as fwi_objective for residual
	ENV["OMP_NUM_THREADS"]=1
	Jls, gradls = lsrtm_objective(model0, q, dobs, 0f0.*dm; options=opt, nlind=true)
	Jfwi, gradfwi = fwi_objective(model0, q, dobs; options=opt)
	@test isapprox(Jls, Jfwi; rtol=0f0, atol=0f0)
	@test isapprox(gradls, gradfwi; rtol=0f0, atol=0f0)
end

# Test if lsrtm_objective produces the same value/gradient as is done by the correct algebra
cases = [(true, false, true), (false, false, true), (true, true, false), (true, false, false), (false, true, false), (false, false, false)]	# DFT and optimal_checkpointing normally don't co-exist
@testset "LSRTM gradient linear algebra test with $(nlayer) layers and tti $(tti) and freesurface $(fs) and isic $(isic) and optimal_checkpointing $(optchk) and DFT $(dft)" for (isic,optchk,dft) = cases

	ftol = 5f-4
	freq = dft ? [[2.5, 4.5],[3.5, 5.5],[10.0, 15.0], [30.0, 32.0]] : []
	J.options.free_surface = fs
	J.options.isic = isic
	J.options.optimal_checkpointing = optchk
	J.options.frequencies = freq

	dm1 = 2f0*circshift(dm, 10)
	d_res = dobs0 + J*dm1 - dobs
	Jm0_1 = 0.5f0 * norm(d_res)^2f0
	grad_1 = J'*d_res

	opt = J.options
	Jm0, grad = lsrtm_objective(model0, q, dobs, dm1; options=opt, nlind=true)
	Jm01, grad1 = lsrtm_objective(model0, q, dobs-dobs0, dm1; options=opt, nlind=false)

	@test isapprox(grad, grad_1; rtol=ftol)
	@test isapprox(Jm0, Jm0_1; rtol=ftol)
	@test isapprox(grad, grad1; rtol=ftol)
	@test isapprox(Jm0, Jm01; rtol=ftol)

end

ext: jl


text:
"""

```
mddf(trajectory::Trajectory, options::Options)
```

Function that computes the minimum-distance distribution function, atomic contributions, and
KB integrals, given the `Trajectory` structure of the simulation and, optionally, parameters
given as a second argument of the `Options` type. This is the main function of the `ComplexMixtures`
package.

### Examples

```julia-repl
julia> trajectory = ComplexMixtures.Trajectory("./trajectory.dcd",solute,solvent);

julia> results = ComplexMixtures.mddf(trajectory);

```

or

```julia-repl
julia> options = Options(lastframe=1000);

julia> results = ComplexMixtures.mddf(trajectory,options);

```

"""
mddf(trajectory::Trajectory) = mddf(trajectory,Options())

# Choose among serial or parallel version, and self and non-self versions

# In both the self and non-self cases, the number of random samples is
# n_random_samples*solvent.nmols. However, in the non-self distribution
# the sampling of solvent distances is proportional to the number of solute
# molecules, thus the md count has to be averaged by the solute.nmols. In the
# case of the self-distribution, we compute n(n-1)/2 distances, and we will
# divide this by n random samples, which is the sampling of the random
# distribution. Therefore, we must weight the self-distance count by dividing
# it by (n-1)/2, so that we have a count proportional to n as well, leading
# to the correct weight relative to the random sample.

function mddf(trajectory::Trajectory, options::Options)

  # Set random number generator
  RNG = init_random(options)

  if options.nthreads < 0
    nthreads = Threads.nthreads()
  else
    nthreads = options.nthreads
  end
  # If the solute and the solvent are the same
  if trajectory.solute.index == trajectory.solvent.index
    samples = Samples(md=(trajectory.solvent.nmols-1)/2,
                      random=options.n_random_samples)
    if nthreads == 1
      mddf_linkedcells(trajectory,options,samples,RNG,mddf_frame_self!)
    else
      mddf_linkedcells_parallel(trajectory,options,samples,RNG,mddf_frame_self!)
    end
  # If solute and solvent are different subsets of the simulation
  else
    samples = Samples(md=trajectory.solute.nmols,
                      random=options.n_random_samples)
    if nthreads == 1
      mddf_linkedcells(trajectory,options,samples,RNG,mddf_frame!)
    else
      mddf_linkedcells_parallel(trajectory,options,samples,RNG,mddf_frame!)
    end
  end
end


text:
# data for Constrained Delaunay Triangulation (CDT)
points = convert(Array{Float64,2}, V')
points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
edges_list = convert(Array{Int64,2}, hcat(EV...)')
edge_boundary = [true for k=1:size(edges_list,1)]
triangles = Triangle.constrained_triangulation(points,points_map,edges_list)

    # edges of the triangulation
    ev = map(sort,cat([[[u,v], [v,w], [w,u]] for (u,v,w) in triangles]))
    # remove duplicated edges from triangulation
    ev_nodups = collect(Set(ev))
    ##Plasm.view(Plasm.numbering(0.35)((V,[[[k] for k=1:size(V,2)], ev_nodups])))
    # dictionary o original edges
    edge_dict = Dict(zip(EV,1:length(EV)))
    triaedges = [edge_dict[[u,v]] for (u,v) in ev if haskey(edge_dict, [u,v]) ]
    # subdivide original edges between inner and outer
    edge_boundary = Vector{Bool}(undef,10)
    counters = zeros(size(edges_list,1))
    for e in triaedges
        counters[e]+=1
    end
    for e in triaedges
        if counters[e] == 1 # outer edge
            edge_boundary[e] = false
        elseif counters[e] == 2 # inner edge
            edge_boundary[e] = true
        else
            error("in computing inner loops")
        end
    end
    # compute inner triangles
    inneredges = Array{Array{Int64,1},1}()
    for (k,value) in enumerate(edge_boundary)
       if value==true
           push!(inneredges, EV[k], reverse(EV[k]))
       end
    end
    # compute hole(s): wheater all triangle edges are inneredges
    holes = Array{Array{Int64,1},1}()
    for (k,(u,v,w)) in enumerate(triangles)
        triangle = [[u,v],[v,w],[w,u]]
        if setdiff(triangle,inneredges)==[]
            push!(holes, triangles[k])
        end
    end

    # all triangles are generated counterclockwise by Triangle.jl
    # Hence, close points to the left of directed hole edges must
    # be inner to the hole, and hence external to polygon, so that
    # the corresponding triangle is cancelled from the triangle array.
    # This must be repeated for each triangle classified as a hole.

    holes_list = Array{Array{Float64,1},1}()
    epsilon = 1e-5
    for (u,v,w) in holes  # each hole is a triangle
        a = normalize(V[:,v] - V[:,u]) #    unit side vector of first edge
        b = [-a[2,:]; a[1,:]]   #   unit orthogonal vector
        m = (V[:,u] + V[:,v])/2 #   middle side point
        pointinhole = m + epsilon * m  # very close to hole's first edge
        push!(holes_list, pointinhole)
    end
    triangles = [tria for tria in triangles if !(tria in holes)]

points = convert(Array{Float64,2}, V')
points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
edges_list = convert(Array{Int64,2}, hcat(EV...)')
edge_boundary = [true for k=1:size(edges_list,1)]
holes_list = convert(Array{Float64,2}, hcat(holes_list...)')
#triangles = Triangle.constrained_triangulation(points,points_map,edges_list,edge_boundary,holes_list)
triangles = [tria for tria in triangles if !(tria in holes)]

Plasm.view(Plasm.numbering(0.35)((V,[[[k] for k=1:size(V,2)], ev_nodups, triangles])))
Plasm.view(V,triangles)

#
# points = [0. 0.; 4. 0.; 2. 3.; 8. 0.; 6. 3.; 4. 6.]
# points_map = Array{Int64,1}(collect(1:1:size(points)[1]))
# edges_list = Array{Int64,2}([1 2; 2 3; 3 1; 2 4; 4 5; 5 2; 3 5; 5 6; 6 3])
# edge_boundary = [false,true,false,false,false,true,true,false,false]
# holes_list = [4. 2.]
# triangles = Triangle.constrained_triangulation(points,points_map,edges_list,edge_boundary,holes_list)


text:
# ] activate .
# in vscode, display(gcf()) or just gcf() after each plot.
# VSCode use won't be interactive.
using Gravityfield
include(joinpath(@__DIR__, "src/plotfield.jl"))

"Plot a sphere for phun"
let
    n = 100
    u = range(0, 2π, length = n)
    v = range(0, π, length = n)
    x = cos.(u) * sin.(v)'m
    y = sin.(u) * sin.(v)'m
    z = ones(n) * cos.(v)'m
    np = nest(x, y, z)
    wrap_surface_and_countour(np, title = "Just a sphere")
end

# 3d space, plot gravity at a plane offset by a moon radius,
# crop to moon and earth
let
    bodies = bodies_3d()
    x, y = xygrids( xmin = -d_moon_earth/5,
                            xmax = d_moon_earth * 6/5,
                            ymin = -d_moon_earth*(6/5 + 1/5)/2,
                            ymax = d_moon_earth*(6/5 + 1/5)/2,
                            n = 300)
    nxyz = nest(x, y, 1.45456e6m)
    np = nest(x, y, gfieldvals(nxyz, bodies)s²)
    wrap_surface_and_countour(np, title = "1 Gravity acceleration")
end

# 3d space, plot gravity in the planet and moon plane,
# cap the plot roughly at moon's gravity.
let
    bodies = bodies_3d()
    x,y = xygrids( xmin = -d_moon_earth/5,
                            xmax = d_moon_earth * 6/5,
                            ymin = -d_moon_earth*(6/5 + 1/5)/2,
                            ymax = d_moon_earth*(6/5 + 1/5)/2,
                            n = 300)
    nxyz = nest(x, y, 0m)
    np = nest(x, y, clamp.(gfieldvals(nxyz, bodies)s², 0m, 2m))
    wrap_surface_and_countour(np, title = "2 Clamped gravitational acceleration")
end

# 3d space, plot gravity in the planet and moon plane,
# cap the plot at 0.5 m/s²
let
    bodies = bodies_3d()
    x,y = xygrids( xmin = -d_moon_earth/5,
                            xmax = d_moon_earth * 6/5,
                            ymin = -d_moon_earth*(6/5 + 1/5)/2,
                            ymax = d_moon_earth*(6/5 + 1/5)/2,
                            n = 300)
    nxyz = nest(x, y, 0m)
    np = nest(x, y, clamp.(gfieldvals(nxyz, bodies)*s², 0m, 0.5m))
    wrap_surface_and_countour(np, title = "3 Clamped gravitational acceleration")
end

# 3d space, plot gravity in the planet and moon plane,
# cap the plot at 0.5 m/s²
let
    bodies = bodies_3d()
    x,y = xygrids( xmin = -d_moon_earth/5,
                            xmax = d_moon_earth * 6/5,
                            ymin = -d_moon_earth*(6/5 + 1/5)/2,
                            ymax = d_moon_earth*(6/5 + 1/5)/2,
                            n = 300)
    nxyz = nest(x, y, 0m)
    np = nest(x, y, clamp.(gfieldvals(nxyz, bodies)s², 0m, 0.5m))
    wrap_surface_and_countour(np, title = "4 Clamped gravitational acceleration")
end

# 3d space, plot gravity in the planet and moon plane,
# cap the plot at 0.05 m/s², fewer levels.
let
    bodies = bodies_3d()
    x,y = xygrids( xmin = -d_moon_earth/5,
                            xmax = d_moon_earth * 6/5,
                            ymin = -d_moon_earth*(6/5 + 1/5)/2,
                            ymax = d_moon_earth*(6/5 + 1/5)/2,
                            n = 300)
    nxyz = nest(x, y, 0m)
    np = nest(x, y, clamp.(gfieldvals(nxyz, bodies)s², 0m, 0.05m))
    wrap_surface_and_countour(np,
                title = "5 Clamped gravitational acceleration",
                levels = 20)
end
# 3d space, plot gravity in the planet and moon plane,
# cap the plot at 0.025 m/s², fewer levels.
# Zoom far out, both earth and moon. Zero gravity points visible.
let
    bodies = bodies_3d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 300)
    nxyz = nest(x, y, 0m)
    np = nest(x, y, clamp.(gfieldvals(nxyz, bodies)s², 0m, 0.025m))
    wrap_surface_and_countour(np,
                            title = "6 Clamped gravity",
                            levels = 20)
end

# 1d space, magnitude of g in one dimension, higher resolution.
# Zoom far out, both earth and moon. Zero gravity points visible,
# but hard to see.
let
    bodies = bodies_1d()
    n = 10000
    nx = nest(range(-d_moon_earth * 1.5, d_moon_earth * 1.5, length = n))
    gmagn = gfieldvals(nx, bodies)
    clf()
    plot(slice_1(nx)/m, gmagn*s²/m)
    grid()
    title("7 1d gravity, moon, earth and zero gravity points")
end

# 1d space, magnitude of g in one dimension, higher resolution.
# Zoom far out, both earth and moon. Zero gravity points visible,
# capped gravity at 0.025 m/s²
let
    bodies = bodies_1d()
    n = 1000
    nx = nest(range(-d_moon_earth * 1.5, d_moon_earth * 1.5, length = n))
    gmagn = clamp.(gfieldvals(nx, bodies)s², 0m, 0.025m)
    clf()
    plot(slice_1(nx)/m, gmagn/m)
    grid()
    title("8 1d gravity, moon, earth and zero gravity points")
    # Add a leader text
    dx = maximum(slice_1(nx)) - minimum(slice_1(nx))
    dy = maximum(gmagn) - minimum(gmagn)
    subsetx = slice_1(nx[1:floor(Int, n/2)])
    subsetg = gmagn[1:floor(Int, n/2)]
    y, xind = findmin(subsetg)
    x = subsetx[xind][1]
    annotate("0 m/s² at x = $(round(x/r_earth)) · r_earth",
        xy = [x; y]/m,
        xytext = [x + 0.1dx; y + 0.2dy]/m,
        xycoords = "data", arrowprops=Dict("facecolor"=>"blue"))
    # And one at the other minimum
    subsetx = slice_1(nx[floor(Int, n/2):end])
    subsetg = gmagn[floor(Int, n/2):end]
    y, xind = findmin(subsetg)
    x = subsetx[xind][1]
    annotate("0 m/s² at x = $(round(x/r_earth)) · r_earth",
        xy = [x; y]/m,
        xytext = [x - 0.3dx; y + 0.1dy]/m,
        xycoords = "data", arrowprops=Dict("facecolor"=>"blue"))
end
# 2d space, quiver plot gravity vectors in the planet and moon plane,
# axes are vectors not grids.
let
    bodies = bodies_2d()
    n = 20
    x = collect(range(-r_earth*1.5, d_moon_earth+ 1.45456e6m * 1.5, length = n))
    dy = maximum(x) -minimum(x)
    y = collect(range(-dy/2, dy/2, length = n))
    mg= fill([0.,0.]m/s², n, n)
    for i = 1:n, j = 1:n
        mg[i, j] = gvec([x[i], y[j]], bodies)
    end
    wrap_quiver(x, y, mg*s², title = "9 Sun, moon, earth")
end

# 2d space, quiver plot gravity vectors in the planet and moon plane,
# axes are vectors not grids. Cap vector length > 0.02.
let
    bodies = bodies_2d()
    n = 50
    x = collect(range(-d_moon_earth * 1.5, d_moon_earth * 1.5, length = n))
    dy = maximum(x) -minimum(x)
    y = collect(range(-dy/2, dy/2, length = n))
    mg = fill([0.,0.]m/s², n, n)
    for i = 1:n, j = 1:n
        g = gvec([x[i], y[j]], bodies)
        if len(g) < 0.02m/s²
            mg[i, j] = g
        else
            mg[i,j] = [NaN, NaN]m/s²
        end
    end
    wrap_quiver(x, y, mg*s², title = "10 Sun, moon, earth gravity")
end

# 2d space, stream plot gravity vectors in the planet and moon plane.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 50)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    wrap_stream_a(x, y, np, title = "11 Sun, moon, earth gravity")
end

# 2d space, stream plot gravity vectors in the planet and moon plane.
# Colour lines based on field strength.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 50)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    wrap_stream_b(x, y, np, title = "12 Sun, moon, earth gravity")
end
# 2d space, stream plot gravity vectors in the planet and moon plane.
# Line widths based on field strength.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 50)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    wrap_stream_c(x, y, np, title = "13 Sun, moon, earth gravity")
end

# 2d space, stream plot gravity vectors in the planet and moon plane.
# Line widths based on field strength.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 50)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    wrap_stream_d(x, y, np, title = "14 Sun, moon, earth gravity")
end

# 2d space, stream plot gravity vectors in the planet and moon plane.
# Cap vector length > 0.02. Much denser grid.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 150)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    np = map(np) do g
        if len(g) < 0.02m
            g
        else
            [NaN, NaN]m
        end
    end
    wrap_stream_a(x, y, np, title = "15 Sun, moon, earth gravity", density = 6.)
end

# 2d space, stream plot gravity vectors in the planet and moon plane.
# Cap vector length > 0.02. Much denser grid.
# Colour lines based on field strength.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 150)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    np = map(np) do g
        if len(g) < 0.02m
            g
        else
            [NaN, NaN]m
        end
    end
    wrap_stream_b(x, y, np, title = "16 Sun, moon, earth gravity", density = 6.)
end
# 2d space, stream plot gravity vectors in the planet and moon plane.
# Cap vector length > 0.02. Much denser grid.
# Line widths based on field strength.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 150)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    np = map(np) do g
        if len(g) < 0.02m
            g
        else
            [NaN, NaN]m
        end
    end
    wrap_stream_c(x, y, np, title = "17 Sun, moon, earth gravity", density = 3.)
end
# 2d space, stream plot gravity vectors in the planet and moon plane.
# Cap vector length > 0.02. Much denser grid.
# Line widths based on field strength.
# Colour lines based on field strength.
let
    bodies = bodies_2d()
    x, y = xygrids( xmin = -d_moon_earth * 1.5,
                    xmax =  d_moon_earth * 1.5,
                    ymin = -d_moon_earth * 1.5,
                    ymax = d_moon_earth * 1.5,
                    n = 150)
    nxy = nest(x, y)
    np = gvec(nxy, bodies)s²
    np = map(np) do g
        if len(g) < 0.02m
            g
        else
            [NaN, NaN]m
        end
    end
    wrap_stream_d(x, y, np, title = "18 Sun, moon, earth gravity", density = 3.)
end

#using ForwardDif

max_stars_repo_name: hustf/Gravityfield

max_issues_repo_name: hustf/Gravityfield

max_forks_repo_name: hustf/Gravityfield


text:
import PyPlot.pcolormesh
import PyPlot.colorbar

function pseudo_color_visualizer(ir_image)
    pcolormesh(sum(ir_image, 1)[1, end:-1:1, 1:end])
    colorbar()
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule cohomCalg_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("cohomCalg")
JLLWrappers.@generate_main_file("cohomCalg", UUID("5558cf25-a90e-53b0-b813-cadaa3ae7ade"))
end  # module cohomCalg_jll


text:
using Mill, Flux, Test, LearnBase
using Mill: LazyNode

@testset "LazyNode" begin
	ss = ["Hello world.", "Make peace.", "Make food.", "Eat penam."]

	function Mill.unpack2mill(ds::LazyNode{:Sentence})
		s = ds.data
		ss = map(x -> split(x, " "),s)
		x = NGramMatrix(reduce(vcat, ss), 3, 256, 2053)
		BagNode(ArrayNode(x), Mill.length2bags(length.(ss)))
	end

	@test LazyNode{:Sentence,Vector{String}}(ss).data == ss
	@test LazyNode(:Sentence, ss).data == ss
	@test LazyNode{:Sentence}(ss).data == ss
	@test nobs(LazyNode{:Sentence}(ss)) == 4
	@test (LazyNode{:Sentence}(ss))[[1,3]].data == ss[[1,3]]
	@test (LazyNode{:Sentence}(ss))[2].data == [ss[2]]
	@test reduce(catobs, [LazyNode{:Sentence}(ss)[i] for i in [1:2,3:4]]).data == ss
	@test reduce(catobs, [LazyNode{:Sentence}(ss)[i] for i in [1,2,3,4]]).data == ss

	ds = LazyNode{:Sentence}(ss)
	m = Mill.reflectinmodel(ds, d -> Dense(d,2), s -> SegmentedMeanMax(s))
	@test m(ds).data ≈ m.m(Mill.unpack2mill(ds)).data
end


text:
"""
    gpupowerdraw()::Float64

The function uses Linux `nvidia-smi` package to sample and get the average electricity
draw of the GPUs.
"""
function gpupowerdraw()
    gpucommand = `nvidia-smi`
    usage = Array{Any}(undef,60)
    cap = Array{Any}(undef,60)
    nogpus = 0

    for count in 1:60
        smis = Array{Any}[]
        smiss = Array{Any}[]
        gpus = Array{Any}[]
        powerdraw = Array{Float64}[]
        powercap = Array{Float64}[]
        smi = read(gpucommand, String);
        smi = split(smi, "\n")
        for s in smi
            push!(smis,split(s, " "))
        end
        for s in smis
            push!(smiss,filter(x->x≠"",s))
        end
        for strings in smiss
            if length(strings) > 5 && strings[6] == "/" && strings[10] == "/"
                push!(gpus,strings)
            end
        end
        nogpus = length(gpus)
        for g in gpus
            usagestr = ""
            capstr = ""
            if g[5] == "N/A"
                usagestr = "0W"
            else
                usagestr = usagestr * g[5]
            end
            if g[7] == "N/A"
                capstr = "0W"
            else
                capstr = capstr * g[7]
            end
            regexw = r"(\d+)"
            wattusg = match(regexw,usagestr)
            wattcap = match(regexw,capstr)
            powerdraw = vcat(powerdraw, parse(Float64,wattusg.match))
            powercap = vcat(powercap, parse(Float64,wattcap.match))
        end
        usage[count] = mean(powerdraw)
        cap[count] = mean(powercap)
        sleep(1)
    end
    return nogpus, mean(usage), mean(cap)
end

"""
    cpupowerdraw()::Float64

This function uses the Linux `powerstat` utility to get the average CPU energy cost.
"""
function cpupowerdraw()
    cpucommand = `powerstat -R -n -d0`
    cpu = read(cpucommand, String);
    cpu = split(cpu,"\n")
    cpu = cpu[66][60:64]
    return parse(Float64,cpu)
end

#TODO: further fine tune the model
"""
    rampowerdraw()::Float64

[Approximate RAM Power Draw](https://www.jedec.org/) the values are provided by the JEDEC we just take the
average of the total RAM with its Watt value to get the approximate power draw.
"""
function rampowerdraw(ramtype="DDR3")
    ramcommand = `free -m`
    powerused = Array{Float64}(undef,60)
    for count in 1:60
        ram = read(ramcommand, String);
        ram = split(ram,"\n")
        ram = split(ram[2]," ")
        filter!(x->x≠"",ram)
        totalram = parse(Float64,ram[2])
        if ramtype == "DDR3"
            powerused[count] = (totalram/1024)*0.3125
        elseif ramtype == "DDR2"
            powerused[count] = (totalram/1024)*0.625
        elseif ramtype == "DDR"
            powerused[count] = (totalram/1024)*6.5
        else
            error("$ramtype unrecognized RAM type.")
        end
        sleep(1)
    end
    return mean(powerused)
end

#TODO: modify the code to work in Windows.
"""
    avgpowerdraw()::Float64

[Average Power Draw](https://arxiv.org/abs/1906.02243) where `pc` is the average power
draw (in watts) from all CPU sockets during training, let `pr` be the average power draw from all
DRAM (main memory) sockets, let `pg` be the average power draw of the GPUs during training and `g`
the number of available gpus.

`apd = 1.58*t*(pc + pr + g*pg)/1000`

returns the average power consumption in kWh.

By default it assumes you use `"DDR3"` memory but you can pass `"DDR2"` or `"DDR"` to get a better estimate.
"""
function avgpowerdraw(freeram="DDR3")
    g, pg, pc, pr = 0.0, 0.0, 0.0, 0.0
    starttime = time()
    try
        g, pg, _ = gpupowerdraw()
    catch ex
        println(ex.msg)
        return 0.0
    end
    try
        pc = cpupowerdraw()
    catch ex
        println(ex.msg)
        return 0.0
    end
    try
        pr = rampowerdraw(freeram)
    catch ex
        println(ex.msg)
        return 0.0
    end
    endtime = time()
    elapsedtime = (endtime - starttime)*0.0002777778
    return 1.58*elapsedtime*(pc + pr + g*pg)/1000
end


text:
"""
    clockwise(line::Union{LineString, Vector{Position}})::Bool

Take a ring and return true or false whether or not the ring is clockwise or counter-clockwise.

# Examples
```jldoctest
julia> line = LineString([[0, 0], [1, 1], [1, 0], [0, 0]])
LineString(Array{Float64,1}[[0.0, 0.0], [1.0, 1.0], [1.0, 0.0], [0.0, 0.0]])

julia> clockwise(line)
true
```
"""
function clockwise(line::Union{LineString, Vector{Position}})::Bool

    let ring end
    # geotype fails with Vector{Vector{...}}
    if geotype(line) === :LineString
        ring = line.coordinates
    else
        ring = line[1]
    end

    sum = 0
    i = 1

    prev = []
    cur = []
    while i < length(ring)
        prev = isempty(cur) ? ring[1] : cur
        cur = ring[i]
        sum += (cur[1] - prev[1]) * (cur[2] + prev[2])
        i += 1
    end

    return sum > 0
end

"""
    concave(poly::Polygon)::Bool

Take a polygon and return true or false as to whether it is concave or not.

# Examples
```jldoctest
julia> poly = Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]])
Polygon(Array{Array{Float64,1},1}[[[0.0, 0.0], [0.0, 1.0], [1.0, 1.0], [1.0, 0.0], [0.0, 0.0]]])

julia> concave(poly)
false
```
"""
function concave(poly::Polygon)::Bool
    coords = poly.coordinates

    length(coords[1]) <= 4 && return false

    sign = false
    n = length(coords[1]) - 1

    for i in 1:n
        j = ((i + 1) % n) === 0 ? 1 : (i + 1) % n
        m = ((i + 2) % n) === 0 ? 1 : (i + 2) % n

        dx1 = coords[1][m][1] - coords[1][j][1]
        dy1 = coords[1][m][2] - coords[1][j][2]
        dx2 = coords[1][i][1] - coords[1][j][1]
        dy2 = coords[1][i][2] - coords[1][j][2]

        cross = (dx1 * dy2) - (dy1 * dx2)

        if i === 0
            sign = cross > 0
        elseif sign !== (cross > 0)
            return true
        end
    end

    return false
end

function equal(geo1::T, geo2::T) where {T <: AbstractGeometry}
    geotype(geo1) !== geotype(geo2) && return false

    geotype(geo1) === :Point && return comparePoints(geo1.coordinates, geo2.coordinates)
    geotype(geo1) === :LineString && return compareLines(geo1.coordinates, geo2.coordinates)
end

function comparePoints(p1::Position, p2::Position)
    length(p1) !== length(p2) && return false

    for i in eachindex(p1)
        round(p1[i]; digits=10) !== round(p2[i]; digits=10) && return false
    end

    return true
end

function compareLines(p1::Vector{Position}, p2::Vector{Position})
    # TODO: complete this
    length(p1[1]) !== length(p2[1]) && return false
end

"""
    parallel(line1::LineString, line2::LineString)::Bool

Return `true` if each segment of `line1` is parallel to the correspondent segment of `line2`

# Examples
```jldoctest
julia> line1 = LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]])
LineString(Array{Float64,1}[[9.17036, 45.478], [9.16443, 45.4826], [9.16664, 45.484]])

julia> line2 = LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]])
LineString(Array{Float64,1}[[9.16936, 45.478], [9.16343, 45.4826], [9.16564, 45.484]])

julia> parallel(line1, line2)
true
```
"""
function parallel(line1::LineString, line2::LineString)::Bool
    seg1 = linesegment(line1)
    seg2 = linesegment(line2)

    for i in eachindex(seg1)
        coors2 = nothing
        coors1 = seg1[i].coordinates

        try
            coors2 = seg2[i].coordinates
        catch e
            isa(e, BoundsError) && break

        end
        isParallel(coors1, coors2) == false && return false
    end

    return true
end

@inline function isParallel(p1::Vector{Position}, p2::Vector{Position})
    slope1 = bearing_to_azimuth(rhumb_bearing(p1[1], p2[1]))
    slope2 = bearing_to_azimuth(rhumb_bearing(p1[2], p2[2]))

    return slope1 === slope2
end

"""
    point_on_line(point::Point, line::LineString, ignoreEndVertices::Bool=false)::Bool

Return true if a point is on a line. Accept a optional parameter to ignore the
start and end vertices of the linestring.

# Examples
```jldoctest
julia> point = Point([1,1])
Point([1.0, 1.0])

julia> line = LineString([[0, 0], [3, 3], [4, 4]])
LineString(Array{Float64,1}[[0.0, 0.0], [3.0, 3.0], [4.0, 4.0]])

julia> point_on_line(point, line)
true
```
"""
function point_on_line(point::Point, line::LineString, ignoreEndVertices::Bool=false)::Bool
    pCoords = point.coordinates
    lCoords = line.coordinates

    ignore = "none"
    for i in 1:length(lCoords) - 1
        if ignoreEndVertices == true
            i === 1 && (ignore = "start")
            i === length(lCoords) - 2 && (ignore = "end")
            (i === 1 && i + 1 === length(lCoords) - 1) && (ignore = "both")
        end
        isPointOnSegment(lCoords[i], lCoords[i + 1], pCoords, ignore) && return true
    end
    return false
end

@inline function isPointOnSegment(start::Position, stop::Position, coords::Position, excludeBoundary::String="none")::Bool
    x, y = coords
    x1, y1 = start
    x2, y2 = stop

    dxc = x - x1
    dyc = y - y1
    dx1 = x2 - x1
    dy1 = y2 - y1

    cross = dxc * dy1 - dyc * dx1
    cross != 0 && return false

    if excludeBoundary === "none"
        if abs(dx1) >= abs(dy1)
            return dx1 > 0 ? x1 <= x && x <= x2 : x2 <= x && x <= x1
        end
        return dy1 > 0 ? y1 <= y && y <= y2 : y2 <= y && y <= y1
    elseif excludeBoundary === "start"
        if abs(dx1) >= abs(dy1)
             return dx1 > 0 ? x1 < x && x <= x2 : x2 <= x && x < x1
        end
        return dy1 > 0 ? y1 < y && y <= y2 : y2 <= y && y < y1
    elseif excludeBoundary === "end"
        if abs(dx1) >= abs(dy1)
            return dx1 > 0 ? x1 <= x && x < x2 : x2 < x && x <= x1
        end
        return dy1 > 0 ? y1 <= y && y < y2 : y2 < y && y <= y1
    elseif excludeBoundary === "both"
        if abs(dx1) >= abs(dy1)
            return dx1 > 0 ? x1 < x && x < x2 : x2 < x && x < x1
        end
        return dy1 > 0 ? y1 < y && y < y2 : y2 < y && y < y1
    end
    return false
end

"""
    point_in_polygon(point::Point, polygon::Union{Polygon, MultiPolygon}, ignoreBoundary::Bool=false)::Bool

Take a Point and a Polygon and determine if the point
resides inside the polygon. The polygon can be convex or concave. The function accounts for holes.

# Examples
```jldoctest
julia> point = Point([-77, 44])
Point([-77.0, 44.0])

julia> poly = Polygon([[ [-81, 41], [-81, 47], [-72, 47], [-72, 41], [-81, 41]]])
Polygon(Array{Array{Float64,1},1}[[[-81.0, 41.0], [-81.0, 47.0], [-72.0, 47.0], [-72.0, 41.0], [-81.0, 41.0]]])

julia> point_in_polygon(point, poly)
true
```
"""
function point_in_polygon(point::Point, polygon::Union{Polygon, MultiPolygon}, ignoreBoundary::Bool=false)::Bool

    pt = point.coordinates
    poly = polygon.coordinates

    inBBox(pt, bbox(polygon)) == false && return false

    geotype(polygon) === :Polygon && (poly = [poly])

    inside = false
    for i in eachindex(poly)
        if inRing(pt, poly[i][1], ignoreBoundary)
            inHole = false
            k = 1

            while k < length(poly[i]) && !inHole
                inRing(pt, poly[i][k], !ignoreBoundary) == true && (inHole = true)
                k += 1
            end

            !inHole && (inside = true)
        end
    end
    return inside
end

function inRing(pt::Position, ring::Vector{Position}, ignoreBoundary::Bool=false)
    inside = false

    (ring[1][1] == ring[length(ring) - 1][1] && ring[1][2] == ring[length(ring) - 1][1]) && (ring = ring[1, length(ring) - 1])

    for i in 1:length(ring) - 1
        j = i + 1

        xi = ring[i][1]
        yi = ring[i][2]
        xj = ring[j][1]
        yj = ring[j][2]

        onBoundary = (pt[2] * (xi - xj) + yi * (xj - pt[1]) + yj * (pt[1] - xi) == 0) &&
            ((xi - pt[1]) * (xj - pt[1]) <= 0) && ((yi - pt[2]) * (yj - pt[2]) <= 0)

        onBoundary && return !ignoreBoundary

        intersect = ((yi > pt[2]) !== (yj > pt[2])) && (pt[1] < (xj - xi) * (pt[2] - yi) / (yj - yi) + xi)

        intersect && (inside =  !inside)
    end

    return inside
end

function inBBox(pt::Position, bbox::Vector{Float64})
    return bbox[1] <= pt[1] &&  bbox[2] <= pt[2] &&
        bbox[3] >= pt[1] && bbox[4] >= pt[2]
end

"""
    contains(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool

Return true if the second geometry is completely contained by the first geometry.
The interiors of both geometries must intersect and, the interior and boundary of the secondary (geometry b)
must not intersect the exterior of the primary (geometry a).
`contains` returns the exact opposite result of `within`.

# Examples
```jldoctest
julia> line = LineString([[1, 1], [1, 2], [1, 3], [1, 4]])
LineString(Array{Float64,1}[[1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0]])

julia> point = Point([1, 2])
Point([1.0, 2.0])

julia> contains(line, point)
true
```
"""
function contains(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool
    type1 = geotype(ft1)
    type2 = geotype(ft2)

    coords1 = ft1.coordinates
    coords2 = ft2.coordinates

    if type1 === :Point
        if type2 === :Point
            return coords1[1] == coords2[1] && coords1[2] == coords2[2]
        end
        throw(error("$(type2) is not a supported type."))
    elseif type1 === :LineString
        if type2 === :Point
            return point_on_line(ft2, ft1, true)
        elseif type2 === :LineString
            return lineOnLine(ft1, ft2)
        else
            throw(error("$(type2) is not a supported type."))
        end

    elseif type1 === :Polygon
        if type2 === :Point
            return point_in_polygon(ft2, ft1, true)
        elseif type2 === :LineString
            return lineInPolygon(ft1, ft2)
        elseif type2 === :Polygon
            return polygonInPolygon(ft1, ft2)
        else
            throw(error("$(type2) is not a supported type."))
        end
    else
        throw(error("$(type1) is not a supported type."))
    end
end

contains(ft1::Feature, ft2::Feature) = contains(ft1.geometry, ft2.geometry)

function lineInPolygon(poly::Polygon, line::LineString)
    out = false

    polybox = bbox(poly)
    linebox = bbox(line)

    !(bboxOverlap(polybox, linebox)) && return false

    coords = line.coordinates

    for i in 1:length(coords) - 1
        mid = [(coords[i][1] + coords[i + 1][1]) / 2, (coords[i][2] + coords[i + 1][2]) / 2]
        if point_in_polygon(Point(mid), poly, true)
            out = true
            break
        end
    end
    return out
end

function lineOnLine(line1::LineString, line2::LineString)
    for i in eachindex(line1.coordinates)
        !(point_on_line(line1.coordinates[i], line2)) && return false
    end
    return true
end

function polygonInPolygon(ft1::Polygon, ft2::Polygon, reverse::Bool=false)
    polybox1 = bbox(ft1)
    polybox2 = bbox(ft2)
    coords = []

    if reverse
        coords = ft1.coordinates
        !(bboxOverlap(polybox2, polybox1)) && return false

        for ring in coords
            for coord in ring
                !(point_in_polygon(Point(coord), ft2)) && return false
            end
        end
    else
        coords = ft2.coordinates
        !(bboxOverlap(polybox1, polybox2)) && return false

        for ring in coords
            for coord in ring
                !(point_in_polygon(Point(coord), ft1)) && return false
            end
        end
    end

    return true
end

function bboxOverlap(box1::Vector{T}, box2::Vector{T}) where {T <: Real}
    box1[1] > box2[1] && return false
    box1[3] < box2[3] && return false
    box1[2] > box2[2] && return false
    box1[4] < box2[4] && return false
    return true
end

"""
    within(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool

Return `true` if the first geometry is completely within the second geometry.
The interiors of both geometries must intersect and, the interior and boundary of the primary (geometry a)
must not intersect the exterior of the secondary (geometry b).
`within` returns the exact opposite result of `contains`.

# Examples
```jldoctest
julia> line = LineString([[1, 1], [1, 2], [1, 3], [1, 4]])
LineString(Array{Float64,1}[[1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0]])

julia> point = Point([1, 2])
Point([1.0, 2.0])

julia> within(point, line)
true
```
"""
function within(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool
    type1 = geotype(ft1)
    type2 = geotype(ft2)

    coords1 = ft1.coordinates
    coords2 = ft2.coordinates

    if type1 === :Point
        if type2 === :LineString
            return point_on_line(ft1, ft2, true)
        elseif type2 === :Polygon
            return point_in_polygon(ft1, ft2, true)
        end
        throw(error("$(type2) is not a supported type."))
    elseif type1 === :LineString
        if type2 === :Polygon
            return lineInPolygon(ft1, ft2)
        elseif type2 === :LineString
            return lineOnLine(ft1, ft2)
        else
            throw(error("$(type2) is not a supported type."))
        end

    elseif type1 === :Polygon
        if type2 === :Polygon
            return polygonInPolygon(ft1, ft2, true)
        else
            throw(error("$(type2) is not a supported type."))
        end
    else
        throw(error("$(type1) is not a supported type."))
    end
end

function lineInPolygon(line::LineString, poly::Polygon)
    polybox = bbox(poly)
    linebox = bbox(line)

    !(bboxOverlap(polybox, linebox)) && return false

    coords = line.coordinates
    inside = false

    for i in 1:length(coords) - 1
        !(point_in_polygon(Point(coords[i]), poly)) && return false
        !inside && (inside = point_in_polygon(Point(coords[i]), poly, true))
        if !inside
            mid = [(coords[i][1] + coords[i + 1][1]) / 2, (coords[i][2] + coords[i + 1][2]) / 2]
            inside = point_in_polygon(Point(mid), poly, true)
        end
    end
    return inside
end

"""
    disjoint(geom1::AbstractGeometry, geom2::AbstractGeometry)::Bool

Return `true` if the intersection of the two geometries is an empty set.

# Examples
```jldoctest
julia> poly = Polygon([[[-1, 2], [3, 2], [3, 3], [-1, 3], [-1, 2]]])
Polygon(Array{Array{Float64,1},1}[[[-1.0, 2.0], [3.0, 2.0], [3.0, 3.0], [-1.0, 3.0], [-1.0, 2.0]]])

julia> point = Point([1, 1])
Point([1.0, 1.0])

julia> disjoint(poly, point)
true
```
"""
function disjoint(geom1::AbstractGeometry, geom2::AbstractGeometry)::Bool
    type1 = geotype(geom1)
    type2 = geotype(geom2)

    coords1 = geom1.coordinates
    coords2 = geom2.coordinates

    if type1 === :Point
        if type2 === :Point
            return !(coords1[1] == coords2[1] && coords1[2] == coords2[2])
        elseif type2 === :LineString
            return !point_on_line(geom1, geom2)
        elseif type2 === :Polygon
            !point_in_polygon(geom1, geom2)
        end

    elseif type1 === :LineString
        if type2 === :Point
            return !point_on_line(geom2, geom1)
        elseif type2 === :LineString
            return !lineOnLine(geom1, geom2)
        elseif type2 === :Polygon
            return !lineInPolygon(geom2, geom1)
        end
    elseif type1 === :Polygon
        if type2 === :Point
            return !point_in_polygon(geom2, geom1)
        elseif type2 === :LineString
            return !lineInPolygon(geom1, geom2)
        elseif type2 === :Polygon
            return !polyInPoly(geom2, geom1)
        end
    end
    return false
end

disjoint(geom1::AbstractFeature, geom2::AbstractFeature) = disjoint(geom1.geometry, geom2.geometry)

function polyInPoly(poly1::Polygon, poly2::Polygon)
    coords1 = poly1.coordinates
    coords2 = poly2.coordinates

    for ring in coords1
        for coord in ring
            (point_in_polygon(Point(coord), poly2)) && return true
        end
    end

    for ring in coords2
        for coord in ring
            (point_in_polygon(Point(coord), poly1)) && return true
        end
    end

    inter = line_intersects(polygon_to_line(poly1), polygon_to_line(poly2))
    inter != nothing && return true

    return false

end

"""
    line_intersects(line1::AbstractLineString, line2::AbstractLineString)

Find a point that intersects LineStrings with two coordinates each.

# Examples
```jldoctest
julia> line1 = LineString([[124.584961,-12.768946],[126.738281,-17.224758]])
LineString(Array{Float64,1}[[124.585, -12.7689], [126.738, -17.2248]])

julia> line2 = LineString([[123.354492,-15.961329],[127.22168,-14.008696]])
LineString(Array{Float64,1}[[123.354, -15.9613], [127.222, -14.0087]])

julia> line_intersects(line1, line2)
Point([125.584, -14.8357])
```
"""
function line_intersects(line1::AbstractLineString, line2::AbstractLineString)
    coords1 = line1.coordinates
    coords2 = line2.coordinates

    (length(coords1) == 2 && length(coords2) == 2) && return intersects(line1, line2)

    result = []

    for i in 1:length(coords1) - 1
        for j in 1:length(coords2) - 1
            inter = intersects(LineString([coords1[i], coords1[i + 1]]), LineString([coords2[j], coords2[j + 1]]))
            inter != nothing && push!(result, inter.coordinates)
        end
    end
    return unique!(result)
end

function intersects(line1::AbstractLineString, line2::AbstractLineString)
    coords1 = line1.coordinates
    coords2 = line2.coordinates

    (length(coords1) != 2 || length(coords2) != 2) && throw(error("Lines must contain only 2 coordinates."))

    x1, y1 = coords1[1]
    x2, y2 = coords1[2]
    x3, y3 = coords2[1]
    x4, y4 = coords2[2]

    d = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1))
    a = ((x4 - x3) * (y1 - y3)) - ((y4 - y3) * (x1 - x3))
    b = ((x2 - x1) * (y1 - y3)) - ((y2 - y1) * (x1 - x3))

    if d == 0
        if a == 0 && b == 0
            return nothing
        end
        return nothing
    end

    ã = a / d
    b̃ = b / d

    if ã >= 0 && ã <= 1 && b̃ >= 0 && b̃ <= 1
        x = x1 + (ã * (x2 - x1))
        y = y1 + (ã * (y2 - y1))
        return Point([x, y])
    end

    return nothing
end

"""
     crosses(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool

Return `true` if the intersection results in a geometry whose dimension is one less than
the maximum dimension of the two source geometries and the intersection set is interior to
both source geometries.

# Examples
```jldoctest
julia> line = LineString([[1, 1], [1, 2], [1, 3], [1, 4]])
LineString(Array{Float64,1}[[1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0]])

julia> line2 = LineString([[-2, 2], [4, 2]])
LineString(Array{Float64,1}[[-2.0, 2.0], [4.0, 2.0]])

julia> crosses(line2, line)
true
```
"""
function crosses(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool
    type1 = geotype(ft1)
    type2 = geotype(ft2)

    if type1 === :MultiPoint
        if type2 === :LineString
            return MpCrossLs(ft1, ft2)
        elseif type2 === :Polygon
            return MpCrossPoly(ft1, ft2)
        else
            throw(error("Geometry $(type2) is not supported."))
        end
    elseif type1 === :LineString
        if type2 === :MultiPoint
            return MpCrossLs(ft2, ft1)
        elseif type2 === :Polygon
            return LsCrossPoly(ft1, ft2)
        elseif type2 === :LineString
            return LsCrossLs(ft1, ft2)
        else
            throw(error("Geometry $(type2) is not supported."))
        end
    elseif type1 === :Polygon
        if type2 === :MultiPoint
            return MpCrossPoly(ft2, ft1)
        elseif type2 === :LineString
            return LsCrossPoly(ft2, ft1)
        else
            throw(error("Geometry $(type2) is not supported."))
        end
    end

    throw(error("Geometry $(type1) is not supported."))
end

crosses(ft1::Feature, ft2::Feature) = crosses(ft1.geometry, ft2.geometry)

function MpCrossLs(geom1::MultiPoint, geom2::LineString)
    intPoint = false
    extPoint = false

    pLength = length(geom1.coordinates)
    i = 1

    while i < pLength && !intPoint && !extPoint

        for j in 1:length(geom2.coordinates) - 1
            incVertices = "both"

            (j === 1 || j === length(geom2.coordinates) - 2) && (incVertices = "none")

            if isPointOnSegment(geom2.coordinates[j], geom2.coordinates[j + 1], geom1.coordinates[i], incVertices)
                intPoint = true
            else
                extPoint = true

            end

        end
        i += 1
    end

    return intPoint && extPoint
end

function LsCrossLs(line1::LineString, line2::LineString)
    inter = line_intersects(line1, line2)

    if length(inter) > 0
        for i in 1:length(line1.coordinates) - 1
            for j in 1:length(line2.coordinates) - 1
                incVertices = "both"

                (j === 1 || j === length(line2.coordinates) - 2) && (incVertices = "none")

                isPointOnSegment(line1.coordinates[i], line1.coordinates[i + 1], line2.coordinates[j], incVertices) && return true
            end
        end
    end
    return false
end

function LsCrossPoly(line::LineString, poly::Polygon)
    line2 = polygon_to_line(poly)
    inter = line_intersects(line, line2)

    (length(inter) > 0) && return true

    return false
end

function MpCrossPoly(mp::MultiPoint, poly::Polygon)
    intPoint = false
    extPoint = false
    pLength = length(mp.coordinates[1])
    i = 1

    while i < pLength && intPoint && extPoint
        if point_in_polygon(Point(mp.coordinates[1][i]), poly)
            intPoint = true
        else
            extPoint = true
        end
        i += 1
    end
    return intPoint && extPoint
end

"""
    overlap(geom1::AbstractGeometry, geom2::AbstractGeometry)::Bool

Compare two Geometries of the same dimension and return true if their intersection set results in a geometry
different from both but of the same dimension. It applies to Polygon/Polygon, LineString/LineString,
Multipoint/Multipoint, MultiLineString/MultiLineString and MultiPolygon/MultiPolygon.

# Examples
```jldoctest
julia> poly1 = Polygon([[[0,0],[0,5],[5,5],[5,0],[0,0]]])
Polygon(Array{Array{Float64,1},1}[[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]])

julia> poly2 = Polygon([[[1,1],[1,6],[6,6],[6,1],[1,1]]])
Polygon(Array{Array{Float64,1},1}[[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]])

julia> overlap(poly1, poly2)
true
```
"""
function overlap(geom1::AbstractGeometry, geom2::AbstractGeometry)::Bool
    # TODO: add support for LineStrings and ML

    type1 = geotype(geom1)
    type2 = geotype(geom2)

    !isequal(type1, type2) && throw(error("Gometries must be of the same type."))
    type1 === :Point && throw(error("$(type1) geometry is not supported."))

    if length(geom1.coordinates) === length(geom2.coordinates)
        geom1.coordinates ≈ geom2.coordinates && return false
    end

    overlap = 0
    coords1 = geom1.coordinates
    coords2 = geom2.coordinates

    if isequal(type1, :MultiPoint)
        for i in eachindex(coords1)
            for j in eachindex(coords2)
                (coords1[i][1] == coords2[j][1] && coords1[i][2] == coords2[j][2]) && (overlap += 1)
            end
        end
    elseif isequal(type1, :Polygon)
        line1 = polygon_to_line(geom1)
        line2 = polygon_to_line(geom2)

        (line_intersects(line1, line2) != []) && (overlap += 1)

    elseif isequal(type1, :MultiPolygon)
        for poly1 in coords1
            for poly2 in coords2
                line1 = polygon_to_line(poly1)
                line2 = polygon_to_line(poly2)

                (line_intersects(line1, line2) != []) && (overlap += 1)
            end
        end
    else
        throw(error("Unsupported geometry."))
    end

    return overlap > 0
end

"""
    overlap(feat1::AbstractFeature, feat2::AbstractFeature)::Bool

Compare two Features of the same dimension and return true if their intersection set results in a geometry
different from both but of the same dimension. It applies to Polygon/Polygon, LineString/LineString,
Multipoint/Multipoint, MultiLineString/MultiLineString and MultiPolygon/MultiPolygon.
"""
overlap(feat1::AbstractFeature, feat2::AbstractFeature)::Bool = overlap(feat1.geometry, feat2.geometry)


text:
EMPTY_ARRAY = spzeros(0,1)*[0.0];

#function get_netlib_problem(file_name::String)
#	file = matopen(file_name)
#	A= sparse(read(file, "A"))
#	b = read(file,"b")
#	b = b[:]
#	c = read(file,"c")
#	c = c[:];
#
#	return A, b, c
#end

function get_netlib_problem(file_name::String)
  dir = dirname(@__FILE__) * "/Problems/"
  return get_netlib_problem(dir, file_name)
end

using MAT
function get_netlib_problem(dir::String, file_name::String)
  #= USE JSON??
  @mput dir file_name
	@matlab begin
    cd(dir)
    load(file_name)
  end
  @mget A b c
  =#

  dir = dirname(@__FILE__)
  file = matopen(dir * "/Problems/" * file_name)

  A = read(file, "A") # note that this does NOT introduce a variable ``varname`` into scope
  b = read(file,"b")[:]
  c = read(file,"c")[:]

  close(file)

  return A, b, c
end

function test_homogeneous_algorithm(A::SparseMatrixCSC{Float64,Int64}, b::Array{Float64,1}, c::Array{Float64,1}, settings::class_settings)
	return test_homogeneous_algorithm(A, b, c, spzeros(length(c),length(c)), settings);
end

function test_homogeneous_algorithm(A::SparseMatrixCSC{Float64,Int64}, b::Array{Float64,1}, c::Array{Float64,1}, Q::SparseMatrixCSC{Float64,Int64}, settings::class_settings)
	try
		qp = class_quadratic_program(A, b, c, Q);

		vars = class_variables(n(qp),m(qp));
		if settings.verbose
			println(n(qp) , " variables, ", m(qp) , " constraints")
		end

		return homogeneous_algorithm(qp, vars, settings)
	catch e
		println("ERROR in test_homogeneous_algorithm")
		throw(e)
	end
end

function tridiagonal(k,di,offdi)
	Q = spzeros(k,k)
	for j = 1:(k-1)
		Q[j,j+1] = offdi
		Q[j,j] = di
		Q[j+1,j] = offdi
	end

	return Q;
end

using JuMP, Ipopt

function solve_with_JuMP(A::SparseMatrixCSC{Float64,Int64}, b::Array{Float64,1}, c::Array{Float64,1}, solver=IpoptSolver(max_iter=300))
	solve_with_JuMP(A, b, c, spzeros(length(c),length(c)), solver)
end

function solve_with_JuMP(A::SparseMatrixCSC{Float64,Int64}, b::Array{Float64,1}, c::Array{Float64,1}, Q::SparseMatrixCSC{Float64,Int64}, solver=IpoptSolver(max_iter=300))
	model = Model(solver=solver)

	n, k = size(A)

	@defVar(model, x[1:k] >= 0)
	@setObjective( model, Min, sum{c[i]*x[i], i=1:k} + sum{0.5*x[j]*Q[i,j]*x[i], i=1:k, j=1:k})

	#@addConstraint(model, A*x .== b)
	@defConstrRef constr[1:n]
	for i = 1:n
		constr[i] = @addConstraint(model, sum{A[i,j]*x[j], j=1:k} == b[i])
	end

	status = solve(model)

	if status == :Optimal
		return 1, getValue(x)[:], getDual(constr)[:]
	elseif status == :Infeasible
		return 2
	elseif status == :Unbounded
		return 3
	else
		return 0
	end
end

function trivial_test(A::SparseMatrixCSC{Float64,Int64}, b::Array{Float64,1}, c::Array{Float64,1}, correct_status::Symbol, problem_name = "", settings::class_settings = class_settings(), verbose = false)
	return trivial_test(A, b, c, spzeros(length(c),length(c)), correct_status, problem_name, settings, verbose)
end

function trivial_test(A::SparseMatrixCSC{Float64,Int64}, b::Array{Float64,1}, c::Array{Float64,1}, Q::SparseMatrixCSC{Float64,Int64}, correct_status::Symbol, problem_name = "", settings::class_settings = class_settings(), verbose = false)
	try
		#settings.linear_system_solver = linear_solver_MATLAB();
		#settings.linear_system_solver.options.sym = 0;
		settings.verbose = verbose;

		status, vars, iterations = test_homogeneous_algorithm(A, b, c, Q, settings);

		if status == correct_status
			println("Problem ", problem_name, " solved successfully in ", iterations, " iterations")
		else
			println("ERROR - Problem ", problem_name, " not solved sucessfully.")
			println("Expected status ", correct_status, " got status ", status)
			#println("x_scaled = ", vars.x()/vars.tau())
		end

		return status, vars, iterations
	catch e
		println("ERROR in trivial_test")
		println("===================================")
		throw(e)
	end
end

max_forks_repo_name: ohinder/Everest


text:
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Quadrature,Cubature,Cuba
using Plots

# Physical and numerical parameters (fixed)
nu = 0.07
nx = 10001 #101
x_max = 2.0 * pi
dx = x_max / (nx - 1.0)
nt = 2 #10
dt = dx * nu
t_max = dt * nt

# Analytic function
analytic_sol_func(t, x) = -2*nu*(-(-8*t + 2*x)*exp(-(-4*t + x)^2/(4*nu*(t + 1)))/
                          (4*nu*(t + 1)) - (-8*t + 2*x - 12.5663706143592)*
                          exp(-(-4*t + x - 6.28318530717959)^2/(4*nu*(t + 1)))/
                          (4*nu*(t + 1)))/(exp(-(-4*t + x - 6.28318530717959)^2/
                          (4*nu*(t + 1))) + exp(-(-4*t + x)^2/(4*nu*(t + 1)))) + 4

function burgers(strategy, minimizer)

    @parameters x t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2

    eq = Dt(u(x, t)) + u(x, t) * Dx(u(x, t))  ~ nu * Dxx(u(x, t))

    bcs = [u(x,0.0) ~ analytic_sol_func(x, 0.0),
           u(0.0,t) ~ u(x_max,t)]

    domains = [x ∈ IntervalDomain(0.0, x_max),
               t ∈ IntervalDomain(0.0, t_max)]

    chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
    discretization = PhysicsInformedNN(chain, strategy)

    indvars = [x, t]   #physically independent variables
    depvars = [u]      #dependent (target) variable

    dim = length(domains)

    losses = []
    error = []
    times = []

    dx_err = 0.00005

    error_strategy = GridTraining(dx_err)

    initθ = Float64.(DiffEqFlux.initial_params(chain))
    eltypeθ = eltype(initθ)
    parameterless_type_θ = DiffEqBase.parameterless_type(initθ)
    phi = NeuralPDE.get_phi(chain,parameterless_type_θ)
    derivative = NeuralPDE.get_numeric_derivative()

    _pde_loss_function = NeuralPDE.build_loss_function(eq,indvars,depvars,
                                             phi,derivative,nothing,chain,initθ,error_strategy)

    bc_indvars = NeuralPDE.get_variables(bcs,indvars,depvars)
    _bc_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars,
                                              phi,derivative,nothing,chain,initθ,error_strategy,
                                              bc_indvars = bc_indvar) for (bc,bc_indvar) in zip(bcs,bc_indvars)]

    train_sets = NeuralPDE.generate_training_sets(domains,dx_err,[eq],bcs,eltypeθ,indvars,depvars)
    train_domain_set, train_bound_set = train_sets

    pde_loss_function = NeuralPDE.get_loss_function(_pde_loss_function,
                                          train_domain_set[1],eltypeθ,
                                          parameterless_type_θ,error_strategy)

    bc_loss_functions = [NeuralPDE.get_loss_function(loss,set,
                                                 eltypeθ, parameterless_type_θ,
                                                 error_strategy) for (loss, set) in zip(_bc_loss_functions,train_bound_set)]

    loss_functions = [pde_loss_function; bc_loss_functions]
    loss_function__ = θ -> sum(map(l->l(θ) ,loss_functions))

    function loss_function_(θ,p)
            return loss_function__(θ)
    end

    cb = function (p,l)

        timeCounter = 0.0
        deltaT_s = time_ns() #Start a clock when the callback begins, this will evaluate questo misurerà anche il calcolo degli uniform error

        ctime = time_ns() - startTime - timeCounter #This variable is the time to use for the time benchmark plot
        append!(times, ctime/10^9) #Conversion nanosec to seconds
        append!(losses, l)
        append!(error, loss_function__(p))
        #println(length(losses), " Current loss is: ", l, " uniform error is, ", loss_function__(p))

        timeCounter = timeCounter + time_ns() - deltaT_s #timeCounter sums all delays due to the callback functions of the previous iterations

        return false
    end

    @named pde_system = PDESystem(eq, bcs, domains, indvars, depvars)

    discretization = NeuralPDE.PhysicsInformedNN(chain, strategy, init_params =initθ)
    prob = NeuralPDE.discretize(pde_system,discretization)

    startTime = time_ns() #Fix initial time (t=0) before starting the training

    if minimizer == "both"
        res = GalacticOptim.solve(prob, ADAM(); cb = cb, maxiters=5)
        prob = remake(prob,u0=res.minimizer)
        res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=15)
    else
        res = GalacticOptim.solve(prob, minimizer; cb = cb, maxiters=500)
    end

    phi = discretization.phi

    params = res.minimizer

    return [error, params, times, losses]
end

# Settings:
#maxIters = [(0,0,0,0,0,0,20000),(300,300,300,300,300,300,300)] #iters

strategies = [NeuralPDE.QuadratureTraining()]

strategies_short_name = ["QuadratureTraining"]

minimizers = [ADAM(),
              ADAM(0.000005),
              ADAM(0.0005),
              RMSProp(),
              RMSProp(0.00005),
              RMSProp(0.05),
              BFGS(),
              LBFGS()]

minimizers_short_name = ["ADAM",
                         "ADAM(0.000005)",
                         "ADAM(0.0005)",
                         "RMS",
                         "RMS(0.00005)",
                         "RMS(0.05)",
                         "BFGS",
                         "LBFGS"]

# Run models
error_res =  Dict()
params_res = Dict()
times = Dict()
losses_res = Dict()

print("Starting run \n")

for min in 1:length(minimizers) # minimizer
      for strat in 1:length(strategies) # strategy
            #println(string(strategies_short_name[1], "  ", minimizers_short_name[min]))
            res = burgers(strategies[strat], minimizers[min])
            push!(error_res, string(strat,min)     => res[1])
            push!(params_res, string(strat,min) => res[2])
            push!(times, string(strat,min)        => res[3])
            push!(losses_res, string(strat,min)        => res[4])
      end
end

#PLOT ERROR VS ITER: to compare to compare between minimizers, keeping the same strategy (easily adjustable to compare between strategies)
error_iter = Plots.plot(1:length(error_res["11"]), error_res["11"], yaxis=:log10, title = string("Burger error vs iter"), ylabel = "Error", label = string(minimizers_short_name[1]), ylims = (0.0001,1))
plot!(error_iter, 1:length(error_res["12"]), error_res["12"], yaxis=:log10, label = string(minimizers_short_name[2]))
plot!(error_iter, 1:length(error_res["13"]), error_res["13"], yaxis=:log10, label = string(minimizers_short_name[3]))
plot!(error_iter, 1:length(error_res["14"]), error_res["14"], yaxis=:log10, label = string(minimizers_short_name[4]))
plot!(error_iter, 1:length(error_res["15"]), error_res["15"], yaxis=:log10, label = string(minimizers_short_name[5]))
plot!(error_iter, 1:length(error_res["16"]), error_res["16"], yaxis=:log10, label = string(minimizers_short_name[6]))
plot!(error_iter, 1:length(error_res["17"]), error_res["17"], yaxis=:log10, label = string(minimizers_short_name[7]))
plot!(error_iter, 1:length(error_res["18"]), error_res["18"], yaxis=:log10, label = string(minimizers_short_name[8]))

Plots.plot(error_iter)

#Use after having modified the analysis setting correctly --> Error vs iter: to compare different strategies, keeping the same minimizer
#error_iter = Plots.plot(1:length(error_res["11"]), error_res["11"], yaxis=:log10, title = string("Burger error vs iter"), ylabel = "Error", label = string(strategies_short_name[1]), ylims = (0.0001,1))
#plot!(error_iter, 1:length(error_res["21"]), error_res["21"], yaxis=:log10, label = string(strategies_short_name[2]))
#plot!(error_iter, 1:length(error_res["31"]), error_res["31"], yaxis=:log10, label = string(strategies_short_name[3]))
#plot!(error_iter, 1:length(error_res["41"]), error_res["41"], yaxis=:log10, label = string(strategies_short_name[4]))
#plot!(error_iter, 1:length(error_res["51"]), error_res["51"], yaxis=:log10, label = string(strategies_short_name[5]))
#plot!(error_iter, 1:length(error_res["61"]), error_res["61"], yaxis=:log10, label = string(strategies_short_name[6]))
#plot!(error_iter, 1:length(error_res["71"]), error_res["71"], yaxis=:log10, label = string(strategies_short_name[7]))

#PLOT ERROR VS TIME: to compare to compare between minimizers, keeping the same strategy
error_time = plot(times["11"], error_res["11"], yaxis=:log10, label = string(minimizers_short_name[1]),title = string("Burger error vs time"), ylabel = "Error", size = (1500,500))
plot!(error_time, times["12"], error_res["12"], yaxis=:log10, label = string(minimizers_short_name[2]))
plot!(error_time, times["13"], error_res["13"], yaxis=:log10, label = string(minimizers_short_name[3]))
plot!(error_time, times["14"], error_res["14"], yaxis=:log10, label = string(minimizers_short_name[4]))
plot!(error_time, times["15"], error_res["15"], yaxis=:log10, label = string(minimizers_short_name[5]))
plot!(error_time, times["16"], error_res["16"], yaxis=:log10, label = string(minimizers_short_name[6]))
plot!(error_time, times["17"], error_res["17"], yaxis=:log10, label = string(minimizers_short_name[7]))
plot!(error_time, times["18"], error_res["18"], yaxis=:log10, label = string(minimizers_short_name[7]))

Plots.plot(error_time)

#Use after having modified the analysis setting correctly --> Error vs time: to compare different strategies, keeping the same minimizer
#error_time = plot(times["11"], error_res["11"], yaxis=:log10, label = string(strategies_short_name[1]),title = string("Burger error vs time"), ylabel = "Error", size = (1500,500))
#plot!(error_time, times["21"], error_res["21"], yaxis=:log10, label = string(strategies_short_name[2]))
#plot!(error_time, times["31"], error_res["31"], yaxis=:log10, label = string(strategies_short_name[3]))
#plot!(error_time, times["41"], error_res["41"], yaxis=:log10, label = string(strategies_short_name[4]))
#plot!(error_time, times["51"], error_res["51"], yaxis=:log10, label = string(strategies_short_name[5]))
#plot!(error_time, times["61"], error_res["61"], yaxis=:log10, label = string(strategies_short_name[6]))
#plot!(error_time, times["71"], error_res["71"], yaxis=:log10, label = string(strategies_short_name[7]))

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
# ------------------------------------------------------------------------
# UTILS

# ------------------------------------------------------------------------
# DIRS

# root_dir
root_dir(gl::GitLink) = gl.root_dir

# repo_dir
const _REPO_DIR_NAME = ".gl-repo"
repo_dir(gl::GitLink) = joinpath(root_dir(gl), _REPO_DIR_NAME)

# state dirs
# glob_state_dir
const _GLOBAL_STATE_DIR_NAME = ".gl-glob-state"
global_state_dir(gl::GitLink) = joinpath(repo_dir(gl), _GLOBAL_STATE_DIR_NAME)

# local_state_dir
const _LOCAL_STATE_DIR_NAME = ".gl-local-state"
local_state_dir(gl::GitLink) = joinpath(root_dir(gl), _LOCAL_STATE_DIR_NAME)

# stage_dir
const _STAGE_DIR_NAME = ".gl-stage"
stage_dir(gl::GitLink) = joinpath(root_dir(gl), _STAGE_DIR_NAME)


text:
"""
`checkinput!(data)`: Check and/or initialize all fields in data dictionary.

where:  `data`::Dict{Symbol, Any}

"""
function checkinput!(data::Dict{Symbol, Any})
  #
  # Check and/or set all input datafields that do not depend on the unitsystem
  # for the checks on the provided values.
  #
  if !(:isotropic in keys(data))
    #
    # Might be updated when :newmaterials are added containing an :er
    # in addition to an :ez value. In that case the material is
    # anisotropic and :er is the 'radial' Young modulus for pipes
    # (used for both :ezx and :ezy in the stiffness matrix).
    #
    # The :isotropic setting will only be updated if an anisotropic material
    # is actually used in the definitions for the segments.
    #
    # Note that either a single :ez value needs to be specified or 2 values (:er and :ez).
    #
    push!(data, :isotropic => true)
  end

  if :ratio in keys(data)
    checkratio!(data)
  else
    push!(data, :ratio => 0.6)
  end

  if :maxiters in keys(data)
    checkmaxiters!(data)
  else
    push!(data, :maxiters => 1000)
  end

  if !(:dlerrormax in keys(data))
    push!(data, :dlerrormax => 1e-6)
  end

  if !(:drerrormax in keys(data))
    push!(data, :drerrormax => 1e-6)
  end

  if !(:maxiters in keys(data))
    push!(data, :maxiters => 1000)
  end

  if !(:media in keys(data))
    createmediatable!(data)
  end

  if !(:materials in keys(data))
    creatematerialstable!(data)
  end

  if :newmedia in keys(data)
    updatemediatable!(data)
  end

  if :newmaterials in keys(data)
    updatematerialstable!(data)
  end

  if !(:segments in keys(data))
    println("No segments defined! Exiting.")
    data[:terminated => true]
    return
  else
    checksegmenttable!(data)
  end

  if :medium in keys(data)
    checkmedium!(data)
  else
    println(":medium not defined. Set to :mud.")
    push!(data, :medium => :mud)
  end

  if !(:surveys in keys(data))
    if !(:welltrajectory in keys(data))
      println("Neither surveys or trajectory data defined!")
      data[:terminated => true]
      return
    end
  else
    checksurveytable!(data)
  end

  if !(:wob in keys(data))
    println("No wob defined! Set to 0.0")
    push!(data, :wob => 0.0)
  else
    data[:wob] *= 1000.0
  end

  if :use_geometric_correction in keys(data)
    checkgeometriccorrectionsettings!(data)
  else
    push!(data, :use_geometric_correction => false)
  end

  if :usetortuosity in keys(data)
    checktortuositysettings!(data)
  end

end

"""
# The checkratio! function
`checkratio!(data)`: Check and/or set the ratio field in data.

  where:
  `data`::Dict{Symbol, Any}

"""
function checkratio!(data::Dict{Symbol, Any})
  if data[:ratio] < 0.01 || data[:ratio] > 10.0
    data[:ratio] = 0.6
    println("Setting of :ratio updated to $(data[:ratio])")
  end
end

function checkmaxiters!(data::Dict{Symbol, Any})
  if data[:maxiters] < 1 || data[:maxiters] > Int(1e6)
    data[:maxiters] = 5000
    println("Setting of :maxiters updated to $(data[:maxiters])")
  end
end

function updatemediatable!(data::Dict{Symbol, Any})
  if :newmedia in keys(data)
    if typeof(data[:newmedia]) == Dict{Symbol, Float64}
      merge!(data[:media], data[:newmedia])
      delete!(data, :newmedia)
    end
  end
end

function updatematerialstable!(data::Dict{Symbol, Any})
  if :newmaterials in keys(data)
    if typeof(data[:newmaterials]) == Dict{Symbol,Dict{Symbol,Float64}}
      for key in keys(data[:newmaterials])
        if key in keys(data[:materials])
          merge!(data[:materials][key], data[:newmaterials][key])
        else
          if !(:sg in keys(data[:newmaterials][key]))
            println("No :sg specified for new material $(key)")
            data[:terminated => true]
            return
          end
          if !(:g in keys(data[:newmaterials][key]))
            println("No :g specified for new material $(key)")
            data[:terminated => true]
            return
          end
          if !(:ez in keys(data[:newmaterials][key]))
            println("No :ez specified for new material $(key)")
            data[:terminated => true]
            return
          end
          if !(:er in keys(data[:newmaterials][key]))
            println("No :er specified for new material $(key). Assumed isotropic material.")
            data[:newmaterials][key][:er] = data[:newmaterials][key][:ez]
          end
          push!(data[:materials], key => data[:newmaterials][key])
        end
      end
      delete!(data, :newmaterials)
    end
  end
end

function checksegmenttable!(data::Dict{Symbol, Any})
  println("Checking $(length(data[:segments])) segments")
  if !(:holedepth in keys(data))
    push!(data, :holedepth => 0.0)
  end
  push!(data, :noofelements => 0)
  push!(data, :noofsegments => length(data[:segments]))
  if typeof(data[:segments][1]) != BHAtp.Bit
    println("First segment needs to be a Bit segment.")
    data[:terminated => true]
    return
  end
  for i in 1:length(data[:segments])
    if typeof(data[:segments][i]) == BHAtp.Rig
      if i != data[:noofsegments]
        println("Only last segment can be a Rig segment.")
        data[:terminated => true]
        return
      end
    end
    checkforanisotropiccomponent!(data[:segments][i], data)
    if !(data[:segments][i].material in keys(data[:materials]))
      println("Material not defined. Set to :steel")
      data[:segments][i].material = :steel
    end
    if data[:segments][i].id < 0.0
      println("Element id not set properly (< 0.0) for segment $i")
      data[:terminated => true]
      return
    end
    if data[:segments][i].od < 0.0
      println("Element od not set properly (< 0.0) for segment $i")
      data[:terminated => true]
      return
    end
    if data[:segments][i].id >= data[:segments][i].od
      println("Element id >= od for segment $i")
      data[:terminated => true]
      return
    end
    if data[:segments][i].fc < 0.0
      println("Element fc not set properly (< 0.0) for segment $i")
      data[:terminated => true]
      return
    end
    if typeof(data[:segments][i])  in [Pipe, Collar]
      if data[:segments][i].length < 0.0
        println("Length of segment $i < 0.0")
        data[:terminated => true]
        return
      end
      checksegmentrecord!(data[:segments][i], data)
      data[:noofelements] +=
        floor(Int64, 1.0+data[:ratio]*data[:segments][i].length/data[:segments][i].od)
    end
  end
  data[:noofnodes] = data[:noofelements] + 1
end

function checksurveytable!(data::Dict{Symbol, Any})
  println("Checking $(length(data[:surveys])) survey records")
  if !(:surveylength in keys(data))
    if size(data[:surveys], 1) > 0
      push!(data, :surveylength => data[:surveys][size(data[:surveys], 1)].depth)
    else
      push!(data, :surveylength => 0.0)
    end
  end
  for i in 1:length(data[:surveys])
    if !(0.0 <= data[:surveys][i].inclination <= 360.0)
      println("Survey inclination not set properly for survey record $i")
      data[:terminated => true]
      return
    end
    if !(0.0 <= data[:surveys][i].heading <= 360.0)
      println("Survey heading not set properly for survey record $i")
      data[:terminated => true]
      return
    end
    if data[:surveys][i].diameter < 0.0
      println("Survey radius not set properly (< 0.0) for survey record $i")
      data[:terminated => true]
      return
    end
    if typeof(data[:surveys][i]) == HoleSurvey
      if data[:surveys][i].depth < 0.0
        println("Depth of survey record $i < 0.0")
        data[:terminated => true]
        return
      end
      checksurveyrecord!(data[:surveys][i], data)
    end
  end
end

"""
`checktrajectoryinput!(data)`: Check trajectory data in data dictionary.

where:  `data`::Dict{Symbol, Any}

"""
function checktrajectoryinput!(data::Dict{Symbol, Any})
  if :welltrajectory in keys(data)
    welltraj = data[:welltrajectory]
    if !(:inclination in keys(welltraj))
      println("No inclination specified in well trajectory.")
      push!(data, :terminated => true)
      return
    end
    if !(:heading in keys(welltraj))
      println("No heading specified in well trajectory.")
      push!(data, :terminated => true)
      return
    end
    if !(:diameter in keys(welltraj))
      println("No hole diameter specified in well trajectory.")
      push!(data, :terminated => true)
      return
    end
  else
    println("No well trajectory data in input dictionary.")
    push!(data, :terminated => true)
    return
  end
end

function checkforcestable!(data::Dict{Symbol, Any})
  println("Checking $(length(data[:forces])) external force records")
  if !(:surveylength in keys(data))
    if size(data[:surveys], 1) > 0
      push!(data, :surveylength => data[:surveys][size(data[:surveys], 1)].depth)
    else
      push!(data, :surveylength => 0.0)
    end
  end
  for i in 1:length(data[:surveys])
    if !(0.0 <= data[:surveys][i].inclination <= 360.0)
      println("Survey inclination not set properly for survey record $i")
      data[:terminated => true]
      return
    end
    if !(0.0 <= data[:surveys][i].heading <= 360.0)
      println("Survey heading not set properly for survey record $i")
      data[:terminated => true]
      return
    end
    if data[:surveys][i].radius < 0.0
      println("Survey radius not set properly (< 0.0) for survey record $i")
      data[:terminated => true]
      return
    end
    if typeof(data[:surveys][i]) == HoleSurvey
      if data[:surveys][i].depth < 0.0
        println("Depth of survey record $i < 0.0")
        data[:terminated => true]
        return
      end
      #checksurveyrecord!(data[:surveys][i], data)
    end
  end
end

function checksegmentrecord!(rec::Collar, data::Dict{Symbol, Any})
  data[:holedepth] += rec.length
end

function checksegmentrecord!(rec::Pipe, data::Dict{Symbol, Any})
  data[:holedepth] += rec.length
end

function checksurveyrecord!(rec::HoleSurvey, data::Dict{Symbol, Any})
  #data[:surveylength] += rec.depth
end

function checkmedium!(data::Dict{Symbol, Any})
  if !(data[:medium] in keys(data[:media]))
    println("Unknown :medium specified. Set to :mud.")
    data[:medium] = :mud
  end
end

function checkforanisotropiccomponent!(seg::SegmentType, data::Dict{Symbol, Any})
  if :er in keys(data[:materials][seg.material])
    data[:isotropic] = false
  end
end

function checkgeometriccorrectionsettings!(data::Dict{Symbol, Any})
  if data[:use_geometric_correction]
    println("Geometric correction is being used.")
  end
end

function checktortuositysettings!(data::Dict{Symbol, Any})
  if (:tortuosity in keys(data)) && (:tortuositylength in keys(data))
    println("Tortuosity is injected.")
  else
    data[:usetortuosity] = false
    println("Artificial tortuosity was requested but not properly defined.")
  end
end


text:
"""
Creates a starting population of size `m` using a uniform random distribution,
where `a` is the lower bound, and `b` is the upper bound for the design
variables.
"""
function rand_population_uniform(m, a, b)
    d = length(a)
    return [a+rand(d).*(b-a) for i in 1:m]
end

abstract type SelectionMethod end
struct TruncationSelection <: SelectionMethod
    k # top k to keep
end
function select(t::TruncationSelection, y)
    p = sortperm(y)
    return [p[rand(1:t.k, 2)] for i in y]
end

abstract type CrossoverMethod end
struct SinglePointCrossover <: CrossoverMethod end
function crossover(::SinglePointCrossover, a, b)
    i = rand(1:length(a))
    return vcat(a[1:i], b[i+1:end])
end

abstract type MutationMethod end
struct GaussianMutation <: MutationMethod
    σ
end
function mutate(M::GaussianMutation, child)
    return child + randn(length(child))*M.σ
end

"""
Performs a Genetic Algorithm optimization with optional parallelization.

  Args:
    f (Function): Objective function.
    population (Array[]): An Array of design variable combinations
      representing an initial population.
    k_max (Integer): Maximum iterations
    S (SelectionMethod): The selection method to be used.
    C (CrossoverMethod): The crossover method to be used.
    M (MutationMethod): The mutation method to be used.
    parallel (Boolean, optional): Optionally use multithreading.

  Returns:
    Array: An Array representing the minimizers for the given objective
      function.
"""
function genetic_algorithm(f, population, k_max, S, C, M; parallel=false)
    n = length(population)
    for k in 1 : k_max
        if parallel
            # Perform function evaluations in parallel
            f_pop = zeros(n)
            Threads.@threads for i=1:n
               f_pop[i] = f(population[i])
            end
        else
            f_pop = f.(population)
        end
        # proceed as normal for selection, crossover, mutation
        parents = select(S, f_pop)
        children = [crossover(C, population[p[1]], population[p[2]]) for p in parents]
        population .= mutate.(Ref(M), children)
    end

    if parallel
        # Perform function evaluations in parallel
        f_pop = zeros(n)
        Threads.@threads for i=1:n
            f_pop[i] = f(population[i])
        end
    else
        f_pop = f.(population)
    end
    population[argmin(f_pop)]
end

max_stars_repo_name: camirmas/distributed-ga

max_issues_repo_name: camirmas/distributed-ga

max_forks_repo_name: camirmas/distributed-ga


text:
# for over/underflow in logs
"""
	softplus_safe(x,T=Float32)

A softplus with small additive constant for safe operations.
"""
softplus_safe(x,T=Float32) = softplus(x) .+ T(1e-6)

### Flux chain builders ###
"""
    function layer_builder(d,k,o,n,ftype,lastlayer = "",ltype = "Dense")

Create a chain with `n` layers of with `k` neurons with activation function `ftype`.
Input and output dimension is `d` / `o`.
If lastlayer is not specified, all layers use the same function.
If lastlayer is "linear", then the last layer is forced to be Dense with linear activation.

It is also possible to specify dimensions in a vector.
"""
layer_builder(k::Vector{Int},l::Vector,f::Vector) =
    Flux.Chain(map(i -> i[1](i[3],i[4],i[2]),zip(l,f,k[1:end-1],k[2:end]))...)

layer_builder(d::Int,k::Int,o::Int,n::Int, args...) =
    layer_builder(vcat(d,fill(k,n-1)...,o), args...)

function layer_builder(ks::Vector{Int},ftype::String,lastlayer::String = "",ltype::String = "Dense")
    ftype = (ftype == "linear") ? "identity" : ftype
    ls = Array{Any}(fill(eval(:($(Symbol(ltype)))),length(ks)-1))
    fs = Array{Any}(fill(eval(:($(Symbol(ftype)))),length(ks)-1))
    if !isempty(lastlayer)
        fs[end] = (lastlayer == "linear") ? identity : eval(:($(Symbol(lastlayer))))
        ls[end] = (lastlayer == "linear") ? Dense : ls[end]
    end
    layer_builder(ks,ls,fs)
end

"""
    stack_layers(lsize, activation, [layer=Dense; last=identity])

Construct a Flux chain (e.g. encoder/decoder) consisting of `length(lsize)-1` layers of
type `layer` with  `activation` in between. Default last activation is `identity`.
Width of layers is defined by the vector `lsize`.
"""
stack_layers(lsize::Union{Vector, Tuple}, activation, layer=Dense; last=identity) =
    layer_builder([x for x in lsize],
        Array{Any}(fill(layer, size(lsize,1)-1)),
        Array{Any}([fill(activation, size(lsize,1)-2); last]))

### training ###
"""
    update_params!(model, data, loss, opt)

Basic training step - computation of the loss, backpropagation of gradients and optimisation
of weights. The loss and opt arguments can be arrays/lists/tuples.
"""
function update_params!(model, data, loss, opt)
    ps = Flux.params(model)
    gs = gradient(ps) do
        loss(data...)
    end
    Flux.Optimise.update!(opt, ps, gs)
end

"""
    train!(model, data, loss, optimiser, callback; [usegpu, memory_efficient])

Train the model. Function callback(model, batch, loss, opt) is
called every iteration - use it to store or print training progress, stop training etc.
"""
function train!(model, data, loss, optimiser, callback;
    usegpu = false, memory_efficient = false)
    for _data in data
        try
            if usegpu
             _data = _data |> gpu
            end
            update_params!(model, _data, loss, optimiser)
            # now call the callback function
            # can be an object so it can store some values between individual calls
            callback(model, _data, loss, optimiser)
        catch e
            # setup a special kind of exception for known cases with a break
            rethrow(e)
        end
        if memory_efficient
            # this might be important for conv nets running on gpu
            # large nets might still train but slowly
            GC.gc();
        end
    end
end

### GAN ###
"""
    gen_loss(sg::AbstractArray)

Generator loss for generated data score `sg`.
"""
function generator_loss(sg::AbstractArray{T}) where T<:Real
    - mean(log.(sg) .+ eps(T))
end

"""
    disc_loss(st::AbstractArray, sg::AbstractArray)

Discriminator loss for true data score `st` and generated data score `sg`.
"""
function discriminator_loss(st::AbstractArray{T}, sg::AbstractArray{T}) where T<:Real
    - T(0.5)*(mean(log.(st .+ eps(T))) + mean(log.(1 .- sg) .+ eps(T)))
end


text:
######################################################################
# fov.jl
#
# field of view sensor
######################################################################

struct FOVb <: Sensor
    cone_width::Float64
    alpha::Float64          # mistake rate

    function FOVb(cone_width::Real, alpha::Float64)
        return new(float(cone_width), alpha)
    end
    FOVb() = FOVb(120.0, 0.1)
end

# 1 means it is in the field of view of front antenna
# 0 means it is not
function observe(m::SearchDomain, s::FOVb, p::Pose)

    # ensure bearing is reflected across
    rel_bearing = fit_180(p[3] - true_bearing(p, m.theta))
    if rel_bearing < 0.0
        rel_bearing = -1.0 * rel_bearing
    end

    a1 = s.cone_width / 2.0
    a2 = 180.0 - a1

    prob_in_view = 0.0
    if rel_bearing < a1
        prob_in_view = 1.0 - s.alpha
    elseif rel_bearing < a2
        prob_in_view = 0.5
    else
        prob_in_view = s.alpha
    end

    o = (rand() < prob_in_view) ? true : false
    return o
end

# o technically has to be a float, but we only have two obs:
#  1 = in field of view
#  0 = not in field of view
function O(s::FOVb, theta::LocTuple, p::Pose, o)

    # determine relative bearing and fix it in 0 to 180
    rel_bearing = fit_180(p[3] - true_bearing(p, theta))
    if rel_bearing < 0.0
        rel_bearing = -1.0 * rel_bearing
    end

    a1 = s.cone_width / 2.0
    a2 = 180.0 - a1

    prob_in_view = 0.0  # probability in view of front antenna

    if rel_bearing < a1
        prob_in_view = 1.0 - s.alpha
    elseif rel_bearing < a2
        prob_in_view = 0.5
    else
        prob_in_view = s.alpha
    end

    # prob that it's one (in view of front antenna)
    ret_val = (o == true) ? prob_in_view : (1.0 - prob_in_view)

    return ret_val
end

get_prob(s::FOVb, theta::LocTuple, xp::Pose, o) = O(s, theta, xp, o)
export get_prob


text:
# examples.jl

#Importing Modules
import Quads.PokerCard
import Quads.PokerDeck
import Quads.Lookup
import Quads.Evaluator

# creating a deck
deck = PokerDeck.make_deck()

# there are 52 cards
@assert length(deck) == 52
# cards are integers
@assert typeof(deck[1]) == Int
#draw some cards
@assert length(PokerDeck.draw!(5, PokerDeck.make_deck())) == 5

# working wth cards
Ah = PokerCard.newCard("Ah") #268447017
@assert Ah == 268447017
@assert typeof(Ah) == Int
# there is a bitstring rep for each card that makes sense, check out PokerCard
@assert bitstring(Ah % UInt32) == "00010000000000000010110100101001"

# working with a hand, eg evaluating a Straight Flush
hand_string = ["2h", "3h", "4h", "5h", "Ah"]
@assert length(PokerCard.newCard.(hand_string)) == 5
hand_ints = PokerCard.newCard.(hand_string)
@assert Evaluator.class_to_string(Evaluator.get_rank_class(Evaluator.evaluate(hand_ints))) == "Straight Flush"

# seven cards from the deck and eval
hand_ints = PokerDeck.draw!(7, PokerDeck.make_deck())
@assert typeof(PokerCard.int_to_pretty_str.(hand_ints)) == Array{String,1}
@assert typeof(Evaluator.get_rank_class(Evaluator.evaluate(hand_ints))) == Int

max_stars_repo_name: HigginTown/Quads

max_issues_repo_name: HigginTown/Quads

max_forks_repo_name: HigginTown/Quads


text:
"""
   inverse

Inverse a permutation.
If n > 0, fill with zeros. Can also use `invperm`
"""
function inverse(p::Vector{Int}, n::Int = 0)
    if n <= 0
        n = length(p)
    end
    q = zeros(Int, n)
    for (i, j) in enumerate(p)
        q[j] = i
    end
    q
end

"""
    concat(a, b)

TODO
"""
concat(a::Vector{T}, i::Vector{Int}) where T =
    T[a[j] for j in i if j > 0]

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

add_source = load_op_and_grad("../build/libADSeismic","add_source", multiple=true)
################## End Load Operator ##################
# tf.debugging.set_log_device_placement(true)
# TODO: specify your input parameters
nx = 20
ny = 20
sigmaxx = constant(ones((nx+2)*(ny+2)))
sigmayy = constant(ones((nx+2)*(ny+2)))
sigmaxy = constant(ones((nx+2)*(ny+2)))
vx = constant(ones((nx+2)*(ny+2)))
vy = constant(ones((nx+2)*(ny+2)))
srci = constant(rand(1:5, 100))
srcj = constant(rand(1:5, 100))
srctype = constant(rand(0:4,100))
srcv = constant(rand(100))
u = add_source(sigmaxx,sigmayy,sigmaxy,vx,vy,srci,srcj,srctype,constant(nx),constant(ny),srcv)
# error()
sess = Session()
init(sess)
run(sess, u)
# error()

# TODO: change your test parameter to `m`
# gradient check -- v
function scalar_function(m)
    c = sum(m^2)
    return sum(sum.(add_source(sigmaxx + c,sigmaxy + c,sigmayy + c,vx + c,vy + c,srci,srcj,srctype,constant(nx),constant(ny),m).^2))
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(100))
v_ = rand(100)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")
savefig("gradtest.png")


text:
function createLogDatabase(dbpath::String)

    db = SQLite.DB(dbpath)

    DBInterface.execute(db, """
        CREATE TABLE IF NOT EXISTS indexes(
            id TEXT PRIMARY KEY,
            timestamp INTEGER,
            numrecords INTEGER
        )""")

    DBInterface.execute(db, """
        CREATE TABLE IF NOT EXISTS records(
            id TEXT PRIMARY KEY,
            timestamp INTEGER,
            content BLOB NOT NULL
        )""")

    return (
        DBInterface.prepare(db, "INSERT INTO indexes VALUES(?,?,?);"),
        DBInterface.prepare(db, "INSERT INTO records VALUES(?,?,?);"),
        DBInterface.prepare(db, "SELECT id FROM indexes;"),
        DBInterface.prepare(db, "SELECT id FROM records;"),
        DBInterface.prepare(db, "BEGIN TRANSACTION;"),
        DBInterface.prepare(db, "COMMIT;")
    )
end

function getLogContent(m::RegexMatch)

    logid::AbstractString = m[:id]

    datetime = Dates.datetime2epochms(DateTime(
            parse(Int, logid[1:4]),
            parse(Int, logid[5:6]),
            parse(Int, logid[7:8]),
            parse(Int, logid[9:10]),
            parse(Int, m[:min])
        ))

    buffer = IOBuffer()
    retry(Downloads.download; delays = [0.2, 0.5, 1, 2, 3])(
            "http://tenhou.net/0/log/?" * logid, buffer
        )

    content = transcode(LZ4FrameCompressor, take!(buffer))

    return logid, datetime, content
end

function buildLogDatabase(indexpath::String, dbpath::String, gametype::String)

    @assert startswith(indexpath, "scraw") && endswith(indexpath, "zip")

    # open database and compile sql statements
    insertind, insertrec, selectind, selectrec, beginsql, endsql =
            createLogDatabase(dbpath)
    println("* Database connected - ", dbpath)

    # get list of subindexes already appended to database
    appended::Vector{String} = Tables.columntable(DBInterface.execute(selectind)).id

    # create regular expression for specified gametype
    rx = Regex(":(?<min>\\d\\d)[\\W\\d]+?$gametype[^?]+\\?log=(?<id>[^\"]+)")

    # iterate over all unique html subindexes in index archive
    zipfile = ZipFile.Dir(indexpath)
    for file in zipfile.files
        if (occursin("scc", file.name) && !(file.name in appended))

            # read subindex (unzip if necessary)
            subindex = read(file, String)
            if occursin("gz", file.name)
                subindex = String(transcode(GzipDecompressor, subindex))
            end

            # download and append log content to database
            matches = collect(eachmatch(rx, subindex))
            DBInterface.execute(beginsql)
            @sync for match in matches
                @async DBInterface.execute(insertrec, getLogContent(match))
            end

            # append subindex information
            DBInterface.execute(insertind, (
                file.name,
                Dates.datetime2epochms(now(UTC)),
                length(matches)
            ))

            # write available data into database
            DBInterface.execute(endsql)
            println(file.name, ": ", length(matches), " logs inserted")
        end
    end

    # after all work is done
    close(zipfile)
    println("* Database complete - ", dbpath)
end

function downloadProgress()

    local lastprinted = Ref(0)
    (total, current) -> begin

        if (current > lastprinted[] + 200000)   # print every 200kb
            lastprinted[] = current
            write(stdout, "\e[1G" * "< Downloading: $current/$total\t")
        end

        return nothing
    end
end

function downloadLogIndex(year::String)

    @assert 2009 <= parse(Int, year) <= Dates.year(now(UTC))

    result = Downloads.download("http://tenhou.net/sc/raw/scraw$year.zip",
            "scraw$year.zip", progress = downloadProgress(), verbose = true)
    write(stdout, "* Download complete: scraw$year.zip\n\n")

    return result
end

const S4P_GAME = "四.南"
const S3P_GAME = "三.南"
const E4P_GAME = "四.東"
const E3P_GAME = "三.東"


text:
using PGFPlotsX, SpectralKit, Colors, ForwardDiff

"""
Plot selected basis functions from a basis.

# Arguments

- `basis`: the basis we use
- `ns`: the polynomials which are picked
- `order`: 0 (value) or 1 (derivative)
- `title`: plot title

# Keyword arguments

- `xmin`, `xmax`: boundaries of the grid for calculating points
- `len`: number of points
- `rel_offsets`: extends the plotting area by a fraction (one value for each side)
"""
function plot_basis(basis, ns, order, title;
                    xmin = domain(basis)[1], xmax = domain(basis)[2],
                    rel_offsets = [0.01, 0.01], len = 500)
    @assert order ∈ 0:1
    axis = @pgf Axis({ xmin = xmin, xmax = xmax, xlabel = "y", title = title,
                       legend_pos = "outer north east"})
    xgrid = range(xmin, xmax; length = len)
    mi, ma = domain(basis)
    colors = distinguishable_colors(length(ns); lchoices = 0:60)
    offsets = [1, -1] .* rel_offsets .* (xmax - xmin)
    ev0(n, x) = first(Iterators.drop(basis_at(basis, x), n - 1)) # NOTE horribly inefficient …
    ev1(n, x) = ForwardDiff.derivative(x -> ev0(n, x), x) # NOTE … but does not matter
    ev(n, x) = order == 0 ? ev0(n, x) : ev1(n, x)
    for (n, col) in zip(ns, colors)
        y = map(x -> ev(n, x), xgrid)
        @pgf push!(axis,
                   Plot({ no_marks, thick, color = col }, Table(xgrid, y)),
                   LegendEntry(string(n)))
    end
    for (i, (n, col)) in enumerate(zip(ns, colors))
        y = [ev(n, mi), ev(n, ma)]
        @pgf push!(axis, Plot({ only_marks, color = col },
                              Table([xmin, xmax] .+ offsets .* i, y)))
    end
    axis
end

N = 5
ns = 1:N                        # polynomials we generate

###
### Chebyshev
###

p = plot_basis(Chebyshev(InteriorGrid(), N), ns, 0, "Chebyshev polynomials (value)")
pgfsave("chebyshev.png", p; dpi = 600)

p = plot_basis(Chebyshev(InteriorGrid(), N), ns, 1, "Chebyshev polynomials (derivative)";
                rel_offsets = [0, 0])
pgfsave("chebyshev_deriv.png", p; dpi = 600)

###
### Rational Chebysev [0,∞)
###

basis = univariate_basis(Chebyshev, InteriorGrid(), N, SemiInfRational(0.0, 1.0))
p = plot_basis(basis, ns, 0, raw"Rational Chebyshev functions $[0,\infty)$ (value)";
               xmax = 1)
pgfsave("semiinf.png", p; dpi = 600)

p = plot_basis(basis, ns, 0, raw"Rational Chebyshev functions $[0,\infty)$ (value)";
               xmax = 200)
pgfsave("semiinf_birdseye.png", p; dpi = 600)

p = plot_basis(basis, ns, 1, raw"Rational Chebyshev functions $[0,\infty)$ (derivative)";
               xmax = 1, rel_offsets = [0, 0.01])
pgfsave("semiinf_deriv.png", p; dpi = 600)

###
### Rational Chebysev (-∞,∞)
###

basis = univariate_basis(Chebyshev, InteriorGrid(), N, InfRational(0.0, 1.0))
p = plot_basis(basis, ns, 0,
               raw"Rational Chebyshev functions $(-\infty,\infty)$ (value)";
               xmin = -1, xmax = 1)
pgfsave("inf.png", p; dpi = 600)

p = plot_basis(basis, ns, 0,
               raw"Rational Chebyshev functions $(-\infty,\infty)$ (value)";
               xmin = -30, xmax = 30, len = 1000)
pgfsave("inf_birdseye.png", p; dpi = 600)

p = plot_basis(basis, ns, 1,
               raw"Rational Chebyshev functions $(-\infty,\infty)$ (derivative)";
               xmin = -5, xmax = 5, rel_offsets = [0.01, 0.01])
pgfsave("inf_deriv.png", p; dpi = 600)


text:
function set_font_matrix(ctx, matrix)
    ccall((:cairo_set_font_matrix, Cairo.libcairo), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), ctx.ptr, Ref(matrix))
end

function get_font_matrix(ctx)
    matrix = Cairo.CairoMatrix()
    ccall((:cairo_get_font_matrix, Cairo.libcairo), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), ctx.ptr, Ref(matrix))
    return matrix
end

function cairo_font_face_destroy(font_face)
    ccall(
        (:cairo_font_face_destroy, Cairo.libcairo),
        Cvoid, (Ptr{Cvoid},),
        font_face
    )
end

function set_ft_font(ctx, font)

    font_face = ccall(
        (:cairo_ft_font_face_create_for_ft_face, Cairo.libcairo),
        Ptr{Cvoid}, (Makie.FreeTypeAbstraction.FT_Face, Cint),
        font, 0
    )

    ccall((:cairo_set_font_face, Cairo.libcairo), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), ctx.ptr, font_face)

    return font_face
end

"""
Finds a font that can represent the unicode character!
Returns Makie.defaultfont() if not representable!
"""
function best_font(c::Char, font = Makie.defaultfont())
    if Makie.FreeType.FT_Get_Char_Index(font, c) == 0
        for afont in Makie.alternativefonts()
            if Makie.FreeType.FT_Get_Char_Index(afont, c) != 0
                return afont
            end
        end
        return Makie.defaultfont()
    end
    return font
end