text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advection_velocity = (0.2, -0.7)
equations = LinearScalarAdvectionEquation2D(advection_velocity)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# Deformed rectangle that looks like a waving flag,
# lower and upper faces are sinus curves, left and right are vertical lines.
f1(s) = SVector(-1.0, s - 1.0)
f2(s) = SVector( 1.0, s + 1.0)
f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
f4(s) = SVector(s,  1.0 + sin(0.5 * pi * s))

# Create P4estMesh with 3 x 2 trees,
# approximate the geometry with a smaller polydeg for testing.
cells_per_dimension = (3, 2)
mesh = P4estMesh(cells_per_dimension, polydeg=3,
                 faces=(f1, f2, f3, f4))

# Refine bottom left quadrant of each tree to level 4
function refine_fn(p4est, which_tree, quadrant)
  if quadrant.x == 0 && quadrant.y == 0 && quadrant.level < 4
    # return true (refine)
    return Cint(1)
  else
    # return false (don't refine)
    return Cint(0)
  end
end

# Refine recursively until each bottom left quadrant of a tree has level 4
# The mesh will be rebalanced before the simulation starts
refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))
Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL)
Trixi.balance!(mesh)

# Refine everything again for convergence tests
refine_fn2(p4est, which_tree, quadrant) = Cint(1)
refine_fn_c2 = @cfunction(refine_fn2, Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))

# Abuse the variable `initial_refinement_level` for this because this will be increased
# during convergence tests
initial_refinement_level = 0
for i in 1:initial_refinement_level
  Trixi.refine_p4est!(mesh.p4est, false, refine_fn_c2, C_NULL)
end

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0));

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100,
                                     solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=1.2)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

# Print the timer summary
summary_callback()


text:
module EnsembleKalmanProcesses

using Distributions, Statistics, LinearAlgebra, DocStringExtensions

#auxilliaries
include("ParameterDistribution.jl")
include("DataStorage.jl")
include("Observations.jl")

#updates:
include("EnsembleKalmanProcess.jl")

end # module

ext: jl


text:
# This file includes an example file by breaking algebraic loops by solving loop equation numerically.

using Causal
using Plots

# Simulation parameter
α = 3.

# Construct model with algebraic loop
@defmodel model begin
    @nodes begin
        gen = RampGenerator()
        adder = Adder(signs=(+,-))
        gain = Gain(gain=α)
        writer = Writer(input=Inport(2))
    end
    @branches begin
        gen[1]      =>      adder[1]
        adder[1]    =>      gain[1]
        gain[1]     =>      adder[2]
        gen[1]      =>      writer[1]
        gain[1]     =>      writer[2]
    end
end

# Simulate the model
simulate!(model, 0., 1., 100.)

# Plot the results
t, y = read(getnode(model, :writer).component)
yt = α / (α + 1) * getnode(model, :gen).component.readout.(t)
err = yt - y[:, 2]
p1 = plot(t, y[:, 1], label=:u)
    plot!(t, y[:, 2], label=:y)
    plot!(t, yt, label=:true)
p2 = plot(t, err, label=:err)
plot(p1, p2, layout=(2, 1))


text:
################################################################################
#
#  Parent constructor
#
################################################################################

function SRowSpace(R::Ring; cached = true)
  T = elem_type(R)
  return SRowSpace{T}(R)
end

base_ring(A::SRow) = parent(A.values[1])

==(x::SRow{T}, y::SRow{T}) where {T} = (x.pos == y.pos) && (x.values == y.values)

################################################################################
#
#  Sparse row creation
#
################################################################################

function SRow(R::Ring)
  return SRow{elem_type(R)}()
end

################################################################################
#
#  Hashing
#
################################################################################

function hash(A::SRow, h::UInt)
  return hash(A.pos, hash(A.values, h))
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::SRow{T}) where T
  print(io, "sparse row with positions $(A.pos) and values $(A.values)\n")
end

################################################################################
#
#  Copy
#
################################################################################

function copy(A::SRow{T}) where T
  sr = SRow{T}()
  for (p, v) = A
    push!(sr.pos, p)
    push!(sr.values, v)
  end
  return sr
end

################################################################################
#
#  Zero row testing
#
################################################################################

function iszero(A::SRow)
  return length(A.pos) == 0
end

################################################################################
#
#  Modular reduction
#
################################################################################

doc"""
***
    mod!(A::SRow{fmpz}, n::Integer)

> Inplace reduction of all entries of $A$ modulo $n$ to the positive residue
> system.
"""
function mod!(A::SRow{fmpz}, n::Integer)
  i=1
  while i<= length(A.pos)
    v = mod(A.values[i], n)
    if v == 0
      deleteat!(A.pos, i)
      deleteat!(A.values, i)
    else
      A.values[i] = v
      i += 1
    end
  end
end

doc"""
***
    mod!(A::SRow{fmpz}, n::fmpz)

> Inplace reduction of all entries of $A$ modulo $n$ to the positive residue
> system.
"""
function mod!(A::SRow{fmpz}, n::fmpz)
  i=1
  while i<= length(A.pos)
    v = mod(A.values[i], n)
    if v == 0
      deleteat!(A.pos, i)
      deleteat!(A.values, i)
    else
      A.values[i] = v
      i += 1
    end
  end
end

# Todo: Do not convert to fmpz
doc"""
***
    mod_sym!(A::SRow{fmpz}, n::Integer)

> Inplace reduction of all entries of $A$ modulo $n$ to the symmetric residue
> system.
"""
function mod_sym!(A::SRow{fmpz}, n::Integer)
  mod_sym!(A, fmpz(n))
end

doc"""
***
    mod_sym!(A::SRow{fmpz}, n::fmpz)

> Inplace reduction of all entries of $A$ modulo $n$ to the symmetric residue
> system.
"""
function mod_sym!(A::SRow{fmpz}, n::fmpz)
  i=1
  while i<= length(A.pos)
    v = mod_sym(A.values[i], n)
    if v == 0
      deleteat!(A.pos, i)
      deleteat!(A.values, i)
    else
      A.values[i] = v
      i += 1
    end
  end
end

################################################################################
#
#  Change ring
#
################################################################################

function SRow(A::SRow{fmpz}, n::Int)
  R = ResidueRing(FlintZZ, n, cached=false)
  return SRow(A, R)
end

function SRow(A::SRow{fmpz}, R::T) where T <: Ring
  z = SRow{elem_type(R)}()
  for (i, v) in A
    nv = R(v)
    if iszero(nv)
      continue
    else
      push!(z.pos, i)
      push!(z.values, nv)
    end
  end
  return z
end

doc"""
    SMat(A::SMat{fmpz}, n::Int) -> SMat{Generic.Res{fmpz}}
    SRow(A::SMat{fmpz}, n::Int) -> SRow{Generic.Res{fmpz}}

> Converts $A$ to ba a sparse matrix (row) over $Z/nZ$
"""
function SRow(A::SRow{fmpz}, n::fmpz)
  R = ResidueRing(FlintZZ, n, cached=false)
  return SRow(A, R)
end

################################################################################
#
#  Make sparse row iteratable
#
################################################################################

function length(A::SRow)
  return length(A.pos)
end

function endof(A::SRow)
  return length(A.pos)
end

function start(A::SRow)
  return 1
end

function next(A::SRow, st::Int)
  return (A.pos[st], A.values[st]), st + 1
end

function done(A::SRow, st::Int)
  return st > length(A.pos)
end

################################################################################
#
#  Multiplication
#
################################################################################

function mul(A::SRow{T}, B::SRow{T}) where T
  @assert length(A) != 0
  v = 0*A.values[1]
  b = 1
  for a=1:length(A.pos)
    while b<=length(B.values) && B.pos[b] < A.pos[a]
      b += 1
    end
    if b>length(B.values)
      return v
    end
    if B.pos[b] == A.pos[a]
      v += B.values[b] * A.values[a]
    end
  end
  return v
end

#in-place scaling
function scale_row!(a::SRow{fmpz}, b::fmpz)
  @assert !iszero(b)
  if isone(b)
    return
  end
  for i=1:length(a.pos)
    a.values[i] *= b
  end
end

################################################################################
#
#  Addition
#
################################################################################

function +(A::SRow{T}, B::SRow{T}) where T
  if length(A.values) == 0
    return B
  elseif length(B.values) == 0
    return A
  end
  return add_scaled_row(A, B, base_ring(A)(1))
end

################################################################################
#
#  Scalar multiplication
#
################################################################################

# The following functions runs into an infinite recursion in case T = fmpz
#function *{T}(b::fmpz, A::SRow{T})
#  r = base_ring(A)(b)
#  return r*A
#end

function *(b::T, A::SRow{T}) where T
  B = SRow{T}()
  if iszero(b)
    return B
  end
  for (p,v) = A
    nv = v*b
    if !iszero(nv)  # there are zero divisors - potentially
      push!(B.pos, p)
      push!(B.values, nv)
    end
  end
  return B
end

function *(b::Integer, A::SRow{T}) where T
  if length(A.values) == 0
    return deepcopy(A)
  end
  return base_ring(A)(b)*A
end

function div(A::SRow{T}, b::T) where T
  B = SRow{T}()
  if iszero(b)
    return error("Division by zero")
  end
  for (p,v) = A
    nv = div(v, b)
    if !iszero(nv)
      push!(B.pos, p)
      push!(B.values, nv)
    end
  end
  return B
end

function div(A::SRow{T}, b::Integer) where T
  if length(A.values) == 0
    return deepcopy(A)
  end
  return div(A, base_ring(A)(b))
end

function divexact(A::SRow{T}, b::T) where T
  B = SRow{T}()
  if iszero(b)
    return error("Division by zero")
  end
  for (p,v) = A
    nv = divexact(v, b)
    @assert !iszero(nv)
    push!(B.pos, p)
    push!(B.values, nv)
  end
  return B
end

function divexact(A::SRow{T}, b::Integer) where T
  if length(A.values) == 0
    return deepcopy(A)
  end
  return divexact(A, base_ring(A)(b))
end

################################################################################
#
#  Maximum and minimum
#
################################################################################

doc"""
***
    maximum(A::SRow{fmpz}) -> fmpz

> Finds the largest entry of $A$.
"""
function maximum(A::SRow{fmpz})
  return maximum(A.values)
end

doc"""
***
    minimum(A::SRow{fmpz}) -> fmpz

> Finds the smallest entry of $A$.
"""
function minimum(A::SRow{fmpz})
  return minimum(A.values)
end

################################################################################
#
#  Lifting
#
################################################################################

doc"""
    lift(a::SRow{nmod}) -> SRow{fmpz}

> Lifts all entries in $a$.
"""
function lift(a::SRow{nmod})
  b = SRow{fmpz}()
  for (p,v) = a
    push!(b.pos, p)
    push!(b.values, lift(v))
  end
  return b
end


text:
module EnsembleKombine

using ..Stats

using LinearAlgebra

import Base:
    rand
using Random

using Statistics

"""
    rescale(pts)

Transform ``pts`` (of size ``(ndim, npts)``) into zero-mean, unit covaraince.
"""
function rescale(pts)
    nd, np = size(pts)

    mu = vec(mean(pts, dims=2))
    cv = cov(pts, dims=2)

    repts = zeros(nd, np)
    for j in 1:np
        for i in 1:nd
            repts[i,j] = pts[i,j] - mu[i]
        end
    end

    F = cholesky(cv)

    # cv = L*L^T
    # so if x is N(0,1) then
    # y = Lx has covariance cv: y y^T = L x x^T L^T = cv

    F.L \ repts
end

