text:
@testset "ObsView" begin
    @test ObsView <: AbstractVector
    @test ObsView <: LearnBase.DataView
    @test ObsView <: LearnBase.AbstractObsView
    @test ObsView <: LearnBase.AbstractObsIterator
    @test ObsView <: LearnBase.AbstractDataIterator
    @test obsview === ObsView

    @testset "constructor" begin
        @test_throws DimensionMismatch ObsView((rand(2,10),rand(9)))
        @test_throws DimensionMismatch ObsView((rand(2,10),rand(4,9,10),rand(9)))
        @test_throws MethodError ObsView(EmptyType())
        @test_throws MethodError ObsView(EmptyType(), ObsDim.Last())
        @test_throws MethodError ObsView(EmptyType(), ObsDim.Undefined())
        @test_throws MethodError ObsView(EmptyType(), ObsDim.Last())
        @test_throws MethodError ObsView((EmptyType(),EmptyType()))
        @test_throws MethodError ObsView(CustomType(), ObsDim.Last())
        @test_throws MethodError ObsView(EmptyType(), obsdim=1)
        @test_throws MethodError ObsView((EmptyType(),EmptyType()))
        for var in (vars..., Xs, ys)
            A = @inferred(ObsView(var))
            @test A.obsdim == ObsDim.Last()
        end
        for var in tuples
            A = @inferred(ObsView(var))
            @test A.obsdim == (fill(ObsDim.Last(),length(var))...,)
        end
        for var in (vars..., tuples..., Xs, ys)
            A = ObsView(var)
            @test @inferred(parent(A)) === var
            @test @inferred(ObsView(A)) == A
            @test @inferred(ObsView(var)) == A
        end
        A = ObsView(X',obsdim=1)
        @test A == @inferred(ObsView(X',ObsDim.First()))
        @test A == ObsView(X',obsdim=:first)
    end

    println("<HEARTBEAT>")

    @testset "typestability" begin
        @test_throws Exception @inferred(ObsView(X, obsdim=2))
        for var in (vars..., tuples..., Xs, ys)
            @test typeof(@inferred(ObsView(var))) <: ObsView
            @test typeof(@inferred(ObsView(var, ObsDim.Last()))) <: ObsView
            @test typeof(@inferred(ObsView(var))) <: ObsView
            @test_throws Exception @inferred(ObsView(var, obsdim=:last))
        end
        for tup in tuples
            @test typeof(@inferred(ObsView(tup))) <: ObsView
            @test typeof(@inferred(ObsView(tup,(fill(ObsDim.Last(),length(tup))...,)))) <: ObsView
            @test typeof(@inferred(ObsView(tup,(fill(ObsDim.Last(),length(tup))...,)))) <: ObsView
            @test_throws Exception @inferred(ObsView(tup, obsdim=:last))
        end
        @test typeof(@inferred(ObsView(CustomType()))) <: ObsView
        @test typeof(@inferred(ObsView(CustomType(), ObsDim.Undefined()))) <: ObsView
    end

    println("<HEARTBEAT>")

    @testset "AbstractArray interface" begin
        for var in (vars..., tuples..., Xs, ys)
            A = ObsView(var)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[151]
            @test @inferred(nobs(A)) == 150
            @test @inferred(length(A)) == 150
            @test @inferred(size(A)) == (150,)
            @test @inferred(A[2:3]) == ObsView(datasubset(var, 2:3))
            @test @inferred(A[[1,3]]) == ObsView(datasubset(var, [1,3]))
            @test @inferred(A[1]) == datasubset(var, 1)
            @test @inferred(A[111]) == datasubset(var, 111)
            @test @inferred(A[150]) == datasubset(var, 150)
            @test A[end] == A[150]
            @test @inferred(getobs(A,1)) == getobs(var, 1)
            @test @inferred(getobs(A,111)) == getobs(var, 111)
            @test @inferred(getobs(A,150)) == getobs(var, 150)
            @test typeof(@inferred(collect(A))) <: Vector
        end
        for var in vars
            A = ObsView(var)
            @test @inferred(getobs(A)) == [getobs(var,i) for i in 1:nobs(var)]
        end
        for var in tuples
            A = ObsView(var)
            @test @inferred(getobs(A)) == [map(x->getobs(x,i), var) for i in 1:nobs(var)]
        end
        A = ObsView(X',obsdim=1)
        @test @inferred(length(A)) == 150
        @test @inferred(size(A)) == (150,)
        @test @inferred(A[2:3]) == ObsView(datasubset(X', 2:3, obsdim=1),obsdim=1)
        @test @inferred(A[1]) == datasubset(X', 1, obsdim=1)
        @test @inferred(A[111]) == datasubset(X', 111, obsdim=1)
        @test @inferred(A[150]) == datasubset(X', 150, obsdim=1)
        @test A[end] == A[150]
        @test @inferred(getobs(A,1)) == getobs(X', 1, obsdim=1)
        @test @inferred(getobs(A,111)) == getobs(X', 111, obsdim=1)
        @test @inferred(getobs(A,150)) == getobs(X', 150, obsdim=1)
    end

    println("<HEARTBEAT>")

    @testset "subsetting" begin
        for var_raw in (vars..., tuples..., Xs, ys)
            for var in (var_raw, DataSubset(var_raw))
                A = ObsView(var)
                @test getobs(@inferred(datasubset(A))) == @inferred(getobs(A))
                S = @inferred(datasubset(A, 1:5))
                @test typeof(S) <: ObsView
                @test @inferred(length(S)) == 5
                @test @inferred(size(S)) == (5,)
                @test @inferred(A[1:5]) == S
                @test @inferred(getobs(A,1:5)) == getobs(S)
                @test @inferred(getobs(S)) == getobs(ObsView(datasubset(var,1:5)))
                S = @inferred(DataSubset(A, 1:5))
                @test typeof(S) <: ObsView
                @test typeof(S.data) <: Union{DataSubset,Tuple}
                @test @inferred(length(S)) == 5
                @test @inferred(size(S)) == (5,)
                @test @inferred(getobs(S)) == getobs(ObsView(DataSubset(var,1:5)))
            end
        end
        A = ObsView(X)
        @test typeof(A.data) <: Array
        S = @inferred(datasubset(A))
        @test typeof(S) <: ObsView
        @test @inferred(length(S)) == 150
        @test typeof(S.data) <: SubArray
    end

    @testset "iteration" begin
        count = 0
        for (i,x) in enumerate(ObsView(X1))
            @test all(i .== x)
            count += 1
        end
        @test count == 150
    end
end

@testset "_compute_batch_settings" begin
    @test MLDataPattern._compute_batch_settings(X) === (30,5)
    @test MLDataPattern._compute_batch_settings(X',-1,-1,ObsDim.First()) === (30,5)
    @test MLDataPattern._compute_batch_settings(X',32,-1,ObsDim.First()) === (32,4)
    @test MLDataPattern._compute_batch_settings(X',32,-1,ObsDim.First(),true) === (30,5)
    @test MLDataPattern._compute_batch_settings(Xv) === (30,5)
    @test MLDataPattern._compute_batch_settings(Xs) === (30,5)
    @test MLDataPattern._compute_batch_settings(DataSubset(X)) === (30,5)
    @test MLDataPattern._compute_batch_settings((X,y)) === (30,5)
    @test MLDataPattern._compute_batch_settings((X',y),-1,-1,ObsDim.First()) === (30,5)
    @test MLDataPattern._compute_batch_settings((Xv,yv)) === (30,5)
    @test MLDataPattern._compute_batch_settings(X,160,-1,ObsDim.Last(),true) === (150,1)

    @test_throws ArgumentError MLDataPattern._compute_batch_settings(X, 160)
    @test_throws ArgumentError MLDataPattern._compute_batch_settings(X, 1, 160)
    @test_throws ArgumentError MLDataPattern._compute_batch_settings(X, 10, 20)

    for inner in (Xs, ys, vars...), var in (inner, DataSubset(inner))
        @test MLDataPattern._compute_batch_settings(var,10) === (10,15)
        @test MLDataPattern._compute_batch_settings(var,0,10) === (15,10)
        @test MLDataPattern._compute_batch_settings(var,-1,10) === (15,10)
        @test MLDataPattern._compute_batch_settings(var,10,10) === (10,15)
        @test MLDataPattern._compute_batch_settings(var,150,1) === (150,1)
        @test MLDataPattern._compute_batch_settings(var,150) === (150,1)
        @test MLDataPattern._compute_batch_settings(var,0,150) === (1,150)
        @test MLDataPattern._compute_batch_settings(var,-1,150) === (1,150)
    end
end

@testset "BatchView" begin
    @test BatchView <: AbstractVector
    @test BatchView <: LearnBase.DataView
    @test BatchView <: LearnBase.AbstractBatchView
    @test BatchView <: LearnBase.AbstractBatchIterator
    @test BatchView <: LearnBase.AbstractDataIterator
    @test batchview == BatchView
    @test_throws MethodError oversample(BatchView(X))
    @test_throws MethodError undersample(BatchView(X))
    @test_throws MethodError stratifiedobs(BatchView(X))

    @testset "constructor" begin
        @test_throws DimensionMismatch BatchView((rand(2,10),rand(9)))
        @test_throws DimensionMismatch BatchView((rand(2,10),rand(9)))
        @test_throws DimensionMismatch BatchView((rand(2,10),rand(4,9,10),rand(9)))
        @test_throws MethodError BatchView(EmptyType())
        @test_throws MethodError BatchView(EmptyType(), ObsDim.Last())
        @test_throws MethodError BatchView(EmptyType(), 10, ObsDim.Last())
        @test_throws MethodError BatchView(EmptyType(), 10, 15, ObsDim.Undefined())
        @test_throws MethodError BatchView(EmptyType(), 5, ObsDim.Last())
        @test_throws MethodError BatchView((EmptyType(),EmptyType()))
        @test_throws MethodError BatchView(CustomType(),5, ObsDim.Last())
        @test_throws MethodError BatchView(EmptyType(), obsdim=1)
        for var in (vars..., tuples..., Xs, ys)
            @test_throws MethodError BatchView(var...)
            @test_throws MethodError BatchView(var..., obsdim=:last)
            @test_throws ArgumentError BatchView(var, 151)
            A = BatchView(var, maxsize=151)
            @test length(A) == 1
            @test batchsize(A) == 150
            @test nobs(A) == nobs(var)
            @test @inferred(parent(BatchView(var))) === var
            A = @inferred(BatchView(var))
            @test @inferred(BatchView(A)) == A
            @test typeof(BatchView(A)) <: typeof(A)
            @test @inferred(BatchView(var)) == A
            @test nobs(A) == nobs(var)
            @test batchsize(A) == 30
        end
        @test BatchView((X,X)) == @inferred(BatchView((X,X), (ObsDim.Last(),ObsDim.Last())))
        @test BatchView((X,X)) == @inferred(BatchView((X,X), -1, (ObsDim.Last(),ObsDim.Last())))
        A = BatchView(X',obsdim=1)
        @test A == @inferred(BatchView(X',ObsDim.First()))
        @test A == @inferred(BatchView(X',-1,ObsDim.First()))
        @test A == BatchView(X',obsdim=:first)
    end

    println("<HEARTBEAT>")

    @testset "typestability" begin
        for var in (vars..., tuples..., Xs, ys)
            @test typeof(@inferred(BatchView(var))) <: BatchView
            @test typeof(@inferred(BatchView(var, 10))) <: BatchView
            @test typeof(@inferred(BatchView(var, 10, 15))) <: BatchView
            @test typeof(@inferred(BatchView(var, -1, 10))) <: BatchView
            @test typeof(@inferred(BatchView(var, 10, ObsDim.Last()))) <: BatchView
            @test typeof(@inferred(BatchView(var))) <: BatchView
        end
        for tup in tuples
            @test typeof(@inferred(BatchView(tup))) <: BatchView
            @test typeof(@inferred(BatchView(tup, 5))) <: BatchView
            @test typeof(@inferred(BatchView(tup, -1, 10))) <: BatchView
        end
        @test typeof(@inferred(BatchView(CustomType()))) <: BatchView
        @test typeof(@inferred(BatchView(CustomType(), 10, ObsDim.Undefined()))) <: BatchView
    end

    println("<HEARTBEAT>")

    @testset "AbstractArray interface" begin
        for var in (vars..., tuples..., Xs, ys)
            A = BatchView(var, 15)
            @test_throws BoundsError A[-1]
            @test_throws BoundsError A[11]
            @test @inferred(nobs(A)) == 150
            @test @inferred(length(A)) == 10
            @test @inferred(batchsize(A)) == 15
            @test @inferred(size(A)) == (10,)
            @test @inferred(getobs(A[2:3])) == getobs(BatchView(datasubset(var, 16:45), 15))
            @test @inferred(getobs(A[[1,3]])) == getobs(BatchView(datasubset(var, [1:15..., 31:45...]), 15))
            @test @inferred(A[1]) == datasubset(var, 1:15)
            @test @inferred(A[2]) == datasubset(var, 16:30)
            @test @inferred(A[3]) == datasubset(var, 31:45)
            @test A[end] == A[10]
            @test @inferred(getobs(A,1)) == getobs(var, 1:15)
            @test @inferred(getobs(A,2)) == getobs(var, 16:30)
            @test @inferred(getobs(A,3)) == getobs(var, 31:45)
            @test typeof(@inferred(collect(A))) <: Vector
        end
        for var in (vars..., tuples...)
            A = BatchView(var, 15)
            @test @inferred(getobs(A)) == A
            @test @inferred(A[2:3]) == BatchView(datasubset(var, 16:45), 15)
            @test @inferred(A[[1,3]]) == BatchView(datasubset(var, [1:15..., 31:45...]), 15)
        end
        A = BatchView(X',size=15,obsdim=1)
        @test_throws ArgumentError BatchView(X',size=15, maxsize=15, obsdim=1)
        @test A == BatchView(X',maxsize=16,obsdim=1)
        @test A == BatchView(X',count=10,obsdim=1)
        @test @inferred(length(A)) == 10
        @test @inferred(size(A)) == (10,)
        @test @inferred(A[1]) == datasubset(X', 1:15, obsdim=1)
        @test @inferred(A[2]) == datasubset(X', 16:30, obsdim=1)
        @test @inferred(A[3]) == datasubset(X', 31:45, obsdim=1)
        @test A[end] == A[10]
        @test @inferred(getobs(A,1)) == getobs(X', 1:15, obsdim=1)
        @test @inferred(getobs(A,2)) == getobs(X', 16:30, obsdim=1)
        @test @inferred(getobs(A,3)) == getobs(X', 31:45, obsdim=1)
    end

    println("<HEARTBEAT>")

    @testset "subsetting" begin
        for var in (vars..., tuples..., Xs, ys)
            A = BatchView(var)
            @test getobs(@inferred(datasubset(A))) == @inferred(getobs(A))
            @test_throws ArgumentError datasubset(A,1:5)
            S = @inferred(datasubset(A, 1:60))
            @test typeof(S) <: BatchView
            @test @inferred(nobs(S)) == 60
            @test @inferred(length(S)) == 2
            @test @inferred(size(S)) == (2,)
            @test getobs(@inferred(A[1:2])) == getobs(S)
            @test @inferred(getobs(A,1:2)) == getobs(S)
            @test @inferred(getobs(S)) == getobs(BatchView(datasubset(var,1:60),30))
            S = @inferred(DataSubset(A, 1:60))
            @test typeof(S) <: BatchView
            @test typeof(S.data) <: Union{DataSubset,Tuple}
            @test @inferred(length(S)) == 2
            @test @inferred(size(S)) == (2,)
            @test @inferred(getobs(S)) == getobs(BatchView(DataSubset(var,1:60),30))
        end
        A = BatchView(X)
        @test typeof(A.data) <: Array
        S = @inferred(datasubset(A))
        @test typeof(S) <: BatchView
        @test @inferred(length(S)) == 5
        @test typeof(S.data) <: SubArray
    end

    println("<HEARTBEAT>")

    @testset "nesting with ObsView" begin
        for var in vars
            @test eltype(@inferred(BatchView(ObsView(var)))[1]) <: Union{SubArray,String}
            @test @inferred(BatchView(BatchView(var))) == BatchView(var)
        end
        for var in tuples
            @test eltype(@inferred(BatchView(ObsView(var)))[1]) <: Tuple
            @test @inferred(BatchView(BatchView(var))) == BatchView(var)
        end
        for var in (Xs, ys)
            @test eltype(@inferred(BatchView(ObsView(var)))[1]) <: DataSubset
            @test @inferred(BatchView(BatchView(var))) == BatchView(var)
        end
        @test ObsView(BatchView(X)) == ObsView(X)
        @test ObsView(BatchView(X', obsdim=1)) == ObsView(X', obsdim=1)
        A = ObsView(X', obsdim=1)
        @test nobs(BatchView(A)) == 150
        @test size(BatchView(A)[1]) == (30,)
        @test typeof(BatchView(A)[1]) <: ObsView
        @test BatchView(A).obsdim == ObsDim.First()
    end
end


text:
using SpecialFunctions: loggamma
using ArgCheck
import Statistics
import Distributions
using StatsFuns: log2π
using ArgCheck: @argcheck
using NNlib: softplus

################################################################################
##### NormalInverseGamma
################################################################################
export NIG
export NIGRegressor
export posterior_predictive
export mean, var, logpdf
export var_aleatoric, var_epistemic, var_predict
export std_aleatoric, std_epistemic, std_predict

"""
    NIG(μ, ν, α, β)

[Normal-inverse-gamma distribution](https://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution).
"""
struct NIG{T}
    "Sample mean."
    μ::T
    "Number of pseudo observations used to estimate μ. Wikipedia uses λ instead of ν."
    ν::T
    "Half number of pseudo observations used to estimate β."
    α::T
    "Half sum of square deviations from the mean."
    β::T
    function NIG(μ, ν, α, β)
        mu, nu, alpha, beta = promote(μ, ν, α, β)
        T = typeof(nu)
        new{T}(mu, nu, alpha, beta)
    end
end
NIG(;μ, ν, α, β) = NIG(μ, ν, α, β)

function Base.show(io::IO, nig::NIG)
    print(io, "NIG(μ=$(nig.μ), ν=$(nig.ν), α=$(nig.α), β=$(nig.β))")
end

"""
NormalInverseGamma regression layer. Takes a tensor of floats of size (4c, n) and
outputs a tensor of NormalInverseGamma distributions of size (c,n).
"""
struct NIGRegressor end

function (o::NIGRegressor)(arr::AbstractMatrix)
    @argcheck mod(size(arr,1), 4) == 0
    nc,nb = size(arr)
    b = nc ÷ 4
    cmu    = view(arr, (   1):(1b),:)
    cnu    = view(arr, (1b+1):(2b),:)
    calpha = view(arr, (2b+1):(3b),:)
    cbeta  = view(arr, (3b+1):(4b),:)
    map(NIG_from_4channels, cmu, cnu, calpha, cbeta)
end

function NIG_from_4channels(cμ, cν, cα, cβ)
    μ = cμ
    ν = softplus(cν)
    α = softplus(cα) + 1
    β = softplus(cβ)
    NIG(μ,ν,α,β)
end

function to_nig_sampletype(pt)
    @argcheck length(pt) == 2
    x,σ²= pt
    return (;x,σ²)
end
function to_nig_sampletype(nt::NamedTuple)
    return (;nt.x, nt.σ²)
end

function Distributions.logpdf(nig::NIG, pt)
    s = to_nig_sampletype(pt)
    x  = s.x
    σ² = s.σ²
    μ  = nig.μ
    λ  = nig.ν
    α  = nig.α
    β  = nig.β
    l = -(2β + λ*(x-μ)^2) / (2*σ²)
    logσ² = log(σ²)
    #return √(λ)       / √(2π*σ²)          * β^α / SPF.gamma(α)       * inv(σ²)^(α+1) * exp(l)
    return 1/2*log(λ) - 1/2*(log2π +logσ²) + α*log(β) - loggamma(α) - logσ²*(α+1) + l
end

Statistics.mean(nig::NIG) = (x=nig.μ, σ²=nig.β/(nig.α - 1))
function Statistics.var(nig::NIG)
    μ = nig.μ
    λ = nig.ν
    α = nig.α
    β = nig.β
    Eσ² = mean(nig).σ²
    x  = Eσ²/ λ
    σ² =  Eσ² / (α-2)
    return (;x,σ²)
end

"""
    posterior_predictive(nig::NIG)

Return the posterior_predictive distribution, where the statistical
model is the two parameter family `Normal(μ,σ²)` and the posterior is given by `nig`.
"""
function posterior_predictive(nig::NIG)
    # https://en.wikipedia.org/wiki/Conjugate_prior
    μ = nig.μ # mean
    α = nig.α # 2α = nobs var
    ν = nig.ν # nobs mean
    β = nig.β # 2β = sum square dev
    var = β * (ν+1) / (ν*α)  # taken from Wikipedia
    studentνμσ(2α, μ, sqrt(var))
end

"""
    studentνμσ(ν,μ,σ)

Create a student t distributions with
* ν degrees of freedom
* μ location (= mean)
* σ scale (≈ std deviation)
"""
function studentνμσ(ν,μ,σ)
    d0 = Distributions.TDist(ν)
    Distributions.LocationScale(μ,σ,d0)
end

function evidence(o::NIG)
    2*o.ν + o.α
end
function predict(o::NIG)
    # This is E[μ]
    o.μ
end

"""
    var_aleatoric(d)

[Alleatoric](https://en.wikipedia.org/wiki/Uncertainty_quantification#Aleatoric_and_epistemic_uncertainty) variance of `predict(d)`. Also known as statistical uncertainty.

See also [var_epistemic](@ref).
"""
function var_aleatoric end
function var_aleatoric(o::NIG)
    # E(σ²)
    return o.β / (o.α - 1)
end

"""
    std_aleatoric(o) = √(var_aleatoric(o))
"""
std_aleatoric(o) = √(var_aleatoric(o))

"""
    var_epistemic(d)

[Epistemic](https://en.wikipedia.org/wiki/Uncertainty_quantification#Aleatoric_and_epistemic_uncertainty) variance of `predict(d)`. Also known as systematic uncertainty.

See also [var_aleatoric](@ref).
"""
function var_epistemic end

function var_epistemic(o::NIG)
    # Var(μ)
    o.β / ((o.α - 1) * o.ν)
end

"""
    std_epistemic(o) = √(var_epistemic(o))
"""
std_epistemic(o) = √(var_epistemic(o))

std_predict(o::NIG) = Statistics.std(posterior_predictive(o))
var_predict(o::NIG) = var(posterior_predictive(o))

ext: jl


text:
module Stats

# Backend module for statistics functionality

using Statistics, Distributions, DataFrames

export

    runaboutmedian, eaboutmedian, runupdown, eupdown,

    cluster, mixture, trend, oscillation,

    cpk, ppk

# ------------------------------------------------------------------------------------------------
# Run chart analysis
function runaboutmedian(points::Vector{Real})
    n = 1
    center = median(points)
    for (i,v) in enumerate(points[2:end])
       (v - center) * (points[i] - center) < 0 && (n += 1)
    end
    n
 end

function eaboutmedian(points::Vector{Real})
    N = length(points)
    m = length(findall(x->x > median(points), points))
    n = N - m
    1 + 2 * m * n / N
end

function cluster(points::Vector{Real})
    N = length(points)
    R = runaboutmedian(points)
    E = eaboutmedian(points)
    Z = (R - E) / sqrt((E - 1) * (E - 2) / (N - 1))
    cdf(Normal(0, 1), Z)
end

mixture(data::Vector{Real}) = 1 - cluster(data)

function runupdown(points::Vector{Real})
    n = 1
    df = diff(points)
    prod = df[1] > 0
    for i in df[2:end]
        i  != 0 && (prod * i < 0) && (prod *= -1 ; n += 1)
    end
    n
end

eupdown(points::Vector{Real}) = 2 / 3 * (2 * length(points) - 1)

function trend(points::Vector{Real})
    N = length(points)
    v = runupdown(points)
    Z = (v - (2 * N - 1) / 3) / sqrt((16 * N  - 29) / 90)
    cdf(Normal(0, 1), Z)
end

oscillation(points::Vector{Real}) = 1 - trend(points)

# ---------------------------------------------------------------------------------------------------------
# Process capability/performance index
# cp/pp/z.bench/z.usl/z.lsl/ppm

function cpk(x::Vector{Real}; s::Real = 3, USL::Real = 110, LSL::Real = 90, correct::Bool = true)
    d = length(x) - 1
    μ = mean(x)
    correct ? (sd = std(x) / (sqrt(2 / d) * gamma(d / 2 + 0.5) / gamma(d / 2))) : (sd = std(x))
    min((USL - μ) / s / sd, (μ - LSL) / s / sd)
 end

function cpk(x::Matrix{Real}; s::Real = 3 USL::Real = 110, LSL::Real = 90, correct::Bool = true)
    sd = map(1:size(x, 2)) do i
       std(x[:,i])
    end
    d = size(x, 1) * size(x, 2) - size(x, 2)
    μ = mean(x)
    sd = sum(sd .* (size(x, 1) - 1)) / d
    correct && (sd /= (sqrt(2 / d) * gamma(d / 2 + 0.5) / gamma(d / 2)))
    min((USL - μ) / s / sd, (μ - LSL) / s  / sd)
end

cpk(x::DataFrame; s::Real = 3 USL::Real = 110, LSL::Real = 90, correct::Bool = true) =
    cpk(Matrix{Real}(x); s = s, USL = USL, LSL = LSL, correct = correct)

function ppk(x::Matrix{Real}; s::Real = 3, USL::Real = 110, LSL::Real = 90, correct::Bool = false)
    d = size(x, 1) * size(x, 2) - 1
    μ = mean(x)
    correct ? (sd = std(x) / (sqrt(2 / d) * gamma(d / 2 + 0.5) / gamma(d / 2))) : (sd = std(x))
    min((USL - μ) / s / sd, (μ - LSL) / s / sd)
end

ppk(x::DataFrame; s::Real, USL::Real = 110, LSL::Real = 90, correct::Bool = false) =
    ppk(Matrix{Real}(x), s = s, USL = USL, LSL = LSL, correct = correct)

end


text:
using Documenter
using Twinkle, Twinkle.FlexUI

makedocs(
    sitename = "Twinkle.jl",
    format = Documenter.HTML(
        # prettyurls = get(ENV, "CI", nothing) == "true",
        assets = [asset("assets/TwinkleJulia.png", class = :ico, islocal = true)],
    ),
    modules = [Twinkle],
    pages = [
        "Home" => "index.md",
    ],
    expandfirst = [])

println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
println("@ SANITY CHECK")
println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")

# deploydocs(
#     repo = "github.com/galran/Twinkle.jl.git",
#     devbranch = "main",
#     push_preview = true,
#     target = "build",
#     deps = nothing,
#     make = nothing,
# )

deploydocs(
    repo = "github.com/galran/Twinkle.jl.git",
    branch = "gh-pages",
    devbranch = "main",
    push_preview = true,
)


text:
#=
    Interpolate using cubic Hermite splines. The breakpoints in arrays xbp and ybp are assumed to be sorted.
    Evaluate the function in all points of the array xeval.
    Methods:
        "Linear"                yuck
        "FiniteDifference"      classic cubic interpolation, no tension parameter
                                Finite difference can overshoot for non-monotonic data
        "Cardinal"              cubic cardinal splines, uses tension parameter which must be between [0,1]
                                cubin cardinal splines can overshoot for non-monotonic data
                                (increasing tension decreases overshoot)
        "FritschCarlson"        monotonic - tangents are first initialized, then adjusted if they are not monotonic
                                can overshoot for non-monotonic data
        "FritschButland"        monotonic - faster algorithm (only requires one pass) but somewhat higher apparent "tension"
        "Steffen"               monotonic - also only one pass, results usually between FritschCarlson and FritschButland
    Sources:
        Fritsch & Carlson (1980), "Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.
        Fritsch & Butland (1984), "A Method for Constructing Local Monotone Piecewise Cubic Interpolants", doi:10.1137/0905021.
        Steffen (1990), "A Simple Method for Monotonic Interpolation in One Dimension", http://adsabs.harvard.edu/abs/1990A%26A...239..443S

    Implementation based on http://bl.ocks.org/niclasmattsson/7bceb05fba6c71c78d507adae3d29417
=#

export
    LinearMonotonicInterpolation,
    FiniteDifferenceMonotonicInterpolation,
    CardinalMonotonicInterpolation,
    FritschCarlsonMonotonicInterpolation,
    FritschButlandMonotonicInterpolation,
    SteffenMonotonicInterpolation

"""
    MonotonicInterpolationType

Abstract class for all types of monotonic interpolation.
"""
abstract type MonotonicInterpolationType <: InterpolationType end

"""
    LinearMonotonicInterpolation

Simple linear interpolation.
"""
struct LinearMonotonicInterpolation <: MonotonicInterpolationType
end

"""
    FiniteDifferenceMonotonicInterpolation

Classic cubic interpolation, no tension parameter.
Finite difference can overshoot for non-monotonic data.
"""
struct FiniteDifferenceMonotonicInterpolation <: MonotonicInterpolationType
end

"""
    CardinalMonotonicInterpolation(tension)

Cubic cardinal splines, uses `tension` parameter which must be between [0,1]
Cubin cardinal splines can overshoot for non-monotonic data
(increasing tension reduces overshoot).
"""
struct CardinalMonotonicInterpolation{TTension<:Number} <: MonotonicInterpolationType
    tension :: TTension # must be in [0, 1]
end

"""
    FritschCarlsonMonotonicInterpolation

Monotonic interpolation based on Fritsch & Carlson (1980),
"Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.

Tangents are first initialized, then adjusted if they are not monotonic
can overshoot for non-monotonic data
"""
struct FritschCarlsonMonotonicInterpolation <: MonotonicInterpolationType
end

"""
    FritschButlandMonotonicInterpolation

Monotonic interpolation based on  Fritsch & Butland (1984),
"A Method for Constructing Local Monotone Piecewise Cubic Interpolants",
doi:10.1137/0905021.

Faster than FritschCarlsonMonotonicInterpolation (only requires one pass)
but somewhat higher apparent "tension".
"""
struct FritschButlandMonotonicInterpolation <: MonotonicInterpolationType
end

"""
    SteffenMonotonicInterpolation

Monotonic interpolation based on Steffen (1990),
"A Simple Method for Monotonic Interpolation in One Dimension",
http://adsabs.harvard.edu/abs/1990A%26A...239..443S

Only one pass, results usually between FritschCarlson and FritschButland.
"""
struct SteffenMonotonicInterpolation <: MonotonicInterpolationType
end

"""
    MonotonicInterpolation

Monotonic interpolation up to third order represented by type, knots and
coefficients.
"""
struct MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType,
    TKnots<:AbstractVector{<:Number}, TACoeff <: AbstractArray{TEl,1}} <: AbstractInterpolation{T,1,DimSpec{TInterpolationType}}

    it::TInterpolationType
    knots::TKnots
    A::TACoeff # constant parts of piecewise polynomials
    m::Vector{TCoeffs} # coefficients of linear parts of piecewise polynomials
    c::Vector{TCoeffs} # coefficients of quadratic parts of piecewise polynomials
    d::Vector{TCoeffs} # coefficients of cubic parts of piecewise polynomials
end

size(A::MonotonicInterpolation) = size(A.knots)
axes(A::MonotonicInterpolation) = axes(A.knots)

function MonotonicInterpolation(::Type{TWeights}, it::TInterpolationType, knots::TKnots, A::AbstractArray{TEl,1},
    m::Vector{TCoeffs}, c::Vector{TCoeffs}, d::Vector{TCoeffs}) where {TWeights, TCoeffs, TEl, TInterpolationType<:MonotonicInterpolationType, TKnots<:AbstractVector{<:Number}}

    isconcretetype(TInterpolationType) || error("The spline type must be a leaf type (was $TInterpolationType)")
    isconcretetype(TCoeffs) || warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")

    check_monotonic(knots, A)

    cZero = zero(TWeights)
    if isempty(A)
        T = Base.promote_op(*, typeof(cZero), eltype(A))
    else
        T = typeof(cZero * first(A))
    end

    MonotonicInterpolation{T, TCoeffs, TEl, TInterpolationType, TKnots, typeof(A)}(it, knots, A, m, c, d)
end

function interpolate(::Type{TWeights}, ::Type{TCoeffs}, knots::TKnots,
    A::AbstractArray{TEl,1}, it::TInterpolationType) where {TWeights,TCoeffs,TEl,TKnots<:AbstractVector{<:Number},TInterpolationType<:MonotonicInterpolationType}

    check_monotonic(knots, A)

    # first we need to determine tangents (m)
    n = length(knots)
    m, Δ = calcTangents(TCoeffs, knots, A, it)
    c = Vector{TCoeffs}(undef, n-1)
    d = Vector{TCoeffs}(undef, n-1)
    for k ∈ 1:n-1
        if TInterpolationType == LinearMonotonicInterpolation
            c[k] = d[k] = zero(TCoeffs)
        else
            xdiff = knots[k+1] - knots[k]
            c[k] = (3*Δ[k] - 2*m[k] - m[k+1]) / xdiff
            d[k] = (m[k] + m[k+1] - 2*Δ[k]) / (xdiff * xdiff)
        end
    end

    MonotonicInterpolation(TWeights, it, knots, A, m, c, d)
end

function interpolate(knots::AbstractVector{<:Number}, A::AbstractArray{TEl,1},
    it::TInterpolationType) where {TEl,TInterpolationType<:MonotonicInterpolationType}

    interpolate(tweight(A), tcoef(A), knots, A, it)
end

function (itp::MonotonicInterpolation)(x::Number)
    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
    k = searchsortedfirst(itp.knots, x)
    if k > 1
        k -= 1
    end
    xdiff = x - itp.knots[k]
    return itp.A[k] + itp.m[k]*xdiff + itp.c[k]*xdiff*xdiff + itp.d[k]*xdiff*xdiff*xdiff
end

function gradient(itp::MonotonicInterpolation, x::Number)
    return SVector(gradient1(itp, x))
end

function gradient1(itp::MonotonicInterpolation, x::Number)
    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
    k = searchsortedfirst(itp.knots, x)
    if k > 1
        k -= 1
    end
    xdiff = x - itp.knots[k]
    return itp.m[k] + 2*itp.c[k]*xdiff + 3*itp.d[k]*xdiff*xdiff
end

function hessian(itp::MonotonicInterpolation, x::Number)
    return SVector(hessian1(itp, x))
end

function hessian1(itp::MonotonicInterpolation, x::Number)
    @boundscheck (checkbounds(Bool, itp, x) || Base.throw_boundserror(itp, (x,)))
    k = searchsortedfirst(itp.knots, x)
    if k > 1
        k -= 1
    end
    xdiff = x - itp.knots[k]
    return 2*itp.c[k] + 6*itp.d[k]*xdiff
end

@inline function check_monotonic(knots, A)
    axes(knots) == axes(A) || throw(DimensionMismatch("knot vector must have the same axes as the corresponding array"))
    issorted(knots) || error("knot-vector must be sorted in increasing order")
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
    y::AbstractVector{TEl}, method::LinearMonotonicInterpolation) where {TCoeffs, TEl}

    n = length(x)
    Δ = Vector{TCoeffs}(undef, n-1)
    m = Vector{TCoeffs}(undef, n)
    for k in 1:n-1
        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
        Δ[k] = Δk
        m[k] = Δk
    end
    m[n] = Δ[n-1]
    return (m, Δ)
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
    y::AbstractVector{TEl}, method::FiniteDifferenceMonotonicInterpolation) where {TCoeffs, TEl}

    n = length(x)
    Δ = Vector{TCoeffs}(undef, n-1)
    m = Vector{TCoeffs}(undef, n)
    for k in 1:n-1
        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
        Δ[k] = Δk
        if k == 1   # left endpoint
            m[k] = Δk
        else
            m[k] = (Δ[k-1] + Δk) / 2
        end
    end
    m[n] = Δ[n-1]
    return (m, Δ)
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
    y::AbstractVector{TEl}, method::CardinalMonotonicInterpolation{TTension}) where {TTension, TCoeffs, TEl}

    n = length(x)
    Δ = Vector{TCoeffs}(undef, n-1)
    m = Vector{TCoeffs}(undef, n)
    for k in 1:n-1
        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
        Δ[k] = Δk
        if k == 1   # left endpoint
            m[k] = Δk
        else
            m[k] = (oneunit(TTension) - method.tension) * (y[k+1] - y[k-1]) / (x[k+1] - x[k-1])
        end
    end
    m[n] = Δ[n-1]
    return (m, Δ)
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
    y::AbstractVector{TEl}, method::FritschCarlsonMonotonicInterpolation) where {TCoeffs, TEl}

    n = length(x)
    Δ = Vector{TCoeffs}(undef, n-1)
    m = Vector{TCoeffs}(undef, n)
    for k in 1:n-1
        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
        Δ[k] = Δk
        if k == 1   # left endpoint
            m[k] = Δk
        else
            # If any consecutive secant lines change sign (i.e. curve changes direction), initialize the tangent to zero.
            # This is needed to make the interpolation monotonic. Otherwise set tangent to the average of the secants.
            if Δk <= zero(Δk)
                m[k] = zero(TCoeffs)
            else
                m[k] = Δ[k-1] * (Δ[k-1] + Δk) / 2.0
            end
        end
    end
    m[n] = Δ[n-1]
    #=
    Fritsch & Carlson derived necessary and sufficient conditions for monotonicity in their 1980 paper. Splines will be
    monotonic if all tangents are in a certain region of the alpha-beta plane, with alpha and beta as defined below.
    A robust choice is to put alpha & beta within a circle around origo with radius 3. The FritschCarlson algorithm
    makes simple initial estimates of tangents and then does another pass over data points to move any outlier tangents
    into the monotonic region. FritschButland & Steffen algorithms make more elaborate first estimates of tangents that
    are guaranteed to lie in the monotonic region, so no second pass is necessary.
    =#

    # Second pass of FritschCarlson: adjust any non-monotonic tangents.
    for k in 1:n-1
        Δk = Δ[k]
        if Δk == zero(TCoeffs)
            m[k] = zero(TCoeffs)
            m[k+1] = zero(TCoeffs)
            continue
        end
        α = m[k] / Δk
        β = m[k+1] / Δk
        τ = 3.0 / sqrt(α^2 + β^2)
        if τ < 1.0 # if we're outside the circle with radius 3 then move onto the circle
            m[k] = τ * α * Δk
            m[k+1] = τ * β * Δk
        end
    end
    return (m, Δ)
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
    y :: AbstractVector{TEl}, method :: FritschButlandMonotonicInterpolation) where {TCoeffs, TEl}

    # based on Fritsch & Butland (1984),
    # "A Method for Constructing Local Monotone Piecewise Cubic Interpolants",
    # doi:10.1137/0905021.

    n = length(x)
    Δ = Vector{TCoeffs}(undef, n-1)
    m = Vector{TCoeffs}(undef, n)
    for k in 1:n-1
        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
        Δ[k] = Δk
        if k == 1   # left endpoint
            m[k] = Δk
        elseif Δ[k-1] * Δk <= zero(TCoeffs)
            m[k] = zero(TCoeffs)
        else
            α = (1.0 + (x[k+1] - x[k]) / (x[k+1] - x[k-1])) / 3.0
            m[k] = Δ[k-1] * Δk / (α*Δk + (1.0 - α)*Δ[k-1])
        end
    end
    m[n] = Δ[n-1]
    return (m, Δ)
end

function calcTangents(::Type{TCoeffs}, x::AbstractVector{<:Number},
    y::AbstractVector{TEl}, method::SteffenMonotonicInterpolation) where {TCoeffs, TEl}

    # Steffen (1990),
    # "A Simple Method for Monotonic Interpolation in One Dimension",
    # http://adsabs.harvard.edu/abs/1990A%26A...239..443S

    n = length(x)
    Δ = Vector{TCoeffs}(undef, n-1)
    m = Vector{TCoeffs}(undef, n)
    for k in 1:n-1
        Δk = (y[k+1] - y[k]) / (x[k+1] - x[k])
        Δ[k] = Δk
        if k == 1   # left endpoint
            m[k] = Δk
        else
            p = ((x[k+1] - x[k]) * Δ[k-1] + (x[k] - x[k-1]) * Δk) / (x[k+1] - x[k-1])
            m[k] = (sign(Δ[k-1]) + sign(Δk)) *
                min(abs(Δ[k-1]), abs(Δk), 0.5*abs(p))
        end
    end
    m[n] = Δ[n-1]
    return (m, Δ)
end

# How many non-NoInterp dimensions are there?
count_interp_dims(::Type{<:MonotonicInterpolationType}) = 1

lbounds(itp::MonotonicInterpolation) = (first(itp.knots),)
ubounds(itp::MonotonicInterpolation) = (last(itp.knots),)


text:
module PiezoSerial

import Base.position

using SerialPorts
using ImagineHardware
using Sockets

export find_piezo_serial,
        kp, ki, kd, set_kp, set_ki, set_kd,
        status,
        runtime_minutes,
        accepts_mod, mod_on, mod_off,
        monsrc, set_monsrc,
        is_fan_on, is_fan_off, set_fan_on, set_fan_off,
        position, set_position,
        is_open_loop, is_closed_loop, set_closed_loop, set_open_loop,
        is_notch_on, is_notch_off, set_notch_on, set_notch_off, notchf, set_notchf, notchb, set_notchb,
        is_lowpass_on, is_lowpass_off, lowpassf, set_lowpassf,
        slewrate, set_slewrate

ENTER = "\r\n" #on Windows

SUPPORTED_AMPLIFIERS = ["30DV300";]

has_serial_control(amp::SimpleAmplifier) = has_serial(name(amp))
has_serial_control(ampname::String) = in(ampname, SUPPORTED_AMPLIFIERS)

#SerialPorts.jl just uses PySerial. Here's the constructor from the PySerial docs:
#__init__(port=None, baudrate=9600, bytesize=EIGHTBITS, parity=PARITY_NONE, stopbits=STOPBITS_ONE, timeout=None, xonxoff=False, rtscts=False, write_timeout=None, dsrdtr=False, inter_byte_timeout=None)
#30DV300 amp serial port properties: 115200 baud, 8 bit, no parity, 1 stop bit, software handshake (XON / XOFF)
#...so the only extra parameters that we need to set different from the default are baudrate and xonxoff (set to true)

function piezo_serialport(portname::String)
    #modified from SerialPorts source code:
    py_ptr = SerialPorts.PySerial[:Serial](portname, 115200, xonxoff=true)
    sp = SerialPort(portname,
               115200,
               py_ptr[:bytesize],
               py_ptr[:parity],
               py_ptr[:stopbits],
               py_ptr[:timeout],
               py_ptr[:xonxoff],
               py_ptr[:rtscts],
               py_ptr[:dsrdtr], py_ptr)
    atexit(()->close(sp))
    return sp
end

#try communication
function has_piezo(p::SerialPort)
    try
        rmin = runtime_minutes(p)
        return true
    catch
        return false
    end
end

#Iterate through serial ports and return an open SerialPort to the first found with a piezo amp attached
function find_piezo_serial()
    allports = list_serialports()
    for p in allports
	prt = piezo_serialport(p)
        if has_piezo(prt)
            return prt
        else
            close(prt)
        end
    end
    error("No supported serial-capable piezo was found")
end

function wait_available(p::SerialPort, timeout = 0.8) #timeout in seconds
    t = time()
    while time() - t < timeout
        if nb_available(p) > 0
            return true
        end
        sleep(0.1)
    end
    error("Timed out")
end

function Sockets.listen(p::SerialPort, timeout = 1.0; prefix = "") #timeout in seconds
	resp = listen_raw(p, timeout)
	temp =  split(resp, ENTER)
	if length(temp) !=2
		error("Expected only one press of the ENTER key in the response string.  Retry the command or try listen_raw.")
	end
	resp =  temp[1] #remove trailing newline
	if !isempty(prefix)
		if resp[1:length(prefix)] != prefix
			error("The expected prefix wasn't returned.  This was returned: $resp")
		end
		resp = resp[length(prefix)+2:end] #one extra index for the comma
	end
	return resp
end
function listen_raw(p::SerialPort, timeout= 1.0)
    wait_available(p, timeout)
    resp =  readavailable(p)
end

function status_register(piezoport::SerialPort)
	comstr = "stat"
    write(piezoport, comstr * ENTER)
    dig_str = listen(piezoport; prefix = comstr)
	return reverse(bits(parse(UInt16, dig_str))) #bstr elements 1 to 16 correspond with 0 to 15 in the 30DV300 manual
end

function status(piezoport::SerialPort)
	bstr = status_register(piezoport)
    bstr[1] == '0' ? print("Piezo is unplugged\n") : print("Piezo is plugged in\n")
    bstr[2] == '0' && bstr[3] == '0' ? print("Piezo has no measurement system\n") : begin
        bstr[2] == '1' && bstr[3] == '0' ? print("Piezo has strain gauge measurement system\n") : print("Piezo has capacitive measurement system\n") end
    bstr[4] == '0' ? print("Piezo system is closed-loop capable\n") : print("Piezo system is open-loop only\n")
    bstr[7] == '0' ? print("Piezo voltage is not enabled\n") : print("Piezo voltage is enabled\n")
    bstr[8] == '0' ? print("Piezo is set to open-loop mode\n") : print("Piezo is set to closed-loop mode\n")
    all(collect(bstr[10:12]).== '0') ? print("Built-in function generator is disabled\n") : print("Built-in function generator is enabled\n")
    bstr[13] == '0' ? print("Input notch filter is off\n") : print("Input notch filter is on\n")
    bstr[14] == '0' ? print("Input lowpass filter is off\n") : print("Input lowpass filter is on\n")
    bstr[16] == '0' ? print("Fan is off\n") : print("Fan is on\n")
    return bstr
end

function available_commands(piezoport::SerialPort)
    write(piezoport, "s" * ENTER)
	@show resp = listen_raw(piezoport)
    temp = split(resp, ['\n', '\r', ' '])
	return temp[find(x->length(x) > 0, temp)]
end

function runtime_minutes(piezoport::SerialPort)
    comstr = "rohm"
    write(piezoport, comstr * ENTER)
    return parse(Int, listen(piezoport; prefix = comstr))
end

function accepts_mod(piezoport::SerialPort)
    write(piezoport, "modon" * ENTER)
    return parse(Bool, first(listen(piezoport; prefix="modon")))
end
mod_on(piezoport::SerialPort) = write(piezoport, "modon,1" * ENTER)
mod_off(piezoport::SerialPort) = write(piezoport, "modon,0" * ENTER)

is_fan_off(piezoport::SerialPort) = status_register(piezoport)[16] == '0'
is_fan_on(piezoport::SerialPort) = !is_fan_off(piezoport)

function set_fan_on(piezoport::SerialPort)
    if is_fan_on(piezoport)
        warn("Fan is already on")
    end
    write(piezoport, "fan,1" * ENTER)
    sleep(0.1)
    if !is_fan_on(piezoport)
        error("Failed to turn on fan")
    end
end

function set_fan_off(piezoport::SerialPort)
    if is_fan_off(piezoport)
        warn("Fan is already off")
    end
    write(piezoport, "fan,0" * ENTER)
    sleep(0.1)
    if !is_fan_off(piezoport)
        error("Failed to turn off fan")
    end
end

function monsrc(piezoport::SerialPort)
    write(piezoport, "monsrc" * ENTER)
    return parse(Int, listen(piezoport; prefix="monsrc"))
end

#0 = position in closed loop (default).  6 = position in open loop (requires rescale). See manual for the rest.
function set_monsrc(piezoport::SerialPort, src_id::Int)
    @assert src_id >= 0 && src_id <= 6
    write(piezoport, "monsrc,$src_id " * ENTER)
end

function position(piezoport::SerialPort)
    write(piezoport, "mess" * ENTER)
    return listen(piezoport; prefix="mess")
end

#Note that tp should be in microns when in closed-loop mode and in volts (-20.0 to 130.0) when in open-loop mode
function set_position(piezoport::SerialPort, p::Float64)
    tp = trunc(p, 3)
	@assert tp >= 0.0 && tp <= 800.0
    write(piezoport, "set,$tp " * ENTER)
end

#kp, ki, and kd much each lie in the range 0 to 999 according to the manual.  Yet our kd was set to -0.3 when shipped(?)
function set_kp(piezoport::SerialPort, val::Float64)
    tval = trunc(val, 3)
    write(piezoport, "kp,$(tval)" * ENTER)
end
#Getting isn't mentioned in the manual, but this works to get rather than set the value
function kp(piezoport::SerialPort)
    comstr = "kp"
    write(piezoport, comstr * ENTER)
    return parse(Float64, listen(piezoport; prefix=comstr))
end
function set_ki(piezoport::SerialPort, val::Float64)
    @assert val >=0.0 && val <= 999.0
    @show tval = trunc(val, 3)
    write(piezoport, "ki,$(tval)" * ENTER)
end
function ki(piezoport::SerialPort)
    comstr = "ki"
    write(piezoport, comstr * ENTER)
    return parse(Float64, listen(piezoport; prefix=comstr))
end
function set_kd(piezoport::SerialPort, val::Float64)
    @assert val >=0.0 && val <= 999.0
    @show tval = trunc(val, 3)
    write(piezoport, "kd,$(tval)" * ENTER)
end
function kd(piezoport::SerialPort)
    comstr = "kd"
    write(piezoport, comstr * ENTER)
    return parse(Float64, listen(piezoport; prefix=comstr))
end

is_open_loop(piezoport::SerialPort) = status_register(piezoport)[8] == '0'
is_closed_loop(piezoport::SerialPort) = !is_open_loop(piezoport)
function set_open_loop(piezoport::SerialPort)
    write(piezoport, "cl,0" * ENTER)
    sleep(0.1)
    if !is_open_loop(piezoport)
        error("Command did not succeed.")
    end
end
function set_closed_loop(piezoport::SerialPort)
    write(piezoport, "cl,1" * ENTER)
    sleep(0.1)
    if !is_closed_loop(piezoport)
        error("Command did not succeed.")
    end
end

#NOTE: was on when shipped
is_notch_off(piezoport::SerialPort) = status_register(piezoport)[13] == '0'
is_notch_on(piezoport::SerialPort) = !is_notch_off(piezoport)

function set_notch_off(piezoport::SerialPort)
    write(piezoport, "notchon,0" * ENTER)
    sleep(0.1)
    if !is_notch_off(piezoport)
        error("Command did not succeed.")
    end
end

function set_notch_on(piezoport::SerialPort)
    write(piezoport, "notchon,1" * ENTER)
    sleep(0.1)
    if !is_notch_on(piezoport)
        error("Command did not succeed.")
    end
end

#NOTE: the value when shipped was 200
function notchf(piezoport::SerialPort)
    comstr = "notchf"
    write(piezoport, comstr * ENTER)
    return parse(Int, listen(piezoport; prefix=comstr))
end
function set_notchf(piezoport::SerialPort, freq::Int) #Hz units
    @assert freq>=0 && freq <= 20000
    write(piezoport, "notchf,$(freq)" * ENTER)
    sleep(0.1)
    f = notchf(piezoport)
    if f != freq
        error("Notch frequency setting failed.  Current value is $f")
    end
end

#NOTE: the value when shipped was 100
function notchb(piezoport::SerialPort)
    comstr = "notchb"
    write(piezoport, comstr * ENTER)
    return parse(Int, listen(piezoport; prefix=comstr))
end
function set_notchb(piezoport::SerialPort, bw::Int) #Hz units
    @assert bw >=0 && bw <= 2 * notchf(piezoport)
    write(piezoport, "notchb,$(bw)" * ENTER)
    sleep(0.1)
    f = notchb(piezoport)
    if f != bw
        error("Notch bandwidth setting failed.  Current value is $f")
    end
end

#sr set/get slewrate 0.0000002 to 500.0 V/s (applied to MOD)
function slewrate(piezoport::SerialPort)
    comstr = "sr"
    write(piezoport, comstr * ENTER)
    return parse(Float64, listen(piezoport; prefix=comstr))
end
#Note: when shipped by the value was 40.0 V/s, seems a bit low?
function set_slewrate(piezoport::SerialPort, sr::Float64) #V/s units
    @assert sr >=0.0000002 && sr<=500.0
    write(piezoport, "sr,$(sr)" * ENTER)
    sleep(0.1)
    s = slewrate(piezoport)
    if s != sr
        error("Slewrate setting failed.  Current value is $s")
    end
end

#Note: this was off when shipped from the company
is_lowpass_off(piezoport::SerialPort) = status_register(piezoport)[14] == '0'
is_lowpass_on(piezoport::SerialPort) = !is_lowpass_off(piezoport)

#lpon set/get (1 on, 0 off)
function set_lowpass_on(piezoport::SerialPort)
    if is_lowpass_on(piezoport)
        warn("Lowpass filter is already on")
    end
    write(piezoport, "lpon,1" * ENTER)
    sleep(0.1)
    if !is_lowpass_on(piezoport)
        error("Failed to turn on lowpass filter")
    end
end

function set_lowpass_off(piezoport::SerialPort)
    if is_lowpass_off(piezoport)
        warn("Lowpass filter is already off")
    end
    write(piezoport, "lpon,0" * ENTER)
    sleep(0.1)
    if !is_lowpass_off(piezoport)
        error("Failed to turn off lowpass filter")
    end
end

#lpf set/get 0 to 20000
#NOTE: the value when shipped was 45 (but it was turned off)
function lowpassf(piezoport::SerialPort)
    comstr = "lpf"
    write(piezoport, comstr * ENTER)
    return parse(Int, listen(piezoport; prefix=comstr))
end
function set_lowpassf(piezoport::SerialPort, fnew::Int) #Hz units
    @assert fnew >=0 && fnew <= 20000
    write(piezoport, "lpf,$(fnew)" * ENTER)
    sleep(0.1)
    f = lowpassf(piezoport)
    if f != fnew
        error("Lowpass frequency setting failed.  Current value is $f")
    end
end

########################
#BELOW IS NOT YET TESTED
########################

#sets default values for kp, ki, kd stored in the device
set_pid_defaults(piezoport::SerialPort) = write(piezoport, "sstd" * ENTER)
#Set it back to the way it was shipped to us. (We asked for custom calibration)
#This may or may not be different from the defaults stored in the controller
function set_pid_factory(piezoport::SerialPort)
    set_kp(piezoport, -0.3)
    set_ki(piezoport, 50.0)
    set_kd(piezoport, 0.1)
end

#TODO: support the function generator commands

end #module


text:
## Conditional Expectations E[X̄_{l}|Y;ϕ] and E[X̄_{l}X̄_{l+1}|Y;ϕ]
using CpelNano
using LinearAlgebra
using BenchmarkTools

# Parameters
L = 300; αs = fill(-1.0, L); βs = fill(0.0, L - 1);

# Data
obs = [CpelNano.MethCallCpgGrp(-20.0, 0.0) for _ in 1:L];

# Get u's
u1 = CpelNano.get_uc(αs[1], obs[1]); uN = CpelNano.get_uc(αs[end], obs[end]);
logu1 = CpelNano.get_log_uc(αs[1], obs[1]); loguL = CpelNano.get_log_uc(αs[end], obs[end]);

# Get W's
Ws = [CpelNano.get_Wc(αs[l], αs[l + 1], βs[l], obs[l], obs[l + 1]) for l = 1:(L - 1)]
logWs = [CpelNano.get_log_Wc(αs[l], αs[l + 1], βs[l], obs[l], obs[l + 1]) for l = 1:(L - 1)]

# Scaling
Ws = UniformScaling(1.0 / maximum(maximum.(Ws))) * Ws

# Get Zc
Zc = CpelNano.get_Zc(u1, uN, Ws)
logZc = CpelNano.get_log_Zc(logu1, loguL, logWs)

# Ex
CpelNano.get_Ec_X(u1,uN,Ws,Zc)
CpelNano.get_Ec_X_log(logu1,loguL,logWs,logZc)

# Exx
CpelNano.get_Ec_XX(u1,uN,Ws,Zc)
CpelNano.get_Ec_XX_log(logu1,loguL,logWs,logZc)

# E[logp(y|x)|y]
CpelNano.get_Ec_logpyx(u1,uN,Ws,Zc,obs)
CpelNano.get_Ec_logpyx_log(logu1,loguL,logWs,logZc,obs)

## Expectations E[X̄_{l};ϕ] and E[X̄_{l}X̄_{l+1};ϕ]
using CpelNano
using LinearAlgebra
using BenchmarkTools

# Parameters
L = 300; αs = fill(1.0, L); βs = fill(1.0, L - 1);

# Get u's
u1 = CpelNano.get_u(αs[1]); uL = CpelNano.get_u(αs[end]);
logu1 = CpelNano.get_log_u(αs[1]); loguL = CpelNano.get_log_u(αs[end]);

# Get W's
Ws = [CpelNano.get_W(αs[l], αs[l + 1], βs[l]) for l = 1:length(βs)];
logWs = [CpelNano.get_log_W(αs[l], αs[l + 1], βs[l]) for l = 1:length(βs)];

# Scaling
Ws = UniformScaling(1.0 / maximum(maximum.(Ws))) * Ws

# Get Z
Z = CpelNano.get_Z(u1, uL, Ws)
logZ = CpelNano.get_log_Z(logu1, loguL, logWs)

# E[X]
CpelNano.get_E_X(u1,uL,Ws,Z)
CpelNano.get_E_X_log(logu1,loguL,logWs,logZ)

# E[XX]
CpelNano.get_E_XX(u1,uL,Ws,Z)
CpelNano.get_E_XX_log(logu1,loguL,logWs,logZ)

# Marginal p(x)
CpelNano.get_marg_px(u1,uL,Ws,Z)
CpelNano.get_marg_px_log(logu1,loguL,logWs,logZ)


text:
struct GeoTable{T<:Array}
    features::T
end

function GeoTable(features::T) where T<:Array
    geom = features[1].geometry
    if(geom === nothing)
        error("Not a valid geojson")
    end
    GeoTable{T}(features)
end

Tables.istable(::Type{<:GeoTable}) = true
Tables.rowaccess(::Type{<:GeoTable}) = true
Tables.rows(gt::GeoTable) = gt

Base.IteratorSize(::Type{<:GeoTable}) = Base.HasLength()
Base.length(gt::GeoTable) = (length(gt.features))
Base.IteratorEltype(::Type{<:GeoTable}) = Base.HasEltype()

# read only AbstractVector
Base.size(gt::GeoTable) = length(gt.features)
Base.getindex(gt::GeoTable, i) = iterate((gt.features), i)
Base.IndexStyle(::Type{<:GeoTable}) = IndexLinear()

struct GeoTableRow
    UUID   # TODO : fetch property names paticular to the geojson
    FID
    Geometry
end

function Tables.schema(features::Array)

    names_fields = (keys(features[1].properties))
    geom_names = ("Geometry")
    names = (names_fields..., geom_names)

    #todo types
    # types_fields = Tuple(typeofvalu(features[1].properties))
    # geom_types = #todo
    # types = (types_fields..., geom_types...)
    types= (1,2,3)
    Tables.Schema(names, types)
end

function Base.iterate(gt::GeoTable, st = 1)
    features = gt.features

    d = Dict{String, Vector}()

    geom = ""
    prop_1 = ""
    prop_2 = ""

    namesz =  String.(Tables.schema(gt.features).names)
    for i in namesz
        if(i == "Geometry")
            d[i] = []
        else
            d[i] = Any[]
            end
    end

    st > length(features) && return nothing
    for (k, v) in pairs(d)
        if k == "Geometry"
            val = features[st].geometry
            geom = val
        elseif k == namesz[1]
            val = features[st].properties[k]
            prop_1 = val
        elseif k == namesz[2]
        # else
            val = features[st].properties[k]
            prop_2 = val
        end

    end

    return GeoTableRow(prop_1, prop_2, geom), st + 1
end

# struct GeoTable_sp{T<:Any}
#     sp::T
# end

# function GeoTable_sp(sp::T) where T<:Any
#     if(sp === nothing)
#         error("Not a valid geojson")
#     end
#     GeoTable_sp{T}(sp)
# end

# Tables.istable(::Type{<:GeoTable_sp}) = true
# Tables.rowaccess(::Type{<:GeoTable_sp}) = true
# Tables.rows(gt_sp::GeoTable_sp) = gt_sp

# Base.IteratorSize(::Type{<:GeoTable_sp}) = Base.HasLength()
# Base.length(gt_sp::GeoTable_sp) = (length(gt_sp.sp))
# Base.IteratorEltype(::Type{<:GeoTable_sp}) = Base.HasEltype()

# # read only AbstractVector
# Base.size(gt_sp::GeoTable_sp) = length(gt_sp.sp)  #use getter fxn here
# Base.getindex(gt_sp::GeoTable_sp, i) =  GeoTableRow_sp(sp(gt_sp)[i])
# Base.IndexStyle(::Type{<:GeoTable_sp}) = IndexLinear()

# # miss(x) = ifelse(x === nothing, missing, x)

# struct GeoTableRow_sp{T}
#     sp_f::T
# end

# Base.propertynames(gtr_sp::GeoTableRow_sp) = propertynames(properties(gtr_sp))

# # "Access the properties JSON3.Object of a Feature"
# # properties(gtr_sp::GeoTableRow_sp) = getfield(sp_f(gtr_sp), :record) ; println("prop got")
# # "Access the JSON3.Object that represents the Feature"
# # sp_f(gtr_sp::GeoTableRow_sp) = (gtr_sp.sp_f)[1] ; println("fetched")
# # "Access the JSON3.Array that represents the FeatureCollection"
# # # sp(gt_sp::GeoTable_sp) = getfield(gt_sp, :sp)
# "Access the JSON3.Object that represents the Feature's geometry"
# geometry(gtr_sp::GeoTableRow_sp) = Shapefile.shape(gtr_sp.sp_f)

# function Base.getproperty(gtr_sp::GeoTableRow_sp, sym::Symbol)
#     nm = sym
#     # println(nm);
#     props = properties(gtr_sp) #todo
#     println("got it")
#     val = props.nm
#     println("passed all")
# end

# @inline function Base.iterate(gt_sp::GeoTable_sp)
#     st = iterate(gt_sp.sp)
#     st === nothing && return nothing
#     val, state = st
#     return GeoTableRow_sp(val), state
# end

# @inline function Base.iterate(gt_sp::GeoTable_sp, st)
#     st = iterate(gt_sp.sp, st)
#     st === nothing && return nothing
#     val, state = st
#     return GeoTableRow_sp(val), state
# end

# struct FeatureCollection{T} <: AbstractVector{eltype(T)}
#     json::T
# end

# function read(source)
#     fc = JSON3.read(source)
#     features = get(fc, :features, nothing)
#     println(features)
#     if get(fc, :type, nothing) == "FeatureCollection" && features isa JSON3.Array
#         FeatureCollection{typeof(features)}(features)
#     else
#         throw(ArgumentError("input source is not a GeoJSON FeatureCollection"))
#     end
# end

# Tables.istable(::Type{<:FeatureCollection}) = true
# Tables.rowaccess(::Type{<:FeatureCollection}) = true
# Tables.rows(fc::FeatureCollection) = fc

# Base.IteratorSize(::Type{<:FeatureCollection}) = Base.HasLength()
# Base.length(fc::FeatureCollection) = length(json(fc))
# Base.IteratorEltype(::Type{<:FeatureCollection}) = Base.HasEltype()

# # read only AbstractVector
# Base.size(fc::FeatureCollection) = size(json(fc))
# Base.getindex(fc::FeatureCollection, i) = Feature(json(fc)[i])
# Base.IndexStyle(::Type{<:FeatureCollection}) = IndexLinear()

# miss(x) = ifelse(x === nothing, missing, x)

# struct Feature{T}
#     json::T
# end

# # these features always have type="Feature", so exclude that
# # the keys in properties are added here for direct access
# Base.propertynames(f::Feature) = propertynames(properties(f))

# "Access the properties JSON3.Object of a Feature"
# properties(f::Feature) = json(f).properties
# "Access the JSON3.Object that represents the Feature"
# json(f::Feature) = getfield(f, :json)
# "Access the JSON3.Array that represents the FeatureCollection"
# json(f::FeatureCollection) = getfield(f, :json)
# "Access the JSON3.Object that represents the Feature's geometry"
# geometry(f::Feature) = json(f).geometry

# """
# Get a specific property of the Feature

# Returns missing for null/nothing or not present, to work nicely with
# properties that are not defined for every feature. If it is a table,
# it should in some sense be defined.
# """
# function Base.getproperty(f::Feature, nm::Symbol)
#     props = properties(f)
#     val = get(props, nm, missing)
#     miss(val)
# end

# @inline function Base.iterate(fc::FeatureCollection)
#     st = iterate(json(fc))
#     st === nothing && return nothing
#     val, state = st
#     return Feature(val), state
# end

# @inline function Base.iterate(gt::GeoTable_sp, st)
#     st = iterate(json(fc), st)
#     st === nothing && return nothing
#     val, state = st
#     return Feature(val), state
# end


text:
function Helmholtz{T<:AbstractFloat}(isz::Array{Int,1},isx::Array{Int,1},ot::Array{T,1},vp::Array{T,2},wav::Array{T,1},fmin::T,fmax::T,nf::Int,nt::Int,dz::T,dx::T,dt::T,ext::Int=10,atten_max::T=1.,flag=2)

# This function solves the acoustic Helmholtz equation for a prescribed
# source wavelet and range of frequencies. A complex-numbered attenuating
# boundary is used to damp the boundary reflections. This code is intended
# for simultaneous sources.
#
# INPUTS:     isz       - z-index for each source location
#             isx       - x-index for each source location
#             ot        - Activation time for each source
#             vp        - Acoustic velocity model
#             wav       - Source wavelet (time domain)
#             fmin      - Minimum frequency to model
#             fmax      - Maximum frequency to model
#             nf        - Number of frequency bins; 'nf' must be greater than or equal to 'length(wav)' and 'nt'
#             nt        - Number of time steps to perform modeling for
#             dz        - Grid spacing in z-direction
#             dx        - Grid spacing in x-direction
#             dt        - Sampling interval (time step)
#             ext       - Thickness (number of grid points) of absorbing boundary region
#             atten_max - Maximum complex amplitude in the attenuating boundary layer
#             flag      - Output frequency-domain data if flag==1, Output time-domain data if flag==2
#
# OUTPUTS:    u         - Solution to Helmhotz equation; stored as a 3D volume with dimensions (w,z,x) or (t,z,x)

    nz = size(vp,1) + 2*ext
    nx = size(vp,2) + 2*ext

    L = Laplacian(nz,nx,dz,dx)
    M = MassMatrix(vp,nz,nx,ext) .* Attenuation(nz,nx,ext,atten_max)

    WAV = fft([wav;zeros(nf-length(wav))])
    fs = 1/dt
    df = fs/nf
    faxis = fftshift(-fs/2:df:fs/2-df)
    waxis = 2*pi*faxis

    _,iwmin = findmin(abs(fmin-faxis))
    _,iwmax = findmin(abs(fmax-faxis))
    U = complex(zeros(nf,nz,nx))
    for iw = iwmin:iwmax
        w = waxis[iw]
        H = L + w^2*M
        s = Source(isz,isx,ot,WAV,waxis,w,nz,nx,ext)
        U_tmp = H\s
        U_tmp = reshape(U_tmp,nz,nx)
        U[iw,:,:] = U_tmp
        U[nf-iw+2,:,:] = conj(U_tmp)
    end

    if flag==1
        return U
    elseif flag==2
        return real(ifft(U,1))[1:nt,ext+1:end-ext,ext+1:end-ext]
    end

end

function Helmholtz{T<:AbstractFloat}(isz::Int,isx::Int,ot::T,vp::Array{T,2},wav::Array{T,1},fmin::T,fmax::T,nf::Int,nt::Int,dz::T,dx::T,dt::T,ext::Int=10,atten_max::T=1.;flag=2)

# This function solves the acoustic Helmholtz equation for a prescribed
# source wavelet and range of frequencies. A complex-numbered attenuating
# boundary is used to damp the boundary reflections. This code is intended
# for a single source.
#
# INPUTS:     isz       - z-index for source location
#             isx       - x-index for source location
#             ot        - Activation time for the source
#             vp        - Acoustic velocity model
#             wav       - Source wavelet (time domain)
#             fmin      - Minimum frequency to model
#             fmax      - Maximum frequency to model
#             nf        - Number of frequency bins; 'nf' must be greater than or equal to 'length(wav)' and 'nt'
#             nt        - Number of time steps to perform modeling for
#             dz        - Grid spacing in z-direction
#             dx        - Grid spacing in x-direction
#             dt        - Sampling interval (time step)
#             ext       - Thickness (number of grid points) of absorbing boundary region
#             atten_max - Maximum complex amplitude in the attenuating boundary layer
#             flag      - Output frequency-domain data if flag==1, Output time-domain data if flag==2
#
# OUTPUTS:    U         - Solution to Helmhotz equation; stored as a 3D volume with dimensions (w,z,x) or (t,z,x)
#
# Contributors: Landon Safron
# Institution: University of Alberta
# Date: March 1, 2017

    nz = size(vp,1) + 2*ext
    nx = size(vp,2) + 2*ext

    L = Laplacian(nz,nx,dz,dx)
    M = MassMatrix(vp,nz,nx,ext) .* Attenuation(nz,nx,ext,atten_max)

    WAV = fft([wav;zeros(nf-length(wav))])
    fs = 1/dt
    df = fs/nf
    faxis = fftshift(-fs/2:df:fs/2-df)
    waxis = 2*pi*faxis

    _,iwmin = findmin(abs(fmin-faxis))
    _,iwmax = findmin(abs(fmax-faxis))
    U = complex(zeros(nf,nz,nx))
    for iw = iwmin:iwmax
        w = waxis[iw]
        H = L + w^2*M
        s = Source(isz,isx,ot,WAV,waxis,w,nz,nx,ext)
        U_tmp = H\s
        U_tmp = reshape(U_tmp,nz,nx)
        U[iw,:,:] = U_tmp
        U[nf-iw+2,:,:] = conj(U_tmp)
    end

    if flag==1
        return U[:,ext+1:end-ext,ext+1:end-ext]
    elseif flag==2
        return real(ifft(U,1))[1:nt,ext+1:end-ext,ext+1:end-ext]
    end

end


text:
using DataStructures

using PyCall
p_VIEW = p["VIEW"]
p_STRUCT = p["STRUCT"]
p_MKPOL = p["MKPOL"]

using LinearAlgebraicRepresentation
using Plasm

import Base.cat

L = LinearAlgebraicRepresentation

View = Plasm.view

"""
	apply(affineMatrix::Matrix)(larmodel::LAR)::LAR

Apply the `affineMatrix` parameter to the vertices of `larmodel`.

# Example

```
julia> square = LinearAlgebraicRepresentation.cuboid([1,1])
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]])

julia> Plasm.apply(LinearAlgebraicRepresentation.t(1,2))(square)
([1.0 1.0 2.0 2.0; 2.0 3.0 2.0 3.0], Array{Int64,1}[[1, 2, 3, 4]])
```
"""
function apply(affineMatrix)
	function apply0(larmodel)
		return L.struct2lar(L.Struct([ affineMatrix,larmodel ]))
	end
	return apply0
end

"""
	comp(funs::Array)

Standard mathematical composition.

Pipe execution from right to left on application to actual parameter.
"""
function comp(funs)
    function compose(f,g)
	  return x -> f(g(x))
	end
    id = x->x
    return reduce(compose, funs; init=id)
end

"""
	cons(funs::Array)(x::Any)::Array

*Construction* functional of FL and PLaSM languages.

Provides a *vector* functional that returns the array of
applications of component functions to actual parameter.

# Example

```
julia> Plasm.cons([cos,sin])(0)
2-element Array{Float64,1}:
 1.0
 0.0
```
"""
function cons(funs)
	return x -> [f(x) for f in funs]
end

"""
	k(Any)(x)

*Constant* functional of FL and PLaSM languages.

Gives a constant functional that always returns the actual parameter
when applied to another parameter.

#	Examples

```
julia> Plasm.k(10)(100)
10

julia> Plasm.k(sin)(cos)
sin
```
"""
function k(Any)
	x->Any
end

"""
	aa(fun::Function)(args::Array)::Array

AA applies fun to each element of the args sequence

# Example

```
julia> Plasm.aa(sqrt)([1,4,9,16])
4-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0
```
"""
function aa(fun)
	function aa1(args::Array)
		map(fun,args)
	end
	return aa1
end

"""
	id(x::Anytype)

Identity function.  Return the argument.

"""
id = x->x

"""
	distr(args::Union{Tuple,Array})::Array

Distribute right. The parameter `args` must contain a `list` and an element `x`.
Return the `pair` array with the elements of `args` coupled with `x`

# Example

```
julia> Plasm.distr(([1,2,3],10))
3-element Array{Array{Int64,1},1}:
 [1, 10]
 [2, 10]
 [3, 10]

julia> Plasm.distr([[1,2,3],10])
3-element Array{Array{Int64,1},1}:
 [1, 10]
 [2, 10]
 [3, 10]
```
"""
function distr(args)
	list,element = args
	return [ [e,element] for e in list ]
end

"""
	distl(args::Union{Tuple,Array})::Array

Distribute right. The parameter `args` must contain an element `x` and a `list`.
Return the `pair` array with `x` coupled with the elements of `args`.

# Example

```
julia> Plasm.distl((10, [1,2,3]))
3-element Array{Array{Int64,1},1}:
 [10, 1]
 [10, 2]
 [10, 3]

julia> Plasm.distl([10, [1,2,3]])
3-element Array{Array{Int64,1},1}:
 [10, 1]
 [10, 2]
 [10, 3]
```
"""
function distl(args)
	element, list = args
	return [ [element, e] for e in list ]
end

# Vector definition of printable ASCII codes as one-dimensional LAR models.
# Font design for *Geometric Programming for Computer-Aided Design*, Wiley, 2003.

ascii32 = ([0.0; 0.0],Array{Int64,1}[[1]])
ascii33 = ([1.75 1.75 2.0 2.0 1.5 1.5; 1.75 5.5 0.5 1.0 1.0 0.5],Array{Int64,1}[[1,2],[3,4],[4,5],[5,6],[6,3]])
ascii34 = ([1.0 2.0 2.0 2.0 1.5 1.5 2.0 3.0 3.0 3.0 2.5 2.5; 4.0 5.0 5.5 6.0 6.0 5.5 4.0 5.0 5.5 6.0 6.0 5.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,3],[7,8],[8,9],[9,10],[10,11],[11,12],[12,9]])
ascii35 = ([1.0 3.0 1.0 3.0 1.25 1.75 2.25 2.75; 2.5 2.5 3.5 3.5 1.75 4.0 1.75 4.0],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8]])
ascii36 = ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 2.0 2.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0 -0.5 6.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[13,14]])
ascii37 = ([2.5 2.5 2.0 2.0 2.5 2.5 2.0 2.0 1.5 3.5; 0.0 0.5 0.5 0.0 5.5 6.0 6.0 5.5 5.5 11.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,5],[9,10]])
ascii38 = ([4.0 3.0 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 5.0 4.0 2.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13],[13,14]])
ascii39 = ([1.0 2.0 2.0 2.0 1.5 1.5; 4.0 5.0 5.5 6.0 6.0 5.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,3]])
ascii40 = ([2.0 1.0 0.5 1.0 2.0; 0.0 1.0 3.0 5.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
ascii41 = ([2.0 3.0 3.5 3.0 2.0; 0.0 1.0 3.0 5.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
ascii42 = ([1.0 3.0 2.0 2.0 1.0 3.0 1.0 3.0; 3.0 3.0 2.0 4.0 2.0 4.0 4.0 2.0],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8]])
ascii43 = ([1.0 3.0 2.0 2.0; 3.0 3.0 2.0 4.0],Array{Int64,1}[[1,2],[3,4]])
ascii44 = ([1.0 2.0 2.0 2.0 1.5 1.5; -1.0 0.0 0.5 1.0 1.0 0.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,3]])
ascii45 = ([1.0 3.0; 3.0 3.0],Array{Int64,1}[[1,2]])
ascii46 = ([2.0 2.0 1.5 1.5 2.0; 0.0 0.5 0.5 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
ascii47 = ([1.0 3.0; 0.0 6.0],Array{Int64,1}[[1,2]])
ascii48 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[4,8]])
ascii49 = ([0.0 2.0 2.0 0.0 4.0; 4.0 6.0 0.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii50 = ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 4.0 5.0 6.0 6.0 5.0 4.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
ascii51 = ([0.0 4.0 2.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 2.0 1.0 0.0 0.0 1.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
ascii52 = ([4.0 0.0 4.0 4.0; 1.0 1.0 6.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
ascii53 = ([4.0 0.0 0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10]])
ascii54 = ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 6.0 6.0 5.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11]])
ascii55 = ([0.0 0.0 4.0 0.0; 5.0 6.0 6.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
ascii56 = ([1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 4.0 4.0 3.0 1.0 0.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 4.0 5.0 6.0 6.0 5.0 4.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[6,5],[5,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,6]])
ascii57 = ([0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 5.0 6.0 6.0 5.0 3.0 2.0 2.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11]])
ascii58 = ([2.0 2.0 1.5 1.5 2.0 2.0 1.5 1.5; 1.0 1.5 1.5 1.0 3.0 3.5 3.5 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,5]])
ascii59 = ([2.0 2.0 1.5 1.5 1.0 2.0 2.0 2.0 1.5 1.5; 3.0 3.5 3.5 3.0 -0.5 0.5 1.0 1.5 1.5 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,9],[9,10],[10,7]])
ascii60 = ([3.0 0.0 3.0; 6.0 3.0 0.0],Array{Int64,1}[[1,2],[2,3]])
ascii61 = ([1.0 3.0 1.0 3.0; 2.5 2.5 3.5 3.5],Array{Int64,1}[[1,2],[3,4]])
ascii62 = ([1.0 4.0 1.0; 6.0 3.0 0.0],Array{Int64,1}[[1,2],[2,3]])
ascii63 = ([2.0 2.0 1.5 1.5 1.75 1.75 3.0 3.0 2.0 1.0 0.0 0.0; 0.0 0.5 0.5 0.0 1.0 2.75 4.0 5.0 6.0 6.0 5.0 4.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,1],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12]])
ascii64 = ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 2.0 1.0 2.0 3.0 2.0; 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 1.0 2.0 3.0 2.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13]])
ascii65 = ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 0.0 5.0 6.0 6.0 5.0 0.0 2.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[7,8]])
ascii66 = ([0.0 0.0 3.0 4.0 4.0 3.0 4.0 4.0 3.0 0.0 0.0 3.0; 0.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[11,12]])
ascii67 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
ascii68 = ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 1.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
ascii69 = ([4.0 0.0 0.0 4.0 0.0 3.0; 0.0 0.0 6.0 6.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[5,6]])
ascii70 = ([0.0 0.0 4.0 0.0 3.0; 0.0 6.0 6.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii71 = ([2.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10]])
ascii72 = ([0.0 0.0 4.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0 3.0 3.0],Array{Int64,1}[[1,2],[3,4],[5,6]])
ascii73 = ([2.0 2.0 1.0 3.0 1.0 3.0; 0.0 6.0 0.0 0.0 6.0 6.0],Array{Int64,1}[[1,2],[3,4],[5,6]])
ascii74 = ([0.0 1.0 2.0 3.0 3.0 2.0 4.0; 1.0 0.0 0.0 1.0 6.0 6.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[6,7]])
ascii75 = ([4.0 0.0 4.0 0.0 0.0; 6.0 3.0 0.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii76 = ([4.0 0.0 0.0; 0.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3]])
ascii77 = ([0.0 0.0 2.0 4.0 4.0; 0.0 6.0 4.0 6.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5]])
ascii78 = ([0.0 0.0 4.0 4.0 4.0; 0.0 6.0 2.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii79 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
ascii80 = ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
ascii81 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0 1.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[10,11]])
ascii82 = ([0.0 0.0 3.0 4.0 4.0 3.0 0.0 3.0 4.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0 2.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[8,9]])
ascii83 = ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12]])
ascii84 = ([2.0 2.0 0.0 4.0; 0.0 6.0 6.0 6.0],Array{Int64,1}[[1,2],[3,4]])
ascii85 = ([0.0 0.0 1.0 3.0 4.0 4.0; 6.0 1.0 0.0 0.0 1.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6]])
ascii86 = ([0.0 2.0 4.0; 6.0 0.0 6.0],Array{Int64,1}[[1,2],[2,3]])
ascii87 = ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 6.0 3.0 0.0 3.0 0.0 3.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
ascii88 = ([0.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0],Array{Int64,1}[[1,2],[3,4]])
ascii89 = ([0.0 2.0 4.0 2.0 2.0; 6.0 2.0 6.0 2.0 0.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii90 = ([0.0 4.0 0.0 4.0; 6.0 6.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
ascii91 = ([2.0 1.0 1.0 2.0; 0.0 0.0 6.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
ascii92 = ([1.0 3.0; 6.0 0.0],Array{Int64,1}[[1,2]])
ascii93 = ([2.0 3.0 3.0 2.0; 0.0 0.0 6.0 6.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
ascii94 = ([1.0 2.0 3.0; 5.0 6.0 5.0],Array{Int64,1}[[1,2],[2,3]])
ascii95 = ([1.0 4.0; 0.0 0.0],Array{Int64,1}[[1,2]])
ascii96 = ([2.0 2.0 3.0 2.5 2.5 2.0; 4.5 5.0 6.0 4.0 4.5 4.0],Array{Int64,1}[[1,2],[2,3],[4,5],[5,1],[1,6],[6,4]])
ascii97 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10]])
ascii98 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0 1.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10],[10,11]])
ascii99 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
ascii100 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0 3.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10],[10,11]])
ascii101 = ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
ascii102 = ([4.0 4.0 3.0 2.0 1.0 1.0 0.0 2.0; 3.0 4.0 5.0 5.0 4.0 0.0 1.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[7,8]])
ascii103 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 -1.0 -1.0 0.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[1,9],[9,10],[10,11],[11,12]])
ascii104 = ([4.0 4.0 3.0 1.0 0.0 0.0 0.0 1.0; 0.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[6,7],[7,8]])
ascii105 = ([1.0 3.0 1.0 3.0 2.0 2.0 2.25 2.25 1.75 1.75; 0.0 0.0 3.0 3.0 0.0 3.0 3.75 4.25 4.25 3.75],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8],[8,9],[9,10],[10,7]])
ascii106 = ([1.0 3.0 2.0 2.0 1.0 0.0 2.25 2.25 1.75 1.75; 3.0 3.0 3.0 0.0 -1.0 0.0 3.75 4.25 4.25 3.75],Array{Int64,1}[[1,2],[3,4],[4,5],[5,6],[7,8],[8,9],[9,10],[10,7]])
ascii107 = ([0.0 1.0 1.0 0.0 4.0 2.0 1.0 3.0 4.0; 0.0 0.0 3.0 3.0 0.0 0.0 1.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[5,6],[6,7],[7,8],[8,9]])
ascii108 = ([2.0 2.0 1.0 1.0 3.0; 0.0 5.0 5.0 0.0 0.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii109 = ([4.0 4.0 2.0 2.0 2.0 0.0 0.0 0.0; 0.0 3.0 2.0 0.0 3.0 2.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5],[5,6],[7,8]])
ascii110 = ([3.0 3.0 1.0 1.0 1.0; 0.0 3.0 2.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii111 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]])
ascii112 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10]])
ascii113 = ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,1],[9,10]])
ascii114 = ([0.0 2.0 1.0 1.0 1.0 2.0 3.0 4.0; 0.0 0.0 0.0 3.0 2.0 3.0 3.0 2.0],Array{Int64,1}[[1,2],[3,4],[5,6],[6,7],[7,8]])
ascii115 = ([0.0 4.0 3.0 1.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 1.0 2.0 3.0 3.0 2.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]])
ascii116 = ([1.0 3.0 2.0 2.0 2.0 3.0; 0.0 0.0 0.0 5.0 4.0 4.0],Array{Int64,1}[[1,2],[3,4],[5,6]])
ascii117 = ([0.0 1.0 1.0 2.0 3.0 4.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
ascii118 = ([0.0 1.0 3.0 4.0; 3.0 0.0 3.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4]])
ascii119 = ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 3.0 2.0 0.0 2.0 0.0 2.0 3.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
ascii120 = ([0.0 1.0 4.0 1.0 4.0; 3.0 3.0 0.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5]])
ascii121 = ([0.0 1.0 2.5 0.0 1.0 4.0; 3.0 3.0 1.5 0.0 0.0 3.0],Array{Int64,1}[[1,2],[2,3],[4,5],[5,6]])
ascii122 = ([0.0 0.0 3.0 0.0 3.0 4.0; 2.0 3.0 3.0 0.0 0.0 1.0],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6]])
ascii123 = ([2.5 2.0 2.0 1.5 2.0 2.0 2.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
ascii124 = ([2.0 2.0; 0.0 5.0],Array{Int64,1}[[1,2]])
ascii125 = ([1.5 2.0 2.0 2.5 2.0 2.0 1.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5],Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]])
ascii126 = ([1.0 1.75 2.75 3.5; 5.0 5.5 5.0 5.5],Array{Int64,1}[[1,2],[2,3],[3,4]])

