text:
"""
Estimate a linear model with high dimensional categorical variables / instrumental variables

### Arguments
* `f` : Formula,
* `df` : AbstractDataFrame
* `vcov_method` : An object of type AbstractVcovMethod. Default to VcovSimple(). For now, `VcovSimple()` (default), `VcovWhite()` and `VcovCluster(cols)` are implemented.
* `weight` : Symbol for weight variables. Corresponds to analytical weights
* `subset` : AbstractVector{Bool} for subsample
* `save` : SHould residuals and eventual estimated fixed effects saved in a dataframe?
* `maxiter` : Maximum number of iterations
* `tol` : tolerance
* `method` : A symbol for the method. Default is :lsmr (akin to conjugate gradient descent). Other choices are :qr and :cholesky (factorization methods)

### Returns
* `::AbstractRegressionResult` : a regression results

### Details
A typical formula is composed of one dependent variable, exogeneous variables, endogeneous variables, instruments, and high dimensional fixed effects
```
depvar ~ exogeneousvars + (endogeneousvars = instrumentvars) |> absorbvars
```
Categorical variable should be of type PooledDataArray.  See the following to create PooledDataArray:
* `pool` : transform one variable into a `PooledDataArray`.
* `group` : combine multiple variables into a `PooledDataArray`.
Models with instruments variables are estimated using 2SLS. `reg` tests for weak instruments by computing the Kleibergen-Paap rk Wald F statistic, a generalization of the Cragg-Donald Wald F statistic for non i.i.d. errors. The statistic is similar to the one returned by the Stata command `ivreg2`.

### Examples
```julia
using DataFrames, RDatasets, FixedEffectModels
df = dataset("plm", "Cigar")
df[:StatePooled] =  pool(df[:State])
df[:YearPooled] =  pool(df[:Year])
reg(Sales ~ Price |> StatePooled + YearPooled, df)
reg(Sales ~ NDI |> StatePooled + StatePooled&Year, df)
reg(Sales ~ NDI |> StatePooled*Year, df)
reg(Sales ~ (Price = Pimin), df)
reg(Sales ~ Price, df, weight = :Pop)
reg(Sales ~ NDI, df, subset = df[:State] .< 30)
reg(Sales ~ NDI, df, VcovWhite())
reg(Sales ~ NDI, df, VcovCluster([:StatePooled]))
reg(Sales ~ NDI, df, VcovCluster([:StatePooled, :YearPooled]))
```
"""

# TODO: minimize memory
function reg(f::Formula, df::AbstractDataFrame,
             vcov_method::AbstractVcovMethod = VcovSimple();
             weight::Union{Symbol, Void} = nothing,
             subset::Union{AbstractVector{Bool}, Void} = nothing,
             maxiter::Integer = 10000, tol::Real= 1e-8, df_add::Integer = 0,
             save::Bool = false,
             method::Symbol = :lsmr)

    ##############################################################################
    ##
    ## Parse formula
    ##
    ##############################################################################

    rf = deepcopy(f)
    (has_absorb, absorb_formula, absorb_terms,
        has_iv,iv_formula,iv_terms,endo_formula,endo_terms) = decompose!(rf)
    rt = Terms(rf)
    has_weight = weight != nothing

    # check depth 1 symbols in original formula are all PooledDataArray
    if has_absorb
        if isa(f.rhs.args[3], Symbol)
            x = f.rhs.args[3]
            !isa(df[x], PooledDataArray) && error("$x should be PooledDataArray")
        elseif f.rhs.args[3].args[1] == :+
            x = f.rhs.args[3].args
            for i in 2:length(x)
                isa(x[i], Symbol) && !isa(df[x[i]], PooledDataArray) && error("$(x[i]) should be PooledDataArray")
            end
        end
    end

    ##############################################################################
    ##
    ## Construct new dataframe after removing missing values
    ##
    ##############################################################################

    # create a dataframe without missing values & negative weights
    vars = allvars(rf)
    absorb_vars = allvars(absorb_formula)
    iv_vars = allvars(iv_formula)
    endo_vars = allvars(endo_formula)
    vcov_vars = allvars(vcov_method)

    # create a dataframe without missing values & negative weights
    all_vars = vcat(vars, vcov_vars, absorb_vars, endo_vars, iv_vars)
    all_vars = unique(convert(Vector{Symbol}, all_vars))
    esample = complete_cases(df[all_vars])

    if has_weight
        esample &= isnaorneg(df[weight])
    end
    if subset != nothing
        if length(subset) != size(df, 1)
            error("df has $(size(df, 1)) rows but the subset vector has $(length(subset)) elements")
        end
        esample &= convert(BitArray, subset)
    end
    nobs = sum(esample)
    (nobs > 0) || error("sample is empty")

    # Compute weight
    sqrtw = get_weight(df, esample, weight)

    # remove unusused levels
    subdf = df[esample, all_vars]
    main_vars = unique(convert(Vector{Symbol}, vcat(vars, endo_vars, iv_vars)))
    for v in main_vars
        # in case subdataframe, don't construct subdf[v] if you dont need to do it
        if typeof(df[v]) <: PooledDataArray
            dropUnusedLevels!(subdf[v])
        end
    end

    # Compute pfe, a FixedEffectProblem
    has_intercept = rt.intercept
    if has_absorb
        fixedeffects = FixedEffect(subdf, absorb_terms, sqrtw)
        # in case some FixedEffect does not have interaction, remove the intercept
        if any([typeof(f.interaction) <: Ones for f in fixedeffects])
            rt.intercept = false
            has_intercept = true
        end
        pfe = FixedEffectProblem(fixedeffects, Val{method})
    else
        pfe = nothing
    end

    # Compute data for std errors
    vcov_method_data = VcovMethodData(vcov_method, subdf)

    ##############################################################################
    ##
    ## Dataframe --> Matrix
    ##
    ##############################################################################

    # initialize iterations and converged
    iterations = Int[]
    converged = Bool[]

    mf = ModelFrame2(rt, subdf, esample)

    # Obtain y
    py = model_response(mf)[:]
    if eltype(py) != Float64
        y = convert(Vector{Float64}, py)
    else
        y = py
    end
    yname = rt.eterms[1]
    broadcast!(*, y, y, sqrtw)
    # old y will be used if fixed effects
    if has_absorb
        oldy = deepcopy(y)
    else
        oldy = y
    end
    residualize!(y, pfe, iterations, converged; maxiter = maxiter, tol = tol)

    # Obtain X
    coef_names = coefnames(mf)
    if isempty(mf.terms.terms) && mf.terms.intercept == false
        Xexo = Matrix{Float64}(sum(mf.msng), 0)
    else
        Xexo = ModelMatrix(mf).m
        if size(Xexo, 2) == 1
            # See pull request #1017 in DataFrames Package
            Xexo = deepcopy(Xexo)
        end
    end
    broadcast!(*, Xexo, Xexo, sqrtw)
    residualize!(Xexo, pfe, iterations, converged; maxiter = maxiter, tol = tol)

    # Obtain Xendo and Z
    if has_iv
        mf = ModelFrame2(endo_terms, subdf, esample)
        coef_names = vcat(coef_names, coefnames(mf))
        Xendo = ModelMatrix(mf).m
        broadcast!(*, Xendo, Xendo, sqrtw)
        residualize!(Xendo, pfe, iterations, converged; maxiter = maxiter, tol = tol)

        mf = ModelFrame2(iv_terms, subdf, esample)
        Z = ModelMatrix(mf).m
        broadcast!(*, Z, Z, sqrtw)
        residualize!(Z, pfe, iterations, converged; maxiter = maxiter, tol = tol)
    end

    # iter and convergence
    if has_absorb
        iterations = maximum(iterations)
        converged = all(converged)
    end

    ##############################################################################
    ##
    ## Regression
    ##
    ##############################################################################

    # Compute Xhat
    if has_iv
        if size(Z, 2) < size(Xendo, 2)
            error("Model not identified. There must be at least as many ivs as endogeneneous variables")
        end
        # get linearly independent columns
        baseall= basecol(Z, Xexo, Xendo)
        basecolXexo = baseall[(size(Z, 2)+1):(size(Z, 2) + size(Xexo, 2))]
        basecolXendo = baseall[(size(Z, 2) + size(Xexo, 2) + 1):end]
        Z = getcols(Z, baseall[1:size(Z, 2)])
        Xexo = getcols(Xexo, basecolXexo)
        Xendo = getcols(Xendo, basecolXendo)
        basecoef = vcat(basecolXexo, basecolXendo)

        # Build
        X = hcat(Xexo, Xendo)
        newZ = hcat(Xexo, Z)
        crossz = cholfact!(At_mul_B(newZ, newZ))
        Pi = crossz \ At_mul_B(newZ, Xendo)
        Xhat = hcat(Xexo, newZ * Pi)
        X = hcat(Xexo, Xendo)

        # prepare residuals used for first stage F statistic
        ## partial out Xendo in place wrt (Xexo, Z)
        Xendo_res = BLAS.gemm!('N', 'N', -1.0, newZ, Pi, 1.0, Xendo)
        ## partial out Z in place wrt Xexo
        Pi2 = cholfact!(At_mul_B(Xexo, Xexo)) \ At_mul_B(Xexo, Z)
        Z_res = BLAS.gemm!('N', 'N', -1.0, Xexo, Pi2, 1.0, Z)

        # free memory (not sure it helps)
        Xexo = nothing
    else
        # get linearly independent columns
        basecolXexo = basecol(Xexo)
        Xexo = getcols(Xexo, basecolXexo)
        Xhat = Xexo
        X = Xexo
        basecoef = basecolXexo
    end

    # Compute coef and residuals
    crossx =  cholfact!(At_mul_B(Xhat, Xhat))
    coef = crossx \ At_mul_B(Xhat, y)
    residuals = y - X * coef

    ##############################################################################
    ##
    ## Optionally save some vectors in a new dataframe
    ##
    ##############################################################################

    # save residuals in a new dataframe
    augmentdf = DataFrame()
    if save
        broadcast!(/, residuals, residuals, sqrtw)
        if all(esample)
            augmentdf[:residuals] = residuals
        else
            augmentdf[:residuals] =  DataArray(Float64, length(esample))
            augmentdf[esample, :residuals] = residuals
        end
        if has_absorb
            mf = ModelFrame2(rt, subdf, esample)
            oldX = ModelMatrix(mf).m
            if !all(basecoef)
                oldX = oldX[:, basecoef]
            end
            broadcast!(*, oldX, oldX, sqrtw)
            BLAS.gemm!('N', 'N', -1.0, oldX, coef, 1.0, oldy)
            axpy!(-1.0, residuals, oldy)

            augmentdf = hcat(augmentdf, getfe!(pfe, oldy, esample; tol = tol, maxiter = maxiter))
        end
    end

    ##############################################################################
    ##
    ## Test Statistics
    ##
    ##############################################################################

    # Compute degrees of freedom
    df_intercept = 0
    if has_absorb || rt.intercept
        df_intercept = 1
    end
    df_absorb = 0
    if has_absorb
        ## poor man adjustement of df for clustedered errors + fe: only if fe name != cluster name
        for fe in fixedeffects
            if typeof(vcov_method) == VcovCluster && in(fe.factorname, vcov_vars)
                df_absorb += 0
                else
                df_absorb += sum(fe.scale .!= zero(Float64))
            end
        end
    end
    nvars = size(X, 2)
    df_residual = max(1, nobs - nvars - df_absorb - df_add)

    # Compute ess, tss, r2, r2 adjusted
    ess = sumabs2(residuals)
    if has_absorb
        tss = compute_tss(y, rt.intercept, sqrtw)
        r2_within = 1 - ess / tss
    end
    tss = compute_tss(oldy, has_intercept, sqrtw)
    r2 = 1 - ess / tss
    r2_a = 1 - ess / tss * (nobs - has_intercept) / df_residual

    # Compute standard error
    vcov_data = VcovData(Xhat, crossx, residuals, df_residual)
    matrix_vcov = vcov!(vcov_method_data, vcov_data)

    # Compute Fstat
    (F, p) = compute_Fstat(coef, matrix_vcov, nobs, rt.intercept, vcov_method_data, vcov_data)

    # Compute Fstat of First Stage
    if has_iv
        Pip = Pi[(size(Pi, 1) - size(Z_res, 2) + 1):end, :]
        (F_kp, p_kp) = ranktest!(Xendo_res, Z_res, Pip,
                                  vcov_method_data, nvars, df_absorb)
    end

    ##############################################################################
    ##
    ## Return regression result
    ##
    ##############################################################################

    # add omitted variables
    if !all(basecoef)
        newcoef = fill(zero(Float64), length(basecoef))
        newmatrix_vcov = fill(NaN, (length(basecoef), length(basecoef)))
        newindex = [searchsortedfirst(cumsum(basecoef), i) for i in 1:length(coef)]
        for i in 1:length(coef)
            newcoef[newindex[i]] = coef[i]
            for j in 1:length(coef)
                newmatrix_vcov[newindex[i], newindex[j]] = matrix_vcov[i, j]
            end
        end
        coef = newcoef
        matrix_vcov = newmatrix_vcov
    end

    # return
    if !has_iv && !has_absorb
        return RegressionResult(coef, matrix_vcov, esample, augmentdf,
                                coef_names, yname, f, nobs, df_residual,
                                r2, r2_a, F, p)
    elseif has_iv && !has_absorb
        return RegressionResultIV(coef, matrix_vcov, esample, augmentdf,
                                  coef_names, yname, f, nobs, df_residual,
                                  r2, r2_a, F, p, F_kp, p_kp)
    elseif !has_iv && has_absorb
        return RegressionResultFE(coef, matrix_vcov, esample, augmentdf,
                                  coef_names, yname, f, nobs, df_residual,
                                  r2, r2_a, r2_within, F, p, iterations, converged)
    elseif has_iv && has_absorb
        return RegressionResultFEIV(coef, matrix_vcov, esample, augmentdf,
                                   coef_names, yname, f, nobs, df_residual,
                                   r2, r2_a, r2_within, F, p, F_kp, p_kp,
                                   iterations, converged)
    end
