text:
using ChaosTools
using DynamicalSystemsBase
using Test
using OrdinaryDiffEq

@testset "Uncertainty exponent / fractal boundaries" begin

@testset "Test uncertainty orginal paper" begin
    ds = Systems.grebogi_map(rand(2))
    θg = range(0, 2π, length = 251)
    xg = range(-0.5, 0.5, length = 251)
    bsn, att = basins_of_attraction((θg, xg), ds; show_progress = false)
    e, f, α = uncertainty_exponent(bsn; range_ε = 3:15)
    # In the paper the value is roughly 0.2
    @test (0.2 ≤ α ≤ 0.3)
end

@testset "Test uncertainty Newton map" begin

function newton_map(dz,z, p, n)
    f(x) = x^p[1]-1
    df(x)= p[1]*x^(p[1]-1)
    z1 = z[1] + im*z[2]
    dz1 = f(z1)/df(z1)
    z1 = z1 - dz1
    dz[1]=real(z1)
    dz[2]=imag(z1)
    return
end

# dummy function to keep the initializator happy
function newton_map_J(J,z0, p, n)
   return
end

ds = DiscreteDynamicalSystem(newton_map,[0.1, 0.2], [3] , newton_map_J)
xg = yg = range(-1.,1.,length=300)
bsn,att = basins_of_attraction((xg, yg), ds; show_progress = false)
e,f,α = uncertainty_exponent(bsn; range_ε = 5:30)

# Value (published) from the box-counting dimension is 1.42. α ≃ 0.6
@test (0.55 ≤ α ≤ 0.65)

end

@testset "Basin entropy and Fractal test" begin
    ds = Systems.grebogi_map()
    θg=range(0,2π,length = 300)
    xg=range(-0.5,0.5,length = 300)
    basin, attractors = basins_of_attraction((θg,xg), ds; show_progress = false)
    Sb, Sbb = basin_entropy(basin, 6)
    @test 0.4 ≤ Sb ≤ 0.42
    @test 0.6 ≤ Sbb ≤ 0.61

    test_res, Sbb = basins_fractal_test(basin; ε = 5)
    @test test_res == :fractal

    ds = Systems.henon(zeros(2); a = 1.4, b = 0.3)
    xg = yg = range(-2.,2.,length = 300)
    basin, attractors = basins_of_attraction((xg,yg), ds; show_progress = false)
    test_res, Sbb = basins_fractal_test(basin; ε = 5)
    @test test_res == :smooth
end

end


text:
using DataFrames

#cd("/Users/VGupta/Documents/Research/Julia Stuff/Demand")
include("BertrandNash2.jl");  # does the generation of the prices
include("fitNash2.jl");  # does the fitting

thetas1 = [-1.2, .5, 1, -10]
thetas2 = [.3, -1, 1, -10]
pmax = .45; N = 200;

function genErrors( ; mu1 = 5, mu2 = 5, mu3 = 5, sigma1 = 1.5, sigma2 = 1.5, sigma3 = 1.5)
    eps1 = rand(Normal(mu1, sigma1)); eps2 = rand(Normal(mu2, sigma2)); eps3 = rand(Normal(mu3, sigma3))
    eps3 = (eps1 + eps2 + eps3)/3
	return eps1, eps2, eps3
end

srand(8675309)
inPrices = Array(Float64, N, 2)
inErrors = Array(Float64, N, 3)
for iRun = 1:N
    inErrors[iRun, :] = [genErrors()...]
    pstar = solveNashPrices(thetas1, thetas2, inErrors[iRun, :], pmax, prices=[.3, .3])
	inPrices[iRun, :] = pstar'
end
indx = find( abs(inErrors[:, 3] - median(inErrors[:, 3]) ) .<= 1e-2 )[1];

function fit(inPrices, inErrors, c, lambda_reg)
    ### Thin out the grids
    # p1_grid = sort(inPrices[:, 1]); p2_grid = sort(inPrices[:, 2])
    # p1_grid = thin(p1_grid, .45/50, pmax)
    # p2_grid = thin(p2_grid, .45/50, pmax)
    p1_grid = p2_grid = linspace(.1, pmax, 50)

    #Augment the data as necessary
    N1 = length(p1_grid); N2 = length(p2_grid); totN = N + N1 + N2
    data = [inPrices inErrors[:, 3]]
    data = [data; p1_grid repmat([inPrices[indx, 2]], N1, 1) repmat([inErrors[indx, 3]], N1, 1)]
    data = [data; repmat([inPrices[indx, 1]], N2, 1) p2_grid repmat([inErrors[indx, 3]], N2, 1)]
    data = [data; p1_grid[1] inPrices[indx, 2] inErrors[indx, 3];
                  inPrices[indx, 1] p2_grid[1] inErrors[indx, 3]]

    #build a kernel, the graham matrix, and it's cholesky
    #only use significant eigenvalues
    k(x,y) = exp( -1/c * norm(x-y)^2 )
    K = [k(data[i, :], data[j, :]) for i =1:N+N1+N2+2, j=1:totN+2]
    D, V = eig(K)
    tailN = find( 1 - (cumsum(D) / sum(D)) .<= 1-1e-6)[1]
    Dbar = copy(D)
    Dbar[1:tailN-1] = 0;
    Vbar = V * diagm(sqrt(Dbar));

    m, vZ1, vZ2, vf1, vf2, vIntercept = setUpFitting(totN+2, Vbar[:, tailN:totN+2])

    normalize(m, vf1[totN + 1], margRev(thetas1, data[totN + 1, 1:2], data[totN + 1, 3], :1))
    normalize(m, vf2[totN + 2], margRev(thetas2, data[totN + 2, 1:2], data[totN + 2, 3], :2))
    addDecreasing(m, vf1, (N+1):N+N1, TOL=1e-6)
    addDecreasing(m, vf2, (N+N1+1):totN, TOL=1e-6)

    resids = Variable[]
    @defVar(m, maxVal >=0)
    for iRun = 1:N
        resid = addResids(m, vf1, vf2, iRun, data, pmax)
        push!(resids, resid)
        @addConstraint(m, maxVal >= resid)
    end

    @defVar(m, reg_term[1:2])
    HNorm1 = QuadExpr(vZ1[:], vZ1[:], ones(length(vZ1)), AffExpr())
    HNorm2 = QuadExpr(vZ2[:], vZ2[:], ones(length(vZ2)), AffExpr())
    addConstraint(m, HNorm1 <= reg_term[1])
    addConstraint(m, HNorm2 <= reg_term[2])

    @setObjective(m, Min, sum{resids[i], i=1:N}/N + lambda_reg * (reg_term[1] + reg_term[2]))
    println(solve(m))

    #extract solution
    fZ1 = [getValue(vZ1[i])::Float64 for i =1:length(vZ1)]
    fZ2 = [getValue(vZ2[i])::Float64 for i = 1:length(vZ2)]
    fIntercept = [getValue(vIntercept[i])::Float64 for i=1:2]
    fMR1_fun = createfun(fZ1, fIntercept[1], Vbar[:, tailN:totN+2], data, k)
    fMR2_fun = createfun(fZ2, fIntercept[2], Vbar[:, tailN:totN+2], data, k)

    return fMR1_fun, fMR2_fun
end

srand(516746226)
function testOut(fMR1_fun, fMR2_fun, Nout, METHOD)
    if METHOD == :pdiffs
        fPrices = Array(Float64, Nout, 2)
        oPrices = Array(Float64, Nout, 2)
        oErrors = Array(Float64, Nout, 3)
        for iRun = 1:Nout
            oErrors[iRun, :] = [genErrors()...]
            pstar = solveNashPrices(thetas1, thetas2, oErrors[iRun, :], pmax)
            oPrices[iRun, :] = pstar'

           println(iRun)

            fPrices[iRun, :] = solveNashPrices(fMR1_fun, fMR2_fun, oErrors[iRun, 3], pmax,
                                                   prices=copy(pstar), trace=false)
        end
        price_diffs = oPrices - fPrices
        return norm(price_diffs) / sqrt(Nout)
    else
        eps = 0.
        oErrors = Array(Float64, Nout, 3)
        for iRun = 1:Nout
            oErrors[iRun, :] = [genErrors()...]
            pstar = solveNashPrices(thetas1, thetas2, oErrors[iRun, :], pmax)
            MR1 = fMR1_fun(pstar..., oErrors[iRun, 3])
            MR2 = fMR2_fun(pstar..., oErrors[iRun, 3])

            y1 = max(0, MR1); y2 = max(0, MR2)
            eps += max(0, pmax * (y1 + y2) - pstar[1] * MR1 -pstar[2] * MR2 )
        end
        return eps/Nout
    end
end

lamb_grid = 10. .^(-6:.5:1)
c_grid = 10. .^(-3:.25:1)

#1 REdo this with tighter lambda, larger c
# lamb_grid = 10. .^(-6:.5:-2)
# c_grid = 10. .^(1:.5:3)
fold_dict = Dict{(Float64, Float64), Array{Float64, 1}}()

NSplits = 10
for (ix, l) in enumerate(lamb_grid)
    for (jx, c) in enumerate(c_grid)
        if haskey(fold_dict, (l,c))
            continue
        end

        fMR1_fun, fMR2_fun = fit(inPrices, inErrors, c, l)
        errors = zeros(NSplits)
        good_indices = Int[]
        for iSplit = 1:NSplits
            try
                errors[iSplit] = testOut(fMR1_fun, fMR2_fun, 200, :pdiffs)
                push!(good_indices, iSplit)
            catch y
                println("Error for $l $c \t", y)
            end
        end
        println("Good_indices: $(length(good_indices))")
        println(good_indices)
        println(errors)
        println(errors[good_indices])
        if length(good_indices) > 0
            s = length(good_indices) > 1 ? std(errors[good_indices]) : -1.
            fold_dict[(l, c)] = [mean(errors[good_indices]), s]
        end
    end
end

#Create an array of the relevant values
fold_vals = Array(Float64, 0, 4)
for l = lamb_grid
    for c = c_grid
        if haskey(fold_dict, (l, c))
            fold_vals = [fold_vals ; [l c fold_dict[(l,c)][1] fold_dict[(l,c)][2] ] ]
            println("$l \t $c \t $(fold_dict[(l,c)][1]) \t $(fold_dict[(l,c)][2])")
        end
    end
end

df = DataFrame(fold_vals)
colnames!(df, ["Lambda", "C", "Mean", "SD"])
writetable("crossVal_pdiffs.csv", df)

ext: jl

max_stars_repo_name: nlperic/InverseVIs

max_issues_repo_name: nlperic/InverseVIs


text:
include("polyagamma/polyagamma_mcmc.jl")
include("mmle/mmle.jl")

function assess_examinee!(
    examinee::AbstractExaminee,
    items::Vector{<:AbstractItem},
    responses::Vector{Response};
    method = "mmle",
    kwargs...
    )
    if method == "pg"
        return assess_examinee_pg!(item, examinees, responses; kwargs...)
    elseif method =="mmle"
        return assess_examinee_mmle!(item, examinees, responses; kwargs...)
    end
end


text:
__precompile__(true)

# main module for `LazySets.jl`
module LazySets

using Requires, SparseArrays, LinearAlgebra
using LinearAlgebra: checksquare
import LinearAlgebra: norm, ×
import Random
using Random: AbstractRNG, GLOBAL_RNG, SamplerType, shuffle
import InteractiveUtils: subtypes

export Arrays, Approximations
export ×

# =======================
# Arrays auxiliary module
# =======================
include("Arrays/Arrays.jl")
using .Arrays

# ===================
# Auxiliary functions
# ===================
include("helper_functions.jl")
include("comparisons.jl")
include("macros.jl")
include("samples.jl")

# ==================
# Abstract set types
# ==================
include("LazySet.jl")
include("AbstractPolyhedron.jl")
include("HalfSpace.jl") # must be here to make LinearConstraint available
include("AbstractPolyhedron_functions.jl")
include("AbstractPolytope.jl")
include("AbstractCentrallySymmetric.jl")
include("AbstractCentrallySymmetricPolytope.jl")
include("AbstractZonotope.jl")
include("AbstractHyperrectangle.jl")
include("AbstractPolygon.jl")
include("AbstractSingleton.jl")
include("AbstractHPolygon.jl")

# =============================
# Types representing basic sets
# =============================
include("Ball1.jl")
include("Ball2.jl")
include("BallInf.jl")
include("Ballp.jl")
include("Ellipsoid.jl")
include("EmptySet.jl")
include("HPolygon.jl")
include("HPolygonOpt.jl")
include("HPolytope.jl")
include("HPolyhedron.jl")
include("Hyperplane.jl")
include("Hyperrectangle.jl")
include("Interval.jl")
include("Line.jl")
include("LineSegment.jl")
include("Singleton.jl")
include("Universe.jl")
include("VPolygon.jl")
include("VPolytope.jl")
include("ZeroSet.jl")
include("Zonotope.jl")
include("Complement.jl")

# ==================================
# Types representing non-convex sets
# ==================================
include("PolynomialZonotope.jl")

# ===================================================
# Algorithms to compute the convex hull of polygons
# ===================================================
include("concrete_convex_hull.jl")

# =================================
# Types representing set operations
# =================================
include("CartesianProduct.jl")
include("ConvexHull.jl")
include("ExponentialMap.jl")
include("Intersection.jl")
include("LinearMap.jl")
include("MinkowskiSum.jl")
include("ResetMap.jl")
include("SymmetricIntervalHull.jl")
include("Translation.jl")
include("UnionSet.jl")
include("Rectification.jl")

# =============================
# Conversions between set types
# =============================
include("intersection_helper.jl")
include("convert.jl")

# =====================
# Approximations module
# =====================
include("Approximations/Approximations.jl")

# ===========================
# Concrete operations on sets
# ===========================
include("concrete_intersection.jl")
include("is_intersection_empty.jl")
include("is_subset.jl")
include("difference.jl")

# =======
# Aliases
# =======
include("aliases.jl")

# ==========================
# Parallel algorithms module
# ==========================
include("Parallel/Parallel.jl")

# ============
# Plot recipes
# ============
include("plot_recipes.jl")

# ===================================================
# Load external packages on-demand (using 'Requires')
# ===================================================
include("init.jl")

end # module