"""
	hpcs::Array

"""
hpcs = [
	ascii32,ascii33,ascii34,ascii35,ascii36,ascii37,ascii38,ascii39,ascii40,ascii41,
	ascii42,ascii43,ascii44,ascii45,ascii46,ascii47,ascii48,ascii49,ascii50,ascii51,
	ascii52,ascii53,ascii54,ascii55,ascii56,ascii57,ascii58,ascii59,ascii60,ascii61,
	ascii62,ascii63,ascii64,ascii65,ascii66,ascii67,ascii68,ascii69,ascii70,ascii71,
	ascii72,ascii73,ascii74,ascii75,ascii76,ascii77,ascii78,ascii79,ascii80,ascii81,
	ascii82,ascii83,ascii84,ascii85,ascii86,ascii87,ascii88,ascii89,ascii90,ascii91,
	ascii92,ascii93,ascii94,ascii95,ascii96,ascii97,ascii98,ascii99,ascii100,ascii101,
	ascii102,ascii103,ascii104,ascii105,ascii106,ascii107,ascii108,ascii109,ascii110,ascii111,
	ascii112,ascii113,ascii114,ascii115,ascii116,ascii117,ascii118,ascii119,ascii120,ascii121,
	ascii122,ascii123,ascii124,ascii125,ascii126 ]

