text:
## Ported Joel CFD's example to Julia with extensions: https://github.com/joel-ho/Field-inversion-example
# Currently uses Zygote for autodiff
using LinearAlgebra
using Zygote
using StaticArrays
using SparseArrays
using Base.Iterators

# Newton iteration
prod_size(tensor) = (prod ∘ size)(tensor)

newton_step(u, R, ∂R∂u) = reshape(∂R∂u \ -R[:], size(u))

## Finite differencing
central_2nd(Tm, T, Tp) = Tm - 2T + Tp
central_2nd_diff(Tm, T, Tp, h) = central_2nd(Tm, T, Tp) / h^2
central_1st_diff(Tm, Tp, h) = (Tm - Tp) / 2h

compute_q(x, y) = 5e5 * (1 - (x - 0.5)^2) * (1 - (y - 0.5)^2)

function get_T(T_grid, i, j, T_boundary)
    if i < 1 || j < 1 || i > length(T_grid[:,1])|| j > length(T_grid[1,:])
        T_boundary
    else
        T_grid[i,j]
    end
end

stencil_1(T, i, j, T_boundary) =
    [ get_T(T, i, j, T_boundary)   ;
      get_T(T, i-1, j, T_boundary) ;
      get_T(T, i+1, j, T_boundary) ;
      get_T(T, i, j-1, T_boundary) ;
      get_T(T, i, j+1, T_boundary) ]

# Difference operators specific to stencil_1
Δ²x(T, h) = central_2nd_diff(T[2], T[1], T[3], h)
Δ²y(T, k) = central_2nd_diff(T[4], T[1], T[5], k)
Δx(T, h) = central_1st_diff(T[3], T[2], h)
Δy(T, k) = central_1st_diff(T[5], T[4], k)

compute_k(β₁, β₂, β₃, β₄,  T) = β₁/T + β₂ + β₃ * T + β₄ * T^2
# + β₃ * T + β₄ * T^2
∇compute_k(β₁, β₂, β₃, β₄, T) = gradient(compute_k, β₁, β₂, β₃, β₄, T)

function residual_cell(T, β, q, dx, dy)
    k = compute_k(β..., first(T))
    ∂k∂T = (last ∘ ∇compute_k)(β..., first(T))

    k * (Δ²x(T, dx) + Δ²y(T, dy)) + ∂k∂T * (Δx(T, dx)^2 + Δy(T, dy)^2 ) + q
end

residual(T, β, q, dx, dy, T_boundary) =
        [ let Ts = stencil_1(T, T_ind.I..., T_boundary);
            residual_cell(Ts, β, q[T_ind], dx, dy) end
            for T_ind in CartesianIndices(T) ]

# Derivative stuff
function solve_direct(x, u, ∂R∂x, ∂R∂u)
    ∂R∂u_sq = reshape(∂R∂u, (length(u[:]), length(u[:])))
    reshape(reduce(hcat, ∂R∂u_sq \ -(∂R∂x)[:,:,i][:] for i in eachindex(x)), (size(u)..., length(x)))
end

function solve_adjoint(u, ∂R∂u, dfdu)
    reshape(∂R∂u, (length(u[:]), length(u[:])))' \ -(dfdu)'[:]
end

total_derivative_direct(∂f∂x, ψ, ∂f∂u) = ∂f∂x + [ sum(∂f∂u * ψ[n]) for n in eachindex(∂f∂x) ]

total_derivative_adjoint(∂f∂x, φ, ∂R∂x) = ∂f∂x + [ sum(permutedims(φ) * reshape(∂R∂x, (length(R[:]), length(∂f∂x)))[:,n]) for n in eachindex(∂f∂x) ]

∇residual_cell(Ts, βs, q, dx, dy) = gradient(residual_cell, Ts, βs, q, dx, dy)

function ∇residual_cell_T_update!(∂R∂T, grad, k, l)
    a, b, c, d = size(∂R∂T)
    ∂R∂T[k,l,k,l]            = grad[1]
    if k > 1 ∂R∂T[k,l,k-1,l] = grad[2] end
    if k < a ∂R∂T[k,l,k+1,l] = grad[3] end
    if l > 1 ∂R∂T[k,l,k,l-1] = grad[4] end
    if l < b ∂R∂T[k,l,k,l+1] = grad[5] end
end

∇residual(T, β, q, dx, dy, T_boundary) = gradient(residual, T, β, q, dx, dy, T_boundary)

∇residual_β(T, β, q₀, dx, dy, T_boundary) =
        [ let Ts = stencil_1(T, T_ind.I..., T_boundary);
            ∇residual_cell(Ts, β, q₀[T_ind], dx, dy)[2][i] end for
            T_ind in CartesianIndices(T), i in eachindex(β) ]

function compute_residual_and_grad(T, β, q, dx, dy)
    # Residual
    R = residual(T, β, q, dx, dy, T_boundary)

    # Jacobian
    T_size = size(T)
    ∂R∂T = zeros(T_size..., T_size...)
    for T_ind in CartesianIndices(T)
        Ts = stencil_1(T, T_ind.I..., T_boundary)
        grad = ∇residual_cell(Ts, β, q₀[T_ind], dx, dy)[1]
        ∇residual_cell_T_update!(∂R∂T, grad, T_ind.I...)
    end

    R, sparse(reshape(∂R∂T, size(R) .* reverse(T_size)))
end

## Grid setup
dx = 0.1
dy = 0.2
xs, ys = -1:dx:1, -1:dy:1
grid = [ (x, y) for x in xs, y in ys ]
T_boundary = 200.
βs = [13221., 27., 5., 1.]
ω = 1.0

T₀ = T_boundary * ones(size(grid))
q₀ = map(xs -> compute_q(xs...), grid);

n = length(xs)
T_truth = deepcopy(T₀)

num_iters = 30
ε = zeros(num_iters);

## NEWTON ITERATIONS
#=================================================================#

for i in 1:num_iters
    R, ∂R∂T = compute_residual_and_grad(T_truth, βs, q₀, dx, dy)
    ΔT = newton_step(T_truth, R, ∂R∂T)
    ε[i] = maximum(abs.(ΔT))
    # display(R)
    println("Newton step error L² norm: $(ε[i])")
    T_truth .+= ω * ΔT
end

##
k_truth = compute_k.(βs..., T_truth)

##
using Plots
gr()

##
plot(ε)

##
surf1 = surface(xs, ys, q₀)

##
surf2 = surface(xs, ys, T_truth)

##
surf3 = surface(xs, ys, k_truth)

## ADJOINT ITERATIONS
#=================================================================#

function cost(T, T_truth, indices)
    norm(T[index] - T_truth[index] for index in indices)
end