text:
#= Node and edge types

* Each graph consists of `Node`s that store information about connections with
each other using pointers of type `NodePtr`.

* Based on those pointers and their corresponding node indices, `Edge` instances
are generated dynamically upon user's request.

* If some metadata needs to be stored, `AbstractAttribute` and its children are used.
They map node or edge indices to their corresponding values.
=#

"""
    NodePtr{N,E}

Internal type for storing information about connections in the graph:

* `node` — id of target node,
* `id`   — id of the corresponding edge.

`NodePtr` instances are compared by IDs of nodes they point to.
"""
struct NodePtr{N,E}
    node :: N
    id   :: E
end

isless(e::NodePtr{N,E}, f::NodePtr{N,E}) where {N,E} = <(e.node, f.node)
isless(e::NodePtr{N,E}, f::N) where {N,E} = <(e.node, f)
isless(e::N, f::NodePtr{N,E}) where {N,E} = <(e, f.node)

==(e::N, f::NodePtr{N,E}) where {N,E} = ==(e, f.node)
==(f::NodePtr{N,E}, e::N) where {N,E} = ==(e, f.node)

"""
    Node{N,E}

A single node in the graph.
Its `forward` and `reverse` fields hold vectors of pointers to other nodes.
"""
struct Node{N,E}
    forward :: Vector{NodePtr{N,E}}
    reverse :: Vector{NodePtr{N,E}}
end

Node{N,E}() where {N,E} = Node{N,E}(NodePtr{N,E}[], NodePtr{N,E}[])

fwd_ptrs(n::Node) = n.forward
fwd_ptrs(n::Node, ::Type{Forward}) = n.forward
fwd_ptrs(n::Node, ::Type{Reverse}) = n.reverse
fwd_ptrs(n::Node, ::Type{Both}) = n.forward

rev_ptrs(n::Node) = n.reverse
rev_ptrs(n::Node, ::Type{Forward}) = n.reverse
rev_ptrs(n::Node, ::Type{Reverse}) = n.forward
rev_ptrs(n::Node, ::Type{Both}) = n.reverse

ptr_length(n::Node, ::Type{Forward}) = length(n.forward)
ptr_length(n::Node, ::Type{Reverse}) = length(n.reverse)
ptr_length(n::Node, ::Type{Both}) = length(n.forward) + length(n.reverse)

findfirst_ptr(x::Vector{NodePtr{N,E}}, i::N) where {N,E} = findfirst(x, N(i))

has_ptr(x::Vector{NodePtr{N,E}}, n::N) where {N,E} =
    let i = searchsortedlast(x, n)
        i != 0 && x[i] == n
    end

get_ptr(n::Node, i::Integer, ::Type{Forward}) = n.forward[i]
get_ptr(n::Node, i::Integer, ::Type{Reverse}) = n.reverse[i]
get_ptr(n::Node, i::Integer, ::Type{Both}) =
    let l = length(n.reverse)
        if i > l
            return n.forward[i - l]
        else
            return n.reverse[i]
        end
    end

function add_ptr!(vec::Vector{T}, val::T) where T
    i = searchsortedlast(vec, val) + 1
    insert!(vec, i, val)
end

"""
    add_one_ptr!(vec, val)

Insert `val` into a sorted vector `vec` if it is not there already.
Returns `true` or `false` depending on whether it was inserted.
"""
function add_one_ptr!(vec::Vector{T}, val::T) where T
    i = searchsortedlast(vec, val)
    @inbounds if i == 0 || vec[i] != val
        insert!(vec, i + 1, val)
        return true
    else
        return false
    end
end

"""
    delete_ptr!()

Internal function to remove a matching NodePtr from a Vector.
"""
function delete_ptr!(ptrs::Vector{NodePtr{N,E}}, nid::N, eid::E) where {N,E}
    @inbounds for i = searchsorted(ptrs, nid)
        if ptrs[i].id == eid
            deleteat!(ptrs, i)
            return true
        end
    end
    return false
end

function delete_ptr!(ptrs::Vector{NodePtr{N,E}}, nid::N) where {N,E}
    i = searchsortedlast(ptrs, nid)
    if i > 0 && ptrs[i] == nid
        deleteat!(ptrs, i)
        return true
    end
    return false
end

function swap_ptr!(ptrs::Vector{NodePtr{N,E}}, nid::N, eid::E, new_nid::N) where {N,E}
    @_inline_meta
    delete_ptr!(ptrs, nid, eid)
    add_ptr!(ptrs, NodePtr{N,E}(new_nid, eid))
end

"""
    Edge{N,E}

A single edge in the graph. It has the following fields:

* `source` — id of the source node,
* `target` — id of the target node,
* `id`     — its own unique id,
* `isdir`  — whether it is directed.
"""
struct Edge{N,E}
    source :: N
    target :: N
    id     :: E
    isdir  :: Bool
end

"""
    AbstractAttribute{T}

Abstract parent for graph attribute types subtyping `AbstractVector`.
To achieve maximum performance, attributes check bounds only on element assignment.
"""
abstract type AbstractAttribute{T,N,F<:Function} <: AbstractArray{T,N} end

const AttributeDict = Dict{Symbol,AbstractAttribute}

"""
    defval(::Type{T})

Get a reasonable default value for type `T`.
Works like `Base.zero`, but has implementations for more types.
"""
defval(::Type{T}) where {T<:Number} = zero(T)
defval(::Type{String}) = ""
defval(::Type{Symbol}) = :|
defval(::Type{Char}) = '\0'
defval(::Type{Nothing}) = nothing
defval(::Type{RGB}) = RGB(0, 0, 0)
defval(::Type{Tuple{Vararg{T,N}}}) where {T,N} = tuple([defval(T) for _ = 1:N]...)
defval(::Type{Any}) = nothing
defval(::Type{T}) where {T} = error("no default value is set for type $T")


text:
module Grids

export
    AbstractTopology, Periodic, Bounded, Flat, topology,
    AbstractGrid, RegularCartesianGrid, VerticallyStretchedCartesianGrid

import Base: size, length, eltype, show

using Oceananigans

"""
    AbstractTopology

Abstract supertype for grid topologies.
"""
abstract type AbstractTopology end

"""
    Periodic

Grid topology for periodic dimensions.
"""
struct Periodic <: AbstractTopology end

"""
    Bounded

Grid topology for bounded dimensions. These could be wall-bounded dimensions
or dimensions
"""
struct Bounded <: AbstractTopology end

"""
    Flat

Grid topology for flat dimensions, generally with one grid point, along which the solution
is uniform and does not vary.
"""
struct Flat <: AbstractTopology end