"""
ascii_LAR = DataStructures.OrderedDict{Int,LinearAlgebraicRepresentation.LAR}()

	ascii_LAR::{Int,LinearAlgebraicRepresentation.LAR}

*Ordered dictionary* of printable ASCII codes as one-dimensional *LAR models* in 2D.

`Key` is the `integer` ASCII code between 32 and 126.
Font design: *Geometric Programming for Computer-Aided Design*, Wiley, 2003.

# Example
```
julia> Plasm.ascii_LAR[46]
([2.0 2.0 … 1.5 2.0; 0.0 0.5 … 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])

julia> Plasm.ascii_LAR[126]
([1.0 1.75 2.75 3.5; 5.0 5.5 5.0 5.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
```
"""
ascii_LAR = OrderedDict(zip(32:126,Plasm.hpcs))

# Attributes for `text` 2D graphics primitive.
# Reference to GKS ISO Graphics Standard (ISO/IEC 7942-4:1998)

const textalignment = "centre" #default value
const textangle = pi/4 #default value
const textwidth = 1.0 #default value
const textheight = 1.0 #default value
const textspacing = 0.25 #default value
const fontwidth = 4.0 #default value
const fontheight = 8.0 #default value
const fontspacing = 1.0 #default value

"""
	charpols(mystring::String)::Array{LAR,1}

Return the `array` of `LAR` models, for chars in `mystring`.
"""
function charpols(mystring)
    return [ascii_LAR[code] for code in [Int(char) for char in mystring]]