end

##############################################################################
##
## Fstat
##
##############################################################################

function compute_Fstat(coef::Vector{Float64}, matrix_vcov::Matrix{Float64},
    nobs::Int, hasintercept::Bool,
    vcov_method_data::AbstractVcovMethodData, vcov_data::VcovData)
    coefF = deepcopy(coef)
    # TODO: check I can't do better
    length(coef) == hasintercept && return NaN, NaN
    if hasintercept && length(coef) > 1
        coefF = coefF[2:end]
        matrix_vcov = matrix_vcov[2:end, 2:end]
    end
    F = (diagm(coefF)' * (matrix_vcov \ diagm(coefF)))[1]
    df_ans = df_FStat(vcov_method_data, vcov_data, hasintercept)
    dist = FDist(nobs - hasintercept, max(df_ans, 1))
    return F, ccdf(dist, F)
end

##############################################################################
##
## Weight
##
##############################################################################

function get_weight(df::AbstractDataFrame, esample::BitVector, weight::Symbol)
    out = df[esample, weight]
    # there are no NA in it. DataVector to Vector
    out = convert(Vector{Float64}, out)
    map!(sqrt, out, out)
    return out
end
get_weight(df::AbstractDataFrame, esample::BitVector, ::Void) = Ones{Float64}(sum(esample))

function compute_tss(y::Vector{Float64}, hasintercept::Bool, ::Ones)
    if hasintercept
        tss = zero(Float64)
        m = mean(y)::Float64
        @inbounds @simd  for i in 1:length(y)
            tss += abs2((y[i] - m))
        end
    else
        tss = sumabs2(y)
    end
    return tss
end

function compute_tss(y::Vector{Float64}, hasintercept::Bool, sqrtw::Vector{Float64})
    if hasintercept
        m = (mean(y) / sum(sqrtw) * length(y))::Float64
        tss = zero(Float64)
        @inbounds @simd  for i in 1:length(y)
            tss += abs2(y[i] - sqrtw[i] * m)
        end
    else
        tss = sumabs2(y)
    end
    return tss
end


text:
# Copyright 2022 RelationalAI, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import Base
using ArgParse: ArgName, ArgParseSettings, add_arg_table!, parse_args

struct Args
    _args::Dict{String,Any}
end

Base.getindex(args::Args, key::String) = getfield(args, :_args)[key]
Base.getproperty(args::Args, name::Symbol) = getfield(args, :_args)[String(name)]

function parseargs(table::Union{ArgName,Vector,Dict}...)
    s = ArgParseSettings()
    add_arg_table!(s, table...)
    args = parse_args(ARGS, s)
    return Args(args)
end

ext: jl

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
using FractionalDiffEq, Plots
α = 0.95; ϕ(x) = 0.5
τ = [2, 2.6]
fun(t, y, ϕ1, ϕ2) = 2*ϕ1/(1+ϕ2^9.65)-y
prob = FDDEProblem(fun, ϕ, α, τ, 100)
delayed, y = solve(prob, 0.01, DelayPECE())

p1=plot(delayed[1, :], y)
p2=plot(delayed[2, :], y)
plot(p1, p2, layout=(1, 2))


text:
import Test: @test, @testset
import Luna: Fields, FFTW, Grid, Maths, PhysData, Processing, Modes, Tools, Maths
import Statistics: mean, std
import Random: MersenneTwister

# note that most of the Fields.jl code is tested in many other modules

function getceo(t, Et, It, ω0)
    Δt = t[argmax(It)] - t[argmax(Et)]
    Δt*ω0
end

@testset "Wavelength" begin
    # real
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕ = [0.0, 0.0]
    grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)
    λ0 = 320e-9
    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)
    λ0 = 800e-9
    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)
    λ0 = 320e-9
    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)

    # Envelope
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕ = [0.0, 0.0]
    grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)
    λ0 = 320e-9
    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)
    λ0 = 800e-9
    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)
    λ0 = 320e-9
    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    @test isapprox(PhysData.wlfreq(grid.ω[argmax(abs2.(Eω))]), λ0, rtol=3e-4)
end

@testset "Energy" begin
    # real
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕ = [0.0, 0.0]
    grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    @test isapprox(energy_t(Et), energy, rtol=1e-14)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    @test isapprox(energy_t(Et), energy, rtol=1e-14)

    # Envelope
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕ = [0.0, 0.0]
    grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    @test isapprox(energy_t(Et), energy, rtol=1e-14)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    @test isapprox(energy_t(Et), energy, rtol=1e-14)
end

@testset "Duration" begin
    # real
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕ = [0.0, 0.0]
    grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=1e-5)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=1e-5)

    # Envelope
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕ = [0.0, 0.0]
    grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=2e-5)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(Maths.fwhm(grid.t, It), τfwhm, rtol=3e-5)
end

@testset "Position" begin
    # real
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    τ0 = 5e-15
    # elements of ϕ are [CEP, group delay, GDD, TOD, ...]
    # so [0.0, τ0] is a delay by τ0
    ϕ = [0.0, τ0]
    grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)

    # Envelope
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    τ0 = 5e-15
    # elements of ϕ are [CEP, group delay, GDD, TOD, ...]
    # so [0.0, τ0] is a delay by τ0
    ϕ = [0.0, τ0]
    grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)

    # non zero
    τ0 = -564e-15

    #real
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    τ0 = 5e-15
    # elements of ϕ are [CEP, group delay, GDD, TOD, ...]
    # so [0.0, τ0] is a delay by τ0
    ϕ = [0.0, τ0]
    grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)

    # Envelope
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    τ0 = 5e-15
    # elements of ϕ are [CEP, group delay, GDD, TOD, ...]
    # so [0.0, τ0] is a delay by τ0
    ϕ = [0.0, τ0]
    grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=ϕ)
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(grid.t[argmax(It)], τ0, rtol=1e-15, atol=1e-15)
end

@testset "CEO" begin
    # real
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕCEO = 0.0
    grid = Grid.RealGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(getceo(grid.t, Et, It, PhysData.wlfreq(λ0)), ϕCEO, rtol=1e-15, atol=1e-15)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Maths.hilbert(Et))
    @test isapprox(getceo(grid.t, Et, It, PhysData.wlfreq(λ0)), ϕCEO, rtol=1e-15, atol=1e-15)

    # Envelope
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    ϕCEO = 0.0
    grid = Grid.EnvGrid(1.0, λ0, (160e-9, 3000e-9), 10e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)

    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(getceo(grid.t, real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0)), ϕCEO, rtol=1e-15, atol=1e-15)

    input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
    Eω = input(grid, FT)
    Et = FT \ Eω
    It = abs2.(Et)
    @test isapprox(getceo(grid.t, real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0)), ϕCEO, rtol=1e-15, atol=1e-15)

    # non zero

    #real
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    τ0 = 0.0
    grid = Grid.RealGrid(1.0, λ0, (100e-9, 3000e-9), 1e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)

    # Make CEO exact multiple of one grid point to avoid issues with argmax() in getceo()
    δt = grid.t[2] - grid.t[1]
    for i = 1:10
        ϕCEO = i*δt*PhysData.wlfreq(λ0)

        input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
        Eω = input(grid, FT)
        Et = FT \ Eω
        It = abs2.(Maths.hilbert(Et))
        @test isapprox(abs(getceo(grid.t, Et, It, PhysData.wlfreq(λ0))), ϕCEO, rtol=1e-10)

        input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
        Eω = input(grid, FT)
        Et = FT \ Eω
        It = abs2.(Maths.hilbert(Et))
        @test isapprox(abs(getceo(grid.t, Et, It, PhysData.wlfreq(λ0))), ϕCEO, rtol=1e-10)
    end

    # Envelope
    τfwhm = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    τ0 = 0.0
    grid = Grid.EnvGrid(1.0, λ0, (100e-9, 3000e-9), 1e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)

    # Make CEO exact multiple of one grid point to avoid issues with argmax() in getceo()
    δt = grid.t[2] - grid.t[1]

    for i = 1:10
        ϕCEO = i*δt*PhysData.wlfreq(λ0)

        input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
        Eω = input(grid, FT)
        Et = FT \ Eω
        It = abs2.(Et)
        @test isapprox(
            abs(getceo(grid.t, real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0))),
            ϕCEO,
            rtol=1e-10)

        input = Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy, ϕ=[ϕCEO])
        Eω = input(grid, FT)
        Et = FT \ Eω
        It = abs2.(Et)
        @test isapprox(
            abs(getceo(grid.t,real(Et.*exp.(im .* grid.ω0 .* grid.t)), It, PhysData.wlfreq(λ0))),
            ϕCEO,
            rtol=1e-10)
    end
end

@testset "CW fields" begin
    λ0 = 1064e-9
    Pavg = 20.0
    Δλ = 4e-9
    grid = Grid.EnvGrid(1.0, λ0, (980e-9, 1160e-9), 500e-12)
    energy_t = Fields.energyfuncs(grid)[1]
    x = Array{ComplexF64}(undef, length(grid.t))
    FT = FFTW.plan_fft(x, 1)
    input = Fields.CWSech(λ0=λ0, Pavg=Pavg, Δλ=Δλ, rng=MersenneTwister(0))
    Eω = input(grid, FT)
    Et = FT \ Eω
    I = Fields.It(Et, grid)
    istart = findfirst(isequal(1.0), grid.twin)
    iend = findlast(isequal(1.0), grid.twin)
    # test average power
    @test isapprox(mean(I[istart:iend]), Pavg, rtol=5e-16)
    # test coherence time
    @test isapprox(Processing.coherence_time(grid, Et), 3.35/(PhysData.c*(Δλ)/λ0^2*2π), rtol=1e-2)
    idcs = sortperm(PhysData.wlfreq.(grid.ω))
    # test spectral width
    @test isapprox(Maths.fwhm(PhysData.wlfreq.(grid.ω)[idcs], abs2.(Eω[idcs])), Δλ, rtol=3e-3)
    # now do the same for a number of realisations
    Eωs = hcat([Fields.CWSech(λ0=λ0, Pavg=Pavg, Δλ=Δλ, rng=MersenneTwister(i))(grid, FT) for i = 1:5]...)
    Iωs = abs2.(Eωs)
    Iωav = mean(Iωs, dims=2)[:,1]
    idcs = sortperm(PhysData.wlfreq.(grid.ω))
    # test average spectral width
    @test isapprox(Maths.fwhm(PhysData.wlfreq.(grid.ω)[idcs], Iωav[idcs], minmax=:max), Δλ, rtol=6e-4)
    Ets = FFTW.ifft(Eωs, 1)
    Its = abs2.(Ets)
    Itav = mean(Its[istart:iend,:])
    # test average power
    @test isapprox(Itav, Pavg, rtol=5e-16)
    # test diversity of power fluctuations
    @test mean(std(Its[istart:iend,:], dims=2)[:,1]) > 10