"""
    AbstractGrid{FT, TX, TY, TZ}

Abstract supertype for grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractGrid{FT, TX, TY, TZ} end

eltype(::AbstractGrid{FT}) where FT = FT
topology(::AbstractGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = (TX(), TY(), TZ())
topology(grid, dim) = topology(grid)[dim]

size(grid::AbstractGrid) = (grid.Nx, grid.Ny, grid.Nz)
length(grid::AbstractGrid) = (grid.Lx, grid.Ly, grid.Lz)
halo_size(grid) = (grid.Hx, grid.Hy, grid.Hz)

total_size(a) = size(a) # fallback

total_size(grid::AbstractGrid) = (grid.Nx + grid.Hx,
                                  grid.Ny + grid.Hy,
                                  grid.Nz + grid.Hz)

include("grid_utils.jl")
include("regular_cartesian_grid.jl")
include("vertically_stretched_cartesian_grid.jl")

end


text:
mutable struct JaxArray{T,N} <: AbstractJaxArray{T,N}
  o::PyObject
  dims::NTuple{N,Int}
end

const JaxVector{T} = JaxArray{T,1}
const JaxMatrix{T} = JaxArray{T,2}

JaxArray(arr::AbstractArray) = JaxArray(collect(arr))
function JaxArray(arr::Array{T,N}) where {T,N}
  # When converting an array to Jax Array we transpose it.
  arr = PyReverseDims(arr)
  return numpy.array(arr)
end

JaxNumber(a::Number) = numpy.array(a)

function JaxArray(parr::PyObject)
  if !pyisinstance(parr, jax.interpreters.xla.DeviceArray)
    error("It is not jax.interpreters.xla.DeviceArray but $parr ")
  end
  T = np_to_jl_type(parr.dtype)
  N = length(parr.shape)
  return JaxArray{T,N}(parr, reverse(parr.shape))
end

JaxArray(d::Union{Integer,Dims}...) = JaxArray(Float32, d...)
JaxArray(T::Type{<:Number}, d...) = JaxArray(T, undef, d...)
JaxArray(T::Type{<:Number}, u::UndefInitializer, d::Integer...) =
  JaxArray(T, u, Dims(d))
function JaxArray(T::Type{<:Number}, u::UndefInitializer, d::Dims)
  return numpy.empty(reverse(d), dtype = jl_to_np_type(T))#JaxArray()
end

Base.size(a::JaxArray) = a.dims
Base.length(a::JaxArray) = prod(a.dims)
Base.eltype(a::JaxArray{T}) where {T} = T
Base.axes(a::JaxArray) = map(Base.OneTo, size(a))

PyCall.PyObject(a::JaxArray) = a.o
Base.convert(::Type{<:JaxArray}, o::PyObject) = JaxArray(o)

_transpose(x::Union{AbstractVector,AbstractMatrix}) = transpose(x)
_transpose(x::AbstractArray{T,N}) where {T,N} = permutedims(x, reverse(1:N))

_values(a::JaxArray) = PyArray_ReadOnly(a.o._value) |>_transpose
_values(a::JaxArray{T,0}) where T = PyArray_ReadOnly(a.o._value)
_values(a::JaxArray{T,1}) where T = PyArray_ReadOnly(a.o._value)
Base.collect(a::JaxArray) = copy(_values(a))
Base.convert(::Type{<:Array}, a::JaxArray) = copy(_values(a))

# show
_szstr(d::Tuple{}) = "$(0)-dimensional"
_szstr(d::Tuple{Int}) = "$(first(d))-element"
_szstr(dims::Dims) = prod(["$i×" for i in dims])[1:end-1]
function Base.show(io::IO, m::MIME"text/plain", a::JaxArray{T,N}) where {T,N}
  print(io, "$(_szstr(a.dims)) JaxArray{$T,$N}")
  isempty(a) && return
  print(io, ":")

  if !haskey(io, :compact) && length(axes(a, 2)) > 1
      io = IOContext(io, :compact => true)
  end
  if get(io, :limit, false) && eltype(a) === Method
      # override usual show method for Vector{Method}: don't abbreviate long lists
      io = IOContext(io, :limit => false)
  end

  if get(io, :limit, false) && displaysize(io)[1]-4 <= 0
      return print(io, " …")
  else
      println(io)
  end

  io = IOContext(io, :typeinfo => eltype(a))

  Base.print_array(io, _values(a))
end

function Base.show(io::IO, a::JaxArray{T,N}) where {T,N}
  print(io, "$(_szstr(a.dims)) JaxArray{$T,$N}:")
  vals = _values(a)
  if N == 0
    print(io, "[$(first(vals))]\n")
  else
    print(io, "[$(first(vals)) … $(last(vals))]\n")
  end
end

function Base.show(io::IO, m::MIME"text/plain", a::JaxArray{Bool,0})
  print(io, "JaxScalar{Bool}[$((first(_values(a)))==true)]")
end

function Base.show(io::IO, m::MIME"text/plain", a::JaxArray{T,0}) where T
  print(io, "JaxScalar{$T}[$((first(_values(a))))]")
end

# very slow
# Base.getindex(a::JaxArray, args...) = (getindex(_values(a), args...))

# iterators
Base.eachrow(a::JaxArray) = (get(a.o, i-1) for i in axes(a, 1))

block_until_ready(a::JaxArray) = a.o.block_until_ready()

##
#=
abstract type AbstractJaxNumber end

mutable struct JaxNumber{T} <: AbstractJaxNumber
  o::PyObject
end

function JaxNumber(n::Array{T,0}) where {T}
  # When converting an array to Jax Array we transpose it.
  return JaxNumber(numpy.array(n))
end
=#


text:
import Sockets: send

export getSerialDevices, resolvedSymlink

function resolvedSymlink(port::String)
  if islink(port)
    resolvedPort = joinpath("/dev", readlink(port))
  else
    resolvedPort = port
  end
  return resolvedPort
end

mutable struct SerialDevice
	sp::SerialPort
	pause_ms::Int
	timeout_ms::Int
	delim_read::String
	delim_write::String
end

"""
Set time to wait in ms for command to be send to serial device.
"""
function set_pause_ms(sd::SerialDevice,pause_ms::Int)
	sd.pause_ms = pause_ms
	return nothing
end

"""
Set maximal time to wait for querry answer in ms.
"""
function set_timeout_ms(sd::SerialDevice,timeout_ms::Int)
	sd.timeout_ms = timeout_ms
	return nothing
end

"""
Set character which terminates querry.
"""
function set_delim_write(sd::SerialDevice,delim::String)
	sd.delim_write = delim
	return nothing
end

"""
Set character which terminates querry answer.
"""
function set_delim_read(sd::SerialDevice,delim::String)
	sd.delim_read = delim
	return nothing
end

"""
Send command string to serial device.
"""
function send(sd::SerialDevice,cmd::String)
	@debug cmd
	write(sd.sp,string(cmd,sd.delim_write))
	sleep(sd.pause_ms/1000)
	return nothing
end

"""
Read out current content of the output buffer of the serial devive. Returns a String.
"""
function receive(sd::SerialDevice)
	return readstring(sd.sp)
end

"""
Send querry to serial device and receive device answer. Returns a String
"""
function query(sd::SerialDevice,cmd::String)
	flush(sd.sp)
	send(sd,string(cmd,sd.delim_write))
	out = readuntil(sd.sp, Vector{Char}(sd.delim_read), sd.timeout_ms)
	flush(sd.sp)
	return rstrip(out,Vector{Char}(sd.delim_read))
end

"""
Close the serial port of the serial device `sd`. The optional `delete` keyword
argument triggers a call to `sp_free_port` in the C library if set to `true`.
"""
function Base.close(sd::SerialDevice; delete::Bool=false)
	sd.sp = close(sd.sp,delete=delete)
	return sd
end

"""
Read out current Serial Ports, returns `Array{String,1}`
"""
function getSerialDevices()
    nports_guess=64
    ports = sp_list_ports()
    devices=Array{String,1}()
    portArray=unsafe_wrap(Array, ports, nports_guess, false)
    for k=1:length(portArray)
        if portArray[k] == C_NULL
          return devices
        else
          try
            device=sp_get_port_name(portArray[k])
            push!(devices,device)
          catch e
            @warn "Exception" e
          end
        end
    end
    return devices
end


text:
using LinearAlgebra
using ColorTypes, ColorVectorSpace

# Fenchel-Moreau conjugate function
fmconj(f) = (x_, s_) -> (x_ - f(x_,s_))

# prox of l2-norm
proxL2(x::AbstractArray, λ::Real) = max(1.0 - λ/norm(x), 0) .* x

# prox. of mixed l1- and l2- norm
groupshrink(x, λ::Real) = proxL2.(x, λ)

# prox. of l1-norm
softshrink(x::AbstractArray, λ::Real) = softshrink.(x, λ)
softshrink(x, λ::Real) = max(1.0 - λ / abs(x), 0) * x

# prox. of nuclear norm.
function shrinkSingularValue(A::AbstractMatrix, λ::Real)
    X = svd(A, full=false)
    X.U * Diagonal(softshrink(X.S, λ)) * X.Vt
end

function boxProj(x::AbstractArray, lowerBound::Real, upperBound::Real)
    max.(min.(x,upperBound),lowerBound)
end

function l2ballProj(x::AbstractArray, radius::Real, c::AbstractArray)
    dist = norm(x - centerVec)

    ifelse(dist <= radius, x, (λ / dist) * (x - c) + c)
end

function hardshrink(x::AbstractArray, k::Integer; lt::Function=(lhs,rhs)->isless(abs2(lhs), abs2(rhs)))
    nzids = sortperm(vec(x); lt=lt, rev=true)[1:k]
    foldl(nzids; init=zero(x)) do mtx, idx
        setindex!(mtx, x[idx], idx)
    end
end

function hardshrink(x::AbstractArray, ks::AbstractArray)
    hardshrink.(x, ks)
end

function hardshrink(x::AbstractArray, ks::NTuple)
    hardshrink.(x, ks)
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
module BenchRandomIncrements

using BenchmarkTools
using TSXPlayground:
    parallel_random_increments!,
    parallel_random_increments_adaptive!,
    random_increments_problem_clusters,
    random_increments_problem_uniform,
    serial_random_increments!

function default_ntasks_list()
    ntasks_list = [Threads.nthreads()]
    if Threads.nthreads() > 2
        pushfirst!(ntasks_list, 2)
    end
    return ntasks_list
end

function default_problems(;
    ncounters_list = [2^30],
    nindices_list = [2^20],
    clustersize_list = [50],
    include_uniform = true,
    include_clusters = true,
)
    problems = []
    for ncounters in ncounters_list, nindices in nindices_list
        label_suffix = "ncounters=$ncounters-nindices=$nindices"
        if include_uniform
            p = () -> random_increments_problem_uniform(ncounters, nindices)
            push!(problems, ("uniform-$label_suffix", p))
        end
        if include_clusters
            for clustersize in clustersize_list
                p = () -> random_increments_problem_clusters(ncounters, nindices, clustersize)
                push!(problems, ("cluster-clustersize=$clustersize-$label_suffix", p))
            end
        end
    end
    return problems
end

const CACHE = Ref{Any}(nothing)

function maybe_cached(f)
    c = CACHE[]
    if c isa Tuple{Any,Any} && c[1] === f
        return c[2]
    end
    t = @elapsed y = f()
    if t > 0.5
        CACHE[] = (f, y)
    end
    return y
end

function setup(;
    include_all = true,
    include_seq = include_all,
    include_atomics = include_all,
    include_tsx_static = include_all,
    include_tsx_adaptive = include_all,
    ntries_list = [10],
    nbatches_list = [50, 100],
    ntasks_list = default_ntasks_list(),
    problems = nothing,
    problems_options...,
)
    suite = BenchmarkGroup()

    if problems === nothing
        problems = default_problems(; problems_options...)
    end

    for (label, p) in problems
        s0 = suite[label] = BenchmarkGroup()

        if include_seq
            s0["seq"] = @benchmarkable(
                serial_random_increments!(counters, indices),
                setup = begin
                    counters, indices = maybe_cached($p)
                    fill!(counters, 0)
                end,
            )
        end
        for ntasks in ntasks_list
            s1 = s0["ntasks=$ntasks"] = BenchmarkGroup()
            if include_atomics
                s1["atomics"] = @benchmarkable(
                    parallel_random_increments!(counters, indices, $ntasks),
                    setup = begin
                        counters, indices = maybe_cached($p)
                        fill!(counters, 0)
                    end,
                )
            end
            if include_tsx_static
                s1["tsx_static"] = s2 = BenchmarkGroup()
                for ntries in ntries_list
                    s2["ntries=$ntries"] = s3 = BenchmarkGroup()
                    for nbatches in nbatches_list
                        s3["nbatches=$nbatches"] = @benchmarkable(
                            parallel_random_increments!(
                                counters,
                                indices,
                                $ntasks,
                                Val(true),
                                $nbatches,
                                $ntries,
                            ),
                            setup = begin
                                counters, indices = maybe_cached($p)
                                fill!(counters, 0)
                            end,
                        )
                    end
                end
            end
            if include_tsx_adaptive
                s1["tsx_adaptive"] = @benchmarkable(
                    parallel_random_increments_adaptive!(counters, indices, $ntasks),
                    setup = begin
                        counters, indices = maybe_cached($p)
                        fill!(counters, 0)
                    end,
                )
            end
        end
    end

    return suite
end

function clear()
    CACHE[] = nothing
end

end  # module


text:
# Base code for timestepping, precomputes a solver and tests it on a simple problem
push!(LOAD_PATH, string(pwd(),"/src"))

using Revise
using PyPlot
using ModStokesFastDirect
using Humanize
import ModifiedStokesSolver
import AnalyticDomains
import PUX

mss = ModifiedStokesSolver
include("../src/plottools.jl")
include("../src/timestepping.jl")

Solution = mss.Solution
PrecomputedSolver = mss.PrecomputedSolver
unitcircle = AnalyticDomains.starfish(n_arms=0, amplitude=0.0)

function estimate_fourier_resolution(U)
    Uhat = fftshift(fft(U))
    Ures = max( norm(Uhat[1:2,:], Inf),
                norm(Uhat[end-1:end,:], Inf),
                norm(Uhat[:,1:2], Inf),
                norm(Uhat[:,end-1:end], Inf) ) / maximum(abs.(Uhat))
    return Ures
end

function print_footprint(var)
    println("Footprint: ", datasize(Base.summarysize(var)))
end

# Precompute a solver
function precompute_solver(;curve=unitcircle, alpha=2.0,
                           numpanels=100, panelorder=16,
                           method=:heuristic, fastdirect_TOL=1e-10,
                           Ngrid=150, pux_ep=2, pux_R=0.4,
                           dotest=true)
    println("======= SETUP")
    @show alpha
    # Discretize
    println("* Discretizing")
    dcurve = mss.discretize(curve, numpanels, panelorder, equal_arclength=true)
    clf()
    plot_boundary(dcurve)
    # PUX prep
    println("* PUX precompute")
    @time PUXStor, X, Y, Lgrid, Ngrid, interior = PUX.pux_precompute(curve, dcurve, Ngrid, pux_ep, R=pux_R)
    print_footprint(PUXStor)
    figure()
    plot(dcurve.points[1,:], dcurve.points[2,:], ".k")
    for s=PUXStor
        PUX.plot_partitions(s)
    end
    axis("image")
    if method==:heuristic
        # Heuristic for small alpha and smooth problem
        # Other problems may change fastdirect build speed and GMRES convergence
        if sum(numpanels) <= 500
            method = :dense_gmres
        else
            method = :fastdirect
        end
        println("* Heuristic solver choice: $method")
    end

    if method==:fastdirect
        newmethod = :dense_gmres
        warn("Fast direct solver not available, switching to $newmethod")
        method = newmethod
    end

    if method==:fastdirect
        # Fast direct prep
        println("* Fast direct build")
        @time Ainv = compressed_system_matrix_inverse(dcurve, alpha, 64, fastdirect_TOL)
        print_footprint(Ainv)
    elseif method==:fmm
        # GMRES prep
        println("* GMRES FMM prep")
        @time flhs = mss.system_matvec(dcurve, alpha)
        print_footprint(flhs)
    elseif method==:dense_lu
        println("* Dense matrix assembly")
        @time MAT = mss.system_matrix(dcurve, alpha)
        print_footprint(MAT)
        println("* Dense matrix LU")
        @time MATL, MATU, MATp = lu(MAT)
    elseif method==:dense_gmres
        println("* Dense matrix assembly")
        @time MAT = mss.system_matrix(dcurve, alpha)
        print_footprint(MAT)
    elseif method==:none
        println("* Dummy solver")
        dotest = false
    else
        error("unknown method")
    end
    # Useful vars
    xbdry = dcurve.points[1,:]
    ybdry = dcurve.points[2,:]
    interior_points = [X[interior]' ; Y[interior]']
    numeval = size(interior_points, 2)
    if method == :none
        solve = nothing
    else
        println("* FMM tree precompute")
        @time fmm_precomp = mss.doublelayer_precomp(dcurve, interior_points, alpha)
        print_footprint(fmm_precomp)
        println("* Near eval precompute")
        @time near_weights = mss.doublelayer_near_weights(dcurve, interior_points, alpha; ifgrad=true)
        print_footprint(near_weights)
        # Form function that does the solve required in a timestep
        function solve(F1, F2, ubc1, ubc2)
            # Extend RHS
            #println("* PUX eval")
            tic()
            FEXT1, FEXT2 = PUX.pux_eval_vec(F1, F2, PUXStor)
            time_pux = toq()
            # Compute particular solution
            #println("* Periodic FFT Solve")
            tic()
            U1, U2, upbdry1, upbdry2, dU1dx, dU1dy, dU2dx, dU2dy =
                mss.per_modstokes(FEXT1, FEXT2, Lgrid, alpha, xbdry, ybdry, ifgrad=true)
            time_fft = toq()
            # Estimate particular solution resolution
            UP_resolution =  max(estimate_fourier_resolution(U1),
                                 estimate_fourier_resolution(U2))
            UPgrad_resolution = max( estimate_fourier_resolution(dU1dx),
                                     estimate_fourier_resolution(dU1dy),
                                     estimate_fourier_resolution(dU2dx),
                                     estimate_fourier_resolution(dU2dy) )
            # Integral equation RHS
            rhs1 = ubc1 - upbdry1
            rhs2 = ubc2 - upbdry2
            rhs = vec([rhs1'; rhs2'])
            # Solve integral equation
            tic()
            if method==:fastdirect
                #println("* Fast direct solve")
                dsol = apply_compressed_matrix(Ainv, rhs)
                density = reshape(dsol, 2, :)
            elseif method==:fmm
                #println("* GMRES FMM solve")
                density = mss.gmres_solve(dcurve, rhs, alpha, flhs=flhs)
            elseif method==:dense_lu
                #println("* LU apply")
                dsol = MATU\(MATL\rhs[MATp])
                density = reshape(dsol, 2, :)
            elseif method==:dense_gmres
                #println("* GMRES dense solve")
                dsol, gmlog = IterativeSolvers.gmres(MAT, rhs; tol=1e-15, restart=100, maxiter=100,
                verbose=false, log=true)
                if !gmlog.isconverged
                    warn("GMRES did not converge, residual = $(gmlog.data[:resnorm][end])")
                end
                density = reshape(dsol, 2, :)
            end
            time_solve = toq()
            # Estimate density resolution
            L = mss.legendre_matrix(panelorder)
            density_resolution = 0.0
            for i=1:dcurve.numpanels
                xcoeff = L*density[1, (1:panelorder) + panelorder*(i-1)]
                ycoeff = L*density[2, (1:panelorder) + panelorder*(i-1)]
                xres = max(abs(xcoeff[end]), abs(xcoeff[end-1]))
                yres = max(abs(ycoeff[end]), abs(ycoeff[end-1]))
                res = max(xres, yres)
                density_resolution = max(density_resolution, res)
            end
            # Evaluate layer potential at interior points
            #println("* FMM eval")
            tic()
            u, u1grad, u2grad = mss.doublelayer_fast(dcurve, density, interior_points, alpha; specquad=true, weights=near_weights, ifgrad=true, precomp=fmm_precomp)
            time_fmm = toq()
            @. begin
                # Add hom sol to part sol on grid
                U1[interior]   += u[1,:]
                U2[interior]   += u[2,:]
                dU1dx[interior] += u1grad[1,:]
                dU1dy[interior] += u1grad[2,:]
                dU2dx[interior] += u2grad[1,:]
                dU2dy[interior] += u2grad[2,:]
                # Clear exterior
                U1[!interior] = 0.0
                U2[!interior] = 0.0
                dU1dx[!interior] = 0.0
                dU1dy[!interior] = 0.0
                dU2dx[!interior] = 0.0
                dU2dy[!interior] = 0.0
            end
            println(@sprintf("Timings: pux=%.3f, fft=%.3f, solve=%.3f, fmm=%.3f", time_pux, time_fft, time_solve, time_fmm))
            info(@sprintf("Resolution est.: den=%.1e, uP=%.1e, uP'=%.1e",
                          density_resolution, UP_resolution, UPgrad_resolution))
            return Solution(U1, U2, dU1dx, dU1dy, dU2dx, dU2dy,
            copy(F1), copy(F2), copy(ubc1), copy(ubc2))
        end
    end
    solver = PrecomputedSolver(alpha, dcurve, X, Y, Lgrid, PUXStor, interior, solve)
    if dotest
        test_solver(solver)
    end
    return solver
end

# Useful functions
function pux_gradient(U1, U2, PUXStor, Lgrid)
    U1EXT, U2EXT = PUX.pux_eval_vec(U1, U2, PUXStor)
    dU1dx, dU1dy = mss.per_gradient(U1EXT, Lgrid)
    dU2dx, dU2dy = mss.per_gradient(U2EXT, Lgrid)
    return dU1dx, dU1dy, dU2dx, dU2dy
end
function ugradu(sol::Solution)
    ugradu(sol.U1, sol.U2, sol.dU1dx, sol.dU1dy, sol.dU2dx, sol.dU2dy)
end
function ugradu(U1, U2, dU1dx, dU1dy, dU2dx, dU2dy)
    F1 = U1.*dU1dx + U2.*dU1dy
    F2 = U1.*dU2dx + U2.*dU2dy
    return F1, F2
end

# Simple tester for precomputed solver
function test_solver(precomp)
    # Unpack
    dcurve = precomp.dcurve
    alpha = precomp.alpha
    X = precomp.X
    Y = precomp.Y
    interior = precomp.interior
    xbdry = dcurve.points[1,:]
    ybdry = dcurve.points[2,:]
    interior_points = [X[interior]' ; Y[interior]']
    numeval = size(interior_points, 2)
    # Define ref problem (linear oscillatory)
    k1 = 1.0
    k2 = 2.0
    ksq = k1^2 + k2^2
    phi_k(x, y) = exp(1im*(k1*x + k2*y))
    fhat1 = 3.0
    fhat2 = 5.0
    kdotfhat = k1*fhat1 + k2*fhat2
    # rhs
    ffunc1(x, y) = real(fhat1*phi_k(x,y))
    ffunc2(x, y) = real(fhat2*phi_k(x,y))
    # solution
    ufunc1(x, y) = real( (ksq*fhat1 - k1*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y) )
    ufunc2(x, y) = real( (ksq*fhat2 - k2*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y) )
    # derivatives of solution
    ufunc1x(x, y) = real( (ksq*fhat1 - k1*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k1 )
    ufunc1y(x, y) = real( (ksq*fhat1 - k1*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k2 )
    ufunc2x(x, y) = real( (ksq*fhat2 - k2*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k1 )
    ufunc2y(x, y) = real( (ksq*fhat2 - k2*kdotfhat)/(ksq*(ksq+alpha^2))*phi_k(x,y)*1im*k2 )

    # Setup test problem: BCs and RHS
    F1 = ffunc1.(X,Y)
    F2 = ffunc2.(X,Y)
    ubc1 = ufunc1.(xbdry, ybdry)
    ubc2 = ufunc2.(xbdry, ybdry)

    ### BEGIN STUFF TO PUT IN TIMESTEP LOOP
    println("======= TEST SOLVE")
    @time begin
        # Solve inhom mod Stokes
        sol = precomp.solve(F1, F2, ubc1, ubc2)
        newF1 = sol.U1.*sol.dU1dx + sol.U2.*sol.dU1dy
        newF2 = sol.U1.*sol.dU2dx + sol.U2.*sol.dU2dy
        println("======= TEST SOLVE DONE")
    end
    ### END STUFF TO PUT IN TIMESTEP LOOP
    # Check errors in test problem
    println("* Computing error")
    # Reference
    xint = interior_points[1,:]
    yint = interior_points[2,:]
    uref1 = ufunc1.(xint, yint)
    uref2 = ufunc2.(xint, yint)

    uref = [uref1 uref2]'
    unorm = norm(vec(uref), Inf)
    # Error
    u = [sol.U1[interior] sol.U2[interior]]'
    err = u .- uref
    maxnorm(v) = maximum(abs.(v))
    l2norm(v) = sqrt(sum(vec(v.^2)))
    println("Max rel err on grid: ", maxnorm(err)/maxnorm(uref))
    println("L2 rel err on grid: ", l2norm(err)/l2norm(uref) )

    println("* Derivative errors (max rel)")
    gradnorm = max( norm(sol.dU1dx[interior], Inf),
                    norm(sol.dU1dy[interior], Inf),
                    norm(sol.dU2dx[interior], Inf),
                    norm(sol.dU2dy[interior], Inf) )
    E_du1dx = norm(ufunc1x.(xint, yint) - sol.dU1dx[interior], Inf) / gradnorm
    E_du1dy = norm(ufunc1y.(xint, yint) - sol.dU1dy[interior], Inf) / gradnorm
    E_du2dx = norm(ufunc2x.(xint, yint) - sol.dU2dx[interior], Inf) / gradnorm
    E_du2dy = norm(ufunc2y.(xint, yint) - sol.dU2dy[interior], Inf) / gradnorm
    Egrad_max = max(E_du1dx, E_du1dy, E_du2dx, E_du2dy)
    @show Egrad_max

    println("* Displaying")
    clf()
    # Put on grid
    E = zeros(X)
    # Max error
    maxerr = maximum(abs.(err), 1)
    E[interior] = maxerr / unorm
    # Or vorticity
    vort_comp = sol.dU2dx[interior] - sol.dU1dy[interior]
    vort_ref = ufunc2x.(xint, yint) - ufunc1y.(xint, yint)
    E[interior] = abs.(vort_comp - vort_ref) / maxnorm(vort_ref)
    title("Vorticity error")

    logE = log10.(E+1e-100)
    #interior_pcolor(X, Y, logE, interior, vmin=-16, vmax=0)
    interior_pcolor(X, Y, logE, interior, vmax=0, cmap=PyPlot.cm[:coolwarm])
    cbar = colorbar()
    cbar[:set_label]("Rel. error")
    plot(dcurve.points[1,:], dcurve.points[2,:], ".k")
    for s=precomp.PUXStor
        PUX.plot_partitions(s)
    end

    axis("image")
    xmax = maximum(dcurve.points[1,:])*1.1
    xmin = minimum(dcurve.points[1,:])*1.1
    ymax = maximum(dcurve.points[2,:])*1.1
    ymin = minimum(dcurve.points[2,:])*1.1
    axis([xmin, xmax, ymin, ymax])
end

function timestep(Nstep, dt, sdc_order, hook, solver, t0, sol, ubc1, ubc2; istart=0)
    dt_substep = sdc_substep(dt, sdc_order)
    Re = solver.alpha^2*dt_substep
    ti = t0 + istart*dt
    for i=istart+1:Nstep
        println("=== STEP ", i, "/", Nstep)
        @time begin
            ti = t0 + i*dt
            newsol = step_sdc(dt_substep, sdc_order, sol, solver, ubc1, ubc2)
            prevsol, sol = sol, newsol
            hook(solver, sol, ti, i)
        end
    end
    return ti, sol
end

hook_null(a1, a2, a3, a4) = nothing

;


text:
################################################################################

"""
Keep track of scores & indexes.