end

"""
	charseq(mystring::String)::Array{Char,1}

# Example
```
julia> Plasm.charseq("PLaSM")
5-element Array{Char,1}:
 'P'
 'L'
 'a'
 'S'
 'M'
```
"""
function charseq(mystring)
	return [char for char in mystring]
end

"""
	text(mystring::String)::LAR

Compute the one-dim *LAR model* drawing the contents of `mystring`

# Example
```
julia> model = Plasm.text("PLaSM")
# output
([0.0 0.0 3.0 4.0 4.0 3.0 0.0 9.0 5.0 5.0 14.0 13.0 11.0 10.0
10.0 11.0 13.0 14.0 14.0 14.0 15.0 16.0 18.0 19.0 19.0 18.0 16.0 15.0 15.0 16.0 18.0
19.0 20.0 20.0 22.0 24.0 24.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0 0.0 0.0 6.0 1.0 0.0 0.0 1.0
2.0 3.0 3.0 2.0 0.0 3.0 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0 0.0 6.0 4.0
6.0 0.0],
Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[8,9],[9,10],[11,12],[12,13],[13,14],
[14,15],[15,16],[16,17],[17,18],[18,11],[19,20],[21,22],[22,23],[23,24],[24,25],[25,26],
[26,27],[27,28],[28,29],[29,30],[30,31],[31,32],[33,34],[34,35],[35,36],[36,37]])

julia> Plasm.view(model)
```
"""
function text(mystring)
	out = comp([ L.struct2lar, L.Struct, cat, distr,
			cons([ charpols, k(L.t(fontspacing+fontwidth,0)) ]),charseq ])(mystring)
	return out