end

@testset "Propagation" begin
    λ0 = 800e-9
    τfwhm = 2.5e-15
    grid = Grid.RealGrid(1, λ0, (400e-9, 1200e-9), 500e-15)
    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=1e-6)
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)
    Eω = input(grid, FT)
    Eωβ1 = Fields.prop_taylor(Eω, grid, [0, 10e-15], λ0)
    Et = FT \ Eωβ1
    @test isapprox(Maths.moment(grid.t, abs2.(Maths.hilbert(Et))), 10e-15, rtol=1e-6)

    # Test sign of dispersion
    Eωβ2 = Fields.prop_taylor(Eω, grid, [0, 0, 15e-30], λ0) # positive chirp
    Et = FT \ Eωβ2
    gab = Maths.gabor(grid.t, Et, [-10e-15, 10e-15], 3e-15) # spectrogram
    ω0 = Maths.moment(grid.ω, abs2.(gab))
    @test ω0[1] < ω0[2] # mean frequency at earlier time should be lower (upchirp)

    # Test pulse stretching for Gaussian pulse
    τfwhm = 30e-15
    τ0 = Tools.τfw_to_τ0(τfwhm, :gauss)
    input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=1e-6)
    Eω = input(grid, FT)
    Eωβ2 = Fields.prop_taylor(Eω, grid, [0, 0, τ0^2], λ0) # should lead to √2 increase
    Et = FT \ Eωβ2
    τfwβ2 = Maths.fwhm(grid.t, abs2.(Maths.hilbert(Et)); method=:spline)
    τ0β2 = Tools.τfw_to_τ0(τfwβ2, :gauss)
    @test τ0β2 ≈ √2 * τ0

    # Test pulse stretching and sign of the dispersion for modal propagation
    ω0 = PhysData.wlfreq(λ0)
    # Artificial mode with τ0^2 2nd order dispersion over 2 m and α=0.1
    β(ω; z=0) = ω0/PhysData.c + 1/(0.999*PhysData.c)*(ω-ω0) + τ0^2/4*(ω-ω0)^2
    α(ω; z=0) = 0.1
    m = Modes.arbitrary(neff=Modes.neff_from_αβ(α, β))
    Eωm = copy(Eω)
    Fields.prop_mode!(Eωm, grid.ω, m, 2, λ0)
    Et = FT \ Eωm
    τfwm = Maths.fwhm(grid.t, abs2.(Maths.hilbert(Et)); method=:spline)
    τ0m = Tools.τfw_to_τ0(τfwm, :gauss)
    @test τ0m ≈ √2 * τ0
    # Check signs are correct:
    # α = 0.1 should give loss
    et, eω = Fields.energyfuncs(grid)
    @test eω(Eω)*exp(-0.2) ≈ eω(Eωm)
    # β2 > 0 should give positive chirp:
    gab = Maths.gabor(grid.t, Et, [-10e-15, 10e-15], 3e-15)
    ω0 = Maths.moment(grid.ω, abs2.(gab))
    @test ω0[1] < ω0[2]

    # Test sign of dispersion for glass
    Eωglass = Fields.prop_material(Eω, grid, :SiO2, 0.5e-3, λ0)
    Et = FT \ Eωglass
    gab = Maths.gabor(grid.t, Et, [-10e-15, 10e-15], 3e-15)
    ω0 = Maths.moment(grid.ω, abs2.(gab))
    @test ω0[1] < ω0[2]

    # Test sign of dispersion for chirped mirrors
    for mirror in (:PC70, :ThorlabsUMC)
        Eωmirr = Fields.prop_mirror(Eω, grid, mirror, 2) # one pair
        Et = FT \ Eωmirr
        gab = Maths.gabor(grid.t, Et, [-10e-15, 10e-15], 3e-15)
        ω0 = Maths.moment(grid.ω, abs2.(gab))
        @test ω0[1] > ω0[2] # negative chirp, so frequency should go down with time
    end
end

@testset "Compression" begin
# Short pulse with 100 fs^2
λ0 = 800e-9
τfwhm = 10e-15
grid = Grid.RealGrid(1, λ0, (400e-9, 1200e-9), 500e-15)
x = Array{Float64}(undef, length(grid.t))
FT = FFTW.plan_rfft(x, 1)
input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=1e-6)
Eω = input(grid, FT)
Et = FT \ Eω
Eωβ2 = Fields.prop_taylor(Eω, grid, [0, 0, 100e-30], λ0)
ϕs, Eωcomp = Fields.optcomp_taylor(Eωβ2, grid, λ0)
Etcomp = FT \ Eωcomp
@test ϕs[3] ≈ -100e-30
@test isapprox(Maths.fwhm(grid.t, abs2.(Maths.hilbert(Etcomp))), τfwhm; rtol=1e-3)

# Long pulse with 40000 fs^2 (stretches 220 fs to ~5 ps)
λ0 = 1030e-9
τfwhm = 220e-15
grid = Grid.RealGrid(1, λ0, (980e-9, 1080e-9), 20e-12)
x = Array{Float64}(undef, length(grid.t))
FT = FFTW.plan_rfft(x, 1)
input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=1e-6)
Eω = input(grid, FT)
Et = FT \ Eω
Eωβ2 = Fields.prop_taylor(Eω, grid, [0, 0, 4e-25], λ0)
ϕs, Eωcomp = Fields.optcomp_taylor(Eωβ2, grid, λ0)
Etcomp = FT \ Eωcomp
@test ϕs[3] ≈ -4e-25
@test isapprox(Maths.fwhm(grid.t, abs2.(Maths.hilbert(Etcomp))), τfwhm; rtol=1e-3)

# Short pulse with GDD and TOD
λ0 = 800e-9
τfwhm = 10e-15
grid = Grid.RealGrid(1, λ0, (400e-9, 1200e-9), 500e-15)
x = Array{Float64}(undef, length(grid.t))
FT = FFTW.plan_rfft(x, 1)
input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=1e-6)
Eω = input(grid, FT)
Et = FT \ Eω
Eωβ2 = Fields.prop_taylor(Eω, grid, [0, 0, 100e-30, 800e-45], λ0)
ϕs, Eωcomp = Fields.optcomp_taylor(Eωβ2, grid, λ0; order=3)
Etcomp = FT \ Eωcomp
@test all(ϕs .≈ [0, 0, -100e-30, -800e-45])
@test isapprox(Maths.fwhm(grid.t, abs2.(Maths.hilbert(Etcomp))), τfwhm; rtol=1e-3)

# Material insertion
λ0 = 800e-9
τfwhm = 10e-15
grid = Grid.RealGrid(1, λ0, (400e-9, 1200e-9), 500e-15)
x = Array{Float64}(undef, length(grid.t))
FT = FFTW.plan_rfft(x, 1)
input = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=1e-6)
Eω = input(grid, FT)
Et = FT \ Eω
EωFS = Fields.prop_material(Eω, grid, :SiO2, 2e-3, λ0)
d, Eωcomp = Fields.optcomp_material(EωFS, grid, :SiO2, λ0, -1e-2, 1e-2)
Etcomp = FT \ Eωcomp
@test d ≈ -2e-3
@test isapprox(Maths.fwhm(grid.t, abs2.(Maths.hilbert(Etcomp))), τfwhm; rtol=1e-3)
end

@testset "Gaussian beam initialisation" begin
    a = 16e-6
    gas = :Kr
    pres = 17.2
    τfwhm = 230e-15
    λ0 = 1030e-9
    energy = 5.2e-6
    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=2, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=3, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=4, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=3, m=1, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0, loss=false)
    )
    inputs = Fields.gauss_beam_init(modes, 2π/λ0, a*0.64, Fields.GaussField, λ0=λ0, τfwhm=τfwhm, energy=energy)
    inputs = (inputs..., ((mode=i, fields=(Fields.ShotNoise(),)) for i=1:length(modes))...)
    @test inputs[1].fields[1].energy/energy ≈ 0.9807131210817726
    @test inputs[2].fields[1].energy/energy ≈ 0.006182621678046407
    @test inputs[3].fields[1].energy/energy ≈ 0.0013567813790567626
    @test inputs[4].fields[1].energy/energy ≈ 0.0008447236094573648
    @test inputs[5].fields[1].energy/energy < 1e-20
    @test inputs[6].fields[1].energy/energy < 2e-20
    @test inputs[7].fields[1].energy/energy < 1e-20
    @test inputs[8].fields[1].energy/energy < 1e-20
end

@testset "DataField" begin
    # create a Gaussian pulse with a delay in the frequency domain, write it to a file,
    # load it back in a check it produces the correct pulse
    import DelimitedFiles: writedlm, readdlm
    import FFTW
    λlims = (200e-9, 4e-6)
    trange = 200e-15
    τfwhm = 10e-15
    τ0 = 20e-15
    λ0 = 800e-9
    f0 = PhysData.c/λ0
    σt = Maths.fwhm_to_σ(τfwhm)
    σf = 1/(4π*σt)
    f = collect(range(PhysData.c/λlims[2], PhysData.c/λlims[1]; length=2048))
    If = Maths.gauss.(f, σf; x0=f0)
    ϕ = @. -2π*τ0*(f-f0) # Fourier transform in the maths convention here--pos. delay = neg. slope
    dat = [f If ϕ]
    grid = Grid.RealGrid(1, λ0, λlims, trange)
    FT = FFTW.plan_rfft(copy(grid.t), 1)
    field = mktempdir() do td
        tf = joinpath(td, tempname())
        writedlm(tf, dat, ' ')

        @test readdlm(tf, ' ') == dat

        Fields.DataField(tf; energy=1e-6)
    end
    Eω = field(grid, FT)
    t, Et = Processing.getEt(grid, Eω)
    @test isapprox(Maths.fwhm(t, abs2.(Et)), τfwhm, rtol=1e-5)
    @test isapprox(Maths.moment(t, abs2.(Et)), τ0, rtol=1e-5)
end


text:
# ---
# title: 447. Number of Boomerangs
# id: problem447
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/number-of-boomerangs/description/>
# hidden: true
# ---
#
# You are given `n` `points` in the plane that are all **distinct** , where
# `points[i] = [xi, yi]`. A **boomerang**  is a tuple of points `(i, j, k)` such
# that the distance between `i` and `j` equals the distance between `i` and `k`
# **(the order of the tuple matters)**.
#
# Return  _the number of boomerangs_.
#
#
#
# **Example 1:**
#
#
#
#     Input: points = [[0,0],[1,0],[2,0]]
#     Output: 2
#     Explanation: The two boomerangs are [[1,0],[0,0],[2,0]] and [[1,0],[2,0],[0,0]].
#
#
# **Example 2:**
#
#
#
#     Input: points = [[1,1],[2,2],[3,3]]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: points = [[1,1]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `n == points.length`
#   * `1 <= n <= 500`
#   * `points[i].length == 2`
#   * `-104 <= xi, yi <= 104`
#   * All the points are **unique**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module LaplaceCoefficients

VERSION < v"0.4-dev" && using Docile

export laplace_coefficients_initialize
export laplace_wisdom
include("laplace_wisdom.jl")

