text:
export calc_k, calc_λ, calc_κ
export calc_NA
export fresnel_number

"""
    calc_NA(focal_length, diameter[, n])

Calculate the numerical aperture of a system with `focal_length`
and `diameter` of the lens. Per default `n=1` is the refractive index.

 # Examples
```julia-repl
julia> calc_NA(100e-3, 200e-3)
1.0

julia> calc_NA(100e-3, 200e-3, 1.33)
1.33
```
"""
function calc_NA(focal_length, diameter, n=1)
    return n * diameter / 2 / focal_length
end

"""
    calc_k(λ, n=1)

Calculate the value of the (angular) wave number with vacuum wavelength `λ` in
medium with refractive index `n`.
It holds: \$ k = \\kappa \\cdot 2 \\pi \$
"""
function calc_k(λ, n=1)
    return 2 * oftype(λ, π) * n / λ
end

"""
    calc_κ(λ, n=1)

Calculate the value of the (non angular) wave number with vacuum wavelength `λ` in
medium with refractive index `n`.
It holds: \$ \\kappa = \\frac{k}{2 \\pi} \$
"""
function calc_κ(λ, n=1)
    return calc_k(λ, n) / 2 / oftype(λ, π)
end

"""
    calc_λ(k, n=1)

Calculate the vacuum wavelength `λ` from the angular wave number `k` in
medium with refractive index `n`.
It holds: \$ \\lambda = \\frac{n \\cdot 2 \\pi}{k} \$
"""
function calc_λ(k, n=1)
    return 2 * oftype(k, π) * n / k
end

"""
    fresnel_number(a, L, λ=λ0)

Calculate the Fresnel number where `a` is the characteristic size,
`L` the distance from screen to aperture and `λ` the wavelength.
"""
function fresnel_number(a, L, λ=λ0)
    return a^2 / L / λ
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
include("./helpers.jl")

include("./ExploreExpression.jl")
include("./React.jl")
include("./Notebook.jl")
include("./WorkspaceManager.jl")

# TODO: test PlutoRunner functions like:
# - from_this_notebook

# TODO: test HTTP
# We could use NodeJS.jl to run these tests without a browser
# We could also do real browser testing, but that seems difficult

ext: jl


text:
## load
using Test, SpecialFunctions, VortexDistributions
using FourierGPE, Plots

## Initialize simulation
L = (18.0,18.0)
N = (256,256)
sim = Sim(L,N)
@unpack_Sim sim

## make initial state
x,y = X
ψi = 100*one.(x.*y') |> complex
ϕi = kspace(ψi,sim)
@pack_Sim! sim
dx,dy = diff(x)[1],diff(y)[1]

## plot
heatmap(x,y,abs2.(ψi))

## call
k = LinRange(0.01,10.0,100) |> collect
Eki = ikespectrum(k,ψi,X,K)

## plot
plot(k,Eki)

## call

Ek = kespectrum(k,ψi,X,K)

## plot
plot(k,Ek)

## what does autocorrelate look like?
x,y = X; kx,ky = K
dx,dy = diff(x)[1],diff(y)[1]
psiac = autocorrelate(ψi,X,K)
psiac2 = autocorrelate(ψi,X,K,periodic=true)

## plot
heatmap(abs2.(psiac2))

##
psiac2


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

JLLWrappers.@generate_main_file_header("Eigen")
JLLWrappers.@generate_main_file("Eigen", UUID("bc6bbf8a-a594-5541-9c57-10b0d0312c70"))
end  # module Eigen_jll


text:
mutable struct MWSingleLayerField3D{K}
    wavenumber::K
end

mutable struct MWDoubleLayerField3D{K}
    wavenumber::K
end

"""
    MWSingleLayerField3D(wavenumber = error())

Create the single layer near field operator, for use with `potential`.
"""
MWSingleLayerField3D(;wavenumber = error("Missing arg: `wavenumber`")) = MWSingleLayerField3D(wavenumber)
MWDoubleLayerField3D(;wavenumber = error("Missing arg: `wavenumber`")) = MWDoubleLayerField3D(wavenumber)

const MWField3D = Union{MWSingleLayerField3D,MWDoubleLayerField3D}
quaddata(op::MWField3D,rs,els) = quadpoints(rs,els,(2,))
quadrule(op::MWField3D,refspace,p,y,q,el,qdata) = qdata[1,q]
function kernelvals(op::MWField3D,y,p)

        k = op.wavenumber
        r = y - cartesian(p)
        R = norm(r)

        kr = k * R
        expn = exp(-im * k * R)
        green = expn / (4pi*R)
        gradgreen = -(im*k + 1/R) * green / R * r

        krn = (trans=r, dist=R, green=green, gradgreen=gradgreen)
end

function integrand(op::MWSingleLayerField3D, krn, y, fp, p)

    γ = im * op.wavenumber

    j = fp.value
    ρ = fp.divergence

    G  = krn.green
    ∇G = krn.gradgreen

    -γ*G*j + ∇G*ρ/γ
end

function integrand(op::MWDoubleLayerField3D, krn, y, fp, p)

    j = fp.value
    ∇G = krn.gradgreen

    ∇G × j
end


text:
# ---
# title: 749. Contain Virus
# id: problem749
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search
# link: <https://leetcode.com/problems/contain-virus/description/>
# hidden: true
# ---
#
# A virus is spreading rapidly, and your task is to quarantine the infected area
# by installing walls.
#
# The world is modeled as a 2-D array of cells, where `0` represents uninfected
# cells, and `1` represents cells contaminated with the virus. A wall (and only
# one wall) can be installed **between any two 4-directionally adjacent cells**
# , on the shared boundary.
#
# Every night, the virus spreads to all neighboring cells in all four directions
# unless blocked by a wall. Resources are limited. Each day, you can install
# walls around only one region -- the affected area (continuous block of
# infected cells) that threatens the most uninfected cells the following night.
# There will never be a tie.
#
# Can you save the day? If so, what is the number of walls required? If not, and
# the world becomes fully infected, return the number of walls used.
#
# **Example 1:**
#
#
#
#     Input: grid =
#     [[0,1,0,0,0,0,0,1],
#      [0,1,0,0,0,0,0,1],
#      [0,0,0,0,0,0,0,1],
#      [0,0,0,0,0,0,0,0]]
#     Output: 10
#     Explanation:
#     There are 2 contaminated regions.
#     On the first day, add 5 walls to quarantine the viral region on the left. The board after the virus spreads is:
#
#     [[0,1,0,0,0,0,1,1],
#      [0,1,0,0,0,0,1,1],
#      [0,0,0,0,0,0,1,1],
#      [0,0,0,0,0,0,0,1]]
#
#     On the second day, add 5 walls to quarantine the viral region on the right. The virus is fully contained.
#
#
# **Example 2:**
#
#
#
#     Input: grid =
#     [[1,1,1],
#      [1,0,1],
#      [1,1,1]]
#     Output: 4
#     Explanation: Even though there is only one cell saved, there are 4 walls built.
#     Notice that walls are only built on the shared boundary of two different cells.
#
#
# **Example 3:**
#
#
#
#     Input: grid =
#     [[1,1,1,0,0,0,0,0,0],
#      [1,0,1,0,1,1,1,1,1],
#      [1,1,1,0,0,0,0,0,0]]
#     Output: 13
#     Explanation: The region on the left only builds two new walls.
#
#
# **Note:**
#
#   1. The number of rows and columns of `grid` will each be in the range `[1, 50]`.
#   2. Each `grid[i][j]` will be either `0` or `1`.
#   3. Throughout the described process, there is always a contiguous viral region that will infect **strictly more** uncontaminated squares in the next round.
#
#
## @lc code=start
using LeetCode

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


text:
using TLA+
using Test

@testset "TLA+.jl" begin
    # Write your tests here.
end


text:
"""
    plot(arg1::Array; kwargs...)

reads (x,y) pairs from files [or standard input] and generates PostScript code that will plot lines,
polygons, or symbols at those locations on a map.

Full option list at [`psxy`]($(GMTdoc)plot.html)

Parameters
----------

- **A** | **steps** | **straight_lines** :: [Type => Str]

	By default, geographic line segments are drawn as great circle arcs.
	To draw them as straight lines, use this option.
    ($(GMTdoc)plot.html#a)
- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **D** | **shift** | **offset** :: [Type => Str]

    Offset the plot symbol or line locations by the given amounts dx/dy.
    ($(GMTdoc)plot.html#d)
- **E** | **error** | **error_bars** :: [Type => Str]

    Draw symmetrical error bars.
    ($(GMTdoc)plot.html#e)
- **F** | **conn** | **connection** :: [Type => Str]

    Alter the way points are connected
    ($(GMTdoc)plot.html#f)
- **G** | **fill** | **markerfacecolor** | **MarkerFaceColor** | **mc** :: [Type => Str]

    Select color or pattern for filling of symbols or polygons. BUT WARN: the alias 'fill' will set the
    color of polygons OR symbols but not the two together. If your plot has polygons and symbols, use
    'fill' for the polygons and 'markerfacecolor' for filling the symbols. Same applyies for W bellow
    ($(GMTdoc)plot.html#g)
- **I** | **intens** :: [Type => Str | number]

    Use the supplied intens value (in the [-1 1] range) to modulate the fill color by simulating illumination.
    ($(GMTdoc)plot.html#i)
- **L** | **close** | **polygon** :: [Type => Str]

    Force closed polygons.
    ($(GMTdoc)plot.html#l)
- **N** | **no_clip** | **noclip** :: [Type => Str or []]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)plot.html#n)
- $(GMT.opt_P)
- **S** | **symbol** | **marker** | **Marker** :: [Type => Str]

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    ($(GMTdoc)plot.html#s)
    Alternatively select a sub-set of symbols using the aliases: **marker** or **Marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
	+ **y**, **y_dash**

    and select their sizes with the **markersize** or **size** keyword [default is 7p].
    The marker size can be a scalar or a vector with same size numeber of rows of data. Units are
	points unless specified otherwise with (for example for cm) *par=(PROJ_LENGTH_UNIT="c")*

- **W** | **pen** | **markeredgecolor** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)
    WARNING: the pen attributes will set the pen of polygons OR symbols but not the two together.
    If your plot has polygons and symbols, use **W** or **pen** for the polygons and
    **markeredgecolor** for filling the symbols. Similar to S above.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)

- **Z** | **level** :: [Type => Str | NamedTuple]	`Arg = value|file[+f|+l] | (data=Array|Number, outline=_, fill=_)`

    Paint polygons after the level given as a cte or a vector with same size of number of polygons. Needs a color map.
    ($(GMTdoc)plot.html#z)

- **axis** | **aspect** :: [Type => Str]

    When equal to "equal" makes a square plot.
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function plot(arg1; first=true, kw...)
	common_plot_xyz("", cat_1_arg(arg1), "plot", first, false, kw...)
end
plot!(arg1; kw...) = plot(arg1; first=false, kw...)

function plot(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	common_plot_xyz("", cat_2_arg2(rang, [f(x) for x in rang]), "plot", first, false, kw...)
end
plot!(f::Function, rang=nothing; kw...) = plot(f, rang; first=false, kw...)

function plot(f1::Function, f2::Function, range_t=nothing; first=true, kw...)	# Parametric version
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; is3D=false, kw...), "plot", first, false, kw...)
end
plot!(f1::Function, f2::Function, range_t=nothing; kw...) = plot(f1, f2, range_t; first=false, kw...)

plot(cmd0::String="", arg1=nothing; kw...)  = common_plot_xyz(cmd0, arg1, "plot", true, false, kw...)
plot!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "plot", false, false, kw...)
plot(arg1, arg2; kw...)  = common_plot_xyz("", cat_2_arg2(arg1, arg2), "plot", true, false, kw...)
plot!(arg1, arg2; kw...) = common_plot_xyz("", cat_2_arg2(arg1, arg2), "plot", false, false, kw...)
plot(arg1::Number, arg2::Number; kw...)  = common_plot_xyz("", cat_2_arg2([arg1], [arg2]), "plot", true, false, kw...)
plot!(arg1::Number, arg2::Number; kw...) = common_plot_xyz("", cat_2_arg2([arg1], [arg2]), "plot", false, false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
plot3d(arg1::Array; kwargs...)

reads (x,y,z) triplets from files [or standard input] and generates PostScript code that will plot lines,
polygons, or symbols at those locations in 3-D.

Full option list at [`plot3d`]($(GMTdoc)plot3d.html)

Parameters
----------

- **A** | **steps** | **straight_lines** :: [Type => Str]

    By default, geographic line segments are drawn as great circle arcs. To draw them as straight lines, use this option.
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **C** | **color** :: [Type => Str]

    Give a CPT or specify -Ccolor1,color2[,color3,...] to build a linear continuous CPT from those colors automatically.
    ($(GMTdoc)plot3d.html#c)
- **D** | **offset** :: [Type => Str]

    Offset the plot symbol or line locations by the given amounts dx/dy.
    ($(GMTdoc)plot3d.html#d)
- **E** | **error_bars** :: [Type => Str]

    Draw symmetrical error bars.
    ($(GMTdoc)plot3d.html#e)
- **F** | **conn** | **connection** :: [Type => Str]

    Alter the way points are connected
    ($(GMTdoc)plot3d.html#f)
- **G** | **fill** | **markerfacecolor** | **MarkerFaceColor** :: [Type => Str]

    Select color or pattern for filling of symbols or polygons. BUT WARN: the alias 'fill' will set the
    color of polygons OR symbols but not the two together. If your plot has polygons and symbols, use
    'fill' for the polygons and 'markerfacecolor' for filling the symbols. Same applyies for W bellow
    ($(GMTdoc)plot3d.html#g)
- **I** | **intens** :: [Type => Str or number]

    Use the supplied intens value (in the [-1 1] range) to modulate the fill color by simulating illumination.
    ($(GMTdoc)plot3d.html#i)
- **L** | **closed_polygon** :: [Type => Str]

    Force closed polygons.
    ($(GMTdoc)plot3d.html#l)
- **N** | **no_clip** :: [Type => Str | []]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)plot3d.html#n)
- $(GMT.opt_P)
- **S** | **symbol** | **marker** | **Marker** :: [Type => Str]

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    ($(GMTdoc)plot3d.html#s)
    Alternatively select a sub-set of symbols using the aliases: **marker** or **Marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
    + **y**, **y_dash**
- **W** | **pen** | **line_attribs** | **markeredgecolor** | **MarkerEdgeColor** :: [Type => Str]
    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot3d.html#w)
    WARNING: the pen attributes will set the pen of polygons OR symbols but not the two together.
    If your plot has polygons and symbols, use **W** or **line_attribs** for the polygons and
    **markeredgecolor** or **MarkerEdgeColor** for filling the symbols. Similar to S above.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)

Example:

    plot3d(x -> sin(x)*cos(10x), y -> sin(y)*sin(10y), z -> cos(z), 0:pi/100:pi, show=true, aspect3=:equal)
"""
plot3d(arg1; kw...)  = common_plot_xyz("", cat_1_arg(arg1), "plot3d", true, true, kw...)
plot3d!(arg1; kw...) = common_plot_xyz("", cat_1_arg(arg1), "plot3d", false, true, kw...)

plot3d(cmd0::String="", arg1=nothing; kw...)  = common_plot_xyz(cmd0, arg1, "plot3d", true, true, kw...)
plot3d!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "plot3d", false, true, kw...)