end

"""
	a2a(mat::Matrix)(models::Array{LAR,1})::Struct.body

Local function, service to `textWithAttributes` implementation.
"""
function a2a(mat)
	function a2a0(models)
		assembly = []
		for model in models
			push!( assembly, LinearAlgebraicRepresentation.Struct([ mat,model ]) )
		end
		assembly
	end
	return a2a0
end

"""
	translate(c::Number)(lar::LAR)::LAR

Local function, service to `textWithAttributes` implementation.
"""
function translate(c)
	function translate0(lar)
		xs = lar[1][1,:]
		width = maximum(xs) - minimum(xs)
		apply(LinearAlgebraicRepresentation.t(width/c,0))(lar)
	end
	return translate0
end

"""
	align(textalignment="centre"::String)(model::LAR)::LAR

Local function, service to `textWithAttributes` implementation.
"""
function align(textalignment="centre")
	function align1(model)
		if ( textalignment == "centre" ) out=translate(-2)(model)
		elseif ( textalignment == "right" ) out=translate(-1)(model)
		elseif ( textalignment == "left" ) out=model
		end
		return out
	end
end

"""
 	cat(args)

Redefined locally, as service to `textWithAttributes` implementation.
"""
function cat(args)
	return reduce( (x,y) -> append!(x,y), args; init=[] )