"""
# This computes the Laplace coefficients via recursion.
"""
function initialize(jmax::Integer,alpha::Number)
const nmax=7
b0=Array{eltype(alpha)}(nmax,jmax+1) # Array to hold the coefficients
# Compute the highest two Laplace coefficients using Wisdom's series approach:
for j=0:jmax
  for i=0:nmax-1
    b0[i+1,j+1]=laplace_wisdom(1//2,j,i,alpha)/alpha^i
  end
end
return b0
end

end # module

max_stars_repo_name: nwtuchow/TTVFaster

max_issues_repo_name: nwtuchow/TTVFaster

max_forks_repo_name: nwtuchow/TTVFaster


text:
using Luna

a = 50e-6
b = 10e-6
gas = :Ar
pres = 5
L = 15e-2

τfwhm = 30e-15
λ0 = 800e-9
energy = 5e-6

grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12)

m = RectModes.RectMode(a, b, gas, pres, :Al)
aeff = let m=m
    z -> Modes.Aeff(m, z=z)
end

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),
             Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)

statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)
Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
Plotting.spec_1D(output, [5e-2, 10e-2, 11e-2])


text:
"""
    ProjectedPointDistribution(m::Manifold, d, proj!, x)

Generates a random point in ambient space of `m` and projects it to `m`
using function `proj!`. Generated arrays are of type `TResult`, which can be
specified by providing the `x` argument.
"""
struct ProjectedPointDistribution{TResult, TM<:Manifold, TD<:Distribution, TProj} <: MPointDistribution{TM}
    manifold::TM
    d::TD
    proj!::TProj
end

function ProjectedPointDistribution(M::Manifold, d::Distribution, proj!, x::TResult) where TResult
    return ProjectedPointDistribution{TResult, typeof(M), typeof(d), typeof(proj!)}(M, d, proj!)
end

function support(d::ProjectedPointDistribution)
    return MPointSupport(d.manifold)
end

function rand(rng::AbstractRNG, d::ProjectedPointDistribution{TResult}) where TResult
    x = convert(TResult, rand(rng, d.d))
    d.proj!(d.manifold, x)
    return x
end

function _rand!(rng::AbstractRNG, d::ProjectedPointDistribution, x::AbstractArray{<:Number})
    _rand!(rng, d.d, x)
    d.proj!(d.manifold, x)
    return x
end

"""
    ProjectedFVectorDistribution(type::VectorBundleFibers, x, d, project_vector!)

Generates a random vector from ambient space of manifold `type.manifold`
at point `x` and projects it to vector space of type `type` using function
`project_vector!`.
Generated arrays are of type `TResult`.
"""
struct ProjectedFVectorDistribution{TResult, TSpace<:VectorBundleFibers, TX, TD<:Distribution, TProj} <: FVectorDistribution{TSpace, TX}
    type::TSpace
    x::TX
    d::TD
    project_vector!::TProj
end

function ProjectedFVectorDistribution(type::VectorBundleFibers, x, d::Distribution, project_vector!, xt::TResult) where TResult
    return ProjectedFVectorDistribution{TResult, typeof(type), typeof(x), typeof(d), typeof(project_vector!)}(type, x, d, project_vector!)
end

function support(tvd::ProjectedFVectorDistribution)
    return FVectorSupport(tvd.type, tvd.x)
end

function rand(rng::AbstractRNG, d::ProjectedFVectorDistribution{TResult}) where TResult
    v = convert(TResult, reshape(rand(rng, d.d), size(d.x)))
    d.project_vector!(d.type, v, d.x, v)
    return v
end

function _rand!(rng::AbstractRNG, d::ProjectedFVectorDistribution, v::AbstractArray{<:Number})
    # calling _rand!(rng, d.d, v) doesn't work for all arrays types
    copyto!(v, rand(rng, d))
    return v
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

# Returns the block expression inside a :let that holds the code to be run.
# The other block (not returned) is for declaring variables in the scope of the
# let.
function _let_code_block(ex::Expr)
    @assert isexpr(ex, :let)
    return ex.args[2]
end

function _error_curly(x)
    Base.error("The curly syntax (sum{},prod{},norm2{}) is no longer supported. Expression: $x.")
end

# generates code which converts an expression into a NodeData array (tape)
# parent is the index of the parent expression
# values is the name of the list of constants which appear in the expression
function _parse_NL_expr(m, x, tapevar, parent, values)
    if isexpr(x, :call) &&
       length(x.args) >= 2 &&
       (isexpr(x.args[2], :generator) || isexpr(x.args[2], :flatten))
        header = x.args[1]
        if _is_sum(header)
            operatorid = operator_to_id[:+]
        elseif _is_prod(header)
            operatorid = operator_to_id[:*]
        else
            error("Unrecognized expression $header(...)")
        end
        codeblock = :(
            let
            end
        )
        block = _let_code_block(codeblock)
        push!(block.args, :(push!($tapevar, NodeData(CALL, $operatorid, $parent))))
        parentvar = gensym()
        push!(block.args, :($parentvar = length($tapevar)))

        code = _MA.rewrite_generator(
            x.args[2],
            t -> _parse_NL_expr(m, t, tapevar, parentvar, values),
        )
        push!(block.args, code)
        return codeblock
    end

    if isexpr(x, :call)
        if isexpr(x.args[1], :.)
            # Functions like foo.bar cannot possibly be registered, because you
            # can register only with a symbol name.
            errorstring =
                "Unexpected function $(x.args[1]). See the " *
                "documentation on how to register a function."
            return :(error($errorstring))
        end
        if _is_sum(x.args[1]) || _is_prod(x.args[1])
            opname = x.args[1]
            errorstring =
                "$opname() can appear in nonlinear expressions " *
                " only if the argument is a generator statement, for example, " *
                "$opname(x[i] for i in 1:N)."
            return :(error($errorstring))
        end
        if length(x.args) == 2 && !isexpr(x.args[2], :...) # univariate
            code = :(
                let
                end
            )
            block = _let_code_block(code)
            @assert isexpr(block, :block)
            if haskey(univariate_operator_to_id, x.args[1])
                operatorid = univariate_operator_to_id[x.args[1]]
                push!(
                    block.args,
                    :(push!($tapevar, NodeData(CALLUNIVAR, $operatorid, $parent))),
                )
            else
                opname = quot(x.args[1])
                errorstring = "Unrecognized function \"$(x.args[1])\" used in nonlinear expression."
                errorstring2 = "Incorrect number of arguments for \"$(x.args[1])\" in nonlinear expression."
                lookupcode = quote
                    if $(esc(m)).nlp_data === nothing
                        error($errorstring)
                    end
                    if !haskey(
                        $(esc(m)).nlp_data.user_operators.univariate_operator_to_id,
                        $opname,
                    )
                        if haskey(
                            $(esc(m)).nlp_data.user_operators.multivariate_operator_to_id,
                            $opname,
                        )
                            error($errorstring2)
                        else
                            error($errorstring)
                        end
                    end
                    operatorid =
                        $(esc(m)).nlp_data.user_operators.univariate_operator_to_id[$opname] +
                        _Derivatives.USER_UNIVAR_OPERATOR_ID_START - 1
                end
                push!(
                    block.args,
                    :($lookupcode;
                    push!($tapevar, NodeData(CALLUNIVAR, operatorid, $parent))),
                )
            end
            parentvar = gensym()
            push!(block.args, :($parentvar = length($tapevar)))
            push!(block.args, _parse_NL_expr(m, x.args[2], tapevar, parentvar, values))
            return code
        else
            code = :(
                let
                end
            )
            block = _let_code_block(code)
            @assert isexpr(block, :block)
            if haskey(operator_to_id, x.args[1]) # fast compile-time lookup
                operatorid = operator_to_id[x.args[1]]
                push!(block.args, :(push!($tapevar, NodeData(CALL, $operatorid, $parent))))
            elseif haskey(comparison_operator_to_id, x.args[1])
                operatorid = comparison_operator_to_id[x.args[1]]
                push!(
                    block.args,
                    :(push!($tapevar, NodeData(COMPARISON, $operatorid, $parent))),
                )
            else # could be user defined
                opname = quot(x.args[1])
                errorstring = "Unrecognized function \"$(x.args[1])\" used in nonlinear expression."
                errorstring2 = "Incorrect number of arguments for \"$(x.args[1])\" in nonlinear expression."
                lookupcode = quote
                    if $(esc(m)).nlp_data === nothing
                        error($errorstring)
                    end
                    if !haskey(
                        $(esc(m)).nlp_data.user_operators.multivariate_operator_to_id,
                        $opname,
                    )
                        if haskey(
                            $(esc(m)).nlp_data.user_operators.univariate_operator_to_id,
                            $opname,
                        )
                            error($errorstring2)
                        else
                            error($errorstring)
                        end
                    end
                    operatorid =
                        $(esc(m)).nlp_data.user_operators.multivariate_operator_to_id[$opname] +
                        _Derivatives.USER_OPERATOR_ID_START - 1
                end
                push!(
                    block.args,
                    :($lookupcode; push!($tapevar, NodeData(CALL, operatorid, $parent))),
                )
            end
            parentvar = gensym()
            push!(block.args, :($parentvar = length($tapevar)))
            for i = 1:length(x.args)-1
                arg = x.args[i+1]
                if isexpr(arg, :...)
                    if !isa(arg.args[1], Symbol)
                        error(
                            "Unexpected expression in $x. JuMP supports " *
                            "splatting only symbols. For example, x... is " *
                            "ok, but (x + 1)..., [x; y]... and g(f(y)...) " *
                            "are not.",
                        )
                    end
                    push!(
                        block.args,
                        quote
                            for val in $(esc(arg.args[1]))
                                _parse_NL_expr_runtime(
                                    $(esc(m)),
                                    val,
                                    $tapevar,
                                    $parentvar,
                                    $values,
                                )
                            end
                        end,
                    )
                else
                    push!(block.args, _parse_NL_expr(m, arg, tapevar, parentvar, values))
                end
            end
            return code
        end
    end
    if isexpr(x, :comparison)
        code = :(
            let
            end
        )
        block = _let_code_block(code)
        op = x.args[2]
        operatorid = comparison_operator_to_id[op]
        for k = 2:2:length(x.args)-1
            @assert x.args[k] == op # don't handle a <= b >= c
        end
        parentvar = gensym()
        push!(block.args, :(push!($tapevar, NodeData(COMPARISON, $operatorid, $parent))))
        push!(block.args, :($parentvar = length($tapevar)))
        for k = 1:2:length(x.args)
            push!(block.args, _parse_NL_expr(m, x.args[k], tapevar, parentvar, values))
        end
        return code
    end
    if isexpr(x, :&&) || isexpr(x, :||)
        code = :(
            let
            end
        )
        block = _let_code_block(code)
        op = x.head
        operatorid = logic_operator_to_id[op]
        parentvar = gensym()
        push!(block.args, :(push!($tapevar, NodeData(LOGIC, $operatorid, $parent))))
        push!(block.args, :($parentvar = length($tapevar)))
        push!(block.args, _parse_NL_expr(m, x.args[1], tapevar, parentvar, values))
        push!(block.args, _parse_NL_expr(m, x.args[2], tapevar, parentvar, values))
        return code
    end
    if isexpr(x, :curly)
        _error_curly(x)
    end
    if isexpr(x, :...)
        error("Unexpected splatting expression $x.")
    end
    # at the lowest level?
    return :(_parse_NL_expr_runtime($(esc(m)), $(esc(x)), $tapevar, $parent, $values))

end

function _parse_NL_expr_runtime(m::Model, x::Real, tape, parent, values)
    push!(values, x)
    push!(tape, NodeData(VALUE, length(values), parent))
    nothing
end

function _parse_NL_expr_runtime(m::Model, x::VariableRef, tape, parent, values)
    if owner_model(x) !== m
        error(
            "Variable in nonlinear expression does not belong to the " *
            "corresponding model",
        )
    end
    push!(tape, NodeData(MOIVARIABLE, x.index.value, parent))
    nothing
end

function _parse_NL_expr_runtime(m::Model, x::NonlinearExpression, tape, parent, values)
    push!(tape, NodeData(SUBEXPRESSION, x.index, parent))
    nothing
end

function _parse_NL_expr_runtime(m::Model, x::NonlinearParameter, tape, parent, values)
    push!(tape, NodeData(PARAMETER, x.index, parent))
    nothing