"""
    update_means(pts, mus, assigns)

Part of the K-means iteration, returns updated means representing the
centroid of each cluster.
"""
function update_means(pts, mus, assigns)
    ncluster = maximum(assigns)

    new_mus = zeros(size(pts, 1), ncluster)
    for i in 1:ncluster
        ps = pts[:, assigns.==i]
        new_mus[:,i] = mean(ps, dims=2)
    end

    new_mus
end

"""
    update_assigns(pts, mus, assigns)

Part of the K-means iteration, returns updated cluster assignments
based on the nearst mean position to each point.
"""
function update_assigns(pts, mus, assigns)
    nd, np = size(pts)

    new_assigns = zeros(Int, size(pts,2))
    for j in 1:np
        p = pts[:,j]

        dmin = Inf
        imin = 0
        for k in 1:size(mus,2)
            d2 = sum((p.-mus[:,k]).*(p.-mus[:,k]))
            if d2 < dmin
                dmin = d2
                imin = k
            end

            new_assigns[j] = imin
        end
    end

    new_assigns
end

"""
    kmeans_init(pts, n)

Returns ``mus, assigns`` to initialise the k-means iteration, using
the locations of randomly chosen points for the ``n`` cluster centres.

The argument ``pts`` should have shape ``(ndim, npts)``.
"""
function kmeans_init(pts, n)
    nd, np = size(pts)

    mus = zeros(nd, n)
    assigns = ones(Int, np)

    inds = randperm(np)
    for j in 1:n
        mus[:,j] = pts[:,inds[j]]
    end

    mus, update_assigns(pts, mus, assigns)
end

"""
    kmeans(pts, n, [maxiter=10000])

Returns ``(mus, assigns)`` after iterating the two K-means update
steps to convergence.  Because the initial assignment of cluster
centres is random, running the code again may lead to a different
solution.
"""
function kmeans(pts, n, maxiter=10000)
    mus, assigns = kmeans_init(pts, n)

    new_mus = nothing
    new_assigns = nothing
    i = 0
    while true
        new_mus = update_means(pts, mus, assigns)
        new_assigns = update_assigns(pts, new_mus, assigns)

        if new_mus != mus
            mus = new_mus
            assigns = new_assigns
        elseif i > maxiter
            break
        else
            break
        end

        i += 1
    end

    new_mus, new_assigns
end

"""
    ClusteredKDE(pts, n)

Returns an object representing the a clustered KDE with ``n``
clusters.  The ``pts`` argument should have shape ``(ndim, npts)``.
Methods defined on this object include

* ``npts(ck)``, ``ndim(ck)``, and ``ncl(ck)`` giving the number of
  points, dimensions, and clusters of the object.

* ``logpdf(ck, x)`` which returns the log of the KDE probability
  density at the point ``x``.

* ``rand(ck[, dims...])`` which returns a random draw from the KDE PDF
  (or an array of such draws).
"""
struct ClusteredKDE
    pts::Array{Float64, 2}
    assigns::Array{Int, 1}
    cholfacts::Array{LinearAlgebra.Cholesky{Float64,Array{Float64,2}}, 1}
end

"""
    ClusteredKDE(pts, n)

Construct a ``ClusteredKDE`` object with ``n`` clusters from the given
``pts``.
"""
function ClusteredKDE(pts, n)
    nd, np = size(pts)

    _, assigns = kmeans(pts, n)

    cfacts = LinearAlgebra.Cholesky{Float64,Array{Float64,2}}[]
    for j in 1:n
        sel = assigns.==j
        ps = pts[:,sel]

        @assert size(ps,2)>n "one cluster has too few points"

        cv = cov(ps, dims=2)
        cv *= (1.0/size(ps, 2)).^(2.0/(4.0 + nd)) # Scott's rule

        push!(cfacts, cholesky(cv))
    end

    ClusteredKDE(copy(pts), assigns, cfacts)
end

"""The number of dimensions in the KDE's points."""
function ndim(ck::ClusteredKDE)
    size(ck.pts, 1)
end

"""The number of points used to construct the KDE."""
function npts(ck::ClusteredKDE)
    size(ck.pts, 2)
end

"""The number of clusters in the KDE."""
function ncl(ck::ClusteredKDE)
    size(ck.cholfacts, 1)
end

"""
    rand(ck::ClusteredKDE, [n::Int])

Return a random point or array of such points from the clustered KDE.
"""
function rand(ck::ClusteredKDE)
    j = rand(1:npts(ck))
    p = ck.pts[:, j]
    F = ck.cholfacts[ck.assigns[j]]

    p + F.L*randn(ndim(ck))
end

function rand(ck::ClusteredKDE, n::Int)
    pts = zeros(ndim(ck), n)

    for j in 1:n
        pts[:,j] = rand(ck)
    end

    pts
end

"""
    logpdf(ck::ClusteredKDE, x)

Return the log of the PDF defined by the clustered KDE at the point
``x``.
"""
function logpdf(ck, x)
    logdets = [sum(log.(diag(ck.cholfacts[i].L))) for i in 1:ncl(ck)]

    lpdf = -Inf
    for j in 1:npts(ck)
        p = ck.pts[:,j]
        a = ck.assigns[j]
        F = ck.cholfacts[a]

        dx = x - p
        r2 = sum(dx .* (F \ dx))

        l = -0.5*ndim(ck)*log(2.0*pi) -logdets[a] - 0.5*r2

        lpdf = logsumexp(lpdf, l)
    end
    lpdf - log(npts(ck))
end

"""
    build_proposal_kde(pts, [ncmax])

Produces a KDE with at most ``ncmax`` clusters that is optimised to
represent the distribution of ``pts``.
"""
function build_proposal_kde(pts, ncmax=nothing)
    # We need to have at least ndim+1 points in each cluster, so that
    # limits how many clusters we can have.  IF all points divided
    # equally, then this is the maximum number of clusters

    nd, np = size(pts)

    @assert np%2==0 "need even number of points!"

    if ncmax == nothing
        ncmax = floor(Int, np / (2*(nd+1)))
    end

    perm = randperm(np)
    ihalf = round(Int, np/2)
    kde_pts = pts[:, perm[1:ihalf]]
    test_pts = pts[:, perm[ihalf+1:end]]

    best = -Inf
    bestck = nothing
    for nc in 1:ncmax
        for j in 1:5 # Five trials at each nc
            try
                ck = ClusteredKDE(kde_pts, nc)
                logprob = sum([logpdf(ck, test_pts[:,k]) for k in 1:size(test_pts,2)])

                if logprob > best
                    best = logprob
                    bestck = ck
                end
            catch x
                if isa(x, AssertionError) || isa(x, LinearAlgebra.PosDefException)
                    # Do nothing
                else
                    rethrow()
                end
            end
        end
    end

    bestck
end

"""
    mcmc_step(pts, lnprobs, lnprops, logpost, proposal)

Returns ``(new_pts, new_lnprobs, new_lnprops)`` for a single MCMC step
using the clustered KDE ``proposal`` as a Metropoli-Hastings proposal.

# Arguments

- ``pts`` of shape ``(ndim, npts)`` is the current state of the MCMC
  chain.

- ``lnprobs`` is the log-probability at the current points.

- ``lnprops`` is the log of the proposal density at the current
  points.

- ``logpost`` is the function that computes the log-probability for
  the MCMC.

- ``proposal::ClusteredKDE`` is the clustered KDE proposal.
"""
function mcmc_step(pts, lnprobs, lnprop, logpost, proposal)
    new_pts = rand(proposal, size(pts, 2))

    new_lnprobs = [logpost(new_pts[:,j]) for j in 1:size(new_pts,2)]

    new_lnprop = [logpdf(proposal, new_pts[:,j]) for j in 1:size(new_pts,2)]

    log_pacc = new_lnprobs .+ lnprop .- lnprobs .- new_lnprop

    out_pts = zeros(size(pts)...)
    out_lnprobs = zeros(size(lnprobs,1))
    out_lnprop = zeros(size(lnprop, 1))

    for j in 1:size(pts, 2)
        if log(rand()) < log_pacc[j]
            out_pts[:,j] = new_pts[:,j]
            out_lnprobs[j] = new_lnprobs[j]
            out_lnprop[j] = new_lnprop[j]
        else
            out_pts[:,j] = pts[:,j]
            out_lnprobs[j] = lnprobs[j]
            out_lnprop[j] = lnprop[j]
        end
    end

    (out_pts, out_lnprobs, out_lnprop)
end

"""
    run_mcmc(pts, lnprobs, lnprops, logpost, proposal, steps; [thin=1])

Repeatedly calls ``mcmc_step`` (see arguments above) ``steps`` times,
saving every ``thin`` iterations into a chain, and returning ``(chain,
chain_lnprobs, chain_lnprops)``.
"""
function run_mcmc(pts, lnprobs, lnprop, logpost, proposal, steps; thin=1)
    nsave = div(steps, thin)

    chain = zeros(size(pts,1), size(pts,2), nsave)
    chainlnprob = zeros(size(pts,2), nsave)
    chainlnprop = zeros(size(pts,2), nsave)
    for i in 1:steps
        pts, lnprobs, lnprop = mcmc_step(pts, lnprobs, lnprop, logpost, proposal)

        if i % thin == 0
            chain[:,:,div(i,thin)] = pts
            chainlnprob[:,div(i,thin)] = lnprobs
            chainlnprop[:,div(i,thin)] = lnprop
        end
    end

    (chain, chainlnprob, chainlnprop)
end

end