∇cost(T, T_truth, indices) = gradient(cost, T, T_truth, indices)

# indices = rand(1:length(T_truth[:]), (Int ∘ floor)(length(T_truth[:]) / 2))
indices = eachindex(T_truth)

##
β =  [13000., 10., 10., 5.]
γs = [100., 5., 5., 5.]

num_adjoint = 100
num_iters   = 30

T₀   = T_boundary * ones(size(grid))
q₀   = map(xs -> compute_q(xs...), grid);

T    = copy(T₀)
ε    = zeros(num_adjoint, num_iters)
εa   = zeros(num_adjoint)
β_l  = [β'; zeros(num_adjoint - 1, length(β)) ];
R    = zeros(size(T))
∂R∂T = zeros(size(R)..., size(T)...);

##
for i in 1:num_adjoint

    println("\nAdjoint Iteration — $i")
    display(β)

    for j in 1:num_iters
        R, ∂R∂T = compute_residual_and_grad(T, β, q₀, dx, dy)
        ΔT      = newton_step(T, R, ∂R∂T)
        ε[i,j]  = maximum(abs.(ΔT))
        # println("Newton step error: $(ε[i,j])")
        T     .+= ω * ΔT
    end

    εa[i] = cost(T, T_truth, indices)
    println("\nCost function value: $(εa[i])")

    # println("Computing df/dT")
    dfdT = ∇cost(T, T_truth, indices)[1]

    # println("Computing ∂f/∂x")
    ∂f∂x = zeros(length(β))

    # println("Computing ∂R/∂x")
    ∂R∂x = ∇residual_β(T, β, q₀, dx, dy, T_boundary)

    # display(∂R∂x)

    # println("Solving adjoint problem: (∂R/∂T)ᵗ φ = -(df/dT)ᵗ")
    # φ = solve_adjoint(T, ∂R∂T, dfdT)
    # dfdβ = total_derivative_adjoint(∂f∂x, φ, ∂R∂x)

    dTdx = solve_direct(β, T, ∂R∂x, ∂R∂T)
    dfdβ = total_derivative_direct(∂f∂x, dTdx, dfdT)

    # display(dfdβ)

    # println("Updating design variables β")
    @. β -= γs * dfdβ
    β_l[i,:] = β'
end

##
display(β)

##
k = compute_k.(β..., T)

##
using Plots
plotlyjs()

##
plot(εa, xlabel = "Iterations", ylabel = "Grid Error")

##
surface(xs, xs, T)
surface!(xs, xs, T_truth)

##
surface(xs, xs, k)
surface!(xs, xs, k_truth)

##
# residual_discrete_k(T, k, q, h) = first(k) * (Δ²x(T, h) + Δ²y(T, h)) + Δx(k, h) * Δx(T, h) + Δy(k, h) * Δy(T, h) + q
# ∇residual_discrete_k(Ts, βs, q, h) = gradient(residual_discrete_k, Ts, βs, q, h)

# function update_k_grid_ghost_points!(k_grid)
#     k_grid[1,:] = 2 * k_grid[2,:] .- k_grid[3,:]
#     k_grid[end,:] = 2 * k_grid[end-1,:] .- k_grid[end-2,:]
#     k_grid[:,1] = 2 * k_grid[:,2] .- k_grid[:,3]
#     k_grid[:,end] = 2 * k_grid[:,end-1] .- k_grid[:,end-2]

#     k_grid
# end

max_stars_repo_name: GodotMisogi/AdjointCFD

max_issues_repo_name: GodotMisogi/AdjointCFD

max_forks_repo_name: GodotMisogi/AdjointCFD


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 4inch)

births = dataset("HistData", "Arbuthnot")

plot(births, x=Col.value(:Year), y=Col.value(:Males, :Females),
     color=Col.index(:Males, :Females), Geom.line)

ext: jl


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_tabs() do
        dcc_tab() do
            html_button(id="btn", children = "Update Input"),
            html_div(id = "output", children = ["Hello"])
        end,
        dcc_tab() do
            html_div("test")
        end
    end
end

callback!(app, Output("output","children"), Input("btn","n_clicks")) do n_clicks
    isnothing(n_clicks) && throw(PreventUpdate())
    return "Bye"
end

run_server(app)


text:
function to_real_field_univariate_vector_function(obj::FieldUnivariateVectorFunction, arg0::Field)
    return jcall(obj, "toRealFieldUnivariateVectorFunction", RealFieldUnivariateVectorFunction, (Field,), arg0)
end

function value(obj::FieldUnivariateVectorFunction, arg0::RealFieldElement)
    return jcall(obj, "value", Vector{RealFieldElement}, (RealFieldElement,), arg0)
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using IntervalTrees
using SparseArrays
using NearestNeighbors

#---------------------------------------------------------------------
#	2D containment test
#---------------------------------------------------------------------

"""

Half-line crossing test. Utility function for `pointInPolygonClassification` function.
Update the `count` depending of the actual crossing of the tile half-line.
"""
function crossingTest(new::Int, old::Int, count::T, status::Int)::Number where T <: Real
    if status == 0
        status = new
        count += 0.5
    else
        if status == old
        	count += 0.5
        else
        	count -= 0.5
        end
        status = 0
    end
end

"""
	setTile(box)(point)

Set the `tileCode` of the 2D bbox `[b1,b2,b3,b4]:=[ymax,ymin,xmax,xmin]:= x,x,y,y`
including the 2D `point` of `x,y` coordinates.
Depending on `point` position, `tileCode` ranges in `0:15`, and uses bit operators.
Used to set the plane tiling depending on position of the query point,
in order to subsequently test the tile codes of edges of a 2D polygon, and determine
if the query point is either internal, external, or on the boundary of the polygon.
Function to be parallelized ...

```julia
c1,c2 = tilecode(p1),tilecode(p2)
c_edge, c_un, c_int = c1 ⊻ c2, c1 | c2, c1 & c2
```

"""
function setTile(box)
	tiles = [[9,1,5],[8,0,4],[10,2,6]]
	b1,b2,b3,b4 = box
	function tileCode(point)
		x,y = point
		code = 0
		if y>b1 code=code|1 end
		if y<b2 code=code|2 end
		if x>b3 code=code|4 end
		if x<b4 code=code|8 end
		return code
	end
	return tileCode
end