end

"""
	textWithAttributes(textalignment='centre', textangle=0,
		textwidth=1.0, textheight=2.0, textspacing=0.25)(strand::String)::LAR

Partial implementation of the GKS's graphics primitive `text`.

# Example

```
Plasm.view(Plasm.textWithAttributes("left", pi/4)("PLaSM"))
```
"""
function textWithAttributes(textalignment="centre", textangle=0,
							textwidth=1.0, textheight=2.0, textspacing=0.25)
	function textWithAttributes(strand)
		id = x->x
		mat = LinearAlgebraicRepresentation.s(textwidth/fontwidth,textheight/fontheight)
		comp([
		   apply(LinearAlgebraicRepresentation.r(textangle)),
		   align(textalignment),
		   L.struct2lar,
		   L.Struct,
		   cat,
		   distr,
		   cons([ a2a(mat) ∘ charpols,
				k(LinearAlgebraicRepresentation.t(textwidth+textspacing,0)) ]),
		   charseq ])(strand)
	end
end

"""
	embed(n::Int)(model::LAR)::LAR

Embed a `LAR` `model` of dimension ``d`` in a space ``R^{d+n}``.
The embedding is done by adding ``d`` zero coordinates to each vertex.

# Example

```
julia> square = LinearAlgebraicRepresentation.cuboid([1,1])
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]])

julia> embed(1)(square)
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0; 0.0 0.0 0.0 0.0], Array{Int64,1}[[1, 2, 3, 4]])
```
"""
function embed(n)
	function embed0( larmodel )
		return [larmodel[1]; zeros(n,size(larmodel[1],2))], larmodel[2]
	end
	return embed0