# Fields
    score::Float64: score
    index::Int64: index

See also: [`HMM`](@ref), [`HMMParams`](@ref)
"""
struct ScorePair
  score::Float64
  index::Int64
end

"""
Hidden Markov model object.

# Fields
    data::Array: hold data
    tb::Array: collect traceback

See also: [`HMMParams`](@ref), [`setup`](@ref),
"""
mutable struct HMM
  data::Array
  tb::Array
end

# TODO: write a function signature verifier
"""
Hidden Markov model parameters.

# Fields
    penalty::Float64: penalty on hidden Markov model. Defualt = 200.
    minimumFrequency::Int64: minimun frequency on model state to split during hypothesis generation. Default = 20.
    verbosity::Bool: increase verbosity. Default = false.
    distance::Function: calculate distance. Possible options provided by this package inlcude: `euclDist`, `bhattDist`. Alternatively, a predefined or an anonymous / λ functions with the form `functionDist(arr::Array{T, 1}, h::Array{T, 1}) where T <: Number` can be passed.

See also: [`HMM`](@ref), [`setup`](@ref), [`euclDist`](@ref), [`bhattDist`](@ref),
"""
@with_kw struct HMMParams
  penalty::Float64           = 200.      # model penalty
  minimumFrequency::Int64    = 20        # minimum frequency for state split
  verbosity::Bool            = false     # verbosity level
  distance::Function         = distance  # distance function
end

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

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_forks_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), Apr 21, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames # load package

# ------------------------------------------------------------------------------------------
# ## Split-apply-combine
# ------------------------------------------------------------------------------------------

x = DataFrame(id=[1,2,3,4,1,2,3,4], id2=[1,2,1,2,1,2,1,2], v=rand(8))

gx1 = groupby(x, :id)

gx2 = groupby(x, [:id, :id2])

vcat(gx2...) # back to the original DataFrame

x = DataFrame(id = [missing, 5, 1, 3, missing], x = 1:5)

showall(groupby(x, :id)) # by default groups include mising values and are not sorted

showall(groupby(x, :id, sort=true, skipmissing=true)) # but we can change it :)

x = DataFrame(id=rand('a':'d', 100), v=rand(100));
by(x, :id, y->mean(y[:v])) # apply a function to each group of a data frame

by(x, :id, y->mean(y[:v]), sort=true) # we can sort the output

by(x, :id, y->DataFrame(res=mean(y[:v]))) # this way we can set a name for a column - DataFramesMeta @by is better

x = DataFrame(id=rand('a':'d', 100), x1=rand(100), x2=rand(100))
aggregate(x, :id, sum) # apply a function over all columns of a data frame in groups given by id

aggregate(x, :id, sum, sort=true) # also can be sorted

# ------------------------------------------------------------------------------------------
# *We omit the discussion of of map/combine as I do not find them very useful (better to use
# by)*
# ------------------------------------------------------------------------------------------

x = DataFrame(rand(3, 5))

map(mean, eachcol(x)) # map a function over each column and return a data frame

foreach(c -> println(c[1], ": ", mean(c[2])), eachcol(x)) # a raw iteration returns a tuple with column name and values

colwise(mean, x) # colwise is similar, but produces a vector

x[:id] = [1,1,2]
colwise(mean,groupby(x, :id)) # and works on GroupedDataFrame

map(r -> r[:x1]/r[:x2], eachrow(x)) # now the returned value is DataFrameRow which works similarly to a one-row DataFrame

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Printf
using PyPlot

using Oceananigans
using Oceananigans.Advection
using JULES

const km = 1000.0
const hPa = 100.0

L= 300km
H = 10km
Δ = 0.5km  # grid spacing [m]

Nx = Int(L/Δ)
Ny = 1
Nz = Int(H/Δ)

grid = RegularCartesianGrid(size=(Nx, Ny, Nz), halo=(4, 4, 4),
                            x=(0, L), y=(0, L), z=(0, H))

model = CompressibleModel(
                      grid = grid,
                     gases = DryEarth(),
    thermodynamic_variable = Entropy(),
                   closure = IsotropicDiffusivity(ν=75.0, κ=75.0)
)

#####
##### Dry thermal bubble perturbation
#####

gas = model.gases.ρ
R, cₚ, cᵥ = gas.R, gas.cₚ, gas.cᵥ
g  = model.gravity
pₛ = 1000hPa
θₛ = 300.0
N² = 1e-4
Γ  = - N² * θₛ / g
U  = 0

# Define an approximately hydrostatic background state
θ₀(x, y, z) = θₛ + Γ * z
p₀(x, y, z) = pₛ * (1 - g / (cₚ*Γ) * log(θ₀(x, y, z)/θₛ))^(cₚ/R)
T₀(x, y, z) = θₛ * (p₀(x, y, z)/pₛ)^(R/cₚ)
ρ₀(x, y, z) = p₀(x, y, z) / (R*T₀(x, y, z))

# Define both energy and entropy
uᵣ, Tᵣ, ρᵣ, sᵣ = gas.u₀, gas.T₀, gas.ρ₀, gas.s₀  # Reference values
ρe₀(x, y, z) = ρ₀(x, y, z) * (uᵣ + cᵥ * (T₀(x, y, z) - Tᵣ) + g*z)
ρs₀(x, y, z) = ρ₀(x, y, z) * (sᵣ + cᵥ * log(T₀(x, y, z)/Tᵣ) - R * log(ρ₀(x, y, z)/ρᵣ))

# Define the initial density perturbation
Δθ₀ = 1e-2
xᶜ = 100km
a = 5km

θ′(x, y, z) = Δθ₀ * sin(π*z/H) / (1 + (x-xᶜ)^2/a^2)
ρ′(x, y, z) = ρ₀(x, y, z) * θ′(x, y, z) / θ₀(x, y, z)

# Define initial state
ρᵢ(x, y, z) = ρ₀(x, y, z) + ρ′(x, y, z)
pᵢ(x, y, z) = p₀(x, y, z)
Tᵢ(x, y, z) = pᵢ(x, y, z) / (R * ρᵢ(x, y, z))

ρuᵢ(x, y, z) = ρᵢ(x, y, z) * U
ρeᵢ(x, y, z) = ρᵢ(x, y, z) * (uᵣ + cᵥ * (Tᵢ(x, y, z) - Tᵣ) + g*z)
ρsᵢ(x, y, z) = ρᵢ(x, y, z) * (sᵣ + cᵥ * log(Tᵢ(x, y, z)/Tᵣ) - R * log(ρᵢ(x, y, z)/ρᵣ))

# Set initial state (which includes the thermal perturbation)
set!(model.tracers.ρ, ρᵢ)
set!(model.momenta.ρu, ρuᵢ)
set!(model.tracers.ρs, ρsᵢ)
update_total_density!(model)

T = 1000
Δt = 0.1
Nt = Int(T/Δt)
for n in 1:Nt
    @info "iteration = $n/$Nt, ρ=$(model.total_density[200, 1, 10])"
    time_step!(model, Δt)
end


text:
function send_command_queue(command_queue::Vector{String})
    print.(command_queue)
    print("\n")
end

dock(ship::Ship, planet::Planet) = "d $(ship.id) $(planet.id)"
undock(ship::Ship)               = "u $(ship.id)"
thrust(ship::Ship, speed, angle) = "t $(ship.id) $(floor(Int, speed)) $(round(Int, angle))"

# function thrust(group::Group, speed, angle)
#     for ship in group.ships
#         thrust(ship::Ship, speed, angle)
#     end
# end

"""
navigate(game_map, ship, target)

Move a `ship` to a specific `target` position (Entity). It is recommended to place the position itself here, else navigate will crash into the target. If `avoid_obstacles` is set to true (default) will avoid obstacles on the way, with up to `max_corrections` corrections. Note that each correction accounts for `angular_step` degrees difference, meaning that the algorithm will naively try `max_correction` degrees before giving up (and returning empty string). The navigation will only consist of up to one command; call this method again in the next turn to continue navigating to the position.