"""
	pointInPolygonClassification(V,EV)(pnt)

Point in polygon classification.

# Example

```julia
result = []
classify = pointInPolygonClassification(V,EV)
```
"""
function pointInPolygonClassification(V,EV)
    function pointInPolygonClassification0(pnt)
        x,y = pnt
        xmin,xmax,ymin,ymax = x,x,y,y
        tilecode = setTile([ymax,ymin,xmax,xmin])
        count,status = 0,0

        for (k,edge) in enumerate(EV)
            p1,p2 = V[:,edge[1]],V[:,edge[2]]
            (x1,y1),(x2,y2) = p1,p2
            c1,c2 = tilecode(p1),tilecode(p2)
            c_edge, c_un, c_int = c1⊻c2, c1|c2, c1&c2

            if (c_edge == 0) & (c_un == 0) return "p_on"
            elseif (c_edge == 12) & (c_un == c_edge) return "p_on"
            elseif c_edge == 3
                if c_int == 0 return "p_on"
                elseif c_int == 4 count += 1 end
            elseif c_edge == 15
                x_int = ((y-y2)*(x1-x2)/(y1-y2))+x2
                if x_int > x count += 1
                elseif x_int == x return "p_on" end
            elseif (c_edge == 13) & ((c1==4) | (c2==4))
                    crossingTest(1,2,status,count)
            elseif (c_edge == 14) & ((c1==4) | (c2==4))
                    crossingTest(2,1,status,count)
            elseif c_edge == 7 count += 1
            elseif c_edge == 11 count = count
            elseif c_edge == 1
                if c_int == 0 return "p_on"
                elseif c_int == 4 crossingTest(1,2,status,count) end
            elseif c_edge == 2
                if c_int == 0 return "p_on"
                elseif c_int == 4 crossingTest(2,1,status,count) end
            elseif (c_edge == 4) & (c_un == c_edge) return "p_on"
            elseif (c_edge == 8) & (c_un == c_edge) return "p_on"
            elseif c_edge == 5
                if (c1==0) | (c2==0) return "p_on"
                else crossingTest(1,2,status,count) end
            elseif c_edge == 6
                if (c1==0) | (c2==0) return "p_on"
                else crossingTest(2,1,status,count) end
            elseif (c_edge == 9) & ((c1==0) | (c2==0)) return "p_on"
            elseif (c_edge == 10) & ((c1==0) | (c2==0)) return "p_on"
            end
        end
        if (round(count)%2)==1
        	return "p_in"
        else
        	return "p_out"
        end
    end
    return pointInPolygonClassification0
end

#---------------------------------------------------------------------
#	Refactoring pipeline
#---------------------------------------------------------------------

"""
	input_collection(data::Array)::Tuple

*Facet selection*. Construction of a ``(d-1)``-dimensional collection from a ``(d-1)``-
or ``d``-dimensional one. ``0-chain`` of `LAR` type are used as *input*.

*Output* is ``admissible input`` for algorithms of the *2D/3D arrangement* pipeline.

# Example 2D

An assembly of geometric objects is generated, and their assembly, including rotated
and translated chains, is built producing a collection of input LAR models.

```julia
V,(_,EV,FV) = Lar.cuboidGrid([4,4],true);
W,(_,EW,FW) = Lar.cuboidGrid([3,5],true);
mycircle(r,n) = Lar.circle(r)(n)

data2d1 = (V,EV)
data2d2 = Lar.Struct([ Lar.t(2,2), Lar.r(pi/3), Lar.t(-1.5,-2.5), (W,EW) ])
data2d3 = Lar.Struct([ Lar.t(2,2), mycircle(2.5,16) ])
data2d4 = Lar.Struct([ Lar.t(3.5,3.5), mycircle(.25,16) ])
data2d5 = Lar.Struct([ Lar.t(5,3.5), mycircle(.5,16) ])
data2d6 = Lar.Struct([ Lar.t(5,3.5), mycircle(.25,16) ])

model2d = input_collection( [ data2d1, data2d2, data2d3, data2d4, data2d5, data2d6 ] )
V,EV = model2d
VV = [[k] for k in 1:size(V,2)];
using Plasm
Plasm.view( Plasm.numbering(.5)((V,[VV,EV])) )
```
Note that `V,EV` is not a cellular complex, since 1-cells intersect out of 0-cells.

# Example 3D

```julia
V,FV = Lar.sphere(2)([3,4])
EV = Lar.simplexFacets(FV)
mysphere = V,FV,EV

data3d1 = mysphere
data3d2 = Lar.Struct([ Lar.t(0,1,0), mysphere ])
data3d3 = Lar.Struct([ Lar.t(0,0.5,0), Lar.s(0.4,0.4,0.4), mysphere ])
data3d4 = Lar.Struct([ Lar.t(4,0,0), Lar.s(0.8,0.8,0.8), mysphere ])
data3d5 = Lar.Struct([ Lar.t(4,0,0), Lar.s(0.4,0.4,0.4), mysphere ])

model3d = input_collection([ data3d1, data3d2, data3d3, data3d4, data3d5 ])
V,FV,EV = model3d
VV = [[k] for k in 1:size(V,2)];
using Plasm
Plasm.view( Plasm.numbering(1)((V,[VV, EV])) )
```

Note that `V,FV,EV` is not a cellular complex, since 1-cells and
2-cells intersect out of 0-cells.

"""
function input_collection(data::Array)::Lar.LAR
	assembly = Lar.Struct(data)
	return Lar.struct2lar(assembly)
end

function boundingbox(vertices::Lar.Points)
   minimum = mapslices(x->min(x...), vertices, dims=2)
   maximum = mapslices(x->max(x...), vertices, dims=2)
   return minimum, maximum
end

function coordintervals(coord,bboxes)
	boxdict = OrderedDict{Array{Float64,1},Array{Int64,1}}()
	for (h,box) in enumerate(bboxes)
		key = box[coord,:]
		if haskey(boxdict,key) == false
			boxdict[key] = [h]
		else
			push!(boxdict[key], h)
		end
	end
	return boxdict
end

function boxcovering(bboxes, index, tree)
	covers = [[] for k=1:length(bboxes)]
	for (i,boundingbox) in enumerate(bboxes)
		extent = bboxes[i][index,:]
		iterator = IntervalTrees.intersect(tree, tuple(extent...))
		for x in iterator
			append!(covers[i],x.value)
		end
	end
	return covers
end

