text:
###### rename, rename! ####################

"""
    rename(ta::TimeArray, colnames::Vector{Symbol})
    rename(ta::TimeArray, colname::Symbol)
    rename(ta::TimeArray, orig => new, ...)
    rename(f::Base.Callable, ta, colnametyp)

Rename the columns of a `TimeArray`.

See also [`rename!`](@ref).

# Arguments
- `colnametyp` is the input type for the function `f`. The valid value is
  `Symbol` or `String`.
"""
function rename(ta::TimeArray, colnames::Vector{Symbol})
    length(colnames) == size(ta, 2) || throw(ArgumentError("Colnames length mismatch"))
    TimeArray(ta, colnames = colnames, unchecked = true)
end

"""
    rename!(ta::TimeArray, colnames::Vector{Symbol})
    rename!(ta::TimeArray, colname::Symbol)
    rename!(ta::TimeArray, orig => new, ...)
    rename!(f::Base.Callable, ta, colnametyp)

In-place rename the columns of a `TimeArray`.

See also [`rename`](@ref).

# Arguments
- `colnametyp` is the input type for the function `f`. The valid value is
  `Symbol` or `String`.
"""
function rename!(ta::TimeArray, colnames::Vector{Symbol})
    length(colnames) == size(ta, 2) ? (_colnames(ta)[:] = colnames) :
    length(colnames) > 0 && throw(ArgumentError("Colnames length mismatch"))
    ta
end

# shared interfaces
for f ∈ [:rename, :rename!]
    @eval begin
        $f(ta::TimeArray, colnames::Symbol) = $f(ta, [colnames])
        $f(ta::TimeArray) = throw(MethodError($f, (ta,)))
        $f(ta::TimeArray, pairs::Pair{Symbol,Symbol}...) =
            $f(ta, _mapcol(colnames(ta), pairs))
        $f(f::Base.Callable, ta::TimeArray, colnametyp::Type{Symbol} = Symbol) =
            $f(ta, map(f, colnames(ta)))
        $f(f::Base.Callable, ta::TimeArray, colnametyp::Type{String}) =
            $f(Symbol ∘ f ∘ string, ta)
    end
end

"""
utility function for `rename` and `rename!`
"""
function _mapcol(col, pairs)
    col′ = copy(col)
    for (key, val) ∈ pairs
        i = findfirst(isequal(key), col)
        if i ≡ nothing  # FIXME: `isnothing(i)` in Julia 1.1
            throw(ArgumentError("Unknown column `$key`"))
        end
        col′[i] = val
    end
    col′
end


text:
export BrickletDustDetectorDustDensityCallbackThreshold
struct BrickletDustDetectorDustDensityCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

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

export BrickletDustDetector
"""
Measures dust density
"""
mutable struct BrickletDustDetector <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_dust_density
"""
    $(SIGNATURES)

Returns the dust density.

If you want to get the dust density periodically, it is recommended
to use the :cb:`Dust Density` callback and set the period with
:func:`Set Dust Density Callback Period`.
"""
function get_dust_density(device::BrickletDustDetector)
    return device.deviceInternal.get_dust_density()
end

export set_dust_density_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Dust Density` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Dust Density` callback is only triggered if the dust density has
changed since the last triggering.
"""
function set_dust_density_callback_period(device::BrickletDustDetector, period)
    device.deviceInternal.set_dust_density_callback_period(period)
end

export get_dust_density_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Dust Density Callback Period`.
"""
function get_dust_density_callback_period(device::BrickletDustDetector)
    return device.deviceInternal.get_dust_density_callback_period()
end

export set_dust_density_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Dust Density Reached` callback.

The following options are possible:

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

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the dust density value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the dust density value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the dust density value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the dust density value is greater than the min value (max is ignored)"
"""
function set_dust_density_callback_threshold(device::BrickletDustDetector, option, min, max)
    device.deviceInternal.set_dust_density_callback_threshold(option, min, max)
end

export get_dust_density_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Dust Density Callback Threshold`.
"""
function get_dust_density_callback_threshold(device::BrickletDustDetector)
    return device.deviceInternal.get_dust_density_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Dust Density Reached`

is triggered, if the threshold

* :func:`Set Dust Density Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletDustDetector, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletDustDetector)
    return device.deviceInternal.get_debounce_period()
end

export set_moving_average
"""
    $(SIGNATURES)

Sets the length of a `moving averaging <https://en.wikipedia.org/wiki/Moving_average>`__
for the dust density.

Setting the length to 0 will turn the averaging completely off. With less
averaging, there is more noise on the data.
"""
function set_moving_average(device::BrickletDustDetector, average)
    device.deviceInternal.set_moving_average(average)
end

export get_moving_average
"""
    $(SIGNATURES)

Returns the length moving average as set by :func:`Set Moving Average`.
"""
function get_moving_average(device::BrickletDustDetector)
    return device.deviceInternal.get_moving_average()
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::BrickletDustDetector)
    return device.deviceInternal.get_identity()
end

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

ext: jl


text:
using Test
using LoopVectorization
using LinearAlgebra

# const START_TIME = time()
# exceeds_time_limit() = (time() - START_TIME) > 35 * 60

function clenshaw(x,coeff)
    len_c = length(coeff)
    tmp = zero(x)
    ret = zero(x)
    for i in len_c:-1:2
        ret     = muladd(x,2tmp,coeff[i]-ret)
        ret,tmp = tmp,ret
    end
    ret = muladd(x,tmp,coeff[1]-ret)
    return ret
end

"""
Causes `check_args` to fail.
"""
struct FallbackArrayWrapper{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
end
Base.size(A::FallbackArrayWrapper) = size(A.data)
Base.@propagate_inbounds Base.getindex(A::FallbackArrayWrapper, i::Vararg{Int, N}) where {N} = getindex(A.data, i...)
Base.@propagate_inbounds Base.setindex!(A::FallbackArrayWrapper, v, i::Vararg{Int, N}) where {N} = setindex!(A.data, v, i...)
Base.IndexStyle(::Type{<:FallbackArrayWrapper}) = IndexLinear()

@show LoopVectorization.REGISTER_COUNT

@time @testset "LoopVectorization.jl" begin

    @test isempty(detect_unbound_args(LoopVectorization))

    @time include("printmethods.jl")

    @time include("fallback.jl")

    @time include("utils.jl")

    @time include("arraywrappers.jl")

    @time include("check_empty.jl")

    @time include("offsetarrays.jl")

    @time include("tensors.jl")

    @time include("map.jl")

    @time include("filter.jl")

    @time include("mapreduce.jl")

    @time include("ifelsemasks.jl")

    @time include("dot.jl")

    @time include("special.jl")

    @time include("gemv.jl")

    @time include("miscellaneous.jl")

    @time include("copy.jl")

    @time include("broadcast.jl")

    # I test  locally on master; times out on Travis.
    if isnothing(get(ENV, "TRAVIS_BRANCH", nothing)) || LoopVectorization.REGISTER_COUNT != 32 || VERSION.minor == 4
        @time include("gemm.jl")
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractDensity

Subtypes of `AbstractDensity` must implement the function

* `BAT.eval_logval_unchecked`

For likelihood densities this is typically sufficient, since shape, and
variate bounds will be inferred from the prior.

Densities with a known variate shape may also implement

* `ValueShapes.varshape`

Densities with known variate bounds may also implement

* `BAT.var_bounds`

!!! note

    The function `BAT.var_bounds` is not part of the stable public BAT-API,
    it's name and arguments may change without deprecation.
"""
abstract type AbstractDensity end
export AbstractDensity

"""
    BAT.eval_logval_unchecked(density::AbstractDensity, v::Any)

Compute log of the value of a multivariate density function for the given
variate/parameter-values.

Input:

* `density`: density function
* `v`: argument, i.e. variate / parameter-values

Note: If `eval_logval_unchecked` is called with an argument that is out of bounds,
the behavior is undefined. The result for arguments that are not within
bounds is *implicitly* `-Inf`, but it is the caller's responsibility to handle
these cases.
"""
function eval_logval_unchecked end

"""
    var_bounds(
        density::AbstractDensity
    )::Union{AbstractVarBounds,Missing}

*BAT-internal, not part of stable public API.*

Get the parameter bounds of `density`. See `eval_logval_unchecked` and
`eval_logval` for the implications and handling of bounds.
If bounds are missing, `eval_logval_unchecked` must be prepared to
handle any parameter values.
"""
var_bounds(density::AbstractDensity) = missing

"""
    ValueShapes.totalndof(density::AbstractDensity)::Union{Int,Missing}

Get the number of degrees of freedom of the variates of `density`. May return
`missing`, if the shape of the variates is not fixed.
"""
function ValueShapes.totalndof(density::AbstractDensity)
    bounds = var_bounds(density)
    ismissing(bounds) ? missing : ValueShapes.totalndof(bounds)
end

"""
    ValueShapes.varshape(
        density::AbstractDensity
    )::Union{ValueShapes.AbstractValueShape,Missing}

    ValueShapes.varshape(
        density::DistLikeDensity
    )::ValueShapes.AbstractValueShape

Get the shapes of the variates of `density`.

For prior densities, the result must not be `missing`, but may be `nothing` if
the prior only supports unshaped variate/parameter vectors.
"""
ValueShapes.varshape(density::AbstractDensity) = missing

"""
    eval_logval(
        density::AbstractDensity,
        v::Any,
        T::Type{:Real} = density_logval_type(v, density);
        use_bounds::Bool = true,
        strict::Bool = false
    )::T

*BAT-internal, not part of stable public API.*

Evaluates density log-value via `eval_logval_unchecked`.

Throws an exception on any of these conditions:

* The variate shape of `density` (if known) does not match the shape of `v`.
* The return value of `eval_logval_unchecked` is `NaN`.
* The return value of `eval_logval_unchecked` is an equivalent of positive
  infinity.

Options:

* `use_bounds`: Return an equivalent of negative infinity if `v` is out of
  bounds. `use_bounds` should only be set to `false` when it is known that
  the `v` is within bounds.

* `strict`: Throw an exception if `v` is out of bounds.
"""
function eval_logval end

function eval_logval(
    density::AbstractDensity,
    v::Any,
    T::Type{<:Real} = density_logval_type(v, density);
    use_bounds::Bool = true,
    strict::Bool = false
)
    v_shaped = reshape_variate(varshape(density), v)
    if use_bounds && !variate_is_inbounds(density, v_shaped, strict)
        return log_zero_density(T)
    end

    # ToDo: Make Zygote-compatible, by wrapping the following exception
    # augmentation mechanism in a function `get_density_logval_with_rethrow`
    # with a custom pullback:
    logval::T = try
        # ToDo: Mechanism to allow versions of eval_logval_unchecked for
        # wrapped distributions and similar that avoid stripscalar:
        eval_logval_unchecked(density, stripscalar(v_shaped))
    catch err
        rethrow(_density_eval_error(density, v, err))
    end

    _check_density_logval(density, v, logval, strict)

    return logval
end

function _density_eval_error(density::AbstractDensity, v::Any, err::Any)
    ErrorException("Density evaluation failed at v = $(variate_for_msg(v)) due to exception $err, density has type $(typeof(density))")
end

function _check_density_logval(density::AbstractDensity, v::Any, logval::Real, strict::Bool)
    if isnan(logval)
        throw(ErrorException("Log-density must not evaluate to NaN, v = $(variate_for_msg(v)) , density has type $(typeof(density))"))
    end

    if !(logval < typeof(logval)(+Inf))
        throw(ErrorException("Log-density must not evaluate to posivite infinity, v = $(variate_for_msg(v)), density has type $(typeof(density))"))
    end

    nothing
end

function variate_is_inbounds(density::AbstractDensity, v::Any, strict::Bool)
    unshaped_v = unshaped_variate(varshape(density), v)
    bounds = var_bounds(density)
    if !ismissing(bounds) && !(unshaped_v in bounds)
        if strict
            throw(ArgumentError("Parameter(s) out of bounds, density has type $(typeof(density))"))
        else
            return false
        end
    else
        true
    end
end

_strip_duals(x) = x
_strip_duals(x::AbstractVector{<:ForwardDiff.Dual}) = ForwardDiff.value.(x)

function variate_for_msg(v::Any)
    # Strip dual numbers to make errors more readable:
    shape = valshape(v)
    v_real_unshaped = _strip_duals(unshaped(v))
    v_real_shaped = shape(v_real_unshaped)
    stripscalar(v_real_shaped)
end