end

function _parse_NL_expr_runtime(m::Model, x::AbstractArray, tape, parent, values)
    error("Unexpected array $x in nonlinear expression. Nonlinear expressions may contain only scalar expressions.")
end

function _parse_NL_expr_runtime(m::Model, x::GenericQuadExpr, tape, parent, values)
    error(
        "Unexpected quadratic expression $x in nonlinear expression. " *
        "Quadratic expressions (e.g., created using @expression) and " *
        "nonlinear expressions cannot be mixed.",
    )
end

function _parse_NL_expr_runtime(m::Model, x::GenericAffExpr, tape, parent, values)
    error(
        "Unexpected affine expression $x in nonlinear expression. " *
        "Affine expressions (e.g., created using @expression) and " *
        "nonlinear expressions cannot be mixed.",
    )
end

function _parse_NL_expr_runtime(m::Model, x, tape, parent, values)
    error("Unexpected object $x (of type $(typeof(x)) in nonlinear expression.")
end

function _expression_complexity(ex::Expr)
    return isempty(ex.args) ? 1 : sum(_expression_complexity, ex.args)
end
_expression_complexity(other) = 1

# This is separated from the macro version to make it available for other @NL*
# macros.
function _process_NL_expr(model, ex)
    # This is an arbitrary cutoff. See issue #1355.
    if _expression_complexity(ex) > 5000
        @warn "Processing a very large nonlinear expression with " *
              "@NLexpression/@NLconstraint/@NLobjective. This may be very " *
              "slow. Consider using setNLobjective() and addNLconstraint() " *
              "instead of the macros or reformulating the expressions using " *
              "sum() and prod() to make them more compact. The macros are " *
              "designed to process smaller, human-readable expressions."
    end
    parsed = _parse_NL_expr(model, ex, :tape, -1, :values)
    return quote
        tape = NodeData[]
        values = Float64[]
        $parsed
        _NonlinearExprData(tape, values)
    end
end

macro _process_NL_expr(model, ex)
    return _process_NL_expr(model, ex)
end

function _Derivatives.expr_to_nodedata(
    ex::VariableRef,
    nd::Vector{NodeData},
    values::Vector{Float64},
    parentid,
    r::_Derivatives.UserOperatorRegistry,
)
    push!(nd, NodeData(MOIVARIABLE, ex.index.value, parentid))
    nothing
end

function _Derivatives.expr_to_nodedata(
    ex::NonlinearExpression,
    nd::Vector{NodeData},
    values::Vector{Float64},
    parentid,
    r::_Derivatives.UserOperatorRegistry,
)
    push!(nd, NodeData(SUBEXPRESSION, ex.index, parentid))
    nothing
end

function _Derivatives.expr_to_nodedata(
    ex::NonlinearParameter,
    nd::Vector{NodeData},
    values::Vector{Float64},
    parentid,
    r::_Derivatives.UserOperatorRegistry,
)
    push!(nd, NodeData(PARAMETER, ex.index, parentid))
    nothing
end

# Construct a _NonlinearExprData from a Julia expression.
# VariableRef objects should be spliced into the expression.
function _NonlinearExprData(m::Model, ex::Expr)
    _init_NLP(m)
    _check_expr(m, ex)
    nd, values = _Derivatives.expr_to_nodedata(ex, m.nlp_data.user_operators)
    return _NonlinearExprData(nd, values)
end
_NonlinearExprData(m::Model, ex) = _NonlinearExprData(m, :($ex + 0))

# Error if:
# 1) Unexpected expression
# 2) VariableRef doesn't match the model
function _check_expr(m::Model, ex::Expr)
    if ex.head == :ref # if we have x[1] already in there, something is wrong
        error("Unrecognized expression $ex. JuMP variable objects and input coefficients should be spliced directly into expressions.")
    end
    for e in ex.args
        _check_expr(m, e)
    end
    return
end
function _check_expr(m::Model, v::VariableRef)
    owner_model(v) === m || error("Variable $v does not belong to this model.")
    return
end
_check_expr(m::Model, ex) = nothing


text:
module CellDataTests

using Gridap
using GridapDistributed
using PartitionedArrays
using Test

function main(parts)
  @assert length(size(parts)) == 2

  output = mkpath(joinpath(@__DIR__,"output"))

  domain = (0,4,0,4)
  cells = (4,4)
  model = CartesianDiscreteModel(parts,domain,cells)
  Ω = Triangulation(model)
  Γ = Boundary(model,tags="boundary")

  f = CellField(x->x[1]+x[2],Ω)
  g = CellField(4.5,Γ)
  v = CellField(x->x,Ω)

  writevtk(Ω,joinpath(output,"Ω"),cellfields=["f"=>f])
  writevtk(Γ,joinpath(output,"Γ"),cellfields=["f"=>f,"g"=>g])

  x_Γ = get_cell_points(Γ)
  @test isa(f(x_Γ),AbstractPData)

  h = 4*f
  h = f*g
  h = (v->3*v)∘(f)
  h = ∇(f)
  h = ∇⋅v
  h = ∇×v

  dΩ = Measure(Ω,1)
  dΓ = Measure(Γ,1)

  @test sum( ∫(1)dΩ ) ≈ 16.0
  @test sum( ∫(1)dΓ ) ≈ 16.0
  @test sum( ∫(g)dΓ ) ≈ 4.5*16.0

  x_Γ = get_cell_points(dΓ)
  @test isa(f(x_Γ),AbstractPData)

end

end # module


text:
using Documenter
using MarketTechnicals

# https://juliadocs.github.io/Documenter.jl/stable/man/doctests/#Module-level-metadata
DocMeta.setdocmeta!(MarketTechnicals, :DocTestSetup, :(using MarketTechnicals);
                    recursive = true)

makedocs(
  format = Documenter.HTML(mathengine = Documenter.Writers.HTMLWriter.MathJax3()),
  sitename = "MarketTechnicals.jl",
  modules = [MarketTechnicals],
  pages = [
    "index.md",
    "getting_started.md",
    "ma.md",
    "levels.md",
    "momentum.md",
    "volatility.md",
    "volume.md",
    "candlesticks.md",
    "utils.md",
  ],
)

deploydocs(
  repo         = "github.com/JuliaQuant/MarketTechnicals.jl.git",
  devbranch    = "master",
  push_preview = true,
)


text:
using OnlineStats

@testset "cubefittable" begin

data = collect(reshape(1:20.,4,5))
axlist = [RangeAxis("XVals",1.0:4.0), CategoricalAxis("YVals",[1,2,3,4,5])]
props = Dict("att1"=>5, "att2"=>"Hallo")
a = YAXArray(axlist, data, props)

cta = CubeTable( data=a, include_axes=("XVals", "YVals"))
meancta = cubefittable(cta,Mean(),:data, by=(:YVals,))
@test meancta.data == [2.5, 6.5, 10.5, 14.5, 18.5]

ashcta = cubefittable(cta, Ash(KHist(3)), :data, by=(:YVals,))
@test all(ashcta[Hist="Frequency"][1,:] .== 0.2222222222222222)

khistcta = cubefittable(cta, KHist(3), :data, by=(:YVals,))
@test all(khistcta[Hist="Frequency"][1,:] .== 1.0)

end


text:
using LinearAlgebraicRepresentation
using Plasm
View = Plasm.view

V = [9. 13 15 17 14 13 11 9 7 5 3 0 2 2 5 7 4 12 6 8 3 5 7 8 10 11 10 13 14 13 11 9 7 4 2 12 12; 0 2 4 8 9 10 11 10 9 9 8 6 3 1 0 1 2 10 3 3 5 5 6 5 5 4 2 4 6 7 9 7 7 7 6 7 5]

EV = Array{Int64,1}[[1, 2], [1, 16], [1, 27], [2, 3], [2, 27], [2, 28], [3, 4], [3, 29], [4, 5], [4, 29], [5, 6], [5, 29], [5, 30], [6, 7], [6, 18], [7, 8], [7, 18], [8, 9], [8, 31], [8, 32], [9, 10], [9, 33], [10, 11], [10, 34], [11, 12], [11, 34], [12, 13], [12, 21], [12, 35], [13, 14], [13, 17], [13, 21], [14, 15], [14, 17], [15, 16], [15, 17], [16, 20], [16, 27], [17, 19], [18, 31], [19, 20], [19, 22], [19, 23], [19, 24], [20, 24], [20, 25], [21, 22], [21, 34], [21, 35], [22, 23], [22, 33], [23, 24], [23, 33], [24, 25], [24, 32], [25, 26], [25, 27], [25, 32], [25, 36], [26, 27], [26, 28], [26, 37], [28, 29], [28, 37], [29, 30], [29, 37], [30, 36], [31, 32], [31, 36], [32, 33], [33, 34], [34, 35], [36, 37]]

FV = Array{Int64,1}[[1, 2, 27], [6, 7, 18], [7, 8, 18, 31], [8, 31, 32], [8, 9, 32, 33], [25, 31, 32, 36], [4, 5, 29], [5, 6, 18, 30, 31, 36], [5, 29, 30], [29, 30, 36, 37], [9, 10, 33, 34], [10, 11, 34], [12, 11, 34, 35], [12, 21, 35], [21, 34, 35], [21, 22, 33, 34], [19, 22, 23], [22, 23, 33], [23, 24, 32, 33], [24, 25, 32], [12, 13, 21], [13, 17, 19, 21, 22], [19, 23, 24], [19, 20, 24], [20, 24, 25], [16, 20, 25, 27], [1, 16, 27], [25, 26, 27], [25, 26, 36, 37], [3, 4, 29], [2, 3, 28, 29], [13, 14, 17], [14, 15, 17], [15, 16, 17, 19, 20], [2, 26, 27, 28], [26, 28, 37], [37, 28, 29]]

VV = [[k] for k=1:size(V,2)]

model = (V, (VV,EV,FV))

View(Plasm.numbering(2.)(model))

View(Plasm.hpc_exploded( model )(1.5,1.5,1.5))


text:
#####################################################################
##                      Data Type Declarations                     ##
#####################################################################