"""
	spaceindex(model::Lar.LAR)::Array{Array{Int,1},1}

Generation of *space indexes* for all ``(d-1)``-dim cell members of `model`.

*Spatial index* made by ``d`` *interval-trees* on
bounding boxes of ``sigma in S_{d−1}``. Spatial queries solved by
intersection of ``d`` queries on IntervalTrees generated by
bounding-boxes of geometric objects (LAR cells).

The return value is an array of arrays of `int`s, indexing cells whose
containment boxes are intersecting the containment box of the first cell.
According to Hoffmann, Hopcroft, and Karasick (1989) the worst-case complexity of
Boolean ops on such complexes equates the total sum of such numbers.

# Examples 2D

```
julia> V = hcat([[0.,0],[1,0],[1,1],[0,1],[2,1]]...);

julia> EV = [[1,2],[2,3],[3,4],[4,1],[1,5]];

julia> Sigma = Lar.spaceindex((V,EV))
5-element Array{Array{Int64,1},1}:
 [4, 5, 2]
 [1, 3, 5]
 [4, 5, 2]
 [1, 3, 5]
 [4, 1, 3, 2]
```

From `model2d` value, available in `?input_collection` docstring:

```julia
julia> Sigma =  spaceindex(model2d);
```

# Example 3D

```julia
model = model3d
Sigma =  spaceindex(model3d);
Sigma
```
"""
function spaceindex(model::Lar.LAR)::Array{Array{Int,1},1}
	V,CV = model[1:2]
	dim = size(V,1)
	cellpoints = [ V[:,CV[k]]::Lar.Points for k=1:length(CV) ]
	#----------------------------------------------------------
	bboxes = [hcat(Lar.boundingbox(cell)...) for cell in cellpoints]
	xboxdict = Lar.coordintervals(1,bboxes)
	yboxdict = Lar.coordintervals(2,bboxes)
	# xs,ys are IntervalTree type
	xs = IntervalTrees.IntervalMap{Float64, Array}()
	for (key,boxset) in xboxdict
		xs[tuple(key...)] = boxset
	end
	ys = IntervalTrees.IntervalMap{Float64, Array}()
	for (key,boxset) in yboxdict
		ys[tuple(key...)] = boxset
	end
	xcovers = Lar.boxcovering(bboxes, 1, xs)
	ycovers = Lar.boxcovering(bboxes, 2, ys)
	covers = [intersect(pair...) for pair in zip(xcovers,ycovers)]

	if dim == 3
		zboxdict = Lar.coordintervals(3,bboxes)
		zs = IntervalTrees.IntervalMap{Float64, Array}()
		for (key,boxset) in zboxdict
			zs[tuple(key...)] = boxset
		end
		zcovers = Lar.boxcovering(bboxes, 3, zs)
		covers = [intersect(pair...) for pair in zip(zcovers,covers)]
	end
	# remove each cell from its cover
	for k=1:length(covers)
		covers[k] = setdiff(covers[k],[k])
	end
	return covers
end

"""
	intersection(line1,line2)

Intersect two line segments in 2D, by computing the two line parameters of the intersection point.

The line segments intersect if both return parameters `α,β` are contained in the interval `[0,1]`.

# Example

```
julia> line1 = [[0.,0], [1,2]]
2-element Array{Array{Float64,1},1}:
 [0.0, 0.0]
 [1.0, 2.0]

julia> line2 = [[2.,0], [0,3]]
2-element Array{Array{Float64,1},1}:
 [2.0, 0.0]
 [0.0, 3.0]

julia> Lar.intersection(line1,line2)
(0.8571428571428571, 0.5714285714285714)
```
"""
function intersection(line1,line2)
	x1,y1,x2,y2 = vcat(line1...)
	x3,y3,x4,y4 = vcat(line2...)

	# intersect lines e1,e2
	# http://www.cs.swan.ac.uk/~cssimon/line_intersection.html
	det = (x4-x3)*(y1-y2)-(x1-x2)*(y4-y3)
	if det != 0.0
		a = 1/det
		b = [y1-y2 x2-x1; y3-y4 x4-x3]  # x1-x2 => x2-x1 bug in the source link !!
		c = [x1-x3; y1-y3]
		(β,α) = a * b * c
	else
		if (y1==y2) == (y3==y4) || (x1==x2) == (x3==x4) # segments collinear
			 return nothing
		else
			 # segments parallel: no intersection
			 return nothing
		end
	end
	return α,β
end

"""
	linefragments(V,EV,Sigma)

Compute the sequences of ordered parameters fragmenting each input lines.

Extreme parameter values (`0.0` and `1.0`) are included in each output line.
`Sigma` is the spatial index providing the subset of lines whose containment boxes intersect the box of each input line (given by `EV`).

```
julia> V = hcat([[0.,0],[1,0],[1,1],[0,1],[2,1]]...);

julia> EV = [[1,2],[2,3],[3,4],[4,1],[1,5]];

julia> Sigma = Lar.spaceindex((V,EV))
5-element Array{Array{Int64,1},1}:
 [4, 5, 2]
 [1, 3, 5]
 [4, 5, 2]
 [1, 3, 5]
 [4, 1, 3, 2]

julia> Lar.linefragments(V,EV,Sigma)
5-element Array{Any,1}:
 [0.0, 1.0]
 [0.0, 0.5, 1.0]
 [0.0, 1.0]
 [0.0, 1.0]
 [0.0, 0.5, 1.0]
```
"""
function linefragments(V,EV,Sigma)
	# remove the double intersections by ordering Sigma
	m = length(Sigma)
	sigma = map(sort,Sigma)
	reducedsigma = sigma ##[filter(x->(x > k), sigma[k]) for k=1:m]
	# pairwise parametric intersection
	params = Array{Float64,1}[[] for i=1:m]
	for h=1:m
		if sigma[h] ≠ []
			line1 = V[:,EV[h]]
			for k in sigma[h]
				line2 = V[:,EV[k]]
				out = Lar.intersection(line1,line2) # TODO: w interval arithmetic
				if out ≠ nothing
					α,β = out
					if 0<=α<=1 && 0<=β<=1
						push!(params[h], α)
						push!(params[k], β)
					end
				end
			end
		end
	end
	# finalize parameters of fragmented lines
	fragparams = []
	for line in params
		push!(line, 0.0, 1.0)
		line = sort(collect(Set(line)))
		push!(fragparams, line)
	end
	return fragparams
end

"""
	fragmentlines(model::Lar.LAR)::Lar.LAR

Pairwise *intersection* of 2D *line segments*.

# Example 2D

```julia
V,EV = model2d
W, EW = Lar.fragmentlines(model2d) # OK
using Plasm
Plasm.viewexploded(W,EW)(1.2,1.2,1.2)
```
"""
function fragmentlines(model)
	V,EV = model
	# acceleration via spatial index computation
	Sigma = Lar.spaceindex(model)
	# actual parametric intersection of each line with the close ones
	lineparams = Lar.linefragments(V,EV,Sigma)
	# initialization of local data structures
	vertdict = OrderedDict{Array{Float64,1},Array{Int,1}}()
	pairs = collect(zip(lineparams, [V[:,e] for e in EV]))
	vertdict = OrderedDict{Array{Float64,1},Int}()
	W = Array[]
	EW = Array[]
	k = 0
	# generation of intersection points
	for (params,linepoints) in pairs
		v1 = linepoints[:,1]
		v2 = linepoints[:,2]
		points = [ v1 + t*(v2 - v1) for t in params]   # !!!! loved !!
		vs = zeros(Int64,1,length(points))
		PRECISION = 8
		# identification via dictionary of points
		for (h,point) in enumerate(points)
			point = map(Lar.approxVal(PRECISION), point)
			if haskey(vertdict, point) == false
				k += 1
				vertdict[point] = k
				push!(W, point)
			end
			vs[h] = vertdict[point]
		end
		[push!(EW, [vs[k], vs[k+1]]) for k=1:length(vs)-1]
	end
	# normalization of output
	W,EW = hcat(W...),convert(Array{Array{Int64,1},1},EW)
	V,EV = Lar.congruence((W,EW))
	return V,EV