"""
    eval_gradlogval(density::AbstractDensity, v::AbstractVector{<:Real})
    grad
*BAT-internal, not part of stable public API.*

Compute the log of the value of a multivariate density function, as well as
the gradient of the log-value for the given variate/parameter-values.

Input:

* `density`: density function
* `v`: argument, i.e. variate / parameter-values

Returns a tuple of the log density value and it's gradient.

Note: This function should *not* be specialized for custom density types!
"""
function eval_gradlogval end

function eval_gradlogval(
    density::AbstractDensity,
    v::Any,
)
    log_f = logvalof(density)

    shape = valshape(v)
    v_unshaped = unshaped(v)

    n = length(eachindex(v_unshaped))
    P = eltype(v_unshaped)
    T = density_logval_type(v_unshaped, density)

    grad_logd_unshaped = Vector{P}(undef, n)
    chunk = ForwardDiff.Chunk(v_unshaped)
    config = ForwardDiff.GradientConfig(log_f, v_unshaped, chunk)
    result = DiffResults.MutableDiffResult(zero(T), (grad_logd_unshaped,))
    ForwardDiff.gradient!(result, log_f, v_unshaped, config)
    logd = DiffResults.value(result)

    gradshape = replace_const_shapes(ValueShapes.const_zero_shape, shape)

    grad_logd = gradshape(grad_logd_unshaped)
    (logd = logd, grad_logd = grad_logd)
end

"""
    BAT.density_logval_type(v::Any, density::AbstractDensity, T::Type{<:Real} = Float32)

*BAT-internal, not part of stable public API.*

Determine a suitable return type of log-density functions, given a variate
`v` and an optional additional default result type `T`.
"""
function density_logval_type end

@inline function density_logval_type(v::AbstractArray{<:Real}, density::AbstractDensity, T::Type{<:Real} = Float32)
    U = float(eltype(v))
    promote_type(T, U)
end

@inline function density_logval_type(v::Any, density::AbstractDensity, T::Type{<:Real} = Float32)
    density_logval_type(unshaped_variate(varshape(density), v), density, T)
end

"""
    BAT.log_zero_density(T::Type{<:Real})

log-density value to assume for regions of implicit zero density, e.g.
outside of variate/parameter bounds/support.

Returns an equivalent of negative infinity.
"""
log_zero_density(T::Type{<:Real}) = float(T)(-Inf)

"""
    BAT.is_log_zero(x::Real, T::Type{<:Real} = typeof(x)}

*BAT-internal, not part of stable public API.*

Check if x is an equivalent of log of zero, resp. negative infinity,
in respect to type `T`.
"""
function is_log_zero(x::Real, T::Type{<:Real} = typeof(x))
    U = typeof(x)

    FT = float(T)
    FU = float(U)

    x_notnan = !isnan(x)
    x_isinf = !isfinite(x)
    x_isneg = x < zero(x)
    x_notgt1 = !(x > log_zero_density(FT))
    x_notgt2 = !(x > log_zero_density(FU))
    x_iseq1 = x ≈ log_zero_density(FT)
    x_iseq2 = x ≈ log_zero_density(FU)

    x_notnan && ((x_isinf && x_isneg) | x_notgt1 | x_notgt2 | x_iseq1 | x_iseq1)
end

struct LogValOfDensity{D<:AbstractDensity} <: Function
    density::D
end

(lvd::LogValOfDensity)(v::Any) = eval_logval(lvd.density, v)

"""
    logvalof(density::AbstractDensity)::Function

Returns a function that computes the logarithmic value of `density` at a
given point:

```julia
f = logvalof(density)
log_density_at_v = f(v)
```
"""
logvalof(density::AbstractDensity) = LogValOfDensity(density)

struct LogValGradOfDensity{D<:AbstractDensity} <: Function
    density::D
end

(lvdg::LogValGradOfDensity)(v::Any) = eval_gradlogval(lvdg.density, v)

"""
    eval_gradlogval(density::AbstractDensity)::Function

Returns a function that is equivalent to

```julia
    v -> eval_gradlogval(density, v)
```
"""
eval_gradlogval(density::AbstractDensity) = LogValGradOfDensity(density)

"""
    abstract type DistLikeDensity <: AbstractDensity

A density that implements part of the `Distributions.Distribution` interface.
Such densities are suitable for use as a priors.

Typically, custom priors should be implemented as subtypes of
`Distributions.Distribution`. BAT will automatically wrap them in a subtype of
`DistLikeDensity`.

Subtypes of `DistLikeDensity` are required to support more functionality
than an [`AbstractDensity`](@ref), but less than a
`Distribution{Multivariate,Continuous}`.

A `d::Distribution{Multivariate,Continuous}` can be converted into (wrapped
in) an `DistLikeDensity` via `conv(DistLikeDensity, d)`.

The following functions must be implemented for subtypes:

* `BAT.eval_logval_unchecked`

* `ValueShapes.varshape`

* `Distributions.sampler`

* `Statistics.cov`

* `BAT.var_bounds`

!!! note

    The function `BAT.var_bounds` is not part of the stable public BAT-API
    and subject to change without deprecation.
"""
abstract type DistLikeDensity <: AbstractDensity end
export DistLikeDensity

"""
    var_bounds(density::DistLikeDensity)::AbstractVarBounds

*BAT-internal, not part of stable public API.*

Get the parameter bounds of `density`. Must not be `missing`.
"""
function var_bounds end

"""
    ValueShapes.totalndof(density::DistLikeDensity)::Int

Get the number of degrees of freedom of the variates of `density`. Must not be
`missing`, a `DistLikeDensity` must have a fixed variate shape.
"""
ValueShapes.totalndof(density::DistLikeDensity) = totalndof(var_bounds(density))

"""
    BAT.AnyDensityLike = Union{...}

Union of all types that BAT will accept as a probability density, resp. that
`convert(AbstractDenstiy, d)` supports:

* [`AbstractDensity`](@ref)
* `Distributions.Distribution`
"""
const AnyDensityLike = Union{
    AbstractDensity,
    Distributions.ContinuousDistribution
}
export AnyDensityLike

"""
    BAT.AnySampleable = Union{...}

Union of all types that BAT can sample from:

* [`AbstractDensity`](@ref)
* [`DensitySampleVector`](@ref)
* `Distributions.Distribution`
"""
const AnySampleable = Union{
    AbstractDensity,
    DensitySampleVector,
    Distributions.Distribution
}
export AnySampleable

"""
    BAT.AnyIIDSampleable = Union{...}

Union of all distribution/density-like types that BAT can draw i.i.d.
(independent and identically distributed) samples from:

* [`DistLikeDensity`](@ref)
* `Distributions.Distribution`
"""
const AnyIIDSampleable = Union{
    DistLikeDensity,
    Distributions.Distribution
}
export AnyIIDSampleable


text:
using Match
using Base.Filesystem
using CSV
using Statistics
# using JuliaDB

function download_and_extract(set::String)

    url = @match set begin
        "destillation_column" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/destill.dat.gz"
        "glass_furnace" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/glassfurnace.dat.gz"
        "power_plant" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/powerplant.dat.gz"
        "evaporator" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/evaporator.dat.gz"
        "stirring_tank" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/pHdata.dat.gz"
        "fractional_distillation_column" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/distill2.dat.gz"
        "industrial_dryer" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/dryer2.dat.gz"
        "heat_exchanger" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/exchanger.dat.gz"
        "winding_process" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/winding.dat.gz"
        "stirred_tank_reactor" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/cstr.dat.gz"
        "power_plant_2" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/steamgen.dat.gz"
        "ball_and_beam" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/ballbeam.dat.gz"
        "hair_dryer" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/dryer.dat.gz"
        "cd_arm" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/CD_player_arm.dat.gz"
        "wing_flutter" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/flutter.dat.gz"
        "flexible_arm" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/robot_arm.dat.gz"
        "steel_subframe" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/flexible_structure.dat.gz"
        "pregnant_woman" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/biomedical/foetal_ecg.dat.gz"
        "tongue_displacement" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/biomedical/tongue.dat.gz"
        "lake_erie" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/environmental/erie.dat.gz"
        "two_layer_wall" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/thermic/thermic_res_wall.dat.gz"
        "heating_system" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/thermic/heating_system.dat.gz"
        "internet_traffic" => "ftp://ftp.esat.kuleuven.be/pub/SISTA/data/timeseries/internet_traffic.dat.gz"
        "silverbox" => "http://www.nonlinearbenchmark.org/FILES/BENCHMARKS/SILVERBOX/SilverboxFiles.zip"
        "cascaded_tanks" => "http://www.nonlinearbenchmark.org/FILES/BENCHMARKS/CASCADEDTANKS/CascadedTanksFiles.zip"
        "f16" => "http://www.nonlinearbenchmark.org/FILES/BENCHMARKS/F16/F16GVT_Files.zip"
        "emps" => "http://www.nonlinearbenchmark.org/FILES/BENCHMARKS/EMPS/EMPS.zip"
        "wiener_hammerstein" => "http://www.nonlinearbenchmark.org/FILES/BENCHMARKS/WIENERHAMMERSTEINPROCESS/WienerHammersteinFiles.zip"
    end
    try
        print("Making directory ", "./data", "...")
        mkdir("./data")
        println("Done!")
    catch
        println("Directory already exists.")
    end

    # paths and file names
    set_name = split(url, "/")[end]
    folder = string("./data/", set, "/")
    file_path = string(folder, set_name)
    unzipped_file_name = string(split(set_name, ".")[1], "")
    unzipped_file_path = string(folder, unzipped_file_name)

    try
        print("Making directory ", folder, "...")
        mkdir(folder)
        println("Done!")
    catch
        println("Directory already exists.")
    end

    # Download data. If data exists, return
    print("Downloading data from ", url, "...")
    if isfile(string(unzipped_file_path, ".zip")) || isdir(unzipped_file_path) || isfile(string(unzipped_file_path, ".dat"))
        println("File already exists.")
        return unzipped_file_path
    end

    download(url, file_path);
    println("Done!")

    # Unzip based on file type
    file_type = split(set_name, ".")[end]
    cmd = @match file_type begin
        "gz" => `gzip -d $file_path`
        "zip" => `unzip $file_path -d $folder`
    end
    run(cmd)

    # remove the .gz to get the atual file name
    return file_path
end

"""
Loads the data, filters missing values and casts.
"""
function load_io_data(fp::String, dtype=Float64; header=false)

    filter_missing(row) = Iterators.filter(r -> !ismissing(r), row)
    cast(row) = map(r -> dtype(r), row)
    collate(row) = Iterators.reduce(vcat, row)
    process_row(row) = row |> filter_missing |> cast |> collate

    data = map(process_row, CSV.File(fp, header=header))
    return data
end

# data pre-processing
# returns a whitener and dewhitener
function standardize(signal)
    μ = mean(signal)
    σ = std(signal)

    process(zi) = (zi - μ) ./ σ
    return process
end

function unstandardize(signal)
    μ = mean(signal)
    σ = std(signal)

    unprocess(zi) = σ .* zi .+ μ
    return unprocess
end

function load_silverbox()
    fp = "./data/silverbox/SilverboxFiles/SNLS80mV.csv"
    data = load_io_data(fp; header=true)
    u_raw = map(d -> d[1:1, :], data)
    y_raw = map(d -> d[2:2, :], data)

    u = standardize(u_raw).(u_raw)
    y = standardize(y_raw).(y_raw)

    uv = u[1:40400]
    yv = y[1:40400]

    ut = u[40400:end]
    yt = y[40400:end]

    stats = (μu = 0.0, μy = 0.0, σu = 1.0, σy = 1.0)

    return (ut, yt), (uv, yv), stats
end

function load_cascaded_tanks()
    fp = "./data/cascaded_tanks/CascadedTanksFiles/dataBenchmark.csv"
    raw_data = load_io_data(fp; header=true)
    raw_data[1] = raw_data[1][1:4]
    # data = standardize(raw_data).(raw_data)

    ut = map(d -> d[1:1, :], raw_data)
    yt = map(d -> d[2:2, :], raw_data)

    μu = mean(ut)
    σu = std(ut)
    μy = mean(yt)
    σy = std(yt)

    uv = standardize(ut).(map(d -> d[3:3, :], raw_data))
    yv = standardize(yt).(map(d -> d[4:4, :], raw_data))

    ut = standardize(ut).(ut)
    yt = standardize(yt).(yt)

    return (ut, yt), (uv, yv), (μu = μu, μy = μy, σu = σu, σy = σy)
end

function load_wing_flutter()
    fp = "./data/wing_flutter/flutter.dat"
    data = load_io_data(fp)

    u_raw = map(d -> d[1:1, :], data)
    y_raw = map(d -> d[2:2, :], data)

    u = standardize(u_raw).(u_raw)
    y = standardize(y_raw).(y_raw)

    uv = u
    yv = y
    ut = u
    yt = y

    # Currently not bothering accounting for stats...
    stats = (μu = 0.0, μy = 0.0, σu = 1.0, σy = 1.0)
    return (ut, yt), (uv, yv), stats