type StateSpace <: LTISystem
    A::Matrix{Float64}
    B::Matrix{Float64}
    C::Matrix{Float64}
    D::Matrix{Float64}
    Ts::Float64
    nx::Int
    nu::Int
    ny::Int
    statenames::Vector{String}
    inputnames::Vector{String}
    outputnames::Vector{String}

    function StateSpace(A::Matrix{Float64}, B::Matrix{Float64},
            C::Matrix{Float64}, D::Matrix{Float64}, Ts::Float64,
            statenames::Vector{String}, inputnames::Vector{String},
            outputnames::Vector{String})
        nx = size(A, 1)
        nu = size(B, 2)
        ny = size(C, 1)

        if size(A, 2) != nx && nx != 0
            error("A must be square")
        elseif size(B, 1) != nx
            error("B must have the same row size as A")
        elseif size(C, 2) != nx
            error("C must have the same column size as A")
        elseif nu != size(D, 2)
            error("D must have the same column size as B")
        elseif ny != size(D, 1)
            error("D must have the same row size as C")
        end
        # Validate names of state, input, and output
        if size(statenames, 1) != nx
            error("Must have same number of statenames as states")
        elseif size(inputnames, 1) != nu
            error("Must have same number of inputnames as inputs")
        elseif size(outputnames, 1) != ny
            error("Must have same number of outputnames as outputs")
        end
        # Validate sampling time
        if Ts < 0 && Ts != -1
            error("Ts must be either a positive number, 0
                   (continuous system), or -1 (unspecified)")
        end
        new(A, B, C, D, Ts, nx, nu, ny, statenames, inputnames, outputnames)
    end
end

function StateSpace(A::Array, B::Array, C::Array, D::Array, Ts::Real,
        statenames::Vector{String}, inputnames::Vector{String},
        outputnames::Vector{String})
    return StateSpace(float64mat(A), float64mat(B), float64mat(C),
            float64mat(D), map(Float64, Ts), statenames, inputnames, outputnames)
end

#####################################################################
##                      Constructor Functions                      ##
#####################################################################

@doc """`ss(A,B,C,D[, Ts, statenames=..., inputnames=..., outputnames=...]) -> sys`

Create a state-space model.
This is a continuous-time model if Ts is omitted or set to 0.
Otherwise, this is a discrete-time model with sampling period Ts.
Set Ts=-1 for a discrete-time model with unspecified sampling period.

State, input and output names: each can be either a vector of strings (one string per dimension),
or a single string (e.g., "x"). In the latter case, an index is automatically appended to identify
the coordinates for each dimension (e.g. "x1", "x2", ...).

`sys = ss(D[, Ts, ...])` specifies a static gain matrix D.""" ->
function ss(A::Array, B::Array, C::Array, D::Array, Ts::Real=0; kwargs...)
    # Check the kwargs for metadata
    nu = size(B, 2)
    ny, nx = size(C, 1, 2)
    kvs = Dict(kwargs)
    statenames = validate_names(kvs, :statenames, nx)
    inputnames = validate_names(kvs, :inputnames, nu)
    outputnames = validate_names(kvs, :outputnames, ny)
    return StateSpace(A, B, C, D, Ts, statenames, inputnames, outputnames)
end

# Function for accepting scalars
function ss(A::Union{Real,Array}, B::Union{Real,Array}, C::Union{Real,Array}, D::Union{Real,Array}, args...; kwargs...)
    A = float64mat(A)
    B = float64mat(B)
    C = float64mat(C)
    if D == 0
        D = zeros(size(C,1),size(B,2))
    else
        D = float64mat(D)
    end
    ss(A, B, C, D, args..., kwargs...)
end

# Function for creation of static gain
function ss(D::Array, Ts::Real=0; kwargs...)
    ny, nu = size(D, 1, 2)
    A = zeros(0, 0)
    B = zeros(0, nu)
    C = zeros(ny, 0)
    return ss(A, B, C, D, Ts, kwargs...)
end
ss(d::Real, Ts::Real=0; kwargs...) = ss([d], Ts, kwargs...)

# ss(sys) converts to StateSpace
ss(sys::LTISystem) = convert(StateSpace, sys)

# Create a random statespace system
function rss(nx::Int, nu::Int=1, ny::Int=1, feedthrough::Bool=true)
    Q = randn(nx, nx)
    A = Q*diagm(-100*abs(randn(nx)))*Q'
    B = randn(nx, nu)
    C = randn(ny, nx)
    if feedthrough
        D = randn(ny, nu)
    else
        D = zeros(ny, nu)
    end
    return ss(A, B, C, D)
end

#####################################################################
##                         Math Operators                          ##
#####################################################################

## EQUALITY ##
function ==(s1::StateSpace, s2::StateSpace)
    fields = [:Ts, :nx, :ny, :nu, :A, :B, :C, :D, :inputnames, :outputnames,
            :statenames]
    for field in fields
        if getfield(s1, field) != getfield(s2, field)
            return false
        end
    end
    return true
end

## Approximate ##
function isapprox(s1::StateSpace, s2::StateSpace)
    fieldsApprox = [:Ts, :nx, :ny, :nu, :A, :B, :C, :D]
    fieldsEqual = [:inputnames, :outputnames, :statenames]
    for field in fieldsApprox
        if !(getfield(s1, field) ≈ getfield(s2, field))
            return false
        end
    end
    for field in fieldsEqual
        if getfield(s1, field) != getfield(s2, field)
            return false
        end
    end
    return true
end

## ADDITION ##
function +(s1::StateSpace, s2::StateSpace)
    #Ensure systems have same dimensions
    if size(s1) != size(s2)
        error("Systems have different shapes.")
    elseif s1.Ts != s2.Ts
        error("Sampling time mismatch")
    end

    A = [s1.A zeros(s1.nx, s2.nx);
         zeros(s2.nx, s1.nx) s2.A]
    B = [s1.B ; s2.B]
    C = [s1.C s2.C;]
    D = [s1.D + s2.D;]

    # Naming strategy: If only one sys is named, use that. If the names are the
    # same, use them. If the names conflict, then they are ignored, and the
    # default "" is used.
    statenames = [s1.statenames; s2.statenames]
    if all(s1.inputnames .== "")
        inputnames = s2.inputnames
    elseif all(s2.inputnames .== "") || (s1.inputnames == s2.inputnames)
        inputnames = s1.inputnames
    else
        inputnames = fill(String(""),s1.ny)
    end
    if all(s1.outputnames .== "")
        outputnames = s2.outputnames
    elseif all(s2.outputnames .== "") || (s1.outputnames == s2.outputnames)
        outputnames = s1.outputnames
    else
        outputnames = fill(String(""),s1.nu)
    end
    return StateSpace(A, B, C, D, s1.Ts, statenames, inputnames, outputnames)
end

+(s::StateSpace, n::Real) = StateSpace(s.A, s.B, s.C, s.D .+ n, s.Ts,
        s.statenames, s.inputnames, s.outputnames)
+(n::Real, s::StateSpace) = +(s, n)

## SUBTRACTION ##
-(s1::StateSpace, s2::StateSpace) = +(s1, -s2)
-(s::StateSpace, n::Real) = +(s, -n)
-(n::Real, s::StateSpace) = +(-s, n)

## NEGATION ##
-(s::StateSpace) = StateSpace(s.A, s.B, -s.C, -s.D, s.Ts, s.statenames,
        s.inputnames, s.outputnames)

## MULTIPLICATION ##
function *(s1::StateSpace, s2::StateSpace)
    #Check dimension alignment
    #Note: s1*s2 = y <- s1 <- s2 <- u
    if s1.nu != s2.ny
        error("s1*s2: s1 must have same number of inputs as s2 has outputs")
    elseif s1.Ts != s2.Ts
        error("Sampling time mismatch")
    end

    A = [s1.A    s1.B*s2.C;
         zeros(s2.nx, s1.nx)  s2.A]
    B = [s1.B*s2.D ; s2.B]
    C = [s1.C   s1.D*s2.C;]
    D = [s1.D*s2.D;]

    statenames = [s1.statenames; s2.statenames]
    return StateSpace(A, B, C, D, s1.Ts, statenames, s2.inputnames,
            s1.outputnames)
end

*(s::StateSpace, n::Real) = StateSpace(s.A, s.B, s.C*n, s.D*n, s.Ts,
        s.statenames, s.inputnames, s.outputnames)
*(n::Real, s::StateSpace) = *(s, n)

## DIVISION ##
/(s1::StateSpace, s2::StateSpace) = s1*(1/s2)

function /(n::Real, s::StateSpace)
    # Ensure s.D is invertible
    Dinv = try
        inv(s.D)
    catch
        error("D isn't invertible")
    end
    return StateSpace(s.A - s.B*Dinv*s.C, s.B*Dinv, -n*Dinv*s.C, n*Dinv, s.Ts,
            s.statenames, s.outputnames, s.inputnames)
end

/(s::StateSpace, n::Real) = StateSpace(s.A, s.B, s.C/n, s.D/n, s.Ts,
        s.statenames, s.inputnames, s.outputnames)

#####################################################################
##                       Indexing Functions                        ##
#####################################################################
Base.ndims(::StateSpace) = 2
Base.size(s::StateSpace) = (s.ny, s.nu)
Base.size(s::StateSpace, d) = d <= 2 ? size(s)[d] : 1

function Base.getindex(s::StateSpace, inds...)
    if size(inds, 1) != 2
        error("Must specify 2 indices to index statespace model")
    end
    rows, cols = ControlSystems.index2range(inds...)
    return StateSpace([s.A;], [s.B[:, cols];], [s.C[rows, :];], [s.D[rows, cols];],
            s.Ts, [s.statenames;], [s.inputnames[cols];], [s.outputnames[rows];])
end

#####################################################################
##                        Display Functions                        ##
#####################################################################

# TODO : this is a very hacky way of handling StateSpace printing.
function _string_mat_with_headers(X::Matrix, cols::Vector{String},
                                rows::Vector{String})
    mat = [[""] reshape(cols,1,length(cols));
           rows X]
    p = (io, m) -> Base.showarray(io, m, false, header=false)
    return replace(sprint(p, mat), "\"", " ")
end

Base.print(io::IO, s::StateSpace) = show(io, s)

function Base.show(io::IO, s::StateSpace)
    # Compose the name vectors
    inputs = format_names(s.inputnames, "u", "?")
    outputs = format_names(s.outputnames, "y", "?")
    println(io, "StateSpace:")
    if s.nx > 0
        states = format_names(s.statenames, "x", "?")
        println(io, "A = \n", _string_mat_with_headers(s.A, states, states))
        println(io, "B = \n", _string_mat_with_headers(s.B, inputs, states))
        println(io, "C = \n", _string_mat_with_headers(s.C, states, outputs))
    end
    println(io, "D = \n", _string_mat_with_headers(s.D, inputs, outputs), "\n")
    # Print sample time
    if s.Ts > 0
        println(io, "Sample Time: ", s.Ts, " (seconds)")
    elseif s.Ts == -1
        println(io, "Sample Time: unspecified")
    end
    # Print model type
    if s.nx == 0
        print(io, "Static gain")
    elseif iscontinuous(s)
        print(io, "Continuous-time state-space model")
    else
        print(io, "Discrete-time state-space model")
    end
end


text:
using Test

@info "CANalyze.Utils tests..."
@testset "endian" begin
    using CANalyze.Utils

    @testset "is_little_or_big_endian" begin
        is_little = is_little_endian()
        is_big = is_big_endian()

        @test (is_little || is_big) == true
    end

    @testset "is_little_and_big_endian" begin
        is_little = is_little_endian()
        is_big = is_big_endian()

        @test (is_little && is_big) == false
    end
end

@testset "convert" begin
    using CANalyze.Utils

    @testset "to_bytes_1" begin
        value = 1337
        new_value = from_bytes(typeof(value), to_bytes(value))

        @test value == new_value
    end

    @testset "from_bytes_1" begin
        types = [UInt8, UInt16, UInt32, UInt64, UInt128]
        for (i, type) in enumerate(types)
            array = [UInt8(j) for j=1:2^(i-1)]
            new_array = to_bytes(from_bytes(type, array))

            @test array == new_array
        end
    end

    @testset "from_bytes_2" begin
        types = [Int8, Int16, Int32, Int64, Int128]
        for (i, type) in enumerate(types)
            array = [UInt8(j) for j=1:2^(i-1)]
            new_array = to_bytes(from_bytes(type, array))

            @test array == new_array
        end
    end

    @testset "from_bytes_4" begin
        types = [Float16, Float32, Float64]
        for (i, type) in enumerate(types)
            array = [UInt8(j) for j=1:2^i]
            new_array = to_bytes(from_bytes(type, array))

            @test array == new_array
        end
    end
end

@testset "mask" begin
    using CANalyze.Utils

    @testset "zero_mask" begin
        @test zero_mask(UInt8) == zero(UInt8)
        @test zero_mask(UInt16) == zero(UInt16)
        @test zero_mask(UInt32) == zero(UInt32)
        @test zero_mask(UInt64) == zero(UInt64)
        @test zero_mask(UInt128) == zero(UInt128)
    end

    @testset "full_mask_1" begin
        @test full_mask(UInt8) == UInt8(0xFF)
        @test full_mask(UInt16) == UInt16(0xFFFF)
        @test full_mask(UInt32) == UInt32(0xFFFFFFFF)
        @test full_mask(UInt64) == UInt64(0xFFFFFFFFFFFFFFFF)
        @test full_mask(UInt128) == UInt128(0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF)
    end

    @testset "full_mask_2" begin
        @test mask(UInt8) == UInt8(0xFF)
        @test mask(UInt16) == UInt16(0xFFFF)
        @test mask(UInt32) == UInt32(0xFFFFFFFF)
        @test mask(UInt64) == UInt64(0xFFFFFFFFFFFFFFFF)
        @test mask(UInt128) == UInt128(0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF)
    end

    @testset "mask_1" begin
        @test mask(UInt8, UInt8(0)) == UInt8(0b0)
        @test mask(UInt8, UInt8(1)) == UInt8(0b1)
        @test mask(UInt8, UInt8(2)) == UInt8(0b11)
        @test mask(UInt8, UInt8(3)) == UInt8(0b111)
        @test mask(UInt8, UInt8(4)) == UInt8(0b1111)
        @test mask(UInt8, UInt8(5)) == UInt8(0b11111)
        @test mask(UInt8, UInt8(6)) == UInt8(0b111111)
        @test mask(UInt8, UInt8(7)) == UInt8(0b1111111)
        @test mask(UInt8, UInt8(8)) == UInt8(0b11111111)
    end

    @testset "mask_2" begin
        @test mask(UInt16, UInt8(0)) == UInt16(0b0)
        value = UInt16(2)
        for i in 1:16
            @test mask(UInt16, UInt8(i)) == (value - 1)
            value *= 2
        end
    end

    @testset "mask_3" begin
        @test mask(UInt32, UInt8(0)) == UInt32(0b0)
        value = UInt16(2)
        for i in 1:32
            @test mask(UInt32, UInt8(i)) == (value - 1)
            value *= 2
        end
    end

    @testset "mask_4" begin
        @test mask(UInt64, UInt8(0)) == UInt64(0b0)
        value = UInt64(2)
        for i in 1:64
            @test mask(UInt64, UInt8(i)) == (value - 1)
            value *= 2
        end
    end

    @testset "mask_5" begin
        @test mask(UInt128, UInt8(0)) == UInt128(0b0)
        value = UInt128(2)
        for i in 1:16
            @test mask(UInt128, UInt8(i)) == (value - 1)
            value *= 2
        end
    end

    @testset "shifted_mask_1" begin
        for i in 0:8
            @test mask(UInt8, i, 0) == mask(UInt8, i)
        end
    end

    @testset "shifted_mask_2" begin
        for i in 0:16
            @test mask(UInt16, i, 0) == mask(UInt16, i)
        end
    end

    @testset "shifted_mask_3" begin
        for i in 0:32
            @test mask(UInt32, i, 0) == mask(UInt32, i)
        end
    end

    @testset "shifted_mask_4" begin
        for i in 0:64
            @test mask(UInt64, i, 0) == mask(UInt64, i)
        end
    end

    @testset "shifted_mask_5" begin
        for i in 0:128
            @test mask(UInt128, i, 0) == mask(UInt128, i)
        end
    end

    @testset "bit_mask_1" begin
        for T in [UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128]
            s = 8*sizeof(T) - 1
            @test bit_mask(T, 0:s) == full_mask(T)
        end
    end

    @testset "bit_mask_2" begin
        for T in [UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128]
            s = 8*sizeof(T) - 1
            @test bit_mask(T, s) == mask(T, 1, s)
        end
    end
end


text:
@testset "Quadruped" begin
    mech = get_mechanism(:quadruped,
        timestep=0.05,
        gravity=-9.81,
        friction_coefficient=0.8,
        damper=1000.0,
        spring=30.0)

    initialize!(mech, :quadruped)

    storage = simulate!(mech, 5.0,
            record=true,
            verbose=false)

    res = get_sdf(mech, storage) # distance from floor to each contact
    @test minimum(minimum([min.(0.0, r) for r in res])) >= 0.0
end

@testset "Box toss" begin
    for timestep in [0.10, 0.05, 0.01, 0.005]
        mech = get_mechanism(:block,
            timestep=timestep,
            gravity=-9.81,
            friction_coefficient = 0.1)

        initialize!(mech, :block,
            x=[0.0, 0.0, 0.5],
            v=[1.0, 1.5, 1.0],
            ω=[5.0, 4.0, 2.0] .* timestep)
        storage = simulate!(mech, 5.0,
            record=true,
            opts=SolverOptions(btol=1e-6, rtol=1e-6,
            verbose=false))

        @test norm(storage.v[1][end], Inf) < 1.0e-8
        @test norm(storage.x[1][end][3] - 0.25, Inf) < 1.0e-3
    end
end

@testset "Four-bar linkage" begin
    for timestep in [0.10, 0.05, 0.01, 0.005]
        mech = Dojo.get_mechanism(:fourbar,
            model="fourbar",
            timestep=timestep)
        Dojo.initialize!(mech, :fourbar,
            angle=0.25)
        loopjoints = mech.joints[end:end]
        Dojo.root_to_leaves_ordering(mech) == [2, 7, 3, 6, 1, 8, 4, 9]

        # Simulation
        function ctrl!(m, t)
            Dojo.set_input!(m, 1.0 * m.timestep * SVector(rand(), -rand(), 0.0, 0.0, 0.0))
            return nothing
        end
        storage = Dojo.simulate!(mech, 5.0, ctrl!, verbose=false, record=true)

        min_coords = Dojo.get_minimal_coordinates(mech)
        @test norm(min_coords[5] - +min_coords[4], Inf) < 1.0e-5
        @test norm(min_coords[5] - -min_coords[3], Inf) < 1.0e-5
        @test norm(min_coords[5] - (min_coords[2] - min_coords[1]), Inf) < 1.0e-5
    end
end

@testset "Tennis racket" begin
    # Simulation
    timestep=0.01
    gravity=0.0
    mech = Dojo.get_mechanism(:dzhanibekov,
            timestep=timestep,
            gravity=gravity);

    # Simulate
    Dojo.initialize!(mech, :dzhanibekov,
        angular_velocity=[15.0; 0.01; 0.0])
    storage = Dojo.simulate!(mech, 4.0,
        record=true,
        verbose=false)

    # The x position of the side body oscillate between a positive and negative
    # value when we observe the Dhanibekov effect. Otherwise, it always remains in the positive position.
    @test Dojo.minimum([x[1] for x in storage.x[2]]) < -0.05
end


text:
"""
    discretize(domain, N)

Discretize domain with `N` points.
"""
function discretize end

discretize(domain::ClosedIntervalDomain, N) = LinRange(domain.left, domain.right, N + 1)

function discretize(domain::PeriodicIntervalDomain, N)
    LinRange(domain.left + (domain.right - domain.left) / N, domain.right, N)
end


text:
#=
file-io.jl
test AVSfldIO through FileIO
=#

using FileIO: File, @format_str
using FileIO: add_format, load, save, UUID
#using FileIO: info, unknown, query
#using AVSfldIO

data1 = reshape(Int32.(1:20), 4, 5)
file = tempname() * ".fld"

#ff = File{format"AVSfld"}(file)
#query(ff)

add_format(format"AVSfld", "# AVS", ".fld",
    [:AVSfldIO => UUID("b6189060-daf9-4c28-845a-cc0984b81781")])

save(file, data1)
data2 = load(file)
@assert data2 == data1

rm(file, force=true)

max_stars_repo_name: JeffFessler/AVSfldIO

max_issues_repo_name: JeffFessler/AVSfldIO


text:
using SimpleDiffEq, Test

@time begin
  @time @testset "Discrete Tests" begin include("discrete_tests.jl") end
  @time @testset "SimpleEM Tests" begin include("simpleem_tests.jl") end
end


text:
using QMTK
using Compat.Test

@testset "Sites" begin
    include("Sites.jl")
end

@testset "SubSites" begin
    include("SubSites.jl")
end

@testset "Tags" begin
    include("Tag/Tags.jl")
end

@testset "Ket" begin
    include("Ket.jl")
end

@testset "Space" begin
    include("Space/Space.jl")
end


text:
const InitializedReference{T} = Union{Base.RefValue{T}, Pointer{T}}

@inline deref(x::Number) = x
@inline deref(x::InitializedReference) = x[]
@inline RESERVED_INCREMENT_SEED_RESERVED!(c::AbstractUninitializedReference, a) = (c[] =  deref(a); nothing)
@inline function RESERVED_INCREMENT_SEED_RESERVED!(c::InitializedReference, a)
    c[] = SIMDPirates.vadd(c[], deref(a)); nothing
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(c::AbstractUninitializedReference, a, b)
    c[] = SIMDPirates.vmul(deref(a), deref(b)); nothing
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(c::InitializedReference, a, b)
    c[] = SIMDPirates.vmuladd(deref(a), deref(b), c[]); nothing
end
# @inline RESERVED_DECREMENT_SEED_RESERVED!(c::Uninitialized, a) = (c[] = -a; nothing)
# @inline function RESERVED_DECREMENT_SEED_RESERVED!(c, a)
    # c[] = SIMDPirates.vsub(c[], a); nothing
# end
# @inline function RESERVED_DECREMENT_SEED_RESERVED!(c::Uninitialized, a, b)
    # c[] = - SIMDPirates.vmul(a, b); nothing
# end
# @inline function RESERVED_DECREMENT_SEED_RESERVED!(c, a, b)
    # c[] = SIMDPirates.fnmadd(a, b, c[]); nothing
# end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::AbstractUninitializedReference, A::AbstractArray
)
    c[] = sum(A)
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::InitializedReference, A::AbstractArray
)
    c[] = Base.FastMath.add_fast(c[], sum(A))
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::AbstractUninitializedReference, A::AbstractArray, B::AbstractArray
)
    c[] = dot(A, B)
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::InitializedReference, A::AbstractArray, B::AbstractArray
)
    c[] = Base.FastMath.add_fast(c[], dot(A, B))
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::AbstractUninitializedReference, A, B::AbstractArray
)
    c[] = Base.FastMath.mul_fast(deref(A), sum(B))
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::InitializedReference, A, B::AbstractArray
)
    c[] = Base.FastMath.add_fast(c[], Base.FastMath.mul_fast(deref(A), sum(B)))
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::AbstractUninitializedReference, A::AbstractArray, B
)
    c[] = Base.FastMath.mul_fast(deref(B), sum(A))