# ------------------------------------------------------------------------------------------------------
function plot3d(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; first=true, kw...)
	common_plot_xyz("", hcat(arg1[:], arg2[:], arg3[:]), "plot3d", first, true, kw...)
end
plot3d!(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; kw...) = plot3d(arg1, arg2, arg3; first=false, kw...)

function plot3d(f1::Function, f2::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; kw...), "plot3d", first, true, kw...)
end
plot3d!(f1::Function, f2::Function, range_t=nothing; kw...) = plot3d(f1, f2, range_t; first=false, kw...)

function plot3d(f1::Function, f2::Function, f3::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_3f(f1, f2, f3, range_t; kw...), "plot3d", first, true, kw...)
end
plot3d!(f1::Function, f2::Function, f3::Function, range_t=nothing; kw...) = plot3d(f1, f2, f3, range_t; first=false, kw...)

const plot3  = plot3d			# Alias
const plot3! = plot3d!
# ------------------------------------------------------------------------------------------------------

"""
    scatter(cmd0::String="", arg1=nothing; kwargs...)

Reads (x,y) pairs and plot symbols at those locations on a map.
This module is a subset of ``plot`` to make it simpler to draw scatter plots. So many of
its (fine) controling parameters are not listed here. For a finer control, user should
consult the ``plot`` module.

Parameters
----------

- **G** | **fill** | **markerfacecolor** :: [Type => Str]

    Select color or pattern for filling of symbols or polygons.
    ($(GMTdoc)plot.html#g)
- **N** | **no_clip** :: [Type => Str | []]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)plot.html#n)
- $(GMT.opt_P)
- **S** :: [Type => Str]

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    ($(GMTdoc)plot.html#s)

    Alternatively select a sub-set of symbols using the aliases: **symbol** or **marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
	+ **y**, **y_dash**

    and select their sizes with the **markersize** or **size** keyword [default is 8p].
    The marker size can be a scalar or a vector with same size numeber of rows of data. Units are
    points unless specified otherwise with (for example for cm) *par=(PROJ_LENGTH_UNIT=:c,)*
- **W** | **pen** | **markeredgecolor** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)

[`Full man page`](https://genericmappingtools.github.io/GMT.jl/latest/scatter/)
[`GMT man page`]($(GMTdoc)plot.html)
"""
function scatter(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	common_plot_xyz("", cat_2_arg2(rang, [f(x) for x in rang]), "scatter",  first, false, kw...)
end
scatter!(f::Function, rang=nothing; kw...) = scatter(f, rang; first=false, kw...)

function scatter(f1::Function, f2::Function, range_t=nothing; first=true, kw...)	# Parametric version
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; is3D=false, kw...), "scatter",  first, false, kw...)
end
scatter!(f1::Function, f2::Function, range_t=nothing; kw...) = scatter(f1, f2, range_t; first=false, kw...)

scatter(cmd0::String="", arg1=nothing; first=true, kw...)  = common_plot_xyz(cmd0, arg1, "scatter",  first, false, kw...)
scatter!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "scatter",  false, false, kw...)

scatter(arg; kw...)  = common_plot_xyz("", cat_1_arg(arg), "scatter", true, false, kw...)
scatter!(arg; kw...) = common_plot_xyz("", cat_1_arg(arg), "scatter", false, false, kw...)

scatter(arg1, arg2; kw...)  = common_plot_xyz("", cat_2_arg2(arg1, arg2), "scatter", true, false, kw...)
scatter!(arg1, arg2; kw...) = common_plot_xyz("", cat_2_arg2(arg1, arg2), "scatter", false, false, kw...)

# ------------------------------------------------------------------------------------------------------
scatter3(cmd0::String="", arg1=nothing; kw...)  = common_plot_xyz(cmd0, arg1, "scatter3",  true, true, kw...)
scatter3!(cmd0::String="", arg1=nothing; kw...) = common_plot_xyz(cmd0, arg1, "scatter3",  false, true, kw...)

scatter3(arg; kw...)  = common_plot_xyz("", cat_1_arg(arg), "scatte3", true, false, kw...)
scatter3!(arg; kw...) = common_plot_xyz("", cat_1_arg(arg), "scatte3", false, false, kw...)

function scatter3(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; kw...)
	common_plot_xyz("", hcat(arg1, arg2, arg3), "scatter3", true, true, kw...)
end
function scatter3!(arg1::AbstractArray, arg2::AbstractArray, arg3::AbstractArray; kw...)
	common_plot_xyz("", hcat(arg1, arg2, arg3), "scatter3", false, true, kw...)
end

function scatter3(f1::Function, f2::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_2f(f1, f2, range_t; kw...), "scatter3", first, true, kw...)
end
scatter3!(f1::Function, f2::Function, range_t=nothing; kw...) = scatter3(f1, f2, range_t; first=false, kw...)

function scatter3(f1::Function, f2::Function, f3::Function, range_t=nothing; first=true, kw...)
	common_plot_xyz("", help_parametric_3f(f1, f2, f3, range_t; kw...), "scatter3", first, true, kw...)
end
scatter3!(f1::Function, f2::Function, f3::Function, range_t=nothing; kw...) = scatter3(f1, f2, f3, range_t; first=false, kw...)

const scatter3d  = scatter3			# Alias
const scatter3d! = scatter3!
# ------------------------------------------------------------------------------------------------------

"""
    bar(cmd0::String="", arg1=nothing; kwargs...)

Reads a file or (x,y) pairs and plots vertical bars extending from base to y.

- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **fill** :: [Type => Str --

    Select color or pattern for filling the bars
    ($(GMTdoc)plot.html#g)
- **base** | **bottom** :: [Type => Str | Num]		``key=value``

    By default, base = ymin. Use this option to change that value. If base is not appended then we read it.
    from the last input data column.
- **size** | **width** :: [Type => Str | Num]		``key=value``

    The size or width is the bar width. Append u if size is in x-units. When *width* is used the default is plot-distance units.

Example:

    bar(sort(randn(10)), fill=:black, axis=:auto, show=true)
"""
function bar(cmd0::String="", arg=nothing; first=true, kw...)
	d = KW(kw)
	do_cat = (haskey(d, :stack) || haskey(d, :stacked) && isvector(arg) && length(arg) > 2) ? false : true
	if (cmd0 == "" && do_cat) arg = cat_1_arg(arg)  end	# If ARG is a vector, prepend it with a 1:N x column
	GMT.common_plot_xyz(cmd0, arg, "bar", first, false, kw...)
end
bar!(cmd0::String="", arg=nothing; kw...) = bar(cmd0, arg; first=false, kw...)