end

function load_flexible_arm()
    fp = "./data/flexible_arm/robot_arm.dat"
    data = load_io_data(fp)

    u_raw = map(d -> d[1:1, :], data)
    y_raw = map(d -> d[2:2, :], data)

    u = standardize(u_raw).(u_raw)
    y = standardize(y_raw).(y_raw)

    uv = u
    yv = y
    ut = u
    yt = y

    # Currently not bothering accounting for stats...
    stats = (μu = 0.0, μy = 0.0, σu = 1.0, σy = 1.0)
    return (ut, yt), (uv, yv), stats
end

function load_liquid_saturated_heat_exchanger()
    fp  = "./data/heat_exchanger/exchanger.dat"
    data = load_io_data(fp; header=false)

    u_raw = map(d -> d[2:2, :], data)
    y_raw = map(d -> d[3:3, :], data)

    μu = mean(u_raw)
    σu = std(u_raw)
    μy = mean(y_raw)
    σy = std(y_raw)

    u = standardize(u_raw).(u_raw)
    y = standardize(y_raw).(y_raw)

    uv = u[1000:3000]
    yv = y[1000:3000]

    ut = u[1:1000]
    yt = y[1:1000]
    stats = (μu = μu, μy = μy, σu = σu, σy = σy)
    return (ut, yt), (uv, yv), stats

end

function load_f16()
    fp1 = "./data/f16/F16GVT_Files/BenchmarkData/F16Data_FullMSine_Level1.csv"
    fp2 = "./data/f16/F16GVT_Files/BenchmarkData/F16Data_FullMSine_Level2_Validation.csv"
    fp3 = "./data/f16/F16GVT_Files/BenchmarkData/F16Data_FullMSine_Level3.csv"
    fp4 = "./data/f16/F16GVT_Files/BenchmarkData/F16Data_FullMSine_Level4_Validation.csv"
    fp5 = "./data/f16/F16GVT_Files/BenchmarkData/F16Data_FullMSine_Level5.csv"
    fp6 = "./data/f16/F16GVT_Files/BenchmarkData/F16Data_FullMSine_Level6_Validation.csv"
    fp7 = "./data/f16/F16GVT_Files/BenchmarkData/F16Data_FullMSine_Level7.csv"

    split_data(xi) = (xi[1:2, :], xi[3:5, :])
    row_filt(D) = map(x -> x[1:5], D)
    stand = standardize(row_filt(load_io_data(fp4, header=true)))  # use this dset for standarization

    loader(fp) = split_data.(stand.(row_filt(load_io_data(fp, header=true))))

    u1, y1 = unzip(loader(fp1))
    u2, y2 = unzip(loader(fp2))
    u3, y3 = unzip(loader(fp3))
    u4, y4 = unzip(loader(fp4))
    u5, y5 = unzip(loader(fp5))
    u6, y6 = unzip(loader(fp6))
    u7, y7 = unzip(loader(fp7))

    ut = hcat.(u1, u3, u5, u7)
    yt = hcat.(y1, y3, y5, y7)
    uv = hcat.(u2, u4, u6)
    yv = hcat.(y2, y4, y6)

    stats = (μu = 0.0, μy = 0.0, σu = 1.0, σy = 1.0)
    return (ut, yt), (uv, yv), stats
end

#
# Identification is model (r, u) -> y
#
function load_wiener_hammerstein()
    fp1 = "./data/wiener_hammerstein/WienerHammersteinFiles/WH_MultisineFadeOut.csv"
    fp2 = "./data/wiener_hammerstein/WienerHammersteinFiles/WH_TestDataset.csv"
    # fp3 = "./data/wiener_hammerstein/WienerHammersteinFiles/WH_SineSweepInput_meas.csv"

    row_filt(D) = map(x -> x[1:6], D)
    raw_train_data = row_filt(load_io_data(fp1, header=true))

    # Creater a standardizer for the inputs and outputs
    # D = [[di[1], di[2], di[5]] for di in raw_train_data]
    # stand = standardize(raw_train_data)

    # split into inputs and outputs, for two experiments
    split_data(data) = map(d -> ([d[1] d[2]; d[3] d[4]], [d[5] d[6]]), data)
    ut, yt = unzip((split_data(raw_train_data)))

    σu = std(ut)[:, 1][:, :]
    σy = std(yt)[:, 1][:, :]
    μu = mean(ut)[:, 1][:, :]
    μy = mean(yt)[:, 1][:, :]

    S(di) = ((di[1].-μu) ./ σu, (di[2].-μy) ./ σy)

    test_data = row_filt(load_io_data(fp2, header=true))

    ut, yt = unzip(S.(split_data(raw_train_data)))
    uv, yv = unzip(S.(split_data(test_data)))

    # statistics = (σy = σ[5:6], σu = σ[1:4], μu = μ[1:4], μy=μ[5:6])
    statistics = (σy = σy, σu = σu, μu = μu, μy = μy)
    return (ut, yt), (uv, yv), statistics
end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
module TestSetfield
include("test_core.jl")
include("test_settable.jl")
include("test_staticarrays.jl")
include("test_kwonly.jl")
include("test_quicktypes.jl")
include("perf.jl")
end  # module


text:
@testset "LaTeXEngine" begin
    @test begin
        texengine!(PDFLaTeX)
        texengine() == PDFLaTeX()
    end

    @test LatexSVG.runlatex(XeLaTeX(), "test.tex", raw"C:\path to\out"; extra_args=["--shell-escape", "-synctex=1"], runcommand=false) == `xelatex --shell-escape -synctex=1 -no-pdf -output-directory="C:\\path to\\out" -quiet -halt-on-error "test.tex"`
    @test LatexSVG.runlatex(PDFLaTeX(), "test.tex", raw"C:\path to\out"; extra_args=["--shell-escape", "-synctex=1"], runcommand=false) == `pdflatex --shell-escape -synctex=1 -output-format=dvi -output-directory="C:\\path to\\out" -quiet -halt-on-error "test.tex"`
    @test LatexSVG.runlatex(LuaLaTeX(), "test.tex", raw"C:\path to\out"; extra_args=["--shell-escape", "--synctex=1"], runcommand=false) == `lualatex --shell-escape --synctex=1 --output-format=dvi --output-directory="C:\\path to\\out" --interaction=batchmode --halt-on-error "test.tex"`

    @test_throws ErrorException begin
        struct _MyEngine <: LatexSVG.LaTeXEngine end
        LatexSVG.runlatex(_MyEngine(), "test.tex", raw"C:\path to\out"; extra_args=["--shell-escape", "--synctex=1"])
    end
end


text:
# # Diagrams
#
# Catlab supports the expression of equations in an arbitrary category via
# diagrams, or functors from a finitely presented category. This view of
# diagrammatic equations builds on the Catlab approach to functorial semantics
# by asserting that a system of equations of shape J is a functor F: J → 𝐂. We
# can draw these equations with Graphviz, where each node is an object in
# X:Ob(J) labeled with its type F(X) and each edge is a morphism in J labeled
# with its image under F.

using Catlab, Catlab.Theories, Catlab.CategoricalAlgebra
using Catlab.Graphs, Catlab.Graphics, Catlab.Programs
using Catlab.Graphs.BasicGraphs: TheoryGraph
using Catlab.Graphics.Graphviz

draw(D::FinFunctor) = to_graphviz(D, node_labels=true, edge_labels=true, prog="neato")

# ## Present a diagram in a given category
#
# Recall that a *diagram* in a category ``C`` is a functor ``F: J → C`` from a
# small category ``J`` into ``C``. Given the category ``C``, this macro presents a
# diagram in ``C``, i.e., constructs a finitely presented indexing category ``J``
# together with a functor ``F: J → C``. This method of simultaneous definition is
# often more convenient than defining ``J`` and ``F`` separately.
# often more convenient than defining ``J`` and ``F`` separately, as could be
# accomplished by calling `@fincat` and then `@finfunctor`.
# For example, the limit of the following diagram consists of the paths of length
# two in a graph:

D₂ = @free_diagram TheoryGraph begin
  v::V
  (e₁, e₂)::E
  tgt(e₁) == v
  src(e₂) == v
end

# Morphisms in the indexing category can be left unnamed, which is convenient for
# defining free diagrams (see also `@free_diagram`. Of course, unnamed morphisms cannot
# be referenced by name within the `@diagram` call or in other settings, which
# can sometimes be problematic.
#
# We can also draw the resulting diagrams.

draw(D₂)

# For small equations the point-free notation commonly employed in functional programming is very convenient; however, there is a reason it is not the standard approach to presenting equations in mathematical writing. Variables are just too useful! As the size of the system of equations grows, it becomes more and more convenient to use variable names. This is why `Catlab.Programs.@program` exists to help people write SMC morphisms with the point-ful notation they are familiar with in imperative or procedural programming.

# We can describe a triangle in a graph using the vertex variables v₁, v₂, v₃ and edge variables e₁, e₂, e₃. Then we use the equation notation to assert the `src` and `tgt` relationships between the edges and vertices.
D₃ = @free_diagram TheoryGraph begin
  (v₁, v₂, v₃)::V
  (e₁, e₂, e₃)::E
  src(e₁) == v₁
  tgt(e₁) == v₂
  src(e₂) == v₂
  src(e₃) == v₁
  tgt(e₂) == v₃
  tgt(e₃) == v₃
end

# You can see the shape of a triangle when you draw this diagram. This coincidence can be systematically understood (by an experienced category theorist) via the Grothendieck construction and representable functors.
draw(D₃)

# ## Sequences in
# In any category with an endomorphism f: A → A, we can think of recurrence equations as aₙ = f(aₙ₋₁)) starting at a₀.
@present 𝗖(FreeCartesianCategory) begin
  (A,B)::Ob
  f::Hom(A,A)
  g::Hom(A,B)
end

seq₃ = @free_diagram 𝗖 begin
  (a₀,a₁,a₂,a₃)::A
  a₁ == f(a₀)
  a₂ == f(a₁)
  a₃ == f(a₂)
end

draw(seq₃)

# Inspired by linear recurrence relations like the Fibonacci sequence we can think of A and B as vector spaces. An endomorphism f: A → A and morphism g: A → B defines a linear recurrence by applying f to update the system state and g to compute the current term in the sequence. This can be visualized as follows.

obs_seq₃ = @free_diagram 𝗖 begin
  (a₀,a₁,a₂,a₃)::A
  (b₁, b₂, b₃ )::B
  a₁ == f(a₀)
  b₁ == g(a₁)
  a₂ == f(a₁)
  b₂ == g(a₂)
  a₃ == f(a₂)
  b₃ == g(a₃)
end
draw(obs_seq₃)

# In the case of the Fibonacci sequence A is ℝ² and B is ℝ with f = [1 1; 1 0] and g = π₁.

@present 𝐃(FreeCartesianCategory) begin
  (ℝ²,ℝ)::Ob
  f::Hom(ℝ²,ℝ²)
  π₁::Hom(ℝ²,ℝ)
end

fib_seq₃ = @free_diagram 𝐃 begin
  (a₀,a₁,a₂,a₃)::ℝ²
  (b₁, b₂, b₃ )::ℝ
  a₁ == f(a₀)
  b₁ == π₁(a₁)
  a₂ == f(a₁)
  b₂ == π₁(a₂)
  a₃ == f(a₂)
  b₃ == π₁(a₃)
end

draw(fib_seq₃)

# A lift of this diagram is a choice of vectors `a₀,…,a₃` and scalars `b₁,…,b₃` that satisfy the linear recurrence formulation of the Fibonacci sequence. The arrows in 𝐃 are treated as fixed for the purpose of the lifting problem. Suppose you knew `b₁,…,b₃`, and wanted to infer the state vectors `a₁,…,a₃`. This can be seen as a boundary value problem where you have observations on a subspace playing the role of boundary conditions and want to infer the system state on the interior of the domain. As you can see from this diagram, solving the lifting problem with only partial information can require some creativity. You need to use the structure of the codomain category, in this case 𝐃 which is a subcategory of finite dimensional vector spaces over ℝ and the combinatorics of 𝐉 to set up some equations you can solve. In order to obtain an efficient algorithm for the lifting problem, you would need to exploit the structure of the morphisms `f,π₁` to create a system of equations.