end
@inline function RESERVED_INCREMENT_SEED_RESERVED!(
    c::InitializedReference, A::AbstractArray, B
)
    c[] = Base.FastMath.add_fast(c[], Base.FastMath.mul_fast(deref(B), sum(A)))
end


text:
@doc raw"""
    PowerKernel([γ=1])

The Power Kernel is a negative definite kernel given by:

```math
\kappa(\mathbf{x},\mathbf{y}) =
\|\mathbf{x} - \mathbf{y} \|^{2\gamma}
\qquad \gamma \in (0,1]
```
"""
struct PowerKernel{T<:AbstractFloat} <: NegativeDefiniteKernel{T}
    gamma::HyperParameter{T}
    PowerKernel{T}(γ::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
    )
end
PowerKernel(γ::T1 = 1.0) where {T1<:Real} = PowerKernel{floattype(T1)}(γ)

@inline powerkernel(z::T, γ::T) where {T<:AbstractFloat} = z^γ

@inline basefunction(::PowerKernel) = SquaredEuclidean()
@inline kappa(κ::PowerKernel{T}, z::T) where {T} = powerkernel(z, getvalue(κ.gamma))


text:
# is caching useful?  only in extreme cases, apparently
# (but we may need to define equality on states!)

function slow(n)
#    matcher = Repeat(Repeat(Equal("a"), 0, n), n, 0)
#    matcher = Seq(Repeat(Equal("a"), 0, ), Repeat(Equal("a"), 0, n))
#    for greedy in (true, false)
    for greedy in (true,)  # false just eats memory
        println("greedy $greedy")
        matcher = Repeat(Equal("a"), 0, n; greedy=greedy)
        for i in 1:n
            matcher = Seq(Repeat(Equal("a"), 0, n; greedy=greedy), matcher)
        end
        source = repeat("a", n)
        for config in (NoCache, Cache)
            println("$(config)")
            all1 = make_all(config)
            @time collect(all1(source, matcher))
            @time n = length(collect(all1(source, matcher)))
            println("n results: $n")
            debug, all2 = make(Debug, source, matcher; delegate=config)
            collect(all2)
            println("max depth: $(debug.max_depth)")
            println("max iter: $(debug.max_iter)")
            println("n calls: $(debug.n_calls)")
        end
    end
end
slow(3)
#slow(7)  # not for travis!

println("slow ok")


text:
abstract type AbstractZeroInflatedLink end

"""
    encprob(link::AbstractZeroInflatedLink, proc1, proc2)

Probability of encounter (non-zero observation) given linear predictors proc1
and proc2. Depending on the link, proc2 may not be required as an
argument.
"""
function encprob end

