text:
mutable struct LinearityExpr
    c
    linearity::Symbol
    coeff::Float64
end

LinearityExpr(c, linearity) = LinearityExpr(c, linearity, 1)

eval(c::LinearityExpr) = eval(get_expr(c))
Base.print(io::IO, c::LinearityExpr) = print(io::IO, "($(c.c),$(c.linearity))")
Base.show(io::IO, c::LinearityExpr) = print(io::IO, "($(c.c),$(c.linearity))")

LinearityExpr(c) = throw("Couldn't determine linearity of expression:\n$c")
LinearityExpr(c::Real) = LinearityExpr(c, :const)
LinearityExpr(c::LinearityExpr) = c  # Already processed; ignore!
function LinearityExpr(c::Expr)
    if c.head == :call
        for i = 2:length(c.args)
            c.args[i] = LinearityExpr(c.args[i])
        end

        args = c.args[2:end]
        if c.args[1] in [:+, :-]
            if check_for_linearity(:nonlinear, args)
                linearity = :nonlinear
            elseif check_for_linearity(:linear, args)
                linearity = :linear
            else
                linearity = :const
            end
        elseif c.args[1] == :*
            if check_for_linearity(:nonlinear, args) ||
               length(filter((c) -> check_linearity(:linear, c), args)) > 1
                linearity = :nonlinear
            elseif check_for_linearity(:linear, args)
                linearity = :linear
            else
                linearity = :const
            end
        elseif c.args[1] == :/
            if c.args[3].linearity != :const
                linearity = :nonlinear
            else
                linearity = c.args[2].linearity
            end
        elseif c.args[1] == :^
            if c.args[3].linearity != :const
                linearity = :nonlinear
            elseif c.args[2].linearity == :linear
                if c.args[3].c == 1
                    linearity = :linear
                else
                    linearity = :nonlinear
                end
            else
                linearity = c.args[2].linearity
            end
        elseif c.args[1] == :ifelse
            if c.args[2].linearity == :const
                # We know which branch to take already - simplify
                c.args[2] = pull_up_constants(c.args[2])
                if c.args[2].c != 0
                    linearity = c.args[3].linearity
                    c = c.args[3].c
                else
                    linearity = c.args[4].linearity
                    c = c.args[4].c
                end
            else
                linearity = :nonlinear
            end
        else
            if check_for_linearity(:linear, args) ||
               check_for_linearity(:nonlinear, args)
                linearity = :nonlinear
            else
                linearity = :const
            end
        end

    elseif c.head == :ref
        linearity = :linear
    elseif c.head == :comparison
        indices = 1:2:length(c.args)
        for i in indices
            c.args[i] = LinearityExpr(c.args[i])
        end
        if check_for_linearity(:linear, c.args[indices]) ||
           check_for_linearity(:nonlinear, c.args[indices])
            linearity = :nonlinear
        else
            linearity = :const
        end
    elseif c.head in [:&&, :||]
        for i = 1:length(c.args)
            c.args[i] = LinearityExpr(c.args[i])
        end
        if check_for_linearity(:linear, c.args) ||
           check_for_linearity(:nonlinear, c.args)
            linearity = :nonlinear
        else
            linearity = :const
        end
    end
    return LinearityExpr(c, linearity)
end

check_linearity(linearity::Symbol, c::LinearityExpr) = c.linearity == linearity

function check_for_linearity(linearity::Symbol, args::Array)
    return !isempty(filter((c) -> check_linearity(linearity, c), args))
end

get_expr(c) = c
get_expr(c::LinearityExpr) = get_expr(c.c)
function get_expr(c::Expr)
    map!(get_expr, c.args, c.args)
    return c
end

pull_up_constants(c) = c
function pull_up_constants(c::LinearityExpr)
    if c.linearity == :const
        c.c = eval(c)
    elseif isa(c.c, Expr)
        map!(pull_up_constants, c.c.args, c.c.args)
    end
    return c
end

function prune_linear_terms!(c::LinearityExpr, lin_constr::Dict{Int, Float64},
                             constant::Float64=0.0, negative_tree::Bool=false)
    if c.linearity != :nonlinear
        constant = add_linear_tree!(c, lin_constr, constant, negative_tree)
        c = LinearityExpr(:(0), :const)
        return true, c, constant
    else
        expr = c.c
        if expr.head == :call
            if expr.args[1] == :+
                n = length(expr.args)
                pruned = falses(n - 1)
                for i in 2:n
                    pruned[i - 1], expr.args[i], constant = prune_linear_terms!(
                        expr.args[i], lin_constr, constant, negative_tree)
                end
                if sum(.!pruned) > 1
                    inds = vcat([1], collect(2:n)[.!pruned])
                    c.c.args = expr.args[inds]
                else
                    c = expr.args[findfirst(.!pruned) + 1]
                end
            elseif expr.args[1] == :-
                if length(expr.args) == 3
                    pruned_first, expr.args[2], constant = prune_linear_terms!(
                        expr.args[2], lin_constr, constant, negative_tree)
                    pruned_second, expr.args[3], constant = prune_linear_terms!(
                        expr.args[3], lin_constr, constant, !negative_tree)
                    if pruned_first
                        new_expr = Expr(:call, :-, expr.args[3])
                        c = LinearityExpr(new_expr, :nonlinear)
                    elseif pruned_second
                        c = expr.args[2]
                    end
                end
            end
        end
        return false, c, constant
    end
end

function add_linear_tree!(c::LinearityExpr, lin_constr::Dict{Int, Float64},
                          constant::Float64=0.0, negative_tree::Bool=false)
    c = collate_linear_terms(c)
    negative_tree && negate(c)
    constant = add_tree_to_constr!(c, lin_constr, constant)
    return constant
end

function collate_linear_terms(c::LinearityExpr)
    if isa(c.c, Expr) && c.c.head == :call
        for i in 2:length(c.c.args)
            c.c.args[i] = collate_linear_terms(c.c.args[i])
        end
    end
    if c.linearity == :linear && isa(c.c, Expr) && c.c.head == :call
        func = c.c.args[1]
        if func == :-
            if length(c.c.args) == 2
                c = negate(c.c.args[2])
            else
                c.c.args[3] = negate(c.c.args[3])
                c.c.args[1] = :+
            end
        elseif func == :*
            # This can be n-ary multiplication, but there is only one variable
            # First, find index of variable in the args
            x = findfirst(ex -> ex.linearity == :linear, c.c.args[2:end]) + 1
            # Multiply variable term by each of the constants in the args
            x_expr = c.c.args[x]
            for i in setdiff(2:length(c.c.args), x)
                x_expr = multiply(x_expr, c.c.args[i].c)
            end
            c = x_expr
        elseif func == :/
            c = multiply(c.c.args[2], 1 / c.c.args[3].c)
        elseif func == :^
            c = c.c.args[2]
        end
    end
    return c
end

negate(c::LinearityExpr) = multiply(c::LinearityExpr, -1)
function multiply(c::LinearityExpr, a::Real)
    if isa(c.c, Expr) && c.c.head == :call
        @assert c.c.args[1] == :+
        map!(arg -> multiply(arg, a), c.c.args[2:end], c.c.args[2:end])
    elseif c.linearity == :const
        c.c *= a
    else
        c.coeff *= a
    end
    return c
end

function add_tree_to_constr!(c::LinearityExpr, lin_constr::Dict{Int, Float64},
                             constant::Float64)
    if isa(c.c, Expr) && c.c.head == :call
        @assert c.c.args[1] == :+
        for i in 2:length(c.c.args)
            constant = add_tree_to_constr!(c.c.args[i], lin_constr, constant)
        end
    elseif c.linearity == :const
        constant += c.c
    else
        lin_constr[c.c.args[2]] += c.coeff
    end
    return constant
end


text:
using RankingDifferences
using Test

@testset "RankingDifferences" begin
    include("ranking.jl")
end

ext: jl


text:
module TestIM

using IterationManagers
using Base.Test

include("newton.jl")

end  # module


text:
module TorchMLP
using PyCallChainRules.Torch: TorchModuleWrapper, torch, functorch, ispysetup

using Statistics
using BenchmarkTools

using Zygote

using CUDA

if !ispysetup[]
    return
end

include("pytorch_utils.jl")

batchsizes = [1, 8, 16, 32]
indim = 8
outdim = 4
hiddendim = 64

suite = BenchmarkGroup()

for device in ["cpu"]
    ss = suite["$device"] = BenchmarkGroup()
    for batchsize in batchsizes
        mlp = torch.nn.Sequential(torch.nn.Linear(indim, hiddendim), torch.nn.ReLU(), torch.nn.Linear(hiddendim, outdim))

        s = ss["bs=$batchsize"] = BenchmarkGroup()
        # Forward pass
        s["forward"] = BenchmarkGroup()
        forward_pass!(s["forward"], mlp, randn(Float32, indim), batchsize, device)

        # Gradient pass
        s["backward"] = BenchmarkGroup()
        backward_pass!(s["backward"], mlp, randn(Float32, indim), batchsize, device)
    end
end

for device in ["cuda"]
    ss = suite["$device"] = BenchmarkGroup()
    for batchsize in batchsizes
        mlp = torch.nn.Sequential(torch.nn.Linear(indim, hiddendim), torch.nn.ReLU(), torch.nn.Linear(hiddendim, outdim))

        s = ss["bs=$batchsize"] = BenchmarkGroup()
        # Forward pass
        s["forward"] = BenchmarkGroup()
        forward_pass!(s["forward"], mlp, randn(Float32, indim), batchsize, device)

        # Gradient pass
        s["backward"] = BenchmarkGroup()
        backward_pass!(s["backward"], mlp, randn(Float32, indim), batchsize, device)
    end
end

end
TorchMLP.suite


text: include("neldermead.jl")


text:
#
# Computes the density function of a list of values
#

function density(v;nbins=nothing,step=nothing,steptype="absolute",vmin=nothing,vmax=nothing)

  ndata = length(v)

  if vmin == nothing
    vmin = minimum(v)
  end
  if vmax == nothing
    vmax = maximum(v)
  end
  if nbins == nothing
    nbins = 100
  end
  if step == nothing
    step = (vmax - vmin)/nbins
  else
    # By default, the step size is absolute
    if steptype == "relative"
      step = step*(vmax-vmin)/nbins
    elseif steptype != "absolute"
      error(" steptype must be \"relative\" or \"absolute\"")
    end
  end

  x = Vector{Float64}(undef,nbins)
  df = Vector{Float64}(undef,nbins)

  binstep = (vmax - vmin)/nbins
  for i in 1:nbins
    x[i] = vmin + (i-1)*binstep + binstep/2
    nv = 0
    for j in 1:ndata
      if ( v[j] > x[i] - step/2 ) && ( v[j] <= x[i] + step/2 )
        nv = nv + 1
      end
    end
    binsize = min(vmax,x[i]+step/2) - max(vmin,x[i]-step/2)
    df[i] = nv/(binsize*ndata)
  end

  return x, df

end


text:
include("wasm-vector.jl")

using .WASMVector

using ExportWebAssembly

function myfun(p)
    x = WVector{Float64}(p, 4)
    x[1] + x[2]
end

@code_warntype myfun(Int32(1))
irgen(myfun, Tuple{Int32}, optimize = true)
write_wasm("myfun.wasm", myfun, Tuple{Int32})
wasm2wast("myfun.wasm")

function myfun2(i)
    x = WVector{Int8}(WASMVector.S"hello", 6)
    x[i]
end

@code_warntype myfun2(Int32(1))
irgen(myfun2, Tuple{Int32}, optimize = true)
write_wasm("myfun2.wasm", myfun2, Tuple{Int32})
wasm2wast("myfun2.wasm")

function myfun3(i)
    x = WVector{Float64}(Int32(6))
    x[i] = 3.0
    x[1] + x[i]