# An observation that `[π₁f;π₁] = f` lets you set up the equation `[b₃ ; b₂] = f a₂` and then you can use the fact that `f` is invertable to solve for a unique `a₂`. Once you have a single value of `aᵢ` you can get the rest from the application of `f` or `f⁻¹`. A general purpose algorithm for solving lifting problems in arbitrary categories would be a wildly powerful tool.

# ## Newton's Method
# The equations that we have seen aren't particularly interesting, so we turn to a classic of numerical methods. Newton's method for root finding. For an overview of Netwon's method see [Fundamentals of Numerical Computation](https://fncbook.github.io/fnc/nonlineqn/newton.html). The following presentation doesn't know that f′ is the derivative of f, they are just two functions that are evocatively named. We could use [CombinatorialSpaces.jl](https://github.com/AlgebraicJulia/CombinatorialSpaces.jl) to formulate this in a richer categorical setting where we could assert `f′ := d/dx(f)`.
@present Analytic(FreeCartesianCategory) begin
  (ℝ,ℝ²)::Ob
  π₁::Hom(ℝ², ℝ)
  π₂::Hom(ℝ², ℝ)
  plus ::Hom(ℝ², ℝ)
  times::Hom(ℝ², ℝ)
  f    ::Hom(ℝ,ℝ)
  f′   ::Hom(ℝ,ℝ)
end
# According to the standard formula xₖ₊₁ = xₖ - f(xₖ)/f′(xₖ). The standard presentation of Newton's method relies on the fact that ℝ is a field to use division in the definition of the iterative procedure. Because of the constraint that you can't divide by 0 in a field, fields are not models of any algebraic theory. Because of this, we can multiply both sides by f′(xₖ) and define a Newton's method iteration without reference to division. We also can avoid negation by adding the f(xₖ) term on both sides.

# Once we have put Newton's method into this relational form we can use a prefix notation for operators which makes the structure closer to the diagrammatic approach `plus(times(f′(xₖ), xₖ₊₁), f(xₖ)) == times(f′(xₖ), xₖ)`.
newtons = @free_diagram Analytic begin
  (xₖ, xₖ₊₁, dₖ, fx, v, ∏)::ℝ
  (p₁, p₂, p₃)::ℝ²
  dₖ  == f′(xₖ)
  π₁(p₁) == dₖ
  π₂(p₁) == xₖ₊₁
  ∏ == times(p₁)
  fx == f(xₖ)
  π₁(p₂) == ∏
  π₂(p₂) == fx
  plus(p₂) == v
  π₁(p₃) == dₖ
  π₂(p₃) == xₖ
  times(p₃) == v
end

draw(newtons)

# This is definitely a case where standard mathematical notation wins for brevity and clarity. To compete with traditional notation, we need a few ergonamic improvements to the diagrammatic equation approach.
# 1. The ability to assert equations without introducing temporary variables,
# 2. The ability to represent morphisms f: A×B → C as bivariate functions like `f(a:A, b:B)`,
# While encoding simple equations can be more verbose than the traditional notation, the diagrammatic encoding provides a direct route to creating a category whose objects are systems of equations and whose morphisms are relationships between systems of equations. That category is the first step towards leveraging the constructive approach in Catlab to making hierarchical representations of complex systems of equations.


text:
# not exported
abstract type ImageEncoder end
struct BigBlocks   <: ImageEncoder end
struct SmallBlocks <: ImageEncoder end

const alpha_chars = ('⋅', '░', '▒', '▓', '█')
function _charof(alpha)
    idx = round(Int, alpha * (length(alpha_chars)-1))
    alpha_chars[clamp(idx + 1, 1, length(alpha_chars))]
end

"""
    encodeimg(enc::ImageEncoder, colordepth::TermColorDepth, img, [maxheight], [maxwidth])

Transforms the pixel of the given image `img`, which has to be an
array of `Colorant`, into a string of unicode characters using
ansi terminal colors.

- The encoder `enc` specifies which kind of unicode represenation
  should be used.

- The `colordepth` can either be `TermColor256()` or `TermColor24bit()`
  and specifies which terminal color codes should be used.

- `maxheight` and `maxwidth` specify the maximum numbers of
  string characters that should be used for the resulting image.
  Larger images are downscaled automatically using `restrict`.

The function returns a tuple with three elements:

1. A vector of strings containing the encoded image.
   Each element represent one line. The lines do not contain
   newline characters.

2. Number of lines in the vector.

3. Number of visible characters per line (the remaining are colorcodes).
"""
function encodeimg(
        ::SmallBlocks,
        colordepth::TermColorDepth,
        img::AbstractMatrix{<:Colorant},
        maxheight::Int = 50,
        maxwidth::Int = 80)
    maxheight = max(maxheight, 5)
    maxwidth  = max(maxwidth,  5)
    h, w = map(length, axes(img))
    while ceil(h/2) > maxheight || w > maxwidth
        img = restrict(img)
        h, w = map(length, axes(img))
    end
    yinds, xinds = axes(img)
    io = IOBuffer()
    for y in first(yinds):2:last(yinds)
        print(io, Crayon(reset = true))
        for x in xinds
            fgcol = _colorant2ansi(img[y,x], colordepth)
            bgcol = if y+1 <= last(yinds)
                _colorant2ansi(img[y+1,x], colordepth)
            else
                # if reached it means that the last character row
                # has only the upper pixel defined.
                Symbol("nothing")
            end
            print(io, Crayon(foreground=fgcol, background=bgcol), "▀")
        end
        println(io, Crayon(reset = true))
    end
    replace.(readlines(seek(io,0)), Ref("\n" => ""))::Vector{String}, length(1:2:h), w
end

function encodeimg(
        ::BigBlocks,
        colordepth::TermColorDepth,
        img::AbstractMatrix{<:Colorant},
        maxheight::Int = 50,
        maxwidth::Int = 80)
    maxheight = max(maxheight, 5)
    maxwidth  = max(maxwidth,  5)
    h, w = map(length, axes(img))
    while h > maxheight || 2w > maxwidth
        img = restrict(img)
        h, w = map(length, axes(img))
    end
    yinds, xinds = axes(img)
    io = IOBuffer()
    for y in yinds
        print(io, Crayon(reset = true))
        for x in xinds
            color = img[y,x]
            fgcol = _colorant2ansi(color, colordepth)
            chr = _charof(alpha(color))
            print(io, Crayon(foreground = fgcol), chr, chr)
        end
        println(io, Crayon(reset = true))
    end
    replace.(readlines(seek(io,0)), Ref("\n" => ""))::Vector{String}, h, 2w
end

# colorant vector
function encodeimg(
        enc::SmallBlocks,
        colordepth::TermColorDepth,
        img::AbstractVector{<:Colorant},
        maxwidth::Int = 80)
    maxwidth  = max(maxwidth, 5)
    w = length(axes(img, 1))
    if w > maxwidth
        img = imresize(img, maxwidth)
        w = length(axes(img, 1))
    end
    io = IOBuffer()
    print(io, Crayon(reset = true))
    for i in axes(img, 1)
        color = img[i]
        fgcol = _colorant2ansi(color, colordepth)
        chr = _charof(alpha(color))
        print(io, Crayon(foreground = fgcol), chr)
    end
    println(io, Crayon(reset = true))
    replace.(readlines(seek(io,0)), Ref("\n" => ""))::Vector{String}, 1, w
end

function encodeimg(
        enc::BigBlocks,
        colordepth::TermColorDepth,
        img::AbstractVector{<:Colorant},
        maxwidth::Int = 80)
    maxwidth  = max(maxwidth, 5)
    inds = axes(img, 1)
    w = length(inds)
    n = 3w > maxwidth ? floor(Int,maxwidth/6) : w
    io = IOBuffer()
    print(io, Crayon(reset = true))
    for i in (0:n-1) .+ first(inds)
        color = img[i]
        fgcol = _colorant2ansi(color, colordepth)
        chr = _charof(alpha(color))
        print(io, Crayon(foreground = fgcol), chr, chr, " ")
    end
    if n < w
        print(io, Crayon(reset = true), " … ")
        for i in last(inds)-n+1:last(inds)
            color = img[i]
            fgcol = _colorant2ansi(color, colordepth)
            chr = _charof(alpha(color))
            print(io, Crayon(foreground = fgcol), chr, chr, " ")
        end
    end
    println(io, Crayon(reset = true))
    replace.(readlines(seek(io,0)), Ref("\n" => ""))::Vector{String}, 1, n < w ? 3*(length(1:n) + 1 + length(w-n+1:w)) : 3w
end


text:
using MLStyle
using YaoLocations
using YaoCompiler
using YaoCompiler.Intrinsics
using OpenQASM
using YaoTargetQASM
using YaoTargetQASM.QASM
using Test

@device function test_basic(theta, phi)
    # syntax sugar
    1 => X
    2 => Z
    4 => Rx(sin(theta)+2)
    4 => Ry(phi)
    @ctrl 1 4 => X
    a = @measure 3
    # direct intrinsic
    apply(Y, 3)
    apply(X, 1, 4)
    c = measure((2, 3))
    if c == 1
        1 => X
    end
    return (a = a, b = c)
end

@device function test_pure_quantum()
    ret = @gate 1:4 => test_basic(1.0, π)
    @barrier 1:4
    @ctrl 2 1 => X
    return ret
end

test_qasm = code_qasm(Operation{typeof(test_pure_quantum), Tuple{}}; include_qelib1=true)

s = """
OPENQASM 2.0;
include "qelib1.inc";
creg a[1];
creg b[2];
qreg qreg_2[2];
qreg qreg_1[2];
x qreg_1[0];
z qreg_2[0];
rx(2.8414709848078967) qreg_1[1];
ry(pi) qreg_1[1];
CX qreg_1[0], qreg_1[1];
measure qreg_2[1] -> a[0];
y qreg_2[1];
CX qreg_1[1], qreg_1[0];
measure qreg_2[0] -> b[0];
measure qreg_2[1] -> b[1];
if (b == 1) x qreg_1[0];
barrier qreg_2, qreg_1;
CX qreg_2[0], qreg_1[0];
"""

target_qasm = OpenQASM.parse(s)
@test test_qasm ≈ target_qasm

@device function test_gate(theta, phi)
    1 => X
    2 => Z
    4 => Rx(-sin(theta)+2)
    4 => Rx(sin(theta)+tan(phi))
    4 => Ry(cos(phi)-log(phi))
    4 => Ry(QASM.cos(phi)*QASM.sqrt(phi))
    @ctrl 1 4 => X
    return
end

ft = Operation{typeof(test_gate), Tuple{Float64, Float64}}
test_qasm = code_qasm(ft; toplevel=false)

s = """
gate test_gate(theta, phi) qreg_2, qreg_3, qreg_1 {
  x qreg_1;
  z qreg_2;
  rx(-(sin(theta))+2.0) qreg_3;
  rx(sin(theta)+tan(phi)) qreg_3;
  ry(cos(phi)-ln(phi)) qreg_3;
  ry(cos(phi)*sqrt(phi)) qreg_3;
  CX qreg_1, qreg_3;
}
"""
target_qasm = OpenQASM.parse_gate(s)
@test test_qasm ≈ target_qasm

test_qasm = code_qasm(ft; toplevel=false, inline_intrinsics=true)

s = """
gate test_gate(theta, phi) qreg_2, qreg_3, qreg_1 {
  U(pi, 0, pi) qreg_1;
  U(0, 0, pi) qreg_2;
  U(-(sin(theta))+2.0, -(pi/2), pi/2) qreg_3;
  U(sin(theta)+tan(phi), -(pi/2), pi/2) qreg_3;
  U(0, cos(phi)-ln(phi), 0) qreg_3;
  U(0, cos(phi)*sqrt(phi), 0) qreg_3;
  CX qreg_1, qreg_3;
}
"""
target_qasm = OpenQASM.parse_gate(s)
@test test_qasm ≈ target_qasm


text:
using Test
using QuantumOptics

@testset "schroedinger" begin

N = 3
Ncutoff = 2
T = [0.:0.1:1.;]

Ω = [0 2 3;
     2 0 1;
     3 1 0]

ω = [1., 1.2, 1.5]

basis_fock = FockBasis(Ncutoff)
basis = tensor([basis_fock for i=1:N]...)

a = destroy(basis_fock)
at = create(basis_fock)
I = identityoperator(basis_fock)

psi0 = tensor([coherentstate(basis_fock, i%Ncutoff) for i=1:N]...)

# Interaction picture
Hrot = SparseOpType[]
for i=1:N, j=1:N
    if i==j
        continue
    end
    h = embed(basis, [i,j], [a, Ω[i,j]*at])
    push!(Hrot, h)
end
Hrot = sum(Hrot)