end
function fraglines(sx::Float64=1.2,sy::Float64=1.2,sz::Float64=1.2)
	function fraglines0(model)
		V,EV = Lar.fragmentlines(model)

		W = zeros(Float64, size(V,1), 2*length(EV))
		EW = Array{Array{Int64,1},1}()
		for (k,(v1,v2)) in enumerate(EV)
			if size(V,1)==2
				x,y = (V[:,v1] + V[:,v2]) ./ 2
				scx,scy = x*sx, y*sy
				t = [scx-x, scy-y]
			elseif size(V,1)==3
				x,y,z = (V[:,v1] + V[:,v2]) ./ 2
				scx,scy,scz = x*sx, y*sy, z*sz
				t = [scx-x, scy-y, scz-z]
			end
			W[:,2*k-1] = V[:,v1] + t
			W[:,2*k] = V[:,v2] + t
			push!(EW, [2*k-1, 2*k])
		end
		return W,EW
	end
	return fraglines0
end

"""
	congruence(model::Lar.LAR)::Lar.LAR
Graded bases of equivalence classes Ck (Uk ), with Uk = Xk /Rk for 0 ≤ k ≤ 2.

# Example

```julia
julia>
```
"""
function congruence(model)
	W,EW = model
	# congruent vertices
	balltree = NearestNeighbors.BallTree(W)
	r = 0.0000000001
	near = Array{Any}(undef, size(W,2))
	for k=1:size(W,2)
		near[k] = union([NearestNeighbors.inrange(balltree, W[:,k], r, true)]...)
	end
	near = map(sort,near)  # check !!!
	for k=1:size(W,2)
		W[:,k] = W[:,near[k][1]]
	end
	pointidx = [ near[k][1] for k=1:size(W,2) ]  # check !!
	invidx = OrderedDict(zip(1:length(pointidx), pointidx))
	V = [W[:,k] for k=1:length(pointidx)]
	# congruent edges
	EV = []
	for e in (EW)
		newedge = [invidx[e[1]],invidx[e[2]]]
		if newedge[1] !== newedge[2]
			push!(EV,newedge)
		end
	end
	EV = [EV[h] for h=1:length(EV) if length(EV[h])==2]
	EV = convert(Lar.Cells, EV)
	#W,EW = Lar.simplifyCells(V,EV)
	return hcat(V...),EV
end

#=
"""
	decomposition()::

Pairwise *intersection* in ``z = 0`` of *line segments* in ``σ ∪ I(σ)``, for each ``σ ∈ Sd−1``.

# Example 3D

```julia
V,FV,EV = model3d
model = model3d

```
"""
function decomposition(model::Lar.LAR)
	V,EV = model
	dim = size(V,1)
	spatialindex = Lar.spaceindex(model)

	function submanifoldmap(vs)
		centroid = [sum(vs[k,:]) for k=1:size(vs,1)]/size(vs,2)
		# u1, u2 always independent
		u1 = normalize( centroid - vs[:,1] )
		u2 = normalize( vs[:,2] - vs[:,1] )
		u3 = normalize(cross(u1, u2))
		# u1, u2, u3 orthonormal
		u1 = cross(u2, u3)
		T = Matrix{Float64}(LinearAlgebra.I, 4, 4)
		T[1:3,4] = - vs[:,1]
		R = Matrix{Float64}(LinearAlgebra.I, 4, 4)
		R[1:3, 1:3] = [u1 u2 u3]'
		return R*T  # roto-translation matrix
	end

	for (k,Sigma) in spatialindex
		sigma = Sigma[k]
		if dim == 3
			# transform Sigma s.t. Sigma[1], i.e. sigma, -> z=0
			vs = V[:, CV[sigma]]
			Q = submanifoldmap(vs)
			vq = Q * [vs; ones(1, size(vs,2))]
			v2d = vq[1:2,:]
		end
	end

end

"""
	Connection()::
Extraction of (X p , ∂p ), maximal connected components of Xd −1 (0 ≤ p ≤ h). d−1 d−1 +p

# Example

```julia
julia>
```
"""
function Connection

end

"""
	Bases()::
Computation of redundant cycle basis [∂d ] for each p-component, via TGW.

# Example

```julia
julia>
```
"""
function Bases

end

"""
	Boundaries()::
Accumulation into H += [o]p (hole-set) of outer boundary cycle from each [∂d+]p .

# Example

```julia
julia>
```
"""
function Boundaries

end

"""
	Containment()::
Computation of antisymmetric containment relation S between [o]p holes in H.

# Example

```julia
julia>
```
"""
function Containment

end

"""
	Reduction()::
Transitive R reduction of S and generation of forest of flat trees ⟨[od ]p , [∂d ]p ⟩.

# Example

```julia
julia>
```
"""
function Reduction

end

"""
	Adjoining()::
of roots [od ]r to (unique) outer cell, and non-roots [∂d+]q to container cells.

# Example

```julia
julia>
```
"""
function Adjoining

end

"""
	Assembling()::
Quasi-block-diagonal assembly of matrices relatives to isolated components [∂d ]p .

# Example

```julia
julia>
```
"""
function Assembling

end

"""
	Output()::
Global boundary map [∂d ] of A(Sd−1), and reconstruction of 0-chains of d-cells in Xd .

# Example

```julia
julia>
```
"""
function Output

end

=#


text:
"""
    fasta_records(filename::AbstractString)
    fasta_records(fasta_file::IO)

Reads records from a file in FASTA format.

If `filename::AbstractString`, first opens the file.

# Examples
```jldoctest
julia> fasta_file = IOBuffer(
\"\"\"
>seqA some description
QIKDLLVSSSTDLDTTLKMK
ILELPFASGDLSM
\"\"\");

julia> records = fasta_records(fasta_file)
1-element Array{BioSequences.FASTA.Record,1}:
 BioSequences.FASTA.Record:
   identifier: seqA
  description: some description
     sequence: QIKDLLVSSSTDLDTTLKMKILELPFASGDLSM
```
"""
function fasta_records(fasta_file::IO)
    records = Vector{BioSequences.FASTA.Record}()

    reader = FASTA.Reader(fasta_file)
    for record in reader
        push!(records, record)
    end
    close(reader)

    records