"""
    posrate(link::AbstractZeroInflatedLink, proc1, proc2)

Expected positive observation conditional on an encounter given linear
predictors proc1 and proc2. Depending on the link, proc1 may be dropped as an
argument.
"""
function posrate end

"""
    LogitLogLink <: AbstractZeroInflatedLink

A link function where for observation y modeled by linear predictors p₁ and p₂,

p(y > 0) = logit⁻¹(p₁)

𝔼[y ∣ y > 0] = exp(p₂) - b
"""
struct LogitLogLink <: AbstractZeroInflatedLink end

encprob(lll::LogitLogLink, proc1) = logistic(proc1)
encprob(lll::LogitLogLink, proc1, proc2) = encprob(lll, proc1)
posrate(lll::LogitLogLink, proc2) = exp(proc2)
posrate(lll::LogitLogLink, proc1, proc2) = posrate(lll, proc2)

"""
    PoissonLink <: AbstractZeroInflatedLink

The Poisson link function of Thorson (2018), where for observation y modeled by
linear predictors n (log numbers-density) and w (log weight-per-group),

p(y > 0) = 1 - exp(-a exp(n))

𝔼[y ∣ y > 0] = exp(n) exp(w) / p(y > 0)

where a is an optional offset (e.g. area swept).

Encounter probability uses the complementary log-log link. This link
approximates a compound Poisson-gamma (Tweedie with 1 < p < 2) distribution.
This implies that number of groups sampled is drawn from a Poisson distribution
with mean exp(n) and each group has expected weight exp(w).
"""
struct PoissonLink{T} <: AbstractZeroInflatedLink
    offset::T

    function PoissonLink(offset::T) where T
        offset > 0 || throw(DomainError("offset must be positive"))
        new{T}(offset)
    end
end
function PoissonLink()
    # Default to integer offset so that it doesn't inadvertantly convert e.g.
    # Float32 to Float64
    PoissonLink(one(Int))
end

encprob(pl::PoissonLink, logn) = ccdf(Poisson(pl.offset * exp(logn)), 0)
encprob(pl::PoissonLink, logn, logw) = encprob(pl, logn)
posrate(pl::PoissonLink, logn, logw) = exp(logn + logw) / encprob(pl, logn)

"""
    IdentityLink <: AbstractZeroInflatedLink

A simple link function where the first process is the probability of encounter
and the second is expected non-zero observation.
"""
struct IdentityLink <: AbstractZeroInflatedLink end

function encprob(il::IdentityLink, proc1)
    0 < proc1 < 1 || throw(DomainError("proc1 must be between zero and one."))
    proc1
end
encprob(il::IdentityLink, proc1, proc2) = encprob(il, proc1)
posrate(il::IdentityLink, proc2) = proc2
posrate(il::IdentityLink, proc1, proc2) = posrate(il, proc2)


text:
## winston.jl  Winston plotting routines for ROC and DET
## (c) 2014--2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

@require Winston begin
    import Winston

    function rocplot(r::Roc, nr=1; ch=false, title="ROC", xlabel="Pfa", ylabel="Pmiss")
        col = string("krgmcb"[(nr-1) % 6 + 1])
        lty = ["-", "--", ";"][div(nr-1,6) + 1]
        if ch
            chi = find(r.ch)
            x, y = r.pfa[chi], r.pmiss[chi]
        else
            x, y = r.pfa, r.pmiss
        end
        if nr==1
            p = Winston.plot(x, y, lty*col)
            Winston.xlabel(xlabel)
            Winston.ylabel(ylabel)
            Winston.title(title)
        else
            p = Winston.oplot(x, y, lty * col)
        end
        return p
    end
    Winston.plot(r::Roc, nr=1; ch=false, title="ROC", xlabel="Pfa", ylabel="Pmiss") = rocplot(r::Roc, nr; title=title, xlabel=xlabel, ylabel=ylabel)

    function detplot(r::Roc, nr=1; xlabel="Pfa (%)", ylabel="Pmiss (%)", title="DET plot")
        if nr==1
            grid = [0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 40]
            ticklabels = AbstractString[]
            for x in grid
                if x<1
                    push!(ticklabels, @sprintf("%3.1f", x))
                else
                    push!(ticklabels, @sprintf("%1.0f", x))
                end
            end
            grid ./= 100
            p = Winston.plot()
            Winston.setattr(p, :aspect_ratio, 1) # essential for DET plots!
            Winston.title(title)
            Winston.xlim(qnorm(0.8e-3), qnorm(0.55))
            Winston.ylim(qnorm(0.8e-3), qnorm(0.55))
            for axis in (p.x1, p.y1)
                Winston.setattr(axis, "draw_subticks", false)
                Winston.setattr(axis, "ticks", qnorm(grid))
                Winston.setattr(axis, "ticklabels", ticklabels)
                Winston.setattr(axis, "tickdir", 1)
                Winston.setattr(axis, "draw_grid", true)
            end
            for axis in (p.x2, p.y2)
                Winston.setattr(axis, "draw_ticks", false)
            end
            Winston.xlabel(xlabel)
            Winston.ylabel(ylabel)
        end
        col = string("krgmcb"[(nr-1) % 6 + 1])
        lty = ["-", "--", ";"][div(nr-1,6) + 1]
        p = Winston.oplot(qnorm(r.pfa), qnorm(r.pmiss), lty * col)
        return p
    end

    function llrplot(r::Roc; title="LLR plot", xlabel="score", ylabel="LLR")
        mi = max(minimum(r.llr), minimum(r.θ))
        ma = min(maximum(r.llr), maximum(r.θ))
        ran = [mi, ma]
        p = Winston.plot(r.θ, r.llr)
        Winston.xlim(ran)
        Winston.ylim(ran)
        Winston.title(title)
        Winston.xlabel(xlabel)
        Winston.ylabel(ylabel)
        return p
    end

    function apeplot(r::Roc; xmin=-7, xmax=7, title="Applied Probability of Error", xlabel="prior log odds", ylabel="Bayes error rate")
        lo = collect(xmin:0.01:xmax)
        be = ber(r, lo)
        mbe = minber(r, lo)
        defbe = 1 ./ normfactor(lo)
        p = Winston.plot(lo, be, "-r")
        Winston.oplot(lo, mbe, "-g")
        Winston.oplot(lo, defbe, "--k")
        Winston.ylim(0, 1.1*maximum(be))
        Winston.title(title)
        Winston.xlabel(xlabel)
        Winston.ylabel(ylabel)
        return p
    end

    function nbeplot(r::Roc; xmin = -10, xmax = 5, title="Normalized Bayes' Error", xlabel="prior log odds", ylabel="normalized Bayes error")
        lo = collect(xmin:0.01:xmax)
        norm = normfactor(lo)
        bfa, bmiss = [be .* norm for be in ber_famiss(r, lo)]
        nbe = bfa + bmiss
        mbfa, mbmiss = [mbe .* norm for mbe in minber_famiss(r, lo)]
        mnbe = minber(r, lo) .* norm
        p = Winston.plot(lo, nbe, "-r")
        Winston.oplot(lo, mnbe, "-g")
        Winston.oplot(lo, bfa, "--r")
        Winston.oplot(lo, bmiss, "-.r")
        Winston.oplot(lo, mbfa, "--g")
        Winston.oplot(lo, mbmiss, "-.g")
        Winston.ylim(0, 1.1)
        Winston.title(title)
        Winston.xlabel(xlabel)
        Winston.ylabel(ylabel)
        return p
    end

end


text:
export Exponent, show

"""
    Exponent(expo::SortedDict{Variable, Degree})

Define a mathematical exponent, that is a product of `Variable` and conjugated
`Variable`.

### Attributes
- `expo` : a dictionary associating `Variable` objects to a `Degree` objects,
their exponent and the exponent of their conjugates.
- `degree`: a `Degree` object indicating the global degree of the exponent :
the sum of the explicit variables degrees and of the conjugated variables
degrees.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > Exponent(SortedDict(a=>Degree(1,0), b=>Degree(1,1))) == a*b*conj(b)
```
"""
struct Exponent <: AbstractPolynomial
    expo::SortedDict{Variable, Degree}
    degree::Degree

    function Exponent(expo::SortedDict{Variable, Degree})
        degexpl, degconj = 0,0
        for (var, degree) in expo
            ((degree.explvar < 0) || (degree.conjvar < 0)) && error(LOGGER, "Exponent(): Expected non negative exponent for variable $var (got $degree)")
            (isreal(var) && degree.conjvar != 0) && error(LOGGER, "Exponent(): Expected nul conj exponent for real variable $var (got $degree)")
            (isbool(var) && degree.explvar ∉ SortedSet([0,1])) && error(LOGGER, "Exponent(): Expected boolean exponent for bool $var (got $degree)")
            if degree != Degree(0,0)
                degexpl = max(degexpl, degree.explvar)
                degconj = max(degconj, degree.conjvar)
            else
                delete!(expo, var)
            end
        end
        return new(expo, Degree(degexpl, degconj))
    end
end

#############################
## Constructors
#############################
Exponent() = Exponent(SortedDict{Variable, Degree}())
Exponent(x::Variable) = Exponent(SortedDict{Variable, Degree}(x=>Degree(1,0)))

#############################
## Equality, hash
#############################
function ==(exp1::Exponent, exp2::Exponent)
  (length(exp1.expo) == length(exp2.expo)) || return false
  return isequal(exp1.expo, exp2.expo)
  # for (varname, deg) in exp2.expo
  #   if !haskey(exp1.expo, varname) || exp1.expo[varname] != deg
  #     return false
  #   end
  # end
  # return true
end
!=(exp1::Exponent, exp2::Exponent) = !(exp1 == exp2)
hash(expo::Exponent, h::UInt) = hash(expo.degree, hash(expo.expo, h))

#############################
## Show
#############################
function Base.show(io::IO, exp::Exponent)
  if exp == Exponent()
    print(io, "1")
  else
    expo = exp.expo
    i=length(expo)
    for (var, deg) in expo
      if var.kind <: Complex && deg.conjvar>0
        print(io, "conj(", var, ")")
        if deg.conjvar > 1
          print(io, "^", deg.conjvar)
        end
        if deg.explvar > 0
          print(io, " * ")
        end
      end
      if deg.explvar == 1
        print(io, var)
      elseif deg.explvar > 1
        print(io, var, "^", deg.explvar)
      end
      if i > 1
        print(io, " * ")
      end
      i -= 1
    end
  end
end


text:
using ExtensibleScheduler
using Test

include("test_priority.jl")
include("test_job.jl")
include("test_jobstores.jl")
include("trigger/test_triggers.jl")
include("test_clock.jl")
include("samples/test_samples.jl")
include("test_action.jl")


text:
# I want to try a pure MC simulation that doesn't rely on a neural network and
# see how things
using Statistics: mean, quantile
using StatsBase: countmap
using Game2048
using Game2048: bitboard_to_array, DIRS
using Serialization: serialize, deserialize
using Base.Threads
using BSON
using BenchmarkTools

### load all results are analysis
full_results = mapreduce(deserialize, vcat, readdir("c:/data/game2048/raw/", join=true));

ls = length.(full_results)

quantile(ls, 0.90)

normalize_bitboard(bb) = Float32.(bb ./ 11)

function make_matrix_x(grvec::Vector{GameRecorder})
    res = zeros(Float32, 4, 4, sum(length.(grvec)))
    i = 1
    for gr in grvec
        for state in gr.states
            @inbounds res[:, :, i] .= normalize_bitboard(bitboard_to_array(state))
            i += 1
        end
    end
    reshape(res, 4, 4, 1, :)
end

function make_matrix_y(grvec::Vector{GameRecorder})
    moves = mapreduce(gr->gr.moves, vcat, grvec)
    onehotbatch(moves, DIRS)
end

# filter out the best results
good_results = filter(x->x.final_state |> maximum >= 11,  full_results);

length(good_results)

@time x = make_matrix_x(good_results);
@time y = make_matrix_y(good_results);

serialize("x-data.jls", x)
serialize("y-data.jls", y)