function bar(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	bar("", cat_2_arg2(rang, [f(x) for x in rang]); first=first, kw...)
end
bar!(f::Function, rang=nothing; kw...) = bar(f, rang; first=false, kw...)

bar(arg1, arg2; first=true, kw...)  = common_plot_xyz("", cat_2_arg2(arg1, arg2), "bar", first, false, kw...)
bar!(arg1, arg2; kw...) = common_plot_xyz("", cat_2_arg2(arg1, arg2), "bar", false, false, kw...)
bar(arg; kw...)  = bar("", arg; kw...)
bar!(arg; kw...) = common_plot_xyz("", cat_1_arg(arg), "bar", false, false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    bar3(cmd0::String="", arg1=nothing; kwargs...)

Read a grid file, a grid or a MxN matrix and plots vertical bars extending from base to z.

- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **fill** :: [Type => Str]		``key=color``

    Select color or pattern for filling the bars
    ($(GMTdoc)plot.html#g)
- **base** :: [Type => Str | Num]		``key=value``

    By default, base = ymin. Use this option to change that value. If base is not appended then we read it.
- $(GMT.opt_p)

Example:

    G = gmt("grdmath -R-15/15/-15/15 -I0.5 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL =");
    bar3(G, lw=:thinnest, show=true)
"""
function bar3(cmd0::String="", arg=nothing; first=true, kwargs...)
	# Contrary to "bar" this one has specific work to do here.
	d = KW(kwargs)
	opt_z = ""

	if (isa(arg, Array{<:GMTdataset,1}))  arg1 = arg[1]	# It makes no sense accepting > 1 datasets
	else                                  arg1 = arg	# Make a copy that may or not become a new thing
	end

	if (isa(arg1, Array))
		ny, nx = size(arg1)
		if ((nx >= 3 && ny > 3))  arg1 = mat2grid(arg1)  end		# Assume it is a 'bare grid'
	elseif (cmd0 != "")
		if ((val = find_in_dict(d, [:grd :grid])[1]) !== nothing)
			arg1 = gmtread(cmd0, grd=true)
		elseif ((val = find_in_dict(d, [:dataset :table])[1]) !== nothing)
			arg1 = gmtread(cmd0, dataset=true);		arg1 = arg1[1]
		else
			error("BAR3: When first arg is a name, must also state its type. e.g. grd=true or dataset=true")
		end
	end

	opt_base = add_opt(d, "", "", [:base])	# No need to purge because base is not a psxy option

	if (isa(arg1, GMTgrid))
		if (haskey(d, :bar))
			opt_S = parse_bar_cmd(d, :bar, "", "So")
		else
			# 0.85 is the % of inc width of bars
			w1::Float64, w2::Float64 = arg1.inc[1]*0.85, arg1.inc[2]*0.85
			opt_S = " -So$(w1)u/$(w2)u"
			if     (haskey(d, :nbands))  opt_z = string("+z", d[:nbands]);	delete!(d, :nbands)
			elseif (haskey(d, :Nbands))  opt_z = string("+Z", d[:Nbands]);	delete!(d, :Nbands)
			end
		end
		opt, = parse_R(d, "", !first)
		if (opt == "" || opt == " -R")			# OK, no R but we know it here so put it in 'd'
			if (arg1.registration == 1)			# Fine, grid is already pixel reg
				push!(d, :R => arg1.range)
			else								# Need to get a pix reg R
				range = deepcopy(arg1.range)
				range[1] -= arg1.inc[1] / 2;	range[2] += arg1.inc[1] / 2;
				range[3] -= arg1.inc[2] / 2;	range[4] += arg1.inc[2] / 2;
				push!(d, :R => range)
			end
			z_min = arg1.range[5]
		elseif (opt_base == "")					# Shit, need to get zmin out of the opt_R string
			t = split(opt, '/')
			(length(t) == 6) ? z_min = t[5] : error("For 3D cases, region must have 6 selements")
		end
		if (opt_base == "")  push!(d, :base => z_min)  end	# Make base = z_min
		arg1 = gmt("grd2xyz", arg1)[1]			# Now arg1 is a GMTdataset
	else
		opt_S = parse_inc(d, "", [:S :width], "So", true)
		if (opt_S == "")
			opt_S = parse_bar_cmd(d, :bar, "", "So", true)
		end
		if (opt_S == "")
			if ((isa(arg1, Array) && size(arg1,2) < 5) || (isa(arg1, GMTdataset) && size(arg1.data,2) < 5))
				error("BAR3: When NOT providing *width* data must contain at least 5 columns.")
			end
		end
		#if (opt_S != "" && !isletter(opt_S[end]))   opt_S = opt_S * 'u'  end
		if     (haskey(d, :nbands))  opt_z = string("+z", d[:nbands]);	delete!(d, :nbands)
		elseif (haskey(d, :Nbands))  opt_z = string("+Z", d[:Nbands]);	delete!(d, :Nbands)
		end
	end

	opt_base = add_opt(d, "", "", [:base])		# Do this again because :base may have been added above
	if (opt_base == "")
		_z_min::Float32 = (isa(arg1, Array)) ? minimum(view(arg1, :, 3)) : minimum(view(arg1.data, :, 3))
		opt_S *= "+b$_z_min"
	else
		opt_S *= "+b" * opt_base
	end

	common_plot_xyz("", arg1, "bar3|" * opt_S * opt_z, first, true, d...)
end

bar3(arg1; kw...) = bar3("", arg1; first=true, kw...)
bar3!(cmd0::String="", arg1=nothing; kw...) = bar3(cmd0, arg1; first=false, kw...)
bar3!(arg1; kw...) = bar3("", arg1; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    arrows(cmd0::String="", arg1=nothing; arrow=(...), kwargs...)

Plots an arrow field. When the keyword *arrow=(...)* or *vector=(...)* is used, the direction (in degrees
counter-clockwise from horizontal) and length must be found in columns 3 and 4, and size, if not specified
on the command-line, should be present in column 5. The size is the length of the vector head. Vector stem
width is set by option *pen* or *line_attrib*.

The *vecmap=(...)* variation is similar to above except azimuth (in degrees east of north) should be
given instead of direction. The azimuth will be mapped into an angle based on the chosen map projection.
If length is not in plot units but in arbitrary user units (e.g., a rate in mm/yr) then you can use the
*input_col* option to scale the corresponding column via the +sscale modifier.

The *geovec=(...)* or *geovector=(...)* keywords plot geovectors. In geovectors azimuth (in degrees east from north) and geographical length must be found in columns 3 and 4. The size is the length of the vector head. Vector width is set by *pen* or *line_attrib*. Note: Geovector stems are drawn as thin filled polygons and hence pen attributes like dashed and dotted are not available. For allowable geographical units, see the *units=()* option.

The full *arrow* options list can be consulted at [Vector Attributes](@ref)

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)

Example:

	arrows([0 8.2 0 6], limits=(-2,4,0,9), arrow=(len=2,stop=1,shape=0.5,fill=:red), axis=:a, pen="6p", show=true)
"""
function arrows(cmd0::String="", arg1=nothing; first=true, kwargs...)
	# A arrows plotting method of plot

	d = KW(kwargs)
	cmd = helper_arrows(d, true)	# Have to delete to avoid double parsing in -W
	if (cmd == "")  cmd = " -Sv0.5+e+h0.5"	# Minimalist defaults
	else            cmd = " -S" * cmd
	end

	GMT.common_plot_xyz(cmd0, arg1, cmd, first, false, d...)
end

function helper_arrows(d::Dict, del::Bool=true)
	# Helper function to set the vector head attributes
	cmd = ""
	val, symb = find_in_dict(d, [:arrow :vector :arrow4 :vector4 :vecmap :geovec :geovector], del)
	if (val !== nothing)
		code = 'v'
		if (symb == :geovec || symb == :geovector)
			code = '='
		elseif (symb == :vecmap)	# Uses azimuth and plots angles taking projection into account
			code = 'V'
		end
		if (isa(val, String))		# An hard core GMT string directly with options
			if (val[1] != code)    cmd = code * val
			else                   cmd = val		# The GMT string already had vector flag char
			end
		elseif (isa(val, Number))  cmd = code * "$val"
		elseif (symb == :arrow4 || symb == :vector4)  cmd = code * vector4_attrib(val)
		else                       cmd = code * vector_attrib(val)
		end
	end
	return cmd
end
arrows!(cmd0::String="", arg1=nothing; kw...) = arrows(cmd0, arg1; first=false, kw...)
arrows(arg1; kw...)  = arrows("", arg1; first=true, kw...)
arrows!(arg1; kw...) = arrows("", arg1; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    lines(cmd0::String="", arg1=nothing; decorated=(...), kwargs...)

Reads a file or (x,y) pairs and plots a collection of different line with decorations

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for lines or the outline of symbols
    ($(GMTdoc)plot.html#w)

Examples:

    lines([0, 10]; [0, 20], limits=(-2,12,-2,22), proj="M2.5", pen=1, fill=:red,
		  decorated=(dist=(val=1,size=0.25), symbol=:box), show=true)

    lines(x -> cos(x) * x, y -> sin(y) * y, linspace(0,2pi,100), region=(-4,7,-5.5,2.5), lw=2, lc=:sienna,
          decorated=(quoted=true, const_label=" In Vino Veritas  - In Aqua, Rãs & Toads", font=(25,"Times-Italic"),
                     curved=true, pen=(0.5,:red)), aspect=:equal, fmt=:png, show=true)
"""
function lines(cmd0::String="", arg1=nothing; first=true, kwargs...)
	# A lines plotting method of plot
	d = KW(kwargs)
	if ((val = find_in_dict(d, [:decorated])[1]) !== nothing)
		cmd = (isa(val, String)) ? val : decorated(val)
	else
		cmd = "lines"
	end

	common_plot_xyz(cmd0, arg1, cmd, first, false, d...)
end
lines!(cmd0::String="", arg=nothing; kw...) = lines(cmd0, arg; first=false, kw...)

function lines(f::Function, range_x=nothing; first=true, kw...)
	rang = gen_coords4funs(range_x, "x"; kw...)
	lines("", cat_2_arg2(rang, [f(x) for x in rang]); first=first, kw...)
end
lines!(f::Function, rang=nothing; kw...) = lines(f, rang; first=false, kw...)

function lines(f1::Function, f2::Function, range_t=nothing; first=true, kw...)	# Parametric version
	lines("", help_parametric_2f(f1, f2, range_t; is3D=false, kw...); first=first, kw...)
end
lines!(f1::Function, f2::Function, range_t=nothing; kw...) = lines(f1, f2, range_t; first=false, kw...)

lines(arg1, arg2; kw...)  = lines("", cat_2_arg2(arg1, arg2); first=true, kw...)
lines!(arg1, arg2; kw...) = lines("", cat_2_arg2(arg1, arg2); first=false, kw...)
lines(arg; kw...)  = lines("", cat_1_arg(arg); first=true, kw...)
lines!(arg; kw...) = lines("", cat_1_arg(arg); first=false, kw...)
# ------------------------------------------------------------------------------------------------------

"""
    hlines(arg; decorated=(...), kwargs...)

Plots one or a collection of horizontal lines with eventual decorations

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for the horizontal lines
    ($(GMTdoc)plot.html#w)

Example:

    plot(rand(5,3))
    hlines!([0.2, 0.6], pen=(1, :red), show=true)
"""
function hlines(arg1=nothing; first=true, kwargs...)
	# A lines plotting method of plot
	d = KW(kwargs)
	(arg1 === nothing && ((arg1 = find_in_dict(d, [:data])[1]) === nothing)) && error("No input data")
	if ((val = find_in_dict(d, [:decorated])[1]) !== nothing)
		cmd = (isa(val, String)) ? val : decorated(val)
	else
		cmd = "lines"
	end
	mat = ones(2, length(arg1))
	[mat[1,k] = mat[2,k] = arg1[k] for k = 1:length(arg1)]
	if ((opt_R = parse_R(d, "")[2]) != "")  x = vec(opt_R2num(opt_R)[1:2])
	else                                    x = [-1e50, 1e50];
	end
	D = mat2ds(mat, x=x, multi=true)

	common_plot_xyz("", D, cmd, first, false, d...)
end
hlines!(arg=nothing; kw...) = hlines(arg; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------------------
"""
    vlines(arg; decorated=(...), kwargs...)

Plots one or a collection of vertical lines with eventual decorations

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **W** | **pen** | **line_attrib** :: [Type => Str]

    Set pen attributes for the horizontal lines
    ($(GMTdoc)plot.html#w)

Example:

    plot(rand(5,3), region=[0,1,0,1])
    vlines!([0.2, 0.6], pen=(1, :red), show=true)
"""
function vlines(arg1=nothing; first=true, kwargs...)
	# A lines plotting method of plot
	d = KW(kwargs)
	(arg1 === nothing && ((arg1 = find_in_dict(d, [:data])[1]) === nothing)) && error("No input data")
	if ((val = find_in_dict(d, [:decorated])[1]) !== nothing)
		cmd = (isa(val, String)) ? val : decorated(val)
	else
		cmd = "lines"
	end
	mat = ones(2, length(arg1))
	mat[1,:] = mat[2,:] = arg1
	if ((opt_R = parse_R(d, "")[2]) != "")  x = vec(opt_R2num(opt_R)[3:4])
	else                                    x = [-1e50, 1e50];
	end
	D = mat2ds(mat, x=x, multi=true)
	# Now we need tp swapp x / y columns because the vlines case is more complicated to implement.
	for k = 1:length(arg1)
		D[k].data[:,1], D[k].data[:,2] = D[k].data[:,2], D[k].data[:,1]
	end

	common_plot_xyz("", D, cmd, first, false, d...)
end
vlines!(arg=nothing; kw...) = vlines(arg; first=false, kw...)
# ------------------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------------------
function ternary(cmd0::String="", arg1=nothing; first=true, kwargs...)
	# A wrapper for psternary
	common_plot_xyz(cmd0, arg1, "ternary", first, false, kwargs...)
end
ternary!(cmd0::String="", arg1=nothing; kw...) = ternary(cmd0, arg1; first=false, kw...)
ternary(arg1;  kw...)  = ternary("", arg1; first=true, kw...)
ternary!(arg1; kw...)  = ternary("", arg1; first=false, kw...)
const psternary  = ternary            # Aliases
const psternary! = ternary!           # Aliases

"""
    events(cmd0::String, arg1=nothing; kwargs...)

Plot event symbols and labels for a moment in time

Full option list at [`events`]($(GMTdoc)events.html)

Parameters
----------

- **T** | **now** :: [Type => Int | Str]

    Set the current plot time. If absolute times are used you must also use -fT.
    ($(GMTdoc)events.html#t)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **D** | **shift** | **offset** :: [Type => Str]	``Arg = [j|J]dx[/dy][+v[pen]]``

    Offset the text from the projected (x,y) point by dx,dy [0/0].
    ($(GMTdoc)events.html#d)
- **E** | **knots** :: [Type => Str]	``Arg = s|t[+o|Odt][+rdt][+pdt][+ddt][+fdt][+ldt]``

    Set the time knots for the symbol or text time-functions.
    ($(GMTdoc)events.html#e)
- **G** | **fill** :: [Type => Str | Int | Touple]

    Set constant shade or color for all symbols.
    ($(GMTdoc)events.html#g)
- $(GMT.opt_J)
- **L** | **duration** :: [Type => Bool | Number | Str]		``Arg = [length|t]``

    Specify the length (i.e., duration) of the event.
    ($(GMTdoc)events.html#l)
- **M** | **rise** :: [Type => Str]		``Arg = i|s|t[val1][+cval2]``

    Modify the initial intensity, size magnification, or transparency of the symbol during the rise interval.
    ($(GMTdoc)events.html#m)
- **Q** | **save** :: [Type => Number]

    Save the intermediate event symbols and labels to permanent files instead of removing them when done.
    ($(GMTdoc)events.html#q)
- $(GMT.opt_R)
- **W** | **pen** | **markeredgecolor** :: [Type => Str]

    Specify symbol outline pen attributes [Default is no outline].
    ($(GMTdoc)events.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_swap_xy)
"""
# ------------------------------------------------------------------------------------------------------
function events(cmd0::String="", arg1=nothing; kwargs...)
	# events share a lot of options with plot
	d = KW(kwargs)
	cmd = add_opt(d, "", "T", [:T :now])
	if (!occursin("-T", cmd))  error("The 'now' (T) option is mandatory")  end
	cmd = add_opt(d, cmd, "E", [:E :knots],
		(symbol=("s", nothing, 1), text=("t", nothing, 1), shift_startEnd = "+o", shift_start="+O", raise="+r", plateau="+p", decay="+d", fade="+f", text_duration="+l"))
	cmd = add_opt(d, cmd, "M", [:M :rise],
		(intensity=("i", arg2str, 1), size=("s", arg2str, 1), transparency=("t", arg2str, 1), coda="+c"))
	cmd = add_opt(d, cmd, "L", [:L :duration])
	cmd = add_opt(d, cmd, "Q", [:Q :save])
	cmd = add_opt(d, cmd, 'D', [:D :offset],
		(away=("j", nothing, 1), corners=("J", nothing, 1), shift="", line=("+v",add_opt_pen)))
	cmd = add_opt(d, cmd, 'F', [:F :attrib],
		(angle="+a", Angle="+A", font=("+f", font), justify="+j", region_justify="+c", header="_+h", label="_+l", rec_number="+r", text="+t", zvalues="+z"), false)
	common_plot_xyz(cmd0, arg1, "events|" * cmd, true, false, d...)
end
const psevents = events            # Alias

# ------------------------------------------------------------------------------------------------------
function cat_1_arg(arg)
	# Add a first column with 1:n to all args that are not GMTdatasets
	(isa(arg, Array{<:GMTdataset,1}) || isa(arg, GMTdataset))  &&  return arg
	if (isa(arg, Vector) || typeof(arg) <: AbstractRange)
		#arg = hcat(collect(1:size(arg,1)), arg)
		arg = hcat(collect(eltype(arg), 1:size(arg,1)), arg)
	#elseif (isa(arg, Array) && size(arg,1) == 1)		# Accept also row arrays. CAN'T IT BREAKS plot([1 1])
		#arg = hcat(collect(eltype(arg), 1:length(arg)), arg')
	elseif (isa(arg, NTuple))
		arg = hcat(collect(eltype(arg), 1:length(arg)), collect(arg))
	end
	return arg
end

# ------------------------------------------------------------------------------------------------------
function cat_2_arg2(arg1, arg2)
	# Cat two vectors (or tuples) or a vector (or tuple) and a matrix in a Mx2 matrix

	!((isa(arg1, Vector) || typeof(arg1) <: AbstractRange || isa(arg1, NTuple) || isa(arg1, Matrix)) && (isa(arg2, Vector) || typeof(arg2) <: AbstractRange || isa(arg2, NTuple) || isa(arg2, Matrix)) ) &&
		error("Unknown types ($(typeof(arg1))) and ($(typeof(arg2))) in cat_2_arg2() function")

	if (isa(arg1, NTuple))  arg1 = collect(arg1)  end
	if (isa(arg2, NTuple))  arg2 = collect(arg2)  end
	if (size(arg1,1) == 1 && size(arg1,2) != 1)  arg1 = arg1[:]  end
	if (size(arg2,1) == 1 && size(arg2,2) != 1)  arg2 = arg2[:]  end
	arg = hcat(arg1, arg2)
	if (size(arg,2) > 2)  global multi_col[1] = true  end
	return arg
end

# ------------------------------------------------------------------------------------------------------
function gen_coords4funs(rang=nothing, axis="x"; kw...)
	# Generate axes coordenates to use when plot functions

	if (rang === nothing)
		symb = (axis == "x") ? [:xlim] : [:ylim]
		if ((val = find_in_dict(KW(kw), symb)[1]) !== nothing)
		(!isa(val, Tuple) && !isa(val, Array{<:Real})) && error("$(string(symb[1])) must be a tuple or array, not '$(typeof(val))'")
			(length(val) != 2) && error("$(string(symb[1])) must have 2 elements")
			rang = linspace(val[1], val[2], 200)
		else
			rang = linspace(-5, 5)
		end
	elseif (isa(rang, Real))
		rang = linspace(-rang, rang)
	end
	return rang
end

# ------------------------------------------------------------------------------------------------------
# Common code shared by the functions that accept parametric equations
function help_parametric_2f(f1::Function, f2::Function, range_t=nothing; is3D=true, kw...)
	# This function is shared by both the 2D & 3D cases
	t = collect(gen_coords4funs(range_t, "x"; kw...));
	x = [f1(x) for x in t];		y = [f2(y) for y in t]
	out = (is3D) ? hcat(x[:], y[:], t[:]) : hcat(x[:], y[:])
	return out
end

function help_parametric_3f(f1::Function, f2::Function, f3::Function, range_t=nothing; kw...)
	t = collect(gen_coords4funs(range_t, "x"; kw...));
	x = [f1(x) for x in t];		y = [f2(y) for y in t];		z = [f3(z) for z in t]
	return hcat(x[:], y[:], z[:])
end
# ------------------------------------------------------------------------------------------------------


text:
@doc raw"""
    HessianProblem <: Problem

specify a problem for hessian based algorithms.

# Fields
* `M`            : a manifold $\mathcal M$
* `cost` : a function $F:\mathcal M→ℝ$ to minimize
* `gradient`     : the gradient $\operatorname{grad}F:\mathcal M
  → \mathcal T\mathcal M$ of the cost function $F$
* `hessian`      : the hessian $\operatorname{Hess}F(x)[⋅]: \mathcal T_{x} \mathcal M
  → \mathcal T_{x} \mathcal M$ of the cost function $F$
* `precon`       : the symmetric, positive definite
    preconditioner (approximation of the inverse of the Hessian of $F$)

# See also
[`truncated_conjugate_gradient_descent`](@ref), [`trust_regions`](@ref)
"""
struct HessianProblem{T,mT<:AbstractManifold,C,G,H,Pre} <: AbstractGradientProblem{T}
    M::mT
    cost::C
    gradient!!::G
    hessian!!::H
    precon::Pre
    function HessianProblem(
        M::mT,
        cost::C,
        grad::G,
        hess::H,
        pre::Pre;
        evaluation::AbstractEvaluationType=AllocatingEvaluation(),
    ) where {mT<:AbstractManifold,C,G,H,Pre}
        return new{typeof(evaluation),mT,C,G,H,Pre}(M, cost, grad, hess, pre)
    end
end

@doc raw"""
    AbstractHessianOptions <: Options

An [`Options`](@ref) type to represent algorithms that employ the Hessian.
These options are assumed to have a field (`gradient`) to store the current gradient ``\operatorname{grad}f(x)``
"""
abstract type AbstractHessianOptions <: AbstractGradientOptions end

@doc raw"""
    TruncatedConjugateGradientOptions <: AbstractHessianOptions

describe the Steihaug-Toint truncated conjugate-gradient method, with

# Fields
a default value is given in brackets if a parameter can be left out in initialization.

* `x` : a point, where the trust-region subproblem needs
    to be solved
* `stop` : a function s,r = @(o,iter,ξ,x,xnew) returning a stop
    indicator and a reason based on an iteration number, the gradient and the
    last and current iterates
* `gradient` : the gradient at the current iterate
* `η` : a tangent vector (called update vector), which solves the
    trust-region subproblem after successful calculation by the algorithm
* `δ` : search direction
* `trust_region_radius` : the trust-region radius
* `residual` : the gradient
* `randomize` : indicates if the trust-region solve and so the algorithm is to be
        initiated with a random tangent vector. If set to true, no
        preconditioner will be used. This option is set to true in some
        scenarios to escape saddle points, but is otherwise seldom activated.
* `project_vector!` : (`copyto!`) specify a projection operation for tangent vectors
    for numerical stability. A function `(M, Y, p, X) -> ...` working in place of `Y`.
    per default, no projection is perfomed, set it to `project!` to activate projection.

# Constructor

    TruncatedConjugateGradientOptions(x, stop, eta, delta, Delta, res, uR)

construct a truncated conjugate-gradient Option with the fields as above.

# See also
[`truncated_conjugate_gradient_descent`](@ref), [`trust_regions`](@ref)
"""
mutable struct TruncatedConjugateGradientOptions{P,T,R<:Real,SC<:StoppingCriterion,Proj} <:
               AbstractHessianOptions
    x::P
    stop::SC
    gradient::T
    η::T
    Hη::T
    δ::T
    Hδ::T
    δHδ::R
    ηPδ::R
    δPδ::R
    ηPη::R
    z::T
    z_r::R
    residual::T
    trust_region_radius::R
    model_value::R
    new_model_value::R
    κ::R
    randomize::Bool
    project_vector!::Proj
    initialResidualNorm::Float64
    function TruncatedConjugateGradientOptions(
        p::HessianProblem,
        x::P,
        η::T,
        trust_region_radius::R,
        randomize::Bool;
        project_vector!::Proj=copyto!,
        θ::Float64=1.0,
        κ::Float64=0.1,
        stop::StoppingCriterion=StopWhenAny(
            StopAfterIteration(manifold_dimension(p.M)),
            StopWhenAll(
                StopIfResidualIsReducedByPower(θ), StopIfResidualIsReducedByFactor(κ)
            ),
            StopWhenTrustRegionIsExceeded(),
            StopWhenCurvatureIsNegative(),
            StopWhenModelIncreased(),
        ),
    ) where {H,G,P,T,R<:Real,Proj}
        o = new{
            typeof(x),
            typeof(η),
            typeof(trust_region_radius),
            typeof(stop),
            typeof(project_vector!),
        }()
        o.x = x
        o.stop = stop
        o.η = η
        o.trust_region_radius = trust_region_radius
        o.randomize = randomize
        o.project_vector! = project_vector!
        o.model_value = zero(trust_region_radius)
        o.κ = zero(trust_region_radius)
        return o
    end
end

@doc raw"""
    TrustRegionsOptions <: AbstractHessianOptions

describe the trust-regions solver, with

# Fields
a default value is given in brackets if a parameter can be left out in initialization.

* `x` : a point as starting point
* `stop` : a function s,r = @(o,iter) returning a stop
    indicator and a reason based on an iteration number and the gradient
* `trust_region_radius` : the (initial) trust-region radius
* `max_trust_region_radius` : the maximum trust-region radius
* `randomize` : indicates if the trust-region solve is to be initiated with a
        random tangent vector. If set to true, no preconditioner will be
        used. This option is set to true in some scenarios to escape saddle
        points, but is otherwise seldom activated.
* `project_vector!` : (`copyto!`) specify a projection operation for tangent vectors
    for numerical stability. A function `(M, Y, p, X) -> ...` working in place of `Y`.
    per default, no projection is perfomed, set it to `project!` to activate projection.
* `ρ_prime` : a lower bound of the performance ratio for the iterate that
        decides if the iteration will be accepted or not. If not, the
        trust-region radius will have been decreased. To ensure this,
        ρ'>= 0 must be strictly smaller than 1/4. If ρ' is negative,
        the algorithm is not guaranteed to produce monotonically decreasing
        cost values. It is strongly recommended to set ρ' > 0, to aid
        convergence.
* `ρ_regularization` : Close to convergence, evaluating the performance ratio ρ
        is numerically challenging. Meanwhile, close to convergence, the
        quadratic model should be a good fit and the steps should be
        accepted. Regularization lets ρ go to 1 as the model decrease and
        the actual decrease go to zero. Set this option to zero to disable
        regularization (not recommended). When this is not zero, it may happen
        that the iterates produced are not monotonically improving the cost
        when very close to convergence. This is because the corrected cost
        improvement could change sign if it is negative but very small.

# Constructor

    TrustRegionsOptions(x, stop, delta, delta_bar, uR, rho_prime, rho_reg)

construct a trust-regions Option with the fields as above.

# See also
[`trust_regions`](@ref)
"""
mutable struct TrustRegionsOptions{
    P,T,SC<:StoppingCriterion,RTR<:AbstractRetractionMethod,R<:Real,Proj
} <: AbstractHessianOptions
    x::P
    gradient::T
    stop::SC
    trust_region_radius::R
    max_trust_region_radius::R
    retraction_method::RTR
    randomize::Bool
    project_vector!::Proj
    ρ_prime::R
    ρ_regularization::R

    tcg_options::TruncatedConjugateGradientOptions{P,T,R}

    x_proposal::P
    f_proposal::R
    timer::Float64
    # Random
    Hgrad::T
    η::T
    Hη::T
    η_Cauchy::T
    Hη_Cauchy::T
    τ::R
    function TrustRegionsOptions{P,T,SC,RTR,R,Proj}(
        x::P,
        grad::T,
        trust_region_radius::R,
        max_trust_region_radius::R,
        ρ_prime::R,
        ρ_regularization::R,
        randomize::Bool,
        stopping_citerion::SC,
        retraction_method::RTR;
        project_vector!::Proj=copyto!,
    ) where {P,T,SC<:StoppingCriterion,RTR<:AbstractRetractionMethod,R<:Real,Proj}
        o = new{P,T,SC,RTR,R,Proj}()
        o.x = x
        o.gradient = grad
        o.stop = stopping_citerion
        o.retraction_method = retraction_method
        o.trust_region_radius = trust_region_radius
        o.max_trust_region_radius = max_trust_region_radius::R
        o.ρ_prime = ρ_prime
        o.ρ_regularization = ρ_regularization
        o.randomize = randomize
        o.project_vector! = project_vector!
        o.timer = 0.0
        return o
    end
end
function TrustRegionsOptions(
    x::P,
    grad::T,
    trust_region_radius::R,
    max_trust_region_radius::R,
    ρ_prime::R,
    ρ_regularization::R,
    randomize::Bool,
    stopping_citerion::SC;
    retraction_method::RTR=ExponentialRetraction(),
    project_vector!::Proj=copyto!,
) where {P,T,R<:Real,SC<:StoppingCriterion,RTR<:AbstractRetractionMethod,Proj}
    return TrustRegionsOptions{P,T,SC,RTR,R,Proj}(
        x,
        grad,
        trust_region_radius,
        max_trust_region_radius,
        ρ_prime,
        ρ_regularization,
        randomize,
        stopping_citerion,
        retraction_method;
        (project_vector!)=project_vector!,
    )
end

@doc raw"""
    get_hessian(p::HessianProblem{T}, q, X)
    get_hessian!(p::HessianProblem{T}, Y, q, X)

evaluate the Hessian of a [`HessianProblem`](@ref) `p` at the point `q`
applied to a tangent vector `X`, i.e. ``\operatorname{Hess}f(q)[X]``.

The evaluation is done in place of `Y` for the `!`-variant.
The `T=`[`AllocatingEvaluation`](@ref) problem might still allocate memory within.
When the non-mutating variant is called with a `T=`[`MutatingEvaluation`](@ref)
memory for the result is allocated.
"""
function get_hessian(p::HessianProblem{AllocatingEvaluation}, q, X)
    return p.hessian!!(p.M, q, X)
end
function get_hessian(p::HessianProblem{MutatingEvaluation}, q, X)
    Y = zero_vector(p.M, q)
    return p.hessian!!(p.M, Y, q, X)
end
function get_hessian!(p::HessianProblem{AllocatingEvaluation}, Y, q, X)
    return copyto!(p.M, Y, p.hessian!!(p.M, q, X))
end
function get_hessian!(p::HessianProblem{MutatingEvaluation}, Y, q, X)
    return p.hessian!!(p.M, Y, q, X)
end

@doc raw"""
    get_preconditioner(p,x,ξ)

evaluate the symmetric, positive definite preconditioner (approximation of the
inverse of the Hessian of the cost function `F`) of a
[`HessianProblem`](@ref) `p` at the point `x`applied to a
tangent vector `ξ`.
"""
get_preconditioner(p::HessianProblem, x, X) = p.precon(p.M, x, X)

@doc raw"""
    ApproxHessianFiniteDifference{T, mT, P, G}

A functor to approximate the Hessian by a finite difference of gradient evaluations

# Constructor

    ApproxHessianFiniteDifference(M, x, gradF)

Initialize the approximate hessian to combute ``\operatorname{Hess}F`` based on the gradient
`gradF` of a function ``F`` on ``\mathcal M``.
The value `x` is used to initialize a few internal fields.

## Optional Keyword arguments

* `steplength` - (`2^-14`) default step size for the approximation
* `evaluation` - ([`AllocatingEvaluation`](@ref)`()`) specify whether the gradient is allocating or mutating.
* `retraction_method` – (`default_retraction_method(M)`) a `retraction(M, p, X)` to use in the approximation.
* `vector_transport_method` - (`default_vector_transport_method(M)`) a vector transport to use
"""
mutable struct ApproxHessianFiniteDifference{E,P,T,G,RTR,VTR,R<:Real}
    x_dir::P
    gradient!!::G
    grad_tmp::T
    grad_tmp_dir::T
    retraction_method::RTR
    vector_transport_method::VTR
    steplength::R
end
function ApproxHessianFiniteDifference(
    M::mT,
    x::P,
    grad::G;
    steplength::R=2^-14,
    evaluation=AllocatingEvaluation(),
    retraction_method::RTR=default_retraction_method(M),
    vector_transport_method::VTR=default_vector_transport_method(M),
) where {
    mT<:AbstractManifold,
    P,
    G,
    R<:Real,
    RTR<:AbstractRetractionMethod,
    VTR<:AbstractVectorTransportMethod,
}
    X = zero_vector(M, x)
    Y = zero_vector(M, x)
    return ApproxHessianFiniteDifference{typeof(evaluation),P,typeof(X),G,RTR,VTR,R}(
        x, grad, X, Y, retraction_method, vector_transport_method, steplength
    )
end
function (f::ApproxHessianFiniteDifference{AllocatingEvaluation})(M, x, X)
    norm_X = norm(M, x, X)
    (norm_X ≈ zero(norm_X)) && return zero_vector(M, x)
    c = f.steplength / norm_X
    f.grad_tmp .= f.gradient!!(M, x)
    f.x_dir .= retract(M, x, c * X, f.retraction_method)
    f.grad_tmp_dir .= f.gradient!!(M, f.x_dir)
    f.grad_tmp_dir .= vector_transport_to(
        M, f.x_dir, f.grad_tmp_dir, x, f.vector_transport_method
    )
    return (1 / c) * (f.grad_tmp_dir - f.grad_tmp)
end
function (f::ApproxHessianFiniteDifference{MutatingEvaluation})(M, Y, x, X)
    norm_X = norm(M, x, X)
    (norm_X ≈ zero(norm_X)) && return zero_vector!(M, X, x)
    c = f.steplength / norm_X
    f.gradient!!(M, f.grad_tmp, x)
    retract!(M, f.x_dir, x, c * X, f.retraction_method)
    f.gradient!!(M, f.grad_tmp_dir, f.x_dir)
    vector_transport_to!(
        M, f.grad_tmp_dir, f.x_dir, f.grad_tmp_dir, x, f.vector_transport_method
    )
    Y .= (1 / c) .* (f.grad_tmp_dir .- f.grad_tmp)
    return Y
end

@doc raw"""
    StopIfResidualIsReducedByFactor <: StoppingCriterion

A functor for testing if the norm of residual at the current iterate is reduced
by a factor compared to the norm of the initial residual, i.e.
$\Vert r_k \Vert_x \leqq κ \Vert r_0 \Vert_x$.
In this case the algorithm reached linear convergence.

# Fields
* `κ` – the reduction factor
* `initialResidualNorm` - stores the norm of the residual at the initial vector
    ``η`` of the Steihaug-Toint tcg mehtod [`truncated_conjugate_gradient_descent`](@ref)
* `reason` – stores a reason of stopping if the stopping criterion has one be
  reached, see [`get_reason`](@ref).

# Constructor

    StopIfResidualIsReducedByFactor(iRN, κ)

initialize the StopIfResidualIsReducedByFactor functor to indicate to stop after
the norm of the current residual is lesser than the norm of the initial residual
iRN times κ.

# See also
[`truncated_conjugate_gradient_descent`](@ref), [`trust_regions`](@ref)
"""
mutable struct StopIfResidualIsReducedByFactor <: StoppingCriterion
    κ::Float64
    reason::String
    StopIfResidualIsReducedByFactor(κ::Float64) = new(κ, "")
end
function (c::StopIfResidualIsReducedByFactor)(
    p::P, o::O, i::Int
) where {P<:HessianProblem,O<:TruncatedConjugateGradientOptions}
    if norm(p.M, o.x, o.residual) <= o.initialResidualNorm * c.κ && i > 0
        c.reason = "The algorithm reached linear convergence (residual at least reduced by κ=$(c.κ)).\n"
        return true
    end
    return false
end

@doc raw"""
    StopIfResidualIsReducedByPower <: StoppingCriterion

A functor for testing if the norm of residual at the current iterate is reduced
by a power of 1+θ compared to the norm of the initial residual, i.e.
$\Vert r_k \Vert_x \leqq  \Vert r_0 \Vert_{x}^{1+\theta}$. In this case the
algorithm reached superlinear convergence.

# Fields
* `θ` – part of the reduction power
* `initialResidualNorm` - stores the norm of the residual at the initial vector
    $η$ of the Steihaug-Toint tcg mehtod [`truncated_conjugate_gradient_descent`](@ref)
* `reason` – stores a reason of stopping if the stopping criterion has one be
    reached, see [`get_reason`](@ref).

# Constructor

    StopIfResidualIsReducedByPower(iRN, θ)

initialize the StopIfResidualIsReducedByFactor functor to indicate to stop after
the norm of the current residual is lesser than the norm of the initial residual
iRN to the power of 1+θ.

# See also
[`truncated_conjugate_gradient_descent`](@ref), [`trust_regions`](@ref)
"""
mutable struct StopIfResidualIsReducedByPower <: StoppingCriterion
    θ::Float64
    reason::String
    StopIfResidualIsReducedByPower(θ::Float64) = new(θ, "")
end
function (c::StopIfResidualIsReducedByPower)(
    p::P, o::O, i::Int
) where {P<:HessianProblem,O<:TruncatedConjugateGradientOptions}
    if norm(p.M, o.x, o.residual) <= o.initialResidualNorm^(1 + c.θ) && i > 0
        c.reason = "The algorithm reached superlinear convergence (residual at least reduced by power 1 + θ=$(1+(c.θ))).\n"
        return true
    end
    return false
end

@doc raw"""
    StopWhenTrustRegionIsExceeded <: StoppingCriterion

A functor for testing if the norm of the next iterate in the  Steihaug-Toint tcg
mehtod is larger than the trust-region radius, i.e. $\Vert η_{k}^{*} \Vert_x
≧ trust_region_radius$. terminate the algorithm when the trust region has been left.

# Fields
* `reason` – stores a reason of stopping if the stopping criterion has one be
    reached, see [`get_reason`](@ref).
* `storage` – stores the necessary parameters `η, δ, residual` to check the
    criterion.

# Constructor

    StopWhenTrustRegionIsExceeded([a])

initialize the StopWhenTrustRegionIsExceeded functor to indicate to stop after
the norm of the next iterate is greater than the trust-region radius using the
[`StoreOptionsAction`](@ref) `a`, which is initialized to store
`:η, :δ, :residual` by default.

# See also
[`truncated_conjugate_gradient_descent`](@ref), [`trust_regions`](@ref)
"""
mutable struct StopWhenTrustRegionIsExceeded <: StoppingCriterion
    reason::String
end
StopWhenTrustRegionIsExceeded() = StopWhenTrustRegionIsExceeded("")
function (c::StopWhenTrustRegionIsExceeded)(
    ::HessianProblem, o::TruncatedConjugateGradientOptions, i::Int
)
    if o.ηPη >= o.trust_region_radius^2 && i >= 0
        c.reason = "Trust-region radius violation (‖η‖² = $(o.ηPη)) >= $(o.trust_region_radius^2) = trust_region_radius²). \n"
        return true
    end
    return false
end

@doc raw"""
    StopWhenCurvatureIsNegative <: StoppingCriterion

A functor for testing if the curvature of the model is negative, i.e.
$\langle \delta_k, \operatorname{Hess}[F](\delta_k)\rangle_x \leqq 0$.
In this case, the model is not strictly convex, and the stepsize as computed
does not give a reduction of the model.

# Fields
* `reason` – stores a reason of stopping if the stopping criterion has one be
    reached, see [`get_reason`](@ref).

# Constructor

    StopWhenCurvatureIsNegative()

# See also
[`truncated_conjugate_gradient_descent`](@ref), [`trust_regions`](@ref)
"""
mutable struct StopWhenCurvatureIsNegative <: StoppingCriterion
    reason::String
end
StopWhenCurvatureIsNegative() = StopWhenCurvatureIsNegative("")
function (c::StopWhenCurvatureIsNegative)(
    ::HessianProblem, o::TruncatedConjugateGradientOptions, i::Int
)
    if o.δHδ <= 0 && i > 0
        c.reason = "Negative curvature. The model is not strictly convex (⟨δ,Hδ⟩_x = $(o.δHδ))) <= 0).\n"
        return true
    end
    return false
end

@doc raw"""
    StopWhenModelIncreased <: StoppingCriterion

A functor for testing if the curvature of the model value increased.

# Fields
* `reason` – stores a reason of stopping if the stopping criterion has one be
    reached, see [`get_reason`](@ref).

# Constructor

    StopWhenModelIncreased()

# See also
[`truncated_conjugate_gradient_descent`](@ref), [`trust_regions`](@ref)
"""
mutable struct StopWhenModelIncreased <: StoppingCriterion
    reason::String
end
StopWhenModelIncreased() = StopWhenModelIncreased("")
function (c::StopWhenModelIncreased)(
    ::HessianProblem, o::TruncatedConjugateGradientOptions, i::Int
)
    if i > 0 && (o.new_model_value > o.model_value)
        c.reason = "Model value increased from $(o.model_value) to $(o.new_model_value).\n"
        return true
    end
    return false
end

function (d::DebugGradientNorm)(p::HessianProblem, o::TrustRegionsOptions, i::Int)
    (i < 1) && return nothing
    Printf.format(d.io, Printf.Format(d.format), norm(p.M, o.x, o.gradient))
    return nothing
end

function (r::RecordGradientNorm)(
    p::P, o::O, i::Int
) where {P<:HessianProblem,O<:TrustRegionsOptions}
    return record_or_reset!(r, norm(p.M, o.x, o.gradient), i)
end


text:
#!/usr/sbin/julia
import Primes: nextprime

struct Data
  g::BigInt
  q::BigInt
  sks::Array{Array{BigInt}}
  pks::Array{BigInt}
  secret::Array{UInt8}
end

function create_key(len::Int, data::Data)
  sk::Array{BigInt} = []
  pk::BigInt = 1

  if length(data.pks) ≥ 10
    println("Enough of key creation!")
    return
  end

  for i in 1:len
    push!(sk, rand(1:data.q))
    pk = powermod(data.g, sk[i], data.q) * pk % data.q
  end
  push!(data.sks, sk)
  push!(data.pks, pk)
  println("Key created sucessfully!")
end

function create_secret(len::Int, data::Data)
  println("Tell me your secret: ")
  readbytes!(stdin, data.secret, len)
end

function show_data(data::Data)
  println("g = ", data.g)
  println("q = ", data.q)
  println("secret = ", data.secret)
  println("pks = ", data.pks)
end

function encrypt(flag::String, data::Data)
  i::Int = 0
  c::BigInt = 1
  d::BigInt = 1
  r::BigInt = 1
  s::Int = 0

  if length(data.pks) < 4
    println("Error: You don't have enough keys to encrypt.")
    return
  end

  if length(data.secret) ≠ length(flag)
    println("Error: Invalid secret. Make sure you create a secret first!")
    return
  end

  println("Which char index you want to encrypt?")
  try i = parse(Int, readline()) catch e i = -1 end
  if !(1 ≤ i ≤ length(flag))
    println("Error: Invalid char index!")
    return
  end

  for ki in 1:2:(length(data.pks)-1)
    c = powermod(data.pks[ki], data.sks[ki + 1][i], data.q) * c % data.q
    d = powermod(data.pks[ki + 1], data.sks[ki][i], data.q) * d % data.q
  end

  r = rand(1:data.q)
  s = data.secret[i] ⊻ Int(flag[i])
  c = powermod(data.g, s + r, data.q) * c % data.q
  d = powermod(data.g, r, data.q) * d % data.q

  println("Here we go: ", (c, d))
end

function menu()::Int
  println()
  println("1- Create Key")
  println("2- Create Secret")
  println("3- Show data")
  println("4- Encrypt Flag")
  println("5- Exit")
  println()
  option::Int = 0
  try option = parse(Int, readline()) catch e return -1 end
  return option
end

function main()
  flag::String = strip(open(f->read(f, String), "flag.txt"))
  data::Data = Data(2, nextprime(rand(0: big(2)^1024)), [], [], [])
  @assert length(flag) == 40

  println("\nWelcome! What's your plan for today?")
  while true
    option::Int = menu()
    if option == 1
      create_key(length(flag), data)
    elseif option == 2
      create_secret(length(flag), data)
    elseif option == 3
      show_data(data)
    elseif option == 4
      encrypt(flag, data)
    elseif option == 5
      println("Bye!")
      return
    else
      println("Invalid option!")
    end
  end
end

main()

max_stars_repo_name: ruhan-islam/ctf-archives

max_issues_repo_name: ruhan-islam/ctf-archives

max_forks_repo_name: ruhan-islam/ctf-archives


text:
using Plots
using Statistics

function Leveling(surface)
    maxlen = 0
    surf = surface
    for i in surf
        if length(i) > maxlen
            maxlen = length(i)
        end
    end
    for i in surf
        append!(i, [3 for i=1:maxlen - length(i)])
    end
    return surf
end

function Deposition_2DG(len, tot_time, rate, color_step)
    surf = [[] for i=1:len + 2]
    for n in 0:tot_time
        randsurf = rand(2:len + 1, rate)
        for i in randsurf
            if i == 2 && length(surf[i]) >= length(surf[i-1]) - 2
                append!(surf[i-1], n%color_step)
            elseif i == len + 1 && length(surf[i]) >= length(surf[i+1]) - 2
                append!(surf[i+1], n%color_step)
            end
            lens_ = Dict(
                length(surf[i-1])=>i-1,
                length(surf[i])=>i,
                length(surf[i+1])=>i+1,
                    )
            minlen = min(length(surf[i-1]),length(surf[i]),length(surf[i+1]))
            append!(surf[lens_[minlen]], n%color_step)
        end
    end
    return Leveling(surf)
end

# theme(:dark)
# plotlyjs()

surf = Deposition_2DG(200, 100, 1000, 10)

heatmap(hcat(surf...), c = :solar, legend = false)
savefig("../../Figs/Q1/Graphic.pdf")


text:
include("fnUtil.jl")
include("clsUtil.jl")

#= ClsMethods = Dict{Symbol, Dict{Expr, Expr}}() =#
ClsMethods = Dict(:Any=>Dict{Expr, Expr}())
ClsFields = Dict(:Any=>Vector{Expr}())

macro class(ClsName, Cbody)
    ClsName, ParentClsNameLst = getCAndP(ClsName)
    AbsClsName = getAbstractCls(ClsName)
    AbsParentClsName = getAbstractCls(ParentClsNameLst)

    ClsFields[ClsName] = fields = copyFields(ParentClsNameLst, ClsFields)
    ClsMethods[ClsName] = methods = Dict{Expr,Expr}()

    cons = Any[]
    hasInit = false

    # record fields and methods separately
    for (i, block) in enumerate(Cbody.args)
        if isa(block, Symbol)
            union!(fields, [:($block::Any)])
        elseif isa(block, LineNumberNode)
            continue
        elseif block.head == :(::)
            union!(fields, [block])
        elseif block.head == :line
            continue
        elseif block.head == :(=) || block.head == :function
            fname = getFnName(block, withoutGeneric=true)
            if fname == ClsName
                append!(cons, [block])
            elseif fname == :__init__
                hasInit = true
                setFnName!(block, ClsName)
                self = findFnSelfArgNameSymbol(block, ClsName)
                deleteFnSelf!(block)
                prepend!(block.args[2].args, [:($self = $ClsName(()))])
                append!(block.args[2].args, [:($self)])
                append!(cons, [block])
            else
                fn = copy(block)
                setFnSelfArgType!(fn, ClsName)
                methods[findFnCall(fn)] = fn
            end
        else
            error("@class: Case not handled")
        end
    end

    ClsFnCalls = Set(keys(methods))
    for parent in ParentClsNameLst
        for pfn in values(ClsMethods[parent])
            fn = copy(pfn)
            setFnSelfArgType!(fn, ClsName)
            fnCall = findFnCall(fn)
            if haskey(methods, fnCall)
                fName = getFnName(fn, withoutGeneric=true)
                if !(fnCall in ClsFnCalls)
                    error("Ambiguious Function Definition: Multiple definition of function $fName while $ClsName does not overwtie this function!!")
                end
                setFnName!(fn, Symbol(string("super_", parent, fName)), withoutGeneric=true)
                methods[fnCall] = fn
            else
                methods[fnCall] = fn
            end
        end
    end

    cons_str = join(cons,"\n") * "\n"
    if hasInit
        cons_str *= "$ClsName(::Tuple{}) = new()\n"
    end

    clsDefStr = """
              mutable struct $ClsName
                  $(join(fields,"\n"))
              """ * cons_str * """
              end
              """
    # Escape here because we want ClsName and the methods be defined in user scope instead of OOPMacro module scope.
    esc(Expr(:block, Meta.parse(clsDefStr), values(methods)...))
end

macro super(ParentClsName, FCall)
    fname = getFnName(FCall, withoutGeneric=true)
    setFnName!(FCall, Symbol(string("super_", ParentClsName, fname)), withoutGeneric=true)
    esc(FCall)
end


text:
"""A `Bank2D` is a two-dimensional wavelet filter bank, parametrized by
* `T`: numeric type of input, e.g. Float32, Float64.
* `D`: transform domain. Either `FourierDomain{2}` or `SpatialDomain{2}`.
* `G`: point group. Either `TrivialGroup` or `RotationGroup`.
* `W`: wavelet class, e.g. `Morlet` or `MexicanHat`.
Its fields are
* `ϕ`: low-pass filter, also called scaling function.
* `ψs`: 3d array of wavelets, indexed by spin, chroma, and octave.
* `behavior`: mutable behavior, e.g. oversampling.
* `spec`: immutable specifications, e.g. number of filters per octave.
To create a `Bank2D`
1. define a `Spec2D`,
2. if needed, provide behavior characteristics as keyword arguments.
Example:
spec = Spec2D(n_orientations = 8)
bank = Bank2D(spec)"""
immutable Bank2D{
        T<:Number,
        D<:PlaneDomains,
        G<:PlaneGroups,
        W<:RedundantWaveletClass} <: AbstractBank{T,D,G,W}
    ϕ::AbstractFilter{T,D}
    ψs::Array{AbstractFilter{T,D},3}
    behavior::Behavior{T}
    spec::Spec2D{T,D,G,W}

    function (::Type{Bank2D}){T,D,G,W}(
            spec::Spec2D{T,D,G,W},
            pathkey::PathKey = PathKey(:space) ;
            is_ϕ_applied::Bool = false,
            j_range::UnitRange{Int} = 0:(spec.n_octaves-1),
            log2_oversampling::Int = 0,
            max_log2_stride::Int = spec.n_octaves-1,
            weighting::AbstractWeighting = EqualWeighting())
        (nΘs, nΧs, nJs) = size(spec.ψmetas)
        ψs = Array(AbstractFilter{T,D}, (nΘs, nΧs, nJs))
        # TODO: write AbstractFilter{T}(y::Vector{T}, spec::AbstractSpec{T,FourierDomain{1}})
        #ψs[1, :, :] =
        #    pmap(AbstractFilter, spec.ψmetas[1, :, :], fill(spec, nΧs, nJs))
        #(nΘs > 1) && spin!(ψs)
        #ϕ = AbstractFilter(spec.ϕmeta, spec)
        #renormalize!(ϕ, ψs, spec)
        #behavior = Behavior(ϕ, ψs, spec, is_ϕ_applied, j_range,
        #    log2_oversampling, max_log2_stride, pathkey, weighting)
        #new{T,D,G,W}(ϕ, ψs, behavior, spec)
    end
end

Base.ndims(::Bank2D) = 2


text:
function bundle_adjustment!(
    cache, camera::Camera; iterations::Int = 10,
    show_trace::Bool = false, repr_ϵ::Real = 5.0,
)
    fx, fy, cx, cy = camera.fx, camera.fy, camera.cx, camera.cy
    n_observations = length(cache.observations)
    n_poses, n_points = length(cache.poses_remap), length(cache.points_remap)
    poses_shift = n_poses * 6

    ignore_outliers = false
    Y = zeros(Float64, n_observations * 2)

    function residue!(Y, X)
        poses = reshape(@view(X[1:poses_shift]), 6, n_poses)
        points = reshape(@view(X[(poses_shift + 1):end]), 3, n_points)

        @simd for i in 1:n_observations
            id = (i - 1) * 2
            if ignore_outliers && cache.outliers[i]
                Y[id + 1] = 0.0
                Y[id + 2] = 0.0
            else
                pt = @view(points[:, cache.points_ids[i]])
                T = @view(poses[:, cache.poses_ids[i]])
                pt = RotZYX(T[1:3]...) * pt .+ T[4:6] # (x, y, z) format.
                px = @view(cache.pixels[:, i]) # (y, x) format.

                inv_z = 1.0 / pt[3]
                Y[id + 1] = px[1] - (fy * pt[2] * inv_z + cy)
                Y[id + 2] = px[2] - (fx * pt[1] * inv_z + cx)
            end
        end
    end

    optimizer = LevenbergMarquardt(LeastSquaresOptim.LSMR())

    # Fast run, to detect outliers.
    θ0 = cache.θ
    sparsity, g! = _get_jacobian_sparsity(cache, residue!, ignore_outliers)

    Y1 = optimize!(
        LeastSquaresProblem(θ0, Y, residue!, sparsity, g!), optimizer;
        iterations=5, show_trace)
    θ1 = Y1.minimizer
    n_outliers = _ba_detect_outliers!(cache, θ1, camera; repr_ϵ)
    @debug "BA N Outliers $n_outliers."

    ignore_outliers = true
    sparsity2, g2! = _get_jacobian_sparsity(cache, residue!, ignore_outliers)

    Y2 = optimize!(
        LeastSquaresProblem(θ1, Y, residue!, sparsity2, g2!), optimizer;
        iterations, show_trace)
    copy!(cache.θ, Y2.minimizer)
end

function _get_jacobian_sparsity(cache, residue_f, ignore_outliers)
    n_observations = length(cache.observations)
    n_poses = length(cache.poses_remap)
    n_points = length(cache.points_remap)
    poses_shift = n_poses * 6
    n_parameters = poses_shift + n_points * 3

    sparsity = spzeros(Float64, n_observations * 2, n_parameters)
    @simd for i in 1:n_observations
        # If mappoint is outlier for that observation,
        # then leave zero Jacobian for it, e.g. no updates.
        if !(ignore_outliers && cache.outliers[i])
            id = 2 * (i - 1)

            # Set Jacobians for point's coordinates.
            pid = poses_shift + (cache.points_ids[i] - 1) * 3
            sparsity[(id+1):(id+2), (pid+1):(pid+3)] .= 1.0

            # Set Jacobians for poses if they are not constant.
            eid = cache.poses_ids[i]
            if !cache.θconst[eid]
                eid = (eid - 1) * 6
                sparsity[(id+1):(id+2), (eid+1):(eid+6)] .= 1.0
            end
        end
    end

    colorvec = matrix_colors(sparsity)
    grad! = (j, x) -> forwarddiff_color_jacobian!(
        j, residue_f, x; colorvec, sparsity)
    sparsity, grad!
end

function _ba_detect_outliers!(cache, θ, camera::Camera; repr_ϵ, depth_ϵ = 1e-6)
    n_observations = length(cache.observations)
    n_poses = length(cache.poses_remap)
    n_points = length(cache.points_remap)
    poses_shift = n_poses * 6

    poses = reshape(@view(θ[1:poses_shift]), 6, n_poses)
    points = reshape(@view(θ[(poses_shift + 1):end]), 3, n_points)

    n_outliers = 0
    @simd for i in 1:n_observations
        T = @view(poses[:, cache.poses_ids[i]])
        pt = @view(points[:, cache.points_ids[i]])
        pt = RotZYX(T[1:3]...) * pt .+ T[4:6]
        Δ = @view(cache.pixels[:, i]) .- project(camera, pt) # (y, x) format

        outlier = pt[3] < depth_ϵ || (Δ[1]^2 + Δ[2]^2) > repr_ϵ
        cache.outliers[i] = outlier
        outlier && (n_outliers += 1;)
    end
    n_outliers
end

function pnp_bundle_adjustment(
    camera::Camera, pose::SMatrix{4, 4, Float64}, pixels, points;
    iterations::Int = 10, show_trace::Bool = false,
    depth_ϵ::Real = 1e-6, repr_ϵ::Real = 5.0,
)
    R = RotZYX(pose[1:3, 1:3])
    X0 = [R.theta1, R.theta2, R.theta3, pose[1:3, 4]...]

    Y = zeros(Float64, length(pixels) * 2)
    outliers = fill(false, length(points))
    ignore_outliers = false

    function residue!(Y, X)
        @simd for i in 1:length(points)
            id = (i - 1) * 2
            if ignore_outliers && outliers[i]
                Y[id + 1] = 0.0
                Y[id + 2] = 0.0
            else
                pt = RotZYX(@view(X[1:3])...) * points[i] .+ @view(X[4:6])
                Y[(id + 1):(id + 2)] .= pixels[i] .- project(camera, pt)
            end
        end
    end

    residue!(Y, X0)
    initial_error = mapreduce(abs2, +, Y)
    Y .= 0.0

    fast_result = optimize!(
        LeastSquaresProblem(;x=X0, y=Y, f! = residue!),
        LevenbergMarquardt(); iterations=5, show_trace)
    X1 = fast_result.minimizer

    # Detect outliers using `fast_result` minimizer.
    n_outliers = 0
    @simd for i in 1:length(points)
        pt = RotZYX(@view(X1[1:3])...) * points[i] .+ @view(X1[4:6])
        r = pixels[i] .- project(camera, pt)
        outlier = pt[3] < depth_ϵ || (r[1]^2 + r[2]^2) > repr_ϵ
        outliers[i] = outlier
        outlier && (n_outliers += 1;)
    end

    if length(points) - n_outliers < 5
        return (
            SMatrix{4, 4, Float64}(I), initial_error, fast_result.ssr,
            outliers, n_outliers)
    end

    Y .= 0.0
    ignore_outliers = true
    result = optimize!(
        LeastSquaresProblem(;x=X1, y=Y, f! = residue!),
        LevenbergMarquardt(); iterations, show_trace)

    new_pose = to_4x4(RotZYX(result.minimizer[1:3]...), result.minimizer[4:6])
    new_pose, initial_error, result.ssr, outliers, n_outliers
end


text:
@symbol_func function Q_E(cur_reactor::AbstractReactor)
  cur_Q = cur_reactor.Q

  cur_Q *= 1 + E_Li / E_F

  cur_Q += 1

  cur_Q *= cur_reactor.eta_T

  cur_Q *= cur_reactor.eta_RF

  cur_Q -= 1

  cur_Q
end


text:
using ParticleFilters

pomdp = RockSample(4, 4)

filter = SIRParticleFilter(pomdp, 500)

custom_bounds = RockSampleBounds(pomdp) # custom lower bound to use with DESPOT solver

solver = DESPOTSolver{RockSampleState,
                      RockSampleAction,
                      RockSampleObs,
                      RockSampleBounds,
                      RandomStreams}(bounds = custom_bounds)

policy = solve(solver, pomdp)

ro = RolloutSimulator(max_steps=10, rng=MersenneTwister(1))

simulate(ro, pomdp, policy, filter)


text:
function value(obj::ExtendedUnivariateDifferentiableVectorFunction, arg0::Derivative)
    return jcall(obj, "value", Vector{Derivative}, (Derivative,), arg0)
end

function value(obj::UnivariateDifferentiableVectorFunction, arg0::DerivativeStructure)
    return jcall(obj, "value", Vector{DerivativeStructure}, (DerivativeStructure,), arg0)
end

function value(obj::UnivariateVectorFunction, arg0::jdouble)
    return jcall(obj, "value", Vector{jdouble}, (jdouble,), arg0)
end


text:
################################################################################
#
#                 LatEnum.jl : Basic lattice enumeration
#
# This file is part of hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# (C) 2015 Claus Fieker
# (C) 2015 Tommy Hofmann
#
################################################################################
#  TODO:
#   - (sh/c)ould be indexed by the type of G and C
#     in fact, since G is not used after the creation, maybe we should drop it.
#     realisticly, x (hence L, U) can be small types mostly
#   - support for lower bound
#   - reference vector
#   - support for all points
#   - if all points are wanted, the spiraling out, enumerating coordinates from
#     the center outwards is counter-productive
#   - lower bounds should be non-trivial speed-up by effectively generating the
#     L, U for the other bound as well and using this for exclusion.
#     (see other comment below)

add_assert_scope(:LatEnum)
add_verbose_scope(:LatEnum)

function show(io::IO, E::enum_ctx)
  println(io, "EnumCtx")
  if isdefined(E, :c)
    println(io, "curr. length ", E.c, " elt ", E.x, "(", (typeof(E.x), typeof(E.C), typeof(E.U)), ")")
  end
end

#need to only compute the top l x l submatrix when using limited enum
function pseudo_cholesky(G::fmpz_mat, den=1;
                 TC::Type=Rational{BigInt}, limit = nrows(G))
  n = ncols(G)
  @hassert :LatEnum 1 nrows(G) == n
  limit = min(limit, n)
  t = fmpz()
  C = Array{TC}(undef, limit, limit)
  for i=1:limit
    for j=1:limit
      getindex!(t, G, i, j)
      C[i,j] = TC(t//den)
    end
  end
  for i = 1:limit-1
    for j = i+1:limit
      C[j,i] = C[i,j]
      C[i,j] = C[i,j]/C[i,i]
    end
    for k = i+1:limit
      for l = i+1:limit
        C[k,l] = C[k,l] - C[k,i]*C[i,l]
      end
    end
  end
  for j = 1:limit-1
    if C[j,j] <= 0
      throw(LowPrecisionCholesky())
    end
    for i=j+1:limit
      C[i,j] = 0
    end
  end
  if C[limit, limit] <= 0
    throw(LowPrecisionCholesky())
  end
  return C
end

function enum_ctx_from_basis(B::FakeFmpqMat; Tx::Type = BigInt, TC::Type = Rational{BigInt}, TU::Type = Rational{BigInt}, limit = nrows(B))
  G = gram(numerator(B))
  return enum_ctx_from_gram(G, denominator(B)^2, Tx=Tx, TC=TC, TU=TU, limit = limit)
end

function enum_ctx_from_gram(G::FakeFmpqMat; Tx = BigInt, TC = Rational{BigInt}, TU = Rational{BigInt}, limit = nrows(G))
  return enum_ctx_from_gram(numerator(G), denominator(G), Tx=Tx, TC=TC, TU=TU, limit = limit)
end

function enum_ctx_from_basis(B::fmpz_mat, den::fmpz = fmpz(1); Tx::Type = BigInt, TC::Type = Rational{BigInt}, TU::Type = Rational{BigInt}, limit = nrows(B))
  G = gram(B)
  return enum_ctx_from_gram(G, den*den, Tx=Tx, TC=TC, TU=TU, limit = limit)
end

function enum_ctx_from_gram(G::fmpz_mat, den = 1; Tx = BigInt, TC = Rational{BigInt}, TU = Rational{BigInt}, limit = nrows(G))
  E = enum_ctx{Tx, TC, TU}()
  E.G = G
  n = nrows(G)
  E.n = n
  limit = min(limit, n)
  E.limit = limit
  E.d = den
  E.C = pseudo_cholesky(E.G, den, TC = TC, limit = limit)
  E.x = zero_matrix(FlintZZ, 1, n)
    #coeffs limit+1:n are going to be zero, always
  E.L = Vector{TU}(undef, limit) #lower and
  E.U = Vector{TU}(undef, limit) #upper bounds for the coordinates

  E.l = Vector{TU}(undef, limit) #current length
  E.tail = Vector{TU}(undef, limit)
  return E
end

function enum_ctx_local_bound(a::Rational{T}, b::Rational{T}) where T
  #return L <= U sth.
  #L = ceil(a-sqrt(b)), U = floor(a+sqrt(b))
  #solves (gives bounds) for (a-x)^2 <= b
  b >= 0 || return a, a-1
  @hassert :LatEnum 1 b >= 0
  d = denominator(b)
  i = isqrt(numerator(b*d*d))
  L = Base.ceil(a-(i+1)//d)
  U = Base.floor(a+(i+1)//d)
#  @show L, U, Base.ceil(BigFloat(a) - sqrt(BigFloat(b))), Base.floor(BigFloat(a) + sqrt(BigFloat(b)))
  if (a-L)^2 >b
    L +=1
  end
  if (a-U)^2>b
    U -= 1
  end

#println("local bound for ", a, " and ", b)
#println("is ", L, " and ", U)
  return L, U
end

function enum_ctx_local_bound(a::Number, b::Number) where Number
  #return L <= U sth.
  #L = ceil(a-sqrt(b)), U = floor(a+sqrt(b))
  #solves (gives bounds) for (a-x)^2 <= b
  b >= 0 || return a, a-1
  @hassert :LatEnum 1 b >= 0
  i = sqrt(b)
  L = Base.ceil(a-i)
  U = Base.floor(a+i)
#println("local bound for ", a, " and ", b)
#println("is ", L, " and ", U)
  return L, U
end

function enum_ctx_start(E::enum_ctx{A,B,C}, c::fmpz) where {A,B,C}
  E.c = c
  zero!(E.x)
  for i=1:E.limit
    E.l[i] = C(E.c//E.d)
    E.tail[i] = 0
    L, U = enum_ctx_local_bound(C(0), C(B(E.c//E.d)/E.C[i,i]))
    @hassert :LatEnum 1 typeof(L) == C
    @hassert :LatEnum 1 typeof(U) == C
    @hassert :LatEnum 1 typeof(E.L) == Vector{C}
    E.U[i] = U
    E.L[i] = L
  end
  E.U[1] = max(E.U[1], 1)
  E.L[1] = -E.U[1]
  E.last_non_zero = 1
  E.cnt = 0
end

function enum_ctx_start(E::enum_ctx{A,B,C}, x::fmpz_mat; eps::Float64=1.0) where {A,B,C}
  E.x = x
  for i=E.limit-1:-1:1
    E.tail[i] = sum(E.C[i, j]*C(E.x[1,j]) for j=i+1:E.limit)
  end
  b = sum(E.C[i,i]*(C(E.x[1,i]) + E.tail[i])^2 for i=1:E.limit) #curr. length
  #@show b, C((x*E.G*x')[1,1])
  #@assert b == C((x*E.G*x')[1,1])
  E.c = ceil(fmpz, b*C(E.d)*eps)
  E.l[E.limit] = C(E.c//E.d)
  for i=E.limit-1:-1:1
    E.l[i] = E.l[i+1] - E.C[i+1,i+1]*(C(E.x[1,i+1]) + E.tail[i+1])^2
  end
  for i=1:E.limit
    E.L[i], E.U[i] = enum_ctx_local_bound(-E.tail[i], E.l[i]/E.C[i,i])
  end
  E.last_non_zero = maximum(findall(i->E.x[1, i] != 0, 1:E.limit))+1
end

#for pol-red-abs we need s.th. else I think
#
#
#missing: lower bound
#         reference vector (eventually)
#         length
#         proper lattice type

@inline function fmpz_mat_entry(a::fmpz_mat, r::Int, c::Int)
  return ccall((:fmpz_mat_entry, libflint), Ptr{fmpz},
               (Ref{fmpz_mat}, Int, Int), a, r - 1, c - 1)
#@inline function fmpz_mat_entry(a::fmpz_mat, r::Int, c::Int)
    #return unsafe_load(reinterpret(Ptr{Ptr{fmpz}}, a.rows), r)+(c-1)*sizeof(Ptr)
end

@inline function fmpz_mat_entry_incref!(a::fmpz_mat, r::Int, c::Int)
  z = ccall((:fmpz_mat_entry, libflint), Ptr{fmpz},
               (Ref{fmpz_mat}, Int, Int), a, r - 1, c - 1)
  ccall((:fmpz_add_ui, libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}, Int), z, z, 1)
  #z = fmpz_mat_entry(a, r, c)
  #unsafe_store!(reinterpret(Ptr{Int}, z), unsafe_load(reinterpret(Ptr{Int}, z), 1)+1, 1)
  ##ccall((:fmpz_add_ui, libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}, Int), z, z, 1)
end

function fmpz_mat_entry_add_ui!(a::fmpz_mat, r::Int, c::Int, v::UInt)
  z = ccall((:fmpz_mat_entry, libflint), Ptr{fmpz},
               (Ref{fmpz_mat}, Int, Int), a, r - 1, c - 1)
  ccall((:fmpz_add_ui, libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}, Int), z, z, v)
end

function enum_ctx_advance_level(E::enum_ctx{A,B,C}, i::Int) where {A,B,C}
#  println("i: ", i, "                                   "[1:2*i], "|")
  t = fmpz()
  if i == E.last_non_zero-1
    fmpz_mat_entry_incref!(E.x, 1, i)
#    E.x[1, i] = getindex!(t, E.x, 1, i) + 1
  elseif i == E.last_non_zero
#    @hassert :LatEnum 1 E.x[1, i] == 0
    E.last_non_zero += 1
    fmpz_mat_entry_incref!(E.x, 1, i)
#    E.x[1, i] = getindex!(t, E.x, 1, i) + 1
  else
    s = E.U[i] + E.L[i]
    getindex!(t, E.x, 1, i)
    x = A(t)
    if s < 2*x  # larger than the middle
      setindex!(E.x, -x + A(Base.floor(s)), 1, i)
    else
      setindex!(E.x, -x + A(Base.floor(s)) + 1, 1, i)
    end
  end
end

_next = 0.0

function enum_ctx_next(E::enum_ctx{A,B,C}) where {A,B,C}
  global _next
  @v_do :ClassGroup_time 2 rt = time_ns()
  E.cnt += 1
  n::Int = 1
  n = E.limit
  i=1
  t = fmpz()
  while true
    enum_ctx_advance_level(E, i)
    getindex!(t, E.x, 1, i)
    if E.L[i] <= A(t) <= E.U[i] #coordinate is valid
      if i==1
        @v_do :ClassGroup_time 2 _next += time_ns()-rt
        return true
      else
        i -= 1
      end
    else # coordinate is invalid
      i += 1
      if i>n
        @v_do :ClassGroup_time 2 _next += time_ns()-rt
        return false
      end
      continue
    end

    @hassert :LatEnum 1 i < n
    while true
      getindex!(t, E.x, 1, i+1)
      t1 = A(t)
      E.tail[i] = E.C[i, i+1]*t1
      for j = i+2:n
        getindex!(t, E.x, 1, j)
        E.tail[i] += E.C[i, j] * A(t)
      end
      #E.tail[i] = sum C[i,j]x[j]
      E.l[i]    = E.l[i+1] - E.C[i+1, i+1]*(t1 + E.tail[i+1])^2
      #l[i] = l[i+1] - C[i+1, i+1]^2*(x[i+1] + tail[i+1])

      if E.l[i] >= 0
        Z = C(B(E.l[i])/E.C[i,i])
#        @hassert :LatEnum 1 E.l[i] >= 0
#        @hassert :LatEnum 1 E.C[i,i] > 0
#        @hassert :LatEnum 1 Z >= 0
        L, U = enum_ctx_local_bound(-E.tail[i], Z)
#        @hassert :LatEnum 1 typeof(L) == C
        E.L[i] = L
        E.U[i] = U

        x = A(Base.ceil((E.L[i] +E.U[i])/2))
        E.x[1, i] = x
        if -E.L[i] == E.U[i] && E.last_non_zero == i+1
          E.last_non_zero = i
#          @hassert :LatEnum 1 x == 0
        end
        if x <= E.U[i] # coordinate is valid
          i -= 1       # go further up
          if i==0
            @v_do :ClassGroup_time 2 _next += time_ns()-rt
            return true
          end
          continue
        else  # coordinate invalid, need to try next on i+1
          i += 1
          if i>n
            @v_do :ClassGroup_time 2 _next += time_ns()-rt
            return false
          end
          break
        end
      else # intervall too short
        i += 1
        if i>n
          @v_do :ClassGroup_time 2 _next += time_ns()-rt
          return false
        end
        break
      end
    end
  end
  @v_do :ClassGroup_time 2 _next += time_ns()-rt
  return true
end

function enum_ctx_short_elements(E::enum_ctx{A,B,C}, c::T, limit=-1) where {A,B,C} where T <: Union{Integer, fmpz}
  enum_ctx_start(E, fmpz(c))
  if enum_ctx_next(E)
    l = deepcopy(E.x) # else the 1st element is not returned....
  else
    l = matrix(FlintZZ, 0, E.n, fmpz[])
  end
  while enum_ctx_next(E) && (limit == -1 || limit >= Base.size(l, 1))
    l = vcat(l, E.x)
  end
  return l
end

################################################################################
#
#  Enumerating using arb objects
#
################################################################################

function pseudo_cholesky(G::arb_mat)
  n = ncols(G)
  C = deepcopy(G)
  for i = 1:n-1
    for j = i+1:n
      C[j,i] = C[i,j]
      C[i,j] = C[i,j]//C[i,i]
      if !isfinite(C[i,j])
          error("Precision not high enough")
      end
    end
    for k = i+1:n
      for l = i+1:n
        C[k,l] = C[k,l] - C[k,i]*C[i,l]
        if !isfinite(C[k,l])
          error("Precision not high enough")
        end
      end
    end
  end
  for j = 1:n-1
    for i=j+1:n
      C[i,j] = parent(C[i,j])(0)
    end
  end
  for i in 1:n
    #println(C[i,i])
    if contains_zero(C[i,i])
      error("Precision not high enough")
    end
  end

  return C
end

function enumerate_using_gram(G::arb_mat, c::arb)
  E = EnumCtxArb(pseudo_cholesky(G))
  return _enumerate(E, c, nrows(G), zero_matrix(FlintZZ, 1, nrows(G)))
end

function _enumerate(E::EnumCtxArb, c::arb, i::Int, x::fmpz_mat)
  # assume x_n,x_n-1,...,x_i+1 are set
  # compute the bound for x_i
  G = E.G
  n = nrows(G)
  p = E.p
  #recprint(n-i)
  #println("I am at level $i with $x")

  @vprint :LatEnum "$(recprint(n - i)) Gii: , $(G[i,i])\n"
  @vprint :LatEnum "$(recprint(n - i)) c: $(ArbField(p)(c))\n"

  C = ArbField(p)(c)//G[i,i]

  @vprint :LatEnum "$(recprint(n - i)) C: $C\n"

  if !contains_zero(C)
    C = sqrt(C)
  else
    # Element C contains zero, this is a bad when taking square root
    # (this would ive NaN).
    # The elements abs(C) also contains zero.
    # The only(?) way to do it is to invoke arb_get_abs_ubound_arf.
    # This should work, potentially we enumerate more elements

    tm = arf_struct(0, 0, 0, 0)
    ccall((:arf_init, libarb), Nothing, (Ref{arf_struct}, ), tm)

    ccall((:arb_get_abs_ubound_arf, libarb), Nothing, (Ref{arf_struct}, Ref{arb}), tm, C)
    ccall((:arb_set_arf, libarb), Nothing, (Ref{arb}, Ref{arf_struct}), C, tm)

    ccall((:arf_clear, libarb), Nothing, (Ref{arf_struct}, ), tm)

    @hassert :LatEnum !contains_zero(C)
  end

  @vprint :LatEnum "$(recprint(n - i)) C: $C\n"

  CC = ArbField(p)(0)

  for j in i+1:n
    CC = CC + G[i,j]*x[1,j]
  end

  lb = -CC - C
  ub = -CC + C

  tr_ptr = ccall((:arb_rad_ptr, libarb), Ptr{Nemo.mag_struct}, (Ref{arb}, ), lb)

  tm_ptr = ccall((:arb_mid_ptr, libarb), Ptr{arf_struct}, (Ref{arb}, ), lb)
  u = arf_struct(0, 0, 0, 0)
  ccall((:arf_init, libarb), Nothing, (Ref{arf_struct}, ), u)

  ccall((:arf_set_mag, libarb), Nothing, (Ref{arf_struct}, Ptr{Nemo.mag_struct}), u, tr_ptr)
  ccall((:arf_sub, libarb), Nothing, (Ref{arf_struct}, Ptr{arf_struct}, Ref{arf_struct}, Int, Cint), u, tm_ptr, u, p, 4) # 4 is round to -infty
  lbfmpz = fmpz()
  ccall((:arf_get_fmpz, libarb), Nothing, (Ref{fmpz}, Ref{arf_struct}, Cint), lbfmpz, u, 4)

  tr = ccall((:arb_rad_ptr, libarb), Ptr{Nemo.mag_struct}, (Ref{arb}, ), ub)
  tm = ccall((:arb_mid_ptr, libarb), Ptr{arf_struct}, (Ref{arb}, ), ub)

  ccall((:arf_set_mag, libarb), Nothing, (Ref{arf_struct}, Ptr{Nemo.mag_struct}), u, tr)
  ccall((:arf_sub, libarb), Nothing, (Ref{arf_struct}, Ptr{arf_struct}, Ref{arf_struct}, Int, Cint), u, tm, u, p, 3) # 3 is round to +infty
  ubfmpz = fmpz()
  ccall((:arf_get_fmpz, libarb), Nothing, (Ref{fmpz}, Ref{arf_struct}, Cint), ubfmpz, u, 3)

  ccall((:arf_clear, libarb), Nothing, (Ref{arf_struct}, ), u)

  @vprint :LatEnum "$(recprint(n - i)) Coordinate $i between $lbfmpz and $ubfmpz\n"

  A = Vector{Vector{fmpz}}()

  if i == 1
    @vprint :LatEnum "$(recprint(n - i)) his is depth $i\n"

    for j in Int(lbfmpz):Int(ubfmpz)
      @vprint :LatEnum "$(recprint(n - i)) j is $j\n"
      if isnegative(c - G[i,i]*(j + CC)^2)
        @vprint :LatEnum "$(recprint(n - i)) Something is negative\n"
        continue
      end
      push!(A, [ fmpz(j) ])
    end
  else
    for j in Int(lbfmpz):Int(ubfmpz)
      @vprint :LatEnum "$(recprint(n - i)) $x $i $j: $c $(G[i,i]) $CC\n"
      @vprint :LatEnum "$(recprint(n - i)) $x $i $j $(G[i,i]*j^2)\n"
      @vprint :LatEnum "$(recprint(n - i)) $x $i $j: $(c - G[i,i]*(j + CC)^2)\n"
      t = c - G[i,i]*(j + CC)^2

      if isnegative(t)
        continue
      end

      x[1,i] = j
      l = _enumerate(E, t, i - 1, x)
      for k in 1:length(l)
        if n == length(l[k]) + 1 && iszero(j) && all(iszero, l[k])
          continue
        end
        push!(A, push!(l[k], fmpz(j)))
      end
    end
  end

  return A
end

function recprint(n::Int)
  s = "  "^n
  return s
end

@doc Markdown.doc"""
    hermite_constant(n::Int, R=ArbField(100))

The $n$-th Hermite constant (for lattices) to the power $n$.
Exact value if $1 \le n \le 8$ and Blichfeld's upper bound
otherwise.
"""
function hermite_constant(n::Int, R=ArbField(20))
  if n <= 8
    return R([1, 4/3, 2, 4, 8, 64/3, 64, 256][n])
  end
  return (2/pi)^n*gamma(R(2+n/2))^2
end

#using semi_factorial and gamma_half we can write the Blichfeld bound
#as rational * pi^1-n (if we want)
#Blichfeld = pi^(1-n) /8 * semi_factorial(n+2)^2
#
function semi_factorial(n::Int)
  if isodd(n)
    return prod([fmpz(2*i-1) for i=1:div(n+1, 2)])
  else
    return prod([fmpz(2*i) for i=1:div(n, 2)])
  end
end

function gamma_half(n::Int)
  #gamma(n/2)
  return sqrt(pi)*semi_factorial(n-2)/2^((n-1)/2)
end