# Arguments
- `game_map::GameMap`: The map of the game, from which obstacles will be extracted
- `ship::Ship`: ship to navigate
- `target::Entity`: The entity to which you will navigate
- `speed::Int`: The (max) speed to navigate. If the obstacle is nearer, will adjust accordingly.
- `avoid_obstacles::Bool`: Whether to avoid the obstacles in the way (simple pathfinding).
- `max_corrections::Int`: The maximum number of degrees to deviate per turn while trying to pathfind. If exceeded returns empty string.
- `angular_step::Float64`: The degree difference to deviate if the original destination has obstacles
- `ignore_ships::Bool`: Whether to ignore ships in calculations (this will make your movement faster, but more precarious)
- `ignore_planets::Bool`: Whether to ignore planets in calculations (useful if you want to crash onto planets)
"""
function navigate(game_map::GameMap, mobile, target;
            speed::Int = 0,
            avoid_obstacles = true,
            max_corrections::Int = 90,
            angular_step::Float64 = 1.0,
            ignore_ships = false,
            ignore_planets = false)
    max_corrections <= 0 && return ""

    dist = distance_between(mobile, target)
    angle = angle_between(mobile, target)
        if avoid_obstacles && !isempty(obstacles_between(game_map, mobile, target, ignore_ships, ignore_planets))
            new_target_dx = cos(deg2rad(angle + angular_step)) * dist
            new_target_dy = sin(deg2rad(angle + angular_step)) * dist
            new_target = Position(mobile.x + new_target_dx, mobile.y + new_target_dy)
            return navigate(game_map, mobile, new_target,
                            speed = speed,
                            avoid_obstacles = avoid_obstacles,
                            max_corrections = max_corrections - 1,
                            angular_step = angular_step,
                            ignore_ships = ignore_ships,
                            ignore_planets = ignore_planets)
        end
    speed = dist >= speed ? speed : floor(Int, dist)
    return thrust(mobile, speed, angle)
end

max_stars_repo_name: Cerbes/Halite-II

max_issues_repo_name: Cerbes/Halite-II

max_forks_repo_name: Cerbes/Halite-II


text:
using Gadgetron
using Test

@testset "ImageSerialization" begin

    meta = MRD.MetaDict()
    meta["penguin"] = 2.0
    meta["dinosaur"] = [2,3,4]
    img = MRD.Image(MRD.ImageHeader(version=4, acquisition_time_stamp=42), zeros(Float32, 2,3,4,1),meta)
    io = IOBuffer()
    MRD.write(io,img)
    seekstart(io)
    img3 = MRD.read(io,MRD.Image)

    @test img.header == img3.header
    @test img.data == img3.data
    @test img.meta == img3.meta
end

@testset "ImageSize" begin

    img = MRD.Image(MRD.ImageHeader(),zeros(Float32,1,1,1,1))
    header = MRD.RawImageHeader(img)

    io = IOBuffer()

    MRD.write(io,header)

    @test position(io) == 198
    meta = MRD.MetaDict()
    MRD.write(io,meta)
    @test position(io) == (198+8)

    io = IOBuffer()

end


text:
include("inner_problem.jl")

# qcqp version of relaxation: it's faster than the socp version
# please see the raw.jl file for the dual problem to this relaxation
function portfolios_socp(thePortfolio::SparsePortfolioData)

  n =thePortfolio.f
  m =thePortfolio.m
  f =thePortfolio.n

  socpProblem = Model(solver=MosekSolver(MSK_DPAR_INTPNT_QO_TOL_PFEAS=1e-5, MSK_DPAR_INTPNT_QO_TOL_DFEAS=1e-5, MSK_IPAR_LOG=0, MSK_IPAR_MAX_NUM_WARNINGS=0))

  @variable(socpProblem, α[1:n])
  @variable(socpProblem, λ)
  @variable(socpProblem, t>=0)
  @variable(socpProblem, v[1:f]>=0.0)
  @variable(socpProblem, w[1:f])
  @variable(socpProblem, βl[1:m]>=0.0)
  @variable(socpProblem, βu[1:m]>=0.0)
  @variable(socpProblem, ρ[1:n]>=0.0)

  @constraint(socpProblem, w.>=thePortfolio.X'*α+thePortfolio.A'*(βl-βu)+λ*ones(f)+ρ-thePortfolio.d)
  @constraint(socpProblem, knorm[i=1:f], v[i]+t+ρ[i]*thePortfolio.min_investment[i]>=(thePortfolio.γ[i]/2.0)*w[i]^2)

  @objective(socpProblem, Max, (-0.5*α'*α+thePortfolio.Y'*α+λ+βl'*thePortfolio.l-βu'*thePortfolio.u-sum(v)-thePortfolio.k*t)') #+π'*thePortfolio.min_investment

  status=solve(socpProblem)

  return Dual(getvalue(α), getvalue(λ), getvalue(βl), getvalue(βu), getvalue(ρ), getvalue(w), getobjectivevalue(socpProblem), status), getobjectivevalue(socpProblem)
end

# SOCP version of relaxation. Not in use, as this version is currently slower.
function portfolios_socp2(thePortfolio::SparsePortfolioData)

  n = size(thePortfolio.X, 1)
  m = size(thePortfolio.A, 1)
  f = size(thePortfolio.X, 2)

  socpProblem = Model(solver=MosekSolver( MSK_DPAR_INTPNT_CO_TOL_PFEAS=1e-8, MSK_DPAR_INTPNT_CO_TOL_DFEAS=1e-8))

  @variable(socpProblem, α[1:n])
  @variable(socpProblem, λ)
  @variable(socpProblem, t>=0)
  @variable(socpProblem, v[1:f]>=0.0)
  @variable(socpProblem, w[1:f])
  @variable(socpProblem, βl[1:m]>=0.0)
  @variable(socpProblem, βu[1:m]>=0.0)
  @variable(socpProblem, ρ[1:n]>=0.0)
  @variable(socpProblem, τ>=0.0)

  @constraint(socpProblem, w.>=thePortfolio.X'*α+thePortfolio.A'*(βl-βu)+λ*ones(f)+ρ-thePortfolio.d)
  @constraint(socpProblem, knorm[i=1:f], v[i]+t+ρ[i]*thePortfolio.min_investment[i]+1.0>=norm([sqrt(2.0*thePortfolio.γ[i])*w[i];v[i]+t+ρ[i]*thePortfolio.min_investment[i]-1.0]))
  @constraint(socpProblem, τ+1.0>=norm([2.0*α;τ-1.0]))

  @objective(socpProblem, Max, (-0.5*τ+thePortfolio.Y'*α+λ+βl'*thePortfolio.l-βu'*thePortfolio.u-sum(v)-thePortfolio.k*t)') #+π'*thePortfolio.min_investment

  status=solve(socpProblem)

  return Dual(getvalue(α), getvalue(λ), getvalue(βl), getvalue(βu), getvalue(ρ), getvalue(w), status, getobjectivevalue(socpProblem)), getobjectivevalue(socpProblem)
end


text:
function closer(ps::ParseState)
    (ps.closer.newline && kindof(ps.ws) == NewLineWS && kindof(ps.t) != Tokens.COMMA) ||
    (ps.closer.semicolon && kindof(ps.ws) == SemiColonWS) ||
    (isoperator(ps.nt) && precedence(ps.nt) <= ps.closer.precedence) ||
    (kindof(ps.nt) == Tokens.WHERE && ps.closer.precedence == LazyAndOp) ||
    (ps.closer.inwhere && kindof(ps.nt) == Tokens.WHERE) ||
    (ps.closer.inwhere && ps.closer.ws && kindof(ps.t) == Tokens.RPAREN && isoperator(ps.nt) && precedence(ps.nt) < DeclarationOp) ||
    (ps.closer.precedence > WhereOp && (
        (kindof(ps.nt) == Tokens.LPAREN && !(ps.t.kind === Tokens.EX_OR)) ||
        kindof(ps.nt) == Tokens.LBRACE ||
        kindof(ps.nt) == Tokens.LSQUARE ||
        (kindof(ps.nt) == Tokens.STRING && isemptyws(ps.ws)) ||
        ((kindof(ps.nt) == Tokens.RPAREN || kindof(ps.nt) == Tokens.RSQUARE) && isidentifier(ps.nt))
    )) ||
    (kindof(ps.nt) == Tokens.COMMA && ps.closer.precedence > 0) ||
    kindof(ps.nt) == Tokens.ENDMARKER ||
    (ps.closer.comma && iscomma(ps.nt)) ||
    (ps.closer.tuple && (iscomma(ps.nt) || isassignment(ps.nt))) ||
    (kindof(ps.nt) == Tokens.FOR && ps.closer.precedence > -1) ||
    (ps.closer.block && kindof(ps.nt) == Tokens.END) ||
    (ps.closer.paren && kindof(ps.nt) == Tokens.RPAREN) ||
    (ps.closer.brace && kindof(ps.nt) == Tokens.RBRACE) ||
    (ps.closer.square && kindof(ps.nt) == Tokens.RSQUARE) ||
    kindof(ps.nt) == Tokens.ELSEIF ||
    kindof(ps.nt) == Tokens.ELSE ||
    kindof(ps.nt) == Tokens.CATCH ||
    kindof(ps.nt) == Tokens.FINALLY ||
    (ps.closer.ifop && isoperator(ps.nt) && (precedence(ps.nt) <= 0 || kindof(ps.nt) == Tokens.COLON)) ||
    (ps.closer.range && (kindof(ps.nt) == Tokens.FOR || iscomma(ps.nt) || kindof(ps.nt) == Tokens.IF)) ||
    (ps.closer.ws && !isemptyws(ps.ws) &&
        !(kindof(ps.nt) == Tokens.COMMA) &&
        !(kindof(ps.t) == Tokens.COMMA) &&
        !(!ps.closer.inmacro && kindof(ps.nt) == Tokens.FOR) &&
        !(kindof(ps.nt) == Tokens.DO) &&
        !(
            (isbinaryop(ps.nt) && !(isemptyws(ps.nws) && isunaryop(ps.nt) && ps.closer.wsop)) ||
            (isunaryop(ps.t) && kindof(ps.ws) == WS)
        )) ||
    (ps.closer.unary && (kindof(ps.t) in (Tokens.INTEGER, Tokens.FLOAT, Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE) && kindof(ps.nt) == Tokens.IDENTIFIER))
end

"""
    @closer ps rule body

Continues parsing closing on `rule`.
"""
macro closer(ps, opt, body)
    quote
        local tmp1 = getfield($(esc(ps)).closer, $opt)
        setfield!($(esc(ps)).closer, $opt, true)
        out = $(esc(body))
        setfield!($(esc(ps)).closer, $opt, tmp1)
        out
    end
end

"""
    @nocloser ps rule body

Continues parsing not closing on `rule`.
"""
macro nocloser(ps, opt, body)
    quote
        local tmp1 = getfield($(esc(ps)).closer, $opt)
        setfield!($(esc(ps)).closer, $opt, false)
        out = $(esc(body))
        setfield!($(esc(ps)).closer, $opt, tmp1)
        out
    end
end

macro closeparen(ps, body)
    quote
        local tmp1 = $(esc(ps)).closer.paren
        $(esc(ps)).closer.paren = true
        push!($(esc(ps)).closer.cc, :paren)
        out = $(esc(body))
        pop!($(esc(ps)).closer.cc)
        $(esc(ps)).closer.paren = tmp1
        out
    end
end

macro closesquare(ps, body)
    quote
        local tmp1 = $(esc(ps)).closer.square
        $(esc(ps)).closer.square = true
        push!($(esc(ps)).closer.cc, :square)
        out = $(esc(body))
        pop!($(esc(ps)).closer.cc)
        $(esc(ps)).closer.square = tmp1
        out
    end
end
macro closebrace(ps, body)
    quote
        local tmp1 = $(esc(ps)).closer.brace
        $(esc(ps)).closer.brace = true
        push!($(esc(ps)).closer.cc, :brace)
        out = $(esc(body))
        pop!($(esc(ps)).closer.cc)
        $(esc(ps)).closer.brace = tmp1
        out
    end
end

"""
    @precedence ps prec body

Continues parsing binary operators until it hits a more loosely binding
operator (with precdence lower than `prec`).
"""
macro precedence(ps, prec, body)
    quote
        local tmp1 = $(esc(ps)).closer.precedence
        $(esc(ps)).closer.precedence = $(esc(prec))
        out = $(esc(body))
        $(esc(ps)).closer.precedence = tmp1
        out
    end
end

# Closer_TMP and ancillary functions help reduce code generation
struct Closer_TMP
    newline::Bool
    semicolon::Bool
    inmacro::Bool
    tuple::Bool
    comma::Bool
    insquare::Bool
    range::Bool
    ifop::Bool
    ws::Bool
    wsop::Bool
    unary::Bool
    precedence::Int
end

@noinline function create_tmp(c::Closer)
    Closer_TMP(c.newline,
        c.semicolon,
        c.inmacro,
        c.tuple,
        c.comma,
        c.insquare,
        c.range,
        c.ifop,
        c.ws,
        c.wsop,
        c.unary,
        c.precedence)
end

@noinline function update_from_tmp!(c::Closer, tmp::Closer_TMP)
    c.newline = tmp.newline
    c.semicolon = tmp.semicolon
    c.inmacro = tmp.inmacro
    c.tuple = tmp.tuple
    c.comma = tmp.comma
    c.insquare = tmp.insquare
    c.range = tmp.range
    c.ifop = tmp.ifop
    c.ws = tmp.ws
    c.wsop = tmp.wsop
    c.unary = tmp.unary
    c.precedence = tmp.precedence
end

@noinline function update_to_default!(c::Closer)
    c.newline = true
    c.semicolon = true
    c.inmacro = false
    c.tuple = false
    c.comma = false
    c.insquare = false
    c.range = false
    c.ifop = false
    c.ws = false
    c.wsop = false
    c.unary = false
    c.precedence = -1
end

"""
    @default ps body