end

@code_warntype myfun3(Int32(1))
irgen(myfun3, Tuple{Int32}, optimize = true)
write_wasm("myfun3.wasm", myfun3, Tuple{Int32}, linkwith = ("../memory/memory-utils.o", "../memory/emmalloc.o"))
#write_wasm("myfun3.wasm", myfun3, Tuple{Int32})


text:
export BrickletMotorizedLinearPotiPositionCallbackConfiguration
struct BrickletMotorizedLinearPotiPositionCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletMotorizedLinearPotiMotorPosition
struct BrickletMotorizedLinearPotiMotorPosition
    position::Integer
    drive_mode::Integer
    hold_position::Bool
    position_reached::Bool
end

export BrickletMotorizedLinearPotiSPITFPErrorCount
struct BrickletMotorizedLinearPotiSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

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

export BrickletMotorizedLinearPoti
"""
Motorized Linear Potentiometer
"""
mutable struct BrickletMotorizedLinearPoti <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_position
"""
    $(SIGNATURES)

Returns the position of the linear potentiometer. The value is
between 0 (slider down) and 100 (slider up).

If you want to get the value periodically, it is recommended to use the
:cb:`Position` callback. You can set the callback configuration
with :func:`Set Position Callback Configuration`.
"""
function get_position(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_position()
end

export set_position_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Position` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Position` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_position_callback_configuration(device::BrickletMotorizedLinearPoti, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_position_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_position_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Position Callback Configuration`.
"""
function get_position_callback_configuration(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_position_callback_configuration()
end

export set_motor_position
"""
    $(SIGNATURES)

Sets the position of the potentiometer. The motorized potentiometer will
immediately start to approach the position. Depending on the chosen drive mode,
the position will either be reached as fast as possible or in a slow but smooth
motion.

The position has to be between 0 (slider down) and 100 (slider up).

If you set the hold position parameter to true, the position will automatically
be retained. If a user changes the position of the potentiometer, it will
automatically drive back to the original set point.

If the hold position parameter is set to false, the potentiometer can be changed
again by the user as soon as the set point was reached once.
"""
function set_motor_position(device::BrickletMotorizedLinearPoti, position, drive_mode, hold_position)
    device.deviceInternal.set_motor_position(position, drive_mode, hold_position)
end

export get_motor_position
"""
    $(SIGNATURES)

Returns the last motor position as set by :func:`Set Motor Position`. This is not
the current position (use :func:`Get Position` to get the current position). This
is the last used set point and configuration.

The position reached parameter is true if the position has been reached at one point.
The position may have been changed again in the meantime by the user.
"""
function get_motor_position(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_motor_position()
end

export calibrate
"""
    $(SIGNATURES)

Starts a calibration procedure. The potentiometer will be driven to the extreme
points to calibrate the potentiometer.

The calibration is saved in flash, it does not have to be called on every start up.

The Motorized Linear Poti Bricklet is already factory-calibrated during
testing at Tinkerforge.
"""
function calibrate(device::BrickletMotorizedLinearPoti)
    device.deviceInternal.calibrate()
end

export set_position_reached_callback_configuration
"""
    $(SIGNATURES)

Enables/Disables :cb:`Position Reached` callback.
"""
function set_position_reached_callback_configuration(device::BrickletMotorizedLinearPoti, enabled)
    device.deviceInternal.set_position_reached_callback_configuration(enabled)
end

export get_position_reached_callback_configuration
"""
    $(SIGNATURES)

Returns the :cb:`Position Reached` callback configuration
as set by :func:`Set Position Reached Callback Configuration`.
"""
function get_position_reached_callback_configuration(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_position_reached_callback_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletMotorizedLinearPoti, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletMotorizedLinearPoti, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletMotorizedLinearPoti, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletMotorizedLinearPoti, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletMotorizedLinearPoti)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletMotorizedLinearPoti, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletMotorizedLinearPoti)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

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

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

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

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


text:
using LinearAlgebra

function Cn(v, n::Integer)
    θ = 2*π / n
    return rotation_matrix(v, θ)
end

"""
    Molecules.rotation_matrix(v, θ::AbstractFloat)

Returns a rotation matrix by θ about a rotation axis v
"""
function rotation_matrix(V, θ::AbstractFloat)
    cosθ = cos(θ)
    sinθ = sin(θ)
    v = normalize(V)
    a = [1,2,3]
    O = zeros(eltype(v), (3,3))
    O .+= 1 - cosθ
    for i = 1:3, j = 1:3
        if i == j
            O[i,i] *= v[i]^2
            O[i,i] += cosθ
        else
            O[i,j] *= v[i]*v[j]
            b = [i,j]
            C = [i for i in a if i ∉ b][1]
            if i < j
                O[i,j] += (-1)^(i+j) * v[C]*sinθ
            else
                O[i,j] += (-1)^(i+j-1) * v[C]*sinθ
            end
        end
    end
    return O
end

function rotation_matrixd(v, θ::AbstractFloat)
    r = deg2rad(θ)
    return rotation_matrix(v, r)
end

"""
    Molecules.reflection_matrix(v::Vector)

Returns a reflection matrix through the plane with normal vector v
"""
function reflection_matrix(V)
    O = zeros(eltype(V), (3,3))
    v = normalize(V)
    for i = 1:3, j = i:3
        if i == j
            O[i,i] = 1 - 2*v[i]^2
        else
            O[i,j] = -2 * v[i] * v[j]
            O[j,i] = O[i,j]
        end
    end
    return O
end

function σ(v)
    return reflection_matrix(v)
end

"""
    Molecules.Sn(v::Vector, n::int)

Returns Sn improper rotation about vector v
"""
function Sn(v, n)
    return Cn(v, n) * σ(v)
end

"""
    Molecules.inversion_matrix()

Returns a diagonal matrix with -1 along the diagonal
"""
function inversion_matrix()
    a = zeros((3,3))
    for i = 1:3
        a[i,i] = -1
    end
    return a
end

function i()
    return inversion_matrix()
end


text:
# https://github.com/felixphew/pokepaste/blob/v3/syntax.go

const STATS_NAMES = ("HP", "Atk", "Def", "SpA", "SpD", "Spe")

function _parse_stats!(stats, text)
    m = match(
        r"^((\d+) HP)?( / )?((\d+) Atk)?( / )?((\d+) Def)?( / )?((\d+) SpA)?( / )?((\d+) SpD)?( / )?((\d+) Spe)?( *)$",
        text
    )
    for (i, field) in enumerate(STATS_FIELDS)
        t = i * 3 - 1
        if !isnothing(m.captures[t])
            setfield!(stats, field, parse(Int, m.captures[t]))
        end
    end
    stats
end

_stats_from_iv_ev_hp(iv, ev, base, level) = if base == 1
    1
else
    (2base + iv + ev ÷ 4) * level ÷ 100 + 10 + level
end
function _stats_from_iv_ev(
    iv::Int64, ev::Int64, base::Int64, level::Int64,
    dex_field = nothing, up = nothing, down = nothing
)
    dex_field == :hp && return _stats_from_iv_ev_hp(iv, ev, base, level)
    value = (2base + iv + ev ÷ 4) * level ÷ 100 + 5
    if dex_field ≡ up
        value = value * 11 ÷ 10
    elseif dex_field ≡ down
        value = value * 9 ÷ 10
    end
    value
end

function calculate_stats(poke::DexPokemon, level::Int64, evs::Stats, ivs::Stats, nature::Missable{Nature})
    stats = Stats()
    up, down = VsRecorder.get_nature_effect(nature)
    for (field, dex_field) in zip(STATS_FIELDS, DEX_STATS_FIELDS)
        base = getfield(poke.base_stats, dex_field)
        ev = getfield(evs, field)
        if ismissing(ev)
            ev = 0
        end
        iv = getfield(ivs, field)
        if ismissing(iv)
            iv = 31
        end
        value = _stats_from_iv_ev(iv, ev, base, level, dex_field, up, down)
        setfield!(stats, field, value)
    end
    stats
end

function import_poke(input::AbstractString)
    input = strip(input)
    lines = split(input, r"[\r\n]", keepempty = false)
    head = match(
        r"^(?:(.* \()([A-Z][a-z0-9:']+\.?(?:[- ][A-Za-z][a-z0-9:']*\.?)*)(\))|([A-Z][a-z0-9:']+\.?(?:[- ][A-Za-z][a-z0-9:']*\.?)*))(?:( \()([MF])(\)))?(?:( @ )([A-Z][a-z0-9:']*(?:[- ][A-Z][a-z0-9:']*)*))?( *)$",
        lines[1]
    )
    @assert !isnothing(head)
    poke = Pokemon(level = 100)  # PokePaste's default level is 100
    pdex = poke_dex()
    if !isnothing(head.captures[2])
        pid = head.captures[2]
        poke.id = search_dex(pdex, pid)
        poke.name = split(head.captures[1], '-')[1] |> strip
    elseif !isnothing(head.captures[4])
        pid = head.captures[4]
        poke.id = search_dex(pdex, pid)
    end
    dex_poke = pdex[poke.id]
    poke.ability = dex_poke.abilities[1]
    if !isnothing(head.captures[6])
        gender = head.captures[6]
        poke.gender = gender == "M" ? GENDER_MALE : gender == "F" ? GENDER_FEMALE : GENDER_NULL
    end
    if !isnothing(head.captures[9])
        iid = head.captures[9]
        poke.item = search_dex(item_list(VsI18n.EN), iid)
    end
    mdex = move_dex()
    adex = ability_list(VsI18n.EN)
    evs = Stats()
    ivs = Stats()
    for line in lines[2:end]
        m = match(
            r"^(-)( ([A-Z][a-z\']*(?:[- ][A-Za-z][a-z\']*)*)(?: \[([A-Z][a-z]+)\])?(?: / [A-Z][a-z\']*(?:[- ][A-Za-z][a-z\']*)*)* *)$",
            line
        )
        if !isnothing(m)
            mid = m.captures[3]
            push!(poke.moves, search_dex(mdex, mid))
            continue
        end
        m = match(
            r"^(.+) Nature.*$",
            line
        )
        if !isnothing(m)
            poke.nature = search_dex(nature_list(VsI18n.EN), strip(m.captures[1]))
            continue
        end
        tmp = split(line, ':', limit = 2)
        length(tmp) == 2 || continue
        key, values = tmp
        values = strip(values)
        if key == "EVs"
            _parse_stats!(evs, values)
        elseif key == "IVs"
            _parse_stats!(ivs, values)
        elseif key == "Ability"
            if length(dex_poke.abilities) == 1
                poke.ability = dex_poke.abilities[1]
                continue
            end
            for id in dex_poke.abilities
                if adex[id] == values
                    poke.ability = id
                    continue
                end
            end
        elseif key == "Level"
            poke.level = parse(Int, values)
        elseif key == "Shiny"
            poke.shiny = values == "Yes"
        else
        end
    end
    poke.stats = calculate_stats(dex_poke, poke.level, evs, ivs, poke.nature)
    poke
end

"""
    import_team(pokepaste_input) -> Team

Import a string formatted in PokePaste syntax as a team. Only support English.
"""
function import_team(
    input::AbstractString;
    title = "",
    author = missing,
    notes = ""
)
    team = Team(title = title, author = author, notes = notes)
    for poke in split(input, r"\r?\n\r?\n", keepempty = false)
        push!(team.pokemons, import_poke(poke))
    end
    team
end

function _stats_to_iv_ev_hp(value, base, level)
    if _stats_from_iv_ev_hp(31, 0, base, level) == value
        return 31, 0
    end
    t = cld((value - level - 10) * 100, level) - 2base
    if t ≤ 31
        t, 0
    else
        31, (t - 31) * 4
    end
end
function _stats_to_iv_ev(
    value::Int64, base::Int64, level::Int64,
    dex_field = nothing, up = nothing, down = nothing
)
    dex_field == :hp && return _stats_to_iv_ev_hp(value, base, level)
    if _stats_from_iv_ev(31, 0, base, level, dex_field, up, down) == value
        return 31, 0
    end
    if dex_field ≡ up
        value = cld(value * 10, 11)
    elseif dex_field ≡ down
        value = cld(value * 10, 9)
    end
    t = cld((value - 5) * 100, level) - 2base
    if t ≤ 31
        t, 0
    else
        31, (t - 31) * 4
    end
end

_get_total_evs(evs::Stats) = sum(getfield(evs, field) for field in STATS_FIELDS)

function _try_natures(poke::Pokemon, level::Int64, up = nothing, down = nothing)
    if isnothing(up) ≠ isnothing(down)
        first = isnothing(up) ? down : up
        for second in DEX_STATS_FIELDS
            first == second && continue
            nature = VsRecorder.get_nature(
                isnothing(up) ? second : first,
                isnothing(up) ? first : second
            )
            nature, ivs, evs = inv_calculate_stats(poke, level, nature)
            _get_total_evs(evs) ≤ 510 && return nature, ivs, evs
        end
    end
    for up in DEX_STATS_FIELDS
        for down in DEX_STATS_FIELDS
            up == down && continue
            nature = VsRecorder.get_nature(up, down)
            nature, ivs, evs = inv_calculate_stats(poke, level, nature)
            _get_total_evs(evs) ≤ 510 && return nature, ivs, evs
        end
    end
    @warn "No valid nature found for $(poke.id)"
    inv_calculate_stats(poke, level, VsRecorder.NATURE_NEUTURAL_DEFAULT)
end

function inv_calculate_stats(poke::Pokemon, level::Int64, nature = poke.nature)
    dex_poke = poke_dex()[poke.id]
    base_stats = dex_poke.base_stats
    stats = poke.stats
    ivs = Stats()
    evs = Stats()
    ismissing(stats) && return nature, ivs, evs
    if !ismissing(stats.hp)
        ivs.hp, evs.hp = _stats_to_iv_ev(stats.hp, base_stats.hp, level, :hp)
    end
    up, down = VsRecorder.get_nature_effect(nature)
    for (field, dex_field) in zip(STATS_FIELDS, DEX_STATS_FIELDS)
        v = getfield(stats, field)
        ismissing(v) && continue
        base = getfield(base_stats, dex_field)
        iv, ev = _stats_to_iv_ev(v, base, level, dex_field, up, down)
        if ismissing(nature)
            if iv < 0
                down = dex_field
                iv, ev = _stats_to_iv_ev(v, base, level, dex_field, up, down)
            elseif ev > 252
                up = dex_field
                iv, ev = _stats_to_iv_ev(v, base, level, dex_field, up, down)
            end
        end
        setfield!.((ivs, evs), (field,), (iv, ev))
    end
    if isnothing(up) ≠ isnothing(down) || ismissing(nature) && _get_total_evs(evs) > 510
        # Just brute force it
        return _try_natures(poke, level, up, down)
    end
    nature = VsRecorder.get_nature(up, down)
    nature, ivs, evs
end

function _print_stats(io::IO, prefix, stats::Stats, default)
    tio = IOBuffer()
    printed = false
    for (i, field) in enumerate(STATS_FIELDS)
        v = getfield(stats, field)
        (ismissing(v) || v == default) && continue
        print(
            tio, printed ? " / " : "", v, ' ',
            STATS_NAMES[i]
        )
        printed = true
    end
    if printed
        println(io, prefix, String(take!(tio)))
    end
end

function export_poke(poke::Pokemon; language = "en")
    if !(language isa VsRecorder.GameLanguage)
        language = VsRecorder.get_game_language(language)
    end
    lang = lowercase(string(language))
    io = IOBuffer()
    if !ismissing(poke.name)
        print(io, poke.name)
        print(io, " (", i18n(poke.id, language = lang), ')')
    else
        print(io, i18n(poke.id, language = lang))
    end
    if !ismissing(poke.gender) && poke.gender ≢ GENDER_NULL
        print(io, " (", poke.gender ≡ GENDER_MALE ? 'M' : 'F', ')')
    end
    if !ismissing(poke.item)
        print(io, " @ ", i18n(poke.item, language = lang))
    end
    println(io)
    if !ismissing(poke.shiny)
        println(io, "Shiny: ", poke.shiny ? "Yes" : "No")
    end
    if !ismissing(poke.ability)
        println(io, "Ability: ", i18n(poke.ability, language = lang))
    end
    level = ismissing(poke.level) ? 50 : poke.level
    println(io, "Level: ", level)
    nature, ivs, evs = inv_calculate_stats(poke, level)
    _print_stats(io, "EVs: ", evs, 0)
    if !ismissing(nature)
        println(io, i18n(nature, language = lang), " Nature")
    end
    _print_stats(io, "IVs: ", ivs, 31)
    for move in poke.moves
        println(io, "- ", i18n(move, language = lang))
    end
    String(take!(io))
end

"""
    export_team(team; language = "en") -> String

Export a team to a string formatted in PokePaste syntax.
"""
function export_team(team::Team; language = "en")
    if !(language isa VsRecorder.GameLanguage)
        language = VsRecorder.get_game_language(language)
    end
    join(
        (export_poke(poke; language = language) for poke in team.pokemons),
        '\n'
    )
end


text:
"""
    conformal_cubed_sphere_mapping(x, y)

Conformal mapping of a cube onto a sphere. Maps `(x, y)` on the north-pole face of a cube
to (X, Y, Z) coordinates in physical space. The face is oriented normal to Z-axis with
X and Y increasing with x and y.

The input coordinates must lie within the range -1 <= x <= 1,  -1 <= y <= 1.

This numerical conformal mapping is described by Rančić et al. (1996).

This is a Julia translation of MATLAB code from MITgcm [1] that is based on
Fortran 77 code from Jim Purser & Misha Rančić.

[1] http://wwwcvs.mitgcm.org/viewvc/MITgcm/MITgcm_contrib/high_res_cube/matlab-grid-generator/map_xy2xyz.m?view=markup

Rančić et al., (1996): Quarterly Journal of the Royal Meteorological Society, A global shallow-water model
    using an expanded spherical cube - Gnomonic versus conformal coordinates
"""
function conformal_cubed_sphere_mapping(x, y)
    X = xᶜ = abs(x)
    Y = yᶜ = abs(y)

    kxy = yᶜ > xᶜ

    xᶜ = 1 - xᶜ
    yᶜ = 1 - yᶜ

    kxy && (xᶜ = 1 - Y)
    kxy && (yᶜ = 1 - X)

    Z = ((xᶜ + im * yᶜ) / 2)^4
    W = W_Rancic(Z)

    im³ = im^(1/3)
    ra = √3 - 1
    cb = -1 + im
    cc = ra * cb / 2

    W = im³ * (W * im)^(1/3)
    W = (W - ra) / (cb + cc * W)
    X, Y = reim(W)

    H = 2 / (1 + X^2 + Y^2)
    X = X * H
    Y = Y * H
    Z = H - 1

    if kxy
        X, Y = Y, X
    end

    y < 0 && (Y = -Y)
    x < 0 && (X = -X)

    # Fix truncation for x = 0 or y = 0.
    x == 0 && (X = 0)
    y == 0 && (Y = 0)

    return X, Y, Z
end

"""
   conformal_cubed_sphere_inverse_mapping(X, Y, Z)

Inverse mapping for conformal cube sphere for quadrant of north-pole face in which X and Y are both
positive. All other mappings to other cube face coordinates can be recovered from rotations of this map.
There a 3 other quadrants for the north-pole face and five other faces for a total of twenty-four quadrants.
Because of symmetry only the reverse for a single quadrant is needed. Because of branch cuts and the complex
transform the inverse mappings are multi-valued in general, using a single quadrant case allows a simple
set of rules to be applied.

The mapping is valid for the cube face quadrant defined by 0 < x < 1 and 0 < y < 1, where a full cube
face has extent -1 < x < 1 and -1 < y < 1. The quadrant for the mapping is from a cube face that has
"north-pole" at its center (x=0,y=0) i.e has X,Y,Z = (0,0,1) at its center. The valid ranges of X and Y
for this mapping and convention are a quadrant defined be geodesics that connect the points A, B, C and D,
on the shell of a sphere of radius R with X,Y coordinates as follows

A = (0,0)
B = (√2,0)
C = (√3/3,√3/3)
D = (0,√2)

"""
function conformal_cubed_sphere_inverse_mapping(X, Y, Z)
    H  = Z + 1
    Xˢ = X/H
    Yˢ = Y/H
    ω  = Xˢ + im*Yˢ

    ra = √3 - 1
    cb = -1 + im
    cc = ra * cb / 2
    ω⁰ = (ω*cb + ra)/(1-ω*cc)
    W⁰ = im*ω⁰^3*im
    Z  = Z_Rancic(W⁰)
    z  = (Z^0.25)*2
    x, y = reim(z)

    kxy = abs(y) > abs(x)
    xx = x
    yy = y
    !kxy && ( x = 1 - abs(yy) )
    !kxy && ( y = 1 - abs(xx) )

    xf = x
    yf = y
    ( X < Y ) && ( xf = y  )
    ( X < Y ) && ( yf = x  )
    x = xf
    y = yf

    return x, y
end

W_Rancic(Z) = sum(A_Rancic[k] * Z^(k-1) for k in 1:length(A_Rancic))
Z_Rancic(W) = sum(B_Rancic[k] * W^(k-1) for k in 1:length(B_Rancic))


text:
ds = Systems.standardmap()
p_index = 1

ds, pmin, pmax = Systems.standardmap(), 0.6, 1.2

integ = integrator(ds)

for i in 1:2
od, xmin, xmax = InteractiveChaos.minimal_normalized_od(integ, i, p_index, pmin, pmax,
                 100, 100, 100, get_state(ds))

@test xmin ≥ 0
@test xmax ≤ 2π
end

od, xmin, xmax = InteractiveChaos.minimal_normalized_od(integ, 1, p_index, pmin, pmax,
                 100, 100, 100, get_state(ds))

# Now to simulate selecting a rectangle:
rect = Observable((origin = [0., 0], widths = [1., 1.]))
OD = Observable(od)
od0 = copy(od)

on(rect) do r
    pmin, xmin = r.origin
    pmax, xmax = r.origin + r.widths

    OD[] = InteractiveChaos.minimal_normalized_od(
        integ, 1,  p_index, pmin, pmax,
        100, 100, 100, get_state(ds), xmin, xmax
    )
end

# change rectangle:
rect[] = (origin = [0.8, 0.2], widths = [1., 5.0])

@test od0 != OD[]


text:
export generate_doc,org_string, create_documented_item_array, create_documented_item

#+OrgString
#
# Generates header for Org export
#
# *Parameters:*
# - *uuid_link:* add a link target <<uuid_link>> (used by index for instance)
# - *extra_link:* add an extra link target <<extra_link>> (used by index for instance)
# - *header_level:* [[(test)]]
#   - *-1*: no header (this function returns "")
#   - *0*: header -
#   - *n*: header *...* with n stars
#
function org_string_documented_item_header(di::Documented_Item;
                                           extra_link::String="", # (ref:test)
                                           uuid_link::String="",
                                           header_level::Int=0)::String
    s=""
    if header_level<0
        return s
    end

    if header_level==0
        header_str="- @@latex:\\phantomsection@@"
    else
        header_str=String(fill('*',header_level))
    end
    identifier_=identifier(di)
    s=s*"$(header_str) *=$(identifier(di))=* "

    if !isempty(uuid_link)
        s=s*" <<$(uuid_link)>>"
    end

    return s
end

#+OrgString
#
# Generates footer for Org export
#
function org_string_documented_item_footer(di::Documented_Item;
                                           index_link::String="")::String
    s=""
    generate_index_link = (index_link!="")

    if !isempty(filename(di))
        s=s*create_file_org_link(di)
        if generate_index_link
            s=s*", "
        end
    end

    if  generate_index_link
        s=s*"[[$(index_link)][back to index]]"
    end
    s=s*"\n"

    return s
end

#+OrgString
# Prints a =Documented_Item= structure
#
# This function is reserved for internal uses
#
# *Note:* we need all the =di_array::Array{Documented_Item,1}= and not simply =di::Documented_Item=
#         this global array is use to generate and check extra link (L:link)
#
# *Optional arguments:*
# - souce_link: uses item filename to create a link to source
# - uuid_link: if != "" add a link target <uuid_link> (for used to create index for instance)
function org_string_documented_item(di::Documented_Item,
                                    di_array::Array{Documented_Item,1},
                                    di_array_universe::Array{Documented_Item,1};
                                    link_prefix::String=randstring(),
                                    source_link::Bool=true,
                                    uuid_link::String="",
                                    index_link::String="",
                                    header_level::Int=0,
                                    boxingModule::String="BoxingModule",
                                    with_body::Bool=false)::String
    # check links for error
    check_for_link_error(di_array,di_array_universe)

    # Header
    s=""
    s=s*org_string_documented_item_header(di,
                                          uuid_link=uuid_link,
                                          header_level=header_level)
    # Header, need completion in case of extra link (L:something)
    if !isempty(link(di))
        if header_level<0
            s=s*" @@latex:\\phantomsection@@ "
        end
        s=s*" <<$(link_prefix*link(di))>>\n"
    else
        s=s*"\n"
    end

    # Code
    s=s*org_string_code(di)

    # Doc
    s=s*org_string_comment(di,
                           di_array,
                           di_array_universe,
                           link_prefix,
                           boxingModule)

     if with_body
        s=s*org_string_code_with_body(di)
     end

    # Footer
    s=s*org_string_documented_item_footer(di,
                                          index_link=index_link)

    return s
end

#+Internal
#
# Given a transformation t(x) creates groups (a vector of UnitRange) where
# t(x_i), i in group_k is constant
#
# *Caveat:* to make sense, the input array must be sorted.
function create_group(array::Array{T,1},t::Function)::Array{UnitRange{Int},1} where {T}
    groups=Array{UnitRange{Int},1}()
    n = length(array)
    if length(array)==0
        return groups
    end
    i,j=1,2
    current_identifier = t(array[1])
    while(j<=n)
        if current_identifier!=t(array[j])
            push!(groups,i:j-1)
            i=j
            current_identifier=t(array[j])
        end
        j=j+1
    end
    push!(groups,i:n)
    return groups
end

#+Internal
# Creates groups sharing the same identifier with correct lexical order
#
# This version is *case nonsensitive*
#
function group_by_identifier_case_nonsensitive(di_array::Array{Documented_Item,1})::Array{UnitRange{Int},1}

    # Group without taking into account uppercase/lowercase
    #
    t(x)=uppercase(identifier(x))
    sort!(di_array, by = t)
    groups = create_group(di_array,t)

    # Sort taking account uppercase/lowercase for each group
    #
    for group_i in groups
        sort!(view(di_array,group_i), by = x -> identifier(x))
    end
    # Rebuild group taking account uppercase/lowercase
    #
    groups = create_group(di_array,x->identifier(x))

    return groups
end

#+Internal
# Creates groups sharing the same identifier with correct lexical order
#
# This version is *case sensitive* (capitals are in front)
#
function group_by_identifier_case_sensitive(di_array::Array{Documented_Item,1})::Array{UnitRange{Int},1}

    # Group without taking into account uppercase/lowercase
    #
    t(x)=identifier(x)
    sort!(di_array, by = t)
    groups = create_group(di_array,t)

    return groups
end

#+OrgString
#
# *Arguments:* see [[print_org_doc_API][]]
#
function org_string_documented_item_array(di_array::Array{Documented_Item,1},
                                          predicate::Function;
                                          header_level::Int=0,
                                          link_prefix::String=randstring(),
                                          complete_link::Bool=false,
                                          case_sensitive::Bool=true,
                                          boxingModule::String="BoxingModule",
                       with_body::Bool=false)::String

    # extract what we need
    di_array_copy=filter(predicate,di_array)

    # try to complete links
    if complete_link
        links = extract_links(di_array_copy)
        di_array_copy_complement=filter(x->!predicate(x),di_array)
        @assert length(di_array_copy_complement) + length(di_array_copy) == length(di_array)
        while length(links)>0
            link_i=pop!(links)
            idx=get_item_idx_from_link_target(link_i,di_array_copy)
            if isempty(idx)
                # We have not found the link target,
                # try with complementary
                #
                idx=get_item_idx_from_link_target(link_i,di_array_copy_complement)

                if !isempty(idx)
                    # The target exists in complementary, we must:
                    # 1/ add all the new links
                    # 2/ move these new items from complement to our current item list
                    #    note:we can delete them as there is only _one_ link target per item

                    # 1/
                    links=vcat(links,extract_links(di_array_copy_complement[idx]))
                    links=remove_link_duplicate(links)
                    # 2/
                    di_array_copy=vcat(di_array_copy,di_array_copy_complement[idx])
                    di_array_copy_complement=deleteat!(di_array_copy_complement,idx)
                end
            end
        end
    end

    # lexical grouping
    n_di = length(di_array_copy)
    s=""
    if n_di==0
        return s
    end
    #    groups = group_by_identifier_case_nonsensitive(di_array_copy)
    # sort!(di_array_copy,by=x->identifier(x))
    # groups = create_group(di_array_copy,x->identifier(x))

    if case_sensitive
        groups = group_by_identifier_case_sensitive(di_array_copy)
    else
        groups = group_by_identifier_case_nonsensitive(di_array_copy)
    end

    n_groups = length(groups)

    # generates (hidden) uuid
    uuid = Array{String,1}()
    for i in 1:n_di
        push!(uuid,randstring())
    end

    # Creates index string
    if header_level>=0

        function get_identifier_first_letter(idx)
            first_letter = identifier(di_array_copy[idx])[1]

            if !case_sensitive
                first_letter = uppercase(first_letter)
            end
            return first_letter
        end

        # if case_sensitive
        #     get_identifier_first_letter(idx) = identifier(di_array_copy[idx])[1]
        # else
        #     get_identifier_first_letter(idx) = uppercase(identifier(di_array_copy[idx])[1])
        # end

        alpha_group = create_group(groups,r->get_identifier_first_letter(first(r)))

        index_label = randstring()
        s="<<$(index_label)>> *Index:* "

        for alpha in alpha_group
            s=s*"*[$(get_identifier_first_letter(first(groups[first(alpha)])))]* "
            for group_j in groups[alpha]
                first_item_idx = first(group_j)
                s=s*"[[$(uuid[first_item_idx])][$(identifier(di_array_copy[first_item_idx]))]]"
                if  group_j!=last(groups[alpha])
                    s=s*","
                end
                s=s*" "
            end
        end
        s=s*"\n"
    else
        index_label = ""
    end

    # Generate documented items
    for group_j in groups
        for j in group_j
            is_first = j==first(group_j)
            s=s*org_string_documented_item(di_array_copy[j],
                                           di_array_copy,
                                           di_array, # universe (mainly use to find external links)
                                           header_level = (is_first ? header_level : -1),
                                           index_link   = index_label,
                                           uuid_link    = uuid[j],
                                           link_prefix  = link_prefix,
                                           boxingModule=boxingModule,
                                           with_body=with_body)
        end
    end

    return s
end

#+OrgString
#
# Generates org doc string
#
# This function filters output according to these options:
#
# *tag:* filter according to this tag, "" <-> take all tags
#
function org_string_documented_item_array(di_array::Array{Documented_Item,1};
                                          tag::Union{String,Array{String,1}}="",
                                          tag_to_ignore::Union{String,Array{String,1}}="",
                                          identifier::String="",
                                          header_level::Int=0,
                                          link_prefix::String=randstring(),
                                          complete_link::Bool=false,
                                          case_sensitive::Bool=true,
                                          boxingModule::String="BoxingModule",
                                          with_body::Bool=false)::String

    function predicate(di::Documented_Item)
        ok = true
        ok = ok && ( (isempty(tag))||contains_tag(di,tag) )
        ok = ok && ( (isempty(tag_to_ignore))||(!contains_tag(di,tag_to_ignore)) )
        ok = ok && ( (identifier=="")||(identifier==J4Org.identifier(di)) )

        return ok
    end

    return org_string_documented_item_array(di_array,
                                            predicate,
                                            header_level=header_level,
                                            link_prefix=link_prefix,
                                            complete_link=complete_link,
                                            case_sensitive=case_sensitive,
                                            boxingModule=boxingModule,
                                           with_body=with_body)
end

#+API                                                                L:print_org_doc_API
#
# Prints generated documentation to be exported by OrgMode, this is the main function of the =J4Org= package.
#
# *Org-Mode Usage example:*
# #+BEGIN_SRC julia :eval never :exports code
# ,#+BEGIN_SRC julia :results output drawer :eval no-export :exports results
# documented_items = create_documented_item_array_dir("~/GitLab/MyPackage.jl/src/");
# print_org_doc(documented_items,tag="API",header_level=0)
# ,#+END_SRC
# #+END_SRC
#
# *Arguments:*
# - =tag=: tags to collect when generating the documentation
# - =tag_to_ignore=: tags to ignore when generating the documentation
# - =identifier=: generates documentation for this "identifier". Can be a function name, a structure name, etc…
# - =link_prefix=: allows to add a prefix to extra link (#+tag L=extra_link). this is can be useful to avoid link name conflict when performing local doc extraction.
# - =complete_link=: if true, try to fix link without target by adding extra items
# - =case_sensitive=: case sensitive index.
# - =boxingModule=: specifies the context in which "#!" code will be executed. See [[initialize_boxing_module][]] for details.
# - =with_body::Bool=: if true include code body
#
function print_org_doc(di_array::Array{Documented_Item,1};
                       tag::Union{String,Array{String,1}}="",
                       tag_to_ignore::Union{String,Array{String,1}}="",
                       identifier::String="",
                       header_level::Int=0,
                       link_prefix::String=randstring(),
                       complete_link::Bool=false,
                       case_sensitive::Bool=true,
                       boxingModule::String="BoxingModule",
                       with_body::Bool=false)
    print(org_string_documented_item_array(di_array,
                                           tag=tag,
                                           tag_to_ignore=tag_to_ignore,
                                           identifier=identifier,
                                           header_level=header_level,
                                           link_prefix=link_prefix,
                                           complete_link=complete_link,
                                           case_sensitive=case_sensitive,
                                           boxingModule=boxingModule,
                                           with_body=with_body))
end


text:
module Readinputfile
    include("data_module.jl")
    include("readinputfile_module.jl")
    using Base.Unicode
    using Printf
    using .Readinputfile_module

    function construct(inputfilename, usethread)
        lines = nothing
        try
            lines = readlines(inputfilename)
        catch e
            @printf "%s。プログラムを終了します。\n" sprint(showerror, e)
            exit(-1)
        end

        data = Readinputfile_module.Data_module.Data_val(
            usethread,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing)

        return Readinputfile_module.Readinputfile_variables(data, lines, 1)
    end

    errormessage(s) = let
        @printf "インプットファイルに%sの行が見つかりませんでした。\n" s
    end

    errormessage(line, s1, s2) = let
        @printf "インプットファイルの[%s]の行が正しくありません。\n" s1
        @printf "%d行目, 未知のトークン:%s\n" line s2
    end

    gettokens(rif_val, word) = let
        line = rif_val.lines[rif_val.lineindex]

        # 読み込んだ行が空、あるいはコメント行でないなら
        if !isempty(line) && line[1] != '#'
            # トークン分割
            tokens = map(line -> lowercase(line), split(line))
            w = lowercase(word)

            if tokens[1] != w
                errormessage(rif_val.lineindex, word, tokens[1])
                return -1, nothing
            end

            return 0, tokens
        else
            return 1, nothing
        end
    end

    function readchemical_number!(rif_val)
        # 原子の種類を読み込む
        chemnum = readdata!(rif_val, Data_module.CHEMICAL_NUMBER)
        if chemnum == nothing
            throw(ErrorException("インプットファイルが異常です"))
        end

        rif_val.data.Z = parse(Float64, chemnum)
    end

    function readdata!(rif_val, word)
        while true
            ret_val, tokens = gettokens(rif_val, word)

            if ret_val == -1
                return nothing
            elseif ret_val == 0
                # 読み込んだトークンの数がもし2個以外だったら
                if length(tokens) != 2
                    @printf "インプットファイル%d行目の、[%s]の行が正しくありません。\n" rif_val.lineindex word
                    return nothing
                end

                rif_val.lineindex += 1
                return tokens[2]
            end

            rif_val.lineindex += 1
        end
    end

    function readdata!(default_value, rif_val, word)
        # グリッドを読み込む
        while true
            ret_val, tokens = gettokens(rif_val, word)

            if ret_val == -1
                return nothing
            elseif ret_val == 0
                rif_val.lineindex += 1

                len = length(tokens)
                # 読み込んだトークンの数をはかる
                if len == 1
                    # デフォルト値を返す
                    return default_value
                elseif len == 2
                    if tokens[2] == "default"
                        # デフォルト値を返す
                        return default_value
                    else
                        val = nothing
                        if typeof(default_value) != String
                            val = tryparse(typeof(default_value), tokens[2])
                        else
                            val = tokens[2]
                        end

                        if val != nothing
                            return val
                        else
                            errormessage(rif_val.lineindex - 1, word, tokens[2])
                            return nothing
                        end
                    end
                else
                    str = tokens[2]

                    if str == "default" || str[1] == '#'
                        return default_value
                    elseif tokens[3][1] != '#'
                        errormessage(rif_val.lineindex - 1, word, tokens[3]);
                        return nothing
                    else
                        val = nothing
                        if typeof(default_value) != String
                            val = tryparse(typeof(default_value), tokens[2])
                        else
                            val = tokens[2]
                        end

                        if val != nothing
                            return val
                        else
                            errormessage(rif_val.lineindex - 1, word, str)
                            return nothing
                        end
                    end
                end
            end

            rif_val.lineindex += 1
        end
    end

    function readdataauto!(rif_val, word)
        while true
            ret_val, tokens = gettokens(rif_val, word)

            if ret_val == -1
                return nothing
            elseif ret_val == 0
                rif_val.lineindex += 1

                # 読み込んだトークンの数をはかる
                len = length(tokens)

                if len == 1
                    return nothing
                elseif len == 2
                    return (tokens[2] == "default" || tokens[2] == "auto") ? "" : tokens[2]
                else
                    if tokens[2] == "default" || tokens[2] == "auto" || tokens[2][0] == '#'
                        # デフォルト値を返す
                        return ""
                    elseif tokens[3][0] != '#'
                        errormessage(rif_val.lineindex - 1, word, tokens[3])

                        # エラー
                        return nothing
                    end

                    return tokens[2]
                end
            end

            rif_val.lineindex += 1
        end
    end

    function readfile(rif_val)
        # 対象の原子番号を読み込む
        readchemical_number!(rif_val)

        # グリッドの最小値を読み込む
        rif_val.data.xmin = readvalue!(Data_module.XMIN_DEFAULT, rif_val, "grid.xmin")

        # グリッドの最大値を読み込む
        rif_val.data.xmax = readvalue!(Data_module.XMAX_DEFAULT, rif_val, "grid.xmax")

        # グリッドのサイズを読み込む
        rif_val.data.grid_num = readvalue!(Data_module.GRID_NUM_DEFAULT, rif_val, "grid.num")

        # 許容誤差を読み込む
        rif_val.data.eps = readvalue!(Data_module.EPS_DEFAULT, rif_val, "eps")

        # マッチングポイントを読み込む
        rif_val.data.x_matching_point = readvalue!(Data_module.X_MATCHING_POINT_DEFAULT, rif_val, "x.matching.point")

        # 積分に使うGauss-Legendreの分点を読み込む
        rif_val.data.gauss_legendre_integ_num  = readvalue!(Data_module.GAUSS_LEGENDRE_INTEG_DEFAULT, rif_val, "gauss.legendre.integ")

        # SCFの最大ループ回数を読み込む
        rif_val.data.iteration_maxiter = readvalue!(Data_module.ITERATION_MAXITER_DEFAULT, rif_val, "iteration.maxIter")

        # SCFの一次混合の重みを読み込む
        readiterationmixingweight!(rif_val)

        # SCFの収束判定条件の値を読み込む
        rif_val.data.iteration_criterion = readvalue!(Data_module.ITERATION_CRITERION_DEFAULT, rif_val, "iteration.criterion")

        return rif_val.data
    end

    function readiterationmixingweight!(rif_val)
        rif_val.data.iteration_mixing_weight = readvalue!(Data_module.ITERATION_MIXING_WEIGHT_DEFAULT, rif_val, "iteration.Mixing.Weight")
        if rif_val.data.iteration_mixing_weight <= 0.0 || rif_val.data.iteration_mixing_weight > 1.0
            @printf "インプットファイルの[scf.Mixing.Weight]の行が正しくありません。\n"
            throw(ErrorException("インプットファイルが異常です"))
        end
    end

    function readvalue!(default_value, rif_val, word)
        val = readdata!(default_value, rif_val, word)
        if val != nothing
            return val
        else
            throw(ErrorException("インプットファイルが異常です"))
        end
    end
end


text:
#=
    We use two sets of infixable comparatives.

    The conventional symbols {==,!=,<,<=,>,>=} are defined using
       the Arb C library's implementation of eq,ne,lt,le,gt,ge.

    The predecessor/successor symbols {≃,≄,≺,≼,≻,≽} are defined from
       Hend Dawood's non-strict total ordering for interval values.
       (q.v. Hend's Master's thesis:
        Interval Mathematics Foundations, Algebraic Structures, and Applications)

=#

for (op,cop) in ((:(==), :(arb_eq)), (:(!=), :(arb_ne)),
                 (:(<=), :(arb_le)), (:(>=), :(arb_ge)),
                 (:(<), :(arb_lt)),  (:(>), :(arb_gt))  )
  @eval begin
    function ($op){T<:ArbFloat}(a::T, b::T)
        return Bool(ccall(@libarb($cop), Cint, (Ptr{T}, Ptr{T}), &a, &b) )
    end
    ($op){P,Q}(a::ArbFloat{P}, b::ArbFloat{Q}) = ($op)(promote(a,b)...)
    ($op){T<:ArbFloat,R<:Real}(a::T, b::R) = ($op)(promote(a,b)...)
    ($op){T<:ArbFloat,R<:Real}(a::R, b::T) = ($op)(promote(a,b)...)
  end
end

function (≃){T<:ArbFloat}(a::T, b::T)
    return Bool(ccall(@libarb(arb_eq), Cint, (Ptr{T}, Ptr{T}), &a, &b))
end
simeq{T<:ArbFloat}(a::T, b::T) = (≃)(a,b)

function (≄){T<:ArbFloat}(a::T, b::T)
    return !Bool(ccall(@libarb(arb_eq), Cint, (Ptr{T}, Ptr{T}), &a, &b))
end
nsime{T<:ArbFloat}(a::T, b::T) = (≄)(a,b)

function (⪰){T<:ArbFloat}(a::T, b::T)
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo < blo) || ((alo == blo) & (ahi <= bhi))
end
succeq{T<:ArbFloat}(a::T, b::T) = (⪰)(a,b)

function (≻){T<:ArbFloat}(a::T, b::T) # (a ≼ b) & (a ≄ b)
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo < blo) || ((alo == blo) & (ahi < bhi))
end
succ{T<:ArbFloat}(a::T, b::T) = (≻)(a,b)

function (⪯){T<:ArbFloat}(a::T, b::T)
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo > blo) || ((alo == blo) & (ahi >= bhi))
end
preceq{T<:ArbFloat}(a::T, b::T) = (⪯)(a,b)

function (≺){T<:ArbFloat}(a::T, b::T)
    alo, ahi = bounds(a)
    blo, bhi = bounds(b)
    return (alo > blo) || ((alo == blo) & (ahi > bhi))
end
prec{T<:ArbFloat}(a::T, b::T) = (≺)(a,b)

# for sorted ordering
isequal{T<:ArbFloat}(a::T, b::T) = !(a != b)
isless{ T<:ArbFloat}(a::T, b::T) = b == max(a,b) # !(a >= b)
isequal{T<:ArbFloat}(a::Void, b::T) = false
isequal{T<:ArbFloat}(a::T, b::Void) = false
isless{ T<:ArbFloat}(a::Void, b::T) = true
isless{ T<:ArbFloat}(a::T, b::Void) = true

# experimental ≖ ≗
(eq){T<:ArbFloat}(x::T, y::T) = !(x != y)
(eq){P,Q}(x::ArbFloat{P}, y::ArbFloat{Q}) = (eq)(promote(x,y)...)
(eq){T1<:ArbFloat,T2<:Real}(x::T1, y::T2) = (eq)(promote(x,y)...)
(eq){T1<:ArbFloat,T2<:Real}(x::T2, y::T1) = (eq)(promote(x,y)...)
(≗){T<:ArbFloat}(x::T, y::T) = !(x != y)
(≗){P,Q}(x::ArbFloat{P}, y::ArbFloat{Q}) = (≗)(promote(x,y)...)
(≗){T1<:ArbFloat,T2<:Real}(x::T1, y::T2) = (≗)(promote(x,y)...)
(≗){T1<:ArbFloat,T2<:Real}(x::T2, y::T1) = (≗)(promote(x,y)...)
(neq){T<:ArbFloat}(x::T, y::T) = donotoverlap(x, y)
(neq){P,Q}(x::ArbFloat{P}, y::ArbFloat{Q}) = (donotoverlap)(promote(x,y)...)
(neq){T1<:ArbFloat,T2<:Real}(x::T1, y::T2) = (donotoverlap)(promote(x,y)...)
(neq){T1<:ArbFloat,T2<:Real}(x::T2, y::T1) = (donotoverlap)(promote(x,y)...)


text:
"""
AP <: ProjectionsMethod
Classical Alternating projection method

Project on one set, than on another, stop after fixed number of iterations or if the desired accuracy is achieved
"""
abstract type AP <: ProjectionsMethod end

Base.@kwdef struct APparam <: AP
    x⁰ = missing
    maxϵ::Union{Float64, Missing} = missing
    maxit::Union{Missing,Int64} = missing
    keephistory::Bool = false
    snapshots::Array{Int64} = Int64[]
end

# APparam() = APparam(missing,missing,missing, false, Int64[])

initial(alg::APparam) = alg.x⁰
tolerance(alg::APparam) = alg.maxϵ
maxit(alg::APparam) = alg.maxit
keephistory(alg::APparam) = alg.keephistory
snapshots(alg::APparam) = alg.snapshots

function solve(p::TwoSetsFP, alg::APparam, x⁰, maxϵ, maxit, keephistory::Bool, snapshots::Vector{Int64})
    A = p.A
    B = p.B

    # process the default parameters
    !ismissing(x⁰) || ( x⁰ = getelement(A) )
    !ismissing(maxϵ) || (maxϵ = 1e-15)
    !ismissing(maxit) || (maxit = 100)

    k = 0
    ϵ = Inf

    xᵏ = copy(x⁰)
    xᵏ⁺¹ = similar(xᵏ)
    x̃ᵏ⁺¹ = similar(xᵏ)
    ỹᵏ = similar(xᵏ)
    yᵏ = similar(xᵏ)

    #  preallocated error vector
    err = similar(xᵏ)

    if keephistory
        errhist = Vector{Float64}(undef, maxit)
        disthist = Vector{Float64}(undef, maxit)
    else
        errhist = Float64[]
        disthist = Float64[]
    end

    if length(snapshots) != 0
        xhist = [copy(x⁰) for i in 1:length(snapshots)]
    else
        xhist = typeof(x⁰)[]
    end
    j = 1

    while k < maxit && ϵ > maxϵ

        project!(yᵏ, xᵏ, B)
        project!(xᵏ⁺¹, yᵏ, A)

        err .= xᵏ⁺¹ .- xᵏ # This doesn't say much in infeasible case, but is OK in case of binary aperture
        # dist .= xᵏ⁺¹ .- yᵏ # this calculates true error but can stay large in case of infeasible case
        ϵ = LinearAlgebra.norm(err)
        xᵏ .= xᵏ⁺¹
        k += 1

    #         println(ϵ)
        if keephistory
            errhist[k] = ϵ
            disthist[k] = LinearAlgebra.norm(xᵏ⁺¹ .- yᵏ)
        end

        if k ∈ snapshots
            println("Saving snapshot # $j, iteration # $k")
            xhist[j] .= xᵏ
             j += 1
        end

    end

    println("Using $(supertype(typeof(alg))):  to converge with $ϵ accuracy, it took me $k iterations")
    # if keephistory
    #     if length(snapshots) != 0
    #         return xᵏ, errhist, xhist
    #     else
    #         return xᵏ, errhist, Vector{T, (0,)}
    #     end
    # else
    #     return xᵏ, Vector{Float64, (0,)}, Vector{T, (0,)}
    # end
    return xᵏ, (lasty = yᵏ, errhist = errhist, xhist = xhist[1:j - 1], disthist = disthist, k= k)

end

export AP, APparam


text:
module test_julia_iterators_peel

using Test

(a, rest) = Iterators.peel("abc")
@test a == 'a'
@test collect(rest) == ['b', 'c']

@test collect(Iterators.product((1,2), (3,4))) == [(1, 3) (1, 4);  (2, 3)  (2, 4)]

end # module test_julia_iterators_peel

module test_julia_iterators_repeated

using Test
using Base.Iterators: repeated

@test collect(repeated((1,2), 3)) == [(1, 2), (1, 2), (1, 2)]

end # module test_julia_iterators_repeated


text:
# Convert a prolongation matrix P to a tentative operator S.
# P is piecewise constant over some number of aggregates.
# S has a n x n for each aggregate where n is the size of the aggregate.

using PETScBinaryIO
using SparseArrays
using LinearAlgebra

function p_to_s(P)
    rows = rowvals(P)
    vals = nonzeros(P)
    m, n = size(P)
    is = Vector{Int}()
    js = Vector{Int}()
    ks = Vector{Float64}()
    j = 0
    bs = 16
    for i = 1:bs:n
        # get row indices of the aggregate
        inds = sort(unique(vcat(map(x -> rows[nzrange(P, i+x)], 0:(bs-1))...)))
        agg_size = length(inds)
        if agg_size > 10000
            println("Aggregate is too large ($agg_size)")
            exit(-1)
        end
        if agg_size == 1
            @warn "Ignoring aggregate of size 1"
            continue
        end
        B = Matrix(P[inds, i:(i+bs-1)])

        # I - B(BᵀB)⁻¹Bᵀ
        block = I - B * inv(B' * B) * B'
        # block = I - fill(1/agg_size, (agg_size, agg_size)) # B(BᵀB)⁻¹Bᵀ for piecewise constant

        U, S, _ = svd(block)
        l = findlast(x -> x > 1e-14, S)
        for x = 1:agg_size
            for y = 1:l
                push!(is, inds[x])
                push!(js, j + y)
                push!(ks, U[x,y])
            end
        end
        j += l

        # verify that S projects out the coarse space
        @assert norm(U[:, 1:l]' * B, Inf) < 1e14
    end
    sparse(is,js,ks,m,j)
end

if joinpath(pwd(), PROGRAM_FILE) == @__FILE__
    P = readpetsc(ARGS[1])[1]
    writepetsc(ARGS[2], p_to_s(P))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.LineEdit
isdefined(Main, :TestHelpers) || @eval Main include(joinpath(dirname(@__FILE__), "TestHelpers.jl"))
using TestHelpers

function run_test(d,buf)
    global a_foo, b_foo, a_bar, b_bar
    a_foo = b_foo = a_bar = b_bar = 0
    while !eof(buf)
        LineEdit.match_input(d, nothing, buf)(nothing,nothing)
    end
end

a_foo = 0

const foo_keymap = Dict(
    'a' => (o...)->(global a_foo; a_foo += 1)
)

b_foo = 0

const foo2_keymap = Dict(
    'b' => (o...)->(global b_foo; b_foo += 1)
)

a_bar = 0
b_bar = 0

const bar_keymap = Dict(
    'a' => (o...)->(global a_bar; a_bar += 1),
    'b' => (o...)->(global b_bar; b_bar += 1)
)

test1_dict = LineEdit.keymap([foo_keymap])

run_test(test1_dict,IOBuffer("aa"))
@test a_foo == 2

test2_dict = LineEdit.keymap([foo2_keymap, foo_keymap])

run_test(test2_dict,IOBuffer("aaabb"))
@test a_foo == 3
@test b_foo == 2

test3_dict = LineEdit.keymap([bar_keymap, foo_keymap])

run_test(test3_dict,IOBuffer("aab"))
@test a_bar == 2
@test b_bar == 1

# Multiple spellings in the same keymap
const test_keymap_1 = Dict(
    "^C" => (o...)->1,
    "\\C-C" => (o...)->2
)

@test_throws ErrorException LineEdit.keymap([test_keymap_1])

a_foo = a_bar = 0

const test_keymap_2 = Dict(
    "abc" => (o...)->(global a_foo = 1)
)

const test_keymap_3 = Dict(
    "a"  => (o...)->(global a_foo = 2),
    "bc" => (o...)->(global a_bar = 3)
)

function keymap_fcn(keymaps)
    d = LineEdit.keymap(keymaps)
    f = buf->(LineEdit.match_input(d, nothing, buf)(nothing,nothing))
end

let f = keymap_fcn([test_keymap_3, test_keymap_2])
    buf = IOBuffer("abc")
    f(buf); f(buf)
    @test a_foo == 2
    @test a_bar == 3
    @test eof(buf)
end

# Eager redirection when the redirected-to behavior is changed.

a_foo = 0

const test_keymap_4 = Dict(
    "a" => (o...)->(global a_foo = 1),
    "b" => "a",
    "c" => (o...)->(global a_foo = 2),
)

const test_keymap_5 = Dict(
    "a" => (o...)->(global a_foo = 3),
    "d" => "c"
)

let f = keymap_fcn([test_keymap_5, test_keymap_4])
    buf = IOBuffer("abd")
    f(buf)
    @test a_foo == 3
    f(buf)
    @test a_foo == 1
    f(buf)
    @test a_foo == 2
    @test eof(buf)
end

# Eager redirection with cycles

const test_cycle = Dict(
    "a" => "b",
    "b" => "a"
)

@test_throws ErrorException keymap([test_cycle])

# Lazy redirection with Cycles

const level1 = Dict(
    "a" => LineEdit.KeyAlias("b")
)

const level2a = Dict(
    "b" => "a"
)

const level2b = Dict(
    "b" => LineEdit.KeyAlias("a")
)

@test_throws ErrorException keymap([level2a,level1])
@test_throws ErrorException keymap([level2b,level1])

# Lazy redirection functionality test

a_foo = 0

const test_keymap_6 = Dict(
    "a" => (o...)->(global a_foo = 1),
    "b" => LineEdit.KeyAlias("a"),
    "c" => (o...)->(global a_foo = 2),
)

const test_keymap_7 = Dict(
    "a" => (o...)->(global a_foo = 3),
    "d" => "c"
)

let f = keymap_fcn([test_keymap_7, test_keymap_6])
    buf = IOBuffer("abd")
    f(buf)
    @test a_foo == 3
    a_foo = 0
    f(buf)
    @test a_foo == 3
    f(buf)
    @test a_foo == 2
    @test eof(buf)
end

# Test requiring postprocessing (see conflict fixing in LineEdit.jl )

global path1 = 0
global path2 = 0
global path3 = 0

const test_keymap_8 = Dict(
    "**" => (o...)->(global path1 += 1),
    "ab" => (o...)->(global path2 += 1),
    "cd" => (o...)->(global path3 += 1),
    "d" => (o...)->(error("This is not the key you're looking for"))
)

let f = keymap_fcn([test_keymap_8])
    buf = IOBuffer("bbabaccd")
    f(buf)
    @test path1 == 1
    f(buf)
    @test path2 == 1
    f(buf)
    @test path1 == 2
    f(buf)
    @test path3 == 1
    @test eof(buf)
end

global path1 = 0
global path2 = 0

const test_keymap_9 = Dict(
    "***" => (o...)->(global path1 += 1),
    "*a*" => (o...)->(global path2 += 1)
)

let f = keymap_fcn([test_keymap_9])
    buf = IOBuffer("abaaaa")
    f(buf)
    @test path1 == 1
    f(buf)
    @test path2 == 1
    @test eof(buf)
end

## edit_move{left,right} ##
buf = IOBuffer("a\na\na\n")
seek(buf, 0)
for i = 1:6
    LineEdit.edit_move_right(buf)
    @test position(buf) == i
end
@test eof(buf)
for i = 5:0
    LineEdit.edit_move_left(buf)
    @test position(buf) == i
end

# skip unicode combining characters
buf = IOBuffer("ŷ")
seek(buf, 0)
LineEdit.edit_move_right(buf)
@test eof(buf)
LineEdit.edit_move_left(buf)
@test position(buf) == 0

## edit_move_{up,down} ##

buf = IOBuffer("type X\n    a::Int\nend")
for i = 0:6
    seek(buf,i)
    @test !LineEdit.edit_move_up(buf)
    @test position(buf) == i
    seek(buf,i)
    @test LineEdit.edit_move_down(buf)
    @test position(buf) == i+7
end
for i = 7:17
    seek(buf,i)
    @test LineEdit.edit_move_up(buf)
    @test position(buf) == min(i-7,6)
    seek(buf,i)
    @test LineEdit.edit_move_down(buf)
    @test position(buf) == min(i+11,21)
end
for i = 18:21
    seek(buf,i)
    @test LineEdit.edit_move_up(buf)
    @test position(buf) == i-11
    seek(buf,i)
    @test !LineEdit.edit_move_down(buf)
    @test position(buf) == i
end

buf = IOBuffer("type X\n\n")
seekend(buf)
@test LineEdit.edit_move_up(buf)
@test position(buf) == 7
@test LineEdit.edit_move_up(buf)
@test position(buf) == 0
@test !LineEdit.edit_move_up(buf)
@test position(buf) == 0
seek(buf,0)
@test LineEdit.edit_move_down(buf)
@test position(buf) == 7
@test LineEdit.edit_move_down(buf)
@test position(buf) == 8
@test !LineEdit.edit_move_down(buf)
@test position(buf) == 8

## edit_delete_prev_word ##

buf = IOBuffer("type X\n ")
seekend(buf)
@test LineEdit.edit_delete_prev_word(buf)
@test position(buf) == 5
@test buf.size == 5
@test String(buf.data[1:buf.size]) == "type "

buf = IOBuffer("4 +aaa+ x")
seek(buf,8)
@test LineEdit.edit_delete_prev_word(buf)
@test position(buf) == 3
@test buf.size == 4
@test String(buf.data[1:buf.size]) == "4 +x"

buf = IOBuffer("x = func(arg1,arg2 , arg3)")
seekend(buf)
LineEdit.char_move_word_left(buf)
@test position(buf) == 21
@test LineEdit.edit_delete_prev_word(buf)
@test String(buf.data[1:buf.size]) == "x = func(arg1,arg3)"
@test LineEdit.edit_delete_prev_word(buf)
@test String(buf.data[1:buf.size]) == "x = func(arg3)"
@test LineEdit.edit_delete_prev_word(buf)
@test String(buf.data[1:buf.size]) == "x = arg3)"

# Unicode combining characters
let buf = IOBuffer()
    LineEdit.edit_insert(buf, "â")
    LineEdit.edit_move_left(buf)
    @test position(buf) == 0
    LineEdit.edit_move_right(buf)
    @test nb_available(buf) == 0
    LineEdit.edit_backspace(buf)
    @test String(buf.data[1:buf.size]) == "a"
end

## edit_transpose ##
let buf = IOBuffer()
    LineEdit.edit_insert(buf, "abcde")
    seek(buf,0)
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "abcde"
    LineEdit.char_move_right(buf)
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "bacde"
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "bcade"
    seekend(buf)
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "bcaed"
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "bcade"

    seek(buf, 0)
    LineEdit.edit_clear(buf)
    LineEdit.edit_insert(buf, "αβγδε")
    seek(buf,0)
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "αβγδε"
    LineEdit.char_move_right(buf)
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "βαγδε"
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "βγαδε"
    seekend(buf)
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "βγαεδ"
    LineEdit.edit_transpose(buf)
    @test String(buf.data[1:buf.size]) == "βγαδε"
end

let
    term = TestHelpers.FakeTerminal(IOBuffer(), IOBuffer(), IOBuffer())
    s = LineEdit.init_state(term, ModalInterface([Prompt("test> ")]))
    buf = LineEdit.buffer(s)

    LineEdit.edit_insert(s,"first line\nsecond line\nthird line")
    @test String(buf.data[1:buf.size]) == "first line\nsecond line\nthird line"

    ## edit_move_line_start/end ##
    seek(buf, 0)
    LineEdit.move_line_end(s)
    @test position(buf) == sizeof("first line")
    LineEdit.move_line_end(s) # Only move to input end on repeated keypresses
    @test position(buf) == sizeof("first line")
    s.key_repeats = 1 # Manually flag a repeated keypress
    LineEdit.move_line_end(s)
    s.key_repeats = 0
    @test eof(buf)

    seekend(buf)
    LineEdit.move_line_start(s)
    @test position(buf) == sizeof("first line\nsecond line\n")
    LineEdit.move_line_start(s)
    @test position(buf) == sizeof("first line\nsecond line\n")
    s.key_repeats = 1 # Manually flag a repeated keypress
    LineEdit.move_line_start(s)
    s.key_repeats = 0
    @test position(buf) == 0

    ## edit_kill_line, edit_yank ##
    seek(buf, 0)
    LineEdit.edit_kill_line(s)
    s.key_repeats = 1 # Manually flag a repeated keypress
    LineEdit.edit_kill_line(s)
    s.key_repeats = 0
    @test String(buf.data[1:buf.size]) == "second line\nthird line"
    LineEdit.move_line_end(s)
    LineEdit.edit_move_right(s)
    LineEdit.edit_yank(s)
    @test String(buf.data[1:buf.size]) == "second line\nfirst line\nthird line"
end

# Issue 7845
# First construct a problematic string:
# julia> is 6 characters + 1 character for space,
# so the rest of the terminal is 73 characters
#########################################################################
let
    buf = IOBuffer(
    "begin\nprint(\"A very very very very very very very very very very very very ve\")\nend")
    seek(buf,4)
    outbuf = IOBuffer()
    termbuf = Base.Terminals.TerminalBuffer(outbuf)
    term = TestHelpers.FakeTerminal(IOBuffer(), IOBuffer(), IOBuffer())
    s = LineEdit.refresh_multi_line(termbuf, term, buf,
        Base.LineEdit.InputAreaState(0,0), "julia> ", indent = 7)
    @test s == Base.LineEdit.InputAreaState(3,1)
end

@testset "function prompt indentation" begin
    term = TestHelpers.FakeTerminal(IOBuffer(), IOBuffer(), IOBuffer(), false)
    # default prompt: PromptState.indent should not be set to a final fixed value
    s = LineEdit.init_state(term, ModalInterface([Prompt("julia> ")]))
    ps::LineEdit.PromptState = s.mode_state[s.current_mode]
    @test ps.indent == -1
    # the prompt is modified afterwards to a function
    ps.p.prompt = let i = 0
        () -> ["Julia is Fun! > ", "> "][mod1(i+=1, 2)] # lengths are 16 and 2
    end
    buf = LineEdit.buffer(ps)
    write(buf, "begin\n    julia = :fun\nend")
    outbuf = IOBuffer()
    termbuf = Base.Terminals.TerminalBuffer(outbuf)
    LineEdit.refresh_multi_line(termbuf, term, ps)
    @test String(take!(outbuf)) ==
        "\r\e[0K\e[1mJulia is Fun! > \e[0m\r\e[16Cbegin\n" *
        "\r\e[16C    julia = :fun\n" *
        "\r\e[16Cend\r\e[19C"
    LineEdit.refresh_multi_line(termbuf, term, ps)
    @test String(take!(copy(outbuf))) ==
        "\r\e[0K\e[1A\r\e[0K\e[1A\r\e[0K\e[1m> \e[0m\r\e[2Cbegin\n" *
        "\r\e[2C    julia = :fun\n" *
        "\r\e[2Cend\r\e[5C"
end

@testset "tab/backspace alignment feature" begin
    term = TestHelpers.FakeTerminal(IOBuffer(), IOBuffer(), IOBuffer())
    s = LineEdit.init_state(term, ModalInterface([Prompt("test> ")]))
    function bufferdata(s)
        buf = LineEdit.buffer(s)
        String(buf.data[1:buf.size])
    end
    move_left(s, n) = for x = 1:n
        LineEdit.edit_move_left(s)
    end

    bufpos(s::Base.LineEdit.MIState) = position(LineEdit.buffer(s))

    LineEdit.edit_insert(s, "for x=1:10\n")
    LineEdit.edit_tab(s)
    @test bufferdata(s) == "for x=1:10\n    "
    LineEdit.edit_backspace(s, true, false)
    @test bufferdata(s) == "for x=1:10\n"
    LineEdit.edit_insert(s, "  ")
    @test bufpos(s) == 13
    LineEdit.edit_tab(s)
    @test bufferdata(s) == "for x=1:10\n    "
    LineEdit.edit_insert(s, "  ")
    LineEdit.edit_backspace(s, true, false)
    @test bufferdata(s) == "for x=1:10\n    "
    LineEdit.edit_insert(s, "éé=3   ")
    LineEdit.edit_tab(s)
    @test bufferdata(s) == "for x=1:10\n    éé=3    "
    LineEdit.edit_backspace(s, true, false)
    @test bufferdata(s) == "for x=1:10\n    éé=3"
    LineEdit.edit_insert(s, "\n    1∉x  ")
    LineEdit.edit_tab(s)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n    1∉x     "
    LineEdit.edit_backspace(s, false, false)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n    1∉x    "
    LineEdit.edit_backspace(s, true, false)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n    1∉x "
    LineEdit.edit_move_word_left(s)
    LineEdit.edit_tab(s)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n        1∉x "
    LineEdit.move_line_start(s)
    @test bufpos(s) == 22
    LineEdit.edit_tab(s, true)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n        1∉x "
    @test bufpos(s) == 30
    LineEdit.edit_move_left(s)
    @test bufpos(s) == 29
    LineEdit.edit_backspace(s, true, true)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n    1∉x "
    @test bufpos(s) == 26
    LineEdit.edit_tab(s, false) # same as edit_tab(s, true) here
    @test bufpos(s) == 30
    move_left(s, 6)
    @test bufpos(s) == 24
    LineEdit.edit_backspace(s, true, true)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n    1∉x "
    @test bufpos(s) == 22
    LineEdit.edit_kill_line(s)
    LineEdit.edit_insert(s, ' '^10)
    move_left(s, 7)
    @test bufferdata(s) == "for x=1:10\n    éé=3\n          "
    @test bufpos(s) == 25
    LineEdit.edit_tab(s, true, false)
    @test bufpos(s) == 32
    move_left(s, 7)
    LineEdit.edit_tab(s, true, true)
    @test bufpos(s) == 26
    @test bufferdata(s) == "for x=1:10\n    éé=3\n    "
    # test again the same, when there is a next line
    LineEdit.edit_insert(s, "      \nend")
    move_left(s, 11)
    @test bufpos(s) == 25
    LineEdit.edit_tab(s, true, false)
    @test bufpos(s) == 32
    move_left(s, 7)
    LineEdit.edit_tab(s, true, true)
    @test bufpos(s) == 26
    @test bufferdata(s) == "for x=1:10\n    éé=3\n    \nend"
end

max_stars_repo_name: nextjournal/julia-lang

max_issues_repo_name: nextjournal/julia-lang

max_forks_repo_name: nextjournal/julia-lang


text:
#=
Recursively instantiate dependent objects and propagate value in hierarchical NamedTuples

* Developer: Hilding Elmqvist, Mogram AB (initial version)
*            Martin Otter, DLR (instantiation of dependent objects added)
* First version: March 2021
* License: MIT (expat)

=#

using OrderedCollections: OrderedDict

subst(ex, environment, modelModule) = ex
subst(ex::Symbol, environment, modelModule) = if length(environment) == 0; ex elseif ex in keys(environment[end]); environment[end][ex] else subst(ex, environment[1:end-1], modelModule) end
subst(ex::Expr, environment, modelModule) = Expr(ex.head, [subst(a, environment, modelModule) for a in ex.args]...)
subst(ex::Vector{Symbol}, environment, modelModule) = [subst(a, environment, modelModule) for a in ex]
subst(ex::Vector{Expr}  , environment, modelModule) = [Core.eval(modelModule, subst(a, environment, modelModule)) for a in ex]

#=
function propagate(model, environment=[])
    #println("\n... Propagate: ", model)
    current = OrderedDict()
    for (k,v) in zip(keys(model), model)
        if k in [:_class, :_Type]
            current[k] = v
        elseif typeof(v) <: NamedTuple
            current[k] = propagate(v, vcat(environment, [current]))
        else
            evalv = v
            try
                evalv = eval(subst(v, vcat(environment, [current])))
            catch e
            end

            if typeof(evalv) <: NamedTuple
                current[k] = v
            else
                current[k] = evalv
            end
        end
    end
    return (; current...)
end
=#

appendKey(path, key) = path == "" ? string(key) : path * "." * string(key)

"""
    map = propagateEvaluateAndInstantiate!(FloatType, unitless::Bool, modelModule::Module, parameters,
                   eqInfo::Modia.EquationInfo; log=false)

Recursively traverse the hierarchical collection `parameters` and perform the following actions:

- Propagate values.
- Evaluate expressions in the context of `modelModule`.
- Instantiate dependent objects.
- Return the evaluated `parameters` if successfully evaluated, and otherwise
  return nothing, if an error occurred (an error message was printed).
"""
function propagateEvaluateAndInstantiate!(FloatType, unitless::Bool, modelModule, parameters, eqInfo, previous_dict, previous, pre_dict, pre, hold_dict, hold; log=false)
    x_found = fill(false, length(eqInfo.x_info))
    map = propagateEvaluateAndInstantiate2!(FloatType, unitless, modelModule, parameters, eqInfo, x_found, previous_dict, previous, pre_dict, pre, hold_dict, hold, [], ""; log=log)

    if isnothing(map)
        return nothing
    end

    # Check that all values of x_start are set:
    #x_start_missing = []
    #for (i, found) in enumerate(x_found)
    #    if !found
    #        push!(x_start_missing, eqInfo.x_info[i].x_name)
    #    end
    #end

    # Check that all previous values are set:
    missingInitValues = false
    namesOfMissingValues = ""
    first = true
    for (name,index) in previous_dict
        if ismissing(previous[index])
            missingInitValues = true
            if first
                first = false
                namesOfMissingValues *= "\n  Variables from previous(..):"
            end
            namesOfMissingValues *= "\n    " * name
        end
    end

    # Check that all pre values are set:
    first = true
    for (name,index) in pre_dict
        if ismissing(pre[index])
            missingInitValues = true
            if first
                first = false
                namesOfMissingValues *= "\n  Variables from pre(..):"
            end
            namesOfMissingValues *= "\n    " * name
        end
    end

    # Check that all hold values are set:
    first = true
    for (name,index) in hold_dict
        if ismissing(hold[index])
            missingInitValues = true
            if first
                first = false
                namesOfMissingValues *= "\n  Variables from hold(..):"
            end
            namesOfMissingValues *= "\n    " * name
        end
    end

    if missingInitValues
        printstyled("Model error: ", bold=true, color=:red)
        printstyled("Missing start/init values for variables: ", namesOfMissingValues,
                    bold=true, color=:red)
        print("\n\n")
        return nothing
    end

    #if length(x_start_missing) > 0
    #    printstyled("Model error: ", bold=true, color=:red)
    #    printstyled("Missing start/init values for variables: ", x_start_missing,
    #                bold=true, color=:red)
    #    print("\n\n")
    #    return nothing
    #end
    return map
end

"""
    firstName(ex::Expr)

If ex = :(a.b.c.d) -> firstName(ex) = :a
"""
function firstName(ex::Expr)
    if ex.head == :(.)
        if typeof(ex.args[1]) == Expr
            firstName(ex.args[1])
        else
            return ex.args[1]
        end
    end
end

function changeDotToRef(ex)
    if ex.head == :(.)
        ex.head = :ref
        if typeof(ex.args[1]) == Expr
            changeDotToRef(ex.args[1])
        end
    end
    return nothing
end

function propagateEvaluateAndInstantiate2!(FloatType, unitless::Bool, modelModule, parameters, eqInfo::Modia.EquationInfo, x_found::Vector{Bool},
                                           previous_dict, previous, pre_dict, pre, hold_dict, hold,
                                           environment, path::String; log=false)
    if log
        println("\n 1: !!! instantiate objects of $path: ", parameters)
    end
    current = OrderedDict{Symbol,Any}()   # should be Map()

    # Determine, whether "parameters" has a ":_constructor" key and handle this specially
    constructor = nothing
    usePath     = false
    if haskey(parameters, :_constructor)
        # For example: obj = (_class = :Par, _constructor = :(Modia3D.Object3D), _path = true, kwargs...)
        #          or: rev = (_constructor = (_class = :Par, value = :(Modia3D.ModiaRevolute), _path=true), kwargs...)
        v = parameters[:_constructor]
        if typeof(v) <: OrderedDict
            constructor = v[:value]
            if haskey(v, :_path)
                usePath = v[:_path]
            end
        else
            constructor = v
            if haskey(parameters, :_path)
                usePath = parameters[:_path]
            end
        end

    elseif haskey(parameters, :value)
        # For example: p1 = (_class = :Var, parameter = true, value = 0.2)
        #          or: p2 = (_class = :Var, parameter = true, value = :(2*p1))
        v = parameters[:value]
        veval = Core.eval(modelModule, subst(v, vcat(environment, [current]), modelModule))
        return veval
    end

    for (k,v) in parameters
        if log
            println(" 2:    ... key = $k, value = $v")
        end
        if k == :_constructor || k == :_path || (k == :_class && !isnothing(constructor))
            if log
                println(" 3:    ... key = $k")
            end
            nothing

        elseif !isnothing(constructor) && (k == :value || k == :init || k == :start)
            error("value, init or start keys are not allowed in combination with a _constructor:\n$parameters")

        elseif typeof(v) <: OrderedDict
            if length(v) > 0
                if haskey(v, :_class) && v[:_class] == :Par && haskey(v, :value)
                    # For example: k = (_class = :Par, value = 2.0) -> k = 2.0
                    #          or: k = (_class = :Par, value = :(2*Lx - 3))   -> k = eval( 2*Lx - 3 )
                    #          or: k = (_class = :Par, value = :(bar.frame0)) -> k = ref(bar.frame0)
                    if log
                        println(" 4:    v[:value] = ", v[:value], ", typeof(v[:value]) = ", typeof(v[:value]))
                        println("        vcat(environment, [current]) = ", vcat(environment, [current]))
                    end
                    subv = subst(v[:value], vcat(environment, [current]), modelModule)
                    if log
                        println(" 5:    _class & value: $k = $subv  # before eval")
                    end
                    if typeof(subv) == Expr && subv.head == :(.)
                        if typeof(firstName(subv)) <: AbstractDict
                            changeDotToRef(subv)
                            if log
                                println(" 5b:    _class & value: $k = $subv  # before eval")
                            end
                        end
                    end
                    current[k] = Core.eval(modelModule, subv)
                    if log
                        println(" 6:                   $k = ", current[k])
                    end
                else
                    if log
                        println(" 7:    ... key = $k, v = $v")
                    end
                    # For example: k = (a = 2.0, b = :(2*Lx))
                    value = propagateEvaluateAndInstantiate2!(FloatType, unitless, modelModule, v, eqInfo, x_found, previous_dict, previous, pre_dict, pre, hold_dict, hold,
                                                              vcat(environment, [current]), appendKey(path, k); log=log)
                    if log
                        println(" 8:    ... key = $k, value = $value")
                    end
                    if isnothing(value)
                        return nothing
                    end
                    current[k] = value
                end
            end

        else
            if log
                println(" 9:    else: typeof(v) = ", typeof(v))
            end
            subv = subst(v, vcat(environment, [current]), modelModule)
            if log
                println(" 10:          $k = $subv   # before eval")
            end
            if typeof(subv) == Expr && subv.head == :(.)
                if typeof(firstName(subv)) <: AbstractDict
                    changeDotToRef(subv)
                    if log
                        println(" 10b:    _class & value: $k = $subv  # before eval")
                    end
                end
            end
            subv = Core.eval(modelModule, subv)
            if unitless && eltype(subv) <: Number
                # Remove unit
                subv = stripUnit(subv)
            end
            current[k] = subv
            if log
                println(" 11:          $k = ", current[k])
            end

            # Set x_start
            full_key = appendKey(path, k)
            if haskey(eqInfo.x_dict, full_key)
                #if log
                #    println(" 12:              (is stored in x_start)")
                #end
                j = eqInfo.x_dict[full_key]
                xe_info = eqInfo.x_info[j]
                x_value = current[k]
                len = hasParticles(x_value) ? 1 : length(x_value)
                if j <= eqInfo.nxFixedLength && len != xe_info.length
                    printstyled("Model error: ", bold=true, color=:red)
                    printstyled("Length of ", xe_info.x_name, " shall be changed from ",
                                xe_info.length, " to $len\n",
                                "This is not possible because variable has a fixed length.", bold=true, color=:red)
                    return nothing
                end
                x_found[j] = true
                xe_info.startOrInit = deepcopy(x_value)

                # Strip units from x_start
                #if xe_info.length == 1
                #    x_start[xe_info.startIndex] = deepcopy( convert(FloatType, stripUnit(x_value)) )
                #else
                #    ibeg = xe_info.startIndex - 1
                #    for i = 1:xe_info.length
                #        x_start[ibeg+i] = deepcopy( convert(FloatType, stripUnit(x_value[i])) )
                #    end
                #end

            elseif haskey(previous_dict, full_key)
                previous[ previous_dict[full_key] ] = current[k]

            elseif haskey(pre_dict, full_key)
                pre[ pre_dict[full_key] ] = current[k]

            elseif haskey(hold_dict, full_key)
                hold[ hold_dict[full_key] ] = current[k]
            end
        end
    end

    if isnothing(constructor)
        return current # (; current...)
    else
        if usePath
            obj = Core.eval(modelModule, :(FloatType = $FloatType; $constructor(; path = $path, $current...)))
        else
            obj = Core.eval(modelModule, :(FloatType = $FloatType; $constructor(; $current...)))
        end
        if log
            println(" 13:    +++ Instantiated $path: typeof(obj) = ", typeof(obj), ", obj = ", obj, "\n\n")
        end
        return obj
    end
end

"""
    splittedPath = spitPlath(path::Union{Symbol, Expr, Nothing})::Vector{Symbol}

# Examples
```
splitPath(nothing)     # = Symbol[]
splitPath(:a)          # = Symbol[:a]
splitPath(:(a.b.c.d))  # = Symbol[:a, :b, :c, :d]
```
"""
function splitPath(path::Union{Symbol, Expr, Nothing})::Vector{Symbol}
    splittedPath = Symbol[]
    if typeof(path) == Symbol
        push!(splittedPath, path)
    elseif typeof(path) == Expr
        while typeof(path.args[1]) == Expr
            pushfirst!(splittedPath, path.args[2].value)
            path = path.args[1]
        end
        pushfirst!(splittedPath, path.args[2].value)
        pushfirst!(splittedPath, path.args[1])
    end
    return splittedPath
end

"""
    modelOfPath = getModelFromSplittedPath(model, splittedPath::Vector{Symbol})

Return reference to the sub-model characterized by `splittedPath`.

# Examples

```
mymodel = Model(a = Model(b = Model(c = Model)))
model_b = getModelFromSplittedPath(mymodel, Symbol["a", "b"])  # = mymodel[:a][:b]
```
"""
function getModelFromSplittedPath(model, splittedPath::Vector{Symbol})
    for name in splittedPath
        model = model[name]
    end
    return model
end


text:
using Mimi
using Base.Test

tests = ["main", "references", "units", "model_structure"]

for t in tests
    fp = joinpath("test_$t.jl")
    println("$fp ...")
    include(fp)
end


text:
using Test
using kNN

@testset "MajorityVote" begin
	votes = ["A", "B", "B", "A", "C", "B"]
    @test isequal(kNN.majority_vote(votes), "B")

    ["A", "B", "B", "A", "C", "B", "C", "C"]
    @test isequal(kNN.majority_vote(votes), "B")

    votes = ["A", "B", "B", "A", "C", "B", "C", "C", "C"]
    @test isequal(kNN.majority_vote(votes), "C")
end


text:
module NaNStatistics

    # AVX vectorziation tools
    using IfElse: ifelse
    using LoopVectorization
    using VectorizationBase: Vec

    using Static
    _dim(::Type{StaticInt{N}}) where {N} = N::Int

    include("ArrayStats/ArrayStats.jl")
    include("ArrayStats/nanmean.jl")
    include("ArrayStats/nansum.jl")
    include("ArrayStats/nanvar.jl")
    include("ArrayStats/nanstd.jl")
    include("ArrayStats/nancov.jl")
    include("Sorting/quicksort.jl")
    include("Sorting/nanmedian.jl")
    include("Sorting/nanpctile.jl")
    include("Histograms.jl")
    include("Binning.jl")

end