# Schroedinger picture
function f(t, psi)
    H = SparseOpType[embed(basis, i, ω[i]*at*a) for i=1:N]
    for i=1:N, j=1:N
        if i==j
            continue
        end
        h = embed(basis, [i,j], [a, exp(1im*(ω[i]-ω[j])*t)*Ω[i,j]*at])
        push!(H, h)
    end
    sum(H)
end

tout, psi_rot_t = timeevolution.schroedinger(T, psi0, Hrot)
tout, psi_t = timeevolution.schroedinger_dynamic(T, psi0, f)

n_op = dense(at*a)
for (i, t) in enumerate(tout)
    R = prod([embed(basis, i, exp(1im*ω[i]*t*n_op)) for i=1:N])
    psi_rot = psi_rot_t[i]
    psi = psi_t[i]
    # @test abs(dagger(psi_rot)*R*psi) < 1e-5
    rho = dm(psi)
    rho_rot = dm(psi_rot)
    @test tracedistance(rho_rot, dense(R)*rho*dagger(dense(R))) < 1e-5
end

function fout(t, psi)
 deepcopy(psi)
end
t_fout, psi_fout = timeevolution.schroedinger(T, psi0, Hrot; fout=fout)
@test t_fout == tout && psi_fout == psi_rot_t

# test integration of propagator using 2 level system
basis = SpinBasis(1//2)
su = spinup(basis)
u0 = identityoperator(basis)
sx = sigmax(basis)
sz = sigmaz(basis)

# for the time dependent equation
f(t, psi) = sx * π
tspan = 0:1.0
t, u = timeevolution.schroedinger(tspan, u0, π * sx)

# I think the tolerance on the differential equation is 1e-6, we expect the operator to be essentially the identity
@test abs(expect(sz, u[end] * su)) - abs(expect(sz, u0 * su)) < 1e-6

t, u = timeevolution.schroedinger_dynamic(tspan, u0, f)
@test abs(expect(sz, u[end] * su)) - abs(expect(sz, u0 * su)) < 1e-6

end # testset


text:
function substitution(code::Union{Expr,Symbol},p::Program_;op=nothing,arg=false,pre=nothing,post=nothing,sym=nothing, #Type of symbol between operator, pre/post script and arguments
    nargs = nothing, addArgs = nothing,                       #Args for the functions that depend on arguments
    only=["Local","Global","Identity","GlobalArray","Medium"],           #Subset of parameters to modify
    addName=nothing, expl=nothing,                            #Add something to the vectorized name of fully substitute by something
    ic=[:ic1_,:ic2_,:ic3_], index=nothing, factor=1,          #Change loop symbols and additional indexes
    dims=nothing,                                             #Operators modifying the behaviour of cells
    opF = nothing                                             #Form of the new index
    )

    if dims === nothing
        dims = p.agent.dims
    end

    if op !== nothing

        if nargs === nothing
            code = postwalk(x->@capture(x, opAux_(argAux_)) && opAux == op ?
                substitution(argAux, p;arg = true, nargs=nargs, only=only, addName=addName, expl=expl, ic=ic, index=index, factor=factor, dims=dims, opF=opF)
            : x, code)
        else
            code = postwalk(x->@capture(x, opAux_(argAux__)) && opAux == op ?
            if length(argAux) != nargs
                error("Operator ", op, " should be declared with just ", nargs, " arguments. ", length(argAux), " where passed.")
            else
                if opF !== nothing
                    opAux = opF
                end
                if addArgs !== nothing
                    :($opAux($(argAux...),$(addArgs...)))
                else
                    :($opAux($(argAux...)))
                end
            end
            : x, code)
        end

    elseif arg

        add = substitution(code, p, only=only, addName=addName, expl=expl, ic=ic, index=index[1], dims=dims)
        if factor[1] != 1
            add = :($(factor[i])*$add)
        end
        for i in 2:length(index)
            add2 = substitution(code, p, only=only, addName=addName, expl=expl, ic=ic, index=index[i], dims=dims)
            if factor[i] != 1
                add = :($add+$(factor[i])*$add2)
            else
                add = :($add+$add2)
            end
        end

        if opF !== nothing
            code = postwalk(x->@capture(x, vAux_) && vAux == :X ? add : x, opF)
        else
            code = add
        end

    else

        if "Local" in only
            for (i,v) in enumerate(p.agent.declaredSymbols["Local"])

                if pre !== nothing
                    v = Meta.parse(string(pre,v))
                elseif post !== nothing
                    v = Meta.parse(string(v,post))
                end

                if expl !== nothing
                    name = expl
                elseif addName !== nothing
                    name = Meta.parse(string(:localV,addName))
                else
                    name = :localV
                end

                if index !== nothing
                    if index[1] != 0
                        ic_ = :($(ic[1])+$(index[1]))
                    else
                        ic_ = :($(ic[1]))
                    end
                else
                    ic_ = ic[1]
                end

                if opF !== nothing
                    subs = postwalk(x->@capture(x, vAux_) && vAux == X ? :($name[$ic_,$i]) : x, opF)
                else
                    subs = :($name[$ic_,$i])
                end

                code = postwalk(x->@capture(x, vAux_) && vAux == v ? subs : x, code)

            end
        end

        if "Identity" in only
            for (i,v) in enumerate(p.agent.declaredSymbols["Identity"])

                if pre !== nothing
                    v = Meta.parse(string(pre,v))
                elseif post !== nothing
                    v = Meta.parse(string(v,post))
                end

                if expl !== nothing
                    name = expl
                elseif addName !== nothing
                    name = Meta.parse(string(:identityV,addName))
                else
                    name = :identityV
                end

                if index !== nothing
                    if index[1] != 0
                        ic_ = :($(ic[1])+$(index[1]))
                    else
                        ic_ = :($(ic[1]))
                    end
                else
                    ic_ = ic[1]
                end

                if opF !== nothing
                    subs = postwalk(x->@capture(x, vAux_) && vAux == X ? :($name[$ic_,$i]) : x, opF)
                else
                    subs = :($name[$ic_,$i])
                end

                code = postwalk(x->@capture(x, vAux_) && vAux == v ? subs : x, code)
            end
        end

        if "Global" in only
            for (i,v) in enumerate(p.agent.declaredSymbols["Global"])

                if pre !== nothing
                    v = Meta.parse(string(pre,v))
                elseif post !== nothing
                    v = Meta.parse(string(v,post))
                end

                if expl !== nothing
                    name = expl
                elseif addName !== nothing
                    name = Meta.parse(string(:globalV,addName))
                else
                    name = :globalV
                end

                if opF !== nothing
                    subs = postwalk(x->@capture(x, vAux_) && vAux == X ? :($name[$i]) : x, opF)
                else
                    subs = :($name[$i])
                end

                code = postwalk(x->@capture(x, vAux_) && vAux == v ? subs : x, code)
            end
        end

        if "GlobalArray" in only
            for (i,v) in enumerate(p.agent.declaredSymbols["GlobalArray"])

                if pre !== nothing
                    v2 = Meta.parse(string(pre,v))
                elseif post !== nothing
                    v2 = Meta.parse(string(v,post))
                else
                    v2 = v
                end

                if expl !== nothing
                    name = expl
                elseif addName !== nothing
                    name = Meta.parse(string(v,addName))
                else
                    name = v
                end

                code = postwalk(x->@capture(x, vAux_) && vAux == v2 ? name : x, code)
            end
        end

        if "Medium" in only
            for (i,v) in enumerate(p.agent.declaredSymbols["Medium"])

                if pre !== nothing
                    v = Meta.parse(string(pre,v))
                elseif post !== nothing
                    v = Meta.parse(string(v,post))
                end

                if expl !== nothing
                    name = expl
                elseif addName !== nothing
                    name = Meta.parse(string(:mediumV,addName))
                else
                    name = :mediumV
                end

                if index !== nothing
                    if index[1] != 0
                        ic1_ = :($(ic[1])+$(index[1]))
                    else
                        ic1_ = :($(ic[1]))
                    end
                    if index[2] != 0
                        ic2_ = :($(ic[2])+$(index[2]))
                    else
                        ic2_ = :($(ic[2]))
                    end
                    if index[3] != 0
                        ic3_ = :($(ic[3])+$(index[3]))
                    else
                        ic3_ = :($(ic[3]))
                    end
                else
                    ic1_ = ic[1]
                    ic2_ = ic[2]
                    ic3_ = ic[3]
                end

                if opF !== nothing
                    if dims == 1
                        subs = postwalk(x->@capture(x, vAux_) && vAux == X ? :($name[$ic1_,$i]) : x, opF)
                    elseif dims == 2
                        subs = postwalk(x->@capture(x, vAux_) && vAux == X ? :($name[$ic1_,$ic2_,$i]) : x, opF)
                    elseif dims == 3
                        subs = postwalk(x->@capture(x, vAux_) && vAux == X ? :($name[$ic1_,$ic2_,$ic3_,$i]) : x, opF)
                    end
                else
                    if dims == 1
                        subs = :($name[$ic1_,$i])
                    elseif dims == 2
                        subs = :($name[$ic1_,$ic2_,$i])
                    elseif dims == 3
                        subs = :($name[$ic1_,$ic2_,$ic3_,$i])
                    end
                end

                code = postwalk(x->@capture(x, vAux_)  && vAux == v  ? subs : x, code)

            end
        end

    end

    return code
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
"""
	flowpayoffs!(u0,u1,m::DDCbasic)

Updates flowpayoffs in matrices `u0` and `u1` based on parameters in model `m`.
"""
function flowpayoffs!(u0,u1,m::DDCbasic)
	# Fill in u0
	u0[:,2] .= -m.δ[1]
	# Fill in u1
	u1 .= m.β[1] .+ m.β[2] .* m.supportX
	u1[:,1] .-= m.δ[2]
	return nothing
end

# NOTES:
# - functions with "!" (e.g. f!(a)) mutate the input argument. The "!"
#	is not required, but it's convention to use it for functions that
# 	mutate its inputs.
# - Using mutating functions in inner functions that are run many times
#	avoids allocations, which can lead to efficiency gains.
# 	(see https://docs.julialang.org/en/v1/manual/performance-tips/)

"""
	bellman!(U0,U1,u0,u1,m::DDCbasic)

Bellman update, i.e. update expected payoffs in `U0` and `U1` given flowpayoffs `u0` and `u1`, discount factor `m.ρ` and transition matrix `m.Π`.
"""
function bellman!(U0,U1,u0,u1,m::DDCbasic)
	r0 = log.(exp.(U0[:,1])+exp.(U1[:,1]))
	U0 .= u0 .+ m.ρ .* m.Π * repeat(r0,1,2)

	r1 = log.(exp.(U0[:,2])+exp.(U1[:,2]))
	U1 .= u1 .+ m.ρ .* m.Π * repeat(r1,1,2)
	return nothing
end

"""
	solve_fixedpoint!(U0,U1,u0,u1,m::DDCbasic)

Solve fixed point and update expected payoffs in `U0` and `U1`, given flowpayoffs `u0` and `u1`, and parameters in model `m`.
"""
function solve_fixedpoint!(U0,U1,u0,u1,m::DDCbasic)
	tol = m.options.tol_fixed_point
	_U0 = U0 .+ 2 .* tol
	_U1 = U1 .+ 2 .* tol
	while maximum(abs.(_U0 .- U0)) > tol || maximum(abs.(_U1 .- U1)) > tol
		_U0 .= U0
		_U1 .= U1
		bellman!(U0,U1,u0,u1,m::DDCbasic)
	end
	return nothing
end

"""
	init_flowpayoffs(m::DDCbasic)

Initialize flowpayoff matrices `u0` and `u1` based on model `m`.
"""
function init_flowpayoffs(m::DDCbasic)
	n = length(m.supportX)
	u0 = zeros(eltype(m.β),n,2)
	u1 = zeros(eltype(m.β),n,2)
	flowpayoffs!(u0,u1,m)
	return u0,u1
end

"""
	calc_ΔU(m::DDCbasic)

Calculate and return matrix `ΔU` for model `m`.
"""
function calc_ΔU(m::DDCbasic)

	u0,u1 = init_flowpayoffs(m)

	U0 = zeros(eltype(u0),length(m.supportX),2)
	U1 = copy(U0)
	solve_fixedpoint!(U0,U1,u0,u1,m)
	ΔU = U1 .- U0
	return ΔU
end

"""
	calc_pinf(Π)

Calculate and return P∞ from transition matrix `Π`.
"""
function calc_pinf(Π)
	nk = size(Π,1)
	A = vcat(Matrix(I,nk,nk)[1:nk-1,:] .- Π[:,1:nk-1]',ones(1,nk))
	return A \ vcat(zeros(nk-1),1)
end

"""
	loglikelihood(m::DDCbasic,X,choices)

Calculates and returns loglikelihood for model `m`, given observed payoff states `X` and choices `choices`.
"""
function loglikelihood(m::DDCbasic,X,choices)
	nperiod,nfirm = size(choices)
	ΔU = calc_ΔU(m::DDCbasic)
	pexit = 1 ./ (1 .+ exp.(ΔU))
	ll = zeros(eltype(ΔU),Threads.nthreads())
	Threads.@threads for threadid in 1:Threads.nthreads()
		for f in getrange(nfirm)
			# Period 1
			ll[threadid] += log(choices[1,f] + (1-2*choices[1,f]) *
								pexit[getindX(m.supportX,X[1,f]),1])
			# Loop through other periods
			for t in 2:nperiod
				ll[threadid] += log(choices[t,f] + (1-2*choices[t,f]) *
									pexit[getindX(m.supportX,X[t,f]),
											choices[t-1,f]+1])
			end
		end
	end
	return sum(ll)
end

"""
	gen_data(m::DDCbasic,nperiod,nfirm)

Generates and returns data on payoff states and firms choices `X,choices` for model `m`, `nperiod` periods and `nfirm` firms.
"""
function gen_data(m::DDCbasic,nperiod,nfirm)
	P∞ = calc_pinf(m.Π)
	dε = GeneralizedExtremeValue(0,1,0)
	ΔU = calc_ΔU(m::DDCbasic)

	X = zeros(Int,nperiod,nfirm)

	Δε = zeros(nperiod,nfirm)
	choices = zeros(Bool,nperiod,nfirm)

	dX0 = DiscreteNonParametric(m.supportX,P∞)
	dX = [DiscreteNonParametric(m.supportX,m.Π[i,:])
			for i in eachindex(m.supportX)]

	Threads.@threads for threadid in 1:Threads.nthreads()
		Random.seed!(m.options.seed+threadid)
		for f in getrange(size(X,2))
			# Period 1
			X[1,f] = rand(dX0)
			Δε[1,f] = rand(dε) - rand(dε)
			choices[1,f] = ΔU[getindX(m.supportX,X[1,f]),1] >
								Δε[1,f]
			# Loop through other periods
			for t in 2:size(X,1)
				X[t,f] = rand(dX[getindX(m.supportX,X[t-1,f])])
				Δε[t,f] = rand(dε) - rand(dε)
				choices[t,f] = ΔU[getindX(m.supportX,X[t,f]),
									choices[t-1,f] + 1 ] > Δε[t,f]
			end

		end
	end

	return X,choices
end

"""
	estimate_model(m::DDCbasic,X,choices;
						startvals = zeros(3),
						estimateΠ=false)

Estimates parameters in model `m`, given observed payoff states `X` and choices `choices`. Returns new model with estimate parameter values.

## Optional keyword arguments
`startvals`: Set 3x1 vector of starting values
`estimateΠ`: Currently not implemented, would first estimate `Π` and then use that in estimation.
"""
function estimate_model(m::DDCbasic,X,choices;
							startvals = zeros(3),
							estimateΠ=false)

	# For parameter vector θ, return loglikelihood by creating new model struct
	# with these values
	function f(θ,m::DDCbasic,X,choices)
		# Get parameters from parameter vector
		β = θ[1:2]
		δ = [m.δ[1];θ[3]]
		# Construct model with parameters
		_m = typeof(m)(;supportX=m.supportX,Π=m.Π,β,δ,ρ=m.ρ)
		# Return neg loglikelihood
		return - loglikelihood(_m,X,choices)
	end

	# Wrap into objective function. Note, p is unused,
	# but required by Optimization.jl
	function f(θ,p)
		f(θ,m,X,choices)
	end

	# Define objective function and problem for Optimization.jl
	of = OptimizationFunction(f,m.options.AD)
	prob = OptimizationProblem(of,startvals)

	# Find optimum
	sol = solve(prob,m.options.algo)

	# Put together in new estimated model
	mout = deepcopy(m)
	mout.β .= sol.u[1:2]
	mout.δ[2] = sol.u[3]

	return mout
end

"""
	calc_stderr(m::DDCbasic,X,choices)

For model `m`, calculate standard errors for parameters given data `X,choices`.
"""
function calc_stderr(m::DDCbasic,X,choices)

	# Define function for parameter values in vector θ
	function f(θ,m::DDCbasic,X,choices)
		β = θ[1:2]
		δ = [m.δ[1];θ[3]]
		_m = typeof(m)(;supportX=m.supportX,Π=m.Π,β,δ,ρ=m.ρ)
		return - loglikelihood(_m,X,choices)
	end

	# Wrap into function depending only on parameter vector θ
	function f(θ)
		f(θ,m,X,choices)
	end

	# Return standard errors (diagonal elements of
	# inverse hessian are variance estimates)
   sqrt.(diag(inv(ForwardDiff.hessian(f,vcat(m.β,m.δ[2])))))
end


text:
# Autogenerated wrapper script for alsa_jll for powerpc64le-linux-gnu
export libasound, libatopology

JLLWrappers.@generate_wrapper_header("alsa")
JLLWrappers.@declare_library_product(libasound, "libasound.so.2")
JLLWrappers.@declare_library_product(libatopology, "libatopology.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libasound,
        "lib/libasound.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libatopology,
        "lib/libatopology.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DelimitedFiles

include("solver.jl")

function hilbert(m,n)
    return [1/(i+j-1) for i in 1:m, j in 1:n]
end

table2 = [(1000,1000),(500,500),(200,200),(100,100),(1000,2000),(500,2000),
(200,2000),(100,2000),(50,2000),(20,2000),(10,2000),(200,10000),(200,20000),
(200,50000),(100,20000),(100,100000),(100,200000),(50,100000),(50,200000),
(50,400000),(20,200000),(20,400000),(20,1000000),(10,100000),(10,200000),
(10,400000),(10,1000000),(10,2000000)]

for (m, n) in table2
    A = rand(m, n) .- 0.5
    z = ones(n)
    x0 = zeros(n)
    b = A * z
    local t0 = time_ns()
    y = systemSolver(A,b,x0)
    elapsed = round((time_ns()-t0)/1e9, digits=2)
    open("Table2.csv", "a") do io
        writedlm(io, [m n norm(y-x0,Inf) norm(b-A*y,Inf)/norm(b,Inf) elapsed], ',')
    end
end

table3 = [(200,200),(100,100),(50,2000),(10,2000),(200,10000),(200,20000),
(100,20000),(100,100000),(50,100000),(50,200000),(20,200000),(20,400000),
(10,100000),(10,200000),(10,400000),(10,1000000),(10,2000000)]

for (m, n) in table3
    A = hilbert(m, n)
    z = ones(n)
    x0 = zeros(n)
    b = A * z
    local t0 = time_ns()
    y = systemSolver(A,b,x0)
    elapsed = round((time_ns()-t0)/1e9, digits=2)
    open("Table3.csv", "a") do io
        writedlm(io, [m n norm(y-x0,Inf) norm(b-A*y,Inf)/norm(b,Inf) elapsed], ',')
    end
end

table4 = [(100,100),(200,200),(1000,1000),(2000,2000),(100,100000),(10,1000000)]

for (m, n) in table4
    A = rand(m, n) .- 0.5
    z = ones(n)
    x0 = zeros(n)
    b = A * z
    local t0 = time_ns()
    y = systemSolver(A,b,x0)
    elapsed = round((time_ns()-t0)/1e9, digits=2)
    open("Table4.csv", "a") do io
        writedlm(io, [m n norm(y-x0,Inf) norm(b-A*y,Inf)/norm(b,Inf) elapsed], ',')
    end
end

table5 = [(10,10),(20,20),(50,50),(100,100),(200,200),(100,100000),(10,1000000)]

for (m, n) in table5
    A = hilbert(m, n)
    z = ones(n)
    x0 = zeros(n)
    b = A * z
    local t0 = time_ns()
    y = systemSolver(A,b,x0)
    elapsed = round((time_ns()-t0)/1e9, digits=2)
    open("Table5.csv", "a") do io
        writedlm(io, [m n norm(y-x0,Inf) norm(b-A*y,Inf)/norm(b,Inf) elapsed], ',')
    end
end

max_stars_repo_name: amir-ni/ALA-extra

max_issues_repo_name: amir-ni/ALA-extra

max_forks_repo_name: amir-ni/ALA-extra


text:
"""
    struct CompilerInstance <: Any
Hold a pointer to a `clang::CompilerInstance` object.
"""
struct CompilerInstance
    ptr::CXCompilerInstance
end

Base.unsafe_convert(::Type{CXCompilerInstance}, x::CompilerInstance) = x.ptr
Base.cconvert(::Type{CXCompilerInstance}, x::CompilerInstance) = x


text:
## port methods of simplify module of SymPy
## from http://docs.sympy.org/0.7.2/modules/simplify/simplify.html

## simple methods (x, args) -> y (y coercion happens via PyCall)
simplify_sympy_meths = (:collect, :rcollect, :separate,
                        :separatevars,
                        :radsimp, :ratsimp, :trigsimp, :besselsimp,
                        :powsimp, :combsimp, :hypersimp,
                        :fraction,
                        :simplify, :nsimplify,
                        :posify, :powdenest, :sqrtdenest,
                        :logcombine, :hyperexpand)

expand_sympy_meths = (:expand_trig,
                      :expand_power_base, :expand_power_exp,
                      :expand_log,
                      :expand_func,
                      :hyperexpand
                      )

## special case cse due to output
## If this pattern is common, will need to collect and use metaprogramming.
@doc """
`$(cse)`: Perform common subexpression elimination on an expression.

[cf.](http://docs.sympy.org/latest/modules/simplify/simplify.html)

Example: (from man page)
```
@vars x w y z x0 x1
cse(((w + x + y + z)*(w + y + z))/(w + x)^3), ([(x0, y + z), (x1, w + x)], [(w + x0)*(x0 + x1)/x1^3]) # tuple of replacements and reduced expressions.

""" ->
function cse{T<:SymbolicObject}(ex::Union{T, AbstractArray{T}}, args...; kwargs...)
    a, b = sympy_meth(:cse, ex, args...; kwargs...)
    a, oftype(ex, b[1])
end
export(cse)

## didn't do traversal tools, EPath tools


text:
using UnicodePlots
import Pkg
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
using SpikeSynchrony
using Statistics
using JLD
using Distributed
using SharedArrays
using Plots
using UnicodePlots
using Evolutionary
using Distributions
using LightGraphs
using Metaheuristics

##
# Override to function to include a state.
##
SNN.@load_units
#unicodeplots()

#outputs = simulate!(pop, T; cb = cb, inputs=input)

global Ne = 200;
global Ni = 50
global σee = 1.0
global pee = 0.5
global σei = 1.0
global pei = 0.5

function make_net_from_graph_structure(xx)#;

    xx = Int(round(xx))
    #@show(xx)
    #h = turan_graph(xx, xx)#, seed=1,cutoff=0.3)

    h = circular_ladder_graph(xx)#, xx)#, seed=1,cutoff=0.3)
    #hi = circular_ladder_graph(xx)#, seed=1,cutoff=0.3)
    E = SNN.IZ(; N = Ne, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
    I = SNN.IZ(; N = Ni, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
    #EE = SNN.SpikingSynapse(E, E, :v; σ = σee, p = 1.0)
    EI = SNN.SpikingSynapse(E, I, :v; σ = σei, p = 1.0)
    IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1.0)
    II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 1.0)
    # PINningSynapse
    P = [E, I]#, EEA]
    C = [EI, IE, II]#, EEA]
    #EE = SNN.PINningSynapse(E, E, :v; σ=0.5, p=0.8)
    #for n in 1:(N - 1)
    #    SNN.connect!(EE, n, n + 1, 50)
    #end
    #for (i,j) in enumerate(h.fadjlist) println(i,j) end
    EE = SNN.SpikingSynapse(E, E, :v; σ = 0.5, p = 0.8)

    @inbounds for (i, j) in enumerate(h.fadjlist)
        @inbounds for k in j
            SNN.connect!(EE, i, k, 10)
        end
    end

    @inbounds for (i, j) in enumerate(hi.fadjlist)
        @inbounds for k in j
            if i < Ni && k < Ni

                SNN.connect!(EI, i, k, 10)
                SNN.connect!(IE, i, k, 10)
                SNN.connect!(II, i, k, 10)
            end
        end
    end

    #for (i,j) in enumerate(h.fadjlist)
    #    for k in j
    #        SNN.connect!(EI, i, k, 50)
    #    end
    #end
    P = [E, I]#, EEA]
    C = [EE, EI, IE, II]#, EEA]
    return P, C

end

function make_net(Ne, Ni; σee = 1.0, pee = 0.5, σei = 1.0, pei = 0.5)
    Ne = 200
    Ni = 50

    E = SNN.IZ(; N = Ne, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
    I = SNN.IZ(; N = Ni, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
    EE = SNN.SpikingSynapse(E, E, :v; σ = σee, p = pee)
    EI = SNN.SpikingSynapse(E, I, :v; σ = σei, p = pei)
    IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 0.5)
    II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.5)
    P = [E, I]#, EEA]
    C = [EE, EI, IE, II]#, EEA]
    return P, C
end
function get_trains(p)
    fire = p.records[:fire]
    x, y = Float32[], Float32[]
    for time in eachindex(fire)
        for neuron_id in findall(fire[time])
            push!(x, time)
            push!(y, neuron_id)
        end
    end
    cellsa = Array{Union{Missing,Any}}(undef, 1, Int(findmax(y)[1]))
    nac = Int(findmax(y)[1])
    for (inx, cell_id) in enumerate(1:nac)
        cellsa[inx] = []
    end
    @inbounds for cell_id in unique(y)
        @inbounds for (time, cell) in collect(zip(x, y))
            if Int(cell_id) == cell
                append!(cellsa[Int(cell_id)], time)

            end

        end
    end

    cellsa

end
global E
global spkd_ground

#P, C = make_net(Ne, Ni, σee = 0.5, pee = 0.8, σei = 0.5, pei = 0.8, a = 0.02)
#sggcu =[ CuArray(convert(Array{Float32,1},sg)) for sg in spkd_ground ]

#Flux.SGD
#Flux.gpu

function rmse(spkd)
    error = Losses(mean(spkd), spkd; agg = mean)
end

function rmse_depr(spkd)
    total = 0.0
    @inbounds for i = 1:size(spkd, 1)
        total += (spkd[i] - mean(spkd[i]))^2.0
    end
    return sqrt(total / size(spkd, 1))
end

function raster_difference(spkd0, spkd_found)
    maxi0 = size(spkd0)[2]
    maxi1 = size(spkd_found)[2]
    mini = findmin([maxi0, maxi1])[1]
    spkd = ones(mini)
    maxi = findmax([maxi0, maxi1])[1]

    if maxi > 0
        if maxi0 != maxi1
            return sum(ones(maxi))

        end
        if isempty(spkd_found[1, :])
            return sum(ones(maxi))
        end
    end
    spkd = ones(mini)
    @inbounds for i in eachindex(spkd)
        if !isempty(spkd0[i]) && !isempty(spkd_found[i])
            maxt1 = findmax(spkd0[i])[1]
            maxt2 = findmax(spkd_found[i])[1]
            maxt = findmax([maxt1, maxt2])[1]
            if maxt1 > 0.0 && maxt2 > 0.0
                t, S = SpikeSynchrony.SPIKE_distance_profile(
                    unique(sort(spkd0[i])),
                    unique(sort(spkd_found[i]));
                    t0 = 0.0,
                    tf = maxt,
                )
                spkd[i] = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1])
            end
        end
    end
    spkd
end

#=
using Base.Threads

function raster_difference_threads(spkd0, spkd_found)
    maxi0 = size(spkd0)[2]
    maxi1 = size(spkd_found)[2]
    mini = findmin([maxi0, maxi1])[1]
    spkd = ones(mini)
    maxi = findmax([maxi0, maxi1])[1]
    if maxi > 0
        if maxi0 != maxi1
            return sum(ones(maxi))

        end
        if isempty(spkd_found[1, :])
            return sum(ones(maxi))
        end
    end
    spkd = ones(mini)
    N = length(spkd)
    Threads.foreach(
        Channel(nthreads()) do chnl
            for i in 1:N
                val = 1.0
                if !isempty(spkd0[i]) && !isempty(spkd_found[i])
                    maxt1 = findmax(spkd0[i])[1]
                    maxt2 = findmax(spkd_found[i])[1]
                    maxt = findmax([maxt1, maxt2])[1]
                    if maxt1 > 0.0 && maxt2 > 0.0
                        t, S = SpikeSynchrony.SPIKE_distance_profile(
                            unique(sort(spkd0[i])),
                            unique(sort(spkd_found[i]));
                            t0 = 0.0,
                            tf = maxt,
                        )
                        val = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1])
                    end
                end
                put!(chnl, val)
                println("thread ", threadid(), ": ", "put ", val)
            end
        end
    end
    spkd
end
=#

function loss(model)
    #σee = model[1]
    #pee = model[2]
    #σei = model[3]
    #pei = model[4]
    #P1, C1 = make_net(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    #@show(model)
    println("best candidate ", 26)
    println(" ")
    #println("found ", model[1])
    P1, C1 = make_net_SNN(model[1])#,a=a)

    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])#vec(E_stim[t,:])#[i]#3randn(Ne)
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    #println("Ground Truth \n")
    #SNN.raster([E]) |> display
    #println("Best Candidate \n")

    #SNN.raster([E1]) |> display

    error = raster_difference(spkd_ground, spkd_found)
    error = sum(error)
    #@show(error)
    #println("Broke here")

    error

end

function eval_best(params)
    xx = Int(round(params[1]))
    @show(xx)
    P1, C1 = make_net_SNN(xx)#,a=a)
    println("found ", xx)

    #σee = params[1]
    #pee = params[2]
    #σei = params[3]
    #pei = params[4]
    #P1, C1 = make_net(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])#vec(E_stim[t,:])#[i]#3randn(Ne)
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth: \n")
    SNN.raster([E]) |> display
    println("candidate: \n")

    SNN.raster([E1]) |> display
    #error = raster_difference(spkd_ground,spkd_found)
    E1, spkd_found

end

function init_b(lower, upper)
    gene = []

    for (i, (l, u)) in enumerate(zip(lower, upper))
        p1 = rand(l:u, 1)
        append!(gene, p1)
    end
    gene
end

function initf(n)
    genesb = []
    for i = 1:n
        genes = init_b(lower, upper)
        append!(genesb, [genes])
    end
    genesb
end

function initd()
    population = initf(10)
    garray = zeros((length(population)[1], length(population[1])))
    for (i, p) in enumerate(population)
        garray[i, :] = p
    end
    garray[1, :]
end
function Evolutionary.trace!(
    record::Dict{String,Any},
    objfun,
    state,
    population,
    method::GA,
    options,
)
    idx = sortperm(state.fitpop)
    record["fitpop"] = state.fitpop[:]#idx[1:last(idx)]]
    record["pop"] = population[:]
    #record["σ"] = state.
end

#lower = Float32[0.0 0.0 0.0 0.0]# 0.03 4.0]
#upper = Float32[1.0 1.0 1.0 1.0]# 0.2 20.0]

#lower = Int32[3]# 0.0 0.0 0.0]# 0.03 4.0]
#upper = Int32[40]# 1.0 1.0 1.0]# 0.2 20.0]

#lower = Float32[0.0 0.0 0.0 0.0]# 0.03 4.0]
#upper = Float32[1.0 1.0 1.0 1.0]# 0.2 20.0]
#lower = vec(lower)
#upper = vec(upper)
#using Evolutionary, MultivariateStats
#range = Any[lower,upper]

# mutation = domainrange(fill(0.5,4))
#
#function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::GA, options)
#    idx = sortperm(state.fitpop)
#record["population"] = population
#    state.fitpop[idx[1:5]]
#    record["fitpop"] = state.fitpop[idx[1:5]]
#end
#using Evolutionary

#=
function Evolutionary.value!(::Val{:serial}, fitness, objfun, population::AbstractVector{IT}) where {IT}
    @show(typeof(fitness))
    fitness = SharedArrays.SharedArray{Float32}(fitness)
    @time @sync @distributed for i in 1:length(population)

    #pmap(f, [::AbstractWorkerPool], c...; distributed=true, batch_size=) -> collection
        fitness[i] = value(objfun, population[i])
        #println("I'm worker $(myid()), working on i=$i")
    end
    @show(typeof(fitness))

    fitness
    #fitness = Array(fitness)
    #@show(fitness)
end
=#