Parses the next expression using default closure rules.
"""
macro default(ps, body)
    quote
        TMP = create_tmp($(esc(ps)).closer)
        update_to_default!($(esc(ps)).closer)
        out = $(esc(body))
        update_from_tmp!($(esc(ps)).closer, TMP)
        out
    end
end

isidentifier(t::AbstractToken) = kindof(t) == Tokens.IDENTIFIER
isidentifier(x::EXPR) = typof(x) === IDENTIFIER

isliteral(t::AbstractToken) = Tokens.begin_literal < kindof(t) < Tokens.end_literal
isliteral(x::EXPR) = typof(x) === LITERAL

isbool(t::AbstractToken) =  Tokens.TRUE ≤ kindof(t) ≤ Tokens.FALSE
iscomma(t::AbstractToken) =  kindof(t) == Tokens.COMMA

iskw(t::AbstractToken) = Tokens.iskeyword(kindof(t))
iskw(x::EXPR) = typof(x) === KEYWORD

isinstance(t::AbstractToken) = isidentifier(t) ||
                       isliteral(t) ||
                       isbool(t) ||
                       iskw(t)

ispunctuation(t::AbstractToken) = kindof(t) == Tokens.COMMA ||
                            kindof(t) == Tokens.END ||
                            Tokens.LSQUARE ≤ kindof(t) ≤ Tokens.RPAREN ||
                            kindof(t) == Tokens.AT_SIGN
ispunctuation(x::EXPR) = typof(x) === PUNCTUATION

isstring(x) = typof(x) === StringH || (isliteral(x) && (kindof(x) == Tokens.STRING || kindof(x) == Tokens.TRIPLE_STRING))
is_integer(x) = isliteral(x) && kindof(x) == Tokens.INTEGER
is_float(x) = isliteral(x) && kindof(x) == Tokens.FLOAT
is_number(x) = isliteral(x) && (kindof(x) == Tokens.INTEGER || kindof(x) == Tokens.FLOAT)
is_nothing(x) = isliteral(x) && kindof(x) == Tokens.NOTHING

isajuxtaposition(ps::ParseState, ret::EXPR) = ((is_number(ret) && (kindof(ps.nt) == Tokens.IDENTIFIER || kindof(ps.nt) == Tokens.LPAREN || kindof(ps.nt) == Tokens.CMD || kindof(ps.nt) == Tokens.STRING || kindof(ps.nt) == Tokens.TRIPLE_STRING)) ||
        ((typof(ret) === UnaryOpCall && is_prime(ret.args[2]) && kindof(ps.nt) == Tokens.IDENTIFIER) ||
        ((kindof(ps.t) == Tokens.RPAREN || kindof(ps.t) == Tokens.RSQUARE) && (kindof(ps.nt) == Tokens.IDENTIFIER || kindof(ps.nt) == Tokens.CMD)) ||
        ((kindof(ps.t) == Tokens.STRING || kindof(ps.t) == Tokens.TRIPLE_STRING) && (kindof(ps.nt) == Tokens.STRING || kindof(ps.nt) == Tokens.TRIPLE_STRING)))) || ((kindof(ps.t) in (Tokens.INTEGER, Tokens.FLOAT) || kindof(ps.t) in (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE)) && kindof(ps.nt) == Tokens.IDENTIFIER)

# When using the FancyDiagnostics package, Base.parse, is the
# same as CSTParser.parse. Manually call the flisp parser here
# to make sure we test what we want, even when people load the
# FancyDiagnostics package.
function flisp_parse(str::AbstractString, pos::Int; greedy::Bool = true, raise::Bool = true)
    # pos is one based byte offset.
    # returns (expr, end_pos). expr is () in case of parse error.
    bstr = String(str)
    ex, pos = ccall(:jl_parse_string, Any,
                    (Ptr{UInt8}, Csize_t, Int32, Int32),
                    bstr, sizeof(bstr), pos - 1, greedy ? 1 : 0)
    if raise && isa(ex, Expr) && ex.head === :error
        throw(Meta.ParseError(ex.args[1]))
    end
    if ex === ()
        raise && throw(Meta.ParseError("end of input"))
        ex = Expr(:error, "end of input")
    end
    return ex, pos + 1 # C is zero-based, Julia is 1-based
end

function flisp_parse(str::AbstractString; raise::Bool = true)
    ex, pos = flisp_parse(str, 1, greedy = true, raise = raise)
    if isa(ex, Expr) && ex.head === :error
        return ex
    end
    if !(pos > ncodeunits(str))
        raise && throw(Meta.ParseError("extra token after end of expression"))
        return Expr(:error, "extra token after end of expression")
    end
    return ex
end

function flisp_parse(stream::IO; greedy::Bool = true, raise::Bool = true)
    pos = position(stream)
    ex, Δ = flisp_parse(read(stream, String), 1, greedy = greedy, raise = raise)
    seek(stream, pos + Δ - 1)
    return ex
end

using Base.Meta
norm_ast(a::Any) = begin
    if isa(a, Expr)
        for (i, arg) in enumerate(a.args)
            a.args[i] = norm_ast(arg)
        end
        if a.head === :line
            return Expr(:line, a.args[1], :none)
        end
        if a.head === :macrocall
            fa = a.args[1]
            if fa === Symbol("@int128_str")
                return Base.parse(Int128, a.args[3])
            elseif fa === Symbol("@uint128_str")
                return Base.parse(UInt128, a.args[3])
            elseif fa === Symbol("@bigint_str")
                return  Base.parse(BigInt, a.args[3])
            elseif fa == Symbol("@big_str")
                s = a.args[3]
                n = tryparse(BigInt, s)
                if !(n === nothing)
                    return (n)
                end
                n = tryparse(BigFloat, s)
                if !(n === nothing)
                    return isnan((n)) ? :NaN : (n)
                end
                return s
            end
        elseif length(a.args) >= 2 && isexpr(a, :call) && a.args[1] == :- && isa(a.args[2], Number)
            return -a.args[2]
        end
        return a
    elseif isa(a, QuoteNode)
        return Expr(:quote, norm_ast(a.value))
    elseif isa(a, AbstractFloat) && isnan(a)
        return :NaN
    end
    return a
end

function flisp_parsefile(str, display = true)
    io = IOBuffer(str)
    failed = false
    x1 = Expr(:file)
    try
        while !eof(io)
            push!(x1.args, flisp_parse(io))
        end
    catch er
        isa(er, InterruptException) && rethrow(er)
        if display
            Base.showerror(stdout, er, catch_backtrace())
            println()
        end
        return x1, true
    end
    if length(x1.args) > 0  && x1.args[end] === nothing
        pop!(x1.args)
    end
    x1 = norm_ast(x1)
    remlineinfo!(x1)
    return x1, false
end

function cst_parsefile(str)
    x, ps = CSTParser.parse(ParseState(str), true)
    sp = check_span(x)
    # remove leading/trailing nothings
    if length(x.args) > 0 && is_nothing(x.args[1])
        popfirst!(x.args)
    end
    if length(x.args) > 0 && is_nothing(x.args[end])
        pop!(x.args)
    end
    x0 = norm_ast(Expr(x))
    x0, ps.errored, sp
end

function check_file(file, ret, neq)
    str = read(file, String)
    x0, cstfailed, sp = cst_parsefile(str)
    x1, flispfailed = flisp_parsefile(str)

    print("\r                             ")
    if !isempty(sp)
        printstyled(file, color = :blue)
        @show sp
        println()
        push!(ret, (file, :span))
    end
    if cstfailed
        printstyled(file, color = :yellow)
        println()
        push!(ret, (file, :errored))
    elseif !(x0 == x1)
        cumfail = 0
        printstyled(file, color = :green)
        println()
        c0, c1 = CSTParser.compare(x0, x1)
        printstyled(string("    ", c0), bold = true, color = :ligth_red)
        println()
        printstyled(string("    ", c1), bold = true, color = :light_green)
        println()
        push!(ret, (file, :noteq))
    end
end

function check_base(dir = dirname(Base.find_source_file("essentials.jl")), display = false)
    N = 0
    neq = 0
    err = 0
    aerr = 0
    fail = 0
    bfail = 0
    ret = []
    oldstderr = stderr
    redirect_stderr()
    for (rp, d, files) in walkdir(dir)
        for f in files
            file = joinpath(rp, f)
            if endswith(file, ".jl")
                N += 1
                try
                    print("\r", rpad(string(N), 5), rpad(string(round(fail / N * 100, sigdigits = 3)), 8), rpad(string(round(err / N * 100, sigdigits = 3)), 8), rpad(string(round(neq / N * 100, sigdigits = 3)), 8))

                    check_file(file, ret, neq)
                catch er
                    isa(er, InterruptException) && rethrow(er)
                    if display
                        Base.showerror(stdout, er, catch_backtrace())
                        println()
                    end
                    fail += 1
                    printstyled(file, color = :red)
                    println()
                    push!(ret, (file, :failed))
                end
            end
        end
    end
    redirect_stderr(oldstderr)
    if bfail + fail + err + neq > 0
        println("\r$N files")
        printstyled("failed", color = :red)
        println(" : $fail    $(100 * fail / N)%")
        printstyled("errored", color = :yellow)
        println(" : $err     $(100 * err / N)%")
        printstyled("not eq.", color = :green)
        println(" : $neq    $(100 * neq / N)%", "  -  $aerr     $(100 * aerr / N)%")
        printstyled("base failed", color = :magenta)
        println(" : $bfail    $(100 * bfail / N)%")
    end
    ret
end

"""
    compare(x,y)

Recursively checks whether two Base.Expr are the same. Returns unequal sub-
expressions.
"""
compare(x, y) = x == y ? true : (x, y)

function compare(x::Expr, y::Expr)
    if x == y
        return true
    else
        if x.head != y.head
            return (x, y)
        end
        if length(x.args) != length(y.args)
            return (x.args, y.args)
        end
        for i = 1:length(x.args)
            t = compare(x.args[i], y.args[i])
            if t != true
                return t
            end
        end
    end
end

"""
check_span(x, neq = [])

Recursively checks whether the span of an expression equals the sum of the span
of its components. Returns a vector of failing expressions.
"""
function check_span(x::EXPR, neq = [])
    (ispunctuation(x) || isidentifier(x) || iskw(x) || isoperator(x) || isliteral(x) || typof(x) == StringH) && return neq

    s = 0
    for a in x.args
        check_span(a, neq)
        s += a.fullspan
    end
    if length(x.args) > 0 && s != x.fullspan
        push!(neq, x)
    end
    neq
end

function speed_test()
    dir = dirname(Base.find_source_file("essentials.jl"))
    println("speed test : ", @timed(for i = 1:5
        parse(read(joinpath(dir, "essentials.jl"), String), true);
        parse(read(joinpath(dir, "abstractarray.jl"), String), true);
    end)[2])
end

"""
    check_reformat()

Reads and parses all files in current directory, applys formatting fixes and checks that the output AST remains the same.
"""
function check_reformat()
    fs = filter(f->endswith(f, ".jl"), readdir())
    for (i, f) in enumerate(fs)
        f == "deprecated.jl" && continue
        str = read(f, String)
        x, ps = parse(ParseState(str), true);
        cnt = 0
        for i = 1:length(x) - 1
            y = x[i]
            sstr = str[cnt + (1:y.span)]

            y1 = parse(sstr)
            @assert Expr(y) == Expr(y1)
            cnt += y.span
        end
    end
end

function trailing_ws_length(x)
    x.fullspan - x.span
end

Base.iterate(x::EXPR) = length(x) == 0 ? nothing : (x.args[1], 1)
Base.iterate(x::EXPR, s) = s < length(x) ? (x.args[s + 1], s + 1) : nothing
Base.length(x::EXPR) = x.args isa Nothing ? 0 : length(x.args)

@inline val(token::RawToken, ps::ParseState) = String(ps.l.io.data[token.startbyte + 1:token.endbyte + 1])

function str_value(x)
    if isidentifier(x) || typof(x) === LITERAL
        return valof(x)
    elseif typof(x) === OPERATOR || typof(x) === MacroName
        return string(Expr(x))
    else
        return ""
    end
end

_unescape_string(s::AbstractString) = sprint(_unescape_string, s, sizehint = lastindex(s))
function _unescape_string(io, s::AbstractString)
    a = Iterators.Stateful(s)
    for c in a
        if !isempty(a) && c == '\\'
            c = popfirst!(a)
            if c == 'x' || c == 'u' || c == 'U'
                n = k = 0
                m = c == 'x' ? 2 :
                    c == 'u' ? 4 : 8
                while (k += 1) <= m && !isempty(a)
                    nc = Base.peek(a)
                    n = '0' <= nc <= '9' ? n << 4 + nc - '0' :
                        'a' <= nc <= 'f' ? n << 4 + nc - 'a' + 10 :
                        'A' <= nc <= 'F' ? n << 4 + nc - 'A' + 10 : break
                    popfirst!(a)
                end
                if k == 1
                    # throw(ArgumentError("invalid $(m == 2 ? "hex (\\x)" :
                    #                         "unicode (\\u)") escape sequence used in $(repr(s))"))
                    # push error to ParseState?
                    n = 0
                end
                if m == 2 # \x escape sequence
                    write(io, UInt8(n))
                else
                    print(io, Char(n))
                end
            elseif '0' <= c <= '7'
                k = 1
                n = c - '0'
                while (k += 1) <= 3 && !isempty(a)
                    c  = Base.peek(a)
                    n = ('0' <= c <= '7') ? n << 3 + c - '0' : break
                    popfirst!(a)
                end
                if n > 255
                    # throw(ArgumentError("octal escape sequence out of range"))
                    # push error to ParseState?
                    n = 255
                end
                write(io, UInt8(n))
            else
                print(io, c == 'a' ? '\a' :
                          c == 'b' ? '\b' :
                          c == 't' ? '\t' :
                          c == 'n' ? '\n' :
                          c == 'v' ? '\v' :
                          c == 'f' ? '\f' :
                          c == 'r' ? '\r' :
                          c == 'e' ? '\e' : c)
            end
        else
            print(io, c)
        end
    end
end

function match_closer(ps::ParseState)
    length(ps.closer.cc) == 0 && return false
    kind = kindof(ps.nt)
    lc = last(ps.closer.cc)
    return (kind === Tokens.RPAREN && lc == :paren) ||
           (kind === Tokens.RSQUARE && lc == :square) ||
           (kind === Tokens.RBRACE && lc == :braces) ||
           (kind === Tokens.END && (lc == :begin || lc == :if))
end


text:
""" [vanZanten2011](@cite) """
function Rico_u(::Type{FT}) where {FT}
    return z -> FT(-9.9) + FT(2.0e-3) * z
end

""" [vanZanten2011](@cite) """
function Rico_v(::Type{FT}) where {FT}
    return z -> FT(-3.8)
end

""" [vanZanten2011](@cite) """
Rico_θ_liq_ice(::Type{FT}) where {FT} =
    z -> if z <= 740.0
        FT(297.9)
    else
        FT(297.9) + (317 - FT(297.9)) / (4000 - 740) * (z - 740)
    end

""" [vanZanten2011](@cite) """
Rico_q_tot(::Type{FT}) where {FT} =
    z -> if z <= 740.0
        (16 + (FT(13.8) - 16) / 740 * z) / 1000
    elseif z > 740.0 && z <= 3260.0
        (FT(13.8) + (FT(2.4) - FT(13.8)) / (3260 - 740) * (z - 740)) / 1000
    else
        (FT(2.4) + (FT(1.8) - FT(2.4)) / (4000 - 3260) * (z - 3260)) / 1000
    end

""" [vanZanten2011](@cite) """
function Rico_geostrophic_ug(::Type{FT}) where {FT}
    return z -> FT(-9.9) + FT(2.0e-3) * z
end

""" [vanZanten2011](@cite) """
function Rico_geostrophic_vg(::Type{FT}) where {FT}
    return z -> FT(-3.8)
end

""" [vanZanten2011](@cite) """
function Rico_dTdt(::Type{FT}) where {FT}
    return (Π, z) -> (FT(-2.5) / (3600 * 24)) * Π
end

""" [vanZanten2011](@cite) """
function Rico_dqtdt(::Type{FT}) where {FT}
    return z -> if z <= 2980.0
        (-1 + FT(1.3456) / FT(2980.0) * z) / FT(86400.0) / 1000   #kg/(kg * s)
    else
        FT(0.3456) / 86400 / 1000
    end
end

""" [vanZanten2011](@cite) """
function Rico_subsidence(::Type{FT}) where {FT}
    return z -> if z <= 2260.0
            -(FT(0.005) / FT(2260.0)) * z
        else
            FT(-0.005)
        end
end


text:
#=
TODO: consider implementing Algorithm to obtain M-Δ Form for Robust Control
Lt Mary K Manningt and Siva S Banda$
=#

"""
    δ(N=32)