end

"""
	numbering(model::LARmodel)

Display a *wireframe* of `model` cells, of any dimension, with their *ordinal number*.
Different `colors` and size are used for the various dimensional cells.

# Examples

```
model = LinearAlgebraicRepresentation.cuboidGrid([3,4,2], true)
Plasm.view(Plasm.numbering()(model))

model = LinearAlgebraicRepresentation.cuboidGrid([10,10], true)
Plasm.view(Plasm.numbering(1.5)(model))
```
"""
function numbering(numberSizeScaling=1)
	function numbering0(model)
		V,cells = model
		if size(V,1)==2
			V = Plasm.embed(1)(model)[1]
		end
		wireframe = Plasm.lar2hpc(V,cells[2])
		ns = numberSizeScaling
		gcode = Plasm.textWithAttributes("centre", 0, 0.1ns, 0.2ns, 0.025ns)
		scene = [wireframe]
		for (h,skel) in enumerate(cells)
			colors = [p.GREEN,p.YELLOW,p.CYAN,p.ORANGE]
			nums = []
			for (k,cell) in enumerate(skel)
				center = sum([V[:,v] for v in cell])/length(cell)
				code = Plasm.embed(1)( gcode(string(k)) )
				scaling = (0.6+0.1h,0.6+0.1h,1)
				push!(nums, LinearAlgebraicRepresentation.struct2lar( LinearAlgebraicRepresentation.Struct([
					LinearAlgebraicRepresentation.t(center...), LinearAlgebraicRepresentation.s(scaling...), code ]) ))
			end
			hpc = Plasm.lar2hpc(nums)
			push!( scene, p.COLOR(colors[h])(hpc) )
		end
		p.STRUCT( scene )
	end
	return numbering0
end


text:
"""
``(N-1)``-dimensional manifold in ``N``-dimensional space.
"""
struct Hypersurface{N,SS<:SearchSpace}
    manifold::Function
    parameter_space::SS

    function Hypersurface(manifold::Function, parameter_space::SS) where {SS<:SearchSpace}
        N = numdims(parameter_space)+1
        int_pt = manifold(0.5*(mins(parameter_space)+maxs(parameter_space)), Val{1})
        length(int_pt) == N || throw(DimensionMismatch())
        new{N,SS}(manifold, parameter_space)
    end
end