text:
module TestCodeMath
  using FactCheck
  import Code
  import Creature
  import Creature.f16
  import Helper
  import Mutator
  import Config

  include("Helper.jl")
  #
  # Tests specified operator in a simple form var = var operator var
  #
  function testOperator(opSym::Symbol, opFn::Symbol, result::Float16)
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, opFn, org)
    @fact Helper.getArg(org.code, [2,Creature.VAR_AMOUNT + 1,2,1]) --> opSym
    @fact eval(org.code)(d.cfg, org) --> result

    Manager.destroy(d.man)
  end
  #
  # Tests specified operator two times
  #
  function testOperatorTwo(opSym::Symbol, opFn::Symbol, result::Float16)
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, opFn, org)
    code(d, opFn, org)
    @fact Helper.getArg(org.code, [2,Creature.VAR_AMOUNT + 1,2,1]) --> opSym
    @fact Helper.getArg(org.code, [2,Creature.VAR_AMOUNT + 2,2,1]) --> opSym
    @fact eval(org.code)(d.cfg, org) --> result

    Manager.destroy(d.man)
  end
  #
  # Tests specified operator inside if operator
  #
  function testOperatorInIf(opFn::Symbol, result::Float16)
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :condition, org)
    code(d, opFn, org, [1, 2, 1])
    #
    # Updates compare operator to "===" to make condition true all the time
    #
    Helper.getArg(org.code, [2,2,1]).args[1] = :(===)
    @fact eval(org.code)(d.cfg, org) --> result

    Manager.destroy(d.man)
  end
  #
  # Tests specified operator inside for operator
  #
  function testOperatorInFor(opFn::Symbol, result::Float16)
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :loop, org)
    code(d, opFn, org, [1, 2, 1])
    @fact eval(org.code)(d.cfg, org) --> result

    Manager.destroy(d.man)
  end
  #
  # plus
  #
  facts("Testing CodeMath.plus()") do; testOperator(:(+), :plus, Float16(10.0)) end
  facts("Testing CodeMath.plus() two times") do; testOperatorTwo(:(+), :plus, Float16(20.0)) end
  facts("Testing CodeMath.plus() inside if operator") do; testOperatorInIf(:plus, Float16(10.0)) end
  facts("Testing CodeMath.plus() inside for operator") do; testOperatorInFor(:plus, Float16(1280.0)) end
  #
  # minus
  #
  facts("Testing CodeMath.minus()") do; testOperator(:(-), :minus, Float16(0.0)) end
  facts("Testing CodeMath.minus() two times") do; testOperatorTwo(:(-), :minus, Float16(0.0)) end
  facts("Testing CodeMath.minus() inside if operator") do; testOperatorInIf(:minus, Float16(0.0)) end
  facts("Testing CodeMath.minus() inside for operator") do; testOperatorInFor(:minus, Float16(0.0)) end
  #
  # multiply
  #
  facts("Testing CodeMath.multiply()") do; testOperator(:(*), :multiply, Float16(25.0)) end
  facts("Testing CodeMath.multiply() two times") do; testOperatorTwo(:(*), :multiply, Float16(625.0)) end
  facts("Testing CodeMath.multiply() inside if operator") do; testOperatorInIf(:multiply, Float16(25.0)) end
  facts("Testing CodeMath.multiply() inside for operator") do; testOperatorInFor(:multiply, Float16(65504.0)) end
  #
  # divide
  #
  facts("Testing CodeMath.divide()") do; testOperator(:(/), :divide, Float16(1.0)) end
  facts("Testing CodeMath.divide() two times") do; testOperatorTwo(:(/), :divide, Float16(1.0)) end
  facts("Testing CodeMath.divide() inside if operator") do; testOperatorInIf(:divide, Float16(1.0)) end
  facts("Testing CodeMath.divide() inside for operator") do; testOperatorInFor(:divide, Float16(1.0)) end
  #
  # not
  #
  facts("Testing CodeMath.not()") do; testOperator(:(*), :not, Float16(-5.0)) end
  facts("Testing CodeMath.not() two times") do; testOperatorTwo(:(*), :not, Float16(5.0)) end
  facts("Testing CodeMath.not() inside if operator") do; testOperatorInIf(:not, Float16(-5.0)) end
  facts("Testing CodeMath.not() inside for operator") do; testOperatorInFor(:not, Float16(5.0)) end
  #
  # and
  #
  facts("Testing CodeMath.and()") do; testOperator(:(&), :and, Float16(5.0)) end
  facts("Testing CodeMath.and() two times") do; testOperatorTwo(:(&), :and, Float16(5.0)) end
  facts("Testing CodeMath.and() inside if operator") do; testOperatorInIf(:and, Float16(5.0)) end
  facts("Testing CodeMath.and() inside for operator") do; testOperatorInFor(:and, Float16(5.0)) end
  #
  # or
  #
  facts("Testing CodeMath.or()") do; testOperator(:(|), :or, Float16(5.0)) end
  facts("Testing CodeMath.or() two times") do; testOperatorTwo(:(|), :or, Float16(5.0)) end
  facts("Testing CodeMath.or() inside if operator") do; testOperatorInIf(:or, Float16(5.0)) end
  facts("Testing CodeMath.or() inside for operator") do; testOperatorInFor(:or, Float16(5.0)) end
  #
  # xor
  #
  facts("Testing CodeMath.xor()") do; testOperator(:($), :xor, Float16(0.0)) end
  facts("Testing CodeMath.xor() two times") do; testOperatorTwo(:($), :xor, Float16(0.0)) end
  facts("Testing CodeMath.xor() inside if operator") do; testOperatorInIf(:xor, Float16(0.0)) end
  facts("Testing CodeMath.xor() inside for operator") do; testOperatorInFor(:xor, Float16(0.0)) end
  #
  # rshift
  #
  facts("Testing CodeMath.rshift()") do; testOperator(:(>>), :rshift, Float16(0.0)) end
  facts("Testing CodeMath.rshift() two times") do; testOperatorTwo(:(>>), :rshift, Float16(0.0)) end
  facts("Testing CodeMath.rshift() inside if operator") do; testOperatorInIf(:rshift, Float16(0.0)) end
  facts("Testing CodeMath.rshift() inside for operator") do; testOperatorInFor(:rshift, Float16(0.0)) end
  #
  # lshift
  #
  facts("Testing CodeMath.lshift()") do; testOperator(:(<<), :lshift, Float16(160.0)) end
  facts("Testing CodeMath.lshift() two times") do; testOperatorTwo(:(<<), :lshift, Float16(0.0)) end
  facts("Testing CodeMath.lshift() inside if operator") do; testOperatorInIf(:lshift, Float16(160.0)) end
  facts("Testing CodeMath.lshift() inside for operator") do; testOperatorInFor(:lshift, Float16(0.0)) end
  #
  # reminder
  #
  facts("Testing CodeMath.reminder()") do; testOperator(:(%), :reminder, Float16(0.0)) end
  facts("Testing CodeMath.reminder() two times") do; testOperatorTwo(:(%), :reminder, Float16(0.0)) end
  facts("Testing CodeMath.reminder() inside if operator") do; testOperatorInIf(:reminder, Float16(0.0)) end
  facts("Testing CodeMath.reminder() inside for operator") do; testOperatorInFor(:reminder, Float16(0.0)) end
  #
  # sqrt
  #
  facts("Testing CodeMath.sqrt()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :sqrt, org)
    @fact eval(org.code)(d.cfg, org) --> Float16(2.236)

    Manager.destroy(d.man)
  end
  facts("Testing CodeMath.sqrt() with negative value") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :sqrt, org)
    val(org, Float16(-5.0))
    @fact eval(org.code)(d.cfg, org) --> Float16(2.236)

    Manager.destroy(d.man)
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
const VERTICES_1 = [[1]]
const VERTICES_2 = [[1,2], [1], [2], [3]]
const VERTICES_3 = [
    [1,2,4], [1,4], [1,2], [2,4], [2,5], [1,6], [3,4], [2], [3,5,6], [1], [4],
    [3,6], [3,5], [5,6], [5], [3], [6], [7]
]
const VERTICES_4 = [
    [1,2,4,8], [1,4,8], [1,2,8], [1,2,4], [2,4,8], [1,6,8], [1,4,10], [1,2,12],
    [2,5,8], [2,4,9], [3,4,8], [1,4], [1,2], [2,8], [2,5,9,12], [2,4],
    [1,6,10,12], [3,5,6,8], [3,4,9,10], [1,8], [4,8], [2,9,12], [2,5,12],
    [2,5,9], [1,6,10], [1,6,12], [1,10,12], [3,6,8], [3,5,8], [3,5,6,9,10,12],
    [3,4,10], [3,4,9], [4,9,10], [5,6,8], [3,8], [3,6,9,10,12], [2,12],
    [3,5,9,10,12], [1,6], [3,5,6,10,12], [3,5,6,9,12], [3,5,6,9,10], [1,10],
    [3,4], [2,5], [1,12], [4,10], [4,9], [2,9], [5,8], [5,6,9,10,12], [6,8],
    [3,6,10,12], [3,6,9,12], [3,6,9,10], [3,5,10,12], [3,5,9,12], [3,5,9,10],
    [3,5,6,12], [3,5,6,10], [3,5,6,9], [1,14], [3,9,10,12], [4,11], [2,13],
    [5,9,10,12], [5,6,10,12], [5,6,9,12], [5,6,9,10], [6,9,10,12], [7,8],
    [3,5,12], [3,10,12], [3,5,10], [3,9,12], [3,9,10], [3,6,9], [1], [2], [4],
    [3,5,9,14], [5,10,12], [5,9,12], [5,9,10], [3,6,12], [5,6,11,12], [5,6,10],
    [5,6,9], [6,10,12], [6,9,12], [6,9,10], [3,5,6], [3,6,10,13], [7,9,10,12],
    [8], [5,10], [3,9,14], [5,9,14], [3,10,13], [3,5,9], [5,6,12], [5,6,11],
    [3,6,13], [3,12], [3,6,10], [6,11,12], [6,10,13], [3,5,14], [6,9],
    [7,10,12], [7,9,12], [7,9,10], [5,11,12], [9,10,12], [6,12], [6,11,13],
    [6,10], [3,6], [3,13,14], [5,12], [7,11,12], [7,10,13], [5,9], [7,9,14],
    [5,6], [3,10], [5,11,14], [3,5], [3,9], [9,12], [9,10], [10,12], [7,12],
    [7,11,13,14], [7,10], [3,14], [6,11], [7,9], [6,13], [5,14], [3,13],
    [5,11], [9,14], [10,13], [11,12], [6], [9], [7,11,14], [3], [10], [5],
    [7,13,14], [11,13,14], [7,11,13], [12], [7,13], [11,14], [11,13], [7,11],
    [7,14], [13,14], [11], [13], [7], [14], [15]
]


text:
## Do not precompile. The pointer to gmp functions needs to be found at run time.

#module PolynomialSequences

using Polynomials
using Memoize

export fibpoly, lucaspoly

fibpoly(n::Int) = _fib_poly(n)
fibpoly(n::BigInt) = _fib_poly_big(n)

lucaspoly(n::Int) = _lucas_poly(n)
lucaspoly(n::BigInt) = _lucas_poly_big(n)

## TODO: refactor code

## This requires Dict maybe ?
let myzero = big(0), myone = big(1), xvar = Polynomials.Poly([myzero,myone]), zerovar = Polynomials.Poly([myzero]),  onevar = Polynomials.Poly([myone])
    global _fib_poly_big
    @memoize Dict function _fib_poly_big(n::BigInt)
        if n == 0
            return zerovar
        elseif n == 1
            return onevar
        else
            return xvar *  _fib_poly_big(n-1) +  _fib_poly_big(n-2)
        end
    end
end

let myzero = 0, myone = 1, xvar = Polynomials.Poly([myzero,myone]), zerovar = Polynomials.Poly([myzero]),  onevar = Polynomials.Poly([myone])
    global _fib_poly
    @memoize function _fib_poly(n::Int)
        if n == 0
            return zerovar
        elseif n == 1
            return onevar
        else
            return xvar *  _fib_poly(n-1) +  _fib_poly(n-2)
        end
    end
end

let  myzero = big(0), myone = big(1), mytwo = big(2), xvar = Polynomials.Poly([myzero,myone]),  twovar = Polynomials.Poly([mytwo])
    global _lucas_poly_big
    @memoize Dict function _lucas_poly_big(n::BigInt)
        if n == 0
            return twovar
        elseif n == 1
            return xvar
        else
            return xvar *  _lucas_poly_big(n-1) +  _lucas_poly_big(n-2)
        end
    end
end

let myzero = 0, myone = 1, mytwo = 2, xvar = Polynomials.Poly([myzero,myone]),  twovar = Polynomials.Poly([mytwo])
    global _lucas_poly
    @memoize function _lucas_poly(n::Int)
        if n == 0
            return twovar
        elseif n == 1
            return xvar
        else
            return xvar * _lucas_poly(n-1) +  _lucas_poly(n-2)
        end
    end
end

#end # module


text:
"""
    scan_hdf5(hdf)

Scan HDF5 archive `hdf` and return station names with waveform data contained
therein as a list of Strings formatted "nn.sssss" (network.station).

    scan_hdf5(hdf, level="trace")

Scan HDF5 archive `hdf` and return channel names with waveform data contained
therein as a list of Strings formatted "nn.sssss.ll.ccc" (network.station.location.channel).
"""
function scan_hdf5( hdf::String;
  fmt::String="asdf",
  level::String="station")

  f = h5open(hdf, "r")
  if fmt =="asdf"
    if level == "station"
      str = names(f["Waveforms"])
    elseif level == "trace"
      D = get_datasets(f)
      str = String[]
      for i in D
        push!(str, name(i))
      end
    else
      error("unsupported level!")
    end
  else
    error("unknown format or NYI!")
  end
  close(f)
  unique!(str)
  return str
end


text:
using DiscreteFunctions, SimpleGraphs

"""
`tree2function(T::SimpleGraph)` converts a tree `T` into a
`DiscreteFunction`. We assume the vertex set of `T` is of the
form `{1,2,...,n}`. The resulting function has each vertex pointing
to the next vertex on a path to the vertex `1`.

**Warning**: I don't check if the input graph is valid. That is
it *must* be a tree with vertex set `1:n`.
"""
function tree2function(G::SimpleGraph)::DiscreteFunction
    n = NV(G)
    data = zeros(Int,n)
    data[1] = 1
    for v = 2:n
        P = find_path(G,v,1)
        data[v] = P[2]
    end
    return DiscreteFunction(data)
end

"""
`RandomTreeFunction(n)` creates a random function whose digraph
is a tree (with a single selfloop).
"""
function RandomTreeFunction(n::Int)::DiscreteFunction
    T = RandomTree(n)
    return tree2function(T)
end

"""
`is_tree_function(f::DiscreteFunction)` returns `true` exactly when the
graph of `f` is a tree (with a single loop at one vertex).
"""
function is_tree_function(f::DiscreteFunction)::Bool
    clist = cycles(f)
    if length(clist) != 1
        return false
    end
    if length(clist[1]) != 1
        return false
    end
    return true
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of filters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Filters" begin
    expected = """
┌─────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 3 │
├─────┼────────┼────────┤
│   2 │      2 │    2.0 │
│   4 │      4 │    4.0 │
│   6 │      6 │    6.0 │
└─────┴────────┴────────┘
"""

    result = pretty_table(
        String,
        data;
        filters_row     = ((data,i) -> i%2 == 0,),
        filters_col     = ((data,i) -> i%2 == 1,),
        formatters      = ft_printf("%.1f", 3),
        show_row_number = true
    )
    @test result == expected

    expected = """
┌─────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 3 │
├─────┼────────┼────────┤
│   2 │   2    │ 2.0    │
│   4 │   4    │ 4.0    │
│   6 │   6    │ 6.0    │
└─────┴────────┴────────┘
"""

    result = pretty_table(
        String,
        data;
        filters_row     = ((data,i) -> i%2 == 0,),
        filters_col     = ((data,i) -> i%2 == 1,),
        formatters      = ft_printf("%.1f", 3),
        show_row_number = true,
        alignment       = [:c, :l, :l, :c]
    )
    @test result == expected

    # No data after filtering
    # --------------------------------------------------------------------------

    result = pretty_table(
        String,
        data;
        filters_row     = ((data,i) -> false,),
        filters_col     = ((data,i) -> i % 2 == 1,),
        formatters      = ft_printf("%.1f", 3),
        show_row_number = true,
        alignment       = [:c, :l, :l, :c]
    )
    @test result == ""

    result = pretty_table(
        String,
        data;
        filters_row     = ((data,i) -> false,),
        filters_col     = ((data,i) -> i % 2 == 1,),
        formatters      = ft_printf("%.1f", 3),
        show_row_number = true,
        alignment       = [:c, :l, :l, :c],
        title           = "Empty"
    )
    @test result == "Empty\n"
end