Create an uncertain element of `N` uniformly distributed samples ∈ [-1, 1]
"""
δ(N=32) = StaticParticles(N, Uniform(-1, 1))

function ControlSystems.lft(M11::AbstractMatrix, M12::AbstractMatrix, M21::AbstractMatrix, M22::AbstractMatrix, d::AbstractMatrix)
    M11 + M12 * d*((I - M22*d)\M21)
end

wass(a, b) = wasserstein(a,b,2)
l2(a,b) = sum(abs2.(a.particles .- b.particles))

wass(a::AbstractArray, b::AbstractArray) = sum(wass(a,b) for (a,b) in zip(a,b)) # shouldn't be allowed to permute all particles individually?
l2(a::AbstractArray, b::AbstractArray) = sum(l2(a,b) for (a,b) in zip(a,b)) # shouldn't be allowed to permute all particles individually?

function pars2lft(M11, p, δ)
    @unpack M12, M21, M22 = p
    Phat = lft(M11, M12, M21, M22, δ)
end

struct LFT
    M11
    M12
    M21
    M22
    d
end

Base.Matrix(l::LFT) = [l.M11 l.M12; l.M21 l.M22]

function Base.show(io::IO, l::LFT)
    display(Matrix(l))
end

function ControlSystems.lft(l::LFT)
    @unpack M11, M12, M21, M22, d = l
    lft(M11, M12, M21, M22, d)
end

function ControlSystems.ss(l::LFT, P0::AbstractStateSpace)
    P = lft(l)
    A = P[:, 1:P0.nx]
    B = P[:, P0.nx+1:end]
    ss(A, B, P0.C, P0.D)
end

find_lft(sys::StateSpace{<:Any, <:StaticParticles{<:Any, N}}, n_uncertain::Int)  where N = find_lft(sys, [δ(N) for _ in 1:n_uncertain], wass)

"""
    l = find_lft(sys::StateSpace{<:Any, <:StaticParticles{<:Any, N}}, δ) where N

Given an systems `sys` with uncertain coefficients in the form of `StaticParticles`, find a lower linear fractional transformation `M` such that `lft(M, δ) ≈ sys`.

`δ` can be either the source of uncertainty in `sys`, i.e., a vector of the unique uncertain parameters that were used to create `sys`. These should be constructed as uniform randomly distributed particles for most robust-control theory to be applicable.
`δ` can also be an integer, in which case a numer of `δ` sources of uncertainty are automatically created. This could be used for order reduction if the number of uncertainty sources in `sys` is large.

Note, uncertainty in `sys` is only supported in `A` and `B`, `C` and `D` must be deterministic.

Returns `l::LFT` that internaly contains all four blocks of `M` as well as `δ`. Call `ss(l,sys)` do obtain `lft(M, δ) ≈ sys`.

Call `Matrix(l)` to obtain `M = [M11 M12; M21 M22]`
"""
function find_lft(sys::StateSpace{<:Any, <:AbstractParticles{<:Any, N}}, delta, dist::F = l2) where {N, F}
    n_uncertain = length(delta)
    P = [sys.A sys.B]
    M11 = pmean.(P)
    δ = Diagonal(delta)
    function loss(p)
        Phat = pars2lft(M11, p, δ)
        dist(P, Phat) #+ 1e-5*sum(abs, p)
    end
    p0 = ComponentVector(
        M12 = 0.00001randn(sys.nx, n_uncertain),
        M21 = 0.00001randn(n_uncertain, sys.nx+sys.nu),
        M22 = 0.00001randn(n_uncertain, n_uncertain),
    )
    res = Optim.optimize(
        loss,
        p0,
        BFGS(),
        Optim.Options(
            store_trace       = false,
            show_trace        = true,
            show_every        = 2,
            iterations        = 40000,
            allow_f_increases = true,
            time_limit        = 45,
            x_tol             = 1e-8,
            f_reltol          = 0,
            g_tol             = 1e-16,
        ),
    )
    p = res.minimizer
    LFT(M11, p.M12, p.M21, p.M22, δ), res

end


text:
# SMS EMOA stands for S Metric Selection Evolutionary Multi-Objective Algorithm
# S metric is the measure of hypervolume
# Simply put, each iteration, the algorithm generate a random feasible offspring and insert it in the popultion,
# and then removes the individual that is the less contributing to the overall hypervolume.
# Thus: |P_k| = |P_{k+1}| (stationary population)
# The algorithm tend to focus more on the "eye" of the pareto front.
# source: https://doi.org/10.1007/978-3-540-31880-4_5

function SMS_EMOA_update(P::Vector{Sol}, config::Config)

    debug && DEBUG_feasible_solutions(P)
    debug && DEBUG_correct_evaluations(P, hv = true)

    N = length(P)

    #--------- generate offspring ----------#

    offspring_has_been_generated = false
    offspring = 0
    while !offspring_has_been_generated

        # choose parents
        parent_1 = SELECTION_binary_tournament(P, criteria = HV)
        parent_2 = SELECTION_binary_tournament(P, criteria = HV)
        while parent_2 == parent_1
            parent_2 = SELECTION_binary_tournament(P, criteria = HV)
        end

        # recombine parents
        offspring = config.crossover(parent_1, parent_2)

        # mutate offspring
        if rand() < config.mutation_rate
            MUTATION_flip!(offspring)
        end

        # repair offspring
        config.reparation!(offspring)

        # check if the offspring is new to the population
        offspring_has_been_generated = true
        for sol in P
            if sol.z == offspring.z
                offspring_has_been_generated = false
                break
            end
        end
    end

    #-------------- add offspring to current population ---------------#

    # add offspring to the population
    push!(P, offspring)

    #- sort population by rank, and compute hypervolume for each rank -#

    # evaluate population
    F = INDICATOR_fast_non_dominated_sort!(P)
    for front in F
        INDICATOR_hypervolume_contribution!(front)
    end

    R_I = F[length(F)]

    # find the solution in R_I that contributes the less to the hypervolume
    idx_min = 1
    sol_min = R_I[1]
    hv_min = R_I[1].hv_contrib
    for idx in 1:length(R_I)
        sol = R_I[idx]
        if sol.hv_contrib < hv_min
            idx_min = idx
            sol_min = sol
            hv_min = sol.hv_contrib
        end
    end

    deleteat!(R_I, idx_min) # remove the worst contributing solution in the front

    # evaluate the new contributions (for next iteration)
    if length(R_I) > 0
        INDICATOR_hypervolume_contribution!(R_I)
    end

    # mark the solution that contributes the less to the hypervolume
    sol_min.hv_contrib = -1

    # find the solution marked in P
    idx_worst_contribution_sol = 0
    idx = 1
    while idx <= N+1 && idx_worst_contribution_sol == 0
        if P[idx].hv_contrib == -1 # we found the marked solution
            idx_worst_contribution_sol = idx
        end
        idx += 1
    end

    @assert idx_worst_contribution_sol != 0 "The marked solution isn't found"

    #------------- Remove the worst contribution solution from P -----------#

    # remove the marked solution (contributing less to the hypervolume)
    deleteat!(P, idx_worst_contribution_sol)

    debug && DEBUG_feasible_solutions(P)
    debug && DEBUG_correct_evaluations(P, hv = true)

    return P
end

function SMS_EMOA_solve(prob::_MOMKP, config::Config; seeding_sols = Vector{Sol}(undef, 0), max_t = 10.0, verbose = true)

    st = time()

    # create population
    P = INIT_random_feasible_population(config.size_pop_init, prob, seeding_sols = seeding_sols)

    # evaluate population
    F = INDICATOR_fast_non_dominated_sort!(P)
    for front in F
        INDICATOR_hypervolume_contribution!(front)
    end

    # iterate
    it = 0
    while !UTIL_termination_criteria_met(time()-st, max_t)
        it += 1
        P = SMS_EMOA_update(P, config)
    end

    # evaluate population
    F = INDICATOR_fast_non_dominated_sort!(P)

    last_front = UTIL_remove_doublons_from_front(F[1])

    sort!(last_front, by = sol -> sol.z[1])

    verbose && println("total time spent : ",time() - st)
    verbose && println("total nb of iter : ",it)
    verbose && println("total nb of sols : ",length(last_front),"\n")

    return last_front
end


text:
using Hiccup

view(x) =
  Dict(:type    => :html,
       :content => stringmime(MIME"text/html"(), x))

render(e::Editor, ::Void) =
  render(e, Atom.icon("check"))

render(::Console, ::Void) = nothing

render(::Inline, x::AbstractFloat) =
  isnan(x) || isinf(x) ?
    view(span(".constant.number", string(x))) :
    Dict(:type => :number, :value => Float64(x), :full => string(x))

@render Inline x::Expr begin
  text = string(x)
  length(split(text, "\n")) == 1 ?
    Model(Dict(:type => :code, :text => text)) :
    Tree(Text("Code"),
         [Model(Dict(:type => :code, :text => text))])
end

render(::Console, x::Expr) =
  Atom.msg("result", Dict(:type => :code, :text => string(x)))

@render Inline x::Text begin
  ls = split(chomp(string(x)), "\n")
  length(ls) > 1 ?
    Tree(Model(ls[1]), c(Model(join(ls[2:end], "\n")))) :
    span(ls[1])
end

getfield′(x, f) = isdefined(x, f) ? getfield(x, f) : fade("#undef#")

@render Inline x begin
  fields = fieldnames(typeof(x))
  if isempty(fields)
    span(c(render(Inline(), typeof(x)), "()"))
  else
    LazyTree(typeof(x), () -> [SubTree(Text("$f → "), getfield′(x, f)) for f in fields])
  end
end

@render Inline x::Type span(".support.type", string(x))

@render Inline x::Module span(".keyword.other", string(x))

@render Inline x::Symbol span(".constant.other.symbol", ":$x")

@render Inline x::Char span(".string.quoted.single", "'$x'")

@render Inline _::Void span(".constant", "nothing")

import Base.Docs: doc

isanon(f) = contains(string(f), "#")

@render Inline f::Function begin
  isanon(f) ? span(".support.function", "λ") :
    Tree(span(".support.function", string(typeof(f).name.mt.name)),
         [(Atom.CodeTools.hasdoc(f) ? [doc(f)] : [])..., methods(f)])
end

@render i::Inline xs::Vector begin
  length(xs) <= 25 ? children = handleundefs(xs) :
                     children = [handleundefs(xs, 1:10); span("..."); handleundefs(xs, length(xs)-9:length(xs))]
    Tree(span(c(render(i, eltype(xs)), Atom.fade("[$(length(xs))]"))),
         children)
end

@render Inline xs::Array begin
  Text(sprint(io -> show(IOContext(io, limit=true), MIME"text/plain"(), xs)))
end

@render i::Inline d::Dict begin
  j = 0
  st = Array{Atom.SubTree}(0)
  for (key, val) in d
    push!(st, SubTree(span(c(render(i, key), " → ")), val))
    j += 1
    j > 25 && (push!(st, SubTree(span("... → "), span("..."))); break)
  end
  Tree(span(c(strong("Dict"),
            Atom.fade(" $(eltype(d).parameters[1]) → $(eltype(d).parameters[2]) with $(length(d)) entries"))), st)
end

@render Inline x::Number span(".constant.number", sprint(show, x))

@render i::Inline x::Complex begin
  re, ima = reim(x)
  span(c(render(i, re), signbit(ima) ? " - " : " + ", render(i, abs(ima)), "im"))
end

@render Inline p::Ptr begin
  Row(Atom.fade(string(typeof(p))), Text(" @"),
       span(".constant.number", c("0x$(hex(UInt(p), Sys.WORD_SIZE>>2))")))
end

@render i::Inline x::AbstractString begin
  span(".string", c(render(i, Text(stringmime("text/plain", x)))))
end

render{sym}(i::Inline, x::Irrational{sym}) =
  render(i, span(c(string(sym), " = ", render(i, float(x)), "...")))

handleundefs(X::Vector) = handleundefs(X, 1:length(X))

function handleundefs(X::Vector, inds)
  Xout = Vector{Union{String, eltype(X)}}(length(inds))
  j = 1
  for i in inds
    Xout[j] = isdefined(X, i) ? X[i] : "#undef"
    j += 1
  end
  Xout
end

@render i::Inline xs::Tuple begin
  span(c("(", interpose(map(x->render(i, x), xs), ", ")..., ")"))
end

@render i::Inline md::Base.Markdown.MD begin
  mds = Atom.CodeTools.flatten(md)
  length(mds) == 1 ? Text(chomp(sprint(show, MIME"text/markdown"(), md))) :
                     Tree(Text("MD"), [HTML(sprint(show, MIME"text/html"(), md))])
end

include("methods.jl")


text:
module MechanicalSketch
import Luxor
import Luxor:
# Layout
Drawing, SVGimage,
paper_sizes,
Tiler, Partition,
rescale, currentdrawing,

finish, preview, snapshot,
origin, rulers, background,

@png, @pdf, @svg, @eps, @draw, @drawsvg, @savesvg,

# Turtle
Turtle, Pencolor, Penwidth, Forward, Turn, HueShift,

# Play

newpath, closepath, newsubpath,

BezierPath, BezierPathSegment, bezier, bezier′, bezier′′, makebezierpath, drawbezierpath, bezierpathtopoly, beziertopoly, pathtobezierpaths,
bezierfrompoints, beziercurvature, bezierstroke, setbezierhandles, beziersegmentangles,
shiftbezierhandles, brush,

strokepath, fillpath,

rect, box, cropmarks,

setantialias, setline, setlinecap, setlinejoin, setdash,

move, rmove, line, rule, rline, arrow, arrowhead, dimension, tickline,

BoundingBox, boxwidth, boxheight, boxdiagonal, boxaspectratio,
boxtop, boxbottom, boxtopleft, boxtopcenter, boxtopright, boxmiddleleft,
boxmiddlecenter, boxmiddleright, boxbottomleft, boxbottomcenter,
boxbottomright,

intersectboundingboxes, boundingboxesintersect, pointcrossesboundingbox,

BoxmapTile, boxmap,

circle, circlepath, crescent, ellipse, hypotrochoid, epitrochoid,
squircle, center3pts, curve,
arc, carc, arc2r, carc2r, isarcclockwise, arc2sagitta, carc2sagitta,
spiral, sector, intersection2circles,
intersection_line_circle, intersectionlinecircle, intersectioncirclecircle, ispointonline,
ispointonpoly,
intersectlinepoly, polyintersect, polyintersections, circlepointtangent,
circletangent2circles, pointinverse, pointcircletangent, circlecircleoutertangents,
circlecircleinnertangents, ellipseinquad,

ngon, ngonside, star, pie, polycross,
do_action, paint, paint_with_alpha, fillstroke, setstrokescale,
setblendextend,

Point, O, randompoint, randompointarray, midpoint, between, slope, intersectionlines,
pointlinedistance, getnearestpointonline, isinside,
perpendicular, crossproduct, dotproduct, distance,
prettypoly, polysmooth, polysplit, poly, simplify,  polycentroid,
polysortbyangle, polysortbydistance, offsetpoly, polyfit,
currentpoint, hascurrentpoint, getworldposition, anglethreepoints,

polyperimeter, polydistances, polyportion, polyremainder, nearestindex,
polyarea, polysample, insertvertices!,

polymove!, polyscale!, polyrotate!, polyreflect!,

@polar, polar,

strokepreserve, fillpreserve,
gsave, grestore, @layer,
scale, rotate, translate, rotate_point_around_point,
clip, clippreserve, clipreset,

getpath, getpathflat, pathtopoly,

fontface, fontsize, text, textpath, label,
textextents, textoutlines, textcurve, textcentred, textcentered, textright,
textcurvecentred, textcurvecentered,
textwrap, textlines, splittext, textbox, texttrack,

setcolor, setopacity, sethue, setgrey, setgray,
randomhue, randomcolor, @setcolor_str,
getmatrix, setmatrix, transform,

setfont, settext,

Blend, setblend, blend, addstop, blendadjust,
blendmatrix, rotationmatrix, scalingmatrix, translationmatrix,
cairotojuliamatrix, juliatocairomatrix, getrotation, getscale, gettranslation,

setmode, getmode,

GridHex, GridRect, nextgridpoint,

Table, highlightcells,

readpng, readsvg, placeimage, svgstring,

julialogo, juliacircles,

barchart,

mesh, setmesh, add_mesh_patch, mask,

# animation
Movie, Scene, animate,

lineartween, easeinquad, easeoutquad, easeinoutquad, easeincubic, easeoutcubic,
easeinoutcubic, easeinquart, easeoutquart, easeinoutquart, easeinquint, easeoutquint,
easeinoutquint, easeinsine, easeoutsine, easeinoutsine, easeinexpo, easeoutexpo,
easeinoutexpo, easeincirc, easeoutcirc, easeinoutcirc, easingflat, easeinoutinversequad, easeinoutbezier,

# noise
noise, initnoise,

# experimental polygon functions
polyremovecollinearpoints, polytriangulate!, polytriangulate,
ispointinsidetriangle, ispolyclockwise, polyorientation, ispolyconvex,

# triangles
trianglecircumcenter, triangleincenter, trianglecenter, triangleorthocenter,

# Turtle
Circle, Rectangle, Penup, Pendown, Message, Reposition, Orientation,

# misc
image_as_matrix, @imagematrix, image_as_matrix!, @imagematrix!,

# experimental
layoutgraph, Style, applystyle,
tidysvg,

# internals and etc (not exported by Luxor).

get_current_redvalue, get_current_greenvalue, get_current_bluevalue, get_current_alpha,
weighted_color_mean, Colors, parse, comp1, comp2, comp3, comp4

using MechanicalUnits
import MechanicalUnits.Unitfu: Power, oneunit, numtype, AbstractQuantity
@import_expand kW # Don't use ~ since that would also import WI, which we use elsewhere.
import Latexify, LaTeXStrings, NodeJS
using Latexify: @latexify, latexify, @latexrecipe
using LaTeXStrings: LaTeXString, @L_str
"""
Rotations given with units occur around the positive z axis, when y is up and x is to the right.
We don't dispatch on the dimension of angles, because the dimension is NonDims (for good reason).
In the context of this package, defining Angle is considered harmless:
"""
const Angle = Union{typeof(1.0°), typeof(1°), typeof(1.0f0°), typeof(1.0rad), typeof(1rad), typeof(1.0f0*rad)}
# TODO move Angle to numeric types...
import ColorSchemes
using ColorSchemes: getinverse, get
using ColorSchemes: RGB, RGBA
using ColorSchemes: HSV, HSVA, HSL, HSLA, LCHab, LCHabA, LCHuv, LCHuvA # Cylindrical hue colorspaces
using ColorSchemes: isoluminant_cgo_70_c39_n256, leonardo, PuOr_8, Greys_9 # We also define ColSchemeNoMiddle in 'colors.jl'
import Base: -, +, *, /, hypot, product, show
using Colors: @colorant_str
using FileIO: @format_str, File, save
#import ForwardDiff
import StaticArrays
using StaticArrays: SA, SVector
#import Interpolations:   interpolate, Linear, Flat, extrapolate, Extrapolation, Gridded
const CHORD_ZERO = 0.25
const FOIL_CHORD_POS = [0, 0.0125, 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1]
const FOIL_HALF_T = [0, 0.117, 0.1575, 0.2179, 0.2649, 0.3042, 0.4146, 0.4764, 0.5, 0.4816, 0.4149, 0.3159, 0.1989, 0.0811, 0.0306, 0]
const FOIL_CAMBER = [0, 0.0494, 0.0975, 0.19, 0.2775, 0.36, 0.64, 0.84, 0.96, 1, 0.96, 0.84, 0.64, 0.36, 0.19, 0]
const CHORD_LENGTH_DEFAULT = 1226mm
const THICKNESS_DEFAULT = 30.2mm
const CAMBER_DEFAULT = 18.15mm
const PALETTE = ColorSchemes.seaborn_bright

"""
Font size FS for Luxor, corresponding to Word points.
    fsize_pt = 11
    fsize_pt * ( 300 / 72 ) * 81 / (81 + 8)