end
fasta_records(filename::AbstractString) = fasta_records(open(filename, "r"))

"""
    signature(record, k, n)

Generates `n` minhash signatures from `k`-shingles of a biological sequence
in a FASTA record.
"""
function signature(record::BioSequences.FASTA.Record, k::Int, n::Int)
    id = FASTA.identifier(record)
    sig = signature(String(sequence(record)), k, n)
    Signature(id, sig)
end

"""
    fasta_signatures(filename::AbstractString, k, n)
    fasta_signatures(fasta_file::IO, k, n)

Reads protein sequences from `filename` in FASTA format and generates `n`
minhash signatures from `k`-shingles.

# Examples
```jldoctest
julia> fasta_file = IOBuffer(
\"\"\"
>seqA some description
QIKDLLVSSSTDLDTTLKMK
ILELPFASGDLSM
>seqB
VLMALGMTDLFIPSANLTG*
\"\"\");

julia> sigs = fasta_signatures(fasta_file)
2-element Array{Signature,1}:
 Signature{String,Array{UInt32,1}}("seqA", UInt32[0x05094006, 0x0afea36a])
 Signature{String,Array{UInt32,1}}("seqB", UInt32[0x1380e49b, 0x17854bad])
```
"""
function fasta_signatures(fasta_file::IO, k::Int, n::Int)
    signatures = Vector{Signature}()

    reader = FASTA.Reader(fasta_file)
    for record in reader
        push!(signatures, LSMH.signature(record, k, n))
    end
    close(reader)

    signatures
end
function fasta_signatures(filename::AbstractString, k::Int, n::Int)
    fasta_signatures(open(filename, "r"), k, n)
end


text:
using LinearAlgebra, LinearMaps

export BlockMap

# we probably don't need T for block map either, but need some upstream
# modification in LinearMaps.jl
struct BlockMap{T, GT <: AbstractBlock} <: LinearMap{T}
    block::GT

    BlockMap(::Type{T}, block::GT) where {T, GT <: AbstractBlock} = new{T, GT}(block)
end

BlockMap(block::AbstractBlock) = BlockMap(ComplexF64, block)

function Base.show(io::IO, A::BlockMap{T}) where T
    println(io, "Quantum Circuit Block as LinearMap{$T}")
    print(io, A.block)
end

function LinearMaps.A_mul_B!(y::AbstractVecOrMat{T1}, op::BlockMap{T2}, x::AbstractVecOrMat{T3}) where {T1, T2, T3}
    error("Do not use BlockMap on different precisions, this is not supported.")
end

function LinearMaps.A_mul_B!(y::AbstractVecOrMat{T}, op::BlockMap{T}, v::AbstractVecOrMat{T}) where T
    copyto!(y, v)
    apply!(ArrayReg(y), op.block)
    return y
end

content(x::BlockMap) = x.block

# NOTE: do not overload operator on different element type
#       we don't support this to error when there is performance
#       issue due to conversion.
Base.:(*)(op::BlockMap{T}, v::VT) where {T, VT <: AbstractVecOrMat{T}} =
    mul!(similar(v), op, v)

"""
    opnorm(A::BlockMap, p::Real=2)

[`opnorm`](@ref) for quantum circuit blocks.
"""
LinearAlgebra.opnorm(A::BlockMap, p::Real=2) = opnorm(mat(A.block), p)
LinearAlgebra.ishermitian(A::BlockMap) = ishermitian(A.block)
LinearAlgebra.issymmetric(A::BlockMap) = issymmetric(A.block)
Base.size(bm::BlockMap) = (1 << nqubits(bm.block), 1 << nqubits(bm.block))


text:
"""
    dircos(incl, decl)

Compute the three cartesian components of the unit vector aligned with a field vector defined by the two angles inclination `incl` and declination `decl`.

# Example
```julia-repl
julia> dircos(56.0, 0.0)
(0.5591929034707468, 0.0, 0.8290375725550417)
```
"""
function dircos(incl::T, decl::T) where T<:Real
    xincl = incl * pi / 180.0
    xdecl = decl * pi / 180.0
    a = cos(xincl) * cos(xdecl)
    b = cos(xincl) * sin(xdecl)
    c = sin(xincl)
    return a, b, c
end

max_stars_repo_name: ruboerner/GravMag

max_issues_repo_name: ruboerner/GravMag

max_forks_repo_name: ruboerner/GravMag


text:
using DiffEqFlux, Flux, Test, OrdinaryDiffEq, CuArrays
using Statistics
#= using Plots =#

## True Solution
u0 = [2.; 0.] |> gpu
datasize = 30
tspan = (0.0,25.0)