text:
export BrickletPiezoSpeakerIdentity
struct BrickletPiezoSpeakerIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletPiezoSpeaker
"""
Creates beep with configurable frequency
"""
mutable struct BrickletPiezoSpeaker <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletPiezoSpeaker(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_piezo_speaker")
        deviceInternal = package.BrickletPiezoSpeaker(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export beep
"""
    $(SIGNATURES)

Beeps with the given frequency for the given duration.

.. versionchanged:: 2.0.2\$nbsp;(Plugin)
   A duration of 0 stops the current beep if any, the frequency parameter is
   ignored. A duration of 4294967295 results in an infinite beep.

The Piezo Speaker Bricklet can only approximate the frequency, it will play
the best possible match by applying the calibration (see :func:`Calibrate`).
"""
function beep(device::BrickletPiezoSpeaker, duration, frequency)
    device.deviceInternal.beep(duration, frequency)
end

export morse_code
"""
    $(SIGNATURES)

Sets morse code that will be played by the piezo buzzer. The morse code
is given as a string consisting of "." (dot), "-" (minus) and " " (space)
for *dits*, *dahs* and *pauses*. Every other character is ignored.

For example: If you set the string "...---...", the piezo buzzer will beep
nine times with the durations "short short short long long long short
short short".
"""
function morse_code(device::BrickletPiezoSpeaker, morse, frequency)
    device.deviceInternal.morse_code(morse, frequency)
end

export calibrate
"""
    $(SIGNATURES)

The Piezo Speaker Bricklet can play 512 different tones. This function
plays each tone and measures the exact frequency back. The result is a
mapping between setting value and frequency. This mapping is stored
in the EEPROM and loaded on startup.

The Bricklet should come calibrated, you only need to call this
function (once) every time you reflash the Bricklet plugin.

Returns *true* after the calibration finishes.
"""
function calibrate(device::BrickletPiezoSpeaker)
    return device.deviceInternal.calibrate()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletPiezoSpeaker)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletPiezoSpeaker, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c4e34a90-77b8-11eb-108e-93d91efc7480
md"# Solving Tic Tac Toe With MiniMax Trees!"

# ╔═╡ a6e0b518-4f39-45d0-b2a3-a95cb1de7114
md"First we need to build out our minimax tree for the game"

# ╔═╡ 72fc2712-61b8-43a6-b68e-b3c5c87ed0a1
struct GameState
	turn::Int8
	tiles::Array{Int8, 2}
end

# ╔═╡ e6e3a756-6526-483d-af9e-9e53cf6d9f52
Base.:(==)(x::GameState, y::GameState) = x.turn == y.turn && all(x.tiles .== y.tiles)

# ╔═╡ b3d7a96e-f3fc-4dbd-94a9-552b21426396
Base.hash(x::GameState) = hash(x.tiles)

# ╔═╡ 82323b7f-1e92-46b0-8cff-5d85152d7d5c
md"0's are empty, 1's are player 1 (x), and 2's are player 2 (o)

Our objective is to optimize the policy of player 1"

# ╔═╡ a51583a3-0faa-4bac-a60b-e763644f674d
EMPTY = 0x00

# ╔═╡ eb2f779f-c651-4444-a34d-42681de9f3c6
PLAYER_1 = Int8(1)

# ╔═╡ f0e3fa72-9acf-4325-9173-622f93497b68
PLAYER_2 = Int8(-1)

# ╔═╡ 2a78b1ad-63cb-4cfa-98f4-a11ab05a9022
next_player(p) = p == PLAYER_1 ? PLAYER_2 : PLAYER_1

# ╔═╡ d1efdbe8-1a18-4a1c-84fc-f80dc3b9272d
initial_tiles = fill(EMPTY, 3, 3)

# ╔═╡ e0e9cbb2-01a8-4523-a6a9-ec7cfdb351d5
initial_state = GameState(PLAYER_1, initial_tiles)

# ╔═╡ e9dd6424-1031-44fb-970c-15aa3a4a1385
struct MinimaxNode
	state::GameState
	children::Vector{MinimaxNode}
	winner::Int8
	leaf::Bool
end

# ╔═╡ 5feafefc-740c-4450-838d-9e2ac47ceb32
function get_winner(state::GameState)
	# Really ugly check for game winners
	# Rows
	rows = [all(state.tiles[i, :] .== state.tiles[i, 1]) ? state.tiles[i, 1] : 0 for i ∈ 1:size(state.tiles)[1]]
	cols = [all(state.tiles[:, j] .== state.tiles[1, j]) ? state.tiles[1, j] : 0 for j ∈ 1:size(state.tiles)[2]]
	diag1 = state.tiles[1, 1] == state.tiles[2, 2] && state.tiles[1, 1] == state.tiles[3, 3] ? state.tiles[1, 1] : 0
	diag2 = state.tiles[3, 1] == state.tiles[2, 2] && state.tiles[3, 1] == state.tiles[1, 3] ? state.tiles[3, 1] : 0

	max_rows = max(rows..., cols..., diag1, diag2)
	min_rows = min(rows..., cols..., diag1, diag2)

	return max_rows > 0 ? max_rows : min_rows
end

# ╔═╡ 7b143cf1-920d-4626-a927-bf33885ae7a9
get_winner(GameState(PLAYER_1, Int8[
	1 -1 1;
	-1 -1 1;
	-1 -1 0;
]))

# ╔═╡ 31ff651f-b275-46d2-97e9-206c9926dfb8
function possible_moves(state::GameState)
	moves = []
	for i ∈ 1:length(state.tiles)
		if state.tiles[i] == EMPTY
			new_tiles = copy(state.tiles)
			new_tiles[i] = state.turn
			push!(moves, GameState(next_player(state.turn), new_tiles))
		end
	end

	moves
end

# ╔═╡ 258017ed-4bf0-4fdd-8239-2fc381e94d0c
function build_minimax(state::GameState)
	winner = get_winner(state)
	nodes = []
	for pm ∈ possible_moves(state)
		push!(nodes, build_minimax(pm))
	end

	if length(nodes) > 0 && winner == EMPTY
		return MinimaxNode(state, nodes, (state.turn == PLAYER_1 ? max : min)((x->x.winner).(nodes)...), true)
	else
		return MinimaxNode(state, [], winner, true)
	end
end

# ╔═╡ 45d2cddd-adfe-470a-89bc-56fb5d5982dc
begin
	tree = build_minimax(initial_state)
	nothing
end

# ╔═╡ d607a5e3-78f5-4308-999f-bf49858004a7
(x->x.winner).(tree.children[1].children[4].children[4].children[2].children)

# ╔═╡ c9254d06-75b9-47dc-8d29-495c80c9215d
tree.children[1].children[2].state

# ╔═╡ 15fd5bd7-d2c4-4337-a2ad-0fc43ad2b4e0
md"## Flattening our Tree
Tic Tac Toe has the nice property that if a game reaches a particular state, regardless of how that state was obtained, the minimax tree from that state will be the same. Therefore we can build up a dictionary of game states for keys, and the corresponding value is yet another game state representing the best move."

# ╔═╡ 248c158c-41bc-478f-be12-b7e4ab43544d
function flatten_minimax!(flat_tree, tree)
	flat_tree[tree.state] = tree.children[argmin((x->x.winner).(tree.children))].state
	for node ∈ tree.children
		if get_winner(node.state) == EMPTY && length(node.children) > 0
			flatten_minimax!(flat_tree, node)
		end
	end
end

# ╔═╡ 72c6bdcd-6d6b-4bf4-a759-1c4256384352
begin
	flat_tree = Dict()
	flatten_minimax!(flat_tree, tree)
end

# ╔═╡ c0634c74-9228-4d05-8146-3492aaf3a8ef
flat_tree[initial_state]

# ╔═╡ decac4e3-d5d3-4a51-a5bb-fdb658dec971
md"Fantastic! As can be seen above, the minimax tree has determined that the best move for player 1 is the center square!"

# ╔═╡ 2bb59351-cb34-45a7-bd58-72409adf48fa
html"

"

# ╔═╡ 4aa39de7-db17-42d7-924f-abb486c03eed
md"## Building the backend"

# ╔═╡ d5ba2170-d687-4ef8-a7fc-e8fef48b3a2c
md"Set **tiles** and **player** from the frontend"

# ╔═╡ b385f71f-c608-4be7-a0b6-30a39397451e
player = PLAYER_2

# ╔═╡ c1480936-bc64-481d-a754-31d5a35a2d6d
md"Note that tiles will be a vector of vectors, not an array"

# ╔═╡ c1cd29c9-b3dc-497b-b200-e6ad08356b39
tiles = [[0, 0, 0], [0, 1, 0], [0, 0, 0]]

# ╔═╡ 7d6a905a-2445-4569-8fe1-8e703b635218
md"And request **next_state**"

# ╔═╡ 4b879529-0359-462e-b032-3bb95a11ac36
next_state = flat_tree[GameState(player, hcat(tiles...))].tiles

# ╔═╡ 3e54ccec-d960-4cc2-a34a-7fce0964addb
md"Viola! Our \"backend\" is complete!"

# ╔═╡ Cell order:
# ╟─c4e34a90-77b8-11eb-108e-93d91efc7480
# ╟─a6e0b518-4f39-45d0-b2a3-a95cb1de7114
# ╠═72fc2712-61b8-43a6-b68e-b3c5c87ed0a1
# ╠═e6e3a756-6526-483d-af9e-9e53cf6d9f52
# ╠═b3d7a96e-f3fc-4dbd-94a9-552b21426396
# ╟─82323b7f-1e92-46b0-8cff-5d85152d7d5c
# ╠═a51583a3-0faa-4bac-a60b-e763644f674d
# ╠═eb2f779f-c651-4444-a34d-42681de9f3c6
# ╠═f0e3fa72-9acf-4325-9173-622f93497b68
# ╠═2a78b1ad-63cb-4cfa-98f4-a11ab05a9022
# ╠═d1efdbe8-1a18-4a1c-84fc-f80dc3b9272d
# ╠═e0e9cbb2-01a8-4523-a6a9-ec7cfdb351d5
# ╠═e9dd6424-1031-44fb-970c-15aa3a4a1385
# ╠═5feafefc-740c-4450-838d-9e2ac47ceb32
# ╠═7b143cf1-920d-4626-a927-bf33885ae7a9
# ╠═31ff651f-b275-46d2-97e9-206c9926dfb8
# ╠═258017ed-4bf0-4fdd-8239-2fc381e94d0c
# ╠═45d2cddd-adfe-470a-89bc-56fb5d5982dc
# ╠═d607a5e3-78f5-4308-999f-bf49858004a7
# ╠═c9254d06-75b9-47dc-8d29-495c80c9215d
# ╟─15fd5bd7-d2c4-4337-a2ad-0fc43ad2b4e0
# ╠═248c158c-41bc-478f-be12-b7e4ab43544d
# ╠═72c6bdcd-6d6b-4bf4-a759-1c4256384352
# ╠═c0634c74-9228-4d05-8146-3492aaf3a8ef
# ╟─decac4e3-d5d3-4a51-a5bb-fdb658dec971
# ╟─2bb59351-cb34-45a7-bd58-72409adf48fa
# ╟─4aa39de7-db17-42d7-924f-abb486c03eed
# ╟─d5ba2170-d687-4ef8-a7fc-e8fef48b3a2c
# ╠═b385f71f-c608-4be7-a0b6-30a39397451e
# ╟─c1480936-bc64-481d-a754-31d5a35a2d6d
# ╠═c1cd29c9-b3dc-497b-b200-e6ad08356b39
# ╟─7d6a905a-2445-4569-8fe1-8e703b635218
# ╠═4b879529-0359-462e-b032-3bb95a11ac36
# ╟─3e54ccec-d960-4cc2-a34a-7fce0964addb

max_stars_repo_path: Unbeatable Tic Tac Toe.jl

max_stars_repo_name: ctrekker/pluto-static-test

max_issues_repo_path: Unbeatable Tic Tac Toe.jl

max_issues_repo_name: ctrekker/pluto-static-test

max_forks_repo_path: Unbeatable Tic Tac Toe.jl

max_forks_repo_name: ctrekker/pluto-static-test


text:
covarianceest_insts(ext::MatSpecExt) = [
    [(d, true, ext) # complex
    for d in vcat(3, 25:25:250)] # includes compile run
    ]

insts = OrderedDict()
insts["logdet"] = (nothing, covarianceest_insts(MatLogdetCone()))
insts["sepspec"] = (nothing, covarianceest_insts(MatNegLog()))
insts["direct"] = (nothing, covarianceest_insts(MatNegLogDirect()))
return (CovarianceEstJuMP, insts)


text:
@testset "GDAL" begin

#=
Gdal.GDALAllRegister()

driver = Gdal.GDALGetDriverByName("GTiff")
srs = Gdal.OSRNewSpatialReference(C_NULL)
Gdal.OSRImportFromEPSG(srs, 4326) 			# fails if GDAL_DATA is not set correctly
Gdal.GDALDestroyDriverManager()
=#

	version = Gdal.GDALVersionInfo("--version")
	n_gdal_driver = Gdal.GDALGetDriverCount()
	n_ogr_driver = Gdal.OGRGetDriverCount()
	@info """$version
	$n_gdal_driver GDAL drivers found
	$n_ogr_driver OGR drivers found
	"""

	dataset = Gdal.create("", driver = getdriver("MEM"), width=241, height=181, nbands=1, dtype=Float64)
	crs = toWKT(importPROJ4("+proj=latlong"), true);
	crs = toWKT(importPROJ4("+proj=latlong"));
	importWKT(crs)
	Gdal.importWKT!(importWKT(crs), toWKT(importPROJ4("+proj=latlong")))
	Gdal.write!(dataset, rand(181,241), 1)
	setproj!(dataset, crs)
	setgeotransform!(dataset, [-4.016666666666667, 0.03333333333333333, 0.0, -3.01666666666, 0.0, 0.03333333333333333])
	Gdal.listcapability(dataset)

	#show(dataset)
	sprint(print, dataset)

	band = getband(dataset);
	Gdal.read(band);
	getproj(dataset)
	band = getband(dataset)
	Gdal.width(dataset)
	Gdal.width(band)
	Gdal.height(dataset)
	Gdal.height(band)
	Gdal.nlayer(dataset)
	Gdal.nraster(dataset)
	Gdal.filelist(dataset)
	Gdal.accessflag(band)
	Gdal.indexof(band)
	Gdal.pixeltype(band)
	Gdal.getcolorinterp(band)
	Gdal.getgeotransform(dataset)

	getdriver(dataset)
	getdriver(5)

	toPROJ4(importEPSG(4326))
	importWKT(crs)

	Gdal.shortname(getdriver("GTiff"));
	Gdal.longname(getdriver("GTiff"));
	Gdal.driveroptions("MEM");

	Gdal.GDALError();
	Gdal.CPLGetLastErrorNo();
	Gdal.CPLGetLastErrorMsg();
	Gdal.GDALGetColorInterpretationName(1);
	Gdal.GDALGetPaletteInterpretationName(1);
	Gdal.OCTDestroyCoordinateTransformation(C_NULL);

	gdalinfo("utmsmall.tif");
	ds_small = Gdal.read("utmsmall.tif");
	Gdal.getlayer(ds_small, 1);
	gdalinfo(ds_small, [""]);
	#gdaldem("utmsmall.tif", "hillshade", ["-q"], save="lixo.nc");
	#rm("lixo.nc")
	gdaldem(ds_small, "hillshade", ["-q"]);
	gdaltranslate(ds_small, [""]);
	gdaltranslate("utmsmall.tif", R="442000/445000/3747000/3750000");
	try		# Stupid Macports gdal does not have AAIGrid driver
	ds_tiny = gdaltranslate(ds_small, ["-of","AAIGrid","-r","cubic","-tr","1200","1200"], gdataset=true) # resample to a 5×5 ascii grid
	@test Gdal.read(ds_tiny, 1) == [128  171  127   93   83; 126  164  148  114  101;
									161  175  177  164  140; 185  206  205  172  128;
									193  205  209  181  122]

	ds_vrt = gdalbuildvrt([ds_tiny])
	@test Gdal.read(ds_vrt, 1) == [128  171  127   93   83; 126  164  148  114  101;
									161  175  177  164  140; 185  206  205  172  128;
									193  205  209  181  122]
	catch
	end

	gdalwarp(ds_small, [""]);
	ds_warped = gdalwarp("utmsmall.tif", ["-of","MEM","-t_srs","EPSG:4326"], I=0.005, gdataset=true)
	ds_warped = gdalwarp("utmsmall.tif", ["-of","MEM","-t_srs","EPSG:4326"], gdataset=true)
	@test Gdal.width(ds_warped) == 109
	@test Gdal.height(ds_warped) == 91

	ds_point = Gdal.read("point.geojson");
	ds_grid = gdalgrid(ds_point, ["-of","MEM","-outsize","3", "10","-txe","100","100.3","-tye","0","0.1"], gdataset=true);
	@test getgeotransform(ds_grid) ≈ [100.0,0.1,0.0,0.0,0.0,0.01]
	ds_point = Gdal.read("point.geojson");
	layer = Gdal.getlayer(ds_point, 0);
	feature = Gdal.nextfeature(layer)
	@test Gdal.asdouble(feature, 0) ≈ 2.0
	@test Gdal.asstring(feature, 1) == "point-a"
	sprint(print, feature)
	sprint(print, Gdal.getfielddefn(feature, 1))
	sprint(print, Gdal.getgeomdefn(feature, 0))
	sprint(print, ds_point)

	ds_rasterize = gdalrasterize(ds_point, ["-of","MEM","-tr","0.05","0.05"], gdataset=true)
	@test getgeotransform(ds_rasterize) ≈ [99.975,0.05,0.0,0.1143,0.0,-0.05]

	ds_csv = gdalvectortranslate(ds_point, ["-f","CSV","-lco", "GEOMETRY=AS_XY"], dest = "point.csv");
	#=
	@test replace(read("point.csv", String), "\r" => "") == """
	X,Y,FID,pointname
	100,0,2,point-a
	100.2785,0.0893,3,point-b
	100,0,0,a
	100.2785,0.0893,3,b
	"""
	=#

	Gdal.createmultipoint([(1251243.7361610543, 598078.7958668759), (1250318.7031934808, 606404.0925750365)]);

	@test GMT.R_inc_to_gd([0.01], " -R-11/-1/33/45")[1] == "-txe";

	dataset = Gdal.create(Gdal.getdriver("MEMORY"))
	layer = Gdal.createlayer(name = "point_out", dataset = dataset, geom = Gdal.wkbPoint)
	Gdal.addfielddefn!(layer, "Name", Gdal.OFTString, nwidth = 32)
	featuredefn = Gdal.layerdefn(layer)
	@test Gdal.getname(featuredefn) == "point_out"
	@test Gdal.nfeature(layer) == 0
	Gdal.createfeature(layer) do feature
		Gdal.setfield!(feature, Gdal.findfieldindex(feature, "Name"), "myname")
		Gdal.setgeom!(feature, Gdal.createpoint(100.123, 0.123))
	end
	@test Gdal.nfeature(layer) == 1

	ds_src = Gdal.read("utmsmall.tif")
	isa(ds_src, Gdal.RasterDataset{UInt8})
	Gdal.write(ds_src, "/vsimem/utmsmall.tif")
	ds_copy = Gdal.read("/vsimem/utmsmall.tif")
	@test Gdal.read(ds_src) == Gdal.read(ds_copy)
	Gdal.metadata(ds_src)
	Gdal.GDALGetDescription(ds_src.ptr)
	imshow(ds_src, Vd=dbg2)

	rb = Gdal.getband(ds_src, 1)
	@test Gdal.getnodatavalue(rb) === nothing
	Gdal.setnodatavalue!(rb, -100)
	@test Gdal.getnodatavalue(rb) ≈ -100
	Gdal.getcolorinterp(rb)

	ct = Gdal.createcolortable(1)		# GDAL.GPI_RGB = 1
	Gdal.createcolorramp!(ct, 128, Gdal.GDALColorEntry(0,0,0,0), 255, Gdal.GDALColorEntry(0,0,255,0))

	line = Gdal.createlinestring()
	Gdal.addpoint!(line, 1116651.439379124,  637392.6969887456)
	Gdal.OGR_G_SetPoints(line.ptr, 3, [1.,2,3], sizeof(Float64), [4.,5,6], sizeof(Float64), [7.,8,9], sizeof(Float64))
	xx = Gdal.getpoint(line, 1)
	@test xx == (2.0, 5.0, 8.0)

	G = GMT.peaks()
	gdalshade(G, C=makecpt(T=(-7,8,1)), zfactor=2);
	GMT.gammacorrection(mat2img(rand(UInt8, 3,4)), 1.1)
	ds = gmt2gd(G)
	G  = gd2gmt(ds)
	G = gd2gmt("utmsmall.tif");
	ds = gmt2gd(G)
	gdalinfo(ds);
	D = mat2ds([0 0; 10 0; 10 10; 11 10], geom=4);
	Dv = gd2gmt(gmt2gd(D));
	@test D == Dv

	GMT.get_FillValue("_FillValue=9999")
	GMT.helper_gmt2gd_xyz(GMT.GMTdataset(rand(Float32,3,3), String[], "", String[], "", "", 0),3)
	@test GMT.gd2gmt_helper_scalefac([1  1; 1 1], 2, 0, false, 0) == [2 2; 2 2]
	@test GMT.gd2gmt_helper_scalefac([1. 1; 1 1], 2, 0, false, 0) == [2. 2.; 2. 2.]

	# Test that we recover the original.
	G = mat2grid(reshape(collect(1.0:12), 3, 4));
	Gv = gd2gmt(gmt2gd(G), layout="BCB");
	@test G == Gv

	I = grdcut("utmsmall.tif", R="442000/445000/3747000/3750000", img=1);
	grdcut("utmsmall.tif", R="442000/445000/3747000/3750000", img=1, save="lixo.tif");
	grdcut("utmsmall.tif", R="442000/445000/3747000/3750000", img=1, save="lixo.tif");

	GMT.Gdal.getmetadata("rgbsmall.tif");
	I = Gdal.dither("rgbsmall.tif");
	gmt2gd(I);
	Gdal.dither("rgbsmall.tif", save="lixo.tif");
	Gdal.gdal_opts2vec("aa bb 'vv pp' aa \"ad uuu\"");

	Gdal.GDALGetDataTypeByName("GTiff");
	Gdal.IFieldDefnView(C_NULL);
	Gdal.IGeomFieldDefnView(C_NULL);
	Gdal.GeomFieldDefn(C_NULL);
	Gdal.RasterBand(C_NULL);
	Gdal.destroy(Gdal.Driver(C_NULL));
	Gdal.destroy(Gdal.Feature(C_NULL));
	Gdal.destroy(Gdal.CoordTransform(C_NULL));
	Gdal.destroy(Gdal.FeatureDefn(C_NULL));
	Gdal.destroy(Gdal.IFieldDefnView(C_NULL));
	Gdal.destroy(Gdal.IGeomFieldDefnView(C_NULL));

	ds = gmt2gd(mat2ds([-8. 37.0; -8.1 37.5; -8.5 38.0]))
	Gdal.getx(Gdal.getgeom(Gdal.unsafe_getfeature(Gdal.getlayer(ds, 0),0)),1)
	Gdal.gety(Gdal.getgeom(Gdal.unsafe_getfeature(Gdal.getlayer(ds, 0),0)),1)
	Gdal.getz(Gdal.getgeom(Gdal.unsafe_getfeature(Gdal.getlayer(ds, 0),0)),1)
	Gdal.buffer(Gdal.getgeom(Gdal.unsafe_getfeature(Gdal.getlayer(ds, 0),0)), 0.2)
	Gdal.GDALGetSpatialRef(ds.ptr)
	Gdal.GDALDatasetGetLayerByName(ds.ptr, "0")
	Gdal.GDALGetPaletteInterpretation(ds.ptr)
	GMT.Gdal.getmetadata(ds)

	#Gdal.identifydriver("lixo.gmt")
	D = mat2ds([-8. 37.0; -8.1 37.5; -8.5 38.0], proj="+proj=longlat");
	ds = gmt2gd(D[1])
	ds = gmt2gd(D, geometry="Polygon")
	gmt2gd(mat2ds([0 10; 1 11; 0 10], x=[0, 1, 2], multi=true), geometry="line")
	gmt2gd(mat2ds([0 10; 1 11; 0 10; 0 10], x=[0, 1, 2, 0], multi=true), geometry="polyg")
	ogr2ogr(D, dest="lixo1.gmt")
	gmt2gd(D, save="lixo2.gmt")
	ds = gmt2gd(D)
	ds2=ogr2ogr(ds, ["-t_srs", "+proj=utm +zone=29", "-overwrite"])
	gd2gmt(ds2)
	ogr2ogr(D, "-t_srs '+proj=utm +zone=29' -overwrite")	# Returns a GMT datset directly
	ogr2ogr(D[1], "-t_srs '+proj=utm +zone=29' -overwrite")

	D1 = mat2ds([0.0 0.0; 1.0 1.0; 1.0 0.0; 0.0 0.0]);
	gmt2gd(D1);		gmt2gd(D1, geometry="line");	gmt2gd(D1, geometry="point")
	D2 = mat2ds([0.0 0.0 1.; 1.0 1.0 2.; 1.0 0.0 3.; 0.0 0.0 1.]);
	gmt2gd(D2);		gmt2gd(D2, geometry="line");	gmt2gd(D2, geometry="point")

	wkt = "POLYGON ((1179091. 712782.,1161053. 667456.,1214705. 641092.,1228580. 682719.,1218405. 721108.,1179091. 712782.))"
	@test Gdal.getgeomtype(Gdal.forceto(Gdal.fromWKT(wkt), Gdal.wkbMultiPolygon)) == Gdal.wkbMultiPolygon
	wkt = "POINT (1198054.34 648493.09)";
	bf = Gdal.buffer(Gdal.fromWKT(wkt), 500)
	@test Gdal.getgeomtype(bf) == Gdal.wkbPolygon
	show(bf)

	@testset "Calculate the Area of a Geometry" begin
		wkt = "POLYGON ((1162440. 672081., 1162440. 647105., 1195279. 647105., 1195279. 672081., 1162440. 672081.))"
		poly = Gdal.fromWKT(wkt)
		@test Gdal.geomarea(poly) ≈ 8.20186864e8
	end

	@testset "Calculate the Length of a Geometry" begin
		wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
		line = Gdal.fromWKT(wkt)
		@test Gdal.geomlength(line) ≈ 76121.94397805972
	end

	mp = Gdal.createmultipolygon_noholes(Vector{Tuple{Float64,Float64}}[
		[(1204067., 634617.), (1204067., 620742.), (1215167., 620742.), (1215167., 634617.), (1204067., 634617.)],
		[(1179553., 647105.), (1179553., 626292.), (1194354., 626292.), (1194354., 647105.), (1179553., 647105.)] ])
	Gdal.wrapgeom(mp)

	gdalinfo("poly_spatialite.sqlite")
	Gdal.read("poly_spatialite.sqlite")
	Gdal.inspect("SELECT HEX(GeomFromText('POINT(10 20)'))", "poly_spatialite.sqlite")

	I1 = mat2img(reshape(collect(UInt8(1):UInt8(20)), 4, 5))	#  layout = TCBa
	I2 = mat2img(reshape(collect(UInt8(11):UInt8(30)), 4, 5))
	GMT.blendimg!(I1, I2)
	I2.layout = "BCBa"
	GMT.blendimg!(I1, I2)
	I2.layout = "BRBa"
	GMT.blendimg!(I1, I2)
	I2.layout = "TRBa"
	GMT.blendimg!(I1, I2)
	I1 = mat2img(reshape(collect(UInt8(1):UInt8(60)), 4, 5, 3))
	I1.layout = "TRPa"
	I2.image[end] = UInt8(200)		# to make tests visit another if branch
	GMT.blendimg!(I1, I2)
end


text:
using ConicBenchmarkUtilities
using ECOS
using SCS
using MathProgBase

# julia testexpsolvers.jl ECOS ../instancedata/minlplib2/syn/*
# julia testexpsolvers.jl SCS ../instancedata/minlplib2/syn/*

# see if the solver can solve the continuous relaxations

if ARGS[1] == "ECOS"
    solver = ECOSSolver()
elseif ARGS[1] == "SCS"
    solver = SCSSolver()
else
    error("Unrecognized solver $(ARGS[1])")
end

solved = 0
suboptimal = 0

results = []

for F in ARGS[2:end]
    dat = readcbfdata(F)

    c, A, b, con_cones, var_cones, vartypes, sense, objoffset = cbftompb(dat)

    @assert dat.sense == :Min
    @assert dat.objoffset == 0.0

    m = MathProgBase.ConicModel(solver)
    MathProgBase.loadproblem!(m, c, A, b, con_cones, var_cones)
    #MathProgBase.setvartype!(m, vartypes)
    MathProgBase.optimize!(m)
    status = MathProgBase.status(m)
    if status == :Optimal
        solved += 1
    elseif status == :Suboptimal
        suboptimal += 1
    end
    push!(results,(F,status))

end

for (name,status) in results
    println("$(basename(name)) $status")
end

println("SOLVED $solved SUBOPTIMAL $suboptimal TOTAL $(length(ARGS)-1)")

max_stars_repo_name: mlubin/PajaritoSupplement

max_issues_repo_name: mlubin/PajaritoSupplement

max_forks_repo_name: mlubin/PajaritoSupplement


text:
dfg = testDFGAPI{NoSolverParams}()

#add types for softtypes
struct TestInferenceVariable1 <: InferenceVariable end
struct TestInferenceVariable2 <: InferenceVariable end

v1 = DFGVariable(:a, TestInferenceVariable1())
v2 = DFGVariable(:b, TestInferenceVariable2())
f1 = DFGFactor{Int, :Symbol}(:f1)

#add tags for filters
append!(v1.tags, [:VARIABLE, :POSE])
append!(v2.tags, [:VARIABLE, :LANDMARK])
append!(f1.tags, [:FACTOR])

st1 = TestInferenceVariable1()
st2 = TestInferenceVariable2()

v1.solverDataDict[:default].softtype = deepcopy(st1)
v2.solverDataDict[:default].softtype = deepcopy(st2)

# NOTE: Just for testing
# Override compareAllSpecial for testing our DFGFactor{Int, :Symbol}
# as we don't have IIF which overloads with the proper compareSpecial(::CommonConvWrapper)
import DistributedFactorGraphs: compareAllSpecial
function compareAllSpecial(A::T1, B::T2;
                    skip=Symbol[], show::Bool=true) where {T1 <: Int, T2 <: Int}
  #
  T1 != T2 && return false
  return compareAll(A, B, skip=skip, show=show)
end

# @testset "Creating Graphs" begin
global dfg,v1,v2,f1
addVariable!(dfg, v1)
addVariable!(dfg, v2)
addFactor!(dfg, [v1, v2], f1)
@test_throws Exception addFactor!(dfg, DFGFactor{Int, :Symbol}("f2"), [v1, DFGVariable("Nope")])
# end

@testset "Adding Removing Nodes" begin
    dfg2 = testDFGAPI{NoSolverParams}()
    v1 = DFGVariable(:a, TestInferenceVariable1())
    v2 = DFGVariable(:b, TestInferenceVariable1())
    v3 = DFGVariable(:c, TestInferenceVariable1())
    f1 = DFGFactor{Int, :Symbol}(:f1)
    f2 = DFGFactor{Int, :Symbol}(:f2)
    # @testset "Creating Graphs" begin
    @test addVariable!(dfg2, v1)
    @test addVariable!(dfg2, v2)
    @test_throws ErrorException updateVariable!(dfg2, v3)
    @test addVariable!(dfg2, v3)
    @test_throws ErrorException addVariable!(dfg2, v3)
    @test addFactor!(dfg2, [v1, v2], f1)
    @test_throws ErrorException addFactor!(dfg2, [v1, v2], f1)
    @test_throws ErrorException updateFactor!(dfg2, f2)
    @test addFactor!(dfg2, [:b, :c], f2)
    @test deleteVariable!(dfg2, v3) == v3
    @test symdiff(ls(dfg2),[:a,:b]) == []
    @test deleteFactor!(dfg2, f2) == f2
    @test lsf(dfg2) == [:f1]
end

@testset "Listing Nodes" begin
    global dfg,v1,v2,f1
    @test length(ls(dfg)) == 2
    @test length(lsf(dfg)) == 1
    @test symdiff([:a, :b], getVariableIds(dfg)) == []
    @test getFactorIds(dfg) == [:f1]
    #
    @test lsf(dfg, :a) == [f1.label]
    # Tags
    @test ls(dfg, tags=[:POSE]) == [:a]
    @test symdiff(ls(dfg, tags=[:POSE, :LANDMARK]), ls(dfg, tags=[:VARIABLE])) == []
    # Regexes
    @test ls(dfg, r"a") == [v1.label]
    @test lsf(dfg, r"f*") == [f1.label]
    # Accessors
    @test getAddHistory(dfg) == [:a, :b] #, :f1
    @test getDescription(dfg) != nothing
    @test getLabelDict(dfg) != nothing
    # Existence
    @test exists(dfg, :a) == true
    @test exists(dfg, v1) == true
    @test exists(dfg, :nope) == false
    # Sorting of results
    # TODO - this function needs to be cleaned up
    unsorted = [:x1_3;:x1_6;:l1;:april1] #this will not work for :x1x2f1
    @test sortDFG(unsorted) == sortVarNested(unsorted)
    @test_skip sortDFG([:x1x2f1, :x1l1f1]) == [:x1l1f1, :x1x2f1]
end

# Gets
@testset "Gets, Sets, and Accessors" begin
    global dfg,v1,v2,f1
    @test getVariable(dfg, v1.label) == v1
    @test getVariable(dfg, v2.label) != v1
    @test getFactor(dfg, f1.label) == f1
    f2 = deepcopy(f1)
    f2.label = :something
    @test f2 != f1
    @test_throws Exception getVariable(dfg, :nope)
    @test_throws Exception getVariable(dfg, "nope")
    @test_throws Exception getFactor(dfg, :nope)
    @test_throws Exception getFactor(dfg, "nope")

    # Sets
    v1Prime = deepcopy(v1)
    #updateVariable! returns the variable updated, so should be equal
    @test updateVariable!(dfg, v1Prime) == v1
    f1Prime = deepcopy(f1)
    #updateFactor! returns the factor updated, so should be equal
    @test updateFactor!(dfg, f1Prime) == f1

    # Accessors
    @test label(v1) == v1.label
    @test tags(v1) == v1.tags
    @test timestamp(v1) == v1.timestamp
    @test estimates(v1) == v1.estimateDict
    @test estimate(v1, :notfound) == nothing
    @test solverData(v1) === v1.solverDataDict[:default]
    @test getData(v1) === v1.solverDataDict[:default]
    @test solverData(v1, :default) === v1.solverDataDict[:default]
    @test solverDataDict(v1) == v1.solverDataDict
    @test internalId(v1) == v1._internalId

    @test softtype(v1) == Symbol(typeof(st1))
    @test softtype(v2) == Symbol(typeof(st2))
    @test getSofttype(v1) == st1

    @test label(f1) == f1.label
    @test tags(f1) == f1.tags
    @test solverData(f1) == f1.data
    # Deprecated functions
    @test data(f1) == f1.data
    @test getData(f1) == f1.data
    # Internal function
    @test internalId(f1) == f1._internalId

    @test getSolverParams(dfg) != nothing
    @test setSolverParams(dfg, getSolverParams(dfg)) == getSolverParams(dfg)

    #solver data is initialized
    @test !isInitialized(dfg, :a)
    @test !isInitialized(v2)

    @test !isInitialized(v2, key=:second)

    # Session, robot, and user small data tests
    smallUserData = Dict{Symbol, String}(:a => "42", :b => "Hello")
    smallRobotData = Dict{Symbol, String}(:a => "43", :b => "Hello")
    smallSessionData = Dict{Symbol, String}(:a => "44", :b => "Hello")
    setUserData(dfg, deepcopy(smallUserData))
    setRobotData(dfg, deepcopy(smallRobotData))
    setSessionData(dfg, deepcopy(smallSessionData))
    @test getUserData(dfg) == smallUserData
    @test getRobotData(dfg) == smallRobotData
    @test getSessionData(dfg) == smallSessionData

end

@testset "BigData" begin
    oid = zeros(UInt8,12); oid[12] = 0x01
    de1 = MongodbBigDataEntry(:key1, NTuple{12,UInt8}(oid))

    oid = zeros(UInt8,12); oid[12] = 0x02
    de2 = MongodbBigDataEntry(:key2, NTuple{12,UInt8}(oid))

    oid = zeros(UInt8,12); oid[12] = 0x03
    de2_update = MongodbBigDataEntry(:key2, NTuple{12,UInt8}(oid))

    #add
    v1 = getVariable(dfg, :a)
    @test addBigDataEntry!(v1, de1)
    @test addBigDataEntry!(dfg, :a, de2)
    @test addBigDataEntry!(v1, de1)

    #get
    @test deepcopy(de1) == getBigDataEntry(v1, :key1)
    @test deepcopy(de2) == getBigDataEntry(dfg, :a, :key2)
    @test_throws Any getBigDataEntry(v2, :key1)
    @test_throws Any getBigDataEntry(dfg, :b, :key1)

    #update
    @test updateBigDataEntry!(dfg, :a, de2_update)
    @test deepcopy(de2_update) == getBigDataEntry(dfg, :a, :key2)
    @test !updateBigDataEntry!(dfg, :b, de2_update)

    #list
    entries = getBigDataEntries(dfg, :a)
    @test length(entries) == 2
    @test symdiff(map(e->e.key, entries), [:key1, :key2]) == Symbol[]
    @test length(getBigDataEntries(dfg, :b)) == 0

    @test symdiff(getBigDataKeys(dfg, :a), [:key1, :key2]) == Symbol[]
    @test getBigDataKeys(dfg, :b) == Symbol[]

    #delete
    @test deepcopy(de1) == deleteBigDataEntry!(v1, :key1)
    @test getBigDataKeys(v1) == Symbol[:key2]
    #delete from dfg
    @test deepcopy(de2_update) == deleteBigDataEntry!(dfg, :a, :key2)
    @test getBigDataKeys(v1) == Symbol[]
end

@testset "Updating Nodes" begin
    global dfg
    #get the variable
    var = getVariable(dfg, :a)
    #make a copy and simulate external changes
    newvar = deepcopy(var)
    estimates(newvar)[:default] = MeanMaxPPE(:default, [150.0], [100.0], [50.0])
    #update
    mergeUpdateVariableSolverData!(dfg, newvar)

    #For now spot check
    @test solverDataDict(newvar) == solverDataDict(var)
    @test estimates(newvar) == estimates(var)

    # Delete :default and replace to see if new ones can be added
    delete!(estimates(newvar), :default)
    estimates(newvar)[:second] = MeanMaxPPE(:second, [15.0], [10.0], [5.0])

    # Persist to the original variable.
    mergeUpdateVariableSolverData!(dfg, newvar)
    # At this point newvar will have only :second, and var should have both (it is the reference)
    @test symdiff(collect(keys(estimates(var))), [:default, :second]) == Symbol[]
    @test symdiff(collect(keys(estimates(newvar))), [:second]) == Symbol[]
    # Get the source too.
    @test symdiff(collect(keys(estimates(getVariable(dfg, :a)))), [:default, :second]) == Symbol[]
end

# Connectivity test
@testset "Connectivity Test" begin
    global dfg,v1,v2,f1
    @test isFullyConnected(dfg) == true
    @test hasOrphans(dfg) == false
    addVariable!(dfg, DFGVariable(:orphan))
    @test isFullyConnected(dfg) == false
    @test hasOrphans(dfg) == true
end

# Adjacency matrices
@testset "Adjacency Matrices" begin
    global dfg,v1,v2,f1
    # Normal
    adjMat = getAdjacencyMatrix(dfg)
    @test size(adjMat) == (2,4)
    @test symdiff(adjMat[1, :], [nothing, :a, :b, :orphan]) == Symbol[]
    @test symdiff(adjMat[2, :], [:f1, :f1, :f1, nothing]) == Symbol[]
    #sparse
    adjMat, v_ll, f_ll = getAdjacencyMatrixSparse(dfg)
    @test size(adjMat) == (1,3)

    # Checking the elements of adjacency, its not sorted so need indexing function
    indexOf = (arr, el1) -> findfirst(el2->el2==el1, arr)
    @test adjMat[1, indexOf(v_ll, :orphan)] == 0
    @test adjMat[1, indexOf(v_ll, :a)] == 1
    @test adjMat[1, indexOf(v_ll, :b)] == 1
    @test symdiff(v_ll, [:a, :b, :orphan]) == Symbol[]
    @test symdiff(f_ll, [:f1, :f1, :f1]) == Symbol[]
end

# Deletions
@testset "Deletions" begin
    deleteFactor!(dfg, :f1)
    @test getFactorIds(dfg) == []
    deleteVariable!(dfg, :b)
    @test symdiff([:a, :orphan], getVariableIds(dfg)) == []
    #delete last also for the LightGraphs implementation coverage
    deleteVariable!(dfg, :orphan)
    @test symdiff([:a], getVariableIds(dfg)) == []
    deleteVariable!(dfg, :a)
    @test getVariableIds(dfg) == []
end

# Now make a complex graph for connectivity tests
numNodes = 10
dfg = testDFGAPI{NoSolverParams}()
verts = map(n -> DFGVariable(Symbol("x$n"), TestInferenceVariable1()), 1:numNodes)
#change ready and backendset for x7,x8 for improved tests on x7x8f1
verts[7].ready = 1
verts[8].backendset = 1

# Can't change the softtypes now.
# #force softytypes to first 2 vertices.
# verts[1].solverDataDict[:default].softtype = deepcopy(st1)
# verts[2].solverDataDict[:default].softtype = deepcopy(st2)

map(v -> addVariable!(dfg, v), verts)
map(n -> addFactor!(dfg, [verts[n], verts[n+1]], DFGFactor{Int, :Symbol}(Symbol("x$(n)x$(n+1)f1"))), 1:(numNodes-1))

@testset "Getting Neighbors" begin
    global dfg,verts
    # Trivial test to validate that intersect([], []) returns order of first parameter
    @test intersect([:x3, :x2, :x1], [:x1, :x2]) == [:x2, :x1]
    # Get neighbors tests
    @test getNeighbors(dfg, verts[1]) == [:x1x2f1]
    neighbors = getNeighbors(dfg, getFactor(dfg, :x1x2f1))
    @test neighbors == [:x1, :x2]
    # Testing aliases
    @test getNeighbors(dfg, getFactor(dfg, :x1x2f1)) == ls(dfg, getFactor(dfg, :x1x2f1))
    @test getNeighbors(dfg, :x1x2f1) == ls(dfg, :x1x2f1)

    # ready and backendset
    @test getNeighbors(dfg, :x5, ready=1) == Symbol[]
    @test getNeighbors(dfg, :x5, ready=0) == [:x4x5f1,:x5x6f1]
    @test getNeighbors(dfg, :x5, backendset=1) == Symbol[]
    @test getNeighbors(dfg, :x5, backendset=0) == [:x4x5f1,:x5x6f1]
    @test getNeighbors(dfg, :x7x8f1, ready=0) == [:x8]
    @test getNeighbors(dfg, :x7x8f1, backendset=0) == [:x7]
    @test getNeighbors(dfg, :x7x8f1, ready=1) == [:x7]
    @test getNeighbors(dfg, :x7x8f1, backendset=1) == [:x8]
    @test getNeighbors(dfg, verts[1], ready=0) == [:x1x2f1]
    @test getNeighbors(dfg, verts[1], ready=1) == Symbol[]
    @test getNeighbors(dfg, verts[1], backendset=0) == [:x1x2f1]
    @test getNeighbors(dfg, verts[1], backendset=1) == Symbol[]

end

@testset "Getting Subgraphs" begin
    # Subgraphs
    dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 2)
    # Only returns x1 and x2
    @test symdiff([:x1, :x1x2f1, :x2], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []
    # Test include orphan factorsVoid
    dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 1, true)
    @test symdiff([:x1, :x1x2f1], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []
    # Test adding to the dfg
    dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 2, true, dfgSubgraph)
    @test symdiff([:x1, :x1x2f1, :x2], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []
    #
    dfgSubgraph = getSubgraph(dfg,[:x1, :x2, :x1x2f1])
    # Only returns x1 and x2
    @test symdiff([:x1, :x1x2f1, :x2], [ls(dfgSubgraph)..., lsf(dfgSubgraph)...]) == []

    # DFG issue #95 - confirming that getSubgraphAroundNode retains order
    # REF: https://github.com/JuliaRobotics/DistributedFactorGraphs.jl/issues/95
    for fId in getVariableIds(dfg)
        # Get a subgraph of this and it's related factors+variables
        dfgSubgraph = getSubgraphAroundNode(dfg, verts[1], 2)
        # For each factor check that the order the copied graph == original
        for fact in getFactors(dfgSubgraph)
            @test fact._variableOrderSymbols == getFactor(dfg, fact.label)._variableOrderSymbols
        end
    end
end

@testset "Summaries and Summary Graphs" begin
    factorFields = fieldnames(DFGFactorSummary)
    variableFields = fieldnames(DFGVariableSummary)

    summary = getSummary(dfg)
    @test symdiff(collect(keys(summary.variables)), ls(dfg)) == Symbol[]
    @test symdiff(collect(keys(summary.factors)), lsf(dfg)) == Symbol[]

    summaryGraph = getSummaryGraph(dfg)
    @test symdiff(ls(summaryGraph), ls(dfg)) == Symbol[]
    @test symdiff(lsf(summaryGraph), lsf(dfg)) == Symbol[]
    # Check all fields are equal for all variables
    for v in ls(summaryGraph)
        for field in variableFields
            if field != :softtypename
                @test getfield(getVariable(dfg, v), field) == getfield(getVariable(summaryGraph, v), field)
            else
                @test softtype(getVariable(dfg, v)) == softtype(getVariable(summaryGraph, v))
            end
        end
    end
    for f in lsf(summaryGraph)
        for field in factorFields
            @test getfield(getFactor(dfg, f), field) == getfield(getFactor(summaryGraph, f), field)
        end
    end
end

@testset "Producing Dot Files" begin
    # create a simpler graph for dot testing
    dotdfg = testDFGAPI{NoSolverParams}()
    v1 = DFGVariable(:a, TestInferenceVariable1())
    v2 = DFGVariable(:b, TestInferenceVariable1())
    f1 = DFGFactor{Int, :Symbol}(:f1)
    addVariable!(dotdfg, v1)
    addVariable!(dotdfg, v2)
    addFactor!(dotdfg, [v1, v2], f1)

    @test toDot(dotdfg) == "graph graphname {\n2 [\"label\"=\"b\",\"shape\"=\"ellipse\",\"fillcolor\"=\"red\",\"color\"=\"red\"]\n2 -- 3\n3 [\"label\"=\"f1\",\"shape\"=\"box\",\"fillcolor\"=\"blue\",\"color\"=\"blue\"]\n1 [\"label\"=\"a\",\"shape\"=\"ellipse\",\"fillcolor\"=\"red\",\"color\"=\"red\"]\n1 -- 3\n}\n"
    @test toDotFile(dotdfg, "something.dot") == nothing
    Base.rm("something.dot")

end


text:
export @get, @set, @reset, @molecule, @lookup, @freezecore

using Formatting
using PrettyTables

"""
    Fermi.Options

Module to manage options in Fermi.

# Functions

    Fermi.Options.set(option, value)     Set an <option> to a given <value>
    Fermi.Options.get(option)            Return the current value of an <option>
    Fermi.Options.reset()                Reset all options to default values
    Fermi.Options.reset(option)          Reset a specific <option> to its default value
    Fermi.Options.molecule(molstring)    Read in a String for the `molstring` option

Alternatively, at global scope, one can use the corresponding macros that create shortcuts
for the commands above

# Macros

    @set <option> <value>   Set an <option> to a given <value>
    @get <option>           Return the current value of an <option>
    @reset                  Reset all options to default values
    @reset <option>         Reset a specific <option> to its default value
    @molecule               Read in a String for the `molstring` option
    @lookup <string>        Search for keywords containing <string>
"""
module Options
export FermiException

"""
    Fermi.Options.Default

Dictionary containing default options for Fermi. Any information not given
explicitly to Methods is obtained from here.
"""
const Default = Dict{String,Union{Float64,Int,String,Bool}}(
                                  "molstring" => """
                                  O        1.2091536548      1.7664118189     -0.0171613972
                                  H        2.1984800075      1.7977100627      0.0121161719
                                  H        0.9197881882      2.4580185570      0.6297938832
                                  """,
                                  "printstyle" => "repl",
                                  "output" => "fermi.out",
                                  "basis" => "sto-3g",
                                  "jkfit" => "auto",
                                  "rifit" => "auto",
                                  "charge" => 0,
                                  "multiplicity" => 1,
                                  "unit" => "angstrom",
                                  "reference" => "rhf",
                                  "scf_max_iter" => 50,
                                  "scf_max_rms" => 10^-9,
                                  "scf_e_conv" => 10^-10,
                                  "scf_alg" => 1,
                                  "mp2_alg" => 1,
                                  "df" => false,
                                  "oda" => true,
                                  "oda_cutoff" => 1E-1,
                                  "oda_shutoff" => 20,
                                  "scf_guess" => "gwh",
                                  "quiet" => true,
                                  "precision" => "double",
                                  "cc_alg" => 1,
                                  "ci_alg" => "aci",
                                  "det_size" => 64,
                                  "σ"      => 0.001,
                                  "γ"      => 1.0,
                                  "ζ"      => 0.95,
                                  "ζsize"  => 0,
                                  "d_conv" => 8,
                                  "cc_max_iter" => 50,
                                  "cc_max_rms" => 10^-10,
                                  "cc_e_conv" => 10^-10,
                                  "bcc_max_t1" => 1^-7,
                                  "drop_occ" => 0,
                                  "drop_vir" => 0,
                                  "diis" => true,
                                  "cc_diis" => true,
                                  "ndiis" => 8,
                                  "cc_ndiis" => 3,
                                  "cc_diis_relax"=>3,
                                  "bcc_tol" => 1.0e-6,
                                  "diis_prec" => "single",
                                  "diis_start" => 3,
                                  "cc_damp_ratio" => 0.0,
                                  "pt_alg" => 1,
                                  "num_frozen" => 0,
                                  "aci_print_screen" => 0,
                                  "cas_frozen" => 0,
                                  "cas_active" => -1,
                                  "cas_cutoff" => 10^-9,
                                  "cas_nroot" => 1,
                                  "min_matrix_elem" => 10^-9,
                                  "precision_override" => false,
                                  "tblis" => true,
                                  "eri_cutoff" => 10^-12,
                                  "findif_disp_size" => 0.005
                                 )
"""
    Fermi.Options.Current

Dictionary containing user defined options for Fermi. Unspecified options are obtained from
Fermi.Default.
"""
Current = Dict{String,Union{Float64,Int,String,Bool}}()

"""
    Fermi.Options.get(key)

Return the current options set for the given `key`.
See also `@get`
"""
function get(key::String)

    key = lowercase(key)
    if haskey(Current, key)
        return Current[key]
    elseif haskey(Default, key)
        return Default[key]
    else
        throw(FermiException(key*" is not a valid option."))
    end
end

"""
    Fermi.Options.set(key, val)

Set the option `key` to the value `val`.
See also `@set`
"""
function set(key::String, val::Union{String, Bool, Float64, Int})
    # Make all string lowercase
    key = lowercase(key)
    if key != "molstring" && val isa String
        val = lowercase(val)
    end
    if !(haskey(Default, key))
        throw(FermiException(key*" is not a valid option."))
    end

    # Check if the type of the variable passed matches the type
    # saved in the Default dictionary
    dtype = typeof(Default[key])
    if !isa(val, dtype)
        # If it does not match...

        # If the expected type is not a String we try to make a conversion
        if dtype !== String
            try newval = dtype(val)
                # Recall the function with the new variable with adjusted type
                return set(key, newval)
            catch InexactError
                # If the conversion is not possible
                throw(FermiException("Impossible to convert data type for $key. Expected: $dtype, Got: $(typeof(val))"))
            end
        end
        # Otherwise, throw an error
        throw(FermiException("Invalid data type for $key. Expected: $dtype, Got: $(typeof(val))"))
    else
        Current[key] = val
    end
    get(key)
end

"""
    Fermi.Options.reset(key...)

Reset all given `key` to the default values.
See also `@reset`
"""
function reset(key::String="all")
    if key == "all"
        for k in keys(Current)
            delete!(Current, k)
        end
    elseif haskey(Current,key)
        delete!(Current, key)
    elseif !haskey(Default, key)
        throw(FermiException(key*" is not a valid option."))
    end
end

"""
    Fermi.Options.FermiException

Error flag used when Fermi encounters an error. Use to give constructive feedback.
"""
struct FermiException <: Exception
    msg::String
end
Base.showerror(io::IO, e::FermiException) = print(io, "Fermi found an Error: ", e.msg)
end #module

# Aux function to convert Symbols to Strings
function symbol_to_string(S)
    S = repr(S)
    return String(strip(filter(c->!occursin(c," {}():"),S)))
end

"""
    Fermi.@get

Returns the current value of an option.

# Examples

```
julia> @get basis
"sto-3g"      # Default

julia> @set basis cc-pVDZ
julia> @get basis
"cc-pVDZ"
```
Without arguments, it returns all user defined keywords
```
julia> @set {
    basis cc-pVDZ
    scf_max_rms 10^-8
    diis true
}
julia> @get
┌─────────────┬───────────────┐
│     Keyword │ Current Value │
├─────────────┼───────────────┤
│ scf_max_rms │   1.00000e-08 │
│        diis │          true │
│       basis │       cc-pvdz │
└─────────────┴───────────────┘
```
"""
macro get(opt)
    A = symbol_to_string(opt)
    quote
        Fermi.Options.get($A)
    end |> esc
end

macro get()
    quote
        Keys = collect(keys(Fermi.Options.Current))
        l = length(Keys)
        if l == 0
            println("No user defined keywords found.")
        else
            data = Array{Union{Nothing,String},2}(nothing,l,2)
            for i in eachindex(Keys)
                k = Keys[i]
                data[i,1] = k
                val = Fermi.Options.get(k)
                if val isa Bool || val isa Int
                    data[i,2] = "$val"
                elseif val isa String
                    # If multiple lines (e.g. molstring) display only the first line
                    if occursin("\n", val)
                        lines = split(val, "\n")
                        data[i,2] = lines[1] * " [+$(length(lines)-1) lines]"
                    else
                        data[i,2] = "$val"
                    end
                else
                    data[i,2] = format("{:5.5e}", val)
                end
            end
            pretty_table(data; header=["Keyword", "Current Value"])
        end
    end
end

"""
    Fermi.@set

Set options for Fermi computations. It saves the options into Fermi.CurrentOptions.

*Usage:*  @set A B

A is set to B. By default A is taken as a string.

# Examples

```
@set basis cc-pVDZ
@set cc_max_iter 100
@set e_conv 10^-9
```
One can also use the block syntax
```
@set {
    basis cc-pVDZ
    cc_max_iter 100
    e_conv 10^-9
}
```
Note for basis set: Having * at the end of a line is considered an incomplete expression by Julia. Thus, for basis such as 6-31g*
you should use quotes
```
@set basis "6-31g"
```
"""
macro set(A::Symbol, B::Symbol)
    clean_up(s) = String(filter(c->!occursin(c," ():"),s))
    key = clean_up(repr(A))
    val = clean_up(repr(B))
    return quote
        Fermi.Options.set($key,$val)
    end |> esc
end

macro set(A::Symbol, B::Union{Number, Bool, String})
    clean_up(s) = String(filter(c->!occursin(c," ():"),s))
    key = clean_up(repr(A))
    quote
        Fermi.Options.set($key, $B)
    end
end

macro set(A::Symbol, B::Expr)
    clean_up(s) = String(filter(c->!occursin(c," ():"),s))
    key = clean_up(repr(A))

    if all(x->typeof(x)<:Number, B.args[2:end])
        quote
            Fermi.Options.set($key, $B)
        end
    else
        val = clean_up(repr(B))
        quote
            Fermi.Options.set($key, $val)
        end
    end
end

macro set(block::Expr)

    # Create an empty quote for the output
    out = quote end

    # Each argument in the expression block is a line
    # e.g.
    # @set {
    #    A B   -> First line (Expr object)
    #    C D   -> Second line (Expr object)
    #}
    lines = block.args

    for line in lines
        # args represent the arguments of each line. In the example above
        # args = {:A, :B} and {:C, :D}
        args = line.args

        # We must assert there are only two arguments
        length(args) == 2 || throw(FermiException("Too many arguments for @set: $line"))

        # A is the dictionary key, which is always taken as a String
        key = symbol_to_string(args[1])

        # B is the val
        val = args[2]
        valtype = typeof(val)

        # First we check the easy case: B is a Number or Bool
        if valtype <: Union{Number, Bool, String}
            push!(out.args, quote Fermi.Options.set($key, $val) end)

        # Now the case where B is a Symbol.
        elseif valtype === Symbol
            # Turn it into a String
            val = symbol_to_string(val)
            push!(out.args, quote Fermi.Options.set($key,$val) end)

        # Finally, if the value is an Expression
        elseif valtype === Expr

            # If the expression is purely numerical, we just used it as it is
            if all(x->typeof(x)<:Number, val.args[2:end])
                push!(out.args, quote Fermi.Options.set($key,$val) end)

            # If not, String it out
            else
                val = symbol_to_string(val)
                push!(out.args, quote Fermi.Options.set($key,$val) end)
            end
        else
            # Outside those cases the user did something sketchy
            throw(FermiException("Invalid input for $key: $val"))
        end
    end

    return out
end

"""
    Fermi.@reset(x="all", y...)

Restores any number of option to default setting. If one argument is "all", restores all options to their default values.

# Examples
```
julia> @set basis 6-31g
julia> @get basis
"6-31g"
julia> @reset basis
julia> @get basis
"sto-3g"

julia> @set {
            scf_guess core
            diis false
        }
julia> @reset
julia> @get scf_core
"gwh"
julia> @get diis
true
```
"""
macro reset(x="all", y...)
    Keys = [String(x)]
    for i in y
        push!(Keys, String(i))
    end

    quote
        for k in $Keys
            Fermi.Options.reset(k)
        end
    end
end

"""
    Fermi.@molecule

Set the molecule to be used in computations.

# Example
```
@molecule {
C                 -0.00000000     0.00000000    -0.00000000
O                  0.00000000     0.00000000    -2.19732722
O                 -0.00000000    -0.00000000     2.19732722
}
```
"""
macro molecule(block)
    mol = repr(block)
    mol = replace(mol, ";"=>"\n")
    mol = strip(filter(c->!occursin(c,"{}():"), mol))
    mol = String(mol)
    quote
        Fermi.Options.set("molstring", $mol)
    end
end

"""
    Fermi.@lookup

Look up valid keywords for Fermi.

# Example
```
julia> @lookup conv
┌────────────┬───────────────┐
│    Keyword │ Current Value │
├────────────┼───────────────┤
│  cc_e_conv │   1.00000e-10 │
│     d_conv │             8 │
│ scf_e_conv │   1.00000e-01 │
└────────────┴───────────────┘
```
"""
macro lookup(A::Symbol)
    A = strip(filter(c->!occursin(c,"{}():"), repr(A)))
    quote
        Keys = collect(keys(Fermi.Options.Default))
        filter!(s->occursin($A, s), Keys)
        l = length(Keys)
        if l == 0
            println("No keywords found containing: $($A)")
        else
            data = Array{Union{Nothing,String},2}(nothing,l,2)
            for i in eachindex(Keys)
                k = Keys[i]
                data[i,1] = k
                val = Fermi.Options.get(k)
                if val isa Bool || val isa Int
                    data[i,2] = "$val"
                elseif val isa String
                    # If multiple lines (e.g. molstring) display only the first line
                    if occursin("\n", val)
                        lines = split(val, "\n")
                        data[i,2] = lines[1] * " [+$(length(lines)-1) lines]"
                    else
                        data[i,2] = "$val"
                    end
                else
                    data[i,2] = format("{:5.5e}", val)
                end
            end
            pretty_table(data; header=["Keyword", "Current Value"])
        end
    end
end

macro freezecore()
    quote
        mol = Fermi.Molecule()
        atoms = mol.atoms
        # Dcore maps ranges of atomic numbers to number of core electrons
        Dcore = Fermi.PhysicalConstants.core_elec_perrow
        Ncore = 0
        for a in atoms
            for k in keys(Dcore)
                a.Z in k ? Ncore += Dcore[k] : nothing
            end
        end
        # Divide by two to get number of dropped occupied
        Ncore = Ncore >> 1
        Fermi.Options.set("drop_occ", Ncore)
    end
end


text:
using LinearAlgebra, SparseArrays, Laplacians, Statistics
include("src/cmg_precondition.jl")


text:
using ProgressMeter

# energy calculation using linked lists method
function utotal(box, DATA, first_atom, next_atom, nc)
  Ut  = 0.
@showprogress 1 "Computing potential energy"  for iat in 1:DATA.N
    icell = trunc(Int64,box[iat][1]/DATA.cutoff) + 1
    jcell = trunc(Int64,box[iat][2]/DATA.cutoff) + 1
    for i in icell-1:icell+1
      for j in jcell-1:jcell+1
        iw, jw = wrapcell(nc, i, j)
        jat = first_atom[iw,jw]
        while jat > 0
          if jat > iat
            rij = pbcseparation(box[jat],box[iat],DATA.side)
            if rij <= 2.
              Ut += upair(rij,DATA)
            end
          end
          jat = next_atom[jat]
        end
      end
    end
  end
  return Ut
end


text:
export ReorderedBasis, reorder, invorder

using LuxurySparse, LinearAlgebra

"""
    ReorderedBasis{N, T}

Lazy reorderd basis.
"""
struct ReorderedBasis{N, T <: Integer}
    orders::NTuple{N, T}
    takers::NTuple{N, T}
    differ::NTuple{N, T}
end

"""
    ReorderedBasis(orders::NTuple{N, <:Integer})

Returns a lazy set of reordered basis.
"""
ReorderedBasis(orders::NTuple{N, T}) where {N, T <: Integer} =
    ReorderedBasis{N, T}(orders, bmask.(orders), unsafe_sub(1:N, orders))

Base.length(::ReorderedBasis{N}) where N = 1 << N
Base.eltype(::ReorderedBasis{N, T}) where {N, T} = T
Base.getindex(it::ReorderedBasis, k::Int) = next_reordered_basis(k - 1, it.takers, it.differ)

function Base.iterate(it::ReorderedBasis{N}, state = 1) where N
    if state - 1 == 1 << N
        return nothing
    else
        return next_reordered_basis(state - 1, it.takers, it.differ), state + 1
    end
end

"""
    next_reordered_basis(basis, takers, differ)

Returns the next reordered basis accroding to current basis.
"""
function next_reordered_basis(basis::T, takers::NTuple{N, T}, differ::NTuple{N, T}) where {N, T}
    out = zero(T)
    # NOTE: do not add simd here, this will cause other loop become less efficient
    #       have some confidence with the compiler :-)
    for i in 1:N
        @inbounds out += (basis & takers[i]) << differ[i]
    end
    return out
end

"""
    unsafe_reorder(X::AbstractArray, orders)

Reorder `X` according to `orders`.

!!! warning

    `unsafe_reorder` won't check whether the length of `orders` and
    the size of first dimension of `X` match, use at your own risk.
"""
function unsafe_reorder end

function unsafe_reorder(v::AbstractVector, orders::NTuple{N, <:Integer}) where N
    nv = similar(v)
    for (k, b) in enumerate(ReorderedBasis(orders))
        @inbounds nv[k] = v[b + 1]
    end
    return nv
end

function unsafe_reorder(A::AbstractMatrix, orders::NTuple{N, <:Integer}) where N
    od = Vector{Int}(undef, 1 << length(orders))
    for (i, b) in enumerate(ReorderedBasis(orders))
        @inbounds od[i] = 1 + b
    end

    od = invperm(od)
    return A[od, od]
end

unsafe_reorder(A::IMatrix, orders::NTuple{N, <:Integer}) where N = A

function unsafe_reorder(A::PermMatrix, orders::NTuple{N, <:Integer}) where N
    od = Vector{Int}(undef, 1 << length(orders))
    for (i, b) in enumerate(ReorderedBasis(orders))
        @inbounds od[i] = 1 + b
    end

    perm = similar(A.perm)
    vals = similar(A.vals)

    @simd for i = 1:length(perm)
        @inbounds perm[od[i]] = od[A.perm[i]]
        @inbounds vals[od[i]] = A.vals[i]
    end

    return PermMatrix(perm, vals)
end

function unsafe_reorder(A::Diagonal, orders::NTuple{N, <:Integer}) where N
    diag = similar(A.diag)
    for (i, b) in enumerate(ReorderedBasis(orders))
        @inbounds diag[b+1] = A.diag[i]
    end
    return Diagonal(diag)
end

function _reorder(v::AbstractArray, orders::NTuple{N, <:Integer}) where N
    length(orders) == log2dim1(v) || throw(DimensionMismatch("size of array not match length of order"))
    return unsafe_reorder(v, orders)
end

"""
    reorder(X::AbstractArray, orders)

Reorder `X` according to `orders`.

!!! tip
    Although `orders` can be any iterable, `Tuple` is preferred inorder to
    gain as much performance as possible. But the conversion won't take much
    anyway.
"""
function reorder end

reorder(v::AbstractArray, itr) = _reorder(v, Tuple(itr))
reorder(v::AbstractArray, orders::Tuple) = _reorder(v, promote(orders...))
reorder(v::AbstractArray, orders::NTuple{N, <:Integer}) where N = _reorder(v, orders)

"""
    invorder(X::AbstractVecOrMat)

Inverse the order of given vector/matrix `X`.
"""
invorder(X::AbstractVecOrMat) = reorder(X, Tuple(log2dim1(X):-1:1))


text:
function _writevalue(io::IO, value::AbstractString, delim, quotechar::Nothing, escapechar, nastring)
    print(io, value)
end

function _writevalue(io::IO, value::AbstractString, delim, quotechar, escapechar, nastring)
    print(io, quotechar)
    for c in value
        if c==quotechar || c==escapechar
            print(io, escapechar)
        end
        print(io, c)
    end
    print(io, quotechar)
end

function _writevalue(io::IO, value, delim, quotechar, escapechar, nastring)
    print(io, value)
end

function _writevalue(io::IO, value::DataValue{T}, delim, quotechar, escapechar, nastring) where {T}
    if isna(value)
        print(io, nastring)
    else
        _writevalue(io, get(value), delim, quotechar, escapechar, nastring)
    end
end

@generated function _writecsv(io::IO, it, ::Type{T}, delim, quotechar, escapechar, nastring) where {T}
    col_names = fieldnames(T)
    n = length(col_names)
    push_exprs = Expr(:block)
    for i in 1:n
        push!(push_exprs.args, :( _writevalue(io, i.$(col_names[i]), delim, quotechar, escapechar, nastring) ))
        if i<n
            push!(push_exprs.args, :( print(io, delim ) ))
        end
    end
    push!(push_exprs.args, :( println(io) ))

    quote
        for i in it
            $push_exprs
        end
    end
end

function _save(io, data; delim=',', quotechar='"', escapechar='"', nastring="NA", header=true)
    isiterabletable(data) || error("Can't write this data to a CSV file.")

    it = getiterator(data)
    colnames = collect(eltype(it).parameters[1])

    if header
        if quotechar===nothing
            join(io,[string(colname) for colname in colnames],delim)
        else
            join(io,["$(quotechar)" * replace(string(colname), quotechar => "$(escapechar)$(quotechar)") * "$(quotechar)" for colname in colnames],delim)
        end
        println(io)
    end
    _writecsv(io, it, eltype(it), delim, quotechar, escapechar, nastring)
end

function _save(filename::AbstractString, data; delim=',', quotechar='"', escapechar='"', nastring="NA", header=true)
    isiterabletable(data) || error("Can't write this data to a CSV file.")

    ext = last(split(filename, '.'))

    if ext == "gz" # Gzipped
        open(GzipCompressorStream, filename, "w") do io
            _save(io, data, delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring,  header=header)
        end
    else
        open(filename, "w") do io
            _save(io, data, delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring,  header=header)
        end
    end
end

function fileio_save(f::FileIO.File{FileIO.format"CSV"}, data; delim=',', quotechar='"', escapechar='"', nastring="NA", header=true)
    return _save(f.filename, data, delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring, header=header)
end

function fileio_save(f::FileIO.File{FileIO.format"TSV"}, data; delim='\t', quotechar='"', escapechar='"', nastring="NA", header=true)
    return _save(f.filename, data, delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring, header=header)
end

function fileio_save(s::FileIO.Stream{FileIO.format"CSV"}, data; delim=',', quotechar='"', escapechar='"', nastring="NA", header=true)
    return _save(s.io, data, delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring, header=header)
end

function fileio_save(s::FileIO.Stream{FileIO.format"TSV"}, data; delim='\t', quotechar='"', escapechar='"', nastring="NA", header=true)
    return _save(s.io, data, delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring, header=header)
end

#
# Streaming version writes header (if any) on first call, then appends on subsequent calls.
#
const CSV_or_TSV = Union{FileIO.format"CSV", FileIO.format"TSV"}

_delim(T) = T <: FileIO.format"CSV" ? ',' : '\t'

mutable struct CSVFileSaveStream{T}
    io::T
    first_data_written::Bool
    delim::Char
    quotechar::Char
    escapechar::Char
    nastring::AbstractString
    header::Bool
end

function fileio_savestreaming(f::FileIO.File{T}, data=nothing; delim=_delim(T), quotechar='"', escapechar='"', nastring="NA",
                              header=true) where T <: CSV_or_TSV
    io = open(f.filename, "w")

    if data!==nothing
        _save(io, data; delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring, header=header)
    end

    return CSVFileSaveStream(io, data!==nothing, delim, quotechar, escapechar, nastring, header)
end

function fileio_savestreaming(s::FileIO.Stream{T}, data=nothing; delim=_delim(T), quotechar='"', escapechar='"', nastring="NA",
                              header=false) where T <: CSV_or_TSV

    if data!==nothing
        _save(s.io, data; delim=delim, quotechar=quotechar, escapechar=escapechar, nastring=nastring, header=header)
    end

    return CSVFileSaveStream(s.io, data!==nothing, delim, quotechar, escapechar, nastring, header)
end

function Base.write(s::CSVFileSaveStream, data)
    _save(s.io, data; delim=s.delim, quotechar=s.quotechar, escapechar=s.escapechar, nastring=s.nastring, header=s.first_data_written ? false : header)
end

function Base.close(s::CSVFileSaveStream)
    close(s.io)
end