"""
const FS = 9 * ( 300 / 72 )
"""
The unit EM, as in .css, corresponds to text + margins above and below
"""
const EM = Int(round(FS * 1.16))
"""
Text position is given from the baseline, but text extends below that.
"""
const TXTOFFV = Int(round(FS * 0.16))
"""
For line thickness, we are used to points.
"""
const PT = Int(round(EM / 12))
"""
WI is the pixel width of the figure. See HE and origo, O.

We want to make quality figures for A4 with 5 cm total margin. Width and height are
    a4_w_72dpi, _ = Luxor.paper_sizes["A4"]
    marginfrac = 50 / 210
    w_300f = a4_w_72dpi * (1 - marginfrac) * 300 / 72
    WI = Int(round(w_300f))
    HE = Int(round(w_300f * 2 / 3))

"""
global WI = 1889
"HE is the pixel height of the figure. See WI and orgo, O"
global HE = 1259
const SCALE = Dict{Symbol, Quantity{Float64}}()

Point(x::Quantity{T, D}, y::Quantity{S, D}) where {T, S, D} = Point(scale_to_pt(x), -scale_to_pt(y))

include("numeric_types.jl")

# Vector subtraction and division using QuantityTuples which don't need to have the exact same type.
+(tup1::QuantityTuple, tup2::QuantityTuple) = .+(tup1, tup2)
-(tup1::QuantityTuple, tup2::QuantityTuple) = .-(tup1, tup2)

# QuantityTuples and scalar operations

*(x::Number, tup::QuantityTuple) = .*(x, tup)
*(tup::QuantityTuple, x::Number) = .*(tup, x)
/(x::Number, tup::QuantityTuple) = ./(x, tup)
/(tup::QuantityTuple, x::Number) = ./(tup, x)
-(tup::QuantityTuple) = .-(tup)

# extending points with quantity tuples - will work for e.g. 3.0m/mm = 3000.0
+(p1::Point, shift::QuantityTuple) = p1 + Point(shift[1], shift[2])
-(p1::Point, shift::QuantityTuple) = p1 - Point(shift[1], shift[2])
*(p1::Point, shift::QuantityTuple) = p1 * Point(shift[1], shift[2])
/(p1::Point, shift::QuantityTuple) = p1 / Point(shift[1], shift[2])

# extending points with complex quantities
+(p1::Point, shift::ComplexQuantity) = p1 + QuantityTuple(shift)

# extend function in base with Point
hypot(p::Point) = hypot(p.x, p.y)
hypot(p::QuantityTuple) = hypot(p[1], p[2])

"""
    empty_figure(filename = "HiThere.png";
    backgroundcolor = color_with_lumin(PALETTE[8], 10),
    hue = PALETTE[8], height = missing, width = missing)

Establish a drawing sized for A4 300 dpi figures (WI, HE),
black on white figure, line width 3 pt  default.

Initial scale is set with set_scale_sketch()
    pt height = 20m
    pt height = 70m/s
    pt height = 20kN
"""
function empty_figure(;filename = :rec,
        backgroundcolor = color_with_lumin(PALETTE[8], 10),
        hue = PALETTE[8], height = missing, width = missing)
    fig = Drawing(WI, HE, filename)
    configure_mechanical(;backgroundcolor, hue, height, width)
    fig
end
"""
configure_mechanical(;
        backgroundcolor = color_with_lumin(PALETTE[8], 10),
        hue = PALETTE[8], height = missing, width = missing)

Establish a drawing sized for A4 300 dpi figures (WI, HE),
black on white figure, line width 3 pt  default.

Initial scale is set with set_scale_sketch()
pt height = 20m
pt height = 70m/s
pt height = 20kN
"""
function configure_mechanical(;
    backgroundcolor = color_with_lumin(PALETTE[8], 10),
    hue = PALETTE[8], height = missing, width = missing)
    # Font for the 'toy' text interface
    # (use e.g. JuliaMono for unicode symbols like ∈. There are no nice fonts for text AND math.
    fontface("Calibri")
    # 1 pt font size = 12/72 inch - by the book.
    # Letter spacing works differently here than in Word, so we adjust a little.
    fontsize(FS)
    setfont("Calibri", FS)
    background(backgroundcolor)
    sethue(hue)
    setline(0.5PT)
    setdash("solid")
    # Origo at centre
    origin()
    # Scale and rotation in pixels - x right, y up ('z' is in.). And rotations
    # are clockwise. Just deal with it. Or stick to using quantities. Functions taking
    # e.g. angles convert to ccw rotations.
    setmatrix([1, 0, 0, 1, WI / 2, HE / 2])
    @assert !ismissing(width) + !ismissing(height) < 2 "Width or height can be specified here."
    if ismissing(width) && ismissing(height)
        set_scale_sketch()
    elseif ismissing(width)
        set_scale_sketch(height, HE)
    else
        set_scale_sketch(width, WI)
    end
end

"""
    text(t, pt::Point, angle::T; kwargs) where {T <: Angle}
For angles with unit, use rotation around z axis.
"""
text(t, pt::Point, angle::T) where {T <: Angle} = text(t, pt; angle = - ustrip( angle |> rad))

polyrotate!(f, ang::Angle) = polyrotate!(f, - ustrip( ang |> rad))

include("lenient_num_analogue.jl")
include("scale.jl")
set_scale_sketch()
include("colors.jl")
include("dimension_aligned.jl")
include("foil.jl")
include("arrow.jl")
include("cart.jl")
include("rope.jl")
include("power.jl")
include("table.jl")
include("curves.jl")
include("labels.jl")
include("flow.jl")
include("runge_kutta.jl")
include("autodiff_unitfu.jl")
include("matrix_interpolation.jl")
include("streamline_convolution.jl")
include("matrix_drawing.jl")
include("colorlegends.jl")
include("colorlegends_vector.jl")
include("place_image.jl")
include("latex.jl")
include("chart.jl")
include("macros.jl")
end # module