function trueODEfunc(du,u,p,t)
    true_A = cu([-0.1 2.0; -2.0 -0.1])
    du .= ((u.^3)'true_A)'
end

true_prob = ODEProblem(trueODEfunc, u0,tspan)

true_sol = solve(true_prob,BS3(),saveat=range(tspan[1],tspan[2],length=datasize))

#= true_sol_plot = solve(true_prob,Tsit5()) =#
#= plot(true_sol_plot) =#

## Neural ODE
dudt = Chain(Dense(2,50,tanh),Dense(50,2)) |> gpu

function ODEfunc(du,u,p,t)
    du .= Flux.data(dudt(u))
end

pred_prob = ODEProblem(ODEfunc, u0,tspan)
pred_sol = solve(pred_prob,BS3(),saveat=range(tspan[1],tspan[2],length=datasize))

## Loss
l1_loss(pred,target) = mean(abs.(pred-target))
l1_loss(pred_sol,true_sol)


text:
using LinearAlgebra: norm

struct QuadDEO end

"""
    quaddeo(f::Function, ω::Real, θ::Real, a::Real, b::Real;
            h0::Real=one(ω)/5, maxlevel::Integer=12,
            atol::Real=zero(ω),
            rtol::Real=atol>0 ? zero(atol) : sqrt(eps(typeof(atol))))

Numerically integrate `f(x)` over an arbitrary interval [a, b] and return the
integral value `I` and an estimated error `E`. The `quaddeo` function is
specialized for the decaying oscillatory integrands,

    f(x) = f₁(x)sin(ωx + θ),

where `f₁(x)` is a decaying algebraic function. `ω` and `θ` are the frequency
and the phase of the oscillatory part of the integrand. If the oscillatory part
is `sin(ωx)`, then `θ = 0.0`; if it is `cos(ωx)` instead, then `θ = π/(2ω)`.
The `E` is not exactly equal to the difference from the true value. However,
one can expect that the integral value `I` is converged if
`E <= max(atol, rtol*norm(I))` is true. Otherwise, the obtained `I` would be
unreliable; the number of repetitions exceeds the `maxlevel` before converged.
Optionally, one can divide the integral interval [a, b, c...], which returns
`∫ f(x)dx in [a, b] + ∫f(x)dx in [b, c[1]] + ⋯`.  It is worth noting that each
endpoint allows discontinuity or singularity.

The integrand `f` can also return any value other than a scalar, as far as
`+`, `-`, multiplication by real values, and `LinearAlgebra.norm`, are
implemented. For example, `Vector` or `Array` of numbers are acceptable
although, unfortunately, it may not be very performant.

# Examples
```jldoctest
julia> using DoubleExponentialFormulas

julia> using LinearAlgebra: norm

julia> f(x) = sin(x)/x;

julia> I, E = quaddeo(f, 1.0, 0.0, 0.0, Inf);

julia> I ≈ π/2
true

julia> E ≤ sqrt(eps(Float64))*norm(I)
true
```
"""
quaddeo(f::Function, ω::Real, θ::Real, a::Real, b::Real;
        h0::Real=one(ω)/5, maxlevel::Integer=12,
        atol::Real=zero(ω), rtol::Real=atol>0 ? zero(atol) : sqrt(eps(typeof(atol)))) =
    quaddeo_entrance(f, ω, θ, a, b, h0, maxlevel, atol, rtol)

function quaddeo(f::Function, ω::Real, θ::Real, a::Real, b::Real, c::Real...;
                 h0::Real=one(ω)/5, maxlevel::Integer=12,
                 atol::Real=zero(ω), rtol::Real=atol>0 ? zero(atol) : sqrt(eps(typeof(atol))))
    bc = (b, c...)
    n = length(bc)
    _atol = atol/n
    I, E = quaddeo_entrance(f, ω, θ, a, b, h0, maxlevel, atol, rtol)
    for i in 2:n
        dI, dE = quaddeo_entrance(f, ω, θ, bc[i-1], bc[i], h0, maxlevel, atol, rtol)
        I += dI
        E += dE
    end
    return I, E
end

function quaddeo_entrance(f, ω, θ, a, b, h0, maxlevel, atol, rtol)
    @assert maxlevel > 0
    @assert h0 > 0
    _ω, _θ, _a, _b, _h0 = float.(promote(ω, θ, a, b, h0))
    if a > b
        I, E = quaddeo_main(f, _ω, _θ, _b, _a, _h0, maxlevel, atol, rtol)
        return -I, E
    else
        return quaddeo_main(f, _ω, _θ, _a, _b, _h0, maxlevel, atol, rtol)
    end
end
function quaddeo_main(f, ω::T, θ::T, a::T, b::T, h0::T, maxlevel, atol, rtol) where {T<:AbstractFloat}
    if a == b
        M = π/h0
        δ = θ/M
        t0 = zero(T) + δ
        x0, w0 = samplepoint(QuadDEO, t0)
        Σ = f(M*x0)*w0
        I = M*h0*Σ
        return zero(typeof(I)), zero(eltype(I))
    end

    if a == -Inf && b == Inf
        _atol = atol/2
        I⁻, E⁻ = quaddeo_main(f, ω, θ, a, zero(a), h0, maxlevel, _atol, rtol)
        I⁺, E⁺ = quaddeo_main(f, ω, θ, zero(b), b, h0, maxlevel, _atol, rtol)
        return I⁺ + I⁻, E⁺ + E⁻
    elseif b == Inf
        if a == 0
            return quaddeo_main(f, ω, θ, h0, maxlevel, atol, rtol)
        else
            return quaddeo_main(u -> f(u + a), ω, θ, h0, maxlevel, atol, rtol)
        end
    elseif a == -Inf
        if b == 0
            return quaddeo_main(u -> f(-u), ω, θ, h0, maxlevel, atol, rtol)
        else
            return quaddeo_main(u -> f(-u + b), ω, θ, h0, maxlevel, atol, rtol)
        end
    else
        _atol = atol/2
        if a < 0 && b ≤ 0
            Ia, Ea = quaddeo_main(f, ω, θ, T(-Inf), a, h0, maxlevel, _atol, rtol)
            Ib, Eb = quaddeo_main(f, ω, θ, T(-Inf), b, h0, maxlevel, _atol, rtol)
            return Ib - Ia, Ea + Eb
        elseif a ≥ 0 && b > 0
            Ia, Ea = quaddeo_main(f, ω, θ, a, T(Inf), h0, maxlevel, _atol, rtol)
            Ib, Eb = quaddeo_main(f, ω, θ, b, T(Inf), h0, maxlevel, _atol, rtol)
            return Ia - Ib, Ea + Eb
        else
            I⁻, E⁻ = quaddeo_main(f, ω, θ, a, zero(a), h0, maxlevel, _atol, rtol)
            I⁺, E⁺ = quaddeo_main(f, ω, θ, zero(b), b, h0, maxlevel, _atol, rtol)
            return I⁺ + I⁻, E⁺ + E⁻
        end
    end
end
function quaddeo_main(f, ω::T, θ::T, h0::T, maxlevel, atol, rtol) where {T<:AbstractFloat}
    I = integrate(QuadDEO, f, ω, θ, h0, atol, rtol)
    E = zero(eltype(I))
    for level in 1:maxlevel
        h = h0/2^level
        prevI = I
        I = integrate(QuadDEO, f, ω, θ, h, atol, rtol)
        E = estimate_error(T, prevI, I)
        tol = max(norm(I)*rtol, atol)
        !(E > tol) && break
    end
    return I, E
end

function integrate(::Type{QuadDEO}, f, ω, θ, h, atol, rtol)
    M = π/(h*ω)
    δ = θ/M
    t0 = zero(δ) + δ
    x0, w0 = samplepoint(QuadDEO, t0)
    Σ = f(M*x0)*w0
    I = M*Σ*h
    chunklen = 10

    k = 1
    while true
        Σ, loop_done = integrate_chunk⁺(QuadDEO, f, M, δ, k, h, Σ, chunklen)
        k += chunklen
        prevI = I
        I = M*h*Σ
        loop_done && break
        dI = norm(prevI - I)
        tol = max(norm(I)*rtol, atol)
        dI ≤ tol && break
    end

    k = -1
    while true
        Σ, loop_done = integrate_chunk⁻(QuadDEO, f, M, δ, k, h, Σ, chunklen)
        k -= chunklen
        prevI = I
        I = M*h*Σ
        loop_done && break
        dI = norm(prevI - I)
        tol = max(norm(I)*rtol, atol)
        dI ≤ tol && break
    end
    return I
end

function integrate_chunk⁺(::Type{QuadDEO}, f, M, δ, k, h, Σ, chunklen)
    loop_done = false
    for _ in 1:chunklen
        t = k*h + δ
        if isinf(sinh(t))
            loop_done = true
            break
        end

        x, w = samplepoint(QuadDEO, t)
        dΣ = f(M*x)*w
        if all(isnan.(dΣ))
            loop_done = true
            break
        end

        Σ += dΣ
        k += 1
    end
    return Σ, loop_done
end

function integrate_chunk⁻(::Type{QuadDEO}, f, M, δ, k, h, Σ, chunklen)
    loop_done = false
    for _ in 1:chunklen
        t = k*h + δ
        if isinf(sinh(t))
            loop_done = true
            break
        end

        x, w = samplepoint(QuadDEO, t)
        if x == 0
            # avoid end-point singularity
            loop_done = true
            break
        end

        dΣ = f(M*x)*w
        if all(isnan.(dΣ))
            loop_done = true
            break
        end

        Σ += dΣ
        k -= 1
    end
    return Σ, loop_done
end

function samplepoint(::Type{QuadDEO}, t::T) where {T<:AbstractFloat}
    # The following equation is used for the variable conversion.
    #   ϕ(t) = (t + sqrt(t² + exp(π - K*cosh(t))))/2
    # Ref. (Japanese) : http://www.kurims.kyoto-u.ac.jp/~ooura/intdefaq-j.html
    #
    # NOTE: might be worth checking this?
    #       https://doi.org/10.1016/S0377-0427(99)00223-X
    K = T(2π)
    A = exp(π - K*cosh(t))
    B = sqrt(t^2 + A)
    ϕ  = (t + B)/2
    ϕ′ = T(0.5) + (2t - K*A*sinh(t))/(4*B)
    return ϕ, ϕ′
end


text:
struct GUI
  id::UUID
  widgets::Vector{Widget}
  expression::Expr
end

macro gui(expr)
  if expr.head != :for
    error(
      "@gui syntax is @gui for ",
      " [<variable>=<domain>,]... <expression> end"
    )
  end

  is_gui = true
  quote_widgets, cur_block, cur_bindings, cur_symbols = _manipulate_outer(expr, is_gui)

  return quote
    cur_widgets = make_gui_list($(quote_widgets...))
    cur_listener = $(esc(make_gui_block(cur_block, cur_symbols)))

    cur_id = _manipulate_inner(cur_widgets, cur_listener, $(is_gui))
    nothing
  end
end

macro manipulate(expr)
  if expr.head != :for
    error(
      "@manipulate syntax is @manipulate for ",
      " [<variable>=<domain>,]... <expression> end"
    )
  end

  is_gui = false
  quote_widgets, cur_block, cur_bindings, cur_symbols = _manipulate_outer(expr, is_gui)

  return quote
    cur_widgets = make_gui_list($(quote_widgets...))
    cur_listener = $(esc(make_gui_block(cur_block, cur_symbols)))

    cur_id = _manipulate_inner(cur_widgets, cur_listener, $(is_gui))
    nothing
  end
end

export @gui, @manipulate


text:
#################################################################
#
#    Common definitions for tests
#
#################################################################

using Base.Test

reload("ReverseDiffSource") # instead of 'using' to retest without exiting Julia
m = ReverseDiffSource

## removes linenumbers from expression to ease comparisons
function striplinenumbers(ex::Expr)
    args = Any[]
    for a in ex.args
        isa(a, LineNumberNode) && continue
        isa(a, Expr) && a.head==:line && continue
        push!(args, isa(a,Expr) ? striplinenumbers(a) : a )
    end
    Expr(ex.head, args...)
end

## first order derivation checking function

#  Error thresholds  #####
DIFF_DELTA = 1e-9
ERROR_THRESHOLD = 2e-2

good_enough(x,y) = isfinite(x) ? (abs(x-y) / max(ERROR_THRESHOLD, abs(x))) < ERROR_THRESHOLD : isequal(x,y)
good_enough(t::Tuple) = good_enough(t[1], t[2])

#  Compares numerical gradient to automated gradient
function compare( ex::Union{Expr, Symbol},
                  x0::Union{Float64, Vector{Float64}, Matrix{Float64}} )
    nx = length(x0)

    if isa(ex, Expr)
        ex2 = m.rdiff( ex, x=typeof(x0) )
        dfunc(x0) = eval( :(let x = $x0 ; $ex2 ; end) )
    else
        dfunc = m.rdiff( eval(ex), (typeof(x0),) )
    end

    l0, grad0 = dfunc(x0)
    if ndims(x0) == 0  # scalar
        grad1 = ( dfunc( x0 + DIFF_DELTA)[1] - l0 ) / DIFF_DELTA
    else # vector and matrices
        grad1 = zeros(size(grad0))
        for i in 1:nx  # i=1
            x1 = copy(x0)
            x1[i] += DIFF_DELTA
            grad1[i] = ( dfunc(x1)[1] - l0 ) / DIFF_DELTA
        end
    end

    if !all(good_enough, zip([grad0;], [grad1;]))
        rg0 = map(x -> round(x,5), grad0)
        rg1 = map(x -> round(x,5), grad1)
        error("\nGradient false at x=$x0, expected $rg1, got $rg0")
    end
end

macro compare(ex, x0)
    if isa(x0, Symbol) || isa(x0, Expr)
        x0 = eval(x0)
    end
    compare(ex, eval(x0))
end

## variables of different dimensions for testing
v0ref = 2.
v1ref = [2., 3, 0.1, 0, -5]
v2ref = [-1. 3 0 ; 0 5 -2]


text:
using Klimakoffer

world = Klimakoffer.read_geography(joinpath(@__DIR__, "..", "input", "world", "The_World128x65.dat"), 128, 65)

rand_long = rand(128:2:2048)

upscaled_world = Klimakoffer.nn_interpolation(world, rand_long)

(nlong,nlat) = size(upscaled_world)

landcover = 100*count(i->(i==1), upscaled_world)/(nlong*nlat)


text:
function sendto(p::Int; args...)
    for (nm, val) in args
        @spawnat(p, Core.eval(Main, Expr(:(=), nm, val)))
    end
end


text:
using Documenter: doctest
doctest((@__DIR__), Module[]; fix = true)


text:
using Nextion
using Test

@testset "scroll" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    nexPage = NexPage(nexSerial, Name("pg_scroll"), pid=PageID(6))
    show(nexPage)

    nexScrollText = NexScrollText(nexSerial, Name("g0"), pid=PageID(1))
    msg = "Hello Nextion!"
    nexScrollText.text = msg
    @test nexScrollText.text == msg

    sleep(10)

    close(nexSerial)
    @test true

end