"""
    generate(surf::Hypersurface, fs::EpsBoxDominanceFitnessScheme, param_step = 0.1*fs.ϵ)

Generate the points of the hypersurface using the
discretization defined by ϵ-box fitness schema.
"""
@generated function generate(surf::Hypersurface{N},
                             fs::EpsBoxDominanceFitnessScheme{N,F},
                             ::Type{Val{NP}},
                             param_step::Vector{F} = 0.1*fs.ϵ) where {N,F,NP}
    quote
        #archive = EpsBoxArchive(fs)
        pf = Dict{NTuple{N,Int}, IndexedTupleFitness{N,F}}()
        param = fill!(Vector{Float64}(undef, N-1), 0.0)
        #hat_compare = HatCompare(fs)
        Base.Cartesian.@nloops $(N-1) t d->range(surf.parameter_space.mins[d], stop=surf.parameter_space.maxs[d],
                                                 length=ceil(Int, surf.parameter_space.deltas[d]/param_step[d])) d->param[d]=t_d begin
            fit = surf.manifold(param, Val{NP})
            if !isnafitness(fit, fs) # NA if given parameters do not correspond to any point on the manifold
                ifit = convert(IndexedTupleFitness, fit, fs)
                pf[ifit.index] = ifit
                #add_candidate!(archive, convert(IndexedTupleFitness, fit, fs), Individual())
            end
        end
        return pf
    end
end

"""
    nondominated(fitnesses, fit_scheme)

Filter `fitnesses` removing all dominated values.
"""
function nondominated(fitnesses, fit_scheme::EpsBoxDominanceFitnessScheme{N,F}) where {N,F}
    arch = EpsBoxArchive(fit_scheme)
    res = sizehint!(Vector{IndexedTupleFitness{N,F}}(), length(fitnesses))
    empty_params = Individual()
    for f in fitnesses
        add_candidate!(arch, f, empty_params)
    end
    IndexedTupleFitness{N,F}[fitness(pt) for pt in pareto_frontier(arch)]
end

"""
    distance(a::NTuple{N,F}, b::NTuple{N,F})

Euclidean distance from `a` to `b`.
"""
@generated function distance(a::NTuple{N,F}, b::NTuple{N,F}) where {N,F}
    quote
        res = 0.0
        Base.Cartesian.@nexprs $N i -> res += (a[i]-b[i])^2
        return sqrt(res)
    end
end

"""
    IGD(A::Vector{NTuple{N,F}}, B::Vector{NTuple{N,F}}, [two_sided=true])

The average Euclidean distance from the points of `A` to the points of `B`.
"""
IGD(A::Vector{NTuple{N,F}}, B::Vector{NTuple{N,F}}) where {N,F<:Number} =
    sum(a -> minimum(b -> distance(a, b), B), A) / length(A)

"""
    IGD(ref::Hypersurface, sol::Vector{FitIndividual}, [two_sided=true])

Average Euclidean distance from the exact Pareto frontier of the problem (`ref`)
to the solution (`sol`) produced by the optimization method.
If `two_sided` is on, returns the maximum of `IGD(sol, ref)` and `IGD(nondominated(ref), sol)`.
"""
function IGD(ref::Hypersurface{N}, sol::AbstractVector{T},
             fit_scheme::EpsBoxDominanceFitnessScheme{N,F},
             ::Type{Val{NP}},
             param_step::Vector{F} = 0.1*fit_scheme.ϵ, two_sided::Bool=true
) where {N,F<:Number,T<:FrontierIndividualWrapper,NP}
    @assert fitness_type(sol[1]) == NTuple{N,F}
    ref_indexed_pts = values(generate(ref, fit_scheme, Val{NP}, param_step))
    nondom_ref_pts = NTuple{N,F}[val.orig for val in nondominated(ref_indexed_pts, fit_scheme)]
    orig_sol_pts = NTuple{N,F}[fitness(candi) for candi in sol]
    # how much `sol` should be inflated to contain all nondom_ref_pts
    res = IGD(nondom_ref_pts, orig_sol_pts)
    if two_sided
        # how much `ref` should be inflated to contain all `sol` points
        return max(res, IGD(orig_sol_pts, NTuple{N,F}[val.orig for val in ref_indexed_pts]))
    else
        return res
    end
end

"""
CEC09 Unconstrained Problem 8 objective function.
"""
function CEC09_UP8(x::Vector{Float64})
    N = length(x)
    N >= 3 || throw(ArgumentError("CEC09_UP8(): insufficient dimensions"))
    sumix = 3
    nterms = zeros(Int, 3)
    sums = zeros(Float64, 3)
    for j in 3:N
        sums[sumix] += (x[j] - 2*x[2]*sin(2π*x[1] + π/N*j))^2
        nterms[sumix] += 1
        sumix += 1
        if sumix > 3
            sumix = 1
        end
    end
    for i in 1:3
        if nterms[i] > 0
            sums[i] *= 2/nterms[i]
        end
    end
    (cos(0.5*π*x[1]) * cos(0.5*π*x[2]) + sums[1],
     cos(0.5*π*x[1]) * sin(0.5*π*x[2]) + sums[2],
     sin(0.5*π*x[1]) + sums[3])
end

"""
CEC09 Unconstrained Problem 8 Pareto Frontier.
Parameterized by t[1]=f[1] and t[2]=f[2].
"""
function CEC09_UP8_PF(t::Vector{Float64}, ::Type{Val{NP}}) where NP
    d=sum(abs2, t)
    (t[1], t[2], d <= 1.0 ? (1.0 - d)^(1/3) : NaN)
end

"""
The collection of CEC09 unconstrained
multi-objective problems.

See http://dces.essex.ac.uk/staff/zhang/MOEAcompetition/cec09testproblem0904.pdf.pdf
"""
const CEC09_Unconstrained_Set = Dict{Int,FunctionBasedProblemFamily}(
    8 => FunctionBasedProblemFamily(CEC09_UP8, "CEC09 UP8",  ParetoFitnessScheme{3}(is_minimizing=true),
                                    (-2.0, 2.0),
                                    Hypersurface(CEC09_UP8_PF, symmetric_search_space(2, (0.0, 1.0))),
                                    symmetric_search_space(2, (0.0, 1.0)))
)

schaffer1(x) = (sum(abs2, x), sum(xx -> abs2(xx - 2.0), x))
schaffer1_PF(t, ::Type{Val{NP}}) where {NP} = (NP*t[1]^2, NP*(2-t[1])^2)

const Schaffer1Family = FunctionBasedProblemFamily(schaffer1, "Schaffer1", ParetoFitnessScheme{2}(is_minimizing=true),
                                (-10.0, 10.0),
                                Hypersurface(schaffer1_PF, symmetric_search_space(1, (0.0, 2.0))))


text:
using ACEbonds
using Test

@testset "ACEbonds.jl" begin
    # Write your tests here.
    @testset "Bonds basics" begin include("test_bonds.jl"); end
    @testset "Invariant Cylindrical" begin include("test_invcyl.jl"); end
    @testset "Calculator" begin include("test_bondpot.jl"); end
end


text:
function get_frequencyweight_tf(responsetype::FilterType)
    designmethod = Butterworth(2)
    digitalfilter(responsetype, designmethod)
end

get_frequencyweight_tf(G::FilterCoefficients) = G

function frequency_weight(system, N)
    ω = range(0, stop = pi, length = N)
    G = get_frequencyweight_tf(system)
    Φ = freqz(G, ω)
    acf = irfft(abs2.(Φ), 2length(Φ) - 1)[1:end÷2+1]
    n = length(acf)

    lastlarge = findlast(x -> abs(x) > sqrt(eps()) * acf[1], acf)
    if lastlarge === nothing || lastlarge == n
        lastlarge = n - 1
    end
    bands = [0 => Fill(acf[1], n)]
    for i = 1:lastlarge
        slice = Fill(acf[i+1], n - i)
        push!(bands, i => slice)
        push!(bands, -i => slice)
    end
    W = BandedMatrix((bands...,), (n, n))

    # W = similar(acf, n, n)
    # for di = -n+1:n-1
    #     W[diagind(W,di)] .= acf[abs(di)+1]
    # end
    # Symmetric(W)
end

weighted_estimator(H::FilterType) = (A, y) -> wls(A, y, frequency_weight(H, size(y, 1)))

"""
    frequency_focus(d::AbstractIdData, responsetype::FilterType)

Filter both input and output of the identification data using zero-phase filtering (`filtfilt`).
Since both input and output is filtered, linear identification will not be affected in any other way than to focus the fit on the selected frequency range, i.e. the range that has high gain in the provided filter. Note, if the system that generated `d` is nonlinear, identification might be severely impacted by this transformation. Verify linearity with, e.g., `coherenceplot`.
"""
function frequency_focus(d::AbstractIdData, responsetype::FilterType)
    H = get_frequencyweight_tf(responsetype)
    y = time1(output(d))
    u = time1(input(d))
    y = filtfilt(H, y)
    u = filtfilt(H, u)
    iddata(y', u', d.Ts)
end

"""
    frequency_focus(d::AbstractIdData, l::Number, u::Number)

Filter input and output with a bandpass filter between `l` and `u` Hz. If `l = 0` a lowpass filter will be used, and if `u = Inf` a highpass filter will be used.
"""
function frequency_focus(d::AbstractIdData, l::Number, u::Number)
    responsetype = if u == Inf
        Highpass(l, fs = d.fs)
    elseif l <= 0
        Lowpass(u, fs = d.fs)
    else
        Bandpass(l, u, fs = d.fs)
    end
    frequency_focus(d, responsetype)
end


text:
using DataFrames
using CSV

products = CSV.read("instacart-market-basket-analysis/products.csv");
departments = CSV.read("instacart-market-basket-analysis/departments.csv");
aisles = CSV.read("instacart-market-basket-analysis/aisles.csv");

function get_department_id(category_string)
    return departments[departments.department .== category_string, :department_id][1]
end
function get_products_in_department(department)
    products[products.department_id .== get_department_id(department), :]
end

excluded_departments = ["frozen",
    "pets", "personal care", "meat seafood", "household", "babies"];

included_products = products[
    indexin(products.department_id,
        [get_department_id(x) for x in excluded_departments]) .== nothing,
    :];

function get_random_product(products)
    maxvalue = size(products)[1]
    return products[rand(1:maxvalue), :]
end

function get_random_products(products, N)
    return DataFrame([get_random_product(products) for i in 1:N])
end

essential_departments = ["produce",
    "dry goods pasta", "dairy eggs", "canned goods"]

essential_products = products[
    indexin(products.department_id,
        [get_department_id(x) for x in essential_departments]) .!= nothing,
    :];

only_produce = ["produce"]

produce_products = products[
    indexin(products.department_id,
        [get_department_id(x) for x in only_produce]) .!= nothing,
    :];

println("Congratulations on yet another day of doing the #randomcookingchallenge!")
println("")

println("""RULES:
    1. You get two options. Choose one.
    2. Each option has 7 ingredients. You need to include at least 5 of them.
    3. You may add other ingredients as you see fit, but the resulting dish
       needs to be based around the random ingredients.""")
println("")

println("Good luck!")
println("")

@show get_random_products(produce_products, 7)

println("")

@show get_random_products(produce_products, 7)


text:
# Generic rand and rand! implementations similar to those in Distributions,
# but with an rng argument:

abstract type BATSampler{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S} end

@inline Random.rand(s::BATSampler, args...) = rand(Random.GLOBAL_RNG, s, args...)
@inline Random.rand!(s::BATSampler, args...) = rand!(Random.GLOBAL_RNG, s, args...)

# To avoid ambiguity with Distributions:
@inline Random.rand(s::BATSampler, dims::Int...) = rand(Random.GLOBAL_RNG, s, dims...)
@inline Random.rand!(s::BATSampler, A::AbstractVector) = rand!(Random.GLOBAL_RNG, s, A)
@inline Random.rand!(s::BATSampler, A::AbstractMatrix) = rand!(Random.GLOBAL_RNG, s, A)

function Random.rand!(rng::AbstractRNG, s::BATSampler{Univariate}, A::AbstractArray)
    @inbounds @simd for i in 1:length(A)
        A[i] = rand(rng, s)
    end
    return A
end

Random.rand(rng::AbstractRNG, s::BATSampler{Univariate}, dims::Dims) =
    rand!(rng, s, Array{eltype(s)}(undef, dims))

Random.rand(rng::AbstractRNG, s::BATSampler{Univariate}, dims::Int...) =
    rand!(s, Array{eltype(s)}(undef, dims))

Random.rand(rng::AbstractRNG, s::BATSampler{Multivariate}) =
    rand!(rng, s, Vector{eltype(s)}(undef, length(s)))

Random.rand(rng::AbstractRNG, s::BATSampler{Multivariate}, n::Integer) =
    rand!(rng, s, Matrix{eltype(s)}(undef, length(s), n))

function Random.rand!(rng::AbstractRNG, s::BATSampler{Multivariate}, A::AbstractMatrix)
    _check_rand_compat(s, A)
    @inbounds for i = 1:size(A,2)
        rand!(rng, s, view(A,:,i))
    end
    return A
end

Random.rand(rng::AbstractRNG, s::BATSampler{Multivariate}, n::Int) = rand!(rng, s, Matrix{eltype(s)}(undef, length(s), n))


text:
using uJSON
using JSON
using Base.Test

println("Running trivial example: ")
input_str = "[1,2,\"hello\", null, 3.2123, {\"a\":1}]"
uresult = uJSON.parse(input_str)
println(input_str, " >> ", uresult)

fn = Pkg.dir("uJSON", "test", "sample.json")
print("\nJSON running simple test:  ")
@time result = JSON.parse(open(fn, "r"))
print("uJSON running simple test: ")
@time uresult = uJSON.parse(open(fn, "r"))
@test uresult == result
println("SUCCESS: results equal")


text:
module TestInputOutput

using ModiaLang
using DifferentialEquations
@usingModiaPlot

FirstOrder = Model(
    T = 0.2,
    u = input | Map(start=0),
    y = output,
    x = Var(init=0.3),
    equations = :[
        T * der(x) + x = u,
        y = 2*x]
)

firstOrder = @instantiateModel(FirstOrder)

simulate!(firstOrder, Tsit5(), stopTime = 2, requiredFinalStates = [1.362130067500907e-5])

plot(firstOrder, ["u", "x", "der(x)", "y"])

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang